;History
; 2016/05/09  merged UV-, VIS-MFRSR, and AERONET k
; 2015/10/06  update plot for large font size and axes                                             
; 2015/08/15  updated plot for matching the version 3                                              
; 2015/03/10  modify OMI k in periodal average (not matchup with MFRSR due to the lack of OMI data)
; 2015/03/04  add writing SSA & AOD 
; 2015/02/03  change sza_min=20 to sza_min=30 
; 2014/12/15  add mean value in the plot as "STAR"
; 2014/12/08  add OMI k in 354 and 388 nm     
; 2014/11/03  output using the calculated KOC in 305, 311,317, 325, 332, and 368      
; 2014/10/17  modify blue curve to use 368/440 k_calc
; 2014/09/24  modify the plot for pulication       
; 2014/04/22  added output var: cvf, cvc, rvf, rvc 
; 2014/04/21  output median K ( and quadratic fit coefficients) and n values 
;2014/04/18   Logarithmic wavelength axis.  Add AERONET 440 and 670nm
;2014/04/11   whisker plot for spectral dependence of SSA
;2014/01/25   whisker plot
;2014/01/21   bin plot
;2014/01/17: plot SSA_MFRSR - SSA_AERONET vs AOT440
;2013/12/09   add MFRSR Error bars                                                                
;2013/12/13   add filters and AERONET diagnstics: sphericity, water, sky error, sun_error
;2013/12/09   add filters and AERONET diagnstics: sphericity, water, sky error, sun_error
;2013/07/08     color index (FSC_color, cgcolor) are changed to adjust to linux machine
;2013/07/02     To compare the averaged SSA (+/- 30 min) between MFRSR and AERONET

k_scale=1000.0  ; stored k values are k_scaled by 1000.
;kirch_scale=0.12

mode='print'
;mode='screen'
inst = '527'
;inst = '582'
;inst = '582_579'

;filter2
sza_min=30 & sza_max=75 ;We use this range for the paper
aot440_min = 0.20   ; AERONET level 2 criteria
sphericity_min=95     ; minimal accepted sphericity [%]
sky_error_max = 5.0 ; maximal sky error [%]
delta_aot_max=0.01
mw_std_max = 0.02   ; max StDev of MFRSR retrievals of SSA  allowed (not used anymore)
mw_error_max = 0.02  ; max StDev of MFRSR retrievals of SSA  allowed
delta_ssa440_max = 0.01 ; 22 points maximal allowed difference between AERONET SSA440 and MFRSR SSA440

;dir='./' & name='bolivia2' & title=' Santa Cruz, Bolivia  2007' & name_output='Figure2_Bolivia_spec_k_plot_bin_whisker'+'AODge'+string(aot440_min,format='(f4.2)') &  y_min=0./k_scale & y_max=60./k_scale & y_lebel=y_min-2.0/k_scale
dir='./' & name='gsfc' & title='NASA/GSFC' & name_output='fig05a' ;'spec_k_boxplot_'+'AODge'+string(aot440_min,format='(f4.2)') 
y_min=0./k_scale & y_max=35./k_scale & y_lebel=y_min-2.0/k_scale
extension = '.txt' 
;extension = '.it0' ;BC only
 
gsfc_read_matchup_allchan_adderr_v3, extension, inst, dir,name, delta_ssa440_max, sza_min, sza_max, aot440_min, sphericity_min, sky_error_max, delta_aot_max, mw_std_max, mw_error_max, aot_mfrsr_spec, w_mfrsr_spec, k_mfrsr_spec,  aot_aero_spec, w_aero_spec, k_aero_spec, n_aero_spec, cvf, cvc, rvf, rvc,  e_tot_aer,sf,sc,alb440,alb670,alb870,alb1020

;output 
;w_mfrsr_spec 
;305, 311, 325, 332, 340, 380, 415, 440, 441, 500

sample_number=size(w_mfrsr_spec)
print,'Sample=', sample_number[2]

if inst eq '579' then begin
real_chan = [415,440,500,615,673,870,940] 
real_chan2= [440,670,870,1020]
vbin = [415,441,500,615,673,870,940]  ; MFRSR  this just location of the box
vbin2 = [451,670,870,1020];  AERONET
endif

if inst eq '582' then begin
real_chan = [305,311,325,332,340,380,440]	 
real_chan2= [440,670,870,1020]
vbin = [305,311,325,332,340,380,430]  ; MFRSR  this just location of the box
vbin2 = [441,670,870,1020];  AERONET
endif

