;
; NEW_SPEC_CALC.PRO			Matthew DeLand			02/14/07
;
; Beginning with the Thuillier et al. [2004] reference solar UV irradiance 
;   spectrum, calculate a new spectrum for a selected date using the Mg II index
;   and scaling factors;  Plot the difference for 250-380 nm, and write out a 
;   new absolute spectrum if desired
;
@read_columns
PROG_NAME = 'new_spec_calc.pro'
;
; Get the Thuillier reference spectrum, binned to 1 nm resolution and limited 
;   to 115.5-419.5 nm
;
read_columns, 'thuillier_atlas1_1nm.dat', 2, 305, 0, WV_REF, IRRAD_REF
;
; Get a recent version of the NOAA SEC Mg II index data set, reformatted for 
;   easier electronic use
;
read_columns, 'mgii_sec_y2006d299.dat', 4, 10215, 0, YR_MG, DAY_MG, $
              DNUM_MG, INDEX_MG
;
; Get the scaling factors that relate irradiance variability at a specific 
;   wavelength to Mg II index variability
;
read_columns, 'scalfact_1nm.dat', 3, 231, [1,2], WV_SCAL, SCAL_FACT
;
; Extract only the 250-380 nm portion of the irradiance and scale factors data;
;   Interpolate the scale factors to the irradiance 0.5 nm center grid
;
SUB1 = where(((WV_REF GT 250.0) AND (WV_REF LT 380.0)), n1)
IRRAD_PLOT = IRRAD_REF(SUB1) & WV_PLOT = WV_REF(SUB1)
;
SUB2 = where(((WV_SCAL GE 250.0) AND (WV_SCAL LE 380.0)), n2)
SCAL_PLOT = spline(WV_SCAL(SUB2), SCAL_FACT(SUB2), WV_PLOT)
;
; Get the Mg II index values for the reference spectrum (29 March 1992) and the
;   desired date, then calculate the new irradiance spectrum and the change
;
STR_NEXT = ''
WHILE (STR_NEXT NE 'N') DO BEGIN
   SUB_MG_REF = where(((YR_MG EQ 1992) AND (DAY_MG EQ 89)), n_mg_ref)
   INDEX_MG_REF = INDEX_MG(SUB_MG_REF(0))
;
   PRINT, ' '
   date:  READ, '@@@  Enter YEAR, DAY for new spectrum [e.g. 1997,182] ? ', $
                YR_NEW, DAY_NEW
   STR_YR_NEW = STRING(form='(I4)', fix(YR_NEW))
   STR_DAY_NEW = STRING(form='(I3.3)', fix(DAY_NEW))
   SUB_MG_NEW = where(((YR_MG EQ YR_NEW) AND (DAY_MG EQ DAY_NEW)), n_mg_new)
   INDEX_MG_NEW = INDEX_MG(SUB_MG_NEW(0))
   IF (INDEX_MG_NEW LT 0.0) THEN BEGIN
      PRINT, '!!!!  No valid Mg II data for this date;  Please select again'
      GOTO, date
   ENDIF
;
   INDEX_NORM = INDEX_MG_NEW/INDEX_MG_REF
   PCNT_MG = 100*(INDEX_NORM-1.0)
   STR_PCNT_MG = STRTRIM(STRING(form='(F5.2)', PCNT_MG), 2)
   IF (PCNT_MG GT 0.0) THEN STR_SIGN = '+' ELSE STR_SIGN = ''
   FRAC_CHNG = SCAL_PLOT*(INDEX_NORM-1.0)
   IRRAD_NEW = IRRAD_PLOT*(1.0+FRAC_CHNG)
   PCNT_CHNG = 100*FRAC_CHNG
   PRINT, '^^  Mg II index CHANGE for '+STR_YR_NEW+'/'+STR_DAY_NEW+' = '+$
          STR_SIGN+STR_PCNT_MG+'%  ^^'
;
; Plot the results
;
   !P.CHARSIZE = 1.1
   !X.TITLE = '!6Wavelength [!8nm!6]'
   !Y.TITLE = '!6Irradiance CHANGE [!8percent!6]'
   !P.TITLE = '!17Estimated Solar Irradiance Change from !181992/89!17 to !18'+$
              STR_YR_NEW+'/'+STR_DAY_NEW+'!6'
   PLOT, WV_PLOT, PCNT_CHNG, /nodata, /xstyle, xrange=[244,386], xticks=13, $
         xminor=5, xtickv=((INDGEN(14)*10)+250), /ystyle, yrange=[-6,+4], $
         yticks=10, yminor=5
   PLOTS, [244,386], [0.0,0.0], lines=2, thick=2
   OPLOT, WV_PLOT, PCNT_CHNG, lines=0, thick=3
   XYOUTS, 0.85, 0.85, /norm, '!7D!5Mg II = '+STR_SIGN+STR_PCNT_MG+'%!6', $
           chars=1.2, align=1.0
   STR_TIME = SYSTIME(0)
   STR_TIME_PLOT = STRMID(STR_TIME,11, 9)+STRMID(STR_TIME, 0,11)+$
                   STRMID(STR_TIME,20, 4)
   XYOUTS, 0.99, 0.09, /norm, '!3'+STR_TIME_PLOT+'!6', orient=90, chars=0.8
;
; Ask if a new ASCII file should be written out 
;
   STR_WRITE = ''
   READ, '@@@  Write the new spectrum to a file [Y/N] ? ', STR_WRITE
   STR_WRITE = STRUPCASE(STR_WRITE)
   IF (STR_WRITE EQ 'Y') THEN BEGIN
      OPENW, UNIT1, 'irrad_new_y'+STR_YR_NEW+'d'+STR_DAY_NEW+'.'dat', /get_lun
      FORM = '$(1X, F5.1, 4X, e10.4)'
      FOR I = 0, N1-1 DO PRINTF, UNIT1, FORM, WV_PLOT(I), IRRAD_NEW(I)
      free_lun, UNIT1
   ENDIF
   READ, '<><>  Another plot [Y/N] ? ', STR_NEXT
   STR_NEXT = STRUPCASE(STR_NEXT)
ENDWHILE
;
; All done
;
END