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

using namespace std;

int main ( );
void dtable_data_write ( ofstream &output, int m, int n, double table[] );
void dtable_header_write ( string output_filename, ofstream &output, int m, 
  int n );
void dtable_write ( string output_filename, int m, int n, double table[], 
  bool header );
void f ( double a, double b, double t0, double t, int n, double x[], 
  double value[] );
double r8_epsilon ( );
void timestamp ( );
char *timestring ( );
void u0 ( double a, double b, double t0, int n, double x[], double value[] );
double ua ( double a, double b, double t0, double t );
double ub ( double a, double b, double t0, double t );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for FD1D_HEAT_EXPLICIT.
//
//  Discussion:
//
//    FD1D_HEAT_EXPLICIT solves the 1D heat equation with an explicit method.
//
//    This program solves
//
//      dUdT - k * d2UdX2 = F(X,T)
//
//    over the interval [A,B] with boundary conditions
//
//      U(A,T) = UA(T),
//      U(B,T) = UB(T),
//
//    over the time interval [T0,T1] with initial conditions
//
//      U(X,T0) = U0(X)
//
//    The code uses the finite difference method to approximate the
//    second derivative in space, and an explicit forward Euler approximation
//    to the first derivative in time.
//
//    The finite difference form can be written as
//
//      U(X,T+dt) - U(X,T)                  ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) )
//      ------------------  = F(X,T) + k *  ------------------------------------
//               dt                                   dx * dx
//
//    or, assuming we have solved for all values of U at time T, we have
//
//      U(X,T+dt) = U(X,T) 
//                + dt * ( F(X,T) 
//                         + k * ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) / dx / dx )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    30 May 2009
//
//  Author:
//
//    John Burkardt
//
{
  double cfl;
  double *fvec;
  bool header;
  int i;
  int j;
  double k;
  double *t;
  double t_delt;
  string t_file;
  double t_max;
  double t_min;
  int t_num;
  double *u;
  string u_file;
  double *x;
  double x_delt;
  string x_file;
  double x_max;
  double x_min;
  int x_num;

  timestamp ( );
  cout << "\n";
  cout << "FD1D_HEAT_EXPLICIT\n";
  cout << "  C++ version\n";
  cout << "\n";
  cout << "  Finite difference solution of\\n";
  cout << "  the time dependent 1D heat equation\n";
  cout << "\n";
  cout << "    Ut - k * Uxx = F(x,t)\n";
  cout << "\n";
  cout << "  for space interval A <= X <= B with boundary conditions\n";
  cout << "\n";
  cout << "    U(A,t) = UA(t)\n";
  cout << "    U(B,t) = UB(t)\n";
  cout << "\n";
  cout << "  and time interval T0 <= T <= T1 with initial condition\n";
  cout << "\n";
  cout << "    U(X,T0) = U0(X).\n";
  cout << "\n";
  cout << "  A second order difference approximation is used for Uxx.\n";
  cout << "\n";
  cout << "  A first order forward Euler difference approximation\n";
  cout << "  is used for Ut.\n";

  k = 5.0E-07;
//
//  Set X values.
//  Increasing X_NUM to 16 while keeping all other data constant will allow
//  instabilities to begin to show up.
//
  x_min = 0.0;
  x_max = 0.3;
  x_num = 11;
  x_delt = ( x_max - x_min ) / ( double ) ( x_num - 1 );

  x = new double[x_num];

  for ( i = 0; i < x_num; i++ )
  {
    x[i] = ( ( double ) ( x_num - i - 1 ) * x_min   
           + ( double ) (         i     ) * x_max ) 
           / ( double ) ( x_num     - 1 );
  }
// 
//  Set T values.
//
  t_min = 0.0;
  t_max = 22000.0;
  t_num = 51;
  t_delt = ( t_max - t_min ) / ( double ) ( t_num - 1 );

  t = new double[t_num];

  for ( j = 0; j < t_num; j++ )
  {
    t[j] = ( ( double ) ( t_num - j - 1 ) * t_min   
           + ( double ) (         j     ) * t_max ) 
           / ( double ) ( t_num     - 1 );
  }
//
//  Check the CFL condition.
//
  cfl = k * t_delt / x_delt / x_delt;

  cout << "\n";
  cout << "  CFL stability criterion value = " << cfl << "\n";

  if ( 0.5 <= cfl )
  {
    cout << "\n";
    cout << "FD1D_HEAT_EXPLICIT - Fatal error!\n";
    cout << "  CFL condition failed.\n";
    cout << "  0.5 <= K * dT / dX / dX = " << cfl << "\n";
    exit ( 1 );
  }
//
//  Set the initial data, for time T_MIN.
//
  u = new double[x_num*t_num];

  u0 ( x_min, x_max, t_min, x_num, x, u );
//
//  Fill in the U array for time steps after T_MIN, up to T_MAX.
//
  fvec = new double[x_num];

  for ( j = 1; j < t_num; j++ )
  {
    u[0+j*x_num] = ua ( x_min, x_max, t_min, t[j-1] );

    f ( x_min, x_max, t_min, t[j-1], x_num, x, fvec );

    for ( i = 1; i < x_num - 1; i++ )
    {
      u[i+j*x_num] = u[i+(j-1)*x_num] + t_delt * ( fvec[i]
                     + k * (      u[i-1+(j-1)*x_num] 
                            - 2 * u[i  +(j-1)*x_num] 
                            +     u[i+1+(j-1)*x_num] ) / x_delt / x_delt 
                     );
    }

    u[x_num-1+j*x_num] = ub ( x_min, x_max, t_min, t[j-1] );

  }

  x_file = "x.txt";
  header = false;
  dtable_write ( x_file, 1, x_num, x, header );

  cout << "\n";
  cout << "  X data written to \"" << x_file << "\".\n";

  t_file = "t.txt";
  header = false;
  dtable_write ( t_file, 1, t_num, t, header );

  cout << "  T data written to \"" << t_file << "\".\n";

  u_file = "u.txt";
  header = false;
  dtable_write ( u_file, x_num, t_num, u, header );

  cout << "  U data written to \"" << u_file << "\".\n";

  cout << "\n";
  cout << "FD1D_HEAT_EXPLICIT\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  delete [] fvec;
  delete [] t;
  delete [] u;
  delete [] x;

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

void dtable_data_write ( ofstream &output, int m, int n, double table[] )

//****************************************************************************80
//
//  Purpose:
//
//    DTABLE_DATA_WRITE writes data to a DTABLE file.
//
//  Discussion:
//
//    The file should already be open.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 December 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, a pointer to the output stream.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double TABLE[M*N], the table data.
//
{
  int i;
  int j;

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      output << setw(10) << table[i+j*m] << "  ";
    }
    output << "\n";
  }

  return;
}
//****************************************************************************80
 
