← Back
Personal / Independent Project (Computational Astrophysics)

Solar System & Halley's Comet 3D Orbit Simulation

AdvisorSelf-directed (Computational Physics)
PeriodFeb 2024

This project implements a high-fidelity, Python-based 3D dynamic orbit simulation of the Solar System. The simulator models the orbits of the eight major planets (Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune) along with the highly eccentric, highly inclined retrograde trajectory of Halley's Comet (1P/Halley). By inputting a specific astronomical epoch (date and hour in UTC), the tool dynamically calculates the exact 3D heliocentric coordinates and outputs relative Euclidean distances between any two celestial bodies, recreating real-time planetary alignments with high physical accuracy.

  • ▪Keplerian Orbital Parameter Integration: Includes semi-major axis (a), eccentricity (e), orbital inclination (i), longitude of ascending node (Ω), and perihelion longitude (ϖ) derived from standard astrodynamic datasets.
  • ▪3D Coordinate Euler Transformations: Implemented coordinate conversion from 2D orbital planes to 3D ecliptic/heliocentric coordinate systems using rigorous spatial rotation matrices.
  • ▪Ephemeris Distance Calculation: Computes instant relative Euclidean distances between any two planets or celestial bodies for any given date and hour (UTC).
  • ▪Astrophysical Data Verification: Pre-loaded boundary and initial values aligned precisely with actual orbital states starting from February 22, 2024.

Methodology

1. 3D Spatial Rotation Transformations (Euler Rotations)

Celestial coordinates on their native 2D orbital plane are represented by standard parametric ellipse equations:

xorb=acos⁡(θ)+cx_{\text{orb}} = a \cos(\theta) + c
yorb=bsin⁡(θ)y_{\text{orb}} = b \sin(\theta)

where c=a⋅ec = a \cdot e is the focal offset distance.

To map these coordinates onto the actual 3D heliocentric ecliptic coordinate system, the project executes a sequential 3D spatial Euler transformation:

1. First, the coordinates are tilted using the orbital inclination (ii) relative to the ecliptic plane:

z′=xorbsin⁡(i)z' = x_{\text{orb}} \sin(i)
x′=xorbcos⁡(i)x' = x_{\text{orb}} \cos(i)
y′=yorby' = y_{\text{orb}}

2. Second, a custom rotation using the spatial rotation matrix is applied to perform 2D plane rotations around the ascending node (Ω\Omega) and perihelion longitude (ϖ\varpi):

(xfinalyfinal)=R(ϕ)(x′y′)\begin{pmatrix} x_{\text{final}} \\ y_{\text{final}} \end{pmatrix} = R(\phi) \begin{pmatrix} x' \\ y' \end{pmatrix}

This maps the tilted trajectory from its isolated orbital plane into the unified heliocentric ecliptic coordinates.

2. Adaptive Time-Stepping & Keplerian Integration

To maintain high numerical accuracy, particularly for cometary bodies with high eccentricities, two adaptive numerical methods are implemented:

A. Keplerian Midpoint Velocity Approximation:

For major planets with low eccentricities, orbital velocity scales inversely with distance (rr) to conserve angular momentum. A midpoint speed correction is implemented during integration:

vn=vavg⋅0.5⋅(1rn+1rn+1)v_n = v_{\text{avg}} \cdot 0.5 \cdot \left(\frac{1}{r_n} + \frac{1}{r_{n+1}}\right)

This ensures planets realistically accelerate near perihelion and decelerate near aphelion, verifying Kepler's Second Law.

B. Vis-Viva Equation Integration for Halley's Comet:

Halley's Comet has an extreme orbital eccentricity (e≈0.9679e \approx 0.9679) and an inclination of 162.18∘162.18^\circ (retrograde). Gravity changes dramatically near perihelion. The simulation utilizes the energy-conserving Vis-Viva Equation to model this motion:

v=μ(2r−1a)v = \sqrt{\mu \left( \frac{2}{r} - \frac{1}{a} \right)}

where μ=1.327×1011 km3/s2\mu = 1.327 \times 10^{11} \text{ km}^3/\text{s}^2 is the Heliocentric Gravitational Constant. This ensures exceptional physical validity and numerical stability even at perihelion passage.

Results & Analysis

The simulation successfully produces highly accurate 3D trajectory datasets and renders dynamic Matplotlib animations showcasing the inner solar system and cometary crossings.

Ephemeris Relative Distance Calculations

Using the built-in return_distance2(year, month, day, hour, planet1, planet2) utility, the system calculates the real-time Euclidean distance between any two celestial bodies. This serves as an intuitive tool for studying synodic periods, planetary alignments, and cometary approaches.

All simulated trajectory data has been exported to structured text files (Earth.txt, Mars.txt, Jupiter.txt, Halley.txt) and validated against NASA JPL Horizons ephemeris data, exhibiting close alignment in short-to-medium-term orbital predictions.

Figures

Figure 1: 3D Dynamic Orbit Simulation of inner planets and Halley's Comet (Matplotlib animation output)

Figure 1: 3D Dynamic Orbit Simulation of inner planets and Halley's Comet (Matplotlib animation output)

Figure 2: Earth-centred view showing the angular positions of Venus, Earth, and Mars — illustrating synodic angle variation over time

Figure 2: Earth-centred view showing the angular positions of Venus, Earth, and Mars — illustrating synodic angle variation over time