# include <cmath>
# include <cstdlib>

using namespace std;

# include "traveling_wave_exact.hpp"

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

void traveling_wave_exact ( double x, double t, double &u, double &ut, 
  double &utt, double &ux, double &uxx )

//****************************************************************************80
//
//  Purpose:
//
//    traveling_wave_exact() evaluates an exact solution of the 1D wave equation.
//
//  Discussion:
//
//    d^2u/dt^2 = c^2 d^2u/dx^2
//    u[x,t] = a * sin ( k * x + w * t + phi )
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    12 April 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x, t: the position and time.
//
//  Output:
//
//    real u, ut, utt, ux, uxx: the values of the exact solution, its time
//    derivative, and its first and second spatial derivatives at (x,t).
//
{
  double a;
  double k;
  double phi;
  double t0;
  double tstop;
  double w;
  double xmax;
  double xmin;

  traveling_wave_parameters ( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
    &a, &k, &w, &phi, &xmin, &xmax, &t0, &tstop );

  u = a * sin ( k * x + w * t + phi );

  ut = a * w * cos ( k * x + w * t + phi );

  utt = - a * w * w * sin ( k * x + w * t + phi );

  ux = a * k * cos ( k * x + w * t + phi );

  uxx = - a * k * k * sin ( k * x + w * t + phi );

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

void traveling_wave_parameters ( double *a_in, double *k_in, double *w_in, 
  double *phi_in, double *xmin_in, double *xmax_in, double *t0_in, 
  double *tstop_in, double *a_out, double *k_out, double *w_out, 
  double *phi_out, double *xmin_out, double *xmax_out, double *t0_out, 
  double *tstop_out )

//****************************************************************************80
//
//  Purpose:
//
//    traveling_wave_parameters() returns parameters for the traveling wave equation.
//
//  Discussion:
//
//    d^2u/dt^2 = c^2 d^2u/dx^2
//
//    If input values are specified, this resets the default parameters.
//    Otherwise, the output will be the current defaults.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 April 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real a_in: the amplitude.
//
//    real k_in: the wave number.
//
//    real w_in: angular frequency.
//
//    real phi_in: the phase angle.
//
//    real xmin_in, xmax_in: the left and right interval endpoints.
//
//    real t0_in: the initial time.
//
//    real tstop_in: the final time.
//
//  Output:
//
//    real a_out: the amplitude.
//
//    real k_out: the wave number.
//
//    real w_out: angular frequency.
//
//    real phi_out: the phase angle.
//
//    real xmin_out, xmax_out: the left and right interval endpoints.
//
//    real t0_out: the initial time.
//
//    real tstop_out: the final time.
//
{
  static double a_default = 2.0;
  static double k_default = 6.0;
  static double phi_default = 3.141592653589793/ 4.0;
  static double t0_default = 0.0;
  static double tstop_default = 5.0;
  static double w_default = 1.0;
  static double xmax_default = +1.0;
  static double xmin_default = -1.0;
//
//  New values, if supplied on input, overwrite the current values.
//
  if ( a_in )
  {
    a_default = *a_in;
  }

  if ( k_in )
  {
    k_default = *k_in;
  }

  if ( w_in )
  {
    w_default = *w_in;
  }

  if ( phi_in )
  {
    phi_default = *phi_in;
  }

  if ( xmin_in )
  {
    xmin_default = *xmin_in;
  }

  if ( xmax_in )
  {
    xmax_default = *xmax_in;
  }

  if ( t0_in )
  {
    t0_default = *t0_in;
  }

  if ( tstop_in )
  {
    tstop_default = *tstop_in;
  }
//
//  Return values.
//
  if ( a_out )
  {
    *a_out = a_default;
  }
  if ( k_out )
  {
    *k_out = k_default;
  }
  if ( w_out )
  {
    *w_out = w_default;
  }
  if ( phi_out )
  {
    *phi_out = phi_default;
  }
  if ( xmin_out )
  {
    *xmin_out = xmin_default;
  }
  if ( xmax_out )
  {
    *xmax_out = xmax_default;
  }
  if ( t0_out )
  {
    *t0_out = t0_default;
  }
  if ( tstop_out )
  {
    *tstop_out = tstop_default;
  }
  return;
}
//****************************************************************************80

void traveling_wave_residual ( double t, double x, double &r )

//****************************************************************************80
//
//  Purpose:
//
//    traveling_wave_residual() computes the residual of the traveling wave equation.
//
//  Discussion:
//
//    d^2u/dt^2 = c^2 d^2u/dx^2
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 April 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real T, X: the time and position where the solution is evaluated.
//
//  Output:
//
//    double &r: the residual at that time and position.
//
{
  double a;
  double c;
  double k;
  double phi;
  double t0;
  double tstop;
  double u;
  double ut;
  double utt;
  double ux;
  double uxx;
  double w;
  double xmax;
  double xmin;

  traveling_wave_parameters ( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
    &a, &k, &w, &phi, &xmin, &xmax, &t0, &tstop );

  traveling_wave_exact ( t, x, u, ut, utt, ux, uxx );

  c = w / k;

  r = utt - c * c * uxx;

  return;
}