diff --git a/GHA_triaxial/gha1_approx.py b/GHA_triaxial/gha1_approx.py
index 042f54a..80b65e9 100644
--- a/GHA_triaxial/gha1_approx.py
+++ b/GHA_triaxial/gha1_approx.py
@@ -47,7 +47,7 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
             pass
 
     if all_points:
-        return points[-1], alphas[-1], np.array(points)
+        return points[-1], alphas[-1], np.array(points), np.array(alphas)
     else:
         return points[-1], alphas[-1]
 
@@ -82,7 +82,7 @@ if __name__ == '__main__':
     P0 = ell.para2cart(0.2, 0.3)
     alpha0 = wu.deg2rad(35)
     s = 13000000
-    P1_app, alpha1_app, points = gha1_approx(ell, P0, alpha0, s, ds=10000, all_points=True)
+    P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=10000, all_points=True)
     P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=16)
     show_points(points, P0, P1_ana)
     print(np.linalg.norm(P1_app - P1_ana))