Commit 5178974e authored by Bruce Allen's avatar Bruce Allen
Browse files

Cleaned up the code a bit, fixing comments and removing some scaffolding

parent 9a03e45e
......@@ -20,15 +20,16 @@ const double deg_to_rad = M_PI/180.0;
// equator. Longitude is measured East of Greenwich. So longitude
// -90 means 90 degrees West of Greenwich (somewhere in the USA).
// Galaxy NGC 4993 (from Wikipedia) Right ascension 13h 09m 47.2s
// Declination −23° 23′ 4″ DEC and RA in radians. We will need to
// convert this to location over the Earth at the event time. It will
// turn out to be Lattitude −23° 23′ 4″ and Longitude 41.092981
// degrees
// Galaxy NGC 4993 (from Wikipedia):
// Right ascension: 13h 09m 47.2s
// Declination: −23° 23′ 4″
// We will need to convert this to location over the Earth at the
// event time. It will turn out to be:
// Lattitude −23.3844 degrees (south of equator)
// Longitude 41.092981 degrees (east of Greenwich)
// From Chapter 11 of "Astronomical Algorithms" by Jean Meeus, published 1991
void print_galaxy_coordinates() {
// Here is the Julian day of the event, including a fractional part
......@@ -56,13 +57,25 @@ void print_galaxy_coordinates() {
// Note these are in the order lattitude,longitude. See routine
// "print_galaxy_coordinates" to convert location on sky to Earth
// location at event time.
// location at event time. These are set/used in the function
// populate_source() below.
struct Source {
// lattitude north of equator, radians
// longitude east from Greenwich, radians
double location[2];
// Unit vector from Earth to source
double vec[3];
const double galaxy[2] = {
-1.0*deg_to_rad*(23.0 + 23.0/60.0 + 4.0/3600.0),
deg_to_rad*41.092981
// Unit vectors defining plane perpendicular to the source
// direction. U points east, V points north
double u[3];
double v[3];
};
struct Source source;
struct Detector {
// null terminated char string
char name[8];
......@@ -76,6 +89,9 @@ struct Detector {
// Unit vector from center of Earth to detector
double vec[3];
// Unit vectors pointing along the north and east directions at the
// detector site
double north[3];
double east[3];
......@@ -87,27 +103,6 @@ struct Detector {
// LLO, LHO, Virgo, in that order
struct Detector detectors[3];
struct Source {
// lattitude north of equator, radians
// longitude east from Greenwich, radians
double location[2];
// Unit vector from Earth to source
double vec[3];
// Unit vectors defining plane perpendicular to the source
// direction. U points
double u[3];
double v[3];
#if 0
// Angle defining line of source ellipse
double orientation;
#endif
};
struct Source source;
// input is a lattitude/longitude set in radians
// output is unit vectors
void make_unit_vectors(double *out, double *in) {
......@@ -116,13 +111,6 @@ void make_unit_vectors(double *out, double *in) {
out[0] = cos(lon)*cos(lat);
out[1] = sin(lon)*cos(lat);
out[2] = sin(lat);
// sanity check
double len = out[0]*out[0]+out[1]*out[1]+out[2]*out[2];
if (fabs(len-1.0)>0.0001) {
fprintf(stderr,"Squares check failed!\n");
exit(1);
}
}
void make_u_v_vectors(struct Source *src) {
......@@ -140,18 +128,6 @@ void make_u_v_vectors(struct Source *src) {
src->v[0] = -sin(lat)*cos(lon);
src->v[1] = -sin(lat)*sin(lon);
src->v[2] = cos(lat);
// Have incorporated orientation angle psi into main code body!
#if 0
// now make a rotated basis, rotated by oriention angle
double cpsi = cos(src->orientation);
double spsi = sin(src->orientation);
int i;
for (i=0; i<3; i++) {
src->u[i] = cpsi*u[i] + spsi*v[i];
src->v[i] = -spsi*u[i] + cpsi*v[i];
}
#endif
}
// input is a lattitude/longitude set in radians
......@@ -177,23 +153,12 @@ void make_north_east_vectors(struct Detector *det) {
}
}
#if 0
// sets out = a X b
void cross_prod(double *out, double *a, double *b) {
out[0]=a[1]*b[2]-a[2]*b[1];
out[1]=a[2]*b[0]-a[0]*b[2];
out[2]=a[0]*b[1]-a[1]*b[0];
}
#endif
// sets up the different vectors needed to define the source
void populate_source() {
source.location[0]= -1.0*deg_to_rad*(23.0 + 23.0/60.0 + 4.0/3600.0);
source.location[1]=deg_to_rad*41.092981;
source.location[0] = -1.0*deg_to_rad*(23.0 + 23.0/60.0 + 4.0/3600.0);
source.location[1] = deg_to_rad*41.092981;
#if 0
// Initially set to zero, fix later
source.orientation = 0.0;
#endif
make_unit_vectors(source.vec,source.location);
make_u_v_vectors(&source);
}
......@@ -279,15 +244,16 @@ void get_UV_combinations(int det, double *alpha, double *beta) {
double *sv=source.v;
double *dx=detectors[det].lx;
double *dy=detectors[det].ly;
for (i=0;i<3;i++) for (j=0;j<3;j++) a+=(su[i]*su[j]-sv[i]*sv[j])*(dx[i]*dx[j]-dy[i]*dy[j]);
for (i=0;i<3;i++) for (j=0;j<3;j++) b+=(su[i]*sv[j]+sv[i]*su[j])*(dx[i]*dx[j]-dy[i]*dy[j]);
*alpha = a;
*beta = b;
*alpha = a;
*beta = b;
}
// This function requires two floating point arguments on the command
// line, iota and psi, angles in degrees.
int main(int argc, char *argv[]) {
// to loop over detectors
......@@ -311,21 +277,22 @@ int main(int argc, char *argv[]) {
// for (i=0; i<3; i++) print_detector(i);
// get inclination angle, polarization axis
double iota=atof(argv[1]);
double psi=atof(argv[2]);
double iota = atof(argv[1]);
double psi = atof(argv[2]);
printf("Iota = %f degrees\nPsi = %f degrees\n", iota, psi);
iota *= deg_to_rad;
psi *= deg_to_rad;
double ci=cos(iota);
double si=sin(iota);
double c2p=cos(2*psi);
double s2p=sin(2*psi);
// get angles needed to calculate antenna response
double ci = cos(iota);
double si = sin(iota);
double c2p = cos(2*psi);
double s2p = sin(2*psi);
// Now find waveforms. At each site, waveform is w^2 [X sin(2wt) + Y cos(2wt) ]
// loop over detectors
for (i=0; i<3; i++) {
double alpha,beta,X,Y;
double alpha, beta, X, Y;
// get antenna pattern overlap with source plane
get_UV_combinations(i, &alpha, &beta);
......@@ -334,7 +301,7 @@ int main(int argc, char *argv[]) {
X = 2.0*ci*(alpha*s2p - beta*c2p);
Y = -(ci*ci + 1.0)*(alpha*c2p + beta*s2p);
printf("For detector %s the waveform is w^2 [ %f sin(2wt) + %f cos(2wt)] \n", detectors[i].name, X, Y);
printf("For detector %s the waveform is w^2 [ %f sin(2wt) + %f cos(2wt) ] \n", detectors[i].name, X, Y);
}
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment