TwoStarsOnePlanet_convergenceTest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GENERATES DATA FOR TwoStarsOnePlanet_Esplit.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 6)\n",
      "[[0. 0. 0.]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0.]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "x = np.zeros((1,6))\n",
    "print(x.shape)\n",
    "print(x[:,0:3])\n",
    "x[:,0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6,)\n",
      "[[0. 0. 0.]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0.]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "x = np.zeros((6,))\n",
    "print(x.shape)\n",
    "x = np.reshape(x, (1,6))\n",
    "print(x[:,0:3])\n",
    "x[:,0:3]"
   ]
  },
  {
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ellipsefns\n",
    "#getfocus, getsemiminor, theoreticalx, theoreticaly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import initdatuneq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rundir='NonAdaptiveData/'\n",
    "runoutprefix='Run4/twostarsoneplanet'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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=.1\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": 8,
   "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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import getxyuvtwouneq_apper #unequal mass version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "useperihelion=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 6)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(xvecouter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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": 14,
   "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": 15,
   "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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([100., 100.,   1.])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "masses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xvec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "avec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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": 20,
   "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": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[3.141592653589793, 0.0], 0.0, 1000, 0.0, [200.0, 1.0]]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initdatouterlist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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": 23,
   "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.1,\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": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "paramsdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyuvatwostarsplanet=masses,xvec,avec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import RK4implicit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "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": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import OrbitDiffEq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "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 0x7fd755b3bdc0>"
      ]
     },
     "execution_count": 29,
     "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": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(xyuvaequneq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timestep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timestepNonAdapt3star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "import parabfns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ellipsefns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hyperbfns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Estability3,Lstability3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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",
    "    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",
    "        if printoutput:\n",
    "            print(xy,lastvalue,thisvalue)\n",
    "        i+=1\n",
    "    xvec1_0=xvec0[0]\n",
    "    xvec2_0=xvec0[1]\n",
    "    xvec3_0=xvec0[2]\n",
    "    xvec1=xvec[0]\n",
    "    xvec2=xvec[1]\n",
    "    xvec3=xvec[2]\n",
    "    ptot0,angLtot0,Etot0=SimpleE3(xvec1_0,xvec2_0,xvec3_0,masses)\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.),xvec0,xvec,Etot0,Etot,deltaE,angLtot0,angLtot,deltaL,ptot0,ptot,deltap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L0, xvec0, xvec, Etot0, Etot, deltaE, angLtot0, angLtot, deltaL,ptot0,ptot,deltap=RunFirstOrbitTwoStarsPlanet(dt,xyuvatwostarsplanet,False,ODEeq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xvec0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Etot0, Etot, deltaE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angLtot0, angLtot, deltaL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptot0, ptot, deltap\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": [
    "xvec3=xvec[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "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",
    "    E0 = np.sum(np.square(p0))/2/masses[0]\n",
    "    E1=np.sum(np.square(p1))/2/masses[1]\n",
    "    E2=np.sum(np.square(p2))/2/masses[2]\n",
    "    ptot=p0+p1+p2\n",
    "    angLtot=angL0+angL1+angL2\n",
    "    Etot=E0+E1+E2\n",
    "    return ptot,angLtot,Etot"
   ]
  },
  {
   "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(L0err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Energy is superficially conserved after one orbit, and so is p and L. deltap/p = 1e-8, delta L/L=2e-6, delta E/E=1e-4 dt = .1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "993.1723721634779 0.009438546692053466 993.1715175187385 -0.004476112797052286 995.0248756218906 0.0\n",
      "500.1145166883877 500.0995024875622 3.002242703871989e-05\n",
      "[   0.            0.         2679.91455938] [   0.            0.         2679.94831257] [0.00000000e+00 0.00000000e+00 1.25947158e-05]\n",
      "[-6.48009771e-12  4.46099737e-01  0.00000000e+00] [-0.          0.44609974  0.        ] [1.45261187e-11 2.45264613e-12 0.00000000e+00]\n",
      "990.8521423476368 0.017799293123606685 990.8495783376584 -0.009964989368821228 995.0248756218906 0.0\n",
      "500.11606487601006 500.0995024875622 3.311818621193772e-05\n",
      "[   0.            0.         2679.91135039] [   0.            0.         2679.94831257] [0.00000000e+00 0.00000000e+00 1.37921244e-05]\n",
      "[-2.85678425e-12  4.46099737e-01  0.00000000e+00] [-0.          0.44609974  0.        ] [6.40391378e-12 8.77564499e-12 0.00000000e+00]\n",
      "987.59483947827 0.027826467518052032 987.5880015323094 -0.027519421935417523 995.0248756218906 0.0\n",
      "500.162377575554 500.0995024875622 0.00012572515605207326\n",
      "[   0.            0.         2679.80857365] [   0.            0.         2679.94831257] [0.00000000e+00 0.00000000e+00 5.21423939e-05]\n",
      "[-2.69662764e-12  4.46099737e-01  0.00000000e+00] [-0.          0.44609974  0.        ] [6.04489853e-12 8.29295527e-12 0.00000000e+00]\n",
      "983.3999836819877 0.095588854829247 983.3828924156597 -0.01463207700454855 995.0248756218906 0.0\n",
      "501.66015437141414 500.0995024875622 0.00312068273631361\n",
      "[   0.            0.         2676.47007654] [   0.            0.         2679.94831257] [0.         0.         0.00129787]\n",
      "[-7.52936602e-12  4.46099737e-01  0.00000000e+00] [-0.          0.44609974  0.        ] [1.68782122e-11 2.00168775e-12 0.00000000e+00]\n",
      "978.2563595796494 0.16463290913980363 978.21533936196 -0.05465005374473758 995.0248756218906 0.0\n",
      "636.2801733558558 500.0995024875622 0.27230715125872484\n",
      "[   0.            0.         2426.55758733] [   0.            0.         2679.94831257] [0.        0.        0.0945506]\n",
      "[-8.41526848e-12  4.46099737e-01  0.00000000e+00] [-0.          0.44609974  0.        ] [1.88640965e-11 1.24758928e-11 0.00000000e+00]\n",
      "[0.03125, 0.0625, 0.125, 0.25, 0.5] [1.852930780782458, 4.174015279243008, 7.43345511660084, 11.633437573066885, 16.789026151085977]\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",
    "dtinit=0.125/4.\n",
    "for i in np.arange(5):\n",
    "    dt=(2**i)*dtinit\n",
    "    dtarr.append(dt)\n",
    "    L0, 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",
    "print(dtarr, L0arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.002242703871989e-05, 3.311818621193772e-05, 0.00012572515605207326, 0.00312068273631361, 0.27230715125872484]\n"
     ]
    }
   ],
   "source": [
    "print(deltaEarr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L0arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "rundir = \"ConvergenceTest/\"\n",
    "runoutprefix=\"Run1/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ConvergenceTest/Run1/dtarr.csv'"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rundir+runoutprefix+'dtarr.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'dtarr.csv',dtarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'L0arr.csv',L0arr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'deltaEarr.csv',deltaEarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'deltaLarr.csv',deltaLarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'deltaparr.csv',deltaparr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Eintarr.csv',Einitarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Lintarr.csv',Linitarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'pintarr.csv',pinitarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Eorbarr.csv',Eorbarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'Lorbarr.csv',Lorbarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(rundir+runoutprefix+'porbarr.csv',porbarr,delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit to E convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, rval, pval, stderr = stats.linregress(np.log10(dtarr),np.log10(deltaEarr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.2851907499866377 -0.22792849025589357 0.9301704934577777 0.02191746292282565\n"
     ]
    }
   ],
   "source": [
    "print(slope, intercept, rval, pval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4RklEQVR4nO3dd3gU1frA8e9rggSEgCgqgjQFhUCIEMEKWEAEFREVBbH+LnK9Kl5FRbEgXQUBUfGiWLDrFcVeQEERuRqkSBFUOgREegsleX9/nEnc1N0ku5kk+36eZ5/slJ15z+5m3tlzZs4RVcUYY0z0OczvAIwxxvjDEoAxxkQpSwDGGBOlLAEYY0yUsgRgjDFRyhKAMcZEKUsAptwRkWNF5FsR2SUio8O87UEi8loxXr9YRNqHL6KyIZKfiSk6SwAlTER6ikiKiOwWkVQR+UxEzvY7rnKmD/AXEK+qd+dcKCIvi8gB7zPYKiJficgp4Q7C28/QwHmqmqCqM8K8n/oiol55Ah89wrmfYirMZ5L5WFDyYUYXSwAlSETuAsYCw4FjgbrAs0BXH8PKRkRi/Y4hDOoBS7TguxwfV9UqQG1gPTCpRCKLrOqqWiXg8XZeK4lITI7pQn3mRfyOhPyZBDxaFGE/BSon3+/wUVV7lMADqAbsBq4sYJ2KuASxwXuMBSp6y9oD64C7gT+BVOBGb1kbYCMQE7CtbsBC7/lhwADgD2AL8A5Qw1tWH1DgZmAN8C0QA4zGnbGtBG7z1okNKMskL4b1wNDMfQM3ALOAUcA27/UXBcRVA3jJK9824IOAZRcD84HtwGwgsYD36kzgJ2CH9/dMb/7LwEHggPd+X5DHa18GhgZMdwb2BEwfD7wHbPbivyNg2SDgtYDpd733fof33iV48/vkiOMjb/4q4AJvH/syPwdv2anee17Bm74JWOq9T18A9fJ5LzI/w9h8lr8MTAA+BfZ4+18F3AcsBPYDscClwGLv/Z8BNAnYRq71I/mZ5FO+63Hf0b+AgQHLw/b9Bq4E5ubY/13AVL+PIRE5LvkdQLQ8gE7Aofz+Sb11BgNzgGOAmriD4BBvWXvv9YOBCriD1l7gSG/5H0CHgG29CwzwnvfztlsHl2T+A7zpLcv8B5kMHAFUAvoCS7z1jwSmkT0BvO9t4wgv1h+BW7xlN3j/7P/w/tH+iTvYi7f8E+Btb7sVgHbe/FNxia2N97rrcQedinm8TzVwB8Xe3j/tNd70Ud7yfA8mOZd7ZXgVWOBNHwbMBR4GDgcaAiuAC73lg8ieAG4CqvJ38p6f134C5q3COwACXwP/CFj2BPCc97wr8DvQxCvjg8DsfMqT+RkWlAB2AGd55Yvz4pgPnOB95o1xyaGD97nc6+3/8IC4s9aP5GdSQPme92JtgUtCTcL9/fZev5XsyW8e0N3vY0hEjkt+BxAtD6AXsDHIOn8AnQOmLwRWec/b484YYwOW/wmc7j0fCrzoPa/q/TPX86aXAucHvK4W7iAdG/AP0jBg+dd4B3Rv+oKAf5BjvX++SgHLrwG+8Z7fAPwesKyy99rjvP1m4CWtHGWfgJfsAuYtw0sQOeb3Bn7MMe8H4AbveSgHmzTcmW4G7iww0VvWBliTY/37gZe854MISAA51qvulbVafnGQPQH8H/C191yAtUBbb/oz4OaA1x2GS/j18thv5me4PcejSUAck/OI46aA6YeAd3Lsbz3QPq/1I/yZZD5eyVG+OgHr/whcHe7vd8B3cZj3PAGXyHKdiJSHh7UBlJwtwNFB6iCPB1YHTK/25mVtQ1UPBUzvBap4z98ALheRisDlwM+qmrmtesD7IrJdRLbj/mHScQfzTGtzxLE2n2X1cGeIqQHb+w/ul0CmjZlPVHWv97QK7uxxq6puy1106gF3Z27T2+4JOcofGN/qHPNW4+rzQzVKVavjDhD7gJMD4jg+RxwPkP29AlxduoiMFJE/RGQn7iAJcHSIMbwHnCEitYC2uGT0XUAc4wJi2IpLEgWV8WhVrR7wWBqwbG0e6+f8zLPeU1XN8JbXzmf9nML2mQQ8rs+xfGPA88Dvfji/3wCvAD1FRHCJ7R1V3V+IcpQZlgBKzg+4M+fLClhnA+7LnKmuNy8oVV2C+4e7COiJSwiZ1uLq4QP/ueJUdX3gJgKep+J+Hmc6Ice29pP9YBOvqgkhhLkWqCEi1fNZNixHjJVV9c081s35PoF7r9bnsW6BVHUNrgphnIhU8uJYmSOOqqraOY+X98RV1VyAaxep782XzM0H2fc24Eugh7ett9Q77fTiuCVHHJVUdXZhy1hALIHzsr2n3sHvBLK/pwWVJ2yfSRGE8/uNqs7BtVecg/tcXo1Q3L6zBFBCVHUHrl75GRG5TEQqi0gFEblIRB73VnsTeFBEaorI0d76hbnm/A3cwawtrg0g03PAMBGpB+Btv2sB23kH6Ccitb2D9X0B5UjFHbRGi0i8iBwmIieKSLtgwXmv/Qx4VkSO9Mrf1lv8PNBXRNqIc4SIdBGRqnls6lOgsXdJbax3uWNT4ONgMeQT11e4A1gfXNXCLhG5T0QqeWf5zUTktDxeWhWXDLfgqrqG51i+CdeGUJA3gOuAK8ietJ8D7heRBAARqSYiVxayaIXxDtBFRM4XkQq4iw3249qhQhHWz6SQwvb9DjAZeBo4qKqzwh1waWEJoASp6mjcFQUP4q4wWYu7AuEDb5WhQAruSotfgJ+9eaF6E2iHq1f+K2D+OOBD4EsR2YVrMGtTwHaexx3kF+IawD7FNUCne8uvwzWQLsHVj/4XV+8ait64+tlfcW0YdwKoagqu4fhpb5u/49oTclHVLbgrhu7GHXzvBS7OUebCesLbTqy37SRc28BfwAu4M/ycJuN+da3HvRdzciyfBDT1qiY+yGe/HwKNcO1DWde9q+r7wGPAW1710iLcr7uCbM9xHf1dQdbPoqrLgGuB8bgyXwJcoqoHQnx9OD6Te3PEH+prw/39BnfW34zCnYCVOfL3L05j8iYiF+GuTsn5E9+YMi+v77dXHfgn0FJVf/MtuAizXwAmF6/qo7P3U7428Aju0k9jyrwQv9//BH4qzwd/sF8AJg8iUhmYCZyCu0LmE6Cfqu70NTBjwiDY91tEVuEa8i9T1Xl+xVkSLAEYY0yUsiogY4yJUmWqY6Sjjz5a69ev73cYxhhTpsydO/cvVa2Zc36ZSgD169cnJSXF7zCMMaZMEZGcd2kDVgVkjDFRyxKAMcZEKUsAxhgTpcpUG0BeDh48yLp160hLS/M7FGNKTFxcHHXq1KFChQp+h2LKsDKfANatW0fVqlWpX78+rgNDY8o3VWXLli2sW7eOBg0a+B2OKcPKfBVQWloaRx11lB38TdQQEY466ij71RtFUnel0u7ldmzcvTH4yoVQ5hMAYAd/E3XsOx9dhnw7hFlrZjF45uCwbrfMVwEZY0x5VWlYJdIO/f1Lb0LKBCakTCAuNo59A/cVe/vl4heA32JiYkhKSiIhIYEWLVowevRoMjIyCnzNqlWreOONNwpcxxgT3VbcsYKezXpSObYyAJVjK9OreS9W9lsZlu1HZQIId31apUqVmD9/PosXL+arr77is88+49FHHy3wNZYAjDHB1Kpai/iK8aSlpxEXG0daehrxFeM5rspxYdl+VCaASNWnARxzzDFMnDiRp59+GlVl1apVnHPOObRs2ZKWLVsye7YbYW/AgAF89913JCUlMWbMmHzXM8ZEt017NtG3VV/m3DyHvq36hrUhuEx1B52cnKw5+wJaunQpTZo0Cen1OevTMhW3Pq1KlSrs3r0727zq1auzbNkyqlatymGHHUZcXBy//fYb11xzDSkpKcyYMYNRo0bx8cduyNS9e/fmuZ4x+SnMd99ENxGZq6rJOedHVSPwijtW0P/L/nzw6wfsPbSXyrGV6dakG6M6jorYPg8ePMhtt93G/PnziYmJYfny5cVazxhjwiWqEkCk69MyrVixgpiYGI455hgeffRRjj32WBYsWEBGRgZxcXF5vmbMmDEhrWeMMeESdW0AkaxPA9i8eTN9+/bltttuQ0TYsWMHtWrV4rDDDuPVV18lPT0dgKpVq7Jr166s1+W3njHGREpU/QIAmNJjStbzZ7o8E5Zt7tu3j6SkJA4ePEhsbCy9e/fmrrvuAuDWW2+le/fuTJ48mU6dOnHEEUcAkJiYSExMDC1atOCGG27Idz1jjImUqGoENqY8se++CVV+jcBRVwVkjDHGsQRgjDFRyhKAMcZEKUsAxhgTpSwBGGNMlLIEYIwxUcoSQBhkdged+Vi1ahVnnnkmUDp7/Xzuuedo3rw5SUlJnH322SxZsiTXOmvXruXcc8+ladOmJCQkMG7cuKxl8+fP5/TTTycpKYnk5GR+/PFHAGbMmEG1atWy3ofBg//ubO+mm27imGOOoVmzZtn2s2DBAs444wyaN2/OJZdcws6dO7OWjRgxgpNOOomTTz6ZL774AoBly5Zle6/j4+MZO3Zstm2OHj0aEeGvv/4CYNu2bXTr1o3ExERat27NokWLstYdN24czZo1IyEhIdt28itjZjkzu/9u165d1vzt27dzxRVXcMopp9CkSRN++OEHAB566CESExNJSkqiY8eObNiwAYAnnngiqxzNmjUjJiaGrVu3hlRGY8JCVUv8AVwJLAYygORQX9eqVSvNacmSJbnmlbQjjjgi32XffPONdunSpQSjCW7Hjh1Zz6dOnaoXXnhhrnU2bNigc+fOVVXVnTt3aqNGjXTx4sWqqtqhQwf99NNPVVX1k08+0Xbt2qlqwWWdOXOmzp07VxMSErLNT05O1hkzZqiq6qRJk/TBBx9UVdXFixdrYmKipqWl6YoVK7Rhw4Z66NChbK89dOiQHnvssbpq1aqseWvWrNGOHTtq3bp1dfPmzaqq2r9/fx00aJCqqi5dulTPO+88VVX95ZdfNCEhQffs2aMHDx7U888/X3/77bcCy7ht2zZt0qSJrl69WlVVN23alLXv6667Tp9//nlVVd2/f79u27Yt1/s9btw4veWWW3K9Px9++KGee+65uebnVcZMpeG7b8oGIEXzOKb69QtgEXA58K1P+4+4KlWqALm7fS7Ia6+9RuvWrUlKSuKWW24hPT2dn376icTERNLS0tizZw8JCQksWrSIGTNm0LZtW7p06cLJJ59M3759gw5Ckyk+Pj7r+Z49e/IcXrBWrVq0bNkScN1WNGnShPXr1wNuOMLMM/UdO3Zw/PHHB91n27ZtqVGjRq75y5cvp23btgB06NCB9957D4CpU6dy9dVXU7FiRRo0aMBJJ52U7SwcYPr06Zx44onUq1cva96///1vHn/88WxlWrJkCeeddx4Ap5xyCqtWrWLTpk0sXbqUNm3aULlyZWJjY2nXrh1TpkwpsIxvvPEGl19+OXXr1gVc99+Z63z77bfcfPPNABx++OFUr14dCO39fvPNN7nmmmtyzc+rjMaEiy9dQajqUgj/uKZ33gnz54d1kyQlQbBf35ldQQA0aNCA999/P2vZyJEjs3X7nJ+lS5fy9ttv8/3331OhQgVuvfVWXn/9da677jouvfRSHnzwQfbt28e1115Ls2bNmDFjBj/++CNLliyhXr16dOrUiSlTpnDFFVfQo0cPli1blmsfd911F9dddx0AzzzzDE8++SQHDhzg66+/LjC2VatWMW/ePNq0aQPA2LFjufDCC+nfvz8ZGRnZxi744YcfaNGiBccffzyjRo0iISGhwG0nJCQwdepULrvsMt59913Wrl0LwPr16zn99NOz1qtTp05WAsr01ltvZTtoTp06ldq1a9OiRYts67Vo0YIpU6Zwzjnn8OOPP7J69WrWrVtHs2bNGDhwIFu2bKFSpUp8+umnJCcnF1jG5cuXc/DgQdq3b8+uXbvo168f1113HStXrqRmzZrceOONLFiwgFatWjFu3LisLj0GDhzI5MmTqVatGt988022+Pbu3cvnn3/O008/nev9yVlGY8Kp1LcBiEgfEUkRkZTNmzf7HU6eMkcEmz9/fraDf2FMnz6duXPnctppp5GUlMT06dNZsWIFAA8//DBfffUVKSkp3HvvvVmvad26NQ0bNiQmJoZrrrmGWbNmAfD2229nxRP4yDz4A/zrX//ijz/+4LHHHmPo0KH5xrV79266d+/O2LFjs85kJ0yYwJgxY1i7di1jxozJOutt2bIlq1evZsGCBdx+++1cdtllQcv94osv8uyzz9KqVSt27drF4YcfHtL7deDAAT788EOuvPJKwB1Ehw8fnq3dIdOAAQPYvn07SUlJjB8/nlNPPZWYmBiaNGnCfffdR8eOHenUqRNJSUnExMQUWMZDhw4xd+5cPvnkE7744guGDBnC8uXLOXToED///DP//Oc/mTdvHkcccQQjR47MimHYsGGsXbuWXr165TrQf/TRR5x11lm5fiHlLKMxYZdXvVA4HsA0XFVPzkfXgHVmUE7bADLnhdoG8NRTT+mAAQPyXLZhwwZt2LChNmnSRHfv3p213bZt22atM2nSJL3zzjtVVfWqq67SFi1a5Hq88sorubadnp6u8fHxee73wIED2rFjRx09enS2+fHx8ZqRkaGqqhkZGVq1atU8X1+vXr2senhV1ZUrV+ZqAwi0bNkyPe2001RVdfjw4Tp8+PCsZR07dtTZs2dnTX/wwQfaoUOHrOmFCxdqzZo1tV69elqvXj2NiYnRE044QVNTU7PtIyMjQ+vVq5etXj7T/fffr88880yBZRwxYoQ+/PDDWa+56aab9J133tHU1FStV69e1vxvv/1WO3funGsfq1evzvUeXHbZZfr666/nWjdnGXMqDd99UzaQTxuAL43AWTuPggSQkpKS7UC9bt26rEbIQIsXL9aTTjopq1Fxy5YtWQ1/l1xyib7++us6dOhQ/de//qWqLgHExcXpihUrND09XTt27Kj//e9/Q4p3+fLlWc8//PBDzet9zcjI0N69e2u/fv1yLTvllFP0m2++UVXVadOmacuWLVVVNTU1Neug+b///U9POOGErGnVvBNAZnnT09O1d+/eOmnSJFVVXbRoUbZG4AYNGmRrBO7Ro4e++OKL+ZYxMPls27ZN9+/fr6qqEydO1N69e+fa/+rVq/Xkk0/OarjNr4xLlizR8847Tw8ePKh79uzRhIQE/eWXX1RV9eyzz9Zff/1VVVUfeeQR7d+/v6pmf7+feuop7d69e9b09u3b9cgjj8xK7IGClbE0fPdN2WAJIIIKSgAHDhzQc889VxMTE/XJJ5/Un376STt27Jjndt566y1t0aKFNm/eXFu2bKk//PCDvvLKK3r55ZerqrsipHXr1jp9+nT95ptv9JxzztHOnTtr48aN9ZZbbtH09PSQ4r3jjju0adOm2qJFC23fvr0uWrRIVVXXr1+vF110kaqqfvfddwpo8+bNs35BfPLJJ1nLWrZsqYmJidq6dWtNSUlRVdXx48dr06ZNNTExUdu0aaPff/991j6vvvpqPe644zQ2NlZr166tL7zwgqqqjh07Vhs1aqSNGjXS++67L1vCGDp0qDZs2FAbN26cdUWOquru3bu1Ro0aun379nzLGJgAZs+erY0aNdLGjRtrt27ddOvWrVnrnX322dqkSRNNTEzUadOmZc3Pr4yqqo8//rg2adJEExISdMyYMVnz582bp61atdLmzZtr165ds/Zz+eWXa0JCgjZv3lwvvvhiXbduXdZrXnrpJe3Ro0eu+EMpY2n47puyIb8E4Et30CLSDRgP1AS2A/NV9cJgrysP3UE//fTT1K1bl0svvbRY28k5prCJPmXtu2/8U6rGBFbV94GitZaWcbfddpvfIRhjDBCFI4KVF+3bt6d9+/Z+h2GMKcNK/WWgofCjGssYP9l33oRDmU8AcXFxbNmyxf4hTNRQVbZs2UJcXJzfoZgSsHYt3HILbNsW/m2X+SqgOnXqsG7dOkrrTWLGREJcXBx16tTxOwwTQbt2wciR8OSToApdu0LnzuHdR5lPABUqVKBBgwZ+h2GMMWFx6BC88AI88gj8+Sf06gXDhkEkuoMq8wnAGGPKA1X47DO45x5YsgTOOQc+/hhOOy1y+yzzbQDGGFPWLVgAHTtCly5w4ABMmQIzZ0b24A+WAIwxxjcbNsBNN8Gpp8LPP8O4cbB4MXTrBmHuLDlPVgVkjDElbM8eeOIJ9zh0CO6+Gx54AI48smTjsARgjDElJD0dXn4ZHnoIUlPhqqtgxAho2NCfeCwBGGNMCfjqK+jfHxYuhDPOgPfec3/9ZG0AxhgTQYsXw0UXuUbeXbvgnXfg++/9P/iDJQBjjImITZvcHbyJiTBnDowaBUuXwpVXlkwDbyisCsgYY8Jo71539+5jj0FaGtx+u6vzP+oovyPLzRKAMcaEQUYGvPYaDBwI69a5SzkfewwaNfI7svxZFZAxxhTTjBnupq3rr4fjjnM3cU2ZUroP/mAJwBhjiuzXX+HSS+Hcc+Gvv+D11+F//4O2bf2OLDSWAIwxppA2b4bbboNmzdzZ/siRsGwZ9OwJh5Who6q1ARhjTIjS0lx3DcOHu7t5b7kFBg2CmjX9jqxoLAEYY0wQGRnw1ltw//2wZg1ccolr4G3SxO/IiqcM/VgxxpiSN2sWnH6665f/qKNg+nT48MOyf/AHSwDGGJOn33+H7t1dv/wbNrg+fFJS4Lzz/I4sfCwBGGNMgK1b4c47oWlT+OILGDIEli93l3iWpQbeUFgbgDHGAPv3wzPPuAP+zp1w880weLC7rr+8Kmf5zBhjCkcV3n3XnfHffber71+wACZOLN8Hf7AEYIyJYnPmwNlnu375jzjCVfl89pm7vj8aWAIwxkSdlSvh6qtdl8wrVsALL8C8ea7L5mhibQDGmKixfTsMGwZPPQUxMfDww3DPPVClit+R+cMSgDGm3Dt4ECZMgEcfhW3b4IYbXGNv7dp+R+YvqwIyxpRbqvDBB5CQAP36wamnws8/w4sv2sEfLAEYY8qplBRo3971yx8bCx9/7MblTUryO7LSwxKAMaZcWbMGrr3W9c+/dKmr+lm4ELp0KT1DMZYW1gZgjCkXdu503TKPGeOm778fBgyA+Hh/4yrNLAEYY8q0Q4fg+efhkUdcP/3XXuuu9Klb1+/ISj9fEoCIPAFcAhwA/gBuVNXtfsRijCmbVOHTT91lnEuXulG4Pv0UkpP9jqzs8KsN4CugmaomAsuB+32KwxhTBs2fDx06wMUXQ3q6u9Jnxgw7+BeWLwlAVb9U1UPe5Bygjh9xGGPKlvXr4cYboWVLlwSeegoWLYKuXa2BtyhKQxvATcDb+S0UkT5AH4C6VqlnTFTavRsefxxGjXJn/P37wwMPQPXqfkdWtoWUAEQkGTgHOB7YBywCvlLVbQW8ZhqQV196A1V1qrfOQOAQ8Hp+21HVicBEgOTkZA0lXmNM+ZCeDi+9BA89BBs3Qo8eMGIENGjgd2TlQ4EJQERuBG4HVgJzgWVAHHA2cJ+ILAIeUtU1OV+rqhcE2fYNwMXA+apqB3ZjTDZffOHO9BctgjPPhPffd101m/AJ9gugMnCWqu7La6GIJAGNgFwJoCAi0gm4F2inqnsL81pjTPm2aJE78H/xBTRs6Prq797d6vgjocAEoKrPBFk+v4j7fRqoCHwl7lOdo6p9i7gtY0w5sHGj651z0iSoVg2efBJuvRUqVvQ7svIrWBXQO6p6lff8MVW9L2DZl6papN6zVfWkorzOGFP+7N0Lo0fDY4/BgQNwxx2uzr9GDb8jK/+CXQbaKOB5hxzLaoY5FmNMFMnIgFdegcaN3Zl/p06wZInrysEO/iUjWAIoqHHWGm6NMUXy9dfQqpXrl792bfjuO/jvf+EkqxsoUUEbgUXkVFyiqOQ9F+9RKdLBGWPKl6VL4d57XdfM9erBG2+4SzsPs36JfREsAaQCT3rPNwY8z5w2xpig/vwTBg2CiRPd4OsjR7oBWuLi/I4sugVLAL1UdUOJRGKMKXf27YNx42D4cNfY27ev67WzprUglgrBEsALIlIDmAF8DswK6MPHGGPylJEBb77pumtYswYuvdR15XDyyX5HZgIVWPOmqp2B9rgE0A2YIyJTRKSPiFjHPMaYXL77zt2xe+21cPTRrsF36lQ7+JdGQfsCUtU03Nn/5wAi0gC4CHhaRI5T1daRDdEYUxb89hvcd5/rsqFOHZg8GXr1sgbe0izYjWCnqOqv3vOKqrpfVVcCz4rIPFz/QMaYKLZlCwweDM8+6xp1hw6Ff/8bKlf2OzITTLDc/EbA8x9yLHtGVQ+EOR5jTBmxf7/rnvnEE+Hpp+Hmm+H332HgQDv4lxXBqoAkn+d5TRtjooCq66BtwABYuRIuugieeAISEvyOzBRWYe4Eznnnr90JbEyU+eEH1zVzjx5QtSp8+aUbh9cO/mVTsF8AdUTkKdzZfuZzvOnaEY3MGFNqrFjhzvjffRdq1XI9dl5/PcTE+B2ZKY5gCeCegOcpOZblnDbGlDPbtsGwYTB+PMTGupu4+veHKlWKvs3UXalc/d7VvH3F2xxXJa9BA01JCTYewCslFYgxpvQ4cAAmTHBX92zb5gZiHzIEjj+++Nse8u0QZq2ZxeCZg3m2y7PF36ApMgllNEYRqQncBzTFDQkJgKqeF7nQcktOTtaUFPvhYUykqMIHH7gO237/HTp0cFf6JCYWf9uVhlUi7VBarvlxsXHsG5jnoIMmTERkrqom55wf6i0arwNLgQbAo8Aq4KewRWeM8d1PP0G7dnD55XD44a5x94svwnPwB1hxxwp6NutJ5Vh3jWjl2Mr0at6Llf1WhmcHptBCTQBHqeok4KCqzlTVm4ASPfs3xkTG6tXujt3WrWHZMnjuOViwwF3eGc5xeGtVrUV8xXjS0tOIi40jLT2N+Irx1g7go6BdQXgOen9TRaQLsAGwMXuMKcN27IARI2DsWHegHzjQdeVQtWrk9rlpzyb6tupLn1Z9mDh3Iqm7UyO3MxNUqG0AFwPfAScA44F4YJCqfhTZ8LKzNgBjiu/gQXj+eXdFz19/Qe/e7kqfE07wOzITKfm1AYT6C2Cbqu4AdgDnehs8K4zxGWMiTNWNxHXvvfDrr9C+vRuMvWVLvyMzfgm1DWB8iPOMMaXQvHlw/vmuX/6MDNc989df28E/2gXrDfQM4EygpojcFbAoHrB7AI0p5datgwcfdF0zH3WU67StTx+oUMHvyExpEKwK6HCgirdeYNPQTuCKSAVljCmeXbvcCFyjR7sz/nvucaNzVavmd2SmNAl2J/BMYKaIvKyqq0soJmNMER06BC++CA8/DJs2wTXXuPF469f3OzJTGgWrAvoIr9dPyeOCYFW9NDJhGWMK6/PPXT89ixfDWWe5ev42bfyOypRmwaqARpVIFMaYIlu40FXxfPmlG5zlvfegW7fw3sRlyqdQqoAAEJFKQF1VXRbxqIwxQaWmwkMPwUsvubr9MWPg1ltdNw7GhCKky0BF5BJgPn8PDJ8kIh9GMC5jTD727HG9dDZq5K7u6dfPddx255128DeFE+p9AIOA1sB2AFWdj+sYzhhTQtLT3dl+48buLt6LLoKlS+HJJ6GGdcxiiiDUBHDQuxM4kA0JaUwJmT4dWrWCm26CunXh++/d6Fwnnuh3ZKYsCzUBLBaRnkCMiDQSkfHA7AjGZYwBliyBiy+GCy5wnbe99RbMnu3G5TWmuEJNALcDCcB+4E3cjWB3RigmY6Lepk3wz3+6vvhnzYInnnDVPT162NU9JnxC6gxOVfcCA72HMSZC9u1zV/OMHOme33qru6nr6KP9jsyURyHfCJaXot4IJiJDgK5ABvAncIOqbijKtowpDzIy4I03XHcNa9dC166uK4fGjf2OzJRnwaqARgGjgZXAPuB577Eb+KMY+31CVRNVNQn4GHi4GNsypkybOdONxtW7NxxzDMyY4cbltYO/ibSQbgQTkdE5BhP4SESKPDKLqu4MmDwCu6LIRKHly13f/FOnusFYXn0VevaEw0JtmTOmmEIdEOYIEWmoqisARKQB7sBdZCIyDLiOgEFm8lmvD9AHoG7dusXZpTGlwl9/uRu5JkyASpVcZ2133umeG1OSQh0SshMwEVgBCFAP6KOqXxbwmmlAXqM9D1TVqQHr3Q/EqeojweKwISFNWZaWBuPHu+EXd+92/fIPGuSqfYyJpGINCamqn4tII+AUb9avqro/yGsuCDG214FPgaAJwJiySBXefhvuvx9WrYIuXVwDb9Omfkdmol2BtY0icnbmc1Xdr6oLvMd+b3m8iDQr7E69ZJKpK/BrYbdhTFnw/fdwxhmuX/5q1WDaNDcurx38TWkQ7BdAdxF5HNcJ3FxgMxAHnISrt68H3F2E/Y4UkZNxl4GuBvoWYRvGlFp//AH33ee6Zj7+eNeHT+/eEGMDqZpSJNhVQP8WkRpAd+BKoBbuctClwH9UdVZRdqqq3YvyOmNKu61bYehQN/bu4YfDo4/C3XfDEcW6ZMKYyAjaBqCqW/n7+n9jTB4OHIBnn3VX9+zY4TptGzwYatXyOzJj8hc0AYjIhcBlQG1v1npgqqp+HsG4jCkTVGHKFFfd88cf0LGj67cnMdHvyIwJLlhXEGOBxsBkYJ03uw5wh4hcpKr9IhueMaXXjz+66p1ZsyAhAT77DDp18jsqY0IX7BdAZ1XNdUO6iLwNLAcsAZios2qVu6Tzrbfg2GNh4kS48UaIDfW2SmNKiWA3naeJyGl5zD8NSItAPMaUWtu3u64bTj7Zdd/w0EPw22/wj3/Ywd+UTcG+tjcAE0SkKn9XAZ2A677hhsiFZUzpcfAg/Oc/7q7drVvhuuvclT516vgdmTHFE+wy0J+BNiJyHAGNwKq6MeKRGeMzVfjoI3fWv2wZnHsujB4Np57qd2TGhEdI/Q6q6kZVnes9NgKIyCnBXmdMWTV3Lpx3nuuXX8QlgunT7eBvypfidDybb0dwxpRVa9e6Kp7kZFi0CJ55BhYudOPy2lCMprwJdhnoU/ktAqqHPRpjfLJrlxuG8cknXdXPgAHuUa2a35EZEznBGoFvxPX1k1fPn9eEPxxjStahQzBpkht3988/oVcv111zvXp+R2ZM5AVLAD8Bi1R1ds4FIjIoIhEZUwJU3Y1b99wDS5bAOee4XjpPy+uiZ2PKqWBtAFcA8/NaoKoNwh6NMSVgwQLXZUOXLq4PnylT3Li8dvA30abABKCqW1V1b0kFY0wkbdjgOmk79VT4+WcYNw4WL4Zu3ayB10SnkO5fFJFfyD1w+w4gBRiqqlvCHZgx4bJnj+ug7YknXJ3/XXfBwIFw5JF+R2aMv0K9gf0zIB14w5u+GqgMbAReBi4Je2TGFFN6OrzyCjz4IKSmwlVXwYgR0LCh35EZUzqEmgAuUNWWAdO/iMjPqtpSRK6NRGDGFMdXX0H//u4a/tNPdyNznXGG31EZU7qEeiNYjIi0zpzwOojLHNzuUNijMqaIFi+Gzp1dI++uXfDOOzB7th38jclLqL8A/g94UUSq4G4C2wncLCJHACMiFZwxodq0yV3L/8ILEB8Po0bBbbdBxYp+R2ZM6RVSAlDVn4DmIlLNm94RsPidSARmTCj27nV37z72GKSlwe23u26ajzrK78iMKf1CvQqoGvAI0NabngkMzpEIjCkxGRnw2mvuap5169ylnI89Bo0a+R2ZMWVHqG0ALwK7gKu8x07gpUgFZUxBZsxwN21dfz0cd5y7iWvKFDv4G1NYobYBnKiq3QOmHxWR+RGIx5h8/fqr65v/o4+gbl14/XW4+mo4rDh92hoTxUL919knImdnTojIWcC+yIRkTHabN7sG3WbN3Nn/iBEuGfTsaQd/Y4oj1F8AfYHJmY3AwDbg+siEZIyTlua6axg+3N3Ne8st8MgjcMwxfkdmTPkQ6lVAC4AWIhLvTe8UkTuBhRGMzUQpVXjrLbj/fli9Gi65xDXwNmnid2TGlC+F+gGtqjtVdac3eVcE4jFRbtYsd+duz55Qo4YbhvHDD+3gb0wkFKcG1fpPNGHz++/Qvbvrl3/9enj5ZUhJcePyGmMiozgJIGfvoMYU2tatcOed0LQpfPEFDBkCy5e7SzytgdeYyCrwX0xEdonIzjweu4DjSyhGUw7t3+/u4D3xRBg/Hm64wf0KePBBqFzZ7+gKJ3VXKu1ebsfG3Rv9DsWYQgk2IExVVY3P41FVVUO9gsiYLKrw7rvujP/uu119/4IFMHGiu6mrLBry7RBmrZnF4JmD/Q7FmEIR1bJTk5OcnKwpKSl+h2GKaM4cd9CfPRuaN3cdtnXs6HdURVdpWCXSDqXlmh8XG8e+gXabjCk9RGSuqibnnG+1rCbiVq50d+yecQasWOF67Jw3r2wf/AFW3LGCns16UjnW1VlVjq1Mr+a9WNlvpc+RGRMaXxOAiNwtIioiR/sZh4mM7dvhnnvglFPcpZwPPwy//QY33wwxMUFfXurVqlqL+IrxpKWnERcbR1p6GvEV4zmuShmtyzJRx7d6fBE5AegIrPErBhMZBw/ChAnw6KOwbZtr4B0yBGrX9juy8Nu0ZxN9W/WlT6s+TJw7kdTdqX6HZEzIfGsDEJH/AkOAqUCyqv4V7DXWBlC6qcLUqa7Dtt9+g/PPd/X8SUl+R2ZMdCtVbQAi0hVY73UxYcqBlBRo3971yx8bCx9/7MbltYO/MaVXxKqARGQakFdl6EDgAVz1Tyjb6QP0Aahbt27Y4jPhsWYNPPCA65q5Zk1X9fN//+eSgDGmdCvxKiARaQ5MB/Z6s+oAG4DWqlrgnTRWBVR67NwJI0fCmDFu+t//hgED3Hi8xpjSJb8qoBI/T1PVX4CsDn1FZBUhtgEY/x06BM8/77pl3rwZrr0Whg1zA7QYY8oW+6FuQqIKn37qLutcuhTatnXTybnOKYwxZYXvN4Kpan07+y/d5s+HDh3g4oshPR0++MCNzGUHf2PKNt8TgCm91q+HG2+Eli1dEnjqKVi0CLp2BbHOwI0p86wKyOSyezc8/ri7hj89Hfr3d1f6VK/ud2TGmHCyBGCypKfDSy/BQw/Bxo3Qo4cbgL1BA78jM8ZEgiUAA7jBWPr3d1U8Z54J77/vumo2xpRf1gYQ5RYtgk6d3GPvXtdXf+a4vMaY8s0SQJTauBH69IEWLeB//4PRo2HJErjiCmvgNSZaWBVQlNm71x3sH3sMDhyAO+5wdf41avgdmTGmpFkCiBIZGTB5MgwcCBs2QPfuriuHk07yOzJjjF+sCigKfP01tGrlrumvUwe++w7++187+BsT7SwBlGNLl8Ill7h++bdtgzfegB9+gLPP9jsyY0xpYAmgHPrzT7j1Vjfw+rffuvr+X3+Fa66Bw+wTN8Z4rA2gHNm3D8aNg+HDXWNv376u186aNf2OzBhTGtn5YDmQkeEGZDnlFLj/fjj3XFi8GJ5+umgH/9RdqbR7uR0bdxc4PIMxpoyzBFDGffsttGnj+uU/+mjX4Dt1Kpx8ctG3OeTbIcxaM4vBMweHL1BjTKkTFQmgPJ7RLl/uxt9t187d1DV5Mvz0kzv7L6pKwyohjwoTUiaQoRlMSJmAPCpUGlYpfIEbY0qNqEgA5emMdssW6NcPEhJg2jQYOhSWLYPevYvfwLvijhX0bNaTyrGVAagcW5lezXuxst/KMERujCltynUjcKVhlUg7lJY1PSFlAhNSJhAXG8e+gft8jKzw9u+H8ePdAX/XLvjHP+DRR+HYY8O3j1pVaxFfMZ609DTiYuNIS08jvmI8x1U5Lnw7McaUGuX6F0B5OKNVhXfegSZN3HCMZ54JCxfCc8+F9+CfadOeTfRt1Zc5N8+hb6u+5arazBiTXbn+BVDWz2hnz4a774Y5cyAxEb780g3NGElTekzJev5Ml2ciuzNjjK/K9S8AKJtntH/8AVddBWedBatXw6RJ8PPPkT/4G2Oii6iq3zGELDk5WVNSUvwOI2K2bXN1/OPHQ4UKrsqnf3+oUsXvyIwxZZmIzFXV5Jzzy3UVUFlx4ABMmACDB7skcOONMGQIHH+835EZY8qzcl8FVJqpwpQp7pLOO++Eli1h3jxX5WMHf2NMpFkC8MlPP7mbuLp3h8MPh08/dY28LVr4HZkxJlpYAihhq1dDr17QurW7geu552DBArjoIhuK0RhTsqwNoITs2AEjRsDYse5AP3Ag3HcfVK3qd2TGmGhlCSDCDh6EiRNh0CD46y/XZcOwYXDCCX5HZoyJdlYFFCGq8NFHblCW226DZs1g7lzXaZsd/I0xpYElgAj4+Wc3DOOll7pEMHWq66a5ZUu/IzPGmL9ZAgijdevg+ushOdn11zN+PCxa5BKBNfAaY0obawMIg1274PHHYfRoSE93d/A+8ABUq+Z3ZMYYkz9LAMVw6BC8+CI8/DBs2uQGXR8+HOrX9zsyY4wJzhJAEajC55+7M/3Fi12nbVOnuqEZjTGmrLA2gEJauBAuvBA6d4a0NHjvPfjuOzv4G2PKHl8SgIgMEpH1IjLfe3T2I47C2LAB/u//ICkJUlJgzBhYsgQuv9waeI0xZZOfVUBjVHWUj/sPyZ49MGqUa+Q9eNB12vbgg1Cjht+RGWNM8VgbQD7S091NWw8+6M7+r7gCRo6EE0/0OzJjjAkPP9sAbhORhSLyoogc6WMcuUybBq1awU03ubt2Z82Cd9+1g78xpnyJWAIQkWkisiiPR1dgAnAikASkAqML2E4fEUkRkZTNmzdHKlzAXdHTubMbenHHDnjrLfjhB3eVjzHGlDe+DwkpIvWBj1W1WbB1IzUk5KZN8Mgj8PzzrnfOgQPh9tshLi7suzLGmBJXqoaEFJFaqprqTXYDFvkRx7597mqeESPcJZ3/+pe7qevoo/2IxhhjSpZfjcCPi0gSoMAq4JaS3HlGBrz+uuuuYd066NrVXeXTuHFJRmGMMf7yJQGoam8/9gswcybcfbfrmrlVK3jtNTc0ozHGRJuouRN42TJ3pt++Pfz5J7z6Kvz4ox38jTHRKyoSwNChbkCWb75xnbUtWwbXXguHRUXpjTEmb1FxI1iDBq4bh0GD4Nhj/Y7GGGNKh6hIAL16uYcxxpi/WSWIMcZEKUsAxhgTpSwBGGNMlLIEYIwxUcoSgDHGRClLAMYYE6UsARhjTJSyBGCMMVHK9/EACkNENgOr/Y4jBEcDf/kdRBiUl3JA+SlLeSkHWFlKUj1VrZlzZplKAGWFiKTkNfhCWVNeygHlpyzlpRxgZSkNrArIGGOilCUAY4yJUpYAImOi3wGESXkpB5SfspSXcoCVxXfWBmCMMVHKfgEYY0yUsgRgjDFRyhJAGIjIlSKyWEQyRKTAS8FEJEZE5onIxyUVX6hCKYeInCAi34jIEm/dfiUdZyhC/UxEpJOILBOR30VkQEnGGAoRqSEiX4nIb97fI/NZ73GvvEtF5CkRkZKONZhClKWuiHzplWWJiNQv4VCDCrUs3rrxIrJORJ4uyRhDYQkgPBYBlwPfhrBuP2BpZMMpslDKcQi4W1WbAqcD/xKRpiURXCEFLYuIxADPABcBTYFrSmFZBgDTVbURMN2bzkZEzgTOAhKBZsBpQLuSDDJEQcvimQw8oapNgNbAnyUUX2GEWhaAIYR2bChxlgDCQFWXquqyYOuJSB2gC/BC5KMqvFDKoaqpqvqz93wXLpnVLon4CiPEz6Q18LuqrlDVA8BbQNfIR1coXYFXvOevAJflsY4CccDhQEWgArCpJIIrpKBl8RJwrKp+BaCqu1V1b4lFGLpQPhdEpBVwLPBlyYRVOJYAStZY4F4gw+c4wsL7aX4q8D+fQymq2sDagOl1lL5kdqyqpnrPN+IOJtmo6g/AN0Cq9/hCVUvjr8ygZQEaA9tFZIpXVfqE90uttAlaFhE5DBgN9C/JwAojKgaFDwcRmQYcl8eigao6NYTXXwz8qapzRaR9mMMLWXHLEbCdKsB7wJ2qujNc8RVGuMrit4LKETihqioiua7bFpGTgCZAHW/WVyJyjqp+F/ZggyhuWXDHpHNwJxZrgLeBG4BJ4Y00uDCU5VbgU1VdVwqbZABLACFT1QuKuYmzgEtFpDPu53q8iLymqtcWP7rQhaEciEgF3MH/dVWdUvyoiiYMZVkPnBAwXcebV6IKKoeIbBKRWqqaKiK1yLs+vBswR1V3e6/5DDgDKPEEEIayrAPmq+oK7zUf4NqaSjwBhKEsZwDniMitQBXgcBHZraql5mIDqwIqIap6v6rWUdX6wNXA1yV98A8H7+qSScBSVX3S73iK6SegkYg0EJHDcZ/Lhz7HlNOHwPXe8+uBvH7ZrAHaiUisl5zbUTovNAilLD8B1UUks+fK84AlJRBbYQUti6r2UtW63v98f2ByaTr4A6Cq9ijmA3cGtg7Yj2t8+8KbfzzuJ2DO9dsDH/sdd1HKAZyNa3RcCMz3Hp39jr2onwnQGVgO/IGrOvI99hzlOAp3lclvwDSghjc/GXjBex4D/Ad30F8CPOl33EUtizfdwft+/QK8DBzud+xFLUvA+jcAT/sdd86HdQVhjDFRyqqAjDEmSlkCMMaYKGUJwBhjopQlAGOMiVKWAIwxJkpZAjBRQ0R2h2EbR3m9oe7O2bujiLQSkV+8nkWz9cgpImNFpG0e22uf2TOs9/zMgGW3ichNxY3ZmPxYAjCmcNKAh8i7f5cJwD+ARt6jE7ikAZyuqsF6hGwPnBkw/SJwezHjNSZflgBMVBORJBGZIyILReT9zH7dReQ0b958r0OyRQCqukdVZ+ESQeB2agHxqjpH3c01k/m7h8juwOcB63YSkV9F5Gdcl9WZHev1Bf7t7fMcdb1grhKR1pF8D0z0sgRgot1k4D5VTcTdefqIN/8l4BZVTQLSQ9hObdydx5kCexY9C5gLICJxwPPAJUArvM7GVHUV8BwwRlWT9O+O3FJwnaMZE3aWAEzUEpFqQHVVnenNegVoKyLVgarqulkGeKOYu6oFbPaenwKsVNXfvF8KrwV57Z+47iuMCTtLAMaEx3r+7o4Zsvcsug/XA2xRxHmvNybsLAGYqKWqO4BtIpJZxdIbmKmq24FdItLGm391CNtKBXaKyOne1T/X8XcPkUuBk7znvwL1ReREb/qagM3sAqrm2HRj3PCWxoSdJQATTSp7g3NnPu7CdeX7hIgsBJKAwd66NwPPi8h84AhgR+ZGRGQV8CRwg7edzHGEb8UN9/k7rnfRz7z5n+Cu8EFV04A+wCdeI3BgP/IfAd0yG4G9eWcBX4Wn+MZkZ72BGpMHEamifw+wMgCopar9irG9WcDF3q+LUF9zKnCXqvYu6n6NKYglAGPyICI9gPtxo+atBm5Q1c0Fv6rA7bUB9qnqwkK8pgPwm3eFkDFhZwnAGGOilLUBGGNMlLIEYIwxUcoSgDHGRClLAMYYE6UsARhjTJT6f+d8c0aeNjijAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "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": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "npdeltap=np.array(pinitarr)-np.array(porbarr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp=npdeltap[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.09412479e-12, -3.91481292e-12, -3.69948516e-12,  8.92952379e-13,\n",
       "        5.56549251e-12])"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope3, intercept3, rval3, pval3, stderr=stats.linregress(np.log10(dtarr),np.log10(np.abs(dp)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.25611759141932516 -11.389452555243532 0.3398523370742992 0.57576689322897\n"
     ]
    }
   ],
   "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": 181,
   "metadata": {},
   "outputs": [],
   "source": [
    "npdeltaL=np.array(deltaLarr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00000000e+00, 0.00000000e+00, 1.25947158e-05],\n",
       "       [0.00000000e+00, 0.00000000e+00, 1.37921244e-05],\n",
       "       [0.00000000e+00, 0.00000000e+00, 5.21423939e-05],\n",
       "       [0.00000000e+00, 0.00000000e+00, 1.29787430e-03],\n",
       "       [0.00000000e+00, 0.00000000e+00, 9.45506016e-02]])"
      ]
     },
     "execution_count": 182,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "npdeltaL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [],
   "source": [
    "dL=npdeltaL[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope4, intercept4, rval4, pval4, stderr=stats.linregress(np.log10(dtarr),np.log10(np.abs(dL)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.2304266856962545,\n",
       " -0.6734525196810321,\n",
       " 0.9320166057156102,\n",
       " 0.021060037537520546)"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope4, intercept4, rval4, pval4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA14klEQVR4nO3deXgUVfbw8e8hAcOqiAtIZHEAZREjRMANkBFEcAXn5wAjoCyi8oLKDIOiM6DgBioqiKOOC+4OrqMwIgICjqgBARFkEYKAyECUnQAh5/3jVkKn6XQ6Sac73X0+z5Mn3bXcOlXdXbfq3qpToqoYY4xJPBWiHYAxxpjosArAGGMSlFUAxhiToKwCMMaYBGUVgDHGJCirAIwxJkFZBWBKREROFZH5IrJHRB4Nc9ljROTVUsz/vYh0DF9EsUFEzhSRpd5nMiza8fgSkZdEZFy04zAFxWUFICK9RSRDRPaKyFYRmSkiF0U7rjgzGNgB1FDVEf4jvR/8Ie8z+FVEPhWRs8IdRKAdi6o2V9V54V5WDBgJzFXV6qr6ZGETedssR0TqRDC2sBKRTO/7dZLf8G9FREWkQZRCC6i8VoBxVwGIyJ3AJOAB4FSgHvA0cHUUwypARJKjHUMY1AdWavA7CR9R1WpAXWAL8M+IRBZl4lTwG1asz7yE35H6wPdFlFsV6AnsAv5UgmVEVKBt6WMD0Mtn2rOBKhEJLF6oatz8AccDe4E/BJnmOFwF8bP3Nwk4zhvXEdgMjAD+B2wFbvTGtQV+AZJ8yroWWO69rgCMAn4EsoC3gRO9cQ0ABQYAPwHzgSTgUdxR9AZgqDdNss+6/NOLYQswLm/ZQH9gITAR+M2b/3KfuE4EXvTW7zfgfZ9xVwBLgZ3Af4GWQbbVBcA3uJ3FN8AF3vCXgMPAIW97Xxpg3peAcT7vuwH7fN6fBrwDbPfiH+Yzbgzwqs/7f3nbfpe37Zp7wwf7xfFvb3gmcKm3jAN5n4M37lxvm1f03t8ErPK20ydA/SDbo523zXYCy4COPuPmAeOBL7xlNvI+z9uAtcAGb7pBwDrgV+BD4DSfMo6ZPkAMV+F28ju9ZTb1hs8BjgDZ3rZoUsj8fYFNwHBghd+4Mbjv7TRgj7ecdJ/xrYBvvXH/At7K+4zxvpN+5SnQyP/7ANQEPvI++9+816nBtmWA9cgE7gG+8Rk2ERjtLbeBz+9omresjd48FXxi/gJ43Nue63Hf+f7eNvof0M9v3zER9xveBjwDVA5h31HY9zR/+wTYRnnljfQp7xrc72gN7vtzd6n3maUtoDz9AV2BHLydaCHT3AcsAk4BTsb9oO/32eg53jQVvY29H6jpjf8R6OxT1r+AUd7r4V65qd4X5R/AG964Bt6HPQ2oClQGhgArvelrArMpWAG855VR1Yv1a+Bmny/uYdzOJAm4BbezF2/8x7gfZ01vPTp4w8/1vkxtvfn64X5IxwXYTififpw3AMm4I63fgFr+X9ZCtrPvl7kq8AqwzHtfAVgM/A2oBJyB+/Fd5rMj8q0AbgKqc7TyXhpoOX47h0u913OAQT7jJgDPeK+vxu2Mm3rreA/w30LWpy6uYu/mxd/Ze3+yN34ebsfQ3Curovd5fupty8pAJ1zl08pbl6eA+T7LKDB9gBiaAPu8ZVfE7RzWAZV8YhhYxG/kM+AR3NlxDtDaZ9wYXAXSzft+PAgs8sZVwu1Ah3vL7oHboZWkAqiFOwup4n2u/6LgQcox2zLAemTiKvnV3ueXhNth1qdgBTAN+MBbTgPcznOAT8w5wI3e/OO85U7xPp8uuMqumjf947hK+0SvvH8DD4a478hf/0DbJ8A2yivvb155g3CV2OvespvjKseGpdpnRmLHHKk/oA/wSxHT/Ah083l/GZDps9EP4FOB4HaY7bzX44AXvNfVcT/G+t77VcDvfearg9tJJ3O0AjjDZ/wcvB269/5Sb5pk3I/zID47AdwOeK7PF3edz7gq3ry1veXm5n3x/NZ9Kl5l5zNsNV4F4Tf8BuBrv2FfAv0L+0L7TfsSbmey04tnA97ZBq4C+slv+ruAF73XY/CpAPymO8Fb1+MLi4OCFcBAYI73WnBHdu299zPxdgbe+wq4H239AMv9K/CK37BP8I4QcTut+/zGK9DJ5/0/cc1iee+red+RBoGmDxDDvcDbfvFuwTsToYgKANccmguk+cT/hM/4McBsn/fNgAPe6/bessRn/EJKUAEEiCsN+M3n/THbMsA8mbjfzD24iqorrvJM9pbbALdTPwQ085nvZmCeT8xrfcad7c17qs+wLC8+wf3ef+cz7nyOntl1JPi+45j1p+gK4ABHz/qre9O39Zl+MXBNsO1U1F+89QFkAScV0X56Gu5IJs9Gb1h+Gaqa4/N+P+6HCq727SEix+GOgJaoal5Z9YH3RGSniOzEVQhHcDvzPJv84thUyLj6uFp/q095/8CdCeT5Je+Fqu73XlYDTgd+VdXfjl116gMj8sr0yj3db/1949voN2wj7kg4VBNV9QTcj/EAcKZPHKf5xXE3BbcVACKSJCIPiciPIrIb98MHOMl/2kK8A5zvdXi2x+0AF/jE8YRPDL/ifuiB1rE+8Ae/mC/CVbh5NgWYz/8zz9+mqroX952tW8j0/vznz/WmD/UzuQFYpapLvfevAb1FpKLPNL/4vN4PpHi/p9OALerteUKItVAiUkVE/iEiG73PdD5wgogklaDsV4DeuJ35NL9xJ+F+R/6/d9/ttc3n9QEAVfUfVg3XWlAFWOzz+f/HG54n2L6jJLJU9YhvbAHiLU35xENnpK8vcUfO1wDTC5nmZwp2ltXzhhVJVVeKyEbgctyX7nWf0ZuAm1T1C//5fK5I8P3xbMU1/+Q53a+sg8BJfl+oUGwCThSRE1R1Z4Bx41V1fAjl5G0nX/VwX/piUdWfRGQ48LKIfOTFsUFVG4cwe29cU82luJ3/8bimKMkrvohl/yYis4DrcU0Fb/rsxPK2x2shxLEJdwYwKNjiihhWYJt6HbK1cEfWwcrwnf9sn/kF973ZUugcBfUF6olI3k4+2Vt+N1wzSTBbgboiIj7b73TcGTW4o+P8DlgRqR2krBG4g4G2qvqLiKTh+hbEZ5qgn2v+RKobRWSDtw4D/EbvwJ1h1cc1t4L7Doe6vfzLOoDrfyrJ/IHWZz8FO61r45qxIiauzgBUdReuzWyKiFzjHWlUFJHLReQRb7I3gHtE5GTvErK/AcW55vx1XDtoe1zbZZ5ngPEiUh/AK//qIOW8DQwXkboicgKuiSFvPbYCs4BHRaSGiFQQkd+JSIeigvPmnQk8LSI1vfVv741+DhgiIm29qyuqikh3EakeoKgZQBPvktpkEbke1yTwUVExFBLXp7gd2GBcf8YeEfmriFT2jvJbiMh5AWatjqsMs3A/lgf8xm/D9SEE8zpu53cdBSvtZ4C7RKQ5gIgcLyJ/KKSMV4ErReQyL94UEekoIqmFTB/IG8CNIpLmnUU+AHylqpkhzv820F1Efu8dtY/AbZv/FjWjiJwP/A5og2vSSANacHTbFOVL3BntUO/7cLVXVp5lQHNv3VJwzUmFqY7bme4UkROBv4ew/GAG4JrO9vkO9I6e38b9Lqt7v807Kd7vPa+sXNzv53EROQXA++1eFmIRgb6nS3FnYEki0hUo8vcdbnFVAQCo6qO4D/keXKfJJtwVNu97k4wDMoDlwHfAEm9YqN7AfVBzVHWHz/AncB1Es0RkD65DuG2Qcp7D7eSX445+ZuA6ffJO+friOt5W4o54p1OwuSGYG3BHPj/g2iFvB1DVDFxn0mSvzHW4U+djqGoW7oqhEbid70jgCr91Lq4JXjnJXtlpuL6BHcDzuKN7f9Nwp+1bcNtikd/4fwLNvNPy9wtZ7odAY1z/0LK8gar6HvAw8KbXFLECd3Z3DFXdhDsTuZuj36u/UIzfkKrOxrXjv4M7ov4d8MdizL8ad+nmU7htdiVwpaoeCmH2fsAHqvqdqv6S94f73l7h7YiDLfsQrtlzAK5f50+4g4GD3vg1uA7Q2birmBYGKW4SrlN8B+7zLPZZpV9sP3rf7UD+H+7sZL0X0+vACyVc1F9xv5lF3vdlNkebNYsS6Hs6HPcZ7sT1X74fcM4yJKohnWmZMiYil+OuTvFvdjGmXBKRr3Df2RejHYspmbg7A4gVXtNHN+90ui7uNPi9aMdlTGFEpIOI1Pa+s/2AlpTy6N1El1UA0SPAWFxTzLe4q4b+FtWIjAnuTFxb/05c0+B1Xp+TiVHWBGSMMQnKzgCMMSZBxdR9ACeddJI2aNAg2mEYY0xMWbx48Q5VPdl/eExVAA0aNCAjo7CrvYwxxgTi3cB6DGsCMsaYBGUVgDHGJCirAIwxJkHFVB9AIIcPH2bz5s1kZ2dHOxRjIiYlJYXU1FQqVqxY9MTGFCIqFYCXcGsMLjtjmyB5PIq0efNmqlevToMGDXDJEY2Jb6pKVlYWmzdvpmHDhtEOx8SwaDUBrcAllppf2oKys7OpVauW7fxNwhARatWqZWe9CWTrnq10eKkDv+z9peiJiyEqFYCqrvIyG4aF7fxNorHvfGK5f/79LPxpIfd9fl9Yyy33fQAiMhiXQ5569epFORpjjImcyuMrk51z9ExvasZUpmZMJSU5hQOjDwSZMzRldgYgIrNFZEWAv2APSTmGqj6rqumqmn7yycfcyFYuJCUlkZaWRvPmzTnnnHN49NFHyc3NDTpPZmYmr7/+etBpjDGJbf2w9fRu0Zsqye7BYVWSq9Dn7D5sGL4hLOWXWQWgqpeqaosAf0U9eq7Mhbs9rXLlyixdupTvv/+eTz/9lJkzZzJ27Nig81gFYIwpSp3qdahxXA2yj2STkpxC9pFsahxXg9rVgj1xM3QJeR9AWbWnAZxyyik8++yzTJ48GVUlMzOTiy++mFatWtGqVSv++1/39L5Ro0axYMEC0tLSePzxxwudzhiT2Lbt28aQ1kNYNGARQ1oPCWtHcFTSQYvItbjH2p2Myy2+VFWLfLZmenq6+ucCWrVqFU2bNg1puf7taXlK255WrVo19u7dW2DYCSecwOrVq6levToVKlQgJSWFtWvX0qtXLzIyMpg3bx4TJ07ko4/cI3b3798fcDpjClOc775JbCKyWFXT/YdHpRPYexZrxJ9+tX7Yev4868+8/8P77M/ZT5XkKlzb9FomdplYZss8fPgwQ4cOZenSpSQlJbFmzZpSTWeMMeFS7q8CCqeybk/Ls379epKSkjjllFMYO3Ysp556KsuWLSM3N5eUlJSA8zz++OMhTWeMMeGScH0AZdmeBrB9+3aGDBnC0KFDERF27dpFnTp1qFChAq+88gpHjhwBoHr16uzZsyd/vsKmM8aYspJQZwAA717/bv7rKd2nhKXMAwcOkJaWxuHDh0lOTuaGG27gzjvvBODWW2+lZ8+eTJs2ja5du1K1alUAWrZsSVJSEueccw79+/cvdDpjjCkrMfVM4NJ2AhsTT+y7b0JVWCdwwjUBGWOMcawCMMaYBGUVgDHGJCirAIwxJkFZBWCMMQnKKgBjjElQVgGEQV466Ly/zMxMLrjgAqB8Zv185plnOPvss0lLS+Oiiy5i5cqVx0yzadMmLrnkEpo1a0bz5s154okn8sfde++9tGzZkrS0NLp06cLPP/8MuEcVDhs2jEaNGtGyZUuWLFlSoMzdu3eTmprK0KFDAZf/qHv37px11lk0b96cUaNGFZj+7bffzl9+796984f/9NNPdOnShaZNm9KsWTMyMzMB+Oyzz2jVqlX+eq1btw6Al156iZNPPjn/83n++eeLLGvOnDm0atWKFi1a0K9fP3JycvLnmTdvXn767w4dOhS5vcaMGUPdunXzlz9jxowC6/nTTz9RrVo1Jk48mpKkQYMG+Z9RevoxV+8ZEx6qGjN/rVu3Vn8rV648ZlikVa1atdBxc+fO1e7du0cwmqLt2rUr//UHH3ygl1122THT/Pzzz7p48WJVVd29e7c2btxYv//++2Pmf+KJJ/Tmm29WVdWPP/5Yu3btqrm5ufrll19qmzZtCpQ5bNgw7dWrl952222qqrpv3z6dM2eOqqoePHhQL7roIp0xY4aqqq5Zs0bT0tL0119/VVXVbdu25ZfToUMHnTVrlqqq7tmzR/ft26eqqo0bN87/PkyZMkX79eunqqovvvhi/jL9BSrryJEjmpqaqqtXr1ZV1XvvvVeff/55VVX97bfftGnTprpx48YCcQXbXn//+991woQJAZevqtqzZ0+97rrrCkxTv3593b59e6HzqJaP776JDUCGBtin2hlAGalWrRpwbNrnYF599VXatGlDWloaN998M0eOHOGbb76hZcuWZGdns2/fPpo3b86KFSuYN28e7du3p3v37px55pkMGTKkyIfQ5KlRo0b+63379gV8vGCdOnVo1aoV4NJWNG3alC1btgSd/4MPPqBv376ICO3atWPnzp1s3boVgMWLF7Nt2za6dOmSP2+VKlW45JJLAKhUqRKtWrVi8+bNADz33HPcdttt1KxZE3BptgFWrlxJTk4OnTt3zt/OVaq4h2WICLt37wZcao3TTjst6HYorKysrCwqVapEkyZNAOjcuTPvvPMOAK+//jo9evTIfzpdXlzBtlcw77//Pg0bNqR58+ZFTmtMuMVVKojbb4elS8NbZloaTJoUfJq8VBAADRs25L33jiY6feihhwqkfS7MqlWreOutt/jiiy+oWLEit956K6+99hp9+/blqquu4p577uHAgQP86U9/okWLFsybN4+vv/6alStXUr9+fbp27cq7777Lddddx/XXX8/q1cc+cvnOO++kb9++AEyZMoXHHnuMQ4cOMWfOnKCxZWZm8u2339K2bdv8YaNHj2batGkcf/zxzJ07F4AtW7Zw+umn50+TmprKli1bOPXUUxkxYgSvvvoqs2fPDriMnTt38u9//5vhw4cD5GdDvfDCCzly5Ahjxoyha9eurFmzhhNOOIEePXqwYcMGLr30Uh566CGSkpJ4/vnn6datG5UrV6ZGjRosWrQov/x33nmH+fPn06RJEx5//HFOP/30Qss66aSTyMnJISMjg/T0dKZPn86mTZvy4zp8+DAdO3Zkz549DB8+PH+bBttekydPZtq0aaSnp/Poo49Ss2ZN9u7dy8MPP8ynn35aoPkHXGXWpUsXRISbb76ZwYMHB/2MjCkJOwMIg7wngi1durTAzr84PvvsMxYvXsx5551HWloan332GevXrwfgb3/7G59++ikZGRmMHDkyf542bdpwxhlnkJSURK9evVi4cCEAb731Vn48vn++O6rbbruNH3/8kYcffphx48YVGtfevXvp2bMnkyZNKnDkP378eDZt2kSfPn2YPHly0HV7+umn6datG6mpqQHH5+Tk0KtXL4YNG8YZZ5yRP2zt2rXMmzePN954g0GDBrFz505ycnJYsGABEydO5JtvvmH9+vW89NJLgMuoOmPGDDZv3syNN96Yn4/pyiuvJDMzk+XLl9O5c2f69euXv4xAZYkIb775JnfccQdt2rShevXqJCUl5c+zePFiPv74Yz755BPuv//+Aqm7A22vW265hR9//JGlS5dSp04dRowYAbi+gTvuuCP/bNHXwoULWbJkCTNnzmTKlCnMnz8/6DY2piTi6gygqCP18kxV6devHw8++OAx47Kysti7dy+HDx8mOzs7P1Gcf9NN3vtQzgDy/PGPf+SWW24JGNPhw4fp2bMnffr0oUePHgGn6dOnD926dWPs2LHUrVs3/0gZYPPmzdStW5cvv/ySBQsW8PTTT7N3714OHTpEtWrVeOihhwAYPHgwjRs35vbbb8+fNzU1lbZt21KxYkUaNmxIkyZNWLt2LampqaSlpeVXFNdccw2LFi3iqquuYtmyZflH3ddffz1du3YFoFatWvnlDhw4ML8SLaysAQMGcP7557NgwQIAZs2alb+TT01NpVatWlStWpWqVavSvn17li1bRpMmTQrdXqeeemr+60GDBnHFFVcA8NVXXzF9+nRGjhzJzp078x8INHToUOrWrQu4JqZrr72Wr7/+mvbt2wf8DIwpsUAdA+X1L5Y6gfOGZWRkaPv27fOHb968WTt16nTM9N9//702atQov1MxKytLMzMzVVX1yiuv1Ndee03HjRuX35k5d+5cTUlJ0fXr1+uRI0e0S5cuOn369JDiXbNmTf7rDz/8UANt19zcXL3hhht0+PDhQed/8skntWfPnqqq+tFHHxXoBD7vvPOOmde/Q3b06NHao0cPPXLkSIHpZs6cqX379lVV1e3bt2tqaqru2LFDc3JytGXLlvq///1PVVX79++vkydP1sOHD2utWrXyO26ff/557dGjh6q6Dto87777rrZt21ZVtdCyVI927mZnZ2unTp30s88+U1X3fevUqZMePnxY9+3bp82bN9fvvvsu6PbyXf5jjz2m119//THT+HYU7927V3fv3p3/+vzzz9eZM2ceM095+O6b2EAhncBxdQZQHvmnfb744otJTj52szdr1oxx48bRpUsXcnNzqVixIlOmTOHzzz+nYsWK9O7dmyNHjnDBBRcwZ84cKlSowHnnncfQoUNZt24dl1xyCddee21IMU2ePJnZs2dTsWJFatasycsvvwzAzz//zMCBA5kxYwZffPEFr7zySv6liAAPPPAA3bp1Y9SoUaxevZoKFSpQv359nnnmGQC6devGjBkzaNSoEVWqVOHFF18MGsfmzZsZP348Z511Vn4H6tChQxk4cCCXXXYZs2bNolmzZiQlJTFhwoT8I/mJEyfy+9//Pu+ggEGDBpGcnMxzzz1Hz549qVChAjVr1uSFF14A4Mknn+TDDz8kOTmZE088Mb/JKCkpKWBZABMmTOCjjz4iNzeXW265hU6dOgHQtGlTunbtSsuWLalQoQIDBw6kRYsWLFy4sNDtNXLkSJYuXYqI0KBBA/7xj38E3S7btm3L/yxzcnLo3bt3/tmMMeFk6aAjbPLkydSrV4+rrrqqVOX4P1PYJJ5Y++6b6ClXzwROZHk3QRljTLRZBRCjOnbsSMeOHaMdhjEmhsXFZaCx1IxlTDjYdz5x5ObC9Ongk40kbGK+AkhJSSErK8t+ECZhqCpZWVmkpKREOxRTxubNg/POgz/8Abyb0cMq5puAUlNT2bx5M9u3b492KMZETEpKSqE31pnYt2YNjBwJH3wA9erBa6+5SiDcYr4CyLtRyBhjYl1WFtx3Hzz9NFSuDA884FLcVK5cNsuL+QrAGGNi3cGDMHkyjBsHu3fDoEEwdiz43EReJqwCMMaYKFF1bft//SusXw+XXw4TJkCkksPGfCewMcbEoq++gosvdm37VarAJ5/AjBmR2/mDVQDGGBNRGzdC797Qrh2sWwfPPuvS2Ps8KiNirAnIGGMiYPduePBBePxxqFAB7rnHXelTvXr0YrIKwBhjylBODjz/PPztb7B9O9xwA4wfDz7PTooaawIyxpgyoOra9M85B265BZo1g4wMmDatfOz8wSoAY4wJu+XL4bLLoHt3OHwY3nsP5s6F1q2jHVlBVgEYY0yYbN0KAwfCuee6o/1Jk2DFCrjmGvB7gF+5YH0AxhhTSvv3w6OPwsMPw6FDMHy46+Q98cRoRxZcVCoAEZkAXAkcAn4EblTVndGIxRhjSio3F159Fe6+G7ZsgZ494aGHoFGjaEcWmmg1AX0KtFDVlsAa4K4oxWGMMSWSl6mzXz847TSYP9+lbY6VnT9EqQJQ1VmqmpfdehFgaQ2NMTFhzRrXpn/JJbBjh8vUuWiRu6s31pSHTuCbgJmFjRSRwSKSISIZlvLZGBMtWVmubb95c5gzx2Xq/OEHd1dvhfKwJy2BMusDEJHZQO0Ao0ar6gfeNKOBHOC1wspR1WeBZ8E9FL4MQjXGmEJFK1NnJJRZBaCqlwYbLyL9gSuA36s9zssYU85EO1NnJETlxEVEugIjgatUdX80YjDGmMKUh0ydkRCtlqvJQHXgUxFZKiLPRCkOY4zJV54ydUZCVO4DUNUYulDKGBPvfDN1ipSPTJ2RYHcCG2MSVnnO1BkJMXrxkjHGlFwsZOqMBKsAjDEJZfly16Zf3jN1RoJVAMaYhJCXqTMtDRYvLv+ZOiPB+gCMMXHNP1Pn7bfHRqbOSLAKwBgTl2I9U2ckWBOQMSbuxEOmzkiwCsAYEzd8M3Vu3+7OAGI1U2ckWAVgjIl5WVkwbFjBTJ2rV0OfPrGbqTMSrA/AGBOz4jlTZyRYBWCMiTmJkKkzEkI6ORKRmiLSXETOEBE7oTLGRE2iZOqMhELPAETkeOA2oBdQCdgOpACnisgi4GlVnRuRKI0xCW/jRrjrLnjjDdfE8+yzcNNNkJQU7chiV7AmoOnANOBiVd3pO0JEWgM3iMgZqvrPMozPGJPgEjVTZyQUWgGoaucg821V1dvDH44xxjiJnqkzEkraCbwIqBfOQIwxBlwH78yZ8Je/wMqV0L69a+NPT492ZPGnpB26CZo6yRhTlgJl6pw3z3b+ZaWkFYA9xN0YEzaWqTM6gl0F9BSBd/QCnFBWARljEodl6oyuYH0AGSUcZ4wxQVmmzvIhWAVQF5ipqt9GKhhjTPybNw9GjIAlS1zGzjfesGRt0RKsD+BHYLiIfCsiL4nI9SJSM1KBGWPii2XqLH+C3QfwFvAWgIicC3QF3hWRJGA28B9V/ToiURpjYlZWlkvQNnUqVK7sMnXefrt7baIrpKuAVPVbVX1QVS8BrgBygYFlGpkxJqYdPOg6eBs1gilTYMAAWLvWpXOwnX/5UOwbwVR1t4gMVlW7EcwYcwzL1Bk77EYwY0zYWKbO2GI3ghljSm3jRujdG9q1g3XrXKbOpUvdXb2m/Ap2I9i/KfxGsFplFpExJmbs2uUydU6aZJk6Y1GwPoCJJRxnjIlzOTnw3HPw979bps5YFuwy0M8jGYgxpvzLy9T55z/DqlXQoYO70qd162hHZkrCHu9ojAmJb6bOnByXqXPu3OLv/Lfu2UqHlzrwy95fyiZQEzKrAIwxQYU7U+f98+9n4U8Lue/z+8IcqSkuUY2dC3rS09M1I8Py0BkTCf6ZOocOLV2mzsrjK5Odk33M8JTkFA6MPlDKaE0wIrJYVY95qkKRZwAicrKITBSRGSIyJ++vbMI0xkRbbi5MmwZNmrjHMXbt6p7M9dhjpUvTvH7Yenq36E2V5CoAVEmuQp+z+7Bh+IYwRW6KK5QmoNeAVUBDYCyQCXxThjEZY6Jk3jyXobNfPzjtNJg/H6ZPD0+a5jrV61DjuBpkH8kmJTmF7CPZ1DiuBrWr1S594aZEQqkAaqnqP4HDqvq5qt4EdCrNQkXkfhFZLiJLRWSWiJxWmvKMMaWzejVcfbXL1LljB7z2Wtlk6ty2bxtDWg9h0YBFDGk9xDqCo6zIPgARWaSq7UTkE+BJ4Gdguqr+rsQLFamhqru918OAZqo6pKj5rA/AmPDyz9R5112WqTMeFdYHEEoyuHEicjwwAngKqAHcXppg8nb+nqpYagljIurgQZg8GcaNg927YdAgVxGcemq0IzORFEoF8Juq7gJ2AZcAiMiFpV2wiIwH+vqWW8h0g4HBAPXqWQJSY0rDMnUaX6H0ATwV4rACRGS2iKwI8Hc1gKqOVtXTcZ3MQwsrR1WfVdV0VU0/+eSTQwjXGBOIb6bOqlUtU6cJngzufOAC4GQRudNnVA0gqaiCVfXSEGN4DZgB/D3E6Y0xxbBxo2vbf+MN18Tz3HNw442QVOSv2MS7YE1AlYBq3jS+uf12A9eVZqEi0lhV13pvrwZ+KE15xphj+WbqrFDBMnWaYxWVDO5zEXlJVTeGebkPiciZuEdLbgSKvALIGBMay9RpQhXS8wAkQMIPVb2qpAtV1Z4lndcYE5hl6jTFVdLnARhjypHly2HECJg9Gxo3dpk6r766ZMnaTOII6XkAIlIZqKeqqyMSlTEmJFu3wr33wgsvQM2arr3/llugUqVoR2ZiQSjJ4K4ElgL/8d6niciHZRyXMSaI/fvh/vvd0f60aXDHHe5ZvMOH287fhC6U+wDGAG2AnQCquhSXGM4YE2G5ufDyywUzda5a5dr6a9aMdnQm1oRSARz27gT2ZakbjImwuXMhPR3693eZOhcscJk6f1firFwm0YVSAXwvIr2BJBFpLCJPAf8t47iMMZ68TJ2dOrnkbXmZOi+6KNqRmVgXSgXw/4DmwEHgDdyNYLeXYUzGGNzOftgwaNHCHf0/8AD88AP07u1u7DKmtIpMBqeq+4HR3p8xpoz5Z+ocPBjGjLFMnSb8QroRLJDS3AhmjDmWZeo0kRbKjWA9gNrAq977XsC2sgzKmETz1VfuRq4vvoCzz3aZOrt0iXZUJt4VeSOYiDzq9ySZf4uIPZbLmDCwTJ0mmkLpSqoqImfkvRGRhrineBljSmjXLhg1Cs48E95/32XqXLsWBg60nb+JnFCeCHYHME9E1gMC1Md7Qpcxpnj8M3X27esydaamRjsyk4hCuQroPyLSGDjLG/SDqh4s27CMiS+WqdOUR4U2AYlI/m0mqnpQVZd5fwe98TVEpEUkgjQmli1f7jp0u3d3ZwDvv++u67edv4m2YGcAPUXkEVwSuMXAdiAFaIR7iHt9YESZR2hMjPLP1PnEEzBkiCVrM+VHsKuA7hCRE4GewB+AOsABYBXwD1VdGJkQjYkt+/a55p1HHoFDh1ymznvusWRtpvwJ2gegqr8Cz3l/xpggcnPhlVdg9GjYsgV69oSHH7Zkbab8CloBiMhlwDVAXW/QFuADVf1PGcdlTEyZO9fdyPXtt9CmDbz5piVrM+VfsFQQk4AmwDRgszc4FRgmIper6vCyD8+Y8m31ahg5Ej78EOrVg9dfh+uvt2RtJjYEOwPopqpN/AeKyFvAGsAqAJOwsrJg7FiYOhUqV3aZOm+/3b02JlYEO07JFpHzAgw/D8guo3iMKdcOHnQdvI0awZQp7s7ddetcOgfb+ZtYE+wMoD8wVUSqc7QJ6HRglzfOmIRhmTpNPAp2GegSoK2I1ManE1hVf4lIZMaUE5ap08SrIruqVPUXVV3s/f0CICJnFTWfMbEuMxN69YJ27eDHH10On2+/tZ2/iR8lvVZhVlijMKYcycvUedZZ8MEH7m5ey9Rp4lGwy0CfLGwUcEKZRGNMFFmmTpNognUC34jL9RMo82evsgnHmMizTJ0mUQWrAL4BVqjqf/1HiMiYMovImAhavtx18M6eDY0bu0ydV10FItGOzJiyF6wP4DpgaaARqtqwTKIxJkK2bnVt+mlpsGSJy9S5YgVcfbXt/E3iCHYZ6K+RDMSYSLBMncYcVeQTwUTkO0D9Bu8CMoBxqppVFoEZE06WqdOYY4XyTOCZwBHgde/9H4EqwC/AS8CVZRKZMWFimTqNCSyUCuBSVW3l8/47EVmiqq1E5E9lFZgxpWWZOo0JLpSfQpKItMl74yWIy7sdJqc0CxeRESKiInJSacoxxldWFgwbBi1auKP/Bx+EH35wd/Xazt+Yo0I5AxgIvCAi1XA3ge0GBohIVeDBki5YRE4HugA/lbQMY3wdPAiTJ8O4cbB7Nwwe7FI2n3JKtCMzpnwqsgJQ1W+As0XkeO/9Lp/Rb5di2Y8DI4EPSlGGMajC9OkufYNl6jQmdEWeEIvI8SLyGPAZ8JmIPJpXGZSUiFyNyyy6rDTlGPPVV65D9//+D6pWhVmzYMYM2/kbE4pQmoBeAFYA/+e9vwF4EegRbCYRmQ3UDjBqNHA3rvmnSCIyGBgMUK9evVBmMQkgM9M9hOXNN6F2bZfD58YbLVmbMcUhqv6X+PtNILJUVdOKGhbyAkXOxp1N7PcGpQI/A22KetZAenq6ZmRklGSxJk7s2uU6dSdNch26f/6zu9KnWrVoR2ZM+SUii1U13X94KGcAB0TkIlVd6BV0IXCgpIGo6ndAfreciGQC6aq6o6RlmvhnmTqNCb9QKoAhwDSfdv/fgH5lF5IxR1mmTmPKTihPBFumqucALYGWqnou0ClcAahqAzv6N4EsW+aevtW9uzsDeP99d12/7fyNCY+Qb4tR1d2qutt7e2cZxWMMW7fCgAFw7rmWqdOYshRKE1Ag9jM0YWeZOo2JrJJWAMEvHTKmGPwzdV53HTz0kGXqNKasBXsm8B4C7+gFqFxmEZmEYpk6jYmeQvsAVLW6qtYI8FddVUt65mAM4DJ1Xn01dOrkkre9/jp8+WVs7vy37tlKh5c68MveoLexGFPuWG5EE1HxmKnz/vn3s/Cnhdz3+X3RDsWYYinyTuDyxO4Ejl0HD8JTT7lMnXv2xEemzsrjK5Odk33M8JTkFA6MLvG9ksaEXWF3AsfoMZeJFarwr39B06bwl7/ABRfA8uUwdWps7/wB1g9bT+8WvamSXAWAKslV6HN2HzYM3xDlyIwJjVUApsz4ZuqsVi3+MnXWqV6HGsfVIPtINinJKWQfyabGcTWoXS1QDkRjyh+rAEzYZWa6Nv127Vx+/ueec1f5dO4c7cjCb9u+bQxpPYRFAxYxpPUQ6wg2McX6AEzYWKZOY8qn0mQDNSYoy9RpTGyyJiBTYqrw8cfQsiXceis0awYZGfDyy7bzNyYWWAVgSiQvU+cVV1imTmNilVUAplgsU6cx8cP6AExI/DN13nmnS95mmTqNiV1WAZigLFOnMfHLmoBMoebOhfR06N8f6taFBQvcXb228zcmPlgFYI6xZk38ZOo0xhTOKgCTb8cOl6mzefP4ydRpjCmc9QGYuMzUaYwpmlUACUwVpk+Hv/4VNmyAbt1gwgR3Q5cxJv7ZiX2CCpSp8+OPbedvTCKxCiDBJFKmTmNMcNYElCD8M3Xee69l6jQm0VkFEOcsU6cxpjDWBBSnLFOnMaYoVgHEIcvUaYwJhVUAccQydRpjisMqgDiwbx/cdx80buwSt915J6xb5+7qrVSp+OVt3bOVDi91sOfbGhPnrAKIYbm5rk3/zDNdJ+/ll8OqVTBxYunSNN8//34W/rSQ+z6/L3zBGmPKnYSoAOLxiLYsMnVWHl8ZGStMzZhKruYyNWMqMlaoPL5y2OI2xpQfCVEBxNMR7erVZZepc/2w9fRu0ZsqyVUAqJJchT5n92HD8A2lL9wYU+7E9X0AlcdXJjsnO//91IypTM2YSkpyCgdGH4hiZMW3Y4dr5586FSpXdjd1DR/uXodLnep1qHFcDbKPZJOSnEL2kWxqHFeD2tVqh28hxphyI67PAOLhiPbgQdem36gRTJkCgwa5Dt5Ro8K788+zbd82hrQewqIBixjSekhcNZsZYwqKyhmAiIwBBgHbvUF3q+qMcC8nlo9oo5Wp893r381/PaX7lLJdmDEmqqJ5BvC4qqZ5f2Hf+eeJxSNay9RpjImEuO4DgNg6os3MhLvugjffhNq1XQ6fG2+EpKRoR2aMiUfRPAMYKiLLReQFESn0qnURGSwiGSKSsX379sImi2m7drk2/bPOgg8+cJk6166FgQNt52+MKTuiqmVTsMhsIFBj+2hgEbADUOB+oI6q3lRUmenp6ZqRkRHWOKMpJweefdbdxLVjh2XqNMaUDRFZrKrp/sPLrAlIVS8NZToReQ74qKziKI9UYcYM+Mtf3J27HTrAo49asjZjTGRFpQlIROr4vL0WWBGNOKLBMnUaY8qLaHUCPyIiabgmoEzg5ijFETFbt8I998CLL7o8PU88AUOGlCxZmzHGhENUKgBVvSEay42Gfftc884jj8ChQy5T5+jRpUvWZowx4RD3l4FGS26uS808ejRs2QLXXQcPPVS6ZG3GGBNOcZ0KIlrKIlOnMcaEm1UAYVSWmTqNMSbcrAIIgx073NO3WrRwR/8PPgg//AC9ekEF28LGmHLK+gBK4eBBeOopGDcO9uyBm2+GMWPglFOiHZkxxhTNKoASiFamTmOMCSdroCgmy9RpjIkXVgGEKDPTtem3awfr17tMnd9+C507RzsyY4wpGWsCKsKuXa5Td9Ik16F7770wcqQ7+jfGmFhmFUAhLFOnMSbeWROQH1XXpt+yJdx2m7u0c/FiePll2/kbY+KLVQA+AmXqnDMHWrWKdmTGGBN+VgHgMnUOGADnngtLlrhMnStWuLt6RaIdnTHGlI2E7gOwTJ3GmESWkBVAXqbOu++Gn3+2TJ3GmMSUcE1Avpk6U1MtU6cxJnElTAVgmTqNMaaghKgAxo2zTJ3GGOMvIfoAGjaEgQNh7FjL1GmMMXkSogLo08f9GWOMOcoaQYwxJkFZBWCMMQnKKgBjjElQVgEYY0yCsgrAGGMSlFUAxhiToKwCMMaYBGUVgDHGJChR1WjHEDIR2Q5sjHYcITgJ2BHtIMIgXtYD4mdd4mU9wNYlkuqr6sn+A2OqAogVIpKhqunRjqO04mU9IH7WJV7WA2xdygNrAjLGmARlFYAxxiQoqwDKxrPRDiBM4mU9IH7WJV7WA2xdos76AIwxJkHZGYAxxiQoqwCMMSZBWQUQBiLyBxH5XkRyRSTopWAikiQi34rIR5GKL1ShrIeInC4ic0VkpTft8EjHGYpQPxMR6Soiq0VknYiMimSMoRCRE0XkUxFZ6/2vWch0j3jru0pEnhQRiXSsRSnGutQTkVneuqwUkQYRDrVIoa6LN20NEdksIpMjGWMorAIIjxVAD2B+CNMOB1aVbTglFsp65AAjVLUZ0A64TUSaRSK4YipyXUQkCZgCXA40A3qVw3UZBXymqo2Bz7z3BYjIBcCFQEugBXAe0CGSQYaoyHXxTAMmqGpToA3wvwjFVxyhrgvA/YS2b4g4qwDCQFVXqerqoqYTkVSgO/B82UdVfKGsh6puVdUl3us9uMqsbiTiK44QP5M2wDpVXa+qh4A3gavLPrpiuRp42Xv9MnBNgGkUSAEqAccBFYFtkQiumIpcF68CTlbVTwFUda+q7o9YhKEL5XNBRFoDpwKzIhNW8VgFEFmTgJFAbpTjCAvv1Pxc4Ksoh1JSdYFNPu83U/4qs1NVdav3+hfczqQAVf0SmAts9f4+UdXyeJZZ5LoATYCdIvKu11Q6wTtTK2+KXBcRqQA8Cvw5koEVR0I8FD4cRGQ2UDvAqNGq+kEI818B/E9VF4tIxzCHF7LSrodPOdWAd4DbVXV3uOIrjnCtS7QFWw/fN6qqInLMddsi0ghoCqR6gz4VkYtVdUHYgy1CadcFt0+6GHdg8RPwFtAf+Gd4Iy1aGNblVmCGqm4uh10ygFUAIVPVS0tZxIXAVSLSDXe6XkNEXlXVP5U+utCFYT0QkYq4nf9rqvpu6aMqmTCsyxbgdJ/3qd6wiAq2HiKyTUTqqOpWEalD4Pbwa4FFqrrXm2cmcD4Q8QogDOuyGViqquu9ed7H9TVFvAIIw7qcD1wsIrcC1YBKIrJXVcvNxQbWBBQhqnqXqqaqagPgj8CcSO/8w8G7uuSfwCpVfSza8ZTSN0BjEWkoIpVwn8uHUY7J34dAP+91PyDQmc1PQAcRSfYq5w6UzwsNQlmXb4ATRCQvc2UnYGUEYiuuItdFVfuoaj3vN/9nYFp52vkDoKr2V8o/3BHYZuAgrvPtE2/4abhTQP/pOwIfRTvukqwHcBGu03E5sNT76xbt2Ev6mQDdgDXAj7imo6jH7rcetXBXmawFZgMnesPTgee910nAP3A7/ZXAY9GOu6Tr4r3v7H2/vgNeAipFO/aSrovP9P2BydGO2//PUkEYY0yCsiYgY4xJUFYBGGNMgrIKwBhjEpRVAMYYk6CsAjDGmARlFYBJGCKyNwxl1PKyoe71z+4oIq1F5Dsvs2iBjJwiMklE2gcor2NeZljv9QU+44aKyE2ljdmYwlgFYEzxZAP3Eji/y1RgENDY++sKrtIA2qlqURkhOwIX+Lx/Afh/pYzXmEJZBWASmoikicgiEVkuIu/l5XUXkfO8YUu9hGQrAFR1n6ouxFUEvuXUAWqo6iJ1N9dM42iGyJ7Af3ym7SoiP4jIElzK6rzEekOAO7xlXqwuC2amiLQpy21gEpdVACbRTQP+qqotcXee/t0b/iJws6qmAUdCKKcu7s7jPL6ZRS8EFgOISArwHHAl0Bov2ZiqZgLPAI+rapoeTeSWgUuOZkzYWQVgEpaIHA+coKqfe4NeBtqLyAlAdXVplgFeL+Wi6gDbvddnARtUda13pvBqEfP+D5e+wpiwswrAmPDYwtF0zFAws+gBXAbYkkjx5jcm7KwCMAlLVXcBv4lIXhPLDcDnqroT2CMibb3hfwyhrK3AbhFp513905ejGSJXAY281z8ADUTkd977Xj7F7AGq+xXdBPd4S2PCzioAk0iqeA/nzvu7E5fKd4KILAfSgPu8aQcAz4nIUqAqsCuvEBHJBB4D+nvl5D1H+Fbc4z7X4bKLzvSGf4y7wgdVzQYGAx97ncC+eeT/DVyb1wnsDbsQ+DQ8q29MQZYN1JgARKSaHn3AyiigjqoOL0V5C4ErvLOLUOc5F7hTVW8o6XKNCcYqAGMCEJHrgbtwT83bCPRX1e3B5wpaXlvggKouL8Y8nYG13hVCxoSdVQDGGJOgrA/AGGMSlFUAxhiToKwCMMaYBGUVgDHGJCirAIwxJkH9fz05wu7fl9ucAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "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": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope2, intercept2, rval2, pval2, stderr=stats.linregress(np.log10(dtarr),np.log10(L0arr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7838044987670053 1.5179134992167318 0.987255815885828 0.001723732858052314\n"
     ]
    }
   ],
   "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": [
    "# Save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "np.savetxt(rundir+runoutprefix+'t.csv',nptimes,delimiter=',')\n",
    "np.savetxt(rundir+runoutprefix+'x1.csv',npstar1x,delimiter=',')\n",
    "np.savetxt(rundir+runoutprefix+'x2.csv',npstar2x,delimiter=',')\n",
    "np.savetxt(rundir+runoutprefix+'x3.csv',npstar3x,delimiter=',')\n",
    "np.savetxt(rundir+runoutprefix+'a1.csv',npstar1a,delimiter=',')\n",
    "np.savetxt(rundir+runoutprefix+'a2.csv',npstar2a,delimiter=',')\n",
    "np.savetxt(rundir+runoutprefix+'a3.csv',npstar3a,delimiter=',')\n",
    "np.savetxt(rundir+runoutprefix+'t.csv',nptimes,delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rundir='NonAdaptiveData/'\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": [
    "rundir='NonAdaptiveData/'\n",
    "runinprefix='Run2/twostarsoneplanet'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "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": "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+'t.csv', newline='') as csvfile:\n",
    "    nptimeslist = list(csv.reader(csvfile))\n",
    "\n",
    "print(nptimeslist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "nptimes=np.zeros(len(nptimeslist))\n",
    "count =0\n",
    "for row in nptimeslist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        nptimes[count]=number\n",
    "    count+=1\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nptimes[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'x1.csv', newline='') as csvfile:\n",
    "    npstar1xlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(npstar1xlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar1x=np.zeros([len(npstar1xlist),len(npstar1xlist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in npstar1xlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        npstar1x[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npstar1x[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'x2.csv', newline='') as csvfile:\n",
    "    npstar2xlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(npstar2xlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar2x=np.zeros([len(npstar2xlist),len(npstar2xlist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in npstar2xlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        npstar2x[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npstar2x[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'x3.csv', newline='') as csvfile:\n",
    "    npstar3xlist = list(csv.reader(csvfile))\n",
    "\n",
    "print(npstar3xlist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar3x=np.zeros([len(npstar3xlist),len(npstar3xlist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in npstar3xlist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        npstar3x[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npstar3x[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'a1.csv', newline='') as csvfile:\n",
    "    npstar1alist = list(csv.reader(csvfile))\n",
    "\n",
    "print(npstar1alist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar1a=np.zeros([len(npstar1alist),len(npstar1alist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in npstar1alist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        npstar1a[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npstar1a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'a2.csv', newline='') as csvfile:\n",
    "    npstar2alist = list(csv.reader(csvfile))\n",
    "\n",
    "print(npstar2alist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar2a=np.zeros([len(npstar2alist),len(npstar2alist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in npstar2alist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        npstar2a[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npstar2a[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open(rundir+runinprefix+'a3.csv', newline='') as csvfile:\n",
    "    npstar3alist = list(csv.reader(csvfile))\n",
    "\n",
    "print(npstar3alist[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npstar3a=np.zeros([len(npstar3alist),len(npstar3alist[0])])\n",
    "rownum=0\n",
    "col=0\n",
    "for row in npstar3alist:\n",
    "    for item in row:\n",
    "        number = float(item)\n",
    "        npstar3a[rownum][col]=number\n",
    "        col+=1\n",
    "    col=0\n",
    "    rownum+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npstar3a[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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}