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

src/antenna.c

+// Copyright Bruce Allen 2017
+// Compile with gcc -o antenna antenna.c -lm
+
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -6,7 +9,6 @@
 // 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
@@ -18,15 +20,16 @@ const double deg_to_rad = M_PI/180.0;
 // 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
+// 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
@@ -54,13 +57,25 @@ void print_galaxy_coordinates() {
 
 //  Note these are in the order lattitude,longitude.  See routine
 // "print_galaxy_coordinates" to convert location on sky to Earth
-// location at event time.
+// location at event time.  These are set/used in the function
+// populate_source() below.
 
-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 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];
@@ -74,6 +89,9 @@ struct Detector {
 
    // 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];
 
@@ -82,33 +100,9 @@ struct Detector {
    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) {
@@ -117,13 +111,6 @@ void make_unit_vectors(double *out, double *in) {
    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) {
@@ -132,26 +119,17 @@ void make_u_v_vectors(struct Source *src) {
    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];
-   }
+   // 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) {
@@ -175,25 +153,14 @@ void make_north_east_vectors(struct Detector *det) {
    }
 }
 
-// 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];
-}
-
+// 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;
-
-   // Initially set to zero, fix later
-
-   source.orientation = 0.0;
+   
+   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
@@ -239,12 +206,16 @@ void print_detector(int det) {
 	  "lattitude (north): %f\n"
 	  "longitude (east): %f\n"
 	  "orientation: (CCW from North): %f\n"
-	  "earth center to detector %f %f %f\n\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].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]
 	  );
 }
 
@@ -252,57 +223,85 @@ void print_source() {
    printf("source at:\n"
 	  "lattitude: %f\n"
 	  "longitude: %f\n"
-	  "earth center to source %f %f %f\n\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.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_UUUVVV(int det, double *uu, double *uv, double *vv) {
+void get_UV_combinations(int det, double *alpha, double *beta) {
+  
+   double a=0,b=0;
+   int i, j;
    
-      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;
+   double *su=source.u;
+   double *sv=source.v;
+   double *dx=detectors[det].lx;
+   double *dy=detectors[det].ly;
+   
+   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]);
+   
+   *alpha = a;
+   *beta = b;
 }
-      
 
-   
+// 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;
 
-   double iota=atof(argv[1]);
-   printf("Iota = %f degrees\n", iota);
-   iota *= deg_to_rad;
-   double ci=cos(iota);
-   double si=sin(iota);
-
+   // 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();   
-   populate_detectors();
+   // print_galaxy_coordinates();
    populate_source();
+   // print_source();
    
+   populate_detectors();
    // 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;
+   // 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);
+   iota *= deg_to_rad;
+   psi *= deg_to_rad;
+
+   // get angles needed to calculate antenna response
+   double ci = cos(iota);
+   double si = sin(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 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);
+      double alpha, beta, X, Y;
+
+      // get antenna pattern overlap with source plane
+      get_UV_combinations(i, &alpha, &beta);
+
+      // 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);
+      
+      printf("For detector %s the waveform is w^2 [ %f sin(2wt) + %f cos(2wt) ] \n", detectors[i].name, X, Y);
    }
    return 0;
 }