;History
; 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'
;DEVICE, DECOMPOSED=0

;filter
sza_min=20 & sza_max=75
aot440_min = 0.40	; GSFC
;aot440_min = 0.60	; Bolivia
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
mw_error_max = 0.02  ; max StDev of MFRSR retrievals of SSA  allowed
;delta_ssa440_max = 0.03	; 22 points maximal allowed difference between AERONET SSA440 and MFRSR SSA440
delta_ssa440_max = 0.01	; 22 points maximal allowed difference between AERONET SSA440 and MFRSR SSA440

;dir='GSFC20140418/' & name='20140415' & title='GSFC  2005-2013' & name_output='GSFC_spec_k_plot_bin_whisker'+'AODge'+string(aot440_min,format='(f4.2)') &  y_min=0./k_scale & y_max=35./k_scale & y_lebel=y_min-1/k_scale

dir='./' & name='bolivia2' & title=' Santa Cruz, Bolivia  2007' & name_output='Bolivia_spec_kBrC_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
 
;read_matchup_allchan_adderr_v2,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
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

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

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=3
!y.thick=3
;xran = [290, 1050]
xran = [300, 705]
yran = [0.00, 0.55]
;yran = [0.005, 0.70]
plot_oi, [0,1], /nodata, yrange=yran, xrange=xran, background=cgcolor('white'), color=cgcolor('black'), $
  xmargin=[10,10], charsize=1.5, charthick=2, xtitle='Wavelength, nm', $
  xtickformat='(A1)', ytitle='BrC Imaginary Refractive Index',  YStyle=1, XStyle=1, xticks=1 ;xtickformat='(A1)' means 'turnnig labels off'

;******Added plot********
;fill lowY and highY
;Everything is done. Uncomment and chage this line
x = findgen(64)+305.
w_high = 6.5997153
highY = 0.14201148*(368./x)^w_high
w_low = 6.7891956
lowY = 0.0247268*(368./x)^w_low
low = 305 & high= 368
for i=0, n_elements(x)-1 do begin
low=x[i] & high=low+1
loY=lowY[i] & hiY=highY[i]
;xpoly = [low,  high, high,  low]
;ypoly = [loY,loY,hiY,hiY]
;polyfill, xpoly, ypoly, Color=cgcolor('light cyan')
xpoly = [low,  high, high,  low, low]
ypoly = [loY,loY,hiY,hiY,loY]
cgColorFill, xpoly, ypoly, /data, color='light cyan'
endfor ;i
;******Added plot********

xx = findgen(1020-305+1)+305.

index = vbin & index2 =vbin2 
;	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_oc_max = k_oc*1.1
	k_oc_min = k_oc*0.9

    k_bc   =[0.770,0.730,0.710,0.720,0.720,0.710,0.720,0.750]

    lam_koc_up=  [305., 311., 317., 325., 332., 368., 440.] 
;	koc_fit_up = [0.54, 0.379, 0.41, 0.24, 0.297, 0.14201148, 0.] ;newly calculated (Oct. 27, 2014)
	koc_fit_up = [0.46, 0.31, 0.36, 0.180, 0.255, 0.14201148, 0.] ;newly calculated (Aug. 21, 2015)
	koc_fit_up_max = koc_fit_up*1.2
	koc_fit_up_min = koc_fit_up*0.8

