;***********************************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;***********************************************************************


;------------------------------------------------------------------------------------------------------
;      Mission:  To reconstruct a data set that from EOFs and  EOF time series, using function eof2data
;------------------------------------------------------------------------------------------------------

begin

test      = cbinread ("/lz11/gxia/MET522/data/RCP4.5/CanESM2/gxia_test.bin",(/64,128,21/), "float")
;printVarSummary(test)


;****************************************************************************
;  print the orignal data, test(0:5,0:5,0) 
;***************************************************************************

print(avg(test(0,0,:)))

a = test(0,0,:)

;******************************************************************************************
;   performing the EOF analysis and then reconstruct the data, following
;   the instruction on https://www.ncl.ucar.edu/Document/Functions/Built-in/eof2data.shtml
;******************************************************************************************

neval = 20                                   ; number of eofs
eof    = eofunc_Wrap(test,neval,False)             ; ev(neval,nlat,nlon)
eof_ts  = eofunc_ts_Wrap(test,eof,False)             ; evts(neval,ntime)

do n=0,neval-1
eof_ts(n,:) = eof_ts(n,:) + eof_ts@ts_mean(n)  ; add time series mean
end do                 


xEof  = eof2data(eof,eof_ts)       ; reconstruct original array


;*********************************************************************************
; print the xEof(0:5,0:5,0), which should be the same as the orignal dataset
;*********************************************************************************

b = xEof(0,0,:)

wks = gsn_open_wks("X11","")
plot = gsn_csm_y(wks, (/a,b/), False)

;**********************************************************************************
; Question: why does the output differ from the original dataset? 
;**********************************************************************************

end