From 74d2d37de5ab92074504a0a10de6e8a212fa9094 Mon Sep 17 00:00:00 2001 From: Justin Willmert Date: Mon, 13 Dec 2021 23:03:38 -0600 Subject: [PATCH] =?UTF-8?q?Add=20new=204=CF=80=20sphere=20(geodesy)=20norm?= =?UTF-8?q?alization?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Defined such that the integrating square of this normalization times the complex exponential in azimuth equals the surface area of the unit sphere. This complements the spherical normalization `LegendreSphereNorm` which instead integrates over the sphere to unity. Also adds an intermediate `AbstractLegendreSphereNorm <: AbstractLegendreNorm` in the abstract type hierarchy to consolidate the coefficients for all three of `LegendreFourPiNorm`, `LegendreOrthoNorm`, and `LegendreSphereNorm`. The original motivation was to see if replacing the initial condition of 1 / sqrt(4π) which must necessarily be rounded (because π is trancendental) with an exactly-representable value (1.0 for the 4π normalization) could lead to better numerical precision. Namely, I thought that something like summing over many terms could show the accumulated bias, but I haven't actually found that to be true in practice. Nonetheless, having another option for a built-in normalizations is still useful. --- Project.toml | 2 +- docs/src/lib/public.md | 2 ++ docs/src/man/devdocs.md | 2 +- src/Legendre.jl | 7 +++++-- src/aliases.jl | 8 ++++++++ src/calculation.jl | 2 +- src/norm_sphere.jl | 36 ++++++++++++++++++++++++++++-------- test/analytic.jl | 18 +++++++++++++++++- test/norms.jl | 3 +++ 9 files changed, 66 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index d05b98c..a302062 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Legendre" uuid = "7642852e-7f09-11e9-134e-0940411082b6" authors = ["Justin Willmert "] -version = "0.2.2" +version = "0.2.3" [compat] julia = "1.2" diff --git a/docs/src/lib/public.md b/docs/src/lib/public.md index d563ba4..f0ebd4d 100644 --- a/docs/src/lib/public.md +++ b/docs/src/lib/public.md @@ -16,6 +16,7 @@ Nlm ```@docs AbstractLegendreNorm LegendreUnitNorm +LegendreFourPiNorm LegendreOrthoNorm LegendreSphereNorm LegendreNormCoeff @@ -36,6 +37,7 @@ Plm! There are also aliases for pre-computed coefficients of the provided normalizations. ```@docs LegendreUnitCoeff +LegendreFourPiCoeff LegendreOrthoCoeff LegendreSphereCoeff ``` diff --git a/docs/src/man/devdocs.md b/docs/src/man/devdocs.md index 6ef812f..0815262 100644 --- a/docs/src/man/devdocs.md +++ b/docs/src/man/devdocs.md @@ -153,7 +153,7 @@ normalization trait type as the first argument. up a type-stable algorithm, which we'll ignore here for the sake of simplicity.) ```jldoctest λNorm julia> initcond(::λNorm, T::Type) = sqrt(1 / 4π) -initcond (generic function with 5 methods) +initcond (generic function with 6 methods) ``` Finally, we provide methods which encode the cofficients as well: ```jldoctest λNorm diff --git a/src/Legendre.jl b/src/Legendre.jl index ce3b181..4ea0116 100644 --- a/src/Legendre.jl +++ b/src/Legendre.jl @@ -9,8 +9,11 @@ export AbstractLegendreNorm include("interface.jl") # Specific normalizations -export LegendreUnitNorm, LegendreOrthoNorm, LegendreSphereNorm, LegendreNormCoeff, - LegendreUnitCoeff, LegendreOrthoCoeff, LegendreSphereCoeff +export LegendreNormCoeff, + LegendreUnitNorm, LegendreUnitCoeff, + LegendreFourPiNorm, LegendreFourPiCoeff, + LegendreOrthoNorm, LegendreOrthoCoeff, + LegendreSphereNorm, LegendreSphereCoeff include("norm_unit.jl") include("norm_sphere.jl") include("norm_table.jl") diff --git a/src/aliases.jl b/src/aliases.jl index 0880be4..de01f92 100644 --- a/src/aliases.jl +++ b/src/aliases.jl @@ -6,6 +6,14 @@ normalization. Alias for `LegendreNormCoeff{LegendreUnitNorm,T}`. """ LegendreUnitCoeff{T} = LegendreNormCoeff{LegendreUnitNorm,T} +""" + LegendreFourPiCoeff{T} + +Table type of precomputed recursion relation coefficients for the ``4\\pi`` spherical +harmonic normalization. Alias for `LegendreNormCoeff{LegendreFourPiNorm,T}`. +""" +LegendreFourPiCoeff{T} = LegendreNormCoeff{LegendreFourPiNorm,T} + """ LegendreOrthoCoeff{T} diff --git a/src/calculation.jl b/src/calculation.jl index a7ab5f2..afc960d 100644 --- a/src/calculation.jl +++ b/src/calculation.jl @@ -279,7 +279,7 @@ end @static if VERSION < v"1.3.0-alpha" # Prior to julia-1.3.0, abstract types could not have methods attached to them. # Provide for just the defined normalization types. - for N in (LegendreUnitNorm,LegendreSphereNorm,LegendreOrthoNorm,LegendreNormCoeff) + for N in (LegendreUnitNorm,LegendreSphereNorm,LegendreOrthoNorm,LegendreFourPiNorm,LegendreNormCoeff) @eval begin @inline function (norm::$N)(l, m, x) return legendre(norm, l, m, x) diff --git a/src/norm_sphere.jl b/src/norm_sphere.jl index fa15468..59acb45 100644 --- a/src/norm_sphere.jl +++ b/src/norm_sphere.jl @@ -1,5 +1,24 @@ # Implements the legendre interface for the unit normalization case +abstract type AbstractLegendreSphereNorm <: AbstractLegendreNorm end + +""" + struct LegendreFourPiNorm <: AbstractLegendreNorm end + +Trait type denoting the ``4\\pi`` (geodesy) normalization of the associated Legendre +functions ``F_\\ell^m(x)``. +The orthonormal normalization is defined such that +```math + \\int_{-1}^{1} \\left[ F_\\ell^m(x) \\right]^2 \\,dx = 2 +``` +The normalization factor with respect to the standard (unnormalized) Legendre +functions ``P_\\ell^m(x)`` ([`LegendreUnitNorm`](@ref)) is given by +```math + F_\\ell^m(x) \\equiv \\sqrt{2\\pi(2\\ell+1) \\frac{(\\ell-m)!}{(\\ell+m)!}} P_\\ell^m(x) +``` +""" +struct LegendreFourPiNorm <: AbstractLegendreSphereNorm end + """ struct LegendreOrthoNorm <: AbstractLegendreNorm end @@ -15,7 +34,7 @@ functions ``P_\\ell^m(x)`` ([`LegendreUnitNorm`](@ref)) is given by O_\\ell^m(x) \\equiv \\sqrt{\\frac{2\\ell+1}{2} \\frac{(\\ell-m)!}{(\\ell+m)!}} P_\\ell^m(x) ``` """ -struct LegendreOrthoNorm <: AbstractLegendreNorm end +struct LegendreOrthoNorm <: AbstractLegendreSphereNorm end """ struct LegendreSphereNorm <: AbstractLegendreNorm end @@ -32,8 +51,11 @@ functions ``P_\\ell^m(x)`` ([`LegendreUnitNorm`](@ref)) is given by \\lambda_\\ell^m(x) \\equiv \\sqrt{\\frac{2\\ell+1}{4\\pi} \\frac{(\\ell-m)!}{(\\ell+m)!}} P_\\ell^m(x) ``` """ -struct LegendreSphereNorm <: AbstractLegendreNorm end +struct LegendreSphereNorm <: AbstractLegendreSphereNorm end +@inline function initcond(::LegendreFourPiNorm, ::Type{T}) where T + return one(T) +end @inline function initcond(::LegendreOrthoNorm, ::Type{T}) where T return sqrt(inv(T(2))) end @@ -44,8 +66,6 @@ end return inv(sqrt(4 * convert(T, π))) end -const OrthoOrSphereNorm = Union{LegendreOrthoNorm,LegendreSphereNorm} - # Version of sqrt() which skips the domain (x < 0) check for the IEEE floating point types. # For nonstandard number types, just fall back to a regular sqrt() since eliminating the # domain check is probably no longer the dominant contributor to not vectorizing. @@ -54,7 +74,7 @@ unchecked_sqrt(x::T) where {T <: Integer} = unchecked_sqrt(float(x)) unchecked_sqrt(x) = Base.sqrt(x) @inline function -coeff_μ(::OrthoOrSphereNorm, ::Type{T}, l::Integer) where T +coeff_μ(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer) where T # The direct derivation of the normalization constant gives # return sqrt(one(T) + inv(convert(T, 2l))) # but when comparing results for T ∈ (Float64,BigFloat), the Float64 results differ by @@ -83,12 +103,12 @@ coeff_μ(::OrthoOrSphereNorm, ::Type{T}, l::Integer) where T end @inline function -coeff_ν(::OrthoOrSphereNorm, ::Type{T}, l::Integer) where T +coeff_ν(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer) where T return unchecked_sqrt(convert(T, 2l + 1)) end @inline function -coeff_α(::OrthoOrSphereNorm, ::Type{T}, l::Integer, m::Integer) where T +coeff_α(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer, m::Integer) where T lT = convert(T, l) mT = convert(T, m) # Write factors in two pieces to make compiler's job easier. In the case where @@ -100,7 +120,7 @@ coeff_α(::OrthoOrSphereNorm, ::Type{T}, l::Integer, m::Integer) where T end @inline function -coeff_β(::OrthoOrSphereNorm, ::Type{T}, l::Integer, m::Integer) where T +coeff_β(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer, m::Integer) where T lT = convert(T, l) mT = convert(T, m) # Write factors in two pieces to make compiler's job easier. In the case where diff --git a/test/analytic.jl b/test/analytic.jl index 0dbe523..e7b45a6 100644 --- a/test/analytic.jl +++ b/test/analytic.jl @@ -35,8 +35,10 @@ end # P_0^0 is constant. Verify the output is invariant for several inputs. @testset "Constant P_0^0 ($T)" for T in NumTypes + @test @inferred(legendre(LegendreFourPiNorm(), 0, 0, T(0.1))) isa T @test @inferred(legendre(LegendreOrthoNorm(), 0, 0, T(0.1))) isa T @test @inferred(legendre(LegendreSphereNorm(), 0, 0, T(0.1))) isa T + @test all(legendre(LegendreFourPiNorm(), 0, 0, x) == one(T) for x in xrange(T)) @test all(legendre(LegendreOrthoNorm(), 0, 0, x) == sqrt(T(0.5)) for x in xrange(T)) @test all(legendre(LegendreSphereNorm(), 0, 0, x) == sqrt(inv(4T(π))) for x in xrange(T)) end @@ -102,7 +104,7 @@ end end end -# The normal P_ℓ^m should be equal to the spherical harmonic normalized +# The orthonormal O_ℓ^m should be equal to the spherical harmonic normalized # λ_ℓ^m if we manually normalize them. @testset "Equality of (2π)^(-1/2) O_ℓ^m and λ_ℓ^m ($T)" for T in NumTypes LMAX = 5 @@ -116,6 +118,20 @@ end end end +# The spherical and 4π-normalized functions should be the same up to the +# differing factor of 1/√4π +@testset "Equality of (4π)^(-1/2) F_ℓ^m and λ_ℓ^m ($T)" for T in NumTypes + LMAX = 5 + atol = max(eps(4T(π)), eps(4.0π)) + plm_4pi = zeros(T, LMAX+1, LMAX+1) + plm_sphr = zeros(T, LMAX+1, LMAX+1) + for x in xrange(T) + legendre!(LegendreFourPiNorm(), plm_4pi, LMAX, LMAX, x) + legendre!(LegendreSphereNorm(), plm_sphr, LMAX, LMAX, x) + @test all(isapprox.(plm_4pi./sqrt(4T(π)), plm_sphr, atol=atol)) + end +end + ################################ # Complex domain ################################ diff --git a/test/norms.jl b/test/norms.jl index b0f504b..b4f8d78 100644 --- a/test/norms.jl +++ b/test/norms.jl @@ -3,6 +3,9 @@ import ..TestSuite @testset "Unit normalization" begin TestSuite.runtests(LegendreUnitNorm()) end +@testset "4π normalization" begin + TestSuite.runtests(LegendreFourPiNorm()) +end @testset "Orthonormal normalization" begin TestSuite.runtests(LegendreOrthoNorm()) end