{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Pertubation of Neptune on Uranus's Angular Velocity

\n", "

Neptune was discovered as the something that was perturbing the angular velocity of Uranus. Here we solve for and visualize that pertubation assuming the planetary orbits remain circular. When Neptune is in front it speeds up the angular velocity of \n", " Uranus, and when Neptune is behind, it slows down the angular velocity. \n", " \n", " \n", "Calling $\\bf{F_T}({\\bf r_{su}},{\\bf r_{nu}})$ the sum of the force of the Sun on Uranus and the force of Neptune on Uranus, then:\n", "\n", "\\begin{eqnarray*}\n", "\\bf{F_T}({\\bf r_{su}},{\\bf r_{nu}})& =& {\\bf F_{su}}({\\bf r_{su}}) + {\\bf F_{nu}}({\\bf r_{nu}})\n", "\\\\ \n", "&=& \\frac{G M_s m_u}{\\vert {\\bf r_{su}\\vert^3}}\\,{\\bf r_{su}}+\n", " \\frac{G m_n m_u}{\\vert {\\bf r_{nu}\\vert^3}}\\,{\\bf r_{nu}}\n", "\\end{eqnarray*}" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/html": [ "

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n", " by RH Landau, MJ Paez, and CC Bordeianu (deceased).\n", " Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n", " C Bordeianu, Univ Bucharest, 2021. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# UranusNeptune.ipnby: Pertubation of neptune on Uranus's angular velocity\n", "\n", "from vpython import *\n", "import numpy as np\n", "\n", "scene = canvas(width=600,height=600, \n", "\ttitle = 'White Neptune & Black Uranus', range=40)\n", "sun = sphere(pos=vec(0,0,0), radius=2, color=color.yellow)\n", "escenau = graph(x=600,width=400,height=400, \n", "\ttitle='Pertubation of Uranus Angular Velocity')\n", "graphu = gcurve(color=color.cyan)\n", "escenan = graph(x=800,y=400,width=400,height=400)\n", "graphn = gcurve(color=color.white)\n", "rfactor = 1.8e-9\n", "G = 4*pi*pi # in units T in years, R AU, Msun=1\n", "mu = 4.366244e-5 # mass Uranus in solar masses\n", "M = 1.0 # mass Sun\n", "mn = 5.151389e-5 # Neptune mass in solar masses\n", "du = 19.1914 # distance Uranus Sun in AU\n", "dn = 30.0611 # distance Neptune sun in AU\n", "Tur = 84.0110 # Uranus Orbital Period yr\n", "Tnp = 164.7901 # Neptune Orbital Period yr\n", "omeur = 2*pi/Tur # Uranus angular velocity (2pi/T)\n", "omennp = 2*pi/Tnp # Neptune angular velocity\n", "omreal = omeur\n", "urvel = 2*pi*du/Tur # Uranus orbital velocity UA/yr\n", "npvel = 2*pi*dn/Tnp # Neptune orbital velocity UA/yr\n", "# 1 Uranus at lon 2 gr 16 min sep 1821\n", "radur = (205.64)*pi/180. # to radians in 1690 -wrt x-axis\n", "urx = du*cos(radur) # init x- pos ur. in 1690\n", "ury = du*sin(radur) # init y-pos ur in 1690\n", "urvelx = urvel*sin(radur)\n", "urvely = -urvel*cos(radur)\n", "# 1690 Neptune at long. \n", "radnp = (288.38)*pi/180. # 1690 rad neptune wrt x-axis\n", "uranus = sphere(pos=vec(urx,ury,0), radius=0.5,color=vec(.88,1,1), \n", "\tmake_trail=True)\n", "urpert = sphere(pos=vec(urx,ury,0), radius=0.5,color=vec(.88,1,1), \n", "\tmake_trail=True)\n", "fnu = arrow(pos=uranus.pos,color=color.orange,axis=vector(0,4,0))\n", "npx = dn*cos(radnp) #init coord x neptune 1690\n", "npy = dn*sin(radnp) # y\n", "npvelx = npvel*sin(radnp)\n", "npvely = -npvel*cos(radnp)\n", "neptune = sphere(pos=vec(npx,npy,0), radius=0.4,color=color.cyan, \n", "\tmake_trail=True)\n", "fun = arrow(pos=neptune.pos,color=color.orange,axis=vector(0,-4,0))\n", "nppert = sphere(pos=vec(npx,npy,0), radius=0.4,color=color.white, make_trail=True)\n", "velour = vector(urvelx,urvely,0) #initial vector velocity Uranus\n", "velnp = vector(npvelx,npvely,0) #initial vector velocity Neptune\n", "dt = 0.5 # time increment in terrestrial year\n", "r = vector(urx,ury,0) # initial position Uranus wrt Sun\n", "rnp = vector(npx,npy,0) # initial position Neptune wrt Sun\n", "veltot = velour\n", "veltotnp = velnp\n", "rtot = r\n", "rtotnp = rnp\n", "\n", "def ftotal(r,rnp,i): # i==1 Uranus i==2 Neptune\n", " Fus = -G*M*mu*r/(du**3) # Force sun over URANUS\n", " Fns = -G*M*mn*rnp/(dn**3) # Force Sun over NEPTUNE\n", " dnu = mag(rnp-r) # distance Neptune-Uranus\n", " Fnu = -G*mu*mn*(rnp-r)/(dnu**3) # force N on U\n", " Fun = -Fnu # force uranus on Neptune\n", " Ftotur = Fus+Fnu # total force on U (sun + N)\n", " Ftotnp = Fns+Fun # On Neptune F sun +F urn\n", " if i==1: return Ftotur\n", " else: return Ftotnp\n", " \n", "def rkn(r,veltot,rnp,m,i): # on Neptune\n", " k1v = ftotal(r,rnp,i)/m\n", " k1r = veltot\n", " k2v = ftotal(r,rnp+0.5*k1r*dt,i)/m\n", " k2r = veltot+0.5*k2v*dt\n", " k3v = ftotal(r,rnp+0.5*k2r*dt,i)/m\n", " k3r = veltot+0.5*k3v*dt\n", " k4v = ftotal(r,rnp+k3r*dt,i)/m\n", " k4r = veltot+k4v*dt\n", " veltot = veltot+(k1v+2*k2v+2*k3v+k4v)*dt/6.0\n", " rnp = rnp+(k1r+2*k2r+2*k3r+k4r)*dt/6.0\n", " return r,veltot\n", "\n", "def rk(r,veltot,rnp,m,i): # on Uranus\n", " k1v = ftotal(r,rnp,i)/m\n", " k1r = veltot\n", " k2v = ftotal(r+0.5*k1r*dt,rnp,i)/m\n", " k2r = veltot+0.5*k2v*dt\n", " k3v = ftotal(r+0.5*k2r*dt,rnp,i)/m\n", " k3r = veltot+0.5*k3v*dt\n", " k4v = ftotal(r+k3r*dt,rnp,i)/m\n", " k4r = veltot+k4v*dt\n", " veltot = veltot+(k1v+2*k2v+2*k3v+k4v)*dt/6.0\n", " r = r+(k1r+2*k2r+2*k3r+k4r)*dt/6.0\n", " return r,veltot\n", "\n", "for i in np.arange(0,320):# estaba 1240\n", " rate(10)\n", " rnewu,velnewu = rk(r,velour,rnp,mu,1) # uranus\n", " rnewn,velnewn = rkn(rnp,velnp,r,mn,2) # neptune\n", " r = rnewu # uranus position\n", " velour = velnewu # uranus velocity\n", " du = mag(r)\n", " omeur = mag(velour)/du # new abgykar velocity of uranus\n", " degr = 205.64*pi/180- omeur*i*dt # angular position uranus\n", " rnp = rnewn # neptune pos\n", " velnp = velnewn # neptune pos\n", " dn = mag(rnp)\n", " omenp = mag(velnp)/dn\n", " radnp = radnp-dt*omenp # radians neptune\n", " npx = dn*cos(radnp) \n", " npy = dn*sin(radnp)\n", " rnp = vector(npx,npy,0) # neptune position \n", " deltaomgs = -omeur+omreal\n", " graphu.plot(pos=(i,deltaomgs*180/pi*3600))\n", " urpert.pos = r\n", " fnu.pos = urpert.pos # position of arrow on uranus\n", " dnu = mag(rnp-r) # distance Neptune-Uranus\n", " fnu.axis = 75*norm(rnp-r)/dnu # axes the arrow over uranus\n", " neptune.pos = rnp # radiovector Neptune\n", " fun.pos = neptune.pos\n", " fun.axis = -fnu.axis # arrow on neptune" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }