solutions/assignment_02_solution_ae470.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "50aea7eb-817d-4619-b695-da0059bf79c6",
   "metadata": {},
   "source": [
    "## AE 470 - Assignment 2\n",
    "## Due in 1 week (Tuesday, Febuary 13, 2024 at 12:00 noon)\n",
    "\n",
    "Submit your solutions in the this Jupyter notebook.  Email the notebook file to me at [jwalton@clarkson.edu](jwalton@clarkson.edu).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "833c49f2-fb4b-4ebb-a22d-94ce9b2928bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Put your name here -- SOLUTION\n",
    "\n",
    "# Add your collaborators names here\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a3d26faa-b487-4269-a0f4-8cdb4e489a58",
   "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": 2,
   "id": "c0355495-256c-4513-abbc-e7348426b21f",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    from modsim import *\n",
    "except ImportError:\n",
    "    print(\"Download modsim from the course repository at https://github.com/jeffwalton/AE470\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bb51d8f3-1e41-4041-9ca4-0cc8719fc6eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6343bc87",
   "metadata": {},
   "source": [
    "## 1. Curtis, page 78, Example 2.5\n",
    "\n",
    "Calculate the altitude $ z_{GEO} $ and the speed $ v_{GEO} $of a geostationary earth satellite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a4d37e49",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6378.137"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define Constants \n",
    "mu_Earth = 3.98600e+5 # km^3/s^2\n",
    "r_Earth  = 6378137e-3 # km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4facd4b3-4318-4ef9-9a4f-38d08d852bc2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z_Geo = 35785.99 km \n",
      "v_Geo = 3.075 km/s\n"
     ]
    }
   ],
   "source": [
    "T_GeoStationary = 23*60*60+56*60+4 # [Seconds]\n",
    "r_GeoStationary = ((T_GeoStationary*np.sqrt(mu_Earth))/(2*np.pi))**(2/3) # [km] r is magnitude of positon vector (Total height)\n",
    "z_Geo = r_GeoStationary-r_Earth # Altitude Abover Earth [km]\n",
    "v_Geo = (2*np.pi*r_GeoStationary)/T_GeoStationary # [km/s]\n",
    "print(f'z_Geo = {z_Geo:.2f} km \\nv_Geo = {v_Geo:.3f} km/s')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "838a5a36",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d4d17522",
   "metadata": {},
   "source": [
    "## 2. Curtis, page 86, Example 2.7\n",
    "\n",
    "An earth satellite is in an orbit with a perigee altitude $ z_{p} = 400 \\, \\mathrm{km} $ and an apogee altitude $ z_{a} = 4000 \\, \\mathrm{km} $.\n",
    "\n",
    "Find each of the following quantities:\n",
    "  1. eccentricity, $ e $\n",
    "  2. angular momentum, $ h $\n",
    "  3. perigee velocity, $ v_{p} $\n",
    "  4. apogee velocity, $ v_{a} $\n",
    "  5. semi-major axis, $ a $\n",
    "  6. period of the orbit, $ T $\n",
    "  7. true anomaly when $ r = \\bar{r}_{\\theta} = \\sqrt{r_{p} r_{a}} $\n",
    "  8. satellite speed when $ r = \\bar{r}_{\\theta} = \\sqrt{r_{p} r_{a}} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d47afe6f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The eccentricity 0.2098.\n",
      "The Angular Momentum 57172.0\n",
      "Apogee Velocity 5.509 km/s\n",
      "Perogee Velocity 8.435 km/s\n",
      "Semi-major axis 8578.0 km\n",
      "Period Time 2.196 hrs\n"
     ]
    }
   ],
   "source": [
    "mu=3.986000e+5 # earth, km^3/s^2\n",
    "r_earth= 6378\n",
    "theta = 0\n",
    "Ra = 4000+r_earth\n",
    "Rp = 400+r_earth\n",
    "e = (Ra-Rp) / (Ra+Rp)\n",
    "print(f\"The eccentricity {e:.4f}.\")\n",
    "\n",
    "a = (Ra+Rp)/2\n",
    "va = np.sqrt(mu * (2 / Ra - 1 / a ))\n",
    "vp = np.sqrt(mu * (2 / Rp - 1 / a ))\n",
    "r0 = Vector( Ra, 0, 0 ) # km\n",
    "v0 = Vector( 0, va, 0 ) # km/s\n",
    "\n",
    "h = vector_mag(np.cross(r0,v0))\n",
    "print(f\"The Angular Momentum {h:.1f}\")\n",
    "print(f\"Apogee Velocity {va:.3f} km/s\")\n",
    "print(f\"Perogee Velocity {vp:.3f} km/s\")\n",
    "\n",
    "\n",
    "p = h**2 / mu\n",
    "r = p / (1+e* np.cos(theta))\n",
    "\n",
    "ab = p / (1-e**2)\n",
    "T = 2*np.pi*np.sqrt(a**3/mu)\n",
    "print(f\"Semi-major axis {a:.1f} km\")\n",
    "print(f\"Period Time {T/3600:.3f} hrs\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c9049e98",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "96.090585043348\n"
     ]
    }
   ],
   "source": [
    "#7 True Anomaly\n",
    "r = np.sqrt(Rp*Ra)\n",
    "theta_rad = np.arccos((p-r)/(r*e))\n",
    "theta = theta_rad*180/np.pi\n",
    "print(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f9e12ac2-fc60-4c86-8ae7-b463938b5e23",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.970218437447556\n"
     ]
    }
   ],
   "source": [
    "#8 Satellite Speed\n",
    "v = np.sqrt(mu*((2/r)-(1/a)))\n",
    "print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0227c764",
   "metadata": {},
   "source": [
    "## 3. Curtis, page 89, Example 2.8\n",
    "\n",
    "At two points on a geocentric orbit, altitud and true anamoly are $ z_1 = 1545 \\, \\mathrm{km} $, $ \\theta_1 = 126^{\\circ} $ and $ z_2 = 852 \\, \\mathrm{km} $, $ \\theta_2 = 58^{\\circ} $, respectively. Find the (a) eccentricity, (b) the altitude of perigee, (c) the semi-major axis, and (d) the period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ff642949",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Algebra\n",
    "# r = p/(1+ecos(theta)\n",
    "# p = a(1-e^2)\n",
    "# so r = a(1-e^2)/(1+ecos(theta))\n",
    "# We know r1, and r2, and theta 1, and theta 2. we dont know a. \n",
    "# After some algebra that cancels out a\n",
    "# e = (r2-r1)/(r1cos(theta1)-r2cos(theta2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "daed14ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "z1 = 1545+r_Earth # km\n",
    "z2 = 852+r_Earth # km \n",
    "theta1 = 126*np.pi/180; # Radians \n",
    "theta2 = 58*np.pi/180; # Radians\n",
    "theta_per = 0; # Radians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "34434066-5c10-475a-b6ae-47642e12b2a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6585.335431899984"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e = (z2-z1)/(z1*np.cos(theta1)-z2*np.cos(theta2)) \n",
    "a = (z1*(1+e*np.cos(theta1)))/(1-e**2) # km\n",
    "z_per = ((a*(1-e**2))/(1+e*np.cos(theta_per)))-r_Earth # km \n",
    "T = np.sqrt(((4*np.pi**2)/mu_Earth)*(a**3)) # Seconds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "120497e3-0845-4392-93d7-54534fea2a35",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "e = 0.08164 \n",
      "z_per = 595.47 km \n",
      "a = 7593.54 km \n",
      "T = 6585.34 seconds\n"
     ]
    }
   ],
   "source": [
    "print(f'e = {e:.5f} \\nz_per = {z_per:.2f} km \\na = {a:.2f} km \\nT = {T:.2f} seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45baccb5-8bad-408b-8384-3a2e36787a93",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3b635967",
   "metadata": {},
   "source": [
    "## 4. Curtis, page 104, Example 2.11\n",
    "\n",
    "An earth orbit has an eccentricity of 0.3, an angular momentum of $ 60,000 km^2/s $, and a true anamoly of $ 120^{\\circ} $.  What are the position vector $ \\vec{\\mathbf{r}} $ and velocity vector $ \\vec{\\mathbf{v}} $ in the perifocal frame of reference?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "75dc5d66",
   "metadata": {},
   "outputs": [],
   "source": [
    "e = 0.3 \n",
    "h = 60000 # km^2/s \n",
    "theta = 120 * np.pi/180; # Radians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d0d525d1-f8dc-488d-ba14-064e391d60f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = (h**2)/mu_Earth # semi-latus rectum\n",
    "r = p/(1+(e*np.cos(theta))) # km\n",
    "a = p/(1-e**2)\n",
    "# Calculate position vector \n",
    "r_vec = Vector((r*np.cos(theta)),(r*np.sin(theta)),0)\n",
    "\n",
    "# Calculate Velocity Vector\n",
    "# After some product rule\n",
    "v_vec = Vector((-(mu_Earth/h)*np.sin(theta)), ((mu_Earth/h)*(e+np.cos(theta))),0);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "51f12a55-250d-4931-a4a0-600eb8f1cd02",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Position Vector \n",
      "x   -5312.712140\n",
      "y    9201.887352\n",
      "z       0.000000\n",
      "Name: component, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(f'Position Vector \\n{r_vec}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5df7e6c7-99b2-47ee-9a45-0031b8fa2968",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Velocity Vector \n",
      "x   -5.753295\n",
      "y   -1.328667\n",
      "z    0.000000\n",
      "Name: component, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(f'Velocity Vector \\n{v_vec}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7fe5f49",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "616cf779",
   "metadata": {},
   "source": [
    "## 5. Curtis, page 104, Example 2.12\n",
    "\n",
    "An earth satellite has the following position and velocity vectors at a given instant:\n",
    "\n",
    "$ \\vec{\\mathbf{r}} = 7000 \\hat{p} + 9000 \\hat{q} \\; (\\mathrm{km}) $\n",
    "\n",
    "$ \\vec{\\mathbf{v}} = -3.3472 \\hat{p} + 9.1251 \\hat{q} \\; (\\mathrm{km/s}) $\n",
    "\n",
    "Calculate the specific angular momentum $ h $, the true anamoly $ \\theta $, and the eccentricity $ e $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "9381a025-602c-4127-9ca1-ceeacab1cc5e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "398600.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu=3.986000e+5 # earth, km^3/s^2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "b7a70981",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = (7000, 9000); #km\n",
    "v = (-3.3472, 9.1251); # km/s\n",
    "h = vector_cross(r,v); # angular momentum\n",
    "rmag = vector_mag(r) # magnitude of position vector\n",
    "vmag = vector_mag(v) # magnitude of velocity vector\n",
    "hs = vector_mag(h); # specific angular momentum\n",
    "mu = 3.986e5; \n",
    "p = hs**2/mu; \n",
    "a = ((2/rmag)-(vmag**2/mu))**-1; #semi-major axis\n",
    "ecc = np.sqrt(-1*((p/a)-1)); # eccentricity\n",
    "ta = np.arctan(r[1]/r[0]); # true anamoly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "bbbf9981",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "94000.5\n",
      "52.125016348901795\n",
      "1.5380010888329392\n"
     ]
    }
   ],
   "source": [
    "print(hs) # specific angular momentum\n",
    "print(ta*180/np.pi) # true anamoly radians\n",
    "print(ecc) # eccentricity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7406940a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9ed6c87-5338-44cf-9561-edbf6eaa029c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0eec1ad-d58c-4b32-af68-23758505cb53",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a0275a-3d9e-46d3-bb8b-ba37e63dcd61",
   "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
}