boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-28 07:22:12 +00:00
Code Issues Projects Releases Wiki Activity

Apply Carlson's latest relations for the Legendre Elliptic integrals.

Browse Source
This commit is contained in:
jzmaddock
2015-02-16 17:49:56 +00:00
parent f7c55524be
commit 5eb3bbbfcc
2 changed files with 34 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
include/boost/math/special_functions/ellint_1.hpp
View File

@@ -103,10 +103,17 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
}
T sinp = sin(rphi);
sinp *= sinp;
T cosp = cos(rphi);
cosp *= cosp;
T c = 1 / sinp;
BOOST_MATH_INSTRUMENT_VARIABLE(sinp);
BOOST_MATH_INSTRUMENT_VARIABLE(cosp);
result = s * sinp * ellint_rf_imp(T(cosp * cosp), T(1 - k * k * sinp * sinp), T(1), pol);
//
// Use http://dlmf.nist.gov/19.25#E5, note that
// c-1 simplifies to cot^2(rphi) which avoid cancellation:
//
result = rphi == 0 ? 0 : s * ellint_rf_imp(T(cosp / sinp), c - k * k, c, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
if(m != 0)
{

32
include/boost/math/special_functions/ellint_2.hpp
View File

@@ -21,6 +21,7 @@
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/ellint_rf.hpp>
#include <boost/math/special_functions/ellint_rd.hpp>
#include <boost/math/special_functions/ellint_rg.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/tools/workaround.hpp>
@@ -87,11 +88,30 @@ T ellint_e_imp(T phi, T k, const Policy& pol)
}
T sinp = sin(rphi);
T cosp = cos(rphi);
T x = cosp * cosp;
T t = k * k * sinp * sinp;
T y = 1 - t;
T z = 1;
result = s * sinp * (ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3);
T c = 1 / (sinp * sinp);
T cm1 = cosp * cosp / (sinp * sinp); // c - 1
T k2 = k * k;
if(k2 > 1)
{
return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
}
else if(rphi == 0)
{
result = 0;
}
else if(k == 0)
{
result = rphi;
}
else if(k2 == 1)
{
result = sinp;
}
else
{
// http://dlmf.nist.gov/19.25#E10
result = s * ((1 - k2) * ellint_rf_imp(cm1, c - k2, c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, c - k2, pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2))));
}
if(m != 0)
result += m * ellint_e_imp(k, pol);
}
@@ -119,7 +139,7 @@ T ellint_e_imp(T k, const Policy& pol)
T t = k * k;
T y = 1 - t;
T z = 1;
T value = ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3;
T value = 2 * ellint_rg_imp(x, y, z, pol);
return value;
}