Add tests for parity and Legendre polynomials (#8)

* add tests

* Add propagate_inbounds

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "SphericalHarmonics"
 uuid = "c489a379-e885-57ff-9236-bd896d33c250"
-version = "0.1.9"
+version = "0.1.10"
 
 [deps]
 ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
@@ -10,14 +10,16 @@ SphericalHarmonicModes = "0e9554e2-b38b-11e9-16d7-9d9abfec665a"
 
 [compat]
 ElasticArrays = "1.2"
+LegendrePolynomials = "0.3"
 SphericalHarmonicArrays = "0.4.4"
 SphericalHarmonicModes = "0.4.7"
 julia = "1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
+LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "HCubature", "Test"]
+test = ["Aqua", "HCubature", "LegendrePolynomials", "Test"]

src/SphericalHarmonics.jl

@@ -40,7 +40,7 @@ Base.parent(P::AssociatedLegendrePolynomials) = P.P
 Base.size(P::AssociatedLegendrePolynomials) = size(parent(P))
 Base.axes(P::AssociatedLegendrePolynomials) = axes(parent(P))
 Base.IndexStyle(::Type{<:AssociatedLegendrePolynomials{<:Any, PLM}}) where {PLM} = IndexStyle(PLM)
-Base.getindex(P::AssociatedLegendrePolynomials, I...) = getindex(parent(P), I...)
+Base.@propagate_inbounds Base.getindex(P::AssociatedLegendrePolynomials, I...) = getindex(parent(P), I...)
 function Base.summary(io::IO, P::AssociatedLegendrePolynomials)
     print(io, "$(length(P))-element AssociatedLegendrePolynomials{$(eltype(P))} for lmax = $(Int(P.lmax))")
     if P.initialized

test/runtests.jl

@@ -3,6 +3,7 @@ using SphericalHarmonicModes
 using Test
 using HCubature
 using Aqua
+using LegendrePolynomials
 
 import SphericalHarmonics: NorthPole, SouthPole, allocate_y, allocate_p, RealHarmonics, ComplexHarmonics
 
@@ -888,3 +889,58 @@ end
         s_exp = "10-element AssociatedLegendrePolynomials{BigFloat} for lmax = 3 (uninitialized)"
     end
 end
+
+@testset "Parity" begin
+    lmax = 10
+    S1 = SphericalHarmonics.cache(lmax)
+    S2 = SphericalHarmonics.cache(lmax)
+    for θ in LinRange(0, pi, 10)
+        P1 = computePlmcostheta!(S1, θ)
+        P2 = computePlmcostheta!(S2, pi - θ)
+        for l in 0:lmax, m in 0:l
+            @test isapprox(P2[(l,m)], (-1)^(l + m) * P1[(l,m)], atol = 1e-13, rtol = 1e-13)
+        end
+        for ϕ in LinRange(0, 2pi, 10)
+            Y1 = computeYlm!(S1, θ, ϕ)
+            Y2 = computeYlm!(S2, pi - θ, pi + ϕ)
+            for l in 0:lmax, m in -l:l
+                @test isapprox(Y2[(l,m)], (-1)^l * Y1[(l,m)], atol = 1e-13, rtol = 1e-13)
+            end
+        end
+    end
+end
+
+@testset "m = 0 Legendre" begin
+    # for m = 0, the spherical harmonics are normalized Legendre Polynomials and are purely real
+    lmax = 100
+    @testset "ComplexHarmonics" begin
+        S = SphericalHarmonics.cache(lmax)
+        LP = zeros(0:lmax)
+        for θ in LinRange(0, pi, 10)
+            LegendrePolynomials.collectPl!(LP, cos(θ), lmax = lmax)
+            computePlmcostheta!(S, θ)
+            for ϕ in LinRange(0, 2pi, 10)
+                computeYlm!(S, θ, ϕ)
+                for l in 0:lmax
+                    @test isapprox(imag(S.Y[(l,0)]), 0, atol = 1e-13)
+                    @test isapprox(S.Y[(l,0)], √((2l+1)/4pi) * LP[l], atol = 1e-13, rtol = 1e-13)
+                end
+            end
+        end
+    end
+    @testset "RealHarmonics" begin
+        S = SphericalHarmonics.cache(lmax, SHType = RealHarmonics())
+        LP = zeros(0:lmax)
+        for θ in LinRange(0, pi, 10)
+            LegendrePolynomials.collectPl!(LP, cos(θ), lmax = lmax)
+            computePlmcostheta!(S, θ)
+            for ϕ in LinRange(0, 2pi, 10)
+                computeYlm!(S, θ, ϕ)
+                for l in 0:lmax
+                    @test isapprox(imag(S.Y[(l,0)]), 0, atol = 1e-13)
+                    @test isapprox(S.Y[(l,0)], √((2l+1)/4pi) * LP[l], atol = 1e-13, rtol = 1e-13)
+                end
+            end
+        end
+    end
+end