function continuity_exact ( )

%*****************************************************************************80
%
%% continuity_exact() computes a spiral velocity vector field.
%
%  Discussion:
%
%    Note that the continuous velocity field (U,V)(X,Y) that is discretely
%    sampled here satisfies the homogeneous continuity equation, that is,
%    it has zero divergence.  In other words:
%
%      dU/dX + dV/dY = 0.
%
%    This is by construction, since we have
%
%      U(X,Y) =  10 * d/dY ( PHI(X) * PHI(Y) )
%      V(X,Y) = -10 * d/dX ( PHI(X) * PHI(Y) )
%
%    which guarantees zero divergence.
%
%    The underlying function PHI is defined by
%
%      PHI(Z) = ( 1 - cos ( C * pi * Z ) ) * ( 1 - Z )^2
%
%    where C is a parameter.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 December 2013
%
%  Author:
%
%    John Burkardt
%
  output_file = 'continuity_exact.txt';

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'continuity_exact\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Write data describing a "spiral" vector flow\n' );
  fprintf ( 1, '  which satisfies the continuity equation.\n' );
  fprintf ( 1, '  (Actually, the flow does not spiral,\n' );
  fprintf ( 1, '  it''s simply concentric loops.)\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The flow is defined on an N+1 x N+1 grid.\n' );
  fprintf ( 1, '  A typical value might be 10 or 20.\n' );
  n = input ( '  Enter your desired value of N:  ' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The flow depends on a parameter C\n' );
  fprintf ( 1, '  A typical value might be 0.75 or 1.0\n' );
  c = input ( '  Enter your desired value of C:  ' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generating data on N+1 x N+1 grid, N = %d\n', n );
  fprintf ( 1, '  Using parameter C = %g\n', c );

  output_unit = fopen ( output_file, 'wt' );

  for i = 0 : n

    x = i / n;

    for j = 0 : n

      y = j / n;

      u =   10.0 * ( 1.0 - cos ( c * pi * x ) ) ...
               * ( 1.0 - x )^2 ...
               * ( ...
                   c * pi * sin ( c * pi * y ) * ( 1.0 - y )^2 ...
                 - ( 1.0 - cos ( c * pi * y ) ) ...
                   * 2.0 * ( 1.0 - y ) ...
                 );

      v = - 10.0 * ( 1.0 - cos ( c * pi * y ) ) ...
               * ( 1.0 - y )^2 ...
               * ( ...
                   c * pi * sin ( c * pi * x ) * ( 1.0 - x )^2 ...
                 - ( 1.0 - cos ( c * pi * x ) ) ...
                 * 2.0 * ( 1.0 - x ) ...
                 );

      fprintf ( output_unit, '  %10.4f  %10.4f  %10.4f  %10.4f\n', x, y, u, v );

    end
  end

  fclose ( output_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XYUV data written to "%s".\n', output_file );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This data can be plotted by GNUPLOT.\n' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'continuity_exact\n' );
  fprintf ( 1, '  Normal end of execution\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  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
