SHGLQ
SHExpandGLQ
MakeGridGLQ

# 4π normalized
# Orthonormalized
# Schmidt semi-normalized
# Unnormalized
# Utilities

# Equally sampled (N×N) and equally spaced (N×2N) grids

# Gauss-Legendre quadrature grids

end
export SHExpandGLQ!
function SHExpandGLQ!(cilm::AbstractArray{Cdouble,3}, lmax::Integer,
                      gridglq::AbstractArray{Cdouble,2},
                      w::AbstractVector{Cdouble},
                      plx::Optional{AbstractArray{Cdouble,2}},
                      zero::Optional{AbstractVector{Cdouble}}; norm::Integer=1,
                      csphase::Integer=1, lmax_calc::Optional{Integer}=nothing,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    lmax_calc′ = optional(lmax_calc, lmax)
    @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(gridglq) == (lmax + 1, 2 * lmax + 1)
    @assert length(w) == lmax + 1
    @assert (plx !== nothing) + (zero !== nothing) == 1
    if plx !== nothing
        @assert size(plx) == (lmax + 1, (lmax + 1) * (lmax + 2) ÷ 2)
    end
    if zero !== nothing
        @assert length(zero) == lmax + 1
    end
    @assert norm ∈ (1, 2, 3, 4)
    @assert csphase ∈ (1, -1)
    exitstatus′ = Ref{Cint}()
    ccall((:SHExpandGLQ, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}), cilm,
          size(cilm, 2), lmax, gridglq, w, optional(plx, Ptr{Cdouble}()),
          optional(zero, Ptr{Cdouble}()), norm, csphase, lmax_calc′,
          exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 && error("SHExpandGLQ!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    return cilm
end
export SHExpandGLQ
function SHExpandGLQ(lmax::Integer, gridglq::AbstractArray{Cdouble,2},
                     w::AbstractVector{Cdouble},
                     plx::Optional{AbstractArray{Cdouble,2}},
                     zero::Optional{AbstractVector{Cdouble}}; norm::Integer=1,
                     csphase::Integer=1, lmax_calc::Optional{Integer}=nothing,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    lmax_calc′ = optional(lmax_calc, lmax)
    @assert lmax_calc′ ≥ 0
    cilm = Array{Cdouble}(undef, 2, lmax_calc′ + 1, lmax_calc′ + 1)
    SHExpandGLQ!(cilm, lmax, gridglq, w, plx, zero; norm=norm, csphase=csphase,
                 lmax_calc=lmax_calc, exitstatus=exitstatus)
    return cilm
end
export MakeGridGLQ!
function MakeGridGLQ!(gridglq::AbstractArray{Cdouble,2},
                      cilm::AbstractArray{Cdouble,3}, lmax::Integer,
                      plx::Optional{AbstractArray{Cdouble,2}},
                      zero::Optional{AbstractVector{Cdouble}}; norm::Integer=1,
                      csphase::Integer=1, lmax_calc::Optional{Integer}=nothing,
                      extend::Integer=0,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    @assert extend ∈ (0, 1)
    @assert size(gridglq) == (lmax + 1, 2 * lmax + 1 + extend)
    @assert size(cilm, 1) == 2
    @assert size(cilm, 2) ≥ 1
    @assert size(cilm, 3) == size(cilm, 2)
    @assert (plx !== nothing) + (zero !== nothing) == 1
    if plx !== nothing
        @assert size(plx) == (lmax + 1, (lmax + 1) * (lmax + 2) ÷ 2)
    end
    if zero !== nothing
        @assert length(zero) == lmax + 1
    end
    @assert norm ∈ (1, 2, 3, 4)
    @assert csphase ∈ (1, -1)
    lmax_calc′ = optional(lmax_calc, lmax)
    exitstatus′ = Ref{Cint}()
    ccall((:MakeGridGLQ, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble},
           Ptr{Cdouble}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}),
          gridglq, size(gridglq, 1), size(gridglq, 2), cilm, size(cilm, 2),
          lmax, optional(plx, Ptr{Cdouble}()), optional(zero, Ptr{Cdouble}()),
          norm, csphase, lmax_calc′, extend, exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 && error("MakeGridGLQ!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    return gridglq

export MakeGridGLQ
function MakeGridGLQ(cilm::AbstractArray{Cdouble,3}, lmax::Integer,
                     plx::Optional{AbstractArray{Cdouble,2}},
                     zero::Optional{AbstractVector{Cdouble}}; norm::Integer=1,
                     csphase::Integer=1, lmax_calc::Optional{Integer}=nothing,
                     extend::Integer=0,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    @assert extend ∈ (0, 1)
    gridglq = Array{Cdouble}(undef, lmax + 1, 2 * lmax + 1 + extend)
    MakeGridGLQ!(gridglq, cilm, lmax, plx, zero; norm=norm, csphase=csphase,
                 lmax_calc=lmax_calc, extend=extend, exitstatus=exitstatus)
    return gridglq

end

    griddh′′, n′′ = MakeGridDH(cilm′, lmax′)
    @test n′′ == n
    @test griddh′′ == griddh′

    griddh″, n″ = MakeGridDH(cilm′, lmax′)
    @test n″ == n
    @test griddh″ == griddh′

    griddh′′, n′′ = MakeGridDHC(cilm′, lmax′)
    @test n′′ == n
    @test griddh′′ == griddh′

    griddh″, n″ = MakeGridDHC(cilm′, lmax′)
    @test n″ == n
    @test griddh″ == griddh′

@testset "SHGLQ" begin
    zero, w = SHGLQ(nothing, 4)

@testset "SHGLQ and SHExpandGLQ" begin
    lmax = 4
    zero, w = SHGLQ(nothing, lmax)
    # Invent a random grid
    gridglq = randn(lmax + 1, 2 * lmax + 1)
    # Calculate coefficients
    cilm = SHExpandGLQ(lmax, gridglq, w, nothing, zero)
    # Go back to the grid (it'll be low-pass filtered)
    gridglq′ = MakeGridGLQ(cilm, lmax, nothing, zero)
    @test size(gridglq′) == size(gridglq)
    # Go back to coefficients again (they must be the same)
    cilm′ = SHExpandGLQ(lmax, gridglq′, w, nothing, zero)
    @test isapprox(cilm′, cilm)
    # Go back to the grid again (it must the be same this time)
    gridglq″ = MakeGridGLQ(cilm′, lmax, nothing, zero)
    @test size(gridglq″) == size(gridglq′)
    @test isapprox(gridglq″, gridglq′)