eschnett/SHTOOLS.jl - Change 24V7KTFI76AULBAV4MVCC2HMJLLNMWEVB5NN7KT7V55PLW565QFQC

Implement some transforms

Created by Erik Schnetter on February 17, 2021

24V7KTFI76AULBAV4MVCC2HMJLLNMWEVB5NN7KT7V55PLW565QFQC

Dependencies

In channels

main

Change contents

Insertion in README.md at line 10 [3.3301]

[3.3678]

[4.0]


Note that the function arguments are not always the same. For example,
output arguments to mutating functions come first in the argument
list, and are omitted in non-mutating functions.

Replacement in README.md at line 17 [3.3301]
B:BD[4.35] → [2.0:19]
```
Legendre Functions
```
[4.35]
[4.35]
```
### Legendre Functions:
```

Insertion in README.md at line 42 [3.3301]

[4.217]


### Spherical harmonic transforms
```Julia
SHExpandDH
MakeGridDH
```

Replacement in src/SHTOOLS.jl at line 5 [3.421]
B:BD[3.457] → [3.457:491]
```
const Maybe{T} = Union{Nothing,T}
```
[3.457]
[2.20]
```
const Optional{T} = Union{Nothing,T}
```

Replacement in src/SHTOOLS.jl at line 7 [3.421]

B:BD[2.45] → [2.45:125]

cmaybe(::Type{T}, x::Maybe{T}) where {T} = x === nothing ? Ptr{T}() : Ref{T}(x)

[2.45]

[3.680]

optional(x::Optional, x0) = x === nothing ? x0 : x
function optional(::Type{T}, x::Optional{T}) where {T}
    return x === nothing ? Ptr{T}() : Ref{T}(x)
end

Insertion in src/SHTOOLS.jl at line 13 [3.421]
[3.762]
[3.762]
```
# Legendre functions
```

Replacement in src/SHTOOLS.jl at line 16 [3.421]

B:BD[3.763] → [3.763:817]

for op in [:PlmBar, :PlmON, :PlmSchmidt, :PLegendreA]

[3.763]

[2.126]

for op in [:PlmBar, :PlmON, :PlmSchmidt]

Replacement in src/SHTOOLS.jl at line 24 [3.421]

B:BD[2.271] → [2.271:499]

                      z::Union{AbstractFloat,Integer};
                      csphase::Maybe{Integer}=nothing,
                      cnorm::Maybe{Integer}=nothing,
                      exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.271]

[2.499]

                      z::Union{AbstractFloat,Integer}; csphase::Integer=1,
                      cnorm::Integer=1,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert cnorm ∈ (0, 1)

Replacement in src/SHTOOLS.jl at line 31 [3.421]

B:BD[2.562] → [2.562:592]

            exitstatus′ = 0

[2.562]

[3.1170]

            exitstatus′ = Ref{Cint}()

Replacement in src/SHTOOLS.jl at line 34 [3.421]

B:BD[3.1297] → [3.1297:1363]

B:BD[3.1363] → [2.593:647]

                   Ref{Cint}), p, lmax, z, cmaybe(Cint, csphase),
                  cmaybe(Cint, cnorm), exitstatus′)

[3.1297]

[2.647]

                   Ref{Cint}), p, lmax, z, csphase, cnorm, exitstatus′)

Replacement in src/SHTOOLS.jl at line 36 [3.421]

B:BD[2.685] → [2.685:768]

                exitstatus′ ≠ 0 && error("$($op!): Error code $exitstatus′")

[2.685]

[2.768]

                exitstatus′[] ≠ 0 &&
                    error("$($op!): Error code $(exitstatus′[])")

Replacement in src/SHTOOLS.jl at line 39 [3.421]

B:BD[2.785] → [2.785:830]

                exitstatus[] = exitstatus′

[2.785]

[2.830]

                exitstatus[] = exitstatus′[]

Replacement in src/SHTOOLS.jl at line 46 [3.421]

B:BD[2.935] → [2.935:1105]

                     csphase::Maybe{Integer}=nothing,
                     cnorm::Maybe{Integer}=nothing,
                     exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.935]

[3.1715]

                     csphase::Integer=1, cnorm::Integer=1,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0

Replacement in src/SHTOOLS.jl at line 58 [3.421]

B:BD[2.1393] → [2.1393:1634]

                         z::Union{AbstractVector,Integer};
                         csphase::Maybe{Integer}=nothing,
                         cnorm::Maybe{Integer}=nothing,
                         exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.1393]

