Application Center - Maplesoft

Contact Maplesoft Request Quote

Submit your work

Application Center Applications Stochastic Model of Filling a Box with Spheres Preview

App Preview:

Stochastic Model of Filling a Box with Spheres

You can switch back to the summary page by clicking here.

Learn about Maple

Download Application

SpheresInBox.mws

Lower and Upper Bounds on Solution
- Lower Bound
  - Use Cubes Instead of Spheres for the Dimensions
- Upper Bound
  - Totally Fill the Volume
  - Tightly Packed Spheres
Source Code for the Stochastic Model of Filling the Box
Applying the Stochastic Model

Stochastic Model of Filling a Box with Spheres

Lee R. Partin

lpartin@chartertn.net

http://users.chartertn.net/lpartin

February 11, 2002

The task is to determine how many balls will fit into a given size box during the random filling of the box. This Maple document first calculates estimates to establish a lower and upper bound on the solution. Then, a stochastic model is programmed to randomly place balls into the box.

The task was based upon the fizz ball example presented to me by Wayne Pafko. ( http://www.pafko.com/wayne)

The box dimensions are set at 6-3/8" x 3" x 5-3/4". The spheres have a radius of 3/4".

Lower and Upper Bounds on Solution

Box dimensions:

`>`	restart;

`>`	Length:=6+3/8: Width:=3: Height:=5+3/4:

Sphere radius:

`>`	Radius:=3/4:

Include sphere plotting routine,

`>`	with(plottools):

Lower Bound

Use Cubes Instead of Spheres for the Dimensions

Fill the box with cubes that have side lengths equal to 2 times the sphere radius.

`>`	Side:=2Radius; CountX:=floor(Length/Side); CountY:=floor(Width/Side); CountZ:=floor(Height/Side); Count:=CountXCountY*CountZ; # total cubes that fit into box LowerBound:=Count;

Upper Bound

Totally Fill the Volume

The number of spheres in the box can not be more than the amount of spheres that have a total volume equal to the box volume.

`>`	BoxVolume:=LengthWidthHeight; SphereVolume:=4/3Pi(Radius)^3; UpperBound:=floor(BoxVolume/SphereVolume);

Tightly Packed Spheres

Spheres can pack tightly as shown in the following plot. Note that only one sphere is shown for the middle layer in the plot.

> c1:=sphere([0,0,0], Radius):
c2:=sphere([2*Radius,0,0], Radius):
c3:=sphere([2*Radius,2*Radius,0], Radius):
c4:=sphere([0,2*Radius,0], Radius):
DiagonalAngle:=arccos(sqrt(8)*Radius/(4*Radius));
Layer3Height:=4*Radius*sin(DiagonalAngle);
c5:=sphere([0,0,Layer3Height], Radius):
c6:=sphere([2*Radius,0,Layer3Height], Radius):
c7:=sphere([2*Radius,2*Radius,Layer3Height], Radius):
c8:=sphere([0,2*Radius,Layer3Height], Radius):
c9:=sphere([Radius,Radius,Layer3Height/2],Radius):
plots[display]([c1,c2,c3,c4,c5,c6,c7,c8,c9], scaling=constrained,
style=patch, axes=boxed);

[Maple Plot]

In this packing method, there are two complete spheres within a box having sides with lengths of Layer3Height x 2*Radius x 2*Radius. The volume taken up a one sphere is then:

`>`	*SphereSpace:=(Layer3Height2Radius2Radius)/2;*

If this packing method could fully use the box volume, the new upper bound estimate is:

`>`	UpperBound:=floor(BoxVolume/SphereSpace);

Source Code for the Stochastic Model of Filling the Box

The model randomly places equal size spheres into a box starting with the first layer followed by placing new spheres on top of existing ones.

Several geometic calculations are required for the model. Maple is used to derive the formulas. Then, a module called Box is programmed to bundle the calculations into a working package.

Geometric Model 1 (3 Spheres to Support New Sphere)

Given three spheres located in 3D, find the origin of a fourth sphere that rests on top of them. Return FAIL for the z coordinate if it is not feasible.

Let: r = radius of the spheres

(x1,y1,z1), (x2,y2,z2), (x3,y3,z3) = origin of the three given spheres

Then, the origin of the resting sphere is the point where the three given spheres would intersect if their diameters were doubled. The algebra is solved by Maple as follows:

`>`	eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2r)^2: eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2r)^2: eqn3:=(x-x3)^2+(y-y3)^2+(z-z3)^2=(2*r)^2: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):

The solutions are very detailed expressions. They are converted to Maple functions as follows:

`>`	assign(solution): fx:=unapply(x,x1,y1,z1,x2,y2,z2,x3,y3,z3,r): fy:=unapply(y,x1,y1,z1,x2,y2,z2,x3,y3,z3,r): fz:=unapply(z,x1,y1,z1,x2,y2,z2,x3,y3,z3,r):

The functions may return multiple values. Also check for an infeasible solution:

