subroutine doughnut_exact ( nt, t, y )

!*****************************************************************************80
!
!! doughnut_exact(): exact solution of doughnut ODE.
!
!  Discussion:
!
!    The formula assumes that t0 = 0.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 October 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer nt: the number of times.
!
!    real ( kind = rk8 ) t(nt): the times.
!
!  Output:
!
!    real ( kind = rk8 ) y(nt,3): the derivative values.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer nt

  real ( kind = rk8 ) a
  real ( kind = rk8 ) b
  real ( kind = rk8 ) c
  real ( kind = rk8 ) delta
  real ( kind = rk8 ) m
  real ( kind = rk8 ) n
  real ( kind = rk8 ) t(nt)
  real ( kind = rk8 ) y(nt,3)
  real ( kind = rk8 ) y0(3)
!
!  This nonsense is another failure of Fortran:
!
  interface
    subroutine doughnut_parameters ( &
      m_in,  n_in,  t0_in,  y0_in,  tstop_in, &
      m_out, n_out, t0_out, y0_out, tstop_out )
      integer, parameter :: rk8 = kind ( 1.0D+00 )
      real ( kind = rk8 ), optional :: m_in
      real ( kind = rk8 ), optional :: m_out
      real ( kind = rk8 ), optional :: n_in
      real ( kind = rk8 ), optional :: n_out
      real ( kind = rk8 ), optional :: t0_in
      real ( kind = rk8 ), optional :: t0_out
      real ( kind = rk8 ), optional :: y0_in(3)
      real ( kind = rk8 ), optional :: y0_out(3)
      real ( kind = rk8 ), optional :: tstop_in
      real ( kind = rk8 ), optional :: tstop_out
    end subroutine
  end interface

  call doughnut_parameters ( m_out = m, n_out = n, y0_out = y0 )

  a = y0(1)
  b = y0(2)
  c = y0(3)

  delta = 1.0 + a**2 + b**2 + c**2
  
  y(1:nt,1) = ( 2.0 * a * cos ( m * t ) - 2.0 * b * sin ( m * t ) ) &
         / ( delta - 2.0 * c * sin ( n * t ) + ( 2.0 - delta ) * cos ( n * t ) )

  y(1:nt,2) = ( 2.0 * a * sin ( m * t ) + 2.0 * b * cos ( m * t ) ) &
         / ( delta - 2.0 * c * sin ( n * t ) + ( 2.0 - delta ) * cos ( n * t ) )

  y(1:nt,3) = ( 2.0 * c * cos ( n * t ) + ( 2.0 - delta ) * sin ( n * t ) ) &
         / ( delta - 2.0 * c * sin ( n * t ) + ( 2.0 - delta ) * cos ( n * t ) )

  return
end
subroutine doughnut_parameters ( &
  m_in,  n_in,  t0_in,  y0_in,  tstop_in, &
  m_out, n_out, t0_out, y0_out, tstop_out )

!*****************************************************************************80
!
!! 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:
!
!    10 October 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) m_in, n_in: parameters.
!
!    real ( kind = rk8 ) t0_in, y0_in: the initial time and value.
!
!    real ( kind = rk8 ) tstop_in: the final time.
!
!  Persistent:
!
!    real ( kind = rk8 ) m_default, n_default.
!
!    real ( kind = rk8 ) t0_default.
!
!    real ( kind = rk8 ) y0_default.
!
!    real ( kind = rk8 ) tstop_default.
!
!  Output:
!
!    real ( kind = rk8 ) m_out, n_out: parameters.
!
!    real ( kind = rk8 ) t0_out, y0_out: the initial time and value.
!
!    real ( kind = rk8 ) tstop_out: the final time.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), save     :: m_default = 3.0D+00
  real ( kind = rk8 ), optional :: m_in
  real ( kind = rk8 ), optional :: m_out
  real ( kind = rk8 ), save     :: n_default = 5.0D+00
  real ( kind = rk8 ), optional :: n_in
  real ( kind = rk8 ), optional :: n_out
  real ( kind = rk8 ), save     :: t0_default = 0.0D+00
  real ( kind = rk8 ), optional :: t0_in
  real ( kind = rk8 ), optional :: t0_out
  real ( kind = rk8 ), save     :: y0_default(3) = (/ 1.0D+00, 1.0D+00, 3.0D+00 /)
  real ( kind = rk8 ), optional :: y0_in(3)
  real ( kind = rk8 ), optional :: y0_out(3)
  real ( kind = rk8 ), save     :: tstop_default = 10.0D+00
  real ( kind = rk8 ), optional :: tstop_in
  real ( kind = rk8 ), optional :: tstop_out
!
!  New values, if supplied on input, overwrite the current values.
!
  if ( present ( m_in ) ) then
    m_default = m_in
  end if

  if ( present ( n_in ) ) then
    n_default = n_in
  end if

  if ( present ( t0_in ) ) then
    t0_default = t0_in
  end if

  if ( present ( y0_in ) ) then
    y0_default = y0_in
  end if

  if ( present ( tstop_in ) ) then
    tstop_default = tstop_in
  end if
!
!  The current values are copied to the output.
!
  if ( present ( m_out ) ) then
    m_out = m_default
  end if

  if ( present ( n_out ) ) then
    n_out = n_default
  end if

  if ( present ( t0_out ) ) then
    t0_out = t0_default
  end if

  if ( present ( y0_out ) ) then
    y0_out = y0_default
  end if

  if ( present ( tstop_out ) ) then
    tstop_out = tstop_default
  end if

  return
end