ms_hm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate as integ\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "mpl.rc('xtick', direction='in', top=True)\n",
    "mpl.rc('ytick', direction='in', right=True)\n",
    "mpl.rc('xtick.minor', visible=True)\n",
    "mpl.rc('ytick.minor', visible=True)\n",
    "import cython\n",
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "import numpy as np\n",
    "import cython\n",
    "cimport numpy as np\n",
    "\n",
    "cimport cython\n",
    "ctypedef np.double_t DTYPE_t\n",
    "\n",
    "from libc.math cimport exp\n",
    "from libc.math cimport sqrt\n",
    "\n",
    "from cython.parallel import prange\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef dfdA(np.ndarray arr, double deltaL, double [::1] L, double c = 1, double bd = 0):\n",
    "    cdef int size = arr.shape[0]\n",
    "    cdef double [:] f = arr\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    res[0] = bd\n",
    "    #res[0] = (f[1] - f[0]) /( deltaL)  * c \n",
    "    cdef int i\n",
    "    for i in prange(1, size - 1,nogil=True):\n",
    "        res[i] = (f[i+1] - f[i-1]) /( 2*deltaL) * c \n",
    "    #res[size - 1] = (25/12*f[-1] - 4*f[-2] + 3*f[-3] - 4/3*f[-4] + 1/4*f[-5]) / deltaL * exp(-L[size-1]) * c\n",
    "    res[size - 1] = (f[-1] - f[-2]) / deltaL * c\n",
    "    return np.asarray(res)\n",
    "\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef ddfddA(np.ndarray arr, double deltaL, double [::1] L, double c = 1, double bd = 0):\n",
    "    cdef int size = arr.shape[0]\n",
    "    cdef double [:] f = arr\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    #res[0] = bd\n",
    "    res[0] = (f[1] + f[1] - 2 * f[0]) /( deltaL**2)  * c**2 \n",
    "    cdef int i\n",
    "    for i in prange(1, size - 1,nogil=True):\n",
    "        res[i] = (f[i+1] + f[i-1] - 2* f[i]) /( deltaL**2) * c**2 \n",
    "    #res[size - 1] = (25/12*f[-1] - 4*f[-2] + 3*f[-3] - 4/3*f[-4] + 1/4*f[-5]) / deltaL * exp(-L[size-1]) * c\n",
    "    res[size - 1] = (f[-3] + f[-1] - 2*f[-2]) / deltaL**2 * c**2\n",
    "    return np.asarray(res)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MS:\n",
    "    \n",
    "    def __init__(self, R, m, U, w, alpha, A, rho0, deltaL, L):\n",
    "        self.R = R\n",
    "        self.m = m\n",
    "        self.U = U\n",
    "        self.w = w\n",
    "        self.deltaL = deltaL\n",
    "        self.L = L\n",
    "        self.A = A\n",
    "        self.alpha = alpha\n",
    "        \n",
    "        self.t0 = self.alpha * np.sqrt(3 / (8*np.pi*rho0))\n",
    "        self.t = t0\n",
    "        self.xi = 0\n",
    "\n",
    "        self.RH = self.t0 / self.alpha\n",
    "        self.Abar = self.A / self.RH\n",
    "        \n",
    "        # initialize null geodesic and MS variables\n",
    "        \n",
    "        return\n",
    "    def gamma(self, R, m, U): \n",
    "        return np.sqrt(np.exp(2 * (1 - self.alpha) * self.xi) \n",
    "                       + (self.Abar * R)**2 * (U**2 - m))\n",
    "    def P(self, rho) : \n",
    "        return self.w * rho\n",
    "    def rho(self, R, m):\n",
    "        return m + (dfdA(m, self.deltaL, self.L,1,0) * self.Abar * R)  \\\n",
    "            / (3* (self.Abar * dfdA(R, self.deltaL, self.L, 1, 0) \n",
    "                  + R ))\n",
    "    def psi(self, rho, p) : \n",
    "        return 0 + \\\n",
    "        integ.cumtrapz(- dfdA(p, self.deltaL, self.L, 1, 0) / (rho + p), self.A, initial=0)\n",
    "    \n",
    "    def Pprime(self, R, m):\n",
    "        mprime = dfdA(m, self.deltaL, self.L, 1, 0)\n",
    "        Rprime = dfdA(R, self.deltaL, self.L, 1, 0)\n",
    "        mpp = ddfddA(m, self.deltaL, self.L, 1, 0)\n",
    "        Rpp = ddfddA(R, self.deltaL, self.L, 1, 0)\n",
    "        return self.w * (4 * self.Abar**2 * mprime * Rprime**2 \n",
    "                + R**2 * (4 * mprime + self.Abar * mpp)\n",
    "               + self.Abar * R * (self.Abar * Rprime * mpp + mprime * (6 * Rprime - A * Rpp))) \\\n",
    "                / (3 * (R + self.Abar * Rprime)**2)\n",
    "    def k_coeffs(self, R, m, U) :\n",
    "        g = self.gamma(R, m, U)\n",
    "        r = self.rho(R, m)\n",
    "        p = self.P(r)\n",
    "  \n",
    "        ep = np.exp(self.psi(r, p))\n",
    "        #print(g,r,p,ep)\n",
    "        kR = self.alpha * R * (U * ep - 1)\n",
    "        km = 2 * m - 3 * self.alpha * U * ep * (p + m)\n",
    "        #kU = U - self.alpha * ep * \\\n",
    "        #    ( g**2 * dfdA(p, self.deltaL, self.L, 1, 0) / \n",
    "        #     (self.Abar * R * (R + self.Abar * dfdA(R, self.deltaL, self.L, 1, 0)) * (r + p)) \n",
    "        #    + (2 * U**2 + m + 3 * p) / 2) \n",
    "        kU = U - self.alpha * ep * \\\n",
    "            (   g**2 * np.concatenate(([self.w\n",
    "                            * ((m[1] + m[1] - 2 * m[0]) /( self.deltaL**2)) *  (1**2 * 5/3)], self.Pprime(R, m)[1:] / \n",
    "             (self.Abar[1:])) ) / (R * (R + self.Abar * dfdA(R, self.deltaL, self.L, 1, 0)) * (r + p))  \n",
    "            + (2 * U**2 + m + 3 * p) / 2) \n",
    "        #kU = U - self.alpha * ep * \\\n",
    "        #    (   g**2 * np.concatenate(([self.w\n",
    "        #                    * ((m[1] + m[1] - 2 * m[0]) /( np.exp(2*L[0]))) *  (5/3)], dfdA(p, self.deltaL, self.L, 1, 0 )[1:] / \n",
    "        #     (self.Abar[1:])) ) / (R * (R + self.Abar * dfdA(R, self.deltaL, self.L, 1, 0)) * (r + p))  \n",
    "        #    + (2 * U**2 + m + 3 * p) / 2) \n",
    "        return kR, km, kU\n",
    "    \n",
    "    def run_steps(self,n_steps) :\n",
    "        step = 0\n",
    "        \n",
    "        deltau = (self.A[1] - self.A[0]) / 1 * 0.1\n",
    "        while(step < n_steps) :\n",
    "            deltau = self.cfl_deltau(R, m, U) * 0.8\n",
    "            kR1, km1, kU1 = self.k_coeffs(self.R, self.m, self.U)\n",
    "            kR2, km2, kU2 = self.k_coeffs(self.R + deltau/2*kR1, self.m + deltau/2*km1, self.U + deltau/2*kU1)\n",
    "            kR3, km3, kU3 = self.k_coeffs(self.R + deltau/2*kR2, self.m + deltau/2*km2, self.U + deltau/2*kU2)\n",
    "            kR4, km4, kU4 = self.k_coeffs(self.R + deltau*kR3, self.m + deltau*km3, self.U + deltau*kU3)\n",
    "\n",
    "    \n",
    "            self.R = self.R + deltau/6*(kR1 + 2*kR2 + 2*kR3 + kR4)\n",
    "            self.m = self.m + deltau/6*(km1 + 2*km2 + 2*km3 + km4)\n",
    "            self.U = self.U + deltau/6*(kU1 + 2*kU2 + 2*kU3 + kU4)\n",
    "\n",
    "            step+=1\n",
    "            self.xi += deltau\n",
    "  \n",
    "    def cfl_deltau(self,R, m, U):\n",
    "        a = np.exp(self.alpha * self.xi)\n",
    "        H = np.exp(-self.xi) / self.RH\n",
    "        \n",
    "        g = self.gamma(R, m, U)\n",
    "        r = self.rho(R, m)\n",
    "        p = self.P(r)\n",
    "  \n",
    "        ep = np.exp(self.psi(r, p))\n",
    "        el = (a * H * self.RH * g) / (dfdA(a * self.A * self.R, self.deltaL, self.L, 1, 1) / 1) \n",
    "        #print(el, ep.min())\n",
    "        return (np.log(1 + el / ep / np.exp(self.xi)  \n",
    "                       * self.alpha * np.concatenate( ([1e10],(Abar[1:] - Abar[0:-1])) ) / np.sqrt(self.w))).min()\n",
    "        #return (np.log(1 + el / ep / np.exp(self.xi)  \n",
    "        #               * self.alpha * (Abar[2] - Abar[1]) / np.sqrt(self.w))).min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=500\n",
    "\n",
    "L = np.linspace(-6, 1.5, n)\n",
    "A = np.exp(L)\n",
    "A = np.linspace(0, np.exp(1.5), n+1)[1:]\n",
    "A = np.linspace(0, np.exp(1.5), n)\n",
    "#A = np.concatenate( ([0], A))\n",
    "\n",
    "w = 1 / 3\n",
    "alpha = 1/2 # alpha = 1/2 for w = 1/3\n",
    "\n",
    "\n",
    "rho0 = 1\n",
    "t0 = alpha * np.sqrt(3 / (8*np.pi*rho0))\n",
    "RH = t0 / alpha\n",
    "Abar = A / RH\n",
    "\n",
    "deltaL = Abar[2] - Abar[1]\n",
    "\n",
    "deltaA = Abar[1] - Abar[0]\n",
    "\n",
    "amp = 0.3\n",
    "delta0 = amp * np.exp(-2 * Abar**2 / RH / 30)\n",
    "\n",
    "m = 1 + delta0 \n",
    "U = 1 - alpha * delta0 / 2\n",
    "R = 1 - alpha / 2 * (delta0 + w * Abar * dfdA(m, deltaL, L, RH) / (1 + 3 * w) )\n",
    "\n",
    "\n",
    "ms = MS(R,m,U,w,alpha,A,rho0, deltaL,L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.026330496766305125"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ms.cfl_deltau(R, m, U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ms.run_steps(200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 30.0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD1CAYAAACrz7WZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARpUlEQVR4nO3de4xcZ3nH8e/DOm4cJ7VDvRCwU3Bb18HlkjQrV4heQijYoVgOEa0c/uAimm0QkcI/KQ5/tLQVImpKVdSkdSetVZBaIiRC8KopBpVyKaLFMXFwnGC6daJkbUrWRDHEdWWcPP1jZ8N02cuc9ZzZM6+/H2mVmTPP2XnfOdmfz77zzNnITCRJZXnBUg9AktR7hrskFchwl6QCGe6SVCDDXZIKZLhLUoGW9fPJLrzwwrzsssu6qj1x4gSrVq0amFqAyclJhoeHl3QcddXWNbeq9YM2v6a8Fk2YX52vxaDNr+prsX///mcy86KudwDIzL59XXDBBdmtG264YaBqMzOvvPLKJR9HXbV1za1q/aDNrymvRRPmV+drMWjzq/paACezYt42dllm27ZtA1VbVRPGXNf8qn7fkufXlNeiCWOo87Wo6/s2oXaxIvv4CdWVK1fmyZMn+/Z8/TYyMsL999+/1MOoRclzA+c36EqfX0T8T2aurLJPX8/c16xZ08+n67vR0dGlHkJtSp4bOL9BV/r8gMmqO/T1zH3Dhg35+te/nm3btvXl1xJJGmRjY2OMjY1x1113jWfmhir79jXcR0ZGsuRfnSSpDhGxPzNHquzT2DdUJUmLZ7hLUoEMd0kqkOEuSQUy3CWpQIa7JBWor+F+4sQJRkdHGRsb6+fTStJAGhsbm/6AVvdXGWuzz12SGs4+d0kSYLhLUpEMd0kqkOEuSQUy3CWpQIa7JBXIPndJaij73CWpYPa5S5IAw12SimS4S1KBugr3iNgaEYcjYjwids7y+MUR8ZmI+FZEfCMiXtn7oUqSurVguEfEEHAncA2wCbg+IjbNKPsgcCAzXw28A/hYrwcqSepeN2fum4HxzDySmaeBu4HtM2o2Af8CkJnfBl4eES/u6UglSV3rJtzXAk903J9ob+v0IHAdQERsBl4GrJv5jSYnJxkZGXn+q9VqLW7UklS4Vqv1fFYCa6ruv6yLmphl28zm+NuAj0XEAeAg8ABwZuZOw8PD2OcuSQsbHR2d/gATEXG86v7dhPsEcGnH/XXAsc6CzPwB8O72IAJ4tP0lSVoC3SzL7AM2RMT6iFgO7AD2dBZExOr2YwC/C3ylHfiSpCWw4Jl7Zp6JiJuAvcAQsDszD0XEje3HdwGvAD4REc8CDwPvqXHMkqQFdLMsQ2beB9w3Y9uujttfBzb0dmiSpMXyE6qSVCDDXZIK5PXcJamhvJ67JBXM67lLkgDDXZKKZLhLUoEMd0kqkOEuSQUy3CWpQPa5S1JD2ecuSQWzz12SBBjuklQkw12SCmS4S1KBDHdJKpDhLkkFss9dkhrKPndJKph97pIkwHCXpCIZ7pJUIMNdkgpkuEtSgQx3SSqQfe6S1FD2uUtSwexzlyQBhrskFamrcI+IrRFxOCLGI2LnLI+vioixiHgwIg5FxLt7P1RJUrcWDPeIGALuBK4BNgHXR8SmGWXvAx7OzNcAVwEfjYjlPR6rJKlL3Zy5bwbGM/NIZp4G7ga2z6hJ4KKICOBC4CngTE9HKknqWjfhvhZ4ouP+RHtbpzuAVwDHgIPAzZn5XE9GKEmqrJtwj1m2zeyf3AIcAF4KXA7cERE/PXOnyclJRkZGnv9qtVoVhytJ54ZWq/V8VgJrqu6/rIuaCeDSjvvrmDpD7/Ru4Lacapofj4hHgcuAb3QWDQ8PY5+7JC1sdHR0+gNMRMTxqvt3c+a+D9gQEevbb5LuAPbMqHkceEN7EC8GNgJHqg5GktQbC565Z+aZiLgJ2AsMAbsz81BE3Nh+fBfwJ8DfR8RBppZxPpCZlf+lkST1RjfLMmTmfcB9M7bt6rh9DHhTb4cmSVosP6EqSQUy3CWpQIa7JBXI67lLUkN5PXdJKpjXc5ckAYa7JBXJcJekAhnuklQgw12SCmS4S1KB7HOXpIayz12SCmafuyQJMNwlqUiGuyQVyHCXpAIZ7pJUIMNdkgpkn7skNZR97pJUMPvcJUmA4S5JRTLcJalAhrskFchwl6QCGe6SVCD73CWpoexzl6SC2ecuSQIMd0kqUlfhHhFbI+JwRIxHxM5ZHr8lIg60vx6KiGcj4oW9H64kqRsLhntEDAF3AtcAm4DrI2JTZ01m3p6Zl2fm5cCtwJcz86kaxitJ6kI3Z+6bgfHMPJKZp4G7ge3z1F8PfLIXg5MkLU434b4WeKLj/kR720+IiAuArcCnz35okqTFWtZFTcyyba7+yW3A1+ZakpmcnGRk5MfdPKOjo9M9nJKkDq1Wi1arNX13TdX9uwn3CeDSjvvrgGNz1O5gniWZ4eFh7HOXpIV1nvxGxPGq+3ezLLMP2BAR6yNiOVMBvmdmUUSsAn4D+GzVQUiSemvBM/fMPBMRNwF7gSFgd2Yeiogb24/vape+Ffh8Zp6sbbSSpK54+QFJajgvPyBJAgx3SSqS4S5JBfJ67pLUUF7PXZIK5huqkiTAcJekIhnuklQgw12SCmS4S1KBDHdJKpB97pLUUPa5S1LB7HOXJAGGuyQVyXCXpAIZ7pJUIMNdkgpkK6QkNZStkJJUMFshJUmA4S5JRTLcJalAhrskFchwl6QCGe6SVCD73CWpoexzl6SC2ecuSQIMd0kqkuEuSQXqKtwjYmtEHI6I8YjYOUfNVRFxICIORcSXeztMSVIVyxYqiIgh4E7gjcAEsC8i9mTmwx01q4G/ArZm5uMR8aKaxitJ6kI3Z+6bgfHMPJKZp4G7ge0zat4O3JOZjwNk5pO9HaYkqYpuwn0t8ETH/Yn2tk6/CFwcEV+KiP0R8Y7ZvtHk5CQjIyPPf7VarcWNWpIK12q1ns9KYE3V/RdclgFilm0zm+OXAVcCbwBWAF+PiH/PzO90Fg0PD2OfuyQtbHR0dPoDTETE8ar7dxPuE8ClHffXAcdmqTmemSeBkxHxFeA1wHeQJPVdN8sy+4ANEbE+IpYDO4A9M2o+C/xaRCyLiAuAXwEe6e1QJUndWvDMPTPPRMRNwF5gCNidmYci4sb247sy85GI+BzwLeA54G8z86E6By5JmpvXlpGkhvPaMpIkwHCXpCJ5PXdJaiiv5y5JBXPNXZIEGO6SVCTDXZIKZLhLUoEMd0kqkOEuSQWyz12SGso+d0kqmH3ukiTAcJekIhnuklQgw12SCmS4S1KBDHdJKpB97pLUUPa5S1LB7HOXJAGGuyQVyXCXpAIZ7pJUIMNdkgpkuEtSgexzl6SGss9dkgpmn7skCTDcJalIXYV7RGyNiMMRMR4RO2d5/KqIOBERB9pff9D7oUqSurVsoYKIGALuBN4ITAD7ImJPZj48o/SrmfmWGsYoSaqomzP3zcB4Zh7JzNPA3cD2eoclSTob3YT7WuCJjvsT7W0zvTYiHoyIf46IX+rJ6CRJi7LgsgwQs2yb2T/5TeBlmflMRLwZuBfYMHOnyclJRkZ+3M0zOjo63cMpSerQarVotVrTd9dU3X/BPveIeC3woczc0r5/K0BmfmSefR4DRjLzeOd2+9wlqbq6+tz3ARsiYn1ELAd2AHtmPPElERHt25vb3/f7VQYiSeqdBZdlMvNMRNwE7AWGgN2ZeSgibmw/vgt4G/DeiDgDnAJ2ZD8/+ipJBbr3gaPcvvcwyy/5hSur7tvNmjuZeR9w34xtuzpu3wHcUfXJJelcMh3Wx54+xUtXr+CWLRu59orZ+lOmam+95yCnfvTsop6rq3CXpHNJ1RDupnZmWB99+hS33nMQYNb62/ceXnSwg+Eu6RxRRwhXqZ0trE/96Flu33t41nEce/rUImc6xWvLSBpY9z5wlNfd9kXW7/wnXnfbF7n3gaNz1t16z0GOPn2K5MchPFv9fCF8NrVzhfVc21+6esWs27vl9dwlNcpSB3aVEK5SO1dYz7X9li0bOfPoPr7/ub/kuf89OWvNfPoa7qtWraLVarFt27Z+Pq2kHus2gBdTu9SBXSWEq9TesmUjK84b+n/bVpw3xC1bNs76Pa69Yi137LyBV+/4fV5w/spZa+bjsowkoJ4z5iq10IzArhLCVWqvvWItH7nuVaxdvYIA1q5ewUeue9Wcb9RO7/O1nVdz+r/H989ZNAffUJUGzFJ3clR5Y7BXbyLOFdhHuwzyW7Zs/Im2wvlCeHrsC71uVWqn6+cL814y3KUaVAngKvVN6OSoa00amhHY0/XdhnA/A7sKw13ntKU+C65aX9dZc11nzFVqoTmBXQLDXUumjmCtUtuEs+Cq9XV2ctRxxlylFgzsXjLcz0FLHarTtXUE66CdBVfdXtdZc1PWpKf3MbDPXl/D/cH/OsYlm9/Ce97+Nj78/nfNWtOU4KljvbQJtU0IVagvWAftLLhqfV1nzefimvQgGBsbm/5c0Kqq+/a1FfIF56/k/Kvfyz3fv+Ss26aaUNuUcTTh03dVaqG+YB3Efua62umqtt5Nt909ettv8bWdVxvIDbBt27bpP9hxouq+S9Ln3u8wqTOkmjCOJnz6rlcfrT7bYB3EfuY6Q9jAPnct2Zp7P8OkzpBqwjiasGZbZ1dEE5Yi6n6jz6UL9dqShXs/w6TOkGrCOJqwZltnV0RTQtgA1iBZknDvd5jUGVJNGEcTzlbr7oowhKVqFvwD2b30Uy/ZkCM3/02ju0mq1jZlHFXHLGlwLOYPZPc13EdGRvL+++/v2/NJUgkWE+5ez12SGmpsbIzR0VFYRJ+7Z+6S1HCNP3OfnJzs59P1XfvDBkUqeW7g/AZd6fMD1lTdoa/hfvz48a5rqyzdNKEWqv0P1oQxV6mta25V6wdtfk15LZowvzpfi0Gb3yKWpoer7tDYv8TUhBe1zvcGmjDmuubXlEBrwvya8lo0YQx1vhZ1fd8m1C5WX9fcI+KHwOyf4/9Jq+j+egpNqIWpX526/fWkCWOuUlvX3KrWD9r8mvJaNGF+db4Wgza/qq/Fxsy8qEJ9f8NdktQfjV2WkSQtnuEuSQXqS7hHxNaIOBwR4xGxsx/P2U8R8VhEHIyIAxEx8I38EbE7Ip6MiIc6tr0wIr4QEf/Z/u/FSznGszHH/D4UEUfbx/BARLx5Kce4WBFxaUT8a0Q8EhGHIuLm9vYijt888yvl+J0fEd+IiAfb8/uj9vbKx6/2NfeIGAK+A7wRmAD2Addn5sO1PnEfRcRjwEhmdt/r2WAR8evAM8AnMvOV7W1/CjyVmbe1/4G+ODM/sJTjXKw55vch4JnM/LOlHNvZioiXAC/JzG9GxEXAfuBa4F0UcPzmmd/vUMbxC2BlZj4TEecB/wbcDFxHxePXjzP3zcB4Zh7JzNPA3cD2PjyvFikzvwI8NWPzduDj7dsfZ+oHaiDNMb8iZOZ3M/Ob7ds/BB4B1lLI8ZtnfkXIKc+0757X/koWcfz6Ee5rgSc67k9Q0MFoS+DzEbE/IkaXejA1eXFmfhemfsCAFy3xeOpwU0R8q71sM5DLFp0i4uXAFcB/UODxmzE/KOT4RcRQRBwAngS+kJmLOn79CPeYZVtp/Zevy8xfBq4B3tf+tV+D5a+BnwcuB74LfHRJR3OWIuJC4NPA+zPzB0s9nl6bZX7FHL/MfDYzLwfWAZsj4pWL+T79CPcJ4NKO++uAY3143r7JzGPt/z4JfIappajSfK+93jm97vnkEo+npzLze+0fqueAuxjgY9heq/008A+ZeU97czHHb7b5lXT8pmXm08CXgK0s4vj1I9z3ARsiYn1ELAd2AHv68Lx9EREr22/sEBErgTcBD82/10DaA7yzffudwGeXcCw9N/2D0/ZWBvQYtt+Q+zvgkcz8846Hijh+c82voOM3HBGr27dXAL8JfJtFHL++fEK13Zb0F8AQsDszP1z7k/ZJRPwcU2frMPVnC/9x0OcXEZ8ErmLqI93fA/4QuBf4FPCzwOPAb2fmQL4pOcf8rmLqV/oEHgN+b3qNc5BExK8CXwUOAs+1N3+QqXXpgT9+88zveso4fq9m6g3TIaZOvj+VmX8cET9DxePn5QckqUB+QlWSCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUoP8Db5tt0Nt1euEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.arange(n), ms.U[:])\n",
    "plt.xlim(0,30)\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fb97f4a6940>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD1CAYAAABJE67gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa+0lEQVR4nO3df5BdZX3H8fd3kwU3wuyCWTu6EIOOjWJBVrc1NrZF7BiQXlisLSJo20F2mDodsZ3UIFFAocUyVabToc6NMtohYsYSrkSKlFFa6o/QJu5KSCOKQmMWplkKRCsZCcm3f+y94bo5v+7dc889Pz6vmTvZPfe55z7Pwn733O/zfZ5j7o6IiBTPQL87ICIi3VEAFxEpKAVwEZGCUgAXESkoBXARkYJSABcRKailWb7Z8uXLfeXKlQDs37+f4eHhxK+dm5tjdHQ0cftOz99p+yzeQ2Puf3+yeI+8jTmLn5HG3Nn5d+zY8aS7H/1id8/s8cY3vtFbLrvsMu9E+2uT6PT8nbbP4j005vTPrzGn376b12jMnZ0f2O4BMbVvKZRarZar83fTnyzeo5fn15jz8x69PH8ef0a9Pn9Vxmye4UrMiYkJ3759e7evpdvXFpXGXA0aczUsZsxmtsPdJxYez/QKfP/+/UxNTbF169aOXzs1NdWDHuWbxlwNGnM1dDPmrVu3tl4XmGwvzBW4iEhV5eIKXERE0pNpGWE3GtOzXLnlQQ4cPHzkmAEXr17BdZOn9a9jIgtsaOxk07Y9aH9PCTJg8J43pRu3ch3AG9Oz/PnmGQ4vOO7Ardv28OXte/nk75/O5PhYP7onFaCgLGk57PNxC0gtiOc6gN94z8NHBe92v3j+MFdsnmH7fz+lq3HpWtCnPJFeue2Bn1QjgD/+zIFE7dL+qyblpqtq6adDKRaO5DqAv3xkiFkFcVkkXWFLniwxS+1cua4DX7d2VUcdvHXbHhrTs911TkqlMT3Laz96NyvX38UVm2cUvCU3LnrTyYnbFr4OvNOrp6HBAXZ/4pxuuiclsKGx88inMZE8WUwVSlgdeK5TKACT42NHqkyS/HIeOHiYxvSsKlMqJG8pkhcfs4TrLzhN/w9Kz+U+gLdr/eWKC+JXbnlQvzwVkGXg7kUNr8hiFSqAw3wQn3jFiXxo80xoFcGBg4fZ0NipX7aS6lXg1gIxKZrCBXDgyNX1FZtnQtvcum0PE684UVfiJXPxxu/wrR89lcq5dFUtRVfIAA7zQfyqO3by8+cOhba55s5dCuAl0ZieZd2XZ1jMRbeusKVsChvAAa6/4LTIq/BnDhzMrjPSM4u96r5EQVtKKtd14HEmx8e4ZPWKyDYbGjtTeS/JXmN6lld/5K6ugvexSwe46cIzeOyGcxW8pbAKXweexOs+9rXQVIoBn77wDKVSCqbbem6V8EkZlXo/8OsvCL/CcuY3xZLi6CZ4v/iYJdx04Rns+vjZCt5SGYXOgbdMjo9x7dZdPP1scM476X4q0n+d5ruPXTqgLYWlskpxBQ5wde11kc9rj5T86zR4X7J6BQ9fd46Ct1RWbAA3s5PN7D4z221mu8zsgwFtLjazB5uPb5vZ63vT3XBxv8RKo+TbhsbOxMG7NUGpyUmpuiQplOeBv3D375rZ8cAOM7vX3f+rrc2jwO+4+9Nmdg5QB97Ug/5GGovYflZplPzq5MpbJYEiL4i9Anf3J9z9u82vfwbsBsYWtPm2uz/d/HYbcFLaHU1i3dpVkc8rjZI/Ct4i3esoB25mK4Fx4IGIZpcCdwc9MTc3x8TExJFHvV7v5O1jxaVRrrlzV6rvJ4uTNHgrZSJVVK/Xj8RKYHlQm8R14GZ2HPBvwPXuviWkzVuBm4G3uPv/Lny+V3Xg7dbc8I3IdMlNqgnPhaTBe3AAfvhX52bQI5H8WlQduJkNArcDmyKC9+nAZ4Hzg4J3VtatXUXUDYs0mdl/SScsBwxu/IMzet8hkYJKUoViwOeA3e7+qZA2K4AtwHvd/QfpdrEzk+NjXByxvF6Tmf2VdJHOsUsH+NQf6tOSSJQkVShrgPcCO81spnnsI8AKAHf/DPAx4CXAzfPxnueDLvezct3kaXzxgT0cDsgOpXc7UelUY3o2UfBe86oT2XTZmzPokUixxQZwd/8mMXHP3d8PvD+tTqUhKHjD/NJ63XKtP666I35jMQVvkeRKsxJzobGRodDnlAfP3oZG9N7tMD9hqeAtklxpA3hUTbjy4NlKkvfWhKVI5wq9H3iUyfExBkISP8qDZydJ3ttAE5YiASqxH3iYlevvCn1O9eDZiNqrvUX/LUSilXo/8DBReXCtyuy9xvRsbPAeGhxQ8BbpUqkDeFQe/JkDB7U3So/FVZ0MGPz1O0/PqDci5VPqAD45PsYJywZDn1c1Su/EVZ0o7y2yeKUO4BB9owdVo/RGkolL3adUZPFKH8BVjZK9uNSJ8t4i6Sh9AIf4VZmSniQLdpT3FklHaevA22lVZjaSpE4uWb1CV98iCVW6DrylMT3LFZtnAp8z4NEbtN90GuJqvocGB9j9iXMy7JFIOVSyDrxlcnyMZYPBQx0eCq9SkeSS1HwrdSKSrkoEcIBjB5cEHn/u+eigI8nETVwqdSKSvsoE8GeePRh4/NmDhzWRuUhxE5dDgwO6n6VID1QmgL9cy+p7IsnEpVInIr1RmQCuZfW9ce3W6D9+qvkW6Z3KBHAtq++Np0NSUy26+hbpnUrUgbdELat/XMvqOxb3qUUTlyKLozrwBU796N08e/DwUcdHhgaZufrtfehRcUXVfavmWyQ9la4DbxdWTmjaGKUjcZUnSp2I9F7lAnhYOWFcLlde0JieZVNM5YlSJyK9V7kAHlZOaGhjq6RuvOdhohJvI1rdKpKJ2ABuZieb2X1mttvMdpnZBwPamJn9nZk9YmYPmtkbetPdxVu3dlXgNrKOKlGSittH/ZrzwieLRSQ9Sa7Anwf+wt1fC6wGPmBmpy5ocw7w6uZjCviHVHuZosnxsdCrR93gIZ4qT0TyIzaAu/sT7v7d5tc/A3YDC39Dzwf+0edtA0bM7GWp9zYlS0JmLMOOywviFu5oybxIdjrKgZvZSmAceGDBU2PAT9q+38vRQZ65uTkmJiaOPOr1eofdTcehkNLJsOPygqjJ3qh910WkM/V6/UisBJYHtVma9GRmdhxwO3CFu/904dMBLzkqGo6OjtLvOnCYDzRB6ZLWRKZSAME2NKJ3HIzarkBEOjM1NdVaxIOZPRnUJtEVuJkNMh+8N7n7loAme4GT274/CXi8o95mSBOZnYsrHdSeJyLZS1KFYsDngN3u/qmQZncC72tWo6wG9rv7Eyn2M1WayOxcXOmgFu6IZC/JFfga4L3AWWY203y8w8wuN7PLm23+Gfgx8AiwEfjT3nQ3PWH5WtWDB4v6wzYyNKirb5E+iM2Bu/s3Cc5xt7dx4ANpdSoL69au4kObZ466qmylURSQXhD3B0113yL9UbmVmC1KoyQXNy+gP3Yi/VHZAA6qB08q6g+aSgdF+qdS+4EvpHrweHHpE5UOivSO9gOPsOaGb4TWg3/6wjOUGgDGP/4vkYt3Hrvh3Ax7I1JN2g88gOrBozWmZ7XyUiTHKh3ANZEZLWrfE0PpE5F+q3QAB9WDR4m6+r5Yuw6K9F3lA7jSKMHi/nhp10GR/qt8AI9Ko1T5TvVR6RPdcUckHyofwCE8IA1XNFDFTV5q5aVIPlS6DrwlbN1OVdfzxF19K/ctko24OvDE+4GnYXh4uG83cYiiO9X/Ml19i+RDrVajVquxcePG/UHPK4WC7lTfLm68uvoWyQ8FcFSJ0i5qvJq8FMkXBXBUidIuagGT0ici+aIA3hS2oCcsvVJWA2ETuih9IpI3CuBNb33NaEfHy6gxPcvhkI8i2p9RJH8UwJvu+/5cR8fLKKp8UBtXieSP6sCbwnLdVdnUKm7xjjauEsme9gNPqOp7g4eNH+arT2aufnvGPRKRFu0HHqPqpYSqPhEpHgXwpirvDR61eEfVJyL5pQDepqp7g0d9wlD1iUh+xQZwM7vFzPaZ2UMhzw+b2VYz+56Z7TKzP0m/m9moahpFd50XKaYkV+CfB86OeP4DwH+5++uBM4G/NbNjFt+17FV1RWbY4h1Q9YlInsUGcHe/H3gqqglwvJkZcFyz7fPpdC97VdsbPGrxDij/LZJnaWwn+/fAncDjwPHAhe5+OKjh3NwcExMvVMJMTU21ahxzo2p7g2vxjkg+1ev19u23lwe1SSOArwVmgLOAVwH3mtm/u/tPFzYcHR0lr3XgLWF7g4cdLzIt3hHJr/YLXDN7MqhNGlUofwJs8XmPAI8Cr0nhvH0RtnlVGVMocVvHKn0ikm9pBPA9wNsAzOxXgFXAj1M4b1+sW7uKwYBZvZ8/93zpSgm1eEek2JKUEd4GfAdYZWZ7zexSM7vczC5vNvkE8JtmthP4OvBhdw+83C+CyfExjnvR0Zmlg4e8dKWE2jpWpNhic+DuflHM848DpdooIyzfXaZSQm0dK1J8WokZICzfXaY8eNSnCVWfiBSDAniAKpQSRn2aUPWJSDFoP/AAYSmUqJK7ohlZFvxpYtnggPLfIjkRtx94pgF8eHiYer1OrVbL8m07FlZKWKZNrX5x8FDg8WOWLsm4JyISplartRbz7A96XimUAGXf1KoxPcuzBwMXy7L/QHk+ZYiUnQJ4gLJvahW1fD7s04eI5I8CeIiwSoyiV6Jo+bxIeSiAhyjrikwtnxcpDwXwEGVdkanl8yLloQAeoYwrMrV8XqQ8VAceoWwrMrV8XqRY4urA09gPPLFWHXhRlG1FppbPixRLrVajVquxceNG1YF3qmw3d4jKf6v6RKR4FMAjlOnmDo3p2cDFSaDqE5GiUgCPUKZSwhvveTgwz22o+kSkqBTAI5SplDAsfeKo+kSkqBTAY5SllDCsfHBJUWdkRUQBPE4ZSgmjygcPuQoIRYpKdeAxylBKqPJBkWJSHfgilaGUUHffESkm1YEvUhlKCXX3HZFyUgCPUYZSQt19R6ScFMBjFL2UUHffESmv2ABuZreY2T4zeyiizZlmNmNmu8zs39LtYv+F5bujlqbnhe6+I1JeSa7APw+cHfakmY0ANwPnufvrgD9IpWc5UtSbHOvuOyLlFhvA3f1+4KmIJu8Btrj7nmb7fSn1LTeKepNj3X1HpNzSyIH/KnCCmf2rme0ws/eFNZybm2NiYuLIoyglhUW9yXFU37T/iUi+1ev1I7ESWB7UJo068KXAG4G3AUPAd8xsm7v/YGHD0dFRtm/fnsJbZm9sZCgw553nPPLIssHAFIrKB0Xyb2pqqrWIBzN7MqhNGlfge4GvufvP3f1J4H7g9SmcN1fe+prRjo7ngcoHRcotjQD+FeC3zGypmS0D3gTsTuG8uXLf9+c6Ot5vKh8UKb/YFIqZ3QacCSw3s73A1cAggLt/xt13m9nXgAeBw8Bn3T205LCowvLJec2BR01g5jntIyLJxQZwd78oQZsbgRtT6VFOvTwkB57XJfXa/0Sk/LQSM6GiLanX/ici5acAnlDRltRrAlOk/LQfeAeKsqReE5gi5RC3H3imAby1H3itVsvybVNTlCX1msAUKYdardZa8Kj9wBerKEvqoz4RaAJTpDwUwDtQhCX1jenZwD8yoP1PRMpGAbxDYfeQzEtq4sZ7Hg78I2No/xORslEA71Del9SHpU8cdPUtUjIK4B3K+5L6JRacQAk7LiLFpQDeobBcd15KCQ95cJY+7LiIFJfqwDuU51LCqAnMsNy9iORXXB24eYZXZhMTE17U/cBbGtOzfGjzTOBE4djIEN9af1bmfWpZc8M3Aj8JGPDpC89QDlykoMxsh7tPLDyuFEqH8lxKqAlMkWpRAO9CXksJNYEpUi0K4F3IaymhJjBFqkUBvAt5LCXUBKZI9SiAdyGPpYRRKzC1/4lIOSmAdyGPpYSawBSpHtWBdyFvuxIqfSJSTqoD75GV6+8KPG7Aozecm2lfVP8tUm6qA09Z2JVtP25yrPSJSDUpgHcpTzc5Vv23SDUpgHcpTzc5Vv23SDXFBnAzu8XM9pnZQzHtft3MDpnZu9LrXr6F3eQ4yyX1msAUqa4kV+CfB86OamBmS4BPAvek0KfCCMt3Z5kHV/23SHXFBnB3vx94KqbZnwG3A/vS6FRRhKWYs0w9h13tawJTpPwWnQM3szHgAuAzcW3n5uaYmJg48qjX64t9+74KS6GEHe+FkWXBV/snhBwXkWKo1+tHYiWwPKjN0bNwnbsJ+LC7H7KYS8/R0VHKUgcO8ysyg0r4skyh/OLgocDjmr8UKbapqanWIh7M7MmgNmlUoUwAXzKzx4B3ATeb2WQK5829fpcSNqZnefbg4cDn9h/I7lOAiPTHogO4u5/i7ivdfSXwT8Cfuntjsectgn6XEka9R7/3JheR3otNoZjZbcCZwHIz2wtcDQwCuHts3rvs+llKGLX7oSpQRMovNoC7+0VJT+buf7yo3hTQ8NAgzwSkK3qdB2/VfwelukeGBlWBIlIBWom5SP0qJYyq/77mvNf19s1FJBcUwBcpLIXydI9LCbWBlYhoP/BF6tfNHbSBlUj5aT/wHmtMz/KhzTOB6YyxkSG+tf6snrxv2H7kAI9lvB+5iPSW9gPvkcnxscDgDb2rRNEGViICCuCpyPrmDtdu3aUNrEREATwNWa7IbEzPhk6QagJTpFoUwFOQ5YrMqPMpfSJSLQrgKclqRWbU+ZQ+EakWBfCUZHVzh7DtY5cNDih9IlIxqgNPSVYrMsO2jz1m6ZJ030hE+i6uDjyN/cATGx4eLvxNHMJksSJT28eKVEutVqNWq7Fx48b9Qc8rhZKSLFZkavtYEWmnAJ6SdWtXBS6ucaIDbyc0gSki7RTAU5LFikxNYIpIOwXwFI30uBJFE5gi0k4BPEW9rETRBKaILKQAnqJeVqJoAlNEFlIdeIp6WYmi+1+KVI/2A89Q1N7gI0ODzFz99lydV0SKQfuBZyCqEuWZAwe7vgrX/S9FJIgCeMqidgTsth5c978UkSAK4CmLykd3Ww8esNU4oPtfilRdbAA3s1vMbJ+ZPRTy/MVm9mDz8W0ze3363SyOyfExlg0G/1i7qQdvTM9yOCQvcyjD+QsRyZ8kV+CfB86OeP5R4Hfc/XTgE0A5d6vqwLGDwQtrurlg1g0cRCRM7G6E7n6/ma2MeP7bbd9uA05KoV+FlmY9uMoHRSRM2jnwS4G7w56cm5tjYmLiyKOsW8umWQ8elv82NIEpUmb1ev1IrASWB7VJbT9wM3sr8wH8LWFtRkdHKXMdeMu6tasC67ZbOxMmDbxR+W9lv0XKbWpqqrWIBzN7MqhNKlfgZnY68FngfHf/3zTOWWRR9eBRKZGFrt26K/Q55b9FZNEB3MxWAFuA97r7DxbfpXIIK/FLOo/ZmJ6NzJkr/y0iScoIbwO+A6wys71mdqmZXW5mlzebfAx4CXCzmc2YWflzJAmElfg5yfLgUdUnI0ODyn+LSKIqlItinn8/8P7UelQSYyNDoemSa+7cFRuAo1ItWj4vIqCVmD0TleJIsi+Kqk9EJI4CeI9Mjo9xQsgt0CA6RaLqExFJQvuB99DVtfBUR1SKRKsvRQS0H3jfvfLKuwKvpg149IZzA1+zcv1doee76cIzlEIRqRjtB94nUamQDY2dRx2Pyo0r/y0i7RTAeywq5bFp256jAnbU4h3lv0WknQJ4j0VVo7SW1rfELd5R/ltE2imA91hcNUr7ZOZVdxydUmkxtPpSRH6ZAngGoqpRYP7Ke0NjJz9/7lBom4tXr1D+W0R+iQJ4BuIC77ovz3Drtj2Rba6bPC3NLolICagOPCNR+euDh6NfO9LFrdhEpPhUB54TjelZrtg809VrVfstUm2qA++zyfExXnxM8L0yo1yi3LeIhFAAz9D1F3SWxx4aHFDuW0RCKYBnaHJ8jEtWr0jc/q/feXoPeyMiRacAnrHrJk9jzatOjG235lUnKnUiIpEUwPtg02VvDr0SH7D5vPemy96cca9EpGhUhSIiknO5qEKpch24iEinVAcuIlJwubgCX4x6vd7vLmROY64GjbkaejHmvgXwTtMonQ6+0/N3k9bp9XtozOmfX2NOv303r9GY0zl/YQJ4r8+fx1/sTmnM6bfP6j16ef48/ox6ff6qjDnTHLiZzQH/3fx2GNjfwcuXA0920L7T83faPov30Jj7358s3iNvY87iZ6Qxd3b+V7j76MIGmQZwERFJT2EmMUVE5JcpgIuIFFTuA7iZnW1mD5vZI2a2vt/9SYuZ3WJm+8zsobZjJ5rZvWb2w+a/J7Q9d2XzZ/Cwma3tT68Xx8xONrP7zGy3me0ysw82j5d23Gb2IjP7DzP7XnPM1zaPl3bMLWa2xMymzeyrze9LPWYze8zMdprZjJltbx7r7ZjdPbcPYAnwI+CVwDHA94BT+92vlMb228AbgIfajv0NsL759Xrgk82vT22O/VjglObPZEm/x9DFmF8GvKH59fHAD5pjK+24mb8f9XHNrweBB4DVZR5z29j/HPgi8NXm96UeM/AYsHzBsZ6OOe9X4L8BPOLuP3b354AvAef3uU+pcPf7gacWHD4f+ELz6y8Ak23Hv+Tuv3D3R4FHmP/ZFIq7P+Hu321+/TNgNzBGicft8/6v+e1g8+GUeMwAZnYScC7w2bbDpR5ziJ6OOe8BfAz4Sdv3e5vHyupX3P0JmA92wEubx0v3czCzlcA481ekpR53M5UwA+wD7nX30o8ZuAn4S6D9jq9lH7MD/2JmO8xsqnmsp2NeuojOZsECjlWx7rFUPwczOw64HbjC3X9qFjS8+aYBxwo3bnc/BJxhZiPAHWb2axHNCz9mM/s9YJ+77zCzM5O8JOBYocbctMbdHzezlwL3mtn3I9qmMua8X4HvBU5u+/4k4PE+9SUL/2NmLwNo/ruvebw0PwczG2Q+eG9y9y3Nw6UfN4C7PwP8K3A25R7zGuA8M3uM+bTnWWZ2K+UeM+7+ePPffcAdzKdEejrmvAfw/wRebWanmNkxwLuBO/vcp166E/ij5td/BHyl7fi7zexYMzsFeDXwH33o36LY/KX254Dd7v6ptqdKO24zG21eeWNmQ8DvAt+nxGN29yvd/SR3X8n87+w33P0SSjxmM3uxmR3f+hp4O/AQvR5zv2duE8zsvoP5aoUfAVf1uz8pjus24AngIPN/jS8FXgJ8Hfhh898T29pf1fwZPAyc0+/+dznmtzD/MfFBYKb5eEeZxw2cDkw3x/wQ8LHm8dKOecH4z+SFKpTSjpn5SrnvNR+7WrGq12PWUnoRkYLKewpFRERCKICLiBSUAriISEEpgIuIFJQCuIhIQSmAi4gUlAK4iEhBKYCLiBTU/wNpvxKtO3PGdQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.arange(n), ms.gamma(ms.R, ms.m, ms.U))\n",
    "#plt.xlim(0,30)\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fb98031caf0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD1CAYAAACrz7WZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYGElEQVR4nO3dfZBddX3H8fdnNxsNWFlt1hYSaFB3bNeHEtyhsXQ6Po0SfAh9BqG2TIc7zEAr2tKJhRa1dkqL49MMQ+aidGSgYrWy3di0wal1OnWMBdxAiDF1RQrZpGUZJbYmYwL59o97Nl5v7sO5u+fee865n9fMnew953fu/n4L+ezJ7/weFBGYmVm5jAy6AmZmlj2Hu5lZCTnczcxKyOFuZlZCDnczsxJyuJuZldCqQVcAYO3atbFhwwYADh8+zBlnnJH62sXFRSYmJlKX7/bzl3NNr8u7zdl//nKucZuzL+82d/f5Dz744FMR0fziiBj466UvfWlcddVVMTs7G1dddVV049WvfnVX5bv9/OVc0+vybnP2n7+ca9zm7Mu7zek+fykrgW9Fi1zNxZ37GWecQbVa7cv3etvb3tbza3pdvltuc2+ucZuzL9+tYW3z0uv2228/3KqcIgczVKenp+OBBx5Y7rUs99qicpuHg9s8HFbSZkkPRsR0s3OFf6BaqVQGXYW+c5uHg9s8HHrV5sLfuZuZDatS37mbmdmpcvFA1fLjxpk93L3rcVby77kRwTt+4Rw+eMkrM6uXmXXH4T7kZuYWeO/nH+bo8ROZfeaJgLt2Pc5dux7/seMCLt/k0DfrB4f7kJqZW+D6z+4mw0zvKDg19H2Xb9YbuXigOjk5Ga973etOjt203unFnXrWHPhm7W3fvp3t27dz++23z0fEZLMyuQh3j5bpj8tv/ypf+fZ3B12NZTl99Sh/8Suv5JKN6wZdFbPcaDdaxuE+BAbRBdMvvsu3YdYu3N3nXnI3zuw55cFmmbR6eAu+27fh5jv3EltJsC8nGIvQn9/IvwCsyNwtM4SWE+wXvuSF3H3VazKvSxFDPy3/crBBcrgPmZm5Ba77zO7U5QcRUGUOfLNuveC0MW5628u7/jvoPvchMjO3wLtTBvtzVo3wV7/2qoHcdV6ycd0p3zeL2bFmRfS9I8e5/nMPAWT299F37iUyM7fAe/5uNydS/Ce9oiAzRR34NkzWja/hK1tfn7p87u/cDx8+TKVS8SSmFbrh3j0dg13AR37rvML0EX/wklc2/SXk0LcyOvj00VTlliYxAS338/Ode0mkfYD60QIF+3LMzC3wvtm9PH30+KCrYta10t2528qkDfYrNp1T6mCH5n35zZR9/L8Vz9iouP7NL8vs83znXnBpR8YUpY+9CPyvA8uaR8vYKW64d0/HMg72bKX914HZIHknpgKbmVvgB8eebVtmzdiIg91sCKUKd0kXSdovaV7S1ibnf1bSVyX9UNIfdXOtLV+nu/YRwV/+6qv6VBszy5OO4S5pFLgV2AxMAZdJmmoo9l3gD4APLeNaW4YbZ/a0vWsX8OHfLPfIGDNrLc2d+wXAfEQ8GhHHgHuALfUFIuLJiLgfaHzC1PFa697M3AJ3dxjpUaSx7GaWvTThvg54ou79geRYGqmuXVxcZHp6+uSrWq2m/PjhdMvO/W0n76wZG3Gwm5VUtVo9mZXA2lbl0oyWUZNjacdPprp2YmICD4VMb6HDLDb3s5uVV6VSoVKpACDpqVbl0ty5HwDOrnu/HjiYsh4rudaamJlbaHt+GCYqmVlnacL9fmBS0rmSVgOXArMpP38l11oTnUbIeNijmUGKbpmIeEbStcBOYBS4IyL2Sro6Ob9N0k8DDwDPB05Iug6YiojvN7u2R20pvU7j2teNr+ljbcwsz1LNUI2IHcCOhmPb6r7+b2pdLqmuteV5//b2vxezXJfCzIrNM1QLYmZuge8dab2WiUfImFk9h3tBdOpr9wgZM6vncC+ATn3tHiFjZo1yEe5LOzElO4tYg0597R4hYzZctm/fvjTW3TsxFVWn9drH14yx+6Y39a9CZpYb7dZzz8Wdu7XW6a79fW9/eZ9qYmZF4nDPuXYjZNzXbmatONxzrNNSA+5rN7NWHO451q5LZnzNWB9rYmZF43DPsXZdMu5rN7N2HO451alLxn3tZtaOwz2n3CVjZiuRi3D3JKYf12kdGXfJmA03T2IqqI0fuK9luHvSkpkt8SSmgvFdu5mtlMM9Z/wg1cyy4HDPmVt27m95zg9SzSwth3vOLDx9tOU5d8mYWVoO95wZUfPjwl0yZpaewz1HZuYWONFi8NLgxzSZWZE43HOk3cSldeNr+lgTMyu6XIS7JzF1nrh0/Ztf1sfamFmeeRJTgVx485daPkz1xCUza8aTmArAo2TMLEsO95zwKBkzy5LDPQc8SsbMspYq3CVdJGm/pHlJW5ucl6SPJ+cflnR+3bl3S9or6RFJn5b03CwbUAbtZqV6lIyZLUfHcJc0CtwKbAamgMskTTUU2wxMJq8KcFty7TrgD4DpiHgFMApcmlntS6Jdf7tHyZjZcqS5c78AmI+IRyPiGHAPsKWhzBbgzqjZBYxLOjM5twpYI2kVcBpwMKO6l4b7280sa2nCfR3wRN37A8mxjmUiYgH4EPA4cAg4HBH3Lb+65eP+djPrhTTh3uy+sjF3mpaR9AJqd/XnAmcBp0u6orHg4uIi09PTJ1/VajVFtcrBs1LNrBvVavVkVgJrW5VbleKzDgBn171fz6ldK63KvBH4TkQsAkj6PPCLwF31F09MTDCMk5g8K9XMulWpVJZmpyLpqVbl0ty53w9MSjpX0mpqD0RnG8rMAu9MRs1sotb9cohad8wmSadJEvAGYF/3zSmnTmu3u7/dzJar4517RDwj6VpgJ7XRLndExF5JVyfntwE7gIuBeeAIcGVy7muSPgd8HXgGmAOGp8+lA89KNbNeSdMtQ0TsoBbg9ce21X0dwDUtrr0JuGkFdSytEdH0YapHyZjZSnmG6oB4lIyZ9ZLDfUA8SsbMesnhPgAeJWNmveZwHwCPkjGzXstFuA/bTkwHPUrGzFbAOzHl1MYP3Ne0W+a0sRG+8eebB1AjMysi78SUMz88/mzT46tXjfa5JmZWVg73PpuZW+DI8RNNzx0+2vohq5lZNxzufdbuYepZHgJpZhlxuPeZN+Yws35wuPfZqJrvzOElB8wsSw73Pnu2xeikwY9ZMrMycbj30czcQtNdTcBLDphZtnIR7sMyien92/c2vUMX7m83s/Q8iSlHZuYWuO4zu1uef+zmt/SvMmZWCp7ElAPthkC6S8bMsuZw7xMPgTSzfnK494mHQJpZPznc+8RDIM2snxzufeAhkGbWbw73Prhl534PgTSzvnK490Grh6mB+9vNrDdyEe5ln8TU6mFqq+NmZu14ElNObNj6jy3PefKSmS2XJzENkB+mmtkgONx7zOvJmNkgONx7aGZuoelG2OCHqWbWW6nCXdJFkvZLmpe0tcl5Sfp4cv5hSefXnRuX9DlJ35S0T9JrsmxAnnk9GTMblI7hLmkUuBXYDEwBl0maaii2GZhMXhXgtrpzHwP+OSJ+Fvh5YF8G9S4ErydjZoOS5s79AmA+Ih6NiGPAPcCWhjJbgDujZhcwLulMSc8Hfhn4JEBEHIuIp7Orfr55PRkzG5Q04b4OeKLu/YHkWJoyLwYWgb+RNCfpE5JOX0F9C8XryZjZoKQJ92a3n4351KrMKuB84LaI2Aj8ADilz35xcZHp6emTr2q1mqJa+eYhkGbWC9Vq9WRWAmtblVuV4rMOAGfXvV8PHExZJoADEfG15PjnaBLuExMTlG0Sk9eTMbNeqFQqS7NTkfRUq3Jp7tzvByYlnStpNXApMNtQZhZ4ZzJqZhNwOCIORcR/A09IWkqzNwDf6LItheT1ZMxskDreuUfEM5KuBXYCo8AdEbFX0tXJ+W3ADuBiYB44AlxZ9xG/D9yd/GJ4tOFcaY1KTfvcvZ6MmfVDmm4ZImIHtQCvP7at7usArmlx7W6g6doHZdbqYWqr42ZmWfIM1R7ww1QzGzSHew/4YaqZDZrDvQf8MNXMBs3hnjF3yZhZHuQi3Mu0E5O7ZMys17wT0wB41yUz6xfvxNRH3i/VzPLA4Z4xj283szxwuGfID1PNLC8c7hnyfqlmlhcO94x4v1QzyxOHe0a8X6qZ5YnDPSMHvV+qmeVILsK9DJOYxk8ba3r8tLERd8mYWaY8iamPpv70nzhy/MQpx8fXjLH7pjcNoEZmVnaexNRjM3MLTYMd4PDR5g9Zzcx6yeGegXYPU8/yw1QzGwCHewb8MNXM8sbhngE/TDWzvHG4Z+CHx59tenz1qtE+18TMrMbhvkJ+mGpmeeRwXyE/TDWzPMpFuBd5ElOr/VLBD1PNrDc8ianHZuYWePdndjddCdKTl8ys1zyJqUfa7Zf6vre/vN/VMTM7yeG+Aq26ZLzEr5kNmsN9BbxfqpnlVapwl3SRpP2S5iVtbXJekj6enH9Y0vkN50clzUn6QlYVzwPvl2pmedUx3CWNArcCm4Ep4DJJUw3FNgOTyasC3NZw/l3AvhXXNke8X6qZ5VmaO/cLgPmIeDQijgH3AFsaymwB7oyaXcC4pDMBJK0H3gJ8IsN6D1y7h6keAmlmg5Ym3NcBT9S9P5AcS1vmo8AfA82ncRaUH6aaWZ6lCfdmvQ+NN61Ny0h6K/BkRDzY7hssLi4yPT198lWtVlNUa7BGWvTJ+GGqmfVStVo9mZXA2lblVqX4rAPA2XXv1wMHU5b5deDtki4Gngs8X9JdEXFF/cUTExMUaRLTzNwCJ1o8M/XDVDPrpUqlsjQ7FUlPtSqX5s79fmBS0rmSVgOXArMNZWaBdyajZjYBhyPiUES8NyLWR8SG5LovNQZ7EbVbT8YPU80sDzreuUfEM5KuBXYCo8AdEbFX0tXJ+W3ADuBiYB44AlzZuyoPnteTMbO8S9MtQ0TsoBbg9ce21X0dwDUdPuPLwJe7rmHOLA2BbLWejB+mmlkeeIZql7yejJkVgcO9Sx4CaWZF4HDvkteTMbMicLh3yevJmFkR5CLci7ITk9eTMbM88E5MGbvw5i817XMX8JHfOs997mbWV96JKSN+mGpmReFw74LXkzGzonC4p+T1ZMysSBzuKb1/+96W5/ww1czyxuGewszcAt87crzlea8nY2Z543BPod0qkF5PxszyyOGewsE2q0B6PRkzy6NchHveJzGNnzbW9PhpYyO+azezvvMkpoxM/ek/ceT4qVvAjq8ZY/dNbxpAjczMPIlpRWbmFpoGO8Dho60fspqZDZLDvYN2D1PP8hBIM8sph3sH3lLPzIrI4d5BqyUHhNeTMbP8cri30W7JgcE/hjYza83h3oaXHDCzonK4t+AlB8ysyHIR7nmcxOQlB8wsrzyJaQU2bP3Hluc+6l2XzCwHPIlpGTxKxsyKzOHehEfJmFnROdyb8CgZMyu6VOEu6SJJ+yXNS9ra5LwkfTw5/7Ck85PjZ0v6V0n7JO2V9K6sG5A1j5IxszLoGO6SRoFbgc3AFHCZpKmGYpuByeRVAW5Ljj8D/GFE/BywCbimybW54lEyZlYGae7cLwDmI+LRiDgG3ANsaSizBbgzanYB45LOjIhDEfF1gIj4X2AfkOt0bLeWjDfmMLOiSBPu64An6t4f4NSA7lhG0gZgI/C1rmvZRx4lY2ZlsCpFmWZx1zhopG0ZSc8D/h64LiK+31hwcXGR6ekfDdWsVCpLA/T7yqNkzCzvqtUq1Wp16e3aVuXShPsB4Oy69+uBg2nLSBqjFux3R8Tnm32DiYkJ8jCJqV1/u0fJmFke1N/8SnqqVbk03TL3A5OSzpW0GrgUmG0oMwu8Mxk1swk4HBGHJAn4JLAvIj68nIb0k9duN7Oy6HjnHhHPSLoW2AmMAndExF5JVyfntwE7gIuBeeAIcGVy+YXAbwN7JO1Ojv1JROzItBUZmJlbaHnO/e1mVjRpumVIwnhHw7FtdV8HcE2T6/6d5v3xudNu4pL7282saDxDlc4Tl9zfbmZF43Cn/YNU4f52MysehzvtH6Revukc97ebWeE43Gk/cemDl7yyr3UxM8tCLsJ9kDsxeeKSmRWNd2JKYeMH7mv5MHXd+Bq+svX1fa6RmVk63ompDS/va2ZlNNThfuPMnrbn/SDVzIpqaMN9Zm6Bu3Y93vL8+JqxPtbGzCxbQxvu7WakgtduN7NiG9pwb9fX7h2XzKzohjLc2y0SBr5rN7PiG8pwv+He1g9S14yN+K7dzAovF+Hez0lMM3ML/ODYsy3P/+WvvqrndTAzWwlPYmqi3aQlgMdufktf6mFmtlKexJTotLSvhz+aWVkMVbi362sHP0g1s/IYmnC/cWZP2772K7y0r5mVyFCEe6fZqOClfc2sXIYi3Dt1x7iv3czKpvTh3mnoI7iv3czKJxfh3stx7p3u2t3XbmZFM/Tj3C+//at85dvfbXl+zdgI+/58c+bf18ysH4ZynHunYAfPRjWz8ipluN84s6djsHsNGTMrs1WDrkDW0tyxg+/azazcShPuM3MLvOczuzmRoqwfoppZ2aUKd0kXAR8DRoFPRMTNDeeVnL8YOAL8bkR8Pc21yzUzt8D7Zvfy9NHWa8U0c+FLXugJS2ZWeh373CWNArcCm4Ep4DJJUw3FNgOTyasC3NbFtV2bmVvg+s8+1HWwj43A3Ve9ZqXf3sws99I8UL0AmI+IRyPiGHAPsKWhzBbgzqjZBYxLOjPltV27Zed+jp/obgjniOCW3zhvpd/azKwQ0oT7OuCJuvcHkmNpyqS5lsXFRaanp0++qtVq2wodfPpoimr/yHNWjfDh3zzP/exmVnjVavVkVgJrW5VL0+euJscab5tblUlzLRMTE3Qziems8TUspAz4yRedzhff89rUn21mlmeVSmVpdiqSnmpVLs2d+wHg7Lr364GDKcukubZr17/5ZYyNNPu98eOu2HSOg93MhlKaO/f7gUlJ5wILwKXAOxrKzALXSroH+AXgcEQckrSY4tquLXWvNBstI+DyTed4RIyZDbWO4R4Rz0i6FthJbTjjHRGxV9LVyfltwA5qwyDnqQ2FvLLdtVlU/JKN69yHbmbWQqkXDjMzK7OhXDjMzGyYFT7cOw2bLCO3eTi4zcOhV23OXbh3u2FHtz+Y5WwI0u01vS7vNmf/+cu5xm3OvrzbnN3n5yLc63di6sVuTPXK8D9Dt9zm3lzjNmdfvlvD2ubC7MSUDJn8r+TtGcDhLi5fC7QcyN9Et5+/nGt6Xd5tzv7zl3ON25x9ebe5u8//mYiYaFYoF+FuZmbZykW3jJmZZcvhbmZWQoUOd0kXSdovaV7S1kHXJyuS7pD0pKRH6o69UNIXJX0r+fMFdefem/wM9kt682BqvXySzpb0r5L2Sdor6V3J8TK3+bmS/kPSQ0mb358cL22bl0galTQn6QvJ+2Fo82OS9kjaLemB5Fhv2x0RhXxRW87g28CLgdXAQ8DUoOuVUdt+GTgfeKTu2F8DW5OvtwJ/lXw9lbT9OcC5yc9kdNBt6LK9ZwLnJ1//BPCfSbvK3GYBz0u+HgO+Bmwqc5vr2v4e4G+BLyTvh6HNjwFrG471tN1FvnPvyUYgeRAR/wY07vK9BfhU8vWngEvqjt8TET+MiO9QW9/ngn7UMysRcSiSbRkj4n+BfdTW/S9zmyMi/i95O5a8ghK3GUDSeuAtwCfqDpe6zW30tN1FDvdUG4GUyE9FxCGohSHwouR4qX4OkjYAG6ndyZa6zUn3xG7gSeCLEVH6NgMfBf4Yfmwv+7K3GWq/uO+T9KCkSnKsp+1OtUF2TqXaCGQIlObnIOl5wN8D10XE92v7rjcv2uRY4docEc8C50kaB+6V9Io2xQvfZklvBZ6MiAclvTbNJU2OFarNdS6MiIOSXgR8UdI325TNpN1FvnPvyUYgOfY/yb60JH8+mRwvxc9B0hi1YL87Ij6fHC51m5dExNPAl4GLKHebLwTeLukxat2or5d0F+VuMwARcTD580ngXmrdLD1td5HD/eQmIpJWU9sIZHbAdeqlWeB3kq9/B/iHuuOXSnpOsinKJPAfA6jfsql2i/5JYF9EfLjuVJnbPJHcsSNpDfBG4JuUuM0R8d6IWB8RG6j9ff1SRFxBidsMIOl0ST+x9DXwJuARet3uQT9FXuET6Iupjaz4NnDDoOuTYbs+DRwCjlP7Lf57wE8C/wJ8K/nzhXXlb0h+BvuBzYOu/zLa+0vU/tn5MLA7eV1c8ja/CphL2vwI8GfJ8dK2uaH9r+VHo2VK3WZqI/oeSl57l7Kq1+328gNmZiVU5G4ZMzNrweFuZlZCDnczsxJyuJuZlZDD3cyshBzuZmYl5HA3Myshh7uZWQn9P7AIqqsF2kYVAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.arange(n), ms.psi(ms.rho(ms.R, ms.m), ms.w * ms.rho(ms.R, ms.m) ) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 30.0)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD1CAYAAACrz7WZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARGElEQVR4nO3df6xkZXnA8e/TZY3rQljb3RZcUGh7g4nYsnCDmm2axdjyoxKosY38odY/HGswwcTQon/4o7FxW1vTGhLoEImSWIwR3bKooSRo1Cag7LKw4HbrpqVlf8S9SLmysFXAp3/cuTLcvXdnzuycM2fe/X6Syc6cec/M+56z9znveed9zonMRJJUll+ZdAUkSeNncJekAhncJalABndJKpDBXZIKZHCXpAKd0uSXnXrqqfna1752qLLz8/OcfvrpU1MWYG5ujg0bNky0HnWVrattVctPW/vasi3a0L46t8W0ta/qttixY8eRzDxt6BUAMrOxxyte8Yoc1nvf+96pKpuZedFFF028HnWVrattVctPW/vasi3a0L46t8W0ta/qtgCeyYrxtrXDMldeeeVUla2qDXWuq31VP7fk9rVlW7ShDnVui7o+tw1lRxXZYIbq2rVr85lnnmns+5o2OzvLAw88MOlq1KLktoHtm3alty8ins3MtVXWabTnvn79+ia/rnGdTmfSVahNyW0D2zftSm8fMFd1hUZ77jMzM3nJJZdw5ZVXNnJaIknTbPv27Wzfvp1bbrllX2bOVFm30eA+OzubJZ86SVIdImJHZs5WWafRqZC7D8yzeeu9XH/peVy9aeOyZbY9eIBP372Xg08d5VXr1rS+rCS1UaPBHeDAU0f58Fd3AxwTMLc9eIAPf3U3R597YSrK9q/jgUNSm0xkKuTR517g03fvPWb5p+/e+8ugOg1l4cWDwYGnjpK8eDDY9uCBEyq7WH7z1ns594avs3nrvSuWk6SlBgb3iHh5RHw/Ih6KiEcj4hPLlImI+GxE7IuIhyPiwkGfe/Cpo0Mta3NZqO/AUeeBwIOGVL5heu4/A96cmb8LXABcFhFvXFLmcmCm9+gANw360FetWzPUsjaXhfoOHHUdCDx7kE4OA4N7L/v1SO/l6t5j6RSbq4DbemXvA9ZFxJkrfeaa1au4/tLzjll+/aXnsWb1qqkpC/UdOOo6ELTl7EFSvYYac4+IVRGxCzgM3JOZ9y8pshF4vO/1/t6yY2xct4ZPve31y/6IePWmjXzqba9n47o1xBSUhfoOHHUdCNpw9tC/jgcDqR5DzZbJzBeACyJiHfC1iDg/Mx/pKxLLrbZ0wRkv+zk/+8pf8MmvwCdZyCpbmll29aaNQ88eaUtZYKgZMFXKXn/peS+ZtQPHPxAcWCY4r3T2MGzZcR0Ilmtf1VlJzjLSyabb7dLtdhdfVk7vrzQVMjOfiohvA5cB/cF9P3B23+uzgINL19+wYUOR13+o48BR14GgroPGOH+IHsf0VGna9Xd+I+KJqusPDO4RsQF4rhfY1wBvAf5mSbE7gQ9ExJeANwDzmXmoamX0UnUcCNpw9gD1nxXYy9fJbpie+5nAFyJiFQtj9F/OzLsi4s8BMvNm4BvAFcA+4FngPTXVVyuYprMHqO+swOEeacHA4J6ZDwOblll+c9/zBK4db9XUBnUcCKC+s4I6h3s8EGiaNH75AZWrrh+iqxwI6hrucdxf08bgromp46ygTbOB7OVrkgzumgrDHgjaMBuo7ovPScNo7T1UpVFUSUKrK6mszovPScNqtOc+Pz9Pp9PxTkyq1aRnA9U55x/s5Z9MFu/EBJxedV3vxCQNadigunnrvcsO92xct4Z/u+HNxyw/94avH5vOzULa939t/aNj6rDcQeZ4l8jQ9Gv9nZikaVbHuD/UN9UT7OWfzBxzl8aszovPjfLDrlf1PDnZc5dqUNec/7YkdKn9DO5SC9Qx5OM8/pObwV2aIm1I6LKXPx0M7tKUmXRCl7386eAPqlKh6kroqvNHXY2PSUxSwSZ9/Z6qvXywp9/PJCZJjamSSFUlQavqZ59MRkliclhGUiVVhnuqXJMHRrsuj3Pzl+cPqpIqqytb11k742PPXVJtqmbr1n31zZOpl2/PXVKtqmTr1pWkdTL28u25S2qNusbzT8Zevj13Sa0y6UsxlNLLt+cuaSq1pZffViYxSZpak+7lQ71JVyYxSdIAddxJq6mkK+/EJEkrqKOX3+aLqBncJalPlevstPmHWoO7JC0xbC+/qYuoveyM375oyKr/krNlJGlEdV0qGV56ueRRGNwlaURNX0StCodlJOkENH0RtWHZc5ekBozrImrDanSe+8zMTF5yySUmMUnSANsePMAHtt7C/N77+L//fpjn/vdgVFnfJCZJaqnF2TIP/OP7+NmhHxncJakk3mZPkgQY3CWpSAODe0ScHRHfiog9EfFoRFy3TJktETEfEbt6j4/WU11J0jCGmef+PPChzNwZEacBOyLinsz84ZJy383Mt46/ipKkqgb23DPzUGbu7D1/GtgDTM/tSCTpJFRpzD0izgE2Afcv8/abIuKhiPhmRLxuHJWTJI1m6MsPRMSpwB3ABzPzp0ve3gm8JjOPRMQVwDZgZulnzM3NMTv74myeTqdDp9MZpd6SVLRut0u32118ub7q+kPNc4+I1cBdwN2Z+Zkhyj8GzGbmE/3LnecuSdXVMs89IgL4HLBnpcAeEWf0yhERF/c+9ydVKiJJGp9hhmU2A+8EdkfErt6yjwCvBsjMm4G3A++PiOeBo8A7ssnUV0nSSwwM7pn5PeC41zTIzBuBG8dVKUnSiTFDVZIKZHCXpAIZ3CWpQAZ3SSpQo8F9fn6eTqfD9u3bm/xaSZpK27dvX0z0PL3qut6sQ5Jazpt1SJIAg7skFcngLkkFMrhLUoEM7pJUIIO7JBXI4C5JBTKJSZJayiQmSSqYSUySJMDgLklFMrhLUoEM7pJUIIO7JBXI4C5JBTK4S1KBTGKSpJYyiUmSCmYSkyQJMLhLUpEM7pJUIIO7JBXI4C5JBTK4S1KBDO6SVCCTmCSppUxikqSCmcQkSQIM7pJUJIO7JBVoYHCPiLMj4lsRsSciHo2I65YpExHx2YjYFxEPR8SF9VRXkjSMU4Yo8zzwoczcGRGnATsi4p7M/GFfmcuBmd7jDcBNvX8lSRMwsOeemYcyc2fv+dPAHmDjkmJXAbflgvuAdRFx5thrK0kaSqUx94g4B9gE3L/krY3A432v93PsAUCS1JBhhmUAiIhTgTuAD2bmT5e+vcwqx0ygn5ubY3b2xamanU5ncYK+JKlPt9ul2+0uvlxfdf2hkpgiYjVwF3B3Zn5mmff/Cfh2Zt7ee70X2JKZh/rLmcQkSdXVksQUEQF8DtizXGDvuRN4V2/WzBuB+aWBXZLUnGGGZTYD7wR2R8Su3rKPAK8GyMybgW8AVwD7gGeB94y9ppKkoQ0M7pn5PZYfU+8vk8C146qUJOnEmKEqSQUyuEtSgQzuklQgg7skFcg7MUlSS3knJkkqmHdikiQBBndJKpLBXZIKZHCXpAIZ3CWpQAZ3SSqQwV2SCmQSkyS1lElMklQwk5gkSYDBXZKKZHCXpAIZ3CWpQAZ3SSqQwV2SCmRwl6QCmcQkSS1lEpMkFcwkJkkSYHCXpCIZ3CWpQAZ3SSqQwV2SCmRwl6QCGdwlqUAmMUlSS5nEJEkFM4lJkgQY3CWpSAZ3SSrQwOAeEbdGxOGIeGSF97dExHxE7Oo9Pjr+akqSqjhliDKfB24EbjtOme9m5lvHUiNJ0gkb2HPPzO8ATzZQF0nSmIxrzP1NEfFQRHwzIl43ps+UJI1omGGZQXYCr8nMIxFxBbANmFmu4NzcHLOzL07V7HQ6ixP0JUl9ut0u3W538eX6qusPlcQUEecAd2Xm+UOUfQyYzcwnlr5nEpMkVTeRJKaIOCMiovf84t5n/uREP1eSNLqBwzIRcTuwBVgfEfuBjwGrATLzZuDtwPsj4nngKPCObPKaBpKkYwwM7pl5zYD3b2RhqqQkqSXMUJWkAhncJalABndJKpDBXZIK5J2YJKmlvBOTJBXMOzFJkgCDuyQVyeAuSQUyuEtSgQzuklQgg7skFcjgLkkFMolJklrKJCZJKphJTJIkwOAuSUUyuEtSgQzuklQgg7skFcjgLkkFMrhLUoFMYpKkljKJSZIKZhKTJAkwuEtSkQzuklQgg7skFcjgLkkFMrhLUoEM7pJUIJOYJKmlTGKSpIKZxCRJAgzuklQkg7skFWhgcI+IWyPicEQ8ssL7ERGfjYh9EfFwRFw4/mpKkqoYpuf+eeCy47x/OTDTe3SAm068WpKkEzEwuGfmd4Anj1PkKuC2XHAfsC4izhxXBSVJ1Y1jzH0j8Hjf6/29ZZKkCTllDJ8RyyxbdvL83Nwcs7MvTtXsdDqLE/QlSX263S7dbnfx5fqq648juO8Hzu57fRZwcLmCGzZswCQmSRqsv/MbEU9UXX8cwzJ3Au/qzZp5IzCfmYfG8LmSpBEN7LlHxO3AFmB9ROwHPgasBsjMm4FvAFcA+4BngffUVVlJ0nAGBvfMvGbA+wlcO7YaSZJOmBmqklQgg7skFcjgLkkFMrhLUoG8E5MktZR3YpKkgnknJkkSYHCXpCIZ3CWpQAZ3SSqQwV2SCmRwl6QCGdwlqUAmMUlSS5nEJEkFa30S09zcXJNf17i++x0Wp+S2ge2bdqW3jxHuodpocH/iieFvA1hl6KYNZaHaf7A21LlK2braVrX8tLWvLduiDe2rc1tMW/tGGJreUHWF1v6g2oaNWudvA22oc13ta0tAa0P72rIt2lCHOrdFXZ/bhrKjanTMPSKeBvYOWfx0YH6KysLCqdOwpydtqHOVsnW1rWr5aWtfW7ZFG9pX57aYtvZV3RbnZeZpFco3G9wlSc1o7bCMJGl0BndJKlAjwT0iLouIvRGxLyJuaOI7mxQRj0XE7ojYFRFTP5E/Im6NiMMR8Ujfsl+NiHsi4ke9f185yTqeiBXa9/GIONDbh7si4opJ1nFUEXF2RHwrIvZExKMRcV1veRH77zjtK2X/vTwivh8RD/Xa94ne8sr7r/Yx94hYBfwH8AfAfuAHwDWZ+cNav7hBEfEYMJuZw8/1bLGI+H3gCHBbZp7fW/a3wJOZubV3gH5lZv7lJOs5qhXa93HgSGb+3STrdqIi4kzgzMzcGRGnATuAq4E/o4D9d5z2/Sll7L8A1mbmkYhYDXwPuA54GxX3XxM994uBfZn5n5n5c+BLwFUNfK9GlJnfAZ5csvgq4Au9519g4Q9qKq3QviJk5qHM3Nl7/jSwB9hIIfvvOO0rQi440nu5uvdIRth/TQT3jcDjfa/3U9DO6EngXyNiR0R0Jl2ZmvxGZh6ChT8w4NcnXJ86fCAiHu4N20zlsEW/iDgH2ATcT4H7b0n7oJD9FxGrImIXcBi4JzNH2n9NBPdYZllp8y83Z+aFwOXAtb3Tfk2Xm4DfAi4ADgF/P9HanKCIOBW4A/hgZv500vUZt2XaV8z+y8wXMvMC4Czg4og4f5TPaSK47wfO7nt9FnCwge9tTGYe7P17GPgaC0NRpflxb7xzcdzz8ITrM1aZ+ePeH9UvgFuY4n3YG6u9A/hiZn61t7iY/bdc+0raf4sy8yng28BljLD/mgjuPwBmIuLciHgZ8A7gzga+txERsbb3ww4RsRb4Q+CR4681le4E3t17/m7gXyZYl7Fb/MPp+WOmdB/2fpD7HLAnMz/T91YR+2+l9hW0/zZExLre8zXAW4B/Z4T910iGam9a0j8Aq4BbM/Ova//ShkTEb7LQWwc4BfjnaW9fRNwObGEhpfvHwMeAbcCXgVcD/wP8SWZO5Y+SK7RvCwun9Ak8BrxvcYxzmkTE7wHfBXYDv+gt/ggL49JTv/+O075rKGP//Q4LP5iuYqHz/eXM/KuI+DUq7j8vPyBJBTJDVZIKZHCXpAIZ3CWpQAZ3SSqQwV2SCmRwl6QCGdwlqUAGd0kq0P8DiwbFA//+ZloAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.arange(n), ms.rho(ms.R, ms.m))\n",
    "plt.xlim(0,30)\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 10.0)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD7CAYAAACCEpQdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWj0lEQVR4nO3df3Bd5Z3f8fdXV/4l/7YljH8FwwIujoGQeKHdZJNQIGuyCJhOt1toM7s7NJp0hp3++KNLOp3ZznQ6TTudTpsubUYhlGQ2A6XszgY7ZByabZYlQ3YJBBYbL1hrCAgZLGFbBtsgS/fbP3RlXQvZvvKVdK983q8Zzb3nuec559Gx9bnPfc5zz4nMRJJULC2NboAkafYZ/pJUQIa/JBWQ4S9JBWT4S1IBGf6SVECGvyQVkOEvSQU0q+EfEYsj4rmIuG029ytJOl1d4R8RD0bEwYjYPaF8e0S8EhE9EXFf1Uu/Bzxazz4lSfWLei7vEBGfBd4HvpOZWytlJeBV4BagF3gWuAtYB7QDC4GBzNxZX9MlSeertZ7KmflURGyaUHw90JOZ+wEi4hHgDmAJsBjYApyIiCcys1zP/iVJ56eu8D+D9cCbVcu9wA2ZeS9ARPw2oz3/jwT/woULs1QqnVpub2+no6NjSjsfHBxk+fLl59Hs6dvGdLShv79/yr/7TLTDYzF92/BYTG8bin4s+vv7GRgYAOD48eMfZubCKW0gM+v6ATYBu6uWfwN4oGr5S8B/r2VbbW1tWa8vf/nLDd/GdLThU5/6VN3b8FhMbzs8FtO3DY/F9LYBOJZTzO6ZmO3TC2ysWt4A9NVScWRkhK6uLnbs2HHeO+/s7DzvutO1jelow3TwWIzzWIzzWIyb68dix44ddHV1AZTOte5EdZ3wBaiM+e/M8RO+rYye8L0JeIvRE753Z+aec21r8eLFeezYsbrac6HYtm0bP/vZzxrdjKbgsRjnsRjnsRgXEcczc/FU6tQ71fNh4Blgc0T0RsQ9mTkM3AvsAvYCj9YS/ACLFy+uu+d/oai8mwuPRTWPxTiPxWk9/yn3muvu+U+nbdu2pe/kkjQ1EfFcZm6bSh0v7yBJBdRU4T84OOiwjyTVqGrYZ8pzRR32kaQ5zmEfSVJNmir8HfaRpNo57CNJBeawjySpJk0V/g77SFLtHPaRpAJz2EeSVBPDX5IKqKnC3zF/SaqdY/6SVGCO+UuSamL4S1IBGf6SVECGvyQVUFOFv7N9JKl2zvaRpAJzto8kqSaGvyQVkOEvSQVk+EtSARn+klRATRX+TvWUpNo51VOSCsypnpKkmhj+klRAhr8kFZDhL0kFZPhLUgEZ/pJUQIa/JBXQrIV/RFwVEd+IiMci4p/O1n4lSR9VV/hHxIMRcTAidk8o3x4Rr0RET0TcB5CZezPzK8A/AKb0ZQRJ0vSqt+f/ELC9uiAiSsD9wK3AFuCuiNhSee124GngR3XuV5JUh7rCPzOfAg5NKL4e6MnM/Zk5BDwC3FFZ//HM/BXgH9WzX0lSfVpnYJvrgTerlnuBGyLi88DfAxYAT0xWsb+/n23bxkeEurq6xi5aJEmq0t3dTXd399hi+1Trz0T4xyRlmZk/Bn58toodHR14YTdJOrfqznFEDEy1/kzM9ukFNlYtbwD6aqnoJZ0lqXYNvaRzRGwCdmbm1spyK/AqcBPwFvAscHdm7jnXtryksyRN3axf0jkiHgaeATZHRG9E3JOZw8C9wC5gL/BoLcEvSZo9dY35Z+ZdZyh/gjOc1D2bsWGfzs5OOjs762maJF3wduzYMTZM7p28JKlo5vydvDzhK0m18x6+klRgc77nL0maHU0V/g77SFLtHPaRpAJz2EeSVJOmCn+HfSSpdg77SFKBOewjSaqJ4S9JBWT4S1IBNVX4e8JXkmrnCV9JKjBP+EqSamL4S1IBGf6SVEBNFf6e8JWk2nnCV5IKzBO+kqSaGP6SVECGvyQVkOEvSQVk+EtSARn+klRATRX+zvOXpNo5z1+SCsx5/pKkmhj+klRAhr8kFZDhL0kFZPhLUgEZ/pJUQIa/JBXQrIV/RNwZEd+MiO9FxBdma7+SpI+qK/wj4sGIOBgRuyeUb4+IVyKiJyLuA8jMP8nMLwO/DfxmPfuVJNWn3p7/Q8D26oKIKAH3A7cCW4C7ImJL1Sr/pvK6JKlB6gr/zHwKODSh+HqgJzP3Z+YQ8AhwR4z6j8APMvP5evYrSapP6wxscz3wZtVyL3AD8LvAzcDyiLg8M78xsWJ/fz/bto1fnqKrq2vsokWSpCrd3d10d3ePLbZPtf5MhH9MUpaZ+XXg62er2NHRgRd2k6Rzq+4cR8TAVOvPxGyfXmBj1fIGoK+Wil7SWZJq19BLOkfEJmBnZm6tLLcCrwI3AW8BzwJ3Z+aec23LSzpL0tTN+iWdI+Jh4Blgc0T0RsQ9mTkM3AvsAvYCj9YS/GDPX5Kmwpu5SFKBeTMXSVJNmir8HfaRpNo57CNJBTbnh33s+UtS7ez5S1KBzfmevyRpdhj+klRATRX+jvlLUu0c85ekAnPMX5JUE8NfkgqoqcLfMX9Jqp1j/pJUYI75S5JqYvhLUgEZ/pJUQIa/JBVQU4W/s30kqXbO9pGkAnO2jySpJoa/JBWQ4S9JBWT4S1IBGf6SVEBNFf5O9ZSk2jnVU5IKzKmekqSaGP6SVECGvyQVkOEvSQVk+EtSARn+klRAhr8kFdCshX9EXBYR34qIx2Zrn5KkydUV/hHxYEQcjIjdE8q3R8QrEdETEfcBZOb+zLynnv1JkqZHvT3/h4Dt1QURUQLuB24FtgB3RcSWOvcjSZpGdYV/Zj4FHJpQfD3QU+npDwGPAHfUsx9J0vSaiTH/9cCbVcu9wPqIWB0R3wCui4ivTlaxv7+fbdu2nfrp7u6egeZJ0tzX3d19KiuB9qnWb53+JhGTlGVmvgt85WwVOzo68MJuknRuXV1dY1f0JCIGplp/Jnr+vcDGquUNQF8tFb2ksyTVrqGXdI6ITcDOzNxaWW4FXgVuAt4CngXuzsw959qWl3SWpKmb9Us6R8TDwDPA5ojojYh7MnMYuBfYBewFHq0l+CVJs6euMf/MvOsM5U8AT0x1e2PDPp2dnXR2dtbTNEm64O3YsWNsmNw7eUlS0cz5O3l5wleSauc9fCWpwOZ8z1+SNDuaKvwd9pGk2jnsI0kF5rCPJKkmTRX+DvtIUu0c9pGkAnPYR5JUE8NfkgrI8JekAmqq8PeEryTVzhO+klRgnvCVJNXE8JekApqJG7ift919g1z9b3cxv9RCaymYV2phfmsLyxbOY/mieSxbNI/li1pZsWg+Fy1bwMXLFrJ2+SIuXr6Q1Yvn09Iy2b3jJUkTNVX4L2mFlUde4WObLmPDxo8xNFLmw+EyR0+c5MiJk7xx6DiDJ04yeOIkI+XTz1XMKwUbVrbxSx2LuaxjCZe1jz5euWYJK9rmN+g3kqSZU7g7eZXLybvHhnh78AMODJ7g7aMf0HfkA37x7jH29x/jtXePMTRcPrX++hWL2Lp+GVvXLWfr+tGfjqULZvJXkaRZcz4nfJuq51+rlpagY+kCOpYu4OoNH33DGyknfUdO0NP/Pq++/R4vvTXInr6j7Nrzzql1Nq1u45c3reKXL13F9ZtWccnqNiIcNpJUDHMy/M+l1BJsXNXGxlVt3Lj5olPl731wkpf7jvLCm0d49vXDPLn3Hf7Pc70AXLR0AZ+5vJ3Pbe7gM5e3s3qJnwwkXbjm5LDPdCmXk7/pf5+/fP0QP91/iKf39XP4+EkiYOu65Xz2ynY+d+VFfPJjK2gtOTFKUnM6n2GfQof/ROVysrtvkD97pZ+n9vXz/BtHGCknK9vmcfNVa9i+9WI+fXk7C+eVGtZGSZrI8J9mRz84ydP7Bvjhnrf50d6DvPfhMIvnl/j837qIX/v4xdy4uYOlC+c1upmSCq4wJ3xny7KF8/ji1Wv54tVrGRou88z+d9m1521+uOcdvv9XB5hfauFzmzu4/dp13HzVGhbN9xOBpLmhqXr+V1xxRd544410dnbS2dnZ6Oac0Ug5+fkbh3nipbf5/kt9vHP0Q9rml/jCljXc/ol1fObyDua3eo5A0swam+f/zW9+syczr5hK3aYK/2Yb9qnFSDn5y9cO8fiLffxg9wGOHD/JirZ53Lp1Lbdfu44bLl3lN48lzSjH/BtsaLjM0z39fO+FPp58+R2OD42wZtkCbrtmHbdfu45rNiz3uwSSpp3h30SODw3zo70HefzFPv7slX6GRspcsrqNzmvWcfsn1nHlmqWNbqKkC4Th36QGT5xk15632fFiHz/pGaCcsHnNUm7/xDo6r1nHx1a3NbqJkuYww38O6H/vQ36w+wCPv9DHz35xGIBrN66g85q1dF67jjXLFja4hZLmGsN/jnnryAl2vtjH4y/2safvKBFw/aZVbN96MbdsWcOGlX4ikHRuhv8c9jf977PjxT6+/1cH2HfwfQC2rF3GFz6+hlu2rGHL2mWeLJY0KcP/AvHawDGefHn0y2TPvXGYzNHLUt+yZQ2fu7KDGy5bRdt8v58naVRTh39ELAb+BzAE/DgzvztxHcP/o/rf+5A//et3ePLld/jzfQN8OFxmXin41CUr+dUrOvjVK9r5+LrllPwugVRYsx7+EfEgcBtwMDO3VpVvB/4bUAIeyMyvRcSXgCOZuSMi/ndm/ubE7Rn+Z/fByRGeff0QT+8b4Kl9A+w9cBSAFW3z2HbJSj55yUq2XbKKazYs9+JzUoE04to+DwF/AHynqhEl4H7gFqAXeDYiHgc2AC9VVhupc7+FtHBeqdLb7+CrjH4q+EnPAD/pGeC5Nw7zf/ceBKC1Jfj4+uVsu2QlW9cv46q1y/iljiXM87LUkirqCv/MfCoiNk0ovh7oycz9ABHxCHAHo28EG4AXAFNoGnQsXcCd163nzuvWA3Do2BDP/+Iwz71xmOdeP8wf/vQXfFi5neX8UgtXrFnCVWtH3wwuv2gJm1a3sX7FIu9VIBXQTJw1XA+8WbXcC9wAfB34g4j4dWDHZBX7+/vZtm38k0tXVxddXV0z0MQL06rF87l5yxpu3rIGgOGRMvsHjrH3wFFe7jvKyweO8uNXDvJY5e5lMHrj+40r27hkdRub2hezfsUiLl6+kLXLF7Jm2UIuWrrQi9RJTai7u5vu7u6xxfap1p+J8J/szGNm5jHgd85WsaOjA8f8p09rqYUr1yzlyjVLueMT60+V97/3Ia8NHOP1gdGb3b8+cIzXBo7x0/2HOHHyoyNy7UsW0L5kPiva5rGybfRxRdt8VlYely5oZdH8Em3zW2mbX6o8L9E2b7R8XilmfZpqZjJSTobLSTlHH0dGKo/lZLhcHn+9XF2ejJTLjJQ54zojVeuWT6tz+nZOqzPWnpGx9pRPLY/k5HXKmZRz9HfJhPKEx2T89VPrUb3ehLpVr52qQ1IuT1636mhWHdfJSseP+WSvnVbnDOtM3Njp9WvY7hnaeK79zG0b4df+3ejT5744MNXaMxH+vcDGquUNQF8tFQcHB+nq6mr6SzrPdR1LF9CxdAHXX7rqtPLMZPDESQ4MfsDbRz/gnbHHox8w8P4Qg8dP0nPwfQ4fP8mR40MMl2v/Syq1BK1jP6WWymPQ2tJCaykoRZz6uxwLotHn43/Yo6E1vs2xAJ0YumNlzaIloLWlhVJLnPpprXo+ttxyqryFUgu0RNASQcTo84Dx5RZoiZbx1069DhFx6vG0OpXHU69TVTcqdQlaWgBGy8ZUv3WfXn76m/qZ60z+5n9+25p8u6fVnrC76m1fSF+X6dnXQ0/PPvbB8qnWrXuqZ2XMf+fYbJ+IaAVeBW4C3gKeBe7OzD3n2pazfeaOzOTY0AiHjw3x/ofDHB8a4cTQCMeHhjlxcoTjQyOVsmFOjoz2dIfHerrl5ORIufI42mseLuepsILRP9Dx51V/ujH+h9zaEpRK40E6FpytkwZq5bHUMmH9qnUjmFdqOa3OZGF9KqBjfP8tEae159Ryy+x/6lHxzPpsn4h4GPg80B4RvcDvZ+a3IuJeYBejUz0frCX4wZ7/XBIRLFnQypIFftlMapSxm7nQiJ7/dLLnL0lTdz49f6dxSFIBNVX4jw37VD7GSJLOYseOHWPT4R32kaSimfPDPvb8Jal29vwlqcDmfM9fkjQ7DH9JKqCmCn/H/CWpdo75S1KBOeYvSaqJ4S9JBdRU4e+YvyTVzjF/SSowx/wlSTUx/CWpgAx/SSogw1+SCqipwt/ZPpJUO2f7SFKBOdtHklQTw1+SCsjwl6QCMvwlqYAMf0kqoKYKf6d6SlLtnOopSQXmVE9JUk0Mf0kqIMNfkgrI8JekAjL8JamADH9JKiDDX5IKaNbCPyIui4hvRcRjs7VPSdLkagr/iHgwIg5GxO4J5dsj4pWI6ImI+862jczcn5n3nG2d/v7+WppTCN3d3Y1uQtPwWIzzWIzzWJymfaoVau35PwRsry6IiBJwP3ArsAW4KyK2RMTVEbFzws9FtexkYGBgCk2f3HRcGqLebUxHG6bjP7bHYnrb4bGYvm14LKa3DUDHVCvUFP6Z+RRwaELx9UBPpUc/BDwC3JGZL2XmbRN+Dk61YefrAvrHrJvHYpzHYpzHYlyRj0XN1/aJiE3AzszcWln++8D2zPwnleUvATdk5r1nqL8a+PfALcADmfkfJlnnA2CkqqgfmOrHgeXA4BTrTPc2pqMN7Uz9d5+Jdngspm8bHovpbUPRj0U74z3+UmYunErl1vPY4ZiYpOyM7ySZ+S7wlbNtcKqNlySdn3pm+/QCG6uWNwB99TVHkjQb6gn/Z4ErIuLSiJgP/EPg8elpliRpJtU61fNh4Blgc0T0RsQ9mTkM3AvsAvYCj2bmnvNtyFSmjV7IImJjRPy/iNgbEXsi4p81uk2NFBGliPh5ROxsdFsaLSJWRMRjEfHXlf8ff6fRbWqEiPgXlb+N3RHxcEQUarh4sqn3EbEqIp6MiH2Vx5Xn3E4z3MylMm30VUZPBvcy+qnirsx8uaENa4CIWAuszcznI2Ip8BxwZxGPBUBE/EtgG7AsM29rdHsaKSK+Dfx5Zj5Q+bTdlplHGtysWRUR64GngS2ZeSIiHgWeyMyHGtuy2RMRnwXeB75TNQHnPwGHMvNrlc7zysz8vbNtp1ku7zDptNEGt6khMvNAZj5fef4eo5+q1je2VY0RERuAXwceaHRbGi0ilgGfBb4FkJlDRQv+Kq3AoohoBdoo2LnGM0y9vwP4duX5t4E7z7WdZgn/9cCbVcu9FDTwqlWm114H/EWDm9Io/xX4V0C5we1oBpcxOvX5f1WGwR6IiMWNbtRsy8y3gP8MvAEcAAYz84eNbVVTWJOZB2C0Awmc84u1zRL+U5o2WgQRsQT4I+CfZ+bRRrdntkXEbcDBzHyu0W1pEq3AJ4H/mZnXAceAwp0bq4xl3wFcCqwDFkfEP25sq+amZgl/p41WiYh5jAb/dzPzjxvdngb5NHB7RLzO6DDg342IP2xskxqqF+jNzLFPgY8x+mZQNDcDr2Vmf2aeBP4Y+JUGt6kZvFM5Xzh23vCcV1VolvB32mhFRASj47p7M/O/NLo9jZKZX83MDZm5idH/D3+amYXt4WXm28CbEbG5UnQTUMRJAG8Afzsi2ip/Kzcxel6s6B4Hfqvy/LeA752rQj3f8J02mTkcEWPTRkvAg/VMG53jPg18CXgpIl6olP3rzHyicU1Sk/hd4LuVDtJ+4Hca3J5Zl5l/Ubks/PPAMPBzoFCX96xMvf880B4RvcDvA18DHo2Iexh9g/yNc26nGaZ6SpJmV7MM+0iSZpHhL0kFZPhLUgEZ/pJUQIa/JBWQ4S9JBWT4S1IBGf6SVED/H27YWCjiJfjQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Abar, ms.rho(ms.R, ms.m) )\n",
    "plt.yscale('log')\n",
    "plt.ylim(10**-1,10**4)\n",
    "plt.xlim(0,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb9805062e0>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD1CAYAAABJE67gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbkklEQVR4nO3df3xU9Z3v8dcnv4FAIJCAEH4T+SEKlRgFagv+xNrIbrWuttVtuzUPb/Vue7e3u+6228dtd3vr1etudx+10mjV7e2urm7dCpaKrkrRROWHgPIrOAlIQgIkAQJE8mvme//IgEkMMgyTOXNm3s/HYx5mZr4n55Mv8ubkO59zjjnnEBER/0nzugAREYmOAlxExKcU4CIiPqUAFxHxKQW4iIhPKcBFRHwqI547GzNmjBs1ahR5eXnnvG1TUxMFBQXnvF1ra+s57y+abeK9r2jmI54/V7TbJet8xHvuNR99+X0+Nm3a1Oyc+/gP4JyL22PBggXurrvuctFYsGBBVNtFs79oa4znvqKZj3j+XNFul6zzEe+513z05ff5ADa6ATI17ksoZWVlCb+/aGuM576iEc+fK9rtknU+4j338dyX5uP8t4t2X+bieCZmSUmJ27hxY7TbEu22yUjz0Zfmoy/NR19+nw8z2+ScK+n/um8+xCwvL/e6hISi+ehL89GX5qOvZJ2PuB6BFxcXu6VLl1JWVhb3pRQREb9ZtWoVq1at4tFHHw0454r7v++bJRQRkVTl+yUUERHpK6594CLSV+/fgHVl5+RmBmYW0++pAJe4C4YcRz7s5HBbJ8fbuzjREeREezdtHd2cCD/aOro52RWkKxiioztEZ/jRFQzRGfzoeWfQEQyFCDkIOUco5Ai5nn045wi6nufOOYLh93rG9LzXOzT75Gef188csr2f9gnjPq+fz2xJstj1d8vIyUyP6fdUgEtMdQdD1B85GX58ePq/ja3ttLT1hPaRDzvPGmpZGWkMzUonMz2NrPQ0sjPSyMpI63me0fPa0KEZZKankZFmpKcZZpBmH32dbkaaGWlpRtpA74W/7s346IXe7/UeFs02vd/ofwzWd5vYHqFJ4shIi/2frQJcotbW0c2WuqNs299K9YHj7DpwnMChE3QGQ6fHpKcZF+TlcEFeDsWFueQPy2J0bjZjcrMYNTSLEUMyyc1OJzc7k2HZ6eRmZzAsuyeYReSTKcAlYu1dQd6qbWHd7mY2fnCY7Q3HCIZ6DqXHjchh5rjhXFk8humFuUzKH0rRqCGMG5FDhsJYZFDENcBbW1spLy9XH7iPnOwM8tKOA7zwbiNvvN/Mya4g2RlpfGrSSL65ZDolU/KZV5THyKFZXpcqknRO9YEDA17pSn3gMqAtdUf517c+YPV7jbR1Bhmfl8M1c8Zy1axCrpg2OuYfxojImZ2pD1xLKHJaKOR4rfoQv1hXy/o9hxmWlc7nLr6AmxcUUToln7RB+BBGRKKnABcAKgPN/O/VO9necIzxeTl8/8bZ3FY6idxs/S8ikqj0tzPF1R3+kB88v43XqpuYMHIID31xHjfNH68uEBEfUICnqGDI8WTVXv7vmmrSDP7mc7O4c+EUrW2L+IgCPAU1He/gW09vpqqmhaUzC/j7P76YCSOHeF2WiJwjBXiKebu2hXuf2szx9i4euPkSvlhSFPPrM4hIfMR1ofNUH3i4r1Hi7Ll36vnyY28zPCeD396zmFsvm6jwFklgq1atOnUzCvWBpyrnHD9fW8ODa6pZNH00K+5YwIicTK/LEpEIqQ88RTnneGBNNY+srWH5/PE8eMs8sjLUYSKSDBTgSe6fXwnwyNoabi+dxI//aK5OxhFJIjoUS2KPrqvlH/9rN7csKFJ4iyShiALczJaZWbWZBczsvgHeH2Vm/2lm75rZejObG/tS5Vy8uK2RH6/eyY0XX8D/ufkShbdIEjprgJtZOvAwcAMwB7jdzOb0G/Y3wBbn3CXAncA/xbpQidy2/a38j3/fyvyJI3no1nmkK7xFklIkR+ClQMA5V+uc6wSeBpb3GzMHeAXAObcLmGJmY2NaqUSk+UQHd/1qI6OGZlJx5wKdWSmSxCIJ8AlAXa/n9eHXetsKfAHAzEqByUBR/2/U1NRESUnJ6UdFRUV0VcuAnHN899mttLR1UnFnCYXDc7wuSUSiVFFRcTorgTEDjYmkC2Wg37/7N4/fD/yTmW0B3gM2A939NyooKEB94IPnyaq9vFbdxA9vuoi5Ewbs+xcRnygvLz91Eg9m1jzQmEgCvB6Y2Ot5EdDQe4Bz7hjwtfCODNgTfkic7Gg4xk9W7+LqWYXcuXCy1+WISBxEsoSyASg2s6lmlgXcBqzsPcDMRobfA/gGsC4c6hIHXcEQ33l2K3lDM3nglkt0erxIijjrEbhzrtvM7gXWAOnA48657WZ2d/j9FcBs4FdmFgR2AH82iDVLP09U7mFn4zFWfOVSRudme12OiMRJRGdiOudWA6v7vbai19dvAsWxLU0iUXf4Q/7x5fe5ZnYh1180zutyRCSOdCamjznn+Nvnt2EGP1w+V0snIilGAe5jr1UfYm11E39x7YW6IYNICtL1wH2qOxjiJ6t3MWX0UO5cOMXrckRkEOh64Enq6fX7uO+593jky5dyw8UXeF2OiAyiM10PXEsoPvRhZzf/8PJuLp00kmVz9cGlSKpSgPvQE5V7OXS8g+/dOFsfXIqkMAW4z7R1dPPY67UsnVnAgsn5XpcjIh5SgPvMv729jyMfdnHvVWq7F0l1CnAfae8KUvF6LYtnjGbB5FFelyMiHlOA+8gzG+toOt7BvUt19C0i6gP3ja5giF/8oZaSyaO4YprWvkVSgfrAk8Tv3m3knn97h0fvLOHaObrZkUgqUR+4zz1ZtYdJ+UO5alah16WISIJQgPvAtv2tbNh7hDsXTtYNikXkNAW4DzxRuZehWel8sWTi2QeLSMpQgCe45hMdrNrawM2XFpE3JNPrckQkgSjAE9wzG+voDIb400W6z6WI9KUAT2DOOZ7dWE/p1HxmFA73uhwRSTDqA09g6/ccZk9zG7dq7VskJakP3Mf+4pktvLT9IOu/dzVDsyK6famIJCH1gfvM8fYuVr/XSNm88QpvERmQAjxBrdraSHtXiFtLirwuRUQSlAI8QT2zsY4Lx+Yyf+JIr0sRkQSlAE9AH7S0saXuKF+4tEh33BGRM1KAJ6CVWxoAKJs33uNKRCSRKcATjHOOlVsbKJ2Sz4SRQ7wuR0QSmPrAE8yuA8d5/9AJyubr6Fsk1akP3Gfu//0uHnu9lvXfu4b8YVlelyMiCUB94D7gnGPV1gY+XTxG4S0iZ6UATyDv7DvK/qMnuUkfXopIBCIKcDNbZmbVZhYws/sGeD/PzFaZ2VYz225mX4t9qcnvxW2NZKWn6ZZpIhKRswa4maUDDwM3AHOA281sTr9h9wA7nHPzgCXAQ2amNYBz4Jzjxe0HWDxjNMNzdN1vETm7SI7AS4GAc67WOdcJPA0s7zfGAcOt56yTXOAw0B3TSpPcjsZj1B0+ybK547wuRUR8IpIAnwDU9XpeH36tt58Bs4EG4D3gW865UEwqTBFrth0gzeCa2Vo+EZHIRBLgA53L3b/38HpgCzAemA/8zMxG9N+oqamJkpKS04+KiopzLDd5vbj9AKVT8xmdm+11KSKSACoqKk5nJTBmoDGRXKe0Huh9R4Eieo60e/sacL/raSoPmNkeYBawvveggoIC1Af+cTVNJ9h98AT/q6z/RwsikqrKy8tPncSDmTUPNCaSI/ANQLGZTQ1/MHkbsLLfmH3A1eEdjQVmArVR1p1y1mw/AMB1F2n9W0Qid9YjcOdct5ndC6wB0oHHnXPbzezu8PsrgL8DnjSz9+hZcvkr59yA/2LIx7284yCXFOUxXtc+EZFzENGtXpxzq4HV/V5b0evrBuC62JaWGlpOdLCl7ijfvvpCr0sREZ/RmZgeW1vdhHNw1axCr0sREZ9RgHvs1V2HKByezUXjP9a0IyLyiRTgHuoKhli3u4mlMwtJS9Odd0Tk3Oh64B7auPcIxzu6WarlExEZgK4HnsB+/Lsd/EvVB7zzg2vJzY7o82QRSUG6HngCenXXIS6flq/wFpGoKMA98kFLGzVNbSydqeUTEYmOAtwjr+46BMDVsxXgIhIdBbhHXqtuYlrBMCaPHuZ1KSLiUwpwD3R0B1m/p4XPFBd4XYqI+JgC3AOb9x2lvSvEoumjvS5FRHxMfeAeqAo0k2Zw+TQFuIicmfrAE9Atj1TRFXI8f89ir0sRER9QH3iCaOvoZkvdURZr+UREzpMCPM7W7zlMd8ixaPqAd0gSEYmYAjzOKgPNZGWkUTJllNeliIjPKcDjrLKmhQWTRpGTme51KSLicwrwODrc1snOxmMsnqH1bxE5fwrwOHqzpgWAhVr/FpEYUB94HFXWNJObncG8ogFbOkVE+lAfeAJZ8uBrTC/I5ZdfvczrUkTER9QH7rH9R0+yt+VDFs3Q8omIxIYCPE6qAs0A+gBTRGJGAR4nVTUtjB6WxYWFw70uRUSShAI8DpxzVAaaWTh9tO4+LyIxowCPg5qmNg4d72Cx1r9FJIYU4HFQVRNe/1b/t4jEkPrA46Ay0MyEkUOYmD/E61JExEfUB+6xYMjxqR+9xLK543jglnlelyMiPqQ+cI/saDjGsfZurX+LSMwpwAdZZXj9e6Fu4CAiMRZRgJvZMjOrNrOAmd03wPvfNbMt4cc2MwuaWX7sy/WfykAzxYW5FA7P8boUEUkyZw1wM0sHHgZuAOYAt5vZnN5jnHMPOufmO+fmA38N/ME5d3gQ6vWVju4gG/Ye1vKJiAyKSI7AS4GAc67WOdcJPA0s/4TxtwNPxaI4v9uy7yjtXSEWaflERAZBJAE+Aajr9bw+/NrHmNlQYBnwm/Mvzf8qa1pIM7h8mgJcRGIvI4IxA537fabewzKg8kzLJ01NTZSUfNQJU15efqrHMSlVBZq5eEIeeUMyvS5FRHymoqKCioqKU08HXIeNJMDrgYm9nhcBDWcYexufsHxSUFBAqvSBt3V0s6XuKHd9ZprXpYiID/U+wDWz5oHGRLKEsgEoNrOpZpZFT0iv7D/IzPKAzwLPR11xElm/9zDdIafT50Vk0Jz1CNw5121m9wJrgHTgcefcdjO7O/z+ivDQPwZecs61DVq1PlIVaCYrI42SKaO8LkVEklQkSyg451YDq/u9tqLf8yeBJ2NVmN9VBlpYMGkUOZnpXpciIklKZ2IOgsNtnexoPKb2QREZVArwQfBWbQuA7n8pIoNKAT4IKgPN5GZnMK9owCtAiojEhK4HPgiqalq4fGo+Gen691FEoqfrgcdZw9GTLLr/Vb5/42y+caV6wEXk/Ol64HFSVdOz/q0LWInIYFOAx1hVoJnRw7KYOXa416WISJJTgMeQc47KmmYWTh9NWtpAl5AREYkdBXgM1TS1cfBYB4t0+ryIxIECPIaqwrdPWzxDJ/CIyOBTgMdQVaCFCSOHMCl/qNeliEgKUB94jARDjjdrW1g8YzRmWv8WkfOnPvA4ea++lbKfvcFP/2Q+f/SpAW9YJCISFfWBD7LK8Pq3LmAlIvGiAI+RqpoWigtzKRyR43UpIpIiFOAx0NkdYsOewzr7UkTiSgEeA5v3HeFkV5CFWj4RkThSgMdAZU0LaQZXTFOAi0j8qI0wBqoCzVw8IY+8IZlelyIiSURthIPseHsX83/0Mnd/dhrfvX6W1+WISBJSG+Egebv2MMGQ49MzCrwuRURSjAL8PL0RaCYnM41LJ4/0uhQRSTEK8PP0RqCZ0qmjyc5I97oUEUkxCvDzcKC1ncChE1yp/m8R8YAC/DxUBk5dPlYBLiLxpwA/D2+Eb582a5xunyYi8ac+8Cg553gj0MziGWN0+zQRGRTqAx8k1QeOc/1P1/HAzZdw62UTvS5HRJKY+sBj7I1T69/FWv8WEW8owKNUGWhm2phhTBg5xOtSRCRFKcCj0Nkd4q3aFnWfiIinIgpwM1tmZtVmFjCz+84wZomZbTGz7Wb2h9iWmVg27zvCh51BBbiIeCrjbAPMLB14GLgWqAc2mNlK59yOXmNGAj8Hljnn9plZ4SDVmxDW7m4iI81YPEOXjxUR70RyBF4KBJxztc65TuBpYHm/MV8CnnPO7QNwzh2KbZmJ5bVdhyiZMorhObp8rIh4J5IAnwDU9XpeH36ttwuBUWa21sw2mdmdA32jpqYmSkpKTj8qKiqiq9pDB1rb2XXgOEtmJvUvGSLisYqKitNZCQy4XnvWJRRgoLNU+jePZwALgKuBIcCbZvaWc25370EFBQX4vQ98bXXPLxdLFeAiMojKy8tPncSDmTUPNCaSAK8Hep+pUgQ0DDCm2TnXBrSZ2TpgHrCbJLO2uokL8nK4cGyu16WISIqLZAllA1BsZlPNLAu4DVjZb8zzwJVmlmFmQ4HLgZ2xLdV7nd0h3gg0s2RmIWY6fV5EvHXWI3DnXLeZ3QusAdKBx51z283s7vD7K5xzO83sReBdIAQ85pzbNpiFe2HTB0c40dHNkpm6+46IeC+SJRScc6uB1f1eW9Hv+YPAg7ErLfGsrT5EZrqp/1tEEoLOxDwHa6ubuGxKPrnZEf27JyIyqBTgEWo4epLqg8fVfSIiCUPXA4/QK7vC7YOztP4tIvGh64HHyB2/fJv9R07yync+qw4UEYkrXQ/8PLSe7OLNmhauvWiswltEEoYCPAJrqw/RHXJcN2es16WIiJymAI/ASzsOMiY3m/kTR3ldiojIaQrws+joDrJ21yGunVNIum5eLCIJRAF+Fm/WtNDWGeS6OeO8LkVEpA8F+Fms2X6QYVnpLJyumzeISGJRH/gn6A6GeGn7AZbMKiQnM93rckQkxagP/Dy8/n4Td/xyPSu+cinL5l7gdTkikqLUBx6FF7Y2MiwrXXffEZGEpAA/g87uEL/f1sh1F43T8omIJCQF+Bm8EWjiWHs3n79ESycikpgU4GfwwtZGRuRkcGWxLl4lIolJAT6A9q4gL+04yLK548jK0BSJSGJSOg3gv3Ye5ERHN2XzxntdiojIGakPfADPbqzngrwcFk3XrdNExDvqAz9HB1rbWXT/K3xzyQz+5/UzvS5HRER94JH6zTv1hBzcsqDI61JERD6RArwX5xz/same0in5TBkzzOtyREQ+kQK8l00fHGFPcxu3lOjoW0QSnwK8l3/fUMfQrHRuvFgn74hI4lOAhx1p62Tl1gaWz5/AsOwMr8sRETkrBXjY0xvq6OgO8aeLJntdiohIRNQHTs91v3/91gdcMS2fWeNGeF2OiAigPvCIvLjtAHf/epOu+y0iCUl94J/gX6r2Mj4vh2tmj/W6FBGRiKV8gG/b38qbtS3csXAKGekpPx0i4iMRJZaZLTOzajMLmNl9A7y/xMxazWxL+PGD2Jc6OB5ZW8Pw7Ay+fMUkr0sRETknZ+2XM7N04GHgWqAe2GBmK51zO/oNfd059/lBqHHQ1DSdYPW2Rr65ZDojcjK9LkdE5JxEcgReCgScc7XOuU7gaWD54JYVHyvW1pCdkcbXFk/1uhQRkXMWSYBPAOp6Pa8Pv9bfQjPbama/N7OLYlLdINrT3MZzm/dze+kkxuRme12OiMg5i+SUQxvgtf69h+8Ak51zJ8zsc8BvgeL+GzU1NVFS8lEnTHl5+akex7h76KVqsjPS+OaSGZ7sX0Tkk1RUVFBRUXHq6YA3J4gkwOuBib2eFwENvQc45471+nq1mf3czMY455p7jysoKCAR+sC37W/lhXcb+fOrZlAwXEffIpJ4eh/gmlnzQGMiWULZABSb2VQzywJuA1b2HmBm48zMwl+Xhr9vy3nUPmicc/zk9zsZOTSTb3xmmtfliIhE7axH4M65bjO7F1gDpAOPO+e2m9nd4fdXALcA/83MuoGTwG0unqd4noMXtx2gMtDCD2+6SJ0nIuJrKXUq/cnOIFc/tJYRQzJ54b9/WifuiIgv6FR64J9ffZ+G1nZ+tHyuwltEfC9lUmxr3VF+8YcavrigiNKp+V6XIyJy3lIiwNu7gnzn2a2MHZHD35bN8bocEZGYSInrgf9k9U4Ch05w/82X6INLEfGNlL8e+MqtDfz5U5v5xqen8v3P6+hbRPwnJT/E3Nl4jPt+8y4LJo/ir26Y5XU5IiIxlbQBvv/oSb76xHpG5GTy8JcuJVNdJyKSZJIy1ZpPdPDVx9fzYUeQJ79+GePycrwuSUQk5iK5FoqvHGht58uPvcX+oyd54qulukmxiCStpArwHQ3HKP9/GznS1smvvn65+r1FJKklxRKKc47/2FTPFx6ppDvoeKr8CoW3iCS9uB6Bn+oDLysro6ysLCbfc29zGz9YuZ11u5sonZrPw1+6VJeIFZGksGrVqlPnzSRPH7hzjs11R/n1Wx/w2837yclM5y+vn8kdC6eQnjbQ/SdERPzrTH3gvlgDX7e7iQPH2jnc1smuxmNs2HuE/UdPkpOZxtcXT6X8s9MoHK5OExFJLb4I8L//3Q52HzwBQOHwbOZPHMm3rylm2dxxDNep8SKSonwR4I98ZQFZ6WmMHJqpwBYRCfNFgE8vyPW6BBGRhJMUbYQiIqlIAS4i4lMpcT1wERE/SvnrgYuI+J3vrwdeUVHhdQkJRfPRl+ajL81HX8k6H3EP8GiXT6L9A4hmf9HWGM99RTMf8fy5ot0uWecj3nOv+egrWefDNwEez/354Q8uGskc4NHww8+l+YjNdvHcVzznI65r4GbWBBwFWqPYfAzQHMV2eVHsL5pt4r2vaOYjnj9XtNsl63zEe+41H335fT4mO+cK+r8Y1wAXEZHY8c2HmCIi0pcCXETEpxI+wM1smZlVm1nAzO7zuh6vmdnjZnbIzLZ5XYvXzGyimb1mZjvNbLuZfcvrmrxkZjlmtt7Mtobn44de15QIzCzdzDab2Qte1xJrCR3gZpYOPAzcAMwBbjezOd5W5bkngWVeF5EguoHvOOdmA1cA96T4/x8dwFXOuXnAfGCZmV3hbUkJ4VvATq+LGAwJHeBAKRBwztU65zqBp4HlHtfkKefcOuCw13UkAudco3PunfDXx+n5SzrB26q843qcCD/NDD9SukvBzIqAG4HHvK5lMCR6gE8A6no9ryeF/4LKmZnZFOBTwNsel+Kp8HLBFuAQ8LJzLqXnA/gp8JdAyOM6BkWiB/hAN7hM6SMK+TgzywV+A3zbOXfM63q85JwLOufmA0VAqZnN9bgkz5jZ54FDzrlNXtcyWBI9wOuBib2eFwENHtUiCcjMMukJ7391zj3ndT2Jwjl3FFhLan9eshi4ycz20rP8epWZ/drbkmIr0QN8A1BsZlPNLAu4DVjpcU2SIMzMgF8CO51z/+B1PV4zswIzGxn+eghwDbDL06I85Jz7a+dckXNuCj3Z8apz7iselxVTCR3gzrlu4F5gDT0fUD3jnNvubVXeMrOngDeBmWZWb2Z/5nVNHloM3EHPkdWW8ONzXhfloQuA18zsXXoOfl52ziVd65x8RKfSi4j4VEIfgYuIyJkpwEVEfEoBLiLiUwpwERGfUoCLiPiUAlxExKcU4CIiPqUAFxHxqf8PMCzA2SKJ1LwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(A[:], ms.U[:] )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1.0)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD1CAYAAACvOmWkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATlklEQVR4nO3dbYylZ13H8e//PMzMzuyy3bptXSqkilhUIsWOSnwgRURrk00hUSMabAxhNApi9AWEF0LiG0zEp6iQE2ioLwSJiOw2KDZVrIQHoVja4oogklrZ2EFw6c7szsw55++L+z5zzmy3nTPnYebs3e8nOTn3w3Wuc+2V3d/ce5//XCcyE0lSNdQOegCSpMkx1CWpQgx1SaoQQ12SKsRQl6QKMdQlqUIa0+z88OHD+bznPW+sPs6dO8fRo0ev+D5WV1e55pprDnQMs9KHc9HnXPQ5F33333//+cw8MtKLM3Nqj8XFxRzXa17zmkr0cfPNNx/4GGalD+eiz7nocy76gLUcMXenevul0+mwsrLC6dOnR+7j5MmTY49jVvqYhTHMSh+zMIZZ6WMWxjArfczCGA6yj9OnT7OysgJQH/W9I6f4G6VLS0u5trY2tf6vJMvLy3z6058+6GHMBOeiz7nocy76ImI9M5dGee1Ur9SPHz8+ze6vKOVPX+FcDHIu+pyLHVZHfeGuV+oRsQDcB8xTfLD6l5n55oi4GvgL4Abgy8DPZObXB1+7vLyc/uSVpL2JiPszc3mU1w5zpb4B/GhmvgC4Cbg1Il4EvBG4NzOfC9xb7kuSDtCuoV5+GHu+3G2WjwRuB+4qj98FvHwaA5QkDW+oe+oRUY+IB4DHgHsy85PAdZl5FqB8vvbS162urrK8vLz9aLVaExy6JFVHq9Xazkpg5A8k91T9EhFXAR8AXgd8NDOvGjj39cw8Ntjee+qStHfTvqe+LTP/D/gIcCvwPxFxohzACYqreEnSAdo11CPimvIKnYg4BPwY8G/AKeCOstkdwAenNEZJ0pCGWfvlBHBXRNQpfgi8LzPvjoiPA++LiFcDjwA/PcVxSpKGsGuoZ+aDwAsvc/x/gZdOY1CSpNG49K4kVchUQ/3cuXNjL+glSU8XAwt6jbxu71QX9LKkUZL2bt9KGiVJs81Ql6QKMdQlqUIMdUmqEENdkirEkkZJmhGWNEpSBVnSKEkCDHVJqhRDXZIqxFCXpAqx+kWSZoTVL5JUQVa/SJIAQ12SKsVQl6QKMdQlqUIMdUmqEEsaJWlGWNIoSRVkSaMkCRgi1CPiWRHxDxFxJiI+FxGvL4+/JSL+OyIeKB+3TX+4kqSn0hiiTRv4zcz8TEQcAe6PiHvKc7+fmb87veFJkvZi11DPzLPA2XL78Yg4A1w/7YFJkvZuT/fUI+IG4IXAJ8tDr42IByPizog4NunBSZL2ZuhQj4jDwPuBX8/MbwBvB54D3ERxJf+2S1+zurrK8vLy9qPVak1m1JJUMa1WazsrgeOj9jNUSWNENIG7gQ9n5u9d5vwNwN2Z+fzB45Y0StLeTbWkMSICeBdwZjDQI+LEQLNXAA+PMgBJ0uQMU/3yQ8CrgIci4oHy2JuAV0bETUACXwZ+aQrjkyTtwTDVLx8F4jKnPjT54UiSxuFvlEpShbiglyTNCBf0kqQKckEvSRJgqEtSpRjqklQhhrokVYihLkkVYkmjJM0ISxolqYIsaZQkAYa6JFWKoS5JFWKoS1KFWP0iSTPC6hdJqiCrXyRJgKEuSZViqEtShRjqklQhhrokVYgljZI0IyxplKQKsqRRkgQMEeoR8ayI+IeIOBMRn4uI15fHr46IeyLiC+XzsekPV5L0VIa5Um8Dv5mZ3wm8CPjViPgu4I3AvZn5XODecl+SdIB2DfXMPJuZnym3HwfOANcDtwN3lc3uAl4+pTFKkoa0p3vqEXED8ELgk8B1mXkWiuAHrp346CRJe9IYtmFEHAbeD/x6Zn4jInZ9zerqKsvL/Q9wV1ZWeuU6kqQBrVaLVqvV2z0+aj9DlTRGRBO4G/hwZv5eeezzwC2ZeTYiTgAfycwbB19nSaMk7d1USxqjuCR/F3CmF+ilU8Ad5fYdwAdHGYAkaXKGuf3yQ8CrgIci4oHy2JuAtwLvi4hXA48APz2VEUqShrZrqGfmR4Enu4H+0skOR5I0Dn+jVJIqxAW9JGlGuKCXJFWQC3pJkgBDXZIqxVCXpAox1CWpQgx1SaoQSxolaUZY0ihJFWRJoyQJMNQlqVIMdUmqEENdkirEUJekCrGkUZJmhCWNklRBljRKkgBDXZIqxVCXpAox1CWpQqx+kaQZYfWLJFWQ1S+SJGCIUI+IOyPisYh4eODYWyLivyPigfJx23SHKUkaxjBX6u8Gbr3M8d/PzJvKx4cmOyxJ0ih2DfXMvA/42j6MRZI0pnHuqb82Ih4sb88cm9iIJEkjGzXU3w48B7gJOAu87XKNVldXWV5e3n60Wq0R306Sqq3Vam1nJXB81H6GKmmMiBuAuzPz+Xs5Z0mjJO3dvpc0RsSJgd1XAA8/WVtJ0v5p7NYgIt4D3AIcj4hHgTcDt0TETUACXwZ+aXpDlCQNa9dQz8xXXubwu6YwFknSmPyNUkmqEENdkirEVRolaUa4SqMkVZCrNEqSAENdkirFUJekCjHUJalCrH6RpBlh9YskVZDVL5IkwFCXpEox1CWpQgx1SaoQQ12SKsSSRkmaEZY0SlIFWdIoSQIMdUmqFENdkirEUJekCrH6RZJmhNUvklRBVr9IkoAhQj0i7oyIxyLi4YFjV0fEPRHxhfL52HSHKUkaxjBX6u8Gbr3k2BuBezPzucC95b4k6YDtGuqZeR/wtUsO3w7cVW7fBbx8ssOSJI1i1Hvq12XmWYDy+drJDUmSNKrGNDtfXV1lebn/Ae7KykqvXEeSNKDVatFqtXq7x0ftZ6iSxoi4Abg7M59f7n8euCUzz0bECeAjmXnjpa+zpFGS9u4gShpPAXeU23cAHxyxH0nSBA1T0vge4OPAjRHxaES8Gngr8LKI+ALwsnJfknTAdr2nnpmvfJJTL53wWCRJY/I3SiWpQgx1SaoQV2mUpBnhKo2SVEGu0ihJAgx1SaoUQ12SKsRQl6QKsfpFkmaE1S+SVEFWv0iSAENdkirFUJekCjHUJalCDHVJqhBLGiVpRljSKEkVZEmjJAkw1CWpUgx1SaoQQ12SKsRQl6QKsaRRkmaEJY2SVEGWNEqSAGiM8+KI+DLwONAB2qP+ZJEkTcZYoV56SWZ+dQL9SJLG5O0XSaqQcUM9gb+LiPsjYuXSk6urqywvL28/Wq3WmG8nSdXUarW2sxI4Pmo/Y1W/RMQzM/MrEXEtcA/wusy8r3fe6hdJ2rsDq37JzK+Uz48BHwC+f5z+JEnjGTnUI2IpIo70toEfBx6e1MAkSXs3TvXLdcAHIqLXz59n5t9OZFSSpJGMHOqZ+SXgBRMciyRpTJY0SlKFuKCXJM0IF/SSpApyQS9JEmCoS1KlGOqSVCGGuiRViKEuSRViSaMkzQhLGiWpgixplCQBhrokVYqhLkkVYqhLUoWMs576rnrVLydPnuTkyZPTfCtJuiJstruc32hz/mKbxze2OH+xzdpmm8cvtvnYpz7DZx46A2NUv0w11I8ePeqXTUu64nW7ydpmezuMz2+0B4K5eF4rjz1+mTaD25ud7lO8Uw0Wvxvg3KhjnWqoS9JByUwubHU4v9FmfaPzhJB9fKMM4t7+E4J5i/MbbdbK1w5jvlHjyEKDw/MNluaL52detcDh+QaHF4pjR8rjhxeaxXN57vB8gyO9Nr8z+p/bUJc0E7rdIoR7wbq20WFts71jf32zt91mbbNo2zu/vtnpnytfO8yv4dSCJ4Tr0UNNvuWqQyzN1zk83+TwQhnGg8Fctu09luYbzDUO/mNKQ13SSHq3JHYN340yfDfbnN/osN5rtzlwbqPN+lZnqBAGWGjWtoN0ca7B4fk6Vy/N8axjiyzN17evknvneu2esfDEYD7UrFN+13IlGOrS00BmstHusl5e3a5vFkF8YbMfuIMBe/5ywby589z6Zmfo91+cGwzaYvvaIwssflN9O5yXyuOD7bbPzffPLc01qNeqE8KTZqhLM2az3WV9s7i9cGHgSnh9o8P6VnGlu7ZZPO/YL4N2vde+d6zc7w55FRwBS3NlkM71Q/XE0QUW53qhWy+vggdCd24gkOf7gbzYrFMzhPeNJY3SiNqdbhmqO696e1fB65s7A3dto1O0Kc+tbbS37yEP7m91hl+Paa5eY7EM1ENzdZbmirA9cXSBQ3ON7f2lMoQX5+rbV8qLc/1jvWA+PN9goWEIH5TTp0/3FkB0QS/pUludLhe2iiC9sNlhfbPDha02FzaLK+HeueJ42W6r2L+41b/yvbi1s01ve7P9VKVpO9Vr0Q/UueJKdnEgdHcG7c4gXpqvc6i5c3+xWYT4LHwwp8kbZ0Evb79o33W7xf3di1sdLrY7XNwqti9sFQG6sdU/d2GzF8zt7cB9YvjuDOLe+b1c8UJx22GxWefQXINDczUWmw0W5uosNut88zOa29uH5orHYvPS4O0H86VXwvONWqU+jNPsGivUI+JW4A+BOvDOzHzrREalfZGZbHWSzU6Xja0Om50um+3isdHusjEQuBe3utuhe3Gr0w/lrYE2A8c2trplYA/2UbTZyxXuoHotWGzWi3CdK0L0ULl9bLHJobnGjtA91CzbDWwvNIugHXxtb9vgVRWMHOoRUQf+BHgZ8CjwqYg4lZn/OqnBXekyk3Y3aXeSrW6XrXaXdjfZ6nSLY50uW52k3e32ty/TdrPdLYO3e0nwdvrnBsK4d7443unvD/Sz0Rk9XActNGssNOssNOr97WaxffXS3BOOzzdr5bHB4/1j85f017siPjRXp1kPQ1faxThX6t8PfDEzvwQQEe8Fbge2Q/3chS1OffYr9O7bZ0KSxXNCUgRfAiR0y+0d7coXbh8vt7vl9mC/nW7RR7uTdDLpdvvP7W7S6Sbd7D/vbAedbpdOt99PpzvwGOiv080ygIvg3Wx3aXeLoN4sA7sI6ul9XlELmG8U91TnGjXm6jXmy+3e80KzxjMWGjvbDZyfrw8eq2/3M7f9+joLjZ1B3Qvc+WbNK1tpBo0T6tcD/zWw/yjwA4MNHvnaOr/2nn8Z4y3GU4viv+z1WlCPoDaw3TteG9je2Q7qtRr1so9aBHON2na7Rq1Gsx406zUa9aBZK5/rxfFGvVZs13rbT9G2VqPZ6Ldt1IO58rlRe2JYF+f8gEyqklarNbhW1vFR+xkn1C93ibbj0vQ7rj3CB37jxUAQUbwgIspniPI4QK32xONRvktc8vpaeYxgx/EikKFRqxVtvIqUdIVYWVnpfT8pEfHVUfsZJ9QfBZ41sP8twFcGG8w3a3z7tUfGeAtJ0l6M83/4TwHPjYhvjYg54GeBU5MZliRpFCNfqWdmOyJeC3yYoqTxzsz83MRGJknas7Hq1DPzQ8CHJjQWSdKYLKGQpAqZaqg/8sgjrKys9BaoeVrza/36nIs+56LPuSgW9CorYJ49cieZObXH4uJijuvUqVOV6OPmm28+8DHMSh/ORZ9z0edc9AFrOWLuzvztl0lc5c9KH7MwhlnpYxbGMCt9zMIYZqWPWRjDrPQxqqkuvRsRjwOfH7Obo4zxzdoz1MdxYORfKJjQGGalD+eiz7nocy76bszMkX7JZ6qhLknaXzN/+0WSNDxDXZIqZCKhHhG3RsTnI+KLEfHGy5yPiPij8vyDEfG9k3jfWTTEXPx8OQcPRsTHIuIFBzHO/bDbXAy0+76I6ETET+3n+PbTMHMREbdExAMR8bmI+Mf9HuN+GeLfyNGIOB0Rny3n4hcPYpzTFhF3RsRjEfHwk5wfLTdHLZvpPSiWCPgP4NuAOeCzwHdd0uY24G8oFlR8EfDJcd93Fh9DzsUPAsfK7Z98Os/FQLu/p/jN5J866HEf4N+Lqyi+i+DZ5f61Bz3uA5yLNwG/U25fA3wNmDvosU9hLl4MfC/w8JOcHyk3J3Glvv1lGZm5CfS+LGPQ7cCfZeETwFURcWIC7z1rdp2LzPxYZn693P0ExeqWVTTM3wuA1wHvBx7bz8Hts2Hm4ueAv8rMRwAys6rzMcxcJHAkirWzD1OEent/hzl9mXkfxZ/tyYyUm5MI9ct9Wcb1I7Spgr3+OV9N8ZO4inadi4i4HngF8I59HNdBGObvxXcAxyLiIxFxf0T8wr6Nbn8NMxd/DHwnxVLeDwGvz8zxv3vxyjNSbo61oFdp1y/LGLJNFQz954yIl1CE+g9PdUQHZ5i5+APgDZnZqfgXmgwzFw3gZuClwCHg4xHxicz892kPbp8NMxc/ATwA/CjwHOCeiPinzPzGlMc2a0bKzUmE+q5fljFkmyoY6s8ZEd8DvBP4ycz8330a234bZi6WgfeWgX4cuC0i2pn51/sywv0z7L+Rr2bmGrAWEfcBLwCqFurDzMUvAm/N4sbyFyPiP4HnAf+8P0OcGSPl5iRuvwzzZRmngF8oP819EXAuM89O4L1nza5zERHPBv4KeFUFr8IG7ToXmfmtmXlDZt4A/CXwKxUMdBju38gHgR+JiEZELFJ83++ZfR7nfhhmLh6h+B8LEXEdcCPwpX0d5WwYKTfHvlLPJ/myjIj45fL8OygqG24DvgisU/wkrpwh5+K3gG8C/rS8Qm1n5vJBjXlahpyLp4Vh5iIzz0TE3wIPAl3gnZl52VK3K9mQfy9+G3h3RDxEcQviDZk57vIBMyci3gPcAhyPiEeBNwNNGC83XSZAkirE3yiVpAox1CWpQgx1SaoQQ12SKsRQl6QKMdQlqUIMdUmqEENdkirk/wHwSvd0EWwQogAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(A, ms.R**2 * ms.m * ms.Abar**2 * np.exp(2 * (ms.alpha-1) * ms.xi))\n",
    "\n",
    "plt.xlim(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb97f4e24f0>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD1CAYAAABJE67gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAchUlEQVR4nO3de3zU9Z3v8dcnkxsQCAQSEMKdCCIKhRQFesE7aindarvaVnvZNg9P69n2tKe7dtvt47S7PXX1dLe7j1pptNZ2u63bbi+CpaKrUjWoXASUW3ACaEICuQABIrnNfM8fGXBIB5iEmfnNb+b9fDzmwcxvvjPz+U7rmx/f+fx+P3POISIi/pPjdQEiIjI4CnAREZ9SgIuI+JQCXETEpxTgIiI+pQAXEfGp3FR+2JgxY9yUKVMAaG9vp7i4OO7XtrS0UFpaGvf4gb7/QMen4jM0Z+/rScVnpNucU/Edac4De//Nmze3Ouf+/MXOuZTdFixY4E753Oc+5wYi+rXxGOj7D3R8Kj5Dc078+2vOiR8/mNdozgN7f2CTi5Gpni2hLF++PK3efzD1pOIzkvn+mnP6fEYy3z8dv6Nkv3+2zNlcCo/ErKysdJs2bRrsaxnsa/1Kc84OmnN2uJA5m9lm51xl/+2++RGzqqrK6xJSTnPODppzdkjGnFO6B15RUeGuuuoqli9fnvR/foiI+N3q1atZvXo1Dz30UNA5V9H/ed8soYiIZCvfL6GIiMiZUtoHLiJniv4XsM7snNnMwMwS+p4KcEm5UNhx5O1uDnd0c7yzhxNdIU509tLR1cuJyK2jq5eTPSF6QmG6esN0R249oTDdoXced4ccoXCYsIOwc4TDjrDr+wznHCHX99g5RyjyXN+YvueiQ/OM/Dxj+9lDNvrhGWF8xvYL+bYkU+z+h2UU5gUS+p4KcEmo3lCYhiMnI7e3T//Z1N5JW0dfaB95u/u8oZafm8PQ/AB5gRzyAzkU5OaQn5vT9zi3b9vQobnkBXLIzTECOYYZ5Ng79wNm5JiRk2PkxHoucj+a8c6G6Oeihw3mNdFP9N8HO/M1id1Dk/SRm5P4/20V4DJoHV29bK0/yvYD7dQePM7ug8cJNp+gOxQ+PSaQY1xUXMhFxYVUlBVRMiyf0UUFjCnKZ9TQfEYMyaOoIEBRQR7DCgIUFeQyrKAvmEXk3BTgErfOnhAv723j+T2tbHrzMDsajxEK9+1KjxtRyMxxw3lvxRimlxUxqWQo5aOGMG5EIbkKY5GkSGmAt7e3U1VVpT5wHznZHeKpnQd54rUmXnyjlZM9IQpyc3jXpJF8ful0KqeUMLe8mJFD870uVSTjnOoDB2KeNUt94BLT1vqj/MfLb7Lm9SY6ukOMLy7k2tljuXpWGVdOG53wH2NE5OzO1geuJRQ5LRx2PFfbzI+e38uGfYcZlh/gpssu4pYF5SycUkJOEn6EEZHBU4ALADXBVv7vml3saDzG+OJCvnHzJdy2cBJFBfq/iEi60n+dWa7+8Nt88/HtPFfbwoSRQ/jeR+bywXnj1QUi4gMK8CwVCjseXb+f/7e2lhyDv7tpFncumqK1bREfUYBnoZbjXXzxsS2sr2vjqpml/ONfXMaEkUO8LktEBkgBnmVe2dvG3b/cwvHOHu675XI+Ulme8PMziEhqpHSh81QfeKSvUVLst6828PGHX2F4YS6//8ISPvruiQpvkTS2evXqUxeCUB94tnLO8cN1ddy/tpbF00ez8o4FjCjM87osEYmT+sCzlHOO+9bW8uC6OlbMG8/9t84lP1cdJiKZQAGe4f7tmSAPrqvj9oWT+M6H5uhgHJEMol2xDPbQ83v5l//ew60LyhXeIhkorgA3s2VmVmtmQTO7J8bzo8zsd2b2mpltMLM5iS9VBuLJ7U18Z80ubr7sIv7plssV3iIZ6LwBbmYB4AHgRmA2cLuZze437O+Arc65y4E7gX9NdKESv+0H2vlf/7mNeRNH8r2PziWg8BbJSPHsgS8Egs65vc65buAxYEW/MbOBZwCcc7uBKWY2NqGVSlxaT3TxuZ9tYtTQPKrvXKAjK0UyWDwBPgGoj3rcENkWbRvwYQAzWwhMBsr7v1FLSwuVlZWnb9XV1YOrWmJyzvHVX2+jraOb6jsrKRte6HVJIjJI1dXVp7MSGBNrTDxdKLH+/d2/efxe4F/NbCvwOrAF6O3/otLSUtQHnjyPrt/Pc7UtfOuDlzJnQsy+fxHxiaqqqlMH8WBmrbHGxBPgDcDEqMflQGP0AOfcMeDTkQ8yYF/kJimys/EY312zm2tmlXHnoslelyMiKRDPEspGoMLMpppZPnAbsCp6gJmNjDwH8Fng+UioSwr0hMJ85dfbKB6ax323Xq7D40WyxHn3wJ1zvWZ2N7AWCACPOOd2mNldkedXApcAPzOzELAT+Ksk1iz9/KRmH7uajrHyE/MZXVTgdTkikiJxHYnpnFsDrOm3bWXU/ZeAisSWJvGoP/w2//L0G1x7SRk3XDrO63JEJIV0JKaPOef4+8e3YwbfWjFHSyciWUYB7mPP1TazrraFL193sS7IIJKFdD5wn+oNhfnumt1MGT2UOxdN8bocEUkCnQ88Qz224S3u+e3rPPjx+dx42UVelyMiSXS284FrCcWH3u7u5Z+f3sP8SSNZNkc/XIpkKwW4D/2kZj/Nx7v4+s2X6IdLkSymAPeZjq5eHn5hL1fNLGXB5BKvyxERDynAfeYXr7zFkbd7uPtqtd2LZDsFuI909oSofmEvS2aMZsHkUV6XIyIeU4D7yK821dNyvIu7r9Let4ioD9w3ekJhfvSnvVROHsWV07T2LZIN1AeeIf7wWhNf+MWrPHRnJdfN1sWORLKJ+sB97tH1+5hUMpSrZ5V5XYqIpAkFuA9sP9DOxv1HuHPRZF2gWEROU4D7wE9q9jM0P8BHKieef7CIZA0FeJprPdHF6m2N3DK/nOIheV6XIyJpRAGe5n61qZ7uUJhPLtZ1LkXkTArwNOac49ebGlg4tYQZZcO9LkdE0oz6wNPYhn2H2dfawUe19i2SldQH7mNf/tVWntpxiA1fv4ah+XFdvlREMpD6wH3meGcPa15vYvnc8QpvEYlJAZ6mVm9rorMnzEcry70uRUTSlAI8Tf1qUz0Xjy1i3sSRXpciImlKAZ6G3mzrYGv9UT48v1xX3BGRs1KAp6FVWxsBWD53vMeViEg6U4CnGeccq7Y1snBKCRNGDvG6HBFJY+oDTzO7Dx7njeYTLJ+nvW+RbKc+cJ+594+7efiFvWz4+rWUDMv3uhwRSQPqA/cB5xyrtzXynooxCm8ROS8FeBp59a2jHDh6kg/qx0sRiUNcAW5my8ys1syCZnZPjOeLzWy1mW0zsx1m9unEl5r5ntzeRH4gR5dME5G4nDfAzSwAPADcCMwGbjez2f2GfQHY6ZybCywFvmdmWgMYAOccT+44yJIZoxleqPN+i8j5xbMHvhAIOuf2Oue6gceAFf3GOGC49R11UgQcBnoTWmmG29l0jPrDJ1k2Z5zXpYiIT8QT4BOA+qjHDZFt0X4AXAI0Aq8DX3TOhRNSYZZYu/0gOQbXXqLlExGJTzwBHutY7v69hzcAW4HxwDzgB2Y2ov+LWlpaqKysPH2rrq4eYLmZ68kdB1k4tYTRRQVelyIiaaC6uvp0VgJjYo2J5zylDUD0FQXK6dvTjvZp4F7X11QeNLN9wCxgQ/Sg0tJS1Af+5+paTrDn0An+z/L+Py2ISLaqqqo6dRAPZtYaa0w8e+AbgQozmxr5YfI2YFW/MW8B10Q+aCwwE9g7yLqzztodBwG4/lKtf4tI/M67B+6c6zWzu4G1QAB4xDm3w8zuijy/EvgH4FEze52+JZe/dc7F/BtD/tzTOw9xeXkx43XuExEZgLgu9eKcWwOs6bdtZdT9RuD6xJaWHdpOdLG1/ihfuuZir0sREZ/RkZgeW1fbgnNw9awyr0sREZ9RgHvs2d3NlA0v4NLxf9a0IyJyTgpwD/WEwjy/p4WrZpaRk6Mr74jIwOh84B7atP8Ix7t6uUrLJyISg84Hnsa+84ed/HT9m7z6zesoKojr92QRyUI6H3gaenZ3M1dMK1F4i8igKMA98mZbB3UtHVw1U8snIjI4CnCPPLu7GYBrLlGAi8jgKMA98lxtC9NKhzF59DCvSxERn1KAe6CrN8SGfW28r6LU61JExMcU4B7Y8tZROnvCLJ4+2utSRMTH1AfugfXBVnIMrpimABeRs1MfeBq69cH19IQdj39hideliIgPqA88TXR09bK1/ihLtHwiIhdIAZ5iG/YdpjfsWDw95hWSRETipgBPsZpgK/m5OVROGeV1KSLicwrwFKupa2PBpFEU5gW8LkVEfE4BnkKHO7rZ1XSMJTO0/i0iF04BnkIv1bUBsEjr3yKSAOoDT6GaulaKCnKZWx6zpVNE5AzqA08jS+9/jumlRfz4U+/2uhQR8RH1gXvswNGT7G97m8UztHwiIomhAE+R9cFWAP2AKSIJowBPkfV1bYwels/FZcO9LkVEMoQCPAWcc9QEW1k0fbSuPi8iCaMAT4G6lg6aj3exROvfIpJACvAUWF8XWf9W/7eIJJD6wFOgJtjKhJFDmFgyxOtSRMRH1AfusVDY8a5vP8WyOeO479a5XpcjIj6kPnCP7Gw8xrHOXq1/i0jCKcCTrCay/r1IF3AQkQSLK8DNbJmZ1ZpZ0MzuifH8V81sa+S23cxCZlaS+HL9pybYSkVZEWXDC70uRUQyzHkD3MwCwAPAjcBs4HYzmx09xjl3v3NunnNuHvA14E/OucNJqNdXunpDbNx/WMsnIpIU8eyBLwSCzrm9zrlu4DFgxTnG3w78MhHF+d3Wt47S2RNmsZZPRCQJ4gnwCUB91OOGyLY/Y2ZDgWXAby68NP+rqWsjx+CKaQpwEUm83DjGxDr2+2y9h8uBmrMtn7S0tFBZ+U4nTFVV1akex4y0PtjKZROKKR6S53UpIuIz1dXVVFdXn3oYcx02ngBvACZGPS4HGs8y9jbOsXxSWlpKtvSBd3T1srX+KJ973zSvSxERH4rewTWz1lhj4llC2QhUmNlUM8unL6RX9R9kZsXA+4HHB11xBtmw/zC9YafD50Ukac67B+6c6zWzu4G1QAB4xDm3w8zuijy/MjL0L4CnnHMdSavWR9YHW8nPzaFyyiivSxGRDBXPEgrOuTXAmn7bVvZ7/CjwaKIK87uaYBsLJo2iMC/gdSkikqF0JGYSHO7oZmfTMbUPikhSKcCT4OW9bQC6/qWIJJUCPAlqgq0UFeQytzzmGSBFRBJC5wNPgvV1bVwxtYTcgP5+FJHB0/nAU6zx6EkW3/ss37j5Ej77XvWAi8iF0/nAU2R9Xd/6t05gJSLJpgBPsPXBVkYPy2fm2OFelyIiGU4BnkDOOWrqWlk0fTQ5ObFOISMikjgK8ASqa+ng0LEuFuvweRFJAQV4Aq2PXD5tyQwdwCMiyacAT6D1wTYmjBzCpJKhXpciIllAfeAJEgo7XtrbxpIZozHT+reIXDj1gafI6w3tLP/Bi3z/L+fxoXfFvGCRiMigqA88yWoi6986gZWIpIoCPEHW17VRUVZE2YhCr0sRkSyhAE+A7t4wG/cd1tGXIpJSCvAE2PLWEU72hFik5RMRSSEFeALU1LWRY3DlNAW4iKSO2ggTYH2wlcsmFFM8JM/rUkQkg6iNMMmOd/Yw79tPc9f7p/HVG2Z5XY6IZCC1ESbJK3sPEwo73jOj1OtSRCTLKMAv0IvBVgrzcpg/eaTXpYhIllGAX6AXg60snDqagtyA16WISJZRgF+Ag+2dBJtP8F71f4uIBxTgF6AmeOr0sQpwEUk9BfgFeDFy+bRZ43T5NBFJPfWBD5JzjheDrSyZMUaXTxORpFAfeJLUHjzODd9/nvtuuZyPvnui1+WISAZTH3iCvXhq/btC698i4g0F+CDVBFuZNmYYE0YO8boUEclSCvBB6O4N8/LeNnWfiIin4gpwM1tmZrVmFjSze84yZqmZbTWzHWb2p8SWmV62vHWEt7tDCnAR8VTu+QaYWQB4ALgOaAA2mtkq59zOqDEjgR8Cy5xzb5lZWZLqTQvr9rSQm2MsmaHTx4qId+LZA18IBJ1ze51z3cBjwIp+Yz4G/NY59xaAc645sWWml+d2N1M5ZRTDC3X6WBHxTjwBPgGoj3rcENkW7WJglJmtM7PNZnZnrDdqaWmhsrLy9K26unpwVXvoYHsnuw8eZ+nMjP5Hhoh4rLq6+nRWAjHXa8+7hALEOkqlf/N4LrAAuAYYArxkZi875/ZEDyotLcXvfeDravv+cXGVAlxEkqiqqurUQTyYWWusMfEEeAMQfaRKOdAYY0yrc64D6DCz54G5wB4yzLraFi4qLuTisUVelyIiWS6eJZSNQIWZTTWzfOA2YFW/MY8D7zWzXDMbClwB7Epsqd7r7g3zYrCVpTPLMNPh8yLirfPugTvnes3sbmAtEAAecc7tMLO7Is+vdM7tMrMngdeAMPCwc257Mgv3wuY3j3Ciq5elM3X1HRHxXjxLKDjn1gBr+m1b2e/x/cD9iSst/ayrbSYvYOr/FpG0oCMxB2BdbQvvnlJCUUFcf++JiCSVAjxOjUdPUnvouLpPRCRt6HzgcXpmd6R9cJbWv0UkNXQ+8AS548evcODISZ75yvvVgSIiKaXzgV+A9pM9vFTXxnWXjlV4i0jaUIDHYV1tM71hx/Wzx3pdiojIaQrwODy18xBjigqYN3GU16WIiJymAD+Prt4Q63Y3c93sMgK6eLGIpBEF+Hm8VNdGR3eI62eP87oUEZEzKMDPY+2OQwzLD7Boui7eICLpRX3g59AbCvPUjoMsnVVGYV7A63JEJMuoD/wCvPBGC3f8eAMrPzGfZXMu8rocEclS6gMfhCe2NTEsP6Cr74hIWlKAn0V3b5g/bm/i+kvHaflERNKSAvwsXgy2cKyzlw9crqUTEUlPCvCzeGJbEyMKc3lvhU5eJSLpSQEeQ2dPiKd2HmLZnHHk5+orEpH0pHSK4b93HeJEVy/L5473uhQRkbNSH3gMv97UwEXFhSyerkuniYh31Ac+QAfbO1l87zN8fukM/vcNM70uR0REfeDx+s2rDYQd3Lqg3OtSRETOSQEexTnHf21uYOGUEqaMGeZ1OSIi56QAj7L5zSPsa+3g1krtfYtI+lOAR/nPjfUMzQ9w82U6eEdE0p8CPOJIRzertjWyYt4EhhXkel2OiMh5KcAjHttYT1dvmE8unux1KSIicVEfOH3n/f75y29y5bQSZo0b4XU5IiKA+sDj8uT2g9z1880677eIpCX1gZ/DT9fvZ3xxIddeMtbrUkRE4pb1Ab79QDsv7W3jjkVTyA1k/dchIj4SV2KZ2TIzqzWzoJndE+P5pWbWbmZbI7dvJr7U5HhwXR3DC3L5+JWTvC5FRGRAztsvZ2YB4AHgOqAB2Ghmq5xzO/sNfcE594Ek1Jg0dS0nWLO9ic8vnc6IwjyvyxERGZB49sAXAkHn3F7nXDfwGLAiuWWlxsp1dRTk5vDpJVO9LkVEZMDiCfAJQH3U44bItv4Wmdk2M/ujmV2akOqSaF9rB7/dcoDbF05iTFGB1+WIiAxYPIccWoxt/XsPXwUmO+dOmNlNwO+Biv4vamlpobLynU6YqqqqUz2OKfe9p2opyM3h80tnePL5IiLnUl1dTXV19amHMS9OEE+ANwATox6XA43RA5xzx6LurzGzH5rZGOdca/S40tJS0qEPfPuBdp54rYm/vnoGpcO19y0i6Sd6B9fMWmONiWcJZSNQYWZTzSwfuA1YFT3AzMaZmUXuL4y8b9sF1J40zjm++8ddjByax2ffN83rckREBu28e+DOuV4zuxtYCwSAR5xzO8zsrsjzK4Fbgf9hZr3ASeA2l8pDPAfgye0HqQm28a0PXqrOExHxtaw6lP5kd4hrvreOEUPyeOJ/vkcH7oiIL+hQeuDfnn2DxvZOvr1ijsJbRHwva1JsW/1RfvSnOj6yoJyFU0u8LkdE5IJlRYB39oT4yq+3MXZEIX+/fLbX5YiIJERWnA/8u2t2EWw+wb23XK4fLkXEN7L+fOCrtjXy17/cwmffM5VvfEB73yLiP1n5I+aupmPc85vXWDB5FH974yyvyxERSaiMDfADR0/yqZ9sYERhHg98bD556joRkQyTkanWeqKLTz2ygbe7Qjz6mXczrrjQ65JERBIunnOh+MrB9k4+/vDLHDh6kp98aqEuUiwiGSujAnxn4zGq/n0TRzq6+dlnrlC/t4hktIxYQnHO8V+bG/jwgzX0hhy/rLpS4S0iGS+le+Cn+sCXL1/O8uXLE/Ke+1s7+OaqHTy/p4WFU0t44GPzdYpYEckIq1evPnXcTOb0gTvn2FJ/lJ+//Ca/33KAwrwAf3PDTO5YNIVATqzrT4iI+NfZ+sB9sQb+/J4WDh7r5HBHN7ubjrFx/xEOHD1JYV4On1kylar3T6NsuDpNRCS7+CLA//EPO9lz6AQAZcMLmDdxJF+6toJlc8YxXIfGi0iW8kWAP/iJBeQHchg5NE+BLSIS4YsAn15a5HUJIiJpJyPaCEVEspECXETEp7LifOAiIn6U9ecDFxHxO9+fD7y6utrrElJOc84OmnN2SMacPQvwgS6jDHTyA33/wSzrJPszNOfEv7/mnPjxg3mN5pyY9/dNgCf7/dPxP+yB0pwTPz5Vn5HM90/H7yjZ758tc07pGriZtQBvRh4WA+0DePkYoHUA4wf6/gMdn4rP0Jy9rycVn5Fuc07Fd6Q5D+z9JzvnSvsPSGmAi4hI4vjmR0wRETmTAlxExKfSPsDNbJmZ1ZpZ0Mzu8bqeRDGzR8ys2cy2R20rMbOnzeyNyJ+jop77WuQ7qDWzG7yp+sKY2UQze87MdpnZDjP7YmR7xs7bzArNbIOZbYvM+VuR7Rk751PMLGBmW8zsicjjjJ6zme03s9fNbKuZbYpsS+6cnXNpewMCQB0wDcgHtgGzva4rQXN7HzAf2B617T7gnsj9e4B/ityfHZl7ATA18p0EvJ7DIOZ8ETA/cn84sCcyt4ydN2BAUeR+HvAKcGUmzzlq7l8GfgE8EXmc0XMG9gNj+m1L6pzTfQ98IRB0zu11znUDjwErPK4pIZxzzwOH+21eAfw0cv+nwIeitj/mnOtyzu0DgvR9N77inGtyzr0auX8c2AVMIIPn7fqciDzMi9wcGTxnADMrB24GHo7anNFzPoukzjndA3wCUB/1uCGyLVONdc41QV/YAWWR7Rn3PZjZFOBd9O2RZvS8I0sJW4Fm4GnnXMbPGfg+8DdAOGpbps/ZAU+Z2WYzq4psS+qc0/184LEucJmNfY8Z9T2YWRHwG+BLzrljZme9jmlGzNs5FwLmmdlI4HdmNuccw30/ZzP7ANDsnNtsZkvjeUmMbb6ac8QS51yjmZUBT5vZ7nOMTcic030PvAGYGPW4HGj0qJZUOGRmFwFE/myObM+Y78HM8ugL7/9wzv02sjnj5w3gnDsKrAOWkdlzXgJ80Mz207fsebWZ/ZzMnjPOucbIn83A7+hbEknqnNM9wDcCFWY21czygduAVR7XlEyrgE9G7n8SeDxq+21mVmBmU4EKYIMH9V0Q69vV/jGwyzn3z1FPZey8zaw0sueNmQ0BrgV2k8Fzds59zTlX7pybQt9/s8865z5BBs/ZzIaZ2fBT94Hrge0ke85e/3Ibxy+7N9HXrVAHfN3rehI4r18CTUAPfX8b/xUwGngGeCPyZ0nU+K9HvoNa4Eav6x/knN9D3z8TXwO2Rm43ZfK8gcuBLZE5bwe+GdmesXPuN/+lvNOFkrFzpq9TblvktuNUViV7zjqUXkTEp9J9CUVERM5CAS4i4lMKcBERn1KAi4j4lAJcRMSnFOAiIj6lABcR8SkFuIiIT/1/PjCo0c6KDMkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot( ms.U[:])\n",
    "#plt.xlim(0,50)\n",
    "#plt.ylim(0.87,0.92)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb97f4cefa0>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD1CAYAAABJE67gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbXUlEQVR4nO3de3SV9Z3v8fc3OwkhXGIwAZFwlYukUKREQLFVvBVtqb05o1ad8XSasUe7Ol2dnjozXdMZXXNOz/SyznS0Q6NjnU6dsdOlbaFFKbVKnYAKcg8I7kQukUB2AoRrIMn+nj8SMI2J2SR772dfPq+19iLP3s+T/c0j+fjjt7/P8zN3R0RE0k9O0AWIiMjAKMBFRNKUAlxEJE0pwEVE0pQCXEQkTSnARUTSVG4y36ykpMQnTZpEJBKhtLT0go9vaWmhqKhoQO890GOTfVyyz02yjxvMsZl+bpJ9XgbznulyXCZkTUlJCatWrVrl7kve86K7J+0xb948d3c/9+eF+sIXvjCg4wZzbLKPS/a5SfZxgzk2089Nss/LYN4zXY7LlKwBNngvmZpWUyhLly5N+rHJPm6g0unn07mJ73GDkS4/Y7r8nRnMsQM5zjyJV2JWVFT4hg0bqKioYMOGDUl733Sic9M3nZve6bz0LVPOjZm94e4VPZ8PZAReWVkZxNumBZ2bvunc9E7npW+Zfm6SOgKfNm2aL168mKVLlwbyz0URkXSyYsUKVqxYweOPPx5292k9Xw9kCkVERGKXUlMoIiIyeEntAxeRgYn1X8q6O3TqMgMzi+v3VIBL3Lk7rW1Rjp4+y9FTbRw91UbL6TZOt7Vz6mwHp7sep9re/fp0WwdtHVHaOpz2aJT2DqetI0p71Gk//6fT1vVaR7SrFxaIuuMOUaeX5xy6/uz+fOfTfv6YmH6umH/+AZ44yWhvPrKEgrxQXL+nAlwuiLvTePwMtZET1B85zaGWVhqOtXb+2dJK88kzHDnVxtn2aL/fKy9kFOSFKMwPUZAXIj+UQ24oh7yQkZtj5IZyKMjLIXdIbtd253N5OUZOjpFjRo6BYeTkAHRu55hhXX9C9+3OEZCdO8Y4/3WsA6OYx08xfsMLGY/FXmN8R3kSH7k58f/vogCXPh1vbaPmwDG21bew/UALtZETvB05ycmzHX+w36hh+YwZWcDYogJmjRtJcWE+RYV5FBfmc9HQPIoK8ygamsew/FyG5oc6H3kh8kL6CEZkMBTgcl7L6TZerWumOtzEutpmwpET56cDxhYVMHX0cCoqRjGldBhTSoYzftRQxowsiPs/C0UkNgrwLNfQcpqV2w6yclsDm/YdIeowNC/E/MmjWDrnUmaXFTF7XBElw4cEXaqI9JDUAG9paaGyslIX8gSsta2D57c38B+v7WP9niMAlI8dyYOLp7JoaglzJxSTn6vpDZGgnbuQB+j19oa6kCeLNB5r5cnqPfzXhv0cPnmWKSXD+My8Mm6dPZbJJcOCLk9E+tDXhTyaQskCjcdbWfZyHU+/tpe2jig3lY/hnoWTWDT14rj3pYpI8ijAM9jZ9ihPVr/NP7/4Fq3tUT41dxwPLp7KJI22RTKCAjxDvf72YR56dit1TSe5ceZo/uZj5ZomEckwCvAMc6a9g+/9ZjdVr9QxvriQH913JYtnjA66LBFJAAV4Btl/+BR//u9vsKPhGHfOH883PlbOsCH6TyySqfTbnSHW1TbzwH9spK0jyuP3VnBT+ZigSxKRBEtqs++5PvCuvkaJk59t2M/d//oaxYV5/PKBRQpvkQyxYsWKc6sKqQ88E/3b2j18c3kN10wt4Qd3f4iRBXlBlyQicaY+8Az0g5fD/OMLu7i5fAz/fNdchuTqniQi2UQBnqb+fd0e/vGFXdx2xaV85/Y5urOfSBaK6bfezJaY2S4zC5vZQ728XmxmPzezrWb2upnNin+pcs6vtzbwt8truHHmGL6r8BbJWv3+5ptZCHgMuAUoB+40s/Ieu/01sNndPwjcC/xTvAuVTq/WNfOVn26mYmIxj941l1yFt0jWiuW3fz4Qdvc6dz8LPAPc1mOfcuBFAHd/E5hkZmqFiLMDR0/zwNMbGT9qKE/ce6Xuwy2S5WIJ8HHA/m7b9V3PdbcF+DSAmc0HJgJl8ShQOp1p7+CLT2+kta2DH95TQVGhuk1Esl0sAd7b7ep69h5+Cyg2s83Al4BNQHvPgyKRCBUVFecfVVVVF1pv1np4xQ627D/Kd26fw9TRw4MuR0QSrKqq6nxWAiW97RNLF0o9ML7bdhlwoPsO7n4MuA/AOu9P+nbX4w+UlpaiPvALt6rmIE+/to8//8gUbpk9NuhyRCQJKisrz13Eg5k19bZPLCPw9cA0M5tsZvnAHcDy7juY2UVdrwH8GfD7rlCXQWo+cYa/+fk2yseO5Ks3zwi6HBFJIf2OwN293cweBFYBIeBJd68xs/u7Xl8GzAR+bGYdwA7g8wmsOWu4O9/4xXaOnW7nJ382R8ucicgfiOlCHndfCazs8dyybl+vA6bFtzRZsbWB57cf5OtLLufyS0YGXY6IpBgN6VLU8dY2HvnVDuaUFVH5kSlBlyMiKUiX0qeo77/4Fk0nzvDEvRWEcrRupYi8l0bgKSjceJwfVe/hjyvGM2f8RUGXIyIpSvcDTzHuzt8t30FhfoivfVRdJyLZTPcDTzMv7Wrkvh+t5++WlvOniyYHXY6IpIC+7geuKZQUEo0631m1i/GjhnLXgolBlyMiKU4BnkJeqDlIzYFjfOXG6er5FpF+KSVSREfU+d7q3UwbPZzbruh5rzARkfdSgKeIX2x6h3DjCb5683S1DYpITBTgKSAadR57KcwHLh3JRz9wSdDliEiaUICngNU7D1HXdJIvXncZnTdzFBHpn/rAA+buLFtTy4RRhSzR6FtEulEfeIp7/e3D/NEP1/HIbR/gnqsmBV2OiKQg9YGnqB+uqWXUsHw+O298/zuLiHSjAA9QuPEEL77ZyJ9cNYmh+VqgWEQujAI8QD95dS95IeNzCycEXYqIpCEFeEBOnW3n2TfquXX2WEqGDwm6HBFJQwrwgCzffIDjZ9q5e6HueSIiA6MAD4C785PX9jJjzAgqJhYHXY6IpCn1gQdgS30L2985xt0LJ+jCHRHpk/rAU9Bf/mwLK7c18Npf38CIgrygyxGRFKc+8BRx8kw7K7c18Ik5lyq8RWRQFOBJ9vz2g5w628Fn5pUFXYqIpDkFeJI9t7GeiRcX6sNLERk0BXgSvXP0NOvqmvn03DJ9eCkig6YAT6Kfb6zHHT79Ia24IyKDpwBPEnfn2Y3vsGDyKMaPKgy6HBHJAArwJNm0/yhvN53kMx/Sh5ciEh+6kCdJfrWlgfxQDktma9EGEYmNLuRJAdGoc/W3fsfssiIev/c9vfgiIu9LF/IE6I19Rzh4rJWPf3Bs0KWISAZRgCfBr7c2MCQ3hxtmjgm6FBHJIDEFuJktMbNdZhY2s4d6eb3IzFaY2RYzqzGz++JfanrqiDortzWweMZohg/JDbocEckg/Qa4mYWAx4BbgHLgTjMr77HbA8AOd58DXAd818zy41xrWlq/5zCNx8/wMU2fiEicxTICnw+E3b3O3c8CzwC39djHgRHWeXnhcOAw0B7XStPUr7c2UJCXww0zRwddiohkmFgCfBywv9t2fddz3T0KzAQOANuAL7t7NC4VprFo1Hmh5iDXXz6awnxNn4hIfMUS4L3dtKNn7+FHgc3ApcAVwKNmNrLnQZFIhIqKivOPqqqqCyw3vWypP0rk+BluLlfvt4hcmKqqqvNZCZT0tk8sw8J6YHy37TI6R9rd3Qd8yzubysNm9jZwOfB6951KS0vJpj7w1TsOEcoxFs/Q9ImIXJjKyspzF/FgZk297RPLCHw9MM3MJnd9MHkHsLzHPvuAG7reaAwwA6gbYN0ZY/WOQyyYPIqiQi3cICLx1+8I3N3bzexBYBUQAp509xozu7/r9WXAI8BTZraNzimXr7t7r//HyBZ7mk7yVuMJ7lowIehSRCRDxfTJmruvBFb2eG5Zt68PADfHt7T0tnrHIQBu1MU7IpIguhIzQVbvOMTll4zQrWNFJGEU4Alw+ORZNuw9zM3lGn2LSOIowBPgxZ2HiDrcpPZBEUkg3Q88AX678xBjiwqYNe49rfAiIjHT/cCTrLWtg7kPr+az88p45JOzgi5HRDKA7geeJOv3HOZ0WwfXX66Ld0QksRTgcfbyrgj5uTksnHJx0KWISIZTgMfZy7saWTB5FEPzQ0GXIiIZTgEeR/sPn6I2cpLrdO8TEUkCBXgcrdkdAeDa6aUBVyIi2UABHkcv74pQVjyUy0qHBV2KiGQB9YHHyZn2DtbWNnHdjFI6FyYSERkc9YEnSXW4ic898RqP31vBTbqEXkTiSH3gCbZmd4T8UA5XX6b2QRFJDgV4nLy8q5ErJxczbIjWvhSR5FCAx8GBo6fZfegE101X+6CIJI8CPA6qw52LD10zrdd1R0VEEkIBHgdra5u5eFg+M8aMCLoUEckiCvBBcneqw01cddnF5OSofVBEkkd94INUGzlB4/EzLJqq6RMRiS/1gSfYv63dwzeX1/D7ry1mwsVa/1JE4k994AlSHW5i3EVDGT9qaNCliEiWUYAPQkfUebWumUVTL9bl8yKSdArwQag50MKx1nbNf4tIIBTgg1AdbgbgKl0+LyIBUIAPwtraJqaPGc7oEQVBlyIiWUgBPkBn2jtYv+cwV1+m6RMRCYYCfIA27TtKa1tUdx8UkcDoQp4BWhtuIsdggVafF5EE0YU8CfKZf1lLe9T55QOLgi5FRDKcLuSJoxNn2tmy/yiLNH0iIgFSgA/A62830x519X+LSKBiCnAzW2Jmu8wsbGYP9fL618xsc9dju5l1mNmo+JebGqrDzeTn5jBvYnHQpYhIFus3wM0sBDwG3AKUA3eaWXn3fdz92+5+hbtfAfwVsMbdDyeg3pSwtraZeROKKcgLBV2KiGSxWEbg84Gwu9e5+1ngGeC299n/TuA/41FcKmo+cYadDcdYNFXz3yISrFgCfBywv9t2fddz72FmhcAS4NnBl5aa1tV1Xj5/tea/RSRgsSyh3ttt9vrqPVwKVPc1fRKJRKioeLcTprKy8lyPY9qoDjczYkguHxzXa1umiEhcVFVVUVVVdW6z1xFjLAFeD4zvtl0GHOhj3zt4n+mT0tJS0r0PfG1tEwumjCI3pAYeEUmc7gNcM2vqbZ9YUmg9MM3MJptZPp0hvbznTmZWBFwL/HLAFae4+iOn2Nt8iqt0/xMRSQH9jsDdvd3MHgRWASHgSXevMbP7u15f1rXrp4DfuPvJhFUbsLW1nfPf+gBTRFJBLFMouPtKYGWP55b12H4KeCpehaWiteEmSobnM2PMiKBLERHRlZixcneqa5u56rISLZ8mIilBAR6jcOMJIsfP6P4nIpIyFOAxqg53fgis+5+ISKrQ/cBjtLa2mbLioYwfVRh0KSKSJXQ/8DjoiDpzH/4Nt8way//97AeDLkdEsozuBz4I299p4VhrO1erfVBEUogCPAbVtZ3z31rAWERSiQI8BmvDzcwYM4LSEUOCLkVE5DwFeD9a2zpYv+cwV6l9UERSjAK8H5v2HeVMe1TtgyKSchTg/Vhb20SOwYIpGbtCnIikKfWB96M63MQHyy5iZEFe0KWISJZRH/ggHG9t44qHV3P/tVP42kcvD7ocEclS6gMfgNffPkxH1NU+KCIpSQH+PqrDzQzJzWHexOKgSxEReQ8F+PuoDjdx5aRRFOSFgi5FROQ9FOB9aDzeyq5Dx9U+KCIpSwHeh7XhzuXTrlGAi0iKUoD3oTrcxEWFeZRfOjLoUkREeqU+8F64O9XhJq6+7GJCOVo+TUSCoT7wAaiLnOD6767hHz41i88tmBh0OSKS5dQHfgHOLZ+m+W8RSWUK8F78d7iJsuKhTNDyaSKSwhTgPXREnXW1zSy6rAQzzX+LSOpSgPdwbvm0RdM0fSIiqU0B3sN/h88tn6YFHEQktSnAe6gONzFz7EhKhmv5NBFJbeoD76a1rYMNe49wjVafF5EUoD7wC7Bmd4Q/efJ1nrrvSq6bMTrockREAPWBx2TNrghDcnNYOEUjcBFJfQrwbtbsbmThlIt1+1gRSQsK8C77D5+iNnKSa6eXBl2KiEhMYgpwM1tiZrvMLGxmD/Wxz3VmttnMasxsTXzLTLw1uyMAXDtDAS4i6SG3vx3MLAQ8BtwE1APrzWy5u+/ots9FwA+AJe6+z8zS7hPANbsjlBUPZUrJsKBLERGJSSwj8PlA2N3r3P0s8AxwW4997gKec/d9AO7eGN8yE+tse5S14SaunV6qy+dFJG3EEuDjgP3dtuu7nutuOlBsZi+b2Rtmdm+8CkyGN/Ye4eTZDs1/i0ha6XcKBehtSNqzeTwXmAfcAAwF1pnZq+6+u/tOkUiEiop3WxkrKyvPNakHas3uCLk5xtW6fayIpIiqqiqqqqrObfYaTrEEeD0wvtt2GXCgl32a3P0kcNLMfg/MAf4gwEtLS0nFC3nW7I5QMamY4UNiOR0iIonXfYBrZk297RPLFMp6YJqZTTazfOAOYHmPfX4JfNjMcs2sEFgA7Bxw5Ul06FgrOxuOce30tPvcVUSyXL9DTndvN7MHgVVACHjS3WvM7P6u15e5+04zewHYCkSBJ9x9eyILj5eXd3V+3qr5bxFJNzHNGbj7SmBlj+eW9dj+NvDt+JWWHL/d2cilRQXMHDsi6FJERC5IVl+J2drWwStvRbixfIzaB0Uk7WR1gK+tbaK1LcoNM8cEXYqIyAXL6gBfvaORYfkhFk4ZFXQpIiIXLGsXdIhGnd+9eYhrZ5QyJFd3HxSR1KMFHfqwtf4on3i0mu/ePofPzCsLuhwRkT5pQYcefruzkRyDxZer/1tE0lPWBvjqHYeYN7GYUcPygy5FRGRAsjLA9zafZGfDMW4uvyToUkREBiwrA3zltoMA3DJbAS4i6SsrA/z57Q3MGX8RZcWFQZciIjJgWRfg+w+fYmt9C7fO0uhbRNJb1vWBv7C9a/pk1tjAahARiYX6wHv45GPVtEej/OpLHw60DhGRWKkPHHjn6Gk27z+q0beIZISsCvDntzUAcOtsBbiIpL+sCvDnNr7DnLIiJpcMC7oUEZFBy5oA33XwODsajvGpueOCLkVEJC6yJsCf21RPKMdYOufSoEsREYmLrAjwjqjzy00HuHZ6KRcPHxJ0OSIicZEVfeCv1jVz8Firpk9EJK2oDxz4y59tYdX2g6z/xo0U5GnxBhFJL1nbB37iTDvPb2vgltmXKLxFJKNkfIAv33yAk2c7+OMrJwRdiohIXGV0gLs7T7+2l8svGcGHJlwUdDkiInGV0QG+pb6FmgPH+NzCiZhZ0OWIiMRVRgf406/upTA/xCevUO+3iGSejA3wltNtrNh6gNuuuJQRBXlBlyMiEncZ2wf+3MZ6Wtui3DV/YsLfS0QkEbKyD7wj6lz3nZcoGT6En//PRQl/PxGRRMqqPvBVNQfZf/g0lR+eEnQpIiIJk5EB/vgrdUwYVcjNH9C6lyKSuTIuwN/Ye5hN+47y+WsmE8pR66CIZK6YAtzMlpjZLjMLm9lDvbx+nZm1mNnmrsffxr/U2Hz/xTDFhXncXlEWVAkiIkmR298OZhYCHgNuAuqB9Wa23N139Nj1FXf/eAJqjNmmfUdYszvC15dcTmF+vz+aiEhai2UEPh8Iu3udu58FngFuS2xZA/NPL75FcWEe916l1kERyXyxBPg4YH+37fqu53q6ysy2mNnzZvaBuFR3ATbuO8LLuyJUfuQyhg3R6FtEMl8sSdfbJ4E9m8c3AhPd/YSZ3Qr8ApjW86BIJEJFxbutjJWVleea1AfF3fnfv95JyfAhGn2LSEaoqqqiqqrq3GZJb/vEEuD1wPhu22XAge47uPuxbl+vNLMfmFmJuzd136+0tJREXMizquYgG/Ye4f98erZG3yKSEboPcM2sqbd9YplCWQ9MM7PJZpYP3AEs776DmV1iXbf7M7P5Xd+3eRC1x+xse5RvPf8m08cM5/Z56jwRkezR73DV3dvN7EFgFRACnnT3GjO7v+v1ZcBngS+aWTtwGrjDk3SN/uOv1LGn+RQ/uu9KckMZ19YuItKnmOYb3H0lsLLHc8u6ff0o8Gh8S+vf3uaTfP/Ft7hl1iUsnjE62W8vIhKotB2yujvf+MV28kI5fHNp0pteREQCl7YB/sz6/bzyVhNf++gMLikqCLocEZGkS8sAr4uc4OEVO7hmagn3LFTboIhkp7Rb0OFMewd/8dPNDMnL4Tu3zyFHN6wSkQyVUQs6uDsPPbuNn27Yz7K757Fklm4XKyKZLyMWdPjxur38dMN+vnT9VIW3iGS9tAnwldsa+PsVNdw4czRfuXF60OWIiAQuLQJ8ze4IX35mE3MnFPP9O+dq3ltEhDQI8LPtUb7xi21MGz2CJ//0St3nW0SkS8qnYX5uDj/+HwsYUZBL0dC8oMsREUkZKR/gAJNLhgVdgohIykm7PnARkWyRUX3gIiLZKCP6wEVE5F0KcBGRNBVIgHdb50160Lnpm85N73Re+pbp5yatAnwwH34O9NhkH5fsc5Ps4wZzbKafm2Sfl8G8Z7ocl+lZk1ZTKOlyUgdz3ECl08+ncxPf4wYjXX7GdPk7M5hjB3JcUrtQzCwC7AVKgF5XWe5HEdAywLcf6LHJPi7Z5ybZxw3m2Ew/N8k+L4N5z3Q5LhOypgnA3Zf0fCGpAS4iIvGTVlMoIiLyLgW4iEiaSnqAm9kSM9tlZmEzeyjZ75+qzOxJM2s0s+1B15JKzGy8mb1kZjvNrMbMvhx0TanCzArM7HUz29J1bv4+6JpSiZmFzGyTmf0q6FoSJakBbmYh4DHgFqAcuNPMypNZQwp7CnjPhxRCO/BVd58JLAQe0N+Z884A17v7HOAKYImZLQy2pJTyZWBn0EUkUrJH4POBsLvXuftZ4BngtiTXkJLc/ffA4aDrSDXu3uDuG7u+Pk7nL+S4YKtKDd7pRNdmXtdDXQmAmZUBHwOeCLqWREp2gI8D9nfbrke/jBIjM5sEzAVeC7iUlNE1TbAZaARWu7vOTaf/B/wvIBpwHQmV7ADvbS00jRikX2Y2HHgW+At3PxZ0PanC3Tvc/QqgDJhvZrMCLilwZvZxoNHd3wi6lkRLdoDXA+O7bZcBB5Jcg6QZM8ujM7yfdvfngq4nFbn7UeBl9DkKwCLgE2a2h85p2uvN7CfBlpQYyQ7w9cA0M5tsZvnAHcDyJNcgacTMDPhXYKe7fy/oelKJmZWa2UVdXw8FbgTeDLSoFODuf+XuZe4+ic6M+Z273x1wWQmR1AB393bgQWAVnR9G/Ze71ySzhlRlZv8JrANmmFm9mX0+6JpSxCLgHjpHUZu7HrcGXVSKGAu8ZGZb6RwcrXb3jG2Zk/fSpfQiImlKV2KKiKQpBbiISJpSgIuIpCkFuIhImlKAi4ikKQW4iEiaUoCLiKQpBbiISJr6/wkrjT93hSHnAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(A[20:], ms.U[20:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb97efad7f0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD1CAYAAACrz7WZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAApsklEQVR4nO3de3SU533g8e9Po/sd0AgGIcxNBmTAYAQmcW6+pAXbCm2PT2unjZu0RaW206Tbc7LJnt327J7dJt1us417XDsi8WZ9EsebJk4LDY6dxNc6xkZcjAEhkAQGgSQkQDd0l377x8zg8aBBr6TRvHP5fc6Zg+Z9n3fe3/ta/s2j53ne5xFVxRhjTHJJczsAY4wx0WfJ3RhjkpAld2OMSUKW3I0xJglZcjfGmCRkyd0YY5JQutsBAJSUlOiSJUscl+/o6MDr9U75PN3d3RQVFcXkuFieC2J7T2J9bdM5LhF+R6Z7nP2ORO+46dyTeLq2AwcOdKrqxBegqq6/VqxYoTt27NDdu3erExs3bnRULtyOHTtidlwsz6Ua23sS62ubznGJ8Dsy3ePsdyR6x03nnsTDte3evVt37NihwCmNkFfjouZeVFREbW3trJ+nuro6ZsfF8lwzkQjXFst7kgjXZr8j0TsulueK5nHV1dVUV1eza9eu7kjHicbBE6pVVVVaV1c3lfJMpXwqsHvyYXY/rmf35HqJfk9E5ICqVk20LyE7VGtqatwOIe7YPfkwux/Xs3tyvWS+JwlZczfGGJMANffu7m5qamrYs2eP26EYY0zc27NnT/CvjojDbxzV3EVkK/AtwAN8R1W/EbZfAvvvBfqBz6vqwcC+vwD+BFDgPeALqjoYerzV3I0xZupmVHMXEQ/wBLANqAQeEpHKsGLbgIrAqwZ4MnBsGfDnQJWqrsH/5fDgNK8jYQWHJhljTKw4GQq5GWhU1WYAEXkO2A4cDymzHXhG/Rlsn4gUi4gv5Bw5IjIC5AIXohZ9nOsfHuV//KyenxxsITvDwxc+upTH7lqBJ03cDs0Yk+SctLmXAedC3rcEtk1aRlXPA/8LOAu0At2q+tL0w00cY+NKzTMHePads1SvW8imJXP53788yX/516Nuh2aMSQFOau4TVTPD2xgmLCMic/DX6pcCXcA/i8gfqOr3Qwt2dHRQVfVBs1FNTU3CD1H63q/P8O+NnXzjd9by4ObFAHz9hXq+/Vozdywv4b51vkk+wRhjrldbWxv60GdJpHJOknsLUB7yfhHXN61EKnMPcFpVOwBE5Hngo8CHkrvX603oBwnCtXYP8PcvNXDnSi+/t+mD2/KV31zFGyc7+foL9dxTWUpWusfFKI0xiSi08isinZHKOWmW2Q9UiMhSEcnE3yG6O6zMbuBh8duCv/mlFX9zzBYRyQ2MqLkbqJ/65SSW775xmqHRcf7b9jX4L9vPkyZ87d5VtFwZ4Pv7zroYoTEm2U2a3FV1FHgMeBF/Yv6Rqh4TkZ0isjNQbC/QDDQCu4BHAse+DfwYOIh/GGQaMPuTyLioZ3CE5/af4761Psrn5l63/+MVXrYsm8vT/36asXEbQWOMmR2OJg5T1b34E3jotqdCflbg0QjH/jXw1zOIMaH8uK6FvqFRdnx8WcQyD39kCY/84CCvNlzk7tXzYxidMSZV2BOqUfYvh8+ztqyItYsiz9v86cr5lBZk8f1978cwMmNMsnDyhGpcJPfglL+xnso02k53XuVISzefuXXhDctleNJ4YOMiXj/VyeWrwzGKzhiTLKqrq4MjZiJO+RsXyT1Z7D58ARG4/9bJhznet87H2Ljy4rG2GERmjEk1ltyj6OfH2th001x8RTmTlq30FbJkXi4/O9Iag8iMManGknuUtHYPUN/aw92rSx2VFxHuXevjreZLXLGmGWNMlMVFck+GDtVXGzoAuHOVs+QO/o7VsXHl9VMdsxWWMSYJWYdqDL1y4iJlxTlUlOY7PmbdomLm5GbwWoMld2OMc9ahGiMjY+O82djJJ1d6P/RE6mQ8acInbvby2skOxu2BJmNMFFlyj4L3zndzdXiMj62IOIdPRJ9a6eXS1WGOXoj4BWyMMVNmyT0K9jVfAmDz0rlTPvZjK7wA/LrpUlRjMsaktrhI7oneobqv+TI3z8+nJD9rysd6C7JY7s3j7WZL7sYYZ6xDNQZGxsapO3OZLcvmTfszbl82j7ozV2wiMWOMI9ahGgNHz3fTPzzG7UtnkNyXzqV3aJTjF3qiGJkxJpVZcp+hQ2e7ANh405xpf0bwi+Ht09Y0Y4yJDkvuM3T4XBe+omwWFGVP+zMWFGVz07xc3j59OYqRGWNSWVwk90TuUD107gobFhfP+HM2L5nL/jOXbby7MWZSUetQFZGtItIgIo0i8tUJ9ouIPB7Yf0REbgtsXykih0NePSLy5fDjE7VDtaN3iHOXB9hQPv0mmaDbl82jq3+Ekxd7oxCZMSaZRaVDVUQ8wBPANqASeEhEKsOKbQMqAq8a4EkAVW1Q1fWquh7YCPQDP53ylcSpw+e6AKJSc78t8BnBNnxjjJkJJzX3zUCjqjar6jDwHLA9rMx24Bn12wcUi0j4pOZ3A02qmjTLDx06e4X0NGFNWeRVl5xaWpJHUU4Ghy25G2OiwElyLwPOhbxvCWybapkHgR9OdIKOjg6qqqquvQJ/bsS9Q2e7qFxYSHaGZ8afJSKsLy++9teAMcZMpLa29lquBCLOeeJkgeyJZsIK7/W7YRkRyQQ+A3xtohN4vV7q6uochBI/xsaVIy1dPLBxUdQ+c315MY+fOkXf0Cj5WY7WLjfGpJiamppgZyoi0hmpnJOaewtQHvJ+EXBhimW2AQdVtd3B+RLC6c4+rg6PsXZRcdQ+c8PiYlThSEtX1D7TGJOanCT3/UCFiCwN1MAfBHaHldkNPBwYNbMF6FbV0PXjHiJCk0yiOnre/zTp2ii0twetLy8GrFPVGDNzk/7tr6qjIvIY8CLgAZ5W1WMisjOw/ylgL3Av0Ih/RMwXgseLSC7waeBPox++e46e7yYrPY3l3ryofWZxbiZLS/Ks3d0YM2OOGnZVdS/+BB667amQnxV4NMKx/cANJ14JPsRUXV2dMGPdj17oZpWvkHRPdJ8D21BezBuNnajqlBb+MMakjj179gQf+rRZIaNJVTl2oYc1Cwuj/tm3lhfT0TtEW89g1D/bGJMcbFbIWXLu8gC9g6NRGd8ebk2Z/wsj2KZvjDHTYcl9GoJL4t0yCzX31b5C0sTfpm+MMdNlyX0ajp7vJj1NuHl+QdQ/OzczneXefI7ZmqrGmBmIi+SeaLNCHrvQQ8X8gqg8mTqRNWVF1ixjjInIltmbJccu9FDpi36TTNAtCwtp6xmko3do1s5hjElc1qE6Cy71DdHZN8RqX/SbZIKCHbXWNGOMmS5L7lN0sr0PYFba24MqFwZHzFhyN8ZMjyX3KTrZ7l9MY+WC2UvuhdkZLJmXa+3uxphpi4vknkgdqifbeynMTqe0IGtWz7OmrOjakEtjjAllHaqz4GR7LysXFMz61ABryopouTJAV//wrJ7HGJN4rEM1ylSVhrbeWW1vD1qzMNipak0zxpips+Q+Be09Q/QMjs5qe3tQ8OnX96xT1RgzDZbcp6Ah0Jkai5r7nLxMFhZlU99qNXdjzNTFRXJPlA7Vk22xS+7gn2fGkrsxJlzUOlRFZKuINIhIo4h8dYL9IiKPB/YfEZHbQvYVi8iPReSEiNSLyEfCj0+UDtWG9l68BVnMzcuMyflW+wpp6rjK4MhYTM5njEkMUelQFREP8AT+dVArgYdEpDKs2DagIvCqAZ4M2fct4Oequgq4FaifwjXElVPtvayMUa0dYJWvgLFxpfFiX8zOaYxJDk5q7puBRlVtVtVh4Dlge1iZ7cAz6rcPKBYRn4gUAp8AvgugqsOq2hW98GNnfFw52d5Hxfz8mJ1zdWD+GmuaMcZMlZPkXgacC3nfEtjmpMwyoAP4PyJySES+IyLRW3Q0hlquDDAwMhbTmvuSeXlkZ6RR39obs3MaY5KDk+Q+0dM66rBMOnAb8KSqbgCuAte12Xd0dFBVVXXtFWhLiivXRsrEYBhkkCdNWLnAOlWNMR+ora29liuBkkjlnCyQ3QKUh7xfBFxwWEaBFlV9O7D9x0yQ3L1eL3V1dQ5CcU9wTpmK0tg1ywBU+gp44WibLZhtjAGgpqYmOFIGEemMVM5JzX0/UCEiS0UkE3gQ2B1WZjfwcGDUzBagW1VbVbUNOCciKwPl7gaOT/Fa4kJDWy9lxTkUZGfE9LyrfYV09Y/YgtnGmCmZtOauqqMi8hjwIuABnlbVYyKyM7D/KWAvcC/QCPQDXwj5iC8CPwh8MTSH7UsYwTllYi20U9VXlBPz8xtjEpOTZhlUdS/+BB667amQnxV4NMKxh4Gq6YfovpGxcZo7rvKplaUxP/eqwBdKfWsvd62aH/PzG2MSkz2h6sD7l64yPDbOzTEcBhlUkJ1B+dwcjlunqjEmwMkTqo5q7rMt+IRqvGpom/3Vl25k9YJCTlhyN8YEVFdXU11dza5du2zK35loaO8lTWBFjEfKBK3yFXK606YhMMY4Z8ndgZNtvYEHijyunL/SV8C4+kfsGGOME5bcHTjZHpsFOiKxaQiMMVMVF8k9njtUB0fGOHPpakyfTA1XPieXvEyPJXdjDGAdqlHR1NHHuBLTOWXCpaUJq3yFNseMMQawDtWoCE47sHKBO52pQat9BdS39eB/pMAYY27MkvskGtr6yPAIN81zdzLL1b5CegdHabky4GocxpjEYMl9Eifbe1nuzSfD4+6tCnaqnrARM8YYB+Iiucdzh2pDm7sjZYJWzi9AxEbMGGOiuIbqbIvXNVT7hkY53zXgyoRh4fKy0rlpbq4ld2NMdNZQTWWnggt0xEHNHaByoS3cYYxxxpL7DVwbKRMnyX31gkLev9xP39Co26EYY+KcJfcbaGjrIyfDw6I58TGP+mpfIarQ0Ga1d2PMjcVFco/XDtWT7b1UzM8nLS0+lrerXOgfMXPcHmYyJqVFrUNVRLaKSIOINIrIdWugBpbXezyw/4iI3Bay74yIvCcih0VkwoVS47VDtcHlOWXC+YqyKcrJ4PgFq7kbk8qcdKhOOv2AiHiAJ4BP418Ie7+I7FbV0LVQtwEVgdftwJOBf4PuVNWIC7nGo8tXh+noHYqb9nYAEfE/qWqdqsaYSTipuW8GGlW1WVWHgeeA7WFltgPPqN8+oFhEfFGONaaCnaluThg2kUpfEQ1tvYyN2zQExpjInCT3MuBcyPuWwDanZRR4SUQOiEjNdAONtVNxNlImaLWvgIHATJXGGBOJk1khJ+pNDK823qjMHap6QURKgV+IyAlVfT20YEdHB1VVH6yhXVNTE+wscE1Dey+F2enML8xyNY5woXO7L/e6O5mZMSb2amtrQ2fRLYlUzklybwHKQ94vAi44LaOqwX8vishP8TfzfCi5e71e6uom7Gt1zcm2PlYuKEAkPkbKBFXMzyc9Tahv7eH+dQvdDscYE2OhlV8RidiX6aRZZj9QISJLRSQTeBDYHVZmN/BwYNTMFqBbVVtFJE9ECgJB5AG/ARyd+uXElqrS0N5LRZw1yQBkpXtYUZpvI2aMMTc0ac1dVUdF5DHgRcADPK2qx0RkZ2D/U8Be4F6gEegHvhA4fD7w00DtNx14VlV/HvWriLKLvUN0D4zEXXt70GpfIW81XXI7DGNMHHO0EpOq7sWfwEO3PRXyswKPTnBcM3DrZJ8ffIgpuLqI24ILUcfTGPdQlb5CfnroPJevDjM3L9PtcIwxMbZnz57gQ582K+RUXBsGOT8+OyxtwWxjUpvNCjlNJ9t7KcnPYl5+fI2UCVrt8/9FYcndGBOJJfcJNLT3ub5m6o3My89ifmGWdaoaYyKy5B5mfFw5FWdzykxkta+Q41ZzN8ZEEBfJPZ5mhTzfNUD/8FjcjpQJWu0rpKmjj+HRcbdDMcbEmC2zNw3BkTLxOMY9VKWvkJEx5dRFm/7XmFRjHarT0BDnI2WCPhgxY8ndGHM9S+5hGtp6KSvOoSA7w+1QbmhpSR7ZGWnWqWqMmZAl9zANbb2sjLNpfifiSRNWLrAFs40xE4uL5B4vHarDo+M0dfQlRHIHqPQVUN/Wg/8BYWNMqrAO1Slq7uxjdFxZlTDJvZCu/hEudA+6HYoxJoasQ3WKgiNlEqbmvtD/pX30fMT/vsaYFGXJPcSJtl7S04RlJfE9Uiao0ldImlhyN8Zcz5J7iIa2XpZ788lMT4zbkpPpoaK0gPcsuRtjwsRFFouXDtVEGSkTak1ZEUfPd1unqjEpJGodqiKyVUQaRKRRRL46wX4RkccD+4+IyG1h+z0ickhE/m2iz4+HDtWewRHOdw0kXHJfW1ZIZ98w7T1DbodijImRqHSoiogHeALYBlQCD4lIZVixbUBF4FUDPBm2/0tAvePIXXAy0JmaKCNlgtYu8n9xW9OMMSaUk5r7ZqBRVZtVdRh4DtgeVmY78Iz67QOKRcQHICKLgPuA70Qx7qg7kWAjZYIqfUWkiSV3Y8yHOUnuZcC5kPctgW1Oy/wD8BUgrqcvbGjrpSArnbLiHLdDmZKcTP+C2TZixhgTyklylwm2hffeTVhGRO4HLqrqgRudoKOjg6qqqmuvQFtSTDW09XLzggICi3knlDVlRVZzNyZF1NbWXsuVQEmkck4WyG4BykPeLwIuOCzzAPAZEbkXyAYKReT7qvoHoQd7vV7q6uochDI7VJUTbT3cf+tC12KYiTULi3j+4HnaewaZX5jtdjjGmFlUU1MTHCmDiHRGKuek5r4fqBCRpSKSCTwI7A4rsxt4ODBqZgvQraqtqvo1VV2kqksCx70cntjjQVvPID2DownXmRp0rVO1xWrvxhi/SZO7qo4CjwEv4h/x8iNVPSYiO0VkZ6DYXqAZaAR2AY/MUryz4lpnapwv0BFJpa8QsU5VY0wIJ80yqOpe/Ak8dNtTIT8r8Ogkn/Eq8OqUI4yBD4ZBFrocyfTkZaWz3GudqsaYD9gTqvg7UxcUZlOUG98LdNzI2rIijl6w5G5MKrApfx06kYDTDoRbU1ZEe88QF3ts+l9jkp1N+evA6Ng4jR19CduZGrQu0Kn6rnWqGmOw5M6ZS1cZHh1P/Jr7wiI8acLhc1fcDsUYEwdSPrkfCywwXbkwMTtTg3IyPaz2FXD4XJfboRhj4kBcJHc3O1SPt/aQ6UljuTcxFui4kfXlxbx7rpuxcZv+15hkZh2qDhy/0MPNC/LJ8MTFrZiRDeVz6Bsapamjz+1QjDGzyDpUJ6GqHL/Qwy2+iF9+CWX94mIADp/tcjUOY4z7Ujq5X+wd4tLV4YRvbw9aOi+PopwMDlmnqjEpL6WT+7HAQz/JktzT0oRby4s5ZDV3Y1JeXCR3tzpUjwdGyiT6GPdQ68uLOdney9WhUbdDMcbMEutQncTx1h6WzMulIDtxpx0It2FxMeMKR+xhJmOSlnWoTuLYhZ6kaZIJWr+oGMDGuxuT4lI2ufcOjvD+pX4qfcmV3OfkZbK0JI+DZ61T1ZhUlrLJvb7VP83vLQuTYxhkqI03zaHuzGX8MzEbY1JRXCR3NzpUjyfZSJlQm5fM5Ur/iD3MZEySilqHqohsFZEGEWkUka9OsF9E5PHA/iMicltge7aIvCMi74rIMRH5rxN9vhsdqsdbe5iXl0lpQVbMzhkrm5bOBeCd09Y0Y0wyikqHqoh4gCeAbUAl8JCIVIYV2wZUBF41wJOB7UPAXap6K7Ae2BpYY9V1wc5UEXE7lKhbMi+XkvxM6s5cdjsUY4xLnNTcNwONqtqsqsPAc8D2sDLbgWfUbx9QLCK+wPtg20BG4OV6Q/Dw6Din2vuSskkGQETYtGQu71hyNyZlOUnuZcC5kPctgW2OyoiIR0QOAxeBX6jq29OONkqaOvoYHhtPupEyoTYtmUvLlQFauwfcDsUY4wInyX2idovw2nfEMqo6pqrrgUXAZhFZE16wo6ODqqqqa69AW9KseS+wkHQyjpQJ2rTE3+6+/4y1uxuTTGpra6/lSqAkUrl0B5/VApSHvF8EXJhqGVXtEpFXga3A0dB9Xq+Xuro6B6FEx3st3eRnpbOsJC9m54y11b4C8jI97D99mc/cutDtcIwxUVJTUxMcKYOIdEYq56Tmvh+oEJGlIpIJPAjsDiuzG3g4MGpmC9Ctqq0i4hWR4kAQOcA9wIkpX02UHTnfzZqyQtLSkq8zNSjdk8ZtN81hv7W7G5OSJk3uqjoKPAa8CNQDP1LVYyKyU0R2BortBZqBRmAX8Ehguw94RUSO4P+S+IWq/luUr2FKhkfHqW/tYV3gMf1ktnnJXBrae+nuH3E7FGNMjDlplkFV9+JP4KHbngr5WYFHJzjuCLBhhjFG1cn2XoZHx1lblrzt7UG3L5uHKrzVfImtaxa4HY4xJoZS7gnV4GyJ6xYlf3JfX15MToaHXzdFbJYzxiQgJ0+oOqq5z7bgE6qx8N75LopyMlg8Nzcm53NTZnoaty+by5uNltyNSSbV1dVUV1eza9cum/I36EhLN+sWFSXlk6kTuWN5CU0dV2nrHnQ7FGNMDKVUch8cGaOhrTcl2tuD7ljhHwZrtXdjUktKJff61h5GxzUl2tuDVi0oYG5eJm9au7sxKSUuknusOlSDT6auTYFhkEFpacJHls/jzcZOm9/dmCRha6iGOdLSTUl+JguLsmf1PPHmYytKaO8ZoqnjqtuhGGOiwNZQDfNeSzdry1KnMzXojuXW7m5MqkmZ5N4/PMqpi70p1SQTtHheLovn5vL6yQ63QzHGxEhcJPdYtLkfu9DDuMK6FBopE+quVaW82dTJ4MiY26EYY2bI2txDvHuuC0iNJ1MncueqUgZHxnmr6ZLboRhjZsja3EMcOttFWXEOpYWp1ZkadPvSueRmenj5xEW3QzHGxEDKJPeDZ6+wYXGx22G4JjvDwx0rSnj5xEUbEmlMCkiJ5N7aPUBr9yC3LZ7jdiiuuntVKee7Bmho73U7FGPMLIuL5D7bHaqHz3YBpHTNHfzt7oA1zRiT4KxDNeDg2Stkpqcl9ZqpTswvzGZNWSEv11tyNyaRRa1DVUS2ikiDiDSKyFcn2C8i8nhg/xERuS2wvVxEXhGRehE5JiJfmu7FzMShs12sWVhIZnpcfJe56q5V8zl49gqX+obcDsUYM4smzXYi4gGeALYBlcBDIlIZVmwbUBF41QBPBraPAn+pqquBLcCjExw7q4ZHxzlyvjvl29uDtt6ygHGFF4+1ux2KMWYWOanKbgYaVbVZVYeB54DtYWW2A8+o3z6gWER8qtqqqgcBVLUX/xqsZVGMf1L1rT0Mj46zwZI7AKt9BSwryeNn711wOxRjzCxyktzLgHMh71u4PkFPWkZEluBfT/Xt8BN0dHRQVVV17RXNVZkOnr0CWGdqkIhw3zofbzVdotOaZoxJOLW1tddyJVASqZyTZfYmmmUrfKD0DcuISD7wE+DLqtoTXtDr9VJXV+cglKk7dLaLBYXZLCzOmZXPT0T3rvXxjy838uKxNn7/9pvcDscYMwU1NTXBkTKISMTZAJ3U3FuA8pD3i4Dwv+kjlhGRDPyJ/Qeq+ryD80VVqj+8NJFVCwpY5s3jZ0da3Q7FGDNLnCT3/UCFiCwVkUzgQWB3WJndwMOBUTNbgG5VbRX/3LrfBepV9ZtRjdyBi72DtFwZsM7UMCLCfWt97Gu2phljktWkyV1VR4HHgBfxd4j+SFWPichOEdkZKLYXaAYagV3AI4HtdwCfA+4SkcOB173RvohI7OGlyO5b52Nc4YWjbW6HYoyZBU7a3FHVvfgTeOi2p0J+VuDRCY77dyZuj/+Q4BOq1dXVUX2Q6cD7V8jwCGtSdJrfG1k5v4CK0nx+erCFz22xdndjEsmePXuCT/RHTG4SD5NIVVVV6Wx0qP7OP72JiPCTP/to1D87GXz7tSa+/sIJfvWXn2S5N9/tcIwxUyQiB1S1aqJ9SfvI5uDIGO+d76ZqibW3R/LbG8pIE/jJgRa3QzHGRFnSJvdDZ7sYGVM2L5nrdihxq7Qwm0/e7OX5g+cZG3f/LzhjTPQkbXKvO3MZgKqbLLnfyAMby2nrGeSNU7a+qjHJJC6S+2xM+fvOmcusnF9AUW5G1D4zGd1TWcrcvEyeffus26EYYxxK2Sl/R8fGOfj+FTYttfb2yWSle/i9TeX8sr6d810DbodjjHEgZddQPdHWy9XhMTZZe7sjn928GAV+aLV3Y5JGUib3d07729stuTtTPjeXu1eV8tz+swyPjrsdjjEmCpIyude9f5my4hybLGwKPveRJXT2DbP7XZsK2JhkEBfJPZodqqrKO6evsMnGt0/JJypKWLWggNrXm4iHB9uMMZGlZIfq+5f66ewbYtNSa5KZChHhTz+5jJPtfbzSYGusGhPPUrJD9Z3A+HZ7eGnq7l+3kIVF2Tz1WrPboRhjZij5kvvpy8zJzbC5UqYhw5PGn3x8Ge+cvsy+5ktuh2OMmYGkS+5vNV3i9qXzSEubdDJKM4HP3r6Y0oIsvvnSSWt7NyaBxUVyj1aH6rnL/ZzvGuAjy+dFKbLUk53h4Yt3reCdM5d541TEFbyMMS6KWoeqiGwVkQYRaRSRr06wX0Tk8cD+IyJyW8i+p0XkoogcjfT50epQfSvQlLBlmSX3mfjdTeWUFefw9y81MG4TihkTd6LSoSoiHuAJYBtQCTwkIpVhxbYBFYFXDfBkyL7vAVunEvh07Wu6xLy8TG6eb+3tM5GV7uE/fPpm3m3p5l/fPe92OMaYaXBSc98MNKpqs6oOA88B28PKbAeeUb99QLGI+ABU9XXgcjSDnoiqsq/5EluWzcO/dKuZid/eUMati4r4xgsnuDo06nY4xpgpcpLcy4BzIe9bAtumWmZWnb3cz4XuQbYssyGQ0ZCWJvxV9S209wzx5KtNbodjjJkiJ8l9ompweEOskzIRdXR0UFVVde0VaEuakuDQPetMjZ6NN83htzeU8e3XmzjZ3ut2OMYYoLa29lquBEoilXOyQHYLUB7yfhEQPgGJkzIReb1eZrqG6ltNlyjJz7Lx7VH2n+9bzasNF/mPPznCj3d+FI8NMTXGVTU1NcGRMohIxCFtTmru+4EKEVkqIpnAg8DusDK7gYcDo2a2AN2q2jq90KfO395+mS3L5lp7e5TNy8/ir6orOXS2i+/9+ozb4RhjHJo0uavqKPAY8CJQD/xIVY+JyE4R2RkothdoBhqBXcAjweNF5IfAW8BKEWkRkT+O8jVw5lI/bT2DNgRylvzW+jLuXlXK3/78BPWtPW6HY4xxwEmzDKq6F38CD932VMjPCjwa4diHZhKgE281WXv7bBIR/vaBdWz71ht88YeH2PPYx8jJ9LgdljHmBpLiCdV9zZcoLchiWUlelCMzQSX5WXzzd2+l8WIf/+mn79nUBMa4yMkTqhIP/5NWVVXpdDtUVZXNf/MrPrJsHo8/tCHKkZlwj//qFN/8xUm+snUlj3xqhdvhGJPSROSAqlZNtM9Rs0w8a+68SkfvkDXJxMgX71rBqYt9/N2LDazw5vMbtyxwOyRjzATiollmJoLt7daZGhsiwt89sI51ZUX8+XOHeNumBjYmLiV8ct/XfIkFhdksmZfrdigpIzvDw3c/v4my4hz+6Hv7OfD+FbdDMsaEiYvkPpMO1eOtPWxYXGzj22OsJD+LZ3dswVuQxeeffocD78/69EHGmICU6FC9/W9+yaduLuVvH1gX5aiMExe6Bvjsrn1c6B7kW7+3nm1rfW6HZEzKuFGHalzU3Geif2jMxly7aGFxDs8/cgdrFhbyyLMH+fZrTTZM0pg4kNDJXVXpHxkjL8uSu5vm5mXy7I4tbFuzgK+/cIIdzxzgytVht8MyJqUldHIfHhtnbFzJzUz4EZ0JLzvDwxOfvY2/ur+S105e5N7H3+C1kx1uh2VMyoqL5D7dDtWB4TEAcjKs5h4PRIQ/+thSnv+zO8jJ9PCHT7/Dn//wEBd7B90OzZikkvQdqhe6BvjoN17mG7+zlgc3L56FyMx0DY2O8eSrTfzTK01keIQ//vgydnx8KQXZGW6HZkzSSNoO1f5h//JvuVnWLBNvstI9fPmem/n5lz/OJ1d6efxXp/jE/3yFb77UwMUeq8kbM9sSOiv2B5plcq1ZJm4t8+bzT7+/kfdauvnWr07yj6808uRrTdy31sfvVpVz+7J5tgCIMbMgOZK7DYWMe2sXFfGdP9zEmc6r/N+3zvDPdS38y+ELlBZkcf+6hdxTWcrGm+aQlW7/LY2JhrhI7sEO1erqaqqrqx0fd61D1ZJ7wlhSksdfV9/CV35zFS+fuMi/Hj7P9/e9z9NvniY308OWZfO4Y0UJGxYXU+krJNv+KjPmOnv27AkOQJlZh6qIbAW+BXiA76jqN8L2S2D/vUA/8HlVPejkWJh+h+rPjrTy6LMHeekvPsHN8wumfLyJD31Do+xrusTrpzp441QnpzuvApCeJqzyFXCLr4gVpfksL81jhbeAsjk51pRjDDOc8ldEPMATwKfxL4S9X0R2q+rxkGLbgIrA63bgSeB2h8dOW7BD1YZCJrb8rHTuqZzPPZXzAWjrHuTwuS6OtHRx+FwXv6xv5//VnbtWPtOTxvyiLHyFOSwoysZXlE1pYTZzcjOYk5tJUeDfObkZFGRn2BeBSUlOmmU2A42q2gwgIs8B24HQBL0deCaw3N4+ESkWER+wxMGx02Zt7slpQVE2W4sWsHXNB3PFX7k6THNnH00Xr3L60lVauwZoDXwJ/PzoIMNj4xE/LyfDQ06mh5wMD7mZ1/+cnpZGukfICPybniakez68LcOTRnqaIAJpgUnq0sT/XoC0NEEAREgTEIJlP/hZRAJlP9gWa9OdYG+6oc7kGmXaZ00sa8oKuWle9FeRc5Lcy4BzIe9b8NfOJytT5vBYOjo6qKr64C+Lmpqa4AD9G/ogucdF14GZRXPyMtmYN5eNN829bp+q0tU/QtfACFf6h+nuH6FrYJgrV/3bBkfG6B8epX94jIHhMQZGxugfHqOzb5iBkTFGx8YZGVNGx8cZHVNGAk8+j4wro2PjjLv/KIhJYv/9t9ZMKbnX1tZSW1sbfFsSqZyTrDjR12f4r3ukMk6Oxev1Mp0294c2l3PnKi/ZGQk9XN/MkIgwJy+TOXmZLCX6NaDxcWUkkPgV/5fJuAIK4/rhbYoGtvt/1mAZxf8K2TYdM/memf7zitOMdQbBptL3aWlB1pTKh1Z+RaQzUjknyb0FKA95vwi44LBMpoNjp604N5Pi3MxofZwxE0pLE7LSPNizciaROKny7gcqRGSpiGQCDwK7w8rsBh4Wvy1At6q2OjzWGGNMlE1aF1HVURF5DHgR/3DGp1X1mIjsDOx/CtiLfxhkI/6hkF+40bGzciXGGGOuiYuJwyoqKvTOO++c8kNMxhiTioIPMe3atatRVSsmKhMXyX0my+wZY0yqSrpZIUOGAZkAuycfZvfjenZPrpfM9ySlkvtUFwOZyXGxPBfE9p7E+tqmc1wi/I5M9zj7HYnecdO5J4lybQmZ3KcrEf4HnO5x05UI1xbLe5II12a/I9E7LpbnivVxcdHmLiLdwDDQBXQ7OKQEiDh4/waKHH5+NI6L5bkgtvck1tc2neMS4XdkusfZ70j0jpvOPYmHaysCioFMVZ1wZsi4SO7GGGOiK6WaZYwxJlVYcjfGmCSUcMldRLaKSIOINIrIV92Ox20i8rSIXBSRo27HEg9EpFxEXhGRehE5JiJfcjsmt4lItoi8IyLvBu7Jf3U7pnggIh4ROSQi/+Z2LLMhoZJ7yOIf24BK4CERqXQ3Ktd9D9jqdhBxZBT4S1VdDWwBHrXfEYaAu1T1VmA9sDUwB1Sq+xJQ73YQsyWhkjshC4eo6jAQXPwjZanq68Blt+OIF6raGlziUVV78f/PW+ZuVO5Sv77A24zAK6VHUojIIuA+4DtuxzJbEi25R1oUxJjriMgSYAPwtsuhuC7QBHEYuAj8QlVT/Z78A/AVIPISXgku0ZK7o8U/jBGRfOAnwJdVtcfteNymqmOquh7/mgqbRWSNyyG5RkTuBy6q6gG3Y5lNiZbcnSwcYlKciGTgT+w/UNXn3Y4nnqhqF/Aqqd1PcwfwGRE5g79p9y4R+b67IUVfoiV3W/zD3JD4V4D+LlCvqt90O554ICJeESkO/JwD3AOccDUoF6nq11R1kaouwZ9DXlbVP3A5rKhLqOSuqqNAcPGPeuBHqb74h4j8EHgLWCkiLSLyx27H5LI7gM/hr40dDrzudTsol/mAV0TkCP4K0i9UNSmH/5kP2PQDxhiThBKq5m6MMcYZS+7GGJOELLkbY0wSsuRujDFJyJK7McYkIUvuxhiThCy5G2NMErLkbowxSej/A084hgEkSaZrAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(A[5:], (dfdA(ms.U[5:], ms.deltaL, ms.L, ms.RH) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'HM' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-5106c39e61ec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mtest\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mHM\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mR\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvs2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdeltaL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mL\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'HM' is not defined"
     ]
    }
   ],
   "source": [
    "test = HM(R, m, U, vs2, deltaL, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title ICs that work\n",
    "\n",
    "n=500\n",
    "\n",
    "L = np.linspace(-2, 0.3, n)\n",
    "A = 10**L\n",
    "eps = 0 # some small number for numerical stability\n",
    "\n",
    "deltaL = L[1] - L[0]\n",
    "deltaA = A[1] - A[0]\n",
    "deltau = np.min(deltaA)/5\n",
    "\n",
    "\n",
    "vs2 = 0.05\n",
    "# #np.sinh(A) + 1 maybe did something\n",
    "g = 1 + 0.01*np.exp(A)  # g -> 1 at A=0\n",
    "c = 0.2*A - 0.01*A**2 # 0.1*np.sqrt(A) - 0.025*A + 0.0005*A*np.exp(-(A-0.2)**2/.2**2) # c -> 0 at A=0 ✓\n",
    "\n",
    "U = -np.sqrt(g**2 - 1 + 2*c)\n",
    "m = integ.cumtrapz(np.concatenate(([1], g - vs2*U)), np.concatenate(([0], A))) # m=0 at A=0\n",
    "R = m/c\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gamma(R,m,U): \n",
    "    return np.sqrt(1+U**2-2*m / R)\n",
    "\n",
    "def P(R,m,U,vs2) : \n",
    "    return vs2*rho(R,m,U)\n",
    "\n",
    "def rho(R,m,U) : \n",
    "    return (gamma(R,m,U) + U) / (4*np.pi*R**2*dfdA(R, deltaL, L))\n",
    "\n",
    "def psi(R,m,U,g,r,p) : \n",
    "  return np.log(g[-1]+U[-1]) + integ.cumtrapz(dfdA(U, deltaL, L)/g + m/(4*np.pi*r*(R**4)*g) + p/(r*g*R), A, initial=0)\n",
    "\n",
    "def k_coeffs(R,m,U,vs2) :\n",
    "  g = gamma(R,m,U)\n",
    "  r = rho(R,m,U)\n",
    "  p = P(R,m,U,vs2)\n",
    "  \n",
    "  ep = np.exp(psi(R,m,U,g,r,p))\n",
    "\n",
    "  #print(ep)\n",
    "  kR = ep*U\n",
    "  km = -1/ep*4*np.pi*R**2*p*U\n",
    "  kU = -ep/(1-vs2)*( 4*np.pi*g*R**2*r/(r + p)*dfdA(p, deltaL, L) + (m + 4*np.pi*R**3*p)/R**2 ) \\\n",
    "    - ep*vs2/(1-vs2)*(4*np.pi*r*R**2*dfdA(U, deltaL, L) + 2*U*g/R)\n",
    "  #kU = 1\n",
    "  return kR, km, kU\n",
    "\n",
    "def cfl_deltau(R,m,U,vs2) :\n",
    "    g = gamma(R,m,U)\n",
    "    r = rho(R,m,U)\n",
    "    p = P(R,m,U,vs2)\n",
    "    ep = np.exp(psi(R,m,U,g,r,p))\n",
    "    Rp = ep * U\n",
    "    return np.abs(Rp / (U + g) / ep).min() * (1/ np.sqrt(1/3) - 1 )\n"
   ]
  }
 ],
 "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": 5
}