mesh_generator_3d_salt.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt #debug\n",
    "from scipy import interpolate\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def triangle_weight(Nx, begin=None, end=None, peak=None):\n",
    "    if(begin == None):\n",
    "        begin = 0\n",
    "    if(end   == None):\n",
    "        end   = Nx\n",
    "    if(peak  == None):\n",
    "        peak  = int((end+begin)/2)\n",
    "    if(peak < begin or peak > end):\n",
    "        raise ValueError(\"Peak outside triangle bounds\")\n",
    "    if(begin >= end):\n",
    "        raise ValueError(\"Incorrect bounds order: begin >= end\")\n",
    "    triangle = np.hstack((np.zeros(begin),\n",
    "                          np.linspace(0,1,peak-begin, endpoint=False),\n",
    "                          np.linspace(1,0,end-peak),\n",
    "                          np.zeros(Nx-end)))\n",
    "    return triangle\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def asym_triangle(Nx, begin=None, end=None, peak=None):\n",
    "    if(begin == None):\n",
    "        begin = 0\n",
    "    if(end   == None):\n",
    "        end   = Nx\n",
    "    if(peak  == None):\n",
    "        peak  = int((end+begin)/2)\n",
    "    triangle = triangle_weight(Nx, begin=begin, end=end, peak=peak-1)\n",
    "    triangle[peak:] = np.hstack((0,-triangle[peak:-1]))\n",
    "    return triangle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "448.7007\n",
      "683.9578777522539\n",
      "1128.0434\n"
     ]
    }
   ],
   "source": [
    "N_file = 321\n",
    "model1_X = np.reshape(np.fromfile(\"model_2-X.txt\", sep=' '), (N_file,N_file)) #read data from file and convert it to 2d array\n",
    "model1_Y = np.reshape(np.fromfile(\"model_2-Y.txt\", sep=' '), (N_file,N_file))\n",
    "model1_Z = np.reshape(np.fromfile(\"model_2-Z.txt\", sep=' '), (N_file,N_file))\n",
    "\n",
    "print(np.amin(model1_Z))\n",
    "print(model1_Z.mean())\n",
    "print(np.amax(model1_Z))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mesh size\n",
    "#(Nx,Ny,Nz) = (161,161,61)\n",
    "(Nx,Ny,Nz) = (51,51,51)\n",
    "#mesh SIZE\n",
    "Lxy = 8000\n",
    "Lz_rock  = 2000\n",
    "Lz = Lz_rock+model1_Z.mean()\n",
    "\n",
    "#boundaries of the curved region\n",
    "z_start = 0\n",
    "z_stop  = Nz\n",
    "\n",
    "#ocean floor\n",
    "k_border = round(Nz * model1_Z.mean() / Lz )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generate an empty array for vertices\n",
    "vertices = np.zeros((Nx,Ny,Nz,3))\n",
    "\n",
    "#generate xy values\n",
    "x = np.linspace(0,1,Nx)\n",
    "y = np.linspace(0,1,Ny)\n",
    "z = np.linspace(0,1,Nz)\n",
    "\n",
    "x = x * Lxy - Lxy / 2\n",
    "y = y * Lxy - Lxy / 2\n",
    "z = z * Lz\n",
    "\n",
    "#populate with rect mesh\n",
    "for k in range(0,Nz):\n",
    "    for j in range(0,Ny):\n",
    "        for i in range(0,Nx):\n",
    "            vertices[i,j,k,:] = (x[i], y[j], z[k])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n"
     ]
    }
   ],
   "source": [
    "ot_file = \"Saltf@@\"\n",
    "(ot_Nx,ot_Ny,ot_Nz)=(676,676,210)\n",
    "\n",
    "ot_x = np.linspace(-Lxy/2,Lxy/2,ot_Nx)\n",
    "ot_y = np.linspace(-Lxy/2,Lxy/2,ot_Ny)\n",
    "ot_z = np.linspace(0,Lz_rock,ot_Nz)\n",
    "\n",
    "ot_c=np.reshape(np.fromfile(ot_file,dtype='>f4'),(ot_Nz,ot_Ny,ot_Nx))\n",
    "f_ot_c=interpolate.RegularGridInterpolator((ot_z,ot_y,ot_x),ot_c)\n",
    "interp_ot_c = np.empty((Nz-k_border,Ny,Nx))\n",
    "for k in range(0,Nz-k_border):\n",
    "    print(k)\n",
    "    for j in range(0,Ny):\n",
    "        for i in range(0,Nx):\n",
    "            interp_ot_c[k,j,i] = f_ot_c((z[k], y[j], x[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#plt.figure(figsize=(16,8))\n",
    "#plt.imshow(ot_c[:,-1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#plt.figure(figsize=(16,8))\n",
    "#plt.imshow(interp_ot_c[:,-1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_dz = interpolate.interp2d(model1_X[0,1:-3], model1_Y[1:-3,0], model1_Z[1:-3,1:-3])\n",
    "\n",
    "#calculate offsets\n",
    "#xygrid = np.transpose(np.array(np.meshgrid(x,y)), (2,1,0)) #iirc this makes a 2d array of (x,y) pairs\n",
    "#dz = f(xygrid[:,:,0], xygrid[:,:,1], 0.7)/2 #feed that array into the bulge function\n",
    "interp_dz = f_dz(x, y)  #interpolate heightmap to target mesh size\n",
    "dz_mean = np.mean(interp_dz)\n",
    "#k_border = np.array(np.where(z < dz_mean)).flatten()[-1] #closest* mesh z-plane to dz_mean\n",
    "#print(k_border)\n",
    "interp_dz -= z[k_border]\n",
    "weights = triangle_weight(Nz, z_start, z_stop, k_border)\n",
    "#weights_condense = asym_triangle(Nz, z_start, z_stop, k_border)\n",
    "weights_condense = 1-weights\n",
    "#weights_condense[[0,-1]] = 0\n",
    "weights_condense[k_border:] = -weights_condense[k_border:]\n",
    "weights_condense = (np.abs(weights_condense)**4)*np.sign(weights_condense)\n",
    "\n",
    "#print(weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6eba752828>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkyUlEQVR4nO3deXxU9b3/8ddnJisJIStrwo4iKoJEkKKt4lKsrWjd0C7Uarm3rdXW28Uuv9tevb3V9l6199YuVK20teJWFe2iqLhVEQKCshNAgbCFfU9I5vP7Yw52iAnbTDKZzPv5eMxjzvmeM3M+B2PeOed7zveYuyMiIukrlOwCREQkuRQEIiJpTkEgIpLmFAQiImlOQSAikuYykl3A8SgtLfW+ffsmuwwRkZQyZ86cze5e1rQ9JYOgb9++VFVVJbsMEZGUYmbvN9euU0MiImlOQSAikuYUBCIiaU5BICKS5hQEIiJpLiFBYGYPmNkmM1vQwnIzs/81s2oze8fMTo9ZNtHMlgeviYmoR0REjl6ijggeBMYdZvlFwKDgNQn4FYCZFQM/BEYBI4EfmllRgmoSEZGjkJAgcPdXga2HWWU88HuPmgkUmlkP4OPAdHff6u7bgOkcPlDi8vS8Gv44s9nLaEVE0lZb9RH0AtbEzK8N2lpq/xAzm2RmVWZWVVtbe1xFPLdwA795dcVxfVZEpKNKmc5id5/s7pXuXllW9qE7pI/K6b2LWLN1H5t27U9wdSIiqautgqAGqIiZLw/aWmpvFaf3iXY/zH1/e2ttQkQk5bRVEEwDPh9cPXQmsMPd1wPPAReaWVHQSXxh0NYqTu5ZQFY4xNzV21prEyIiKSchg86Z2cPAOUCpma0leiVQJoC7/xr4K/AJoBrYC1wXLNtqZrcDs4Ovus3dD9fpHJfsjDCnlndh7vsKAhGRgxISBO5+zRGWO/DVFpY9ADyQiDqOxum9C5ny5vvUN0TIykiZLhIRkVaTdr8JR/Qpor4hwsJ1O5JdiohIu5B2QXB672iH8RydHhIRAdIwCLoW5FBelMvbq7cnuxQRkXYh7YIAokcFOiIQEYlKyyAY0aeIDTv3s277vmSXIiKSdGkZBOonEBH5p7QMgsE9OpObGVYQiIiQpkGQGQ4xtLwLb+sOYxGR9AwCiPYTLFy3k/0HGpNdiohIUqV1EDREnHfW6sYyEUlvaRsEw9VhLCICpHEQFOdl0b80TyORikjaS9sggOhRwdz3txEdE09EJD2ldRCM6FPElj31rN66N9mliIgkTVoHwel9CgH1E4hIekvrIBjUtTOdszPUTyAiaS0hQWBm48xsqZlVm9mtzSy/28zmBa9lZrY9ZlljzLJpiajnaIVDxrDehczRM4xFJI3F/YQyMwsD9wIXAGuB2WY2zd0XHVzH3b8Rs/7XgOExX7HP3YfFW8fxOr13Ef/30nJ21zWQn52QB7aJiKSURBwRjASq3X2lu9cDU4Hxh1n/GuDhBGw3IU7vU0TEYf6a7ckuRUQkKRIRBL2ANTHza4O2DzGzPkA/4KWY5hwzqzKzmWZ2aUsbMbNJwXpVtbW1CSg7alhFIWbqMBaR9NXWncUTgMfdPXaAnz7uXglcC9xjZgOa+6C7T3b3SnevLCsrS1hBXXIzGdQ1Xx3GIpK2EhEENUBFzHx50NacCTQ5LeTuNcH7SuBlDu0/aBOVfYupem8bBxojbb1pEZGkS0QQzAYGmVk/M8si+sv+Q1f/mNlgoAh4M6atyMyyg+lSYAywqOlnW9vHTihjd12DTg+JSFqKOwjcvQG4EXgOWAw86u4Lzew2M7skZtUJwFQ/dDyHk4AqM5sPzADuiL3aqK2MGVhKZtiYsXRTW29aRCTpLBXH2amsrPSqqqqEfue1v53Jlt31PPeNjyb0e0VE2gszmxP0yR4ire8sjnXuiV1ZunEXNXqgvYikGQVB4NzB0SuRXtbpIRFJMwqCwICyfMqLcpmxJHH3KIiIpAIFQcDMOPfErvyjejN1DXqOsYikDwVBjHMHl7HvQCOzVm1NdikiIm1GQRBjdP9SsjJCOj0kImlFQRAjNyvM6P4lvLxMHcYikj4UBE2cc2IZK2v38P6WPckuRUSkTSgImjj3xK4AvLxUp4dEJD0oCJroW5pHv9I8DTchImlDQdCMc04s480VW9hXr8tIRaTjUxA049wTu1LXEGHmyi3JLkVEpNUpCJoxsl8xuZlhnR4SkbSgIGhGTmaYMQNLeGnJJlJxdFYRkWOhIGjBOSd2Ze22fayo1WWkItKxKQhacM6JGo1URNJDQoLAzMaZ2VIzqzazW5tZ/gUzqzWzecHrhphlE81sefCamIh6EqG8qBODuuarn0BEOry4g8DMwsC9wEXAEOAaMxvSzKqPuPuw4HVf8Nli4IfAKGAk8EMzK4q3pkQZO7grb63cytY99ckuRUSk1STiiGAkUO3uK929HpgKjD/Kz34cmO7uW919GzAdGJeAmhLi0uG9aIg4T8+rSXYpIiKtJhFB0AtYEzO/Nmhr6nIze8fMHjezimP8LGY2ycyqzKyqtrZthn84qUcBp/Qq4LGqtW2yPRGRZGirzuJngL7uPpToX/1TjvUL3H2yu1e6e2VZWVnCC2zJlSMqWLR+JwtqdrTZNkVE2lIigqAGqIiZLw/aPuDuW9y9Lpi9DxhxtJ9NtvHDepIVDvH4HB0ViEjHlIggmA0MMrN+ZpYFTACmxa5gZj1iZi8BFgfTzwEXmllR0El8YdDWbhR2yuKCk7vx1LwaPcJSRDqkuIPA3RuAG4n+Al8MPOruC83sNjO7JFjtJjNbaGbzgZuALwSf3QrcTjRMZgO3BW3typUjytm+9wAvLtalpCLS8VgqDqFQWVnpVVVVbba9xogz5o6XOKlHZ3533cg2266ISCKZ2Rx3r2zarjuLj0I4ZHz69F68sqyWjTv3J7scEZGEUhAcpStGlBNx+PPcdtWXLSISNwXBUepflk9lnyIem7NGI5KKSIeiIDgGV1aWs7J2D3NXb092KSIiCaMgOAYXD+1JbmaYx+esOfLKIiIpQkFwDPKzM7jo1O48M3+9nmcsIh2GguAYXTmigt11Dfx94fpklyIikhAKgmM0ql8xFcW5GohORDoMBcExCoWMK0dU8MaKLSzdsCvZ5YiIxE1BcBw+d2Yf8rMzuOeFZckuRUQkbgqC41CUl8UXx/Tlbws2sHCdhqcWkdSmIDhO15/dn4KcDO6erqMCEUltCoLj1CU3ky+d3Z8XFm9i3prtyS5HROS4KQjicN1Z/SjqlMldOioQkRSmIIhDfnYG//KxAby6rJaq99rdYxRERI5KQoLAzMaZ2VIzqzazW5tZfouZLQoeXv+imfWJWdZoZvOC17Smn23vPj+6D6X5WfzP8zoqEJHUFHcQmFkYuBe4CBgCXGNmQ5qs9jZQGTy8/nHgpzHL9rn7sOB1CSmmU1YGXz5nIG+u3MIbKzYnuxwRkWOWiCOCkUC1u69093pgKjA+dgV3n+Hue4PZmUQfUt9hfGZUb7oVZHPX88s0RLWIpJxEBEEvIHY4zrVBW0uuB/4WM59jZlVmNtPMLm3pQ2Y2KVivqra2Nq6CEy0nM8yN5w6k6v1tvLpcRwUiklratLPYzD4LVAI/i2nuEzxD81rgHjMb0Nxn3X2yu1e6e2VZWVkbVHtsrjqjgl6Fudz1/FIdFYhISklEENQAFTHz5UHbIczsfOD7wCXuXnew3d1rgveVwMvA8ATU1OayM8LcdN5A5q/dwZ9mrU52OSIiRy0RQTAbGGRm/cwsC5gAHHL1j5kNB35DNAQ2xbQXmVl2MF0KjAEWJaCmpLhyRAVnDyrl9mcXUb1JA9KJSGqIOwjcvQG4EXgOWAw86u4Lzew2Mzt4FdDPgHzgsSaXiZ4EVJnZfGAGcIe7p2wQhELGf195GrmZYW56eB51DXp4jYi0f5aK57MrKyu9qqoq2WW0aPqijXzp91VM+mh/vveJk5JdjogIAGY2J+iTPYTuLG4FFwzpxmdG9Wbyqyt5XVcRiUg7pyBoJT+4eAgDyvK45dF5bNtTn+xyRERapCBoJblZYX4+YTjb9tbznSfe0SWlItJuKQha0Sm9uvDtjw/m+UUbeXjWmiN/QEQkCRQErez6s/px1sBSbnt2oZ5bICLtkoKglYVCxl1XnUbXzjl89r63mK3hqkWknVEQtIGuBTk8+i+j6VqQzefvn8U/qnUlkYi0HwqCNtK9Sw6PTBpNn5JOXPfgbGYs2XTkD4mItAEFQRsq65zNw186kxO65TPpD1X8fcGGZJckIqIgaGtFeVk8dMOZnNKrC1/901yenveh8flERNqUgiAJuuRm8ofrR1HZp4ivPzKPO/62hH31GpdIRJJDQZAk+dkZPHjdSK4cUc6vX1nBhfe8wivL2tcDd0QkPSgIkig3K8xPrziNh790JpmhEBMfmMVND79N7a66I39YRCRBFATtwOgBJfzt62dz83mD+PuCDZx/1ytMnbWaSETDUohI69Mw1O1M9abdfO/Jd5m1aiu9CnO5fEQ5V44op6K4U7JLE5EU19Iw1AqCdigScf6+cANTZ6/hteW1uMPo/iVcdUY5407uQW5WONklikgKatUgMLNxwM+BMHCfu9/RZHk28HtgBLAFuNrd3wuWfRe4HmgEbnL35460vY4eBLHWbd/HE3PW8tictazeupfO2RmM6l/CsIouDC0v5LTyQrp0ykx2mSKSAlotCMwsDCwDLgDWEn2G8TWxj5w0s68AQ939X81sAnCZu19tZkOAh4GRQE/gBeAEdz/stZTpFAQHRSLOrPe28uTcGma/v5WVtXs+WNavNI/TyrvQtzSPbgU5dO2c/cF7SX424ZAlsXIRaS9aCoKMBHz3SKDa3VcGG5oKjOfQh9CPB34UTD8O/MLMLGif6u51wCozqw6+780E1NWhhELGmf1LOLN/CQA79h1gQc0O5q3Zzvw125m5citPzVv34c9Z9FLV3KwwOZlhcjLC5GSFyckIkZURImRGyIi+h6LThmEJyo6Wvsew2JlDJs3+udQs2hYKJoygxqBmMyMcCuoPXuFQ9N8rI2SEQyHCZmSEjXDQlhkOkRGOvmd+8B4iOyNEdkaYrIzodE5mdD43K0xuZvTfT6EqHVEigqAXEDvY/lpgVEvruHuDme0ASoL2mU0+26u5jZjZJGASQO/evRNQdmrrkpvJmIGljBlY+kFbfUOEzbvr2LhzP5t21UVfO/eza38D+w80su9AY/AeYX99I7vrGog4uDuNEf9gOpKgfqOWvsYPWccPbfd/Lnd3PPgex4lEou0Rd9yj7x/U704k2IfovkT3qSHBV15lZYTolBWmU2aYvOwM8nMyyM+OeeVkUJibRZfcDLp0yqQwN4uC3EwKO2VSmpdNQW4GlqiUFUmQRARBm3D3ycBkiJ4aSnI57VJWRoiehbn0LMxNdinthvs/w6EhEqEh4hxoiL7XB+8HGiPUN0Sob4xQdyBCXUMj9Q0R6hoi7P8gPBvZWx+EaX0ju+sa2VPXwJ76Bnbtb2DDjv3srotO765raLGezLBRnJdFcV42pflZlOZn07Ugm+4FOXQvyKFbl+h7WedsMsO6ulvaRiKCoAaoiJkvD9qaW2etmWUAXYh2Gh/NZ0WOm5kRNgiHjKw2um2moTHCzv0NbN9bz459B9i+7wDb99azZXc9W/dE37fsqWPz7npWbd7Dpp111DdGDvmOkEGPLrmUF+VSUdyJiqJOVBTn0ru4E/1K8yjOy9KRhSRMIoJgNjDIzPoR/SU+Abi2yTrTgIlEz/1fAbzk7m5m04A/mdldRDuLBwGzElCTSNJkhEPBX/1ZR7W+u7N1Tz0bdu5n4879rN+xnw079rN22z7WbN3La8tr2bjz0LvNCztlMqAsnwFlecF7Pid270x5Ua4CQo5Z3EEQnPO/EXiO6OWjD7j7QjO7Dahy92nA/cAfgs7grUTDgmC9R4l2LDcAXz3SFUMiHY2ZUZIfvcLr5J5dml1n/4FGarbvY/XWvays3cOK2t2srN3NjKW1PFq19oP1OudkcFKPAob0KGBIz+j7Cd06k5Wh00zSMt1QJpLiduw7QPWm3SzZsJPF63eyaN1OlmzYxd5gRNvsjBCn9urC8N6FDO9dxPDehfToon6kdKQ7i0XSSCTivL91LwvX7WD+mu28vXo779TsoL4h2hfRvSCHM/oV85EBJXxkQAm9izvplFIaaM37CESknQmFjH6lefQrzeOTQ3sC0cuLF6/fydurtzFn9XbeWrmFZ+ZH7z3pVZjLmf2joXD2oFK6FuQks3xpYzoiEElT7s6K2j28uWIzb67cwpsrtrBt7wEATu3VhfNO6sp5g7txSq8CHS10EDo1JCKHFYk4izfs5OWltby4eCNvr9mOO3QryGbs4K5cOKQ7YwaWquM5hSkIROSYbNldx4yltby0ZCOvLtvM7roGCjtlctEp3fnU0J6M6l+iITdSjIJARI5bXUMjry/fzDPz1/H8oo3srW+krHM2F5/ag0uH9+K08i46fZQCFAQikhD76ht5ackmnpm/jpeWbqK+IcLg7p2ZcEYFlw0v17Do7ZiCQEQSbuf+Azw7fz1TZ6/mnbU7yM4I8YlTe3DNyN6c0bdIRwntjIJARFrVwnU7mDprDU+9XcOuugYGds3n+rP6cdnwXuRk6ql67YGCQETaxL76Rv7y7noefGMVC2p2UpKXxedG9+FzZ/ahJD872eWlNQWBiLQpd2fmyq3c99pKXlyyieyMEJ8+vZwbzu7HgLL8ZJeXlnRnsYi0KTNj9IASRg8ooXrTLu5/fRVPzF3L1NmrueS0ntx83iD6KxDaBR0RiEib2by7jt++tpLfv/E+dQ2NXDa8nJvPG0Tvkk7JLi0t6NSQiLQbtbvq+PUrK/jjzPdpjDhXjCjnxrEDKS9SILQmBYGItDsbd+7nlzOqeXhW9LHn153VlxvPHUjnHN2L0BoUBCLSbq3bvo//fn4pf55bQ2l+Nt/++IlcMaKckIawSKiWgiCu0aPMrNjMppvZ8uC9qJl1hpnZm2a20MzeMbOrY5Y9aGarzGxe8BoWTz0ikpp6FuZy11XDeOqrY+hdnMu3n3iHS+59ndnvbU12aWkh3mEEbwVedPdBwIvBfFN7gc+7+8nAOOAeMyuMWf4tdx8WvObFWY+IpLBhFYU88eWP8PMJw9iyu54rf/0mN/5pLpt27k92aR1avEEwHpgSTE8BLm26grsvc/flwfQ6YBNQFud2RaSDMjPGD+vFi//2MW4+bxDTF23kvLteYeqs1aTiqexUEG8QdHP39cH0BqDb4VY2s5FAFrAipvnHwSmju82sxdsOzWySmVWZWVVtbW2cZYtIe9cpK4NvXHACf//6RxnSo4Bb//wu1/72Ld7bvCfZpXU4R+wsNrMXgO7NLPo+MMXdC2PW3ebuH+onCJb1AF4GJrr7zJi2DUTDYTKwwt1vO1LR6iwWSS+RiDN19hp+8tfFHIhEuOWCE/jimH5khPWQnGNx3HcWu/v5h/nSjWbWw93XB7/UN7WwXgHwF+D7B0Mg+O6DRxN1ZvY74JtHqkdE0k8oZFw7qjdjB3flB08t4L/+uoRn5q/nrqtOY1C3zskuL+XFG6fTgInB9ETg6aYrmFkW8CTwe3d/vMmyHsG7Ee1fWBBnPSLSgXXvksNvPz+Ce689nZrt+/jUL17nT2+p7yBe8QbBHcAFZrYcOD+Yx8wqzey+YJ2rgI8CX2jmMtGHzOxd4F2gFPjPOOsRkQ7OzLh4aA/+fvPZnNG3mO89+S5f/uNctu+tT3ZpKUs3lIlIyopEnPtfX8VPn1tCSV42d189jNEDSpJdVrvVKjeUiYgkUyhkfOmj/XnyK2PolBXm2vtm8t/PLeVAYyTZpaUUBYGIpLxTenXhma+dxVUjKvjFjGo+d/9bbNldl+yyUoaCQEQ6hLzsDO68Yih3X30ab6/eziW/+AcLanYku6yUoCAQkQ7lsuHlPP6vH8HdufxXb/DU2zXJLqndUxCISIdzankXpn3tLIZVFPL1R+Zx+7OLaFC/QYsUBCLSIZXmZ/PHG0bxhY/05f7XV/H5B2axbY8uMW2OgkBEOqzMcIgfXXIyP7tiKFXvb+PyX73Bmq17k11Wu6MgEJEO78rKCh66YRRb9tRz2S/fUCdyEwoCEUkLZ/Qt5okvjyY7I8TVv3mTV5dpFOODFAQikjYGdu3Mn7/yESqKO/HFB2fz57lrk11Su6AgEJG00q0gh0f/dTSj+hdzy6PzuXdGddoPWqcgEJG0U5CTye++MJJLh/XkZ88t5fZnF6d1GBzxeQQiIh1RVkaIu64aRlFeFg/8YxUNkQj/ccnJREfFTy8KAhFJW6GQ8e+fHEJmOMTkV1fSGHFuH38KoVB6hYGCQETSmpnx3YsGkxEyfvnyChojzn9ddmpahUFcQWBmxcAjQF/gPeAqd9/WzHqNRB8+A7Da3S8J2vsBU4ESYA7wOXfXrX8i0qbMjG99/EQyQsb/vlRNQ8S58/KhhNMkDOLtLL4VeNHdBwEvBvPN2efuw4LXJTHtdwJ3u/tAYBtwfZz1iIgcFzPjlgtP5Bvnn8Djc9byzcfm0xhJjw7keINgPDAlmJ5C9LnDRyV4TvFY4OBzjI/p8yIireHm8wfxzQtP4Mm3a/jWY/OJpEEYxNtH0M3d1wfTG4BuLayXY2ZVQANwh7s/RfR00HZ3bwjWWQv0amlDZjYJmATQu3fvOMsWEWnZjWMH4Q7/M30ZBbmZ/PBTQzr01URHDAIzewHo3syi78fOuLubWUvR2cfda8ysP/BS8MD6Yxrsw90nA5Mh+sziY/msiMixunHsQLbvO8D9r6+iOC+Lm84blOySWs0Rg8Ddz29pmZltNLMe7r7ezHoAm1r4jprgfaWZvQwMB54ACs0sIzgqKAf0BAkRaRfMjO9/4iS27a3nrunLKOqUyedG9012Wa0i3j6CacDEYHoi8HTTFcysyMyyg+lSYAywyKO38c0Arjjc50VEkiUUMu68fCjnn9SVf5+2kGnz1yW7pFYRbxDcAVxgZsuB84N5zKzSzO4L1jkJqDKz+UR/8d/h7ouCZd8BbjGzaqJ9BvfHWY+ISEJlhkP84trTOaNvMbc8Mo+XlzZ74iOlWSqOr1FZWelVVVXJLkNE0sjO/QeY8JuZrNy8m4duGMWIPsXJLumYmdkcd69s2q5B50REjkJBTiZTvjiS7gU53DClive37El2SQmjIBAROUplnbP53XUjiThcP6WKnfsPJLukhFAQiIgcg36lefz6syN4b/MevvrQXBoaI8kuKW4KAhGRYzR6QAn/eekpvLZ8M7c/u+jIH2jnNPqoiMhxmDCyN9WbdnPf66sY2DU/pe8xUBCIiByn737iJFZt3sOPnllEn5I8PnpCWbJLOi46NSQicpzCIePn1wxnUNd8vvrQXKo37Up2ScdFQSAiEof87Azum1hJdmaIG1L0SiIFgYhInMqLOvGrz45gzbZ9fPuxd0i1G3UVBCIiCXBG32JuHTeYvy/cwP2vr0p2OcdEQSAikiA3nN2Pj5/cjZ/8bQmz39ua7HKOmoJARCRBzIyfXXkaFUW53PinuWzeXZfsko6KgkBEJIEKcjL55WdGsH3vAW56+O2UeO6xgkBEJMGG9Czg9ktP4Y0VW7h7+rJkl3NECgIRkVZwVWUFV1WW84sZ1by0ZGOyyzksBYGISCu5bfwpDOlRwC2Pzmfjzv3JLqdFcQWBmRWb2XQzWx68FzWzzrlmNi/mtd/MLg2WPWhmq2KWDYunHhGR9iQnM8z/XTuc/Qca+eZj84m00/6CeI8IbgVedPdBwIvB/CHcfYa7D3P3YcBYYC/wfMwq3zq43N3nxVmPiEi7MqAsnx9cPITXlm9mypvvJbucZsUbBOOBKcH0FODSI6x/BfA3d98b53ZFRFLGZ0b1Zuzgrvzkb0tYtrH9jUcUbxB0c/f1wfQGoNsR1p8APNyk7cdm9o6Z3W1m2S190MwmmVmVmVXV1tbGUbKISNsyM+68fCidszP4+tR51DU0JrukQxwxCMzsBTNb0MxrfOx6Hh1co8UTYGbWAzgVeC6m+bvAYOAMoBj4Tkufd/fJ7l7p7pVlZak51KuIpK+yztnceflQFq3fyV3t7JLSIz6PwN3Pb2mZmW00sx7uvj74Rb/pMF91FfCku38wNF/M0USdmf0O+OZR1i0iknLOH9KNa0b2ZvKrKznnhK6MHlCS7JKA+E8NTQMmBtMTgacPs+41NDktFIQHZmZE+xcWxFmPiEi79v8+eRJ9S/L4t0fnsWNf+xiyOt4guAO4wMyWA+cH85hZpZndd3AlM+sLVACvNPn8Q2b2LvAuUAr8Z5z1iIi0a52yMrj76mFs3FXHvz/dPv72jetRle6+BTivmfYq4IaY+feAXs2sNzae7YuIpKJhFYXcNHYQd7+wjItP7cGFJ3dPaj26s1hEJAm+cu4ABnfvzP97ekHSn2qmIBARSYLMcIifXjGU2l11/OSvS5Jai4JARCRJhpYXcv1Z/Xh41mpmrtyStDoUBCIiSXTLBSfSu7gTtz7xDvsPJOdGMwWBiEgS5WaF+cmnT+W9LXu554XlSalBQSAikmRjBpZydWUFv31tJQtqdrT59hUEIiLtwPc+cRLFeVl854l3aGiMtOm2FQQiIu1Al06Z3D7+ZBau28lvX1vVpttWEIiItBPjTunBuJO7c88Ly1i9pe1G61cQiIi0I/8x/mQyQsZtzy5qs20qCERE2pFuBTl87bxBvLB4IzOWHm5A58RREIiItDNfHNOP/qV53PbMojZ5iI2CQESkncnKCPGjS05m1eY93P9663ccKwhERNqhj55QxsdP7sb/vVjN+h37WnVbCgIRkXbqBxcPIeLOj/+yuFW3oyAQEWmnKoo78ZVzBvLsO+t5Y8XmVttOXEFgZlea2UIzi5hZ5WHWG2dmS82s2sxujWnvZ2ZvBe2PmFlWPPWIiHQ0//Kx/pQX5fKjaQs50Ep3HMd7RLAA+DTwaksrmFkYuBe4CBgCXGNmQ4LFdwJ3u/tAYBtwfZz1iIh0KDmZYf79k0NYtnE3f3jz/VbZRlxB4O6L3X3pEVYbCVS7+0p3rwemAuODB9aPBR4P1ptC9AH2IiIS44Ih3fjYCWXcPX0ZtbvqEv79bdFH0AtYEzO/NmgrAba7e0OT9maZ2SQzqzKzqtra2lYrVkSkvTEzfvipIYzoW9Qqzyw44sPrzewFoLknK3/f3Z9OeEUtcPfJwGSAyspKb6vtioi0B/3L8nnwupGt8t1HDAJ3Pz/ObdQAFTHz5UHbFqDQzDKCo4KD7SIi0oba4tTQbGBQcIVQFjABmObuDswArgjWmwi02RGGiIhExXv56GVmthYYDfzFzJ4L2nua2V8Bgr/2bwSeAxYDj7r7wuArvgPcYmbVRPsM7o+nHhEROXYW/cM8tVRWVnpVVVWyyxARSSlmNsfdP3TPl+4sFhFJcwoCEZE0pyAQEUlzCgIRkTSXkp3FZlYLHO+gG6VA6w3j1z5pn9OD9rnji3d/+7h7WdPGlAyCeJhZVXO95h2Z9jk9aJ87vtbaX50aEhFJcwoCEZE0l45BMDnZBSSB9jk9aJ87vlbZ37TrIxARkUOl4xGBiIjEUBCIiKS5tAoCMxtnZkvNrNrMbk12Pa3BzB4ws01mtiCmrdjMppvZ8uC9KJk1JpKZVZjZDDNbZGYLzezmoL0j73OOmc0ys/nBPv9H0N7PzN4Kfr4fCYZ971DMLGxmb5vZs8F8h95nM3vPzN41s3lmVhW0JfxnO22CwMzCwL3ARcAQ4BozG5LcqlrFg8C4Jm23Ai+6+yDgxWC+o2gA/s3dhwBnAl8N/rt25H2uA8a6+2nAMGCcmZ0J3Anc7e4DgW3A9ckrsdXcTHQ4+4PSYZ/PdfdhMfcPJPxnO22CABgJVLv7SnevB6YC45NcU8K5+6vA1ibN44EpwfQU4NK2rKk1uft6d58bTO8i+kuiFx17n93ddwezmcHLgbHA40F7h9pnADMrBy4G7gvmjQ6+zy1I+M92OgVBL2BNzPzaoC0ddHP39cH0BqBbMotpLWbWFxgOvEUH3+fgFMk8YBMwHVgBbA8eBAUd8+f7HuDbQCSYL6Hj77MDz5vZHDObFLQl/Gf7iM8slo7F3d3MOtw1w2aWDzwBfN3dd0b/WIzqiPvs7o3AMDMrBJ4EBie3otZlZp8ENrn7HDM7J8nltKWz3L3GzLoC081sSezCRP1sp9MRQQ1QETNfHrSlg41m1gMgeN+U5HoSyswyiYbAQ+7+56C5Q+/zQe6+neizv0cDhWZ28I+7jvbzPQa4xMzeI3padyzwczr2PuPuNcH7JqKBP5JW+NlOpyCYDQwKrjLIAiYA05JcU1uZBkwMpicCTyexloQKzhPfDyx297tiFnXkfS4LjgQws1zgAqJ9IzOAK4LVOtQ+u/t33b3c3fsS/X/3JXf/DB14n80sz8w6H5wGLgQW0Ao/22l1Z7GZfYLoecYw8IC7/zi5FSWemT0MnEN0uNqNwA+Bp4BHgd5Eh+++yt2bdiinJDM7C3gNeJd/njv+HtF+go66z0OJdhKGif4x96i732Zm/Yn+tVwMvA181t3rkldp6whODX3T3T/Zkfc52Lcng9kM4E/u/mMzKyHBP9tpFQQiIvJh6XRqSEREmqEgEBFJcwoCEZE0pyAQEUlzCgIRkTSnIBARSXMKAhGRNPf/AbvWB9Y65eBmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(weights_condense)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#plt.imshow(model1_Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "214.7166302201803\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6eba614b70>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAAD4CAYAAABLyM+7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABKnElEQVR4nO3dd3hUVf7H8fdJI0ACJJAESCEhJPQeOgquomLDir2wKurKWrfolp9r27WtZVeERcWOiAVFRRGVoiC9d0ISSAIkgVBSSJ3z+yODm4WAlCR3JvN5PQ/PzJy5mfvNJTN3Pveee46x1iIiIiIiIiLeyc/pAkREREREROTUKdSJiIiIiIh4MYU6ERERERERL6ZQJyIiIiIi4sUU6kRERERERLxYgNMFnIhWrVrZ+Ph4p8sQEZF6sHz58j3W2gin6/AW2keKiPiG4+0fvSLUxcfHs2zZMqfLEBGRemCM2e50Dd5E+0gREd9wvP2jul+KiIiIiIh4MYU6ERERERERL6ZQJyIiIiIi4sUU6kRERERERLyYQp2IiIiIiIgXU6gTERERERHxYgp1IiIiIiIiXswr5qkTETms0mXZX1xGflEZe4v+e7uvqAyXtQT6+xHobwjwq7qteuxH6+bBxIY1oU2LYAL9dTxLRERE6sfy7ftYlpHPHcMS62wdCnUi4pFKyitJzS1k8+4CNucUsGl3AVt2F5BbUILLnvrr+vsZ2rYIJi68CXHhTYgNb0Lv2DB6x7UgONC/9n4BERER8WnWWt7+aTtPfLmBti0ac8PAdjRtVDfxS6FORDxC9v5DLEjdw0/b9rImaz8Ze4updKe3oAA/OkSEMDixJTFhjQlvGkRY0yBaNm1EeNMgWoYEEdYkCH8/Q3mliwqXpaLSRVmli4pKS1mFi10HSsjML2aH+1/mvmJmb8hhT2EZAI0C/OgTF8agxJYMSmxJz5gWBAXojJ6IiIicvOKyCv70yVo+XbWTX3WK5IXRveos0IFCnYg4JL+ojIXb9rBw214Wpu4hY28xAK1CgugdF8YF3dvQsXUonVqHEt+yKQEn2GXS36/ms23xrZoyKLHlUe0HistZkpHPT9v28lPaXp6fvQVmQ+NAf1Liw7ikZ1su7NGGJkH6uBQREZFflrGniDvfXc7mnAIeHJHM3Wd1wM/P1Ok69S1FROpN9v5DzFyziy/W7mJ15n4AQhoFMLB9ODcNimdIh1YkR4VgTN1+8FXXvEkgI7pEMaJLFAD7ispYnJ7PorS9zN2cy+8/WsPfZqzn4p5tGd0vlt6xLeq1PhEREfEeszfk8MC0Vfj7Gd4c059hyRH1sl6FOhGpUzkHS/hyzS6+WLOTFTv2A9A9ujkPjkhmSFIrekQ3P+GzcPUhrGkQ53drzfndWmNtF5Zt38cHSzP5bNVOpi7NJCkyhNEpsVzWJ5pWIY2cLldEREQ8QKXL8sLsLbw8J5Vu0c2YcH1fYsOb1Nv6jbWnMeJAPUlJSbHLli1zugwROUHFZRXMWLWTT1Zks3R7PtZC5zbNuKhHGy7s3ob4Vk2dLvGkFZZW8OWanXywNJMVO/YT6G+4tn8c437VgcjQYKfLa1CMMcuttSlO1+EttI8UEXFWXkEp932wkgWpe7k6JZZHR3Wtk8HXjrd/1Jk6Eak1m3cX8N7i7UxfkU1BaQVJkSHcf04yF/ZoQ2JEiNPlnZaQRgFc3S+Oq/vFkZpbwOs/ZvDe4h18uCyLMUPiuePMRJo3CXS6TBEREalHi9P28tv3V3LgUDlPX9Gdq/vFOVKHQp2InJbSikq+XrebdxdtZ2nGPoIC/LiwextuGBhHn7iwBnn9WYfIUP5xeXfuOLM9L3y7hQnztvHuou3cMSyRMUPiNaiKiIhIA+dyWSbO38ZzszbTrmVT3vp1fzq3aeZYPfrmISKnZG9hKW8urDpblV9URruWTfjTBZ24sm8s4U2DnC6vXsS3aspL1/TmzmGJPDdrM8/O2swbCzK45+wOXD+gHf51PNKViIiI1L99RWU8MG0VczbncVGPNvzj8u6EBjvbW0ehTkROSvb+Q7w6P42pS3dQUu5iRJcobhrUjiGJrep8uF5P1blNM16/pR/Lt+fzzNeb+b/P1vPZqp08d1VPErzw+kERERGp2Yod+xj33gr2FJbx+Kiu3DCwnUf0SlKoE5ETsjWngAnztjFj1U4ALu0dzZ3D2tMhMtThyjxH33bhTB07kM9W7eT/PlvHyJfm89D5nbhpULzPBl4REZGGwOWyTF6QzlNfbaJNi2A+vmsw3WOaO13WzxTqROS41mUf4F/fbeWbDTkEB/pxw8B23H5me6JbNHa6NI9kjOHS3tEMbN+Shz5Zw98+38Cs9Tk8c2WPeh3aWOqPMSYWeBuIAiwwyVr7kjHmb8DtQJ570T9Za2e6f+Zh4FagErjHWjvL3X4+8BLgD7xmrX2qPn8XERE5Wn5RGQ+6u1ue1zWKZ67o6XGDoynUiUiN0vIK+efsLXy5ZhfNggO451cduGVIgs9cL3e6WjcP5o1b+vHB0kwe/2ID5784n79e1IWr+8V6RDcNqVUVwIPW2hXGmFBguTFmtvu5F6y1z1Vf2BjTBbgG6Aq0Bb41xiS7nx4PjACygKXGmBnW2g318luIiMhRFqXt5d6pK9lXVM5jo7pyo4d0tzySQp2I/I9dBw7x0rdb+XB5Fo0C/Pjtrzpw+5ntaebwBcDeyBjDNf3jGNKhFX/4aA0PfbKWWet38+LVvT3uCJ+cOmvtLmCX+36BMWYjEH2cHxkFTLXWlgLpxphUoL/7uVRrbRqAMWaqe1mFOhGRelbpsvzru638+/utxLdsyuRb+tG1red0tzySX228iDFmsjEm1xizrlpbuDFmtjFmq/s2zN1ujDH/MsakGmPWGGP61EYNInJ68ovKeOKLDQx7di6frMjmxoHtmPf7s3jw3I4KdKcpNrwJ7902gL9d3IUfU/dw2YQFpO8pcrosqQPGmHigN7DY3TTOva+bfHg/SFXgy6z2Y1nutmO117SescaYZcaYZXl5eTUtIiIip2j3gRKue3URL323lUt7R/P5b4d6dKCDWgp1wJvA+Ue0PQR8Z61NAr5zPwYYCSS5/40FJtRSDSJyCkrKKxk/J5Uzn5nD5AXpXNKzLd//bhh/u6QrEaGNnC6vwfDzM9wyJIH3bhvIvqIyLh2/gJ+27XW6LKlFxpgQ4GPgPmvtQar2b4lAL6rO5P2zttZlrZ1krU2x1qZERETU1suKiPi87zbmcMG/fmBt9gH+eVVPnh/di6aNPL9zY62EOmvtfCD/iOZRwFvu+28Bl1Zrf9tWWQS0MMa0qY06ROTEWWv5et1uRrwwj2dnbWZg+5bMuu9MnruqJzFhGtCjrvRPCOfTu4cQEdqIG19fzAdLdzhdktQCY0wgVYHuPWvtJwDW2hxrbaW11gW8yn+7WGYDsdV+PMbddqx2ERGpYyXllfzfZ+u49a1ltG4WzOe/HcoVfWOcLuuE1WXsjHJfZwCwm6pRweDY3Ut2VWvDGDOWqjN5xMXF1WGZIr5n0+6DPPb5BhZu20tyVAjv3jqAoUmtnC7LZ7Rr2ZRPfjOYu99bwR8/XktqbiEPjeysycq9lKm6Yv51YKO19vlq7W2q7QcvAw5fojADmGKMeZ6qgVKSgCWAAZKMMQlUhblrgOvq57cQEfFdG3cd5N6pK9mSU8jtZyTwu/M60ijA3+myTkq9nEu01lpjjD3Jn5kETAJISUk5qZ8VkZrlF5Xx/OzNTFm8g2aNA3lsVFeu6x9HgH9t9cSWE9UsOJA3bunH419s4NUf0knfU8SL1/QmxAu6eMhRhgA3AmuNMavcbX8CrjXG9KJqmoMM4A4Aa+16Y8w0qgZAqQDuttZWAhhjxgGzqJrSYLK1dn39/RoiIr7FWsubCzP4x1ebaN44kHdu7c8ZSd7Zpb0uvz3kHD5K6e5emetuV/cSkXrmclneW7KDZ7/eRFFZJTcNiue+c5Jo0UTTEzgpwN+PR0d1IzEyhEc/38DoiT/x7m0DNG2El7HW/kjVWbYjzTzOzzwJPFlD+8zj/ZyIiNSOvIJSfv/RauZuzuOczpE8fUUPWoZ471gCdRnqZgA3A0+5bz+r1j7OPVTzAOBAte4pIlLLNuw8yJ+mr2VV5n4GtW/Jo6O6khwV6nRZUs1Ng+KJC2/CHe8s54bXFjPl9gEK3CIiInXku405/PHjNRSUVPD4qK7c4KFzz52MWgl1xpj3geFAK2NMFvAIVWFumjHmVmA7MNq9+EzgAiAVKAbG1EYNIvK/ikorePHbLUxekEGLxoE8P7onl/WO9voPrYZqeMdIJt2Uwu1vLeOG1xfz3q0DNZediIhILSoqreCJLzfy/pIddGodypTbBzaYA921Euqstdce46mza1jWAnfXxnpFpGbfbsjhkRnryd5/iGv6xfLQyE468+MFhiVH8J8b+zL2nWXcNHkx79w2QHMEioiI1IIVO/bxwAer2J5fzB3D2vPAiGSvGwzleDQ6gkgDknOwhDveWcZtby8jpFEAH905iKeu6KFA50XO6hTJhOv7smHXQW6evISCknKnSxIREfFa5ZUunp+9hasm/kR5peX92wfy8MjODSrQQT2Nfikidctay4fLs3j8iw2UVbj44/mduO2MBAI1qqVXOqdLFC9f14e731vBLW8s5a1f99eomCIiIidpW14h93+wijVZB7i8TzR/u6Rrg+0Bo28JIl4ue/8hHv5kLfO35NE/IZynr+hBQqumTpclp+m8rq3597W9Gff+Ssa8sYQ3x/SnqYKdiIjIL3K5LG//lMFTX28iONCfV67vwwXd2zhdVp3SNwQRL2WtZcqSHfxj5iZc1vLYqK7cMKAdfprAusEY2b0NL1nLPe+v5Pa3l/HmmP4EBejsq4iIyLFk7SvmDx+tYeG2vQxLjuCZK3sQ1SzY6bLqnEKdiBfasbeYhz6p+sAa0qElT13eg9jwJk6XJXXgoh5tKS138eCHq3lkxjr+fll3jWAqIiJyBGstHy7L4rEvNmCt5anLu3N1v1if2Wcq1Il4kcNn5578ciN+xvCPy7tzjQ99YPmqK/rGsC2vkFfmbiMpMpRfD01wuiQRERGPkXuwhIc+Wcv3m3IZ2D6cZ6/s6XMHuxXqRLxE7sES/vDxGuZuzuOMpFY8fUUP2rZo7HRZUk9+d25HUnMLeeLLDbSPaMrwjpFOlyQiIuIoay2fr9nFXz9dR2lFJY9c3IWbB8X75KUoujhDxAt8uWYX5744n0Vpe3n0kq68Naa/Ap2P8fMzvHB1Lzq2bsZvp6wkNbfA6ZJEREQck1dQyl3vruCe91fSPqIpM+85gzFDEnwy0IFCnYhHO3ConPumruTuKStoF96EL+85g5sH++YRKIGmjQJ47eYUGgX6cetby9hXVOZ0SSIiIvXKWstnq7IZ8cI8vt+cy8MjO/HhHYNoHxHidGmOUqgT8VALU/dw/ovz+XzNLu47J4mP7hpMoo9/YAlEt2jMf25MYdf+Eu56bznllS6nSxIREakXuQdLuP3t5dw7dRUJrarOzt0xLJEAzcurUCfiacoqXPx95kaue20xjYP8+eSuwdx3TrImEpef9W0XxlNXdGdRWj6PzFiPtdbpkkREROqMtZZPVmRxzvPz+GFrHn+5sDMf3TmYDpE62H2YBkoR8SBpeYXcM3Ul67IPcv2AOP5yYRcaB/k7XZZ4oMv7xLAlp5CJ87bRMSqUmwfHO12SiIhIrdt14BB/nr6O7zflktIujGeu7OHzXS1rolAn4gEOz63yyIz1NAr04z839uW8rq2dLks83B/O60hqbgFPfLmBvu3C6Bbd3OmSREREaoXLVTWN01NfbaLSZfnrRV24ZXA8/hpXoEbqzyXisAPF5YybspI/fLyGXrEt+PreMxXo5IT4+RmevbInLZs24p6pKykuq3C6JBERkdOWllfINa8u4i+frqNnbHNm3Xcmtw5NUKA7DoU6EQctSc9n5EvzmbV+N388vxPv3jaA1s2DnS5LvEhY0yCeH92T9D1FPP7FRqfLEREROWUVlS4mztvGyJd+YNOugzxzZQ/evXUAcS19ayLxU6HulyIOqHRZXpmTygvfbiE2vAkf3zWYnrEtnC5LvNTgDq0Ye2Z7/jMvjeEdI3SmV0REvM76nQf448drWJd9kPO7tuaxUV2JbKYD3SdKoU6knuUWlHD/B6tYkLqXUb3a8uRl3QlppLeinJ4HR3RkQeoeHnJ3443SjlBERLzAobJKXvxuC6/9kE5YkyAmXN+Hkd3bOF2W11H3S5F69OPWPVzw0g8s376PZ67owYtX91Kgk1oRFODHS9f0pqTcxQPTVuFyaZoDERHxbPO35HHui/P4z7w0ruobw7cPnKlAd4oU6kTqQUWli+dmbebGyYsJaxLEjHFDGd0vFmN0wa/UnsSIEP7v4i4sSN3Laz+mOV2OiIhIjfYWlnL/B6u4afISAv39mDp2IE9d0YMWTYKcLs1r6RSBSB3bdeAQ97y/kqUZ+xidEsOjl3TT3HNSZ67pF8vczbk8O2szgxNbaZoDERHxGNZaPl6RzRNfbqCotIJ7zk7iN8MTCQ7U96LTpTN1InVo3pY8LnjpB9bvPMgLV/fkmSt7KtBJnTLG8NTlPTTNgYiIeJRteYVc/9pifvfhajpEhDDznjN4YESyAl0tUagTqQOVLsvz32zmljeWENUsmM9/O5TLesc4XZb4iOrTHPx9pqY5EBER55SUV/L87C2MfPEH1mYf4MnLujHtjkEkRYU6XVqDou6XIrVsT2Ep905dyYLUvVzVN4bHRqm7pdS/wR1aMWZwApMXpHNZ7xj6tgtzuiQREfExP2zN46+friNjbzGX9mrLny/sQkRoI6fLapAU6kRq0ZL0fMZNWcGBQ+U8c0UPRveLdbok8WEPnJvMV+t28efpa/n8t0MJ9FfnDBERqXu5BSU88cVGZqzeSUKrprx76wCGJrVyuqwGTXt4kVpgrWXivG1c++oimgT5M/03QxToxHEhjQJ45OKubNpdwOQf050up0EzxsQaY+YYYzYYY9YbY+51t4cbY2YbY7a6b8Pc7cYY8y9jTKoxZo0xpk+117rZvfxWY8zNTv1OIiInq9JleWfRds7+5zy+Xrebe89O4qt7z1Cgqwc6Uydymg4cKud3H65m9oYcLujemqev6EFocKDTZYkAcF7XKM7pHMmL327lwh5tiAlr4nRJDVUF8KC1doUxJhRYboyZDdwCfGetfcoY8xDwEPBHYCSQ5P43AJgADDDGhAOPACmAdb/ODGvtvnr/jURETsLqzP389bN1rMk6wODEljxxaTfaR4Q4XZbP0Jk6kdOwcddBLnn5R+ZsyuWvF3Vh/HV9FOjEoxhj+NslXQF45LP1WKtJyeuCtXaXtXaF+34BsBGIBkYBb7kXewu41H1/FPC2rbIIaGGMaQOcB8y21ua7g9xs4Pz6+01ERE7O/uIy/jR9LZe+soDdB0r417W9ee+2AQp09Uxn6kRO0ScrsvjT9LU0Cw5k6tiBpMSHO12SSI1iwppw/4gk/j5zE7PW53B+t9ZOl9SgGWPigd7AYiDKWrvL/dRuIMp9PxrIrPZjWe62Y7WLiHgUl8vy0fIsnvp6EwcOlTNmcAL3j0jSwW2HKNSJnKTSikoe/2ID7y7awYCEcP59XW8iQ4OdLkvkuMYMSeCTFdn8bcZ6hia1IqSRPv7rgjEmBPgYuM9ae9AY8/Nz1lprjKmVU6XGmLHAWIC4uLjaeEkRkRO2YedB/vrZOpZv30dKuzAev7Qbnds0c7osn6bulyInYef+Q1z9n0W8u2gHY89sz3u3DVCgE68Q6O/H3y/vTk5BCc9/s8XpchokY0wgVYHuPWvtJ+7mHHe3Sty3ue72bKD6aEox7rZjtf8Pa+0ka22KtTYlIiKidn8REZFjOFBcziOfreOif/9Axp4inr2yB9PuGKRA5wF0qFbkBC1I3cNv319JaXklr1zfhwu6t3G6JJGT0icujOv6x/HmwnQu7xNNt+jmTpfUYJiqU3KvAxuttc9Xe2oGcDPwlPv2s2rt44wxU6kaKOWAtXaXMWYW8PfDo2QC5wIP18fvICJyLC6X5cPlmTz99Wb2F5dx48B2PDCiI82bqKulp1CoE/kF1lomzU/j6a830T4ihIk39KVDpC7+Fe/0h/M7MWt9Dn+avpbpvxmCv5/55R+SEzEEuBFYa4xZ5W77E1Vhbpox5lZgOzDa/dxM4AIgFSgGxgBYa/ONMY8DS93LPWatza+X30BEpAarM/fzfzPWszpzP/3iw3j0kgF0aaszc56mzkOdMSYDKAAqgQprbYp7yOYPgHggAxit4ZrFExWVVvCHj9bw5dpdXNC9Nc9c2VPXIolXa944kL9e1Jl7p67ivcXbuWlQvNMlNQjW2h+BYyXks2tY3gJ3H+O1JgOTa686EZGTl19UxrOzNjF1aSatQhrxwtU9ubRXNNWvFRbPUV/fTs+y1u6p9vghap63R8RjpOUVcsc7y9mWV8jDIzsx9sz2+iCTBuGSnm35YGkmL8zewqW9o2mmkcpERMStvNLFu4u288LsLRSXVXLb0ATuOVujWno6pwZKOda8PSIeYfaGHEa9vIA9haW8/esB3DEsUYFOGgxjDH+6oDP7isuZOHeb0+WIiIiH+HHrHi546Qce/XwDPWNb8NW9Z/DnC7so0HmB+jhTZ4Fv3MM4/8daO4ljz9vzMw3XLE5wuSwvfreVf323lW7RzZh4Q19iwpo4XZZIresW3ZxRvdry+o/p3DioHW2aN3a6JBERcUhmfjFPfLmBWetziA1vzKQb+zKiS5QOaHuR+gh1Q6212caYSGC2MWZT9SePNW+PO/xNAkhJSamVeX1EjufAoXLu/2AV32/K5cq+MTxxaTeCA/2dLkukzvzu3I58tXY3L8zewjNX9nS6HBERqWfFZRVMmLuN/8xPw98Yfn9eR24dmqDvP16ozkOdtTbbfZtrjJkO9Mc9b497+Obq8/aIOGJrTgFj31lOZn4xj4/qyg0D2+nolDR4seFNuHFQO95YkM6tQ9vTsXWo0yWJiEg9cLksn63O5umvNrP7YAmjerXloZGd1GvDi9XpNXXGmKbGmNDD96mab2cd/523B/533h6Revf1ul1cOn4BBSUVvD92IDcOilegE58x7qwONG0UwNNfb/rlhUVExOut2LGPyycs5P4PVhMR2oiP7hzES9f0VqDzcnV9pi4KmO7+ghwATLHWfm2MWUrN8/aI1JtKl+WF2Vt4eU4qPWNbMPGGPvpAE58T1jSI3wzvwNNfb+KnbXsZlNjS6ZJERKQO7DpwiKe/2sSnq3YSEdqIZ6/swRV9YvDTfKUNQp2GOmttGnDUhRrW2r3UMG+PSH05cKic+6auZM7mPEanxPDYKF0/J75rzJB43v4pg398tZFPfzNEO3gRkQbkUFklk+anMXHeNiqt5e6zErlreAfNu9vA6H9TfI6unxP5X8GB/jwwIpnff7SGL9fu4uKebZ0uSURETtPh6+ae+Xozuw6UcEH31jw8sjOx4RrVuyFSqBOfMmv9bh74YBWNg/yZcvtA+ieEO12SiEe4vE8Mr/+YzrOzNnNe19YEBTg1jamIiJyupRn5PPHFBlZnHaB7dHNevLoXA9qre31DplAnPsHlsrz03VZe+m4rPWOaM/HGvrp+TqQafz/DH0d2YswbS3lv8XbGDElwuiQRETlJO/YW8/TXm/hy7S5aNwvm+dE9ubRXtLrV+wCFOmnwCksreOCDVXyzIYfL+0Tz98u66/o5kRoMT45gcGJL/v19Klf0jaFZcKDTJYmIyAk4WFLO+O9TeWNBBv5+hvvPSeb2MxNoEqSv+r5C/9PSoGXsKeL2t5eRtqeI/7uoC2OGaLoCkWMxxvDwyM5c/PKPvDY/jQfO7eh0SSIichzllS6mLN7Bi99uYf+hcq7oE8Pvz+tIVLNgp0uTeqZQJw3W3M253PP+Svz9DO/8uj+DO7RyuiQRj9c9pjnnd23NGwszuO3M9jpbJyLigay1zN6Qw1NfbSJtTxGD2rfkzxd2plt0c6dLE4co1EmDY61l0vw0nv56E8lRobx6U4pGehI5CXef1YGv1+/mnZ+2c/dZHZwuR0REqlmTtZ8nvtzIkvR8OkSGMPmWFM7qGKmeSD5OoU4alENllfzx4zXMWL2TC7u34dmreqg/uchJ6h7TnOEdI3j9x3TGDInXe0hExANk7SvmuVmb+XTVTlo2DeKJS7txTb9YAvw1WrEo1EkDkr3/EHe8s4z1Ow/y+/M68pvhiTpqJXKKxp3VgSsn/sSUxTu47Yz2TpcjIuKzDhSXM35uKm8uyMAYuPusRO4clkiousdLNQp10iAsSc/nN+8tp6TcxWs3pXB25yinSxLxainx4QxsH86rP6Rx46B2NArQiLEiIvWptKKStxdu5+U5qRwsqRoE5cFzkzUlk9RIoU683nuLt/PIZ+uJDW/C1LF96RAZ6nRJIg3CuLOSuOH1xXy0PIvrB7RzuhwREZ/gclk+X7OTZ2dtJmvfIYYlR/DQyE50btPM6dLEgynUidcqq3Dx6OfreW/xDoYlR/Cva3vTvLG6IojUliEdWtIrtgUT5m5jdEosgbpuQ0SkTi1M3cM/vtrE2uwDdGnTjHdv7cHQJI3eLb9MoU680p7CUn7z7gqWZORz57BEfn9eR/z9dP2cSG0yxjDurA7c9vYyZqzayRV9Y5wuSUSkQVq/8wBPf72Z+VvyiG7RmOdH9+TSXtH46buNnCCFOvE667IPcMc7y9lTWMpL1/RiVK9op0sSabDO7hxJ5zbNGD83lUt7R+vgiYhILcrML+b52Vv4dFU2zYID+fMFnblxUDuCA3Uds5wchTrxKl+s2cnvPlxNWJMgPrpzMN1jNMmmSF0yxnD3WYmMm7KSr9ft5sIebZwuSUTE6+0rKuPlOam889N2jIE7zkzkruGJuoxETplCnXgFl8vy/OwtvDwnlb7twphwQx8iQ4OdLkvEJ4zs1ob2EVXvvwu6t9ZUISIip6i4rII3FmQwcd42ikoruLJvDPeP0IiWcvoU6sTjFZSUc/8Hq/h2Yy7X9Ivl0VFdNby6SD3y9zPcPbwDD364mu835WrKEBGRk1Re6WLq0kz+9d1W8gpKOadzFH84vyPJURqxW2qHQp14tIw9Rdz+9jLS9hTx2Kiu3Diwnc4SiDjgkl5teeHbLfz7+1R+1SlS70MRkRPgclm+XLuLf36zmYy9xfSLD2PC9X1IiQ93ujRpYBTqxGP9sDWPcVNWYgy8c2t/BidqSF8RpwT6+3HX8ET+PH0dC7ftZUgHvR9FRI7FWsuPqXt4+utNrMs+SKfWoUy+JYWzOuqgmNQNhTrxONZaJi/I4MkvN5AUGcprN6cQG97E6bJEfN6VfWN4YfZWXv8xXaFOROQYVmXu55mvN7Fw296fpycY1UujB0vd0kyy4lFKKyr5/UdrePyLDYzoEsUnvxmsQCfiIRoF+HPDwDi+35RL+p4ip8vxKMaYycaYXGPMumptfzPGZBtjVrn/XVDtuYeNManGmM3GmPOqtZ/vbks1xjxU37+HiJy6rTkF3PHOMi4dv4DNuwt45OIufP+7YVzeJ0aBTuqcztSJx8g9WMId7y5n5Y793Ht2EveenaRJN0U8zHUD4hg/J5W3Fmbwt0u6Ol2OJ3kTeBl4+4j2F6y1z1VvMMZ0Aa4BugJtgW+NMcnup8cDI4AsYKkxZoa1dkNdFi4ipydrXzEvfruVT1Zk0SQogAdGJPProQmENNLXbKk/+msTj7A6cz93vLOcA4fKmXB9H0Z211xYIp4oMjSYi3q05aPlWTx4bjKhwZpTCcBaO98YE3+Ci48CplprS4F0Y0wq0N/9XKq1Ng3AGDPVvaxCnYgH2lNYyvg5qby3aAcYuHVoAncN70B40yCnSxMfpFAnjvt0ZTZ//HgNrUIa8fFdg+nStpnTJYnIcdwyOJ7pK7P5eHkWtwxJcLocTzfOGHMTsAx40Fq7D4gGFlVbJsvdBpB5RPuAml7UGDMWGAsQFxdX2zWLyHEcOFTOq/PTmLwgndIKF1f2ieHec5Jo20JzzYlzFOrEMZUuyzOzNvGfeWn0TwhnwvV9aBnSyOmyROQX9IxtQe+4Frz103ZuGhSvbtLHNgF4HLDu238Cv66NF7bWTgImAaSkpNjaeE0ROb7isgreXJjBxLnbOFhSwUU92nD/iGQSI0KcLk1EoU6ccbCknHvfX8mczXncMDCORy7uSqC/xu0R8Ra3DI7n3qmrmLclj7M6RTpdjkey1uYcvm+MeRX4wv0wG4ittmiMu43jtIuIQ8oqXLy/ZAf//j6VPYWl/KpTJA+em0zXts2dLk3kZwp1Uu/S8gq57e1l7NhbzOOXduPGge2cLklETtIF3dvw95kbeWNhhkLdMRhj2lhrd7kfXgYcHhlzBjDFGPM8VQOlJAFLAAMkGWMSqApz1wDX1W/VInJYRaWLT1Zm89K3W8nef4j+CeFMvEETh4tnUqiTejVvSx7jpqwgwM/w7m0DGNi+pdMlicgpCPT344YB7fjn7C2k5hbSIdK3ux8ZY94HhgOtjDFZwCPAcGNML6q6X2YAdwBYa9cbY6ZRNQBKBXC3tbbS/TrjgFmAPzDZWru+fn8TEXG5LF+u3cUL324hLa+I7tHN+fvl3TkzqZUmDhePpVAn9cJay+s/pvP3mRtJjgrl1Zs0obiIt7t2QBz//j6Vt3/K4LFR3Zwux1HW2mtraH79OMs/CTxZQ/tMYGYtliYiJ8hay3cbc/nn7C1s3HWQ5KgQJt7Ql/O6RinMicdTqJM6V1JeyZ+nr+PjFVmc37U1/xzdk6aau0XE67UKacTFPaumN/jdeR1ppukNRMRLLUjdw3PfbGbljv20a9mEF6/uxcU922rScPEa+mYtdSr3YAlj31nOqkxNKC7SEN0yOJ6PV2QxbWkmt53R3ulyREROyrKMfJ77ZjOL0vJp0zyYf1zenSv7xmjwNvE6CnVSZ1Zn7mfsO8soKKnQhOIiDVT3mOaktAvj7Z+2M2ZIgo5qi4hXWJO1n39+s4V5W/JoFdKIRy7uwrX94wgO9He6NJFT4lioM8acD7xE1cXgr1lrn3KqFql9n67M5g8fryEytGpC8c5tNKG4SEM1ZkgCd09ZwZxNuZzTJcrpckREjmnT7oO8MHsLs9bn0KJJIA+N7MRNg9rRJEjnOcS7OfIXbIzxB8YDI4AsYKkxZoa1doMT9UjtqT6h+ICEcCbc0JfwpkFOlyUidejcrlG0aR7MmwszFOpExCNtyyvkxW+38sWanYQEBXD/Ocn8emg8oboWWBoIpw5L9AdSrbVpAMaYqcAoqoZ3rnXFZRU8OmMDZ3eO5NyuretiFUJVoBv79jK+25TLjQPb8X8Xd1GfdBEfEOjvxw0D2/HsrM1sySkgOSrU6ZJERADYvreIl77byqcrs2kU4M+dwxK548z2tGiiA87SsDgV6qKBzGqPs4AB1RcwxowFxgLExcWd1sqCA/z5fM1OGgf5K9TVoW835vDdplz+dEEnxp6Z6HQ5IlKPru0fx5sLM9iWW6hQJyKOy95/iJe/38qHy7Lw9zP8ekgCdw5PpFVII6dLE6kTHtuB2Fo7CZgEkJKSYk/ntfz8DEmRIWzJKaiV2qRmby7IILpFY349JMHpUkSknoU3DeKnh35FgM7Oi4iDcg6WMH5OKlOXVJ07uH5AHL85qwNRzYIdrkykbjkV6rKB2GqPY9xtdSYpKpR5W/LqchU+bdPug/yUtpeHR3bSlzoRH6X3vog4Ja+glInztvHuou1UuixXpcQy7lcdiG7R2OnSROqFU6FuKZBkjEmgKsxdA1xXlytMjgrho+VZ7C8uUz/qOvDWwgyCA/24ul/sLy8sIiIiUgvyi8r4z7xtvPVTBuWVlst7R/PbXyUR17KJ06WJ1CtHQp21tsIYMw6YRdWUBpOttevrcp1J7ms8tuQU0j8hvC5X5XP2FZUxfWU2l/WOUWAWERGROre/uIxXf0jjzQUZFJdXcmmvaO45O4mEVk2dLk3EEY5dU2etnQnMrK/1Jf8c6goU6mrZ1KWZlJS7uGVwvNOliIiISAN24FA5k39MZ/KP6RSWVXBh9zbcd04SHSI1QJP4No8dKKW2tW0eTEijALZqsJRaVVHp4p2fMhic2JKOrfWBKiIiIrWvoKScNxZk8NoPaRwsqeD8rq25b0QSnVo3c7o0EY/gM6HOGEOHyBC25BQ6XUqDMntDDjsPlPC3S7o6XYqIiIg0MIWlFby1MINXf0hjf3E5I7pEcd85SXRt29zp0kQ8is+EOoCkyBDmbM51uowG5Y2FGcSENebszlFOlyIiIiINRHFZBW//tJ1J89PILyrjV50iuf+cZLrHKMyJ1MSnQl1yVCgfLs8iv6iM8KYa0ON0rd95gCXp+fz5gs74+xmnyxEREREvd6iskncXbWfivG3sLSpjWHIE949IpldsC6dLE/FoPhXqkqJCgKrBUga2b+lwNd7vrYUZNA70Z3SKpjEQERGRU1dSfjjMpbGnsJShHVpx/4gk+rbT4HYiJ8KnQt3hETC3KtSdtvyiMj5dtZOr+sbQvEmg0+WIiIiIFyopr+T9JTt4Ze428gpKGZzYkgk39KFfvMKcyMnwqVDXpnkwoY0CNFhKLXh/yQ7KKjSNgYiIiJy8kvJKPliayStzU8k5WMrA9uG8fG1vBuigu8gp8alQZ4yhQ1QIWzStwWkpr3Tx7qLtDO3Q6udJ3UVERER+SWlFJdOWZjJ+zjZ2HyyhX3wYL1zdi8GJrZwuTcSr+VSoA0iODOXbjTlOl+HVvlmfw64DJTw+qpvTpYiIiIgXKK2oZNqyLF6Zk8quAyWktAvjn6N7MjixJcZosDWR0+VzoS4pKoQPlmWyt7CUliGNnC7HK725MJ248Cac1SnS6VJERETEg5VVuPhweSbjv09l54ES+rYL49krezKkg8KcSG3yuVB3eLCULTmFDFKoO2nrsg+wNGMff7lQ0xiIiIhIzcoqXHy0PIvxc1LJ3n+I3nEteOqKHpyR1EphTqQO+Gyo25pbwKBEXYx7st5YkEGTIH+u0jQGIiIicoSyChcfr8ji5e+rwlyv2BY8eVk3hiVHKMyJ1CGfC3VRzRoRGhygwVJOQV5BKZ+v3sk1/WNp3ljTGIiIiEiV8sqqM3MKcyLO8LlQZ4whOSpU0xqcgimLd1BW6eJmTWMgIiIiVIW5T1Zk8e/vU8nad4iesS144rJuDFeYE6lXPhfqAJKjQvh63W6stfrAOUFlFS7eXbydYckRJEaEOF2OiIiIOOioMBfTnMcvVZgTcYpPhrqkyFDeL85kT2EZEaEaLOVEzFy7i7yCUsZcGe90KSIiIuKQGsPcqG4M76gwJ+Iknwx1Pw+WklOgUHcCrLW8sSCd9hFNOTMpwulyREQ8jjFmMnARkGut7eZuCwc+AOKBDGC0tXafqfrm+xJwAVAM3GKtXeH+mZuBv7hf9glr7Vv1+XuIHIvCnIhn83O6ACckR1V1H9RgKSdmxY79rM46wJjB8fhpGgMRkZq8CZx/RNtDwHfW2iTgO/djgJFAkvvfWGAC/BwCHwEGAP2BR4wxYXVeuchxlFe6mLpkB2c9N5c/fryWlk2DeOOWfnx69xDO6hSpQCfiIXzyTF1EaCOaBQewJVeDpZyINxdmEBocwOV9YpwuRUTEI1lr5xtj4o9oHgUMd99/C5gL/NHd/ra11gKLjDEtjDFt3MvOttbmAxhjZlMVFN+v6/pFjlRe6eLj5Vm8POe/Z+YeG9WVszoqyIl4Ip8MdYdHwNyqM3W/aPeBEr5au4tbBsfTtJFP/rmIiJyqKGvtLvf93UCU+340kFltuSx327Haj2KMGUvVWT7i4uJqsWTxdWUVVd0sq4c5dbMU8Xw++y09KSqUmWt3aQTMX/DOogxc1moaAxGR02CttcYYW4uvNwmYBJCSklJrryu+68hJwxXmRLyLz4a65KgQ3l9STl5hKZGhwU6X45FKyiuZsngH53SOIja8idPliIh4mxxjTBtr7S5398pcd3s2EFttuRh3Wzb/7a55uH1uPdQpPuyoMBfbgicuVZgT8TY+HOoOj4BZqFB3DDNW7WRfcTljhiQ4XYqIiDeaAdwMPOW+/axa+zhjzFSqBkU54A5+s4C/Vxsc5Vzg4XquWXxEWYWLj5ZnMX5OtTCnScNFvJbPhrqkaiNgDunQyuFqPI+1lskL0unUOpSB7cOdLkdExKMZY96n6ixbK2NMFlWjWD4FTDPG3ApsB0a7F59J1XQGqVRNaTAGwFqbb4x5HFjqXu6xw4OmiNSWI8Ncr9gWPHlZN4YpzIl4NZ8NdREhjWjRJJAtORoBsyaL0vLZtLuAp6/org95EZFfYK299hhPnV3Dsha4+xivMxmYXIuliQBVYe7D5Zm8MmebwpxIA+Szoc4YQ3KkRsA8ljcXphPWJJBRvWoceE1ERES8wJFhrndcC/5+eXfOTGqlMCfSgPhsqIOqLpifr96pETCPsGNvMd9syOE3wxMJDvR3uhwRERE5SaUVlXy4LItX5qSy80AJfeJa8I/Lu3OGwpxIg+TToS45KpSDJRXkFpQS1UyDpRz2xsJ0AvwMNw2Kd7oUEREROQmlFZVMW5bFBHeY69sujKev7MHQDgpzIg2ZT4e66oOlKNRVOVhSzrSlmVzUo622iYiIiJcorahk2tJMXpm7jV0KcyI+x6dD3eFpDbbkFHJGUoTD1XiGD5ZkUlRWya1DNY2BiIiIpyspr2TaskwmuMNcSrswnr2yJ0M6tFSYE/EhPh3qWoU0IrxpkAZLcauodPHmwgwGJITTLbq50+WIiIjIMZSUV/LB0qowt/tgCf3iw3juqp4MTlSYE/FFPh3qAJIiQ9iiUAfArPU5ZO8/xCMXd3G6FBEREalBSXklU5fsYMK8beQcLKVffBjPj+7JIIU5EZ9WZ6HOGPM34HYgz930J2vtTPdzDwO3ApXAPdbaWXVVxy9Jjgrl01XZGgETeO3HNNq1bMLZnaOcLkVERESqKSmv5P0lO5gwdxu5BaX0Twjnhat7Mai9wpyI1P2Zuhestc9VbzDGdAGuAboCbYFvjTHJ1trKOq6lRklRIRSUVJBzsJTWzX13YJAVO/axcsd+Hr2kK/5+2jmIiIh4gpLySqYsrjozl1dQyoCEcF66pjeDEls6XZqIeBAnul+OAqZaa0uBdGNMKtAf+MmBWkiKPDxYSoFPh7rXf0ynWXAAV/aNcboUERERn1dSXsl7i3cw0R3mBrYP518KcyJyDHUd6sYZY24ClgEPWmv3AdHAomrLZLnb/ocxZiwwFiAuLq7OCuzUuirUrd95kDOTfXMEzKx9xXy1dhe3n9Gepo18/jJLERERxxwqq+S9xduZOC+NPYWlDGrfkn9f25uB7RXmROTYTusbvDHmW6B1DU/9GZgAPA5Y9+0/gV+f6GtbaycBkwBSUlLs6dR5PGFNg0iMaMrSjHzuIrGuVuPR3lqYgTGGmwfHO12KiIiITzoyzA1ObMn463ozQGFORE7AaYU6a+05J7KcMeZV4Av3w2wgttrTMe42x/RPaMkXa3ZS6bI+dz1ZYWkFU5dkckH3NrRt0djpckRERHxKcVkF7y3awX/mb2NPYRlDOrTklbP70D8h3OnSRMSL1OXol22stbvcDy8D1rnvzwCmGGOep2qglCRgSV3VcSIGJITz/pIdbNp9kK5tfWt+tmlLMykordBk4yIiIvWouKyCdxdt5z/z0thbVBXmJpyTTL94hTkROXl1eQHVM8aYXlR1v8wA7gCw1q43xkwDNgAVwN1OjXx52OGjYUvS830q1FW6LG8sTKdvuzB6xbZwuhwREZEGr7isgnd+2s6k+VVh7oykVtx7dhIpCnMichrqLNRZa288znNPAk/W1bpPVtsWjYkJa8yS9HzGDPGdM1azN+SQmX+Ih0d2droUERGRBq2gpJz3Fu9g0vw08t1h7r5zkujbTmFORE6fhjp0658QzrzNeT41CfnrP6YRE9aYc7tosnEREZHaVl7pYv6WPKavzGb2hhxKK1wMS47gnrOT6NsuzOnyRKQBUahzG5AQzicrstmWV0SHyBCny6lzy7fnszRjH/93URcC/P2cLkdERKRBsNayJusA01dm8/nqnewtKiOsSSBX94vlqr6xdI/xncs8RKT+KNS59U+oGjJ4SXq+T4S6CXPTaNEkkGv6x/7ywiIiInJcmfnFfLoym+mrsknLKyIowI8RnaO4rHc0wzpGEKgDqCJShxTq3OJbNiEitBFL0vdy3YC6m+zcE2zJKeDbjTnce3YSTYL0JyAiInIqDhwqZ+baXUxfkc2SjHygqufPHWe2Z2T3NjQLDnS4QhHxFfpG72aMoX9COIvT8xv8dXX/mZdG40B/TTYuIiJyksorXczb7L5ObmMOZRUuEiOa8vvzOjKqV1tiwpo4XaKI+CCFumoGJITz5ZpdZO07RGx4w/xQzt5/iM9WZXPDwHaENw1yuhwRERGPZ61lXfZBPl6RxYzVO8kvKqNl0yCu6x/H5X2i6R7dvEEfDBYRz6dQV031+eoaaqh7/Yd0AG47w3embhARETkVuw+UMH1lNp+syGJrbiFB/n6M6BLF5X2iOTNZ18mJiOdQqKsmOTKU5o0DWZKezxV9Y5wup9btKyrj/SU7uKSnuoeIiIjUpKS8klnrd/PR8ix+TN2DtdC3XRh/v6w7F3ZvQ/Mmuk5ORDyPQl01fn6GfvHhP1/s3NC8/dN2DpVXcsewRKdLERER8RjWWlZl7ufD5Vl8vnonBSUVRLdozG/P6sDlfWKIb9XU6RJFRI5Loe4IAxLC+XZjDrkHS4hsFux0ObWmuKyCNxemc07nSDq2DnW6HBEREcflFpQwfUU2Hy7PIjW3kOBAP0Z2a8NVfWMY2L4lfn66Tk5EvINC3RF+vq4uI5+LerR1uJraM21pJvuKy7lTZ+lERMSHVbos87fm8f7iHXy3KZdKl6VPXAv+cXl3LuyhaQhExDsp1B2ha9tmNAnyZ0l6wwl15ZUuXv0hnX7xYaTEhztdjoiISL3bdeAQ05ZmMW1ZJtn7D9EqJIjbzkjgqr6xdIgMcbo8EZHTolB3hAB/P/q2C2NJesO5ru6LNTvJ3n+Ix0Z1dboUERGfY4zJAAqASqDCWptijAkHPgDigQxgtLV2n6kaF/8l4AKgGLjFWrvCibobgkqXZe7mXN5fsoPvN+XisnBGUiv+fGFnzukcRVCARq8UkYZBoa4GAxLCee6bLewvLqNFE++ey81ay8S5aXSMCuWsjpFOlyMi4qvOstbuqfb4IeA7a+1TxpiH3I//CIwEktz/BgAT3LdyEg4Ul/PBsh28s2g7mfmHiAhtxF3DE7k6JY64lhr9WUQaHoW6GvRPaAnA0ox9jOgS5XA1p2fO5lw25xTw/OieuuBbRMRzjAKGu++/BcylKtSNAt621lpgkTGmhTGmjbV2lyNVeplNuw/y1sIMpq/MpqTcRf+EcB4e2ZkRXaI0p5yINGgKdTXoEdOcoAA/lqTv9epQZ61l/JxtRLdozMU9G8b1gSIiXsgC3xhjLPAfa+0kIKpaUNsNHN7ZRAOZ1X42y932P6HOGDMWGAsQFxdXh6V7vkqXZfaGHN5cmM6itHwaBfhxWe9obhoUT5e2zZwuT0SkXijU1SA40J9esS28/rq6uVvyWL59H09e1k1HKEVEnDPUWpttjIkEZhtjNlV/0lpr3YHvhLmD4SSAlJSUk/rZhqK0opLpK7KZND+NtD1FRLdozEMjO3F1SixhTb370gkRkZOlUHcMAxLCeWXuNgpLKwhp5H2byVrLP7/ZTGx4Y67qG+t0OSIiPstam+2+zTXGTAf6AzmHu1UaY9oAue7Fs4HqH9ox7jZxKygpZ8riHbz+Yzq5BaV0i27G+Ov6cF7XKAJ0AFNEfJQ+/Y6hf0I4lS7Liu37nC7llMxan8O67IPce3ayRvcSEXGIMaapMSb08H3gXGAdMAO42b3YzcBn7vszgJtMlYHAAV1PVyWvoJRnZ21i8FPf84+vNpEUFcK7tw7g83FDubBHGwU6EfFp3ncKqp70iQvD38+wJD2fM5MjnC7npFS6LM/P3kxiRFMu6x3tdDkiIr4sCpheNVMBAcAUa+3XxpilwDRjzK3AdmC0e/mZVE1nkErVlAZj6r9kz7K/uIyJ89J4c2E6pRUuRnZrzZ3DEukR08Lp0kREPIZC3TE0bRRAt+jmXnld3RdrdrIlp5CXr+uNv0a8FBFxjLU2DehZQ/te4Owa2i1wdz2U5vEKSyt448d0Js1Po7Csgkt7RfPbX3WgfYQmChcROZJC3XEMSAjnzQUZlJRXEhzo73Q5J6Si0sWL326lU+tQLujWxulyRERETkpJeSXvLd7BK3NS2VtUxrldonjw3I50bB3qdGkiIh5Loe44+sWHM2l+Git27GNwYiunyzkhn6zIJn1PEa/elKJ56URExGu4XJaPVmTxwuwt7DpQwtAOrXjw3GR6x4U5XZqIiMdTqDuOQYktaRzoz+erd3pFqCutqOSl77bSM6Y553SOdLocERGRE7Imaz9//Ww9qzP30yu2Bf+8qieDO3j+fldExFMo1B1HSKMALujehs9X7+KvF3WhSZBnb65pSzPJ3n+If1zeHfdF+SIiIh4rv6iMZ2dtYurSTFqFNOKFq3tyaa9o7cNERE6SZ6cUD3B1v1g+XpHFl2t2cVWK5873VlJeyb+/T6V/QjhnJOnopoiIeK5Kl+X9JTt47pvNFJRUcOuQBO49J4nQ4ECnSxMR8UoKdb+gX3wY7Vs15cNlWR4d6t5dtJ3cglL+fW1vHeEUERGPtSpzP3/5dC3rsg8ysH04j43qRnKUBkERETkdCnW/wBjDVSmxPP31JtLyCj1yKOXC0gpembuNM5JaMaB9S6fLEREROUpZhYuXvtvChLnbiAhtxL+u7c3FPdroQKSISC3wc7oAb3BFn2j8/QzTlmU5XUqN3vgxnfyiMh48t6PTpYiIiBxlS04Bl72ygPFztnFFnxi+fWAYl/Rsq0AnIlJLFOpOQGSzYM7qGMHHK7KoqHQ5Xc7/yMwvZvzcVM7v2ppesS2cLkdERORnLpfl1flpXPTvH9l9oIRJN/bl2at66to5EZFaplB3gkanxJJXUMqczXlOl/Izay2PzFiPnzH838VdnC5HRETkZ5n5xVz76iKenLmRM5MimHX/mZzbtbXTZYmINEi6pu4EndUpklYhjZi2LJMRXaKcLgeAWetz+H5TLn+5sDNtWzR2uhwREREAPl2ZzV8+XQfAM1f24Kq+MepqKSJSh07rTJ0x5ipjzHpjjMsYk3LEcw8bY1KNMZuNMedVaz/f3ZZqjHnodNZfnwL9/biibzTfb8olt6DE6XIoLK3gbzPW07lNM24ZHO90OSIiIpRXunj08/Xc98EqurRpxlf3nsHolFgFOhGROna63S/XAZcD86s3GmO6ANcAXYHzgVeMMf7GGH9gPDAS6AJc617WK4xOiaXSZflkRbbTpfDC7C3kFJTw5GXdCPBXL1oREXFWflEZN72+hDcWZDBmSDzv3T6A2PAmTpclIuITTqv7pbV2I1DTEbhRwFRrbSmQboxJBfq7n0u11qa5f26qe9kNp1NHfUmMCKFffBjTlmZyx5ntHTvyuC77AG8sSOe6/nH0iQtzpAYREZHD1mUf4I53lpNXWMpzV/Xkyr4xTpckIuJT6uoUTzSQWe1xlrvtWO1HMcaMNcYsM8Ysy8vznMFJRqfEkraniGXb9zmy/kqX5c+friO8aRB/OK+TIzWIiIgc9tmqbK6cuBCXtXx4xyAFOhERB/xiqDPGfGuMWVfDv1F1WZi1dpK1NsVamxIREVGXqzopF/ZoQ9Mgfz5YmvnLC9eBKUt2sDpzP3+5sAvNm2hIaBERcUaly/KPmRu5d+oqukc3Z8a4ofTU1DoiIo74xe6X1tpzTuF1s4HYao9j3G0cp90rNAkK4OKebfls1U7+dklXQhrV3wCiuQUlPPP1JoZ0aMmoXm3rbb0iIiLVlZRXMm7KCr7dmMuNA9vx14u6EBSg67tFRJxSV5/AM4BrjDGNjDEJQBKwBFgKJBljEowxQVQNpjKjjmqoM6P7xXKovJIvVu+s1/U+8cVGSstdPD6qm0YSExERRxSWVjDmjaV8tymXx0Z15fFLuynQiYg47HSnNLjMGJMFDAK+NMbMArDWrgemUTUAytfA3dbaSmttBTAOmAVsBKa5l/UqvWNbkBQZwgfL6q8L5vwtecxYvZO7hifSPiKk3tYrIiJy2P7iMq5/bTFLMvJ5fnRPbhoU73RJIiLC6Y9+OR2YfoznngSerKF9JjDzdNbrNGMMV/eL5YkvN7Jix746H4EyY08R905dSWJEU+4anlin6xIREalJbkEJN762hPQ9RUy4vg/ndm3tdEkiIuKm/hKn6KqUWNo2D+ae91dyoLi8ztazr6iMMW8uxRjD6zf3IzjQv87WJSIiUpOsfcWMnvgTO/KLmXxLPwU6EREPo1B3ipo3DuTl6/uQc7CE3320Gmttra+jpLySse8sI3v/IV69qS/xrZrW+jpERESOZ1teIVdN/In8ojLevW0AQ5NaOV2SiIgcQaHuNPSJC+PhkZ2ZvSGH135Ir9XXdrksf/hoDUsz9vH86J70bRdeq68vIiLySzbuOsjoiT9RVuHi/bED6duubi83EBGRU6NQd5rGDIlnZLfWPPX1JpZl5Nfa6z4/ewszVu/kD+d35KIemr5ARETq1469xdz4+hIC/f2YducgurZt7nRJIiJyDAp1p8kYw9NX9iAmrDHjpqxkb2Hpab/mtKWZvDwnlWv7x3LXMA2MIiIi9WtPYSk3TV5MeaWLd2/rT6JGXRYR8WgKdbWgWXAgr1zfh/ziMu77YBWVrlO/vu7HrXv40/S1nJHUisc0H52IiNSzotIKfv3mUnYfLGHyLSl0iAx1uiQREfkFCnW1pGvb5jx6SVd+2LqH8XNST+k1lmbkc9e7y+kQGcIr1/ch0F//PSIiUn/KKlzc+e5y1u88yMvX9tH13CIiXuK05qmT/3VNv1iWpOfzwrdb6NsujCEdTmyEsC05BTw3azPfbMihbfNgJt/Sj9DgwDquVkRE5L+qBuhazQ9b9/D0Fd05p0uU0yWJiMgJUqirRcYYnrysG+uyDzBuygquGxDHsORI+sS1IKCGs25Z+4p58dutfLIii6ZBAfzu3GTGDEmgaSP9t4iI+CpjzPnAS4A/8Jq19qn6WO8/vtrIp6t28rtzk7m6X1x9rFJERGqJ0kMtaxIUwMQb+/Ln6WuZOC+N8XO2ERocwBlJrRieHMmZyREE+hvGz9nGu4u2g4Fbhybwm+EdCGsa5HT5IiLiIGOMPzAeGAFkAUuNMTOstRvqcr2vzk/j1R/SuWlQO+4+q0NdrkpEROqAQl0dSIwIYerYQRwsKWfB1j3M25LH3M15zFy7G4Agfz8qXC6u6hvLveck0bZFY4crFhERD9EfSLXWpgEYY6YCo4A6C3WfrszmyZkbuaB7ax65uKsG6BIR8UIKdXWoWXAgI7u3YWT3Nlhr2ZxTwLzNeezcf4gbB7XTiGIiInKkaCCz2uMsYMCRCxljxgJjAeLiTq+rZOvmwZzbJYrnR/fC30+BTkTEGynU1RNjDJ1aN6NT62ZOlyIiIl7OWjsJmASQkpJy6vPoAAPbt2Rg+5a1UpeIiDhDY+aLiIh4jmwgttrjGHebiIjIMSnUiYiIeI6lQJIxJsEYEwRcA8xwuCYREfFw6n4pIiLiIay1FcaYccAsqqY0mGytXe9wWSIi4uEU6kRERDyItXYmMNPpOkRExHuo+6WIiIiIiIgXU6gTERERERHxYgp1IiIiIiIiXkyhTkRERERExIsZa09rztJ6YYzJA7bXwku1AvbUwus0JNomR9M2OZq2ydG0TWpWG9ulnbU2ojaK8QW1tI/U3/PRtE1qpu1yNG2To2mbHK1O949eEepqizFmmbU2xek6PIm2ydG0TY6mbXI0bZOaabt4J/2/HU3bpGbaLkfTNjmatsnR6nqbqPuliIiIiIiIF1OoExERERER8WK+FuomOV2AB9I2OZq2ydG0TY6mbVIzbRfvpP+3o2mb1Ezb5WjaJkfTNjlanW4Tn7qmTkREREREpKHxtTN1IiIiIiIiDYpCnYiIiIiIiBfziVBnjDnfGLPZGJNqjHnI6XqcYoyZbIzJNcasq9YWboyZbYzZ6r4Nc7LG+maMiTXGzDHGbDDGrDfG3Otu99ntYowJNsYsMcasdm+TR93tCcaYxe730QfGmCCna61vxhh/Y8xKY8wX7sc+vU2MMRnGmLXGmFXGmGXuNp9973gr7SO1f6yJ9o9H0/7x2LR/PFp97yMbfKgzxvgD44GRQBfgWmNMF2ercsybwPlHtD0EfGetTQK+cz/2JRXAg9baLsBA4G7334cvb5dS4FfW2p5AL+B8Y8xA4GngBWttB2AfcKtzJTrmXmBjtcfaJnCWtbZXtbl3fPm943W0j/zZm2j/eCTtH4+m/eOxaf9Ys3rbRzb4UAf0B1KttWnW2jJgKjDK4ZocYa2dD+Qf0TwKeMt9/y3g0vqsyWnW2l3W2hXu+wVUfSBF48PbxVYpdD8MdP+zwK+Aj9ztPrVNAIwxMcCFwGvuxwYf3ybH4LPvHS+lfSTaP9ZE+8ejaf9YM+0fT0qdvX98IdRFA5nVHme526RKlLV2l/v+biDKyWKcZIyJB3oDi/Hx7eLuRrEKyAVmA9uA/dbaCvcivvg+ehH4A+ByP26JtokFvjHGLDfGjHW3+fR7xwtpH3ls+lt20/7xv7R/rNGLaP9Yk3rdRwbU1guJ97PWWmOMT85xYYwJAT4G7rPWHqw6yFTFF7eLtbYS6GWMaQFMBzo5W5GzjDEXAbnW2uXGmOEOl+NJhlprs40xkcBsY8ym6k/64ntHGiZf/lvW/vF/af/4v7R/PK563Uf6wpm6bCC22uMYd5tUyTHGtAFw3+Y6XE+9M8YEUrXDes9a+4m72ee3C4C1dj8wBxgEtDDGHD4Q5GvvoyHAJcaYDKq6p/0KeAnf3iZYa7Pdt7lUfbnpj9473kb7yGPz+b9l7R+PTfvHn2n/eAz1vY/0hVC3FEhyj8ITBFwDzHC4Jk8yA7jZff9m4DMHa6l37n7frwMbrbXPV3vKZ7eLMSbCfQQSY0xjYARV11LMAa50L+ZT28Ra+7C1NsZaG0/VZ8j31trr8eFtYoxpaowJPXwfOBdYhw+/d7yU9pHH5tN/y9o/Hk37x6Np/1gzJ/aRxtqGf9bcGHMBVf19/YHJ1tonna3IGcaY94HhQCsgB3gE+BSYBsQB24HR1tojLxZvsIwxQ4EfgLX8ty/4n6i6bsAnt4sxpgdVF+/6U3XgZ5q19jFjTHuqjsKFAyuBG6y1pc5V6gx395LfWWsv8uVt4v7dp7sfBgBTrLVPGmNa4qPvHW+lfaT2jzXR/vFo2j8en/aP/+XEPtInQp2IiIiIiEhD5QvdL0VERERERBoshToREREREREvplAnIiIiIiLixRTqREREREREvJhCnYiIiIiIiBdTqBMREREREfFiCnUiIiIiIiJe7P8BZIdKpAw7oTEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "delta_z=(Lz/(Nz-1))*4\n",
    "print(delta_z)\n",
    "condense_dz = delta_z*weights_condense - delta_z\n",
    "unshrink_coef = Lz/(Lz-2*delta_z)\n",
    "condense_dz = (z+condense_dz)*unshrink_coef-z\n",
    "condense_dz -= weights*condense_dz[k_border]\n",
    "ax = plt.figure(figsize=(15,4)).subplots(1,2)\n",
    "ax[0].plot(condense_dz)\n",
    "ax[1].plot((z-condense_dz))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#apply offsets\n",
    "for k in range(0,Nz):\n",
    "    vertices[:,:,k,2] += (interp_dz*weights[k] - condense_dz[k])\n",
    "#todo: make condense offset and curve offset aware of each other or act one after another or something. \n",
    "#Right now they can sometimes clip through each other and result in negative cell height."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_heights = vertices[:,:,1:,2]-vertices[:,:,:-1,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.976833179364064"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.amin(cell_heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n"
     ]
    }
   ],
   "source": [
    "\n",
    "N_NODES = Nx*Ny*Nz\n",
    "N_CELLS = (Nx-1)*(Ny-1)*(Nz-1)\n",
    "\n",
    "with open('salt_index.vtk', 'w') as f:\n",
    "    # write header\n",
    "    f.write('# vtk DataFile Version 3.0\\n')\n",
    "    f.write('Generated by mfem-seismo\\n')\n",
    "    f.write('ASCII\\n')\n",
    "    f.write('DATASET UNSTRUCTURED_GRID\\n')\n",
    "    \n",
    "    # write point coordinates\n",
    "    f.write(' '.join(['POINTS', str(N_NODES), 'double\\n']))\n",
    "    #for k in range(Nz-1, -1, -1):\n",
    "    for k in range(0, Nz):\n",
    "        for j in range(0, Ny):\n",
    "            for i in range(0, Nx):\n",
    "                f.write(' '.join([str(vertices[i][j][k][0]), str(vertices[i][j][k][1]), str(vertices[i][j][k][2]), '\\n']))\n",
    "                \n",
    "    #fill cell vertices\n",
    "    f.write(' '.join(['CELLS', str(N_CELLS), str(9 * N_CELLS), '\\n']))\n",
    "    Nx -= 1\n",
    "    Ny -= 1\n",
    "    Nz -= 1\n",
    "    for i in range(Nz):\n",
    "        for j in range(Ny):\n",
    "            for k in range(Nx):\n",
    "                f.write(' '.join(['8', str((Nx+1)*(Ny+1)*i+j*(Nx+1)+k), str((Nx+1)*(Ny+1)*i+j*(Nx+1)+k+1), # 8 means num_of_nodes_per_cell\n",
    "                    str((Nx+1)*(Ny+1)*i+(j+1)*(Nx+1)+k+1), str(((Nx+1)*(Ny+1)*i+(j+1)*(Nx+1)+k)),\n",
    "                    str((Nx+1)*(Ny+1)*(i+1)+j*(Nx+1)+k), str((Nx+1)*(Ny+1)*(i+1)+j*(Nx+1)+k+1),\n",
    "                    str((Nx+1)*(Ny+1)*(i+1)+(j+1)*(Nx+1)+k+1), str(((Nx+1)*(Ny+1)*(i+1)+(j+1)*(Nx+1)+k)),\n",
    "                    '\\n'])) #\n",
    "    \n",
    "    #fill cell types for somr reason\n",
    "    f.write(' '.join(['CELL_TYPES', str(N_CELLS), '\\n']))\n",
    "    for _ in range(N_CELLS):\n",
    "        f.write('12\\n') # magic constant for VTK_QUAD\n",
    "    \n",
    "    f.write(' '.join(['CELL_DATA', str(N_CELLS), '\\n']))\n",
    "    f.write('SCALARS index int\\n')\n",
    "    f.write('LOOKUP_TABLE default\\n')\n",
    "    \n",
    "    \n",
    "    #for _ in range(0, k_border * (Nx * Ny)):\n",
    "    #    f.write('1\\n')\n",
    "    #for _ in range(k_border * (Nx * Ny), N_CELLS):\n",
    "    #    f.write('2\\n')\n",
    "    for k in range(0, Nz):\n",
    "        print(k)\n",
    "        for j in range(0, Ny):\n",
    "            for i in range(0, Nx):\n",
    "                #f.write(''.join([str(interp_ot_c[k-k_border,j,i]) if k >= k_border else '1500','\\n']))\n",
    "                f.write(''.join([str(k*Ny*Nx+j*Nx+i+1),'\\n']))\n",
    "    \n",
    "    Nx += 1\n",
    "    Ny += 1\n",
    "    Nz += 1 #so that i don't have to restart kernel every time i rerun this block\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n"
     ]
    }
   ],
   "source": [
    "with open('salt_values.txt', 'w') as f:\n",
    "    Nx -= 1\n",
    "    Ny -= 1\n",
    "    Nz -= 1\n",
    "    for k in range(0, Nz):\n",
    "        print(k)\n",
    "        for j in range(0, Ny):\n",
    "            for i in range(0, Nx):\n",
    "                f.write(''.join([str(interp_ot_c[k-k_border,j,i]) if k >= k_border else '1500','\\n']))\n",
    "                #f.write(''.join([str(k*Ny*Nx+j*Nx+i+1),'\\n']))\n",
    "    Nx += 1\n",
    "    Ny += 1\n",
    "    Nz += 1 #so that i don't have to restart kernel every time i rerun this block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}