;
;	plot_see.pro
;
;	plot TIMED SEE L3 time series or spectrum
;
;	USAGE: plot_see, 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, planet=planet
;	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
;
;       planet = planet name or number (1-9) to time shift SEE data to planet
;
;		data = option to return data that was plotted
;			the "data" can be passed back to plot_see as PROXY data
;
;
;	HISTORY:
;	2/1/03		Tom Woods	Original Code - based on plot_sol.pro
;   8/11/04     Karen Turk  Modified for IDL version 5.0.3
;   5/26/05     Tom Woods   Updated with "planet" option
;   10/25/07    Don Woodraska Fixed /undo1au option, 1au is interpolated for 
;                             missing values
;
;idver='$Id: plot_see.pro,v 9.2 2008/10/20 14:57:18 see_sw Exp $'
;

pro plot_see, wave_or_day, type=type, xrange=xrange, yrange=yrange, file=file, $
	proxy=proxy, ptitle=ptitle, pcolor=pcolor, list_lines=list_lines, $
	line_num=line_num, xps_num=xps_num, data=data, undo1au=undo1au, $
	planet=planet, debug=debug
; , 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_see, wave_or_day, type=type, file=file,  $'
  print, '                 line_num=line_num, /list_lines, xps_num=xps_num, $'
  print, '                 xrange=xrange, yrange=yrange, data=data, /undo1au, $'
  print, '                 proxy=proxy, ptitle=ptitle, pcolor=pcolor, $'
  print, '                 planet=planet'
  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, '    planet      = planet name or number (1-9) to time shift SEE data to planet'
  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 L3 data in common block
;	so this procedure doesn't have to read in data for every call
;
common  plot_see_data, see, fname, windowset
ssee = size(see)
if (ssee[0] ne 1) and (ssee[2] ne 8) then begin
read_file:
	windowset = 0
	if keyword_set(file) then fname = file $
	else fname = 'see__L3_merged_*.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_see 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 L3 data from ', fname
	read_netcdf, fname, see, seeatt, status
    ;interpolate days with bad cor_1au values
    bad_cor1au = where(see.cor_1au lt .5 or see.cor_1au gt 1.5,comp=comp,$
                       n_bad_cor1au)
    if n_bad_cor1au gt 0 then begin
        cor1au=see.cor_1au
        jd=yd2jd(see.date)
        cor1au[bad_cor1au]=interpol(cor1au[comp],jd[comp],jd[bad_cor1au])
        see.cor_1au = cor1au
    endif
	print, 'The "see" 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

;
;	keep Planet information in common block 
;	so don't have to re-calculate planet position on every call
;	unless the planet changes
;
common  plot_see_planet_info, lastplanet, planetdist2, planetshift, planetstr
psize = size(lastplanet)
if (psize[psize[0]+1] eq 0) then begin
  lastplanet=-1
  planetstr = '?'
endif
if keyword_set(planet) then begin
  planetnum = 0L
  psize = size(planet)
  if (psize[1] eq 7) then begin
	  ;  Figure out which Planet string was given
	  pstr = strupcase(strmid(planet,0,3))
	  pall = [ "SUN", "MER", "VEN", "EAR", "MAR", "JUP", "SAT", "URA", "NEP", "PLU" ]
	  planetnum = where( pall eq pstr )
  endif else planetnum = planet
  if (planetnum ne lastplanet) then begin
    print, 'plot_see: calculating planet position...'
    if (planetnum lt 1) or (planetnum gt 9) then planetnum = 3	; Earth
    lastplanet = planetnum
    case planetnum of
      1: planetstr = 'Mercury'
      2: planetstr = 'Venus'
      3: planetstr = 'Earth'
      4: planetstr = 'Mars'
      5: planetstr = 'Jupiter'
      6: planetstr = 'Saturn'
      7: planetstr = 'Uranus'
      8: planetstr = 'Neptune'
      9: planetstr = 'Pluto'
      else: planetstr = '?'
    endcase
    ndates = n_elements(see.date)
    planetdist2 = fltarr(ndates)
    planetshift = fltarr(ndates,2,2)
    for k=0,ndates-1 do begin
      earth2planet, planetnum, yd2jd(see.date[k]), pdist, pshift
      ; division correction is R^2
      planetdist2[k] = pdist^2.    
      ;  force time shifts to be integers
      pshift[0,0] = round(pshift[0,0])
      pshift[1,0] = round(pshift[1,0])
      planetshift[k,*,*] = pshift
    endfor
    ;
    ;  fix end points so shifting will not crash program
    ;
    for k=0,27 do begin
      if (planetshift[k,0,0] lt (-1*k)) then planetshift[k,0,0] = -1*k
      if (planetshift[ndates-1-k,0,0] gt k) then planetshift[ndates-1-k,0,0] = k
      if (planetshift[k,1,0] lt (-1*k)) then planetshift[k,1,0] = -1*k
      if (planetshift[ndates-1-k,1,0] gt k) then planetshift[ndates-1-k,1,0] = k
    endfor
  endif
