Commit 8c3f7e60 authored by Bruce Allen's avatar Bruce Allen
Browse files

Preliminary code to calculate waveforms. HAS A MISTAKE THAT I NEED TO FIX

parent 15aaa3f0
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
// 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
// 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″ 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
// 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.
const double galaxy[2] = {
-1.0*deg_to_rad*(23.0 + 23.0/60.0 + 4.0/3600.0),
deg_to_rad*41.092981
};
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];
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];
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
double u[3];
double v[3];
// Angle defining line of source ellipse
double orientation;
};
struct Source source;
// 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);
// 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) {
// lattitude and longitude
double lat = src->location[0];
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;
// 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];
}
}
// 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 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];
}
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;
// Initially set to zero, fix later
source.orientation = 0.0;
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*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*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*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\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]
);
}
void print_source() {
printf("source at:\n"
"lattitude: %f\n"
"longitude: %f\n"
"earth center to source %f %f %f\n\n",
source.location[0],
source.location[1],
source.vec[0], source.vec[1],source.vec[2]
);
}
// This dots the mass quadrupole with the antenna functions
void get_UUUVVV(int det, double *uu, double *uv, double *vv) {
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]);
*uu = UU;
*uv = UV;
*vv = VV;
}
int main(int argc, char *argv[]) {
double iota=atof(argv[1]);
printf("Iota = %f degrees\n", iota);
iota *= deg_to_rad;
double ci=cos(iota);
double si=sin(iota);
// 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;
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);
}
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