diff options
Diffstat (limited to 'src/path.c')
| -rw-r--r-- | src/path.c | 25 |
1 files changed, 13 insertions, 12 deletions
@@ -35,7 +35,7 @@ double path_get_coefficient(unsigned int k, double *s, double *x, double theta0, path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m) { size_t i, j; - double E, mom, beta, theta, phi; + double T, E, mom, theta, phi; double normal[3], k[3], tmp[3], tmp2[3]; path *p = malloc(sizeof(path)); @@ -56,7 +56,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 double *x = calloc(N,sizeof(double)); double *y = calloc(N,sizeof(double)); double *z = calloc(N,sizeof(double)); - double *T = calloc(N,sizeof(double)); + double *beta = calloc(N,sizeof(double)); double *t = calloc(N,sizeof(double)); double *dx = calloc(N,sizeof(double)); double *dy = calloc(N,sizeof(double)); @@ -69,7 +69,12 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 theta1[i] += z1[j]*foo(s[i],range,theta0,j+1); theta2[i] += z2[j]*foo(s[i],range,theta0,j+1); } - T[i] = fun(s[i],params); + + T = fun(s[i],params); + E = T + m; + mom = sqrt(E*E - m*m); + beta[i] = mom/E; + if (i > 0) { theta = sqrt(theta1[i]*theta1[i] + theta2[i]*theta2[i]); phi = atan2(theta2[i],theta1[i]); @@ -79,11 +84,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 dz[i] = (s[i]-s[i-1])*cos(theta); /* Calculate total energy */ - E = T[i] + m; - mom = sqrt(E*E - m*m); - beta = mom/E; - - t[i] = t[i-1] + (s[i]-s[i-1])/(beta*SPEED_OF_LIGHT); + t[i] = t[i-1] + (s[i]-s[i-1])/(beta[i]*SPEED_OF_LIGHT); x[i] = x[i-1] + dx[i]; y[i] = y[i-1] + dy[i]; @@ -168,7 +169,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 p->x = x; p->y = y; p->z = z; - p->T = T; + p->beta = beta; p->t = t; p->dx = dx; p->dy = dy; @@ -180,13 +181,13 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 return p; } -void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0) +void path_eval(path *p, double s, double *pos, double *dir, double *beta, double *t, double *theta0) { pos[0] = interp1d(s,p->s,p->x,N); pos[1] = interp1d(s,p->s,p->y,N); pos[2] = interp1d(s,p->s,p->z,N); - *T = interp1d(s,p->s,p->T,N); + *beta = interp1d(s,p->s,p->beta,N); *t = interp1d(s,p->s,p->t,N); /* Technically we should be interpolating between the direction vectors @@ -214,7 +215,7 @@ void path_free(path *p) free(p->x); free(p->y); free(p->z); - free(p->T); + free(p->beta); free(p->t); free(p->dx); free(p->dy); |