function reactor_simulation_test ( )

%*****************************************************************************80
%
%% reactor_simulation_test() tests reactor_simulation().
%
%  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:
%
%    18 November 2022
%
%  Author:
%
%    Original FORTRAN77 version by Kahaner, Moler, Nash.
%    This 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:
%
%    real AZM, the azimuthal angle of the particle's direction.
%
%    real D, the distance that the particle can
%    travel through the slab, given its current energy.
%
%    real E, the energy of the particle.
%
%    real EA, energy absorbed by the slab.
%
%    real ER, energy reflected by the slab.
%
%    real ET, energy transmitted through the slab.
%
%    real MU, the cosine of the angle between the
%    particle's direction and the X axis.
%
%    integer NA, number of particles absorbed by the slab.
%
%    integer NPART, the index of the current particle.
%
%    integer NR, number of particles reflected by the slab.
%
%    integer NT, number of particles transmitted by the slab.
%
%    integer NTOT, the total number of particles to be
%    emitted from the neutron source.
%
%    real SA, standard deviation of absorbed energy.
%
%    real SR, standard deviation of reflected energy.
%
%    real ST, standard deviation of transmitted energy.
%
%    real THICK, the thickness of the slab that is
%    intended to absorb most of the particles.
%
%    real X, Y, Z, the current position of the particle.
%
  addpath ( '../reactor_simulation' );

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'reactor_simulation_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test reactor_simulation()\n' );

  ntot = 100000;
  test_num = 5;
  thick = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The reactor shielding simulation.\n' );
  fprintf ( 1, '  Shield thickness is THICK = %g\n', thick );
  fprintf ( 1, '  Number of simulated particles is NTOT = %d\n', ntot );
  fprintf ( 1, '  Number of tests TEST_NUM = %d\n', test_num );

  for test = 1 : test_num

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Test # %d\n', test );
%
%  Initialize.
%
    ea = 0.0;
    er = 0.0;
    et = 0.0;
    na = 0;
    nr = 0;
    nt = 0;
    sa = 0.0;
    sr = 0.0;
    st = 0.0;
%
%  Main loop over NTOT particles.
%
    for npart = 1 : ntot
%
%  Generate a new particle.
%
      [ e, mu, azm, x, y, z ] = neutron_source ( );

      while ( true )
%
%  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.
%
        [ x, y, z ] = update ( mu, azm, d, x, y, z );
%
%  The particle was reflected by the shield, and this path is complete.
%
        if ( x < 0.0 )

          nr = nr + 1;
          er = er + e;
          sr = sr + e * e;
          break
%
%  The particle was transmitted through the shield, and this path is complete.
%
        elseif ( thick < x )

          nt = nt + 1;
          et = et + e;
          st = st + e * e;
          break
%
%  The particle collided with the shield, and was absorbed.  This path is done.
%
        else

          absorbed = absorb ( );

          if ( absorbed )

            na = na + 1;
            ea = ea + e;
            sa = sa + e * e;
            break
%
%  The particle collided with the shield and was scattered.
%  Find the scattering angle and energy, and continue along the new path.
%
          else

            [ mu, azm, e ] = scattering ( e );

          end

        end

      end

    end
%
%  Print the results of the simulation.
%
    output ( na, ea, sa, nr, er, sr, nt, et, st, ntot );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'reactor_simulation_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../reactor_simulation' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

