;History
; 2016/11/15  AMT paper for GSFC data
; 2015/10/12  changed Figure caption: Figure 4a -> Figure 4b
; 2015/09/08  Instead of extrapolation, use the direct retrievals from AERONET only inversion
; 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'
inst = '527'

;filter
sza_min=30 & sza_max=75 ;relaxed filter we suggested in this paper
;sza_min=20 & sza_max=75 ;BrC paper
aot440_min = 0.20	; 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
delta_ssa440_max = 0.01	; 22 points maximal allowed difference between AERONET SSA440 and MFRSR SSA440


dir='./' & name='gsfc' & title=' NASA/GSFC' & name_output='fig05c';'spec_aaod_boxplot_'+'AODge'+string(aot440_min,format='(f4.2)') 
y_min=0.001 & y_max=0.1 & y_lebel=y_min-0.0002

extension = '.txt'
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

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 = [441,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=5
!y.thick=5
!P.Thick=5
plot_oo, [0,1], /nodata, yrange=[y_min,y_max], xrange=[290, 1050], background=cgcolor('white'), color=cgcolor('black'), $
  xtickformat='(A1)', ytitle='Aerosol Absorption Optical Depth',  YStyle=1, XStyle=1, xticks=1, $ ;xtickformat='(A1)' means 'turnnig labels off'
  charsize=2.7, charthick=5, xtitle='Wavelength, nm', font=1
index = vbin & index2 =vbin2 

;    MFRSR
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
width = 8 ;nm
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

; 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: ', string(med_data, format='(7(F0.3,2x))')

num_chan=7
;calculation of AAE in the UV wavelength
n_lamda = 127. ;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=305. + findgen(n_lamda)*0.5
xf = alog(w_cont)
yf = exp(result[0]+result[1]*xf)
oplot,w_cont,yf,line=2,color=cgcolor('blue'), thick=5.0
mfr_UV_AAE = string(result[1]*(-1.),format='(f0.1)')
print, 'AAE(305 to 368 nm)='+mfr_UV_AAE

;**********************
;AAE = 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

;AAE=1 (flat)
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 = 9 ;nm  AERONET
cgboxplot_nowhisker, data2, /overplot, xlocation=index2, width=width, $
  color=cgcolor('tan6'), boxcolor='tan6', /fillboxes, outlinecolor='red',outliercolor='red', stat=stat_aero, thick=5
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[1:3],yy2[1:3],1,v_fit2,/double) ;670, 870, and 1020

;calculation of AAE using AERONET wavelength
n_lamda2 = 1161. 
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=5.0 
aero_AAE = string(result2[1]*(-1.),format='(f0.1)')
print, 'AAE(440 to 1020 nm)='+aero_AAE
print, 'AAE(670 to 1020 nm)='+string(result2_limit[1]*(-1.),format='(f0.1)')

;;add later if I receive OMI data fron Hiren
;;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, aaod388, /overplot, xlocation=vbin3, width=width, $
;  color='Pale Green', /fillboxes,  boxcolor='Pale Green', outlinecolor='Olive', outliercolor='Olive', stat=stat_omi1, thick=5
;
;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=5
;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=5,color=cgcolor('Olive')
;plots,[354,354],[omi2_min,omi2_max],thick=5,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=5.0
;omi_AAE = string(result3[1]*(-1.),format='(f0.1)')
;print, 'AAE(354 and 388 nm)='+omi_AAE


;  X labels wavelengths
;index3 = [305,332,368,440,670,870,1020]	
index3 = [305,368,440,670,870]
labels = string(index3,format='(I4)')
cgText, index3, y_lebel, labels, Alignment=0.5, color=cgcolor('black'), charsize=2.1, charthick=5, font=1
print, 'AOD440 > '+string(aot440_min,format='(f4.2)')+'  SZA > '+string(sza_min,format='(i2)')+'    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.87, 0.82, 'c', ALIGNMENT=0.5, /NORMAL, charsize=4, charthick=3, font=1, color=cgcolor('Black')
cgText, 340, 0.05,'AAE='+mfr_UV_AAE,color=cgcolor('blue'),charsize=3.0, charthick=3, font=1
cgText, 500, 0.02,'AAE='+aero_AAE,color=cgcolor('red'),charsize=3.0, charthick=3, font=1
;cgText, 380, 0.38,'OMI AAE='+omi_AAE,color=cgcolor('Olive'),charsize=2.5, charthick=3, font=1

if mode eq 'print' then device,/close

;if you want to calculate the MAE, please comment out 'stop'
stop  ;will be changed later

extension = '.it0' ;Inversion from AERONET only
final_read_matchup_allchan_adderr_v3,extension, 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
data_aero = aot_mfrsr_spec*(1.-w_mfrsr_spec)   ;mfrsr_aaod
data2_aero= aot_aero_spec*(1.-w_aero_spec)    ;aeronet_aaod

