program main
!*****************************************************************************80
!
!! reactor_simulation() simulations the reactor shielding problem.
!
! Discussion:
!
! This is a Monte Carlo simulation, using
! uniform random numbers, which investigates the
! effectiveness of a shield intended to absorb the
! neutrons emitted from a nuclear reactor.
!
! The reactor is modeled as a point source,
! located at (0,0,0).
!
! A particle emitted from the reactor has a random
! initial direction, and an energy selected from
! [Emin,Emax] with a 1/Sqrt(E) distribution.
!
! The shield is modeled as a wall of thickness THICK,
! extending from 0 to THICK in the X direction, and
! extending forever in the Y and Z directions.
!
! Based on the particle energy, a distance D is computed
! which measures how far the particle could travel through
! the shield before colliding.
!
! Based on the particle direction, the position is updated
! by D units.
!
! If the particle is now to the left of the shield, it is
! counted as being REFLECTED.
!
! If the particle is to the right of the shield, it is
! counted as being ABSORBED.
!
! If the particle is inside the shield, it has COLLIDED.
! A particle that collides is either absorbed (end of story)
! or SCATTERED with a new random direction and a new (lower)
! energy.
!
! Every particle is followed from origin to its final fate,
! which is reflection, transmission, or absorption.
! At the end, a summary is printed, giving the number of
! particles with each fate, and the average energy of each
! group of particles.
!
! Increasing NTOT, the number of particles used, will improve the
! expected reliability of the results.
!
! Increasing THICK, the thickness of the shield, should
! result in more absorptions and reflections.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Local Parameters:
!
! Local, real ( kind = rk ) AZM, the azimuthal angle of the particle's
! direction.
!
! Local, real ( kind = rk ) D, the distance that the particle can
! travel through the slab, given its current energy.
!
! Local, real ( kind = rk ) E, the energy of the particle.
!
! Local, real ( kind = rk ) EA, energy absorbed by the slab.
!
! Local, real ( kind = rk ) ER, energy reflected by the slab.
!
! Local, real ( kind = rk ) ET, energy transmitted through the slab.
!
! Local, real ( kind = rk ) MU, the cosine of the angle between the
! particle's direction and the X axis.
!
! Local, integer NA, number of particles absorbed by the slab.
!
! Local, integer NPART, the index of the current particle.
!
! Local, integer NR, number of particles reflected by the slab.
!
! Local, integer NT, number of particles transmitted
! by the slab.
!
! Local, integer NTOT, the total number of particles to be
! emitted from the neutron source.
!
! Local, real ( kind = rk ) SA, standard deviation of absorbed energy.
!
! Local, real ( kind = rk ) SR, standard deviation of reflected energy.
!
! Local, real ( kind = rk ) ST, standard deviation of transmitted energy.
!
! Local, real ( kind = rk ) THICK, the thickness of the slab that is
! intended to absorb most of the particles.
!
! Local, real ( kind = rk ) X, Y, Z, the current position of the particle.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
logical absorb
real ( kind = rk ) azm
real ( kind = rk ) d
real ( kind = rk ) dist2c
real ( kind = rk ) e
real ( kind = rk ) ea
real ( kind = rk ) er
real ( kind = rk ) et
real ( kind = rk ) mu
integer na
integer npart
integer nr
integer nt
integer, parameter :: ntot = 100000
real ( kind = rk ) sa
real ( kind = rk ) sr
real ( kind = rk ) st
integer test
integer, parameter :: test_num = 5
real ( kind = rk ), parameter :: thick = 2.0D+00
real ( kind = rk ) x
real ( kind = rk ) y
real ( kind = rk ) z
call timestamp ( )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'reactor_simulation():'
write ( *, '(a)' ) ' FORTRAN90 version'
write ( *, '(a)' ) ' The reactor shielding simulation.'
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' Shield thickness is THICK = ', thick
write ( *, '(a,i8)' ) ' Number of simulated particles is NTOT = ', ntot
write ( *, '(a,i8)' ) ' Number of tests TEST_NUM = ', test_num
do test = 1, test_num
write ( *, '(a)' ) ' '
write ( *, '(a,i2)' ) ' Test # ', test
!
! Initialize.
!
ea = 0.0D+00
er = 0.0D+00
et = 0.0D+00
na = 0
nr = 0
nt = 0
sa = 0.0D+00
sr = 0.0D+00
st = 0.0D+00
!
! Loop over the particles.
!
do npart = 1, ntot
!
! Generate a new particle.
!
call source ( e, mu, azm, x, y, z )
do
!
! Compute the distance that the particle can travel through the slab,
! based on its current energy.
!
d = dist2c ( e )
!
! Update the particle's position by D units.
!
call update ( mu, azm, d, x, y, z )
!
! The particle was reflected by the shield, and this path is complete.
!
if ( x < 0.0D+00 ) then
nr = nr + 1
er = er + e
sr = sr + e * e
exit
!
! The particle was transmitted through the shield, and this path is complete.
!
else if ( thick < x ) then
nt = nt + 1
et = et + e
st = st + e * e
exit
!
! The particle collided with the shield, and was absorbed. This path is done.
!
else if ( absorb ( ) ) then
na = na + 1
ea = ea + e
sa = sa + e * e
exit
!
! The particle collided with the shield and was scattered.
! Find the scattering angle and energy, and continue along the new path.
!
else
call scatter ( e, mu, azm )
end if
end do
end do
!
! Print the results of the simulation.
!
call output ( na, ea, sa, nr, er, sr, nt, et, st, ntot )
end do
!
! Terminate.
!
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'REACTOR_SIMULATION:'
write ( *, '(a)' ) ' Normal end of execution.'
write ( *, '(a)' ) ' '
call timestamp ( )
stop
end
function absorb ( )
!*****************************************************************************80
!
!! ABSORB determines if a colliding particle is absorbed.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Output, logical ABSORB, is TRUE if the particle is absorbed.
!
! Local:
!
! real ( kind = rk ) PA, the probability of absorption.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
logical absorb
real ( kind = rk ), parameter :: pa = 0.1D+00
real ( kind = rk ) u
call random_number ( harvest = u )
if ( u <= pa ) then
absorb = .true.
else
absorb = .false.
end if
return
end
function cross ( e )
!*****************************************************************************80
!
!! CROSS returns the "cross section" of a particle based on its energy.
!
! Discussion:
!
! The particle's cross section is a measure of its likelihood to collide
! with the material of the slab. This quantity typically depends on both
! the particle's energy and the kind of medium through which it is traveling.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Input, real ( kind = rk ) E, the energy of the particle.
!
! Output, real ( kind = rk ) CROSS, the cross section.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) cross
real ( kind = rk ) e
real ( kind = rk ) s
real ( kind = rk ) y
s = abs ( sin ( 100.0D+00 * ( exp ( e ) - 1.0D+00 ) ) &
+ sin ( 18.81D+00 * ( exp ( e ) - 1.0D+00 ) ) )
y = max ( 0.02D+00, s )
cross = 10.0D+00 * exp ( -0.1D+00 / y )
return
end
function dist2c ( e )
!*****************************************************************************80
!
!! DIST2C returns the distance to collision.
!
! Discussion:
!
! Assuming the particle has a given energy, and assuming it is currently
! somewhere inside the shield, it is possible to determine a typical distance
! which the particle can travel before it collides with the material of
! the shield.
!
! The computation of the collision distance is made by estimating a
! "cross section" (as though having more energy made the particle "bigger"
! and hence more likely to collide) and then randomly selecting a distance
! that is logarithmically distributed.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Input, real ( kind = rk ) E, the energy of the particle.
!
! Output, real ( kind = rk ) DIST2C, the distance the particle can travel
! through the slab before colliding.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) cross
real ( kind = rk ) dist2c
real ( kind = rk ) e
real ( kind = rk ) u
call random_number ( harvest = u )
dist2c = - log ( u ) / cross ( e )
return
end
function energy ( )
!*****************************************************************************80
!
!! ENERGY assigns an energy to an emitted particle.
!
! Discussion:
!
! The energy E is in the range [EMIN,EMAX], with distribution
! const/sqrt(energy).
!
! An inverse function approach is used to compute this.
!
! The energies are measured in MeV.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Output, real ( kind = rk ) ENERGY, a randomly chosen energy that is
! distributed as described above.
!
! Local:
!
! real ( kind = rk ) EMIN, EMAX, the minimum and maximum
! energies.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) c
real ( kind = rk ), parameter :: emax = 2.5D+00
real ( kind = rk ), parameter :: emin = 1.0D-03
real ( kind = rk ) energy
real ( kind = rk ) u
call random_number ( harvest = u )
c = 1.0D+00 / ( 2.0D+00 * ( sqrt ( emax ) - sqrt ( emin ) ) )
energy = ( u / ( 2.0D+00 * c ) + sqrt ( emin ) )
energy = energy * energy
return
end
subroutine output ( na, ea, sa, nr, er, sr, nt, et, st, ntot )
!*****************************************************************************80
!
!! OUTPUT prints the results of the reactor shielding simulation.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Input, integer NA, number of particles absorbed by the slab.
!
! Input, real ( kind = rk ) EA, energy absorbed by the slab.
!
! Input, real ( kind = rk ) SA, the sum of the squares of the
! absorbed energies.
!
! Input, integer NR, number of particles reflected by the slab.
!
! Input, real ( kind = rk ) ER, energy reflected by the slab.
!
! Input, real ( kind = rk ) SR, the sum of the squares of the
! reflected energies.
!
! Input, integer NT, number of particles transmitted
! by the slab.
!
! Input, real ( kind = rk ) ET, energy transmitted through the slab.
!
! Input, real ( kind = rk ) ST, the sum of the squares of the
! transmitted energies.
!
! Input, integer NTOT, the total number of particles.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) ea
real ( kind = rk ) ea_ave
real ( kind = rk ) er
real ( kind = rk ) er_ave
real ( kind = rk ) et
real ( kind = rk ) et_ave
real ( kind = rk ) etot
integer na
integer nr
integer nt
integer ntot
real ( kind = rk ) pa
real ( kind = rk ) pr
real ( kind = rk ) pt
real ( kind = rk ) ptot
real ( kind = rk ) sa
real ( kind = rk ) sr
real ( kind = rk ) st
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' The Reactor Shielding Problem:'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) &
' Total Average'
write ( *, '(a)' ) ' # Energy ' &
// 'Percent Energy StDev'
write ( *, '(a)' ) ' '
etot = ea + er + et
if ( 0 < na ) then
ea_ave = ea / real ( na, kind = rk )
sa = sqrt ( sa / real ( na, kind = rk ) - ea_ave * ea_ave )
else
ea_ave = 0.0D+00
end if
pa = real ( na * 100, kind = rk ) / real ( ntot, kind = rk )
write ( *, '(a,2x,i8,2x,g14.6,2x,f6.2,2x,g14.6,2x,g14.6)' ) &
'Absorbed ', na, ea, pa, ea_ave, sa
if ( 0 < nr ) then
er_ave = er / real ( nr, kind = rk )
sr = sqrt ( sr / real ( nr, kind = rk ) - er_ave * er_ave )
else
er_ave = 0.0D+00
end if
pr = real ( nr * 100, kind = rk ) / real ( ntot, kind = rk )
write ( *, '(a,2x,i8,2x,g14.6,2x,f6.2,2x,g14.6,2x,g14.6)' ) &
'Reflected ', nr, er, pr, er_ave, sr
if ( 0 < nt ) then
et_ave = et / real ( nt, kind = rk )
st = sqrt ( st / real ( nt, kind = rk ) - et_ave * et_ave )
else
et_ave = 0.0D+00
end if
pt = real ( nt * 100, kind = rk ) / real ( ntot, kind = rk )
write ( *, '(a,2x,i8,2x,g14.6,2x,f6.2,2x,g14.6,2x,g14.6)' ) &
'Transmitted', nt, et, pt, et_ave, st
ptot = 100.0D+00
write ( *, '(a)' ) ' '
write ( *, '(a,2x,i8,2x,g14.6,2x,f6.2,2x,g14.6,2x,g14.6)' ) &
'Total ', ntot, etot, ptot
return
end
subroutine scatter ( e, mu, azm )
!*****************************************************************************80
!
!! SCATTER returns the new direction and energy of a particle that is scattered.
!
! Discussion:
!
! The scattering direction is chosen uniformly on the sphere.
!
! The energy of the scattered particle is chosen uniformly in
! [ 0.3*E, E ].
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Input/output, real ( kind = rk ) E. On input, the particle energy
! before collision. On output, the particle energy after collision
! and scattering.
!
! Output, real ( kind = rk ) MU, the cosine of the angle between the
! particle's direction and the X axis.
!
! Output, real ( kind = rk ) AZM, the azimuthal angle of the particle's
! direction.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) azm
real ( kind = rk ) e
real ( kind = rk ) mu
real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
real ( kind = rk ) u
call random_number ( harvest = u )
mu = - 1.0D+00 + 2.0D+00 * u
call random_number ( harvest = u )
azm = u * 2.0D+00 * pi
call random_number ( harvest = u )
e = ( u * 0.7D+00 + 0.3D+00 ) * e
return
end
subroutine source ( e, mu, azm, x, y, z )
!*****************************************************************************80
!
!! SOURCE generates a new particle from the neutron source.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Output, real ( kind = rk ) E, the initial energy of the particle.
!
! Output, real ( kind = rk ) MU, the cosine of the angle between the
! particle's direction and the X axis.
!
! Output, real ( kind = rk ) AZM, the azimuthal angle of the particle's
! direction.
!
! Output, real ( kind = rk ) X, Y, Z, the initial coordinates of the particle.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) azm
real ( kind = rk ) e
real ( kind = rk ) energy
real ( kind = rk ) mu
real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
real ( kind = rk ) u
real ( kind = rk ) x
real ( kind = rk ) y
real ( kind = rk ) z
call random_number ( harvest = u )
mu = u
call random_number ( harvest = u )
azm = u * 2.0D+00 * pi
x = 0.0D+00
y = 0.0D+00
z = 0.0D+00
e = energy ( )
return
end
subroutine timestamp ( )
!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
! Example:
!
! 31 May 2001 9:45:54.872 AM
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 18 May 2013
!
! Author:
!
! John Burkardt
!
implicit none
character ( len = 8 ) ampm
integer d
integer h
integer m
integer mm
character ( len = 9 ), parameter, dimension(12) :: month = (/ &
'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', &
'September', 'October ', 'November ', 'December ' /)
integer n
integer s
integer values(8)
integer y
call date_and_time ( values = values )
y = values(1)
m = values(2)
d = values(3)
h = values(5)
n = values(6)
s = values(7)
mm = values(8)
if ( h < 12 ) then
ampm = 'AM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Noon'
else
ampm = 'PM'
end if
else
h = h - 12
if ( h < 12 ) then
ampm = 'PM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Midnight'
else
ampm = 'AM'
end if
end if
end if
write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )
return
end
subroutine update ( mu, azm, d, x, y, z )
!*****************************************************************************80
!
!! UPDATE determines the position of the particle after it has traveled D units.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 April 2008
!
! Author:
!
! Original FORTRAN77 version by Kahaner, Moler, Nash.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters:
!
! Input, real ( kind = rk ) MU, the cosine of the angle between the
! particle's direction and the X axis.
!
! Input, real ( kind = rk ) AZM, the azimuthal angle of the particle's
! direction.
!
! Input, real ( kind = rk ) D, the distance the particle traveled.
!
! Input/output, real ( kind = rk ) X, Y, Z. On input, the previous
! coordinates of the particle. On output, the updated coordinates of the
! particle.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) azm
real ( kind = rk ) d
real ( kind = rk ) mu
real ( kind = rk ) s
real ( kind = rk ) x
real ( kind = rk ) y
real ( kind = rk ) z
s = sqrt ( 1.0D+00 - mu * mu )
x = x + d * mu
y = y + d * s * cos ( azm )
z = z + d * s * sin ( azm )
return
end