function [ h_array, uh_array, x, t ] = shallow_water_1d ( nx, nt, x_max, ...
  t_max, g )

%*****************************************************************************80
%
%% shallow_water_1d() approximates the 1D shallow water equations.
%
%  Discussion:
%
%    This code can be considered a 1D version of Cleve Moler's shallow
%    water equation solver.
%
%    The version of the shallow water equations being solved here is in
%    conservative form, and omits the Coriolis force.  The state variables
%    are H (the height) and UH (the mass velocity).
%
%    The equations have the form
%
%      dH/dt + d UH/dx = 0
%
%      d UH/dt + d ( U^2 H + 1/2 g H^2 )/dx = 0
%
%    Here U is the ordinary velocity, U = UH/H, and g is the gravitational
%    acceleration.
%
%    The initial conditions are used to specify ( H, UH ) at an equally
%    spaced set of points, and then the Lax-Wendroff method is used to advance
%    the solution through a number of equally spaced points in time, with 
%    boundary conditions supplying the first and last spatial values.
%
%
%    Some input values will result in an unstable calculation that
%    quickly blows up.  This is related to the Courant-Friedrichs-Lewy
%    condition, which requires that DT be small enough, relative to DX and
%    the velocity, that information cannot cross an entire cell.
%
%    A "reasonable" set of input quantities is
%
%      shallow_water_1d ( 41, 101, 1.0, 0.2, 9.8 )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    "The Shallow Water Equations",
%    Experiments with MATLAB.
%
%  Input:
%
%    integer NX, the number of spatial nodes, including two boundary nodes.
%
%    integer NT, the number of times steps, including the starting time.
%
%    real x_max, the length of the region.
%
%    real t_max, the time extent.
%
%    real G, the gravity constant.  G = 9.8 meters per second^2.
%
%  Output:
%
%    real H_ARRAY(NX,NT), the height for all space and time points.
%
%    real UH_ARRAY(NX,NT), the mass velocity for all space and time points.
%
%    real X(NX), the X coordinates.
%
%    real T(NT), the T coordinates.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'shallow_water_1d:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Solve the shallow water equations in 1D.\n' );

  if ( nargin < 5 )
    fprintf ( 1, '\n' );
  end
  if ( nargin < 1 )
    nx = 41;
    fprintf ( 1, '  Missing input NX defaults to %d\n', nx );
  end
  if ( nargin < 2 )
    nt = 101;
    fprintf ( 1, '  Missing input NT defaults to %d\n', nt );
  end
  if ( nargin < 3 )
    x_max = 1.0;
    fprintf ( 1, '  Missing input x_max defaults to %g\n', x_max );
  end
  if ( nargin < 4 )
    t_max = 0.2;
    fprintf ( 1, '  Missing input t_max defaults to %g\n', t_max );
  end
  if ( nargin < 5 )
    g = 9.8;
    fprintf ( 1, '  Missing input G defaults to %g\n', g );
  end
%
%  Set up arrays to collect all data.
%
  h_array = zeros ( nx, nt );
  uh_array = zeros ( nx, nt );

  hm = zeros ( nx, 1 );
  uhm = zeros ( nx, 1 );
%
%  Define the locations of the nodes and time steps and the spacing.
%
  x_min = 0.0;
  x = ( linspace ( x_min, x_max, nx ) )';
  dx = x_max / ( nx - 1 );

  t_min = 0.0;
  t = ( linspace ( t_min, t_max, nt ) )';
  dt = t_max / ( nt - 1 );
%
%  Time loop.
%
  for it = 1 : nt
%
%  Initial condition.
%
    if ( it == 1 )

      [ h, uh ] = initial_conditions ( x );
      h_min = 0.0;
      h_max = max ( h );
      uh_min = min ( uh );
      uh_max = max ( uh );

    else

%     for i = 1 : 10
%       fprintf ( 1, '%d %g %g %g\n', i, x(i), h(i), uh(i) );
%     end
%     fprintf ( 1, '\n' );
%
%  Take a half time step, estimating H and UH at the NX-1 spatial midpoints.
%
      hm(1:nx-1) = ( h(1:nx-1) + h(2:nx) ) / 2.0 ...
        - ( dt / 2.0 ) * ( uh(2:nx) - uh(1:nx-1) ) / dx;

      uhm(1:nx-1) = ( uh(1:nx-1) + uh(2:nx) ) / 2.0 ...
        - ( dt / 2.0 ) * ( ...
          uh(2:nx).^2    ./ h(2:nx)   + 0.5 * g * h(2:nx).^2 ...
        - uh(1:nx-1).^2  ./ h(1:nx-1) - 0.5 * g * h(1:nx-1).^2 ) / dx;
