#include "Elliptic_integrals.h"

/**
 * Purpose: Compute elliptic integrals F(phi, k) and E(phi, k). Take from Shanjie Zhang, Jian-Ming Jin, "Computation of Special Functions", page 664-665
 * 
 * @param k - modulus k (-1 <= k <= 1)
 * @param phi - argument (in degrees), 0 <= phi <= 90
 * @return Vector with values of F(phi, k)) (the first place) and E(phi, k) (the second place)
 * @throws "Error! k should be from -1 to 1", "Error! phi should be from 0 to 90" 
 */

std::vector<double> Elliptic(const double phi, const double k)
{
    std::vector<double> result;

    if ((phi > 90) || (phi < 0))
    {
        std::cerr << "Parameter phi must be in the range [0, 90]" << std::endl;
        return result;
    }


    if ((k > 1) || (k < -1))
    {
        std::cerr << "Parameter k must be in the range [-1, 1]" << std::endl;
        return result;
    }

    double G = 0.0;
    double pi = M_PI; 
    double A0 = 1.0;
    double B0 = sqrt(1.0 - k * k);
    double D0 = pi / 180.0 * phi;
    double R = k * k;

    if (k == 1 and phi == 90)
    {
        double FE = std::numeric_limits<double>::infinity();
        double EE = 1.0;
        result.push_back(FE);
        result.push_back(EE);
    }
    else if (k == 1)
    {
        double FE = log((1.0 + sin(D0)) / cos(D0));
        double EE = sin(D0);
        result.push_back(FE);
        result.push_back(EE);
    }
    else
    {
        double FAC = 1.0;
        double A;
        double B;
        double C;
        double D;
        for (int k = 1; k <= 1000; k++)
        {
            A = (A0 + B0) / 2.0;
            B = sqrt(A0 * B0);
            C = (A0 - B0) / 2.0;
            FAC = 2.0 * FAC;
            R = R + FAC * C * C;
            if (phi != 90)
            {
                D = D0 + atan((B0 / A0) * tan(D0));
                G = G + C * sin(D);
                D0 = D + pi * int(D / pi + 0.5);
            }
            A0 = A;
            B0 = B;
            if (C < 1e-7) goto stop;
            continue;
        }
        stop:
        double CK = pi / (2.0 * A);
        double CE = pi * (2.0 - R) / (4.0 * A);
        if (phi == 90)
        {
            double FE = CK;
            double EE = CE;
            result.push_back(FE);
            result.push_back(EE);
        }
        else
        {
            double FE = D / (FAC * A);
            double EE = FE * CE / CK + G;
            result.push_back(FE);
            result.push_back(EE);
        }
    }
    return result;
}