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

# include "stiff_exact.h"

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

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

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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 = ( double * ) malloc ( n * sizeof ( double ) );

  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;
}
/******************************************************************************/

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 )

/******************************************************************************/
/*
  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;
}