# Further On An Ethereal Orrery

Last time we met[1] we spoke of my fellow students' and my interest in constructing a model of the motion of heavenly bodies using mathematical formulae in the place of brass. In particular we have sought to do so from first principals using Sir N-----'s law of universal gravitation, which states that the force attracting two bodies is proportional to the product of their masses divided by the square of the distance between them, and his laws of motion, which state that a body will remain at rest or in constant motion unless a force acts upon it, that it will be accelerated in the direction of that force at a rate proportional to its magnitude divided the body's mass and that a force acting upon it will be met with an equal force in the opposite direction.
Whilst Sir N----- showed that a pair of bodies traversed conic sections under gravity, being those curves that arise from the intersection of planes with cones, the general case of several bodies has proved utterly resistant to mathematical reckoning. We must therefore approximate the equations of motion and I shall now report on our first attempt at doing so.

### Approximating Equations Of Motion

Recall that we experimented with a variety of schemes by which we might do so for small steps in time, finally settling upon one whose errors grow in proportion to the sixth power of their lengths[2][3]. Specifically, given the coefficients
\begin{align*} c_0 &= \phantom{-}0.784513610477560\\ c_1 &= \phantom{-}0.235573213359357\\ c_2 &= -1.17767998417887\\ c_3 &= \phantom{-}1.31518632068391\\ c_4 &= \phantom{-}c_2\\ c_5 &= \phantom{-}c_1\\ c_6 &= \phantom{-}c_0 \end{align*}
defining an initial state of the motion at time $$t_i$$ as
\begin{align*} x_{i,0} &= x_i\\ v_{i,0} &= v_i\\ a_{i,0} &= a_i \end{align*}
and the step in time as
$\Delta t_i = t_{i+1} - t_i$
then updating it with the sub-steps
\begin{align*} \Delta t_{i,j} &= c_j \Delta t_i\\ x_{i,j+1} &= x_{i,j} + v_{i,j} \Delta t_{i,j} + \tfrac12 a_{i,j} \Delta t_{i,j}^2\\ v_{i,j+1} &= v_{i,j} + \tfrac12 \left(a_{i,j}+a_{i,j+1}\right) \Delta t_{i,j} \end{align*}
which require that would should be able to figure the acceleration $$a_{i,j+1}$$ from the position $$x_{i,j+1}$$, result in the state at time $$t_{i+1}$$ being well approximated by
\begin{align*} x_{i+1} &= x_{i,7}\\ v_{i+1} &= v_{i,7}\\ a_{i+1} &= a_{i,7} \end{align*}
We tested this scheme with an exponential motion whose position at time $$t$$ was defined by the function
$x(t) = e^t$
which has a velocity $$v(t)$$ and an acceleration $$a(t)$$ defined by the derivatives
\begin{align*} v(t) & = \frac{\mathrm{d}x(t)}{\mathrm{d}t} = e^t\\ a(t) & = \frac{\mathrm{d}v(t)}{\mathrm{d}t} = e^t \end{align*}
so that we could approximate $$a_{i+1}$$ with $$x_{i+1}$$ and we assembled a deck to plot the positions yielded by our chosen integration algorithm in red against the exact positions in black using time steps of one tenth, presented here as deck 1.

Deck 1: The Exact And Approximate Motions

Since the approximate and exact positions are indistinguishable at its scale we put together another deck to plot the difference between them at a scale that was several orders of magnitude finer

Deck 2: The Difference Between Them

which most certainly convinced us that it would be an appropriate engine for our ethereal orrery.

### The Gravity Of The Situation

