pressure_loop_thermal_boundary_layer.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "64c675b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import interpolate\n",
    "from matplotlib.tri import Triangulation\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ef0fb949",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=np.loadtxt('./data_tomo/V_mean.txt',comments='#')\n",
    "X = np.unique(data[:,0])\n",
    "Y = np.unique(data[:,1])\n",
    "Z = np.unique(data[:,2])\n",
    "nx,ny,nz = np.meshgrid(X,Y,Z)\n",
    "V= np.zeros_like(nz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "59dca96f",
   "metadata": {},
   "outputs": [],
   "source": [
    "m,n,p = np.shape(nz)\n",
    "for i in range(m):\n",
    "    for j in range(n):\n",
    "        for k in range(p):\n",
    "            ind = np.where( (data[:,0] == nx[i][j][k]) & (data[:,1] == ny[i][j][k]) & (data[:,2] == nz[i][j][k]))\n",
    "            #print(nx[i][j][k],ny[i][j][k],nz[i][j][k])\n",
    "            try:\n",
    "                V[i][j][k] = data[ind[0],3]\n",
    "            except:\n",
    "                print(ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ebf2665f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "300.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nz[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e368d394",
   "metadata": {},
   "outputs": [],
   "source": [
    "import V2RhoT_gibbs_lib as lib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5abc2f84",
   "metadata": {},
   "outputs": [],
   "source": [
    "DMM_no_atten = np.loadtxt('./databases/DMM_HP',comments='#')\n",
    "# correction using grain size = 10 mm and oscillatio period of 75 seconds.\n",
    "# Attenuation model of Jackson and Faul 2010\n",
    "# Function: lib.atten_correction (T (oC),P (Pascal),VP (km/s),Vs (km/s),oscilation period (s), grain size (mm))\n",
    "DMM_atten_corrected = np.copy(DMM_no_atten)\n",
    "for i in range(len(DMM_atten_corrected)):\n",
    "    DMM_atten_corrected[i,3],DMM_atten_corrected[i,4] = lib.atten_correction(DMM_atten_corrected[i,0],DMM_atten_corrected[i,1]*100000,\n",
    "                                                         DMM_atten_corrected[i,3],DMM_atten_corrected[i,4],75,10)\n",
    "\n",
    "# correction for melts\n",
    "# These are relations from lab experiments. More details in Afonso et al., 2016 III\n",
    "# Function: lib.velocity_melt_correction_mantle (T (oC),P (GPa),VP (km/s),Vs (km/s),oscilation period (s), grain size (mm))\n",
    "DMM_atten_melt_corrected = np.copy(DMM_atten_corrected)\n",
    "melt = np.zeros_like(DMM_atten_melt_corrected[:,0])\n",
    "for i in range(len(DMM_atten_melt_corrected)):\n",
    "    DMM_atten_melt_corrected[i,3],DMM_atten_melt_corrected[i,4],melt[i] = lib.velocity_melt_correction_mantle(DMM_atten_melt_corrected[i,0]-273.15,\n",
    "                                                                                             DMM_atten_melt_corrected[i,1]/10000,\n",
    "                                                             DMM_atten_melt_corrected[i,3],DMM_atten_melt_corrected[i,4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "a87a3b20",
   "metadata": {},
   "outputs": [],
   "source": [
    "### lithostatic pressure from ak135\n",
    "def lithostatic_pressure(depth,density):\n",
    "    Pressure_lith = np.zeros_like(depth)\n",
    "    for j in range(len(depth)-1): \n",
    "        DepthCalculated = False\n",
    "        i = 0\n",
    "        Pressure = 0.\n",
    "        d = np.absolute(depth[j])*1000.\n",
    "        while DepthCalculated is False:\n",
    "            if i > len(depth):\n",
    "                print('Maximum depth reached.')\n",
    "                break\n",
    "            z1   = depth[i]*1000.\n",
    "            Rho1 = density[i]\n",
    "            z2   = depth[i+1]*1000.\n",
    "            Rho2 = density[i+1]\n",
    "            if d < depth[i+1]*1000.:\n",
    "                Rho2 = Rho1 + (d - z1)/(z2 - z1)*(Rho2 - Rho1)\n",
    "                dz = d - z1\n",
    "                DepthCalculated = True\n",
    "            else:\n",
    "                dz = z2 - z1\n",
    "                Pressure = Pressure + dz*(Rho1 + (Rho2 - Rho1)/2)*9.8\n",
    "                i += 1\n",
    "        Pressure_lith[j]=Pressure*1e-5 # in bar\n",
    "        print(Pressure_lith[j],depth[j],density[j])\n",
    "    p_func = interpolate.interp1d(depth,Pressure_lith)\n",
    "    return p_func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5cae951d",
   "metadata": {},
   "outputs": [],
   "source": [
    "### lithostatic pressure \n",
    "def lithostatic_pressure(depth,density):\n",
    "    Pressure_lith = np.zeros_like(depth)\n",
    "    # Check is the first index is at zero depth or not\n",
    "    # If it is not then pressures is set equal rho*g*h\n",
    "    # else is it set to 0\n",
    "    if depth[0] != 0:\n",
    "        Pressure_lith[0] = density[0] * depth[0]* 9.8 * 1e3 * 1e-5\n",
    "    else:\n",
    "        #thickness = (depth[1]-depth[0]) #*1e3\n",
    "        Pressure_lith[0] = 0\n",
    "    # Now looping through depths\n",
    "    for j in range(len(depth)-2):\n",
    "        thickness = (depth[j+2] - depth [j+1]) *1e3\n",
    "        Pressure_lith[j+1]=Pressure_lith[j] + thickness*density[j+1] *9.8*1e-5\n",
    "        \n",
    "    # fixing last index\n",
    "    Pressure_lith[-1] = Pressure_lith[-2] + (depth[-1] - depth [-2])*1e3 * density[-2]*9.8*1e-5 \n",
    "    #print(depth[:],density[:],Pressure_lith[:])\n",
    "\n",
    "    p_func = interpolate.interp1d(depth,Pressure_lith)\n",
    "    return p_func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "2d3f370f",
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh_hold=1\n",
    "def vel_to_temp_P_loop(nx,ny,nz,Vs,Table):\n",
    "    Temperature_out = np.zeros_like(nx)\n",
    "    Density_out     = np.zeros_like(nx)#+3300\n",
    "    diff_Vs         = np.zeros_like(nx)\n",
    "    P_out           = np.zeros_like(nx)\n",
    "    melt_out        = np.zeros_like(nx)\n",
    "    #Vp_out          = []#np.zeros_like(tomo[:,1])\n",
    "    #Vs_out          = []#np.zeros_like(tomo[:,1])\n",
    "    m,n,p=np.shape(nx)\n",
    "    # looping over all the geospatial nodes\n",
    "    for i in range(1):\n",
    "        for j in range(1):\n",
    "            #P_func=lithostatic_pressure(nz[i][j][:],Density_out[i][j][:])\n",
    "            P_func=lib.pressure_inter\n",
    "            threshhold = 10\n",
    "            p_counter  = 0\n",
    "            while threshhold > 0:\n",
    "                #print(f'{p_counter} Pressure counter {threshhold}.')\n",
    "                p_counter = p_counter+1\n",
    "                for k in range(p):\n",
    "                    P  = P_func(nz[i][j][k])\n",
    "                    Vs_in = Vs[i][j][k]\n",
    "                    P_table,temp,dens,vp,vs=lib.lookup_vs_P_accurate(Vs_in,P.tolist(),Table)\n",
    "                    Temperature_out[i][j][k]=temp\n",
    "                    Density_out[i][j][k]=dens\n",
    "                P_func_new=lithostatic_pressure(nz[i][j][:],Density_out[i][j][:])\n",
    "                threshhold = np.sqrt(((P_func_new(nz[i][j][:]) - P_func(nz[i][j][:])) ** 2).mean())\n",
    "                #threshhold = P_func_new(nz[i][j][-1]) - P_func(nz[i][j][-1])\n",
    "                \n",
    "                #thresh_hold = abs(P_func_new(nz[i][j][:]) - P_func(nz[i][j][-1]))\n",
    "                print(f'{p_counter} Pressure updated {threshhold}.')\n",
    "                #print(threshhold)\n",
    "                P_func = P_func_new\n",
    "                #Temperature_out[i][j][k]=temp\n",
    "                #Temperature_out[i][j][k]=temp\n",
    "    return Temperature_out,Density_out \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "16ae4ee9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 Pressure updated 913.6671263240239.\n",
      "2 Pressure updated 79.69359081098095.\n",
      "3 Pressure updated 8.501500000010129.\n",
      "4 Pressure updated 0.0.\n"
     ]
    }
   ],
   "source": [
    "T_P_loop,D_P_loop=vel_to_temp_P_loop(nx,ny,nz,V,DMM_atten_melt_corrected)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "a6e37fb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vel_to_temp(nx,ny,nz,Vs,Table):\n",
    "    Temperature_out = np.zeros_like(nx)\n",
    "    Density_out     = np.zeros_like(nx) #+3300\n",
    "    diff_Vs         = np.zeros_like(nx)\n",
    "    P_out           = np.zeros_like(nx)\n",
    "    melt_out        = np.zeros_like(nx)\n",
    "    #Vp_out          = []#np.zeros_like(tomo[:,1])\n",
    "    #Vs_out          = []#np.zeros_like(tomo[:,1])\n",
    "    m,n,p=np.shape(nx)\n",
    "    # looping over all the geospatial nodes\n",
    "    for i in range(1):\n",
    "        for j in range(1):\n",
    "            #P_func=lithostatic_pressure(nz[i][j][:],Density_out[i][j][:])\n",
    "            P_func=lib.pressure_inter\n",
    "            for k in range(p):\n",
    "                    P  = P_func(nz[i][j][k])\n",
    "                    print(P)\n",
    "                    Vs_in = Vs[i][j][k]\n",
    "                    print(Vs_in)\n",
    "                    P_table,temp,dens,vp,vs=lib.lookup_vs_P_accurate(Vs_in,P.tolist(),Table)\n",
    "                    #print(dens)\n",
    "                    Temperature_out[i][j][k]=temp\n",
    "                    Density_out[i][j][k]=dens\n",
    "                    print(Density_out[i][j][k])\n",
    "    return Temperature_out,Density_out "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "5a05d4d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65324.44800000001\n",
      "4.4008613\n",
      "3385.64\n",
      "82301.478\n",
      "4.5544051\n",
      "3450.565\n",
      "100055.6088\n",
      "4.6854635\n",
      "3525.44\n",
      "118315.36080000002\n",
      "4.7774704\n",
      "3591.6400000000003\n",
      "137120.1104\n",
      "4.8259661\n",
      "3642.66\n"
     ]
    }
   ],
   "source": [
    "T,D=vel_to_temp(nx,ny,nz,V,DMM_atten_melt_corrected)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "ddab981c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2e241782188>]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgcUlEQVR4nO3deXzU9b3v8deHbOx7wpIBwr6JQIKIu7IIahUltNfe01atvba1nrq0Lj09j3t6zr33cdWeqj1dTmtbezynttpLUNQKCAhaVyQBZAsQ9kkISYAEkpD9e/+YQSMNQyAz85uZvJ+PBw9ny3w//HDemfnN/N5jzjlERCQxdfF6ABERiRyFvIhIAlPIi4gkMIW8iEgCU8iLiCSwZK8HaG3gwIEuKyvL6zFEROJKfn5+hXMuva3rYirks7Ky2LBhg9djiIjEFTM7cLbrtLtGRCSBhS3kzSzJzDaa2evB8yPN7CMzKzKzl8wsNVxriYhI+4Tzmfz9wI5W558AnnbOjQGOA3eHcS0REWmHsIS8mfmAm4DfBs8bMBtYErzJ88Ct4VhLRETaL1zP5J8BHgFagucHAJXOuabgeT+Q2dYPmtk9ZrbBzDaUl5eHaRwREYEwhLyZfQEoc87lX8jPO+eedc7NcM7NSE9v8xNAIiJygcLxEcorgFvM7EagK9Ab+CnQ18ySg8/mfUBxGNYSEZHz0OFn8s65HzjnfM65LOB24C3n3N8Ba4HFwZvdASzr6FoiIonomdW7yD9wPCL3HcnPyT8KPGRmRQT20f8ugmuJiMSlgoPHeWb1bt7dXRGR+w/rEa/OuXXAuuDpvcDMcN6/iEgicc7x+PJCBvZM5RtXjYzIGjriVUTEI+t2lrN+3zHunzOWHmmRaZlRyIuIeKC5xfHEikKyBnTn9pnDI7aOQl5ExAOvbCymsPQk358/npSkyEWxQl5EJMrqGpt5atUuLvb14caLhkR0LYW8iEiU/eHDAxRXnuKxBRPo0sUiupZCXkQkiqpONfLztUVcPS6dy8cMjPh6CnkRkSj69dt7qKxt5NEF46OynkJeRCRKSqvqeO69fdw6bSiTh/aJypoKeRGRKPnpml00tzi+d310nsWDQl5EJCqKyqp56eNDfGXWCIb17x61dRXyIiJR8OOVhXRPTea+68ZEdV2FvIhIhOUfOM7KbUe45+pRDOiZFtW1FfIiIhHknOOJ5YUM7JnG3VdGpoQsFIW8iEgEvVVYxvr9x7h/buRKyEJRyIuIRMjnSsguGebJDAp5EZEIeXljMbuOVPPw/AkRLSELRSEvIhIBdY3NPPXmTqb6+nDjlMGezaGQFxGJgP/64AAlVXU8esMEzCJbQhaKQl5EJMxOl5BdMy6dy0dHvoQsFIW8iEiY/ertPZyoa+TRBRO8HkUhLyISTqVVdTz37j5unZbJpKG9vR5HIS8iEk7PrN6Fc/DQvHFejwIo5EVEwqao7CR/3hD9ErJQFPIiImHy5IqdgRKy2dEtIQtFIS8iEgb5B47x5vYjfOuaUfTvker1OJ9SyIuIdJBzjseXF5LeK42ve1BCFopCXkSkg9bsKOPj/cd5YO5YuqdGv4QslA6HvJl1NbP1ZrbZzLaZ2T8HLx9pZh+ZWZGZvWRmsfP6RUQkTE6XkI0a2IMvzfCmhCyUcDyTrwdmO+emAtOABWY2C3gCeNo5NwY4DtwdhrVERGJKXoGf3WXVPDx/vGclZKF0eCIXUB08mxL844DZwJLg5c8Dt3Z0LRGRWFLX2MzTq3YxdVhfFlzkXQlZKGH5tWNmSWa2CSgDVgF7gErnXFPwJn4gMxxriYjEiuff38/hqjp+4HEJWShhCXnnXLNzbhrgA2YC7S5sMLN7zGyDmW0oLy8PxzgiIhFXVdvIL9YWcd34dGaNGuD1OGcV1h1IzrlKYC1wGdDXzE6/zewDis/yM88652Y452akp6eHcxwRkYj55dtFnKxv4pEYKCELJRyfrkk3s77B092AecAOAmG/OHizO4BlHV1LRCQWlFSe4vfv7ee2aZlMHOJ9CVko4fhA5xDgeTNLIvBL48/OudfNbDvwopn9b2Aj8LswrCUi4rlnVu8CBw/GSAlZKB0OeefcJ8D0Ni7fS2D/vIhIwth15CRL8v3cdcXImCkhCyX2PtQpIhLDnlyxkx6pyXznutgpIQtFIS8i0k4f7z/G6h1H+Na1o2OqhCwUhbyISDs453hieSEZvdK464osr8dpN4W8iEg7rN5RxoYDx3lg7riYKyELRSEvInIOTc0tPPlpCZnP63HOi0JeROQclhYUs7usmkcWjCc5BkvIQomvaUVEoqyusZmnVu1i2rC+zJ8cmyVkoSjkRURC+I/391N6oo7HYriELBSFvIjIWVTWNvDLtUXMnpAR0yVkoSjkRUTO4t/X7QmWkI33epQLppAXEWlDSeUpfv/+fhZN9zFhcGyXkIWikBcRacPTq3YB8ND1sV9CFopCXkTkDDtLT5JX4OeOy0aQ2beb1+N0iEJeROQMP15ZSI+0ZO69Nj5KyEJRyIuItLJ+3zFW7yjj29eOpl+clJCFopAXEQlyzvH48h0M6p3GXZeP9HqcsFDIi4gEvbn9CAUHK3lw7ji6pSZ5PU5YKORFRPishGx0eg8W58RXCVkoCnkREWBJvp895TU8PH9C3JWQhZI4fxMRkQt0qqGZp1fvYvrwvsyfPMjrccJKIS8ind7v39/HkRP1PLYgPkvIQlHIi0indrymgX9ft4c5EzK4NE5LyEJRyItIp/bLdUVU1zfxyIIJXo8SEQp5Eem0/Mdref79A+Rm+xg/uJfX40SEQl5EOq2nV+0GgwfnxXcJWSgKeRHplApLT7B0o587L8+K+xKyUBTyItIp/XjFTnqmJXPvtaO9HiWiFPIi0ul8tPcoawrLuPfaMfTtHv8lZKF0OOTNbJiZrTWz7Wa2zczuD17e38xWmdnu4H/7dXxcEZGOcc7x+IpCBvfuyl1XZHk9TsSF45l8E/A959wkYBbwHTObBDwGrHHOjQXWBM+LiHhq5bYjbDxYyYPzxtI1JTFKyELpcMg75w475wqCp08CO4BMYCHwfPBmzwO3dnQtEZGOaGpu4cmVhYzJ6EluduKUkIUS1n3yZpYFTAc+AgY55w4HryoF2iyEMLN7zGyDmW0oLy8P5zgiIp/z//L97C2v4ZH54xOqhCyUsP0tzawnkAc84Jw70fo655wDXFs/55x71jk3wzk3Iz09PVzjiIh8zqmGZp5etYucEf2YNymxSshCCUvIm1kKgYB/wTm3NHjxETMbErx+CFAWjrVERC7Ec+/to+xkPY/dkHglZKGE49M1BvwO2OGce6rVVa8CdwRP3wEs6+haIiIX4nhNA79at4e5EwdxSVZ/r8eJquQw3McVwFeBLWa2KXjZPwCPA382s7uBA8CXwrCWiMh5+8XaImoamnhkwXivR4m6Doe8c+5d4GyvfeZ09P5FRDrCf7yW//zgAItzfIwblJglZKF0jreXRaTTemrVLszggbmJW0IWikJeRBLWjsMneHljMXdekcXQBC4hC0UhLyIJ68kVhfRKS+bea8Z4PYpnFPIikpA+2HOUtTvLufe6MfTpnuL1OJ5RyItIwmldQnbn5Vlej+MphbyIJJwVW0vZfKiSh+aN6xQlZKEo5EUkoTQ2t/DjlTsZm9GTRdmZXo/jOYW8iCSUP284xN6KGh5ZMKHTlJCFoi0gIgmjtqGJZ1bvZsaIfsydmOH1ODFBIS8iCeO5d/dR3glLyEJRyItIQjhW08Cv3t7LvEmDmNHJSshCUciLSEL4+VtF1DY08cj8zldCFopCXkTi3qFjtfzhwwN8MWcYYzthCVkoCnkRiXtPny4hmzfW61FijkJeROLa9pITvLypmLuuGMmQPp2zhCwUhbyIxLUnVxbSu2sK375mtNejxCSFvIjErff3VLBuZznfuW50py4hC0UhLyJxyTnHE8sLGdqnK1+7LMvrcWKWQl5E4tLyraVs9lfxoErIQlLIi0jcOV1CNm5QTxZl+7weJ6Yp5EUk7rz08SH2VdTw6IIJJHVRfUEoCnkRiSs19YESsplZ/Zk9QSVk56KQF5G48ty7+6iorudRlZC1i0JeROLG0ep6fv3OXuZPHkTOiH5ejxMXFPIiEjd+vjZQQvawSsjaTSEvInHhdAnZl2YMY0yGSsjaSyEvInHhJ2/upIsZD8wd5/UocSUsIW9mz5lZmZltbXVZfzNbZWa7g//VDjQRuSBbi6t4ZVMJX79yJIP7dPV6nLgSrmfy/wEsOOOyx4A1zrmxwJrgeRGR8/bkyp306ZbCt1RCdt7CEvLOuXeAY2dcvBB4Pnj6eeDWcKwlIp3Le0UVvLOrnPuuG0OfbiohO1+R3Cc/yDl3OHi6FBjU1o3M7B4z22BmG8rLyyM4jojEm5YWx+PBErKvXjbC63HiUlTeeHXOOcCd5bpnnXMznHMz0tPTozGOiMSJN7YeZktxFQ9dP14lZBcokiF/xMyGAAT/WxbBtUQkwZwuIRs/qBe3Tc/0epy4FcmQfxW4I3j6DmBZBNcSkQTz4vqDHDhay6M3jFcJWQeE6yOUfwI+AMabmd/M7gYeB+aZ2W5gbvC8iMg51dQ38dM1u5k5sj/XjVcJWUckh+NOnHNfPstVc8Jx/yLSufzu3X1UVDfw7NdUQtZROuJVRGLK0ep6fv32HhZMHkz2cB1D2VEKeRGJKT97q4i6phYeXqASsnBQyItIzDh4tJYXPgqUkI1O7+n1OAlBIS8iMeMnq3aS1MV4YO5Yr0dJGAp5EYkJW4urWLaphLuvHMmg3iohCxeFvIjEhCdWFNK3ewrfVAlZWCnkRcRz7+6u4K+7K7jvujH07qoSsnBSyIuIp1paHE+sKCSzbzeVkEWAQl5EPPWXLYESsu9dP460ZJWQhZtCXkQ809DUwr++uZMJg3uxcJpKyCJBIS8innnx42AJ2YIJKiGLEIW8iHiiur6Jf1uzm0tH9ufa8fouiUhRyIuIJ377171UVDfw2A0qIYskhbyIRF35yXp+885ebrhoMNNVQhZRCnkRibqfv7WbuqYWvj9fJWSRppAXkag6cLSGFz46yH+7RCVk0aCQF5Go+tc3d5GS1IUH5qiELBoU8iISNVv8Vby2OVBClqESsqhQyItI1DyxopB+3VO455pRXo/SaSjkRSQq/rq7nHeLKrhv9liVkEWRQl5EIq6lxfH48kAJ2VdmDfd6nE5FIS8iEffaJyVsKznB9+erhCzaFPIiElENTS385M1dTBzSm4VTVUIWbQp5EYmoP60/yMFjtTy6YDxdVEIWdQp5EYmY0yVkl40awDXjVELmBYW8iETMb97Zy9EalZB5SSEvIhFRfrKe3/x1LzdNGcLUYX29HqfTinjIm9kCM9tpZkVm9lik1xOR2PCzt3ZTrxIyz0U05M0sCfgFcAMwCfiymU2K5Joi4r39FTX88aODfHnmMEYO7OH1OJ1apJ/JzwSKnHN7nXMNwIvAwgivKSIeampu4YkVhaQkdeG7KiHzXHKE7z8TONTqvB+4tPUNzOwe4B6A4cN1JJxIvNpx+AR5+X5e2VRCRXU9988ZS0YvlZB5LdIhf07OuWeBZwFmzJjhPB5HRM5DRXU9yzaVkJfvZ/vhEyR3MWZPyCA3x8e8iYO8Hk+IfMgXA8NanfcFLxOROFXf1MxbO8rIK/Czbmc5TS2OKZl9+NHNk7hlWib9e6R6PaK0EumQ/xgYa2YjCYT77cB/j/CaIhJmzjk2+6vIy/fz2iclVNY2ktErjbuvHElujo9xg3p5PaKcRURD3jnXZGb3ASuBJOA559y2SK4pIuFTWlXH0o1+lhYUU1RWTVpyF66fPJjc7EyuHDOQ5CQdahPrIr5P3jn3BvBGpNcRkfA41dDMm9tLWZLv592iCpyDGSP68X8XTeHGKUPo001d8PHE8zdeRcR7zjk+3n+cvHw/f9lymOr6JjL7duPvrxvDomwfWfqse9xSyIt0YoeO1ZJXENgdc/BYLd1Tk7hxyhBys31cOrK/WiMTgEJepJM5WdfI8i2lLCnws37fMczg8tEDeGDuWOZPHkyPNMVCItG/pkgn0NzieH9PBXn5flZsK6WusYVRA3vw8Pzx3Do9k8y+3bweUSJEIS+SwIrKqskr8PNyQTGlJ+ro3TWZRdk+Fuf4mD6sr+p/OwGFvEiCqaxt4LXNJSwpKGbzoUqSuhhXjx3IP35hInMnDqJrir5jtTNRyIskgMbmFt7eWU5egZ81O8poaG5hwuBe/PDGiSycPlQdMp2YQl4kjm0rqSIvv5hXNxdTUd3AgB6p/N2s4eRm+5g8tLd2x4hCXiTelJ+sZ9mmYpbk+yksPUlKkjFnwiByc3xcOz6dFB2FKq0o5EXiQF1jM2uCpWBv7yqnucUx1deHf1k4mZsvHko/lYLJWSjkRWKUc46NhyrJy/fz+ieHqTrVyKDeafyPq0axOCeTMRkqBZNzU8iLxJjDVadYWlBMXoGfveU1dE3pwvzJg8nN9nHFmIEk6ShUOQ8KeZEYUNvQxMptpeTlF/PenkAp2Mys/nzz6lHcOGUIvbqqFEwujEJexCMtLY71+4+Rl+/njS2HqWloZlj/bnx39lhys30MH9Dd6xElASjkRaLswNEa8gqKWVrgx3/8FD3Tkrnp4kAp2CVZKgWT8FLIi0TBibpG3vjkMHkFfj7efxwzuHLMQL5//XjmTx5Mt1QdhSqRoZAXiZDmFse7RYFSsJXbSqlvamF0eg8eWTCe26ZnMqSPSsEk8hTyImG2+8hJlhT4eWVjMUdO1NOnWwpfnOEjN9vHNJWCSZQp5EXC4HhNA69uLiGvwM8n/iqSuhjXjkvnn272MWdiBmnJ2h0j3lDIi1ygxuYW1hYGjkJ9q7CMxmbHxCG9+cebJrJwWibpvdK8HlFEIS9yPpxzbCs5wZJ8P69uLuFYTQMDe6bytcuyyM32MWlob69HFPmchAj5hqYWjtbU640siZiyk3Us2xjYHVNYepLUpC7MnZRBbraPq8epFExiV0KE/FuFR/j2CwVcMXoguTmZzJ88mO6pCfFXEw/VNTazescR8vL9vLO7guYWx7Rhfflft17EzRcPoW93lYJJ7EuIJJzi68vfzx7L0gI/D760mR6pWwPfOJ/jY6YOLpHz4Jyj4GAleQV+Xt9cwom6Job06co3rx7FomwfYzJ6ej2iyHkx55zXM3xqxowZbsOGDRf8820dJu7r141F2T5yszMZMaBHGKeVRFJceYqXC/zkFRSzr6KGbilJLLgoUAp22egBKgWTmGZm+c65GW1el0gh31pbhU+XZPUjN9vHjRcPobcKnzq9mvomVmwtJa/Azwd7j+IcXDqyP7k5Pm6cMoSeaQnxQlc6gU4Z8q2VVJ7i5Y2fVbemJQerW3N8XKnq1k6lpcXx4d6jLCnws2JrKbUNzYwY0J1F030sys5kWH+Vgkn86fQhf5pzjk2HAvtbX9v82Zcw3Do9k8XZPsYO0pcwJKp9FTXk5ft5eWMxxZWn6HW6FCzHx4wR/XQUqsS1iIW8mX0R+BEwEZjpnNvQ6rofAHcDzcB3nXMrz3V/kQ751uqbgl+nlu9nXauvU8vN8enr1BJE1alGXv+khLx8PwUHK+licOXYdHKzA5/A6pqio1AlMUQy5CcCLcCvge+fDnkzmwT8CZgJDAVWA+Occ82h7i+aId/a6S9GzisoZsfhE/pi5DjW1NzCX3dXsKTAz6rtR2hoamFsRk9yc3zcNj2TQb27ej2iSNiFCvkOvbPknNsRXODMqxYCLzrn6oF9ZlZEIPA/6Mh6kZLeK41vXDWKb1w1iu0lJ8gr8LNsUzErtpUyoEcqt0wbSm62j8lDe+tlfYwqLD1BXr6fVzaVUH6ynn7dU/jyJcPIzfExJbOP/t2k04rUxwcygQ9bnfcHL/sbZnYPcA/A8OHDIzRO+00a2ptJQyfx2A0TeGdXOXkFfl748CC/f28/Ewb3Ijfbx8LpQ8nopWeEXjtaXc+yTYGjULeVnCC5izF7QgaLsn3MnpBBarJegYmcM+TNbDUwuI2rfuicW9bRAZxzzwLPQmB3TUfvL1xSkrowZ+Ig5kwcRGVtA699cpi8fD//540dPL6ikKvHDiQ3x8fciYO0bzeK6puaWVtYxpL8YtbtLKOpxTElsw8/unkSt0zLpL/eSxH5nHOGvHNu7gXcbzEwrNV5X/CyuNS3eypfnTWCr84aQVFZNUsL/CwtKOa+P26kd9dkvjA1sDsne7i6wiPBOccn/iryCgKlYJW1jWT0SuPuK0eyKNvH+MH6VJTI2YTlI5Rmto7Pv/E6Gfgjn73xugYYG6tvvF6I5hbH+3sC3/qzYlspdY0tjBrYg0XZmdyW7SOzr8rSOqq0qu7T4xuKyqpJS+7C9ZMHk5udyZVjBpKsN8RFgMh+uuY24GdAOlAJbHLOzQ9e90Pg60AT8IBzbvm57i+eQr61k3WNLN9SypICP+v3HcMMLhs1gNxsHzdMUVna+TjV0Myb20tZku/nvaIKWhzMGNHv06NQ+3TTkcoiZ9LBUFF06FgtecHdOQeP1dI9NYkbLhpCbk4ms0YOUFlaG5xzfLz/OHn5fv6y5TDV9U1k9u1GbnYmi7J9ZA1U55BIKAp5D5wtuBZlZ5Kr4ALg4NHgL8SNfg4dO0X31KRAe2i2j0tHqj1UpL0U8h5raxdEzohAWdpNF3euXRAn6xp5Y8th8vKLWb8/sGvr8tGBXVsLLtKuLZELoZCPIW29mThvUuDo2qsS9M3E5hbHe0UV5BX4WdnqTercHB+3Ts/Um9QiHaSQj0Fn+1jgbdMzyc3xMS4BytKKyk6yJL+YVzYWU3qijt5dk7l56lByc3xMH6aPm4qEi0I+xp3tAJ/c7My4O8DneE0DrwVLwTb7q0jqYlwzLp3cbB9zJmbowDGRCFDIx5GK6npebXWofkqScd34DHJzfFw3PjYP1W9sbmHdznLy8v2sKTxCY7NjwuBeLM7xccs0VUCIRJpCPk6dLt16eWMJFdX19O+Ryi3Bo2svyvS2LM05x7Zgmdurm0o4WtPAgB6pLJyWSW5OJpOH9vFsNpHORiEf59qqzx03qCe52YH63Iwo1ueWnaxj2cbAK43C0pOkJnVhzsQMcrN9XKNaZhFPKOQTSFVtI699UsLSgs++COOqsenk5vi4flJkytLqGptZveMIefl+3tldEfiClWF9WZydyc1Th9K3e/y8ZyCSiBTyCWpveTVLC4pZWuCnpKqOXl2T+cLFgYOJcjr4lXbOOQoOBr4q8fXNJZyoa2Jw767cFjyYa0xGzzD+TUSkIxTyCe7TL6fO97N8aymnGpvJGtCdRdmBL6f29Wv/l1MXV57i5WAtw96KGrqmdGFB8EvPLx+tLz0XiUUK+U6kur6J5VsOk1fg58O9xwCYNao/udmBgq8eaX97RGlNfRMrtpaSV+Dng71HcQ5mjuzP4mDBWq+uneeIXJF4pJDvpA4dq+XljYHdOfuP1tItJYkbLgo8K790ZH/W7zvGkgI/K7aWUtvQzPD+3T/t1hnWv/3P/kXEWwr5Ts45R/6B48H964c5Wd9EWnIX6pta6JWWzE0XDyE3x8eMDu7HFxFvKOTlU3WNzby5/Qgf7Klg1qgBXD9pMN1SdRSqSDwLFfKq/OtkuqYkccvUodwydajXo4hIFOjIFRGRBKaQFxFJYAp5EZEEppAXEUlgCnkRkQSmkBcRSWAKeRGRBKaQFxFJYDF1xKuZlQMHLvDHBwIVYRwnXGJ1Lojd2TTX+dFc5ycR5xrhnEtv64qYCvmOMLMNZzus10uxOhfE7mya6/xorvPT2ebS7hoRkQSmkBcRSWCJFPLPej3AWcTqXBC7s2mu86O5zk+nmith9smLiMjfSqRn8iIicgaFvIhIAou7kDezBWa208yKzOyxNq5PM7OXgtd/ZGZZMTLXnWZWbmabgn++EaW5njOzMjPbepbrzcz+LTj3J2aWHSNzXWtmVa221/+MwkzDzGytmW03s21mdn8bt4n69mrnXFHfXsF1u5rZejPbHJztn9u4TdQfk+2cy6vHZJKZbTSz19u4LvzbyjkXN3+AJGAPMApIBTYDk864zb3Ar4KnbwdeipG57gR+7sE2uxrIBrae5fobgeWAAbOAj2JkrmuB16O8rYYA2cHTvYBdbfw7Rn17tXOuqG+v4LoG9AyeTgE+AmadcRsvHpPtmcurx+RDwB/b+veKxLaKt2fyM4Ei59xe51wD8CKw8IzbLASeD55eAsyxyH87dXvm8oRz7h3gWIibLAT+0wV8CPQ1syExMFfUOecOO+cKgqdPAjuAzDNuFvXt1c65PBHcDtXBsynBP2d+miPqj8l2zhV1ZuYDbgJ+e5abhH1bxVvIZwKHWp3387f/s396G+dcE1AFDIiBuQBygy/xl5jZsAjP1F7tnd0LlwVfbi83s8nRXDj4Mnk6gWeArXm6vULMBR5tr+Duh01AGbDKOXfWbRbFx2R75oLoPyafAR4BWs5yfdi3VbyFfDx7Dchyzl0MrOKz39bStgICfRxTgZ8Br0RrYTPrCeQBDzjnTkRr3XM5x1yebS/nXLNzbhrgA2aa2UXRWjuUdswV1cekmX0BKHPO5UdynTPFW8gXA61/2/qCl7V5GzNLBvoAR72eyzl31DlXHzz7WyAnwjO1V3u2adQ5506cfrntnHsDSDGzgZFe18xSCATpC865pW3cxJPtda65vNpeZ8xQCawFFpxxlRePyXPO5cFj8grgFjPbT2CX7mwz+8MZtwn7toq3kP8YGGtmI80slcAbE6+ecZtXgTuCpxcDb7nguxheznXGfttbCOxXjQWvAl8LfmpkFlDlnDvs9VBmNvj0vkgzm0ng/9WIBkNwvd8BO5xzT53lZlHfXu2Zy4vtFVwr3cz6Bk93A+YBhWfcLOqPyfbMFe3HpHPuB845n3Mui0BGvOWc+8oZNwv7tkruyA9Hm3OuyczuA1YS+ETLc865bWb2L8AG59yrBB4M/2VmRQTe2Ls9Rub6rpndAjQF57oz0nMBmNmfCHzyYqCZ+YF/IvAmFM65XwFvEPjESBFQC9wVI3MtBr5tZk3AKeD2KPyyvgL4KrAluC8X4B+A4a3m8mJ7tWcuL7YXBD7587yZJRH4xfJn59zrXj8m2zmXJ4/JM0V6W6nWQEQkgcXb7hoRETkPCnkRkQSmkBcRSWAKeRGRBKaQFxFJYAp5EZEEppAXEUlg/x9bikqpA22x9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(T[0][0][:]-T_P_loop[0][0][:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "277a0e81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3389.165, 3453.32 , 3526.99 , 3591.64 , 3647.18 ])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_P_loop[0][0][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "cbae8729",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 65324.448   82301.478  100055.6088 118315.3608 137120.1104] [200. 250. 300. 350. 400.]\n",
      "[ 66427.634  83348.902 100631.153 118230.189 135829.225] [200. 250. 300. 350. 400.]\n"
     ]
    }
   ],
   "source": [
    "P=lithostatic_pressure(nz[0][0][:],D[0][0][:])\n",
    "print(lib.pressure_inter(nz[0][0][:]),nz[0][0][:])\n",
    "print(P(nz[0][0][:]),nz[0][0][:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "345335d6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([66427.634])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_new([200])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "514e9054",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "A value in x_new is above the interpolation range.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-20-015230ef9ae4>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp_ini\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mD\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mp_new\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mD\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mc:\\users\\kumar\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\scipy\\interpolate\\polyint.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, x)\u001b[0m\n\u001b[0;32m     76\u001b[0m         \"\"\"\n\u001b[0;32m     77\u001b[0m         \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx_shape\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_prepare_x\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 78\u001b[1;33m         \u001b[0my\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_evaluate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     79\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_finish_y\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx_shape\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     80\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\users\\kumar\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\scipy\\interpolate\\interpolate.py\u001b[0m in \u001b[0;36m_evaluate\u001b[1;34m(self, x_new)\u001b[0m\n\u001b[0;32m    675\u001b[0m         \u001b[0my_new\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_call\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx_new\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    676\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_extrapolate\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 677\u001b[1;33m             \u001b[0mbelow_bounds\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mabove_bounds\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_check_bounds\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx_new\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    678\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my_new\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    679\u001b[0m                 \u001b[1;31m# Note fill_value must be broadcast up to the proper size\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\users\\kumar\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\scipy\\interpolate\\interpolate.py\u001b[0m in \u001b[0;36m_check_bounds\u001b[1;34m(self, x_new)\u001b[0m\n\u001b[0;32m    707\u001b[0m                              \"range.\")\n\u001b[0;32m    708\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbounds_error\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mabove_bounds\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0many\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 709\u001b[1;33m             raise ValueError(\"A value in x_new is above the interpolation \"\n\u001b[0m\u001b[0;32m    710\u001b[0m                              \"range.\")\n\u001b[0;32m    711\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: A value in x_new is above the interpolation range."
     ]
    }
   ],
   "source": [
    "plt.plot(p_ini(D[0][0][:])-p_new(D[0][0][1:-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "70d17ca1",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=np.loadtxt('./data_tomo/NA_age_vel_stack.dat',comments='#')\n",
    "X = np.unique(data[:,0])\n",
    "#Y = np.unique(data[:,1])\n",
    "Z = np.unique(data[:,1])\n",
    "nx,nz = np.meshgrid(X,Z)\n",
    "V= np.zeros_like(nz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "af228b0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "m,n = np.shape(nz)\n",
    "for i in range(m):\n",
    "    for k in range(p):\n",
    "        ind = np.where( (data[:,0] == nx[i][k]) & (data[:,1] == nz[i][k]))\n",
    "        #print(nx[i][j][k],ny[i][j][k],nz[i][j][k])\n",
    "        #try:\n",
    "        V[i][k] = data[ind[0],2]\n",
    "        #except:\n",
    "        #    print(ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "cdebd632",
   "metadata": {},
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "Error in qhull Delaunay triangulation calculation: singular input data (exitcode=2); use python verbose option (-v) to see original qhull error.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-27-c39bacff7070>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mtri\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTriangulation\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mnz\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msharex\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msharey\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mm_v_in\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtricontourf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtri\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mV\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcmap\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'RdBu'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlevels\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m3.8\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m5\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0.05\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_ylim\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m400\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_xlabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Age (Ma)'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\users\\kumar\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\matplotlib\\tri\\triangulation.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, x, y, triangles, mask)\u001b[0m\n\u001b[0;32m     52\u001b[0m             \u001b[1;31m# No triangulation specified, so use matplotlib._qhull to obtain\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     53\u001b[0m             \u001b[1;31m# Delaunay triangulation.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 54\u001b[1;33m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtriangles\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_neighbors\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_qhull\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdelaunay\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     55\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mis_delaunay\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     56\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mRuntimeError\u001b[0m: Error in qhull Delaunay triangulation calculation: singular input data (exitcode=2); use python verbose option (-v) to see original qhull error."
     ]
    }
   ],
   "source": [
    "tri = Triangulation(nx[1][:],-nz[1][:])\n",
    "fig, ax = plt.subplots(1, 1, sharex=True, sharey=True,figsize=(10,5))\n",
    "m_v_in = ax.tricontourf(tri, V[1][:],cmap='RdBu',levels=np.arange(3.8,5,0.05))\n",
    "ax.set_ylim([-400,0])\n",
    "ax.set_xlabel('Age (Ma)')\n",
    "ax.set_ylabel('Depth (km)')\n",
    "fig.colorbar(m_v_in, ax=ax, label='CSEM Vs(km/s)', orientation='vertical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9e43a6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "eb8c4ce6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAADqCAYAAAB3L3NsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABSAUlEQVR4nO29eXRb9Zk+/lxJlmzLlmVZ3u3YcRzHjh3biZOQaSlLWBq2JJQ1pQVKCxQGDvD9lZYWmIHpsLSdztADLaUMDG0pdMjSsoTSAUrYSYjBa7zF+67dlrVL9/P7w3wuV7J2Xy129JyTcxJHy7Wt577v532f93kZQghSSCGFlQ9Roi8ghRRSEAYpMqeQwipBiswppLBKkCJzCimsEqTInEIKqwQpMqeQwiqBJMT/p/pWKaQQezBCvEgqMqeQwipBiswppLBKkCJzCimsEqTInEIKqwQpMqeQwipBiswppLBKkCJzCimsEqTInEIKqwQpMqeQwipBiswppLBKkCJzCimsEqTInEIKqwQpMqeQwipBiswppLBKkCJzCimsEoSaZ04hBmBZFna7HSKRCBKJBGKxGAwjyEhrCqcwUmSOIwgh8Hg8cLlccLlcYFmWI7FEIkFaWhokEglEIlGK3ClEDCaECX7KaUQgEELgdDo5ArtcLq//4/8xmUzIzc1FRkZGitynBgT55aYicxzAsiycTicIIX5JyTCM19enp6eRkZHh9f8SiYT7kyJ3Cv6QInMMQQiB2+2G2+0GwzAQiUTc1wMRG/iS3GKx2Ot1aDRPkTsFf0iROUbwTasjJRv/+OP7fEIId+6m/0/P3GKxOEXuUxQpMscAvlE0UmKFejw/agP+yU2LaSlynzpIkVlA+KbVoQgU7AwdyXZOf+R2Op1wOBwAAJFIxJFbIpFEdYNJIfmRIrNAYFkWU1NTkEqlUCgUQclCyRSrdbrByE3fOy0tjUvLU+ReHUgpwJYJGo0dDgfm5uZgs9mWTQyhiU7JzS+WOZ1OWCwWdHR0wGw2w2q1cr3vFFYmUpF5GaBnVY/Hw1Wrk315Pb3R0Mit1+tRVVUFp9MJp9MJAJwyjV9QSyH5kSJzlPDtHQuZOscyBff3Xr5tMJqW88ntW1BLIfmQInOECNQ7BuJLwljBXxuMEAKHw+G3oJYid/IgReYIEKp3LFSanUw3hXDILRaLuZScVstTiD9SZA4TNBrz02pfMAyz6gtI/shNp8Doz4aSOzURFl+kyBwCwdJqXwhF5mSKzKEQjNwUKXLHBykyBwHLsl6jiuEos1YKCWOFYOSenp6GQqGAQqFIkTsGSJHZD2g0Hh8fR3FxcdgFnpVYzY41+OS2WCzIysoCy7Je/fhU5BYGKTL7gPaO3W43RkdHUVpaGvZzV0KfOZEghEAkEnF/6Nd8yc2fCEuRO3ykyMwDv3ccTbslkohqt9u5iLSc11lJ8KdFp5GbT27qxuLrwiIWi1PkDoIUmRFZkSsYwimAEUIwNDSE2dlZAIvRXKlUQqVSQaFQcHLL1Uhmvk1SIAQit9vt5h6Tsljyj1OezMudO+YjFAmdTic6OzuRlZWFlpYWAIDL5YLJZMLMzAz6+/shlUrBsiykUilUKtWq+qBGk/H4K6jxyZ0yavgSpzSZaToXrHccCYKR2Wg04sSJE1i/fj0KCgq4KrlUKkVBQQEKCgoALKbffX190Gq1mJqaQmZmJlQqFecJtpI/qMHcVcKFP3KnXFgWcUqSWai02hf+yEwIwcjICDQaDbZs2eLl7eUP6enpyM7OhkKhQF5eHqxWKwwGA06ePAm73Y6srCyO3DKZTJDrjheEILMvgrmwaLVaZGdnIzs7+5RwYTnlyBxp7zgS+FazaVotl8uxbdu2iFtcDMNALpdDLpejvLwcLMtiYWEBBoMBJ06cgNvtRk5ODnJzc6FUKv0W05IJsSCzL/hDI3Nzc8jIyDhlXFhOGTLzq6RAdHY+ocAvgJlMJnR3d6O6uhqFhYWCvL5IJOJEF5WVlfB4PJibm4PRaMTY2BgIIcjNzUVubi5ycnK8DAr8Id5FtniQmQ+WZbkKOP8a/Bk1rAYXllOCzEIWuYKBknlkZAQzMzPYvHkzMjMzo3qdcIgmFouhUqmgUqkALOrHjUYjdDodBgcHIZFIOHJnZ2cnfLopnGq20O/n+z2Ha7G0El1YVj2Zae/4o48+wle+8pWIfzGRRBOPxwO9Xg+JRILt27fHnTwSiQT5+fnIz88HADgcDhiNRkxNTcFsNiM9PZ0jt1wuj+u1AYmJzKF+B76z3ABWrFHDqiWzr/ggmjusSCTiUrVQmJubQ3d3NzIyMlBXVxftZQMQTjQik8lQVFSEoqIiEEJgs9lgNBoxMjICi8UCuVwOp9MJm80WsjAnBKIV40SLcMjMh68Ly0oj96oks1BpdTikIoRgbGwM09PTqK2txdTUVNivHU8wDIPMzExkZmaitLQUhBCYzWZ0dnaiv78fDocDCoWCi9xSqVTwa4h3ZPZ4PMsi20oj96ojsz87n1DomTYjSyZGucr7fEsjcyC4XC50d3dDKpVi27ZtcDgcYY9ABruueCjAGIZBVlYWZDIZmpqawLIs5ufnYTQaMTk5CY/HA6VSyVXKJRJhPirJlmZHAn/kDmSxpNFoUFhYiPT0dMHePxRWDZnD6R0Higw//78BfDRkQE2BHOfUFuCyzcUoV2UGJfP8/Dy6urqwdu1aFBcXA4iMhMkg1+RfA5WVKpVKrF27Fh6PByaTiUvLGYbxqpQnU3oZCLHOBIK5sNx11114+OGHUVtbG7P398WqIHM4aXWw8++/XVKLt3u1eLtPi6feH8bTH4zgsi0l+Jrqy82MwJdkHR8fx+TkJJqamrwKSStxBDLQh10sFiMvLw95eXkAFrMQo9EIjUaDgYEBSKVS5ObmQqVSISsrK2krvvG8Lv5nz2q1IisrK27vDawCMoe7CiYYmctVmbj+KxW4/isVmJ2346n3R/BS6yT+10MgevNDsF/wSiJiIBEBF6/PxL9dsX3JawlJ5mRDWlraEtmp0WjE+Pg4zGYzMjMzOXKvdNmpEKCz2/HEiiVzpKtgxGJxWOfZQkU6/uWiWnz3qxV46u9tyMpRIkMmhcPhxOTMLD7VAON2qd+bQqTzzMHSwGRIw4MhPT0dxcXFKC4uBiEEVqsVRqMRg4ODsNlsyMrK4si90mSnQsBqtca9/bciyRyNJFMkEsHj8YT9HqXKDFyxMQtr1pRgfn4eExMabDp9K+766wAMFpff54QzAhkO4pVmC/UefNlpWVkZVymnwyVutxsKhYKTVSa77FQIsCwrWNEwXKwoMtNe6cLCQsSKpmiJNjAwAJlMhu3bF9NquVSCcaMt4Hske0T1RSzSYYZhONlpRUUFWJbF3NwctFotOjo6QAjhZrjDkZ2uNPDrLPHEiiEznYaZn5/H1NQU6uvrI3q+WCyO6AdsNpuh1WpRVlaGmpoa7utymRgWh/8IL+SZeTVZ9opEIm7Kq6WlBW63GyaTCXq9HoODgxCLxVxKLpTsNBluqvGuG6wIMvN7x2KxOKJ0mSKSNHtiYgJjY2NQq9VcNZdCLhXD4nT7fd5Ki8yJulaJRAK1Wg21Wg1gUYTBl53KZDJuzFMul0dFingLVJIBSU1mf73jcAtZvgglAAEWFUMnTpwAIQTbt2/H0NDQkufIZRJYnZ6AflZCIBlaU/GEVCpFYWEhN11ms9lgMBi8ZKd8g4ZwILRgJBI4nc6YKOhCIWnJHKh3vJzIHIzMCwsL6OzsRHl5OUpLS7n39CVVplQMQgCbi0WmNPqzXjKQKFmRkZGB0tJSTnZqsVhgMBg42Wl2djZH7kCkSSSZLRZLVNNyy0VSkjnYKphwIqw/BHve1NQURkZGsGnTJmRnZwd9jly6+COzONzLInMwrLR0PRxE+/1Q2WlWVhbWrFkDlmVhNpthMBiCyk4TSeaFhYW495iBJCNzOJJMISOzx+NBT08PPB4Ptm/fvqSV4K9vLJctEtji9CA/4qtYhNvtxokTJ2Cz2ZCXlweVShX12XA5iOe5Uqj3EolEyMnJQU5OTlDZaSKFK4lQfwFJROZwe8dCRWaLxYKOjg6UlZWhrKzM7/v5qypzZHb4L4KFAv9916xZA5PJxJ0NafpIs5LVhFjdOILJTo1GIz7//HMuJc/Ozo4LwU/ZNNvXziecYfJowCfz9PQ0hoeH0dDQAIVCEdZzKLg02xl5djA7O4vBwUE0NDRALpfD5XKhpKQEJSUlnNDCYDBgdnYWbrcbTqcTKpUKSqVyxfdi45UFUNlpeno6JBIJKisrYTQaMTEx4SU7zc3NRWZmZkyu6ZRMs2nv2OPxxNyeRSQSwe12o7u7Gy6Xy29a7Qt/Z1cama0RkJllWQwMDGBhYQHbtm1DWlrakqMCX2iRmZmJhYUF5OTkwGAwYGhoCBKJhLMIEmqwYSWm2eGCnpkDyU6HhoZgtVqRnZ3NkVuoccVTjswsy2J2dhZKpTIuPksulwtjY2OoqqpCeXl52BJQmjFQ8Atg4cDpdKK9vR25ubnYsmVL2N+nSCTySh8dDgcMBgPGxsa4Dwsl90rQPieDy0gw2Wlvby+cTifndpqbmxu17PSUOTPzi1w9PT346le/GvP3nJmZwfj4OAoKCrBmzZqwn+e3ACalZ+bQkZk6dNbU1HC+XOHAX0Ygk8m8IszCwgKMRiN6enrgcrmQk5PDnQ2TMSVPBjM/XwSSndJpMCo7pZXycH+utDceb8SVzPFyyaRgWRa9vb1wOBxYv349rFZrRM/3XwCjZ+bAkZkQgomJCUxOTkbl0BmqNcUwDGfuvmbNGs5ylwotRCIRF7WDFX1OhTQ7ElDZaW5uLgBwslN61AlXdrqwsCCYvXIkiBuZhVwFQ0kW7JdltVrR0dGBoqIi1NXVQafTRdzS8lcAywwRmQkh6OrqAgBs27YtLlHS13LX6XTCYDBwRZ9oFFRCI1Ge2ctBINnp9PQ0+vr6IJPJOHLzW4tCzDIzDJMO4D0AMizy9AAh5F+DXu+y3jEMhOodR/NLpiQLRObZ2VmcPHkS9fX1UCqVAMKfZ+bDX4QUixhkpIn8Rmar1Qqr1Yry8vKQ5/JYeoBJpVIvV05fBRVNHeNJ7JUQmUPBn+zUaDRidHQUCwsLkMvl0Ov10Ov1QpyZHQB2EkIWGIZJA/ABwzB/I4R8EugJMSVzqN4x/dBG+kumwhHfajTLsujr64PNZsP27du9ChjRLEIP1NOWyyRLWlNarRb9/f1IT08Pu8AWD/hTUNGUfHh4GHa7HcPDw4JOLPnDaiCzLzIyMpCRkcG1Fi0WC9ra2vD+++/jzTffxMGDB3H//fejsrIy4tcmix/WhS/+mfbFn6Af4JiQOdxVMJSU0ZxtfElms9nQ0dGBgoIC1NbWLnm/SM0JAr0PsJhq0zSbEILBwUEYjUZs27YNn3322bI/uLGUc/LPhSUlJejv70dmZiampqYwPz/PbZ2k9j9CIRFkjqc5AL1pfutb38KHH36I22+/HYQQ5OTkLOc1xQBaAVQD+DUh5Giwxwv+3UZi50PJHGkLwFfSSU3mNm7cyBUvfBGNcixQNKdjkC6XCx0dHdy+ZbqELJG64EghEom41JH2YfkpuRCtGiD+WmmPx5Owqr7VaoVSqcSGDRuW9TqEEA+AZoZhlAD+wjBMAyGkK9DjBSVzpHY+YrGYW5odCej5l2VZ9Pf3w2KxYNu2bUHHzqIhcyCTALlMgnmrA59++inWrVvnVbkUIqomyjaI34elWyf5rRoAXMFHoVBERM7VmGYHgtCiEUKIiWGYdwDsAhB7MtO2ExBakkmxnKEJq9WKnp4e5OfnY8OGDSE/KEJGZjHrwuycDU1Npy3pJ660iadgPzffVg3VPc/MzHD1AX5KHuy1TiUyU539csAwTD4A1xdEzgBwHoCfBXuOYGSOpt0UrdGAw+FAX18fGhsbA6bVvhAiMtO+tYh1ARKZX2FANIU2f++bjDcEvt0u9WPjL4KngyIqlWpJSn4qkdlmswkxaFEM4PdfnJtFAF4ihLwW7AmCptmREibSNJtqnM1mM6qrq8MmcjTX5vscu92O9vZ2FBYWoqxQjcEBvd/nrDb/rkBgmC93V5WVlXnNGU9MTIAQwqXkOTk5pxSZqb3VMl+jA8DmSJ6T0EELiUQSdpptt9vR0dGBvLw8lJaWRvzDijYyE0Kg1+vR29uLuro6qFQqZPb0xdQHbKVZ7QJL54zprmhanGSYxaXmdDwwHuq/RJA5Uc6cQILJHO6ZWafToa+vjyPTyMhIxGftaAjCMAwsFgtOnjyJlpYWbqomlA9YMqbIgRArUvnuip6YmIBer/eaVqIpeSz8shLdUUiEzkBQMkej5ApGSkIITp48CZPJhK1bt3LTQdEUziK9Njou6fF4sG3bNq8PhvwLHzCr08NptfnvEy6ZA13TSrshhIO0tDQolUpUVFRw00p6vR5dXV1gWVZwH+1ERuZECYYSnmY7HA6//+dwONDR0YHc3Fxs3bp1iQ+Y72iikFhYWEBHRwcqKirgdDqXfCi+HLZYSmahCmDxQDxvGPwPOX9aiabkJpMJOp0OJ0+ehFQq5aJ2tHZKiSJzopw5gSRNs+kZtba2dolvNX2e3W6PyTXNzMxgaGgImzZtQlZWFkZHR5c8xss6KNt7llioAthqs9oNFrF8BxrsdruX1W40s9uJIjPVaCcCCU2zfcnMl0byz6i+iNYHLBioAMVqtXJuIIFAJ6f8uY2spAJYPBFJ+pmenu5lp7SwsAC9Xs/traIpebAZ40SSORHGBEASRWaaViuVyiVpdbDnCQH63iqVCps3bw75oQvmAxYOEelAiFar5T6YKpUq7ovG4oloz5L82e3KykrOjTMcO6VEnF0T5TICJAmZDQYDenp6sGHDBi7VCgYhIzPdVBjuewPBHTpDkdnpdHI3rW3btnG92bGxMc5UIDMzMy696nhGfyHmi4Glbpz+7JRyc3MT1utPlDMnkATV7Lm5OdhstqBptS+ijcz8IQhCCMbHxzE1NYUtW7ZENCEULDIHK4CZzWZ0dnaiuroa+fn5cDqdUCqVUCqVqKqq4kwFZmZmOMsh6qsdq6JKPM/MsUh7A9kp2e12fPrpp5ydEt8gP5Y4JdNsp9OJEydOwOPxYOvWrRH9oqMZZ6TPY1kWhBB0d3dDJBJF5QaSFSIy+4sK1Ga3sbERWVlZfglPTQUUCgUGBwexZs0ar/YNTScjHXJIBsSjZcNPyWdnZ7Fly5ao7JSWg1OOzDS1raqqwvj4eMQfzOUsj1tYWEBPTw/Ky8tRVlYW8WsAfIfO0Gdm2iufn58PWVjzfR3+WZEqquiQQ0ZGBhe1o7WITVRrKl4Ix06JSk6Fmt0+Zc7MhBAMDw9Dq9Viy5YtkMlkGBkZifh1ok2zXS4Xurq6sGnTpmUNjWeEWc12u93o6OiAXC6PyGbX37mbr6iic8e0hedyuaBUKpGXlxex6CIZWlPxgj87JaPR6DW7Tb3Sok3JT4nWlNPpRGdnJ7KysrwUVdFEh0gLYDQ62u12tLS0LIvIQHAfMEpEi8WC9vZ2rF27FsXFxct6P3/vQeeOqTsnFV0MDg5yoou8vLyE7lziI55kDuezwbdT4s9uGwwGjI6OgmEYr5Q83OzRYrGgtLR0ud9CVIhLZKbFnPXr16OgoGDZrxdJZKbVY3rXjaai6u+DKJdJ/KbZtKg3ODiITZs2BV1/EwiR9pl9K7y+o4k5OTnIy8uLWxHIH+Lpmx1Nsc13dpu/AJ5vp0TX2gSCzWZbHZHZF4QQjIyMQKPRRFwxDoZwPxRzc3Po6uribiKdnZ1RT04tIfMX1kF8EEJgMBg4Q8FEyfr4+435EWdkZIQ7R+bl5SU9waKFEIIRvhMn306J3iAVCgVHbn4dRIgCGMMw5QD+AKAQiyZ+vyOE/CrU82KWZrtcLnR2diIjI2PJoEI8MDExgfHxcTQ3N3N3yuXMNPvTZ/Mjs8fjQVdXF9xuNyoqKpZFZCEVYL4Rh/ZlR0ZGYDabwTAMZmdn/RoKCIl4p9lCft782SnNz8/DYDB42Sm53W6hqtluAP8fIeQzhmGyAbQyDPMmIeREsCcJHpkZhoHJZEJXVxeqq6vDcvYX8hfNsizX8tq2bZtXWhkLUz9gMbVqb2/nhvSXi1h+6Pl9WZPJhPHxcVitVkxMTAAAF7WFbt2sZDL7QiQScfoA4Es7pd/97nf429/+hoGBAVxxxRX49re/HVV9hhAyDWD6i7+bGYbpAVAKIH5kpmn19PR02GtZKMGEUAdRUhUXF2PNmjV+7XaFMvXLlIqhW3By6jXqDEp3FC0X8WobyWQyrF27FmvXroXL5YLBYMDk5CTm5+chl8u59tdyl9OtJjL7gtop3Xfffejr68Ptt9+OgYEBoW7slVh0HAlqswvEIDKLxWJs3749YlO/5ZKZtmniZrcrk2BQM4+BgQEv9ZoQU1OJqj6npaV5nRMtFgs34ODxeLxsgCIly2omMx8LCwvYuHEjzjnnnGW/FsMwWQAOAriTEDIf6vGCn5np1oRwsZyhCUqy4eFh6HQ6LwODQO8lRGRmWRZOyzwsjkX1Gv9GRPdALxeJnprit24qKiq4mWNqA5Sens5F7XAKm/EstgmV6UUDIZw5AYBZXElzEMCfCCGHwnlOwsd0lmO363A40NPTg4yMjLAkocs19QMWC0htbW1QZErhYJ1LPjQraQQykmjJnzmmzpx6vd5rfxVtf/kjUjyr2dFsSREKdrt92V0bZvGX8gyAHkLIf4b7vBVLZkIIWltbUVVVFbYoIxpNN59YtNVVW1uLYosJ1u45sCyBSOQ9drfa3Tn5zpzl5eXcSlnq8ZWWlsYV0qh536mSZgt00/oqgG8D6GQYpu2Lr/2EEPJ6sCetSDJPT0/DYrFg8+bNfp1IAmE5+6YmJycxNjbGFfbkssWdXlaXB1k866CVFJmFgq8GmjqFDA0NwWazQaFQwG63CzqDHgyJ9P8S6HU+ABDxnS8mralIEAmZ6VC/3W6PasAgGu8whmEwPDwMQohXq4u/p1loMscLsbpOvlMI7ckajUb09PR4iVZ8zQSEwqnozAmsoMjscDjQ3t6O/Px81NbWRqXmivTM7HQ6odVqoVarUV9f7/VLogYFVqcbi/uwF7HSInOsP3i0J5uRkYGNGzcCgJeZQCwsd09FZ05ghZCZjkzyDf6iSc8jITM1EsjJyUFRUZEfOad/gwIh3DlXI2g1Oy0tzWtyiTqt8Ge2qWglWkImisx2uz3qcVQhkNRpNiEEo6Oj3KA5v0oo5L5lX1CHzsbGRkxPTwfoM/s3KBCqz7zSNlqE816+BONb7tKZbYPBgOnpafT19UU9sx3v3cwUiRx/BJIkMvvzzqYm9BKJxK+2O5qecSgy+zMSCETOQJE5XCK63W6Mj48jJycHCoUiYelZMs0zSyQSr8V0vjPbVLSiVCqDRt5T0ZkTSBIy+0ZYi8WCjo4OrFmzJuBsqNBpdiAjgUDP+TIyR05mq9WKtrY2qNVqTE1Nobe3F1lZWdwYYywHHhKFSM+TwWa2T548CZlMxkVtX9lwoshstVpXV2RebppNvbIaGhqCzgILIQChCGYkEGzQAvCfZgcjM9Vy19XVIT09nXssHbHr6OgAsFjwM5vNMav4JgLL+T78zWzr9Xq/M9uJIrPFYlldZI4UlMyEEG5dazheWUJFZq1Wi/7+/oBGAqHSbF/roGAFsImJCUxMTKC5uRlSqRSEEM5gkIow6Axya2srV/FVKBRQq9XLsrPxh5VcqMvIyEBZWRk3qUZFKyMjI3A4HPB4PJBKpVGvt4kGp3yaLZFI4HQ60draCqVSGbZXVjQaaD6Z6YSXTqfDtm3bArZFAkVzrs/sDB2ZCSHo7++HzWbjvj+GYSASiTjpIyW1x+OBxWKBVCrF+vXrwTAMFhYWODsbfoQSYjXqaoj6vjPb3d3dSEtL49bbKBQK5OXlLTESEBqJNPMDkiDNtlqt0Ol0aGpq4tZ/hgOxWAyn0xnRe1FiUiOBtLQ0tLS0BE3JAkVmkYhBplTs98zMfzw9i2dnZ6O+vp57jL/xTGBx+mt4eBibNm2CRCIBy7LcwMOaNWu42VmqrqLpZW5ubsKGC5INDMOgoKCAszT2NRKI1cz2KV3NHh8fx/j4OBQKRUREBqJrTYnFYrhcLnz66adcihbO+wRSjS0aFAQugNlsNrS1taGiooKr0PojMvBlG85gMGDLli1eEYRlWe5PWloa8vPzOS81+kEdHh6GVCrlonY4Yv+VnGYHA//MzDCM1xJ4OrNN7XajWUoXCHTJXaKQEDJ7PB709PSAZVm0tLSgra0t4teIpjU1NzcHs9mMrVu3Bpx59kWwM/BiZHb7fTwVumzcuJGLEIEyAJZl0dvbCwBobm5e8jiRSMR9jabjlNwKhYLz1qYmdP39/dy2DFoUCvTeqyHN9kWwAlgsZ7YtFktEswJCI+5ktlqt6OjoQElJCcrLywGEZ43qi0giM38VTUZGRthEBoKLQOQyid/IbLPZ0NvbyxW6AkVjYNFypqOjA2q12q87ii/oB4x/1qZHB18Turm5Oa6VQwUYeXl5y45AyY5wq9n+ZraNRmPUM9unVGuKVo7r6+s5/yQgunQv3MhMPcHooMTRoyHdV7wQrAUml4q/0GYvghCCoaEhOBwOnH766RyJA/1M6I2tqqoqagtiGrXp+ZpGbY/H4+VTZbfbuWyBRiCxWLwqU+1oW1O+iwb4M9s00wm2SnZVVrP9rWgZGhqCwWAI6QYSLsJpTVEjgcLCQlRUVESVUgZLs+UyCbTmRfWa2+1GZ2cn0tPTIZVK4XQ6kZ6eHvA9jUYjent70dDQIIgzBb1WYPFnk5aWxpGaEIL09HQUFxdzmmiTyYSpqSksLCzAZrPFfEFdPCFEnzmcmW2+aIVhGEHOzAzDPAvgYgAaQkhDJM+NeZpNLXflcnnIynEkCCUa4RsJLOccEzTNloox4vRwha41a9agoKAAMpkM3d3dYFkWeXl5yM/P96qcTk1NYWJiAps3b46pMN/fWdvj8XADDR6Phys+GgwGdHZ2AohdtTdeiIVFkb+ZbUpsm82G48ePQ6PRCPFWzwF4Aou+2REhpmSmk0dVVVUoKioS9LWDRWZfI4HlIFRkNttc+Oyzz7hCFwBUVFSgoqICLpcLer0eo6OjnPjD7XZzhb94tpJ8z9rz8/MYHx9HTU2Nl2m+x+OB0WjE+Pg4d800Aq2UZfDBjjZCIT093WvRgN1ux8GDB3H77bcjJycHBw8e5IgfCQgh733hyBkxYpZmT05OYmRkhFthGgqRanf9RWZqXuBwOJZ4ZkeLYJGZddqw4HChqakFMplsyYeIP+7ndrvR1tbGRcj29nao1Wrk5+cLtukjXBgMBvT396OpqQlyudwragOAWq1GXl4etzWTvwyeFtEiUVatxnM5HyKRCGeccQays7PxwgsvAEBERVahEBMy9/b2wmq1Yvv27WERikbZSMjnG5mdTifa29uhUqlQW1sr2J3Z302DTlexDiscHiAtRMXa4XBwFXw6OGKz2aDT6dDb2wuHw4G8vDyo1eqoWiKRYGZmhstaaO3CX4WcpuP8TQ50RHFkZARWqzVswUqih/bjBWq2sNxsMFrEhMxlZWURbR+Mhsz8IhtN56urq8OqCkfy4fIt5nk8Hq7QVbWmBOgfgt3FIivd/4fZbDaju7sbNTU1XmlXRkYGysvLueIKneOlE1Q0OgpZkBodHYVer8eWLVuC/qx9z9r0D6320p8xNRYYHh7mCkJUZsrHqUJmp9OZ0LZfTMicnZ0dkTormqEJ+uHgGwmEk85H6hTJj8x2ux1tbW0oKytDUVERWk2TAACri0WWnzqWVqvltkEG6z+KxWKvlsjCwgK0Wi3a29sBgEvHox0aoNpwl8vlV5QSDMEEK/xl8FSwQqeYcnNzOcHKqULmeNoJ+0NSVDSiITMhBHa7HZOTk2FNWVEEWgQX6vG0Ol5XV4fs7GwQQpCVvvieFocH4HWXqEiFLpWPJLoyDMORpKqqCk6nEzqdDkNDQ7BYLFAqlcjPzw9bi+3xeNDd3Y3MzEzU1NQsi1TB0nG6ooVGbb5gRSaTwel0xsVWJ1Hnc6Hel2GYFwGcBUDNMMwEgH8lhDwTznNXJJnp8AKAsKesKCKdgxaJRLDb7Thx4gSampq8Cl1yP5NTtAjHsiw2b94syGpRvtOlyWSCVqvlSJKfnw+1Wu2XJC6XC+3t7SgqKgpLhx4pgrW+qB6aYRiYzWb09fVxjiG09aVQKASPZCvdmZMQsi/a5644MvONBOx2e1RmCOGSmY5J2u12nH766RCJRF6FLr7dLvBlTz03NxeVlZWCp5Yikcir12mxWKDT6dDd3Q232831tKlPNVWXRTrEEu21AUujNiEEbrcbGRkZqKurAwCYTCbMzMygv78fmZmZ3FlbiPpAPNfg+L5vIm8iQAxbU5EgXDL7GglQP+tYRGb+mKRUKoXb7eaiMgVnHfSFcKSjowOVlZVhrbEVArTSTHXFer0eExMTMJlMcDqdqKysTEiLBPgyatPqd0NDA3cjzc3NhVKp5HTsRqPRy51TrVZHLVhJ1J4pm80W9xajL1ZEZCaEYHh4GHq93stIIJq1J+GQ2W63o729HSUlJSgqKoJUKkVHRwdEIpFXMYq6jehM82gzDGDjxo1R7eMVAhKJBIWFhZBIJJifn8fGjRuxsLCAzz77DGKxmNsTFc9BAI1Gg5GRkSV1A/5ZOyMjAxkZGdwxwmg0cuOJ2dnZnGAl3JpIIs38EjlkAawAMtNWkFQqXSIHpc+LtDob7MYxPz+Pzs5O1NbWIicnB4QQbn+xw+GAVqvFwMAA7HY7kLFI3KHxKey+eGtCPZOBRZno5OQkWlpauAmqdevWwW63Q6fTcdedm5vL2RDF6oM/OTnJ7en2JWKgszYALuWmDiuU3AzDcFE7WFU/kWZ+iRyyAJI8zaaa5/Lycr8FHPq8SKxggkVmaibY2NjIEZP/wZDJZJypgdvtRmt79+J1ulgMDAwgPz8/Ie6a9GxvMpmwZcuWJWlmeno6d91UrkmPLJmZmVwRTaieNjVZ2Lx5c8iUl3/WJoR4FdHoMYJeN03ZLRaLl2CF3zNPReYEw58FEHWx9B2X5EMoh0461WU0GjkyBFN0UXMFhXxRIKAuKkNFhRparRZjY2Nc31itVsdcDUQIQW9vLwghaGpqCvlB5qfcdDhfp9Oho6MDhBBOiRbNmZUQgsHBQdhstrCuxRd8bzTAW7ACgOvF0wq5P1+0RK1zTbQzJ5BEZKaRmRCCsbExzMzMoKWlJWjqKoRDJ+3DSiQSNDU1AQgu1Hc6nejo6EBhYSHKy8uRkaaB1enhNjPQtFar1XI6cVplpq0aoUCPIAqFAmvXro34tfnD+ZWVlXC5XNDpdNxgSE5ODqdECxVhCSHo6emBSCRCQ0ODIN9nIMGKx+PxMhXg+6ItLCxALBZDr9cHnDuOBRI9ywwkWZrNNxLYunVrWCnaciIznXcuKipCaWlpUI8uYPEX1tXVhfXr13NjlXLZUh+w9PR0L6mmXq/H5OQkenp6kJ2dzaXjyxkEoVr00tJSlJSURP06fKSlpaG4uBjFxcWcUEan03FyTZpt+FZtWZZFV1cX5HI5qqqqYtIaCjarzTcVMBgM0Gq13ErZSH3RosWqPTNHCppmf/rppygqKgrLPoc+L9rIbDab0dHRgQ0bNnhJDgO9r16vx8DAABoaGrx+aYtuI4GvQSwWe61cmZ+fh1arxcjICNLS0qKanKIOJdXV1VCr1eF/8xHA176WDob09PRwwg+1Wo2srCx0dnZytkfxgr+oTW1/srOzUVJSAoZh4HA4vHzRqMxU6IGW1Jn5C1itVmg0mrgtTzeZTJydLSVRsF/s+Pg4ZmZm/Eoz/UXmQOA7RVZXVy8hCF/0EeimMjc3hxMnTqC+vj7oxg+h4W8wZHJyErOzs8jOzkZaWhpcLldCVuvQmfORkRGIRCJUVFRwBKdVfb7MlHp8CemLlmhnTiAJ0mw690ydJCNBpA6dhBAYDAZYLBZs27YtZKGLDig4nU6/VWJgcbOFr0NnuOAThIo+xsfHYTabOQcQ/nmVDm40NzcnVKAgFos50Q69IWq1WrS1tYFhGO6mFK9tErQIKBKJuPFXf8sFWJaFUqnkahe+vmh8mWk0/u/xEgsFQsIiM99IYPPmzThx4kTErxFJZGZZlpM9lpaWhiQy9fRSKBRBBxTkUjFmzUu3WEYKKvrgO2tqtVoMDQ1BJpNBIpFwGzES7dNF0/wNGzZwabjvYMjg4CCsVqtXTzsWxShCCE6cOAGpVIrq6uqAywX8DYdQXzRaI6C+aNEs8lu1kRkIvkDN10jA4/FEnC4D4Udmp9OJtrY2FBQUoLi4GD09PVhYWFgS+Siorrm8vHzJIjlfZPoxwl8uGIbhnDWrq6vR19cHo9GItLQ0tLW1+fUVixfMZjO6uroCpvnBBkPS09O5tpgQAhtaMM3IyAi78BZoVhsAp3un5nxGo5Eb6KHEDrTI75RsTdHC0/r167lzTDSFLCC8yEzfr6amBrm5uSCEYMeOHVwhamhoCOnp6Vw11Gaz4cSJE6irqwvY3+ZDLlu6okYosCyLnp4eiMVi7NixAwzDLPEVy8nJQX5+PlQqVczbMCaTCb29vWhsbAzrgxtsMISmtaFqBIFAK+jZ2dlYu3ZtVN+PL7EBcFHbd5Gf0Wj0WuTn64smRGuKYZhdAH4FQAzgvwkhj0by/LiSmRoJNDU1eX3j0UYXsVgcdHkclV7Scx1/eJxGvvXr18NisUCr1eL48eOw2+0oLy8Pu5CTKZV4eWcLBZrm5+bmetkE833FaPuInqVp5MvPzxfc8YKuT21ubo46qvIHQ/hrYubn57lNl+G07FiWRWdnJ3JyclBZWRnVtfiCfi58ozYVofBlpmazmSO3WCzG6OgozGbzsiIzwzBiAL8GcB6ACQCfMgzzCiEk7PNnXNLsSNe1hotAkZnubaIVcolEEvR8nJmZyXlLNzU1wWQycTrmUIIPuVQMm4uFhyUQi4RJeR0OB9rb20Om+b7tIxr5Ojs7wbIsR+zl7nienZ3F6OgoNm/eLNh53XdNDM2UqKKLXruvgs7j8aCjowN5eXkxbYUFWy7gu8jvxIkT6OrqwlVXXYWzzz4bP//5z6OJ0tsBnCSEDAEAwzB/BrAHQOLJTEHXr2RnZ0dsJBAK/s7MfOHJ5s2bAQRXdNFUViQScZY6crmcs53lCz78VZjpGOR/vDWIc2vVaC7LWRapLRYLOjs7l3iGhQPfyEcFH3yHEpVKFVF/dWJiArOzsyF9w5YDfssOADcY0t/fzw2G0BpBZ2cnCgoKYmK2EAihBCtXXnklnnzySbzyyisYGRmJNkKXAhjn/XsCwGmRvEBMyUyNBGLhmw0sjcy0sKZWq1FeXh5SCOJ0OtHZ2Yn8/HyUl5cveZyv4INfYabn7NPX5uCjqlz86dgEfv/JOHIz03Dm+jzs3KDGV6pUnIFBODCZTOjp6RFky4WvmotfiMrIyAg5XEH7tvPz82hubo7rjLC/wZCZmRm0tbVBLpdDJBLB6XQmrKrve9Z+/vnnQQiBUqnEmWeemZBrAmJIZq1Wi97eXs5IIBxEOpvMj8wLCwucKkqlUoUkMo2A69atC8uJg19h5p+zTVP9uGkDcOf2Neg3S/DR6ALe7tXhr+0zkElE+Ke1udi5QY0za/KQnxX4HKvRaDA8PByTLRf8QpS/4Qpfw0B6LHK5XNi0aVNCHTTEYjGUSiVGRkZQW1sLhULhZXZIj0HLPUpEi5deegn/+7//i48//ni51exJAOW8f5d98bWwwYQwIovapWxychIKhSLsu+fRo0fR0tISUSpnNpsxPDyMkpIS9PX1ha3oMhgM6OvrE2zPk9PphFarhVarhd1uR06uCpPODBybtOEffXpMzdnBAGgsVeDsDWrs3KDGOnUm9+EbHx+HRqNBY2Nj3BVUtC+s1WphtVqhVCphs9mQmZmJDRs2JNxV0+Vycat/fEUZTqcTer0eOp0u7pV9ADh06BCefvppvPbaa8v9HDEMw0gA9AM4B4sk/hTANwkh3WG/SKzI7Ha7I2o3tba2oqGhIaIqrMViQVtbGyQSCRobG0MWuoDFm8zU1BQaGxtj4nFMpY4ajQbz8/PIzs6GWazA5xo3jgwY0DVlBgCsUWXg7Jo8bJA7UCH3oKkxsREQWPydff7555xiSi6Xc33hRKS0VB9QWVkZ0g+dPxhiMBiCDoYIgVdffRWPP/44XnvttbBamCHAAADDMBcCeAyLralnCSEPRfQisSKzx+MJ2jbyRVtbG2pqasKe/6V9Rr1ej9NPPx1A8EIX3UJhtVo5P6pYg3/O1uv1i1shM5XoMjJ496QRnwwb4WaBnAwJd87+6joVZ0cUT7jdbrS3t6OgoICrNywsLECn00Gn0wH4cp6Ybj2MJZxOJz7//HOsW7cuqmESm80GrVYLnU7HDYYINYb6xhtv4Be/+AUOHz4c1T4pPxDkh5k0ZO7o6MDatWvDSldo6qVUKrlKK11qHuha6IjeunXrEpY60nO2RqOBxWKBLCsHOkk+Pho140i/HvN2N9LEDHaszcXZNYvpeEF27Dck0MJheXl5wEKlw+HgiG2z2bgKs1KpFDyjoKOp1dXVy9rgSUHX6uh0OszNzSErKytqV5i3334bP/3pT/H6668LObG2ush84sQJlJSUhExZaIV83bp1UKlUGB0dhVar9doKwU+rqDSzrKxMsLnf5YBeD905xZ2zlbmYdmfi2KQd7/TrMG60AwA2lWQvnrNr1FhfIPzgAjUvjCQC8q2HTCYTt05HrVYv+8xPt4bwdd9CghACs9nM3ZioSWMobzEAeO+993Dffffh8OHDQg9VrC4y9/X1cQqgQNDr9VyF3LfQRd09NBoNPB4P98sZHBxEXV1dwixn+VhYWODMAvnX4++cbZEo0Kbx4J0BAzom5wEAZcp0jthb1uQgTby8iEgr+rW1tVGf+2g6rtFooNfrIRKJuLNqpNVdm82G9vb2ZV1PpHA4HNDr9VwBMNBgyIcffogf/ehHeO2112IRFJKbzCzLwuVyhf34kydPQqFQBCx0jI2NYXp6Gps2bUJaWlrQ87HL5cLQ0BCmpqYgk8mgVqtRUFAguG1PJKBrVH3NDXzh75wtkuei28jgvUETPh4ywulhoUiX4Iz1eTi7Ro2vVauQJYvsnE0HJoSq6FNQwYdWq4XD4fA6qwZLx61WK9rb2xNqV0w12DqdDkajEenp6RgYGEB6ejoeffRRvPLKKygvLw/9QpFjdZF5eHgYMplsyV2Pjko6nU7U1dVxJA5W6BobG4NOp0NjYyNnxE6jXk5ODgoKCiJWQi0HdI0qXW8TCeg5W6fTgRCCbGUeTlrS8NHoAo7062GyuSARfXnOPqsmD8U5wfvURqMRfX19aGxsjKnhIM04tFot5ubmkJ2dzWVf/HTcYrGgo6ND8BvLcmGxWPDQQw/hpZdeQlFREfbu3Yv7778/FgFhdZGZLvPmy/ToriQ6bBBKCMKyLOdUWVdXt4SstH2h0WhgMBggl8tRUFAQU3tcukaVts6WgyX9bGUuZtyZ+HTKjnf69Rg12AAAG4uzsLNGjbM3qFFb6C2moAq2pqamuPp8U/21TqeDXq/nahyZmZncMEyi54F90d7ejptvvhmHDh1Cfn4+jh07hvPOOy8Wb5XcZCaELLHPDYbJyUm4XC5uCoYvBaW2sKFS687OTqhUKq8po2DXR896Op2O22Io1MQRf43qxo0bBc8C/EU9W1oOPtd48O6AHm0T8yAAinNkXGW8XGbHzNQkmpubE2Lvw4fdbsf4+DjGx8e9RlATeRTio7u7G9/97nexf/9+bNiwIdZvt7rIPDMzA4vFgnXr1nGe2Q0NDVwaGOq81dnZibVr14a1bN0fbDYbNBoNtFotCCHchysaiR5/jWo8WmG+52yZTIa0bNUX5+w5fDxkgN3NIlMCnFGTj3Nq8/G1ahUU6YkjNPUya2pqglQq5W5MkY5DxgK9vb24/vrr8eKLL6K+vj4eb7m6yKzVamE0GpGZmYmJiQlO2hhK0UWHE4Q0uKPprEajgcPh4LTL4QzRx3qNajig52ytVrt4TYRBt86DEZcC7w7oYbAunrO3Viixc4MaZ9eoUaqMX8pNTQ6ampqWqLPojYmm47FWcvliYGAA3/72t/HHP/6R81GPA5KbzMBi2T9c0LaTXC7Hxo0bQxa6AGB6ehrj4+Ne62SEBjXa02q1MJvNnFjC354mugkyXmtUQ4H6Y1mtVqSlpXHnbC0rx7FJO44M6DGkswIANhRmYecXba+NxbEbWqDFt3BNDnyVXOE4mEaLkZER7Nu3D88++yxaWloEfe0QWD1kdrvdaG1thcfjwfbt27l9QcEq1kNDQzCbzWhoaIhbKkZbFzSLoIb2arUaVqsVXV1dCW2t+F7riRMnIJPJOKM733O2QqGAPU2BNq0HRwYM+Hx8DiwBCrNlOHtDHnbWqLG9MhdSiTDnfb5bSTR1Cd8bayQbN0JhfHwcV111FZ566imcdlpEY8RCIPnJ7HQ6A5r6UdD+YlFREbRaLRoaGkJKM7u7u5Geno7169cnrFhCq7MajQazs7NwuVyoqKhAWVlZwt0z6doapVIZ0FbH7zk7Kxc9c2K8P2TCh4MG2Fws5FIxTq9W4ewaNc5cn4ecjOjO2dSxUyi3En/XT2+skWZpU1NTuOKKK/D4449zOv84Y+WTmXoW19fXIzMzE6Ojo9DpdJBIJH4ryw6HAx0dHSguLk7YedQX09PTmJiYwPr16znBAcMwyM/PR0FBQdz9renARGFhYUQ/o6X9bBWGbTJ8PLqAIwN66BacEDMMtlbkcNXxstzwvje6p7m5uTlmNzqr1cpdv8fjCdvBdGZmBldccQV++ctf4qyzzorJtYWBlU3miYkJTExMYNOmTVwkpj90ek7SaDQghKCgoAByuRwnT5702vOUSPDXqDY2NnqleVRaqtVq4XK5OAVarAfo6chgRUXFsrTD/Blnm80GZW4udKwcx6cc+Ee/Hie1FgDA+gI5189uKMmGyM/3Njs7i7Gxsbi2w6iDqVarDTrnrNVqcdlll+GRRx6JVf84XCQ/mV0ul9/1qf39/bBaraivrw9Z6HI4HBgeHsbU1BQyMjK4VSOJFBjw16jW1tYGbZtRLy6tVguLxQKVSoWCggIolUpBiU11zULf7Pyds53SHO6c/dnYHDyEID9LirO+iNg71iohk4gxPT2NycnFvnYiWkwAOMskOuecnp6Oubk55OXl4bbbbsMDDzyACy+8MCHXxsPKI7Pb7UZHRweysrJQVVUVUggCLBYmZmdn0djYCIZhoNPpoNFoYLPZuIgXTzP45axR9UcMKi1dTgGHyiFjXXzzd06VZavQMyfC+0Nz+GDQAKvTg4w0MVpKM1Cf48Y1ZzdDrUjcKh1fWCwWPPvss3jyySeRkZGBa665Bv/8z/+c6EGclUVmm83GpYDUIC+Uxrqvrw9ut9uvgsrj8XDEXlhYiFnE44PO/ZaUlHAjjNGCEoNKSzMyMlBQUBDxGCEVXyRCDul7zs7JzcOITYo3T8zioxEzjA4CEQNsKc9Z7GdvUKNCFdvl86EwPz+Pyy67DHfddRd27tyJv/3tb7jkkkviuoTPD5KfzNQ6iBa6Nm7ciOzs7JBEpgbwOTk5YUU/lmW5YYq5uTnk5OSgsLDQby84WsRyjSo12aPSUuoKmp+fH7QySyex/Ikv4g16zh4dHYXNZkNRcTFMjAKfTjlwZECHvtnFc/a6/MzFAlqNGo1lCr/n7FhhYWEBV1xxBb7//e9j3759Ub3GDTfcgNdeew0FBQXo6uoCsPh7uOqqqzAyMoLKykq89NJLASP9/Pw8Nm7ciL179+KJJ56gX14ZZB4fH8fo6CgaGxuXFLr8gQovKioqorLnJYTAaDRCo9FwvWA6TBFtKhvvNaq0AKjVauHxeLjKOF9aSivE0UxixQojIyOYm5vDxo0bOWtfepxwy3LQpmVxZECP46OL5+w8uRRn1yyOcf5TVS7S02Jn5WS1WnHllVfiuuuuw3XXXRf167z33nvIysrCtddey5H5hz/8IVQqFe655x48+uijMBqN+NnPfub3+XfccQe0Wi1UKtXKIvPJkyeh1WpRX1/PiUCCEZmSRqizH78XrNfro0pldTodTp48mbDoRyOeRqPhNmwwDAOj0ZgUAxMUw8PDnIiHnw35O2enK1Tomxfj/aE5vDegh8XpQUaaCF+pUuHsDYv97Dy5cC0sm82Gq6++GldeeSVuvPHGZb/eyMgILr74Yo7MGzZswJEjR1BcXIzp6WmcddZZ6OvrW/K81tZW/OIXv8CuXbtw/PjxlUVmq9UKQkhYha7Z2VmMjIygsbExJqThp7JarZabkiooKAjY+6ROnnQYINHweDzo7e2FXq+HRCJBbm4uCgoKBD1ORAqqxrPZbGFNh1mtVu44Qc/ZY3YZPh5bwDv9eszMO8AA2Fyeg7Nr8rBzQz7WqqM/ZzscDlxzzTW4+OKLccsttwhST/Els1KphMlkArD488jNzeX+TcGyLHbu3Innn38eb731VkzIHNN+wYMPPoiSkhLs3r074DQTIQTDw8OYm5uL2Dc7EjAMw+0Iqqqq4j5U7e3tYBiGI3Z6ejr3AV1YWAi4ZD3eoD8nlmU5lZLJZIJGo0F/fz+ysrK440S82kDU8dTpdHJtxlDIzMxEZWUlKisruawj36LFeSobrrwgDyYmC59OOfBOvx6/fHsIv3x7CJV5Gdj5RdurKYL1P06nE9dddx3OP/98wYgcCoGC1m9+8xtceOGFMRU7xTQyDw4O4sCBA3j55ZchlUpxySWXYO/evSgqKgLDMJx+WCKRoKamJmHRxdc/jGVZyOVyNDQ0JMVsLe1rA0Btbe2Sa6ImdfQ4IZVKuXN2rDIKqhdgWdbvNUUKf/PZJEOJNi2Ld08acGzEBDdLsG9rKe6/sCbk67lcLtxwww047bTTcPfddwv6e4wmzb7mmmvw/vvvQyQSYWFhAU6nE7feeiseffRRYCWk2dyLfGHlc/DgQfz1r38Fy7I455xz8M477+Cxxx5DbW2tEG+zbFApJN3853Q646beCgTqDx7JbDTNOugIJCW2UBZB9OYiEolQU1Mj+M/Ft9Yhk8mQmaNC37wEa/Ozsak0eBHS7XbjpptuQn19Pe67776YTFfxyXz33XcjLy+PK4AZDAb8/Oc/D/j85557buWdmf2+ICF4//33cc0116CqqgoOhwMXX3wx9uzZg6qqqoRFQn9rVN1uN1d8slgsyMvLi6sxIF1fSt1TooHD4eAq4/TmFI5mORDoWKVUKuWmsWINqrumxhG+u7H48Hg8uPXWW1FZWYl/+7d/E/z69u3bhyNHjkCn06GwsBAPPvgg9u7diyuvvBJjY2OoqKjASy+9BJVKhePHj+O3v/0t/vu//9vrNVYNmQHgySefxNe+9jU0NDRAo9HgL3/5Cw4ePAij0YgLL7wQe/bsieueI2o5G0wK6WuHS4tPsTCBB740OSgpKRHM2pXenKhmmTpnhvs90GNRRkZGwm68vrpxvvsnANx5551QqVR49NFHE77uJwKsXDIHgl6vx8svv4yDBw9iZmYGX//617F3796YeGhRRLNGlep9aS9bKFkmBc0SwtmxFC18hTb+dk/7Pr6rq4srICYD+Ofsf/3Xf8Xs7CyKiorw4osvJpXLZxhYfWTmw2Qy4dVXX8XBgwcxMjKC8847D3v37kVTU5NgxKZrVJfTDuPLMvV6Pef4qVaro6oq04GJaJatRwvfXjDd35yfn88tF6eKvEDz0YkEy7L4yU9+Ao1Gg6qqKvzjH//A22+/vdwVq/HE6iYzH2azGYcPH8bBgwfR19eHnTt3Yu/evdi6dWvUxI7FGlVfx0+pVIrCwsKwtyjSjRfxUpr5A+3H0zMqXWxeWFiIdevWJeSagoFlWTz44IMwGAx4+umno/o8RCvRbGtrwy233IL5+XmIxWLce++9uOqqq6L5Nk4dMvNhtVrxxhtv4ODBg+jo6MAZZ5yBvXv3YseOHWGluLQ3arPZlqiVhAZfpEL11gUFBX7llzTdb2xsTJqI4vF48Nlnn0Emk8HlcnFrf6i0NNFtO0IIHn74YYyNjeG5556L+ogTrUSzv78fDMNg/fr1mJqaQktLC3p6eqJZrXNqkpkPu92ON998EwcOHMDx48fx1a9+FZdeeim++tWv+k1xaQEnLS0tJi2VUNdKic2yrFe7SK/XY2BgICkGJihom664uJgrwNHZbDqCSifVEuF1TQjBf/zHf6C3txd//OMfly2UiVaiyUdTUxMOHDiA9evXR/r2KTLz4XQ68c477+DAgQP46KOPcNppp2HPnj0488wzIZVK4XA40N3djby8vKjbPEJeq0ajgUajgdVqBcuy2LRpU0zHNyOB2+1GW1sbSktLuTadL3yr+/Fc+0MIweOPP45PP/0Uf/7znwU5JkUj0eTj2LFjuO6669Dd3R3N958icyC43W6899572L9/P95//33U1taiq6sLzzzzDDZv3pzoy+MwMTGB6elplJSUcHuPaS87Flay4YDuvl6zZk3Y1kOB1v5EWwQMBkIInnrqKRw5cgQHDhwQTOEWjMwAkJubC6PR6Pe5NHL//ve/x44dO6J5++TXZicKEokEO3fuxM6dO9HZ2YlvfOMbaGlpwfe//33U19djz549OO+882K6NC0UqH8Y1X6XlpbC4/FAr9djfHyc8+guLCyMW8SmHmKRtsREIhFyc3ORm5vrVQQcHR0VdO0PIQTPPPMM3nrrLRw6dCimwy+FhYWYnp7m0uxAP4/5+XlcdNFFeOihh6IlsmBYlZGZjzfffBOlpaXYuHEjWJbFsWPHsH//frz55ptYv3499u7di69//etxc+mgBTiHwxG0f049umdnZznDhVimsU6nE59//nlES9fDga96azlrf/7whz/g0KFDePnllwWvLUQj0XQ6nbjgggtwySWX4M4771zO26fS7OWAZVl8/vnnOHDgAN544w2Ul5djz549uPDCC2Pmo0UIQU9PD8RicUQFOEIIJ1IxGAzchJRarRZMpNLW1obq6uqYOp8uZ+3Piy++iD/96U949dVXBa/2RyvRfP755/Gd73zHax/Vc889h+bm5kgvIUVmoUAIQVdXF/bv34/XX38d+fn52LNnDy6++GLBhBtUeJGdnR2xEaDvtdIJKZ1OF7V3GIXdbkdbWxs2bNgQV1M7up2CergFW/tz4MABPPPMMzh8+HDSrX0VCCkyxwJ0IujAgQN47bXXoFAosHv3blxyySXIz8+PioTUlTQ/Px/l5eWCXi9fpEKXB4Q7+kjVZrW1tdH0RgUDPVJoNBqYTCZkZ2cjJycHOTk5OHLkCJ544gkcPnw4Kdb+xAgpMscahBBuJvuVV16BTCbDJZdcgj179nAz2aFAq8NlZWUB2zxCga6l1Wg0nOFCfn6+3/MlXQuULLuxKOj44yeffIIf/OAHWFhYwH333Yd9+/YJbqSYREiROZ4ghGB0dJSbySaEcGYLZWVlfoltt9u5hfHx3grpcDg4YvuaAlKf7UiGS+KNt956Cw899BAee+wxvPfee1AoFLjlllsSfVmxQnKT+Y033sAdd9wBj8eD733ve7jnnnuifamkAyEEU1NTOHjwIP7yl7/AbrdzM9n0PEzHKuN9FvUHl8vFFZ6sViucTifq6upQUFCQFCIVX7z77ru4//778frrr0c9NbYcS9zf//73+Pd//3cAwH333bcsN88wkbxk9ng8qKmpwZtvvomysjJs27YNL774IjZu3BjdVSYxCCHcTPahQ4dgMpmwY8cOHD16FH/961+TKoU1m83o7OxEaWkp5ufnuZnmwsLChEgy/eGDDz7Aj3/8Y7z22mvLOpZEq7c2GAzYunUrjh8/DoZh0NLSgtbW1ljfkJOXzB9//DEeeOAB/P3vfwcAPPLIIwCAH//4x9G83IrC66+/jltuuQX19fWYnZ3Frl27sHfvXtTV1SV0WJ7aGDc1NXFiGZZluYry/Pw8lEplQt0+jx49iv/3//4fXn31VUGM76LRW7/44os4cuQInnrqKQDAzTffjLPOOitq0/wwkbwKsMnJSa+qbVlZGY4ePRqLt0o6DAwM4KOPPkJpaSlMJhNeeeUVPPTQQxgdHcW5556LSy+9FI2NjXEli8lkQm9vL5qbm72KYSKRiBNx8A0X+vv7BVkeEAlaW1tx11134eWXX46Zg+Xs7CwX7YuKijA7O7vkMf4+u5OTkzG5HqGxKuWcicQdd9zB/V2pVOLaa6/Ftddei/n5eRw+fBj/9V//hb6+PpxzzjnYu3cvWlpaYkpso9GIvr4+NDc3B111IxKJoFKpoFKpvAz1hoaGkJmZGTOtNQC0t7fjtttuw8GDB+M2BBPKx30lIiZkLi0txfj4OPfviYmJZS9aW+lQKBTYt28f9u3bB6vVitdffx1PPfUUOjs7ceaZZ2Lv3r047bTTBI2Cer0eJ0+exObNmyPSRTMMw/V5q6urvbTWUqmUa3kJoY3u7u7GzTffjP3796O6unrZrxcM4eitS0tLceTIEe7fExMTiVzCHhFicmZ2u92oqanB22+/jdLSUmzbtg0vvPCCl+wthUXQmez9+/ejtbUVp59+Oi699FJ85StfWVYU1Ol0GBwcxObNmwUdSODb+NI0nS4PiBS9vb24/vrr8eKLL8bksxGN3tpgMKClpQWfffYZAGDLli1obW2NtYVT8hbAgMVC0J133gmPx4MbbrgB9957b7QvdcrA6XTiH//4Bw4cOICPP/4YO3bswJ49e3DGGWdEREitVovh4WE0NzfHdLKIvzzA13AhFAYGBvDtb38bzz//PBobGwW/tuVY4j777LN4+OGHAQD33nsvvvOd7wh+fT5IbjKnsDy43W68++672L9/Pz744ANs2bIFe/bswc6dO4OmzLOzsxgbG4v7Ujn+EEWo5QEjIyPYt28f/ud//gdbtmyJ2zUmMU5tMldWViI7OxtisRgSiQTHjx+PaE/uSoLH48EHH3yAgwcP4p133kF9fT327t2Lc8891ysKTk9PY3JyEs3NzXHbN+UP/OUBVquVM1zIzs7G5OQkrrrqKvzud7/D9u3bE3aNSYYUmY8fP+6l141kT+5KBcuyOHr0KA4cOMDNZF966aUYHx9Hfn4+rrzyyoQS2RfUcGFgYAA33XQTGIbBXXfdhVtuuSWprjPBSJHZl8zRmLCtZLAsi88++wz3338/Ojo60NLSgr179+KCCy5IKuUZAMzMzOCKK67Avn37MDg4CKfTiWeeeSbRl5UsSF7RSDzAMAzOP/98MAyDm2++GTfddFNYooDVBJFIhMzMTKSlpWFgYAAnT57E/v37cfHFF6OwsBB79uzBRRddFDcz/UDQaDS44oor8LOf/QznnnuuIK/5q1/9Ck8//TQIIbjxxhuXOH3Mzc3hW9/6FsbGxuB2u/GDH/wgHoWshGLFRubJyUmUlpZCo9HgvPPOw+OPP47du3eHbcK2msCyrJfwhDqa0JlspVLJzWSr1eq4iiX0ej0uu+wyPPjgg7jgggsEec2uri5cffXVOHbsGKRSKXbt2oXf/va3Xn3qhx9+GHNzc/jZz34GrVaLDRs2YGZmJqbV/WVAkF/Iitms5QsqQikoKMCll16KY8eOcaIAAEFN2FYbfBVkDMNg48aN+Jd/+RccPXoUv/71r2E2m7Fv3z5ccskleOqppzAzM4MQN/Jlw2g04oorrsB9990nGJEBoKenB6eddhoyMzMhkUhw5pln4tChQ16PYRgGZrOZMxhUqVSr/oy+IslssVhgNpu5v//f//0fGhoasHv3bvz+978HsDjGtmfPnkReZlKAblz48Y9/jA8//BDPPvssWJbFtddei127duGJJ57AxMSE4MSen5/HlVdeibvvvhu7d+8W9LUbGhrw/vvvQ6/Xc2o6vuIQAG677Tb09PSgpKQEmzZtwq9+9auVtBUyKqzINHtoaAiXXnopgMU2yDe/+U3ce++90Ov1fkUBKSwFfyb70KFDcDgcnItKZWXlslLxhYUFXH755bjllltiNm30zDPP4De/+Q3kcjnq6+shk8nw2GOPcf9/4MABfPjhh/jP//xPDA4O4rzzzkN7e3vCdniFwKldzU5BOBBCMDs7y81kz83N4aKLLsKePXuwfv36iIhtsVhw1VVX4frrr8e1114bw6v+Ej/5yU9QVlaGW2+9lfvaRRddhHvuuQdf+9rXAAA7d+7Eo48+mqy97RSZU4gNdDod/vrXv+LgwYPQarW44IILsGfPHtTV1QUlts1mw9VXX42rrroK3/ve92J6jRqNBgUFBRgbG8P555+PTz75xMuU8JZbbkFhYSEeeOABzM7OYsuWLWhvb09WHzFhKpKEkGB/Tll85zvfIfn5+aS+vp77ml6vJ+eeey6prq4m5557LjEYDIQQQliWJbfffjtZt24d2bRpE2ltbU3UZQsOg8FAnnvuObJ7927S3NxM7r77bvLRRx8Rs9lMLBYL98dgMJBdu3aRX//614Rl2Zhf1+mnn07q6upIY2MjeeuttwghhDz55JPkySefJIQQMjk5Sc477zzS0NBA6uvryR//+MeYX9MyEIqHYf1JkTkA3n33XdLa2upF5rvvvps88sgjhBBCHnnkEfLDH/6QEELI4cOHya5duwjLsuTjjz8m27dvT8g1xxpzc3PkT3/6E7nssstIY2Mjueuuu8i7775LDAYDufjii8ljjz0WFyKvQqTIHGsMDw97kbmmpoZMTU0RQgiZmpoiNTU1hBBCbrrpJvLCCy/4fdxqhcViIfv37ydXX301UavV5Cc/+UmKyNFDEDKv7sabwAikMAtkNRNrn+xEIjMzE5dffjkuv/xyGAyGpDEEPJWRInOUWI22M9Ei1f5LDqzuLrrACKQwS9kkCYNf/epXaGhoQH19vVfPmI8jR46gubkZ9fX1OPPMM+N7gUmOFJkjQCCF2e7du/GHP/wBhBB88sknyMnJWdUpdizQ1dWFp59+GseOHUN7eztee+01nDx50usxJpMJt956K1555RV0d3dj//79Cbra5ESKzAGwb98+/NM//RP6+vpQVlaGZ555Bvfccw83Q/zWW29xWzouvPBCVFVVobq6GjfeeCN+85vfJPjqVx7C0Vu/8MIL+MY3voE1a9YAwCmjvQ8XKdFICkmBnp4e7NmzBx9//DEyMjJwzjnnYOvWrXj88ce5x9x5551wuVzo7u6G2WzGHXfcETeVWYxxak9NrUbccMMNKCgoQENDA/e1Bx54AKWlpWhubkZzczNef/117v8eeeQRVFdXY8OGDdz2kJWKuro6/OhHP8L555+PXbt2obm5eYntsNvtRmtrKw4fPoy///3v+OlPf4r+/v4EXXHyIUXmJML111+PN954Y8nX77rrLrS1taGtrQ0XXnghAODEiRP485//jO7ubrzxxhu49dZb4fF44n3JguK73/0uWltb8d577yE3Nxc1NTVe/19WVoavf/3rkMvlUKvVOOOMM9De3p6gq00+pMicRDjjjDPCbvO8/PLLuPrqqyGTybB27VpUV1fj2LFjMb7C2EKj0QAAxsbGcOjQIXzzm9/0+v89e/bggw8+gNvthtVqxdGjR1FXV5eIS01KpPrMKwBPPPEE/vCHP2Dr1q345S9/idzcXExOTmLHjh3cY1bSTqRAuOyyy6DX65GWloZf//rXUCqV+O1vfwsA+P73v4+6ujrs2rWL29X1ve99z+tIcqojVQBLMvhuYZidneWsfu6//35MT0/j2WefxW233YYdO3bgW9/6FoDFFPWCCy7A5ZdfnsjLTyE6pApgpwIKCwshFoshEolw4403cql0SqiSgi9SZE5yUMUZAPzlL3/h0srdu3fjz3/+MxwOB4aHhzEwMJCsg/cpxAmh0uwU4giGYV4EcBYANYBZAP/6xb+bsXjkGQFwMyFk+ovH3wvgBgBuAHcSQv4W72tOIXmQInMKKawSpNLsFFJYJUiROYUUVglSZE4hhVWCFJlTSGGVIEXmFFJYJUiROYUUVglSZE4hhVWCFJlTSGGV4P8HTxgNgTvQWSIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot(nx[1],nz[1],V[1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "9e7dd8c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh_hold=1\n",
    "def vel_to_temp_2D(nx,nz,Vs,Table):\n",
    "    Temperature_out = np.zeros_like(nx)\n",
    "    Density_out     = np.zeros_like(nx) #+3300\n",
    "    diff_Vs         = np.zeros_like(nx)\n",
    "    P_out           = np.zeros_like(nx)\n",
    "    melt_out        = np.zeros_like(nx)\n",
    "    #Vp_out          = []#np.zeros_like(tomo[:,1])\n",
    "    #Vs_out          = []#np.zeros_like(tomo[:,1])\n",
    "    m,p=np.shape(nx)\n",
    "    # looping over all the geospatial nodes\n",
    "    for i in range(m):\n",
    "        #P_func=lithostatic_pressure(nz[i][:],Density_out[i][:])\n",
    "        P_func=lib.pressure_inter\n",
    "        threshhold = 12000\n",
    "        while threshhold >1000:\n",
    "            for k in range(p):\n",
    "                P  = P_func(nz[i][k])\n",
    "                Vs_in = Vs[i][k]\n",
    "                P_table,temp,dens,vp,vs=lib.lookup_vs_P_accurate(Vs_in,P.tolist(),Table)\n",
    "                Temperature_out[i][k]=temp\n",
    "                Density_out[i][k]=dens\n",
    "\n",
    "            P_func_new=lithostatic_pressure(nz[i][:],Density_out[i][:])\n",
    "            threshhold= np.sqrt(((P_func_new(nz[i][:]) - P_func(nz[i][:])) ** 2).mean())\n",
    "            #thresh_hold = abs(P_func_new(nz[i][j][:]) - P_func(nz[i][j][-1]))\n",
    "            print(threshhold)\n",
    "            P_func = P_func_new\n",
    "            #Temperature_out[i][j][k]=temp\n",
    "            #Temperature_out[i][j][k]=temp\n",
    "    return Temperature_out,Density_out "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "5a331035",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\kumar\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\ipykernel_launcher.py:18: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "c:\\users\\kumar\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\ipykernel_launcher.py:18: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1274.0000000000002\n",
      "0.0\n",
      "2861.6000000000004\n",
      "0.0\n",
      "4292.4\n",
      "0.0\n",
      "7115.8293920000015\n",
      "0.0\n",
      "8848.592872000001\n",
      "0.0\n",
      "10581.356352000003\n",
      "0.0\n",
      "12314.119832\n",
      "0.0\n",
      "14046.883312000002\n",
      "0.0\n",
      "15755.143486486488\n",
      "0.0\n",
      "17426.648702702707\n",
      "0.0\n",
      "19098.15391891892\n",
      "0.0\n",
      "20769.659135135138\n",
      "0.0\n",
      "22441.16435135135\n",
      "0.0\n",
      "24112.669567567566\n",
      "0.0\n",
      "25784.17478378379\n",
      "0.0\n",
      "27455.680000000008\n",
      "0.0\n",
      "29061.116\n",
      "0.0\n",
      "30666.552\n",
      "0.0\n",
      "32271.98800000001\n",
      "0.0\n",
      "33877.42400000001\n",
      "0.0\n",
      "35482.860000000015\n",
      "0.0\n",
      "37088.296\n",
      "0.0\n",
      "38693.731999999996\n",
      "0.0\n",
      "40299.168000000005\n",
      "0.0\n",
      "41878.22566666667\n",
      "0.0\n",
      "43457.28333333333\n",
      "0.0\n",
      "45036.34100000002\n",
      "0.0\n",
      "46615.39866666667\n",
      "0.0\n",
      "48194.45633333336\n",
      "0.0\n",
      "49773.51400000002\n",
      "0.0\n",
      "51352.57166666668\n",
      "0.0\n",
      "52931.62933333334\n",
      "0.0\n",
      "54510.68700000002\n",
      "0.0\n",
      "56055.51000000002\n",
      "0.0\n",
      "57600.33300000001\n",
      "0.0\n",
      "59145.15600000002\n",
      "0.0\n",
      "60689.979000000014\n",
      "0.0\n",
      "62234.80200000001\n",
      "0.0\n",
      "63779.625000000015\n",
      "0.0\n",
      "65324.448000000004\n",
      "0.0\n",
      "66869.27100000001\n",
      "0.0\n",
      "68414.09400000003\n",
      "0.0\n",
      "70150.01699999999\n",
      "0.0\n",
      "71885.94\n",
      "0.0\n",
      "73621.863\n",
      "0.0\n",
      "75357.78600000002\n",
      "0.0\n",
      "77093.709\n",
      "0.0\n",
      "78829.63200000003\n",
      "0.0\n",
      "80565.55500000002\n",
      "0.0\n",
      "82301.47800000002\n",
      "0.0\n",
      "84037.40100000003\n",
      "0.0\n",
      "85773.32400000001\n",
      "0.0\n",
      "87558.6096\n",
      "0.0\n",
      "89343.89520000003\n",
      "0.0\n",
      "91129.18080000002\n",
      "0.0\n",
      "92914.4664\n",
      "0.0\n",
      "94699.75200000002\n",
      "0.0\n",
      "96485.0376\n",
      "0.0\n",
      "98270.32320000001\n",
      "0.0\n",
      "100055.6088\n",
      "0.0\n",
      "101840.8944\n",
      "0.0\n",
      "103626.18\n",
      "0.0\n",
      "105462.32759999998\n",
      "0.0\n",
      "107298.4752\n",
      "0.0\n",
      "109134.62280000003\n",
      "0.0\n",
      "110970.7704\n",
      "0.0\n",
      "112806.91800000003\n",
      "0.0\n",
      "114643.0656\n",
      "0.0\n",
      "116479.21320000001\n",
      "0.0\n",
      "118315.36080000004\n",
      "0.0\n",
      "120151.50840000002\n",
      "0.0\n",
      "121987.65600000005\n",
      "0.0\n",
      "123879.21280000005\n",
      "0.0\n",
      "125770.7696\n",
      "0.0\n",
      "127662.3264\n",
      "0.0\n",
      "129553.88320000001\n",
      "0.0\n",
      "131445.44000000003\n",
      "0.0\n",
      "133336.9968\n",
      "0.0\n",
      "135228.5536\n",
      "0.0\n",
      "137120.11039999998\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "T,D=vel_to_temp_2D(nx,nz,V,DMM_atten_melt_corrected)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "2156d897",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAADtCAYAAABH/TtHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABsQElEQVR4nO2deXxcZdn+v2dmsu97s7VJm6ZbmnSnUChKKUUoXQRZXBARFRXeVxEQfvqK+oqg8Aq+riggvsoipWVHoCCIQLEF2mxNky5p0ixNJvs+6/P7Y/IcTiaznElmkpTO9fn003Zy5syZyTzXuZ/7vu7rVoQQhBFGGGH4gmG6LyCMMMKY+QgTRRhhhOEXYaIII4ww/CJMFGGEEYZfhIkijDDC8AuTn5+HSyJhhBF6KNN9Af4QjijCCCMMvwgTRRhhhOEXYaIII4ww/CJMFGGEEYZfhIkijDDC8IswUYQRRhh+ESaKMMIIwy/CRBFGGGH4RZgowggjDL8IE0UYYYThF2GiCCOMMPwiTBRhhBGGX4SJIowwwvCLMFGEEUYYfhEmijDCCMMvwkQxDRBCYLPZcDgchF3QwzgV4M+4Jowgw+l0YrPZGB4eBkBRFEwmExERERiNRgwGA4oy431MwjjNoPi5o4Vvd0GCEAKHw4HNZgPAbrerjzudTvU4SRJRUVFERESEieP0wIz/BYcjiimAdqvhvugVRcFoNI45trKyksLCQmJjY1EUhYiICEwmEyaTCUVRwsQRxpQjTBQhhtPppKWlBbvdzqxZs1AUBSEEQgiPC15RFAwGA0ajEaPRiBACq9WKxWIBwGAwEBERoW5VwsQRxlQgTBQhghACu92O3W7HarVis9kmtKC1EYfcJlqtVqxWK+AiDvccRxhhBBthoggBZBTgdDrVO742F9TZ2UltbS2xsbGkpKSQkpJCTEyM+nP347WPA2HiCGPKESaKIEMmLOXWQksUTqeTI0eO0Nvby9KlS7FarfT09FBXV4fFYiEhIYGUlBQcDoeu1woTRxhThXDVI0jQbjXc8wYnT56kp6eH3t5e0tPTmTt37hgyAVcuY2BggK6uLpqamjCZTCQnJ6sRR2Rk5ISuSf6RkDkOk8kUJo6ZgxmfZApHFEGA1EZotxpa9Pb20tLSwvLly0lJSfF4DoPBQGJiIomJiQwNDZGXl4fT6aS7u5vm5mYcDgdJSUmkpKSQnJxMRESE3+tyvxZJGhaLRU2O9vb2kpWVhdFoVKsqYYThjjBRTALu2gj3helwOKitraW3t5ecnByvJOEOeZ7k5GSSk5MpLCzE4XDQ29tLd3c3jY2NCCHGEIfJ5P9X6Yk4jhw5QlJSkvqY0WgcE3GEiSMMCBPFhCGEoKuri8jISCIiIsYtqMHBQSoqKsjJySE9PZ2enp5JvZ7RaCQ1NZXU1FTAJdiSxHH8+HGVWFJSUkhKShqjzfAGSRzaHIfT6WRkZGTM64aJI4wwUUwATqcTq9XKsWPHmDNnzrj8QXNzMw0NDZSUlJCYmEhHR0dAPR3eqh5amEwm0tLSSEtLA8Bms9HT00NHRwdHjx7FaDSOIQ49uQhPEUeYOMKAMFEEBO1WQwqjtLDb7Rw8eBCANWvWqNsBPQtfi4ksvoiICDIyMsjIyABQKyrt7e0cPnyYiIgINTGakJAwKeIYHh4eU3EJE8fHH2Gi0AlP2giDwaD2afT19VFVVcWcOXPIzc0d89xAiUK+3mQQGRlJZmYmmZmZAFgsFrq7u2lpaaG/v5+oqCi1FOtNJeoO7fuW1xgmjtMDYaLQAbnV0GojwLVwnE4nDQ0NtLS0UFpaSnx8/LjnTySiCHb7eVRUFLNmzWLWrFkADA8P09PTg9VqZe/evcTExKgRR1xcXNCIQ/aohInj1EaYKHzAXRvhHq47nU4OHz5MYmIia9as8ZpADMXCnyxiYmKIiYmhqamJVatWMTw8rCZGBwcHiYuLUysqsjnNHzwRh8PhUDtl+/r6iIuLIz4+HpPJFO6MPYUQJgov8KeN6O7upq2tjby8PIqLi32eayZEFP5eLzY2ltjYWHJzcxFCMDQ0RHd3N8eOHWNoaIj4+HiPcnN/59V+bu3t7WRmZqqEKr045J8wccxchInCDf60EUIIjh07RkdHB7NmzSI5OdnvOWdiROELiqIQFxdHXFwceXl5CCEYGBigu7t7jNxcVlWio6N1nVcIoXbFyv/b7fYxn3WYOGYmwkShgbyTtrS0kJ+fP+5LarFYqKysJCEhgdWrV3Ps2DFdBOCJKHyRx0wjFkVRSEhIICEhgdmzZ4+Rmx86dAir1UpiYqJfubl70tQTCYeJY2YiTBSjkFsNq9VKR0cHs2fPHvPzjo4OamtrWbBgAenp6YBLdj1RojiVoZWbg+uz6+vr8ys391dd8UQcNpttHHGEbQOnHqc9UbhrI6RZjIS243PVqlVERUWpP5NVD3+Y6TmKycJgMOiSm4+MjOjujAXP7l/uxBF2/5oanNZE4U8bMTw8TEVFBenp6axatcqjjd3pGFH4gze5ufThcFeN6pGbg2fiCLt/TQ1OW6Lwp41oa2vjyJEjLF682GfH50QjCnln9WaH93EiFik3b25uZv78+ZhMpknLzcE3cSiKwvDwsFrNCbfUTw6nHVH400YIIRgcHKSlpYXVq1f79IGYaERx8uRJDh8+DKDqFVJTU3WXHU9VCCHUu36w5eYw3jawubmZzMxM9bMPm/hMHKcVUXjaamgxMDBAZWUlBoOBZcuW+Q1bA81ROJ1OamtrGR4eZuXKlRgMBlWvcPjwYUZGRkhISMBut08Jabib2kzF63n6TPXKzZOTk0lISNAt/nI6naoiNOz+NTmcNkThXnZz/7JpOz6rqqp0fRkNBoN6Tl9QFAW73c6+ffvIzMxk4cKF2O12nE4n8fHxxMfHk5+fj9PppL+/n6NHj9LY2EhTU5Makuv1nAgUU7mH19tT4k1u3tTURH9/v265udPpVBe/thdFXou8cWiJI+z+5Rkfe6KQW43KykoWL1487kvlreNTD/RuPTo7O+nv72fVqlU+zWsMBgNJSUmkpqYSGxtLWloaPT09qrTaYDCoCyQxMfGU+yLrJQp3SLl5dnY2QgiPcnNJqFq5uZYo3OGpFOvu/qUljtPd/etjTRRaGXZvb++4X7Svjk898EcUshekr69PlUAHcl6j0TjOc0JKx+vq6tSQPDU1VXcjlxZTnTCdKFFoEYjc3GazBZTf8EccsjP2dLQN/FgShbsM21PCsrGx0WfHpx74qnqMjIxQUVFBWloay5cv5/3335/Qa2gRERExZi8v76wNDQ0MDAwQFxdHampqwP0YUwWZGwomfMnNBwcHOXDgAImJiQHLzcMmPmPxsSMK9/F97r88q9VKVVUV0dHRPjs+9cBbRNHZ2cmhQ4dYuHAhaWlpOJ3Occf5Ihm9WxoZkufk5KjVGm0/hh5Z9VQj1ItJKzc3m82UlZUxNDQUsNzc03lPZ+L4WBGFN22ERFdXFzU1NRQVFZGVlTXp13Nf0EIIjh49Snd39xgVp6eFH+ywX1EUj4nRrq4umpubcTqdav4jOTl5UgQ52eucKjidToxG44Tk5v7gjzi+//3v861vfYuioqLgv7FpwMeCKPRoI6xWK4cPH2bFihVBKz1qy6NWq5WKigqSkpLU0udkzjtZIpGJ0aSkJFVWLROj9fX16s8dDofPpN+pDndi0is3D7TS5E4c9fX1IalSTRdO+XfiTxthsVioqKgAYNmyZWN6Nfyd19/dTzaFdXd3c/DgQYqLi1URkRYzIQR1T4xarVa6urpobW3l/fffJyoqSs1vTCQxeqoiFO7mgJpU/bjglCYKT+P7tNB2fB4/flz3eeUd3ddicTgF/7GrjpHhIUrSO/n0WYtJT08N6Pq9vcZUSLgjIyPJyMigubmZlStXjis5TsSo5uOAYLmby8/w44JTkijkVqOmpobi4mKvFnVSuxAVFUVjY6MuFSXoax/vGxphb2M/UUaoMDt4rGY/mQlRrJuXylnzUjlrbirp8fqil+mC9j16Sox2dXWNSYzK/MZMSYxOBSYqN7darbqj11MBpxxRaLURnZ2d4+7IQ0NDVFZWkpGRwcqVK9Wf623ggo9yD97CzN7eXvbtrwTgqiXxXL2hjHePdvHO0S7eqOvg6QOtABRnxbNubirritKwO2Zmm7m3iEYmRqVRjUwANjU14XQ6xyhGpysxOh3QIzd/4403iIiIwOFweM1TXHvttbzwwgtkZmZSVVUFgKIo9wCXAFbgKPAlIUTP6M9uB74MOID/EEK8Mvr4hcAvASPwoBDi7lC871OGKDxZ1Lnj5MmTHD161GPHZyBE4e1Yqb9obW1lSclSeGMfigK5yTF8ZmUun1mZi9MpOHiyf5Q4Ovnr3hP8aU8jJgVWHf2AdfPSOGteKsUZp044754AlPv4rq4uNTEqhV+BNHF9HOBJbv7OO+/Q3t7OihUrWLJkCY899ti47+s111zDDTfcwNVXX619eDdwuxDCrijKz4Dbge8qirIYuBJYAuQArymKIo1afwNsBJqAfYqiPCeEOBjs93lKEIU/bYSc8WmxWFizZo3HElegROF+R7fb7VRVVREZGcnq1avpHnY5Syvj9BEKJTmJlOQk8tVzChi2Oni/sYcd/6rk+JCN/3ntCP/zGiTHmFhTkMxZhSmcOTeFnKSPhEB6m80mi4lGLe77eKvVSnd3N62trdTW1s6IxOhUfH6eEBMTw9e//nUef/xx9u/fT0tLi8f3v379+nF5MyHEq5r/vgdcNvrvrcATQggLUK8oyhFgzejPjgghjgEoivLE6LGnH1H400bIjs/c3FwWLVrk9Us5ka2HhJR6FxYWkp2dDYDDOUoUiu/FFhNp5JyiNIztkZx11lo6Biy8e6yLtw938O6xbl6t6QBgdko0Z85N4czCFGZHOYgxTc3iCsYijoyMJCsrS9WmDA8P09XVNS4xOpWLV7a0TwccDodq0zeR1oBRXAv8bfTfubiIQ6Jp9DGAE26PnzHRF/SFGUsU/sb3ybJoRUWFOuPTFyay9ZCeBidOnBgn9XZKj4MA31d6fBRbSrO5aHEGNpuN+q4R9tR3s+dYD89VtPG3D1oxKLAgPYpzF9g5c24KS3MSiDCeOuF8TEwMubm5ai+GlFSPjIywd+9eVaeQkpKiW+AUKKZTGzLZ0qiiKN8D7MCjQbuoSWJGEoU/bYTs+HQ4HKxatUpXFj5QopAdp4qieJR6O5wykgg8fBdC0NTUxPDwMOlpaXxuVQ5fWJOHzeHkQFMfuytP8EHzEH94p5Hfv91IXKSR1XOSWFvoijgK02JOGZ2DlFTHx8fT1tbGypUr6evro6urK6SJUV/J6FBjYGCA2NjYCT1XUZRrgM3ABvHR3rAZyNccljf6GD4eDypmHFH422poOz5HRkZ03zUCIQqHw0FlZSWFhYXk5eV5vk41ogiMKCQBRUVFER8fz8mTJ6mrqyM6OprU1FSWZKSStzKdz5faSMrIZm9DD3uO9bCnvps3D3cBMCsxirUFyZw5N4W1Bcmkxk2sXDkdpjXaxCi4Po+enh66uro4duwYRqMxKInR6YwoJqqhGK1g3AqcK4QY0vzoOeAxRVF+gSuZOR/YCyjAfEVRCnERxJXAZ/28hkkIYQ/02mYMUQgh6O/vV4fL+Or4LCsrIy4ujtbWVp8lKC2MRqMuomhpaaGzs5OFCxf63F/KiCKQ+/rAwAAVFRUUFhaSmZmJzWZj1qxZqsdCV1cXR44cob+/n8jISKKiovjEvFQ2LnTV8E90D7Onvpv36nt4va6TZyraAFiYFafmN1bkJxFl0r9Apioy8SYuM5lMpKenqyMQZGJUlhujo6MDnokK00sUAwMDfoniqquu4s0336SjowNFUZqAO3BVOaKA3aPv8z0hxPVCiGpFUZ7ElaS0A98UQjgAFEW5AXgFV3n0YSFEtafXUxQlHzgfsAF/VRQlDTAKIdr1vKcZQRRSG9HV1cXg4OC4fIO3jk+j0ajb/l1RFJ/HOhwODh06hM1mIzs72287srrz0BlR2Gw2KioqWLp0KQkJCWOuReuxkJeXx8mTJ+nt7WVoaIimpiaEEOpd9rJls7h8RQ4Op+Bgaz976l3Rxl/+3cyf9jQRZTKwIj+RM0e3KcVZcRhmwDZFrxeFNjHqyaQmPj5eraj4+h1Nd0QRFxfn85jHH39c+18Ztj7k7XghxJ3AnR4efwl4yddrKYqSBXwfWISLUP4KnAlcD2xWFMUghPB5F51WonDXRphMpnGLWfZReOr4DIY2AlzJp/LycnJzc8nPz+fw4cN+w3K7U9/Ww+l0UldXh81mY926dbqiH4PBQFRUFAUFBapuobu7W1UDSsOawtRUStbl89WzZzNkdfB+Y+9oxNHNL/5RD9STGhvBGaPblDMLU5iV+JFaMBhGMnoxkdfyZFIjE6OyZdxbYnSmE8VUQFEUZTTPMR9IBzbwEam0APrMOZhGonDv+FRGHZTlYhaaGZ/eOj71bifksZ78LaVIq6SkhKSkJECfjsEpicLHd99isVBeXk56ejoxMTET7iY0mUxjZMTud9mEhARSU1NZOzuV9UWufpO2Pgv/Pu6KNvbUd/P3g2YACtNi1GijNHvqRF/BcrdyH20oG7hOnDiBEILk5GRSU1OnnSgSEhKm5bW9QAB9wCcAaaAxH+gd/bffX8y0EIW3SeFyKzEyMkJlZSVJSUmsXr3a6y/cYDDo3nq4RxRaR2x3kZaeSMUhfBOFjISkec3Jkyd1LxZ/Em73vgzpO1FVVYXD4VDvsBcvSWdLqSuEP2weGi3DdrPrwEkee78Fo0FhXrKB8/saOLMwhZKcBEy+mG8SCEX0ovUQhY8So52dnXR0dCCEwGQyBWz7P1nMlIhCUzU5CPwL+BpgUhTlP4BNwO9Hf+73bjulROFPG2E0GhkaGuKDDz4YM+PTGya69ZATwKQjtifPAn9bj48iirHPFULQ0NBAW1vbmEhIT0fqRKAoimrMUlBQgMPhoLu7W60imEwmUlNTyU5N5eo1uXzxjDwsdicHmnp5+3AH/6xt53dvNfDbtxpIiDKypiBZjTjyU6KDdr2hsMFzhzYxmpSURF9fH9HR0eMSo9K8OFTXMzg4qEanMwSRQC2uROZZuCon3xVCVGm2Jz4xZUThTxvhdDpV78d169bp6rwLJJkpicJsNlNXV+dzApierYcaUWges9vtVFdXExERMS4SCqTRazJNYUajcUwVwWKx0NXVRWNjo5qNT01NZVl2KmWz8tmYNczsokW8d7yH90a3Ka/XdgKQmxTl0m7MTeGMOckkx05cHDXVSkmn00lkZKTah6FNjEoT3oSEBDUi0eulqQeDg4OTUWQGDZok5dXAXCHE14G/aI/RQxIwRUThTxshOz5TU1OJj4/X3Z4bqCy7s7OTvr4+vxPAdG09Rn8sI/XBwUEqKiqYPXu2xy/JVHaEahEVFUV2drZqdS+TgQcPHsRqteJwOLAP9XJ+cQoXLs5wlaG7R9hzzEUar9SY2XngJAqwaFY8Z40mRZflJRIZQBl2KhOnMD6Z6S0xqvXS1FoFTkYxOgO9KJqAVYqirAdOAsOjf3qFEP4H0xBiovBnUQdjOz4TExPZt2+f7vPrjShGRkaoq6tDUZQxrefeoIcopOBKUaCtrU1NiHqTkrsTha9rCBWpuCcDBwYGOHToED09PercEOn2dMXKbK5clYPdKahq6VfzG3/ac4IH3z1BTISBFflJajVlfobvUH66icId2s9izpw5YxKj0hJPRhuBOFvBjHS3GgTKgF8Dx3ElN7OA+4C/TWt5VM5E2L9/P8uXL/fY8SmZXCYTpUGpXuhZ0NIROz8/n76+Pl1fVinh9gUpuLLbrDQ1NbF69Wqfd6Hpiih8wWg0EhUVpRrASnu8pqYmBgYGiImJITU1lQVpqSzLm8PXz5nDgMXOvoZetZpy72vHAEiPi2RtoSu/sbYwmcyEsVHhdBBFIFUmb4lR6Wwlk6JSMerrvcygZKZcHC8BrwIRQDwuUVcC0Op2nFeEhCi04/tGRkZ8dnzm5+erPw/0i+St5AkfOWJ3dXWxcuVK1dJMD/TkKCxW1+saFYUVK1b4vfapylFMBu57emlzL12uZGh+ztwUPlnsajE/2Teiir7ePdbNC1UuoV9RRqzam7JqdtKMiyj8wV0xKg1qmpub6evrU0nUfToZ6FNmTgU0icrVuERdbcAArrKoGejWe66gEoUnbYT7z5ubm2lsbNTV8ekP3iIK2VWamJjIqlWr1DJqsKzwenp6OFRbC0BUVERQSp4zDYry0WAdaf8vzWoaGhpQFEW9w25dmsn2slk4haC2bVCNNp78oIW/7m3GZFAomRXLohS4JKGPxdkJGENUhpUIto5Ca1Cjldy7J0ajo6P96ii8uFt9BvghLvXkGiHE+6OPFwA1uKoWMCrrHv3ZSuARIAZX1PCfbslJBdc241PAVbg0FAlAxuhztgPPTunWQ0/HZ3V1NQaDIeAZn97gKUfhzRE7WA5XJ06coKWlhcK5c+HDQ7rbzGdiRBHIa7iH5nK8oXtT25zUVBauzePaM/MZsTn48EQf79V38/aRDh6vHuHx6gMkRptcatHCFM6cm0xecvCFX6EUXLlL7qWWpbu7m9tvv50DBw7w4x//mM2bN3PRRReNS5x7cbeqAj4NPODhJY8KIZZ5ePx3wFeAf+MiiguBv8sfCiGcoyTwA+AHmutPAb4FHJLH+XvPQSMKSQyefjm9vb1UV1dTUFBATk5OsF5yzIIWQnD8+HHVgsxdyRmIitPT1sPhcFBdXY2iKKxevVrt5DTo1EcEShRThYm+lna8oXtT28jICAkJCaSlpbE6P4Wz5qZwzfJkjrWYaVdS2HOsm3fru9l9yGXak58SrWo31hQkkxg9+a/lVCoztVqWRx55hHPPPZfLL7+cf/3rX1x00UXjjvfiblUjz6XzNbOBRCHEe6P//z9gGxqiGD2vU1GUVFw9HrbRx7pHKyDSd9OvliKoWw/3O7EYtbA7ePCg2vGpB3r3szKisNlsVFZWEhsb61XJOZmIQvaC5OXlkZ/vav9X28wN+q430CjhVNumaO+w2illUlodHR1NjKJw4aJ0LlriIpf6zmG1N+WFqnae/NBl2lOSnTCq30imLDdxQqY90ynhdjqdbNy4kQsuuCBYpyxUFGU/Lhn294UQ/8LlcNWkOUbregWMIYDPALNwVT+GRztJDbhKpbq0FCGresiOTyEEy5cv1y1okYtfb/OUdE3yNyYwULm3/OykQEvbCwIfVT2MBn0E4E4Ug4ODdHZ2kpaWNk43MhO3HoHAfUqZ3W6noaGBrq6uMcOGMlNS+OyqHD63Ohebw0llS78abTz4biN/eKeRmAgDq+dItWgyc9P1KSqnkyiCHBG2ArOFEJ2jOYlnFEVZEuA5InDlKxJxNYd1AlcKIVr1niAkRCFnfM6fP5/m5uaAvpB6twhCCNra2uju7ubMM8/06ygUaOjvcDg4fPgwvb29HgVaMqIw6jyv9vWldiQ9PZ2DBw9it9vVxKA0dJkqTMU2x2QykZiYiMFgoLCwcJynpmxqK8lKZUV+Et88t4C+ETv7NE1tbx1xbfUyEyJZW5jCWYUpbFqc4bU3ZbqIQggR1AqPcBnqWkb//YGiKEeBYlxGNVpXpXHuVppIoUsI8WvtzxRFuVBRlK7R8/tF0Ini6NGjdHZ2qnmCkydP6r6Tgz4RlXTElpl3PbZjgfziHA4HHR0d5ObmehVoSWWm0ahvSyPzHocPH6avr49Vq1apj8v+jI6ODo4cOYLRaFTLkzExp47tnS9oF4+7p6anprbU1FQ+MT+FDQtd5cmmnmHeGy3D/vNwJ+8d6+aiJePHN0p8XOapKoqSgWuhOxRFmYur6/OYEKJLUZQ+RVHW4kpmXg38ysNzo4G7FEUpB6RJTS9wPy4n76knirq6OoQQakkSXAvfn3hJC39bhP7+fiorKykoKCA1NZXqao+GPhOGTLzGxMRQXFzs9TjZ66F36yGE4NChQ6SkpLBixQo1fwPj+zPMZjONjY0cPXqU4eFhEhMTSUtLIyUl5ZQdfOtt4bo3tXkSOkm16KXLZnHZ8mwcTkFbv8UngU4XUdjtdr+/I627VV5eHs3NzV8GunAt9AzgRUVRDgghNgHrgR8rimLD1eV5vRCia/RU3+Cj8ujfcUtkAsuAC4AUXK7e8hwmXNuPQb3vK6jfuuLi4nF310Aat/wd39TUNMYRW876CBbk+RcvXkxjY6PPY2X3qJ6tx8DAAGazmYKCAubNmwf4zg9ER0cTGxvLokWL1CldnZ2dNDQ0jJFZ+1MI+sNUiqD0vpYnoZOnpra0VN9zXqeLKCbgbgUfOVs97f4DIcROYKen84xqLUp8vFQ98C4wB3gbV3QRhSuK+JkYtdPTg5BWPSA4RCFdtwFWr16tMnYglQxfcDgc1NTU4HQ6WbNmjeqV4fM5o0Rh8rP1kH0gaWlp6sTsQOBuRquVWff396sLJzU1dUbPupwoKXlqauvq6lJzO9LBOyUlZUw/xnQRxcDAwIyQbwMIIY4AR4CnR8up8u40HOi5Qh7HmkymgLYe7kQh5d75+fnjHLGDQRTDw8OUl5eTnZ3N7Nmz1ZyBXqLwtvUQQnDkyBE1GXrkyJFxx3lbPL4Sr+4ya+3CcTgcJCcnk5aW5nPS9nQgGNGLeyOXw+FQ1aLuTW3TGVHMBPm2hKIoUbjmmW7B1RjWhavycRy4ZMb4UUwmomhpaeH48eOqIa07JvvF6+jooLa2dpw3hR7jGpmjMHmo0tjtdioqKoiLi1OToZ4Wf7AXjtzfa701ZSTjKeE7lVoNEQI/CqPRqBIDjI22hoaGqKqqUhOjnqwUQ4GZ0hCmwSzgO8BPgQXANcDNjGowZowfha/GLW/HS7m3zWYLmtxbQn4ux44do6uri1WrVnnUMej1zHSPKKQl/9y5c9XBtfKcoW4zd9/fy6YuqZaUTV3apOhMy1FMBtpoa2BggMLCQo9NbaFMCs+UhjANknFFEf/AJdZqUBTll8CfgJ/p6fOAIBOFpy+CyWRiZGTEw9GeIcU5BQUF6lYgWDAYDFgsFg4ePKje7Seq4nSMrmMtUbS3t3PkyBGPEdB0NIW5qyW1TV0Gg4G4uDhsNtuULOKp7h7V09SmTQoHK9qZgRGFFSjHVRk5pijKZ3EpOHUb68IURRR6cxRtbW00NTWRlZXFnDlzgn4tTqeTDz74gHnz5o2527tDn/Lvo2Smw+HgyJEj9PT0eB1xON1NYe5NXVarlebmZjo7O9m7dy/x8fHqNkXPiMZAMd26Bm9NbVo/TUkck9GuzLQcBa4y6N9H//4V8DPgKHAjgN7Kx4zIUWgdsYuKirBYdGlAAkJLSwtDQ0OsWLFiQtUHd8gchYLg8OHDJCcns2LFCq+LYaY1hUVGRpKWlobFYmHBggUMDAzQ2dlJVVUVTqdT3dsHKyk61RGFP/hraktMTFS3KYHY4s0UdytNknIxsHK0P+Rd4JyJnG9Kth6+iMLdEbujo4OhoSGvx3t6TV93K6fTqTpppaSkBC2pJSOKDnM7hbPzfYqz5HXO1KYwbVJUip60A4fk3TYtLW3Cn99UazYCgaemNjlIWTa16SXOgYEB0tLSJvsWgokEYKHi+vBNuARXQk9eQotpjSg8OWIHWiWR+QRPv7yRkRHKy8vJyspi0aJFlJeXB0V3AdA3MABARnq6VzdvLWZaROEL2oFD2rutNimYlpZGcnKy7qTgVBLFZLc5ngYpu09qk9sUd3erwcHBkGybJwBpWjMILAEeBN5gVEOhKEqVEKLW+9PHYlpyFE6nkyNHjnh0xJ4oUbhDemUuWrRI3WoEQ3chhGuCWVd3DwBxsTET6h7V8zqhht72ePe7rZw+fvz48TElyvj4eK/nO5WIwh2eJrVJd6vh4WFV9JacnOw3RxGIu9Xoz24Hvgw4gP8QQkgPiQuBX+LymXhQCHG320vJL5AD2Af04xpSHA1k46p61M6Yqof7wpdTwFJSUjw2XE1Wyak1sFm5cuWY9vZAzGs8wW63U1lZSUxMDFmzsqHuOEYdmguY/mRmsKAVNcF4ibXsBHVPip7KROEOb01td9xxBy+99BLNzc1ERUWxfv36caX3QNytFEVZDFyJKyLIAV5TFEXucX8DbMSlh9inKMpzQoiD8rlCCDFKAv8E/qkoyhwhRIP7e9G7BZnSrYe8y8sxe54Q6F1fe7xcyNHR0R4NbALxpHDH4OAg5eXlqkvX35uOYDR4FlJ5w0xd/JOBu8Ra2wkqk6JpaWk4HI4pq3pMl7vVz372M3p6eli9ejUvvPAChYWFqsO5RIDuVluBJ0ZbwesVRTmCq+MT4IgQ4tjo854YPVYlitFkpnOUbLYD1ymK8v+EEI8rinIjLu9N3bMxgk4UnoRFYtQRu7Ozc9xdftwF+Ul+ukMSkRQ6+bLbm+jWQ+ZSli5dqhoCO4TAoOg/p7vac2hoiBMnTqhzIyY6VWwyCPZd3lMnqPTVbG9vZ2hoiKysrJArJaezFDsyMsKFF17IokWLgnG6XOA9zf+1LlYn3B4/w+25BlzbjhuB/bjyEzK8WY3LjXvftGw9PMFqtTI8PIzD4RjTfu4Ngd71DQYDZrOZtrY2r1Jv7bGBEIXT6eT48eN0dnaOy6U4nEw4opCRVW5urpogi4mJIS0tjbS0tGlPZgYL2r290+kkPT0di8WiJkXl5HH3hq7JYjqJQjpyzyBk4BpQPA+XlgJcjle6rfohxEQhHbEjIiL8lg8lAslRSLVdf3+/3wE8EPgIwvLycqKjoz0qOJ1CYBjt4dAbUTidThobG2ltbVXzM3Ii+dDQEJ2dnRw8eBCbzaYKgmZac9dkEBsbS2ZmJvn5+WMauurr6zGZTKrgKy4ublJkOZ1EEWTBVTOQr/m/1sXK2+MS8kv5Lq7BxJ8A6hRFWY1Lpdky+vPp6fWQC0friL1//37dYa7eL4jFYqG8vByTyURBQYEuUYxeohgaGmJgYGCMma47HE6B0aDoaiADVzTR0tJCbGysGllZrVZgrNx49uzZalnXU7QRzFbyqW4K0/5u3Ru6LBYLnZ2dY+zxJHEEOgd0uokiiBLu54DHFEX5Ba5k5nxgL67S53xFUQpxEcSVwGe1T9Q0e/0vcCuurcY24Abg/wkh9rod5xNBJwqbzUZ5eTkxMTFqQlEu0GCFlzJSWbhwIT09PZOe16GFzEfEx8f7NOt1CoFRZ0RhtVppaGggJiaGpUuX+t2umEwmIiIiWLBgwbhoQ/prBquVfKY0hUVFRZGTk6NGWP39/XR2dtLU5DKalu9ZT1/GdDtw+9KWuLtb/ehHP+K6667bjgd3KyFEtaIoT+JKUtqBb0rJtaIoN+Cy2zcCDwshPFq9CSHswE8VRXkZ1zjBf+v1ydQi6ERx+PBhcnJyxvRSyO3EZIlCCEFDQwNtbW2qJ2dfX19QiEIIQX19PR0dHaxatUrN2nuDwykwGFzn9NUdK6370tPTiYqKCnhhukcb7q3kcrRdWlqabqfz6YAcDKUH2qRoYWHhuL4Mf+95Oo11/b1HD+5WfPnLX34aD+5Wo+e8E7jTw+Mv4Rr64xOKopThai3PwZXgvEJRlF8JIQ75e64WQSeKJUuWjMsxBFrJ8ATZem4ymcaUPgPJaXhb1NKsNzIyUt0W+NtSOAVqROHtONlNWlpaysDAgG5puq9zalvJtdFGTU1NwNHGqeJH4d6XIdvnDx06hM1mU5OiycnJqlZmuiKKmdLToun1+DXwGi5lpgW4DHhIUZSLhBC9vs6hxZQ4tQZqsAtjP/DBwUEqKiqYPXs2ubljZpwEPK/DPUqQw33cz+1vS2F3OjF4qXrI6KSzs1PtJh0YGAi6hDsY0cZM2XrohXv7uMPhUJWix44dIyIiApPJRHR09IxZtNMBTe7BJoT4keZHdyuKchk63bclpowoJiLLNhqNquekt6HG/kJ/9+vQLn7pcCWH+3QNWkmNixxzDd7gdLoiCvfjHA4HVVVVREREjKmW6E16SkzWuEYbbcg7bzBzG4EiVIvWaDSqiV5w6RiOHj1Kd3c3+/btU5OigXaBTgRWqzUkLfqThF1RlHuBF3F5UJQBx4AYRVEcQghdiyckVQ93TESWbbPZOHLkCAMDAz5Ln4GcW0YJUuZtNptVh6unPmzm+8/V8M+bziYrMdrvwnYIMU5HMTIywoEDB8jNzfVYLZnKpjD3aEPODtFGG9HR0UF1MfeFqbq7R0dHqy3iWVlZalL0xAmXPklWWhITE4N+PTPNtEZRFAPQgWt48Tm4+jxygcPAW7gEWcv0nGtKIopAcxRSw5CWlsaKFSt8/kID3XrYbDYqKiqIiIgYIwBbnp+MEPBajZnPnZGvI6IYSxS9vb1UVVWNaUJzf0/T2RSmnR0io40TJ07Q1dXFvn37Qh5tTEevh3a0Ibgqcl1dXbS0tHDo0CHi4uJU4ghGIngmOXCD2sfxWXDlLHBVSIy4CMNEAOt/xuUoenp66Onpobi4mNmzZ+s6t96qh81mo62tjeLi4nGO3vMy4pibHsvuQy6i8JejkBJuRVHU/obly5d7nVo2k9rMZbSRkZFBREQEBQUF4/wnZDgfrEpKIFWPYLyWJ7KLiIggKyuLrKwshBAMDg76TIoGipnmbjVKDiW4+kOcuPISVlyTx/4RyLlmzNZDCMGJEydobm4mPT19zEBgX9Arours7KSuro6kpKRxJCGxcVEmD77TQM+QzX/Vw+nSUTQ1NTE4OMi6det81s89EYX8TOTdb7rgKdrQLqBgRRvTTRTu1xIfH098fLy6NXNPimodzPVc+9DQ0IyKKIBMXK3oI7jGCUbjMrJpA/6hqYz4xZRtPaQK0RMcDgfV1dUoisKaNWuoq6vTHSX4IyGt9qKkpETdq3rC+YsyeOBfx3mjzkxJnG8CsjucWCwjGAwur0l/Bi5aohBCqN2UTqcTh8MxhjSm0wnKU1UhWNHGTCIKd3hKinZ2dqqeE3qs8Wba1gNXPiJSCHGepx/qJQmYAVUPWZ7My8sjLy8PZbSSoHer4iuikARkNBpZvXo1IyMjPhf/0pxEZiVGsbvGzNLVcV6PHRoaorO7m6jISGbPnu2TfCQkUcjhQoqiqF84mWCVP9P+LT+PUMHf4tUTbchwfab0pARDBRwdHa16Trhb4wFq2Vk71nEGWvU3AS8rirId18CfYVzRRY8QoieQE01rjkLKpWV5Unv8ZCMKOQFMW4Hwt01RFIXzF2aw48MWbljhmShkuTEuLgGHVQTUaKYlCe0C1YrHRkZGqKqqYvbs2QghVMKQP9ceP9XwFm2YzWaOHDkSktzGRCAJNljwNNaxu7tbHesom9haW1sDdrfq6uoiLS1tN1CAazFfLoToVhTlE8CzuOaHAuwSQvwYdLlbSThwDSi+H1ePiAnXnI93gO/rbTGHacpRCLdxe+6150DVlu4LVS5mTxPA/C3qjYsz+eveJj5sHeG8hLGJSW3n58NHD2I02HUlKYUQREZG0tPTQ1VVFRkZGaSlpY173319fVRXV7Nw4UL1up1Op/pHfiZy2xLqaMMf3KON4eHhcboNGW1MJUKtzIyMjByXFH3jjTd4+OGH6evro6enh+uuu25cx7Qnd6u7774b4HUhxN2KotwG3AZ8d/TH/xJCbNaeQ1EUI37crTQEsAk4G1iFq8/DhMuToh/0u1vBNJRHrVYrFRUVJCYmerTCg8CIQnusNh/hySBHD1Gsmp1McmwE7zYO8omiZMD1xaupqVE9NYxGIw6BrjZzuaWIjIxk7dq1DA4O0tHRQXl5OQBpaWlkZGQwNDREfX09ZWVlYyon2kSn3KK4k4ckjEAXRzBLsIrGW9NTtDEyMkJzc/OURBuTkYsHCpkUveSSS2hoaCA5OZnZs2d7/E54crd69tlnAf48+t8/A2/yEVF4whr8uFvxUet4N/BPIYQZMAf41sZgSrceUmswf/58MjMzfR6vV22pTQhWV1djMBg82uBpj/UFk9HAecXpvHKwDYvNgdVqpby8nPT0dAoKClRi09NmLklCaggUjS1+YWEhVquVjo4OqqqqGB4eJjMzk4GBASIjIz0mR7VbFGBcPsM9IapnsYQqwegebezdu1cdnSBLkdLJO9iLOpidyoFgcHCQefPmcfHFF+t+TltbG0KI1tH/ngS0LctnKopSjss74ubRDtFc/LtbSQduO3C2oigPAHtw5ScU4IC039OLKdt6DA4OcvDgQZYtW+Y3Myz36npfz+l0sm/fPnJycnxqL/RqGc5flMmuA63sa+jhZEsTxcXFqgOzhNa4xtM5tVsFbwvBZDLR09NDUlISq1evpr+/n46ODurr64mIiFAXmjdthqdoQzuJfTLRRjAhryE/P39KchtT6c+pxWSH/wghhKIo8sv0ITBHCDGgKMpFwDO4/Ch0nWr0bwsu45oBXFuQGFxDix8BaqY1RwFjF6TD4aCurg6r1cr69et1zYEIZOvR1dXF0NAQq1ev9jtfQ+/dc928VKJNCq/XdnDf587w+Mt3OAUREYZxUYp2a+CetNRCKkTT09PVGavakXfDw8NqL4rFYiElJYX00Rki3qIlGB9tSPKQ/zYajerCncruUS385TYmG21MV/foRARXWVlZKIqSLYRoVRQlG5feASFEnzxGCPGSoii/VRQlHd+uV/J4oSiKUQjxlqIobwNzhRBH3F97xuQoZOUhOzubvr4+3cNi9MiyhRA0NjZy8uRJYmNjdQ3h0QMhBM2Nx1mSqlDVrRAb6zn6cRnXjNdH6CGJoaEhdeK5ty1YTEyMxztwXV0dsbGx6kLz5njlL9qw2+1T1ufhC95yGx0dHWOijUAMeaera3QivR5btmzh3nvv/SJwN/BFXJUOFEWZBbSNLvo1uLwkOoEe/LhbjZKEQ1GUUlz5i2sVRbldCPGEoijfAt4Vow5XehEyopB3Q1l5aG52t/TzcVEmk2/5tMPBwYMHURSF1atX895773k9NhBoOz8/VZrDT15rory5l+X5yeOPdTKmzVwvSXR3d3Po0CFKSkp0m7C634FlQrSqqgqHw0FaWhrp6eleG508RRtWq5X29nYyMjLUfNBMq6QAagdsbW1tQNHGdBGFr4jCk7vVbbfdxr333rtRUZQvAw3A5aOHXwZ8XVEUOy79w5WjAim74sPdalRtKe8AN+Jy8X4TkCW2VbiSnHunfetx9OhR1SlqIh6PviIK2aGZk5NDfn5+0L4Q7p2f0U0nMSpN7K4xeyQKaYWnTaZqk5ae0NLSQlNTE8uXL5/wXlwrPS4oKFAbnU6cOEF/fz8JCQlkZGT49JqURj2zZs1SrefccxvBJI3JbHGCHW2EEtLr0xs8uVsBCCE2eHjs17hMZzwd79XdajQCuRyX32YirhzFIlxdpOAijIAcuCFERJGUlERBQcGEv2TechTSK9NTh+Zkws2enh6qq6vHnDcpJoLSrCh217Rzy8aiced2WeG5HrNYLBw7dozMzEyPdxQxOtdkcHCQlStXBjUj797o1NfXR0dHBw0NDRgMBvXuLEVBcjs4b968MQlarchNkoa72Guiv89gbQX0RhvTlXuZQU1hNwK1wE5c5dSzgCpFUc7AVfVoGj1ueiXcco6DO/R+YTwJtE6cOEFLS4tPfcREFmBLSwsNDQ3jOj8NBgNn5Ebz+w96qWsfZEHW2C+AwylcE1YcDtasWaPazg8ODpKcnEx6erpKOtXV1cTExFBaWhrSkFhRFLWtet68eVgsFjo6Ojh69KjasNTb20tJSYnHnI7WZAfGi70mUn6F0OUM3KONnp4eOjo61BzQVEcbQ0NDXitUU4xBYDEuBeaPcKkzv4rLN/MWIcSHMAN7PSCwxawlCqfTycGDBxFCsHr1ao/PD8S8V5ZTFUWhrq5OrZi4J1oNBgOrc6J44EPYXdM+hiiEEC6iGE1masfqySG+MjweGRkhIyMjqNskvYiKilL7Fdrb26mrqyM1NZXa2lqioqJIT08nIyPD6zbIm9hLqxDVU36diuSitqmrp6eHoqKiCeU2JgNZVZoBSAcuweVmVYtLjGXCRSArFEV5QwjR4eP54xCy8qg7AlnM8lg532LWrFlqCdETAum3kPM0qqurSUxMZNmyZR7PqygKSVEKy/KSeK2mnRs+MRf4qLIhtx7uzzUYDGp+oKuriwULFmCz2aiurtaVeAwFWlpaaG5uZs2aNapsfGhoiI6ODnXgUGpqqtre76/8GhER4bH8Ko9zjzamsgohb5Leog2Z25BNXcGKNqar0uIFDuBlXHmJJFx5iUg+mhCmT6SkwZRFFCaTCbvdrstTUFEU7HY7H3zwgVfHKC0C0V0IIfjggw+YO3cu2dnZXo+T5LNxUSY/f/UwJ7qHyUuOVheHc9QKzxPMZjNHjx5l6dKlarlszpw54xKPiYmJauJRb+k4EIhRk9++vj5WrFgxhqRjY2OZPXu2asrb1dVFa2srhw4dIj4+nvT0dI/9KNrPR6/Ya6qJwp3o3FvIZW5D6nuCFW3MILJIwTXc2Lu3Q4CYMqIIZDGfOHECi8XC2WefrYvxAzGv6e/vZ9myZeOUlp7OKYRg46IMfv7qYV6raefzq3PUnzmEwORGFDKXYjabWbFixbhF5p547O3txWw2j1FiZmRkBOUuJ+XSiqJQWlrqcwGYTKYxdvgDAwOYzeZx/Sjx8fG6y6/aaMNisaj/DnX5VY/YKhTRxnQlUL3gIVxJS+mbqV5cIHkJLaZ86+ELMh/hdDqJjY3V/Yvy15auFWelpKToEsXIXEZ+SgzFmXHsPtjOF9bkqu/N6XQ1hWmvvba2FqfTyfLly3U5LGlbl6USs6amBqvVqi7OpKSkgO9SDoeDiooKkpOTx/Sn6IG2H2Xu3LlYrVY6OzvVRG1SUpKaqPUWBWmjjZGREWpra8nPzx8j9pLHhKLPI5BzBivasFqtM2YAkxDiZ5p/65/K7QNTuvXwRRQyH5GVlcWcOXPYs2eP7nP70l1oOz9Xr15NZWVlQEOFnU4n5y1I54F/NdA5aCM93hUluFy4XcfabDYqKytJSUkJeGFKaJWYcivQ0tJCTU0NCQkJ6lbAn+W8bGLLzc0lJycn4OtwR2Rk5JhEbW9vr+5+FFmKnT9/vroQg9HI5guTlW/7ijaioqJUUnG/iQ0MDMyUikdIMCO2HlLHsHDhQvULBfr3fN6Iwmq1cuDAATIyMtQFrHeboigKFosFq9XKxoUZ/P5fDbxR18FnVrgWn9PpagobHh6moqKCgoICn7NKA4H7VqC/vx+z2UxjYyMGg4GMjAxVG6GFLAtqF2YwYTAY/PajZGRkkJycrH4uixcvHmNK5J7b0P6Rv+/JkEYw+zw8RRtdXV0eow1/7lbeTGuuuOIKXnvttcOMNa1RcBnTXAQMAdfIkqaiKF8Evj962p8IIf7s/lqhwJRuPTy5XElDXTlLVCLQcqr74u/v76eiomJc56ceopBf2MzMTA4cOIDBYGBWvIlXqk+qROEQArvNRUTuiyGYUJSP5nBKbYTZbObw4cOMjIyoi1NRlICl4ZOFp36U9vZ2ampqsFgsFBQU+AzH/Xlt2O0uY6BAxF6hbAiT0UZeXt6YaOOdd97hl7/8JSaTifr6egoLC8c915tpzYYNG9i9e/d8N9OaT+HqFJ2Pq4X8d8AZiqKkAnfgkmEL4INR05qAlZaBYtoiCvctgTshBFpO1S5+OV2srKxsHMv7IwohhEpoc+fOZe7cuYyMjHBWQw3P1fTwxtt7yMtKx2530NvTw/ILV0/p3jQqKkr1F5WLs6Ghga6uLlJTU+nv7ycqKmrKJ1ZJ1aTBYKC7u5uysjJ1QLPT6ZxQP4p2e6J3izJVnaPaaKO4uBi73c4f/vAHbrjhBs4991xuvfXWMcd7M61588035X+1pjVbgf8bTTy+pyhK8mhn6SeA3UKILgBFUXbjGu7jWRseRExpjkI6cVssFg4cOKDmI7xFIIFoI2SW/dixY3R3d3udLubPjNdTU1d0dDTbVxeyq3o/AwlzGB7uwO4UKLik2aEscfqC9O1wOBycc845qhJTVitk/sBbtSLYaG9v5/jx46xYsULdz8t+FDmtS28/iow2ZIOgXrHXdLWYJycns2rVKu6//37dz2lra9OW6LWmNZ7MaXJ9PB5yTPnWw1s+wh0GQ2BO3Ha7nfLycqKiolixYoXXL4snotDT+bksL4m0uAieeb+eW89OR1EM5OXmkJubPqbEmZGR4VPtGCxo+0eWL1+O0WgkMjJynHuWrFZo/SxCoR5saWmhpaWF5cuXj1v8ERERzJo1i1mzZunqR3GHHrGXVEXKispUY7IO3G6mNTMOU7r16OnpwWw2+5yopT0+EHfr48ePU1RU5HW4j4Q3oxltx6TH17DbWJoK/261UTCvCIdoxeihxGk2m8eoMDMyMsZYugcDsowcERHhtX8kMjKSnJwccnJyVFm51lFKajYm0t3rDrn1kYTlC/76UfQQmi+x1/DwMDC2BDsVmIi7VVZWFq2trWRnZ6M1rcG7OU0zru2H9vE3J3rNgWBKiMLpdHLixAmGhoZYt25dwP0evtDT08Px48dJT0/3SxLyvJIotJ6WvkhiYGCAqqoqLi7L5c0Tx3nvWLdLR+EmuIqJiVHVjjLcbmhoYGBggKSkJDXcnswd3W63qw1Pc+bM0fUcKSuXCldPfhYTITQZ1QwPD1NWVjahRantR3E6nWPayPX2o4Dr99rR0aEOegJCUn71holEFFu2bOHPf/4zt912G2hMa3C1iN+guIxzzwB6Rx2wXgF+qiiK7Oi7ALg9KG/AD0JOFBaLhfLychITEzEYDLoXiR6iaG5uprGxkaKiIoaGhnSdV0YUnoxvPaGzs5PDhw9TUlJCZHQs8btP8Fpth0tH4WNNacNtbaPY0aNHiY6OVkucgdzR5Wc5e/ZsZs2apft57pCzOaSsvLOzk8bGRvr7+1VBVVpams/flRCCQ4cOAVBSUhKUiMlgMIwrR+rtR5H5keXLl6uJXG/OXlo7wGBhcHCQ3Fzv6QJvpjWXX345t99++2HGmta8hKs0egRXefRLAEKILkVR/hvYN3rcj2ViM9QIaY5Cum4vWLCA2NhYamtrdZ/DF1EIIcZ0fvb29jIwMKDrvDKfIUtvvr4sJ06c4OTJk2Pk2OfOT+MftR2j59K3ODzd0c1mM5WVlQghdCUdBwcHqayspLi42G/vSyBwzx+4y8oloWlL106nU22dnzdvXsgSpXr7Ubq7u2lsbByXH/EmLffktaE9fiLwN3fUm2nN66+/Dm6muaPVjm96Ol4I8TDw8ESvc6IIWUQh7/YyH2G1WgPyaPRGFNKUVtv5qXebIhNeTU1NREZGkp6e7vGuKYnIarWOa6basDCDF6tcW0njBBeIvKMXFBR4TDpmZGSMMdHt6elRNRKhNEbxJCs3m83U1NRgs9lUf4fjx48HtPUJBrz1o+zbtw+r1Up+fj4jIyOYTCafXcbevDYmK/by5251qiMkRCGdlbU+D4E0hXk7fnBwkPLycubOnTsm9NYropL1/IiICHVxRkVFkZmZqW4D7HY7lZWVJCYmUlxcPO5Ld05RKpFGA1aHU3dE4QvuSUetiW5cXByRkZF0d3dPyj5votDmXOx2O+3t7VRWVqIoCpGRkbS1tfkscYYKsh+lv7+f6Oholi9fTm9v74T7UYIxWGkixrqnEkK29Zg/f/6Yjjo9ztpauJdHpZ370qVLSUxMHHOsnonm2qSllCHPnz9/zDbA4XCod6eCggKP54qLNHHW3BTePNzpM0cxEWj36GJ07GJbWxsRERFUVVWpib3p+EI6HA6amppYtGgRGRkZY0qcUmyVkZExZf0OTU1NtLe3s2zZMoxGI7GxsRPqR4HgDFaaQTZ4IUHIth7SnVr7/0AgBVrazk9vZr2+Igp/SUu5DUhJSaG6uprs7Gx6enp47733SE1NVXsXtM/bsDCDNw93BiWi8HbNhw8fxmKxcNZZZ2EwGNQyopRuy2vzZjQTTMjmLm1+RFviHBkZoaOjg7q6OvXa0tPTQ+Ym1djYSGdnJ2VlZeO2joH0o/jS2uj12pDHhbce0wQZUVRXV6s2eN5+sd4iCq3S0tcXtq2tTc2Yy6Sdw+EYkzyTJjNpaWmctyCNn79qIith8hoEd8hEYXR09JhqgraM6H5tgXSXBoqBgQEqKytZsmTJuEhOIjo6eoysvKuri7a2Nmpra4mLi/M6lHkiOH78OL29vbrLsd76UWpra3XPRwHvXhsy6u3t7Q1HFNMBp9NJS0sLBQUFXmXeEt5EVP5mbAghOH78ON3d3axYsWLMIjMajarKUlsNOHbsGFFRUTxxZQHZWcFtBJOJ2szMTPLz870e535tchvQ2Ng45meTNcDp7e3l4MGDlJaW6t7uuF/bwMDAOFm53D4FGmUeO3aMgYEBli5dOqFIxdt8FL39KOA52qipqaG6unqmuFuFBIofw5sJS0ptNtu47cC7777LmWee6fcD7e/vZ//+/cTGxrJq1Sq/r+V0Ovn3v//NmWeeqZskZFOawWBgwYIFAX3xZF6jo6NDLW9mZmYSGxs74S+L9OMoLCz0OcBZz3nMZjNms1mtVEjtQSDXJvUjZWVlQfOVlBUes9msqjDdKzyeIIVdIyMjLF68OCTbGakn6ejo0N2PAlBXV8fVV1/No48+ytKlSyf68jOeYUJGFJ5G1r333ntenbQlZOfnvHnzaGtro7S0VNfrSRLSI6KyWq1UVlYGxR1b++UfHh72mtfwBRneL1q0SC1NBgNSe2A2m+nr6xuzffL3O2hoaGDZsmUh60LVqjC7u7uJiYnxuA2QSV2r1crixYun5K4tozSz2UxXV5fXfpT6+no++9nP8sgjj7B8+fLJvOSMJ4op3XpIlytv2gVt56fNZqOlpSWg8+shCSlcch+AM1Foy5u+8hreFqY0QgkkvNcLd+2B+/ZJiqm0Zdempiba2tpYsWJFSLth3Ss8Q0NDavXJ6XSqOZfW1laEEFNGEjC2HwUY149y/PhxhoaG+MMf/sDDDz88WZI4JTClEcX+/ftVlaYWDoeDyspKoqKi1G2AxWKhqqqKlStX+n0tIQR79uwhISGBzMxMUlNTPYanclEuWbIk5Blq7cLs7OxUF6a2EevkyZM0NjZSVlYWlOasQCDl0WazWe33sNlsDA8PU1paOq3zKWw2m7ow7XY7mZmZQemTCQacTifPP/88d999NyMjIxQXF/Ob3/yG2bNnT+a0Mz6iCBlRaDPCEhUVFRQWFo5ZpLL0JrPmEna7nQ8//JA1a9Z4vzhNPkK7MLu6uoiLi1OFVCaTiebmZlpaWigtLZ3yRQnj8xqyJdpTW/ZUQ845GRgYwGQykZycrOYOpmNhCiE4ePAgkZGRzJ07V03WSsL11ygWSpw8eZLLLruM+++/n/Xr13PkyBFyc3Mnm8cJE4UWBw8eJCcnR92Hy1micuK5FtoEpccL85G0lNn29vZ2Ojo6sNlsmEwmSktLp90AVS6CoaEhIiIixmgiAslrBPN6ampqMBqNFBcXI4RQW9K1uYOMjIwpcc2SLfQxMTHMnTt33OchIyH5e51osnYiaG9v59JLL+XnP/85GzaMmys8GYSJQova2lp17yl7QZYtW+aVjd99913OOuus8Rels7LhcDioqqoiMjKS2NhYzGYzQgh1CzDVCkd5PfHx8eoikHkNs9lMb2+v7oRjMOB0OtXrKSwsHPc5anMHZrMZmFx5U+/1SPMdf7Db7WqlQiZr5fcr2PmVjo4OPv3pT/OTn/yECy+8MKjn5nQmCqfTic1mG/PYkSNHiI+Pp6enh+HhYZYuXerzF+qJKPS2h4+MjFBRUUFeXt4Y23pZpWhvb2dkZIS0tDQyMzNDPuLPZrOp4xG9+WboyWsEC9LXQlZ+9MC9vCkVmP7Km3rgdDqprKwkKSnJq3zeF7SVis7OTkwmk5qsnWwU2d3dzac//Wn+67/+i82bN0/qXF4QJgotZO9CVlaWrvZkd6LQGt/6+mL29/dTVVXFwoULPU7tlnA4HHR2dtLe3k5/fz/JyclkZmYG5YuvhczDBFpp8aTXCMbd3GazceDAAfLy8nyOVfQFbQNbd3e3qsBMT08POOciBxalpaVNNimoQsrKzWYzFotlwpL33t5eLr30Um655Ra2b98elGvzgDBRSAwODrJv3z4yMjJYsmSJrnNIotC71YCP5n4Gmo+Q5jLt7e10d3cTHx9PZmbmpMPYvr4+qqurJ23pHwy9Bnwk7Jo7d25QysPwUU5Ikpqv2SPucDgclJeXk5mZqcuhbCKQ27uOjg56enpULwt/pNbf389ll13GDTfcwBVXXBGSaxvF6UsUQgjVdVs25cyaNQuDwaBr/wkuoli7dq1KFP7k2I2NjXR0dFBaWjqpSoJ26E5HRweRkZET2gJIdWOwk6gTzWvIAUELFizwGWlNFnL2iPvd3J3UpCFydnZ2UKaa6YGW1Do7O1EUxaOYanBwkMsvv5wvf/nLfP7znw/1ZZ3eRGGxWGhoaKCtrY1ly5bR09NDf38/RUVFus6xZ88eVcLtT44tbdkWLlwYdImvNqGnNxna0tJCc3MzZWVlIa0W6M1ryO3YVA4IgvGkJqXRSUlJVFVVkZubO+HtTzBgsVjo7OxUI7X+/n56enr4v//7Pz73uc9x7bXXTsVlnL5E4XA4OHDgAEIIlixZgsFgUDPUCxYs8Pt8IQQffvghkZGRZGVlec0byLmfqampfpvHggGr1YrZbKa9vR2LxaLmDWQyVAhBfX09fX19LF26dMp1CENDQ2pZWOY1oqOjOX78eEjUn4FAJhzb2tpoamoiJiaG3NzcoE1wnyycTid79uzh+9//Ps3Nzaxdu5b/9//+HytWrAj1S5++RGG322lsbCQnJ0ddvD09PbS0tLB48WKfz5VJS2lC0t7eTk9PD4mJiary0mg0MjQ0RGVl5aQbqSYKWZ4zm81qMnR4eJioqCgWLVo0LfMltLBardTX19PS0jJGqDQdeg3tNR04cICCggISEhLUvIvURGhJd6phsVj4whe+wAUXXMANN9ygyvBDlTvR4PQlCnB98Fr09/dTX1/vs9HL27QuGWK3t7fT2dlJREQEQ0NDLF26NKT7bb2w2Wzs378fcL2HYCVDJ4OTJ09y4sQJ1eDFfQsgr2+qoh6r1cr+/fuZN28e6enpY37mTrp63cCDBZvNxjXXXMO6dev4zne+M9VEdXoThXSokhgaGqK2ttZjE00glY2WlhYaGhpITU2lp6eHyMhItR9gqmdugut9lpeXk5ubS05OjpoMlaQ2Hdd34sQJ2tvbKSsrG0dUcgsgry+Ueg0JOUayqKjI76R1b3kX9wa2YMFut/PlL3+Z5cuXc/vtt09HNBMmCu35rVYrFRUV4zwmAjGaOXr0KAMDA5SUlKgLQO7LzWazOoU8MzNzSnoBZCVh/vz5XheA9voAMjIyVP+KYEPmSPr7+3UbvLgna4OtvhwZGeHAgQMTrrbI6+vo6MDhcKjXF4yZqg6Hg+uvv5758+dzxx13TNeWLEwU2vM7HA727dvH2rVrP3oBnUpLh8OhWsTNnz/fpyJTJhsdDoe6KEORxJMOUIFUEmTLcnt7O1arVVWGBmPsoBwzYLfbJ5wjCZZeQ0KKzRYuXBgUrw3ZWWo2mxkcHFQb2Lx1DPuCw+HgxhtvJDs7m5/+9KeT+vyvvfZaXnjhBTIzM6mqqgLglltu4fnnnycyMpJ58+bxpz/9Sf0MKioq+NrXvkZfXx8HDx6sAlYLIUYmfAEhRkiJwt3lSraDS7WlXqWlxWKhoqKC7OzsgBJLNptNJY2RkZFxFYrJQPo6lJaWTjhjL/fl7e3tDAwM6HZ88gTp2BUREeGTSAOBp9JmIHmNoaEhysvLJy028wbtTFXZwCa3KP62eE6nk29/+9skJiZyzz33TDrx/NZbbxEfH8/VV1+tEsWrr77Keeedh8lk4rvf/S4AP/vZz7Db7axYsYK//OUvlJWVoShKOtAjhNBvUz/FmNIsm/bL6y1p6Q4599NXaO8NERERY0xlpK+kXJRSrh3oompqalIniE1G2GUymcjKyiIrK2vcTA+pN/A2pEgL6eeRlJSkW8ymB568Odvb28cY33jLawwODlJRURFS3YZ2Apv0wDSbzX79OZ1OJ7fddhtRUVFBIQmA9evXc/z48TGPXXDBBeq/165dy1NPPQW4CKS0tJSysjIAhBCdk76AEGNa0vF6SUIOq126dOmktw5Go3HMouzq6uLkyZPU1taqZde0tDRd3o2Dg4O6JncHAnfHJ9ngdPz4cZ/JUKluzMrKCmkZT+v6NH/+/DGOVO55DekitnTp0ilzplYUhfj4eLUTVm6h5BDllJQUNeK48847sdls/O53v5uyEvbDDz+sysDr6upQFIVNmzZhNpvZv3//rUKIn0/JhUwQISUKT23LcpH62/OeOHFCtWQLdqVA64EoM+xtbW0cOXJknOGNhAztjUYjpaWlIU16aRdlUVHRmDuloijqndxkMnHgwIFJDy2eCGJjY5kzZw5z5swZsygHBgaw2WwUFxdPq/eHu0Vhd3c3O3fu5L777sNkMvHf//3fahk21LjzzjsxmUx87nOfA1zk/vbbb7Nv3z5iY2OJi4vbrijKB0KI10N+MRPElEUUMmm5YMECWlpaqK2tVcN/7TAWp9OpJuRWrFgRcsbXztvUGt40NDSod/KUlBQOHTo0ZepPd2hnlco+ioMHD9Lb26tWT2QyeDogF2VcXBzV1dUUFRXR09NDQ0PDtOg13GE0GtVhxmeffTbf+c53ePHFF6mtrfXpoBYMPPLII7zwwgu8/vrr6u8nLy+P9evXa7UkLwErgBlLFCFNZkrfTE+VDbknl6rLpKQk1dAmOTnZo5HKVGNwcJDW1lYaGxvVITdTVXb1d12VlZXMnz9fnQmqzbuEakKXL8hBylp7f3ePiIk2100WQgjuv/9+ysvLeeyxx0ImgDt+/DibN29Wk5kvv/wyN910E//85z/HdOp2d3ezYcMG3n77bSIjI4mIiHgduE8I8WJILiwICDlRWK1WXfoIOVlKjoSb7rsQfLQgZRg9VWVXX5Bt6+77f3finco7eXd3N7W1tSxbtswniYZar+EJQgh++9vf8u677/K3v/0tZIK3q666ijfffJOOjg6ysrL40Y9+xF133YXFYlGT8GvXruX3v/89AH/961+56667UBSF6urqe4QQt4bkwoKEkBLFSy+9RFFRETk5OT7vcFKPsHjxYhITE8eoBmNjYz3mDEKNnp4eampqPCbkPJVdg6WF8AWttb+v/b+78jI6Olr9DIO9UDo7Ozly5AjLli0LKEpw12vIm0MwZ6kKIXjwwQfZvXs3O3funBZTZZ04vQVXDz74IH/5y19wOBxs3ryZbdu2jRu4c/LkSRoaGjzqEWTOoK2tjY6ODvULn5GREVLn6vb2durr6ykrK/O7zZBlVxn+p6amquF/MElD6jYCXZDwkVOWVK7KaGiyHZsygbl8+fJJEZC3VvTJ9MkIIfjzn//Mc889xzPPPDPp7WKggipwDVNevHgxP/zhD7n55pt9nf70Jgpw/cJOnjzJrl272LVrF4ODg2zevJlLLrmE5557jnXr1rFq1SpdX4jBwUFVCq0dbhPMu2RjYyNms3lC5jeyotPe3k5vby9JSUnqF34yd8nW1laamppYtmzZpAlSJkPb29ux2WxqNBSoHLq9vZ3jx48HfZqYNq/R0dEx4T6Uv/71rzz55JM899xzQam+BCKokrjssstQFIUzzjgjTBSBwmw28+STT3L33XeTnp7O5s2b2b59OwsWLAjoizo8PExbWxtmsxmDwaBGGhO9c8jRdSMjI6p/xmQgbe/b29vVOSNZWVkB3yWla5fsAA0m7Ha7Gg0NDg6OkWv7m/4uHdRDPZNkInmNJ598kkceeYQXX3wxqDkk92SlFk8//TRPPfUUjz76KADPPPMM77zzDnFxccTHx4eJYiK47bbbyM/P56qrruLZZ59l165dNDc3s2nTJrZv3x7wINqRkRE10nA6nWRkZJCVlaU7tHY6nWofSVFRUdDzDBPpJpUjFgcHBykpKQl5FUNGQ2azWfX+8GSv19raSnNzM8uWLZvy9nk9eY2nn36aBx54gBdeeIHExMSgvr4vorjkkku44oor+PznP8/AwAAbN25k9+7d3HvvvWGimCicTue4L35vby/PP/88u3bt4tixY2zcuJFt27ZRVlYW0CKxWq20t7fT3t6O3W73W52w2WyqbX2wHKD9QW6htEa02rKrEILa2lqEECxcuHBahgJp27xlbshut2M2mz22rk813PMar7zyCk6nk3/961+8/PLLIfEo8UYUd955J++//z67du1CURRuvvlm1qxZw+WXX84Pf/jDMFGECgMDA7z00ks89dRTHDp0iPPOO4+tW7eyevXqgEhDW52QtnXa/bh0pC4oKCArKyuE78g73Ltd09PT1YSenpEGU4HBwUEOHz6supNnZWXNGPs6cBHb73//e/70pz+p/T1PPvlk0EvXnojikUce4YEHHuD1119XcyHnnHMOJ06cAFzVM4PBwI9//GNuuOEGb6ee/l+yH8xIotBieHiYV155haeeeory8nLWr1/P1q1bOfPMMwPas2v340NDQyQkJNDd3c2SJUtmhEMWuEhj//79CCFUmflUlF39obGxkc7OTkpLS1XylfZ1E02GBhP/+Mc/+NGPfsRLL71ERkYGx48fn9AQIX/QK6jSIhxRTAMsFguvvfYaO3bs4P333+ess85i27ZtrFu3LqCkWkdHBzU1NSQkJKieC6EoaQYCOUksOzub3NzccS3o03WNx48fp7e316MJjrs3xHRc41tvvcX3vvc9XnzxxZD2uwQqqJIIE8U0w2az8cYbb/DUU0/x7rvvsnr1arZt28a5557rs1x38uRJGhsbKSsrIyoqSk3itbW10dfXF7JpYb5gsVjULZAnk2BPZVdpMhzKa5ROWXqSqe7XOBUzVPfs2cPNN9/MCy+8QG5ubkheY4oQJoqpgN1u51//+hdPPfUUb731FmVlZWzbto3zzjtvTIKwoaGBrq4uSktLPSbj3KeF6W0/nwykA5Revw33smsoTHxlxWV4eDjgCpR8vjYZKtu7gymU27dvH//5n//Jc889N+kkdCBiqt27d3PbbbdhtVqJjIzknnvu4bzzzpvs2wkTxVTD4XDw7rvvsnPnTl5//XUWL17Mli1b+Mc//sFnPvMZzj77bF1ffG37uXZB6jGS0YuBgQEqKysn7AClLbtqlauTkWpLPYnVamXx4sVBsefTKkO1ZjgTTYbu37+fb3zjGzzzzDNBMeoJREy1f/9+srKyyMnJoaqqik2bNtHc3DzZSwgTxXTC6XTyzjvvcN111xEVFUVRURFbt27lwgsvDMh1SS7ItrY29Q4pdRATvYvL/pZgmrtolatGozFgEZr03HQ6nSEry3qq8gRilFtZWclXvvIVdu7cyfz584N2XYGIqSSEEKSlpdHa2jrZPpIwUUw3/vCHPwBw3XXXUV5ezo4dO3j55ZfJy8tjy5YtXHTRRQGZvso7pOw/kZPMAgmr5UxSbUt2sOEuQpPVCW8lQyEEhw4dwmAwUFxcPCXJyECToQcPHuTaa6/lb3/7G4sWLQrqtegVU2nx1FNP8fvf/57XXnttsi8fJoqZCCEEVVVVPPXUU7z00kukp6ezdetWLr744oB9OT31n/jqSwhVj4QveBqDqC27CiE4ePAgkZGRIVGm6oGnwcvaqXB1dXVcffXVPProoyxdujTor69XTCVRXV3Nli1bePXVV5k3b95kXz5MFDMdUgX51FNP8cILL5CQkMCWLVu45JJLyMjICLj/RKpCPc0XaW5uprW1lbKyspD3SHiDLLu2tbWpd/HBwUESExNnjMBLOxXuzTffZNeuXbS0tPDwww+zfv36kLymXjEVuMyVzzvvPP70pz+xbt26YLz89H/ofnDaE4UWMtu/c+dOnnnmGaKiorjkkkvYunUrs2bNCmgRWSwWlTScTicmkwmHwxF0U97JwG63s3//ftV9bKrKroGgvr6e6667jmXLllFeXs5nP/tZXwrHCUOvmKqnp4dzzz2XO+64g09/+tPBevkwUZyqEELQ2NjIrl27ePrpp3E6nVxyySVs27aNvLw83aQhE4Q9PT1ERERgt9vHqBmnC06nU7X4Lygo8Fp2DWaVJ1A0Nzfzmc98ht/97neceeaZgGuLEuzrCURM9ZOf/IS77rprTCL11VdfneyQ7DBRfBwghKC1tVUljaGhIS6++GK2bt3K3LlzfVr81dTUYDAY1DZ6mcBrb29neHh4WmTaDoeDiooK0tLSPGoQvJVdQ20YpEVrayuXXXYZv/zlL0O23ZhBCBPFxw1CCMxmM08//TS7du2iu7ubT33qU2zbtm1MtcDpdFJVVUVcXJxXMpHuWDJfIMcLJiUlhYw0HA4H5eXlZGZm6p4DMtmya6Boa2vj0ksv5d577w2GmOlUwOlLFC+//DL/+Z//icPh4LrrruO2226b6KlmNDo7O1VPjZMnT3LBBRdwwQUX8Oijj3LTTTcxZ84cXeeRWf+2tjb6+/snNcnMG+SwoOzsbHJyciZ0juHhYbWCIoQI+sDljo4OPv3pT3PnnXeyadOmSZ8vUAu7u+66i4ceegij0cj//u//BuUadOD0JAqHw0FxcTG7d+8mLy+P1atX8/jjj7N48eKJXeUpgp6eHp544gl+8IMfMH/+fM466yy2b99OaWlpQMlB6ajd1tam9nZkZWVNqv/Ebrdz4MABcnNzyc7OntA53KEtu1qt1kl3knZ1dXHppZfyX//1X2zevDko1xiI6vLgwYNcddVV7N27l5aWFs4//3zq6uqmIkcz44kiJO4je/fupaioiLlz5wJw5ZVX8uyzz37siSI5OZkDBw7w0EMP8YlPfIKXXnqJ++67j7q6OtVTY9WqVX4Xu/t4QWnDL2eSZmVlqfoCPbDZbOpEsWD6bkRGRpKbm6t2u3Z0dFBfX6+WXbOysnRvo3p7e7n88su57bbbgkYSENhM0GeffZYrr7ySqKgoCgsLKSoqYu/evWoi9XRGSIiiubmZ/Px89f95eXn8+9//DsVLzTj87ne/UxfGFVdcwRVXXMHw8DAvv/wyf/zjH7nxxhs599xz2bp1K2vXrvW72BVFGTOIV+oL5PjDrKwsn5UJq9XKgQMHvHamBgsmk4lZs2Yxa9YsdRvV3NxMTU2N34iov7+fyy+/nG9961ts3749ZNfoCdqZoM3Nzaxdu1b9WV5eXsB9HFr3tumc3hZsTK+f2ccQnr4YMTExbN++ne3bt2OxWNi9ezd//etf+fa3v826detUTw1/fSPu4w9lZaK+vl7tP0lPT1crE1arlf379zNv3jzt+LqQQ9v4pe3I9TSlfXBwkCuvvJKvfe1rXH755VN2jTB+JmgwIEnCZrNNm6guFAgJUeTm5qpWYOBSsp3ifgFBQ1RUFJs3b2bz5s1YrVbVU+OWW25hzZo1bNu2jfXr1/uVdyuKQmJiIomJiRQVFakzUz/88EMiIyNJSUmhtbWV4uLigGXpwYTBYBgTEUkr/gMHDnDfffdhsVj4whe+MK6PItTwNBN0Mt/b559/nksuuQRw6TJ6e3u54IIL2LBhQ0gk51ONkCQz7XY7xcXFvP766+Tm5rJ69Woee+wxlixZMrGrPA1gt9t56623VE+N5cuXs23bNj75yU8GXIbs7u6moqKCiIgIoqKiVCn5TJqUNTIywuc//3mioqLUqfW//vWvQ/JaelWX1dXVfPazn1WTmRs2bODw4cN+t4f9/f1s3LiRsrIySkpKOHz4MJs2beLNN9/EZDJx2WWXsXLlSl+nmPH7k5CVR1966SW+9a1v4XA4uPbaa/ne97430VOddnA4HLzzzjvs3LmTf/zjHyxZsoRt27Zx/vnn+y1DSiOchQsXkpycrHaRynKmJI3pNMaVUcSmTZu44YYbUBSFoaGhoJVYtQjUwu7OO+/k4YcfxmQycf/99/OpT31K1+u0tLRwxx138Oabb/LGG2+Ql5fHvn37ePHFFxkZGWHz5s2cffbZ3p5++hJFGMGB0+lk37597Nixg927d1NUVMS2bdvYtGnTOAn40NAQ5eXlXo1wtFPC9IwyCAVsNhtf/OIXOeecc7jppptO6WTfoUOHyMvLIz4+nieeeIKzzz6bSy+9lKysLJ577jnAZbLz6KOPUlJSwjXXXOPtVDP+QzhliaKgoICEhASMRiMmk4n333+frq4urrjiCtWF+cknn5wxDtvBgNPp5MCBA+zYsYNXXnmF/Px81VOjoaGBlpYW1q1bp8uUR7ppt7W1qRqIrKyskE4Wt9vtXHvttaxcuZLbbrvtlCYJgFdeeYXrr7+eM888k5iYGB566CE6Ojq4/vrriYyM5LHHHgNcjW1+nLhm/AdxShPF+++/Pyabf+utt5Kamsptt93G3XffTXd395hZkB8nSE+NHTt2sHPnTrq6uvjGN77Bl770JVJTUwM6lxzsI/tP0tLSyMrKCmr/icPh4Gtf+xrFxcXccccdkz6vJ8Xljh07+OEPf0hNTQ179+5l1apVgIsUr7vuOj788EPsdjtXX301t99++4RfW1v2vPjii3n99dfZt28fS5cuxW6309fXxze/+U1OnjzJG2+84fF5bpjxRDEzeomDhGeffZYvfvGLAHzxi1/kmWeemd4LCiEURWHp0qV84QtfICIiggcffBCASy+9lG3btvHwww+r8zr9wWQykZ2dTVlZGatXryYxMZGGhgbee+89amtr6enp0XUeb3A4HNx4440UFBQEhSQArrnmGl5++eUxj5WUlLBr165xTWQ7duzAYrFQWVnJBx98wAMPPDBOhKUXTqdzzPXffPPN3HvvvXziE5/gvffew2QykZqayv/8z/+wZcuWMc89lSOoU1ZHoSgKF1xwAYqi8LWvfY2vfvWrtLW1qfLkWbNm0dbWNs1XGXrMmTNHdaK++OKL+f73v8/Ro0fZuXMnV111FVFRUWzZsoWtW7eSlZXl98tqNBrJysoiKysLp9NJZ2enKpxKTk4mKyvL7xBjLZxOJzfddBNpaWn85Cc/Cdpi8aS49GaPpygKg4OD2O12hoeHiYyMnNBcUq2Y6v7776e7u5tly5Zxww03kJqayqZNm3j77bd5+umniYuL4zvf+c64552qOGWJ4u233yY3N5f29nY2btzIwoULx/xcUZRTmsH1IjIyckyruKIoFBUV8d3vfpdbb72VxsZGdu7cyTXXXIMQQvXUyM3N9fv5yLmoUjgl+09qa2tJTExUpeTeFoHT6eS73/0u0dHR3HPPPdO2WC677DKeffZZsrOzGRoa4r777gt4ewYfial++ctf8vzzz7N+/XpefvllXn75ZX7xi19gt9u57rrrmDVrFjt27Bj3vFMZpyxRSCFMZmYm27dvZ+/evWRlZdHa2kp2djatra0hlSyfClAUhTlz5nDTTTfx7W9/m9bWVnbu3Mn111/P8PAwmzdvZuvWrRQWFuoiDW3/iVRbHj58mISEBHW2iNQcOJ1OfvCDH+BwOPjVr341rYtl7969GI1GWlpa6O7u5pxzzuH8889Xe5ECwXPPPcdDDz3EW2+9RXJyMjU1NTz++OM88sgjfPOb32TDhg3MmjULo9EYEpOd6cIpSXWDg4P09/er/3711VcpKSlhy5Yt/PnPfwbgz3/+M1u3bp3Oy5xRUBSFnJwcbrzxRl5//XWeeeYZ0tLSuOmmm/jkJz/Jz3/+c+rq6nTlIhRFISUlhQULFrB27Vry8/Pp7e1l3759vPPOO/zxj3/kjjvuoLu7m9/85jfTfkd97LHHuPDCC4mIiCAzM5N169bx/vvv63qu++dhMBgwm8385Cc/AVzbnZKSEv75z38CrhuY0WhECPGxIQk4RYmira2Ns88+m7KyMtasWcPFF1/MhRdeyG233cbu3buZP38+r7322sfWA2OyUBSFrKwsrr/+el599VX+/ve/k5+fz/e+9z3OPfdc7rzzTg4ePKibNJKSkpg/fz5nnHEGubm5vPrqq/zlL3+ho6OD559/fgrekW/Mnj2bf/zjH4DrxvLee++N26p6grZKUVNTQ319PZs3b+bRRx+lvr6eH/3oR4Dr+2ixWLDZbOpzP27b3lO2PBpGaNDT08Nzzz3Hrl27aGhoYOPGjWzfvt3jkGJ3CCG47777qKio4LHHHuPo0aOcOHGC888/P+jX6UlxmZqayo033ojZbCY5OZlly5bxyiuvMDAwwJe+9CWV/L70pS9xyy236H6tu+++mzfffBObzUZpaSk33ngjjY2NfOUrXyE5OZl58+bxhz/8gcTExIl2jM54VgkTRRhe0d/fz4svvsjOnTupq6tjw4YNbN26lZUrV44jDSEEv/nNb9izZw9/+9vfpmxmSbDhvtDfeOMNfv7zn/P3v/+d66+/no6ODp566insdjt79uzh4YcfZvbs2Wp0MUHMeKI4JbceU4Frr72WzMxMSkpK1Me6urrYuHEj8+fPZ+PGjXR3dwOuL9d//Md/UFRURGlpKR9++OF0XXZQkZCQwJVXXsmOHTvYs2cP69at44EHHuCss87iu9/9Lu+++65q9f/HP/6Rt956iyeeeOKUJQkYv2WIiIhg48aN/OAHP6ChoYG//e1vAPz73//m7LPP5itf+Qp79uzhT3/603Rc7tRBCOHrz2mLf/7zn+KDDz4QS5YsUR+75ZZbxF133SWEEOKuu+4St956qxBCiBdffFFceOGFwul0ij179og1a9ZMyzVPFYaHh8Vzzz0nrr76alFSUiLWr18vPvnJT4rh4eGgnP9LX/qSyMjIGPPZP/nkk2Lx4sVCURSxb9++MceXl5eLtWvXisWLF4uSkpIJX0d5ebm48847xS233CJuvfVWYbFYRF1dndi0aZM4//zzRX9/vxBCiF/96lfi7LPPFsPDw8JisYiqqirhdDon/oZ9r8EZ8SdMFD5QX18/5staXFwsWlpahBBCtLS0iOLiYiGEEF/96lfFY4895vG4jzssFov4xS9+Icxmc9DO6YmkDx48KA4dOiTOPffcMURhs9nE0qVLxYEDB4QQQnR0dAi73R7wa5aXl4sFCxaIe+65R9x5551i+/btoqSkRDQ0NIjf/va34oorrhA//vGPxX//93+LkpIScejQocm/0Y8w7UTg788pq6OYDnhTfnqy/mtubg6aie1MRmRkJN/+9reDes5AVJevvvoqpaWllJWVAUzIpEdOSL/33nvH+HVef/31bN68mffff5+ioiKqq6vp7e1l586dFBcXf6x0Ev4QJooJ4nRRfs501NXVoSgKmzZtwmw2c+WVV3LrrbcGdI633nqL2tpalSSkN8bvf/97rrzySr73ve9xzz33sHHjRvU5TqfztCEJCBNFQPCm/Axb/00f7HY7b7/9Nvv27SM2NpYNGzawcuVKNmzYoPsc3/zmNxkcHGTdunW88847xMbGMjw8TExMDOvXrx8zuFiMVkWmW0Q21Ti93u0k4U35uWXLFv7v//4PIQTvvfceSUlJp8W2YyYgLy+P9evXk56eTmxsLBdddNGEqk633norl19+OWeccQaA6gAWGxtLfHw8TqfzY+WqHTD8JDFOW1x55ZVi1qxZwmQyidzcXPHggw+Kjo4Ocd5554mioiKxYcMG0dnZKYQQwul0im984xti7ty5oqSkZFxWPozA4Z5IlnBPZnZ1dYnly5eLwcFBYbPZxIYNG8QLL7ww4df97W9/K5YvXy6EEOKDDz4QxcXF4u23357w+XRi2pOV/v6EiSKMGQdPJL1r1y6Rm5srIiMjRWZmprjgggvU4//yl7+IxYsXiyVLlohbbrll0q//8MMPi5SUFFFaWipee+01IYSYbPnTH6adCPz9CRPFDIIn/cAdd9whcnJyRFlZmSgrKxMvvvii+rOf/vSnYt68eaK4uFi8/PLL03HJH1vs2LFDPPPMM0KIkJOEEDOACPz9CUu4ZxA8zcn84Q9/SHx8PDfffPOYY6dxTuZpBSGmJC8x4xMf4WTmDML69et1G6p4m5MZRnBx2iYv3RAmilMAv/71ryktLeXaa69V+0u8ibxmMjz1z+zYsYMlS5ZgMBg8ekQ0NjYSHx/PvffeO5WXGoYbTiuieOaZZ1AUhUOHDk33pejG17/+dY4ePcqBAwfIzs5WfRhPRQRiiCtx00036R7CE0bocFoRxeOPP87ZZ5/N448/Pt2XohtZWVkYjUYMBgNf+cpX1O3FqSjy8rS1WrRoEQsWLPB4/DPPPENhYWF4FOUMwGlDFAMDA7z99ts89NBDPPHEE4BLhvuNb3yDhQsXsnHjRi666CKeeuopAD744APOPfdcVq5cyaZNm2htbZ2W69a+7tNPP62G7Vu2bOGJJ57AYrFQX1/P4cOHWbNmzbRcYygwMDDAz372M+64447pvpQwwG959GPzB/gc8NDov98FVgKXAS/hIsxZQPfoYxGjx2SMHn8F8PAUXOPjQCtgA5qALwN/ASqBCuA5IFtz/PeAo0At8Knp/ox1vscCoMrD428CqzT/vxe4fPTfPwRunu5rP53/nE69HlcBvxz99xOj/zcBO4QQTuCkoihyrNMCoATYPZr1NuJawCGFEOIqDw8/5OP4O4E7Q3dF04ozgMsURfk5kAw4FUUZEUKEZuR5GD5xWhCFoiipwHnAUkVRBK6FL4CnvT0FqBZCnDlFlxiGG4QQ58h/K4ryQ2AgTBLTh9MlR3EZ8BchxBwhRIEQIh+oB7qASxVFMSiKkgV8YvT4WiBDUZQzARRFiVAUJZxRmyQURXkc2AMsUBSlSVGULyuKsl1RlCbgTOBFRVFemd6rDMMT/CkzPxYY3VL8TAjxsuax/wAW4YoePgGcGP33z4QQuxVFWQb8L5CEK/K6Xwjxxym+9DDCmBE4LYjCFxRFiRdCDCiKkgbsBdYJIU5O93WFEcZMwmmRo/CDFxRFSQYigf8Ok0QYYYzHaR9RhBFGGP5xuiQzwwgjjEkgTBRhhBGGX4SJIowwwvCLMFGEEUYYfhEmijDCCMMv/j8aQGkrCA7JswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "index=23\n",
    "ax.plot(nx[:][index],nz[:][index],T[:][index])\n",
    "ax.set_xlabel('Age')\n",
    "ax.set_ylabel('Depth')\n",
    "ax.set_zlabel('Temperature')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "eb617d57",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,\n",
       "         55.,  60.,  65.,  70.,  75.,  80.,  85.,  90.,  95., 100., 105.,\n",
       "        110., 115., 120., 125., 130., 135., 140., 145., 150., 155., 160.,\n",
       "        165., 170., 175., 180., 185., 190.]),\n",
       " array([55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55.,\n",
       "        55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55.,\n",
       "        55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55., 55.]),\n",
       " array([4.3982029 , 4.40572834, 4.41502047, 4.41996193, 4.42710114,\n",
       "        0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "        0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "        0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "        0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "        0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "        0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "        0.        , 0.        , 0.        , 0.        ]),\n",
       " array([ 774.37  ,  751.02  ,  727.6765,  712.1095,  688.762 , 1210.195 ,\n",
       "        1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 ,\n",
       "        1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 ,\n",
       "        1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 ,\n",
       "        1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 ,\n",
       "        1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 , 1210.195 ,\n",
       "        1210.195 , 1210.195 , 1210.195 ]))"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nx[:][10],nz[:][10],V[:][10],T[:][10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "64b90ee8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x233eed4b908>]"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVK0lEQVR4nO3df4zkdX3H8dfr9vbu9kBDzW0Nctg79KQeaBTXKyRUiY31MFp+VJo702hTIlAlbZMaBZtUbCSptlWrUturnqfVQomVSgL0pK3JtQaUPUS448QuoLIneotUKc7szezMu3/Md85x3V/Hzvf7+dzO85Fsdubz/e7se75hX3zu8/3M5+OIEABgMKxKXQAAoDqEPgAMEEIfAAYIoQ8AA4TQB4ABsjp1AYvZsGFDbNq0KXUZAHDC2L9//xMRMTrXsexDf9OmTRofH09dBgCcMGx/d75jDO8AwAAh9AFggBD6ADBACH0AGCCEPgAMEEIfAAYIoQ8AAyT7efor1d6DP9DBwz9JXQaATK1fu1pXvfoFfX9dQj+R93zxAf3opw3ZqSsBkKMNJ68l9FeSp4/O6MpXnaFrX//i1KUAGCCM6SfQboeOzrS1bngodSkABgyhn0C92ZIkjawh9AFUi9BPoBv66wl9ABUj9BOoNzqhz/AOgKoR+gkcG94h9AFUjNBPoNvTZ3gHQNUI/QRqDXr6ANIg9BOYZvYOgEQI/QSYsgkgFUI/AYZ3AKRC6CdATx9AKoR+AtP09AEkQugnwPAOgFQI/QTqzZbWDK3S6iEuP4BqLZo6tnfbPmL7QE/bX9r+lu37bd9i+5SeY9fanrD9kO3X9bRvL9ombF/T93dyAqk3ZrRumMAHUL2lJM8eSdtntd0p6eyIeKmkb0u6VpJsb5W0Q9JZxc/8re0h20OSbpB0oaStknYW5w6kerOl9WvYygBA9RYN/YjYJ+nJWW1fjoiZ4undkjYWjy+SdFNEHI2IRyVNSNpWfE1ExCMR0ZB0U3HuQKo328zcAZBEP8YYfl/SHcXj0yQ91nNssmibr31Otq+wPW57fGpqqg8l5qUzvEPoA6jeskLf9p9KmpH0+f6U0xERuyJiLCLGRkdH+/nSWegM7xD6AKr3jAeWbf+epDdI+o2IiKL5sKTTe07bWLRpgfaBU28wpg8gjWfU07e9XdK7JP1WRNR6Dt0qaYfttbY3S9oi6euS7pG0xfZm22vUudl76/JKP3HVGi2GdwAksWh30/aNki6QtMH2pKT3qjNbZ62kO21L0t0RcVVEHLR9s6QH1Rn2eUdEtIrXuVrSXklDknZHxMES3s8JYZrhHQCJLBr6EbFzjuZPLXD+9ZKun6P9dkm3H1d1K1St0eLTuACS4BNCCdSbLaZsAkiC0E9gmtAHkAihX7Fmq61mKxjeAZAEoV+x7lr63MgFkAKhX7HuWvpM2QSQAqFfMdbSB5ASoV8xhncApEToV6zb019H6ANIgNCv2HS3p8/wDoAECP2K1btj+vT0ASRA6Fes1uRGLoB0CP2KTdPTB5AQoV+xOj19AAkR+hWr0dMHkBChX7FuT3/dakIfQPUI/Yp1NkVfpVWrnLoUAAOI0K9YZ1N09scFkAahX7F6o81NXADJEPoVqzc7wzsAkALpU7F6g+EdAOkQ+hWrN9kUHUA6hH7F6o0WK2wCSIbQr1i92WKFTQDJEPoVqzVafBoXQDKEfsWmmy32xwWQDKFfsc7sHUIfQBqEfoUiQjVm7wBIaNHQt73b9hHbB3raLrN90Hbb9lhP+7Dtz9h+wPYh29f2HNtu+yHbE7av6f9byd/RmbYiWGETQDpL6envkbR9VtsBSZdK2jer/TJJayPiJZJeIelK25tsD0m6QdKFkrZK2ml763IKPxFNs5Y+gMQWDf2I2CfpyVlthyLioblOl3SS7dWSRiQ1JD0laZukiYh4JCIakm6SdNFyiz/RsJY+gNT6Pab/BUk/lfS4pO9J+quIeFLSaZIe6zlvsmibk+0rbI/bHp+amupziel019LnRi6AVPod+tsktSQ9T9JmSX9i+4zjfZGI2BURYxExNjo62ucS06kXPX2mbAJIpd+h/2ZJ/xYRzYg4IumrksYkHZZ0es95G4u2gcL+uABS63fof0/SayTJ9kmSzpX0LUn3SNpie7PtNZJ2SLq1z787e92ePsM7AFJZypTNGyXdJelM25O2L7d9ie1JSedJus323uL0GySdbPugOkH/6Yi4PyJmJF0taa+kQ5JujoiDZbyhnNUY3gGQ2KILu0fEznkO3TLHuU+rM21zrte5XdLtx1XdCjPNjVwAifGJ3AodG9Mn9AEkQuhX6Ng8fYZ3ACRC6Fdomp4+gMQI/QrVGjNaZWnNEJcdQBqkT4XqjbZGhodkO3UpAAYUoV+herOlkTWLTpgCgNIQ+hWqN2Y0soZLDiAdEqhCnU3R6ekDSIfQr1C92dY6Zu4ASIjQr1C9MaORYS45gHRIoArVmy2t50YugIQI/QrVGmyKDiAtQr9C040WK2wCSIrQr1BneIfQB5AOoV+hWqPFujsAkiL0K9Juh47OtBnTB5AUoV+R6RlW2ASQHqFfEdbSB5ADQr8i3U3R6ekDSInQr8ixDVTo6QNIiNCvCMM7AHJA6Fekuyk68/QBpEToV6Q7ps8qmwBSIvQrQk8fQA4I/YrUGdMHkAFCvyI1Zu8AyAChX5Fp5ukDyMCioW97t+0jtg/0tF1m+6Dttu2xWee/1PZdxfEHbK8r2l9RPJ+w/VHb7v/byVd3TJ+llQGktJSe/h5J22e1HZB0qaR9vY22V0v6nKSrIuIsSRdIahaHPyHpbZK2FF+zX3NFqzVaGh6yhof4xxWAdBZNoIjYJ+nJWW2HIuKhOU7/TUn3R8Q3i/N+FBEt26dKenZE3B0RIemzki5edvUnkOkmu2YBSK/f3c4XSQrbe23fa/tdRftpkiZ7zpss2uZk+wrb47bHp6am+lxiGrXGDOP5AJLr9y7dqyWdL+mVkmqS/sP2fkk/OZ4XiYhdknZJ0tjYWPS5xiTqzTabogNIrt89/UlJ+yLiiYioSbpd0jmSDkva2HPexqJtYNTZHxdABvod+nslvcT2+uKm7qslPRgRj0t6yva5xaydt0j6Up9/d9bqzRmNDHMTF0BaS5myeaOkuySdaXvS9uW2L7E9Kek8SbfZ3itJEfG/kj4k6R5J90m6NyJuK17q7ZI+KWlC0sOS7uj3m8lZvdFieAdAcoumUETsnOfQLfOc/zl1pm3Obh+XdPZxVbeC1JttPeckhncApMV4Q0XqzN4BkAFCvyL1ZkvruZELIDFCvyK1RouePoDkCP2KTDcJfQDpEfoVaLbaaraCZRgAJEfoV6DOWvoAMkHoV4C19AHkgtCvAD19ALkg9CtQo6cPIBOEfgWO9fQJfQCJEfoVqDcY3gGQB0K/AoQ+gFwQ+hXoDu+sZ3gHQGKEfgW6PX02UQGQGqFfAXr6AHJB6FeA2TsAckHoV6A7T3/dakIfQFqEfgWmmy2tG16lVaucuhQAA47Qr0CtMcN0TQBZIPQrUG+0CX0AWSD0K8AGKgByQehXoMam6AAyQehXoLMp+urUZQAAoV+FerOtdfT0AWSA0K9AvTGjkWEuNYD0SKIK1JstrV/D8A6A9Aj9CtQbLRZbA5CFRUPf9m7bR2wf6Gm7zPZB223bY3P8zPNtP237nT1t220/ZHvC9jX9ewv5qzdazNMHkIWl9PT3SNo+q+2ApEsl7ZvnZz4k6Y7uE9tDkm6QdKGkrZJ22t56vMWeiCKiGN4h9AGkt+hAc0Tss71pVtshSbJ/cS0Z2xdLelTST3uat0maiIhHinNuknSRpAefYd0njKMzbbWDFTYB5KGvY/q2T5b0bknvm3XoNEmP9TyfLNpWvOkmWyUCyEe/b+ReJ+nDEfH0cl7E9hW2x22PT01N9aeyRFhLH0BO+j2P8Nckvcn2ByWdIqlte1rSfkmn95y3UdLh+V4kInZJ2iVJY2Nj0ecaK1VjU3QAGelr6EfEr3cf275O0tMR8XHbqyVtsb1ZnbDfIenN/fzdueruj0tPH0AOljJl80ZJd0k60/ak7cttX2J7UtJ5km6zvXeh14iIGUlXS9or6ZCkmyPi4PLLz1+dMX0AGVnK7J2d8xy6ZZGfu27W89sl3b7kylYIevoAcsIncktGTx9ATgj9ktHTB5ATQr9k3Z4+n8gFkANCv2R1pmwCyAihX7JuT59VNgHkgNAvWb3R0ipLa1dzqQGkRxKVrN7sLKs81+J0AFA1Qr9ktUaLmTsAskHol2y6SegDyAehX7JaY4aZOwCyQeiXrN5sa4RN0QFkgtAv2XSjpZFhLjOAPJBGJas1Gd4BkA9Cv2T1RkvrGd4BkAlCv2TTzTafxgWQDUK/ZLXGjEbWcJkB5IE0Klm9yfAOgHwQ+iVqt4PhHQBZIfRLND3DWvoA8kLol4i19AHkhtAvUY3QB5AZQr9E0032xwWQF0K/RN1ds+jpA8gFoV+iY8M79PQBZILQL1Gd4R0AmSH0S8TsHQC5IfRL1A195ukDyAWhXyJu5ALIzaKhb3u37SO2D/S0XWb7oO227bGe9tfa3m/7geL7a3qOvaJon7D9Udvu/9vJS7env46ePoBMLKWnv0fS9lltByRdKmnfrPYnJL0xIl4i6a2S/rHn2CckvU3SluJr9muuOPT0AeRm0dCPiH2SnpzVdigiHprj3G9ExPeLpwcljdhea/tUSc+OiLsjIiR9VtLFy64+c/VmS8ND1vAQo2gA8lBmGv22pHsj4qik0yRN9hybLNrmZPsK2+O2x6empkossVz1RosVNgFkpZTQt32WpA9IuvKZ/HxE7IqIsYgYGx0d7W9xFepslUjoA8hH30Pf9kZJt0h6S0Q8XDQflrSx57SNRduKVmu2GM8HkJW+hr7tUyTdJumaiPhqtz0iHpf0lO1zi1k7b5H0pX7+7hzVGy2NsGsWgIwsZcrmjZLuknSm7Unbl9u+xPakpPMk3WZ7b3H61ZJeKOnPbN9XfP1yceztkj4paULSw5Lu6Pebyc10s6WRYW7iAsjHot3QiNg5z6Fb5jj3/ZLeP8/rjEs6+7iqO8F1NkVneAdAPuiGlqjebGtkmOEdAPkg9Es03WzR0weQFUK/RLXGDGP6ALJCIpWoM0+f4R0A+SD0S1Rv8olcAHkh9EvSbLXVbAUfzgKQFUK/JNNNNlABkB9CvySspQ8gR4R+Sbpr6a9neAdARgj9khzbQIWePoCMEPolqTXYNQtAfgj9kkw36OkDyA+hXxJ6+gByROiXhDF9ADki9EtyLPTp6QPICKFfkjpj+gAyROiXpM4ncgFkiNAvybFP5K4m9AHkg9AvSb3Z0trVq7RqlVOXAgDHEPol6aylTy8fQF4I/ZLUGi1m7gDIDqFfkulmixU2AWSH0C9JvcnwDoD8EPol6WyKTugDyAuhX5J6s60RNkUHkBlCvyTTjZZGhrm8APJCKpWk1mR4B0B+Fh1/sL1b0hskHYmIs4u2yyRdJ+nFkrZFxHjP+ddKulxSS9IfRsTeon27pL+RNCTpkxHxF/19Kz/vjR/772Obk6fw/R9P6/wXEvoA8rKUQec9kj4u6bM9bQckXSrp73tPtL1V0g5JZ0l6nqR/t/2i4vANkl4raVLSPbZvjYgHl1X9Al4wepIarXZZL7+oFz33Wbr0nI3Jfj8AzGXR0I+IfbY3zWo7JEn2LywxcJGkmyLiqKRHbU9I2lYcm4iIR4qfu6k4t7TQ/8iOl5f10gBwwur3mP5pkh7reT5ZtM3XPifbV9getz0+NTXV5xIBYHBleSM3InZFxFhEjI2OjqYuBwBWjH5PJD8s6fSe5xuLNi3QDgCoSL97+rdK2mF7re3NkrZI+rqkeyRtsb3Z9hp1bvbe2uffDQBYxFKmbN4o6QJJG2xPSnqvpCclfUzSqKTbbN8XEa+LiIO2b1bnBu2MpHdERKt4nasl7VVnyubuiDhYxhsCAMzPEZG6hgWNjY3F+Pj44icCACRJtvdHxNhcx7K8kQsAKAehDwADJPvhHdtTkr77DH98g6Qn+lhOv1Hf8lDf8lDf8uRc369ExJzz3bMP/eWwPT7fuFYOqG95qG95qG95cq9vPgzvAMAAIfQBYICs9NDflbqARVDf8lDf8lDf8uRe35xW9Jg+AODnrfSePgCgB6EPAANkRYa+7e22H7I9Yfua1PXMZvs7th+wfZ/tLNaYsL3b9hHbB3ranmP7Ttv/U3z/pczqu8724eI63mf79YlqO932V2w/aPug7T8q2rO4fgvUl8X1K2pZZ/vrtr9Z1Pi+on2z7a8Vf8v/XCzYmFN9e2w/2nMNX5aivuMSESvqS50F3R6WdIakNZK+KWlr6rpm1fgdSRtS1zGrpldJOkfSgZ62D0q6pnh8jaQPZFbfdZLemcG1O1XSOcXjZ0n6tqStuVy/BerL4voVdVnSycXjYUlfk3SupJsl7Sja/07SH2RW3x5Jb0p9/Y7nayX29Lep2JoxIhqSulszYgERsU+d1VN7XSTpM8Xjz0i6uMqaes1TXxYi4vGIuLd4/H+SDqmzM1wW12+B+rIRHU8XT4eLr5D0GklfKNpTXsP56jvhrMTQP66tGRMJSV+2vd/2FamLWcBzI+Lx4vEPJD03ZTHzuNr2/cXwT7Lhp65iP+mXq9MTzO76zapPyuj62R6yfZ+kI5LuVOdf7D+OiJnilKR/y7Pri4juNby+uIYftr02VX1LtRJD/0RwfkScI+lCSe+w/arUBS0mOv+uza1n8wlJL5D0MkmPS/rrlMXYPlnSv0j644h4qvdYDtdvjvqyun4R0YqIl6mzs942Sb+asp7ZZtdn+2xJ16pT5yslPUfSu9NVuDQrMfQX2rIxCxFxuPh+RNIt6vwHnqMf2j5VkorvRxLX83Mi4ofFH2Jb0j8o4XW0PaxOoH4+Ir5YNGdz/eaqL6fr1ysifizpK5LOk3SK7e5mT1n8LffUt70YOouIOCrp08rkGi5kJYZ+1lsz2j7J9rO6jyX9pqQDC/9UMrdKemvx+K2SvpSwll/QDdTCJUp0HW1b0qckHYqID/UcyuL6zVdfLtevqGXU9inF4xFJr1Xn3sNXJL2pOC3lNZyrvm/1/E/d6txvyPVv+ZgV+YncYurZR/SzrRmvT1vRz9g+Q53evdTZrvKfcqjPPdtiSvqhOtti/qs6syeer87y1r8TEUlups5T3wXqDE2EOjOiruwZQ6+ytvMl/ZekByS1i+b3qDNunvz6LVDfTmVw/STJ9kvVuVE7pE5n9OaI+PPi7+UmdYZOviHpd4tedS71/ac628Za0n2Sruq54ZulFRn6AIC5rcThHQDAPAh9ABgghD4ADBBCHwAGCKEPAAOE0AeAAULoA8AA+X/2IekGSmxXhQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(T[1][:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff50998b",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}