diff --git a/polyrat/basis.py b/polyrat/basis.py
index 3353f8f..c0fd5ed 100644
--- a/polyrat/basis.py
+++ b/polyrat/basis.py
@@ -110,7 +110,6 @@ def _scale(self, X):
 
 
 class LaguerrePolynomialBasis(TensorProductPolynomialBasis):
-	# TODO: Change scaling to match orthogonality conditions
 	def _vander(self, *args, **kwargs):
 		return lagvander(*args, **kwargs)
 	def _der(self, *args, **kwargs):
@@ -119,6 +118,9 @@ def _roots(self, *args, **kwargs):
 		return lagroots(*args, **kwargs)
 
 	def _set_scale(self, X):
+		r""" Laguerre polynomial expects x[i] to be distributed like exp[-lam*x] on [0,infty)
+		so we shift so that all entries are positive and then set a scaling
+		"""
 		self._lb = np.min(X, axis = 0)
 		self._a = 1./np.mean(X - self._lb[None,:], axis = 0)
 		
@@ -243,8 +245,8 @@ def vandermonde_arnoldi_eval(X, R, indices, mode, weight = None):
 	W[:,0] = weight/R[0,0]
 	
 	# Now work on the remaining columns
-	for k, ids in iteridx:
-		i, j = _update_vec(idx[:k], ids)
+	for k, ids in iter_indices:
+		i, j = update_rule(indices[:k], ids)
 		# Form new column
 		w = X[:,i] * W[:,j]
 
@@ -286,5 +288,5 @@ def basis(self):
 		return self._Q   
 
 	def vandermonde(self, X, weight = None):
-		return vanderonde_arnoldi_eval(X, self._R, self._indices, self.mode, weight = weight)
+		return vandermonde_arnoldi_eval(X, self._R, self._indices, self.mode, weight = weight)
 	
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 59cd982..0debaf5 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 from polyrat import *
-
+from scipy.linalg import subspace_angles
 
 
 
@@ -75,8 +75,97 @@ def test_monomial_max(M, degree, Basis, seed):
 		assert err_norm < 1e-8, "Vector not contained in space"
 
 
+@pytest.mark.parametrize("n_grid", [10])
+@pytest.mark.parametrize("degree",[
+	5,
+	(7,2),
+	(2,7),
+	(3,0,4)
+	])
+@pytest.mark.parametrize("dim", [1,2,3])
+@pytest.mark.parametrize("Basis", 
+	[MonomialPolynomialBasis,
+	 ChebyshevPolynomialBasis,
+	 HermitePolynomialBasis, 
+	 LaguerrePolynomialBasis,
+	 ArnoldiPolynomialBasis])
+def test_subspace_angles_basis(n_grid, degree, dim,  Basis):
+
+	# Exit without error if testing a total degree problem
+	try:
+		degree = int(degree)
+	except (TypeError, ValueError):
+		if len(degree) != dim: return
+
+	# Construct grid 
+	Xs = np.meshgrid(*[np.linspace(-1,1,n_grid) for i in range(dim)])
+	X = np.vstack([Xi.flatten() for Xi in Xs]).T
+	
+	
+	# We use Legendre polynomials as a reference solution
+	leg = LegendrePolynomialBasis(X, degree)
+	basis = Basis(X, degree)
+	
+	V1 = leg.basis()
+	V2 = basis.basis()
+
+	tol = 1e-10
+	# This Laguerre polynomials are not particularly well conditioned
+	if Basis is LaguerrePolynomialBasis: tol = 1e-7
+
+	for k in range(1,V1.shape[1]):
+		ang = subspace_angles(V1[:,:k], V2[:,:k])
+		print(basis._indices[k], np.max(ang))
+		assert np.max(ang) < tol
+
+
+
+@pytest.mark.parametrize("seed", [0])
+@pytest.mark.parametrize("degree",[
+	5,
+	(7,2),
+	(2,7),
+	(3,0,4)
+	])
+@pytest.mark.parametrize("dim", [1,2,3])
+@pytest.mark.parametrize("Basis", 
+	[MonomialPolynomialBasis,
+	 ChebyshevPolynomialBasis,
+	 HermitePolynomialBasis, 
+	 LaguerrePolynomialBasis,
+	 ArnoldiPolynomialBasis])
+def test_subspace_angles(seed, degree, dim,  Basis):
+
+	# Exit without error if testing a total degree problem
+	try:
+		degree = int(degree)
+	except (TypeError, ValueError):
+		if len(degree) != dim: return
+
+	# Construct grid 
+	X = np.random.randn(500, dim)
+
+	X2 = np.random.randn(200, dim)	
+	# We use Legendre polynomials as a reference solution
+	leg = LegendrePolynomialBasis(X, degree)
+	basis = Basis(X, degree)
+	
+	V1 = leg.vandermonde(X2)
+	V2 = basis.vandermonde(X2)
+
+	tol = 1e-10
+	# This Laguerre polynomials are not particularly well conditioned
+	if Basis is LaguerrePolynomialBasis: tol = 1e-7
+
+	for k in range(1,V1.shape[1]):
+		ang = subspace_angles(V1[:,:k], V2[:,:k])
+		print(basis._indices[k], np.max(ang))
+		assert np.max(ang) < tol
+
+
 
 
 if __name__ == '__main__':
 	#test_monomial_total(100, 10, 3, MonomialPolynomialBasis,0) 
-	test_monomial_max(100, [4], MonomialPolynomialBasis,0) 
+	#test_monomial_max(100, [4], MonomialPolynomialBasis,0) 
+	test_subspace_angles(10, [1,4], 2, ArnoldiPolynomialBasis)