Conference rusure::math

1 !Spheredot.bas Thomas E. Janzen ; calculates points on unit sphere with
!random distribution per spherical surface area.  For use with 3d renderer
!called render.c by the author
2 option angle = degrees ! do trigonometry in degrees instead of radians
5 randomize ! seed pseudo-random number generator
7 sign = 1	! initialize sign; this will toggle -1, 1 to get 2 hemispheres.
10 open 'spheredot.dat' FOR OUTPUT AS FILE #1%,&
	ORGANIZATION SEQUENTIAL,ALLOW WRITE,ACCESS WRITE
20 restore #1 ! open a file for the 3-D triples in text form
40 input "numpoints: ",numpoints ! input the number of points to find on sphere
90 print #1,numpoints ! tell the rendering program the number of 3d points
600 for n=1 to numpoints ! loop for finding the random location points
601	sign = -sign ! toggle sign of lattitude to get 2 hemispheres
605     longitude = (rnd * 360) - 180 !  find random longitude in degrees
608     axisradius=sqrt(1-((rnd-1)**2)) !find a random radius, dropped from axis
609     ! and bias the flat distribution with 90 degrees of a circle with 
	! center at (1,0) [y = sqrt(1-(x-1)**2)]
610	a = sqrt(1-(axisradius**2)) ! do arc cosine of axisradius
620     lat = atn(a/axisradius) * sign ! finish arc cosine and factor in sign
660     x = sin(longitude) * cos(lat) ! covert polar coordinates to rectilinear
670     y = cos(longitude) * cos(lat)
675     z = sin(lat)
680     print #1 using " #.##^^^^ ",x; ! print in text form to data file
685     print #1 using " #.##^^^^ ",y;
690     print #1 using " #.##^^^^",z
695 next n
1230 close #1%
1830 end

\

