Only use approx for x > 1.0 in elliptic integrals

ProjectTorreyPines · Dec 12, 2024 · adfe191 · adfe191
1 parent 957514c
commit adfe191
Showing 1 changed file with 18 additions and 6 deletions.
diff --git 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