// 楕円積分と Jacobi の楕円関数
#include <cmath>
using namespace std;

double ellpC(double x);
double ellpI(double x, double m);
int sncndn(double u, double k, double &sn, double &cn, double &dn);

// 第１種完全楕円積分
// 　算術幾何平均による方法を利用
//   山内他編：“電子計算機のための数値計算法Ⅲ”，p.258，培風館，1972 参照
//   1-x*x に対する値を計算
double ellpC(double x)
{
    int j;
    double a, at, b;

    a = 1.0;
	b = sqrt(x);
    for (j=0; j<20; j++)
    {
        at = (a + b)/2.0;
        b = sqrt(a*b);
        a = at;
        if ((a - b)/a < 1.0E-8) break;
    }
    if (j >= 20)
    {
        ShowMessage("ellpC: 収束しません．");
        return 0;
    }
    else
        return M_PI_2/a;
}

// 第１種不完全楕円積分
//   ただし，kc = sqrt(1-k*k)
double ellpI(double x, double kc)
{
    double a, b, y, i, j, bt;
    a = 1.0;
    b = kc;
    y = 1.0/x;
    i = 0;
    for (j=0; j<15; j++)
    {
        bt = a*b;
        a = a + b;
        b = 2.0*sqrt(bt);
        y = y - bt/y;
        if (y == 0.0) y = sqrt(bt)*1.0E-10;
        if (fabs(a - b) < a*1.0E-10) break;
        i = 2*i;
        if (y < 0.0) i++;
    }
    if (j >= 15)
    {
        ShowMessage("ellpI: 収束しません．");
        return 0;
    }
    if (y < 0.0) i++;
    return ( (atan(a/y) + M_PI*i)/a );
}

// Jacobi の楕円関数，sn(u,kc), cn(u,kc), dn(u,kc) の計算
//   ただし，kc = 1-k*k
int sncndn(double u, double kc, double &sn, double &cn, double &dn)
{
    double a, at, b, c, d, e, aa[16], bb[16];
    int i;

    if (u == 0.0)
    {
        sn = 0.0;
        cn = 1.0;
        dn = 1.0;
        return 0;
    }
    a = 1.0;
    b = sqrt(kc);
    for (i=0; i<16; i++)
    {
        aa[i] = a;
        bb[i] = b;
        at = (a + b)/2.0;
        b = sqrt(a*b);
        a = at;
        if ((a-b)/a < 1.0E-8) break;
     }
     if (i >= 16)
     {
        ShowMessage("sncndn: 収束しません．");
        return -1;
     }
     c = a/tan(u*a);
     d = 1.0;
     do
     {
        e = c*c/a;
        c = c*d;
        a = aa[i];
        d = (e + bb[i])/(e + a);
     } while (i--);
     sn = 1.0/sqrt(1.0 + c*c);
     cn = sn*c;
     dn = d;
     return 0;
}

