;History
; 2015/03/10  use OMI periodal avearage instead of matchup (due to the lack of data)
; 2015/02/04  add OMI AAOD w.r.t. wavelength (hard-coded)
; 2014/11/24  add 2014 data and 1020 channel       
; 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 *(-1.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

mode='print'
;mode='screen'

;filter
;sza_min=30 & sza_max=75 ;do not use (we miss one more measurement)
sza_min=20 & sza_max=75 ;BrC paper
aot440_min = 0.40	; GSFC
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 (currently it is not used)
mw_error_max = 0.02  ; max StDev of MFRSR retrievals of SSA  allowed
;mw_error_max = 0.05  ; 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='GSFC20140418/' & name='2005t2014' & title='GSFC  2005-2014' & name_output='GSFC_spec_AAOD_plot_bin_whisker'+'AODge'+string(aot440_min,format='(f4.2)') &  y_min=0.0001 & y_max=0.3 & y_lebel=y_min-(y_min/5.)

dir='./' & name='bolivia2' & title=' Santa Cruz, Bolivia  2007' & name_output='Bolivia_spec_AAOD_whisker'+'AODge'+string(aot440_min,format='(f4.2)') &  y_min=0.01 & y_max=0.6 & y_lebel=y_min-0.0015

read_matchup_allchan_adderr_v3,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
;read_matchup_allchan,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 

data = aot_mfrsr_spec*(1.-w_mfrsr_spec)  ;mfrsr_aaod
data2= aot_aero_spec*(1.-w_aero_spec)    ;aeronet_aaod

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

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 = [440,670,870,1020];  AERONET

;     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

!x.thick=3
!y.thick=3
!P.Thick=3
plot_oo, [0,1], /nodata, yrange=[y_min,y_max], xrange=[290, 1050], background=cgcolor('white'), color=cgcolor('black'), $
  xtickformat='(A1)', ytitle='Absorption Aerosol Optical Depth',  YStyle=1, XStyle=1, xticks=1, $ ;xtickformat='(A1)' means 'turnnig labels off'
  charsize=1.5, charthick=2, xtitle='Wavelength, nm'
;cgPlot, [0,1], /nodata, yrange=[y_min,1], xrange=[300, 700], $
;  xtickformat='(A1)', ytitle='k', title=title, xtitle='wavelength', YStyle=1, XStyle=1, xticks=1 ;xtickformat='(A1)' means 'turnnig labels off'
index = vbin & index2 =vbin2 

;    MFRSR
width = 4 ;nm
;cgBoxPlot_noout, data, /overplot, xlocation=index, width=width, $
cgBoxPlot, 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=2
width = 10 ;nm
cgBoxPlot, data[6,*], /overplot, xlocation=index[6], width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfrsr2, thick=2

; fitting 
med_data = [stat_mfrsr1.median,stat_mfrsr2.median]
tot_aaod = stat_mfrsr1.median ;in UV wavelength
print, '305, 311, 317, 325, 332, 368, 440 nm'
print, 'MFRSR AAOD', med_data, format='(7(F0.3,2x))'

num_chan=7
;calculation of AAE in the UV wavelength
n_lamda = 136. ;number of wavelength
xx = alog([real_chan[0:num_chan-2]]) 
yy = alog([med_data[0:num_chan-2] ]) 
result = poly_fit(xx,yy,1,v_fit,/double)
w_cont=300. + findgen(n_lamda)*0.5
xf = alog(w_cont)
yf = exp(result[0]+result[1]*xf)
oplot,w_cont,yf,line=0,color=cgcolor('blue'), thick=2.0
print, 'AAE(305 to 368 nm)='+string(result[1]*(-1.),format='(f0.1)')
cgText, 317., 0.06,'AAE='+string(result[1]*(-1.),format='(f0.1)'),color=cgcolor('blue'),charsize=2.5, charthick=3

;**********************
;-3
b3 = 3.*xf(n_lamda-1)+ alog(yf(n_lamda-1))
yf3 = exp(b3 + (-3.)*xf)
;oplot,w_cont,yf3,line=4,color=cgcolor('Gold'), thick=2.0

;-1
b1 = 1.*xf(n_lamda-1)+ alog(yf(n_lamda-1))
yf1 = exp(b1 + (-1.)*xf)
;oplot,w_cont,yf1,line=4,color=cgcolor('Gold'), thick=2.0
;**********************

width = 10 ;nm  AERONET
;cgboxplot_noout, data2, /overplot, xlocation=index2, width=width, $
cgboxplot, data2, /overplot, xlocation=index2, width=width, $
  color=cgcolor('tan6'), boxcolor='tan6', /fillboxes, outlinecolor='red',outliercolor='red', stat=stat_aero, thick=2
med_data2 = stat_aero.median
xx2 = alog(real_chan2)
yy2 = alog(med_data2)
result2 = poly_fit(xx2,yy2,1,v_fit2,/double)
result2_limit = poly_fit(xx2[0:2],yy2[0:2],1,v_fit2,/double) ;670, 870, and 1020
;calculation of AAE using AERONET wavelength
n_lamda2 = 1160. 
w_cont=440. + findgen(n_lamda2)*0.5
xf = alog(w_cont)
yf = exp(result2[0]+result2[1]*xf)
oplot,w_cont,yf,line=2,color=cgcolor('red'), thick=2.0 
print, 'AAE(440 to 1020 nm)='+string(result2[1]*(-1.),format='(f0.1)')
print, 'AAE(670 to 1020 nm)='+string(result2_limit[1]*(-1.),format='(f0.1)')
cgText,500, 0.1,'AAE='+string(result2[1]*(-1.),format='(f0.1)'),color=cgcolor('red'),charsize=2.5, charthick=3

;BC AAOD: bc_aaod
bc_aaod = exp(result2[0]+result2[1]*xx)
;total AAOD: tot_aaod
brc_aaod = tot_aaod - bc_aaod
;MAE for BrC [m^2/g]
print, 'wavelength: 305, 311, 317, 325, 332, 368 nm'
print, 'total AAOD', tot_aaod
print, 'BC AAOD', bc_aaod
mae_saleh = brc_aaod/0.05657
mae_kirch = brc_aaod/0.01024
print, 'MAE [Saleh et al., 2014] '+string(mae_saleh, format='(6(F0.3,2x))')
print, 'MAE [Kirchsetter et al., 2004] '+string(mae_kirch, format='(6(F0.3,2x))')

;;how about other wavelength?
;target_wave = 350.
;bc_aaod = exp(result2[0]+result2[1]*target_wave)
;tot_aaod = exp(result[0]+result[1]*target_wave)
;brc_aaod = tot_aaod - bc_aaod
;print, 'MAE for highest fBrC', brc_aaod/0.05657
;print, 'MAE for lowest fBrC', brc_aaod/0.01024

;OMI AAE
;In fact, OMI retrieve only 388 nm
;354 nm is estimated value assuming spectral dependence.
;This is just for visualization. 

;388
vbin3 = 388.
;aaod388 = [0.0645650, 0.129901, 0.200024, 0.177037]
aaod388 = [0.200, 0.139, 0.159, 0.279, 0.177]

;354
vbin4 = 354.
;aaod354 = [0.0710648, 0.143086, 0.252849, 0.227573]
aaod354 = [0.253, 0.179, 0.198, 0.357, 0.228] 

;width = 6 ;nm
width = 10 ;nm
;cgBoxPlot_mod
cgBoxPlot, aaod388, /overplot, xlocation=vbin3, width=width, $
  color='Pale Green', /fillboxes,  boxcolor='Pale Green', outlinecolor='Olive', outliercolor='Olive', stat=stat_omi1, thick=2

width = 8 ;nm ;354 nm is estimated value in OMI
cgBoxPlot, aaod354, /overplot, xlocation=vbin4, width=width, $
  color='Pale Green', /fillboxes,  boxcolor='Pale Green', outlinecolor='Olive', outliercolor='Olive', stat=stat_omi2, thick=2
stat_omi = [stat_omi2.median, stat_omi1.median]

omi1_max = stat_omi1.median*1.3 ;+30%
omi1_min = stat_omi1.median*0.7 ;-30%
omi2_max = stat_omi2.median*1.3 ;+30%
omi2_min = stat_omi2.median*0.7 ;-30%
plots,[388,388],[omi1_min,omi1_max],thick=2,color=cgcolor('Olive')
plots,[354,354],[omi2_min,omi2_max],thick=2,color=cgcolor('Olive')

;calculation of AAE in OMI wavelength
real_chan3 = [354,388]
n_lamda3 = 388-354.+1. ;number of wavelength
w_cont3=354. + findgen(n_lamda3)

xx3 = alog(real_chan3) 
yy3 = alog(stat_omi) 
result3 = poly_fit(xx3,yy3,1,v_fit3,/double)
xf3 = alog(w_cont3)
yf3 = exp(result3[0]+result3[1]*xf3)
oplot,w_cont3,yf3,line=0,color=cgcolor('Olive'), thick=2.0
print, 'AAE(354 and 388 nm)='+string(result3[1]*(-1.),format='(f0.1)')
cgText, 380, 0.38,'AAE='+string(result3[1]*(-1.),format='(f0.1)'),color=cgcolor('Olive'),charsize=2.5, charthick=3

;  X labels wavelengths
index3 = [305,332,368,440,670,870,1020]	
labels = string(index3,format='(I4)')
cgText, index3, y_lebel, labels, Alignment=0.5, color=cgcolor('black'), charsize=1.5, charthick=2
print, 'AOD440'+string(aot440_min,format='(f4.2)')+'  SZA>'+string(sza_min,format='(i2)')+'!uo!n'+'    Number = '+string(sample_number[2],format='(i3)')
print, 'Delta_AOD440 < '+string(delta_aot_max,format='(f4.2)')+'    '+'Delta_SSA440 < '+string(delta_ssa440_max,format='(f4.2)')
;cgText,350,y_max-0.12, '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,350,y_max-0.17, 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')

Device, Set_Font='Helvetica Bold', /TT_Font
cgText, 0.9, 0.85, 'a', ALIGNMENT=0.5, /NORMAL, charsize=4, charthick=3, font=1
if mode eq 'print' then device,/close

stop
;Add MAE plot
mae_xran = [300, 380]
;mae_xran = [300, 700]
mae_yran = [0.1, 16]
;mae_yran = [0, 14]
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+'_MAE.ps'
endelse

!x.thick=3
!y.thick=3
!P.Thick=3
plot_oi, [0,1], /nodata, yrange=mae_yran, xrange=mae_xran, background=cgcolor('white'), color=cgcolor('black'), $
  ytitle='MAE (m!U2!N/g)',  YStyle=1, XStyle=1, $ 
  charsize=1.5, charthick=2, xtitle='Wavelength, nm';, xtickformat='(A1)', xticks=1 ;xtickformat='(A1)' means 'turnnig labels off'

mae_x = [305, 311, 317, 325, 332, 368]
labels = string(mae_x,format='(I4)')
oplot,mae_x,mae_saleh,color=cgcolor('Black'), thick=2.0, psym=cgSymcat(16), symsize=2
oplot,mae_x,mae_kirch,color=cgcolor('Black'), thick=2.0, psym=cgSymcat(9), symsize=2