Given a pair of particles $$P_i$$ and $$P_j$$, having masses $$m_i$$ and $$m_j$$ and vector positions $$\mathbf{x}_i$$ and $$\mathbf{x}_j$$, Sir N-----'s law of gravitation states that the first will be attracted toward the second with a force of
$\mathbf{F}_{i,j} = \frac{G m_i m_j}{d_{i,j}^2} \times \mathbf{u}_{j,i}$
where $$d_{i,j}$$ is the distance between them, $$\mathbf{u}_{j,i}$$ is a vector of length one in the direction of the first toward the second and $$G$$ is known as the gravitational constant.
Noting that the distance between a pair of points is simply the magnitude of their difference and that dividing a vector by its magnitude yields a vector with magnitude one pointing in the same direction, we have
$\mathbf{F}_{i,j} = \frac{G m_i m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^2} \times \frac{\mathbf{x}_j - \mathbf{x}_i}{\left|\mathbf{x}_j - \mathbf{x}_i\right|} = \frac{G m_i m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^3} \times \left(\mathbf{x}_j - \mathbf{x}_i\right)$
Similarly, the force attracting the second toward the first will be
$\mathbf{F}_{j,i} = \frac{G m_j m_i}{\left|\mathbf{x}_i - \mathbf{x}_j\right|^3} \times \left(\mathbf{x}_i - \mathbf{x}_j\right) = \frac{G m_i m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^3} \times -\left(\mathbf{x}_j - \mathbf{x}_i\right) = -\mathbf{F}_{i,j}$
which is wholly consistent with his third law of motion. Furthermore, the second law of motion implies that the acceleration of the first particle toward the second will be
$\mathbf{a}_{i,j} = \frac{\alpha \mathbf{F}_{i,j}}{m_i} = \frac{\alpha \, G m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^3} \times \left(\mathbf{x}_j - \mathbf{x}_i\right)$
for some positive constant $$\alpha$$, and consequently, if there are several particles, the acceleration of the $$i^\mathrm{th}$$ will be the vector sum of the accelerations caused by the others, so that
$\mathbf{a}_i = \sum_{j \neq i} \mathbf{a}_{i,j} = \alpha \, G \sum_{j \neq i} \frac{m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^3} \times \left(\mathbf{x}_j - \mathbf{x}_i\right)$
where $$\sum$$ is the summation sign.

Under the, quite frankly revolutionary, French system of units which measures lengths in metres, masses in kilograms, time in seconds and forces in newtons, $$G$$ is roughly equal to $$6.67408 \times 10^{-11}$$ and $$\alpha$$ is exactly equal to one.
Whilst these units are supremely well suited to earthly measurements they are rather insignificant at the heavenly scale and so my fellow students and I decided to instead build our orrery with lengths measured in terametres, masses in hellagrams, time in megaseconds and forces in yottanewtons.
The prefixes kilo, mega, tera, yotta and hella represent multiplicative factors of ten to the power of three, six, twelve, twenty four and twenty seven respectively so that, using $$\mathrm{m}$$ to represent one metre, $$\mathrm{g}$$ one gram, $$\mathrm{s}$$ one second and $$\mathrm{N}$$ one newton, with the symbols $$\mathrm{k}$$ representing kilo, $$\mathrm{T}$$ tera, $$\mathrm{M}$$ mega, $$\mathrm{Y}$$ yotta and $$\mathrm{H}$$ hella, we have
\begin{align*} 1 \, \mathrm{Tm} &= 10^{12} \, \mathrm{m}\\ 1 \, \mathrm{kg} &= 10^{3} \, \mathrm{g}\\ 1 \, \mathrm{Hg} &= 10^{27} \, \mathrm{g} = 10^{24} \, \mathrm{kg}\\ 1 \, \mathrm{Ms} &= 10^6 \, \mathrm{s}\\ 1 \, \mathrm{YN} &= 10^{24} \, \mathrm{N} \end{align*}
The benefits of these units are that, very roughly, the distance from the Earth to the Sun is one seventh of a terametre, the mass of the Earth is six hellagrams and one megasecond is eleven and a half days, so we shall not have to reckon with very large or very small numbers.
Serendipitously, the constants $$\alpha$$ and $$G$$ remain unchanged since the unit of force is equal to the unit of mass times the unit of acceleration, which is in turn equal to the unit of length divided by the square of the unit of time, and so
\begin{align*} \mathrm{N} &= \mathrm{kg} \times \frac{\mathrm{m}}{\mathrm{s}^2}\\ 10^{24} \, \mathrm{N} &= 10^{24} \, \mathrm{kg} \times \frac{10^{12} \, \mathrm{m}}{10^{12} \, \mathrm{s}^2}\\ 10^{24} \, \mathrm{N} &= 10^{24} \, \mathrm{kg} \times \frac{10^{12} \, \mathrm{m}}{\left(10^6 \, \mathrm{s}\right)^2}\\ \mathrm{YN} &= \mathrm{Hg} \times \frac{\mathrm{Tm}}{\left(\mathrm{Ms}\right)^2} \end{align*}

