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

# include "doughnut_exact.h"

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

double *doughnut_exact ( int nt, double t[] )

/******************************************************************************/
/*
  Purpose:

    doughnut_exact(): exact solution of doughnut ODE.

  Discussion:

    The formula assumes that t0 = 0.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    12 October 2025

  Author:

    John Burkardt

  Input:

    int nt: the number of times.

    double t[nt]: the times.

  Output:

    double doughnut_exact[nt*3]: the derivative values.
*/
{
  double a;
  double b;
  double c;
  double delta;
  int i;
  double m;
  double n;
  double *y;
  double y0[3];

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

  a = y0[0];
  b = y0[1];
  c = y0[2];

  delta = 1.0 + a * a + b * b + c * c;
  
  y = ( double * ) malloc ( nt * 3 * sizeof ( double ) );

  for ( i = 0; i < nt; i++ )
  {
    y[i*3+0] = ( 2.0 * a * cos ( m * t[i] ) - 2.0 * b * sin ( m * t[i] ) ) 
           / ( delta - 2.0 * c * sin ( n * t[i] ) + ( 2.0 - delta ) * cos ( n * t[i] ) );

    y[i*3+1] = ( 2.0 * a * sin ( m * t[i] ) + 2.0 * b * cos ( m * t[i] ) ) 
           / ( delta - 2.0 * c * sin ( n * t[i] ) + ( 2.0 - delta ) * cos ( n * t[i] ) );

    y[i*3+2] = ( 2.0 * c * cos ( n * t[i] ) + ( 2.0 - delta ) * sin ( n * t[i] ) ) 
           / ( delta - 2.0 * c * sin ( n * t[i] ) + ( 2.0 - delta ) * cos ( n * t[i] ) );
  }

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

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;
  }
/*
  How is this to be handled?
*/
  if ( y0_out )
  {
    for ( i = 0; i < 3; i++ )
    {
      y0_out[i] = y0_default[i];
    }
  }
  if ( tstop_out )
  {
    *tstop_out = tstop_default;
  }

  return;
}