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

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

struct plot3d 
{
  double a;             // a < x < b
  double b;             // a < x < b
  double T;             // 0 < t < T
  int n;                // number of x subintervals
  int m;                // number of t subintervals
  double **u;           // (m+1)x(n+1) data array
  char *maple_out;      // output file for maple graphics 
  char *matlab_out;     // output file for matlab graphics 
  char *geomview_out;   // output file for geomview graphics 
};

double error_vs_exact ( struct heat_solve *prob );
void heat_solve ( struct heat_solve *prob );
void heat_solve_crank_nicolson ( struct heat_solve *prob );
void heat_solve_explicit ( struct heat_solve *prob );
void heat_solve_implicit ( struct heat_solve *prob );
void heat_solve_seidman_sweep ( struct heat_solve *prob );
void plot_with_geomview ( struct plot3d *data );
void plot_with_maple ( struct plot3d *data );
void plot_with_matlab ( struct plot3d *data );
void plot3d ( struct plot3d *data );
void show_usage_and_exit ( char *progname );
void trisolve ( int n, const double *a, double *d, const double *c, double *b, 
  double *x );
void write_plotting_script ( struct heat_solve *prob );

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

double error_vs_exact ( struct heat_solve *prob )

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

    error_vs_exact() compares computed solution to formula for exact solution.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 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:

    struct heat_solve *prob: data defining the problem.

  Output:

    double diff_max: the maximum value of the absolute value of the
    difference between the computed and exact solutions at grid points.
*/
{
  double diff_max;
  int i;
  int j;
  double t;
  double u_comp;
  double u_exact;
  double x;

  diff_max = 0.0;

  for ( i = 0; i <= prob->m; i++ )
  {
    t = i * prob->T / prob->m;
    for ( j = 0; j <= prob->n; j++ )
    {
      x = ( ( prob->n - j ) * prob->a + j * prob->b ) / prob->n;
      u_exact = prob->exact_sol ( x, t );
      u_comp = prob->u[i][j];
      diff_max = fmax ( diff_max, fabs ( u_exact - u_comp ) );
    }
  }

  return diff_max;
}
/******************************************************************************/

void heat_solve ( struct heat_solve *prob )

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

    heat_solve() chooses a method to solve the heat equation.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 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:

    struct heat_solve *prob: data defining the problem.

  Output:

    struct heat_solve *prob: updated data which now includes a solution.
*/
{
  switch ( prob->method )
  {
    case FD_explicit:
      heat_solve_explicit ( prob );
      break;
    case FD_implicit:
      heat_solve_implicit ( prob );
      break;
    case FD_crank_nicolson:
      heat_solve_crank_nicolson ( prob );
      break;
    case FD_seidman_sweep:
      heat_solve_seidman_sweep ( prob );
      break;
    default:
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "heat_solve():\n" );
      fprintf ( stderr, "  Missing method specification\n" );
      exit ( EXIT_FAILURE );
  }
  write_plotting_script ( prob );

  if ( prob->exact_sol != NULL )
  {
    prob->error = error_vs_exact ( prob );
  }

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

void heat_solve_crank_nicolson ( struct heat_solve *prob )

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

    heat_solve_crank_nicolson() uses Crank-Nicolson to solve the heat equation.

  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:

    struct heat_solve *prob: data defining the problem.

  Output:

    struct heat_solve *prob: updated data which now includes a solution.
