Application Center - Maplesoft

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

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:=2*Radius;
CountX:=floor(Length/Side);
CountY:=floor(Width/Side);
CountZ:=floor(Height/Side);
Count:=CountX*CountY*CountZ; # total cubes that fit into box
LowerBound:=Count;

Side := 3/2

CountX := 4

CountY := 2

CountZ := 3

Count := 24

LowerBound := 24

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:=Length*Width*Height;
SphereVolume:=4/3*Pi*(Radius)^3;
UpperBound:=floor(BoxVolume/SphereVolume);

BoxVolume := 3519/32

SphereVolume := 9/16*Pi

UpperBound := 62

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);

DiagonalAngle := 1/4*Pi

Layer3Height := 3/2*2^(1/2)

[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:=(Layer3Height*2*Radius*2*Radius)/2;

SphereSpace := 27/16*2^(1/2)

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

>    UpperBound:=floor(BoxVolume/SphereSpace);

UpperBound := 46

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=(2*r)^2:
eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^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));

x4 := .1297169812

y4 := .6957547170

z4 := .8531544192

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=(2*r)^2:
eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^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));

y3 := 2.375000000

z3 := 1.430690393

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=(2*r)^2:
eqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^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));

x3 := 1.803333334

z3 := 1.806835785

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+2*x1*x2-x1^2-y2^2+2*y1*y2-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);

z2 := 1.870828693

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);

Distance := proc (x1, y1, z1, x2, y2, z2) options operator, arrow; sqrt((-x2+x1)^2+(y1-y2)^2+(-z2+z1)^2) end proc

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;

module Box () local i; export width, depth, height, radius, makeSphere, SphereSet, plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres, fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzb...
module Box () local i; export width, depth, height, radius, makeSphere, SphereSet, plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres, fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzb...
module Box () local i; export width, depth, height, radius, makeSphere, SphereSet, plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres, fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzb...
module Box () local i; export width, depth, height, radius, makeSphere, SphereSet, plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres, fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzb...
module Box () local i; export width, depth, height, radius, makeSphere, SphereSet, plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres, fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzb...
module Box () local i; export width, depth, height, radius, makeSphere, SphereSet, plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres, fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzb...
module Box () local i; export width, depth, height, radius, makeSphere, SphereSet, plot, firstLayer, Spheres, checkDistances, checkDistances2, moreSpheres, fxbox, fybox, fzbox, gybox, gzbox, hxbox, hzb...

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);

7

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

9

2, S4, S7, S10

10

2, S10, S3, S11

11

2, S2, S1, S12

12

"corner3", S10, S13

"corner4", S10, S14

14

2, S14, S8, S15

15

2, S9, S12, S16

16

2, S14, S15, S17

17

2, S17, S13, S18

18

2, S12, S11, S19

19

2, S19, S16, S20

20

"corner2", S16, S21

21

2, S9, S15, S22

22

2, S17, S22, S23

23

"corner1", S20, S24

24

2, S19, S11, S25

25

2, S25, S18, S26

26

2, S23, S21, S27

27

2, S26, S20, S28

28

28

Plot the Spheres Within the Box

Plotting all spheres within the box,

>    Box:-plot();

[Maple Plot]

>