# include <cmath>
# include <cstdlib>
# include <ctime>
# include <fstream>
# include <iomanip>
# include <iostream>
# include <string>

using namespace std;

int main ( int argc, char *argv[] );
double cpu_time ( );
void gnuplot_fxy ( int m, int n, double *x, double *y, double *z, string header );
double *r8vec_linspace_new ( int n, double a_first, double a_last );
void timestamp ( );

//****************************************************************************80

int main ( int argc, char *argv[] )

//****************************************************************************80
//
//  Purpose:
//
//    heated_plate() solves the steady state heat equation on a rectangular region.
//
//  Discussion:
//
//    The sequential version of this program needs approximately
//    18/epsilon iterations to complete.
//
//
//    The physical region, and the boundary conditions, are suggested
//    by this diagram;
//
//                   W = 0
//             +------------------+
//             |                  |
//    W = 100  |                  | W = 100
//             |                  |
//             +------------------+
//                   W = 100
//
//    The region is covered with a grid of M by N nodes, and an N by N
//    array W is used to record the temperature.  The correspondence between
//    array indices and locations in the region is suggested by giving the
//    indices of the four corners:
//
//                  I = 0
//          [0][0]-------------[0][N-1]
//             |                  |
//      J = 0  |                  |  J = N-1
//             |                  |
//        [M-1][0]-----------[M-1][N-1]
//                  I = M-1
//
//    The steady state solution to the discrete heat equation satisfies the
//    following condition at an interior grid point:
//
//      W[Central] = (1/4) * ( W[North] + W[South] + W[East] + W[West] )
//
//    where "Central" is the index of the grid point, "North" is the index
//    of its immediate neighbor to the "north", and so on.
//
//    Given an approximate solution of the steady state heat equation, a
//    "better" solution is given by replacing each interior point by the
//    average of its 4 neighbors - in other words, by using the condition
//    as an ASSIGNMENT statement:
//
//      W[Central]  <=  (1/4) * ( W[North] + W[South] + W[East] + W[West] )
//
//    If this process is repeated often enough, the difference between successive
//    estimates of the solution will go to zero.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 April 2025
//
//  Author:
//
//    Original C version by Michael Quinn.
//    This version by John Burkardt.
//
//  Reference:
//
//    Michael Quinn,
//    Parallel Programming in C with MPI and OpenMP,
//    McGraw-Hill, 2004,
//    ISBN13: 978-0071232654,
//    LC: QA76.73.C15.Q55.
//
//  Local:
//
//    double DIFF, the norm of the change in the solution from one iteration
//    to the next.
//
//    double MEAN, the average of the boundary values, used to initialize
//    the values of the solution in the interior.
//
//    double U[M][N], the solution at the previous iteration.
//
//    double W[M][N], the solution computed at the latest iteration.
//
{
# define M 100
# define N 500

  double ctime;
  double ctime1;
  double ctime2;
  double diff;
  double epsilon;
  int i;
  int iterations;
  int iterations_print;
  int j;
  int m;
  double mean;
  int n;
  double u[M*N];
  double w[M*N];
  double *x;
  double *y;

  timestamp ( );
  cout << "\n";
  cout << "heated_plate():\n";
  cout << "  C++ version\n";
  cout << "  A program to solve for the steady state temperature distribution\n";
  cout << "  over a rectangular plate.\n";

  m = M;
  n = N;
  cout << "\n";
  cout << "  Spatial grid of " << M << " by " << N << " points.\n";

  epsilon = 0.0001;
  cout << "\n";
  cout << "  The iteration will be repeated until the change is <= "
       << epsilon << "\n";
  diff = epsilon;

  x = r8vec_linspace_new ( n, 0.0, 1.0 );
  y = r8vec_linspace_new ( m, 0.0, 1.0 );
//
//  Set the boundary values, which don't change.
//
  for ( i = 1; i < M - 1; i++ )
  {
    w[i+0*m] = 100.0;
  }
  for ( i = 1; i < M - 1; i++ )
  {
    w[i+(N-1)*m] = 100.0;
  }
  for ( j = 0; j < N; j++ )
  {
    w[M-1+j*m] = 100.0;
  }
  for ( j = 0; j < N; j++ )
  {
    w[0+j*m] = 0.0;
  }
//
//  Average the boundary values, to come up with a reasonable
//  initial value for the interior.
//
  mean = 50.0;
//
//  Initialize the interior solution to the mean value.
//
  for ( i = 1; i < M - 1; i++ )
  {
    for ( j = 1; j < N - 1; j++ )
    {
      w[i+j*m] = mean;
    }
  }
//
//  iterate until the  new solution W differs from the old solution U
//  by no more than EPSILON.
//
  iterations = 0;
  iterations_print = 1;
  cout << "\n";
  cout << " Iteration  Change\n";
  cout << "\n";
  ctime1 = cpu_time ( );

  while ( epsilon <= diff )
  {
//
//  Save the old solution in U.
//
    for ( i = 0; i < M; i++ )
    {
      for ( j = 0; j < N; j++ )
      {
        u[i+j*m] = w[i+j*m];
      }
    }
//
//  Determine the new estimate of the solution at the interior points.
//  The new solution W is the average of north, south, east and west neighbors.
//
    diff = 0.0;
    for ( i = 1; i < M - 1; i++ )
    {
      for ( j = 1; j < N - 1; j++ )
      {
        w[i+j*m] = ( u[i-1+j*m] + u[i+1+j*m] + u[i+(j-1)*m] + u[i+(j+1)*m] ) / 4.0;

        if ( diff < fabs ( w[i+j*m] - u[i+j*m] ) )
        {
          diff = fabs ( w[i+j*m] - u[i+j*m] );
        }
      }
    }
    iterations++;
    if ( iterations == iterations_print )
    {
      cout << "  " << setw(8) << iterations
           << "  " << diff << "\n";
      iterations_print = 2 * iterations_print;
    }
  }
  ctime2 = cpu_time ( );
  ctime = ctime2 - ctime1;

  cout << "\n";
  cout << "  " << setw(8) << iterations
       << "  " << diff << "\n";
  cout << "\n";
  cout << "  Error tolerance achieved.\n";
  cout << "  CPU time = " << ctime << "\n";
//
//  Write the solution to a file.
//  GNUPLOT requires a regular grid, and a blank line before a new 
//  Y value.
//
  gnuplot_fxy ( m, n, x, y, w, "heated_plate" );
//
//  Terminate.
//
  cout << "\n";
  cout << "heated_plate():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  delete [] x;
  delete [] y;

  return 0;

# undef M
# undef N
}
//****************************************************************************80

