pro IRIS_PREP_TREND_AIA, $ win = win, png = png, inpath = inpath, daymedian = daymedian, $ pstop = pstop, loud = loud, makepngs = makepngs, $ save = save, outpath = outpath, chisqlim = chisqlim, numlim = numlim ;+ ; ; Read the output of various IRIS_PREP AIAcorr logs ; ;- tt0 = SYSTIME(/sec) if N_ELEMENTS(outpath) eq 0 then outpath = '/net/xema/Volumes/disk2/data/iris/iris_prep/trend/' if N_ELEMENTS(inpath) eq 0 then inpath = '/irisa/data/prep/' if N_ELEMENTS(chisqlim) eq 0 then chisqlim = 100 if N_ELEMENTS(numlim) eq 0 then numlim = 100 if N_ELEMENTS(win) eq 0 then win = 10 ;+++++++++++++++++++++++++++++++++++++++++++++++++ ; Read the results of the per-OBS sine wave fits ; to the AIA cross-correlation into a structure ; array (obsdata) ;------------------------------------------------- RESTGEN, file = CONCAT_DIR(outpath, 'iris_prep_filter_aia_save.genx'), aiafits numobs = N_ELEMENTS(aiafits) all_obsdata = IRIS_READ_META() for i = 0, numobs - 1 do begin metamatch = WHERE(all_obsdata.path eq aiafits[i].filename, nummatch) if nummatch gt 0 then begin this_obsdata = CREATE_STRUCT(aiafits[i], all_obsdata[metamatch]) if N_ELEMENTS(obsdata) eq 0 then obsdata = this_obsdata else obsdata = [obsdata, this_obsdata] endif endfor taistart = ANYTIM2TAI(FILE2TIME(obsdata.path)) obstag = TAG_NAMES(obsdata) ;+++++++++++++++++++++++++++++++++++++++++++++++++ ; Read the results of the individual image cross- ; correlation into a structure array (imgdata) ;------------------------------------------------- RESTGEN, file = CONCAT_DIR(outpath, 'iris_prep_all_aia_save.genx'), imgdata imgtags = TAG_NAMES(imgdata) xtag = WHERE(imgtags eq 'X1') if KEYWORD_SET(loud) then PRINT, 'IRIS_PREP_TREND_AIA loaded data. Elapsed time : ', $ STRING(SYSTIME(/sec) - tt0, form = '(f10.1)') if KEYWORD_SET(pstop) then STOP ;+++++++++++++++++++++++++++++++++++++++++++++++++ ; Set up plotting ;------------------------------------------------- if KEYWORD_SET(loud) then begin if N_ELEMENTS(win) eq 0 then win = 10 TVLCT, rr, gg, bb, /get PB_SET_LINE_COLOR oldpmulti = !p.multi endif ;+++++++++++++++++++++++++++++++++++++++++++++++++ ; Generate a best-overall sine wave by looking at ; only the best-fit OBS and locally smoothing them ;------------------------------------------------- chans = ['SJI_1400', 'SJI_2832'] axes = ['X', 'Y'] params = obsdata[0].paramnames cols = [255, 2, 6] for i = 0, N_ELEMENTS(chans) - 1 do begin ; Loop through img_path (1400, 2832) ; First just plot the raw (per-image) fit results for this channel WDEF, (win+3) mod 32, 1000, 700 !p.multi = [0, 1, 2] thistype = WHERE(imgdata.wave_iris eq chans[i], numtype) UTPLOT, imgdata[thistype].t_obs, imgdata[thistype].x_off, $ psym = 3, yrange = [-10, 10] + [-3,3]*i, /xstyle, chars = 1.5, $ title = chans[i] + ' X offset', /ysty, ytit = 'Pixels [X]' UTPLOT, imgdata[thistype].t_obs, imgdata[thistype].y_off, $ psym = 3, yrange = [-10, 10] + [-3,3]*i, /xstyle, chars = 1.5, $ title = chans[i] + ' Y offset', /ysty, ytit = 'Pixels [Y]' if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_aia_data.png') ; Now generate the average fit for j = 0, N_ELEMENTS(axes) - 1 do begin ; Loop through axes numtag = WHERE(obstag eq chans[i] + '_' + axes[j] + '_NUM') chisqtag = WHERE(obstag eq chans[i] + '_' + axes[j] + '_CHISQ') partag = WHERE(obstag eq chans[i] + '_' + axes[j] + '_PARAMS') goodobs = WHERE(obsdata.(numtag) ge numlim and obsdata.(chisqtag) lt chisqlim $ and obsdata.(chisqtag) gt 0, numgood) if numgood gt 0 then for k = 0, numgood - 1 do begin thisdata = CREATE_STRUCT('date_obs', obsdata[goodobs[k]].date_obs, $ 'imgpath', chans[i], $ 'axis', axes[j], $ 'num', obsdata[goodobs[k]].(numtag), $ 'chisq', obsdata[goodobs[k]].(chisqtag)) for m = 0, N_ELEMENTS(params) - 1 do begin thisdata = CREATE_STRUCT(thisdata, params[m], obsdata[goodobs[k]].(partag)[m]) endfor if N_TAGS(chandat) eq 0 then chandat = thisdata else chandat = [chandat, thisdata] endfor endfor numpoints = N_ELEMENTS(chandat) ; Generate running median (smoothed) corrections if numpoints gt 1 then begin if not KEYWORD_SET(daymedian) then daymedian = 8 taiaxis = daymedian * 86400d numdays = (MAX(taistart) - MIN(taistart)) / 86400. tais = DINDGEN(numdays)*86400 + taistart[0] days = TIME2FILE(TAI2UTC(tais), /date) tais = ANYTIM2TAI(FILE2TIME(days)) for j = 0, numdays - 1 do begin useobs = WHERE(ABS(ANYTIM2TAI(chandat.date_obs) - tais[j]) le taiaxis, numuse) if numuse gt 1 then begin usespec = chandat[useobs] thissmooth = CREATE_STRUCT('tai', tais[j], 'date', days[j], $ 'amplitude', MEDIAN(usespec.amplitude), $ 'phase', MEDIAN(usespec.phase), $ 'offset', MEDIAN(usespec.offset) ) if N_TAGS(smoothstr) eq 0 then smoothstr = thissmooth $ else smoothstr = [smoothstr, thissmooth] endif endfor endif if KEYWORD_SET(pstop) then STOP ;+++++++++++++++++++++++++++++++++++++++++++++++++ ; Generate plots ;------------------------------------------------- if KEYWORD_SET(loud) then begin ; Plot amplitude WDEF, win, 900, 600 !p.multi = 0 UTPLOT, chandat.date_obs, chandat.amplitude, psym = 4, chars = 1.5, $ ytitle = 'Amplitude [pixels]', title = chans[i] + ' Amplitude', $ yrange = [0, 3] for j = 1, N_ELEMENTS(axes) - 1 do begin rdat = chandat[WHERE(chandat.axis eq axes[j])] OUTPLOT, rdat.date_obs, rdat.amplitude, col=cols[j], psym = 4 endfor OUTPLOT, FILE2TIME(smoothstr.date), smoothstr.amplitude, thick = 2, col = 7, line = 2 if i gt 0 then begin OUTPLOT, FILE2TIME(smooths.nuv.date), smooths.nuv.amplitude, $ thick = 2, col = 6, line = 3 LEGEND, pos = 12, ['FUV trend', 'NUV trend'], col = [7, 6], chars = 1.5, linesty = [2,3] endif LEGEND, pos = 10, chars = 1.5, axes, col = cols, psym = 4 if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_wave_amplitude.png') ; Plot offset WDEF, win+1, 900, 600 pdb = IRIS_MK_POINTDB(version = pdb_ver) case i of 0 : xnom = pdb.cpx1_nuv - 1 1 : xnom = pdb.cpx1_fu1 - 1 ; 2 : xnom = pdb.cpx1_fu2 - 1 2 : xnom = 3023 endcase UTPLOT, chandat.date_obs, chandat.offset, psym = 4, chars = 1.5, $ ytitle = 'Offset [pixels]', title = chans[i] + ' average line position', $ yrange = xnom + [-5,5] OUTPLOT, FILE2TIME(smoothstr.date), smoothstr.offset, thick = 2, col = 7, line = 2 if i gt 0 then begin OUTPLOT, FILE2TIME(smooths.nuv.date), thick = 2, col = 6, line = 3, $ 0 - smooths.nuv.offset + ( MEAN(smoothstr.offset) + MEAN(smooths.nuv.offset) ) LEGEND, pos = 6, col = [7, 6], linesty = [2,3], $ ['FUV trend', 'NUV trend (sign flipped)'], chars = 1.5 endif for j = 1, N_ELEMENTS(axes) - 1 do begin rdat = chandat[WHERE(chandat.axis eq axes[j])] OUTPLOT, rdat.date_obs, rdat.offset, col=cols[j], psym = 4 endfor OUTPLOT, chandat.date_obs, chandat.offset * 0. + xnom, col = 5, line = 1 LEGEND, pos = 4, chars = 1.5, axes, col = cols, psym = 4 if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_wave_offset.png') ; Plot phase WDEF, win+2, 900, 600 UTPLOT, chandat.date_obs, (((chandat.phase/!pi)+1) mod 2) - 1, psym = 4, chars = 1.5, $ ytitle = 'Phase [pi radians]', title = chans[i] + ' Phase', yrange = [-1, 1] OUTPLOT, FILE2TIME(smoothstr.date), (((smoothstr.phase/!pi)+1) mod 2) - 1, thick = 2, col = 7, line = 2 if i gt 0 then begin OUTPLOT, FILE2TIME(smooths.nuv.date), (((smooths.nuv.phase/!pi)+1) mod 2) - 1, $ thick = 2, col = 6, line = 3 LEGEND, pos = 12, ['FUV trend', 'NUV trend'], col = [7, 6], chars = 1.5, linesty = [2,3] endif for j = 1, N_ELEMENTS(axes) - 1 do begin rdat = chandat[WHERE(chandat.axis eq axes[j])] OUTPLOT, rdat.date_obs, (((rdat.phase/!pi)+1) mod 2) - 1, col=cols[j], psym = 4 endfor LEGEND, pos = 6, chars = 1.5, axes, col = cols, psym = 4 if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_wave_phase.png') endif if i eq 0 then begin result = CREATE_STRUCT(chans[i], chandat) smooths = CREATE_STRUCT(chans[i], smoothstr) endif else begin result = CREATE_STRUCT(result, chans[i], chandat) smooths = CREATE_STRUCT(smooths, chans[i], smoothstr) endelse chandat = 0 smoothstr = 0 endfor if KEYWORD_SET(save) then begin output = CREATE_STRUCT('fitdat', result, 'smoothdat', smooths) SAVEGEN, file = CONCAT_DIR(outpath, 'wave_fit.genx'), str = output endif ;+++++++++++++++++++++++++++++++++++++++++++++++++ ; Done generating overall best-fit sine wave ;------------------------------------------------- if KEYWORD_SET(loud) then PRINT, 'IRIS_PREP_TREND_WAVE generated average sine. Elapsed time : ', $ STRING(SYSTIME(/sec) - tt0, form = '(f10.1)') if KEYWORD_SET(pstop) then STOP ;+++++++++++++++++++++++++++++++++++++++++++++++++ ; Now go back and plot the global best-fit wave ; over each obs (relies on the genx file written ; by IRIS_PREP_FILTER_WAVE, /write_logdat) ;------------------------------------------------- if KEYWORD_SET(makepngs) then begin WDEF, win, 900, 700 !p.multi = 0 bottag = WHERE(imgtags eq 'LINE_BOT') imtai = ANYTIM2TAI(imgdata.t_obs) CD, inpath, current = old_dir ff = FILE_SEARCH('*/*/*/*/wavecorr/*.txt') numf = N_ELEMENTS(ff) dates = STRMID(ff, 11, 15) years = STRMID(dates, 0, 4) months = STRMID(dates, 4, 2) days = STRMID(dates, 6, 2) tais = ANYTIM2TAI(FILE2TIME(dates)) for i = 0, numf - 2 do begin ; Loop through OBS for j = 0, N_ELEMENTS(chans)-1 do begin ; Loop through IMG_PATH usefiles = WHERE(imtai gt tais[i] and imtai lt tais[i+1] and $ imgdata.img_path eq chans[j], numuse) if numuse gt 1 then begin ; Generate a plot for this OBS and wavelength channel plotname = years[i] + '/' + months[i] + '/' + days[i] + '/' + $ dates[i] + '/wavecorr/' + dates[i] + '_' + chans[j] + '.png' usedat = imgdata[usefiles] thissum = usedat[0].sumsptrl reset_flag = 0 for k = 0, N_ELEMENTS(axes) - 1 do begin ; goodrange = WHERE(FINITE(usedat.(bottag + 3 + k)), numfin) goodrange = WHERE(usedat.(bottag + 3 + k) le chisqlim, numfin) if numfin gt 0 then begin if reset_flag eq 0 then begin flatt = usedat[goodrange].t_obs flatwave = usedat[goodrange].(bottag + k) * thissum endif else begin flatt = [flatt, usedat[goodrange].t_obs] flatwave = [flatwave, usedat[goodrange].(bottag + k) * thissum] endelse reset_flag = 1 endif endfor if N_ELEMENTS(flatt) lt 2 then goto, skipobs UTPLOT, flatt, flatwave, /ynoz, ytitle = 'Pixel location', $ chars = 1.5, psym = 4, title = plotname, $ yrange = MEDIAN(flatwave) + [-8, 8], /ystyle for k = 0, 2 do OUTPLOT, usedat.t_obs, usedat.(bottag+k), psym = 4, col = cols[k] ; Compute the smoothed sine wave for this wavelength and time interval nparams = smooths.nuv[WHERE(smooths.nuv.date eq TIME2FILE(/date, flatt[0]))] fparams = smooths.fuvs[WHERE(smooths.fuvs.date eq TIME2FILE(/date, flatt[0]))] thisphase = (nparams.phase + imgdata[usefiles[0]].roll / 180. * !pi) mod (2 * !pi) case j of 0 : begin thisoffset = nparams.offset thisamp = nparams.amplitude end 1 : begin thisoffset = fparams.offset thisamp = 0 - thisamp end 2 : begin thisoffset = 3020 thisamp = 0 - thisamp end endcase ; Look up the per-range fits and note the average pixel deviation fitind = (WHERE(obsdata.filename eq dates[i]))[0] chisqstrings = STRARR(3) for k = 0, 2 do begin thischisq = obsdata[fitind].(WHERE(obstag eq chans[j] + '_BOT_CHISQ') + k*3) if thischisq eq 0 then chisqstrings[k] = 'NaN' else $ chisqstrings[k] = STRING(thischisq, '(f6.2)') chisqstrings[k] = axes[k] + ' : ' + chisqstrings[k] endfor ; include the SINE wave fit if numuse gt 10 then begin xtai = ANYTIM2TAI(flatt) taigrid = MIN(xtai) + DINDGEN(MAX(xtai) - MIN(xtai)) ysine = thisoffset + thisamp * SIN( ((2 * !pi) / 5856d * taigrid) + thisphase ) OUTPLOT, TAI2UTC(taigrid), ysine, thick = 2 sinepred = INTERPOL(ysine, taigrid, ANYTIM2TAI(flatt)) fitdev = MEAN(ABS(flatwave - sinepred), /nan) chisqstrings = ['Fit : ' + STRING(fitdev, form = '(f6.2)'), chisqstrings] chisqpsym = [-3,4,4,4] ccol = [255, cols] clin = [0,1,1,1] for k = 0, 2 do begin thisparams = obsdata[fitind].(WHERE(obstag eq chans[j] + '_BOT_PARAMS') + k*3) ysine = thisparams[0] * SIN( ((2 * !pi) / 5856d * taigrid) + thisparams[2] ) ysine = ABS(thisparams[3]) + ysine OUTPLOT, TAI2UTC(taigrid), ysine, col = cols[k], line = 1 if KEYWORD_SET(pstop) then STOP endfor endif else begin cpsym = [-4,-4,-4] ccol = cols clin = [1,1,1] endelse LEGEND, pos = 6, chisqstrings, col = ccol, chars = 1.5, psym = cpsym PB_WIN2PNG, plotname skipobs : if KEYWORD_SET(pstop) then STOP endif endfor endfor CD, old_dir endif if KEYWORD_SET(loud) then begin TVLCT, rr, gg, bb !p.multi = oldpmulti PRINT, 'IRIS_PREP_TREND_WAVE elapsed time : ', STRING(SYSTIME(/sec) - tt0, form = '(f10.1)') endif if KEYWORD_SET(pstop) then STOP end