diff --git a/test/side_tests.jl b/test/side_tests.jl
index a1665c7..f786a15 100644
--- a/test/side_tests.jl
+++ b/test/side_tests.jl
@@ -8,7 +8,7 @@ function fit_recursion_oneterm(l1,m1,μ1,l2,m2,μ2)
    local theta, dset
    for i = 1:N_samp
       theta = rand() * 2pi
-      dset = [ Ctran(L) * wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
+      dset = [ Ctran(L) * WignerD.wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
 
       D(l,m,μ) = dset[l+1][m+l+1,μ+l+1]
       for (j,L) in enumerate(λset)
@@ -26,7 +26,7 @@ function fit_recursion_oneterm(l1,m1,μ1,l2,m2,μ2)
    
    for i = 1:N_samp
       theta = rand() * 2pi
-      dset = [ Ctran(L) * wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
+      dset = [ Ctran(L) * WignerD.wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
 
       D(l,m,μ) = dset[l+1][m+l+1,μ+l+1]
       for (j,L) in enumerate(λset)
@@ -50,7 +50,7 @@ function fit_recursion(l1,m1,μ1,l2,m2,μ2)
 
    for i = 1:N_samp
       theta = rand() * 2pi
-      dset = [ Ctran(L) * wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
+      dset = [ Ctran(L) * WignerD.wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
 
       D(l,m,μ) = dset[l+1][m+l+1,μ+l+1]
       for (j,L) in enumerate(λset)
@@ -111,7 +111,7 @@ function fit_recursion(l1,m1,μ1,l2,m2,μ2)
    
    for i = 1:N_samp
       theta = rand() * 2pi
-      dset = [ Ctran(L) * wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
+      dset = [ Ctran(L) * WignerD.wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
 
       D(l,m,μ) = dset[l+1][m+l+1,μ+l+1]
       for (j,L) in enumerate(λset)
@@ -203,9 +203,9 @@ end
 @info("1. Only when abs(m) == abs(μ) can the corresponding element be nonzero")
 for i = 1:50
    theta = rand() * 2pi
-   dset = [ Ctran(L) * wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
+   dset = [ Ctran(L) * WignerD.wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
    D(l,m,μ) = dset[l+1][m+l+1,μ+l+1]
-   # dset = [ wignerD(L, 0, 0, theta) for L = 0:Lmax ]
+   # dset = [ WignerD.wignerD(L, 0, 0, theta) for L = 0:Lmax ]
    
    l = rand(0:Lmax)
    m = rand(-l:l)
@@ -217,8 +217,8 @@ println()
 @info("2. D(l,m,μ) = D(l,-μ,-m)")
 for i = 1:50
    theta = rand() * 2pi
-   dset = [ Ctran(L) * wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
-   # dset = [ wignerD(L, 0, 0, theta) for L = 0:Lmax ]
+   dset = [ Ctran(L) * WignerD.wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
+   # dset = [ WignerD.wignerD(L, 0, 0, theta) for L = 0:Lmax ]
    D(l,m,μ) = dset[l+1][m+l+1,μ+l+1]
 
    l = rand(0:Lmax)
@@ -231,8 +231,8 @@ println()
 @info("3. D(l,m,μ) = -(-1)^{m==μ} * D(l,-m,-μ) - this is the most important part to simplify the recursion")
 for i = 1:50
    theta = rand() * 2pi
-   dset = [ Ctran(L) * wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
-   # dset = [ wignerD(L, 0, 0, theta) for L = 0:Lmax ]
+   dset = [ Ctran(L) * WignerD.wignerD(L, 0, 0, theta) * Ctran(L)' for L = 0:Lmax ]
+   # dset = [ WignerD.wignerD(L, 0, 0, theta) for L = 0:Lmax ]
    D(l,m,μ) = dset[l+1][m+l+1,μ+l+1]
 
    l = rand(0:Lmax)
diff --git a/test/test_obsolete_cg.jl b/test/test_obsolete_cg.jl
index 5687671..58fecc4 100644
--- a/test/test_obsolete_cg.jl
+++ b/test/test_obsolete_cg.jl
@@ -48,7 +48,6 @@ cg = ClebschGordan()
 # end
 
 ##
-index_y(l, m) = m + l + (l*l) + 1
 @info("Checking the SphH expansion in terms of CG coeffs")
 # expansion coefficients of a product of two spherical harmonics in terms a
 # single spherical harmonic
@@ -65,7 +64,8 @@ for ntest = 1:200
       R = SVector( cos(φ)*sin(θ), sin(φ)*sin(θ), cos(θ) )
       # evaluate all relevant Ylms (up to l1 + l2)
       
-      local Ylm =eval_cY(SphericalHarmonics(maximum(l1+l2)),R)
+      basis = complex_sphericalharmonics(maximum(l1+l2))
+      local Ylm = basis(R)
       # local Ylm = evaluate(CYlmBasis(l1+l2), R)
       # evaluate the product p = Y_l1_m1 * Y_l2_m2
       p = Ylm[index_y(l1,  m1)] * Ylm[index_y(l2,m2)]
diff --git a/test/test_obsolete_o3.jl b/test/test_obsolete_o3.jl
index 12ee575..2511119 100644
--- a/test/test_obsolete_o3.jl
+++ b/test/test_obsolete_o3.jl
@@ -35,7 +35,7 @@ for ntest = 1:30
    for L = 0:Lmax
       YL = Ylm[L^2+1:(L+1)^2]
       YrL = Ylm_r[L^2+1:(L+1)^2]
-      D = transpose(wignerD(L, θ...)) 
+      D = transpose(WignerD.wignerD(L, θ...)) 
       # Here transpose is needed because of the different convention (between ours and WignerD.jl's) of the D-matrix; 
       # Same for the cases below.
       print_tf(@test norm(D * YrL- collect(YL)) < 1e-12)
@@ -46,7 +46,7 @@ println()
 @info("The cSH D-matrix recursion")
 Lmax = 4
 θ = rand(3) * 2pi
-Dset = [ transpose(wignerD(L, θ...)) for L = 0:Lmax ]
+Dset = [ transpose(WignerD.wignerD(L, θ...)) for L = 0:Lmax ]
 
 for ntest = 1:100
    l1 = rand(0:Lmax)
@@ -80,7 +80,7 @@ for ntest = 1:30
    for L = 0:Lmax
       YL = Ylm[L^2+1:(L+1)^2]
       YrL = Ylm_r[L^2+1:(L+1)^2]
-      D = transpose(wignerD(L, θ...))
+      D = transpose(WignerD.wignerD(L, θ...))
       print_tf(@test norm(Ctran(L) * D * Ctran(L)' * YrL- collect(YL)) < 1e-12)
    end
 end
@@ -95,8 +95,8 @@ include("side_tests.jl")
 @info("The rSH D-matrix recursion")
 Lmax = 4
 θ = rand(3) * 2pi
-Dset = [ Ctran(L) * transpose(wignerD(L, θ...)) * Ctran(L)' for L = 0:Lmax ]
-DDset = [ transpose(wignerD(L, θ...)) for L = 0:Lmax ]
+Dset = [ Ctran(L) * transpose(WignerD.wignerD(L, θ...)) * Ctran(L)' for L = 0:Lmax ]
+DDset = [ transpose(WignerD.wignerD(L, θ...)) for L = 0:Lmax ]
 
 for ntest = 1:100
    l1 = rand(0:Lmax)
@@ -159,7 +159,7 @@ for L = 0:2
          if L == 0
             print_tf(@test norm(B1 - B2)<1e-12)
          else
-            D = transpose(wignerD(L, θ...))
+            D = transpose(WignerD.wignerD(L, θ...))
             print_tf(@test norm(B1 - Ref(D) .* B2)<1e-12)
          end
       end
@@ -190,7 +190,7 @@ for L = 0:0
          if L == 0
             print_tf(@test norm(B1 - B2)<1e-12)
          else
-            D = Ctran(L) * transpose(wignerD(L, θ...)) * Ctran(L)'
+            D = Ctran(L) * transpose(WignerD.wignerD(L, θ...)) * Ctran(L)'
             print_tf(@test norm(B1 - Ref(D) .* B2)<1e-12)
          end
       end
@@ -220,7 +220,7 @@ for L = 0:2
 
          B1 = eval_basis(ll, Ure, Mll, X; Real = false)
          B2 = eval_basis(ll, Ure, Mll, Ref(Q) .* X; Real = false)
-         D = BlockDiagonal([ transpose(wignerD(l, θ...)) for l = 0:L] )
+         D = BlockDiagonal([ transpose(WignerD.wignerD(l, θ...)) for l = 0:L] )
          print_tf(@test norm(B1 - Ref(D) .* B2)<1e-12)
       end
       println()
@@ -249,7 +249,7 @@ for ntest = 1:30
    for l = 0:Lmax
       B1l = [ B1[i][Val(l)] for i = 1:length(B1) ]
       B2l = [ B2[i][Val(l)] for i = 1:length(B2) ]
-      D = transpose(wignerD(l, θ...))
+      D = transpose(WignerD.wignerD(l, θ...))
       print_tf(@test norm(B1l - Ref(D) .* B2l)<1e-12)
    end
 end