      PROGRAM MCARLO
C
C     MONTE CARLO EVALUATION OF INTEGRALS
C     
C     -----------------------------------------------
C     AUTHOR: ROBERT EHRLICH (GEORGE MASON U    1973)
C     Revised: Peter Signell (Mich. St. U.      3/86)
C     Revised: Gary Herzenstiel (Mich. St. U.   3/86)
C     Revised: E.J.D.Kales (Mich. St. U.       11/87)
C
C     Calculates the mass and center of mass coordinates of
C     and object of specified shape and density using
C     randomly located points
C
C     SUBROUTINES: GETINF  - get N & NRUNS from the user
C                  CALC    - evaluate and print results
C                  RANGEN  - generate 3 random numbers for
C                            CALC
C                  QUERY   - ask if User wants to run the
C                            program again
C                  
C                          M C A R L O 
C                               :
C                  ------------------------------------
C                  :                    :             :
C               GETINF                CALC           QUERY
C          (N,NRUNS,IX,KOUNT)   (N,NRUNS,IX,KOUNT)  (AGAIN)
C                                       :
C                                     RANGEN
C                                     (IX,R)
C
C     ***** DECLARATIONS *****
      REAL N,NRUNS
      INTEGER IX,KOUNT
      CHARACTER AGAIN
      OPEN(UNIT=3,FILE='m355p1f.out')
C                    
C     ***** M A I N   B O D Y *****
C
      WRITE(3,'(/)')
 10   CALL GETINF (          N,NRUNS,IX,KOUNT)
      CALL CALC   (N,NRUNS,IX,KOUNT          )
      CALL QUERY  (                     AGAIN)
C     if user is done, END
      IF(AGAIN.EQ.'Y') THEN
          WRITE(3,'(A1)') 12
          GOTO 10
      ENDIF
      IF(AGAIN.EQ.'y') THEN
          WRITE(3,'(A1)') 12
          GOTO 10
      ENDIF
      CLOSE(UNIT=3)
      STOP
      END
C
C     ***** G E T I N F *****
C
      SUBROUTINE GETINF(   N,NRUNS,IX,KOUNT)
C   
C     Prompt user for input of N and NRUNS
C
C     ***** DECLARATIONS *****
      REAL N,NRUNS
      INTEGER IX,KOUNT
C
C     IX    = odd integral value needed to initialize RanGen
C     N     = number of randomly located points to use
C     NRUNS = number of runs to make using this value of N
C
      KOUNT=0
      IX=13
      PRINT *,'Please supply N & NRUNS values:'
      WRITE(*,20) '    N='
      READ *,N
      WRITE(*,20) 'NRUNS='
      READ *,NRUNS
  20  FORMAT(/A6)
      PRINT *
      PRINT *,'Montecarlo Evaluation of Integrals:'
      WRITE(*,'(A10,F5.0/A10,F5.0)') ' Points=',N,
     +'  No.Runs=',NRUNS
      WRITE(3,24) N,NRUNS
  24  FORMAT('m355p1f.for: MONTECARLO EVALUATION OF INTEGRALS'/
     +'   Points=',F5.0/
     +'  No.Runs=',F5.0/
     +' Coordinates of first 50 random points'//
     +'           X              Y              Z'/)
      RETURN
      END
C                              
C     ***** C A L C *****
C
      SUBROUTINE CALC(N,NRUNS,IX,KOUNT     )
C
C     does main calculations and prints results
C
C     ***** DECLARATIONS *****
      INTEGER IX,MRUNS,M,I,J,K,KOUNT
      REAL N,NRUNS,MASS,AVM,AVM2,AVCM(3),AVCM2(3),RCM(3)
      REAL RAN,R(3),DRCM(3),R1,R2,RHO,DM,DMASS
C
      MRUNS=NRUNS
      AVM=0.0
      AVM2=0.0
      AVCM(1)=0.0
      AVCM(2)=0.0
      AVCM(3)=0.0
      AVCM2(1)=0.0
      AVCM2(2)=0.0
      AVCM2(3)=0.0
