{ Particle Distribution in a Debye Cell © Charles Chandler http://qdl.scs-inc.us/?top=12797 } PROCEDURE Main; {This generates a random distribution of points around a center, for use in calculating the electric and gravitational forces between Debye cells. To make it realistic, we put more points nearer the center, and fewer further away (i.e., 80 of them within 25% of the radius of the cell, 40 of them within 50%, 20 of them within 75%, and another 10 around the outside.} VAR spherePts :ARRAY [1..80] OF VECTOR; tmpVect :VECTOR; cnt :INTEGER; str :STRING;     FUNCTION Rnd(input :REAL; places :INTEGER) :REAL; {Returns a real which is rounded to the number of places past the decimal point.} BEGIN Rnd := Round(input * (10 ^ places)) / (10 ^ places); END;     PROCEDURE GoldenSpiral(numPts :INTEGER; orbit, fuzz :REAL); {Get the 3D points on the surface of the sphere, using the Golden Spiral method. Note that this misses the last one, so we tell it to create (numPts - 1), then we create the last one manually, bringing the total up to numPts.} VAR inc, off, x, y, z, r, phi :REAL; BEGIN inc := PI * (3 - sqrt(5)); off := 2 / (numPts - 1); FOR cnt := 1 TO (numPts - 1) DO BEGIN y := (cnt * off) - 1 + (off / 2); r := SqRt(Abs(1 - (y * y))); phi := cnt * inc; spherePts[cnt].x := cos(phi) * r; spherePts[cnt].y := y; spherePts[cnt].z := sin(phi) * r; END; {Throw in the one that it misses.} spherePts[numPts].x := .05; spherePts[numPts].y := -1; spherePts[numPts].z := 0;   {Write the points to a text file, and create little spheres in the drawing. We multiply the unit vector by .3686 to scale the whole thing to be r = 1 m.} FOR cnt := 1 TO numPts DO BEGIN tmpVect := UnitVec(spherePts[cnt]) * (orbit + (fuzz * Random)) * .3686; WriteLn(Rnd(tmpVect.x, 8), Chr(9), Rnd(tmpVect.y, 8), Chr(9), Rnd(tmpVect.z, 8)); SetFillBack(CreateSphere(tmpVect.x, tmpVect.y, tmpVect.z, .02), 0, 0, 65535); END; END;     PROCEDURE FindSpherePts(h :HANDLE); BEGIN Get3DCntr(h, tmpVect.x, tmpVect.y, tmpVect.z); WriteLn(Rnd(tmpVect.x, 8), Chr(9), Rnd(tmpVect.y, 8), Chr(9), Rnd(tmpVect.z, 8)); END;     BEGIN {Call the GoldenSpiral procedure multiple times, to get varying densities of points as we move away from the center. The "orbit" argument sets the base distance from the center, and then the "fuzz" argument determines the range of randomization.} ReWrite('DebyeCellPoints.txt'); GoldenSpiral(80, 1.00, .25); GoldenSpiral(40, 1.50, .25); GoldenSpiral(20, 2.00, .25); GoldenSpiral(10, 2.50, .25); Close('DebyeCellPoints.txt');   {Put a sphere in the middle, for the dust grain.} SetFillBack(CreateSphere(0, 0, 0, .1), 65535, 0, 0);   {The VectorWorks drawing contains the HCP closest-packed arrangement of spheres, which was done just for the documentation of this code, but here we just query the drawing to get the centers of the spheres, so we don't have to calculate them.} ReWrite('ClosestPackedCenters.txt'); str := Concat('((L=''DebyeHCP'') & (T=95))'); ForEachObject(FindSpherePts, str); Close('ClosestPackedCenters.txt'); END; RUN(Main);