10   'Program to play 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
  130   Iter=1:'                        how many times we arranged the points
  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-1772.5: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   '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!
 1080   end