endif

;
;	List the lines (name, wavelength) if requested
;
if keyword_set(list_lines) then begin
  num_lines = n_elements(see[0].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[0].linename[k-1], see[0].linewave[k-1], $
      k+num2, see[0].linename[k+num2-1], see[0].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[0].linename[k-1], see[0].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

version = long( string(see.ver.sdp) ) ;8/11/04 KBT added string() to correct bug in IDL version 5.0.3

;
;	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'

mtitle = 'TIMED SEE'
if keyword_set(planet) then begin
  mtitle = planetstr + ' - ' + mtitle
endif

if (ptype eq 'SP') then begin
	;
	;	Plot SPECTRUM
	;
    if (windowset eq 0) then begin
        window, retain=2
        windowset = 1
	endif
    temp=min( abs(see.date - wave_or_day[0]), wgood )
    theDay = see.date[wgood]
    dayStr = strtrim( long(theDay/1000L), 2 )+'/'+strtrim(long(theDay mod 1000L),2)
    factor = 1.0
    if keyword_set(planet) then begin
        factor = 1. / planetdist2[wgood]
	    ;  do weighted average for planet time shifting
        spflux = see.sp_flux[*,wgood+long(planetshift[wgood,0,0])] * planetshift[wgood,0,1] $
                 + see.sp_flux[*,wgood+long(planetshift[wgood,1,0])] * planetshift[wgood,1,1]
        spflux = reform(spflux) * 1000. * factor
        if (version lt 7) then begin
            spfluxfill = see.sp_flux_filled[*,wgood+long(planetshift[wgood,0,0])] * planetshift[wgood,0,1] $
                         + see.sp_flux_filled[*,wgood+long(planetshift[wgood,1,0])] * planetshift[wgood,1,1]
            spfluxfill = reform(spfluxfill) * 1000. * factor
        endif
        splines = see.line_flux[*,wgood+long(planetshift[wgood,0,0])] * planetshift[wgood,0,1] $
                  + see.line_flux[*,wgood+long(planetshift[wgood,1,0])] * planetshift[wgood,1,1]
        splines = splines * 1000. * factor
    endif else begin
        if keyword_set(undo1au) then factor = see.cor_1au[wgood]
        spflux = reform(see.sp_flux[*,wgood]) * 1000. * factor
        if (version lt 7) then spfluxfill = reform(see.sp_flux_filled[*,wgood]) * 1000. * factor
        splines = see.line_flux[*,wgood] * 1000. * factor
    endelse
    if (n_elements(wave_or_day) gt 1) then begin
	    ;  over-plot more days and include average into "data"
        temp=min( abs(see.date - wave_or_day[1]), wgood2 )
        theDay2 = see.date[wgood2]
        dayStr = dayStr + ' - ' + $
                 strtrim( long(theDay2/1000L), 2 )+'/'+strtrim(long(theDay2 mod 1000L),2)
        if (wgood2 lt wgood) then kstep = -1L else kstep = 1L

        cnt = fltarr( n_elements(spflux) )
        wgd = where( spflux gt 0 )
        if (wgd[0] ne -1) then cnt[wgd] = 1
        if (version lt 7) then begin
            cntfill = fltarr( n_elements(spfluxfill) )
            wgd = where( spfluxfill gt 0 )
            if (wgd[0] ne -1) then cntfill[wgd] = 1
        endif
        cntline = fltarr( n_elements(splines) )
	    wgd = where( splines 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(planet) then begin
                    factor = 1. / planetdist2[k]
				    ;  do weighted average for planet time shifting
                    spflux2 = see.sp_flux[*,k+long(planetshift[k,0,0])] * planetshift[k,0,1] $
                              + see.sp_flux[*,k+long(planetshift[k,1,0])] * planetshift[k,1,1]
                    spflux2 = reform(spflux2) * 1000. * factor
                    if (version lt 7) then begin
                        spfluxfill2 = see.sp_flux_filled[*,k+long(planetshift[k,0,0])] * planetshift[k,0,1] $
                                      + see.sp_flux_filled[*,k+long(planetshift[k,1,0])] * planetshift[k,1,1]
                        spfluxfill2 = reform(spfluxfill2) * 1000. * factor
                    endif
                    splines2 = see.line_flux[*,k+long(planetshift[k,0,0])] * planetshift[k,0,1] $
                               + see.line_flux[*,k+long(planetshift[k,1,0])] * planetshift[k,1,1]
                    splines2 = splines2 * 1000. * factor
                endif else begin
                    if keyword_set(undo1au) then factor = see.cor_1au[k]
                    spflux2 = reform(see.sp_flux[*,k]) * 1000. * factor
                    if (version lt 7) then spfluxfill2 = reform(see.sp_flux_filled[*,k]) * 1000. * factor
                    splines2 = see.line_flux[*,k] * 1000. * factor
                endelse
                wgd = where( spflux2 gt 0 )
                if (wgd[0] ne -1) then begin
                    spflux[wgd] = spflux[wgd] + spflux2[wgd]
                    cnt[wgd] = cnt[wgd] + 1
                endif
                if (version lt 7) then begin
                    wgd = where( spfluxfill2 gt 0 )
                    if (wgd[0] ne -1) then begin
                        spfluxfill[wgd] = spfluxfill[wgd] + spfluxfill2[wgd]
                        cntfill[wgd] = cntfill[wgd] + 1
                    endif
                endif
                wgd = where( splines2 gt 0 )
                if (wgd[0] ne -1) then begin
                    splines[wgd] = splines[wgd] + splines2[wgd]
                    cntline[wgd] = cntline[wgd] + 1
                endif
            endfor
            wgd = where( cnt gt 0 )
            if (wgd[0] ne -1) then spflux[wgd] = spflux[wgd] / cnt[wgd]
            if (version lt 7) then begin
                wgd = where( cntfill gt 0 )
                if (wgd[0] ne -1) then spfluxfill[wgd] = spfluxfill[wgd] / cntfill[wgd]
            endif
            wgd = where( cntline gt 0 )
            if (wgd[0] ne -1) then splines[wgd] = splines[wgd] / cntline[wgd]
        endif
    endif
    theflux = spflux
    if (version lt 7) then theflux = spfluxfill
	if keyword_set(xrange) then theXrange=xrange else theXrange=[0,200]
    if keyword_set(yrange) then theYrange=yrange else begin
        theYmax = 10.^round(alog10(max(theflux)))
        if (theYmax lt max(theflux)) then theYmax = theYmax * 3.
        theYrange = [theYmax*1E-4, theYmax]
    endelse
	plot_io, see.sp_wave, theflux, xrange=theXrange, xstyle=1, $
             yrange=theYrange, ystyle=1, xtitle='Wavelength [nm]', $
             ytitle='Irradiance [mW/m!U2!N/nm]', $
             title=mtitle+': '+dayStr, thick=2, charthick=1.5
    if (version lt 7) then oplot, see.sp_wave, spflux, line=2
	if keyword_set(line_num) then begin
        oplot, see.linewave, splines, psym=4
		data = dblarr(2,n_elements(splines))
		data[0,*] = see.linewave
		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,*] = theflux[wwv]
		endif
	endelse
endif else begin
	;
	;	Time Series Plot OR Comparison Plot with PROXY data
	;	Both have Time Series Plot
	;
    if (windowset eq 0) then begin
        window, retain=2
        windowset = 1
    endif
    if (ptype eq 'CP') then !p.multi=[0,1,2]
	;
	;	set "theTime" array as fraction of year
	; 		leap year is either (divisible by 4 and not divisible by 100) or divisible by 400
	;
    yrdays = 365.D0
	year = long(see.date / 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 + (long((see.date-1) mod 1000L)+0.5D0) / (yrdays + yrisleap)
	;
	;	first plot (top) is TIME SERIES
	;
	factor = 1.0
	if keyword_set(planet) then begin
        factor = 1. / planetdist2
	endif else begin
        if keyword_set(undo1au) then factor = see.cor_1au
	endelse
    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
        if keyword_set(planet) then begin
            wpos = where(see.line_flux[i_param,*] gt 0)
            allFlux = interpol(see.line_flux[i_param,wpos],theTime[wpos], $
                               theTime+long(planetshift[*,0,0])/yrdays) * planetshift[*,0,1] $
                      + interpol(see.line_flux[i_param,wpos],theTime[wpos], $
                                 theTime+long(planetshift[*,1,0])/yrdays) * planetshift[*,1,1]
        endif else begin
            allFlux = see.line_flux[i_param,*]	    
        endelse
        allFlux = allFlux * 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
        if keyword_set(planet) then begin
            wpos = where(see.xp_flux[i_param,*] gt 0)
            allFlux = interpol(see.xp_flux[i_param,wpos], theTime[wpos], $
                               theTime+long(planetshift[*,0,0])/yrdays) * planetshift[*,0,1] $
                      + interpol(see.xp_flux[i_param,wpos], theTime[wpos], $
                                 theTime+long(planetshift[*,1,0])/yrdays) * planetshift[*,1,1]
        endif else begin
            allFlux = see.xp_flux[i_param,*]	    
	    endelse
	    allFlux = allFlux * 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)' )
	    if (version lt 7) then begin
            if keyword_set(planet) then begin
                wpos = where( see.sp_flux_filled[wgood,*] gt 0 )
                allFlux = interpol(see.sp_flux_filled[wgood,wpos], theTime[wpos], $
                                   theTime+long(planetshift[*,0,0])/yrdays) * planetshift[*,0,1] $
                          + interpol(see.sp_flux_filled[wgood,wpos], theTime[wpos], $
                                     theTIme+long(planetshift[*,1,0])/yrdays) * planetshift[*,1,1]
            endif else begin
                allFlux = see.sp_flux_filled[wgood,*]	    
            endelse
 	    endif else begin
            if keyword_set(planet) then begin
                wpos = where( see.sp_flux[wgood,*] gt 0 )
                allFlux = interpol(see.sp_flux[wgood,wpos], theTime[wpos], $
                                   theTime+long(planetshift[*,0,0])/yrdays) * planetshift[*,0,1] $
                          + interpol(see.sp_flux[wgood,wpos], theTime[wpos], $
                                     theTime+long(planetshift[*,1,0])/yrdays) * planetshift[*,1,1]
            endif else begin
                allFlux = see.sp_flux[wgood,*]	    
            endelse
        endelse
        allFlux = reform(allFlux) * 1000. * factor
        theYtitle = 'Irradiance [mW/m!U2!N/nm]'
        theTitle = mtitle+' 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(see.sp_wave) )
                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	           
                   if (version lt 7) then begin
                       if keyword_set(planet) then begin
                           wpos = where( see.sp_flux_filled[k,*] gt 0 )
                           tempFlux = interpol(see.sp_flux_filled[k,wpos], theTime[wpos], $
                                               theTime+long(planetshift[*,0,0])/yrdays) * planetshift[*,0,1] $
                                      + interpol(see.sp_flux_filled[k,wpos], theTime[wpos], $
                                                 theTime+long(planetshift[*,1,0])/yrdays) * planetshift[*,1,1]
                       endif else begin
                           tempFlux = see.sp_flux_filled[k,*]	    
                       endelse
                   endif else begin
                       if keyword_set(planet) then begin
                           wpos = where( see.sp_flux[k,*] gt 0 )
                           tempFlux = interpol(see.sp_flux[k,wpos], theTime[wpos], $
                                               theTime+long(planetshift[*,0,0])/yrdays) * planetshift[*,0,1] $
                                      + interpol(see.sp_flux[k,wpos], theTime[wpos], $
                                                 theTime+long(planetshift[*,1,0])/yrdays) * planetshift[*,1,1]
                       endif else begin
                           tempFlux = see.sp_flux[k,*]	    
                       endelse
                   endelse
                   tempFlux = reform(tempFlux) * 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 = mtitle+': '+waveStr+'-'+waveStr2+' nm'
        endif
    endelse
	;  FIX bad data
    wbad = finite(allFlux)
	wfix = where(wbad eq 0)
	if (wfix[0] ne -1) then allFlux[wfix] = 0.0
	wthere = where( allFlux gt 0. )
    if (wthere[0] eq -1) then begin
        print, 'WARNING: No valid data found.'
        return
	endif
	if keyword_set(xrange) then theXrange=xrange else theXrange=[2002,long(2.*(max(theTime)+0.5))/2.]
	if keyword_set(yrange) then theYrange=yrange else begin
        theYrange = [ min(allFlux[wthere])*errormult[0], $
                      max(allFlux[wthere])*errormult[1] ]
	endelse
	plot, theTime[wthere], allFlux[wthere],  xrange=theXrange, xstyle=1, $
          ystyle=1, thick=2, xtitle='Time [year]', ytitle=theYtitle, $
          title=theTitle, xmargin=[12,3], yrange=theYrange
	; 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
	    
          plot, proxyFlux[wthere2], allFlux[wthere2], psym=3, thick=3, $
                xstyle=1, ystyle=1, xtitle=ptitle, $
                ytitle=theYtitle, title='', xmargin=[12,3]
	    ; 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

if keyword_set(debug) then stop, 'Check out results at end of plot_see...'

return
end