void dtable_header_write ( string output_filename, ofstream &output, int m, 
  int n )
 
//****************************************************************************80
//
//  Purpose:
//
//    DTABLE_HEADER_WRITE writes the header of a DTABLE file.
//
//  Discussion:
//
//    The file should already be open.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string OUTPUT_FILENAME, the output filename.
//
//    Input, ofstream &OUTPUT, the output stream.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
{
  char *s;

  s = timestring ( );

  output << "#  " << output_filename << "\n";
  output << "#  created by dtable_write in table_io.C" << "\n";
  output << "#  at " << s << "\n";
  output << "#\n";
  output << "#  Spatial dimension M = " << m << "\n";
  output << "#  Number of points N =  " << n << "\n";
  output << "#  EPSILON (unit roundoff) = " << r8_epsilon ( ) << "\n";
  output << "#\n";

  delete [] s;

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

void dtable_write ( string output_filename, int m, int n, double table[], 
  bool header )

//****************************************************************************80
//
//  Purpose:
//
//    DTABLE_WRITE writes information to a DTABLE file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string OUTPUT_FILENAME, the output filename.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double TABLE[M*N], the table data.
//
//    Input, bool HEADER, is TRUE if the header is to be included.
//
{
  ofstream output;

  output.open ( output_filename.c_str ( ) );

  if ( !output )
  {
    cerr << "\n";
    cerr << "DTABLE_WRITE - Fatal error!\n";
    cerr << "  Could not open the output file.\n";
    return;
  }

  if ( header )
  {
    dtable_header_write ( output_filename, output, m, n );
  }

  dtable_data_write ( output, m, n, table );

  output.close ( );

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

void f ( double a, double b, double t0, double t, int n, double x[], 
  double value[] )

//****************************************************************************80
//
//  Purpose:
//
//    F returns the right hand side of the heat equation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    15 May 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double A, B, the left and right endpoints.
//
//    Input, double T0, the initial time.
//
//    Input, double T, the current time.
//
//    Input, int N, the number of points.
//
//    Input, double X[N], the current spatial positions.
//
//    Output, double VALUE[N], the prescribed value of U(X(:),T0).
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    value[i] = 0.0;
  }
  return;
}
//****************************************************************************80

double r8_epsilon ( )

//****************************************************************************80
//
//  Purpose:
//
//    R8_EPSILON returns the R8 roundoff unit.
//
//  Discussion:
//
//    The roundoff unit is a number R which is a power of 2 with the 
//    property that, to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but 
//      1 = ( 1 + R / 2 )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    01 July 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_EPSILON, the double precision round-off unit.
//
{
  double r;

  r = 1.0;

  while ( 1.0 < ( double ) ( 1.0 + r )  )
  {
    r = r / 2.0;
  }

  return ( 2.0 * r );
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    02 October 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}
//****************************************************************************80

char *timestring ( )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTRING returns the current YMDHMS date as a string.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *TIMESTRING, a string containing the current YMDHMS date.
//
{
# define TIME_SIZE 40

  const struct tm *tm;
  size_t len;
  time_t now;
  char *s;

  now = time ( NULL );
  tm = localtime ( &now );

  s = new char[TIME_SIZE];

  len = strftime ( s, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  return s;
# undef TIME_SIZE
}
//****************************************************************************80

void u0 ( double a, double b, double t0, int n, double x[], double value[] )

//****************************************************************************80
//
//  Purpose:
//
//    U0 returns the initial condition at the starting time.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    15 May 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double A, B, the left and right endpoints
//
//    Input, double T0, the initial time.
//
//    Input, double T, the current time.
//
//    Input, int N, the number of points where initial data is needed.
//
//    Input, double X[N], the positions where initial data is needed.
//
//    Output, double VALUE[N], the prescribed value of U(X,T0).
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    value[i] = 100.0;
  }
  return;
}
//****************************************************************************80

double ua ( double a, double b, double t0, double t )

//****************************************************************************80
//
//  Purpose:
//
//    UA returns the Dirichlet boundary condition at the left endpoint.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    15 May 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double A, B, the left and right endpoints
//
//    Input, double T0, the initial time.
//
//    Input, double T, the current time.
//
//    Output, double UA, the prescribed value of U(A,T).
//
{
  double value;

  value = 20.0;

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

double ub ( double a, double b, double t0, double t )

//****************************************************************************80
//
//  Purpose:
//
//    UB returns the Dirichlet boundary condition at the right endpoint.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    15 May 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double A, B, the left and right endpoints
//
//    Input, double T0, the initial time.
//
//    Input, double T, the current time.
//
//    Output, double UB, the prescribed value of U(B,T).
//
{
  double value;

  value = 20.0;

  return value;
}