feat: TOI-5398 against PyORBIT and ARoME with matched jitter and offset priors, instantaneous RVs

# PyORBIT + ARoME with the dataset priors of PyAETNA's runs (review M4.4, 2026-10-02): jitter log-uniform in

# [0.01, 100] m/s (PyAETNA: log_jitter uniform in [-4.6, 4.6]) and offset uniform in [-9600, -9200] m/s (PyAETNA:

# normalization uniform in [-9.6, -9.2] km/s). PyORBIT's defaults, used by toi5398_pyorbit_arome.yaml, are a

# uniform jitter in [0.045, 697] m/s and a uniform offset in [-19470, 643] m/s. Everything else as in

# toi5398_pyorbit_arome.yaml; compared with PyAETNA's runs/ccf_instant (EXPTIME 0: ARoME does not integrate over

# the exposures).

# TOI-5398 b with PyORBIT's closed-form RM model (rossitermclaughlin_pyarome: ARoME, Boue et al. 2013, the

# centre of a Gaussian fitted to the CCF), for a like-for-like comparison with PyAETNA's rm_observable: ccf

# (runs/ccf). Same data and priors as configs/toi5398_pyorbit.yaml, with the model of PyAETNA's run:

# - circular orbit and no convective blueshift (12 free parameters, as PyAETNA);

# - line model: non-rotating line of standard deviation 1.28 km/s (PyAETNA's rm_line_beta = 1.81 km/s in

#   Hirano's convention, sqrt(2) x the standard deviation), no macroturbulence;

# - ARoME takes the width of the Gaussian fitted to the out-of-transit CCF as a fixed input (PyAETNA refits it

#   at every v sin i): FWHM 9.235 km/s, the width PyAETNA's line model gives at v sin i = 6.92 km/s (its

#   posterior median) with the prior-mean limb darkening.

# Needs pyarome (github.com/andres-jordan/PyARoME @9cf60eb) and ARoME 1.0.0 (www.astro.up.pt/resources/arome),

# built by ../environment/setup_envs.sh (component arome).

inputs:

  RV_HARPSN:

    file: toi5398_pyorbit.dat

    kind: RV

    ancillary: toi5398_pyorbit_exptime.dat

    boundaries:

      jitter: [0.01, 100.0]

      offset: [-9600.0, -9200.0]

    spaces:

      jitter: Logarithmic   # PyORBIT's log2 space: a flat prior in it is log-uniform, as PyAETNA's

    models:

      - radial_velocities

      - rm_model

common:

  planets:

    b:

      orbit: circular

      parametrization: Standard

      use_time_inferior_conjunction: True

      boundaries:

        P: [10.58, 10.60]

        Tc: [2459616.48, 2459616.50]

        K: [0.0, 50.0]

        R_Rs: [0.05, 0.15]

        b: [0.0, 1.0]

        lambda: [-180.0, 180.0]

      spaces:

        P: Linear

        K: Linear

      priors:

        P: ['Gaussian', 10.590547, 0.000012]

        Tc: ['Gaussian', 2459616.49232, 0.00022]

        K: ['Gaussian', 15.7, 1.5]

        R_Rs: ['Gaussian', 0.0899, 0.0007]

        b: ['Gaussian', 0.27, 0.07]

  star:

    star_parameters:

      convective_order: 0

      boundaries:

        v_sini: [0.0, 20.0]

      priors:

        density: ['Gaussian', 0.971, 0.05]    # solar units

        v_sini: ['Gaussian', 7.5, 0.6]        # km/s

    limb_darkening:

      model: ld_quadratic

      priors:

        ld_c1: ['Gaussian', 0.5760, 0.0149]

        ld_c2: ['Gaussian', 0.1377, 0.0318]

models:

  radial_velocities:

    planets:

      - b

  rm_model:

    kind: rossitermclaughlin_pyarome

    limb_darkening: limb_darkening

    measurement_technique: ccf

    instrumental_broadening: 1.2799     # km/s, standard deviation of the non-rotating line

    measured_ccf_width: 9.2346          # km/s, FWHM of the Gaussian fitted to the CCF

    macroturbulence: 0.0

    planets:

      - b

