rotate and translate the path in path_init to speed things up

1 files changed, 63 insertions, 26 deletions

diff --git a/src/path.c b/src/path.c
index dff1dc7..1e53d79 100644
--- a/src/path.c
+++ b/src/path.c
@@ -33,6 +33,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
 {
     size_t i, j;
     double E, mom, beta, theta, phi;
+    double normal[3], k[3], tmp[3], tmp2[3];
 
     path *p = malloc(sizeof(path));
 
@@ -87,6 +88,62 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
         }
     }
 
+    /* Now, we rotate and translate the path so that it starts at `pos` and
+     * points in the direction `dir`. */
+
+    /* We need to compute an arbitrary vector which is orthogonal to the
+     * direction vector. To do this, all we need is another vector not parallel
+     * to the direction. */
+    k[0] = 0.0;
+    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);
+
+    /* Make sure it's normalized. */
+    normalize(normal);
+
+    /* Compute the angle required to rotate the unit vector to `dir`. */
+    phi = acos(DOT(k,dir));
+
+    /* Now we rotate and translate all the positions and rotate all the
+     * directions. */
+    for (i = 0; i < N; i++) {
+        tmp[0] = x[i];
+        tmp[1] = y[i];
+        tmp[2] = z[i];
+
+        rotate(tmp2,tmp,normal,phi);
+
+        ADD(tmp2,tmp2,pos);
+
+        x[i] = tmp2[0];
+        y[i] = tmp2[1];
+        z[i] = tmp2[2];
+
+        tmp[0] = dx[i];
+        tmp[1] = dy[i];
+        tmp[2] = dz[i];
+
+        normalize(tmp);
+
+        rotate(tmp2,tmp,normal,phi);
+
+        dx[i] = tmp2[0];
+        dy[i] = tmp2[1];
+        dz[i] = tmp2[2];
+    }
+
     for (i = 0; i < 8; i++) {
         p->acc[i] = gsl_interp_accel_alloc();
         p->spline[i] = gsl_spline_alloc(gsl_interp_linear,N);
@@ -118,39 +175,19 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
 
 void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0)
 {
-    double normal[3], k[3], tmp[3], phi;
-
     if (s > p->spline[0]->x[p->spline[0]->size-1])
         s = p->spline[0]->x[p->spline[0]->size-1];
 
-    tmp[0] = gsl_spline_eval(p->spline[0],s,p->acc[0]);
-    tmp[1] = gsl_spline_eval(p->spline[1],s,p->acc[1]);
-    tmp[2] = gsl_spline_eval(p->spline[2],s,p->acc[2]);
+    pos[0] = gsl_spline_eval(p->spline[0],s,p->acc[0]);
+    pos[1] = gsl_spline_eval(p->spline[1],s,p->acc[1]);
+    pos[2] = gsl_spline_eval(p->spline[2],s,p->acc[2]);
 
     *T = gsl_spline_eval(p->spline[3],s,p->acc[3]);
     *t = gsl_spline_eval(p->spline[4],s,p->acc[4]);
 
-    k[0] = 0.0;
-    k[1] = 0.0;
-    k[2] = 1.0;
-
-    CROSS(normal,k,p->dir);
-
-    normalize(normal);
-
-    phi = acos(DOT(k,p->dir));
-
-    rotate(pos,tmp,normal,phi);
-
-    ADD(pos,pos,p->pos);
-
-    tmp[0] = gsl_spline_eval(p->spline[5],s,p->acc[5]);
-    tmp[1] = gsl_spline_eval(p->spline[6],s,p->acc[6]);
-    tmp[2] = gsl_spline_eval(p->spline[7],s,p->acc[7]);
-
-    normalize(tmp);
-
-    rotate(dir,tmp,normal,phi);
+    dir[0] = gsl_spline_eval(p->spline[5],s,p->acc[5]);
+    dir[1] = gsl_spline_eval(p->spline[6],s,p->acc[6]);
+    dir[2] = gsl_spline_eval(p->spline[7],s,p->acc[7]);
 
     /* FIXME: This should be the *residual* scattering RMS. */
     *theta0 = p->theta0*sqrt(s);