%
%  Take a full time step, evaluating the derivative at the half time step,
%  to estimate the solution at the NX-2 interior nodes.
%
      h(2:nx-1) = h(2:nx-1) ...
        - dt * ( uhm(2:nx-1) - uhm(1:nx-2) ) / dx;

      uh(2:nx-1) = uh(2:nx-1) ...
        - dt * ( ...
          uhm(2:nx-1).^2  ./ hm(2:nx-1) + 0.5 * g * hm(2:nx-1).^2 ...
        - uhm(1:nx-2).^2  ./ hm(1:nx-2) - 0.5 * g * hm(1:nx-2).^2 ) / dx;

%     for i = 1 : 10
%       fprintf ( 1, '%d %g %g %g\n', i, x(i), hm(i), uhm(i) );
%     end
%     fprintf ( 1, '\n' );

%     for i = 1 : 10
%       fprintf ( 1, '%d %g %g %g\n', i, x(i), h(i), uh(i) );
%     end
%     return

    end
%
%  Apply the boundary conditions.
%
    [ h, uh ] = boundary_conditions ( nx, nt, h, uh, t(it) );
%
%  Copy data into the big arrays.
%
    h_array(1:nx,it) = h(1:nx);
    uh_array(1:nx,it) = uh(1:nx);

    figure ( 1 );
    clf ( );

    subplot ( 2, 1, 1 );
    axis ( [ x_min, x_max, h_min, h_max ] );
    area ( x, h_array (:,it) );
    title_string = sprintf ( 'H(T), Step %3d, Time = %f', it, t(it) );
    title ( title_string );
    xlabel ( 'X' );
    ylabel ( 'H(X,T)' );
    grid ( 'on' );

    subplot ( 2, 1, 2 );
    area ( x, uh_array(:,it), 0.0, 'FaceColor', [ 0.6, 0.1, 0.1 ] );
    title_string = sprintf ( 'UH(T), Step %3d, Time = %f', it, t(it) );
    title ( title_string );
    xlabel ( 'X' );
    ylabel ( 'UH(X,T)' );
    grid ( 'on' );
    
  end
%
%  Write data for possible animation.
%
  if ( false )
    filename = 'shallow_water_1d.mat';
    save ( filename, 'x', 't', 'h_array', 'uh_array' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  data saved in file "%s"\n', filename );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'shallow_water_1d:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function [ h, uh ] = boundary_conditions ( nx, nt, h, uh, t )

%*****************************************************************************80
%
%% boundary_conditions() sets the boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NX, the number of spatial nodes.
%
%    integer NT, the number of times steps.
%
%    real H(NX), the height for all space.
%
%    real UH(NX), the mass velocity for all space.
%
%    real T, the current time.
%
%  Output:
%
%    real H(NX), the height, with H(1) and H(NX) adjusted for
%    boundary conditions.
%
%    real UH(NX), the mass velocity, with UH(1) and UH(NX)
%    adjusted for boundary conditions.
%
  bc = 1;
%
%  Periodic boundary conditions on H and UH.
%
  if ( bc == 1 )
    h(1) = h(nx-1);
    h(nx) = h(2);
    uh(1) = uh(nx-1);
    uh(nx) = uh(2);
%
%  Free boundary conditions on H and UH.
%
  elseif ( bc == 2 )
    h(1) = h(2);
    h(nx) = h(nx-1);
    uh(1) = uh(2);
    uh(nx) = uh(nx-1);
%
%  Reflective boundary conditions on UH, free boundary conditions on H.
%
  elseif ( bc == 3 )
    h(1) = h(2);
    h(nx) = h(nx-1);
    uh(1) = - uh(2);
    uh(nx) = - uh(nx-1);
  end

  return
end
function [ h, uh ] = initial_conditions ( x )

%*****************************************************************************80
%
%% initial_conditions() sets the initial conditions.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 June 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X(NX,1), the coordinates of the nodes.
%
%  Output:
%
%    real H(NX,1), the initial height for all space.
%
%    real UH(NX,1), the initial mass velocity for all space.
%
  h = 2.0 + sin ( 2.0 * pi * x );
  uh = zeros ( size ( x ) );

  return
end