*/
{
  int m = prob->m;
  int n = prob->n;

  double *b;
  double *c;
  double *d;
  double dt = prob->T / m;
  double dx = ( prob->b - prob->a ) / n;
  double r = dt / ( dx * dx );
  double **u = prob->u;

  printf ( "  Crank Nicolson linear equation solver:\n");
  printf ( "  r = %g\n", r );

  make_vector ( d, n-1 );
  make_vector ( c, n-2 );
  make_vector ( b, n+1 );
/*
  Set the lower and upper diagonals.
*/
  for ( int j = 0; j < n - 2; j++)
  {
    c[j] = -r;
  }
/*
  Set the initial condition.
*/
  for ( int j = 0; j <= n; j++ )
  {
    double x = prob->a + j * dx;
    u[0][j] = prob->ic ( x );
  }
/*
  Take m timesteps.
*/
  for ( int i = 1; i <= m; i++ )
  {
/*
  Set boundary values.
*/
    double t = i * dt;
    u[i][0] = prob->bcL(t);
    u[i][n] = prob->bcR(t);
/*
  Set the right hand side.
*/
    for ( int j = 1; j <= n - 1; j++ )
    {
      b[j] = 2.0 * ( 1.0 - r ) * u[i-1][j];
    }
    for ( int j = 2; j <= n - 1; j++ )
    {
      b[j] = b[j] + r * u[i-1][j-1];
    }
    for ( int j = 1; j <= n - 2; j++ )
    {
      b[j] = b[j] + r * u[i-1][j+1];
    }
    b[1]   = b[1]   + r * ( u[i-1][0] + u[i][0] );
    b[n-1] = b[n-1] + r * ( u[i-1][n] + u[i][n] );
/*
  Set the diagonal of the matrix.
*/
    for ( int k = 0; k < n - 1; k++ )
    {
      d[k] = 2.0 + 2.0 * r;
    }
/*
  Solve for interior values.
  The "+1" on b and u[i] allows us to skip the first equation which is
  a boundary condition.
*/
    trisolve ( n-1, c, d, c, b+1, u[i]+1 );
  }

  free_vector ( b );
  free_vector ( c );
  free_vector ( d );

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

void heat_solve_explicit ( struct heat_solve *prob )

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

    heat_solve_explicit() uses an explicit method to solve the heat equation.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 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:

    struct heat_solve *prob: data defining the problem.

  Output:

    struct heat_solve *prob: updated data which now includes a solution.
*/
{
  int m = prob->m;
  int n = prob->n;

  double dt = prob->T / m;
  double dx = ( prob->b - prob->a ) / n;
  double r = dt / ( dx * dx );
  double **u = prob->u;

  printf ( "  Explicit finite difference scheme:\n");
  printf ( "  r = %g\n", r );
/*
  Set the initial condition.
*/
  for ( int j = 0; j <= n; j++ )
  {
    double x = prob->a + j * dx;
    u[0][j] = prob->ic ( x );
  }
/*
  Take m timesteps.
*/
  for ( int i = 1; i <= m; i++ )
  {
/*
  Set boundary values.
*/
    double t = i * dt;
    u[i][0] = prob->bcL(t);
    u[i][n] = prob->bcR(t);
/*
  Apply explicit formula.
*/
    for ( int j = 1; j <= n - 1; j++ )
    {
      u[i][j] = u[i-1][j] + r * ( u[i-1][j-1] - 2.0 * u[i-1][j] + u[i-1][j+1] );
    }
  }

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

void heat_solve_implicit ( struct heat_solve *prob )

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

    heat_solve_implicit() uses an implicit method to solve the heat equation.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 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:

    struct heat_solve *prob: data defining the problem.

  Output:

    struct heat_solve *prob: updated data which now includes a solution.
*/
{
  int m = prob->m;
  int n = prob->n;

  double *b;
  double *c;
  double *d;
  double dt = prob->T / m;
  double dx = ( prob->b - prob->a ) / n;
  double r = dt / ( dx * dx );
  double **u = prob->u;

  printf ( "  Implicit finite difference scheme:\n");
  printf ( "  r = %g\n", r );

  make_vector ( d, n-1 );
  make_vector ( c, n-2 );
  make_vector ( b, n+1 );

  for ( int j = 0; j < n - 2; j++)
  {
    c[j] = -r;
  }
/*
  Set the initial condition.
*/
  for ( int j = 0; j <= n; j++ )
  {
    double x = prob->a + j * dx;
    u[0][j] = prob->ic ( x );
  }
/*
  Take m timesteps.
*/
  for ( int i = 1; i <= m; i++ )
  {
/*
  Set boundary values.
*/
    double t = i * dt;
    u[i][0] = prob->bcL(t);
    u[i][n] = prob->bcR(t);
/*
  Set up matrix for interior values.
*/
    for ( int j = 1; j <= n - 1; j++ )
    {
      b[j] = u[i-1][j];
    }
    b[1]   = b[1]   + r * u[i][0];
    b[n-1] = b[n-1] + r * u[i][n];

    for ( int k = 0; k < n - 1; k++ )
    {
      d[k] = 1.0 + 2.0 * r;
    }
/*
  Solve for interior values.
*/
    trisolve ( n-1, c, d, c, b+1, u[i]+1 );
  }

  free_vector ( b );
  free_vector ( c );
  free_vector ( d );

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

void heat_solve_seidman_sweep ( struct heat_solve *prob )

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

    heat_solve_seidman_sweep() uses Seidman sweep to solve the heat equation.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 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:

    struct heat_solve *prob: data defining the problem.

  Output:

    struct heat_solve *prob: updated data which now includes a solution.
*/
{
  int m = prob->m;
  int n = prob->n;

  double dt = prob->T / m;
  double dx = ( prob->b - prob->a ) / n;
  double r = 0.5 * dt / ( dx * dx );
  double t;
  double **u = prob->u;
  double *v;

  make_vector ( v, n + 1 );

  printf ( "  Seidman sweep scheme:\n");
  printf ( "  r = %g\n", r );
/*
  Set the initial condition.
*/
  for ( int j = 0; j <= n; j++ )
  {
    double x = prob->a + j * dx;
    u[0][j] = prob->ic ( x );
  }
/*
  Take m timesteps.
*/
  for ( int i = 1; i <= m; i++ )
  {
/*
  LEFT HALF TIMESTEP
  Set boundary values.
*/
    t = ( 2 * i - 1 ) * dt / 2.0;
    v[0] = prob->bcL(t);
    v[n] = prob->bcR(t);
/*
  Apply explicit formula, including just updated value at j-1.
*/
    for ( int j = 1; j < n; j++ )
    {
      v[j] = r * v[j-1] + ( 1.0 - r ) * u[i-1][j] + r * u[i-1][j+1];
      v[j] = v[j] / ( 1.0 + r );
    }
/*
  RIGHT HALF TIMESTEP
  Set boundary values.
*/
    t = i * dt;
    u[i][0] = prob->bcL(t);
    u[i][n] = prob->bcR(t);
/*
  Apply explicit formula, including just updated value at j+1.
*/
    for ( int j = n - 1; 1 <= j; j-- )
    {
      u[i][j] = r * v[j-1] + ( 1.0 - r ) * v[j] + r * u[i][j+1];
      u[i][j] = u[i][j] / ( 1.0 + r );
    }
  }

  free_vector ( v );

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

void plot_with_geomview ( struct plot3d *data )

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

    plot_with_geomview() plots ppc_fd1() data using geomview().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by 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:

    struct plot3d *data: the plot data.
*/
{
  FILE *fp;
  int m = data->m;
  int n = data->n;
  double dx = ( data->b - data->a ) / n;
  double dt = data->T / m;

  fp = fopen ( data->geomview_out, "w" );

  if ( fp == NULL )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "plot_with_geomview(): Fatal error!\n" );
    fprintf ( stderr, "  Cannot open script file '%s'\n", data->geomview_out );
    return;
  }

  fprintf ( fp, "# Geomview script for ppc_fd1() finite difference solver\n");

  fprintf ( fp, "{ appearance { +edge }\n" );
  fprintf ( fp, "MESH %d %d\n", n+1, m+1 );

  for ( int i = 0; i <= m; i++ )
  {
    double t = i * dt;
    for ( int j = 0; j <= n; j++ )
    {
      double x = data->a + j * dx;
      fprintf ( fp, "%g %g %g\n", x, t, data->u[i][j] );
    }
  }

  fprintf ( fp, "}\n" );

  fclose ( fp );

  printf ( "  Geomview script written to '%s'\n", data->geomview_out );

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

void plot_with_maple ( struct plot3d *data )

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

    plot_with_maple() plots ppc_fd1() data using maple().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by 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:

    struct plot3d *data: the plot data.
*/
{
  FILE *fp;
  int m = data->m;
  int n = data->n;
  double dx = ( data->b - data->a ) / n;
  double dt = data->T / m;

  fp = fopen ( data->maple_out, "w" );

  if ( fp == NULL )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "plot_with_maple(): Fatal error!\n" );
    fprintf ( stderr, "  Cannot open script file '%s'\n", data->maple_out );
    return;
  }

  fprintf ( fp, "# Maple script for finite difference solver\n" );

  fprintf ( fp, "plots:-display(MESH([\n");

  for ( int i = 0; i <= m; i++ )
  {
    double t = i * dt;
    fprintf ( fp, "\t[" );
    for ( int j = 0; j <= n; j++ )
    {
      double x = data->a + j * dx;
      fprintf ( fp, "[%g,%g,%g],", x, t, data->u[i][j] );
    }
    fprintf ( fp, "NULL],\n" );
  }
  fprintf ( fp, "NULL]), labels=[x,t,u]);" );

  fclose ( fp );

  printf ( "  Maple script written to '%s'\n", data->maple_out );

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

void plot_with_matlab ( struct plot3d *data )

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

    plot_with_matlab() plots ppc_fd1() data using matlab().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by 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:

    struct plot3d *data: the plot data.
*/
{
# define XVEC	"PLOT_WITH_MATLAB_FD_XVEC"
# define TVEC	"PLOT_WITH_MATLAB_FD_TVEC"
# define ZMAT	"PLOT_WITH_MATLAB_FD_ZMAT"

  FILE *fp;
  int m = data->m;
  int n = data->n;
  double dx = ( data->b - data->a ) / n;
  double dt = data->T / m;

  fp = fopen ( data->matlab_out, "w" );

  if ( fp == NULL )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "plot_with_matlab(): Fatal error!\n" );
    fprintf ( stderr, "  Cannot open script file '%s' \n", data->matlab_out );
    return;
  }

  fprintf ( fp, "%% Matlab plot script for ppc_fd1() finite difference solver.\n");

  fprintf ( fp, "%s = [\n", XVEC );
  for ( int j = 0; j <= n; j++ )
  {
    double x = data->a + j * dx;
    fprintf ( fp, "%g ", x );
  }
  fprintf ( fp, "\n];\n" );

  fprintf ( fp, "%s = [\n", TVEC );
  for ( int i = 0; i <= m; i++ )
  {
    double t = i * dt;
    fprintf ( fp, "%g ", t );
  }
  fprintf ( fp, "\n];\n" );

  fprintf ( fp, "%s = [\n", ZMAT );
  for ( int i = 0; i <= m; i++ )
  {
    for ( int j = 0; j <= n; j++ )
    {
      fprintf ( fp, "%g ", data->u[i][j] );
    }
    fprintf ( fp, "\n" );
  }
  fprintf ( fp, "];\n" );

  fprintf ( fp, "surf ( %s, %s, %s )\n", XVEC, TVEC, ZMAT );
  fprintf ( fp, "xlabel('x')\n" );
  fprintf ( fp, "ylabel('t')\n" );
  fprintf ( fp, "zlabel('u')\n" );
  fprintf ( fp, "colormap jet\n" );

  fclose ( fp );

  printf ( "  Matlab script written to '%s'\n", data->matlab_out );

  return;

#undef XVEC
#undef TVEC
#undef ZMAT
}
/******************************************************************************/

