Phyllotaxis and Fibonacci

Lagrange points

Imagine you were a hamster in a hamster ball living on a hilly surface.  Base camp is surrounded on all sides by high summits, but you have a powerful catapult there that can fire you to the top of any of them.  Once fired you can influence your trajectory, but it's hard work and you don't have much energy in your little legs.  Suppose you know where you want to end up beyond the hills.  What's the best strategy for getting there?


A smart plan would be to aim the catapult at the top of one of the surrounding hills and supply it with exactly the energy needed to get you to the summit.  Then you should arrive at the summit with zero velocity, at which point only a tiny amount of energy is required to completely determine the direction you take from there, meaning that you can take your time to reorient and plan the remainder of your journey.

Planning a journey in a rocket ship is a bit like this.  The amount of rocket fuel needed to escape Earth's potential well is an exponential function of the escape velocity, so, like the hamster, you can't afford to take very much fuel for controlling your trajectory once you're beyond Earth's pull.

[To see this suppose that $m_p$ is the mass of the payload, $m_0$ is the takeoff mass, $m$ is the total mass including the remaining fuel, $u$ is the relative speed at which burnt fuel is ejected, and $v$ is the speed of the craft, then
$$
\begin{align}
v &= \int \frac{F}{m} dt \\&= - \int  u \frac{dm}{dt}\frac{dt}{m} \\
&= -u \int_{m_0}^{m_p} \frac{dm}{m} \\
&= u \, ln\left(\frac{m_0}{m_p}\right)
\end{align}
$$
So each time you want to add u to the speed the $\frac{\text{takeoff mass}}{\text{payload mass}}$ ratio needs to go up by 2.718.]


The hamster ball analogy seems to be inappropriate since, in terms of gravitational potential, it doesn't look like Earth is surrounded by summits.  Instead it looks like we are stuck in a deep pit of despair and we need to climb indefinitely to proceed in any direction!  This is because our gravitational potential is $-\frac{GM_e}{r}$ where $M_e$ is the mass of the Earth and $r$ the radial distance from it.  (Note that the moon also contributes a term like this, but that causes a saddle point to appear between us and it, not a local maxima.)

But this is only because we are looking at it in the wrong way.  The Lagrangian for a spaceship of unit mass is as follows (where $M_e$ and $M_m$ are the masses of the Earth and moon, and $r_e$ and $r_m$ are the distances to the centres of the Earth and moon, respectively)
$$
\begin{align}
L &= T - V \\
&= \frac{1}{2}\left(\dot{r}^2 + r^2\dot{\theta}^2 \right) - V\\
&= \frac{1}{2}\left(\dot{r}^2 + r^2\dot{\theta}^2 \right) + \frac{GM_e}{r_e} +\frac{GM_m}{r_m}
\end{align}
$$
What if we rewrite this for the rotating frame of reference where the origin is the centre of mass of the Earth and moon and the two bodies are fixed.  Then, letting $R_m$ and $R_e$ be, repectively, the distances of the moon and Earth from the origin, and letting $\omega$ be the rate of rotation of the new frame relative to the original one (i.e. $2\pi$ per month):
$$
\begin{align}
L &= \frac{1}{2}\left(\dot{r}^2 + r^2(\dot\theta+\omega)^2 \right) + \frac{GM_e}{r_e} +\frac{GM_m}{r_m}\\
&= \frac{1}{2}\left(\dot{r}^2 + r^2\dot\theta^2 + 2r^2\dot\theta\omega\right) + \frac{GM_e}{r_e} +\frac{GM_m}{r_m} + \frac{1}{2}r^2\omega^2 \\
r_e &= \sqrt{r^2 + R_e^2 + 2R_e r cos\theta}\\
r_m &= \sqrt{r^2 + R_m^2 - 2R_m r cos\theta}\\
\end{align}
$$
Now that $\theta$ refers to an angle in the rotating reference frame, a "fictitious" additional potential of $-\frac{1}{2}r^2\omega^2$ appears in $V^{(\dagger)}$.  The image at the top shows this new potential $V$.  This picture makes it clear that in this new frame of reference we are surrounded by hilltops.  Actually there are 2 hilltops and 3 saddle points, but at all of these the spacecraft is weightless.  These 5 locations are called Lagrange points and if you manage to use up all your kinetic energy getting there you'll remain there indefinitely.

