;
;	plot_whi_ref.pro
;
;	Plot WHI reference solar spectral irradiance
;
;	Option to plot ratio instead of spectra
;		Ratio is Mar 29 to Apr 13 if neither ATLAS or VUV2002 option is given
;		or Ratio is to the other reference / model if that option is given
;
;	Options to overplot previous reference spectrum or model
;
;	Options to set Xrange and Yrange on plot
;
;	Options to set Log scale for x or y axis
;
;	Option to not plot March spectrum
;
;	Option to normalize spectra to TSI (removed /normalize option as Ver 2 is normalized to TSI now)
;
;	Option to overplot Planck derivation (temperature change) if doing variability ratio
;
pro plot_whi_ref, atlas=atlas, vuv2002=vuv2002, xrange=xrange, yrange=yrange, ratio=ratio, $
		activeratio=activeratio, xlog=xlog, ylog=ylog, noMarch=noMarch, debug=debug, title=title, $
		tempchange=tempchange 

;  read the WHI reference spectra
ref = read_dat('ref_solar_irradiance_whi-2008_ver2.dat')

; if keyword_set(normalize) then begin
;	ref[1,*] = ref[1,*] / 1.00921
;	ref[2,*] = ref[2,*] / 1.00925
;	ref[3,*] = ref[3,*] / 1.00909
; endif

;  setup default parameters for plotting
if keyword_set(xrange) then xr=xrange else xr=[1,2400]
if keyword_set(yrange) then yr=yrange else yr=[1E-6,1E1]

if (n_elements(xlog) gt 0) then myxlog=xlog else myxlog=0
if (n_elements(ylog) gt 0) then myylog=ylog else myylog=1

if keyword_set(atlas) then begin
 	atlas = read_dat('atlas3_1-nm_sim_res.dat')
endif

if keyword_set(vuv2002) then begin
	vuv2002 = read_dat('vuv2002_whi-2008.dat')
endif

mtitle='LASP SSI Reference for WHI-2008'
if keyword_set(title) then mtitle=title

setplot
cc=rainbow(7)

cs = 1.8		; charsize for xyouts
ct = 1.6		; charthick for xyouts

plo = 0.9	; factor to lower plot min
phi = 1.1	; factor to raise plot max

if keyword_set(ratio) then begin
  ;
  ;	plot ratio - compare to ATLAS or VUV2002 if that option is given
  ;		or else compare March to April
  ;
  if (keyword_set(atlas) or keyword_set(vuv2002)) then begin
    maxbin = long(n_elements(ref[3,*])/10)*10L
    ref1nm = rebin(reform(ref[3,0:maxbin-1]),maxbin/10L)
    wv1nm = rebin(reform(ref[0,0:maxbin-1]),maxbin/10L)
    if keyword_set(vuv2002) then begin
      awave = reform(vuv2002[0,*])
      aratio = interpol(ref1nm,wv1nm,awave) / reform(vuv2002[2,*])
      otherstr = 'VUV2002 Model'
    endif
    if keyword_set(atlas) then begin
      ; ATLAS has higher priority for plotting ratio
      awave = reform(atlas[0,*])
      aratio = interpol(ref1nm,wv1nm,awave) / reform(atlas[1,*])
      otherstr = 'ATLAS-3'
    endif
    nsm = 5  ;  smooth ratio to 5-nm
    aratio = smooth(aratio,nsm)
    if (n_elements(ylog) eq 0) then myylog=0
    if (not keyword_set(yrange)) then begin
      wx = where( (awave ge xr[0]) and (awave le xr[1]), numx )
      if (numx gt 1) then yr = [ min(aratio[wx])*plo, max(aratio[wx])*phi ] else yr=[0,2]
    endif
    plot, awave, aratio, $
		xrange=xr, xstyle=1, yrange=yr, ystyle=1, xlog=myxlog, ylog=myylog, $
		xtitle='Wavelength (nm)', ytitle='Ratio WHI / '+otherstr, $
		title=mtitle
	oplot, xr, [1,1], line=2, color=cc[3]
  endif else begin
    nsm = 11  ;  smooth to 1-nm before doing ratio
    if (xr[0] ge 300) then nsm = 51	; smooth to 5-nm before doing ratio
    ratio1 = smooth(ref[1,*],nsm)/smooth(ref[3,*],nsm) - 1.
    ratio1 = smooth(ratio1,nsm)
    ratio2 = smooth(ref[2,*],nsm)/smooth(ref[3,*],nsm) - 1.
    ratio2 = smooth(ratio2,nsm)
    ratio3 = smooth(ref[1,*],nsm)/smooth(ref[2,*],nsm) - 1.
    ratio3 = smooth(ratio3,nsm)
    aratio = (smooth(ref[1,*],nsm)+smooth(ref[2,*],nsm))/2./smooth(ref[3,*],nsm) - 1.
    aratio = smooth(aratio,nsm)
    if (not keyword_set(yrange)) then begin
      wx = where( (ref[0,*] ge xr[0]) and (ref[0,*] le xr[1]), numx )
      if (numx gt 1) then yr = [ min(abs(ratio1[wx]))*plo, max(abs(ratio1[wx]))*phi ] else yr=[1E-5,1E0]
    endif
    if keyword_set(activeratio) then begin
      ;  if activeratio option given, then plot Active-1 / Active-2 ratio (Sunspot / Plage)
      plot, ref[0,*], ratio3, min_value=yr[0], $
	  	  xrange=xr, xstyle=1, yrange=yr, ystyle=1, xlog=myxlog, ylog=myylog, $
		  xtitle='Wavelength (nm)', ytitle='Variability (Ratio-1)', $
		  title=mtitle
	  if (myylog ne 0) then begin
	    oplot, ref[0,*], -1. * ratio3, $
	  		  min_value=yr[0], color=cc[0]
	  endif
    endif else begin
       ;  plot average (Active-1 and Active-2) ratio to QS
       plot, ref[0,*], aratio, min_value=yr[0], $
	  	  xrange=xr, xstyle=1, yrange=yr, ystyle=1, xlog=myxlog, ylog=myylog, $
		  xtitle='Wavelength (nm)', ytitle='Variability (Ratio-1)', $
		  title=mtitle
	  if (myylog ne 0) then begin
	    oplot, ref[0,*], -1. * aratio, $
	  		  min_value=yr[0], color=cc[0]
	  endif
	endelse
	if keyword_set(tempchange) then begin
		;
		;	overplot Planck derivative with temperature change as example of variability
		;
		wvbb = findgen(3000)+0.5   ; in nm, but use Angstrom in planck.pro calls
		temp = 5700.
		bb = planck(wvbb*10., temp)
		bbdev = planck_dev(wvbb*10., temp)
		bblow = max(bb) * 1D-16
		bbratio1 = (bbdev / (bb > bblow))
		;  keyword tempchange = temperature change for Planck function (e.g. 0.1 K)
		oplot, wvbb, bbratio1*tempchange, color=cc[3]
	endif
 endelse
  
