"Pstrang Rzekle" wrote in message <itd7aa$hl2$1@newscl01ah.mathworks.com>...
> Can anyone come up with code equivalent to the following for loop?
>
> m1 should be a 100x3 matrix.
>
> for n = 1:100
> d(n,:) = (acos((m1(n,1).*m1(:,1)) ...
> + (m1(n,2).*m1(:,2)) ...
> + (m1(n,3).*m1(:,3))))';
> end
>
> The result d should be a 100x100 matrix. Each row of d is the distance to all the other indices of m1. Obviously the distance to itself is 0, so I would expect diagonal marching zeros down the matrix.
>
> Big thanks!
         
There is a numerical weakness in this method of computing angles. When the unit vectors are very nearly in the same or opposite directions, there is a loss of accuracy. Just look at a plot of 'acos' between 1 and +1 to see why.
For that reason I tend to recommend a more robust formula using the 'atan2' function, even though it is more expensive computationwise:
a = atan2(norm(cross(v1,v2)),dot(v1,v2));
It has the added advantage that neither v1 nor v2 need be normalized.
To apply this to your case you could either use for loops or do the following:
n = size(m1,1);
[I,J] = ndgrid(1:n);
d = reshape(sqrt(sum(cross(m1(I,:),m1(J,:),2).^2,2)),n,n);
d = atan2(d,m1*m1.');
Note that you cannot expect the diagonal to have exact zeros since there will inevitably be rounding errors in the calculations.
Roger Stafford
