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

using namespace std;

# include "doughnut_exact.hpp"

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

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

//****************************************************************************80
//
//  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:
//
//    13 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 = new double[nt*3];

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