C       REACTOR SHIELDING PROBLEM
C
      REAL MU,E,AZM,D,DIST2C,X,Y,Z,EA,ER,ET,SA,SR,ST,THICK
      INTEGER EXIT,I,NA,NT,NR,NTOT,NPART
      LOGICAL ABSORB
      COMMON THICK
C
C
  100 WRITE(*,*)'ENTER TOTAL NUMBER OF PARTICLES, AND SLAB THICKNESS: '
      READ(*,*) NTOT,THICK
      IF(NTOT.LE.0)STOP
C       INITIALIZE PARTICLE COUNTS AND ENERGY TALLIES
      CALL INIT(NA,NT,NR,NPART,EA,ET,ER,SA,SR,ST)
C
C       MAIN LOOP OVER NTOT PARTICLES....................................
    1 CONTINUE
      NPART=NPART+1
      IF(NPART.EQ.NTOT)THEN
C               FINISHED, COMPUTE AND PRINT AVERAGES, STANDARD DEVIATIONS
        CALL OUTPUT(NA,EA,SA,NR,ER,SR,NT,ET,ST,NTOT)
        GOTO 100
      ENDIF
C
C       SOURCE GENERATES A NEW PARTICLE WITH POSITION, DIRECTION, ENERGY
      CALL SOURCE(E,MU,AZM,X,Y,Z)
C
    2 CONTINUE
C
C       COMPUTE-DISTANCE-TO-COLLISION
      D=DIST2C(E)
C
C        UPDATE POSITION TO DECIDE IF PARTICLE EXITS OR COLLIDES
      CALL UPDATE(MU,AZM,D,X,Y,Z)
C
      I=EXIT(X,Y,Z)
C       RETURNS -1, 0 , 1 FOR 'OUT ON LEFT', 'IN', 'OUT ON RIGHT'
      IF(I.LT.0)THEN
C               EXIT ON LEFT (REFLECTED), TALLY AND GOTO SOURCE
           NR=NR+1
           ER=ER+E  
           SR=SR+E*E
           GOTO 1
      ELSEIF(I.GT.0)THEN
C               EXIT ON RIGHT (TRANSMITTED), TALLY AND GOTO SOURCE
           NT=NT+1
           ET=ET+E
           ST=ST+E*E
           GOTO 1
      ELSEIF(I.EQ.0 .AND. ABSORB())THEN
C               COLLISION & ABSORBED, TALLY AND GOTO SOURCE
           NA=NA+1
           EA=EA+E
           SA=SA+E*E
           GOTO 1
      ELSE
C               COLLISION & SCATTERED, FIND SCATTERING ANGLE AND ENERGY
             CALL SCATT(MU,AZM,E)
C               GOTO COMPUTE-DISTANCE-TO-COLLISION
             GOTO 2
      ENDIF
      END
C
      SUBROUTINE INIT(NA,NT,NR,NPART,EA,ET,ER,SA,SR,ST)
         REAL EA,ER,ET,SA,SR,ST
         INTEGER NA,NT,NR,NPART
         NA=0
         NT=0
         NR=0
         NPART=0
         EA=0.
         ET=0.
         ER=0.
         SA=0.
         SR=0.
         ST=0.
         RETURN
      END
C
      SUBROUTINE SOURCE(E,MU,AZM,X,Y,Z)
C          EMITTER  
C          RETURNS MU UNIFORM ON [0,1] AS THE COSINE OF AN ANGLE, 
C              AND AZM UNIFORM ON [0,2*PI] AS ASIMUTHAL ANGLE 
        REAL E,MU,AZM,X,Y,Z,PI,UNI,ENERGY
        PI=4*ATAN(1.0)
        MU=UNI()
        AZM=UNI()*(2.0*PI)
        X=0.0
        Y=0.0
        Z=0.0
C          RETURNS ENERGY FROM SOURCE ENERGY DISTRIBUTION
        E=ENERGY()
        RETURN
      END
C
      SUBROUTINE UPDATE(MU,AZM,D,X,Y,Z)
         REAL MU,AZM,D,X,Y,Z,R      
         X=X+D*MU
         R=D*SQRT(1.-MU*MU)
         Y=Y+R*COS(AZM)
         Z=Z+R*SIN(AZM)
         RETURN
      END
