[antonsoo]

For multi-planet fits with fit_secondary_mass=True, batching identical parameter vectors changes the predicted astrometry and likelihood:

import os
import numpy as np
import orbitize
from orbitize import read_input, sampler, system

data = read_input.read_file(os.path.join(orbitize.DATADIR, "test_val_multi.csv"))
s = system.System(2, data, 1.47, 24.3, mass_err=.11, plx_err=.7,
                  fit_secondary_mass=True, tau_ref_epoch=50000)
p = np.array([72.77401, .041116819, .56322372, 3.5251172, 4.2904768,
              .094234377, 45.418411, .0014317369, .56322372, 3.1016846,
              4.2904768, .34033456, 24.589758, .007, .008, 1.47])
mc = sampler.MCMC(s, num_walkers=64, num_threads=1)
print(mc._logl(p))
print(mc._logl(np.repeat(p[None, :], 3, axis=0)))

Before: -23.46248 for the scalar and [-32.23987, -32.23987, -32.23987] for the batch. After: all are -23.46248 (rounded).

np.where on the planet-by-sample ordering mask returns the same body index once per matching sample. Select each affected body once, then apply its correction only to samples where that companion is interior. Deduplicating alone still gives incorrect results when orbital ordering varies within the batch.

Regressions cover two/three companions, identical/mixed/equal semimajor axes, and singleton/three-sample batches. Native checks on Python 3.10 and 3.12 compare the bundled astrometry with separately assembled scalar Kepler predictions and Gaussian likelihoods, including HDF5 roundtrips. All working controls and scalar predictions stay exact; a clean native wheel reproduces the results. The change composes with #451 and preserves its checkpoint workloads.

Related to #392; this fixes the multi-planet astrometric batch case, without claiming every configuration in that report is resolved. #426's relative-RV fix still reproduces this separate error.

@sblunt @semaphoreP, could you review this change?