PROGRAM NEUTRONY 
 
   IMPLICIT NONE
 
   ! ------  PRZECHODZENIE NEUTRONOW PRZEZ PLYTE ROWNOLEGLOSCIENNA -----
 
   ! ------  DANE:
   real, parameter :: h=1.            ! grubosc plytki 
   real, parameter :: slam=3e-1*h     ! srednia droga swobodna 
   real, parameter :: pcap=0.2        ! prawdopodobienstwo pochlaniania
 
   integer, parameter :: ntraj=10**8  ! liczba trajektorii 
 
   real N0, Np, Nm, Nt, gam, zlam, xcos, z, sigma
   integer k
 
   ! pseudorandom seed
   CALL set_seed()
 
   sigma= 1.0/slam
 
   ! liczniki
   N0=0.; Np=0.; Nm=0.;
 
   ! trajektorie
   do k=1, ntraj
 
      ! cos kata rozproszenia dla wchodzacej trajektorii
      xcos=1.
 
      ! droga od brzegu plytki; glebokosc wnikania
      z=0.
 
      ! zdarzenia
1     continue
 
      ! przebyta droga (z rozkladu exp(-sigma*x))
      call random_number(gam)
      zlam=-slam * Log(gam)
      ! droga przebyta wzdluz osi z
      z=z+zlam*xcos
 
      ! przejscie na druga strone plytki
      if (z.ge.h) then
         Np=Np+1.
         cycle
      end if
 
      ! przejscie na lewa strone plytki
      if (z.le.0.) then
         Nm=Nm+1.
         cycle
      end if
 
      ! pochlanianie wewnatrz plytki
      call random_number(gam)
      if (gam.lt.pcap) then
         N0=N0+1.
         cycle
      end if
 
      ! rozpraszanie; nowy cosinus kata
      call random_number(gam)
      xcos=2.*gam-1.
      go to 1
 
   end do
 
   Nt=N0+Nm+Np
   write(*,"('Liczba n przechodzacych Np =',F10.0,f7.2,'%')") Np, Np/Nt*100
   write(*,"('Liczba n odbitych       Nm =',F10.0,f7.2,'%')") Nm, Nm/Nt*100
   write(*,"('Liczba n pochlonietych  N0 =',F10.0,f7.2,'%')") N0, N0/Nt*100
 
END PROGRAM NEUTRONY
 
 
 
SUBROUTINE set_seed()
  
  IMPLICIT NONE
 
  ! ----- variables for portable seed setting -----
  INTEGER :: i_seed
  INTEGER, DIMENSION(:), ALLOCATABLE :: a_seed
  INTEGER, DIMENSION(1:8) :: dt_seed
  ! ----- end of variables for seed setting -----
  
  ! ----- Set up random seed portably -----
  CALL RANDOM_SEED(size=i_seed)
  ALLOCATE(a_seed(1:i_seed))
  CALL RANDOM_SEED(get=a_seed)
  CALL DATE_AND_TIME(values=dt_seed)
  a_seed(i_seed)=dt_seed(8);
  a_seed(1)=dt_seed(8)*dt_seed(7)*dt_seed(6)
  CALL RANDOM_SEED(put=a_seed)
  DEALLOCATE(a_seed)
  ! ----- Done setting up random seed -----
   
END SUBROUTINE set_seed
