Compute Jacobi elliptic functions for a complex parameter by using `sicpy.special.ellipj`. This is Julia implementation of [ELLIPJI by Moiseev Igor](https://www.mathworks.com/matlabcentral/fileexchange/17747-jacobi-elliptic-functions-sn-cn-and-dn-of-complex-phase).

ellipj.jl

using PyCall

special = pyimport("scipy.special")

ellipj(u::Float64, m::Float64) = special.ellipj(u, m)[1:3]


"""
    ellipj(u::ComplexF64, m::Float64)

Compute Jacobi elliptic functions for complex u by using `sicpy.special.ellipj`.
This is Julia implementation of [ELLIPJI by Moiseev Igor](https://www.mathworks.com/matlabcentral/fileexchange/17747-jacobi-elliptic-functions-sn-cn-and-dn-of-complex-phase).
"""
function ellipj(u::ComplexF64, m::Float64)
    if !(isreal(m))
        error("m must be real.")
    end

    if m < 0 || m > 1
        error("m must be in range 0 <= m <= 1.")
    end

    ϕ = real(u)
    ψ = imag(u)

    s, c, d = ellipj(ϕ, m)
    s1, c1, d1 = ellipj(ψ, 1-m)

    δ = c1^2 + m*s^2*s1^2
    sn = (s*d1 + im*c*d*s1*c1) / δ
    cn = (c*c1 - im*s*d*s1*d1) / δ
    dn = (d*c1*d1 - im*m*s*c*s1) / δ

    return sn, cn, dn
end