## [x,y,z] = ellipsoidrand (a,b,c,n)
## Generate up to n random points uniformly distributed on the
## surface of an ellipsoid with axes a,b,c. The number of points
## returned may be less than n/2 due to the rejection technique
## used.
##
## Reference:
## Herman Rubin, Dept. of Statistics, Purdue University.
##
http://groups.google.com/group/sci.math/msg/24356105b2f3bc41?hl=en## Author: Paul Kienzle <
pkienzle@...>
## This program is public domain.
function [ ex, ey, ez ] = ellipsoidrand (a,b,c,n)
## Normalize the ellipsoid to have a smallest axis of 1.
ai = [a,b,c]./min([a,b,c]);
## Generate points with uniform density on a sphere surface
u = rand(n,1); theta = 2*pi*u;
v = rand(n,1); phi = asin(2*v-1);
x = cos(theta).*cos(phi); y=sin(theta).*cos(phi); z=sin(phi);
## Select points with probability inversely related to how
## far they are from the surface of the original sphere.
## The acceptance rate appears to be about 50% in the worst
## case of an elongated cigar shape.
t = rand(size(u));
idx = find((x/ai(1)).^2 + (y/ai(2)).^2 + (z/ai(3)).^2 >= t.^2);
acceptance_rate = 100*size(idx)/size(t)
## stretch the selected points onto the ellipse
ex = x(idx)*a; ey=y(idx)*b; ez = z(idx)*c;
endfunction
_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources
