From adfe191ac2c396c1ef85d7f7ded9c11b8647b9e0 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Wed, 11 Dec 2024 17:25:12 -0800 Subject: [PATCH] Only use approx for x > 1.0 in elliptic integrals --- src/elliptic.jl | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/elliptic.jl b/src/elliptic.jl index 01f09a8..9b0768c 100644 --- a/src/elliptic.jl +++ b/src/elliptic.jl @@ -89,11 +89,15 @@ As suggested in this paper, the domain is restricted to ``(-\infty,1]``. if x == 0.0 return pi / 2 - elseif x ≈ 1.0 + elseif x == 1.0 return Inf elseif x > 1.0 - throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]")) + if x ≈ 1.0 + return Inf + else + throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]")) + end elseif 0.0 <= x < 0.1 #Table 2 from paper t = x - 0.05 @@ -265,11 +269,15 @@ function ellipe(m::Real) if x == 0.0 return pi / 2 - elseif x ≈ 1.0 + elseif x == 1.0 return 1.0 elseif x > 1.0 - throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]")) + if x ≈ 1.0 + return 1.0 + else + throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]")) + end elseif 0.0 <= x < 0.1 #Table 2 from paper t = x - 0.05 @@ -403,11 +411,15 @@ end if x == 0.0 return km, pi / 2 - elseif x ≈ 1.0 + elseif x == 1.0 return km, 1.0 elseif x > 1.0 - throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]")) + if x ≈ 1.0 + return km, 1.0 + else + throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]")) + end elseif 0.0 <= x < 0.1 #Table 2 from paper t = x - 0.05