API reference

def cli():
    """
    Read the settings-file path from the command line, then run

    Kept separate from main() so that importing Morrigan and calling
    main() from another program never inspects that program's own
    command line.
    """
    parser = argparse.ArgumentParser(description='Run the Morrigan giant-impact model')
    parser.add_argument(
        '-c', '--config', default=DEFAULT_CONFIG, help='path to the .toml settings file'
    )
    main(parser.parse_args().config)

def main(config_path):
    """
    Run every system described by a settings file

    Parameters
    ----------
    config_path : str
        Path to the .toml settings file

    Notes
    -----
    A batch of more than one system runs on a process pool. Python starts
    those workers by re-importing the module that launched them, so a
    caller that reaches this function from an unguarded top level would
    have each worker launch the batch again. Put the call behind
    ``if __name__ == '__main__':``; a batch that cannot start a pool runs
    its systems one after another instead, which gives the same numbers
    because each system seeds itself from its own index.
    """
    # a worker re-importing an unguarded caller arrives back here, where it
    # would launch the batch a second time; it has its own task to get on with
    if current_process().name != 'MainProcess':
        return

    # each disk is initialised with the same conditions
    config = read_config(config_path)

    # number of systems to run (defaults to a single run, unless specified)
    ndisk = config.get('batch', {}).get('ndisk', 1)

    # a relative save_directory is taken from the working directory, which is
    # not necessarily where the settings file lives, so say where results land
    save_directory = config['run_simulation']['save_directory']
    os.makedirs(save_directory, exist_ok=True)
    print(f'Writing results to {os.path.abspath(save_directory)}')

    if ndisk <= 1:
        # single run so skip multiprocessing entirely, allows for debugging
        start = time.time()
        result = run_once(0, config)
        end = time.time()

        summary = Table(rows=[result], names=['run_idx', 'runtime_s', 'n_survivors'])
        ascii.write(
            summary,
            os.path.join(config['run_simulation']['save_directory'], 'batch_summary.csv'),
            format='fixed_width',
            overwrite=True,
        )
        print(f'Ran 1 system in {round(end - start, 3)}s')

    else:
        # number of cores to use, leaves one free by default
        nproc = config.get('batch', {}).get('nproc', max(1, cpu_count() - 1))

        worker = partial(run_once, config=config)

        start = time.time()
        try:
            with Pool(processes=nproc) as pool:
                results = pool.map(worker, range(ndisk))
        except RuntimeError:
            # the caller's top level is unguarded, so a worker cannot re-import
            # it safely; run the systems in turn instead of failing the batch
            print(
                'Could not start a process pool, running the systems one at a time. '
                "Put the call behind if __name__ == '__main__': to run them in parallel."
            )
            results = [run_once(i, config) for i in range(ndisk)]
        end = time.time()

        # high-level statistics for each system (remaining planets, ids, etc)
        summary = Table(rows=results, names=['run_idx', 'runtime_s', 'n_survivors'])
        ascii.write(
            summary,
            os.path.join(config['run_simulation']['save_directory'], 'batch_summary.csv'),
            format='fixed_width',
            overwrite=True,
        )
        print(f'Ran {ndisk} systems in {round(end - start, 3)}s')

def read_config(config_path=DEFAULT_CONFIG):
    """
    Read a settings file

    Parameters
    ----------
    config_path : str
        Path to the .toml settings file

    Returns
    -------
    config : dict
        Parsed settings

    Raises
    ------
    FileNotFoundError
        If no settings file is at the given path. The default is a bare
        file name, so it only resolves when the command runs from a
        checkout root; the message names the path that was tried.
    """
    if not os.path.isfile(config_path):
        raise FileNotFoundError(
            f'No settings file at {os.path.abspath(config_path)}. '
            'Pass one with -c, for example: morrigan -c initialise.toml'
        )
    with open(config_path, 'r') as f:
        return toml.load(f)

def run_system(
    seed,
    masses,
    eccentricity,
    inner_edge,
    spacing,
    density,
    impact_angle,
    evolution_time,
    inner_cutoff,
    stellar_mass,
):
    """Evolve one system and return its survivors and impact history in memory.

    This is the entry point an in-process caller uses instead of the
    file-writing command line. Inputs are SI apart from the stellar mass,
    in solar masses, and the evolution time, in Gyr; masses are in kg and
    lengths in metres. The returned quantities are SI with times in years.

    Parameters
    ----------
    seed : int
        Seed for this system's random draws.
    masses : sequence of float
        Initial embryo masses [kg].
    eccentricity : float
        Initial eccentricity shared by the embryos.
    inner_edge : float
        Orbit of the innermost embryo [m].
    spacing : float
        Spacing between adjacent embryos [mutual Hill radii].
    density : float
        Bulk density shared by the embryos [kg m-3].
    impact_angle : float
        Impact angle [deg]; its sine is the impact parameter.
    evolution_time : float
        Length of the dynamical evolution [Gyr].
    inner_cutoff : float
        Orbit inside which a body is taken to have fallen into the star [m].
    stellar_mass : float
        Host stellar mass [Msun].

    Returns
    -------
    dict
        ``survivors``: one record per surviving body with ``id``,
        ``mass_initial`` [kg], ``a_initial`` [m], ``mass_final`` [kg] and
        ``a_final`` [m]. ``impacts``: a dict keyed by surviving-body id,
        each value the list of impacts that body experienced in time
        order, every impact a dict in the shared impact schema (SI, time
        in years, the merged mass reported as the plain sum of the two
        bodies).
    """
    n = len(masses)
    config = {
        'run_simulation': {
            't': 0.0,
            't_ref': 0.0,
            't_event': 0.0,
            'flag_event': 1,
            'a_min': inner_cutoff / au2m,  # [AU]
            'max_time': evolution_time,  # [Gyr]
            'random_seed': seed,
            'save_directory': '',  # unused: collect mode writes nothing
        },
        'init_par': {
            'N': n,
            'e': eccentricity,
            'impact_angle': impact_angle,
            'Mp': [m / M_earth for m in masses],  # [M_earth]
            'Ms': stellar_mass,  # [Msun]
            'rho_p': density,
            'inner_edge': inner_edge / au2m,  # [AU]
            'spacing': spacing,
        },
    }

    raw = run_once(0, config, collect=True)

    year_sec = gyr2sec / 1.0e9  # match the model's own 365-day year
    b = np.sin(np.deg2rad(impact_angle))  # constant impact parameter of the run

    impacts = {}
    for m in raw['mergers']:
        M_t, M_i = m['M_target_before'], m['M_impactor_before']
        R_t, R_i = m['R_target_before'], m['R_impactor']
        target = int(m['id_target'])
        record = {
            'time': m['t'] / year_sec,  # [yr]
            'M_target_before': M_t,
            'M_impactor': M_i,
            # perfect-merger mass: the plain sum, leaving atmospheric loss to the caller
            # Carry the model's own merged mass rather than re-adding the two
            # fields above: recomputing the sum here would make any consumer
            # check of mass closure a tautology and hide a mass sink in the
            # dynamics behind a record that always looks exact.
            'M_merged_after': m['M_merged_after'],
            'v_impact': m['v_c'],
            # mutual escape velocity of the pair, the floor the collision speed sits above
            'v_esc': np.sqrt(2.0 * G * (M_t + M_i) / (R_t + R_i)),
            'impact_parameter': b,
            'R_target_before': R_t,
            'R_impactor': R_i,
            # shared bulk density recovered from mass and radius
            'rho_target': M_t / ((4.0 / 3.0) * np.pi * R_t**3),
            'rho_impactor': M_i / ((4.0 / 3.0) * np.pi * R_i**3),
            'a_before': m['a_before'],  # [m]
            'a_after': m['a_final_AU'] * au2m,  # [m]
            # both elements before and after, so a consumer following a planet on a
            # different orbit can apply the change rather than the absolute value
            'e_before': m['e_before'],
            'e_after': m['e_after'],
            'id_target': target,
            'id_impactor': int(m['id_impactor']),
        }
        impacts.setdefault(target, []).append(record)
    for chain in impacts.values():
        chain.sort(key=lambda r: r['time'])

    survivors = []
    for pid, mf, af in zip(raw['planet_id'], raw['masses'], raw['a']):
        pid = int(pid)
        mass_initial, a_initial = raw['initial_state'][pid]
        survivors.append(
            {
                'id': pid,
                'mass_initial': mass_initial,
                'a_initial': a_initial,
                'mass_final': float(mf),
                'a_final': float(af),
            }
        )

    # a body that was a target early and then died carries a partial history that
    # belongs to no survivor, and can even hold an unbound post-merge orbit, so keep
    # only the survivors' chains; each survivor is present, with an empty list if it
    # never merged, so a caller can always look one up
    survivor_ids = {s['id'] for s in survivors}
    impacts = {pid: impacts.get(pid, []) for pid in survivor_ids}

    return {'survivors': survivors, 'impacts': impacts}

def system_seed(base_seed, run_idx):
    """Return the seed for one system of a batch.

    Mixes the settings-file seed with the run index rather than adding
    them. Adding makes neighbouring ensembles overlap: run k of seed s
    draws the same stream as run k-1 of seed s+1, so two ensembles one
    apart share all but one of their systems, and any scatter measured
    across seeds comes out far too small. Mixing decorrelates them while
    keeping each system exactly reproducible from its settings file,
    which is the property that matters for rerunning a result.

    Parameters
    ----------
    base_seed : int
        Seed from the settings file, shared by the whole batch.
    run_idx : int
        Index of this system within the batch.

    Returns
    -------
    int
        Seed for this system, in numpy's valid seed range.
    """
    return int(np.random.SeedSequence([base_seed, run_idx]).generate_state(1)[0])

def crossing_pair(
    ap, Mp, Rp, Ms, ecc, ecc_vec, g, beta, interact, N, t, t_ref
):  # identify crossing pair from triplet, return pair and t_event
    """
    Identify interacting pair of planets and calculate when their orbits will cross

    Parameters
    ----------
    ap : list
        Semi major axes of all planets in system [m]
    Mp : list
        Masses of all planets in system [kg]
    Rp : list
        Radii of all planets in system [m]
    ecc : list
        Eccentriticies of all planets
    ecc_vec : array
        Eccentricity vector components per planet
    g : array
        Secular eigenfrequency
    beta : array
        Phase angle
    interact : boolean list
        Flags for which planets are live and interacting
    N : int
        Number of planets
    t : float
        Time in simulation [s]
    t_ref : float

    Returns
    -------
    icross : int
        Index of the inner planet in the crossing pair
    t_event : float
        Time of next crossing event
    """

    # evaluate tau_cross for all planets
    Tcross = np.full(
        N - 1, 1e20
    )  # initial array to eventually store the predicted crossing times for every pair, initialised at 'never'
    Naffect = np.ones(N, dtype=int)  # the interacting planets, used to compute K factor in eq 5
    interacting_pair = np.arange(N - 1)  # stores which pair of the triplet interact

    # calculate current orbital separation
    bKmin = 2.5  # result from N-body simulation, Kokubo et al 2025
    bKij = np.empty(N - 1)
    for i in range(N - 1):
        a_mean = 0.5 * (ap[i] + ap[i + 1])
        # maximum separation before ejection
        rKij = (
            esc_ecc(Ms, Mp[i], Mp[i + 1], Rp[i], Rp[i + 1], a_mean) * a_mean
        )  # physical separation
        bKij[i] = (ap[i + 1] - ap[i]) / rKij  # current gap relative to physical separation

    # 2-body case
    if N == 2:
        aM = (Mp[0] * ap[0] + Mp[1] * ap[1]) / (Mp[0] + Mp[1])  # mean semi-major axis
        h = hill_sphere(aM, Mp[0] + Mp[1], Ms) / aM
        # checking for stability again here
        EJbef = (
            5.0 / 8.0 * (ecc[0] ** 2 + ecc[1] ** 2) / h**2
            - 3.0 / 8.0 * ((ap[0] - ap[1]) / (h * aM)) ** 2
            + 4.5
        )
        if EJbef > 0.0:  # pair is unstable
            return (
                0,
                1.5 * t,
            )  # schedule a prompt event; a stable pair instead falls through to the 'never' sentinel below

    # calculate orbit crossing time for each pair with the original eccentricities
    for i in range(N - 2):
        if not all(interact[i : i + 3]):
            continue  # skip triplets where a planet is not live or interacting

        # determine 'packed planets'
        group = np.zeros(N, dtype=bool)
        # true if conditions are satisfied, otherwise false
        group[: i + 1] = (bKij[: i + 1] < bKmin) & interact[: i + 1]
        group[i + 1 :] = (bKij[i : N - 1] < bKmin) & interact[i + 1 :]

        #'let up to 2 neighboring planets inside and outside the triplet affect the number of resonances'
        # find indices of planets to the left (i_in) and right (i_out) of the triplet
        non_group_in = np.where(~group[: i + 1])[
            0
        ]  # indices if planets not in the triplet, flips boolean
        if len(non_group_in) == 0 or all(
            group[: i + 1]
        ):  # if there are no planets, or all of them are packed
            i_in = 0
        else:  # otherwise, i_in is the index of the planet closest to inner edge of triplet
            i_in = min(non_group_in[-1] + 1, i)

        # similarly for outer planet
        non_group_out = np.where(~group[i + 1 :])[0]
        if len(non_group_out) == 0 or all(group[i + 1 :]):
            i_out = N - 1  # if everything is packed, the group extends to the end of the system
        else:  # index of planet closest to outer edge of the triplet
            i_out = max(non_group_out[0] + i, i + 1)

        # numerical factor to account for the assumption that resonance density in a N-body system is K times larger than in the 3-planet case
        # scales up 3-planet system to N-planet system
        Naffect_val = max(i_out - i_in + 1, 3)  # eq 5
        Naffect[i] = Naffect_val

        # chooses interacting pair from the triplet
        d1 = (1.0 - ecc[i + 1]) * ap[i + 1] - (1.0 + ecc[i]) * ap[
            i
        ]  # eq 4, closest physical distance so not normalised
        d2 = (1.0 - ecc[i + 2]) * ap[i + 2] - (1.0 + ecc[i + 1]) * ap[i + 1]
        pair_index = i if d1 <= d2 else i + 1
        interacting_pair[i] = pair_index

        # eccentricity at overlap eq6
        ecc_cross = [
            np.sqrt(Mp[pair_index + 1])
            * abs(ap[pair_index + 1] - ap[pair_index])
            / (
                np.sqrt(Mp[pair_index + 1]) * ap[pair_index]
                + np.sqrt(Mp[pair_index]) * ap[pair_index + 1]
            ),
            np.sqrt(Mp[pair_index])
            * abs(ap[pair_index + 1] - ap[pair_index])
            / (
                np.sqrt(Mp[pair_index + 1]) * ap[pair_index]
                + np.sqrt(Mp[pair_index]) * ap[pair_index + 1]
            ),
        ]

        # as before, the eccentricity can't be lower than eccentricity actually required for an overlap
        ecc_cross[0] = max(ecc[pair_index], ecc_cross[0])  # eq 23
        ecc_cross[1] = max(ecc[pair_index + 1], ecc_cross[1])

        # calculate interaction timescales for every triplet, calling timescale functions from above
        viscous_timescale = tau_vis(
            ap[pair_index : pair_index + 2],
            Mp[pair_index : pair_index + 2],
            Rp[pair_index : pair_index + 2],
            Ms,
            ecc_cross,
        )
        collision_timescale = tau_col(
            ap[pair_index : pair_index + 2],
            Mp[pair_index : pair_index + 2],
            Rp[pair_index : pair_index + 2],
            Ms,
            ecc_cross,
        )

        # correction with secular perturbations, eq 21
        ecc_dmy = [
            np.sqrt(np.sum(ecc_vec[i, :] ** 2)),
            np.sqrt(np.sum(ecc_vec[i + 1, :] ** 2)),
            np.sqrt(np.sum(ecc_vec[i + 2, :] ** 2)),
        ]

        # passes the TRIPLETS ap, Mp, ecc to the wrapper function to calculate crossing timescale
        crossing_timescale = interaction_wrapper(
            ap[i : i + 3], Mp[i : i + 3], Ms, ecc_dmy, Naffect_val
        )
        # predicted crossing timescale is 'current' time t + predicted next interaction + duration of the interaction
        Tcross[i] = t + crossing_timescale + min(viscous_timescale, collision_timescale)

    idmy_min = np.argmin(Tcross)  # index of minimum crossing time
    t_event = Tcross[idmy_min]
    icross = interacting_pair[idmy_min]

    return icross, t_event  # inner planet index, event time

def interaction_wrapper(
    ap, Mp, Ms, ecc, N_affect
):  # determine if system is stable, and if not, calculate timescale to instability
    """
    Function to determine if triplet is stable, and if not, calculate timescale of instability

    Parameters
    ----------
    ap : list
        Semi-major axes of triplet [m]
    Mp : list
        Masses of triplet [kg]
    Ms : float
        Stellar mass [kg]
    ecc : list
        Eccentricities of triplet
    N_affect : int
        Number of planets interacting

    Returns
    -------
    float
        1e20 if the system is stable, otherwise the instability timescale
        from ``tau_cross_petit`` [s]
    """
    # N_affect is planets participating in crossing event
    aM = (Mp[0] * ap[0] + Mp[1] * ap[1]) / (Mp[0] + Mp[1])
    h = hill_sphere(aM, Mp[0] + Mp[1], Ms) / aM
    # stability criterion
    EJbef = (
        5.0 / 8.0 * (ecc[0] ** 2 + ecc[1] ** 2) / h**2
        - 3.0 / 8.0 * ((ap[0] - ap[1]) / (h * aM)) ** 2
        + 4.5
    )  # eq 28

    if (
        N_affect <= 2 and EJbef < 0.0
    ):  # if system is stable return 'infinite' stability, 2 planets can still interact
        return 1e20

    # otherwise return the crossing timescale from Petit 2020
    return tau_cross_petit(ap, Mp, Ms, ecc, N_affect)

def tau_col(ap, Mp, Rp, Ms, ecc):
    """
    Function to calculate duration of a collision event

    Parameters
    ----------
    ap : list
        Semi-major axes of interacting pair [m]
    Mp : list
        Masses of interacting pair [kg]
    Rp : list
        Radii of interacting pair [m]
    Ms : float
        Stellar mass [kg]
    ecc : list
        Eccentricities of interacting pair

    Returns
    -------
    timescale : float
        Duration of collisional event [s]
    """
    mu_a = sum(ap) / 2  # average semi-major axis of interacting pair
    M_T = sum(Mp)  # sum of masses
    R_T = sum(Rp)  # sum of radii
    impact_parameter = abs(ap[1] - ap[0])

    # eccentricities at the onset of crossing
    ecross_i = (np.sqrt(Mp[1]) * impact_parameter) / (
        (np.sqrt(Mp[1]) * ap[0]) + np.sqrt(Mp[0]) * ap[1]
    )  # eq 6
    ecross_j = (np.sqrt(Mp[0]) * impact_parameter) / (
        (np.sqrt(Mp[1]) * ap[0]) + np.sqrt(Mp[0]) * ap[1]
    )  # eq 6, same denominator for both planets
    ecross = [ecross_i, ecross_j]

    rep_e = 0.5 * max(sum(ecross), sum(ecc))
    kep_vel = np.sqrt((G * Ms) / mu_a)
    ran_vel = rep_e * kep_vel
    esc_vel = np.sqrt(2.0 * G * M_T / R_T)
    n = 1.0 / (2.0 * np.pi * rep_e * mu_a**2 * impact_parameter)

    # only difference from tau_vis here is the timescale
    timescale = n * np.pi * (R_T) ** 2 * (1 + esc_vel**2 / ran_vel**2) * ran_vel  # eq 11

    return 1 / timescale

def tau_cross_petit(
    a, Mp, Ms, ecc, N_affect
):  # evaluates every planetary triplet for instability
    """
    Function to calculate timescale to instability as a result of 3-body mean motion resonances

    Parameters
    ----------
    a : list
        Semi-major axes of planetary triplet [m]
    Mp : list
        Masses of planetary triplet [kg]
    Ms : float
        Stellar mass [kg]
    ecc : list
        Eccentricities of planetary triplet
    N_affect : int
        Number of planets interacting

    Returns
    -------
    tau_cross : float
        Instability timescale for planetary triplet [s]
    """
    K = min(0.5 * (N_affect - 3) + 1, 3)
    alpha_01, alpha_12 = a[0] / a[1], a[1] / a[2]

    nu_01, nu_12 = alpha_01**1.5, alpha_12**1.5  # different from paper appendix B

    eta = (nu_01 * (1 - nu_12)) / (1 - nu_01 * nu_12)
    M = (
        np.sqrt(
            Mp[0] * Mp[2]
            + Mp[1] * Mp[2] * eta**2 * alpha_01 ** (-2)
            + Mp[0] * Mp[1] * alpha_12**2 * (1 - eta) ** 2
        )
        / Ms
    )
    delta_01, delta_12 = (
        1.0 - ecc[1] - (1.0 + ecc[0]) * alpha_01,
        1.0 - ecc[2] - (1.0 + ecc[1]) * alpha_12,
    )  # also different from paper, but noted in the code

    if delta_01 < 0.0 or delta_12 < 0.0:
        return 0.0

    delta = (delta_01 * delta_12) / (delta_01 + delta_12)
    delta_ov = (6.55 * K * M) ** (1 / 4) * (eta * (1 - eta)) ** (3 / 8)

    if delta >= delta_ov:  # tau cross approaches infinity
        return 1e20

    log_arg = (
        -np.log10((32 * np.sqrt(19) * M * np.sqrt(eta * (1 - eta))) / (3 * np.sqrt(np.pi)))
        + np.log10(delta**6 / (delta_ov**6 * (1 - (delta / delta_ov) ** 4)))
        + np.sqrt(-np.log(1 - (delta / delta_ov) ** 4))
    )
    tau_cross = 10 ** (log_arg) * kepler_period(Mp[0], Ms, a[0])

    return tau_cross

