;----------------------------------------------------------------------------
;+
;  FILE:
;       hypirsno_prime_rac_plot.pro
;
;  PURPOSE:
;	IDL procedure to plot Primary GSICS Re-Analysis Correction (RAC)
;	calculated by hypirsno_prime_rac
;
;  REFERENCE: 
;	None
;
;  COMMON BLOCKS:
;       none
;
;  DEPENDENCIES:
;	msg_rad2tb		;in /tcc1/home/timothyh/idl/
;	msg_tb2rad		;in /tcc1/home/timothyh/idl/
;
;  INPUTS:
;	x1	Structure for coefficients of GSICS Correction with Primary Reference
;	x2	Coefficients of GSICS Correction with Secondary Reference
;	x3	Coefficients of GSICS Correction with Secondary Reference after Delta Correction
;               n.b. x3 is an abbreviation of x{2,1/2} in the ATBD
;	x4	RAC from Primary Ref, merged with Secondary Reference after Delta Correction 
;	x0	Coefficients of Prime GSICS Correction from all References merged cumulatively
;
;  OUTPUTS:
;     	delta 		;NetCDF file of output Prime GSICS Correction
;
;  OPTIONS:
;	Tb=Tb		;scene brightness temperature(s) (default = standard scene)
;	plotvar		;Select what to plot: 
;				dTb=biases, dTba=all biases
;				dd=double differences of Tb biases, 
;				delta=biases from delta correction,
;				ddelta=dd+delta
;				unc=uncertainties,
;				wt=weights, 
;				offset(_se), slope(_se), covar, 
;				ncol
;
;  MODIFICATION HISTORY:
;   Written by:	Tim Hewison 2014-09-03 - Based on ~/satcal/msg-iasi/delta_applied3.pro
;   Modified by	Tim Hewison 2014-11-25 - to include more plot options
;   Modified by	Tim Hewison 2015-07-23 - changed date matching for 'dd' to account for gaps
;-
;----------------------------------------------------------------------------
pro hypirsno_prime_rac_plot, x0, x1, x2, x3, delta, Tb=Tb, crash=crash, plotvar=plotvar

if (n_elements(tb) eq 8) then x1[0].std_scene_tb=tb	;Override standard
if (n_elements(tb) eq 8) then x2[0].std_scene_tb=tb	;Override standard

nch=N_ELEMENTS(x1.offset[*,0])

chname=STRCOMPRESS(STRING(x1.central_wavelength*1e6,FORMAT='("IR",f4.1)'),/REMOVE_ALL)

!p.multi=[0,2,4]
date_label = LABEL_DATE(DATE_FORMAT=['%N','%Y'],/ROUND_UP)
!X.TICKFORMAT = ['LABEL_DATE', 'LABEL_DATE']
!X.TICKUNITS = ['Month','Year'] 
!X.TICKINTERVAL = 1  
!X.MINOR=1;4	;30
!X.TICKLAYOUT = 2  
!X.STYLE=1
!P.CHARSIZE=2;4	;Settings for maximised plotting window
!P.THICK=2
ymax=KEYWORD_SET(tb)?3.0:1.0

