source: issm/oecreview/Archive/21337-21723/ISSM-21588-21589.diff@ 21726

Index: ../trunk-jpl/externalpackages/pcatool/install.sh
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install.sh    (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install.sh    (revision 21589)
@@ -0,0 +1,12 @@
+#!/bin/bash
+
+#Some cleanup
+rm -rf install pcatool
+mkdir install
+
+#Download from ISSM server
+$ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/pcatool.tar.gz' 'pcatool.tar.gz'
+
+#Untar  into install
+cd install 
+tar -zxvf  ../pcatool.tar.gz

Property changes on: ../trunk-jpl/externalpackages/pcatool/install.sh
___________________________________________________________________
Added: svn:executable
   + *

Index: ../trunk-jpl/externalpackages/pcatool/pcatool.tar.gz
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream

Property changes on: ../trunk-jpl/externalpackages/pcatool/pcatool.tar.gz
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream

Index: ../trunk-jpl/externalpackages/pcatool/install/license.txt
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/license.txt   (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/license.txt   (revision 21589)
@@ -0,0 +1,24 @@
+Copyright (c) 2008, Guillaume MAZE
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in
+      the documentation and/or other materials provided with the distribution
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
Index: ../trunk-jpl/externalpackages/pcatool/install/dispEEOF_v2.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispEEOF_v2.m (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispEEOF_v2.m (revision 21589)
@@ -0,0 +1,76 @@
+% dispEEOF_v2(CHP,PC,EXPVAR,DT,NLAG,MOD) Display few EEOFs with their
+%                                     PC time evolution
+%
+% => DISPLAY FEW EEOFs with their PC time evolution.
+% CHP contains all the EEOFs as EOF*LAG*X*Y.
+% PC contains all the PCs.
+% EXPVAR is a matrix with the explained variance of each
+%  EEOFs in %. This is just for title.
+% DT is time step between maps.
+% NLAG is the number of LAG to display.
+% MOD contains explicitly the EEOFs to display.
+%
+% Rq: it's using the plotm function
+% See also: dispEEOF
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - June 2004
+%  gmaze@univ-brest.fr
+
+function [] = dispEEOF(CHP,PC,EXPVAR,DT,NLAG,MOD)
+
+% Number of EEOF to display     
+NMOD = length(MOD);
+
+% Open figure and first guest
+figur;
+clf;hold on
+typg=8;
+width  = .9/NMOD;
+height = .9/(NLAG+1);
+dleft = (.95-width*NMOD)/4;
+left  = dleft;
+load mapanom
+
+% Let'go :
+
+% We choose an EEOF
+for imod = 1 : NMOD
+    mod = MOD(imod);
+    
+    % and plot maps
+    for iplot=1:NLAG+1
+      bottom = .95-height*iplot;
+      subplot('Position',[left bottom width height]);
+      if iplot==1 % Plot the pc's time serie
+         plot(PC(imod,:));
+         grid on;box on;
+         set(gca,'XAxisLocation','top');axis tight;
+         titre=strcat('EEOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)');
+         titre=strcat(titre,'; DT=',num2str(DT),' ');
+         title(titre);
+         if (imod==1),ax1=get(gca,'YLim');else,set(gca,'YLim',ax1);end
+      else % Plot the lag maps
+         lag=iplot-1;
+         C = squeeze(real(CHP(mod,lag,:,:)));
+%         C = C./xtrm(C); % Eventually normalise field
+         plotm(C,typg);
+         if(lag==1),caxis([-abs(xtrm(C)) abs(xtrm(C))]);cx=caxis;end;
+         caxis(cx);
+         colormap(mapanom);
+         %ylabel(strcat('{\bf',num2str((lag-1)*DT),'}'))
+         m_text(10,-70,strcat('{\bf',num2str((lag-1)*DT),'}'));
+      end %if
+
+    end
+    left = left + width + dleft;
+
+    % Adjust plot width via the plot Aspect ratio
+    ch=get(gcf,'Children');
+    posiPC=get(ch(iplot),'PlotBoxAspectRatio');
+    posiMAP=get(ch(iplot-1),'PlotBoxAspectRatio');
+    set(ch(iplot),'PlotBoxAspectRatio',posiMAP);
+
+end
+
Index: ../trunk-jpl/externalpackages/pcatool/install/confexpvar.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/confexpvar.m  (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/confexpvar.m  (revision 21589)
@@ -0,0 +1,67 @@
+% [CONF,CONFLEVEL] = CONFEXPVAR(A,B,nsamples,N)
+%
+% A(Ntime,MmapA), B(Ntime,MmapB) are time-synchronised 
+%  and time-constant sampled data sets.
+% MmapA and MmapB are not required to be identic.
+% nsamples is what it is (given at least a 100 factor 
+%  will be nice to provide a trustable result).
+% N is the max number of eigen vectors to be computed.
+%
+% CONFEXPVAR Compute the confidence limits for significant 
+% eigenvalues in a SVD decomposition of data sets A and B.
+% This generates nsamples random-in-time data set of B 
+% and finds percent variance of N SVD modes computed 
+% with data set A. 
+% This is exactly the monte-carlo method but based on 
+% a random-in-time only data set. This considerably 
+% decrease the cputime.
+%
+% 2005/11. Guillaume Maze
+% gmaze@univ-brest.fr
+
+function [conf,lev,dsumCF] = confexpvar(A,B,nsamples,N)
+  
+% We try to find the p confidence levels of expvar.
+% Keeping as it is the first field, we determine SVDs with the 
+% second field arbitrary sorted in time.
+% 1st field : A
+% 2nd field : B
+
+n = size(A,1);
+if size(B,1)~=n
+  disp('Error, A and B required to be of the same time length !');
+  disp('See help confexpvar');
+  stop
+end
+  
+timer = cputime;
+disp('Compute confidence levels of explained variance...');
+for iR=1:nsamples
+    if iR==fix(nsamples/2)
+       disp(strcat(num2str(iR),'/',num2str(nsamples)));
+       disp(['Please, still wait a little time ... around '...
+             num2str(cputime-timer) ' cpu seconds']);
+    end
+    
+    % Generate a random new time index:
+    t=randperm(n);
+    % and sort the 2nd field with it:
+    Br=B(t,:);
+
+    % Now we compute SVD and expvar:
+    C=A'*Br;
+    [U,Lambda,V] = svds(C,N);
+    L2 = Lambda.^2;   dsum = diag(L2)/trace(L2);
+    L  = Lambda;      dsum =  diag(L)/trace(L);
+    dsumCF(iR,:) = dsum';
+
+end %for iR
+disp(['OK, this took ' num2str(cputime-timer) ' cpu seconds'])
+
+% CF(nsamples,N) contains expvar of the nsamples realisations:
+CF=sort(dsumCF,1);
+
+% Now we just take the decadal significant levels:
+lev=fix([nsamples/10:nsamples/10:nsamples]);
+conf = CF(lev,:) ; 
+lev=(lev./max(lev))*100;
Index: ../trunk-jpl/externalpackages/pcatool/install/calsvd.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/calsvd.m      (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/calsvd.m      (revision 21589)
@@ -0,0 +1,54 @@
+% [SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = CALSVD(A,B,N) Compute SVDs
+%
+% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - 
+%      Analyses of climatic Data" 1997
+%================================================================
+%
+%  Guillaume MAZE - LPO/LMD - March 2004
+%  gmaze@univ-brest.fr
+
+
+   function [e1,e2,pc1,pc2,expvar,Lambda] = calsvd(A,B,N);
+
+
+disp('====> Let''go for SVDs and pc computing...')
+
+% Assume that A is (time*map) matrix
+[n p]=size(A);
+
+% Remove the mean of each column (ie the time mean in each station records)
+S=detrend(A,'constant');
+P=detrend(B,'constant');
+
+% Form the covariance matrix:
+C = S'*P;
+
+% Find eigenvectors and singular values
+[U,Lambda,V] = svds(C,N);
+
+% PC
+a = S*U;
+b = P*V;
+
+% Make them clear for output
+for iN=1:N
+    e1(iN,:) = squeeze( U(:,iN) )';
+   pc1(iN,:) = squeeze( a(:,iN) )';
+    e2(iN,:) = squeeze( V(:,iN) )';
+   pc2(iN,:) = squeeze( b(:,iN) )';
+end
+
+% Amount of variance explained a 0.1 pres et en %
+L2=Lambda.^2;
+dsum=diag(L2)/trace(L2);
+
+for iN=1:N
+   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
+end
+
+
+disp('====> Finished !')
+
+
+
+
Index: ../trunk-jpl/externalpackages/pcatool/install/dispSVD.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispSVD.m     (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispSVD.m     (revision 21589)
@@ -0,0 +1,86 @@
+% dispSVD(typg,CHP1,CHP2,N,expvar) Display SVDs
+%
+% => Display SVDs
+% typg is plot type you want (see plotm.m).
+% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
+% N is number of SVD to display
+% expvar contained the explained variance of SVDs
+%
+% Rq: it''s using the plotm function
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - December 2004
+%  gmaze@univ-brest.fr
+
+function [] = dispSVD(typg,CHP1,CHP2,N,expvar)
+
+load mapanom2
+global ffile1 ffile2
+
+if length(N)>1
+  mod=N;
+else
+  mod=[1:N];
+end  
+Nx=length(mod);
+% --------------------------------------------------
+% AFFICHAGE
+% --------------------------------------------------
+if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
+   CBAROR='horiz';
+   iw=Nx;jw=2;
+else
+   CBAROR='vert';
+   iw=Nx;jw=2;
+end
+
+
+f1=figur;clf; hold on
+for iN=1:Nx
+     subplot(iw,jw,2*iN-1)
+C = squeeze(real(CHP1(mod(iN),:,:)));
+%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+plotm(C,typg);
+caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+%if(iN==1);cx=caxis;end;caxis(cx);
+colormap(mycolormap(mapanom,21));
+c=colorbar(CBAROR);
+titre1=strcat('SVD',num2str(mod(iN)),':',blanks(1),ffile1);
+titre2=strcat(' (',num2str(expvar(mod(iN))),'%)');
+tt=title(strvcat(strcat(titre1,titre2),char(1)));
+set(tt,'fontweight','bold')
+
+     subplot(iw,jw,2*iN)
+C = squeeze(real(CHP2(mod(iN),:,:)));
+%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+plotm(C,typg);
+caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+%if(iN==1);cx=caxis;end;caxis(cx);
+colormap(mycolormap(mapanom,21));
+c=colorbar(CBAROR);
+titre1=strcat('SVD',num2str(mod(iN)),':',blanks(1),ffile2);
+titre2=strcat(' (',num2str(expvar(mod(iN))),'%)');
+tt=title(strvcat(strcat(titre1,titre2),char(1)));
+set(tt,'fontweight','bold')
+
+end %for iN
+
+
+% $$$ for iN=1:2*N
+% $$$   figure(f1);
+% $$$   subplot(iw,jw,iN);
+% $$$ %  colorbar(CBAROR);
+% $$$ end
+% $$$ refresh
+
+
+%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
+titre2=strcat(' Total explained variance:',num2str(sum(expvar(mod))),'%');
+%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')');
+%titre=strvcat(titre2,titre3);
+suptitle(titre2);
+
+set(f1,'Position',[378 39 313 647]); % Laptop screen
+%set(f1,'Position',[369 55 316 899]); % Desktop screen
+%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));
Index: ../trunk-jpl/externalpackages/pcatool/install/caleof.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/caleof.m      (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/caleof.m      (revision 21589)
@@ -0,0 +1,234 @@
+% [EOFs,PC,EXPVAR] = CALEOF(M,N,METHOD) Compute EOF
+%
+% => Compute the Nth first EOFs of matrix M(TIME,MAP).
+% EOFs is a matrix of the form EOFs(N,MAP), PC is the principal
+% components matrix ie it has the form PC(N,TIME) and EXPVAR is
+% the fraction of total variance "explained" by each EOF ie it has
+% the form EXPVAR(N).
+% Differents method can be used:
+% 1 - The "classic" one, ie we compute eigenvectors of the 
+%     temporal covariance matrix with the eig Matlab function.
+% 2 - A faster "classic" one, same as method 1 but we use the
+%     eigs Matlab function.
+% 3 - We compute eigenvectors by the singular value decomposition,
+%     by used of the svd Matlab function.
+% 4 - Same as method 3 but faster by used of the svds Matlab function
+%
+% See also EIG, EIGS, SVD, SVDS
+%
+% Ref: L. Hartmann: "Objective Analysis" 2002
+% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - 
+%      Analyses of climatic Data" 1997
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - March 2004
+%  Revised July 2006
+%  gmaze@univ-brest.fr
+
+
+function [e,pc,expvar,L] = caleof(M,N,method);
+
+% Get dimensions
+[n p]=size(M);
+
+% Temporal covariance is p*p matrix, that why max EOF computable is p, 
+% so we perform a test on parameter N:
+if(N>p)
+ disp('Warning: N is larger than possible so it''s modified to perform')
+ disp('         EOFs computing...');
+ N = p; 
+end
+
+
+% Eventualy time filtering of data
+if 0==1
+   disp('====> Time filtering...')
+   Fc  = 1/20; Fc2 = 1/1;
+   Fc  = 1/7 ; Fc2 = 1/3;
+   SIGNAL = M(:,1);
+   nj = fix(length(SIGNAL)/10); % Nombre de points du filtre
+   for ipt = 1 : p
+       SIGNAL = M(:,ipt);
+       SIGNALF = lanczos(SIGNAL,Fc2,nj);
+       SIGNALF = SIGNALF - lanczos(SIGNALF,Fc,nj);
+       Y(:,ipt) = SIGNALF;
+   end
+   M = Y;
+end
+
+
+disp('====> Let''go for EOFs and pc computing...')
+switch method
+    case 1 % CLASSIC METHOD
+%================================================================
+% Transform the data matrix in the correct form (map*time) for eig
+M = M';
+
+% Remove the time mean (ie the mean of each rows of M)
+% Rq: detrend remove the mean of columns ie we take M'.
+F = detrend(M','constant')';
+
+% Covariance Matrix (inner product over space = covariance in time)
+R = F * F';
+
+% Eigenanalysis of the covariance matrix R
+[E,L] = eig(R);
+
+% Get PC by projecting eigenvectors on original data
+Z = E'*F;
+
+% Make them clear for output
+for iN=1:N
+    e(iN,:) = squeeze( E(:,p-(iN-1)) )';
+   pc(iN,:) = squeeze( Z(p-(iN-1),:) );
+end
+
+% Amount of explained variance (at 0.1%)
+dsum = diag(L)./trace(L);
+for iN=1:N
+   expvar(iN)=fix((dsum(p-(iN-1))*100/sum(dsum))*10)/10;
+end
+
+% Plots Original field and reconstructed one
+if 0==1
+figure;
+subplot(1,2,1);imagesc(abs(M));title('ORIGINAL');cx=caxis;
+%subplot(1,2,2);imagesc((E*Z));title('RECONSTRUCTED')
+subplot(1,2,2);imagesc(abs(e'*pc));title('RECONSTRUCTED');caxis(cx);
+end
+
+    case 2 % RAPID CLASSIC METHOD 
+%================================================================
+% Remove the time mean of each column
+F = detrend(M,'constant');
+
+% Covariance Matrix
+if n >= p
+   R = F' * F;
+else 
+   R = F * F';
+end
+
+% Eigen analysis of the square covariance matrix
+[E,L] = eigs(R,N);
+if n < p
+  E = F' * E;
+  sq = [sqrt(diag(L))+eps]';
+  sq = sq(ones(1,p),:);
+  E = E ./ sq;
+end
+
+% Get PC by projecting eigenvectors on original data
+if n >= p
+   Z = (F*E)';
+else
+   Z =  E'*F';
+end
+
+
+% Make them clear for output
+for iN=1:N
+    e(iN,:) = squeeze( E(:,iN) )';
+   pc(iN,:) = squeeze( Z(iN,:) );
+end
+
+% Amount of variance explained a 0.1 pres et en %
+dsum=diag(L)./trace(L);
+for iN=1:N
+   expvar(iN)=fix((dsum(iN)*100/sum(dsum))*10)/10;
+end
+
+
+    case 3 % SVD METHOD
+%================================================================
+% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - 
+% Analyses of climatic Data" 1997 => p18
+
+% Assume that M is (time*map) matrix
+[n p]=size(M);
+
+% Remove the mean of each column (ie the time mean in each station records)
+F=detrend(M,'constant');
+
+% Form the covariance matrix:
+R = F'*F;
+
+% Find eigenvectors and singular values
+[C,L,CC] = svd(R);
+% Eigenvectors are in CC and the squared diagonal values of L
+% are the eigenvalues of the temporal covariance matrix R=F'*F
+
+% find the PC corresponding to eigenvalue
+PC = F*CC;
+
+% Make them clear for output
+for iN=1:N
+    e(iN,:) = squeeze( CC(:,iN) )';
+   pc(iN,:) = squeeze( PC(:,iN) )';
+end
+
+if 0
+figure;
+subplot(1,2,1);imagesc(F);title('ORIGINAL');cx=caxis;
+subplot(1,2,2);imagesc(C*L*CC');title('RECONSTRUCTED');caxis(cx);
+end
+
+% Amount of variance explained at 0.1%
+dsum=diag(L)./trace(L);
+if length(dsum)<N % L was not squared
+  dsum = [dsum ;zeros(N-length(dsum),1)];
+end
+for iN = 1 : N
+   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
+end
+
+
+    case 4 % FAST SVD METHOD
+%================================================================
+% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - 
+% Analyses of climatic Data" 1997 => p18
+
+% Assume that M is (time*map) matrix
+[n p]=size(M);
+
+% Remove the mean of each column (ie the time mean in each station records)
+F=detrend(M,'constant');
+
+% Form the covariance matrix:
+R = F' * F;
+
+% Find eigenvectors and singular values
+[C,L,CC,flag] = svds(R,N);
+% Eigenvectors are in CC and the squared diagonal values of L
+% are the eigenvalues of the temporal covariance matrix R=F'*F
+% (Sometimes, CC stops for nul eigenvector, then we need to fill to reach N)
+if size(CC,2)<N
+  CC = [CC  zeros(size(CC,1),N-size(CC,2)+1)];
+end
+
+% find the PC corresponding to eigenvalue
+PC = F*CC;
+% Which is similar to: C*L
+
+% Make them clear for output
+for iN=1:N
+    e(iN,:) = squeeze( CC(:,iN) )';
+   pc(iN,:) = squeeze( PC(:,iN) )';
+end
+
+% Amount of variance explained a 0.1 pres et en %
+dsum=diag(L)./trace(L);
+if length(dsum)<N % L was not squared
+  dsum = [dsum ;zeros(N-length(dsum),1)];
+end
+for iN=1:N
+   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
+end
+
+%figure;
+%subplot(1,2,1);imagesc(M);title('ORIGINAL');cx=caxis;
+%subplot(1,2,2);imagesc((e'*pc)');title('RECONSTRUCTED');caxis(cx);
+
+
+end % switch method
+disp('====> Finished !')
Index: ../trunk-jpl/externalpackages/pcatool/install/dispEOF.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispEOF.m     (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispEOF.m     (revision 21589)
@@ -0,0 +1,130 @@
+% dispEOF(typg,CHP,N,expvar) Display EOFs
+%
+% => Display EOFs
+% typg is plot type you want (see plotm.m).
+% CHP contains EOFs and is of the form: CHP(EOF,X,Y)
+% N is number of EOF to display
+% expvar contained the explained variance of EOFs
+%
+% Rq: it''s using the plotm function
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - March 2004
+%  gmaze@univ-brest.fr
+
+function [] = dispEOF(typg,CHP,N,expvar)
+
+load mapanom2
+global ffile
+% --------------------------------------------------
+% AFFICHAGE
+% --------------------------------------------------
+if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
+
+   CBAROR='horiz';
+if N>12,NN=12;else,NN=N;end
+switch NN
+  case 1,tblmagique=[1];                         iw=1;jw=1;
+  case 2,tblmagique=[1 2];                       iw=2;jw=1;
+  case 3,tblmagique=[1 2 3];                     iw=3;jw=1;
+  case 4,tblmagique=[1 2 3 4];                   iw=2;jw=2;
+  case 5,tblmagique=[1 2 3 4 5];                 iw=3;jw=2;
+  case 6,tblmagique=[1 2 3 4 5 6];               iw=3;jw=2;
+  case 7,tblmagique=[1 2 3 4 5 6 7];             iw=4;jw=2;
+  case 8,tblmagique=[1 2 3 4 5 6 7 8];           iw=4;jw=2;
+  case 9,tblmagique=[1 2 3 4 5 6 7 8 9];         iw=3;jw=3;
+ case 10,tblmagique=[1 2 3 4 5 6 7 8 9 10];      iw=4;jw=3;
+ case 11,tblmagique=[1 2 3 4 5 6 7 8 9 10 11];   iw=4;jw=3;
+ case 12,tblmagique=[1 2 3 4 5 6 7 8 9 10 11 12];iw=4;jw=3;
+end
+iw=N;jw=1;
+
+else
+
+   CBAROR='vert';
+if N>12,NN=12;else,NN=N;end
+switch NN
+  case 1,tblmagique=[1];                         iw=1;jw=1;
+  case 2,tblmagique=[1 2];                       iw=2;jw=1;
+  case 3,tblmagique=[1 2 3];                     iw=3;jw=1;
+  case 4,tblmagique=[1 2 3 4];                   iw=2;jw=2;
+  case 5,tblmagique=[1 2 3 4 5];                 iw=3;jw=2;
+  case 6,tblmagique=[1 2 3 4 5 6];               iw=3;jw=2;
+  case 7,tblmagique=[1 2 3 4 5 6 7];             iw=4;jw=2;
+  case 8,tblmagique=[1 2 3 4 5 6 7 8];           iw=4;jw=2;
+  case 9,tblmagique=[1 2 3 4 5 6 7 8 9];         iw=3;jw=3;
+ case 10,tblmagique=[1 2 3 4 5 6 7 8 9 10];      iw=4;jw=3;
+ case 11,tblmagique=[1 2 3 4 5 6 7 8 9 10 11];   iw=4;jw=3;
+ case 12,tblmagique=[1 2 3 4 5 6 7 8 9 10 11 12];iw=4;jw=3;
+end
+
+
+
+end %if
+
+
+
+f1=figur;clf; hold on
+
+for iN=1:NN
+     subplot(iw,jw,tblmagique(iN)); hold on
+C = squeeze(real(CHP(iN,:,:)));
+%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+if(iN==1);cx=caxis;end;caxis(cx);
+colormap(mycolormap(mapanom,21));
+
+%c=colorbar(CBAROR);
+titre1=strcat('EOF',num2str(iN));
+titre2=strcat(' (',num2str(expvar(iN)),'%)');
+title(strcat(titre1,titre2));
+
+end %for iN
+
+%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
+titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
+%titre3=strcat('(Analysed file: ',ffile,')');
+%titre=strvcat(titre2,titre3);
+suptitle(titre2);
+
+set(f1,'Position',[378 39 313 647]); % Laptop screen
+%set(f1,'Position',[369 55 316 899]); % Desktop screen
+%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));
+set(f1,'Name','<> REAL PART <>');
+
+
+% -------------------------------------------
+if isreal(CHP)==0
+
+f2=figur;clf; hold on
+
+for iN=1:NN
+     subplot(iw,jw,tblmagique(iN)); hold on
+
+C = squeeze(imag(CHP(iN,:,:)));
+plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+%if(iN==1);cx=caxis;end;
+caxis(cx);
+colormap(mycolormap(mapanom,21));
+
+c=colorbar(CBAROR);
+titre1=strcat('EOF',num2str(iN));
+titre2=strcat(' (',num2str(expvar(iN)),'%)');
+title(strcat(titre1,titre2));
+
+end %for iNN
+
+%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
+titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
+%titre3=strcat('(Analysed file: ',ffile,')');
+%titre=strvcat(titre2,titre3);
+suptitle(titre2);
+
+set(f2,'Position',[378 39 313 647]); % Laptop screen
+%set(f2,'Position',[369 55 316 899]); % Desktop screen
+set(f2,'Name','<> IMAGINARY PART <>');
+
+end
+
+
Index: ../trunk-jpl/externalpackages/pcatool/install/plotm.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/plotm.m       (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/plotm.m       (revision 21589)
@@ -0,0 +1,214 @@
+%PLOTM  Trace un champ de donnee geographique
+% [CS,H]=PLOTM(C,[ZONE]) 
+%   Trace le champ C(LONGITUDE,LATITUDE) sur une carte 
+%   du monde en projection Mercator par defaut.
+%   Le parametre optionnel ZONE peut prendre les valeurs
+%     ZONE = 1 : Projection Stereographique centree sur le pole sud
+%     ZONE = 2 : Projection Mercator de l'hemisphere sud
+%     ZONE = 3 : Projection Mercator du globe
+%     ZONE = 4 : Pas de projection (imagesc)
+%     ZONE = 5 : Projection Mercator de l'hemisphere sud complet
+%     ZONE = 6 : Projection Mercator du globe complet
+%     ZONE = 7 : Projection Stereographique centree sur le pole nord
+%     ZONE = 8 : Idem 2 mais sans label sur les axes
+%     ZONE = 9 : Projection Mercator de l'hemisphere sud aux limites du Fmask
+%   Returns contour matrix CS as described in CONTOURC and a vector H
+%   of handles to PATCH objects (for use by CLABEL).
+%
+%   La definition des grilles (longitude,latitude) est celle du
+%   modele par defaut, mais elle peut etre modifiee en passant en
+%   global les nouvelles variables nommees lon et lat.
+%
+% NECESSITE LE PACKAGE M_MAP
+%
+%
+% (nov. 2003 -> gmaze@univ-brest.fr)
+
+function varargout = plotm(C,varargin)
+
+if (nargin>2)|(nargin<1)
+     help plotm.m
+     error('plotm.m : Wrong number or bad parameter(s)')
+     return
+elseif (nargin==2)
+   arg = varargin{:};
+   if (arg>14)|(arg<1)
+     help plotm.m
+     error('plotm.m : Wrong number or bad parameter(s)')
+     return
+   end
+end %if
+
+%--------------------------------------------------
+% Variabes
+%--------------------------------------------------
+[lx,ly]=size(C);
+% Si lon et lat ne sont pas definie comme dans le modele
+% on prend en global n'importe quel definition
+if ((lx~=64)|(ly~=32)) 
+   global lon lat
+else
+   lon= linspace(0.,360.,64);
+   lat= [85.76 80.27 74.74 69.21 63.68 58.14 52.61 47.07 41.53 ...
+      36.00 30.46 24.92 19.38 13.84 8.31 2.77 -2.77 -8.31 ...
+      -13.84 -19.38 -24.92 -30.46 -36.00 -41.53 -47.07 -52.61 ...
+      -58.14 -63.68 -69.21 -74.74 -80.27 -85.76];
+end
+
+% Le nombre de contour est par defaut a 10, mais on peut le
+% modifier en declarant nbcont en global
+global nbcont
+if isempty(nbcont)==1
+   nbcont=10;
+end
+
+%-----------------------------------------------
+% Projection
+%-----------------------------------------------
+if (nargin==2)
+   arg = varargin{:};
+   switch arg(1)
+     case 1 % Projection stereo 
+       m_proj('stereo','lon',0,'lat',-90,'rad',71);
+     case 2 % Projection mercator (HS)
+       m_proj('miller','lon',[0 360],'lat',[-80 -20]);
+       m_proj('miller','lon',[0 360],'lat',[-60 -20]);
+     case 3 % Projection mercator (Tout)  -> Defaut
+       m_proj('miller','lon',[0 360],'lat',[-80 80]);
+     case 5 % Projection mercator (HS-complet)
+       m_proj('miller','lon',[0 360],'lat',[-85 -5]);
+     case 6 % Projection mercator (Tout-complet)
+       m_proj('miller','lon',[0 360],'lat',[-85 85]);
+     case 7 % Projection stereo 
+       m_proj('stereo','lon',0,'lat',90,'rad',71);
+     case 8 % Projection mercator (HS) sans Label sur les axes
+       m_proj('miller','lon',[0 360],'lat',[-80 -20]);
+       m_proj('miller','lon',[0 360],'lat',[-75 -30]);
+     case 9 % Projection mercator (HS) sans Label sur les axes aux limtes du FMASK
+       m_proj('miller','lon',[0 360],'lat',[-69 -25]);
+     case 10 % Projection stereo recentree sur le pole sud
+       m_proj('stereo','lon',0,'lat',-90,'rad',51);
+     case 11 % Indien
+       m_proj('lambert','long',[11 125],'lat',[-60 -30]);
+     case 12 % Pacifique
+       m_proj('lambert','long',[120 290],'lat',[-60 -30]);
+     case 13 % Atlantique
+       m_proj('lambert','long',[-70 20],'lat',[-60 -30]);
+     case 14 % Centré Atlantique        
+       m_proj('lambert','long',[-90 90],'lat',[-70 -20],'rectbox','on');
+   end
+else
+% Projection mercator (Tout)  -> Defaut
+  m_proj('mercator','lon',[0 360],'lat',[-80 80]);
+end
+
+%------------------------------------------------
+% Affichage
+%------------------------------------------------
+load MODELCOAST
+if length(lon)==65
+  modelcoast=[modelcoast;modelcoast(1,:)];
+end
+
+hold on
+if (nargin==2)
+   arg = varargin{:};
+   switch arg(1)
+     case 1 % Projection stereo (Par defaut)
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-30 -40 -50 -60 -70],'yticklabels',[],...
+             'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top')
+     case 2 % Projection mercator (HS)
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-80 -70 -60 -50 -40 -30 -20]);
+     case 3 % Projection mercator (Tout)
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]);
+     case 4 % Pas de projection
+      pcolor(lon,lat,C');
+      axis tight;box on
+%      axis([lon(1) lon(64) lat(32) lat(1)]);
+%      gc=gca;set(gc,'YDir','reverse');
+     case 5 % Projection mercator (Hs-complet)
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      COASTYPE='color';COASTCOLOR=[0 0 0];
+%      COASTYPE='patch';COASTCOLOR=[0 0.5 0.6];
+%      m_coast(COASTYPE,COASTCOLOR);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-80 -60 -40 -20 0]);
+     case 6 % Projection mercator (Tout-complet)
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]);
+     case 7 % Projection stereo
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[30 40 50 60 70],'yticklabels',[],...
+             'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top')
+     case 8 % Projection mercator (HS) sans label sur les axes
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_pcolor(lon,lat,C');shading interp
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick','','xtick','');
+     case 9 % Projection mercator (HS) sans label sur les axes aux limtes du FMASK
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick','','xtick','');
+     case 10 % Projection stereo recentree sur le pole Sud
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+%      m_coast('color',[0 0 0]);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
+             'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top')
+     case 11 % Indien
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
+             'xtick',[30 60 90 120 150 180],'XaxisLocation','top')
+     case 12 % Pacifique
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
+             'xtick',[120 150 180 210 240 270 300],'XaxisLocation','top')
+     case 13 % Atlantique
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
+             'xtick',[-60 -30 0 30 60],'XaxisLocation','top')
+     case 14 % Centré Atlantique
+      [cs,h]=m_contourf(lon,lat,C',nbcont);
+      m_contour(lon,lat,modelcoast',1,'k');
+      m_grid('yticklabels',[],'XaxisLocation','top',...
+             'fontsize',6,...
+             'xtick',[0 60 120 180 -120 -60],'xticklabels',[180 -120 -60 0 60 120]);
+            
+   end % switch
+else
+  %  Projection Mercator (Par defaut)
+     [cs,h]=m_contourf(lon,lat,C',nbcont);
+%     m_coast('color',[0 0 0]);
+     m_contour(lon,lat,modelcoast',1,'k');
+     m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]);
+end
+
+
+
+% OUTPUT VARIABLES
+switch nargout
+  case 1
+   varargout(1) = {cs} ;
+  case 2
+   varargout(1) = {cs} ;
+   varargout(2) = {h};
+end
Index: ../trunk-jpl/externalpackages/pcatool/install/calCeof.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/calCeof.m     (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/calCeof.m     (revision 21589)
@@ -0,0 +1,11 @@
+% [e,pc,expvar] = calCeof(M,N,METHOD)
+%   Compute Complex EOF
+
+    function [e,pc,expvar] = calCeof(M,N,METHOD)
+
+
+% Compute Hilbert transform Mh and create M as: M = M + i*Mh
+M = hilbert(M);
+
+% Compute EOFs
+ [e,pc,expvar] = caleof(M,N,METHOD);
Index: ../trunk-jpl/externalpackages/pcatool/install/dispCEOF.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispCEOF.m    (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispCEOF.m    (revision 21589)
@@ -0,0 +1,65 @@
+% dispCEOF(CEOF,EXPVAR,PHI,MOD) Display few CEOFs.
+%
+% => DISPLAY FEW CEOFs.
+% CEOF contains all the CEOFs as nCEOF*X*Y.
+% EXPVAR is a matrix with the explained variance of each
+%  CEOFs in %. This is just for title.
+% PHI is phase when ploting maps.
+% MOD contains explicitly the CEOFs to display.
+%
+% Rq: it's using the plotm function
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - June 2004
+%  gmaze@univ-brest.fr
+
+function [] = dispCEOF(CEOF,EXPVAR,PHI,MOD)
+
+% Number of EEOF to display     
+NMOD = length(MOD);
+NANG = length(PHI);
+
+% Open figure and first guest
+figur;clf;hold on
+typg=8;
+width  = .9/NMOD;
+height = .9/NANG;
+dleft = (.95-width*NMOD)/4;
+left  = dleft;
+load mapanom2
+
+% Let'go :
+
+% We select a CEOF:
+for imod = 1 : NMOD
+    mod = MOD(imod);
+    chp=squeeze(CEOF(mod,:,:));
+    Rchp=real(chp);
+    Ichp=-imag(chp);
+
+    % 
+    for ang=1:NANG
+      C = cos(PHI(ang))*Rchp + sin(PHI(ang))*Ichp;
+      m(ang)=abs(xtrm(C));
+    end
+    m=max(m);
+    
+    % and plot maps
+    for ang=1:NANG
+      C = cos(PHI(ang))*Rchp + sin(PHI(ang))*Ichp;     
+      bottom = .95-height*ang;
+      subplot('Position',[left bottom width height]);
+      plotm(C,typg);
+      m_text(10,-70,strcat('\phi=',num2str(PHI(ang)*180/pi),'^o'));
+%      ylabel(strcat('\phi=',num2str(PHI(ang)*180/pi),'^o'),'Rotation',0);
+      caxis([-m m]);
+      colormap(mapanom);
+      if(ang==1)
+         titre=strcat('CEOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)');
+         title(titre);
+      end;
+    end
+    left = left + width + dleft;
+
+end
Index: ../trunk-jpl/externalpackages/pcatool/install/dispSVDnb.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispSVDnb.m   (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispSVDnb.m   (revision 21589)
@@ -0,0 +1,63 @@
+% dispSVDnb(typg,CHP1,CHP2,N,expvar) Display SVDs
+%
+% => Display SVDs in white and black
+% typg is plot type you want (see plotm.m).
+% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
+% N is number of SVD to display
+% expvar contained the explained variance of SVDs
+%
+% Rq: it''s using the plotmNB function
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - December 2004
+%  gmaze@univ-brest.fr
+
+function [] = dispSVDnb(typg,CHP1,CHP2,N,expvar)
+
+global ffile1 ffile2
+% --------------------------------------------------
+% AFFICHAGE
+% --------------------------------------------------
+if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
+   CBAROR='horiz';
+   iw=N;jw=2;
+else
+   CBAROR='vert';
+   iw=N;jw=2;
+end
+
+
+f1=figur;clf; hold on
+for iN=1:N
+     subplot(iw,jw,2*iN-1)
+C = squeeze(real(CHP1(iN,:,:)));
+%[cs,h]=plotmNB(C,typg); whos cs h
+[cs,h]=plotmNB(C,typg);
+%clabel(cs,h,'labelspacing',100);
+titre1=strcat('SVD',num2str(iN),':',ffile1);
+titre2=strcat(' (',num2str(expvar(iN)),'%)');
+title(strcat(titre1,titre2));
+
+     subplot(iw,jw,2*iN)
+C = squeeze(real(CHP2(iN,:,:)));
+[cs,h]=plotmNB(C,typg);
+%clabel(cs,h,'labelspacing',100);
+titre1=strcat('SVD',num2str(iN),':',ffile2);
+titre2=strcat(' (',num2str(expvar(iN)),'%)');
+title(strcat(titre1,titre2));
+
+end %for iN
+
+
+
+
+%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
+titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
+%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')');
+%titre=strvcat(titre2,titre3);
+suptitle(titre2);
+
+set(f1,'Position',[378 39 313 647]); % Laptop screen
+%set(f1,'Position',[369 55 316 899]); % Desktop screen
+%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));
Index: ../trunk-jpl/externalpackages/pcatool/install/Contents.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/Contents.m    (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/Contents.m    (revision 21589)
@@ -0,0 +1,32 @@
+% PRINCIPAL COMPONENTS ANALYSIS TOOLBOX
+% Last Modified: Dec 2007
+%
+% This software is provided "as is" without express 
+% or implied warranty...
+%
+%
+% EOF and SVD 
+%    caleof       - Compute EOFs
+%    calEeof      - Compute Extended EOFs
+%    calCeof      - Compute Complex  EOFs
+%    calsvd       - Compute SVDs
+%    confexpvar   - Compute statistical confidence of explained variance
+%
+% VISUALISATION
+%    dispEOF      - Display EOFs
+%    dispEEOF     - Display EEOFs as time sequence map
+%    dispEEOF_v2  - Display EEOFs as time sequence map with
+%                   the PC time series
+%    dispCEOF     - Display CEOFs
+%    dispSVD      - Display SVDs
+%    plotm        - Plot a 2D geo field
+%
+% Matrix manipulation
+%    map2mat      - Transform C(TIME,LON,LAT) in D(TIME,PT)
+%                   under a geographic points selection
+%    mat2map      - Invert transform done by map2mat
+%
+%
+%================================================================
+% Author: Guillaume MAZE - LPO/LMD - March 2004
+%         gmaze@univ-brest.fr
Index: ../trunk-jpl/externalpackages/pcatool/install/calEeof.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/calEeof.m     (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/calEeof.m     (revision 21589)
@@ -0,0 +1,128 @@
+% [EOFs,PC,EXPVAR] = calEeof(M,N,METHOD,NLAG,DT) Compute EEOF
+%
+% => Compute N Extended EOF of a M(TIME*MAP) matrix.
+% METHOD is the EOF method computing (see caleof.m)
+% NLAG is the number of map you want
+% DT is the increment between each map (not the time step
+% of your data)
+%
+% The function can filter datas with a Lanczos band-pass 
+% filter. Edit the calEeof.m file and change option in it. 
+%
+% See also: CALEOF
+%
+% Ref: Weare and Nasstrom 1982
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - March 2004
+%  gmaze@univ-brest.fr
+
+    function varargout = calEeof(M,N,METHOD,NLAG,DT)
+
+% ---------------------------------
+% Divers
+tic
+LAG = NLAG*DT;
+foll = 1; % Turn this value to 0 if you don t want commentary along computing
+
+% ---------------------------------
+% LAG must be >0, LAG=1 correspond to initial position
+if LAG<=0
+   LAG = 1;
+end
+
+% ---------------------------------
+% Usely for us, M is a TIME*MAP matrix so:
+% We get dimensions
+  [n p]=size(M);
+% and put M in the correct form for computing
+  M = M';
+
+% ---------------------------------
+% Eventualy time filtering of data
+if 0==1
+   if (foll~=0),disp('==> Time filtering...');end
+   dt = 1;         % Real time step of your data (diff from DT)
+   Fc  = 1/dt/8;   % This the low cut off frequency
+   Fc2 = 1/dt/2;   % and this the high cut off frequency
+   SIGNAL = M(1,:);
+   nj = fix(length(SIGNAL)/10); % Nb of points of the window
+
+   for ipt = 1 : p
+       SIGNAL = M(ipt,:)';
+       SIGNALF = lanczos(SIGNAL,Fc2,nj);
+       SIGNALF = SIGNALF - lanczos(SIGNALF,Fc,nj);
+       Y(:,ipt) = SIGNALF;
+       if mod(ipt,10)==0 % We display a plot of filtered data
+          clf;
+          plot((0:n-1),SIGNAL,'k',(0:n-1),SIGNALF,'r');          
+          drawnow;
+       end
+   end
+   M = Y';
+end, clear Fc Fc2 SIGNAL nj SIGNALF Y ipt, close
+
+
+% ---------------------------------
+% This is the matrix where we concat all submatrices from M.
+% Size of each submatrix is : NLAG*p rows vs n-LAG+1 colums 
+F = zeros(NLAG*p,n-LAG+1);
+
+% ---------------------------------
+% Let's go for F:
+if (foll~=0),disp('==> Forming concatenated matrix...');end
+if DT==1
+% CLASSIC CASE
+for ilag = 1 : LAG
+       % Extract submatrix ilag from M
+       Msub = M(:,ilag:n-LAG+ilag);
+       % Concat matrix
+       F( (ilag-1)*p+1 : ilag*p , : ) = Msub;
+end
+else
+% DT>1
+it=0;
+for ilag = 1 : DT : LAG
+       % Extract submatrix ilag from M
+         Msub = M( : , ilag : n-LAG+ilag);
+       % Concat matrix
+         it = it + 1;
+         F( (it-1)*p+1 : it*p , : ) = Msub;
+         % imagesc(Msub);pause % for debugging
+end
+end
+
+% ---------------------------------
+% Compute EOFs by normal way
+% (Don't forget that caleof is taking a TIME*MAP argument matrix,
+% that's why we have F' !)
+if (foll~=0),disp('==> Computing EOFs ...');end
+[e,pc,expvar] = caleof(F',N,METHOD);
+
+% ---------------------------------
+% e is the matrix with EOFs inside. We should extract each map
+% for different lags from it.
+% e is NEOF*(LAG*MAP) and we construct 3D matrix CHP as NEOF*LAG*MAP
+if (foll~=0),disp('==> Reshaping vectors ...');end
+for ilag = 1 : NLAG
+    E = e( : , (ilag-1)*p+1 : ilag*p );
+    CHP(:,ilag,:) = E;
+end
+
+
+% ---------------------------------
+% OUTPUT VARIABLES
+switch nargout
+  case 1
+   varargout(1) = {CHP} ;
+  case 2
+   varargout(1) = {CHP} ;
+   varargout(2) = {pc};
+  case 3
+   varargout(1) = {CHP} ;
+   varargout(2) = {pc};
+   varargout(3) = {expvar};
+end
+
+% ---------------------------------
+if (foll~=0),toc,disp('==> That''s all folks !');end
Index: ../trunk-jpl/externalpackages/pcatool/install/calsvd2.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/calsvd2.m     (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/calsvd2.m     (revision 21589)
@@ -0,0 +1,49 @@
+% [SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = CALSVD2(A,B,N) Compute SVDs
+%
+% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - 
+%      Analyses of climatic Data" 1997
+%================================================================
+%
+%  Guillaume MAZE - LPO/LMD - March 2004
+%  gmaze@univ-brest.fr
+
+
+   function [e1,e2,pc1,pc2,expvar,Lambda,dsumCF] = calsvd2(A,B,N);
+
+
+%================================================================
+% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - 
+% Analyses of climatic Data" 1997 => p18
+
+% Assume that A is (time*map) matrix
+[n p]=size(A);
+
+% Remove the mean of each column (ie the time mean in each station records)
+S=detrend(A,'constant');
+P=detrend(B,'constant');
+
+% Form the covariance matrix:
+C=S'*P;
+
+% Find eigenvectors and singular values
+[U,Lambda,V] = svds(C,N);
+
+% PC
+a=S*U;
+b=P*V;
+
+% Make them clear for output
+for iN=1:N
+    e1(iN,:) = squeeze( U(:,iN) )';
+   pc1(iN,:) = squeeze( a(:,iN) )';
+    e2(iN,:) = squeeze( V(:,iN) )';
+   pc2(iN,:) = squeeze( b(:,iN) )';
+end
+
+% Amount of variance explained a 0.1 pres et en %
+L2=Lambda.^2;
+dsum=diag(L2)/trace(L2);
+for iN=1:N
+   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
+end
+
Index: ../trunk-jpl/externalpackages/pcatool/install/dispEEOF.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispEEOF.m    (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispEEOF.m    (revision 21589)
@@ -0,0 +1,57 @@
+% dispEEOF(CHP,EXPVAR,DT,NLAG,MOD) Display few EEOFs.
+%
+% => DISPLAY FEW EEOFs.
+% CHP contains all the EEOFs as EOF*LAG*X*Y.
+% EXPVAR is a matrix with the explained variance of each
+%  EEOFs in %. This is just for title.
+% DT is time step between maps.
+% NLAG is the number of LAG to display.
+% MOD contains explicitly the EEOFs to display.
+%
+% Rq: it's using the plotm function
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - March 2004
+%  gmaze@univ-brest.fr
+
+function [] = dispEEOF(CHP,EXPVAR,DT,NLAG,MOD)
+
+% Number of EEOF to display     
+NMOD = length(MOD);
+
+% Open figure and first guest
+figur;
+clf;hold on
+typg=8;
+width  = .9/NMOD;
+height = .9/NLAG;
+dleft = (.95-width*NMOD)/4;
+left  = dleft;
+load mapanom
+
+% Let'go :
+
+% We choose an EEOF
+for imod = 1 : NMOD
+    mod = MOD(imod);
+    
+    % and plot maps
+    for lag=1:NLAG
+      bottom = .95-height*lag;
+      subplot('Position',[left bottom width height]);
+      C = squeeze(real(CHP(mod,lag,:,:)));
+%      C = C./xtrm(C); % Eventually normalise field
+      plotm(C,typg);
+      if(lag==1),caxis([-abs(xtrm(C)) abs(xtrm(C))]);cx=caxis;end;
+      caxis(cx);
+      colormap(mapanom);
+      if(lag==1)
+         titre=strcat('EOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)');
+         titre=strcat(titre,'; DT=',num2str(DT),' ');
+         title(titre);
+      end;
+    end
+    left = left + width + dleft;
+
+end
Index: ../trunk-jpl/externalpackages/pcatool/install/dispSVDnbc.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispSVDnbc.m  (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispSVDnbc.m  (revision 21589)
@@ -0,0 +1,66 @@
+% dispSVDnb(typg,CHP1,CHP2,N,expvar) Display SVDs
+%
+% => Display SVDs in white and black with correlation form
+% typg is plot type you want (see plotm.m).
+% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
+% N is number of SVD to display
+% expvar contained the explained variance of SVDs
+%
+% Rq: it''s using the plotmNB function
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - December 2004
+%  gmaze@univ-brest.fr
+
+function [] = dispSVDnbc(typg,CHP1,CHP2,N,expvar)
+
+global ffile1 ffile2 nbcont
+nbcont=[-1:0.2:1];
+nbcont2=[-1:0.4:1];
+
+% --------------------------------------------------
+% AFFICHAGE
+% --------------------------------------------------
+if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
+   CBAROR='horiz';
+   iw=N;jw=2;
+else
+   CBAROR='vert';
+   iw=N;jw=2;
+end
+
+
+f1=figur;clf; hold on
+for iN=1:N
+     subplot(iw,jw,2*iN-1)
+C = squeeze(real(CHP1(iN,:,:)));
+[cs,h]=plotmNB(C,typg);
+clabel(cs,h,nbcont2,'labelspacing',100,'fontsize',9);
+titre1=strcat('SVD',num2str(iN),':',ffile1);
+titre2=strcat(' (',num2str(expvar(iN)),'%)');
+title(strcat(titre1,titre2));
+
+     subplot(iw,jw,2*iN)
+C = squeeze(real(CHP2(iN,:,:)));
+[cs,h]=plotmNB(C,typg);
+clabel(cs,h,nbcont2,'labelspacing',100,'fontsize',9);
+%clabel(cs,h,'labelspacing',100);
+titre1=strcat('SVD',num2str(iN),':',ffile2);
+titre2=strcat(' (',num2str(expvar(iN)),'%)');
+title(strcat(titre1,titre2));
+
+end %for iN
+
+
+
+
+%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
+titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
+%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')');
+%titre=strvcat(titre2,titre3);
+suptitle(titre2);
+
+set(f1,'Position',[378 39 313 647]); % Laptop screen
+%set(f1,'Position',[369 55 316 899]); % Desktop screen
+%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));
Index: ../trunk-jpl/externalpackages/pcatool/install/test_pcatool.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/test_pcatool.m        (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/test_pcatool.m        (revision 21589)
@@ -0,0 +1,160 @@
+% Here is a simple script to test caleof.m from PCATOOL package.
+% 
+
+% gmaze@mit.edu
+% 2006/07/31
+%
+
+clear
+
+field = 1; % The best one
+
+% Generate the data
+x = -10 : 1 : 10 ;
+y = x;
+[X,Y] = meshgrid(x,y);
+
+switch field
+ case 1   % H. Bjornson and S.A. Venegas
+   HH = 12*( 1.2 - 0.35*sqrt( (X-0).^2 + (Y-0).^2 ) );
+   HH = HH - mean(mean(HH));
+   H(1,:,:) = 1012 + HH;
+   H(2,:,:) = 1012 + (1012-squeeze(H(1,:,:)));;
+   H(3,:,:) = 1022.8 - 3.6*Y;
+   H(4,:,:) = 1012 + (1012-squeeze(H(3,:,:)));
+   H(5,:,:) = 1001.2 + 3.6*X;
+   H(6,:,:) = 1012 + (1012-squeeze(H(5,:,:)));
+  cont = [980:5:1050];
+  N = 2;
+  
+ case 11   % H. Bjornson and S.A. Venegas modified
+   % Here it must be only 2 non-zero EOF (45degre oriented ) of similar variance 50%
+   H = ones(4,length(x),length(y));
+   H(1,:,:) =  Y;
+   H(2,:,:) = -Y;
+   H(3,:,:) = X;
+   H(4,:,:) = -X;
+   N = 2;
+
+ case 2 % Hartmann eg:  analysis in place
+  a = [2 4 -6 8 ; 1 2 -3 4];
+  H(1,:,:) = a(1,:)';
+  H(2,:,:) = a(2,:)';
+  
+  % Get EOFs:
+  N = 4;
+  for method = 2 : -1 : 1
+    [E,pc,expvar] = caleof(H,N,method);
+    %E'
+  end
+
+  return
+  
+ case 3  % My signal
+  np = 4; % Numer of random signal (...and eof)
+  nt = 10; % Number of timestep
+  H = zeros(nt,length(x),length(y));
+  for ip = 1 : np
+    xc(ip) = abs(fix(rand(1)*10));
+    yc(ip) = abs(fix(rand(1)*10));
+    if ip>=fix(np/2)
+      xc(ip) = -xc(ip);
+      yc(ip) = -yc(ip);
+    end
+    dd(ip) = fix(rand(1)*10);
+  end %for ip
+  f = 1/nt;
+  for it = 1 : nt 
+    H2 = zeros(length(x),length(y));
+    for ip = 1 : ip
+        if it==1,[xc(ip) yc(ip) dd(ip)],end
+        HH = 12*( 1.2 - 0.35*sqrt(  ((X-xc(ip)).^2 + (Y-yc(ip)).^2 )/dd(ip) ));
+        H2 = HH - mean(mean(HH));
+       H(it,:,:) = squeeze(H(it,:,:)) + H2.*cos(pi*it/dd(ip)) ;
+    end %for ip
+  end %for it
+  H = 1012 + H;
+  cont = [980:5:1050];
+  N = 3;
+  
+ case 4 % My signal 2
+  x = -pi : pi/6 : pi ;
+  y = x;
+  [X,Y] = meshgrid(x,y);
+  HH  = cos(X) + cos(Y);
+  HH2 = cos(Y);
+  nt = 12;
+  for it = 1 : nt
+%     H(it,:,:) = cos(pi*it/nt)*HH; cont=[-2 2];
+     H(it,:,:) = 2*cos(pi*it/nt)*HH + 3*cos(pi*it/nt/2)*HH2; cont=[-10 10];
+     xtrm(squeeze(H(it,:,:)));
+  end
+  cont=[-2 2];
+  N = 3;
+  
+end % switch field
+
+%return
+
+% Plot field time serie:
+if 1
+figure;iw=2;jw=size(H,1)/iw;
+set(gcf,'position',[11 533 560 420]);
+
+for i=1:iw*jw
+  C = squeeze(H(i,:,:));
+  
+  subplot(iw,jw,i);
+  pcolor(X,Y,C);
+  title(num2str(i));
+  if i==1,cx=caxis;end
+  axis square
+  caxis(cx);
+
+end %for i
+end %fi
+%return
+
+% Get EOFs:
+G = map2mat(ones(size(H,2),size(H,3)),H);
+
+for method = 1 : 4
+  [E,pc,expvar] = caleof(G,N,method);  
+  eof = mat2map(ones(size(H,2),size(H,3)),E);
+  
+
+figure;iw=1;jw=N+1;
+set(gcf,'MenuBar','none');
+posi = [576 0 710 205];
+set(gcf,'position',[posi(1) (method-1)*250+50 posi(3) posi(4)]);
+
+for i=1:iw*jw
+  if i<= iw*jw-1
+    
+  C = squeeze(eof(i,:,:));
+  subplot(iw,jw,i);
+  cont = 12;
+  [cs,h] = contourf(X,Y,C,cont);
+  clabel(cs,h); 
+  title(strcat('EOF:',num2str(i),'/',num2str(expvar(i)),'%'));
+  axis square;
+  %caxis([cont(1) cont(end)]);
+  
+  else
+  subplot(iw,jw,iw*jw);
+  plot(pc');
+  grid on
+  xlabel('time')
+  title('PC')
+  legend(num2str([1:N]'),2);
+  box on
+  
+
+  end %if
+  
+end %for i
+%suptitle(strcat('METHODE:',num2str(method)));
+
+end %for method
+
+
Index: ../trunk-jpl/externalpackages/pcatool/install/dispSVDlag.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/dispSVDlag.m  (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/dispSVDlag.m  (revision 21589)
@@ -0,0 +1,62 @@
+% dispSVDlag(typg,CHP1,CHP2,N,expvar,LAGS,dt) Display SVDs lagged
+%
+% => Display SVDs lagged
+% typg is plot type you want (see plotm.m).
+% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
+% N is number of SVD to display
+% expvar contained the explained variance of SVDs
+%
+% Rq: it''s using the plotm function
+%
+%================================================================
+
+%  Guillaume MAZE - LPO/LMD - December 2004
+%  gmaze@univ-brest.fr
+
+   function [] = dispSVD(typg,CHP1,CHP2,N,expvar,LAGS,dt)
+
+load mapanom2
+global ffile1 ffile2
+% --------------------------------------------------
+% AFFICHAGE
+% --------------------------------------------------
+if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
+   CBAROR='horiz';
+else
+   CBAROR='vert';
+end
+
+   CHP1=squeeze(CHP1(:,N,:,:));
+   CHP2=squeeze(CHP2(:,N,:,:));
+
+Nlag=size(CHP1,1);
+
+iw=round(Nlag/dt);jw=2;
+
+
+f1=figur;clf; hold on
+ilag=0;
+for lag=1:dt:Nlag
+  ilag=ilag+1;
+
+    subplot(iw,jw,2*ilag-1)
+C = squeeze(real(CHP1(lag,:,:)));
+plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+if(ilag==1);cx=caxis;end;caxis(cx);
+colormap(mycolormap(mapanom,21));
+c=colorbar(CBAROR);
+titre1=strcat('lag:',num2str(LAGS(lag)),'y: ',ffile1);
+title(titre1);
+
+     subplot(iw,jw,2*ilag)
+C = squeeze(real(CHP2(lag,:,:)));
+plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
+if(ilag==1);cx=caxis;end;caxis(cx);
+colormap(mycolormap(mapanom,21));
+c=colorbar(CBAROR);
+titre1=strcat('lag:',num2str(LAGS(lag)),'y: ',ffile2);
+title(titre1);
+
+end %for iN
+
+suptitle(strcat('Lagged SVD, mode:',num2str(N)));
Index: ../trunk-jpl/externalpackages/pcatool/install/map2mat.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/map2mat.m     (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/map2mat.m     (revision 21589)
@@ -0,0 +1,40 @@
+% D = MAP2MAT(F,C) Reshaping matrix
+%
+% => Conversion of a 'map' matrix C(TIME,LON,LAT) into a D(TIME,PT) matrix
+% under the mask F(LON,LAT).
+% F is a matrix contenaing 1 where you would like to keep the point and
+%  0 elsewhere.
+%
+% Rq: No check is done about the input.
+%
+% See also: mat2map
+%================================================================
+
+% March 2004
+% gmaze@univ-brest.fr
+
+function [D] = map2mat(F,C);
+
+% Get dimensions
+[tps nolon nolat] = size(C);
+
+% So output matrix will be:
+D = zeros(tps,nolon*nolat);
+
+% point indice
+ipt = 0;
+
+% 'Un-mapping' :
+for iy=1:nolat
+  for ix=1:nolon
+      if F(ix,iy)>0
+         ipt = ipt + 1;
+         D(:,ipt)=squeeze(C(:,ix,iy));
+      end % if
+  end % for
+end %for
+
+% OUTPUT:
+D = squeeze(D(:,1:ipt));
+
+
Index: ../trunk-jpl/externalpackages/pcatool/install/mat2map.m
===================================================================
--- ../trunk-jpl/externalpackages/pcatool/install/mat2map.m     (revision 0)
+++ ../trunk-jpl/externalpackages/pcatool/install/mat2map.m     (revision 21589)
@@ -0,0 +1,41 @@
+% C = MAP2MAT(F,D) Reshaping matrix
+%
+% => Creation of a 'map' matrix C(TIME,LON,LAT) from D(TIME,PT)
+% following mask F(LON,LAT).
+% F is a matrix contenaing 1 where you would like to keep the point and
+%  0 elsewhere (see mat2map).
+% 
+% Rq: No check is done about the input.
+%
+% See also: map2mat
+%================================================================
+
+% March 2004
+% gmaze@univ-brest.fr
+
+function [C] = mat2map(F,D);
+
+% Get dimensions
+[nolon nolat] = size(F);
+[time npt] = size(D);
+
+% So output 'map' matrix has the form:
+C = zeros(time,nolon,nolat);
+
+% Variables
+nul = NaN.*ones(time,1);
+ipt = 0 ;
+
+
+% 'mapping' :
+  for iy=1:nolat
+      for ix=1:nolon
+         if F(ix,iy)>0
+            ipt = ipt + 1;
+            C(:,ix,iy) = D(:,ipt);
+         else
+            C(:,ix,iy) = nul;
+         end %if
+      end %for ix
+  end %for iy
+