/*************************************************************************
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier

Contributors:
    * Sergey Bochkanov (ALGLIB project). Translation from C to
      pseudocode.

See subroutines comments for additional copyrights.

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/

#include <stdafx.h>
#include "jacobianelliptic.h"

/*************************************************************************
Jacobian Elliptic Functions

Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
and dn(u|m) of parameter m between 0 and 1, and real
argument u.

These functions are periodic, with quarter-period on the
real axis equal to the complete elliptic integral
ellpk(1.0-m).

Relation to incomplete elliptic integral:
If u = ellik(phi,m), then sn(u|m) = sin(phi),
and cn(u|m) = cos(phi).  Phi is called the amplitude of u.

Computation is by means of the arithmetic-geometric mean
algorithm, except when m is within 1e-9 of 0 or 1.  In the
latter case with m close to 1, the approximation applies
only for phi < pi/2.

ACCURACY:

Tested at random points with u between 0 and 10, m between
0 and 1.

           Absolute error (* = relative error):
arithmetic   function   # trials      peak         rms
   IEEE      phi         10000       9.2e-16*    1.4e-16*
   IEEE      sn          50000       4.1e-15     4.6e-16
   IEEE      cn          40000       3.6e-15     4.4e-16
   IEEE      dn          10000       1.3e-12     1.8e-14

 Peak error observed in consistency check using addition
theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
the above relation to the incomplete elliptic integral.
Accuracy deteriorates when u is large.

Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*************************************************************************/
void jacobianellipticfunctions(double u,
     double m,
     double& sn,
     double& cn,
     double& dn,
     double& ph)
{
    double ai;
    double b;
    double phi;
    double t;
    double twon;
    ap::real_1d_array a;
    ap::real_1d_array c;
    int i;

    ap::ap_error::make_assertion(ap::fp_greater_eq(m,0)&&ap::fp_less_eq(m,1), "Domain error in JacobianEllipticFunctions: m<0 or m>1");
    a.setbounds(0, 8);
    c.setbounds(0, 8);
    if( ap::fp_less(m,1.0e-9) )
    {
        t = sin(u);
        b = cos(u);
        ai = 0.25*m*(u-t*b);
        sn = t-ai*b;
        cn = b+ai*t;
        ph = u-ai;
        dn = 1.0-0.5*m*t*t;
        return;
    }
    if( ap::fp_greater_eq(m,0.9999999999) )
    {
        ai = 0.25*(1.0-m);
        b = cosh(u);
        t = tanh(u);
        phi = 1.0/b;
        twon = b*sinh(u);
        sn = t+ai*(twon-u)/(b*b);
        ph = 2.0*atan(exp(u))-1.57079632679489661923+ai*(twon-u)/b;
        ai = ai*t*phi;
        cn = phi-ai*(twon-u);
        dn = phi+ai*(twon+u);
        return;
    }
    a(0) = 1.0;
    b = sqrt(1.0-m);
    c(0) = sqrt(m);
    twon = 1.0;
    i = 0;
    while(ap::fp_greater(fabs(c(i)/a(i)),ap::machineepsilon))
    {
        if( i>7 )
        {
            ap::ap_error::make_assertion(false, "Overflow in JacobianEllipticFunctions");
            break;
        }
        ai = a(i);
        i = i+1;
        c(i) = 0.5*(ai-b);
        t = sqrt(ai*b);
        a(i) = 0.5*(ai+b);
        b = t;
        twon = twon*2.0;
    }
    phi = twon*a(i)*u;
    do
    {
        t = c(i)*sin(phi)/a(i);
        b = phi;
        phi = (asin(t)+phi)/2.0;
        i = i-1;
    }
    while(i!=0);
    sn = sin(phi);
    t = cos(phi);
    cn = t;
    dn = t/cos(phi-b);
    ph = phi;
}




