Mathematics & Computation

This page explains (in student-friendly terms) how the simulation computes planet positions using Orbital Elements and Kepler's laws.

1. Orbital Elements

Most planets move around the Sun on an ellipse (an oval shape). Keplerian orbital elements (Orbital Elements) are 6 numbers that describe the size, shape, tilt, and direction of that ellipse:

• a — Semi-major axis (AU): half the longest diameter of the ellipse.
• e — Eccentricity (0 = circle, <1 = ellipse): how elongated the orbit is.
• I — Inclination (°): tilt of the orbital plane relative to the ecliptic.
• Omega — Longitude of the ascending node (°): where the orbit crosses the ecliptic going north.
• wBar — Longitude of perihelion (°): direction of the closest approach to the Sun.
• L0 — Mean longitude at epoch (°): where the planet was at the reference time J2000.0.

These values change slowly over time because planets pull on each other. This simulation uses a simple linear approximation for that slow change.

2. Julian Date

In astronomy, time is often written as a single day counter called Julian Date (JD). A common reference time is J2000.0 (noon on 2000-01-01), which equals JD 2451545.0.

A JavaScript Date object stores milliseconds from the Unix epoch (1 Jan 1970 00:00 UTC), so the conversion to JD is:

3. Mean Anomaly

If a planet moved around its orbit at a constant speed, its angle would be the mean longitude L:

If we subtract the direction of perihelion (wBar), we get the mean anomaly M:

Real planets move faster when they are closer to the Sun (Kepler's second law). So the real angle is different from M. That real angle is called the true anomaly (here we write it as nu).

4. Kepler's Equation

To go from M to nu, we first solve for the eccentric anomaly E. (It is another angle that is easier to compute.) Kepler's equation is:

We cannot solve this equation for E in one step, so we use repeated steps (numerical iteration). This simulation uses the Newton-Raphson method:

For solar-system planet eccentricities (e < 0.21 for Mercury, ≤ 0.09 for all others) this converges to machine precision in fewer than 5 iterations.

5. True Anomaly & Heliocentric Distance

Once we know E, we can compute the true anomaly (nu) and the heliocentric distance r (AU):

atan2 helps keep nu in the correct direction (quadrant).

6. Heliocentric Ecliptic Coordinates

The planet's position in its own orbital plane is:

Three successive rotations transform these to the heliocentric ecliptic frame (J2000.0):

• Rotate by omega (argument of perihelion = wBar - Omega) about the z-axis.
• Rotate by I (inclination) about the x-axis.
• Rotate by Omega (longitude of ascending node) about the z-axis.

The resulting (x, y, z) are in AU with the Sun at the origin, x pointing toward the vernal equinox, and z pointing north of the ecliptic.

7. Drawing the Orbit Path

To draw the orbit line, we sample nu from 0 to 360 degrees and compute a point each time. For each nu, the distance from the Sun is:

Applying the same rotation matrices as above gives 360 points that are fed directly into a Three.js LineLoop geometry.

8. Accuracy & Limitations

• VSOP87 (Bretagnon & Francou, 1987): analytical series with thousands of terms giving ~1-arcsecond accuracy over several millennia.
• JPL DE440: numerical integration of full equations of motion, valid from 1550–2650, sub-kilometre accuracy.
• Neither model includes asteroid perturbations on inner planets, which accumulate over very long timescales.

9. Step-by-Step Algorithm Summary

1. Convert the chosen date to Julian Date (JD).
2. Compute T = (JD - 2451545) / 36525.
3. Evaluate all six orbital elements at time T.
4. Compute mean anomaly M = L - wBar.
5. Solve Kepler's equation for eccentric anomaly E (Newton-Raphson).
6. Compute true anomaly nu and distance r.
7. Rotate to heliocentric ecliptic (x, y, z) in AU.
8. Scale to screen coordinates and render in Three.js.