diff options
| -rw-r--r-- | calculate_limits.c | 67 |
1 files changed, 67 insertions, 0 deletions
diff --git a/calculate_limits.c b/calculate_limits.c index 8815324..47ed1c9 100644 --- a/calculate_limits.c +++ b/calculate_limits.c @@ -17,6 +17,10 @@ double sno_depth = 2072640; /* Fiducial volume in mm. */ double radius_fiducial = 5000; +/* Cartesian coordinates of SNO in earth frame. They need to be global since + * they are used in some functions. */ +double x_sno[3]; + double deg2rad(double deg) { return deg*M_PI/180.0; @@ -47,6 +51,67 @@ void cartesian2sphere(double x, double y, double z, double *r, double *theta, do *phi = atan2(y,x); } +void cross(double *a, double *b, double *c) +{ + c[0] = a[1]*b[2] - a[2]*b[1]; + c[1] = a[2]*b[0] - a[0]*b[2]; + c[2] = a[0]*b[1] - a[1]*b[0]; +} + +double dot(double *a, double *b) +{ + return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; +} + +double norm(double *a) +{ + return sqrt(dot(a,a)); +} + +void normalize(double *a) +{ + double n = norm(a); + a[0] /= n; + a[1] /= n; + a[2] /= n; +} + +/* Rotate a vector x around the vector dir by an angle theta. */ +void rotate(double *result, double *x, double *dir, double theta) +{ + double a = dot(dir,x); + double b[3]; + + double sin_theta = sin(theta); + double cos_theta = cos(theta); + + /* Make sure the direction vector is normalized. */ + normalize(dir); + + cross(x,dir,b); + + result[0] = x[0]*cos_theta + dir[0]*a*(1-cos_theta) + b[0]*sin_theta; + result[1] = x[1]*cos_theta + dir[1]*a*(1-cos_theta) + b[1]*sin_theta; + result[2] = x[2]*cos_theta + dir[2]*a*(1-cos_theta) + b[2]*sin_theta; +} + +/* Rotate a vector in earth centered coordinates to SNO coordinates (doesn't do + * the translation). */ +void rotate_earth_to_sno(double *x_earth, double *x_sno) +{ + double dir[3]; + double z[3] = {0,0,1}; + + cross(x_sno, z, dir); + + /* Normalize. */ + normalize(dir); + + double theta = acos(dot(x_sno,z)/norm(x_sno)); + + rotate(x_sno, x_earth, dir, theta); +} + int main(int argc, char **argv) { /* Spherical angles for the SNO detector in the earth frame which has z @@ -55,5 +120,7 @@ int main(int argc, char **argv) double sno_theta = deg2rad(latitude + 90.0); double sno_phi = deg2rad(longitude); + sphere2cartesian(radius_earth - sno_depth, sno_theta, sno_phi, x_sno, x_sno+1, x_sno+2); + return 0; } |