function [ sn, cn, dn ] = sncndn ( u, m )

%*****************************************************************************80
%
%% SNCNDN evaluates Jacobi elliptic functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2018
%
%  Author:
%
%    Original ALGOL version by Roland Bulirsch.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Roland Bulirsch,
%    Numerical calculation of elliptic integrals and elliptic functions,
%    Numerische Mathematik,
%    Volume 7, Number 1, 1965, pages 78-90.
%
%  Parameters:
%
%    Input, real U, M, the arguments.
%
%    Output, real SN, CN, DN, the value of the Jacobi
%    elliptic functions sn(u,m), cn(u,m), and dn(u,m).
%
  m_comp = 1.0 - m;
  u_copy = u;

  if ( m_comp == 0.0 )
    cn = 1.0 / cosh ( u_copy );
    dn = cn;
    sn = tanh ( u_copy );
    return
  end

  if ( 1.0 < m )
    d = 1.0 - m_comp;
    m_comp = - m_comp / d;
    d = sqrt ( d );
    u_copy = d * u_copy;
  end

  ca = sqrt ( eps );

  a = 1.0;
  dn = 1.0;
  l = 25;

  m_array = zeros ( 25, 1 );
  n_array = zeros ( 25, 1 );

  for i = 1 : 25

    m_array(i) = a;
    m_comp = sqrt ( m_comp );
    n_array(i) = m_comp;
    c = 0.5 * ( a + m_comp );

    if ( abs ( a - m_comp ) <= ca * a )
      l = i;
      break;
    end

    m_comp = a * m_comp;
    a = c;

  end

  u_copy = c * u_copy;
  sn = sin ( u_copy );
  cn = cos ( u_copy );

  if ( sn ~= 0.0 )

    a = cn / sn;
    c = a * c;

    for i = l : -1 : 1
      b = m_array(i);
      a = c * a;
      c = dn * c;
      dn = ( n_array(i) + a ) / ( b + a );
      a = c / b;
    end

    a = 1.0 / sqrt ( c * c + 1.0 );

    if ( sn < 0.0 )
      sn = - a;
    else
      sn = a;
    end

    cn = c * sn;

  end

  if ( 1.0 < m )
    a = dn;
    dn = cn;
    cn = a;
    sn = sn / d;
  end

  return
end