diff --git a/src/ModalInterlace.jl b/src/ModalInterlace.jl
index 34d106c..4ca21da 100644
--- a/src/ModalInterlace.jl
+++ b/src/ModalInterlace.jl
@@ -1,7 +1,11 @@
+abstract type AbstractModalInterlace{T} <: AbstractBandedBlockBandedMatrix{T} end
+
+axes(Z::AbstractModalInterlace) = blockedrange.(oneto.(Z.MN))
+
 """
 ModalInterlace
 """
-struct ModalInterlace{T, MMNN<:Tuple} <: AbstractBandedBlockBandedMatrix{T}
+struct ModalInterlace{T, MMNN<:Tuple} <: AbstractModalInterlace{T}
     ops
     MN::MMNN
     bandwidths::NTuple{2,Int}
@@ -10,7 +14,7 @@ end
 ModalInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = ModalInterlace{T,typeof(MN)}(ops, MN, bandwidths)
 ModalInterlace(ops::AbstractVector{<:AbstractMatrix{T}}, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = ModalInterlace{T}(ops, MN, bandwidths)
 
-axes(Z::ModalInterlace) = blockedrange.(oneto.(Z.MN))
+
 
 blockbandwidths(R::ModalInterlace) = R.bandwidths
 subblockbandwidths(::ModalInterlace) = (0,0)
@@ -56,3 +60,69 @@ end
 
 # act like lazy array
 Base.BroadcastStyle(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()
+
+
+
+
+"""
+ShiftModalInterlace
+
+for operators that shift the mode
+"""
+struct ShiftModalInterlace{T, MMNN<:Tuple} <: AbstractModalInterlace{T}
+    ops
+    MN::MMNN
+    bandwidths::NTuple{2,Int}
+    subbandwidth::Int
+end
+
+ShiftModalInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}, subbandwidth::Int) where T = ShiftModalInterlace{T,typeof(MN)}(ops, MN, bandwidths, subbandwidth)
+ShiftModalInterlace(ops::AbstractVector{<:AbstractMatrix{T}}, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}, subbandwidth::Int) where T = ShiftModalInterlace{T}(ops, MN, bandwidths, subbandwidth)
+
+
+
+blockbandwidths(R::ShiftModalInterlace) = R.bandwidths
+subblockbandwidths(R::ShiftModalInterlace) = (-R.subbandwidth,R.subbandwidth)
+
+
+function Base.view(R::ShiftModalInterlace{T}, KJ::Block{2}) where T
+    K,J = KJ.n
+    dat = Matrix{T}(undef,1,J)
+    l,u = blockbandwidths(R)
+    λ = R.subbandwidth
+    if isodd(J-K) && -l ≤ J - K ≤ u
+        sh = (J-K-1)÷2
+        if iseven(K)
+            k = K÷2+1
+            dat[1,1] = R.ops[1][k,k+sh]
+        end
+        for m in range(2-iseven(K); step=2, length=J÷2-max(0,sh))
+            k = K÷2-m÷2+isodd(K)
+            dat[1,m] = dat[1,m+1] = R.ops[m+1][k,k+sh]
+        end
+    else
+        fill!(dat, zero(T))
+    end
+    _BandedMatrix(dat, K, 0, 0)
+end
+
+getindex(R::ShiftModalInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
+
+struct ShiftModalInterlaceLayout <: AbstractBandedBlockBandedLayout end
+struct LazyShiftModalInterlaceLayout <: AbstractLazyBandedBlockBandedLayout end
+
+MemoryLayout(::Type{<:ShiftModalInterlace}) = ShiftModalInterlaceLayout()
+MemoryLayout(::Type{<:ShiftModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyShiftModalInterlaceLayout()
+sublayout(::Union{ShiftModalInterlaceLayout,LazyShiftModalInterlaceLayout}, ::Type{<:NTuple{2,BlockSlice{<:BlockOneTo}}}) = ShiftModalInterlaceLayout()
+
+
+function sub_materialize(::ShiftModalInterlaceLayout, V::AbstractMatrix{T}) where T
+    kr,jr = parentindices(V)
+    KR,JR = kr.block,jr.block
+    M,N = Int(last(KR)), Int(last(JR))
+    R = parent(V)
+    ShiftModalInterlace{T}([R.ops[m][1:(M-m+2)÷2,1:(N-m+2)÷2] for m=1:min(N,M)], (M,N), R.bandwidths)
+end
+
+# act like lazy array
+Base.BroadcastStyle(::Type{<:ShiftModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()
diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
index 98d42d3..c993888 100644
--- a/src/MultivariateOrthogonalPolynomials.jl
+++ b/src/MultivariateOrthogonalPolynomials.jl
@@ -4,7 +4,7 @@ using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Block
       LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
       HarmonicOrthogonalPolynomials
 
-import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size, oneto
+import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size, oneto, conj
 import DomainSets: boundary
 
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
@@ -26,8 +26,8 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        UnitTriangle, UnitDisk,
        JacobiTriangle, TriangleWeight, WeightedTriangle,
        DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
-       Zernike, ZernikeWeight, zerniker, zernikez,
-       PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum,
+       Zernike, ComplexZernike, ZernikeWeight, zerniker, zernikez,
+       PartialDerivative, ComplexDerivative, Laplacian, AbsLaplacianPower, AngularMomentum,
        RadialCoordinate, Weighted, Block
 
 include("ModalInterlace.jl")
diff --git a/src/disk.jl b/src/disk.jl
index d5b00f0..8e6502f 100644
--- a/src/disk.jl
+++ b/src/disk.jl
@@ -101,20 +101,43 @@ function getindex(w::ZernikeWeight, xy::StaticVector{2})
 end
 
 
+abstract type AbstractZernike{T} <: BivariateOrthogonalPolynomial{T} end
+
 """
     Zernike(a, b)
 
 is a quasi-matrix orthogonal `r^(2a) * (1-r^2)^b`
 """
-struct Zernike{T} <: BivariateOrthogonalPolynomial{T}
+struct Zernike{T} <: AbstractZernike{T}
     a::T
     b::T
     Zernike{T}(a::T, b::T) where T = new{T}(a, b)
 end
-Zernike{T}(a, b) where T = Zernike{T}(convert(T,a), convert(T,b))
-Zernike(a::T, b::V) where {T,V} = Zernike{float(promote_type(T,V))}(a, b)
-Zernike{T}(b) where T = Zernike{T}(zero(b), b)
-Zernike{T}() where T = Zernike{T}(zero(T))
+
+
+"""
+    ComplexZernike(a, b)
+
+is a complex-valued quasi-matrix orthogonal `r^(2a) * (1-r^2)^b`
+"""
+struct ComplexZernike{T} <: AbstractZernike{T}
+    a::T
+    b::T
+    ComplexZernike{T}(a::T, b::T) where T = new{T}(a, b)
+end
+
+for Zer in (:Zernike, :ComplexZernike)
+    @eval begin
+        $Zer{T}(a, b) where T = $Zer{T}(convert(T,a), convert(T,b))
+        $Zer(a::T, b::V) where {T,V} = $Zer{float(promote_type(T,V))}(a, b)
+        $Zer{T}(b) where T = $Zer{T}(zero(b), b)
+        $Zer{T}() where T = $Zer{T}(zero(T))
+
+        $Zer() = $Zer{Float64}()
+
+        ==(w::$Zer, v::$Zer) = w.a == v.a && w.b == v.b
+    end
+end
 
 """
     Zernike(b)
@@ -122,15 +145,18 @@ Zernike{T}() where T = Zernike{T}(zero(T))
 is a quasi-matrix orthogonal `(1-r^2)^b`
 """
 Zernike(b) = Zernike(zero(b), b)
-Zernike() = Zernike{Float64}()
 
-axes(P::Zernike{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞)))
+"""
+    ComplexZernike(b)
 
-==(w::Zernike, v::Zernike) = w.a == v.a && w.b == v.b
+is a complex quasi-matrix orthogonal `(1-r^2)^b`
+"""
+ComplexZernike(b) = ComplexZernike(zero(b), b)
 
-copy(A::Zernike) = A
+axes(P::AbstractZernike{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞)))
+copy(A::AbstractZernike) = A
 
-orthogonalityweight(Z::Zernike) = ZernikeWeight(Z.a, Z.b)
+orthogonalityweight(Z::AbstractZernike) = ZernikeWeight(Z.a, Z.b)
 
 zerniker(ℓ, m, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/π) * r^m * normalizedjacobip((ℓ-m) ÷ 2, b, m+a, 2r^2-1)
 zerniker(ℓ, m, b, r) = zerniker(ℓ, m, zero(b), b, r)
@@ -142,9 +168,19 @@ function zernikez(ℓ, ms, a, b, rθ::RadialCoordinate{T}) where T
     zerniker(ℓ, m, a, b, r) * (signbit(ms) ? sin(m*θ) : cos(m*θ))
 end
 
-zernikez(ℓ, ms, a, b, xy::StaticVector{2}) = zernikez(ℓ, ms, a, b, RadialCoordinate(xy))
-zernikez(ℓ, ms, b, xy::StaticVector{2}) = zernikez(ℓ, ms, zero(b), b, xy)
-zernikez(ℓ, ms, xy::StaticVector{2,T}) where T = zernikez(ℓ, ms, zero(T), xy)
+function complexzernikez(ℓ, ms, a, b, rθ::RadialCoordinate{T}) where T
+    r,θ = rθ.r,rθ.θ
+    m = abs(ms)
+    zerniker(ℓ, m, a, b, r) * exp(im*m*θ)
+end
+
+for func in (:zernikez, :complexzernikez)
+    @eval begin
+        $func(ℓ, ms, a, b, xy::StaticVector{2}) = $func(ℓ, ms, a, b, RadialCoordinate(xy))
+        $func(ℓ, ms, b, xy::StaticVector{2}) = $func(ℓ, ms, zero(b), b, xy)
+        $func(ℓ, ms, xy::StaticVector{2,T}) where T = $func(ℓ, ms, zero(T), xy)
+    end
+end
 
 function getindex(Z::Zernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
     ℓ = Int(block(B))-1
@@ -153,10 +189,17 @@ function getindex(Z::Zernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where
     zernikez(ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
 end
 
+function getindex(Z::ComplexZernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
+    ℓ = Int(block(B))-1
+    k = blockindex(B)
+    m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
+    complexzernikez(ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
+end
+
 
-getindex(Z::Zernike, xy::StaticVector{2}, B::BlockIndex{1}) = Z[RadialCoordinate(xy), B]
-getindex(Z::Zernike, xy::StaticVector{2}, B::Block{1}) = [Z[xy, B[j]] for j=1:Int(B)]
-getindex(Z::Zernike, xy::StaticVector{2}, JR::BlockOneTo) = mortar([Z[xy,Block(J)] for J = 1:Int(JR[end])])
+getindex(Z::AbstractZernike, xy::StaticVector{2}, B::BlockIndex{1}) = Z[RadialCoordinate(xy), B]
+getindex(Z::AbstractZernike, xy::StaticVector{2}, B::Block{1}) = [Z[xy, B[j]] for j=1:Int(B)]
+getindex(Z::AbstractZernike, xy::StaticVector{2}, JR::BlockOneTo) = mortar([Z[xy,Block(J)] for J = 1:Int(JR[end])])
 
 
 
@@ -220,7 +263,11 @@ end
 
 factorize(S::FiniteZernike{T}) where T = TransformFactorization(grid(S), ZernikeTransform{T}(blocksize(S,2), parent(S).a, parent(S).b))
 
-# gives the entries for the Laplacian times (1-r^2) * Zernike(1)
+"""
+    WeightedZernikeLaplacianDiag{T}()
+
+gives the entries for the Laplacian times (1-r^2) * Zernike(1)
+"""
 struct WeightedZernikeLaplacianDiag{T} <: AbstractBlockVector{T} end
 
 axes(::WeightedZernikeLaplacianDiag) = (blockedrange(oneto(∞)),)
@@ -243,7 +290,7 @@ end
 
 getindex(W::WeightedZernikeLaplacianDiag, k::Integer) = W[findblockindex(axes(W,1),k)]
 
-@simplify function *(Δ::Laplacian, WZ::Weighted{<:Any,<:Zernike})
+@simplify function *(Δ::Laplacian, WZ::Weighted{<:Any,<:AbstractZernike})
     @assert WZ.P.a == 0 && WZ.P.b == 1
     WZ.P * Diagonal(WeightedZernikeLaplacianDiag{eltype(eltype(WZ))}())
 end
@@ -288,22 +335,55 @@ end
 # 2 dimensional special case, again without the 2^(2*β) factor
 fractionalcfs2d(l::Integer, m::Integer, β) = fractionalcfs(l,m,β,2)
 
-function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
-    TV = promote_type(T,V)
-    A.a == B.a && A.b == B.b && return Eye{TV}(∞)
-    @assert A.a == 0 && A.b == 1
-    @assert B.a == 0 && B.b == 0
-    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(1,0:∞)) .\ Normalized.(Jacobi{TV}.(0,0:∞))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (0,2))
+for Zer in (:Zernike, :ComplexZernike)
+    @eval begin
+        function \(A::$Zer{T}, B::$Zer{V}) where {T,V}
+            TV = promote_type(T,V)
+            A.a == B.a && A.b == B.b && return Eye{TV}(∞)
+            @assert A.a == 0 && A.b == 1
+            @assert B.a == 0 && B.b == 0
+            ModalInterlace{TV}((Normalized.(Jacobi{TV}.(1,0:∞)) .\ Normalized.(Jacobi{TV}.(0,0:∞))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (0,2))
+        end
+
+        function \(A::$Zer{T}, B::Weighted{V,$Zer{V}}) where {T,V}
+            TV = promote_type(T,V)
+            A.a == B.P.a == A.b == B.P.b == 0 && return Eye{TV}(∞)
+            if A.a == A.b == 0
+                @assert B.P.a == 0 && B.P.b == 1
+                ModalInterlace{TV}((Normalized.(Jacobi{TV}.(0, 0:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(1, 0:∞)))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (2,0))
+            else
+                Z = $Zer{TV}() 
+                (A \ Z) * (Z \ B)
+            end
+        end
+    end
 end
 
-function \(A::Zernike{T}, B::Weighted{V,Zernike{V}}) where {T,V}
-    TV = promote_type(T,V)
-    A.a == B.P.a == A.b == B.P.b == 0 && return Eye{TV}(∞)
-    if A.a == A.b == 0
-        @assert B.P.a == 0 && B.P.b == 1
-        ModalInterlace{TV}((Normalized.(Jacobi{TV}.(0, 0:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(1, 0:∞)))) ./ sqrt(convert(TV, 2)), (ℵ₀,ℵ₀), (2,0))
-    else
-        Z = Zernike{TV}() 
-        (A \ Z) * (Z \ B)
+
+#########
+# ComplexDerivative
+# is the complex differential (∂ˣ - im*∂ʸ)/2
+#########
+
+
+for Der in (:ComplexDerivative, :ConjComplexDerivative) 
+    @eval begin
+        struct $Der{T,Ax<:Inclusion} <: LazyQuasiMatrix{Complex{T}}
+            axis::Ax
+        end
+
+        $Der{T}(axis::Inclusion) where {k,T} = $Der{T,typeof(axis)}(axis)
+        $Der{T}(domain) where {k,T} = $Der{T}(Inclusion(domain))
+        $Der(axis) where k = $Der{eltype(eltype(axis))}(axis)
+
+        axes(D::$Der) = (D.axis, D.axis)
+        ==(a::$Der, b::$Der) where k = a.axis == b.axis
+        copy(D::$Der) where k = $Der(copy(D.axis))
+
+        ^(D::$Der, k::Integer) = ApplyQuasiArray(^, D, k)
     end
-end
\ No newline at end of file
+end
+
+conj(D::ComplexDerivative{T}) where T = ConjComplexDerivative{T}(D.axis)
+conj(D::ConjComplexDerivative{T}) where T = ComplexDerivative{T}(D.axis)
+
diff --git a/test/test_disk.jl b/test/test_disk.jl
index ac9528c..8e1d8ae 100644
--- a/test/test_disk.jl
+++ b/test/test_disk.jl
@@ -1,488 +1,501 @@
-using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, RecipesBase, Test, SpecialFunctions
-import MultivariateOrthogonalPolynomials: DiskTrav, grid, ZernikeTransform, ZernikeITransform, *
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, 
+        BandedMatrices, FastTransforms, LinearAlgebra, RecipesBase, Test, SpecialFunctions, ForwardDiff
+import MultivariateOrthogonalPolynomials: DiskTrav, grid, ZernikeTransform, ZernikeITransform, jacobian
 import ClassicalOrthogonalPolynomials: HalfWeighted
 
 
 @testset "Disk" begin
-    @testset "Transform" begin
-        N = 5
-        T = ZernikeTransform{Float64}(N, 0, 0)
-        Ti = ZernikeITransform{Float64}(N, 0, 0)
+    @testset "Real" begin
+        @testset "Transform" begin
+            N = 5
+            T = ZernikeTransform{Float64}(N, 0, 0)
+            Ti = ZernikeITransform{Float64}(N, 0, 0)
 
-        v = PseudoBlockArray(randn(sum(1:N)),1:N)
-        @test T * (Ti * v) ≈ v
-        
+            v = PseudoBlockArray(randn(sum(1:N)),1:N)
+            @test T * (Ti * v) ≈ v
+            
 
-        @test_throws MethodError T * randn(15)
-    end
-    @testset "Basics" begin
-        @test ZernikeWeight(1)[SVector(0.1,0.2)] ≈ (1 - 0.1^2 - 0.2^2)
-        @test ZernikeWeight(1) == ZernikeWeight(1)
-
-        @test Zernike() == Zernike()
-        @test Zernike(1) ≠ Zernike()
-        @test Zernike() ≡ copy(Zernike())
-
-        @test ZernikeWeight() == ZernikeWeight() == ZernikeWeight(0,0) == 
-                ZernikeWeight(0) == ZernikeWeight{Float64}() == 
-                ZernikeWeight{Float64}(0) == ZernikeWeight{Float64}(0, 0)
-        @test ZernikeWeight(1) ≠ ZernikeWeight()
-        @test ZernikeWeight() ≡ copy(ZernikeWeight())
-    end
+            @test_throws MethodError T * randn(15)
+        end
+        @testset "Basics" begin
+            @test ZernikeWeight(1)[SVector(0.1,0.2)] ≈ (1 - 0.1^2 - 0.2^2)
+            @test ZernikeWeight(1) == ZernikeWeight(1)
+
+            @test Zernike() == Zernike()
+            @test Zernike(1) ≠ Zernike()
+            @test Zernike() ≡ copy(Zernike())
+
+            @test ZernikeWeight() == ZernikeWeight() == ZernikeWeight(0,0) == 
+                    ZernikeWeight(0) == ZernikeWeight{Float64}() == 
+                    ZernikeWeight{Float64}(0) == ZernikeWeight{Float64}(0, 0)
+            @test ZernikeWeight(1) ≠ ZernikeWeight()
+            @test ZernikeWeight() ≡ copy(ZernikeWeight())
+        end
 
-    @testset "Evaluation" begin
-        r,θ = 0.1, 0.2
-        rθ = RadialCoordinate(r,θ)
-        xy = SVector(rθ)
-        @test Zernike()[rθ,1] ≈ Zernike()[xy,1] ≈ inv(sqrt(π)) ≈ zernikez(0, 0, rθ)
-        @test Zernike()[rθ,Block(1)] ≈ Zernike()[xy,Block(1)] ≈ [inv(sqrt(π))]
-        @test Zernike()[rθ,Block(2)] ≈ [2r/sqrt(π)*sin(θ), 2r/sqrt(π)*cos(θ)] ≈ [zernikez(1, -1, rθ), zernikez(1, 1, rθ)]
-        @test Zernike()[rθ,Block(3)] ≈ [sqrt(3/π)*(2r^2-1),sqrt(6/π)*r^2*sin(2θ),sqrt(6/π)*r^2*cos(2θ)] ≈ [zernikez(2, 0, rθ), zernikez(2, -2, rθ), zernikez(2, 2, rθ)]
-        @test Zernike()[rθ,Block(4)] ≈ [zernikez(3, -1, rθ), zernikez(3, 1, rθ), zernikez(3, -3, rθ), zernikez(3, 3, rθ)]
-
-        @test zerniker(5, 0, norm(xy)) ≈ zernikez(5, 0, xy)
-    end
+        @testset "Evaluation" begin
+            r,θ = 0.1, 0.2
+            rθ = RadialCoordinate(r,θ)
+            xy = SVector(rθ)
+            @test Zernike()[rθ,1] ≈ Zernike()[xy,1] ≈ inv(sqrt(π)) ≈ zernikez(0, 0, rθ)
+            @test Zernike()[rθ,Block(1)] ≈ Zernike()[xy,Block(1)] ≈ [inv(sqrt(π))]
+            @test Zernike()[rθ,Block(2)] ≈ [2r/sqrt(π)*sin(θ), 2r/sqrt(π)*cos(θ)] ≈ [zernikez(1, -1, rθ), zernikez(1, 1, rθ)]
+            @test Zernike()[rθ,Block(3)] ≈ [sqrt(3/π)*(2r^2-1),sqrt(6/π)*r^2*sin(2θ),sqrt(6/π)*r^2*cos(2θ)] ≈ [zernikez(2, 0, rθ), zernikez(2, -2, rθ), zernikez(2, 2, rθ)]
+            @test Zernike()[rθ,Block(4)] ≈ [zernikez(3, -1, rθ), zernikez(3, 1, rθ), zernikez(3, -3, rθ), zernikez(3, 3, rθ)]
+
+            @test zerniker(5, 0, norm(xy)) ≈ zernikez(5, 0, xy)
+        end
 
-    @testset "DiskTrav" begin
-        @test DiskTrav(reshape([1],1,1)) == [1]
-        @test DiskTrav([1 2 3]) == 1:3
-        @test DiskTrav([1 2 3 5 6;
-                        4 0 0 0 0]) == 1:6
-        @test DiskTrav([1 2 3 5 6 9 10;
-                        4 7 8 0 0 0  0]) == 1:10
-
-        @test DiskTrav([1 2 3 5 6 9 10; 4 7 8 0 0 0  0])[Block(3)] == 4:6
-
-        @test_throws ArgumentError DiskTrav([1 2])
-        @test_throws ArgumentError DiskTrav([1 2 3 4])
-        @test_throws ArgumentError DiskTrav([1 2 3; 4 5 6])
-        @test_throws ArgumentError DiskTrav([1 2 3 4; 5 6 7 8])
-
-        for N = 1:10
-            v = PseudoBlockArray(1:sum(1:N),1:N)
-            if iseven(N)
-                @test DiskTrav(v) == [v; zeros(N+1)]
-            else
-                @test DiskTrav(v) == v
+        @testset "DiskTrav" begin
+            @test DiskTrav(reshape([1],1,1)) == [1]
+            @test DiskTrav([1 2 3]) == 1:3
+            @test DiskTrav([1 2 3 5 6;
+                            4 0 0 0 0]) == 1:6
+            @test DiskTrav([1 2 3 5 6 9 10;
+                            4 7 8 0 0 0  0]) == 1:10
+
+            @test DiskTrav([1 2 3 5 6 9 10; 4 7 8 0 0 0  0])[Block(3)] == 4:6
+
+            @test_throws ArgumentError DiskTrav([1 2])
+            @test_throws ArgumentError DiskTrav([1 2 3 4])
+            @test_throws ArgumentError DiskTrav([1 2 3; 4 5 6])
+            @test_throws ArgumentError DiskTrav([1 2 3 4; 5 6 7 8])
+
+            for N = 1:10
+                v = PseudoBlockArray(1:sum(1:N),1:N)
+                if iseven(N)
+                    @test DiskTrav(v) == [v; zeros(N+1)]
+                else
+                    @test DiskTrav(v) == v
+                end
             end
         end
-    end
 
-    @testset "Transform" begin
-        for (a,b) in ((0,0), (0.1, 0.2), (0,1))
-            Zn = Zernike(a,b)[:,Block.(Base.OneTo(3))]
-            for k = 1:6
-                @test factorize(Zn) \ Zernike(a,b)[:,k] ≈ [zeros(k-1); 1; zeros(6-k)]
+        @testset "Transform" begin
+            for (a,b) in ((0,0), (0.1, 0.2), (0,1))
+                Zn = Zernike(a,b)[:,Block.(Base.OneTo(3))]
+                for k = 1:6
+                    @test factorize(Zn) \ Zernike(a,b)[:,k] ≈ [zeros(k-1); 1; zeros(6-k)]
+                end
+
+                Z = Zernike(a,b);
+                xy = axes(Z,1); x,y = first.(xy),last.(xy);
+                u = Z * (Z \ exp.(x .* cos.(y)))
+                @test u[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2))
             end
-
-            Z = Zernike(a,b);
-            xy = axes(Z,1); x,y = first.(xy),last.(xy);
-            u = Z * (Z \ exp.(x .* cos.(y)))
-            @test u[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2))
         end
-    end
-
-    @testset "Laplacian" begin
-        # u = r^m*f(r^2) * cos(m*θ)
-        # u_r = (m*r^(m-1)*f(r^2) + 2r^(m+1)*f'(r^2)) * cos(m*θ)
-        # u_rr = (m*(m-1)*r^(m-2)*f(r^2) + (4m+2)*r^m*f'(r^2) + 2r^(m+1)*f''(r^2)) * cos(m*θ)
-        # u_rr + u_r/r + u_θθ/r^2 = (4*(m+1)*f'(r^2) + 2r*f''(r^2)) * r^m * cos(m*θ)
-        # t = r^2, dt = 2r * dr, 4*(m+1)*f'(t) + 2sqrt(t)*f''(t) = 4 t^(-m) * d/dt * t^(m+1) f'(t)
-        # d/ds * (1-s) * P_n^(1,m)(s) = -n*P_n^(0,m+1)(s)
-        # use L^6 and L^6'
-        # 2t-1 = s, 2dt = ds
-
-
-        ℓ, m, b = 6, 2, 1
-        x,y = 0.1,0.2
-        r = sqrt(x^2+y^2)
-        θ = atan(y,x)
-        t = r^2
-        u = xy -> zernikez(ℓ, m, b, xy)
-        ur = r -> zerniker(ℓ, m, b, r)
-
-        f = t -> sqrt(2^(m+b+2-iszero(m))/π) * normalizedjacobip((ℓ-m) ÷ 2, b, m, 2t-1)
-        ur = r -> r^m*f(r^2)
-        @test ur(r) ≈ zerniker(ℓ, m, b, r)
-        @test f(r^2) ≈ r^(-m) * zerniker(ℓ, m, b, r)
-        # u = xy -> ((x,y) = xy; ur(norm(xy)) * cos(m*atan(y,x)))
-        # t = r^2; 4*(m+1)*derivative(f,t) + 4t*derivative2(f,t)
-        
-        # @test derivative(ur,r) ≈  m*r^(m-1)*f(r^2) + 2r^(m+1)*derivative(f,r^2)
-        # @test derivative2(ur,r) ≈ m*(m-1)*r^(m-2)*f(r^2) + (4m+2)*r^m * derivative(f,r^2) + 4r^(m+2)*derivative2(f,r^2)
-        # @test lapr(ur, m, r) ≈ 4*((m+1) * derivative(f,r^2) + r^2*derivative2(f,r^2)) * r^m
-        # @test lapr(ur, m, r) ≈ 4*((m+1) * derivative(f,t) + t*derivative2(f,t)) * t^(m/2)
-        
-
-        ℓ, m, b = 1, 1, 1
-
-        f = t -> sqrt(2^(m+b+2-iszero(m))/π) * (1-t) * normalizedjacobip((ℓ-m) ÷ 2, b, m, 2t-1)
-        ur = r -> r^m*f(r^2)
-        @test ur(r) ≈ (1-r^2) * zerniker(ℓ, m, b, r)
-
-        D = Derivative(axes(Chebyshev(),1))
-        D1 = Normalized(Jacobi(0, m+1)) \ (D * (HalfWeighted{:a}(Normalized(Jacobi(1, m)))))
-        D2 = HalfWeighted{:b}(Normalized(Jacobi(1, m))) \ (D * (HalfWeighted{:b}(Normalized(Jacobi(0, m+1)))))
-
-        @test (D1 * D2)[band(0)][1:10] ≈ -((1:∞) .* ((1+m):∞))[1:10]
-
-                
-        ℓ = m = 0; b= 1
-        d = -4*((1:∞) .* ((m+1):∞))
-        xy = SVector(0.1,0.2)
-        # ℓ = m = 0; b= 1
-        # u = xy -> (1 - norm(xy)^2) * zernikez(0, 0, 1, xy)
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,1] * (-4)
-
-        # u = xy -> (1 - norm(xy)^2) * zernikez(1 , -1, 1, xy)
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,2] * (-4) * 1 * 2
-        # u = xy -> (1 - norm(xy)^2) * zernikez(1 , 1, 1, xy)
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,3] * (-4) * 1 * 2
-
-        # u = xy -> (1 - norm(xy)^2) * zernikez(2 , 0, 1, xy) # eval at 2
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,4] * (-4) * 2^2
-        # u = xy -> (1 - norm(xy)^2) * zernikez(2 , -2, 1, xy) # eval at 1
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,5] * (-4) * 3
-        # u = xy -> (1 - norm(xy)^2) * zernikez(2 , 2, 1, xy) # eval at 1
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,6] * (-4) * 3
-
-        # u = xy -> (1 - norm(xy)^2) * zernikez(3 , -1, 1, xy) # eval at 2
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,7] * (-4) * 2 * 3
-        # u = xy -> (1 - norm(xy)^2) * zernikez(3 , 1, 1, xy)
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,8] * (-4) * 2 * 3
-        # u = xy -> (1 - norm(xy)^2) * zernikez(3 , -3, 1, xy) # eval at 1
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,9] * (-4) * 4 * 1
-        # u = xy -> (1 - norm(xy)^2) * zernikez(3 , 3, 1, xy) # eval at 1
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,10] * (-4) * 4 * 1
-
-        # u = xy -> (1 - norm(xy)^2) * zernikez(4 , 0, 1, xy) # eval at 3
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,11] * (-4) * 3^2
-        # u = xy -> (1 - norm(xy)^2) * zernikez(4 , -2, 1, xy) # eval at 2
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,12] * (-4) * 4 * 2
-        # u = xy -> (1 - norm(xy)^2) * zernikez(4 , 2, 1, xy) # eval at 2
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,13] * (-4) * 4 * 2
-        # u = xy -> (1 - norm(xy)^2) * zernikez(4 , -4, 1, xy) # eval at 1
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,14] * (-4) * 5 * 1
-        # u = xy -> (1 - norm(xy)^2) * zernikez(4 , 4, 1, xy) # eval at 1
-        # @test lap(u, xy...) ≈ Zernike(1)[xy,15] * (-4) * 5 * 1
-
-        WZ = Weighted(Zernike(1)) # Zernike(1) weighted by (1-r^2)
-        Δ = Laplacian(axes(WZ,1))
-        Δ_Z = Zernike(1) \ (Δ * WZ)
-        @test exp.(Δ_Z)[1:10,1:10] == exp.(Δ_Z[1:10,1:10])
-
-        xy = axes(WZ,1)
-        x,y = first.(xy),last.(xy)
-        u = @. (1 - x^2 - y^2) * exp(x*cos(y))
-        Δu = @. (-exp(x*cos(y)) * (4 - x*(-5 + x^2 + y^2)cos(y) + (-1 + x^2 + y^2)cos(y)^2 - 4x*y*sin(y) + x^2*(x^2 + y^2-1)*sin(y)^2))
-        @test (WZ * (WZ \ u))[SVector(0.1,0.2)] ≈ u[SVector(0.1,0.2)]
-        @test (Δ_Z * (WZ \ u))[1:100]  ≈ (Zernike(1) \ Δu)[1:100] 
-    end
 
