05_ae470_keplerian_orbits_tof.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d972813-d41d-4d6f-87fa-578a4ec5adfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e0be755-e2c3-423b-aed3-7e1f53e3a545",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import modsim and install if necessary\n",
    "try:\n",
    "    from modsim import *\n",
    "except ImportError:\n",
    "    print(\"Download modsim from the course repository at https://github.com/jeffwalton/AE470\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3f7b7a4-6dbc-41b8-a3c0-04ad84d24c21",
   "metadata": {},
   "source": [
    "# Keplerian Orbits (continued)\n",
    "\n",
    "From the last notebook, the period of the orbit, $ T $, is\n",
    "\n",
    "\\begin{equation}\n",
    "T = 2 \\pi \\sqrt{ \\frac{a^3}{\\mu} }\n",
    "\\end{equation}\n",
    "\n",
    "The *mean orbital motion*, $ n $ is defined as \n",
    "\n",
    "\\begin{equation}\n",
    "n \\triangleq \\frac{ 2 \\pi }{ T } = \\sqrt{ \\frac{\\mu}{a^3} }\n",
    "\\end{equation}\n",
    "\n",
    "This is the average angular rate of the orbit.\n",
    "\n",
    "## Time of Flight\n",
    "\n",
    "To solve the orbit problem, the position of $ m_2 $ with repsect to time must be known, specifically $ \\theta(t) $, true anomaly $ \\theta $ as a function of time, $ t $.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e9c2678-7778-4e78-90f3-f77c8f0ce5f8",
   "metadata": {},
   "source": [
    "### Circular Orbits\n",
    "\n",
    "Since $ h = r^2 \\dot{ \\theta } $, rearranging \n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{ \\theta } = \\frac{h}{r^2}\n",
    "\\end{equation}\n",
    "\n",
    "For a circular orbit, $ r = a $ and $ h = \\sqrt{ \\mu  a } $.\n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{ \\theta } = \\sqrt{ \\frac{\\mu}{r^3} } = \\sqrt{ \\frac{\\mu}{a^3} } = n\n",
    "\\end{equation}\n",
    "\n",
    "**For a circular orbit, the angular rate is the mean orbital motion.**  \n",
    "\n",
    "Note: $ \\dot{ \\theta } = \\frac{h}{r^2} $ is only constant if $ r $ is constant (a circular orbit).\n",
    "\n",
    "Integrating this gives the relationship between $ \\theta $ and $ t $.\n",
    "\n",
    "\\begin{equation}\n",
    "\\theta(t) - \\theta(t_0) = n(t - t_0)\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "027bdfe7-f35e-407a-97f6-f11a22d8e69e",
   "metadata": {},
   "source": [
    "### Elliptical Orbits\n",
    "\n",
    "#### Eccentric Anomaly\n",
    "\n",
    "A direct relationship between time, $ t $, and true anamoly, $ \\theta $, is difficult to obtain.\n",
    "\n",
    "An indirect realtionship with time is created by defining an auxiliary variable called *eccentric anamoly*, $ E $.  The geometric interpertation of $ E $ is shown in the following diagram. \n",
    "\n",
    "![eccentric anomaly](https://upload.wikimedia.org/wikipedia/commons/thumb/4/48/EccentricAnomaly.svg/437px-EccentricAnomaly.svg.png)\n",
    "\n",
    "Eccentric anomaly is the angle between the apse line and a fictitious replica of the orbiting body that moves on a circle at a vertical projection from the major axis. $ E $ is measured at the center of the circle.\n",
    "\n",
    "The relationship between true anomaly, $ \\theta $, and eccentric amomaly, $ E $, is\n",
    "\n",
    "\\begin{equation}\n",
    "\\boxed{ \\tan{\\frac{\\theta}{2}} = sqrt{ \\frac{1+e}{1-e} } \\tan{\\frac{E}{2}} }\n",
    "\\end{equation}\n",
    "\n",
    "A lengthy derivation can be found in most texts, see deRuiter, pages 84-86, or [Eccentric Anomaly at wikipedia.org](https://en.wikipedia.org/wiki/Eccentric_anomaly). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d806fab8-bd10-40fa-890d-8c2e5b139e36",
   "metadata": {},
   "source": [
    "#### Mean Anomaly\n",
    "\n",
    "A second variable, $ M $ is defined as\n",
    "\n",
    "\\begin{equation}\n",
    "M = n(t - t_0)\n",
    "\\end{equation}\n",
    "\n",
    "![mean anomaly](https://upload.wikimedia.org/wikipedia/commons/thumb/b/bc/Mean_Anomaly.svg/522px-Mean_Anomaly.svg.png)\n",
    "\n",
    "Mean anomaly is the angle between the apse line and a fictitious replica of the orbiting body that moves uniformly on a circle with the same period as the elliptical orbit (measured at the center of the circle).\n",
    "\n",
    "A lengthy derivation can be found in most texts, see deRuiter, pages 86-87, or [Mean Anomaly at wikipedia.org](https://en.wikipedia.org/wiki/Mean_anomaly)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9324fa23-ffe3-45f8-89f7-0d2d64026701",
   "metadata": {},
   "source": [
    "#### Kelper's Equation\n",
    "\n",
    "*Kepler's equation* defines the relationship between mean anomaly, $ M $, and eccentric anomaly, $ E $.\n",
    "\n",
    "\\begin{equation}\n",
    "\\boxed{  M = E - e \\sin{ E }  }\n",
    "\\end{equation}\n",
    "\n",
    "A lengthy derivation can be found in most texts (see deRuiter, pages 84-88 or [Kepler's Equation at wikipedia.org](https://en.wikipedia.org/wiki/Kepler%27s_equation)).\n",
    "\n",
    "We cannot solve this transcendental equation directly for $ E $. A solution requires an iterative, numerical procedure such as Newton's method.\n",
    "\n",
    "\\begin{equation}\n",
    "E_{k+1} = E_{k} - \\frac{f(E_{k})}{f'(E_{k})}\n",
    "\\end{equation}\n",
    "\n",
    "where\n",
    "\\begin{align}\n",
    "f(E) &= E - e \\sin{E} - M \\\\\n",
    "f'(E) &= 1 - e \\cos{E} \\\\\n",
    "\\end{align}\n",
    " \n",
    "Iterate until the desired tolerance, $ tol $, is reached, e.g. $ 10^{-8} $.\n",
    "\n",
    "\\begin{equation}\n",
    "| E_{k} - E_{k+1} | = \\frac{f(E_{k})}{f'(E_{k})} \\leq tol\n",
    "\\end{equation}\n",
    "\n",
    "An initial guess $ E_0 = M $ is usually a good starting point.\n",
    "\n",
    "There are versions of Kepler's equation for parabolic and hyperbolic orbits."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c333127-d657-4ff5-9ca7-00000278cee4",
   "metadata": {},
   "source": [
    "---\n",
    "+ For an elliptical orbit, to find true anomaly, $ \\theta $, given time since periapsis passage, $ t - t_0 $, \n",
    "\n",
    " 1. Compute mean anomaly, $ M $, from\n",
    "\n",
    "\\begin{equation}\n",
    "M = n(t - t_0)\n",
    "\\end{equation}\n",
    "\n",
    " 2. Solve Kepler's equation numerically for the eccentric anomaly, $ E $\n",
    "\n",
    "\\begin{equation}\n",
    "M = E - e \\sin{ E }\n",
    "\\end{equation}\n",
    "\n",
    " 3. Compute the true anomaly, $ \\theta $ from\n",
    "\n",
    "\\begin{equation}\n",
    "\\tan{\\frac{\\theta}{2}} = sqrt{ \\frac{1+e}{1-e} } \\tan{\\frac{E}{2}}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3b308bd-aa86-482e-901a-c49337580846",
   "metadata": {},
   "source": [
    "---\n",
    "+ The procedure can be reversed. Given true anomaly, $ \\theta $, the time since periapsis passage, $ t - t_0 $, can be determined by \n",
    "\n",
    " 1. Compute the eccentric anomaly, $ E $, from\n",
    "\n",
    "\\begin{equation}\n",
    "\\tan{\\frac{\\theta}{2}} = sqrt{ \\frac{1+e}{1-e} } \\tan{\\frac{E}{2}}\n",
    "\\end{equation}\n",
    "\n",
    " 2. Compute mean anomaly, $ M $, from\n",
    "\n",
    "\\begin{equation}\n",
    "M = E - e \\sin{ E }\n",
    "\\end{equation}\n",
    "\n",
    " 3. Compute the true anomaly, $ \\theta $ from\n",
    "\n",
    "\\begin{equation}\n",
    "M = n(t - t_0)\n",
    "\\end{equation}\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2f92a51-914a-4d52-b9c2-48c31a31a943",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b018415c-19fc-4262-b85b-cd601e182ed5",
   "metadata": {},
   "source": [
    "### 1. Bazzocchi 2023\n",
    "\n",
    "If a spacecraft is in orbit with an eccentricity $ e = 0.5 $, period $ T = 7.82 h $ and a true anomaly $ \\theta = \\pi/2 $, how long has it been since periapsis passage (assume $ 𝑡_0 = 0 $)?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d64fe2d-a4a8-4829-8bfa-dd6bd0281f8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ecc = 0.5\n",
    "T = 7.82 * 60 * 60 # s\n",
    "theta = np.pi / 2\n",
    "\n",
    "E = 2 * np.arctan( np.sqrt((1-ecc)/(1+ecc))*np.tan(theta/2) )\n",
    "M = E - ecc*np.sin(E)\n",
    "n = 2*np.pi/T\n",
    "\n",
    "t = M/n\n",
    "\n",
    "print(f\"time after perigee passage = {t:.1f} seconds, or {t/60:.1f} minutes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bd2022d-8712-40f7-a3bc-8f9744632fb1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7459e55f-2318-484d-bcbe-e1017d7dd5eb",
   "metadata": {},
   "source": [
    "### 2. Curtis, Example 3.1, page 149\n",
    "\n",
    "A geocentric elliptical orbit has a perigee radius of 9600 km and an apogee radius of 21,000 km.  Calculate the time to fly from perigee to true anomaly of $ 120^{\\circ} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "661c2d91-9585-45e4-b1c9-757ec9cfe5ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3a5b69c-5e32-405d-9311-c427fb0211a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_a = 21000 # km\n",
    "r_p = 9600 # km\n",
    "theta = 120 / 180 * np.pi\n",
    "mu = 3.986000e+5 # km^3/s^2\n",
    "\n",
    "ecc = ( r_a - r_p )/( r_a + r_p ) \n",
    "h = np.sqrt( r_p*( mu * ( 1 + ecc*np.cos(0) ) ) )\n",
    "# a = (h**2/mu)/(1-ecc**2)\n",
    "a = (r_a + r_p)/2\n",
    "T = (2*np.pi/np.sqrt(mu))*np.sqrt(a**3)\n",
    "E = 2 * np.arctan( np.sqrt((1-ecc)/(1+ecc))*np.tan(theta/2) )\n",
    "M = E - ecc*np.sin(E)\n",
    "n = 2*np.pi/T\n",
    "# n = np.sqrt(mu/a**3)\n",
    "\n",
    "t = M/n\n",
    "\n",
    "print(f\"time after perigee passage = {t:.1f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d2f1a7f-1b89-4a50-9451-8f80f587563c",
   "metadata": {},
   "source": [
    "### 3. Curtis, Example 3.2, page 150\n",
    "\n",
    "In the previous example, find the true anomaly at 3 h after perigee passage.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba16e69a-5fc3-4237-83e9-644a55734e06",
   "metadata": {},
   "source": [
    "First, let's write a python function to solve Kelper's equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84b36a8b-6018-4a14-acae-9a367885bb98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solve Kepler's equation\n",
    "def keplers_equation(ecc, M, err_tol=1.0e-6, max_iter=10):\n",
    "    if M < np.pi:\n",
    "        E = M + ecc/2\n",
    "    else:\n",
    "        E = M - ecc/2\n",
    "    ratio = 0.1\n",
    "    # print(f\"init: {ecc=}, {M=}, {err_tol=}, {max_iter=}, {E=}, {ratio=}\")\n",
    "    for iter in range(0,max_iter+1):\n",
    "        # print(f\"{iter=}, {E=}\" )\n",
    "        f = E - ecc*np.sin(E) - M\n",
    "        f_prime = 1 - ecc*np.cos(E)\n",
    "        ratio = f/f_prime\n",
    "        E -= ratio\n",
    "        # print(f\"{f=}, {f_prime=}, {ratio=}, {E=}\")\n",
    "        if np.abs(ratio) < err_tol:\n",
    "            break\n",
    "    return E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2991ec06-1763-4ad0-af8a-4c25b7001503",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 10800 # s\n",
    "print(f\"eccentricity, ecc (from above): {ecc:.5f}\")\n",
    "print(f\"Orbit period, T (from above): {T:.1f} seconds\")\n",
    "M = 2*np.pi*t/T\n",
    "print(f\"Mean Anomaly, M: {M:.4f} radians\")\n",
    "\n",
    "E = keplers_equation(ecc, M, err_tol=1e-6)\n",
    "\n",
    "theta = 2*np.pi + 2 * np.arctan( np.sqrt((1+ecc)/(1-ecc))*np.tan(E/2) )\n",
    "theta_deg = theta * 180/np.pi\n",
    "\n",
    "print(f\"true anomaly, theta = {theta_deg:.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95e42025-5353-46aa-b2b3-7cd567e473d0",
   "metadata": {},
   "source": [
    "### 4. Curtis, Example 3.3, page 151\n",
    "\n",
    "A satellite is in a 500 km by 500 km orbit with its apse line parallel to the line from the earth to the sun.  Find the time that the satellite is in the earth's shadow if\n",
    "  1. the apogee is toward the sun\n",
    "  2. the perigee is toward the sun.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d73d3478-949c-45d1-aa42-9184a19e69dc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cab220a1-8989-41da-af42-c27fe821157e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0e9cd065-9036-48bc-a0f9-f2362555d10e",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "    R.R. Bate, D.D. Mueller, J.E. White, W.W. Saylor; Fundamentals of Astrodynamics, Dover, 2nd Ed, 2020.  \n",
    "\n",
    "    H.D. Curtis; Orbital Mechanics for Engineering Students; Butterworth-Heinemann, 4th Ed, 2021.  \n",
    "\n",
    "    A.H.J. de Ruiter, C.J. Damaren, J.R. Forbes; Spacecraft Dynamics and Control: An Introduction; Wiley, 2013.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aff3fc6-3d5a-40c3-8214-660ce6111bf3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}