{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 06: Environmental Perturbations\n",
    "\n",
    "This tutorial demonstrates how to model and simulate non-gravitational and third-body environmental perturbations in `OpenGNC`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--- \n",
    "## Theory Prerequisites\n",
    "\n",
    "### 1. Atmospheric Drag\n",
    "The drag acceleration opposes the spacecraft's velocity relative to the rotating atmosphere:\n",
    "$$\n",
    "\\mathbf{a}_{drag} = -\\frac{1}{2} \\rho v_{rel}^2 \\frac{C_D A}{m} \\hat{\\mathbf{v}}_{rel}\n",
    "$$\n",
    "where $\\rho$ is local density, $C_D$ is drag coefficient, $A$ is area, and $m$ is mass.\n",
    "\n",
    "### 2. Solar Radiation Pressure (SRP)\n",
    "The force exerted by solar photons. It pushes away from the Sun:\n",
    "$$\n",
    "\\mathbf{a}_{srp} = - \\nu P_{sun} C_R \\frac{A}{m} \\hat{\\mathbf{u}}_{sun}\n",
    "$$\n",
    "where $\\nu$ is the shadow function (1 in sunlight, 0 in eclipse), $P_{sun}$ is radiation pressure at distance, and $C_R$ is reflectivity coefficient.\n",
    "\n",
    "### 3. Third-Body Gravity\n",
    "Gravitational pull from other celestial bodies (Sun, Moon):\n",
    "$$\n",
    "\\mathbf{a}_{3B} = \\sum \\mu_k \\left( \\frac{\\mathbf{s}_k - \\mathbf{r}}{|\\mathbf{s}_k - \\mathbf{r}|^3} - \\frac{\\mathbf{s}_k}{|\\mathbf{s}_k|^3} \\right)\n",
    "$$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Environment modules loaded.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from opengnc.disturbances.drag import LumpedDrag\n",
    "from opengnc.environment.density import HarrisPriester\n",
    "from opengnc.disturbances.srp import Canonball\n",
    "from opengnc.disturbances.gravity import ThirdBodyGravity, TwoBodyGravity\n",
    "from opengnc.propagators.cowell import CowellPropagator\n",
    "\n",
    "print(\"Environment modules loaded.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Atmospheric Drag & Density\n",
    "\n",
    "Let's evaluate the atmospheric density and resulting drag acceleration at different altitudes and compare them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drag at 400km: Accordance Vector: [-0.00000000e+00 -2.74518523e-06 -0.00000000e+00] m/s^2\n",
      "Magnitude: 2.745185e-06 m/s^2\n"
     ]
    }
   ],
   "source": [
    "density_model = HarrisPriester()\n",
    "drag_model = LumpedDrag(density_model)\n",
    "\n",
    "# Spacecraft Properties\n",
    "mass = 4.0 # 4 kg (3U CubeSat)\n",
    "area = 0.03 # 0.03 m^2 cross section\n",
    "cd = 2.2 # Drag coefficient\n",
    "\n",
    "# State: 400km Altitude\n",
    "r_400 = np.array([6778137.0, 0.0, 0.0])\n",
    "v_400 = np.array([0.0, 7670.0, 0.0])\n",
    "jd = 2451545.0\n",
    "\n",
    "a_drag = drag_model.get_acceleration(r_400, v_400, jd, mass, area, cd)\n",
    "print(f\"Drag at 400km: Accordance Vector: {a_drag} m/s^2\")\n",
    "print(f\"Magnitude: {np.linalg.norm(a_drag):.6e} m/s^2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Solar Radiation Pressure (SRP)\n",
    "\n",
    "Evaluating are drag acceleration vs SRP at equivalent conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SRP Accel Vector: [-8.00476320e-09  4.01254285e-08  1.73962601e-08] m/s^2\n",
      "Magnitude: 4.446073e-08 m/s^2\n",
      "\n",
      "Ratio (Drag/SRP): 61.74x\n",
      "At 400km LEO, Atmospheric Drag dominates heavily. SRP becomes relevant at higher altitudes (MEO/GEO).\n"
     ]
    }
   ],
   "source": [
    "srp_model = Canonball()\n",
    "cr = 1.3 # Reflectivity\n",
    "\n",
    "a_srp = srp_model.get_acceleration(r_400, jd, mass, area, cr)\n",
    "print(f\"SRP Accel Vector: {a_srp} m/s^2\")\n",
    "print(f\"Magnitude: {np.linalg.norm(a_srp):.6e} m/s^2\")\n",
    "\n",
    "print(f\"\\nRatio (Drag/SRP): {np.linalg.norm(a_drag)/np.linalg.norm(a_srp):.2f}x\")\n",
    "print(\"At 400km LEO, Atmospheric Drag dominates heavily. SRP becomes relevant at higher altitudes (MEO/GEO).\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Third-Body Gravity\n",
    "\n",
    "Lastly, let's see the acceleration from the Moon and Sun."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Third-Body Accel Vector: [5.58127159e-08 5.88545936e-07 1.51971346e-07] m/s^2\n",
      "Magnitude: 6.104070e-07 m/s^2\n"
     ]
    }
   ],
   "source": [
    "third_body = ThirdBodyGravity()\n",
    "a_3b = third_body.get_acceleration(r_400, jd)\n",
    "print(f\"Third-Body Accel Vector: {a_3b} m/s^2\")\n",
    "print(f\"Magnitude: {np.linalg.norm(a_3b):.6e} m/s^2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Total Disturbance Propagation\n",
    "\n",
    "We can aggregate these and feed into a propagator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Propagated position after 1 hour: [-4.06313714e+06 -5.43038765e+06  1.11069242e-01] m\n"
     ]
    }
   ],
   "source": [
    "propagator = CowellPropagator()\n",
    "\n",
    "# Combined Perturbations Function\n",
    "def environmental_perturbations(t, r, v):\n",
    "    # Note: t is local inside integrator loop relative to integration start\n",
    "    curr_jd = jd + (t / 86400.0) \n",
    "    \n",
    "    acc_drag = drag_model.get_acceleration(r, v, curr_jd, mass, area, cd)\n",
    "    acc_srp  = srp_model.get_acceleration(r, curr_jd, mass, area, cr)\n",
    "    acc_3b   = third_body.get_acceleration(r, curr_jd)\n",
    "    \n",
    "    return acc_drag + acc_srp + acc_3b\n",
    "\n",
    "# Propagate short window\n",
    "r_f, v_f = propagator.propagate(r_400, v_400, 3600, perturbation_acc_fn=environmental_perturbations, dt_step=60)\n",
    "print(f\"\\nPropagated position after 1 hour: {r_f} m\")"
   ]
  }
 ],
 "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
}