function [p,F,df1,df2]=pc2(Gr1,Gr2)
% [p,F,df1,df2]=pc2(Gr1,Gr2)    or    [p,F,df1,df2]=pc2(Data,n1)
% In the first form, Gr1 and Gr2 are data from two independent groups of cases and have identical dimensions, except possibly for the number of
% levels of the last dimension, which correspond to the number of cases in each group
% In the second form, all the data are in Data (with cases as last index) and n1 indicates that the first n1 cases from Data belong to one group, 
% the remainder belonging to the other group
% if Data or Gr1 and Gr2 have more than two dimensions, they are made two dimensional by merging all but their last
% dimension
% p is the probability of the F ratio produced by Hotelling T^2 for the difference between the means of independent samples on two variables
% The two variables are the weights of the the individual cases on the first two principal components of all cases taken together
% (the methods saves computations by extracting the two PCs in the case domains, which differs from ordinary
% factor scores only by a multiplicative constant)
% The theory is that independent replications within two groups contain at most one systematic
% pattern for each group, therefore at most two systematic patterns in two groups.  If the groups do
% not differ, their means on the two most prominent patterns should not differ.
% The F test has (df1,df2) degrees of freedom with df1 normally 2 if k>1 (or if j*k>1 for 3 dimensional data),
% but df1=1 if this function is used for a simple independent group t- test with a single variable
% Reference:
% Achim, A.  Statistical detection of between groups differences in event-related potentials.  Clinical Neurophysiology, 2001, 112: 1023-1034.
a=size(Gr1);
n1=a(end);
k=length(a);
if length(Gr2)>1,
    b=size(Gr2);
    if length(b)~=k | any(a(1:k-1)~=b(1:k-1)), error('The two matrices passed to PC2 are not compatible.'); end 
    if k>2,  %if 3 dimensions, transform to 2 dimensions
        D=[reshape(Gr1,prod(a(1:k-1)),n1) reshape(Gr2,prod(b(1:k-1)),b(k))];
    else
        D=[Gr1 Gr2];    
    end
else
    if length(a)>2, D=reshape(Gr1,prod(a(1:k-1)),n1); else D=Gr1; end
    n1=Gr2;
end
n=size(D,1);
if n>2,
%	opt.disp=0;    % pour éviter des rapports de eigs sur l'écran
%   [V,d,flg]=eigs(C'*C,2,'LM',opt);
	V=extrais_n_cp(D'*D,2);  % plutôt que par eigs qui a une composante aléatoire
    [p,t2,F,df1,df2]=hotelling2gr(V',n1);
else
    [p,t2,F,df1,df2]=hotelling2gr(D,n1);
end