view release on metacpan or  search on metacpan

src/boost/math/special_functions/ellint_3.hpp  view on Meta::CPAN

#include <boost/math/special_functions/ellint_rf.hpp>
#include <boost/math/special_functions/ellint_rj.hpp>
#include <boost/math/special_functions/ellint_1.hpp>
#include <boost/math/special_functions/ellint_2.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/tools/workaround.hpp>

// Elliptic integrals (complete and incomplete) of the third kind
// Carlson, Numerische Mathematik, vol 33, 1 (1979)

namespace boost { namespace math { 
   
namespace detail{

template <typename T, typename Policy>
T ellint_pi_imp(T v, T k, T vc, const Policy& pol);

// Elliptic integral (Legendre form) of the third kind
template <typename T, typename Policy>
T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
{
    // Note vc = 1-v presumably without cancellation error.
    T value, x, y, z, p, t;

    BOOST_MATH_STD_USING
    using namespace boost::math::tools;
    using namespace boost::math::constants;

    static const char* function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)";

    if (abs(k) > 1)
    {
       return policies::raise_domain_error<T>(function,
            "Got k = %1%, function requires |k| <= 1", k, pol);
    }

    T sphi = sin(fabs(phi));

    if(v > 1 / (sphi * sphi))
    {
        // Complex result is a domain error:
       return policies::raise_domain_error<T>(function,
            "Got v = %1%, but result is complex for v > 1 / sin^2(phi)", v, pol);
    }

    // Special cases first:
    if(v == 0)
    {
       // A&S 17.7.18 & 19
       return (k == 0) ? phi : ellint_f_imp(phi, k, pol);
    }
    if(phi == constants::pi<T>() / 2)
    {
       // Have to filter this case out before the next
       // special case, otherwise we might get an infinity from
       // tan(phi).
       // Also note that since we can't represent PI/2 exactly
       // in a T, this is a bit of a guess as to the users true
       // intent...
       //
       return ellint_pi_imp(v, k, vc, pol);
    }
    if(k == 0)
    {
       // A&S 17.7.20:
       if(v < 1)
       {
          T vcr = sqrt(vc);
          return atan(vcr * tan(phi)) / vcr;
       }
       else if(v == 1)
       {
          return tan(phi);
       }
       else
       {
          // v > 1:
          T vcr = sqrt(-vc);
          T arg = vcr * tan(phi);
          return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr);
       }
    }

    if(v < 0)
    {
       //
       // If we don't shift to 0 <= v <= 1 we get
       // cancellation errors later on.  Use
       // A&S 17.7.15/16 to shift to v > 0:
       //
       T k2 = k * k;
       T N = (k2 - v) / (1 - v);
       T Nm1 = (1 - k2) / (1 - v);
       T p2 = sqrt(-v * (k2 - v) / (1 - v));
       T delta = sqrt(1 - k2 * sphi * sphi);
       T result = ellint_pi_imp(N, phi, k, Nm1, pol);

       result *= sqrt(Nm1 * (1 - k2 / N));
       result += ellint_f_imp(phi, k, pol) * k2 / p2;
       result += atan((p2/2) * sin(2 * phi) / delta);
       result /= sqrt((1 - v) * (1 - k2 / v));
       return result;
    }
#if 0  // disabled but retained for future reference: see below.
    if(v > 1)
    {
       //
       // If v > 1 we can use the identity in A&S 17.7.7/8
       // to shift to 0 <= v <= 1.  Unfortunately this
       // identity appears only to function correctly when
       // 0 <= phi <= pi/2, but it's when phi is outside that
       // range that we really need it: That's when
       // Carlson's formula fails, and what's more the periodicity
       // reduction used below on phi doesn't work when v > 1.
       //
       // So we're stuck... the code is archived here in case
       // some bright spart can figure out the fix.
       //