fix a bug in path_init() when the direction was equal to (0,0,1)

2 files changed, 28 insertions, 21 deletions

diff --git a/src/path.c b/src/path.c
index fd5d2a3..3904b44 100644
--- a/src/path.c
+++ b/src/path.c
@@ -104,20 +104,21 @@ path *path_init(double *pos, double *dir, double T0, double range, size_t len, d
     k[1] = 0.0;
     k[2] = 1.0;
 
-    /* If k is approximately equal to the unit vector in the z direction, we
-     * switch to the unit vector in the x direction. */
-    if (allclose(k,dir,3,1e-9,1e-9)) {
-        k[0] = 1.0;
-        k[1] = 0.0;
-        k[2] = 0.0;
-    }
-
-    /* Take the cross product between `k` and `dir` to get a vector orthogonal
-     * to `dir`. */
-    CROSS(normal,k,dir);
+    if (allclose(k,dir,3,0,0)) {
+        /* If the direction vector is already in the z direction, we don't need
+         * to rotate it, so the normal direction is arbitrary since `phi` will
+         * be 0. */
+        normal[0] = 1.0;
+        normal[1] = 0.0;
+        normal[2] = 0.0;
+    } else {
+        /* Take the cross product between `k` and `dir` to get the axis of
+         * rotation. */
+        CROSS(normal,k,dir);
 
-    /* Make sure it's normalized. */
-    normalize(normal);
+        /* Make sure it's normalized. */
+        normalize(normal);
+    }
 
     /* Compute the angle required to rotate the unit vector to `dir`. */
     phi = acos(DOT(k,dir));
diff --git a/src/test.c b/src/test.c
index 5c44243..3bd8c5c 100644
--- a/src/test.c
+++ b/src/test.c
@@ -488,6 +488,7 @@ int test_path(char *err)
     double E, mom, beta0;
     double s;
     double theta0;
+    path *p;
 
     T0 = 1.0;
     range = 1.0;
@@ -495,17 +496,17 @@ int test_path(char *err)
 
     init_genrand(0);
 
-    for (i = 0; i < 100; i++) {
+    for (i = 0; i < 100000; i++) {
         pos0[0] = genrand_real2();
         pos0[1] = genrand_real2();
         pos0[2] = genrand_real2();
 
         rand_sphere(dir0);
 
-        path *p = path_init(pos0,dir0,T0,range,1000,0.1,getKineticEnergy,NULL,NULL,NULL,0,m);
+        p = path_init(pos0,dir0,T0,range,1000,0.1,getKineticEnergy,NULL,NULL,NULL,0,m);
 
-        for (j = 0; j < 100; j++) {
-            s = range*j/99;
+        for (j = 0; j < p->len; j++) {
+            s = range*j/(p->len-1);
             pos_expected[0] = pos0[0] + dir0[0]*s;
             pos_expected[1] = pos0[1] + dir0[1]*s;
             pos_expected[2] = pos0[2] + dir0[2]*s;
@@ -522,25 +523,25 @@ int test_path(char *err)
             for (k = 0; k < 3; k++) {
                 if (!isclose(pos[k], pos_expected[k], 1e-9, 1e-9)) {
                     sprintf(err, "path_eval(%.2g) returned pos[%i] = %.5g, but expected %.5g", s, k, pos[k], pos_expected[k]);
-                    return 1;
+                    goto err;
                 }
             }
 
             for (k = 0; k < 3; k++) {
                 if (!isclose(dir[k], dir0[k], 1e-9, 1e-9)) {
                     sprintf(err, "path_eval(%.2g) returned dir[%i] = %.5g, but expected %.5g", s, k, dir[k], dir0[k]);
-                    return 1;
+                    goto err;
                 }
             }
 
             if (!isclose(beta, beta0, 1e-9, 1e-9)) {
                 sprintf(err, "path_eval(%.2g) returned beta = %.5g, but expected %.5g", s, beta, beta0);
-                return 1;
+                goto err;
             }
 
             if (!isclose(t, t_expected, 1e-9, 1e-9)) {
                 sprintf(err, "path_eval(%.2g) returned t = %.5g, but expected %.5g", s, t, t_expected);
-                return 1;
+                goto err;
             }
         }
 
@@ -548,6 +549,11 @@ int test_path(char *err)
     }
 
     return 0;
+
+err:
+    path_free(p);
+
+    return 1;
 }
 
 int test_interp1d(char *err)