antenna_lib.cpp

// 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 <string.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;
   
#if 0
   // for testing purposes, put source directly over a detector
   source.location[0] = 46.45*deg_to_rad;
   source.location[1] = -119.41*deg_to_rad;
#endif

   
   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);

      // 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 = 180.0*atan2(Y, X)/M_PI;
      
      // pass outputs In this code the relevant phases are the
      // DIFFERENCES between those at the different detectors.  So to
      // make for a less confusing visual display in the GUI, I have
      // offset the output phase as a function of psi, in a way that
      // tries to reduce confusion about the meaining of the orbital
      // coalescence phase.
      out->amp[i]= amp;
      out->phase[i] = ang-2*psi*ci/deg_to_rad;
      out->dt[i] = dt;

      fprintf(stderr, "ang = %f psi = %f  ci =% f\n", ang, psi/deg_to_rad, ci);
      
#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;
}