P1q#1;1;0;50;100#2;1;0;100;0#3;1;240;50;100#4;1;60;49;59#5;1;300;49;59#6;1;186;49;59#7;1;196;50;100#8;1;240;39;34#9;1;0;46;28#10;1;120;42;38#11;1;240;46;28#12;1;60;46;28#13;1;300;46;28#14;1;180;46;28#15;1;300;50;100------------#2!255?!11?_!16?G???OC!6?OG
G???C??_?C??C????_!6?G???O_??_!116?_!16?G?O?C!5?OG???G?C????_C??C!6?_???G???OO???_$-#2!254?_G??_G!5?@_????O`????C@!9?@??@?@_!10?C!12?_?C?A??G??_C?G?C?@?KAC!8?g??O_!79?_???GG??_???@???_????@O_???C@!5?@!6?@?@?_!9?C!11?A?_CG!7?c?GC@?ACAC!5?_G??_O$-#2!233?_
O?O???_A???C!11?@?KA???A???AC_!8?O????@????@???@?C??A!9?G!11?C_A!13?@!16?HA@_??A_??O?_!47?_oO??_??A!5?C!8?K??@?A?A?C?_??A????O!5?@???@???@?C!9?A!10?G!6?A???C?_!12?@!13?@?H@??_A??_O?_$-#2!223?_???_??@A???A!8?C??J???S???O???CA?O?AO!10?@!9?@???O??G???O?@??
?_!17?O!7?G!19?O!9?@A????G?@G??WG!30?_!6?@?A????A???C??B??OC?G!7?O????CAO??AO!9?@!9?@????O??G_??O?@!12?O!7?G!27?O!8?@@???GG@GG?G$-#2!218?O??A?_!7?AA?B??_!13?C?_G????O???A??@!5?C????@!15?@!10?C!14?@?A???O!13?A?QOG!6?O????c??`!17?@@?C!18?O??a!6?@A????a?@C
!17?C?_G????O???A?@!6?C????@!6?@!19?C!11?O??@?A!8?A??AG!6?PO????_?O?_??C??@!14?@@??C$-#2!215?A?a?O??_O!6?A!9?H!7?G!5?G??C!8?@???@!8?A!12?_!5?G?_!6?C!10?O!9?O??G!14?OO!8?O!23?FC?CG!11?AOA_!5?O?A!15?GA????G??K!12?@????@!12?A!9?G???_???C??_!7?O!20?O??G????
O!9?O!8?O!17?F?C???KG$-#2!210?O?K??CC?C_????O!7?A@!25?C@?C!9?G!16?G!19?O!6?O!15?C!12?C??_??C???G??C??O??C!7?A????A?O????_oE?C?OG?CC!5?C_???@?O?A!6?@!18?@!7?C??GC!26?G!19?O!6?O!5?C!15?_???C?CG??C???O!10?CA????A?O!5?__OA?CC$-#2!211?C???A!5?_????A?G!7?@@!1
4?@!7?C!5?_?GA??@!13?@A????G!12?G!24?g!6?_G????OA?O????O???O!8?G?@??A???G!14?@_C!6?A!6?_????A??H!8?@!5?@!13?_???C?@!5?GA!11?G????@A!17?G!24?g???A?O_G????O!7?O??O!7?G??H?A!16?`$-#2!210?@?@??_??AG?G!13?G???A!9?a?@!40?A!16?A!5?C!11?G???C??A!10?G!5?_!9?@?@!
7?A?O!8?B????AH@?@!5?_??A?G?G!13?G????A?_@!6?A!43?A!6?A????C!16?C??A??G!16?G!5?_?@?@!8?AO!14?B?A@G$-#2!214?OB??G??A????_??_C????_????_?A!11?A???_!9?A?@!9?_?A??OO!9?O!7?C!7?G????__!10?@@!5?@!8?G!5?@!10?O??A???G??O????O!5?pC??AA!7?OBA?G!6?C?_??_!5?a!5?_!1
3?A???_A!11?@?A???O???_????O!7?C??O!10?_????G????_?@!6?@??@!14?G!5?@!5?A???O???O??G?O!7?_?@PCAA$-#2!218?_O!5?C??C???PA???@?@G!18?_!14?A!20?A????@!8?_O!17?O!8?O??C!15?A!10?G???@?_!5?@G?C??@!11?_???O???C?C?@???@?OA!6?H!10?_!23?A!15?@!5?A!12?_O!17?OO?C!16?
A!15?@G!7?@???CG??@$-#2!219?F??CCC!6?G!5?C!10?O????@C@???P!12?@!12?O!7?GG???_!8?C!20?W!7?__??_!5?CC!7?G!13?@??G??G?@G?A@!16?@AECA?C!8?CC!13?@?CO??@O???@!17?@???O!6?G????_????G???G!20?G!8?O???_???__C!7?C?G!13?@????@?GG??G?@A@$-#2!229?C??OA??_?DO??G!13?A?
?_!11?C!8?C!6?C!30?_O???O!10?c???G!5?_O!6?O?G!6?oF???A??OGE???@!32?C???O??A?__GCP!17?A??_???C!9?C!5?C!36?O?_O!14?_C?_G!5?_OO?G!6?O?a??C?@??AW?E???@$-#2!236?A????GGO????_????_G?O????G!5?C!9?O!19?G_????A!5?G!10?O!10?O???O???@A??C!5?A!5?D?A?C?CA?@!52?A???G
O!5?__??g!5?O????CG!15?O!12?G?_???A!5?G!11?O!10?O!7?@O??C?A????A!5?C???HC??CA?@$-#2!247?@@?@????E?C?G!6?@_??_I?@!12?_?O!5?O!10?G???C!5?_???B???AO??O?@O??G???Q????@!5?A@!77?@?@????AC??CG!6?_?@g???A?@!8?_????O!6?O!5?G???C!9?_????AP??C@?o??G!5?O?A????@???A
?@$-#2!255?!18?@???A?A??A?A@!12?A?F!6?C@????A!6?A!132?@???A?A?A?A??@!11?@AC!6?D????C!7?A$-------------\

    A more efficient way to do this is via a random walk on the sphere.

    At each step, make a random plane rotation (the distribution of
    theta is not critical as long as it is not too pathological) and
    cyclically permute the XYZ coordinates of a unit vector.

    For example:

	x = y = z = 1.0/sqrt(3);		/* unit vector */

    loop:

	c = rand();				/* random number in [-1,1] */
	s = rand();
	t = 1.0/sqrt(c*c+s*s);
	c *= t;
	s *= t;					/* random sine and cosine */

	t = x;					/* cyclic permutation and */
	x = c*y - s*z;				/* random rotation	  */
	y = s*y + c*z;
	z = t;

	goto loop;

    This avoids the expense of transcendental function evaluation in the 
    loop.  You can obtain uniform distributions on any compact, symmetric
    manifold by similar techniques.

    You may want to renormalize the unit vector once in a while, say
    every few hundred points or so.

    Another simple alternative is to use rejection.  Generate 3 randoms
    uniform in the cube [-1,1]^3.  If the point lies outside the unit sphere
    try again (probability of sucess is a little over 1/2.)

    Then simply normalize to a unit vector (project it onto the
    surface of a unit sphere.)  But the random walk is faster.

    - Jim

    That's very interesting.  When I ran my old program in batch, it took
    3.44 seconds CPU time total (I don't know but I assume that includes
    this big fat setup file the sys mgr set up for batch files).  A BASIC
    version of your C code tool 4.44, so it seems less efficient.  Also, it
    does not return lattitude and longitude, which necessary in this
    problem; I apologize if that was not clear.  This is about random
    points on the earth.
    I will post the new program.  The picture from the new program
    looks about the same, but I
    will post it anyway.
    Tom

1 !Spheredot.bas Thomas E. Janzen ; calculates points on unit sphere with
!random distribution per spherical surface area.  For use with 3d renderer
!called render.c by the author
2 option angle = degrees ! do trigonometry in degrees instead of radians
5 randomize ! seed pseudo-random number generator
7 sign = 1	! initialize sign; this will toggle -1, 1 to get 2 hemispheres.
x = 1.0/sqrt(3)
y = 1.0/sqrt(3)
z = 1.0/sqrt(3)
10 open 'spheredot.dat' FOR OUTPUT AS FILE #1%,&
	ORGANIZATION SEQUENTIAL,ALLOW WRITE,ACCESS WRITE
20 restore #1 ! open a file for the 3-D triples in text form
40 input "numpoints: ",numpoints ! input the number of points to find on sphere
90 print #1,numpoints ! tell the rendering program the number of 3d points
600 for n=1 to numpoints ! loop for finding the random location points
	c = rnd				! random number in [-1,1] */
	s = rnd
	t = 1.0/sqrt(c*c+s*s)
	c = c * t
	s = s * t					! random sine and cosine */
	t = x					! cyclic permutation and */
	x = c*y - s*z				! random rotation	  */
	y = s*y + c*z
	z = t
680     print #1 using " #.##^^^^ ",x; ! print in text form to data file
685     print #1 using " #.##^^^^ ",y;
690     print #1 using " #.##^^^^",z
695 next n
1230 close #1%
1830 end

\

P1q#1;1;0;50;100#2;1;0;100;0#3;1;240;50;100#4;1;60;49;59#5;1;300;49;59#6;1;186;49;59#7;1;196;50;100#8;1;240;39;34#9;1;0;46;28#10;1;120;42;38#11;1;240;46;28#12;1;60;46;28#13;1;300;46;28#14;1;180;46;28#15;1;300;50;100------------#2!255?!10?_!5?O??OO?G???_O
!5?C???OC????cC_C???C????C?CG?G???OoW???O????_!119?O???OO?W?_O!7?C???CO??_?CC?C_??C????CC?GG??O_?OG???O????_????_$-#2!245?__??_!5?CCC??P???A@O????_O?A!10?@???@???W???@A!9?O??A?C!11?O?J!7?@??_??K????_OO!76?_??_!5?C?CC??@?O?@A??OO????_?A!6?@!7?@??@GO???A?
??O!9?AC!6?O!5?@I!10?_?CG???_OO$-#2!236?_?C?A??C_?P??@??O!5?C!5?_?A!20?C?A??G!8?Q???_!5?O!9?O!6?O?C???C!9?G!7?O!6?_????__O!48?_??CA?C!5?`OO?@???C????_!7?A!16?A!5?C????GA???_???O!9?O?O!8?C????S!14?G??O!11?OO_?O$-#2!229?Q??@_C??G?@D!19?C!5?O!22?G???@?A?D!
24?_????_S??`!8?A????G!5?_??@c!12?CG???G!35?AO_?@GG?C???@@!13?C!13?O!22?G??@??A?D!21?OC?_????_??@_A????G!6?_!7?@?c!9?GC??G$-#2!219?_?_@@??AA@??A!11?O??A!14?C!5?_!6?A????_!6?@!13?G!6?@!37?G!7?A!17?GO?O?Q@A??O!17?__?@@?A??AA@???A????O!9?A!13?_?C!7?_????A!
11?@!11?@??G!43?G!6?A!16?GG?OBA??O$-#2!215?GC!5?A??_!8?C?O??OO!41?O!17?G????A?H!20?G!10?_!14?C????@!9?_!5?_!6?IA!10?G?GC!5?_A!8?O?O??C!5?O!31?O!17?G!7?@!7?AG!21?G_!23?C???@!9?__!8?AI$-#2!212?E?_??O!15?@?A????O!8?G!11?_!8?G??@????_!12?C!13?G!14?@??C!5?O!
8?A!10?_!10?_!8?G?@???AAG??O!7?__??cA!5?O!10?@A????O!8?G!18?G?_!5?_????@!18?C!13?G!7?C????O?@!7?A!20?_?_!8?G??@????AH??Q!11?_$-#2!209?_K?GA??O??@?C!13?O????a!10?O!9?OP!14?_!16?_!15?C!8?A!7?AA???G!11?A_!5?_!21?@!5?_!5?_?OokGG???A?@CO!11?O????a!9?O!10?P!8
?O!15?_!16?_!13?A?C!16?AA??G??A!9?_!5?_!20?@!5?_??O?_Ko$-#2!213?C!8?@?O??A?o!8?O!11?A_???G!5?A!6?I???_!5?C!10?A???G?@!10?_@????A!46?@!12?C!13?ACpO@?C!6?@?O??A_!7?O!8?O??A?_!8?A??G!12?KA???_!6?A!5?@!7?G!12?_?@???A!45?@!5?C!14?A?@Osp$-#2!212?C?_???G!18?O!
7?D????I!6?C!18?@!10?O!9?G!12?CC?_G?@?C_A??G?A?A?O!15?_!6?G??`O?_???G???O??_@A???_@???@O??@@OT!5?c!6?G!12?O!6?@????GA??C???C!27?@!10?O!9?G???C????@?cAG?GgA!12?A?O!14?_?@????G?_?O?_??G@??O??_?A???@???@P@?CO@$-#2!214?@!10?C!6?C????@!16?O?O!5?A?_?G?@!6?OO?
??C@!18?_A!7?G!9?C!17?@!5?@O!5?O????O!7?W!12?P?@??S??C!8?@!13?C!7?C!5?@!9?O!6?O??G????A??_???@!6?O????C@!8?a!7?G!10?C!26?@!5?@O??O??O????OG!12?@????OC@?OCC$-#2!218?@A?KGO?c!6?A!7?gO????A!12?C????A????Q!10?G_!6?G!11?G!9?C?@!6?O!17?O!7?_?GG!14?CO@_A?O????
G?G?CCA!17?B?AG?S???_!7?A?GG????_!5?A!5?C!8?A???@????O??_!6?GG!10?G!18?S??@!23?O!9?GG!9?@?AC?O_??O???G?K?CA$-#2!227?@?GC?A?G_????C@o???AA????O??C!11?C!19?G?C!17?A!7?_!8?_?EKG!11?O!11?O_O!6?@????C???D?B!35?@KA!5?g?@?_?C?Q?_?OA??C!17?C!12?GC!17?A!16?_?A?G
G??_?C?C????O????@!7?O_O!7?@!7?C@??E?@$-#2!239?@????@???AO!12?A?G??A!9?C!5?a?@????G???C!22?_G!15?O?O?OO!9?_c????AI????C!56?@!9?O??A!13?A?G??A!8?AC!5?g?@?C!23?G!5?_!10?O?O??O????O!5?c_????A???I??C$-#2!254?@?A??G?H??OA??_?o!5?C???C????g@?g?I!11?G??G????_!
5?G!6?_OG?AO??_????K@OG??@???AA?A!85?B???H?G?AO??_?O_??C!9?D??ag??G?G!9?G?G???_!9?G??_O???OG?A??_?CG?GOH???@??AA@A$-#2!255?!17?@???A????A???AD!9?C????@@????CA!6?A!8?@!126?@????A????A???CA@!8?C??@@!6?E!7?A!7?@$-------------\

    Your program indeed finds random points on the surface of a sphere with 
equal probability of the point falling within any given equal area on the
surface of the sphere.

    Let's first change parameters.  Choose a random u uniformly within 
[-1,+1] and a random longitude v uniformly within [0,360[, or, equivalently, 
within [-180,+180[.  Then compute the point M(u,v) of coordinates:

	x = sqrt(1-u*u) * cos(v)
	y = sqrt(1-u*u) * sin(v)
	z = u

Note that sqrt(1-u*u) is the variable called axisradius in your
program, and that the latitude is the arc sinus of z=u=a.  So both 
methods are the same, up to notations.

    Now the density of points falling around the point M(u,v) is the area 
of the parallelogram generated by the variations of M, that is the vectors 
dM/du and dM/dv (where I intend to mean the partial derivatives).  This is
easily computed with a cross-product and gives:

         d� = || dM/du x dM/dv || = du x dv

Which is the uniform measure on the sphere since du and dv are the uniform
measure density.

    Alternatively, cut a small "rectangle" on the sphere of the form 
[u,u+du] x [v,v+dv].  Its area is allmost d� = du*dv provided that du 
and dv are small enough compared to the radius of the sphere (i.e. 
if you can neglect the curvature).  When generating a large number of 
points,the probability that a point has a u falling in [u,u+du] is P=du, 
the probability that a point has a v falling in [v,v+dv] is P=dv, and 
the probability that a point falls in the rectangle [u,u+du] x [v,v+dv] 
is P=du*dv.  Thus the number of points fallining within the small rectangle 
is proportional to its area.  As any measurable surface of the sphere is the 
limit of a covevering of such triangles, the property extend to any area, 
and the points' density is uniform on the sphere.

    One of the wierdnesses of BASIC I guess, the method in .3 *should* be
    much more efficient.  If what is desired is to get a whole bunch of
    points uniformly distributed over the unit sphere its about as good
    as you can do, I would say.  Its weakness is that there is a strong
    correlation between adjacent points (somewhat disguised by the
    "permutation" of the axises), which might be a problem for some uses
    (including graphical if only a few points are to be generated for each
    view).

    The standard way for generating uniform points on a three-sphere is
    rejection on the unit three-cube.  There are special methods but they
    only marginally improve on rejection in performance.  For dimensions
    greater than 3, rejection becomes more and more expensive (the ratio of
    the n-volume of the n-sphere to that of the n-volume of the n-cube goes
    to 0 as n goes to infinity).  For n above 5, the standard method is to
    generate n normally distributed random values, then scale them by the
    square-root of their sum of squares, to get the n coordinates needed.

    There might be a way to efficiently generate a correct latitude and
    longitude by using the normal method for n=3, and generating the
    normal variates by means of the polar method (which is really the
    inverse of this process for n=2).  Things may cancel out to produce
    something simple and efficient, but I haven't checked yet.

    In any case, given the BASIC wierdness you are dealing with, I can only
    suggest a simplification to your program.  If I haven't screwed up,
    you should be able to use:

	    sign = -sign
	    longitude = (rnd - .5)*360
	    latitude = arccos(rnd)

    which I'm pretty sure will come out to be uniform across the two
    hemispheres.  Probably easier for you to try it than for me to run
    sanity checks.

					Topher

    I am afraid that there must have been something wrong with the way that
    you measured these two programs, Tom.  I removed all of the output
    statements and measured them with Numpoints = 10000.  The resulsts were
    that your original program took 7.89 seconds on our system (8350) and
    the new version took only 1.02 seconds.(!)  So that is an increase of
    over 7x.
    
    --  Barry