;	 koc_saleh_min_int = [0.089, 0.0608, 0.0665, 0.0384, 0.0474, 0.0223,  0.] 		;foc=2.6800E-01 (Saleh et al., supplementary file: #7 internal mixing in Table 1) 
;	koc_saleh_min_ext = [0.0983, 0.0672, 0.0735, 0.0425, 0.0524, 0.0247268, 0.] 	;foc=2.4300E-01 (Saleh et al., supplementary file: #7 external mixing in Table 1) 
	koc_saleh_min_ext = [0.0820, 0.0547, 0.0645, 0.0313, 0.0450, 0.0247268, 0.] 	;foc=2.4300E-01 (Saleh et al., supplementary file: #7 external mixing in Table 1) 
    koc_saleh_min_ext_max = koc_saleh_min_ext*1.2
    koc_saleh_min_ext_min = koc_saleh_min_ext*0.8

;Kirchstetter plot
test_k = interpol(k_oc,lam_kir,[368,440,670,870,1020])
print, 'linearly interpolated k: ', test_k
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, symsize=1.4

;our kBrC calculation 
lam_koc2 = findgen(136*5)/5+305.
koc_fit2 = interpol(koc_fit_up, lam_koc_up, lam_koc2)
USERSYM, COS(A), SIN(A), color=cgcolor('blue')
oplot,lam_koc_up[0:5],koc_fit_up[0:5],psym=8, linestyle=1, color=cgcolor('blue'),thick=2, symsize=1.4
;oplot,lam_koc_up[5:6],koc_fit_up[5:6],psym=-8, linestyle=1, color=cgcolor('blue'),thick=2, symsize=1.4

;   PLOT error bars
koc_fit_up_max = koc_fit_up_max < yran[1]
for i=0, n_elements(lam_kir)-2 do plots,[lam_kir[i],lam_kir[i]],[k_oc_min[i],k_oc_max[i]],thick=0.1,color=cgcolor('Black')
for i=0, n_elements(lam_koc_up)-2 do plots,[lam_koc_up[i],lam_koc_up[i]],[koc_fit_up_min[i], koc_fit_up_max[i]],thick=0.1,color=cgcolor('blue')

;  X labels wavelengths
;index3 = [310,350,400,440,500,550,650]
index3 = [305,332,368,440,550,650]
labels = string(index3,format='(I3.3)')
cgText, index3, y_lebel-0.015, labels, Alignment=0.5, color=cgcolor('black'), charsize=1.2, charthick=2
;xyouts,0.5, 0.9, title,color=0,charthick=2, charsize=2, /normal, alignment=0.5 

;Overplot the line
;1. fitting (Kirchstetter) 350-700
;305-1020 nm
x = lam_kir
y = k_oc
koc550 = k_oc[4]
expr = '0.03*(550./x)^p[0]' ;Kirchstetter k in 550 nm = 0.03
;expr = 'koc550*(550./x)^p[0]'
start=1
result = MPFITEXPR(expr,x,y,start)
print, 'w_550 Kirchstetter et al. (2004)'
print, 'w in Fitting Kirchstetter Koc: ', result[0]
;xx = findgen(396)+305.
oplot, xx, 0.03*(550./xx)^result[0], color=cgcolor('black'), thick=3, symsize=1.4, linestyle=2
print, 'koc_368: ', 0.03*(550./368.)^result[0]
print, 'koc_440: ', 0.03*(550./440.)^result[0]
print, 'koc_670: ', 0.03*(550./670.)^result[0]
print, 'koc_870: ', 0.03*(550./870.)^result[0]
w_kirch = result[0]

;2. inferred k using k368 from Kirchstetter et al. (2004) 
;305-368 nm
xx = findgen(64)+305.
x = lam_koc_up
y = koc_fit_up
expr = '0.14201148*(368./x)^p[0]' ;interpolated Kirchstetter k in 368 nm = 0.14201148
start=1
result = MPFITEXPR(expr,x,y,start)
;result = MPFITEXPR(expr,x[0:5],y[0:5],start) ;305 to 368 nm: w=6.58
print, 'w_368 Kirchstetter et al. (2004)'
print, 'w in Fitting Koc from Kirchstetter (highY): ', result[0]
highY = 0.14201148*(368./xx)^result[0]
oplot, xx, highY, color=cgcolor('blue'), thick=3, symsize=1.4
w_kirch_mfr = result[0] 

;3. inferred k using k368 from Saleh et al.(2014)
;305-368 nm
xx = findgen(64)+305.
A = FINDGEN(17) * (!PI*2/16.)
USERSYM, COS(A), SIN(A), color=cgcolor('blue')
x = lam_koc_up
y = koc_saleh_min_ext
oplot,x[0:5],y[0:5],psym=8, linestyle=1, color=cgcolor('blue'),thick=2, symsize=1.4
;oplot,x[5:6],y[5:6],psym=-8, linestyle=1, color=cgcolor('blue'),thick=2, symsize=1.4
for i=0, n_elements(lam_koc_up)-2 do plots,[x[i],x[i]],[koc_saleh_min_ext_min[i], koc_saleh_min_ext_max[i]],thick=0.1,color=cgcolor('blue')
expr = '0.0247268*(368./x)^p[0]' ;interpolated Saleh et al.(2014) k in 368 nm
start=1
result = MPFITEXPR(expr,x,y,start)
;result = MPFITEXPR(expr,x[0:5],y[0:5],start) ;305 to 368 nm: w=6.77
print, 'w_368 Saleh et al.(2014)'
print, 'w in Fitting Koc from Saleh et al.(2014) (lowY): ', result[0]
lowY = 0.0247268*(368./xx)^result[0]
oplot, xx, lowY, color=cgcolor('blue'), thick=3, symsize=1.4
w_saleh_mfr = result[0]

;4. the parameterized value from Saleh et al. (2014) 
;305-1020 nm
xx = 305.+findgen(1020-305+1.)
kOA1_lamda=[550]
kOA1 = [0.013] & w1 = [1.6]
print, 'w_550 Saleh et al.(2014): ', w1[0]
yy_saleh_2014 = kOA1[0]*(550./xx)^w1[0]
oplot, xx, yy_saleh_2014, color=cgcolor('Forest Green'), thick=3, linestyle=2, symsize=1.4
w_saleh = w1
cgText, 325, 0.08, 'w='+string(w_saleh_mfr,format='(F0.1)'), color=cgcolor('blue'), charsize=2.5, charthick=3
cgText, 450, 0.1,  'w='+string(w_kirch,format='(F0.1)'), color=cgcolor('black'), charsize=2.5, charthick=3
cgText, 325, 0.4,  'w='+string(w_kirch_mfr,format='(F0.1)'), color=cgcolor('blue'), charsize=2.5, charthick=3

;;background = cgSnapShot()
;;fill lowY and highY
;x = findgen(64)+305.
;low = 305 & high= 368
;for i=0, n_elements(x)-1 do begin
;low=x[i] & high=low+1
;loY=lowY[i] & hiY=highY[i]
;;xpoly = [low,  high, high,  low]
;;ypoly = [loY,loY,hiY,hiY]
;;polyfill, xpoly, ypoly, Color=cgcolor('light cyan')
;xpoly = [low,  high, high,  low, low]
;ypoly = [loY,loY,hiY,hiY,loY]
;cgColorFill, xpoly, ypoly, /data, color='blue'
;endfor ;i
;foreground = cgSnapshot()
;;cgDisplay, Title='Transparent Polygon Window', /FREE
;;cgBlendImage, foreground, background, ALPHA=0.25


;;fitting (combined) 305-700
;x = [lam_koc_up[0:5],lam_kir]
;y = [koc_fit_up[0:5],k_oc]
;print, x, y
;idx = sort(x)
;x = x(idx)
;y = y(idx)
;expr = '0.03*(550./x)^p[0]'
;start=1
;result = MPFITEXPR(expr,x,y,start)
;print, 'w in Fitting the combined Koc: ', result[0]
;;xx = findgen(396)+305.
;;oplot, xx, 0.03*(550./xx)^result[0], color=cgcolor('Olive'), thick=3


if mode eq 'print' then device,/close

stop
end