diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
index 5d1f42c..b8ffc22 100644
--- a/src/MultivariateOrthogonalPolynomials.jl
+++ b/src/MultivariateOrthogonalPolynomials.jl
@@ -4,7 +4,7 @@ using QuasiArrays: AbstractVector
 using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets,
       QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra,
       LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
-      HarmonicOrthogonalPolynomials, RecurrenceRelationships
+      HarmonicOrthogonalPolynomials, RecurrenceRelationships, FillArrays
 
 import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary, diff
 import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
@@ -14,6 +14,7 @@ import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain
 import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis
 
 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
+import FillArrays: SquareEye
 import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle
 import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, _BandedMatrix, blockbandwidths, subblockbandwidths
 import LinearAlgebra: factorize
diff --git a/src/disk.jl b/src/disk.jl
index c7e4a0c..2622d3f 100644
--- a/src/disk.jl
+++ b/src/disk.jl
@@ -68,7 +68,7 @@ summary(io::IO, P::Zernike) = print(io, "Zernike($(P.a), $(P.b))")
 
 orthogonalityweight(Z::Zernike) = 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, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/π) * r^m * normalizedjacobip(ℓ, b, m+a, 2r^2-1)
 zerniker(ℓ, m, b, r) = zerniker(ℓ, m, zero(b), b, r)
 zerniker(ℓ, m, r) = zerniker(ℓ, m, zero(r), r)
 
@@ -86,7 +86,7 @@ function getindex(Z::Zernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where
     ℓ = 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θ)
+    zernikez((ℓ-m) ÷ 2, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
 end
 
 
@@ -94,6 +94,7 @@ getindex(Z::Zernike, xy::StaticVector{2}, B::BlockIndex{1}) = Z[RadialCoordinate
 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])])
 
+basis_axes(::Inclusion{<:Any,<:UnitDisk}, v) = Zernike()
 
 ###
 # Jacobi matrices
@@ -192,6 +193,11 @@ end
 # Transforms
 ###
 
+function grammatrix(Z::Zernike{T}) where T
+    @assert Z.a == Z.b == 0
+    SquareEye{T}((axes(Z,2),))
+end
+
 function grid(S::Zernike{T}, B::Block{1}) where T
     N = Int(B) ÷ 2 + 1 # matrix rows
     M = 4N-3 # matrix columns
diff --git a/test/test_diskvector.jl b/test/test_diskvector.jl
new file mode 100644
index 0000000..0778aef
--- /dev/null
+++ b/test/test_diskvector.jl
@@ -0,0 +1,69 @@
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test, ForwardDiff, StaticArrays
+using ForwardDiff: gradient
+
+
+k = 0; m = 0; n = 2
+
+Z_x = (n,m) -> (𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[1])
+Z_y = (n,m) -> (𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[2])
+
+
+𝐱 = SVector(0.1,0.2)
+r = norm(𝐱); θ = atan(𝐱[2], 𝐱[1])
+z = 2r^2 - 1
+zernikez(n,m,𝐱)
+
+W = (n,a,b) -> 2^(a+b+1)/(2n+a+b+1) * gamma(n+a+1)gamma(n+b+1)/(gamma(n+a+b+1)factorial(n))
+
+@test jacobip(n,k, m,z) / sqrt(W(n,k,m)) ≈ normalizedjacobip(n,k,m,z)
+
+r^m * cos(m*θ) * jacobip(n,k, m,z) / sqrt(W(n,k,m) / 2^(2+k+m))
+
+sqrt(π) * zernikez(n,m,𝐱)
+
+@time expand(Zernike(), Z_x(3,2))
+
+# vector OPs
+o = expand(Zernike(), _ -> 1)
+
+v = [[o,0*o], [0*o,o]]
+ip = (v,w) -> dot(v[1],w[1]) + dot(v[2],w[2])
+ip(v[1],v[2])
+
+expand(Zernike(), Z_x(3,2))
+
+
+W_x = (n,m) -> (𝐱 -> gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), 𝐱)[1])
+W_y = (n,m) -> (𝐱 -> gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), 𝐱)[2])
+
+∇W = (n,m) -> [expand(Zernike(), W_x(n,m)),expand(Zernike(), W_y(n,m))]
+
+ip(∇W(2,3), [expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3)),expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3))])
+
+ip(,[expand(Zernike(), W_x(3,2)),expand(Zernike(), W_y(3,2))])
+
+w = [expand(Zernike(), splat((x,y)->1-y^2)) expand(Zernike(), splat((x,y)->x*y)); expand(Zernike(), splat((x,y)->x*y)) expand(Zernike(), splat((x,y)->1-x^2))]
+
+wiW1 = (n,m) -> expand(Zernike()[:,Block.(1:20)], splat((x,y) -> [1-x^2,-x*y]' * gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), SVector(x,y))/(1-x^2-y^2)))
+wiW2 = (n,m) -> expand(Zernike()[:,Block.(1:20)], splat((x,y) -> [-x*y,1-y^2]' * gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), SVector(x,y))/(1-x^2-y^2)))
+
+[wiW1(3,4),wiW2(3,4)], [expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3)),expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3))]
+
+v = [wiW1(3,4),wiW2(3,4)]
+[ip(v, ∇W(n,m)) for n=0:5, m=0:5]
+dot(∇W(0,6)[1], ∇W(4,6)[1])
+dot(∇W(0,6)[2], ∇W(4,6)[2])
+
+
+∇W(0,6)[1][SVector(0.1,0.2)]
+gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(0,6,1,𝐱), SVector(0.1,0.2))
+ip(∇W(8,4), ∇W(9,4))
+[ip(∇W(8,4), ∇W(n,m)) for n=0:10, m=0:6]
+v = [wiW1(3,4),wiW2(3,4)]
+[ip(v, ∇W(n,m)) for n=0:10, m=0:6]
+
+zernikez(4,6,1,0.1*SVector(cos(0.2),sin(0.2)))
+v[1][SVector(0.1,0.2)]
+
+(1-x^2) *P^(1,1) * (1-x^2) *P^(1,1)
+(𝐱 -> (1-norm(𝐱)^2)*zernikez(3,4,1,𝐱))(0.1*SVector(cos(0.2),sin(0.2)))
\ No newline at end of file