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

# include "doughnut_ode.h"

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

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

/******************************************************************************/
/*
  doughnut_deriv() evaluates the right hand side of the doughnut ODE.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

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

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 )

/******************************************************************************/
/*
  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:

    11 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;
}