diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org
index f997913..93d9b92 100644
--- a/src/qmckl_ao.org
+++ b/src/qmckl_ao.org
@@ -485,7 +485,6 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
   Y(:) = X(:) - R(:)
 
   lmax = 4;
-  n = 0;
   ldl = 3;
   ldv = 100;
 
@@ -495,23 +494,26 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
 
   test_qmckl_ao_polynomial_vgl = &
        qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)
-  if (test_qmckl_ao_polynomial_vgl /= 0) return
 
-  test_qmckl_ao_polynomial_vgl = -1
+  if (test_qmckl_ao_polynomial_vgl /= QMCKL_INVALID_ARG_5) return
 
-  if (n /= d) return 
+  n=d
+  test_qmckl_ao_polynomial_vgl = &
+       qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)
+
+  if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return
 
   do j=1,n
-     test_qmckl_ao_polynomial_vgl = -11
+     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
      do i=1,3
         if (L(i,j) < 0) return
      end do
-     test_qmckl_ao_polynomial_vgl = -12
+     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
      if (dabs(1.d0 - VGL(1,j) / (&
           Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j)  &
           )) > epsilon ) return
 
-     test_qmckl_ao_polynomial_vgl = -13
+     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
      if (L(1,j) < 1) then
         if (VGL(2,j) /= 0.d0) return
      else
@@ -520,7 +522,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
              )) > epsilon ) return
      end if
 
-     test_qmckl_ao_polynomial_vgl = -14
+     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
      if (L(2,j) < 1) then
         if (VGL(3,j) /= 0.d0) return
      else
@@ -529,7 +531,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
              )) > epsilon ) return
      end if
 
-     test_qmckl_ao_polynomial_vgl = -15
+     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
      if (L(3,j) < 1) then
         if (VGL(4,j) /= 0.d0) return
      else
@@ -538,7 +540,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
              )) > epsilon ) return
      end if
 
-     test_qmckl_ao_polynomial_vgl = -16
+     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
      w = 0.d0
      if (L(1,j) > 1) then
         w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j) 
@@ -552,7 +554,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
      if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
   end do
 
-  test_qmckl_ao_polynomial_vgl = 0
+  test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS
 
   deallocate(L,VGL)
 end function test_qmckl_ao_polynomial_vgl
