3I ATLAS trajectory Bayes factor Calculation

3iatlas_combined_odds.py

import math
from dataclasses import dataclass


def calculate_unbiased_odds() -> None:
    """
    Calculates the combined odds of orbital anomalies summing probabilities
    across all valid candidates to avoid p-hacking.
    """

    @dataclass
    class Planet:
        name: str
        a: float  # Semi-major axis in AU

    # Solar System Data
    planets = [
        Planet("Mercury", 0.387098),
        Planet("Venus", 0.728213),
        Planet("Earth", 1.000),
        Planet("Mars", 1.52368055),
        Planet("Jupiter", 5.2038),
        Planet("Saturn", 9.5826),
        Planet("Uranus", 19.19126),
        Planet("Neptune", 30.07),
        Planet("Pluto", 39.482),
    ]

    # Interstellar Object Parameters
    Q_PERIHELION = 1.36

    def get_encounter_prob(planet: Planet, limit_au: float) -> float:
        """Calculates probability of encounter within limit_au."""
        circumference = 2 * math.pi * planet.a

        # Case 1: Crossing (Planet orbit is larger than object perihelion)
        if planet.a > Q_PERIHELION:
            target_width = 2 * limit_au
            # Two crossings (inbound/outbound)
            return 2 * (target_width / circumference)

        # Case 2: Grazing (Planet orbit is inside object perihelion)
        else:
            min_dist = abs(Q_PERIHELION - planet.a)
            if min_dist < limit_au:
                # Chord length calculation
                # L = 2 * sqrt(R^2 - d^2)
                chord = 2 * math.sqrt(limit_au**2 - min_dist**2)
                return chord / circumference
            else:
                return 0.0

    print(f"{'Condition':<40} | {'Prob %':<10}")
    print("-" * 55)

    # 1. Inner Planet Condition (< 0.25 AU)
    # Candidates: Mercury, Venus, Earth, Mars
    inner_candidates = planets[:4]
    prob_inner = sum(get_encounter_prob(p, 0.25) for p in inner_candidates)
    print(f"Any Inner Planet (<0.25 AU)             | {prob_inner * 100:.2f}%")

    # 2. Outer Planet Condition (< 0.5 AU)
    # Candidates: Jupiter through Pluto
    outer_candidates = planets[4:]
    prob_outer = sum(get_encounter_prob(p, 0.50) for p in outer_candidates)
    print(f"Any Outer Planet (<0.50 AU)             | {prob_outer * 100:.2f}%")

    # 3. Any Other Planet Condition (< 0.75 AU)
    # Candidates: All 9
    # Note: This assumes the events are statistically independent or
    # we are looking for the compound probability of matching all 3 conditions.
    prob_other = sum(get_encounter_prob(p, 0.75) for p in planets)
    print(f"Any Planet (<0.75 AU)                   | {prob_other * 100:.2f}%")

    # 4. Angular Constraints
    # Plane < +-5 degrees = 10 deg (Isotropic)
    prob_plane = 2 * (1 - math.cos(math.radians(10)))

    # Alignment with ANY Inner Planet (4 targets) (+-5 deg = 10 deg)
    # 4 targets * (10 degrees / 360 degrees)
    prob_align = 4 * (10 / 360)

    print(f"Orbital Plane (<5 deg)                  | {prob_plane * 100:.3f}%")
    print(f"Perihelion Align (Any Inner Planet)     | {prob_align * 100:.2f}%")

    # Total Calculation
    total_p = prob_inner * prob_outer * prob_other * prob_plane * prob_align

    print("-" * 55)
    print(f"Combined Probability given random: {total_p:.4e}")
    if total_p > 0:
        print(f"Null Odds: 1 in {1 / total_p:,.0f}")

    # ========================== Alternate explanation: Aliens ========================
    print(f"{'Condition':<40} | {'Prob %':<10}")
    print("-" * 55)

    # Voyager 2 encounters:
    # Neptune    4950 km
    # Galatea   18360 km
    # Larissa   60180 km
    # Proteus   97860 km
    # Triton    39800 km
    # Uranus   107000 km
    # Miranda   29000 km
    # Ariel    127000 km
    # Saturn   161000 km
    # Calypso  151590 km
    # Pandora  107000 km
    # Enceladus 87010 km
    # Tethys    93010 km
    # Jupiter  721670 km
    # Ganymede  62130 km
    # Clearly 100,000 km is a reasonable upper bound target = 0.0006684587 AU

    # 1. Inner Planet Condition (< 100k km)
    # Candidates: Mercury, Venus, Earth, Mars
    inner_candidates = planets[:4]
    prob_inner2 = sum(get_encounter_prob(p, 0.0006684587) for p in inner_candidates)
    print(f"Any Inner Planet (<100k km)             | {prob_inner2 * 100:.2f}%")

    # 2. Outer Planet Condition (< 100k km)
    # Candidates: Jupiter through Pluto
    outer_candidates = planets[4:]
    prob_outer2 = sum(get_encounter_prob(p, 0.0006684587) for p in outer_candidates)
    print(f"Any Outer Planet (<100k km)             | {prob_outer2 * 100:.2f}%")

    # 3. Any Other Planet Condition (< 100k km)
    # Candidates: All 9
    # Note: This assumes the events are statistically independent or
    # we are looking for the compound probability of matching all 3 conditions.
    prob_other2 = sum(get_encounter_prob(p, 0.0006684587) for p in planets)
    print(f"Any Other Planet (<100k km)             | {prob_other2 * 100:.2f}%")

    # 4. Angular Constraints
    # Plane < +-1 degree = 2 deg (Isotropic)
    prob_plane2 = 2 * (1 - math.cos(math.radians(2)))

    # Alignment with ANY Inner Planet (4 targets) (+-1 deg = 2 deg)
    # 4 targets * (2 degrees / 360 degrees)
    prob_align2 = 4 * (2 / 360)

    print(f"Orbital Plane (<1 deg)                  | {prob_plane2 * 100:.3f}%")
    print(f"Perihelion Align (Any Inner Planet)     | {prob_align2 * 100:.2f}%")

    # Total Calculation
    total_p2 = prob_inner2 * prob_outer2 * prob_other2 * prob_plane2 * prob_align2

    print("-" * 55)
    print(f"Combined Probability given aliens: {total_p2:.4e}")
    if total_p > 0:
        print(f"Alien Odds: 1 in {1 / total_p2:,.0f}")

    print(f"Odds ratio: {total_p2 / total_p:.4e}")
    print(f"Bayes factor: {(total_p2 / total_p) / total_p}")


if __name__ == "__main__":
    calculate_unbiased_odds()