Commit fa2da70f authored by Bruce Allen's avatar Bruce Allen
Browse files

I have now extracted the antenna dependence on iota and psi explictly.

I am calculating the second derivative of the mass quadrupole in the 2-dimensional plane
orthogonal to the line of site.  The mass quadrupole is traceless w.r.t. to the two-dimensional
metric in that plane
parent d1759ca6
......@@ -6,7 +6,6 @@
// For converting degrees to radians
const double deg_to_rad = M_PI/180.0;
// Event time
// GPS 1187008882
// UTC: Thu 2017 Aug 17 12:41:04
......@@ -82,8 +81,6 @@ struct Detector {
double ly[3];
};
// LLO, LHO, Virgo, in that order
struct Detector detectors[3];
......@@ -96,15 +93,14 @@ struct Source {
double vec[3];
// Unit vectors defining plane perpendicular to the source
// direction
// direction. U points
double u[3];
double v[3];
#if 0
// Angle defining line of source ellipse
double orientation;
#endif
};
struct Source source;
......@@ -132,14 +128,18 @@ void make_u_v_vectors(struct Source *src) {
double lon = src->location[1];
// construct unit vectors perpendicular to line of sight
double u[3],v[3];
u[0] = -sin(lat)*cos(lon);
u[1] = -sin(lat)*sin(lon);
u[2] = cos(lat);
v[0] = -sin(lon);
v[1] = cos(lon);
v[2] = 0.0;
// u points east
src->u[0] = -sin(lon);
src->u[1] = cos(lon);
src->u[2] = 0.0;
// v points north, so u,v are a right-handed pair like x,y
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);
......@@ -148,10 +148,9 @@ void make_u_v_vectors(struct Source *src) {
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
// output is unit vectors in the north and east directions
void make_north_east_vectors(struct Detector *det) {
......@@ -175,25 +174,25 @@ 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
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;
#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);
//
}
// initially take detector locations and orientations from
......@@ -239,12 +238,16 @@ void print_detector(int det) {
"lattitude (north): %f\n"
"longitude (east): %f\n"
"orientation: (CCW from North): %f\n"
"earth center to detector %f %f %f\n\n",
"earth center to detector %f %f %f\n"
"vector along X arm %f %f %f\n"
"vector along Y arm %f %f %f\n\n",
detectors[det].name,
detectors[det].location[0],
detectors[det].location[1],
detectors[det].orientation,
detectors[det].vec[0], detectors[det].vec[1], detectors[det].vec[2]
detectors[det].vec[0], detectors[det].vec[1], detectors[det].vec[2],
detectors[det].lx[0], detectors[det].lx[1], detectors[det].lx[2],
detectors[det].ly[0], detectors[det].ly[1], detectors[det].ly[2]
);
}
......@@ -252,57 +255,74 @@ void print_source() {
printf("source at:\n"
"lattitude: %f\n"
"longitude: %f\n"
"earth center to source %f %f %f\n\n",
"earth center to source: %f %f %f\n"
"U vector: %f %f %f\n"
"V vector: %f %f %f\n\n",
source.location[0],
source.location[1],
source.vec[0], source.vec[1],source.vec[2]
source.vec[0], source.vec[1],source.vec[2],
source.u[0],source.u[1],source.u[2],
source.v[0],source.v[1],source.v[2]
);
}
// This dots the mass quadrupole with the antenna functions
void get_UUUVVV(int det, double *uu, double *uv, double *vv) {
void get_UV_combinations(int det, double *alpha, double *beta) {
double a=0,b=0;
int i, j;
double UU=0,VV=0,UV=0;
int i, j;
for (i=0;i<3;i++) for (j=0;j<3;j++) UU+=source.u[i]*source.u[j]*(detectors[det].lx[i]*detectors[det].lx[j]-detectors[det].ly[i]*detectors[det].ly[j]);
for (i=0;i<3;i++) for (j=0;j<3;j++) UV+=(source.u[i]*source.v[j]+source.v[i]*source.u[j])*(detectors[det].lx[i]*detectors[det].lx[j]-detectors[det].ly[i]*detectors[det].ly[j]);
for (i=0;i<3;i++) for (j=0;j<3;j++) VV+=source.v[i]*source.v[j]*(detectors[det].lx[i]*detectors[det].lx[j]-detectors[det].ly[i]*detectors[det].ly[j]);
double *su=source.u;
double *sv=source.v;
double *dx=detectors[det].lx;
double *dy=detectors[det].ly;
*uu = UU;
*uv = UV;
*vv = VV;
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;
}
int main(int argc, char *argv[]) {
// to loop over detectors
int i;
// setup and test
// print_galaxy_coordinates();
populate_source();
// print_source();
populate_detectors();
// for (i=0; i<3; i++) print_detector(i);
// get inclination angle, polarization axis
double iota=atof(argv[1]);
printf("Iota = %f degrees\n", iota);
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);
// setup and test
// print_galaxy_coordinates();
populate_detectors();
populate_source();
// for (i=0; i<3; i++) print_detector(i);
// print_source();
// Now find waveforms. At each site, waveform is A sin^2(wt) + B cos^2(wt) + C sin(wt)cos(wt)
double UU, UV, VV;
int i;
// 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 A=0,B=0,C=0;
get_UUUVVV(i, &UU, &UV, &VV);
A = ci*ci*(2.0*VV/3.0-UU/3.0);
B = 2.0*(UU+VV)/3.0;
C = ci*UV;
printf("For detector %s the waveform is %f sin^2(wt) + %f cos^2(wt) + %f sin(wt)cos(wt)\n", detectors[i].name, A, B, C);
double alpha,beta,X,Y;
// get antenna pattern overlap with source plane
get_UV_combinations(i, &alpha, &beta);
// form combinations as functions of inclination and polarization angles
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);
}
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