;
;	 plot of TIMED SEE XPS Level4 time series
;
;	Tom Woods
;	July 12, 2006
;
;XPS Level 4 Structure for TIMED XPS and SORCE XPS (shared code)
; ---------------------
;   DATE            LONG           YYYYDOY
;   TIME            FLOAT          Seconds of Day
;   XPS_QS          FLOAT          QS scale factor (usually 1.0)
;   XPS_AR          FLOAT          AR scale factor (usually < 0.1)
;   XPS_FLARE       FLOAT          Flare scale factor
;   GOES_QS         FLOAT          GOES QS scale factor
;   GOES_AR         FLOAT          GOES AR scale factor
;   GOES_FLARE      FLOAT          GOES Flare scale factor
;   FMTEMP          FLOAT          Flare Model temperature
;   FMINDEX         INT            Flare Model index into temperature array
;   FMWEIGHT        FLOAT          Flare Model weight at "index"
;   ERR_ABS         FLOAT          Error Absolute (accuracy)
;   ERR_MEAS        FLOAT          Error Measurement (precision)
;   MODELFLUX       FLOAT Array[400]  Irradiance in W/m^2/nm from 0.05 to 39.95 nm
;
;
;	optmized to use with merged L4A data set which has only 3-min orbit averages
;	Run make_xps4_merged.pro (if necessary) to create this file
;
;	Changed 1/31/07 (TNW) to have MINUTE, ORBIT, and DAY averaged merged data sets
;   Changed 4/3/07 (TNW) to have FILE input option
;   Changed 4/20/07 (DLW) Added clearcommon keyword to clear all common blocks
;   Changed 12/10/07 (TNW) to only select good (positive) irradiance values
;
;$Id: plot_xps4.pro,v 8.3 2007/04/21 00:17:09 dlwoodra Exp $
;
pro plot_xps4, waverange, data, log=log, $
	minute=minute, orbit=orbit, day=day, file=file, $
	xrange=xrange, yrange=yrange, ignoreflare=ignoreflare, $
	ddata=ddata, f107plot=f107plot, $
	nodayplot=nodayplot, debug=debug, clearcommon=clearcommon

if (n_params() lt 1) then begin
  print, 'USAGE: plot_xps4, waverange, data, file=file, /minute, /orbit, /day, '
  print, '                    /log, xrange=xrange, yrange=yrange, '
  print, '                    /ignoreflare, /f107plot, ddata=ddata'
  print, ' '
  print, '    waverange     Wavelength range [in nm] for plot: 0-40 nm valid range'
  print, '    data          output data array: data[0,*] is time, [1,*] is irradiance'
  print, ' '
  print, '    /day          Option for DAY average file in $see_data_merged'
  print, '    /orbit        Option for ORBIT average file in $see_data_merged [DEFAULT]'
  print, '    /minute       Option for 1-MINUTE average file in $see_data_merged'
  print, '    file=file     Option to specify input data file'
  print, ' '
  print, '    /log          Option to plot Y-axis as log scale'
  print, '    xrange=xrange Option to specify X-range for plot'
  print, '    yrange=yrange Option to specify Y-range for plot'
  print, '    /ignoreflare  Option to set Y-range for daily results instead of flares'
  print, '    /nodayplot    Option to NOT plot the day average (for ORBIT, MINUTE option)'
  print, ' '
  print, '    /f107plot     Option to over plot F10.7 [if you have get_index()]'
  print, ' '
  print, '    ddata=ddata   Optional output for daily average for ORBIT or MINUTE option'
  print, '    /clearcommon  Option to allow switching between ORBIT or DAY data that clears all the common blocks prior to proceeding'
  print, ' '
endif

common xpsL4ts_common, xps, xps_num, ver_num
;   xps is the data structure read from the NetCDF file
;   xps_num is the number of array elements of "xps"
;   ver_num is the version number of the XPS data

common xpsL4extra_common, time, timedays, uniqdays, uniqtime, datatype
;  time is the fractional year time of each data point
;  timedays is not used
;  uniqdays is YYYYDOY dates for each unique date found in the XPS data set
;  uniqtime is "time" associated with the "uniqdays"
;  datatype is 3 for MINUTE and 1 for DAY average

if keyword_set(clearcommon) then begin
    t=temporary(xps)      & t=temporary(xps_num)   & t=temporary(ver_num)
    t=temporary(time)     & t=temporary(timedays)  & t=temporary(uniqdays)
    t=temporary(uniqtime) & t=temporary(datatype)
    t=0
endif

;
;	minuteavg has priority over orbitavg
;	and orbitavg has priority over day average
;	DEFAULT is orbitavg
;
orbitavg = 1
minuteavg = 0
if keyword_set(day) then orbitavg = 0
if keyword_set(minute) then minuteavg = 1


dotimecalc = 0    ; flag to set if read in new file - forces "time" calculation

;
;	Read FILE if keyword option "file" is given
;	ELSE
;	  check to see if need to read a different XPS Level 4 merged file
;		from changing between minute, orbit, or day averaged
;	  if so, then clear "xps"
;
if keyword_set(file) then begin
    ;
    ;	Read FILE if keyword option "file" is given
    ;
    fname = findfile( file, count=count )
    if (count gt 0) then begin
       fname=fname[count-1]
       print, 'Reading XPS L4 merged data from ', fname
       read_netcdf, fname, xps, xpsattr, status
       xps_num = n_elements(xps)
       rdot = strpos( fname, '.ncdf', /reverse_search )
       if (rdot gt 4) then ver_num = long(strmid(fname,rdot-3,3)) $
       else ver_num = -1L
       rbase = strpos( fname, 'xps_L4' )
       datatype = 0
       if (rbase ge 0) then begin
         atype = strupcase(strmid(fname,rbase+6,1))
         if (atype eq '_') then datatype = 1
         if (atype eq 'M') then datatype = 3
         if (atype eq 'A') then datatype = 2
       endif 
    endif else status = -1L
    if (status ne 0) then begin
       print, 'ERROR finding XPS L4 merged data file'
       if keyword_set(debug) then stop, 'Check out error with fname, file, count...'
       return
    endif
    dotimecalc = 1
endif else begin
  ;
  ;	Read Merged File (if necessary) based on the MINUTE, ORBIT, and DAY keywords
  ;
 if (n_elements(xps) ge 2) then begin
  if (not keyword_set(minute)) and (not keyword_set(orbit)) and (not keyword_set(day)) then begin
    orbitavg = 0
    minuteavg = 0
    if (datatype eq 3) then minuteavg = 1
    if (datatype eq 2) then orbitavg = 1
  endif else begin
    if (minuteavg ne 0) and (datatype ne 3) then xps=0
    if (orbitavg ne 0) and (datatype ne 2) then xps=0
    if (orbitavg eq 0) and (minuteavg eq 0) and (datatype ne 1) then xps=0
  endelse
 endif
 ;
 ;	this procedure only knows how to read merged data set (in IDL saveset)
 ;
 if (n_elements(xps) lt 2) then begin
    if ( strupcase(!version.os_family) ne 'UNIX' ) then begin
      data_dir = '' ; for non Sun (X) users
    endif else begin
      data_dir = getenv( 'see_data_merged' )  ; for SUN OS / Unix
      if (strlen(data_dir) gt 0) then data_dir = data_dir + '/' 
    endelse
    print, 'Checking for Merged Save Set...'

    fnamesearch = 'see_xps_L4'
    if (minuteavg ne 0) then begin
      fnamesearch = fnamesearch+'M'
      datatype = 3
    endif else if (orbitavg ne 0) then begin
      fnamesearch = fnamesearch+'A'
      datatype = 2
    endif else begin
      fnamesearch = fnamesearch   ; no letter for daily average  ; +'D'
      datatype = 1
    endelse
    fnamesearch = fnamesearch + '_merged*.ncdf'

    fname = findfile( data_dir+fnamesearch, count=count )
    if (count gt 0) then begin
       fname=fname[count-1]
       print, 'Reading XPS L4 merged data from ', fname
       read_netcdf, fname, xps, xpsattr, status
       xps_num = n_elements(xps)
       rdot = strpos( fname, '.ncdf', /reverse_search )
       if (rdot gt 4) then ver_num = long(strmid(fname,rdot-3,3)) $
       else ver_num = -1L
    endif else status = -1L
    if (status ne 0) then begin
       print, 'ERROR finding XPS L4 merged data file'
       if keyword_set(debug) then stop, 'Check out error with fname, data_dir, fnamesearch, count...'
       return
    endif
    dotimecalc = 1
  endif
