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 when viewed from
   // infinity.  So U x V points outwards from Earth center.
   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.  Note that East x North points outwards from
   // center of Earth
   double east[3];
   double north[3];
   

   // Unit vectors along the two arms
   double lx[3];
   double ly[3];

   // response tensor with same conventions as LAL
   // meaning 1/2(L_x L_x - L_y L_y)
   double response[3][3];
};

// LLO, LHO, Virgo, KAGRA, LIGO-Infia in that order
struct Detector detectors[5];

// 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 (viewed from infinity)
   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];

   // vector pointing Eastwards (viewed from infinity)
   det->east[0] = -sin(lon);
   det->east[1] = cos(lon);
   det->east[2] = 0.0;

   // orthogonal vector which projects onto Northwards direction
   // (viewed from infinity)
   det->north[0] = -sin(lat)*cos(lon);
   det->north[1] = -sin(lat)*sin(lon);
   det->north[2] = cos(lat);

   // rotate these through offset angle if desired. Increasing
   // orientation angle is CCW viewed from infinity
   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];
   }
}

// makes a detector response tensor with the same normalization
// conventions as LAL, meaning 1/2(Lx Lx - Ly Ly)
void make_detector_response(struct Detector *det) {
   int i,j;

   for (i=0;i<3; i++)
      for (j=0;j<3;j++)
	 det->response[i][j] =
	    0.5*(det->lx[i]*det->lx[j]-det->ly[i]*det->ly[j]);
   return;
}

// sets up the different vectors needed to define the source
void populate_source() {
   // this is the GW170817 location
   //    Lattitude −23.3844 degrees (south of equator)
   //    Longitude 41.092981 degrees (east of Greenwich)
   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 LHO
   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;

#if DEBUG
   fprintf(stderr, "populate_detectors(): LHO orientation %f\n", inputdata.orientation[1]);
#endif
   
   // 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;

   // KAGRA
   strcpy(detectors[3].name, "KAGRA");
   detectors[3].location[0] = 77.76*deg_to_rad;
   detectors[3].location[1] = 10.5*deg_to_rad;
   detectors[3].orientation = (25+inputdata.orientation[3])*deg_to_rad;

   // LIGO-India
   // Based on fig. 7 in http://iopscience.iop.org/article/10.1088/0264-9381/30/15/155004#cqg460254s5
   strcpy(detectors[4].name, " LI  ");
   detectors[4].location[0] = 77.76*deg_to_rad;
   detectors[4].location[1] = 10.02*deg_to_rad;
   detectors[4].orientation = ((90 + 58.2)+inputdata.orientation[4])*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<5;i++) {
      make_unit_vectors(detectors[i].vec, detectors[i].location);
      make_north_east_vectors(detectors+i);
      make_detector_response(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]
	  );

   int i;
   printf("Response tensor:\n");
   for (i=0;i<3; i++)
      printf("% 0.4f % 0.4f % 0.4f\n",
	     detectors[det].response[i][0],
	     detectors[det].response[i][1],
	     detectors[det].response[i][2]
	     );
   printf("\n");
   return;
}

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 *u=source.u;
   double *v=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 += (u[i]*u[j] - v[i]*v[j]) * (dx[i]*dx[j] - dy[i]*dy[j]);
   
   for (i=0; i<3; i++)
      for (j=0; j<3; j++)
	 b += (u[i]*v[j] + v[i]*u[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;
}


void print_angles(int i) {
   // prints angles between line in between detector i arms and direction to source

   double line1[3];
   double line2[3];
   double sqrt2 = sqrt(2.0);

   int j;
   // loop over coordinates
   for (j=0; j<3; j++) {
      line1[j]=(detectors[i].lx[j]+detectors[i].ly[j])/sqrt2;
      line2[j]=(detectors[i].lx[j]-detectors[i].ly[j])/sqrt2;
   }
   float cosang1 = 0.0;
   float cosang2 = 0.0;
   for (j=0; j<3; j++) {
      cosang1 += line1[j]*source.vec[j];
      cosang2 += line2[j]*source.vec[j];
   }	 
   
   printf("%s: Cos ang1 = %f, ang1 = %f deg\n", detectors[i].name, cosang1, acos(cosang1)/deg_to_rad);
   printf("%s: Cos ang2 = %f, ang2 = %f deg\n", detectors[i].name, cosang2, acos(cosang2)/deg_to_rad);

   return;
}

// prints antenna pattern functions for detector i
void print_patterns() {
   int i;
   for (i=0; i<3; i++) {
      
      int j;
      
      // vectors defining source frame
      double *u=source.u;
      double *v=source.v;
      double *z=source.vec;
      
      double *x = detectors[i].lx;
      double *y = detectors[i].ly;
      
      printf("Detector %s:\n", detectors[i].name);
      
      // dotted into X and Y detector arms
      double UdotX=0.0;
      double UdotY=0.0;
      
      double VdotX=0.0;
      double VdotY=0.0;
      
      double ZdotX=0.0;
      double ZdotY=0.0;
      
      for (j=0; j<3; j++) {
	 UdotX += u[j]*x[j];
	 UdotY += u[j]*y[j];
	 VdotX += v[j]*x[j];
	 VdotY += v[j]*y[j];
	 ZdotX += z[j]*x[j];
	 ZdotY += z[j]*y[j];
      }
      
      double Fp = UdotX*UdotX+VdotY*VdotY - UdotY*UdotY-VdotX*VdotX;
      double Fc = 2.0*UdotX*VdotX - 2.0*UdotY*VdotY;
      double Fl = 2.0*ZdotX*UdotX - 2.0*ZdotY*UdotY;
      double Fr = 2.0*ZdotX*VdotX - 2.0*ZdotY*VdotY;
      double mag1 = sqrt(Fp*Fp + Fc*Fc);
      double mag2 = sqrt(Fl*Fl + Fr*Fr);
      double ang1 = 180*atan2(Fc, Fp)/M_PI;
      double ang2 = 180*atan2(Fr, Fl)/M_PI;
      
      
      printf("F+ = % 0.4f Fx = % 0.4f  |F| = % 0.4f, phi = % 0.4f degrees\n", Fp, Fc, mag1, ang1);
      printf("Fl = % 0.4f Fr = % 0.4f  |F| = % 0.4f, psi = % 0.4f degrees\n", Fl, Fr, mag2, ang2);
   }
   return;
}
   
   
// 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();

   // 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.

   // I also change the sign of psi to correspond to LAL.  In my
   // notes, increasing psi means rotating the source ellipse
   // CCW as viewed from infinity.  But in LAL,
   // increasing psi means rotating the source ellipse CCW as viewed
   // from earth.

   double iota = 180-inputdata.iota;
   double psi = -inputdata.psi;
   
#ifdef DEBUG
   fprintf(stderr, "get_antenna(): iota = %f degrees\nPsi = %f degrees\n", iota, psi);
#endif
   iota *= deg_to_rad;
   psi *= deg_to_rad;

   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<5; 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
      
      // this is h_cross * (antenna pattern)
      X =         2.0*ci*(alpha*s2p - beta*c2p);
      
      // this is h_plus * (antenna pattern)
      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->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;
}