Introduction
The initial intention was to use Newton's laws as the sole starting point to simulate planetary motion in our solar system, using stepwise simulation of orbits:
However, the method turned out to be useful to try a whole bunch of other stuff. The possibilities of playing and experimenting with this simulation software are immense, and the number of remarkable phenomena in our solar system is impressive.
The journey took us from the orbital resonance of the Galilean moons, over the Sun's orbit revealing the presence of its planets, up to the Lissajous-like orbits of Trojan asteroids around Jupiter's Lagrange points.
In what follows we will start by introducing the (simple) mathematics required to implement the simulation. If you are interested in the results, you may skip to the section on Solar System Simulation.
Newton's Law of Universal Gravitation
To develop the simulation, we start off with Newton's law of universal gravitation, defining the strength of the gravitational force between two bodies as proportional to the mass of both bodies and inversely proportional to the square of the distance between the bodies, with a factor G:
$$F=G{\frac {m_{1}m_{2}}{r^{2}}}.$$ Here, G is the gravitational constant, whose value has been measured to be:
$$G = 6.6740831\,10^{11} \frac{m^3}{kg\,s^2}.$$ By the way, thanks to MathJax for their TeX support featuring beautiful math in all browsers.
The law also states that the force acts along the line intersecting the two bodies, in the direction of the other body. Therefore, to obtain the force vector, we multiply the unit vector in this direction with the strength of the force:
$${\vec F}=G{\frac {m_{1}m_{2}}{r^{2}}}{\frac {{\vec p_2}-{\vec p_1}}{r}},$$ where p represents the position of a body, and knowing that
$$r = \|{\vec p_2}-{\vec p_1}\|.$$ This allows us to calculate, for a specific body, as a vector, the net force exerted on that body by all other bodies, as the vector sum of the force vectors associated with all other bodies: $${\vec F}_T = {\sum_{i=1}^{N} {\vec F}_i},$$ where N is the number of bodies involved in the simulation.
The calculation of the net force vector is illustrated with Java code at the end of the post, in the section on Celestial Body Representation. The implementation of the vector calculations involved is presented in the section on Vector Representation and Calculations.
Newton's Laws of Motion
We will then be using Newton's laws of motion:
- I - A body either remains at rest or continues to move at a constant velocity, unless acted upon by a force.
- II - The vector sum of the forces on a body is equal to the mass of that body multiplied by the acceleration.
- III - When one body exerts a force on a second body, the second body simultaneously exerts a force equal in magnitude and opposite in direction on the first body.
$${\vec F}=m{\vec a}\;\Rightarrow\;{\vec a} = \frac{\vec F}{m}.$$ This allows us to calculate the acceleration of each body out of the net force vector and its mass.
We can then revert to the definitions of acceleration and velocity, again rewritten to our needs:
$${\vec a} = \frac{d{\vec v}}{dt}\;\Rightarrow\;d{\vec v} = {\vec a}\,dt,\;{\vec v} = \frac{d{\vec p}}{dt}\;\Rightarrow\;d{\vec p} = {\vec v}\,dt,$$ where v represents the velocity and p represents the position of a body.
For a single simulation step, we will be using the discrete equivalent of these equations:
$$\Delta {\vec v} = {\vec a} \Delta t,\;\Delta {\vec p} = {\vec v} \Delta t.$$ This allows us to calculate the change in velocity out of the acceleration and the time interval of the simulation step, and to calculate the change in position out of the velocity and the time interval of the simulation step.
The implementation of a single simulation step is presented at the end of this post, in the section on The Simulation Step. We explicitly detailed the use of the second law of motion. The first and third laws of motion are equally used, implicitly wired in the control structure of the Java code snippet.
Celestial Body Characteristics
The next step, before we can simulate anything, is to obtain the characteristics of the celestial bodies we wish to include in our simulation. Although not strictly necessary for the simulation, we prefer working with the real positions and velocities of the celestial bodies.
The mass and average distance to the sun of all planets can easily be found on the Internet. However, synchronous information on the position and velocity of the planets is a lot harder to find.
Looking for this information, we stumbled upon the Jet Propulsion Laboratory Horizons System. It offers detailed information for a vast number of celestial bodies, accessible using different interfaces. Its use is documented in the section Jet Propulsion Laboratory to the Rescue.
For this exercise, we requested positions and velocities at midnight on 19-07-2019, expressed as vector coordinates against the ecliptic with the sun as the origin.
As an example, the position of the earth on the aforementioned date and time is obtained as approximately: $$(6.61850\,10^{10} m, -1.36870\,10^{11} m, 6.33303\,10^6 m),$$ and the velocity is obtained as approximately: $$(2.63245\,10^4 \frac {m}{s}, 1.28468\,10^4 \frac {m}{s}, 2.89760\,10^{-1} \frac {m}{s}).$$
Elimination of Net Impulse
With all the collections of celestial bodies we used in a simulation, we systematically observed a net impulse.
We suspect this is caused by the lack in velocity for the sun in our data, as the sun is used as a reference point in the coordinate system. However, in our simulations, we will let the sun move around as any other celestial body, subject to the gravitational forces of all other bodies.
Impulse is defined as the product of mass and velocity: $${\vec I} = m{\vec v}.$$ The net impulse for all bodies involved in a simulation is given by $${\vec I}_T = {\sum_{i=1}^{N} m_i {\vec v}_i},$$ where N is the number of bodies involved in the simulation.
This net impulse causes a common drift of the bodies involved in the simulation, with constant velocity. This velocity can be calculated using the net impulse and the total mass of all bodies: $${\vec v}_T = \frac {{\vec I}_T}{m_T}, m_T = \sum_{i=1}^{N} m_i.$$ For example, when we calculate the net impulse for a simulation involving the Sun and eight planets, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus and Venus, it is approximately: $$(1.45843\,10^1 {\frac {m}{s}}, 2.53651\,10^{-1} {\frac {m}{s}}, -3.81053\,10^{-1} {\frac {m}{s}}).$$ To avoid drift, for each simulation, and consequently for each collection of celestial bodies, we systematically correct all initial velocities to compensate for this drift, by subtracting the velocity associated with the net impulse and the total mass.
Projection on the Ecliptic
Still we are not ready to start simulating. The simulations we will be running have to be visualised.
We will take the path of least resistance for this. As we requested information in reference of the ecliptic, and planets are orbiting the sun largely parallel with that plane, we can merely set the Z-coordinate to zero, implementing a parallel projection onto the ecliptic.
Scaling the projection onto a drawing canvas typically involves flipping the Y-coordinate and scaling both the X-axis and the Y-axis with the same scaling factor. We will leave this transformation as an exercise.
Solar System Simulation
Now we are ready to run our first simulation, finally. We will simulate the Sun, with the eight planets mentioned before. For nostalgic reasons, we will add Pluto to the mix. And Halley's Comet. The result can be observed as an animated gif:
In the simulation, we observe the orbits of the eight planets. The time is indicated at the bottom, measured in earth years. Pluto completes its orbit in approximately 248 years. Its orbital eccentricity is clearly visible. During part of the orbit, Pluto even moves closer to the Sun than Neptune, luckily without ever coming too close to Neptune itself. If you look carefully, you can observe the eccentricity of Mars' orbit as well.
We can also observe Halley's comet has a retrograde orbit, moving clockwise rather than counterclockwise, from our perspective, and that it completes somewhat over 3 orbits during the 248 year period.
All these observations are largely consistent with reality.
A Note on Precision
The simulation technique we are using is simplistic. It is expected to generate significant errors in simulations. In addition, errors will accumulate over time. More precise orbit models have been designed and implemented. The Jet Propulsion Laboratory Development Ephemeris can give you a good impression, with precision of 25 meters for the position of any planet in the solar system.
The step size we have been using for the previous simulation is 12 hours. This means that we implicitly assume all movements to be linear for a period of 12 hours, after which new velocities, with a different direction, are calculated. Clearly, this will have a negative impact on the precision of the simulation.
All other simulations have been done with a step size of 1 hour. This step size has been chosen pragmatically, finding the balance between simulation time and the results we are looking for. Comparing simulated positions with data provided by the Jet Propulsion Laboratory Horizons System over a period of 10 earth years gives accumulated errors below 0.1% for the Jupiter orbit and below 0.35% for the orbit of comet Halley.
The Venus Pentagram
As a frequent reader of Science Alert, we recently stumbled upon the article Invisible Dance of Earth And Venus Forms a Stunning Pentagrammic Pattern in Space.
The periods of the orbits of Earth and Venus are very close to having a relatively small common multiple. As the original post by Matthew Henderson on Tumblr states, eight Earth years are roughly equal to thirteen Venus years, meaning the two planets approximately trace out this pattern with 5-fold symmetry as they orbit the Sun.
Being low-hanging fruit, we could not resist replicating this idea with our simulation software:
Indeed, drawing a line connecting the centers of mass of the earth and Venus every five days, a similar pentagrammic animation appears. The resulting spirographic image is a little less regular than the original, but that is because we did not simplify the scenario to have the orbits perfectly in sync.
Galilean Moons
We will now investigate whether we can add some satellite bodies, by taking the Galilean moons into the simulation. The Galilean moons are the four largest moons of Jupiter: Io, Europa, Callisto and Ganymede.
Simulating one orbit of Jupiter accompanied by those four moons, and zooming in on a fragment of that orbit, we get an interesting result:
One can clearly see the moons engaged in a kind of leapfrog game. The attentive reader will also notice the orbital resonance of the first three moons, Io, Europa and Ganymede, with ratio's 4-to-2-to-1 respectively.
The Sun's Orbit
We are commonly talking about the orbit of a planet around the sun. However, in a two-body system, both bodies are orbiting around the center of mass (barycenter) of the system. This obviously extrapolates to N-body systems.
Therefore, as mentioned before, during our simulation, we will calculate and update the position and velocity of the Sun, just as for any other celestial body. This naturally raises our curiosity on the orbit of the Sun.
Obviously, the major impact on the Sun's orbit will come from Jupiter and Saturn, being by far the heaviest of the planets. Simulating a system featuring only the Sun and these two gas giants gives a nice picture:
The impact of two different frequencies in the perturbation of the Sun's position is clearly visible. Given the ratio of the periods of the individual oscillations that make up the composite oscillation, the orbit of the Sun kind of resembles a spirograph creation.
The Sun's volume is drawn in dark grey at the start of the simulation, to get an idea of the relative size of the Sun's orbit. The animation however involves only the center of mass of the Sun.
Repeating this experiment with all eight planets and Pluto on board, gives a different picture:
Here, the composition of many oscillations with different frequencies and amplitudes creates a picture of the Sun's orbit that seems rather chaotic and unpredictable. But notice the similarity of the orbit with the image on the wikipedia page on astrometry showing the motion of the barycenter of the solar system relative to the Sun, both in terms of the relative sizes of the Sun and the orbit and in terms of the orbit's shape.
An Excursion into Astrometry
Although somewhat off-topic to the story being told, we could not resist an excursion into the realm of exoplanet discovery. Skip to the section on Lagrange Points for more simulations.
There are many different methods of detecting exoplanets. Technology related to this theme is evolving rapidly. One of the older methods is based on Astrometry, where changes in position of a star may reveal a companion in orbit.
Although this method of detection has been unsuccessful in the past, the space-based observatory Gaia, launched in 2013, is expected to find thousands of planets via astrometry. An interesting follow-up of the previous simulation is to check whether we would be able to detect the contributions of individual planets in perturbations of the Sun's orbit.
A Fourier transform can be used to decompose a function of time into its constituent frequencies, and it is just the technique we would use for this exercise. We recommend the Visual Introduction to Fourier Transforms by 3Blue1Brown as an excellent read.
Rerunning the last simulation, with eight planets and Pluto in orbit around the Sun, and registering the position of the Sun on the X-axis every 8 days over a period of approximately 718 years, we get our input data for the Fourier transform: a time series with 32,678 (=2^15) data points.
The resulting frequency spectrum coming out of the Fourier transformation confirms our hypothesis: there are 8 planets orbiting this star.
Up to grid point 200, representing orbits with a period exceeding 3.59 year, we observe Neptune, Uranus, Saturn and Jupiter. On grid point 121, representing and orbit period of 5.94 year, we observe a small peak that does not correspond to any of the known planets. Might there be something out there ?
Up to grid point 2000, representing orbits with a period exceeding 0.36 year, we observe Mars, Earth and Venus.
Up to grid point 3000, representing orbits with a period exceeding 0.24 year, we observe Mercury.
The only planet we cannot not observe in the frequency diagram is Pluto. Its influence is expected to be very weak, and our data set covers less that three Pluto orbits. Anyway, that is in line with the 2006 regulations of the International Astronomical Union, that excluded Pluto as a planet and reclassified it as a dwarf planet.
We hereby provide an overview of the peaks observed in the frequency diagram and the matching planets in our simulation:
| # | Planet | Estimated Period | Official Period | Error |
|---|---|---|---|---|
| 4 | Neptune | 179.55 | 163.84 | +9.59% |
| 9 | Uranus | 79.80 | 83.81 | -4.78% |
| 24 | Saturn | 29.92 | 29.44 | +1.63% |
| 61 | Jupiter | 11.77 | 11.87 | -0.84% |
| 381 | Mars | 1.89 | 1.88 | +0.53% |
| 719 | Earth | 1.00 | 1.00 | 0.00% |
| 1,166 | Venus | 0.62 | 0.61 | +1.64% |
| 2,979 | Mercury | 0.24 | 0.24 | 0.00% |
Note that the rather large errors for more remote planets here are mainly introduced by the fixed frequency bucket size imposed by the Fourier transform, giving a variable period bucket size, with large buckets, and consequently large errors on the side with long periods.
The implementation of the Fourier transform is detailed in the section Winding Time Series.
Lagrange Points
Lagrange points are the points near two large bodies in orbit where a smaller object will maintain its position relative to the large orbiting bodies.
In a two-body system, there are five Lagrange points, labeled L1 to L5. Lagrange points L1, L2 and L3 are all located somewhere on the axis connecting the two bodies. They are unstable, in the sense that eventually, given enough time, a small body will drift away from the Lagrange point.
Lagrange points L4 and L5 each form an equilateral triangle with the large bodies, more or less, respectively ahead and behind in orbit. They are stable under certain conditions, in the sense that a small body will maintain a stable orbit around the Lagrange point.
Observational spacecraft are often put in orbit around a Lagrange point, with the main advantage that they have to spend less resources to maintain a stable orbit. Even L1 and L2, although not stable, offer this advantage, as their instability is only in the radial direction, while they offer stability in directions orthogonal to the direction of the radius. The effective potential in such a Lagrange point is a saddle point.
See the Document Created by Neil Cornish for WMAP for an analytical approach on Lagrange points. For now, we will continue our experiment trying to simulate a body orbiting a Lagrange point.
In order to experiment with Lagrange points and their effects, we will have to extend the visualisation logic, requiring some additional (simple) mathematics. If you are interested in the results, you may skip to the section on Pluto's Dance.
A Rotating Frame
For the visualisation of an orbit around a Lagrange point, we need to use a rotating frame, where the perspective we take follows the rotation of the large bodies in orbit at the same angular velocity. If not, we would just observe a wobbly orbit around the sun, obfuscating the fine details of all motion relative to the large bodies and the Lagrange points.
With respect to our current visualisation logic, we basically need to add a rotation that compensates for the rotation of the two large bodies. We will insert this after the projection of the ecliptic, such that it becomes a 2D problem, and before the scaling to the drawing canvas.
In order to do that, we first need to know the angle of the line connecting the two large bodies. We would use the inverse sine, cosine or tangent function and the values of the X- and Y-coordinates for this. However, we need to be careful to avoid discontinuities in our function. In fact, the angle should span the range from 0 to 2π radians when rotating the unit vector on the X-axis 360 degrees around the Z-axis.
We will be using a formula proposed by Theodore Panagos on StackExchange for this purpose, choosing the tangent variant, as this is the shortest: $$\alpha = \pi-{\frac {\pi}{2}}(1+sgn(x))(1-sgn(y^2))-{\frac {\pi}{4}}(2+sgn(x))sgn(y)\\-sgn(xy)atan({\frac {\vert x\vert-\vert y\vert}{\vert x\vert+\vert y\vert}}).$$ Do not let the sheer size of this monster frighten you. It is merely an arithmetical trick to create a kind of if-then-else construct. The sign functions will activate the right part of the formula, depending on the quadrant we are in with particular values for x and y, more or less.
Once we know the angle of the line connecting the two large bodies, we will use it to rotate back. We will rotate back and additional 90 degrees to position the line vertically on the drawing canvas.
As we are working in the Ecliptic, we need to consider only a 2D rotation: $$x_2 = x_1 cos(\alpha) - y_1 sin(\alpha),\\y_2 = x_1 sin(\alpha) + y_1 cos(\alpha).$$ The implementation of these calculations is detailed in the section Getting Oriented in 2D.
Pluto's Dance
Now that we have our rotating frame of reference, let's continue playing.
In the section on the Galilean Moons we encountered a first example of orbital resonance. Another example of this phenomenon can be found in the orbits of the planets Neptune and Pluto, with a period ratio of 3-to-2, as mentioned in Dr. Renu Malhotra's talk on The Search for Planet Nine.
Observed from a fixed perspective, this does not look very spectacular. However, using a rotating perspective, Pluto's orbit takes a much more interesting shape:
We see that the position of Neptune is fixed on the drawing canvas. This is the effect of our rotating frame. We are basically following the rotation of Neptune around the sun, by rotating our camera with the same angular velocity. Consequently, Neptune appears fixed on our drawing canvas.
Looking at the night sky on Neptune, meticulously measuring Pluto's position night after night, one can imagine its orbit going back and forth, much like in the dancing procession of Echternach. At least if you take enough time for your observations.
Playing with Lagrange and Lissajous
Back to the Lagrange points ...
We will put a 1 kg mass in the vicinity of L4, precisely on Jupiter's orbit, with the same orbital velocity, but at 66 degrees ahead of the orbit, rather than 60 degrees, where the Lagrange point is.
There are a few things to explain here.
First of all, we see that the position of Jupiter is fixed on the drawing canvas, as was the case with Neptune before, except for some vertical movement. This is due to the eccentricity of Jupiter's orbit: while orbiting the Sun, the distance between Jupiter and the Sun varies slightly.
Secondly, we see our fictitious body following an elliptical, fairly stable orbit around L4. The orbit as a whole is shifting back and forth along Jupiter's orbit. We also see that the total area covered by the orbit of our body is shaped along Jupiter's orbit as well.
Finally, we see that the orbit period is rather long, something in the order of around 12 years. Probably, the orbit is in 1-to-1 resonance with Jupiter's orbit.
Next, we will put our 1 kg mass exactly at 60 degrees, but 0.5% farther away from the Sun. As before, the orbital velocity is kept the same, giving the mass a kind of deficit in velocity for a nice elliptical orbit around the Sun.
Again, the situation is dynamic, but seems to be stable. This time, the frequency of oscillation in the radial direction is much higher than the frequency in the direction of the orbit. As a composition of two orthogonal oscillations with different frequency, we see a Lissajous-like curve coming to view.
We assume the radial oscillation is due to the radial oscillation of Jupiter itself, as it is consistent in many of these experiments. The oscillation in the direction of the orbit will be due to the deficit in velocity.
But hey, nothing beats the real thing ...
Trojan Asteroids
The Trojan asteroids are a large group of asteroids that share planet Jupiter's orbit around the Sun. They are orbiting one of Jupiter's two stable Lagrange points, L4 or L5. Shortly after the discovery of the first few of these asteroids, and inspired by the legend of the Trojan war, they were divided into a Greek camp, at L4, ahead of Jupiter in orbit, and a Trojan camp, at L5, behind of Jupiter in orbit.
The first Trojan asteriod to be discovered, in 1906, was Achilles, the fourth largest known so far, and part of the Greek camp. The largest three known until now are Hektor, Patroclus and Agamemnon.
Agamemnon is the third largest Trojan asteroid, with a diameter of approximately 168 kilometers. It was discovered on 19 March 1919, as part of the Greek camp. However, that is a historical anomaly, because it was discovered before the naming convention for Trojan asteroids was in place.
Using its identifier 911, we obtained its position on midnight 19-07-2019 from the Jet Propulsion Laboratory Horizons System as approximately: $$(6.62428\,10^{14} m, -4.22598\,10^{14} m, -5.73893\,10^{13} m),$$ and its velocity as approximately: $$(6.06120\,10^{6} \frac {m}{s}, 1.04530\,10^{7} \frac {m}{s}, 4.77569\,10^{6} \frac {m}{s}).$$ When we add this body to the simulation, we observe yet another variation in orbit shape.
Nowadays, Agamemnon is suspected to be a binary asteroid, meaning it has a (much) smaller companion. However, we will leave that simulation for another post.
Horsing around the Solar System
One thing that came up while playing with fictitious bodies orbiting Jupiter's Lagrange points is that, with specific starting conditions, a body may jump from an orbit around L4 to an orbit around L5 over the far end, away from Jupiter, kind of orbiting both points alternatingly. The starting conditions are delicate in terms of orbit stability, and the easiest way to reproduce this phenomenon is to start near L3, with small deviations in distance and velocity.
Consider this simulation, over a 400 year period, with a 1 kg body on Jupiter's opposite side of the Sun, at a distance increased by 0.1% and with an orbital velocity increased by 1.4%.
Seen from the rotating frame perspective, the body falls away in the direction of L4, oscillating, only to round L4, return, cross over, and eventually fall away in the direction of L5, where the same thing happens. After 400 years of simulation, the body has rounded L4 and L5 and is back where it started.
Some investigation teaches us that this phenomenon is known as a horseshoe orbit.
Although not shown in the animated gif, the simulation continues as follows:
- After 1,200 years, the body is not able to cross over from L5 to L4 again and falls back to L5, as if it lacks kinetic energy.
- After 1,500 years, coming back from L5, the body does cross over to L4 again, as if the lack in kinetic energy was due to an increase in potential energy, and the energy conversion inverted along the way.
- After 3,000 years, after the body has rounded both L4 and L5 eight times, the stability of the orbit is broken, due to a close encounter with Jupiter, and the body starts rotating the Sun like crazy in its own spirographic way.
The phenomenon has been observed with real asteroids though, around the Sun-Earth Lagrange points for example, with asteroids 2010 SO16 and 2002 AA29. This suggests stable horseshoe orbits, excuse the wording, exist, at least for some configurations.
The Mischievous Bee-Zed
We recently bumped into an article on the Bee-Zed asteroid, designation 2015 BZ509.
The asteroid follows a retrograde orbit around the sun, co-orbital with Jupiter. From the rotating frame perspective, the orbit is shaped like a trisectrix. Adding that body to the simulation was straightforward, with an immediate result:
One can observe the asteroid rotating clockwise, with a orbital period of around 12 earth years, perfectly in orbital resonance with Jupiter. However, a single orbit consists of a longer loop (in diameter) around the sun, followed by a second, shorter loop around the sun, before returning to the starting configuration, slightly shifted. Running a (much) longer simulation, the shorter loop gradually swings to the right-hand side, only to swing back to the left-hand side.
The asteroid is also referred to as Ka'epaoka'awela. This is Hawaiian for "the mischievous opposite-moving companion of Jupiter".
Conclusions
The problem of planetary dynamics is an N-body problem. Consequently, with N>2, the solar system is a chaotic system, where long term predictions or simulations will be inaccurate.
The simulation approach taken in this post is a very simplistic one, using nothing more than Newton's law of universal gravitation and Newton's laws of motion, to evolve a system in discrete steps using first-order relations only. It does not adapt the simulation step size to the circumstances and ignores higher-order effects. Relative errors will be larger for strongly curved orbits and generally speaking errors will accumulate and get magnified over time.
Its implementation however is very compact, with, given the right building blocks, 10 lines of code for the calculation of the net gravitation force, and 6 lines of code for the whole of a simulation step. And it provides an acceptable prediction accuracy over short (planetary) time scales. As such, it is a suitable tool to play with some non-trivial aspects of planetary motion:
- It resulted in a visually accurate simulation of the planetary orbits in the solar system over the time frame in which Pluto completes its orbit.
- It visually demonstrated the orbital resonance of some of Jupiter's large moons.
- It allowed us to observe the symmetric nature of gravitation by demonstrating the perturbations of the Sun's orbit induced by the planets orbiting it.
- It could be used to demonstrate asteroids orbiting Jupiter's Lagrange points.
- And it could be used to demonstrate asteroids in a horseshoe orbit, at least for some length in time.
Code Fragments
In what follows, we will detail the Java implementation of some of the subjects in more depth.
Vector Representation and Calculations
We start the implementation with an immutable representation of 3D vectors. Immutable in principle, because we will declare the properties public to avoid the verbosity of getters and setters.
public class Vector {
public double x;
public double y;
public double z;
public Vector(double x, double y, double z) {
this.x = x;
this.y = y;
this.z = z;
}
public Vector() {
this(0, 0, 0);
}
}
We add the basic vector functions and operations to the class, because we anticipate we will need those as we go along:
public class Vector {
...
public double norm() {
return Math.sqrt(x*x+y*y+z*z);
}
public Vector add(Vector v) {
return new Vector(x+v.x, y+v.y, z+v.z);
}
public Vector sub(Vector v) {
return add(v.mul(-1));
}
public Vector mul(double r) {
return new Vector(r*x, r*y, r*z);
}
public Vector div(double r) {
return mul(1/r);
}
}
Note that none of these operations will modify the vector acted upon.
Being an immutable type, potentially destructive operations will create new instances of the class.
This technique greatly simplifies using the vector class.
Celestial Body Representation
We continue the implementation with the representation of celestial bodies and the function to calculate the net force vector.
We do not need a lot of information on a celestial body. The implementation is restricted to hold the name of the body, if any, the mass, position, velocity and the force.
public class Body {
public String name;
public double mass;
public Vector position;
public Vector velocity;
public Vector force;
public Body(String name, double mass, Vector position, Vector velocity) {
this.name = name;
this.mass = mass;
this.position = position;
this.velocity = velocity;
}
}
The position, velocity and force are vector values.
The force can be seen as the force that is currently exerted on the body.
It is used for temporary storage during the simulation.
Note that we will consider all celestial bodies to be point masses for the whole of this exercise.
We could have created all required bodies dynamically, but opted for an explicit implementation of the bodies at play. Consider the implementation of our home planet:
public class EarthBody extends Body {
public EarthBody() {
super(
"Earth", 5.97219e24,
new Vector(6.618496041458324e10, -1.368702026239706e11, 6.333028190493584e6),
new Vector(2.632453093013447e4, 1.284679136106616e4, 2.897601110483095e-1)
);
}
}
We can calculate the force exerted on a body as described in the section on Newton's Law of Universal Gravitation.
public class EarthBody extends Body {
...
public Vector calculateForce(List<Body> bodies) {
Vector forceVector = new Vector();
for (Body body : bodies) {
if (this!=body) {
Vector deltaVector = body.position.sub(position);
double delta = deltaVector.norm();
Vector deltaUnitVector = deltaVector.div(delta);
double force = Constants.G*mass*body.mass/(delta*delta);
forceVector = forceVector.add(deltaUnitVector.mul(force));
}
}
return forceVector;
}
}
The function initialises the net force vector as the origin.
It then iterates over all known bodies, a collection passed as a parameter.
For each body, it calculates the strength of the force exerted by that particular body on this body,
uses it to multiple the unit vector in the direction towards that body and adds the resulting vector to the net force vector.
This way, all individual forces are accumulated in the net force vector.
It skips the subject body itself to avoid division by zero.
The Simulation Step
Now we are ready to write the code for the simulation, as a function that performs a single step in the simulation, representing an elapse time of dt, on all bodies that participate in the simulation:
public void simulationStep(List<Body> bodies, long dt) {
// calculate forces
for (Body body : bodies) {
body.force = body.calculateForce(bodies);
}
// move bodies
for (Body body : bodies) {
Vector acceleration = body.force.div(body.mass);
body.velocity = body.velocity.add(acceleration.mul(dt));
body.position = body.position.add(body.velocity.mul(dt));
}
}
We are sure you are pleasantly surprised by the brevity and simplicity of this code fragment, just as we were.
The function first updates the force exerted on all bodies, and then updates the position and velocity of all bodies. It calculates the acceleration, the change in velocity and the change in position of each body accordingly, using the techniques described in the section on Newton's Laws of Motion.
Getting Oriented in 2D
This section provides the implementation of the building blocks to use a rotating frame of reference.
The calculation of the angle of a 2D vector against the X-axis is implemented as a static utility function:
public static double angle(double x, double y) {
return Math.PI
-Math.PI/2*(1+Math.signum(x))*(1-Math.signum(y*y))
-Math.PI/4*(2+Math.signum(x))*Math.signum(y)
-Math.signum(x*y)*Math.atan((Math.abs(x)-Math.abs(y))/(Math.abs(x)+Math.abs(y)));
}
The 2D rotation of a vector is implemented as an operation on the 3D vector data type we already introduced, for convenience:
public class Vector {
...
public Vector rotate(double angle) {
return new Vector(
x*Math.cos(angle)-y*Math.sin(angle),
x*Math.sin(angle)+y*Math.cos(angle),
0
);
}
}
These two building blocks will allow us to compensate the rotation of bodies and keep the bodies fixed in our drawing canvas.
Winding Time Series
If you are new to the concept of a Fourier tranform, as stated earlier, take a look at the Visual Introduction to Fourier Transforms by 3Blue1Brown.
For our purpose, we will be using Apache Commons Math to implement the Fourier transform:
double[] data = new double[size]; // populate data[] with time series ... FastFourierTransformer transformer = new FastFourierTransformer(DftNormalization.STANDARD); Complex[] complex = transformer.transform(data, TransformType.FORWARD); // use real part (amplitude) of complex[] ...
The data array is populated with the deviations of the Sun's position along our X-axis. Note that the implementation provided by Apache Commons Math requires the size of the input array to be a power of two. For our purpose, we provided the deviation every 8 days, over a period of around 718 years, 2^15 values in total.
The output provided by the transformation in the complex array should be interpreted as follows:
- The index in the array represents the angular frequency of the data point relative to the size of the array. Therefore, divide the index by the size of the array. Invert this frequency to obtain the period. In our case, one grid point represents 8 days, so multiply the period by 8 to get the period in days. Divide the period by 365.25, the average number of days in a year, to get the period in year fractions.
- The real part of the complex number represents the amplitude for that frequency.
- The imaginary part of the complex number represents the phase for that frequency.
Jet Propulsion Laboratory to the Rescue
The mass and average distance to the sun of all planets can easily be found on the Internet. However, synchronous information on the position and velocity of the planets is a lot harder to find.
The Jet Propulsion Laboratory Horizons System comes to the rescue. It offers detailed information for a vast number of celestial bodies, accessible using different interfaces.
For the exercise at hand, we opted to use the mail interface. To use it, just send an e-mail to horizons@ssd.jpl.nasa.gov, with subject "job" and a body like:
!$$SOF EMAIL_ADDR='' START_TIME = '2019-Jul-19 00:00:00' STOP_TIME = '2019-Jul-19 00:00:01' TABLE_TYPE = 'Vector' REF_PLANE = 'Ecliptic' CENTER = '@010' COMMAND='399' !$$EOFHere, we specify the date and time for the orbit information, namely midnight 19-07-2019, and the coordinate system to be used, namely vector coordinates against the ecliptica. We specify the sun (identifier 010) as the origin of the vector space and request information for the earth (identifier 399). The identifier of a celestial body can easily be found using the Horizon's Search Function.
You will get a reply in a few seconds, with rich, structured information on the celestial body. For our exercise, we extract the mass of the body, and the position and velocity vectors at the requested date and time:
... Revised: July 31, 2013 Earth 399 GEOPHYSICAL PROPERTIES (revised Aug 15, 2018): Vol. Mean Radius (km) = 6371.01+-0.02 Mass x10^24 (kg)= 5.97219+-0.0006 ... $$SOE 2458683.500000000 = A.D. 2019-Jul-19 00:00:00.0000 TDB X = 6.618496041458324E+07 Y =-1.368702026239706E+08 Z = 6.333028190493584E+03 VX= 2.632453093013447E+01 VY= 1.284679136106616E+01 VZ= 2.897601110483095E-04 LT= 5.071260561131307E+02 RG= 1.520325668780014E+08 RR=-1.056016984269061E-01 $$EOE ...After a session of careful mailing and processing replies, our local database of celestial body characteristics is populated.
Note that the position and velocity information obtained are expressed in kilometers. We converted everything into meters before using the information in the simulation.
