Skip to content

Commit

Permalink
Update ComplexElliptic.jl
Browse files Browse the repository at this point in the history
  • Loading branch information
tinatorabi authored Feb 22, 2024
1 parent 78a8780 commit 01ddf31
Showing 1 changed file with 14 additions and 15 deletions.
29 changes: 14 additions & 15 deletions src/ComplexElliptic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

2 comments on commit 01ddf31

@tinatorabi
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Error while trying to register: "Tag with name v1.0.3 already exists and points to a different commit"

Please sign in to comment.