lab13.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import pi\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "global g, d, m, i, v15, P0, T0, w, zt, cd, m3, S, Tv, Tz, v0\n",
    "g = 9.80665 #ускорение свободного падения\n",
    "d = 12.7 * 10**-3 #диаметр пули\n",
    "m = 48.3 * 10**-3 #масса пули\n",
    "i = 1.0629 #коэффиециен формы\n",
    "v15 = 820.0 #начальная скорость пули T = 15 c\n",
    "P0 = 750 #начальное давление\n",
    "T0 = 15 #начальная температура\n",
    "w = 0.5 #влажность воздуха\n",
    "zt = 1.35 * 10**-3 #коэффиециент температуры заряда\n",
    "cd = 0.0423 #коэффиециент деривации\n",
    "m3 = 0.1744 #коэффиециент убывания угловой скорости\n",
    "S = pi*d*d/4\n",
    "Tv = (T0 + 273.15)/(1 - 3*12.7/(8.0*P0) * w)\n",
    "\n",
    "def T(y):\n",
    "    global Tv\n",
    "    return Tv - 0.0065*y\n",
    "\n",
    "def p(y):\n",
    "    global Tv,P0\n",
    "    return 1.225 * np.power(T(y)/Tv, 4.256) * P0*288.15/(Tv*760.0)\n",
    "\n",
    "def s_speed(y):\n",
    "    return 340.294 * (T(y)/288.15)**(1/2.0)\n",
    "\n",
    "def cx(x):\n",
    "    pa = np.array([0.0525, -0.9476, 8.9342, -9.4610, 0.3207, 4.2980, -1.9382])\n",
    "    pb = np.array([1.0000, -15.4071, 178.6690, -580.8643, 985.5873, -853.9492, 296.9213])\n",
    "    pc = np.array([0.0531, 0.9449, 90.5063, 0.1639])\n",
    "    return np.polyval(pa,x**2)/np.polyval(pb,x**2) + pc[0] / (1+np.exp(-(x-pc[1])*pc[2]))+pc[3]\n",
    "\n",
    "Tz = T0\n",
    "v0 = v15*(1+zt*(Tz-15))\n",
    "\n",
    "def Fs(i,S,p,v,a):\n",
    "    return i*S*p*v*v/2*cx(v/a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "butcher = [np.array([0, 0, 0, 0, 0, 0, 0], dtype=np.longdouble), \n",
    "           np.array([4/7, 0, 0, 0, 0, 0, 0], dtype=np.longdouble), \n",
    "           np.array([115/112, -5/16, 0, 0, 0, 0, 0], dtype=np.longdouble), \n",
    "           np.array([589/630, 5/18, -16/45, 0, 0, 0, 0], dtype=np.longdouble),\n",
    "           np.array([229/1200-29/6000*5**0.5, 119/240-187/1200*5**0.5, -14/75+34/375*5**0.5, -3/100*5**0.5, 0, 0, 0], dtype=np.longdouble),\n",
    "           np.array([71/2400-587/12000*5**0.5, 187/480-391/2400*5**0.5, -38/75+26/375*5**0.5, 27/80-3/400*5**0.5, (1+5**0.5)/4, 0, 0], dtype=np.longdouble),\n",
    "           np.array([-49/480+43/160*5**0.5, -425/96+51/32*5**0.5, 52/15-4/5*5**0.5, -27/16+3/16*5**0.5, 5/4-3/4*5**0.5, 5/2-0.5*5**0.5, 0], dtype=np.longdouble) \n",
    "          ]\n",
    "a = np.array([0, 4/7, 5/7, 6/7, (5-5**0.5)/10, (5+5**0.5)/10, 1],dtype=np.longdouble)\n",
    "b = np.array([1/12, 0, 0, 0, 5/12, 5/12, 1/12], dtype=np.longdouble)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "butcher = [[0, 0, 0, 0], [0.5, 0, 0, 0], [0, 0.5, 0, 0], [0, 0, 1, 0]]\n",
    "a = np.matrix('0; 0.5; 0.5; 1')\n",
    "v = 1/6\n",
    "b = np.matrix([[v], [2*v], [2*v], [v]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mult(vect,list_):\n",
    "    res = np.matrix(np.zeros(vect[0].shape))\n",
    "    for v,l in zip(vect,list_):\n",
    "        res = res + v*l\n",
    "    return res\n",
    "def universal_runge_kutta(f ,u0, t, butcher, a, b):\n",
    "    u = np.array([u0] * len(t))\n",
    "    w = [np.matrix(np.zeros(u0.shape))] * len(butcher)\n",
    "    for i,ui in enumerate(u[:-1]):\n",
    "        tau = t[i+1] - t[i]\n",
    "        for j,wj in enumerate(w):\n",
    "            w[j] = f(u[i] + tau * mult(w, butcher[j]), t[i] + tau * a[j])\n",
    "        u[i+1] = u[i] + tau * mult(w,b)\n",
    "    return u, t\n",
    "def f(ui, x):\n",
    "    #print(ui,'\\n')\n",
    "    global m, g, i, S\n",
    "    u = ui[0].item(0)\n",
    "    gamma = ui[1].item(0)\n",
    "    y = ui[2].item(0)\n",
    "    t = ui[3].item(0)\n",
    "    v = u * np.sqrt(1+gamma**2)\n",
    "    return np.matrix([\n",
    "    [-Fs(i, S, p(y), v, s_speed(y)) / (m*v)],\n",
    "    [-g/u**2],\n",
    "    [gamma],\n",
    "    [1/u]\n",
    "    ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "print('T0 = 15\\n')\n",
    "theta0 = pi/300\n",
    "theta1 = pi/50\n",
    "eps = 10**-15\n",
    "n = 0 \n",
    "while True:\n",
    "    print(\"iteration\",n)\n",
    "    theta = (theta0+theta1)/2\n",
    "    print(theta)\n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta)], #u\n",
    "        [np.tan(theta)],    #gamma\n",
    "        [0.0],              #y\n",
    "        [0.0]               #t\n",
    "        ])\n",
    "    ans = universal_runge_kutta(f, u0, np.linspace(0,2000,16), butcher, a, b)[0][-1][2].item(0)\n",
    "    print(ans)\n",
    "    if abs(ans) < 10**-5:\n",
    "        break\n",
    "    if ans > 0:\n",
    "        theta1 = theta\n",
    "    else:\n",
    "        theta0 = theta\n",
    "    n+=1\n",
    "theta*180/pi/0.06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T0 = 15\n",
      "\n",
      "iteration 0\n",
      "0.019999999999999997\n",
      "iteration 1\n",
      "0.005527825819328758\n",
      "iteration 2\n",
      "4.0500266431584e-06\n",
      "iteration 3\n",
      "4.514686541279289e-10\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "33.93044250040451"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('T0 = 15\\n')\n",
    "theta0 = 0.01\n",
    "theta1 = 0.03\n",
    "eps = 10**-15\n",
    "n = 0 \n",
    "while abs(theta1 - theta0) > eps:\n",
    "    print(\"iteration\",n)\n",
    "    print(abs(theta1 - theta0))\n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta1)], #u\n",
    "        [np.tan(theta1)],   #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0]                  #t\n",
    "        ])\n",
    "    ans1 = universal_runge_kutta(f, u0, np.linspace(0,2000,8), butcher, a, b)\n",
    "    \n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta0)],   #u\n",
    "        [np.tan(theta0)],      #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0]                  #t\n",
    "        ])\n",
    "    ans0 = universal_runge_kutta(f, u0, np.linspace(0,2000,8), butcher, a, b)\n",
    "    new_theta = theta1 - (theta1 - theta0) * ans1[0][-1][2].item(0) / (ans1[0][-1][2].item(0) - ans0[0][-1][2].item(0))\n",
    "    theta0 = theta1\n",
    "    theta1 = new_theta\n",
    "    n+=1\n",
    "theta1*180/pi/0.06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T0 = 5\n",
      "\n",
      "iteration 0\n",
      "0.09000000000000001\n",
      "iteration 1\n",
      "0.06213691282389146\n",
      "iteration 2\n",
      "1.5325823330733512e-05\n",
      "iteration 3\n",
      "4.0579831162013136e-10\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "36.141949438033194"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "global g, d, m, i, v15, P0, T0, w, zt, cd, m3, S, Tv, Tz, v0\n",
    "g = 9.80665 #ускорение свободного падения\n",
    "d = 12.7 * 10**-3 #диаметр пули\n",
    "m = 48.3 * 10**-3 #масса пули\n",
    "i = 1.0629 #коэффиециен формы\n",
    "v15 = 820.0 #начальная скорость пули T = 15 c\n",
    "P0 = 750 #начальное давление\n",
    "T0 = 5 #начальная температура\n",
    "w = 0.5 #влажность воздуха\n",
    "zt = 1.35 * 10**-3 #коэффиециент температуры заряда\n",
    "cd = 0.0423 #коэффиециент деривации\n",
    "m3 = 0.1744 #коэффиециент убывания угловой скорости\n",
    "S = pi*d*d/4\n",
    "Tv = (T0 + 273.15)/(1 - 3*12.7/(8.0*P0) * w)\n",
    "\n",
    "def T(y):\n",
    "    global Tv\n",
    "    return Tv - 0.0065*y\n",
    "\n",
    "def p(y):\n",
    "    global Tv,P0\n",
    "    return 1.225 * np.power(T(y)/Tv, 4.256) * P0*288.15/(Tv*760.0)\n",
    "\n",
    "def s_speed(y):\n",
    "    return 340.294 * (T(y)/288.15)**(1/2.0)\n",
    "\n",
    "def cx(x):\n",
    "    pa = np.array([0.0525, -0.9476, 8.9342, -9.4610, 0.3207, 4.2980, -1.9382])\n",
    "    pb = np.array([1.0000, -15.4071, 178.6690, -580.8643, 985.5873, -853.9492, 296.9213])\n",
    "    pc = np.array([0.0531, 0.9449, 90.5063, 0.1639])\n",
    "    return np.polyval(pa,x**2)/np.polyval(pb,x**2) + pc[0] / (1+np.exp(-(x-pc[1])*pc[2]))+pc[3]\n",
    "\n",
    "Tz = T0\n",
    "v0 = v15*(1+zt*(Tz-15))\n",
    "\n",
    "def Fs(i,S,p,v,a):\n",
    "    return i*S*p*v*v/2*cx(v/a)\n",
    "\n",
    "print('T0 = 5\\n')\n",
    "T0 = 5\n",
    "theta0 = 0.01\n",
    "theta1 = 0.1\n",
    "eps = 10**-15\n",
    "n = 0 \n",
    "while abs(theta1 - theta0) > eps:\n",
    "    print(\"iteration\",n)\n",
    "    print(abs(theta1 - theta0))\n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta1)], #u\n",
    "        [np.tan(theta1)],   #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0]                  #t\n",
    "        ])\n",
    "    ans1 = universal_runge_kutta(f, u0, np.linspace(0,2000,8), butcher, a, b)\n",
    "    \n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta0)],   #u\n",
    "        [np.tan(theta0)],      #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0]                  #t\n",
    "        ])\n",
    "    ans0 = universal_runge_kutta(f, u0, np.linspace(0,2000,8), butcher, a, b)\n",
    "    new_theta = theta1 - (theta1 - theta0) * ans1[0][-1][2].item(0) / (ans1[0][-1][2].item(0) - ans0[0][-1][2].item(0))\n",
    "    theta0 = theta1\n",
    "    theta1 = new_theta\n",
    "    n+=1\n",
    "theta1*180/pi/0.06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 0\n",
      "Iteration: 1\n",
      "Iteration: 2\n",
      "Iteration: 3\n",
      "Theta: 33.93076144422984\n",
      "Alpha: 0.1903695942140304\n"
     ]
    }
   ],
   "source": [
    "global g, d, m, i, v15, P0, T0, w, zt, cd, m3, S, Tv, Tz, v0, V\n",
    "g = 9.80665 #ускорение свободного падения\n",
    "d = 12.7 * 10**-3 #диаметр пули\n",
    "m = 48.3 * 10**-3 #масса пули\n",
    "i = 1.0629 #коэффиециен формы\n",
    "v15 = 820.0 #начальная скорость пули T = 15 c\n",
    "P0 = 750 #начальное давление\n",
    "T0 = 15 #начальная температура\n",
    "w = 0.5 #влажность воздуха\n",
    "zt = 1.35 * 10**-3 #коэффиециент температуры заряда\n",
    "cd = 0.0423 #коэффиециент деривации\n",
    "m3 = 0.1744 #коэффиециент убывания угловой скорости\n",
    "S = pi*d*d/4\n",
    "V = -10\n",
    "Tv = (T0 + 273.15)/(1 - 3*12.7/(8.0*P0) * w)\n",
    "\n",
    "def T(y):\n",
    "    global Tv\n",
    "    return Tv - 0.0065*y\n",
    "\n",
    "def p(y):\n",
    "    global Tv,P0\n",
    "    return 1.225 * np.power(T(y)/Tv, 4.256) * P0*288.15/(Tv*760.0)\n",
    "\n",
    "def s_speed(y):\n",
    "    return 340.294 * (T(y)/288.15)**(1/2.0)\n",
    "\n",
    "def cx(x):\n",
    "    pa = np.array([0.0525, -0.9476, 8.9342, -9.4610, 0.3207, 4.2980, -1.9382])\n",
    "    pb = np.array([1.0000, -15.4071, 178.6690, -580.8643, 985.5873, -853.9492, 296.9213])\n",
    "    pc = np.array([0.0531, 0.9449, 90.5063, 0.1639])\n",
    "    return np.polyval(pa,x**2)/np.polyval(pb,x**2) + pc[0] / (1+np.exp(-(x-pc[1])*pc[2]))+pc[3]\n",
    "\n",
    "Tz = T0\n",
    "v0 = np.sqrt((v15*(1+zt*(Tz-15)))**2 + V**2)\n",
    "\n",
    "def Fs(i,S,p,v,a):\n",
    "    return i*S*p*v*v/2*cx(v/a)\n",
    "\n",
    "def f(ui, x):\n",
    "    global m, g, i, S, cd, m3,V\n",
    "    u = ui[0].item(0)\n",
    "    gamma = ui[1].item(0)\n",
    "    y = ui[2].item(0)\n",
    "    t = ui[3].item(0)\n",
    "    z = ui[4].item(0)\n",
    "    q = ui[5].item(0)\n",
    "    v = np.sqrt((u * np.sqrt(1+gamma**2))**2 + V**2)\n",
    "    return np.matrix([\n",
    "    [-Fs(i, S, p(y), v, s_speed(y)) / (m*v)],\n",
    "    [-g/u**2],\n",
    "    [gamma],\n",
    "    [1/u],\n",
    "    [q*pi*v0*cd+V/u],\n",
    "    [np.exp(-m3*t)/(v**2)/u]\n",
    "    ])\n",
    "\n",
    "theta = 0.01\n",
    "alpha = pi/4\n",
    "\n",
    "def NewtonSolver(f, u0, eps=10**-12, h = 10**-5):\n",
    "    u = u0\n",
    "    dx = np.matrix(np.ones(u0.shape))\n",
    "    i = 0\n",
    "    fu = np.zeros( (u0.shape[0], u0.shape[0]) )\n",
    "    while abs(min(dx).item(0)) > eps :\n",
    "        print(\"Iteration:\", i)\n",
    "        i+=1\n",
    "        du = np.matrix(np.zeros(u0.shape), dtype = np.float64)\n",
    "        for j in range( u0.shape[0] ):\n",
    "            du[j-1] = 0 \n",
    "            du[j] = h\n",
    "            val = ( f(u+du) - f(u-du) )/(2*h)\n",
    "            for k in range( u0.shape[0] ):\n",
    "                fu[j][k]= val[k] \n",
    "        fu = fu.T\n",
    "        dx = -1 * np.linalg.inv(fu) * f(u)\n",
    "        u = u + dx\n",
    "    return u\n",
    "\n",
    "def newton_f(u):\n",
    "    global v0\n",
    "    theta = u[0].item(0)\n",
    "    alpha = u[1].item(0)\n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta)], #u\n",
    "        [np.tan(theta)],   #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0],                  #t\n",
    "        [0.0],                 #z\n",
    "        [np.sin(alpha)]        #q\n",
    "        ])\n",
    "    ans0 = universal_runge_kutta(f, u0, np.linspace(0,2000,8), butcher, a, b)\n",
    "    return np.matrix( [[ans0[0][-1][2].item(0)],[ans0[0][-1][4].item(0)]])\n",
    "\n",
    "res = NewtonSolver(newton_f, np.matrix( [[theta],[alpha]]))\n",
    "res*=180/pi/0.06\n",
    "print('Theta:', res[0].item(0))\n",
    "print('Alpha:', res[1].item(0))\n",
    "res/=180/pi/0.06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.012194517448251816"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20.297385080487224"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta = res[0].item(0)\n",
    "alpha = np.arctan(V/(v15*(1+zt*(Tz-15))))\n",
    "u0 = np.array([          \n",
    "    [v0*np.cos(theta)], #u\n",
    "    [np.tan(theta)],   #gamma\n",
    "    [0.0],                 #y\n",
    "    [0.0],                  #t\n",
    "    [0.0],                 #z\n",
    "    [0.0]                  #q\n",
    "    ])\n",
    "ans0 = universal_runge_kutta(f, u0, np.linspace(0,2000,1024), butcher, a, b)\n",
    "2000*V/(v15*(1+zt*(Tz-15))) - V*ans0[0][-1][3].item(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: overflow encountered in exp\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\numpy\\lib\\polynomial.py:680: RuntimeWarning: overflow encountered in double_scalars\n",
      "  y = y * x + p[i]\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in double_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:23: RuntimeWarning: invalid value encountered in power\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in true_divide\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:49: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:13: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  del sys.path[0]\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:38: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\numpy\\lib\\polynomial.py:680: RuntimeWarning: overflow encountered in cdouble_scalars\n",
      "  y = y * x + p[i]\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\numpy\\lib\\polynomial.py:680: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "  y = y * x + p[i]\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "global g, d, m, i, v15, P0, T0, w, zt, cd, m3, S, Tv, Tz, v0\n",
    "g = 9.80665 #ускорение свободного падения\n",
    "d = 12.7 * 10**-3 #диаметр пули\n",
    "m = 48.3 * 10**-3 #масса пули\n",
    "i = 1.0629 #коэффиециен формы\n",
    "v15 = 820.0 #начальная скорость пули T = 15 c\n",
    "P0 = 750 #начальное давление\n",
    "T0 = 5 #начальная температура\n",
    "w = 0.5 #влажность воздуха\n",
    "zt = 1.35 * 10**-3 #коэффиециент температуры заряда\n",
    "cd = 0.0423 #коэффиециент деривации\n",
    "m3 = 0.1744 #коэффиециент убывания угловой скорости\n",
    "S = pi*d*d/4\n",
    "Tv = (T0 + 273.15)/(1 - 3*12.7/(8.0*P0) * w)\n",
    "\n",
    "def T(y):\n",
    "    global Tv\n",
    "    return Tv - 0.0065*y\n",
    "\n",
    "def p(y):\n",
    "    global Tv,P0\n",
    "    return 1.225 * np.power(T(y)/Tv, 4.256) * P0*288.15/(Tv*760.0)\n",
    "\n",
    "def s_speed(y):\n",
    "    return 340.294 * (T(y)/288.15)**(1/2.0)\n",
    "\n",
    "def cx(x):\n",
    "    pa = np.array([0.0525, -0.9476, 8.9342, -9.4610, 0.3207, 4.2980, -1.9382])\n",
    "    pb = np.array([1.0000, -15.4071, 178.6690, -580.8643, 985.5873, -853.9492, 296.9213])\n",
    "    pc = np.array([0.0531, 0.9449, 90.5063, 0.1639])\n",
    "    return np.polyval(pa,x**2)/np.polyval(pb,x**2) + pc[0] / (1+np.exp(-(x-pc[1])*pc[2]))+pc[3]\n",
    "\n",
    "Tz = T0\n",
    "v0 = v15*(1+zt*(Tz-15))\n",
    "\n",
    "def Fs(i,S,p,v,a):\n",
    "    return i*S*p*v*v/2*cx(v/a)\n",
    "\n",
    "def f(ui, x):\n",
    "    #print(ui,'\\n')\n",
    "    global m, g, i, S\n",
    "    u = ui[0].item(0)\n",
    "    gamma = ui[1].item(0)\n",
    "    y = ui[2].item(0)\n",
    "    t = ui[3].item(0)\n",
    "    v = u * np.sqrt(1+gamma**2)\n",
    "    return np.matrix([\n",
    "    [-Fs(i, S, p(y), v, s_speed(y)) / (m*v)],\n",
    "    [-g/u**2],\n",
    "    [gamma],\n",
    "    [1/u]\n",
    "    ])\n",
    "\n",
    "\n",
    "x = np.linspace(0.0001,8000,256)\n",
    "phi = np.linspace(0.0001,pi/2-0.0001,60)\n",
    "dist = []\n",
    "for theta in phi:\n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta)],    #u\n",
    "        [np.tan(theta)],       #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0]                  #t\n",
    "        ])\n",
    "    ans0 = universal_runge_kutta(f, u0, x, butcher, a, b)\n",
    "    arr = [n[2].item(0)>=0 for n in ans0[0]].index(False)\n",
    "    dist.append(x[arr-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2f4483b6978>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VFX+//HXJx1CGiTUBBIgdAglUkTEChFdsaCLroIVXXX3q65992ddy7q2VVddsWFZyyoKNhARbBQJLRBaQguhpEASEiB1zu+PubCRDWQSJnOnfJ6PRx4z98ydue8JZD5zzz33XDHGoJRSKvAE2R1AKaWUPbQAKKVUgNICoJRSAUoLgFJKBSgtAEopFaC0ACilVIDSAqCUUgFKC4BSSgUoLQBKKRWgQlxZSURigdeAAYABrgE2Ah8CycA24FJjTImICPAPYAJwELjKGLPCep2pwF+sl/2rMWbG8bYbHx9vkpOTm/aOlFIqwC1fvrzYGJPQ2HriylQQIjID+NEY85qIhAGtgfuAfcaYJ0TkHiDOGHO3iEwA/oCzAIwA/mGMGSEibYFMIB1nEVkODDPGlBxru+np6SYzM7PRfEoppf5LRJYbY9IbW6/RLiARiQZOBV4HMMZUG2NKgYnA4W/wM4ALrPsTgbeN0xIgVkQ6AeOBecaYfdaH/jwgo4nvSymllJu4cgygO1AEvCkiK0XkNRGJBDoYY3YDWLftrfW7ADvqPT/fajtWu1JKKRu4UgBCgKHAy8aYIcAB4J7jrC8NtJnjtP/6ySLTRCRTRDKLiopciKeUUqo5XCkA+UC+MWaptfwxzoJQYHXtYN0W1ls/qd7zE4Fdx2n/FWPMq8aYdGNMekJCo8cwlFJKNVOjBcAYswfYISK9raYzgXXAbGCq1TYVmGXdnw1MEaeRQJnVRTQXGCcicSISB4yz2pRSStnApWGgOEf1vGeNANoCXI2zeHwkItcCecAl1rpf4RwBlItzGOjVAMaYfSLyCLDMWu9hY8w+t7wLpZRSTebSMFC76DBQpZRqOleHgbq6B6CUx1TV1rF+dzlr8kspKq9y74uL0D0+kkGJMSS3iyQoqKGxCUoFBi0AylZ1DkNuYQWr80vJyi8lK7+M9bv3U1P33z1TceNndP0d3qiIEAYlxjAoMZY067ZTTATizg0q5cW0ACiPqq518N2GQjK37SMrv4y1u8o4WF0HQJvwEAZ2ieGaU1JIS4wlLSmWzm7+QK6tc5BTWEFWfimr88vIyi9l+g9bqHU4K0NCVDiDulhFIcl52zYyzG3bV8qb6DEA5RE7Sw/x/tI8PliWR3FFNWEhQfTvHE1aYuyRb+Hd4+3pkqmsqWPd7v1k7XDugazOL2VL8YEjewtJbVv9ai9hYJcYIsP1u5PyXnoMQNnO4TD8lFvMO0u2M399AQY4s097fjeyG6f0jCc02Dsmo40IDWZo1ziGdo070lZeWcOanWVkWXsJq/JK+TJrN+DskuqZ0ObIXsKwbnH06xStXUfK5+gegHK70oPVfLw8n3eXbGfb3oO0iwzj0pOSuHx4V5LatrY7XrPtrag6sodwuDAUV1QDMKBLNFeO7Mb5aV1oFRZsc1IV6FzdA9ACoNwmK7+UdxZvZ/bqXVTVOkjvFseVo7qRMaAj4SH+96FojGFXWSXfrS/g3SV5bCwoJzoihEnDkrhiZFe6J7SxO6IKUFoAlEdU1tTx+epdvLtkO6vzy2gdFswFQ7pwxYhu9OscbXc8jzHGsGxbCe8s2c6ctbupqTOc0jOeK0Z246y+7Qnxku4uFRi0AKgWta34AO8t3c5HmfmUHaqhZ/s2XDGiKxcNSyQ6ItTueLYqLK/ko2U7+PfSPHaVVdIxOoLLR3Rl8klJtI+OsDueCgBaAFSLKK6o4rEv1zNz5U5CgoTx/TtyxchujOzeVg+CHqW2zjnk9Z0l2/kxp5iQIOHq0cncelYvHUWkWpSOAlJu5XAYPsrcweNfb+BgdS03ju3B1aOT6aDfaI8pJDiIcf07Mq5/R7YWH+CVhZuZ/uNWvszazUMTB3B2vw52R1QBTvcAVKM2FZRz38w1ZG4vYXhKWx67cAA920fZHcsnZW7bx58/XcvGgnLG9evAg+f3p3NsK7tjKT+jXUDqhB2qruP573KY/sMWoiJCuG9CXyYNS9SunhNUU+fgtR+38o/5mwgW4baze3HVycl6oFi5jRYAdUIWbCzk/llr2bHvEJcMS+TeCX11SgQ327HvIPfPWsuCjUX07xzNYxcOJC0p1u5Yyg9oAVDNUrC/koc/X8eXa3bTIyGSxy4cyIju7eyO5beMMXy9dg8PfZ5NYXkVV47sxh3jewf8SCp1YvQgsGqSOofh3SXbeWruRqrqHPzp7F5MG9vdL0/g8iYiwoSBnRiTGs/T32xixuJtzFm7hwd+058JAztqd5tqUboHoFi7s4w/f7qG1flljEmN55GJA0iOj7Q7VkDKyi/lvk/XsHbnfk7rncAjEwf49PQZyh7aBaQaVVFVyzPfbOKtRVtpGxnG/zuvH+enddZvnTarrXMwY/F2nvlmI3XG8MczU7l+THevmTxPeT8tAOq45q0r4P5Za9ldVsnlI7py9/g+xLTWfmdvsrvsEA/OzmZudgG9OrTh0QsHclJyW7tjKR/gagHQrxQBxuEw/H3uBq5/O5OYVqF88vuTeezCgfrh74U6xbTiX1em89qUdA5U1XHJK4u555MsSg9W2x1N+Qk9CBxADlTVcuuHq5i3roDLhifx0PkDCAvR7wDe7qx+HRjVox3/mJ/D6z9tZd66Av58bl8uHNJFu+vUCdG//gCRX3KQi19exPz1BTz4m348duFA/fD3IZHhzhPxPr/lFLq2a83tH63md68tZUtRhd3RlA/TT4AAsGzbPia++DM7Sw/x1tXDuWp0in5z9FH9OkfzyY0n89cLBrBmZxkTnv+RL7J22R1L+SgtAH7uo8wdXD59CTGtQvns5tGc2ivB7kjqBAUFCVeM7Mb828cyoHMMt/x7Jc/M24TD4b0DOpR30gLgp+ochke/XMddH2cxsns7Pr1pND30ClV+pX10BO9dP4JLhiXy/Pwcbv73Cg5W19odS/kQlwqAiGwTkTUiskpEMq22tiIyT0RyrNs4q11E5HkRyRWRLBEZWu91plrr54jI1JZ5S2p/ZQ3XzljG9B+3ctXJybx51Uk6ysdPhYcE8+SkQfzl3L7Mzd7DpJcXs6v0kN2xlI9oyh7A6caYwfXGlt4DzDfGpALzrWWAc4BU62ca8DI4CwbwADACGA48cLhoKPfZVnyAi15axE85xTx24UAePL+/zjLp50SE68Z05/WpJ7Fj30HOf/Fnlm8vsTuW8gEn8skwEZhh3Z8BXFCv/W3jtASIFZFOwHhgnjFmnzGmBJgHZJzA9tVRFm0u5oKXfqa4oop3rh3B5SO62h1JedDpfdoz86aTiQwP5rJXlzBzRb7dkZSXc7UAGOAbEVkuItOstg7GmN0A1m17q70LsKPec/OttmO1Kzd4d8l2prz+Cwltwpl182hG9dAZPANRaocoPrtpNMO6xXH7R6t54usNenBYHZOrJ4KNNsbsEpH2wDwR2XCcdRsaX2iO0/7rJzsLzDSArl31G2xjauocPPz5Ot5Zsp0z+rTnH5MHE6VTCQe0uMgw3r52OA/OzuaV7zeTW1jOc5OH0EavQ6yO4tIegDFml3VbCHyKsw+/wOrawbottFbPB5LqPT0R2HWc9qO39aoxJt0Yk56QoEMWj6f0YDVT3/iFd5Zs54ZTuzN9Srp++CsAQoODePTCgTwysT8LNhZx8UuL2F2mB4fVrzVaAEQkUkSiDt8HxgFrgdnA4ZE8U4FZ1v3ZwBRrNNBIoMzqIpoLjBOROOvg7zirTTXD9r0HuOCfP5O5rYSnLknj3gl9CQ7Sk7vUr105Kpm3rxnOrtJDXD59KYX7K+2OpLyIK3sAHYCfRGQ18AvwpTFmDvAEcLaI5ABnW8sAXwFbgFxgOnATgDFmH/AIsMz6edhqU020v7KGq99aRtmhGt6fNoJJwxLtjqS82Oie8bx1zUkU7K/kd68tZW9Fld2RlJfQ6aB9jMNhmPZOJgs2FvHedSMYqZdrVC5avHkvV735Cz0S2vD+9SP13BA/ptNB+6nn5ufw7fpC7j+vn374qyYZ1aMdr05JJ7ewgilv/kJ5ZY3dkZTNtAD4kDlr9/D8/BwuGZbIlFHd7I6jfNDYXgn883dDyd5ZxtVvLuNAlU4dEci0APiITQXl/OmjVaQlxfLIBQN0Nk/VbGf368A/Jg9hRV4J17+dSWVNnd2RlE20APiAsoM1THs7k1ZhIfzrimFEhAbbHUn5uHMHdeLpS9NYvGUvN7yznKpaLQKBSAuAlztUXccN72ays/QQr1wxlI4xEXZHUn7iwiGJPH7hQL7fVMQt/15JTZ3D7kjKw7QAeLGq2jqmvZPJ0q37eOqSNNL1guDKzSYP78pD5/dn3roCbvtwFXU6bURA0XPDvVRNnYOb31vJjznFPHnxICYO1mmTVMuYenIyVbV1PPbVBsJCgnhqUhpBelJhQNAC4IVq6xzc+sEqvl1fwMMT+3PpSUmNP0mpEzDt1B5U1jh4Zt4mIkKDeVQHGgQELQBeps5huPPjLL5cs5s/T+jLlFHJdkdSAeIPZ/SksqaOlxZuJjwkiPvP66dFwM9pAfAyj3yxjk9X7uSOcb24/tTudsdRAUREuHN8byprHLzx81bCQ4K5O6O3FgE/pgXAi3y7roC3Fm3jmtEp3HJGqt1xVAASEf7feX2pqq3jle83ExEaxK1n9bI7lmohWgC8xN6KKu6ZmUXfTtHcc04fu+OoACYiPDJxAFW1Dp77NoeI0GBuHNvD7liqBWgB8ALGGO77dA37D9Xy7nVphIXo6Fxlr6Ag4W8XD6Kq1sETX28gPCSIq0en2B1LuZkWAC8wc8VO5mYXcO85fejTMdruOEoBEBwkPHNpGtW1dTz0+TrCQ4L1OtN+Rr9q2mxn6SEenJ3N8OS2XDdGD/oq7xIaHMTzlw3h9N4J/PmzNXyyXC8070+0ANjI4TDc8dFqHMbw9KVpekUv5ZXCQ4J5+YphnNyjHXd+vJrPV//PlVyVj9ICYKM3F21j8Za93P+bfiS1bW13HKWOKSI0mOlT0knv1pZbP1zFN9l77I6k3EALgE1yCsr525wNnNW3PZem65m+yvu1Dgvh9avSGdglhpv/vYIFGwvtjqROkBYAG9TUObj9o9W0CQ/h8YsG6Yk2ymdERYQy4+rh9OoQxe/fXU72rjK7I6kToAXABi/Mz2HNzjIeu3AgCVHhdsdRqkliWofy1tXDiW0VxrS3l7PvQLXdkVQzaQHwsJV5Jfxz4WYuGtqFjAEd7Y6jVLMkRIXzryuHUVRRxc3vraBWryXgk7QAeNCh6jr+9NFqOkSF8+D5/e2Oo9QJSUuK5fELB7J4y14e+2qD3XFUM+iJYB70xNfr2VJ8gH9fN4LoiFC74yh1wi4elsjaXWW88fNW+neO5uJhiXZHUk2gewAe8mNOETMWb+ea0Smc3DPe7jhKuc19E/oyqns77v10DVn5pXbHUU2gBcAD9lfWcOd/suiREMldGb3tjqOUW4UGB/Hi5UNIaBPODe8sp6i8yu5IykUuFwARCRaRlSLyhbWcIiJLRSRHRD4UkTCrPdxazrUeT673Gvda7RtFZLy734y3enbeJgrKK3n60sFEhAbbHUcpt2vXJpxXpwyj5GA1N723nOpaPSjsC5qyB/B/wPp6y38DnjXGpAIlwLVW+7VAiTGmJ/CstR4i0g+YDPQHMoCXRMTvPw3X7drPjEXbuHx4VwYnxdodR6kW079zDH+7eBDLtpXwyBfr7I6jXOBSARCRROBc4DVrWYAzgI+tVWYAF1j3J1rLWI+faa0/EfjAGFNljNkK5ALD3fEmvJXDYbh/1lpiW4dx53jt+lH+b+LgLtxwanfeWbKdD37JszuOaoSrewDPAXcBh/fr2gGlxphaazkf6GLd7wLsALAeL7PWP9LewHP80syVO8ncXsI9GX2IbR1mdxylPOKujD6MSY3n/lnZLN9eYnccdRyNFgAROQ8oNMYsr9/cwKqmkceO95z625smIpkikllUVNRYPK9VdqiGx79az9CusUzSoXEqgAQHCS9cNoSOMRH8/t3lFOyvtDuSOgZX9gBGA+eLyDbgA5xdP88BsSJy+DyCRODwHLH5QBKA9XgMsK9+ewPPOcIY86oxJt0Yk56QkNDkN+Qtnv5mIyUHq3l44gCCdJpnFWBiW4fx6pRhVFTVcuO7y6mqrbM7kmpAowXAGHOvMSbRGJOM8yDud8aY3wELgEnWalOBWdb92dYy1uPfGWOM1T7ZGiWUAqQCv7jtnXiRtTvLeHfJdq4c2Y0BXWLsjqOULfp0jObpS9JYmVfK/Z9l4/wYUN7kRM4DuBu4XURycfbxv261vw60s9pvB+4BMMZkAx8B64A5wM3GGL/7WuBwGP7y2VraRoZx+zg98KsC2zkDO3HL6T35MHMH7y7Vg8LepklTQRhjFgILrftbaGAUjzGmErjkGM9/FHi0qSF9yWerdrJqRylPXZJGTCud7kGp287uxbrd+3lodja9O0QxPKWt3ZGURc8EdqOq2jqe/mYTA7pEc9EQvx7gpJTLgoOEZ387mKS2rbnpveXsLjtkdyRl0QLgRu8uyWNn6SHuzuijB36VqiemVSjTpwyjssbBDe8sp7LG73p/fZIWADfZX1nDi9/lMLpnO8ak+u7oJaVaSs/2UTxzaRpZ+WX8+dO1elDYC2gBcJPpP2yh5GANd2f0sTuKUl5rXP+O3HpWKp+syOc9PShsOy0AblBYXslrP27l3EGdGJSo8/0odTx/PCOVManxPPrlerYWH7A7TkDTAuAGL8zPpabOwR067FOpRgUFCX+flEZYSBC3f7RKLydpIy0AJ2hb8QHe/yWPycOTSImPtDuOUj6hY0wEj1wwgJV5pbzy/Wa74wQsLQAn6KlvNhIaHMQfz0y1O4pSPuX8tM6cN6gTz32bw9qdZXbHCUhaAE7AmvwyvsjazbWnpNA+KsLuOEr5nL9eMIB2bcK47cNVOjTUBloATsArP2wmKiKEaWO72x1FKZ8U2zqMJyelkVNYwVNzN9odJ+BoAWimnaWHmLN2D5cN70p0hE75oFRzje2VwBUju/L6z1tZvHmv3XECihaAZnp78TaMMUwZ1c3uKEr5vPsm9KVb29bc8Z/VlFfW2B0nYGgBaIYDVbW8vzSPjAEdSYxrbXccpXxe67AQnvntYHaXHeLhz/V6wp6iBaAZZq7IZ39lLdeekmJ3FKX8xtCucdx0Wk/+szyfb7L32B0nIGgBaCKHw/Dmz9tIS4xhaNc4u+Mo5Vf+eGYq/TtHc+/MNRRXVNkdx+9pAWii7zcVsaX4ANeckoKIzviplDuFhQTx7G8HU15Vy70z1+iEcS1MC0ATvf7TVjpEhzNhYCe7oyjll3p1iOKu8b2Zt66A/yzPtzuOX9MC0AQb95TzU24xU0YlExqsvzqlWso1o1MYkdKWhz9fx459B+2O47f0U6wJ3vx5KxGhQVw+vKvdUZTya0FBwtOXpgFw9ydZ2hXUQrQAuGhvRRUzV+7koqGJxEWG2R1HKb+XGNeau8/pw6LNe5m9epfdcfySFgAXvf9LHtW1Dq4ZnWx3FKUCxuXDuzIoMYZHv1yvJ4i1AC0ALjDG8J/l+Yzq3o6e7aPsjqNUwAgOEh6ZOICiiiqe+zbH7jh+RwuACzK3l7B970EuHpZodxSlAk5aUiyTT+rKW4u2sWHPfrvj+BUtAC74ZHk+rcOCOWdAR7ujKBWQ7hrfm+iIEO7/LFsPCLuRFoBGVNbU8WXWbs4Z0InI8BC74ygVkOIiw7g7ow+/bNvHzBU77Y7jNxotACISISK/iMhqEckWkYes9hQRWSoiOSLyoYiEWe3h1nKu9Xhyvde612rfKCLjW+pNudPc7D2UV9Vy8bAudkdRKqBdmp7E4KRYHv96PWWH9ICwO7iyB1AFnGGMSQMGAxkiMhL4G/CsMSYVKAGutda/FigxxvQEnrXWQ0T6AZOB/kAG8JKIBLvzzbSET1bspEtsK0amtLM7ilIBLShI+OsFA9h7oJpnvtGLx7hDowXAOFVYi6HWjwHOAD622mcAF1j3J1rLWI+fKc5JcyYCHxhjqowxW4FcYLhb3kUL2VNWyU85RVw0tAtBQTrvj1J2G9AlhitHduPtJdtZukUvHnOiXDoGICLBIrIKKATmAZuBUmNMrbVKPnC4j6QLsAPAerwMaFe/vYHneKVPV+7EYeCioTr6RylvcXdGH5LiWvMnvXjMCXOpABhj6owxg4FEnN/a+za0mnXb0Fdlc5z2XxGRaSKSKSKZRUVFrsRrEcYYPlmRz7BucaTER9qWQyn1a5HhITxzaRq7Sg/x1y/W2x3HpzVpFJAxphRYCIwEYkXk8LCYRODwudr5QBKA9XgMsK9+ewPPqb+NV40x6caY9ISEhKbEc6us/DJyCyu4WL/9K+V10pPbcsPYHnyYuYNv1xXYHcdnuTIKKEFEYq37rYCzgPXAAmCStdpUYJZ1f7a1jPX4d8Y5cHc2MNkaJZQCpAK/uOuNuNsnK/IJCwni3EE67bNS3ujWs1Lp0zGKe2ZmsVcvHtMsruwBdAIWiEgWsAyYZ4z5ArgbuF1EcnH28b9urf860M5qvx24B8AYkw18BKwD5gA3G2Pq3Plm3KWqto7Zq3cxvn9HYlqF2h1HKdWA8JBgnv3tYPYfquW+T/XiMc3R6JlNxpgsYEgD7VtoYBSPMaYSuOQYr/Uo8GjTY3rWgg2FlB6s4eKhXn2MWqmA17dTNLeP68UTX29g5oqdOl1LE+mZwA34ZMVO2keFc0rPeLujKKUacf2Y7pyUHMeDs7Mp2F9pdxyfogXgKBVVtXy/qYhzB3UiRK/6pZTXCw4S/j4pjapaB499paOCmkI/4Y6ycGMh1bUOMvrrxG9K+Yrk+EhuHNudWat2sWhzsd1xfIYWgKPMWbuH+DZhpCe3tTuKUqoJbjq9J0ltW3H/rGxq6hx2x/EJWgDqqaypY8GGQs7u15FgnfpBKZ8SERrMg7/pT25hBW/8tNXuOD5BC0A9P+UUc6C6jgyd918pn3Rm3w6c1bc9/5ifw+6yQ3bH8XpaAOqZk72H6IgQRnXXmT+V8lUP/KY/dQ7DX7/UA8KN0QJgqalzMG9dAWf17UBYiP5alPJVSW1bc/PpPfkyazc/5egB4ePRTzrL0i37KDtUw3jt/lHK5007tTvJ7Vpz/6y1VNV65YQDXkELgGVO9m5ahQYztpd9E9AppdwjIjSYB8/vz5biA7y3JM/uOF5LCwDgcBjmZhdwep8EIkK9/iJlSikXnNa7Paf0jOfFBbl63YBj0AIArMgroai8iowBOvOnUv7k7ow+7DtQzfQfttgdxStpAcB58ldYcBCn99buH6X8ycDEGM4d1InXftpKYbnOE3S0gC8Axhi+XruHU1LjiYrQqZ+V8jd3jOtNda2DF+bn2h3F6wR8AcjetZ+dpYf05C+l/FRKfCSThyfx/i95bCs+YHccrxLwBWDO2j0EBwln9e1gdxSlVAv545mphAYH8dQ3G+2O4lUCvgB8s24PI1La0jYyzO4oSqkW0j4qguvGpPBF1m7W5JfZHcdrBHQBKDlQzaaCCk5J1Qu/KOXvpp3anbjWoTw5d4PdUbxGQBeAVfmlAAxJirM5iVKqpUVFhHLz6T35MaeYHzYV2R3HKwR0AVi5vYQggUGJMXZHUUp5wJWjutGtXWse/Dxbp4gg0AvAjlL6dIwmMjzE7ihKKQ8ID7GmiCg6wOt6zYDALQAOh2FVXilDusbaHUUp5UGn927PuH4deGF+LjtLA/uaAQFbADYXVVBeVcuQrtr/r1Sguf83/TAYHvl8nd1RbBWwBWBlnnUAWPcAlAo4iXGt+cMZqczJ3sPCjYV2x7FNwBaAFXklxLQKJaVdpN1RlFI2uG5MCt3jI3lgdjaVNYF5QLjRAiAiSSKyQETWi0i2iPyf1d5WROaJSI51G2e1i4g8LyK5IpIlIkPrvdZUa/0cEZnacm+rcSvzShmcFEuQXvxdqYAUHhLMQxP7s33vQV4N0NlCXdkDqAX+ZIzpC4wEbhaRfsA9wHxjTCow31oGOAdItX6mAS+Ds2AADwAjgOHAA4eLhqeVV9awqbCcodr/r1RAG5OawLkDO/HPBbns2HfQ7jge12gBMMbsNsassO6XA+uBLsBEYIa12gzgAuv+ROBt47QEiBWRTsB4YJ4xZp8xpgSYB2S49d24KCu/DGO0/18pBX85ry/BQcLjXwfeReSbdAxARJKBIcBSoIMxZjc4iwTQ3lqtC7Cj3tPyrbZjtXvcyrwSANKStAAoFeg6xbTiujHd+WrNHlbvKLU7jke5XABEpA3wCXCrMWb/8VZtoM0cp/3o7UwTkUwRySwqapnTtVfkldKzfRtiWun8/0opuH5MCu0iw3ji6w0Y8z8fS37LpQIgIqE4P/zfM8bMtJoLrK4drNvDY6nygaR6T08Edh2n/VeMMa8aY9KNMekJCe6/QpcxhpV5JQzRb/9KKUtURCi3nNGTxVv28kNOsd1xPMaVUUACvA6sN8Y8U++h2cDhkTxTgVn12qdYo4FGAmVWF9FcYJyIxFkHf8dZbR61fe9BSg7W6AlgSqlfuXxEVxLjWvG3rzfgcATGXoArewCjgSuBM0RklfUzAXgCOFtEcoCzrWWAr4AtQC4wHbgJwBizD3gEWGb9PGy1edTKHc7+/6HddA9AKfVf4SHB3DGuN+t27+fzrP/pnPBLjc6CZoz5iYb77wHObGB9A9x8jNd6A3ijKQHdbWVeKZFhwaS2j7IzhlLKC52f1pl//bCFp7/ZxDkDOhEW4t/nyvr3u2vAirwS0pJiCdYTwJRSRwkKEu7K6E3evoO8/0ue3XFaXEAVgEPVdazfXa7j/5VSx3RarwRGdm/LC9/lUFFVa3ecFhVQBWDNzjLqHEavAKaUOiYR4e6MPhRXVPPaj/49RURAFYDDJ4AN1j0ApdRxDOkaR0b/jkz/YQvFFVX9Qq48AAANfklEQVR2x2kxAVYASunWrjXxbcLtjqKU8nJ3ZvSmstbBi9/l2h2lxQRMATDGsEJPAFNKuahHQhsuTU/kvaXbydvrnxPFBUwBKD1YQ2F5Ff076wXglVKu+b8zexEkwjPzNtodpUUETAEosvrxOsRE2JxEKeUrOsZEcM0pKcxavYt1u443BZpvCpgCUFzuLADxbcJsTqKU8iU3ju1BdEQoT87dYHcUtwuYAnB4DyBBDwArpZogplUoN53Wg4Ubi1i8ea/dcdwqcAqAtQeQEKUFQCnVNFNPTqZTTARPzPGv6aIDpgAUV1QTGix6DQClVJNFhAZz21m9WL2jlLnZe+yO4zYBVACqaBcZjnN2a6WUapqLhnYhtX0bnpy7kdo6h91x3CJgCkBReZV2/yilmi0kOIg7x/dmS9EB/rM83+44bhEwBaC4okpHACmlTsjZ/TowrFscz327iUPVdXbHOWEBVgB0D0Ap1XyHJ4or2F/FW4u22R3nhAVEAXA4DMUV1doFpJQ6YcNT2nJmn/a8tDCX0oPVdsc5IQFRAEoP1VDnMLoHoJRyizszelNRVcvLCzfbHeWEBEQBODyda7zuASil3KBPx2guGpLIm4u2sav0kN1xmi0gCsCRk8B0D0Ap5Sa3nZ0KBp77dpPdUZotIArA4T2AhCgdBaSUco/EuNZMGdWNj5fnk1NQbnecZgmIAlB0ZCI43QNQSrnPzaf3JDIshCfn+uZ00YFRACqqdBoIpZTbxUWGMe3U7sxbV8Ca/DK74zRZQBSA4vJq4tvoNBBKKfe7anQy0REh/HOB7106MjAKgJ4EppRqIVERoVx1cjJzsvewyceOBTRaAETkDREpFJG19draisg8EcmxbuOsdhGR50UkV0SyRGRovedMtdbPEZGpLfN2GlZcofMAKaVaztWjU2gdFuxzewGu7AG8BWQc1XYPMN8YkwrMt5YBzgFSrZ9pwMvgLBjAA8AIYDjwwOGi4QlF5ToPkFKq5cRFhnHFyG58vnoX24oP2B3HZY0WAGPMD8C+o5onAjOs+zOAC+q1v22clgCxItIJGA/MM8bsM8aUAPP436LSIhwOw94D1doFpJRqUdeNSSEkOMinzg5u7jGADsaY3QDWbXurvQuwo956+Vbbsdpb3OFpILQLSCnVktpHRTD5pCRmrsxnp4+cHezug8ANDbMxx2n/3xcQmSYimSKSWVRUdMKB9BwApZSn3DC2B8bAq9/7xl5AcwtAgdW1g3VbaLXnA0n11ksEdh2n/X8YY141xqQbY9ITEhKaGe+/jswDpAVAKdXCusS24uKhiXywbAeF5ZV2x2lUcwvAbODwSJ6pwKx67VOs0UAjgTKri2guME5E4qyDv+Osthan00AopTzp96f1oKbOwes/brU7SqNcGQb6PrAY6C0i+SJyLfAEcLaI5ABnW8sAXwFbgFxgOnATgDFmH/AIsMz6edhqa3H/nQguwhObU0oFuOT4SH6T1pl3lmxn3wHvvl5ASGMrGGMuO8ZDZzawrgFuPsbrvAG80aR0blBUUUVYcBDRrRp9q0op5Ra3nN6Tz1fv4qUFufzlvH52xzkmvz8TuLi8mnZtwnQaCKWUx6R2iOLioYm8vXi7V48I8vsCUKRnASulbHDb2b1A4Nl53nu9AL8vAMXlOg+QUsrzOse2YuqobnyyIp+Ne7xzjiD/LwAVOg2EUsoeN53WkzbhIfx97ga7ozTIrwvA4WkgtAtIKWWHuMgwbhzbg2/XF7Jsm0cGPjaJXxeAkoPV1DmMdgEppWxzzegU2keF88TXG3AOlPQefl0AiiucY3C1ACil7NIqLJhbz+rF8u0lfLu+sPEneJBfF4AjJ4FpF5BSykaXpifSPT6SJ+dsoM7hPXsBfl0AdB4gpZQ3CAkO4o7xvckprOCzlTvtjnNEQBSABC0ASimbnTOgI306RvHSwlyv2Qvw6wJQVK7TQCilvIOIcMsZPdlcdIA5a/fYHQfw9wJgnQOg00AopbzBOQM60T0hkhcX5HrFiCC/LgDFFdXE6wFgpZSXCA4SbjqtJ+t37+e7DfaPCPLrAlBUXqX9/0oprzJxcGcS41rxwnf27wX4dQFwTgOhBUAp5T1Cg4O4cWwPVu0oZdHmvbZm8dsC4HAY9h2oJl6vBKaU8jKThiXSITqcF77LsTWH3xaAw9NAaBeQUsrbRIQGc/2Y7izZso9MG+cI8tsCUHT4JDA9CKyU8kKXj+hK28gwXlyQa1sGvy0AxeU6D5BSynu1Dgvh2lNSWLixiDX5ZbZk8N8CUKHzACmlvNuVo7oRHRHCM/M22rJ9vy0AhyeC0z0ApZS3io4I5fen9WTBxiKWbvH8iCC/LQDFFdY0EBE6DYRSyntdPTqZjtERPDHH89cL8NsCcPhi8DoNhFLKm0WEBnPrWamszCtlbnaBR7ftvwWgXK8FrJTyDZOGJdIjIZK/z91AbZ3DY9v12wJQXFGt/f9KKZ8QEhzEneP7sLnoAJ+syPfYdj1eAEQkQ0Q2ikiuiNzTUtvRaSCUUr5kfP8ODOkay7PzcqisqfPINj1aAEQkGPgncA7QD7hMRPq5ezt1DsNe6xiAUkr5AhHh7ow+7NlfyVuLtnlkm57eAxgO5BpjthhjqoEPgInu3kjJwWocBj0GoJTyKSO7t+P03gm8tCCXsoM1Lb49TxeALsCOesv5VptbFes0EEopH3VXRh/Kq2p56fuWnyLC0wWgoTGZvxr4KiLTRCRTRDKLioqatZHQ4CDOHdiJlPjIZj1fKaXs0rdTNFNHJdM5plWLb0s8eeKBiIwCHjTGjLeW7wUwxjze0Prp6ekmMzPTY/mUUsofiMhyY0x6Y+t5eg9gGZAqIikiEgZMBmZ7OINSSinAo/MkGGNqReQWYC4QDLxhjMn2ZAallFJOHp8oxxjzFfCVp7erlFLq1/z2TGCllFLHpwVAKaUClBYApZQKUFoAlFIqQGkBUEqpAOXRE8GaSkSKgO0n8BLxQLGb4riTt+YCzdZc3prNW3OBZmsuV7J1M8YkNPZCXl0ATpSIZLpyNpyneWsu0GzN5a3ZvDUXaLbmcmc27QJSSqkApQVAKaUClL8XgFftDnAM3poLNFtzeWs2b80Fmq253JbNr48BKKWUOjZ/3wNQSil1DD5fABq7yLyIhIvIh9bjS0Uk2Yuy3S4i60QkS0Tmi0g3b8lWb71JImJExCMjIlzJJSKXWr+3bBH5tydyuZJNRLqKyAIRWWn9m07wUK43RKRQRNYe43ERkeet3FkiMtQTuVzM9jsrU5aILBKRNG/JVm+9k0SkTkQmeVM2ETlNRFZZfwffN2tDxhif/cE5pfRmoDsQBqwG+h21zk3AK9b9ycCHXpTtdKC1df/33pTNWi8K+AFYAqR7Qy4gFVgJxFnL7b3ld4azb/b31v1+wDYPZTsVGAqsPcbjE4CvcV6RbySw1BO5XMx2cr1/y3O8KVu9f/fvcM5gPMlbsgGxwDqgq7XcrL8DX98DcOUi8xOBGdb9j4EzRaShS1N6PJsxZoEx5qC1uARI9EAul7JZHgGeBCq9KNf1wD+NMSUAxphCL8pmgGjrfgywyxPBjDE/APuOs8pE4G3jtASIFZFO3pDNGLPo8L8lnv0bcOX3BvAH4BPAU//PAJeyXQ7MNMbkWes3K5+vFwBXLjJ/ZB1jTC1QBrTzkmz1XYvzW5onNJpNRIYAScaYLzyUyaVcQC+gl4j8LCJLRCTDi7I9CFwhIvk4vzH+wTPRGtXU/4t28eTfQKNEpAtwIfCK3Vka0AuIE5GFIrJcRKY050U8fkEYN2v0IvMurtMSXN6uiFwBpANjWzRRvU020HYkm4gEAc8CV3koz5FNN9B29O8sBGc30Gk4vy3+KCIDjDGlXpDtMuAtY8zT1vWv37GyOVo4W2Ps+htwmYicjrMAnGJ3lnqeA+42xtR5ptOgSUKAYcCZQCtgsYgsMcZsauqL+LJ8IKneciL/u9t9eJ18EQnBuWve2G6fp7IhImcBfwbGGmOqPJDLlWxRwABgofUfvyMwW0TON8Zk2pjr8DpLjDE1wFYR2YizICxrwVyuZrsWyAAwxiwWkQic87Z4tPugAS79X7SLiAwCXgPOMcbstTtPPenAB9bfQDwwQURqjTGf2RsLcP6bFhtjDgAHROQHIA1oUgHwyAGNFjxQEgJsAVL474G5/ketczO/Pgj8kRdlG4LzwGKqt/3ejlp/IZ45COzK7ywDmGHdj8fZtdHOS7J9DVxl3e+L80NWPPRvmsyxDxiey68PAv/i4f9vx8vWFcgFTvZkJleyHbXeW3jwILALv7e+wHzr/2VrYC0woKnb8Ok9AHOMi8yLyMNApjFmNvA6zl3xXJzf/Cd7Uba/A22A/1jfMvKMMed7STaPczHXXGCciKwD6oA7jQe+NbqY7U/AdBG5DWcXy1XG+mttSSLyPs4usXjr+MMDQKiV+xWcxyMm4PygPQhc3dKZmpDtfpzH5F6y/gZqjYcmYXMhm20ay2aMWS8ic4AswAG8Zow57nDWBrfjgf+fSimlvJCvjwJSSinVTFoAlFIqQGkBUEqpAKUFQCmlApQWAKWUClBaAJRSKkBpAVBKqQClBUAppQLU/weLr23BM7jQRwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(phi,dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5803.921596078431"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: overflow encountered in exp\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\numpy\\lib\\polynomial.py:680: RuntimeWarning: overflow encountered in double_scalars\n",
      "  y = y * x + p[i]\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in double_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:23: RuntimeWarning: invalid value encountered in power\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:13: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  del sys.path[0]\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in true_divide\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:49: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:38: RuntimeWarning: invalid value encountered in cdouble_scalars\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5803.921596078431"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOXd9/HPLxthCWFJCBC2IBGJQEZIEYsbKohUiVYEanvX1rY8bbUV3OrWqrVV6wba2sW23rf2uSuLiqSiIFpbV5SgAwQQiOyEfd8Tkuv5I0OfoIRMSDJnZs73/Xrllck1ZzLfi4T5Zs6cOZc55xAREf9J8DqAiIh4QwUgIuJTKgAREZ9SAYiI+JQKQETEp1QAIiI+pQIQEfEpFYCIiE+pAEREfCrJ6wAnk5GR4Xr06OF1DBGRmLJgwYLtzrnMuraL6gLo0aMHxcXFXscQEYkpZrY2nO20C0hExKdUACIiPqUCEBHxKRWAiIhPqQBERHwqrAIwszVmttjMgmZWHBoLmNm8Y2NmNig0bmb2lJmVmtkiMxtQ4/tcZ2YrQx/XNc2UREQkHPU5DHSoc257ja8fAe53zr1uZiNDX18IXAbkhj7OBv4AnG1m7YB7gQLAAQvMrMg5t6vh0xARkfpqyPsAHNA6dDkdKAtdLgSed9VrTc4zszZm1onqcpjrnNsJYGZzgRHACw3IIFHicEUl//vROvYcLPc6SlRonpLEtWd3I715stdRRGoVbgE44A0zc8CfnHPPABOAOWb2GNW7kr4a2jYbWF/jthtCY7WNSxx48LVlPP/hWsy8ThIdnINP1+3iT/81ENM/ikSpcAtgiHOuzMw6AHPN7DNgNDDROfeSmY0B/gpcApzot92dZPw4ZjYeGA/QrVu3MOOJl97+bCvPf7iW752bw88vz/M6TlT4y7ur+NWsZUwrXs/Yr+j3WKJTWC8CO+fKQp+3AjOAQcB1wMuhTaaHxqD6L/uuNW7eherdQ7WNf/G+nnHOFTjnCjIz6zyVhXhsx/4j3PbiIs7omMZtl/b2Ok7UuH5IDkN6tef+fyxl9fYDXscROaE6C8DMWppZ2rHLwHCghOoH7wtCm10ErAxdLgK+HToaaDCwxzm3CZgDDDeztmbWNvR95jTqbCSinHP87KXF7D1UweRxAVKTE72OFDUSEozHrsknOTGBCVODVFRWeR1J5EvC2QWUBcwI7cdMAv7unJttZvuBJ80sCThMaLcN8BowEigFDgLfBXDO7TSzB4D5oe1+eewFYYlNU+av581lW7jna304o2Prum/gM53Sm/PgVf244e+f8Nt/lnLzsNO9jiRyHKs+WCc6FRQUOJ0NNDqt3n6AkU++y8DubXn++kEkJOiFztrcPC3IK59uZPoPv8rA7m29jiM+YGYLnHMFdW0X1aeDFu9VVTlWbN1HVY09GA7HXTNKSElK4LFr8vXgX4f7R53Jx6t3MnFqkKevHUDiF/69crNakZyoN+VL5KkApFbOOX465VNeXbTphNf//psD6JieGuFUsSctNZnJYwOM+dOHXPG79750/aAe7fj7D84mSSUgEaYCkFq99MlGXl20ie8O6cHZOe2Puy6rdTPO6qbdGeEq6NGOWT89j7U7Dh43vnLLPh6fu4I//vtzbrwo16N04lcqADmhdTsOcu/MEgbltOOer+V9abeF1F+fTq3p0+n4F8svPTOLFVv3M/nNlZyXm0l+1zYepRM/0nNO+ZKjlVVMnBYkIcGYNDagB/8mZGb8qrAvHdKaMWFqkIPlR72OJD6iApAv+cO/PmfB2l386sq+ZLdp7nWcuJfeIpnHxuSzZscBfjVrmddxxEdUAHKc4PrdTH5rJYWBzhQGdKqmSPnqaRmMP68nf/9oHXOXbvE6jviECkD+48CRo0yY8ikdW6fyy8K+XsfxnZuHn05ep9b87KVFbN132Os44gN6EdiHKqsct7+4iPlrjn8j9sHySnYcOMILPxis0xh7oFlSIk+OC3D5b9/j0knvkJZa98+gQ1ozfnetDseVU6MC8KE/v7uKlz7ZwCV9skhLPf5X4PzTMxjcs30tt5SmlpuVxh+/NZCihV86T+KXOOeYs2QLt724kOe+q3djS/2pAHymZOMeHn9jOSP7deTpawfoXPVRaOgZHRh6Roewtv3fj9Zy94wS/vuDNXzv3JwmTibxRq8B+Mih8kpumvIp7Vqm8Osr++nBPw5cO6gbl/TpwG9mf8Znm/d6HUdijArARx5+fRmfbzvA49cEaNsyxes40gjMjIev7k/r1CQmTAlyuKLS60gSQ1QAPvH28q08F1q169zcDK/jSCPKaNWMR0b357PN+3hsznKv40gMUQH4wI79R7htulbtimcXnZHFtwZ34y/vreb90u1ex5EYoReB48jhikr+54M17DxQftx48Zqd7D1Uwf/9/iCt2hXH7h6Zxwef7+CWaQuZPeE82rTQbj45ORVAHHl0znL++t5qmn/hQT4pwbi/8Eyt2hXnmqck8uTYs7jq9+9z94wSfnftWXqhX05KBRAn3lu5nb++t5pvn9Nd7+L1sX5d0pk47HQenbOciz7pwNUDu3gdSaKYXgOIA7sOlHPL9CC9OrTizsv6eB1HPPbDC05jUI923Fu0hPU7D9Z9A/EtFUCMc85x14zF7DxQzuSxAZqnaB+/3yUmGE+MzceAiVODHK2sqvM24k8qgBj34oINvF6ymZuH9aZvdrrXcSRKdGnbggeu7Evx2l384V+fex1HopQKIIat23GQ+4qWcHZOO8af39PrOBJlCgOduSK/M5PfWklw/W6v40gU0ovAMaCyylGycQ9Hq9xx4w++toyEBOMJrdolJ3BstbEFa3YycWqQR0f3P+6oIDPo2zmdlCT9HehXKoAYcMu0IK8ET3x2yCfHBbRql9QqvUUyj48JcO1f5jH6jx9+6fpzerbn/37/bP0B4VMqgCg3M7iRV4JlXD8khwt6Zx53XfuWKdrvL3U657T2zJ14Pht3H7/IzOINu3nsjRX8+d1V/PCC0zxKJ15SAUSxDbsOcs8rJQzs3pa7Rp5BUqKeqsup6dUhjV4d0o4bOz83gyVle3n8jeWc2ytDf0z4kB5RolRlleOWaQupqnJMGhPQg780OjPjwav60a5lChOm6kyifqRHlSj1zDur+Gj1Tu4bdSbd2rfwOo7EqbYtU3jsmnxKt+7nodeWeR1HIkwFEIVKNu7hibnVq3aN1lv5pYmdl5vJ9UNyeO7Dtby9fKvXcSSCVABRRqt2iRduH9Gb3llp3P7iInbsP+J1HIkQFUCUeSi0atdj1+Rr1S6JmNTkRCaPC7DnYAV3vLwY51zdN5KYpwKIIm9/tpXnP1zL9UNyOC83s+4biDSiPp1ac/uI3sxduoWp89d7HUciQAUQJbbvP8JtLy6qfho+Qqt2iTeuH5LDkF7tuf8fS1m9/YDXcaSJqQCigHOOO15azN5DFUweF9CqXeKZhATjsWvySUlKYMLUIBU6k2hcUwFEgSnz1/Pmsi3cPqI3fTpp1S7xVqf05jx4VT8Wrt/Nb/9Z6nUcaUIqAI+t2rafX/5jKef2yuD6ITlexxEB4Gv9O3H1gC787p8rWbB2l9dxpImoADxUUVnFxKlBUpISeOyafBJ0Qi6JIveNyiO7bXMmTg2y/8hRr+NIE9C5gDz01FsrWbhhD7//5gA6pqd6HUfkOGmpyUwaE2DMnz7k+v+Zzxkdjz+XUKf05vzgvBydpiSGhVUAZrYG2AdUAkedcwWh8Z8ANwJHgVnOudtD43cC3wtt/1Pn3JzQ+AjgSSAR+Itz7uFGnU0MKV6zk6ffLuXqAV0Y2a+T13FETqigRzvuGtmHP/zrc1Zu2fefcQfsPlhBZVUVN16U611AaZD6PAMY6pzbfuwLMxsKFAL9nXNHzKxDaDwPGAecCXQG3jSz00M3exoYBmwA5ptZkXNuaSPMI6bsO1zBxGlBsts2575ReV7HETmp75/Xk++f9+UV5376wqdMfnMl5+Vmkt+1jQfJpKEa8tztR8DDzrkjAM65YycRKQSmOOeOOOdWA6XAoNBHqXNulXOuHJgS2tZ37itaysZdh5g8NkBaarLXcUROyQOFfemQ1owJU4McLNdrBLEo3AJwwBtmtsDMxofGTgfOM7OPzOzfZvaV0Hg2UPNthBtCY7WN+8qsRZt46ZMN3Di0FwO7t/M6jsgpO7ba2JodB3jgVZ1JNBaFWwBDnHMDgMuAG8zsfKp3H7UFBgO3AdOs+sxlJzqUxZ1k/DhmNt7Mis2seNu2bWHGiw2b9xzmrhmLye/ahp9crP2mEvvOOa0948/vyQsfr2Pu0i1ex5F6CqsAnHNloc9bgRlU787ZALzsqn0MVAEZofGuNW7eBSg7yfgX7+sZ51yBc64gMzN+zodTVeW4ZXqQ8qNVTB4bIFlHTkicuHnY6eR1as3PXlrE1n2H676BRI06H4XMrKWZpR27DAwHSoBXgItC46cDKcB2oAgYZ2bNzCwHyAU+BuYDuWaWY2YpVL9QXNT4U/Leyi37+KB0+3Efj89dzvulO/jFFXnkZLT0OqJIo2mWlMiT4wIcOHKUn724SGcSjSHhHAWUBcwInZc+Cfi7c2526EH8WTMrAcqB61z1T36JmU0DllJ9eOgNzrlKADO7EZhD9WGgzzrnljT6jDz2xpLNjP/bghNeNywvi3Ff6XrC60RiWW5WGneN7MO9RUv427y1fPucHl5HkjBYNLd1QUGBKy4u9jpG2LbuO8yIye/SKT2Vn1+ed9yLHokJRqBrG71pRuKWc47v/Pd85q3awayfnvulReglcsxswbH3a52M3gncSJxz3DZ9EQeOHOXJcQH98ovvmBmPju7PiCff5aYpQWb8eAgpSfqDJ5rpp9NI/jZvLf9esY27v9ZHD/7iWx1ap/LQ1/uxpGwvT8xd4XUcqYMKoBGs3LKPX89axoW9M/mvwd29jiPiqUvP7Mi4r3TlT+98zrxVO7yOIyehAmig8qNV3DQlSMtmSTwyur8WcRcBfn55Ht3bteCWaQvZc6jC6zhSCxVAAz0xdwVLN+3lN1f3p0OazugpAtCyWRKTxgbYvPcwv5hZ4nUcqYUKoAHmrdrBn975nG8M6sawvCyv44hElbO6teWmi3OZGSxjZnCj13HkBFQAp2jPoQpunhqkR/uW3PO1Pl7HEYlKP77wNAZ2b8s9r5Swcfchr+PIF6gATtEvZpawZd8RJo0N0LKZjqYVOZGkxAQmjQlQVeW4eWqQyqrofd+RH6kATsHM4EZmBsuYcHEuAZ0HXeSkurVvwX2jzuSj1Tv587urvI4jNagA6mnDroPc80oJA7u35UcXnuZ1HJGYMHpgFy7r25HH31hOycY9XseREBVAPVRWOW6ZthDnYNKYgE7rIBImM+PBq/rRrmUKE6YGOVRe6XUkQQVQL8+8s4qPVu/kvlFn0q19C6/jiMSUti1TeOyafEq37ufh17WATDRQAYSpZOMenpi7nJH9OnL1AN8tZCbSKM7LzeR75+bw3IdreXv51rpvIE1KBRCGQ+WV3DTlU9q1TOHXV/bTu31FGuC2S3vTOyuN26YvYsf+I17H8TUVQBgeen0Zn287wOPXBGjbMsXrOCIxLTU5kcnjAuw9VMEdLy/WAjIeUgHU4e3PtvL8h2v53rk5nJub4XUckbjQp1Nrbh/Rm7lLtzBl/nqv4/iWCuAktu8/wm0vLuKMjmncdmlvr+OIxJXrh+QwpFd7fvmPpazefsDrOL6kAqiFc447XlrM3sMVTB4XIDU50etIInElIcF47Jp8UpISmDA1SEVlldeRfEcFUIsXPl7Pm8u28LMRZ3BGx9ZexxGJS53Sm/PQ1/uxcP1ufvvWSq/j+I4K4ARWbdvPA68u5dxeGXz3qz28jiMS10b268TVA7rwu7dLWbB2p9dxfEUF8AUVlVVMnBokJSmBx67JJyFBh3yKNLX7RuWR3bY5E6YG2XdYC8hEigrgC556ayULN+zhoa/3o2O6FngRiYS01GQmjQmwcdch7v/HUq/j+IYKoIb5a3by9NuljB7YhZH9OnkdR8RXCnq044ahvXhxwQZeW7zJ6zi+oAII2Xe4golTg3RpW33qWhGJvJ9enEt+l3TufHkxm/cc9jpO3FMBhNxbtISy3YeYNDafVlrgRcQTyYkJTBoboPxoFbdOX0iVFpBpUioA4NVFZbz8yUZuvCiXgd3beR1HxNd6ZrbiF1fk8V7pdp59f7XXceKa7wtg055D3D2jhPyubfjJRb28jiMiwLivdGVYXhaPzF7Osk17vY4Tt3xdAFVVjlunL6T8aBWTxwZI1gIvIlHBzHj46/1o3TyZiVODHK7QAjJNwdePeM++v5r3S3dw7xV55GS09DqOiNTQvlUzHr2mP59t3sejc5Z7HScu+bYAlm3ayyOzlzMsL4uxX+nqdRwROYGhvTvw7XO689f3VvPeyu1ex4k7viyAwxWVTJgSJL1FMg9/XQu8iESzOy/rQ68OrbhlepDdB8u9jhNXfFkAj8xezvIt+3h0dH/at2rmdRwROYnmKYlMHhtg54Fy7pqhBWQak+8K4N2V23j2/dVcd053Luzdwes4IhKGvtnp3DysN68t3sxLn2z0Ok7c8FUB7DpQzq3TF9KrQyvuHNnH6zgiUg/jz+/JoJx23DuzhHU7DnodJy74qgB+PrOEnQfKmTxWC7yIxJrEBOOJMdVn6J04LchRLSDTYL4pgI27D/Hqok38n/NPo292utdxROQUdGnbgl9d2ZcFa3fxh3997nWcmOebAigKlgHokE+RGFcYyGZUfmcmv7WS4PrdXseJab4pgJnBjQzo1oau7Vp4HUVEGuiBK/uSldaMCVM+5cCRo17HiVlhFYCZrTGzxWYWNLPiL1x3q5k5M8sIfW1m9pSZlZrZIjMbUGPb68xsZejjusadSu2Wb97HZ5v3URjIjtRdikgTSm+ezONjAqzdeZBfzVrmdZyYVZ9nAEOdcwHnXMGxATPrCgwD1tXY7jIgN/QxHvhDaNt2wL3A2cAg4F4za9uw+OEpWriRxATTIi8iceSc09oz/vyevPDxOt5YstnrODGpobuAJgG3AzXfmVEIPO+qzQPamFkn4FJgrnNup3NuFzAXGNHA+6+Tc46ZwTKG9MogM01v+hKJJ7cM601ep9bc8fJitu7TAjL1FW4BOOANM1tgZuMBzGwUsNE5t/AL22YD62t8vSE0Vtt4k/pk3S427DrElYHOTX1XIhJhKUkJPPWNAAeOHOW26Yv0LuF6CrcAhjjnBlC9e+cGMzsfuBv4xQm2PdGJddxJxo+/sdl4Mys2s+Jt27aFGa92M4NlNEtKYPiZHRv8vUQk+vTqkMbdX+vDv1ds42/z1nodJ6aEVQDOubLQ563ADOACIAdYaGZrgC7AJ2bWkeq/7Gsea9kFKDvJ+Bfv6xnnXIFzriAzM7PeE6qporKKWYs2cUlelpZ5FIlj/zW4Oxf2zuTXs5ZRunWf13FiRp0FYGYtzSzt2GVgODDfOdfBOdfDOdeD6gf3Ac65zUAR8O3Q0UCDgT3OuU3AHGC4mbUNvfg7PDTWZN4v3c6OA+UU5mv3j0g8MzMeGd2fls2SuGlKkPKjepdwOMJ5BpAFvGdmC4GPgVnOudkn2f41YBVQCvwZ+DGAc24n8AAwP/Txy9BYk5kZLCO9ebJO+ibiAx3SUvnN1f1ZUraXJ+au8DpOTKhzv4hzbhWQX8c2PWpcdsANtWz3LPBs/SKemkPllcxZspnCQGdSknzzfjcRXxuWl8U3BnXjT+98zoW9Mxncs73XkaJa3D4yvrlsCwfLKxmVrzd/ifjJzy/vQ4/2Lbl5apA9hyq8jhPV4rYAZgbL6Ng6lbNz2nkdRUQiqEVKEpPHBtiy7wi/mFnidZyoFpcFsPtgOf9esZVRgc4kJGi5RxG/ye/ahgkX5zIzWMbMoBaQqU1cFoBhTLjkdK4e0MXrKCLikR9deBoF3dtyzyslbNilBWROJC4LIL1FMjcM7UXvjmleRxERjyQlJjBpbADn4JZpC6ms0ruEvyguC0BEBKBruxbcN+pMPlq9k2feWeV1nKijAhCRuHb1gGxG9uvIE3OXU7Jxj9dxoooKQETimpnx4FX9aNcyhZumfMqh8kqvI0UNFYCIxL02LVJ4/JoAn287wEOvawGZY1QAIuIL5+Zm8L1zc3j+w7W8s6LhZxqOByoAEfGN2y7tTef0VJ59f7XXUaKCCkBEfCM1OZHCs7J5d+V2duw/4nUcz6kARMRXCgOdqaxyvLZ4k9dRPKcCEBFfOaNja3pnpTEz+KX1qHxHBSAivjMq0JnitbtYv9Pfp4hQAYiI74wKrRJYtNDfzwJUACLiO13btaCge1uKfL4bSAUgIr5UGOjM8i37+GzzXq+jeEYFICK+NLJfJxITzNcvBqsARMSX2rdqxnm5GRQFy6jy6amiVQAi4ltXBrLZuPsQC9bt8jqKJ1QAIuJbw/KySE1O8O2ykSoAEfGtls2SGJbXkVmLNlFRWeV1nIhTAYiIrxXmd2bXwQreW7nd6ygRpwIQEV87//RM2rRI5r8/WINz/noxWAUgIr6WkpTATy7K5Z0V23jh4/Vex4koFYCI+N53v9qDc3tl8MCrS1m1bb/XcSJGBSAivpeQYDx2TT7NkhOYODXomxeEVQAiIkDH9FQevKofCzfs4am3VnodJyJUACIiISP7dWL0wC48/XYpxWt2eh2nyakARERquPeKPLLbNmfitCD7Dld4HadJqQBERGpIS01m8tgAG3cd4r6ipV7HaVIqABGRLxjYvR03Du3FS59sYNai+F07WAUgInICP7k4l/yubbhrxmI27TnkdZwmoQIQETmB5MQEJo8NUH60ilunL4zLU0arAEREapGT0ZJfXJHH+6U7ePb91V7HaXQqABGRkxj3la4My8vikdnLWbYpvpaPVAGIiJyEmfHw1/uR3iKZCVOCHK6o9DpSowmrAMxsjZktNrOgmRWHxh41s8/MbJGZzTCzNjW2v9PMSs1suZldWmN8RGis1MzuaPzpiIg0vvatmvHo6P4s37KP38z+zOs4jaY+zwCGOucCzrmC0Ndzgb7Ouf7ACuBOADPLA8YBZwIjgN+bWaKZJQJPA5cBecA3QtuKiES9C3t34Jtnd+O5D9awZe9hr+M0ilPeBeSce8M5dzT05TygS+hyITDFOXfEObcaKAUGhT5KnXOrnHPlwJTQtiIiMeG7Q3KocvBqnLw3INwCcMAbZrbAzMaf4PrrgddDl7OBmifV3hAaq238OGY23syKzax427ZtYcYTEWl6vTq0om92a4riZA3hcAtgiHNuANW7b24ws/OPXWFmdwNHgf89NnSC27uTjB8/4NwzzrkC51xBZmZmmPFERCKjMD+bhRv2sHr7Aa+jNFhYBeCcKwt93grMoHp3DmZ2HXA58E33/9dS2wB0rXHzLkDZScZFRGLGFfmdMYOZcfAsoM4CMLOWZpZ27DIwHCgxsxHAz4BRzrmDNW5SBIwzs2ZmlgPkAh8D84FcM8sxsxSqXyguatzpiIg0rY7pqQzOaU9RsCzm1xAO5xlAFvCemS2k+oF8lnNuNvA7IA2YGzo89I8AzrklwDRgKTAbuME5Vxl6wfhGYA6wDJgW2lZEJKYUBjqzavsBSjbG9hvDLJobrKCgwBUXF3sdQ0TkOHsOVlDw67lcd04P7rk8+o5mN7MFNQ7Zr5XeCSwiUk/pLZK5sHcHihaWURnDJ4lTAYiInIIrA9ls3XeEj1bt8DrKKVMBiIicgov7dKBlSiIzg7F7MKMKQETkFKQmJ3Jp3468VrKJI0dj8wRxKgARkVNUGMhm3+Gj/Gt5bJ61QAUgInKKhpzWnoxWKTH7pjAVgIjIKUpKTODy/p15c9lWdh4o9zpOvakAREQa4Nqzu+Gc466XF8fcO4NVACIiDXB6Vhq3Du/N7CWbmb5gg9dx6kUFICLSQN8/ryeDe7bj/qIlrN0RO2cJVQGIiDRQYoLx+JgACQnGxKlBjlZWeR0pLCoAEZFGkN2mOb++qh+frNvN029/7nWcsKgAREQayaj8zlwZ6MxT/1zJp+t2eR2nTioAEZFGdH9hXzq2TmXi1CAHjhyt+wYeUgGIiDSi9ObJPD4mn7U7D/LAq0u9jnNSKgARkUY2uGd7fnjBaUyZv545SzZ7HadWKgARkSYw8ZLT6ZvdmjteWsTWvYe9jnNCKgARkSaQkpTA5LFncaiikltfXBSV7xJWAYiINJFeHVpx98g+vLNiG899sMbrOF+iAhARaULfGtydob0zeej1z1ixZZ/XcY6jAhARaUJmxiOj82nVLImbpgSjavEYFYCISBPLTGvGb67uz7JNe3n+g7Vex/kPFYCISARckpfFWd3a8PKn0bN4jApARCRCCvM7s2zT3qh5LUAFICISIV/r35nEBIuaJSRVACIiEZKZ1owhvTKYGSyLivcFqABERCKoML8zG3Yd4pN1u72OogIQEYmk4Wdm0SwpISp2A6kAREQiKC01mUv6ZDFr0SYqPF45TAUgIhJhowKd2XGgnPdLt3uaQwUgIhJhF/bOpHVqEkXBMk9zqABERCKsWVIiI/t1Ys6SzRwq9+7UECoAEREPjAp05kB5JW8u2+JZBhWAiIgHzs5pT1brZsz0cDeQCkBExAOJCcYV/Tvz7xVb2X2w3JMMKgAREY8UBrKpqHS8XuLNusEqABERj/TNbk3PjJa84tEZQsMqADNbY2aLzSxoZsWhsXZmNtfMVoY+tw2Nm5k9ZWalZrbIzAbU+D7XhbZfaWbXNc2URERig5lRGMjm4zU7Kdt9KOL3X59nAEOdcwHnXEHo6zuAt5xzucBboa8BLgNyQx/jgT9AdWEA9wJnA4OAe4+VhoiIX40KdMY5eHVR5F8MbsguoELgudDl54Ara4w/76rNA9qYWSfgUmCuc26nc24XMBcY0YD7FxGJeTkZLcnvku7J0UDhFoAD3jCzBWY2PjSW5ZzbBBD63CE0ng2sr3HbDaGx2sZFRHxtVCCbJWV7Kd0a2YViwi2AIc65AVTv3rnBzM4/ybZ2gjF3kvHjb2w23syKzax427ZtYcYTEYmpX1NJAAAGHUlEQVRdV/TvRIIR8VNDhFUAzrmy0OetwAyq9+FvCe3aIfR5a2jzDUDXGjfvApSdZPyL9/WMc67AOVeQmZlZv9mIiMSgDq1T+eppGcxcGNmFYuosADNraWZpxy4Dw4ESoAg4diTPdcDM0OUi4Nuho4EGA3tCu4jmAMPNrG3oxd/hoTEREd8bFejM2h0HCa6P3EIx4TwDyALeM7OFwMfALOfcbOBhYJiZrQSGhb4GeA1YBZQCfwZ+DOCc2wk8AMwPffwyNCYi4nsj+nYkJSkhoi8GWzSsS1mbgoICV1xc7HUMEZGI+OHfFlC8dhfz7ryIpMRTP0jTzBbUOGS/VnonsIhIlLjyrM5s33+ED1ftiMj9qQBERKLEhb07kNYsKWK7gVQAIiJRIjU5kRF9OzK7ZDOHK5p+oRgVgIhIFCkMZLP/yFH++dnWujduIBWAiEgUOee09mSmNWNmsOnPEJrU5PcgIiJhS0wwvvPVHhFZK1gFICISZW4Y2isi96NdQCIiPqUCEBHxKRWAiIhPqQBERHxKBSAi4lMqABERn1IBiIj4lApARMSnono9ADPbBqxtwLfIALY3Uhwvxcs8QHOJVvEyl3iZBzRsLt2dc3WuqRvVBdBQZlYczqII0S5e5gGaS7SKl7nEyzwgMnPRLiAREZ9SAYiI+FS8F8AzXgdoJPEyD9BcolW8zCVe5gERmEtcvwYgIiK1i/dnACIiUouYLwAzG2Fmy82s1MzuOMl2o83MmVnUHiFQ11zM7Dtmts3MgqGP73uRMxzh/FzMbIyZLTWzJWb290hnDFcYP5dJNX4mK8xstxc56xLGPLqZ2dtm9qmZLTKzkV7kDEcYc+luZm+F5vEvM+viRc66mNmzZrbVzEpqud7M7KnQPBeZ2YBGDeCci9kPIBH4HOgJpAALgbwTbJcGvAPMAwq8zn2qcwG+A/zO66yNNJdc4FOgbejrDl7nbsjvWI3tfwI863XuU/yZPAP8KHQ5D1jjde4GzGU6cF3o8kXA37zOXctczgcGACW1XD8SeB0wYDDwUWPef6w/AxgElDrnVjnnyoEpQOEJtnsAeAQ4HMlw9RTuXGJBOHP5AfC0c24XgHOu6VfAPjX1/bl8A3ghIsnqJ5x5OKB16HI6UBbBfPURzlzygLdCl98+wfVRwTn3DrDzJJsUAs+7avOANmbWqbHuP9YLIBtYX+PrDaGx/zCzs4CuzrlXIxnsFNQ5l5CrQ08FXzSzrpGJVm/hzOV04HQze9/M5pnZiIilq59wfy6YWXcgB/hnBHLVVzjzuA/4lpltAF6j+tlMNApnLguBq0OXrwLSzKx9BLI1trB//05FrBeAnWDsP4c1mVkCMAm4JWKJTt1J5xLyD6CHc64/8CbwXJOnOjXhzCWJ6t1AF1L9V/NfzKxNE+c6FeHM5ZhxwIvOuaZfzbv+wpnHN4D/cc51oXrXw99C/4eiTThzuRW4wMw+BS4ANgJHmzpYE6jP71+9ReMPtz42ADX/Cu7C8U9b04C+wL/MbA3V+9CKovSF4LrmgnNuh3PuSOjLPwMDI5StvuqcS2ibmc65CufcamA51YUQbcKZyzHjiM7dPxDePL4HTANwzn0IpFJ9PppoE87/lTLn3Nedc2cBd4fG9kQuYqOpz+9fvcV6AcwHcs0sx8xSqP4PWHTsSufcHudchnOuh3OuB9UvAo9yzhV7E/ekTjoXgC/s+xsFLItgvvqocy7AK8BQADPLoHqX0KqIpgxPOHPBzHoDbYEPI5wvXOHMYx1wMYCZ9aG6ALZFNGV4wvm/klHj2cudwLMRzthYioBvh44GGgzscc5taqxvntRY38gLzrmjZnYjMIfqIwOedc4tMbNfAsXOuS/9R41WYc7lp2Y2iuqnsjupPioo6oQ5lznAcDNbClQCtznndniX+sTq8Tv2DWCKCx26EW3CnMctwJ/NbCLVuxm+E43zCXMuFwIPmZmj+gjAGzwLfBJm9gLVWTNCr73cCyQDOOf+SPVrMSOBUuAg8N1Gvf8o/PmKiEgExPouIBEROUUqABERn1IBiIj4lApARMSnVAAiIj6lAhAR8SkVgIiIT6kARER86v8BbkcFzzsu5w0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0.0001,8000,256)\n",
    "phi = np.linspace(0.4,1,60)\n",
    "dist = []\n",
    "for theta in phi:\n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta)],    #u\n",
    "        [np.tan(theta)],       #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0]                  #t\n",
    "        ])\n",
    "    ans0 = universal_runge_kutta(f, u0, x, butcher, a, b)\n",
    "    arr = [n[2].item(0)>=0 for n in ans0[0]].index(False)\n",
    "    dist.append(x[arr-1])\n",
    "plt.figure()\n",
    "plt.plot(phi,dist)\n",
    "max(dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: overflow encountered in exp\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\numpy\\lib\\polynomial.py:680: RuntimeWarning: overflow encountered in double_scalars\n",
      "  y = y * x + p[i]\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in double_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:23: RuntimeWarning: invalid value encountered in power\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:13: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  del sys.path[0]\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in true_divide\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:49: RuntimeWarning: invalid value encountered in cdouble_scalars\n",
      "D:\\Programs\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:38: RuntimeWarning: invalid value encountered in cdouble_scalars\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5823.874782583171"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt0VOd57/HvozuIm6QRBAPiIoFtbGMMAhtNHAfXITg9rdM0dXBOY7tNSuvaWSs5p6dJmnXaxGlP0zSXk5y6aWnq1EmWb01MS2viS5vGiSVsEDYYg7GRhAwC20LiJhmQkPScP2YLDyBZo+ue0f591tLSnndf5nm1YZ7Z77v3+5q7IyIi0ZQVdgAiIhIeJQERkQhTEhARiTAlARGRCFMSEBGJMCUBEZEIUxIQEYkwJQERkQhTEhARibCcsAMYSCwW83nz5oUdhohIxti+fXuLu5emsm3aJ4F58+ZRW1sbdhgiIhnDzF5PdVs1B4mIRJiSgIhIhCkJiIhEmJKAiEiEKQmIiESYkoCISIQpCYiIRFjaPycgo2Nb41F++dqRsMMQ4LryEqrKY2GHIRGlJBBRf/LYLvY1t2MWdiTR5g4/eeEQz35uNaaTISFIKQmYWSPQBnQDXe5eaWZLgb8DCoAu4A/dfauZ/Xfgc8Gu7cBd7r4zOM5a4NtANvA9d//qSFZGUtPcdoZ9ze18/ubL+IMbysMOJ9IeqGnkzzbt5uDR05SVTAw7HImgwfQJrHb3pe5eGbz+GvBld18K/GnwGmA/cIO7LwG+AmwAMLNs4D7gZmAxcJuZLR6BOsggbalvBSCuJojQxStKAKiubwk5Eomq4XQMOzAlWJ4KHAZw9xp3PxaUPwfMDpZXAnXu3uDuncDDwC3DeH8Zouq6FqZOyGXxJVMG3lhGVXnpJKZPzqe6TklAwpFqn4ADT5mZA3/v7huAzwBPmtnXSSSTqj72+yTw02B5FnAwaV0TcO2QopYhc3eq61q5bkEx2Vlqgw6bmRGviPGL147Q0+Nk6ZzIGEv1SiDu7stINOXcbWbvA+4CPuvuc4DPAv+YvIOZrSaRBHr7B/r61+19vZmZrTezWjOrPXJEd7CMpANHT3Ho+GniFWoKSheryktofbuTV99qCzsUiaCUkoC79zb1NAMbSTTt3AE8Fmzyz0EZAGa2BPgecIu7twbFTcCcpMPOJmhC6uP9Nrh7pbtXlpamNCS2pKgm6A/QLYnpozch954bkbE0YBIws0Izm9y7DKwBXibxAX5DsNmNwL5gmzISyeET7v5a0qG2AQvNbL6Z5QHrgE0jVRFJTXVdCzOm5FNeWhh2KBKYNW0C80omUqN+AQlBKn0CM4CNwT3MOcCD7v6EmbUD3zazHOAMsD7Y/k+BEuBvg326gm/1XWZ2D/AkiVtE73f33SNbHXk3PT3OlvpWblhUqnvS00xVRYxNOw7T1d1DTrYe5JexM2AScPcG4Oo+yp8FlvdR/ingU/0cazOwefBhykh49a02Wt/uZFV5SdihyAXi5TEefP4AO5tOsHxuUdjhSIToK0eE9N6GqE7h9NObmNUkJGNNSSBCttS3Mj9WyCXTJoQdilyguDCPxTOnqHNYxpySQER0dffw/P6jVKkpKG1VlZew/cAxzpztDjsUiRAlgYjY2XSC9o4uNQWlsXhFjM6uHmobjw28scgIURKIiN625usW6EogXa2cX0xOlmkcIRlTSgIRUV3fwuKZUyguzAs7FOlHYX4OS+dMU+ewjCklgQg4c7abF14/fm7ESklfVRUxdh06wYnTZ8MORSJCSSACahuP0dndQ5X6A9JevLyEHofnG3SXkIwNJYEIqK5vISfLWDmvOOxQZABLy6ZRkJulW0VlzGh6yTRTU9fC3/68nh7vc4DVIXnljZMsnTONwnyd7nSXn5PNinnFPPZCE6+N4KiiZvCp6xew+tLpI3ZMGR/0qZBm/qmmkR0Hj3P5zMkjdsyK6ZP45Hvnj9jxZHT9bnw+3/15PWe7e0bsmHvfbOP+Z/crCchFlATSSHeP81xDK7961Uz+6qNLwg5HQrL6sumsvmxkP6y/tGk3D287QEdXN/k52SN6bMls6hNIIy8fOsHJM11U6S4eGWHxihhnzvbwwuvHww5F0oySQBrpfUhIE77ISLt2QTFZBjV6EE0uoCSQRrbUt3LpjMmUTs4POxQZZ6YU5LJk9jTddSQXURJIEx1d3WxrPKqmIBk18YoSdh48TntHV9ihSBpREkgTL7x+nDNne9QUJKOmqjxGV4+zdb+uBuQdKSUBM2s0s11mtsPMaoOypWb2XG+Zma0Mys3MvmNmdWb2kpktSzrOHWa2L/i5Y3SqlJlq6lvIskTbrchoWD63iLycLKrrlATkHYO5RXS1uyf3Kn0N+LK7/9TMPhS8fj9wM7Aw+LkW+C5wrZkVA38GVAIObDezTe6ucXNJzPq1ZPY0phTkhh2KjFMFudlUzi06N8OcCAyvOciBKcHyVOBwsHwL8ANPeA6YZmYzgQ8CT7v70eCD/2lg7TDef9xo7+hiZ9MJDfAmoy5eEWPvm220tneEHYqkiVSTgANPmdl2M1sflH0G+GszOwh8HfhCUD4LOJi0b1NQ1l/5RcxsfdDEVHvkyJEUQ8xcW/e30t3jxNUfIKOsd2a5LRqgTgKpJoG4uy8j0dRzt5m9D7gL+Ky7zwE+C/xjsK31sb+/S/nFhe4b3L3S3StLS0tTDDFzVde1kpeTxbK5RWGHIuPcVbOmMjk/R/0Cck5KScDdDwe/m4GNwErgDuCxYJN/Dsog8Q1/TtLus0k0FfVXHnnVdS1Uzi2iIFeP88voysnO4toFxXpoTM4ZMAmYWaGZTe5dBtYAL5P4AL8h2OxGYF+wvAm4PbhL6DrghLu/ATwJrDGzIjMrCo7z5IjWJgO1tHew9802zf0rY6aqPMbrradoOnYq7FAkDaRyd9AMYKOZ9W7/oLs/YWbtwLfNLAc4A/T2FWwGPgTUAaeA3wFw96Nm9hVgW7Ddve5+dMRqkqGeC9pme9tqRUZb7xeOmvpWbq2cGHI0ErYBk4C7NwBX91H+LLC8j3IH7u7nWPcD9w8+zPGruq6Vyfk5XDVratihSEQsmjGJ2KQ8aupauLVyzsA7yLimJ4ZDVlPfwrULisnJ1qmQsWFmrCqPUV3fio/g5EWSmfTJE6KmY6d4vfWUhoqQMRcvL+FIWwd1ze1hhyIhUxIIUU1wm546hWWs9f6b09PDopnFhqH55Bk+9J1nOXn67JD27+rpITYpj0UzJo1wZCLvbk7xROYUT+Def9/D/9m8d0jHKMzPZtM972VOsTqXM5mSwDD8/NUjtLR3cPuquUOexH3l/GKCO69ExtRXP7KEZ4d4JdBxtof7q/fzs73N3FE1b2QDkzGlJDAMNfUtxCbl8eVfv0If5JJx4hWxYTVFPrXnTWrqW5QEMpz6BIbI3amub2VVeUwJQCIpXh5jS31i3CvJXEoCQ1TX3M6Rtg7ieshLIqqqooSTZ7rYffhE2KHIMCgJDFHvXRW6s0eialXwBUiD0WU2JYEhqq5vZXbRBN0ZIZE1fXIBi2ZM0mB0GU5JYAi6e5znGlo1/r9EXlV5jG2NR+no6g47FBkiJYEhePnQCdrOdFGlmcAk4uIVMc6c7eHFA8fDDkWGSElgCKqDy18N9yBRd+2CYrIMavTkccZSEhiCmrpWLp0xmdLJ+WGHIhKqKQW5XDV7GtX16hzOVEoCg3TmbDfbGo+qKUgkEC8vYefB47R3dIUdigyBksAgvXjgOB1dPeoUFgnEK2J09Tjb9kd+jqiMlFISMLNGM9tlZjvMrDYoeyR4vSNYvyMozzWzB4LtXzGzLyQdZ62ZvWpmdWb2+dGp0uiqqW8hy2DlguKwQxFJC8vnFpGXk6URSTPUYMYOWu3u586yu3+sd9nMvgH0Pjb4W0C+u19lZhOBPWb2EHAQuA/4AIlJ57eZ2SZ33zPcSoyl6roWlsyexpSC3LBDEUkLBbnZVM4tUr9Ahhp2c5AlBs65FXgoKHKgMJh7eALQCZwEVgJ17t7g7p3Aw8Atw33/sdR25iw7m04QV3+AyHmqykt45Y2TtLZ3hB2KDFKqScCBp8xsu5mtv2Dd9cBb7r4veP1j4G3gDeAA8PVgQvlZJK4GejUFZRlj6/6jdPe4+gNELlAVDJ+ypUFXA5km1SQQd/dlwM3A3Wb2vqR1t/HOVQAkvvF3A5cA84H/aWYLgL6G2uxz+EEzW29mtWZWe+TIkRRDHH019a3k52SxbG5R2KGIpJUls6YyOT+HGjUJZZyUkoC7Hw5+NwMbSXzQEzT5fAR4JGnzjwNPuPvZYPtqoJLEN/85SdvNBg73834b3L3S3StLS0sHV6NRVF3XQuW8Igpys8MORSSt5GRnce2CYj00loEG7Bg2s0Igy93bguU1wL3B6puAve7elLTLAeBGM/sRMBG4Dvi/wB5goZnNBw4B60gkjFCd7uxmX3PbgNu93dHN3jfb+F8fvHQMohLJPFXlMf7jlWaeee0IRRMHvnGivHTSkGfkk5GTyhmYAWwMJk7JAR509yeCdes4vykIEncAfR94mUQT0Pfd/SUAM7sHeBLIBu53993DrsEwffFfdvHYC4dS3v76heoPEOlL7/+NO+7fmtL2a694D3/3ieWjGZKkYMAk4O4NwNX9rLuzj7J2EreJ9rX9ZmDz4EIcPe7OM68e4X2LSrlj1dwBt58yIZcls6eNQWQimWfhjMn85K4qjp/qHHDbh7cd5Nm6Frq6e8jJ1jOrYYr0tdirb7XR+nYnv7ZkJr9y+YywwxHJeMtTvGni9Nlunt7zFi8dOsGyMt1oEaZIp+DeGZGqNDuYyJhatSDxrI06ksMX6SRQU9fCvJKJzJo2IexQRCKlZFI+l8+coqkp00Bkk0BXdw/P7z+qqwCRkMTLS9h+4BhnzmpWsjBFNgnsbDpBe0eXnv4VCUm8IkZnVw/bXz8WdiiRFtkk0NsWuapc4wCJhGHF/GJyskyjj4YsukmgvpXFM6dQXJgXdigikTQpP4er52hWsrBFMgmcOdvN9gPHqNJVgEio4uUl7Go6zonTZ8MOJbIimQRqG4/R2dVDXJ3CIqGqqojR4/C8Rh8NTSSTQHV9CzlZxsr5mh1MJEzXlE2jIDdLo4+GKJJJoKauhaVzpmnwKpGQ5edks2JeMTX16hwOS+SSwInTZ9l16ISeDxBJE1XlMV57q53mtjNhhxJJkUsCzze00uOJDikRCV/vdK1b1CQUisglgZr6Vgpys7hGg1aJpIUrLpnKlIIcPS8Qksglgeq6FlbMKyYvJ3JVF0lL2VnGqvISdQ6HJFKfhM0nz7CvuV23hoqkmaryGE3HTnOg9VTYoUROSknAzBrNbJeZ7TCz2qDskeD1jmD9jqTtl5jZFjPbHexXEJQvD17Xmdl3LJiubKz0ftPQeEEi6aW3X6BadwmNucHcI7na3c+dIXf/WO+ymX0DOBEs5wA/Aj7h7jvNrATofRzwu8B64DkSM4ytBX46rBr044GaRrp6/LyyJ3e/ydQJuSy+ZMpovKWIDFF56SSmT87n0dqDnOoceFTRqRNy+c1lsxjj75Hj0rBvlA++zd8K3BgUrQFecvedAO7eGmw3E5ji7luC1z8APswoJYGv/nQvp/sYova3ls8mO0v/cETSiZlx85Xv4YEtr/PigeMp7TO3ZCIr5umBz+FKNQk48JSZOfD37r4had31wFvuvi94vQhwM3sSKAUedvevAbOApqT9moKyUfHcn/xKn+VTCvSAmEg6+tKvX8H/WHPpgNu1d3Tx3r/62bmbPGR4Uv1EjLv7YTObDjxtZnvd/RfButuAhy445nuBFcAp4D/NbDtwso/jeh9lmNl6Es1GlJWVpRji+aZOyB3SfiISDjNL6f/t1Am5XHnJVGrqWvnMTWMQ2DiXUsewux8OfjcDG4GVcK79/yPAI0mbNwHPuHuLu58i0fa/LCifnbTdbOBwP++3wd0r3b2ytLR0cDUSkXGvqqKEFw8e41RnV9ihZLwBk4CZFZrZ5N5lEm3+LwerbwL2untyM8+TwBIzmxgkiRuAPe7+BtBmZtcF/Qi3A/86gnURkYiIl8c42+1sa9SsZMOVypXADOBZM9sJbAUed/cngnXrOL8pCHc/BnwT2AbsAF5w98eD1XcB3wPqgHpGqVNYRMa3FfOKyc22czMEytAN2Cfg7g3A1f2su7Of8h+RuE30wvJa4MrBhSgicr4JedlcU1ak5wpGQKSeGBaR8SNeHmP34ZMcP9UZdigZTUlARDJSvKIEd3hOs5INi5KAiGSkq+dMozAvm+o6JYHhUBIQkYyUm53FyvnF6hcYJiUBEclYVeUxGo68zZsnNCvZUCkJiEjGquodfVS3ig6ZkoCIZKzL3zOF4sI8NQkNg5KAiGSsrCxj1YISttS34t7nUGQyACUBEcloVRUlvHHiDPtb3g47lIykJCAiGa0qmCmwWnMUD4mSgIhktHklE7lkaoHGERoiJQERyWhmRlVFjC0NrfT0qF9gsDTNlohkvHhFCT/e3sQnH9hGQW72ufKsLOMP31/OFZdMDTG69KYkICIZ7/2LprN8bhGHjp8+r7yx5RSFedl87aN9DoQsKAmIyDhQVJjHT+6quqj8939YS406jN+V+gREZNyKV8RoOnaaA62nwg4lbaWUBMys0cx2mdkOM6sNyh4JXu8I1u+4YJ8yM2s3sz9KKltrZq+aWZ2ZfX5kqyIicr53bh/VnUP9GUxz0Gp3P/eXdPeP9S6b2TeAExds/y2Spo80s2zgPuADJCad32Zmm9x9z1ACFxEZSHlpIdMn51Nd18JtK8vCDictDbtPIJg0/lbgxqSyDwMNQPIjfCuBumC6SszsYeAWQElAREaFmRGviPGL147Q0+NkZVnYIaWdVPsEHHjKzLab2foL1l0PvOXu+wDMrBD4HPDlC7abBRxMet0UlImIjJqq8hJa3+7ktea2sENJS6kmgbi7LwNuBu42s/clrbsNeCjp9ZeBb7l7+wXH6CsF9/lkh5mtN7NaM6s9cuRIiiGKiFwsXhH0C2gGsj6llATc/XDwuxnYSKJpBzPLAT4CPJK0+bXA18ysEfgM8Cdmdg+Jb/5zkrabDRzu5/02uHulu1eWlpYOqkIiIskumTaB+bFCDSvRjwH7BILmnSx3bwuW1wD3BqtvAva6e1Pv9u5+fdK+XwLa3f1vgoSx0MzmA4eAdcDHR6wmIiL9WFVewqYdh+nq7iEnW3fGJ0vlrzEDeNbMdgJbgcfd/Ylg3TrObwrql7t3AfcATwKvAI+6++7BhywiMjjx8hjtHV3sbLrwJkYZ8EoguJunz2eu3f3OAfb90gWvNwObUw9PRGT4VpUnpqHcUt/C8rlFIUeTXnRdJCLjXnFhHotnTlHncB+UBEQkEuIVJWw/cIwzZ7vDDiWtKAmISCRUlcfo7OqhtvFY2KGkFSUBEYmElfOLyckyjSN0ASUBEYmEwvwcls6ZpqGlL6AkICKRUVURY1fTcU6cPht2KGlDSUBEIiNeXkKPw/MNuhropZnFRCQylpZNoyA3i88+soOJ+ed//K294j185cNXhhRZeJQERCQy8nOy+fMPX8X218+/Q2j34RM8WnuQL/7q5edNVB8FSgIiEikfXT6bjy6ffV7Z03ve4vd+UMsLB46dm40sKtQnICKRd+2CYrIMaiL4RLGSgIhE3pSCXJbMnhbJZwiUBERESAwr8VLTCdrOROv2USUBERESw0139zhb9x8NO5QxpSQgIgIsm1tEfk5W5EYaVRIQEQEKcrOpnFdETcT6BZQEREQCVeUx9r7ZRkt7R9ihjJmUkoCZNZrZLjPbYWa1QdkjwesdwfodQfkHzGx7sP12M7sx6TjLg/I6M/uOmdnoVEtEZPCqzs1AFp0mocFcCax296XuXgng7h8LXi8FfgI8FmzXAvyau18F3AH8MOkY3wXWAwuDn7XDrYCIyEi5atZUJufnRKpJaNjNQcG3+VsJJpx39xfd/XCwejdQYGb5ZjYTmOLuW9zdgR8AHx7u+4uIjJSc7CyuXVASqc7hVJOAA08FzTvrL1h3PfCWu+/rY7/fBF509w5gFtCUtK4pKLuIma03s1ozqz1y5EiKIYqIDF+8ooQDR09x8OipsEMZE6kmgbi7LwNuBu42s/clrbuN4CogmZldAfwV8Pu9RX0c1/t6M3ff4O6V7l5ZWlqaYogiIsMXr0iMHRSVJqGUkkBv8467NwMbgZUAZpYDfAR4JHl7M5sdbHe7u9cHxU1A8qhNs4HDiIikkYXTJxGblB+ZGcgGTAJmVmhmk3uXgTXAy8Hqm4C97t6UtP004HHgC+5e3Vvu7m8AbWZ2XdCPcDvwryNWExGREWBmVJWXUFPfSqL7cnxL5UpgBvCsme0EtgKPu/sTwbp1XNwUdA9QAfzvpFtIpwfr7gK+B9QB9cBPh1sBEZGRFq8o4UhbB/ua28MOZdQNOJ+AuzcAV/ez7s4+yv4c+PN+tq8Fojd1j4hklN45BarrWlg0Y3LI0YwuTSojInKBOcUTKSueyH+9eoQbL5t+3rrJBbkUF+aFFNnIUxIQEenDexfGePD5A9zw1z8/rzwvO4tffm41M6YUhBPYCFMSEBHpwx+tuZQV84pI7htube/kLza/wi/3tVw0RWWmUhIQEelDcWEev3HN+R/0PT3Od5+pp6Z+/CQBjSIqIpKirCxjVXkJNXXj5/ZRJQERkUGIl8d48+QZGlreDjuUEaEkICIyCL3DTdfUjY9hJZQEREQGYW7JRGZNmzBuRhpVEhARGYTeYSW2NLTS05P5/QJKAiIigxSviHHi9Fn2vHEy7FCGTUlARGSQevsFqsdBv4CSgIjIIE2fUkDF9ElUj4PhppUERESGIF5ewrb9R+ns6gk7lGFREhARGYKqihinz3az4+DxsEMZFiUBEZEhuG5BCVmW+f0CSgIiIkMwdUIuV82amvFzEaeUBMys0cx2BbOE1QZljyTNHNZoZjuStv+CmdWZ2atm9sGk8rVBWZ2ZfX7kqyMiMnZWlcd48cBx3u7oCjuUIRvMlcBqd1/q7pUA7v6x4PVS4CfAYwBmtpjEtJNXAGuBvzWzbDPLBu4DbgYWA7cF24qIZKR4RQldPc7WxqNhhzJkw24OCiaNv5V35hq+BXjY3TvcfT+J+YRXBj917t7g7p3Aw8G2IiIZqXJuMXnZWWzJ4FtFU00CDjxlZtvNbP0F664H3nL3fcHrWcDBpPVNQVl/5SIiGWlCXjbL5k7L6M7hVCeVibv7YTObDjxtZnvd/RfButt45yoAwPrY3+k74fQ58EaQaNYDlJWVpRiiiMjYi5fH+OZ/vMYj2w6Qk/XOx1xOtvErl89gUn56z92VUnTufjj43WxmG0k07fzCzHKAjwDLkzZvAuYkvZ4NHA6W+yu/8P02ABsAKisrM3+EJhEZt1ZfNp1vPP0an/vJrovW/dGaRdxz48IQokrdgEnAzAqBLHdvC5bXAPcGq28C9rp7U9Ium4AHzeybwCXAQmAriSuEhWY2HzhEovP44yNWExGREFw5ayrbvngTpzu7zytf/8Nanq1ryfwkAMwANib6f8kBHnT3J4J16zi/KQh3321mjwJ7gC7gbnfvBjCze4AngWzgfnffPSK1EBEJUenk/IvKrl8Y44Ga1znd2c2EvOwQokrNgEnA3RuAq/tZd2c/5X8B/EUf5ZuBzYMLUUQk81RVxPiHX+6n9vWjXL+wNOxw+qUnhkVERsHKecXkZFnaz0CmJCAiMgoK83O4pmxa2g8roSQgIjJKqspjvHzoBCdOnQ07lH4pCYiIjJKq8hJ6HJ7bn75NQkoCIiKj5JqyIibkZlOTxk8UKwmIiIySvJwsVswvTutpKJUERERGUby8hLrmdppPngk7lD4pCYiIjKJ4RQyAmjS9GlASEBEZRZfPnMLUCblpO9KokoCIyCjKzjJWLSihpr4V9/QbD1NJQERklMUrSjh0/DSvt54KO5SLKAmIiIyyqjTuF1ASEBEZZQtihbxnSgHVaTiEhJKAiMgoMzOqykvYUt9KT0969QsoCYiIjIGqihhH3+5k75ttYYdynvSe/FJEZJyIV5QA8OV/282C0sLz1t142Qw+sHhGGGGllgTMrBFoA7qBLnevDMo/DdxDYgaxx939j80sF/gesCw4/g/c/S+D7dcC3yYxs9j33P2rI1sdEZH0NHPqBNYsnsGOg8fZ3/L2ufK2M11U17WmdxIIrHb3c70aZrYauAVY4u4dZjY9WPVbQL67X2VmE4E9ZvYQcBC4D/gAicnot5nZJnffMyI1ERFJcxtur7yo7PvV+/nyv+3h4NFTzCmeOOYxDadP4C7gq+7eAeDuzUG5A4VmlgNMADqBk8BKoM7dG9y9E3iYRBIREYmsd4aVCOfOoVSTgANPmdl2M1sflC0Crjez583sGTNbEZT/GHgbeAM4AHzd3Y8Cs0hcDfRqCspERCJr4fRJlE7OD+0ZglSbg+Lufjho8nnazPYG+xYB1wErgEfNbAGJb/zdwCXB+l+a2X8A1sdx+7xXKkg06wHKysoGUR0RkczSe/to77ASZn19VI6elK4E3P1w8LsZ2Ejig74JeMwTtgI9QAz4OPCEu58Ntq8GKoPt5yQddjZwuJ/32+Dule5eWVpaOrSaiYhkiHh5jCNtHexrbh/z9x4wCZhZoZlN7l0G1gAvA/8C3BiULwLygBYSTUA3WkIhiSuFvcA2YKGZzTezPGAdsGnkqyQikllWlSduHw1jpNFUrgRmAM+a2U5gK4lbQZ8A7gcWmNnLJDp57/DEEHn3AZNIJIptwPfd/SV37yJxO+mTwCvAo+6+e8RrJCKSYeYUT6SseCLVdWPfLzBgn4C7NwBX91HeCfx2H+XtJG4T7etYm4HNgw9TRGR8i1eU8O8736Cru4ec7LEbzEHDRoiIpIGq8hhtHV28fPjkmL6vkoCISBqoCqlfQElARCQNlEzK57L3TB7zh8aUBERE0kRVeYzaxmOcOds9Zu+pJCAikibiFSV0dPXwwoFjY/aeSgIiImli5fxisrOMmjG8VVRJQEQkTUwuyOXq2VPHdBp8m/mNAAAGAklEQVRKJQERkTQSr4jxUtMJ2s6cHZP3UxIQEUkjq8pL6O5xnm84OibvpyQgIpJGlpUVkZ+TNWZDSysJiIikkYLcbFbMKx6z5wWUBERE0sx/WzKTa8qK6O7pc8qVETWYOYZFRGQMrFtZxrqVY/NeuhIQEYkwJQERkQhTEhARibCUkoCZNZrZLjPbYWa1SeWfNrNXzWy3mX0tqXyJmW0JyneZWUFQvjx4XWdm37GxnlFZRETOM5iO4dXufu6eJTNbDdwCLHH3DjObHpTnAD8CPuHuO82sBOh99O27wHrgORIzjK0Ffjr8aoiIyFAMpznoLuCr7t4B4O7NQfka4CV33xmUt7p7t5nNBKa4+5ZgLuIfAB8exvuLiMgwpZoEHHjKzLab2fqgbBFwvZk9b2bPmNmKpHI3syfN7AUz++OgfBbQlHTMpqBMRERCkmpzUNzdDwdNPk+b2d5g3yLgOmAF8KiZLQjK3xuUnQL+08y2A31NnNnnkxBBolkPUFZWNojqiIjIYKSUBNz9cPC72cw2AitJfJN/LGja2WpmPUAsKH+mt//AzDYDy0j0E8xOOuxs4HA/77cB2BDsf8TMXh9C3QjiGdu52sKnOo9/UasvqM6DNTfVDQdMAmZWCGS5e1uwvAa4F2gHbgR+bmaLgDwSAT8J/LGZTQQ6gRuAb7n7G2bWZmbXAc8DtwP/b6D3d/fSVCvTR+y17l451P0zkeo8/kWtvqA6j6ZUrgRmABuDuzlzgAfd/QkzywPuN7OXSXzY3xFcFRwzs28C20g092x298eDY90F/BMwgcRdQbozSEQkRAMmAXdvAK7uo7wT+O1+9vkRieafC8trgSsHH6aIiIyG8f7E8IawAwiB6jz+Ra2+oDqPGku04IiISBSN9ysBERF5FxmZBMxsbTBmUZ2Zfb6P9XcGt5buCH4+lbSuO6l809hGPnQD1TnY5lYz2xOM2fRgUvkdZrYv+Llj7KIenmHWeVyeZzP7VlK9XjOz40nrxuV5HqDO4/U8l5nZf5nZi2b2kpl9KGndF4L9XjWzDw47GHfPqB8gG6gHFpC4LXUnsPiCbe4E/qaf/dvDrsMo1Xkh8CJQFLyeHvwuBhqC30XBclHYdRrNOo/n83zB9p8G7h/v57m/Oo/n80yiP+CuYHkx0Ji0vBPIB+YHx8keTjyZeCWwEqhz9wZP3KH0MImB7MazVOr8e8B97n4MzhvL6YPA0+5+NFj3NImB+9LdcOqcqQb7b/s24KFgeTyf52TJdc5UqdTZgSnB8lTeebD2FuBhd+9w9/1AXXC8IcvEJDALOJj0ur8xiH4zuIz6sZnNSSovMLNaM3vOzDJlALtU6rwIWGRm1UHd1g5i33Q0nDrD+D3PAJjZXBLfBH822H3TzHDqDOP3PH8J+G0zayIx4vKnB7HvoGTiHMN9zUFw4S1O/wY85Ikhrv8AeIDE080AZZ4YB2kB8DMz2+Xu9aMY70hIpc45JJpH3k9iSI5fmtmVKe6bjoZcZ3c/zvg9z73WAT929+4h7JtOhlNnGL/n+Tbgn9z9G2a2CvjhaP1/zsQrgSYg+Zv9RWMQeWL46o7g5T8Ay5PW9Y6D1AD8HLhmNIMdIQPWOdjmX939bHCZ+CqJD8hU9k1Hw6nzeD7PvdZxfrPIeD7PvS6s83g+z58EHgVw9y1AAe+MzTay5znsTpIhdKrkkOj0ms87nSpXXLDNzKTl3wCeC5aLgPxgOQbs4106odLlJ8U6rwUeSKrbQaCEREfh/qDuRcFycdh1GuU6j9vzHGx3KdBI8JxPUDZuz/O71HncnmcSQ+rcGSxfTuKD3oArOL9juIFhdgyH/gcZ4h/xQ8BrJHrGvxiU3Qv8erD8l8Du4I/1X8BlQXkVsCso3wV8Muy6jGCdDfgmsCeo27qkfX+XRAdSHfA7YddltOs8ns9z8PpLJCZ0unDfcXme+6vzeD7PJO4Cqg7qtgNYk7TvF4P9XgVuHm4semJYRCTCMrFPQERERoiSgIhIhCkJiIhEmJKAiEiEKQmIiESYkoCISIQpCYiIRJiSgIhIhP1/Rx2NkxAGSjQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0.0001,8000,512)\n",
    "phi = np.linspace(0.55,0.8,60)\n",
    "dist = []\n",
    "for theta in phi:\n",
    "    u0 = np.array([          \n",
    "        [v0*np.cos(theta)],    #u\n",
    "        [np.tan(theta)],       #gamma\n",
    "        [0.0],                 #y\n",
    "        [0.0]                  #t\n",
    "        ])\n",
    "    ans0 = universal_runge_kutta(f, u0, x, butcher, a, b)\n",
    "    arr = [n[2].item(0)>=0 for n in ans0[0]].index(False)\n",
    "    dist.append(x[arr-1])\n",
    "plt.figure()\n",
    "plt.plot(phi,dist)\n",
    "max(dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}