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

Corrected sign of iota to correspond to findchirp paper

Included effects of arrival time delay in the formulae for the waveforms (assuming a round earth)

Expressed output X sin(phi) + Y cos(phi) in the alternative form r sin(phi+offset)

Currently offset is printed in degrees.  Printing code assumes UTF-8 on the terminal. This may be fragile, I'm not sure.
parent 39647a64
......@@ -9,6 +9,11 @@
// For converting degrees to radians
const double deg_to_rad = M_PI/180.0;
// For computing arrival time delays, need mean radius of the earth in
// seconds. Get this by dividing radius in km by speed of light in
// km/s.
const double radius_earth = 6371.0/299792.458;
// Event time
// GPS 1187008882
// UTC: Thu 2017 Aug 17 12:41:04
......@@ -235,21 +240,41 @@ void print_source() {
}
// This dots the mass quadrupole with the antenna functions
void get_UV_combinations(int det, double *alpha, double *beta) {
void get_UV_combinations(int det, double *alpha, double *beta, double *rdotn) {
double a=0,b=0;
double a=0,b=0,dot=0;
int i, j;
double *su=source.u;
double *sv=source.v;
double *n=source.vec;
double *dx=detectors[det].lx;
double *dy=detectors[det].ly;
double *r=detectors[det].vec;
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]);
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]);
// The quantity dot is positive if the vector from the earth center
// to a detector has the SAME direction as the vector from the
// earth center to the source. This means that the signal arrives
// EARLIER at that detector. Hence the signal waveform seen at the
// detector is of the form h(t) = wave(t + R*dot) where R is the
// (positive) earth radius in units of time and wave(t) is the
// emitted source waveform.
for (i=0; i<3; i++)
dot += n[i]*r[i];
*alpha = a;
*beta = b;
*rdotn = dot;
}
// This function requires two floating point arguments on the command
......@@ -263,7 +288,7 @@ int main(int argc, char *argv[]) {
if (argc != 3) {
fprintf(stderr,
"Wrong argument count! Correct usage:\n"
"%s float_iota float_psi\n",
"%s float_iota_in_degrees float_psi_in_degrees\n",
argv[0]);
exit(1);
}
......@@ -283,25 +308,48 @@ int main(int argc, char *argv[]) {
iota *= deg_to_rad;
psi *= deg_to_rad;
// get angles needed to calculate antenna response
double ci = cos(iota);
double si = sin(iota);
// get angles needed to calculate antenna response. The minus sign
// in the definition of cos(iota) is because my calculational notes
// assume that iota=0 corresponds to orbital motion CCW in the
// (right-handed) UV coordinate system. That corresponds to having
// the orbital angular momentum pointing radially AWAY from earth.
// But the standard convention is that iota=0 is face-on, meaning
// orbital angular momentum pointing TO the earth. So I change the
// sign of cos(iota) below to correct for this.
double ci = -cos(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, dt;
// get antenna pattern overlap with source plane
get_UV_combinations(i, &alpha, &beta);
get_UV_combinations(i, &alpha, &beta, &dt);
// 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);
// compute time delay
dt *= radius_earth;
printf("For detector %s the waveform is w^2 [ %.3f sin(2w(t%+.3f s))%+.3f cos(2w(t%+.3f s)) ]\n", detectors[i].name, X, dt, Y, dt);
// compute alternative form of output
// X sin(phi) + Y cos(phi) = sqrt(X^2+Y^2) sin(phi + ang) where ang=atan2(Y,X)
double amp = sqrt(X*X + Y*Y);
double ang = atan2(Y, X);
// degree character in UTF-8 character set (most likely terminal type!)
int deg1=0xC2, deg2=0xB0;
printf(" = %.3f w^2 cos[2w(t%+.3f s) %+.1f%c%c]\n", amp, dt, ang, deg1, deg2);
}
return 0;
}
......@@ -8,16 +8,19 @@ x,y,z coordinate system for Earth is as follows:
(a) The equator lies in the XY plane
(b) The north pole lies on the positive Z axis
(c) The line of longitude 0 lies in the XZ plane and passes through
(c) The line of longitude 0 lies in the XZ plane, and passes through
the positive X axis
#THE FOLLOWING STATEMENT MAY NEED A SIGN CHANGE
The orbital inclination angle iota (from 0 to 180 degrees) is the angle
between the orbital angular momentum vector (defined with the
right-hand rule) and the vector from earth to source.
According to the findchirp paper, the orbital inclination angle iota
(from 0 to 180 degrees) is the angle between the orbital angular
momentum vector (defined with the right-hand rule) and the vector from
source to earth.
Face-on is iota=0 and face-off is iota=180.
Face-on is iota=0 (orbital angular momentum pointing to earth) and
face-off is iota=180 (orbital angular momentum pointing away from
earth).
## CHECK THIS!
The orientation angle psi is the angle (from 0 to 180 degrees,
measured counterclockwise from North) between the long angle of the
ellipse (projected into the plane of the sky) and the northerly
......
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