From e81ce17bac4c16d010af711a18ed706b96cff912 Mon Sep 17 00:00:00 2001 From: Shuhua Date: Sun, 30 Aug 2020 16:09:16 +0800 Subject: [PATCH 1/3] Make ellipj and ellipke more MATLAB-like --- src/Elliptic.jl | 42 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/src/Elliptic.jl b/src/Elliptic.jl index f1ff1e7..d481790 100644 --- a/src/Elliptic.jl +++ b/src/Elliptic.jl @@ -36,11 +36,14 @@ function _E(sinphi::Float64, m::Float64) end E(phi::Real, m::Real) = E(Float64(phi), Float64(m)) - """ -`ellipke(m::Real)` -returns `(K(m), E(m))` for scalar `0 ≤ m ≤ 1` + ellipke(M) + +Complete elliptic integrals of first and second kind. Each element of `M`, say `m`, should satisfy `m ≤ 1.0`. +A tuple `(K, E)` is returned. If `M` is an array, then `K` and `E` share the same size with `M`. """ +function ellipke end + function ellipke(m::Float64) if m < 1. y = 1. - m @@ -58,6 +61,15 @@ function ellipke(m::Float64) end ellipke(x::Real) = ellipke(Float64(x)) +function ellipke(M::AbstractArray{T}) where T <: Real + K = similar(M, Float64) + E = similar(M, Float64) + for i in eachindex(M) + @inbounds K[i], E[i] = ellipke(M[i]) + end + return K, E +end + E(m::Float64) = ellipke(m)[2] E(x::Float32) = Float32(E(Float64(x))) E(x::Real) = E(Float64(x)) @@ -122,6 +134,18 @@ end Pi(n::Real, phi::Real, m::Real) = Pi(Float64(n), Float64(phi), Float64(m)) Π = Pi +""" + ellipj(U, M) + +Compute the Jacobi elliptic functions `SN`, `CN`, and `DN` evaluated for corresponding +elements of argument `U` and parameter `M`. +Note that `M` can only take values in range `[0, 1]` (both ends closed). +If both `U` and `M` are arrays, they must be the same size, and each element of the returned tuple `(SN, CN, DN)` +also has the same size. If either `U` or `M` is an array, then the other (a scalar) is broadcasted, and +each element of the returned tuple `(SN, CN, DN)` has the same size as the array. +""" +function ellipj end + function ellipj(u::Float64, m::Float64, tol::Float64) phi = Jacobi.am(u, m, tol) s = sin(phi) @@ -132,4 +156,16 @@ end ellipj(u::Float64, m::Float64) = ellipj(u, m, eps(Float64)) ellipj(u::Real, m::Real) = ellipj(Float64(u), Float64(m)) +function ellipj(U::AbstractArray{TU}, M::AbstractArray{TM}) where {TU <: Real,TM <: Real} + @assert size(U) == size(M) "U and M must be the same size" + SN = similar(U, Float64) + CN = similar(U, Float64) + DN = similar(U, Float64) + for i in eachindex(U) + @inbounds SN[i], CN[i], DN[i] = ellipj(U[i], M[i]) + end + return SN, CN, DN +end +ellipj(U::AbstractArray{TU}, m::Real) where TU <: Real = ellipj(U, fill(Float64(m), size(U))) +ellipj(u::Real, M::AbstractArray{TM}) where TM <: Real = ellipj(fill(Float64(u), size(M)), M) end # module From 988c992bd97eba8028de8163f7120a03973b5ccf Mon Sep 17 00:00:00 2001 From: Shuhua Date: Sun, 30 Aug 2020 16:10:56 +0800 Subject: [PATCH 2/3] Test improved ellipj and ellipke --- test/runtests.jl | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 9fe1008..b89e95f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -717,4 +717,32 @@ end end end +@testset "MATLAB compatibility" begin + # compare the results with the MATLAB results + Base.isapprox(x::Tuple, y::Tuple; kws...) = isapprox(collect(x), collect(y); kws...) + @testset "ellipke" begin + @test ellipke(0.5) ≈ (1.854074677301372, 1.350643881047675) + @test ellipke([0, 0.3, 0.7, 1.0]) ≈ ([1.570796326794897, 1.713889448178791, 2.075363135292469, Inf], + [1.570796326794897, 1.445363064412665, 1.241670567945823, 1.0]) + end + @testset "ellipj" begin + @test ellipj(1.2, 0.5) ≈ (0.887715488619278, 0.460392453527896, 0.778447561260692) + S, C, D = ellipj(-1.0, [0.0, 0.7, 1.0]) + @test S ≈ [-0.841470984807897, -0.786762346175161, -0.761594155955765] + @test C ≈ [0.540302305868140, 0.617256033296522, 0.648054273663885] + @test D ≈ [1.0, 0.752797122370078, 0.648054273663885] + S, C, D = ellipj([-1.5, 1.5], 0.23) + @test S ≈ [-0.988195152460524, 0.988195152460524] + @test C ≈ [0.153200328503307, 0.153200328503307] + @test D ≈ [0.880566964148843, 0.880566964148843] + S, C, D = ellipj([-1.5 1.5; 0.1 9.9], [0.1 0.3; 0.5 0.9]) + @test S ≈ [-0.994313629062765 0.983958114459644 + 0.099750685474625 -0.391507114124818] + @test C ≈ [0.106491347348199 0.178399632816390 + 0.995012462609058 0.920175080943652] + @test D ≈ [0.949280801821044 0.842346679637729 + 0.997509348514424 0.928466456922753] + end +end + include("jacobi_tests.jl") From dbcc8ac05b88269ed62591df534029723804eeb2 Mon Sep 17 00:00:00 2001 From: Shuhua Date: Sun, 30 Aug 2020 16:13:23 +0800 Subject: [PATCH 3/3] Update README for improved ellipj and ellipke --- README.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/README.md b/README.md index b88bc9e..8e53778 100644 --- a/README.md +++ b/README.md @@ -76,6 +76,10 @@ also provided. `ellipj(u,m)` is equivalent to `sn(u,m), cn(u,m), dn(u,m)`, but faster if you want all three. Likewise, `ellipke(m)` is equivalent to `K(m), E(m)`, but faster if you want both. +Additionally, you may also input arrays of the compatible size to calculate elliptic integrals +for multiple values together. Please check the documentation of `ellipke` and `ellipj` for details. +Example usages are also available in `test/runtests.jl" (see the "MATLAB compatibility" testset). + ```jlcon julia> import Elliptic @@ -84,6 +88,17 @@ julia> k,e = Elliptic.ellipke(0.5) julia> sn,cn,dn = Elliptic.ellipj(0.672, 0.36) (0.6095196917919022,0.792770928653356,0.9307281387786907) + +julia> using Elliptic + +julia> ellipke([0, 0.3, 0.7, 1.0]) +([1.5707963267948968, 1.7138894481787914, 2.0753631352924686, Inf], + [1.5707963267948968, 1.4453630644126656, 1.2416705679458224, 1.0]) + +julia> ellipj([-1.5, 1.5], 0.23) +([-0.9881951524605244, 0.9881951524605244], + [0.15320032850330667, 0.15320032850330667], + [0.8805669641488431, 0.8805669641488431]) ``` Installation