C
C     LOOP OVER THE NUMBER OF RUNS
C
      DO 40,I=1,MRUNS
           M=N
           MASS=0.0
           RCM(1)=0.0
           RCM(2)=0.0
           RCM(3)=0.0
C
C          LOOP OVER THE NUMBER OF POINTS
C
           DO 32,J=1,M
C
C               3 coord's per point; get 3 random numbers
C
                CALL GENRAN(   IX,   R)
C
C               If < 51 points found, write out coordinates
C
                IF(KOUNT.LT.50)THEN
                     WRITE(3,'(3F14.6)') R(1),R(2),R(3)
                     KOUNT=KOUNT+1
                ELSEIF(KOUNT.EQ.50)THEN
                     WRITE(3,'(1X)')
                     KOUNT=100
                ENDIF
C
                R1=R(1)**2+R(2)**2+R(3)**2
C
C               see if point is < cylinder radius 0.5
C
                IF(R1.LE.0.25)THEN
                     R2=R(1)**2+R(2)**2
C
C                    see if point is > cylinder radius 0.3
C
                     IF(R2.GE.0.09)THEN
C
C                         if both conditions satisfied, use
C                         RHO=1.0 and add this point's
C                         contribution to mass and center of
C                         mass cooridinates (RCM(K)).
C
                          RHO=1.0
                          DM=RHO/N
                          MASS=MASS+DM
                          DO 30,K=1,3
                               RCM(K)=RCM(K)+R(K)*DM
  30                        CONTINUE
                     ENDIF
                ENDIF
  32       CONTINUE
C
C          compute average  MASS, MASS*2, RCM, RCM**2:
C
           AVM=AVM+MASS/NRUNS
           AVM2=AVM2+MASS**2/NRUNS
           DO 34,K=1,3
                RCM(K)=RCM(K)/MASS
                AVCM(K)=AVCM(K)+RCM(K)/NRUNS
                AVCM2(K)=AVCM2(K)+RCM(K)**2/NRUNS
 34        CONTINUE
           IF(I.EQ.1) WRITE(3,36) 12,N,NRUNS
 36        FORMAT(A1//'--------------------------------------------',
     +'--------------------------'/'  Points=',
     +F7.0/' No.Runs=',F7.0)
           WRITE(3,38) I,MASS,RCM(1),RCM(2),RCM(3)
 38        FORMAT(//' #',I3,'  Mass=',F10.5/
     +     7X,'C.M. coordinates  X=',F12.8,' Y=',F12.8,
     +' Z=',F12.8)
C
 40   CONTINUE
C
C     compute RMS deviations
C
      WRITE(3,'(//1X)')
      DMASS=SQRT(AVM2-AVM**2)
      DO 44,K=1,3
 44   DRCM(K)=SQRT(AVCM2(K)-AVCM(K)**2)
      WRITE(3,46) 'Averages  ',AVM,AVCM(1),AVCM(2),AVCM(3)
      WRITE(3,46) 'Deviations',DMASS,DRCM(1),DRCM(2),DRCM(3)
 46   FORMAT(/1X,
     +    A10,'  M=',F9.5,' X=',F12.8,' Y=',F12.8,
     +' Z=',F12.8)
      RETURN
      END
C
C     ***** G E N R A N *****
C
      SUBROUTINE GENRAN(   IX,   R)
C  
C     Generate 3 random numbers
C
C     ***** DECLARATIONS *****
      REAL R(3),RAN
      INTEGER IX,K
C
      DO 50,K=1,3
           IX=IABS(899*IX)
           IF(IX.GT.32767) IX=IX-(32767+1)*(IX/(32767+1))
           RAN=IX
           RAN=RAN/32767.0
C
C          subtract 0.5 to put coord's between -0.5 and +0.5
C
           R(K)=RAN-0.5
 50   CONTINUE
      RETURN
      END
C
C     ***** Q U E R Y *****
C
      SUBROUTINE QUERY(    AGAIN)
C 
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
C
      PRINT *
      WRITE(*,62)
 62   FORMAT(
     +'Would you like to do this again using new values?'/
     +'(Y/N)? ')
      READ( *,'(A1)')AGAIN
      RETURN
      END
C     ***** P R O G R A M    E N D *****