if inst eq '582_579' then begin
real_chan = [305,311,325,332,340,380,440,415,441,500] 
real_chan2= [440,670,870,1020]
vbin = [305,311,325,332,340,380,430,415,440,500]  ; MFRSR  this just location of the box
vbin2 = [450,670,870,1020];  AERONET
endif

if inst eq '527' then begin
real_chan = [305,311,317,325,332,368,440]	 
real_chan2= [440,670,870,1020]
vbin = [305,311,317,325,332,368,430]  ; MFRSR  this just location of the box
vbin2 = [441,670,870,1020];  AERONET
endif

;     Set plotting mode
if mode eq 'screen' then begin
set_plot,'X' ;&  device,pseudo_color=8
endif else begin
set_plot,'PS' & device,/color,/landscape,filename=name_output+'.ps'
endelse

;k
data = 1./k_scale*k_mfrsr_spec  ; 1000.*k
data2= 1./k_scale*k_aero_spec  ; 1000.*k
data3= 1.      *n_aero_spec  ; Real part of AERONET refractive index

!x.thick=5
!y.thick=5
!P.Thick=5
cgloadct
;y_range = [0.005, 0.055]
y_range = [0.00, 0.030]
plot_oi, [0,1], /nodata, yrange=y_range, xrange=[290, 1050], background=cgcolor('white'), color=cgcolor('black'), $
  charsize=2.7, charthick=5, font=1, xtitle='Wavelength, nm', $
  xtickformat='(A1)', ytitle='Imaginary Refractive Index', YStyle=1, XStyle=1, xticks=1 ;xmargin=[10,10], xtickformat='(A1)' means 'turnnig labels off'
index = vbin & index2 =vbin2 

;UV-MFRSR 305, 311, 325, 332, 340, 380 nm
width = 4 ;nm
cgBoxPlot_nowhisker, data(0:5,*), /overplot, xlocation=index(0:5), width=width, $
  color='light cyan', /fillboxes,  boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfrsr1, thick=5 
med_data1  = stat_mfrsr1.median

;UV-MFRSR 440 nm
width = 8 
cgBoxPlot_nowhisker, data(6,*), /overplot, xlocation=index(6), width=width, $
  color='light cyan', /fillboxes,  boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfrsr2, thick=5 
med_data2  = stat_mfrsr2.median

;AERONET
width = 9 
cgboxplot_nowhisker, data2, /overplot, xlocation=index2, width=width, $
  color=cgcolor('tan6'), boxcolor='tan6', /fillboxes, outlinecolor='red',outliercolor='red', stat=stat_aero, thick=2
med_data4 = stat_aero.median


;#added
;fitting UV regime (BC+BrC) and ViS-NIR regime (BC only)
num_uvchan=6
uv_wave=real_chan[0:num_uvchan-1]
result = poly_fit(uv_wave,stat_mfrsr1.median,1,v_fit,/double)
;w_cont=305. + findgen(127)*0.5
oplot,uv_wave,v_fit,color=cgcolor('blue'),linestyle=2
print, 'UV regime: ', result
print, stat_mfrsr1.median

vis_nir = real_chan2
result2 = poly_fit(vis_nir,stat_aero.median,1,v_fit2,/double)
;w_cont=440. + findgen(1161)*0.5
oplot,vis_nir,v_fit2,color=cgcolor('red'),linestyle=2
print, 'VIS-NIR regime: ', result2
print, stat_aero.median
;	Spectral dependence of OC absorption:	Kirchetter et al., JGR 2004
;                0    1      2    3      4    5     6     7
      lam_kir=[350., 400.,  450.,500. ,550. ,600. ,650. ,700]
      k_oc   =[0.168,0.112,0.063,0.045,0.030,0.017,0.005,0.001]
      k_bc   =[0.770,0.730,0.710,0.720,0.720,0.710,0.720,0.750]
      lam_koc=  [311., 317., 325., 332., 368.]
      koc_fit = [0.620,0.600,0.40, 0.48, 0.1478]

lam_kir2 = findgen(351)+350
k_oc2 = interpol(k_oc, lam_kir, lam_kir2)

; Spectral MG calculated K values
w_mg=[305,311,317,325,332,368,430,670,870,1020]
w_mg2=[305,311,317,325,332,368,440,670,870,1020]

