I am building a model of a simple sun-planet-moon solar system (not the real Sol-Terra-Luna system).
Purpose: to calculate the positions of the sun and moon in the sky of the planet at a given time, and the timing of astronomical events like eclipses and moon phases.
Question: how do I account for the gravity of the sun acting on the planet-moon sub-system?
The starting point for my model is a two-body Keplerian orbital system, the planet-moon system, with the planet and moon in circular orbits about their shared barycenter. This system is then treated as a single object in a larger two-body Keplerian orbital system, the sun-[planet-moon] system, with the sun and the planet-moon barycenter in circular orbits about their shared barycenter. As it currently stands, the sun does not affect the orbit of the moon about the planet. This is what I want to fix.
In reality, the orbits of the planet and moon about their shared barycenter would not be perfectly circular, and the speed of their orbits about their shared barycenter would not be constant, due to the gravitational influence of the sun. How do I calculate the physically-accurate motions?
In the real Sun-Earth-Moon system, the deviation of the moon's orbit from uniform circular motion has several named components:
- The largest is known as the equation of the center (EOC), and accounts for the fact that the moon's orbit is slightly eccentric (shaped like an ellipse, not a perfect circle, with the Earth[-moon barycenter] located at one focus of the ellipse). In my system the orbital eccentricity is zero, so this component is not present.
- The next largest component is evection, caused by the sun, where the moon's longitude varies from its position as predicted by the simple circular orbit and EOC according to the term $+1.27401° \sin(2D - ℓ)$, where $D$ is the mean angular distance of the Moon from the Sun (elongation) and $ℓ$ is the mean angular distance of the Moon from its perigee (mean anomaly). I want to calculate what this term would be for my custom solar system.
- The third largest component in the real moon's perturbation is variation, also caused by the sun, where the moon is closer to Earth and moving faster at full-moon and new-moon phases, and farther from Earth and moving slower at quarter-moon phases, with longitude varying according to the term $+0.6583° \sin(2D)$. This makes the orbit an ellipse, but with the Earth at the center of the ellipse rather than at a focus. I want to calculate what this term would be for my custom solar system, and also how the distance varies from this effect.
The total deviation of the moon's position in the sky caused by evection and variation only gets up to about two degrees, but since the moon is only about half a degree wide in the sky, that amounts to being off by up to four times the width of the moon - a pretty hefty distance. I'd like to get the maximum error in my model down to significantly less than that.
In the model, the sun-planet-moon barycenter is located at the origin of the reference frame, with the plane of the ecliptic as the $xy$-plane and the orbital axes parallel to the $z$-axis. All orbits are circular (or at least have zero eccentricity) and are in the plane of the ecliptic. At $t = 0$ the three bodies are colinear along the $x$-axis, with the moon located between the sun and the planet. The sun and the planet-moon barycenter both orbit the sun-planet-moon barycenter in exactly 360 local solar days; the planet and the moon both orbit the planet-moon barycenter in approximately $\frac{360}{13}$ (≈27.7) local solar days, with one synodic month (time between full/new moons) being exactly 30 local solar days. One local solar day is 75,325.804 Earth seconds (approximately 21 Earth hours), making one local sidereal day 75,117.145 Earth seconds.
In the model as it currently stands:
- The center of the sun has coordinates (−SunBarycenterDistance * cos($t$), −SunBarycenterDistance * sin($t$))
- The center of the planet has coordinates (PlanetOrbitRadius * cos($t$) + PlanetBarycenterDistance * cos(13$t$), PlanetOrbitRadius * sin($t$) + PlanetBarycenterDistance * sin(13$t$))
- The center of the moon has coordinates (PlanetOrbitRadius * cos($t$) − MoonOrbitRadius * cos(13$t$), PlanetOrbitRadius * sin($t$) − MoonOrbitRadius * sin(13$t$))
Where
- SunBarycenterDistance is the distance between the sun-planet-moon barycenter and the center of the sun, calculated based on the masses of the sun and the planet-moon subsystem and the total distance between them: (DistanceSunPlanetMoon (MassPlanet + MassMoon)) / (MassSun + MassPlanet + MassMoon)
- PlanetOrbitRadius is the distance between the sun-planet-moon barycenter and the planet-moon barycenter, calculated based on the masses of the sun and the planet-moon subsystem and the total distance between them: (DistanceSunPlanetMoon * MassSun) / (MassSun + MassPlanet + MassMoon)
- PlanetBarycenterDistance is the distance between the planet-moon barycenter and the center of the planet, calculated based on the masses of the planet and the moon and the total distance between them: (DistancePlanetMoon * MassMoon) / (MassPlanet + MassMoon)
- MoonOrbitRadius is the distance between the planet-moon barycenter and the center of the moon, calculated based on the masses of the planet and the moon and the total distance between them: (DistancePlanetMoon * MassPlanet) / (MassPlanet + MassMoon)
- DistanceSunPlanetMoon is the distance between the center of the sun and the planet-moon barycenter, calculated based on the masses of the sun, planet, and moon and on the orbital period of the planet: cbrt((TimePlanetOrbit² * G * (MassSun + MassPlanet + MassMoon)) / (4π²))
- DistancePlanetMoon is the distance between the center of the planet and the center of the moon, calculated based on the masses of the planet and moon and on the orbital period of the moon: cbrt((TimeMoonOrbit² * G * (MassPlanet + MassMoon)) / (4π²))
- G is the gravitational constant
One thing that was suggested to deal with lunar variation is putting the moon on an additional small epicycle, making the coordinates for the center of the moon something like (PlanetOrbitRadius * cos($t$) − MoonOrbitRadius * cos(13$t$) + EpicycleRadius cos(−11$t$), PlanetOrbitRadius * sin($t$) − MoonOrbitRadius * sin(13$t$) + EpicycleRadius sin(−11$t$)), but to do that I need to figure out the radius of the epicycle and whether or not the speeding-up and slowing-down effect of this epicycle align with the actual speeding-up and slowing-down effect of the sun's gravity. I assume I would also need to put the planet on a similar epicycle, sized and timed to make sure the center of mass of the planet-moon system always tracks with the modeled barycenter point.
Things from the real Sun-Earth-Moon system that I am purposefully leaving out:
- eccentric orbits causing changes in orbital velocity at different points in the orbit; as mentioned above, the model's orbits have zero eccentricity
- apsidal precession; because the orbits are circular, there are no apsides
- nodal precession; because the orbits are in the same plane, there are no orbital nodes
- perturbations from other bodies; the model only considers the one sun, one planet, and one moon, with no other masses acting on these three bodies
- precession of the equinoxes caused by general relativity and the oblateness of the sun and planet and the axial tilt of the planet; I'm still working out how to calculate this and want to include it in the future, but I'm leaving it out for now
ADDENDUM: I received a question as to the validity of having the moon on a circular orbit around the planet-moon barycenter even as a first approximation, from the thought that this would make the moon's path sometimes curve away from the sun. Here is an image showing the moon's path (black) superimposed over the perfect circle of the barycenter path (red). Note that the black path is always curving towards the sun, though the curvature and speed relative to the sun increase and decrease as the moon circles inside and outside the planet:
And here is a GIF of the planet and moon orbiting the sun, centered on the planet-moon barycenter, with the moon's trail shown as a black line:
