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. 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?
 A: I concur with the qualification that you have recreated in software a tellurion
To put your creation in perspective: the Antikythera device offered a representation of the motion of the Moon that is superior to what you created in software.

With a mechanical model the planets move like the hands of a clock; they move mechanically. What that means is that in effect you are not even modeling gravity.
(It is inevitable of course that in your model the Moon-Earth distance is not to scale with the Sun-Earth distance; every mechanical model has to deal with the distances of the solar system in one way or another.)
The animated GIF is about 1.5 MB, and you have bought very little for that.
You did take care to make it a looping animated GIF: when the 240 frames are completed the loop closes seamlessly.
Of course, in order to make the animation loop seamlessly you had to take liberty with the duration of the Moon orbit.
My recommendation:
Look up how the Antikythera mechanism gave the Moon a non-uniform angular velocity and incorporate that in the motion that goes into creating the frames of the animated GIF.

You mention the possibility of including relativistic effects at some future point. But the liberties you are already taking are far, far, far greater than any relativistic effect.

I'm assuming the intented usage is publication in the form of an animated GIF.
I'm assuming that is your constraint.
(Without that constraint you would not be making an animation, but a simulation, with the motions of the celestial bodies being calculated in real time.)
Within the constraint of creating an animated GIF the best you can do is to include the non-uniform angular velocity of the Moon. Anything more is beyond an animated GIF anyway.
A: There are aspects of the solar system that are not well approximated by this model.  For example in the actual sole system the moon never goes backwards in its orbit around the sun.
A: You have created a very simple digital orrery - in fact, since it only includes the Sun, the Moon and the Earth it is technically a tellurion. Since you have intentionally omitted or approximated so many aspects of the actual motion, including scales and time periods, it is impossible to tell if you have accidentally omitted anything.
A: The orbit of the Moon could not be circular. As the Moon moved toward the Sun, in its orbit around Earth, it would accelerate towards the Sun. Then as it moved around Earth and started moving away from the Sun it would be slowed, this would result in an elliptical lunar orbit giving Earth a non-circular wobble around around the Earth - Moon barycenter.
