diff --git a/Project.toml b/Project.toml
index 294ad99..e793401 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.9.11"
+version = "0.10"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -34,15 +34,15 @@ ArrayLayouts = "1.11"
 BandedMatrices = "1"
 BlockArrays = "1.7.3"
 BlockBandedMatrices = "0.13"
-ClassicalOrthogonalPolynomials = "0.15.8"
+ClassicalOrthogonalPolynomials = "0.15.26"
 ContinuumArrays = "0.20.4, 0.21"
 DomainSets = "0.7, 0.8"
 FastTransforms = "0.17"
 FillArrays = "1.0"
-HarmonicOrthogonalPolynomials = "0.7.1, 0.8"
+HarmonicOrthogonalPolynomials = "0.8.4"
 InfiniteArrays = "0.15"
 InfiniteLinearAlgebra = "0.10"
-LazyArrays = "2.10.1"
+LazyArrays = "2.14.3"
 LazyBandedMatrices = "0.11.7"
 QuasiArrays = "0.13.3"
 Random = "1"
diff --git a/examples/cylinder.jl b/examples/cylinder.jl
index 27f636d..f02c042 100644
--- a/examples/cylinder.jl
+++ b/examples/cylinder.jl
@@ -23,7 +23,6 @@ end
 
 @test Zernike()[SVector(0.1,0.2),Block(1):(M+1)]'C*Fourier()[0.3,1:n] ≈ f(0.1,0.2,0.3)
 
-Z_ipl = inv(Z_pl)
 F_ipl = inv(F_pl)
 
 
diff --git a/examples/diskhelmholtz.jl b/examples/diskhelmholtz.jl
index ee3c3ad..af51e90 100644
--- a/examples/diskhelmholtz.jl
+++ b/examples/diskhelmholtz.jl
@@ -12,7 +12,7 @@ pyplot()
 # with the weight (1-|𝐱|^2) = (1-x^2-y^2)
 ####
 
-Z = Zernike(1)
+Z = Normalized(Zernike(1)) # orthonormal so that the mass matrix below is symmetric
 W = Weighted(Z) # w*Z
 x, y = coordinates(W)
 Δ = Z \ laplacian(W)
diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
index e589d63..411a279 100644
--- a/src/MultivariateOrthogonalPolynomials.jl
+++ b/src/MultivariateOrthogonalPolynomials.jl
@@ -21,7 +21,7 @@ import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayo
 import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout, diagtrav, invdiagtrav, ApplyBandedBlockBandedLayout, krontrav
 import InfiniteArrays: InfiniteCardinal, OneToInf
 
-import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw, OPLayout, normalized
+import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw, OPLayout, normalized, normalizationconstant
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
                                           AngularMomentum, angularmomentum, BlockOneTo, BlockRange1, interlace,
                                           MultivariateOPLayout, AbstractMultivariateOPLayout, MAX_PLOT_BLOCKS
@@ -30,7 +30,7 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        UnitTriangle, UnitDisk,
        JacobiTriangle, TriangleWeight, WeightedTriangle,
        DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
-       Zernike, ZernikeWeight, zerniker, zernikez,
+       Zernike, ZernikeWeight, zerniker, zernikez, normalizedzerniker, normalizedzernikez,
        AngularMomentum,
        RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial,
        grammatrix, oneto, coordinates, Laplacian, AbsLaplacian, laplacian, abslaplacian, angularmomentum, weaklaplacian
diff --git a/src/disk.jl b/src/disk.jl
index f3d2a75..c7536f7 100644
--- a/src/disk.jl
+++ b/src/disk.jl
@@ -38,7 +38,12 @@ end
 """
     Zernike(a, b)
 
-is a quasi-matrix orthogonal `r^(2a) * (1-r^2)^b`
+is a quasi-matrix orthogonal `r^(2a) * (1-r^2)^b`. The polynomials are not normalized:
+the entry of degree `ℓ` and Fourier mode `m` is
+
+    r^|m| * jacobip((ℓ-|m|) ÷ 2, b, |m|+a, 2r^2-1) * (signbit(m) ? sin(|m|*θ) : cos(|m|*θ))
+
+Use `Normalized(Zernike(a, b))` for the orthonormal polynomials.
 """
 struct Zernike{T} <: BivariateOrthogonalPolynomial{T}
     a::T
@@ -74,37 +79,109 @@ orthogonalityweight(Z::Zernike) = ZernikeWeight(Z.a, Z.b)
 
 basis_axes(::Inclusion{<:Any,<:UnitDisk}, v) = Zernike()
 
-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)
+###
+# Evaluation
+###
+
+zerniker(ℓ, m, a, b, r::T) where T = r^m * jacobip((ℓ-m) ÷ 2, b, m+a, 2r^2-1)
 zerniker(ℓ, m, b, r) = zerniker(ℓ, m, zero(b), b, r)
 zerniker(ℓ, m, r) = zerniker(ℓ, m, zero(r), r)
 
-function zernikez(ℓ, ms, a, b, rθ::RadialCoordinate{T}) where T
+normalizedzerniker(ℓ, 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)
+normalizedzerniker(ℓ, m, b, r) = normalizedzerniker(ℓ, m, zero(b), b, r)
+normalizedzerniker(ℓ, m, r) = normalizedzerniker(ℓ, m, zero(r), r)
+
+function _zernikez(radial, ℓ, ms, a, b, rθ::RadialCoordinate)
     r,θ = rθ.r,rθ.θ
     m = abs(ms)
-    zerniker(ℓ, m, a, b, r) * (signbit(ms) ? sin(m*θ) : cos(m*θ))
+    radial(ℓ, 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)
+for (z, radial) in ((:zernikez, :zerniker), (:normalizedzernikez, :normalizedzerniker))
+    @eval begin
+        $z(ℓ, ms, a, b, rθ::RadialCoordinate) = _zernikez($radial, ℓ, ms, a, b, rθ)
+        $z(ℓ, ms, a, b, xy::StaticVector{2}) = $z(ℓ, ms, a, b, RadialCoordinate(xy))
+        $z(ℓ, ms, b, xy::StaticVector{2}) = $z(ℓ, ms, zero(b), b, xy)
+        $z(ℓ, ms, xy::StaticVector{2,T}) where T = $z(ℓ, ms, zero(T), xy)
+    end
+end
 
-function getindex(Z::Zernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
+# degree ℓ and signed Fourier mode m of the entry B
+function _zernikelm(B::BlockIndex{1})
     ℓ = Int(block(B))-1
     k = blockindex(B)
     m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
-    zernikez(ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
+    ℓ, (isodd(k+ℓ) ? 1 : -1) * m
 end
 
+const NormalizedZernike{T} = Normalized{T,Zernike{T}}
+
+getindex(Z::Zernike, rθ::RadialCoordinate, B::BlockIndex{1}) = zernikez(_zernikelm(B)..., Z.a, Z.b, rθ)
+getindex(Q::NormalizedZernike, rθ::RadialCoordinate, B::BlockIndex{1}) = normalizedzernikez(_zernikelm(B)..., Q.P.a, Q.P.b, rθ)
+getindex(Q::NormalizedZernike, xy::StaticVector{2}, B::BlockIndex{1}) = Q[RadialCoordinate(xy), B]
 
 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])])
 
+###
+# Normalization constants
+###
+
+"""
+    zernikenormalizationconstant(T, n, m, a, b)
+
+gives `c` such that `normalizedzerniker(ℓ, m, a, b, r) == c * zerniker(ℓ, m, a, b, r)` where `n == (ℓ-m) ÷ 2`.
+"""
+function zernikenormalizationconstant(::Type{T}, n::Integer, m::Integer, a, b) where T
+    α, β = convert(T, b), convert(T, m+a)
+    # ratio of 2^(α+β+1) to the norm squared of the Jacobi polynomial
+    lr = iszero(n) ? loggamma(α+β+2) - loggamma(α+1) - loggamma(β+1) :
+                     log(2n+α+β+1) + loggamma(n+α+β+1) + loggamma(n+one(T)) - loggamma(n+α+1) - loggamma(n+β+1)
+    sqrt(convert(T,2)^(1-iszero(m)) / convert(T,π) * exp(lr))
+end
+
+# the scaling such that Normalized(Z) == Z * Diagonal(ZernikeNormalizationConstant(Z.a, Z.b))
+struct ZernikeNormalizationConstant{T} <: AbstractBlockVector{T}
+    a::T
+    b::T
+end
+
+ZernikeNormalizationConstant{T}(Z::Zernike) where T = ZernikeNormalizationConstant{T}(Z.a, Z.b)
+
+axes(::ZernikeNormalizationConstant) = (blockedrange(oneto(∞)),)
+copy(c::ZernikeNormalizationConstant) = c
+
+MemoryLayout(::Type{<:ZernikeNormalizationConstant}) = LazyLayout()
+Base.BroadcastStyle(::Type{<:ZernikeNormalizationConstant}) = LazyArrayStyle{1}()
+Base.BroadcastStyle(::Type{<:Diagonal{<:Any,<:ZernikeNormalizationConstant}}) = LazyArrayStyle{2}()
+
+function getindex(c::ZernikeNormalizationConstant{T}, Kk::BlockIndex{1}) where T
+    ℓ, ms = _zernikelm(Kk)
+    m = abs(ms)
+    zernikenormalizationconstant(T, (ℓ-m) ÷ 2, m, c.a, c.b)
+end
+getindex(c::ZernikeNormalizationConstant, k::Integer) = c[findblockindex(axes(c,1),k)]
+Base.view(c::ZernikeNormalizationConstant, K::Block{1}) = [c[K[j]] for j = 1:Int(K)]
+
+normalizationconstant(Z::Zernike{T}) where T = ZernikeNormalizationConstant{real(T)}(Z)
+
+# the normalization constants laid out as the coefficient matrix of a ModalTrav:
+# row i corresponds to n = i-1 and column j to Fourier mode j ÷ 2
+_zernikemodalconstant(::Type{T}, (M,N)::NTuple{2,Int}, a, b) where T =
+    [zernikenormalizationconstant(T, i-1, j ÷ 2, a, b) for i = 1:M, j = 1:N]
+
+# coefficients in Normalized(Z) to coefficients in Z
+_normalized2zernike(Z::Zernike, c::ModalTrav) = ModalTrav(_zernikemodalconstant(real(eltype(Z)), size(c.matrix), Z.a, Z.b) .* c.matrix)
+# coefficients in Z to coefficients in Normalized(Z)
+_zernike2normalized(Z::Zernike, c::AbstractVector) = c ./ normalizationconstant(Z)[axes(c,1)]
+
 
 ###
 # Jacobi matrices
 ###
-function jacobimatrix(::Val{1}, Z::Zernike{T}) where T
+function jacobimatrix(::Val{1}, Q::NormalizedZernike{T}) where T
+    Z = Q.P
     if iszero(Z.a)
         α = Z.b     # extract second basis parameter
 
@@ -153,7 +230,8 @@ function jacobimatrix(::Val{1}, Z::Zernike{T}) where T
     end
 end
 
-function jacobimatrix(::Val{2}, Z::Zernike{T}) where T
+function jacobimatrix(::Val{2}, Q::NormalizedZernike{T}) where T
+    Z = Q.P
     if iszero(Z.a)
         α = Z.b     # extract second basis parameter
 
@@ -194,6 +272,15 @@ function jacobimatrix(::Val{2}, Z::Zernike{T}) where T
     end
 end
 
+# Normalized(Z) == Z * D so x * Z == Z * D * X * inv(D) where X is the Jacobi matrix of Normalized(Z)
+function jacobimatrix(v::Union{Val{1},Val{2}}, Z::Zernike)
+    s = Normalized(Z).scaling
+    _bandedblockbandeddiagonal(s) * jacobimatrix(v, Normalized(Z)) * _bandedblockbandeddiagonal(inv.(s))
+end
+
+# a Diagonal would lose the block-banded structure in the product
+_bandedblockbandeddiagonal(d::AbstractVector) = _BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d)', axes(d,1), (0,0), (0,0))
+
 ###
 # Transforms
 ###
@@ -226,7 +313,7 @@ function plotvalues(u::ApplyQuasiVector{T,typeof(*),<:Tuple{Zernike, AbstractVec
     B = findblock(axes(Z,2), last(colsupport(c)))
     N = Int(B) ÷ 2 + 1 # polynomial degree
     F = ZernikeITransform{T}(min(2N, MAX_PLOT_BLOCKS), Z.a, Z.b)
-    C = F * c[Block.(OneTo(min(2N, MAX_PLOT_BLOCKS)))] # transform to grid
+    C = F * _zernike2normalized(Z, c[Block.(OneTo(min(2N, MAX_PLOT_BLOCKS)))]) # transform to grid
     [permutedims(u[x[1,:]]); # evaluate on edge of disk
      C C[:,1];
      fill(u[x[end,1]], 1, size(x,2))] # evaluate at origin and repeat
@@ -267,7 +354,9 @@ deduceeltype(T, f) = all(isreal, f) ? T : Complex{T}
 
 inv(P::ZernikeTransform) = ZernikeITransform(P.N, P.disk2cxf, inv(P.analysis))
 
-plan_transform(Z::Zernike{T}, (N,)::Tuple{Block{1}}, dims=1) where T = ZernikeTransform{real(T)}(Int(N), Z.a, Z.b)
+plan_transform(Q::NormalizedZernike{T}, (N,)::Tuple{Block{1}}, dims=1) where T = ZernikeTransform{real(T)}(Int(N), Q.P.a, Q.P.b)
+# transforms for Zernike(a,b) rescale those for Normalized(Zernike(a,b))
+plan_transform(Z::Zernike{T}, (N,)::Tuple{Block{1}}, dims=1) where T = ApplyPlan(Base.Fix1(_normalized2zernike, Z), ZernikeTransform{real(T)}(Int(N), Z.a, Z.b))
 
 ##
 # Laplacian
@@ -281,10 +370,10 @@ end
 function laplacian(Z::Zernike{T}; dims...) where T
     a,b = Z.a,Z.b
     @assert a == 0
-    D = Derivative(Inclusion(ChebyshevInterval{T}())) 
-    Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 4(HalfWeighted{:b}(C)\(D*HalfWeighted{:b}(B)))*(B\(D*A)), Normalized.(Jacobi.(b+2,a:∞)), Normalized.(Jacobi.(b+1,(a+1):∞)), Normalized.(Jacobi.(b,a:∞)))
+    D = Derivative(Inclusion(ChebyshevInterval{T}()))
+    Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 8(HalfWeighted{:b}(C)\(D*HalfWeighted{:b}(B)))*(B\(D*A)), Jacobi{T}.(b+2,a:∞), Jacobi{T}.(b+1,(a+1):∞), Jacobi{T}.(b,a:∞))
     Δ = ModalInterlace(Δs, (ℵ₀,ℵ₀), (-2,2))
-    Zernike(a,b+2) * Δ
+    Zernike{T}(a,b+2) * Δ
 end
 
 ###
@@ -329,18 +418,18 @@ 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}((axes(A,2),))
+    A == B && return Eye{TV}((axes(A,2),))
     st = Int(A.a - B.a + A.b - B.b)
-    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(A.b,A.a:∞)) .\ Normalized.(Jacobi{TV}.(B.b,B.a:∞))) .* convert(TV, 2)^(-st/2), (ℵ₀,ℵ₀), (0,2st))
+    ModalInterlace{TV}(Jacobi{TV}.(A.b,A.a:∞) .\ Jacobi{TV}.(B.b,B.a:∞), (ℵ₀,ℵ₀), (0,2st))
 end
 
 function \(A::Zernike{T}, wB::Weighted{V,Zernike{V}}) where {T,V}
     TV = promote_type(T,V)
     B = wB.P
     A.a == B.a == A.b == B.b == 0 && return Eye{TV}((axes(A,2),))
-    c = Int(B.a - A.a + B.b - A.b)
     @assert iszero(B.a)
-    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(A.b, A.a:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(B.b, B.a:∞)))) .* convert(TV, 2)^(-c/2), (ℵ₀,ℵ₀), (2Int(B.b), 2Int(A.a+A.b)))
+    # (1-r^2)^b == 2^(-b) * (1-s)^b where s = 2r^2-1
+    ModalInterlace{TV}((Jacobi{TV}.(A.b, A.a:∞) .\ HalfWeighted{:a}.(Jacobi{TV}.(B.b, B.a:∞))) ./ convert(TV, 2)^B.b, (ℵ₀,ℵ₀), (2Int(B.b), 2Int(A.a+A.b)))
 end
 
 
@@ -351,5 +440,5 @@ end
 function Base._sum(P::Zernike{T}, dims) where T
     @assert dims == 1
     @assert P.a == P.b == 0
-    Hcat(sqrt(convert(T, π)), Zeros{T}(1,∞))
-end
\ No newline at end of file
+    Hcat(convert(T, π), Zeros{T}(1,∞))
+end
diff --git a/test/test_disk.jl b/test/test_disk.jl
index a94a3e1..1b57d7c 100644
--- a/test/test_disk.jl
+++ b/test/test_disk.jl
@@ -42,13 +42,40 @@ import ForwardDiff: hessian
         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θ,1] ≈ Zernike()[xy,1] ≈ 1 ≈ zernikez(0, 0, rθ)
+        @test Zernike()[rθ,Block(1)] ≈ Zernike()[xy,Block(1)] ≈ [1]
+        @test Zernike()[rθ,Block(2)] ≈ [r*sin(θ), r*cos(θ)] ≈ [zernikez(1, -1, rθ), zernikez(1, 1, rθ)]
+        @test Zernike()[rθ,Block(3)] ≈ [2r^2-1, r^2*sin(2θ), 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)
+        @test zerniker(5, 1, 0.3, 0.2, r) ≈ r * jacobip(2, 0.2, 1.3, 2r^2-1)
+
+        Q = Normalized(Zernike())
+        @test Q[rθ,1] ≈ Q[xy,1] ≈ inv(sqrt(π)) ≈ normalizedzernikez(0, 0, rθ)
+        @test Q[rθ,Block(1)] ≈ Q[xy,Block(1)] ≈ [inv(sqrt(π))]
+        @test Q[rθ,Block(2)] ≈ [2r/sqrt(π)*sin(θ), 2r/sqrt(π)*cos(θ)] ≈ [normalizedzernikez(1, -1, rθ), normalizedzernikez(1, 1, rθ)]
+        @test Q[rθ,Block(3)] ≈ [sqrt(3/π)*(2r^2-1),sqrt(6/π)*r^2*sin(2θ),sqrt(6/π)*r^2*cos(2θ)] ≈ [normalizedzernikez(2, 0, rθ), normalizedzernikez(2, -2, rθ), normalizedzernikez(2, 2, rθ)]
+        @test Q[rθ,Block(4)] ≈ [normalizedzernikez(3, -1, rθ), normalizedzernikez(3, 1, rθ), normalizedzernikez(3, -3, rθ), normalizedzernikez(3, 3, rθ)]
+        @test Q[xy,1:6] ≈ Q[xy,Block.(1:3)]
+        @test Q[xy,Block.(2:3)] ≈ Q[xy,2:6]
+
+        @test normalizedzerniker(5, 0, norm(xy)) ≈ normalizedzernikez(5, 0, xy)
+    end
+
+    @testset "Normalized" begin
+        xy = SVector(0.1,0.2)
+        for (a,b) in ((0,0), (0.1,0.2), (0,1))
+            Z = Zernike(a,b)
+            Q = Normalized(Z)
+            @test Q[xy,Block.(1:10)] ≈ Z[xy,Block.(1:10)] .* Q.scaling[1:55]
+        end
+        @test Normalized(Zernike()) == Normalized(Zernike())
+        @test Normalized(Zernike()) ≠ Normalized(Zernike(1))
+        @test Normalized(Zernike()) ≠ Zernike()
+        @test Zernike() ≠ Normalized(Zernike())
+        @test copy(Normalized(Zernike())) ≡ Normalized(Zernike())
+        @test stringmime("text/plain", Normalized(Zernike(1))) == "Normalized(Zernike(0.0, 1.0))"
     end
 
     @testset "ModalTrav" begin
@@ -79,70 +106,84 @@ import ForwardDiff: hessian
     @testset "expand" begin
         @test expand(Zernike(), splat((x,y) -> exp(x*cos(y))))[SVector(0.1,0.2)] ≈ expand(Zernike{ComplexF64}(), splat((x,y) -> exp(x*cos(y))))[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2))
         @test expand(Zernike{ComplexF64}(), splat((x,y) -> exp(x*cos(y)+im*y)))[SVector(0.1,0.2)] ≈ expand(Zernike(), splat((x,y) -> exp(x*cos(y)+im*y)))[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2)+im*0.2)
+        @test expand(Normalized(Zernike()), splat((x,y) -> exp(x*cos(y))))[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2))
+        @test expand(Normalized(Zernike()), splat((x,y) -> exp(x*cos(y)+im*y)))[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2)+im*0.2)
     end
 
     @testset "Jacobi matrices" begin
         # Setup
         α = 10 * rand()
-        Z = Zernike(α)
-        x, y = coordinates(Z)
-        n = 150
-
-        # X tests
-        JX = zeros(n,n)
-        for j = 1:n
-            JX[1:n,j] = (Z \ (x .* Z[:,j]))[1:n]
-        end 
-        X = Z \ (x .* Z)
-        # The Zernike Jacobi matrices are symmetric for this normalization
-        @test issymmetric(X)
-        # Consistency with expansion
-        @test X[1:150,1:150] ≈ JX
-        # Multiplication by x
-        f = Z \ (sin.(x.*y) .+ x.^2 .- y)
-        xf = Z \ (x.*sin.(x.*y) .+ x.^3 .- x.*y)
-        @test X[Block.(1:20),Block.(1:20)]*f[Block.(1:20)] ≈ xf[Block.(1:20)]
-
-        # Y tests
-        JY = zeros(n,n)
-        for j = 1:n
-            JY[1:n,j] = (Z \ (y .* Z[:,j]))[1:n]
-        end 
-        Y = Z \ (y .* Z)
-        # The Zernike Jacobi matrices are symmetric for this normalization
-        @test issymmetric(Y)
-        # Consistency with expansion
-        @test Y[1:150,1:150] ≈ JY
-        # Multiplication by y
-        f = Z \ (sin.(x.*y) .+ x.^2 .- y)
-        yf = Z \ (y.*sin.(x.*y) .+ y .* x.^2 .- y.^2)
-        @test Y[Block.(1:20),Block.(1:20)]*f[Block.(1:20)] ≈ yf[Block.(1:20)]
+        for Z in (Normalized(Zernike(α)), Zernike(α))
+            x, y = coordinates(Z)
+            n = 150
+
+            # X tests
+            JX = zeros(n,n)
+            for j = 1:n
+                JX[1:n,j] = (Z \ (x .* Z[:,j]))[1:n]
+            end 
+            X = Z \ (x .* Z)
+            # The Zernike Jacobi matrices are symmetric for the orthonormal polynomials
+            @test issymmetric(X) == (Z isa Normalized)
+            # Consistency with expansion
+            @test X[1:150,1:150] ≈ JX
+            # Multiplication by x
+            f = Z \ (sin.(x.*y) .+ x.^2 .- y)
+            xf = Z \ (x.*sin.(x.*y) .+ x.^3 .- x.*y)
+            @test X[Block.(1:20),Block.(1:20)]*f[Block.(1:20)] ≈ xf[Block.(1:20)]
+
+            # Y tests
+            JY = zeros(n,n)
+            for j = 1:n
+                JY[1:n,j] = (Z \ (y .* Z[:,j]))[1:n]
+            end 
+            Y = Z \ (y .* Z)
+            # The Zernike Jacobi matrices are symmetric for the orthonormal polynomials
+            @test issymmetric(Y) == (Z isa Normalized)
+            # Consistency with expansion
+            @test Y[1:150,1:150] ≈ JY
+            # Multiplication by y
+            f = Z \ (sin.(x.*y) .+ x.^2 .- y)
+            yf = Z \ (y.*sin.(x.*y) .+ y .* x.^2 .- y.^2)
+            @test Y[Block.(1:20),Block.(1:20)]*f[Block.(1:20)] ≈ yf[Block.(1:20)]
             
-        # Multiplication of Jacobi matrices
-        @test (X*X)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:10),Block.(1:10)]*X[Block.(1:10),Block.(1:10)])[Block.(1:6),Block.(1:6)]
-        @test (X*Y)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:10),Block.(1:10)]*Y[Block.(1:10),Block.(1:10)])[Block.(1:6),Block.(1:6)]
+            # Multiplication of Jacobi matrices
+            @test (X*X)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:10),Block.(1:10)]*X[Block.(1:10),Block.(1:10)])[Block.(1:6),Block.(1:6)]
+            @test (X*Y)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:10),Block.(1:10)]*Y[Block.(1:10),Block.(1:10)])[Block.(1:6),Block.(1:6)]
+
+            # Addition of Jacobi matrices
+            @test (X+Y)[Block.(1:6),Block.(1:6)] ≈ X[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]
+            @test (Y+Y)[Block.(1:6),Block.(1:6)] ≈ Y[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]
+        end
 
-        # Addition of Jacobi matrices
-        @test (X+Y)[Block.(1:6),Block.(1:6)] ≈ X[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]
-        @test (Y+Y)[Block.(1:6),Block.(1:6)] ≈ Y[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]
-            
         # for now, reject non-zero first parameter options
         @test_throws ErrorException("Implement for non-zero first basis parameter.") jacobimatrix(Val(1),Zernike(1,1))  
         @test_throws ErrorException("Implement for non-zero first basis parameter.") jacobimatrix(Val(2),Zernike(1,1))
+        @test_throws ErrorException("Implement for non-zero first basis parameter.") jacobimatrix(Val(1),Normalized(Zernike(1,1)))
     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 (a,b) in ((0,0), (0.1, 0.2), (0,1)), Z in (Zernike(a,b), Normalized(Zernike(a,b)))
+            Zn = Z[:,Block.(Base.OneTo(3))]
             for k = 1:6
-                @test factorize(Zn) \ Zernike(a,b)[:,k] ≈ [zeros(k-1); 1; zeros(6-k)]
+                @test factorize(Zn) \ Z[:,k] ≈ [zeros(k-1); 1; zeros(6-k)]
             end
 
-            Z = Zernike(a,b);
             x,y = coordinates(Z)
             u = Z * (Z \ exp.(x .* cos.(y)))
             @test u[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2))
         end
+
+        # coefficients differ by the normalization constants
+        Z = Zernike(0.1, 0.2)
+        x,y = coordinates(Z)
+        f = exp.(x .* cos.(y))
+        @test (Z \ f)[1:55] ≈ Normalized(Z).scaling[1:55] .* (Normalized(Z) \ f)[1:55]
+
+        P = plan_transform(Z, Block(5))
+        c = BlockedArray(randn(sum(1:5)), 1:5)
+        V = [(Z * [c; zeros(∞)])[SVector(𝐱)] for 𝐱 in grid(Z, Block(5))]
+        @test P * V ≈ c
     end
 
     @testset "Laplacian" begin
@@ -166,8 +207,8 @@ import ForwardDiff: hessian
 
         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)
+        @test ur(r) ≈ normalizedzerniker(ℓ, m, b, r)
+        @test f(r^2) ≈ r^(-m) * normalizedzerniker(ℓ, 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)
 
@@ -181,7 +222,7 @@ import ForwardDiff: hessian
 
         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)
+        @test ur(r) ≈ (1-r^2) * normalizedzerniker(ℓ, m, b, r)
 
         D = Derivative(Chebyshev())
         D1 = Normalized(Jacobi(0, m+1)) \ (D * (HalfWeighted{:a}(Normalized(Jacobi(1, m)))))
@@ -229,114 +270,160 @@ import ForwardDiff: hessian
         # 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(WZ)
-        Δ_Z = Zernike(1) \ (Δ * WZ)
-        @test exp.(Δ_Z)[1:10,1:10] == exp.(Δ_Z[1:10,1:10])
+        for Z in (Zernike(1), Normalized(Zernike(1)))
+            WZ = Weighted(Z) # Z weighted by (1-r^2)
+            Δ = Laplacian(WZ)
+            Δ_Z = Z \ (Δ * WZ)
+            @test exp.(Δ_Z)[1:10,1:10] == exp.(Δ_Z[1:10,1:10])
 
-        x,y = coordinates(WZ)
-        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]
+            x,y = coordinates(WZ)
+            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]  ≈ (Z \ Δu)[1:100]
+        end
 
         @testset "Unweighted" begin
             c = [randn(100); zeros(∞)]
-            Z = Zernike()
-            Δ = Zernike(2) \ (Laplacian(Z) * Z)
-            @test tr(hessian(xy -> (Zernike{eltype(xy)}()*c)[xy], SVector(0.1,0.2))) ≈ (Zernike(2)*(Δ*c))[SVector(0.1,0.2)]
-
-            b = 0.2
-            Z = Zernike(b)
-            Δ = Zernike(b+2) \ (Laplacian(Z) * Z)
-            @test tr(hessian(xy -> (Zernike{eltype(xy)}(b)*c)[xy], SVector(0.1,0.2))) ≈ (Zernike(b+2)*(Δ*c))[SVector(0.1,0.2)]
+            for b in (0, 0.2), normalize in (identity, Normalized)
+                Z = normalize(Zernike(b))
+                Z₂ = normalize(Zernike(b+2))
+                Δ = Z₂ \ (Laplacian(Z) * Z)
+                @test tr(hessian(xy -> (normalize(Zernike{eltype(xy)}(b))*c)[xy], SVector(0.1,0.2))) ≈ (Z₂*(Δ*c))[SVector(0.1,0.2)]
+            end
         end
     end
 
     @testset "Conversion" begin
+        xy = SVector(0.1,0.2)
+
+        # orthonormal
         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))
 
-        xy = SVector(0.1,0.2)
-        @test Zernike()[xy,Block(1)[1]] ≈ Zernike(1)[xy,Block(1)[1]]/sqrt(2)
+        Q, Q₁ = Normalized(Zernike()), Normalized(Zernike(1))
+        @test Q[xy,Block(1)[1]] ≈ Q₁[xy,Block(1)[1]]/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 Q[xy,Block(2)[1]] ≈ Q₁[xy,Block(2)[1]]*R1[1,1]/sqrt(2)
+        @test Q[xy,Block(2)[2]] ≈ Q₁[xy,Block(2)[2]]*R1[1,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 Q[xy,Block(3)[1]] ≈ R0[1:2,2]'*Q₁[xy,getindex.(Block.(1:2:3),1)]/sqrt(2)
+        @test Q[xy,Block(3)[2]] ≈ R2[1,1]*Q₁[xy,Block(3)[2]]/sqrt(2)
+        @test Q[xy,Block(3)[3]] ≈ R2[1,1]*Q₁[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)
+        @test Q[xy,Block(4)[1]] ≈ R1[1:2,2]'*Q₁[xy,getindex.(Block.(2:2:4),1)]/sqrt(2)
+        @test Q[xy,Block(4)[2]] ≈ R1[1:2,2]'*Q₁[xy,getindex.(Block.(2:2:4),2)]/sqrt(2)
+        @test Q[xy,Block(4)[3]] ≈ R3[1,1]*Q₁[xy,Block(4)[3]]/sqrt(2)
+        @test Q[xy,Block(4)[4]] ≈ R3[1,1]*Q₁[xy,Block(4)[4]]/sqrt(2)
 
-        @test Zernike()[xy,Block(5)[1]] ≈ R0[2:3,3]'*Zernike(1)[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
+        @test Q[xy,Block(5)[1]] ≈ R0[2:3,3]'*Q₁[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
 
+        # unnormalized
+        R0 = Jacobi(1, 0) \ Jacobi(0, 0)
+        R1 = Jacobi(1, 1) \ Jacobi(0, 1)
+        R2 = Jacobi(1, 2) \ Jacobi(0, 2)
 
-        R = Zernike(1) \ Zernike()
+        @test Zernike()[xy,Block(1)[1]] ≈ Zernike(1)[xy,Block(1)[1]]
+        @test Zernike()[xy,Block(2)[1]] ≈ Zernike(1)[xy,Block(2)[1]]*R1[1,1]
+        @test Zernike()[xy,Block(3)[1]] ≈ R0[1:2,2]'*Zernike(1)[xy,getindex.(Block.(1:2:3),1)]
+        @test Zernike()[xy,Block(3)[2]] ≈ R2[1,1]*Zernike(1)[xy,Block(3)[2]]
+        @test Zernike()[xy,Block(4)[1]] ≈ R1[1:2,2]'*Zernike(1)[xy,getindex.(Block.(2:2:4),1)]
+        @test Zernike()[xy,Block(5)[1]] ≈ R0[2:3,3]'*Zernike(1)[xy,getindex.(Block.(3:2:5),1)]
 
+        R = Zernike(1) \ Zernike()
         @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)]
 
-        R = Zernike(2) \ Zernike()
-        @test Zernike()[xy,Block.(1:6)]' ≈ Zernike(2)[xy,Block.(1:6)]'*R[Block.(1:6),Block.(1:6)]
+        for (A,B) in ((Zernike(1), Zernike()), (Zernike(2), Zernike()), (Zernike(0.1,1.2), Zernike(0.1,0.2)),
+                      (Normalized(Zernike(1)), Normalized(Zernike())), (Normalized(Zernike(2)), Normalized(Zernike())),
+                      (Zernike(), Normalized(Zernike())), (Normalized(Zernike()), Zernike()),
+                      (Zernike(1), Normalized(Zernike())), (Normalized(Zernike(2)), Zernike()))
+            R = A \ B
+            @test B[xy,Block.(1:6)]' ≈ A[xy,Block.(1:6)]'*R[Block.(1:6),Block.(1:6)]
+        end
+        @test Zernike() \ Zernike() isa Eye
+        @test Normalized(Zernike()) \ Normalized(Zernike()) isa Eye
+        @test (Zernike() \ Normalized(Zernike()))[1:10,1:10] ≈ Diagonal(Normalized(Zernike()).scaling[1:10])
     end
 
     @testset "Lowering" begin
+        xy = SVector(0.1,0.2)
+        r = norm(xy)
+        w = 1 - r^2
+
+        # orthonormal
         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)))
 
-        xy = SVector(0.1,0.2)
-        r = norm(xy)
-        w = 1 - r^2
+        Q, Q₁ = Normalized(Zernike()), Normalized(Zernike(1))
+        @test w*Q₁[xy,Block(1)[1]] ≈ L0[1:2,1]'*Q[xy,getindex.(Block.(1:2:3),1)] / sqrt(2)
 
-        @test w*Zernike(1)[xy,Block(1)[1]] ≈ L0[1:2,1]'*Zernike()[xy,getindex.(Block.(1:2:3),1)] / sqrt(2)
+        @test w*Q₁[xy,Block(2)[1]] ≈ L1[1:2,1]'*Q[xy,getindex.(Block.(2:2:4),1)]/sqrt(2)
+        @test w*Q₁[xy,Block(2)[2]] ≈ L1[1:2,1]'*Q[xy,getindex.(Block.(2:2:4),2)]/sqrt(2)
 
-        @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 w*Q₁[xy,Block(3)[1]] ≈ L0[2:3,2]'*Q[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
+        @test w*Q₁[xy,Block(3)[2]] ≈ L2[1:2,1]'*Q[xy,getindex.(Block.(3:2:5),2)]/sqrt(2)
+        @test w*Q₁[xy,Block(3)[3]] ≈ L2[1:2,1]'*Q[xy,getindex.(Block.(3:2:5),3)]/sqrt(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)
+        # unnormalized: (1-r^2) == (1-s)/2 where s = 2r^2-1
+        L0 = Jacobi(0, 0) \ HalfWeighted{:a}(Jacobi(1, 0))
+        L1 = Jacobi(0, 1) \ HalfWeighted{:a}(Jacobi(1, 1))
+        L2 = Jacobi(0, 2) \ HalfWeighted{:a}(Jacobi(1, 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(1)[1]] ≈ L0[1:2,1]'*Zernike()[xy,getindex.(Block.(1:2:3),1)]/2
+        @test w*Zernike(1)[xy,Block(2)[1]] ≈ L1[1:2,1]'*Zernike()[xy,getindex.(Block.(2:2:4),1)]/2
+        @test w*Zernike(1)[xy,Block(3)[1]] ≈ L0[2:3,2]'*Zernike()[xy,getindex.(Block.(3:2:5),1)]/2
+        @test w*Zernike(1)[xy,Block(3)[2]] ≈ L2[1:2,1]'*Zernike()[xy,getindex.(Block.(3:2:5),2)]/2
 
+        L = Zernike() \ Weighted(Zernike(1))
         @test exp.(L)[1:10,1:10] == exp.(L[1:10,1:10])
 
-        L = Zernike(1) \ Weighted(Zernike(1))
-        @test 2w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike(1)[xy,Block.(1:7)]'*L[Block.(1:7),Block.(1:5)]
+        for normalize in (identity, Normalized)
+            A, B = normalize(Zernike()), normalize(Zernike(1))
+            L = A \ Weighted(B)
+            @test w*B[xy,Block.(1:5)] ≈ transpose(L[Block.(1:7),Block.(1:5)])*A[xy,Block.(1:7)]
 
-        L = Zernike() \ Weighted(Zernike(2))
-        @test w^2*Zernike(2)[xy,Block.(1:5)]' ≈ Zernike()[xy,Block.(1:9)]'*L[Block.(1:9),Block.(1:5)]
+            L = B \ Weighted(B)
+            @test w*B[xy,Block.(1:5)] ≈ transpose(L[Block.(1:7),Block.(1:5)])*B[xy,Block.(1:7)]
+
+            B = normalize(Zernike(2))
+            L = A \ Weighted(B)
+            @test w^2*B[xy,Block.(1:5)] ≈ transpose(L[Block.(1:9),Block.(1:5)])*A[xy,Block.(1:9)]
+        end
+
+        L = Zernike() \ Weighted(Normalized(Zernike(1)))
+        @test w*Q₁[xy,Block.(1:5)] ≈ transpose(L[Block.(1:7),Block.(1:5)])*Zernike()[xy,Block.(1:7)]
+        @test Zernike() \ Weighted(Zernike()) isa Eye
     end
 
     @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 all(rep[1].args .≈ (first.(g),last.(g),u[g]))
-
-        W = Weighted(Zernike(1))
-        u = W * [1; 2; zeros(∞)];
-        rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
-        g = MultivariateOrthogonalPolynomials.plotgrid(W[:,1:3])
-        @test all(rep[1].args .≈ (first.(g),last.(g),u[g]))
+        for Z in (Zernike(), Normalized(Zernike()))
+            u = Z * [1; 2; zeros(∞)];
+            rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u);
+            g = MultivariateOrthogonalPolynomials.plotgrid(Z[:,1:3])
+            @test all(rep[1].args .≈ (first.(g),last.(g),u[g]))
+        end
+
+        for W in (Weighted(Zernike(1)), Weighted(Normalized(Zernike(1))))
+            u = W * [1; 2; zeros(∞)];
+            rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
+            g = MultivariateOrthogonalPolynomials.plotgrid(W[:,1:3])
+            @test all(rep[1].args .≈ (first.(g),last.(g),u[g]))
+        end
     end
 
     @testset "sum" begin
-        P = Zernike()
-        x,y = coordinates(P)
-        @test sum(expand(P, 𝐱 -> 1)) ≈ π
-        @test sum(expand(P, 𝐱 -> let (x,y) = 𝐱; exp(x*cos(y)) end)) ≈ sum(exp.(x.*cos.(y))) ≈ sum(exp.(x.*cos.(y)) for (x,y) in UnitDisk()) ≈ 3.4898933353782744
-        @test [sum(P[:,k] .* P[:,j]) for k=1:10, j=1:10] ≈ I
+        for P in (Zernike(), Normalized(Zernike()))
+            x,y = coordinates(P)
+            @test sum(expand(P, 𝐱 -> 1)) ≈ π
+            @test sum(expand(P, 𝐱 -> let (x,y) = 𝐱; exp(x*cos(y)) end)) ≈ sum(exp.(x.*cos.(y))) ≈ sum(exp.(x.*cos.(y)) for (x,y) in UnitDisk()) ≈ 3.4898933353782744
+        end
+        Q = Normalized(Zernike())
+        @test [sum(Q[:,k] .* Q[:,j]) for k=1:10, j=1:10] ≈ I
+        Z = Zernike()
+        @test [sum(Z[:,k] .* Z[:,j]) for k=1:10, j=1:10] ≈ Diagonal(inv.(Q.scaling[1:10]).^2)
     end
 
     @testset "Show" begin