10   '
20   '        MONTE CARLO EVALUATION OF INTEGRALS
30   '
40   '------------------- 1/5/88 -----------------------
50   '
60   '   Program to calculate the mass and center of mass 
70   '   coordinates of an object of soecified shape and
80   '   density using randomly located points.
90   '
100  '         --> LIST OF PROCEDURES <--
110  '
120        GOSUB  230  'INITIALIZATIONS
130        GOSUB  390  'INPUT
140        GOSUB  600  'SUM OVER NUM RANDOM-POINTS-PER-RUN
                       'FOR NRUNS
150        GOSUB  930  'OUTPUT
160        '
170        PRINT
180        INPUT "WANT TO RUN ANOTHER (Y/N)? ", AGAIN$
190        IF AGAIN$ = "Y" OR AGAIN$ = "y" THEN 130
200        '
210        END
220  '--------------------------------------------------
230  'INITIALIZATIONS
240  '
250  DIM X(3),INCOORD(3),RCM(3),DRCM(3),CMAV(3),AVCM2(3)
260  COUNT = 0
270  '
280  ' -- An odd integral value of IX is needed to
     '    initialize the random number generator.
300  '
310  IX = 13
320  IC = 1899
330  N  = 32768!
340  NM1 = N - 1
350  RHO = 1
360  '
370  RETURN
380  '--------------------------------------------------
390  'INPUT
400  '
410  CLS
420  PRINT "INPUT VALUES FOR THE FOLLOWING:"
430  PRINT "NUMBER OF RANDOMLY LOCATED POINTS TO USE (NUM)"
440  PRINT "NUMBER OF RUNS TO MAKE USING THIS VALUE OF NUM"
     PRINT "(NRUNS)"
450  PRINT "NUMBER OF RANDOM NUMBERS TO BE PRINTED ON"
     PRINT "SCREEN (NPTS)"
460  PRINT
470  PRINT "SEPARATE THE VALUES WITH A COMMA, PRESS RETURN"
480  PRINT
490  INPUT NUM,NRUNS,NPTS
500  '
510  MAV = 0
520  M2AV = 0
530  FOR J = 1 TO 3
540    CMAV(J) = 0
550    AVCM2(J) = 0
560    NEXT J
570  '
580  RETURN
590  '--------------------------------------------------
600  'SUM OVER NUM RANDOM POINTS/RUN FOR NRUNS
610  '
620  MAV=0: M2AV=0
     FOR K=1 TO 3: CMAV(K)=0: AVCM2(K)=0: NEXT K
630  '
640  FOR NR = 1 TO NRUNS
650    '
660    GOSUB 1310  'COMPUTE MSS AND RCM W. NUM RANDOM POINTS
670    '
680    ' -- Normalize and compute averages:
690    '
700    MSS = MSS * RHO / NUM
710    MAV = MAV + MSS / NRUNS
720    M2AV = M2AV + MSS ^ 2 / NRUNS
730    FOR K = 1 TO 3
740      RCM(K) = RCM (K) / NUM
750      CMAV(K) = CMAV(K) + RCM(K) / NRUNS
760      AVCM2(K) = AVCM2(K) + RCM(K) ^ 2 / NRUNS
770      NEXT K
780      '
790    PRINT USING " ###";NR;
800    PRINT USING " ##.###    ";MSS;RCM(1);RCM(2);RCM(3)
810    '
820    NEXT NR
830    '
840  '
850  DMSS =  SQR (M2AV - MAV ^ 2)
860  FOR K = 1 TO 3
870    DRCM(K) =  SQR (AVCM2(K) - CMAV(K) ^ 2)
880    NEXT K
890    '
900  '
910  RETURN
920  '--------------------------------------------------
930  'OUTPUT
940  '
950  INPUT "PRINTER (Y/N)? ", ANS$
960  IF ANS$ = "y" OR ANS$ = "Y" THEN 
965     OPEN "LPT1:" FOR OUTPUT AS #1
970  IF ANS$ <> "y" AND ANS$ <> "Y" THEN 
980     OPEN "SCRN:" FOR OUTPUT AS #1
990  PRINT #1, NUM;"  RANDOM POINTS"
1000 PRINT #1, NRUNS;"  RUNS"
1010 PRINT #1,
1020 PRINT #1, "              MASS      C.M. X      C.M. Y_
      C.M. Z"
1030 PRINT #1, "AVERAGES:  ";
1040 PRINT #1, USING " ##.###     ";MAV;CMAV(1);CMAV(2);_
CMAV(3)
1050 PRINT #1, "DEVIATIONS:";
1060 PRINT #1, USING " ##.###     ";DMSS;DRCM(1);DRCM(2);_
DRCM(3)
1070 '
1080 CLOSE #1
1090 '
1100 RETURN
1110 '--------------------------------------------------
1120 'GENERATE 3 RANDOM NUMBERS IN THE INTERVAL (-0.5,0.5)
1130 '
1140 FOR K = 1 TO 3
1150   '
1160   ' -- Generate random number RAN in the interval (0,1).
1170   '
1180   IX =  ABS (IC * IX)
1190   IX = IX - N * INT (IX / N)
1200   RAN = IX / NM1
1210   '
1220   ' -- End of random number generator.
1230   '
1240   ' -- Subtract 0.5 to get coordinates in (-0.5,+0.5).
1250   '
1260   X(K) = RAN - .5
1270   NEXT K
1280 '
1290 RETURN
1300 '--------------------------------------------------
1310 'COMPUTE MSS AND RCM USING NUM RANDOM POINTS
1320 '
1330 MSS = 0: FOR K = 1 TO 3: RCM(K) = 0: NEXT K
1340 '
1350 FOR J = 1 TO NUM
1360   '
1370  GOSUB 1120 'GENERATE 3 RANDOM NUMBERS IN THE INTERVAL
          '(-.5,.5)
1380  '
1390  GOSUB 1580 'PRINT THE RANDOM NUMBERS IF IN FIRST NPTS
                 'POINTS
1400  '
1410  GOSUB 1670 'DOES THE RANDOM POINT LIE INSIDE THE
                 'OBJECT?
1420  '
1430   ' -- If both conditions are satisfied, we use RHO=1
       '    and add this point's contribution to the mass,
       '    INSIDE, and to the 3 center-of-mass coordinates,
       '    RCM(K).
1460   '
1470   IF ANS$ = "OUTSIDE" THEN 1540
1480   '
1490   MSS = MSS + 1
1500   FOR K = 1 TO 3
1510     RCM(K) = RCM(K) + X(K)
1520     NEXT K
1530     '
1540   NEXT J
1550   '
1560 RETURN
1570 '--------------------------------------------------
1580 'PRINT THE RANDOM NUMBERS IF IN FIRST NPTS POINTS
1590 '
1600 IF COUNT = NPTS THEN COUNT = 100
1610 IF COUNT > NPTS THEN RETURN
1620 PRINT USING "#.#####       ";X(1);X(2);X(3)
1630 COUNT = COUNT + 1
1640 '
1650 RETURN
1660 '--------------------------------------------------
1670 'DOES THE RANDOM POINT LIE OUTSIDE THE OBJECT?
1680 '
1690 ' -- See if point lies outside a sphere of radius 0.5.
1700 '
1710 R1 = X(1) ^ 2 + X(2) ^ 2 + X(3) ^ 2
1720 IF R1 > .25 THEN ANS$ = "OUTSIDE": RETURN
1730 '
1740 ' -- See if point lies inside a cylinder of radius 0.3.
1750 '
1760 R2 = X(1) ^ 2 + X(2) ^ 2
1770 IF R2 < 9.000001E-02 THEN 
1780    
1790 RETURN
1800 '---------------- END OF PROGRAM LINES ----------------