;History
; change to plot StDev as error bars 
;2016/07/22  Tom's comments: So you may want
;			 to show AERONET points in red in Slides 2&3 for AOD > 0.4 only but without
;			 the SZA criteria applied.
;2016/05/03  update program with version 3
;2014/12/04  modification: add aot440_max
;2014/11/19  modification: 1) add 2014 data 2) add 1020 channel
;2013/12/15  reduce time averaging to +/-5min of the AERONET almucantar measurements
;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'

ssa_offset=0.03  ;MFRSR paper
ssa_offset2=0.05

;	relaxed  filtering criteria (GSFC)
sza_min=30 & sza_max=75
;sza_min=38 & sza_max=75 ;minimum SZA = 38.9
aot440_min = 0.20
;aot440_min = 0.15
aot440_idx='n'
;aot440_max = 0.40 & aot440_idx='y'
sphericity_min=95     ; minimal accepted sphericity [%]
;sphericity_min=00     ; minimal accepted sphericity [%] ;need discussion
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.05; max StDev of MFRSR retrievals of SSA  allowed 
mw_error_max = 0.02; max StDev of MFRSR retrievals of SSA  allowed 
;mw_error_max = 0.2	; max StDev of MFRSR retrievals of SSA  allowed 

;please choose which data you plot
site = 'gsfc'
chan='440'
inst = '527' & label = '(a)'
;inst = '582' & label = '(a)'
;inst = '579' & label = '(b)'
if inst eq '527' then mfr  = 'UV'
if inst eq '582' then mfr  = 'UV'
if inst eq '579' then mfr  = 'VIS'

;dir='CHINA2008/'
;file='matchup_v2_out527.dd440_china_20141120.txt_aotge0.15'                   ; China campaign data 
;file='matchup_out527.dd440_china.txt_aotge0.15'                   ; China campaign data 

;dir='GREECE2006/'
;file='matchup_out527.dd440_greece_nb.txt_aotge0.15'                ; Greece campaign data(use previous month nighttime vbias) 
;file='matchup_out527.dd440_greece.txt_aotge0.15'                   ; Greece campaign data 

;dir='BOLIVIA2007/'
;file='matchup_out527.dd440_bolivia2.txt_aotge0.15'                   ; Bolivia campaign data 

if site eq 'gsfc'   then dir='./'        
if site eq 'yonsei' then dir= 'Yonsei/' 

file= 'matchup_out'+inst+'.dd'+chan+'_'+site+'_sfc_albedo.txt_v3_aotge0.15'
ctitle=mfr+'-MFRSR-'+inst+' '+chan+' nm'

;************          read matchup file   **************
ifile=dir+file
get_matchup_no2_v3,ifile,  yy,mm,dd, doy,   UT,   SZA,          mAOT440,mn,         mW440,     mk440,     aAOT440,   aW440,   ak440,   an440,    aAOT670,   aW670,   ak670,   an670,$
 aAOT870,   aW870,   ak870,   an870, aAOT1020, aW1020, ak1020, an1020, $
 Trans,    no2,    water,    sphericity,    sky_error,    sun_error,    mk_error,    mw_error,     mk_std,     mw_std,    cvf,     cvc,     rvf,     rvc,  $
 e_tot_aer, sf, sc, alb440, alb670, alb870, alb1020

;		!!!!!!!!!!   need  AERONET error estimates   !!!!!!!!!!!!!!!!
ak_error = 0.005	; AERONET error estimate in k
aw_error = 0.03		; AERONET error estimate in w
;		!!!!!!!!!!   need  AERONET error estimates   !!!!!!!!!!!!!!!!

;jul_date = julday_original(mm,dd,yy,UT)

;		Filter good retrievals 
idx = where(mW440 ge 0. $
;        and abs(aAOT440 - mAOT440) le delta_aot_max  $
 	and aAOT440 ge aot440_min $
;	and aAOT440 lt aot440_max $
        and SZA le sza_max  			$
        and SZA ge sza_min  			$
        and sphericity ge sphericity_min 	$
        and sky_error le sky_error_max  	$
        and mw_error le mw_error_max		$
;       and mw_std le mw_std_max		$
;       and mW440 le 0.935			$
            , count )
print,'Filters:'
print,'sza_min=',sza_min,'   sza_max=',sza_max
print,'delta_aot_max=',delta_aot_max,'   aot440_min=',aot440_min
print,'sky_error <',sky_error_max,'%'
print,'sphericity > ',sphericity_min,'%'
print,'Number of filtered retrievals=',count
for i=0,count-1 do print,'yy=',yy[idx[i]],'  mm=',mm[idx[i]],'  dd=',dd[idx[i]],'  UT=',UT[idx[i]],'  SZA=',SZA[idx[i]],'  AERONET SSA440=', aW440[idx[i]],'  AOT440=',aAOT440[idx[i]]
year = yy & mon = mm & day = dd

