C       REACTOR SHIELDING PROBLEM
C
      DOUBLE PRECISION 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)
         DOUBLE PRECISION EA,ER,ET,SA,SR,ST
         INTEGER NA,NT,NR,NPART
         NA=0
         NT=0
         NR=0
         NPART=0
         EA=0.0D0
         ET=0.0D0
         ER=0.0D0
         SA=0.0D0
         SR=0.0D0
         ST=0.0D0
         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 
        DOUBLE PRECISION E,MU,AZM,X,Y,Z,PI,ENERGY,DUNI
        PI=4*ATAN(1.0D0)
        MU=DUNI()
        AZM=DUNI()*(2.0D0*PI)
        X=0.0D0
        Y=0.0D0
        Z=0.0D0
C          RETURNS ENERGY FROM SOURCE ENERGY DISTRIBUTION
        E=ENERGY()
        RETURN
      END
C
      SUBROUTINE UPDATE(MU,AZM,D,X,Y,Z)
         DOUBLE PRECISION MU,AZM,D,X,Y,Z,R      
         X=X+D*MU
         R=D*SQRT(1.0D0-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
        DOUBLE PRECISION MU,AZM,E,PI,DUNI
C
C              ISOTROPIC SCATTERING, I.E., UNIFORM ON SPHERE 
        PI=4.0D0*ATAN(1.0D0)
        MU=-1.0D0+2.0D0*DUNI()
        AZM=DUNI()*2D0*PI
C
C               FIND SCATTERING ENERGY, UNIFORM IN [.3D0*E,E]
         E=(DUNI()*0.7D0+0.3D0)*E
         RETURN
      END
C
      LOGICAL FUNCTION ABSORB()
C         RETURNS TRUE IF PARTICLE IS ABSORBED. OCCURS WITH PROB PA.
        DOUBLE PRECISION PA,DUNI
        PARAMETER(PA=0.1D0)
        IF(DUNI().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.
        DOUBLE PRECISION X,Y,Z,THICK
        COMMON THICK
        IF(X.GT.THICK)THEN
                EXIT=1
        ELSEIF(X.LT.0.0D0)THEN
                EXIT=-1
        ELSE
                EXIT=0
        ENDIF
        RETURN
      END
C
      DOUBLE PRECISION FUNCTION DIST2C(E)
C       RETURNS DISTANCE TO COLLISION, WITH EXPONENTIAL DISTRIBUTION  
C               WITH PARAMETER  `CROSS SECTION'
        DOUBLE PRECISION LOG,CROSS,E,DUNI
        DIST2C=-LOG(DUNI())/CROSS(E)
        RETURN
      END
C
      DOUBLE PRECISION 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
        DOUBLE PRECISION C,SQRT,EMIN,EMAX,DUNI
        PARAMETER(EMIN=1.0D-3,EMAX=2.5D0)
        C=1.0D0/(2.0D0*(SQRT(EMAX)-SQRT(EMIN)))
        ENERGY=( DUNI()/(2*C)+SQRT(EMIN) )**2
        RETURN
      END
C
C
      DOUBLE PRECISION FUNCTION CROSS(E)
C          RETURNS CROSS SECTION (FICTIONAL) FOR ENERGY IN RANGE [EMIN,EMAX]
        DOUBLE PRECISION E,S,ABS,SIN,Y,EXP
        S=ABS(SIN(100.0D0*(EXP(E)-1.0D0))+SIN(18.81D0*(EXP(E)-1.0D0)))
        Y=MAX(0.02D0,S)
        CROSS=10.0D0*EXP(-0.1D0/Y)
        RETURN
      END
C
      SUBROUTINE OUTPUT(NA,EA,SA,NR,ER,SR,NT,ET,ST,NTOT)
        DOUBLE PRECISION 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(*,'(1X,A,D10.4,2D20.10)') '% 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(*,'(1X,A,D10.4,2D20.10)') '% 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(*,'(1X,A,D10.4,2D20.10)') '% TRANSMITTED, ENERGY, SD: ',
     *       FLOAT(NT)/NTOT*100,ET,ST
        WRITE(*,*)
        RETURN
      END