Add an interface for working with partial derivatives and gradients

These polynomials are analytically simple so there's no need for
numerical differentiation (something which is better suited for the
WavefrontErrors); calculating these amounts to manipulating the radial
coefficient vector and tracking the harmonic's sign.

src/Derivatives.jl

@@ -0,0 +1,59 @@
+struct PartialDerivative{T <: Union{RadialPolynomial, Harmonic}} <: Polynomials
+    order::Int
+    inds::NamedTuple{(:j, :n, :m), NTuple{3, Int}}
+    N::Float64
+    R::RadialPolynomial
+    M::Harmonic
+end
+
+struct Gradient{T <: Polynomial}
+    r::PartialDerivative{RadialPolynomial}
+    t::PartialDerivative{Harmonic}
+    Gradient{Polynomial}(Z::Polynomial) = new(derivatives(Z)...)
+end
+
+"""
+    Zernike.Gradient(Z::Polynomial)
+
+Returns ∇Z(ρ, θ).
+"""
+Gradient(Z::Polynomial) = Gradient{Polynomial}(Z)
+
+(g::Gradient)(ρ::Real, θ::Real = 0) = [g.r(ρ, θ), g.t(ρ, θ)]
+
+p(i, order) = ∏(i-order+1:float(i))
+
+s(m, order)::Int = order % 4 ∈ (0, 1+2(m > 0)) || -1
+
+# TODO: this needs unit tests
+"""
+    Zernike.derivatives(Z::Polynomial, order::Int = 1)
+
+Computes the nth order partial derivatives of Z(ρ, θ).
+"""
+function derivatives(Z::Polynomial, order::Int = 1)
+    order < 1 && throw(DomainError(order))
+    (; inds, N, R, M) = Z
+    (; λ, γ, ν) = R
+    (; m) = M
+    λ′ = zeros(length(λ))
+    i = eachindex(λ) .- 1
+    j = @. i ≥ order
+    @. λ′[j] = λ[j] * p(i[j], order)
+    λ′ = shift(λ′, -order)
+    ν′ = Int[νᵢ - order for νᵢ ∈ ν if νᵢ ≥ order]
+    γ′ = Float64[λ′[νᵢ+1] for νᵢ ∈ ν′]
+    N′ = N * s(m, order) * abs(m) ^ order
+    m *= (-1) ^ order
+    R′ = RadialPolynomial(λ′, γ′, ν′)
+    M′ = Harmonic(m)
+    ∂Z_∂ρ = PartialDerivative{RadialPolynomial}(order, inds, N, R′, M)
+    ∂Z_∂θ = PartialDerivative{Harmonic}(order, inds, N′, R, M′)
+    return ∂Z_∂ρ, ∂Z_∂θ
+end
+
+function show(io::IO, ∂::T) where T <: PartialDerivative
+    print(io, T, " order: ", ∂.order)
+end
+
+show(io::IO, ::T) where T <: Gradient = print(io, T)

src/Zernike.jl

@@ -13,7 +13,9 @@ export zernike, wavefront, transform, Z, W, P, WavefrontError, get_j, get_mn,
 
 const public_names = "public \
     radial_coefficients, wavefront_coefficients, transform_coefficients, \
-    metrics, scale, J, Superposition, Product, sieve, format_strings, valid_fringes"
+    metrics, scale, J, Superposition, Product, \
+    sieve, format_strings, valid_fringes, \
+    derivatives, PartialDerivative, Gradient"
 
 VERSION >= v"1.11.0-DEV.469" && eval(Meta.parse(public_names))
 
@@ -38,6 +40,7 @@ const FloatMat = AbstractMatrix{<:AbstractFloat}
 const SurfacePlot = Surface{Tuple{Matrix{Float64}, Matrix{Float64}, Matrix{Float32}}}
 
 abstract type Phase end
+abstract type Polynomials <: Phase end
 
 struct RadialPolynomial
     λ::Vector{Float64}
@@ -49,7 +52,7 @@ struct Harmonic
     m::Int
 end
 
-struct Polynomial <: Phase
+struct Polynomial <: Polynomials
     inds::NamedTuple{(:j, :n, :m), NTuple{3, Int}}
     N::Float64
     R::RadialPolynomial
@@ -78,6 +81,7 @@ include("Arithmetic.jl")
 include("ZernikePlot.jl")
 include("ScaleAperture.jl")
 include("TransformAperture.jl")
+include("Derivatives.jl")
 include("Docstrings.jl")
 
 function (R::RadialPolynomial)(ρ::Real)
@@ -91,7 +95,7 @@ function (M::Harmonic)(θ::Real)
     m < 0 ? -sin(m * θ) : cos(m * θ)
 end
 
-function (Z::Polynomial)(ρ::Real, θ::Real = 0)
+function (Z::Polynomials)(ρ::Real, θ::Real = 0)
     (; N, R, M) = Z
     N * R(ρ) * M(θ)
 end