ms_hm_universal.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate as integ\n",
    "import scipy.interpolate as interp\n",
    "\n",
    "\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": [
    "# # B. Fornberg, \"Calculation of weights in finite difference formulas\"\n",
    "# https://epubs.siam.org/doi/10.1137/S0036144596322507\n",
    "# m is the the order of the derative; x is your grid spacing domain\n",
    "#This will create an array for  a point \"xbar\" in refrance to all other domain points. for a domain which is equally \n",
    "#spaced x=0,1,2 for first order (m=1) a array will be created for point a \"xbar\" coefficient for each point for the domain.\n",
    "#For example if xbar was 1 then the array [-0.5,0,0.5] \n",
    "\n",
    "def fdcoeff(m, xbar, x, fullC=False):\n",
    "    n = len(x) - 1\n",
    "    c1 = 1.\n",
    "    c4 = x[0] - xbar\n",
    "    C = np.zeros((n+1,m+1))\n",
    "    C[0,0] = 1.\n",
    "    for i in range(1,n+1):\n",
    "        mn = min(i,m)\n",
    "        c2 = 1.\n",
    "        c5 = c4\n",
    "        c4 = x[i] - xbar\n",
    "        for j in range(i):\n",
    "            c3 = x[i] - x[j]\n",
    "            c2 = c2*c3\n",
    "            if j==i-1:\n",
    "                for s in range(mn,0,-1):\n",
    "                    C[i,s] = c1*(s*C[i-1,s-1] - c5*C[i-1,s])/c2\n",
    "                C[i,0] = -c1*c5*C[i-1,0]/c2\n",
    "            for s in range(mn,0,-1):\n",
    "                C[j,s] = (c4*C[j,s] - s*C[j,s-1])/c3\n",
    "            C[j,0] = c4*C[j,0]/c3\n",
    "        c1 = c2\n",
    "    \n",
    "    if fullC:\n",
    "        return C\n",
    "    else:\n",
    "        c = C[:,-1] # last column of C\n",
    "        return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.5  2.  -0.5]\n",
      " [-0.5  0.   0.5]\n",
      " [ 0.5 -2.   1.5]]\n"
     ]
    }
   ],
   "source": [
    "# m is the the order of the derative; x is your grid spacing domain \n",
    "# This function uses the the fdceoff to create a matrix where each row is the coefficients for each point in order\n",
    "#For example if x=0,1,2 then it would create a 3x3 matrix where the first row is the coefficients for 0 then next row \n",
    "# are the coefficient for one and so on. this matrix can be multiplied by the range of a function and this should result in\n",
    "# first derative of a range by the domain.\n",
    "\n",
    "def fde(m,x):   \n",
    "    w=np.zeros((len(x),len(x)))\n",
    "    j=0\n",
    "    for i in x:\n",
    "        q=fdcoeff(m, i, x)\n",
    "        w[j]=q\n",
    "        j=j+1\n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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",
    "# first derivative\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef dfdA(np.ndarray arr, double deltaL, double [::1] A, double bd = 0, int exec_pos = -1):\n",
    "    cdef int size = arr.shape[0]\n",
    "    cdef double [:] f = arr\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    \n",
    "    cdef int i\n",
    "    for i in prange(exec_pos+2, size - 1,nogil=True):\n",
    "        res[i] = (f[i+1] - f[i-1]) / (A[i+1] - A[i-1])\n",
    "    if(exec_pos > -1):\n",
    "        res[exec_pos+1] = (f[exec_pos+2] - f[exec_pos+1]) / (A[exec_pos+2]-A[exec_pos+1]) \n",
    "    elif(bd < 1e99): # if bd value is valid, set 0th component to bd\n",
    "        res[0] = bd\n",
    "    else: # invalid bd value, so use one-side derivative\n",
    "        res[0] = (f[1] - f[0]) / (A[1] - A[0])\n",
    "        #res[0] = (-11 * f[0] / 6 + 3 * f[1] - 3 * f[2] / 2 + f[3] / 3) / (A[1] - A[0])\n",
    "    res[size - 1] = (f[-1] - f[-2]) / (A[size - 1] - A[size-2]) \n",
    "    return np.asarray(res)\n",
    "\n",
    "# first derivative\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef stg_dfdA(np.ndarray arr, double [::1] A):\n",
    "    cdef int size = arr.shape[0]\n",
    "    cdef double [:] f = arr\n",
    "    cdef double [:] res =  np.zeros(size-1, dtype=np.double)\n",
    "    \n",
    "    cdef int i\n",
    "    for i in prange(0, size - 1,nogil=True):\n",
    "        res[i] = (f[i+1] - f[i]) / (A[i+1] - A[i])\n",
    "    \n",
    "    return np.asarray(res)\n",
    "\n",
    "# first derivative\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef WENO_dfdA(np.ndarray arr, double deltaL, double [::1] A, double bd=0):\n",
    "    cdef int size = arr.shape[0]\n",
    "    cdef double [:] f = arr\n",
    "    cdef double [:] st_f =  np.zeros(size, dtype=np.double) # i+1/2 field\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    \n",
    "    cdef int i\n",
    "    cdef double b1, b2, b3, w1, w2, w3, ws, st_f1, st_f2, st_f3\n",
    "    for i in prange(2, size - 2,nogil=True):\n",
    "        b1 = (4 * f[i-2]**2 - 19 * f[i-2] * f[i-1] + 25 * f[i-1]**2 + 11 * f[i-2] * f[i] - 31 * f[i-1] * f[i] + 10 * f[i]**2) / 3\n",
    "        b2 = (4 * f[i-1]**2 - 13 * f[i-1] * f[i] + 13 * f[i]**2 + 5 * f[i-1] * f[i+1] - 13 * f[i] * f[i+1] + 4 * f[i+1]**2) / 3\n",
    "        b3 = (10 * f[i]**2 - 31 * f[i] * f[i+1] + 25 * f[i+1]**2 + 11 *f[i] * f[i+2] - 19 * f[i+1] * f[i+2] + 4 * f[i+2]**2) / 3\n",
    "        w1 = 1/10/(1e-6 + b1)**2\n",
    "        w2 = 6/10/(1e-6 + b2)**2\n",
    "        w3 = 3/10/(1e-6 + b3)**2\n",
    "        ws = w1 + w2 + w3\n",
    "        w1 = w1 / ws\n",
    "        w2 = w2 / ws\n",
    "        w3 = w3 / ws\n",
    "        st_f1 = 3 * f[i-2] / 8 - 5 * f[i-1] / 4 + 15 * f[i] / 8\n",
    "        st_f2 = -f[i-1] / 8 + 3 * f[i] / 4 + 3 * f[i+1] / 8\n",
    "        st_f3 = 3 * f[i] / 8 + 3 * f[i+1] /4 - f[i+2] / 8\n",
    "        st_f[i] = st_f1 * w1 + st_f2 * w2 + st_f3 * w3\n",
    "        \n",
    "    for i in prange(3, size - 2,nogil=True):\n",
    "        res[i] = (st_f[i] - st_f[i-1]) / ((A[i+1] - A[i-1]) / 2)\n",
    "    \n",
    "    res[2] = (f[3] - f[1]) / (A[3] - A[1])\n",
    "    res[1] = (f[2] - f[0]) / (A[2] - A[0])\n",
    "    res[0] = (f[1] - f[0]) / (A[1] - A[0])\n",
    "    \n",
    "    i = size-1\n",
    "    res[size - 1] = (f[size-1] - f[size-2]) / (A[size-1] - A[size-2])\n",
    "    i = size-2\n",
    "    res[size-2] = (f[size-1] - f[size-3]) / (A[size-1] - A[size-3])\n",
    "    \n",
    "    return np.asarray(res)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef inv_derv_phi(np.ndarray rho, np.ndarray P, double off_set=0, int exec_pos = -1):\n",
    "    cdef int size = rho.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    cdef int i\n",
    "    res[size - 1] = off_set\n",
    "    res[size - 2] =  (P[size-1] - P[size-2]) / (P[size-2] + rho[size-2]) + res[size-1]\n",
    "    for i in range(size-3, -1, -1):\n",
    "        res[i] = (P[i+2] - P[i]) / (P[i+1] + rho[i+1]) + res[i+2]\n",
    "    return np.asarray(res)\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "# for exp(\\psi) in Eq. 165e\n",
    "cpdef derv_psi( np.ndarray xi, np.ndarray c, double off_set=0):\n",
    "    cdef int size = xi.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    cdef int i\n",
    "    res[0] = off_set\n",
    "    res[1] = - (xi[2] - xi[1]) * c[0] + res[0]\n",
    "\n",
    "    for i in range(2, size):\n",
    "        res[i] = - (xi[i] - xi[i-2]) * c[i - 1] + res[i-2]\n",
    "    return np.asarray(res)\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "# for exp(\\psi) in Eq. 165e\n",
    "cpdef inv_derv_psi( np.ndarray xi, np.ndarray c, double off_set=0):\n",
    "    cdef int size = xi.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    cdef int i\n",
    "    res[size - 1] = off_set\n",
    "    res[size - 2] =  (xi[size-1] - xi[size-2]) * c[size-2] + res[size-1]\n",
    "    #for i in range(size-3, -1, -1):\n",
    "    #    res[i] = (xi[i+1] - xi[i]) * c[i] + res[i+1]\n",
    "    for i in range(size-3, -1, -1):\n",
    "        res[i] = (xi[i+2] - xi[i]) * c[i+1] + res[i+2]\n",
    "    return np.asarray(res)\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef find_exec_pos(np.ndarray arr):\n",
    "    cdef int size = arr.shape[0]\n",
    "    cdef double [:] f = arr\n",
    "    \n",
    "    cdef int i\n",
    "    for i in range(size - 1, 1,-1):\n",
    "        if(arr[i] < 1 and arr[i-1] > 1): # found a horizon\n",
    "            return i-1\n",
    "    return -1\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef ms_rho_term(np.ndarray R_in, double [:] m, double [:] A):\n",
    "    cdef int size = R_in.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    cdef double [:] R = R_in\n",
    "    cdef int i\n",
    "    for i in prange(1, size - 1,nogil=True):\n",
    "        res[i] = A[i] * R[i] * (m[i+1] - m[i-1]) \\\n",
    "                / (3 * (A[i+1] * R[i+1] - A[i-1] * R[i-1]))\n",
    "    i = 0\n",
    "    res[i] = A[i] * R[i] * (m[i+1] - m[i]) \\\n",
    "                / (3 * (A[i+1] * R[i+1] - A[i] * R[i]))\n",
    "    i = size - 1\n",
    "    res[i] = A[i] * R[i] * (m[i] - m[i-1]) \\\n",
    "                / (3 * (A[i] * R[i] - A[i-1] * R[i-1]))\n",
    "    return np.asarray(res)\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef ms_rho_term_stg(np.ndarray R_in, double [:] m, double [:] A,\n",
    "                     double [:] R_stg, double [:] A_stg):\n",
    "    cdef int size = R_in.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size-1, dtype=np.double)\n",
    "    cdef double [:] R = R_in\n",
    "    cdef int i\n",
    "    for i in prange(0, size - 1,nogil=True):\n",
    "        res[i] = A_stg[i] * R_stg[i] * (m[i+1] - m[i]) \\\n",
    "                / (3 * (A[i+1] * R[i+1] - A[i] * R[i]))\n",
    "    return np.asarray(res)\n",
    "    \n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef hm_rho_term(np.ndarray R_in, double [:] m, double [:] A, double [:] xi,\n",
    "               double a):\n",
    "    cdef int size = R_in.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    cdef double [:] R = R_in\n",
    "    cdef int i\n",
    "    for i in prange(1, size - 1,nogil=True):\n",
    "        res[i] = (m[i+1] - m[i-1] - 2 * m[i] * (xi[i+1] - xi[i-1])) \\\n",
    "                / (a * A[i] * R[i] * (xi[i+1] - xi[i-1]) + (A[i+1] * R[i+1] - A[i-1] * R[i-1]))\n",
    "        \n",
    "    i = 0\n",
    "    res[0] = (m[i+1] - m[i] - 2 * m[i] * (xi[i+1] - xi[i])) \\\n",
    "                / (a * A[i] * R[i] * (xi[i+1] - xi[i]) + (A[i+1] * R[i+1] - A[i] * R[i]))\n",
    "    i = size - 1\n",
    "    res[size-1] = (m[size-1] - m[i-1] - 2 * m[i] * (xi[size-1] - xi[i-1])) \\\n",
    "                / (a * A[i] * R[i] * (xi[size-1] - xi[i-1]) + (A[size-1] * R[size-1] - A[i-1] * R[i-1]))\n",
    "    return np.asarray(res)\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef rho_term_stg(np.ndarray R_in, double [:] m, double [:] A, double [:] xi, \n",
    "               double [:] R_stg, double [:] m_stg, double [:] A_stg, double a):\n",
    "    cdef int size = R_in.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size-1, dtype=np.double)\n",
    "    cdef double [:] R = R_in\n",
    "    cdef int i\n",
    "    for i in prange(0, size - 1,nogil=True):\n",
    "        res[i] = (m[i+1] - m[i] - 2 * m_stg[i] * (xi[i+1] - xi[i])) \\\n",
    "                / (a * A_stg[i] * R_stg[i] * (xi[i+1] - xi[i]) + (A[i+1] * R[i+1] - A[i] * R[i]))\n",
    "        \n",
    "    return np.asarray(res)\n",
    "\n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef ephi_term(np.ndarray R_in, double [:] U, double [:] A, double [:] xi, double [:] g,\n",
    "               double a):\n",
    "    cdef int size = R_in.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    cdef double [:] R = R_in\n",
    "    cdef int i\n",
    "    for i in prange(1, size - 1,nogil=True):\n",
    "        res[i] = (a * (xi[i+1] - xi[i-1]) * A[i] * R[i] + (A[i+1] * R[i+1] - A[i-1] * R[i-1])) /\\\n",
    "                (a * (xi[i+1] - xi[i-1]) * (g[i] + A[i] * R[i] * U[i]))\n",
    "    \n",
    "    i = 0\n",
    "    res[0] = (a * (xi[i+1] - xi[i]) * A[i] * R[i] + (A[i+1] * R[i+1] - A[i] * R[i])) /\\\n",
    "                (a * (xi[i+1] - xi[i]) * (g[i] + A[i] * R[i] * U[i]))\n",
    "    i = size - 1\n",
    "    res[size-1] = (a * (xi[size-1] - xi[i-1]) * A[i] * R[i] + (A[size-1] * R[size-1] - A[i-1] * R[i-1])) /\\\n",
    "                (a * (xi[size-1] - xi[i-1]) * (g[i] + A[i] * R[i] * U[i]))\n",
    "    return np.asarray(res)\n",
    "    \n",
    "@cython.boundscheck(False)  # Deactivate bounds checking\n",
    "cpdef rho_prime(np.ndarray R_in, double [:] m, double [:] A, double [:] xi,\n",
    "               double [:] Rp, double [:] mp, double [:] xip,\n",
    "               double [:] Rpp, double [:] mpp, double [:] xipp, double a):\n",
    "    cdef int size = R_in.shape[0]\n",
    "    cdef double [:] res =  np.zeros(size, dtype=np.double)\n",
    "    cdef double [:] R = R_in\n",
    "    cdef int i\n",
    "    for i in prange(0, size-1 ,nogil=True):\n",
    "        res[i] = ((R[i] + A[i] * Rp[i] + A[i] * a * R[i] * xip[i]  ) \n",
    "                  * (-2 * mp[i] * xip[i] + mpp[i] - 2 * m[i] * xipp[i]) \n",
    "                 - (mp[i] - 2 * m[i] * xip[i] ) * (Rp[i] * (2 + A[i] * a * xip[i]) \n",
    "                        + A[i] * Rpp[i] + a * R[i] * (xip[i] + A[i] * xipp[i])) ) \\\n",
    "                / (R[i] + A[i] * Rp[i] + A[i] * a * R[i] * xip[i])**2\n",
    "    return np.asarray(res)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MS:\n",
    "    \n",
    "    def __init__(self, R, m, U, w, alpha, A, rho0, deltaL, L, trace_ray=False):\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",
    "        self.N = R.shape[0]\n",
    "        self.exec_pos = -1\n",
    "\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",
    "        self.Abar_stg = self.to_stg(self.Abar)\n",
    "        \n",
    "        self.q = 1\n",
    "        self.deltau_i = self.cfl_deltau(self.R, self.m, self.U) * 0.05\n",
    "        \n",
    "        # initialize the poton\n",
    "        self.trace_ray = trace_ray\n",
    "        self.Abar_p = self.Abar[0]\n",
    "        self.U_p = self.U[0]\n",
    "        self.m_p = self.m[0]\n",
    "        self.R_p = self.R[0]\n",
    "        \n",
    "        self.U_hm = np.zeros(self.N-1)\n",
    "        self.m_hm = np.zeros(self.N-1)\n",
    "        self.R_hm = np.zeros(self.N-1)\n",
    "        self.xi_hm = np.zeros(self.N-1)\n",
    "        \n",
    "        self.U_hm[0] = self.U[0]\n",
    "        self.m_hm[0] = self.m[0]\n",
    "        self.R_hm[0] = self.R[0]\n",
    "        \n",
    "        self.xi_hm[0] = 0\n",
    "        self.r_old = 0\n",
    "        return\n",
    "    # convert to half grid\n",
    "    def to_stg(self,arr):\n",
    "        return (arr[0:-1] + arr[1:]) / 2\n",
    "    def to_cubic_stg(self,arr):\n",
    "        a1 = arr[0:-1] ** 3 \n",
    "        a2 = arr[1:] ** 3\n",
    "        return (a1 + a2) / (np.abs(a1 + a2)) * ( np.abs(a1 + a2)/ 2) **(1/3)\n",
    "    def to_idx(self, pos):\n",
    "        return np.searchsorted(self.Abar, pos, \"right\") - 1\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 + ms_rho_term(R, m, self.Abar)\n",
    "        return m + (dfdA(m, self.deltaL, self.Abar,0, self.exec_pos) * self.Abar * R)  \\\n",
    "            / (3* (self.Abar * dfdA(R, self.deltaL, self.Abar, 0, self.exec_pos) \n",
    "                  + R ))\n",
    "    def psi(self, rho, p, Pprime):\n",
    "        #return np.log(rho ** (-3 * self.alpha * self.w / 2))\n",
    "        offset = + np.log(rho[-1]**(-3 * self.alpha * self.w / 2))\n",
    "        return inv_derv_phi(rho, p, offset) \n",
    "    \n",
    "    def Pprime(self, R, m):\n",
    "        R_stg = self.to_stg(R)\n",
    "        m_stg = self.to_stg(m)\n",
    "        rho_stg = m_stg + ms_rho_term_stg(R, m, self.Abar, R_stg, self.Abar_stg)\n",
    "        return self.w * np.concatenate( ([0], stg_dfdA(rho_stg, self.Abar_stg) ,[0] ))\n",
    "        \n",
    "        mprime = dfdA(m, self.deltaL, self.Abar, 0, self.exec_pos)\n",
    "        Rprime = dfdA(R, self.deltaL, self.Abar, 0, self.exec_pos)\n",
    "        mpp = ddfddA(m, self.deltaL, self.Abar, 0, self.exec_pos)\n",
    "        Rpp = ddfddA(R, self.deltaL, self.Abar, 0, self.exec_pos)\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 - self.Abar * Rpp))) \\\n",
    "                / (3 * (R + self.Abar * Rprime)**2)\n",
    "    def k_coeffs(self, R, m, U, Abar_p) :\n",
    "        g = self.gamma(R, m, U)\n",
    "        r = self.rho(R, m)\n",
    "        p = self.P(r)\n",
    "        \n",
    "        Pprime = self.Pprime(R, m)\n",
    "        ep = np.exp(self.psi(r, p, Pprime))\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",
    "        \n",
    "        AR_prime = R + self.Abar * dfdA(R, self.deltaL, self.Abar, 0, self.exec_pos)\n",
    "        \n",
    "        kU = U - self.alpha * ep * \\\n",
    "            (   g**2 * np.concatenate( ([self.w\n",
    "                            * ((m[1] + m[1] - 2 * m[0]) /( self.Abar[1]**2)) *  (1**2 * 5/3)] ,Pprime[1:] / \n",
    "             (self.Abar[1:]) ))  / (R * (AR_prime) * (r + p))  \n",
    "            + (2 * U**2 + m + 3 * p) / 2) \n",
    "\n",
    "        kA_p = self.alpha * interp.griddata(self.Abar, ep * g / AR_prime, \n",
    "                                            Abar_p, method='cubic') \n",
    "        #kA_p.fill(0)\n",
    "        return kR, km, kU, kA_p\n",
    "    \n",
    "\n",
    "    \n",
    "    def run_steps(self,n_steps, exc_intv=0) :\n",
    "        step = 0\n",
    "        \n",
    "        deltau = self.deltau_i\n",
    "        if(self.trace_ray == True):\n",
    "            print(\"Tracing ray is enabled and excision will be performed!\")\n",
    "        else:\n",
    "            print('Not Tracing ray and NO excision will be performed!')\n",
    "        while(step < n_steps) :\n",
    "            #if(step % 100 == 0):\n",
    "                #deltau = self.cfl_deltau(self.R, self.m, self.U) * 0.2\n",
    "            \n",
    "            if(self.to_idx(self.Abar_p) > 50 and self.to_idx(self.Abar_p) < self.N * 0.8):\n",
    "                self.exec_pos = np.max([self.exec_pos, self.to_idx(self.Abar_p) - 10])\n",
    "            exec_arr = np.concatenate(([0] * (self.exec_pos+1),[1] * (self.N - self.exec_pos - 1)))\n",
    "\n",
    "            kR1, km1, kU1, kA_p1 = self.k_coeffs(self.R, self.m, self.U, self.Abar_p)\n",
    "            kR2, km2, kU2, kA_p2 = self.k_coeffs(self.R + deltau/2*kR1, self.m + deltau/2*km1, \n",
    "                                                 self.U + deltau/2*kU1, self.Abar_p + deltau/2*kA_p1)\n",
    "            kR3, km3, kU3, kA_p3 = self.k_coeffs(self.R + deltau/2*kR2, self.m + deltau/2*km2, \n",
    "                                                 self.U + deltau/2*kU2, self.Abar_p + deltau/2*kA_p2)\n",
    "            kR4, km4, kU4, kA_p4 = self.k_coeffs(self.R + deltau*kR3, self.m + deltau*km3, \n",
    "                                                 self.U + deltau*kU3, self.Abar_p + deltau*kA_p3)\n",
    "            \n",
    "            self.R = self.R + (deltau/6*(kR1 + 2*kR2 + 2*kR3 + kR4)) * exec_arr\n",
    "            self.m = self.m + (deltau/6*(km1 + 2*km2 + 2*km3 + km4)) * exec_arr\n",
    "            self.U = self.U + (deltau/6*(kU1 + 2*kU2 + 2*kU3 + kU4)) * exec_arr\n",
    "            \n",
    "            if(self.trace_ray == True):\n",
    "                Abar_p_new = self.Abar_p + deltau/6*(kA_p1 + 2*kA_p2 + 2*kA_p3 + kA_p4)\n",
    "                idx_p_new = self.to_idx(Abar_p_new)\n",
    "                U_p_new = interp.griddata(self.Abar, self.U, Abar_p_new, method='cubic')\n",
    "                m_p_new = interp.griddata(self.Abar, self.m, Abar_p_new, method='cubic')\n",
    "                R_p_new = interp.griddata(self.Abar, self.R, Abar_p_new, method='cubic')\n",
    "\n",
    "                diff = idx_p_new - self.to_idx(self.Abar_p)\n",
    "                if (diff > 1): ##move across more than two grid pints!\n",
    "                    print('Warning!' + str(self.Abar_p) + ' ' + str(Abar_p_new))\n",
    "                    #break\n",
    "                if ( diff > 0): # move across one grid point            \n",
    "                    interp_w = (self.Abar[idx_p_new] - self.Abar_p) / (Abar_p_new - self.Abar_p)\n",
    "                    # linear interpolation\n",
    "                    self.U_hm[idx_p_new] = U_p_new * interp_w + self.U_p * (1 - interp_w)\n",
    "                    self.m_hm[idx_p_new] = m_p_new * interp_w + self.m_p * (1 - interp_w)\n",
    "                    self.R_hm[idx_p_new] = R_p_new * interp_w + self.R_p * (1 - interp_w)\n",
    "                    self.xi_hm[idx_p_new] = self.xi + deltau * interp_w\n",
    "\n",
    "                if(self.xi >= self.xi_hm[0]): #only start advancing photon when the system time is large enough\n",
    "                    self.Abar_p = Abar_p_new\n",
    "                    self.U_p = U_p_new\n",
    "                    self.m_p = m_p_new\n",
    "                    self.R_p = R_p_new\n",
    "                \n",
    "                if(idx_p_new == self.N-2): #going out of the boundary\n",
    "                    break\n",
    "                    \n",
    "            r = self.rho(self.R, self.m)\n",
    "            #if(step %100 == 0):\n",
    "            #    print(deltau, (r - self.r_old).max())  \n",
    "            self.r_old = r\n",
    "            \n",
    "            step+=1\n",
    "            self.xi += deltau\n",
    "            \n",
    "            if(step % 10 == 0):\n",
    "                if(find_exec_pos(self.R**2 * self.m * self.Abar**2 * np.exp(2 * (self.alpha-1) * self.xi)) > 0):\n",
    "                    break\n",
    "            \n",
    "        if step < n_steps:\n",
    "            print(\"Photon has gone out of the outter boundary or horizon is identified!\")\n",
    "  \n",
    "    def run_steps_rk2(self,n_steps) :\n",
    "        step = 0\n",
    "        \n",
    "        #deltau = self.cfl_deltau(R, m, U) * 0.9\n",
    "        deltau = self.deltau_i\n",
    "        while(step < n_steps) :\n",
    "            deltau = self.cfl_deltau(R, m, U) * 0.9\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",
    "\n",
    "\n",
    "            self.R = self.R + deltau * kR2 \n",
    "            self.m = self.m + deltau * km2 \n",
    "            self.U = self.U + deltau * kU2\n",
    "            \n",
    "            step+=1\n",
    "            self.xi += deltau\n",
    "            \n",
    "    def adap_run_steps(self,n_steps, adjust_steps=100, tol=1e-7) :\n",
    "        step = 0\n",
    "        deltau = self.deltau_i\n",
    "        #q = 1\n",
    "        while(step < n_steps) :\n",
    "            #if (step % adjust_steps == 0): # reset_stepping size\n",
    "            #    deltau = self.cfl_deltau(R, m, U) * 0.9\n",
    "            #    q = 1\n",
    "            \n",
    "            kR1, km1, kU1 = self.k_coeffs(self.R, self.m, self.U)\n",
    "            \n",
    "            kR2, km2, kU2 = self.k_coeffs(self.R + deltau/2*kR1, self.m + deltau/2*km1, self.U + deltau/2*kU1)\n",
    "            \n",
    "            kR3, km3, kU3 = self.k_coeffs(self.R + 3*deltau/4*kR2, self.m + 3*deltau/4*km2, self.U + 3*deltau/4*kU2)\n",
    "\n",
    "            \n",
    "\n",
    "            R_new = self.R + deltau/9*(2*kR1 + 3*kR2 + 4*kR3 )\n",
    "            m_new = self.m + deltau/9*(2*km1 + 3*km2 + 4*km3 )\n",
    "            U_new = self.U + deltau/9*(2*kU1 + 3*kU2 + 4*kU3 )\n",
    "            \n",
    "            kR4, km4, kU4 = self.k_coeffs(R_new , m_new , U_new )\n",
    "            \n",
    "            E_R = np.max( np.abs((deltau * (-5*kR1/72 + kR2/12 + kR3/9 - kR4/8))))\n",
    "            E_m = np.max( np.abs(deltau * (-5*km1/72 + km2/12 + km3/9 - km4/8)))\n",
    "            E_U = np.max( np.abs(deltau * (-5*kU1/72 + kU2/12 + kU3/9 - kU4/8)))\n",
    "    \n",
    "            max_err_R = np.max(np.abs(self.R)) * tol\n",
    "            max_err_m = np.max(np.abs(self.m)) * tol\n",
    "            max_err_U = np.max(np.abs(self.U)) * tol\n",
    "            \n",
    "            if(E_R < max_err_R and E_m < max_err_m and E_U < max_err_U):\n",
    "                self.R = R_new\n",
    "                self.m = m_new\n",
    "                self.U = U_new\n",
    "                \n",
    "                kR1 = kR4\n",
    "                km1 = km4\n",
    "                kU1 = kU4\n",
    "                \n",
    "                step+=1\n",
    "                self.xi += deltau\n",
    "                         \n",
    "            # Adjust step size.\n",
    "            self.q = 0.8*np.min((max_err_R/E_R, max_err_m/E_m, max_err_U/E_U) )**(1/3)   # conservative optimal step factor\n",
    "            self.q = min(self.q,10)               # limit stepsize growth\n",
    "            deltau *= self.q\n",
    "            #print(q, E_R, max_err_R)\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",
    "        Pprime = self.Pprime(R, m)\n",
    "                           \n",
    "        ep = np.exp(self.psi(r, p, Pprime))\n",
    "        el =  (dfdA(a * self.A * self.R, self.deltaL, self.Abar, 1) / self.RH) /(a * H * self.RH * g) \n",
    "\n",
    "        return (np.log(1 + el / ep / np.exp(self.xi)  \n",
    "                       * self.alpha * np.concatenate( ([1e10],(self.Abar[1:] - self.Abar[0:-1])) ) / np.sqrt(self.w))).min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HM:\n",
    "    \n",
    "    def __init__(self, MS):\n",
    "        self.R = MS.R_hm\n",
    "        self.m = MS.m_hm\n",
    "        self.U = MS.U_hm\n",
    "        self.xi = MS.xi_hm\n",
    "        self.N = MS.R_hm.shape[0]\n",
    "        \n",
    "        self.w = MS.w\n",
    "        self.deltaL = MS.deltaL\n",
    "        self.L = MS.L\n",
    "        self.A = MS.A[:self.N]\n",
    "        self.alpha = MS.alpha\n",
    "\n",
    "\n",
    "        self.t0 = MS.t0\n",
    "        self.t = t0\n",
    "        self.u = 0\n",
    "\n",
    "        self.RH = self.t0 / self.alpha\n",
    "        self.Abar = MS.Abar[:self.N]\n",
    "        \n",
    "        self.kappa = 2\n",
    "        self.Q = np.zeros(self.N)\n",
    "        self.Q_du = np.zeros(self.N)\n",
    "        self.Q_old = np.zeros(self.N)\n",
    "        self.Qprime = np.zeros(self.N)\n",
    "        \n",
    "        self.deltau_i = self.cfl_deltau(self.R, self.m, self.U, self.xi) * 0.02\n",
    "        \n",
    "\n",
    "        return\n",
    "    # convert to half grid\n",
    "    def to_stg(self,arr):\n",
    "        return (arr[0:-1] + arr[1:]) / 2\n",
    "    def to_cubic_stg(self,arr):\n",
    "        a1 = arr[0:-1] ** 3 \n",
    "        a2 = arr[1:] ** 3\n",
    "        return (a1 + a2) / (np.abs(a1 + a2)) * ( np.abs(a1 + a2)/ 2) **(1/3)\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, U, xi, g, xiprime, Rprime, mprime):\n",
    "        temp = (g + self.Abar * R * U) / (g - (self.w + self.Q ) * self.Abar * R * U ) \\\n",
    "            * (m + self.Abar * R * hm_rho_term(R, m, self.Abar, xi, self.alpha) / 3)\n",
    "        return temp\n",
    "\n",
    "    def rho_stg(self, R, m, U, xi, g, xiprime, Rprime, mprime):\n",
    "        R_stg = self.to_stg(R)\n",
    "        m_stg = self.to_stg(m)\n",
    "        U_stg = self.to_stg(U)\n",
    "        g_stg = self.to_stg(g)\n",
    "        A_stg = self.to_stg(self.Abar)\n",
    "        \n",
    "        if self.Q.shape[0] > 1:\n",
    "            Q_stg = self.to_stg(self.Q)\n",
    "        else:\n",
    "            Q_stg = self.Q\n",
    "        \n",
    "        temp = (g_stg + A_stg * R_stg * U_stg) / (g_stg - (self.w + Q_stg ) * A_stg * R_stg * U_stg ) \\\n",
    "            * (m_stg + A_stg * R_stg * rho_term_stg(R, m, self.Abar, xi, R_stg, m_stg, A_stg, self.alpha) / 3)\n",
    "        return temp\n",
    "        \n",
    "    def ephi(self, R, U, g, xi, xiprime, Rprime):\n",
    "        return ephi_term(R, U, self.Abar, xi, g, self.alpha)\n",
    "    \n",
    "        return (self.alpha * xiprime * self.Abar * R + self.Abar * Rprime + R) \\\n",
    "            / (self.alpha * xiprime* (g + self.Abar * R * U))\n",
    "    \n",
    "    def elambda(self, ephi, exi, xiprime):\n",
    "        return self.alpha * ephi * exi * xiprime\n",
    "    \n",
    "    def epsi(self, R, U, g, xi, rho, ephi):\n",
    "        c = self.alpha - 1 + ephi * self.Abar * R * rho * (1 + self.w + self.Q) /\\\n",
    "            ((g + self.Abar * R * U) * (1+self.w))\n",
    "        #offset = 0\n",
    "        #temp = derv_psi(xi, c, offset)\n",
    "        offset = np.log(1/ephi[-1] * (g[-1] + self.Abar[-1] * R[-1] * U[-1])) \n",
    "        temp = inv_derv_psi(xi, c, offset)\n",
    "        return  (g + self.Abar * R * U) / np.exp(temp)\n",
    "    \n",
    "    def drho(self, R, m, U, g, xi, Rp, mp, xip):\n",
    "        rho_stg = self.rho_stg(R, m, U, xi, g, xip, Rp, mp)\n",
    "        return np.concatenate( ([0], stg_dfdA(rho_stg, self.to_stg(self.Abar)) ,[0]) )\n",
    "        \n",
    "        \n",
    "    \n",
    "            \n",
    "    def k_coeffs(self, R, m, U, xi) :\n",
    "        xiprime = WENO_dfdA(xi, self.deltaL, self.Abar, 1e100)\n",
    "        Rprime = WENO_dfdA(R, self.deltaL, self.Abar, 1e100)\n",
    "        mprime = WENO_dfdA(m, self.deltaL, self.Abar, 1e100)\n",
    "        Uprime = WENO_dfdA(U, self.deltaL, self.Abar, 1e100)\n",
    "                \n",
    "        g = self.gamma(R, m, U)\n",
    "        r = self.rho(R, m, U, xi, g, xiprime, Rprime, mprime)\n",
    "        p = self.P(r)\n",
    "                \n",
    "        exi = np.exp(xi)\n",
    "        ephi = self.ephi(R, U, g, xi, xiprime, Rprime)\n",
    "        elambda = self.elambda(ephi, exi, xiprime)\n",
    "        epsi = self.epsi(R, U, g, xi, r, ephi)\n",
    "        \n",
    "\n",
    "        drho = self.drho(R, m, U, g, xi, Rprime, mprime, xiprime)\n",
    "        drho[0] = 3 * elambda[0] / exi[0] * ((1 + w) * r[0] / ephi[0] - (r[0] + p[0]) * U[0])\n",
    "        drho[-1] = 0\n",
    "        \n",
    "        kxi = epsi / ephi / np.exp(xi) / self.alpha\n",
    "        \n",
    "        kR = epsi / exi * R * (U - 1/ephi)\n",
    "        \n",
    "        km = 3 * epsi / exi * (1/ephi * m * (1+w) - U * (p +m)) \n",
    "                \n",
    "        kU = - epsi / exi / (1 - w - self.Q) * (\n",
    "            (m + 3 * p) / 2 + U**2 - U / self.alpha / ephi \n",
    "        + (w + self.Q) * exi / elambda * Uprime + g * (w + self.Q) / ( np.concatenate( ([1], self.Abar[1:]) ) * R * (1+w + self.Q)) * (\n",
    "            3 * self.Q * U + exi * (self.Qprime / elambda - self.Q_du / epsi) / (self.w + self.Q)\n",
    "         + 3 * (1 + w) * (U - 1 / ephi) + exi / elambda * drho / r) )\n",
    "        \n",
    "        # boundary conditions \n",
    "        kxi[0] = epsi[0] / elambda[0] * (xi[1] - xi[0]) / ( (self.Abar[1] - self.Abar[0]) )\n",
    "        kR[0] = epsi[0] / elambda[0] * (R[1] - R[0]) / ( (self.Abar[1] - self.Abar[0]) )\n",
    "        km[0] = epsi[0] / elambda[0] * (m[1] - m[0]) / ( (self.Abar[1] - self.Abar[0]) )\n",
    "        kU[0] = epsi[0] / elambda[0] * (U[1] - U[0]) / ( (self.Abar[1] - self.Abar[0]) )\n",
    "        \n",
    "        return kxi, kR, km, kU\n",
    "    \n",
    "\n",
    "    \n",
    "    def run_steps(self,n_steps, exc_intv=0) :\n",
    "        step = 0\n",
    "        \n",
    "        deltau = self.deltau_i\n",
    "        while(step < n_steps) :\n",
    "            #if(step % 1000 == 0):\n",
    "            #    deltau = self.cfl_deltau(self.R, self.m, self.U, self.xi) * 0.02\n",
    "            \n",
    "            #der_U = dfdA(np.exp(self.xi * (self.alpha-1)) * self.Abar * self.R * self.U , self.deltaL, self.Abar, 1e100)\n",
    "            \n",
    "            #self.Q = self.kappa * (self.Abar[1])**2 * der_U**2\n",
    "            #self.Qprime = dfdA(self.Q , self.deltaL, self.Abar, 1e100)\n",
    "            #self.Q_du = (self.Q - self.Q_old) / deltau\n",
    "            \n",
    "            #self.Q[der_U > 0] = 0\n",
    "            #self.Qprime[der_U>0] = 0\n",
    "            #self.Q_du[der_U>0] = 0\n",
    "            \n",
    "            kxi1, kR1, km1, kU1 = self.k_coeffs(self.R, self.m, self.U,  self.xi)\n",
    "            kxi2, kR2, km2, kU2 = self.k_coeffs(self.R + deltau/2*kR1, self.m + deltau/2*km1, \n",
    "                                                 self.U + deltau/2*kU1, self.xi)\n",
    "            kxi3, kR3, km3, kU3 = self.k_coeffs(self.R + deltau/2*kR2, self.m + deltau/2*km2, \n",
    "                                                 self.U + deltau/2*kU2, self.xi)\n",
    "            kxi4, kR4, km4, kU4 = self.k_coeffs(self.R + deltau*kR3, self.m + deltau*km3, \n",
    "                                                 self.U + deltau*kU3,  self.xi)\n",
    "            \n",
    "            self.xi = self.xi + (deltau/6*(kxi1 + 2*kxi2 + 2*kxi3 + kxi4)) \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",
    "            self.Q_old = self.Q\n",
    "            \n",
    "            step+=1\n",
    "            self.u += deltau\n",
    "            if self.u > 1100:\n",
    "                print(\"u is larger than 1100!\" )\n",
    "                #break\n",
    "              \n",
    "    def cfl_deltau(self,R, m, U, xi):\n",
    "        xiprime = dfdA(xi, self.deltaL, self.Abar, 1e100)\n",
    "        Rprime = dfdA(R, self.deltaL, self.Abar, 1e100)\n",
    "        mprime = dfdA(m, self.deltaL, self.Abar, 1e100)\n",
    "        Uprime = dfdA(U, self.deltaL, self.Abar, 1e100)\n",
    "        #print((xiprime ))\n",
    "        g = self.gamma(R, m, U)\n",
    "        r = self.rho(R, m, U, xi, g, xiprime, Rprime, mprime)\n",
    "        p = self.P(r)\n",
    "        exi = np.exp(xi)\n",
    "        ephi = self.ephi(R, U, g, xi, xiprime, Rprime)\n",
    "        elambda = self.elambda(ephi, exi, xiprime)\n",
    "        epsi = self.epsi(R, U, g, xi, r, ephi)\n",
    " \n",
    "        return ((1 / np.sqrt(w) - 1)  * (self.Abar[1] - self.Abar[0]) * elambda / epsi).min()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating a mix of uniform and logarithmic grids\n",
    "def mix_grid(left, right, n):\n",
    "    A = np.exp(np.linspace(left, right, n))\n",
    "    dA = A[1] - A[0]\n",
    "    A = np.concatenate( (np.linspace(0, A[0],int(np.ceil(A[0] / dA)), endpoint=False), A))\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating uniform grids\n",
    "def uni_grid( right, n):\n",
    "    A = np.linspace(0, np.exp(right), n)\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=1500\n",
    "\n",
    "lower = 1\n",
    "upper = 2.2\n",
    "#L = np.linspace(-6, 1.5, n)\n",
    "#A = np.exp(L)\n",
    "#A = np.linspace(0, np.exp(1.5), n)\n",
    "A = mix_grid(lower,upper,n)\n",
    "A = uni_grid(upper, n)\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 = 0\n",
    "\n",
    "L = Abar[-1]\n",
    "deltaA = Abar[1] - Abar[0]\n",
    "\n",
    "amp = 0.3\n",
    "delta0 = amp * np.exp(-Abar**2 / 2 /(1.6)**2)\n",
    "delta0P = amp * delta0 * 2 * (-1 / 2 / (1.6)**2 ) * Abar\n",
    "\n",
    "\n",
    "#delta0 = amp * (1 - Abar**2 / (Abar[-1] * RH)**2) *  np.exp(-3 * Abar**2 / 2 / (Abar[-1] * RH)**2)\n",
    "#delta0P = amp * delta0 * 2 * (-2 / (Abar[-1] * RH)**2 ) * Abar\n",
    "\n",
    "m = 1 + delta0 \n",
    "U = 1 - alpha * delta0 / 2\n",
    "R = 1 - alpha / 2 * (delta0 + w * Abar * delta0P / (1 + 3 * w) )\n",
    "\n",
    "\n",
    "ms = MS(R,m,U,w,alpha,A,rho0, deltaL,L, trace_ray=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "26.122044362383328"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(ms.Abar*ms.R * np.exp((ms.alpha-1) * ms.xi))[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1500,)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Abar.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01742631378411197"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Abar[350]- Abar[349]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.017426313784111625"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Abar[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.718281828459045"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.1, 10000.0)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD7CAYAAACMlyg3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAR8UlEQVR4nO3dbYxc113H8e9/d1vcuJBW9RaBncqUB1MrtKJZkqqCqiUU7KqLoS0PTtWqlckqIFe8bEAIXiCgfQeFlGhIo6gSirFC1XatlIgHBYPqQlrUBxsrkTFSs0klrxvkCgeI7P3zYmZ278zOeu/eOzsz3vv9SCvvfTjnnjnend/ec87ejcxEktRMU+NugCRpfAwBSWowQ0CSGswQkKQGMwQkqcEMAUlqMENAkhrMEJCkBhtpCETE7oj4SkS8e5TXlSQNVisEIuLhiLgUEWf79h+KiKcj4kJE3F849FHgZJ1rSpKGJ+o8NiIi3gb8N/DpzLy9s28aeAZ4J7AEPAUcBb4f2APsAi5n5ql6TZck1TVTp3Bmno6I/X277wQuZOZFgIg4ARwBXgnsBg4C/xMRj2fmSp3rS5LqqRUCG9gLPFvYXgLuyszjABHxIdp3AusCYNeuXTk9Pb26vWfPHmZnZ7d08StXrnDrrbdWaPbw6hhGG5aXl7f82rejHfbF8OqwL4bbhqb3xfLyMpcvXwbgxRdf/L/M3FWpAZlZ6wPYD5wtbP8S8FBh+wPAn5ap65Zbbsm67r333rHXMYw23HHHHbXrsC+G2w77Ynh12BfDbQNwNSu+h2/H6qAl4LbC9j7g+TIFr1+/zsLCAouLi5UvPj8/X7nssOoYRhuGwb5YY1+ssS/W3Ox9sbi4yMLCAsD0ZudupNbEMEBnTuBUrk0Mz9CeGL4beI72xPA9mXlus7p2796dV69erdWenWJubo4vf/nL427GRLAv1tgXa+yLNRHxYmburlK27hLRR4EzwIGIWIqIY5l5DTgOPAGcB06WCQCA3bt3174T2Ck66S7siyL7Yo190XMnUPmn59p3AsM0NzeXJrskbU1EfCUz56qU9bERktRgExUCV65ccThIkkoqDAdVXqPqcJAk3eQcDpIkVTJRIeBwkCSV53CQJMnhIElSNRMVAg4HSVJ5DgdJkhwOkiRVYwhIUoNNVAg4JyBJ5TknIElyTkCSVI0hIEkNZghIUoMZApLUYBMVAq4OkqTyXB0kSXJ1kCSpGkNAkhrMEJCkBjMEJKnBDAFJarCJCgGXiEpSeS4RlSS5RFSSVI0hIEkNZghIUoMZApLUYIaAJDWYISBJDWYISFKDjSwEIuINEfFgRDwWEb8+qutKkjZWKwQi4uGIuBQRZ/v2H4qIpyPiQkTcD5CZ5zPzPuCXgUq/1CBJGq66dwKPAIeKOyJiGngAOAwcBI5GxMHOsZ8H/hn4+5rXlSQNQa0QyMzTwAt9u+8ELmTmxcx8CTgBHOmc//nMfCvw/jrXlSQNx8w21LkXeLawvQTcFRFvB94DfBfw+KCCy8vLzM2tjRQtLCx0H44kSSpotVq0Wq3u5p6q9WxHCMSAfZmZTwJP3qjg7OwsPkBOkjZX/CE5Ii5XrWc7VgctAbcVtvcBz5cp6KOkJam8iXiUdETsB05l5u2d7RngGeBu4DngKeCezDy3WV0+SlqStm5sj5KOiEeBM8CBiFiKiGOZeQ04DjwBnAdOlgkASdLo1ZoTyMyjG+x/nA0mf2+kOxw0Pz/P/Px8naZJ0o63uLjYHT73L4tJUlPtmL8s5sSwJJU3ERPDw+SdgCRt3Y65E5AkjdZEhYDDQZJUnsNBkiSHgyRJ1UxUCDgcJEnlORwkSXI4SJJUjSEgSQ1mCEhSg01UCDgxLEnlOTEsSXJiWJJUjSEgSQ1mCEhSg01UCDgxLEnlOTEsSXJiWJJUjSEgSQ1mCEhSgxkCktRghoAkNZghIEkNNlEh4O8JSFJ5/p6AJMnfE5AkVWMISFKDGQKS1GCGgCQ1mCEgSQ1mCEhSgxkCktRgIwuBiPiFiPiLiPhcRPzsqK4rSdpYrRCIiIcj4lJEnO3bfygino6ICxFxP0BmfjYz7wU+BPxKnetKkoaj7p3AI8Ch4o6ImAYeAA4DB4GjEXGwcMrvdI5LksasVghk5mnghb7ddwIXMvNiZr4EnACORNvHgS9k5r/Vua4kaThmtqHOvcCzhe0l4C7gI8DPALdGxA9l5oP9BZeXl5mbW3v8xcLCQvfhSJKkglarRavV6m7uqVrPdoRADNiXmfkJ4BM3Kjg7O4sPkJOkzRV/SI6Iy1Xr2Y7VQUvAbYXtfcDzZQr6KGlJKm8iHiUdEfuBU5l5e2d7BngGuBt4DngKuCczz21Wl4+SlqStG9ujpCPiUeAMcCAiliLiWGZeA44DTwDngZNlAgC8E5CkrZiIO4Fh8k5AkrbOPyojSapkokLA4SBJKs/hIEnSzhkO8k5AksrzTkCStHPuBCRJo2UISFKDTVQIOCcgSeU5JyBJck5AklSNISBJDTZRIeCcgCSV55yAJMk5AUlSNYaAJDWYISBJDWYISFKDTVQIuDpIkspzdZAkydVBkqRqDAFJajBDQJIazBCQpAYzBCSpwSYqBFwiKknluURUkuQSUUlSNYaAJDWYISBJDWYISFKDGQKS1GCGgCQ1mCEgSQ02shCIiNdHxKci4rFRXVOSdGO1QiAiHo6ISxFxtm//oYh4OiIuRMT9AJl5MTOP1bmeJGm46t4JPAIcKu6IiGngAeAwcBA4GhEHa15HkrQNaoVAZp4GXujbfSdwofOT/0vACeBInetIkrbHdswJ7AWeLWwvAXsj4jUR8SDw4xHxW4MKLi8vMzc3t/rRarW2oXmSdPNrtVqr75XAnqr1zAyvSatiwL7MzG8D992o4OzsLD5ATpI2t7Cw0H2CKBFxuWo923EnsATcVtjeBzxfpqCPkpak8ibiUdIRsR84lZm3d7ZngGeAu4HngKeAezLz3GZ1+ShpSdq6sT1KOiIeBc4AByJiKSKOZeY14DjwBHAeOFkmACRJo1drTiAzj26w/3Hg8a3W1x0Omp+fZ35+vk7TJGnHW1xc7A6f+5fFJKmpdsxfFnNiWJLKm4iJ4WHyTkCStm7H3AlIkkZrokLA4SBJKs/hIEmSw0GSpGomKgQcDpKk8hwOkiQ5HCRJqsYQkKQGMwQkqcEmKgScGJak8pwYliQ5MSxJqsYQkKQGMwQkqcEmKgScGJak8pwYliQ5MSxJqsYQkKQGMwQkqcEMAUlqMENAkhrMEJCkBpuoEPD3BCSpPH9PQJLk7wlIkqoxBCSpwQwBSWowQ0CSGswQkKQGmxl3A4rOf+s73PWHf1fq3CBK1xsbnDpodww4ecPyhf3d9vTu6603+g70Hy8ejtVz1tfbX9/quQPKrF0zeuqOvnp6yvYcb+8rXqO4XWx7se717YvCsbU6up+zQZt66iqWG1DPhtfo1lF4jcXzpiJ66p7qK0tPewrndF93Z/9U9Lcpeq5RrK//GmttKJw7oK5iPd3XNNU9NnXj8+k5NwbXUbj2VN+xQfVOdfqgW9+68p193dfX7bee1z/w+ht802noJioEvmfXy3jHgdduet5WVrUmg08eVEcOOL5R+eLu7qfF5bb9dfWf01Nr37VWr73R/p72rz/Wvc7a58VyhTZmu+7i9Xq2gVyBZKXndRRfQ891O9fs1rNapr9c4Trd/cWXU9zXc+5qOwccW23L+nozYaXQ0EH1r2T3WPvz4nGNx6DA6Ammvn+nNgipGwXN1Oo5ayHes81aiG12/W6oFusZfK2+gJzqPb7unKkBZVbrbH9ex8hCICJ2A58EXgKezMy/7D9n76tfwcfe+8ZRNUkqLXN9uKzk+gDtD5BiuZXVz9s7VvpCbGVlLSRXVuso/NtfJnvrXFeuUG9Pe/qu0a6jr86+cuvrTVZWWF+++DpzrT9WCkG8srq/+FrW9hWvNejclZXs6Yfu9VY652/8Orrletu8MqCd/X17faX9+fWV7Om3dX3UV25Qff3/Xu/pr06Zld6vsZ5/+7726qoVAhHxMPBu4FJm3l7Yfwj4E2AaeCgzPwa8B3gsMxcj4q+AdSEgTaricNLggURptIoh9rKPV6+n7sTwI8Ch4o6ImAYeAA4DB4GjEXEQ2Ac82zntes3rSlKjRQRTU8HMdL238VqlM/M08ELf7juBC5l5MTNfAk4AR4Al2kFQ+7qSpOHYjjfjvaz9xA/tN/+9wGeA90bEnwMDnxC3vLzM3Nzc6ker1dqG5knSza/Vaq2+VwJ7qtazHRPDgwZMMzOvAh++UcHZ2Vl8gJwkbW5hYaH7BFEi4nLVerbjTmAJuK2wvQ94vkxBHyUtSeVNxKOkI2I/cKq7OigiZoBngLuB54CngHsy89xmdfkoaUnaurE9SjoiHgXOAAciYikijmXmNeA48ARwHjhZJgDAOwFJ2oqJuBMYJu8EJGnr/KMykqRKJioEHA6SpPIcDpIk7ZzhIO8EJKk87wQkSTvnTkCSNFqGgCQ12ESFgHMCklSecwKSJOcEJEnVGAKS1GATFQLOCUhSec4JSJKcE5AkVWMISFKDGQKS1GCGgCQ12ESFgKuDJKk8VwdJklwdJEmqxhCQpAYzBCSpwQwBSWowQ0CSGmyiQsAlopJUnktEJUkuEZUkVWMISFKDGQKS1GCGgCQ1mCEgSQ1mCEhSgxkCktRgIwuBiHh9RHwqIh4b1TUlSTdWKgQi4uGIuBQRZ/v2H4qIpyPiQkTcf6M6MvNiZh670TnLy8tlmtMIrVZr3E2YGPbFGvtijX3RY0/VgmXvBB4BDhV3RMQ08ABwGDgIHI2IgxHxYxFxqu/jtWUucvny5S00fbBhPHKibh3DaMMwvsDti+G2w74YXh32xXDbAMxWLVgqBDLzNPBC3+47gQudn/BfAk4ARzLzG5n57r6PS1UbuFU76D+1NvtijX2xxr5YY19s4dlBEbEfOJWZt3e23wccysxf62x/ALgrM49vUP41wB8A7wQeysw/GnDO/wLXC7uWga3eHtwKXNlimWHXMYw27GHrr3072mFfDK8O+2K4bWh6X+xh7Q5gOjN3Vbn4TJVCHTFg34aJkpnfBu67UYVVX4QkqZo6q4OWgNsK2/uA5+s1R5I0SnVC4CnghyPiByLi5cCvAp8fTrMkSaNQdonoo8AZ4EBELEXEscy8BhwHngDOAycz81zJ+m64tDTaPtE5/vWIeHP5l3RzKdEX7+/0wdcj4osR8aZxtHMUyi45joifiIjrnXmpHalMX0TE2yPiqxFxLiL+cdRtHJUS3yO3RsRiRHyt0xcfHkc7t9tGS/ULx6u9b2bmSD+AaeA/gNcDLwe+BhzsO+ddwBdozzu8BfiXUbdzgvrircCrO58fbnJfFM77B+Bx4H3jbvcYvy5eBfw78LrO9mvH3e4x9sVvAx/vfD5LeyXjy8fd9m3oi7cBbwbObnC80vvmOB4bMXBpad85R4BPZ9uXgFdFxPeNuqEjsGlfZOYXM/O/Optfoj33shOV+boA+Ajw18DIlh2PQZm+uAf4TGZ+EyBHuAx7xMr0RQLfHREBvJJ2CFwbbTO3Xw5eql9U6X1zHCGwF3i2sL3U2bfVc3aCrb7OY7STfifatC8iYi/wi8CDI2zXOJT5uvgR4NUR8WREfCUiPjiy1o1Wmb74M+ANtBemfAP4zcxcGU3zJkql9806S0SrKrO0dEvLT29ipV9nRLyDdgj85La2aHzK9MUfAx/NzOvtH/p2rDJ9MQPcAdwNvAI4ExFfysxntrtxI1amL34O+Crw08APAn8bEf+Umd/Z7sZNmErvm+MIgTJLS5uy/LTU64yINwIPAYez/fsWO1GZvpgDTnQCYA/wroi4lpmfHU0TR6bs98jlzLwKXI2I08CbgJ0WAmX64sPAx7I9MH4hIv4T+FHgX0fTxIlR6X1zHMNBZZaWfh74YGe2+y3Alcz81qgbOgKb9kVEvA74DPCBHfhTXtGmfZGZP5CZ+zNzP/AY8Bs7MACg3PfI54CfioiZiLgFuIv2Kr2dpkxffJP2HRER8b3AAeDiSFs5GSq9b478TiAzr0VEd2npNPBwZp6LiPs6xx+kvfLjXcAF4EXaSb/jlOyL3wVeA3yy8xPwtcycG1ebt0vJvmiEMn2Rmecj4m+ArwMrtB/FMnDp4M2s5NfF7wOPRMQ3aA+JfDQz6z+NcsJ0luq/HdgTEUvA7wEvg3rvm6WfHSRJ2nn8y2KS1GCGgCQ1mCEgSQ1mCEhSgxkCktRghoAkNZghIEkNZghIUoP9PzHtN3X5MlMRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot((ms.Abar*ms.R * np.exp((ms.alpha-1) * ms.xi)), ms.rho(ms.R,ms.m))\n",
    "plt.xlim(0,1)\n",
    "plt.yscale('log')\n",
    "plt.ylim(1e-1,1e4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20e8b9745b0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD1CAYAAABJE67gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQiUlEQVR4nO3df5Bd5V3H8feXDZkpPwTabJEJqaQdaifTKVrWoJZWaqc0MINpHbWAlhaFLSNBmVEH9A/tDOOMWquMU2i8YFoYW6JOsU1sLHacCrUtNQQQCBTcSZGEIGyKBoHRGPL1j70L995zs3tvcu/uee6+XzOZ7Dn3uSffZ0/uJyfPPs85kZlIkspzzGIXIEk6Mga4JBXKAJekQhngklQoA1ySCmWAS1Khli3kH7ZixYo844wzANi/fz8nnXRSz++dnp5mfHy85/b9HL/fWobZfpj97Lf9sL8v/fTVfi58ez+j3S3GZ3THjh37MrP6h2bmgv06++yzc9aVV16Z/Wh9by/6OX6/tQyz/TD72W/7YX9f+umr/Vz49n5Gu1uMzyhwX3bJ1EUbQrnoootqc/x+axl2+2Eeu07fl2Ee234efXs/o4MxzFoiF3Al5sTERN53331H+l6O9L0lWSr9hKXTV/s5WhajnxGxIzMnOvcv6BX4/v37mZycZOvWrX2/d3JycggV1c9S6Scsnb7az9GykP3cunXr7J/XdRC9mCtwSVqqanEFLkkaHANckgpVRIAfOpQcOuRtbyWpVREB/oEb72HDHfcvdhmSVCtFBLgkqcoAl6RCFRPgPvlNktoVEeARi12BJNVPEQEuSaoqZim9JC01I7GU/vw/vZs3rziBjR85ewhVSVK9Fb2UPnAQXJI6FRHgkqQqA1ySClVMgCdOBJekVkUEuPPAJamqiACXJFUZ4JJUqGIC3HuhSFK7YgJcktTOpfSSVFMjsZR+3Y33sOr1x3HLZZWVpJI08opeSi9JqioiwMOJ4JJUUUSAS5KqiglwpxFKUrsiAtwBFEmq6inAI2JdRDweEVMRcX2X10+KiK0R8a8RsTMiLh98qZKkVvMGeESMATcBFwBrgEsiYk1Hs6uBRzPzLOA84FMRsXzAtUqSWvRyBb4WmMrMXZl5ANgMrO9ok8CJMTNd5ATgeeDgQCv1drKS1GZZD21WArtbtvcA53S0+TSwBdgLnAh8ODMPdR5oenqaiYnX5qJPTk7OrjKak7MIJS01jUaDRqMxu7miW5teArxbfHZeDn8AeBD4aeAtwNci4huZ+UJro/HxcY5kJaYkLTWtF7gRsa9bm16GUPYAq1q2T2fmSrvV5cCdOWMK+B7wtr4rliT1rJcA3w6cGRGrmz+YvJiZ4ZJWTwHvA4iIU4EfBnYNslDngUtSu3mHUDLzYERsAO4CxoBNmbkzIq5qvr4RuAH4XEQ8zMyQy3WZ2fWS/0g4Bi5JVb2MgZOZ24BtHfs2tny9Fzh/sKVJkuZSxEpMSVJVMQHuELgktSsiwMO7oUhSRREBLkmqMsAlqVDFPNR4IZ/dKUl1MBIPNf6ZT/8zbzh+OZ+9fO0QqpKkevOhxpI0YooJcAdQJKldEQHuJEJJqioiwCVJVQa4JBWqmAB3FqEktSsjwL2frCRVlBHgkqQKA1ySClXOUvoh1CNJdTYSS+k/eNM3+YHXHcvtv+xSeklLj0vpJWnEGOCSVKhiAtzbyUpSuyIC3GngklRVRIBLkqoMcEkqlAEuSYUqIsAdApekqiICXJJUVc5SemcRSlpiRmIp/c/e/E2OW76Mv7zinCFUJUn1VvRS+nAiuCRVFBHgkqSqYgI8vaGsJLUpIsAdQJGkqiICXJJU1VOAR8S6iHg8IqYi4vrDtDkvIh6MiJ0Rcfdgy5QkdVo2X4OIGANuAt4P7AG2R8SWzHy0pc3JwM3Ausx8KiLeOOhCnQcuSe16uQJfC0xl5q7MPABsBtZ3tLkUuDMznwLIzOcGWaSzCCWpqpcAXwnsbtne09zX6q3AKRHxTxGxIyIu63ag6elpJiYmXv3VaDSOrGpJGnGNRuPVrARWdGsz7xAK3SeBdA5oLAPOBt4HvA74dkTcm5lPtDYaHx/nSFZiStJSMzk5ObuMnojY161NLwG+B1jVsn06sLdLm32Z+RLwUkTcA5wFPMGAOAYuSe16GULZDpwZEasjYjlwMbClo82XgXdHxLKIOA44B3hsUEWGM8ElqWLeK/DMPBgRG4C7gDFgU2bujIirmq9vzMzHIuKrwEPAIeDWzHxkmIVL0lLXyxAKmbkN2Naxb2PH9ieBTw6uNEnSXIpZiem9UCSpXRkB7hC4JFWUEeCSpIpiAtxphJLUrpgAlyS1K+Khxg6BS1qKRuKhxh/+828D8Fcf/4lBlyRJtVf0Q42hevMVSVrqighwbycrSVVFBLgkqcoAl6RClRPgDoJLUpsiAtzbyUpSVREBLkmqMsAlqVDFBLi3k5WkdmUspXcIXNISNBJL6S+95V7+75VD/M1VPzmEqiSp3opfSi9JaldMgHs/cElqV0SAOwYuSVVFBLgkqaqYAHcERZLaFRPgkqR2RQS490KRpKoiAlySVFXESkyAhVxwJEl1MBIrMX/p1u/w8oGD3Pmr7xpCVZJUb0WvxHQeuCRVFRHgkqSqYgLcEXBJaldMgEuS2hngklSongI8ItZFxOMRMRUR18/R7sci4pWI+LnBlShJ6mbeAI+IMeAm4AJgDXBJRKw5TLs/BO4adJHg7WQlqVMvV+BrganM3JWZB4DNwPou7a4Bvgg8N8D6AAjnEUpSRS8BvhLY3bK9p7nvVRGxEvgQsHFwpUmS5rKshzbdLn87BzRuBK7LzFfmulqenp5mYuK1xUSTk5Ozy0QlSS0ajQaNRmN2c0W3Nr0E+B5gVcv26cDejjYTwOZmeK8ALoyIg5n5pdZG4+PjHMlSenAeuKSlpfUCNyL2dWvTS4BvB86MiNXA08DFwKWtDTJz9ezXEfE54O86w/toOAIuSVXzBnhmHoyIDczMLhkDNmXmzoi4qvm6496StAh6uQInM7cB2zr2dQ3uzPzY0ZfV9cBDOawklcqVmJJUqCIC3GngklRVRIBLkqqKCXBHwCWpXTEBLklqV0SAOwQuSVXFPJVekpaakXgq/eWf/Rf2vXiArdecO4SqJKneCn8qvYMoktSpiACXJFUZ4JJUqGICPJ0JLkltighwR8AlqaqIAJckVRngklSoYgLc24FLUrsiVmI6DVzSUjQSKzGvuG07z+z/H77ya+8eQlWSVG9Fr8QEh1AkqVMxAS5JaldIgDsILkmdCglwSVKnYgLcIXBJaldMgEuS2hUR4M4Dl6SqIgJcklRVTIAv5IIjSSpBGUvph1SPJNXZSCyln7z9Pp56/mW+eu17hlCVJNVb8UvpJUntDHBJKlQRAe40QkmqKiLAJUlVBrgkFaqnAI+IdRHxeERMRcT1XV7/xYh4qPnrWxFx1qALdRq4JLWbN8AjYgy4CbgAWANcEhFrOpp9D/ipzHwHcAPQGGSR4UxwSaro5Qp8LTCVmbsy8wCwGVjf2iAzv5WZ/9ncvBc4fbBlSpI69RLgK4HdLdt7mvsO51eAvz+aorpJbygrSW2W9dCm2/hF1zSNiPcyE+Dndnt9enqaiYnXFhNNTk7OLhOVJLVoNBo0Gq+ORq/o1qaXAN8DrGrZPh3Y29koIt4B3ApckJnf73ag8fFxjmQpvfPAJS01rRe4EbGvW5tehlC2A2dGxOqIWA5cDGxpbRARbwLuBD6SmU8cVdWSpJ7MewWemQcjYgNwFzAGbMrMnRFxVfP1jcDvAm8Abo6Zy+WD3W68cjScRihJ7XoZQiEztwHbOvZtbPn6CuCKwZYmSZpLESsxHQOXpKoiAlySVFVMgDsELkntiglwSVK7IgLce6FIUlURDzWWpKVoJB5qfPXn7+e7//EC//gb5w2+KEmqubIfauwIiiRVlBHgkqQKA1ySClVMgDsPXJLaFRHgDoFLUlURAS5JqionwB1DkaQ25QS4JKlNEQEe3k9WkipcSi9JNTUSS+mvueMBHnl6P1//zfMGX5Qk1VzZS+klSRVFBLgj4JJUVUSAS5KqignwhRyrl6QSFBPgkqR2RQS408AlqaqIAJckVRUT4I6AS1K7IgLcERRJqnIpvSTV1Egspb928wM8sPu/uPu33juEqiSp3opfSu80cElqV0SAeztZSaoqIsAlSVXFBHg6kVCS2hQT4JKkdj0FeESsi4jHI2IqIq7v8npExJ81X38oIt45yCIdAZekqnkDPCLGgJuAC4A1wCURsaaj2QXAmc1fk8BnBlynJKnDsh7arAWmMnMXQERsBtYDj7a0WQ/cnjOTyu+NiJMj4rTMfGZQhb78v6/wlYcGdjhJWlDr3v6DjB0z2PGEXgJ8JbC7ZXsPcE4PbVYCbYk7PT3NxMRrc9EnJydnVxnN6eTjlvP9lw5w9Rfu76FcSaqf796wjrFjxnpu32g0aDQas5srurXpJcC7/ZPROSWklzaMj49zJCsxf/vCt3Hx2lV9v0+S6mL5WH9zRlovcCNiX7c2vQT4HqA1PU8H9h5BmyN27NgxvPXUEwd1OEkaCb38k7AdODMiVkfEcuBiYEtHmy3AZc3ZKD8O7B/k+LckqWreK/DMPBgRG4C7gDFgU2bujIirmq9vBLYBFwJTwMvA5cMrWZIEvQ2hkJnbmAnp1n0bW75O4OrBliZJmosrMSWpUMUEeMt0mpG2VPoJS6ev9nO01Kmfixbg/T6Vp99vWj/H77eWYbYfZj/7bT/s70s/fbWfC9/ez2h3dfqMFhPgwzx+nf5y9KtOH2D7efTq9HfRz+hgDLOWBX2kWkRMA//e3DwJ2N/H21cAXSezH0Y/x++3lmG2H2Y/+20/7O9LP321nwvf3s9od4vxGf2hzBzvbLigAS5JGpxifogpSWpngEtSoWof4PM9TGKURMSTEfFwRDwYEf3f9aumImJTRDwXEY+07Ht9RHwtIv6t+fspi1njoBymr5+IiKeb5/XBiLhwMWs8WhGxKiK+HhGPRcTOiPj15v6RO6dz9LUW57TWY+DNh0k8AbyfmRtmbQcuycxH53xjoSLiSWAiM/v5AUntRcR7gBeZuWf825v7/gh4PjP/oPkP8ymZed1i1jkIh+nrJ4AXM/OPF7O2QYmI04DTMvP+iDgR2AF8EPgYI3ZO5+jrL1CDc1r3K/BXHyaRmQeA2YdJqCCZeQ/wfMfu9cBtza9vY+ZDUbzD9HWkZOYzmXl/8+v/Bh5j5v7/I3dO5+hrLdQ9wA/3oIhRlcA/RMSOiJj/SRdlO3X2jpXN39+4yPUM24bm82I3jcLQwqyIOAP4UeA7jPg57egr1OCc1j3Ae3pQxAh5V2a+k5lnjF7d/O+4yvcZ4C3AjzDzlKpPLW45gxERJwBfBK7NzBcWu55h6tLXWpzTugf4UB8UUTeZubf5+3PA3zIzhDSqnm2OL86OMz63yPUMTWY+m5mvZOYh4BZG4LxGxLHMBNrnM/PO5u6RPKfd+lqXc1r3AO/lYRIjISKOb/6QhIg4HjgfeGTudxVtC/DR5tcfBb68iLUM1WyoNX2Iws9rRATwF8BjmfknLS+N3Dk9XF/rck5rPQsFoDk950Zee5jE7y9ySUMREW9m5qobZu7T/oVR6WtE3AGcx8wS5GeB3wO+BPw18CbgKeDnM7P4H/4dpq/nMfNf7QSeBD5e8hOrIuJc4BvAw8Ch5u7fYWZseKTO6Rx9vYQanNPaB7gkqbu6D6FIkg7DAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVD/DxrslTCODd1SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ms.Abar[0:ms.N-1], ms.U_hm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20e8ba12340>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD1CAYAAACrz7WZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfZRV9X3v8feHgQF5EBBGVEBARRGpGp0QmwfzoImYdmrjbbqkN01qNFNXgze9XW2u16a1z/X2aYXeZUtPEnrjSqorSUOKKQ0xNpHcGxUfAHnWAQSGpxlAYAaEYZjv/ePsIcfJDHOGOWfOPud8XmudNWfv/dt7f3/nMN/Z/H6//duKCMzMrLIMK3UAZmZWeE7uZmYVyMndzKwCObmbmVUgJ3czswrk5G5mVoGGlzoAgMmTJ8fMmTPPLh89epTx48fntW9rayt1dXV5n2sgxx5o+WIeGwZW12LHUszPxfUc+liKWd6/o8WL5eWXXz4YEb2fMCJK/rr55psj12c+85nIV899+zOQYw+0fDGPHTGwuhY7lmJ+Lq7n0MdSzPL+HS1eLMBL0UdeTWWzTENDQ2qOPZDyxTz2QBU7lmJ+LsWKo9jli1nPgR4/TZ/LQKXpOyrXegIoUnCHan19fbz00kvnuy/nu2+5qZa6up6VpVrqCUNfV0kvR0R9b9tSceV+9OhRGhsbeeqppwa8b2NjYxEiSqdqqavrWVmqpZ4wdHV96qmnus/VZ6N92V+5m5lVq0FduUtaKqlF0oY+tkvS30tqkvSqpJtyti2QtDXZ9tD5V8HMzAYin2aZ/wMsOMf2O4HZyasR+EcASTXAY8n2ucBCSXMHE6yZmeWn33HuEbFK0sxzFLkLeDwZlvO8pAmSLgVmAk0RsR1A0pNJ2U2DDdqKp5TNdCloITQbchJIKvhxC3ET01Rgd85yc7Kut/XvKsD5qlpE0NJ2ip2HTrDz0HFa2k5x+HgHb57o4M3jHRzvOMOpzi5Onc7+PJn87IqgqyuIIPs++fnT5ew6Mxtaf/bL8/jELTMKftxCJPfe/uTEOdb/jNbWVurrf9on0NjYWFU97Odypit4Zdeb/Pi1VtbsPsK63Uc4drLzbWVG19YwcXQtE0aPYOzI4Yy/YASjxo1k5IgaRg0fRu3wYQwfJiQxTGKYYNgwIfHTZWW3i+yVRKmo1382ZpXr+mn536UKkMlkyGQy3YuT+ypXiOTeDEzPWZ4G7AVq+1j/M+rq6qpmHGy+tu5v42vP7+Tf1+/j8PEOhgmuueRCfuH6S5lzyYXMmDSaGZPGcOn4UYwaUVPqcM1siORe/Eo62Fe5QiT35cCipE39XcDRiNgnqRWYLWkWsAe4B/i1Apyvom3ce5S/+t5Wnn2tldrhw7jjuku447op3Hp1HReOGlHq8MysTPSb3CU9AXwAmCypGXgEGAEQEUuAFcBHgSbgBHBvsq1T0iJgJVADLI2IjUWoQ0U4fqqTv/yPzXzt+V2Mv2AEv3fHNSycfzkXjaktdWhmVobyGS2zsJ/tAXy2j20ryCZ/O4eNe4/yW19/hV2HT/Dp98zic7fPZvwFvko3s/NX9tMPlLunNx3g40ueo6Oziyc/cwt/2DDXid3MzsnTD6Tcf6zfx6In1nDdZRfy5U/Wc/GFo0odkpmVkXNNP5CKh3VUox9uaeHBJ9Zw4/QJPP7p+YwZ6a/CzAonFc0y1aappY0Hn1jDnEvH8c/3vtOJ3cwKzsl9iLWdPM1nHn+ZUSOGkfn1eg9vNLOiSEVyr6YO1T/97iZ2HjrOP/zXm7lswgWlDsfMypA7VFPmmc0HuO+rL/FbH7iSzy+YU+pwzKzMpf5JTNXgREcnX/jOBuZcMo7P3T671OGYWYVzT94QWfLsdvYdPcn/XvgORg73XDBmVly+ch8Ce468xT89u42GGy6jfuZFpQ7HzKqAk/sQ+PsfvE4AD93pdnYzGxqpSO6VPFpm9+ET/OsrzSx853SmenSMmRWAR8ukwO8vW883X2rm2c9/gEvHO7mbWeF4tEyJtLSd5JsvNfMr9dOc2M1sSDm5F9ETL+ym40wX9793VqlDMbMq4+ReJB2dXXz9hZ28/+o6rqgbW+pwzKzKpCK5V2KH6sqN+2lpO8Wn3l34p5qbWXUrWIeqpAXAYrKPy/tyRDzaY/tEYClwJXAS+HREbEi2fQ74DCDgSxHxxZ7Hr8QO1YWZ52k+coIf/e4HqRmmUodjZhVoUB2qkmqAx4A7gbnAQklzexR7GFgbEdcDnyT7hwBJ88gm9vnADcAvSqr4e+/3HHmL57Yf4ldumu7EbmYlkU+zzHygKSK2R0QH8CRwV48yc4FnACJiCzBT0hTgWuD5iDgREZ3As8DHChZ9Sn1nzR4APvaOqSWOxMyqVT7JfSqwO2e5OVmXax1wN4Ck+cAMYBqwAbhV0iRJo4GPAtMHG3SaRQTffqWZd86cyOWTRpc6HDOrUvlMHNZbu0LPhvpHgcWS1gLrgTVAZ0RslvS/gKeBdrJ/BDp7Hqy1tZX6+p82GzU2NnZ3FpSd9XuOsq31OPe/74pSh2JmFSiTyZDJZLoXJ/dVLp/k3szbr7anAXtzC0TEMeBeAEkCdiQvIuIrwFeSbX+RHO9t6urqqJQO1eVr91JbM4yP/tylpQ7FzCpQ7sWvpIN9lcunWeZFYLakWZJqgXuA5bkFJE1ItgHcD6xKEj6SLk5+Xk626eaJAdalbEQE39u4n/fOnsz4C/z4PDMrnX6v3COiU9IiYCXZoZBLI2KjpAeS7UvIdpw+LukMsAm4L+cQ/yppEnAa+GxEvFnoSqTFpn3HaH7zLR780FWlDsXMqlxeD+uIiBXAih7rluS8fw7odYhjRLxvMAGWk5UbDzBMcNu1U0odiplVOd+hWkDf37if+hkXMXnsyFKHYmYVzFP+DqGdh47z/r/+EV/4hWs9UsbMhoSn/B0Cz2xuAeAjcy8pcSRmZk7uBbPq9VZmTR7jG5fMLBWc3Avg5OkzPL/9ELfO7vN+AjOzIeXkXgAvvfEmJ093cevVdaUOxcwMSElyL/fRMqteb2VEjbjlikmlDsXMqoBHywyRBV9cxcTRtTzReEupQzGzKuLRMkXU0naSLfvbeN/Vbm83s/Rwch+k57cfBuA9Vzq5m1l6OLkP0uodhxg7cjjXXXZhqUMxMzsrFcm9nDtUV+84zM0zJjK8JhUfpZlVAXeoFtnh4x3c9KdP83t3XMNnP+iZIM1saLlDtUhW78i2t99yxUUljsTM7O2c3AfhhR2HGDViGD83dUKpQzEzexsn90FYveMwN10+kdrh/hjNLF1SkZXKsUP12MnTbNp3jPmz3CRjZkPLHapF9OxrrXxq6Wq+fv+7eM9VHuNuZkNv0B2qkhZI2iqpSdJDvWyfKGmZpFclrZY0L2fbf5e0UdIGSU9IGnX+VUmPNbveRIIbpru93czSp9/kLqkGeAy4E5gLLJQ0t0exh4G1EXE98ElgcbLvVOC/AfURMY/sA7bvKVz4pbNm1xGuvngcY0fm9RhaM7Mhlc+V+3ygKSK2R0QH8CRwV48yc4FnACJiCzBTUvdToocDF0gaDowG9hYk8hKKCNbuPsI7LvdVu5mlUz7JfSqwO2e5OVmXax1wN4Ck+cAMYFpE7AH+BtgF7AOORsT3e56gtbWV+vr6s69MJjPwmgyhHQePc/St007uZjbkMpnM2VwJ9Nnhl0+bgnpZ17MX9lFgsaS1wHpgDdApaSLZq/xZwBHgm5I+ERFfy925rq6OcupQXbv7CAA3Tp9Y4kjMrNo0NjZ2j5RB0sG+yuWT3JuB6TnL0+jRtBIRx4B7k5MJ2JG87gB2RERrsu3bwLuBtyX3crNm1xHGjhzOVRePLXUoZma9yqdZ5kVgtqRZkmrJdoguzy0gaUKyDeB+YFWS8HcBt0ganST924DNhQu/NNbuPsL108ZTM6y3/9SYmZVev8k9IjqBRcBKson5GxGxUdIDkh5Iil0LbJS0heyoms8l+74AfAt4hWxzzTAg3Q3q/Th5+gyb9x1ze7uZpVpe4/giYgWwose6JTnvnwNm97HvI8Ajg4gxVTbuPUpnV3DDNCd3M0svTz8wQOubjwJwvZO7mZWIpx8ogt/75jp+uLWFF3//drLdCGZmpeH53Atow95jXHfZeCd2M0s1J/cBOHn6DK8faGPeVD8v1czSzcl9AF470EZnVzDvsj6buczMUiEVyb1cOlQ37DkGwLypTu5mVjruUC2wh5et57vr9rLukY+4zd3MSs4dqgWycc9R5k11Z6qZpZ+Te55On+li8/42N8mYWVlwcs9TU0s7HZ1dXHeZR8qYWfqlIrmXQ4fqhj3ZO1N95W5mpeYO1QL6k6c28cTqXWz44zs8G6SZpYI7VAtg64FjXD1lrBO7mZUFJ/c8bdnXxjWXjCt1GGZmeXFyz0Nr2ykOHe9gziXuTDWz8uDknoet+9sAmOMrdzMrE6lI7mkfLbNlf3baATfLmFka5DNaJq/kLmmBpK2SmiQ91Mv2iZKWSXpV0mpJ85L110ham/M6Jum3e+4/fvx4MpkMDQ0N+dduCG3Z30bduJFMGjuy1KGYmdHQ0EAmkwE42leZfh+zJ6kGeAz4MNAMvChpeURsyin2MLA2Ij4maU5S/raI2ArcmHOcPcCy86xPyWzd3+YmGTMrK/lcuc8HmiJie0R0AE8Cd/UoMxd4BiAitgAzJU3pUeY2YFtE7BxkzEPqTFfw2oE2rpni5G5m5SOf5D4V2J2z3Jysy7UOuBtA0nxgBjCtR5l7gCfOL8zSeePQcU51drm93czKSj7Jvbe7dnre1vooMFHSWuBBYA3QefYAUi3wS8A3eztBa2sr9fX1Z19JW1IqbNmXHSlz7aUeBmlmpZfJZM7mSmByX+X6bXMne6U+PWd5GrA3t0BEHAPuBVB2PtwdyavbncArEXGgtxPU1dWR1ukHtu4/xjDBVRePLXUoZmY0NjZ2j5RB0sG+yuVz5f4iMFvSrOQK/B5geW4BSROSbQD3A6uShN9tIWXYJAPZkTIzJ49h1IiaUodiZpa3fq/cI6JT0iJgJVADLI2IjZIeSLYvAa4FHpd0BtgE3Ne9v6TRZEfa/GYR4i+6Lfvb+DnPBGlmZSafZhkiYgWwose6JTnvnwNm97HvCWDSIGIsmeOnOtl1+AS/cnPPvmEzs3TzHarnsK21HYCrp7i93czSw/O5D9K3X2nmd76xjh/8zq1cdbGHQppZung+9/O0rbWdmmHi8ovGlDoUM7MBcXI/h6aWdmZMGk3tcH9MZlZenLXOoamlnavq3N5uZuUnFck9jR2qp890sfPQCd+8ZGap4w7VQWhqaef2v3uWv/34DfwXD4U0sxRyh+p5aGrJDoP0lbuZlSMn9z50j3G/0sndzMqQk3sftrW0c+n4UYwdmddNvGZmqeLk3oem1nau9EgZMytTqUjuaRstExFsa2l3e7uZpZJHy5ynfUff4uf/8j/501+ex6/fMqPU4ZiZ9cqjZQaoe6TMlXWedsDMypOTey88DNLMyp2Tey+aWtq5cNRw6saOLHUoZmbnJRXJPW0dqttas52p2cfBmpmliztUz1P9n/2AD15Tx19//IZSh2Jm1qdBd6hKWiBpq6QmSQ/1sn2ipGWSXpW0WtK8nG0TJH1L0hZJmyX9/PlXpfiOnjjNwfZTbm83s7LWb3KXVAM8BtwJzAUWSprbo9jDwNqIuB74JLA4Z9ti4HsRMQe4AdhciMCLpanVnalmVv7yuXKfDzRFxPaI6ACeBO7qUWYu8AxARGwBZkqaIulC4FbgK8m2jog4UrDoi2CbR8qYWQXIJ7lPBXbnLDcn63KtA+4GkDQfmAFMA64AWoF/lrRG0pcl/czg8dbWVurr68++MpnMeVSlMJpa26kdPoxpE0eXLAYzs75kMpmzuRKY3Fe5fGbF6m3ISM9e2EeBxZLWAuuBNUAnMAK4CXgwIl6QtBh4CPiD3J3r6upIS4dqU0s7V0weQ80wj5Qxs/RpbGzsHimDpIN9lcsnuTcD03OWpwF7cwtExDHg3uRkAnYkr9FAc0S8kBT9FtnknlrbWtuZN7XP0UVmZmUhn2aZF4HZkmZJqgXuAZbnFkhGxNQmi/cDqyLiWETsB3ZLuibZdhuwqUCxF9zJ02fYffiEZ4M0s7LX75V7RHRKWgSsBGqApRGxUdIDyfYlwLXA45LOkE3e9+Uc4kHg60ny305yhZ9GOw4epyvcmWpm5S+vJ1FExApgRY91S3LePwfM7mPftUCvg+y7dd+h2tDQQENDQz4hFUX305eu8pW7maXYU0891X1Hv+9QzccXf/Aai595nc1/soBRI2pKHY6Z2Tl5yt88NbW0M23iBU7sZlb2nNxzNLW0u0nGzCqCk3viTFew4+Bxd6aaWUVwck/sefMtTnV2ObmbWUVIRXJPw3zuTa1tAB7jbmap5/ncB+BLq7bz5ys2s/YPP8yE0bX972BmVmIeLZOHppZ2Jo+tdWI3s4rg5J5oam3nCjfJmFmFcHIHIiI7DNKdqWZWIVKR3EvdoXroeAdH3zrtMe5mVhbcoZqn57cf4p7M8zz+6fncenVdyeIwMxsId6j2oyl5tN6VbpYxswrh5E52NsjRtTVcNn5UqUMxMysIJ3eyV+5X1o0l+xApM7Pyl4rkXuoO1W0t7VxZ9zPP7TYzSyV3qObh+KlOrntkJb/7katZ9KFenzdiZpZKg+5QlbRA0lZJTZJ+5gHXkiZKWibpVUmrJc3L2faGpPWS1koq/RM5etjeehzwo/XMrLL0+5g9STXAY8CHgWbgRUnLIyL3QdcPA2sj4mOS5iTlb8vZ/sGIOFjAuAume8IwJ3czqyT5XLnPB5oiYntEdABPAnf1KDMXeAYgIrYAMyVNKWikRdLU0k7NMHH5RW5zN7PKkU9ynwrszlluTtblWgfcDSBpPjADmJZsC+D7kl6W1Di4cAuvqaWdGZNGUzs8FX3LZmYF0W+zDNDb+MCevbCPAoslrQXWA2uAzmTbeyJir6SLgaclbYmIVbk7t7a2Ul//0z6BxsbG7p7gomtqaWe2m2TMrExkMhkymUz34uS+yuWT3JuB6TnL04C9uQUi4hhwL4Cyg8V3JC8iYm/ys0XSMrLNPG9L7nV1dZRitExHZxc7D51gwbxLhvzcZmbnI/fiV1KffZn5tEW8CMyWNEtSLXAPsDy3gKQJyTaA+4FVEXFM0hhJ45IyY4CPABsGXJsi2XnoOJ1d4c5UM6s4/V65R0SnpEXASqAGWBoRGyU9kGxfAlwLPC7pDLAJuC/ZfQqwLLnzczjwLxHxvcJX4/x0zylzVd24EkdiZlZY+TTLEBErgBU91i3Jef8c8DN3AEXEduCGQcZYND+dMMwjZcyssqRiiEipph94vaWdqRMuYHRtXn/jzMxSwdMP9OOji39M3biRfPXT84f83GZmg+X53HvR1RVsP+hhkGZWmao2ue858hYnT3d5pIyZVaSqTe6vt3hOGTOrXKlI7qXoUD07DNLJ3czKjDtUz+Hz31rHf25p5aUv3D6k5zUzKxR3qPbi9ZZ2rvL4djOrUFWZ3COCppZ2N8mYWcWqyuTe2naKtpOdXFXn5G5mlSkVyX2oO1RfTzpTZ0/xnDJmVn7codqHr/7kDR5ZvpEXHr6NKReOGrLzmpkVkjtUe2hqaWfcyOFcPG5kqUMxMyuKqkzuWw+0MXvKWJKpiM3MKk7VJfeI4LUDbVxzyYWlDsXMrGiqLrm3tJ3iyInTXDPFI2XMrHKlIrkP5WiZLfuzc8r4yt3MylU+o2XySu6SFkjaKqlJ0kO9bJ8oaZmkVyWtljSvx/YaSWskfbe3448fP55MJkNDQ0M+4QzK1v3HAJhziYdBmll5amhoIJPJABztq0y/yV1SDfAYcCcwF1goaW6PYg8DayPieuCTwOIe2z8HbM4/9OLZsr+Ni8eNZOKY2v4Lm5mVqXyu3OcDTRGxPSI6gCeBu3qUmQs8AxARW4CZkqYASJoG/ALw5YJFPQjZzlRftZtZZcsnuU8FducsNyfrcq0D7gaQNB+YAUxLtn0R+DzQNahIC+BMV/D6gXau8Z2pZlbh8knuvQ0G73lb66PARElrgQeBNUCnpF8EWiLi5XOdoLW1lfr6+rOvpC2p4N44dJxTnV2+cjezspXJZM7mSmByX+WG53GsZmB6zvI0YG9ugYg4BtwLoOydQTuS1z3AL0n6KDAKuFDS1yLiE7n719XVMRTTD2xNRsrM8UgZMytTjY2N3SNlkHSwr3L5XLm/CMyWNEtSLdmEvTy3gKQJyTaA+4FVEXEsIv5nREyLiJnJfv/ZM7EPpa3725D89CUzq3z9XrlHRKekRcBKoAZYGhEbJT2QbF8CXAs8LukMsAm4r4gxn7et+9uYOWkMF9TWlDoUM7OiyqdZhohYAazosW5JzvvngNn9HONHwI8GHGEBbT3Q5s5UM6sKVXOHatvJ0+w4eJy5l7m93czKm+dzz7F6x2F+9Z+eY+lv1POhOVOKei4zs6Hg+dyBjXuzd+nOu6zPP3RmZhWjapL7hj3HmDx2JBf7yUtmVgWqJrlv3HuUeVPd3m5m1SEVyb3YHaonT5/h9ZZ2N8mYWUVwh2pi3e4j3PXY/2PJJ25iwbxLi3YeM7OhVPUdqhuSztTrfOVuZlWiOpL7nmOMv2AE0yZeUOpQzMyGRFUk91ebjzBv6oVk5zQzM6t8FZ/cT3R0smV/GzddPrHUoZiZDZlUJPdijpZ5tfkoZ7qCd1w+oeDHNjMrBY+WAf7hR0381fe2suYPPuznpppZRanq0TKv7DzCFZPHOLGbWVWp6OQeEazd/SY3uknGzKpMRSf33Yff4mB7hztTzazqpCK5F6tD9eVdhwHcmWpmFaXqO1Q//611rNx4gDV/8GGGDfMYdzOrLIPuUJW0QNJWSU2SHupl+0RJyyS9Kmm1pHnJ+lHJ8jpJGyX98eCqMjA/2XaIW664yIndzKpOv8ldUg3wGHAnMBdYKGluj2IPA2sj4nrgk8DiZP0p4EMRcQNwI7BA0i2FCv5cdh8+QfObb/HuKycPxenMzFIlnyv3+UBTRGyPiA7gSeCuHmXmAs8ARMQWYKakKZHVnpQZkbyGpB3oJ9sOAvDuKycNxenMzFIln+Q+Fdids9ycrMu1DrgbQNJ8YAYwLVmukbQWaAGejogXep6gtbWV+vr6s69MJjPwmvTwk22HmDx2JFddPHbQxzIzS4tMJnM2VwJ9Nk0Mz+NYvTVY97z6fhRYnCTx9cAaoBMgIs4AN0qaACyTNC8iNuTuXFdXRyE7VM90Batea+X9V9d5sjAzqyiNjY3dI2WQdLCvcvkk92Zges7yNGBvboGIOAbcm5xMwI7klVvmiKQfAQuAtyX3Qntl15u8eeI0t8+dUszTmJmlVj7NMi8CsyXNklQL3AMszy0gaUKyDeB+YFVEHJNUl1yxI+kC4HZgS+HC790PNh1gRI249eq6Yp/KzCyV+r1yj4hOSYuAlUANsDQiNkp6INm+BLgWeFzSGWATcF+y+6XAV5MRN8OAb0TEd4tQj9x4eXrTAW65YhIXjhpRzFOZmaVWPs0yRMQKYEWPdUty3j8HzO5lv1eBdwwyxgHZuPcY2w8e5973zhrK05qZpUrFTT/w7Vf2MKJGNFzvB2GbWWWquukHOjq7ePejz1A/4yKW/PrNBYjMzCy9qmY+939bu4eD7R382rsuL3UoZmYlVTHJvasr+NKPtzPnknG8b7anHDCz6lYxyf2bL+/mtQPt/NYHr/KNS2ZW9VKR3Afbobrz0HH+/N83886ZE92RamYVr+I7VN84eJz/23SQx37YxFunz/DUovcy/aLRRYjQzCx9KrZDdeXG/XzhOxsYM3I4/3L/LU7sZmaJvG5iSqu7b5rGHdddwuUXjfYDOczMcpR1cq8bN5K6cSNLHYaZWeqkolmmWA/INjOrRBXfoWpmVs0qtkPVzMx65+RuZlaByj65F+J5q+WiWurqelaWaqknpKuuqUjuPTtUB9KxOtAPc6CdtgMpX8xjw8DqWuxYivm5uJ5DH0sxy/t3tPCx5NOhmorkPn78eDKZDA0NDcDAKz0Q5fwPp5jHTtPnUqw4il2+2KO90vRvsVy/00qpZ0NDQ/cfkqN9lU/FaBlJrcDOnFXjOUfQPUwG+nwCeC8GcuyBli/msWFgdS12LMX8XFzPoY+lmOX9O1q8WGZERK8Pi05Fcjczs8JKRbOMmZkVlpO7mVkFKuvkLmmBpK2SmiQ9VOp4ikXSG5LWS1orqaJu5ZW0VFKLpA056y6S9LSk15OfE0sZYyH0Uc8/krQn+V7XSvpoKWMsBEnTJf1Q0mZJGyV9LllfUd/pOeqZmu+0bNvcJdUArwEfBpqBF4GFEbGppIEVgaQ3gPqIGEinVFmQdCvQDjweEfOSdX8FHI6IR5M/2hMj4n+UMs7B6qOefwS0R8TflDK2QpJ0KXBpRLwiaRzwMvDLwG9QQd/pOer5q6TkOy3nK/f5QFNEbI+IDuBJ4K4Sx2QDFBGrgMM9Vt8FfDV5/1WyvzRlrY96VpyI2BcRryTv24DNwFQq7Ds9Rz1To5yT+1Rgd85yMyn7cAsogO9LellSY6mDGQJTImIfZH+JgItLHE8xLZL0atJsU9ZNFT1Jmgm8A3iBCv5Oe9QTUvKdlnNy7+3pHOXZxtS/90TETcCdwGeT/+Jb+ftH4ErgRmAf8LelDadwJI0F/hX47Yg4Vup4iqWXeqbmOy3n5N4MTM9ZngbsLVEsRRURe5OfLcAysk1SlexA0qbZ3bbZUuJ4iiIiDkTEmYjoAr5EhXyvkkaQTXhfj4hvJ6sr7jvtrZ5p+k7LObm/CMyWNEtSLXAPsLzEMRWcpDFJhw2SxgAfATace6+ytxz4VPL+U8C/lTCWoulOdomPUQHfqyQBXwE2R8Tf5WyqqO+0r3qm6Tst29EyAMkwoy8CNcDSiPjzEodUcJKuIHu1DtnHIv5LJdVT0hPAB8jetn0AeAT4DvAN4HJgF/DxiCjrzsg+6vkBsv99D+AN4De726XLlaT3Aj8G1gNdyeqHybZHV8x3eo56LiQl32lZJ3czM+tdOTfLmJlZH5zczcwqkGyKgpIAAAAoSURBVJO7mVkFcnI3M6tATu5mZhXIyd3MrAI5uZuZVSAndzOzCvT/AeFIhkxcfkqWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ms.Abar, ms.U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20e8bafe4f0>]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD1CAYAAACrz7WZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAcnElEQVR4nO3deXAc53nn8e+DiyAOggABghTAU6YOSiJ1QJRsyTYZrxVKWa4sZ9cW5TiWYgvxxqrEW/FWHG/FSmznLNuVdUWWduQwsnd1lFJrxmKF1rGOI9qWZJOUZYm0RBIiKRK8ABIkSIgXjmf/mAEJUTMYHNPTPT2/T9UUZ/rt6X7eaeJB4+2n3zF3R0RE4qUk7ABERCT3lNxFRGJIyV1EJIaU3EVEYkjJXUQkhpTcRURiqCzsAAAaGxt9/vz551739vZSV1c3pvd2d3fT1NQ05n2NZ9vjXT/IbcP4+hp0LEF+Lupn/mMJcn39jAYXy+bNmw+7e/odunvoj+uuu85Huvfee32sLnxvNuPZ9njXD3Lb7uPra9CxBPm5qJ/5jyXI9fUzGlwswCbPkFcjOSyzatWqyGx7POsHue3xCjqWID+XoOIIev0g+zne7UfpcxmvKB2jQu0ngHkE7lBdtGiRr1ixglWrVo27A21tbWzatCmgyKKlWPqqfsZLsfQT8tfXdevWsW7dOh5++OEOd1+Ubp1IjLnX1dWRSCQm9N729vYcRxNdxdJX9TNeiqWfkL++Dp8IP/zww72Z1onEmXtbW5sXy292EZFcMbPN7t6Wri3rmLuZrTGzLjPbkqH9djN71cxeMbNNZnbziLaVZrbNzDrM7IsT74KIiIzHWC6oPgKsHKX9R8BSd78a+D3gOwBmVgo8ANwKLAZWm9niSUUrIiJjkjW5u/sGoGeU9j4/P7ZTDQw/XwZ0uPtOdz8LPAHcPsl43+HU2UHODAzmcpMiIrGQk1JIM7vDzN4A/pXk2TtAC7B3xGqdqWU58+SmvVz9F8/xl//6awaHwr92ICISFTlJ7u6+1t0vAz4CfDW12NKtmu793d3dtLW1nXuMtXJmSWsdv3lFMw//ZBf/80c7JhS7iEghSSQS53Il0JhpvZyWQrr7BjO72MwaSZ6pzxnR3ArsT/e+pqamCdWGXjO3nmvm1gPw4L93cOf1c7ho+tTxBy4iUiDa29vPlVya2eFM6036zN3M3mNmlnp+LVABHAE2AovMbIGZVQB3Ak9Ndn/p/PeVlwHw0PNvBrF5EZGCM5ZSyMeBF4FLzazTzD5tZp81s8+mVvltYIuZvUKyOubjqWkPBoD7gGeA14En3X1rEJ1omT6V/7jkIta+vI/T/brAKiKSdVjG3Vdnaf9b4G8ztK0H1k8stPH5z9e1svaX+3j214f4T0svyscuRUQiK5ITh03EexfOoKl2Cs9uPRh2KCIioYtEcu/t7aW9vZ1169ZNeBslJcbyS5rYsL2bgcGhHEYnIhIt69atG76omnFC+FjNLbP+tQP8waMv88+ffS/Xz2/IQWQiItE1qbllCsnNixopKzGe39YddigiIqGKVXKfVlnOFS11/GJ3xtkSRESKQqySO8D18+r51d5jmnNGRIpa7JJ72/wGzgwMsWVfxjnsRURiLxLJPRfVMsPa5ienI9i4++iktyUiEkVFVy0z7AN/92OubJnGtz9xXc62KSISNUVTLTPsqpY6tuw7HnYYIiKhiWVyv6JlGnt6TtJ7sj/sUEREQhHL5H5VS3IYast+XVQVkeIUy+R+5UWp5K6KGREpUrFM7vXVFbRMn8prSu4iUqQikdxzWQo5bPFF03jj4ImcbU9EJCqKthQS4OvPbOOh599k61d+kyllpTndtohIFBRdKSTAouYaBoacXYffDjsUEZG8i21yv6S5FoDth/pCjkREJP9im9wXNlVTWmLsOKRxdxEpPrFN7lPKSpk/o4rtSu4iUoQikdyDqJaB5NCMhmVEJG6KuloG4JvPbecf/m0Hv/7KSirLVTEjIvEyqWoZM1tjZl1mtiVD+yfM7NXU4wUzWzqibbeZvWZmr5hZ7rN3Fpc01zDksLNbFTMiUlzGMizzCLBylPZdwAfdfQnwVSBxQfsKd78602+XIC1srAFQOaSIFJ2ybCu4+wYzmz9K+wsjXr4EtE4+rNyY31gFwK7DGncXkeKS6wuqnwZ+OOK1A8+a2WYza8/xvrKqqihjdl0lO3XmLiJFJuuZ+1iZ2QqSyf3mEYtvcvf9ZjYTeM7M3nD3DRe+t7u7m7a286M27e3tw1eCJ23+jGp2K7mLSEwkEgkSiXOj342Z1stJcjezJcB3gFvd/cjwcnffn/q3y8zWAsuAdyX3pqYmgqiWAVjQVM0PXzsQyLZFRPJt5MmvmR3OtN6kh2XMbC7wfeCT7r59xPJqM6sdfg7cAqStuAnSwsZqjp7s5+jbZ/O9axGR0GQ9czezx4HlQKOZdQL3A+UA7v4Q8GVgBvBtMwMYSFXGNANrU8vKgMfc/ekA+jCqBY3VAOw68jb11RX53r2ISCjGUi2zOkv7Z4DPpFm+E1j67nfk17nk3v02186tDzkaEZH8iMT0A0Ga01BFaYmp1l1EikokkntQc8sAlJeWMLehSsldRGKj6OeWGXbPP/2CQ8fPsP6P3h/YPkRE8q0ov4lppHkzqtnbc5Io/CITEcmHokjucxqqOHFmgGMn+8MORUQkL4oiuc9tSM4xs6fnZMiRiIjkh5K7iEgMFUVyn9MwFVByF5HiEYnkHmQpJCRnh2ysmcJeJXcRiQGVQo5wx7d/xtTyUh6798ZA9yMiki9FXwoJyXF3DcuISLEoquS+/9gp+geHwg5FRCRwRZPc5zRUMeSw/9ipsEMREQlc0SR3lUOKSDGJRHIPuloGlNxFJD5ULTPC4JBz+Z89zT03z+dPb7080H2JiOSDqmWA0hKjtX6qat1FpCgUTXKH5EVVDcuISDEoquTeWj+VzqOqlhGR+Cuy5F7FsZP99J0ZCDsUEZFARSK556NaBqClPjmB2D6dvYtIAVO1zAVe3nOUj377Bdbc3cZvXNYc+P5ERIKkapmU1uk6cxeR4pA1uZvZGjPrMrMtGdo/YWavph4vmNnSEW0rzWybmXWY2RdzGfhENNZMoaK0RBdVRST2xnLm/giwcpT2XcAH3X0J8FUgAWBmpcADwK3AYmC1mS2eVLSTVFJitNRPpVPzy4hIzGVN7u6+AegZpf0Fdz+aevkS0Jp6vgzocPed7n4WeAK4fZLxTlrLdJVDikj85XrM/dPAD1PPW4C9I9o6U8vepbu7m7a2tnOPRCKR47DOa62fqjF3ESlYiUTiXK4EGjOtV5arHZrZCpLJ/ebhRWlWS1ua09TURD6qZSB55n647wyn+wepLC/Nyz5FRHKlvb19uAwSMzucab2cnLmb2RLgO8Dt7n4ktbgTmDNitVZgfy72NxmtqS/L3qdxdxGJsUkndzObC3wf+KS7bx/RtBFYZGYLzKwCuBN4arL7m6yW6cmpfzU0IyJxlnVYxsweB5YDjWbWCdwPlAO4+0PAl4EZwLfNDGDA3dvcfcDM7gOeAUqBNe6+NZBejMPwXaq6qCoicZY1ubv76iztnwE+k6FtPbB+YqEFo7l2CmUlxr5jmh1SROKrqO5QBSgrLWFWXaXO3EUk1iKR3PM1cdgwlUOKSCHTxGEZ/PGTv+JnHYd56Usfyts+RURyTROHXaC1fiqHTpzm7MBQ2KGIiASiKJN7S/1U3OFAr4ZmRCSeijK5t+pLO0Qk5oozuaduZNLskCISV5FI7vmulplVV4mZbmQSkcKkaplRvPevf8T7Lm7kGx9bmn1lEZEIUrVMGsl53XWXqojEU9Em99Z6fWmHiMRX0Sb3lvqpHDx+moFB1bqLSPwUbXJvra9icMg5dOJM2KGIiORcJJJ7vqtlIDnmDqp1F5HCo2qZUbzZ3ceHvvE83/zYUj56bWv2N4iIRIyqZdLQmbuIxFnRJvfK8lIaa6aoYkZEYqlokzuk5nXXFAQiEkNFndxblNxFJKYikdzDqJaB89/INDQU/kVlEZGxUrVMFv/7xd382Q+28osvfYiZ0yrzvn8RkclQtUwGLal53ffqoqqIxEzW5G5ma8ysy8y2ZGi/zMxeNLMzZvaFC9p2m9lrZvaKmeX/1DyL1vrkvO4adxeRuBnLmfsjwMpR2nuAPwS+nqF9hbtfnelPhzCp1l1E4iprcnf3DSQTeKb2LnffCPTnMrB8qJ5SRn1Vuab+FZHYCXrM3YFnzWyzmbUHvK8JUTmkiMRRWcDbv8nd95vZTOA5M3sj9ZfAO3R3d9PWdn7Upr29fbjMJ3Ct06vo6O7Ly75ERCYrkUiQSCSGXzZmWi/Q5O7u+1P/dpnZWmAZ8K7k3tTURBilkJA8c39+ezfujpmFEoOIyFiNPPk1s8OZ1gtsWMbMqs2sdvg5cAuQtuImTC3Tp3Kqf5CjJwvukoGISEZZz9zN7HFgOdBoZp3A/UA5gLs/ZGazgE3ANGDIzD4PLCb558La1NlwGfCYuz8dRCcmozVV69559CQN1RUhRyMikhtZk7u7r87SfhBINyH6cWDpBOPKm+EbmfYdPcWS1ukhRyMikhtFfYcqnL+RSVP/ikicRCK5hzVxGEDd1HJqp5SpHFJECoYmDhujlX+/gdb6Kr7zqcjdRCsikpEmDsuitX6q7lIVkVhRcidZDqlhGRGJEyV3khdVT5weoPeUat1FJB6U3HlnOaSISBxEIrmHWS0D529k0tCMiBQCVcuM0ZG+M1z3tf/H/asWc89NC0KLQ0RkPFQtk0VDdQWV5SUalhGR2FByB8yM1voq3aUqIrGh5J6ickgRiRMl95RWfSOTiMSIkntKS/1Uet4+y8mzA2GHIiIyaZFI7mGXQsL52SF1UVVEok6lkOOw+a2j/PaDL/BPd1/PistmhhqLiMhYqBRyDM59I5PG3UUkBpTcU5pqplBRWqLZIUUkFpTcU0pKjIumV2rMXURiQcl9hDkNVexVcheRGIhEco9CtQykknuPhmVEJNpULTNO/+v5N/nrH77Ba39+C7WV5WGHIyIyqklVy5jZGjPrMrMtGdovM7MXzeyMmX3hgraVZrbNzDrM7IsTCz9/5jYka93fOqKzdxEpbGMZlnkEWDlKew/wh8DXRy40s1LgAeBWYDGw2swWTyzM/Jg7I5nc92hoRkQKXNbk7u4bSCbwTO1d7r4RuPA76pYBHe6+093PAk8At08m2KANn7kruYtIoQvygmoLsHfE687UssiqrSynobpCwzIiUvDKAty2pVmW9uptd3c3bW3nrwm0t7cPXwnOu7kNVezpeTuUfYuIZJNIJEgkEsMvGzOtF2Ry7wTmjHjdCuxPt2JTUxNRqJYBmDejis1vHQ07DBGRtEae/JrZ4UzrBTkssxFYZGYLzKwCuBN4KsD95cS8hir2HztF/+BQ2KGIiExY1jN3M3scWA40mlkncD9QDuDuD5nZLGATMA0YMrPPA4vd/biZ3Qc8A5QCa9x9azDdyJ05DVUMeXLq3/mN1WGHIyIyIVmTu7uvztJ+kOSQS7q29cD6iYUWjnkzkgn9rZ6TSu4iUrAiMf1AlMxTrbuIxEAkkntU5paB5NS/U8pK2HNEFTMiEk2aW2aCPvzN51nQWE3id9NO2SAiEgn6JqZxmjejSsMyIlLQlNzTmNOQTO5R+KtGRGQilNzTmNdQxcmzgxzuOxt2KCIiE6Lknsa8VAnkbl1UFZECpeSexsJUct/VreQuIoUpEsk9SqWQAK31VVSUlvDm4b6wQxEReReVQk7Ch7/5PPMbq3lY5ZAiElEqhZyAhU3V7OzWmbuIFCYl9wwWNNawp+ckA5odUkQKkJJ7BgubqukfdDqPngo7FBGRcVNyz+DipmTFzE5dVBWRAhSJ5B61ahmAhY01AOxUOaSIRIyqZSbp6q88y21Xzeav7rgq7FBERN5F1TITtLBRFTMiUpiU3EexsKlGwzIiUpCU3EexsKmarhNnOHG6P+xQRETGRcl9FLqoKiKFKhLJPYrVMgCXNCeT+44ujbuLSHSoWmaSBoecy7/8NHe/bz5fuu3ysMMREXkHVctMUGmJsWhmDdsOngg7FBGRccma3M1sjZl1mdmWDO1mZt8ysw4ze9XMrh3RttvMXjOzV8wseqfmY3Bpcy3bDym5i0hhGcuZ+yPAylHabwUWpR7twIMXtK9w96sz/ekQdZfMquVA72l6T6liRkQKR9bk7u4bgJ5RVrkd+J4nvQRMN7PZuQowbJc21wKwQ2fvIlJAcjHm3gLsHfG6M7UMwIFnzWyzmbVn2kB3dzdtbW3nHolEIgdh5cYls5LJfZuSu4hEQCKROJcrgcZM65XlYF+WZtlwCc5N7r7fzGYCz5nZG6m/BN6hqamJKFbLAFxUV0nNlDK266KqiERAe3v7cBkkZnY403q5OHPvBOaMeN0K7Adw9+F/u4C1wLIc7C+vzIxLmmt05i4iBSUXyf0p4HdTVTM3Ar3ufsDMqs2sFsDMqoFbgLQVN1F36axath08QRTuCRARGYuswzJm9jiwHGg0s07gfqAcwN0fAtYDtwEdwEngntRbm4G1Zja8n8fc/ekcx58XlzbX8vgv9tJ14gzN0yrDDkdEJKusyd3dV2dpd+BzaZbvBJZOPLTouLIleYfvln29Su4iUhB0h+oYXD57Gmbw2r7esEMRERmTSCT3qE4cNqx6ShkXN9WwZd/xsEMREdHEYbn0+Sd+yc939fDin34o7FBERABNHJYTV7bUcaD3NIf7zoQdiohIVkruY3TFRecvqoqIRJ2S+xhd0TINUHIXkcKg5D5G0yrLWdBYza86ldxFJPoikdyjXi0z7Jq50/nlnqO6U1VEQqVqmRx79Odv8T/WbuHfv7Cc+Y3VYYcjIkVO1TI50javAYBNbx0NORIRkdEpuY/Dopk1TKssY7OSu4hEnJL7OJSUGNfOq2fzW6N9MZWISPiU3Mfpurn1bD/UR+9JfaeqiERXJJJ7oVTLALTNHx5319m7iIRD1TIBON0/yNK/eJa7bpjL/auuCDscESliqpbJocryUpYtaOCnOzJ+daGISOiU3Cfg/Ysa2dHVx8He02GHIiKSlpL7BNz8niYAftqhs3cRiSYl9wm4bFYtjTUVbNjeHXYoIiJpRSK5F1K1DCTr3VdcOpMfb+vi7MBQ2OGISJFRtUyAfvT6IT793U08cs/1LL90ZtjhiEgRUrVMAG56TyM1U8p4ZuvBsEMREXmXrMndzNaYWZeZbcnQbmb2LTPrMLNXzezaEW0rzWxbqu2LuQw8bJXlpay4bCbPbj3EwKCGZkQkWsZy5v4IsHKU9luBRalHO/AggJmVAg+k2hcDq81s8WSCjZrfumoWR94+q6oZEYmcrMnd3TcAo91rfzvwPU96CZhuZrOBZUCHu+9097PAE6l1Y+M3LmumobqCJzftDTsUEZF3KMvBNlqAkdmtM7Us3fIbcrC/yKgoK+GOa1r43ou7OdJ3hhk1U0KN58Tpfvb0nOTQ8dMc7D3DweOnOXbyLH2nB+g7c/5xpn+I/qEh+geHGBh0+gedgaHk88Ehx0leZHeHc5fbnbTL3X3E83z2ViQevvaRK/mdG+flfLu5SO6WZpmPsvxduru7aWs7f8G3vb19uMwn8j5+/Rz+8ae7eGLjXj634j1522/vyX427u5h41s9vH7gBDsOneDABXfMmiW/+7W2soyaKclHfVUFleUllJcmH2UlRllpCeWlRnlpCaUlljxwqaOXeoXZ+QNqI9rOP39nY7qDLyLvtqQ1YzVjWolEgkQiMfyyMdN6uUjuncCcEa9bgf1ARYbl79LU1EShlUIOu6S5luWXNrHmp7v4vZsWMLWiNJD9uDs7uvp4estBntl6kF8fOI47lJcalzTXcuPCGSxqrmHBjGpm1VXSPK2SptoplJeqIEokTkae/JpZxgt+uUjuTwH3mdkTJIddet39gJl1A4vMbAGwD7gTuCsH+4uc//rBi/l44iWe2LiHe25akNNtHz/dz9qX9/HYz/ew7dAJzJJzyv+3/3AJyxY0cPWc6VSWB/MLRUQKV9bkbmaPA8uBRjPrBO4HygHc/SFgPXAb0AGcBO5JtQ2Y2X3AM0ApsMbdtwbQh9AtW9DAjQsb+NaPdvDRa1qpqyqf9Db39pzkweffZO3L+zjVP8iS1jq+cvsVrLxiFjOnVeYgahGJM92hmiOvHzjOb33rJ9x1w1y+9pGrJrydt468zQM/7uD7L++jxIw7rmnhEzfOZUnr9BxGKyJxMNodqrkYlhHg8tnTuPt9C1jzs118YFETt1wxa1zv33X4bf7h3zr4l1f2UVpi/M6N8/jsBy9mVp3O0kVk/CKR3IcnDlu1ahWrVq0KO5wJ+5NbL2Xj7h7++Mlf8X8+U8nSOdnPtju6+njgxx384JV9VJSVcPf75vP7H1iooRcRyWjdunXDEy1q4rB82X/sFB9PvMjRt/v5q49exaolszF7Z2Ggu/PiziM8+tIe1m85QGVZKZ987zzuff9CmmrDrZUXkcIx2rCMknsA9h87xR88+jKv7D3GVS113LK4mZb6qfSdGeD1AyfYsL2bfcdOMa2yjLtumMe9718Q+g1QIlJ4lNxD0D84xD9v6uTRn7/F1v3Hzy2vmVLGjQsbuO2q2dx21WyVMYrIhOmCagjKS0u464a53HXDXI6f7udI31lqppQxo7qCkhLdvykiwVJyz4NpleVMq5x87buIyFhF4t70QvuaPRGRMOlr9kREYkxfsyciUmSU3EVEYqjgk/uIeY1jr1j6qn7GS7H0E6LV10gm9/FcWB3vhznei7bjWT/IbcP4+hp0LEF+Lupn/mMJcn39jIYTS8En96C3HaX/OEFuO0qfS1BxBL1+0NVeUfq/WKjHtFj6CRGpljGzXuAscAzoJVne0zvGtzcCGb+NJI3xbHu86we5bRhfX4OOJcjPRf3MfyxBrq+f0dzHUgdMByrcPW05ZCSSu4iI5FYkh2VERGRylNxFRGKooJO7ma00s21m1mFmXww7nqCY2W4ze83MXjGzWN3Ka2ZrzKzLzLaMWNZgZs+Z2Y7Uv/VhxpgLGfr552a2L3VcXzGz28KMMRfMbI6Z/djMXjezrWb2R6nlsTqmo/QzMse0YMfczawU2A58GOgENgKr3f3XoQYWADPbDbS5+3guShUEM/sA0Ad8z92vTC37O6DH3f8m9Uu73t3/JMw4JytDP/8c6HP3r4cZWy6Z2Wxgtru/bGa1wGbgI8DdxOiYjtLPjxGRY1rIZ+7LgA533+nuZ4EngNtDjknGyd03AD0XLL4d+G7q+XdJ/tAUtAz9jB13P+DuL6eenwBeB1qI2TEdpZ+RUcjJvQXYO+J1JxH7cHPIgWfNbLOZtYcdTB40u/sBSP4QATNDjidI95nZq6lhm4IeqriQmc0HrgF+ToyP6QX9hIgc00JO7um+8aIwx5iyu8ndrwVuBT6X+hNfCt+DwMXA1cAB4BvhhpM7ZlYD/F/g8+5+PNv6hSpNPyNzTAs5uXcCc0a8bgX2hxRLoNx9f+rfLmAtySGpODuUGtMcHtvsCjmeQLj7IXcfdPch4GFiclzNrJxkwnvU3b+fWhy7Y5qun1E6poWc3DcCi8xsgZlVAHcCT4UcU86ZWXXqgg1mVg3cAmwZ/V0F7yngU6nnnwJ+EGIsgRlOdil3EIPjamYG/CPwurt/c0RTrI5ppn5G6ZgWbLUMQKrM6O+BUmCNu/9lyCHlnJktJHm2DsmvRXwsTv00s8eB5SRv2z4E3A/8C/AkMBfYA/wXdy/oi5EZ+rmc5J/vDuwGfn94XLpQmdnNwE+A14Ch1OIvkRyPjs0xHaWfq4nIMS3o5C4iIukV8rCMiIhkoOQuIhJDSu4iIjGk5C4iEkNK7iIiMaTkLiISQ0ruIiIxpOQuIhJD/x9h0uAp3LcORwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ms.Abar, ms.rho(ms.R,ms.m))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tracing ray is enabled and excision will be performed!\n"
     ]
    }
   ],
   "source": [
    "ms.run_steps(53000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "ms_mix = copy.copy(ms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "ms_uni = copy.copy(ms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "hm_uni = HM(ms_uni)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix = HM(ms_mix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#hm_mix.run_steps(2000)\n",
    "hm=hm_mix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_uni.run_steps(100)\n",
    "hm=hm_uni"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hm_mix.run_steps(40000)\n",
    "hm=hm_mix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_uni=copy.copy(hm_uni_bak)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_uni_bak=copy.copy(hm_uni)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix_bak=copy.copy(hm_mix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix = copy.copy(hm_mix_bak)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(hm_mix_bak.u -hm_mix.u)/ hm_mix.deltau_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix.u/hm_mix.deltau_i*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm.u=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_uni.u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.exp(hm.alpha * hm.xi)\n",
    "H = np.exp(-hm.xi) / hm.RH\n",
    "rho_b = a**(1+hm.w)\n",
    "Rb = a * hm.A\n",
    "plt.plot(hm.R * Rb ,  (np.exp(-hm.xi/2) * hm.R **3 * hm.Abar**3 * hm.m ) / 2 )\n",
    "plt.xlim(0,10)\n",
    "plt.ylim(0,20)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix.U.argmin()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix.Abar[600] - hm_mix.Abar[599]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm_mix.Abar), hm_mix.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "plt.xlim(5,7)\n",
    "plt.ylim(0.2,1)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm.Abar), hm.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "#plt.xlim(20,28)\n",
    "#plt.ylim(-0.95,1.2)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm.Abar), hm.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "#plt.xlim(20,28)\n",
    "#plt.ylim(-0.95,1.2)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "#plt.xlim(0,5)\n",
    "#plt.ylim(-60,10)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm_mix.Abar), hm_mix.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "#plt.xlim(0,5)\n",
    "#plt.ylim(-60,10)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm_mix.Abar), hm_mix.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "#plt.xlim(0,5)\n",
    "#plt.ylim(-60,10)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix.U[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter((hm.Abar*hm.R * np.exp((hm.alpha) * hm.xi)), hm.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "plt.xlim(0,50)\n",
    "plt.ylim(-60,10)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix_bak_low.Abar[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix.Abar[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm_mix.Abar*hm_mix.R * np.exp((hm_mix.alpha) * hm_mix.xi))*hm_mix.RH, hm_mix.U)\n",
    "#plt.plot((hm_uni.Abar), hm_uni.U)\n",
    "\n",
    "\n",
    "\n",
    "plt.xlim(0,54)\n",
    "#plt.ylim(-60,10)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms=ms_mix_cubic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((ms.Abar*ms.R * np.exp((ms.alpha-1) * ms.xi)), ms.rho(ms.R,ms.m))\n",
    "\n",
    "plt.ylim(1e-1,1e4)\n",
    "plt.xlim(0,1)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm_mix.Abar*hm_mix.R * np.exp((hm_mix.alpha-1) * hm_mix.xi)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(np.arange(hm.N), r)\n",
    "\n",
    "#plt.ylim(0.5,1)\n",
    "plt.xlim(1470,1500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(np.arange(hm.N), WENO_dfdA(r, hm.deltaL, hm.Abar, 1e100)[:])\n",
    "\n",
    "#plt.ylim(0.5,1)\n",
    "plt.xlim(0,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(np.arange(hm.N), hm.U)\n",
    "\n",
    "#plt.ylim(0.5,1)\n",
    "plt.xlim(0,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        xiprime = dfdA(hm.xi, hm.deltaL, hm.Abar, 1e100)\n",
    "        Rprime = dfdA(hm.R, hm.deltaL, hm.Abar, 1e100)\n",
    "        mprime = dfdA(hm.m, hm.deltaL, hm.Abar, 1e100)\n",
    "\n",
    "\n",
    "        \n",
    "        #xiprime.fill(1)\n",
    "        \n",
    "        g = hm.gamma(hm.R, hm.m, hm.U)\n",
    "        r = hm.rho(hm.R, hm.m, hm.U, hm.xi, g, xiprime, Rprime, mprime)\n",
    "        drho = hm.drho(hm.R, hm.m, hm.U, g, hm.xi, Rprime, mprime, xiprime)\n",
    "        exi = np.exp(hm.xi)\n",
    "        ephi = hm.ephi(hm.R, hm.U, g, hm.xi, xiprime, Rprime)\n",
    "        elambda = hm.elambda(ephi, exi, xiprime)\n",
    "        epsi = hm.epsi(hm.R, hm.U, g, hm.xi, r, ephi)\n",
    "        \n",
    "        plt.plot(hm.Abar[2:], ddfddA(hm.U, hm.deltaL, hm.Abar, 0)[2:])\n",
    "        #plt.plot(np.arange(hm.N), dfdA(r, hm.deltaL, hm.Abar, 1e100))\n",
    "        #plt.plot((hm_mix.Abar*hm_mix.R * np.exp((hm_mix.alpha) * hm_mix.xi)), r* (np.exp((hm_mix.alpha) * hm_mix.xi))**(-3*(1+w)))\n",
    "        #plt.plot(np.arange(hm.N), ephi)\n",
    "        #plt.plot(np.arange(hm.N), epsi)\n",
    "        #plt.plot(np.arange(hm.N),elambda / epsi)\n",
    "        #plt.xlim(0.1,3)\n",
    "        #plt.ylim(-0.1,0)\n",
    "        #plt.yscale('log')\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " epsi[0] / elambda[0] * (hm.U[1] - hm.U[0]) / ( (hm.Abar[1] - hm.Abar[0]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm=hm_mix_bak\n",
    "hm_mix = hm_mix_bak"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        xiprime = dfdA(hm.xi, hm.deltaL, hm.Abar, 1e100)\n",
    "        Rprime = dfdA(hm.R, hm.deltaL, hm.Abar, 1e100)\n",
    "        mprime = dfdA(hm.m, hm.deltaL, hm.Abar, 1e100)\n",
    "\n",
    "\n",
    "        \n",
    "        #xiprime.fill(1)\n",
    "        \n",
    "        g = hm.gamma(hm.R, hm.m, hm.U)\n",
    "        r = hm.rho(hm.R, hm.m, hm.U, hm.xi, g, xiprime, Rprime, mprime)\n",
    "        drho = hm.drho(hm.R, hm.m, hm.U, g, hm.xi, Rprime, mprime, xiprime)\n",
    "        exi = np.exp(hm.xi)\n",
    "        ephi = hm.ephi(hm.R, hm.U, g, hm.xi, xiprime, Rprime)\n",
    "        elambda = hm.elambda(ephi, exi, xiprime)\n",
    "        epsi = hm.epsi(hm.R, hm.U, g, hm.xi, r, ephi)\n",
    "        \n",
    "        #plt.plot(np.arange(hm.N-1), ddfddA(hm.U, hm.deltaL, hm.Abar, 0)[1:])\n",
    "        #plt.plot(np.arange(hm.N), dfdA(r, hm.deltaL, hm.Abar, 1e100))\n",
    "        plt.plot((hm_mix.Abar*hm_mix.R * np.exp((hm_mix.alpha) * hm_mix.xi)), r* (np.exp((hm_mix.alpha) * hm_mix.xi))**(-3*(1+w)))\n",
    "        #plt.plot(np.arange(hm.N), ephi)\n",
    "        #plt.plot(np.arange(hm.N), epsi)\n",
    "        #plt.plot(np.arange(hm.N),elambda / epsi)\n",
    "        plt.xlim(0,15)\n",
    "        plt.ylim(1e-8,1)\n",
    "        plt.yscale('log')\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        xiprime = dfdA(hm.xi, hm.deltaL, hm.Abar, 1e100)\n",
    "        Rprime = dfdA(hm.R, hm.deltaL, hm.Abar, 1e100)\n",
    "        mprime = dfdA(hm.m, hm.deltaL, hm.Abar, 1e100)\n",
    "\n",
    "\n",
    "        \n",
    "        #xiprime.fill(1)\n",
    "        \n",
    "        g = hm.gamma(hm.R, hm.m, hm.U)\n",
    "        r = hm.rho(hm.R, hm.m, hm.U, hm.xi, g, xiprime, Rprime, mprime)\n",
    "        drho = hm.drho(hm.R, hm.m, hm.U, g, hm.xi, Rprime, mprime, xiprime)\n",
    "        exi = np.exp(hm.xi)\n",
    "        ephi = hm.ephi(hm.R, hm.U, g, hm.xi, xiprime, Rprime)\n",
    "        elambda = hm.elambda(ephi, exi, xiprime)\n",
    "        epsi = hm.epsi(hm.R, hm.U, g, hm.xi, r, ephi)\n",
    "        \n",
    "        #plt.plot(np.arange(hm.N-1), ddfddA(hm.U, hm.deltaL, hm.Abar, 0)[1:])\n",
    "        #plt.plot(np.arange(hm.N), dfdA(r, hm.deltaL, hm.Abar, 1e100))\n",
    "        plt.plot((hm_mix.Abar*hm_mix.R * np.exp((hm_mix.alpha) * hm_mix.xi)), r* (np.exp((hm_mix.alpha) * hm_mix.xi))**(-3*(1+w)))\n",
    "        #plt.plot(np.arange(hm.N), ephi)\n",
    "        #plt.plot(np.arange(hm.N), epsi)\n",
    "        #plt.plot(np.arange(hm.N),elambda / epsi)\n",
    "        plt.xlim(0,25)\n",
    "        plt.ylim(1e-8,1)\n",
    "        plt.yscale('log')\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        xiprime = dfdA(hm.xi, hm.deltaL, hm.Abar, 1e100)\n",
    "        Rprime = dfdA(hm.R, hm.deltaL, hm.Abar, 1e100)\n",
    "        mprime = dfdA(hm.m, hm.deltaL, hm.Abar, 1e100)\n",
    "\n",
    "\n",
    "        \n",
    "        #xiprime.fill(1)\n",
    "        \n",
    "        g = hm.gamma(hm.R, hm.m, hm.U)\n",
    "        r = hm.rho(hm.R, hm.m, hm.U, hm.xi, g, xiprime, Rprime, mprime)\n",
    "        drho = hm.drho(hm.R, hm.m, hm.U, g, hm.xi, Rprime, mprime, xiprime)\n",
    "        exi = np.exp(hm.xi)\n",
    "        ephi = hm.ephi(hm.R, hm.U, g, hm.xi, xiprime, Rprime)\n",
    "        elambda = hm.elambda(ephi, exi, xiprime)\n",
    "        epsi = hm.epsi(hm.R, hm.U, g, hm.xi, r, ephi)\n",
    "        \n",
    "        #plt.plot(np.arange(hm.N-1), ddfddA(hm.U, hm.deltaL, hm.Abar, 0)[1:])\n",
    "        #plt.plot(np.arange(hm.N), dfdA(r, hm.deltaL, hm.Abar, 1e100))\n",
    "        plt.plot((hm.Abar*hm.R * np.exp((hm.alpha) * hm.xi)), epsi)\n",
    "        #plt.plot(np.arange(hm.N), ephi)\n",
    "        #plt.plot(np.arange(hm.N), epsi)\n",
    "        #plt.plot(np.arange(hm.N),elambda / epsi)\n",
    "        plt.xlim(0,25)\n",
    "        #plt.ylim(1e-8,1)\n",
    "        plt.yscale('log')\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_uni.u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mix.u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(hm_mix.R**2 * hm_mix.m * hm_mix.Abar**2 * np.exp(2*(hm_mix.alpha-1)*hm_mix.xi)).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(hm_uni.R**2 * hm_uni.m * hm_uni.Abar**2 * np.exp(2*(hm_uni.alpha-1)*hm_uni.xi)).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(hm_mix.Abar[:], hm_mix.U[:])\n",
    "#plt.plot(hm_bak.Abar[:], hm_bak.U[:])\n",
    "\n",
    "\n",
    "#plt.xlim(0,30)\n",
    "#plt.ylim(-10,-19.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(hm.R**2 * hm.m * hm.Abar**2 * np.exp(2*(hm.alpha-1)*hm.xi)).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((hm.Abar*hm.R * np.exp((hm.alpha) * hm.xi)), hm.R**2 * hm.m * hm.Abar**2 * np.exp(2*(hm.alpha-1)*hm.xi))\n",
    "#plt.plot(hm.Abar, hm_bak.R**2 * hm_bak.m * hm_bak.Abar**2 * np.exp(2*(hm_bak.alpha-1)*(hm_bak.xi)))\n",
    "\n",
    "plt.xlim(0,60)\n",
    "plt.ylim(0.,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        xiprime = dfdA(hm.xi, hm.deltaL, hm.Abar, 0)\n",
    "        Rprime = dfdA(hm.R, hm.deltaL, hm.Abar, 0)\n",
    "        mprime = dfdA(hm.m, hm.deltaL, hm.Abar, 0)\n",
    "        Uprime = dfdA(hm.U, hm.deltaL, hm.Abar, 0)\n",
    "        #print((xiprime ))\n",
    "        g = hm.gamma(hm.R, hm.m, hm.U)\n",
    "        r = hm.rho(hm.R, hm.m, hm.U, hm.xi, g, xiprime, Rprime, mprime)\n",
    "        p = hm.P(r)\n",
    "        exi = np.exp(hm.xi)\n",
    "        ephi = hm.ephi(hm.R, hm.U, g, xiprime, Rprime)\n",
    "        \n",
    "        epsi = hm.epsi(hm.R, hm.U, g, hm.xi, r, ephi)\n",
    "        plt.plot(hm.Abar, epsi)\n",
    "        plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(ms.Abar, ms_bak.R**2 * ms_bak.m * ms_bak.Abar**2 * np.exp(2 * (ms_bak.alpha-1) * ms_bak.xi))\n",
    "plt.plot(ms.Abar, np.array([1]*ms.N))\n",
    "#plt.xlim(0,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(Abar[:], ms.U[:])\n",
    "#plt.xlim(0,30)\n",
    "#plt.ylim(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(ms_mix.Abar[:ms_mix.N-1], ms_mix.U_hm[:ms_mix.N-1])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(ms.Abar[:ms.N-1], ms.U_hm[:ms.N-1])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm.deltau_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "ms_std = copy.copy(ms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms = copy.copy(ms_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hm = HM(ms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hm.run_steps(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "hm_bak = copy.copy(hm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "hm_bak2 = copy.copy(hm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm = copy.copy(hm_bak2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "hm_bak3 = copy.copy(hm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm.u / hm.deltau_i / 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.exp(hm.alpha * hm.xi)\n",
    "H = np.exp(-hm.xi) / hm.RH\n",
    "rho_b = a**(1+hm.w)\n",
    "Rb = a * hm.A\n",
    "plt.plot(hm.R * Rb ,  (np.exp(-hm.xi/2) * hm.R **3 * hm.Abar**3 * hm.m ) / 2 )\n",
    "plt.xlim(0,15)\n",
    "plt.ylim(0,10)\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(hm.Abar, hm.m)\n",
    "plt.plot(hm_bak.Abar, hm_bak.m)\n",
    "\n",
    "\n",
    "#plt.xlim(0,300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        xiprime = dfdA(hm.xi, hm.deltaL, hm.Abar, 0)\n",
    "        Rprime = dfdA(hm.R, hm.deltaL, hm.Abar, 0)\n",
    "        mprime = dfdA(hm.m, hm.deltaL, hm.Abar, 0)\n",
    "\n",
    "\n",
    "        \n",
    "        #xiprime.fill(1)\n",
    "        \n",
    "        g = hm.gamma(hm.R, hm.m, hm.U)\n",
    "        r = hm.rho(hm.R, hm.m, hm.U, hm.xi, g, xiprime, Rprime, mprime)\n",
    "        exi = np.exp(hm.xi)\n",
    "        ephi = hm.ephi(hm.R, hm.U, g, xiprime, Rprime)\n",
    "        elambda = hm.elambda(ephi, exi, xiprime)\n",
    "        epsi = hm.epsi(hm.R, hm.U, g, hm.xi, r, ephi)\n",
    "        plt.scatter(np.arange(hm.N), ephi)\n",
    "        plt.scatter(np.arange(hm.N), epsi)\n",
    "        plt.xlim(0,300)\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(hm.Abar,  hm.xi-0)\n",
    "\n",
    "\n",
    "#plt.xlim(0,300)\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "( hm.R**2 * hm.m * hm.Abar**2 * np.exp(2*(hm.alpha-1)*hm.xi)).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(hm.Abar, hm.R**2 * hm.m * hm.Abar**2 * np.exp(2*(hm.alpha-1)*hm.xi))\n",
    "plt.plot(hm.Abar, hm_bak.R**2 * hm_bak.m * hm_bak.Abar**2 * np.exp(2*(hm_bak.alpha-1)*(hm_bak.xi)))\n",
    "\n",
    "\n",
    "#plt.xlim(0,300)\n",
    "plt.ylim(0.,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(hm.Abar, hm.R**2 * hm.m * hm.A**2 * np.exp(2*(hm.alpha-1)*hm.xi))\n",
    "plt.plot(hm.Abar, hm_bak.R**2 * hm_bak.m * hm_bak.A**2 * np.exp(2*(hm_bak.alpha-1)*hm_bak.xi))\n",
    "\n",
    "\n",
    "#plt.xlim(0,300)\n",
    "#plt.ylim(0.915,0.93)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        xiprime = dfdA(hm.xi, hm.deltaL, hm.Abar, 0)\n",
    "        Rprime = dfdA(hm.R, hm.deltaL, hm.Abar, 0)\n",
    "        mprime = dfdA(hm.m, hm.deltaL, hm.Abar, 0)\n",
    "        Uprime = dfdA(hm.U, hm.deltaL, hm.Abar, 0)\n",
    "        #print((xiprime ))\n",
    "        g = hm.gamma(hm.R, hm.m, hm.U)\n",
    "        r = hm.rho(hm.R, hm.m, hm.U, hm.xi, g, xiprime, Rprime, mprime)\n",
    "        p = hm.P(r)\n",
    "        exi = np.exp(hm.xi)\n",
    "        ephi = hm.ephi(hm.R, hm.U, g, xiprime, Rprime)\n",
    "        \n",
    "        epsi = hm.epsi(hm.R, hm.U, g, hm.xi, r, ephi)\n",
    "        plt.plot(hm.Abar, epsi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(hm.Abar, hm.U[:])\n",
    "#plt.xlim(0,30)\n",
    "#plt.ylim(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(Abar[:n-1], ms.xi_hm[:])\n",
    "#plt.xlim(0,30)\n",
    "#plt.ylim(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(Abar[:n-1], ms_std.xi_hm[:])\n",
    "#plt.xlim(0,30)\n",
    "#plt.ylim(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}