FOR i=0,nch-1 DO BEGIN
  Tb0=N_ELEMENTS(tb)?tb:x1[0].std_scene_tb[i]
  Rad=(msg_tb2rad(Tb0,ich=i+1))[0]
  dRad0=TRANSPOSE(x0.offset[i,*]+(x0.slope[i,*]-1)*Rad)
  dRad1=TRANSPOSE(x1.offset[i,*]+(x1.slope[i,*]-1)*Rad)
  dRad2=TRANSPOSE(x2.offset[i,*]+(x2.slope[i,*]-1)*Rad)
  dRad3=TRANSPOSE(x3.offset[i,*]+(x3.slope[i,*]-1)*Rad)
  sdRad0=TRANSPOSE(sqrt(x0.offset_se[i,*]^2+x0.slope_se[i,*]^2*Rad^2+2*x0.covariance[i,*]*Rad))
  sdRad1=TRANSPOSE(sqrt(x1.offset_se[i,*]^2+x1.slope_se[i,*]^2*Rad^2+2*x1.covariance[i,*]*Rad))
  sdRad2=TRANSPOSE(sqrt(x2.offset_se[i,*]^2+x2.slope_se[i,*]^2*Rad^2+2*x2.covariance[i,*]*Rad))
  sdRad3=TRANSPOSE(sqrt(x3.offset_se[i,*]^2+x3.slope_se[i,*]^2*Rad^2+2*x3.covariance[i,*]*Rad))
  dTb0=REFORM(msg_rad2tb(Rad+dRad0,ich=i+1)-Tb0)
  dTb1=REFORM(msg_rad2tb(Rad+dRad1,ich=i+1)-Tb0)
  dTb2=REFORM(msg_rad2tb(Rad+dRad2,ich=i+1)-Tb0)
  dTb3=REFORM(msg_rad2tb(Rad+dRad3,ich=i+1)-Tb0)
  sdTb0=REFORM(msg_rad2tb(Rad+dRad0+sdRad0,ich=i+1)-Tb0-dTb0)
  sdTb1=REFORM(msg_rad2tb(Rad+dRad1+sdRad1,ich=i+1)-Tb0-dTb1)
  sdTb2=REFORM(msg_rad2tb(Rad+dRad2+sdRad2,ich=i+1)-Tb0-dTb2)
  sdTb3=REFORM(msg_rad2tb(Rad+dRad3+sdRad3,ich=i+1)-Tb0-dTb3)
  !P.TITLE=x1.monitored_instrument+' '+chname[i]
  !Y.TITLE=plotvar
  IF plotvar eq 'unc' THEN BEGIN
    PLOT,x1.jd,sdtb1,YTITLE='Uncertainty [K]'
    OPLOT,x2.jd,sdtb2,COLOR=2
    OPLOT,x3.jd,sdtb3,COLOR=3
    OPLOT,x3.jd,sdtb3-sdtb2,COLOR=5
    OPLOT,x0.jd,sdtb0,COLOR=4
  ENDIF
  IF STRMID(plotvar,0,3)  EQ 'dTb' THEN BEGIN
    PLOT,x0.jd,dTb0,YTITLE='Tb Bias [K]',/NODATA,$
      YRANGE=[-1,+1]*ymax*(1+(i EQ 7))-(i EQ 7)
    OPLOT,x1.jd,dTb1,PSYM=-1
    OPLOT,x2.jd,dTb2,COLOR=2,PSYM=-4
    OPLOT,x0.jd,dTb0*0.,LINESTYLE=2 
  END
  IF plotvar EQ 'dTba' THEN BEGIN
    OPLOT,x3.jd,dTb3,COLOR=3;,PSYM=-2
    OPLOT,x0.jd,dTb0,COLOR=4;,PSYM=-2
    OPLOT,x0.jd,dTb0*0.,LINESTYLE=2 
  END
  IF STRMID(plotvar,0,2) EQ 'dd' THEN BEGIN
    ok1=WHERE(x1.jd GE MIN(x2.jd) AND x1.jd LE MAX(x2.jd))
    ok2=WHERE(x2.jd GE MIN(x1.jd) AND x2.jd LE MAX(x1.jd))
;but this is not enough if the time series are discontinuous - so replaced with the following:
    MATCH,x1.jd,x2.jd,ok1,ok2				;http://idlastro.gsfc.nasa.gov/ftp/pro/misc/match.pro
    mdd=dTb1[ok1]-dTb2[ok2]
    sdd=SQRT((sdTb1[ok1]^2-sdTb2[ok2]^2)>0)			;Catch NANs->0
    dmax=max([1.0*abs([mdd+sdd*2,mdd-sdd*2]),0.2*ymax])		;Re-scale if needed
    PLOT,x2.jd[ok2],mdd,YTITLE='Double Difference [K]',YRANGE=[-1,+1]*dmax
    POLYFILL,[x2.jd[ok2],REVERSE(x2.jd[ok2])],[mdd+2*sdd,REVERSE(mdd-2*sdd)],color=15
    POLYFILL,[x2.jd[ok2],REVERSE(x2.jd[ok2])],[mdd+sdd,REVERSE(mdd-sdd)],color=14
    OPLOT,x2.jd[ok2],mdd
    OPLOT,x2.jd[ok2],dTb2[ok2]*0.,LINESTYLE=2 
;    caldat,x2.jd[ok2],mm,dd,yy					;Decorrelate time series 
;    ok2=ok2[where(dd eq 01)]					;by only taking 1st day of each month
;    mdd=mdd[ok2]
    coeff=POLY_FIT(x2.jd[ok2]/365.25,mdd,1,yft,SIGMA=sigma)
    OverSamplingFactor=SQRT((1+A_CORRELATE(mdd,1))/(1-A_CORRELATE(mdd,1)))
    PRINT,'DD Trend=',coeff[1],'+/-',sigma[1]*OverSamplingFactor,' K/yr',FORMAT='(a9,f8.3,a3,f8.3,a5)'
  END
  IF plotvar EQ 'ddelta' THEN BEGIN
    mdd=dTb3-dTb2
    sdd=SQRT((sdTb3^2-sdTb2^2)>0)				;Catch NANs->0
    POLYFILL,[x3.jd,REVERSE(x3.jd)],[mdd+2*sdd,REVERSE(mdd-2*sdd)],color=7
    POLYFILL,[x3.jd,REVERSE(x3.jd)],[mdd+sdd,REVERSE(mdd-sdd)],color=8
    OPLOT,x3.jd,mdd,COLOR=2
