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