X5ISFPUVKYET33KOY3IN3FSI54BWY2JDPLGBAK5ENTAGBH54NK4QC
74KULQ54O2O72SWNQ7K3BNMGD2IF4ODPA7LOAIEINQUTDBVXFPLAC
DKEAEDXCBBRU2D42EQRQ6EQ7NEXJMR76VLF3RIORFMUYYPAXWJLQC
24V7KTFI76AULBAV4MVCC2HMJLLNMWEVB5NN7KT7V55PLW565QFQC
WG44JQLOE3BQOS6DZMCTCRWSNDLU2CA4DP7ZTWADEOP2GQ5TQDBQC
XEWLD2DYMEMKHCAWJZ7XEHAZTAYSTWCZMTRTEWMBPYJ4EE3NBSXAC
MEY5CG3E57E42TBUCYJ4URBEHYQV2NL645WZWXUGR37Y3R4IBUIQC
7VIMUPV7XRRR3NTHQBBAFKUGI3Y23PSCOEY4CXQ3DVFKC2ACHN5AC
end
# Other routines
export SHExpandLSQ!
function SHExpandLSQ!(cilm::AbstractArray{Cdouble,3},
d::AbstractVector{Cdouble}, lat::AbstractVector{Cdouble},
lon::AbstractVector{Cdouble}, nmax::Integer,
lmax::Integer; norm::Integer=1, csphase::Integer=1,
weights::Vector{Cdouble},
exitstatus::Optional{Ref{<:Integer}}=nothing)
@assert lmax ≥ 0
@assert size(cilm, 1) == 2
@assert size(cilm, 2) ≥ lmax + 1
@assert size(cilm, 3) == size(cilm, 2)
@assert nmax ≥ 0
@assert nmax ≥ (lmax + 1)^2
@assert length(d) ≥ nmax
@assert length(lat) ≥ nmax
@assert length(lon) ≥ nmax
@assert norm ∈ (1, 2, 3, 4)
@assert length(weights) == nmax
chi2′ = Ref{Cdouble}()
exitstatus′ = Ref{Cint}()
ccall((:SHExpandLSQ, libSHTOOLS), Cvoid,
(Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
Cint, Ref{Cint}, Ref{Cdouble}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
cilm, size(cilm, 2), d, lat, lon, nmax, lmax, norm, chi2′, csphase,
weights, exitstatus′)
if exitstatus === nothing
exitstatus′[] ≠ 0 && error("SHExpandLSQ!: Error code $(exitstatus′[])")
else
exitstatus[] = exitstatus′[]
end
return cilm, Cdouble(chi2′[])
end
export SHExpandLSQ
function SHExpandLSQ(d::AbstractVector{Cdouble}, lat::AbstractVector{Cdouble},
lon::AbstractVector{Cdouble}, nmax::Integer, lmax::Integer;
norm::Integer=1, csphase::Integer=1,
weights::Vector{Cdouble},
exitstatus::Optional{Ref{<:Integer}}=nothing)
@assert lmax ≥ 0
cilm = Array{Cdouble}(undef, 2, lmax + 1, lmax + 1)
_, chi2 = SHExpandLSQ!(cilm, d, lat, lon, nmax, lmax; norm=norm,
csphase=csphase, weights=weights, exitstatus)
return cilm, chi2
end
export MakeGrid2d!
function MakeGrid2d!(grid::AbstractArray{Cdouble,2},
cilm::AbstractArray{Cdouble,3}, lmax::Integer,
interval::Union{AbstractFloat,Integer}; norm::Integer=1,
csphase::Integer=1,
f::Optional{Union{AbstractFloat,Integer}}=nothing,
a::Optional{Union{AbstractFloat,Integer}}=nothing,
north::Union{AbstractFloat,Integer}=90.0,
south::Union{AbstractFloat,Integer}=-90.0,
east::Union{AbstractFloat,Integer}=360.0,
west::Union{AbstractFloat,Integer}=0.0,
dealloc::Bool=false,
exitstatus::Optional{Ref{<:Integer}}=nothing)
@assert interval > 0
@assert size(grid, 1) ≥ floor(Int, 180 / interval + 1)
@assert size(grid, 2) ≥ floor(Int, 360 / interval + 1)
@assert size(cilm, 1) == 2
@assert lmax ≥ 0
@assert size(cilm, 2) ≥ lmax + 1
@assert size(cilm, 3) == size(cilm, 2)
@assert (f !== nothing) == (a !== nothing)
nlat = Ref{Cint}()
nlong = Ref{Cint}()
exitstatus′ = Ref{Cint}()
ccall((:MakeGrid2d, libSHTOOLS), Cvoid,
(Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Cdouble,
Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cdouble},
Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
Ref{Cint}, Ref{Cint}), grid, size(grid, 1), size(grid, 2), cilm,
size(cilm, 2), lmax, interval, nlat, nlong, norm, csphase,
optional(f, Ptr{Cdouble}()), optional(a, Ptr{Cdouble}()), north,
south, east, west, dealloc, exitstatus′)
if exitstatus === nothing
exitstatus′[] ≠ 0 && error("MakeGrid2d!: Error code $(exitstatus′[])")
else
exitstatus[] = exitstatus′[]
end
return grid, Int(nlat[]), Int(nlong[])
end
export MakeGrid2d
function MakeGrid2d(cilm::AbstractArray{Cdouble,3}, lmax::Integer,
interval::Union{AbstractFloat,Integer}; norm::Integer=1,
csphase::Integer=1,
f::Optional{Union{AbstractFloat,Integer}}=nothing,
a::Optional{Union{AbstractFloat,Integer}}=nothing,
north::Union{AbstractFloat,Integer}=90.0,
south::Union{AbstractFloat,Integer}=-90.0,
east::Union{AbstractFloat,Integer}=360.0,
west::Union{AbstractFloat,Integer}=0.0, dealloc::Bool=false,
exitstatus::Optional{Ref{<:Integer}}=nothing)
@assert interval > 0
grid = Array{Cdouble}(undef, floor(Int, 180 / interval + 1),
floor(Int, 360 / interval + 1))
MakeGrid2d!(grid, cilm, lmax, interval; norm=norm, csphase=csphase, f=f,
a=a, north=north, south=south, east=east, west=west,
dealloc=dealloc, exitstatus=exitstatus)
return grid
end
export MakeGridPoint
function MakeGridPoint(cilm::AbstractArray{Cdouble,3}, lmax::Integer,
lat::Union{AbstractFloat,Integer},
lon::Union{AbstractFloat,Integer}; norm::Integer=1,
csphase::Integer=1, dealloc::Bool=false)
@assert lmax ≥ 0
@assert size(cilm, 1) == 2
@assert size(cilm, 2) ≥ lmax + 1
@assert size(cilm, 3) ≥ size(cilm, 2)
@assert norm ∈ (1, 2, 3, 4)
@assert csphase ∈ (1, -1)
value = ccall((:MakeGridPoint, libSHTOOLS), Cdouble,
(Ptr{Cdouble}, Cint, Cint, Cdouble, Cdouble, Ref{Cint},
Ref{Cint}, Ref{Cint}), cilm, size(cilm, 2), lmax, lat, lon,
norm, csphase, dealloc)
return Float64(value)
export MakeGridPointC
function MakeGridPointC(cilm::AbstractArray{Complex{Cdouble},3}, lmax::Integer,
lat::Union{AbstractFloat,Integer},
lon::Union{AbstractFloat,Integer}; norm::Integer=1,
csphase::Integer=1, dealloc::Bool=false)
@assert lmax ≥ 0
@assert size(cilm, 1) == 2
@assert size(cilm, 2) ≥ lmax + 1
@assert size(cilm, 3) ≥ size(cilm, 2)
@assert norm ∈ (1, 2, 3, 4)
@assert csphase ∈ (1, -1)
value = ccall((:MakeGridPointC, libSHTOOLS), Complex{Cdouble},
(Ptr{Complex{Cdouble}}, Cint, Cint, Cdouble, Cdouble,
Ref{Cint}, Ref{Cint}, Ref{Cint}), cilm, size(cilm, 2), lmax,
lat, lon, norm, csphase, dealloc)
return Complex{Float64}(value)
export SHMultiply!
function SHMultiply!(cilmout::AbstractArray{Cdouble,3},
cilm1::AbstractArray{Cdouble,3}, lmax1::Integer,
cilm2::AbstractArray{Cdouble,3}, lmax2::Integer;
precomp::Bool=false, norm::Integer=1, csphase::Integer=1,
exitstatus::Optional{Ref{<:Integer}}=nothing)
@assert lmax1 ≥ 0
@assert lmax2 ≥ 0
@assert size(cilmout, 1) == 2
@assert size(cilmout, 2) ≥ lmax1 + lmax2 + 1
@assert size(cilmout, 3) == size(cilmout, 2)
@assert size(cilm1, 1) == 2
@assert size(cilm1, 2) ≥ lmax1 + 1
@assert size(cilm1, 3) == size(cilm1, 2)
@assert size(cilm2, 1) == 2
@assert size(cilm2, 2) ≥ lmax2 + 1
@assert size(cilm2, 3) == size(cilm2, 2)
@assert norm ∈ (1, 2, 3, 4)
@assert csphase ∈ (1, -1)
exitstatus′ = Ref{Cint}()
ccall((:SHMultiply, libSHTOOLS), Cvoid,
(Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint,
Cint, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}), cilmout,
size(cilmout, 2), cilm1, size(cilm1, 2), lmax1, cilm2, size(cilm2, 2),
lmax2, precomp, norm, csphase, exitstatus′)
if exitstatus === nothing
exitstatus′[] ≠ 0 && error("SHMultiply!: Error code $(exitstatus′[])")
else
exitstatus[] = exitstatus′[]
end
return cilmout
end
export SHMultiply
function SHMultiply(cilm1::AbstractArray{Cdouble,3}, lmax1::Integer,
cilm2::AbstractArray{Cdouble,3}, lmax2::Integer;
precomp::Bool=false, norm::Integer=1, csphase::Integer=1,
exitstatus::Optional{Ref{<:Integer}}=nothing)
@assert lmax1 ≥ 0
@assert lmax2 ≥ 0
cilmout = Array{Cdouble}(undef, 2, lmax1 + lmax2 + 1, lmax1 + lmax2 + 1)
SHMultiply!(cilmout, cilm1, lmax1, cilm2, lmax2; precomp=precomp, norm=norm,
csphase=csphase, exitstatus=exitstatus)
return cilmout
end
end
end
@testset "SHExpandLSQ" begin
Random.seed!(0)
lmax = 2
n = (lmax + 1)^2
lat = π * rand(n)
lon = 2π * rand(n)
d = randn(n)
weights = Float64[1 for i in 1:n]
cilm, chi2 = SHExpandLSQ(d, lat, lon, n, lmax; weights=weights)
@test size(cilm) == (2, lmax + 1, lmax + 1)
@test abs(chi2) < 10 * eps(Float64)
end
@testset "MakeGrid2d" begin
Random.seed!(0)
lmax = 4
cilm = randn(2, lmax + 1, lmax + 1)
interval = 15 # grid spacing in degrees
grid = MakeGrid2d(cilm, lmax, interval)
@test size(grid) == (180 ÷ interval + 1, 360 ÷ interval + 1)
end
@testset "MakeGridPoint" begin
Random.seed!(0)
# Choose random values
lmax = 2
nmax = (lmax + 1)^2
lat = π * rand(nmax)
lon = 2π * rand(nmax)
values = randn(nmax)
weights = Float64[1 for n in 1:nmax]
# Expand
cilm, chi2 = SHExpandLSQ(values, lat, lon, nmax, lmax; weights=weights)
@test size(cilm) == (2, lmax + 1, lmax + 1)
@test abs(chi2) < 10 * eps(Float64)
# Convert back to values (they will be low-pass filtered)
values′ = Float64[MakeGridPoint(cilm, lmax, lat[n], lon[n]) for n in 1:nmax]
# Expand again (result must be the same)
cilm′, chi2 = SHExpandLSQ(values′, lat, lon, nmax, lmax; weights=weights)
@test abs(chi2) < 10 * eps(Float64)
@test isapprox(cilm′, cilm)
# Convert back to values (result must be the same this time)
values″ = Float64[MakeGridPoint(cilm′, lmax, lat[n], lon[n])
for n in 1:nmax]
@test isapprox(values″, values′)
@testset "MakeGridPointC" begin
Random.seed!(0)
# Choose random values
lmax = 2
cilm = randn(Complex{Float64}, 2, lmax + 1, lmax + 1)
lat = π * rand()
lon = 2π * rand()
value = MakeGridPointC(cilm, lmax, lat, lon)
@test value isa Complex{Float64}
end
@testset "SHMultiply" begin
Random.seed!(0)
# Choose random values
lmax1 = 2
lmax2 = 3
cilm1 = randn(2, lmax1 + 1, lmax1 + 1)
cilm2 = randn(2, lmax2 + 1, lmax2 + 1)
cilmout = SHMultiply(cilm1, lmax1, cilm2, lmax2)
@test size(cilmout) == (2, lmax1 + lmax2 + 1, lmax1 + lmax2 + 1)
end