Vector Approach to Circular Orbits with RK4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bokeh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import math\n",
    "def InitialDataReducedMass(orbitalangle,r0):\n",
    "    phi=orbitalangle\n",
    "    radiusseparation=r0 #separation one hundred\n",
    "    eccentricity=np.zeros(2)\n",
    "    mass=np.ones(2)\n",
    "    masssun=1.989*10**30\n",
    "    lsun=3.828*10**26\n",
    "    massearth=5.9722*10**24\n",
    "    masses=mass #*masssun (natural units)\n",
    "    return phi, radiusseparation,eccentricity, masses\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getxyuveqreducedmass(initdat): #using reduced mass cooordinate system\n",
    "    phi,radiusseparation,eccentricity, masses=initdatreducedmass\n",
    "    #print(orbitalradius, phi, np.cos(phi), np.sin(phi))\n",
    "    metersperAU=1\n",
    "    Gconstant=1\n",
    "    x0=np.zeros(2)\n",
    "    y0=np.zeros(2)\n",
    "    x0[1]=radiusseparation*np.cos(phi)*metersperAU\n",
    "    y0[1]=radiusseparation*np.sin(phi)*metersperAU\n",
    "    if phi==0:\n",
    "        print(\"zero\")\n",
    "        x0[1]=radiusseparation\n",
    "        y0[1]=0\n",
    "    if phi==math.pi:\n",
    "        print(\"pi\")\n",
    "        x0[1]=-radiusseparation\n",
    "        y0[1]=0\n",
    "    z0=np.zeros(2)\n",
    "    print(\"x0\", x0)\n",
    "    print(\"y0\", y0)\n",
    "    \n",
    "    #start at perihelion for both (eliptical, doesn't generalize to three body)\n",
    "    #actually start with circular orbit\n",
    "    ux0=np.zeros(2) #*149597870700\n",
    "    #centrepital force balances gravitational force\n",
    "    reducedmass=np.zeros(2)\n",
    "    masstotal=np.sum(masses)\n",
    "    reducedmass[0]=masstotal\n",
    "    print(masses)\n",
    "    for i in np.arange(1):\n",
    "        j=(i+1)%2 #reverse masses\n",
    "        reducedmass[i+1]=masses[i]*masses[j]/np.sum(masses)\n",
    "    print(reducedmass)\n",
    "    metersperAU=1 #natural units\n",
    "    rphys=np.zeros(len(masses))\n",
    "    for i in np.arange(2):\n",
    "        #rphys is r in the reduced mass system\n",
    "        rphys[i]=radiusseparation #*metersperAU*reducedmass[i]/masses[i]\n",
    "    x0[0]=0.0\n",
    "    x0phys=rphys*x0/radiusseparation\n",
    "    x0phys[0]=0.0\n",
    "    y0phys=rphys*y0/radiusseparation\n",
    "    z0phys=np.zeros(2)\n",
    "    r0=np.sqrt(x0**2+y0**2)\n",
    "    print(rphys)\n",
    "    #G=1\n",
    "    Gconstant=1\n",
    "    F=(Gconstant*masstotal**2/(2*rphys)**2)\n",
    "    print(F)\n",
    "    #centF=reducedmass*v**2/rphys\n",
    "    #centF=accel\n",
    "    v=np.sqrt(Gconstant*masstotal/(rphys))\n",
    "    print(v)\n",
    "    print(\"r0\", r0)\n",
    "    ux0=-v*y0/radiusseparation #[1]\n",
    "    uy0=v*x0/radiusseparation #[1]\n",
    "        \n",
    "    #ux0=-v*y0/r0\n",
    "    #uy0=v*x0/r0 #initial data in y only\n",
    "    #evolve in plane only\n",
    "    #circular orbit\n",
    "    uz0=np.zeros(2)\n",
    "    \n",
    "    a=Gconstant*masstotal/((2*rphys)**2)\n",
    "    ax0=-a*x0/radiusseparation\n",
    "    ay0=-a*y0/radiusseparation\n",
    "    az0=np.zeros(2)\n",
    "    \n",
    "    statevec=[]\n",
    "    avec=[]\n",
    "    for i in np.arange(len(x0phys)):\n",
    "        stateveci=np.array([x0phys[i],y0phys[i],z0phys[i],ux0[i],uy0[i],uz0[i]])\n",
    "        aveci=np.array([ax0[i],ay0[i],az0[i]])\n",
    "        statevec.append(stateveci)\n",
    "        avec.append(aveci)\n",
    "    statevecnp=np.array(statevec)\n",
    "    avecnp=np.array(avec)\n",
    "    return masstotal, reducedmass,statevecnp,avecnp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random,numpy as np\n",
    "initdatreducedmass=InitialDataReducedMass(0,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zero\n",
      "x0 [ 0. 50.]\n",
      "y0 [0. 0.]\n",
      "[1. 1.]\n",
      "[2.  0.5]\n",
      "[50. 50.]\n",
      "[0.0004 0.0004]\n",
      "[0.2 0.2]\n",
      "r0 [ 0. 50.]\n",
      "(2.0, array([2. , 0.5]), array([[ 0. ,  0. ,  0. , -0. ,  0. ,  0. ],\n",
      "       [50. ,  0. ,  0. , -0. ,  0.2,  0. ]]), array([[-0.    , -0.    ,  0.    ],\n",
      "       [-0.0002, -0.    ,  0.    ]]))\n"
     ]
    }
   ],
   "source": [
    "xyuvaeqreducedmass=getxyuveqreducedmass(initdatreducedmass)\n",
    "print(xyuvaeqreducedmass)\n",
    "masstotal, reducedmass,xvec,avec=xyuvaeqreducedmass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0.])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xvec[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0., 50.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "axii=np.zeros(len(xvec[:,0]))\n",
    "rii=np.sqrt(xvec[:,0]**2+xvec[:,1]**2+xvec[:,2]**2)\n",
    "rii=np.sqrt(xvec[:,0]**2+xvec[:,1]**2+xvec[:,2]**2)\n",
    "rii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK4implicit(h,t,xvec,f): #not a finite difference so no step in y\n",
    "    k1= h*f(t,xvec)\n",
    "    k2=h*f(t+h/2, xvec+h*k1/2)\n",
    "    k3=h*f(t+h/2,xvec+h*k2/2)\n",
    "    k4=h*f(t+h,xvec+h*k3)\n",
    "    return t+h, xvec+1/6.*(k1+2.*k2+2.*k3+k4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class OrbitDiffEq:\n",
    "    def __init__(self,reducedmass,xvec,avec,t0):\n",
    "        self.reducedmass=reducedmass\n",
    "        self.xvec=xvec\n",
    "        self.ti=t0\n",
    "        self.avec=avec\n",
    "    def dxidt(self,t,xvec):\n",
    "        return xvec[:,3].transpose()\n",
    "    def dyidt(self,t,xvec):\n",
    "        return xvec[:,4].transpose()\n",
    "    def dzidt(self,t,xvec):\n",
    "        return xvec[:,5].transpose()\n",
    "    def dvxidt(self,t,xvec):\n",
    "        #return axi[m]\n",
    "        axii=np.zeros(len(xvec[:,0]))\n",
    "        rii=np.sqrt(xvec[:,0]**2+xvec[:,1]**2+xvec[:,2]**2)\n",
    "        Gconstant=1 #6.408*10**-11\n",
    "        for k in np.arange(len(rii)):\n",
    "            for j in np.arange(1,len(rii)):\n",
    "                if j==0:\n",
    "                    axii[j]=0.0\n",
    "                if j!=k:\n",
    "                    rreljk=np.abs((xvec[j,0] - xvec[k,0])**2+(xvec[j,1]-xvec[k,1])**2+(xvec[j,2]-xvec[k,2])**2)**(1./2.)\n",
    "                    axii[j]-=Gconstant*self.reducedmass[k]*(xvec[j,0]  - xvec[k,0])/rreljk**3\n",
    "        return axii\n",
    "    def dvyidt(self,t,xvec):\n",
    "        #return axi[m]\n",
    "        ayii=np.zeros(len(xvec[:,0]))\n",
    "        rii=np.sqrt(xvec[:,0]**2+xvec[:,1]**2+xvec[:,2]**2)\n",
    "        Gconstant=1 #6.408*10**-11\n",
    "        for k in np.arange(len(rii)):\n",
    "            for j in np.arange(1,len(rii)):\n",
    "                if j==0:\n",
    "                    ayii[j]=0.0\n",
    "                if j!=k:\n",
    "                    rreljk=np.abs((xvec[j,0] - xvec[k,0])**2+(xvec[j,1]-xvec[k,1])**2+(xvec[j,2]-xvec[k,2])**2)**(1./2.)\n",
    "                    ayii[j]-=Gconstant*self.reducedmass[k]*(xvec[j,1]  - xvec[k,1])/rreljk**3\n",
    "        return ayii\n",
    "    def dvzidt(self,t,xvec):\n",
    "        #return axi[m]\n",
    "        azii=np.zeros(len(xvec[:,0]))\n",
    "        rii=np.sqrt(xvec[:,0]**2+xvec[:,1]**2+xvec[:,2]**2)\n",
    "        Gconstant=1 #6.408*10**-11\n",
    "        for k in np.arange(len(rii)):\n",
    "            for j in np.arange(1,len(rii)):\n",
    "                if j==0:\n",
    "                    azii[j]=0.0\n",
    "                if j!=k:\n",
    "                    rreljk=np.abs((xvec[j,0] - xvec[k,0])**2+(xvec[j,1]-xvec[k,1])**2+(xvec[j,2]-xvec[k,2])**2)**(1./2.)\n",
    "                    azii[j]-=Gconstant*self.reducedmass[k]*(xvec[j,2]  - xvec[k,2])/rreljk**3\n",
    "        return azii\n",
    "    \n",
    "    def dvecdt(self,t,xvec):\n",
    "        avec=np.array([self.dvxidt(self,xvec),self.dvyidt(self,xvec),self.dvzidt(self,xvec)])\n",
    "        avecT=avec.transpose()\n",
    "        self.avec=avecT\n",
    "        dvec2=np.array([self.dxidt(t,xvec),self.dyidt(t,xvec),self.dzidt(self,xvec),avec[0,:],avec[1,:],avec[2,:]])\n",
    "                       \n",
    "        dvec2T=dvec2.transpose()\n",
    "        return dvec2T\n",
    "    def updateINTERNAL(self,xvecii,tii):\n",
    "        self.xvec=xvecii\n",
    "        self.ti=tii\n",
    "        return self\n",
    "    def update(self,xvecii,avecii,tii):\n",
    "        self.xvec=xvecii\n",
    "        self.avec=avecii,\n",
    "        self.ti=tii\n",
    "    def print2D(self):\n",
    "        print(self.reducedmass,self.xvec[1], self.ti)\n",
    "        return self\n",
    "    def list2D(self):\n",
    "        return self.reducedmass,self.xvec,self.avec,self.ti\n",
    "    def timestepRK4ODE(self,step,dt):\n",
    "\n",
    "    \n",
    "        h=dt\n",
    "        #tnew,ynew, intval=RK4(h,t,y,f)\n",
    "        #m represents choices of mass\n",
    "        i=step\n",
    "        \n",
    "        tnew,intvalxvec=RK4implicit(h,self.ti,self.xvec,self.dvecdt)\n",
    "        #xii = intvalx\n",
    "        #tnew,intvalvx=RK4implicit(h,self.ti,self.vxi,self.dvxidt)\n",
    "        #vxii=intvalvx\n",
    "        #tnew,intvaly=RK4implicit(h,self.ti,self.yi,self.dyidt)\n",
    "        #yii = intvaly\n",
    "        #tnew,intvalvy=RK4implicit(h,self.ti,self.vyi,self.dvyidt)\n",
    "        #vyii=intvalvy\n",
    "        #tnew,intvalz=RK4implicit(h,self.ti,self.zi,self.dzidt)\n",
    "        #zii = intvalz\n",
    "        #tnew,intvalvz=RK4implicit(h,self.ti,self.vzi,self.dvzidt)\n",
    "        #vzii=intvalvz\n",
    "        #axii=self.dvxidt(self.ti,np.zeros(np.shape(self.xi)))\n",
    "        #ayii=self.dvyidt(self.ti,np.zeros(np.shape(self.yi)))\n",
    "        #azii=self.dvzidt(self.ti,np.zeros(np.shape(self.zi)))\n",
    " \n",
    "        #print(xii)\n",
    "        self.updateINTERNAL(intvalxvec,tnew)\n",
    "        return reducedmass, self.xvec,self.avec,self.ti"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.  0.5] [50.   0.   0.  -0.   0.2  0. ] 0.0\n",
      "[2.  0.5] [50.   0.   0.  -0.   0.2  0. ] 0.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<__main__.OrbitDiffEq at 0x7f4d24f8fee0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mtotal,reducedmass,xvec,avec=xyuvaeqreducedmass\n",
    "ODEeq= OrbitDiffEq(reducedmass,xvec,avec,0.0)\n",
    "ODEeq.print2D()\n",
    "ODEeq.print2D()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.  0.5] [50.   0.   0.  -0.   0.2  0. ] 0.0\n",
      "[2.  0.5] [ 4.99996000e+01  1.99999467e-01  0.00000000e+00 -7.99997867e-04\n",
      "  1.99998400e-01  0.00000000e+00] 1.0\n",
      "[2.  0.5] [ 4.99984000e+01  3.99995733e-01  0.00000000e+00 -1.59998293e-03\n",
      "  1.99993600e-01  0.00000000e+00] 2.0\n",
      "[2.  0.5] [ 4.9996400e+01  5.9998560e-01  0.0000000e+00 -2.3999424e-03\n",
      "  1.9998560e-01  0.0000000e+00] 3.0\n",
      "[2.  0.5] [ 4.99936001e+01  7.99965867e-01  0.00000000e+00 -3.19986347e-03\n",
      "  1.99974401e-01  0.00000000e+00] 4.0\n",
      "[2.  0.5] [ 4.99900003e+01  9.99933335e-01  0.00000000e+00 -3.99973334e-03\n",
      "  1.99960001e-01  0.00000000e+00] 5.0\n",
      "[2.  0.5] [ 4.99856007e+01  1.19988480e+00  0.00000000e+00 -4.79953921e-03\n",
      "  1.99942403e-01  0.00000000e+00] 6.0\n",
      "[2.  0.5] [ 4.99804013e+01  1.39981707e+00  0.00000000e+00 -5.59926830e-03\n",
      "  1.99921605e-01  0.00000000e+00] 7.0\n",
      "[2.  0.5] [ 4.99744022e+01  1.59972695e+00  0.00000000e+00 -6.39890779e-03\n",
      "  1.99897609e-01  0.00000000e+00] 8.0\n",
      "[2.  0.5] [ 4.99676035e+01  1.79961123e+00  0.00000000e+00 -7.19844490e-03\n",
      "  1.99870414e-01  0.00000000e+00] 9.0\n"
     ]
    }
   ],
   "source": [
    "dt=1. #natural units #*31556926 #seconds per year\n",
    "numsteps=20000\n",
    "mtotal, reducedmass,xvec0,avec0=xyuvaeqreducedmass\n",
    "ODEeq= OrbitDiffEq(reducedmass,xvec0,avec0,0.0)\n",
    "ODEeq.print2D()\n",
    "t=0.0\n",
    "star1x=[]\n",
    "star2x=[]\n",
    "star1a=[]\n",
    "star2a=[]\n",
    "times=[]\n",
    "for i in np.arange(1,numsteps):\n",
    "    reducedmass,xvec,avec,t=ODEeq.timestepRK4ODE(i,dt)\n",
    "    if i< 10:\n",
    "        ODEeq.print2D()\n",
    "    star1x.append(xvec[0])\n",
    "    star2x.append(xvec[1])\n",
    "    star1a.append(avec[0])\n",
    "    star2a.append(avec[1])\n",
    "    times.append(t)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar2x=np.array(star2x)\n",
    "npstar1x=np.array(star1x)\n",
    "npstar2a=np.array(star2a)\n",
    "npstar1a=np.array(star1a)\n",
    "nptimes=np.array(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.io import output_notebook, show\n",
    "from bokeh.plotting import figure\n",
    "output_notebook()\n",
    "# create a new plot with default tools, using figure\n",
    "p = figure(title=\"Circular orbit, RK4\", plot_width=400, plot_height=400)\n",
    "\n",
    "# add a circle renderer with x and y coordinates, size, color, and alpha\n",
    "p.circle(npstar1x[:,0],npstar1x[:,1], size=15, line_color=\"purple\",fill_color=\"purple\", fill_alpha=0.5, legend=\"M\")\n",
    "p.line(npstar2x[:,0],npstar2x[:,1], line_color=\"blue\", legend=\"mu\")\n",
    "\n",
    "\n",
    "p.legend.location = \"bottom_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "show(p) # show the"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%%html\n",
    "<a href=\"https://github.com/sdorsher/ManyBodyProblem/blob/master/reduced%20mass%20circular%20orbit%20rk4.png\"> Circular orbit in the reduced mass system </a> \n",
    "Note that the central heavier mass is stationary in this reference frame with the smaller lighter mass orbits about it in a circular orbit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.io import output_notebook, show\n",
    "from bokeh.plotting import figure\n",
    "output_notebook()\n",
    "# create a new plot with default tools, using figure\n",
    "p = figure(title=\"X-acceleration, Circular orbit, RK4\", plot_width=400, plot_height=400)\n",
    "\n",
    "# add a circle renderer with x and y coordinates, size, color, and alpha\n",
    "p.circle(times,npstar1a[:,0], size=1, line_color=\"purple\",fill_color=\"purple\", fill_alpha=0.5, legend=\"M\")\n",
    "p.circle(times,npstar2a[:,0], size=1, line_color=\"blue\", legend=\"mu\")\n",
    "\n",
    "p.legend.location = \"bottom_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "show(p) # show the"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%%html\n",
    "<a href=\"https://github.com/sdorsher/ManyBodyProblem/blob/master/reduced%20mass%20acceleration%20circular%20RK4.png\">Acceleration in the reduced mass reference frame for a circular orbit </a> \n",
    "\n",
    "In the reduced mass reference frame, the acceleration of the larger mass is zero and the acceleration of the x componenet of the smaller mass is sinusiodal. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RunFirstOrbit(dt,initcondition,printoutput):\n",
    "    masstotal,mass0,xvec, avec=initcondition\n",
    "    x0=xvec[1,0]\n",
    "    y0=xvec[1,1]\n",
    "    lastvalue=[1,1]\n",
    "    thisvalue=[xvec[1,0],xvec[1,1]]\n",
    "    i=1\n",
    "    t=0\n",
    "    while not (lastvalue[0]>0 and lastvalue[1]<0 and thisvalue[0]> 0 and thisvalue[1]>=0):\n",
    "        reducedmass,xvec,avec,t=ODEeq.timestepRK4ODE(i,dt)\n",
    "        #if printoutput:\n",
    "        #    print(xvec[1,0])\n",
    "        #print(ay) #forces should be equal and opposite, but in reduced mass framework accelerations are also equal and opposite\n",
    "        #accelerations should evolve from y to x with time in a sinusoidal manner even in reduced mass framework\n",
    "        #print(ax)\n",
    "        xy=[xvec[1,0],xvec[1,1]]\n",
    "        lastvalue,thisvalue=thisvalue, xy\n",
    "        #lastvalue=thisvalue\n",
    "        #thisvalue[0]=x[1]\n",
    "        #thisvalue[1]=y[1]\n",
    "        #thisvalue[0]=x[1]\n",
    "        #thisvalue[1]=y[1]\n",
    "        if printoutput:\n",
    "            print(xy,lastvalue,thisvalue)\n",
    "        i+=1\n",
    "    #slope=(thisvalue[1]-lastvalue[1])/(thisvalue[0]-lastvalue[0])\n",
    "    #xcrossing=1./slope*(0.-lastvalue[1])+lastvalue[0] #(y1-y2)=m(x1-x2)\n",
    "    xcrossing=(lastvalue[0]+thisvalue[0])/2.\n",
    "    print(dt, thisvalue[0], thisvalue[1], lastvalue[0], lastvalue[1], x0, y0,xcrossing,np.abs(x0-xcrossing))\n",
    "    #interpolate where it crosses zero\n",
    "    #slope=(thisvalue[1]-lastvalue[1])/(thisvalue[0]-lastvalue[0])\n",
    "    #xcrossing=1./slope*(0.-lastvalue[1])+lastvalue[0] #(y1-y2)=m(x1-x2)\n",
    "    return np.abs(x0-xcrossing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt=.1 #*31556926 #seconds per year\n",
    "numsteps=2000\n",
    "mtotal,mass0,xvec,avec=xyuvaeqreducedmass\n",
    "L0arr=[]\n",
    "dtarr=[]\n",
    "dtinit=0.125/4.\n",
    "for j in np.arange(6):\n",
    "    dt=(2**j)*dtinit\n",
    "    dtarr.append(dt)\n",
    "    L0=RunFirstOrbit(dt,xyuvaeqreducedmass,False)\n",
    "    L0arr.append(L0)\n",
    "print(dtarr, L0arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xdat=np.log2(np.array(dtarr))\n",
    "ydat=np.log2(np.array(L0arr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitparam=np.polyfit(xdat,ydat,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitparam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errfn=np.poly1d(fitparam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errfn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.io import output_notebook, show\n",
    "from bokeh.plotting import figure\n",
    "output_notebook()\n",
    "# create a new plot with default tools, using figure\n",
    "p = figure(title=\"L0 error after one orbit, Circular, Vector, RK4\", x_axis_label=\"dt\", y_axis_label=\"L0\", plot_width=400, plot_height=400)\n",
    "\n",
    "# add a circle renderer with x and y coordinates, size, color, and alpha\n",
    "\n",
    "p.circle(np.log2(np.array(dtarr)),np.log2(np.array(L0arr)), line_color=\"blue\", legend=\"L0\")\n",
    "p.line(xdat,errfn(xdat),line_color=\"green\", legend=\"Fit, dt^4.1\")\n",
    "p.legend.location = \"bottom_right\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "show(p) # show the"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%%html \n",
    "<a href=\"https://github.com/sdorsher/ManyBodyProblem/blob/master/L0%20error%20Circular%20Vector%20RK4%20dt%204p1.png\"> L0 error after one circular orbit in the reduced mass sytem </a>\n",
    "\n",
    "L0 error is the absolute value of the difference between the x coordinate at t=0 and the x coordinate after one orbit. I have plotted the behavior of the L0 error as a function of the timestep dt. For the RK4 numerical differential equation solution method, the error should go as dt^4 or better. It doe, since it goes as dt^4.1. "
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}