Calculate the period for tidal locking of the Earth and Moon

Calculate Moon synchronous period.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the synchronous period between Earth's rotation and the Moon's orbit\n",
    "\n",
    "The Moon is currently gradually moving away from the Earth, at a rate of about 3.8 cm per year. At the same time (and for the same reason), the Earth's rotation is gradually slowing down.  This is due to exchange of angular momentum between the Earth's rotation and the Moon's orbit, due to the Moon's pull on the Earth's tidal bulge, as explained [here](http://www.astronomy.ohio-state.edu/%7Epogge/Ast161/Unit4/tides.html).  \n",
    "\n",
    "While the process is too slow to completely bring the Earth's rotation and the Moon's orbit to the same period (at which they would be tidally locked) before the Sun becomes a red giant, likely swallowing the Earth and certainly evaporating Earth's oceans, it is still an interesting physics problem to work out what period the Moon would have if there were time for the process to go to completion. \n",
    "\n",
    "This notebook calculates what the synchronous rotation / orbital period would be for the Earth-Moon system, assuming the current angular momentum in the system is conserved. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.constants as c\n",
    "import astropy.units as u\n",
    "from math import pi\n",
    "from sympy import symbols, Eq, solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the masses, radii, and periods we need: \n",
    "\n",
    "# Moon constants from\n",
    "# https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html\n",
    "R_Moon = 1738.1 * u.km\n",
    "M_Moon = 0.07346e24 * u.kg \n",
    "R_Earth_Moon = 0.3844e6 * u.km\n",
    "P_Moon = 27.3217 * u.day\n",
    "M_Earth = c.M_earth\n",
    "R_Earth = c.R_earth\n",
    "P_Earth = u.day"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Angular momentum for each rotating or orbiting object is $L = I \\omega$, where $I$ is the moment of inertia, and $\\omega$ is the angular frequency, related to the rotational or orbital period $P$ by $\\omega = \\frac{2\\pi}{P}$. \n",
    "\n",
    "A uniform density sphere has moment of inertia $I = 0.4 MR^2$, but the Moon and (especially) the Earth are somewhat centrally condensed, and thus have slightly smaller moments of inertia.  Values below are from the Wikipedia entry for [\"moment of inertia factor\"](https://en.wikipedia.org/wiki/Moment_of_inertia_factor). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_Earth = 0.33 * M_Earth * R_Earth**2\n",
    "L_Earth = (2 * pi * I_Earth)/P_Earth\n",
    "I_Moon = 0.39 * M_Moon * R_Moon**2\n",
    "L_Moon_rot = (2 * pi * I_Moon)/P_Moon\n",
    "I_Moon_orbit = M_Moon * R_Earth_Moon**2\n",
    "L_Moon_orbit = (2 * pi * I_Moon_orbit)/P_Moon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Total angular momentum now\n",
    "L_1 = L_Earth + L_Moon_orbit + L_Moon_rot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.8320805 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 0.8320805>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Most of the angular momentum is in the Moon's orbit: \n",
    "(L_Moon_orbit/L_1).decompose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After setting the angular momentum now ($L_1$ above) to the eventual angular momentum (where all three bodies have the same period), and using Kepler's third law to re-write the new (unknown) semimajor axis of the Moon's orbit in terms of its orbital period, we end up with an equation to solve for the new period that involves the terms below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partial terms in eventual angular momentum equation to solve:\n",
    "A = 2 * pi * (I_Earth + I_Moon)/L_1\n",
    "B = 2 * pi * I_Moon_orbit / (L_1 * P_Moon**(4./3.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check out the values of these terms, and make sure we convert them from quantities into scalar numbers, both in units based on days: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.16809413 \\; \\mathrm{d}$"
      ],
      "text/plain": [
       "<Quantity 0.16809413 d>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.27626727 \\; \\mathrm{\\frac{1}{d^{1/3}}}$"
      ],
      "text/plain": [
       "<Quantity 0.27626727 1 / d(1/3)>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Want B also in day-like units, make the unit up: \n",
    "inverse_day_cuberoot = (u.day)**(-1/3)\n",
    "B.to(inverse_day_cuberoot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save scalars of these in day-based units that we can use in our equation: \n",
    "a = A.to(u.d).value\n",
    "b = B.to(inverse_day_cuberoot).value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define our symbol to solve for, and the equation: \n",
    "P_new = symbols('P_new')\n",
    "eq = Eq(P_new,  a + b*P_new**(4/3))\n",
    "sol = solve(eq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.200518234104636, 46.9175571571954]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the solution - there is a longer-period solution, \n",
    "# but also one at a very short period if the Moon were to \n",
    "# spiral *in*!   Both of these values are in days, since that's\n",
    "# how we specified a and b. \n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve for new semimajor axis of Moon: \n",
    "P_new = sol[1] * u.day\n",
    "# This is Kepler's third law, solved for P: \n",
    "R_Earth_Moon_new = ( P_new**2 * (M_Moon + M_Earth) * c.G /(4*pi**2))**(1/3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1.4353082 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 1.43530823>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# New orbit is farther away as expected: \n",
    "(R_Earth_Moon_new/R_Earth_Moon).decompose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1.0018068 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 1.00180679>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check angular momentum; should be the same as before: \n",
    "L_2 = (2 * pi/P_new) * (I_Earth + I_Moon + M_Moon * R_Earth_Moon_new**2)\n",
    "\n",
    "L_2/L_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the bottom line is, the tidal locking would be achieved when the Earth's day, the month (orbital period of the Moon), and the Moon's rotational period all reached a value of 46.9 days, at which point the Moon would orbit 43.5% farther from Earth than it does today. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}