Let's see if we can work out what these Lagrange points are.  Well the Euler-Lagrange equation is $\frac{d}{dt}\frac{\partial L}{\partial \dot q} = \frac{\partial L}{\partial q}$ and our $q$s are $r$ and $\theta$, so
$$
\begin{align}
\frac{d}{dt}\frac{\partial L}{\partial \dot r} = \frac{\partial L}{\partial r}\tag{1}\\
\frac{d}{dt}\frac{\partial L}{\partial \dot \theta} = \frac{\partial L}{\partial \theta}\tag{2}\\
\end{align}
$$
If we calculate each side of this equation and then set all time derivatives to zero (because we're looking for solutions to Euler-Lagrange corresponding to no movement or coordinate acceleration) equations (1) and (2) become
$$
\begin{align}
0 &= - \frac{GM_e(r+R_e cos\theta)}{\left(r^2+R_e^2+2R_e r cos\theta\right)^{\frac{3}{2}}} - \frac{GM_m(r-R_m cos\theta)}{\left(r^2+R_m^2-2R_m r cos\theta\right)^{\frac{3}{2}}} + r\omega^2 \tag{3}\\
0 &= - \frac{GM_e R_e r sin\theta}{\left( r^2 + R_e^2 + 2 r R_e cos\theta \right)^{\frac{3}{2}}} + \frac{GM_m R_m r sin\theta}{\left( r^2 + R_m^2 - 2 r R_m cos\theta \right)^{\frac{3}{2}}}\tag{4}\\
\end{align}
$$
However, by the definition of $R_e$ and $R_m$ we know that $M_e R_e = M_m R_m$ so equation (4) simplifies to
$$
\begin{align}
r^2 + R_e^2 + 2 r R_e cos\theta &=  r^2 + R_m^2 - 2 r R_m cos\theta \\\implies r cos\theta &= \frac{1}{2}\frac{R_m^2-R_e^2}{R_e+R_m} \\
\implies r cos\theta&\approx \frac{1}{2}{R_m}\tag{4a} \\
\text{- OR -} \\
r sin \theta &= 0 \tag{4b}
\end{align}
$$
Now let's try to break this down into cases

Case I

Assume that equation (4b) holds because $\theta=0$.  Then equation (3) becomes
$$
\frac{GM_e(R_e+r)}{\left( r^2 + R_e^2 +2 r R_e \right)^{\frac{3}{2}}} = \frac{GM_m(R_m-r)}{\left( r^2 + R_m^2 -2 r R_m \right)^{\frac{3}{2}}} + r\omega^2
$$
Putting in $r=R_m$ and letting $d = R_m + R_e$ (i.e. the distance from the centre of Earth to the centre of the moon) we get
$$
\frac{GM_e}{d^2}= R_m \omega^2
$$
which is approximately correct since the LHS is the gravitational pull of the Earth on the moon divided by the moon's mass, and the RHS is the acceleration of the moon towards the origin.  (The slight inaccuracy is due to the fact that $\omega$ is relative to an origin that is not quite at the centre of the Earth.)

So we know that if $\theta = 0$ then $r \approx R_m$.  In fact it is obvious from the picture that there are two such solutions: $r = R_m - \epsilon_1$ and $r = R_m + \epsilon_2$.

Case II

Assume that equation (4b) holds because $\theta=\pi$.  Then equation (3) becomes
$$
\begin{align}
r\omega^2 &= \frac{GM_e(r-R_e)}{\left( r^2 + R_e^2 -2 r R_e \right)^{\frac{3}{2}}} + \frac{GM_m(r+R_m)}{\left( r^2 + R_m^2 +2 r R_m \right)^{\frac{3}{2}}}\\
&=\frac{GM_e}{(r-R_e)^2}+\frac{GM_m}{(r+R_m)^2}\\
&\approx  \frac{GM_e}{(r-R_e)^2}\\
&\approx  \frac{GM_e}{r^2}
\end{align}
$$
which, can be rearranged to give
$$
\begin{align}
r &\approx \sqrt[3]{\frac{GM_e}{\omega^2}}\\
&\approx R_m
\end{align}
$$
So we know that if $\theta = \pi$ then $r \approx R_m$ and it is clear from the picture that there is exactly one such solution.

Case III

Assume that equation (4a) holds and let's try to work out what $\theta$ could be.  Substituting (4a) into (3) gives:
$$
\begin{align}
r\omega^2 &= \frac{GM_e(r+R_e cos\theta)}{\left(r^2+R_e^2+R_e R_m\right)^{\frac{3}{2}}} + \frac{GM_m(r-R_m cos\theta)}{\left(r^2+R_m^2-R_m^2\right)^{\frac{3}{2}}}\\
&\approx \frac{GM_e r}{r^3} + \frac{GM_m(r-2 r cos^2\theta)}{r^3}\\
&\approx \frac{GM_e}{r^2}\\
\end{align}
$$
which, as in Case II, leads to $r \approx R_m$, combining this with equation (4a) gives $\theta = \pm \frac{\pi}{3}$

So, the final 2 of the 5 lagrange points are approximately $r = R_m$ and  $\theta = \pm \frac{\pi}{3}$

Upshot

So the instructions to our hamster astronaut on its mission to explore the solar system should be:

Stay in the Sun/Earth/Moon plane. Start off with the moon at 2 o'clock and keep the Earth behind you. Keep adding extra boosts to prevent stalling and falling back until the Moon is at 4 o'clock. Then stay there and await further instructions!  Good luck captain!

Alternatively, we could build a Lagrangian Guidance system with the familiar look and feel of being in a hamster tunnel.  The circles centre on the locations where $V$ is maximal at different distances, and the ones corresponding to all our current kinetic energy being consumed and the top of the hill are marked.  All the hamster captain needs to do is to keep aiming for the crosshairs and keep the All K.E. converted circle in front of the P.E. Max one.
This should help the captain counter the Coriolis force which would otherwise make the spacecraft tend to the right - a bit like a shopping trolley with a skew whiff wheel.

 Footnotes

  • ($\dagger$) And also a fictitious kinetic term $r^2\dot\theta\omega$ which depends on angular velocity $\dot \theta$.  This is responsible for the Coriolis Force.

POSTSCRIPT

The following short python script plots out the 3D graph above:


import matplotlib.pyplot as plt
import matplotlib.cm
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from numpy import pi, arctan2, sin, cos

# Constants
R = distance_to_moon = 384400e3
radius_of_earth = 6371e3
radius_of_moon = 1737e3
Me = mass_of_earth = 5.972e24
Mm = mass_of_moon = 7.34767309e22
omega = 2*pi / (27*24*3600)
G = 6.67408e-11

# Create 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Create square matrices x and y representing n**2 grid points d metres apart where 0,0 is centre of Earth
n = 30;
d = (distance_to_moon / n) * 1.3
x = (np.ones((2*n,2*n))*(np.arange(-n,n)+.5))
y = (np.ones((2*n,2*n))*(np.arange(-n,n)+.5)).transpose()
x *= d
y *= d

# Calculate distance from centre of Earth
r = (x**2 + y**2)**.5

# Calculate angle subtended at earth and distance from centre of moon
theta = arctan2(x,y)
r_m = ((r*sin(theta))**2 + (R - r*cos(theta))**2)**.5

# replace r and r_m near centres of earth and moon to avoid divide by zero
r = np.vectorize(lambda _r: max(_r, 40*radius_of_earth))(r)
r_m = np.vectorize(lambda _r: max(_r, 15*radius_of_moon))(r_m)

# Calculate potential energy in frame where moon and earth are stationary
V = -((r*omega)**2)/2. - G*Me/r - G*Mm/r_m

# Plot out V for disc
plt.ion()
ax = fig.gca(projection='3d')
ax.set_aspect('equal')       
ax.set_axis_off()
ax.plot_surface(x, y, V, rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm,  linewidth=0, shade=True)
#ax.plot_wireframe(x, y, V)
plt.show()

Comments