endelse

if (dotimecalc ne 0) then begin
   ;
   ;  calculate time as fractional year (and save in common block)
   ;
   timeyd = xps.date + xps.time/3600.D0/24.D0
   yrdays = 365.D0
   year = long(xps.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
   time = year + (long((xps.date-1) mod 1000L) + (xps.time/3600.D0/24.D0)) / (yrdays + yrisleap)
   
   ;  calculate timedays as days since 2002/001
   ; timejd = yd2jd(timeyd)
   ; timedays = timejd  - yd2jd(2002001L)
   
   ;  find uniq days for making daily average later
   uniqindex = uniq(xps.date)
   uniqdays = xps[uniqindex].date
   uniqtime = time[uniqindex]
   ; stop, 'Check out timedays, uniqdays ...'
endif

if (n_elements(xps) lt 2) then stop, 'STOP: not enough XPS Level 4 data...'

if (n_params() lt 1) then begin
	waverange = [0., 7]
	read, 'Enter wavelength range [0, 40] ? ', waverange
endif
if n_elements(waverange) lt 2 then waverange = [ long(waverange-0.5), long(waverange+0.5)]
waverange[0] = long(waverange[0])
if (waverange[0] lt 0) then waverange[0] = 0.0
waverange[1] = long(waverange[1]+0.5)
if (waverange[1] gt 40) then waverange[1] = 40.0

; !fancy=4
; setplot		; font='Helvetica Bold', thick=2, white background
cc = rainbow(7)

orig_xtitle=!xtitle
orig_ytitle=!ytitle
orig_mtitle=!mtitle
!xtitle='Time (year)'
!ytitle='Irradiance (mW/m!U2!N)'
!mtitle='TIMED XPS ' + strtrim(long(waverange[0]),2) + '-' + $
	strtrim(long(waverange[1]),2) + ' nm'

;
;	Make irradiance time series
;
nwv = n_elements(xps[0].modelflux)
if (nwv le 40) then begin
  modelwave = findgen(nwv) + 0.5
  dwave = 1.0
endif else begin
  modelwave = findgen(nwv)/10. + 0.05
  dwave = 0.1
endelse
wwave = where( (modelwave ge waverange[0]) and (modelwave le waverange[1]) )
xdata = fltarr(n_elements(xps))
numday = lonarr(n_elements(xps))
numwave = n_elements(wwave)
for k=0,numwave-1 do begin
  temp = reform(xps.modelflux[wwave[k]])
  wgd = where( temp gt 0, numgd )
  if (numgd gt 0) then begin
  	xdata[wgd] = xdata[wgd] + temp[wgd]
  	numday[wgd] = numday[wgd] + 1
  endif
endfor
;  estimate missing spectral coverage
xdata = xdata * numwave / float(numday > 1)
xdata = xdata * dwave	; convert W/m^2/nm to W/m^2
xdata = xdata * 1000.	; convert W to mW

wgood = where( (xdata gt 0) and (numday gt (numwave/2)) )
ptime = time[wgood]
xdata = xdata[wgood]

;
;  make daily averaged array (only for minute and orbit averaged)
;
if (datatype ne 1) then begin
  ;days = long(timedays)
  ;day1 = min(days)
  ;ndays = max(days) - day1 + 1
  ;adays = indgen(ndays) + day1
  ;afy = yd_to_yfrac(jd2yd(yd2jd(2002001L) + double(adays)))
  
  days = xps[wgood].date
  ndays = n_elements(uniqdays)
  afy = uniqtime
  avgdata = fltarr(ndays)
  for k=0,ndays-1 do begin
    ; wd = where( days eq (day1+k) )
    wd = where( days eq uniqdays[k] )
    if (wd[0] ne -1) then avgdata[k] = total(xdata[wd])/n_elements(wd)
  endfor
  wgdavg = where( avgdata gt 0 )
  extraplots = 1
endif else extraplots = 0

;
;	prepare for plotting
;
if keyword_set(xrange) then xr=xrange else xr=[2002.,long((max(ptime)+0.1)*10.)/10.]
wactive = where((ptime ge xr[0]) and (ptime le xr[1]))

if keyword_set(yrange) then yr=yrange else begin
  yfactor=1.3
  if (keyword_set(ignoreflare)) then ymax=avgdata[0]*yfactor $
  else ymax=max(xdata[wactive])
  yr=[0.7*min(xdata[wactive]),ymax*1.1]
endelse
extralimit = 0.2 ;  fractional year
extraplots = ((xr[1]-xr[0]) gt extralimit)
if keyword_set(nodayplot) then extraplots = 0
if (datatype eq 1) then extraplots = 0

if ((xr[1]-xr[0]) le extralimit) then psym=4 else psym=0

if keyword_set(log) then begin
  plot_io, ptime, xdata, yrange=yr, ys=1, psym=psym, xrange=xr, xs=1
endif else begin
  plot, ptime, xdata, yrange=yr, ys=1, psym=psym, xrange=xr, xs=1
endelse

ans='N'

if (extraplots) then begin
 oplot, afy[wgdavg], avgdata[wgdavg], thick=3, psym=10, color=cc[0]
 y2 = !y.crange[1] - (!y.crange[1] - !y.crange[0])/10.
 y1 = max(avgdata[wgdavg]) * 1.05
 if (y1 gt y2) then y1 = y2
 x1 = xr[0] + 0.05 * (xr[1]-xr[0])
 xyouts, x1, y1, 'Daily Average', color=cc[0]
 
 ans='Y'
endif
; else begin
; oplot, ptime, xdata, color=cc[0], psym=0
; endelse

if keyword_set(f107plot) then ans = 'Y' else ans = 'N'
; read, 'Overplot F10.7 ? (Y/N) ', ans
ans = strupcase(strmid(ans,0,1))
if (ans eq 'Y') then begin
  ii = get_index()
  iitime = ii.tyr
  wlap = where( (iitime ge min(ptime)) and (iitime le max(ptime)) )
  factor = 0.8 * mean(xdata) / mean(ii.ten7[wlap])
  oplot, iitime, ii.ten7 * factor, color=cc[3]
  temp = min( abs(iitime - 46), wmin )
  factor2=factor * 0.40
  x1 = xr[0] + (xr[1]-xr[0])*0.05
  xyouts, x1, ii.ten7[wmin] * factor2, 'F10.7', color=cc[3]
endif

;
;	save the best data in an array
;
if (ans eq 'Y') then ndata=3 else ndata=2
data=dblarr(ndata,n_elements(ptime))
data[0,*] = ptime
data[1,*] = xdata
if (ans eq 'Y') then begin
  wgd7 = where( ii.ten7 gt 0 )
  data[2,*] = interpol(ii.ten7[wgd7], ii.tyr[wgd7], ptime ) > 65.
endif

;
;	save daily values in array
;
if (datatype ne 1) then begin
  if (ans eq 'Y') then ddata = fltarr(3,ndays) else ddata = fltarr(2,ndays)
  ddata[0,*] = afy
  ddata[1,*] = avgdata
  if (ans eq 'Y') then ddata[2,*] = interpol( ii.ten7, iitime, afy )
endif

exit_ts:
; !p.multi=0

if (keyword_set(debug)) then stop, 'STOPPED at end of plot_xps4...'
!xtitle=orig_xtitle
!ytitle=orig_ytitle
!mtitle=orig_mtitle

return
end