def tau_vis(
    ap, Mp, Rp, Ms, ecc
):  # viscous relaxation timescale for an interacting planetary pair
    """
    Viscous relaxation timescale for an interacting planetary pair

    Parameters
    ----------
    ap : list
        Semi-major axes of interacting pair [m]
    Mp : list
        Masses of interacting pair [kg]
    Rp : list
        Radii of interacting pair [m]
    Ms : float
        Stellar mass [kg]
    ecc : list
        Eccentricities of interacting pair

    Returns
    -------
    timescale : float
        Viscous relaxation timescale (when system settles down after scattering)
    """
    mu_a = sum(ap) / 2  # average semi-major axis of interacting pair
    M_T = sum(Mp)  # sum of masses
    impact_parameter = abs(ap[1] - ap[0])

    # eccentricities at the onset of crossing
    ecross_i = (np.sqrt(Mp[1]) * impact_parameter) / (
        (np.sqrt(Mp[1]) * ap[0]) + np.sqrt(Mp[0]) * ap[1]
    )  # eq 6
    ecross_j = (np.sqrt(Mp[0]) * impact_parameter) / (
        (np.sqrt(Mp[1]) * ap[0]) + np.sqrt(Mp[0]) * ap[1]
    )  # eq 6, same denominator for both planets
    ecross = [ecross_i, ecross_j]

    rep_e = 0.5 * max((sum(ecross)), sum(ecc))  # eq 23 used to calculate lambda in eq 12
    kep_vel = np.sqrt((G * Ms) / mu_a)  # eq 10
    ran_vel = rep_e * kep_vel  # before eq 9
    n = 1 / (2 * np.pi * rep_e * mu_a**2 * impact_parameter)  # eq 8

    timescale = ran_vel**3 / (n * np.pi * G**2 * M_T**2)  # matches Fortran (no factor of 3)
    return timescale

def secular_solution(ap, Mp, ecc, Rp, Ms, N):
    """
    Solution to secular perturbation theory of orbital mechanics. Only recalculates after an event

    Parameters
    ----------
    ap : list
        Semi-major axes of all planets in system [m]
    Mp : list
        Masses of all planets in system [kg]
    ecc : list
        Eccentricities of all planets in system
    Rp : list
        Radii of all planets in system [m]
    Ms : float
        Stellar mass [kg]
    N : int
        Number of planets

    Returns
    -------
    ecc_vec : array
        Scaled eigenvectors - columns represent eccentricity contribution from each mode per planet
    g : array
        Secular eigenfrequencies
    beta : array
        Phase angle for each planet in system

    """
    varpi = np.random.uniform(0.0, 2.0 * np.pi, N)

    h0 = ecc * np.sin(varpi)
    k0 = ecc * np.cos(varpi)

    mean_motion = np.sqrt(G * Ms / ap**3)
    A = np.zeros((N, N))  # empty interaction matrix
    apo = (1.0 + ecc) * ap
    peri = (1.0 - ecc) * ap

    laplace = LaplaceCoefficient(method='Brute')  # to calculate laplace coefficients
    for i in range(N):
        for j in range(N):
            if i == j:
                continue  # skip self-interactions
            if min(apo[i], apo[j]) > max(peri[i], peri[j]):
                continue  # skip overlapping orbits
            if ap[i] < ap[j]:
                alpha = ap[i] / ap[j]
                alpha_bar = alpha
            else:
                alpha = ap[j] / ap[i]
                alpha_bar = 1

            # calculate laplacian coefficients, b1,3/2 and b2,3/2
            # (a,s,m,p,q)
            coeff_m1 = laplace(alpha, 3 / 2, 1, 1, 1)  # m = 1 for A_ii
            coeff_m2 = laplace(alpha, 3 / 2, 2, 1, 1)  # m = 2 for A_ij

            factor = mean_motion[i] * 0.25 * Mp[j] / (Ms + Mp[i]) * alpha * alpha_bar
            A[i, i] += factor * coeff_m1
            A[i, j] = -factor * coeff_m2

    # solve the eigenvalue problem
    g, S = np.linalg.eig(A)

    g = np.real(g)
    S = np.real(S)

    # solve for integration constants using S
    Csinb = np.linalg.solve(S, h0)
    Ccosb = np.linalg.solve(S, k0)

    # calculate amplitudes (C) and phase angles (beta)
    C = np.sqrt(Csinb**2 + Ccosb**2)
    beta = np.arctan2(Csinb, Ccosb)

    # scale eigenvectors by the amplitudes (columns of ecc_vec)
    ecc_vec = S * C

    return ecc_vec, g, beta

