; Given a single AR id #: ; - finds all MDI images containing that AR ; - copies from MDI database and renames fits files ; - uses pe and fastrack modules to make projection and extract region ; - applies find_nl program to get a neutral line for the AR images ; - export data files of qualitative characteristics: total flux, effective distance, primary and total neutral line length ; - use make_movie to make movie of results @/auto/home2/jmason86/Forecast/pros/find_nl_n.pro @/auto/home2/jmason86/Forecast/pros/posplot.pro @/auto/home2/jmason86/Forecast/pros/gifplot.pro @/auto/home2/jmason86/Forecast/pros/sel_AR.pro @/auto/home2/jmason86/Forecast/pros/WLSG.pro @/auto/home2/jmason86/Forecast/pros/plotj.pro @/auto/home2/jmason86/Forecast/pros/mapout.pro @/auto/home2/jmason86/Forecast/pros/clean.pro @/auto/home2/jmason86/Forecast/pros/bzcorrect.pro @/auto/home2/jmason86/Forecast/pros/mdiquality.pro @/auto/home2/jmason86/Forecast/pros/count_1min5minmag.pro @/auto/home2/jmason86/Forecast/pros/mktable.pro @/auto/home2/jmason86/Forecast/pros/tiltonly.pro pro solignis, no_copy=no_copy, out_typ=out_typ, dat_out=dat_out, gif=gif, ps=ps, loud=loud, typ=typ, lat=lat, lon=lon, area=area, lng_ext=lng_ext, lbl=lbl, src=src, z=z, long_cut=long_cut, fix_dir=fix_dir, yn=yn, nar=nar, sorc=sorc, seldate=seldate, debug=debug IF n_elements(debug) eq 0 then debug = 0 ; create and set location of working directory IF debug then goto, debuger ;print, "Where would you like your output saved?(full directory name)" ;fix_dir=STRING('') & read, fix_dir fix_dir='/auto/home2/jmason86' cd, fix_dir & spawn, 'mkdir -p solignis/data' working = fix_dir + STRING('/solignis/data') cd, working loc = STRING('/auto/home2/jmason86/Forecast/pros') ; Ask if user wants movie output ;print, "Would you like a movie of the output? y/n" ;yn = STRING('') & read, yn yn='y' ;determine NAR to run on if n_elements(noaa) eq 0 then begin print, "Select ARs by: 1) NOAA AR#'s, 2) Date Range (>1 week ago), 3) Most Recent ARs 4) Range of ARs" choice=0 & read, choice if choice eq 1 then goto, by_AR & if choice eq 2 then goto, by_date & if choice eq 3 then goto, by_yester & if choice eq 4 then goto, by_rng by_AR: iii=0 print, 'How many ARs?:' totar = 0 & read, totar noaa_all = intarr(totar) print, 'Input NOAA AR#:' read, noaa_all data_loc = STRING('/mag') goto, done by_date: iii=0 print, 'Input start date. Format = yyyy-MON-dd' seldate = STRING('') & read, seldate print, 'Input end date. Format = yyyy-MON-dd' seldate2 = STRING('') & read, seldate2 data_loc = STRING('/mag') narset = get_nar(seldate,seldate2) noaa_all = narset.noaa goto, done by_yester: print, 'Recent data compiling...' seldate = STRING('') data_loc = STRING('/synoptic/eof') narset = get_nar(seldate) noaa_all = narset.noaa goto, done by_rng: iii=0 print, 'Input first NOAA AR#:' firar = 0 & read, firar print, 'Input last NOAA AR#;' lasar = 0 & read, lasar totar = lasar - firar + 1 noaa_all = intarr(totar) for alpha = 0, totar-1 DO BEGIN noaa_all(alpha) = firar + alpha ENDFOR data_loc = '/mag' goto, done debuger: working = '/auto/home2/jmason86/solignis/solignis/data' loc = '/auto/home2/jmason86/Forecast/pros' yn = 'n' noaa_all=7999 data_loc = '/mag' cd, working done: endif ;NAR loop print, 'NOAA Active regions to run:', noaa_all last = (n_elements(noaa_all)-1) for iii = 0,last do begin noaa = noaa_all[iii] print,'Program running on AR:', noaa defaults: ; set call parameter defaults if n_elements(no_copy) eq 0 then no_copy=0 if n_elements(out_typ) eq 0 then out_typ=1 if n_elements(dat_out) eq 0 then dat_out=1 if n_elements(gif) eq 0 then gif=1 if n_elements(ps) eq 0 then ps=0 if n_elements(loud) eq 0 then loud=0 if n_elements(typ) eq 0 then typ=0 if n_elements(z) eq 0 then z=0 if n_elements(long_cut) eq 0 then long_cut=40 if typ eq 1 and n_elements(src) eq 0 then no_copy=1 if gif eq 0 and ps eq 0 then begin & gif=1 & endif lbl=strcompress(string(fix(noaa)),/remove_all) spawn, 'rm -rf '+lbl spawn, 'mkdir '+lbl ; clean up leftover temporary files if (no_copy eq 0) then clean else clean, no_src=1 ; create directories to store temporary files in if (no_copy eq 0) then spawn, 'mkdir '+lbl+'/fits' spawn, 'mkdir out' if n_elements(src) ne 0 then begin & spawn, 'cp '+src+'/*.fits ./'+lbl+'/fits/' & no_copy=1 & endif spawn, 'mkdir '+lbl+'/calc' ; find the range of mdi images that the AR appears in ; loud turns on feedback messages in sel_AR program ; long_cut sets the boundaries of longitude the AR must be inside if (typ eq 0) then begin sel_AR, noaa, srt_dt, subsq, edpt, errdat, dims, loud=loud, long_cut=long_cut if errdat eq 1 then begin error = 'Not contained in MDI data' goto, fin2 endif endif else if (typ eq 1) then begin dims=fltarr(6) if n_elements(lat) eq 0 then begin & print, "latitude:" & lat=0.0 & read, lat & endif dims(2)=lat if n_elements(lon) eq 0 then begin & print, "carrington longitude:" & lon=0.0 & read, lon & endif dims(3)=lon if n_elements(area) eq 0 then begin & print, "area (in 1e-6 * FD):" & area=0.0 & read, area & endif dims(1)=area if n_elements(lng_ext) eq 0 then begin & print, "longitudinal extent in degrees:" & lng_ext=0.0 & read, lng_ext & endif dims(0)=lng_ext if n_elements(lbl) eq 0 then begin & print, "output folder label" & lbl='' & read, lbl & endif errdat=0 endif if errdat eq 1 then goto, fin ; using time_index module, converts the start date found to a SOHO era day if (no_copy eq 0) then begin dt=strcompress(string(srt_dt(0)),/remove_all)+'.'+strcompress(string(srt_dt(1)),/remove_all)+'.'+strcompress(string(srt_dt(2)),/remove_all) spawn, 'rm -f tmp' spawn, 'time_index -d time='+dt+' > tmp' close, 5 openr,5,'tmp' srt_num=0 & readf,5,srt_num close,5 num=srt_num idx=0 & ctr=1 & tst=0 subsq=[subsq, 0, 0] if subsq(0) eq 0 then begin & idx=2 & endif ; copies mdi fits files in whole day blocks from database to temporary folder src repeat begin number=strcompress(string(num),/remove_all) spawn, 'ls /synoptic/eof/fd_M_96m_01d.00'+number+'/fd_M_96m_01d.'+number+'.00*.fits >'+lbl+'dlist_mvp' log_name=+lbl+'dlist_mvp' close, 50 & openr,50,log_name while not eof(50) do begin data_loc=STRING('/synoptic/eof') endwhile close, 50 print, 'cp '+data_loc+'/fd_M_96m_01d.00'+number+'/fd_M_96m_01d.'+number+'.00*.fits ./'+lbl+'/fits/' spawn, 'cp '+data_loc+'/fd_M_96m_01d.00'+number+'/fd_M_96m_01d.'+number+'.00*.fits ./'+lbl+'/fits/' if ctr lt subsq(idx) then begin & ctr=ctr+1 & num=num+1 & endif else begin ctr=1 & idx=idx+1 if (n_elements(subsq) gt idx) then num=num+subsq(idx) & idx=idx+1 & endelse if n_elements(subsq) le idx then tst=1 else if subsq(idx) eq 0 then tst=1 endrep until tst eq 1 ; cuts out files where AR is outside boundaries in first and last days' files id=-1 & if edpt(0) gt 0 then repeat begin id=id+1 if id lt 10 then idn='0'+strcompress(string(id),/remove_all) else idn=strcompress(string(id),/remove_all) ; print, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(srt_num),/remove_all)+'.00'+idn+'.fits' ;spawn, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(srt_num),/remove_all)+'.00'+idn+'.fits' endrep until id eq edpt(0)-1 id=15 & if edpt(1) lt 14 then repeat begin id=id-1 if id lt 10 then idn='0'+strcompress(string(id),/remove_all) else idn=strcompress(string(id),/remove_all) print, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(num),/remove_all)+'.00'+idn+'.fits' spawn, 'rm -f ./'+lbl+'/fits/fd_M_96m_01d.'+strcompress(string(num),/remove_all)+'.00'+idn+'.fits' endrep until id eq edpt(1)+1 endif ;if no_src eq 1 then begin & del_stuff=1 & goto, no_files & endif ; makes directory to store output in spawn, 'mkdir '+lbl ;Determine good data/drop bad data ;list fits files to get mdi_id from file name spawn, 'ls ./'+lbl+'/fits/*.fits >> ./'+lbl+'/fits/list' openr, 6, './'+lbl+'/fits/list' array = strarr(100000) & ar=0L & ar_end=ar while not eof(6) do begin file = '' & readf, 6, file array(ar) = file fits = readfits(array(ar),hd) DATAFILE = strarr(27) d=fxpar(hd,'DATAFILE') mdi_id = strmid(d, 13, 4) fits_id = strmid(d, 18, 4) t=fxpar(hd,'T_REC') time = strmid(t,0,23) ;determine mdi quality record='/auto/PAS/D136007/mdi_log/lev1.5/fd_M_96m_01d/00'+mdi_id+'.record.rdb' info=strarr(17) mag_qual=strarr(15) & mag_vld=strarr(15) mag_qual(0)='NoFile' openr, 1, record, error=errout if errout eq 0 then begin readf, 1, info quality = 0L & qualvld = 0L & mag_qual=0L & mag_vld=0L pos1=strpos(info(0),'QUALITY') val=strpos(info(0),'QUAL_VLD') timepos=strpos(info(0),time) for j=2, 16 do begin time2=strmid(info(j), timepos, 23) if time EQ time2 then begin mag_qual=strmid(info(j), pos1, 10) reads, mag_qual, quality, format='(%"%x")' mag_vld=strmid(info(j), val, 10) reads, mag_vld, qualvld, format='(%"%x")' endif endfor ; j loop maskstring = '0x402c0c04' & qualmask = 0L reads, maskstring, qualmask, format='(%"%x")' endif bad = 0 if ((quality AND qualmask AND qualvld) NE 0) then bad = 1 if bad eq 1 then begin spawn, 'mkdir ./'+lbl+'/fits/bad_data' spawn, 'mv ./'+lbl+'/fits/fd_M_96m_01d.'+mdi_id+'.'+fits_id+'.fits ./'+lbl+'/fits/bad_data' print, 'Removing bad data: ', d endif close,1 ar=ar+1L & ar_end=ar-1L endwhile & close, 6 spawn, 'rm ./'+lbl+'/fits/list' ;end determine good data/drop bad data ; gets list of files in src folder no_src=1 spawn,'\rm '+lbl+'dlist_mvp' spawn,'ls ./'+lbl+'/fits/*.* > '+lbl+'dlist_mvp' close, 5 openr,5,+lbl+'dlist_mvp' dfname = strarr(10000) & dfp=0L & dfp_end=dfp while not eof(5) do begin no_src=0 data_file='' & readf, 5, data_file dfname(dfp) = data_file dfp=dfp+1L & dfp_end=dfp-1L endwhile & close,5 del_stuff = 0 & if no_src eq 1 then begin & del_stuff=1 & goto, no_files & endif ; starts output of analysis data close, 77 & openw, 77, ('./'+lbl+'/'+lbl+'_NL_length.dat') & close, 77 close, 78 & openw,78,('./'+lbl+'/'+lbl+'_pos_flux_sum.dat') & close, 78 close, 79 & openw,79,('./'+lbl+'/'+lbl+'_neg_flux_sum.dat') & close, 79 close, 80 & openw,80,('./'+lbl+'/'+lbl+'_av_NL_grad.dat') & close, 80 close, 81 & openw,81,('./'+lbl+'/'+lbl+'_total_flux.dat') & close, 81 close, 82 & openw,82,('./'+lbl+'/'+lbl+'_eff_dist.dat') & close, 82 close, 83 & openw,83,('./'+lbl+'/'+lbl+'_ang_on_disc.dat') & close, 83 ; analyses the size of the AR over its existence and picks a good size ; for the region of analysis to extract, maintains consistent sample size ; over analysis run data = readfits(dfname(0), hd) xscale=fxpar(hd,'XSCALE') yscale=fxpar(hd,'YSCALE') if (n_elements(where(dims))/5-1) gt 0 then if n_elements(dims(0,*)) gt 1 then dims=dims(*,0:n_elements(where(dims))/5-1) ;if (typ eq 0) then begin ;place=0 & for i=0, n_elements(subsq)-1 do begin ;if (i mod 2 eq 1) and (i ne 1) and (subsq(i) ne 0) then dims(5,place)=1 ;if i mod 2 eq 0 then place=place+subsq(i) ;endfor ;endif sub=0 & if n_elements(dims(0,*)) gt 1 then long_ext=max(dims(0,0:n_elements(dims(0,*))-2), sub) else long_ext=dims(0) area=dims(1,sub) & lat=dims(2,sub) width=(long_ext>1) * 8.55 * cos(lat*3.14159/180) * xscale/2 + 50 if area le 40 then begin err_dat = 1 print, "skipping: size too small" error = '' & error = 'size too small' endif else err_dat = 0 if err_dat eq 1 then goto, fin2 area=(area > 10)*xscale*yscale height=area/width+100 day_ctr=intarr(2) & day_ctr(0)=-1 & day_ctr(1)=0 ; begins loops that processes each fits file separately for ip=0L, dfp_end do begin ; reads fits file data = readfits(dfname(ip), hd) ; reads date and time from header d=fxpar(hd,'T_REC') if day_ctr(1) ne fix(strmid(d,8,2)) then begin & day_ctr(0)=day_ctr(0)+1 & day_ctr(1)=fix(strmid(d,8,2)) & endif dt=strmid(d, 0, 4)+'-'+strmid(d,5,2)+'-'+strmid(d,8,2) tm=fix(strmid(d,11,2))*60+fix(strmid(d,14,2)) xscale=fxpar(hd,'XSCALE') yscale=fxpar(hd,'YSCALE') tmpname=fxpar(hd,'DATAFILE') ; gets NOAA AR data for date ;if (dims(5,day_ctr(0)) eq 1) or (dims(4,day_ctr(0)) eq 2) or ((dims(4,day_ctr(0)) eq 1) and (tm eq 0)) then begin ;lat=dims(2,day_ctr(0)) & lon=dims(3,day_ctr(0)) ;endif else begin ;lat=dims(2,day_ctr(0))+((dims(2,day_ctr(0)+1)-dims(2,day_ctr(0)))*(tm/1440.)) ;lon=dims(3,day_ctr(0))+(dims(3,day_ctr(0)+1)-dims(3,day_ctr(0)))*(tm/1440.) ;endelse lat=median(dims(2,0:max(where(dims(2,*))))) lon=median(dims(3,0:max(where(dims(3,*))))) ; assigns file number if ip lt 10 then nm='000'+strcompress(string(ip),/remove_all) else if ip lt 100 then nm='00'+strcompress(string(ip),/remove_all) else if ip lt 1000 then nm='0'+strcompress(string(ip),/remove_all) else nm=strcompress(string(ip),/remove_all) ; corrects Bz values of pixels in image using helioscopic angle bzcorrect, +working+'/'+lbl+'/fits/'+tmpname spawn, 'mkdir '+working+'/'+lbl+'out/' ; creates tst.map file for running pe module mapout, fn=+lbl+'tst', in= +working+'/'+lbl+'/fits/'+tmpname, out= +working+'/'+lbl+'out/', rows=fix(height), cols=fix(width), lat=lat, lon=lon ; runs project module using pe and tst.map file spawn, 'pe map='+lbl+'tst.map' ; renames pe output so reiteration will not overwrite cd, working spawn, 'mv ./'+lbl+'out/0000.fits ./'+lbl+'out/'+lbl+'PD'+nm+'.fits' ; reads data out of pe output file data=readfits('./'+lbl+'out/'+lbl+'PD'+nm+'.fits', hdr) ; reads filename and time fn=fxpar(hd, 'DATAFILE') if typ eq 0 then fn=strmid(fn, 13, 9) else fn=strmid(fn, 9, 10) t_obs=fxpar(hd, 'T_OBS') ; puts data in mdi_src mdi_src=data spawn, 'rm '+lbl+'tst.map' ; first smoothed run of find_nl program if loud then begin & print, "first run" & endif find_nl_n, hd, mdi, mdi_src, no_line_1, nl_xy_1, nl_info_1, edpt_1, coords_1, detail=6, sep=30, minsize=15, thresh=20, pbt_thresh=120, isth_cut=5, loud=loud, z=z ; second run: compares first smoothed run to less smoothed iteration if loud then begin & print, "second run" & endif if no_line_1 eq 0 then find_nl_n, hd, mdi, mdi_src, no_line_2, nl_xy_2, nl_info_2, edpt_2, coords_2, nl_xy_3, nl_info_3, detail=3, sep=30, minsize=5, thresh=8, pbt_thresh=100, isth_cut=3, comp=1, nl_xy_comp=nl_xy_1, nl_info_comp=nl_info_1, loud=loud, z=z else no_line_2=1 ; third run: compares iterative output of second to unsmoothed data ; to create best guess for nl if loud then begin & print, "third run" & endif if no_line_2 eq 0 then find_nl_n, hd, mdi, mdi_src, no_line_3, nl_xy_4, nl_info_4, edpt_3, coords_3, nl_xy_5, nl_info_5, nl1, pbt, detail=1, sep=30, minsize=5, thresh=6, pbt_thresh=100, isth_cut=2, comp=1, nl_xy_comp=nl_xy_3, nl_info_comp=nl_info_3, loud=loud, z=z else no_line_3=1 ; runs an isthmus cutting routine on nl guess if (no_line_1 eq 0) and (no_line_2 eq 0) and (no_line_3 eq 0) then begin isth, nl_info_5, nl_xy_5, nl1, isth_cut=1, thresh=10, pbt_thresh=100 ; analysis data output routines ; 77. Neutral line length ; 78. Positive sum of flux above noise level (100 Gauss) ; 79. Negative sum of flux above noise level ; 80. Average gradient across neutral line ; 81. Absolute value of total flux above noise level ; 82. Effective distance between flux weighter centers ; 83. Central angle from 0 degrees longitude and latitude in degrees if (dat_out eq 1) then begin openu,77,('./'+lbl+'/'+lbl+'_NL_length.dat'),/append openu,78,('./'+lbl+'/'+lbl+'_pos_flux_sum.dat'),/append openu,79,('./'+lbl+'/'+lbl+'_neg_flux_sum.dat'),/append openu,80,('./'+lbl+'/'+lbl+'_av_NL_grad.dat'),/append openu,81,('./'+lbl+'/'+lbl+'_total_flux.dat'),/append openu,82,('./'+lbl+'/'+lbl+'_eff_dist.dat'),/append openu,83,('./'+lbl+'/'+lbl+'_ang_on_disc.dat'),/append openu,84,('./'+lbl+'/'+lbl+'_central_latitude.dat'),/append openu,85,('./'+lbl+'/'+lbl+'_central_longitude.dat'),/append openu,86,('./'+lbl+'/'+lbl+'_tilt_angle.dat'),/append printf, 84, fn, ': ', lat & close, 84 printf, 85, fn, ': ', lon & close, 85 close, 99 &openw, 99, +lbl+'/calc/'+lbl+'_xy_coords',/append printf, 99, nl_xy_4 close, 99 sz=size(mdi_src) dxmin=0 & dymin=0 & dxmax=sz(1)-1 & dymax=sz(2)-1 ; Total length of primary neutral line length=0. for o=0, n_elements(nl_info_5)-1 do begin if nl_info_5(o).typ ne 2 then begin for p=nl_info_5(o).offset, nl_info_5(o).offset+nl_info_5(o).n-2 do begin if (nl_xy_5(0,p)+nl_xy_5(0,p+1))/2 ge dxmin and (nl_xy_5(0,p)+nl_xy_5(0,p+1))/2 le dxmax and (nl_xy_5(1,p)+nl_xy_5(1,p+1))/2 ge dymin and (nl_xy_5(1,p)+nl_xy_5(1,p+1))/2 le dymax then $ ;James: average positions are within fits image data length=length+((nl_xy_5(0,p)-nl_xy_5(0,p+1))^2 + (nl_xy_5(1,p)-nl_xy_5(1,p+1))^2)^.5 ;James: solid trig, nothing wrong with method here endfor endif else if (nl_xy_5(0,nl_info_5(o).offset)+nl_xy_5(0,nl_info_5(o).n))/2 ge dxmin and (nl_xy_5(0,nl_info_5(o).offset)+nl_xy_5(0,nl_info_5(o).n))/2 le dxmax and (nl_xy_5(1,nl_info_5(o).offset)+nl_xy_5(1,nl_info_5(o).n))/2 ge dymin and (nl_xy_5(1,nl_info_5(o).offset)+nl_xy_5(1,nl_info_5(o).n))/2 le dymax then $ length=length+((nl_xy_5(0,nl_info_5(o).offset)-nl_xy_5(0,nl_info_5(o).n))^2 + (nl_xy_5(1,nl_info_5(o).offset)-nl_xy_5(1,nl_info_5(o).n))^2)^.5 endfor printf,77,fn,': ',length ; Average gradient across neutral line sum=0. & length=0. & for q=0, n_elements(nl_info_5)-1 do begin if nl_info_5(q).typ ne 2 then begin for p=nl_info_5(q).offset, nl_info_5(q).offset+nl_info_5(q).n-2 do begin med=fltarr(2) & nmed=fltarr(2) & smed=fltarr(2);James: med = mid-point of hypotenuse med(0)=(nl_xy_5(0,p)+nl_xy_5(0,p+1))/2. ;James: x-average position med(1)=(nl_xy_5(1,p)+nl_xy_5(1,p+1))/2. ;James: y-average position if med(0) ge dxmin and med(0) le dxmax and med(1) ge dymin and med(1) le dymax then begin ;make sure med is inside fits data image angle=atan((abs(nl_xy_5(0,p)-nl_xy_5(0,p+1)))/(abs(nl_xy_5(1,p)-nl_xy_5(1,p+1)))) ;James: trig to determine the angle between x and y directions angle=abs(3.14159/2-angle) ;James: rotate angle ccw (so that nmed,smed point across NL instead of along) nmed=med+[(.1*sin(angle)),(.1*cos(angle))] & smed=med-[(.1*sin(angle)),(.1*cos(angle))];James: points a distance away from med a=interpolate(mdi_src,nmed(0),nmed(1)) & b=interpolate(mdi_src,smed(0),smed(1));James: get B values at points on either side of med if ((a ge 0) and (b le 0)) or ((a le 0) and (b ge 0)) then begin ;James: if a and b are opposite (bipolar flux) len=((nl_xy_5(0,p)-nl_xy_5(0,p+1))^2 + (nl_xy_5(1,p)-nl_xy_5(1,p+1))^2)^.5 ;James: hypotenuse sum=sum+len*(abs(a-b)) ;James: length of hypotenuse * difference in B field length=length+len ;James: total the lengths of each small hypotenuse endif endif endfor endif else begin ;James: same process as above med=fltarr(2) & nmed=fltarr(2) & smed=fltarr(2) med(0)=(nl_xy_5(0,nl_info_5(q).offset)+nl_xy_5(0,nl_info_5(q).n))/2. med(1)=(nl_xy_5(1,nl_info_5(q).offset)+nl_xy_5(1,nl_info_5(q).n))/2. if med(0) ge dxmin and med(0) le dxmax and med(1) ge dymin and med(1) le dymax then begin angle=atan((abs(nl_xy_5(0,nl_info_5(q).offset)-nl_xy_5(0,nl_info_5(q).n)))/(abs(nl_xy_5(1,nl_info_5(q).offset)-nl_xy_5(1,nl_info_5(q).n)))) angle=abs(3.14159/2-angle) nmed=med+[(.1*sin(angle)),(.1*cos(angle))] & smed=med-[(.1*sin(angle)),(.1*cos(angle))] a=interpolate(mdi_src,nmed(0),nmed(1)) & b=interpolate(mdi_src,smed(0),smed(1)) if ((a ge 0) and (b le 0)) or ((a le 0) and (b ge 0)) then begin len=((nl_xy_5(0,nl_info_5(q).offset)-nl_xy_5(0,nl_info_5(q).n))^2 + (nl_xy_5(1,nl_info_5(q).offset)-nl_xy_5(1,nl_info_5(q).n))^2)^.5 sum=sum+len*(abs(a-b)) length=length+len endif endif endelse endfor av=sum/length printf,80, fn,": ", av ; Summed flux ;James method using mask ;read pixel size from fits file (in arcsec) angarc = fxpar(hd, 'CDELT1') pixel = angarc * 725e5 ;conversion to cm pos_mask = mdi_src GE 100 neg_mask = mdi_src LE -100 pos = (mdi_src * pos_mask) neg = (mdi_src * neg_mask) ;Natty method using more junk ;pos=mdi_src(dxmin:dxmax,dymin:dymax) * (mdi_src(dxmin:dxmax,dymin:dymax) gt 100.0) ;neg=mdi_src(dxmin:dxmax,dymin:dymax) * (mdi_src(dxmin:dxmax,dymin:dymax) lt -100.0) possum=total(pos) * pixel^2 negsum=total(neg) * pixel^2 totalsum=(total(pos)+abs(total(neg)))*pixel^2 printf,81,fn,': ',totalsum printf,78,fn,': ',possum printf,79,fn,': ',negsum ; Effective distance between flux weighted centers cen=fltarr(4) & cen(*)=0. for tmx=dxmin, dxmax do begin for tmy=dymin, dymax do begin cen(0)=pos(tmx,tmy)*tmx+cen(0) cen(1)=pos(tmx,tmy)*tmy+cen(1) cen(2)=neg(tmx,tmy)*tmx+cen(2) cen(3)=neg(tmx,tmy)*tmy+cen(3) endfor endfor cen(0:1)=cen(0:1)/possum & cen(2:3)=cen(2:3)/possum effd_tmp=((cen(0)-cen(2))^2+(cen(1)-cen(3))^2)^(.5) effd = effd_tmp * 2.1025E16 ;this number puts effd back into units of pixels printf,82,fn,': ',effd ; location of center of active region bin on visible solar surface in deg dtt=dt+strmid(d,10,9) centlon=TIM2CARR(dtt) longitude=lon-centlon if longitude gt 90 then longitude=longitude-360. if longitude lt -90 then longitude=longitude+360. centang=(acos( cos(longitude*0.0174532925) * cos(lat*0.0174532925) ))/0.0174532925 printf,83,fn,': ',centang ;Yang's Tilt Angle IF lat GT 0 THEN tiltonly, mdi_src, tiltangle, /north ELSE $ tiltonly, mdi_src, tiltangle printf,86,fn,': ', tiltangle close, 77 close, 78 close, 79 close, 80 close, 81 close, 82 close, 83 close, 86 endif ; outputs images in AR designated folder if (out_typ) eq 0 then begin if ps then posplot, mdi_src, nl_xy_5, nl_info_5, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl), t_obs=t_obs if gif then gifplot, mdi_src, nl_xy_5, nl_info_5, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl+'_sub'),t_obs=t_obs endif else if out_typ eq 1 then begin if ps then posplot, mdi_src, nl_xy_4, nl_info_4, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl), nl_xy_2=nl_xy_3, nl_info_2=nl_info_3, nl_xy_3=nl_xy_5, nl_info_3=nl_info_5,t_obs=t_obs if gif then gifplot, mdi_src, nl_xy_4, nl_info_4, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl), nl_xy_2=nl_xy_3, nl_info_2=nl_info_3, nl_xy_3=nl_xy_5, nl_info_3=nl_info_5,t_obs=t_obs endif else if out_typ eq 2 then begin if ps then posplot, mdi_src, nl_xy_4, nl_info_4, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl+'_1'), nl_xy_2=nl_xy_3, nl_info_2=nl_info_3, nl_xy_3=nl_xy_5, nl_info_3=nl_info_5,t_obs=t_obs if gif then gifplot, mdi_src, nl_xy_4, nl_info_4, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl+'_1'), nl_xy_2=nl_xy_3, nl_info_2=nl_info_3, nl_xy_3=nl_xy_5, nl_info_3=nl_info_5,t_obs=t_obs if ps then posplot, mdi_src, nl_xy_2, nl_info_2, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl+'_2'), nl_xy_2=nl_xy_1, nl_info_2=nl_info_1, nl_xy_3=nl_xy_3, nl_info_3=nl_info_3,t_obs=t_obs if gif then gifplot, mdi_src, nl_xy_2, nl_info_2, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl+'_2'), nl_xy_2=nl_xy_1, nl_info_2=nl_info_1, nl_xy_3=nl_xy_3, nl_info_3=nl_info_3,t_obs=t_obs endif endif else begin print, "no nl found in ", fn+'.AR'+lbl if (dat_out eq 1) then begin openu,77,('./'+lbl+'/'+lbl+'_NL_length.dat'),/append openu,78,('./'+lbl+'/'+lbl+'_pos_flux_sum.dat'),/append openu,79,('./'+lbl+'/'+lbl+'_neg_flux_sum.dat'),/append openu,80,('./'+lbl+'/'+lbl+'_av_NL_grad.dat'),/append openu,81,('./'+lbl+'/'+lbl+'_total_flux.dat'),/append openu,82,('./'+lbl+'/'+lbl+'_eff_dist.dat'),/append openu,83,('./'+lbl+'/'+lbl+'_ang_on_disc.dat'),/append sz=size(mdi_src) dxmin=0 & dymin=0 & dxmax=sz(1)-1 & dymax=sz(2)-1 ; Total length of primary neutral line length=0. printf,77,fn,': ',length ; Average gradient across neutral line printf,80, fn," : N/A" ; Summed flux ;James method using mask ;read pixel size from fits file (in arcsec) angarc = fxpar(hd, 'CDELT1') pixel = angarc * 725e5 ;conversion to cm pos_mask = mdi_src GE 100 neg_mask = mdi_src LE -100 pos = (mdi_src * pos_mask) neg = (mdi_src * neg_mask) possum=total(pos) * pixel^2 negsum=total(neg) * pixel^2 totalsum=(total(pos)+abs(total(neg)))*pixel^2 printf,81,fn,': ',totalsum printf,78,fn,': ',possum printf,79,fn,': ',negsum ; Effective distance between flux weighted centers cen=fltarr(4) & cen(*)=0. for tmx=dxmin, dxmax do begin for tmy=dymin, dymax do begin cen(0)=pos(tmx,tmy)*tmx+cen(0) cen(1)=pos(tmx,tmy)*tmy+cen(1) cen(2)=neg(tmx,tmy)*tmx+cen(2) cen(3)=neg(tmx,tmy)*tmy+cen(3) endfor endfor cen(0:1)=cen(0:1)/possum & cen(2:3)=cen(2:3)/possum effd=((cen(0)-cen(2))^2+(cen(1)-cen(3))^2)^(.5) printf,82,fn,': ',effd ; location of center of active region bin on visible solar surface in deg dtt=dt+strmid(d,10,9) centlon=TIM2CARR(dtt) longitude=lon-centlon if longitude gt 90 then longitude=longitude-360. if longitude lt -90 then longitude=longitude+360. centang=(acos( cos(longitude*0.0174532925) * cos(lat*0.0174532925) ))/0.0174532925 printf,83,fn,': ',centang close, 77 close, 78 close, 79 close, 80 close, 81 close, 82 close, 83 nl_xy_4=0 & nl_info_4=0 & coords=0 & edpt=0 & no_lines=1 if ps then posplot, mdi_src, nl_xy_4, nl_info_4, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl), no_lines=no_lines, t_obs=t_obs if gif then gifplot, mdi_src, nl_xy_4, nl_info_4, coords, edpt, ('./'+lbl+'/'+fn+'.AR'+lbl),no_lines=no_lines, t_obs=t_obs endif endelse endfor no_files: if no_src eq 1 then begin print, "No files found for ", lbl error = 'No files found AR' err_dat = 1 goto, fin2 endif else err_dat = 0 if del_stuff eq 1 then begin spawn, 'rm -f ./'+lbl+'/*.*' spawn, 'rm -rf ./'+lbl endif if yn EQ 'no' then goto, fin & if yn EQ 'n' then goto, fin ; From make_movie.pro ; makes collection of gif files into a movie set_plot, 'z' ; gets list of files in target folder spawn,'rm '+loc+'/'+lbl+'dlist_mvp' spawn, 'rm -rf make_movie_tmp/*' spawn, 'mkdir make_movie_tmp' spawn, 'cp '+lbl+'/*.gif make_movie_tmp/' spawn,'ls make_movie_tmp/*.gif > '+loc+'/'+lbl+'dlist_mvp' close, 5 openr,5, +loc+'/'+lbl+'dlist_mvp' dfname = strarr(10000) & dfp=0L & dmax=0 while not eof(5) do begin data_file='' & readf, 5, data_file dfname(dfp) = data_file dmax=dfp dfp=dfp+1L endwhile & close,5 if dmax eq 0 then goto, fin dfname=dfname(0:dmax) labels=strarr(n_elements(dfname)) if n_elements(sorc) eq 0 then begin ; add label for i=0L, dmax do begin if n_elements(sorc) eq 0 then begin nm1=strmid(dfname(i),15,4) nm2=strmid(dfname(i),20,4) path=data_loc+'/fd_M_96m_01d.00'+nm1+'/fd_M_96m_01d.'+nm1+'.'+nm2+'.fits' endif else begin ; assign labels from other source? ; print, dfname(i) ; wt='' & read, wt endelse data=readfits(path,hd) labels(i)=fxpar(hd,'T_OBS') read_gif,dfname(i),image sz=size(image) r=indgen(256) & g=r & b=g tvlct, r, g, b device, set_resolution=[sz(1),sz(2)] ;window, 1, xsize=sz(1), ysize=sz(2) tv,image label=strmid(labels(i),0,10)+' '+strmid(labels(i),11,8)+' '+strmid(labels(i),20,3) xyouts,sz(1)-400, 30, label, color=255, charsize=3, charthick=2.5, /device loadct, 3 tvlct, r,g,b,/get g(255)=190 & b(255)=100 & r(255)=0 tvlct,r,g,b write_gif, dfname(i), tvrd(), r,g,b endfor device, /close set_plot, 'x' ;wdelete, 1 endif spawn, 'rm -f /scr/*' spawn, 'mkdir /scr' image2movie, dfname, movie_name=lbl+'.mpg', /gif, scratchdir=''+working+'/scr/' spawn, 'mv '+lbl+'.mpg '+lbl+'/'+lbl+'.mpg' print, 'movie moved to: '+working+'/' +lbl+'/'+lbl+'.mpg' spawn, 'rm -rf make_movie_tmp/*' ;end make_movie if (no_copy eq 0) then clean else clean, no_src=1 fin: WLSG, lbl plotj, lbl, loc, dt mktable, lbl, lat, lon, height, width fin2: if n_elements(err_dat) eq 0 then err_dat = 0 & if err_dat eq 1 then begin cd, working spawn, 'rm -rf '+lbl close, 10 & if iii eq 0 then openw, 10, 'skipped_regions' else $ openw, 10, 'skipped_regions', /append printf, 10, lbl, ' ', error & print, lbl, ' ', error close, 10 endif spawn, 'pwd' spawn, 'rm -rf ../fits' cd, working spawn, 'mkdir junk' & spawn, 'mv *.param *.gif *tmp *mvp junk' spawn, 'rm -rf junk' finito: if errdat eq 1 then print, 'No AR data found for '+lbl+' in NOAA database' else print, '************Program normal completion************' endfor ;iii loop END