add schmidt normalization (#15)

* add schmidt normalization

* update docs with norm details

* mention Ambix norm

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "LegendrePolynomials"
 uuid = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
-version = "0.4.2"
+version = "0.4.3"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

docs/src/index.md

@@ -5,12 +5,14 @@ DocTestSetup = quote
 end
 ```
 
-# Introduction
+## Definition
 
-Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials), [Associated Legendre polynomials](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials), and their derivatives using a 3-term recursion relation (Bonnet’s recursion formula).
+### Legendre polynomials
+
+[Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) satisfy the differential equation
 
 ```math
-P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
+\frac{d}{dx}\left(\left(1-x^{2}\right)\frac{d}{dx}\right)P_{\ell}\left(x\right)=-\ell\left(\ell+1\right)P_{\ell}\left(x\right)
 ```
 
 By default this package evaluates the standard polynomials that satisfy ``P_\ell(1) = 1`` and ``P_0(x) = 1``. These are normalized as
@@ -19,43 +21,29 @@ By default this package evaluates the standard polynomials that satisfy ``P_\ell
 \int_{-1}^1 P_m(x) P_n(x) dx = \frac{2}{2n+1} \delta_{mn}.
 ```
 
-Optionally, normalized polynomials may be evaluated that have an L2 norm of `1`.
-
-Analogous to Legendre polynomials, one may evaluate associated Legendre polynomials using a 3-term recursion relation. This is evaluated by iterating over the normalized associated Legendre functions, and multiplying the norm at the final stage. Such an iteration avoids floating-point overflow.
+### Associated Legendre polynomials
 
-The relation used in evaluating normalized associated Legendre polynomials ``\bar{P}_{\ell}^{m}\left(x\right)`` is
-
-```math
-\bar{P}_{\ell}^{m}\left(x\right)=\alpha_{\ell m}\left(x\bar{P}_{\ell-1}^{m}\left(x\right)-\frac{1}{\alpha_{\ell-1m}}\bar{P}_{\ell-2}^{m}\left(x\right)\right),
-```
+[Associated Legendre polynomials](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials) satisfy the differential equation
 
-where
 ```math
-\alpha_{\ell m}=\sqrt{\frac{\left(2\ell+1\right)\left(2\ell-1\right)}{\left(\ell-m\right)\left(\ell+m\right)}},
+\left[\frac{d}{dx}\left(\left(1-x^{2}\right)\frac{d}{dx}\right)-\frac{m^{2}}{1-x^{2}}\right]P_{\ell}^{m}\left(x\right)=-\ell\left(\ell+1\right)P_{\ell}^{m}\left(x\right),
 ```
 
-The starting value for a specific ``m`` is obtained by setting ``\ell = m``, in which case we use the explicit relation
+and are related to the Legendre polynomials through
 
 ```math
-\bar{P}_{m}^{m}\left(x\right)=\left(-1\right)^{m}\left(2m-1\right)!!\sqrt{\frac{\left(2m+1\right)}{2\left(2m\right)!}}\left(1-x^{2}\right)^{\left(m/2\right)}.
+P_{\ell}^{m}\left(x\right)=\left(-1\right)^{m}\left(1-x^{2}\right)^{m/2}\frac{d^{m}}{dx^{m}}\left(P_\ell\left(x\right)\right).
 ```
 
-The polynomials, thus defined, include the Condon-Shortley phase ``(-1)^m``, and may be evaluated using the standard normalization as
+For ``m=0``, we obtain
 
 ```math
-P_{\ell}^{m}\left(x\right)=\sqrt{\frac{2\left(\ell+m\right)!}{\left(2\ell+1\right)\left(\ell-m\right)!}}\bar{P}_{\ell}^{m}\left(x\right).
+P_{\ell}^{0}\left(x\right)=P_\ell\left(x\right).
 ```
 
-There are six main functions:
-
-* [`Pl(x,l; [norm = Val(:standard)])`](@ref Pl): this evaluates the Legendre polynomial for a given degree `l` at the argument `x`. The argument needs to satisfy `-1 <= x <= 1`.
-* [`collectPl(x; lmax, [lmin = 0], [norm = Val(:standard)])`](@ref collectPl): this evaluates all the polynomials for `l` lying in `0:lmax` at the argument `x`. As before the argument needs to lie in the domain of validity. Functionally this is equivalent to `Pl.(x, lmin:lmax)`, except `collectPl` evaluates the result in one pass, and is therefore faster. There is also the in-place version [`collectPl!`](@ref) that uses a pre-allocated array.
-* [`Plm(x, l, m; [norm = Val(:standard)])`](@ref Plm): this evaluates the associated Legendre polynomial ``P_\ell^m(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
-* [`collectPlm(x; m, lmax, [lmin = abs(m)], [norm = Val(:standard)])`](@ref collectPlm): this evaluates the associated Legendre polynomials with coefficient `m` for `l = lmin:lmax`. There is also an in-place version [`collectPlm!`](@ref) that uses a pre-allocated array.
-* [`dnPl(x, l, n)`](@ref dnPl): this evaluates the ``n``-th derivative of the Legendre polynomial ``P_\ell(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
-* [`collectdnPl(x; n, lmax)`](@ref collectdnPl): this evaluates the ``n``-th derivative of all the Legendre polynomials for `l = 0:lmax`. There is also an in-place version [`collectdnPl!`](@ref) that uses a pre-allocated array.
+The Condon-Shortley phase is included by default, and may be toggled by using the keyword argument `csphase`.
 
-# Quick Start
+## Quick Start
 
 Evaluate the Legendre polynomial for one `l` at an argument`x` as `Pl(x, l)`:
 
@@ -120,7 +108,75 @@ julia> collectdnPl(0.5, lmax = 5, n = 3)
  65.625
 ```
 
-# Reference
+## Recursion relations
+
+Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials), [Associated Legendre polynomials](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials), and their derivatives using a 3-term recursion relation (Bonnet’s recursion formula).
+
+```math
+P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
+```
+
+Analogous to Legendre polynomials, one may evaluate associated Legendre polynomials using a 3-term recursion relation. This is evaluated by iterating over the normalized associated Legendre functions, and multiplying the norm at the final stage. Such an iteration avoids floating-point overflow.
+
+The relation used in evaluating normalized associated Legendre polynomials ``\bar{P}_{\ell}^{m}\left(x\right)`` is
+
+```math
+\bar{P}_{\ell}^{m}\left(x\right)=\alpha_{\ell m}\left(x\bar{P}_{\ell-1}^{m}\left(x\right)-\frac{1}{\alpha_{\ell-1m}}\bar{P}_{\ell-2}^{m}\left(x\right)\right),
+```
+
+where
+```math
+\alpha_{\ell m}=\sqrt{\frac{\left(2\ell+1\right)\left(2\ell-1\right)}{\left(\ell-m\right)\left(\ell+m\right)}},
+```
+
+The starting value for a specific ``m`` is obtained by setting ``\ell = m``, in which case we use the explicit relation
+
+```math
+\bar{P}_{m}^{m}\left(x\right)=\left(-1\right)^{m}\left(2m-1\right)!!\sqrt{\frac{\left(2m+1\right)}{2\left(2m\right)!}}\left(1-x^{2}\right)^{\left(m/2\right)}.
+```
+
+The polynomials, thus defined, include the Condon-Shortley phase ``(-1)^m``, and may be evaluated using the standard normalization as
+
+```math
+P_{\ell}^{m}\left(x\right)=\sqrt{\frac{2\left(\ell+m\right)!}{\left(2\ell+1\right)\left(\ell-m\right)!}}\bar{P}_{\ell}^{m}\left(x\right).
+```
+
+## Main functions
+
+There are six main functions:
+
+* [`Pl(x,l; [norm = Val(:standard)])`](@ref Pl): this evaluates the Legendre polynomial for a given degree `l` at the argument `x`. The argument needs to satisfy `-1 <= x <= 1`.
+* [`collectPl(x; lmax, [lmin = 0], [norm = Val(:standard)])`](@ref collectPl): this evaluates all the polynomials for `l` lying in `0:lmax` at the argument `x`. As before the argument needs to lie in the domain of validity. Functionally this is equivalent to `Pl.(x, lmin:lmax)`, except `collectPl` evaluates the result in one pass, and is therefore faster. There is also the in-place version [`collectPl!`](@ref) that uses a pre-allocated array.
+* [`Plm(x, l, m; [norm = Val(:standard)], [csphase = true])`](@ref Plm): this evaluates the associated Legendre polynomial ``P_\ell^m(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
+* [`collectPlm(x; m, lmax, [lmin = abs(m)], [norm = Val(:standard)], [csphase = true])`](@ref collectPlm): this evaluates the associated Legendre polynomials with coefficient `m` for `l = lmin:lmax`. There is also an in-place version [`collectPlm!`](@ref) that uses a pre-allocated array.
+* [`dnPl(x, l, n)`](@ref dnPl): this evaluates the ``n``-th derivative of the Legendre polynomial ``P_\ell(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
+* [`collectdnPl(x; n, lmax)`](@ref collectdnPl): this evaluates the ``n``-th derivative of all the Legendre polynomials for `l = 0:lmax`. There is also an in-place version [`collectdnPl!`](@ref) that uses a pre-allocated array.
+
+## Normalization options
+
+By default, the Legendre and associated Legendre polynomials are not normalized.
+One may specify the normalization through the keyword argument `norm`.
+The normalization options accepted are
+
+* `norm = Val(:standard)`: standard, unnormalized polynomials. This is the default option. Choosing this option will lead to the standard Legendre or associated Legendre polynomials that satisfy ``P_\ell(0)=1`` and ``P_{\ell0}(0)=1``
+* `norm = Val(:normalized)`: fully normalized polynomials, with an L2 norm of 1. These are defined as
+```math
+  \bar{P}_{\ell m}\left(x\right)=\sqrt{\frac{2\ell+1}{2}\frac{\left(\ell-m\right)!}{\left(\ell+m\right)!}}P_{\ell m}\left(x\right)
+```
+* `norm = Val(:schmidtquasi)`: Schmidt quasi-normalized polynomials, also known as Schmidt semi-normalized polynomials. These have an L2 norm of ``\sqrt{2(2-\delta_{m0})/(2\ell+1)}``. These are defined as
+```math
+  \bar{P}_{\ell m}\left(x\right)=\sqrt{\left(2-\delta_{m0}\right)\frac{\left(\ell-m\right)!}{\left(\ell+m\right)!}}P_{\ell m}\left(x\right)
+```
+* `norm = Val(:schmidt)`: Schmidt normalized polynomials. These have an L2 norm of ``\sqrt{2(2-\delta_{m0})}``. These are defined as
+```math
+  \bar{P}_{\ell m}\left(x\right)=\sqrt{\left(2-\delta_{m0}\right)\left(2\ell+1\right)\frac{\left(\ell-m\right)!}{\left(\ell+m\right)!}}P_{\ell m}\left(x\right)
+```
+
+!!! note
+    Irrespective of the norm specified, the 3-term recursion relations used are stable ones,
+    so polynomials of high degrees may be computed without overflow.
+
+## Reference
 
 ```@autodocs
 Modules = [LegendrePolynomials]

docs/src/tutorial.md

@@ -4,16 +4,7 @@ DocTestSetup = :(using LegendrePolynomials)
 
 ## Computing normalized Legendre polynomials
 
-By default, the Legendre and associated Legendre polynomials are not normalized.
-One may specify the normalization through the keyword argument `norm`.
-The normalization options accepted are
-
-* `norm = Val(:standard)`: standard, unnormalized polynomials. This is the default option.
-* `norm = Val(:normalized)`: fully normalized polynomials with an L2 norm of 1
-
-!!! note
-	Irrespective of the norm specified, the 3-term recursion relations used are stable ones,
-	so polynomials of high degrees may be computed without overflow.
+The norm of the polynomials may be specified through the keyword argument `norm`. The various normalization options are listed in the [Normalization options](@ref) section.
 
 ```jldoctest
 julia> l = 3;
@@ -33,6 +24,16 @@ julia> Plm(0.5, l, m, norm = Val(:normalized))
 2.172276347346834e-187
 ```
 
+Starting from these, other normalization such as the [Ambix SN3D format](https://en.wikipedia.org/wiki/Ambisonic_data_exchange_formats#SN3D) may be constructed as
+
+```julia
+julia> Plmambix(x, l, m) = Plm(x, l, m, norm=Val(:schmidtquasi)) / √(4pi)
+Plmambix (generic function with 1 method)
+
+julia> Plmambix(0.5, 2, 1)
+-0.21157109383040862
+```
+
 ## Condon-Shortley phase
 
 By default, the associated Legendre polynomials are computed by including the Condon-Shortley phase ``(-1)^m``.

src/LegendrePolynomials.jl

@@ -200,25 +200,31 @@ function Base.iterate(iter::LegendrePolynomialIterator{T,<:Integer}, state) wher
     return Pl, (l+1, nextstate)
 end
 
+throw_normerr(norm) =
+    throw(ArgumentError("norm = $norm undefined,"*
+        " valid norms are Val(:standard), Val(:normalized), Val(:schmidtquasi) or Val(:schmidt)"))
+
 maybenormalize(P, l, ::Val{:normalized}) = oftype(P, P * √(l+1/2))
-maybenormalize(P, l, ::Val{:standard}) = P
-maybenormalize(P, l, norm) = throw(ArgumentError("norm = $norm undefined, valid norms are Val(:standard) and Val(:normalized)"))
+maybenormalize(P, l, ::Union{Val{:standard}, Val{:schmidtquasi}}) = P
+maybenormalize(P, l, ::Val{:schmidt}) = oftype(P, P * √(2l+1))
+maybenormalize(P, l, norm) = throw_normerr(norm)
 
-function maybenormalize(P, l, m, norm::Val{:normalized}, csphase = true)
-    if m == 0
-        return maybenormalize(P, l, norm)
-    else
-        return P
-    end
+plm_norm(::Val{:normalized}, lm...) = 1.0
+plm_norm(::Val{:standard}, lm...) = plm_norm(lm...)
+plm_norm(::Val{:schmidtquasi}, l, m) = √(2*(2 - iszero(m))/(2l + 1))
+function plm_norm(::Val{:schmidt}, l, m)
+    two = first(promote(2, l, m))
+    √(two*(two - iszero(m)))
 end
-function maybenormalize(P, l, m, norm::Val{:standard}, csphase = true)
+plm_norm(norm, args...) = throw_normerr(norm)
+
+function maybenormalize(P, l, m, norm, csphase = true)
     if m == 0
         return maybenormalize(P, l, norm)
     else
-        return (m < 0 ? neg1pow(m) : 1) * plm_norm(l, m) * P
+        return (m < 0 ? neg1pow(m) : 1) * plm_norm(norm, l, m) * P
     end
 end
-maybenormalize(P, l, m, norm, args...) = throw(ArgumentError("norm = $norm undefined, valid norms are Val(:standard) and Val(:normalized)"))
 
 """
     Pl(x, l::Integer; [norm = Val(:standard)])
@@ -268,8 +274,10 @@ For `m == 0` this function returns Legendre polynomials.
 
 # Optional keyword arguments
 * `norm`: The normalization used in the function. Possible
-    options are `Val(:standard)` (default) and `Val(:normalized)`.
-    The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
+    options are `Val(:standard)` (default), `Val(:normalized)`, `Val(:schmidtquasi)` or `Val(:schmidt)`.
+    * The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
+    * The functions obtained with `norm = Val(:schmidtquasi)` have an L2 norm of ``\\sqrt{\\frac{2(2-\\delta_{m0})}{2l+1}}``.
+    * The functions obtained with `norm = Val(:schmidt)` have an L2 norm of ``\\sqrt{2(2-\\delta_{m0})}``.
 * `csphase::Bool`: The Condon-shortley phase ``(-1)^m``, which is included by default.
 
 # Examples
@@ -474,8 +482,10 @@ Returns `v` with indices `lmin:lmax`, with `v[l] == Plm(x, l, m; kwargs...)`.
 
 # Optional keyword arguments
 * `norm`: The normalization used in the function. Possible
-    options are `Val(:standard)` (default) and `Val(:normalized)`.
-    The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
+    options are `Val(:standard)` (default), `Val(:normalized)`, `Val(:schmidtquasi)` or `Val(:schmidt)`.
+    * The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
+    * The functions obtained with `norm = Val(:schmidtquasi)` have an L2 norm of ``\\sqrt{\\frac{2(2-\\delta_{m0})}{2l+1}}``.
+    * The functions obtained with `norm = Val(:schmidt)` have an L2 norm of ``\\sqrt{2(2-\\delta_{m0})}``.
 * `csphase::Bool`: The Condon-shortley phase ``(-1)^m``, which is included by default.
 * `Tnorm`: Type of the normalization factor, which is dynamicall inferred by default,
     and is used in allocating an appropriate container.
@@ -525,8 +535,10 @@ At output, `v[l + firstindex(v)] == Plm(x, l, m; kwargs...)` for `l = lmin:lmax`
 
 # Optional keyword arguments
 * `norm`: The normalization used in the function. Possible
-    options are `Val(:standard)` (default) and `Val(:normalized)`.
-    The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
+    options are `Val(:standard)` (default), `Val(:normalized)`, `Val(:schmidtquasi)` or `Val(:schmidt)`.
+    * The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
+    * The functions obtained with `norm = Val(:schmidtquasi)` have an L2 norm of ``\\sqrt{\\frac{2(2-\\delta_{m0})}{2l+1}}``.
+    * The functions obtained with `norm = Val(:schmidt)` have an L2 norm of ``\\sqrt{2(2-\\delta_{m0})}``.
 * `csphase::Bool`: The Condon-shortley phase ``(-1)^m``, which is included by default.
 
 # Examples