PROGRAM FISSCH
C
C     CRITICAL MASS FOR FISSION CHAIN REACTION
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     Revised: Peter Signell (Mich. St. U.      3/02)
C
C     Compute the survival fraction F for a rectangular slab
C     of fissionable material.  F is defined as the number
C     of induced fissions resulting from one spontaneous
C     fission.  For F > 1.0  a chain reaction occurs. The
C     slab has the 'critical mass' when F = 1.  It is
C     assumed that two neutrons are emitted during each
C     fission.
C
C     SUBROUTINES: GETINF  - get input of M, S, & N
C                  CALC    - calculate survival fraction
C                  RANGEN  - generate 9 random nums for CALC
C                  DISPLA  - display the survival fraction
C                  QUERY   - ask if user wants to repeat
C                  
C                              F I S S C H 
C                                   :
C                  -----------------------------------
C                  :            :           :        :
C               GETINF         CALC      DISPLA    QUERY
C               (M,S,N)      (M,S,N)       (F)    (AGAIN)
C                               :
C                             RANGEN
C                             (IX,R)
C
C     ***** DECLARATIONS *****
      INTEGER IX
      REAL M,S,N,F
      CHARACTER AGAIN
      OPEN(UNIT=3,FILE='m356p1f.out')
      WRITE(3,'(/)')
      IX=13
C                    
C     ***** M A I N   B O D Y *****
C
 10   CALL GETINF (    M,S,N)
      CALL CALC   (M,S,N,IX,  F)
      CALL DISPLA (F        )
      CALL QUERY  (    AGAIN)
      IF(AGAIN.EQ.'Y') THEN
          WRITE(3,'(//)')
          GOTO 10
      ENDIF
      IF(AGAIN.EQ.'y') THEN
          WRITE(3,'(//)')
          GOTO 10
      ENDIF
      CLOSE(UNIT=3)
      STOP
      END
C
C     ***** G E T I N F *****
C
      SUBROUTINE GETINF(   M,S,N)
C    
C     Prompt user for M, S, & N
C
C     ***** DECLARATIONS *****
      REAL M,S,N
C
C     M  = mass of the rectangular slab 
C     S  = shape of slab (ratio of length to thickness)
C     N  = number of randomly generated neutrons
C
      PRINT *
      PRINT *,'Please supply M, S, & N values:'
      WRITE(*,20) 'M='
      READ *,M
      WRITE(*,20) 'S='
      READ *,S
      WRITE(*,20) 'N='
      READ *,N
  20  FORMAT(/A4)
      WRITE(*,24) M,S,N
      WRITE(3,24) M,S,N
  24  FORMAT(//'m356p1.for: FISSION CHAIN REACTION:'//
     +'               Mass=',F8.3/
     +'   Length/Thickness=',F8.3/
     +'    Random Fissions=',F8.3/)
      RETURN
      END
C                              
C     ***** C A L C *****
C
      SUBROUTINE CALC(M,S,N,IX,  F)
C  
C     calculate survival fraction
C
C     ***** DECLARATIONS *****
      INTEGER NUM,J,IX
      REAL M,S,N,R(9),A,B,NIN,PI,X0,Y0,Z0,PHI,COSTH
      REAL SINTH,D,X1,Y1,Z1,F
C
      PI=3.14159265
C
C     find dimensions of slab from M and S
C
      A=(M*S)**(1.0/3.0)
      B=(M/S**2)**(1.0/3.0)
C
C     NIN, the number of cases in which an emitted neutron
C     reacts with a neucleus inside the boundaries of the
C     slab is initialized to zero
C
      NIN=0.0
      NUM=N
C
C     generate N random fissions giving rise to 2N random
C     neutrons
C
      DO 40,J=1,NUM
C
C          each fission needs 9 random numbers
C
           CALL GENRAN(   IX,   R)
C
C          X, Y, and Z coordinates of the nucleus emitting
C          a neutron are random numbers inside the
C          boundaries of the slab.  X and Y range from -A/2
C          to +A/2, and Z ranges from -B/2 to +B/2
C
           X0=A*(R(1)-0.5)
           Y0=A*(R(2)-0.5)
           Z0=B*(R(3)-0.5)
C
C          two random neutrons from each fission:
C
           DO 30,K=1,2
C
C               two angles define direction in space
C
                PHI=2.0*PI*R(2*K+2)
                COSTH=2.0*(R(2*K+3)-0.5)
                SINTH=SQRT(1.0-COSTH**2)
C
C               distanct traveled by a neutron is assumed
C               to be a random number between 0 and 1
C
                D=R(K+7)
C
C               coordinates of interaction point for neutron
C
                X1=X0+D*SINTH*COS(PHI)
                Y1=Y0+D*SINTH*SIN(PHI)
                Z1=Z0+D*COSTH
C
C               if interaction pt is in dide slab, add 1 to
C                   NIN
C
                IF(ABS(X1).LE.A/2.0.AND.ABS(Y1).LE.A/2.0
     +                   .AND.ABS(Z1).LE.B/2.0) NIN=NIN+1.0
C
 30         CONTINUE
 40   CONTINUE
C
      F=NIN/N
      RETURN
      END
C
C     ***** D I S P L A *****
C
      SUBROUTINE DISPLA(F)
C
C     write the survival fraction F
C
C     ***** DECLARATIONS *****
      REAL F
C
      WRITE(*,44) F
      WRITE(3,44) F
 44   FORMAT(/'Survival Fraction F=',F8.3/)
      RETURN
      END
C
C     ***** G E N R A N *****
      SUBROUTINE GENRAN(   IX,   R)
C  
C     generate 9 random numbers
C
C     ***** DECLARATIONS *****
      REAL R(9),RAN
      INTEGER IX,K
C
      DO 50,K=1,9
           IX=IABS(899*IX)
           IX=IX-(32767+1)*(IX/(32767+1))
           RAN=IX
           R(K)=RAN/32767.0
 50   CONTINUE
      RETURN
      END
C
C     ***** Q U E R Y *****
      SUBROUTINE QUERY(    AGAIN)
C 
C     asks if user would like to try another
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
      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 *****