# include <cmath>
# include <cstdlib>

using namespace std;

# include "kdv_exact.hpp"

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

void kdv_exact_rational ( double x, double t, double *u, double *ut, double *ux, 
  double *uxx, double *uxxx )

//****************************************************************************80
//
//  Purpose:
//
//    kdv_exact_rational(): exact solution of KDV PDE using a rational function.
//
//  Discussion:
//
//    This solution u(x,t) satisfies the Korteweg-Devries PDE:
//      u' - 6 u ux + uxxx = 0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John D Cook,
//    Rational solution to the Korteweg-De Vries equation,
//    13 November 2023,
//    https://www.johndcook.com/blog/2023/11/13/rational-kdv/
//
//  Input:
//
//    double x: the position.
//
//    double t: the time.
//
//  Output:
//
//    double *u, *ut, *ux, *uxx, *uxxx:
//    the values of the solution and its derivatives.
//
{
  double t2;
  double t3;
  double x2;
  double x3;
  double x5;
  double x6;
  double x8;
  double x9;

  t2 = pow ( t, 2 );
  t3 = pow ( t, 3 );
  x2 = pow ( x, 2 );
  x3 = pow ( x, 3 );
  x5 = pow ( x, 5 );
  x6 = pow ( x, 6 );
  x8 = pow ( x, 8 );
  x9 = pow ( x, 9 );

  *u = 6.0 * x * ( x3 - 24.0 * t ) / pow ( x3 + 12.0 * t, 2 );

  *ut = - 288.0 * x * ( x3 - 6.0 * t ) / pow ( x3 + 12.0 * t, 3 );

  *ux = - 12.0 * ( x6 - 84.0 * t * x3 + 144.0 * t2 ) / pow ( x3 + 12.0 * t, 3 );

  *uxx =  36.0 * ( x8 - 192.0 * t * x5 + 1440.0 * t2 * x2 ) 
    / pow ( x3 + 12.0 * t, 4 );

  *uxxx = - 144.0 * x * ( x9 - 360.0 * t * x6 + 6480.0 * t2 * x3 - 8640.0 * t3 ) 
    / pow ( x3 + 12.0 * t, 5 );

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

void kdv_exact_sech ( double x, double t, double *u, double *ut, double *ux, 
  double *uxx, double *uxxx )

//****************************************************************************80
//
//  Purpose:
//
//    kdv_exact_sech(): exact solution of KDV PDE using sech().
//
//  Discussion:
//
//    This solution u(x,t) satisfies the Korteweg-Devries PDE:
//      u' - 6 u ux + uxxx = 0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John D Cook,
//    Solitons and the KdV equation,
//    03 November 2023,
//    https://www.johndcook.com/blog/2023/11/03/solitons-and-the-kdv-equation/
//
//  Input:
//
//    double x: the position.
//
//    double t: the time.
//
//  Output:
//
//    double *u, *ut, *ux, *uxx, *uxxx:
//    the values of the solution and its derivatives.
//
{
  double a;
  double argument;
  double sa;
  double ta;
  double v;
//
//  Retrieve parameters a and v.
//
  kdv_parameters ( NULL, NULL, NULL, NULL, &a, &v, NULL, NULL );

  argument = 0.5 * sqrt ( v ) * ( x - v * t - a );
  sa = 1.0 / cosh ( argument );
  ta = tanh ( argument );

  *u =    - 0.5  *       v        * pow ( sa, 2 );

  *ut =   - 0.5  * pow ( v, 2.5 ) * pow ( sa, 2 ) * ta;

  *ux =     0.5  * pow ( v, 1.5 ) * pow ( sa, 2 ) * ta;

  *uxx =  + 0.25 * pow ( v, 2 )   * pow ( sa, 4 ) 
          - 0.5  * pow ( v, 2 )   * pow ( sa, 2 ) * pow ( ta, 2 );

  *uxxx =        - pow ( v, 2.5 ) * pow ( sa, 4 ) * ta 
          + 0.5  * pow ( v, 2.5 ) * pow ( sa, 2 ) * pow ( ta, 3 );

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

void kdv_parameters ( double *a_in, double *v_in, double *t0_in, 
  double *tstop_in, double *a_out, double *v_out, double *t0_out, 
  double *tstop_out )

//****************************************************************************80
//
//  Purpose:
//
//    kdv_parameters(): parameters for the Korteweg-Devries PDE.
//
//  Discussion:
//
//    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:
//
//    02 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double *a_user: the phase
//
//    double *v_user: the velocity
//
//    double *t0_user: the initial time
//
//    double *tstop_user: the final time.
//
//  Output:
//
//    double *a_out: the phase.
//
//    double *v_out: the velocity
//
//    double *t0_out: the initial time
//
//    double *tstop_out: the final time.
//
{
  static double a_default = 0.0;
  static double t0_default = 0.0;
  static double tstop_default = 10.0;
  static double v_default = 1.0;
//
//  New values, if supplied on input, overwrite the current values.
//
  if ( a_in )
  {
    a_default = *a_in;
  }

  if ( v_in )
  {
    v_default = *v_in;
  }

  if ( t0_in )
  {
    t0_default = *t0_in;
  }

  if ( tstop_in )
  {
    tstop_default = *tstop_in;
  }
//
//  The current values are copied to the output.
//
  if ( a_out )
  {
    *a_out = a_default;
  }

  if ( v_out )
  {
    *v_out = v_default;
  }

  if ( t0_out )
  {
    *t0_out = t0_default;
  }

  if ( tstop_out )
  {
    *tstop_out = tstop_default;
  }
  
  return;
}
//****************************************************************************80

double kdv_residual ( double u, double ut, double ux, double uxxx )

//****************************************************************************80
//
//  Purpose:
//
//    kdv_residual(): evaluate the residual of KDV PDE.
//
//  Discussion:
//
//    The residual of the Korteweg-Devries PDE:
//      u' - 6 u ux + uxxx
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John D Cook,
//    Solitons and the KdV equation,
//    03 November 2023,
//    https://www.johndcook.com/blog/2023/11/03/solitons-and-the-kdv-equation/
//
//  Input:
//
//    double u, ut, ux, uxxx:
//    the values of the solution and its derivatives.
//
//  Output:
//
//    double r: the residual.
//
{
  double r;

  r = ut - 6.0 * u * ux + uxxx;

  return r;
}