double cpu_time ( )

//****************************************************************************80
//
//  Purpose:
//
//    cpu_time() returns the current reading on the CPU clock.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    06 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double CPU_TIME, the current reading of the CPU clock, in seconds.
//
{
  double value;

  value = ( double ) clock ( ) / ( double ) CLOCKS_PER_SEC;

  return value;
}
//****************************************************************************80

void gnuplot_fxy ( int m, int n, double *x, double *y, double *z, 
  string header )

//****************************************************************************80
//
//  Purpose:
//
//    gnuplot_fxy() plots z=f(x,y).
//
//  Discussion:
//
//    Actually, we simply create two files for processing by gnuplot().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 April 2025
//
//  Author:
//
//    John Burkardt
//
{
  string command_filename;
  ofstream command_unit;
  string data_filename;
  ofstream data_unit;
  int i;
  int j;
  string png_filename;
//
//  Create a graphics data file.
//
  data_filename = header + "_data.txt";
  data_unit.open ( data_filename.c_str ( ) );
  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      data_unit << x[j] << "  " << y[i] << "  " << z[i+j*m] << "\n";
    }
    if ( i < m - 1 )
    {
      data_unit << "\n";
    }
  }
  data_unit.close ( );
  cout << "  Created graphics data file '" << data_filename << "'\n";

  png_filename = header + ".png";

  png_filename = header + ".png";
//
//  Create graphics command file.
//
  command_filename = header + "_commands.txt";
  command_unit.open ( command_filename.c_str ( ) );
  command_unit << "# " << command_filename << "\n";
  command_unit << "#\n";
  command_unit << "# Usage:\n";
  command_unit << "#  gnuplot < " << command_filename << "\n";
  command_unit << "#\n";
  command_unit << "set term png\n";
  command_unit << "set output '" << png_filename << "'\n";
  command_unit << "set xlabel '<--- X --->'\n";
  command_unit << "set ylabel '<--- Y --->'\n";
  command_unit << "set title '" << header << "'\n";
  command_unit << "set grid\n";
  command_unit << "unset surface\n";
  command_unit << "set contour base\n";
  command_unit << "set view map\n";
  command_unit << "set pm3d\n";
  command_unit << 
    "splot '" << data_filename << "' with pm3d\n";
  command_unit.close ( );
  cout << "  Created command file '" << command_filename << "'\n";

  return;
}
//****************************************************************************80

double *r8vec_linspace_new ( int n, double a_first, double a_last )

//****************************************************************************80
//
//  Purpose:
//
//    r8vec_linspace_new() creates a vector of linearly spaced values.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//    4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
//
//    In other words, the interval is divided into N-1 even subintervals,
//    and the endpoints of intervals are used as the points.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2011
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of entries in the vector.
//
//    double A_FIRST, A_LAST, the first and last entries.
//
//  Output:
//
//    double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data.
//
{
  double *a;
  int i;

  a = new double[n];

  if ( n == 1 )
  {
    a[0] = ( a_first + a_last ) / 2.0;
  }
  else
  {
    for ( i = 0; i < n; i++ )
    {
      a[i] = ( ( double ) ( n - 1 - i ) * a_first 
             + ( double ) (         i ) * a_last ) 
             / ( double ) ( n - 1     );
    }
  }
  return a;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    19 March 2018
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}