# Further Still On An Ethereal Orrery

Recently, my fellow students and I constructed a mathematical orrery which modelled the motion of heavenly bodies employing Sir N-----'s laws of gravitation and motion, rather than clockwork, as its engine[1][2]. Those laws state that bodies are attracted toward each other with a force proportional to the product of their masses divided by the square of the distance between them, that a body will remain at rest or in constant motion unless a force acts upon it, that if a force acts upon it then it will be accelerated in the direction of that force at a rate proportional to its strength divided by its mass and that, if so, it will reciprocate with an opposing force of equal strength.
Its operation was most satisfactory, which set us to wondering whether we might use its engine to investigate the motions of entirely hypothetical arrangements of heavenly bodies and I should now like to report upon our progress in doing so.

You will no doubt recall that the motion of more than two bodies under gravity defies mathematical reckoning and so we were compelled to approximate the motion of the planets with a scheme for solving their equations of motion that admitted errors on the order of the sixth power of the magnitude of the steps in time by which we advanced them[3][4], given again in listing 1.

Listing 1: Our Orrery's Engine
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;
}
}
}

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)]);
}


Using the units of terameter, hellagram, megasecond and yottanewton given by
\begin{align*} 1 \, \mathrm{Tm} &= 10^{12} \, \mathrm{m}\\ 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*}
for distance, mass, time and force respectively, we finally used this to model the solar system upon a quadratic scale of distance from the Sun's initial position, as shown by deck 1 which plots the Sun in orange, Mercury in black, Venus in grey, Earth in green, Mars in red, Jupiter in light brown, Saturn in reddish brown, Uranus in cyan and Neptune in blue and which you can stop running by pressing the halt button.

Deck 1: Our Ethereal Orrery

### Uniformly Distributed Bodies

We commenced our experimentation by modelling systems of bodies having uniformly distributed masses, positions and velocities, as demonstrated by deck 2.

Deck 2: Uniformly Distributed Bodies

I must confess that my fellow students and I were quite taken aback by the outcome; our ethereal orrery yielded supremely regular motion of the heavenly bodies, but applying its engine to such systems casts many of them asunder at tremendous velocities!

To resolve this apparently contradictory behaviour we returned to Sir N-----'s laws which, expressed in the language of mathematics, state that the force that a body $$P_i$$, located at a vector position $$\mathbf{x}_i$$ and having a mass of $$m_i$$, is subjected to by another $$P_j$$, having a position $$\mathbf{x}_j$$ and mass $$m_j$$, is given by
$\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)$
for some positive constant $$G$$, where the vertical bars stand for the length of the vector between them, that this force is reciprocated with
$\mathbf{F}_{j,i} = -\mathbf{F}_{i,j}$
and that they consequently undergo accelerations of
\begin{align*} \mathbf{a}_{i,j} &= \phantom{-}\frac{\alpha \mathbf{F}_{i,j}}{m_i} = \phantom{-}\frac{\alpha \, G m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^3} \times \left(\mathbf{x}_j - \mathbf{x}_i\right)\\ \mathbf{a}_{j,i} &= -\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) \end{align*}
for some positive constant $$\alpha$$.
That the distance between two bodies grows ever smaller as they approach each other is clearly tautological and therefore $$\left|\mathbf{x}_j - \mathbf{x}_i\right|$$ might be arbitrarily close, or even equal, to zero. Since the magnitude of the acceleration $$\mathbf{a}_{i,j}$$ is given by
$\left|\mathbf{a}_{i,j}\right| = \frac{\alpha \, G m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^3} \times \left|\mathbf{x}_j - \mathbf{x}_i\right| = \frac{\alpha \, G m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^2}$
it is without bound which obviously, in retrospect, implies that $$P_i$$ might reach an exceedingly large velocity.

