Merge branch 'update-constants-and-gmst-conversion' into 'master'

Resolve issues identified in !408

See merge request !597

bilby/core/utils.py

@@ -10,6 +10,7 @@ import subprocess
 import multiprocessing
 from importlib import import_module
 import json
+import warnings
 
 import numpy as np
 from scipy.interpolate import interp2d
@@ -18,12 +19,11 @@ import pandas as pd
 
 logger = logging.getLogger('bilby')
 
-
-# Constants: values taken from astropy v3.0.4
+# Constants: values taken from LAL cd65f38ce43cef3a1dec217c060de25caf99bf14
 speed_of_light = 299792458.0  # m/s
-parsec = 3.0856775814671916e+16  # m
-solar_mass = 1.9884754153381438e+30  # Kg
-radius_of_earth = 6378100.0  # m
+parsec = 3.085677581491367e+16  # m
+solar_mass = 1.9885469549614615e+30  # Kg
+radius_of_earth = 6378136.6  # m
 
 _TOL = 14
 
@@ -254,6 +254,10 @@ def gps_time_to_gmst(gps_time):
     float: Greenwich mean sidereal time in radians
 
     """
+    warnings.warn(
+        "Function gps_time_to_gmst deprecated, use "
+        "lal.GreenwichMeanSiderealTime(time) instead",
+        DeprecationWarning)
     omega_earth = 2 * np.pi * (1 / 365.2425 + 1) / 86400.
     gps_2000 = 630720013.
     gmst_2000 = (6 + 39. / 60 + 51.251406103947375 / 3600) * np.pi / 12

bilby/gw/utils.py

 from __future__ import division
 import os
 import json
+from math import fmod
 
 import numpy as np
 from scipy.interpolate import interp1d
 import matplotlib.pyplot as plt
 
-from ..core.utils import (gps_time_to_gmst, ra_dec_to_theta_phi,
+from ..core.utils import (ra_dec_to_theta_phi,
                           speed_of_light, logger, run_commandline,
                           check_directory_exists_and_if_not_mkdir,
                           SamplesSummary)
@@ -87,7 +88,7 @@ def time_delay_geocentric(detector1, detector2, ra, dec, time):
     float: Time delay between the two detectors in the geocentric frame
 
     """
-    gmst = gps_time_to_gmst(time)
+    gmst = fmod(lal.GreenwichMeanSiderealTime(time), 2 * np.pi)
     theta, phi = ra_dec_to_theta_phi(ra, dec, gmst)
     omega = np.array([np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)])
     delta_d = detector2 - detector1
@@ -120,8 +121,8 @@ def get_polarization_tensor(ra, dec, time, psi, mode):
     array_like: A 3x3 representation of the polarization_tensor for the specified mode.
 
     """
-    greenwich_mean_sidereal_time = gps_time_to_gmst(time)
-    theta, phi = ra_dec_to_theta_phi(ra, dec, greenwich_mean_sidereal_time)
+    gmst = fmod(lal.GreenwichMeanSiderealTime(time), 2 * np.pi)
+    theta, phi = ra_dec_to_theta_phi(ra, dec, gmst)
     u = np.array([np.cos(phi) * np.cos(theta), np.cos(theta) * np.sin(phi), -np.sin(theta)])
     v = np.array([-np.sin(phi), np.cos(phi), 0])
     m = -u * np.sin(psi) - v * np.cos(psi)

test/gw_likelihood_test.py

@@ -42,7 +42,7 @@ class TestBasicGWTransient(unittest.TestCase):
         """Test log likelihood matches precomputed value"""
         self.likelihood.log_likelihood()
         self.assertAlmostEqual(self.likelihood.log_likelihood(),
-                               -4055.236283345252, 3)
+                               -4055.25243177871, 3)
 
     def test_log_likelihood_ratio(self):
         """Test log likelihood ratio returns the correct value"""
@@ -111,7 +111,7 @@ class TestGWTransient(unittest.TestCase):
         """Test log likelihood matches precomputed value"""
         self.likelihood.log_likelihood()
         self.assertAlmostEqual(self.likelihood.log_likelihood(),
-                               -4055.236283345252, 3)
+                               -4055.25243177871, 3)
 
     def test_log_likelihood_ratio(self):
         """Test log likelihood ratio returns the correct value"""