Commit c3acd36e authored by Oliver Bock's avatar Oliver Bock
Browse files

Merge branch 'antenna-patterns' of gitlab:einsteinathome/pulsatingscience into antenna-patterns

parents afa95f7c 24745a7c
antenna_main: antenna_lib.o antenna_main.c antenna_lib.h
cc -o antenna_main antenna_main.c antenna_lib.o -lm
antenna.o: antenna_lib.c antenna_lib.h
cc -c antenna_lib.c
clean:
rm -f antenna_lib.o antenna_main
// Copyright Bruce Allen 2017
// Compile with gcc -o antenna antenna.c -lm
// REMAINING THINGS TO CHECK:
// (a) sign conventions for h, which arm is positive
// (b) direction and origin conventions for the polarization axis
// (c) double check the hand calculations
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
// Structure for passing information about the GW source and detectors
struct InputStruct {
// orbital orientation in degrees, 0 to 180. Zero degrees has
// orbital angular momentum pointing to the earth
double iota;
// orientation of long axis of ellipse, 0 to 360, in degrees CCW
// from North
double psi;
// orientation of detector arms away from actual, 0 to 360, in
// degrees CCW when viewed from directly overhead
// Ordering is LLO, LHO, VIRGO
double orientation[3];
};
// Structure for returning information about the response. The array
// elements are in the order LLO, LHO, VIRGO
struct OutputStruct {
// amplitude of sin function, dimensionless
double amp[3];
// phase of sin function, in degrees from 0 to 360
double phase[3];
// time delays, in milliseconds
double dt[3];
};
// Global variables for passing information. Ugly style but doesn't
// matter for this!
struct InputStruct inputdata;
// 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
// milliseconds. Get this by dividing radius in km by speed of light in
// km/s.
const double radius_earth = 6371.0*1000.0/299792.458;
// Event time
// GPS 1187008882
// UTC: Thu 2017 Aug 17 12:41:04
......@@ -171,24 +211,24 @@ void populate_detectors(){
// LLO
// Sanity checked using Google Earth!
strcpy(detectors[0].name, "LLO");
strcpy(detectors[0].name, " LLO ");
detectors[0].location[0] = 30.56*deg_to_rad;
detectors[0].location[1] = -90.77*deg_to_rad;
detectors[0].orientation = 198.0*deg_to_rad;
detectors[0].orientation = (198.0+inputdata.orientation[0])*deg_to_rad;
// LHO
// Sanity checked using Google Earth!
strcpy(detectors[1].name, "LHO");
strcpy(detectors[1].name, " LHO ");
detectors[1].location[0] = 46.45*deg_to_rad;
detectors[1].location[1] = -119.41*deg_to_rad;
detectors[1].orientation = 126.8*deg_to_rad;
detectors[1].orientation = (126.8+inputdata.orientation[1])*deg_to_rad;
// VIRGO
// Sanity checked using Google Earth!
strcpy(detectors[2].name, "Virgo");
strcpy(detectors[2].name, "VIRGO");
detectors[2].location[0] = 43.63*deg_to_rad;
detectors[2].location[1] = 10.5*deg_to_rad;
detectors[2].orientation = 71.5*deg_to_rad;
detectors[2].orientation = (71.5+inputdata.orientation[2])*deg_to_rad;
int i;
......@@ -235,73 +275,148 @@ 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++) 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++)
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
// line, iota and psi, angles in degrees.
int main(int argc, char *argv[]) {
// to loop over detectors
// library function that can be called either from the GUI code or
// from a stand-alone terminal program. This function does not modify
// the input struct, but does populate/modify the output strut!
void get_antenna(struct OutputStruct *out, struct InputStruct *in) {
int i;
// check syntax crudely, issue usage message
if (argc != 3) {
fprintf(stderr,
"Wrong argument count! Correct usage:\n"
"%s float_iota float_psi\n",
argv[0]);
exit(1);
}
// setup and test
// print_galaxy_coordinates();
// set up the source and detectors
inputdata=*in;
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]);
double psi = atof(argv[2]);
printf("Iota = %f degrees\nPsi = %f degrees\n", iota, psi);
double iota = inputdata.iota;
double psi = inputdata.psi;
#ifdef DEBUG
fprintf(stderror, "Iota = %f degrees\nPsi = %f degrees\n", iota, psi);
#endif
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);
// compute time delay in milliseconds
dt *= radius_earth;
// printf("For detector %s the waveform is w^2 [ %.3f sin(2w(t%+.1f ms))%+.3f cos(2w(t%+.1f ms)) ]\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)
printf("For detector %s the waveform is w^2 [ %f sin(2wt) + %f cos(2wt) ] \n", detectors[i].name, X, Y);
double amp = sqrt(X*X + Y*Y);
double ang = 180.0*atan2(Y, X)/M_PI;
// pass outputs
out->amp[i]= amp;
out->phase[i] = ang;
out->dt[i] = dt;
#ifdef DEBUG
// degree character in UTF-8 character set (most likely terminal type!)
int deg1=0xC2, deg2=0xB0;
fprintf(stderr, "For detector %s the waveform is %.3f w^2 sin(2w[t%+.1f ms]%+.1f%c%c)\n", detectors[i].name, amp, dt, ang, deg1, deg2);
#endif
}
return;
}
// 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
int i;
// to pass data in and out
struct InputStruct myinput;
struct OutputStruct myoutput;
// check syntax crudely, issue usage message
if (argc != 3) {
fprintf(stderr,
"Wrong argument count! Correct usage:\n"
"%s float_iota_in_degrees float_psi_in_degrees\n",
argv[0]);
exit(1);
}
// pass inclination angle, polarization axis, orientation offsets
myinput.iota = atof(argv[1]);
myinput.psi = atof(argv[2]);
for (i=0;i<3;i++) myinput.orientation[i]=0.0;
printf("Iota = %f degrees\nPsi = %f degrees\n", myinput.iota, myinput.psi);
// now compute responses
get_antenna(&myoutput, &myinput);
// degree character in UTF-8 character set (most likely terminal type!)
int deg1=0xC2, deg2=0xB0;
// Now display waveforms
for (i=0; i<3; i++)
printf("For detector %s the waveform is %.3f w^2 sin(2w[t%+.1f ms]%+.1f%c%c)\n", detectors[i].name, myoutput.amp[i], myoutput.dt[i], myoutput.phase[i], deg1, deg2);
return 0;
}
// Copyright Bruce Allen 2017
// Compile with gcc -o antenna antenna.c -lm
// REMAINING THINGS TO CHECK:
// (a) sign conventions for h, which arm is positive
// (b) direction and origin conventions for the polarization axis
// (c) double check the hand calculations
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include "antenna_lib.h"
// Global variables for passing information. Ugly style but doesn't
// matter for this!
struct InputStruct inputdata;
// 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
// milliseconds. Get this by dividing radius in km by speed of light in
// km/s.
const double radius_earth = 6371.0*1000.0/299792.458;
// Event time
// GPS 1187008882
// UTC: Thu 2017 Aug 17 12:41:04
// Julian Day 2457983.02852
// Modified Julian Day 57982.52852
// The convention in this code is that Lattitude is measured north of
// the equator. So Lattitude -10 means 10 degrees SOUTH of the
// 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″
// 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
double JD = 2457983.02852;
// This is the Julian day (note previous integer part!)
double JD0 = 2457982.5;
// Compute Julian century t
double day = JD0-2451545.0;
double t = day/36525.0;
// The following is equation 11.4 of the previous reference
double degrees = 280.46061837 + 360.98564736629*(JD-2451545.0)+t*t*0.000387933 - t*t*t/38710000.0;
degrees = fmod(degrees,360.0);
printf("Greenwich Mean Sidereal time is %f degrees\n", degrees);
// Now the longitude of our source (with sign conventions above) is given by L = RA - MST, where
// MST (in degrees) is given above and RA is the RA of our source
// Here a positive number means "East of Greenwich"
double longitude = (13*3600.0+9*60+47.2)/240 - degrees;
printf("Longitude of source is %f degrees\n", longitude);
}
// Note these are in the order lattitude,longitude. See routine
// "print_galaxy_coordinates" to convert location on sky to Earth
// 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];
// 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];
// lattitude north of equator, radians
// longitude east from Greenwich, radians
double location[2];
// orientation of Y arm CCW from North, radians
double orientation;
// 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];
// Unit vectors along the two arms
double lx[3];
double ly[3];
};
// LLO, LHO, Virgo, in that order
struct Detector detectors[3];
// input is a lattitude/longitude set in radians
// output is unit vectors
void make_unit_vectors(double *out, double *in) {
double lat = in[0];
double lon = in[1];
out[0] = cos(lon)*cos(lat);
out[1] = sin(lon)*cos(lat);
out[2] = sin(lat);
}
void make_u_v_vectors(struct Source *src) {
// lattitude and longitude
double lat = src->location[0];
double lon = src->location[1];
// construct unit vectors perpendicular to line of sight
// 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);
}
// 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) {
double lat = det->location[0];
double lon = det->location[1];
det->north[0] = -sin(lat)*cos(lon);
det->north[1] = -sin(lat)*sin(lon);
det->north[2] = cos(lat);
det->east[0] = -sin(lon);
det->east[1] = cos(lon);
det->east[2] = 0.0;
double cpsi = cos(det->orientation);
double spsi = sin(det->orientation);
int i;
for (i=0; i<3; i++) {
det->lx[i] = cpsi*det->east[i] + spsi*det->north[i];
det->ly[i] = -spsi*det->east[i] + cpsi*det->north[i];
}
}
// 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;
make_unit_vectors(source.vec,source.location);
make_u_v_vectors(&source);
}
// initially take detector locations and orientations from
// https://arxiv.org/pdf/gr-qc/9607075.pdf
// Be sure to check this later!!
void populate_detectors(){
// LLO
// Sanity checked using Google Earth!
strcpy(detectors[0].name, " LLO ");
detectors[0].location[0] = 30.56*deg_to_rad;
detectors[0].location[1] = -90.77*deg_to_rad;
detectors[0].orientation = (198.0+inputdata.orientation[0])*deg_to_rad;
// LHO
// Sanity checked using Google Earth!
strcpy(detectors[1].name, " LHO ");
detectors[1].location[0] = 46.45*deg_to_rad;
detectors[1].location[1] = -119.41*deg_to_rad;
detectors[1].orientation = (126.8+inputdata.orientation[1])*deg_to_rad;
// VIRGO
// Sanity checked using Google Earth!
strcpy(detectors[2].name, "VIRGO");
detectors[2].location[0] = 43.63*deg_to_rad;
detectors[2].location[1] = 10.5*deg_to_rad;
detectors[2].orientation = (71.5+inputdata.orientation[2])*deg_to_rad;
int i;
// Coordinate system has x/y plane through the equator, north pole
// along positive z axis, and Greenwich passing through x axis (ie
// y=0).
for (i=0;i<3;i++) {
make_unit_vectors(detectors[i].vec, detectors[i].location);
make_north_east_vectors(detectors+i);
}
}
void print_detector(int det) {
printf("name: %s\n"
"lattitude (north): %f\n"
"longitude (east): %f\n"
"orientation: (CCW from North): %f\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].lx[0], detectors[det].lx[1], detectors[det].lx[2],
detectors[det].ly[0], detectors[det].ly[1], detectors[det].ly[2]
);
}
void print_source() {
printf("source at:\n"
"lattitude: %f\n"
"longitude: %f\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.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_UV_combinations(int det, double *alpha, double *beta, double *rdotn) {
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]);
// 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;
}
// library function that can be called either from the GUI code or
// from a stand-alone terminal program. This function does not modify
// the input struct, but does populate/modify the output strut!
void get_antenna(struct OutputStruct *out, struct InputStruct *in) {
int i;
// set up the source and detectors
inputdata=*in;
populate_source();
populate_detectors();
double iota = inputdata.iota;
double psi = inputdata.psi;
#ifdef DEBUG
fprintf(stderror, "Iota = %f degrees\nPsi = %f degrees\n", iota, psi);
#endif
iota *= deg_to_rad;
psi *= deg_to_rad;
// 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, dt;
// Ugly, should pass as argument to populate_detectors()
strcpy(out->name[i], detectors[i].name);
// get antenna pattern overlap with source plane
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);
// compute time delay in milliseconds
dt *= radius_earth;
// printf("For detector %s the waveform is w^2 [ %.3f sin(2w(t%+.1f ms))%+.3f cos(2w(t%+.1f ms)) ]\n", detectors[i].name, X, dt, Y, dt);