This is a BASIC program to accompany the FAQ file on spheres. It demonstrates the method for placing points roughly uniformly on the sphere (using O(N) steps) and the method of improving the distribution based on the method of repulsion of charged particles (using O(N^2) steps per iteration -- far from maximally efficient!) This sample assumes you want to place 46 points on the sphere. I found really very little change after 100-200 iterations when eps=0.01 and K=6. We expect D to drop by an amount roughly equal to (the variable called Totgrad) * (eps). I never got D to drop below 1772.5, so I suspect that's close to the (local) minimum. I find eps=.01 to be appropriate, at least at first, since we are adding in a vector as large as length 30 or so when the initial points are only about pi/6 units apart. Observe that this program never stops, nor produces any output of coordinates! You'll have to insert some print statements and some statements allowing control either by the operator or on the basis of some measure of accuracy. 10 'Program to place some points evenly on the sphere. 20 'Written 1/21/95 by dave rusin, rusin@math.niu.edu 30 'Runs OK on UBASIC Basic compiler 40 ' 50 'Method: start with points evenly spaced on K evenly spaced rings 60 '(Here K=6 makes N=46 points). Then have points repel with inverse 70 ' square law. 80 ' 90 dim X(46),Y(46),Z(46):' points' cartesian coordinates 100 dim X1(46),Y1(46),Z1(46):' improved cartesian coordinates 110 K=6:' number of intervals (K-1 real rings, plus two poles) 120 Eps=+0.01:' well, initial spacing about pi/K, so... 130 Iter=1:' how many times points rearranged so far 140 ' 150 'set up the points in the first place 160 N=1:' will be number of points 170 X(N)=0:Y(N)=0:Z(N)=1 180 for I=1 to K-1 190 M=round(2*K*sin(#pi*(I/K))) 200 for J=0 to M-1 210 N=N+1 220 X(N)=sin(#pi*(I/K))*cos(#pi*(J/M)) 230 Y(N)=sin(#pi*(I/K))*sin(#pi*(J/M)) 240 Z(N)=cos(#pi*(I/K)) 250 next J 260 next I 270 N=N+1 280 X(N)=0:Y(N)=0:Z(N)=-1 290 'N is now the number of points (46) 300 ' 310 print "Current value of D is ";:' D=potential energy function 320 D=0 330 for I=1 to 46 340 for J=1 to 46 350 if J=I then 370 360 D=D+((X(I)-X(J))^2+(Y(I)-Y(J))^2+(Z(I)-Z(J))^2)^(-1/2) 370 next J 380 next I 390 print D:print 400 ' 410 print "Begining iteration #";Iter 420 Iter=Iter+1 430 ' 440 'compute total gradient in tangent planes to spheres 450 Totgrad=0 460 for I=1 to 46 470 'compute gradient (force) vector at i-th point 480 A=0:B=0:C=0 490 for J=1 to 46 500 if J=I then 550 510 F=((X(I)-X(J))^2+(Y(I)-Y(J))^2+(Z(I)-Z(J))^2)^(3/2) 520 A=A-(X(I)-X(J))/F:' first minus sign is part of derivative formula 530 B=B-(Y(I)-Y(J))/F 540 C=C-(Z(I)-Z(J))/F 550 next J 560 'subtract portion of gradient vector sticking out 570 F=A*X(I)+B*Y(I)+C*Z(I) 580 X1(I)=A-F*X(I) 590 Y1(I)=B-F*Y(I) 600 Z1(I)=C-F*Z(I) 610 Totgrad=Totgrad+X1(I)*X1(I)+Y1(I)*Y1(I)+Z1(I)*Z1(I) 620 next I 630 print "Total gradient size is ";sqrt(Totgrad) 640 'compute changed coordinates. First add in tangential change 650 for I=1 to 46 660 X1(I)=X(I)-X1(I)*Eps:' minus sign because we want D to DEcrease 670 Y1(I)=Y(I)-Y1(I)*Eps 680 Z1(I)=Z(I)-Z1(I)*Eps 690 'Next scale so we're back on the sphere 700 D=(X1(I)^2+Y1(I)^2+Z1(I)^2)^(1/2) 710 X1(I)=X1(I)/D 720 Y1(I)=Y1(I)/D 730 Z1(I)=Z1(I)/D 740 next I 750 'Now twist back so first two points are right. Pt 1= N pole; E pole on 760 'line thru first two (=(a,b,c)); third pole (d,e,f) is a cross product 770 F=X1(1)*X1(2)+Y1(1)*Y1(2)+Z1(1)*Z1(2) 780 A=X1(2)-X1(1)*F:B=Y1(2)-Y1(1)*F:C=Z1(2)-Z1(1)*F 790 F=sqrt(A*A+B*B+C*C) 800 A=A/F:B=B/F:C=C/F 810 D=Y1(1)*C-Z1(1)*B:E=Z1(1)*A-X1(1)*C:F=X1(1)*B-Y1(1)*A 820 for I=1 to 46 830 U=A*X1(I)+B*Y1(I)+C*Z1(I) 840 V=D*X1(I)+E*Y1(I)+F*Z1(I) 850 W=X1(1)*X1(I)+Y1(1)*Y1(I)+Z1(1)*Z1(I) 860 X1(I)=U:Y1(I)=V:Z1(I)=W 870 next I 880 'Find out how much we moved, THEN copy new coordinates in 890 Maxchg=0 900 for I=1 to 46 910 D=sqrt((X(I)-X1(I))^2+(Y(I)-Y1(I))^2+(Z(I)-Z1(I))^2) 920 if D>Maxchg then Maxchg=D:Changed=I 930 next I 940 'Just for fun, reorder points as you copy them back 950 for I=1 to 46 960 Test=-2 970 for J=1 to 46 980 if Z1(J)>Test then Test=Z1(J):Best=J 990 next J 1000 X(I)=X1(Best):Y(I)=Y1(Best):Z(I)=Z1(Best):Z1(Best)=-2 1010 next I 1020 ' 1030 'report on the changes, then loop back 1040 print "Largest change: vertex #";Changed;" changed by ";Maxchg 1050 goto 310 1060 ' 1070 'Should never get here! But you should alter the program so it does 1080 end