-    @testset "Conversion" begin
-        R0 = Normalized(Jacobi(1, 0)) \ Normalized(Jacobi(0, 0))
-        R1 = Normalized(Jacobi(1, 1)) \ Normalized(Jacobi(0, 1))
-        R2 = Normalized(Jacobi(1, 2)) \ Normalized(Jacobi(0, 2))
-        R3 = Normalized(Jacobi(1, 3)) \ Normalized(Jacobi(0, 3))
+        @testset "Laplacian" begin
+            # u = r^m*f(r^2) * cos(m*θ)
+            # u_r = (m*r^(m-1)*f(r^2) + 2r^(m+1)*f'(r^2)) * cos(m*θ)
+            # u_rr = (m*(m-1)*r^(m-2)*f(r^2) + (4m+2)*r^m*f'(r^2) + 2r^(m+1)*f''(r^2)) * cos(m*θ)
+            # u_rr + u_r/r + u_θθ/r^2 = (4*(m+1)*f'(r^2) + 2r*f''(r^2)) * r^m * cos(m*θ)
+            # t = r^2, dt = 2r * dr, 4*(m+1)*f'(t) + 2sqrt(t)*f''(t) = 4 t^(-m) * d/dt * t^(m+1) f'(t)
+            # d/ds * (1-s) * P_n^(1,m)(s) = -n*P_n^(0,m+1)(s)
+            # use L^6 and L^6'
+            # 2t-1 = s, 2dt = ds
+
+
+            ℓ, m, b = 6, 2, 1
+            x,y = 0.1,0.2
+            r = sqrt(x^2+y^2)
+            θ = atan(y,x)
+            t = r^2
+            u = xy -> zernikez(ℓ, m, b, xy)
+            ur = r -> zerniker(ℓ, m, b, r)
+
+            f = t -> sqrt(2^(m+b+2-iszero(m))/π) * normalizedjacobip((ℓ-m) ÷ 2, b, m, 2t-1)
+            ur = r -> r^m*f(r^2)
+            @test ur(r) ≈ zerniker(ℓ, m, b, r)
+            @test f(r^2) ≈ r^(-m) * zerniker(ℓ, m, b, r)
+            # u = xy -> ((x,y) = xy; ur(norm(xy)) * cos(m*atan(y,x)))
+            # t = r^2; 4*(m+1)*derivative(f,t) + 4t*derivative2(f,t)
+            
+            # @test derivative(ur,r) ≈  m*r^(m-1)*f(r^2) + 2r^(m+1)*derivative(f,r^2)
+            # @test derivative2(ur,r) ≈ m*(m-1)*r^(m-2)*f(r^2) + (4m+2)*r^m * derivative(f,r^2) + 4r^(m+2)*derivative2(f,r^2)
+            # @test lapr(ur, m, r) ≈ 4*((m+1) * derivative(f,r^2) + r^2*derivative2(f,r^2)) * r^m
+            # @test lapr(ur, m, r) ≈ 4*((m+1) * derivative(f,t) + t*derivative2(f,t)) * t^(m/2)
+            
+
+            ℓ, m, b = 1, 1, 1
+
+            f = t -> sqrt(2^(m+b+2-iszero(m))/π) * (1-t) * normalizedjacobip((ℓ-m) ÷ 2, b, m, 2t-1)
+            ur = r -> r^m*f(r^2)
+            @test ur(r) ≈ (1-r^2) * zerniker(ℓ, m, b, r)
+
+            D = Derivative(axes(Chebyshev(),1))
+            D1 = Normalized(Jacobi(0, m+1)) \ (D * (HalfWeighted{:a}(Normalized(Jacobi(1, m)))))
+            D2 = HalfWeighted{:b}(Normalized(Jacobi(1, m))) \ (D * (HalfWeighted{:b}(Normalized(Jacobi(0, m+1)))))
+
+            @test (D1 * D2)[band(0)][1:10] ≈ -((1:∞) .* ((1+m):∞))[1:10]
+
+                    
+            ℓ = m = 0; b= 1
+            d = -4*((1:∞) .* ((m+1):∞))
+            xy = SVector(0.1,0.2)
+            # ℓ = m = 0; b= 1
+            # u = xy -> (1 - norm(xy)^2) * zernikez(0, 0, 1, xy)
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,1] * (-4)
+
+            # u = xy -> (1 - norm(xy)^2) * zernikez(1 , -1, 1, xy)
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,2] * (-4) * 1 * 2
+            # u = xy -> (1 - norm(xy)^2) * zernikez(1 , 1, 1, xy)
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,3] * (-4) * 1 * 2
+
+            # u = xy -> (1 - norm(xy)^2) * zernikez(2 , 0, 1, xy) # eval at 2
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,4] * (-4) * 2^2
+            # u = xy -> (1 - norm(xy)^2) * zernikez(2 , -2, 1, xy) # eval at 1
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,5] * (-4) * 3
+            # u = xy -> (1 - norm(xy)^2) * zernikez(2 , 2, 1, xy) # eval at 1
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,6] * (-4) * 3
+
+            # u = xy -> (1 - norm(xy)^2) * zernikez(3 , -1, 1, xy) # eval at 2
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,7] * (-4) * 2 * 3
+            # u = xy -> (1 - norm(xy)^2) * zernikez(3 , 1, 1, xy)
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,8] * (-4) * 2 * 3
+            # u = xy -> (1 - norm(xy)^2) * zernikez(3 , -3, 1, xy) # eval at 1
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,9] * (-4) * 4 * 1
+            # u = xy -> (1 - norm(xy)^2) * zernikez(3 , 3, 1, xy) # eval at 1
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,10] * (-4) * 4 * 1
+
+            # u = xy -> (1 - norm(xy)^2) * zernikez(4 , 0, 1, xy) # eval at 3
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,11] * (-4) * 3^2
+            # u = xy -> (1 - norm(xy)^2) * zernikez(4 , -2, 1, xy) # eval at 2
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,12] * (-4) * 4 * 2
+            # u = xy -> (1 - norm(xy)^2) * zernikez(4 , 2, 1, xy) # eval at 2
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,13] * (-4) * 4 * 2
+            # u = xy -> (1 - norm(xy)^2) * zernikez(4 , -4, 1, xy) # eval at 1
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,14] * (-4) * 5 * 1
+            # u = xy -> (1 - norm(xy)^2) * zernikez(4 , 4, 1, xy) # eval at 1
+            # @test lap(u, xy...) ≈ Zernike(1)[xy,15] * (-4) * 5 * 1
+
+            WZ = Weighted(Zernike(1)) # Zernike(1) weighted by (1-r^2)
+            Δ = Laplacian(axes(WZ,1))
+            Δ_Z = Zernike(1) \ (Δ * WZ)
+            @test exp.(Δ_Z)[1:10,1:10] == exp.(Δ_Z[1:10,1:10])
 
