Ex_3.5_Solve_Universal_Kepler_Equation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0bf517a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1045349349268365"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Solving Universal Kepler Equation\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def S(z):\n",
    "    return (1/6) - z/120 + z**2/5040 - z**3/362880 + z**4/39916800\n",
    "\n",
    "def C(z):\n",
    "    return (1/2) - z/24 + z**2/720 - z**3/40320 + z**4/3628800\n",
    "    \n",
    "\n",
    "def Universal_F(Dt,r0,vr0,alpha,mu):\n",
    "    \n",
    "    X = np.sqrt(mu)*np.abs(alpha)*Dt\n",
    "    z = alpha*X\n",
    "    Universal_f= (r0*vr0)/(np.sqrt(mu)) * X**2*C(z)+(1-alpha*r0)*X**3 * S(z)+r0*X-np.sqrt(mu)*Dt\n",
    "    Universal_f_dot = (r0*vr0)/(np.sqrt(mu)) * X*(1- (alpha*X**2* S(z)))+(1-alpha*r0)*X**2 * C(z)+r0\n",
    "\n",
    "    return(Universal_f,Universal_f_dot)\n",
    "\n",
    "def Universal_Kepler(Dt,r0,vr0,alpha,mu):\n",
    "    ratio=1\n",
    "    X = np.sqrt(mu)*np.abs(alpha)*Dt\n",
    "    z = alpha*X\n",
    "    while ratio > 10**(-12):\n",
    "\n",
    "        ratio=Universal_F(Dt,r0,vr0,alpha,mu)[0]/Universal_F(Dt,r0,vr0,alpha,mu)[1]\n",
    "        X = X-ratio\n",
    "    return(X)\n",
    "\n",
    "Universal_Kepler(1,0.5,0.5,0.5,0.5)"
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}