diff --git a/src/AstroBasis.jl b/src/AstroBasis.jl
index 104066a..3c8b23e 100644
--- a/src/AstroBasis.jl
+++ b/src/AstroBasis.jl
@@ -27,9 +27,11 @@ export getUln,getDln,tabUl!,tabDl!
 
 # bring in the Clutton-Brock (1973) spherical basis
 include("CB73.jl")
+export CB73Basis
 
 # bring in the Hernquist & Ostriker (1992) disc basis
 include("Hernquist.jl")
+export HernquistBasis
 
 # @IMPROVE, add Bessel basis
 
@@ -37,9 +39,11 @@ include("Hernquist.jl")
 
 # bring in the Clutton-Brock (1972) disc basis
 include("CB72.jl")
+export CB72Basis
 
 # bring in the Kalnajs (1976) disc basis
 include("Kalnajs76.jl")
+export K76Basis
 
 
 end # module
diff --git a/src/Hernquist.jl b/src/Hernquist.jl
index 977a399..34951b6 100644
--- a/src/Hernquist.jl
+++ b/src/Hernquist.jl
@@ -134,12 +134,12 @@ end
 
 
 
-"""ClnHernquist(α,n,rho)
+"""_ClnHernquist(α,n,rho)
 Definition of the Gegenbauer polynomials
 These coefficients are computed through an upward recurrence
 @IMPROVE compute all the basis elements (n)_{1<=n<=nmax} at once
 """
-function ClnHernquist(α::Float64,n::Int64,xi::Float64)
+function _ClnHernquist(α::Float64,n::Int64,xi::Float64)
 
     v0 = 1.0 # Initial value for n=0
     if (n == 0)
@@ -173,7 +173,7 @@ function UlnpHernquist(l::Int64,np::Int64,r::Float64,tabPrefHernquist_Ulnp::Matr
     x    = r/rb                                    # Dimensionless radius
     rho  = rhoHernquist(x)                              # Value of the rescaled parameter rho
     valR = ((x/(1.0+x^(2)))^(l))/(sqrt(1.0+x^(2))) # Value of the multipole factor
-    valC = ClnHernquist(l+1.0,np,rho)                 # Value of the Gegenbauer polynomials
+    valC = _ClnHernquist(l+1.0,np,rho)                 # Value of the Gegenbauer polynomials
     res  = pref*valR*valC                          # Value of the radial function
     return res
 end
@@ -193,7 +193,7 @@ function DlnpHernquist(l::Int64,np::Int64,r::Float64,tabPrefHernquist_Dlnp::Matr
     x    = r/rb                                      # Dimensionless radius
     rho  = rhoHernquist(x)                                # Value of the rescaled parameter rho
     valR = ((x/(1.0+x^(2)))^(l))/((1.0+x^(2))^(5/2)) # Value of the multipole factor
-    valC = ClnHernquist(l+1.0,np,rho)                   # Value of the Gegenbauer polynomials
+    valC = _ClnHernquist(l+1.0,np,rho)                   # Value of the Gegenbauer polynomials
     res  = pref*valR*valC
     return res
 end
diff --git a/test/test_razorthin.jl b/test/test_razorthin.jl
index ed91490..b0da5c1 100644
--- a/test/test_razorthin.jl
+++ b/test/test_razorthin.jl
@@ -5,9 +5,10 @@ ltest, nradial = 2, 5
 
 @testset "razorthinbases" begin
     @testset "CB72" begin
-        basis = AstroBasis.CB72Basis(lmax=ltest,nradial=nradial,G=G, rb=rb)
+        basis = CB72Basis(lmax=ltest,nradial=nradial,G=G, rb=rb)
         @test dimension(basis) == 2
         @test getparameters(basis)["name"] == "CB72"
+        @test getDln(basis,ltest,1,rb) == 0.0
         @test getDln(basis,ltest,nradial-1,rb) ≈ 0.33126698841066865 atol=1e-6
         @test getUln(basis,ltest,nradial-1,rb) ≈ -0.3202172114362374 atol=1e-6
         tabUl!(basis,0,rb)
@@ -16,9 +17,10 @@ ltest, nradial = 2, 5
         @test basis.tabDl[2] == 0
     end
     @testset "Kalnajs76" begin
-        basis = AstroBasis.K76Basis(lmax=ltest,nradial=nradial,G=G, rb=rb)
+        basis = K76Basis(lmax=ltest,nradial=nradial,G=G, rb=rb)
         @test dimension(basis) == 2
         @test getparameters(basis)["name"] == "K76"
+        @test getDln(basis,ltest,1,rb) == 0.0
         @test getDln(basis,ltest,nradial-1,rb) == 0.0
         @test getUln(basis,ltest,nradial-1,rb) ≈ -0.3167679231608087 atol=1e-6
         tabUl!(basis,0,rb)
diff --git a/test/test_spherical.jl b/test/test_spherical.jl
index 07f83cc..3e84471 100644
--- a/test/test_spherical.jl
+++ b/test/test_spherical.jl
@@ -6,11 +6,11 @@ ltest, nradial = 2, 5
 
 @testset "sphericalbases" begin
     @testset "CB73" begin
-        # build bases with unique cases
-        basis = AstroBasis.CB73Basis(lmax=ltest,nradial=1,G=G, rb=rb)
-        basis = AstroBasis.CB73Basis(lmax=ltest,nradial=2,G=G, rb=rb)
+        # test Gegenbauer polynomial (non)recursions
+        @test AstroBasis._ClnCB73(1.0,0,1.0) == 1.0
+        @test AstroBasis._ClnCB73(1.0,1,1.0) == 2.0
         # build a standard basis
-        basis = AstroBasis.CB73Basis(lmax=ltest,nradial=nradial,G=G, rb=rb)
+        basis = CB73Basis(lmax=ltest,nradial=nradial,G=G, rb=rb)
         @test dimension(basis) == 3
         @test getparameters(basis)["name"] == "CB73"
         # backward compatibility check
@@ -25,7 +25,10 @@ ltest, nradial = 2, 5
         AstroBasis.WriteParameters("tmp.h5",basis,"w")
     end
     @testset "Hernquist" begin
-        basis = AstroBasis.HernquistBasis(lmax=ltest,nradial=nradial,G=G, rb=rb)
+        # test Gegenbauer polynomial (non)recursions
+        @test AstroBasis._ClnHernquist(1.0,0,1.0) == 1.0
+        @test AstroBasis._ClnHernquist(1.0,1,1.0) == 2.0
+        basis = HernquistBasis(lmax=ltest,nradial=nradial,G=G, rb=rb)
         @test dimension(basis) == 3
         @test getparameters(basis)["name"] == "Hernquist"
         @test getDln(basis,ltest,nradial-1,rb) ≈ 1.552003244163087 atol=1e-6