-        xy = SVector(0.1,0.2)
-        @test Zernike()[xy,Block(1)[1]] ≈ Zernike(1)[xy,Block(1)[1]]/sqrt(2)
+            xy = axes(WZ,1)
+            x,y = first.(xy),last.(xy)
+            u = @. (1 - x^2 - y^2) * exp(x*cos(y))
+            Δu = @. (-exp(x*cos(y)) * (4 - x*(-5 + x^2 + y^2)cos(y) + (-1 + x^2 + y^2)cos(y)^2 - 4x*y*sin(y) + x^2*(x^2 + y^2-1)*sin(y)^2))
+            @test (WZ * (WZ \ u))[SVector(0.1,0.2)] ≈ u[SVector(0.1,0.2)]
+            @test (Δ_Z * (WZ \ u))[1:100]  ≈ (Zernike(1) \ Δu)[1:100] 
+        end
 
-        @test Zernike()[xy,Block(2)[1]] ≈ Zernike(1)[xy,Block(2)[1]]*R1[1,1]/sqrt(2)
-        @test Zernike()[xy,Block(2)[2]] ≈ Zernike(1)[xy,Block(2)[2]]*R1[1,1]/sqrt(2)
+        @testset "Conversion" begin
+            R0 = Normalized(Jacobi(1, 0)) \ Normalized(Jacobi(0, 0))
+            R1 = Normalized(Jacobi(1, 1)) \ Normalized(Jacobi(0, 1))
+            R2 = Normalized(Jacobi(1, 2)) \ Normalized(Jacobi(0, 2))
+            R3 = Normalized(Jacobi(1, 3)) \ Normalized(Jacobi(0, 3))
 
-        @test Zernike()[xy,Block(3)[1]] ≈ R0[1:2,2]'*Zernike(1)[xy,getindex.(Block.(1:2:3),1)]/sqrt(2)
-        @test Zernike()[xy,Block(3)[2]] ≈ R2[1,1]*Zernike(1)[xy,Block(3)[2]]/sqrt(2)
-        @test Zernike()[xy,Block(3)[3]] ≈ R2[1,1]*Zernike(1)[xy,Block(3)[3]]/sqrt(2)
+            xy = SVector(0.1,0.2)
+            @test Zernike()[xy,Block(1)[1]] ≈ Zernike(1)[xy,Block(1)[1]]/sqrt(2)
 
-        @test Zernike()[xy,Block(4)[1]] ≈ R1[1:2,2]'*Zernike(1)[xy,getindex.(Block.(2:2:4),1)]/sqrt(2)
-        @test Zernike()[xy,Block(4)[2]] ≈ R1[1:2,2]'*Zernike(1)[xy,getindex.(Block.(2:2:4),2)]/sqrt(2)
-        @test Zernike()[xy,Block(4)[3]] ≈ R3[1,1]*Zernike(1)[xy,Block(4)[3]]/sqrt(2)
-        @test Zernike()[xy,Block(4)[4]] ≈ R3[1,1]*Zernike(1)[xy,Block(4)[4]]/sqrt(2)
+            @test Zernike()[xy,Block(2)[1]] ≈ Zernike(1)[xy,Block(2)[1]]*R1[1,1]/sqrt(2)
+            @test Zernike()[xy,Block(2)[2]] ≈ Zernike(1)[xy,Block(2)[2]]*R1[1,1]/sqrt(2)
 
-        @test Zernike()[xy,Block(5)[1]] ≈ R0[2:3,3]'*Zernike(1)[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
+            @test Zernike()[xy,Block(3)[1]] ≈ R0[1:2,2]'*Zernike(1)[xy,getindex.(Block.(1:2:3),1)]/sqrt(2)
+            @test Zernike()[xy,Block(3)[2]] ≈ R2[1,1]*Zernike(1)[xy,Block(3)[2]]/sqrt(2)
+            @test Zernike()[xy,Block(3)[3]] ≈ R2[1,1]*Zernike(1)[xy,Block(3)[3]]/sqrt(2)
 
+            @test Zernike()[xy,Block(4)[1]] ≈ R1[1:2,2]'*Zernike(1)[xy,getindex.(Block.(2:2:4),1)]/sqrt(2)
+            @test Zernike()[xy,Block(4)[2]] ≈ R1[1:2,2]'*Zernike(1)[xy,getindex.(Block.(2:2:4),2)]/sqrt(2)
+            @test Zernike()[xy,Block(4)[3]] ≈ R3[1,1]*Zernike(1)[xy,Block(4)[3]]/sqrt(2)
+            @test Zernike()[xy,Block(4)[4]] ≈ R3[1,1]*Zernike(1)[xy,Block(4)[4]]/sqrt(2)
 
-        R = Zernike(1) \ Zernike()
+            @test Zernike()[xy,Block(5)[1]] ≈ R0[2:3,3]'*Zernike(1)[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
 
-        @test R[Block.(Base.OneTo(6)), Block.(Base.OneTo(7))] == R[Block.(1:6), Block.(1:7)]
-        @test Zernike()[xy,Block.(1:6)]' ≈ Zernike(1)[xy,Block.(1:6)]'*R[Block.(1:6),Block.(1:6)]
-    end
 
-    @testset "Lowering" begin
-        L0 = Normalized(Jacobi(0, 0)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 0)))
-        L1 = Normalized(Jacobi(0, 1)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 1)))
-        L2 = Normalized(Jacobi(0, 2)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 2)))
-        L3 = Normalized(Jacobi(0, 3)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 3)))
+            R = Zernike(1) \ Zernike()
 
-        xy = SVector(0.1,0.2)
-        r = norm(xy)
-        w = 1 - r^2
+            @test R[Block.(Base.OneTo(6)), Block.(Base.OneTo(7))] == R[Block.(1:6), Block.(1:7)]
+            @test Zernike()[xy,Block.(1:6)]' ≈ Zernike(1)[xy,Block.(1:6)]'*R[Block.(1:6),Block.(1:6)]
+        end
 
-        @test w*Zernike(1)[xy,Block(1)[1]] ≈ L0[1:2,1]'*Zernike()[xy,getindex.(Block.(1:2:3),1)] / sqrt(2)
+        @testset "Lowering" begin
+            L0 = Normalized(Jacobi(0, 0)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 0)))
+            L1 = Normalized(Jacobi(0, 1)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 1)))
+            L2 = Normalized(Jacobi(0, 2)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 2)))
+            L3 = Normalized(Jacobi(0, 3)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 3)))
 
-        @test w*Zernike(1)[xy,Block(2)[1]] ≈ L1[1:2,1]'*Zernike()[xy,getindex.(Block.(2:2:4),1)]/sqrt(2)
-        @test w*Zernike(1)[xy,Block(2)[2]] ≈ L1[1:2,1]'*Zernike()[xy,getindex.(Block.(2:2:4),2)]/sqrt(2)
+            xy = SVector(0.1,0.2)
+            r = norm(xy)
+            w = 1 - r^2
 
-        @test w*Zernike(1)[xy,Block(3)[1]] ≈ L0[2:3,2]'*Zernike()[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
-        @test w*Zernike(1)[xy,Block(3)[2]] ≈ L2[1:2,1]'*Zernike()[xy,getindex.(Block.(3:2:5),2)]/sqrt(2)
-        @test w*Zernike(1)[xy,Block(3)[3]] ≈ L2[1:2,1]'*Zernike()[xy,getindex.(Block.(3:2:5),3)]/sqrt(2)
+            @test w*Zernike(1)[xy,Block(1)[1]] ≈ L0[1:2,1]'*Zernike()[xy,getindex.(Block.(1:2:3),1)] / sqrt(2)
 
-        L = Zernike() \ Weighted(Zernike(1))
-        @test w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike()[xy,Block.(1:7)]'*L[Block.(1:7),Block.(1:5)] 
+            @test w*Zernike(1)[xy,Block(2)[1]] ≈ L1[1:2,1]'*Zernike()[xy,getindex.(Block.(2:2:4),1)]/sqrt(2)
+            @test w*Zernike(1)[xy,Block(2)[2]] ≈ L1[1:2,1]'*Zernike()[xy,getindex.(Block.(2:2:4),2)]/sqrt(2)
 
-        @test exp.(L)[1:10,1:10] == exp.(L[1:10,1:10])
+            @test w*Zernike(1)[xy,Block(3)[1]] ≈ L0[2:3,2]'*Zernike()[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
+            @test w*Zernike(1)[xy,Block(3)[2]] ≈ L2[1:2,1]'*Zernike()[xy,getindex.(Block.(3:2:5),2)]/sqrt(2)
+            @test w*Zernike(1)[xy,Block(3)[3]] ≈ L2[1:2,1]'*Zernike()[xy,getindex.(Block.(3:2:5),3)]/sqrt(2)
 
-        L2 = Zernike(1) \ Weighted(Zernike(1))
-        @test w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike(1)[xy,Block.(1:7)]'*L2[Block.(1:7),Block.(1:5)] 
-    end
+            L = Zernike() \ Weighted(Zernike(1))
+            @test w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike()[xy,Block.(1:7)]'*L[Block.(1:7),Block.(1:5)] 
 
-    @testset "plotting" begin
-        Z = Zernike()
-        u = Z * [1; 2; zeros(∞)];
-        rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
-        g = MultivariateOrthogonalPolynomials.plotgrid(Z[:,1:3])
-        @test rep[1].args == (first.(g),last.(g),u[g])
-    end
-end
-
-@testset "Fractional Laplacian on Unit Disk" begin
-    @testset "Fractional Laplacian on Disk: (-Δ)^(β) == -Δ when β=1" begin
-        WZ = Weighted(Zernike(1.))
-        Δ = Laplacian(axes(WZ,1))
-        Δ_Z = Zernike(1) \ (Δ * WZ)
-        Δfrac = AbsLaplacianPower(axes(WZ,1),1.)
-        Δ_Zfrac = Zernike(1) \ (Δfrac * WZ)
-        @test Δ_Z[1:100,1:100] ≈ -Δ_Zfrac[1:100,1:100]
-    end
+            @test exp.(L)[1:10,1:10] == exp.(L[1:10,1:10])
 
-    @testset "Fractional Laplacian on Disk: Computing f where (-Δ)^(β) u = f" begin
-        @testset "Set 1 - Explicitly known constant f" begin
-            # set up basis
-            β = 1.34
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^β
-            # explicit and computed solutions
-            fexplicit0(d,α) = 2^α*gamma(α/2+1)*gamma((d+α)/2)/gamma(d/2) # note that here, α = 2*β
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit0(2,2*β) ≈ f[(0.1,0.4)] ≈ f[(0.1137,0.001893)] ≈ f[(0.3721,0.3333)]
-
-            # again for different β
-            β = 2.11
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^β
-            # computed solution
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit0(2,2*β) ≈ f[(0.14,0.41)] ≈ f[(0.1731,0.091893)] ≈ f[(0.3791,0.333333)]
-
-            # again for different β
-            β = 3.14159
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^β
-            # computed solution
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit0(2,2*β) ≈ f[(0.14,0.41)] ≈ f[(0.1837,0.101893)] ≈ f[(0.37222,0.2222)]
-        end
-        @testset "Set 2 - Explicitly known radially symmetric f" begin
-            β = 1.1
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^(β+1)
-            # explicit and computed solutions
-            fexplicit1(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2)/gamma(d/2)*(1-(1+α/d)*norm(x)^2) # α = 2*β
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit1(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
-            @test fexplicit1(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
-            @test fexplicit1(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
-
-            # again for different β
-            β = 2.71999
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^(β+1)
-            # explicit and computed solutions
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit1(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
-            @test fexplicit1(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
-            @test fexplicit1(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+            L2 = Zernike(1) \ Weighted(Zernike(1))
+            @test w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike(1)[xy,Block.(1:7)]'*L2[Block.(1:7),Block.(1:5)] 
         end
-        @testset "Set 3 - Explicitly known f, not radially symmetric" begin
-            # dependence on x
-            β = 2.71999
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^(β)*x
-            # explicit and computed solutions
-            fexplicit2(d,α,x) = 2^α*gamma(α/2+1)*gamma((d+α)/2+1)/gamma(d/2+1)*x[1] # α = 2*β
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit2(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
-            @test fexplicit2(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
-            @test fexplicit2(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
-
-            # different β, dependence on y
-            β = 1.91239
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^(β)*y
-            # explicit and computed solutions
-            fexplicit3(d,α,x) = 2^α*gamma(α/2+1)*gamma((d+α)/2+1)/gamma(d/2+1)*x[2] # α = 2*β
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit3(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
-            @test fexplicit3(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
-            @test fexplicit3(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
-        end
-        @testset "Set 4 - Explicitly known f, different non-radially-symmetric example" begin
-            # dependence on x
-            β = 1.21999
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^(β+1)*x
-            # explicit and computed solutions
-            fexplicit4(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[1] # α = 2*β
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit4(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
-            @test fexplicit4(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
-            @test fexplicit4(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
-
-            # different β, dependence on y
-            β = 0.141
-            Z = Zernike(β)
-            WZ = Weighted(Z)
-            xy = axes(WZ,1)
-            x,y = first.(xy),last.(xy)
-            # generate fractional Laplacian
-            Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-            Δ_Zfrac = Z \ (Δfrac * WZ)
-            # define function whose fractional Laplacian is known
-            u = @. (1 - x^2 - y^2).^(β+1)*y
-            # explicit and computed solutions
-            fexplicit5(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[2] # α = 2*β
-            f = Z*(Δ_Zfrac*(WZ \ u))
-            # compare
-            @test fexplicit5(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
-            @test fexplicit5(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
-            @test fexplicit5(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+
+        @testset "plotting" begin
+            Z = Zernike()
+            u = Z * [1; 2; zeros(∞)];
+            rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
+            g = MultivariateOrthogonalPolynomials.plotgrid(Z[:,1:3])
+            @test rep[1].args == (first.(g),last.(g),u[g])
         end
 
-        @testset "Fractional Poisson equation on Disk: Comparison with explicitly known solutions" begin
-            @testset "Set 1 - Radially symmetric solution" begin
-                # define basis
-                β = 1.1812
-                Z = Zernike(β)
-                WZ = Weighted(Z)
-                xy = axes(WZ,1)
-                x,y = first.(xy),last.(xy)
-                # generate fractional Laplacian
-                Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-                Δ_Zfrac = Z \ (Δfrac * WZ)
-                # define function whose fractional Laplacian is known
-                uexplicit = @. (1 - x^2 - y^2).^(β+1)
-                uexplicitcfs = WZ \ uexplicit
-                # RHS
-                RHS(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2)/gamma(d/2)*(1-(1+α/d)*norm(x)^2) # α = 2*β
-                RHScfs = Z \ @. RHS.(2,2*β,xy)
-                # compute solution
-                ucomputed = Δ_Zfrac \ RHScfs
-                @test uexplicitcfs[1:100] ≈ ucomputed[1:100]
+        @testset "Fractional Laplacian on Unit Disk" begin
+            @testset "Fractional Laplacian on Disk: (-Δ)^(β) == -Δ when β=1" begin
+                WZ = Weighted(Zernike(1.))
+                Δ = Laplacian(axes(WZ,1))
+                Δ_Z = Zernike(1) \ (Δ * WZ)
+                Δfrac = AbsLaplacianPower(axes(WZ,1),1.)
+                Δ_Zfrac = Zernike(1) \ (Δfrac * WZ)
+                @test Δ_Z[1:100,1:100] ≈ -Δ_Zfrac[1:100,1:100]
             end
-            @testset "Set 2 - Non-radially-symmetric solutions" begin
-                # dependence on y
-                β = 0.98812
-                Z = Zernike(β)
-                WZ = Weighted(Z)
-                xy = axes(WZ,1)
-                x,y = first.(xy),last.(xy)
-                # generate fractional Laplacian
-                Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-                Δ_Zfrac = Z \ (Δfrac * WZ)
-                # define function whose fractional Laplacian is known
-                uexplicit = @. (1 - x^2 - y^2).^(β+1)*y
-                uexplicitcfs = WZ \ uexplicit
-                # RHS
-                RHS2(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[2] # α = 2*β
-                RHS2cfs = Z \ @. RHS2.(2,2*β,xy)
-                # compute solution
-                ucomputed = Δ_Zfrac \ RHS2cfs
-                @test uexplicitcfs[1:100] ≈ ucomputed[1:100]
-        
-                # different β, dependence on x
-                β = 0.506
-                Z = Zernike(β)
-                WZ = Weighted(Z)
-                xy = axes(WZ,1)
-                x,y = first.(xy),last.(xy)
-                # generate fractional Laplacian
-                Δfrac = AbsLaplacianPower(axes(WZ,1),β)
-                Δ_Zfrac = Z \ (Δfrac * WZ)
-                # define function whose fractional Laplacian is known
-                uexplicit = @. (1 - x^2 - y^2).^(β+1)*x
-                uexplicitcfs = WZ \ uexplicit
-                # RHS
-                RHS3(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[1] # α = 2*β
-                RHS3cfs = Z \ @. RHS3.(2,2*β,xy)
-                # compute solution
-                ucomputed = Δ_Zfrac \ RHS3cfs
-                @test uexplicitcfs[1:100] ≈ ucomputed[1:100]
+
+            @testset "Fractional Laplacian on Disk: Computing f where (-Δ)^(β) u = f" begin
+                @testset "Set 1 - Explicitly known constant f" begin
+                    # set up basis
+                    β = 1.34
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^β
+                    # explicit and computed solutions
+                    fexplicit0(d,α) = 2^α*gamma(α/2+1)*gamma((d+α)/2)/gamma(d/2) # note that here, α = 2*β
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit0(2,2*β) ≈ f[(0.1,0.4)] ≈ f[(0.1137,0.001893)] ≈ f[(0.3721,0.3333)]
+
+                    # again for different β
+                    β = 2.11
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^β
+                    # computed solution
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit0(2,2*β) ≈ f[(0.14,0.41)] ≈ f[(0.1731,0.091893)] ≈ f[(0.3791,0.333333)]
+
+                    # again for different β
+                    β = 3.14159
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^β
+                    # computed solution
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit0(2,2*β) ≈ f[(0.14,0.41)] ≈ f[(0.1837,0.101893)] ≈ f[(0.37222,0.2222)]
+                end
+                @testset "Set 2 - Explicitly known radially symmetric f" begin
+                    β = 1.1
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^(β+1)
+                    # explicit and computed solutions
+                    fexplicit1(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2)/gamma(d/2)*(1-(1+α/d)*norm(x)^2) # α = 2*β
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit1(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
+                    @test fexplicit1(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
+                    @test fexplicit1(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+
+                    # again for different β
+                    β = 2.71999
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^(β+1)
+                    # explicit and computed solutions
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit1(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
+                    @test fexplicit1(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
+                    @test fexplicit1(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+                end
+                @testset "Set 3 - Explicitly known f, not radially symmetric" begin
+                    # dependence on x
+                    β = 2.71999
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^(β)*x
+                    # explicit and computed solutions
+                    fexplicit2(d,α,x) = 2^α*gamma(α/2+1)*gamma((d+α)/2+1)/gamma(d/2+1)*x[1] # α = 2*β
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit2(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
+                    @test fexplicit2(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
+                    @test fexplicit2(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+
+                    # different β, dependence on y
+                    β = 1.91239
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^(β)*y
+                    # explicit and computed solutions
+                    fexplicit3(d,α,x) = 2^α*gamma(α/2+1)*gamma((d+α)/2+1)/gamma(d/2+1)*x[2] # α = 2*β
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit3(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
+                    @test fexplicit3(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
+                    @test fexplicit3(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+                end
+                @testset "Set 4 - Explicitly known f, different non-radially-symmetric example" begin
+                    # dependence on x
+                    β = 1.21999
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^(β+1)*x
+                    # explicit and computed solutions
+                    fexplicit4(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[1] # α = 2*β
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit4(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
+                    @test fexplicit4(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
+                    @test fexplicit4(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+
+                    # different β, dependence on y
+                    β = 0.141
+                    Z = Zernike(β)
+                    WZ = Weighted(Z)
+                    xy = axes(WZ,1)
+                    x,y = first.(xy),last.(xy)
+                    # generate fractional Laplacian
+                    Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                    Δ_Zfrac = Z \ (Δfrac * WZ)
+                    # define function whose fractional Laplacian is known
+                    u = @. (1 - x^2 - y^2).^(β+1)*y
+                    # explicit and computed solutions
+                    fexplicit5(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[2] # α = 2*β
+                    f = Z*(Δ_Zfrac*(WZ \ u))
+                    # compare
+                    @test fexplicit5(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
+                    @test fexplicit5(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
+                    @test fexplicit5(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
+                end
+
+                @testset "Fractional Poisson equation on Disk: Comparison with explicitly known solutions" begin
+                    @testset "Set 1 - Radially symmetric solution" begin
+                        # define basis
+                        β = 1.1812
+                        Z = Zernike(β)
+                        WZ = Weighted(Z)
+                        xy = axes(WZ,1)
+                        x,y = first.(xy),last.(xy)
+                        # generate fractional Laplacian
+                        Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                        Δ_Zfrac = Z \ (Δfrac * WZ)
+                        # define function whose fractional Laplacian is known
+                        uexplicit = @. (1 - x^2 - y^2).^(β+1)
+                        uexplicitcfs = WZ \ uexplicit
+                        # RHS
+                        RHS(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2)/gamma(d/2)*(1-(1+α/d)*norm(x)^2) # α = 2*β
+                        RHScfs = Z \ @. RHS.(2,2*β,xy)
+                        # compute solution
+                        ucomputed = Δ_Zfrac \ RHScfs
+                        @test uexplicitcfs[1:100] ≈ ucomputed[1:100]
+                    end
+                    @testset "Set 2 - Non-radially-symmetric solutions" begin
+                        # dependence on y
+                        β = 0.98812
+                        Z = Zernike(β)
+                        WZ = Weighted(Z)
+                        xy = axes(WZ,1)
+                        x,y = first.(xy),last.(xy)
+                        # generate fractional Laplacian
+                        Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                        Δ_Zfrac = Z \ (Δfrac * WZ)
+                        # define function whose fractional Laplacian is known
+                        uexplicit = @. (1 - x^2 - y^2).^(β+1)*y
+                        uexplicitcfs = WZ \ uexplicit
+                        # RHS
+                        RHS2(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[2] # α = 2*β
+                        RHS2cfs = Z \ @. RHS2.(2,2*β,xy)
+                        # compute solution
+                        ucomputed = Δ_Zfrac \ RHS2cfs
+                        @test uexplicitcfs[1:100] ≈ ucomputed[1:100]
+                
+                        # different β, dependence on x
+                        β = 0.506
+                        Z = Zernike(β)
+                        WZ = Weighted(Z)
+                        xy = axes(WZ,1)
+                        x,y = first.(xy),last.(xy)
+                        # generate fractional Laplacian
+                        Δfrac = AbsLaplacianPower(axes(WZ,1),β)
+                        Δ_Zfrac = Z \ (Δfrac * WZ)
+                        # define function whose fractional Laplacian is known
+                        uexplicit = @. (1 - x^2 - y^2).^(β+1)*x
+                        uexplicitcfs = WZ \ uexplicit
+                        # RHS
+                        RHS3(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[1] # α = 2*β
+                        RHS3cfs = Z \ @. RHS3.(2,2*β,xy)
+                        # compute solution
+                        ucomputed = Δ_Zfrac \ RHS3cfs
+                        @test uexplicitcfs[1:100] ≈ ucomputed[1:100]
+                    end
+                end
             end
         end
     end
+
+    @testset "Complex" begin
+        @testset "∂" begin
+            A = ComplexZernike()
+            B = ComplexZernike(1)
+            xy = axes(A,1)
+            ∂ = ComplexDerivative(xy)
+            conj(∂) * A
+        end
+    end
 end
\ No newline at end of file