function [maxR,pcind,RT,T,P,T_v,maxR_v] = rotation(R_c,grps,maxrank,R_v)
%%  obtaining the resolutions of all the PC subspaces and finding the PC
%%  subspace with maximal resolution
% Input
%      R_c: the scores of the training set after PCA
%      grps: the structure of each class
%      marank: the number of PCs used in PCA
%      R_v:  the scores of the test set after projection
% Output
%       maxR: the training set with maximal Rs after the rotation and projection 
%       pcind: the PC subspace with maximal Rs
%       RT1: the resolutions of all the PCˇˇsubspaces
%       T: The scores of the training set after the rotation and  projection
%       P: the index of all PC subspaces
%       T_v: The scores of the test set after the rotation and  projection
%       maxR_v: the test set with maximal Rs after the rotation and projection

grpnumber = size(grps,2);

for i = 1:grpnumber
    mm(i) = mean(R_c((grps(i).index),1));
end

[tmp, ind] = sort(mm); %Sort in ascending order
RT = []; T = []; P = []; T_v = [];

sumR = zeros(size(R_c(:,1))); sumR_v = zeros(size(R_v(:,1))); maxRs = 0; 
for d1 = 1:maxrank-1;
 for d2 = d1+1:maxrank
    for i = 1:grpnumber-1  
        for j = i+1:grpnumber  
            x1 = [R_c(grps(ind(i)).index,d1) R_c(grps(ind(i)).index,d2)];
            x2 = [R_c(grps(ind(j)).index,d1) R_c(grps(ind(j)).index,d2)];
            X = [R_c(:,d1) R_c(:,d2)];
            X_v = [R_v(:,d1) R_v(:,d2)];%
            [XR1, XR2, XR1_v, XR2_v] = rotpcs(x1, x2, X, X_v); %
            R(:,d1) = (XR1-min(XR1))/(max(XR1)-min(XR1));  
            R(:,d2) = (XR2-min(XR2))/(max(XR2)-min(XR2));  
            RR_v(:,d1) = (XR1_v-min(XR1_v))/(max(XR1_v)-min(XR1_v)); % 
            RR_v(:,d2) = (XR2_v-min(XR2_v))/(max(XR2_v)-min(XR2_v)); %
            %disprotpcsallpoints(R, grps,2,8);
            T = [T R(:,d1)];P = [P;[d1 d2]];T_v = [T_v RR_v(:,d1)]; %
            tR = (R(:,d1)-min(R(:,d1)))./(max(R(:,d1))-min(R(:,d1)));
            DR = mean(tR(grps(i).index))-mean(tR(grps(j).index)); 
            STDR = std(tR(grps(i).index))+std(tR(grps(j).index)); 
            Rs = abs(DR/STDR); 
            [rt1 meanR11 meanR22 stdR11 stdR22] = dispsumRs(tR,grps);
            RT = [RT;rt1]; 
            if Rs>maxRs, maxRs = Rs; maxR = R(:,d1); pcind = [d1,d2]; end 
        end         
    end
 end
end
m1 = find(RT==maxRs);
maxR_v = (T_v(:,m1)-min(T_v(:,m1)))./(max(T_v(:,m1))-min(T_v(:,m1)));