;later update 
aero_k_mg= [1.2592E-02,  1.2592E-02,  1.2592E-02,  1.2592E-02,  1.2592E-02, 1.2592E-02,  1.2592E-02,  1.2592E-02,  1.2592E-02,  1.2592E-02] ;using 670, 870, 1020
;oplot,w_mg2,aero_k_mg,line=0,color=cgcolor('red'), thick=3 

;** outputt k, n  ****
openw, 1,name_output+'.txt'
printf,1,title
printf,1,'K values'
printf,1,'Coefficients of spectral fit polynomial: a0+a1*wavelength'
printf,1,'a0,a1:'
printf,1,result
printf,1,'Polynomial fitted values at MFRSR UV channels (excluding 440nm):'
for i=0, num_uvchan-1 do printf,1,real_chan[i],v_fit[i]

printf,1,'Retrieved k values for each MFRSR channel:'
merge_chan= ['mfrsr'+string(real_chan(*),format='(I3)')+'_nm','aeronet'+string(real_chan2(*),format='(I4)')+'_nm']
merge_med = [med_data1,med_data2,med_data4]

for k=0, n_elements(merge_chan)-1 do printf, 1, merge_chan[k], merge_med[k], format='(A,2x,F0.9)'
printf,1
printf,1,'Real part of AERONET Refractive Index:  n '
printf,1,'AERONET wavelengths:'
printf,1,real_chan2[0],real_chan2[1],real_chan2[2],real_chan2[3]
printf,1,'Median AERONET n for all retrrievals:'
printf,1,median(n_aero_spec[0,*]),median(n_aero_spec[1,*]),median(n_aero_spec[2,*]),median(n_aero_spec[3,*])
printf,1
printf,1, 'Median values of cvf[micron^3/micron^2], cvc[micron^3/micron^2], rvf[micron], rvc[micron]'
printf,1, Median(cvf, /EVEN, /DOUBLE), Median(cvc, /EVEN, /DOUBLE), Median(rvf, /EVEN, /DOUBLE), Median(rvc, /EVEN, /DOUBLE), format='(4(F0.4,2x))'

;  X labels wavelengths
index3 = [305,368,440,670,870]	
;index3 = [305,340,380,440,670,870]	
labels = string(index3,format='(I3.3)')
cgText, index3, y_lebel+0.0003, labels, Alignment=0.5, color=cgcolor('black'), charsize=2.1, charthick=5, font=1

Device, Set_Font='Helvetica Bold', /TT_Font
cgtext,340, 20/k_scale,'UV-MFRSR',color=cgcolor('blue'), charsize=4, charthick=4, font=1
;cgtext,430, 22/k_scale,'VIS-MFRSR',color=cgcolor('Gold'), charsize=3, charthick=3, font=1
cgtext,550, 11/k_scale,'AERONET',color=cgcolor('red'), charsize=4, charthick=4, font=1
print, 'Sample='+string(sample_number[2],format='(i3)')

if mode eq 'print' then device,/close

print, 'imaginary part of the refractive index'
print, med_data1
print, med_data2

;The following contents are the inputs of a 1-D box chemical model simulation
n440 = reform(n_aero_spec [0,*])
print, 'median n(real part of the refractive index):' , median(n440,/even)
print, 'median surface albedo (440 nm):' , median(alb440,/even)
print, 'median sf(sigma from volume median radius for fine mode):' , median(sf,/even) 
print, 'median sc(sigma from volume median radius for coarse mode):' , median(sc,/even)

fine_vc = median(cvf,/even)
coarse_vc = median(cvc,/even)
tot_vc = fine_vc + coarse_vc
print, 'median cvf(volume concentration for fine mode):' , fine_vc
print, 'median cvc(volume concentration for coarse mode):' , coarse_vc  
print, 'total volume concentration (fine+coarse): ', tot_vc 

print, 'median rvf(effective radius for fine mode):' , median(rvf,/even)
print, 'median rvc(effective radius for coarse mode):' , median(rvc,/even)

stop
;for the supplementary table (and chemical model input)
;checking plot
cgBoxPlot, n440, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=real_n, thick=2, ytitle='real part of the refractive index'
med_real_n = real_n.median
print, 'mean   n:' , real_n.mean
print, 'median n:' , real_n.median
print, 'stdev  n:' , real_n.SDEV
print, 'num.p  n:' , real_n.N  ;number of points used to construct the box plot

