{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Target in a solenoid\n", "\n", "In this example, we set up a solenoid source and position a target, which is a conductive, permeable ellipsoid in the center. We simulate in the frequency domain. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Install packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Basic python packages\n", "import numpy as np \n", "import matplotlib.pyplot as plt\n", "from scipy.constants import mu_0\n", "\n", "import ipywidgets\n", "from matplotlib.colors import LogNorm\n", "\n", "# Modules of SimPEG we will use for forward modelling\n", "from SimPEG import Mesh, Utils, Maps\n", "from SimPEG.EM import FDEM\n", "from pymatsolver import Pardiso as Solver\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## model parameters\n", "\n", "Here, we define a model consisting of a uniform background, and a target" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# conductivities (S/m) of the wholespace and target\n", "sigma_wholespace = 1e-3\n", "sigma_target = 1\n", "\n", "# relative permeability of target and wholespace\n", "mur_wholespace = 1\n", "mur_target = 1" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# semi-axes of the target\n", "radius_a = 40 # horizontal axis\n", "radius_c = 20 # vertical axis " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.00000000e+00 3.72759372e+00 1.38949549e+01 5.17947468e+01\n", " 1.93069773e+02 7.19685673e+02 2.68269580e+03 1.00000000e+04]\n" ] } ], "source": [ "# frequencies at which to simulate\n", "freqs = np.logspace(0, 4, 8)\n", "print(freqs)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Minimum skin depth (in target): 5.00e+00 m \n", "Maximum skin depth (in background): 1.58e+04 m \n" ] } ], "source": [ "# print the min and max skin depths to make sure mesh is fine enough and \n", "# extends far enough \n", "\n", "def skin_depth(sigma, f):\n", " return 500./np.sqrt(sigma * f)\n", "\n", "print(\n", " 'Minimum skin depth (in target): {:.2e} m '.format(skin_depth(sigma_target, freqs.max()))\n", ")\n", "print(\n", " 'Maximum skin depth (in background): {:.2e} m '.format(skin_depth(sigma_wholespace, freqs.min()))\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set up a mesh" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "cs = 2.5 # core cell size\n", "ncx = np.ceil(500/cs) # number of core cells in the x direction\n", "ncz = np.ceil(1000/cs) # number of core cells in the z direction\n", "\n", "pf = 1.3 # padding factor (how much we expand the cells)\n", "npadx = 29 # number of padding cells in the x direction (need to be > max skin depth)\n", "npadz = 30 # number of padding cells in z-direction\n", "\n", "# vectors of cell sizes in each direction\n", "hx = Utils.meshTensor([(cs, ncx), (cs, npadx, pf)])\n", "hz = Utils.meshTensor([(cs, npadz, -pf), (cs, ncz), (cs, npadz, pf)])\n", "\n", "# create a cylindrical mesh (has geometry, and differential operators)\n", "mesh = Mesh.CylMesh(\n", " [hx, 1, hz], # cylindrically symmetric mesh\n", " x0='00C' # origin at the axis of symmetry in 0, centered in z\n", ")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEKCAYAAADJvIhZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4HNWd6P3vT5sXSd6wLWwgXmIJhi0JGLDJhfSFBAhMwjA3yROSgO9MHpwJcDN5ZjI35IUAw4RJJjOE9w2BJGRgQogJIcmASWxW220Tb3jHC9hqSZZsy7asraVuSa1ezvtHV7er291Su1tdJdu/z/PoUfepU+ec+rWko6o6dY4YY1BKKaWcUOJ2A5RSSp05tNNRSinlGO10lFJKOUY7HaWUUo7RTkcppZRjtNNRSinlGO10lFJKOUY7HaWUUo7RTkcppZRjytxuwGgzdepUM3v27Lz2DQaDVFZWjmyDTmEajxNpTFJpPFKdyvHYsmVLuzFm2nD5tNNJM3v2bDZv3pzXvl6vF4/HM7INOoVpPE6kMUml8Uh1KsdDRJpzyaeX15RSSjnGtU5HRMaKyLsiskNEdovIP1vpc0Rko4jUi8hvRaTCSh9jvfdZ22fbyvqOlb5XRG60pd9kpflE5D6nj1EppVQqN890QsB1xpiPAB8FbhKRBcC/AY8bY2qBLuCrVv6vAl3GmHnA41Y+RORC4IvARcBNwFMiUioipcCTwKeBC4HbrbxKKaVc4lqnY+IC1tty68sA1wG/t9KfA/7Ken2r9R5r+/UiIlb6i8aYkDGmCfABV1pfPmNMozFmEHjRyquUUsolrg4ksM5GtgDziJ+VNADdxpiIleUgcI71+hzgAIAxJiIifuAsK32DrVj7PgfS0q/K0o7FwGKAmpoavF5vXscTCATy3vd0pPE4kcYklcYj1ZkQD1c7HWNMFPioiEwCXgb+IlM267tk2ZYtPdNZXMYV64wxTwNPA8yfP9/kO3rkVB55UgwajxNpTFJpPFKdCfEYFaPXjDHdgBdYAEwSkURneC7Qar0+CJwHYG2fCHTa09P2yZaulFLKJa6d6YjINCBsjOkWkXHAJ4kPDlgFfI74PZhFwFJrl1et9+ut7SuNMUZEXgVeEJEfATOBWuBd4mdAtSIyBzhEfLDBl4p1PI8u28Mv3gnC68uy5jln0jg+ct5E5k2rYmxFKcbA3KmV7Dsa4IrZk9l2oJu6mmoqK0rZ3NzFJedMZExZCZGYYXNzFwvmTEmWFY4ZtjZ3cZUtLaGhPcj48lJmTBybsR3BwSj7jvbysfMmDXtcPQMR9ncEufSciTlEIdUf9g1Sfm57xlPRk3Gwq5+oMcyaMr7AkjJr9Q8wEI4yd2rxH8rb0xGlwteefN83GGVvjp+FW1o6+xCB8yaPfPzT4+GUD470MrV6DFMrKxyveyhuxQNgQ2MH/+f6WspLi3su4ubltRnAc9Z9nRLgJWPMn0RkD/CiiHwP2AY8Y+V/BnheRHzEz3C+CGCM2S0iLwF7gAhwj3XZDhG5F3gDKAWeNcbsLtbB/OKdpmHzHOru51B3f8ZtC+eexfrGDqZWVVA7vZr1jR3MOms8MyeOIzgY4b2DfjbNPSuZ398fZs/hHhba0hLWN3Yky8ykvq2X9sBg1u122w50MRCO5ZT3xHaEOWrqKZXCup3hjqdQxS7frrt7kHc6fMn3vmMBjvWGHKk7X8WMT3o8nOLkZ34y3IoHxGNy17VzT99OxxjzHvCxDOmNxEeepacPAJ/PUtajwKMZ0pcDywtu7AgbU1bC3GlVhMJRvnfbxXzzxe38ZvECPvWj1Tz55cuoq6nm6u+v4IW7FnDOpHEEQhGuevRtfrN4QbIMf1+Ya364MiUt4f6Xd3LBjAncsWBWxvr3He3lniVbM+6bbmtLF4/8cU9OedPNuW8Zv7lrAaUlhXU6T6yoJxSJ8a0bzy+onGz+a20TzR19PPzZi4pSvl38mv3xWPraevna81vyiq9THntzL+WlJXzj+toRLzs9Hk65e8kWbrlkJrdcOsPxuofiVjwA6h54regdDoySezpKKaXODNrpKKWUcox2OkoppRyjnY5SSinHaKejlFLKMdrpKKWUcox2OkoppRyjnY5SSinHaKejlFLKMdrpKKWUcox2OkoppRyjnY5SSinHaKejlFLKMdrpKKWUcox2OkoppRyjnY5SSinHaKejlFLKMdrpKKWUcox2OkoppRyjnY5SSinHiDHG7TaMKvPnzzebN28+6f1m37esoHrnTa/C1xZIeT2+opSZk8YRDEU47B9g3vSqZH5/f5hjvaGUtAR7OZkMtz3fvJn2nTu1kpISOel9R6oNo6F8u75gkPGVla7Una9itjE9Hk4ZrXF3Kx4Qj8l7D9/AhLHlee0vIluMMfOHy1eWV+mqYOWlQjhqeOTWi3hw6W5++uXL+NTja3jyS5dRV1PFpx5fw1NfvoxzJo0jEIpw21Pr+OmXL0vu7+8P87mfrU9JS7j/lV38xdnVfGXBrIx17zsa4J4XtmbcN922lm4e+dOenPKm+9Tja3jqK5dRKoV1Ok+s9BGKRPnWDecXVE42/7VuP80dQR7+zEVFKd/u3U2buPKK47H0tQX4+pLcPgu3PPbmPsrLSvjGdfNGvOz0eDjl7iVbueXSGdxyyQzH6x6KW/GA+O9reUnxL35pp+OCMWUlzJ1WRSgcZd70KqZXj6G2ppra6VXU1lRRW1PNzIljqa2pTnY6lRWl1NZUJ8vw94WZMLYsJS2hdnoV82qqM24DMFaebNvtekMR5uWYN50AtdOrKS3wTKd2ehWhSCyvNuRafkVpSdHKtztUlVqPCHx4WqUjdeertqaK8iLFJz0eTqmtqaJ2evbfEbe4FQ+AirISCvz/MCd6T0cppZRjtNNRSinlGNc6HRE5T0RWicj7IrJbRP7eSp8iIm+JSL31fbKVLiLyYxHxich7InKZraxFVv56EVlkS79cRHZa+/xYxImTR6WUUtm4eaYTAf7RGPMXwALgHhG5ELgPWGGMqQVWWO8BPg3UWl+LgZ9CvJMCHgKuAq4EHkp0VFaexbb9bnLguJRSSmXhWqdjjDlsjNlqve4F3gfOAW4FnrOyPQf8lfX6VuBXJm4DMElEZgA3Am8ZYzqNMV3AW8BN1rYJxpj1Jj4u/Fe2spRSSrlgVNzTEZHZwMeAjUCNMeYwxDsmYLqV7RzggG23g1baUOkHM6QrpZRyietDpkWkCvgD8E1jTM8Qt10ybTB5pGdqw2Lil+GoqanB6/UO0+rCxGIxAoEA4ahhx/YdDA4O4vV6Cfb1sendTbRWlzAQCrFh/XrOGldCf8QQjUZT2hUMGyKRSMa2traGKO09gnegKWP9h3pjBPsGcjpOX3eUnp7BPGNiWL3aS0mBt9Ka9g8SjoLXe7igcrLx7Q9ztC+G13usKOXbBQKBlFi2BmL05fhZuGV/8yBlAl7voREvOz0eTmlrG2D3ni4qO/c6XvdQ3IoHxP8urVmzhorS4t76drXTEZFy4h3OEmPMf1vJR0VkhjHmsHWJrM1KPwicZ9v9XKDVSvekpXut9HMz5D+BMeZp4GmIz0jg8XgyZRva67nPSFBSUkJVVfw5nY989GIqPtiOx+OhcutqrrjyMupqqhm7fgULFi5MPqdTuuZt7O3y94UpW7uSTG19q2sntTMm4Mn6cGgvlfVb8Xg+MWxbJ7R08adDe/B4Pp7z8SW9voxPfMJT8HM6O6P1hCIxPJ7iPBzatLaJ8o4+PJ7iPxzq9XpTPjNfWy/j927J+DmOFlsG91JeWoLHUzviZafHwykvHdrCRRfOxHPp6Ho41K14AJS8/RrXXnstY8tLi1tPUUsfgjWS7BngfWPMj2ybXgUSI9AWAUtt6Xdao9gWAH7r8tsbwA0iMtkaQHAD8Ia1rVdEFlh13WkrSymllAvcPNP5OHAHsFNEtltp/w/wA+AlEfkq0AJ83tq2HLgZ8AF9wN8AGGM6ReRfgE1WvkeMMZ3W668DvwTGAa9ZX0oppVziWqdjjPkzme+7AFyfIb8B7slS1rPAsxnSNwMXF9BMpZRSI2hUjF5TSil1ZtBORymllGO001FKKeUY15/TOROFIjHeP9wDwHdf2UVbb4gHXtlJfVuAB17ZRV1NFa3+Ab77yi5mThpLMBQlOBjlgVd2Jsvw90foGYikpCUs2dgCwN4jPRnr33c0YNV14r7ptrV0s7u1J6e86QzwwCu7KC3wX5tfb4gfT3f/YGEFDVN+JBYrSvl2rYdCvN19PJa+tgANx4J5xdcpifi09Q6MeNnp8XDK8p1H2NDYyfrGdsfrHopb8QAYjMSIxIq/qKd2Oi47/+xqGo4FOd9aQ+P8mmrqauKrGdbVVHPOpLEEQtHktgR/f5g/pqWllFtTnXWbILzb1Jl1u11fKMru1p6c8mZuR1XBz+lMrx5DKBLLuw3DmTutkuaOvqKVbyc9JdTZ6ikVYUNjbp+FWyaNL6e8tKQobUyPh5PqaqpGXdzdjAdkH9k1krTTcYF9EbevLJjF5v1d3LFwNr9a38wdC2dRV1PNT70N3LFwVvLh0J+srOeOhbOTZfj7wjy9pjElLeGDI71cMGMCdwzxcOiGxo6M+6a76JwuGtuDOeVN9+DS3dyxcHbBnU53X5hQJJZXG3IRiRmaO/qKVr6dN7Qfj60eX1svf/a1O1J3vtp6Q5SXlhSljenxcMr6xg5uuWQmt4y2h0NdigfAvyx7v+Df1VzoPR2llFKO0U5HKaWUY7TTUUop5RjtdJRSSjlGOx2llFKO0U5HKaWUY7TTUUop5RjtdJRSSjlGOx2llFKO0U5HKaWUY3QaHBfYJ/x8Y9cR2npDvL7rMPVtAV7fdYTGYwFa/QO8vutIcu614GCU13cdTpbh7w/TMxBJSUt4Y/cRmjv6mFZVkbH+xISfmfZNt62lm+0HunPKm84Ar+86UvCEn6/vPkIoEuPicyYUVlC28nfF47Vg7pSilG+360iEAVssExN+5hNfp7y+6wjlpSXJOQFHUno8nBL/uSwp+GdzpLkVD4hP+BkzOuHnae/lbYdSvr+y7RC11i/3K9sOJWeZtueBeKeTnpbQHhjkz752KseUZqyz/mgg677ptrV055w3k5e3HSr4F3t3a09BbRjOxqbOopZv194eoT58vB5fW+6fhVvqi9jG9Hg4JWbgjztaGYxEHa97KG7FIyGqs0yfnuwTfn7vtov55ovb+fkd8/nUj1bz5Jcvo66mmqu/v4Kf3XF5csLPqx59m5/fMT9Zhr8vzDU/XJmSlnD/yzuHnfDzniVbM+6bbmtLF4/8cU9OedPNuW8ZP7/j8oInEXxiRT2hSIxv3Xh+QeVk819rm2ju6OPhz15UlPLtvF4vHs/xWPraevna81vyiq9THntzL+WlJXzj+toRLzs9Hk65e8mW0Tnhp0vxAKh74DXKHTj1G2Unl0oppU5n2ukopZRyjHY6SimlHKOdjlJKKcdop6OUUsox2ukopZRyjKudjog8KyJtIrLLljZFRN4SkXrr+2QrXUTkxyLiE5H3ROQy2z6LrPz1IrLIln65iOy09vmxiBR/AXCllFJZuf2czi+BnwC/sqXdB6wwxvxARO6z3n8b+DRQa31dBfwUuEpEpgAPAfOJPwS/RUReNcZ0WXkWAxuA5cBNwGsOHNeQQpEYh/39DISjHO4eoK03RGt3P4f9A7R291M1poxW/wCHu/sRIBiKEByM0trdnywjMSOBPS3hsH+AiePKM24DUuoazuHuAQ77+3PKm85YdRX6nE6rf4BQJJpXG3KRiEWxyrfr6I+l1NPaPZDzZ+GW1u4BKsqkKG1Mj4dTWgv4uS4mt+IB8RkJnCDGgWkPhmyAyGzgT8aYi633ewGPMeawiMwAvMaY80Xk59br39jzJb6MMV+z0n8OeK2vVcaYC6z02+35spk/f77ZvHnzSR/H7PuWnfQ+djMmjuWwfyDj60AoQu9AhBkTxybz+/vD9A1GU9IS7PtmMtz2fPNm2nd69ZiCO51C2jAayrcLhUKMGTPGlbrzVcw2psfDKaM17m7FA+Ix2fHQDUwcV57X/iKyxRgz7JOtbp/pZFJjjDkMYHU80630c4ADtnwHrbSh0g9mSHddphkJ1n/n+hNmJPjd169OmZFg/XeuT5aRmJHAnpaQ64wEb/3DJ4Zta2JGglfu+fhJH+ec+5ax/jvX64wENvEnzj3J94kZCVb8oyfrPm4r/owEnhEvdzije0YCjyt11z3wGmPKin/HZTR2Otlk+stl8kg/sWCRxcQvw1FTU4PX682zibmJxWIEAgHCUcOO7TsYHBzE6/US7Otj07ubaK0uYSAUYsP69Zw1roT+iCEajaa0Kxg2RCKRjG1tbQ1R2nsE70BTxvoP9cYI9g3kdJy+7ig9PYN5xsSwerWXkgJvpTXtHyQcBa+3OBMh+vaHOdoXw+s9VpTy7QKBQEosWwMx+nL8LNyyv3mQMgGvd+TnBEuPh1Pa2gbYvaeLys69jtc9FLfiAfG/S2vWrKGitLi3vkdjp3NURGbYLq+1WekHgfNs+c4FWq10T1q610o/N0P+Exhjngaehvjltbz+03g998trJSUlVFXFz3Q+8tGLqfhgOx6Ph8qtq7niyviZztj1K1iwcGHyTKd0zdsp/wH5+8KUrV2Z8b+it7p2UjtjAp4hznQq67fi8Qx/pjOhpYs/HdqDx3PyZzq8voxPfMJT8JnOzmj8TMfjKc6ZTtPaJso7+vB43DnTGb93i2v/3eZiy2D8TMfjOX3OdF46tIWLLpyJR890kkrefo1rr72WseWZJwoesXqKWnp+XgUSI9AWAUtt6Xdao9gWAH7rMtwbwA0iMtka6XYD8Ia1rVdEFlij1u60laWUUsoFrp7piMhviJ+lTBWRg8RHof0AeElEvgq0AJ+3si8HbgZ8QB/wNwDGmE4R+Rdgk5XvEWNMp/X668RHyI0jPmrN9ZFrSil1JnO10zHG3J5l0wl3x018mN09Wcp5Fng2Q/pm4OJC2qiUUmrkjMZ7Oqc9+8qh63wdtPWGWOtrp74twFpfO8d6Q7T6B1jra0/e0wkORlnra0+WkXhOx56WsK6hg87gIHOnVmasf9/R3mRdw9nW0sX2A9055U1ngLW+9oLv6axtaCcUibHQd1ZB5WQt39dBc0cwr2M8WXs6opTb6kmsHOpE3fla62unvLSEy2dNHvGy0+PhlLW+DiaNr2DS+PyGBxeLW/EAXTn0jPGk1xf/virxvSG5LPBTq3zMnDSOYCiSkgeOrxxqT0toag/S1B5M5km3z1o5NNO+6RIrh+aSN5MnV/kK7nQ2NHYW1IbhrGvoKGr5dl1dg/y583g9iZVDnag7X1sL/BkYSno8nOLvD/PCxhb2twcdr3sobsUjQVcOPU1lek7nhbsWnPCczpK7FqQ8p/PCXQuSZSSe07GnJeT6nE6mfdMlntPJJW+6Ofct44W7FuhzOjbx0UnHY5l4Tief+Dql+M/pOH/so/s5HXd+FnTlUKWUUqcd7XSUUko5RjsdpZRSjtF7Oi6wj17b1tJNW2+ILc1d1LcF2NrcRe9AhFb/AFubuzjiH0jOMr2luStZRo81es2elrC1pZv+wSgXzpiQsf56a/Rapn3TJUav5ZI3nSF+T6jQaXC2tnQRisTyakNu5XfT3BEsWvl2vq4o1bZ6GqzRa07Una+tLV2Ul5YUpY3p8XDK1uZuZkwcx9mjbMJPt+IB8dFrTsz/7Pos06ONW7NMf+xDk9jW0k3VmDJqa6rY1tJNzYQxydFr+44G+NiHJiXz+/vDNB4LpqQlJEacZdoGUH80QCAUybr9ZMoabt+PnDuRkgIHEhTShtFQvl2Pv4cJE4//M+BrC9A7kNtn4ZZixic9Hk5x8jM/GW7FA+Ixee/hG5gw9sybZfq0l2n02st3f3zYWaZfvvv4/GeJ0Wv2tIRcR69l2jddYvRaLnnTzblvGf9998d19JpNfHTS8VgmRq/lE1+nFH/0mvPHPrpHr7nzs1D3wGtU6Og1pZRSpxPtdJRSSjlGL6+5wD6Q4KGlu2nrDfHg0l3UtwV4cOku6mqqafUP8NDSXcy0TYPz4NJdyTIS0+DY0xKWbGyhROIDBjJJTIOTad9021q62XnIn1PedAZ4cOmugi+v/Wp9MwA9A5lnWChUonwnpgA5dCjESv/xWCamwcknvk5JxKc9EBrxstPj4ZTlO4/wblMnG5s6HK97KG7FA+IDCSI6I8Hp78PTqqhvC/DhaVXJ9/bXiU4n8T4hMcWNPS293GzbojHDhsbOrNvtevrD7DzkzylvtnYU2ulMrRpDKBLNuw3DmTO1kv0dwaKVbxfrLkmpx5j4NDxO1J2vCWPLqCgrKUob0+PhpLlD/I64xc14QOaVL0eadjousA8kuPPqWWxt6WLR1bP59YZmFl09m7qaan6+uoE7r56dHEjw1Cofi66enSzD3xfmP99pTElL2He0d9iBBO82dWbcN90l505kf0dfTnnTPfzqbhZdPbvgTqenP0woEsurDbmIGUNznsd4sryD+/HY6vG19bKuod2RuvPVHghRXlpSlDamx8MpG5s6RudAApfiAfDo8vcL/l3NhXY6LjCAMSa+drZJfW+M/bU5/mW9P17GiWn28rH2y1i/rc5h23oSeTPvbwoe+2+wjrdIl7+SMXfg8lp6PYXG1wnx+BSnjU7F/cR6i/szlS+34hGv3Jlq9DmdNG49p6OUUm7b8dANTBynz+mcdjI9p/Pu/Z8c9jmd3Y/clCwj8ZzOew/feEL5uT6n89Y/fGLYtiae03nlnvye0/H96836nI5N/DkMT/J94jmdFf/oybqP24r/nI5nxMsdzuh+TsfjSt11D7zGmLLiD2jWTscF9tFr6xvii7itsxZxW+drp91axG1d2iJu6zIs4rYuw4JP6xs66Oob5MPDLOKWad902w50s/1Ad0550xlgXUM7pQVOg7OuoYNQJMq6Ii3itq4hvohbPsd4svZ0RKmwL+J2LD56zYm687WuoYPyUmF+kRZxq3Dh2Nc1dDB5fAWTR+Eibm7EA3QRtzPGEyt9Kd9/sspH7fTq5OuZE8cRHIyk5IHjo9fsaQmN7UEa24N0BTMPMa5v6826b7ptB7pyzpvJEyt9BXc66xs7CmqD2+XbdXcP8k6HbRG3YwHH6s5XYs41YeTbmB4Pp3T3hVmysYXGY6NrETe34pGgi7idxsZXlDIQjvLobZfwnf/eyff/+hJueHwN/3rbJdTVVOP5Dy/fv+2S5JDpv3ziz3z/ry9J7u/vD3Prk2tT0hIeeGUXF5xdzVeGuLx27wvbMu6bbtuB+OW1XPKm8/yHl+//9SUFdzpPrPQRikT51g3Fubz2y3X72d8R5OHPFP/y2saNG7nqquOx9LUF+Ltfb8krvk557K19lJcK37hu5C+vpcfDKXcv2cotl87glktG1+U1t+IB8d/XspLiX17TgQRp3BpIMOus8TR39KW8Li+V5ISf7YFBZp01Ppnf3x+muy+ckpZgLyeT4bbnmzfTvudOHlfwPZ1C2jAayrfr7+9n3LhxrtSdr2K2MT0eThmtcXcrHhCPiU74eQa4fNZktjR38ZlLZ/KTVT5uuLCGuppqfrLKx2cunZk80/nluv185tKZyf38/WGe39Cckpbwk1U+Jo0vz7gN4mc6b+45mnW73bYDXaz1deSUN1M7PvORmQWf6fxkVfxyQz5tyMUv1+0nEIoUrXy75uZmZs06Xo+vLcDru484Une+ihn/9Hg45SerfMydWsnNo+xMx614QDwmhS5DkgvtdFyWuF6+aX8nAJubu5L3azY1dzKz+/g9nUQeOH5Px55m190Xzrqtvi0w5L522w5055w3k837O0fsBznfNgwnMeNDscq36+6O0maO19NwLPfPwm3FaGN6PJzU2B4cdXF3Mx7gzPNienktjT6no5Q6U+148AYm5jmiL9fLazrLtFJKKcCZSW9P+05HRG4Skb0i4hOR+9xuj1JKjVaD0VjR6zitOx0RKQWeBD4NXAjcLiIXutsqpZQanRL3GIvptO50gCsBnzGm0RgzCLwI3Opym5RSalT60i82Fr2OnDodEVkhIjenpT1dnCaNqHOAA7b3B600pZRSLsh1yPQc4NsicoUx5p+ttGFHKYwCmcbqnnCnTEQWA4sBampq8Hq9RW6WUkqNTsX++5drp9MNXA/8WET+CHyleE0aUQeB82zvzwVa0zMZY54Gnob4kOm8Znl9XYdMK6VOfcWe5TrXezpijIkYY+4G/gD8GZhevGaNmE1ArYjMEZEK4IvAqy63SSmlRqV/dmB5j1zPdH6WeGGM+aWI7ATuKU6TRo4xJiIi9wJvAKXAs8aY3S43SymlRqUbLzq76HXk1OkYY36e9n4L8LdFadEIM8YsB5a73Q6llBrtKhxYxO10HzKtlFIqR4VOzpsLnfBzlLhqzhQ2NnUypbKC2ulVbGzq5Lwp45KLuO061MNVc6Yk8/v7w3xwpDclLWFjU2eyzEzq2wJ0BgezbrfbdqCbwUgsp7yZ2nHF7MkFT/g53PEUqtjl23V3dzNp0qTk+4ZjAdoDuX0WbilmfNLj4RQnP/OT4VY8IB4TB5bT0U7HbRPHlePvDzN5fAUAk8aXJ1/Hl9OtSJ7yJtKB5B9ye1q6bNsmjS+nMzg45L7HyyjnaE8op7yZ66oYsf+e8m3DaCkfINInKfVMGl9BeyC3z8JtxWhjejycNtri7nY8nKCdjsvC1lxH9u+J14OR+OvBSGqeTPsNVfZwdQ7dPpNz3mx1xUao08m3DaOlfIBI7OQ/x9GiGG1Mj4fTRlvc3Y6HE7TTcdntV36IZ/7cxO1XfogVH7Rx+5Ufom56NSs+aONLV30oeXltY1Mnt1/5oeR+/v4wm5u7UtISdrf2cMGM6ozbAPa19fLD1/dm3W53wYxqnlzVkFPedInjKfRMpz04SCgczasNuQhFYuzvCBatfLudu/xccvHxenzHAvzgtQ8cqTtfrf4BKkqlKG1Mj4dT3vG185eXzBh1i7i5FQ+I/74WuspvLrTTcdnqfcdSv+89xuHugeTrxHLV9jxwfBE3e1rCkZ4QO2i9AAAegklEQVQBjvYOcN7kzEvx7jvam3XfdNsOdOWcN5PV+44V3OnssBaSy7cNw/mzr72o5dsdao/SaavHZy2o50Td+Xr/cA9QnDamx8Mpg5EYq/cdo3LM6PoT6FY8EmIOLK82uiJ+BqqdXoWvLUBtTVX8fU1V8vW8mirOsZarTmxL8PeFT0hLLzfbNoNhY1Nn1u12gVB8EEMuebO1o9D/nqZVjyEUjubdhuHMnVrJ/o5g0cq3M/6SlHpEYH1jhyN152viuHLKS0uK0sb0eDhp3hC/I25xMx6Qed6wkaadjgvGlJUwd1oVoXCUOxbOYktzF3cunM3z65u5c+Fs6mqq+Zm3gTsXzk52Ok+u9HHnwtnJMvx9YX7xTmNKWsLeI71cMGMCdyyYlbH+fUd72djYmXHfdBef00VTezCnvOkeWrqbOxfOLrjT8feFCUViebUhF9GYobmjr2jl23lD+/HY6vG19bLW1+5I3fk61huivLSkKG1Mj4dTNjR2cMslM7nl0tF1ec2teAB8b9n7jlxe0+Wq07i1XPU864zH/np8RWny8tph/wDzptvOdPrDHOsNpaQl2MvJZLjt+ebNtO/cqZWUFPiDXEgbRkP5dn3BIOMrK12pO1/FbGN6PJwyWuPuVjwgHpP3Hr6BCWOLu1y1num4pLxUCEcNj9x6EQ8u3c1Pv3wZn3p8DU9+6TLqaqr41ONreOrLlyXPdG57ah0//fJlyf39/WE+97P1KWkJ97+yi784u5qvZD3TCXDPC1sz7ptuW0s3j/xpT055033q8TU89ZXLCr6n88RKH6FIlG/dcH5B5WTzX+v209wR5OHPFH/eqXc3beLKK47H0tcW4OtLcvss3PLYm/soLyvhG9fNG/Gy0+PhlLuXbOWWS2dwyygbSOBWPCD++1ruwIM62um4wH55bd70KqZXj6G2pjp5H6a2ppqZE8dSW1Od7HQqK0qpralOluHvCzNhbFlKWkLt9Crm1VRn3AbxtR3idWXebtcbiljXvofPm06A2unVBZ+y106vIhSJ5dWGXMuvKC0pWvl2h6pS6xGBD0+rdKTufNXWVFn3dEa+jenxcEptTRW107P/jrjFrXhAfAocByYk0E7HDaFILDkiaJ2vg7beEGt97dS3BVjra+dYb4hW/wBrfe3JTic4GGWtNcoK4mc6PQORlLSEdQ0ddAYHmTs182n6vqO9ybqGs62li+0HunPKm84Aa33tBXc6axvaCUViLPSdVVA5Wcv3ddDcEczrGE/Wno4o5bZ6fG0BGo45U3e+1vraKS8t4fJZk0e87PR4OGWtr4NJ4yuYND6/S0nF4lY8ID6iL+bA7RbtdFz2pNcX/74q8b2BOmv0ylOrfClDphN54PiQaXtaQlN7kKb2YDJPun1HA1n3TbetpTvnvJk8ucpXcKezobGzoDYMZ11DR1HLt+vqGuTPncfrSdxbcKLufG0t8GdgKOnxcIq/P8wLG1vY3x50vO6huBWPhKgDY6a103GB/fLa9267mG++uJ0X7lrAp360mie/fBl1NdVc/f0VLLlrQfJM56pH3+aFuxYky/D3hbnmhytT0hLuf3nnsKPX7lmyNeO+6ba2dPHIH/fklDfdnPuW8cJdCwrudJ5YUU8oEuNbNxbpns7aJpo7+njYgbVEvF4vHs/xWPraevna81vyiq9THntzL+WlJXzj+toRLzs9Hk65e8mW0Tl6zaV4ANQ98BrlpTrLtFJKqdOIdjpKKaUco5fXXGAfSPB/f/8ebb0h/ul3O6hvC/BPv3+PuulVtPoH+L+/35Gcey04GOWffrcjWUZiIIE9LeF3Ww4C8J41fUy6fW2BeF0Z9k237UA3vhzzpjPAP/1+R8FDphPHc7RnoKByhis/ce+smI4cCbHs2PFY+o7FBxLkE1+nJOJzoLNvxMtOj4dTlu88wuq9x/DubXO87qG4FQ+IDySI6D2d098Vs6dwsOsQV8yewu+2HOSKWZOpq6nmd1sOMn/WlOQ9neU7j3DF7NT1dN7cczQlLWHp9lYumFGdcRvEpzXZcaA763a7slLB1xbIKW+63205yBWzpxTc6Wxo6iAUjuXVhlzsau2huSNYtPLtPgi1cYGtnsmVFWxrye2zcMua+mOUl5YUpY3p8XDK77YcZP7sKaMu7m7FA+Ix0WlwTlP2gQSfn38ua33tfOGK8/jFO4184YrzqKup5v99ex9fuOK8ZKfz2Jt7+cIV5yXL8PeFeWJlfUpawo6D3VwwY0LGbRAfSLB637Gs2+3m1VTx/uHenPKm+/Yf3uML888reCDB0Z4BQpFYXm3IRXAwQnNHX9HKt/MGG/DY6vG19bLi/aOO1J2vA119lJeWFKWN6fFwindf2+gcSOBSPAAeWLrLkWlw9J6OUkopx2ino5RSyjHa6SillHKM3tNxgX302kNLd9PWG+LBpbuobwvw4NJd1NVU0+of4KGlu5hpmwbnwaW7kmUkRq/Z0xKWbGyhRKDeWqwtXWIanEz7ptvW0s3OQ/6c8qYzwIMjcJ34V+ubAegZyDzDQqES5TsxBcihQyFW+o/HMjENTj7xdUoiPu2B0IiXnR4PpyzfeYR3mzrZ2NTheN1DcSseoKPXzhgfnlZFfVuAD0+rSr63v55pW8QtkQ7Hp8Gxp6WXm21bNGbY0NiZdbtdT3+YnYf8OeXN1o5CO52pVWMIRaJ5t2E4c6xF3IpVvl2suySlHmPi0/A4UXe+Jowto6KspChtTI+Hk+YO8TviFjfjAbqI22nLPnrtzqtnsbWli0VXz+bXG5pZdHV8Ebefr27gzquPL+L21Cofi66enSzD3xfmP99pTElL2Hd0+EXc3m3qzLhvukvOncj+jr6c8qZ7+NXdLLq68EXcevrji7jl04ZcxEx8EbdilW/nHdyPx1aPr62XdQ3tjtSdr/ZAfBG3YrQxPR5O2dg0ShdxcykeAI8ud2YRN72no5RSyjGudDoi8nkR2S0iMRGZn7btOyLiE5G9InKjLf0mK80nIvfZ0ueIyEYRqReR34pIhZU+xnrvs7bPdur4lFJKZebWmc4u4K+BNfZEEbkQ+CJwEXAT8JSIlIpIKfAk8GngQuB2Ky/AvwGPG2NqgS7gq1b6V4EuY8w84HErn1JKKRe50ukYY943xuzNsOlW4EVjTMgY0wT4gCutL58xptEYMwi8CNwqIgJcB/ze2v854K9sZT1nvf49cL2VXymllEtG20CCc4ANtvcHrTSAA2npVwFnAd3GmEiG/Ock9jHGRETEb+U/YVk+EVkMLAaoqanB6/WOxLFkFYvFCAQChKOGHdt3MDg4iNfrJdjXx6Z3N9FaXcJAKMSG9es5a1wJ/RFDNBpNaVcwbIhEIhnb2toaorT3CN6Bpoz1H+qNEewbyOk4fd1RenoG84yJYfVqLyUF9vVN+wcJR8HrPVxQOdn49oc52hfD6z1WlPLtAoFASixbAzH6cvws3LK/eZAyAa/30IiXnR4Pp7S1DbB7TxeVnZn+93WPW/GA+N+lNWvWUFFa3P/Ni9bpiMjbwNkZNt1vjFmabbcMaYbMZ2RmiPxDlXViojFPA08DzJ8/33g8nizNG8Lry3LOGo7Bgd4YAI9tCRGJGZ5pGEdrIMiPd8aoq6mkcyDIk3tKkqPXBqJ9PNMwLlmGvz9MX4SUtIR3DsRXQ2wKT8hY/76jvRwNmIz7ptvW0k0gFMspbzpDkP/0jSt4RMw79fHj6S6dWFA5w5UfGjOpKOXbdXb2M2XK8Vj62gIcDub2WbjlnYZ4fNpl5OOfHg+nbD4aZPPRKDtrqx2veyhuxQMgEguy8OP/g+qxxV3Cu2idjjHmk3nsdhCwz3Z3LtBqvc6U3g5MEpEy62zHnj9R1kERKQMmAp15tKmo7rp2Lj/1NnDXNXN5p76du66ZS11NNe/Ut7P4mrnJ5aq/vmQrd10zN7mfvz/M//nNtpS0hOaOPi44u5qvDDFk+nvL3s+4b7ptLd08/va+nPKme6e+nbuunVvwLNOhcIxQJJpXG3JRXlrC/o5g0cq3e++9Xi699Hg9vrYAj/xpjyN156tnIEJFqRSljenxcMr2lm5uuXQGN18yuoZMuxUPiP++OjFkerRdXnsVeEFEfgTMBGqBd4mftdSKyBzgEPHBBl8yxhgRWQV8jvh9nkXAUltZi4D11vaVxjjwyHkO7M/pXFM7lT9sOci1ddOonV7FtXXTqKupZubEsVxTNy15plNZUcq1ddOSZfj7wkwYW5aSlnBN7VQumDEh4zaAsyeOpXbTgazb7arGlrFqb1tOedMJcG3ttIJ/kHcc6CYUieXVhlw0HAvwoSnji1a+Xaw19TObOWksH55W6Ujd+dq0v5Py0pKitDE9Hk65pm4q19ROG3VxdyseABVlJQVfCs+FW0OmbxORg8BCYJmIvAFgjNkNvATsAV4H7jHGRK2zmHuBN4D3gZesvADfBv5BRHzE79k8Y6U/A5xlpf8DkBxmrZRSyh2unOkYY14GXs6y7VHg0Qzpy4HlGdIbiY9uS08fAD5fcGOVUkqNGJ2RQCmllGO001FKKeUY7XSUUko5RjsdpZRSjtFORymllGNG23M6ZwT7yqH/8cZe2npD/PsbH1DfFuDf39hLXU0Vrf4B/uONvcycNJZgKEpwMMq/v/FBsozEyqH2tIQlG1uYOK6cI/7+jPXvOxqw6jpx33TbWrrZfqA7p7zpDPDvb+yltMB/bZ5c1WCVV5zHrH65dj/BwSiVY0qLUr5dc/Mgm0LHY5lYOTSf+DolEf9QJDriZafHwynLdx7h/cO97Dnsd7zuobgVD9CVQ88Y4yrif+jGlR//nng91nodjZGSB+I/IOlpKeXaysm0bah9882bbf9CO50xZSWEIrG82zCccRWlBAejRSvfrqI0NZaFxtcJIvFZG4rRxvR4OGnsEL8jbnEzHqArh5627DMS3PM/51F/NMC919WydHsr9143j7qaal7Y2MK9181Lzkjw/Pr93HtdbbIMf1+Y3246kJKWcNg/MOzKoe8f7sm4b7qtLV10BAdzypvusTf3ce918wqekcCY+NlhPm3IReWYMpo7+opWvp3XewiP53g9vrZedh7yO1J3vkKRGOWlJUVpY3o8nLLncM/oXDnUpXgA/HilT1cOVUopdXrRTkcppZRjtNNRSinlGO10lFJKOUY7HaWUUo7RTkcppZRjtNNRSinlGO10lFJKOUY7HaWUUo7RTkcppZRjdBocF9gn/NzQ0EFbb4h1De3UtwVY39BBeyBEq3+A9Q0dKRN+rmtoT5bRY034aU9LWN/YQXdfmA9Pq8xYf7014WemfdMlJvzMJW86A6xv6KCkwH9t1jd2EIrE8mpDTuU3dNDc0Ve08u3e74hSYaunwZrw04m687W+oYPy0hLmz5484mWnx8Mp6xs6mFJZweTKcsfrHopb8YD4fI4xoxN+nvZ+vNIX/76iPvm9tqYq+TrR6djzAPj7IyekJTQeC9J4LEhHMJSxzvqjgaz7ptvW0p1z3kx+vKK+4E5nQ2NnQW1wu3y77u5BVnccr8fXFnSs7nxtbu4CijPLd3o8nNLVF+bXG1rwtQUcr3sobsUjIaqzTJ+e7BN+fu+2i/nmi9t5cfFCPvWj1Tz55cuoq6nm6u+v4DeLFyQn/Lzq0bd5cfHCZBn+vjDX/HBlSlrC/S/vHHbCz3uWbM24b7qtLV088sc9OeVNN+e+Zfxm8YKCJxF8YkU9oUiMb914fkHlZPNfa5to7ujj4c9eVJTy7bxeLx7P8Vj62nr52vNb8oqvUx57cy/lpSV84/piTPiZGg+n3L1kyyid8NOdeADUPfAa5YVOCZ8DvaejlFLKMdrpKKWUcox2OkoppRyjnY5SSinHuNLpiMi/i8gHIvKeiLwsIpNs274jIj4R2SsiN9rSb7LSfCJyny19johsFJF6EfmtiFRY6WOs9z5r+2wnj1EppdSJ3DrTeQu42BhzKbAP+A6AiFwIfBG4CLgJeEpESkWkFHgS+DRwIXC7lRfg34DHjTG1QBfwVSv9q0CXMWYe8LiVTymllItc6XSMMW8aYyLW2w3AudbrW4EXjTEhY0wT4AOutL58xphGY8wg8CJwq4gIcB3we2v/54C/spX1nPX698D1Vn6llFIuGQ33dP4WeM16fQ5wwLbtoJWWLf0soNvWgSXSU8qytvut/EoppVxStIdDReRt4OwMm+43xiy18twPRIAlid0y5Ddk7hzNEPmHKitTWxcDiwFqamrwer2Zso2YWCxGIBAgHDXs2L6DwcFBvF4vwb4+Nr27idbqEgZCITasX89Z40rojxii0WhKu4JhQyQSydjW1tYQpb1H8A40Zaz/UG+MYN9ATsfp647S0zOYZ0wMq1d7KSnwBLNp/yDhKHi9hwsqJxvf/jBH+2J4vceKUr5dIBBIiWVrIEZfjp+FW/Y3D1Im4PUeGvGy0+PhlLa2AXbv6aKyc6/jdQ/FrXhA/O/SmjVrqCgt7gWhonU6xphPDrVdRBYBfwlcb0xywp+DwHm2bOcCrdbrTOntwCQRKbPOZuz5E2UdFJEyYCLQmaWtTwNPA8yfP994PJ5cDjHV68tyzlpSUkJVVXxGgo989GIqPtiOx+OhcutqrrgyPiPB2PUrWLBwYXJGgtI1b2Nvl78vTNnalWRq61tdO6mdMQHPEDMSVNZvxeP5xLBtndDSxZ8O7cHj+XjOx5f0+jI+8QlPwTMS7IzGZyTweIozI0HT2ibKO/rweJyakcCTfO9r62X83i0ZP8fRYstgfEYCj6dYMxJ4Rrzc4bx0aAsXXTgTz6ickcDjSt0lb7/Gtddey9jy0uLWU9TSsxCRm4BvA581xvTZNr0KfNEaeTYHqAXeBTYBtdZItQrigw1etTqrVcDnrP0XAUttZS2yXn8OWGnr3JRSSrnArbnXfgKMAd6y7u1vMMb8nTFmt4i8BOwhftntHmNMFEBE7gXeAEqBZ40xu62yvg28KCLfA7YBz1jpzwDPi4iP+BnOF505NKWUUtm40ulYw5izbXsUeDRD+nJgeYb0RuKj29LTB4DPF9ZSpZRSI2k0jF5TSil1htClDVxgAGOs1UlM6ntj7K/N8S/r/fEyTkyzl4+1X8b6bXUO29aTyJt5f0Ohd9IM1vEW6ZZcMuYO3PJLr6fQ+DohHp/itNGpuJ9Yb3F/pvLlVjzilTtTjYy2oLtt/vz5ZvPmzSe93+z7ch+9ppRSo9GOh25g4rj8VlMVkS3GmPnD5dMzHRdkWsTt3fs/ecIibr/7+tUpi7jtfuSmZBmJRdzee/jGE8rPdRG3t/5h+CHTiUXcXrnn5IdMz7lvGb5/vVkXcbPJNGT6a89vYcU/erLu47biL+LmGfFyhzO6F3HzuFJ33QOvMaZMF3FTSil1GtFORymllGO001FKKeUY7XSUUko5RjsdpZRSjtFORymllGO001FKKeUY7XSUUko5RjsdpZRSjtFORymllGO001FKKeUY7XSUUko5RjsdpZRSjtFORymllGO001FKKeUY7XSUUko5RjsdpZRSjtFORymllGO001FKKeUY7XSUUko5psztBpyJQpEY7x/uAeDBpbtp6w3x3Vd2Ud8W4Luv7KKupppW/wAPvrKLmZPGEQxFCA5G+e4ru5Jl+PvD9AxEUtISlmxsQQT2HenNWP++o73Juoaz7UAXuw715JQ3nQG+u3QXpSInva/d8xuagfgxF0Oi/GjMFKV8u0OtIVZ0H4+lry1Aw7FgXvF1SiI+x3pDI152ejycsnznETY2drKhscPxuofiVjwABiMxIg78DrjS6YjIvwC3AjGgDfjfxphWERHg/wNuBvqs9K3WPouAB6wivmeMec5Kvxz4JTAOWA78vTHGiMgU4LfAbGA/8AVjTFexjmnGxLEc9g/klPf8mmomjisnEotx1dyzGFdeSm1NFZfPmszH501l0vhyPnLuRC6fPZmqMWUMRmLsOdxDbU1VsoxQOIavLZCSlnDzJWczeXxFxm0A06rHEImZrNvtJo0vp2pMWU55082sEuqmV1FSUlin87nLzyWaY3vzcfuVH6JnIFy08u2MvySlnmnVYxiMxhypO19//bFzEJGitDE9Hk65/oLpzDqrktlTxzte91DcigfE/y6Vlxb2u5oTY4zjX8AE2+tvAD+zXt8MvAYIsADYaKVPARqt75Ot15Otbe8CC619XgM+baX/ELjPen0f8G+5tO3yyy83+Vq1alXe+56ONB4n0pik0nikOpXjAWw2OfyNdeWejjGmx/a2kviVGIif/fzKOoYNwCQRmQHcCLxljOk08bOVt4CbrG0TjDHrrYP+FfBXtrKes14/Z0tXSinlEtfu6YjIo8CdgB/4n1byOcABW7aDVtpQ6QczpAPUGGMOAxhjDovI9JE+BqWUUienaJ2OiLwNnJ1h0/3GmKXGmPuB+0XkO8C9wEPEL5GlM3mkn2xbFwOLAWpqavB6vSdbBACBQCDvfU9HGo8TaUxSaTxSnQnxKFqnY4z5ZI5ZXwCWEe90DgLn2badC7Ra6Z60dK+Vfm6G/ABHRWSGdZYzg/iAhWxtfRp4GmD+/PnG4/Fkyzokr9dLvvuejjQeJ9KYpNJ4pDoT4uHKPR0RqbW9/SzwgfX6VeBOiVsA+K1LZG8AN4jIZBGZDNwAvGFt6xWRBdbItzuBpbayFlmvF9nSlVJKucStezo/EJHziQ+Zbgb+zkpfTnwEm4/4kOm/ATDGdFrDrDdZ+R4xxnRar7/O8SHTr1lfAD8AXhKRrwItwOeLeUBKKaWG50qnY4z5X1nSDXBPlm3PAs9mSN8MXJwhvQO4vrCWKqWUGkk6DY5SSinHSPzkQiWIyDHil/zyMRVoH8HmnOo0HifSmKTSeKQ6leMxyxgzbbhM2umMIBHZbIyZ73Y7RguNx4k0Jqk0HqnOhHjo5TWllFKO0U5HKaWUY7TTGVlPu92AUUbjcSKNSSqNR6rTPh56T0cppZRj9ExHKaWUY7TTGSEicpOI7BURn4jc53Z7iklE9ovIThHZLiKbrbQpIvKWiNRb3ydb6SIiP7bi8p6IXGYrZ5GVv95apO+UICLPikibiOyypY3Y8YvI5VZ8fda+Dqyslb8s8XhYRA5ZPyPbReRm27bvWMe2V0RutKVn/B0SkTkistGK029FpMK5ozt5InKeiKwSkfdFZLeI/L2Vfsb+jKTIZdEd/Rp2UbpSoAGYC1QAO4AL3W5XEY93PzA1LS3jonnksTDfaP8CrgUuA3YV4/jJsjDhaP3KEo+HgW9lyHuh9fsxBphj/d6UDvU7BLwEfNF6/TPg624f8zDxmAFcZr2uBvZZx33G/ozYv/RMZ2RcCfiMMY3GmEHgReKLyJ1Jsi2ad1IL8znd6HwYY9YAnWnJI3L8MvTChKNSlnhkcyvwojEmZIxpIj7P4pVk+R2y/oO/Dvi9tf+oX5DRGHPYGLPVet0LvE98na8z9mfETjudkZFtkbnTlQHeFJEtEl+LCNIWzQMSi+ad7MJ8p6qROv6hFiY81dxrXS56NnEpiZOPx1lAtzEmkpZ+ShCR2cDHgI3ozwignc5IGZHF5E4hHzfGXAZ8GrhHRK4dIm9RF+A7BbiyMOEo8FPgw8BHgcPAY1b6GRMPEakC/gB80xjTM1TWDGmnZUxAO52Rkm3xudOSMabV+t4GvEz80shR67QfSV00b6iF+U6nmI3U8Q+1MOEpwxhz1BgTNcbEgF8Q/xmBk49HO/HLTWVp6aOaiJQT73CWGGP+20rWnxG00xkpm4Baa5RNBfBF4ovInXZEpFJEqhOviS+ot4vsi+ad1MJ8Dh7KSBuR4zdDL0x4ykj8cbXcRvxnBOLx+KKIjBGROUAt8ZviGX+HrHsWq4DPWfuP+gUZrc/tGeB9Y8yPbJv0ZwR09NpIfREfgbKP+Aic+91uTxGPcy7xkUU7gN2JYyV+7X0FUG99n2KlC/CkFZedwHxbWX9L/EayD/gbt4/tJGLwG+KXjMLE/+v86kgePzCf+B/pBuAnWA9xj9avLPF43jre94j/UZ1hy3+/dWx7sY26yvY7ZP3MvWvF6XfAGLePeZh4/A/il7veA7ZbXzefyT8j9i+dkUAppZRj9PKaUkopx2ino5RSyjHa6SillHKMdjpKKaUco52OUkopx2ino5RSyjHa6SillHKMdjpKjXIicoU1ceZYa0aI3SJysdvtUiof+nCoUqcAEfkeMBYYBxw0xnzf5SYplRftdJQ6BVjzkW0CBoCrjTFRl5ukVF708ppSp4YpQBXxlSjHutwWpfKmZzpKnQJE5FXiq2nOIT555r0uN0mpvJQNn0Up5SYRuROIGGNeEJFSYJ2IXGeMWel225Q6WXqmo5RSyjF6T0cppZRjtNNRSinlGO10lFJKOUY7HaWUUo7RTkcppZRjtNNRSinlGO10lFJKOUY7HaWUUo75/wGDMqxfh/855wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot the mesh\n", "fig, ax = plt.subplots(1, 1)\n", "mesh.plotGrid(ax=ax)\n", "# ax.set_xlim([-1., 1500.]) # uncomment to zoom in\n", "# ax.set_ylim([-1000., 1000.])" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total number of cells: 105340\n" ] } ], "source": [ "print(\"total number of cells: {}\".format(mesh.nC))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# put the model on the mesh\n", "sigma_background = sigma_wholespace*np.ones(mesh.nC)\n", "mur_background = mur_wholespace*np.ones(mesh.nC)\n", "\n", "# ellipsoid\n", "sigma = sigma_background.copy()\n", "mur = mur_background.copy()\n", "\n", "# indicies in the mesh where the ellipsoid is\n", "inds_ellipsoid = (mesh.gridCC[:,0]**2/radius_a**2 + mesh.gridCC[:,2]**2/radius_c**2) <= 1.\n", "\n", "# models with target\n", "sigma[inds_ellipsoid] = sigma_target\n", "mur[inds_ellipsoid] = mur_target" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the model" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAAEYCAYAAABftDB3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcXFWd///XmwBBNkkIS4AEUOMCLlERGJlRQGT7joZFEBxlcYnMwG/UUQeQGXBAhFEUdXDAqBHCFhAEogTZBBkXNEEDhE3CojSJRHYUjKT78/vjnoaboqq7qrv63ltV7+fjcR+pOnepcyvwyafOPYsiAjMzMzOzXrda2RUwMzMzM6sCJ8ZmZmZmZjgxNjMzMzMDnBibmZmZmQFOjM3MzMzMACfGZmZmZmaAE2NrQNJhkn42BtedKunPksY1cexZkv6z3XVoJ0k3Svpok8eGpFeNdZ3MzGoNFatq43L+WEn/JOmaIuvaKkmfl3Rek8c2HbOtNzkx7jCSPiBpYQpiyyRdJenvy65XI5IelLTb4PuI+ENErBsR/cOdGxFHRMRJ6To7S+obRT0+nxLTf60p/2Qq//xIr21mvSPFtOdSDH5E0vckrVt2vUZjqLgcEedHxO6D70fzAz/F8ZD0g5ryN6XyG0dyXbN2cmLcQST9G/A14IvAJsBU4H+BGWXWq4P8Dji0puyQVG5m1qz3RMS6wFuAtwH/0eoFJK3e9lp1hj8Bb5e0Ya7sUByHrSKcGHcISS8HTgSOjIgfRMRfIuL5iPhhRHw2HTNe0tckLU3b1ySNT/t2ltQn6dOSlqfW5sNz199Q0jxJT0v6NfDK3L6t0q/51XNlqzyOkvQxSXdJekbSnZLeIulcsuT9h6l15d/z15J0kKSFNff5KUnz0uuzJX1B0jrAVcBm6Tp/lrSZpGfzwVXSWyX9SdIaDb7GBcDakrZNx28LvCyV5+vwMUlLJD2evpPNcvveLeluSU9JOgNQzbkfTt/DE5KulrRlo79TM+tsEfEwWWx6PWRxWtJ3U3x9OMWvwe4Jh0n6uaTTJT0OfL6m7ElJ90t6eyp/KMXqF37Mpxh/mqQ/pNbqsyS9LO2bIOlHKQY+kV5vUVPlV0r6dYpfV0iamM59SYzPfeYL3eok3ZSKb01x+P2SFkt6T+74NSQ9Kml6g6/tb8DlwEHp+HHAgcD5NZ/7dkkLUl0XSHp7bt/Wkn6a/r25FphUc+6Okn6RvtNbJe3coC5mL+HEuHP8HbAWcNkQxxwH7AhMB94EbM+qLRmbAi8HNgc+AnxT0oS075vAX4HJwIfT1hRJBwCfJ2t9XR94L/BYRHwI+AOpdSUivlRz6jzgNZKm5co+AFyQPygi/gLsBSxN11k3IpYCN5IF1EEfBOZGxPNDVPfcVE/IWinm1NzLrsAp6bqTgd8Dc9O+ScClZN/pJOA+YKfcufsAnwP2AzYC/g+4cIi6mFkHkzQF2Bv4bSo6B1gJvAp4M7A7kO/PugNwP7AxcHKu7DZgQ7LYN5esFfpVZDHtDL3YVeO/gVeTxfhXkcXy49O+1YDvAVuSNUg8B5xRU+VDyGL7Zqme32jlfiPiHenlm1Icvogshn4wd9jewLKIWDTEpebwYhzeA7gDWDq4MyXsV6b6bQh8Fbgy1xByAXALWRw+idyTQEmbp3O/AEwEPgNcKmmjVu7VepcT486xIfBoRKwc4ph/Ak6MiOUR8Sfgv4AP5fY/n/Y/HxHzgT+TJabjgP2B41NL9GKyAN+sjwJfiogFkVkSEb8f7qSIeBa4AjgYICXIryVLmJtxDikgp3s4mCzxHcp5wMGpVfmg9D7vn4DZEfGbiFgBHAv8naStyAL+nRFxSUq+vwb8MXfux4FTIuKu9Pf0RWC6W43Nus7lkp4Efgb8FPiipE3IfsB/MsXR5cDppJbRZGlE/E9ErIyI51LZAxHxvdS/9yJgClmcXhER15C1sL5KkoCPAZ+KiMcj4hmyGHMQQEQ8FhGXRsSzad/JwDtr6n1uRCxOjQ3/CRyoJgZCD+M8YG9J66f3H2KYOBwRvwAmSnoNWYI8p+aQ/wfcGxHnpu/qQuBu4D2SppL9cPjP9B3dBPwwd+4HgfkRMT8iBiLiWmAhWfw2G5YT487xGDCp3qOunM3IWjgH/T6VvXCNmsT6WWBdstbN1YGHas5t1hSy1tORuICUGJO1Fl+eEuZmXAFsI+kVwLuBpyLi10OdEBF/AJaQ/YNyb0Q8VHPIKt9hRPyZ7LvfPO17KLcvWPU72xL4enp89yTwOFlXi82bvB8z6wz7RMQGEbFlRPxLSnK3BNYAluViwLfIWocH1cYbgEdyr58DiIjassE4vTZwS+76P07lSFpb0rck/V7S08BNwAY1iW9tjF+Dmm4IrUpP734O7C9pA7IfB+cPfRaQJc9HAbvw0iehtf+WDdZ3MA4/kZL7/L5BWwIHDH5H6Xv6e7IngGbD6tXO/53ol2RdHfYBLmlwzFKyoHBHej+V3OOpIfyJ7LHaFLJf5YPnDhoMQGsDT6fXm+b2P0SuT3KNGOazryFL+KeTJcifavY6EfFXSReTtfK+luFbiwfNAWYDh9fZN/gdAqCsf/OGwMPAMrLvaHCf8u/JvoeTI6KZfxTMrLs8BKwAJg3xZG+4eDiUR8mS5G1T3+ZanwZeA+wQEX9MMfW3rDoOIh+vppI9RXy0pnwkziF7crg68MsG9at1LlkjxZyIeDYLpy9YJQ7n6vtjsjg8QdI6ueR4Ki9+tw+RtYx/bER3Yj3PLcYdIiKeIutL9k1J+6TWgTUk7SVpsO/uhcB/SNoo9Yc9npd2Fah37X7gB2SDQdaWtA25PlupW8bDwAcljZP0YVZNhL8DfEbZ4DdJelWu+8AjwCuG+OyVZIn+l8n6g13b4NBHgA2VDULMmwMcRtavual5LMkeV+4OXFxn3wXA4ZKmKxu4+EXgVxHxIFm/tW0l7Zda7v+VVX8gnAUcqxcH97089b82sy4XEcvIfuh/RdL6klaT9EpJtd0ZRnr9AeDbwOmSNoasP62kPdIh65Elzk+mPron1LnMByVtI2ltssHclzQzdWaNejH9crIZOj7BS7tFNLqfB8i6ehxXZ/d84NXKpiddXdL7gW2AH6VueguB/5K0prLpSt+TO/c8si4Xe6R/r9ZSNvi8diCiWV1OjDtIRHwV+DeywV9/IvtlfBRZUIJssMFCsoEctwO/SWXNOIrscd0fgbPJBnHkfQz4LFm3gm2BX+Tq9X2y/mwXAM+k+kxMu08hS9aflPSZBp99AbAb8P1GLS0RcTdZ4n9/utZmqfznwADwm5S8DisinouI63J9/PL7rifre3cpWcvEK3mxD9+jwAHAqel7mEb2CHHw3MvIBsfMTY8yF5M9VjSz3nAIsCZwJ/AE2Y/+dj7CP5qslfXmFGOuI2slhmzMw8vIWoBvJmtdrXUuWXz/I9lg7n+tc8xwPg+ck+LwgZDFVLKYuTVZI0tTIuJnqStGbfljwD+StYI/Bvw78I8pBkPW7W4Hsu5qJ5BLxlP3uBlkA6EH/538LM53rEnKukmadS5JPwEuiIjvlF0XM7NeJOl44NUR8cFhDzarMPcxto4m6W1kj/C8yImZWQlS142PsOosSGYdyY8WrGNJOofsUeIn0/RE1mMk7SnpHmULshxTZ/94SRel/b9K0+6ZWZtI+hhZd4Wr0tRp1mMkzVa2GM3iBvtfK+mXklbUdqlsFMOVLeLyK0n3phi+5ljfxwuf7a4UZtaJ0jRUvyObqq+PbAXDgyPiztwx/wK8MSKOkHQQsG9EvL+UCpuZdSFJ7yBbF2FORLy+zv6NyWYZ2Ydsqr3TUnnDGJ5mnPpBRMyVdBZwa0ScWcT9uMXYzDrV9sCSiLg/Iv5GtmJYbZeaGby4WM0lwLtUMy+UmZmNXHpS8PgQ+5dHxAKy6QHz6sbwFKN35cWpac8hS6oL0dV9jNfU+FiLdcquhllXeoYnHo2IlpdZ3WOXdeKxx5ubIeqW21bcQTZ/96BZETErvd6cVRcs6CMbqZ73wjERsVLSU6RVJFutt43MuHXXidUnThz+QDNr2d8e6is7Do9Goxi+IfBkbpaqPgpcKKurE+O1WIcd9K6yq2HWla6LS1pZHfEFjz3ez6+vnjr8gcC4yff+NSK2a7C7Xstvbd+wZo6xMbT6xIls/ulG6/aY2Wg88MlPlx2HR6NRfC41bnd1Ymxm1RPAAAPtuFQfq67YtQUvXelx8Ji+tCjLyxnikZ+ZWS9oYxwejUYx/FGy5cxXT63G9WL7mHEfYzMrVBA8H/1NbcNYAExLo5fXJFuIZV7NMfN4cRXH9wE/CY84NrMe18Y4PBp1Y3iK0TeQxWzIYvgVY1mRPLcYm1nh2tFSkfoMHwVcDYwDZkfEHZJOBBZGxDzgu8C5kpaQtRQfNOoPNjPrAu1qMZZ0IbAzMElSH9lqhGsARMRZkjYlW5V3fWBA0ieBbSLi6XoxPF32aLJVZL8A/JYslhei1MRY0myyZR+XD07xkSYKvwjYCngQODAinkijFL8O7A08CxwWEb8po95mNnJB0N+mRtuImA/Mryk7Pvf6r2TLeFsDjsNmvafNcfjgYfb/kaw7RL19L4nhqfx+slkrCld2V4qzgT1ryo4Bro+IacD16T3AXsC0tM0ECpnPzszab4BoarNCnI3jsFnPcRyur9TEuMHcd/l5R/Nz180gmzw6IuJmso7Zk4upqZm1SwD9RFObjT3HYbPe4zjcWNktxvVsEhHLANKfG6fyevPdvWReO0kzJS2UtPB5Vox5Zc2sdW6pqLy2xeH+P/9lzCtrZq1zHK6vkwbfNTWvXZp0ehbA+prYe3+jZhUXwPOeGKJTtRyHx0+d4r9ss4pxHG6sionxI5ImR8Sy9IhueSpvZs5SM6u46NHHcx3GcdisizkON1bFrhT5eUfzc9fNAw5RZkfgqcFHfWbWQQL6m9ysNI7DZt3Mcbihsqdrqzf33anAxZI+AvyBF6damk82RdASsmmCDi+8wmY2atmKS1YVjsNmvcdxuLFSE+Mh5r57V51jAzhybGtkZmNP9NftqmplcBw260WOw41UsY+xmXWxAAZ68PGcmVlVOA435sTYzAoVwN8qObzBzKw3OA435sTYzAo3EH6EZ2ZWJsfh+pwYm1mhshWXHJDNzMriONyYE2MzK1Qg+v0Iz8ysNI7DjTkxNrPC+RGemVm5HIfrc2JsZoXyIzwzs3I5DjfmxNjMChWI58Ohx8ysLI7DjflbMbPCuaXCzKxcjsP1OTE2s0JFiP7woA8zs7I4DjfmxNjMCjfglgozs1I5DtfnnwtmVqhs0MdqTW1mZtZ+7YzDkmZLWi5pcYP9kvQNSUsk3SbpLal8F0mLcttfJe2T9p0t6YHcvuntvP+huMXYzArmR3hmZuVqaxw+GzgDmNNg/17AtLTtAJwJ7BARNwDTASRNBJYA1+TO+2xEXNKuSjbLibGZFSqA52Nc2dUwM+tZ7YzDEXGTpK2GOGQGMCciArhZ0gaSJkfEstwx7wOuiohn21KpUXCzjZkVanDFJXelMDMrR4txeJKkhbltZosftznwUO59XyrLOwi4sKbs5NT14nRJ41v8zBFzi7GZFW7AXSnMzErVQhx+NCK2G8VH1RvlFy/slCYDbwCuzu0/FvgjsCYwCzgaOHEUdWhaJRNjSa8BLsoVvQI4HtgA+Bjwp1T+uYiYX3D1zGwUBgd9WHU5Bpt1t4LjcB8wJfd+C2Bp7v2BwGUR8fxgQa6bxQpJ3wM+M+a1TCqZGEfEPbzYIXsc8DBwGXA4cHpEnFZi9cxsFALRH54mqMocg826W8FxeB5wlKS5ZIPvnqrpX3wwWQvxCwb7IEsSsA9Qd8aLsVDJxLjGu4D7IuL32fdjZp1uwC3GncQx2KwLtSsOS7oQ2JmsL3IfcAKwBkBEnAXMB/Ymm3XiWbIf2IPnbkXWmvzTmsueL2kjsm4Yi4Aj2lLZJnRCYlzbIfsoSYcAC4FPR8QT+YNTp/CZAGuxdmGVNLPmRMizUnSWlmIwrBqHx02YUEglzax57YzDEXHwMPsDOLLBvgd56UA8ImLXtlRuBCrdbCNpTeC9wPdT0ZnAK8ke8S0DvlJ7TkTMiojtImK7NShsEKOZNSmA/litqW00JE2UdK2ke9OfdTM0Sf25SeTnjepDu8xIYjCsGofHrbtOIXU1s+YVFYc7UdXveC/gNxHxCEBEPBIR/RExAHwb2L7U2pnZiBQ0XdsxwPURMQ24Pr2v57mImJ629472Q7uMY7BZl/K0mfVV/Y4PJvcIL03pMWhfCuyMbWbtEYiBaG4bpRnAOen1OWQDOKw1jsFmXajAONxxKtvHWNLawLuBj+eKv5TWyw7gwZp9ZtYhWmiFmCRpYe79rIiY1eS5mwyOfE6jmzducNxa6TNWAqdGxOXNVq6bOQabdbdebA1uRmUT47Qs4IY1ZR8qqTpm1iZB+yaWl3QdsGmdXce1UKWpEbFU0iuAn0i6PSLua+H8ruQYbNa9WozDPaWyibGZdaegraOhd2u0T9IjubkwJwPLG1xjafrzfkk3Am8Gej4xNrPu1c443G38c8HMCtePmtpGaR5waHp9KHBF7QGSJkgan15PAnYC7hztB5uZVV1BcbjjuMXYzAoVoaIe4Z0KXCzpI8AfgAMAJG0HHBERHwVeB3xL0gBZQ8GpEeHE2My6WoFxuOM4MTazwhUxN2ZEPEa2altt+ULgo+n1L4A3jHllzMwqphfnKG6GE2MzK1QAAz34eM7MrCochxtzYmxmBZNbKszMSuU43IgTYzMrVIBHQ5uZlchxuDEnxmZWqMEVl8zMrByOw405MTazwg14pkgzs1I5DtfnxNjMChUB/W6pMDMrjeNwY06MzaxwfoRnZlYux+H6nBibWaGyvm1+hGdmVhbH4cacGJtZobLR0A7IZmZlcRxuzN+KmRUsa6loZjMzs7HQvjgsabak5ZIWN9gvSd+QtETSbZLektvXL2lR2ublyreW9CtJ90q6SNKabbntJvhfHjMr3ABqajMzs7HRxjh8NrDnEPv3AqalbSZwZm7fcxExPW3vzZX/N3B6REwDngA+0sq9jYYTYzMr1OBo6GY2MzNrv3bG4Yi4CXh8iENmAHMiczOwgaTJjQ6WJGBX4JJUdA6wT9M3N0ruY2xmhXM3CTOzcrUQhydJWph7PysiZrXwUZsDD+Xe96WyZcBa6dorgVMj4nJgQ+DJiFhZc3whKpsYS3oQeAboB1ZGxHaSJgIXAVsBDwIHRsQTZdXRzFrnFZc6h+OwWXdqMQ4/GhHbjeLj6n1QpD+nRsRSSa8AfiLpduDpIY4fc1Vvttkl9TsZ/As5Brg+9Tm5Pr03sw4SwMpYranNKsFx2KzLFByH+4ApufdbAEsBImLwz/uBG4E3A4+SdbdYvfb4InTavzwzyPqaQMF9TsysfTwrRUdzHDbrAgXG4XnAIWl2ih2BpyJimaQJksYDSJoE7ATcGREB3AC8L51/KHBFOyrSjMp2pSD7QXONpAC+lfqzbBIRywDSl7px7UmSZpKNemQt1i6yvmbWjHBXig4y6jg8bsKEIutrZs1oYxyWdCGwM1lf5D7gBGANgIg4C5gP7A0sAZ4FDk+nvg74lqQBsobaUyPizrTvaGCupC8AvwW+25bKNqHKifFOqd/JxsC1ku5u5qQUuGcBrK+JhfVJMbPmBHgqts4x6jg8fuoUx2GzimlnHI6Ig4fZH8CRdcp/AbyhwTn3A9u3pYItqmxinOt3slzSZWRf0COSJqdWisnA8lIraWYj4hbjzuA4bNa9HIfrq2QnPknrSFpv8DWwO7CYrJ/KoemwQvucmFl7BFlAbmaz8jgOm3Uvx+HGqtpivAlwWTbHM6sDF0TEjyUtAC6W9BHgD8ABJdbRzEYgECsHKvmb3FblOGzWpRyHG6tkYpz6lrypTvljwLuKr5GZtZP7GFef47BZd3Mcrq+SibGZdbFw3zYzs1I5DjfkxNjMCjXYt83MzMrhONyYE2MzK5wDsplZuRyH63NibGaFCnpzpLOZWVU4DjfmxNjMCtfv5Z7NzErlOFyfvxUzK1REMfNnSjpA0h2SBiRtN8Rxe0q6R9ISSceM6kPNzDpAUXG4EzkxNrPCRaipbZQWA/sBNzU6QNI44JvAXsA2wMGSthntB5uZVV1BcbjjuCuFmRWsmFaIiLgLIC1Q0cj2wJI0Zy+S5gIzgDvHvIJmZqXpzdbgZjgxNrPCVagVYnPgodz7PmCHkupiZlaYCsXhSnFibGaFanH+zEmSFubez4qIWYNvJF0HbFrnvOMi4oomrl+vItFs5czMOpHnMW7MibGZFSugv/mA/GhENBw4FxG7jbI2fcCU3PstgKWjvKaZWbW1Fod7ihNjMytUUKlHeAuAaZK2Bh4GDgI+UG6VzMzGVsXicKV4VgozK1hzUwS1Ybq2fSX1AX8HXCnp6lS+maT5ABGxEjgKuBq4C7g4Iu4Y1QebmVVeMXG4E7nF2MwKFwX04o2Iy4DL6pQvBfbOvZ8PzB/7GpmZVUcRcbgTucXYzArn+TPNzMrVrjgsabak5ZIWN9gvSd9IiyjdJuktqXy6pF+mhZhuk/T+3DlnS3pA0qK0TW/bjQ/DLcZmVqgI920zMytTm+Pw2cAZwJwG+/cCpqVtB+DM9OezwCERca+kzYBbJF0dEU+m8z4bEZe0q5LNqlyLsaQpkm6QdFf6FfGJVP55SQ/nfj3sPdy1zKya+gfU1GblcBw2637tisMRcRPw+BCHzADmROZmYANJkyPidxFxb7rGUmA5sFEbbm1UqthivBL4dET8RtJ6ZL8grk37To+I00qsm5m1gVuMK89x2KzLFRiH6y2ktDmwbLBA0vbAmsB9ueNOlnQ8cD1wTESsKKCu1WsxjohlEfGb9PoZspHim5dbKzNrl6C5fm1OnsvjOGzW3VqMw5MkLcxtM1v8uCEXUpI0GTgXODwiBlLxscBrgbcBE4GjW77JEapcYpwnaSvgzcCvUtFRqYP2bEkTSquYmY1KNLlZ+RyHzbpTC3H40YjYLrfNqnvBxhoupCRpfeBK4D9SN4usbtmP80itxN8Dtm/5BkeosomxpHWBS4FPRsTTZJ21XwlMJ2t+/0qD82YO/qp5nkJa3c2sFeFZKTpFO+Jw/5//Ulh9zaxJxcbhecAhaXaKHYGnImKZpDXJptScExHfz5+QWpGRJGAfoO6MF2Ohin2MkbQGWTA+PyJ+ABARj+T2fxv4Ub1z0y+ZWQDra6IbncyqyP9nVl674vD4qVP8t21WRW36P1PShcDOZF0u+oATgDUAIuIssnni9waWkM1EcXg69UDgHcCGkg5LZYdFxCLgfEkbkXXDWAQc0Z7aDq9yiXH6dfBd4K6I+GqufHJEDHbU3pcCfz2YWXsNeMaJSnMcNut+7YrDEXHwMPsDOLJO+XnAeQ3O2bUtlRuByiXGwE7Ah4DbJS1KZZ8DDk4TPAfwIPDxcqpnZqMReFaKDuA4bNbFHIcbq1xiHBE/o/4IRi/ZatYNAnBArjTHYbMu5zjcUOUSY+tcVy+9tbTP3mOzN5X22da6cK9TM7NSOQ7X58TYzIrngGxmVi7H4boqO12bmXWr7ljgQ9LRaVU4M7MO0x1xeCy4xdhaUmZ3iaEMVy93taiQgOiOWSm+CNyWpjVbFhELyq6QmVlTuicOt50TYzMrXnc8whPZPL5PA2tL+iPwzxHhAWpmVn3dEYfbzomxraKqLcKj5RblqumaloqDIuL7ktYCDgC+J+kDEXF92RUzMxta18ThtnJibGbF646Wir8AlwBExF+BcyUtBf4LcGJsZtXWHXG47Tz4zsyKF01u1fYgsGNN2Q3ANsVXxcysRd0Rh9vOLcYGdG8XimYN3r+7VBSgeyaW/wbwfUn/AlwZEf3AfsDj5VbLzGwY3ROH286JsZkVLgbKrsHoRcS3Ja0GfAdYR9JfgA2BmeXWzMxseN0Qh8eCE+Me1uutxPXkvxO3Ho+hLmmpiIhvSfoO8A5gM2BxRPh/LDOrvi6Jw+3mxNjMCqcu6rcWEf2SbgReCWwo6VXAfRFecNXMqqub4nA7OTE2s2J12YAOSbsCpwMPAE8CE4CtJH0qIn5SauXMzOrpsjjcTk6Me4y7TzTPA/LGirrtEd5JwC4R8cKgO0kTgR8CTozNrIK6Lg63jRNjMyted7VUiGz1u7xn8Oz5ZlZl3RWH28aJcY9wS/HIeUDeGOiu0dD/C/xK0k3AE8BE4O+Br5VaKzOzoXRXHAZAkoDJEbF0pNdwYmxmxeqy+TMj4jxJPwR2IJuu7THghIh4qtyamZk10GVxeFBEhKSzgd1Heo2OW/lO0p6S7pG0RNIxZdfHzFqnaG4b1WdIB0i6Q9KApO2GOO5BSbdLWiRp4Ug+KyKeiohrIuLC9GdXJ8WOw2adr11xWNJsScslLW6wX5K+keLFbZLektt3qKR703ZorvytKS4vSee2ksX/RtJrWzh+FR2VGEsaB3wT2Its2dWDJXn5VbNOU8xSpIvJVqK7qYljd4mI6RHRMIFulaQr2nWtKnEcNusS7YvDZwN7DrF/L2Ba2mYCZ8ILg5RPIHvatj1wgqQJ6Zwz07GD5w11/VrvAn4uaa6kT0jaqYVzm+tKIel64CsRMT9XNisiil7haXtgSUTcn+owF5gB3FlwPcys4iLiLoDWGhpaJ+lz9YqB17X5cxyHzaxyIuImSVsNccgMYE6a2/1mSRtImgzsDFw7OKOPpGuBPdO88OtHxC9T+RxgH+CqJuvzNkmrA9sCbwEOBH7e7P0022K8NXC0pBNyZW1rWWnB5sBDufd9qewFkmZKWihp4fOsKLRyZtacFh7hTRr8/zltY5EEBnCNpFtGeP1PkcWih3NbH7Q9AHVkHO7/818KrZyZNafAONwoZgxV3lenvP59SGdI+rCkN6UnWoPWjojvRcQnWqlss4PvniRrmv5GGmTywVY+pI3qNf2s0tAfEbOAWQDra6InIzGrmgAGmm7FfXSo7g2SrgM2rbPruIhotivDThGxVNLGwLWS7o6IZrpfDFoEXFc7ClrSO1q4RjM6Mg6PnzrFcdisatoYh5vQKGa0Wt4v+PDnAAAbAUlEQVTIAPBR4H+A1STdDtwN/D+yAdEtaTYxVkSsBP5F0mHAz8hWdypaHzAl934LYMRTcphZSdqUKkXEbm24xtL053JJl5F1FWg6MY6Idzco/8ho61bDcdjM2qe4n6yNYkYfWXeKfPmNqXyLOsfXFRH/Ci+Mf9iGrD/yMcBZI6lss10pXrh4RJwNHAZcM5IPHKUFwDRJW0taEzgImFdCPcxsFIqYlaKpekjrSFpv8DXZFD91R1bXOXdqg22TMaqu47CZtU2BcXgecEianWJH4KmIWAZcDewuaUIadLc7cHXa94ykHdNsFIcAwz4BjIj+iLg9Ir5MNvD6FSOpbFMtxhHxrZr3twAfHskHjkZErJR0FNmXOQ6YHRF3FF0PMxulYpLefckerW0EXClpUUTsIWkz4DsRsTewCXBZGqC3OnBBRPy4yY94kAZ3ImkFMBf4ZETUroo3Io7DZtZWbYrDki4ka/mdJKmPbKaJNQAi4ixgPrA3sAR4Fjg87Xtc0klkP7YBThwciAf8M9lsFy8jG3TX1MC7nP8DLh3J/XTcAh9pRPb8YQ80s+oqIDGOiMuAy+qULyUL0qSZFUa6nOFHgQ8BJwG/B7YEjgMuJkuaTwJOI5tyqKs4Dpt1gfZ1aTt4mP0BHNlg32xgdp3yhcDrm/l8SYuAX5Ml2AuB24F3Avc2c36tjkuMbWQGlzL20tCt8zLQ7VVUN4kC/Dvw9xHxaHp/X5rg/qaIeK2ke2mhr7KZWVG6KA5D1gAxnaxb16nABqn8Kkn/H1mivDgXq4fkxNjMitf8aOgq2xR4rqbs2VRORNwv6eWF18rMrBndEYeJiPOA8wbfS9oSeDNZsrwb8Fmy6d7G1b1ADSfGPSbf+unW46G5pXjsdElLxf8BZ0v6DNlcnFPJWiv+D0DSG4A/llc9M7PGuiQOv0RE/J6se9vlg2Vplb2mdNSS0GbWJYpZEnqsfRSYCDwAPA/cRzbQ76Npv4CPlVM1M7NhdEccbkpuUN+w3GJsZsXqkr5tEfEI8C5Jgys1PRwRD+f231Za5czMhtIlcXgsODHuYe5W8VLuPlGQ7grIK8lWXlpZdkXMzJrWXXG4bdyVwsyK1wWP8NKk9D8ClpFNFbRU0g9b6ctmZlaaLojDY8GJsQFZS2kvt5b2+v0XrSor343S6enP15JNZv86sn9GvlpajczMmtQlcbjt3JXCzIrXHcF2d+B1EfFUev87SYcCd5ZYJzOz5nRHHG47J8ZmVqzuaoWovZOBUmphZtaK7orDbeXE2FYxXHeCTh2k524SFdMdAfk64FxJnyJbAnor4CvAtSXWycysOd0Rh9vOfYzNrHjdMejjk8B4YAnZPMZLgJcB/1ZmpczMmtIdcbjt3GJsZoUS3fEIL00Yv6ekycAU4KGIWFZytczMhtUtcXgsODG2lgzVJaHMbhbuKtFBAtShPXElfW6Y/QBExBcLqZCZ2Uh0cBwea06Mzax4ndtS8e4mjgnAibGZVVvnxuEx5cTY2sattta0Dg3IEbFL2XUwM2uLDo3DY61Sg+8kfVnS3ZJuk3SZpA1S+VaSnpO0KG1nlV1XMxs5TyxfXY7DZr3Bcbi+SiXGZNMcvT4i3gj8Djg2t+++iJietiPKqZ6ZtYVHQ1eZ47BZL3AcrqtSiXFEXBMRK9Pbm4EtyqyPmY2BZoNxDwbkKnAcNusBbYzDkvaUdI+kJZKOqbN/S0nXp6dQN0raIpXvknsCtUjSXyXtk/adLemB3L7p7bnx4VUqMa7xYeCq3PutJf1W0k8l/UNZlTKz0dNAc5uVznHYrEu1Iw5LGgd8E9gL2AY4WNI2NYedBsxJT6FOBE4BiIgbBp9AAbsCzwLX5M77bO4J1aJ23HMzCh98J+k6YNM6u46LiCvSMccBK4Hz075lwNSIeEzSW4HLJW0bEU/Xuf5MYCbAWqw9FrdgZqPUi/3WqqTIODxuwoSxuAUzG6U2xeHtgSURcT+ApLnADODO3DHbAJ9Kr28ALq9znfcBV0XEs22p1SgUnhhHxG5D7Zd0KPCPwLsiItI5K4AV6fUtku4DXg0srHP9WcAsgPU10f/8mlWR/88sVZFxePzUKf7bNqui9vyfuTnwUO59H7BDzTG3AvsDXwf2BdaTtGFEPJY75iDgqzXnnSzpeOB64JgUg8ZcpbpSSNoTOBp4b/5Xg6SNUnM9kl4BTAPuL6eWZjYq7mNcaY7DZj2gtTg8SdLC3DYzdyU1uHreZ4B3Svot8E7gYbKnUdkFstVD3wBcnTvnWOC1wNuAiWQxqRBVm8f4DGA8cG1aQermNPL5HcCJklYC/cARaTlWM+swon4ktcpwHDbrci3G4UcjYrsG+/qAKbn3WwBL8wdExFJgPwBJ6wL7R8RTuUMOBC6LiOdz5yxLL1dI+h5Zcl2ISiXGEfGqBuWXApcWXB0zGytuDa4sx2GzHtGeOLwAmCZpa7KW4IOAD+QPkDQJeDwiBshagmfXXONgVp0WEkmTI2KZsl/n+wCL21LbJlQqMTaz3uAZJ8zMytWOOBwRKyUdRdYNYhwwOyLukHQisDAi5gE7A6dICuAm4MgX6iBtRdbi/NOaS58vaSOyhu1FQGHzpjsxNrPiucXYzKxcbYrDETEfmF9Tdnzu9SXAJQ3OfZBsAF9t+a7tqV3rKjX4zsx6QJPLkI52KqFGSxvXOW7IyenNzLpOQXG4EzkxNrPiFTMrxVBLGwNNT05vZtZ9PDtQXU6MzaxwRbRUNLm08QuT00fE34DByenNzLqaW4zrc2JsZsUrvqWidmnjQfUmp39Jfzczs67jFuO6PPjOzIoVLY2GniQpv7LarLSqGjDipY3zmpmc3sysu7QWh3uKE2MzK17zqedQE8uPaGnjGsNOTm9m1pXcBFCXE2MzK5Qopt9abmnjd+aXNq4x7OT0Zmbdpqg43Incx9jMildM37YzgPXIljZeJOksAEmbSZoP2eT0wODk9HcBF0fEHaP+ZDOzqnMf47rcYmxmhVPdXg3tNcTSxkuBvXPvXzI5vZlZtysiDnciJ8ZmVqwebYUwM6sMx+GGnBibWeE8GtrMrFyOw/U5MTazwnnQh5lZuRyH63NibGbFc0A2MyuX43BdTozNrFg9usyomVllOA435MTYzIrngGxmVi7H4boqN4+xpM9LejjNO7pI0t65fcdKWiLpHkl7lFlPMxuZwYnlm9msHI7DZt3NcbixqrYYnx4Rp+ULJG1DtirVtsBmwHWSXh0R/WVU0MxGTgM9GG07j+OwWRdzHK6vci3GQ5gBzI2IFRHxALAE2L7kOplZq5pdbckxu4och826geNwQ1VNjI+SdJuk2ZImpLLNgYdyx/SlslVImilpoaSFz7OiiLqaWYs00NxmpWpLHO7/81+KqKuZtchxuL5SEmNJ10laXGebAZwJvBKYDiwDvjJ4Wp1LveS3TETMiojtImK7NRg/ZvdgZqPglorSFRWHx627zpjdg5mNQpvisKQ905iDJZKOqbN/S0nXpx/aN0raIrevPzeWYV6ufGtJv5J0r6SLJK05yrttWil9jCNit2aOk/Rt4EfpbR8wJbd7C2Bpm6tmZgXoxQEdVeM4bNbb2hGHJY0Dvgm8myw+LJA0LyLuzB12GjAnIs6RtCtwCvChtO+5iJhe59L/TTbOYa6ks4CPkP1gH3OV60ohaXLu7b7A4vR6HnCQpPGStgamAb8uun5mNkoBRDS3WSkch826XPvi8PbAkoi4PyL+BswlG4uQtw1wfXp9Q539q5AkYFfgklR0DrBP8zc3OlWcleJLkqaT/bU9CHwcICLukHQxcCewEjjSI6HNOlMv9lvrMI7DZl2uhTg8SdLC3PtZETErva437mCHmvNvBfYHvk72Q3s9SRtGxGPAWunaK4FTI+JyYEPgyYhYmbvmS8YyjJXKJcYR8aEh9p0MnFxgdcyszQbnz7Tqchw2624txuFHI2K7IS5Vq/bKnwHOkHQYcBPwMFkiDDA1IpZKegXwE0m3A083cc0xU7nE2My6nLtJmJmVq31xeNhxBxGxFNgPQNK6wP4R8VRuHxFxv6QbgTcDlwIbSFo9tRoXOpahcn2Mzaz7ecUlM7NytSkOLwCmpVkk1iRbAGhe/gBJkyQN5pvHArNT+QRJ4wePAXYC7oyIIOuL/L50zqHAFaO/4+Y4MTaz4nm6NjOzcrUhDqcW3aOAq4G7gIvTWIQTJb03HbYzcI+k3wGb8GJXrNcBCyXdSpYIn5qbzeJo4N8kLSHrc/zd0d5us9yVwswK59ZgM7NytSsOR8R8YH5N2fG515fw4gwT+WN+AbyhwTXvp6RVNZ0Ym1mxAuh3ZmxmVhrH4YacGJtZ4dxibGZWLsfh+pwYm1nxPCuFmVm5HIfrcmJsZoVzS4WZWbkch+tzYmxmxfKME2Zm5XIcbsiJsZkVKltxyRHZzKwsjsONOTE2s8LJo6HNzErlOFyfE2MzK5Yf4ZmZlctxuCEnxmZWsPBoaDOzUjkON+LE2MwKV8RoaElfBt4D/A24Dzg8Ip6sc9yDwDNAP7AyIrYb+9qZmZXLs1LUt1rZFTCzHhTR3DY61wKvj4g3Ar8Djh3i2F0iYrqTYjPrGcXE4Y7jFmMzK1aABgr4mIhrcm9vBt439p9qZtYBCorDncgtxmZWvIFoboNJkhbmtpkj/MQPA1c12BfANZJuGcX1zcw6S/NxuKdUqsVY0kXAa9LbDYAnI2K6pK2Au4B70r6bI+KI4mtoZu3QwvyZjw7VvUHSdcCmdXYdFxFXpGOOA1YC5ze4zE4RsVTSxsC1ku6OiJuarWC3cRw26w2ex7i+SiXGEfH+wdeSvgI8ldt9X0RML75WZtZ2bQrIEbHbUPslHQr8I/CuiPofGhFL05/LJV0GbA/0bGLsOGzWI5wY11WpxHiQJAEHAruWXRcza7MACujbJmlP4GjgnRHxbINj1gFWi4hn0uvdgRPHvnbV5zhs1sUKisOdqKp9jP8BeCQi7s2VbS3pt5J+KukfGp0oaeZgf8TnWTH2NTWzlohA0dw2SmcA65F1j1gk6SwASZtJmp+O2QT4maRbgV8DV0bEj0f7wV2iLXG4/89/GfuamllLCozDHafwFuNm+gQCBwMX5vYtA6ZGxGOS3gpcLmnbiHi69iIRMQuYBbC+Jvbe36hZJygg2EbEqxqULwX2Tq/vB9405pWpmCLj8PipUxyHzaqoB5PeZhSeGDfRJ3B1YD/grblzVkDW/BsRt0i6D3g1sHAMq2pmYyGAfgfkMjkOm/W4Nsbh1G3t68A44DsRcWrN/i2B2cBGwOPAByOiT9J04ExgfbIFlk6OiIvSOWcD7+TFMQ6HRcSitlR4GFXsY7wbcHdE9A0WSNoIeDwi+iW9ApgG3F9WBc1sdHrx8VyHcRw263LtiMOSxgHfBN4N9AELJM2LiDtzh50GzImIcyTtCpwCfAh4FjgkIu6VtBlwi6SrcyuUfjYiLhl1JVtUxcT4IFZ9fAfwDuBESSvJflUcERGPF14zM2sPJ8ZV5zhs1u3aE4e3B5akbmlImgvMAPKJ8TbAp9LrG4DLs4+P371YlVgqaTlZq/KTlKhyiXFEHFan7FLg0uJrY2bt15vLjHYSx2GzbtdSHJ4kKd9lalYaRwCwOfBQbl8fsEPN+bcC+5N1t9gXWE/ShhHx2OABkrYH1gTuy513sqTjgeuBY1J3rjFXucTYzLpc4MTYzKxMrcXhoRZaUoOr530GOEPSYWRzxD9MtuhSdgFpMnAucGhEDE4idyzwR7JkeRbZ1JuFTKXpxNjMiuf5M83MytWeONwHTMm93wJYmj8gzQS0H4CkdYH9I+Kp9H594ErgPyLi5tw5y9LLFZK+R5ZcF8KJsZkVTgPOjM3MytSmOLwAmCZpa7KW4IOAD6zyOdIksoG7A2QtwbNT+ZrAZWQD875fc87kiFiWFhraB1jcjso2w4mxmRUrgAF3pTAzK02b4nBErJR0FHA12XRtsyPiDkknAgsjYh6wM3CKpCDrSnFkOv1AskG9G6ZuFvDitGznp5lwBCwCjhh1ZZvkxNjMCubBd2Zm5WpfHI6I+cD8mrLjc68vAV4y7VpEnAec1+CapS1F78TYzIrnxNjMrFyOw3U5MTaz4jkgm5mVy3G4LifGZlYs9zE2MyuX43BDTozNrGABA/1lV8LMrIc5DjfixNjMiuWWCjOzcjkON+TE2MyK575tZmblchyuy4mxmRXPAdnMrFyOw3U5MTazgnkeYzOzcjkON+LE2MyKFYCXhDYzK4/jcENOjM2seA7IZmblchyuy4mxmRUsPBrazKxUjsONrFbGh0o6QNIdkgYkbVez71hJSyTdI2mPXPmeqWyJpGOKr7WZtUVAxEBTm40dx2GzHuY43FBZLcaLgf2Ab+ULJW0DHARsC2wGXCfp1Wn3N4F3A33AAknzIuLO4qpsZm3jlooqcBw262WOw3WVkhhHxF0Akmp3zQDmRsQK4AFJS4Dt074lEXF/Om9uOtYB2awTeTR06RyHzXqc43BdVetjvDlwc+59XyoDeKimfId6F5A0E5gJsBZrj0EVzWxUIjzoo9raGofHTZgwBlU0s1FxHG5ozBJjSdcBm9bZdVxEXNHotDplQf2+0HV/6kTELGAWwPqa6J9DZhUU/f1lV6EnVCEOj586xXHYrIIch+sbs8Q4InYbwWl9wJTc+y2Apel1o3Iz6yieWL4ojsNmVp/jcCOlzEoxhHnAQZLGS9oamAb8GlgATJO0taQ1yQaGzCuxnmY2UkE26KOZzcrgOGzW7RyHGyqlj7GkfYH/ATYCrpS0KCL2iIg7JF1MNphjJXBkRPSnc44CrgbGAbMj4o4y6m5mbdCDUwBVjeOwWY9zHK6rlBbjiLgsIraIiPERsUlE7JHbd3JEvDIiXhMRV+XK50fEq9O+k8uot5mNXgAxEE1toyHpJEm3SVok6RpJmzU47lBJ96bt0FF9aAdxHDbrXe2Mw8PNby5pS0nXp3h8o6Qtcvvqxl9Jb5V0e7rmN1Rn+pyxUrWuFGbW7SKylopmttH5ckS8MSKmAz8Cjq89QNJE4ASy2RW2B06Q5GkUzKy7tSkOSxpHNr/5XsA2wMFpLvS804A5EfFG4ETglHTuUPH3TLKZbaalbc923HYznBibWeGiv7+pbVSfEfF07u061J9BYQ/g2oh4PCKeAK6lwABsZlaWNsXh7Unzm0fE34DB+c3ztgGuT69vyO2vG38lTQbWj4hfRkQAc4B9Rn/HzanaPMZt9QxPPHpdXPL7OrsmAY8WXZ8x4nuprm66n3r3suVILvQMT1x9XVwyqcnD15K0MPd+VpoKrCmSTgYOAZ4CdqlzyOa8dG7ezescZyP0t4f6Hn3gk592HO4cvpfqqmIcrhdDa+c3vxXYH/g6sC+wnqQNG5y7edr66pQXoqsT44jYqF65pIURsV3R9RkLvpfq6qb7aee9RETbWmSHm6c3Io4DjpN0LHAU2WO7VS5Rr4rtqp85Dnca30t1VTQONxNDPwOcIekw4CbgYbKBvY3OLTUud3VibGbdrYV5ei8AruSliXEfsHPu/RbAjaOumJlZbxhq3nMAImIpsB+ApHWB/SPiKUmN4m9fet3wmmPJfYzNrCtJmpZ7+17g7jqHXQ3sLmlCGvSxeyozM7PhDTu/uaRJkgbzzWOB2el13fgbEcuAZyTtmGajOARotFJn2/VqYtx0H8UO4Huprm66n068l1MlLZZ0G1nA/QSApO0kfQcgIh4HTiIL7guAE1OZjb1O/G+qEd9LNXXTvUAF7yciVpJ1U7sauAu4OM2FfqKk96bDdgbukfQ7YBPg5HTuUPH3n4HvAEuA+4AXpo0cawovCWhmZmZm1rMtxmZmZmZmq3BibGZmZmZGlyfGkg6QdIekAUnb1ew7Ni01eI+kPXLlQy5tWBWSPi/p4bTc7SJJe+f21b23KuuU770RSQ+m5SsXDc73KGmipGvTUpfXVnVFNUmzJS2XtDhXVrfuynwj/T3dJukt5dXcOoHjsONwURyHrR26OjEGFpNNEXJTvlDZcoUHAduSrXL1v5LGqbmlDavk9IiYnrb50PjeyqzkcDrwe29kl/R3MfiP/zHA9RExjWzVn6r+Q3M2L13trVHd9+LFJTpnki3baTYUx2HH4SI5DtuodHViHBF3RcQ9dXbNAOZGxIqIeIBs1OP2NLe0YdU1urcq64bvvZ4ZwDnp9TkUuKRlKyLiJqB2JoZGdZ9BtuZ9RMTNwAbKlu80q8tx2HG4ZI7D1pKuToyHMNQyhJ20POxR6THK7NzjoU67B+jMOtcK4BpJt0iamco2SfMxkv7cuLTata5R3bvh78qqwXG4WjqxzrUch23UOn7lOw2zJGyj0+qUBfV/KJQ2n91Q90b26OQksvqdBHwF+DCducRtJ9a51k4RsVTSxsC1kuotJtENuuHvytrMcdhxuCIch23UOj4xbmFJ2LyhljAccmnDIjV7b5K+DfwovR12ecYK6sQ6ryIteUlELJd0GdljyUckTY6IZekx1/JSK9maRnXv+L8raz/HYcfhKnActnbo1a4U84CDJI2XtDVZB/Zf08TShlVR059oX7IBLtD43qqsY773eiStI2m9wddkq6wtJruHQ9Nhh1LgkpZt0Kju84BD0qjoHYGnBh/1mbXIcbhaOuZ7r8dx2HG4XTq+xXgokvYF/gfYCLhS0qKI2CMtV3gxcCewEjgyIvrTOYNLG44DZkfEHSVVfzhfkjSd7PHJg8DHAYa6t6qKiJUd9L3XswlwmSTI/p+6ICJ+LGkBcLGkjwB/AA4osY4NSbqQbMnOSZL6gBOAU6lf9/nA3mSDiZ4FDi+8wtZRHIcdhwviOGxt4SWhzczMzMzo3a4UZmZmZmarcGJsZmZmZoYTYzMzMzMzwImxmZmZmRngxNjMzMzMDHBibGZmZmYGODE2MzMzMwOcGFsFSHqbpNskrZVWL7pD0uvLrpeZWa9wHDbLeIEPqwRJXwDWAl4G9EXEKSVXycyspzgOmzkxtoqQtCawAPgr8PaqL59qZtZtHIfN3JXCqmMisC6wHlmLhZmZFctx2HqeW4ytEiTNA+YCWwOTI+KokqtkZtZTHIfNYPWyK2Am6RBgZURcIGkc8AtJu0bET8qum5lZL3AcNsu4xdjMzMzMDPcxNjMzMzMDnBibmZmZmQFOjM3MzMzMACfGZmZmZmaAE2MzMzMzM8CJsZmZmZkZ4MTYzMzMzAyA/x/NZR/UR7AwFwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "xlim = np.r_[-100., 100.]\n", "ylim = np.r_[-100., 100]\n", "\n", "# Plot a cross section of the conductivity model\n", "fig, ax = plt.subplots(1,2, figsize=(10, 4))\n", "\n", "# conductivity model\n", "cb = plt.colorbar(mesh.plotImage(np.log10(sigma), ax=ax[0], mirror=True)[0], ax=ax[0])\n", "\n", "# plot formatting and titles\n", "cb.set_label('$\\log_{10}\\sigma$', fontsize=13)\n", "ax[0].axis('equal')\n", "ax[0].set_title('Conductivity Model')\n", "\n", "# permeability model\n", "cb = plt.colorbar(mesh.plotImage(mur, ax=ax[1], mirror=True)[0], ax=ax[1])\n", "\n", "# plot formatting and titles\n", "cb.set_label('$\\mu_r$', fontsize=13)\n", "ax[1].axis('equal')\n", "ax[1].set_title('Permeability Model')\n", "\n", "[a.set_xlim(xlim) for a in ax]\n", "[a.set_ylim(ylim) for a in ax]\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## set up the source" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# find the indices where we want to put the solenoid - it is on edges \n", "# (electric field is defined on cell edges for EB formulation)\n", "solenoid_x = 500 # radius of the solenoid\n", "solenoid_zlims = np.r_[-2000., 2000.] # 1km long solenoid\n", "\n", "solenoid_z = mesh.vectorNz[(mesh.vectorNz>=solenoid_zlims[0]) & (mesh.vectorNz<=solenoid_zlims[1])]" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "solenoid_inds = Utils.closestPoints(mesh, Utils.ndgrid([np.r_[solenoid_x], np.r_[0], solenoid_z]), 'Ey')\n", "src_vec = np.zeros(mesh.nEy)\n", "src_vec[solenoid_inds] = 1." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-3000, 3000)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAF3CAYAAACIZdGbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8lOWd///XlTMCERFEDm0JWyge6ooLouXQVK0FH92l26+7te63atfqul/11/rdfRS17tJW3a1+d9e2q6tlq10P7arr1gWxIirEiAIKhSAJ4UxMCJMDIZmEkMNkrt8fc0+YzMw95DCZuYe8n49HHjP3dV/3dX/uyfDJxTX3dY2x1iIiIt6Tle4AREQkPiVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj0pbgjbGFBhjPjTGlBljyo0xP3LKi4wxW4wx+4wxLxlj8pzyfGd7v7N/ekRb9znle4wxX0nPFYmIJFc6e9CdwFXW2j8ELgWWGGOuAB4BHrPWzgSOA7c69W8FjltrPws85tTDGHMhcANwEbAE+DdjTHZKr0REZBikLUHbkDZnM9f5scBVwCtO+bPA15zny5xtnP1XG2OMU/6itbbTWnsI2A9cnoJLEBEZVmkdgzbGZBtjdgD1wFvAAaDZWhtwqtQAU53nU4FqAGd/C3BuZHmcY0REMlZOOk9ure0BLjXGjANeBS6IV815NC773MpjGGNuB24HyBpV+Ec5Z5834Ji9JMdAQY7hZMAyZUwW1a1Bphdm0dRhyTZwdr7hRLelvRsmnhV6mbqC0NAeZOqYvn+bj54Ick5BFgVxBod8J4KcnW8YlRPvpQ451mHJzYLCPPc6wWCQrKzQeevbLaNzYXSue303nT1wrCPIlNHJ6V/0WKhtC/KpsfHbi4x7qKr8QT5dmBX3TTsc3GI/7A+9V7wqOu4jbUEmnpVFnndD7rV3795Ga+3EZLSV1gQdZq1tNsaUAFcA44wxOU4veRpQ61SrAT4F1BhjcoCzgaaI8rDIY6LPsxJYCZA/eaadfPNPh+FqUmfquFHMLxpP6b4G1ty9iAWPrGfPP1zHQ2sqmFRYwG2LZ7C6rJZ15T4ev/EyAHYf9XPPSztY+73Ffdq6/skPWL50NvOmj485z7ee3sJti2aweJb7e27Fql0UTRjNLQuKXOuUlJRQXFwMwB3Pb2PZpVNY+vnJA77uHdXNrFi1i1V3LRzwsfH4WjpY9sRGttx/Tdz9kXEP1awH3mDnimspyE3NxyRusRfd9zq7H76OrKxU/akYmOi4l/y0lMe+cSkXTC5MX1D9ZIypSlZb6byLY6LTc8YYMwq4BtgNbACud6rdDKxynq92tnH2r7ehlZ5WAzc4d3kUATOBD1NzFSIiwyedPejJwLPOHRdZwMvW2jXGmArgRWPMQ8B24Gmn/tPA88aY/YR6zjcAWGvLjTEvAxVAALjTGToREcloaUvQ1tqdwJw45QeJcxeGtbYD+DOXth4GHk52jCIi6ZQBQ+4iIiOTErSIiEcpQYuIeJQStIiIRylBi4h4lBK0iIhHKUGLiHiUErSIiEcpQYuIeJQStIiIRylBi4h4lBK0iIhHKUGLiHiUErSIiEcpQYuIeJQStIiIRylBi4h4lBK0iIhHeeJbvWVwjjSfZP2eeprbu1lfWU9P0LK+so71e+qZNLaAPzhvNBsq69lQWc/6yjoAdh9tpdLX2rsdtrXqOOsr62nt6I45z3v7GvmDiWMIBIOusazfU0/RsTF8+tyzXOvsrA8QdM67fk89YwtyyM8deB9hR3ULZTUtMdcwWL6WTur8na7tRcY9VF2BIBsq6wd13YPhFru1sL6yniyPdtGi4w69Z+s52nIyjVGlnhJ0hmtuDyXU5zeHvun9+U1VHGw4wcGGExTkZrFhT0NvOYTe6JHbkV7YVEXlUX/c87ywuYqqYydc46huOkl100myjXusx5oClJ0InbcrEOS/ttXQ2NZ5miuMVVbTAsS/hsHw+TsTthcZdzI8v7mK/JzUZMZEsT+/uYqsBL+vdIoX9wubq5h9/tg0RZQeStAZbOq4UcwvGk/pvgZ+dcs8Fjyynl99+3IeWlPBpMICbls8g9Vltawr9/H4jZcBsPuon3te2sGvvt33e3mvf/IDli+dzbzp42PO862nt3DbohksnjXRNZYVq3ZRNGE0tywocq1TUlJCcXHovHc8v41ll05h6ecnD/i6d1Q3s2LVrphrGCxfSwfLntjo2l5k3EM164E3eOaWeRTkZielvdNxi73ovtf51S3zyPJoho6Oe8lPS3nsG5dyweTCNEbVP//xl8lry6P/wRERESVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8aicdAcgg3ek+SS/3X4EgKfePUBP0PJkyQF+ufEQ547OIxC0vFZWS8VRPxdNOQBApc9Ppa+VJ0sO9Glra9Vxnio5wNbpx2PO896+RnqClvJav2ssz26q4vzCAk52B13rHDzYxW5C511b7uNoy0kOH2sf8HWXVTdTVtMScw2DVefvoM7f6dpeZNxD1RUI8tS7B8jPyU5Ke6fjFru18OS7B8gyJiVxDFR03JW+Vp569wCzzy9MY1SppwR9hmhu7wo9ngw/dtN8souWk919y9v7bvdpwzkmfvvu+04d35WwzonuvudNdL7TnSfycaiiX7to0XEP/Xzd5Of2JK29RBLF3nKyG4/m57hx9+c9eKZJW4I2xnwKeA44HwgCK621PzPGjAdeAqYDh4E/t9YeN8YY4GfAdUA7cIu19vdOWzcDDzhNP2StfTaV15IuU8eNYn7ReEr3NXDv0gt4bedR7lt6AT09lkmFBdy2eAYXTTmbdeU+7lt6AQC7j/qp83f0bodtO3yc5UtnM2/6+JjzVNT6uW3RDBbPmugaS0dXD0UTRnPLgiLXOiUldRQXh85b1djOskunsPTzkwd83Tuqm1mxalfMNQyWr6WDTQePubYXGfdQ/er9w9y7dDYFuanpQbvFvrL0IPcumU1WljczdHTc7+5p4N6ls7lgsvd70Pcnsa10jkEHgL+x1l4AXAHcaYy5ELgXeMdaOxN4x9kGWArMdH5uB54EcBL6CmA+cDmwwhhzTiovRERkOKQtQVtrj4Z7wNbaVmA3MBVYBoR7wM8CX3OeLwOesyGbgXHGmMnAV4C3rLVN1trjwFvAkhReiojIsPDEGLQxZjowB9gCTLLWHoVQEjfGnOdUmwpURxxW45S5lcc7z+2Eet/knf/Z5F1AmnR0dOCrq6OrK8CmTR9graWkpITqmk5a87MoCX5CxdEA9fUBSkpKAKhuDdLW1tG7HdbScpLt27dz4nDsf72bmk6yc2crwVr3t0vNkU4CzVmUdFe51mlra+s9b0NjB7vKjzPq2J4BX/fB5h78rV0x1zBYxzuCdHa6txcZ91AFg0FKS0vJy07N0IJr7BZK3i3x7IeE0XG3tbWzdetW6saOrBvP0p6gjTFjgP8Gvmet9Rv3N0y8HTZBeWyhtSuBlQD5k2fGrZNJCgoKOH/SePb6G7jyyi9gStdTXFzMxrYKJhUWULx4Bv6yWmqsj+Liy4DQGPSvD+yguHhxn7Ye3/0Bc+bEH4N++sAWLrkk8Rj0hpZdFE0YTXHCMegSiouLAXixehsXXzSF4kGMQY+rbmZVzS6KixcO+Nh4fC0d5P9+Y29s0SLjHqqst99g8eLFKRyDdon9zdcp/mKxh8eg+8Y9Zkcpc+demhFj0MmU1j9HxphcQsn519ba3zrFdc7QBc5jvVNeA3wq4vBpQG2CchGRjJa2BO3clfE0sNta+y8Ru1YDNzvPbwZWRZTfZEKuAFqcoZA3gWuNMec4Hw5e65SJiGS0dA5xLAC+BXxsjNnhlN0P/AR42RhzK/AJ8GfOvt8RusVuP6Hb7L4NYK1tMsY8CHzk1PuxtbYpNZcgIjJ80pagrbUbiT9+DHB1nPoWuNOlrWeAZ5IXnYhI+o2sj0RFRDKIErSIiEcpQYuIeJQStIiIRylBi4h4lBK0iIhHKUGLiHhU2tfikP75k/INfL/0Oab4G6ktnMCji29i24Lr0h2WiAwj9aAzwJ+Ub+Anax9nmr+BLCzT/A38ZO3jfHn72+kOTUSGkXrQGeD7pc9xVqCzT9lZgU6+s/aXLCxaAMD9r35MT9By328/5j8//ASAg40nWFNWS2tngLEFHwOnvpPwvt9+3Ke9rVXHuf+3HzM3zmp27+1r5FDjCRbNdF/NLnzOPXVtrnWO1nbyZlPovGvLfWw5dIzSfY2nu/wYZdXNVBz1x1zDYIW/k9Ctvci4h6orEOT+Vz9O2XcSusVubeg9k2D1yLSKjrvS18r9r36s7yQU75nij5/EIssvnlLI+sp6Lp5ayFl52UwqLODiqYUcamxj88EmLp4aemPnZBm2f9Lcux3p4qlnxy0Pte++DyA321A0YXTCOnltPmZF7E90vkR6gkEqjvoHdWw8E8fk9b528UTHPVQXTzmb/NzU/Oc1UewXTT0bj642Gjfui6eczezJY9MUUXooQWeA2sIJTPM3xJTXj5vI1+dMpXRfAzfO/wxPlBzgL+Z/hkMNJ5hUWMBfzP8MYwtymTAmn7+Y/xkgtB70R4eberfDXv39EW6c/+m460Gv3eXjxvmfTrge9F5fK0UTRse0G6nk5CGKnf3v7W0c9HcSXjTlbPb4WhOeayB8LR28tLXatb3IuIfqR69VcOP8T6duPWiX2B/4n138xeWf9u560FFxP7+pihvnfzoj1oP+30lsS2PQGeDRxTfRnpPfp6w9J5+nrr01TRGJSCqoB50BVl/0JYDYuzjmXMP8NMcmIsNHCTpDrL7oS72JOizuFy+KyBlDQxwiIh6lBC0i4lFK0CIiHqUELSLiUUrQIiIepQQtIuJRStAiIh6l+6AzhJYbFRl51IPOAFpuVGRkUg86AyRabnTJ7C/S1hlgXYWPnqDlzXIf6yrqmFSYz6fPPYt1zvab5T4AKo+2Uulr7d0O21p1nHXlPppOdMWc/719jXzm3LM42d3jGuO6ijqKJoxm8rhRrnV21QXodM67rsLHqLzsQS3WU1bdTFlNS8w1DFZ4uVG39iLjHqquQJB1FXXk56Smb+QWu7Wh34FXlxuNjrvS18q68jo+aWpPY1SppwSdARItN9rWGQDglW01vY+fNLXzSVM7486q4a2Kuj77K33+PtuRXtlWw+Fj8f8BvLKthjp/Z9x9AEdbOjja0sHofPe3VGNjgMrO0HmDFl7dfqQ3/oEoq27ujSkZ6vwdCduLjDsZXtlWk7IEnSj2V7bVeDZBx4v7ld9Xaz1o8Z7+LDe68ltzWfDIev79prk8tKaCSYUF3LZ4BqvLallX7uPxGy8DQsuN3vPSDv79prl92rr+yQ9YvnR23OVGv/X0Fm5bNCPhcqMrVu2iaMJobllQ5FqnpKSE4uLQee94ftuglxvdUd3MilW7Yq5hsHwtHSx7YqNre5FxD9WsB95g5bf+KHXLjbrEXnTf66z81lzvLjcaFfeSn5by2DcuzYjlRn95c/La0hh0BtByoyIjk3rQGUDLjYqMTErQGULLjYqMPBriEBHxKCVoERGPUoIWEfEoJWgREY/Sh4QZQmtxiIw86kFnAK3FITIyqQedARKtxbGwaAEAD66poCdo+fFrFTzz/iEgNP36tZ21NLR2ct7YCiA01bvS18qPX6vo097WquM8uKaCuZ+JnUn43r5GfC0dLJrpPpPw2U1VZGcZPmk66VqnpqaT0tbQedeW+9hZ08xHh4/34xXoq6wmtBZH9DUMVngtDrf2IuMeqq5AkAfXVJCfk5qZhG6xWws/XlNBlkenekfHXelr5cE1FZrqLd6TaC2O3ufjCnofjYFJYwuYMq6AKeNG0dDa2bvf39Hdp36f9s4eFbc8VN99X2QMieqcPJbVZ39/2oynoW0U26qOD+rYeLKzwvHEby867qGaMm5UytbiSBT71HGj8Gh+jhv3YN8vmUwJOgP0Zy2OWxfO4Jn3D/OdRTPwtXQwqbCA7yyawXmFBawr9/GdRTOA0Focu4609G6Hrd3l49ZFRXHX4nh3bwO3LixKuBZHdVP76dfi6PmEYue8Ww8fH9JaHJ8cOxFzDYPla+lgdVmta3uRcQ/Vo2/u4daFRalbi8Ml9od/t5tbFxZ5dy2OqLhf2VbDrQuLMmItjtuS2JbGoDOA1uIQGZnUg84AWotDZGRSgs4QWotDZOTREIeIiEepB50hNFFFZORRDzoDaKKKyMiU1h60MeYZ4KtAvbX2YqdsPPASMB04DPy5tfa4CX152s+A64B24BZr7e+dY24GHnCafcha+2wqr2O49Weiyj+t20NP0PJPb+7hlxsPkZ+TRcvJbl7bWUvVsXamn7sHODVR5Z/e3NOnva1Vx/mnN/fEvc3uvX2N+E928+GhJtcYn91Uxdj8HBrbYr90NqyqqoutnaHzri33sa++lfJaf/9ehAjhiSrR1zBYPmeiilt7kXEPVVcgyD+v25OyiSpusVsbes94daJKdNyVvlb+ed0eTVRJsf8AHgeeiyi7F3jHWvsTY8y9zvZyYCkw0/mZDzwJzHcS+gpgLmCBbcaY1dbagU9R86j+TFQJT3yIfAz/9C3P7rMdKT8323UCRX6O+75Tx2clrJOb3fe8/Wkzfix9r2moTtdedNxDP9/grnswEsWen5ONR2+Djht3Kl83r0hrgrbWlhpjpkcVLwOKnefPAiWEEvQy4DlrrQU2G2PGGWMmO3XfstY2ARhj3gKWAP85zOGnTH8mqtx91Uxe/Kiau6+eScvJ7t4vjf3MhNGsK/dx99UzgdBElQMNbb3bYe/ubeDuqz4btwf94eGm035pbGNbZz++NPYIxcWh85bX+oc0UaWhtTPmGgbL19LB27vrXNuLjHuo/nXDfu666rMp/NLY+LH/y9t7ufuqz3p3okpU3K9/fJS7rvpsRkxU+f+S2Fa6e9DxTLLWHgWw1h41xpznlE8FqiPq1ThlbuUxjDG3A7cD5J3/2SSHPXweXXwTP1n7eJ9hjvacfP71qpvw1dXR1RVg06YPsNZSUlJCdU0nrflZlAQ/oeJogPr6ACUlJQBUtwZpa+vo3Q5raTnJ9u3bOXE4NnE0NZ1k585WgrXub5eaI50EmrMo6a5yrdPW1tZ73obGDnaVH2fUsYEPHRxs7sHf2hVzDYN1vCNIZ6d7e5FxD1UwGKS0tJS87NQkRtfYLZS8W+LZIY7ouNva2tm6dSt1Y9WD9qp47ySboDy20NqVwEqA/Mkz49bxIteJKvOWMn/SePb6G7jyyi9gStdTXFzMxrYKJhUWULx4Bv6yWmqsj+Liy4BQD/rXB3ZQXLy4zzke3/0Bc+bMjtuDfvrAFi65JHEPekPLLoomjKY4YQ+6hOLiYgBerN7GxRdNoXgQPehx1c2sqtlFcfHCAR8bj6+lg/zfb+yNLVpk3EOV9fYbLF68OIU9aJfY33yd4i8We7gH3TfuMTtKmTv30ozoQSeTF/8c1TlDFziP9U55DfCpiHrTgNoE5SIiGc2LCXo1cLPz/GZgVUT5TSbkCqDFGQp5E7jWGHOOMeYc4Fqn7Iyh2+xERqZ032b3n4Q+5JtgjKkhdDfGT4CXjTG3Ap8Af+ZU/x2hW+z2E7rN7tsA1tomY8yDwEdOvR+HPzA8U2g96L60HnT/aT3ozJbuuzi+6bLr6jh1LXCnSzvPAM8kMTRP0XrQfWk96P7TetCZLZM+JByxtB50X1oPuv+0HnTqaT3oEebRxTfRld33b2lXdo7WgxY5wylBZwgbtAm3ReTMoyGODPD90ufItz19yvJtD99Z+0tuvvxaGtu6KKtppidoKatuZmdNC+cVdoSeO9tl1c3AqQ8Jw9th4Tp52bF/s0NtNnP2qFzXGMtqWmjtCMS0G+lgSw/nOPt31jQzfcJopowb1e/X4VSsoQ8JE51rIMJrcbi1Fxn3UHUFguysaUnZGLRb7NaGPmz16oeE0XFX+lrZWdNMVyCYxqhSz4Q+ext58ifPtJNv/mm6w+iXg4/8MVlx5t4EMcxY/hoA543Np761k0umnc3OmhaAmOcAlUdb6eoJ9m6HRdfr776B1Gn1tzK2cGy/6w/lXAPha+nofe3iiYx7qHbWtPC5SWPJz01NgnaLfWdNCxdPLfRsgo6Oe2dNC3nZWcyenJzfw3B67e5F26y1c5PRlnrQGaA/HxKuvmshCx5Zz+q7FvLQmoretThWl9WyrtzH4zeemkl4z0s7WH1X31l41z/5AcuXxp9J+K2nt5x2LY4Vq3b1Yy2Okt7Zf3c8v21IHxKuWLWLVXclbybhsic2xrwmYZFxD9WsB95g1V0LUjyTMDb2ovteZ/WdC737IWFU3Et+Wspj38iMmYTm7uS1pTHoDPDOjHkx/WcLvP+5K9IRjoikiHrQGeDqgx/FLDhigCt2b+Z/2jppbOuisa2TnqClobWTxrZOsrNM6Lmz3dAamujS2NZ3O6yxLVQ3uhzobTPevlPHdzG2IDdhnZZOe9o4+iN0TV2DOjZue22J24uMe6i6AkEa2zpTNlHFLXZrQ9dtPDrEER33UN4vmUxj0BmgP2PQYRPG5PUumh/9HIjZDnMrP92+gdTp6uoiLy9xHP0xlGMH015k3Mk419j8nJSNQbvF3tjWxfjReZ5dDzo67mT/zofTtr+7VmPQI0l/xqDX3L2IBY+sZ+sDX+7XGPTa7/VdzS51Y9DFgDfHoLfcf03c/clczW7WA2/w0QPXpH01u6L7XmfrD67x+Bh0ce92Ro1B/13y2tIYdAbQGLTIyKQhjgyw8clvx+1B1xROZOFf/6pP2bizcmlu7477HIjZDnMrP92+gdTp7u4mNzdxHP0xlGMH015k3Mk416jc7JQNcbjF3tzeTWFBjmd70NFxJ/t3PpzKVnxFQxwjydQ4yTlc/qXPTaR0XyNr7l7I0p+9x4a/Keah13czqTCf2xbN4LWdtawrr+NfvzkHgN0+Z4jju1FDHE99wPIlLkMczzhDHAlWs1uxujw0xPGF6a513n//fRYsCK2+d8cL21h26VSWXnz+6S4/xo6aZlasKmfVnQsGfGw8Pn8Hy554nw1/Uxx3f2TcQzXnwbfY8LfFKZuo4hZ7OA6v3gcdHfeSnzlDHBmwmt34FclrSwk6A/SYLHJs7AyqHpPFhj2h5P385qrex//+fQ1jC3IoyM3mtbJa9tW39e6v9Pmp83f2bocdaDjB85urqPS1xpxn1xE/L2yuoupYu2uMq8tqOXd0HtkJemSHPunmYHbovFsONeHvCNDYNvBP5cuqm/mkqT3mGgbL5++gKxB0bS8y7mR4fvPhlN3FkSj25zdXeTZBR8dd5+/khc1VWm5UvCc7TnKOLu92psB295x6DP9ElncFbJ/tSJH1o3X1WNd9p+q4Hw8QsH3Pm+h8pztP+PhkiH7tokXHPeTz9ViyTGqmLCeKvbsn6NkEHS/ursDp34NnGiXoDBA0huw4nxUEI/5x+fwdoceW0GNHdxBfS0dMeV3UdqQ6f2fccoC6lg7XfWGtHYGEdZo7bJ/9/WkzbiwJrmEwol+jaNFxD/l8LR0pG+JIFLuvpdOzt9nFi7vO38E5GTAGnUxK0Bkgy+WD3MjywoLQG3dswak38NiCXMYW5NLR3dlbPrYgJ6beqfo5cctPty8sO8skrNOSEx3f6duMH0vstQ5Fe1dPwvai4x6qsQU5qZuokiD2sQU5nu1Bx4t7sO+XTKa7ODLAoUe+6vrV5UXL1/QpO7+woLdHGP0ciNkOcys/3b6B1Ons7CQ/P7/f9YdyrmS2Fxl3Ms51zlm5KUvQbrH7/B2cNzbfswk6Ou5k/86H05YfXKO7OCS+8Bs50fN426crP92+ftfp7F8c/TGUYwfcXmfyznW8vRvoTlp7p+USe73Xp03HiTvZv3Ov00QVERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfGoMyZBG2OWGGP2GGP2G2PuTXc8IiJDdUYkaGNMNvAEsBS4EPimMebC9EYlIjI0Z0SCBi4H9ltrD1pru4AXgWVpjklEZEj6laCNMe8YY66LKls5PCENylSgOmK7xikTEclY/e1BFwHLjTErIsrmDkM8g2XilNmYSsbcbozZaozZmoKYRESGpL8Juhm4GphkjHnNGHP2MMY0GDXApyK2pwG10ZWstSuttXOttV764yIiEld/E7Sx1gastf8H+G9gI3De8IU1YB8BM40xRcaYPOAGYHWaYxIRGZKcftZ7KvzEWvsfxpiPgTuHJ6SBs9YGjDF3AW8C2cAz1tryNIclIjIk/UrQ1tpfRG1vA/5yWCIaJGvt74DfpTsOEZFkOVNusxMROeMoQYuIeJQStIiIRylBi4h4lBK0iIhHKUGLiHiUErSIiEcpQYuIeJQStIiIRylBi4h4lBK0iIhHKUGLiHiUErSIiEcpQYuIeJQStIiIRylBi4h4lBK0iIhHKUGLiHiUErSIiEcpQYuIeJQStIiIRylBi4h4lBK0iIhHKUGLiHiUErSIiEflpDsASa5JhfnU+TvjPgditsPcyk+3byB1Oju7yM/P63f9oZwrme1Fxp2Mc407K5f8nNT0jdxir/N3MnFsPlkmJWEMWHTcyf6dD6eqJLalBH2GWDxrIqV7G/jirIm8vLUGgC/OmshrZUc52d3DF2dNBKDS18rOmpbe7bCXt9YwY8Jo5k4/J6btl7fWMPnsAhbNnOB6/shzujl61MfkyRN7648tyElY301ZdQt76loHdWw8Pn9n72sXT2TcQ/Xy1hq+OGtiyhK0W+zhOLyaoKPjfnlrDZdMO5vZ549NY1RmHWeeAAAeCklEQVT982ES2zLW2iQ2lznyJ8+0k2/+abrD6JeDj3w17lhUEJixfE2fsjH5ObR1BuI+B2K2w9zKT7dvIHUCgQA5OYnj6I+hHDuY9iLjTsa5crMN+TnZSWnvdNxib+sMcFZeNlnGmxk6Ou5k/86HU/mPl2yz1s5NRlvev1rB7Z+QAS6fPp4PDzfxyh1Xcv1Tm3jju4t4+PXdTCrM5zuLZvDazlrWldfxr9+cA4R60Pe8tIM3vruoT1t/9tQmli/9HHM/Mz7mPDc98yG3LZqRsAf9w9XlFE0Yzc1fmO5aZ/PmzVxxxRUA/PWvt7HsD6ey5OLzE157PGU1zaxYVc7/3LlgwMfGU+fv6H3t4omMe6gWPbqBN767OGU9aLfYFz26gbXfXYxH83NM3Et/9h6PfePSjOhBf/rHyWtLPegMcOiRr8ZN0hYoiupBnz0ql5aT3XGfAzHbYW7lp9s3kDqBQDc5OYnj6I+hHDuY9iLjTsa5CnKzUtiDjh97y8luxhbkeLgH3TfuZP/Oh9POH35FPeiR5ERuPmO6O2PK23ML+PqcqZTua2DN3YtY8Mh6ylZcy0NrKphUWMBti2ewuqyWdeU+Hr/xMgB2H/Vzz0s7WPu9xX3auv7JD1i+dDbzpsf2oL/19BZuWzSDxQnGfFes2kXRhNHcsqDItU5JSQnFxcUA3PH8NpZdOoWln5/cn5egjx3VzaxYtYtVdy0c8LHx+Fo6WPbERrbcf03c/ZFxD9WsB95gx99fS0FuahK0W+xF971O2d9fS5ZHB6Gj417y01Ie+8alXDC5MH1B9ZP5YfLaUoLOAJ3ZeXETdEd2Lr/dfgSAe17aQU/Qcs9LO3jVKas46ue1sloCQUtu9g4glKDDwxyRtlYd556XdsRN0O/ta2T3UT+LZron6PA5y2paXOv46jpYVRc679pyHyV761lXUZfo0uMqq27mYOOJmGsYLF9LB3X+Ttf2IuMeqq5AkP/78o6U9aDdYrcW7nl5h2d70NFxh9+zmZCgk0kJOgOc09F22vJFsyaw6eAxFs2cQOneBs4rDN110dzexYY9Db3jxxPH5lPpa40ZT351+xEWzZzIvDh3cYT3JRqDfnt3HTMmjE5YZ3eggQuc/f1p001hQQ4HG08M6th4fP6O3tcunsi4hyp83akag3aLPRyHRzvQMXG/uv0Ii2dNzIgx6GQOnCpBZ4DmUWMYf7I1ptx/1tjeIY6vz5nGP6/by9cvm0ZFrZ9JhQV8/bJp5GRnMTo/h69fNg0I9aBL9zb0bof9ZssnfP2yqXF70K9uP8KfzpmacIijrLqZogmjY9qNVOLfT7Gzf1153aCHOGZMHMOO6uaE5xoIX0sHz35w2LW9yLiH6t7ffsyfzpmauiEOl9j/5r/K+Pqcqd4d4oiKe2XpQf50zlT1oMV73D7H7Qna3iGOJ0v20xO0/FvJfn658RDjR+fRHQzyWtlRdh/1c+GU/QBUHm2l0tfKv5Xs79PW1qrjPFlygI+mN8Wc5719jQR6LLtq3Ycvnt1UxaTCfNq7e1zrHDzYRQWh864t91HbcpJDx04kvPZ4yqqbKatpibmGwapzhjjc2ouMe6i6AkGeLDlAfm5qetBusVsLT757wLN3cUTHXelr5cmSA8ye7P0edDIpQWeA/gxx+DtC94n6T4Yfu/GfDOB3Pv3uLe/oux0pfEw8/g73faeODySs097d97yJzne680Q+DlX0axctOu6hn6+b/EBqetCJYvef7MZ4NEPHi7s/78EzjRJ0BqgtnMA0f0NMef24ib1DHMuXzGZ1WS33Lp1NoCfYexfHhVMKWVfu496ls4HQEIevpaN3O2zr4SbXuzjKa1tOexfHya5AP+7i8FFcHDrv4cYTQ76LI/oaBsvX0sEHBxpd24uMe6ieef8Qy5fMTuFdHPFj/0XpAZYvme3dIY6ouEv21LN8yeyMGOK4L4ltabGkDPDOjHlEj3JY4P3PJWfyhIh4k3rQGeDqgx/FTFQxwPzdm/h+WS09QcuqHUd6H1eV1TKpMJ/zCvN5rayWt3fX8+ULQ2PVu50x6FU7jvRpb2vVcVbtOEJt88mY87+3r5HzxhZwvL3LNcZVZbUUTRjNOaPdFxWqqA3Q4px3bbkPgK6eYD9egb52OGPQ0dcwWOHb7Nzai4x7qLoCQVbvqE3ZGLRb7NbCqrIjnr3NLjru0Hu2lr11sR+Wn8mUoDPAFH+ja3lPMNS3fqeyPvS4u56G1k4aWjt5Z3d9n3KASp+/z3ak9bvrXcf41lfW0Z0gmTa3d7P9k2Y+dU5su2H1DQGORpx3fWU9eYO43aysphmIfw2D4fN3JGwvOu6heqeyLmX3QSeK/Z3d9Z5N0PHiXl9ZF7cDcSZTgs4A/RmD/vkNc1jwyHp+/s05nDc2v3cM+poLJ7Gu3MfPnbU4wjMJw9u952g+6ToGfby967Rj0OecldvPmYSh83YFgkMeg46+hsEKzyR0ay8y7qFaW+7jZzfMSfFMwtjYX9tZy89vmOPhMei+ce+ta82YmYT/emPy2tIYdAZ4dPFNdGX3/VvalZ3DU9femqaIRCQVlKAzhA3ahNsicuZRgs4A3y99jnzbdwJIvu3hjnVPpykiEUkFLTeaAQ4+8sdkxdxoB0EMM5a/loaIRMRN1SNf1XKjI0l/PiQMLzd64B+u03KjA5Tq5UZ3rvDGcqMHHr7O4x8SFvduZ9Ryo48kr620DHEYY/7MGFNujAkaY+ZG7bvPGLPfGLPHGPOViPIlTtl+Y8y9EeVFxpgtxph9xpiXjDHJ+XZPD3l08U205/T9ssz2nHx9SChyhkvXGPQu4OtAaWShMeZC4AbgImAJ8G/GmGxjTDbwBLAUuBD4plMX4BHgMWvtTOA4cMZlrdUXfYl7l9xFTeFEghhqCidy75K7eGtO/B6fiJwZ0pKgrbW7rbV74uxaBrxore201h4C9gOXOz/7rbUHrbVdwIvAMhNa6eUq4BXn+GeBrw3/FYiIDD+vjUFPBTZHbNc4ZQDVUeXzgXOBZmttIE79GMaY24HbAfLO/2ySQh5+f1K+gZ+sfZyzAqFvVZnmb+Anax/n4TzDobOuoasrwKZNH2CtpaSkhOqaTlrzsygJfkLF0QD19QFKSkoAqG4N0tbW0bsd1tJyku3bt3PicOzYaFPTSXbubCVY6/52qTnSSaA5i5LuKtc6bW1tvedtaOxgV/lxRh2L93c6sYPNPfhbu2KuYbCOdwTp7HRvLzLuoQoGg5SWlpKXnZqxX9fYLZS8W+LZmYTRcbe1tbN161bqxo6sG8+GLUEbY94G4n1l8w+stavcDotTZonf07cJ6sdlrV0JrITQXRxu9bzm+6XP9SbnsLMCndz5znP86GvfoLuhkUv+6HKCJe8y78qFvNkU+lbveVfOwFdWy97OOuZdGZqVNdbnx+zewbwr+37AllO2iVkXfi7uh4T5FVv47OwZzEvwlVdrGsqZMmE0866c7lrnvffe6z3vWQe2UTRzKvMG8a3eedXNZB8sZ96VyflWb5+/g55N78e8JmGRcQ9VYO2bzL1iQcq+UcUtdrv2TeZdudCzCTo6bvNhKRdecimzz/f+h4TJlNbb7IwxJcDfWmu3Otv3AVhr/9HZfhP4oVP9h9bar0TWA34CNADnW2sDxpgrI+slcqbeZndWXjbtXT1xnwMx22Fu5afbN5A6PT09ZGcnjqM/hnLsYNqLjDsZ58rOMilL0G6xt3f1UJCb5dkEHR13sn/nw2n3g0vP2NvsVgO/Mcb8CzAFmAl8SKinPNMYUwQcIfRB4o3WWmuM2QBcT2hc+mbArXeesQZym13Fj5foNrsBGqm32VX8aIlusxsG5sHktZWu2+z+1BhTA1wJvO70lLHWlgMvAxXAWuBOa22PM8Z8F/AmsBt42akLsBz4v8aY/YTGpM+46XW6zU5kZEpLD9pa+yrwqsu+h4GH45T/DvhdnPKDhO7yOGOtvuhLQGgseoq/kdrCCTy6+Ca2zbmG+WmOTUSGj9eGOMTF6ou+1Juow1xvVxGRM8LIumdFRCSDqAedIf6kfEPsEMeC69IdlogMI/WgM0B4oso0fwNZ2N6JKl/e/na6QxORYaQedAZwm6jyV+ue5tGv/znWgsVirfNDnG3nfndrwz9RXwBA/PLeY4i/r8/xp6vjxEM4vtPUT9jOII+N2x7W9dp7z5es+QL2VJup4BZ7f36n6RQdt9v79kyn9aAzgNaDFskcWg96hNF60H1pokr/aT3o1Mv49aBlYDRRRWRkUg86A2iiisjIpASdITRRRWTk0RCHiIhHKUGLiHiUErSIiEcpQYuIeJQ+JMwQWotDZORRDzoDaC0OkZFJPegM4LYWx3fW/pJvzLmGxrYuthw6Rk/QsvngMbYcamJSYT6fn3Y2W5ztzQePAVB51E+lr7V3O2xr1XG2HAy1EW3LwSYumXY2eQm+R2/LoSbqWzuZnWCmV2VTDwXOebccOsb5Zxdwzui8fr8OYWXVzZTVtMRcw2DV+Tuo83e6thcZ91B1BYJsOdSUsu8kdIvdWth86Jhnv5MwOu5KXytbDh6j5WR3GqNKPa3FkQH6sxbHZ849i6pj7VxeNJ4PDzUBxDyHUIL2dwR6t8M+PNREloG5caZ6R7cRT3/qtDQ3c/a4cf2u76asupnOQHBQx8ZT5+/ofe3iiYx7qD481MQffmpcyhK0W+wfHmpi3vRzMB5N0NFxf3ioicKCnIQdAK/4rzu+oLU4RpL+rMXx0u1XsuCR9bz8V1f2ay2Ol//qyj5tpW4tjtB5k7EWR/Q1DFZ4LQ639iLjHqpZD7zBS7dfkeK1OGJjL7rvdV66/UqPr8VxKu6MWovjjuS1pTHoDKC1OERGJvWgM4DW4hAZmZSgM4TW4hAZeTTEISLiUUrQIiIepQQtIuJRStAiIh6lBC0i4lFK0CIiHqUELSLiUUrQIiIepYkqGULrQYuMPOpBZwCtBy0yMmm50Qyw8clvx13NrqZwIgv/+ld9yvJzsugMBOM+B2K2w9zKT7dvIHWCwSBZWYnj6I+hHDuY9iLjHu5zJZtb7J2BIHnZWXh0tdGYuFP9ug3F3oev03KjI8kUf6Nr+XWfP5/SvY2suXshxf9UQtmKa3no9QomjY1cbrSOx2+cA5xabnTt9xb3aev6pz5g+ZKhLDdaTtHE0dzyhemudUpLS1m8OHTeO15wlhu9eLDLjZaz6q4FAz42ntByo++z5f6r4+6PjHuoZv/dWsr+/lryc1OTaNxin/13a9mx4sueXbA/Ou5MWm501MPJa0sJOgMkWg+6ICebgtwsCnKzyc4yFORmO2WRP1m96w9HlkeKPCaa2zF962RRkJOVsE5etjkVR86pOAcq+pqG6nTtRcY9VHk5WeQnMfbTns8ldmNCvwOvrgcdHXd/3oNnIu//f0G0HrTICKUedAbQetAiI5MSdIbQetAiI4+GOEREPEoJWkTEo5SgRUQ8SglaRMSjlKBFRDxKCVpExKOUoEVEPCotCdoY8/+MMZXGmJ3GmFeNMeMi9t1njNlvjNljjPlKRPkSp2y/MebeiPIiY8wWY8w+Y8xLxpi8VF+PiMhwSFcP+i3gYmvtJcBe4D4AY8yFwA3ARcAS4N+MMdnGmGzgCWApcCHwTacuwCPAY9bamcBxQPOfReSMkJYEba1dZ60NOJubgWnO82XAi9baTmvtIWA/cLnzs99ae9Ba2wW8CCwzxhjgKuAV5/hnga+l6jpERIaTF6Z6/yXwkvN8KqGEHVbDqRnN1VHl84FzgeaIZB9ZP4Yx5nbgdoC88z875MDTraOjA19dHV1dATZt+gBrLSUlJVTXdNKan0VJ8BMqjgaorw9QUlICQHVrkLa2jt7tsJaWk2zfvp0Th2NXC2tqOsnOna0Ea93fLjVHOgk0Z1HSXeVap62trfe8DY0d7Co/zqhjewZ83Qebe/C3dsVcw2Ad7wjS2eneXmTcQxUMBiktLSUvOzWryLnGbqHk3RLPLjcaHXdbWztbt26lbuzI+ths2BK0MeZt4Pw4u35grV3l1PkBEAB+HT4sTn1L/J6+TVA/LmvtSmAlhBbsdw0+QxQUFHD+pPHs9Tdw5ZVfwJSup7i4mI1tFUwqLKB48Qz8ZbXUWB/FxZcBofWgf31gB8XFfdcIfnz3B8yZE3896KcPbOGSSxKvB72hZRdFE0ZTvKDItU5JSQnFxcUAvFi9jYsvmkLx5we+HvS46mZW1eyiuHjhgI+Nx9fSQf7vN/bGFi0y7qHKevsNFi9enLJlM11jf/N1ir9Y7NnlRqPjHrOjlLlzM2M96GQatgRtrb0m0X5jzM3AV4Gr7amvdakBPhVRbRpQ6zyPV94IjDPG5Di96Mj6IiIZLS1DHMaYJcBy4IvW2vaIXauB3xhj/gWYAswEPiTUU55pjCkCjhD6IPFGa601xmwAric0Ln0zsCp1V5Je9a0d1Ld20tjWRX1rBz1BS70/VGYMoefOdr2/wzmms8/2qbY6qffHlgOh8jjHRB8/Oj8nYZ3mjmBEHB2nbdP1XP7BHxu3vdPEEhn3UHUFgjS0dqbsq5vcYrc29DvzaAc6Ju7w+/bc0cn5PWSKtHwnoTFmP5APHHOKNltr73D2/YDQuHQA+J619g2n/Drgp0A28Iy19mGnfAah5Dwe2A78b2tt5+liyKTvJByIiWPzaWjtjPsciNkOcys/3b6B1Onq6iIvL6/f9YdyrmS2Fxl3Ms5VWJBDfoqGONxib2jtZMKYPIxHx6Cj407273w4bX3gy0n7TkJ9aWwGmzpuFPOLxlO6r4E1dy9iwSPrOfAP1/HQmtAY9KnvJPTx+I2nxqDjfifhkx+wfOlQvpMwNAZ9Sz/HoO943vlOwkGMQYe+k3AXq+5K3hj0sic2suX++KNyyRyDnvXAG+xccW3ax6CL7nudAw9flzFj0Jn0nYTGmKQl6JH1kaiISAZRghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFhHxKCVoERGPykl3ADJ4R5pP8t7+RhrbunhvXwM9Qct7+xrYuL+R8woLmD15LBud7ff2NQBQebSVSl9r73bY1qrjvLevkY7unpjzvLevkQsmF2KMeyzv7W/kSPNJ/uC8Ma51djX2kO2cd+P+RsaPyWNMwcDfgmXVzZTVtMRcw2D5Wjqo83e6thcZ91B1BYJs3NdIfm5q+kZusVsb+p1lJfidplN03JW+Vjbua6SxrTONUaWeEnSGa2gNvWF/UXow9PjuQSp9oSQcDFo27m/sLQeo9Pn7bEf6xbsH+H3V+LjnWVl6kIpav2scBxtOcLDhBB3dQdc6x493sbk5dN62zgC/2fIJnxxrP90lxiirbnbijb2GwfD5OxK2Fxl3Mvyi9AD5OdlJay+RRLH/4t0DZCX6q5tG8eL+RekBZp9fmKaI0kMJOoNNHTeK+UXjKd3XwAu3zmfBI+t54TvzeWhNBZMKC7ht8QxWl9WyrtzH4zdeBsDuo37ueWkHL3xnfp+2rn/yA5Yvnc286bEJ+ltPb+G2RTNYPGuiaywrVu2iaMJobllQ5FqnpKSE4uLQee94fhvLLp3C0s9PHvB176huZsWqXTHXMFi+lg6WPbHRtb3IuIdq1gNv8Pyt8ynITU2Cdou96L7XeeHW+WR5tAsdHfeSn5by2Dcu5YLJ3k/Qv74teW1pDFpExKOUoEVEPEoJWkTEo5SgRUQ8SglaRMSjlKBFRDxKCVpExKOUoEVEPEoJWkTEo5SgRUQ8SglaRMSjlKBFRDxKCVpExKOUoEVEPEoJWkTEo5SgRUQ8SglaRMSjlKBFRDxKCVpExKPSkqCNMQ8aY3YaY3YYY9YZY6Y45cYY83NjzH5n/2URx9xsjNnn/NwcUf5HxpiPnWN+boxHvwVTRGSA0tWD/n/W2kustZcCa4C/d8qXAjOdn9uBJwGMMeOBFcB84HJghTHmHOeYJ5264eOWpOoiRESGU1oStLXWH7E5GrDO82XAczZkMzDOGDMZ+ArwlrW2yVp7HHgLWOLsK7TWbrLWWuA54GupuxIRkeGTk64TG2MeBm4CWoAvOcVTgeqIajVOWaLymjjlIiIZz4Q6nsPQsDFvA+fH2fUDa+2qiHr3AQXW2hXGmNeBf7TWbnT2vQN8H7gKyLfWPuSU/x3QDpQ69a9xyhcB37fW/rFLTLcTGg4BuBjYNfQrTbkJQGO6gxikTI09U+OGzI09U+MG+Jy1dmwyGhq2HnQ4afbDb4DXCY0x1wCfitg3Dah1youjykuc8mlx6rvFtBJYCWCM2WqtndvPGD0jU+OGzI09U+OGzI09U+OGUOzJaitdd3HMjNj8E6DSeb4auMm5m+MKoMVaexR4E7jWGHOO8+HgtcCbzr5WY8wVzt0bNwGrEBE5A6RrDPonxpjPAUGgCrjDKf8dcB2wn9AQxrcBrLVNxpgHgY+cej+21jY5z/8a+A9gFPCG8yMikvHSkqCttf/LpdwCd7rsewZ4Jk75VkLjyQO1chDHeEGmxg2ZG3umxg2ZG3umxg1JjH3YPiQUEZGh0VRvERGPGnEJ2hizxBizx5kafm+644nHGHPYmb6+I/yJsDFmvDHmLWeq+1vhmZSJpsenIM5njDH1xphdEWUDjtNtGn8aYv+hMeaI87rvMMZcF7HvPif2PcaYr0SUp/T9ZIz5lDFmgzFmtzGm3BjzXafc0697grgz4TUvMMZ8aIwpc2L/kVNeZIzZ4rx+Lxlj8pzyfGd7v7N/+umuyZW1dsT8ANnAAWAGkAeUARemO644cR4GJkSVPQrc6zy/F3jEeX4doQ9GDXAFsCWFcS4GLgN2DTZOYDxw0Hk8x3l+Tppi/yHwt3HqXui8V/KBIuc9lJ2O9xMwGbjMeT4W2OvE5+nXPUHcmfCaG2CM8zwX2OK8li8DNzjlTwF/7Tz/P8BTzvMbgJcSXVOic4+0HvTlwH5r7UFrbRfwIqHp5ZlgGfCs8/xZTk1pd5seP+ystaVAU1TxQOOMO40/TbG7WQa8aK3ttNYeInSX0eWk4f1krT1qrf2987wV2E1o9qynX/cEcbvx0mturbVtzmau82MJTaB7xSmPfs3Dv4tXgKud24DdrsnVSEvQblPGvcYC64wx20xo9iPAJBu67xvn8Tyn3GvXNNA4vRb/Xc5QwDPm1IJcnozd+a/zHEI9uox53aPihgx4zY0x2caYHUA9oT9mB4Bma20gThy9MTr7W4BzBxP7SEvQ8ZYi9eJtLAustZcRWt3vTmPM4gR1M+Wa3OL0UvxPAn8AXAocBf7ZKfdc7MaYMcB/A9+zfRcfi6kapyxtsceJOyNec2ttjw2tvjmNUK/3ggRxJC32kZag3aaSe4q1ttZ5rAdeJfSGqAsPXTiP9U51r13TQOP0TPzW2jrnH2IQ+HdO/ffTU7EbY3IJJblfW2t/6xR7/nWPF3emvOZh1tpmQstMXEFouCg8lyQyjt4Ynf1nExpOG3DsIy1BfwTMdD59zSM0gL86zTH1YYwZbYwZG35OaFr7LkJxhj9pv5lTU9rdpseny0DjjDuNP9VBQ29iC/tTTi2mtRq4wfl0vojQuuMfkob3kzOW+TSw21r7LxG7PP26u8WdIa/5RGPMOOf5KOAaQmPoG4DrnWrRr3n4d3E9sN6GPiV0uyZ3w/nppxd/CH2qvZfQGNIP0h1PnPhmEPqktwwoD8dIaAzrHWCf8zjenvqE+Qnnej4G5qYw1v8k9N/SbkK9g1sHEyfwl4Q+MNkPfDuNsT/vxLbT+cc0OaL+D5zY9wBL0/V+AhYS+m/xTmCH83Od11/3BHFnwmt+CbDdiXEX8PdO+QxCCXY/8F+EVtwEKHC29zv7Z5zumtx+NJNQRMSjRtoQh4hIxlCCFhHxKCVoERGPUoIWEfEoJWgREY9SghYR8SglaBERj1KCFnFhjJnnLOJT4MzwLDfGDObr1UQGRRNVRBIwxjxEaGbYKKDGWvuPaQ5JRhAlaJEEnPUePgI6gC9Ya3vSHJKMIBriEElsPDCG0LeAFKQ5Fhlh1IMWScAYs5rQt3YUEVrI5640hyQjSM7pq4iMTMaYm4CAtfY3xphs4ANjzFXW2vXpjk1GBvWgRUQ8SmPQIiIepQQtIuJRStAiIh6lBC0i4lFK0CIiHqUELSLiUUrQIiIepQQtIuJR/z+B6FR+B6aFogAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot the location of the source\n", "fig, ax = plt.subplots(1,1, figsize=(5, 6))\n", "mesh.plotGrid(ax=ax)\n", "ax.plot(mesh.gridEy[solenoid_inds,0], mesh.gridEy[solenoid_inds, 2], 'ro')\n", "ax.set_xlim([-1, 3e3])\n", "ax.set_ylim([-3000, 3000])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## set up the forward simulation\n", "\n", "Including the survey and problem" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "rxList = [] # see FDEM.Rx for receiver types\n", "srcList = [FDEM.Src.RawVec_e(rxList, f, src_vec) for f in freqs]" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# define a problem - the statement of which discrete pde system we want to solve\n", "wires = Maps.Wires(('sigma', mesh.nC), ('mu', mesh.nC)) # mapping to work with a vector m = [sigma, mu] as the model\n", "prob = FDEM.Problem3D_e(mesh, sigmaMap=wires.sigma, muMap=wires.mu) \n", "prob.solver = Solver\n", "\n", "survey = FDEM.Survey(srcList)\n", "\n", "# tell the problem and survey about each other - so the RHS can be constructed \n", "# for the problem and the resulting fields and fluxes can be sampled by the receiver. \n", "prob.pair(survey) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the simulation" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 52.3 s, sys: 9.88 s, total: 1min 2s\n", "Wall time: 51 s\n" ] } ], "source": [ "%%time\n", "fields_background = prob.fields(np.r_[sigma_background, mu_0*mur_background])\n", "fields = prob.fields(np.r_[sigma, mu_0*mur])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the magnetic fields" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "def streamplot_b(field_to_plot, ax):\n", " # Plot the magnetic flux \n", " max_field = np.abs(field_to_plot).max() #use to set colorbar limits\n", " cb_range = 5e2 # dynamic range of colorbar\n", " \n", " cb = plt.colorbar(mesh.plotImage(\n", " field_to_plot, \n", " vType='F', view='vec', \n", " range_x=xlim*1.1, range_y=ylim*1.1,\n", " pcolorOpts={\n", " 'norm': LogNorm(), \n", " 'cmap': plt.get_cmap('viridis')\n", " },\n", " streamOpts={'color': 'w'}, mirror=True, ax=ax, \n", " clim=[max_field/cb_range, max_field]\n", " )[0], ax=ax)\n", "\n", " ax.set_xlim(xlim)\n", " ax.set_ylim(ylim)\n", " return ax\n", " \n", "def plot_b(\n", " field,\n", " freq_ind=0, # which frequency would you like to look at?\n", " reim='real', # real or imaginary part\n", " xlim=np.r_[-1000, 1000],\n", " ylim=np.r_[-1000, 1000],\n", " ax=None\n", "):\n", " if ax is None: \n", " fig, ax = plt.subplots(1,1, figsize=(6,5))\n", " \n", " field_to_plot = getattr(field[srcList[freq_ind], 'b'], reim)\n", " ax = streamplot_b(field_to_plot, ax)\n", " cb.set_label('|B {reim}|'.format(reim=reim))\n", "\n", " # give it a title\n", " ax.set_title(\n", " 'B {reim}, {freq:10.2f} Hz'.format(\n", " reim=reim, \n", " freq=freqs[freq_ind]\n", " )\n", " )\n", " plt.show()\n", " return ax\n", "\n", "def plot_difference(\n", " field, \n", " field_background,\n", " freq_ind=0, # which frequency would you like to look at?\n", " reim='real', # real or imaginary part\n", " ax=None, \n", " xlim=np.r_[-1000, 1000],\n", " ylim=np.r_[-1000, 1000]\n", "):\n", " \n", " if ax is None: \n", " fig, ax = plt.subplots(1,1, figsize=(6,5))\n", " \n", " field_to_plot = getattr(field[srcList[freq_ind], 'b'], reim)\n", " background = getattr(field_background[srcList[freq_ind], 'b'], reim)\n", " ax = streamplot_b(field_to_plot - background, ax)\n", " \n", " cb.set_label('|B {reim}|'.format(reim=reim))\n", "\n", " # give it a title\n", " ax.set_title(\n", " 'B {reim}, {freq:10.2f} Hz'.format(\n", " reim=reim, \n", " freq=freqs[freq_ind]\n", " )\n", " )\n", " plt.show()\n", " return ax" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2e38b97d0a9a4f0b8bd5c1cbfe4b59a5", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(IntSlider(value=0, description='freq_ind', max=7), ToggleButtons(description='reim', opt…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xlim=np.r_[-200., 200.]\n", "ylim=np.r_[-200., 200.]\n", "\n", "def plot_fields(freq_ind, reim, model): \n", " if model == \"difference\":\n", " return plot_difference(fields, fields_background, freq_ind, reim=reim, xlim=xlim, ylim=ylim)\n", " else:\n", " return plot_b(fields if model == \"with target\" else fields_background, freq_ind, reim, xlim, ylim)\n", "\n", "ipywidgets.interact(\n", " plot_fields,\n", " freq_ind=ipywidgets.IntSlider(min=0, max=len(freqs)-1, value=0), \n", " reim=ipywidgets.ToggleButtons(options=['real', 'imag']), \n", " model=ipywidgets.ToggleButtons(\n", " options=[\n", " 'background', # wholespace\n", " 'with target', # wholespace with target\n", " 'difference' # response due to target (target - background)\n", " ]\n", " ),\n", ")" ] }, { "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.6" } }, "nbformat": 4, "nbformat_minor": 2 }