add weighted partial derivatives

actually test rectdisk
import AngularMomentum

Author: MikaelSlevinsky

src/rectdisk.jl

@@ -57,6 +57,28 @@ end
     DunklXuDisk(β+1) * _BandedBlockBandedMatrix(((k .+ T(2β)) ./ 2)', axes(k,1), (-1,1), (-1,1))
 end
 
+@simplify function *(Dx::PartialDerivative{1}, w_P::WeightedDunklXuDisk)
+    wP, P = w_P.args
+    @assert P.β == wP.β
+    β = P.β
+    n = mortar(Fill.(oneto(∞),oneto(∞)))
+    k = mortar(Base.OneTo.(oneto(∞)))
+    dat = BlockBroadcastArray(hcat,
+        (-4 .* (k .+ (β - 1)).*(n .- k .+ 1)./(2k .+ (2β - 1))),
+        0 .* n,
+        (-k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (k .+ β)) .* (n .+ k .+ 2β))
+        )
+    WeightedDunklXuDisk(β-1) * _BandedBlockBandedMatrix(dat', axes(k,1), (1,-1), (2,0))
+end
+
+@simplify function *(Dy::PartialDerivative{2}, w_P::WeightedDunklXuDisk)
+    wP, P = w_P.args
+    @assert P.β == wP.β
+    k = mortar(Base.OneTo.(oneto(∞)))
+    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
+    WeightedDunklXuDisk(P.β-1) * _BandedBlockBandedMatrix((T(-2).*k)', axes(k,1), (1,-1), (1,-1))
+end
+
 # P^{(β)} ↗ P^{(β+1)}
 function dunklxu_raising(β)
     n = mortar(Fill.(oneto(∞),oneto(∞)))

test/runtests.jl

@@ -2,5 +2,6 @@ using MultivariateOrthogonalPolynomials, Test
 
 include("test_modalinterlace.jl")
 include("test_disk.jl")
+include("test_rectdisk.jl")
 include("test_triangle.jl")
 # include("test_dirichlettriangle.jl")

test/test_rectdisk.jl

@@ -1,6 +1,6 @@
 using MultivariateOrthogonalPolynomials, StaticArrays, BlockArrays, BlockBandedMatrices, ArrayLayouts,
         QuasiArrays, Test, ClassicalOrthogonalPolynomials, BandedMatrices, FastTransforms, LinearAlgebra
-import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering
+import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, AngularMomentum
 
 @testset "Dunkl-Xu disk" begin
     @testset "basics" begin
@@ -34,8 +34,8 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering
 
         @test (L*R)[Block.(1:N), Block.(1:N)] ≈ (I - X^2 - Y^2)[Block.(1:N), Block.(1:N)]
 
-        ∂x = PartialDerivative{1}(P)
-        ∂y = PartialDerivative{2}(P)
+        ∂x = PartialDerivative{1}(axes(P, 1))
+        ∂y = PartialDerivative{2}(axes(P, 1))
 
         Dx = Q \ (∂x * P)
         Dy = Q \ (∂y * P)
@@ -61,6 +61,12 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering
 
         @test A[Block.(1:N), Block.(1:N)] ≈ C
 
+        ∂x = PartialDerivative{1}(axes(WQ, 1))
+        ∂y = PartialDerivative{2}(axes(WQ, 1))
+
+        wDx = WP \ (∂x * WQ)
+        wDy = WP \ (∂y * WQ)
+
         @testset "truncations" begin
             KR,JR = Block.(1:N),Block.(1:N)