function [A,B,cc] = ccakirby(X,Y)
%formulation in Johnson and Wichern but using SVD
%let X, Y be data matrix of data set I, II X = [x1;x2;...]

mx = mean(X);
MX = repmat(mx,size(X,1),1);
my=mean(Y);
MY = repmat(my,size(Y,1),1);

if 0
  S11 = (X-MX)'*(X-MX)/size(X,1);
  S22 = (Y-MY)'*(Y-MY)/size(Y,1);
  S12 = (X-MX)'*(Y-MY)/size(X,1);
else
  S11 = cov(X);
  S22 = cov(Y);
end  
S12 = (X-MX)'*(Y-MY)/(size(X,1)-1);

K = inv(real(sqrtm(S11)))*S12*inv(real(sqrtm(S22)));

%The matrix RHO contains the correlation coefs along its diagonal

[E,RHO,F] = svd(K);
E = real(E);
F = real(F);

A = inv(real(sqrtm(S11)))*E;
if nargout > 1
  B = inv(real(sqrtm(S22)))*F;
end

cc = diag(RHO);