C
      SUBROUTINE SCATT(MU,AZM,E)
C       RETURNS SCATTERING ANGLE AND ENERGY
        REAL MU,AZM,E,UNI,PI
C
C              ISOTROPIC SCATTERING, I.E., UNIFORM ON SPHERE 
        PI=4.0*ATAN(1.)
        MU=-1.+2.*UNI()
        AZM=UNI()*2*PI
C
C               FIND SCATTERING ENERGY, UNIFORM IN [.3*E,E]
         E=(UNI()*0.7+0.3)*E
         RETURN
      END
C
      LOGICAL FUNCTION ABSORB()
C         RETURNS TRUE IF PARTICLE IS ABSORBED. OCCURS WITH PROB PA.
        REAL PA,UNI
        PARAMETER(PA=0.1)
        IF(UNI().LE.PA)THEN
                ABSORB=.TRUE.
        ELSE
                ABSORB=.FALSE.
        ENDIF
        RETURN
      END
C
      INTEGER FUNCTION EXIT(X,Y,Z)
C        RETURNS -1, 0, +1 AS PARTICLE IS   OUTSIDE ON LEFT, INSIDE, 
C                            OR OUTSIDE ON RIGHT.
        REAL X,Y,Z,THICK
        COMMON THICK
        IF(X.GT.THICK)THEN
                EXIT=1
        ELSEIF(X.LT.0.0)THEN
                EXIT=-1
        ELSE
                EXIT=0
        ENDIF
        RETURN
      END
C
      REAL FUNCTION DIST2C(E)
C       RETURNS DISTANCE TO COLLISION, WITH EXPONENTIAL DISTRIBUTION  
C               WITH PARAMETER  `CROSS SECTION'
        REAL UNI,LOG,CROSS,E
        DIST2C=-LOG(UNI())/CROSS(E)
        RETURN
      END
C
      REAL FUNCTION ENERGY()
C      RETURNS ENERGY, WITH DISTRIBUTION CONST/SQRT(ENERGY)
C                               OVER [EMIN,EMAX]
C               USE INVERSE FUNCTION APPROACH TO COMPUTE THIS
C
C      ENERGY MIN, MAX IN MEV
        REAL C,SQRT,UNI,EMIN,EMAX
        PARAMETER(EMIN=1.0E-3,EMAX=2.5)
        C=1.0/(2.0*(SQRT(EMAX)-SQRT(EMIN)))
        ENERGY=( UNI()/(2*C)+SQRT(EMIN) )**2
        RETURN
      END
C
C
      REAL FUNCTION CROSS(E)
C          RETURNS CROSS SECTION (FICTIONAL) FOR ENERGY IN RANGE [EMIN,EMAX]
        REAL E,S,ABS,SIN,Y,EXP
        S=ABS(SIN(100.*(EXP(E)-1.))+SIN(18.81*(EXP(E)-1.)))
        Y=MAX(0.02,S)
        CROSS=10.0*EXP(-0.1/Y)
        RETURN
      END
C
      SUBROUTINE OUTPUT(NA,EA,SA,NR,ER,SR,NT,ET,ST,NTOT)
        REAL EA,SA,ER,SR,ET,ST
        INTEGER NA,NR,NT,NTOT
        WRITE(*,*) 'TALLIES '
        IF(NA.GT.0)THEN
           EA=EA/NA
           SA=SQRT(SA/NA-EA*EA)
        ENDIF
        WRITE(*,*) '% ABSORBED, ENERGY, SD: ',
     *       FLOAT(NA)/NTOT*100,EA,SA
        IF(NR.GT.0)THEN
           ER=ER/NR
           SR=SQRT(SR/NR-ER*ER)
        ENDIF
        WRITE(*,*) '% REFLECTED, ENERGY, SD: ',
     *       FLOAT(NR)/NTOT*100,ER,SR
        IF(NT.GT.0)THEN
           ET=ET/NT
           ST=SQRT(ST/NT-ET*ET)
        ENDIF
        WRITE(*,*) '% TRANSMITTED, ENERGY, SD: ',
     *       FLOAT(NT)/NTOT*100,ET,ST
        WRITE(*,*)
        RETURN
      END