
%Compute distance matrix:

%Initializing:
clear;
Npts=10000;

%Generating a random orbit:
rng('default'); % re-initialize rand so that results are reproducible
x=rand(Npts,1);
y=rand(Npts,1);

%Using internal built-in function:
tic
XY=[x,y]; % collate x's and y's
D1 = pdist2(XY,XY);
T1=toc;
fprintf('It took:\n')
fprintf('%f sec for built-in\n',T1)

%Using loop within loop:
tic
D2=zeros(Npts,Npts);
for jj=1:Npts
 for ii=1:Npts
  D2(ii,jj)=(x(ii)-x(jj))^2 + (y(ii)-y(jj))^2;
 end
end
D2=sqrt(D2);
T2=toc;
fprintf('%f sec for loop-in-loop\n',T2)

%Using meshgrid:
tic
[X1,X2]=meshgrid(x,x);
[Y1,Y2]=meshgrid(y,y);
D3=sqrt((X1-X2).^2+(Y1-Y2).^2);
T3=toc;
fprintf('%f sec for meshgrid\n',T3)

%Using Kronecker product;
tic
X1=kron(x,ones(1,Npts));
Y1=kron(y,ones(1,Npts));
D4=sqrt((X1-X1.').^2+(Y1-Y1.').^2);
T4=toc;
fprintf('%f sec for Kronecker product\n',T4)

%difference between them: (just to check that they are equal)
diff_D1_D2=max(max(abs(D1-D2)));
diff_D1_D3=max(max(abs(D1-D3)));
diff_D1_D4=max(max(abs(D1-D4)));
if(diff_D1_D2>10*eps|diff_D1_D3>10*eps|diff_D1_D4>10*eps) error('Distance matrices not equal!');end

%timings:
fprintf('Built-in is %1.1f times faster than loop-in-loop!!!\n',round(10*T2/T1)/10)
fprintf('Built-in is %1.1f times faster than meshgrid!!!\n',round(10*T3/T1)/10)
fprintf('Built-in is %1.1f times faster than Kronecker!!!\n',round(10*T4/T1)/10)