[2.1634]

                         z::Union{AbstractVector,Integer}; csphase::Integer=1,
                         cnorm::Integer=1,
                         exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert cnorm ∈ (0, 1)

Replacement in src/SHTOOLS.jl at line 66 [3.421]
B:BD[2.1721] → [2.1721:1751]
```
            exitstatus′ = 0
```
[2.1721]
[3.1852]
```
            exitstatus′ = Ref{Cint}()
```

Replacement in src/SHTOOLS.jl at line 69 [3.421]

B:BD[3.1985] → [3.1985:2066]

B:BD[3.2066] → [2.1752:1806]

                   Ref{Cint}, Ref{Cint}), p, dp, lmax, z, cmaybe(Cint, csphase),
                  cmaybe(Cint, cnorm), exitstatus′)

[3.1985]

[2.1806]

                   Ref{Cint}, Ref{Cint}), p, dp, lmax, z, csphase, cnorm,
                  exitstatus′)

Replacement in src/SHTOOLS.jl at line 72 [3.421]

B:BD[2.1844] → [2.1844:1930]

                exitstatus′ ≠ 0 && error("$($op_d1!): Error code $exitstatus′")

[2.1844]

[2.1930]

                exitstatus′[] ≠ 0 &&
                    error("$($op_d1!): Error code $(exitstatus′[])")

Replacement in src/SHTOOLS.jl at line 75 [3.421]

B:BD[2.1947] → [2.1947:1992]

                exitstatus[] = exitstatus′

[2.1947]

[2.1992]

                exitstatus[] = exitstatus′[]

Replacement in src/SHTOOLS.jl at line 82 [3.421]

B:BD[2.2149] → [2.2149:2328]

                        csphase::Maybe{Integer}=nothing,
                        cnorm::Maybe{Integer}=nothing,
                        exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.2149]

[2.2328]

                        csphase::Integer=1, cnorm::Integer=1,
                        exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0

Insertion in src/SHTOOLS.jl at line 89 [3.421]

[2.2574]

[3.2134]

            return p, dp
        end
        #
    end
end
for op in [:PLegendreA]
    op! = Symbol(op, "!")
    op_d1 = Symbol(op, "_d1")
    op_d1! = Symbol(op, "_d1!")
    @eval begin
        #
        export $op!
        function $op!(p::AbstractVector{Cdouble}, lmax::Integer,
                      z::Union{AbstractFloat,Integer}; csphase::Integer=1,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert length(p) ≥ (lmax + 1) * (lmax + 2) ÷ 2
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Cint, Cdouble, Ref{Cint}, Ref{Cint}), p, lmax,
                  z, csphase, exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p
        end
        #
        export $op
        function $op(lmax::Integer, z::Union{AbstractFloat,Integer};
                     csphase::Integer=1,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            $op!(p, lmax, z; csphase=csphase, exitstatus=exitstatus)
            return p
        end
        #
        export $op_d1!
        function $op_d1!(p::AbstractVector{Cdouble},
                         dp::AbstractVector{Cdouble}, lmax::Integer,
                         z::Union{AbstractVector,Integer}; csphase::Integer=1,
                         exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert length(p) ≥ lmax + 1
            @assert length(dp) ≥ lmax + 1
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op_d1)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble, Ref{Cint},
                   Ref{Cint}), p, dp, lmax, z, csphase, exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op_d1!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p, dp
        end
        #
        export $op_d1
        function $op_d1(lmax::Integer, z::Union{AbstractFloat,Integer};
                        csphase::Integer=1,
                        exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            dp = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            $op_d1!(p, dp, lmax, z; csphase=csphase, exitstatus=exitstatus)

Replacement in src/SHTOOLS.jl at line 176 [3.421]

B:BD[2.2775] → [2.2775:2840]

                      exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.2775]

[2.2840]

                      exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0

Replacement in src/SHTOOLS.jl at line 179 [3.421]
B:BD[2.2883] → [2.2883:2913]
```
            exitstatus′ = 0
```
[2.2883]
[3.2517]
```
            exitstatus′ = Ref{Cint}()
```

Replacement in src/SHTOOLS.jl at line 184 [3.421]

B:BD[2.2985] → [2.2985:3068]

                exitstatus′ ≠ 0 && error("$($op!): Error code $exitstatus′")

[2.2985]

