diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile | 2 | ||||
| -rw-r--r-- | src/misc.c | 24 | ||||
| -rw-r--r-- | src/misc.h | 2 | ||||
| -rw-r--r-- | src/path.c | 89 | ||||
| -rw-r--r-- | src/test.c | 9 |
5 files changed, 90 insertions, 36 deletions
diff --git a/src/Makefile b/src/Makefile index 75a75cb..2e9169a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -18,7 +18,7 @@ test-vector: test-vector.o vector.o mt19937ar.o test-likelihood: test-likelihood.o muon.o random.o optics.o quantum_efficiency.o mt19937ar.o pdg.o vector.o solid_angle.o scattering.o -test-path: test-path.o mt19937ar.o vector.o path.o random.o +test-path: test-path.o mt19937ar.o vector.o path.o random.o misc.o test-charge: test-charge.o sno_charge.o misc.o @@ -2,6 +2,30 @@ #include <math.h> #include <stdlib.h> /* for size_t */ +int isclose(double a, double b, double rel_tol, double abs_tol) +{ + /* Returns 1 if a and b are "close". This algorithm is taken from Python's + * math.isclose() function. + * + * See https://www.python.org/dev/peps/pep-0485/. */ + return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol); +} + +int allclose(double *a, double *b, size_t n, double rel_tol, double abs_tol) +{ + /* Returns 1 if all the elements of a and b are "close". This algorithm is + * taken from Python's math.isclose() function. + * + * See https://www.python.org/dev/peps/pep-0485/. */ + size_t i; + + for (i = 0; i < n; i++) { + if (!isclose(a[i],b[i],rel_tol,abs_tol)) return 0; + } + + return 1; +} + double logsumexp(double *a, size_t n) { /* Returns the log of the sum of the exponentials of the array `a`. @@ -3,6 +3,8 @@ #include <stdlib.h> /* for size_t */ +int isclose(double a, double b, double rel_tol, double abs_tol); +int allclose(double *a, double *b, size_t n, double rel_tol, double abs_tol); double logsumexp(double *a, size_t n); double norm(double x, double mu, double sigma); double norm_cdf(double x, double mu, double sigma); @@ -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); @@ -108,15 +108,6 @@ struct solid_angle_results { {2.0,2.0,0.282707} }; -int isclose(double a, double b, double rel_tol, double abs_tol) -{ - /* Returns 1 if a and b are "close". This algorithm is taken from Python's - * math.isclose() function. - * - * See https://www.python.org/dev/peps/pep-0485/. */ - return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol); -} - int test_muon_get_T(char *err) { /* A very simple test to make sure the energy as a function of distance |