Improved transit bounding box calculations for RoadRunnerModel to better handle eccentric orbits, updated dependency to meepmeep >= 0.8.0, and bumped version to 2.6.17.

Author: hpparvi

CHANGELOG.md

@@ -1,5 +1,11 @@
 # Changelog
 
+## [2.6.17] - 2025-11-27
+
+### Changed
+
+- Improved RoadRunnerModel transit bounding box calculation for eccentric orbits.
+
 ## [2.6.16] - 2025-10-29
 
 ### Changed

pytransit/__init__.py

@@ -41,7 +41,7 @@
 
 """
 
-__version__ = '2.6.16'
+__version__ = '2.6.17'
 
 # Generic
 # -------

pytransit/models/roadrunner/model_full.py

@@ -1,4 +1,6 @@
 from math import fabs, floor
+
+from meepmeep.tsorbit import bounding_box
 from numba import njit, prange
 from numpy import zeros, dot, ndarray, isnan, nan, ones, full
 
@@ -54,7 +56,7 @@ def rr_full_serial(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: ndarr
     pv_is_good = full(npv, True)
     ldm = zeros((npv, npb, ng))  # Limb darkening means
     xyc = zeros((npv, 2, 5))     # Taylor series coefficients for the (x, y) position
-    hwws = zeros((npv, npb))     # Half-window widths [d]
+    bbs = zeros((npv, nlc, 2))
 
     for ipv in range(npv):
         if isnan(a[ipv]) or (a[ipv] <= 1.0) or (e[ipv] < 0.0) or (isnan(ldp[ipv, 0, 0])):
@@ -79,12 +81,15 @@ def rr_full_serial(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: ndarr
         # ------------------------------------------------------#
         xyc[ipv, :, :] = solve_xy_p5s(0.0, p[ipv], a[ipv], i[ipv], e[ipv], w[ipv])
 
-        # ---------------------------------#
-        # Calculate the half-window widths #
-        # ---------------------------------#
-        hww = 0.5 * d_from_pkaiews(p[ipv], ks[ipv, 0], a[ipv], i[ipv], e[ipv], w[ipv], 1, 14)
+        # -----------------------------#
+        # Calculate the bounding boxes #
+        # -----------------------------#
+        bt1, bt4 = bounding_box(k, xyc)
+        bbs[ipv, :, 0] = bt1
+        bbs[ipv, :, 1] = bt4
         for ilc in range(nlc):
-            hwws[ipv, ilc] = 0.0015 + _exptimes[ilc] + hww
+            bbs[ipv, ilc, 0] -= 0.003 + _exptimes[ilc]
+            bbs[ipv, ilc, 1] += 0.003 + _exptimes[ilc]
 
     # ---------------------------#
     # Calculate the light curves #
@@ -104,8 +109,8 @@ def rr_full_serial(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: ndarr
 
         epoch = floor((times[ipt] - t0[ipv, iep] + 0.5 * p[ipv]) / p[ipv])
         tc = times[ipt] - (t0[ipv, iep] + epoch * p[ipv])
-        if fabs(tc) > hwws[ipv, ilc]:
-            flux[ipv, ipt] = 1.0
+        if not (bbs[ipv, ilc, 0] <= tc <= bbs[ipv, ilc, 1]):
+            flux[ipt] = 1.0
         else:
             for isample in range(1, nsamples[ilc] + 1):
                 time_offset = exptimes[ilc] * ((isample - 0.5) / nsamples[ilc] - 0.5)
@@ -116,6 +121,7 @@ def rr_full_serial(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: ndarr
             flux[ipv, ipt] /= nsamples[ilc]
     return flux
 
+
 @njit(parallel=True, fastmath=False)
 def rr_full_parallel(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: ndarray, i: ndarray, e: ndarray, w: ndarray,
             nlc: int, npb: int, nep: int,
@@ -146,6 +152,7 @@ def rr_full_parallel(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: nda
     ldm = zeros((npv, npb, ng))  # Limb darkening means
     xyc = zeros((npv, 2, 5))  # Taylor series coefficients for the (x, y) position
     hwws = zeros((npv, npb))  # Half-window widths [d]
+    bbs = zeros((npv, nlc, 2))
 
     for ipv in range(npv):
         if isnan(a[ipv]) or (a[ipv] <= 1.0) or (e[ipv] < 0.0) or (isnan(ldp[ipv, 0, 0])):
@@ -171,12 +178,15 @@ def rr_full_parallel(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: nda
         # ------------------------------------------------------#
         xyc[ipv, :, :] = solve_xy_p5s(0.0, p[ipv], a[ipv], i[ipv], e[ipv], w[ipv])
 
-        # ---------------------------------#
-        # Calculate the half-window widths #
-        # ---------------------------------#
-        hww = 0.5 * d_from_pkaiews(p[ipv], ks[ipv, 0], a[ipv], i[ipv], e[ipv], w[ipv], 1, 14)
+        # -----------------------------#
+        # Calculate the bounding boxes #
+        # -----------------------------#
+        bt1, bt4 = bounding_box(k, xyc)
+        bbs[ipv, :, 0] = bt1
+        bbs[ipv, :, 1] = bt4
         for ilc in range(nlc):
-            hwws[ipv, ilc] = 0.0015 + _exptimes[ilc] + hww
+            bbs[ipv, ilc, 0] -= 0.003 + _exptimes[ilc]
+            bbs[ipv, ilc, 1] += 0.003 + _exptimes[ilc]
 
     # ---------------------------#
     # Calculate the light curves #
@@ -196,8 +206,8 @@ def rr_full_parallel(times: ndarray, k: ndarray, t0: ndarray, p: ndarray, a: nda
 
         epoch = floor((times[ipt] - t0[ipv, iep] + 0.5 * p[ipv]) / p[ipv])
         tc = times[ipt] - (t0[ipv, iep] + epoch * p[ipv])
-        if fabs(tc) > hwws[ipv, ilc]:
-            flux[ipv, ipt] = 1.0
+        if not (bbs[ipv, ilc, 0] <= tc <= bbs[ipv, ilc, 1]):
+            flux[ipt] = 1.0
         else:
             for isample in range(1, nsamples[ilc] + 1):
                 time_offset = exptimes[ilc] * ((isample - 0.5) / nsamples[ilc] - 0.5)

pytransit/models/roadrunner/model_simple.py

@@ -2,7 +2,7 @@
 from numba import njit, prange
 from numpy import zeros, dot, ndarray, isnan, full, nan
 
-from meepmeep.xy.position import solve_xy_p5s, pd_t15sc
+from meepmeep.xy.position import solve_xy_p5s, pd_t15sc, bounding_box
 from meepmeep.utils import d_from_pkaiews
 
 from .common import calculate_weights_2d, interpolate_mean_limb_darkening_s
@@ -54,11 +54,12 @@ def rr_simple_serial(times: ndarray, k: float, t0: float, p: float, a: float, i:
     # -----------------------------------------------------#
     xyc[:, :] = solve_xy_p5s(0.0, p, a, i, e, w)
 
-    # --------------------------------#
-    # Calculate the half-window width #
-    # --------------------------------#
-    hww = 0.5 * d_from_pkaiews(p, k, a, i, e, w, 1, 14)
-    hww = 0.0015 + exptimes + hww
+    # ---------------------------#
+    # Calculate the bounding box #
+    # ---------------------------#
+    bt1, bt4 = bounding_box(k, xyc)
+    bt1 -= 0.003 + exptimes
+    bt4 += 0.003 + exptimes
 
     # --------------------------#
     # Calculate the light curve #
@@ -67,7 +68,7 @@ def rr_simple_serial(times: ndarray, k: float, t0: float, p: float, a: float, i:
     for ipt in range(npt):
         epoch = floor((times[ipt] - t0 + 0.5 * p) / p)
         tc = times[ipt] - (t0 + epoch * p)
-        if fabs(tc) > hww:
+        if not (bt1 <= tc <= bt4):
             flux[ipt] = 1.0
         else:
             for isample in range(1, nsamples+ 1):
@@ -111,11 +112,12 @@ def rr_simple_parallel(times: ndarray, k: float, t0: float, p: float, a: float,
     # -----------------------------------------------------#
     xyc[:, :] = solve_xy_p5s(0.0, p, a, i, e, w)
 
-    # --------------------------------#
-    # Calculate the half-window width #
-    # --------------------------------#
-    hww = 0.5 * d_from_pkaiews(p, k, a, i, e, w, 1, 14)
-    hww = 0.0015 + exptimes + hww
+    # ---------------------------#
+    # Calculate the bounding box #
+    # ---------------------------#
+    bt1, bt4 = bounding_box(k, xyc)
+    bt1 -= 0.003 + exptimes
+    bt4 += 0.003 + exptimes
 
     # --------------------------#
     # Calculate the light curve #
@@ -124,7 +126,7 @@ def rr_simple_parallel(times: ndarray, k: float, t0: float, p: float, a: float,
     for ipt in prange(npt):
         epoch = floor((times[ipt] - t0 + 0.5 * p) / p)
         tc = times[ipt] - (t0 + epoch * p)
-        if fabs(tc) > hww:
+        if not (bt1 <= tc <= bt4):
             flux[ipt] = 1.0
         else:
             for isample in range(1, nsamples+ 1):

requirements.txt

@@ -17,4 +17,4 @@ emcee
 pyopencl
 corner
 celerite
-meepmeep
+meepmeep >= 0.8.0