diff --git a/src/qmckl_distance.org b/src/qmckl_distance.org
index e797393..b239621 100644
--- a/src/qmckl_distance.org
+++ b/src/qmckl_distance.org
@@ -21,19 +21,16 @@ MunitResult test_<<filename()>>() {
 
 * Squared distance
 
-** ~qmckl_distance_sq~
+  ~qmckl_distance_sq~ computes the matrix of the squared distances
+   between all pairs of points in two sets, one point within each set:
 
-   Computes the matrix of the squared distances between all pairs of
-   points in two sets, one point within each set:
    \[
    C_{ij} = \sum_{k=1}^3 (A_{k,i}-B_{k,j})^2
    \]
 
-*** Arguments
-
     | ~context~  | input  | Global state                                 |
-    | ~transa~   | input  | Array ~A~ is ~N~: Normal, ~T~: Transposed   |
-    | ~transb~   | input  | Array ~B~ is ~N~: Normal, ~T~: Transposed   |
+    | ~transa~   | input  | Array ~A~ is ~N~: Normal, ~T~: Transposed    |
+    | ~transb~   | input  | Array ~B~ is ~N~: Normal, ~T~: Transposed    |
     | ~m~        | input  | Number of points in the first set            |
     | ~n~        | input  | Number of points in the second set           |
     | ~A(lda,3)~ | input  | Array containing the $m \times 3$ matrix $A$ |
@@ -74,6 +71,7 @@ qmckl_exit_code qmckl_distance_sq(const qmckl_context context,
 *** Source
     #+begin_src f90 :tangle (eval f)
 integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) result(info)
+  use qmckl
   implicit none
   integer*8  , intent(in)  :: context
   character  , intent(in)  :: transa, transb
@@ -91,18 +89,18 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, L
 
   info = 0
 
-  if (context == 0_8) then
-     info = -1
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
      return
   endif
 
   if (m <= 0_8) then
-     info = -2
+     info = QMCKL_INVALID_ARG_4
      return
   endif
 
   if (n <= 0_8) then
-     info = -3
+     info = QMCKL_INVALID_ARG_5
      return
   endif
 
@@ -123,27 +121,27 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, L
   endif
 
   if (transab < 0) then
-     info = -4
+     info = QMCKL_INVALID_ARG_1
      return 
   endif
 
   if (iand(transab,1) == 0 .and. LDA < 3) then
-     info = -5
+     info = QMCKL_INVALID_ARG_7
      return
   endif
 
   if (iand(transab,1) == 1 .and. LDA < m) then
-     info = -6
+     info = QMCKL_INVALID_ARG_7
      return
   endif
 
   if (iand(transab,2) == 0 .and. LDA < 3) then
-     info = -6
+     info = QMCKL_INVALID_ARG_7
      return
   endif
 
   if (iand(transab,2) == 2 .and. LDA < m) then
-     info = -7
+     info = QMCKL_INVALID_ARG_7
      return
   endif
 
@@ -270,13 +268,19 @@ integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
      end do
   end do
 
-  test_qmckl_distance_sq = qmckl_distance_sq(context, 'X', 't', m, n, A, LDA, B, LDB, C, LDC)
+  test_qmckl_distance_sq = &
+       qmckl_distance_sq(context, 'X', 't', m, n, A, LDA, B, LDB, C, LDC)
+
   if (test_qmckl_distance_sq == 0) return 
 
-  test_qmckl_distance_sq = qmckl_distance_sq(context, 't', 'X', m, n, A, LDA, B, LDB, C, LDC)
+  test_qmckl_distance_sq = &
+       qmckl_distance_sq(context, 't', 'X', m, n, A, LDA, B, LDB, C, LDC)
+
   if (test_qmckl_distance_sq == 0) return 
 
-  test_qmckl_distance_sq = qmckl_distance_sq(context, 'T', 't', m, n, A, LDA, B, LDB, C, LDC)
+  test_qmckl_distance_sq = &
+       qmckl_distance_sq(context, 'T', 't', m, n, A, LDA, B, LDB, C, LDC)
+
   if (test_qmckl_distance_sq /= 0) return 
 
   test_qmckl_distance_sq = -1
@@ -290,7 +294,9 @@ integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
      end do
   end do
   
-  test_qmckl_distance_sq = qmckl_distance_sq(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)
+  test_qmckl_distance_sq = &
+       qmckl_distance_sq(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)
+
   if (test_qmckl_distance_sq /= 0) return 
 
   test_qmckl_distance_sq = -1
@@ -304,7 +310,9 @@ integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
      end do
   end do
 
-  test_qmckl_distance_sq = qmckl_distance_sq(context, 'T', 'n', m, n, A, LDA, B, LDB, C, LDC)
+  test_qmckl_distance_sq = &
+       qmckl_distance_sq(context, 'T', 'n', m, n, A, LDA, B, LDB, C, LDC)
+
   if (test_qmckl_distance_sq /= 0) return 
 
   test_qmckl_distance_sq = -1
@@ -318,7 +326,9 @@ integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
      end do
   end do
 
-  test_qmckl_distance_sq = qmckl_distance_sq(context, 'n', 'N', m, n, A, LDA, B, LDB, C, LDC)
+  test_qmckl_distance_sq = &
+       qmckl_distance_sq(context, 'n', 'N', m, n, A, LDA, B, LDB, C, LDC)
+
   if (test_qmckl_distance_sq /= 0) return 
 
   test_qmckl_distance_sq = -1