From 3c673a5c9da941e516802f0037ab7d966f8a9e8a Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Wed, 21 Feb 2024 16:48:13 -0800 Subject: [PATCH] Update ComplexElliptic.jl --- src/ComplexElliptic.jl | 48 +++++++++++++++++++----------------------- 1 file changed, 22 insertions(+), 26 deletions(-) diff --git a/src/ComplexElliptic.jl b/src/ComplexElliptic.jl index 0eeaf09..a25e10f 100644 --- a/src/ComplexElliptic.jl +++ b/src/ComplexElliptic.jl @@ -5,7 +5,6 @@ module ComplexElliptic Author: Tina Torabi Year: 2023 """ - function ellipkkp(L) # Complete elliptic integral of the first kind, with complement. if L > 10 @@ -26,7 +25,7 @@ function ellipkkp(L) c1 = (a0 - b0) / 2 i1 += 1 w1 = 2^i1 * c1^2 - mm = maximum(w1) + mm = maximum([w1]) # Ensure w1 is in an array for maximum to work correctly s0 += w1 a0 = a1 b0 = b1 @@ -37,30 +36,27 @@ function ellipkkp(L) K = Inf end - if nargout() > 1 - a0 = 1.0 - b0 = sqrt(m) - a1=0 - s0 = 1 - m - i1 = 0 - mm = 1.0 - while mm > eps() - a1 = (a0 + b0) / 2 - b1 = sqrt(a0 * b0) - c1 = (a0 - b0) / 2 - i1 += 1 - w1 = 2^i1 * c1^2 - mm = maximum(w1) - s0 += w1 - a0 = a1 - b0 = b1 - end - Kp = π / (2 * a1) - if m == 0 - Kp = Inf - end - else - Kp = nothing + # Recalculate for Kp using the complement + a0 = 1.0 + b0 = sqrt(m) + a1 = 0 + s0 = 1 - m + i1 = 0 + mm = 1.0 + while mm > eps() + a1 = (a0 + b0) / 2 + b1 = sqrt(a0 * b0) + c1 = (a0 - b0) / 2 + i1 += 1 + w1 = 2^i1 * c1^2 + mm = maximum([w1]) # Ensure w1 is in an array for maximum to work correctly + s0 += w1 + a0 = a1 + b0 = b1 + end + Kp = π / (2 * a1) + if m == 0 + Kp = Inf end return K, Kp