; =============================================================
; test rofunc & eof2data using Davis' data
; =============================================================
; John C. Davis
; Statistics and Data Analysis in Geology
; Wiley, 2002, 3rd Edition, p517

; ================================>  ; PARAMETERS
  ntime  = 25                        ; # time steps  (rows)   
  nsta   = 7                         ; # stations    (columns)

  neval  = nsta                      ; # EOFs to calculate (max=nsta)
  npts   = nsta*ntime
; ================================>  ; READ THE ASCII FILE
  dir    = "./"
  fname  = "davis.BoxData.txt"
  data   = asciiread (dir+fname,(/ntime,nsta/), "float")

  data!0 = "time"                    ; name the dimensions
  data!1 = "station"

  opt = True
  opt@tspace = 5 
  opt@title  = "BOX DATA: DAVIS"
  write_matrix (data, (nsta+"f9.3"), opt)
  print("==")

; assorted statistics on data; **not needed** 

  data_avg   = avg(data)           ; overall mean, variance and std
  data_var   = variance(data)
  data_std   = stddev(data)
  print("data_avg="+data_avg+"   data_var="+data_var+"   data_std="+data_std)
  print("==")
 
  data_avg_dims = dim_avg_n(data, 0)
  data_var_dims = dim_variance_n(data, 0)
  data_std_dims = dim_stddev_n(data, 0)
  print("dim_avg_dims: "+data_avg_dims(0)+"  "+data_avg_dims(1)+"  "+data_avg_dims(2)+"  " \
                        +data_avg_dims(3)+"  "+data_avg_dims(4)+"  "+data_avg_dims(5)+"  " \
                        +data_avg_dims(6) )
  print("dim_var_dims: "+data_var_dims(0)+"  "+data_var_dims(1)+"  "+data_var_dims(2)+"  " \
                        +data_var_dims(3)+"  "+data_var_dims(4)+"  "+data_var_dims(5)+"  " \
                        +data_var_dims(6) )
  print("dim_std_dims: "+data_std_dims(0)+"  "+data_std_dims(1)+"  "+data_std_dims(2)+"  " \
                        +data_std_dims(3)+"  "+data_std_dims(4)+"  "+data_std_dims(5)+"  " \
                        +data_std_dims(6) )
   
  print ("==========> EOFUNC <==========")
   
  x        = data(station|:, time|:)     ; 'eofunc' wants time as rightmost dimension

  evecv    = eofunc_Wrap    (x,neval,False)
  evecv_ts = eofunc_ts_Wrap (x,evecv,False)   

  printVarSummary(evecv)
  print ("NCL: evecv@pcvar= " + evecv@pcvar)
  printVarSummary(evecv_ts)

  opt@title  = "Eigenvector components: evecv"
  write_matrix (evecv(station|:,evn|:), (nsta+"f9.3"), opt) ; roorder to match book
  
  opt@title  = "Eigenvector time series: evecv_ts"
  write_matrix (evecv_ts(time|:,evn|:), (nsta+"f9.3"), opt)
  print("==")
  print("evecv_ts@ts_mean="+evecv_ts@ts_mean)
  print("==")
  print("avg(evecv_ts@ts_mean)="+avg(evecv_ts@ts_mean))
  print("==")
   
; ===== Reconstruct data  =========================================

  x0 = eof2data (evecv,evecv_ts)   ; retuens no meta data
  x0!0 = "station"
  x0!1 = "time"
  printVarSummary(x0)
  
  opt@title  = "Reconstructed data via NCL: default mode: constant offset"
  write_matrix (x0(time|:,station|:), (nsta+"f9.3"), opt)

; add the time series component means 

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

  x1      = eof2data (evecv,evecv_ts)
  x1!0 = "station"
  x1!1 = "time"
  opt@title = "Reconstructed data: x1"
  write_matrix (x1(time|:,station|:), (nsta+"f9.3"), opt)