def orbit_cross_K25(
    ap, Mp, Rp, Ms, impact_parameter, ecc, interact, live_status, N, planet_id, icross
):  # determine outcome of crossing event
    """
    Function to determine what happens once an orbit crossing event has been established

    Parameters
    ----------
    ap : list
        Semi major axes of interacting pair [m]
    Mp : list
        Masses of interacting pair [kg]
    Rp : list
        Radii of interacting pair [m]
    Ms : float
        Stellar mass [kg]
    impact_parameter: float
        Defined as sin(impact angle). Describes angle of contact between target and impactor
    ecc : list
        Eccentricities of interacing pair
    interact : bool
        Interaction status of interacting pair (set to True)
    live_status : bool
        Live status of interacting pair (starts as True and will be modified depending on scattering vs. merge events)
    N : int
        Number of planets
    planet_id : list
        Persistent ID number for interacting planets to track evolutions
    icross: int
        Index of inner planet in interacting pair

    Returns
    -------
    Modified arrays for ap, Mp, Rp, ecc, interact, live_status, N following scattering or merging
    merge_record : dict
        Stores information about targets and impactors for each merging event
    """

    # now working with an interacting pair of planets i,j
    # modifies the arrays of ap,Mp,Rp,ecc,interact,live_status based on what happens
    jcross = (
        icross + 1
    )  # sets indices of interacting pair, aj>ai always, +1 to be able to index a pair later

    # mass-weighted mean semi-major axis of the pair (used for EJbef stability condition and the ejection formula)
    aM = (Mp[icross] * ap[icross] + Mp[jcross] * ap[jcross]) / (Mp[icross] + Mp[jcross])
    h = hill_sphere(aM, Mp[icross] + Mp[jcross], Ms) / aM
    EJbef = (
        5.0 / 8.0 * (ecc[icross] ** 2 + ecc[jcross] ** 2) / h**2
        - 3.0 / 8.0 * ((ap[icross] - ap[jcross]) / (h * aM)) ** 2
        + 4.5
    )

    # 2-body case: if the pair is Jacobi-stable and not currently overlapping, nothing happens
    if N == 2 and (
        EJbef < 0.0 and (1.0 + ecc[icross]) * ap[icross] < (1.0 - ecc[jcross]) * ap[jcross]
    ):
        return

    mean_ap = (ap[icross] + ap[jcross]) / 2  # simple average semi-major axis, used for e_esc
    e_esc = esc_ecc(
        Ms, Mp[icross], Mp[jcross], Rp[icross], Rp[jcross], mean_ap
    )  # escape eccentricity

    ecross_i = (
        (ap[jcross] - ap[icross])
        * np.sqrt(Mp[jcross])
        / (ap[icross] * np.sqrt(Mp[jcross]) + ap[jcross] * np.sqrt(Mp[icross]))
    )  # eq 6 again
    ecross_j = (
        (ap[jcross] - ap[icross])
        * np.sqrt(Mp[icross])
        / (ap[icross] * np.sqrt(Mp[jcross]) + ap[jcross] * np.sqrt(Mp[icross]))
    )

    ecc_cross = [ecross_i, ecross_j]  # store together
    # minimum eccentricity of interacting planets to cause a collision

    ecc_encounter = [max(ecc_cross[0], ecc[icross]), max(ecc_cross[1], ecc[jcross])]  # eq 23
    # enforces that the planets ACTUALLY interact (actual eccentricites from ecc[] are secular, not at exact crossing time)
    # crossing_pair is still just a statistical prediction
    # ensures the eccentricities used during a crossing are geometrically consistent
    eij = np.sqrt(ecc_encounter[0] ** 2 + ecc_encounter[1] ** 2)  # relative eccentricity

    # calculates collosion probability, Pcol
    ln_lambda = 3
    lambdaa = (2 * eij / e_esc) ** 2 * (1 + (eij**2 / e_esc**2)) * (1 / ln_lambda)  # eq 12
    p_col = 1 - np.exp(-lambdaa)  # eq 14, caps between [0,1]

    # use Monte Carlo approach to say whether or not a collision actually happens
    # Bernoulli sampling?
    merge_record = None  # stores when a merge happens, stays None if nothing happens
    draw = np.random.uniform(0, 1)  # draws a random number to compare against p_col
    if draw < p_col:  # merge event
        count = 0  # keeps track of rejection-sampling structure, and ensures while loop doesnt go forever
        while True:  # keeps sampling until the eccentricity is consistent with the orbits actually overlapping
            rayleigh_ecc = rayleigh(1 / np.sqrt(2), eij / e_esc)
            # mutate ecc[icross]/ecc[jcross] in place so each retry's max() compares against the
            # running (already-updated) eccentricity, matching the Fortran goto-loop behaviour
            ecc[icross] = max(
                rayleigh_ecc * np.sqrt(Mp[jcross]) / np.sqrt(Mp[icross] + Mp[jcross]) * e_esc,
                ecc[icross],
            )
            ecc[jcross] = max(
                rayleigh_ecc * np.sqrt(Mp[icross]) / np.sqrt(Mp[icross] + Mp[jcross]) * e_esc,
                ecc[jcross],
            )

            if np.sqrt(ecc_cross[0] ** 2 + ecc_cross[1] ** 2) / e_esc > 2.0 or count > 500:
                # unable to find a random draw that satisfies the condition, defaul to an orbital overlap of 0.1%
                # if orbits will never overlap OR took too many interations.
                # Take the larger of the geometric default and the eccentricity
                # the body already carries, the same convention as every other
                # site: overwriting outright discards the secular value, which
                # can be the larger of the two and then reports a collision
                # slower than the incoming orbits imply.
                ecc[icross] = max(1.001 * ecc_cross[0], ecc[icross])
                ecc[jcross] = max(1.001 * ecc_cross[1], ecc[jcross])
                break
            # check if epicycle amplitudes sum to at least aj-ai (they will overlap)
            # OVERLAP CONDITION
            elif ap[icross] * ecc[icross] + ap[jcross] * ecc[jcross] >= abs(
                ap[jcross] - ap[icross]
            ):
                # yay it worked, eccentricities are already updated above
                break
            else:  # if all else, try again with a new random number
                count += 1

        # velocity at collision, calculated from energy conservation
        v_c = collision_velocity(
            ap[icross : jcross + 1],
            Mp[icross : jcross + 1],
            Rp[icross : jcross + 1],
            Ms,
            ecc[icross : jcross + 1],
        )

        # capture pre-merge state for the impact log, since merge_embryo overwrites Mp[icross:jcross+1] below
        # (target = larger/surviving body, impactor = smaller/destroyed body, as per merge_embryo convention)
        target_idx = icross if Mp[icross] >= Mp[jcross] else jcross
        impactor_idx = jcross if target_idx == icross else icross
        id_target, id_impactor = (
            planet_id[target_idx],
            planet_id[impactor_idx],
        )  # carry index through all other parameters
        M_target_before, M_impactor_before = Mp[target_idx], Mp[impactor_idx]
        # pre-merge geometry, captured before merge_embryo overwrites ap below; these
        # let a consumer rebuild the impact kinematics (escape velocity, densities)
        a_target_before = ap[target_idx]  # [m]
        e_target_before = ecc[target_idx]  # [1], before merge_embryo overwrites it
        R_target_before, R_impactor = Rp[target_idx], Rp[impactor_idx]  # [m]

        # call merge_embryo function to update parameters for interacting pair
        # jcross+1 to include that planet in the interacting pair
        # print(f"[COLLISION] Planets {planet_id[icross]} and {planet_id[jcross]} merged")
        ap_merge, Mp_merge, ecc_merge, live_status_merge = merge_embryo(
            ap[icross : jcross + 1],
            Mp[icross : jcross + 1],
            Rp[icross : jcross + 1],
            Ms,
            ecc[icross : jcross + 1],
            v_c,
            live_status[icross : jcross + 1],
            impact_parameter,
        )
        # update system
        ap[icross : jcross + 1] = ap_merge
        Mp[icross : jcross + 1] = Mp_merge
        ecc[icross : jcross + 1] = ecc_merge
        live_status[icross : jcross + 1] = live_status_merge  # smaller planet dies

        # record impact velocity (v_c). The pre-merge radii, target semi-major axis
        # and target eccentricity, and the post-merge semi-major axis and eccentricity,
        # are kept alongside so the full impact geometry can be reconstructed later to calculate atmospheric mass loss.
        # Both orbital elements are reported before and after, so a consumer can apply the
        # change this collision made rather than the absolute value, which matters when the
        # consumer follows a planet on a different orbit from this body's.
        merge_record = {
            'id_target': id_target,
            'id_impactor': id_impactor,
            'M_target_before': M_target_before,
            'M_impactor_before': M_impactor_before,
            'M_merged_after': Mp[target_idx],
            'v_c': v_c,
            'a_final_AU': ap[target_idx] / au2m,
            'R_target_before': R_target_before,
            'R_impactor': R_impactor,
            'a_before': a_target_before,
            'e_before': e_target_before,
            'e_after': ecc[target_idx],
        }

    else:  # scattering event
        rayleigh_ecc = rayleigh(
            1.0 / np.sqrt(2.0), 0.0
        )  # 0 because no truncation at geometric overlap constraint as for merge case
        # assuming energy equipartition
        # taken from fortran code, why is it not a max() anymore as before with merging?

        ecc[icross] = (
            rayleigh_ecc * np.sqrt(Mp[jcross]) / np.sqrt(Mp[icross] + Mp[jcross]) * e_esc
        )
        ecc[jcross] = (
            rayleigh_ecc * np.sqrt(Mp[icross]) / np.sqrt(Mp[icross] + Mp[jcross]) * e_esc
        )

        # Identify which planet was excited more. Equipartition sets the two
        # excitations in the fixed ratio sqrt(M_other/M_self), so for an
        # equal-mass pair they are identical to the last bit and the choice is
        # a tie that has to be made deliberately: the paper's prescription does
        # not cover it, and it decides whether such an ejection removes one body
        # or two. The outer body escapes, which leaves the survivor inside the
        # escaper, where the apocentre closure below gives a bound orbit.
        if ecc[icross] == ecc[jcross]:
            ilarge, ismall = jcross, icross
        else:
            ilarge = icross if ecc[icross] > ecc[jcross] else jcross
            ismall = jcross if ecc[icross] > ecc[jcross] else icross

        if max(ecc[icross], ecc[jcross]) >= 1.0:  # planet got bumped out
            # print(f"[EJECTION] Planet {planet_id[ilarge]} was ejected")
            # planet with smaller excited eccentricity remains in the system, and orbital parameters are recalculated
            # The survivor absorbs the escaper's binding energy, so its orbit
            # tightens. It must still pass through the place the encounter
            # happened, and since that place is now outside its new orbit it
            # becomes the apocentre: a_new (1 + e_new) = a_old. Capture the old
            # axis before overwriting it.
            a_rem_before = ap[ismall]
            ap[ismall] = Mp[ismall] / (Mp[ismall] / ap[ismall] + Mp[ilarge] / ap[ilarge])
            ecc[ismall] = a_rem_before / ap[ismall] - 1.0
            live_status[ilarge] = False  # the body excited past e = 1 is the one that escapes
            if ecc[ismall] >= 1.0:
                # A comparable-mass pair leaves no bound survivor: the orbit that
                # closes the energy books is too tight to still reach the
                # encounter, so the remaining body is unbound too and the
                # encounter destroys both.
                live_status[ismall] = False
        else:  #'normal' scattering conditions
            #'change in orbital separation is assumed to be equal to the sum of the excited epicycle amplitude'
            # db = delta_a essentially how much the orbit is shifted either in or out
            # print(f"[SCATTERING] Planets {planet_id[icross]} and {planet_id[jcross]} scattered")
            db = ecc[icross] * ap[icross] + ecc[jcross] * ap[jcross]  # delta b_ij eq 18
            ap[icross] = ap[icross] - Mp[jcross] / (Mp[icross] + Mp[jcross]) * db  # eq 19
            ap[jcross] = ap[jcross] + Mp[icross] / (Mp[icross] + Mp[jcross]) * db  # eq 20

    # a negative semi-major axis is a hyperbolic orbit: the scattering unbound
    # the inner body, so it leaves the system rather than falling inward
    if ap[icross] < 0.0:
        live_status[icross] = False  # it's now dead

    return merge_record

