%==========================================================
function [Z, Z0] = rmsd_data()
%==========================================================
clear all
close all
view(25,0)
colors = ['k', 'r', 'g', 'b', 'c', 'm', 'k'];
symbol = ['^', 'o', '*', 's', 'x', '+', 'd'];
[a1,b1,c1]=sphere(10);format long
% Target
Z0 = [0.000   0.000   0.000;...
     1.520   0.000   0.000;...
     2.162   0.055   1.376;...
     1.779   1.262   2.217;...
     0.430   1.159   2.911;...
    -0.612   2.161   2.441;...
    -1.442   1.714   1.249;...
    -0.642   1.378   0.000]';
Z0 = set_mean(Z0); % target shifted to barycenter
% Model set; the one that is closest to the target
% is sought.
Z{1}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162  -0.379   1.324;...
       2.582   0.795   2.194;...
       2.426   2.162   1.547;...
       0.999   2.680   1.478;...
       0.330   2.540   0.120;...
      -0.642   1.378   0.000]';
Z{2}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162  -0.262   1.353;...
       1.695   0.653   2.473;...
       1.221   2.026   2.026;...
      -0.233   2.338   2.342;...
      -1.237   1.809   1.331;...
      -0.642   1.378   0.000]';
Z{3}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162  -0.010   1.378;...
       2.462   1.364   1.956;...
       1.419   1.894   2.926;...
       0.082   2.252   2.299;...
       0.109   2.437   0.790;...
      -0.642   1.378   0.000]';
Z{4}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162  -0.326  -1.338;...
       1.622   0.463  -2.520;...
       1.863   1.962  -2.454;...
       1.042   2.700  -1.410;...
      -0.386   2.206  -1.249;...
      -0.642   1.378   0.000]';
Z{5}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162   0.618  -1.231;...
       1.443   0.338  -2.540;...
       0.227   1.209  -2.805;...
       0.236   2.556  -2.101;...
       0.218   2.487  -0.582;...
      -0.642   1.378   0.000]';
Z{6}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162  -0.106   1.374;...
       2.029   1.134   2.242;...
       0.688   1.288   2.942;...
      -0.233   2.338   2.342;...
      -1.237   1.809   1.331;...
      -0.642   1.378   0.000]';
Z{7}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162   1.010   0.936;...
       2.031   0.689   2.416;...
       0.872   1.374   3.122;...
      -0.501   1.074   2.542;...
      -0.914   1.960   1.378;...
      -0.642   1.378   0.000]';
Z{8}= [0.000   0.000   0.000;...
       1.520   0.000   0.000;...
       2.162   1.375  -0.086;...
       2.927   1.805   1.156;...
       2.067   2.353   2.283;...
       0.589   2.489   1.955;...
       0.259   2.507   0.471;...
      -0.642   1.378   0.000]';
for i = 1:8
    Z{i} = set_mean(Z{i})
end

for i = 1:8
    draw_molecule(8, Z{i}, a1,b1,c1,1)
    draw_molecule(8, Z0,   a1,b1,c1,2)
    %figure
end

%==========================================================

function Xs = set_mean(X)
[M,N] = size(X);
Xm(1:M,1) = zeros(M,1); 
for i = 1:N
    Xm(:) = Xm(:) + X(:,i);
end
    Xm = (1/N)*Xm;
for i = 1:N
    Xs(:,i) = X(:,i) - Xm(:);
end
%end function set_mean

%==========================================================

function draw_molecule(N, atoms, a,b,c, k)
%[a,b,c]=sphere(10);
ats = zeros(3, N+1);
ats(1:3, 1:N) = atoms(1:3,1:N);
ats(1:3,N+1)  = atoms(1:3,1);
hold on
for i = 1:N
   surf(.05*a+ats(2,i),.05*b+ats(3,i),.05*c+ats(1,i));
end
if (k == 1)
   plot3(ats(2,:),ats(3,:),ats(1,:),'g-')
else
   plot3(ats(2,:),ats(3,:),ats(1,:),'r-') 
end
hold off
%==========================================================