function [pcs,eofs,vxpl,lambda,weights]=is_svdcov(lon,lat,data,weight,scaling,trend);
% function [pcs,eofs,vxpl,lambda,weights]=is_svdcov(lon,lat,data,weight,scaling);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
% Computes the EOFs from the covariance matrix by SVD
% Note the difference between the different scalings of the EOFs 
%
% Input: lon      - longitude vector  
%        lat      - latitude vector 
%        data     - data matrix (time,lat,lon)                    
%        weight   - 0: no area weighting (default) 
%                   1: cosine of latitude weighting   
%                   2: sqrt of cosine weighting
%        trend    - 'constant':  detrend using 'constant'
%                   'linear': i  detrend using 'linear'
%
% The scaling uses the 'synthesis' formulation of PCA
% 
%        scaling  - 0: EOFs have unit length
%                      PCs  have unit length     
%                   1: EOFs MULTIPLIED by square root of eigenvalue   
%                      PCs have standard deviation equal to one   
%                      (often used before rotation in sample space) 
%                   2: EOFs DIVIDED by square root of eigenvalue                                             
%                   3: EOFs scaled to original data units    
%                      by regressing the original anomaly data matrix on 
%                      the PCs. PCs have standard deviation equal to one   
%
% Output: pcs  	  - principal components (columns)
%         eofs    - empirical orthogonal functions (columns)
%         lambda  - the eigenvalues of the covariance matrix   
%         vxpl    - variance represented by the EOFs
%         weights - the weight matrix 
%
% Example: [pcs,eofs,v,l,w]=is_svd(lon,lat,time,hgt,1,1);
%
% See also: is_ploteof, is_plotfield, is_rotateeof, is_svdcorr
%
% IS_EOF TOOLBOX
% Ivar Seierstad (C) 2003
% Tested:
% Revised:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Checking input arguments
if nargin<1
help is_svd; return; end

% Preparing data
[d1,d2,d3]=size(data)
data=reshape(data,d1,d2*d3);
matrix=['covariance'];

% Remove mean from time series
%data=detrend(data,'constant');
if nargin<6
data=detrend(data);
else
data=detrend(data,trend);
end;
disp('Data are detrended');
dataold=data;  % keep a copy of the data matrix prior to normalization/weighting

% The input data are weighted by area 
[lon,lat]=meshgrid(lon,lat);     
lon=reshape(lon,1,d2*d3);
lat=reshape(lat,1,d2*d3); 
 if weight==1
     weight=cos(abs(lat)*pi/180); 
     data=data.*meshgrid(weight,1:d1);
     disp('Data are weighted by the cosine of latitude'); 
     weights=weight; % return the weights    
 elseif weight==2   
     %weight=sqrt(cos(abs(lat)*pi/180)); 
     weight=sqrt(cos(abs(double(lat))*pi/180)); 
     data=data.*meshgrid(weight,1:d1);
     disp('Data are weighted by the sqrt of the cosine of latitude');  
     weights=weight; % return the weights
 else 
     weights=NaN; 
     disp('No area weighting');
 end    
% Make mask

mask=data(1,:);
I=find(isnan(mask)==0);

tmp=data;
data=data(:,I);
[d1,d4]=size(data)
     
% EOFs calculated by SVD
%if d1>(d2*d3)
if d1>d4
  [U,L,A]=svd(data,0);
  lambda=L.^2;         % the singular values are the square roots of the eigenvalues
  if scaling~=0
     U=U*L; end       % amplitude information in PCs      
  A=A';                % put the EOFs as columns 
else           
  [A,L,U]=svd(data',0);
  lambda=L.^2;         % the singular values are the square roots of the eigenvalues
  if scaling~=0
     U=L*U; end        % amplitude information in the PCs 
  U=U';                % put the PCs as columns
end;
disp('The eigenvectors are calculated by singular value decomposition')
vxpl=abs(diag(lambda)/sum(diag(lambda))); % Variance represented by the PCs

% Scaling of eigenvectors (EOFs) using the "synthesis" formulation of PCA  
if scaling==0
    disp('EOFs and PCs scaled to unit length'); % default setting                
elseif scaling==1    
    A=A.*sqrt(meshgrid(diag(lambda)/(d1-1),1:size(A,1)));
    U=U./sqrt(meshgrid(diag(lambda)/(d1-1),1:size(U,1)));
    disp('EOFs multiplied by square root of eigenvalue');
    disp('PCs have st.dev 1');
elseif scaling==2
    A=A./sqrt(meshgrid(diag(lambda)/(d1-1),1:size(A,1)));
    U=U.*sqrt(meshgrid(diag(lambda)/(d1-1),1:size(U,1)));
    disp('EOFs divided by square root of eigenvalue');  
end    
% Postprocessing
eofs=A;
pcs=U;
size(A)
size(U)

% Data regressed on principal components, so as to represent
% the amplitude of the data. The EOF is then strictly speaking 
% EOF coefficients

if scaling==3
    disp('EOFs scaled to represent the amplitude of the data');  
    disp('PCs have st.dev 1');
    pcs=pcs./sqrt(meshgrid(diag(lambda)/(d1-1),1:size(pcs,1))); % PCs normalized
    eofs=((pcs*dataold)/(d1-1))'; 
end
      
% Information for plotting. See also is_ploteofs    
mx=max(max(max(eofs(:,1:3))));
mn=min(min(min(eofs(:,1:3))));
disp(['The first 3 EOFs take values in the range: ',num2str(mn),' to ',num2str(mx)]);
