;
; COMP_TIME_PUBLIC.PRO		Matthew DeLand			10/16/07
;
; Plot time series of composite UV irradiance data in any requested wavelength
;   band, using 1 nm averaged data as input
;
;##  05/05/08  Revise to use V3 data
;
PROG_NAME = 'comp_time_public.pro'
;
PRINT,'This program plots a time series of composite UV solar irradiance data'
PRINT,'  averaged over any wavelength region between 120-400 nm, and normalized'
PRINT,'  to the start of the data set.'
;
; Load all DATA;  Irradiances available for 9764 days from 1978/311 - 2005/213
;   ( 7 Nov 1978 - 1 Aug 2005);  Wavelength coverage is 120.5-399.5 nm, binned
;   at 1 nm intervals;  Data source flag values are extracted, but not plotted
;
OPENR, UNIT1, 'comp_dates.unf', /get_lun
OPENR, UNIT2, 'comp_irrad.unf', /get_lun
OPENR, UNIT3, 'comp_flags.unf', /get_lun
INFO1 = fstat(UNIT1) & N_DY = INFO1.size / 4
DNUM_ALL = FLTARR(N_DY) & IRRAD_ALL = FLTARR(N_DY,280) & FLAG_ALL = IRRAD_ALL
READU, UNIT1, DNUM_ALL
READU, UNIT2, IRRAD_ALL
READU, UNIT3, FLAG_ALL
;
; Select wavelength range & extract valid data;  Use loop to allow repeats
;
STR_NEXT = ' '
WHILE (STR_NEXT NE 'N') DO BEGIN
   PRINT, ' '
   READ, '&&&  Enter STARTING, ENDING wavelength to use [e.g. 200,208] ? ', $
         WV_STRT, WV_END
;;   PRINT, ' '
   READ, '%%  Enter number of days for SMOOTHING function '+$
         '(no smoothing = 1) ? ', N_SMTH
   STR_SMTH = STRTRIM(STRING(format='(I2)', fix(N_SMTH)), 2)
;
; Extract appropriate data & create average;  Remember that averaged data are 
;   stored on 0.5 nm centers
;
   NUM_WV = fix(WV_END - WV_STRT)
   SUB_WV1 = fix(WV_STRT-120)
   SUB_WV2 = SUB_WV1 + (NUM_WV-1)
;
   IRRAD_AVG = FLTARR(N_DY)
   FOR I = 0, N_DY-1 DO BEGIN
      DUM_IRRAD = IRRAD_ALL(I, SUB_WV1:SUB_WV2) 
      SUB_DAY = where((DUM_IRRAD GT 0.0), n_day)
      IF (N_DAY GT 0) THEN BEGIN
         IRRAD_AVG(I) = total(DUM_IRRAD(SUB_DAY)) / N_DAY
      ENDIF ELSE BEGIN
         IRRAD_AVG(I) = -999.0
      ENDELSE
   ENDFOR
;
; Screen out dates with no good irradiance values, then apply smoothing if 
;   desired;  Also extract data source flags for first wavelength of band, but 
;   don't plot
;
   SUB_GOOD = where((IRRAD_AVG GT 0.0), n_good)
   IRRAD_GOOD = IRRAD_AVG(SUB_GOOD)
   DNUM_GOOD = DNUM_ALL(SUB_GOOD)
   FLAG_GOOD = FLAG_ALL(SUB_GOOD,SUB_WV1)
;
   IF (N_SMTH GT 1) THEN BEGIN
      IRRAD_SMTH = smooth(IRRAD_GOOD, N_SMTH)
   ENDIF ELSE BEGIN
      IRRAD_SMTH = IRRAD_GOOD
   ENDELSE
;
; Remove any unsmoothed data points, then normalize to first remaining date;  
;   Create approximate fractional time for plotting
;
   IF (N_SMTH GT 1) THEN BEGIN
      S1 = fix(N_SMTH+1)/2 & S2 = N_GOOD-(S1+1)
   ENDIF ELSE BEGIN
      S1 = 0 & S2 = N_GOOD-1
   ENDELSE
;
   IRRAD_SMTH_PLOT = IRRAD_SMTH(S1:S2)
   IRRAD_NORM = IRRAD_SMTH_PLOT/IRRAD_SMTH_PLOT(0)
   PLOT_TIME = 1978.0 + (DNUM_GOOD(S1:S2)/365.25)
;
; PLOTTING SECTION
; ****************
;
   !P.CHARSIZE = 1.2 
   STR_WV1 = STRING(format='(I3)', fix(WV_STRT))
   STR_WV2 = STRING(format='(I3)', fix(WV_END))
   !X.TITLE = '!6DATE'
   !Y.TITLE = '!6NORMALIZED Irradiance [!8relative to November 1978!6]'
   !P.TITLE = '!6Composite UV Irradiance Time Series:  !8'+STR_WV1+'-'+$
              STR_WV2+' nm!6'
;
;;   PRINT, ' '
   PRINT, '%%  Specify NORMALIZED Y-axis minimum, maximum, # of major & '+$
          'minor ticks'
   READ, '%%    [e.g. 0.96,1.12,8,4] ? ', Y_MIN, Y_MAX, Y_MAJOR, Y_MINOR
;
; Create plot, then ask for new plot if desired
;
   PLOT, PLOT_TIME, IRRAD_NORM, /nodata, /xstyle, xrange=[1978,2006], $
         xticks=14, xminor=4, /ystyle, yrange=[Y_MIN,Y_MAX], yticks=Y_MAJOR, $
         yminor=Y_MINOR
   PLOTS, [1978.0,2006.0], [1.0,1.0], lines=2
   OPLOT, PLOT_TIME, IRRAD_NORM
;
   IF (N_SMTH GT 1) THEN BEGIN
      XYOUTS, 0.87, 0.87, /norm, '!5'+STR_SMTH+' day running average!6', $
              chars=1.1, align=1.0
   ENDIF
;
;;   PRINT, ' '
   READ, '@@@  Another plot [Y/N] ? ', STR_NEXT
   STR_NEXT = STRUPCASE(STR_NEXT)
ENDWHILE
;
CLOSE, /all
END