{VERSION 5 0 "IBM INTEL NT" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "Hyperlink" -1 17 "" 0 1 0 128 128 1 2 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1" -1 3 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 2" -1 4 1 {CSTYLE "" -1 -1 "Times " 1 14 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 } {PSTYLE "Heading 3" -1 5 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 1 1 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Plot" -1 13 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Title" -1 18 1 {CSTYLE "" -1 -1 "Times " 1 18 0 0 0 1 2 1 1 2 2 2 1 1 1 1 }3 1 0 0 12 12 1 0 1 0 2 2 19 1 } {PSTYLE "Author" -1 19 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 8 8 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Output" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 3 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 46 "Stochastic Model of Filli ng a Box with Spheres" }}{PARA 19 "" 0 "" {TEXT -1 13 "Lee R. Partin" }}{PARA 0 "" 0 "" {TEXT -1 21 "lpartin@chartertn.net" }}{PARA 0 "" 0 " " {URLLINK 17 "http://users.chartertn.net/lpartin" 4 "" "" }}{PARA 0 " " 0 "" {TEXT -1 17 "February 11, 2002" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 288 "The task is to determine how many ba lls will fit into a given size box during the random filling of the bo x. This Maple document first calculates estimates to establish a lowe r and upper bound on the solution. Then, a stochastic model is progra mmed to randomly place balls into the box." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 80 "The task was based upon the fiz z ball example presented to me by Wayne Pafko. (" }{URLLINK 17 "http: //www.pafko.com/wayne" 4 "" "" }{TEXT -1 1 ")" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 87 "The box dimensions are se t at 6-3/8\" x 3\" x 5-3/4\". The spheres have a radius of 3/4\"." }} }{SECT 0 {PARA 3 "" 0 "" {TEXT -1 34 "Lower and Upper Bounds on Soluti on" }}{PARA 0 "" 0 "" {TEXT -1 15 "Box dimensions:" }}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "Length:=6+3/8:\nWidth:=3:\nHeight:=5+3/4: " }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 14 "Sphere radius:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Radius:=3/4:" }{TEXT -1 0 " " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 32 "Include sphere plotting routi ne," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "with(plottools):" }} }{SECT 0 {PARA 4 "" 0 "" {TEXT -1 11 "Lower Bound" }}{SECT 0 {PARA 5 " " 0 "" {TEXT -1 47 "Use Cubes Instead of Spheres for the Dimensions" } }{PARA 0 "" 0 "" {TEXT -1 82 "Fill the box with cubes that have side l engths equal to 2 times the sphere radius." }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 178 "Side:=2*Radius;\nCountX:=floor(Length/Side);\nCoun tY:=floor(Width/Side);\nCountZ:=floor(Height/Side);\nCount:=CountX*Cou ntY*CountZ; # total cubes that fit into box\nLowerBound:=Count;" }}}}} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 11 "Upper Bound" }}{SECT 0 {PARA 5 " " 0 "" {TEXT -1 23 "Totally Fill the Volume" }}{PARA 0 "" 0 "" {TEXT -1 125 "The number of spheres in the box can not be more than the amou nt of spheres that have a total volume equal to the box volume." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 107 "BoxVolume:=Length*Width*Hei ght;\nSphereVolume:=4/3*Pi*(Radius)^3;\nUpperBound:=floor(BoxVolume/Sp hereVolume);" }}}}{SECT 0 {PARA 5 "" 0 "" {TEXT -1 22 "Tightly Packed \+ Spheres" }}{PARA 0 "" 0 "" {TEXT -1 126 "Spheres can pack tightly as s hown in the following plot. Note that only one sphere is shown for th e middle layer in the plot." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 582 "c1:=sphere([0,0,0], Radius):\nc2:=sphere([2*Radius,0,0], Radius): \nc3:=sphere([2*Radius,2*Radius,0], Radius):\nc4:=sphere([0,2*Radius,0 ], Radius):\nDiagonalAngle:=arccos(sqrt(8)*Radius/(4*Radius));\nLayer3 Height:=4*Radius*sin(DiagonalAngle);\nc5:=sphere([0,0,Layer3Height], R adius):\nc6:=sphere([2*Radius,0,Layer3Height], Radius):\nc7:=sphere([2 *Radius,2*Radius,Layer3Height], Radius):\nc8:=sphere([0,2*Radius,Layer 3Height], Radius):\nc9:=sphere([Radius,Radius,Layer3Height/2],Radius): \nplots[display]([c1,c2,c3,c4,c5,c6,c7,c8,c9], scaling=constrained,\n \+ style=patch, axes=boxed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 175 "In this packing method, there are two complete spheres w ithin a box having sides with lengths of Layer3Height x 2*Radius x 2*R adius. The volume taken up a one sphere is then:" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 48 "SphereSpace:=(Layer3Height*2*Radius*2*Radius )/2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 87 "If this packing method co uld fully use the box volume, the new upper bound estimate is:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "UpperBound:=floor(BoxVolume/ SphereSpace);" }}}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 55 "Source Code \+ for the Stochastic Model of Filling the Box" }}{PARA 0 "" 0 "" {TEXT -1 142 "The model randomly places equal size spheres into a box starti ng with the first layer followed by placing new spheres on top of exis ting ones." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 188 "Several geometic calculations are required for the model. Map le is used to derive the formulas. Then, a module called Box is progr ammed to bundle the calculations into a working package." }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 51 "Geometric Model 1 (3 Spheres to Support N ew Sphere)" }}{PARA 0 "" 0 "" {TEXT -1 153 "Given three spheres locate d 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." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 32 "Let: r = radius o f the spheres " }}{PARA 0 "" 0 "" {TEXT -1 77 " (x1,y1,z1), (x2, y2,z2), (x3,y3,z3) = origin of the three given spheres" }}{PARA 0 "" 0 "" {TEXT -1 174 "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:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 168 "eqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2 :\neqn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^2:\neqn3:=(x-x3)^2+(y-y3)^2+ (z-z3)^2=(2*r)^2:\nsolution:=solve(\{eqn1,eqn2,eqn3\},\{x,y,z\}):" }}} {PARA 0 "" 0 "" {TEXT -1 95 "The solutions are very detailed expressio ns. They are converted to Maple functions as follows:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 152 "assign(solution):\nfx:=unapply(x,x 1,y1,z1,x2,y2,z2,x3,y3,z3,r):\nfy:=unapply(y,x1,y1,z1,x2,y2,z2,x3,y3,z 3,r):\nfz:=unapply(z,x1,y1,z1,x2,y2,z2,x3,y3,z3,r):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "The functions may return multiple values. Also check for an infeasible solution:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 811 "fx2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res; \n try\n res:=[allvalues(fx(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]\n catch : res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi;\n if no ps(res)=2 and Im(res[1])<>0 then \n RETURN(FAIL);\n else RETURN(o p(res)) fi; \nend:\nfy2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local res ;\n try\n res:=[allvalues(fy(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]\n catc h: res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi;\n if n ops(res)=2 and Im(res[1])<>0 then \n RETURN(FAIL);\n else RETURN( op(res)) fi; \nend:\nfz2:=proc (x1,y1,z1,x2,y2,z2,x3,y3,z3,r) local re s;\n try\n res:=[allvalues(fz(x1,y1,z1,x2,y2,z2,x3,y3,z3,r))]\n cat ch: res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi;\n if \+ nops(res)=2 and Im(res[1])<>0 then \n RETURN(FAIL);\n else RETURN (op(res)) fi;\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "Testing th e results," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 149 "x1:=-0.5: y1 :=-1: z1:=0: # sphere 1\nx2:=1.75: y2:=1.5: z2:=0: # sphere 2\nx3:= 0: y3:=2.5: z3:=0: # sphere 3\n r:=1: # radius " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 121 "x4:=fx2(x1,y1,z1,x2,y2 ,z2,x3,y3,z3,r);\ny4:=fy2(x1,y1,z1,x2,y2,z2,x3,y3,z3,r);\nz4:=max(fz2( x1,y1,z1,x2,y2,z2,x3,y3,z3,r));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Viewing the Plot of Spheres" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 199 "c1:=sphere([x1,y1,z1], r):\nc2:=sphere([x2,y2,z2], r):\nc3:=s phere([x3,y3,z3], r):\nc4:=sphere([x4,y4,z4], r):\nplots[display]([c1, c2,c3,c4], scaling=constrained,\n style=patch, axes=boxe d);" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 63 "Geometric Model 2 (2 Sph eres and Fixed x to Support New Sphere)" }}{PARA 0 "" 0 "" {TEXT -1 210 "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 o f the two spheres. Return FAIL for the z coordinate if it is not feas ible." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 32 " Let: r = radius of the spheres " }}{PARA 0 "" 0 "" {TEXT -1 115 " \+ (x1,y1,z1), (x2,y2,z2) = origins of the two given spheres\n x 3 = known x location of the resting sphere" }}{PARA 0 "" 0 "" {TEXT -1 172 "Then, the origin of the resting sphere is the point where the \+ two given spheres would intersect if their diameters were doubled. Th e algebra is solved by Maple as follows:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 210 "unassign('x1','x2','x3','y1','y2','y3','z1','z2','z3 ','r','x','y','z'):\neqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:\neqn2:= (x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^2:\neqn3:=x=x3:\nsolution:=solve(\{eq n1,eqn2,eqn3\},\{x,y,z\}):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 95 "The solutions are very detailed expressions. They are converted to Maple functions as follows:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "a ssign(solution):\ngy:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r):\ngz:=unapply( z,x1,y1,z1,x2,y2,z2,x3,r):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "The functions may return multiple values. Also check for an infeasible s olution:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 518 "gy2:=proc (x1, y1,z1,x2,y2,z2,x3,r) local res;\n try\n res:=[allvalues(gy(x1,y1,z1, x2,y2,z2,x3,r))]\n catch: res:=FAIL;\n end try;\n if res=FAIL then \+ RETURN(FAIL) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n RETUR N(FAIL);\n else RETURN(op(res)) fi; \nend:\ngz2:=proc (x1,y1,z1,x2,y2 ,z2,x3,r) local res;\n try\n res:=[allvalues(gz(x1,y1,z1,x2,y2,z2,x3 ,r))]\n catch: res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL ) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n RETURN(FAIL);\n \+ else RETURN(op(res)) fi; \nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "Testing the results," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 156 "x1:=1.25: y1:=1: z1:=0: # sphere 1\nx2:=2.25: y2:=3: z2:=0: # sphere 2\nx3:=1: # sphere 3 x location\nr:=1: \+ # radius" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "y3:=gy2(x1,y1,z1,x2,y2,z2,x3,r);\nz3:=max(gz2(x1,y1,z1,x2,y2,z2,x3,r) );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Viewing the Plot of Spheres " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 169 "c1:=sphere([x1,y1,z1], r):\nc2:=sphere([x2,y2,z2], r):\nc3:=sphere([x3,y3,z3], r):\nplots[di splay]([c1,c2,c3], scaling=constrained,\n style=patch, a xes=boxed);" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 116 "Th e above solutions fail with a zero divide when y1=y2. Therefore, a ne w set of equations is derived for this case." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 829 "unassign('x1','x2','x3','y1','y2','y3','z1','z2 ','z3','r','x','y','z'):\neqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:\ne qn2:=(x-x2)^2+(y-y1)^2+(z-z2)^2=(2*r)^2:\neqn3:=x=x3:\nsolution:=solve (\{eqn1,eqn2,eqn3\},\{x,y,z\}):\nassign(solution):\ngy3:=unapply(y,x1, y1,z1,x2,y2,z2,x3,r):\ngz3:=unapply(z,x1,y1,z1,x2,y2,z2,x3,r):\ngy4:=p roc (x1,y1,z1,x2,y2,z2,x3,r) local res;\n try\n res:=[allvalues(gy3( x1,y1,z1,x2,y2,z2,x3,r))]\n catch: res:=FAIL;\n end try;\n if res=F AIL then RETURN(FAIL) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n \+ RETURN(FAIL);\n else RETURN(op(res)) fi; \nend:\ngz4:=proc (x1,y1 ,z1,x2,y2,z2,x3,r) local res;\n try\n res:=[allvalues(gz3(x1,y1,z1,x 2,y2,z2,x3,r))]\n catch: res:=FAIL;\n end try;\n if res=FAIL then R ETURN(FAIL) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n RETURN (FAIL);\n else RETURN(op(res)) fi; \nend:" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 62 "Another set is needed for the special case of y1=y2 and z1=z2." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 829 "unassign('x1',' x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'):\neqn1:=(x-x1) ^2+(y-y1)^2+(z-z1)^2=(2*r)^2:\neqn2:=(x-x2)^2+(y-y1)^2+(z-z1)^2=(2*r)^ 2:\neqn3:=x=x3:\nsolution:=solve(\{eqn1,eqn2,eqn3\},\{x,y,z\}):\nassig n(solution):\ngy5:=unapply(y,x1,y1,z1,x2,y2,z2,x3,r):\ngz5:=unapply(z, x1,y1,z1,x2,y2,z2,x3,r):\ngy6:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res ;\n try\n res:=[allvalues(gy5(x1,y1,z1,x2,y2,z2,x3,r))]\n catch: re s:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi;\n if nops(r es)=2 and Im(res[1])<>0 then \n RETURN(FAIL);\n else RETURN(op(re s)) fi; \nend:\ngz6:=proc (x1,y1,z1,x2,y2,z2,x3,r) local res;\n try\n res:=[allvalues(gz5(x1,y1,z1,x2,y2,z2,x3,r))]\n catch: res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi;\n if nops(res)=2 and \+ Im(res[1])<>0 then \n RETURN(FAIL);\n else RETURN(op(res)) fi; \+ \nend:" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 63 "Geometric Model 3 (2 \+ Spheres and Fixed y to Support New Sphere)" }}{PARA 0 "" 0 "" {TEXT -1 210 "Given two spheres located in 3D and fixed y origin location fo r a third sphere, find the origin of the third sphere that rests on to p of the two spheres. Return FAIL for the z coordinate if it is not f easible." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 32 "Let: r = radius of the spheres " }}{PARA 0 "" 0 "" {TEXT -1 115 " (x1,y1,z1), (x2,y2,z2) = origins of the two given spheres\n \+ y3 = known y location of the resting sphere" }}{PARA 0 "" 0 "" {TEXT -1 172 "Then, the origin of the resting sphere is the point wher e the two given spheres would intersect if their diameters were double d. The algebra is solved by Maple as follows:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 210 "unassign('x1','x2','x3','y1','y2','y3','z1','z2 ','z3','r','x','y','z'):\neqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:\ne qn2:=(x-x2)^2+(y-y2)^2+(z-z2)^2=(2*r)^2:\neqn3:=y=y3:\nsolution:=solve (\{eqn1,eqn2,eqn3\},\{x,y,z\}):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 95 "The solutions are very detailed expressions. They are converted t o Maple functions as follows:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "assign(solution):\nhx:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r):\nhz:=u napply(z,x1,y1,z1,x2,y2,z2,y3,r):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "The functions may return multiple values. Also check for an infea sible solution:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 520 "hx2:=pr oc (x1,y1,z1,x2,y2,z2,y3,r) local res;\n try\n res:=[allvalues(hx(x1 ,y1,z1,x2,y2,z2,y3,r))]\n catch: res:=FAIL;\n end try;\n if res=FAI L then RETURN(FAIL) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n \+ RETURN(FAIL);\n else RETURN(op(res)) fi; \nend:\nhz2:=proc (x1,y1, z1,x2,y2,z2,y3,r) local res;\n try\n res:=[allvalues(hz(x1,y1,z1,x2, y2,z2,y3,r))]\n catch: res:=FAIL;\n end try;\n if res=FAIL then RET URN(FAIL) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n RETURN(F AIL);\n else RETURN(op(res)) fi; \nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "Testing the results," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 156 "x1:=1.: y1:=1.3: z1:=0: # sphere 1\nx2:=2.5: y2:=1 .5: z2:=0: # sphere 2\ny3:=1: # sphere 3 y locatio n\nr:=1: # radius" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "x3:=hx2(x1,y1,z1,x2,y2,z2,y3,r);\nz3:=max(hz2(x1,y1,z 1,x2,y2,z2,y3,r));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Viewing the Plot of Spheres" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 169 "c1:=sp here([x1,y1,z1], r):\nc2:=sphere([x2,y2,z2], r):\nc3:=sphere([x3,y3,z3 ], r):\nplots[display]([c1,c2,c3], scaling=constrained,\n \+ style=patch, axes=boxed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 116 "T he above solutions fail with a zero divide when x1=x2. Therefore, a n ew set of equations is derived for this case." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 831 "unassign('x1','x2','x3','y1','y2','y3','z1','z2 ','z3','r','x','y','z'):\neqn1:=(x-x1)^2+(y-y1)^2+(z-z1)^2=(2*r)^2:\ne qn2:=(x-x1)^2+(y-y2)^2+(z-z2)^2=(2*r)^2:\neqn3:=y=y3:\nsolution:=solve (\{eqn1,eqn2,eqn3\},\{x,y,z\}):\nassign(solution):\nhx3:=unapply(x,x1, y1,z1,x2,y2,z2,y3,r):\nhz3:=unapply(z,x1,y1,z1,x2,y2,z2,y3,r):\nhx4:=p roc (x1,y1,z1,x2,y2,z2,y3,r) local res;\n try\n res:=[allvalues(hx3( x1,y1,z1,x2,y2,z2,y3,r))]\n catch: res:=FAIL;\n end try;\n if res=F AIL then RETURN(FAIL) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n \+ RETURN(FAIL);\n else RETURN(op(res)) fi; \nend:\nhz4:=proc (x1,y 1,z1,x2,y2,z2,y3,r) local res;\n try\n res:=[allvalues(hz3(x1,y1,z1, x2,y2,z2,y3,r))]\n catch: res:=FAIL;\n end try;\n if res=FAIL then \+ RETURN(FAIL) fi;\n if nops(res)=2 and Im(res[1])<>0 then \n RETUR N(FAIL);\n else RETURN(op(res)) fi; \nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 62 "Another set is needed for the special case of y1=y2 an d z1=z2." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 831 "unassign('x1', 'x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'):\neqn1:=(x-x1 )^2+(y-y1)^2+(z-z1)^2=(2*r)^2:\neqn2:=(x-x1)^2+(y-y2)^2+(z-z1)^2=(2*r) ^2:\neqn3:=y=y3:\nsolution:=solve(\{eqn1,eqn2,eqn3\},\{x,y,z\}):\nassi gn(solution):\nhx5:=unapply(x,x1,y1,z1,x2,y2,z2,y3,r):\nhz5:=unapply(z ,x1,y1,z1,x2,y2,z2,y3,r):\nhx6:=proc (x1,y1,z1,x2,y2,z2,y3,r) local re s;\n try\n res:=[allvalues(hx5(x1,y1,z1,x2,y2,z2,y3,r))]\n catch: r es:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi;\n if nops( res)=2 and Im(res[1])<>0 then \n RETURN(FAIL);\n else RETURN(op(r es)) fi; \nend:\nhz6:=proc (x1,y1,z1,x2,y2,z2,y3,r) local res;\n try \n res:=[allvalues(hz5(x1,y1,z1,x2,y2,z2,y3,r))]\n catch: res:=FAIL; \n end try;\n if res=FAIL then RETURN(FAIL) fi;\n if nops(res)=2 an d Im(res[1])<>0 then \n RETURN(FAIL);\n else RETURN(op(res)) fi; \+ \nend:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 66 "Geometric Model 4 (1 Sphere and Fixed (x,y) to Support Ne w Sphere)" }}{PARA 0 "" 0 "" {TEXT -1 168 "Given fixed (x,y) origin lo cations 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 ." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 32 "Let: r = radius of the spheres " }}{PARA 0 "" 0 "" {TEXT -1 107 " ( x1,y1,z1) = origin of the given spheres\n (x2,y2) = known (x,y) \+ location of the resting sphere" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 60 "The z-coordinate of the origin is found b y Maple as follows:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 140 "unas sign('x1','x2','x3','y1','y2','y3','z1','z2','z3','r','x','y','z'):\ne qn:=(x2-x1)^2+(y2-y1)^2+(z-z1)^2=(2*r)^2:\nsolution:=solve(eqn,\{z\}); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 "A Maple function is programme d for the calculation." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 159 " iz:=proc (x1,y1,z1,x2,y2,r) local res;\n res:=-x2^2+2*x1*x2-x1^2-y2^2 +2*y1*y2-y1^2+4*r^2;\n if res>0 then RETURN(z1+sqrt(res))\n else \+ RETURN(FAIL) fi; \nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "Testin g the results," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 161 "x1:=0: y 1:=0: z1:=0: # sphere origin\nx2:=0.5: y2:=0.5: # fixed (x,y) o f second sphere\nr:=1: # radius of spheres\nz2:=iz(x1 ,y1,z1,x2,y2,r);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Viewing the P lot of Spheres" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 139 "c1:=sphe re([x1,y1,z1], r):\nc2:=sphere([x2,y2,z2], r):\nplots[display]([c1,c2] , scaling=constrained,\n style=patch, axes=boxed);" }}}} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 48 "Geometric Model 5 (Distance Betwe en Two Spheres)" }}{PARA 0 "" 0 "" {TEXT -1 63 "Given two spheres in 3 D, find the distance between the origins." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 62 "Let (x1,y1,z1), (x2,y2,z2) = orig ins of the two given spheres" }{TEXT 256 1 "." }{TEXT -1 0 "" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "Distance:=(x1,y1,z1,x2,y2,z2 )->sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2);" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 22 "Miscellaneous Routines" }}{PARA 0 "" 0 "" {TEXT -1 23 "Ra ndom number generator" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "ran domize():\nrandomX:=stats[random,uniform[0,1]]('generator'):" }}}} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 10 "Box Module" }}{PARA 0 "" 0 "" {TEXT -1 25 "Module Box() Programming" }}{PARA 0 "" 0 "" {TEXT -1 168 " 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." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 12 "Definitions:" }}{PARA 0 "" 0 "" {TEXT -1 31 " Box:-width = width of the box" }}{PARA 0 "" 0 "" {TEXT -1 31 " \+ Box:-depth = depth of the box" }}{PARA 0 "" 0 "" {TEXT -1 33 " Box:-h eight = height of the box" }}{PARA 0 "" 0 "" {TEXT -1 63 " Box:-radiu s = radius of the spheres to be placed into the box" }}{PARA 0 "" 0 " " {TEXT -1 61 " Box:-SphereSet = set of data for the spheres within t he box" }}{PARA 0 "" 0 "" {TEXT -1 49 " Box:-Spheres = number of sphe res within the box" }}{PARA 0 "" 0 "" {TEXT -1 305 " Box:-fxbox, Box: -fybox, Box:-fzbox = routines that link with the calculations of secti on 'Geometric Model 1' for a sphere resting on three spheres\n Box:- gybox, Box:-gzbox = routines that link with the calculations of sectio n 'Geometric Model 2' for a sphere resting on two spheres and a wall ( fixed x)" }}{PARA 0 "" 0 "" {TEXT -1 155 " Box:-hxbox, Box:-hzbox = r outines that link with the calculations of section 'Geometric Model 3' for a sphere resting on two spheres and a wall (fixed y)" }}{PARA 0 " " 0 "" {TEXT -1 117 " Box:-makeSphere = routine for adding a sphere t o the set. It tests to make sure that the sphere location is valid." }}{PARA 0 "" 0 "" {TEXT -1 129 " Box:-checkDistances = routine for ch ecking the distances between a sphere and other spheres. It must be g reater than 2*radius." }}{PARA 0 "" 0 "" {TEXT -1 58 " Box:-checkDist ances2 = another distance checking routine" }}{PARA 0 "" 0 "" {TEXT -1 91 " Box:-firstLayer = routine for randomly placing the first laye r of spheres within the box." }}{PARA 0 "" 0 "" {TEXT -1 185 " Box:-m oreSpheres = routine for randomly adding spheres to the box. It rando mly picks three spheres from the set of spheres already within the box and tries to place a sphere on them." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26132 "module Box()\n description \"Box for randomly pla cing spheres\";\n local i;\n export width, depth, height, radius, ma keSphere, SphereSet,\n plot, firstLayer, Spheres, checkDistanc es, checkDistances2, moreSpheres,\n fxbox, fybox, fzbox, gybox , gzbox, hxbox, hzbox;\n # initialize parameters\n width:=6.+3./8.; \n depth:=3.;\n height:=5.+3./4.;\n radius:=3./4.;\n SphereSet:=\{ \};\n Spheres:=0;\n fxbox:=proc(T1::`name`,T2::`name`,T3::`name`)\n \+ description \"place a sphere on 3 spheres - calculate x of origin\" ;\n local x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, res, rand omList;\n x1:=T1:-x;\n x2:=T2:-x;\n x3:=T3:-x;\n y1:=T1:-y ;\n y2:=T2:-y;\n y3:=T3:-y;\n z1:=T1:-z;\n z2:=T2:-z;\n \+ z3:=T3:-z;\n try\n res:=fx2(x1,y1,z1,x2,y2,z2,x3,y3,z3,radius );\n catch: res:=FAIL;\n end try;\n if res=FAIL then RETURN (FAIL) fi; \n if nops([res])=1 then \n if type(res,numeric)= true then RETURN([res,res]) else RETURN(FAIL) fi;\n fi;\n # for \+ multiple real solutions, return both \n RETURN([res]);\n end proc: \n fybox:=proc(T1::`name`,T2::`name`,T3::`name`)\n description \"p lace a sphere on 3 spheres - calculate y of origin\";\n local x1, y 1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, res, randomList;\n x1:=T 1:-x;\n x2:=T2:-x;\n x3:=T3:-x;\n y1:=T1:-y;\n y2:=T2:-y; \n y3:=T3:-y;\n z1:=T1:-z;\n z2:=T2:-z;\n z3:=T3:-z;\n \+ try\n res:=fy2(x1,y1,z1,x2,y2,z2,x3,y3,z3,radius);\n catch: \+ res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi;\n i f nops([res])=1 then \n if type(res,numeric)=true then RETURN([re s,res]) else RETURN(FAIL) fi;\n fi;\n # for multiple real soluti ons, return both\n RETURN([res]);\n end proc:\n fzbox:=proc(T1::` name`,T2::`name`,T3::`name`)\n description \"place a sphere on 3 sp heres - calculate z of origin\";\n local x1, y1, z1, x2, y2, z2, x3 , y3, z3, x4, y4, z4, res, randomList;\n x1:=T1:-x;\n x2:=T2:-x; \n x3:=T3:-x;\n y1:=T1:-y;\n y2:=T2:-y;\n y3:=T3:-y;\n \+ z1:=T1:-z;\n z2:=T2:-z;\n z3:=T3:-z;\n try\n res:=fz2(x1 ,y1,z1,x2,y2,z2,x3,y3,z3,radius);\n catch: res:=FAIL;\n end tr y;\n if res=FAIL then RETURN(FAIL) fi;\n if nops([res])=1 then \+ \n if type(res,numeric)=true then RETURN([res,res]) else RETURN(F AIL) fi;\n fi;\n # for multiple solutions, return both\n RETU RN([res]);\n end proc:\n gybox:=proc(T1::`name`,T2::`name`,x3::`nume ric`)\n description \"place a sphere on 2 spheres and a x-wall - ca lculate y of origin\";\n local x1, y1, z1, x2, y2, z2, y3, z3, res, randomList;\n x1:=T1:-x;\n x2:=T2:-x;\n y1:=T1:-y;\n y2:= T2:-y;\n z1:=T1:-z;\n z2:=T2:-z;\n if y1=y2 and z1=z2 then\n \+ try\n res:=gy6(x1,y1,z1,x2,y2,z2,x3,radius);\n catc h: res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) fi; \n elif y1=y2 then\n try\n res:=gy4(x1,y1,z1,x2,y2,z2,x 3,radius);\n catch: res:=FAIL;\n end try;\n if res=FA IL then RETURN(FAIL) fi;\n else\n try\n res:=gy2(x1,y1, z1,x2,y2,z2,x3,radius);\n catch: res:=FAIL;\n end try;\n \+ if res=FAIL then RETURN(FAIL) fi;\n fi;\n if nops([res])=1 t hen\n if resmax(y1,y2) then res:=FAIL fi; \n \+ if type(res,numeric)=true then RETURN([res,res]) else RETURN(FAIL) fi;\n fi;\n # for multiple real solutions, return both\n res :=[res];\n if res[1]max(y1,y2) then res[1]:=F AIL fi;\n if res[2]max(y1,y2) then res[2]:=FA IL fi;\n if res[1]=FAIL and res[2]=FAIL then RETURN(FAIL) fi;\n \+ RETURN(res);\n end proc:\n gzbox:=proc(T1::`name`,T2::`name`,x3::`nu meric`)\n description \"place a sphere on 2 spheres and a x-wall - \+ calculate z of origin\";\n local x1, y1, z1, x2, y2, z2, y3, z3, re s, randomList;\n x1:=T1:-x;\n x2:=T2:-x;\n y1:=T1:-y;\n y2 :=T2:-y;\n z1:=T1:-z;\n z2:=T2:-z;\n if y1=y2 and z1=z2 then \n try\n res:=gz6(x1,y1,z1,x2,y2,z2,x3,radius);\n c atch: res:=FAIL;\n end try;\n if res=FAIL then RETURN(FAIL) \+ fi;\n elif y1=y2 then\n try\n res:=gz4(x1,y1,z1,x2,y2,z 2,x3,radius);\n catch: res:=FAIL;\n end try;\n if res =FAIL then RETURN(FAIL) fi;\n else\n try\n res:=gz2(x1, y1,z1,x2,y2,z2,x3,radius);\n catch: res:=FAIL;\n end try; \n if res=FAIL then RETURN(FAIL) fi;\n fi;\n if nops([res]) =1 then \n if type(res,numeric)=true then RETURN([res,res]) else \+ RETURN(FAIL) fi;\n fi;\n # for multiple real solutions, return b oth\n RETURN([res]);\n end proc:\n hxbox:=proc(T1::`name`,T2::`na me`,y3::`numeric`)\n description \"place a sphere on 2 spheres and \+ a y-wall - calculate x of origin\";\n local x1, y1, z1, x2, y2, z2, x3, z3, res, randomList;\n x1:=T1:-x;\n x2:=T2:-x;\n y1:=T1: -y;\n y2:=T2:-y;\n z1:=T1:-z;\n z2:=T2:-z;\n if x1=yx and \+ z1=z2 then\n try\n res:=hx6(x1,y1,z1,x2,y2,z2,y3,radius); \n catch: res:=FAIL;\n end try;\n if res=FAIL then RE TURN(FAIL) fi;\n elif x1=x2 then\n try\n res:=hx4(x1,y1 ,z1,x2,y2,z2,y3,radius);\n catch: res:=FAIL;\n end try;\n \+ if res=FAIL then RETURN(FAIL) fi;\n else\n try\n r es:=hx2(x1,y1,z1,x2,y2,z2,y3,radius);\n catch: res:=FAIL;\n \+ end try;\n if res=FAIL then RETURN(FAIL) fi;\n fi;\n if n ops([res])=1 then\n if resmax(x1,x2) then res: =FAIL fi; \n if type(res,numeric)=true then RETURN([res,res]) els e RETURN(FAIL) fi;\n fi;\n # for multiple real solutions, return both\n res:=[res];\n if res[1]max(x1,x2) \+ then res[1]:=FAIL fi;\n if res[2]max(x1,x2) t hen res[2]:=FAIL fi;\n if res[1]=FAIL and res[2]=FAIL then RETURN(F AIL) fi;\n RETURN(res);\n end proc:\n hzbox:=proc(T1::`name`,T2:: `name`,y3::`numeric`)\n description \"place a sphere on 2 spheres a nd a y-wall - calculate z of origin\";\n local x1, y1, z1, x2, y2, \+ z2, x3, z3, res, randomList;\n x1:=T1:-x;\n x2:=T2:-x;\n y1:= T1:-y;\n y2:=T2:-y;\n z1:=T1:-z;\n z2:=T2:-z;\n if x1=x2 a nd z1=z2 then\n try\n res:=hz6(x1,y1,z1,x2,y2,z2,y3,radius );\n catch: res:=FAIL;\n end try;\n if res=FAIL then \+ RETURN(FAIL) fi;\n elif x1=x2 then\n try\n res:=hz4(x1, y1,z1,x2,y2,z2,y3,radius);\n catch: res:=FAIL;\n end try; \n if res=FAIL then RETURN(FAIL) fi;\n else\n try\n \+ res:=hz2(x1,y1,z1,x2,y2,z2,y3,radius);\n catch: res:=FAIL;\n \+ end try;\n if res=FAIL then RETURN(FAIL) fi;\n fi;\n i f nops([res])=1 then \n if type(res,numeric)=true then RETURN([re s,res]) else RETURN(FAIL) fi;\n fi;\n # for multiple real soluti ons, return both\n RETURN([res]);\n end proc:\n makeSphere:=proc( f::`name`,xi::`complex`,yi::`complex`,\n zi::`comple x`)\n description \"add a sphere to the sphere set within the box\" ;\n local test;\n # check for valid (x,y) values\n #lprint(f, xi,yi,zi);\n if xi(width-radius) or \n yi(depth-radius) or \n zi(height-radius*0.5) or \n Im(xi)<>0 or Im(yi)<>0 or Im(zi)<>0 then\n RETURN( 'FAIL');\n fi; \n f:=module() export x, y, z; \n x:=xi; y :=yi; z:=zi; end module;\n Spheres:=Spheres+1;\n SphereSet:=Sphe reSet union \{f\};\n test:=checkDistances(Spheres);\n if test='F AIL' then \n SphereSet:=SphereSet minus \{f\};\n unassign(S| |Spheres);\n Spheres:=Spheres - 1;\n RETURN('FAIL');\n fi ;\n Spheres;\n end proc;\n plot:=proc(SList::list)\n descripti on \"create a 3D plot of the spheres in the box\"; \n local i,l1,l 2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12;\n l1:=plottools[line]([0,0,0],[ width,0,0],color=red);\n l2:=plottools[line]([width,0,0],[width,dep th,0],color=red);\n l3:=plottools[line]([width,depth,0],[0,depth,0] ,color=red);\n l4:=plott\\ools[line]([0,depth,0],[0,0,0],color=red) ;\nl5:=plottools[line]([0,0,height],[width,0,height],color=red);\nl6:= plottools[line]([width,0,height],[width,depth,height],color=red);\nl7: =plottools[line]([width,depth,height],[0,depth,height],color=red);\nl8 :=plottools[line]([0,depth,height],[0,0,height],color=red);\nl9:=plott ools[line]([0,0,0],[0,0,height],color=red);\nl10:=plottools[line]([wid th,0,0],[width,0,height],color=red);\nl11:=plottools[line]([width,dept h,0],[width,depth,height],color=red);\nl12:=plottools[line]([0,depth,0 ],[0,depth,height],color=red);\n if nargs=0 then\n for i from \+ 1 to nops(SphereSet) do\n c||i:=plottools[sphere]([SphereSet[i] :-x,SphereSet[i]:-y,\n SphereSet[i]:-z ], radius); \n end do;\n plots[display]([seq(c||i,i=1..nops( SphereSet)),\n l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12],\n \+ title=\"Spheres in Box\", scaling=constrained,\n style=patch, axes=boxed); \n else\n for i from 1 to nops(SList) do\n \+ c||i:=plottools[sphere]([SList[i]:-x,SList[i]:-y,\n \+ SList[i]:-z], radius);\n end do;\n plots[d isplay]([seq(c||i,i=1..nops(SList)),\n l1,l2,l3,l4,l5,l6,l7,l 8,l9,l10,l11,l12],\n title=\"Spheres in Box\", scaling=constr ained,\n style=patch, axes=boxed);\n fi;\n end proc;\n f irstLayer:=proc(times::posint)\n description \"randomly add spheres to the first layer\";\n local i, x, y, z, intList, randomList, ira ndom, Sphere1, movex, movey;\n z:=radius; # origin height for first layer\n for i from 1 to times do\n x:=randomX();\n x:=ra dius+x*(width-2.*radius); \n y:=randomX();\n y:=radius+y*(de pth-2.*radius);\n if modp(i,5)=0 then x:=radius fi;\n if mod p(i,10)=0 then x:=width-radius fi;\n if modp(i,4)=0 then y:=radiu s fi;\n if modp(i,8)=0 then y:=depth-radius fi;\n makeSphere (S||(Spheres+1),x,y,z);\n if modp(i,10)=0 then \n # pick \+ a random sphere from the box\n intList:=[seq(i,i=1..Spheres)]; \n randomList:=combinat[randperm](intList);\n irandom: =randomList[1];\n Sphere1:=S||irandom;\n # try moving \+ it to a new location toward a lower x value.\n movex:=randomX( )*0.5; # find the random move\n movey:=randomX()*0.25-0.125; \n x:=Sphere1:-x; # save the old location\n y:=Sphere 1:-y;\n # move the sphere\n movex:=x-movex;\n \+ if movex(depth- radius) then movey:=depth-radius fi;\n Sphere1:-x:=movex; \n \+ Sphere1:-y:=movey;\n #lprint(Sphere1,movex,movey,checkD istances(irandom),irandom,x,y);\n if checkDistances(irandom)=' FAIL' then # check that it is an ok move\n Sphere1:-x:=x; \n Sphere1:-y:=y;\n #lprint(\"FAIL\",Sphere1:-x, Sphere1:-y);\n fi; \n fi;\n # try each corner\n \+ makeSphere(S||(Spheres+1),radius,radius,radius);\n makeSphere( S||(Spheres+1),width-radius,radius,radius);\n makeSphere(S||(Sphe res+1),width-radius,depth-radius,radius);\n makeSphere(S||(Sphere s+1),radius,depth-radius,radius);\n Spheres; \n end do;\n end proc;\n checkDistances:=proc(SphereNumber::posint)\n description \+ \"check sphere for being too close to another sphere\";\n local i, \+ d;\n for i from 1 to Spheres do\n if i<>SphereNumber then\n \+ d:=Distance(S||i:-x,S||i:-y,S||i:-z,S||SphereNumber:-x,\n \+ S||SphereNumber:-y,S||SphereNumber:-z);\n d:=evalf( d);\n if (d+0.00001)<2.*radius then RETURN('FAIL') fi;\n f i;\n end do;\n RETURN(0); \n end proc;\n checkDistances2:=pro c(x::`numeric`,y::`numeric`,z::`numeric`)\n description \"check sph ere for being too close to another sphere\";\n local i, d;\n for i from 1 to Spheres do\n d:=Distance(S||i:-x,S||i:-y,S||i:-z,S ||SphereNumber:-x,\n S||SphereNumber:-y,S||SphereNu mber:-z);\n d:=evalf(d);\n if (d+0.00001)<2.*radius then RETURN('FAIL', S||i, d) fi;\n end do;\n RETURN(0); \n end pro c; \n moreSpheres:=proc(maxInBox::posint, maxTries::posint, trySpher es::list)\n description \"add more spheres to the box by placing th em on existing spheres and the walls\";\n local intList, randomList , Sphere1, Sphere2, Sphere3, i, imax, i1, i2, i3, \n x1, y1, \+ z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, res, oldSpheres, delta,\n \+ Pt1, Pt2, Pt3, Pt4, line1, line2, anglesum;\n oldSpheres:=Sphe res;\n if nargs=3 then \n imax:=1\n else\n imax:=maxTr ies\n fi; \n for i from 1 to imax do\n if Spheres>oldSph eres then \n lprint(Spheres);\n oldSpheres:=Spheres;\n fi;\n if nargs=2 then\n # randomly pick three sphere s from the box\n intList:=[seq(i,i=1..Spheres)];\n rando mList:=combinat[randperm](intList);\n i1:=randomList[1];\n \+ i2:=randomList[2];\n i3:=randomList[3];\n Sphere1:=S| |i1;\n Sphere2:=S||i2;\n Sphere3:=S||i3;\n #lprin t(Sphere1,Sphere2,Sphere3);\n else\n Sphere1:=trySpheres[1 ];\n Sphere2:=trySpheres[2];\n Sphere3:=trySpheres[3];\n fi;\n # will a sphere rest on the three random spheres?\n \+ z4:=fzbox(Sphere1,Sphere2,Sphere3);\n if z4<>FAIL then\n \+ x4:=fxbox(Sphere1,Sphere2,Sphere3);\n y4:=fybox(Sphere1,Sphe re2,Sphere3);\n fi;\n if z4<>FAIL and x4<>FAIL and y4<>FAIL \+ then\n # try solution #1\n # test that (x4,y4) of the so lution fails within the triangle of the\n # three supporting sp heres. Otherwise, the sphere is not fully supported.\n geometr y[point](Pt1,[Sphere1:-x,Sphere1:-y]);\n geometry[point](Pt2,[S phere2:-x,Sphere2:-y]);\n geometry[point](Pt3,[Sphere3:-x,Spher e3:-y]);\n geometry[point](Pt4,[x4[1],y4[1]]);\n try\n \+ anglesum:=geometry[FindAngle] \n (geome try[line](line1,[Pt4,Pt1]),geometry[line](line2,[Pt4,Pt2]))\n \+ +geometry[FindAngle]\n (geometry[line](line1,[ Pt4,Pt2]),geometry[line](line2,[Pt4,Pt3]))\n +geometry [FindAngle]\n (geometry[line](line1,[Pt4,Pt1]),geometr y[line](line2,[Pt4,Pt3]));\n catch: anglesum:=0;\n end try; \n if abs(anglesum-2.*evalf(Pi))<0.001 then \n \+ try\n res:=makeSphere(S||(Spheres+1),x4[1],y4[1],z4[1]); \n catch: res:=FAIL;\n end try;\n if res< >FAIL then \n lprint(3,Sphere1,Sphere2,Sphere3,S||Spheres) \n fi;\n fi;\n # try solution #2\n # tes t that (x4,y4) of the solution fails within the triangle of the\n \+ # three supporting spheres. Otherwise, the sphere is not fully sup ported.\n geometry[point](Pt4,[x4[2],y4[2]]);\n try\n \+ anglesum:=geometry[FindAngle] \n (geomet ry[line](line1,[Pt4,Pt1]),geometry[line](line2,[Pt4,Pt2]))\n \+ +geometry[FindAngle]\n (geometry[line](line1,[P t4,Pt2]),geometry[line](line2,[Pt4,Pt3]))\n +geometry[ FindAngle]\n (geometry[line](line1,[Pt4,Pt1]),geometry [line](line2,[Pt4,Pt3]));\n catch: anglesum:=0;\n end \+ try;\n if abs(anglesum-2.*evalf(Pi))<0.001 then \n tr y\n res:=makeSphere(S||(Spheres+1),x4[2],y4[2],z4[2]);\n \+ catch: res:=FAIL;\n end try;\n if res<>FAIL then \n lprint(3,Sphere1,Sphere2,Sphere3,S||Spheres)\n \+ fi;\n fi;\n fi;\n # try lots of positions every 3th time\n if nargs=3 or modp(i,3)=0 then\n # will a sphere rest on the first two random spheres and a wall?\n x4:=radius; \+ \n z4:=gzbox(Sphere1,Sphere2,x4);\n if z4<>FAIL then\n \+ y4:=gybox(Sphere1,Sphere2,x4);\n fi;\n if z4<>FAIL and y4< >FAIL then\n try\n res:=FAIL;\n if y4[1]<>FA IL and z4<>FAIL then\n res:=makeSphere(S||(Spheres+1),x4,y 4[1],z4[1]) fi;\n catch: res:=FAIL;\n end try;\n \+ if res<>FAIL then \n lprint(2,Sphere1,Sphere2,S||Sphere s)\n fi;\n try\n res:=FAIL;\n if y4[ 2]<>FAIL and z4[2]<>FAIL then\n res:=makeSphere(S||(Sphere s+1),x4,y4[2],z4[2]) fi;\n catch: res:=FAIL;\n end tr y;\n if res<>FAIL then \n lprint(2,Sphere1,Sphere2, S||Spheres)\n fi;\n fi;\n x4:=width-radius;\n z 4:=gzbox(Sphere1,Sphere2,x4);\n if z4<>FAIL then\n y4:=gyb ox(Sphere1,Sphere2,x4);\n fi;\n if z4<>FAIL and y4<>FAIL the n\n try\n res:=FAIL;\n if y4[1]<>FAIL and z4[ 1]<>FAIL then \n res:=makeSphere(S||(Spheres+1),x4,y4[1],z 4[1]) fi;\n catch: res:=FAIL;\n end try;\n if r es<>FAIL then \n lprint(2,Sphere1,Sphere2,S||Spheres)\n \+ fi;\n try\n res:=FAIL;\n if y4[2]<>FAIL an d z4[2]<>FAIL then\n res:=makeSphere(S||(Spheres+1),x4,y4[2 ],z4[2]) fi;\n catch: res:=FAIL;\n end try;\n i f res<>FAIL then \n lprint(2,Sphere1,Sphere2,S||Spheres)\n \+ fi;\n fi;\n y4:=radius;\n z4:=hzbox(Sphere1,Spher e2,y4);\n if z4<>FAIL then\n x4:=hxbox(Sphere1,Sphere2,y4) ;\n fi;\n if z4<>FAIL and x4<>FAIL then\n try\n \+ res:=FAIL;\n if x4[1]<>FAIL and z4[1]<>FAIL then\n \+ res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;\n cat ch: res:=FAIL;\n end try;\n if res<>FAIL then \n \+ lprint(2,Sphere1,Sphere2,S||Spheres)\n fi;\n fi;\n \+ if z4<>FAIL and x4<>FAIL then\n try\n res:=FAIL; \n if x4[2]<>FAIL and z4[2]<>FAIL then\n res:=make Sphere(S||(Spheres+1),x4[2],y4,z4[2]) fi;\n catch: res:=FAIL; \n end try;\n if res<>FAIL then \n lprint(2 ,Sphere1,Sphere2,S||Spheres)\n fi;\n fi; \n y4:=de pth-radius;\n z4:=hzbox(Sphere1,Sphere2,y4);\n if z4<>FAIL t hen\n x4:=hxbox(Sphere1,Sphere2,y4);\n fi;\n if z4<>F AIL and x4<>FAIL then\n try\n res:=FAIL;\n if x4[1]<>FAIL and z4[1]<>FAIL then\n res:=makeSphere(S||(Sph eres+1),x4[1],y4,z4[1]) fi;\n catch: res:=FAIL;\n end try;\n if res<>FAIL then \n lprint(2,Sphere1,Spher e2,S||Spheres)\n fi;\n try\n res:=FAIL;\n \+ if x4[2]<>FAIL and z4[2]<>FAIL then\n res:=makeSphere (S||(Spheres+1),x4[2],y4,z4[2]) fi;\n catch: res:=FAIL;\n \+ end try;\n if res<>FAIL then \n lprint(2,Spher e1,Sphere2,S||Spheres)\n fi;\n fi;\n # will a sphere rest on the second two random spheres and a wall?\n x4:=radius; \+ \n z4:=gzbox(Sphere2,Sphere3,x4);\n if z4<>FAIL then\n \+ y4:=gybox(Sphere2,Sphere3,x4);\n fi;\n if z4<>FAIL and y4 <>FAIL then\n try\n res:=FAIL;\n if y4[1]<>FA IL and z4[1]<>FAIL then\n res:=makeSphere(S||(Spheres+1),x4 ,y4[1],z4[1]) fi;\n catch: res:=FAIL;\n end try;\n \+ if res<>FAIL then \n lprint(2,Sphere2,Sphere3,S||Sphe res)\n fi;\n try\n res:=FAIL;\n if y 4[2]<>FAIL and z4[2]<>FAIL then\n res:=makeSphere(S||(Spher es+1),x4,y4[2],z4[2]) fi;\n catch: res:=FAIL;\n end t ry;\n if res<>FAIL then \n lprint(2,Sphere2,Sphere3 ,S||Spheres)\n fi;\n fi;\n x4:=width-radius;\n \+ z4:=gzbox(Sphere2,Sphere3,x4);\n if z4<>FAIL then\n y4:=gy box(Sphere2,Sphere3,x4);\n fi;\n if z4<>FAIL and y4<>FAIL th en\n try\n res:=FAIL;\n if y4[1]<>FAIL and z4 [1]<>FAIL then\n res:=makeSphere(S||(Spheres+1),x4,y4[1],z4 [1]) fi;\n catch: res:=FAIL;\n end try;\n if \+ res<>FAIL then \n lprint(2,Sphere2,Sphere3,S||Spheres)\n \+ fi;\n try\n res:=FAIL;\n if y4[2]<>FAI L and z4[2]<>FAIL then\n res:=makeSphere(S||(Spheres+1),x4, y4[2],z4[2]) fi;\n catch: res:=FAIL;\n end try;\n \+ if res<>FAIL then \n lprint(2,Sphere2,Sphere3,S||Spher es)\n fi;\n fi;\n y4:=radius;\n z4:=hzbox(Spher e2,Sphere3,y4);\n if z4<>FAIL then\n x4:=hxbox(Sphere2,Sph ere3,y4);\n fi;\n if z4<>FAIL and x4<>FAIL then\n try \n res:=FAIL;\n if x4[1]<>FAIL and z4[1]<>FAIL then \n res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;\n \+ catch: res:=FAIL;\n end try;\n if res<>FAIL then \n lprint(2,Sphere2,Sphere3,S||Spheres)\n fi;\n \+ try\n res:=FAIL;\n if x4[2]<>FAIL and z4[2]<>F AIL then\n res:=makeSphere(S||(Spheres+1),x4[2],y4,z4[2]) f i;\n catch: res:=FAIL;\n end try;\n if res<>F AIL then \n lprint(2,Sphere2,Sphere3,S||Spheres)\n \+ fi;\n fi; \n y4:=depth-radius;\n z4:=hzbox(Sphere2,Sph ere3,y4);\n if z4<>FAIL then\n x4:=hxbox(Sphere2,Sphere3,y 4);\n fi;\n if z4<>FAIL and x4<>FAIL then\n try\n \+ res:=FAIL;\n if x4[1]<>FAIL and z4[1]<>FAIL then\n \+ res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;\n c atch: res:=FAIL;\n end try;\n if res<>FAIL then \n \+ lprint(2,Sphere2,Sphere3,S||Spheres)\n fi;\n t ry\n res:=FAIL;\n if x4[1]<>FAIL and z4[1]<>FAIL the n\n res:=makeSphere(S||(Spheres+1),x4[1],y4,z4[1]) fi;\n \+ catch: res:=FAIL;\n end try;\n if res<>FAIL the n \n lprint(2,Sphere2,Sphere3,S||Spheres)\n fi;\n \+ fi;\n # will a sphere rest on the first and third random sphe res and a wall?\n x4:=radius; \n z4:=gzbox(Sphere1,Sphere3, x4);\n if z4<>FAIL then\n y4:=gybox(Sphere1,Sphere3,x4);\n fi;\n if z4<>FAIL and y4<>FAIL then\n try\n \+ res:=FAIL;\n if y4[1]<>FAIL and z4[1]<>FAIL then\n \+ res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;\n catch: res:=FAIL;\n end try;\n if res<>FAIL then \n \+ lprint(2,Sphere1,Sphere3,S||Spheres)\n fi;\n try\n \+ res:=FAIL;\n if y4[2]<>FAIL and z4[2]<>FAIL then\n \+ res:=makeSphere(S||(Spheres+1),x4,y4[2],z4[2]) fi;\n \+ catch: res:=FAIL;\n end try;\n if res<>FAIL then \n \+ lprint(2,Sphere1,Sphere3,S||Spheres)\n fi;\n f i;\n x4:=width-radius;\n z4:=gzbox(Sphere1,Sphere3,x4);\n \+ if z4<>FAIL then\n y4:=gybox(Sphere1,Sphere3,x4);\n fi; \n if z4<>FAIL and y4<>FAIL then\n try\n res:=FAI L;\n if y4[1]<>FAIL and z4[1]<>FAIL then\n res:=ma keSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;\n catch: res:=FAI L;\n end try;\n if res<>FAIL then \n lprint (2,Sphere1,Sphere3,S||Spheres)\n fi;\n try\n \+ res:=FAIL;\n if y4[1]<>FAIL and z4[1]<>FAIL then\n \+ res:=makeSphere(S||(Spheres+1),x4,y4[1],z4[1]) fi;\n catch: \+ res:=FAIL;\n end try;\n if res<>FAIL then \n \+ lprint(2,Sphere1,Sphere3,S||Spheres)\n fi;\n fi;\n \+ y4:=radius;\n z4:=hzbox(Sphere1,Sphere3,y4);\n if z4<>FAIL \+ then\n x4:=hxbox(Sphere1,Sphere3,y4);\n fi;\n if z4<> FAIL and x4<>FAIL then\n try\n res:=FAIL;\n i f x4[1]<>FAIL and z4[1]<>FAIL then\n res:=makeSphere(S||(Sp heres+1),x4[1],y4,z4[1]) fi;\n catch: res:=FAIL;\n en d try;\n if res<>FAIL then \n lprint(2,Sphere1,Sphe re3,S||Spheres)\n fi;\n try\n res:=FAIL;\n \+ if x4[2]<>FAIL and z4[2]<>FAIL then\n res:=makeSpher e(S||(Spheres+1),x4[2],y4,z4[2]) fi;\n catch: res:=FAIL;\n \+ end try;\n if res<>FAIL then \n lprint(2,Sphe re1,Sphere3,S||Spheres)\n fi;\n fi; \n y4:=depth-ra dius;\n z4:=hzbox(Sphere1,Sphere3,y4);\n if z4<>FAIL then\n \+ x4:=hxbox(Sphere1,Sphere3,y4);\n fi;\n if z4<>FAIL an d x4<>FAIL then\n try\n res:=FAIL;\n if x4[1] <>FAIL and z4[1]<>FAIL then\n res:=makeSphere(S||(Spheres+1 ),x4[1],y4,z4[1]) fi;\n catch: res:=FAIL;\n end try; \n if res<>FAIL then \n lprint(2,Sphere1,Sphere3,S| |Spheres)\n fi;\n try\n res:=FAIL;\n \+ if x4[2]<>FAIL and z4[2]<>FAIL then\n res:=makeSphere(S||( Spheres+1),x4[2],y4,z4[2]) fi;\n catch: res:=FAIL;\n \+ end try;\n if res<>FAIL then \n lprint(2,Sphere1,Sp here3,S||Spheres)\n fi;\n fi;\n fi;\n # try cor ner positions every 10th time\n if nargs=3 or modp(i,10)=0 then\n # check on resting a sphere on the first random sphere and a cor ner\n x1:=Sphere1:-x;\n y1:=Sphere1:-y;\n z1:=Sphere1:- z;\n x4:=radius;\n y4:=radius; \n z4:=iz(x1,y1,z1,x4,y4 ,radius);\n if type(z4,numeric)=true and (z4-z1)<1.999*radius the n \n res:=makeSphere(S||(Spheres+1),x4,y4,z4);\n if ty pe(res,numeric)=true then \n lprint(\"corner1\",Sphere1,S|| Spheres)\n fi;\n fi;\n x4:=radius;\n y4:=depth- radius; \n z4:=iz(x1,y1,z1,x4,y4,radius);\n if type(z4,numer ic)=true and (z4-z1)<1.999*radius then \n res:=makeSphere(S||( Spheres+1),x4,y4,z4);\n if type(res,numeric)=true then \n \+ lprint(\"corner2\",Sphere1,S||Spheres)\n fi;\n fi; \n x4:=width-radius;\n y4:=radius; \n z4:=iz(x1,y1,z1,x 4,y4,radius);\n if type(z4,numeric)=true and (z4-z1)<1.999*radius then \n res:=makeSphere(S||(Spheres+1),x4,y4,z4);\n i f type(res,numeric)=true then \n lprint(\"corner3\",Sphere1 ,S||Spheres)\n fi;\n fi;\n x4:=width-radius;\n \+ y4:=depth-radius; \n z4:=iz(x1,y1,z1,x4,y4,radius);\n if typ e(z4,numeric)=true and (z4-z1)<1.999*radius then \n res:=makeS phere(S||(Spheres+1),x4,y4,z4);\n if type(res,numeric)=true th en \n lprint(\"corner4\",Sphere1,S||Spheres)\n fi; \n fi;\n if Spheres>=maxInBox then RETURN(Spheres) fi;\n \+ fi;\n end do;\n RETURN(Spheres);\n end proc; \nend module;" } }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 29 "Applying the Stochastic Model" }}{PARA 0 "" 0 "" {TEXT -1 88 "The \+ Box module and supporting functions are used to model the random filli ng of the box." }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 11 "First Layer" }} {PARA 0 "" 0 "" {TEXT -1 292 "The spheres of the first layer are place d into the box by random selection of (x,y) coordinates. It takes man y 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. " }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 22 "Box:-firstLayer(5000);" }}}}{SECT 0 {PARA 4 "" 0 " " {TEXT -1 12 "Fill the Box" }}{PARA 256 "" 1 "" {TEXT -1 412 "The box is filled by placing more spheres on top of the existing spheres. Th ree spheres are randomly selected from the box and tested to see if sp here can be supported on top of them. On every 3rd time, tests are pe rformed to see if a sphere could rest upon two of the random spheres a nd 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." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "Box:-moreSpheres(40,2000);" }{TEXT -1 1 "\n" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 31 "Plot the Spheres Wi thin the Box" }}{PARA 0 "" 0 "" {TEXT -1 36 "Plotting all spheres with in the box," }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Box:-plot(); " }}{PARA 13 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}}{MARK "4" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }