;History
;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'

sza_min=30 & sza_max=65
aot440_min = 0.40
;aot440_min = 0.25
sphericity_min=95     ; minimal accepted sphericity [%]
sky_error_max = 5.0	; maximal sky error [%]
delta_aot_max=0.05


ncolumns=24
chan='440'

;file = 'matchup_out527.t440nm_gsfc20131209_aotge0.25'  ; T: 527 re-process 2007-2013 applying 2007 CUCF calibration for 2007-2009 and 2012 cal to 2010-2012
file = 'matchup_out527.dd440nm_gsfc20131208_aotge0.25'  ; DD:527 re-process 2007-2013 applying 2007 CUCF calibration for 2007-2009 and 2012 cal to 2010-2012
;file = 'matchup_out527.dd440nm_gsfc2005-13.aer025_aotge0.25'        ; 527 re-process 2005-2013 applying old CUCF calibration from 2006/2007
;file = 'matchup_out527.dd440nm_gsfc2005-13.aer025_aotge0.39'        ; 527 re-process 2005-2013 applying old CUCF calibration from 2006/2007
;file = 'matchup_out527.dd440nm_gsfc2010-13.aer025noNO2_aotge0.39'   ; 527 no NO2 absorption applying old CUCF calibration from 2006/2007
;file = 'matchup_out527.dd440nm_gsfc2010-13.aer025.oldcal_aotge0.39'	; 527 applying old CUCF calibration from 2006/2007
;file = 'matchup_out527.dd440nm_gsfc2010-13.aer025.oldcal_aotge0.25'	; 527 applying old CUCF calibration from 2006/2007
;file = 'matchup_out582.dd440nm_gsfc2013.aer025_aotge0.39' ;582
;file = 'matchup_out582.dd440nm_gsfc2013.aer025_aotge0.25' ;582
;file = 'matchup_out579.dd440nm_gsfc2012.aer025_aotge0.25' ;579
;file = './matchup_out579.dd440nm_gsfc2012.mfr025_aotge0.25' ;579

;***********527*****************
hdr=strarr(3)
a=''
close,/all &  openr,1,file &  name_output='Scatter_plot_'+file+'.ps'
readf,1,hdr		;read header
c=0.& data=strarr(ncolumns,5000)
WHILE NOT (EOF(1)) DO BEGIN
    readf,1,a
    x=strsplit (a," ",/extract)
    data(*,c)=x
    c=c+1
ENDWHILE
close, 1
print,'Total number of matchups:',c

data=float(data(*,0:c-1))
yy		=data(0,*) ;year
mm		=data(1,*) ;month
dd		=data(2,*) ;day
doy		=data(3,*) ;day of year
UT		=data(4,*) ;time of day
SZA		=data(5,*) ;solar zenith angle
mAOT440 =data(6,*) ;aot(mfrsr)
mn		=data(7,*) ;number of MFRSR data averaged for one AERONET inversions
mW440	=data(8,*) ;ssa(mfrsr)
mk440	=data(9,*) ;k(mfrsr)
aAOT440	=data(10,*);aot 440nm (AERONET)
aW440	=data(11,*);ssa 440nm (AERONET)
ak440	=data(12,*)
an440	=data(13,*)
aAOT670	=data(14,*)
aW670	=data(15,*)
ak670	=data(16,*)
an670	=data(17,*)
Trans	=data(18,*)
NO2	=data(19,*)
water   =data(20,*)
sphericity	=data(21,*)
sky_error	=data(22,*)
sun_error	=data(23,*)
;mk_error = data(24,*)
mk_error = 0.005	; !!! this needs later estimates
ak_error = 0.005     ;  !!!  need further estimations
;mw_error = data(25,*)
mw_error = 0.02		;  absolute SSA retrieval error: deltaW
aw_error = 0.03       

;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 SZA le sza_max  			$
        and SZA ge sza_min  			$
        and sphericity ge sphericity_min 	$
        and sky_error le sky_error_max, 	$
            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,'%'

;jul_date=jul_date(idx)
mW440 = mW440(idx) & mAOT440 = mAOT440(idx) ;MFRSR
aW440 = aW440(idx) & aAOT440 = aAOT440(idx) ;AERONET

;*******SSA PLOT (AERONET vs. MFRSR 579)*********

xx = aW440 &  xx_error = aw_error &  xx_min = xx - xx_error & xx_max = xx + xx_error  < 1.0 ;AERONET
yy = mW440 & yy_error = mw_error &  yy_min = yy - yy_error & yy_max = yy + yy_error    < 1.0  ;MFRSR
xlim=[0.86,1.02] &  ylim=[0.86,1.02]

;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)
print,'moment(mW440):',moment(mW440) &  print,'moment(aW440):',moment(aW440)

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=900, ys=900
loadct, 39
endif else begin
set_plot,'PS' & device,/color,/landscape,filename='scatter_plot_'+file+'.ps'
endelse

A = FINDGEN(17) * (!PI*2/16.)
USERSYM, COS(A), SIN(A), color=64, /fill
plot, xx, yy, psym=8, background=255, xmargin=[10,5],ymargin=[7,7], thick=1.2, charsize=1.2, charthick=1.5, $
  xtitle='SSA (AERONET)', color=0, symsize=1.5, title=chan+' nm', ytitle='SSA (MFRSR)',xr=xlim,xstyle=1, yr=ylim,ystyle=1
oplot, xxx, yyy, psym=0, symsize=2.0, color=0, thick=1.5
oplot, xx, y_fit, psym=0, symsize=1.2, color=250, thick=1.2
XYOUTS, 0.15, 1.01, 'Total Irradiance Fitting ',color=0
XYOUTS, 0.15, 0.80, 'Number of matchup data: '+strcompress(string(N_elements(xx)),/remove_all),color=0,charthick=1.5, charsize=1.5, /normal
XYOUTS, 0.15, 0.75, 'Correlation coefficient: '+strcompress(string(r),/remove_all),color=0,charthick=1.5, charsize=1.5, /normal
XYOUTS, 0.15, 0.70, 'y=ax+b',color=0,charthick=1.5, charsize=1.5, /normal
XYOUTS, 0.15, 0.65, 'a='+strcompress(string(result(1)),/remove_all)+', b='+strcompress(string(result(0)),/remove_all),color=0,charthick=1.5, charsize=1.5, /normal


;	PLOT error bars
for i=0, n_elements(xx)-1 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


close,/all
if mode eq 'print' then device,/close

;return
end