The reason why our orrery's engine admits such behaviour when it most evidently does not occur in the heavens is a consequence of our treating bodies as infinitesimal points that are entirely defined by their masses, their positions and their velocities whereas, in reality, a pair of them will collide long before their centres are coincident.
To rectify this my fellow students and I gave its bodies radii so that it might determine when worlds collide and to simplify the reckoning of the outcomes of such collisions we decided to make them perfectly inelastic so that a pair of bodies stick together rather than rebound if the distance between their centres is less than the sum of their radii.
To figure the velocity of the combined body we employed Sir N-----'s laws once again with
\begin{align*} \mathbf{F}_{j,i} &= -\mathbf{F}_{i,j}\\ \frac{m_j}{\alpha} \times \mathbf{a}_{j,i} &= -\frac{m_i}{\alpha} \times \mathbf{a}_{i,j}\\ \frac{m_j}{\alpha} \times \frac{\mathrm{d}}{\mathrm{d}t}\mathbf{v}_{j,i} &= -\frac{m_i}{\alpha} \times \frac{\mathrm{d}}{\mathrm{d}t}\mathbf{v}_{i,j} \end{align*}
where $$\frac{\mathrm{d}}{\mathrm{d}t}$$ represents the derivative with respect to time of the term following it, which is to say its rate of change. Rearranging yields
\begin{align*} \frac{1}{\alpha} \times \frac{\mathrm{d}}{\mathrm{d}t} \left(m_j \mathbf{v}_{j,i} + m_i \mathbf{v}_{i,j}\right) &= 0\\ \frac{\mathrm{d}}{\mathrm{d}t} \left(m_j \mathbf{v}_{j,i} + m_i \mathbf{v}_{i,j}\right) &= 0\\ m_j \mathbf{v}_{j,i} + m_i \mathbf{v}_{i,j} &= c \end{align*}
for some constant $$c$$. At the instant of collision this constant remains unchanged and so the velocity $$\mathbf{v}^\prime$$ of the combined body is given by
\begin{align*} \left(m_j + m_i\right) \mathbf{v}^\prime &= m_j \mathbf{v}_{j,i} + m_i \mathbf{v}_{i,j}\\ \mathbf{v}^\prime &= \frac{m_j \mathbf{v}_{j,i} + m_i \mathbf{v}_{i,j}}{m_j + m_i} \end{align*}
The collisions function given in listing 2 merges colliding bodies, giving their unions such velocities and indicating whether or not any collisions took place.

Listing 2: Figuring The Collisions
function collisions(p) {
var c = 0;
var n = p.length;
var i, j, r, dx, dy;

i = 1;
while(i<n) {
j = 0;
while(i<n && j<i) {
r = p[i].r+p[j].r;
dx = p[i].x-p[j].x;
dy = p[i].y-p[j].y;
if(dx*dx+dy*dy<=r*r) {
p[j].r = Math.pow(Math.pow(p[i].r, 3) + Math.pow(p[j].r, 3), 1/3);
p[j].vx = (p[i].m*p[i].vx + p[j].m*p[j].vx) / (p[i].m+p[j].m);
p[j].vy = (p[i].m*p[i].vy + p[j].m*p[j].vy) / (p[i].m+p[j].m);
p[j].m += p[i].m;
p[i] = p[--n];
p.length = n;
++c;
}
else {
++j;
}
}
++i;
}
return c>0;
}


Listing 3 shows how we adapted the substep function to employ the collisions function, recalculating the accelerations that our orrery's bodies are subject to in the event of collisions

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;
}

if(collisions(p)) accelerations(p);
}


and deck 3 demonstrates how our new engine processes uniformly distributed matter, with the size of the tick marks representing the size of the bodies and the width of their lines their mass.

Deck 3: Our New Engine

This certainly behaves in a much more plausible manner, frequently resulting in smaller bodies orbiting larger ones.

### Orbiting Bodies

My fellow students and I were curious as to how our engine might treat many relatively small bodies in approximately the same orbit about a star. We therefore constructed deck 4 to place such bodies roughly in Earth's orbit about the Sun, with the mass of each being between one twentieth and one half that of the Earth and their radii being five hundred times that that they should have had had they its density, so as to significantly increase the chances of collisions.

Deck 4: Orbiting Bodies

We were certainly impressed by the outcome and it set us to wondering what effect other heavenly bodies might have upon their accretion. To that end put together deck 5 which adds Venus, Mars and Jupiter to represent the former whilst putting the latter much closer to Earth's orbit and reducing their masses.

Deck 5: With Venus, Mars And Jupiter

Whilst we have no doubt that these orbital models paint wholly fictional pictures of planetary formation, we are confident that they are at least illustrative of the process.

We are not quite finished exploring the behaviour of our orrery's engine, but for now our studies must take precedence. Of course I shall be sure to report upon our future findings once we have had the time to do so.
$$\Box$$

### References

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

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

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

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