### The Point Of The Matter

In constructing our orrery, my fellow students and I made two significant simplifying assumptions; firstly, that its bodies move upon a plane and, secondly, that they are point masses whose only defining properties are mass, position and velocity. We may therefore fully specify a body with the five figures $$m$$, $$x$$, $$y$$, $$v_x$$ and $$v_y$$ representing its mass, its position on the $$x$$ axis, its position on the $$y$$ axis, its velocity along the $$x$$ axis and its velocity along the $$y$$ axis respectively.
The accelerations function furnished by listing one adds the current accelerations in the directions of the $$x$$ and $$y$$ axes to records defining the masses, positions and velocities of bodies, stored in the array p, exploiting Sir N-----'s third law to halve the number of calculations required to do so.

Listing 1: Figuring The Accelerations
var G = 6.67408e-11;

function accelerations(p) {
var n = p.length;
var i, j, dx, dy, d2, g, ax, ay;

for(i=0;i<n;++i) {
p[i].ax = 0;
p[i].ay = 0;
}

for(i=1;i<n;++i) {
for(j=0;j<i;++j) {
dx = p[j].x-p[i].x;     dy = p[j].y-p[i].y;
d2 = dx*dx+dy*dy;       g  = G/(d2*Math.sqrt(d2));
ax = g*(p[j].x-p[i].x); ay = g*(p[j].y-p[i].y);

p[i].ax += p[j].m*ax;   p[i].ay += p[j].m*ay;
p[j].ax -= p[i].m*ax;   p[j].ay -= p[i].m*ay;
}
}
}


The substep and step functions given in listing 2 apply our chosen scheme to those bodies, updating their velocities in two steps; halfway with the accelerations at the present time and halfway after a time of dt has passed.

Listing 2: Figuring The Motion
function substep(p, dt) {
var n = p.length;
var i;

for(i=0;i<n;++i) {
p[i].x += p[i].vx*dt + 0.5*p[i].ax*dt*dt;
p[i].y += p[i].vy*dt + 0.5*p[i].ay*dt*dt;

p[i].vx += 0.5*p[i].ax*dt;
p[i].vy += 0.5*p[i].ay*dt;
}

accelerations(p);

for(i=0;i<n;++i) {
p[i].vx += 0.5*p[i].ax*dt;
p[i].vy += 0.5*p[i].ay*dt;
}
}

var C = [0.784513610477560, 0.235573213359357, -1.17767998417887, 1.31518632068391];

function step(p, dt) {
var i;
for(i=0;i<7;++i) substep(p, dt*C[Math.min(i, 6-i)]);
}


To put these functions to the test we decided to model the motion of the Earth about the Sun and to do so we have estimated the figures specifying the initial state of the Earth with
\begin{align*} m &= 5.97 \, \mathrm{Hg}\\ x &= 0.152 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.0293 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
and those of the Sun with
\begin{align*} m &= 1,990,000 \, \mathrm{Hg}\\ x &= 0 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
so that all positions and velocities shall be measured relative to those initially held by the Sun.
Now, one year is approximately equal to
$1 \, \mathrm{yr} \approx 365.26 \times 24 \times 60 \times 60 \, \mathrm{s} = 31,558,464 \, \mathrm{s} = 31.558464 \, \mathrm{Ms}$
and so we should expect our model Earth to make three orbits of our model Sun in slightly less than $$94.7 \, \mathrm{Ms}$$.
We therefore assembled deck three to plot the progress of our model for nine hundred and forty seven steps of $$0.1 \, \mathrm{Ms}$$ apiece.

Deck 3: The Sun And The Earth

I must confess that my fellow students and I were most relieved that not only did our Earth follow a very nearly circular path about our Sun but that it did so slightly more than three times!
Furthermore, changing the initial value of Earth's $$v_y$$ to $$0.0294$$ results in it doing so slightly less than three times, as you may confirm for yourself, which was sufficient to convince us that the discrepancy is a consequence of the inexactitude of our initial values for their states.

### The Inner Planets

