Add new 4π sphere (geodesy) normalization

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

@@ -1,7 +1,7 @@
 name = "Legendre"
 uuid = "7642852e-7f09-11e9-134e-0940411082b6"
 authors = ["Justin Willmert <[email protected]>"]
-version = "0.2.2"
+version = "0.2.3"
 
 [compat]
 julia = "1.2"

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
 ```

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

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")

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}
 

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)

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

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
 ################################

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