3i/Atlas Closest

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Orbital elements for 3I/ATLAS (from search results)
# Using JPL solution from July 2025
atlas_e = 6.143  # eccentricity
atlas_q = 1.357  # perihelion distance (AU)
atlas_i = 175.1  # inclination (degrees) - retrograde
atlas_omega = 0  # argument of perihelion (approximated)
atlas_Omega = 0  # longitude of ascending node (approximated)
atlas_v_inf = 58  # velocity at infinity (km/s)

# New Horizons trajectory data (from search results)
# Post-Jupiter trajectory characteristics
nh_inclination = 2.3  # degrees (after Jupiter flyby)
nh_eccentricity = 1.40  # hyperbolic orbit
nh_perihelion = 2.2  # AU

# Convert angles to radians
atlas_i_rad = np.radians(atlas_i)
atlas_omega_rad = np.radians(atlas_omega)
atlas_Omega_rad = np.radians(atlas_Omega)
nh_i_rad = np.radians(nh_inclination)


def calculate_position(a, e, i, omega, Omega, true_anomaly, mu=1.0):
    """
    Calculate 3D position from orbital elements
    Using simplified two-body problem
    """
    # Distance from focus
    r = a * (1 - e**2) / (1 + e * np.cos(true_anomaly))

    # Position in orbital plane
    x_orb = r * np.cos(true_anomaly)
    y_orb = r * np.sin(true_anomaly)
    z_orb = 0

    # Rotation matrices to transform to heliocentric ecliptic frame
    # Rotate by argument of periapsis
    cos_omega = np.cos(omega)
    sin_omega = np.sin(omega)

    # Rotate by inclination
    cos_i = np.cos(i)
    sin_i = np.sin(i)

    # Rotate by longitude of ascending node
    cos_Omega = np.cos(Omega)
    sin_Omega = np.sin(Omega)

    # Combined rotation matrix
    x = x_orb * (cos_omega * cos_Omega - sin_omega * sin_Omega * cos_i) - y_orb * (
        sin_omega * cos_Omega + cos_omega * sin_Omega * cos_i
    )

    y = x_orb * (cos_omega * sin_Omega + sin_omega * cos_Omega * cos_i) + y_orb * (
        cos_omega * cos_Omega * cos_i - sin_omega * sin_Omega
    )

    z = x_orb * sin_omega * sin_i + y_orb * cos_omega * sin_i

    return np.array([x, y, z])


def hyperbolic_trajectory_points(e, q, i, omega, Omega, n_points=100):
    """Generate points along hyperbolic trajectory"""
    a = q / (1 - e)  # semi-major axis (negative for hyperbola)

    # For hyperbolic orbits, true anomaly is limited
    max_true_anomaly = np.arccos(-1 / e)  # Maximum physical true anomaly
    true_anomalies = np.linspace(
        -max_true_anomaly * 0.99, max_true_anomaly * 0.99, n_points
    )

    positions = []
    for ta in true_anomalies:
        pos = calculate_position(abs(a), e, i, omega, Omega, ta)
        positions.append(pos)

    return np.array(positions)


# Generate trajectory points for 3I/ATLAS
print("Calculating 3I/ATLAS trajectory...")
atlas_a = atlas_q / (1 - atlas_e)  # semi-major axis
atlas_positions = hyperbolic_trajectory_points(
    atlas_e, atlas_q, atlas_i_rad, atlas_omega_rad, atlas_Omega_rad
)

# For New Horizons, approximate as a straight line trajectory
# Since it's been traveling for ~20 years since launch
print("Calculating New Horizons trajectory...")

# New Horizons approximate position progression
# Started near Earth, passed Jupiter in 2007, now ~62 AU out
# Simplified trajectory: starts near 1 AU, goes to ~62+ AU

# Time points (years since 2006 launch)
time_years = np.linspace(0, 25, 100)
nh_distances = 1 + time_years * 2.5  # Rough approximation: ~2.5 AU per year average

