From 01ddf31eadc8acb17e430461aa61afb8dfc57119 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Wed, 21 Feb 2024 18:37:26 -0800 Subject: [PATCH] Update ComplexElliptic.jl --- src/ComplexElliptic.jl | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/ComplexElliptic.jl b/src/ComplexElliptic.jl index a25e10f..37b46f4 100644 --- a/src/ComplexElliptic.jl +++ b/src/ComplexElliptic.jl @@ -5,64 +5,63 @@ module ComplexElliptic Author: Tina Torabi Year: 2023 """ -function ellipkkp(L) - # Complete elliptic integral of the first kind, with complement. - if L > 10 + + + function ellipkkpp(L) + if abs(L) > 10 K = π / 2 Kp = π * L + log(4) return K, Kp end - m = exp(-2 * π * L) + m = exp(-2 * π * L) a0 = 1.0 b0 = sqrt(1 - m) + a1=0 s0 = m i1 = 0 mm = 1.0 - while mm > eps() + while abs(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 + mm = w1[argmax(abs.(w1))] s0 += w1 a0 = a1 b0 = b1 end - K = π / (2 * a1) - - if m == 1 - K = Inf - end + K_ = π / (2 * a1) - # 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() + while abs(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 + mm = w1[argmax(abs.(w1))] s0 += w1 a0 = a1 b0 = b1 end Kp = π / (2 * a1) + if m == 0 Kp = Inf end - return K, Kp + return K_, Kp end + function polyval(p::Vector{T}, x::T) where T <: Real result = zero(T) for i in 1:length(p)