Have I accounted for all relevant physics in my simplified three-body solar system model?
EDIT: Re-framed here to be more on-topic.
I am building a model of a simplified, fictional sun-planet-moon solar system (not the real Sun-Earth-Luna system) to calculate astronomical events like eclipses and lunar phases and the positions of the sun and moon in the sky of the planet at a given time, and want to know if there's anything I'm missing that would affect the accuracy of the model over the scale of several hundred years (that isn't being purposefully left out; see below). I want to be physically accurate within the restricted scope of the model.
In the model, all orbits are circular and in the same plane. At $t = 0$ the three bodies are colinear, with the moon located between the sun and the planet. 319
The sun-planet-moon barycenter is located at the origin. The sun and the planet-moon barycenter both orbit the sun-planet-moon barycenter in 360 days; the planet and the moon both orbit the planet-moon barycenter in $\frac{360}{13}$ (≈27.69) days, with one synodic month (time between full/new moons) being exactly 30 days. The following diagram (very much NOT to scale) illustrates these motions:
* 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-plus-moon subsystem: (PlanetBarycenterDistance * (MassPlanet + MassMoon)) / (MassSun + MassPlanet + MassMoon) - PlanetOrbitRadius is the distance between the planet-moon barycenter and the sun-planet-moon barycenter, calculated based on the mass of the sun and the orbital period of the planet: cbrt((TimePlanetOrbit² * G * MassSun) / (4π²)) - PlanetBarycenterDistance is the distance between the planet-moon barycenter and the center of the planet: (MoonOrbitRadius * MassMoon) / (MassPlanetPlusMoon) - MoonOrbitRadius is the distance between the planet-moon barycenter and the center of the moon: cbrt((TimeMoonOrbit² * G * MassPlanet) / (4π²)) - G is the gravitational constant
Things I am purposefully leaving out: - non-circular orbits causing changes in orbital velocity at different points in the orbit; as mentioned above, the model's orbits are circular - 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 - precession of the equinoxes and the difference between sidereal and solar years caused by general relativity and the oblateness of the sun and 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 - 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
Aside from these effects that I am purposefully leaving out of the model, are there any effects or motions that should be considered? 380