plot_oo, [0,1], /nodata, yrange=[y_min,y_max], xrange=[290, 1050], background=cgcolor('white'), color=cgcolor('black'), $
  xtickformat='(A1)', ytitle='AAOD_iter0',  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 inversion (it0)
width = 4 ;nm
cgBoxPlot, data_aero[0:5,*], /overplot, xlocation=index[0:5], width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfrsr1_aero, thick=2
width = 10 ;nm
cgBoxPlot, data_aero[6,*], /overplot, xlocation=index[6], width=width, $
  color='light cyan', /fillboxes, boxcolor='light cyan', outlinecolor='blue', outliercolor='blue', stat=stat_mfrsr2_aero, thick=2
med_data_aero = [stat_mfrsr1_aero.median,stat_mfrsr2_aero.median]

width = 10 ;nm  AERONET
;cgboxplot_noout, data2, /overplot, xlocation=index2, width=width, $
cgboxplot, data2_aero, /overplot, xlocation=index2, width=width, $
  color=cgcolor('tan6'), boxcolor='tan6', /fillboxes, outlinecolor='red',outliercolor='red', stat=stat_aero, thick=2
med_data2_aero = stat_aero.median

;print, med_data - med_data_aero
;print, med_data2- med_data2_aero

bc_aaod = stat_mfrsr1_aero.median ;in UV wavelength (it0)
print, '305, 311, 317, 325, 332, 368 nm'
print, 'MFRSR BC AAOD', bc_aaod

;;extrapolation
;bc_aaod = exp(result2[0]+result2[1]*xx)
;print, '305, 311, 317, 325, 332, 368 nm'
;print, 'MFRSR BC AAOD (extrapolated)', bc_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: '+string(tot_aaod, format='(6(F0.3,2x))')
print, 'BC AAOD:   '+string(bc_aaod , format='(6(F0.3,2x))')
print, 'BrC AAOD:  '+string(brc_aaod, format='(6(F0.3,2x))')

;######
brc_density = 1.2  ;[gm-3] [Turpin and Lim, 2001]
bc_density = 1.8  ;[gm-3] [Bond and Bergstrom, 2006] 

;the following information is exported from Figure 2 program
;fbc      = 0.021
;fbrc_kir = 0.049
;fbrc_sal = 0.273

;updated 04/04/2016
fbc      = 0.019
fbrc_kir = 0.048
fbrc_sal = 0.270

;AERONET
tot_vol_conc_bc = 0.1890   ;[m3/m2]
tot_vol_conc_brc = 0.1890  ;[m3/m2]

;result [g/m^2]
mass_con_brc_kir  = tot_vol_conc_brc*brc_density*fbrc_kir   ;[BrC] (Kirchstetter et al., 2004)
mass_con_brc_sal  = tot_vol_conc_brc*brc_density*fbrc_sal  ;[BrC] (Saleh et al., 2014)
mass_con_bc       = tot_vol_conc_bc *bc_density *fbc

;[m^2/g]
mae_kirch = brc_aaod/mass_con_brc_kir
mae_saleh = brc_aaod/mass_con_brc_sal
print, '********************************************'
print, 'wavelength: 305, 311, 317, 325, 332, 368 nm'
print, 'BrC MAE [Saleh et al., 2014]: '+string(mae_saleh, format='(6(F0.3,2x))')
print, 'BrC MAE [Kirchsetter et al., 2004]: '+string(mae_kirch, format='(6(F0.3,2x))')
print, '********************************************'

;AAE_BrC (AAE for BrC only)
xx3 = alog([real_chan[0:num_chan-2]])
yy3 = alog(brc_AAOD)
result3 = poly_fit(xx3,yy3,1,v_fit,/double)
mfr_UV_AAE_BrC = string(result3[1]*(-1.),format='(f0.1)')
print, 'wavelength: ', [real_chan[0:num_chan-2]]
print, 'BrC AAOD (total AAOD - AAOD (BC only)): ', brc_AAOD
print, 'AAE_BrC(305 to 368 nm)='+mfr_UV_AAE_BrC

;test: total MAE
mae_bcbrc_kirch = tot_aaod/(mass_con_brc_kir+mass_con_bc)
mae_bcbrc_saleh = tot_aaod/(mass_con_brc_sal+mass_con_bc)
print, 'wavelength: 305, 311, 317, 325, 332, 368 nm'
print, 'total MAE [Saleh et al., 2014]: '+string(mae_bcbrc_saleh, format='(6(F0.3,2x))')
print, 'total MAE [Kirchsetter et al., 2004]: '+string(mae_bcbrc_kirch, format='(6(F0.3,2x))')

;Add MAE plot
mae_xran = [300, 380]
;mae_xran = [300, 700]
mae_yran = [0.1, 20]
;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.7, 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)

stop
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