// ell_ag.cpp -- Complete and incomplete elliptic integrals of
//      the 1st and 2nd kinds. Algorithms and coefficient values
//      from "Computation of Special Functions", Zhang and Jin,
//      John Wiley and Sons, 1996.
//
//  (C) 2005, C. Bond. All rights reserved.
//
//
#include <iostream.h>
#include <math.h>
//
//  int ellit(double phi,double k,double *fe,double *ee);
//
//  Input:
//          phi -- argument in radians (pi/2 for complete elliptic integral)
//          k  -- (k) modulus, (0.0 - 1.0)
//
//  Returns: (error code)
//          0 -- no error
//          1 -- bad modulus
//          2 -- bad argument
//
//  Output:
//          fe -- F(k,phi)  (in)complete Elliptic function of 1st kind
//          ee -- E(k,phi)  (in)complete Elliptic function fo 2nd kind
//
//
//  This routine has been modified to accept an argument (phi) in radians,
//  rather than degrees. It has also been modified to take the argument
//  as the first parameter and the modulus as the second.
//
//  Some treatments use the modulus 'm', rather than 'k' where m = k^2.
//
//
int ellit(double phi,double k,double *fe,double *ee)
{
    double g,fac,a,b,c,d,a0,b0,d0,r,ck,ce;

    if ((k < 0.0) || (k > 1.0)) return 1;
    if ((phi < 0.0) || (phi > 0.5*M_PI)) return 2;

    g = 0.0;
    a0 = 1.0;
    r = k*k;
    b0 = sqrt(1.0-r);
    d0 = phi;
    if ((k == 1.0) && (phi == 0.5*M_PI)) {
        *fe = 1.0e300;
        *ee = 1.0;
    }
    else if (k == 1.0) {
        *fe = log((1.0+sin(d0))/cos(d0));
        *ee = sin(d0);
    }
    else {
        fac = 1.0;
        for (int n=1;n<41;n++) {
            a = 0.5*(a0+b0);
            b = sqrt(a0*b0);
            c = 0.5*(a0-b0);
            fac *= 2.0;
            r += (fac*c*c);
            if (phi != 0.5*M_PI) {
                d = d0+atan((b0/a0)*tan(d0));
                g += (c*sin(d));
                d0 = d+M_PI*int(d/M_PI+0.5);
            }
            a0 = a;
            b0 = b;
            if (c < 1.0e-7) break;
        }
        ck = M_PI/(2.0*a);
        ce = M_PI*(2.0-r)/(4.0*a);
        if (phi == 90.0) {
            *fe = ck;
            *ee = ce;
        }
        else {
            *fe = d/(fac*a);
            *ee = *fe*ce/ck+g;
        }

    }
    return 0;
}