;sfc_albedo440
cgBoxPlot, alb440, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=stat_alb, thick=2, ytitle='surface albedo (440 nm)'
print, 'median surface albedo (440 nm):' , stat_alb.median

;sf
cgBoxPlot, sf, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=stat_sf, thick=2, ytitle='sigma from volume median radius (fine mode)'
print, 'median sf:' , stat_sf.median
print, 'mean   sf:' , stat_sf.mean

;sc
cgBoxPlot, sc, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=stat_sc, thick=2, ytitle='sigma from volume median radius (coarse mode)'
print, 'median sc:' , stat_sc.median
print, 'mean   sc:' , stat_sc.mean

;cvf
cgBoxPlot, cvf, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=stat_cvf, thick=2, ytitle='volume concentration (fine mode)'
print, 'median cvf:' , stat_cvf.median

;cvc
cgBoxPlot, cvc, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=stat_cvc, thick=2, ytitle='volume concentration (coarse mode)'
print, 'median cvc:' , stat_cvc.median

;rvf
cgBoxPlot, rvf, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=stat_rvf, thick=2, ytitle='effective radius (fine mode)'
print, 'median rvf:' , stat_rvf.median
print, 'mean   rvf:' , stat_rvf.mean
print, 'stdev  rvf:' , stat_rvf.SDEV
print, 'num.p  rvf:' , stat_rvf.N  ;number of points used to construct the box plot

;rvc
cgBoxPlot, rvc, xtickformat='(A1)', $
  color='black', boxcolor='black', outlinecolor='black', outliercolor='black', stat=stat_rvc, thick=2, ytitle='effective radius (coarse mode)'
print, 'median rvc:' , stat_rvc.median
print, 'mean   rvc:' , stat_rvc.mean
print, 'stdev  rvc:' , stat_rvc.SDEV
print, 'num.p  rvc:' , stat_rvc.N  ;number of points used to construct the box plot

stop
;addional plot and calculation of AOD and SSA
;**********calculate median SSA****************
;     Set plotting mode
if mode eq 'screen' then begin
set_plot,'X' ;&  device,pseudo_color=8
endif else begin
set_plot,'PS' & device,/color,/landscape,filename=name_output+'_ssa.ps'
endelse

!x.thick=3
!y.thick=3
cgloadct

;SSA
data_ssa   = w_mfrsr_spec
data_ssa2  = w_aero_spec

y_range = [0.7,1.0]
plot_oi, [0,1], /nodata, yrange=y_range, xrange=[290, 1050], background=cgcolor('white'), color=cgcolor('black'), $
  xmargin=[10,10], $
  xtickformat='(A1)', ytitle='SSA', YStyle=1, XStyle=1, xticks=1, charsize=1.5, charthick=2
width = 4 ;nm MFRSR
cgBoxPlot, data_ssa(0:5,*), /overplot, xlocation=index(0:5), width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfr_ssa1, thick=2
width = 10 ;nm
cgBoxPlot, data_ssa(6,*), /overplot, xlocation=index(6), width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfr_ssa2, thick=2
stat_mfr_ssa = [stat_mfr_ssa1,stat_mfr_ssa2]
med_ssa1 = stat_mfr_ssa.median

width = 10 ;nm  AERONET
cgboxplot, data_ssa2, /overplot, xlocation=index2, width=width, $
  color=cgcolor('tan6'), boxcolor='tan6', /fillboxes, outlinecolor='red',outliercolor='red', stat=stat_aero_ssa, thick=2
med_ssa2 = stat_aero_ssa.median

;  X labels wavelengths
index3 = [305,368,440,670,870]	
labels = string(index3,format='(I3.3)')+'nm'
cgText, index3, y_range[0]-0.01, labels, Alignment=0.5, color=cgcolor('black');, charsize=1.5, charthick=2
;xyouts,0.5, 0.9, title,color=0,charthick=2, charsize=2, /normal, alignment=0.5
print, 'Sample='+string(sample_number[2],format='(i3)')

if mode eq 'print' then device,/close

;stop


;**********calculate median AOD****************
;     Set plotting mode
if mode eq 'screen' then begin
set_plot,'X' ;&  device,pseudo_color=8
endif else begin
set_plot,'PS' & device,/color,/landscape,filename=name_output+'_aod.ps'
endelse

!x.thick=3
!y.thick=3
cgloadct

;AOD
data_aod  = aot_mfrsr_spec
data_aod2 = aot_aero_spec

