gnuplot - plot random points inside an ellipsoid

Question

I want to create and plot a set of random points within the boundaries of an ellipsoid using gnuplot. Is it possible to do this directly in gnuplot, or would I need to generate my random data points in an external program?

Ultimately I would like to produce an ellipsoid figure similar to this ellipse figure.

There are some examples on the use of rand and cylindrical/spherical coordinates, but I'm not sure how to generate the random points inside the boundaries of an ellipsoid.

Answer 1

Based on @Bernhard's answer, here is how you can do that with gnuplot only. To reuse a random number, you can put the two rand calls and the assignments to variables in the first argument of the using statement, separated by commas. The using statements are evaluated from left to right, so you can access those variables in all following using parameters.

To demonstrate this, see the following examples:

set samples 1000
plot '+' using (x=rand(0), y=rand(0), x):(y)

Applying this to the ellipsoid, this gives the script:

a=3
b=2
phi=30*pi/180

max(x,y) = (x > y ? x : y)
set xrange[-max(a,b):max(a,b)]
set yrange[-max(a,b):max(a,b)]

set offset 0.1,0.1,0.1,0.1
set samples 2000

ex(x, y) = a*(2*x-1)
ey(x, y) = b*(sqrt(1-((2*x-1))**2))*(2*y-1)

unset key    
plot '+' using (x=rand(0), y=rand(0), ex(x,y)*cos(phi)-ey(x,y)*sin(phi)):\
               (ey(x,y)*cos(phi)+ex(x,y)*sin(phi)) pt 7 ps 0.5

with the result:

This, however, leads to an apparently unequal distribution of the points (see the agglomarations at the ellipse's ends, see @andyras's comment). To avoid this, here is how you can filter equally distributed random points to be inside the ellipsoid:

a=3
b=2
set angles degree
phi=30

max(x,y) = (x > y ? x : y)
set xrange[-max(a,b):max(a,b)]
set yrange[-max(a,b):max(a,b)]

set offset 0.1,0.1,0.1,0.1
set samples 2000
set size ratio 1

check(x, y) = (((x/a)**2 + (y/b)**2) <= 1)
unset key
plot '+' using (x=2*a*(rand(0)-0.5), y=2*b*(rand(0)-0.5), \
     check(x,y) ? x*cos(phi)-y*sin(phi) : 1/0):\
     (x*sin(phi)+y*cos(phi)) pt 7 ps 0.5

This gives the much better result:

Extending this to three dimensions:

a=3
b=1
c=1
set angles degree
phi=30

mx(x,y) = (x > y ? x : y)
max(x,y,z) = mx(mx(x,y), mx(x,z))
set xrange[-max(a,b,c):max(a,b,c)]
set yrange[-max(a,b,c):max(a,b,c)]
set zrange[-max(a,b,c):max(a,b,c)]

set offset 0.1,0.1,0.1,0.1
set samples 2000
set size ratio 1

set ticslevel 0
set view 60, 330

check(x, y, z) = (((x/a)**2 + (y/b)**2 + (z/c)**2) <= 1)
unset key
splot '+' using (x = 2*a*(rand(0)-0.5), \
                 y = 2*b*(rand(0)-0.5), \
                 z=2*c*(rand(0)-0.5), \
                 check(x,y,z) ? x*cos(phi)-y*sin(phi) : 1/0):\
                 (x*sin(phi)+y*cos(phi)):(z) pt 7 ps 0.5

with the result: