;-------------------------------------------------------------
;+
; NAME:
;       eofcalc.pro
; PURPOSE:
;       calculate eofs sequentially
; CALLING SEQUENCE:
;       eofcalc,data,weights,neofs,eofs,pcs,var,/timelast
; INPUTS:
;       data: The time dimension of data must be the first or last dimension.
;       If last, set the /timelast flag.
;       weights: same dimensions as data, except no time dimension
;       neofs: number of eofs/pc to return
; OUTPUTS:
;       eofs, pcs, var
; ALGORITHM
;       iterates between time series and spatial pattern  
;       The data matrix is not pre-multiplied by sqrt(w). Rather,
;       w is carried through the iterative equations. Doing this 
;       avoids dividing by w after the EOFs are calculated. Thus,
;       if some elements of w are zero, this algorithm still works.
; NOTES
;       Please cite: Baldwin, M.P., D.B. Stephenson, I.T. Jolliffe, 
;       Spatial weighting and iterative projection methods for EOFs, 
;       J. Climate, revised June 2008.
;
;       Please send comments and suggestions to help improve this code.
;
;       tolerance  and itmax can be adjusted.
;
;       The EOF time series have unit variance.
;
;       The weighting matrix, w, has the same spatial dimensions as 
;       the data array. It consists of the diagonal elements of w, 
;       as described in the JCLIM paper.
;
;       If you run into memory usage issues, email mark@nwra.com
;       and I will send you a version of this code that is longer,
;       but carefully manages memory (doesn't duplicate arrays).
; MODIFICATION HISTORY:
;       Mark Baldwin 9/1/2003
;       MPB 6/8/2007 keeps w (no pre-multiplication)
;-------------------------------------------------------------
;
pro eofcalc,xx,ww,neofs,eofs,pcs,var,timelast=timelast
s=size(xx,/dimensions) & ndims=n_elements(s)-1 & nn=n_elements(xx)
tolerance=0.0001 & itmax=200
if(keyword_set(timelast)) then X=transpose(reform(xx,nn/s(ndims),s(ndims))) $
  else X=reform(xx,s(0),nn/s(0))
s2=size(X,/dimensions) & n=s2(0) & dim1=s2(1)
eofs=fltarr(dim1,neofs) & pcs=fltarr(n,neofs) & var=fltarr(neofs)
w=reform(ww,dim1) 
v0=0 & for i=0,n-1 do v0=v0+total(w*X(i,*)*X(i,*)) & v00=v0
conv=fltarr(itmax,neofs) ;& sd=stdev(x)
for ii=0,neofs-1 do begin
  y=randomn(seed,n) & y_old=y & e_old=randomn(seed,dim1)
  for jj=0,itmax-1 do begin 
    print, 'eof ',ii,' iteration ',jj
    e=(transpose(X) # y)/total(transpose(y) # y) ; "total" of 1-element vector!
    y=(X # (w*e))/total((transpose(e)#(w*e)))
    conv(jj,ii)=mean(((e-e_old)/stdev(e))^2)
    y_old=y & e_old=e
    if(conv(jj,ii) lt tolerance) then goto, done
  endfor
  print, 'did not converge', jj,' iterations. conv= ',conv(jj-1,ii)
  done: X=x-y # transpose(e) 
  v1=0 & for i=0,n-1 do v1=v1+total(w*X(i,*)*X(i,*))
  var(ii)=(v0-v1)/v00 & v0=v1
  pcs(*,ii)=y/stdev(y) & eofs(*,ii)=e*stdev(y) ; normalized
  print,'eofcalc: eof ',ii,' converged in ',jj,' iterations',var(ii)
endfor
pcs=reform(pcs) & eofs=reform(eofs)
if(keyword_set(timelast)) then eofs=reform(eofs,[s(0:ndims-1),neofs])
if(not keyword_set(timelast)) then eofs=reform(eofs,[s(1:*),neofs])
jump: return
end