; ; plot_see3a.pro ; ; plot TIMED SEE L3A time series or spectrum ; L3A is orbit average, for daily average use L3 data and plot_see.pro ; ; USAGE: plot_see3A, wave_or_day, type=type, xrange=xrange, file=file, ; proxy=proxy, ptitle=ptitle, line_num=line_num, xps_num=xps_num, ; data=data, list_lines=list_lines, undo1au=undo1au ; where ; wave_or_day = wavelength (nm) if type='TS' or 'CMP' ; = UARS Day if type='SP' ; ; type = 'TS' for Time Series (default if not given) ; = 'SP' for SPectrum ; = 'CP' for Comparison Plot between SEE and PROXY data ; ; xrange = [ xmin, xmax ] for setting plot's X range ; X range is time for Time Series (TS) ; X range is wavelength for Spectrum (SP) ; X range is time range for Comparison Plot (CP) ; ; file = NAME OF TIMED SEE L3 SAVE SET (not required) ; (IDL save set made using make_see.pro) ; ; line_num = 1-38 for plotting a LINE TS instead of ; wavelength (ignores wave_or_day value) ; ; list_lines = option to list the lines (name, wavelength) ; ; xps_num = 1-12 (exclude 4,8,12) for plotting a XPS TS ; instead of wavelength (ignores wave_or_day value) ; ; proxy = fltarr(2, *) of proxy data (or other SEE data) ; proxy[0,*] = time - fractional year ; proxy[1,*] = proxy data ; ONLY valid if use the type = 'CP' ; ; ptitle= title for Proxy plots ; ; data = option to return data that was plotted ; the "data" can be passed back to plot_see as PROXY data ; ; EXAMPLES: ; IDL> plot_see3a,30 ;plots a time series of the 30.5 nm bin ; IDL> plot_see3a,/list_lines ;lists extracted lines and numbers ; IDL> plot_see3a,line_num=1 ;plots a time series of the extracted ; He II 30.4 nm line ; IDL> plot_see3a,xps_num=1 ;plots a time series of XPS diode #1 ; ; IDL> plot_see3a,30,data=data ;plots a time series of the 30.5 nm bin ; and returns the data in "data" ; IDL> plot_see3a,line_num=1,proxy=data ;plots the time series and ; plots the correlation of 30.5 with the proxy ; ; HISTORY: ; 10/20/03 Tom Woods Original Code - based on plot_see.pro ; 05/10/04 Don Woodraska Modified plot_see3a.pro ; 02/19/07 Don Woodraska Adde Tom's div by zero fix ; ;idver='$Id: plot_see3a.pro,v 8.3 2007/02/19 19:24:02 dlwoodra Exp $' ; pro plot_see3a, wave_or_day, type=type, xrange=xrange, file=file, $ proxy=proxy, ptitle=ptitle, pcolor=pcolor, list_lines=list_lines, $ line_num=line_num, xps_num=xps_num, data=data, undo1au=undo1au, $ yrange=yrange, ylog=ylog ; , errorplot=errorplot ; ; give user HELP if called without parameters ; if (n_params(0) lt 1) and (not keyword_set(line_num)) $ and (not keyword_set(xps_num)) then begin print, ' ' print, 'USAGE: plot_see3a, wave_or_day, type=type, file=file, $' print, ' line_num=line_num, /list_lines, xps_num=xps_num, $' print, ' xrange=xrange, data=data, /undo1au, $' print, ' proxy=proxy, ptitle=ptitle, pcolor=pcolor' print, ' ' print, " wave_or_day = wavelength if type='TS' or 'CP'" print, ' (can be wavelength range)' print, " = Day (YYYYDOY) if type='SP'" print, ' (can be day range)' print, ' ' print, ' /undo1au = undo correction to 1 AU (put irradiance at Earth)' print, ' ' print, ' line_num = 1-38 are valid values for L3 LINEs (emission lines)' print, " or /line_num is option to plot the lines if type='SP'" print, ' /list_lines = option to list the lines (name, wavelength)' print, ' ' print, ' xps_num = 1-12 (exclude 4,8,12) for plotting XPS diode data' print, ' ' print, " proxy = fltarr(2,*) of proxy data used for type='CP'" print, ' ptitle = title for PROXY plots' print, ' pcolor = color for PROXY plots or color table for multiple SP plots' print, ' ' print, ' file = specify the file for the SEE L3 merged file' print, ' (normally not needed as it will search for latest version)' print, ' ' print, " data = output for this plot, fltarr(2,*) for 'TS' or 'SP'" print, " or fltarr(3,*) for 'CP' type" print, ' ' wave_or_day = 121.5 ; use default value to make some kind of plot endif ; ; keep TIMED SEE L3A data in common block ; so this procedure doesn't have to read in data for every call ; common plot_see3a_data, see, fname ssee = size(see) if (ssee[0] ne 1) and (ssee[2] ne 8) then begin read_file: if keyword_set(file) then fname = file $ else fname = 'see__L3A_merged_2*.ncdf' if (float(!version.release) gt 6.0) then $ filelist = file_search( fname, count=fcount ) $ else filelist = findfile( fname, count=fcount ) if (fcount lt 1) then begin fdir = getenv( 'see_data_merged' ) if (strlen(fdir) gt 0) then begin fsearch = fdir + '/' + fname if (float(!version.release) gt 6.0) then $ filelist = file_search( fsearch, count=fcount ) $ else filelist = findfile( fsearch, count=fcount ) endif else fcount = 0 if (fcount lt 1) then begin fname = '' print, 'ERROR: plot_see3a can not find data file - ', fname return endif endif fname = filelist[fcount-1] ; pick latest version of file file = fname print, 'Reading in TIMED SEE L3A data from ', fname read_netcdf, fname, see, seeatt, status print, 'The "SEE L3A" data structure is:' help, see, /str print, ' ' endif ; ; check if need to re-read SEE data if "file" changes ; if keyword_set(file) then begin if (file ne fname) then goto, read_file endif ; ; List the lines (name, wavelength) if requested ; if keyword_set(list_lines) then begin num_lines = n_elements(see.linewave) print, '********************* SEE L3 Line List *********************' print, ' ID Emission Wave (nm) ID Emission Wave (nm)' print, ' -- -------- --------- -- -------- ---------' num2 = num_lines/2 for k=1,num2 do $ print, k, see.linename[k-1], see.linewave[k-1], $ k+num2, see.linename[k+num2-1], see.linewave[k+num2-1], $ format='(I3,A12,F10.1,I12,A12,F10.1)' k = num_lines if ( num_lines gt (num2 * 2)) then $ print, k, see.linename[k-1], see.linewave[k-1], $ format='(I3,A12,F10.1)' print, ' ' if (n_params() lt 1) and (not keyword_set(line_num)) $ and (not keyword_set(xps_num)) then return endif ; ; common stuff for making the plots ; !fancy=4 !psym=10 data = 0.0 ; empty data set errormult = [0.9, 1.1] ; 10% typical error ; ; now do plot based on "type" ; ptype = 'TS' if keyword_set(type) then ptype = strupcase(type) if keyword_set(proxy) then ptype = 'CP' if (not keyword_set(proxy)) and (ptype eq 'CP') then ptype='TS' ; make date a fraction of day ;seedate = double(see.date) + (see.start_time + see.stop_time)/2.D0/24./3600. seedate = double(see.date) + (see.time)/86400.d0 if (ptype eq 'SP') then begin ; ; Plot SPECTRUM ; temp=min( abs(seedate - wave_or_day[0]), wgood ) theDay = seedate[wgood] dayStr = strtrim( long(theDay/1000L), 2 )+'/'+string((theDay mod 1000L),format='(F5.1)') factor = 1.0 if keyword_set(undo1au) then factor = see.cor_1au[wgood] spflux = reform(see.sp_flux[*,wgood]) * 1000. * factor spfluxfill = reform(see.sp_flux[*,wgood]) * 1000. * factor if keyword_set(xrange) then theXrange=xrange else theXrange=[0,200] splines = see.line_flux[*,wgood] * 1000. * factor if (n_elements(wave_or_day) gt 1) then begin ; over-plot more days and include average into "data" temp=min( abs(seedate - wave_or_day[1]), wgood2 ) theDay2 = seedate[wgood2] dayStr = dayStr + ' - ' + $ strtrim( long(theDay2/1000L), 2 )+'/'+string((theDay2 mod 1000L),format='(F5.1)') if (wgood2 lt wgood) then kstep = -1L else kstep = 1L cnt = fltarr( n_elements(see.sp_flux[*,0]) ) wgd = where( see.sp_flux[*,wgood] gt 0 ) if (wgd[0] ne -1) then cnt[wgd] = 1 cntfill = fltarr( n_elements(see.sp_flux[*,0]) ) wgd = where( see.sp_flux[*,wgood] gt 0 ) if (wgd[0] ne -1) then cntfill[wgd] = 1 cntline = fltarr( n_elements(see.line_flux[*,0]) ) wgd = where( see.line_flux[*,wgood] gt 0 ) if (wgd[0] ne -1) then cntline[wgd] = 1 if abs(wgood2-wgood) gt 0 then begin for k=wgood+kstep,wgood2,kstep do begin factor = 1.0 if keyword_set(undo1au) then factor = see.cor_1au[k] wgd = where( see.sp_flux[*,k] gt 0 ) if (wgd[0] ne -1) then begin spflux[wgd] = spflux[wgd] + $ (reform(see.sp_flux[wgd,k]) * 1000. * factor) cnt[wgd] = cnt[wgd] + 1 endif wgd = where( see.sp_flux[*,k] gt 0 ) if (wgd[0] ne -1) then begin spfluxfill[wgd] = spfluxfill[wgd] + $ (reform(see.sp_flux[wgd,k]) * 1000. * factor) cntfill[wgd] = cntfill[wgd] + 1 endif wgd = where( see.line_flux[*,k] gt 0 ) if (wgd[0] ne -1) then begin splines[wgd] = splines[wgd] + $ (reform(see.line_flux[wgd,k]) * 1000. * factor) cntline[wgd] = cntline[wgd] + 1 endif endfor wgd = where( cnt gt 0 ) if (wgd[0] ne -1) then spflux[wgd] = spflux[wgd] / cnt[wgd] wgd = where( cntfill gt 0 ) if (wgd[0] ne -1) then spfluxfill[wgd] = spfluxfill[wgd] / cntfill[wgd] wgd = where( cntline gt 0 ) if (wgd[0] ne -1) then splines[wgd] = splines[wgd] / cntline[wgd] endif endif if size(yrange,/type) eq 0 then yr=[1e-3,3e1] else yr=yrange plot_io, see.sp_wave, spfluxfill, xrange=theXrange, xstyle=1, $ yrange=yr, ystyle=1, xtitle='Wavelength [nm]', $ ytitle='Irradiance [mW/m!U2!N/nm]', $ title='TIMED SEE: '+dayStr, thick=2, charthick=1.5 oplot, see.sp_wave, spflux, line=2 if keyword_set(line_num) then begin oplot, see.line_wave, splines, psym=4 data = dblarr(2,n_elements(splines)) data[0,*] = see.line_wave data[1,*] = splines endif else begin wwv = where( (see.sp_wave ge theXrange[0]) and (see.sp_wave le theXrange[1])) if (wwv[0]) ne -1 then begin data = dblarr(2,n_elements(wwv)) data[0,*] = see.sp_wave[wwv] data[1,*] = spfluxfill[wwv] endif endelse endif else begin ; ; Time Series Plot OR Comparison Plot with PROXY data ; Both have Time Series Plot ; if (ptype eq 'CP') then !p.multi=[0,1,2] ; ; first plot (top) is TIME SERIES ; factor = 1.0 if keyword_set(undo1au) then factor = see.cor_1au if keyword_set(line_num) then begin i_param = line_num[0] if (i_param lt 1) then i_param=1 if (i_param gt n_elements(see.linewave)) then i_param=n_elements(see.linewave) i_param = i_param - 1 ; index is zero based allFlux = see.line_flux[i_param,*]*1000.*factor theYtitle = 'Irr [mW/m!U2!N]' theTitle = 'EGS ' + see.linename[i_param] + $ string( see.linewave[i_param], format='(F6.1)') + ' nm' endif else if keyword_set(xps_num) then begin xpsvalid = [1,2,3,5,6,7,9,10,11] i_param = where( xpsvalid eq xps_num[0] ) i_param = i_param[0] if (i_param lt 0) then i_param=0 allFlux = see.xp_flux[i_param,*] * 1000. * factor theYtitle = 'XP#'+strtrim(xpsvalid[i_param],2) + ' Irr [mW/m!U2!N]' theTitle = 'XPS '+strtrim(long(see.xpwave1[i_param]),2) + '-' + $ strtrim(long(see.xpwave2[i_param]),2) + ' nm' endif else begin theWave = long(wave_or_day[0]) + 0.5 temp=min( abs(see.sp_wave - theWave), wgood ) theWave = see.sp_wave[wgood] waveStr = string( theWave, format='(F5.1)' ) allFlux = reform(see.sp_flux[wgood,*]) * 1000. * factor theYtitle = 'Irradiance [mW/m!U2!N/nm]' theTitle = 'TIMED SEE at '+waveStr+' nm' if (n_elements(wave_or_day) gt 1) then begin ; integrate over more wavelengths theWave2 = long(wave_or_day[1]) - 0.5 temp=min( abs(see.sp_wave - theWave2), wgood2 ) theWave2 = see.sp_wave[wgood2,*] waveStr = strtrim( long(theWave-0.5), 2 ) waveStr2 = strtrim( long(theWave2+0.5), 2 ) if (wgood2 lt wgood) then kstep = -1L else kstep = 1L if abs(wgood2-wgood) gt 0 then begin cnt = fltarr( n_elements(allFlux) ) wgd = where( allFlux gt 0. ) if (wgd[0] ne -1) then cnt[wgd] = 1 cnttotal = abs(wgood2-wgood) + 1 for k=wgood+kstep,wgood2,kstep do begin tempFlux = reform(see.sp_flux[k,*]) * 1000. * factor wgd = where( tempFlux gt 0.) if (wgd[0] ne -1) then begin allFlux[wgd] = allFlux[wgd] + tempFlux[wgd] cnt[wgd] = cnt[wgd] + 1 endif endfor wgd = where( cnt gt 0 ) allFlux[wgd] = allFlux[wgd] * cnttotal / cnt[wgd] endif theTitle = 'TIMED SEE: '+waveStr+'-'+waveStr2+' nm' endif endelse wthere = where( allFlux gt 0. ) if (wthere[0] eq -1) then begin print, 'WARNING: No valid data found.' return endif yrdays = 365.D0 year = long(seedate) / 1000L yrisleap = lonarr(n_elements(year)) wleap=where( ( ((year mod 4) eq 0) and ((year mod 100) ne 0) ) or ((year mod 400) eq 0), n_true) if n_true gt 0 then yrisleap[wleap] = 1 theTime = year + ((seedate-1) mod 1000L) / (yrdays + yrisleap) if keyword_set(xrange) then theXrange=xrange else theXrange=[2002,long(2.*(max(theTime)+0.5))/2.] theYrange = [ min(allFlux[wthere])*errormult[0], $ max(allFlux[wthere])*errormult[1] ] if size(yrange,/type) eq 0 then yr=theYrange else yr=yrange plot, theTime[wthere], allFlux[wthere], xrange=theXrange, xstyle=1, $ ystyle=1, thick=2, xtitle='Time [year]', ytitle=theYtitle, $ title=theTitle, xmargin=[12,3], yrange=yr,ylog=ylog ; if keyword_set(errorplot) then begin ; t1 = 1996.5 ; dt = 0.25 ; temp = min( abs(sol.year[wthere] - t1), wpt ) ; y1 = (smooth( allflux[wthere], 81 ))[wpt] ; oplot, [t1-dt, t1+dt, t1, t1, t1-dt, t1+dt], $ ; [ [y1,y1,y1]*errormult[0], [y1,y1,y1]*errormult[1] ], $ ; thick=3 ; endif ; ; second plot (bottom) is comparison to PROXY data ; if (ptype ne 'CP') then begin data = dblarr(2,n_elements(wthere)) data[0,*] = theTime[wthere] data[1,*] = allFlux[wthere] endif else begin if (not keyword_set(ptitle)) then ptitle='Proxy Data' proxyFlux = interpol( proxy[1,*], proxy[0,*], theTime ) wthere2 = where( (allFlux gt 0.) and (proxyFlux gt 0.) $ and (theTime ge theXrange[0]) and (theTime le theXrange[1]) ) if (wthere2[0]) ne -1 then begin ; overplot proxy data first cfit = poly_fit( proxyFlux[wthere2], allFlux[wthere2], 1, yfit = pfit ) if (keyword_set(pcolor)) then oplot, theTime[wthere2], pfit, color=pcolor $ else oplot, theTime[wthere2], pfit, line=2 rcor = correlate( proxyFlux[wthere2], allFlux[wthere2] ) plot, proxyFlux[wthere2], allFlux[wthere2], psym=3, thick=3, $ xstyle=1, ystyle=1, xtitle=ptitle, xmargin=[12,3], $ ytitle=theYtitle, title='R = '+string(rcor,format="(F5.2)") ; if (keyword_set(pcolor)) then oplot,proxyFlux[wthere2], allFlux[wthere2], $ ; psym=3, color=pcolor data = dblarr(3,n_elements(wthere2)) data[0,*] = theTime[wthere2] data[1,*] = allFlux[wthere2] data[2,*] = proxyFlux[wthere2] endif !p.multi=0 endelse endelse return end