TwoStarsOnePlanet_convergenceTest-Esplitvst_ELvst.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# THIS NOTEBOOK GENERATES DATA FOR CONVERGENCE OVER TIME, PLOTTED IN EsplitEL vs t convergence plot.ipynb."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of a reduced three body problem where there is a central binary star system and an outer planet, I see precession of the inner star system's orbit about a mutual center of mass. I also see energy transfer between the inner star system and the outer planet in the form of loosening of the inner coils and inspiraling of the outer planet. I see less of this when the planet is further out and less massive. I also see that the Jacobi integral is the only conserved quantity, implying that this is in principle possible, but that symplectice methods are usually used to get numerical stability for the two body problem (and this is three). So I am not yet sure what the actual source of this energy and angular momentum transfer is. I believe the source of the precession and center of mass drift is in the initial conditions. The acceleration of the center of mass and wobble in it is likely due to numerical instabliity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ellipsefns\n",
    "#getfocus, getsemiminor, theoreticalx, theoreticaly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import initdatuneq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rundir='ConvergenceTestwithKEPE/'\n",
    "runoutprefix='Run4/twostarsoneplanet'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "rad0=[10,1000]\n",
    "ecc0=[0.,0.] #set ecc = 0 for circle, ecc < 1 for ellepse, = 1 for parabola, >1 for hyperbola\n",
    "theta0= [0.0,0.] #math.pi/6. #math.pi/6.\n",
    "mass0_1=[100.0,100.0]\n",
    "mass0_2=1.0\n",
    "numsteps=200000 #2000000  #200000 #10000 #50000\n",
    "#numsteps=6000*10\n",
    "dt=2*2*.1 #2*2*.1 #/2./2./2.\n",
    "outputevery=100 #1\n",
    "#rad0=[10,5000]\n",
    "#ecc0=[0.,0.] #set ecc = 0 for circle, ecc < 1 for ellepse, = 1 for parabola, >1 for hyperbola\n",
    "#theta0= [0.0,0.] #math.pi/6. #math.pi/6.\n",
    "#mass0_1=[100.0,100.0]\n",
    "#mass0_2=.1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "angle 0.0\n",
      "0.0\n",
      "angle 0.0\n",
      "0.0\n",
      "(array([3.14159265, 0.        ]), 0.0, 10, 0.0, array([100., 100.]))\n",
      "(array([3.14159265, 0.        ]), 0.0, 1000, 0.0, array([200.,   1.]))\n"
     ]
    }
   ],
   "source": [
    "initdatinner=initdatuneq.InitialDataEqualMassConic(rad0[0],ecc0[0],theta0[0],mass0_1[0],mass0_1[0])\n",
    "initdatouter=initdatuneq.InitialDataEqualMassConic(rad0[1],ecc0[1],theta0[1],np.sum(np.array(mass0_1)),mass0_2)\n",
    "print(initdatinner)\n",
    "print(initdatouter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import getxyuvtwouneq_apper #unequal mass version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "useperihelion=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initdat 1000 [3.14159265 0.        ] 0.0 [-1.  1.] [1.2246468e-16 0.0000000e+00]\n",
      "coordsep 1000\n",
      "pi\n",
      "zero\n",
      "[-1000.  1000.]\n",
      "[0. 0.]\n",
      "circle\n",
      "r0 2000.0\n",
      "(array([200.,   1.]), array([[-4.97512438e+00,  0.00000000e+00,  0.00000000e+00,\n",
      "        -0.00000000e+00, -2.23049868e-03,  0.00000000e+00],\n",
      "       [ 9.95024876e+02,  0.00000000e+00,  0.00000000e+00,\n",
      "        -0.00000000e+00,  4.46099737e-01,  0.00000000e+00]]), array([[ 1.e-06, -0.e+00,  0.e+00],\n",
      "       [-1.e-06, -0.e+00,  0.e+00]]))\n"
     ]
    }
   ],
   "source": [
    "xyuvouter=getxyuvtwouneq_apper.getxyuvequneq(initdatouter,useperihelion)\n",
    "print(xyuvouter)\n",
    "massesouter,xvecouter,avecouter=xyuvouter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 6)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(xvecouter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initdat 10 [3.14159265 0.        ] 0.0 [-1.  1.] [1.2246468e-16 0.0000000e+00]\n",
      "coordsep 10\n",
      "pi\n",
      "zero\n",
      "[-10.  10.]\n",
      "[0. 0.]\n",
      "circle\n",
      "r0 20.0\n",
      "(array([100., 100.]), array([[-5.        ,  0.        ,  0.        , -0.        , -2.23606798,\n",
      "         0.        ],\n",
      "       [ 5.        ,  0.        ,  0.        , -0.        ,  2.23606798,\n",
      "         0.        ]]), array([[ 1., -0.,  0.],\n",
      "       [-1., -0.,  0.]]))\n"
     ]
    }
   ],
   "source": [
    "xyuvinner=getxyuvtwouneq_apper.getxyuvequneq(initdatinner,useperihelion)\n",
    "print(xyuvinner)\n",
    "massesinner,xvecinner,avecinner=xyuvinner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def weavexyuvtwostarsoneplanet(massesinner,massesouter,xvecinner,xvecouter, avecinner,avecouter):\n",
    "    masses=np.zeros(3)\n",
    "    xvec=np.zeros([3,6])\n",
    "    avec=np.zeros([3,3])\n",
    "    masses[0:2]=massesinner\n",
    "    masses[2]=(massesouter[1])\n",
    "    xvec[0:2,:]=xvecinner\n",
    "    xvec[2,:]=(xvecouter[1,:])\n",
    "    avec[0:2,:]=avecinner\n",
    "    avec[2,:]=avecouter[1,:]\n",
    "    return masses,xvec,avec\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "masses,xvec,avec=weavexyuvtwostarsoneplanet(massesinner,massesouter,xvecinner,xvecouter,avecinner,avecouter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([100., 100.,   1.])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "masses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-5.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "        -0.00000000e+00, -2.23606798e+00,  0.00000000e+00],\n",
       "       [ 5.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "        -0.00000000e+00,  2.23606798e+00,  0.00000000e+00],\n",
       "       [ 9.95024876e+02,  0.00000000e+00,  0.00000000e+00,\n",
       "        -0.00000000e+00,  4.46099737e-01,  0.00000000e+00]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xvec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.e+00, -0.e+00,  0.e+00],\n",
       "       [-1.e+00, -0.e+00,  0.e+00],\n",
       "       [-1.e-06, -0.e+00,  0.e+00]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "avec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "initdatinnerlist=[]\n",
    "for item in initdatinner:\n",
    "    if np.isscalar(item):\n",
    "        initdatinnerlist.append(item)\n",
    "    else:\n",
    "        initdatinnerlist.append(list(item))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "initdatouterlist=[]\n",
    "for item in initdatouter:\n",
    "    if np.isscalar(item):\n",
    "        initdatouterlist.append(item)\n",
    "    else:\n",
    "        initdatouterlist.append(list(item))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[3.141592653589793, 0.0], 0.0, 1000, 0.0, [200.0, 1.0]]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initdatouterlist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "paramsdict = {}\n",
    "paramsdict['rad0']=list(rad0)\n",
    "paramsdict['ecc0']=list(ecc0)\n",
    "paramsdict['theta0']=list(theta0)\n",
    "paramsdict['masses']=list(masses)\n",
    "paramsdict['numsteps']=numsteps\n",
    "paramsdict['outputevery']=outputevery\n",
    "paramsdict['dt']=dt\n",
    "paramsdict['initdatinner']=initdatinnerlist\n",
    "paramsdict['initdatouter']=initdatouterlist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'rad0': [10, 1000],\n",
       " 'ecc0': [0.0, 0.0],\n",
       " 'theta0': [0.0, 0.0],\n",
       " 'masses': [100.0, 100.0, 1.0],\n",
       " 'numsteps': 200000,\n",
       " 'outputevery': 100,\n",
       " 'dt': 0.4,\n",
       " 'initdatinner': [[3.141592653589793, 0.0], 0.0, 10, 0.0, [100.0, 100.0]],\n",
       " 'initdatouter': [[3.141592653589793, 0.0], 0.0, 1000, 0.0, [200.0, 1.0]]}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "paramsdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "json = json.dumps(paramsdict)\n",
    "jsonfile=open(rundir+runoutprefix+'params.json',\"w\")\n",
    "jsonfile.write(json)\n",
    "jsonfile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyuvatwostarsplanet=masses,xvec,avec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import RK4implicit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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+k1/2)\n",
    "    k3=h*f(t+h/2,xvec+k2/2)\n",
    "    k4=h*f(t+h,xvec+k3)\n",
    "    return t+h, xvec+1/6.*(k1+2.*k2+2.*k3+k4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import OrbitDiffEq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[100. 100.   1.] [[-5.00000000e+00  0.00000000e+00  0.00000000e+00 -0.00000000e+00\n",
      "  -2.23606798e+00  0.00000000e+00]\n",
      " [ 5.00000000e+00  0.00000000e+00  0.00000000e+00 -0.00000000e+00\n",
      "   2.23606798e+00  0.00000000e+00]\n",
      " [ 9.95024876e+02  0.00000000e+00  0.00000000e+00 -0.00000000e+00\n",
      "   4.46099737e-01  0.00000000e+00]] 0.0\n",
      "[100. 100.   1.] [[-5.00000000e+00  0.00000000e+00  0.00000000e+00 -0.00000000e+00\n",
      "  -2.23606798e+00  0.00000000e+00]\n",
      " [ 5.00000000e+00  0.00000000e+00  0.00000000e+00 -0.00000000e+00\n",
      "   2.23606798e+00  0.00000000e+00]\n",
      " [ 9.95024876e+02  0.00000000e+00  0.00000000e+00 -0.00000000e+00\n",
      "   4.46099737e-01  0.00000000e+00]] 0.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<OrbitDiffEq.OrbitDiffEq at 0x7fdfe120f3d0>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#masses, xvec,avec=xyuvaequneq\n",
    "ODEeq= OrbitDiffEq.OrbitDiffEq(masses,xvec,avec,0.0)\n",
    "ODEeq.print2D()\n",
    "ODEeq.print2D()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(xyuvaequneq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timestep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timestepNonAdapt3star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "import parabfns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ellipsefns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hyperbfns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Estability3,Lstability3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SimpleE3(xvec1,xvec2,xvec3,masses):\n",
    "    pos0=xvec1[0:3]\n",
    "    pos1=xvec2[0:3]\n",
    "    pos2=xvec3[0:3]\n",
    "    vel0=xvec1[3:]\n",
    "    vel1=xvec2[3:]\n",
    "    vel2=xvec3[3:]\n",
    "    p0=masses[0]*vel0\n",
    "    p1 = masses[1]*vel1\n",
    "    p2 = masses[2]*vel2\n",
    "    angL0 = np.cross(pos0,p0)\n",
    "    angL1=np.cross(pos1,p1)\n",
    "    angL2=np.cross(pos2,p2)\n",
    "    KE0 = np.sum(np.square(p0))/2/masses[0]\n",
    "    KE1=np.sum(np.square(p1))/2/masses[1]\n",
    "    KE2=np.sum(np.square(p2))/2/masses[2]\n",
    "    r01=np.sqrt(np.sum(np.square(pos0-pos1)))\n",
    "    r02=np.sqrt(np.sum(np.square(pos0-pos2)))\n",
    "    r12=np.sqrt(np.sum(np.square(pos1-pos2)))\n",
    "    PE01=-masses[0]*masses[1]/r01\n",
    "    PE02=-masses[0]*masses[2]/r02\n",
    "    PE12=-masses[1]*masses[2]/r12\n",
    "    PE=PE01+PE02+PE12\n",
    "    ptot=p0+p1+p2\n",
    "    angLtot=angL0+angL1+angL2\n",
    "    Etot=KE0+KE1+KE2+PE01+PE02+PE12\n",
    "    return ptot,angLtot,Etot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def KEbinaryplanet3(xvec1,xvec2,xvec3,masses):\n",
    "    vel0=xvec1[3:]\n",
    "    vel1=xvec2[3:]\n",
    "    vel2=xvec3[3:]\n",
    "    pos0=xvec1[:3]\n",
    "    pos1=xvec2[:3]\n",
    "    pos2=xvec3[:3]\n",
    "    KE0 = 1/2 * masses[0]*np.sum(np.square(vel0))\n",
    "    KE1 = 1/2*masses[1]*np.sum(np.square(vel1))\n",
    "    KE2=1/2*masses[2]*np.sum(np.square(vel2))\n",
    "    KEplanet = KE2\n",
    "    KEbinary = KE0+KE1\n",
    "    #print(pos0)\n",
    "    ##print(pos1)\n",
    "    r01=np.sqrt(np.sum(np.square(pos0-pos1)))\n",
    "    #print(r01)\n",
    "    \n",
    "    PE01 = - masses[0]*masses[1]/np.sqrt(np.sum(np.square(pos0-pos1)))\n",
    "    PE02 = -masses[0]*masses[2]/np.sqrt(np.sum(np.square(pos0-pos2)))\n",
    "    PE12= - masses[1]*masses[2]/np.sqrt(np.sum(np.square(pos1-pos2)))\n",
    "    PEplanet=PE12+PE02\n",
    "    PEbinary = PE01\n",
    "    Etotal = PEplanet+PEbinary+KEplanet+KEbinary\n",
    "    return KEbinary, KEplanet, PEbinary, PEplanet, Etotal\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "xvec1=xvec[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "xvec2=xvec[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "xvec3=xvec[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-0.        ,  0.44609974,  0.        ]),\n",
       " array([   0.        ,    0.        , 2679.94831257]),\n",
       " -500.1015025879416)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SimpleE3(xvec1,xvec2,xvec3,masses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RunFirstOrbitTwoStarsPlanet(dt,initcondition,printoutput,ODEeq):\n",
    "    masses,xvec0, avec=initcondition\n",
    "    x0=xvec0[2,0]\n",
    "    y0=xvec0[2,1]\n",
    "    lastvalue=[1,1]\n",
    "    thisvalue=[xvec0[2,0],xvec0[2,1]]\n",
    "    xvec=xvec0\n",
    "    i=1\n",
    "    t=0\n",
    "    Earr = []\n",
    "    Larr = []\n",
    "    parr =[]\n",
    "    deltaEarr =[]\n",
    "    deltaLarr =[]\n",
    "    deltaparr =[]\n",
    "    xvec1_0=xvec0[0]\n",
    "    xvec2_0=xvec0[1]\n",
    "    xvec3_0=xvec0[2]\n",
    "    KEplanetarr=[]\n",
    "    PEplanetarr=[]\n",
    "    Eplanetarr=[]\n",
    "    ptot0,angLtot0,Etot0=SimpleE3(xvec1_0,xvec2_0,xvec3_0,masses)\n",
    "    KEbinary, KEplanet, PEbinary, PEplanet, Etotal=KEbinaryplanet3(xvec1_0,xvec2_0,xvec3_0,masses)\n",
    "    while not (lastvalue[0]>0 and lastvalue[1]<0 and thisvalue[0]> 0 and thisvalue[1]>=0):\n",
    "        masses,xvec,avec,t=ODEeq.timestepRK4ODE(i,dt)\n",
    "        if printoutput:\n",
    "            print(xvec[2])\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[2,0],xvec[2,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",
    "        xvec1=xvec[0]\n",
    "        xvec2=xvec[1]\n",
    "        xvec3=xvec[2]\n",
    "        ptot,angLtot,Etot=SimpleE3(xvec1,xvec2,xvec3,masses)\n",
    "        KEbinary, KEplanet, PEbinary, PEplanet, Etotal=KEbinaryplanet3(xvec1,xvec2,xvec3,masses)\n",
    "        deltaE=np.abs((Etot-Etot0)/Etot0)\n",
    "        deltaL=np.abs((angLtot-angLtot0)/np.sqrt(np.sum(np.square(angLtot0))))\n",
    "        deltap=np.abs((ptot-ptot0)/np.sqrt(np.sum(np.square(ptot0))))\n",
    "        Earr.append(Etot)\n",
    "        Larr.append(angLtot)\n",
    "        parr.append(ptot)\n",
    "        deltaEarr.append(deltaE)\n",
    "        deltaLarr.append(deltaL)\n",
    "        deltaparr.append(deltap)\n",
    "        KEplanetarr.append(KEplanet)\n",
    "        PEplanetarr.append(PEplanet)\n",
    "        Eplanetarr.append(KEplanet+PEplanet)\n",
    "        if printoutput:\n",
    "            print(xy,lastvalue,thisvalue)\n",
    "        i+=1\n",
    "    #xvec1=xvec[0]\n",
    "    #xvec2=xvec[1]\n",
    "    #xvec3=xvec[2]\n",
    "    #ptot,angLtot,Etot=SimpleE3(xvec1,xvec2,xvec3,masses)\n",
    "    deltaE=np.abs((Etot-Etot0)/Etot0)\n",
    "    deltaL=np.abs((angLtot-angLtot0)/np.sqrt(np.sum(np.square(angLtot0))))\n",
    "    deltap=np.abs((ptot-ptot0)/np.sqrt(np.sum(np.square(ptot0))))\n",
    "    #print(thisvalue[0], thisvalue[1], lastvalue[0], lastvalue[1], x0, y0)\n",
    "    #print(Etot,Etot0,deltaE)\n",
    "    #print(angLtot,angLtot0,deltaL)\n",
    "    #print(ptot,ptot0,deltap)\n",
    "    return np.abs(x0-(thisvalue[0]+lastvalue[0])/2.),Eplanetarr,KEplanetarr,PEplanetarr,xvec0,xvec,Etot0,Earr,deltaEarr,angLtot0,Larr,deltaLarr,ptot0,parr,deltaparr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run Run 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L0, Eplanetarr,KEplanetarr,PEplanetarr,xvec0,xvec,Etot0,Earr,deltaEarr,angLtot0,Larr,deltaLarr,ptot0,parr,deltaparr=RunFirstOrbitTwoStarsPlanet(dt,xyuvatwostarsplanet,False,ODEeq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'dtarr.csv',dt*np.arange(len(Earr)),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'deltaEarr.csv',np.array(deltaEarr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'deltaLarr.csv',np.array(deltaLarr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'deltaparr.csv',np.array(deltaparr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'E0arr.csv',np.array([Etot0,Etot0]),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Etot0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'L0arr.csv',np.array(angLtot0),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'p0arr.csv',np.array(ptot0),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Eorbarr.csv',np.array(Earr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Lorbarr.csv',np.array(Larr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'porbarr.csv',np.array(parr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+\"Eplanet.csv\",np.array(Eplanetarr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+\"KEplanet.csv\",np.array(KEplanetarr),delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+\"PEplanet.csv\",np.array(PEplanetarr),delimiter=',')"
   ]
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": [
    "#xvec0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Etot0,Earr[len(Earr)-1],deltaEarr[len(Earr)-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angLtot0, Larr[len(Larr)-1], deltaLarr[len(Larr)-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptot0, parr[len(parr)-1], deltaparr[len(parr)-1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##xvec1_0=xvec0[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#xvec1_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#xvec2_0=xvec0[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#xvec3_0=xvec0[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#xvec1=xvec[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#xvec2=xvec[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(dt*np.arange(len(Earr)),Earr)\n",
    "plt.show()\n",
    "#not showing the right plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(dt*np.arange(len(parr)),np.array(parr)[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(dt*np.arange(len(Larr)),np.array(Larr)[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eplanetarr[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(dt*np.arange(len(Eplanetarr)),np.array(Eplanetarr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lorbarr\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dtarr is times\n",
    "#Eorbarr is only for one specific orbit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eorbarr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ptot0,angLtot0,Etot0=SimpleE3(xvec1_0,xvec2_0,xvec3_0,masses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ptot0,angLtot0,Etot0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ptotf,angLtotf,Etotf=SimpleE3(xvec1,xvec2,xvec3,masses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ptotf,angLtotf,Etotf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(L0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Energy is superficially conserved after one orbit, and so is p and L. deltap/p <7e-12, delta L/L=1e-5, delta E/E=3e-5 dt = .1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convergence vs t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.03125, 0.0625, 0.125, 0.25, 0.5] [0.4629301289722889, 1.8530160871680437, 4.172217434659274, 7.439584629903379, 11.626947088124211]\n"
     ]
    }
   ],
   "source": [
    "dt=1 #*31556926 #seconds per year\n",
    "numsteps=2000\n",
    "#mtotal,mass0,xvec,avec=xyuvaeqreducedmass\n",
    "L0arr=[]\n",
    "dtarr=[]\n",
    "Einitarr=[]\n",
    "Eorbarr=[]\n",
    "deltaEarr=[]\n",
    "Linitarr=[]\n",
    "Lorbarr=[]\n",
    "deltaLarr=[]\n",
    "pinitarr=[]\n",
    "porbarr=[]\n",
    "deltaparr=[]\n",
    "Eplanetarr=[]\n",
    "KEplanetarr=[]\n",
    "PEplanetarr=[]\n",
    "dtinit=0.125/4.\n",
    "for i in np.arange(5):\n",
    "    dt=(2**i)*dtinit\n",
    "    dtarr.append(dt)\n",
    "    #L0, xvec0,xvec,Etot0,Earr,deltaEarr,angLtot0,Larr,deltaLarr,ptot0,parr,deltaparr\n",
    "    L0, Eplanet,KEplanet,PEplanet, xvec0, xvec, Einit, Eorb, deltaE, Linit, Lorb, deltaL,pinit,porb,deltap=RunFirstOrbitTwoStarsPlanet(dt,xyuvatwostarsplanet,False,ODEeq)\n",
    "    #L0,Einit,Eorb,deltaE,Linit,Lorb,deltaL=RunFirstOrbitTwoStarsPlanet(dt,xyuvatwostarsplanet,False,ODEeq)\n",
    "    L0arr.append(L0)\n",
    "    Einitarr.append(Einit)\n",
    "    Eorbarr.append(Eorb)\n",
    "    deltaEarr.append(deltaE)\n",
    "    Linitarr.append(Linit)\n",
    "    Lorbarr.append(Lorb)\n",
    "    deltaLarr.append(deltaL)\n",
    "    pinitarr.append(pinit)\n",
    "    porbarr.append(porb)\n",
    "    deltaparr.append(deltap)\n",
    "    Eplanetarr.append(Eplanet)\n",
    "    KEplanetarr.append(KEplanet)\n",
    "    PEplanetarr.append(PEplanet)\n",
    "print(dtarr, L0arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "216603"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(deltaEarr[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "#L0arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "rundir = \"ConvergencevstTest/\"\n",
    "runoutprefix=\"Run3/\"\n",
    "subprefix = \"dt\"\n",
    "maxdt=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ConvergencevstTest/Run3/dtarr.csv'"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rundir+runoutprefix+'dtarr.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'dtarr.csv',dtarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.4629301289722889,\n",
       " 1.8530160871680437,\n",
       " 4.172217434659274,\n",
       " 7.439584629903379,\n",
       " 11.626947088124211]"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L0arr[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'L0arr.csv',L0arr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'deltaEarr.csv',deltaEarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'deltaLarr.csv',deltaLarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'deltaparr.csv',deltaparr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Eintarr.csv',Einitarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Lintarr.csv',Linitarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'pintarr.csv',pinitarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'Eorbarr.csv',Eorbarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'Lorbarr.csv',Lorbarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'porbarr.csv',porbarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'Eplanet.csv',Eplanetarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'KEplanet.csv',KEplanetarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(maxdt):\n",
    "    np.savetxt(rundir+runoutprefix+subprefix+str(i)+\"/\"+'PEplanet.csv',PEplanetarr[i],delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DONT NEED A FIT SINCE ITS VS TIME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit to E convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, rval, pval, stderr = stats.linregress(np.log10(dtarr),np.log10(deltaEarr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(slope, intercept, rval, pval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.xlabel('Log10(dt)')\n",
    "plt.ylabel('Log10(deltaE/E)')\n",
    "plt.title('Convergence of Relative Error of Energy')\n",
    "ax1,=plt.plot(np.log10(dtarr),np.log10(deltaEarr),'g*', label=\"Data\")\n",
    "ax2,=plt.plot(np.log10(dtarr),slope*np.log10(dtarr)+intercept,'b',label=\"Fit, exp=\"+str(slope))\n",
    "plt.legend([ax1,ax2],['Data','Fit, exp='+str(slope)])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# P convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npdeltap=np.array(pinitarr)-np.array(porbarr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp=npdeltap[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope3, intercept3, rval3, pval3, stderr=stats.linregress(np.log10(dtarr),np.log10(np.abs(dp)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(slope3,intercept3,rval3,pval3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "looks like p may be in roundoff. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npdeltaL=np.array(deltaLarr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npdeltaL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dL=npdeltaL[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope4, intercept4, rval4, pval4, stderr=stats.linregress(np.log10(dtarr),np.log10(np.abs(dL)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope4, intercept4, rval4, pval4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.xlabel('Log10(dt)')\n",
    "plt.ylabel('Log10(delta L/L)')\n",
    "plt.title('Convergence of Relative error of Angular Momentum')\n",
    "ax1,=plt.plot(np.log10(dtarr),np.log10(dL),'g*', label=\"Data\")\n",
    "ax2,=plt.plot(np.log10(dtarr),slope4*np.log10(dtarr)+intercept4,'b',label=\"Fit, exp=\"+str(slope4))\n",
    "plt.legend([ax1,ax2],['Data','Fit, exp='+str(slope4)])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit to L0 convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope2, intercept2, rval2, pval2, stderr=stats.linregress(np.log10(dtarr),np.log10(L0arr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(slope2,intercept2, rval2, pval2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitparam=np.polyfit(np.log2(dtarr[0:5]),np.log2(L0arr),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": [
    "xvals=2*np.arange(1,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yvals = xvals**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyfitparam=np.polyfit(np.log2(xvals),np.log2(yvals),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyfitparam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(np.log10(xvals),np.log10(yvals),'*')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errfn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.log2(dtarr[0:5]),np.log2(L0arr),'*', label=\"L0 error\")\n",
    "plt.xlabel('Log2(dt)')\n",
    "plt.ylabel('Log2(L0)')\n",
    "plt.title('Convergence of L0')\n",
    "plt.plot(np.log2(dtarr[0:5]),errfn(np.log2(dtarr[0:5])),'g', label=\"Fit, dt^.78\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Slope is 0.8 rather than 4. It looks like I may have interchanged the slope and the intercept of the log-log plot in my circular orbits convergence test as well. Probably there is something wrong with my RK4 routine underlying all of this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "nptimes, npstar1x, npstar2x,npstar3x, npstar1a, npstar2a,npstar3a= timestepNonAdapt3star.timestepNonAdapt3star(numsteps,dt,xyuvatwostarsplanet,outputevery)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar1x[0:10,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rundir='NonAdaptiveData/'\n",
    "#runinprefix='Run3/twostarsoneplanet'\n",
    "rundir='ConvergenceTestwithKEPE/'\n",
    "runinprefix='Run3/twostarsoneplanet'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'dtarr.csv', newline='') as csvfile:\n",
    "    dtarrlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(dtarrlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "dtarr=np.zeros(len(dtarrlist))\n",
    "count =0\n",
    "for row in dtarrlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        dtarr[count]=number\n",
    "    count+=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'L0arr.csv', newline='') as csvfile:\n",
    "    L0arrlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(L0arrlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "L0arr=np.zeros(len(L0arrlist))\n",
    "count =0\n",
    "for row in L0arrlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        L0arr[count]=number\n",
    "    count+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L0arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "paramsdict={}\n",
    "import json \n",
    "with open(rundir+runinprefix+'params.json') as jsonfile:\n",
    "    paramsdict = json.load(jsonfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "paramsdict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rad0 = paramsdict['rad0']\n",
    "ecc0 = paramsdict['ecc0']\n",
    "theta0=paramsdict['theta0']\n",
    "masses = paramsdict['masses']\n",
    "numsteps=paramsdict['numsteps']\n",
    "outputevery=paramsdict['outputevery']\n",
    "dt = paramsdict['dt']\n",
    "initdatinnerlist=paramsdict['initdatinner']\n",
    "initdatouterlist=paramsdict['initdatouter']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initdatinner =[]\n",
    "for item in initdatinnerlist:\n",
    "    if np.isscalar(item):\n",
    "        initdatinner.append(item)\n",
    "    else:\n",
    "        initdatinner.append(np.array(item))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initdatinner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initdatouter =[]\n",
    "for item in initdatouterlist:\n",
    "    if np.isscalar(item):\n",
    "        initdatouter.append(item)\n",
    "    else:\n",
    "        initdatouter.append(np.array(item))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initdatouter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "useperihelion=True\n",
    "import getxyuvtwouneq_apper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyuvouter=getxyuvtwouneq_apper.getxyuvequneq(initdatouter,useperihelion)\n",
    "print(xyuvouter)\n",
    "massesouter,xvecouter,avecouter=xyuvouter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyuvinner=getxyuvtwouneq_apper.getxyuvequneq(initdatinner,useperihelion)\n",
    "print(xyuvinner)\n",
    "massesinner,xvecinner,avecinner=xyuvinner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'deltaEarr.csv', newline='') as csvfile:\n",
    "    deltaElist = list(csv.reader(csvfile))\n",
    "\n",
    "print(deltaElist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaEarr=np.zeros([len(deltaElist),len(deltaElist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in deltaElist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        deltaEarr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(deltaEarr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'deltaLarr.csv', newline='') as csvfile:\n",
    "    deltaLlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(deltaLlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaLarr=np.zeros([len(deltaLlist),len(deltaLlist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in deltaLlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        deltaLarr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(deltaLarr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'deltaparr.csv', newline='') as csvfile:\n",
    "    deltaplist = list(csv.reader(csvfile))\n",
    "\n",
    "print(deltaplist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaparr=np.zeros([len(deltaplist),len(deltaplist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in deltaplist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        deltaparr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(deltaparr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'E0arr.csv', newline='') as csvfile:\n",
    "    Eintlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(Eintlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0arr=np.zeros([len(Eintlist),len(Eintlist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in Eintlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        E0arr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(E0arr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'L0arr.csv', newline='') as csvfile:\n",
    "    Lintlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(Lintlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L0arr=np.zeros([len(Lintlist),len(Lintlist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in Lintlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        L0arr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(L0arr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'p0arr.csv', newline='') as csvfile:\n",
    "    pintlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(pintlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0arr=np.zeros([len(pintlist),len(pintlist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in pintlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        p0arr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(p0arr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'Eorbarr.csv', newline='') as csvfile:\n",
    "    Eorblist = list(csv.reader(csvfile))\n",
    "\n",
    "print(Eorblist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eorbarr=np.zeros([len(Eorblist),len(Eorblist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in Eorblist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        Eorbarr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Eorbarr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'Lorbarr.csv', newline='') as csvfile:\n",
    "    Lorblist = list(csv.reader(csvfile))\n",
    "\n",
    "print(Lorblist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lorbarr=np.zeros([len(Lorblist),len(Lorblist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in Lorblist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        Lorbarr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Lorbarr[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'porbarr.csv', newline='') as csvfile:\n",
    "    porblist = list(csv.reader(csvfile))\n",
    "\n",
    "print(porblist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "porbarr=np.zeros([len(porblist),len(porblist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in porblist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        porbarr[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(porbarr[0:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Orbit Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import ellipsefns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "t= np.linspace(-math.pi,math.pi,100000)\n",
    "xax=-200+np.arange(0,400)\n",
    "plt.ylim(-5000,5000)\n",
    "plt.xlim(-5000,5000)\n",
    "#plt.axis('equal')\n",
    "mtot=masses[0]+masses[1]\n",
    "semimajor0=rad0[0]*masses[1]/mtot\n",
    "semimajor1=rad0[0]*masses[0]/mtot\n",
    "#plt.title(\"Adaptive RK4, Rp=\"+str(rad0)+\", m1=\"+str(masses[0])+\", m2=\"+str(masses[1])+\" e=\"+str(ecc0)+\", theta=\"+str(round(180*theta0/math.pi)))\n",
    "plt.title(\"Two massive stars and a planet, circular orbits\")\n",
    "plt.plot(ellipsefns.getxtheory(theta0[0],ecc0[0],semimajor0,t,1),ellipsefns.getytheory(theta0[0],ecc0[0],semimajor0,t,1),'b')\n",
    "plt.plot(ellipsefns.getxtheory(theta0[0],ecc0[0],semimajor1,t,-1),ellipsefns.getytheory(theta0[0],ecc0[0],semimajor1,t,-1),'g')\n",
    "plt.plot(ellipsefns.getxtheory(theta0[1],ecc0[1],rad0[1],t,1),ellipsefns.getytheory(theta0[1],ecc0[1],rad0[1],t,1),'c')\n",
    "plt.plot(npstar1x[:,0],npstar1x[:,1],'y')\n",
    "plt.plot(npstar2x[:,0],npstar2x[:,1],'r')\n",
    "plt.plot(npstar3x[:,0],npstar3x[:,1],'k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "t= np.linspace(-math.pi,math.pi,100000)\n",
    "xax=-200+np.arange(0,400)\n",
    "plt.ylim(-50,50)\n",
    "plt.xlim(-50,50)\n",
    "#plt.axis('equal')\n",
    "mtot=masses[0]+masses[1]\n",
    "semimajor0=rad0[0]*masses[1]/mtot\n",
    "semimajor1=rad0[0]*masses[0]/mtot\n",
    "print(masses)\n",
    "print(mtot)\n",
    "print(rad0)\n",
    "#circular\n",
    "#plt.title(\"Adaptive RK4, Rp=\"+str(rad0)+\", m1=\"+str(masses[0])+\", m2=\"+str(masses[1])+\" e=\"+str(ecc0)+\", theta=\"+str(round(180*theta0/math.pi)))\n",
    "plt.title(\"Two massive stars and a planet, circular orbits\")\n",
    "plt.plot(ellipsefns.getxtheory(theta0[0],ecc0[0],semimajor0,t,1),ellipsefns.getytheory(theta0[0],ecc0[0],semimajor0,t,1),'b')\n",
    "plt.plot(ellipsefns.getxtheory(theta0[0],ecc0[0],semimajor1,t,-1),ellipsefns.getytheory(theta0[0],ecc0[0],semimajor1,t,-1),'g')\n",
    "plt.plot(ellipsefns.getxtheory(theta0[1],ecc0[1],rad0[1],t,1),ellipsefns.getytheory(theta0[1],ecc0[1],rad0[1],t,1),'c')\n",
    "plt.plot(npstar1x[:,0],npstar1x[:,1],'y')\n",
    "plt.plot(npstar2x[:,0],npstar2x[:,1],'r')\n",
    "plt.plot(npstar3x[:,0],npstar3x[:,1],'k')\n",
    "plt.plot(ellipsefns.getxtheory(theta0[0],ecc0[0],semimajor0,t,1),ellipsefns.getytheory(theta0[0],ecc0[0],semimajor0,t,1),'b')\n",
    "plt.plot(ellipsefns.getxtheory(theta0[0],ecc0[0],semimajor1,t,-1),ellipsefns.getytheory(theta0[0],ecc0[0],semimajor1,t,-1),'g')\n",
    "plt.plot(ellipsefns.getxtheory(theta0[1],ecc0[1],rad0[1],t,1),ellipsefns.getytheory(theta0[1],ecc0[1],rad0[1],t,1),'c')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Precession of the central two stars about the three object center of mass and transfer of binding energy from the central two stars to the outer planet as the outer planet moves inward and the radii of the inner stars orbits gets larger (less bound)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eqn 15 pg 95 of \n",
    "An Introduction to the Classical Three-Body Problem\n",
    "From Periodic Solutions to Instabilities and Chaos\n",
    "Govind S Krishnaswami and Himalaya Senapati\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Restricted Three-Body Problem\n",
    "The restricted three-body problem is a simplified version of the\n",
    "three-body problem where one of the masses m3 is assumed much\n",
    "smaller than the primaries m1 and m2. The Sun–Earth–Moon\n",
    "system provides an\n",
    "example of the restricted\n",
    "three-body problem\n",
    "where we further have\n",
    "m2 = mE \n",
    " m1 = mS.\n",
    "Thus, m1 and m2 move\n",
    "in Keplerian orbits which are not affected by m3. In the planar circular restricted three-body problem, the primaries move in\n",
    "fixed circular orbits around their common CM with angular speed\n",
    "Ω = (G(m1+m2)/d3)\n",
    "1/2 given by Kepler’s third law and m3 moves\n",
    "in the same plane as m1 and m2. Here, d is the separation between the primaries. This system has 2 degrees of freedom associated to the planar motion of m3, and therefore, a 4-dimensional\n",
    "phase space just like the planar Kepler problem for the reduced\n",
    "mass. However, unlike the latter which has three conserved quantities (energy, z-component of angular momentum and direction\n",
    "of LRL vector) and is exactly solvable, the planar restricted threebody problem has only one known conserved quantity, the ‘Jacobi integral’, which is the energy of m3 in the co-rotating (noninertial) frame of the primaries: The ‘Roche’ effective\n",
    "potential Veff, named\n",
    "after the French\n",
    "astronomer Edouard ´\n",
    "Albert Roche, is a sum\n",
    "of centrifugal and\n",
    "gravitational energies\n",
    "due to m1 and m2.\n",
    "E =[1/2 m3 rdot^2  +1/2 m3 r^2 φdot^2]−1/2 m3 Ω^2 r^2 −G m3 [m1 /r1 + m2/r2]≡ T+Veff. (15)\n",
    "Here, (r, φ) are the plane polar coordinates of m3 in the co-rotating\n",
    "frame of the primaries with origin located at their CM while r1\n",
    "and r2 are the distances of m3 from m1 and m2 (see Figure 4)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RK4 Quality Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "t = np.linspace(0, 2*math.pi, 100)\n",
    "xax=np.arange(0,100*len(npstar2x[:,1]))\n",
    "xax=xax-len(xax)/2\n",
    "oneellipse=int(len(nptimes))\n",
    "plt.ylim(-200,200)\n",
    "plt.xlim(-200,200)\n",
    "#plt.ylim(-1,1)\n",
    "#plt.xlim(74.999,75)\n",
    "plt.title(\"Elliptical Orbit, RK4, a=\"+str(rad0)+\", e=\"+ str(ecc0)+\", theta=\"+str(round(180*theta0/math.pi))+ \" degrees\")\n",
    "plt.plot(ellipsefns.getxtheory(theta0,ecc0,rad0,t,1),ellipsefns.getytheory(theta0,ecc0,rad0,t,1),'b')\n",
    "plt.plot(ellipsefns.getxtheory(theta0,ecc0,rad0,t,-1),ellipsefns.getytheory(theta0,ecc0,rad0,t,-1),'g')\n",
    "plt.plot(npstar1x[:oneellipse,0],npstar1x[:oneellipse,1],'y')\n",
    "plt.plot(npstar2x[:oneellipse,0],npstar2x[:oneellipse,1],'r')\n",
    "plt.plot(xax, np.tan(theta0)*xax,'c')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "t = np.linspace(0, 2*math.pi, 100)\n",
    "xax=np.arange(0,100*len(npstar2x[:,1]))\n",
    "xax=xax-len(xax)/2\n",
    "\n",
    "#plt.ylim(-200,200)\n",
    "#plt.xlim(-200,200)\n",
    "plt.ylim(-1,1)\n",
    "plt.xlim(74.999,75)\n",
    "plt.title(\"Elliptical Orbit, RK4, a=\"+str(rad0)+\", e=\"+ str(ecc0)+\", theta=\"+str(round(180*theta0/math.pi))+ \" degrees\")\n",
    "plt.plot(ellipsefns.getxtheory(theta0,ecc0,rad0,t,1),ellipsefns.getytheory(theta0,ecc0,rad0,t,1),'b')\n",
    "plt.plot(ellipsefns.getxtheory(theta0,ecc0,rad0,t,-1),ellipsefns.getytheory(theta0,ecc0,rad0,t,-1),'g')\n",
    "plt.plot(npstar1x[:oneellipse,0],npstar1x[:oneellipse,1],'*y')\n",
    "plt.plot(npstar2x[:oneellipse,0],npstar2x[:oneellipse,1],'*r')\n",
    "plt.plot(xax, np.tan(theta0)*xax,'c')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "t = np.linspace(0, 2*math.pi, 100)\n",
    "xax=np.arange(0,100*len(npstar2x[:,1]))\n",
    "xax=xax-len(xax)/2\n",
    "\n",
    "#plt.ylim(-200,200)\n",
    "#plt.xlim(-200,200)\n",
    "plt.ylim(-10,10)\n",
    "plt.xlim(-25,-20)\n",
    "plt.title(\"Elliptical Orbit, RK4, a=\"+str(rad0)+\", e=\"+ str(ecc0)+\", theta=\"+str(round(180*theta0/math.pi))+ \" degrees\")\n",
    "plt.plot(ellipsefns.getxtheory(theta0,ecc0,rad0,t,1),ellipsefns.getytheory(theta0,ecc0,rad0,t,1),'b')\n",
    "plt.plot(ellipsefns.getxtheory(theta0,ecc0,rad0,t,-1),ellipsefns.getytheory(theta0,ecc0,rad0,t,-1),'g')\n",
    "plt.plot(npstar1x[:,0],npstar1x[:,1],'*y')\n",
    "plt.plot(npstar2x[:,0],npstar2x[:,1],'*r')\n",
    "plt.plot(xax, np.tan(theta0)*xax,'c')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oneellipse=len(nptimes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "plt.xlabel(\"timestep\")\n",
    "plt.ylabel(\"Rel err: 1/30 *|x2-x1|*(h'/h)^5\")\n",
    "plt.title(\"Relative error, RK4\")\n",
    "relerrplt,=plt.plot(np.arange(2,oneellipse),nperr[2:oneellipse],'b',label=\" dtmin = \"+ str(dtmin)+\" delta = \"+str(delta))\n",
    "plt.legend(handles=[relerrplt])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This curve should be constant and about equal to delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "plt.xlabel(\"timestep number\")\n",
    "plt.ylabel(\"timestep size\")\n",
    "plt.title(\"Evolution of timestep with step\")\n",
    "relerrplt,=plt.plot(np.arange(2,oneellipse),npdt[2:oneellipse],'b',label=\" dtmin = \"+ str(dtmin)+\" delta = \"+str(delta))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For appropriately set tolerances, adaptive timestep no longer has feedback issues. Tol=delta=0.0005 or less. The important thing is that this curve is smooth and has no constant regions where the value is pinned to either the maximum or minimum. We also seek reasonable starting values for the timestep. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"relative error\")\n",
    "plt.title(\"horizontal position vs relative error\")\n",
    "plt.plot(npstar1x[2:oneellipse,0], nperr[2:oneellipse],'b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"y\")\n",
    "plt.ylabel(\"rel error\")\n",
    "plt.title(\"y vs rel error\")\n",
    "plt.plot(npstar1x[2:oneellipse,1],nperr[2:oneellipse])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"r\")\n",
    "plt.ylabel(\"rel error\")\n",
    "plt.title(\"r vs rel error\")\n",
    "rstar1x=np.sqrt(npstar1x[:,0]**2+npstar1x[:,1]**2)\n",
    "plt.plot(rstar1x[2:oneellipse],nperr[2:oneellipse],'*b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"r\")\n",
    "plt.title(\"r vs t\")\n",
    "plt.plot(nptimes[2:oneellipse], rstar1x[2:oneellipse], '*b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thetaellipse = np.arctan2(npstar1x[2:oneellipse,1],npstar1x[2:oneellipse,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"theta\")\n",
    "plt.title(\"theta vs t\")\n",
    "plt.plot(nptimes[2:oneellipse],thetaellipse,'b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"theta\")\n",
    "plt.ylabel(\"rel error\")\n",
    "plt.title(\"rel error vs theta\")\n",
    "plt.plot(thetaellipse, nperr[2:oneellipse],'*b')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conservation of angular momentum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Lstability3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaL, Lavg, npstarLtot,npcm,npvcm=Lstability3.Lstability3(npstar1x,npstar2x,npstar3x,masses)\n",
    "deltaL, Lavg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaL*Lavg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel('time')\n",
    "plt.ylabel('Center of Mass X coordinate')\n",
    "plt.title('M1*x1+M2*x2+M3*x3')\n",
    "plt.plot(nptimes,npcm[:,0],'*b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(npcm[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(npcm[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel('time')\n",
    "plt.ylabel('Center of Mass Y coordinate')\n",
    "plt.title('M1*y1+M2*y2+M3*y3')\n",
    "plt.plot(nptimes,npcm[:,1],'*b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel('time')\n",
    "plt.ylabel('Center of Mass Z coordinate')\n",
    "plt.title('M1*z1+M2*z2+M3*z3')\n",
    "plt.plot(nptimes,npcm[:,2],'*b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel('time')\n",
    "plt.ylabel('Total L')\n",
    "plt.title('L1+L2+L3')\n",
    "plt.plot(nptimes,npstarLtot[:,2],'*b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npstarLtot[0:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total L should not be conserved (See Euler problem on wikipedia). It is not. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conservation of Energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Estability3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaE, Eavg, Egrav, Ekinetic, Etot,Ecmkinetic, Eplanetgexact, Eplanetgapprox,Eplanetk= Estability3.Estability3(npstar1x, npstar2x, npstar3x, masses)\n",
    "deltaE, Eavg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaE*Eavg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,Egrav)\n",
    "plt.title('Potential energy')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('U1+U2+U3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,Ekinetic)\n",
    "plt.title('Kinetic Energy')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('K1+K2+K3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,Etot)\n",
    "plt.title('Total Energy')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('E1+E2+E3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,Egrav+2*Ekinetic,'b')\n",
    "plt.title('Virial (binding) energy')\n",
    "plt.ylabel('2*Ekinetic-Egrav')\n",
    "plt.xlabel('time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Eplanetgexact+Eplanetk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Eplanetgapprox+Eplanetk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Eplanetgexact+2*Eplanetk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Ekinetic-Eplanetk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Egrav-Eplanetgexact)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ApproxBindingPlanet=Eplanetgexact+2*Eplanetk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ExactBindingPlanet=Eplanetgexact+2*Eplanetk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BindingSystem=Ekinetic+2*Egrav"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,BindingSystem-ApproxBindingPlanet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,BindingSystem-ExactBindingPlanet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,Etot-0.5*Egrav)\n",
    "plt.title('Residual of binding energy from total energy')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('Etot+V/2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binding energy for a circular orbit should be constant around the orbit. It appears the orbit is not circular-- it is eliptical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,Ecmkinetic)\n",
    "plt.title('Kinetic energy of the center of mass')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Ecm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total E should be conserved. See Euler problem on wikipedia. It is not. This is a problem. Perhaps there is a problem in the numerical integration routine still. Total energy is linearly decreasing, binding energy is oscillating (2K+V). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conservation of linear momentum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pstability3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltap, pavg, pmag,ptot = pstability3.pstability3(npstar1x, npstar2x, npstar3x, masses)\n",
    "deltap, pavg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deltap*pavg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,ptot[:,0])\n",
    "plt.title('Linear momentum x')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('m1*v1x+m2*v2x+m3*v3x')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,ptot[:,1])\n",
    "plt.title('Linear momentum y')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('m1*v1y+m2*v2y+m3*v3y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,ptot[:,2])\n",
    "plt.title('Linear momentum z')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('m1*v1z+m2*v2z+m3*v3z')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes,pmag)\n",
    "plt.title('Magnitude of linear momentum')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('|p|')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conservation of Jacobi Integral for 3 body"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orbit conserves energy to a relative error of deltaE. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "def JacobiStability2(npstar1x,npstar2x,npstar3x,masses,rad0):\n",
    "    npstar1pos=npstar1x[:,0:3]\n",
    "    npstar1v=npstar1x[:,3:]\n",
    "    npstar2pos=npstar2x[:,0:3]\n",
    "    npstar2v=npstar2x[:,3:]\n",
    "    npstar3pos=npstar3x[:,0:3]\n",
    "    npstar3v=npstar3x[:,3:]\n",
    "    cmcoords=(masses[0]*npstar1pos+masses[1]*npstar2pos+masses[2]*npstar3pos)/(masses[0]+masses[1]+masses[2])\n",
    "    \n",
    "    cmprimarycoords=(masses[0]*npstar1pos+masses[1]*npstar2pos)/(masses[0]+masses[1])\n",
    "    #print(cmcoords)\n",
    "    cmprimaryv=(masses[0]*npstar1v+masses[1]*npstar2v)/(masses[0]+masses[1])\n",
    "    r3cm=npstar3pos-cmprimarycoords\n",
    "    v3cm=npstar3v-cmprimaryv\n",
    "    r12=rad0[0] #np.sqrt(np.sum((npstar1pos-npstar2pos)**2,axis=1))\n",
    "    r23=np.sqrt(np.sum((npstar2pos-npstar3pos)**2,axis=1))\n",
    "    r31=np.sqrt(np.sum((npstar3pos-npstar1pos)**2,axis=1))\n",
    "    print(r12,r23,r31)\n",
    "    Omega=((masses[0]+masses[1])/r12**3)**0.5 #angular velocity of inner two stars  \n",
    "    rotmat = np.array([[np.cos(Omega),np.sin(Omega)],[-np.sin(Omega),np.cos(Omega)]])\n",
    "    print(rotmat)\n",
    "    v3shape=v3cm.shape\n",
    "    print(v3shape)\n",
    "    v3cmrot=np.zeros([v3shape[0],2])\n",
    "    for row in np.arange(v3shape[0]):\n",
    "        for i in np.arange(0,2):\n",
    "            for j in np.arange(0,2):\n",
    "                v3cmrot[row,i] = rotmat[i,j]*v3cm[row,j]\n",
    "    print(\"Omegasq=\",Omega**2)\n",
    "    Egrav=-masses[2]*(masses[0]/r31+masses[1]/r23)\n",
    "    Ekinetic=masses[2]*0.5*np.sum(v3cmrot**2,axis=1)\n",
    "    r3cmmagsq=np.sum(r3cm**2,axis=1)\n",
    "    print(\"r3\",r3cmmagsq)\n",
    "    Erot=-0.5*masses[2]*Omega**2*r3cmmagsq\n",
    "    Etot=Ekinetic+Egrav+Erot\n",
    "    #print(Egrav,Ekinetic,Erot,Etot)\n",
    "    deltaE=np.abs(np.std(Etot)/np.mean(Etot))\n",
    "    return deltaE, np.mean(Etot), cmcoords, cmprimarycoords, Ekinetic, Egrav, Erot, Etot\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import JacobiStability\n",
    "deltaE, Javg,cm, cmprimary, Ekinetic, Egrav, Erot, Etot = JacobiStability2(npstar1x,npstar2x,npstar3x,masses,rad0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(deltaE, Javg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line,=plt.plot(nptimes,Egrav,'.y',label=\"Egrav\")\n",
    "plt.plot(nptimes,Ekinetic,'c',label=\"Ekinetic\")\n",
    "#plt.plot(nptimes,Egrav+Ekinetic,'k')\n",
    "plt.plot(nptimes,Erot,'.b',label=\"Erot\")\n",
    "plt.plot(nptimes,Etot,'r',label=\"Etot\")\n",
    "plt.title(\"Jacobi integral energy components\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"Energy\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(cmprimary[:,0]-cm[:,0],cmprimary[:,1]-cm[:,1])\n",
    "plt.title(\"Motion of the center of mass of the binary star system\")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(cm[:,0], cm[:,1])\n",
    "plt.title(\"Motion of the center of the mass of the three bodies\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symmetry of the two stars orbits with respect to eachother "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import orbitSymmetricTol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allcancelsx, allcancelsa=orbitSymmetricTol.howSymmetric(npstar1x,npstar2x,npstar1a,npstar2a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allcancelsx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orbit is symmetric with respect to x and v up to a tolerance of allcancelsx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allcancelsa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar1x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orbit is symmetric with respect to a to up to or better than a tolerance of allcancelsa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing for precision in perihelion and aphelion recovery from knowledge of the initial angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import aeTest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "plt.plot( np.abs(npstar1x[:,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findPerapHelionThetaIsZero(npstar1x,ecc0,rad0,useperihelion):\n",
    "    maxindex=len(npstar1x[:,1])\n",
    "    apindex=np.argmin(np.abs(npstar1x[10:int(maxindex/2),1]))+1\n",
    "    print(apindex)\n",
    "    delta=math.floor(.2*apindex)\n",
    "    print(delta)\n",
    "    perindex=np.argmin(np.abs(npstar1x[delta:apindex-delta,1]))+delta\n",
    "    print(perindex)\n",
    "    print(apindex/perindex)\n",
    "    if np.abs(apindex/perindex-2.)>0.05: \n",
    "        perindex=np.argmin(np.abs(npstar1x[apindex+delta:2*apindex+delta,1]))+delta+apindex\n",
    "        print(perindex)\n",
    "        apindex,perindex=perindex,apindex\n",
    "        print(apindex,perindex)\n",
    "        print(apindex/perindex)\n",
    "    #apindex3=np.argmin(np.abs(npstar1x[100:apindex2-100,1]))+100\n",
    "    #apindex4=np.argmin(np.abs(npstar1x[100:apindex3-100,1]))+100\n",
    "    #perindex=np.argmin(np.abs(npstar1x[100:apindex2-100,1]))+100\n",
    "    if useperihelion:\n",
    "        apindex,perindex=perindex,apindex\n",
    "    coordap=[npstar1x[apindex,0],npstar1x[apindex,1]]\n",
    "    coordper=[npstar1x[perindex,0],npstar1x[perindex,1]]\n",
    "    ra=np.sqrt(np.sum(np.array(coordap)**2))\n",
    "    rp=np.sqrt(np.sum(np.array(coordper)**2))\n",
    "    e=(ra-rp)/(ra+rp)\n",
    "    a=1./2.*(ra+rp)\n",
    "    deltae=np.abs(e-ecc0)/ecc0\n",
    "    deltaa=np.abs(a-rad0)/rad0\n",
    "    return apindex,perindex,coordap,coordper,ra,rp,e,a,deltae,deltaa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if theta0==0:\n",
    "    apindex,perindex,coordap,coordper,ra,rp,e,a,deltae,deltaa=findPerapHelionThetaIsZero(npstar1x,ecc0,rad0,useperihelion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if theta0!=0:\n",
    "    perindex,apindex,coordper,coordap,rp,ra,e,a,deltae,deltaa=aeTest.findPerApHelion(npstar1x,theta0,ecc0,rad0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recovered perihelion radius given by r_per, recovered aphelion radius given by r_ap, a is semimajor axis computed from the orbit, e is eccentricty computed from the orbit, delta a is relative error on the semimajor axis, and delta e is eccenticity relative error. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apindex,perindex,coordap,coordper,ra,rp,e,a,deltae,deltaa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimes[:perindex],npstar1x[:perindex,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit t as a function of the eccentric anomoly (angle from the center of the elipse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findEccentricAnomoly(coordap, coordper):\n",
    "    Ecenter=np.array(coordap)+np.array(coordper)\n",
    "    Ecenter/=2.\n",
    "    return Ecenter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ellipsecenter=findEccentricAnomoly(coordap,coordper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ellipsecenter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getEllipsecenter(ecc,rad,theta, useperihelion):\n",
    "    #ind of use perihelion\n",
    "    rp=rad*(1-ecc)\n",
    "    ra=-rad*(1+ecc)\n",
    "    rcent=(rp+ra)/2\n",
    "    #print(rp,ra,rcent)\n",
    "    rcentx=rcent*np.cos(theta)\n",
    "    rcenty=rcent* np.sin(theta)\n",
    "    return [rcentx,rcenty]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(np.array(coordap)+np.array(coordper))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "getEllipsecenter(ecc0,rad0,theta0, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arctan2(-1.,0)/math.pi*180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getEang(npstar2x,ecc,rad,theta, useperihelion):\n",
    "    ellipsecenter=getEllipsecenter(ecc,rad,theta, useperihelion)\n",
    "    #print(ellipsecenter)\n",
    "    #print(npstar1x)\n",
    "    #print(npstar1x[:,0])\n",
    "    xfortan=npstar2x[:,0]-ellipsecenter[0]\n",
    "    yfortan=npstar2x[:,1]-ellipsecenter[1]\n",
    "    #print(\"xy\",xfortan,yfortan)\n",
    "    Eang=np.arctan2(yfortan,xfortan)\n",
    "    Eangtheta=2*np.arctan2(np.sqrt(1-ecc)*np.sin(theta/2.),np.sqrt(1+ecc)*np.cos(theta/2.))\n",
    "    Eang-=Eangtheta\n",
    "    #if useperihelion:\n",
    "    #    Eang-=math.pi\n",
    "    print(len(Eang))\n",
    "    for i in np.arange(len(Eang)):\n",
    "        #if xfortan[i]<0:\n",
    "        #    Eang[i]+=math.pi\n",
    "        if Eang[i]<0:\n",
    "            Eang[i]+=2*math.pi\n",
    "    print(xfortan,yfortan,Eang/math.pi*180)\n",
    "    print(\"Eang\",Eang)\n",
    "    return xfortan, yfortan,Eang"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arctan2(-6.69327422*10**-6, -0.02240526)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def xofeccanom(rad,Eang,ecc):\n",
    "    return rad*(np.cos(Eang)+ecc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def yofeccanom(rad, Eang, ecc):\n",
    "    b=getsemiminor(rad,ecc)\n",
    "    return b*np.sin(Eang) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tofEccAnom(npstar2x, rad,ecc,theta,mass, useperihelion):\n",
    "    xE0,yE0,Eang=getEang(npstar2x,ecc,rad,theta, useperihelion)\n",
    "    b=getsemiminor(rad,ecc)\n",
    "    #fitval=rad*b*(Eang-ecc*np.sin(Eang))\n",
    "    sgn=-1\n",
    "    if useperihelion:\n",
    "        sgn*=-1\n",
    "    # added a factor of 2 in front of fitval\n",
    "    fitval=rad**(1.5)*mass**(-.5)*(Eang-sgn*ecc*np.sin(Eang))\n",
    "    fitval*=2\n",
    "    print(Eang, ecc,rad,mass,fitval)\n",
    "    return fitval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def residualtEang(pars,nptimes,npstar2x):\n",
    "    residual=np.abs(nptimes-tofEccAnom(npstar2x,pars['rad'],pars['ecc'],pars['theta'],pars['mass'],pars['useperihelion']))\n",
    "    #print(nptimes, residual)\n",
    "    return residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endindex=apindex\n",
    "if useperihelion:\n",
    "    endindex=perindex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endindex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xE0,yE0,Eang0=getEang(npstar2x[0:endindex-1,:],ecc0,rad0,theta0, useperihelion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nptimesoneorbit=nptimes[0:endindex-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Eang0, xE0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Eang0, yE0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(xE0,yE0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,xE0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,yE0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,yE0/xE0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arctan(0./1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yE0[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xE0[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,np.arctan2(yE0,xE0)*180/math.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eangtest=tofEccAnom(npstar2x[0:endindex-1,:],rad0,ecc0,theta0,masses[0], useperihelion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(Eangtest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nptimes[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nptimesoneorbit=nptimes[0:endindex-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tpredoneorbit=Eangtest[0:endindex-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eangoneorbit=Eang0[0:endindex-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,npstar2x[0:endindex-1,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,Eangoneorbit/math.pi*180)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,tpredoneorbit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(nptimesoneorbit,tpredoneorbit, label='times')\n",
    "plt.xlabel(\"T timesteps\")\n",
    "plt.ylabel(\"T computed from angle\")\n",
    "plt.title(\"Times using start parameters\")\n",
    "plt.show()\n",
    "plt.plot(nptimesoneorbit,tpredoneorbit-nptimesoneorbit, label='residual')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(nptimesoneorbit,tpredoneorbit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lmfit import Minimizer,Parameters\n",
    "pfittEang = Parameters()\n",
    "pfittEang.add(name='theta', value=theta0, min=theta0-math.pi, max=theta0+math.pi, vary=True)\n",
    "pfittEang.add(name='rad', value=rad0, min=0, vary=True)\n",
    "pfittEang.add(name='ecc', value=ecc0, min=0, max=1.0, vary=True)\n",
    "pfittEang.add(name='mass', value=masses[0], vary=False)\n",
    "pfittEang.add(name='useperihelion', value=True, vary=False)\n",
    "\n",
    "minit = Minimizer(residualtEang, pfittEang, fcn_args=(nptimesoneorbit,npstar2x[:endindex-1,:]))\n",
    "\n",
    "outt = minit.leastsq()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "outt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Semimajor axis fit is good to about .1 percent at an eccentricity of 0.5. rad=49.886512 for delta=0.0001. rad=49.8863857 for delta = 0.0005. ecc = 0.51764747 for delta = 0.0001. ecc = 0.51765166 for delta = 0.0005. so both ecc and rad improve with a smaller delta, but not quickly. The program runs very slowly for large delta and requires 500000 timesteps for delta=0.0001. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tfinal=nptimesoneorbit[endindex-2]\n",
    "tfinal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.fft import fft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import blackman\n",
    "wndow=blackman(len(nptimesoneorbit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftresiduals=fft(nptimesoneorbit-tpredoneorbit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xf = np.linspace(0.0, 1.0/(2.0*dt), len(nptimesoneorbit)//2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(xf[1:len(nptimesoneorbit)//2], 2.0/len(nptimesoneorbit) * np.abs(ftresiduals[1:len(nptimesoneorbit)//2]), '-b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ftr=fft(np.sqrt(npstar1x[:endindex-1,0]**2+npstar1x[:endindex-1,1]**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"frequency\")\n",
    "plt.ylabel(\"Fourier transform of t(E)\")\n",
    "plt.title(\"Fourier transform of time computed from eccentric anomaly\")\n",
    "plt.semilogy(xf[1:len(nptimesoneorbit)//2], 2.0/len(nptimesoneorbit) * np.abs(ftr[1:len(nptimesoneorbit)//2]), '-b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the fourier transform of the triangle wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.title(\"Fourier transform of residuals of times\")\n",
    "plt.xlabel(\"frequency\")\n",
    "plt.ylabel(\"FT(t(E))\")\n",
    "plt.semilogy(xf[1:50], 2.0/len(nptimesoneorbit) * np.abs(ftresiduals[1:50]), '-b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"bin number\")\n",
    "plt.ylabel(\"FT(t(E))\")\n",
    "plt.title(\"Fourier transform of residuals of times\")\n",
    "plt.semilogy(np.arange(1,30), 2.0/len(nptimesoneorbit[:29]) * np.abs(ftresiduals[1:30]), '*b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Peak at the second harmonic and the sixth harmonic. The second harmonic seems natural due to the symmetry of the orbit but the sixth harmonic is harder to explain. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implementing the Adaptive RK4 hasn't changed the fit, residuals, or harmonics, sadly. There is still a 2% error in the fit for an eccentricity of 0.5"
   ]
  },
  {
   "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": []
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}