!------------------- Shielding for Nuc. Reactors -----------------
!                 vs.2 (WRITE data on HD, many TRIALS)
! Program to simulate neutron behaviour in a shield using MC tech.
! if number of neutron collisions exeed a MAX, it is absorbed in
! the medium.  Otherwise it continues having more collisions.
! NHIST particles are entered into reactor shield, having random
! walks, finally emerged on either side of the slab sheild.
! FORTRAN M.Gharib
!--------------------------------------------------------------

integer :: e,w,n,s
call random_seed
print *," Number of Trials & Histories(particles) for each="
print *,"for example: 25, 1000"
read *,ntrial,nhist
print *,"Number of Max collisions beyond"
print *,"which particle is absorbed"
read *, ncoll
print *,"Shield Thickness (in terms of m.f.p.)"
read *,d
open (1,file="e:\shieldresult2.txt")
print     *,"NEUTRON RANDOM WALK THROUGH A SHIELD: MC Method"
print     *,"Shield Thick(mfp)=",d; print *,"No. of Trials=",ntrial
print     *,"No. of  Histories for each Trial=",nhist
print     *,"Exceeding",ncoll," collisions, particle absorbed"
print     *,"Ntransmit    Prob.   Nref      Prob.     Nabs     Prob."
write (1,*) "NEUTRON RANDOM WALK THROUGH A SHIELD: MC Method"
write (1,*) "Shield Thick(mfp)=",d; write(1,*)"No. of Trials=",ntrial
write (1,*) "No. of Histories for each Trial=",nhist
write (1,*) "Exceeding",ncoll," collisions, particle absorbed"
write (1,*) "Ntransmit    Prob.   Nref      Prob.     Nabs     Prob."
print *,repeat ("=",65) ; write(1,*)repeat ("=",65)

do k=1, ntrial
nt=0; nr=0 ;!-- RESET counters
do j=1, nhist
 e=0; w=0; n=0; s=0 ;!-- East, on penetration
 do i= 1 , ncoll
  call random_number (r)
  if (r < .25 ) then
    e=e+1
  else if(r <.50 ) then
    w=w+1
  else if(r < .75) then
    n=n+1
  else
    s=s+1
  end if
  x=e-w ; y=n-s
  if (x > d) then
    nt=nt+1
    exit
  end if  
  if (x < 0) then
    nr=nr+1
    exit
  end if  
 end do 
end do 
na=nhist-(nt+nr)
ptr=real (nt)/nhist
prf=real (nr)/nhist
pab=real (na)/nhist
print 33,   nt,ptr,nr,prf,na,pab
write (1,33)nt,ptr,nr,prf,na,pab
33 format (1x,3(i6,e13.3))
end do
close (1)
end
      
                
