# include <cmath>
# include <iomanip>
# include <iostream>

using namespace std;

# include "doughnut_ode.hpp"

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

void doughnut_deriv ( double t, double *y, double dydt[] )

//****************************************************************************80
//
//  doughnut_deriv() evaluates the right hand side of the doughnut ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    14 October 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double T, Y[3]: the time and solution value.
//
//  Output:
//
//    double DYDT[3]: the derivative value.
//
{
  double m;
  double n;

  doughnut_parameters ( NULL, NULL, NULL, NULL, NULL,
                        &m,   &n,   NULL, NULL, NULL );

  dydt[0] = - m * y[1] + n * y[0] * y[2];
  dydt[1] =   m * y[0] + n * y[1] * y[2];
  dydt[2] =   0.5 * n * ( 1.0 - y[0] * y[0] - y[1] * y[1] + y[2] * y[2] );

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

void doughnut_parameters ( 
  double *m_in,  double *n_in,  double *t0_in,  double *y0_in,  double *tstop_in, 
  double *m_out, double *n_out, double *t0_out, double *y0_out, double *tstop_out )

//****************************************************************************80
//
//  Purpose:
//
//    doughnut_parameters() returns parameters of the doughnut ODE.
//
//  Discussion:
//
//    If input values are specified, this resets the default parameters.
//    Otherwise, the output will be the current defaults.
//
//    Suggested choices for m, n, (y0) include:
//
//      3,  5, (1, 1, 3.0)
//      4,  5, (1, 1, 0.3)
//      pi, 5, (1, 1, 0.3)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 October 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double m_in, n_in: parameters.
//
//    double t0_in, y0_in[3]: the initial time and value.
//
//    double tstop_in: the final time.
//
//  Persistent:
//
//    double m_default, n_default.
//
//    double t0_default.
//
//    double y0_default[3].
//
//    double tstop_default.
//
//  Output:
//
//    double m_out, n_out: parameters.
//
//    double t0_out, y0_out[3]: the initial time and value.
//
//    double tstop_out: the final time.
//
{
  int i;
  static double m_default = 3.0;
  static double n_default = 5.0;
  static double t0_default = 0.0;
  static double y0_default[3] = { 1.0, 1.0, 3.0 };
  static double tstop_default = 10.0;
//
//  New values, if supplied on input, overwrite the current values.
//
  if ( m_in )
  {
    m_default = *m_in;
  }

  if ( n_in )
  {
    n_default = *n_in;
  }

  if ( t0_in )
  {
    t0_default = *t0_in;
  }
  if ( y0_in )
  {
    for ( i = 0; i < 3; i++ )
    {
      y0_default[i] = y0_in[i];
    }
  }

  if ( tstop_in )
  {
    tstop_default = *tstop_in;
  }
//
//  The current values are copied to the output.
//
  if ( m_out )
  {
    *m_out = m_default;
  }
  if ( n_out )
  {
    *n_out = n_default;
  }
  if ( t0_out )
  {
    *t0_out = t0_default;
  }
  if ( y0_out )
  {
    for ( i = 0; i < 3; i++ )
    {
      y0_out[i] = y0_default[i];
    }
  }
  if ( tstop_out )
  {
    *tstop_out = tstop_default;
  }

  return;
}