using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, RecipesBase, Test, SpecialFunctions
import MultivariateOrthogonalPolynomials: DiskTrav, grid, ZernikeTransform, ZernikeITransform, *
import ClassicalOrthogonalPolynomials: HalfWeighted


@testset "Disk" 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
        

        @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 "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

    @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
    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))

        xy = SVector(0.1,0.2)
        @test Zernike()[xy,Block(1)[1]] ≈ Zernike(1)[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 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)

        @test Zernike()[xy,Block(5)[1]] ≈ R0[2:3,3]'*Zernike(1)[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)


        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)]
    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)))

        xy = SVector(0.1,0.2)
        r = norm(xy)
        w = 1 - r^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*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*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)

        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 exp.(L)[1:10,1:10] == exp.(L[1:10,1:10])

        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 "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

    @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