cgText, mae_x, -0.5, labels, Alignment=0.5, color=cgcolor('black'), charsize=1.5, charthick=2

;Add AOD whisker plot
data = aot_mfrsr_spec
data2= aot_aero_spec
yran = [0,3.2] & y_lebel=yran[0]-0.1
;print, max(data), max(data2)
;print, min(data), min(data2)

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
!P.Thick=3

plot_oi, [0,1], /nodata, yrange=yran, xrange=[290, 1050], background=cgcolor('white'), color=cgcolor('black'), $
  xtickformat='(A1)', ytitle='Aerosol Optical Depth',  YStyle=1, XStyle=1, xticks=1, $ ;xtickformat='(A1)' means 'turnnig labels off'
  charsize=1.5, charthick=2, xtitle='Wavelength, nm'
index = vbin & index2 =vbin2

;    MFRSR
width = 4 ;nm
;cgBoxPlot_noout, data, /overplot, xlocation=index, width=width, $
cgBoxPlot, 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=2
width = 10 ;nm
cgBoxPlot, data[6,*], /overplot, xlocation=index[6], width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfrsr2, thick=2
med_data = [stat_mfrsr1.median,stat_mfrsr2.median]
print, 'AOD: 305, 311, 317, 325, 332, 368, 440 nm'
print, med_data

width = 10 ;nm  AERONET
;cgboxplot_noout, data2, /overplot, xlocation=index2, width=width, $
cgboxplot, data2, /overplot, xlocation=index2, width=width, $
  color=cgcolor('tan6'), boxcolor='tan6', /fillboxes, outlinecolor='red',outliercolor='red', stat=stat_aero, thick=2
med_data2 = stat_aero.median
print, 'AOD: 440, 670, 870, 1020 nm'
print, med_data2

index3 = [305,332,368,440,670,870,1020]
labels = string(index3,format='(I4)')
cgText, index3, y_lebel, labels, Alignment=0.5, color=cgcolor('black'), charsize=1.5, charthick=2


;calculate SSA for BC mixture
SSA_bc = 1. - (bc_aaod/med_data[0:5])
print, 'SSA(BC): 305, 311, 317, 325, 332, 368 nm'
print, SSA_bc
if mode eq 'print' then device,/close
stop
end