# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "ppc_array.h"
# include "ppc_fd1.h"

int main ( int argc, char **argv );
double bc_L ( double t );
double bc_R ( double t );
double exact_sol ( double x, double t );
double ic ( double x );
void timestamp ( );

/******************************************************************************/

int main ( int argc, char **argv )

/******************************************************************************/
/*
  Purpose:

    demo2() solves the heat equation using four different solvers.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

   int argc, char **argv: the usual commandline arguments.
*/
{
  char *endptr;

  timestamp ( );
  printf ( "\n" );
  printf ( "demo2():\n" );
  printf ( "  C version\n" );
  printf ( "  Solve heat equation with exact solution\n" );
  printf ( "  u(x,t) = erf ( 0.5 * x / sqrt ( t ) ) )\n" );

  if ( argc != 4 )
  {
    show_usage_and_exit ( argv[0] );
  }
  double T = strtod ( argv[1], &endptr );

  if ( *endptr != '\0' || T <= 0.0 )
  {
    show_usage_and_exit ( argv[0] );
  }

  int n = strtol ( argv[2], &endptr, 10 );
  if ( *endptr != '\0' || n < 1 )
  {
    show_usage_and_exit ( argv[0] );
  }

  int m = strtol ( argv[3], &endptr, 10 );
  if (* endptr != '\0' || m < 1 )
  {
    show_usage_and_exit ( argv[0] );
  }
/*
  Report input values.
*/
  printf ( "\n" );
  printf ( "  Final time T = %g\n", T );
  printf ( "  Number of spatial nodes n = %d\n", n );
  printf ( "  Number of time steps m = %d\n", m );
/*
  Define problem parameters.
*/
  struct heat_solve prob = 
  {
    .a = -1.0,
    .b = 1.0,
    .T = T,
    .n = n,
    .m = m,
    .ic = ic,
    .bcL = bc_L,
    .bcR = bc_R,
    .method = FD_undefined,
    .exact_sol = exact_sol,
    .maple_out = NULL,
    .matlab_out = NULL,
    .geomview_out = NULL,
  };

  make_matrix ( prob.u, m+1, n+1 );
/*
  Solve the problem through four different algorithms.
*/
  printf ( "\n" );
  prob.method = FD_explicit;
  prob.maple_out = "demo2_explicit.mpl";
  prob.matlab_out = "demo2_explicit.m";
  prob.geomview_out = "demo2_explicit.gv";
  heat_solve ( &prob );
  if ( prob.exact_sol != NULL )
  {
    printf ( "  Absolute error = %g\n", prob.error );
  }

  printf ( "\n" );
  prob.method = FD_implicit;
  prob.maple_out = "demo2_implicit.mpl";
  prob.matlab_out = "demo2_implicit.m";
  prob.geomview_out = "demo2_implicit.gv";
  heat_solve ( &prob );
  if ( prob.exact_sol != NULL )
  {
    printf ( "  Absolute error = %g\n", prob.error );
  }

  printf ( "\n" );
  prob.method = FD_crank_nicolson;
  prob.maple_out = "demo2_crank_nicolson.mpl";
  prob.matlab_out = "demo2_crank_nicolson.m";
  prob.geomview_out = "demo2_crank_nicolson.gv";
  heat_solve ( &prob );
  if ( prob.exact_sol != NULL )
  {
    printf ( "  Absolute error = %g\n", prob.error );
  }

  printf ( "\n" );
  prob.method = FD_seidman_sweep;
  prob.maple_out = "demo2_seidman_sweep.mpl";
  prob.matlab_out = "demo2_seidman_sweep.m";
  prob.geomview_out = "demo2_seidman_sweep.gv";
  heat_solve ( &prob );
  if ( prob.exact_sol != NULL )
  {
    printf ( "  Absolute error = %g\n", prob.error );
  }
/*
  Free memory.
*/
  free_matrix ( prob.u );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "demo2():\n" );
  printf ( "  Normal end of execution.\n" );
  timestamp ( );

  return EXIT_SUCCESS;
}
/******************************************************************************/

double bc_L ( double t )

/******************************************************************************/
/*
  Purpose:

    bc_L() evaluates the left hand boundary condition at any time.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double t: the time.

  Output:

    double bc_L: the solution value at (-1.0,t).
*/
{
  return exact_sol ( -1.0, t );
}
/******************************************************************************/

double bc_R ( double t )

/******************************************************************************/
/*
  Purpose:

    bc_R() evaluates the right hand boundary condition at any time.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double t: the time.

  Output:

    double bc_R: the solution value at (+1.0,t).
*/
{
  return exact_sol ( 1.0, t );
}
/******************************************************************************/

double exact_sol ( double x, double t )

/******************************************************************************/
/*
  Purpose:

    exact_sol() evaluates the exact solution at any place and time.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double x: the position.

    double t: the time.

  Output:

    double exact_sol: the solution value at (x,t).
*/
{
  double value;

  if ( t == 0.0 )
  {
    if ( x < 0.0 )
    {
      value = -1.0;
    }
    else if ( x == 0.0 )
    {
      value = 0.0;
    }
    else
    {
      value = 1.0;
    }
  }
  else
  {
    value = erf ( 0.5 * x / sqrt ( t ) );
  }
  return value;
}
/******************************************************************************/

double ic ( double x )

/******************************************************************************/
/*
  Purpose:

    ic() evaluates the initial condition at any place.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double x: the position.

  Output:

    double ic: the solution value at (x,0.0).
*/
{
  return exact_sol ( x, 0.0 );
}
/******************************************************************************/

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

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

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

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

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}