Correct term in peters_g(n, e)

.github/workflows/python-package.yml

@@ -33,7 +33,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.11]
+        python-version: [3.11, 3.12, 3.13]
 
     steps:
     - uses: actions/checkout@v4

changelog.md

@@ -186,4 +186,11 @@ having two when we could just set `f_dom=2 f_orb`
 
 ## 0.5.2
 *TW 02/06/25*
-- Bug fix: Ensure that `plot_sources_on_sc` employs the sensitivity curve settings that are saved in a `Source` class when plotting y-values of sources (thanks to Kierstin Sorensen for pointing out this issue)
+- Bug fix: Ensure that `plot_sources_on_sc` employs the sensitivity curve settings that are saved in a `Source` class when plotting y-values of sources (thanks to Kierstin Sorensen for pointing out this issue)
+
+## 0.5.3
+*TW 09/07/26*
+
+- Bug fix for [#127](https://github.com/TeamLEGWORK/LEGWORK/issues/127) found by @willcerny. This corrects $g(n, e)$ where one term was cubed instead of squared. This causes slight issues for low harmonic modes of very eccentric sources. See the issue for a plot demonstrating where this is focused.
+- Use decorators for visualisation to avoid changing global matplotlib settings for plot
+- Change to use `np.trapezoid` over `np.trapz` and fix numpy to >= 2.0

helpers/generate_g_data.py

@@ -0,0 +1,11 @@
+import numpy as np
+import legwork as lw
+
+n_range = np.arange(1, 10000 + 1).astype(int)
+e_range = np.linspace(0, 1, 1000)
+
+N, E = np.meshgrid(n_range, e_range)
+
+g_vals = lw.utils.peters_g(N, E)
+
+np.save("../legwork/peters_g.npy", g_vals)

legwork/_version.py

@@ -1 +1 @@
-__version__ = "0.5.2"
+__version__ = "0.5.3"

legwork/snr.py

@@ -228,7 +228,7 @@ def snr_circ_evolving(m_1, m_2, f_orb_i, dist, t_obs, n_step, t_merge=None,
         h_f_lisa_2 = psd.power_spectral_density(f=2 * f_orb_evol, t_obs=t_obs, **kwargs)
     h_c_lisa_2 = (2 * f_orb_evol)**2 * h_f_lisa_2
 
-    snr = np.trapz(y=h_c_n_2 / h_c_lisa_2, x=2 * f_orb_evol, axis=1)**0.5
+    snr = np.trapezoid(y=h_c_n_2 / h_c_lisa_2, x=2 * f_orb_evol, axis=1)**0.5
 
     return snr.decompose().value
 
@@ -340,7 +340,7 @@ def snr_ecc_evolving(m_1, m_2, f_orb_i, dist, ecc, harmonics_required, t_obs, n_
     snr_evol = h_c_n_2 / h_c_lisa_2
 
     # integrate, sum and square root to get SNR
-    snr_n_2 = np.trapz(y=snr_evol, x=f_n_evol, axis=1)
+    snr_n_2 = np.trapezoid(y=snr_evol, x=f_n_evol, axis=1)
 
     if ret_snr2_by_harmonic:
         return snr_n_2.decompose().value

legwork/utils.py

@@ -58,7 +58,7 @@ def peters_g(n, e):
     bracket_2 = jv(n-2, n*e) - 2*jv(n, n*e) + jv(n+2, n*e)
     bracket_3 = jv(n, n*e)
 
-    g = n**4/32 * (bracket_1**2 + (1 - e**2) * bracket_2**2 + 4 / (3 * n**3) * bracket_3**2)
+    g = n**4/32 * (bracket_1**2 + (1 - e**2) * bracket_2**2 + 4 / (3 * n**2) * bracket_3**2)
 
     return g
 

legwork/visualisation.py

@@ -5,23 +5,40 @@
 import legwork.psd as psd
 from astropy.visualization import quantity_support
 
-# set the default font and fontsize
-plt.rc('font', family='serif')
-plt.rcParams['text.usetex'] = False
 fs = 24
-
-# update various fontsizes to match
-params = {'figure.figsize': (12, 8),
-          'legend.fontsize': fs,
-          'axes.labelsize': fs,
-          'xtick.labelsize': 0.7 * fs,
-          'ytick.labelsize': 0.7 * fs}
-plt.rcParams.update(params)
+style_params = {
+    'font.family': 'serif',
+    'text.usetex': False,
+    'figure.figsize': (12, 8),
+    'axes.titlesize': 24,
+    'legend.title_fontsize': 18,
+    'legend.fontsize': 16,
+    'axes.labelsize': 24,
+    'xtick.labelsize': 21,
+    'ytick.labelsize': 21,
+    'axes.linewidth': 1.1,
+    'xtick.major.size': 7,
+    'xtick.minor.size': 4,
+    'ytick.major.size': 7,
+    'ytick.minor.size': 4,
+    'xtick.minor.visible': True,
+    'ytick.minor.visible': True,
+}
 
 __all__ = ['plot_1D_dist', 'plot_2D_dist', 'plot_sensitivity_curve',
            'plot_sources_on_sc']
 
 
+# define a decorator for plotting functions which applies the style params in a context environment so that
+# the style is only applied to plots made by legwork and not to any other plots the user might make
+def legwork_plot_style(func):
+    def wrapper(*args, **kwargs):
+        with plt.rc_context(style_params):
+            return func(*args, **kwargs)
+    return wrapper
+
+
+@legwork_plot_style
 def plot_1D_dist(x, weights=None, disttype="hist", log_scale=False, fig=None, ax=None, show=True,
                  figsize=(10, 7), xlabel=None, ylabel=None, xlim=None, ylim=None, color=None, **kwargs):
     """Plot a 1D distribution of ``x``.
@@ -150,6 +167,7 @@ def plot_1D_dist(x, weights=None, disttype="hist", log_scale=False, fig=None, ax
     return fig, ax
 
 
+@legwork_plot_style
 def plot_2D_dist(x, y, weights=None, disttype="scatter", fig=None, ax=None, show=True, figsize=(12, 7),
                  xlabel=None, ylabel=None, xlim=None, ylim=None, log_scale=False,
                  color=None, scatter_s=20, **kwargs):
@@ -281,6 +299,7 @@ def plot_2D_dist(x, y, weights=None, disttype="scatter", fig=None, ax=None, show
     return fig, ax
 
 
+@legwork_plot_style
 def plot_sensitivity_curve(frequency_range=None, y_quantity="ASD", fig=None, ax=None, show=True,
                            figsize=(10, 7), color="#18068b", fill=True, alpha=0.2, linewidth=1, label=None,
                            **kwargs):
@@ -372,6 +391,7 @@ def plot_sensitivity_curve(frequency_range=None, y_quantity="ASD", fig=None, ax=
     return fig, ax
 
 
+@legwork_plot_style
 def plot_sources_on_sc(f_dom, snr, weights=None, snr_cutoff=0, t_obs="auto",
                        instrument="LISA", custom_psd=None, L="auto", approximate_R=False,
                        confusion_noise="auto", fig=None, ax=None, show=True, sc_vis_settings={}, **kwargs):

setup.cfg

@@ -19,7 +19,7 @@ python_requires = >=3.11
 packages = find: 
 install_requires = 
     numba >= 0.58
-    numpy >= 1.26
+    numpy >= 2.0
     astropy >= 6.1
     scipy >= 1.13
     matplotlib >= 3.8