home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / prgramer / pli / sphere < prev    next >
Text File  |  1989-03-12  |  6KB  |  205 lines
  1. ClrScr;
  2. PRINT(SYSPRINT,'                        ',
  3.   'Evenly Spaced Points on a Sphere');
  4. PUTCRLF(SYSPRINT);
  5. PUTCRLF(SYSPRINT);
  6. PUTCRLF(SYSPRINT);
  7. num_points=0;
  8. DO WHILE(num_points < 2);
  9.   PRINT(SYSPRINT,'   How many points are to be placed?  ');
  10.   num_points=GETINT;
  11.   IF num_points < 2 THEN
  12.     DO;
  13.       PRINT(SYSPRINT,
  14.        '   ? The number of points must be no less 2');
  15.       PUTCRLF(SYSPRINT);
  16.     END;
  17. END;
  18. b=2.5;
  19. c=(b+1.0)*EXP(b*LOG(1.0+1.0/b));
  20. PUTCRLF(sysprint);
  21. Random=0.0;
  22. DO WHILE((Random <= 0.0) | (Random >= 1.0));
  23.   PRINT(SYSPRINT,'   Random number seed?  ');
  24.   Random=GETREAL;
  25.   IF ((Random <= 0.0) | (Random >= 1.0)) THEN
  26.     DO;
  27.       PRINT(SYSPRINT,
  28.        '   ? The seed must be between 0 and 1, exclusively');
  29.       PUTCRLF(SYSPRINT);
  30.     END;
  31. END;
  32. point_num_1=1;
  33. DO WHILE(point_num_1 <= num_points);
  34.   sum=0.0;
  35.   dimension_num=1;
  36.   DO WHILE(dimension_num <= 3);
  37.     tem_real=1.0;
  38.     DO WHILE(tem_real >= 1.0);
  39.       Random=c*Random*EXP(b*LOG(1.0-Random));
  40.       v1=2.0*Random-1.0;
  41.       Random=c*Random*EXP(b*LOG(1.0-Random));
  42.       v2=2.0*Random-1.0;
  43.       tem_real=v1*v1+v2*v2;
  44.     END;
  45.     tem_real=v1*SQRT(-2.0*LOG(tem_real)/tem_real);
  46.     sum=sum+tem_real*tem_real;
  47.     position(dimension_num)=tem_real;
  48.     dimension_num=dimension_num+1;
  49.   END;
  50.   sum=SQRT(sum);
  51.   dimension_num=1;
  52.   DO WHILE(dimension_num <= 3);
  53.     location(point_num_1,dimension_num)
  54.      =position(dimension_num)/sum;
  55.     dimension_num=dimension_num+1;
  56.   END;
  57.   point_num_1=point_num_1+1;
  58. END;
  59. max_delta=num_points;
  60. max_delta=SQRT(4.0*ATAN(1.0)/max_delta);
  61. current_potential=0.0;
  62. point_num_1=1;
  63. DO WHILE (point_num_1 < num_points);
  64.   point_num_2=point_num_1;
  65.   DO WHILE (point_num_2 < num_points);
  66.     point_num_2=point_num_2+1;
  67.     distance=0.0;
  68.     dimension_num=1;
  69.     DO WHILE(dimension_num <= 3);
  70.       distance=distance
  71.        +SQR(location(point_num_1,dimension_num)
  72.        -location(point_num_2,dimension_num));
  73.       dimension_num=dimension_num+1;
  74.     END;
  75.     distance=SQRT(distance);
  76.     current_potential=current_potential+1.0/distance;
  77.   END;
  78.   point_num_1=point_num_1+1;
  79. END;
  80. PUTCRLF(SYSPRINT);
  81. degrees_per_radian=180.0/PI;
  82. point_num_1=1;
  83. DO WHILE(point_num_1 <= num_points);
  84.   x=location(point_num_1,1);
  85.   y=location(point_num_1,2);
  86.   z=location(point_num_1,3);
  87.   latitude=ATAN(z/SQRT(SQR(x)+SQR(y)));
  88.   IF x >= 0.0 THEN
  89.     longitude=-PI/2+ATAN(y/x);
  90.   ELSE
  91.     longitude=PI/2+ATAN(y/x);
  92.   latitude=degrees_per_radian*latitude;
  93.   longitude=degrees_per_radian*longitude;
  94.   PRINT(SYSPRINT,latitude,' ',longitude);
  95.   PUTCRLF(SYSPRINT);
  96.   point_num_1=point_num_1+1;
  97. END;
  98. previous_potential=current_potential;
  99. better_min_found=TRUE;
  100. DO WHILE(better_min_found);
  101.   point_num_1=1;
  102.   DO WHILE(point_num_1 <= num_points);
  103.     dimension_num=1;
  104.     DO WHILE(dimension_num <= 3);
  105.       force(dimension_num)=0.0;
  106.       dimension_num=dimension_num+1;
  107.     END;
  108.     point_num_2=1;
  109.     DO WHILE(point_num_2 <= num_points);
  110.       IF point_num_1 != point_num_2 THEN
  111.         DO;
  112.           distance=0.0;
  113.           dimension_num=1;
  114.           DO WHILE(dimension_num <= 3);
  115.             distance=distance
  116.              +SQR(location(point_num_1,dimension_num)
  117.              -location(point_num_2,dimension_num));
  118.             dimension_num=dimension_num+1;
  119.           END;
  120.           distance=SQRT(distance);
  121.           dimension_num=1;
  122.           DO WHILE(dimension_num <= 3);
  123.             force(dimension_num)=force(dimension_num)
  124.              +(location(point_num_1,dimension_num)
  125.              -location(point_num_2,dimension_num))
  126.              /(distance*distance*distance);
  127.             dimension_num=dimension_num+1;
  128.           END;
  129.         END;
  130.       point_num_2=point_num_2+1;
  131.     END;
  132.     magnitude_of_force=0.0;
  133.     dimension_num=1;
  134.     DO WHILE(dimension_num <= 3);
  135.       magnitude_of_force
  136.        =magnitude_of_force+SQR(force(dimension_num));
  137.       dimension_num=dimension_num+1;
  138.     END;
  139.     magnitude_of_force=SQRT(magnitude_of_force);
  140.     magnitude_of_position=0.0;
  141.     dimension_num=1;
  142.     DO WHILE(dimension_num <= 3);
  143.       coordinate
  144.        =location(point_num_1,dimension_num)
  145.        +max_delta*force(dimension_num)/magnitude_of_force;
  146.       magnitude_of_position=magnitude_of_position
  147.        +coordinate*coordinate;
  148.       position(dimension_num)=coordinate;
  149.       dimension_num=dimension_num+1;
  150.     END;
  151.     magnitude_of_position=SQRT(magnitude_of_position);
  152.     dimension_num=1;
  153.     DO WHILE(dimension_num <= 3);
  154.       location(point_num_1,dimension_num)
  155.        =position(dimension_num)/magnitude_of_position;
  156.       dimension_num=dimension_num+1;
  157.     END;
  158.     point_num_1=point_num_1+1;
  159.   END;
  160.   current_potential=0.0;
  161.   point_num_1=1;
  162.   DO WHILE (point_num_1 < num_points);
  163.     point_num_2=point_num_1;
  164.     DO WHILE (point_num_2 < num_points);
  165.       point_num_2=point_num_2+1;
  166.       distance=0.0;
  167.       dimension_num=1;
  168.       DO WHILE(dimension_num <= 3);
  169.         distance=distance
  170.          +SQR(location(point_num_1,dimension_num)
  171.          -location(point_num_2,dimension_num));
  172.         dimension_num=dimension_num+1;
  173.       END;
  174.       distance=SQRT(distance);
  175.       current_potential=current_potential+1.0/distance;
  176.     END;
  177.     point_num_1=point_num_1+1;
  178.   END;
  179.   PUTCRLF(SYSPRINT);
  180.   point_num_1=1;
  181.   DO WHILE(point_num_1 <= num_points);
  182.     x=location(point_num_1,1);
  183.     y=location(point_num_1,2);
  184.     z=location(point_num_1,3);
  185.     latitude=ATAN(z/SQRT(SQR(x)+SQR(y)));
  186.     IF x >= 0.0 THEN
  187.       longitude=-PI/2+ATAN(y/x);
  188.     ELSE
  189.       longitude=PI/2+ATAN(y/x);
  190.     latitude=degrees_per_radian*latitude;
  191.     longitude=degrees_per_radian*longitude;
  192.     PRINT(SYSPRINT,latitude,' ',longitude);
  193.     PUTCRLF(SYSPRINT);
  194.     point_num_1=point_num_1+1;
  195.   END;
  196.   IF ABS(current_potential - previous_potential)/previous_potential
  197.    < 5.0E-6 THEN
  198.     better_min_found=FALSE;
  199.   ELSE
  200.     DO;
  201.       better_min_found=TRUE;
  202.       previous_potential=current_potential;
  203.     END;
  204. END;
  205.