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: