diff --git a/Project.toml b/Project.toml
index 64654a5..a2e3044 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.9.4"
+version = "0.9.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
diff --git a/src/disk.jl b/src/disk.jl
index c7e4a0c..89574ab 100644
--- a/src/disk.jl
+++ b/src/disk.jl
@@ -57,7 +57,7 @@ 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(∞)))
+axes(P::Zernike{T}) where T = (Inclusion(UnitDisk{real(T)}()),blockedrange(oneto(∞)))
 
 ==(w::Zernike, v::Zernike) = w.a == v.a && w.b == v.b
 
@@ -192,7 +192,8 @@ end
 # Transforms
 ###
 
-function grid(S::Zernike{T}, B::Block{1}) where T
+function grid(S::Zernike, B::Block{1})
+    T = real(eltype(S))
     N = Int(B) ÷ 2 + 1 # matrix rows
     M = 4N-3 # matrix columns
 
@@ -253,10 +254,11 @@ function ZernikeITransform{T}(N::Int, a::Number, b::Number) where T<:Real
 end
 
 *(P::ZernikeTransform{T}, f::AbstractArray) where T = P * convert(Matrix{T}, f)
+*(P::ZernikeTransform{T}, f::Matrix{Complex{T}}) where T = ModalTrav((P.disk2cxf \ (P.analysis * real(f))) + im * (P.disk2cxf \ (P.analysis * imag(f))))
 *(P::ZernikeTransform{T}, f::Matrix{T}) where T = ModalTrav(P.disk2cxf \ (P.analysis * f))
 *(P::ZernikeITransform, f::AbstractVector) = P.synthesis * (P.disk2cxf * ModalTrav(f).matrix)
 
-plan_transform(Z::Zernike{T}, (N,)::Tuple{Block{1}}, dims=1) where T = ZernikeTransform{T}(Int(N), Z.a, Z.b)
+plan_transform(Z::Zernike{T}, (N,)::Tuple{Block{1}}, dims=1) where T = ZernikeTransform{real(T)}(Int(N), Z.a, Z.b)
 
 ##
 # Laplacian
diff --git a/test/test_disk.jl b/test/test_disk.jl
index 12272c4..649d6f2 100644
--- a/test/test_disk.jl
+++ b/test/test_disk.jl
@@ -68,6 +68,11 @@ import ForwardDiff: hessian
             end
         end
     end
+    
+    @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)] ≈ exp(0.1cos(0.2)+im*0.2)
+    end
 
     @testset "Jacobi matrices" begin
         # Setup