\SET MODE VERIFY
! ef_eof_demo.jnl ( acm  9/2000)

! Description: Demonstration of computing EOFs using the
! External Functions EOF_SPACE, EOF_TFUNC, EOF_STAT

! The functions implements Chelton's '82 method for finding EOFs of gappy
! time series. All the functions perform the same computation but return
! different portions of the results.

!   * EOF_SPACE returns EOF eigenvectors: spatial EOF in x and y with the
!     same units as the incoming data
!   * EOF_TFUNC returns EOF time amplitude functions; dimensionless
!   * EOF_STAT returns statistics on the EOF computation: number of EOFs
!     scaled and returned; Percent variance explained by each EOF; and the
!     eigenvalues.
!
! For all functions the arguments are as follows:
!
!   * 1st argument: Input data field, a function of x, y, and time; may be a
!     function of z
!   * 2nd argument: The minimum percent variance explained by the EOF's that
!     are computed and scaled, e.g. to return EOFs that explain at least 2%
!     variance, use 2.0
!   * 3nd argument: The fraction of each time series that must be present to
!     include it in the calculations, e.g. to use all time series that have
!     at least half the data present, use 0.5
!
!
! EOF functions Example 1: define a function of x,y,time using trig
! functions.  Decompose into spatial and time EOF's and display statistics.
!
PAUSE
LET time = t[t=1-jan-1990:10-jan-1990:24]   ! 24 hour resolution time axis

DEFINE AXIS/x=0:10:0.5 x10
DEFINE AXIS/y=0:10:0.5 y10
DEFINE GRID/x=x10/y=y10/t=time g10x10
SET GRID g10x10

LET fcn1 = 15.* sin(omega1*t)*cos(r)/(r+1)
LET fcn2 = 20.* sin(omega2*t)*(sin(s)-.2*sin(q))/(s+1)

LET r = ((6.*(xpts-x0)^2 + 7.*(ypts-y0)^2)^0.5)
LET s = (((xpts-x1)^2 + 2*(ypts-y1)^2)^0.5)
LET q = ((3*(xpts-x0)^2 + (ypts-y1)^2)^0.5)

LET x0 = 2
LET y0 = 4
LET x1 = 5
LET y1 = 7

LET omega1 = 1/10*2*3.14159
LET omega2 = 2 * omega1

LET sample_function = fcn1 + fcn2
LET xpts = x
LET ypts = y

!  Compute the statistics on the EOF solution.  EOF_STAT returns:
!  for J=1  the number of EOFs returned
!  for J=2  the percent variance explained by each EOF
!  for J=3  the eigenvalue for each EOF

LET estat = eof_stat(sample_function, 0.1)

LIST/I=1/J=1 estat

LIST/I=1:4/J=2 estat

LIST/I=1:4/J=3 estat


!  Plot the original function (averaged over time) and its spatial
!  decomposition by EOF_SPACE

PAUSE

DEFINE VIEW/xlim=0.,.33/ylim=.6,1./text=0.2 vul
DEFINE VIEW/xlim=.33,.66/ylim=.6,1./text=0.2 vuc
DEFINE VIEW/xlim=.66,1./ylim=.6,1./text=0.2 vur

DEFINE VIEW/xlim=0.,.33/ylim=.1,.5/text=0.2 vll
DEFINE VIEW/xlim=.33,.66/ylim=.1,.5/text=0.2 vlc
DEFINE VIEW/xlim=.66,1./ylim=.1,.5/text=0.2 vlr

SET VIEW vul; CONTOUR/TITLE="FCN1" fcn1[l=1:10@ave]
SET VIEW vuc; CONTOUR/TITLE="FCN2" fcn2[l=1:10@ave]
SET VIEW vur; CONTOUR/TITLE="FCN1 + FCN2" sample_function[l=1:10@ave]


LET exy = eof_space(sample_function, 0.1)

SET VIEW vll; CONTOUR/L=1/TITLE="EOF 1" exy
SET VIEW vlc; CONTOUR/L=2/TITLE="EOF 2" exy

PAUSE
CANCEL VIEW

! Now plot the time amplitude functions.

LET etim = eof_tfunc(sample_function, 0.1)

SET VIEW ul
PLOT/I=1/TITLE=taf1 etim

SET VIEW ur
PLOT/I=2/TITLE=taf2 etim

! should be all bad flags...
SET VIEW ll
PLOT/I=3/TITLE=taf3 etim

PAUSE
CANCEL DATA/ALL
CANCEL VARIABLE/ALL
CANCEL REGION
CANCEL VIEW

!  A second example, using the COADS climatology data.
!

USE coads_climatology
SET REGION/X=67w:1w/Y=11S:11N

!  Compute and save the spatial EOF functions.  These have the same units
!  as the data.

LET eof_xyfcn = eof_space(sst, 0.5)
SAVE/CLOBBER/FILE=sst_clim_eof_space.cdf eof_xyfcn


CANCEL DATA/ALL
CANCEL VARIABLE/ALL
CANCEL REGION

USE sst_clim_eof_space.cdf
SET VIEW ul; fill/l=1/TITLE="eof 1" eof_xyfcn; go land
SET VIEW ur; fill/l=2/TITLE="eof 2" eof_xyfcn; go land
SET VIEW ll; fill/l=3/TITLE="eof 3" eof_xyfcn; go land

PAUSE
CANCEL DATA/ALL
CANCEL VARIABLE/ALL
CANCEL REGION

!  Compute the statistics on the EOFs: number of EOFs scaled
!  and returned; Percent variance explained by each EOF; and the
!  eigenvalues.


USE coads_climatology
SET REGION/X=67w:1w/Y=11S:11N
LET eofstat = eof_stat(sst[X=67w:1w,Y=11S:11N], 0.5)

LET nout = eofstat[i=1,j=1]
LET pcts = eofstat[i=1:`nout`,j=2]

LET eigenv = eofstat[i=1:`nout`,j=3]

LIST/I=1:10 nout
LIST/I=1:10 pcts
LIST/I=1:10 eigenv

PAUSE

!  Compute and save time amplitude functions: Note they are dimensionless.

USE coads_climatology
LET eoftime = eof_tfunc(sst[X=67w:1w,Y=11S:11N], 0.5)

SAVE/CLOBBER/FILE=sst_clim_eof_tfunc.cdf eoftime[i=1:`nout`]

CANCEL VIEW
CANCEL VARIABLE/ALL
USE sst_clim_eof_tfunc.cdf

SET VIEW ul
PLOT/I=1/TITLE="time function 1"/VLIMITS=-2:2:0.5 eoftime
SET VIEW ur
PLOT/I=2/TITLE="time function 2"/VLIMITS=-2:2:0.5 eoftime
SET VIEW ll
PLOT/I=3/TITLE="time function 3"/VLIMITS=-2:2:0.5 eoftime

\SET MODE/LAST VERIFY