;    OPLOT,x3.jd,dTb3*0.,LINESTYLE=2 
    mdd=dTb1[ok1]-dTb2[ok2]
    sdd=SQRT((sdTb1[ok1]^2-sdTb2[ok2]^2)>0)				;Catch NANs->0
    POLYFILL,[x2.jd[ok2],REVERSE(x2.jd[ok2])],[mdd+2*sdd,REVERSE(mdd-2*sdd)],color=15
    POLYFILL,[x2.jd[ok2],REVERSE(x2.jd[ok2])],[mdd+sdd,REVERSE(mdd-sdd)],color=14
    OPLOT,x2.jd[ok2],mdd
    OPLOT,x2.jd[ok2],dTb2[ok2]*0.,LINESTYLE=2 
  END
  IF plotvar EQ 'delta' THEN BEGIN
    mdd=dTb3-dTb2
    sdd=SQRT((sdTb3^2-sdTb2^2)>0)				;Catch NANs->0
    dmax=max([0.5*abs([mdd+sdd*2,mdd-sdd*2]),0.1*ymax])		;Re-scale if needed
    PLOT,x3.jd,mdd,YTITLE='Delta Correction [K]',YRANGE=[-1,+1]*dmax
    POLYFILL,[x3.jd,REVERSE(x3.jd)],[mdd+2*sdd,REVERSE(mdd-2*sdd)],color=15
    POLYFILL,[x3.jd,REVERSE(x3.jd)],[mdd+sdd,REVERSE(mdd-sdd)],color=14
    OPLOT,x3.jd,mdd
    OPLOT,x3.jd,dTb3*0.,LINESTYLE=2 
  END
  IF plotvar eq 'covar' THEN BEGIN
    plot,x1.jd,x1.covariance[i,*]
    oplot,x2.jd,x2.covariance[i,*],color=2
    oplot,x3.jd,x3.covariance[i,*],color=3
    oplot,x0.jd,x0.covariance[i,*],color=4
  ENDIF
  IF plotvar eq 'offset_se' THEN BEGIN
    plot,x1.jd,x1.offset_se[i,*]
    oplot,x2.jd,x2.offset_se[i,*],color=2
    oplot,x3.jd,x3.offset_se[i,*],color=3
    oplot,x0.jd,x0.offset_se[i,*],color=4
  ENDIF
  IF plotvar eq 'slope_se' THEN BEGIN
    plot,x1.jd,x1.slope_se[i,*]
    oplot,x2.jd,x2.slope_se[i,*],color=2
    oplot,x3.jd,x3.slope_se[i,*],color=3
    oplot,x0.jd,x0.slope_se[i,*],color=4
  ENDIF
  IF plotvar eq 'offset' THEN BEGIN
    plot,x1.jd,x1.offset[i,*]
    oplot,x2.jd,x2.offset[i,*],color=2
    oplot,x3.jd,x3.offset[i,*],color=3
    oplot,x0.jd,x0.offset[i,*],color=4
  ENDIF
  IF plotvar eq 'slope' THEN BEGIN
    plot,x1.jd,x1.slope[i,*]
    oplot,x2.jd,x2.slope[i,*],color=2
    oplot,x3.jd,x3.slope[i,*],color=3
    oplot,x0.jd,x0.slope[i,*],color=4
  ENDIF
  IF plotvar eq 'ncol' THEN BEGIN
    plot,x1.jd,x1.number_of_collocations[i,*]
    oplot,x2.jd,x2.number_of_collocations[i,*],color=2
    oplot,x3.jd,x3.number_of_collocations[i,*],color=3
    oplot,x0.jd,x0.number_of_collocations[i,*],color=4
  ENDIF
  IF plotvar EQ 'wt' THEN BEGIN
    PLOT, x0.jd,x0.w0a[i,*],YTITLE='Weight',/NO_DATA;,LINESTYLE=2
;    OPLOT,x0.jd,x0.w3a[i,*],COLOR=3,LINESTYLE=2
;    OPLOT,x0.jd,x0.w0b[i,*],COLOR=5,LINESTYLE=2
;    OPLOT,x0.jd,x0.w3b[i,*],COLOR=6,LINESTYLE=2
    FOR k=0,4 DO $
      OPLOT,x0.jd,x0.ref_wt[i,k,*],COLOR=2+k
  ENDIF
  poly_prop,reform(delta.coeff[i,*]),reform(delta.covar[i,*,*]),rad,rad12,srad12
  dTb12=REFORM(msg_rad2tb(Rad12,ich=i+1)-Tb0)
  sdTb12=REFORM(msg_rad2tb(Rad12+sRad12,ich=i+1)-Tb0-dTb12)
  PRINT,chname[i],dTb12,sdTb12,$
      FORMAT='(a6," Delta(MetopB-A)=",f6.3," +/-",f6.3)'
ENDFOR

IF KEYWORD_SET(crash) THEN STOP

!P.MULTI=0
!P.TITLE=''
!Y.TITLE=''
!X.TICKFORMAT = ''
!X.TICKUNITS = ''
!X.TICKINTERVAL = 0  
!X.MINOR=0	;30
!X.TICKLAYOUT = 0  
!P.CHARSIZE=0
!X.STYLE=0
!P.THICK=0
!P.CHARTHICK=0

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

RETURN
END