# include <cmath>
# include <cstdlib>

using namespace std;

# include "stiff_exact.hpp"

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

void stiff_deriv ( double t, double y[], double dydt[] )

//****************************************************************************80
//
//  stiff_deriv() evaluates the right hand side of the stiff ODE.
//
//  Discussion:
//
//    y' = lambda * ( cos(t-t0) - y )
//    y(t0) = y0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double T, Y[1]: the time and solution value.
//
//  Output:
//
//    double DYDT[1]: the derivative value.
//
{
  double lambda;
  double t0;

  stiff_parameters ( NULL, NULL, NULL, NULL, &lambda, &t0, NULL, NULL );

  dydt[0] = lambda * ( cos ( t - t0 ) - y[0] );

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

double *stiff_exact ( int n, double t[] )

//****************************************************************************80
//
//  Purpose:
//
//    stiff_exact() evaluates the exact solution of the stiff ODE.
//
//  Discussion:
//
//    y' = lambda * ( cos(t-t0) - y )
//    y(t0) = y0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N: the number of values.
//
//    double T[N]: the evaluation times.
//
//  Output:
//
//    double STIFF_EXACT[N]: the exact solution values.
//
{
  int i;
  double lambda;
  double mu;
  double t0;
  double *y;
  double y0;

  stiff_parameters ( NULL, NULL, NULL, NULL, &lambda, &t0, &y0, NULL );

  mu = y0 - lambda * lambda / ( lambda * lambda + 1.0 );

  y = new double[n];

  for ( i = 0; i < n; i++ )
  {
    y[i] = lambda          * sin ( t[i] - t0 ) / ( lambda * lambda + 1.0 ) 
         + lambda * lambda * cos ( t[i] - t0 ) / ( lambda * lambda + 1.0 )
         + mu * exp ( - lambda * t[i] - t0 );
  }

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

void stiff_parameters ( 
  double *lambda_in,  double *t0_in,  double *y0_in,  double *tstop_in, 
  double *lambda_out, double *t0_out, double *y0_out, double *tstop_out )

//****************************************************************************80
//
//  Purpose:
//
//    stiff_parameters() returns parameters of the stiff ODE.
//
//  Discussion:
//
//    y' = lambda * ( cos(t-t0) - y )
//    y(t0) = y0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double *lambda_in, a parameter.
//
//    double *t0_in, double *y0_in: the initial time and value.
//
//    double *tstop_in: the final time.
//
//  Output:
//
//    double *lambda_out, a parameter.
//
//    double *t0_out, double *y0_out: the initial time and value.
//
//    double *tstop_out: the final time.
//
{
  static double lambda_default = 50.0;
  static double t0_default = 0.0;
  static double y0_default = 0.0;
  static double tstop_default = 1.0;
//
//  New values, if supplied on input, overwrite the current values.
//
  if ( lambda_in )
  {
    lambda_default = *lambda_in;
  }

  if ( t0_in )
  {
    t0_default = *t0_in;
  }

  if ( y0_in )
  {
    y0_default = *y0_in;
  }

  if ( tstop_in )
  {
    tstop_default = *tstop_in;
  }
//
//  The current values are copied to the output.
//
  if ( lambda_out )
  {
    *lambda_out = lambda_default;
  }
  if ( t0_out )
  {
    *t0_out = t0_default;
  }
  if ( y0_out )
  {
    *y0_out = y0_default;
  }
  if ( tstop_out )
  {
    *tstop_out = tstop_default;
  }

  return;
}