;*******SSA SCATTER PLOT of AERONET(idx) vs. MFRSR527(idx) )*********
print, 'RMS error of SSA'
rms_err = mw_error(idx)
print, 'number of sample: ', n_elements(rms_err)
print, max(rms_err), min(rms_err) 
print, mean(rms_err,/nan), median(rms_err,/even)

xx = aW440(idx)  &  xx_error = aw_error(idx) &  xx_min = xx - xx_error & xx_max = xx + xx_error  < 1.0 ;AERONET
yy = mW440(idx) & yy_error = mw_error(idx) ;& yy_error = mw_std(idx) 
yy_min = yy - yy_error & yy_max = yy + yy_error    < 1.0  ;MFRSR

;;test print
;year = year(idx) & mon = mon(idx) & day = day(idx) & UT = UT(idx)
;sphericity = sphericity(idx)
;for h=0, n_elements(year)-1 do begin
;  sdate = string(year[h], format='(I4)')+string(mon[h],format='(I2.2)')+string(day[h],format='(I2.2)')+' '+string(UT[h],format='(F0.2)')
;  print, sdate, sphericity[h]
;endfor ;h

;xlim=[0.80,1.02] &  ylim=[0.80,1.02] ;GSFC
xlim=[0.75,1.01] &  ylim=[0.75,1.01] ;Yonsei

;xx = aAOT440*(1.-aW440) ;AERONET
;yy = mAOT440*(1.-mW440) ;MFRSR
;xlim=[0.,0.15]  & ylim=[0.0,0.15]

result = POLY_FIT(xx, yy ,1,yfit=y_fit,status=status,CHISQ=chi)
;print, result
r = correlate(xx, yy) &  print, 'Correlation R=',r, '  Slope=',result(1), ' Intercept=',result(0), ' sqrt( CHISQ)=',sqrt(chi)
mom_mW440=moment(mW440(idx)) & mom_aW440=moment(aW440(idx))
print,'Mean (mW440)=',mom_mW440(0), '   StDev(mW440)=',sqrt(mom_mW440(1)) 
print,'Mean (aW440)=',mom_aW440(0), '   StDev(aW440)=',sqrt(mom_aW440(1))

;add RMSD
dif = (yy - xx)^2.
tnum= n_elements(xx)
rmsd= sqrt(total(dif)/tnum)
print, 'num= ', tnum
print, 'RMSD SSA= ', rmsd

xxx = [0,1.1] & yyy=[0,1.1]
;			Set plotting mode
if mode eq 'screen' then begin
;set_plot,'X' & device,pseudo_color=8
device, decomposed=0
window, 1, xs=1000, ys=1000
loadct, 39
endif else begin
set_plot,'PS' & device,/color,/landscape,filename='scatter_plot_'+file+'SZAgt'+strcompress(string(sza_min,format='(i2)'),/remove_all)+'_AOT440gt'+strcompress(string(aot440_min,format='(f3.1)'),/remove_all)+'.ps'
endelse

!x.thick = 5.0
!y.thick = 5.0
!p.thick = 5.0
!p.charthick = 5.0

!p.font=0
DEVICE, SET_FONT='Helvetica Bold', /TT_FONT

A = FINDGEN(17) * (!PI*2/16.)
USERSYM, COS(A), SIN(A), color=cgcolor('blue')  ;, /fill
plot, xx, yy, psym=8, background=255, xmargin=[10,5],ymargin=[7,7], thick=2, charsize=1.7, charthick=2, $
  xtitle='SSA!D'+chan+'!N (AERONET)', color=0, symsize=1.3, $
;  title=ctitle, $
  ytitle='SSA!D'+chan+'!N ('+mfr+'-MFRSR)',xr=xlim,xstyle=1, yr=ylim,ystyle=1

;		1:1 line
oplot, xxx, yyy, psym=0, symsize=2.0, color=0, thick=2.5
oplot, xxx, yyy+ssa_offset, linestyle=2, symsize=2.0, color=0, thick=1.5
oplot, xxx, yyy-ssa_offset, linestyle=2, symsize=2.0, color=0, thick=1.5
oplot, xxx, yyy+ssa_offset2, linestyle=1, symsize=2.0, color=0, thick=1.5
oplot, xxx, yyy-ssa_offset2, linestyle=1, symsize=2.0, color=0, thick=1.5

if (aot440_idx eq 'y') then begin
;XYOUTS, 0.20, 0.73, string(aot440_min,format='(f0.2)')+' <= AOD!d440!n < '+strcompress(string(aot440_max,format='(f0.2)'),/remove_all),color=0,charthick=1.5, charsize=1.3, /normal
cgtext, 0.20, 0.73, string(aot440_min,format='(f0.2)')+' $\leq$ AOD!d440!n < '+strcompress(string(aot440_max,format='(f0.2)'),/remove_all),color=0,charthick=1.5, charsize=1.3, /normal
endif else begin
;XYOUTS, 0.20, 0.73, 'AOD!d440!n > '+strcompress(string(aot440_min,format='(f0.2)'),/remove_all),color=0,charthick=1.5, charsize=1.3, /normal
cgText, 0.20, 0.73, 'AOD!d440!n $\geq$ '+strcompress(string(aot440_min,format='(f0.2)'),/remove_all),color=0,charthick=1.5, charsize=1.3, /normal
endelse

