eof <- function(lon,lat,array3d, ...) {

# Compute EOFs and PCs of a three-dimensional p x q x n array, with first
# two space dimensions p (longitude) and q (latitude), and third time 
# dimension n. 
#
# Description:
#
#      Returns principal component time series, spatial patterns (loadings),
#      standard deviation of each principal component, percentage of total 
#      variance accounted by each principal component, a vector of longitudes 
#      and a vector of latitudes
#
# Usage:
#
#      eofs(lon,lat,array3d)
#
# Input:
#
#      array3d: a three-demensional array with p longitude points and q latitude
#               points as the first two dimensions and n as the third time 
#               dimension 
#
#      lon: vector of longitudes
#
#      lat: vector of latitudes
#
#      ...: Additional arguments to be passed for prcomp function such as,
#           retx:   logical value indicating whether the rotated variables 
#                   should be returned.
#           center: logical value indicating whether the variables should 
#                   be zero centered. Alternately, a vector of
#                   length equal the number of columns of x can be supplied. 
#                   The value is passed to scale.
#           scale.: logical value indicating whether the variables should 
#                   be scaled to have unit variance before the analysis takes 
#                   place. The default is FALSE for consistency with S, but
#                   in general scaling is advisable. Alternatively, a vector
#                   of length equal the number of columns of x can be supplied. 
#                   The value is passed to scale.
#              tol: value indicating the magnitude below which components 
#                   should be omitted. (Components are omitted if their 
#                   standard deviations are less than or equal to tol times the
#                   standard deviation of the first component). With the default
#                   null setting, no components are omitted. Other settings for 
#                   tol could be tol = 0 or tol = sqrt(.Machine$double.eps), 
#                   which would omit essentially constant components.
# Output:
#
#               $pcs: matrix of principal components (first column contains 
#                     first principal component and so on)
#             $stdev: vector containing the standard deviation of each principal
#                     component, starting from the first.
#               $eof: three-dimensional array containing the loadings.
#      $accountedvar: vector containing the amount of variance accounted for 
#                     each principal component
#               $lon: vector of longitudes
#               $lat: vector of latitudes
#
#
# Authors:
#
#      Dag Johan Steinskog <dag.johan.steinskog@nersc.no> 10 June 2005
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 
#
# Examples:
#
#      lon <- seq(50,70,20)
#      lat <- seq(-10,10,20)
#      t <- 20
#      x <- array(rnorm(length(lon)*length(lat)*t,sd=10),dim=c(2,2,20))
#      eofs(lon,lat,x)

data <- reshapefield(lon,lat,array3d)$out

aux <- prcomp(data, ...)

eofs<-reshapefield(lon,lat,t(aux$r))$out

varsquared <- (aux$sdev)^2
sumvarsquared <- sum(varsquared)

accountedvar<-round((varsquared/sumvarsquared)*100,2)

invisible(list(pcs=aux$x,stdev=aux$sdev,eofs=eofs,accountedvar=accountedvar,lon=lon,lat=lat))
}