> fx2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res;
try
res:=[allvalues(fx(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
fy2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res;
try
res:=[allvalues(fy(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
fz2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res;
try
res:=[allvalues(fz(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:

Testing the results,

`>`	x1:=-0.5: y1:=-1: z1:=0: # sphere 1 x2:=1.75: y2:=1.5: z2:=0: # sphere 2 x3:=0: y3:=2.5: z3:=0: # sphere 3 r:=1: # radius

`>`	x4:=fx2(x1,y1,z1,x2,y2,z2,x3,y3,z3,r); y4:=fy2(x1,y1,z1,x2,y2,z2,x3,y3,z3,r); z4:=max(fz2(x1,y1,z1,x2,y2,z2,x3,y3,z3,r));

Viewing the Plot of Spheres

`>`	c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): c3:=sphere([x3,y3,z3], r): c4:=sphere([x4,y4,z4], r): plots[display]([c1,c2,c3,c4], scaling=constrained, style=patch, axes=boxed);

[Maple Plot]

Geometric Model 2 (2 Spheres and Fixed x to Support New Sphere)

Given two spheres located in 3D and fixed x origin location for a third sphere, find the origin of the third sphere that rests on top of the two spheres. Return FAIL for the z coordinate if it is not feasible.

Let: r = radius of the spheres

(x1,y1,z1), (x2,y2,z2) = origins of the two given spheres
x3 = known x location of the resting sphere

Then, the origin of the resting sphere is the point where the two given spheres would intersect if their diameters were doubled. The algebra is solved by Maple as follows:

`>`	unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2r)^2: eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2r)^2: eqn3:=x=x3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):

The solutions are very detailed expressions. They are converted to Maple functions as follows:

`>`	assign(solution): gy:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r): gz:=unapply(z,x1,y1,z1,x2,y2,z2,x3,r):

The functions may return multiple values. Also check for an infeasible solution:

> gy2:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;
try
res:=[allvalues(gy(x1,y1,z1,x2,y2,z2,x3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
gz2:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;
try
res:=[allvalues(gz(x1,y1,z1,x2,y2,z2,x3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:

Testing the results,

`>`	x1:=1.25: y1:=1: z1:=0: # sphere 1 x2:=2.25: y2:=3: z2:=0: # sphere 2 x3:=1: # sphere 3 x location r:=1: # radius

`>`	y3:=gy2(x1,y1,z1,x2,y2,z2,x3,r); z3:=max(gz2(x1,y1,z1,x2,y2,z2,x3,r));

Viewing the Plot of Spheres

`>`	c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): c3:=sphere([x3,y3,z3], r): plots[display]([c1,c2,c3], scaling=constrained, style=patch, axes=boxed);

[Maple Plot]

The above solutions fail with a zero divide when y1=y2. Therefore, a new set of equations is derived for this case.

> unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'):
eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:
eqn2:=(x-x2)^2+(y-y1)^2+(z-z2)^2=(2*r)^2:
eqn3:=x=x3:
solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):
assign(solution):
gy3:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r):
gz3:=unapply(z,x1,y1,z1,x2,y2,z2,x3,r):
gy4:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;
try
res:=[allvalues(gy3(x1,y1,z1,x2,y2,z2,x3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
gz4:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;
try
res:=[allvalues(gz3(x1,y1,z1,x2,y2,z2,x3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:

Another set is needed for the special case of y1=y2 and z1=z2.

> unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'):
eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:
eqn2:=(x-x2)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:
eqn3:=x=x3:
solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):
assign(solution):
gy5:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r):
gz5:=unapply(z,x1,y1,z1,x2,y2,z2,x3,r):
gy6:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;
try
res:=[allvalues(gy5(x1,y1,z1,x2,y2,z2,x3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
gz6:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;
try
res:=[allvalues(gz5(x1,y1,z1,x2,y2,z2,x3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:

Geometric Model 3 (2 Spheres and Fixed y to Support New Sphere)

Given two spheres located in 3D and fixed y origin location for a third sphere, find the origin of the third sphere that rests on top of the two spheres. Return FAIL for the z coordinate if it is not feasible.

Let: r = radius of the spheres

(x1,y1,z1), (x2,y2,z2) = origins of the two given spheres
y3 = known y location of the resting sphere

Then, the origin of the resting sphere is the point where the two given spheres would intersect if their diameters were doubled. The algebra is solved by Maple as follows:

`>`	unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2r)^2: eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2r)^2: eqn3:=y=y3: solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):

The solutions are very detailed expressions. They are converted to Maple functions as follows:

`>`	assign(solution): hx:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r): hz:=unapply(z,x1,y1,z1,x2,y2,z2,y3,r):

The functions may return multiple values. Also check for an infeasible solution:

> hx2:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;
try
res:=[allvalues(hx(x1,y1,z1,x2,y2,z2,y3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
hz2:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;
try
res:=[allvalues(hz(x1,y1,z1,x2,y2,z2,y3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:

Testing the results,

`>`	x1:=1.: y1:=1.3: z1:=0: # sphere 1 x2:=2.5: y2:=1.5: z2:=0: # sphere 2 y3:=1: # sphere 3 y location r:=1: # radius

`>`	x3:=hx2(x1,y1,z1,x2,y2,z2,y3,r); z3:=max(hz2(x1,y1,z1,x2,y2,z2,y3,r));

Viewing the Plot of Spheres

`>`	c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): c3:=sphere([x3,y3,z3], r): plots[display]([c1,c2,c3], scaling=constrained, style=patch, axes=boxed);

[Maple Plot]

The above solutions fail with a zero divide when x1=x2. Therefore, a new set of equations is derived for this case.

> unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'):
eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:
eqn2:=(x-x1)^2+(y-y2)^2+(z-z2)^2=(2*r)^2:
eqn3:=y=y3:
solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):
assign(solution):
hx3:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r):
hz3:=unapply(z,x1,y1,z1,x2,y2,z2,y3,r):
hx4:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;
try
res:=[allvalues(hx3(x1,y1,z1,x2,y2,z2,y3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
hz4:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;
try
res:=[allvalues(hz3(x1,y1,z1,x2,y2,z2,y3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:

Another set is needed for the special case of y1=y2 and z1=z2.

> unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'):
eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:
eqn2:=(x-x1)^2+(y-y2)^2+(z-z1)^2=(2*r)^2:
eqn3:=y=y3:
solution:=solve({eqn1,eqn2,eqn3},{x,y,z}):
assign(solution):
hx5:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r):
hz5:=unapply(z,x1,y1,z1,x2,y2,z2,y3,r):
hx6:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;
try
res:=[allvalues(hx5(x1,y1,z1,x2,y2,z2,y3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:
hz6:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;
try
res:=[allvalues(hz5(x1,y1,z1,x2,y2,z2,y3,r))]
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops(res)=2 and Im(res[1])<>0 then
RETURN(FAIL);
else RETURN(op(res)) fi;
end:

Geometric Model 4 (1 Sphere and Fixed (x,y) to Support New Sphere)

Given fixed (x,y) origin locations and the origin of a sphere, find the origin of the new sphere resting on it. Return FAIL for the z coordinate if it is not feasible.

Let: r = radius of the spheres

(x1,y1,z1) = origin of the given spheres
(x2,y2) = known (x,y) location of the resting sphere

The z-coordinate of the origin is found by Maple as follows:

`>`	unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'): eqn:=(x2-x1)^2+(y2-y1)^2+(z-z1)^2=(2*r)^2: solution:=solve(eqn,{z});

$solution := {z = z1+(-x2^2+2*x1*x2-x1^2-y1^2+2*y2*y1-y2^2+4*r^2)^(1/2)}, {z = z1-(-x2^2+2*x1*x2-x1^2-y1^2+2*y2*y1-y2^2+4*r^2)^(1/2)}$

A Maple function is programmed for the calculation.

`>`	iz:=proc (x1,y1,z1,x2,y2,r) local res; res:=-x2^2+2x1x2-x1^2-y2^2+2y1y2-y1^2+4*r^2; if res>0 then RETURN(z1+sqrt(res)) else RETURN(FAIL) fi; end:

Testing the results,

`>`	x1:=0: y1:=0: z1:=0: # sphere origin x2:=0.5: y2:=0.5: # fixed (x,y) of second sphere r:=1: # radius of spheres z2:=iz(x1,y1,z1,x2,y2,r);

Viewing the Plot of Spheres

`>`	c1:=sphere([x1,y1,z1], r): c2:=sphere([x2,y2,z2], r): plots[display]([c1,c2], scaling=constrained, style=patch, axes=boxed);

[Maple Plot]

Geometric Model 5 (Distance Between Two Spheres)

Given two spheres in 3D, find the distance between the origins.

Let (x1,y1,z1), (x2,y2,z2) = origins of the two given spheres .

`>`	Distance:=(x1,y1,z1,x2,y2,z2)->sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2);

Miscellaneous Routines

Random number generator

`>`	randomize(): randomX:=stats[random,uniform[0,1]]('generator'):

Box Module

Module Box() Programming

The Box module performs the calculations of randomly adding balls to the box, keeps track on the balls within the box and provides a means of plotting the results.

Definitions:

Box:-width = width of the box

Box:-depth = depth of the box

Box:-height = height of the box

Box:-radius = radius of the spheres to be placed into the box

Box:-SphereSet = set of data for the spheres within the box

Box:-Spheres = number of spheres within the box

Box:-fxbox, Box:-fybox, Box:-fzbox = routines that link with the calculations of section 'Geometric Model 1' for a sphere resting on three spheres
Box:-gybox, Box:-gzbox = routines that link with the calculations of section 'Geometric Model 2' for a sphere resting on two spheres and a wall (fixed x)

Box:-hxbox, Box:-hzbox = routines that link with the calculations of section 'Geometric Model 3' for a sphere resting on two spheres and a wall (fixed y)

Box:-makeSphere = routine for adding a sphere to the set. It tests to make sure that the sphere location is valid.

Box:-checkDistances = routine for checking the distances between a sphere and other spheres. It must be greater than 2*radius.

Box:-checkDistances2 = another distance checking routine

Box:-firstLayer = routine for randomly placing the first layer of spheres within the box.

Box:-moreSpheres = routine for randomly adding spheres to the box. It randomly picks three spheres from the set of spheres already within the box and tries to place a sphere on them.

> module Box()
description "Box for randomly placing spheres";
local i;
export width, depth, height, radius, makeSphere, SphereSet,
plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres,
fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzbox;
# initialize parameters
width:=6.+3./8.;
depth:=3.;
height:=5.+3./4.;
radius:=3./4.;
SphereSet:={};
Spheres:=0;
fxbox:=proc(T1::`name`,T2::`name`,T3::`name`)
description "place a sphere on 3 spheres - calculate x of origin";
local x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, res, randomList;
x1:=T1:-x;
x2:=T2:-x;
x3:=T3:-x;
y1:=T1:-y;
y2:=T2:-y;
y3:=T3:-y;
z1:=T1:-z;
z2:=T2:-z;
z3:=T3:-z;
try
res:=fx2(x1,y1,z1,x2,y2,z2,x3,y3,z3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops([res])=1 then
if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;
fi;
# for multiple real solutions, return both
RETURN([res]);
end proc:
fybox:=proc(T1::`name`,T2::`name`,T3::`name`)
description "place a sphere on 3 spheres - calculate y of origin";
local x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, res, randomList;
x1:=T1:-x;
x2:=T2:-x;
x3:=T3:-x;
y1:=T1:-y;
y2:=T2:-y;
y3:=T3:-y;
z1:=T1:-z;
z2:=T2:-z;
z3:=T3:-z;
try
res:=fy2(x1,y1,z1,x2,y2,z2,x3,y3,z3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops([res])=1 then
if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;
fi;
# for multiple real solutions, return both
RETURN([res]);
end proc:
fzbox:=proc(T1::`name`,T2::`name`,T3::`name`)
description "place a sphere on 3 spheres - calculate z of origin";
local x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, res, randomList;
x1:=T1:-x;
x2:=T2:-x;
x3:=T3:-x;
y1:=T1:-y;
y2:=T2:-y;
y3:=T3:-y;
z1:=T1:-z;
z2:=T2:-z;
z3:=T3:-z;
try
res:=fz2(x1,y1,z1,x2,y2,z2,x3,y3,z3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
if nops([res])=1 then
if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;
fi;
# for multiple solutions, return both
RETURN([res]);
end proc:
gybox:=proc(T1::`name`,T2::`name`,x3::`numeric`)
description "place a sphere on 2 spheres and a x-wall - calculate y of origin";
local x1, y1, z1, x2, y2, z2, y3, z3, res, randomList;
x1:=T1:-x;
x2:=T2:-x;
y1:=T1:-y;
y2:=T2:-y;
z1:=T1:-z;
z2:=T2:-z;
if y1=y2 and z1=z2 then
try
res:=gy6(x1,y1,z1,x2,y2,z2,x3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
elif y1=y2 then
try
res:=gy4(x1,y1,z1,x2,y2,z2,x3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
else
try
res:=gy2(x1,y1,z1,x2,y2,z2,x3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
fi;
if nops([res])=1 then
if res<min(y1,y2) or res>max(y1,y2) then res:=FAIL fi;
if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;
fi;
# for multiple real solutions, return both
res:=[res];
if res[1]<min(y1,y2) or res[1]>max(y1,y2) then res[1]:=FAIL fi;
if res[2]<min(y1,y2) or res[2]>max(y1,y2) then res[2]:=FAIL fi;
if res[1]=FAIL and res[2]=FAIL then RETURN(FAIL) fi;
RETURN(res);
end proc:
gzbox:=proc(T1::`name`,T2::`name`,x3::`numeric`)
description "place a sphere on 2 spheres and a x-wall - calculate z of origin";
local x1, y1, z1, x2, y2, z2, y3, z3, res, randomList;
x1:=T1:-x;
x2:=T2:-x;
y1:=T1:-y;
y2:=T2:-y;
z1:=T1:-z;
z2:=T2:-z;
if y1=y2 and z1=z2 then
try
res:=gz6(x1,y1,z1,x2,y2,z2,x3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
elif y1=y2 then
try
res:=gz4(x1,y1,z1,x2,y2,z2,x3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
else
try
res:=gz2(x1,y1,z1,x2,y2,z2,x3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
fi;
if nops([res])=1 then
if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;
fi;
# for multiple real solutions, return both
RETURN([res]);
end proc:
hxbox:=proc(T1::`name`,T2::`name`,y3::`numeric`)
description "place a sphere on 2 spheres and a y-wall - calculate x of origin";
local x1, y1, z1, x2, y2, z2, x3, z3, res, randomList;
x1:=T1:-x;
x2:=T2:-x;
y1:=T1:-y;
y2:=T2:-y;
z1:=T1:-z;
z2:=T2:-z;
if x1=yx and z1=z2 then
try
res:=hx6(x1,y1,z1,x2,y2,z2,y3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
elif x1=x2 then
try
res:=hx4(x1,y1,z1,x2,y2,z2,y3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
else
try
res:=hx2(x1,y1,z1,x2,y2,z2,y3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
fi;
if nops([res])=1 then
if res<min(x1,x2) or res>max(x1,x2) then res:=FAIL fi;
if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;
fi;
# for multiple real solutions, return both
res:=[res];
if res[1]<min(x1,x2) or res[1]>max(x1,x2) then res[1]:=FAIL fi;
if res[2]<min(x1,x2) or res[2]>max(x1,x2) then res[2]:=FAIL fi;
if res[1]=FAIL and res[2]=FAIL then RETURN(FAIL) fi;
RETURN(res);
end proc:
hzbox:=proc(T1::`name`,T2::`name`,y3::`numeric`)
description "place a sphere on 2 spheres and a y-wall - calculate z of origin";
local x1, y1, z1, x2, y2, z2, x3, z3, res, randomList;
x1:=T1:-x;
x2:=T2:-x;
y1:=T1:-y;
y2:=T2:-y;
z1:=T1:-z;
z2:=T2:-z;
if x1=x2 and z1=z2 then
try
res:=hz6(x1,y1,z1,x2,y2,z2,y3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
elif x1=x2 then
try
res:=hz4(x1,y1,z1,x2,y2,z2,y3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
else
try
res:=hz2(x1,y1,z1,x2,y2,z2,y3,radius);
catch: res:=FAIL;
end try;
if res=FAIL then RETURN(FAIL) fi;
fi;
if nops([res])=1 then
if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;
fi;
# for multiple real solutions, return both
RETURN([res]);
end proc:
makeSphere:=proc(f::`name`,xi::`complex`,yi::`complex`,
zi::`complex`)
description "add a sphere to the sphere set within the box";
local test;
# check for valid (x,y) values
#lprint(f,xi,yi,zi);
if xi<radius or xi>(width-radius) or
yi<radius or yi>(depth-radius) or
zi<radius or zi>(height-radius*0.5) or
Im(xi)<>0 or Im(yi)<>0 or Im(zi)<>0 then
RETURN('FAIL');
fi;
f:=module() export x, y, z;
x:=xi; y:=yi; z:=zi; end module;
Spheres:=Spheres+1;
SphereSet:=SphereSet union {f};
test:=checkDistances(Spheres);
if test='FAIL' then
SphereSet:=SphereSet minus {f};
unassign(S||Spheres);
Spheres:=Spheres - 1;
RETURN('FAIL');
fi;
Spheres;
end proc;
plot:=proc(SList::list)
description "create a 3D plot of the spheres in the box";
local i,l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12;
l1:=plottools[line]([0,0,0],[width,0,0],color=red);
l2:=plottools[line]([width,0,0],[width,depth,0],color=red);
l3:=plottools[line]([width,depth,0],[0,depth,0],color=red);
l4:=plott\ools[line]([0,depth,0],[0,0,0],color=red);
l5:=plottools[line]([0,0,height],[width,0,height],color=red);
l6:=plottools[line]([width,0,height],[width,depth,height],color=red);
l7:=plottools[line]([width,depth,height],[0,depth,height],color=red);
l8:=plottools[line]([0,depth,height],[0,0,height],color=red);
l9:=plottools[line]([0,0,0],[0,0,height],color=red);
l10:=plottools[line]([width,0,0],[width,0,height],color=red);
l11:=plottools[line]([width,depth,0],[width,depth,height],color=red);
l12:=plottools[line]([0,depth,0],[0,depth,height],color=red);
if nargs=0 then
for i from 1 to nops(SphereSet) do
c||i:=plottools[sphere]([SphereSet[i]:-x,SphereSet[i]:-y,
SphereSet[i]:-z], radius);
end do;
plots[display]([seq(c||i,i=1..nops(SphereSet)),
l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12],
title="Spheres in Box", scaling=constrained,
style=patch, axes=boxed);
else
for i from 1 to nops(SList) do
c||i:=plottools[sphere]([SList[i]:-x,SList[i]:-y,
SList[i]:-z], radius);
end do;
plots[display]([seq(c||i,i=1..nops(SList)),
l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12],
title="Spheres in Box", scaling=constrained,
style=patch, axes=boxed);
fi;
end proc;
firstLayer:=proc(times::posint)
description "randomly add spheres to the first layer";
local i, x, y, z, intList, randomList, irandom, Sphere1, movex, movey;
z:=radius; # origin height for first layer
for i from 1 to times do
x:=randomX();
x:=radius+x*(width-2.*radius);
y:=randomX();
y:=radius+y*(depth-2.*radius);
if modp(i,5)=0 then x:=radius fi;
if modp(i,10)=0 then x:=width-radius fi;
if modp(i,4)=0 then y:=radius fi;
if modp(i,8)=0 then y:=depth-radius fi;
makeSphere(S||(Spheres+1),x,y,z);
if modp(i,10)=0 then
# pick a random sphere from the box
intList:=[seq(i,i=1..Spheres)];
randomList:=combinat[randperm](intList);
irandom:=randomList[1];
Sphere1:=S||irandom;
# try moving it to a new location toward a lower x value.
movex:=randomX()*0.5; # find the random move
movey:=randomX()*0.25-0.125;
x:=Sphere1:-x; # save the old location
y:=Sphere1:-y;
# move the sphere
movex:=x-movex;
if movex<radius then movex:=radius fi;
movey:=y+movey;
if movey<radius then movey:=radius fi;
if movey>(depth-radius) then movey:=depth-radius fi;
Sphere1:-x:=movex;
Sphere1:-y:=movey;
#lprint(Sphere1,movex,movey,checkDistances(irandom),irandom,x,y);
if checkDistances(irandom)='FAIL' then # check that it is an ok move
Sphere1:-x:=x;
Sphere1:-y:=y;
#lprint("FAIL",Sphere1:-x, Sphere1:-y);
fi;
fi;
# try each corner
makeSphere(S||(Spheres+1),radius,radius,radius);
makeSphere(S||(Spheres+1),width-radius,radius,radius);
makeSphere(S||(Spheres+1),width-radius,depth-radius,radius);
makeSphere(S||(Spheres+1),radius,depth-radius,radius);
Spheres;
end do;
end proc;
checkDistances:=proc(SphereNumber::posint)
description "check sphere for being too close to another sphere";
local i, d;
for i from 1 to Spheres do
if i<>SphereNumber then
d:=Distance(S||i:-x,S||i:-y,S||i:-z,S||SphereNumber:-x,
S||SphereNumber:-y,S||SphereNumber:-z);
d:=evalf(d);
if (d+0.00001)<2.*radius then RETURN('FAIL') fi;
fi;
end do;
RETURN(0);
end proc;
checkDistances2:=proc(x::`numeric`,y::`numeric`,z::`numeric`)
description "check sphere for being too close to another sphere";
local i, d;
for i from 1 to Spheres do
d:=Distance(S||i:-x,S||i:-y,S||i:-z,S||SphereNumber:-x,
S||SphereNumber:-y,S||SphereNumber:-z);
d:=evalf(d);
if (d+0.00001)<2.*radius then RETURN('FAIL', S||i, d) fi;
end do;
RETURN(0);
end proc;
moreSpheres:=proc(maxInBox::posint, maxTries::posint, trySpheres::list)
description "add more spheres to the box by placing them on existing spheres and the walls";
local intList, randomList, Sphere1, Sphere2, Sphere3, i, imax, i1, i2, i3,
x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, res, oldSpheres, delta,
Pt1, Pt2, Pt3, Pt4, line1, line2, anglesum;
oldSpheres:=Spheres;
if nargs=3 then
imax:=1
else
imax:=maxTries
fi;
for i from 1 to imax do
if Spheres>oldSpheres then
lprint(Spheres);
oldSpheres:=Spheres;
fi;
if nargs=2 then
# randomly pick three spheres from the box
intList:=[seq(i,i=1..Spheres)];
randomList:=combinat[randperm](intList);
i1:=randomList[1];
i2:=randomList[2];
i3:=randomList[3];
Sphere1:=S||i1;
Sphere2:=S||i2;
Sphere3:=S||i3;
#lprint(Sphere1,Sphere2,Sphere3);
else
Sphere1:=trySpheres[1];
Sphere2:=trySpheres[2];
Sphere3:=trySpheres[3];
fi;
# will a sphere rest on the three random spheres?
z4:=fzbox(Sphere1,Sphere2,Sphere3);
if z4<>FAIL then
x4:=fxbox(Sphere1,Sphere2,Sphere3);
y4:=fybox(Sphere1,Sphere2,Sphere3);
fi;
if z4<>FAIL and x4<>FAIL and y4<>FAIL then
# try solution #1
# test that (x4,y4) of the solution fails within the triangle of the
# three supporting spheres. Otherwise, the sphere is not fully supported.
geometry[point](Pt1,[Sphere1:-x,Sphere1:-y]);
geometry[point](Pt2,[Sphere2:-x,Sphere2:-y]);
geometry[point](Pt3,[Sphere3:-x,Sphere3:-y]);
geometry[point](Pt4,[x4[1],y4[1]]);
try
anglesum:=geometry[FindAngle]
(geometry[line](line1,[Pt4,Pt1]),geometry[line](line2,[Pt4,Pt2]))
+geometry[FindAngle]
(geometry[line](line1,[Pt4,Pt2]),geometry[line](line2,[Pt4,Pt3]))
+geometry[FindAngle]
(geometry[line](line1,[Pt4,Pt1]),geometry[line](line2,[Pt4,Pt3]));
catch: anglesum:=0;
end try;
if abs(anglesum-2.*evalf(Pi))<0.001 then
try
res:=makeSphere(S||(Spheres+1),x4[1],y4[1],z4[1]);
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(3,Sphere1,Sphere2,Sphere3,S||Spheres)
fi;
fi;
# try solution #2
# test that (x4,y4) of the solution fails within the triangle of the
# three supporting spheres. Otherwise, the sphere is not fully supported.
geometry[point](Pt4,[x4[2],y4[2]]);
try
anglesum:=geometry[FindAngle]
(geometry[line](line1,[Pt4,Pt1]),geometry[line](line2,[Pt4,Pt2]))
+geometry[FindAngle]
(geometry[line](line1,[Pt4,Pt2]),geometry[line](line2,[Pt4,Pt3]))
+geometry[FindAngle]
(geometry[line](line1,[Pt4,Pt1]),geometry[line](line2,[Pt4,Pt3]));
catch: anglesum:=0;
end try;
if abs(anglesum-2.*evalf(Pi))<0.001 then
try
res:=makeSphere(S||(Spheres+1),x4[2],y4[2],z4[2]);
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(3,Sphere1,Sphere2,Sphere3,S||Spheres)
fi;
fi;
fi;
# try lots of positions every 3th time
if nargs=3 or modp(i,3)=0 then
# will a sphere rest on the first two random spheres and a wall?
x4:=radius;
z4:=gzbox(Sphere1,Sphere2,x4);
if z4<>FAIL then
y4:=gybox(Sphere1,Sphere2,x4);
fi;
if z4<>FAIL and y4<>FAIL then
try
res:=FAIL;
if y4[1]<>FAIL and z4<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
try
res:=FAIL;
if y4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[2],z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
fi;
x4:=width-radius;
z4:=gzbox(Sphere1,Sphere2,x4);
if z4<>FAIL then
y4:=gybox(Sphere1,Sphere2,x4);
fi;
if z4<>FAIL and y4<>FAIL then
try
res:=FAIL;
if y4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
try
res:=FAIL;
if y4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[2],z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
fi;
y4:=radius;
z4:=hzbox(Sphere1,Sphere2,y4);
if z4<>FAIL then
x4:=hxbox(Sphere1,Sphere2,y4);
fi;
if z4<>FAIL and x4<>FAIL then
try
res:=FAIL;
if x4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
fi;
if z4<>FAIL and x4<>FAIL then
try
res:=FAIL;
if x4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[2],y4,z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
fi;
y4:=depth-radius;
z4:=hzbox(Sphere1,Sphere2,y4);
if z4<>FAIL then
x4:=hxbox(Sphere1,Sphere2,y4);
fi;
if z4<>FAIL and x4<>FAIL then
try
res:=FAIL;
if x4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
try
res:=FAIL;
if x4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[2],y4,z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere2,S||Spheres)
fi;
fi;
# will a sphere rest on the second two random spheres and a wall?
x4:=radius;
z4:=gzbox(Sphere2,Sphere3,x4);
if z4<>FAIL then
y4:=gybox(Sphere2,Sphere3,x4);
fi;
if z4<>FAIL and y4<>FAIL then
try
res:=FAIL;
if y4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if y4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[2],z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
fi;
x4:=width-radius;
z4:=gzbox(Sphere2,Sphere3,x4);
if z4<>FAIL then
y4:=gybox(Sphere2,Sphere3,x4);
fi;
if z4<>FAIL and y4<>FAIL then
try
res:=FAIL;
if y4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if y4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[2],z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
fi;
y4:=radius;
z4:=hzbox(Sphere2,Sphere3,y4);
if z4<>FAIL then
x4:=hxbox(Sphere2,Sphere3,y4);
fi;
if z4<>FAIL and x4<>FAIL then
try
res:=FAIL;
if x4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if x4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[2],y4,z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
fi;
y4:=depth-radius;
z4:=hzbox(Sphere2,Sphere3,y4);
if z4<>FAIL then
x4:=hxbox(Sphere2,Sphere3,y4);
fi;
if z4<>FAIL and x4<>FAIL then
try
res:=FAIL;
if x4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if x4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere2,Sphere3,S||Spheres)
fi;
fi;
# will a sphere rest on the first and third random spheres and a wall?
x4:=radius;
z4:=gzbox(Sphere1,Sphere3,x4);
if z4<>FAIL then
y4:=gybox(Sphere1,Sphere3,x4);
fi;
if z4<>FAIL and y4<>FAIL then
try
res:=FAIL;
if y4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if y4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[2],z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
fi;
x4:=width-radius;
z4:=gzbox(Sphere1,Sphere3,x4);
if z4<>FAIL then
y4:=gybox(Sphere1,Sphere3,x4);
fi;
if z4<>FAIL and y4<>FAIL then
try
res:=FAIL;
if y4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if y4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
fi;
y4:=radius;
z4:=hzbox(Sphere1,Sphere3,y4);
if z4<>FAIL then
x4:=hxbox(Sphere1,Sphere3,y4);
fi;
if z4<>FAIL and x4<>FAIL then
try
res:=FAIL;
if x4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if x4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[2],y4,z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
fi;
y4:=depth-radius;
z4:=hzbox(Sphere1,Sphere3,y4);
if z4<>FAIL then
x4:=hxbox(Sphere1,Sphere3,y4);
fi;
if z4<>FAIL and x4<>FAIL then
try
res:=FAIL;
if x4[1]<>FAIL and z4[1]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
try
res:=FAIL;
if x4[2]<>FAIL and z4[2]<>FAIL then
res:=makeSphere(S||(Spheres+1),x4[2],y4,z4[2]) fi;
catch: res:=FAIL;
end try;
if res<>FAIL then
lprint(2,Sphere1,Sphere3,S||Spheres)
fi;
fi;
fi;
# try corner positions every 10th time
if nargs=3 or modp(i,10)=0 then
# check on resting a sphere on the first random sphere and a corner
x1:=Sphere1:-x;
y1:=Sphere1:-y;
z1:=Sphere1:-z;
x4:=radius;
y4:=radius;
z4:=iz(x1,y1,z1,x4,y4,radius);
if type(z4,numeric)=true and (z4-z1)<1.999*radius then
res:=makeSphere(S||(Spheres+1),x4,y4,z4);
if type(res,numeric)=true then
lprint("corner1",Sphere1,S||Spheres)
fi;
fi;
x4:=radius;
y4:=depth-radius;
z4:=iz(x1,y1,z1,x4,y4,radius);
if type(z4,numeric)=true and (z4-z1)<1.999*radius then
res:=makeSphere(S||(Spheres+1),x4,y4,z4);
if type(res,numeric)=true then
lprint("corner2",Sphere1,S||Spheres)
fi;
fi;
x4:=width-radius;
y4:=radius;
z4:=iz(x1,y1,z1,x4,y4,radius);
if type(z4,numeric)=true and (z4-z1)<1.999*radius then
res:=makeSphere(S||(Spheres+1),x4,y4,z4);
if type(res,numeric)=true then
lprint("corner3",Sphere1,S||Spheres)
fi;
fi;
x4:=width-radius;
y4:=depth-radius;
z4:=iz(x1,y1,z1,x4,y4,radius);
if type(z4,numeric)=true and (z4-z1)<1.999*radius then
res:=makeSphere(S||(Spheres+1),x4,y4,z4);
if type(res,numeric)=true then
lprint("corner4",Sphere1,S||Spheres)
fi;
fi;
if Spheres>=maxInBox then RETURN(Spheres) fi;
fi;
end do;
RETURN(Spheres);
end proc;
end module;

Applying the Stochastic Model

The Box module and supporting functions are used to model the random filling of the box.

First Layer

The spheres of the first layer are placed into the box by random selection of (x,y) coordinates. It takes many tries to find locations that are not already controlled by existing spheres. Every 10th try, a random sphere is randomly moved toward the x-axis to make space for more spheres.

`>`	Box:-firstLayer(5000);

Fill the Box

The box is filled by placing more spheres on top of the existing spheres. Three spheres are randomly selected from the box and tested to see if sphere can be supported on top of them. On every 3rd time, tests are performed to see if a sphere could rest upon two of the random spheres and a wall. On every 10th time, test are performed to see if a sphere could rest on the random sphere and one of the corners.

`>`	Box:-moreSpheres(40,2000);

2, S3, S6, S8

2, S5, S6, S9

2, S4, S7, S10

2, S10, S3, S11

2, S2, S1, S12

"corner3", S10, S13

"corner4", S10, S14

2, S14, S8, S15

2, S9, S12, S16

2, S14, S15, S17

2, S17, S13, S18

2, S12, S11, S19

2, S19, S16, S20

"corner2", S16, S21

2, S9, S15, S22

2, S17, S22, S23

"corner1", S20, S24

2, S19, S11, S25

2, S25, S18, S26

2, S23, S21, S27

2, S26, S20, S28

Plot the Spheres Within the Box

Plotting all spheres within the box,

`>`	Box:-plot();

[Maple Plot]

>

Maple

Maple Add-Ons

MapleSim

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

App Preview:

Stochastic Model of Filling a Box with Spheres

Maple

Powerful math software that is easy to use

Maple Add-Ons

MapleSim

Advanced System Level Modeling

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

App Preview:

Stochastic Model of Filling a Box with Spheres