diff --git a/bilby/core/utils.py b/bilby/core/utils.py index ea66c5614d4eb3aaad0950dd6052096209470a4c..d14c642e928044524ad7105a2774dc092939c0e4 100644 --- a/bilby/core/utils.py +++ b/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 diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index eb2fce9230fd0b1c2e03adc8f6e430f1306673ec..e22641bb0a934e35bbd863ceb25b593f0c17f0d7 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -1,12 +1,13 @@ 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) diff --git a/test/gw_likelihood_test.py b/test/gw_likelihood_test.py index 658044a1f6d6e6e8e068e3274b44674f1df14cbe..2317baaee6baf5919d5d157fb8b9c00583653c33 100644 --- a/test/gw_likelihood_test.py +++ b/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"""