parameters:

  Tref: 2460029.5

  cpu_threads: 16          # as for juliet in the K2-140 benchmark

solver:

  # nested sampling (pyorbit_run dynesty_static), as the other benchmarks: valid here

  # because v sin i is sampled directly (the paper used emcee for its prior on a

  # derived v sin i)

  nested_sampling:

    nlive: 500

    dlogz: 0.1          # bound and sample: PyORBIT defaults (multi, auto)

  pyde:

    ngen: 20000

    npop_mult: 4

  emcee:

    npop_mult: 4

    nsteps: 50000

    nburn: 20000

    nsave: 10000

    thin: 100

  recenter_bounds: True

# TOI-5398: HARPS-N RVs of the 2023-03-25 transit (Mantovan et al. 2024, CDS J/A+A/684/L17); RV in m/s; EXPTIME 0 (instantaneous model)

BJD Y SIGMA EXPTIME VISIT INSTRUMENT

2460029.352911 -9421.096 6.572 0.0 1 1

2460029.360121 -9412.300 6.644 0.0 1 1

2460029.367413 -9412.513 5.391 0.0 1 1

2460029.374681 -9417.624 5.286 0.0 1 1

2460029.381683 -9417.289 5.452 0.0 1 1

2460029.389124 -9413.045 6.132 0.0 1 1

2460029.396369 -9421.042 5.853 0.0 1 1

2460029.403522 -9414.836 5.589 0.0 1 1

2460029.410373 -9423.115 5.836 0.0 1 1

2460029.418058 -9422.045 6.337 0.0 1 1

2460029.425130 -9427.087 6.326 0.0 1 1

2460029.432641 -9409.101 6.465 0.0 1 1

2460029.440071 -9403.621 5.728 0.0 1 1

2460029.446934 -9394.796 5.356 0.0 1 1

2460029.454225 -9373.578 5.009 0.0 1 1

2460029.461586 -9367.374 5.082 0.0 1 1

2460029.468657 -9366.677 5.113 0.0 1 1

2460029.476041 -9357.724 4.787 0.0 1 1

2460029.483171 -9357.447 4.525 0.0 1 1

2460029.490369 -9370.988 4.752 0.0 1 1

2460029.497718 -9378.465 5.103 0.0 1 1

2460029.505125 -9392.763 4.749 0.0 1 1

2460029.511965 -9401.925 4.990 0.0 1 1

2460029.519477 -9415.941 5.200 0.0 1 1

2460029.526768 -9417.712 4.964 0.0 1 1

2460029.533932 -9432.976 4.928 0.0 1 1

2460029.541142 -9450.814 4.877 0.0 1 1

2460029.548341 -9454.765 4.589 0.0 1 1

2460029.555655 -9456.060 4.779 0.0 1 1

2460029.562993 -9465.437 5.111 0.0 1 1

2460029.570330 -9470.478 4.759 0.0 1 1

2460029.577309 -9470.340 4.986 0.0 1 1

2460029.584716 -9470.104 5.230 0.0 1 1

2460029.591845 -9462.444 5.226 0.0 1 1

2460029.599264 -9440.003 5.216 0.0 1 1

2460029.606498 -9432.490 5.451 0.0 1 1

2460029.613627 -9421.743 5.234 0.0 1 1

2460029.621034 -9408.090 5.504 0.0 1 1

2460029.628267 -9419.133 5.371 0.0 1 1

2460029.635466 -9424.711 5.329 0.0 1 1

2460029.642618 -9427.364 5.648 0.0 1 1

2460029.649875 -9421.984 5.580 0.0 1 1

2460029.657143 -9414.513 6.001 0.0 1 1

2460029.664365 -9434.382 6.127 0.0 1 1

2460029.671587 -9411.301 6.967 0.0 1 1

2460029.678971 -9420.501 6.169 0.0 1 1

2460029.686135 -9415.935 6.746 0.0 1 1

Reads runs/pyorbit_ns/toi5398_pyorbit/dynesty_static/results.p (dynesty's

results, pickled by PyORBIT) and the parameter ids printed in runs/pyorbit_ns/run.out

("<name> id: <i> s:<space> ..."). Every parameter of this configuration is in

