      PROGRAM SHIELD
************************************************************************
* This program uses the random number generator  RAND  to simulate     *
* neutrons entering a shield and to determine what percentage reaches  *
* the outside.  The neutrons are assumed to move forward, backward,    *
* left, and right with equal likelihood and to die within the shield   *
* if a certain number of changes of direction have occurred.           *
* Identifiers used are:                                                *
*     THICK  : thickness of shield                                     *
*     LIMDIR : limit on # of direction changes before energy dissipated*
*     RANNUM : a random real number in the range 0 to 1                *
*     NEWDIR : a random integer 1, 2, 3, or 4 representing direction   *
*     OLDDIR : previous direction of neutron                           *
*     NDIRS  : number of changes of direction                          *
*     FORW   : net units forward traveled                              *
*     N      : number of neutrons simulated                            *
*     COUNT  : number of neutrons reaching outside of shield           *
*     RANDOM : random number generator (subroutine)                    *
*     I      : subscript                                               *
*                                                                      *
* Input:   THICK, LIMDIR, N, and seed for random number generator      *
* Output:  Percentage of neutrons that reach the outside               *
************************************************************************

      INTEGER THICK, LIMDIR, NEWDIR, OLDDIR, NDIRS, FORW, N, COUNT, I
      REAL RANNUM

      PRINT *, 'ENTER THICKNESS OF SHIELD, LIMIT ON # OF DIRECTION'
      PRINT *, 'CHANGES, AND THE NUMBER OF NEUTRONS TO SIMULATE'
      READ *, THICK, LIMDIR, N
      PRINT *, 'SEED FOR RANDOM NUMBER GENERATOR'
      READ *, RANNUM
      COUNT = 0

* Begin the simulation

      DO 20 I = 1, N
         FORW = 0
         OLDDIR = 0
         NDIRS = 0

*        Repeat the following until neutron reaches outside of
*        shield, returns inside reactor, or dies within shield

 10      CONTINUE
            CALL RANDOM(RANNUM)
            NEWDIR = 1 + INT(4 * RANNUM)
            IF (NEWDIR.NE. OLDDIR) THEN
               NDIRS = NDIRS + 1
               OLDDIR = NEWDIR
            END IF
            IF (NEWDIR .EQ. 1) THEN
               FORW = FORW + 1
            ELSE IF (NEWDIR .EQ. 2) THEN
               FORW = FORW - 1
            END IF
         IF ((FORW .LT. THICK) .AND. (FORW .GE. 0) .AND.
     +       (NDIRS .LT. LIMDIR)) GO TO 10

         IF (FORW .EQ. THICK) THEN
            COUNT = COUNT + 1
         END IF
 20   CONTINUE

      PRINT 30, 100 * COUNT / REAL(N)
 30   FORMAT (1X, F5.2, '% OF THE NEUTRONS ESCAPED')

      END



**RANDOM****************************************************************
* This subroutine generates a random real number in the interval from  *
* 0 to 1.  The variable M is initially the seed supplied by the        *
* user;  thereafter, it is the random integer generated on the pre-    *
* ceding call to the function and saved using the SAVE statement.      *
* NOTE:  The constants 2147483647 and .4656613E-9 used in this sub-    *
* program are appropriate when it is executed on a machine having      *
* 32-bit memory words.  For a machine having M-bit words, these two    *
* constants should be replaced by the values of  2**M - 1  and         *
* 1/(2**M - 1), respectively.                                          *
************************************************************************

      SUBROUTINE RANDOM(RANNUM)

      INTEGER M, CONST1
      REAL RANNUM, CONST2
      PARAMETER (CONST1 = 2147483647, CONST2 = .4656613E-9)
      SAVE
      DATA M /0/

      IF (M .EQ. 0) M = INT(RANNUM)
      M = M * 65539
      IF (M .LT. 0) M = (M + 1) + CONST1
      RANNUM = M * CONST2

      END
