Merge pull request #108 from TomWagg/angles

Correct position coordinates

changelog.md

@@ -158,4 +158,9 @@ having two when we could just set `f_dom=2 f_orb`
 
 ## 0.4.4
 *TW 14/09/22*
-- [Issues [#106](https://github.com/TeamLEGWORK/LEGWORK/issues/106)] Fix sampling of inclination and polarisation when they are not supplied to `Source`
+- [Issues [#106](https://github.com/TeamLEGWORK/LEGWORK/issues/106)] Fix sampling of inclination and polarisation when they are not supplied to `Source`
+
+## 0.4.5
+*TW 27/09/22*
+- Change positions coordinate system to barycentric rather than heliocentric
+- Ensure `theta` is a co-latitude in `amplitude_modulation`

legwork/_version.py

@@ -1 +1 @@
-__version__ = "0.4.4"
+__version__ = "0.4.5"

legwork/source.py

@@ -136,13 +136,15 @@ def __init__(self, m_1, m_2, ecc, dist, n_proc=1, f_orb=None, a=None, position=N
                                  "modulation is only valued for circular sources")
 
             # ensure position is in the correct coordinate frame
-            position = position.transform_to("heliocentrictrueecliptic")
+            position = position.transform_to("barycentrictrueecliptic")
+
+            # barycentrictrueecliptic
 
             # ensure that the position, polarisation, and inclination
             # quantities are at least 1d for masking later on
             lon, lat, polarisation, inclination = np.atleast_1d(position.lon, position.lat,
                                                                 polarisation, inclination)
-            position = SkyCoord(lon=lon, lat=lat, distance=dist, frame='heliocentrictrueecliptic')
+            position = SkyCoord(lon=lon, lat=lat, distance=dist, frame='barycentrictrueecliptic')
 
         # calculate whichever one wasn't supplied
         f_orb = utils.get_f_orb_from_a(a, m_1, m_2) if f_orb is None else f_orb

legwork/strain.py

@@ -1,5 +1,6 @@
 """Computes several types of gravitational wave strains"""
 
+import astropy.units as u
 import astropy.constants as c
 from legwork import utils
 import numpy as np
@@ -36,7 +37,8 @@ def amplitude_modulation(position, polarisation, inclination):
     modulation : `float/array`
         modulation to apply to strain from detector response
     """
-    theta, phi = position.lat, position.lon
+    # the pi/2 subtraction ensures that theta is a co-latitude to match Cornish+03
+    theta, phi = (np.pi/2 * u.rad) - position.lat, position.lon
 
     # a_plus/a_cross as defined in Robson+19 Eq.15 and Babak+21 Eq. 67
     a_plus = (1 + np.cos(inclination)**2) / 2

legwork/tests/test_sources.py

@@ -141,7 +141,7 @@ def test_amplitude_modulation_h_0_n(self):
         phis = np.random.uniform(0, 2 * np.pi, n_values) * u.rad
         psis = np.random.uniform(0, 2 * np.pi, n_values) * u.rad
 
-        positions = SkyCoord(phis, thetas, distance=dist, frame='heliocentrictrueecliptic')
+        positions = SkyCoord(phis, thetas, distance=dist, frame='barycentrictrueecliptic')
 
         sources = source.Source(m_1=m_1, m_2=m_2, f_orb=f_orb,
                                 ecc=ecc, dist=dist,
@@ -169,7 +169,7 @@ def test_amplitude_modulation_h_c_n(self):
         phis = np.random.uniform(0, 2 * np.pi, n_values) * u.rad
         psis = np.random.uniform(0, 2 * np.pi, n_values) * u.rad
 
-        positions = SkyCoord(phis, thetas, distance=dist, frame='heliocentrictrueecliptic')
+        positions = SkyCoord(phis, thetas, distance=dist, frame='barycentrictrueecliptic')
 
         sources = source.Source(m_1=m_1, m_2=m_2, f_orb=f_orb,
                                 ecc=ecc, dist=dist,