endif else begin
  ;
  ;	plot irradiance spectrum
  ;
  if (not keyword_set(yrange)) then begin
    if (xr[0] lt 1) then pmin = 1 else pmin=xr[0]
    wx = where( (ref[0,*] ge pmin) and (ref[0,*] le xr[1]), numx )
    if (numx gt 1) then yr = [ min(ref[3,wx])*plo, max(ref[3,wx])*phi ] else yr=[1E-5,1E0]
  endif

  plot, ref[0,*], ref[3,*], /nodata, $
	xrange=xr, xstyle=1, yrange=yr, ystyle=1, xlog=myxlog, ylog=myylog, $
	xtitle='Wavelength (nm)', ytitle='Irradiance (W/m!U2!N/nm)', $
	title=mtitle
	
  ; oplot, ref[0,*], smooth(ref[3,*],3), psym=10, color=cc[0]

  dx = (xr[1] - xr[0])/15.
  xx = xr[0] + dx
  if (xr[1] gt 400) then xx = xx + 2*dx
  if (myylog ne 0) then begin
     factor = 10.^((alog10(yr[1])-alog10(yr[0]))/15.)
     yy1 = yr[1] / factor / factor
     yy2 = yy1 / factor
     yy3 = yy2 / factor
     yy4 = yy3 / factor
  endif else begin
     dy = (yr[1]-yr[0])/15.
     yy1 = yr[1] - 2*dy
     yy2 = yy1 - dy
     yy3 = yy2 - dy
     yy4 = yy3 - dy
  endelse
  
  if keyword_set(atlas) then begin
  	oplot, atlas[0,*], atlas[1,*], psym=10, color=cc[5]
  	if keyword_set(noMarch) then ayy = yy2 else ayy = yy3
  	xyouts, xx, ayy, 'ATLAS-3', charsize=cs,charthick=ct, color=cc[5]
  endif

  if keyword_set(vuv2002) then begin
	oplot, vuv2002[0,*], vuv2002[2,*], psym=10, color=cc[6]
  	if keyword_set(noMarch) then ayy = yy3 else ayy = yy4
  	xyouts, xx, ayy, 'VUV2002 Model', charsize=cs,charthick=ct, color=cc[6]
  endif

  ; Option to plot  March 29 - Active Regions
  if (not keyword_set(noMarch)) then begin
    oplot, ref[0,*], ref[1,*], psym=10, color=cc[0]	
    xyouts, xx, yy2, 'WHI Mar 24 - Apr 3', charsize=cs,charthick=ct, color=cc[0]
  endif
  
  if (not keyword_set(noMarch)) or keyword_set(atlas) or keyword_set(vuv2002) then begin
    xyouts, xx, yy1, 'WHI Apr 10-16', charsize=cs,charthick=ct, color=cc[3]
    oplot, ref[0,*], ref[3,*], psym=10, color=cc[3]	; April 13 - Quiet Sun (QS)
  endif else begin
    oplot, ref[0,*], ref[3,*], psym=10		; April 13 - Quiet Sun (QS)    
  endelse
endelse

if keyword_set(debug) then stop, 'STOP: debug data at end of plot_whi_ref.pro...'

return
end