[2.3068]

                exitstatus′[] ≠ 0 &&
                    error("$($op!): Error code $(exitstatus′[])")

Replacement in src/SHTOOLS.jl at line 187 [3.421]

B:BD[2.3085] → [2.3085:3130]

                exitstatus[] = exitstatus′

[2.3085]

[2.3130]

                exitstatus[] = exitstatus′[]

Replacement in src/SHTOOLS.jl at line 194 [3.421]

B:BD[2.3235] → [2.3235:3299]

                     exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.3235]

[2.3299]

                     exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0

Replacement in src/SHTOOLS.jl at line 205 [3.421]

B:BD[2.3645] → [2.3645:3713]

                         exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.3645]

[2.3713]

                         exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0

Replacement in src/SHTOOLS.jl at line 209 [3.421]
B:BD[2.3800] → [2.3800:3830]
```
            exitstatus′ = 0
```
[2.3800]
[3.3028]
```
            exitstatus′ = Ref{Cint}()
```

Replacement in src/SHTOOLS.jl at line 214 [3.421]

B:BD[2.3911] → [2.3911:3997]

                exitstatus′ ≠ 0 && error("$($op_d1!): Error code $exitstatus′")

[2.3911]

[2.3997]

                exitstatus′[] ≠ 0 &&
                    error("$($op_d1!): Error code $(exitstatus′[])")

Replacement in src/SHTOOLS.jl at line 217 [3.421]

B:BD[2.4014] → [2.4014:4059]

                exitstatus[] = exitstatus′

[2.4014]

[2.4059]

                exitstatus[] = exitstatus′[]

Replacement in src/SHTOOLS.jl at line 224 [3.421]

B:BD[2.4216] → [2.4216:4283]

                        exitstatus::Maybe{Ref{<:Integer}}=nothing)

[2.4216]

[2.4283]

                        exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0

Insertion in src/SHTOOLS.jl at line 242 [3.421]

[3.3288]

[5.776]


################################################################################
# Spherical harmonic transforms

Insertion in src/SHTOOLS.jl at line 247 [3.421]

[5.777]

export SHExpandDH!
function SHExpandDH!(cilm::Array{Cdouble,3}, lmax::Optional{Ref{<:Integer}},
                     griddh::Array{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((:SHExpandDH, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Cint, Ptr{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("SHExpandDH!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    if lmax !== nothing
        lmax[] = lmax′[]
    end
    return cilm, Int(lmax′[])

Insertion in src/SHTOOLS.jl at line 281 [3.421]

[5.781]


export SHExpandDH
function SHExpandDH(griddh::Array{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{Cdouble}(undef, 2, lmax_calc′ + 1, lmax_calc′ + 1)
    lmax = Ref{Int}()
    SHExpandDH!(cilm, lmax, griddh, n; norm=norm, sampling=sampling,
                csphase=csphase, lmax_calc=lmax_calc, exitstatus=exitstatus)
    return cilm, lmax[]
end
export MakeGridDH!
function MakeGridDH!(griddh::Array{Cdouble,2}, n::Optional{Ref{<:Integer}},
                     cilm::Array{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((:MakeGridDH, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ref{Cint}, Ptr{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("MakeGridDH!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    if n !== nothing
        n[] = n′
    end
    return griddh, n′
end
export MakeGridDH
function MakeGridDH(cilm::Array{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{Cdouble}(undef, n′ + extend, sampling * n′ + extend)
    _, n = MakeGridDH!(griddh, nothing, cilm, lmax; norm=norm,
                       sampling=sampling, csphase=csphase, lmax_calc=lmax_calc,
                       extend=extend, exitstatus=exitstatus)
    return griddh, n
end
end

Insertion in test/runtests.jl at line 3 [3.22]

[3.48]

[2.4506]


################################################################################
# Legendre Polynomials

Insertion in test/runtests.jl at line 107 [3.22]

[2.6445]

end
################################################################################
# Spherical harmonic transforms
@testset "MakeGridDH and SHExpandDH" begin
    n = 10
    griddh = randn(n, n)
    cilm, lmax = SHExpandDH(griddh, n)
    @test lmax == 4
    griddh′, n′ = MakeGridDH(cilm, lmax)
    @test n′ == n
    
    cilm′, lmax′ = SHExpandDH(griddh, n)
    @test lmax′ == lmax
    @test cilm′ == cilm
    griddh′′, n′′ = MakeGridDH(cilm′, lmax′)
    @test n′′ == n
    @test griddh′′ == griddh′