Algorythms/modules/Elliptic_integrals.cpp

#include "Elliptic_integrals.h"

// std::string Error()
// {
//     std::string mes = "Error! k should be from -1 to 1";
//     std::exception(mes);
// }

/**
 * Purpose: Constructor
 * 
 * @return None
 * @throws None
 */

Elliptic_Integrals::Elliptic_Integrals()
{

}

/**
 * Purpose: Destructor
 * 
 * @return None
 * @throws None
 */

Elliptic_Integrals::~Elliptic_Integrals()
{
    
}

/**
 * Purpose: Compute complete elliptic integrals K(k) and E(k)
 * 
 * @param k - modulus k (-1 <= k <= 1)
 * @return Vector with values of K(k) (the first place) and E(k) (the second place)
 * @throws "Error! k should be from -1 to 1"
 */

std::vector<double> Elliptic_Integrals::CompleteElliptic(const double k)
{
    std::vector<double> result;

    double PK = 1.0 - k * k;

    if ((k > 1) || (k < -1))
    {
        std::cerr << "Parameter k must be in the range [-1, 1]" << std::endl;
        return result;
    }
    
    if (k == 1.0)
    {
        double CK = std::numeric_limits<double>::infinity();
        double CE = 1.0;
        result.push_back(CK);
        result.push_back(CE);
    }
    else
    {
        double AK = (((0.01451196212 * PK + 0.03742563713) * PK + 0.03590092383) * PK + 0.09666344259) * PK + 1.38629436112;
        double BK = (((0.00441787012 * PK + 0.03328355346) * PK + 0.06880248576) * PK + 0.12498593597) * PK + 0.5;
        double CK = AK - BK * log(PK);
        double AE = (((0.01736506451 * PK + 0.04757383546) * PK + 0.06260601220) * PK + 0.44325141463) * PK + 1.0;
        double BE = (((0.00526449639 * PK + 0.04069697526) * PK + 0.09200180037) * PK + 0.24998368310) * PK;
        double CE = AE - BE * log(PK);
        result.push_back(CK);
        result.push_back(CE);
    }
    return result;
}


/**
 * Purpose: Compute elliptic integrals F(phi, k) and E(phi, k)
 * 
 * @param k - modulus k (-1 <= k <= 1)
 * @param phi - argument (in degrees)
 * @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"
 */

std::vector<double> Elliptic_Integrals::Elliptic(const double phi, const double k)
{
    std::vector<double> 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;
}