# New Horizons travels roughly in ecliptic plane with small inclination
nh_positions = []
for i, dist in enumerate(nh_distances):
    # Approximate trajectory direction (roughly toward constellation Sagittarius)
    # Using simplified coordinates
    x = dist * np.cos(np.radians(270))  # Roughly toward galactic center
    y = dist * np.sin(np.radians(270))
    z = dist * np.sin(nh_i_rad) * 0.1  # Small inclination effect
    nh_positions.append([x, y, z])

nh_positions = np.array(nh_positions)

# Calculate closest approach
print("Finding closest approach...")
min_distance = float("inf")
closest_atlas_idx = 0
closest_nh_idx = 0

for i, atlas_pos in enumerate(atlas_positions):
    for j, nh_pos in enumerate(nh_positions):
        distance = np.linalg.norm(atlas_pos - nh_pos)
        if distance < min_distance:
            min_distance = distance
            closest_atlas_idx = i
            closest_nh_idx = j

# Results
closest_atlas_pos = atlas_positions[closest_atlas_idx]
closest_nh_pos = nh_positions[closest_nh_idx]

print(f"\nRESULTS:")
print(f"Minimum distance: {min_distance:.2f} AU")
print(f"Minimum distance: {min_distance * 149.6:.0f} million km")

print(f"\n3I/ATLAS position at closest approach:")
print(f"  X: {closest_atlas_pos[0]:.2f} AU")
print(f"  Y: {closest_atlas_pos[1]:.2f} AU")
print(f"  Z: {closest_atlas_pos[2]:.2f} AU")
print(f"  Distance from Sun: {np.linalg.norm(closest_atlas_pos):.2f} AU")

print(f"\nNew Horizons position at closest approach:")
print(f"  X: {closest_nh_pos[0]:.2f} AU")
print(f"  Y: {closest_nh_pos[1]:.2f} AU")
print(f"  Z: {closest_nh_pos[2]:.2f} AU")
print(f"  Distance from Sun: {np.linalg.norm(closest_nh_pos):.2f} AU")

# Create 3D visualization
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection="3d")

# Plot Sun at origin
ax.scatter([0], [0], [0], color="yellow", s=100, label="Sun")

# Plot trajectories
ax.plot(
    atlas_positions[:, 0],
    atlas_positions[:, 1],
    atlas_positions[:, 2],
    "r-",
    label="3I/ATLAS trajectory",
    alpha=0.7,
)
ax.plot(
    nh_positions[:, 0],
    nh_positions[:, 1],
    nh_positions[:, 2],
    "b-",
    label="New Horizons trajectory",
    alpha=0.7,
)

# Mark closest approach points
ax.scatter(
    closest_atlas_pos[0],
    closest_atlas_pos[1],
    closest_atlas_pos[2],
    color="red",
    s=100,
    label=f"3I/ATLAS closest",
)
ax.scatter(
    closest_nh_pos[0],
    closest_nh_pos[1],
    closest_nh_pos[2],
    color="blue",
    s=100,
    label=f"New Horizons closest",
)

# Draw line between closest points
ax.plot(
    [closest_atlas_pos[0], closest_nh_pos[0]],
    [closest_atlas_pos[1], closest_nh_pos[1]],
    [closest_atlas_pos[2], closest_nh_pos[2]],
    "k--",
    alpha=0.5,
    linewidth=2,
    label=f"Min distance: {min_distance:.1f} AU",
)

ax.set_xlabel("X (AU)")
ax.set_ylabel("Y (AU)")
ax.set_zlabel("Z (AU)")
ax.set_title("3I/ATLAS vs New Horizons Trajectories\nClosest Approach Calculation")
ax.legend()

# Set equal aspect ratio
max_range = 70
ax.set_xlim([-max_range, max_range])
ax.set_ylim([-max_range, max_range])
ax.set_zlim([-10, 10])

plt.tight_layout()
plt.show()

print(f"\nNOTE: This is an approximation based on:")
print(f"- Limited orbital element data for both objects")
print(f"- Simplified New Horizons trajectory model")
print(f"- Two-body orbital mechanics (ignoring planetary perturbations)")
print(f"- The actual minimum distance could vary significantly from this estimate")