Modified the library by changing sign conventions for psi to correspond to LAL.

Also added some additional debugging info and a detector response tensor that
can be directly compared to the one in LAL (they are equal, to good precision).

lib/Makefile

-.DEFAULT_GOAL := libantenna.a
+all: libantenna.a antenna_lib.o antenna_main
 
 antenna_main: antenna_lib.o antenna_main.cpp antenna_lib.h
 	g++ -o antenna_main antenna_main.cpp antenna_lib.o -lm

lib/antenna_lib.cpp

@@ -84,7 +84,8 @@ struct Source {
    double vec[3];
 
    // Unit vectors defining plane perpendicular to the source
-   // direction.  U points east, V points north
+   // 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];
 };
@@ -106,13 +107,19 @@ struct Detector {
    double vec[3];
 
    // Unit vectors pointing along the north and east directions at the
-   // detector site
-   double north[3];
+   // 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]; 
+   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, in that order
@@ -134,7 +141,7 @@ void make_u_v_vectors(struct Source *src) {
    double lon = src->location[1];
 
    // construct unit vectors perpendicular to line of sight
-   // u points east
+   // u points east (viewed from infinity)
    src->u[0] = -sin(lon);
    src->u[1] = cos(lon);
    src->u[2] = 0.0;
@@ -151,14 +158,20 @@ 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);
+
+   // 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;
@@ -168,19 +181,32 @@ void make_north_east_vectors(struct Detector *det) {
    }
 }
 
+// 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 a detector
+   // 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);
 }
@@ -204,6 +230,10 @@ void populate_detectors(){
    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!
@@ -220,7 +250,11 @@ void populate_detectors(){
    for (i=0;i<3;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) {
@@ -239,6 +273,17 @@ void print_detector(int det) {
 	  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() {
@@ -262,8 +307,8 @@ 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 *u=source.u;
+   double *v=source.v;
    double *n=source.vec;
    
    double *dx=detectors[det].lx;
@@ -272,11 +317,11 @@ void get_UV_combinations(int det, double *alpha, double *beta, double *rdotn) {
    
    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]);
+	 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 += (su[i]*sv[j] + sv[i]*su[j]) * (dx[i]*dx[j] - dy[i]*dy[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
@@ -305,13 +350,6 @@ void get_antenna(struct OutputStruct *out, struct InputStruct *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
@@ -321,8 +359,23 @@ void get_antenna(struct OutputStruct *out, struct InputStruct *in) {
    // 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 ci =  cos(iota);
    double c2p = cos(2*psi);
    double s2p = sin(2*psi);
 
@@ -338,8 +391,13 @@ void get_antenna(struct OutputStruct *out, struct InputStruct *in) {
       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;
       
@@ -358,10 +416,10 @@ void get_antenna(struct OutputStruct *out, struct InputStruct *in) {
       // 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-2*psi*ci/deg_to_rad;
+      out->phase[i] = ang;
       out->dt[i] = dt;
 
-#define DEBUG      
 #ifdef DEBUG
       // degree character in UTF-8 character set (most likely terminal type!)
       int deg1=0xC2, deg2=0xB0;
@@ -376,8 +434,6 @@ void get_antenna(struct OutputStruct *out, struct InputStruct *in) {
    return;
 }
 
-
-
 void print_angles(int i) {
    // prints angles between line in between detector i arms and direction to source
 

lib/antenna_lib.h

 // Double inclusion protection
 #ifndef ANTENNA_LIB_H
 #define ANTENNA_LIB_H
+
+void print_detector(int i);
+void print_source();
+
 // 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
+   // from North when viewed from the EARTH
    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

lib/antenna_main.cpp

@@ -2,9 +2,20 @@
 // 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
+
+// (a) sign conventions for h, which arm is positive CONFIRMED.  VIRGO
+//  AND LIGO DEFINE THE POSITIVE ARM AS "X" where X defined by right
+//  hand rule with two arms and Z away from center of earth.
+
+// (b) direction and origin conventions for the polarization axis The
+// convention for LAL is that psi is the angle that the source plane
+// is rotated CCW as viewed from the Earth.  The convention in my
+// calculation is the opposite.  So in my code I have changed the sign
+// inside get_antenna() so that the arguments now correspond to LAL
+// ones.  Thus in this code, the conventions correspond to LAL.
+
+// IOTA is zero when orbital angular momentum points to earth
+// PSI is orbital plane rotation CCW as viewed from Earth
 
 #include <math.h>
 #include <stdio.h>
@@ -21,7 +32,7 @@ int main(int argc, char *argv[]) {
    // to pass data in and out
    struct InputStruct myinput;
    struct OutputStruct myoutput;
-   
+
    // check syntax crudely, issue usage message
    if (argc != 3) {
       fprintf(stderr,
@@ -30,10 +41,11 @@ int main(int argc, char *argv[]) {
 	      argv[0]);
       exit(1);
    }
-
-   // pass inclination angle, polarization axis, orientation offsets
    myinput.iota = atof(argv[1]);
    myinput.psi = atof(argv[2]);
+
+
+   // pass inclination angle, polarization axis, orientation offsets
    for (i=0;i<3;i++) myinput.orientation[i]=0.0;
    
    printf("Iota = %f degrees\nPsi = %f degrees\n", myinput.iota, myinput.psi);
@@ -41,14 +53,19 @@ int main(int argc, char *argv[]) {
    // now compute responses
    get_antenna(&myoutput, &myinput);
 
+   for (i=0; i<3; i++) print_detector(i);
+   print_source();
+
    // 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", myoutput.name[i], myoutput.amp[i], myoutput.dt[i], myoutput.phase[i], deg1, deg2);
+      printf("For detector %s the waveform is %.3f w^2 sin(2w[t%+.1f ms]%+.1f%c%c)\n",
+	     myoutput.name[i], myoutput.amp[i], myoutput.dt[i], myoutput.phase[i], deg1, deg2);
 
-   for (i=0; i<3; i++) print_angles(i);
+   // prints angles between line in between detector i arms and direction to source
+   // for (i=0; i<3; i++) print_angles(i);
 
    return 0;
 }