linear space, so the samples are the parameter values themselves. Writes

linear space, so the samples are the parameter values themselves; parameters in a logarithmic space

(the matched configuration's jitter) are transformed back. Writes

reference/pyorbit_samples.npz (equal-weight resampling, fixed seed) and

reference/pyorbit_evidence.json.

HERE = Path(__file__).parent

# PyORBIT's sampling spaces and their inverse transformations (pyorbit.subroutines.common)

SPACES = {

    "Linear": lambda x: x,

    "Logarithmic": lambda x: 2.0**x,

    "Log_Base2": lambda x: 2.0**x,

    "Log_Natural": np.exp,

    "Log_Base10": lambda x: 10.0**x,

}

def main(run: Path, config: str, prefix: str) -> None:

    ids = {}

    ids: dict[int, str] = {}

    spaces: dict[int, str] = {}

    for line in (run / "run.out").read_text().splitlines():

        m = re.match(r"^(\S+)\s+id:\s+(\d+)\s+s:(\S+)", line)

        if m:

            if m.group(3) != "Linear":

                raise SystemExit(f"{m.group(1)} is not in linear space")

            if m.group(3) not in SPACES:

                raise SystemExit(f"{m.group(1)}: unknown PyORBIT space {m.group(3)}")

            ids[int(m.group(2))] = m.group(1)

            spaces[int(m.group(2))] = m.group(3)

    with open(run / config / "dynesty_static" / "results.p", "rb") as f:

        res = pickle.load(f)

    assert sorted(ids) == list(range(res.samples.shape[1]))

    rng = np.random.default_rng(1)

    n = int(np.round(1 / np.sum((weights / weights.sum()) ** 2)))  # Kish ESS

    idx = rng.choice(len(weights), size=n, p=weights / weights.sum())

    samples = {ids[i]: res.samples[idx, i] for i in ids}

    samples = {ids[i]: SPACES[spaces[i]](res.samples[idx, i]) for i in ids}

    (HERE / "reference").mkdir(exist_ok=True)

    np.savez_compressed(HERE / "reference" / f"{prefix}_samples.npz", **samples)

    (HERE / "reference" / f"{prefix}_evidence.json").write_text(

- data/toi5398.dat, PyAETNA's format: BJD Y SIGMA EXPTIME VISIT INSTRUMENT,

  with Y and SIGMA in m/s (the systemic velocity is then fitted in km/s, as

  PyAETNA expects) and EXPTIME = 600 s;

- data/toi5398_instant.dat, the same with EXPTIME = 0 (for the ARoME comparison: ARoME does not integrate);

- data/toi5398_pyorbit.dat, PyORBIT's RV format: BJD, RV [m/s], error [m/s],

  jitter flag 0, offset flag 0;

- data/toi5398_pyorbit_exptime.dat, PyORBIT's ancillary file with the exposure

        f.write("BJD Y SIGMA EXPTIME VISIT INSTRUMENT\n")

        for t, rv, err in rows:

            f.write(f"{t:.6f} {rv} {err} {EXPTIME:.1f} 1 1\n")

    # the same RVs as instantaneous measurements (EXPTIME = 0), for the comparison with PyORBIT's ARoME model,

    # which does not integrate over the exposures (PyAETNA averages the RM anomaly when EXPTIME > 0)

    with open(HERE / "toi5398_instant.dat", "w") as f:

        f.write(

            "# TOI-5398: HARPS-N RVs of the 2023-03-25 transit (Mantovan et al. 2024, CDS "

            "J/A+A/684/L17); RV in m/s; EXPTIME 0 (instantaneous model)\n"

        )

        f.write("BJD Y SIGMA EXPTIME VISIT INSTRUMENT\n")

        for t, rv, err in rows:

            f.write(f"{t:.6f} {rv} {err} 0.0 1 1\n")

    with open(HERE / "toi5398_pyorbit.dat", "w") as f:

        f.write("# epoch value error jitter_flag offset_flag\n")

        for t, rv, err in rows:

{

 "log_evidence": -160.43848983532055,

 "log_evidence_err": 0.26148149531965226,

 "niter": 7205,

 "ncall": 190645,

 "n_posterior_samples": 2550

}

 No newline at end of file