{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modeling A Particle In the Gravity Field of an Oblate Planet\n",
"### Cee Gould"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.mplot3d import axes3d\n",
"from matplotlib import cm\n",
"from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
"%matplotlib osx \n",
"from matplotlib import cm\n",
"from colorspacious import cspace_converter\n",
"from collections import OrderedDict\n",
"\n",
"cmaps = OrderedDict()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$V(r,u) = \\frac{GM}{r} [1 - \\sum (\\frac{a}{r})^{2n} J_{2n} P_{2n}(\\mu)] $ , only using n = 1\n",
"
\n",
"
\n",
"*$P_{n2} = \\frac{1}{2}(3x^{2} - 1)$\n",
"
\n",
"*$J_{2n} = 2$\n",
"
\n",
"$\\mu = \\frac{z}{r}$\n",
"
\n",
"$V(r,u) = \\frac{GM}{r} [1 - (\\frac{a}{r})^{2} * J_{2} * (\\frac{1}{2} (3x^{2} - 1)*\n",
"(\\mu)] $ , only using n = 1\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"#----------// User Inputs // ----------\n",
"x = 2\n",
"y = -1.5\n",
"z = 0.5\n",
"vx = 0\n",
"vy = 12\n",
"vz = 0\n",
"b = 0.5 #short radius\n",
"a = 2 #long radius\n",
"ms = 3e+6 #planet mass\n",
"mp = 1 #particle mass"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# //---------- Oblate Planet Function ----------\n",
"def V(t, X):\n",
" \n",
" x = X[0]\n",
" y = X[1]\n",
" z = X[2]\n",
" \n",
" r = X[0:3]\n",
" v = X[3:6]\n",
" \n",
" # // PERIODIC ORBITS ABOUT AN OBLATE SPHEROID, MACMILLAN 1910\n",
" # https://www.ams.org/journals/tran/1910-011-01/S0002-9947-1910-1500856-2/S0002-9947-1910-1500856-2.pdf\n",
" \n",
" dVx = - (k**2) * ms * mp / np.linalg.norm(r)**3 * (1 + (3/10)*b**2*(x**2 + y**2 - 4*z**2)/ np.linalg.norm(r)**4 * mu**2 ) * x\n",
" dVy = - (k**2) * ms * mp / np.linalg.norm(r)**3 * (1 + (3/10)*b**2*(x**2 + y**2 - 4*z**2)/ np.linalg.norm(r)**4 * mu**2 ) * y\n",
" dVz = - (k**2) * ms * mp / np.linalg.norm(r)**3 * (1 + (3/10)*b**2*(3*(x**2 + y**2) - 2*z**2)/ np.linalg.norm(r)**4 * mu**2 ) * z\n",
" \n",
" #// First attempt\n",
"# dVx = -((G*M)/(np.linalg.norm(r)**2))*x + 3*((6*M*a**2)/(np.linalg.norm(r)**4))*x\n",
"# dVy = -((G*M)/(np.linalg.norm(r)**2))*y + 3*((6*M*a**2)/(np.linalg.norm(r)**4))*y\n",
"# dVz = -((G*M)/(np.linalg.norm(r)**2))*z - ((30*G*M*a**2*z**2)/np.linalg.norm(r)**6) + 3*((6*M*a**2)/(np.linalg.norm(r)**4))*z\n",
" \n",
" #F = -G * ms * mp /np.linalg.norm(r)**3 * r #Newtons Method\n",
" \n",
" \n",
" F = np.array([dVx, dVy, dVz])\n",
" acceleration = F / mp\n",
" \n",
" return np.concatenate((v, acceleration))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# // ----------Constants ----------//\n",
"year = 3.154e+7 #s\n",
"au = 1.496e+8 #km\n",
"\n",
"k = 0.01720209895 #rad/day\n",
"#k = np.sqrt(G*Msun)\n",
"\n",
"G = k**2 #1e-6\n",
"M = ms + mp\n",
"\n",
"#G = 6.674e-20 #km^3/kg s^2\n",
"\n",
"grav_param = G * M\n",
"\n",
"a_earth = 6371 #km # equatorial radius\n",
"b_earth = 6356.8 #km #short length of earth #polar radius \n",
"m_earth = 5.972e24 #kg\n",
"m_moon = 1 #7.35e22 #kg"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"r = [ 2. -1.5 0.5]\n",
"v = [ 0 12 0]\n",
"mu = 0.75\n"
]
}
],
"source": [
"r = np.array([x,y,z])\n",
"v = np.array([vx,vy,vz])\n",
"\n",
"print('r = ', r)\n",
"print('v = ', v)\n",
"\n",
"mu = 1 - b/a #oblateness of the spheroid 0 < mu < 1. mu = 0 is a sphere. \n",
"print(\"mu = \", mu)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"e = 0.755659949446618 elliptical orbit\n",
"E = -276.19906784137 bound\n",
"a = 1.60705995083017\n",
"p = 2.0372677831379082\n",
"i = 0.24497866312686378\n",
"Omega = 358.4292036732051\n",
"argp = 357.53678973688943\n",
"initial nu = 2.878193519700783 degrees\n"
]
}
],
"source": [
"#/// ----------Calculate orbital elements ----------/// \n",
"#Fundamentals of Astrodynamics and Applications, by Vallado, 2007\n",
"\n",
"h = np.cross(r,v) # angular momentum\n",
"K = np.array([0, 0, 1])\n",
"nhat=np.cross(K,h) #node vector\n",
"\n",
"#eccentricity vector\n",
"evec = ((np.linalg.norm(v)**2 - grav_param /np.linalg.norm(r))*r - np.dot(r,v)*v)/ grav_param \n",
"# evec = (np.linalg.norm(abs(v))**2/mu - 1/np.linalg.norm(abs(r)))*r - (np.dot(r,v)*v / mu) \n",
"\n",
"#evec = np.array([e_p9, 0, 0])\n",
"e = np.linalg.norm(evec)\n",
"\n",
"if e == 0:\n",
" orbit = 'circular orbit'\n",
"if e > 0 and e < 1:\n",
" orbit = 'elliptical orbit'\n",
"if e == 1:\n",
" orbit ='parabolic orbit'\n",
"if e > 1:\n",
" orbit = 'hyperbolic orbit'\n",
"print('e = ', e, orbit)\n",
"\n",
"energy = np.linalg.norm(v)**2/2 - grav_param /np.linalg.norm(r) #Specific mechanical energy\n",
"if energy < 0: #absolute value of potential energy is larger than kinetic\n",
" bound = 'bound'\n",
"if energy > 0: #kinetic energy is larger than the absolute value of the potential energy\n",
" bound = 'unbound'\n",
"print('E = ', energy, bound)\n",
"\n",
"if abs(e) != 1: \n",
" a = - grav_param/(2*energy)\n",
" #p = a*(1-e**2)\n",
" p = np.sqrt(abs(a)**3)\n",
"else:\n",
" #p = np.linalg.norm(h)**2/mu\n",
" a = 'inf'\n",
" p = 'inf'\n",
" \n",
"print('a = ', a) #semi major axis\n",
"\n",
"\n",
"print('p = ', p) #period\n",
"\n",
"i = np.arccos(h[2]/np.linalg.norm(h)) #inclination\n",
"print('i = ', i) \n",
"\n",
"Omega = np.arccos(nhat[0]/np.linalg.norm(nhat)) #swivel: the angle from the principal direction to the ascending node, 0° ≤ Ω < 360°\n",
"if nhat[1] < 0:\n",
" Omega = 360-Omega\n",
"print('Omega = ', Omega)\n",
"\n",
"argp = np.arccos(np.dot(nhat,evec)/(np.linalg.norm(nhat)*e)) # argument of perigee, ω: 0° ≤ ω < 360\n",
"if evec[2]<0:\n",
" argp = 360-argp\n",
"print('argp = ', argp)\n",
"\n",
"#nu_0 = nu_p9\n",
"nu_0 = np.arccos(np.dot(evec,r)/(e*np.linalg.norm(r))) # True anomaly, ν, is the angle along the orbital path from perigee to the spacecraft’s position vector r. 0 <= v < 360°\n",
"if np.dot(r,v)<0:\n",
" nu = 360 - nu_0\n",
"print('initial nu = ', nu_0, 'degrees') # changes with time, location of the spacecraft in its orbit"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"#// ----------Set Time ----------//\n",
"dt = .001\n",
"tmax = 4.0\n",
"tt = []\n",
"t = 0"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# //---------- RUN THE CODE!! ----------// \n",
"xt = []\n",
"yt = []\n",
"zt = []\n",
"vx = []\n",
"vy = []\n",
"vz = []\n",
"\n",
"X = np.concatenate((r,v))\n",
"\n",
"%matplotlib osx \n",
"\n",
"while(t < tmax):\n",
"\n",
" r = X[0:3]\n",
" v = X[3:6]\n",
" \n",
" \n",
" xt.append(X[0])\n",
" yt.append(X[1])\n",
" zt.append(X[2])\n",
" \n",
" \n",
" #X = X + V(X) * dt #EULERS METHOD\n",
" \n",
" f1 = V(t ,X ) #4th order RUNGE KUTTA METHOD\n",
" f2 = V(t+dt/2.0,X+f1*dt/2.0)\n",
" f3 = V(t+dt/2.0,X+f2*dt/2.0)\n",
" f4 = V(t+dt ,X+f3*dt )\n",
"\n",
" X = X + (f1 + 2.0*f2 + 2.0*f3 + f4) / 6.0 * dt\n",
" \n",
" \n",
" tt.append(t)\n",
" t += dt\n",
" \n",
" #// -----------Real Time Plot --------------//\n",
" \n",
" %matplotlib osx \n",
" fig = plt.figure(figsize=(5,5))\n",
" \n",
" plt.subplot(2, 1, 1)\n",
" plt.plot(yt, xt, 'r-')\n",
" plt.xlabel('Y axis')\n",
" plt.ylabel('X axis')\n",
" plt.plot(0,0,'*',mfc='w',ms=10)\n",
" plt.gca().set_aspect('equal', adjustable='box')\n",
" \n",
" \n",
" \n",
" fig.add_subplot(2, 1, 2)\n",
" plt.plot(yt, zt, 'r-')\n",
" plt.xlabel('Y axis')\n",
" plt.ylabel('Z axis')\n",
" plt.plot(0,0,'*',mfc='w',ms=10)\n",
" plt.gca().set_aspect('equal', adjustable='box')\n",
" \n",
" plt.draw()\n",
" plt.show()\n",
" plt.pause(0.0005) \n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"#----------// time plots //----------\n",
"\n",
"plt.subplot(3, 1, 1)\n",
"\n",
"plt.plot(tt,xt, 'r-')\n",
"plt.title('Height of X over time')\n",
"\n",
"fig.add_subplot(3, 1, 2)\n",
"plt.plot(tt,yt, 'r-')\n",
"plt.title('Height of Y over time')\n",
"\n",
"fig.add_subplot(3, 1, 3)\n",
"plt.plot(tt,zt, 'r-')\n",
"plt.title('Height of Z over time')\n",
"\n",
"\n",
"plt.show()\n",
"plt.pause(5) \n",
"plt.close()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"#---------- // 3D Plot // ----------\n",
"\n",
"# // Create 3D plot\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"ax.set_aspect(\"equal\")\n",
"ax.set_xlabel('X axis')\n",
"ax.set_ylabel('Y axis')\n",
"ax.set_zlabel('Z axis')\n",
"\n",
"# ax.set_xlim(-1,1)\n",
"# ax.set_ylim(-1, 1)\n",
"# ax.set_zlim(-1, 1)\n",
"\n",
"#// Orbit //\n",
"XT = np.array(xt)\n",
"YT = np.array(yt)\n",
"ZT = np.array(zt)\n",
"orbit = ax.plot(XT, YT, ZT, color = 'b') \n",
" \n",
"#// Sphere //\n",
"\n",
"u, v = np.mgrid[0:2*np.pi:200j, 0:np.pi:100j] \n",
"x = 0.1 * a * np.cos(u)*np.sin(v)\n",
"y = 0.1* a * np.sin(u)*np.sin(v)\n",
"z = 0.1 * b * np.cos(v)\n",
"planet = ax.plot_surface(x, y, z, color=\"b\", cmap='gist_earth')\n",
"\n",
"\n",
"# // Rotating Plot //\n",
"\n",
"for angle in range(0, 360):\n",
" orbit = ax.plot(XT, YT, ZT, color = 'b') \n",
" planet = ax.plot_surface(x, y, z, color=\"b\", cmap='gist_earth')\n",
" \n",
" ax.autoscale_view(True,True,True) \n",
" ax.view_init(20, angle)\n",
" plt.pause(.00001)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}