void plot3d ( struct plot3d *data )

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

    plot3d() plots ppc_fd1() data using geomview(), maple(), or matlab().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by 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:

    struct plot3d *data: the plot data.
*/
{
  if ( data->geomview_out != NULL )
  {
    plot_with_geomview ( data );
  }

  if ( data->maple_out != NULL )
  {
    plot_with_maple ( data );
  }

  if ( data->matlab_out != NULL )
  {
    plot_with_matlab ( data );
  }

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

void show_usage_and_exit ( char *progname )

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

    show_usage_and_exit() prints a usage message.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 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:

    char *progname: the program name as invoked by the user.
*/
{
  fprintf ( stderr, "\n" );
  fprintf ( stderr, "Usage: %s T n m\n", progname );
  fprintf ( stderr, "  T : solve over 0 ≤ t ≤ T\n" );
  fprintf ( stderr, "  n : number of subintervals in the x direction\n" );
  fprintf ( stderr, "  m : number of subintervals in the t direction\n" );

  exit ( EXIT_SUCCESS );
}
/******************************************************************************/

void trisolve ( int n, const double *a, double *d, const double *c, double *b, 
  double *x )

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

    trisolve() solves a tridiagonal linear system.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 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 n:

    double a[n], double d[n], double c[n]

    double b[n]: the right hand side.

  Output:

    double b[n]: altered during the solution process.

    double d[n]: altered during the solution process.

    double x[n]: the solution.
*/
{
  double m;

  for ( int i = 1; i < n; i++ ) 
  {
    m = a[i-1] / d[i-1];
    d[i] = d[i] - m * c[i-1];
    b[i] = b[i] - m * b[i-1];
  }

  x[n-1] = b[n-1] / d[n-1];
  for ( int i = n - 2; 0 <= i; i-- )
  {
    x[i] = ( b[i] - c[i] * x[i+1] ) / d[i];
  }

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

void write_plotting_script ( struct heat_solve *prob )

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

    write_plotting_script() prepares solution data to be plotted.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 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:

    struct heat_solve *prob: the finite difference solution data.
*/
{
  struct plot3d data;

  data.a = prob->a;
  data.b = prob->b;
  data.T = prob->T;
  data.n = prob->n;
  data.m = prob->m;
  data.u = prob->u;
  data.maple_out = prob->maple_out;
  data.matlab_out = prob->matlab_out;
  data.geomview_out = prob->geomview_out;

  plot3d ( &data );
}