"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plot a cross section of the conductivity model\n",
"fig, ax = plt.subplots(1,1)\n",
"cb = plt.colorbar(mesh.plotImage(np.log10(sigma), ax=ax, mirror=True)[0])\n",
"\n",
"# plot formatting and titles\n",
"cb.set_label('$\\log_{10}\\sigma$', fontsize=13)\n",
"ax.axis('equal')\n",
"ax.set_xlim([-120., 120.])\n",
"ax.set_ylim([-100., 30.])\n",
"ax.set_title('Conductivity Model')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Set up the Survey\n",
"\n",
"Here, we define sources and receivers. For this example, the receivers sample $db/dt$, the time-derivative of the magnetic flux. The source is a vertical magnetic dipole with unit moment. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# Define the receivers, we will sample the real secondary magnetic flux density as well as the \n",
"# imaginary magnetic flux density \n",
"\n",
"dbdt_z = TDEM.Rx.Point_dbdt(locs=rx_loc, times=times, orientation='z') # vertical db_dt\n",
"rxList = [dbdt_z] # list of receivers"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"There is 1 source. Each source has 1 receivers sampling the resulting db/dt-field at 30 times\n"
]
}
],
"source": [
"# Define the list of sources - one source for each frequency. The source is a point dipole oriented\n",
"# in the z-direction\n",
"\n",
"srcList = [\n",
" TDEM.Src.MagDipole(\n",
" rxList, loc=src_loc, orientation='z', waveform=TDEM.Src.StepOffWaveform()\n",
" )\n",
"]\n",
"\n",
"print(\n",
" 'There is {nsrc} source. Each source has {nrx} receivers '\n",
" 'sampling the resulting db/dt-field at {ntimes} times'.format(\n",
" nsrc = len(srcList), \n",
" nrx = len(rxList),\n",
" ntimes = len(times)\n",
" )\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Set up Forward Simulation \n",
"\n",
"A forward simulation consists of a paired SimPEG problem and Survey. For this example, we use the E-formulation of Maxwell's equations, solving the second-order system for the electric field, which is defined on the cell edges of the mesh. This is the `prob` variable below. The `survey` takes the source list which is used to construct the RHS for the problem. The source list also contains the receiver information, so the `survey` knows how to sample fields and fluxes that are produced by solving the `prob`.\n",
"\n",
"\n",
"\n",
"**An Aside: Mappings** \n",
"\n",
"The `sigmaMap` defines a [mapping](http://dev-docs.simpeg.xyz/content/api_core/api_Maps.html) which translates an input vector to a physical property - specifically, $\\sigma$. This comes in handy when you want to invert for $\\log(\\sigma)$, in that case, the model is defined in terms of log-conductivity, and you would need an `ExpMap` to translate $\\log(\\sigma) \\to \\sigma$"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The maximum time is 1.3e-04. \n",
" There are 80 timesteps, 6 of them distinct (which is the same as the number of matrices that need to be factored)\n"
]
}
],
"source": [
"# solve the problem at these times\n",
"timeSteps = [(1e-8, 20), (1e-7, 20), (2e-7, 10), (1e-6, 10), (2e-6, 10), (1e-5, 10)] \n",
"\n",
"# define a problem - the statement of which discrete pde system we want to solve\n",
"prob = TDEM.Problem3D_e(mesh, timeSteps = timeSteps, sigmaMap=Maps.IdentityMap(mesh)) \n",
"prob.solver = Solver\n",
"\n",
"survey = TDEM.Survey(srcList)\n",
"\n",
"# tell the problem and survey about each other - so the RHS can be constructed for \n",
"# the problem and the resulting fields and fluxes can be sampled by the receiver. \n",
"prob.pair(survey) \n",
"\n",
"print(\n",
" 'The maximum time is {:1.1e}. \\n There are {} timesteps, '\n",
" '{} of them distinct (which is the same as the number of matrices that need to be factored)'.format(\n",
" prob.times[-1], prob.nT, (len(timeSteps))\n",
" )\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solve the forward simulation \n",
"\n",
"Here, we solve the problem for the fields everywhere on the mesh, for with and without the sphere. "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"solving with sphere ... \n",
"... done \n",
"CPU times: user 6.31 s, sys: 415 ms, total: 6.72 s\n",
"Wall time: 4.93 s\n"
]
}
],
"source": [
"%%time\n",
"print('solving with sphere ... ')\n",
"fields = prob.fields(sigma)\n",
"print('... done ')"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"solving without sphere ... \n",
"... done \n",
"CPU times: user 6.51 s, sys: 487 ms, total: 7 s\n",
"Wall time: 6.62 s\n"
]
}
],
"source": [
"%%time\n",
"print('solving without sphere ... ')\n",
"fields_background = prob.fields(sigma_background)\n",
"print('... done ')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot the fields\n",
"\n",
"Lets look at the physics!"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# log-scale the colorbar\n",
"from matplotlib.colors import LogNorm \n",
"import ipywidgets"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "bb711db917ac4490b952a8724589df31",
"version_major": 2,
"version_minor": 0
},
"text/html": [
"Failed to display Jupyter Widget of type interactive
.
\n",
"\n",
" If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
" that the widgets JavaScript is still loading. If this message persists, it\n",
" likely means that the widgets JavaScript library is either not installed or\n",
" not enabled. See the Jupyter\n",
" Widgets Documentation for setup instructions.\n",
"
\n",
"\n",
" If you're reading this message in another frontend (for example, a static\n",
" rendering on GitHub or NBViewer),\n",
" it may mean that your frontend doesn't currently support widgets.\n",
"
\n"
],
"text/plain": [
"interactive(children=(IntSlider(value=1, description='time_ind', max=79, min=1), Output()), _dom_classes=('widget-interact',))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def plot_dbdtSphere(\n",
" time_ind=0 # which frequency would you like to look at?\n",
"# ax=ax\n",
"):\n",
" fig, ax = plt.subplots(1,1, figsize=(6,5))\n",
" \n",
" # Plot the magnetic flux\n",
" field_to_plot = fields[srcList, 'dbdt', time_ind]\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=[-120., 120.], range_y=[-180., 80.],\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([-120., 120.])\n",
" ax.set_ylim([-180., 70.])\n",
" cb.set_label('|db/dt|')\n",
"\n",
" # plot the outline of the sphere\n",
" x = np.linspace(-sphere_radius, sphere_radius, 100)\n",
" ax.plot(x, np.sqrt(sphere_radius**2 - x**2) + sphere_z, color='k')\n",
" ax.plot(x, -np.sqrt(sphere_radius**2 - x**2) + sphere_z, color='k')\n",
"\n",
" # plot the source and receiver locations\n",
" ax.plot(src_loc[0],src_loc[2],'co', markersize=6)\n",
" # ax.plot(rx_loc[0,0],rx_loc[0,2],'co', markersize=6)\n",
"\n",
" # plot the surface of the earth\n",
" ax.plot(np.r_[-200, 200], np.r_[0., 0.], 'w:')\n",
"\n",
" # give it a title\n",
" ax.set_title(\n",
" 'db/dt, {time:10.2e} s'.format( \n",
" time=prob.times[time_ind]\n",
" )\n",
" )\n",
" plt.show()\n",
" return ax\n",
"\n",
"ipywidgets.interact(\n",
" plot_dbdtSphere, \n",
" time_ind=ipywidgets.IntSlider(min=1, max=len(prob.timeSteps)-1, value=1) \n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot data with and without sphere\n",
"\n",
"In what follows, we plot db/dt data with and without the sphere"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"dpred = survey.dpred(sigma, f=fields)\n",
"dpred_background = survey.dpred(sigma_background, f=fields_background)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5,0,'time (ms)')"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEaCAYAAAAyinE1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd8zWf7wPHPlUGIRBA7RsggQ4m9Q0ujtFpVpWiVPjx2q4sOPPoU1aE6VFXVVkqVoqW/1mztvWKv2DFiBknu3x8n8qSpkXGS70lyvV+v83LOd145t3Ouc3/v+3vfYoxBKaWUshcnqwNQSimVs2hiUUopZVeaWJRSStmVJhallFJ2pYlFKaWUXWliUUopZVeaWJRKQUQmich/E583FJG9VsfkSETEiIif1XEox6WJRTkMEWkgIn+JSIyIXBCRP0WkppUxGWNWGWMCrYwhM4hINxGJFJErInJGRBaJiIfVcamcwcXqAJQCEBFPYCHQE5gN5AEaAjetjCu7ExEXY0xcimWNgeFAhDFmi4gUBh53hNhUzqA1FuUoAgCMMTONMfHGmBvGmKXGmO0AIuIkIu+IyFEROSsiU0SkYOK6cBGJSn4wETkiIo8kPh8qIrMT97kiIrtEpEaybauJyObEdbMAt2Tr/nbsxOO+JiLbE2tWs0Qk+fZviMgpETkpIi/d77KRiJQSkQWJtbMDIvKvZMtvJH7hJ48xWkRcE193FZE9InJRRJaISLlk2xoR6S0i+4H9dzl1TWCNMWZL4nt+wRgz2RhzJXH/SSIyTkR+S3xPViQ/fqJHRGR/4vm/FBFJdv40xSYilRLPdUFE9opIu7u9Xyr70MSiHMU+IF5EJotICxEplGJ9l8RHE6ACUAD4Ig3HfwL4HvACFtzZV0TyAD8BU4HCwA/A0w84VjsgAvAFqiTGhYhEAAOARwA/oPEDjjMTiAJKAW2B4SLysDHmJLAmRRzPAXOMMbdF5EngLaANUBRYlXis5J4EagNBdznvOuBREfmPiNQXkbx32aYj8B7gDWwFpqdY3wpbgnoI2/vxaOJ7kKbYRMQd+A2YARQDOgBjRST4LjGp7MIYow99OMQDqAxMwvZlG4ctARRPXPc70CvZtoHAbWyXc8OBqBTHOgI8kvh8KPB/ydYFATcSnzcCTgKSbP1fwH8Tn//t2InH7ZTs9ShgXOLzicCIZOv8AAP43eVvLQPEAx7Jlo0AJiU+fwn4I/G5AMeBRomvfwG6JdvPCbgOlEt8bYCmD3ivWwA/A5eAq8AngHPiuknA98m2LZAYa5lkx2+QbP1sYGB6YgOeBValiO1rYIjV/x/1kf6H1liUwzDG7DHGdDHG+AAh2H7Jf5q4uhRwNNnmR7ElleKpPPzpZM+vA24i4pJ43BMm8Rst2bHTcqwCyWI8nmxd8ucplQIumMTLT8nOWzrx+RygroiUwpb8DLZf/wDlgDEicklELgEXsCWf0smOdb9zY4z5xRjzOLZaWmtsta6X7ra/MeZq4jlKJVt/r/cgrbGVA2rf2T5xn45AifvFrxybNt4rh2SMiRSRSUCPxEUnsX0J3VEWW63mDLYvvPx3VoiIM7bLMKlxCigtIpIsuZQFDqYj7FOAT7LXZe6z7UmgsIh4JEsuZYETAMaYSyKyFNtlpsrAzGTxHQfeN8akvDyVXKqGLTfGJAC/i8gf2JL5P2IXkQLYEtDJVBwyrbEdB1YYY5qlJl6VPWiNRTmExAbcV0XEJ/F1GWzX29cmbjITeEVEfBO/6IYDs4ytV9E+bDWQlomN2+8Ad2s3uJs12BJUPxFxEZE2QK10/hmzgRdFpLKI5AcG32tDY8xxbJfcRoiIm4hUAbrx97aMGcDz2NpaZiRbPg4YdKcdQkQKisgzqQ1SRFqLSHsRKSQ2tbC1B61NttljYuv+nQdbW8u6xJgfJK2xLQQCRKSziLgmPmqKSOXU/j3K8WhiUY7iCrYG3XUicg3bl9xO4NXE9ROxNbCvBA4DsUBfAGNMDNALmIDtF/81bO00D2SMuYWtobkLcBHbNf8f0/MHGGN+AT4DlgEHsCUtuHeX6Q5AeWw1gXnY2hV+S7Z+AeAPnDHGbEt2nnnAB8D3InIZ2/vUIg2hXgT+ha1X1mVgGvBhilrGDGAItktZ1bFdnnqgtMaWWFtrDrTH9j6cTtw/tT8MlAOSv19aVkrZS+Kv7p1AXpON7tdIvAQZZYx5x+pYVPakNRal7EhEnhKRPIndpT8Afs5OSUUpe9DEopR99QDOYWv8j8c2koBSuYpeClNKKWVXWmNRSillV5pYlFJK2VWuvEHS29vblC9fPs37JSQkAODkpPnYkWi5OB4tE8eU0XLZtGlTtDHmgTcf58rEUr58eTZu3Jjm/a5csd0g7eGh01Y4Ei0Xx6Nl4pgyWi4i8qDhjgC9FKaUUsrONLEopZSyK00sSiml7CpXtrEopbKH27dvExUVRWxsrNWh5Aipbbx3c3PDx8cHV1fXdJ1HE4tSymFFRUXh4eFB+fLlSTb7sUqn+Ph4AJydne+5jTGG8+fPExUVha+vb7rOo5fClFIOKzY2liJFimhSyUIiQpEiRTJUS9TEkgZHL9xg9cELnLtyr1HQlVL2pkkl62X0PdfEkga/7DpLz+93UvP9/6P28P+j66QNfLx0L7/uPM3xC9fRcdeUyn0ee+wxLl26xKVLlxg7dmzS8uXLl9OqVSu7niszjpkZckQbi4g0xDYRkQsQZIyplxnn6VizNNXLFuTQxdvsPnmZnSdjWL73LAmJ+aRgPleCS3kmPgoSXMqTCkUL4Oykv7iUyqkWL14MwJEjRxg7diy9evWyOKJ7i4uLw8Ul87/2LU8sIjIRaAWcNcaEJFseAYwBnIEJxpiR9zqGMWYVsEpEngQ2ZFasHm4u1CznRdOQ/921Gns7nsjTV9h5IoZdJy+z+2QMk9cc5VacrfeFm6sTlUr8PdkElvDAzfXejWdKKccwatQo3Nzc6NevH6+88grbtm3jjz/+4Pfff+e7775j2rRpSSN5DBw4kIMHD1K1alWaNWtGy5YtuXr1Km3btmXnzp1Ur16dadOm/eMy02effca4ceNwcXEhKCiI77//nqFDh3Lw4EFOnDjB8ePHeeONN/jXv/4FcM9jbtq0iQEDBnD16lW8vb2ZNGkSJUuWJDw8nHr16vHnn3/SqlUrOnfuTO/evTl27BgAn376KfXr17fr+2Z5YgEmAV8AU+4sEBFn4EugGbYpZjeIyAJsSWZEiv27GmPOJj5/DnjpQSdMSEhIGtogLe61T0UvZyp6FaZ1cGEAbscncPj8DSJPX2XPmatEnr7Kgq0nmL7OVpDOAr7e+alUvACVShSgcvECBBZ3p2C+9HXty+3SU5Yqc9mrTBISEpJ6Mr23aA97Tl22y3HvqFzSk3dbVr7n+vr16zN69Gh69+7Nhg0buHnzJrGxsaxcuZL69esnxRYfH8/777/Pzp072bRpE2C7bLVlyxa2b99OqVKlaNiwIStXrqRBgwZ/O8fIkSM5cOAAefPm5dKlS8THx5OQkMD27dv5888/uXbtGjVq1CAiIoL4+Pi7HrN27dr06dOHefPmUbRoUWbPns1bb73FhAkTMMZw8eJF/vjjD+Lj43n++efp168fDRo04NixYzz22GPs3LnzH397er8nwQESizFmpYiUT7G4FnDAGHMIQES+B1obY0Zgq938g4iUBWKMMXf9nyci3YHuAGXKlLFP8Pfg6uxEQDF3Aoq58wTFAVsXvqhLsUSeucqe01eJPH2NdUcusXDn2aT9SnjmJbCYOwHF3QksXoDAYu6ULZwPJ228VMoS1atXZ/PmzVy5coW8efMSFhbGxo0bWb16NZ9++ukD969ZsyY+Pj4AVK1alaNHj/4jsYSGhtK5c2dat25N69atk5Y//vjj5MuXj3z58hEeHs6GDRsoWLDgXY/p5eXFrl27iIiIAGyJrkSJEknHeuaZZ5Ke//HHH0RGRia9vnz5MleuXLHruG6WJ5Z7KA0cT/Y6Cqj9gH26Ad/da6UxZjwwHqBGjRomI29ievcN8vQkqGyxvy07d+Umu09dZk/iY/fJy6w+dJH4xIab/HmcCSzhQeWSnlQu6Yl/sQL4ertTzCOv9pZJQQc8dDwZLRMnJ6ekey6GPhHygK3tz9nZmfLlyzNlyhTq169PlSpVWLlyJYcOHSIkJCTpM+js7JwUZ/J/3dzckl67uLiQkJDwj3tIFi9ezMqVK1mwYAHvv/8+u3btwsnJCRFJ2vbO83sd08nJieDgYNasWfOPv0FE8PT0TNonISGBNWvWkC9fvvv+7U5OTukuP0dNLHf7xrxvlytjzJBMiiVTFfXIS2OPojQO+N9I1LG349l/5qot0SQ+ft52khmJl9LAlnDKF3HHt6g7vkXcKe/tjm/io1B+V006StlJo0aN+Oijj5g4cSKhoaEMGDCA6tWr/+Mz5uHhkeZLRwkJCRw/fpwmTZrQoEEDZsyYwdWrVwGYP38+gwYN4tq1ayxfvpyRI0eyb9++ux4nMDCQc+fOsWbNGurWrcvt27fZt28fwcHB/9i2WbNmfPHFF7z++usAbN26lapVq6Yp7gdx1MQSBSS/XuUDnLQolizn5upMqE9BQn0KJi0zxnAyJpZD565yOPpa0mPXiRh+3Xk6qYYDtt5pFYq6E1TSk6BSngSV9KRSCU/y5dEOA0qlVcOGDXn//fepW7cu7u7uuLm50bBhw39sV6RIEerXr09ISAgtWrSgZcuWDzx2fHw8nTp1IiYmBmMMr7zyCl5eXgDUqlWLli1bcuzYMd59911KlSp1z8SSJ08e5syZQ79+/YiJiSEuLo6XX375ronl008/pV+/flSpUoW4uDgaNWrEuHHj0viu3J9DzHmf2May8E6vMBFxAfYBDwMnsPX0es4Ys8se56tRo4bJSfOx3I5P4PiF6xw5f43D0dc5HH01qcZzOTYOACcBX293gkoV/FvCKeqR1+LoM85RyyU3s1eZ7Nmzh8qV7924nlMNHTqUAgUK8Nprr9n1uKkZ0uWOu733IrLJGFPjQftaXmMRkZlAOOAtIlHAEGPMtyLSB1iCrSfYRHsllZzI1dmJCkULUKFogb8tN8Zw4tINdp9MvKR28jJbjl3k523/q/wV9chLSClPqvh48VCZgoSW9soRyUYpZR3LE4sxpsM9li8GFmdxODmKiOBTKD8+hfLTPPh/PURibtxO6iiw6+Rldp6IYcW+/Uk3epYq6EaoT0FbsvHxIrR0QQrm167QSmWVoUOHWh1ChlieWFTWK5jPlToVilCnQpGkZdduxrHr5GW2R11ie1QM26MusWTXmaT15Yvk56EyXlQr40W1soUIKuWJq7OOCKSU+idNLAoA97wu1PItTC3fwknLYq7fZseJGLZFXWJ71CXWHjrP/K22y2h5XZyo4lOQamULEVbWi7CyhSjm6WZV+EopB6KJRd1TwfyuNPD3poG/d9Kyk5dusPnYRbYcu8TmYxeZ9OcRxq+0DV9T2isf1RKTTM3yhQkq5anjpCmVC2liUWlSyisfpbzy0apKKQBuxsWz6+RlNh9NTDZHL7Jw+ykAPPK6UKN8IWpXKEIt38KEli6ol8+UygU0sagMyeviTFjZQoSVLZS07FTMDdYfvsC6wxdYd+g8y/aeA2w3dVYvV4javoWpXaEIVXwKktdF761RjuvIkSO0atXqrmNppdby5cv56KOPWLhwoR0js48uXbrQqlUr2rZta9fjamJRdleyYD5aVy1N66qlAduwNesPX2D94fOsO3yBj5babvLK6+JEWNlCtsttft6ElC6ol86USsYYgzHmgXPUO5rsFa3Klop65KVllZL8p3UIv77ciC3vNuPrztXpWLscl27c5sMle2n95Z9U/+9v9Jq+iRnrjnH8wnWrw1YKsM1h8sILL1ClShXatm3L9evXGTZsGDVr1iQkJITu3bsnTfJ34MABHnnkER566CHCwsI4ePDg3461YcMGqlWrxqFDhzh37hzNmjUjLCyMHj16UK5cOaKjozly5AiVK1emV69ehIWFcfz4cWbOnEloaCghISG8+eabSccrUOB/967NmTOHLl26ALaaSL9+/ahXrx4VKlRgzpw5gC1R9evXj6CgIFq2bMnZs2fJDFpjUVmukHseHg0uwaOJ99acu3KTvw5Gs2p/NKv3R7N4x2kAyhXJT30/bxr6eVOvorfeS6MIDw//x7J27drRq1cvrl+/zmOPPfaP9V26dKFLly5ER0f/45LP8uXLH3jOvXv38u2331K/fn26du3K2LFj6dOnD4MHDwagc+fOLFy4kMcff5yOHTsycOBAnnrqKWJjY5PGAgP466+/6Nu3L/Pnz6ds2bL06dOHpk2bMmjQIH799VfGjx//t3N+9913jB07lpMnT/Lmm2+yadMmChUqRPPmzfnpp5948skn7xv3qVOnWL16NZGRkTzxxBO0bduWn376ib1797Jjxw7OnDlDUFAQXbt2feB7kFaaWJTlinrkTbp0Zozh4LlrrN5/jtUHolmw1Tb4ppNAtbKFeLhyMR6uVJyA4gV0oE2VJcqUKZM0EVanTp347LPP8PX1ZdSoUVy/fp0LFy4QHBxMeHg4J06c4KmnngLAze1/3e/37NlD9+7dWbp0KaVK2Tq+rF69mnnz5gEQERFBoUL/a6csV64cderUAWy1nPDwcIoWtQ1U27FjR1auXPnAxPLkk0/i5OREUFAQZ87Y7klbtWoV7du3x9nZmVKlStG0aVN7vEX/oIlFORQRwa9YAfyKFaBLfV9uxyew7fglVu6PZlnkWUb9updRv+7Fp1A+HqlcnKaVihFcNA95XPSqbm5wvxpG/vz577ve29s7VTWUlFL+gBERevXqxcaNGylTpgxDhw4lNjaW+427WLJkSWJjY9myZUtSYrnf9u7u7knP77dd8thiY2P/ti5v3v8NzZT8GFnxg0w/jcqhuTo7UaN8YQY0C+Dnvg1YO+hhRrQJpVIJD77fcIznJ66n4eg1vDJnNz9sPE701ZtWh6xymGPHjiXNczJz5sykibq8vb25evVqUvuFp6cnPj4+/PTTTwDcvHmT69dtbYVeXl4sWrSIt956Kym5NWjQgNmzZwOwdOlSLl68eNfz165dmxUrVhAdHU18fDwzZ86kcePGABQvXpw9e/aQkJCQVPu5n4YNGzJr1izi4+M5deoUy5YtS+e7cn9aY1HZSomCbnSoVZYOtcoSezuevw5G88u2E6w4cJ7/2xuNCFQr48VjoSVpEVqS0l73n8xIqQepXLkykydPpkePHvj7+9OzZ08uXrxIaGgo5cuXp2bNmknbTp06lR49ejB48GBcXV354YcfktYVL16cn3/+mRYtWjBx4kSGDBlChw4dmDVrFo0bN6ZkyZJ4eHgkzcdyR8mSJRkxYgRNmjTBGMNjjz2WNNPkyJEjadWqFWXKlCEkJOQf+6b05JNPsmzZMkJDQwkICEhKUPbmEMPmZ7WcNmx+bnflyhWMMRy7Yvgj8ixLdp1m10nbDNVVy3jRMrQkLUJL4FMov8WR5h46bP6D3bx5E2dnZ1xcXFizZg09e/Zk69atmXrOXDNsvlL2ICKElPYkpHRB+j3sz5HoayzeeYrFO07x/uI9vL94Dw+V8aJlaAlahJSkTGFNMspax44do127diQkJJAnTx6++eYbq0OyG00sKkcq7+1Or3A/eoX7cfT8NRbtsCWZ4YsjGb44kod8CvJYaEkef6gUpfRymbKAv78/W7ZssTqMTKGJReV45Yr8L8kcO389KcmM+CWSkb9G0sDPm7bVfXg0uARurjrEjFIZpYlF5Spli+SnZ3hFeoZX5Oj5a8zdfIK5m6Lo//1WPNxcePyhUjxT3YeqZbz0PhkHYYzRsshiGW1718Sicq1yRdwZ0CyAlx/2Z+2h8/ywKYofN0cxY90x/IoVoG11H9pUK63zzFjIzc2N8+fPU6RIEU0uWcQYw/nz5/92g2daaa+wNLhy5Qpz584lNDSU6tWrZ0JkKj3s2VvvcuxtFm8/xQ+both09CJOAo0DivJszTI8Urk4Ljrsf6rYq0xu375NVFTUP27+U+mTkGCbO+lBg1q6ubnh4+ODq+vfh1FKba8wTSxpcPHiRWrWrMnhw4fp168fw4YN067HDiCzuoEfOneVOZui+HHzCU5fjqW0Vz5eqFeOZ2uWpWA+HbfsfrRrvmPKaLmkNrHoz680cHFxYdmyZfTo0YMxY8YQFBTE/PnzrQ5LZZIKRQvwRkQl/hzYlPGdq+NTKB/DF0dSd8TvDJ6/k0Pn7n8zmlK5ldZY0iB5tl+zZg09evRgx44d7Nixg5CQEHuHqVIpK38d7zwRw3d/HuHnbSe5FZ/Aw5WK0bWBL/UqahtAclpjcUxZVWPRxJIGKQvl9u3b/Pbbb0lDda9YsYIGDRqk6q5WZT9WfImdvRLL9LXHmL7uKNFXbxFY3IOuDcrTumpp7bKMJhZHpZfCsgFXV9ekpLJ3716aNm1K7dq12bx5s8WRqcxWzMONV5oFsPrNpnzYtgpOTsKbc3dQb+QffPb7fi7H3rY6RKUso4nFTgICApg5cyYnTpygZs2aDBgw4IEDwqnsz83VmWdqlGFxvwbM/FcdqpXx4pPf9tFg5B98+n/7iLmhCUblPnopLA1SU428dOkSgwYNYty4cfj7+7Nr165/dNlT9uVol112nohhzO/7+W33GTzcXHixvi/d6vvmqhkwHa1MlI22sWSirBjdeM2aNezZs4euXbtijOHs2bMUL148zedUD+aoX2K7Tsbw2e/7WbLrDB55XehSvzzdGvjilT+P1aFlOkctk9xO21iyubp16ybNJT1//nwqVKjAqFGjuH1bL43kFsGlCvJ15xr80r8hDQO8+fyPAzT4YBkfLonk4rVbVoenVKbRxJIFwsLCaN68OW+++SZhYWGsXr3a6pBUFqpc0pOxHauz5OVGNA4sytjlB2nwwR98uCSSazfjrA5PKbvTxJIFypYty7x585g/fz4xMTE0bNiQN9980+qwVBYLLOHBl8+FseTlRjStXJwvlx3k4Y9XsHD7yQwP+qeUI9HEkoWeeOIJdu/ezeuvv85DDz0E2GZ00y+V3CWguAefd6jG3J71KOyehz4zttDp23UcOKu9CFXOoIklixUoUIBRo0bx3HPPAfD5558THh6e6VOSKsdTvVwhfu7bgGGtg9kRFUOLMSsZ+YteHlPZnyYWixUpUoTdu3cTFhbGSy+9xOnTp60OSWUhZyfh+brl+eO1cJ6sWppxKw7yyCcrWLT9lNZkVbalicVinTt3Zv/+/QwYMIApU6bg7+/PlClTrA5LZTHvAnn58JmHmNuzHoXy56H3jM10/na9Xh5T2ZImFgfg5eXFRx99xO7du3nkkUcoV64cALGxsfqrNZdJfnlsW9SlpMtj12/p5TGVfWT7xCIiQSIyW0S+EpG2VseTEX5+fsybN4/GjRsD8Oabb1K/fn3WrVtncWQqK925PLbstXBaJ14eaz56Jav3R1sdmlKpYmliEZGJInJWRHamWB4hIntF5ICIDHzAYVoAnxtjegLPZ1qwFggLC+Pw4cPUqVOHTp06cfz4catDUlnIu0BePnrmIWb3qIursxOdvl3HG3O26fhjyuFZOqSLiDQCrgJTjDEhicucgX1AMyAK2AB0AJyBESkO0TXx3yHAdaCeMab+g84bFhZmVqxYkeZ4rRim4sqVK4wePZrPP/8cJycnJkyYwOOPP55l588OcsPwIbG34xm3+hiT1hynsHse3o7w4+FAb6vDuqfcUCbZUUbLxdPT0/GHdDHGrAQupFhcCzhgjDlkjLkFfA+0NsbsMMa0SvE4m/joDQwE7nmtQES6i8hGEdkYHZ19Lil4eHgwePBgNm3axNNPP42vry9gm/ulb9++/P777zpMTC7g5urMy018mf5iNQrnd+XlObt59cfdRF/VoWGU43GxOoC7KA0kv+YTBdS+18YiUh54C3AHPrzXdsaY8cB4sA1CmZFfUlb8CgsODmbq1KlJr0+ePMncuXOZPHkyhQsX5qmnnqJdu3Y88sgjODll+6azdMkNv47reHiwsH8Jxq88xJj/28/6o5sY8ngQT1Yt7ZAzWOaGMsmOMrtcHPEb6G6fjnterzPGHDHGdDfGdDTG5JpBuLp37865c+f46aefiIiIYNasWbz44otJ6w8ePKg9ynIoV2cnejfxY3H/BlTwdueVWdt4cdIGTly6YXVoSgGOmViigDLJXvsAJy2KxaG5ubnRunVrpk+fzrlz51iyZAlOTk4kJCTQsGFD6tSpw6JFizTB5FB+xTz44d/1GPp4EOsPX6D5JyuYtvaolreynCMmlg2Av4j4ikgeoD2wwOKYHJ6bmxshISEAGGMYNmwYZ8+epVWrVtSqVYuff/5Zv3ByIGcnoUt9X5a83IiwcoV456eddJ20gXNXblodmsrFrO5uPBNYAwSKSJSIdDPGxAF9gCXAHmC2MWaXlXFmN87Ozrz00kvs27ePCRMmcP78eZ544glWrlxpdWgqk5QpnJ8pXWvxnyeC+fPgeVqMWcmyyLNWh6VyKZ1BMg2yaxfK27dvs3DhQp588klEhLFjx1KyZElat26dIxr6s2u5ZJZ9Z67Qb+YWIk9f4YW65Rj0WGXcXJ2zNAYtE8ekM0gqu3F1deWpp55CREhISGDChAm0adOGatWqMWfOHBISEqwOUdlRQHEPfupdn24NfJm85iiPf76a3ScvWx2WykU0seQyTk5OrF+/nqlTpxIbG8szzzxD9erVWb9+vdWhKTtyc3Xm3VZBTOlai0s3bvPkl38yYdUhEhJy3xUKlfU0seRCLi4udOrUid27dzN16lSuXLmCm5sbAHFxOthhTtIooGjSlMj/XbSHF75bz5nLsVaHpXI4TSy5mLOzM506dWLv3r1UqVIFgBdeeIG2bduyd+9ei6NT9lLYPQ/jO1dn+FOhbDhygYhPV7Jkl877ozKPJhaFs7OtYdcYQ6VKlViyZAnBwcF0796dEydOWBydsgcR4bnaZVnYtyGlC+Wjx9RNjPo1Urugq0yhiUUlERHeffddDh48SK9evZg0aRKFFVW2AAAgAElEQVR+fn4sWKC3EeUUfsUK8GPP+nSoVYaxyw8yeP4ubXdRdqeJRf1DsWLF+Oyzz4iMjKR9+/bUrm0bqu3YsWPExur1+ewuj4sTw58KpUfjCkxde5TXfthGXLz2DFT2o4lF3VOFChX47rvvKF68OMYY2rdvT6VKlZgxY4Z2Uc7mRISBEZV4/dFAftxygl7TN3MzLt7qsFQOoYlFpYqI8N///pfChQvTsWNHatWqxfLly60OS2WAiNC7iR9DHw9i6e4zvDR5o06BrOxCE4tKtaZNm7Jx40amTJnC2bNnadKkCT/88IPVYakM6lLfl4+eeYg/D0TT+dv1OkOlyjBNLCpNnJyc6Ny5M3v37uXTTz9Nms1y3bp1nDlzxuLoVHq1re7Dl8+FsT3qEh3GryX6qg5iqdJPE4tKl3z58tG/f3/c3NxISEigU6dO+Pn58d5773Ht2jWrw1Pp0CK0JBNeqMmh6Ku0+3oNJ3V+F5VOmlhUhjk5ObFo0SKaN2/O4MGDCQgIYMqUKdrAnw01DijK1G61OXf5Js+MW8ORaP2RoNJOE4uyi4CAAObOncuqVavw8fHhhRde4Ndff7U6LJUONcsXZmb3Oty4Hc8zX68h8rQOYKnSRhOLsqsGDRqwZs0a5s+fT4sWLQD47bffOH/+vMWRqbQIKV2Q2T3q4CTQYfxaDmvNRaWBJhZld05OTjzxxBOICNeuXePZZ5/F39+fsWPHEh+v90pkF37FPJjVvS4iwovfrefCtVtWh6SyCU0sKlO5u7uzatUqqlatSu/evalevTqrVq2yOiyVSuW93fnm+RqcionlpckbiL2tPwzUg2liUZkuODiY33//ndmzZ3PhwgUaN27M/v37rQ5LpVL1coX49NmqbDl+iVdnb9OxxdQDaWJRWUJEeOaZZ9izZw+zZ8/G398fsLW/3Lqll1gcXYvQkrzVojKLdpzig18jrQ5HOThNLCpLubu707ZtWwAOHjxIREQEVatW1eFhsoGXGvryfN1yfL3yEFPXHrU6HOXANLEoy1SsWJH58+dz48YNmjRpQseOHTl16pTVYal7EBEGtwri4UrFGDJ/J39E6kgL6u40sShLtWrVit27d/Puu+8yZ84cqlWrxo0bese3o3JxduKzDtUIKuVJnxlb2HkixuqQlAPSxKIsly9fPoYNG8aOHTsYPXo0+fLlwxjD7t27rQ5N3YV7XhcmvlCTQvnz0HXSBk7o0C8qBU0symEEBATQoUMHABYtWkRwcDAvvfQS0dHRFkemUirm6cZ3L9bkxq14un63gcuxOiKy+h9NLMohhYeH8/rrrzN58mQCAwP55ptvdOwxBxNQ3INxnatz8NxVek7bxK04LR9lo4lFOaQCBQowatQotm7dSkhICN27d6dNmzZWh6VSqO/nzcinq/DngfO8NW8Hxug9LgpcrA5AqfsJDg5m+fLlTJs2jfz58wMQFxdHbGwsBQoUsDg6Bba5XI5fuM6Y3/cTWNyDfzWqYHVIymJaY1EOT0To3LkzTz/9NABffPEFlStXZt68efoL2UG8/Ig/jwYXZ9SSSO0pplKfWESkXGYGolRq1a5dm8KFC9OmTRsef/xxjhw5YnVIuZ6IMLJNFQq75+HlWVu5oWOK5WppqbHMS7lAROrYMRalUqVu3bps2rSJjz/+mOXLl1O7dm2mTZtmdVi5XiH3PHz8TFUOnL3Kx78fsjocZaEHJhYRaSciIwEPEaksIs7JVo/PvNCUujcXFxcGDBjAnj17aNasGb6+vgDac8xiDfy9+VdDX2ZtOsWK/ToHT26VmhrLn8BuoBDwCbBfRDaLyEJA74xSlipTpgzTpk2jfv36AAwYMIAuXbpw9uxZiyPLvV57NJDAYu4MXriPs1dirQ5HWSA1iaUK8H/AE8aYFsaYCsAjwBCgaWYGp1RaGGPw8PBg+vTpBAYG8tVXX+nEYhbI6+LMyCcrce1WPK//sF07WORCqUksTwMLgFki8puIfAg8BtwE9OeIchgiwnvvvcf27dupVq0avXr1ok6dOuzcudPq0HIdv6LuvPpwBVbsO8fkv45YHY7KYg9MLMaYl4wxNYCPgX3AYaAJsB7QsbOVw6lcuTK///47M2bM4Pz587i5uVkdUq7UvnpJmgQWZfgvkew7c8XqcFQWSkuvsBeNMb2NMWONMd2AhsDqTIrrnkSkgoh8KyJz7rdM5W4iQocOHdi/fz9+fn4AdOvWje+++04b+LOIiDCq7UN4urnQb+YWndY4F0lLYrksItXvvDDGbAIC0nIyEZkoImdFZGeK5REisldEDojIwPsdwxhzKDGx3XeZUgDOzrZOjFevXmXv3r107dqVhg0bsm3bNosjyx2KeuTlw7YPEXn6Ch8u2Wt1OCqLpGVIl27AVBHZDWwCQoG0Dmk6CfgCmHJnQWL35S+BZkAUsEFEFgDOwIgU+3c1xmS4u09CQgJXrqS9ap6efVTmS225LFq0iJkzZ/LOO+8QFhZGjx49ePvtt/H09MzkCHOf5GVSo3Q+2lcvxberD1OrjDv1KhS2MLLcLau+w1JzH0tdERFjzD6gHrAYKA7swdaIn2rGmJXAhRSLawEHEmsdt4DvgdbGmB3GmFYpHulOKiLSXUQ2ishGHYY9d3JycqJjx45s2rSJF198kblz52qvsSzy6sO+VPDOz9sL9nHxug6xn9PJg7oCisg4bF/++4BfgV+NMafTfUKR8sBCY0xI4uu2QIQx5qXE152B2saYPvfYvwjwPrYazgRjzIi7LbtfDDVq1DAbN25Mc+x3sr2Hh0ea91WZJ73lcvnyZTw9PYmLi6NXr1707duX0NDQzAgx17lbmew6GcNTX/5F48CijO9cHRGxKrxcK6PfYSKyKbEz132lplfYv40xYcBQbDdJThKRNSIyXEQapbgTP12x3u2094nnfGJMFe8kkLstU+pB7lwC27dvHz/++CPVqlXjlVdeISZGB1HMDMGlCvJGRCC/7T7D9xuOWx2OykSpbrw3xkQaY0YbYyKw3Ri5GngGWJfBGKKAMsle+wAnM3hMpVItKCiIvXv38tJLLzFmzBgCAwOZNm2a3tiXCbrW96WBnzfvLdzNqRgduCOnSu+w+V8Dy4wxfbFdfsqIDYC/iPiKSB6gPbYbMpXKMkWKFGHcuHGsX7+esmXLMmLECOLi4qwOK8dxchKGPxVKXIJh5C+RVoejMkl6E4sT8JWI5ANeSe1OIjITWAMEikiUiHQzxsQBfYAl2DoEzDbG7EpnXEplSI0aNVi7di1Lly7F1dWVy5cv88477+jlMTsqWyQ/PRpVYP7Wk6w/nLIvj8oJ0ptYDmNrc/kKcE/tTsaYDsaYksYYV2OMjzHm28Tli40xAYltJO+nMyal7MLJyYnSpUsDsHTpUoYPH05gYCBTp07Vy2N20ivcj1IF3RiyYBfxCfqe5jTpTSzfGGOOYEsuEXaLRikH07ZtW9atW0e5cuV4/vnnadiwIVu3brU6rGwvXx5n3m4ZxJ5Tl5mx/pjV4Sg7S1diMcYcS/z3iDEm2L4hKeVYatasyZo1a/j222/Zt28fb731ltUh5QiPhZagboUifLx0Lxev3bI6HGVHaZmauJGIrBSRXSIyQ0RqZWZgSjkSJycnunbtyr59+xg/3ja/3ZEjRxg/frzeZJlOIsKQJ4K4EhvHx7/pcC85SVpqLBOB94BwbEOyfCoi7TIjKKUclZeXFz4+PgBMnDiRHj16ULt2bdauXWtxZNlTpRKedK5TjhnrjrHrpHaQyCnSkliijTG/GWPOGWN+BZoDgzMpLqUc3n/+8x+mT5/OqVOnqFu3Ll27duXMmTNWh5XtvPJIAF758/CfBbu1c0QOkZqxwqaIyMvAahEZLCJ3Bq7Uib5UriYiPPfcc0RGRvLGG28wdepUhg8fbnVY2U7B/K68/mgg649cYME2vTc6J0hNjeVbbEOsFAaeBA6IyP8BkdgGpFQqV/Pw8OCDDz5gx44dDB5sq8SvW7eOH374QX+Bp1K7GmUILV2Q4Yv3cO2m3pia3aUmseQHfjDGdE0cM6witpsihwD5MjM4pbKTSpUqUaRIEQDGjh1Lu3btqF+/PmvWrLE4Msfn7CQMfSKYM5dv8uWyA1aHozIo1XPei8hxEfkNGAk8BGwF7jspl1K51cSJE5kwYQKHDx+mXr16PPvssxw6dMjqsBxa9XKFaBNWmgmrDnMk+prV4agM0DnvlcoEzs7OdOvWjf379zNkyBAWLlzIvHnzrA7L4Q2MqEQeFyfeW7jb6lBUBmS7Oe+Vyk4KFCjA0KFD2b9/P3362KYYmjNnDmPGjOHWLb0pMKVinm70e9iP3yPPsiwyw5PFKotk6Zz3SuVWpUqVIm/evIBtiuSXX36Z4OBgfvrpJ23gT6FLPV8qFHVn2MLd3IzTm0+zo7Qklq7AWBH5TkT6iMjXpH3Oe6VyvYkTJ7J48WLy5MnDU089RdOmTdm2bZvVYTmMPC5ODHk8mMPR15i4+ojV4ah0SMtEX/v535z3VUnHnPdKKdv9Ly1atGDbtm2MHTuWnTt3cuCA9oRKrnFAUZoFFefzP/ZzOkZvl8tu0jQIpTEm3hjzAxBmjPnUGHM+k+JSKsdzcXGhZ8+eHDx4kDZt2gDwySef8N5773H9+nWLo7Peuy2DiEswfLRUxxHLbtI7bP7d5qlXSqWDp6cnIraP1Pbt2xk8eHDS9MgJCQkWR2edskXy82K98szdHMXuk5etDkelQbrnY7FrFEopACZNmsTKlSspXrw4nTt3pk6dOmzevNnqsCzTq4kfBfO5MuKXPVaHotIgvfOxjLV3IEopm4YNG7J+/XomT57MqVOnkmotuXF4/oL5XOnb1J9V+6NZue+c1eGoVEpvjUUplYmcnJx4/vnnOXToEDVq1ACge/futGvXjt27c9fNg53rlKNs4fwMX7xHpzHOJjSxKOXAXF1dATDGULZsWX755RdCQkLo3LlzrulJlsfFiTciAok8fYUfN0dZHY5KBU0sSmUDIsKQIUM4fPgwr732GnPnzqVSpUpMmzbN6tCyRMvQklQt48XHS/dx41buuySY3WhiUSob8fb2ZtSoURw6dIi+ffsSHh4OwN69ezlx4oS1wWUiEeHtlpU5fTmWb1frYJ6OThOLUtlQiRIlGD16dNI0yf3796dixYq8/PLLnDyZMyfLqlm+MM2DijNuxSGir960Ohx1H5pYlMoBxo0bR6dOnfjiiy/w9fWlV69eHD2a8wYff7NFJW7cjmfM/+23OhR1H5pYlMoBypcvz4QJE9i/fz9dunRhwoQJzJgxw+qw7K5i0QI8V6ssM9Yf4+C5q1aHo+5BE4tSOYivry9ff/01Bw8eTBqmf9asWXTu3JnIyEiLo7OP/o/4k8/VmQ9+yRl/T06kiUWpHKhMmTJ4eHgAcPr0aX788UeCgoJ49tln2bFjh8XRZYx3gbz8u3EFlu4+w/rDF6wOR92FJhalcrj+/ftz5MgRBg4cyC+//EKVKlV49dVXrQ4rQ7o1qEAJTzeGL96j89k4IE0sSuUCRYsWZfjw4Rw5coQhQ4ZQp04dAG7cuEFUVPa76TBfHmcGNA9g6/FLLNpxyupwVAqaWJTKRQoXLszQoUN55plnABgzZgz+/v4MGjSImJgYi6NLm6fDfKhUwoNRv+7VmSYdjCYWpXKx9u3b07ZtW0aOHEnFihUZM2YMN29mj3tEnJ2EQY9V5tiF60xbe8zqcFQymliUysXKly/P1KlT2bx5M9WqVePll1/mxRdftDqsVGscUJSG/t58/sd+Ym7oTOmOQhOLUopq1arx22+/sWTJEl577TUATp06xbJlyyyO7MEGtahMzI3bjF2WOwblzA40sSilkjRv3pywsDAARo8eTdOmTWnZsqVDD9UfVMqTNtV8+O6vI5y5HGt1OApNLEqpexg2bBgffPABf/75J1WqVKFnz56cOXPG6rDuqv/D/sQnGMatOGh1KIpsmFhEpIKIfCsic5Itqywi40Rkjoj0tDI+pXIKNzc33njjDQ4cOECvXr2YMGECb731ltVh3VXZIvlpU600M9Yd46zWWiyXpYlFRCaKyFkR2ZlieYSI7BWRAyIy8H7HMMYcMsZ0S7FsjzHm30A7oIb9I1cq9/L29uazzz5j586dDBs2DIAdO3Ywffr0pGmTHUGfpn7EJRjGrdBh9a3mksXnmwR8AUy5s0BEnIEvgWZAFLBBRBYAzsCIFPt3NcacvduBReQJYGDi8e8rISGBK1eupDn49OyjMp+WS9YoVaoUYHu/v/jiC8aPH88nn3zC8OHDqVev3t+2taJMCueBViHFmL7uKJ1rFMe7QJ4sj8HRZVW5ZGmNxRizEkg5uE8t4EBiTeQW8D3Q2hizwxjTKsXjrkkl8dgLjDH1gI53Wy8i3UVko4hsjI6OttefpFSuNGrUKL7++mvOnDlDREQEHTt2dIipkrvXL0tcfAIT1xy3OpRcLatrLHdTGkj+vyAKqH2vjUWkCPA+UE1EBhljRohIONAGyAssvtt+xpjxwHiAGjVqmDsD9KVHRvZVmUfLJWt1796dTp06MXr0aEaMGEGVKlX473//+7dtsrpMgj08eLKaDz9sOUnfZpUo5uGWpefPLjK7XBwhschdlt1zVDljzHng3ymWLQeW2zUqpdQD5c+fn7fffptu3brh7u4OwNKlS9m2bRvdunV7wN6Zo09TP+ZtiWL8ikO80yrIkhhyO0foFRYFlEn22gfImXOrKpVDlShRIulX8Ny5c3njjTeoW7cuixYtyvLRh3293XmyammmrTvKuSvZY3ianMYREssGwF9EfEUkD9AeWGBxTEqpdBo3bhyzZs0iISGBVq1aERERwa5du7I0hj5N/bgVl8A3q7SHmBWyurvxTGANECgiUSLSzRgTB/QBlgB7gNnGmKz9X6iUshsRoUWLFqxdu5bRo0ezfv161q5dm6UxVChagNZVSzNlzRGir2qtJatJbpwkp0aNGmbjxo1p3u9OVz1tJHYsWi6OJ3mZnD9/Hi8vL5ydnfn222+5ePEiffv2JW/evJkaw8FzV2n2yQr+1bACgx6rnKnnyi4y+lkRkU3GmAfeK+gIl8KUUjlYkSJFcHZ2BmDFihW8/vrrBAcH8+OPP2Zq+0vFogV44qFSTFlzlPNaa8lSmliUUllmypQpLFmyBDc3N55++mnCw8PZvn17pp2vT1N/YuPiGa9tLVlKE4tSKks1b96crVu38tVXX7Fnzx4y84Zlv2K2WsvUNUe5cO1Wpp1H/Z0mFqVUlnNxceHf//43R44coWnTpgC8++67DBs2jOvXr9v1XH2b+nHjdrz2EMtCmliUUpbJnz8/AMYYDh06xJAhQwgMDGTatGl2G+DSr5gHraqUYvJfR7TWkkU0sSilLCciTJ8+nRUrVlC8eHE6d+5M3bp17db+0i+x1jJBay1ZQhOLUsphNGrUiPXr1zNp0iTOnDmDm5ttrK+M9h7zL+5By9CSTP7rCBe11pLpNLEopRyKk5MTL7zwAgcPHiQgIACA9u3b079//ww19Pd72J/rt+OZsFprLZlNE4tSyiHdufclLi4OLy8vvvjiCypWrMgHH3zAjRs30ny8gOIePBZaksl/HdVaSybTxKKUcmguLi58/fXX7Nixg0aNGjFw4EACAwNZv359mo/Vr6k/127FMXzxnkyIVN2hiUUplS0EBQXx888/s2zZMgIDA/Hz8wMgJiYm1ccILOFB3yZ+/LApirmbojIr1FxPE4tSKlsJDw/nt99+o3DhwiQkJNCkSRMiIiLYsWNHqvbv/0gAdSoU5p2fdnLgrE5rnRk0sSilsq34+Hg6duzIunXrqFq1Kj169ODs2XvOYA6As5PwWftquOd1ptf0zdy4FZ9F0eYemliUUtmWq6srr776KgcPHqRv375MnDgRf39/Nm3adN/9inm6MfrZquw/e5XB83dmUbS5hyYWpVS2V7hwYT799FN27NhBhw4dCA0NBeDEiRP3vAemoX/RpPaWOdreYleaWJRSOUalSpUYN24cefLk4erVq9SsWZOHH374nnfw32lvefennew/o+0t9qKJRSmVI7m5ufH222+zbds2qlWrdtf2l5TtLddvxVkUbc6iiUUplSO5uLjQu3dvDhw4QL9+/ZLaXw4d+vud93faWw6cu8qQ+Toruj1oYlFK5WiFChVi9OjR7Ny5k/79++Pr6wvAgQMHktpftL3FvjSxKKVyhcDAQIYNG4aIEBUVRZUqVYiIiGD37t1A8vtbdrBP21syRBOLUirXKV68OB988AHr16+nSpUq9OvXj5hLF/msfTUK5HWht7a3ZIgmFqVUruPq6krfvn3Zv38/PXr04Msvv6RSpUrkSYhNam8ZrO0t6aaJRSmVa3l7e/Pll1+ydetW3n77bby8vGjoX5R2FYU5m6KYvfG41SFmSy5WB6CUUlYLDQ1Nuqly48aNjPpXS0pXbcwrlzrz14EwXns0EJ9C+S2OMvvQGotSSiUTEhLC8OHDubR/I6cn/JvvRrxBg0FTGflLJJdjb1sdXragiUUppZJxc3Nj0KBB7Nu3j359+3Br/5+cnvoaY/9vF+EfLmfyX0e4HZ9gdZgOTROLUkrdRalSpRg9ejTHjh5lwby5LB7wCAHFCtCnT29q9R7DrztP3XMcstxOE4tSSt1H0aJFad68OSGlCzLy0ZK4Ht/A1q8H0PrRJjTs/RFbjl6wOkSHo4lFKaVSqXz58pw4dpTPv/iCguYaf371BrVrhPHCpws4dv661eE5DE0sSimVBvny5aNP796cOHqIbyZOwrecD2vPQpOPl/OvzxZw+Oxlq0O0nCYWpZRKB1dXV1568QX2blzNqoHN6VC9JJPe6UZAYCAte/+Ho+dyb4LRxKKUUhlUzNON99pUZcL4ryjqXYTFY4fi5x9Aq15DOJYLE4wmFqWUsgMR4YX2bTmxbweTZ/1IsWLFWPTVMOr2G8N/ft7F2SuxVoeYZTSxKKWUHYkIz7d7iqi92/hh4VKea/M4U9YcJaRNX1r9+22iomOsDjHTaWJRSqlMICK0bdmMD5+pyv+90gj36D0s+no4vhUq0n7A+8Rcy7k1mGyXWESkgoh8KyJzki0LF5FVIjJORMItDE8ppf7Bt2gBDm1exdS5CylcrCSzRr9DCd9KDB7/Y468iz9LE4uITBSRsyKyM8XyCBHZKyIHRGTg/Y5hjDlkjOmWcjFwFXADdPo3pZTDERE6tWnJ6f3b+XDcZFydYcLaEzw6eiULtxzPUXfxS1b+MSLSCFsCmGKMCUlc5gzsA5phSwobgA6AMzAixSG6GmPOJu43xxjTNvG5kzEmQUSKA58YYzreL46wsDCzYsWKNMd/5YptVjkPD48076syj5aL49EyebD4+HhWHrzEp8sOs37Se3g4x/HBf//DU+E1Mu2cGS0XT0/PTcaYBwaYpTUWY8xKIOX4B7WAA4k1kVvA90BrY8wOY0yrFI+z9zjunbrkRSDv3bYRke4islFENkZHR9vpL1JKqfRxdnamSUAR5rwURkS9apzdu4kXWj9CjZadWLPzoNXhZYgjzMdSGkg+m04UUPteG4tIEeB9oJqIDDLGjBCRNsCjgBfwxd32M8aMB8YD1KhRw2Tkl5T+CnNMWi6OR8skdeaM/5ijQ16jS79BrPhpOhGNlvDUq6P46u3uFPNws/v5MrtcHCGxyF2W3fP6nDHmPPDvFMt+BH60c1xKKZVlypUuybK5k9i8cyDd+g9i843ChH+4nGeD3OnfqjpeBbLPRGOO0CssCiiT7LUPcNKiWJRSylJhIZXY8vs8lr37FI0DvBnxWg9KlA/gtVHfEJdNepA5QmLZAPiLiK+I5AHaAwssjkkppSxV3tudsR2r88F/h5DH1ZWP3+xOUf+qjPvhV6tDe6Cs7m48E1gDBIpIlIh0M8bEAX2AJcAeYLYxZldWxqWUUo5IROjf5VnOH9tHn8Efcu1cFD3btSC8zygiTzvuGGRZ2t3YUdSoUcNs3LgxzftpF0rHpOXieLRMMsf5SzF0Hzicfd4NuBbvRMNCVxnUriGVypdK1f4ZLRcRcbzuxkoppdKviFdB5o77gFVvPUqXOmWZOfIVgisH8nj3N7l45ZrV4SXRxKKUUtmMV/48DG4dysL5P+ETUIWF34yiZDk/Xv1gHPEO0MCviUUppbKp5g1rcXTbn3w++QfyuOXnk4E9qdP7Y9YcPG9pXJpYlFIqm+vzfFvOH43k9VHjcC0XRodv1vLogNGs3hppSTyaWJRSKgdwdXVl1Os9+OO1cF5uUp7fxw+jUc2HaPxsD46fydoajCYWpZTKQdxcnXn50WDWrVtL5bqPsHL2eCpU9OeFN0YQe+t2lsSgiUUppXKg6sGB7Fq5iJkLf8erhA9TPnyLiP/MYt2Ri5l+bk0sSimVg7Vv2ZQz+7by2YyFeJYJJD4h8+9ddIRBKJVSSmUiJycn+nZoSaeYy7g43W3cXzufL9PPoJRSyiFkRVIBTSxKKaXsTBOLUkopu9LEopRSyq40sSillLIrTSxKKaXsShOLUkopu9LEopRSyq5y5QySIhID7E+xuCAQc5fNUy73BqIzKbQHuVeMmX2c1G7/oO3utz617/+9lllVLlaVSVr2SW+5ZHS5flbSv52jflbKGWOKPnArY0yuewDjU7PsbsuBjY4Ud1YcJ7XbP2i7+61P7ft/n2WWlItVZZIV5ZLR5fpZsX+ZpLVcrPqs5NZLYT+nctn9llvBXrGk9Tip3f5B291vfVrefy2TtO2T3nKx13Ir6GcldefJFLnyUlhGiMhGY0wNq+NQf6fl4ni0TBxTVpRLbq2xZMR4qwNQd6Xl4ni0TBxTppeL1liUUkrZldZYlFJK2ZUmFqWUUnaliUUppZRdaWKxMxFxF5FNItLK6lgUiEhlERknInNEpKfV8SgbEXlSRBLBpAIAAASrSURBVL4Rkfki0tzqeBSISAUR+VZE5mT0WJpYEonIRBE5KyI7UyyPEJG9InJARAam4lBvArMzJ8rcxR5lYozZY4z5N9AO0K6vdmCncvnJGPMvoAvwbCaGmyvYqUwOGWO62SUe7RVmIyKNgKvAFGNMSOIyZ2Af0AyIAjYAHQBnYESKQ3QFqmAbLsENiDbGLMya6HMme5SJMeasiDwBDAS+MMbMyKr4cyp7lUvifh8D040xm7Mo/BzJzmUyxxjTNiPxuGRk55zEGPP/7d1NiFVlHMfx78/epBe0dFFO0ZSmEZKOi4RaaBG00wqKIAjJFpK2qE1J2K4SLCqDCt3MriiEmEl7IXAyegGjlHGYRaJEg9AbvcxChJh/i/PMcLjN3Obc+1zvofv7rObcee5znjt/Dj+ee+7c/xFJ/Q0P3wacjIhTAJLeATZHxIvAv97qknQncBlwC3BW0qGImOrowv/HctQkzTMEDEk6CDhY2pTpWhGwG/jQodK+XNdKLg6W5vqAH0vHE8D6uQZHxLMAkrZQ7FgcKvlVqomkjcD9wCXAoY6urLdVqgvwBHA3sEjSioh4q5OL61FVr5UlwPPAgKSdKYBa4mBpTrM89p/vHUbEYP6lWFKpJhExAox0ajE2o2pd9gJ7O7cco3pNfgO25Tixb943NwFcVzq+FjjTpbVYwTWpJ9elfrpWEwdLc0eBmyTdIOli4CFgqMtr6nWuST25LvXTtZo4WBJJbwNfAaskTUjaGhF/AzuAj4Fx4N2IGOvmOnuJa1JPrkv91K0m/rixmZll5R2LmZll5WAxM7OsHCxmZpaVg8XMzLJysJiZWVYOFjMzy8rBYmZmWTlYzCqQtFjS46XjZTkaI81xrnslPZdhnpck3ZVjTWbz4X+QNKsgfTX5B9M9Lzp8ri+BTRHxa5vzXA/sjwh3arTzwjsWs2p2A8slHZO0R1L/dNc+SVskvS9pWNJpSTskPSXpO0lfS7oqjVsu6aPUwvpzSTc3nkTSSuDcdKhIGpT0pqTDkk5J2pC6Bo5LGkxjLkjjTkgalfQkQET8ACyRdPX5+RNZr/PX5ptV8wywOiLWwswOpmw1MEDRRfQk8HREDEh6BXgEeBXYB2yLiO8lrQfeABrfqroDaGyAdWUatwkYTmMeA45KWkvRGbCv1EFwcem536bxB1p72Wbz52Axy+twREwCk5L+pAgAgFHgVkmXA7cD7xVNFIGiCVmja4BfGh4bjoiQNAr8FBGjAJLGgH7gM+BGSa8DB4FPSs/9GVjW7oszmw8Hi1le50o/T5WOpyiutwXAH9M7nibOAovmmLs878zcEfG7pDXAPcB24EHg0TRmYZrTrON8j8WsmkngilafHBF/AaclPQBF7/cUBo3GgRVV5pa0FFgQEQeAXcC60q9XAidaW7VZNQ4WswpS+9Yv0g3yPS1O8zCwVdJxYAzYPMuYIxS9x2drLzuXPmBE0jFgENgJIOkiipD6psX1mlXijxub1ZSk1yjuq3za5jz3AesiYleelZk15x2LWX29AFyaYZ4LgZczzGM2L96xmJlZVt6xmJlZVg4WMzPLysFiZmZZOVjMzCwrB4uZmWX1D38NV91h+pOAAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plot\n",
"fig, ax = plt.subplots(1,1)\n",
"\n",
"ax.loglog(1e3*times, -dpred)\n",
"ax.loglog(1e3*times, -dpred_background, '--k')\n",
"ax.grid(True, color='k',linestyle=\"-\", linewidth=0.1)\n",
"ax.legend(['with sphere', 'background'])\n",
"\n",
"ax.set_title('Sounding over Sphere')\n",
"ax.set_ylabel('-$db_z/dt$')\n",
"ax.set_xlabel('time (ms)')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"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.4"
}
},
"nbformat": 4,
"nbformat_minor": 1
}