{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 05: Comparison: Cowell vs Encke\n",
    "\n",
    "This tutorial compares absolute numerical integration (**Cowell's Method**) and perturbative numerical integration (**Encke's Method**) within `OpenGNC`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--- \n",
    "## Theory Prerequisites\n",
    "\n",
    "### 1. Cowell's Method\n",
    "Integrates the absolute equations of motion directly. The state integrated is the absolute position $\\mathbf{r}$ and velocity $\\mathbf{v}$:\n",
    "$$\n",
    "\\ddot{\\mathbf{r}} = -\\frac{\\mu}{r^3}\\mathbf{r} + \\mathbf{a}_{pert}\n",
    "$$\n",
    "\n",
    "### 2. Encke's Method\n",
    "Useful for orbits where perturbations are small. It integrates the **deviation** $\\delta \\mathbf{r}$ from a reference Keplerian orbit $\\mathbf{r}_{ref}$:\n",
    "$$\n",
    "\\mathbf{r} = \\mathbf{r}_{ref} + \\delta \\mathbf{r}\n",
    "$$\n",
    "$$\n",
    "\\ddot{\\delta \\mathbf{r}} = \\frac{\\mu}{r_{ref}^3} \\left( f(q)\\mathbf{r}_{ref} - \\delta\\mathbf{r} \\right) + \\mathbf{a}_{pert}\n",
    "$$\n",
    "where $q = \\frac{\\delta\\mathbf{r} \\cdot (\\delta\\mathbf{r} - 2\\mathbf{r}_{ref})}{r_{ref}^2}$. Encke's method reduces numerical roundoff errors for bound orbits with small disturbances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Propagator modules loaded.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "from opengnc.propagators.cowell import CowellPropagator\n",
    "from opengnc.propagators.encke import EnckePropagator\n",
    "from opengnc.disturbances.gravity import J2Gravity, TwoBodyGravity\n",
    "from opengnc.integrators.rk4 import RK4\n",
    "\n",
    "print(\"Propagator modules loaded.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Setup Simulation Environment\n",
    "\n",
    "We define a LEO orbit and introduce a $J_2$ disturbance to trigger the perturbative integrator loops on both solvers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial Orbit (m, m/s)\n",
    "r_0 = np.array([7000000.0, 0.0, 0.0])\n",
    "v_0 = np.array([0.0, 7546.0 * np.cos(np.radians(28)), 7546.0 * np.sin(np.radians(28))]) \n",
    "\n",
    "duration = 10000.0 # [s] - Approx 2 orbits\n",
    "dt_step = 10.0 # Steps size [s]\n",
    "\n",
    "def j2_perturbation(t, r, v):\n",
    "    return J2Gravity().get_acceleration(r) - TwoBodyGravity().get_acceleration(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Cowell Propagation Benchmarking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cowell Final Pos:  [-1310807.66896629 -6060113.86490409 -3226579.61689067] m\n",
      "Cowell Execution Time: 0.096331 seconds\n"
     ]
    }
   ],
   "source": [
    "cowell = CowellPropagator()\n",
    "\n",
    "t0 = time.perf_counter()\n",
    "r_f_cowell, v_f_cowell = cowell.propagate(r_0, v_0, duration, perturbation_acc_fn=j2_perturbation, dt_step=dt_step)\n",
    "t_cowell = time.perf_counter() - t0\n",
    "\n",
    "print(f\"Cowell Final Pos:  {r_f_cowell} m\")\n",
    "print(f\"Cowell Execution Time: {t_cowell:.6f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Encke Propagation Benchmarking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Encke Final Pos:  [-1311079.24910403 -6060081.37923322 -3226564.18300095] m\n",
      "Encke Execution Time: 0.337163 seconds\n"
     ]
    }
   ],
   "source": [
    "encke = EnckePropagator(rect_tol=1e-6) # Triggers rectification when deviation gets large\n",
    "\n",
    "t0 = time.perf_counter()\n",
    "r_f_encke, v_f_encke = encke.propagate(r_0, v_0, duration, perturbation_acc_fn=j2_perturbation, dt_step=dt_step)\n",
    "t_encke = time.perf_counter() - t0\n",
    "\n",
    "print(f\"Encke Final Pos:  {r_f_encke} m\")\n",
    "print(f\"Encke Execution Time: {t_encke:.6f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Comparison Summary\n",
    "\n",
    "Let's compare the magnitude difference between both methods and speed ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final Position Difference: 273.951264 m\n",
      "Speed Advantage (Cowell/Encke ratio): 0.29x\n",
      "\n",
      "Note: Encke is often faster/more accurate over long continuous orbits because it integrates smaller variables, reducing numerical noise limits on large loops.\n"
     ]
    }
   ],
   "source": [
    "pos_diff = np.linalg.norm(r_f_cowell - r_f_encke)\n",
    "speedup = t_cowell / t_encke if t_encke > 0 else 1.0\n",
    "\n",
    "print(f\"Final Position Difference: {pos_diff:.6f} m\")\n",
    "print(f\"Speed Advantage (Cowell/Encke ratio): {speedup:.2f}x\")\n",
    "print(\"\\nNote: Encke is often faster/more accurate over long continuous orbits because it integrates smaller variables, reducing numerical noise limits on large loops.\")"
   ]
  }
 ],
 "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}