y_range = [0,3.1]
plot_oi, [0,1], /nodata, yrange=y_range, xrange=[290, 1050], background=cgcolor('white'), color=cgcolor('black'), $
  xmargin=[10,10], charsize=1.5, charthick=2, $
  xtickformat='(A1)', ytitle='AOD', YStyle=1, XStyle=1, xticks=1

width = 4 ;nm MFRSR
cgBoxPlot, data_aod(0:5,*), /overplot, xlocation=index(0:5), width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfr_aod1, thick=2
width = 10 ;nm
cgBoxPlot, data_aod(6,*), /overplot, xlocation=index(6), width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfr_aod2, thick=2
stat_mfr_aod = [stat_mfr_aod1,stat_mfr_aod2]
med_aod1 = stat_mfr_aod.median

width = 10 ;nm  AERONET
cgBoxPlot, data_aod2, /overplot, xlocation=index2, width=width, $
  color=cgcolor('tan6'), boxcolor='tan6', /fillboxes, outlinecolor='red',outliercolor='red', stat=stat_aero_aod, thick=2
med_aod2 = stat_aero_aod.median
print, 'MFRSR AOD: ',   med_aod1
print, 'AERONET AOD: ', med_aod2

;  X labels wavelengths
index3 = [305,380,440,670,870]	
labels = string(index3,format='(I3.3)')+'nm'
cgText, index3, y_range[0]-0.1, labels, Alignment=0.5, color=cgcolor('black'), charsize=1.5, charthick=2
xyouts,0.5, 0.9, title,color=0,charthick=2, charsize=2, /normal, alignment=0.5

if mode eq 'print' then device,/close

printf, 1
printf,1,'Retrieved SSA values for each MFRSR channel:'
merge_chan= ['mfrsr'+string(real_chan(*),format='(I3)')+'_nm','aeronet'+string(real_chan2(*),format='(I4)')+'_nm']
merge_med_ssa = [med_ssa1,med_ssa2]
for k=0, n_elements(merge_chan)-1 do printf, 1, merge_chan[k], merge_med_ssa[k], format='(A,2x,F0.3)'
printf, 1

printf,1,'Retrieved AOD values for each MFRSR channel:'
merge_chan= ['mfrsr'+string(real_chan(*),format='(I3)')+'_nm','aeronet'+string(real_chan2(*),format='(I4)')+'_nm']
merge_med_aod = [med_aod1,med_aod2]
for k=0, n_elements(merge_chan)-1 do printf, 1, merge_chan[k], merge_med_aod[k], format='(A,2x,F0.3)'
close, 1
;cgText,380,y_max-10./k_scale, 'AOD!d440!n>'+string(aot440_min,format='(f4.2)')+'  SZA>'+string(sza_min,format='(i2)')+'!uo!n'+'    Number = '+string(sample_number[2],format='(i3)'),charsize=1.5,color=0

;cgText,380,y_max-15./k_scale, cggreek('Delta',/capital)+'AOD!d440!n <'+string(delta_aot_max,format='(f4.2)')+'    '+cggreek('Delta',/capital)+'(SSA!d440!n) < '+string(delta_ssa440_max,format='(f4.2)'),charsize=1.5,color=0

; print number in each spectral bin
   ;cgText, 6,y_min+0.01, 'Sample='+string(sample_number[2],format='(i3)'),Alignment=0.5;, color=cgcolor('black')
;Kirchstetter plot
;cgAxis, /YAxis, Color='Black', charsize=1.5, charthick=2, YRange=[0,0.45], ystyle=1, ytitle='OC Imaginary Refractive Index', /Save
;cgAxis, /YAxis, Color='Black', charsize=1.5, charthick=2, YRange=[0,0.65], ystyle=1, ytitle='OC Imaginary Refractive Index', /Save
;oplot,lam_kir2,k_oc2,line=0,color=cgcolor('Black'),thick=2
;A = FINDGEN(17) * (!PI*2/16.)
;USERSYM, COS(A), SIN(A), color=cgcolor('Black'), /fill
;oplot,lam_kir,k_oc,psym=8, color=cgcolor('Black'),thick=2

;lam_koc2 = findgen(21)+311.
;koc_fit2 = interpol(koc_fit, lam_koc, lam_koc2)
;oplot,lam_koc2,koc_fit2,line=2, color=cgcolor('Black'),thick=2
;oplot,lam_koc,koc_fit,psym=8, color=cgcolor('Black'),thick=2

stop
end