SHExpandDHC
MakeGridDHC

export SHExpandDHC!
function SHExpandDHC!(cilm::Array{Complex{Cdouble},3},
                      lmax::Optional{Ref{<:Integer}},
                      griddh::Array{Complex{Cdouble},2}, n::Integer;
                      norm::Integer=1, sampling::Integer=1, csphase::Integer=1,
                      lmax_calc::Optional{Integer}=nothing,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert n > 0 && iseven(n)
    @assert norm ∈ (1, 2, 3, 4)
    @assert sampling ∈ (1, 2)
    @assert csphase ∈ (1, -1)
    lmax_calc′ = optional(lmax_calc, n ÷ 2 - 1)
    @assert lmax_calc′ ≥ 0
    @assert size(cilm, 1) == 2
    @assert size(cilm, 2) ≥ lmax_calc′ + 1
    @assert size(cilm, 3) == size(cilm, 2)
    @assert size(griddh, 1) ≥ n
    @assert size(griddh, 2) ≥ sampling * n
    lmax′ = Ref{Cint}()
    exitstatus′ = Ref{Cint}()
    ccall((:SHExpandDHC, libSHTOOLS), Cvoid,
          (Ptr{Complex{Cdouble}}, Cint, Cint, Cint, Ptr{Complex{Cdouble}}, Cint,
           Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}),
          griddh, size(griddh, 1), size(griddh, 2), n, cilm, size(cilm, 2),
          lmax′, norm, sampling, csphase, lmax_calc′, exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 && error("SHExpandDHC!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    if lmax !== nothing
        lmax[] = lmax′[]
    end
    return cilm, Int(lmax′[])
end
export SHExpandDHC
function SHExpandDHC(griddh::Array{Complex{Cdouble},2}, n::Integer;
                     norm::Integer=1, sampling::Integer=1, csphase::Integer=1,
                     lmax_calc::Optional{Integer}=nothing,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert n > 0 && iseven(n)
    lmax_calc′ = optional(lmax_calc, n ÷ 2 - 1)
    cilm = Array{Complex{Cdouble}}(undef, 2, lmax_calc′ + 1, lmax_calc′ + 1)
    lmax = Ref{Int}()
    SHExpandDHC!(cilm, lmax, griddh, n; norm=norm, sampling=sampling,
                 csphase=csphase, lmax_calc=lmax_calc, exitstatus=exitstatus)
    return cilm, lmax[]

export MakeGridDHC!
function MakeGridDHC!(griddh::Array{Complex{Cdouble},2},
                      n::Optional{Ref{<:Integer}},
                      cilm::Array{Complex{Cdouble},3}, lmax::Integer;
                      norm::Integer=1, sampling::Integer=1, csphase::Integer=1,
                      lmax_calc::Optional{Integer}=nothing, extend::Integer=0,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    @assert norm ∈ (1, 2, 3, 4)
    @assert sampling ∈ (1, 2)
    @assert csphase ∈ (1, -1)
    @assert extend ∈ (0, 1)
    n′ = 2 * lmax + 2
    @assert size(griddh, 1) ≥ n′ + extend
    @assert size(griddh, 2) ≥ sampling * n′ + extend
    @assert size(cilm, 1) == 2
    @assert size(cilm, 2) ≥ lmax + 1
    @assert size(cilm, 3) == size(cilm, 2)
    lmax_calc′ = optional(lmax_calc, lmax)
    exitstatus′ = Ref{Cint}()
    ccall((:MakeGridDHC, libSHTOOLS), Cvoid,
          (Ptr{Complex{Cdouble}}, Cint, Cint, Ref{Cint}, Ptr{Complex{Cdouble}},
           Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint},
           Ref{Cint}), griddh, size(griddh, 1), size(griddh, 2), n′, cilm,
          size(cilm, 2), lmax, norm, sampling, csphase, lmax_calc′, extend,
          exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 && error("MakeGridDHC!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    if n !== nothing
        n[] = n′
    end
    return griddh, n′
end
export MakeGridDHC
function MakeGridDHC(cilm::Array{Complex{Cdouble},3}, lmax::Integer;
                     norm::Integer=1, sampling::Integer=1, csphase::Integer=1,
                     lmax_calc::Optional{Integer}=nothing, extend::Integer=0,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    @assert sampling ∈ (1, 2)
    @assert extend ∈ (0, 1)
    n′ = 2 * lmax + 2
    griddh = Array{Complex{Cdouble}}(undef, n′ + extend, sampling * n′ + extend)
    _, n = MakeGridDHC!(griddh, nothing, cilm, lmax; norm=norm,
                        sampling=sampling, csphase=csphase, lmax_calc=lmax_calc,
                        extend=extend, exitstatus=exitstatus)
    return griddh, n
end
# TODO: Add MakeGradientDH once there is a C wrapper
end

    @test n′′ == n
    @test griddh′′ == griddh′
end
@testset "MakeGridDHC and SHExpandDHC" begin
    # Invent a random grid
    n = 10
    griddh = randn(Complex{Float64}, n, n)
    # Calculate coefficients
    cilm, lmax = SHExpandDHC(griddh, n)
    @test lmax == 4
    # Go back to the grid (it'll be low-pass filtered)
    griddh′, n′ = MakeGridDHC(cilm, lmax)
    @test n′ == n
    # Go back to coefficients again (they must be the same)
    cilm′, lmax′ = SHExpandDHC(griddh, n)
    @test lmax′ == lmax
    @test cilm′ == cilm
    # Go back to the grid again (it must the be same this time)
    griddh′′, n′′ = MakeGridDHC(cilm′, lmax′)