Bolstered by our initial success, we next built a model of the inner planets by adding Mercury with the initial state
\begin{align*} m &= 0.330 \, \mathrm{Hg}\\ x &= 0.0698 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.0387 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
Venus with
\begin{align*} m &= 4.87 \, \mathrm{Hg}\\ x &= 0.109 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.0348 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
and Mars with
\begin{align*} m &= 0.642 \, \mathrm{Hg}\\ x &= 0.249 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.0219 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
Deck 4 plots the motion of the inner planets, with Mercury in black, Venus in grey, Earth again in green and Mars in red.

Deck 4: The Inner Planets

We were once more most satisfied with the behaviour of our orrery; not least because it shows Mercury's orbit to be highly eccentric, being clearly elliptical rather than very nearly circular, just as it is in the heavens.

### The Outer Planets

Turning to the outer planets we gave Jupiter an initial state of
\begin{align*} m &= 1900 \, \mathrm{Hg}\\ x &= 0.817 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.0125 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
Saturn one of
\begin{align*} m &= 568 \, \mathrm{Hg}\\ x &= 1.51 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.00915 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
Uranus
\begin{align*} m &= 86.8 \, \mathrm{Hg}\\ x &= 3.01 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.00649 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
and finally Neptune
\begin{align*} m &= 102 \, \mathrm{Hg}\\ x &= 4.54 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.00538 \, \mathrm{Tm}/\mathrm{Ms} \end{align*}
Deck 5 plots their motion in isolation from the inner planets, with Jupiter in light brown, Saturn in reddish brown, Uranus in cyan and Neptune in blue.

Deck 5: The Outer Planets

Note that, since the action of the Sun's gravity at the immense distances at which these bodies reside is so weak, they proceed quite ponderously and we are consequently able to figure their motion with the relatively large time step of $$5 \, \mathrm{Ms}$$.
When combining these parts to fashion a whole, the time step must be chosen carefully; too small and it will take a tremendous effort to reckon the motion of the outer planets, too large and the motion of the inner planets will be wholly inaccurate. My fellow students and I put together deck 6 to do so, thereby creating our ethereal orrery, the action of which you may stop by pressing the halt button, and after some tinkering settled upon a time step of $$0.5 \, \mathrm{Ms}$$.

Deck 6: The Solar System

Unfortunately, given the quite staggering difference in scale between the distances of the inner and the outer planets from the Sun, it's rather difficult to discern the motions of the former.
To rectify this, we recast the rectangular coordinates of a body $$(x, y)$$ as polar coordinates $$(r, \theta)$$, representing its distance from the origin $$(0,0)$$ and the angle between the vector from the origin to it and the $$x$$ axis respectively, which are related by the formulae
\begin{align*} x &= r \times \cos \theta\\ y &= r \times \sin \theta \end{align*}
so that
\begin{align*} x^2 + y^2 &= r^2 \times \cos^2 \theta + r^2 \times \sin^2 \theta = r^2 \times \left(\cos^2 \theta + \sin^2 \theta\right) = r^2\\ \frac{y}{x} &= \frac{r \times \sin \theta}{r \times \cos \theta} = \frac{\sin \theta}{\cos \theta} = \tan \theta \end{align*}
To exaggerate the distances from the Sun of the inner planets and to depreciate those of the outer planets we transformed them back to rectangular coordinates using the square root of $$r$$ with
\begin{align*} x^\prime &= \sqrt{r} \times \cos \theta\\ y^\prime &= \sqrt{r} \times \sin \theta \end{align*}
and the results of doing so are plotted by deck 7.

Deck 7: Upon The Transformed Coordinates

In addition to making the motions of the inner planets much plainer, this makes it most evident that the Sun too is moved by Sir N-----'s laws!

Whilst my fellow students and I are fully satisfied with our orrery we cannot help but wonder what its engine might make of hypothetical arrangements of bodies, but unfortunately our curiosity shall have to remain unsatisfied until we have time enough away from our studies to satiate it.
$$\Box$$

### References

[1] On An Ethereal Orrery, www.thusspakeak.com, 2019.

[2] Yoshida, H., Construction of higher order symplectic integrators, Physics Letters A, Vol. 150, No. 5, 6, 7, Elsevier, 1990.

[3] Hut, P. & Makino, J., The 2-Body Problem: Higher-Order Integrators, www.artcompsci.org, 2007.