;Tom Eck's comments
;The minimum SZA of the AERONET almucantars
;for the Yonsei Cimel is 40 degrees, and we are likely to include data >40
;degrees SZA as Level 2 in the upcoming Version 3 database.

;XYOUTS, 0.40, 0.73, 'SZA  > '+strcompress(string(sza_min),/remove_all)+'!uo!n',color=0,charthick=1.5, charsize=1.3, /normal
XYOUTS, 0.20, 0.69, 'Number of matchup data: '+strcompress(string(N_elements(xx)),/remove_all),color=0,charthick=1.5, charsize=1.3, /normal
XYOUTS, 0.20, 0.65, 'Correlation coefficient: '+strcompress(string(r,format='(f5.2)'),/remove_all),color=0,charthick=1.5, charsize=1.3, /normal

;	Average SSA
momSSAMFRSR=moment(mW440) 
xyouts,0.20,0.59,'<SSA!dMFRSR!n>= '+string(mom_mW440(0),format='(f7.3)')     ,color=0,charthick=1.5, charsize=1.3, /normal
xyouts,0.20,0.55,'StDev!dMFRSR!n= '+string(sqrt(mom_mW440(1)),format='(f7.3)'),color=0,charthick=1.5, charsize=1.3, /normal

xyouts,0.60,0.34,'<SSA!dAERONET!n>= '+string(mom_aW440(0),format='(f7.3)')      ,color=0,charthick=1.5, charsize=1.3, /normal
xyouts,0.60,0.30,'StDev!dAERONET!n= '+string(sqrt(mom_aW440(1)),format='(f7.3)'),color=0,charthick=1.5, charsize=1.3, /normal

;	PLOT error bars
;for i=0, n_elements(xx)-1,15 do plots,[xx_min[i],xx_max[i]],[yy[i],yy[i]],thick=0.5,color=64 ; AERONET SSA error bars
for i=0, n_elements(yy)-1 do plots,[xx[i],xx[i]],[yy_min[i],yy_max[i]],thick=0.1,color=64 ; MFRSR   SSA error bars

;if you want to add overplot, please comment out next three lines
;close,/all
;if mode eq 'print' then device,/close
;stop

;  Re-plot relaxed criteria
;sza_min=50 & sza_max=75
;sza_min=39.5 & sza_max=75 ;Tom's comments
aot440_min = 0.40
;sphericity_min = 0.
;		Filter good retrievals 
idx = where(mW440 ge 0. $
;        and abs(aAOT440 - mAOT440) le delta_aot_max  $
 	and aAOT440 ge aot440_min $
;	and aAOT440 lt aot440_max $
        and SZA le sza_max  			$
        and SZA ge sza_min  			$
;        and sphericity ge sphericity_min 	$
        and sky_error le sky_error_max  	$
        and mw_error le mw_error_max		$
;       and mw_std le mw_std_max		$
;       and mW440 le 0.935			$
            , count )
xx = aW440(idx) &  xx_error = aw_error(idx) &  xx_min = xx - xx_error & xx_max = xx + xx_error  < 1.0 ;AERONET
yy = mW440(idx) & yy_error = mw_error(idx) ;& yy_error = mw_std(idx) 
yy_min = yy - yy_error & yy_max = yy + yy_error    < 1.0  ;MFRSR

print, inst
print, 're-plot aeronet criteria'
r = correlate(xx, yy) &  print, 'Correlation R=',r, '  Slope=',result(1), ' Intercept=',result(0), ' sqrt( CHISQ)=',sqrt(chi)
mom_mW440=moment(yy) & mom_aW440=moment(xx)
print,'Mean (mW440)=',mom_mW440(0), '   StDev(mW440)=',sqrt(mom_mW440(1))
print,'Mean (aW440)=',mom_aW440(0), '   StDev(aW440)=',sqrt(mom_aW440(1))

;add RMSD
dif = (yy - xx)^2.
tnum= n_elements(xx)
rmsd= sqrt(total(dif)/tnum)
print, 'num= ', tnum
print, 'RMSD SSA= ', rmsd

color=cgcolor('red') 
A = FINDGEN(17) * (!PI*2/16.)
USERSYM, COS(A), SIN(A), color=color, /fill
oplot, xx, yy, psym=8, thick=2, symsize=1.3
;	oPLOT error bars
for i=0, n_elements(yy)-1 do plots,[xx[i],xx[i]],[yy_min[i],yy_max[i]],thick=0.1,color=color ; MFRSR   SSA error bars

close,/all

Device, Set_Font='Helvetica Bold', /TT_Font
cgText, 0.89, 0.25, label, ALIGNMENT=0.5, /NORMAL, charsize=2.5, charthick=3, font=1, color=cgcolor('Black')

if mode eq 'print' then device,/close

;return
print, 'finished'
end