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

1 files changed, 14 insertions, 13 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));