def collision_velocity(ap, Mp, Rp, Ms, ecc):
    """
    Velocity of collision between target and impactor during merge events

    Parameters:
    ap : list
        Semi-major axes of interacting pair [m]
    Mp : list
        Masses of interacting pair [kg]
    Rp : list
        Radii of interacting pair [m]
    Ms : float
        Stellar mass [kg]
    ecc : list
        Eccentricities of interacting pair

    Returns
    -------
    v_c : float
        Collision velocity between target and impactor during merging [m/s]
    """
    mu_a = sum(ap) / 2
    kep_vel = np.sqrt(G * Ms / mu_a)
    rep_e = np.sqrt(ecc[0] ** 2 + ecc[1] ** 2)  # same as eij in orbit_cross_K25
    v_inf = rep_e * kep_vel
    v_esc = np.sqrt(2 * G * (Mp[0] + Mp[1]) / (Rp[0] + Rp[1]))
    v_c = np.sqrt(v_inf**2 + v_esc**2)
    return v_c

def merge_embryo(
    ap, Mp, Rp, Ms, ecc, v_c, live_status, b
):  # calculate orbital parameters post collision
    """
    Function to compute mass, orbital separation, and eccentricity after a giant impact

    Parameters
    ----------
    ap : list
        Semi major axes of colliding pair [m]
    Mp : list
        Masses of colliding pair [kg]
    Rp : list
        Radii of colliding pair [m]
    ecc : list
        Secular eccentricity of colliding pair just before collision
    v_c : float
        Collision velocity [m/s]
    live_status : list
        Live_status of colliding pair
    b : float
        Impact parameter, defined as sin(beta) where beta is impact angle in [deg]

    Returns
    -------
    ap : float
        Updated semi-major axis of surviving planet [m]
    Mp : float
        Updated mass of surviving planet [kg]
    ecc : float
        Updated eccentricity of surviving planet
    live_status : list
        Sets consumed planet's status to 'False'

    """
    # Mp is a list containing interacting pair masses
    Mp_new = sum(Mp)  # eq 15
    ap_new = Mp_new / (Mp[0] / ap[0] + Mp[1] / ap[1])  # eq 16

    if ap[1] * (1.0 - ecc[1] ** 2) < ap[0] * (1.0 - ecc[0] ** 2):
        min_dvarpi = 0.0
    elif ecc[0] == 0.0 or ecc[1] == 0.0:
        # a circular orbit has no pericentre direction, so no alignment is
        # forbidden and the draw is unrestricted. Handled explicitly because the
        # general expression divides by the product of the two eccentricities;
        # it happens to recover the same answer through an infinity and the
        # clamp below, but only by accident, and it warns on the way.
        min_dvarpi = 0.0
    else:
        cosdvarpi = ((ecc[0] * ap[0]) ** 2 + (ecc[1] * ap[1]) ** 2 - (ap[1] - ap[0]) ** 2) / (
            2.0 * ecc[0] * ecc[1] * ap[0] * ap[1]
        )
        cosdvarpi = max(-1.0, min(1.0, cosdvarpi))
        min_dvarpi = np.arccos(cosdvarpi)

    dvarpi = np.random.uniform(min_dvarpi, 2.0 * np.pi - min_dvarpi)
    ecc_new = np.sqrt(
        (
            (Mp[0] ** 2 * ecc[0] ** 2)
            + (Mp[1] ** 2 * ecc[1] ** 2)
            + 2 * Mp[0] * Mp[1] * ecc[0] * ecc[1] * np.cos(dvarpi)
        )
        / Mp_new**2
    )  # eq 17

    if Mp[0] >= Mp[1]:  # larger planet consumes smaller one
        alive, dead = 0, 1
    else:
        alive, dead = 1, 0

    live_status[dead] = False  # kill consumed planet
    ap[alive] = ap_new
    Mp[alive] = Mp_new
    ecc[alive] = ecc_new

    return ap, Mp, ecc, live_status