diff --git a/docs/source/api_Acous.rst b/docs/source/api_Acous.rst index a3b3791..f11063d 100644 --- a/docs/source/api_Acous.rst +++ b/docs/source/api_Acous.rst @@ -14,17 +14,22 @@ \renewcommand{\dudt}{\frac{\partial\u}{\partial t}} \renewcommand{\dphidt}{\frac{\partial\phi}{\partial t}} \renewcommand{\dsdt}{\frac{\partial s}{\partial t}} + \newcommand{\M}{{\mathbf M}} + \newcommand{\MfMui}{{\M^f_{\mu^{-1}}}} Acoustic wave equation ********************** +Backgrounds +----------- + Time domain acoustic wave equation in 1st order form: .. math:: - \rho\frac{\partial \phi}{\partial t} = \nabla \cdot \vec{u}+\frac{\partial s}{\partial t}\delta(\vec{r}-\vec{r}_s) + \mu^{-1}\frac{\partial \vec{u}}{\partial t} = \nabla \phi - \mu^{-1}\frac{\partial \vec{u}}{\partial t} = \nabla \phi + \rho\frac{\partial \phi}{\partial t} = \nabla \cdot \vec{u}+\frac{\partial s}{\partial t}\delta(\vec{r}-\vec{r}_s) where \\(\\rho\\) is the density, \\(\\mu\\) is the adiabatic compression modulus, \\(\\delta(\\vec{r}-\\vec{r}_s)\\) represents point source, and \\(\\phi\\) is the pressure. @@ -39,13 +44,13 @@ By applying Fourier transform to frequency domain we have: .. math:: - -\rho\omega^2\phi - \div\mu\grad\phi = -\omega ^2 s \delta(\vec{r}-\vec{r}_s) + -\rho\omega^2\phi - \div\mu\grad\phi = -\omega ^2 s \delta(\vec{r}-\vec{r}_s) where \\(\\omega = 2\\pi f\\) is angular frequency. Assuming contant \\(\\mu\\) and \\(\\rho\\) we have Helmholtz equation: .. math:: - (\grad^2 + k^2)\phi = \mu^{-1}\omega^2 s \delta(\vec{r}-\vec{r}_s) + (\grad^2 + k^2)\phi = \mu^{-1}\omega^2 s \delta(\vec{r}-\vec{r}_s) where \\(\ k = \\omega\\sqrt{\\frac{\\rho}{\\mu}}=\\frac{\\omega}{v} \\) is the wave propagation constant. @@ -53,24 +58,311 @@ where \\(\ k = \\omega\\sqrt{\\frac{\\rho}{\\mu}}=\\frac{\\omega}{v} \\) is the Discretization of problem ------------------------- -Artificial Boundary condition ------------------------------ +To compute the solution of time or frequency dependent partial differential equation (PDE) as discussed above, we need to discretize those equations in both time and space. Our domain in real situation (i.e., earth) can be considered as infinite in space. However, in discrete space where we compute the solution of those PDE should be finite. Therefore, we need to implement artificial boundary conditions to get rid of out going wave, which propagates beyond our domain of interest. We first implement sponge boundary condition then perfectly matched layer (PML) boundary condition. + + +.. raw:: html + :file: examples/center_nobc.html Sponge boundary =============== +A fundamental idea of sponge boundary condition is adding damping term to acoustic wave equation: + +.. math:: + + \rho\frac{\partial^2 \phi}{\partial t^2} + c\frac{\partial \phi}{\partial t}- \div\mu\grad\phi = \frac{\partial^2 s}{\partial t^2}\delta(\vec{r}-\vec{r}_s) + +Similarly in frequency domain we have: + +.. math :: + + -\rho\omega^2\phi +\imath\omega c\phi- \div\mu\grad\phi = -\omega ^2 s \delta(\vec{r}-\vec{r}_s) + + +The 1st order form of time domain equation can be written as: + +.. math :: + + \mu^{-1}\frac{\partial \vec{u}}{\partial t} = \nabla \phi + + \rho\frac{\partial \phi}{\partial t} +c\phi- \nabla \cdot \vec{u}= +\frac{\partial s}{\partial t}\delta(\vec{r}-\vec{r}_s) + +We discretize 2nd order acoustic wave equation in time (CITE) with central difference: + +.. math :: + + \rho\triangle t^{-2}(\phi^{n+1}-2\phi^{n}-\phi^{n-1}) +\frac{c}{2}(\phi^{n+1}-\phi^{n-1})- \nabla \cdot \grad \phi^{n} = q + +where + +.. math :: + + q = \triangle t^{-2}(s^{n+1}-2s^{n}-s^{n-1})\delta(\vec{r}-\vec{r}_s) + +By rearranging them, we have: + +.. math :: + + \phi^{n+1} = (1+\frac{c\rho}{2}\triangle t)^{-1}[2\phi^{n}-\phi^{n-1}+\frac{c\rho}{2}\triangle t\phi^{n-1} + +\triangle t^2\rho\div\mu\grad\phi^n + \triangle t^2\rho q] + +By letting \\(\ f = (1+\\frac{\c \\rho}{2}\\triangle t )^{-1}\\), we can substitute \\(\ c = \\frac{1-\f}{\f}\\triangle t^{-1}2 \\rho\\). Range of \\(\ f \\) for time and frequency domain are 0.98-1 and 0.9-1, respectively. And 35 cells are used for the sponge boundary. + +Using staggered grid, we discretize acoustic wave equation in time (CITE): + +.. math :: + + \mu^{-1}\triangle t^{-1}(\u^{n+1}-\u^{n}) = \nabla \phi + + \rho\triangle t^{-1}(\phi^{n+1}-\phi^{n}) +\rho\sigma\phi^{n}- \nabla \cdot \vec{u}^{n+1/2}= +\triangle t^{-1}(s^{n+}-s^{n})\delta(\vec{r}-\vec{r}_s) + + + +where \\(\\sigma = 2 \\frac{1-\f}{\f} \\triangle t^{-1}\\). + +Using mimetic finite volume approach, we spatially discretize above equations: + +.. math :: + + \mathbf{diag}(\Acf^T\mu^{-1})\frac{\mathbf{u}^{n+1/2}-\mathbf{u}^{n-1/2}}{\triangle t} + \mathbf{diag}(\Acf^T\mu^{-1}\sigma)\mathbf{u}^{n+1/2}= \mathbf{Grad} \mathbf{\phi}^{n} + + \mathbf{diag}(\rho)\frac{\mathbf{\phi}^{n+1}-\mathbf{\phi}^{n}}{\triangle t} + \mathbf{diag}(\rho\sigma)\mathbf{\phi}^n -\mathbf{Div} \mathbf{u}^{n+1/2}= +\mathbf{diag}(\mathbf{vol})\frac{\mathbf{s}^{n+1}-\mathbf{s}^{n}}{\triangle t} + +where \\(\\mathbf{Div}\\) and \\(\\mathbf{Grad}\\) and discrete differential operators and \\(\\Acf\\) is averaging operator from cell face to center. \\(\\mu\\), \\(\\rho\\) and \\(\\sigma\\) are defined on the cell center, and \\(\\mathbf{vol} \\) is volume of each cell. When we compute these, we first compute: + +.. math :: + + + \mathbf{u}^{n+1/2} = \mathbf{u}^{n-1/2} - \triangle t\mathbf{diag}(\sigma)\mathbf{u}^{n-1/2}+\triangle t\mathbf{diag}(\Acf^T\mu^{-1})^{-1}\mathbf{Grad} \mathbf{\phi}^{n} + +Then we compute: + +.. math :: + + \mathbf{\phi}^{n+1} = \mathbf{\phi}^{n} - \triangle t\mathbf{diag}(\sigma)\mathbf{\phi}^n + \mathbf{diag}(\rho)^{-1}\triangle t[\mathbf{Div}\mathbf{u}^{n+1/2}+\mathbf{diag}(\mathbf{vol})^{-1}\frac{\mathbf{s}^{n+1}-\mathbf{s}^{n}}{\triangle t}] + +.. note :: + + Choice of \\(\\sigma \\) in sponge boundary condition case can be expanded to PML case. + + .. math :: + + \sigma = 2 \frac{1-f}{f} \triangle t^{-1} + + And range of \\(\ f \\) is 0.98-1, which can be useful reference property. + PML boundary ============ +PML has two fundamental factors: a. matching the impedance and b. damping. These purposes can be realized by considering solution of Helmholtz equation on complex plane: + +.. math :: + + \tilde{x} = x + \frac{1}{\imath \omega}\int^x_0 \sigma^x(\xi) d\xi + + \partial x = \frac{\imath\omega}{\sigma^x+\imath\omega} \partial\tilde{x} + +By parameterizing the physical coordinate as + +.. math :: + + \tilde{x} = \tilde{x}(x) + + x<0 \ : \ \text{real axis} + + x>0 \ : \ \Im[\tilde{x}] < 0 \ (\text{decaying term}) + +Another core treatment of PML is decomposing \\(\\phi \\) as: + +.. math :: + + \phi_d = + \begin{bmatrix} + \phi^x \\[0.3em] + \phi^y \\[0.3em] + \phi^z + \end{bmatrix} + + \phi = [1,1,1]\phi_d = \phi^x + \phi^y +\phi^z + +Substituting those yields: + +.. math :: + + \imath\omega \mu^{-1}\u = + \begin{bmatrix} + \frac{1}{\sigma^x+\imath\omega} \frac{\partial \phi}{\partial x} \\[0.3em] + \frac{1}{\sigma^y+\imath\omega} \frac{\partial \phi}{\partial y} \\[0.3em] + \frac{1}{\sigma^z+\imath\omega} \frac{\partial \phi}{\partial z} + \end{bmatrix} + + \imath\omega\rho + \begin{bmatrix} + \phi^x \\[0.3em] + \phi^y \\[0.3em] + \phi^z + \end{bmatrix} + -\rho\imath\omega + \begin{bmatrix} + \frac{1}{\sigma^x+\imath\omega} \frac{\partial u^x}{\partial x} \\[0.3em] + \frac{1}{\sigma^y+\imath\omega} \frac{\partial u^y}{\partial y} \\[0.3em] + \frac{1}{\sigma^z+\imath\omega} \frac{\partial u^z}{\partial z} + \end{bmatrix} + = \frac{1}{3}\imath\omega \delta(\vec{r}-\vec{r}_s) + \begin{bmatrix} + s \\[0.3em] + s \\[0.3em] + s + \end{bmatrix} + +With some linear algebra: + +.. math :: + + \mu^{-1}(\imath\omega \u + \Sigma\u) = \grad\phi + + \imath\omega\rho + \begin{bmatrix} + \phi^x \\[0.3em] + \phi^y \\[0.3em] + \phi^z + \end{bmatrix} + -\rho + \begin{bmatrix} + \sigma^x\phi^x \\[0.3em] + \sigma^y\phi^y \\[0.3em] + \sigma^z\phi^z + \end{bmatrix} + + + \begin{bmatrix} + u^x \\[0.3em] + u^y \\[0.3em] + u^z + \end{bmatrix} + = \frac{1}{3}\imath\omega \delta(\vec{r}-\vec{r}_s) + \begin{bmatrix} + s \\[0.3em] + s \\[0.3em] + s + \end{bmatrix} + +where + +.. math :: + + \Sigma = + \begin{bmatrix} + \sigma^x & 0 & 0 \\[0.3em] + 0 & \sigma^y & 0 \\[0.3em] + 0 & 0 & \sigma^z + \end{bmatrix} + +In time domain we have: + +.. math:: + + \mu^{-1}(\frac{\partial \vec{u}}{\partial t} + \Sigma \u) = \nabla \phi + + \rho\frac{\partial \phi_d}{\partial t} +\rho\Sigma \phi_d - + \begin{bmatrix} + \frac{\partial u^x}{\partial x} \\[0.3em] + \frac{\partial u^y}{\partial y} \\[0.3em] + \frac{\partial u^z}{\partial z} + \end{bmatrix} + = \frac{1}{3}\imath\omega \delta(\vec{r}-\vec{r}_s) + \begin{bmatrix} + s \\[0.3em] + s \\[0.3em] + s + \end{bmatrix} + +We discretize above equations in both space and time: + +.. math :: + + \MfMui \triangle t^{-1} (\mathbf{u}^{n+1/2}-\mathbf{u}^{n-1/2}) + \MfMui\mathbf{\Sigma}^{f}\mathbf{u}^{n-1/2} - \mathbf{Grad}\mathbf{I}_d\phi_d^{n} = 0 + + \mathbf{\Omega}^{cc} \triangle t^{-1} (\phi_d^{n+1}-\phi_d^{n}) + \mathbf{\Omega}^{cc} \mathbf{\Sigma}^{cc} \phi_d - \mathbf{Div}_{vec} \mathbf{u} + = \triangle t^{-1}\mathbf{diag}(\mathbf{vol})^{-1}(\mathbf{s}_d^{n+1}-\mathbf{s}_d^{n}) + +where + +.. math :: + + \MfMui = \mathbf{diag}(\Acf^T\mu^{-1}), \ + \mathbf{\Sigma}^{f} = \mathbf{diag}(\Acf_{vec}^T + \begin{bmatrix} + \sigma^x \\[0.3em] + \sigma^y \\[0.3em] + \sigma^z + \end{bmatrix} + ) + + \mathbf{I}_d = [\mathbf{I}^{cc}, \mathbf{I}^{cc}, \mathbf{I}^{cc}], \ + \mathbf{\Omega}^{cc} = + \begin{bmatrix} + \rho & 0 & 0 \\[0.3em] + 0 & \rho & 0 \\[0.3em] + 0 & 0 & \rho + \end{bmatrix}, + + \mathbf{\Sigma}^{cc} = + \begin{bmatrix} + \sigma^x & 0 & 0 \\[0.3em] + 0 & \sigma^y & 0 \\[0.3em] + 0 & 0 & \sigma^z + \end{bmatrix} + + \mathbf{s}_d =\frac{1}{3} + \begin{bmatrix} + \mathbf{s} \\[0.3em] + \mathbf{s} \\[0.3em] + \mathbf{s} + \end{bmatrix} + +Similarly we first compute: + +.. math :: + + \mathbf{u}^{n+1/2} = \mathbf{u}^{n-1/2} - \triangle t \mathbf{\Sigma}^{f}\mathbf{u}^{n-1/2} + \triangle t \MfMui^{-1}\mathbf{Grad}\mathbf{I}_d\phi_d^{n} + +Then we compute: + +.. math :: + + \phi_d^{n+1} = \phi_d^{n} - \triangle t \mathbf{\Sigma}^{cc} \phi_d + \mathbf{\Omega}^{cc \ -1} \triangle t\mathbf{Div}_{vec} \mathbf{u} + +\mathbf{\Omega}^{cc \ -1}\mathbf{diag}(\mathbf{vol})^{-1}(\mathbf{s}_d^{n+1}-\mathbf{s}_d^{n}) .. raw:: html - :file: examples\refraction.html + :file: examples/center_pml.html +Stability conditions +==================== -Backgrounds -=========== +Stability of forward modeling have two fundamental factors: cell size and time step size ( \\(\\triangle t\\) ). First, we determine cell size based on the number of cell per wavelength ( \\(\ G\\) ): + +.. math :: + + G = \frac{\lambda}{\triangle x} \approx 16 + +where + +.. math :: + + \lambda = \frac{v_{min}}{f_{main}} + +Second, we determine \\(\\triangle t\\) + +.. math :: + + \triangle t = \frac{\triangle x}{v_{max}}c + +where \\(\ c \\) is a proper constant. Notebooks ========= +1. `How to run acoustic wave modeling in simpegSeis `_ + Examples shown are generated through this ipython notebook. diff --git a/docs/source/api_license.rst b/docs/source/api_license.rst new file mode 100644 index 0000000..7fc38d2 --- /dev/null +++ b/docs/source/api_license.rst @@ -0,0 +1,25 @@ +.. _api_license: + +License +******* + +The MIT License (MIT) + +Copyright (c) 2013-2014 SimPEG Developers + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software is furnished to do so, +subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/docs/source/examples/center_nobc.html b/docs/source/examples/center_nobc.html new file mode 100644 index 0000000..6dab49a --- /dev/null +++ b/docs/source/examples/center_nobc.html @@ -0,0 +1,10799 @@ +
+ + +
+ +
+ +
+ + + + + + + + + +
+ Once + Loop + Reflect +
+
+ + + +
diff --git a/docs/source/examples/center_pml.html b/docs/source/examples/center_pml.html new file mode 100644 index 0000000..869bea1 --- /dev/null +++ b/docs/source/examples/center_pml.html @@ -0,0 +1,9919 @@ +
+ + +
+ +
+ +
+ + + + + + + + + +
+ Once + Loop + Reflect +
+
+ + + +
diff --git a/docs/source/index.rst b/docs/source/index.rst index 632bebf..d90e857 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -4,10 +4,15 @@ :align: center SimPEG (Simulation and Parameter Estimation in Geophysics) is a python package for simulation and gradient based parameter estimation in the -context of geoscience applications. simpegSeis uses SimPEG as the framework for the forward modeling and inversion of seismic geophysical proble +context of geoscience applications. simpegSeis uses SimPEG as the framework for the forward modeling and inversion of seismic geophysical problem. We consider two fundamental governing equations for seismic wave propagation: acoustic and elastic wave equations in both time and frequency domains. We initially discretize those problems to compute forward responses. Obvious next step is to solve inverse problems +We welcome any people who want to contribute this project developing forward modeling and inversion package of seimic data, simpegSeis. + +.. raw:: html + :file: examples/refraction.html + Acoustic wave ============= @@ -19,9 +24,16 @@ Acoustic wave Elastic wave ============ +License +======= -Indices and tables -================== +.. toctree:: + :maxdepth: 2 + + api_license + +Project Index & Search +====================== * :ref:`genindex` * :ref:`modindex` diff --git a/notebooks/SeismicEx-sponge_explicit.ipynb b/notebooks/SeismicEx-sponge_explicit.ipynb new file mode 100644 index 0000000..e955049 --- /dev/null +++ b/notebooks/SeismicEx-sponge_explicit.ipynb @@ -0,0 +1,40153 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:5081835fb75e1b0331697bd22049838b1282f54319cf4228737a059be7421dcf" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from SimPEG import *\n", + "import scipy\n", + "%pylab inline" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "cs = 0.5\n", + "hx = np.ones(500)*cs\n", + "hy = np.ones(500)*cs\n", + "mesh = Mesh.TensorMesh([hx, hy], 'CC')" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Acoustic Wave equation in time domain (1st order form)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\rho\\frac{\\partial p}{\\partial t} = \\nabla \\cdot \\vec{u}+\\frac{\\partial s}{\\partial t}\\delta(\\vec{r}-\\vec{r}_s)$\n", + "\n", + "$ \\mu^{-1}\\frac{\\partial \\vec{u}}{\\partial t} = \\nabla p$" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Add damping term" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\rho[\\frac{\\partial p}{\\partial t} + \\sigma p] = \\nabla \\cdot \\vec{u}+\\frac{\\partial s}{\\partial t}\\delta(\\vec{r}-\\vec{r}_s)$\n", + "\n", + "$ \\mu^{-1}[\\frac{\\partial \\vec{u}}{\\partial t}+\\sigma \\vec{u}] = \\nabla p$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discretized form" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\mathbf{diag}(\\rho)\\frac{\\mathbf{p}^{n+1}-\\mathbf{p}^{n}}{\\triangle t} + \\mathbf{diag}(\\rho\\sigma)\\mathbf{p}^n = \\mathbf{Div} \\mathbf{u}^{n+1/2}+\\mathbf{M}^{cc -1}\\frac{\\mathbf{s}^{n+1}-\\mathbf{s}^{n}}{\\triangle t}$ \n", + "\n", + "\n", + "$ \\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1})\\frac{\\mathbf{u}^{n+1/2}-\\mathbf{u}^{n-1/2}}{\\triangle t} + \\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1}\\sigma)\\mathbf{u}^{n+1/2}= \\mathbf{Grad} \\mathbf{p}^{n}$" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Compute p first" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\mathbf{p}^{n+1} = \\mathbf{p}^{n} - \\triangle t\\mathbf{diag}(\\sigma)\\mathbf{p}^n + \\mathbf{diag}(\\rho)^{-1}\\triangle t[\\mathbf{Div} \\mathbf{u}^{n+1/2}+\\mathbf{M}^{cc -1}\\frac{\\mathbf{s}^{n+1}-\\mathbf{s}^{n}}{\\triangle t}]$ \n", + "\n" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Then compute u" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "$ \\mathbf{u}^{n+1/2} = \\mathbf{u}^{n-1/2} - \\triangle t\\mathbf{diag}(\\sigma)\\mathbf{u}^{n-1/2}+\\triangle t\\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1})^{-1}\\mathbf{Grad} \\mathbf{p}^{n}$" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "mesh.setCellGradBC('dirichlet')\n", + "Grad = mesh.cellGrad\n", + "Div = mesh.faceDiv\n", + "rho = np.ones(mesh.nC)*2.7\n", + "vhalf = 2000\n", + "vblk = 2800\n", + "v = np.ones(mesh.nC)*vhalf\n", + "blkind = np.logical_and(mesh.gridCC[:,1]>-12.5, mesh.gridCC[:,1]<12.5) & np.logical_and(mesh.gridCC[:,0]>-12.5, mesh.gridCC[:,0]<12.5)\n", + "v[blkind] = vblk\n", + "mu = rho*v**2\n", + "AvF2CC = mesh.aveF2CC" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 19 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "dat = mesh.plotImage(v)\n", + "plt.colorbar(dat[0])" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 20, + "text": [ + "" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEPCAYAAABFpK+YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHN1JREFUeJzt3X2wXOV92PHvtbiAjHAwYw/ojV6RiBYUCLKMwJETBDVE\ndFJeyiSYFExaxmYqDAp1zYvcFhE3LiY1INqipIMoCAVSKtkMHl4MZiQMGUCISqAXFIQjMehawhhI\nhAJGL9z+8XvW99y9u3t3de7VnnP3+5k5s2efc3b3ORfxO8/+zrPnB5IkSZIkSZIkSZIkSZIkSZIk\ndaTJwApgA7AeuDq1zwRWAWuAF4FTMq+5AdgMbALOzrTPANalbQtHtNeSpP12NHByWh8H/C1wPLAS\n+L3Ufg5xcgA4AVgLdAM9wOtAV9q2ijhhADwKzKn3oZ8Yjp5LkvbLDiKQA+wCXgUmAtuBX0vtRwC9\naf084AFgD7CVCPynAuOBw4ngD7AEOL/ehx40XL2XJOXSA0wHnifSNc8C/40YoH8h7TMhba/YRpwo\n9qT1it7UXpMjfklqv3HAMmAeMfJfTOT7jwGuAe4ezg8bhSP+f9IHb7S7E5LK4Wlg9v6++FDo+2Vr\nL3kf+FRVWzewHFgKPJTaZgJfSuvLgLvSei9xQbhiEjHS703r2fZe6hiFgf8NYMEB+qwVwBkH6LMO\nlNF4TDA6j2s0HhMc2ONacHqeV/8S+C8t7P8fIw+f1UWM7jcCt2faXwdOJ05MZwKvpfaHgfuBW4lU\nzlQir98H7CTy/auAS4E76vVjFAZ+STpwuvO9fBZwCfAKMXUTYD7wNeB/AocAH6bnECeIB9PjXmAu\nEfRJ6/cAY4lZPY/X+1ADvyTlkDOIPkv9a62n1mn/TlqqvQSc2MyHGvhz6Wl3B0ZAT7s7MEJ62t2B\nEdDT7g6MkJ52d6AlY9vdgf1g4M9lSrs7MAJG4zHB6Dyu0XhMULbjypnqaQsDvyTlUMYgWsY+S1Jh\nOOKXpA5TxiBaxj5LUmE44pekDmPgl6QO43ROSeowZQyiZeyzJBWGqR5J6jBlDKJl7LMkFYYjfknq\nMGUMomXssyQVRhlH/O0uvXg38BawLtN2JPAkUXjgCaLQcMUNRC3KTcDZB6iPklTX2BaWGiYTlWc2\nAOuJcosA/4e4P/8aYAv99+o/C1hN3L9/NQMr1swgYulmYGGjPrc78P9vYE5V2/VE4D8OeCo9BzgB\nuCg9zgHupP39l9ThultYathD1NSdBpwGXAkcT8S66WlZnhaAt4HfB04CLgPuy7zXIuByoirXVAbH\n1l9pd+B8Bnivqu1c4N60fi9wflo/D3iA+ENtJUqTzRz5LkpSfQe1sNSwA1ib1ncBrwITMtu7gD8k\nYh9p3x1pfSPxRaIbGE+UdVyVti2hP3bW7HPRHEWkf0iPR6X1CcDzmf22ETUnJaltuluJonsbbu0h\nRvgvZNp+h4iDP62x/4VE1a09RCzcltnWS4P4WMTAn9VHfz3JettrWJFZ76FshR0kjZQtRMJg+BzU\nIIo+sw+e/biptxkHLAPmESP/iouJ4urVpgE3Ezn/lhUx8L8FHE18nRkP/Dy19xIXQiompbYazqjd\nLKnDTWHgQPDp3O/YPab+tjPHwJmZ5zf/Y+23IHL4S4GHMu0HARcAn6vafxLwfeBS4kwGEQsnVe1T\nJz62P8dfy8PERQvS40OZ9i8DBxP/5abSn8+SpLY46KDmlxq6gMVEvv72qm1fInL+P8u0HQE8AlwH\nPJdp3w7sJAq0dxEnhexJZGCfWzrC4fcAcDrwGeBN4D8TX18eJK5ObyUubED8YR5Mj3uBuTROA0nS\niOs+JNfLZwGXENMzK1M2bwAeJ2b2PFC1/9eBXwduTAtEuucXREy8h7jg+2h6j5q6cnW5mPpgQbv7\nIKkUFkC+ONjXN2HonSq6Yuze9rjb7hG/JJVbCaNoCbssSQVSwihawi5LUoE0mNVTVAZ+ScqjhFG0\nhF2WpALJN6unLQz8kpRHCaNoCbssSQVSwihawi5LUoF4cVeSOkwJo2gJuyxJBVLCKFrCLktSgZQw\nipawy5JUIE7nlKQOU8IoWsT78UtSeYxpYRlsMlEycAOwHrg6s+0q4n7864HvVr3uGKJS1zcybTOA\ndcBmYGGjLpfwXCVJBZIviu4BriGKqI8jaug+SVQhPBc4Ke3z2arX3UoUZMlaRNQxWUXcj38Ode7J\nb+CXpDzyRdEdaYEYwb9KFEn/KvBfiaAP8HbmNecDfwdkCzmOBw6nvyrhkrRfzcBvqkeS8siX6snq\nAaYDLwDHAb8LPA+sBD6f9hkHXMvgalMTgW2Z572prSZH/JKUR4MouvJnsHJ7U+8yDlgGzAPeT+/6\naeA04BSi7OyxRMC/DfiAHJW8DPySlMeh9TfNPjaWipvW1NytG1gOLKW/QPo24Ptp/UXgY6I2+Uzg\nQuAWovD6x8CHad9JmfecRIz6azLwS1Ie+e7V0wUsBjYCt2faHwLOBJ4m0j4HEwXVfzezz43Et4M7\n0/OdwKlEnv9S4I56H2rgl6Q88kXRWcAlwCtA5fvADcDdaVkH7Aa+0sR7zQXuAcYSs3pqXtgFA78k\n5ZMvij5L/Uk2lw7x2puqnr8EnNjMhxr4JSkPb8ssSR2mhFG0hF2WpAIpYRQtYZclqUC8O6ckdZgS\nRtESdlmSCqSEUbSEXZakAnFWjyR1mBJG0RJ2WZIKpIRRtIRdlqQCMdUjSR2mwd05i8rAL0l5lDCK\nlrDLklQgJUz1FLn04lb6b1VaqSN5JFGI+DXgCaIQgSS1z0EtLINNBlYAG4D1wNWpfQFRjGVNWuZk\nXnMS8Fza/xXiXv0AM4jbOG8GFjbqcpEDfx8wm6hBOTO1XU8E/uOAp9JzSWqffIF/D3ANMI0os3gl\ncDwR/24l4t90+u+tfxBwH/A14DeB04G9adsi4HJgalqyJ4sBihz4YXBNyXOBe9P6vUQVeUlqn3zF\n1ncAa9P6LuBV+ouk16qpezYxyl+Xnr9HlF8cDxxOf3ZkCQ3iY5EDfx/wY2A18NXUdhTwVlp/Kz2X\npPY5tIWlsR5idP98en4V8DJRmrGS1p5KxMbHicIr30ztE4nUUEUv/SeQQYp8cXcWsB34LJHe2VS1\nvS8tktQ+w3NxdxywDJhHjPwXAX+atn0b+B6RxukGvgh8niiy/hRxAviHVj6syIF/e3p8G/gBked/\nCzia+Ho0Hvh57ZeuyKz3AFNGqIuSymULMW9kGDWIoitfgpX/b8h36AaWA0uJIuswMLbdBfwwrb8J\n/AR4Nz1/FPhceu2kzGsmEaP+mmrlkIrgk8R59H3gMGIGz03Al4B3gO8SF3aPYPAF3r64IC5JQ1kA\n+eJgX9/q5nfu+jzVn9dFXK98h7jIWzGe/sHvNcApwB8BnyZS4F8kLgw/RlwEfgx4gZgVtAp4BLiD\nOgXXizriP4oY5UP08a+I4L8aeJD4yrMV+MN2dE6SfiVfFJ0FXEL/1HWA+cDFwMlEOnsLcEXa9h4R\n6F9M2x4hgj7AXOAeYCzxTaBm0IfijvjzcMQvqUkLIO+If93QO1V0nUjezxsWRR3xS1I5lDCKlrDL\nklQg1tyVpA5Twihawi5LUoGUMIqWsMuSVCAljKIl7LIkFUdfCW/LbOCXpBz2lTCKlrDLklQcBn5J\n6jAfHXLw0Dv9yu4R60crDPySlMO+MeVL8hv4JSmHfSUsumvgl6Qc9hr4Jamz7CthGC1y6UVJKrx9\njGl6qWEyUTlqA7CeuJ9+1jeImrpHpueHAg8Qt3HeyMB6JDOIWrybgYWN+mzgl6Qccgb+PUShlWnA\nacCVwPFp22TgLOCNzP5fTo8nEYH+CuCY1LaIqFUyNS1z6vXZwC9JOXzEwU0vNewA1qb1XcCrwIT0\n/Fbg2qr9txNVCcekx93ATqJi1+FE9S2AJcD59fpcvuSUJBXIMOb4e4DpRAnF84BtREon60fApcQJ\n4JPAnwB/D/xG2r+iF5hY74MM/JKUwzBN5xwHLAPmETn9+USap6JStesSorTieCLv/wzwVKsfZuCX\npBwaBf7VK/+R1Ss/GOotuoHlwFLgIeBEYvT/cto+CXgJOBX4baIe+T7gbeBviFz/s2k/Mq/prfeB\nBn5JyqHRPP6TZ3+Kk2d/6lfP/9dNv6jepQtYTMzQuT21rQOOyuyzhQju7wKbgDOJk8RhxAXh24hr\nBTuJk8MqIh10R71+eXFXknLYx0FNLzXMItI3ZwBr0nJO1T59mfW/BA4mTg6rgLuJaaAAc4G7iOmc\nrwOP1+uzI35JyiFnjv9Zhh6AH5tZ/4g4UdTyEpEmGpKBX5Jy2F17mmahGfglKQfv1SNJHaaM9+op\nX48lqUC8LbMkdRgDvyR1GHP8ktRhdnNIu7vQMgO/JOVgqkcaYQtYMCo/S+VlqkeSOozTOSWpw5jq\nkaQOY+CXpA5j4JekDvNRCadzlvF+/HOIYgSbgeva3BdJHW4fY5peapgMrAA2EPfVv7pq+zeIUoxH\nZtpuIOLfJuDsTPsM4j79m4GFjfpctsA/BvgfRPA/AbgYOL6tPZLU0XIG/j3ANcA0oprWlfTHtMlE\n3d03MvufAFyUHucAd9Jfj3cRcDkwNS1z6vW5mcB/NfDpJvY7EGYSlWW2En+wvyaq0UtSW+xlTNNL\nDTuAtWl9F/AqMCE9vxW4tmr/84AHiPi3lYiHpxLF1w8nqnIBLAHOr9fnZgL/UcCLwIPEGaSr8e4j\naiLwZub5ttQmSW2Rs/RiVg8wHXiBCPDbgFeq9pmQ2isqMbC6vZcGsbGZi7vfAv4TkUv6YyLV8iBR\nIPinTbx+OPUNvQtEyqyiB5gyAl2RVD5biIHy8Gk0q2fryjd4Y+UbdbdnjAOWAfOInP58Is1TMawD\n7mZn9XxMfCV5C9hHpH6WAT8GvjmcHRpCL5H3qpjMwLNccsYB6o6kcpnCwIHg07nfsVHgnzz7WCbP\n7i+Z+5Obnq21WzewHFgKPETUze0BXk7bJxH1dE9lcAycRMTA3rSebe+t169mAv884CvAO0QF9/9A\n5Jc+QVw9PpCBfzVx0aIH+BlxkePiA/j5kjTAR/lq7nYR2ZONwO2pbR2RYq/YQszYeRd4GLifyP9P\nJOLhKiIbspM4OawCLgXuqPehzQT+I4F/xcAryxDfAv5lE68fTnuBrwM/Imb4LCYuhkhSW+S8V88s\n4BIil78mtc0HHsvsk01xbyRS7RuJeDg3s30ucA8wFngUeLzehzbT4xsbbNvYxOuH22MM/KNIUtvk\n/OXusww9yebYquffSUu1l4g00ZD85a4k5eAtGySpw3g/fknqMN6PX5I6jKkeSeowu/NN52wLA78k\n5WCOX5I6jDl+aYQtYEG7uyANYI5fkjqMgV+SOow5fknqMOb4JanDOJ1TkjpMGVM9ZSu2LkmFkrP0\n4mSiZOAGYD1R4xzg20QhlrXAU/QXXzmLqEvySnrMVp2aQdzLfzOwsFGfDfySlMM+xjS91LAHuAaY\nBpwGXAkcD9wC/BZwMlGVq3J7/LeB3wdOAi4D7su81yLgcqI4y1SiRnpNpnokKYec0zl3pAVgF1FY\nagIDC0yNA36R1tdm2jcSRVe6gc8AhxPVtwCWAOdTpxiLgV+SchjGefw9wHTghfT8z4gSih8Q3waq\nXUgUX9lDlGHM1h/vTW01meqRpBw+4pCmlwbGAcuIGue7Utu3gGOIcoq3Ve0/DbgZuGJ/+uyIX5Jy\naDTi/2Dli3ywcvVQb9ENLAeWEvn8avcTNXQrJgHfJ74NbEltvak9u09vvQ808EtSDo0C/yGzT+OQ\n2f1Zmndv+ovqXbqAxUS+/vZM+1Ridg7AefQXYj8CeAS4Dngus/92YCdwKpHnvxS4o16/DPySlEPO\nefyzgEuI6ZmV4D6fmJ3zT4F9wE+Bf5e2fR34dWKWT2Wmz1nExd+5RFpoLPENoeaFXYizzWjTh3dw\nlNSUBZAvDvZN6ts89F7Jtq6peT9vWDjil6QcvDunJHUYA78kdZiPdnuTNknqKPv2li+Mlq/HklQg\n+/aa6pGkjmLgl6QOs3ePgV+SOsrH+8oXRsvXY0kqElM9ktRhflm+MFq+HktSkextdwdaZ+CXpDwM\n/JLUYUoY+ItYgWsBUUJsTVrOyWy7gbhH9Sbg7APeM0mqtqeFZbDJwApgA7AeuDq1/zlRd/dloujK\nr1W97hiiUtc3Mm0zgHVEjFzYqMtFDPx9wK1E7cnpwGOp/QTgovQ4B7iTYvZfUifZ18Iy2B7gGqKU\n4mnAlcDxwBOp7beA14hBb9atREGWrEXEffynpmVOvS4XNXDWul/1ecADxB9qK/A6MPMA9kmSBtvb\nwjLYDmBtWt9FjPInAE8CH6f2FxhYVvF84O+Iql0V44HDiepbAEvSfjUVNfBfRXzFWUyUGoP4Y2Sr\nyG+jQRV5STogftnC0lgPkeV4oar939Jfc3cccC2Dq01NZGB87KVBfGxX4H+SyEVVL+cSX1emACcT\ndSS/1+B9+ka2m5I0hHwj/opxwDJgHjHyr/gWsJsouA4R8G8DPiBHJa92zeo5q8n97gJ+mNZ7iQsh\nFQ2qyK/IrPcQ5xFJ2kJkiodRo4C+biWsXznUO3QDy4GlwEOZ9j8G/gXwzzNtM4ELgVuIbMjHwIfE\nBeBsOqhBfCxA7ccaxhMjfYiLHqcAf0Rc1L2fOPCJwI+B32DwqN+au5KatABy1txleQuJhwu7qj+v\nC7gXeIeIdxVziGzH6UQh9VpuBN4nLvRCpIiuJvL8jwB3UKfgehHn8X+XSPP0EafnK1L7RuDB9LiX\nqChvqkdSe9WeptmsWcAlwCvE9HWA+UTQPphIiwM8R8S8RuYC9wBjiWsCNYM+FHPEn5cjfklNWgB5\nR/x/1cL4818PGvG3RRFH/JJUHiX85a6BX5LyGHqaZuEY+CUpD0f8ktRhDPyS1GEM/JLUYfJN52wL\nA78k5VH7rpuFZuCXpDyc1SNJHcYcvyR1GHP8ktRhzPFLUocx1SNJHcbAL0kdpoQ5/qLW3JWkcvio\nhWWwyUTJwA3AeqKQCsAfpLZ9wOeqXnMScX/+9cR9/A9O7TOIErabgYWNumzgl6Q88tXc3UNU3poG\nnAZcCRxPBPALgJ9U7X8QcB/wNeA3iQpdlXdeBFwOTE3LnHpdNtUjSXnkS/XsSAtEkfVXgQnAU3X2\nP5sY5a9Lz99Lj+OBw4myiwBLgPOpU4XLEb8k5bGvhaWxHmA6UTu3nqlEydnHgZeAb6b2icC2zH69\nqa0mR/ySlEejWT2/WAnvrGzmXcYBy4B5xMi/nm7gi8DngQ+JbwYvAf/QzIdUGPglKY9Ggf+I2bFU\nvHZTrb26geXAUuChIT7tTSLv/256/ihx8XcpMCmz3yRi1F+TqR5JymNPC8tgXcBiYCNwe51PyBZn\n/xFwIjCWGLifTsz+2QHsBE5N+19Kg5OII35JyqP2NM1mzQIuIS7Yrklt84FDgP8OfAZ4JG07B/h7\n4FbgRSLX/wjwWHrdXOAe4qTwKHUu7MLAM8lo0QcL2t0HSaWwAPLFwT6+0Nf83s915f28YeGIX5Ly\nKOEvdw38kpSHd+eUpA7jTdokqcMY+CWpw5jjl6QOk286Z1sY+CUpD1M9ktRhTPVIUodxOqckdRhT\nPZLUYQz8ktRhzPFLUocp4Yi/Xffjb1RB/gaiSvwmor5kRdMV5CWpJCYDK4h4uB64OrUfCTwJvAY8\nARyR2g8FHiBu47wRuD7zXk3HyHYF/noV5E8ALkqPc4A76b+FadMV5CWpJPYA1wDTgNOAK4HjiYD+\nJHAcUV6xEuC/nB5PIgL9FcAxqa3pGNmuwL+JOJNVO484m+0BtgKvExVl6lWQl6Qy2wGsTeu7gFeJ\nIunnAvem9nvpj3fbgcOAMelxN1F5q6UYWbTSixMYWCl+G/FHqG5vWEFekg6cfLUXM3qA6cALwFHA\nW6n9rfQcovTiTuIEsBX4c6Iq10RaiJEjeXH3SeDoGu3zgR+O4OcSKbOKHmDKyH6cpJLYQsTL4dTo\n6u5PGJzRrmkcUXB9HvB+1ba+tECUaRxLjPCPBJ4hUkEtGcnAf9Z+vKaXuNhRMYk4i/XSQgV5OGM/\nPlrS6DeFgQPBp4fhPRuN5L+Qlorv1Nqpmwj699FfIP0tYuC8gwjyP0/tvw38gJgY8zbwN0Su/1la\niJFFSPVk608+TFy8OJj4rzOVyFm1VEFekg6cD1tYBukCFhMzdG7PtD8MXJbWL6M/3m0CzkzrhxEX\nhDfRYoxsV+C/AHiT6HS2SvxG4MH0+BhRNb7yFWcucBcxVel1GlSQl6QDJ1eOfxaRvjkDWJOWOcDN\nRNbkNSLQ35z2/0tiYLyOGBTfTUwDhRZiZNurvY+APljQ7j5IKoUFkC8O9sV1g2ZNyft5w8Jf7kpS\nLuW7Z4OBX5JyKd89Gwz8kpSLI35J6jA1Z+sUmoFfknIx1SNJHcZUjyR1GEf8ktRhHPFLUodxxC9J\nHcYRvyR1GKdzSlKHccQvSR2mfDn+ItyPX5JKLNdtmScTJQM3ELdXvjq1H0lUMXwNeAI4IvOaG4hb\nL28Czs60zyBu17wZWNioxwb+XFq5HWtZjMZjgtF5XKPxmKB8x7W3hWWQPcA1wDSiPsmVwPHA9UTg\nP44orXh92v8E4KL0OAe4k/7bPC8CLicKWE1N22sy8Oeytd0dGAFb292BEbK13R0YAVvb3YERsrXd\nHWhRrhH/DmBtWt8FvEoUST8XuDe13wucn9bPAx5Ib7aVKLhyKlGe8XCiOAvAksxrBjHHL0m5DFuO\nvweYDrwAHEXU3SU9HpXWJwDPZ16zjThR7EnrFb2pvSYDvyTlMizTOccRBdfnAe9XbeujvwSt6lhJ\n/x/KxcXFpdGyknxa/bydNd6jG/gR8CeZtk3A0Wl9fHoOkeu/PrPf40Sq52giTVRxMfAX+3lMkqQR\n1EXk42+rar8FuC6tX09/sfUTiGsCBxMFfH9K/8XdF4iTQBfwKA0u7kqS2ueLwMdEMF+TljnEdM4f\nU3s653ziou4m4Pcy7ZXpnK8Dd4x0xyVJGnX+gPiBxT7gc1Xbcv+YoiAWELMCKqOOczLb6h1jGcwh\n+r2Z/q/OZbUVeIX471OZttfohz5FdDcxS2Vdpm1/fqwkjbh/RvyQYgUDA38l39ZNTMV6nf582ypg\nZlovQ77tRuDf12ivdYxl+f3HGKK/PUT/1xI/jimrLUSQzLoFuDatX0d/LriofoeYspgN/PWOocz/\n9grNP2JzNhGjkWrD8mOKAumq0VbrGGfW2K+IZhL93Ur0/6+J4ymz6v9G9X7oU1TPAO9VtbXyY6Wy\n/NsrNAN/PhMY+KOJyo8pqtsb/piiQK4CXgYW0/91u94xlsFE4M3M8zL1vZY+4oLfauCrqa3eD33K\npNGPlcr6b6/Q/AFXvyfpnzebNR/44QHuy0ipd4zfIu7z8afp+beB7xH3/ailb/i7NiLK0s9mzQK2\nA58l/ltuqtpemSteZkMdQ9mPrxAM/P3O2o/X9BJ316uYRIxKetN6tr13/7s2bJo9xrvoP9nVOsYi\nHEszqvs+mYEjyLLZnh7fBn5ApD3eIk7mO4gU48/b07Vc6h1Dmf/tFZqpntZlc6wPA1+m/8cUU4m8\n/g7iF3qVH1NcCjx0YLvZsvGZ9Qvov/hW7xjLYDXR3x6i/xcRx1NGnySuGwEcRsxwWUccz2Wp/TKK\n/++slnrHUOZ/exoFLiByxR8SQf2xzLbR8mOKJcRUwZeJ//GyueJ6x1gG5wB/S/T/hjb3JY8pxAyX\ntcR92yvH0uiHPkX0APAzYDfx/9S/Yf9+rCRJkiRJkiRJkiRJkiRJkiRJkiRJkqQD6RTil8iHELc5\nWE/c413qaLXuvy6NJt8GDgXGErcI+G57uyNJGmndxKj/eRzoSIB359To9xkizTOOGPVLHc8RkEa7\nh4H7gWOJW09f1d7uSJJG0leA/5vWP0Gke2a3rTeSJEmSJEmSJEmSJEmSJEmSJEmSJGl0+/9BB90c\njwdxcwAAAABJRU5ErkJggg==\n", + "text": [ + "" + ] + } + ], + "prompt_number": 20 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "time = np.linspace(0, 0.08, 2**10)\n", + "dt = time[1]-time[0]" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 21 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print dt" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "7.82013685239e-05\n" + ] + } + ], + "prompt_number": 22 + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Determine $\\sigma$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\sigma = 2 \\frac{1-f}{f} \\triangle t^{-1} $" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "ax = mesh.vectorCCx[-30]\n", + "ay = mesh.vectorCCy[-30]\n", + "indy = np.logical_or(mesh.gridCC[:,1]<=-ay, mesh.gridCC[:,1]>=ay)\n", + "indx = np.logical_or(mesh.gridCC[:,0]<=-ax, mesh.gridCC[:,0]>=ax)\n", + "tempx = zeros_like(mesh.gridCC[:,0])\n", + "tempx[indx] = (abs(mesh.gridCC[:,0][indx])-ax)**2\n", + "tempx[indx] = tempx[indx]-tempx[indx].min()\n", + "tempx[indx] = tempx[indx]/tempx[indx].max()\n", + "tempy = zeros_like(mesh.gridCC[:,1])\n", + "tempy[indy] = (abs(mesh.gridCC[:,1][indy])-ay)**2\n", + "tempy[indy] = tempy[indy]-tempy[indy].min()\n", + "tempy[indy] = tempy[indy]/tempy[indy].max()\n", + "temp = tempx+tempy\n", + "temp[temp>1.] = 1.\n", + "f = 1-temp*0.1\n", + "sig = (1.-f)/f*2./dt" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 23 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "dat = mesh.plotImage(sig)\n", + "plt.colorbar(dat[0])" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 24, + "text": [ + "" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X98VPWd7/EXJBADBIEAgUA0KD9qFliFSrxwq7Mtssjj\nVqX7UNStpV26u48H96qP3h9W7N0ttrve6m691d0rd7f+AFqh67XVxceiFXwYVFAoLAgYkIBESQKJ\nBBDQGEjI/ePzTWcYZoaZnJnMOTPv5+NxmDPfc+bM9xuSz5z5nO/5fkFERERERERERERERERERERE\nRERERLKoT7YrkG6XQ9dH2a6EiATFBiDU0xdfAl1fpPaS48Cwnr5fuuRc4Ae6libYeDlQCYRKgQlA\nBVAKDEz9jZZugqUzU3+dn+VimyA325WLbYIU29UBfAocBj6E+jqox6J5Uu9lD17iYNffpLDz//T+\nfmlRmO0K9JYS4E+BshuAcuBq4Cq3PgwY1IODPgrcn64a+kQutglys1252CZIrV3dgb8F+BAqd0Nl\nE4Q2w/oPYWPGKhnWrxfeI93yIvBPA74+AZgHuMD/2dV9qS8aRwsjOc4QzlBEJwUpHXfnwFqeG1GV\ngRpnTy62CXKzXbnYJkitXQV0UjL6FKVfaqXi+kOMOdQKTcBmmP1buGotPJXZ6gYyiAaxzimZBcy+\nDrgNmAMfTx5JMyOppYp6LPCfYAjtFKV87NbQJF4kt/7wcrFNkJvtysU2QWrtKqCDAbRRylEqOMT4\nigOMrGjmmqt2UzgOxpTD/avh0c8yV9/izB06Y7Kea8qA3+f4vwuM+S5wJ3z81ZFsYiZbmU4LZXzA\nJOqppKWpDI4WQRvQmbU6i0hPXAIMguJRxykf3MSVHKCMZqawi2o2c/3hLfAb4ClYvgOiO34stQdP\nOf6fp7Dzn3t/v7TI2TP+WcCYecDN0PjVUt4gxCZmso0v00wZDfvGw36gATgKpHhpXkR8oBAYBG2j\nhnKgcihNk8spHdzqvsX3h9Fw/bwt0ALfboL1LenP+3sMohXASmAk0AX8M/AE8C/ARLfPEOAEcI17\nvgT4M+xU9V7gNVc+HViOfRyuBe7LUJ39aTYw6wrgj6HjBtjMDLZQzRaq2dF0DRwqgveAOqAR+ARo\nz2aNRaRHXOCnDBgHbUeH0lA2lDPV/QEo4gxDxp1g6tf2QR3M3gwbP0xvFTxe3D0LfA/YgbVkG7AO\nWBCxz99jgR+gym2rAsYA67H+iV3AMmARsAUL/HOBV2O9aU4G/lmXAdW2bB88mR1cw06msLN5CrxV\nBM3Yj7n7jP8EcCqLFRaRnumHnQ8Px/6uTwOjoIXL2F59hiGcoJwmKmfVM3jXGQBu+DD57p7J8BhE\nj7gFrPZ7sL6Ge1xZH+B24I/c81uA1dgHRj0WxaqxLFYJFvTBvkXcSj4FfuYAs6CxupRNzGQz1Ww5\nVs251wfC21iw3w0cBDgGtKLILxJAHYXQNgwOl0F9P+vPPxo4DQ2d49k0s50STjGyoJkbv/Y2DITQ\nVthQl74qpLE7ZyWWztkcUfYV7CPtgHteDrwbsb0BO/M/69a7NbrymHIz8F9hyyEqOEQFTYzmTN1g\n+3ysx3489WA/z3os+LdhPzsRCY4B2N/vKTg13oJ/O5b6qYSmyeU0DS7nEJfx2bi+DGw5B1fA5XUX\nXujtqTQF/kHAC1he/nRE+Z3AqvS8RVhuBv6RtrRSSivDaWW4ndR3f6k6CnS1YR+KLdgHQFvWqisi\nPVUIDHbrJdA61iJxM3AE2o4M5ejgUosFRcMZOLIFRsNY0hf4E3XnfM8tF9EP+DXwS+CliPJCYD52\nK1K3RuyCcLex2Klso1uPLG+M94a5GfgvteUUJXxOMac+L7E8/mkso3MC4KRbWt2jAr9I8HSfbw8G\njkHX2PDfultOU8Ipt3BpC1ya3sFyEgXR6W7p9osLd+kDPA3UAj+L2jYby/U3RZStwb4BPIalciZg\nef0uLJBVu+d3Y72DUq5zcBUBhdBJAZ0U0tlRYLd2f4F1gDqL+6cjahGR4OkAPnePZ+Fsv/Cf9BfQ\nQYGLBQUW8QrTe9OVx1TPLOCbwE5guytbgl2UXYBdyI1UCzzvHjuAxVjQx60vx5q3ljgXdiFXA7+I\nSC/xGETfBvrG2fadOOUPuyXaNmBKMm+qwC8i4oEGaRMRyTNBDKJBrLOIiG8E8Yw/Xm6ptzyDdbza\nFVE2DLtleR82BsWQiG1LsIEW9mK3aYmIZFVxCotfZDvwP4uNJxHpASzwTwRed8/h/DEq5gJPkv36\ni0ie65fC4hfZDpxvYXNQRroZWOHWV2DjTUDsMSpmZL6KIiLxFaaw+IWf6tKtDEv/4B7L3Hq8MSpE\nRLKmXypR1Ce3C/kx8EfqInxzQrztF1i6CnjTpnBrDU2CafMyUjkRCZaajVCzyUZifz9NxyxU4E+L\nZmAUNqrOaGwwHYg9RkXMsSiW3gVcB8+NqLIp3E5msLYiEhihWRCaCXvesdtf0zE8c7/Upur2hWzn\n+GNZAyx06wsJD1q0BrgD6A+MIzxGhYhI1hQWJr/4Rbarshq4AZtG4RDw18BPsA/jRdhF3NvdvonG\nqBARyYp+RdmuQeqyHfjvjFM+O055vDEqRESyI9tRtAcCWGURER8JYBQNYJVFRHwkgFE0gFUWEfER\n9eoREckz3m7drQDewG4r2A3cG7X9vwHnOH/SsHhjlk3Hxj2rAx6/WJVFRKSnvPXqOQt8D9iBTbi+\nDRurbA/2oXAj508PHDlm2RhgPda1vQtYhvWG3ILNwDWXOLNw6YxfRMQLb2f8R7CgDzZL8B5seBqw\neXXvj9o/1phl1djNriWE721aSXics5hVFhGRnkpfFK0ErgE2YwG+AZuLN1K8McvOuvVujSQYy0yB\nX0TEiwQXd2s+tSUJg4AXgPuwnP6DWJqnW58e1y8GBX4RES8SRNFQqS3dHjoUc7d+wK+BX2JD1EzB\nzv7fc9vHYrn/amKPWdbgysdGlcccywyU4xcR8cZbjr8P8DQ2FM3PXNkubDj6cW5pAKZhA1jGG7Ps\nCDYcZbU75t2ExzmLWWUREekpb1F0FvBNLJe/3ZU9CLwSsU/kmGSJxixbDCzHZnlcS5wePd6rLCKS\n77x153ybi2deroh6Hm/Msm1YmuiiFPhFRLwIYBQNYJVFRHwkgEM2KPCLiHgRwCgawCqLiPhIAKNo\nAKssIuIjSvWIiOSZAEbRAFZZRMRHLsl2BVKnwC8i4oVSPSIieSaAUTSAVRYR8ZEARtEAVllExEeU\n6hERyTMBjKIBrLKIiI8EMIoGsMoiIj7ibXTOrNBELCIiXnibiKUCeAN4H9gN3OvKb3NlndgkLJGW\nAHXAXmBORPl0bBKXOuDxRFVW4BcR8cJb4D8LfA/4A+A64D8DV2EBfD7wZtT+VcAC9zgXeJLwfLzL\ngEXYrFwT3PaYFPhFRLwoSGG50BFgh1s/DewByrGz+X0x9r8FWI19YNQD+7HpFkcDJdg0jAArgVvj\nVVk5fhERL9IXRSuBa4DNCfYpB96NeN4AjME+CBoiyhtdeUwK/CIiXqQnig4CXgDuw878M0qBX0TE\niwQ3cNXsteUi+gG/Bn4JvHSRfRuxC8LdxmJn+o1uPbK8Md5BFPhFRLxIMDpn6Gpbuj30rxfs0gd4\nGqgFfhbnMH0i1tcAq4DHsFTOBCyv3wWcxPL9W4C7gSfi1UuBX0TEC29RdBbwTWAnsN2VPYjdHfAP\nwHDg39y2m7APiOfdYwewGAv6uPXlQDGwFng1M1UWEcl33sbqeZv4vSvjpX0edku0bcCUZN7Uz4G/\nHvvq0oldsZ4BDAP+Bbjcbb8dOJGd6omI4O8oGoef+/F3ASGse9MMV/YAsA6YCLzunouIZI+3G7iy\nws+BH86/qAFwM7DCra8gwQ0KIiK9wtsNXFnho8+gC3QB67FUzz8BPwfKgGa3vdk9FxHJHs25m1az\ngMPACCy9E90btovw1WwRkezw0Zl8svwc+A+7x0+AF7E8fzMwChvfYjTQEuuFS1cBb8LOgbW0hibB\ntHm9UF0R8buajVCzyYLK++k6qJ+jaBx+rfIA7HP0FDAQG3r0IezmhYXAI+4xZnenpXcB18FzI6p4\nkSrrGyQieS80C0IzYc871hl+QzoO6tcomoBfq1yGneWD1fE54DVgK/b/tYhwd04RkezxaxRNwK9V\nPghcHaP8GDC7l+siIhKfcvwiInkmgFE0gFUWEfGRAM65q8AvIuJFAKNoAKssIuIjAYyiAayyiIiP\nBDCKBrDKIiL+0RXAXj1+H6RNRMTXOguTX2J4BhuRYFdE2QxsFq3twO+AayO2LQHqsCFs5kSUT3fH\nqAMev1idFfhFRDzwGPifBeZGlT0K/BU2JP1fu+cAVcAC9zgXeJLwCMbLsBtbJ7gl+pjnUapHRMSD\n9qL+Kex9JrrgLaAyquwwcKlbH0J40vRbgNXYxFT1wH5sjt2PgBLsWwLASmzIek29KCKSCZ0FaU/y\nP4BNyfj3WFbmP7jycuDdiP0asAnXz7r1bo2uPC4FfhERDzoTjNmwsaaDjTWdqR7yaeBebLyy27Dr\nADf2tH6xKPCLiHjQkSDwV4cKqA6Fn//dQ58mc8gZhMckewF4yq03AhUR+43FzvQb3XpkeSMJ6OKu\niIgHnRQmvSRpP3CDW/8qsM+trwHuAPoD47CLuFuw+UlOYvn+PsDdxBmyvpvO+EVEPEiU6knCaizI\nDwcOYb14/gL4P9goQG3uOUAtNix9LdABLCY8C+FiYDlQDKwlwYVdUOAXEfHEY+C/M055dZzyh90S\nbRswJdk3VeAXEfGgnVS6c/qDAr+IiAcp5O59I3g1FhHxEY+pnqxQ4BcR8UCBX0QkzyTqx+9XCvwi\nIh4oxy8ikmeU6hERyTNn1J1TRCS/KMcvIpJnlOMXEckzyvGLiOQZBX4RkTyjHL+ISJ45Q1G2q5Ay\nBX4REQ+CmOrRDFwiIh50UJD0EsMzQDOwK6JsKTal4na33BSxbQlQB+wF5kSUT3fHqAMev1idFfhF\nRDzwOPXis8DcqLIu4DHgGre84sqrgAXucS7wJDbVIsAyYBE2HeOEGMc8jwK/iIgHnRQkvcTwFnA8\nRnmfGGW3YFM1ngXqsbl5q4HRQAk2/y7ASuDWRHVW4BcR8cBj4I/nHuA94GlgiCsrx1JA3RqAMTHK\nG115XLq4KyLiQaKAvq/mMHU1h1M95DLgR279x8BPsTRO2ijwi4h40J6gO+floUouD1X+/vkrD21P\n5pAtEetPAS+79UagImLbWOxMv9GtR5Y3JnqDIKZ65mJXtOuA72e5LiKS5zKQ6hkdsT6fcI+fNcAd\nQH9gHHYRdwtwBDiJ5fv7AHcDLyV6g6Cd8RcA/wjMxj7Rfof9MPZks1Iikr889uNfDdwADAcOAT8E\nQsDVWO+eg8Bfun1rgefdYwew2O2DW18OFANrgVcTvWkygf9e4BfEvvLc22ZgV7Lr3fNfYVe6FfhF\nJCs8DtlwZ4yyZxLs/7Bbom0DpiT7psmkesqwM+vnsTRLrG5GvWUM9qnYrfuqtohIVnjsx58VydTk\nB8BfYXeJfRtLtTyPdTM6kLGaxdZ18V1g6SrgTdg5sJbW0CSYNi/D1RKRIKjZCDWb4BPg/TQdM4hD\nNiT7EXQOu4DQDHQCQ4EXgPXA/8hM1WKKvqpdwfn9VwFYehdwHTw3oooXqbLLHiKS90KzIDQT9rxj\nZ68b0nDMXA389wHfAlqxrkX/HbtzrC/Ws6Y3A/9W7Ep2JdCE3b4cK0cmItIr2nN0zt1hwDeAj6LK\nzwFfT3uNEusA/gvwW6yHz9Powq6IZJGfcvfJSqbGP0ywrTZdFUnBK4QHLRIRyapcTfWIiEgcCvwi\nInlGUy+KiOSZXM3xi4hIHEr1iIjkmTM52p1TRETiUI5fRCTPKMcvIpJngpjjD+JELCIivuFxIpZn\nsDHQdkWU/R02IsF7wG+ASyO2LcGGytmLDZzZbbo7Rh3w+MXqrMAvIuJBBwVJLzE8iw13H+k14A+A\nPwT2YcEeoAobn6zKveZJwsPkL8Pm5Z3gluhjnkeBX0TEA4/j8b/FhZNcrcPGQgPYTHg+3VuwGbvO\nYpNR7cemWxwNlGDTMAKsBG5NVGfl+EVEPMhwd84/w4I9QDnwbsS27omoznL+8PSNXGSCKgV+EREP\nMtid8wfAGWBVug+swC8i4kGi7pyf1uzgZM2Onhz228A84GsRZdETUY3FzvQbCaeDussbEx1cgV9E\nxINE3TkHhaYzKDT9988bHlqZzCHnYhNc3QB8EVG+Bjv7fwxL5UzA8vpd2DyD1e753cATid5AgV9E\nxAOP/fhXYwF+OHAIm/9kCdAfu8gL8A6wGJv/5Hn32OHKuuchXwwsB4qBtcCrid5UgV9ExAOPgT/W\n1LHPJNj/YbdE2wZMSfZNFfhFRDxopyjbVUiZAr+IiAdBHLJBgV9ExAMFfhGRPKNhmUVE8oyGZRYR\nyTNK9YiI5BkFfhGRPNN+RnPuiojklc6O4IXR4NVYRMRHOjuU6hERySsK/CIieabjrAK/iEheOdcZ\nvDAavBqLiPiJUj0iInnmi+CF0b7ZroCISKB1pLDEdh+wC9jt1gGGYROx7ANeA4ZE7L8EqAP2AnN6\nUmUFfhERL7wF/snAd4FrgT8E/hNwJfAAFvgnAq+75wBVwAL3OBd4kh7EcQV+EREvvAX+LwGbsbl1\nO4ENwJ8ANwMr3D4rgFvd+i3YdI1ngXpgPzAj1Sr7MfAvxWaO3+6WmyK2ef6KIyKSVmdTWC60G/gK\nltoZAMwDxgJlQLPbp9k9ByjH4mO3Bmzi9ZT48apEFzaL/GNR5ZFfccYA67GvQed6tXYiIpE6E2z7\n9xrYXpPo1XuBR7A8/mfAjhhH7CI8qXosibbF5MfAD9AnRlm8rzjv9l61RESixL9oC1NDtnR75qFY\nez1DeIL1v8XO4puBUcARYDTQ4rY3AhURrx3rylLix1QPwD3Ae8DThK9mp+UrjohIWn2RwhLbSPd4\nGfANYBWwBljoyhcCL7n1NcAdQH9gHDAB2JJqlbN1xr8O+zSL9gNgGfAj9/zHwE+BRXGOk/JXHBGR\ntEp0xp+cF4BSLJuxGPgU+AnwPBb76oHb3b61rrzWvfNiApTquTHJ/Z4CXnbrSX/FWboKeBN2Dqyl\nNTQJps3rcUVFJHfUbISaTfAJ8H66Duo98F8fo+wYMDvO/g+7pcf8mOMfDRx26/OxGxvAvuKswi76\njiHBV5yldwHXwXMjqniRKjiZ2QqLSDCEZkFoJux5x06bN6TjoN4Df6/zY+B/BLga+/pyEPhLV56W\nrzgiImkVu5umr/kx8H8rwTbPX3FERNIqUXdOn/Jj4BcRCQ6lekRE8kz8bpq+pcAvIuKFzvhFRPKM\nAr+ISJ5R4BcRyTPqzikikmfUnVNEJM+oV4+ISJ5Rjl9EJM8oxy8ikmcCmOP360QsIiLB4G2ydbDJ\npl4A9mCDUFZjc/CuA/Zh0zIOidjf89zjCvwiIl54D/yPA2uBq4CpWEB/AAv8E4HX3XM4f+7xucCT\n9CCOK/CLiHhxNoXlQpcCXyE8524HNgPXzcAKV7YCuNWtx5t7PCUK/CIiXrSnsFxoHDYh2LPAvwM/\nBwYCZdiE67jHMreelrnHFfhFRLzwluopBKZhKZtpwGeE0zrdukg86VRg5twVEckNibpzttTAJzWJ\nXt3glt+55y9gF2+PAKPc42igxW1Peu7xRBT4RUS8SNSdszRkS7c9D0XvcQQ4hF3E3YdNsP6+WxZi\nU9EuBF5y+yc993giCvwiIl54v3P3HuA5oD9wAPgOUIDNMb4Iu4h7u9s3LXOPK/CLiHjhPfC/B1wb\no3x2nP09zz2uwC8i4oWGbBARyTOxu2n6mgK/iIgXGp1TRCTPKNUjIpJnAjg6pwK/iIgXSvWIiOQZ\nBX4RkTyjHL+ISJ5Rd04RkTyjVI+ISJ5RqkdEJM+oO6eISJ5RqkdEJM8o8IuI5JkA5vg1566IiBfe\n5ty9BNgM7MAmV/lfrnwYsA6bles1YEjEa5YAdcBeYE5PqpytwH8bNrVYJzbBcKR4jZoO7HLbHu+F\nOoqIZNoXwB8BVwNT3fp/xCZcX4dNyfg64QnYq4AF7nEuNkl7ynE8W4F/FzAfeDOqPFaj+rhty7Bp\nyCa4ZW6v1FREJLM+d4/9sSkXjwM3Aytc+QrgVrd+C7AaSzDVA/uBGam+YbYC/17sK0y0WI2qxmaZ\nLyE8qfBKwj8IEZEg64ulepqBN7BsSJl7jnssc+vlQEPEaxuwSddT4reLu+XAuxHPuxt1lvMb20gP\nGisikn6Jru5ucEtC57BUz6XAb7F0T6QuEk+o7qvJ1tcBo2KUPwi8nMH3Zekq4E3YObCW1tAkmDYv\nk28nIgFRsxFqNsEn2Gl1eiTqzznLLd3+JtGBPgX+Dbue2YzFzyNYxqPF7dMIVES8ZqwrS0kmA/+N\nPXhNrEY1uPKxUeVxG7v0LuA6eG5EFS9SBSd7UBMRyTmhWRCaCXvegedJ4lw8KZ76cw7HPjlOAMVY\n3HwIWAMsBB5xjy+5/dcAq4DHsKzHBMIp8KT5IdXTJ2I9XqO6sPBd7Z7fDTzRu9UUEYmlzcuLR2MX\nb/u65RdYL57t2GfTIux65+1u/1pXXot9YCzGZ6meROZjgXs49tVmO3ATiRu1GFiOfSquBV7t1RqL\niMTk6Yx/Fxd2aQc4BsyO85qH3dJj2Qr8L7ollniN2gZMyViNRER6JHhjNvgh1SMiEmDBG7NBgV9E\nxBOd8YuI5Bmd8YuI5BlPvXqyQoFfRMQTpXr8wf0/FETPiRaztbn5IxDJL/2iHp2C2HunNzmjVI8/\ntAOnof+Idopop/8l7bQVAUXY6NclwPEBbqXYLSISPIWE/4ZLrKiE8N96MQygjWI+ZwCfw2ngs3Qn\nZ3TG7w+twDEYOu4EQzjBkP4n+HTEKBiB3TL2BXB8GDbXwSn3ojaC+B8okt8KgVLsb7nU4v8Qwn/r\nw9sp4RRDsVjAMaC1B4PbJKQzfn/YAQyEidP3MZEP2M94Dl9dzpkjg21EjO65bA5OwH5pWgl/AIhI\ncBRiQb8MSvrBJGwQhMm2XFm+n0l8wEQ+YPje03af7A57SJ/gnTDmZuDfDVwKYw61Mr7iAOPZT/2w\nSvZNngpHscDfgeX/GodB2zD3wuB9covkt3422lcpUInd2z8KmAzFk49zJQe4kgOM54DFhd1QX5fu\nOgQvbuRk4D+2FYYVAVthSsUu6qmkmTJOTS3h8OlxcAhr+RDs7OAolv6JvjAkIv5WgOX0R2GB/0vY\nlCVfbqdqcC1T2MkUdjLu4GHYCmyFmrRXQt05feE3HfD1DVB2NVw2oYXpk7dxnCF8TjHMhObmkZwb\nOtB+WRqwwbnbCeI3NpH8dgkwiHDgnwx9yz5jatkuprOVL7ONL7MNNgIboPld+CjtldAZvy80Av8X\n+O7jMOYzuO7OHZR/tYnxHKCKWlrKyvhgziTq51TS0lQGR4vsQ7vzIgcWEX9xgb941HHKBzdxJQco\no5kp7KKazVx/eAv8BngKlu/IRNCHIJ4x5mTg7/YUcP9qKC6EyzpaKJnzW0o5SgtlVHCIQ1TQVF7O\nqfIS+zYgIoFSSCeDOMVwWinn/MA/9eA+G8B9NazPWNCHNJzxzwV+hiWunsImX8monA78AI9+Bt9b\nDoM/haEtbdw4621OXtafiQUfUM84WhjJcYbQxoCUj/1RTT2XhyrTXudsysU2QW62KxfbBKm1q4BO\nSjhFKa1U8DHjOUBpZyuDN56x6UzWWNDfmNEaezrjLwD+ERt7vxH4HTYh1R7v9Yov5wM/wP/+Auav\nhqkHgd0wuOIM107ZzbUTdnN09CBOMITPexD4l9UcZUFoePornEW52CbIzXblYpsgtXb9PvC3H2Xg\nwXPWc6cZi/Svw/KWTJ7pd/N0xj8D2I/NsgXwK+AWFPjT40Vgy7vwja0wrByb0/4qGF5xmuGlp+1O\nvxR/GmUfwNSXj6W/slmUi22C3GxXLrYJUmxXOzZFeQvwIRb4m2Djx7A+YzWM5umMfwzWz7BbAzbF\nbEblTeAH+x71Dx3Ax7YUr4Hx2E++mNQ7c34C7PlVeuuYbbnYJsjNduVimyC1drVhN+M20Btn9olq\n0WMpz5ebDn0uvkvg1AA3ZLsSIhIIG4CQh9enGrhPAYMjnl8HLMUu8AIsAc7RCxd4RUQkOwqBA9hd\nCP2xAWeuymaFREQk824CPsAu8i7Jcl1ERETy123A+9i9vdOiti0B6oC9wJyI8unYIIB1wOO9UEev\nlmLXyLa75aaIbfHaGARzsXrXAd/Pcl28qgd2Yv8/W1zZMGAdsA94jfDYs371DNbhMnKAzERtCPLv\nngTcl4CJwBucH/irsJxcPyxHt5/wBfMtWB9dsPsH5+JvPwT+a4zyWG3s23vV8qQAq28lVv+g508P\nYkEy0qPA/W79+8BPerVGqfsKcA3nB/54bQjy756v6YeYnL3Y2Ui0W4DV2B0c9dgvZjU25mcJ4bOy\nlcCtGa+ld7F6ecVq44wY+/lR5M0xZwnfHBNk0f9HNwMr3PoK/P979hZwPKosXhuC/Lvnawr83pRj\n6ZFuDdhtAdHlja7c7+4B3gOeJvx1O14bgyDWzTFBqXssXdh9SVuBP3dlZVjqBPdYloV6eRWvDUH+\n3fO1vLqB6yLWYYO7RnsQeLmX65Ip8dr4A2AZ8CP3/MfAT4FFcY6TlZtOeiAo9UzWLOAwNrHgOuyb\naKQugt/mi7Uh6O3zBQX+sBt78JpGoCLi+VjsrKTRrUeWp3eaz55Jto1PEf6wi9VGP7QlGdF1r+D8\nM8igOeweP8FGIZmBnSGPAo5gKcaW7FTNk3htCPLvnq8p1ZO6yBzrGuAO7MaLccAELK9/BDiJ5fv7\nAHcDL/VuNVM2OmJ9PuGLb/HaGARbsfpWYvVfgLUniAZg140ABmI9XHZh7Vnoyhfi/9+zWOK1Ici/\ne5ID5mO54jYsqL8Sse1B7KLTXuCPI8q7u3PuB57onWp6shLrKvge9ocXmSuO18YgyJWbY8ZhPVx2\nYEORdbdwseShAAAAt0lEQVRlGJb3D0p3ztVAE3AG+5v6DonbEOTfPRERERERERERERERERERERER\nEREREREREREREZHccS12J3IRNszBbmyMd5G8Fmv8dZFc8mPgEqAYGyLgkexWR0REMq0fdtb/LjrR\nEQE0OqfkvuFYmmcQdtYvkvd0BiS5bg2wCrgCG3r6nuxWR0REMulbwP9z632xdE8oa7URERERERER\nERERERERERERERERERERERERERGR3Pb/AarWjtK2Y6UFAAAAAElFTkSuQmCC\n", + "text": [ + "" + ] + } + ], + "prompt_number": 24 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "plt.plot(mesh.vectorCCx, sig[mesh.gridCC[:,1]==mesh.vectorCCy[125]], '.')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 25, + "text": [ + "[]" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE/5JREFUeJzt3X+snXV9wPF3pSBUftwQSLktaDt+BEpMLFUw0cWzRLvi\nElqy8DMjdRKvphsalwjUJetVkw01mtUssBBRyrZ2aTAizIIUwnHERDscPyq1649ZYu/oRR1YSYyW\neffH93s5D8d7es/v83y/z/uVnJznfs9zznm+uef5fM7z/XVAkiRJkiRJkiRJkiRJkiRJFXMy8APg\nGWAP8Hex/ExgJ7APeBQYKzxnI7Af2AusLpSvAnbHxzYP9KglSX2zKN4vBL4PvBf4AnBrLL8NuCNu\nryAkjBOBZcABYEF8bBdwedzeAawZ5EFLkvprEfAfwKWEb/mLY/k58W8IVwG3FZ7zCPBuYBz4caH8\neuAfB3mwkqT5vanNfZ4BpoEngOcJCWA6Pj5NIyEsAQ4XnnsYWDpH+VQslySN0MI29vkd8A7gDOA7\nwB81PT4Tb5KkxLSTBGb9Evg2oYN3mtAMdITQ1PNS3GcKOK/wnHMJVwBTcbtYPtX8Bueff/7MwYMH\nOzgkSRJwELigmyfO1xx0Fo2RP6cAHwCeBh4E1sfy9cADcftBQnv/ScBy4EJCh/AR4ChwBaGj+KbC\ncxq1OHiQmZmZbG+bNm0a+TFYP+tXtbpVoX7A+d0kAJj/SmAc2EJIFm8C/gl4nJAItgM3A4eAa+P+\ne2L5HuA1YAONpqINwL2EZLKD0GksSRqh+ZLAbuCyOcr/F3h/i+f8bbw1+yHw9vYPTZI0aO2MDlKf\n1Gq1UR/CQFm/dOVcN8i/fr1YMP8uQzUT27ckSW1asGABdBnPvRKQpAozCUhShZkEJKnCTAKSVGGV\nTQITE1CrwQc/CK+8MuqjkZSL1GJLZZPAvn3w3e/Cww+Hf5ok9UNqsaWySWBR/JWEd74T7r57tMci\nKR+pxZbKzhN45ZWQpe++G8bG5t9fktoxitjSyzyByiYBScqFk8UkSV0xCUhShZkEJKnCTAKSVGEm\nAUmqMJOAJFWYSUCSKswkIEkVZhKQpAozCUhShZkEJKnCKp0EUlv3W1K5pRhTKp0EUlv3W1K5pRhT\nKp0EUlv3W1K5pRhTKr2UtL8pIKmfRhVT/D0BSaowf09AktSV+ZLAecATwPPAj4CPx/JJ4DDwdLxd\nWXjORmA/sBdYXShfBeyOj23u8bglSX0w3+XDOfH2DHAq8ENgHXAt8Cvgy037rwC2Au8ClgKPARcC\nM8Au4C/j/Q7gK8AjTc+3OUiSOjTI5qAjhAQA8CrwY0Jwb/WGa4FtwDHgEHAAuAIYB04jJACA+wjJ\nRJI0Qp30CSwDVgLfj3/fAjwL3APM9oMvITQTzTpMSBrN5VM0kokkaUTaTQKnAvcDnyBcEdwFLAfe\nAbwIfGkgRydJGqiFbexzIvAN4J+BB2LZS4XHvwo8FLenCJ3Js84lXAFMxe1i+dRcbzY5Ofn6dq1W\no1artXGIklQd9Xqder3el9earyNhAbAF+AXwyUL5OOEKgFj+LuBGGh3Dl9PoGL6A0DH8A8Lool3A\nt7FjWJL6opeO4fmuBN4D/BnwHGEoKMCngRsITUEzwE+Aj8bH9gDb4/1rwIa4D3H7XuAUwuig5gQg\nSRoyZwxLUuKcMdyDFJd+lVROKcaTyieBFJd+lVROKcaTyieBFJd+lVROKcaTyvcJuJy0pH5xKene\n2TEsSR2yY1iS1BWTgCRVmElAkirMJCBJFWYSkKQKMwlIUoWZBCSpwkwCklRhJgFJqrDKJ4EUV/2T\nVD6pxpLKJ4EUV/2TVD6pxpLKJ4EUV/2TVD6pxpLKLyDnKqKS+mGUscRVRCWpwlxFVJLUFZOAJFWY\nSUCSKswkIEkVZhKQpAozCUhShZkEJKnCTAKSVGEmAUmqMJMA6a7+J6k8Uo0j8yWB84AngOeBHwEf\nj+VnAjuBfcCjQHGljI3AfmAvsLpQvgrYHR/b3OuB91Oqq/9JKo9U48h8SeAY8EngUuDdwF8AlwC3\nE5LARcDj8W+AFcB18X4NcCeN9SzuAm4GLoy3Nf2qRK9SXf1PUnmkGkfmSwJHgGfi9qvAj4GlwFXA\nlli+BVgXt9cC2wjJ4xBwALgCGAdOA3bF/e4rPGfktm6Fa66BnTtdSVRSd1KNIws72HcZsBL4AbAY\nmI7l0/FvgCXA9wvPOUxIGsfi9qypWF4KY2Owffuoj0JSylKNI+0mgVOBbwCfAH7V9NhMvPXF5OTk\n69u1Wo1ardavl5akLNTrder1el9eq531p08E/g14GPj7WLYXqBGai8YJnccX0+gbuCPePwJsAl6I\n+1wSy28A3gd8rOm9/D0BSerQIH9PYAFwD7CHRgIAeBBYH7fXAw8Uyq8HTgKWEzqAdxGSxVFC/8AC\n4KbCcyRJIzJf5ngv8O/AczSafDYSAvt24K2EDuBrgdmRsZ8GPgy8Rmg++k4sXwXcC5wC7KAx3LTI\nKwFJ6pA/LylJFebPS0qSumISkKQKMwlIUoWZBEh34SdJ5ZByDDEJkO7CT5LKIeUYYhIg3YWfJJVD\nyjHEIaKEy7eJifDPS2nhJ0nlMOoY4jwBSaow5wlIkrpiEpCkCjMJSFKFmQQkqcJMApJUYSYBSaow\nk4AkVZhJIEp57Q9Jo5N67DAJRCmv/SFpdFKPHSaBKOW1PySNTuqxw2UjolGv/SEpTWWIHa4dJEkV\n5tpBkqSumAQkqcJMApJUYSYBSaowk4AkVZhJQJIqzCQgSRVmEohSX/9D0mikHjvaSQJfA6aB3YWy\nSeAw8HS8XVl4bCOwH9gLrC6Ur4qvsR/Y3PURD0jq639IGo3UY0c7SeDrwJqmshngy8DKeHs4lq8A\nrov3a4A7acxiuwu4Gbgw3ppfc6RSX/9D0mikHjvaSQJPAi/PUT7XFOW1wDbgGHAIOABcAYwDpwG7\n4n73Aes6PNaB2roVrrkGdu507SBJ7Us9dvTSJ3AL8CxwDzBb9SWEZqJZh4Glc5RPxfLSGBuD7dvT\n/CdKGp3UY8fCLp93F/DZuP054EuEpp6eTU5Ovr5dq9Wo1Wr9eFlJyka9Xqder/fltdpddW4Z8BDw\n9nkeuz2W3RHvHwE2AS8ATwCXxPIbgPcBH2t6LVcRlaQOjWIV0fHC9tU0Rg49CFwPnAQsJ3QA7wKO\nAEcJ/QMLgJuAB7p8b0lSn7TTHLSN8K39LOCnhG/2NeAdhFFCPwE+GvfdA2yP968BG+I+xO17gVOA\nHYSrBEnSCPmjMpKUOH9URpLUFZNAQerTvyUNVw4xwyRQkPr0b0nDlUPMMAkUpD79W9Jw5RAz7Bgu\neOWVkM3vvjvd2X+ShqcsMaOXjmGTgCQlztFBkqSumAQkqcJMApJUYSYBSaowk4AkVZhJoCCH2X+S\nhieHmGESKMhh9p+k4ckhZpgECnKY/SdpeHKIGU4WKyjL7D9JaShLzHDGsCRVmDOGJUldMQlIUoWZ\nBCSpwkwCklRhJoEmOUz+kDR4ucQKk0CTHCZ/SBq8XGKFSaBJDpM/JA1eLrHCeQJNyjL5Q1K5lSlW\nOFlMkirMyWKSpK6YBCSpwkwCklRh7SSBrwHTwO5C2ZnATmAf8ChQ7BbZCOwH9gKrC+Wr4mvsBzZ3\nf8iSpH5pJwl8HVjTVHY7IQlcBDwe/wZYAVwX79cAd9LorLgLuBm4MN6aX1OSNGTtJIEngZebyq4C\ntsTtLcC6uL0W2AYcAw4BB4ArgHHgNGBX3O++wnNKJ5eZgJIGJ5c40W2fwGJCExHxfnHcXgIcLux3\nGFg6R/lULC+lXGYCShqcXOLEwj68xky89cXk5OTr27VajVqt1q+XblsuMwElDc4o40S9Xqder/fl\ntdqdXLAMeAh4e/x7L1ADjhCaep4ALqbRN3BHvH8E2AS8EPe5JJbfALwP+FjT+5RisliZZgJKKqcy\nxYlRTBZ7EFgft9cDDxTKrwdOApYTOoB3EZLFUUL/wALgpsJzSmdsDLZvH/0/VlJ55RIn2mkO2kb4\n1n4W8FPgbwjf9LcTRvscAq6N++6J5XuA14ANNJqKNgD3AqcAOwhXCZKkEXLtIElKnGsHSZK6YhKY\nQy7jfyUNRk4xwiQwh1zG/0oajJxihElgDs4TkHQ8OcUIO4bnUKbxv5LKp2wxwl8Wk6QKc3SQJKkr\nJgFJqjCTgCRVmEmghZzGAUvqr5zig0mghZzGAUvqr5zig0mghZzGAUvqr5zig0NEWyjbOGBJ5VG2\n+OA8AUmqMOcJSJK6YhKQpAozCbSQ0xAwSf2TW2wwCbSQ0xAwSf2TW2wwCbSQ0xAwSf2TW2xwdFAL\nZRsCJqkcyhgbHCIqSRXmEFFJUldMApJUYSaB48htKJik3uQYE0wCx5HbUDBJvckxJpgEjiO3oWCS\nepNjTHB00HGUcSiYpNEpa0xwiKgkVdgoh4geAp4DngZ2xbIzgZ3APuBRoJgvNwL7gb3A6h7fW5LU\no16TwAxQA1YCl8ey2wlJ4CLg8fg3wArguni/BrizD+8vSepBP4Jw8yXIVcCWuL0FWBe31wLbgGOE\nK4gDNBJHKeU4HExS93KMCf24EngMeAr4SCxbDEzH7en4N8AS4HDhuYeBpT2+/0DlOBxMUvdyjAkL\ne3z+e4AXgbMJTUB7mx6fibdWSt0LnONwMEndyzEm9JoEXoz3PwO+SWjemQbOAY4A48BLcZ8p4LzC\nc8+NZW8wOTn5+natVqNWq/V4iN3burWcw8EkjUZZYkK9Xqder/fltXoZIroIOAH4FfAWwkigzwDv\nB34BfJ7QKTwW71cAWwmJYimhGekC3ng14BBRSepQL0NEe7kSWEz49j/7Ov9CSARPAduBmwkdwNfG\nffbE8j3Aa8AGSt4cJEm5c7LYPCYmQmfQokXhUtBmIamayhwL/D2BAcpxNICkzuUaC0wC88hxNICk\nzuUaC2wOmkdZF4ySNFxljgUuICdJFWafwADlOE1cUudyjQUmgXnk2hkkqTO5xgKTwDxy7QyS1Jlc\nY4F9AvMoc2eQpOEpcyywY1iSKsyO4QHLtUNIUntyjgEmgTbk2iEkqT05xwCTQBty7RCS1J6cY4B9\nAm0oc4eQpMErewywT2DAxsbCbd26PNsEJR3frbfCSy/BjTfmd/6bBNqUc5ugpOPL+fw3CbQp5zZB\nSceX8/lvn0Cbyt4mKGlwyn7+2ycwBDm3CUpqbWIi9Ae++uqoj2QwTAJtyrlNUFJruZ/7JoE25dwm\nKKm13M99+wTaVPY2QUmDkcK5b5/AEDhXQKqm3PsDTQIdyL1tUNLvy/28Nwl0IPe2QUm/L/fz3iTQ\ngbPPhrPOKm+7oKT+mpiAo0fhnHPg/vvzPPdNAh144QX4+c/hscfyvCyU9Eb79sH3vgdHjsCnPjXq\noxkMk0AHZi8LTz0VXn45z04iScHEBDz3XNheuTLPpiAwCXRk69bQHPTqq14NSLnbty982QN461vz\nbAqC4SeBNcBeYD9w25Dfu2djY3DyyWH79NPhi18c7fFIGpyDB8P9GWfA5s2jPZZBGmYSOAH4B0Ii\nWAHcAFwyxPfvi7e9LdwfPdp5G2G9Xu/78ZSJ9UtXznWD7uo3e67/8pf59gfAcJPA5cAB4BBwDPhX\nYO0Q378vTj893HfTL+CJlrac65dz3aDz+k1MwJ49YTvn/gAYbhJYCvy08PfhWJaUrVvhzW9u9At8\n6EOjPiJJ/fbQQ43+gPHxfPsDABYO8b3KuShQh8bGwiih3/wm/P2tb8GCDlbs+MxnBnNcZWH90pVz\n3aD7+p14Yn+Po2yGuYDcu4FJQp8AwEbgd8DnC/scAM4f4jFJUg4OAheM+iDms5BwoMuAk4BnSLBj\nWJLUvSuB/yJ849844mORJEmSNGzXAM8D/wdcVihfBvwaeDre7iw8tgrYTZhoVvapG63qB+EKaD9h\n0tzqQnlK9SuaJIz0mv2fXVl4rFVdU5P0JMcWDgHPEf5nu2LZmcBOYB/wKJDSmJivAdOEc2jW8eqT\n0mdzrrpNkvh5dzFwEfAEv58Eds/1BMIH9fK4vYNGB3MZtarfCkJfyImEuh6g0TmfUv2KNgF/NUf5\nXHVNcZmSEwjHvoxQl1z6sn5CCJJFXwBujdu3AXcM9Yh684fASt4YP1rVJ7XP5lx169t5N6qK7yVk\n53aNA6fR+MZyH7Cu3wfVR63qtxbYRpgsd4jwD7qC9OrXbK5RZnPV9fI59iu7LCY5ttD8f7sK2BK3\nt5DWZ/BJ4OWmslb1Se2zOVfdoE/nXRmz33LC5U0deG8sW0q49Jk1RYITzYAlvLEesxPmmstTq98t\nwLPAPTQuuVvVNTVZTHKcwwzwGPAU8JFYtpjQ7EC8XzyC4+qnVvXJ5bPZl/NukJPFdgLnzFH+aeCh\nFs/5H+A8Qta7DHgAuHQgR9e7buqXqlZ1/WvgLuCz8e/PAV8Cbm7xOilOGEzxmNvxHuBF4GzC/3dv\n0+Mz5FX3+eqTWl37dt4NMgl8oIvn/DbeAP6TMK/gQsI343ML+50by0apm/pNEZLcrHMJmbqM9Stq\nt65fpZEA56prmerUruZ6nMcbv2ml6sV4/zPgm4Qmg2lCsj9CaKJ8aTSH1jet6pPDZ7P4v+npvCtD\nc1CxXessQkccwB8QEsB/Ez6wRwnt5wuAmwhXCSko1u9B4HrCZLnlhPrtInxIU63feGH7ahqdV63q\nmpqnCMe+jFCX6wh1S9kiQh8UwFsII0h2E+q1PpavJ53PYCut6pPDZzP58+5qQjvrrwkB8OFY/qfA\njwh9Aj8E/qTwnNkhlAeArwztSLvTqn4QmosOEC6//7hQnlL9iu4jDDV8lnCSFduRW9U1NblNclxO\nGEHyDOF8m63TmYR+ghSHiG4jNCf/lnDu/TnHr09Kn83mun2Yapx3kiRJkiRJkiRJkiRJkiRJkiRJ\nkqR2/T9WIn4U8vLtbQAAAABJRU5ErkJggg==\n", + "text": [ + "" + ] + } + ], + "prompt_number": 25 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "sigy = zeros_like(mesh.gridCC[:,1])\n", + "sigy[mesh.gridCC[:,1]<-ay] = (mesh.gridCC[:,1][mesh.gridCC[:,1]<-ay]-ay)**2\n", + "refy = sigy[mesh.gridCC[:,0]<-ay].max()\n", + "sigy[mesh.gridCC[:,1]>ay] = (mesh.gridCC[:,1][mesh.gridCC[:,1]>ay]+ay)**2\n", + "sigy = sigy/refy" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 26 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def ricker(fpeak, t, tlag):\n", + " \"\"\"\n", + " Generating Ricker Wavelet\n", + " \n", + " .. math ::\n", + " \n", + " \n", + " \"\"\"\n", + " return (1-2*np.pi**2*fpeak**2*(t-tlag)**2)*np.exp(-np.pi**2*fpeak**2*(t-tlag)**2)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 29 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "wave = ricker(400, time, 0.0025)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 30 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fig, ax = plt.subplots(1,1, figsize = (7, 5))\n", + "ax.plot(time, wave, '.-')\n", + "ax.set_xlim(0, 0.02)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 31, + "text": [ + "(0, 0.02)" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAE4CAYAAADPUy0vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUHHWd9/H3MJOQK5mES+5s5CYEuQ8Xb+uAcDZGBFbI\n4UkexcsjwccHl12PCj7n+DAehUV3VfTg4xIUNq4GCLCrYTGE6KHVBzCiQEgCgSSahYQQCLmQC7n3\n88evh5lMuqeruqqnq3ver3P6dHV3dVUlRefD93epAkmSJEmSJEmSJEmSJEmSJEmSpEy6E1gPLOll\nne8DK4DFwBl9cVCSJFXD+wlBVir0pgK/LCyfC/y+Lw5KkqRqmUTp0PsX4Mpur5cDo6t9QJIk9XRI\nH+xjPPByt9drgAl9sF9Jkg7QF6EH0NTjdb6P9itJ0tta+mAfa4GJ3V5PKLx3gGOPPTa/atWqPjgc\nSVKdWAUcl+YG+6LSmwdcVVg+D9hMGO15gFWrVpHP5+vyMXt2npaWPJCnrS3Ppk21P6ZKHjfeeGPN\nj8GH5yErD89D7R/AsWkHUhqhdzfwOPBOQt/dp4FrCg8IIzf/DKwEbgc+l8I+M+WBB+Dmm6G5GX7x\nC2htrfURSZKKSaN5c3qEda5NYT+ZtHUr5HIwezY89BA8/TSMG1fro5IkFdNXA1ka1hVXwCGHwIwZ\nIQCvvhqmToXNm2t9ZPG1t7fX+hCE5yErPA+NqeeoylrKF9pw68pxx0Hn+JsRI2DLlrA8bRrMnVu7\n45KketfU1AQp55SVXkL794fntjZ417u6lmfNqt0xSZKK64spCw2trQ2OOAIeeQR27YKxY2HBAgez\nSFIWGXoJvfYa3HJLV8gdfjjs2VPbY5IkFWfzZkJr18L48V2vx42DV16p3fFIkkoz9BLI5w8OvfHj\nw3uSpOwx9BLYtAkGDoRhw7res9KTpOwy9BLoWeWBoSdJWWboJVAs9GzelKTsMvQSWLMGJvS4M6CV\nniRll6GXgM2bklRfDL0EilV6Nm9KUnYZegkUq/SOPDJcbHr37tockySpNC84ncDhh8OkSTB6NMyZ\n03VVlokT4bHH4Oija3p4klTXqnHBaUMvgYEDuy451v2uCkcdFSrAsWMPDENJUnSGXsY0N4e7LLS1\nwcKFB15/c+PGsOwthiSpMt5aKEN27gyhN23agYEHMGhQePYWQ5KULd5loUKbNsGoUcWruIsvhkWL\nDg5DSVJtWelVaNMmGDmy+Gdjx8Kllxp4kpQ1hl6FNm4MlV4xI0eGUJQkZYuhV6HeKj1DT5KyydCr\nkJWeJNUfQ69CVnqSVH8MvQoZepJUfwy9Ctm8KUn1J43QmwIsB1YA1xf5/AjgYeAZYCnwyRT2WXNW\nepJUf5KGXjNwGyH4JgPTgZN6rHMt8DRwOtAOfJsGmBTfOTm9mMGDw/Nbb/Xd8UiSyksaeucAK4HV\nwB7gHuDSHuusAw4rLB8GvAHsTbjfmtu4sXSlB1Z7kpRFSUNvPPByt9drCu91dwdwMvAKsBi4LuE+\nM6G35k0w9CQpi5KGXpTbIvxvQn/eOEIT5w+A4Qn3W3O9NW+CoSdJWZS0b20tMLHb64mEaq+79wA3\nFZZXAX8B3gn8sefGOjo63l5ub2+nvb094eFVRz5vpSdJacvlcuRyuaruI+l9ilqAF4APEpov/0AY\nzPJ8t3W+A2wBvgaMBv4EnAps7LGturmf3rZt4W7p27eXXufjH4eLLoKrruq745KkRlKN++klrfT2\nEkZnLiCM5PwxIfCuKXx+O3AzcBehP+8Q4MscHHh1pbc5ep2s9CQpe9KYOjC/8Oju9m7LG4CPpLCf\nzCjXtAmGniRlkVdkqUC5QSxg6ElSFhl6FbjpJli6FKZOhc2bi69j6ElS9hh6FXjpJXjjDZg/H2bO\nLL6OoSdJ2WPoVaC5OTy3tcGsWcXXueMO+O1ve68GJUl9y9CrwEc/CieeCAsXQmtr8XXWrYM33+y9\nGpQk9S1DrwL79oX5d6UCD2DYsPDcWzUoSepbhl4Ftm3rCrVS7roLBgzovRqUJPUtQ68CW7fC8DJX\nD50wIVyubMSIvjkmSVJ5hl4FolR6AwZASwvs3Nk3xyRJKs/Qq8DWreVDD0I1uHVr9Y9HkhSNoVeB\nbdvKN28CHHZYGMEpScoGQ68CUZo3wUpPkrLG0KtAlIEsYKUnSVlj6FXASk+S6pOhV4GoA1ms9CQp\nWwy9mPbuhd27YfDg8uta6UlSthh6MW3fHqq8pgg3sLfSk6RsMfRiitq0CVZ6kpQ1hl5MUefogZWe\nJGWNoReTlZ4k1S9DL6Y4ld7w4VZ6kpQlhl5MUefoQWjetNKTpOww9GKK27xppSdJ2WHoxRR3IIuV\nniRlh6EXk5WeJNUvQy8mKz1Jql9phN4UYDmwAri+xDrtwNPAUiCXwj5rJs5Als5KL5+v7jFJkqJp\nSfj9ZuA24EJgLfAkMA94vts6rcAPgL8B1gBHJNxnTW3dCsccE23dAQPCY+fOaNfqlCRVV9JK7xxg\nJbAa2APcA1zaY50ZwAOEwAPYkHCfNRWneRPs15OkLEkaeuOBl7u9XlN4r7vjgVHAo8AfgY8n3GdN\nxRnIAvbrSVKWJG3ejNJbNQA4E/ggMAR4Avg9oQ+w7jz5JKxeDbNmwZw50Nra+/pWepKUHUlDby0w\nsdvriXQ1Y3Z6mdCk+Vbh8VvgNIqEXkdHx9vL7e3ttLe3Jzy89G3dCuvWwbPPwsyZMHdu7+tb6UlS\nNLlcjlwuV9V9RLgrXK9agBcIVdwrwB+A6Rw4kOVEwmCXvwEOBRYBVwLP9dhWPl8HwxyHDQv31Gtr\ng4ULy1d6Rx8dgu/oo6NVhpKkoCncuDRpTh0gaZ/eXuBaYAEhxO4lBN41hQeE6QwPA88SAu8ODg68\nutHaCh/+cLTAgzByc9kymD8/VIaSpNpJNUETqotKb+RIWLUKRo2Ktv7EibBmTfTKUJIUZLHS61fy\n+XiT0wEuvhhOO83Ak6QsMPRi2L0bDjkEBg6M/p0jjoDLLzfwJCkLDL0Y4s7Rg7D+tm3VOR5JUjyG\nXgxxr8YCYX1DT5KywdCLwUpPkuqboRdDJZWeoSdJ2WHoxRB35CYYepKUJYZeDDZvSlJ9M/RisHlT\nkuqboReDlZ4k1TdDLwb79CSpvhl6MVTavOmthSQpGwy9GCpp3hw8GHbtgn37qnNMkqToDL0YKqn0\nmppg6NBwDz5JUm0ZejFUUumB/XqSlBWGXgyVDGQBQ0+SssLQi6GS5k0w9CQpKwy9GGzelKT6ZujF\nYKUnSfXN0IvBSk+S6puhF4MDWSSpvhl6Ee3fDzt2hDl3cRl6kpQNhl5EO3bAoEHQ3Bz/u4aeJGWD\noRdRpYNYwNCTpKww9CKqdBALGHqSlBWGXkSVDmKBUCEaepJUe2mE3hRgObACuL6X9c4G9gIfTWGf\nfc7mTUmqf0lDrxm4jRB8k4HpwEkl1vsm8DDQlHCfNXHTTbBsGUydCps3x/uuoSdJ2ZA09M4BVgKr\ngT3APcClRdb7PHA/8HrC/dXMSy/Bpk0wfz7MnBnvu4aeJGVD0tAbD7zc7fWawns917kU+GHhdT7h\nPmuic6pCWxvMmhXvu9/9Ljz1VGVVoiQpPS0Jvx8lwG4Fbiis20QvzZsdHR1vL7e3t9Pe3p7s6FI0\nYwbceScsXAitrfG+u2ZNuIlsZ5U4d251jlGS6lkulyOXy1V1H0n7184DOgh9egBfAfYT+u86/bnb\nfo4AdgBXA/N6bCufz2e3CLzppjBB/aab4n/3ggvg0UdDlVhJaEpSf9TU1AQpjwNJ2rz5R+B4YBIw\nELiSg8PsGOAdhcf9wP8ssk7mJZmnN3s2tLQYeJJUa0lDby9wLbAAeA64F3geuKbwaBhJpiyMGxeu\n3TliRLrHJEmKJ2mfHsD8wqO720us+6kU9lcTSSq95uZw3c5KL1gtSUqHV2SJKMkVWcBpC5KUBYZe\nREmaN8HQk6QsMPQiStK8CYaeJGWBoReRzZuSVP8MvYi2brV5U5LqnaEXkZWeJNU/Qy8iB7JIUv0z\n9CLYuxd27w5z7Spl6ElS7Rl6EXQ2bTYluAKcoSdJtWfoRZB0EAsYepKUBYZeBEkHsYChJ0lZYOhF\nkHQQCxh6kpQFhl4ESa/GAiE0DT1Jqi1DLwKbNyWpMRh6ETiQRZIag6EXgZWeJDUGQy8CB7JIUmMw\n9CJIYyCLoSdJtWfoRWDzpiQ1BkMvgjQGsgwdGkIvn0/nmCRJ8Rl6EaRR6Q0YAM3NsGtXOsckSYrP\n0IsgjYEsYBOnJNWaoRdBGgNZwNCTpFoz9CJIo3kTDD1JqjVDL4I0BrKAoSdJtWboRWClJ0mNIY3Q\nmwIsB1YA1xf5/L8Di4FngceAU1PYZ59avx6mT4epU2Hz5sq3Y+hJUm0lDb1m4DZC8E0GpgMn9Vjn\nz8BfE8Lu68CshPvsU/k87NkDjz8O8+fDzJmVb8vQk6Taakr4/XcDNxJCD+CGwvMtJdYfCSwBJhT5\nLJ/P4MztnTthyJAQfm1tsHAhtLZWtq2TToJ9++C442DOnMq3I0n9QVNTEyTPqQMkrfTGAy93e72m\n8F4p/wP4ZcJ99qmtW0M4TZuWLPAgVHkrViSvGCVJlWlJ+P04pdn5wKeB95ZaoaOj4+3l9vZ22tvb\nKz2u1GzbBocdBnPnJt/WoEHhua0NZtVVI68kVV8ulyOXy1V1H0nLxvOADrqaN78C7Ae+2WO9U4F/\nL6y3ssS2Mtm8uWQJzJgRnpP6+tfhpz+FRYts2pSkcrLYvPlH4HhgEjAQuBKY12OdowmB9zFKB15m\npXU1FoAjj4TzzzfwJKlWkjZv7gWuBRYQRnL+GHgeuKbw+e3A/yEMYPlh4b09wDkJ99tn0pqjB47e\nlKRaSxp6APMLj+5u77b8mcKjLqV1NRYI2zH0JKl2vCJLGVZ6ktQ4DL0yDD1JahyGXhlpNm8aepJU\nW4ZeGVZ6ktQ4DL0yrPQkqXEYemVY6UlS4zD0ykgz9AYODBec3rMnne1JkuIx9MpIs3mzqclqT5Jq\nydArI81KDww9SaolQ6+MNCs9MPQkqZYMvTKs9CSpcRh6ZRh6ktQ4DL0ybN6UpMZh6PVi/37YsQOG\nDk1vm4aeJNWOodeL7dthyBA4JMW/JUNPkmrH0OtF2v15YOhJUi0Zer0w9CSpsRh6vUh7EAsYepJU\nS4ZeL6z0JKmxGHq9sNKTpMZi6PXCSk+SGouh14tqhN7w4YaeJNWKodeLH/0IFiyAqVNh8+Z0tmml\nJ0m1Y+j1Yt06eOUVmD8fZs5MZ5uGniTVjqHXi84rsbS1waxZ6Wzz5pvhhRfSrR4lSdEYer34wAfg\nzDNh4UJobU1nmy+9BLt2pVs9SpKiSSP0pgDLgRXA9SXW+X7h88XAGSnss0/s2gVf+EJ6gQddUyDS\nrB4lSdEkDb1m4DZC8E0GpgMn9VhnKnAccDwwE/hhwn32mTffhMMOS3eb99wDTU3wn/+ZbphKkspL\nGnrnACuB1cAe4B7g0h7rXALMLiwvAlqB0Qn32yeqEXojR8KoUdDcnO52JUnltST8/njg5W6v1wDn\nRlhnArC+58ZGjYKzzoL77stGFVSN0IOwzTffhCOOSH/blbj6anj+eVi5EvbsCY9DD4W9e8Nj4EDY\nty8sH3po+eWBA8O9CIstH3UUbNhw4LZLrZtkua/201/22ch/Nv8+s7vPamhK+P3LCU2bVxdef4wQ\nep/vts6DwC3AY4XXvwK+DDzVY1t5uPHtFwMGtDNiRDsnnAAjRsCcOX0fhMceC488Ep7TdPrp8K//\nGp772syZsGQJLF8eqs0tW6r3H5ckxZMrPDp9DZLnVKrOAx7u9vorHDyY5V+A/9bt9XKKN2/mIZ8f\nMiSfh4MfY8bk85s25fvUEUfk8+vXp7/d978/n//Nb9LfbjmXX57Pt7QU//uFAz8bPrx6y4cd1lj7\n6S/7bOQ/m3+f2dxnyIV0Je1ZepVQns0DdgDfA24GNnRbZz+hEvwZISTbC+v11HHZZR2MHg1//nN4\nY/hw2L07LG/bBnfcAb/5DXz4wzBoUMIjj+CrX4WODmhJ2gjcwwMPwCmnwAknpLvdUq6+Gv7+7+Hx\nx0OzAYQ/0/79oal1167Q1/jEE2FC/kknwbx58Je/pL+8di38x39UZ9u12E9/2Wcj/9n8+8zuPpcv\n/xoUyr20pFE2fgi4lRCgPwb+Ebim8NnthefOEZ7bgU9xcNMmQD6fz7N5M3zyk2GE4623wnnnwauv\nhn+kO5vhpk2DuXNTOPJe7NoVQnfXrnAsaZoxAy6+ODxXWz4PEyaEK8t0GjkScjn4xjfgn/4JvvSl\nMH0iC/2oktSpKfzjm+q/wFlqK83n8wdXsps3h36oTZvgV78K/VBnnx3+4a5mP9/rr4f/09iwofy6\ncX32s6E/77OfTX/bPZ1+OixdGgaZnHJK6J+86y4DTlL2VSP0Um64S19ra6jqOsMvl4Pf/z58NnNm\n9Sq+ao3chK7Rm9X2yCNhVOa+feH1MceEpgZJ6q/q5jJkneF32mnh9aBBofqr1vUr6z30PvGJ0IQ6\neHB43dYWRoxKUn9WN6HX6b77Qmjs3BmaO6t1/cp6D71HHw3z7bZsCX16aV4/VJLqVd2FXmsrvOc9\nYXnSpOpdv7KeQ+/VV8MDQoW3ZImBJ0lQh6EHcPfd0N4eBptcckl1btOzdWv9ht6FF4arvYwZA/ff\nb+BJUqe6DL3W1tB819wMv/tddW7T8+abXXdESFs1Q2/79nC1lXXrQrX3pS9VZz+SVI/qMvQ6/dVf\nhedq3KanXps3770XDj88LHv7Ikk6UF2H3sMPw4AB1Zl3Vq+hN2sWfO97YQK/g1ck6UB1HXrjxsFx\nx4Wh+Wn361Uz9EaMqE7oXXEFPP10mJrgFVYk6WB1HXoQbmnzX/+Vfr9ePVZ6f/pTuFbpggXVm8oh\nSfWs7kNvzJjwfMop6fZfVTP0hgwJ8ww7r5SSlo0bw7N9eZJUXN2H3t13h/l6M2ak25xXzdBragoj\nQ7duTW+bL78MhxwSmjjty5Ok4uo+9Fpbw90YFi5Md7vVDD1Iv4lz3jz4yEeyc9d5Scqiug89gIsu\ngt/+Ft73vvQGtGzdWr15ehBCb8uW9LZ3yy2hT68aE/UlqVE0ROgNGQLDhsFjj6U3oKWeKr09e8L9\n8p57rjoT9SWpUTRE6AEcdVR4TmsQRz2F3pNPwtChYdlBLJJUWsOE3ve/H/qy0hjE8ZnPwLZtcOWV\n1WsqTDP0cjn42MeckC5J5TRM6J1/PuzdC0Vuvh7b88+H54cfrl5T4eLF8NWvptMHl8vBlCnhfoMG\nniSV1jChN3AgnHdeuAB1UgMGhOdqNhVu2wYrViTvg9u9G554Av76r9M7NklqVA0TehDupP65zyWv\nnm68Mb2m0lKGDAnPSYP1iitCdTtjhqM2Jamchgq9vXth7drk1dPevXDWWdVtKpw5E445JnmwLlsW\nbifkqE1JKq+hQq9zBOdZZyWrnjZtgpEj0zmmUsaPh7PPTh6sO3aEZ0dtSlJ5DRV6c+eGUZH//M/J\nwqQvQm/kyLCfpAYODINYHLUpSeU1VOi1tsLll4c7hyfRF6E3alTXBaIr9frr4aouDz1k4ElSFA0V\negDnnAN/+EOybdRLpffkk6GJ9JCGO4uSVB1J/7kcBSwEXgQeAYrVGxOBR4FlwFLg7xLus1dnn91/\nQm/RohDykqRokobeDYTQOwH4deF1T3uAfwBOBs4D/hdwUsL9lnTKKaF5M8nFp/si9FpbQ9Pk/v2V\nb+POO+EXv/Ai05IUVdLQuwSYXVieDVxWZJ1XgWcKy9uA54FxCfdb0sCBMHhwsotP90XotbSE62Um\nuRTZ+vVhyoLTFSQpmqShNxpYX1heX3jdm0nAGcCihPvtVeegjkqH8fdF6EGyJs5XX+1adrqCJEUT\nJfQWAkuKPC7psV6+8ChlGHA/cB2h4qua666Dd7yj8mH89RB6ixeHy655kWlJiq4lwjoX9fLZemAM\noQlzLPBaifUGAA8APwV+XmpjHR0dby+3t7fT3t4e4fAO9p73wL33Vh4E9RB6zzwTKrzvfCfdY5Kk\nWsnlcuRyuaruoynh978FvAF8kzCIpZWDB7M0Efr73iAMaCkln0/jFgmEu56PHh36y1qixHo3+/eH\nfsFdu6C5OZXDKemKK8Lti6ZNi//d6dPhQx+Cq65K/7gkKQuampogeU4dIGmf3i2ESvBF4ILCawgD\nVR4qLL8X+BhwPvB04TEl4X57NXx4uMzXiy/G/+6WLeEu7NUOPEjevHn66ekejyQ1uph10EE2AhcW\nef8V4MOF5f9HDSbBn3ZaCIbJk+N9r6+aNqHy0NuxA1avhhNPTP2QJKmhNey1PDpDL66+DL1KL0U2\nY0a4Cstllzk/T5LiaNjQ+93vwjD+uBO366HSe/55byckSZVo2NB7880QKHGDoR5Cb9eu8Oz8PEmK\nJ2mfXmZ1BtcZZ8QLhnoIvUmTYOJEePBB5+dJUhwNG3p33w3HHgvf+Ea8YKiHPr0XXggX1TbwJCme\nhg29znvrrV4d73tz5oSm0aVLw3I1g6WSSm/DhjB6c8KE6hyTJDWyhu3TA3jXu2DJknjfWb8+BGVf\nDBL52tfgpZfiDbZZtiz8uZpSna4pSf1DQ4feKaeEii2OffvCc18MElm9OlwBJk7ALl0aQk+SFF/D\nNm9CV6WXz0evjMaNC3P87r+/+n1mQ4eG51NPjR6whp4kVS5LjWSpXXuzuyFDQqiMGhWtj27CBHj8\ncTj66NQP5SCbN4e7QTzwAFxwQbTvjBkTris6fnz1+xwlqZayeO3NzDv0UFi0KFoTYj4Pr70GRx3V\nN8fW2hruCLFjR7T18/kwkOXZZ52YLkmVaPjQGzEiPEfpo9u0KVSGgwZV/7g6jR594A1he7NmTdeF\nsJ2YLknxNXzoffGLoakyyo1W168PIdSXxoyJHnpLl8K73+2NYyWpUg0feueeG/rzogRELUJv9Oiw\n3yiWLIEzz4S5cw08SapEw4fe5MnhCiZ795ZfN+uV3pIlYRqGJKkyDR96Q4eGaQgrV5Zftx4qPacr\nSFLlGj70IFRHUa7MkuVKb8+eULGefHL1j0mSGlW/CL2olyPLcqW3cmWYQzhkSPWPSZIaVb8IvaiX\nI6tF6LW2ws6d8NZbva937bVhSkXcm+JKkrr0i9C7774wmbtcYNQi9JqaolV7K1fCG284KV2SkugX\nobd+faimygVGLUIPovXr7dwZnp2ULkmV6xehN2xYeJ48uXRg5PO1C70oV2UZPBimTHFSuiQl0S9C\nb86ccFWWz32udGB88pNhhOQVV/R9n9mKFfDlL5duft22LVwT9MEHDTxJSqJfhF5rawi8VatKr7N0\nafx726XlrbdC8JXa93PPwYknQktD3whKkqqvX4QelB/BuX9/eK5Fn1m5i2J7JRZJSke/Cb1yc/Uu\nvhje+c7a9Jl9+9swcmTpfRt6kpSOJKE3ClgIvAg8AvQWFc3A08CDCfaXyMSJoRlxw4bin69ZE+7I\nUIs+szPPhH37uiq+nubOhZ/8xDl6kpRUktC7gRB6JwC/Lrwu5TrgOSD9W6NH1NQU7pP3wQ8WD4+V\nK+G442pzbIcfHu6T9/rrB3+Wz4f3lyxxjp4kJZUk9C4BZheWZwOXlVhvAjAV+BEp3/Y9rpaW0ncd\nr2XoAZxwQhjM0tOKFV0DWJyjJ0nJJAm90UDndUTWF14X813gS8D+BPtKxZFHhuee4bFtG2zZEu7G\nUCvHHw8vvnjw+088ESpTbxwrScmVC72FwJIij0t6rJeneNPlxcBrhP68mlZ5AD/7GQwcCI88cmB4\nrFoFxxwDh9RwWE+pSu/xx+EDH/DGsZKUhnIzvy7q5bP1wBjgVWAsIdx6eg8hIKcCg4DDgJ8AVxXb\nYEdHx9vL7e3ttLe3lzm8eCZPhrFjw0TvkSO73l+xorZNmxAqvQceOPj9J56Az3ym749HkvpaLpcj\nl8tVdR9Jqq9vAW8A3yQMYmml98EsHwC+CHykxOf5fL7641xmzIALL4RPf7rrvXPOgXXrwrSAOXNq\nU1E99RR86lOweHHXe1u2wPjx4e4KAwb0/TFJUi01NTVByq2ESRr0biFUgi8CFxReA4wDHirxnZqN\n3uz03veGJsPuVq8OUxZqOTry+ONh2bLQlNk5unT69DCq89JLnaogSWmoeT9bN31S6V1+OTz0EJx/\nPtx9NwwdGh579oQBLrUcLDJsGGzfHpanTQvhvHZt1+u5c2tzXJJUC1mr9OrShg2waxc8/HCo6h57\nDE4+ORujIztHj7a1wQ9+EPoeO187VUGSkut3oTd0aHgeNSoEybx58Ld/m43Rkf/2b2EC/YIFYQDL\nmWdmI4wlqVH0u9CbMwcuuQS2bg3TBG69FX75y2z0mZ17bpg68fTT8PnPw8aNYQ6hJCkd/a5Pr9P4\n8fDKK12vs9Jn1tYGf/pTGMCyb194LyvHJkl9yT69FJ12WtfyGWdkp8+ss/m1M/Dsz5Ok9PTb25LO\nmRPult7UBHfdlZ0+s87QO/10mDQpW8cmSfWu3zZvZtXmzWFU6axZhp2k/q0azZuGniQpk+zTkyQp\nAUNPktRvGHqSpH7D0JMk9RuGniSp3zD0JEn9hqEnSeo3DD1JUr9h6EmS+g1DT5LUbxh6kqR+w9CT\nJPUbhp4kqd8w9CRJ/YahJ0nqNww9SVK/YehJkvoNQ0+S1G8kCb1RwELgReARoLXEeq3A/cDzwHPA\neQn2KUlSxZKE3g2E0DsB+HXhdTHfA34JnAScSgg/ZVAul6v1IQjPQ1Z4HhpTktC7BJhdWJ4NXFZk\nnRHA+4E7C6/3AlsS7FNV5I88GzwP2eB5aExJQm80sL6wvL7wuqd3AK8DdwFPAXcAQxLsU5KkipUL\nvYXAkiKPS3qsly88emoBzgT+b+F5O6WbQSVJqqqmBN9dDrQDrwJjgUeBE3usMwZ4glDxAbyPEHoX\nF9neSuDYBMcjSWosq4Dj0txgS4LvzgM+AXyz8PzzIuu8CrxMGOzyInAhsKzE9lL9g0mSlKZRwK84\neMrCOODfTUB8AAACnElEQVShbuudBjwJLAb+nTC4RZIkSZJUb6YQ+vxWANeXWOf7hc8XA2dE+G7U\nyfDqUo3z0AGsAZ4uPKakesSNJ8k5uJMwMnpJj/X9LcRXjfPQgb+FuCo9DxMJ40aWAUuBv+u2fs1/\nD82EQSmTgAHAM4SJ6d1NJUxYBzgX+H2E734L+HJh+XrgltSPvLFU6zzcCHyhSsfcaJKcAwhzXM/g\n4H9s/S3EU63z4G8hniTnYQxwemF5GPACXQMnY/0eqnHtzXMIf7DVwB7gHuDSHut0n9i+iJDMY8p8\nN8pkeHWp1nmAZKN++5Mk5wDgd8CmItv1txBPtc4D+FuIo9LzMJowKPKZwvvbCFf2Gl/kO2V/D9UI\nvfGEEZud1tB1cOXWGdfLd6NMhleXap0HgM8Tmh5+jE1rvUlyDnrjbyGeap0H8LcQR6XnYUKPdSYR\nKu9Fhdexfg/VCL1ik9SLifJ/SE0ltldqMry6pHkeuvshYd7l6cA64Nsxv9+fVHoO4vy37W+hvGqd\nB38L8aRxHoYRbmBwHaHiK7aPXvdTjdBbS+h07DSRkNa9rTOhsE6x99cWltfT1dwwFngtpeNtVGme\nh+7ffY2u/7B+RGiyUHGVnoO19M7fQjzVOg/+FuJJeh4GAA8AP+XAeeE1/z20EGbRTwIGUr6z8jy6\nOit7++636BrtcwN23pdTrfMwttv3/wGYk+5hN5Qk56DTJIoPZPG3EF21zoO/hXiSnIcm4CfAd4ts\nNxO/hw8RRtesBL5SeO+awqPTbYXPFxOuy9nbd6H0ZHiVVo3z8BPg2cL6P8f+pHKSnIO7gVeAXYR+\njk8V3ve3EF81zoO/hfgqPQ/vA/YTgrLnFBF/D5IkSZIkSZIkSZIkSZIkSZIkSZIkSZKkdP1/CX6I\nogOHP4IAAAAASUVORK5CYII=\n", + "text": [ + "" + ] + } + ], + "prompt_number": 31 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "txind = Utils.closestPoints(mesh, [0., 30.], gridLoc='CC')\n", + "q = Utils.sdiag(1/mesh.vol)*np.zeros(mesh.nC)\n", + "q[txind] = 1." + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 32 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "p = np.zeros((mesh.nC, time.size))\n", + "pn = np.zeros(mesh.nC)\n", + "p0 = np.zeros_like(pn)\n", + "un = np.zeros(mesh.nF)\n", + "u0 = np.zeros_like(un)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 34 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Proj = mesh.getInterpolationMat(np.r_[3., 30], 'CC')" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 35 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%%time\n", + "for i in range(time.size-1):\n", + " sn = ricker(400, time[i+1], 0.0025)\n", + " s0 = ricker(400, time[i], 0.0025)\n", + " pn = p0-dt*Utils.sdiag(sig)*p0+Utils.sdiag(1/rho)*dt*(Div*un+(sn-s0)/dt*q)\n", + " p0 = pn.copy()\n", + " un = u0 - dt*Utils.sdiag(AvF2CC.T*sig)*u0 + dt*Utils.sdiag(1/(AvF2CC.T*(1/mu)))*Grad*p0\n", + " u0 = un.copy()\n", + " p[:,i+1] = pn" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Wall time: 2min 55s\n" + ] + } + ], + "prompt_number": 36 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "data = Proj*p" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 37 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import JSAnimation" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 38 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import urllib" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 39 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min(), mesh.vectorCCy.max()]\n", + "extent[:2]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 40, + "text": [ + "[-124.75, 124.75]" + ] + } + ], + "prompt_number": 40 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from JSAnimation import IPython_display\n", + "from matplotlib import animation\n", + "fig, ax = subplots(1,2, figsize = (16, 8))\n", + "ax[0].set_xlabel('Easting (m)')\n", + "ax[0].set_ylabel('Depth (m)')\n", + "extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min(), mesh.vectorCCy.max()]\n", + "ax[0].set_xlim(extent[:2])\n", + "ax[0].set_ylim(extent[2:])\n", + "ax[1].plot(Utils.mkvc(data), time, 'k--',)\n", + "ax[1].set_xlim(data.min(), data.max())\n", + "ax[1].set_ylim(time.min(), 0.04)\n", + "ax[1].invert_yaxis()\n", + "\n", + "line, = ax[1].plot([], [], color=\"black\", lw=2)\n", + "nskip = 20\n", + "def animate(i_id):\n", + " icount = i_id*nskip\n", + " frame = ax[0].imshow(np.flipud(p[:,icount].reshape((500, 500), order = 'F').T), cmap = 'binary', extent=extent)\n", + " \n", + " tx = ax[0].plot(mesh.gridCC[txind,0], mesh.gridCC[txind,1], 'k.', ms = 10)\n", + " rx = ax[0].plot(10, 30, 'r.', ms = 10)\n", + " text_tx = ax[0].text(mesh.gridCC[txind,0]-15., mesh.gridCC[txind,1], 'Tx', fontsize = 18)\n", + " text_tx = ax[0].text(10+5., 30, 'Rx', fontsize = 18, color=\"red\")\n", + " ax[0].plot(np.r_[-12.5, 12.5, 12.5, -12.5, -12.5], np.r_[-12.5, -12.5, 12.5, 12.5, -12.5], 'w-', lw=2)\n", + " line.set_data([Utils.mkvc(data)[:icount]], [time[:icount]])\n", + " return frame, line\n", + "animation.FuncAnimation(fig, animate, frames=40, interval=40, blit=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Once \n", + " Loop \n", + " Reflect \n", + "
\n", + "
\n", + "\n", + "\n", + "\n" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 41, + "text": [ + "" + ] + } + ], + "prompt_number": 41 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# !ipython nbconvert --to html SeismicEx-sponge_explicit.ipynb" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 43 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/notebooks/SeismicEx-sponge_implicit.ipynb b/notebooks/SeismicEx-sponge_implicit.ipynb new file mode 100644 index 0000000..989699a --- /dev/null +++ b/notebooks/SeismicEx-sponge_implicit.ipynb @@ -0,0 +1,39426 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:97e555d042bfefd57c2df46c9164b543d1b3777851e08a95cefdb104e8f7d6e2" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from SimPEG import *\n", + "import scipy\n", + "%pylab inline" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "cs = 0.5\n", + "hx = np.ones(500)*cs\n", + "hy = np.ones(500)*cs\n", + "mesh = Mesh.TensorMesh([hx, hy], 'CC')" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Acoustic Wave equation in time domain (1st order form)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\rho\\frac{\\partial p}{\\partial t} = \\nabla \\cdot \\vec{u}+\\frac{\\partial s}{\\partial t}\\delta(\\vec{r}-\\vec{r}_s)$\n", + "\n", + "$ \\mu^{-1}\\frac{\\partial \\vec{u}}{\\partial t} = \\nabla p$" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Add damping term" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\rho[\\frac{\\partial p}{\\partial t} + \\sigma p] = \\nabla \\cdot \\vec{u}+\\frac{\\partial s}{\\partial t}\\delta(\\vec{r}-\\vec{r}_s)$\n", + "\n", + "$ \\mu^{-1}[\\frac{\\partial \\vec{u}}{\\partial t}+\\sigma \\vec{u}] = \\nabla p$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discretized form" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$ \\mathbf{diag}(\\rho)\\frac{\\mathbf{p}^{n+1}-\\mathbf{p}^{n}}{\\triangle t} + \\mathbf{diag}(\\rho\\sigma)\\mathbf{p}^n = \\mathbf{Div} \\mathbf{u}^{n}+\\mathbf{M}^{cc -1}\\frac{\\mathbf{s}^{n+1}-\\mathbf{s}^{n}}{\\triangle t}$ \n", + "\n", + "\n", + "$ \\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1})\\frac{\\mathbf{u}^{n+1}-\\mathbf{u}^{n}}{\\triangle t} + \\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1}\\sigma)\\mathbf{u}^{n+1}= \\mathbf{Grad} \\mathbf{p}^{n+1}$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "\\begin{eqnarray}\n", + " \\begin{bmatrix}\n", + " \\triangle t^{-1}\\mathbf{diag}(\\rho) & \\tilde{\\mathbf{0}} \\\\[0.3em]\n", + " -\\mathbf{Grad} & \\triangle t^{-1}\\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1}) + \\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1}\\sigma) \\\\[0.3em]\n", + " \\end{bmatrix}\n", + " \\begin{bmatrix}\n", + " \\mathbf{p}^{n+1} \\\\[0.3em]\n", + " \\mathbf{u}^{n+1} \\\\[0.3em]\n", + " \\end{bmatrix} + \n", + " \\begin{bmatrix}\n", + " -\\triangle t^{-1}\\mathbf{diag}(\\rho)+ \\mathbf{diag}(\\rho\\sigma)& -\\mathbf{Div} \\\\[0.3em]\n", + " \\tilde{\\mathbf{0}} & -\\triangle t^{-1}\\mathbf{diag}(\\mathbf{Av}^{f}_{cc}\\mu^{-1}) \\\\[0.3em]\n", + " \\end{bmatrix}\n", + " \\begin{bmatrix}\n", + " \\mathbf{p}^{n} \\\\[0.3em]\n", + " \\mathbf{u}^{n} \\\\[0.3em]\n", + " \\end{bmatrix} = \n", + " \\begin{bmatrix}\n", + " \\triangle t^{-1}\\mathbf{M}^{cc -1}(\\mathbf{s}^{n+1}-\\mathbf{s}^{n}) \\\\[0.3em]\n", + " \\tilde{\\mathbf{0}} \\\\[0.3em]\n", + " \\end{bmatrix}\n", + "\\end{eqnarray}" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "mesh.setCellGradBC('dirichlet')\n", + "Grad = mesh.cellGrad\n", + "Div = mesh.faceDiv\n", + "rho = np.ones(mesh.nC)*2.7\n", + "vhalf = 2000\n", + "vblk = 2800\n", + "v = np.ones(mesh.nC)*vhalf\n", + "blkind = np.logical_and(mesh.gridCC[:,1]>-12.5, mesh.gridCC[:,1]<12.5) & np.logical_and(mesh.gridCC[:,0]>-12.5, mesh.gridCC[:,0]<12.5)\n", + "v[blkind] = vblk\n", + "mu = rho*v**2\n", + "AvF2CC = mesh.aveF2CC" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "dat = mesh.plotImage(v)\n", + "plt.colorbar(dat[0])" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 5, + "text": [ + "" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEPCAYAAABFpK+YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHN1JREFUeJzt3X2wXOV92PHvtbiAjHAwYw/ojV6RiBYUCLKMwJETBDVE\ndFJeyiSYFExaxmYqDAp1zYvcFhE3LiY1INqipIMoCAVSKtkMHl4MZiQMGUCISqAXFIQjMehawhhI\nhAJGL9z+8XvW99y9u3t3de7VnnP3+5k5s2efc3b3ORfxO8/+zrPnB5IkSZIkSZIkSZIkSZIkSZIk\ndaTJwApgA7AeuDq1zwRWAWuAF4FTMq+5AdgMbALOzrTPANalbQtHtNeSpP12NHByWh8H/C1wPLAS\n+L3Ufg5xcgA4AVgLdAM9wOtAV9q2ijhhADwKzKn3oZ8Yjp5LkvbLDiKQA+wCXgUmAtuBX0vtRwC9\naf084AFgD7CVCPynAuOBw4ngD7AEOL/ehx40XL2XJOXSA0wHnifSNc8C/40YoH8h7TMhba/YRpwo\n9qT1it7UXpMjfklqv3HAMmAeMfJfTOT7jwGuAe4ezg8bhSP+f9IHb7S7E5LK4Wlg9v6++FDo+2Vr\nL3kf+FRVWzewHFgKPJTaZgJfSuvLgLvSei9xQbhiEjHS703r2fZe6hiFgf8NYMEB+qwVwBkH6LMO\nlNF4TDA6j2s0HhMc2ONacHqeV/8S+C8t7P8fIw+f1UWM7jcCt2faXwdOJ05MZwKvpfaHgfuBW4lU\nzlQir98H7CTy/auAS4E76vVjFAZ+STpwuvO9fBZwCfAKMXUTYD7wNeB/AocAH6bnECeIB9PjXmAu\nEfRJ6/cAY4lZPY/X+1ADvyTlkDOIPkv9a62n1mn/TlqqvQSc2MyHGvhz6Wl3B0ZAT7s7MEJ62t2B\nEdDT7g6MkJ52d6AlY9vdgf1g4M9lSrs7MAJG4zHB6Dyu0XhMULbjypnqaQsDvyTlUMYgWsY+S1Jh\nOOKXpA5TxiBaxj5LUmE44pekDmPgl6QO43ROSeowZQyiZeyzJBWGqR5J6jBlDKJl7LMkFYYjfknq\nMGUMomXssyQVRhlH/O0uvXg38BawLtN2JPAkUXjgCaLQcMUNRC3KTcDZB6iPklTX2BaWGiYTlWc2\nAOuJcosA/4e4P/8aYAv99+o/C1hN3L9/NQMr1swgYulmYGGjPrc78P9vYE5V2/VE4D8OeCo9BzgB\nuCg9zgHupP39l9ThultYathD1NSdBpwGXAkcT8S66WlZnhaAt4HfB04CLgPuy7zXIuByoirXVAbH\n1l9pd+B8Bnivqu1c4N60fi9wflo/D3iA+ENtJUqTzRz5LkpSfQe1sNSwA1ib1ncBrwITMtu7gD8k\nYh9p3x1pfSPxRaIbGE+UdVyVti2hP3bW7HPRHEWkf0iPR6X1CcDzmf22ETUnJaltuluJonsbbu0h\nRvgvZNp+h4iDP62x/4VE1a09RCzcltnWS4P4WMTAn9VHfz3JettrWJFZ76FshR0kjZQtRMJg+BzU\nIIo+sw+e/biptxkHLAPmESP/iouJ4urVpgE3Ezn/lhUx8L8FHE18nRkP/Dy19xIXQiompbYazqjd\nLKnDTWHgQPDp3O/YPab+tjPHwJmZ5zf/Y+23IHL4S4GHMu0HARcAn6vafxLwfeBS4kwGEQsnVe1T\nJz62P8dfy8PERQvS40OZ9i8DBxP/5abSn8+SpLY46KDmlxq6gMVEvv72qm1fInL+P8u0HQE8AlwH\nPJdp3w7sJAq0dxEnhexJZGCfWzrC4fcAcDrwGeBN4D8TX18eJK5ObyUubED8YR5Mj3uBuTROA0nS\niOs+JNfLZwGXENMzK1M2bwAeJ2b2PFC1/9eBXwduTAtEuucXREy8h7jg+2h6j5q6cnW5mPpgQbv7\nIKkUFkC+ONjXN2HonSq6Yuze9rjb7hG/JJVbCaNoCbssSQVSwihawi5LUoE0mNVTVAZ+ScqjhFG0\nhF2WpALJN6unLQz8kpRHCaNoCbssSQVSwihawi5LUoF4cVeSOkwJo2gJuyxJBVLCKFrCLktSgZQw\nipawy5JUIE7nlKQOU8IoWsT78UtSeYxpYRlsMlEycAOwHrg6s+0q4n7864HvVr3uGKJS1zcybTOA\ndcBmYGGjLpfwXCVJBZIviu4BriGKqI8jaug+SVQhPBc4Ke3z2arX3UoUZMlaRNQxWUXcj38Ode7J\nb+CXpDzyRdEdaYEYwb9KFEn/KvBfiaAP8HbmNecDfwdkCzmOBw6nvyrhkrRfzcBvqkeS8siX6snq\nAaYDLwDHAb8LPA+sBD6f9hkHXMvgalMTgW2Z572prSZH/JKUR4MouvJnsHJ7U+8yDlgGzAPeT+/6\naeA04BSi7OyxRMC/DfiAHJW8DPySlMeh9TfNPjaWipvW1NytG1gOLKW/QPo24Ptp/UXgY6I2+Uzg\nQuAWovD6x8CHad9JmfecRIz6azLwS1Ie+e7V0wUsBjYCt2faHwLOBJ4m0j4HEwXVfzezz43Et4M7\n0/OdwKlEnv9S4I56H2rgl6Q88kXRWcAlwCtA5fvADcDdaVkH7Aa+0sR7zQXuAcYSs3pqXtgFA78k\n5ZMvij5L/Uk2lw7x2puqnr8EnNjMhxr4JSkPb8ssSR2mhFG0hF2WpAIpYRQtYZclqUC8O6ckdZgS\nRtESdlmSCqSEUbSEXZakAnFWjyR1mBJG0RJ2WZIKpIRRtIRdlqQCMdUjSR2mwd05i8rAL0l5lDCK\nlrDLklQgJUz1FLn04lb6b1VaqSN5JFGI+DXgCaIQgSS1z0EtLINNBlYAG4D1wNWpfQFRjGVNWuZk\nXnMS8Fza/xXiXv0AM4jbOG8GFjbqcpEDfx8wm6hBOTO1XU8E/uOAp9JzSWqffIF/D3ANMI0os3gl\ncDwR/24l4t90+u+tfxBwH/A14DeB04G9adsi4HJgalqyJ4sBihz4YXBNyXOBe9P6vUQVeUlqn3zF\n1ncAa9P6LuBV+ouk16qpezYxyl+Xnr9HlF8cDxxOf3ZkCQ3iY5EDfx/wY2A18NXUdhTwVlp/Kz2X\npPY5tIWlsR5idP98en4V8DJRmrGS1p5KxMbHicIr30ztE4nUUEUv/SeQQYp8cXcWsB34LJHe2VS1\nvS8tktQ+w3NxdxywDJhHjPwXAX+atn0b+B6RxukGvgh8niiy/hRxAviHVj6syIF/e3p8G/gBked/\nCzia+Ho0Hvh57ZeuyKz3AFNGqIuSymULMW9kGDWIoitfgpX/b8h36AaWA0uJIuswMLbdBfwwrb8J\n/AR4Nz1/FPhceu2kzGsmEaP+mmrlkIrgk8R59H3gMGIGz03Al4B3gO8SF3aPYPAF3r64IC5JQ1kA\n+eJgX9/q5nfu+jzVn9dFXK98h7jIWzGe/sHvNcApwB8BnyZS4F8kLgw/RlwEfgx4gZgVtAp4BLiD\nOgXXizriP4oY5UP08a+I4L8aeJD4yrMV+MN2dE6SfiVfFJ0FXEL/1HWA+cDFwMlEOnsLcEXa9h4R\n6F9M2x4hgj7AXOAeYCzxTaBm0IfijvjzcMQvqUkLIO+If93QO1V0nUjezxsWRR3xS1I5lDCKlrDL\nklQg1tyVpA5Twihawi5LUoGUMIqWsMuSVCAljKIl7LIkFUdfCW/LbOCXpBz2lTCKlrDLklQcBn5J\n6jAfHXLw0Dv9yu4R60crDPySlMO+MeVL8hv4JSmHfSUsumvgl6Qc9hr4Jamz7CthGC1y6UVJKrx9\njGl6qWEyUTlqA7CeuJ9+1jeImrpHpueHAg8Qt3HeyMB6JDOIWrybgYWN+mzgl6Qccgb+PUShlWnA\nacCVwPFp22TgLOCNzP5fTo8nEYH+CuCY1LaIqFUyNS1z6vXZwC9JOXzEwU0vNewA1qb1XcCrwIT0\n/Fbg2qr9txNVCcekx93ATqJi1+FE9S2AJcD59fpcvuSUJBXIMOb4e4DpRAnF84BtREon60fApcQJ\n4JPAnwB/D/xG2r+iF5hY74MM/JKUwzBN5xwHLAPmETn9+USap6JStesSorTieCLv/wzwVKsfZuCX\npBwaBf7VK/+R1Ss/GOotuoHlwFLgIeBEYvT/cto+CXgJOBX4baIe+T7gbeBviFz/s2k/Mq/prfeB\nBn5JyqHRPP6TZ3+Kk2d/6lfP/9dNv6jepQtYTMzQuT21rQOOyuyzhQju7wKbgDOJk8RhxAXh24hr\nBTuJk8MqIh10R71+eXFXknLYx0FNLzXMItI3ZwBr0nJO1T59mfW/BA4mTg6rgLuJaaAAc4G7iOmc\nrwOP1+uzI35JyiFnjv9Zhh6AH5tZ/4g4UdTyEpEmGpKBX5Jy2F17mmahGfglKQfv1SNJHaaM9+op\nX48lqUC8LbMkdRgDvyR1GHP8ktRhdnNIu7vQMgO/JOVgqkcaYQtYMCo/S+VlqkeSOozTOSWpw5jq\nkaQOY+CXpA5j4JekDvNRCadzlvF+/HOIYgSbgeva3BdJHW4fY5peapgMrAA2EPfVv7pq+zeIUoxH\nZtpuIOLfJuDsTPsM4j79m4GFjfpctsA/BvgfRPA/AbgYOL6tPZLU0XIG/j3ANcA0oprWlfTHtMlE\n3d03MvufAFyUHucAd9Jfj3cRcDkwNS1z6vW5mcB/NfDpJvY7EGYSlWW2En+wvyaq0UtSW+xlTNNL\nDTuAtWl9F/AqMCE9vxW4tmr/84AHiPi3lYiHpxLF1w8nqnIBLAHOr9fnZgL/UcCLwIPEGaSr8e4j\naiLwZub5ttQmSW2Rs/RiVg8wHXiBCPDbgFeq9pmQ2isqMbC6vZcGsbGZi7vfAv4TkUv6YyLV8iBR\nIPinTbx+OPUNvQtEyqyiB5gyAl2RVD5biIHy8Gk0q2fryjd4Y+UbdbdnjAOWAfOInP58Is1TMawD\n7mZn9XxMfCV5C9hHpH6WAT8GvjmcHRpCL5H3qpjMwLNccsYB6o6kcpnCwIHg07nfsVHgnzz7WCbP\n7i+Z+5Obnq21WzewHFgKPETUze0BXk7bJxH1dE9lcAycRMTA3rSebe+t169mAv884CvAO0QF9/9A\n5Jc+QVw9PpCBfzVx0aIH+BlxkePiA/j5kjTAR/lq7nYR2ZONwO2pbR2RYq/YQszYeRd4GLifyP9P\nJOLhKiIbspM4OawCLgXuqPehzQT+I4F/xcAryxDfAv5lE68fTnuBrwM/Imb4LCYuhkhSW+S8V88s\n4BIil78mtc0HHsvsk01xbyRS7RuJeDg3s30ucA8wFngUeLzehzbT4xsbbNvYxOuH22MM/KNIUtvk\n/OXusww9yebYquffSUu1l4g00ZD85a4k5eAtGySpw3g/fknqMN6PX5I6jKkeSeowu/NN52wLA78k\n5WCOX5I6jDl+aYQtYEG7uyANYI5fkjqMgV+SOow5fknqMOb4JanDOJ1TkjpMGVM9ZSu2LkmFkrP0\n4mSiZOAGYD1R4xzg20QhlrXAU/QXXzmLqEvySnrMVp2aQdzLfzOwsFGfDfySlMM+xjS91LAHuAaY\nBpwGXAkcD9wC/BZwMlGVq3J7/LeB3wdOAi4D7su81yLgcqI4y1SiRnpNpnokKYec0zl3pAVgF1FY\nagIDC0yNA36R1tdm2jcSRVe6gc8AhxPVtwCWAOdTpxiLgV+SchjGefw9wHTghfT8z4gSih8Q3waq\nXUgUX9lDlGHM1h/vTW01meqRpBw+4pCmlwbGAcuIGue7Utu3gGOIcoq3Ve0/DbgZuGJ/+uyIX5Jy\naDTi/2Dli3ywcvVQb9ENLAeWEvn8avcTNXQrJgHfJ74NbEltvak9u09vvQ808EtSDo0C/yGzT+OQ\n2f1Zmndv+ovqXbqAxUS+/vZM+1Ridg7AefQXYj8CeAS4Dngus/92YCdwKpHnvxS4o16/DPySlEPO\nefyzgEuI6ZmV4D6fmJ3zT4F9wE+Bf5e2fR34dWKWT2Wmz1nExd+5RFpoLPENoeaFXYizzWjTh3dw\nlNSUBZAvDvZN6ts89F7Jtq6peT9vWDjil6QcvDunJHUYA78kdZiPdnuTNknqKPv2li+Mlq/HklQg\n+/aa6pGkjmLgl6QOs3ePgV+SOsrH+8oXRsvXY0kqElM9ktRhflm+MFq+HktSkextdwdaZ+CXpDwM\n/JLUYUoY+ItYgWsBUUJsTVrOyWy7gbhH9Sbg7APeM0mqtqeFZbDJwApgA7AeuDq1/zlRd/dloujK\nr1W97hiiUtc3Mm0zgHVEjFzYqMtFDPx9wK1E7cnpwGOp/QTgovQ4B7iTYvZfUifZ18Iy2B7gGqKU\n4mnAlcDxwBOp7beA14hBb9atREGWrEXEffynpmVOvS4XNXDWul/1ecADxB9qK/A6MPMA9kmSBtvb\nwjLYDmBtWt9FjPInAE8CH6f2FxhYVvF84O+Iql0V44HDiepbAEvSfjUVNfBfRXzFWUyUGoP4Y2Sr\nyG+jQRV5STogftnC0lgPkeV4oar939Jfc3cccC2Dq01NZGB87KVBfGxX4H+SyEVVL+cSX1emACcT\ndSS/1+B9+ka2m5I0hHwj/opxwDJgHjHyr/gWsJsouA4R8G8DPiBHJa92zeo5q8n97gJ+mNZ7iQsh\nFQ2qyK/IrPcQ5xFJ2kJkiodRo4C+biWsXznUO3QDy4GlwEOZ9j8G/gXwzzNtM4ELgVuIbMjHwIfE\nBeBsOqhBfCxA7ccaxhMjfYiLHqcAf0Rc1L2fOPCJwI+B32DwqN+au5KatABy1txleQuJhwu7qj+v\nC7gXeIeIdxVziGzH6UQh9VpuBN4nLvRCpIiuJvL8jwB3UKfgehHn8X+XSPP0EafnK1L7RuDB9LiX\nqChvqkdSe9WeptmsWcAlwCvE9HWA+UTQPphIiwM8R8S8RuYC9wBjiWsCNYM+FHPEn5cjfklNWgB5\nR/x/1cL4818PGvG3RRFH/JJUHiX85a6BX5LyGHqaZuEY+CUpD0f8ktRhDPyS1GEM/JLUYfJN52wL\nA78k5VH7rpuFZuCXpDyc1SNJHcYcvyR1GHP8ktRhzPFLUocx1SNJHcbAL0kdpoQ5/qLW3JWkcvio\nhWWwyUTJwA3AeqKQCsAfpLZ9wOeqXnMScX/+9cR9/A9O7TOIErabgYWNumzgl6Q88tXc3UNU3poG\nnAZcCRxPBPALgJ9U7X8QcB/wNeA3iQpdlXdeBFwOTE3LnHpdNtUjSXnkS/XsSAtEkfVXgQnAU3X2\nP5sY5a9Lz99Lj+OBw4myiwBLgPOpU4XLEb8k5bGvhaWxHmA6UTu3nqlEydnHgZeAb6b2icC2zH69\nqa0mR/ySlEejWT2/WAnvrGzmXcYBy4B5xMi/nm7gi8DngQ+JbwYvAf/QzIdUGPglKY9Ggf+I2bFU\nvHZTrb26geXAUuChIT7tTSLv/256/ihx8XcpMCmz3yRi1F+TqR5JymNPC8tgXcBiYCNwe51PyBZn\n/xFwIjCWGLifTsz+2QHsBE5N+19Kg5OII35JyqP2NM1mzQIuIS7Yrklt84FDgP8OfAZ4JG07B/h7\n4FbgRSLX/wjwWHrdXOAe4qTwKHUu7MLAM8lo0QcL2t0HSaWwAPLFwT6+0Nf83s915f28YeGIX5Ly\nKOEvdw38kpSHd+eUpA7jTdokqcMY+CWpw5jjl6QOk286Z1sY+CUpD1M9ktRhTPVIUodxOqckdRhT\nPZLUYQz8ktRhzPFLUocp4Yi/Xffjb1RB/gaiSvwmor5kRdMV5CWpJCYDK4h4uB64OrUfCTwJvAY8\nARyR2g8FHiBu47wRuD7zXk3HyHYF/noV5E8ALkqPc4A76b+FadMV5CWpJPYA1wDTgNOAK4HjiYD+\nJHAcUV6xEuC/nB5PIgL9FcAxqa3pGNmuwL+JOJNVO484m+0BtgKvExVl6lWQl6Qy2wGsTeu7gFeJ\nIunnAvem9nvpj3fbgcOAMelxN1F5q6UYWbTSixMYWCl+G/FHqG5vWEFekg6cfLUXM3qA6cALwFHA\nW6n9rfQcovTiTuIEsBX4c6Iq10RaiJEjeXH3SeDoGu3zgR+O4OcSKbOKHmDKyH6cpJLYQsTL4dTo\n6u5PGJzRrmkcUXB9HvB+1ba+tECUaRxLjPCPBJ4hUkEtGcnAf9Z+vKaXuNhRMYk4i/XSQgV5OGM/\nPlrS6DeFgQPBp4fhPRuN5L+Qlorv1Nqpmwj699FfIP0tYuC8gwjyP0/tvw38gJgY8zbwN0Su/1la\niJFFSPVk608+TFy8OJj4rzOVyFm1VEFekg6cD1tYBukCFhMzdG7PtD8MXJbWL6M/3m0CzkzrhxEX\nhDfRYoxsV+C/AHiT6HS2SvxG4MH0+BhRNb7yFWcucBcxVel1GlSQl6QDJ1eOfxaRvjkDWJOWOcDN\nRNbkNSLQ35z2/0tiYLyOGBTfTUwDhRZiZNurvY+APljQ7j5IKoUFkC8O9sV1g2ZNyft5w8Jf7kpS\nLuW7Z4OBX5JyKd89Gwz8kpSLI35J6jA1Z+sUmoFfknIx1SNJHcZUjyR1GEf8ktRhHPFLUodxxC9J\nHcYRvyR1GKdzSlKHccQvSR2mfDn+ItyPX5JKLNdtmScTJQM3ELdXvjq1H0lUMXwNeAI4IvOaG4hb\nL28Czs60zyBu17wZWNioxwb+XFq5HWtZjMZjgtF5XKPxmKB8x7W3hWWQPcA1wDSiPsmVwPHA9UTg\nP44orXh92v8E4KL0OAe4k/7bPC8CLicKWE1N22sy8Oeytd0dGAFb292BEbK13R0YAVvb3YERsrXd\nHWhRrhH/DmBtWt8FvEoUST8XuDe13wucn9bPAx5Ib7aVKLhyKlGe8XCiOAvAksxrBjHHL0m5DFuO\nvweYDrwAHEXU3SU9HpXWJwDPZ16zjThR7EnrFb2pvSYDvyTlMizTOccRBdfnAe9XbeujvwSt6lhJ\n/x/KxcXFpdGyknxa/bydNd6jG/gR8CeZtk3A0Wl9fHoOkeu/PrPf40Sq52giTVRxMfAX+3lMkqQR\n1EXk42+rar8FuC6tX09/sfUTiGsCBxMFfH9K/8XdF4iTQBfwKA0u7kqS2ueLwMdEMF+TljnEdM4f\nU3s653ziou4m4Pcy7ZXpnK8Dd4x0xyVJGnX+gPiBxT7gc1Xbcv+YoiAWELMCKqOOczLb6h1jGcwh\n+r2Z/q/OZbUVeIX471OZttfohz5FdDcxS2Vdpm1/fqwkjbh/RvyQYgUDA38l39ZNTMV6nf582ypg\nZlovQ77tRuDf12ivdYxl+f3HGKK/PUT/1xI/jimrLUSQzLoFuDatX0d/LriofoeYspgN/PWOocz/\n9grNP2JzNhGjkWrD8mOKAumq0VbrGGfW2K+IZhL93Ur0/6+J4ymz6v9G9X7oU1TPAO9VtbXyY6Wy\n/NsrNAN/PhMY+KOJyo8pqtsb/piiQK4CXgYW0/91u94xlsFE4M3M8zL1vZY+4oLfauCrqa3eD33K\npNGPlcr6b6/Q/AFXvyfpnzebNR/44QHuy0ipd4zfIu7z8afp+beB7xH3/ailb/i7NiLK0s9mzQK2\nA58l/ltuqtpemSteZkMdQ9mPrxAM/P3O2o/X9BJ316uYRIxKetN6tr13/7s2bJo9xrvoP9nVOsYi\nHEszqvs+mYEjyLLZnh7fBn5ApD3eIk7mO4gU48/b07Vc6h1Dmf/tFZqpntZlc6wPA1+m/8cUU4m8\n/g7iF3qVH1NcCjx0YLvZsvGZ9Qvov/hW7xjLYDXR3x6i/xcRx1NGnySuGwEcRsxwWUccz2Wp/TKK\n/++slnrHUOZ/exoFLiByxR8SQf2xzLbR8mOKJcRUwZeJ//GyueJ6x1gG5wB/S/T/hjb3JY8pxAyX\ntcR92yvH0uiHPkX0APAzYDfx/9S/Yf9+rCRJkiRJkiRJkiRJkiRJkiRJkiRJkqQD6RTil8iHELc5\nWE/c413qaLXuvy6NJt8GDgXGErcI+G57uyNJGmndxKj/eRzoSIB359To9xkizTOOGPVLHc8RkEa7\nh4H7gWOJW09f1d7uSJJG0leA/5vWP0Gke2a3rTeSJEmSJEmSJEmSJEmSJEmSJEmSJGl0+/9BB90c\njwdxcwAAAABJRU5ErkJggg==\n", + "text": [ + "" + ] + } + ], + "prompt_number": 5 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "time = np.linspace(0, 0.08, 2**10)\n", + "dt = time[1]-time[0]" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 7 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print dt" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "7.82013685239e-05\n" + ] + } + ], + "prompt_number": 8 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "ax = mesh.vectorCCx[-30]\n", + "ay = mesh.vectorCCy[-30]\n", + "indy = np.logical_or(mesh.gridCC[:,1]<=-ay, mesh.gridCC[:,1]>=ay)\n", + "indx = np.logical_or(mesh.gridCC[:,0]<=-ax, mesh.gridCC[:,0]>=ax)\n", + "tempx = zeros_like(mesh.gridCC[:,0])\n", + "tempx[indx] = (abs(mesh.gridCC[:,0][indx])-ax)**2\n", + "tempx[indx] = tempx[indx]-tempx[indx].min()\n", + "tempx[indx] = tempx[indx]/tempx[indx].max()\n", + "tempy = zeros_like(mesh.gridCC[:,1])\n", + "tempy[indy] = (abs(mesh.gridCC[:,1][indy])-ay)**2\n", + "tempy[indy] = tempy[indy]-tempy[indy].min()\n", + "tempy[indy] = tempy[indy]/tempy[indy].max()\n", + "temp = tempx+tempy\n", + "temp[temp>1.] = 1.\n", + "f = 1-temp*0.1\n", + "sig = (1.-f)/f*2./dt" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 9 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "An = sp.vstack((sp.hstack((1/dt*Utils.sdiag(rho), Utils.spzeros(mesh.nC, mesh.nF))),\\\n", + " sp.hstack((-Grad, 1/dt*Utils.sdiag(AvF2CC.T*(1/mu)) + Utils.sdiag(AvF2CC.T*(1/mu*sig))))))\n", + "Bn = sp.vstack((sp.hstack((-1/dt*Utils.sdiag(rho)+Utils.sdiag(rho*sig), -Div)),\\\n", + " sp.hstack((Utils.spzeros(mesh.nF, mesh.nC), -1/dt*Utils.sdiag(AvF2CC.T*(1/mu))))))" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 10 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fig, ax = plt.subplots(1,2, figsize = (10, 5))\n", + "ax[0].spy(An, ms = 1)\n", + "ax[0].grid(True)\n", + "ax[1].spy(Bn, ms = 1)\n", + "ax[1].grid(True)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAEaCAYAAAC7JJutAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2QXNV55/Hv2LJl47AMIl5eBUOtpcSKSQZNGRwn2CNj\nQM7aQP4BuWoLDSbLGuwFbZwYQbxgr3EMbGl5SRU4GL8AZfMSs07WZcDIpuVsZRfJyG6HDSYSkHGh\nwQgsZ1j5JVmIe/84p+nbVz0zt7tP3/ucc3+fqq7pe/vt0dO3nznnnqdHICIiIiIiIiIiIiIiIiIi\nIiIiIiIiIiIiIiIiIiIiIiIiIiVZDzwB7AYuG/A5Pg/sBR7L7FsBbAV2AQ8B45nbLvev9wRwemb/\nlH+O3cCNmf3LgXv8/keA4zK3bfSv8TTwOPB3wP8BLqkojt3AU0DTx/PpiuLYBZwHvBr4HvC1CuN4\nAfhbH8eOCuP4IPAV4Ae49+bkkuP4B+CfgFmfixdxx2lVx0YKUqlfu4BNQAPVsGxONqL6ZaV+7cLV\nsFmfC9WwEXk18CQwAbwG90F88wDPcwpwIt2F7Trgo/76ZcA1/voa/zqv8a/7JDDmb9sBnOSv348r\nugAXAzf76+cCd/vrK3AFZBxYDTzjr/8K8Pf+31J2HOO4QeI4sAx3YP1uRXE8BVwBfAn4H/4+VcTx\nkn/OrCri2A98yN+2DDikojieAg4FfgSsrDCGbPGMUUr1axz3C+8Uf5tqmLv8GLgX1S9r9WsceBWq\nYSPx28CDme3N/jKICboL2xPA4f76EX4b3Og5O7N9EHgbcCRuJtC2AfhM5j4n++vLcDMYgPcDt2Qe\n8xn/OIC/BN5dcRwHAd8BfqOiOO7EvSfr6Mw4q4jj/wJ/QLey4zjEx9E+PqqKA3//TwH/s+IY8rmI\nTcr1C1TDjgH2AP8F1S9r9WsD7uyX2Rr2qoVuiMDRuLNMbXv8vhAOxy0P4H+237Sj/OvkXzO/fy4T\nSzbOl3GnTw9b5LkmcDPg7RXFMQds8a/XXpqoIo4TcaeVf5nZV0UcLwOfAB4F/n1FcRzv9/8h8F3g\ns8AbKsrHHuA04K6KcpF9rpilWr9ANQzgeuCv/P421S8b9eto3KDIbA2LeWDWKvF1ynqt1wL3AZfi\nTv1WEUcL+G+4Gd87cGesyo7jvcBPgWfpnDrOKysft+Ly8R7cqfhTcreXEccy3Af7fwFrgZ9x4NmV\nsvLxatwp/r/ocVuZn5XYpVi/wC1j1r2GvRd4HrdUpvplq36Bq2Hvw3ANi3lgNodbH25bSfeIdBh7\ncac0wZ22fH6B12yfrp7z1/P724851l9vr63v6/Fcx+FG8XfilgGqiqOdxxeBr+OaHMuO4+3AW3Cn\nke8C3oXLSxX5ONQ/5gXgq7i+grLj2IMrZv/bb38FV+CeqyAfv43rtWifoq/yGI1ZavVrJW4gch+q\nYW8HzsT1Lb0f1S9L9Wsl8EZgJ6phI7EM10A3gTvTNGjzLBzYo3EdnbXlzRzYDPha3OnZp+jMiLbj\n1pbHOLAZsL22vIHuZsB2k+qhuNllu2mwqjj+Da6Jdxx4PfDXwKkV5aN9/Z10ejTKjuMo3Dd5xnGn\n3v8G15tQRT5+gfsFA/BxH0MVcfwUuIiOKo+NmKVWv57GTaKuz8VW9xr2b1H9slS/nsZNHjbSoRoW\n2Htw3/55EteoN4i7cEtm/w+3Nnw+LonfpPfXZ6/wr/cEcEZmf/vrs08CN2X2L8d9M6f99dmJzG3n\n+/3P4HqpmnS+yru+gjh+iCtqTdxXrP/Y3152HLvpfHDeSedbTWXH8Q+4nDRxfwKgfYxVkY+P4RqZ\nvw/8d9xMrOw4nsRNIA7O3FblsRG7VOrXbuBqVMPyOdmI6pel+rUb12f3Y1TDRERERCQVIf4Io4hI\nVVTDRCQZof4Io4hIFVTDRKQv1r+VeRKuqM3i/oLx3cBZVQYkItIH1TAR6Yv1gdko/wijiMioqYaJ\nSF+sD8wq/0NvIiJDUA0Tkb4sqzqAJRT4o2zLWu5/PhCRGvk+MFl1EAWoholI3qL1y/oZs0eBVXT+\nCOO5dP4mjPcy7o8Iu0ur1arksnHjxspeW3HYjUFxjCYO4LdKrEPD6KuGxfyeWIgj+7tgmHymko8Q\nF/dnvvQ7NmQcLFG/rA/MXgY+DHwDeBy4h+7/1f0AY2NTi90sIlKmvmqY6tdwWq2dXdvK5/A2bnxf\n17ZyOnrWB2YADwC/BrwJ+HSvO1j4ME5MTJT+mr0oDlsxgOLIsxJHiZasYVmqX8MJ8fsgpXyEiEG/\nYzvKiCOGgVkhVR8409PTpb7eQhSHrRhAceRZicMS1S8nVBzD5jO1fISIQceoU0YcyQzMoPoDR0Rk\nUKpfYSmf4Smn5UhqYAY6cEQkXqpfYSmf4SmnozdWdQABtPy3HLrkD5b8wSQi8RobG4M06hf0qGGq\nX2Epn+Epp4Nbqn4ld8asTaN6EYmV6ldYymd4yunoJDswg3IPnG3bto3sufuhOGzFAIojz0oc1ql+\nhdVPPuuQjxAx6BgdjaQHZqBRvYjES/UrLOUzPOU0vBR6NHr2mOVpPVwkHan3mOWpfoWlfIannBZX\n2x6zPI3qRSRWql9hKZ/hKafh1GZgBqM9cOq0/l2EhTgsxACKI89KHLFR/QprsXzWMR8hYtAxGkat\nBmagUb2IxEv1KyzlMzzldHgp9GgU6jHL03q4SLzq1mOWp/oVlvIZnnK6MPWYLUCjehGJlepXWMpn\neMrp4Go7MIOwB06d1r+LsBCHhRhAceRZiSN2ql9hdedzv4mBhIX3ZZgYdIwOptYDM9CoXkTipfoV\nlvIZnnLavxR6NAbqMcvTerhIPOreY5an+hWW8hmectqhHrOCNKoXkVipfoWlfIannBangVnGMAdO\nnda/i7AQh4UYQHHkWYkjNapf4Wzbts3EQMJCPkLGoGO0GA3Mcix8GEVEBqH6FZbyGZ5yurQUejSC\n9JjlaT1cxC71mC1O9Sss5TO8OudUPWYD0qheRGKl+hWW8hmecrowDcwW0c+BU6f17yIsxGEhBlAc\neVbiSJ3q1+B6xVHFQMJCPkYZg47R3jQwW4JG9SISK9WvsJTP8JTTA6XQozGSHrO8Oq+Hi1ijHrP+\nqH6FpXyGV6ecqscsEI3qRSRWql9hKZ/hKacdRQZmnwf2Ao9l9q0AtgK7gIeA8cxtlwO7gSeA0zP7\np/xz7AZuzOxfDtzj9z8CHJe5baN/jV3AeQViHanFDpw6rX8XYSEOCzGA4sgrOQ7VL0/1q7gicZQx\nkLCQjzJj0DHqFBmYfQFYn9u3GVfYVgPf8tsAa4Bz/c/1wM10TtfdAlwArPKX9nNeAOzz+64HrvX7\nVwBXAif5y1V0F9BKaFQvEhXVrwzVr7CUz/CU0+I9GhPA14AT/PYTwDtxM9EjgG3Ar+Nmm7+kU5we\nBD4O/BB4GHiz378BmAY+6O9zFbAdWAb8CHgj8H7gHcBF/jGf8a9zdy62UnrM8uq0Hi5iTZ89ZhPY\nrV9QQQ1T/QpL+Qwv5ZyOqsfscFxRw/883F8/CtiTud8e4Oge++f8fvzPZ/z1l4EXgcMWeS4TNKoX\niZbql+pXUMpneHXOaYjm/5a/1E73gbPfxIFTp3X4GGIAxZFnJQ5P9QtQ/eo2SByjGEhYyEeVMdT1\nGF024OPaSwDPAUcCz/v9c8DKzP2Owc0U5/z1/P72Y44FnvXxHILr2ZjDLRe0rcQtJxxgZmaGiYkJ\nAMbHx5mcnGR62j20ncRRbTcaW1i37sJXYhkbW02jcWtpr5/fbjabpb7eQtttVb2+pe1ms2kqnqq3\nB8lH+/rs7CwBmKpfUF0Na7V2Mja2Gvg5cDBjY1M0GltG9nqxbA/6me3kE1LJZ9W/U1L4HdtsNpmf\nnwcoVMMG7dG4Dld8rsU1zo77n2uAL+OaXY8Gvgm8CTcj3Q5cAuwAvg7chOvPuNg/70W43o2z/c8V\nwKPAWh/nTn99PhdbJT1meSmvh4tYM2SPmaX6BQZqmOpXWMpneCnldKn6VaSw3YVrlP1V3EzzSuCv\ngHtxM8VZ4Bw6BecK4AO4fotLgW/4/VPAF4HXA/fjihy4r5vfCZyIK5Yb/HMCnO+fD+Bq4PYe8VVe\n1NpSOnBELOtjYGa9foGRGqb6FZbyGV4qOU3sD2T31LKg0Wi0Wq1WC9Z2XaqKo2oW4rAQQ6ulOPJC\nxEFafWHDJ3VIql/dQsUxbD4t5MNCDK1WWscoS9Qv/eX/wOr8TRIRiZvqV1jKZ3h1yGkKp9L8ANSW\nVE65iliU2FKAuRqm+hWW8hlezDnV/5VZkTqM6kUkTapfYSmf4aWcUw3MAsl+rb+tigOnVxxVsBCH\nhRhAceRZiUM6VL+6jSKOQfJpIR8WYoB6HaMamI1YyqN6EUmb6ldYymd4KeY0hR4Nc/0ZvcS8Hi5i\njXrMyqX6FZbyGV5MOVWPmREpjupFpB5Uv8JSPsNLKacamAVSZN25jAPHcj9AHWMAxZFnJQ7pUP3q\nVkYcRfJpIR8WYoB6HaMamJUspVG9iNSL6ldYymd4KeQ0hR4N8/0ZvcS0Hi5ijXrMqqX6FZbyGZ7l\nnKrHzKgURvUiUk+qX2Epn+HFnFMNzAIZZN15FAdOTP0AdYgBFEeelTikQ/WrWxVx9MqnhXxYiAHq\ndYxqYFaxmEf1IlJvql9h5fO5bt2FFUWSjhiP0RR6NKLrz+jF8nq4iDXqMbNF9Sss5TM8Szldqn6l\nUNiiL2ptlg4cEcs0MLNH9SusXmd2lNPhWDlG1fxfkhDrziFOucbcD5BiDKA48qzEIR2qX90sxNFq\n7eSMM07o2lfFMpyFXEC9jlENzIyJcT1cRARUv0LbvHlGOQ0shnymsBSQxDJAnpVTriIWaSnTNtWv\n8JTTsKrMp5YyIxXDqF5EpBfVr/CU07As51MDs0BGse48yIGTUj9ACjGA4sizEod0qH51sxpHFYMJ\nq7kIweoxqoGZcZZH9SIii1H9Ck85DctiPlPo0UiuP6MX9ReIdKjHLC6qX+Epp2GVmU/1mCXC4qhe\nRKQI1a/wlNOwLOVTA7NAylh3LnLgpNwPEGMMoDjyrMQhHapf3WKJo4zBRCy5CMHKMaqBWWQsjepF\nRPqh+hWechqWhXwW6dFYCdwB/GugBdwK3ASsAO4BjgNmgXOAef+Yy4EPAP8CXAI85PdPAV8EXgfc\nD1zq9y/3r7EW2AecC/zQ37YR+BN//Wp/v6zk+zN6UX+B1FkfPWbW6xfUsIapfoWnnIY1ynyG6DF7\nCfhPwG8AbwM+BLwZ2AxsBVYD3/LbAGtwhWkNsB64ORPALcAFwCp/We/3X4AraKuA64Fr/f4VwJXA\nSf5yFTBeIObkWRjVi0RA9csg1a/wlNOwqsxnkYHZc0DTX/8p8APgaOBM4Ha//3bgbH/9LOAuXEGc\nBZ4ETgaOBA4Gdvj73ZF5TPa57gNO9dfPwM1W5/1lK51iaEoV6/C9Dpw69QPEEAMojryS41D9KkD1\nq1uscYxiMBFrLkKo6hjtt8dsAjgR2A4cDuz1+/f6bYCjgD2Zx+zBFcL8/jm/H//zGX/9ZeBF4LBF\nnku8/IGzbt2FFUUiYt4Eql+mqH6FpzNnYVVxjPYzMPsV3GzwUmB/7raWv9TW9PR0Za/dfeAcbOKD\nWGU+LMUAiiOvojhUvxah+tUt9s9KyMFZ7LkIoexjdFnB+70GV9TuBP7S79sLHIFbKjgSeN7vn8M1\n3LYdg5spzvnr+f3txxwLPOtjOgTXszEHTGcesxJ4OB/czMwMExMTAIyPjzM5OfnKm9g+7Zj6dqu1\n0x8s7nfO2NgUrdZOM/FpW9vDbLevz87OMgDT9QtUwxqNLaxb9xGfjf2Mja2m1dplJr4Yt/U7Iez2\nMPlsNpvMz7vvFg1Yww4whuunuD63/zrgMn99M3CNv74G19PxWuB44Ck6zbPbcf0aY7hvNbX7LS7G\nNdYCbADu9tdXAE/jGmYPzVzPalnQaDSqDqHVarVasKoFa1+5VMVCPizE0GopjrwQcVD8DJf1+mWi\nhlk5NlS/uoX5rKwdKqcp5SKEEMfoUvXrVYvd6P0O8O+AdcD3/GW9L2SnAbuAd2UK2+PAvf7nA75o\ntYO4GLgN2I1rqn3Q7/8cridjN7CJzjekfgJ8EvgOrun2E3S+0i49NBq3dm1bWBYQqZDqV0RUv8JT\nz1lYZRyjKfxfc34AKln6mzaSMv1fmWlT/QpPOQ1rmHzq/8qsKc2SRCRWql/hKadhjTKfGpgFkm1S\nrlI2jio/iBbyYSEGUBx5VuKQDivviepXt9BxDJLTVHMxqDKOUQ3MEqdZkojESvUrPOU0rFHkM4Ue\nDfVnFKD+AkmJeszqRfUrPOU0rH7yqR4zATRLEpF4qX6Fp5yGFTKfGpgFYnH9O6/MD6KFfFiIARRH\nnpU4pMPKe6L61W3UcRTJaV1yUVQZx6gGZjWjWZKIxEr1KzzlNKwQ+UyhR0P9GQNQf4HETD1m9ab6\nFZ5yGtZi+VSPmfSkWZKIxEr1KzzlNKxh8qmBWSAxrH/njfKDaCEfFmIAxZFnJQ7psPKeqH51KzuO\nXjmtay4WUsYxqoFZzWmWJCKxUv0KL5/TdesurCiSNAxyjKbQo6H+jADUXyAxUY+ZZKl+haechtWd\nz+9C6j1mmiUNTzNPEYmV6ld4GoiF1U8+kxiYQfUfxBjXv/NCFjcL+bAQAyiOPCtxSIeV90T1q1vV\ncbic7q80hraqc9EW8hhdSDIDM6h+cJYCzTxFJFaqX+E1Grfq7FlARXKZQo9GC9Z279BBNDT1F4hl\n6jGTxah+iWW1+DtmmiWFp5yKSKxUvyRmSQzMoPoPYgrr33nD5NRCPizEAIojz0oc0mHlPVH96mYh\nDgsxQL3iSGZgBtUPzlKknIpIrFS/JEYp9Ggc0J+h/oLwlFOxRD1m0g/VL7GkFj1meZolhaecikis\nVL8kJkkOzKD8D2Id1r/7yamFfFiIARRHnpU4pMPKe6L61c1CHBZigHrFkezADDRLGgXlVERipfol\nMUihR2PJ/gz1F4SnnEqV1GMmw1D9kirVsscsT7Ok8JRTEYmV6pdYttTA7HXAdqAJPA582u9fAWwF\ndgEPAeOZx1wO7AaeAE7P7J8CHvO33ZjZvxy4x+9/BDguc9tG/xq7gPMK/pt6GvUHsU7r322L5dRC\nPizEAIojr+Q4kqlho1THY8N6/QIbcViIAeoVx1IDs38C1gGTwG/6678LbMYVtdXAt/w2wBrgXP9z\nPXAzndN1twAXAKv8Zb3ffwGwz++7HrjW718BXAmc5C9X0V08+6ZZUnjKqRiXVA2TsFS/xKJ+ejQO\nAr4NzAD3Ae8E9gJHANuAX8fNNH9JpzA9CHwc+CHwMPBmv38DMA180N/nKtysdhnwI+CNwPuBdwAX\n+cd8xr/O3bm4+u7PUH9BeMqplGnAHrNkapiEpfolZQrRY/Yq3DLAXqAB/B1wuN/G/zzcXz8K2JN5\n7B7g6B775/x+/M9n/PWXgReBwxZ5rqFplhSeciqGJVfDJCzVL7FkWYH7/BK3DHAI8A3cUkBWy18q\nMzMzw8TEBADj4+NMTk4yPT0NdNaD89ut1k7/4dsPuA9iq7Vzwfsvtd3eN+jjQ23fcMMNhf79o9h2\nOV39Sj7GxqZoNLZUlo/8e1P267e3m80mmzZtquz129sx56N9fXZ2lgEkWcOqfk9SO0a7fyf8fOjf\nCbHno71d5e+U7HZ7X5XvBwyWj2azyfz8PMCgNWxR/xn4I1xT7BF+35F+G1yfxubM/R8ETvb3/UFm\n//tx/Rrt+7zNX18GvOCvb8Cd+m/7c1zvR15rGLC26zKoRqMxVByhWIjD5XLV0DkdloVctFqKIy9E\nHAw+kEquhoWQ0rExLCv1q9WykQ8LMbRaacXBEvVrqR6NX8Wdmp8HXo+bbX4COAPX7HotroiN+59r\ngC/jGl2PBr4JvMkHsR24BNgBfB24yRe0i4ETcH0YG4Cz/c8VwKPAWh/nTn99vkdRW+KfsTj1F4Sn\nnMoo9dFjVosaJmGpfskoLVW/lipsJwC343o0XgXcCfxXXMG5FzgWmAXOoVNsrgA+gCuGl+IKIbiv\nmn8RVxzvxxU4cF81vxM4EVcoN/jnBDjfPx/A1T6WvCBFTR/E8JRTGZU+Bma1qWESluqXjEpifyC7\np6FPK3ZOLw6+rJnSadYQ2nGEWioeJoaqKY5uFS9lWjR8UoeU0rERgoX6lY2jShZiaLXSioMl6lct\n/vJ/UfpmTnjKqYjESvVLqpDCqTQ/AA1Hp7DDU04lpMSWAoLXMAlL9UtC0v+VOQDNksJTTkUkVqpf\nUiYNzBbQ7wcx+7dWqmQ5jrKLm+VcVEFxyEKsvCeW46hicGYhHxZigHrFoYHZIjRLCk85FZFYqX5J\nGVLo0Rh5f4b6C8JTTmUY6jGTKql+yTDUYxaAZknhKaciEivVLxklDcwKWuqDWKf17yKKxDHq4hZT\nLsqgOGQhVt6TmOIoY3BmIR8WYoB6xaGBWR80SwpPORWRWKl+ySik0KNRen+G+gvCU06lH+oxE0tU\nv6Qf6jEbAc2SwlNORSRWql8SkgZmAzrwg7i6oki6xbwOH7q4xZyLUVAcshAr70nMcYxicGYhHxZi\ngHrFoYHZEDRLCk85FZFYqX5JCCn0aFTen6H+gvCUU1mMeszEMtUvWYx6zEqgWVJ4yqmIxEr1S4ah\ngVkgjcaWru2qPogprcMPW9xSykUIikMWYuU9SSmOEIMzC/mwEAPUKw4NzALSLCk85VREYqX6JYNI\noUfDXH+G+gvCU04lSz1mEhPVL8lSj1kFNEsKTzkVkVipfkk/NDALJL/uXNUHMeV1+H5zmnIuBqE4\nZCFW3pOU4xjkd4KFfFiIAeoVhwZmI6RZUnjKqYjESvVLikihR8N8f4b6C8JTTutNPWYSM9WvelOP\nmQGaJYWnnIpIrFS/ZDEamAWy1LpzWR/EOq3DL5XTOuWiCMUhC7HyntQpjiK/Eyzkw0IMUK84ig7M\nXg18D/ia314BbAV2AQ8B45n7Xg7sBp4ATs/snwIe87fdmNm/HLjH738EOC5z20b/GruA8wrGapZm\nSeEpp1KQapiYo/olvRTt0fhDXFE6GDgTuA74sf95GXAosBlYA3wZeCtwNPBNYBXQAnYAH/Y/7wdu\nAh4ELgbe4n+eC/w+sAFXOL/jXxdgp78+n4stuv4M9ReEp5zWywA9ZqphYpbqV72E6DE7Bvg94LbM\nE50J3O6v3w6c7a+fBdwFvATMAk8CJwNH4griDn+/OzKPyT7XfcCp/voZuJnsvL9sBdYXiNc8zZLC\nU05lEaphYprql2QVGZhdD/wx8MvMvsOBvf76Xr8NcBSwJ3O/PbhZZ37/nN+P//mMv/4y8CJw2CLP\nZVK/686j+iDWaR0+78Ccri49hl7q/J70UkEcqmFLqPGx0ZON+jVlIh8WYoB6xbHUwOy9wPO43oyF\nTru1/EX6pFlSeMqp5KiGSTTy9WvdugsrikSqtGyJ29+OO03/e8DrgH8F3ImbYR4BPIc7xf+8v/8c\nsDLz+GNws8Q5fz2/v/2YY4FnfTyHAPv8/unMY1YCD/cKcmZmhomJCQDGx8eZnJxketo9tD26tbrd\naGzxH76DAXeWp9G4deDna++z8u+rYtvl9CO4nO5nbGw1rdauSuNrqzI/09PTJt6frH7uv23bNmZn\nZ+mTapiO0YG228p+/e7fCQcP/Tth2O32vqrfDyvb7X39PL7ZbDI/71pLi9Swfppn3wn8EfA+XMPs\nPuBaXMPsON2NsyfRaZx9E242uh24BNej8XW6G2dPAC7CNcyeTadx9lFgrY9zp7+eZOOsmj/DU07T\nNeAfmFUNk2iofqUr9B+YbVePa4DTcF8Bf5ffBngcuNf/fABXsNqPuRjXfLsb11D7oN//OVw/xm5g\nE644AvwE+CTuW007gE9wYEEzIz/L6leoJbhh4wjFQhyNxpau7aqWNS3kAhSHpxrWg46NbhbicL8T\n9r+yrfq1reoQgHLiWGopM+vb/gKu4Lx7gfv9qb/k7cTNKvP+GThngef6gr/UQqu1s+vDNzY2pVnS\nkJRTyVANk6g0Grf6tgxH9aseUvi/5pJbBtAp7PCU07To/8qUOlH9Sov+r8wI6ZuF4SmnIhIr1a96\n0cAskNDrzoN+EOu0Dt9vDFUVNwu5AMUhC7PyniiObtk4VL+2VR0CUE4cGpgZpllSeMqpiMRK9ase\nUujRSL4/Q/0F4SmncVOPmdSZ6lfc1GOWAM2SwlNORSRWql9p08AskFGvOxf9INZpHX7YGMoqbhZy\nAYpDFmblPVEc3RaLQ/WrGuoxky6aJYWnnIpIrFS/0pRCj0bt+jPUXxCechoX9ZiJdKh+xUU9ZgnS\nLCk85VREYqX6lRYNzAIpe/17oQ9indbhQ8cwquJmIRegOGRhVt4TxdGtnzhUv8qhHjNZlGZJ4Smn\nIhIr1a80pNCjUfv+DPUXhKec2qYeM5GFqX7Zph6zGtAsKTzlVERipfoVNw3MAql6/bvzQdwPVP9B\nrDofIWIIVdws5AIUhyzMynuiOLoNE4fq12iox0z6ollSeMqpiMRK9StOKfRoqD8jR/0F4SmntqjH\nTKQ41S9b1GNWQ5olhaecikisVL/iooFZINbWv6v+IFrIR+gYBs2phVyA4pCFWXlPFEe3kHGofoWh\nHjMZStWDsxQppyISK9WvOKTQo6H+jCWovyA85bRa6jETGZzqV7XUYyaaJY2AcioisVL9sk0Ds0Cs\nr3+X/UG0kI9Rx1A0pxZyAYpDFmblPVEc3UYZh+rXYNRjJkFplhSecioisVL9simFHg31Z/RJ/QXh\nKaflUo+ZSDiqX+UK1WM2C/wt8D1gh9+3AtgK7AIeAsYz978c2A08AZye2T8FPOZvuzGzfzlwj9//\nCHBc5raN/jV2AecVjFcWoVlSeMqpabOofoksSPXLlqIDsxYwDZwInOT3bcYVttXAt/w2wBrgXP9z\nPXAznZHhLcAFwCp/We/3XwDs8/uuB671+1cAV/rXPAm4iu4CakZs69+j/iBayEfZMSyUUwu5gFrH\nofq1hBrKUpirAAAPcUlEQVQfGz3VMQ7Vr2Ks9ZjlT7udCdzur98OnO2vnwXcBbyEm6k+CZwMHAkc\nTGfGekfmMdnnug841V8/AzebnfeXrXSKoQxJs6TwlFOzVL9ElqD6ZUPRHo2ngReBfwH+HPgs8I/A\noZnn+Ynf/jPc6fwv+dtuAx7AFblrgNP8/lOAjwLvwy0PnAE8629rF8MZ4HXAp/z+jwG/ALZkYlN/\nxpDUXxCecjpaffaYWa5foBomxqh+jVaoHrPfwS0DvAf4EK4oZbX8pRIa1Q9Hs6TwlFNTTNcvEWtU\nv6q1rOD9fuR/vgB8FdcvsRc4AngOd5r/eX+fOWBl5rHHAHv8/mN67G8/5ljcjHMZcAiuZ2MO1xvS\nthJ4+MDwZhkbO5KrrvoPjI+PMzk5yfS0e1h7PXjU2+19Zb3eQts33HDDQP/+Vmun//DtJyvmfORj\nKfv1XU5XAz8HDmdsbIpGY0tt89HebjabbNq0qa/Ht6/Pzs4yAOP1C2ZmZpiYmACopIYN8p7oGB39\ndpX56PxO2AscxNjYFK3Wzsry0d5X5fsBg/2ObTabzM/PAwxaww5wEK63AuANwN/gvql0HXCZ378Z\nd5ofXNNsE3gtcDzwFJ1Tdttxp/jHgPvp9FtcjGusBdgA3O2vr8AtQ4zjlhna17NasPaVS1UajUZl\nr52lOGzF0Gq1WrBKx2hGiDgofobLev0Ct5RZqZSOjRAUR4fqV7cy6leRHo3jcbNMcLPBLwGf9kXn\nXtxMcRY4B9fgCnAF8AHgZeBS4Bt+/xTwReD1uMJ2id+/HLgTt9ywD1fcZv1t5/vnA7iaTpNtWwvW\ndu/QergYo56NsProMbNev0A9ZmKc6ldYS9WvFP5A4wEDM9CBI/aouIWjPzArUi7Vr3Bq8Z+Y9zpA\nym5WzK6DV0lx2IoBuns2snSMihVW3hPF0c1CHKpf3cqII4mBGdgYnIksperiJiIyKNWvcqSwFNC1\nDNDrQNEpV7FGywLD0VKmSHVUv4ZTi6XMLJ05kxho5ikisVL9Gq3kBmZQzeCsTuvfRViIw0IMsHAc\nZRc36/mQ6lh5TxRHNwtxqH51U4/ZEHTmTGKgmaeIxEr1azRS6NFYtD9DPWcSA/Vs9Ec9ZiJ2qH71\np3Y9Znk6cyYx0MxTRGKl+hVW8gMzKGdwVqf17yIsxGEhBigex6iLW2z5kPJYeU8URzcLcah+dVOP\nWUA6cyYx0MxTRGKl+hVGCj0affVnqOdMYqCejcWpx0zELtWvxdW+xyxPZ84kBpp5ikisVL+GU7uB\nGYxmcFan9e8iLMRhIQYYPI7QxS32fMjoWHlPFEc3C3GofnVTj9kI6cyZxEAzTxGJlerXYFLo0Riq\nP0M9ZxID9Wx0U4+ZSDxUv7qpx2wJOnMmMdDMU0RipfrVn9oPzCDM4KxO699FWIjDQgwQLo5hi1tq\n+ZBwrLwniqObhThUv7qpx6xEOnMmMdDMU0RipfpVTAo9GkH7M9RzJjGoe8+GesxE4qX6pR6zvujM\nmcRAM08RiZXq1+I0MOthkMFZnda/i7AQh4UYYHRx9FvcUs+HDM7Ke6I4ulmIQ/Wrm3rMKqQzZxID\nzTxFJFaqX72l0KMx0v4M9ZxJDOrWs6EeM5F0qH510xmzJejMmcRAM08RiZXqVzcNzAooMjir0/p3\nERbisBADlBfHUsWtbvmQ4qy8J4qjm4U4VL+6WeoxGwe+AvwAeBw4GVgBbAV2AQ/5+7RdDuwGngBO\nz+yfAh7zt92Y2b8cuMfvfwQ4LnPbRv8au4DzCsYbnM6cSQw08+yp9vVLJAaqX07RHo3bgW8DnweW\nAW8A/gT4MXAdcBlwKLAZWAN8GXgrcDTwTWAV0AJ2AB/2P+8HbgIeBC4G3uJ/ngv8PrABVzy/gyuI\nADv99flMbKX2Z6jnTGKQes9Gnz1mlusXqMdMpEvd61eRM2aHAKfgihrAy8CLwJm4gof/eba/fhZw\nF/ASMAs8iZuhHgkcjCtqAHdkHpN9rvuAU/31M3Cz2Xl/2QqsLxDzyOjMmcRAM89XqH6JRKbu9avI\nwOx44AXgC8B3gc/iZpyHA3v9ffb6bYCjgD2Zx+/BzTzz++f8fvzPZ/z1duE8bJHnqlTvwdnqCiI5\nUJ3W4WOIAaqL48DiVstjVPWrgLp/VvIUR/Ux1Ll+FRmYLQPWAjf7nz/DnfLPavlLbejMmcSg7jNP\nVL9EolXX+rWswH32+Mt3/PZXcM2xzwFH+J9HAs/72+eAlZnHH+MfP+ev5/e3H3Ms8KyP6RBgn98/\nnXnMSuDhfIAzMzNMTEwAMD4+zuTkJNPT7mHt0e0otlutnZlR/MGAG9U3GreW8vq9ttv7qnp9S9vT\n09Nm4mmr4vUbjS2sW/cR3DG6n7Gx1bRauyqLJ6uf+2/bto3Z2Vn6ZL5+QXU1LLvdps+s226rez7a\n+6p6fVe/LiTm37HNZpP5eddaWqSGFW2e/WvgD3DfLPo4cJDfvw+4FjcDHae7efYkOs2zb8LNSLcD\nl+D6NL5Od/PsCcBFuKbZs+k0zz6Km+mO4Zpn11Jh838v+kKAxCClhto+m/8t1y8wUMNErKtT/Sr6\n5zL+I/Al4PvAbwKfAq4BTsMVu3f5bXBfR7/X/3wAV7TaVedi4Dbc18qfxBU1gM/hejJ2A5voLDX8\nBPgkbra7A/gEBxa1yrkDZH/XvqpOueZne1WxEIeFGMBOHI3Glq7tGh2jql9LsHKMKo5uFuKwEAPU\nq34VWcoEV9De2mP/uxe4/5/6S95O3Mwy75+BcxZ4ri/4i2mNxq1+uahjbGwq6lG9pMctv3cKWk2O\nUdUvkQTUpX6l8H/NmVoG0LKmxCD2ZQH9X5ki9ZV6/Sq6lCkF6duaEoO6fttJROKXev3SwCyQ7Lpz\nlYMzK/0AFuKwEAPYjaOq4mYlH9Jh5T1RHN0sxGEhBqhX/dLAbER05kxikPrMU0TSlWr9SqFHw3R/\nhnrOJAax9Wyox0xE2lKrXzpjNmI6cyYxSHXmKSLpS61+aWAWyGLrzmUOzqz2A9Q1BognjrKKm5V8\nSIeV90RxdLMQh4UYoF71SwOzkujMmcQgtZmniNRHKvUrhR6NqPoz1HMmMbDes6EeMxFZSOz1S2fM\nSqYzZxKDVGaeIlI/sdcvDcwC6WfdeZSDs1j6AeoSA8Qbx6iKm5V8SIeV90RxdLMQh4UYoF71SwOz\niujMmcQg9pmniNRXrPUrhR6NqPsz1HMmMbDWs6EeMxEpKrb6pTNmFdOZM4lBrDNPEZHY6pcGZoEM\ns+4ccnAWaz9AqjFAOnGEKm5W8iEdVt4TxdHNQhwWYoB61S8NzIzQmTOJQWwzTxGRtljqVwo9Gkn1\nZ6jnTGJQdc+GesxEZFDW65fOmBmjM2cSg1hmniIiedbrlwZmgYRcdx5mcJZKP0AqMUC6cQxa3Kzk\nQzqsvCeKo5uFOCzEAPWqXxqYGaUzZxID6zNPEZGFWK1fKfRoJN2foZ4ziUHZPRvqMRORUKzVL50x\nM05nziQGVmeeIiJLsVa/NDALZJTrzv0MzlLtB4g1BqhPHEWLm5V8SIeV90RxdLMQh4UYoF71SwOz\nSOjMmcTA2sxTRKQoK/UrhR6NWvVnqOdMYjDqng31mInIqFRdv4qcMfs14HuZy4vAJcAKYCuwC3gI\nGM885nJgN/AEcHpm/xTwmL/txsz+5cA9fv8jwHGZ2zb619gFnFcg3qTpzJnEwMrME9UvEelT1fWr\nyMDs74ET/WUK+DnwVWAzrrCtBr7ltwHWAOf6n+uBm+mMDG8BLgBW+ct6v/8CYJ/fdz1wrd+/ArgS\nOMlfrqK7gJpR5jr8YoOzuvQDxBID1DeOhYpbyXGofhVQ12N0IYrDVgxQr/rVb4/Zu4EngWeAM4Hb\n/f7bgbP99bOAu4CXgFl//5OBI4GDgR3+fndkHpN9rvuAU/31M3Cz2Xl/2UqnGNaazpxJDKqeeeao\nfolIYVXVr357ND4PPIqbRf4jcGjmeX7it/8Mdzr/S/6224AHcEXuGuA0v/8U4KPA+3DLA2cAz/rb\n2sVwBngd8Cm//2PAL4AtmZhq3Z+hnjOJQeiejQF7zCzWL6h5DROxruz61c8Zs9fiitBf9Lit5S9S\nMp05kxgYOHOm+iUiAym7fi3r477vAXYCL/jtvcARwHO40/zP+/1zwMrM444B9vj9x/TY337MsbgZ\n5zLgEFzPxhwwnXnMSuDhfGAzMzNMTEwAMD4+zuTkJNPT7mHt9eBRb7f3lfV62e1GYwvr1n3ER7EX\nOIixsSlarZ2VxFN1PvKvXdXrt7ebzSabNm2q7PXb21Xno9XaydjYalyb1+GMjU3RaGwp9Pj29dnZ\nWQZktn5B9TVMx6jysdD2DTfcUMnv1F41IJuTsl/f1a8pBvkd22w2mZ+fBximhvV0N+4bRm3XAZf5\n65txp/nBNc02cTPU44Gn6Jyy2447xT8G3E+n3+JiXGMtwAb/WuCaZ5/GNcwemrme1bKg0WhUHUIL\n1rZglf/pLlWxkA8LMbRaiiMvxDFK/2e4rNYvEzXMyrGhOLpZiMNCDK2WnTjKqF9FezTeAPzQF6r9\nmaJzL26mOAucg2twBbgC+ADwMnAp8A2/fwr4IvB6XGG7xO9fDtyJ++bUPlxxm/W3ne+fD+BqOk22\nbf7fKaCeM4nDsD0bffaYWa5foBomEpVR168U/kCjilqOBmcSg2GKm/7ArIhUaZT1S/8lUyDZdfAq\nbdu2zcQXAizkw0IMoDjysj0bWfrSSnWsHRtVUxy2YgB7cYyyfmlgligLgzORpWhwJiKxGlX9SmEp\nQMsAi9CypsSg32UBLWWKiBWh65fOmCVOZ84kBjpzJiKxCl2/NDALxNr6d1YVgzML+bAQAyiOvIXi\n0OCsOtaPjbIpDlsxgP04QtYvDcxqQmfOJAYanIlIrELVrxR6NNSf0Qf1nEkMlurZUI+ZiFg1bP3S\nGbOa0ZkziYHOnIlIrIatXxqYBWJ9/TurjMGZhXxYiAEUR17RODQ4K09sx8aoKQ5bMUB8cQxTvzQw\nC6TZbFYdAlA8jlEPzizkw0IMoDjy+olDg7NyxHhsjJLisBUDxBnHoPVLA7NA2v9zfNX6iWOUgzML\n+bAQAyiOvH7j0OBs9GI9NkZFcdiKAeKNY5D6pYFZzannTGKgL6iISKz6rV8amAUyOztbdQjAYHGM\n4peehXxYiAEUR96gcWhwNjqxHxuhKQ5bMUD8cfT1n5wP9Aq2NIHfqjoIESnTG4CfpVC/QDVMpGaS\nql8iIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiFfn/wDFCxka64TcAAAAASUVORK5CYII=\n", + "text": [ + "" + ] + } + ], + "prompt_number": 11 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def ricker(fpeak, t, tlag):\n", + " \"\"\"\n", + " Generating Ricker Wavelet\n", + " \n", + " .. math ::\n", + " \n", + " \n", + " \"\"\"\n", + " return (1-2*np.pi**2*fpeak**2*(t-tlag)**2)*np.exp(-np.pi**2*fpeak**2*(t-tlag)**2)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 12 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "wave = ricker(600, time, 0.0025)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 13 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fig, ax = plt.subplots(1,1, figsize = (7, 5))\n", + "ax.plot(time, wave, '.-')\n", + "ax.set_xlim(0, 0.1)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 14, + "text": [ + "(0, 0.1)" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAE4CAYAAAAtjzZWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG55JREFUeJzt3X2UXGWB5/FvpTuvRNK8hCBJPJEXF/TIi/YGzMhaAzjD\nRg2wo3BwXMVZCXNmUHdVBNcz0npgx4zLujKMbOIwnIyzDA4vMnAkYqPWqrsShA3vL4EgmoQlhIRO\nxoSQ7nTtH09VurpS1V333k73faq/n3Pq3Hurbt17c0+nf/28XpAkSZIkSZIkSZIkSZIkSZIkScrs\n74AtwOMj7HM98BzwKHDaeFyUJElj5UxCeDULuqXAvZX104EHxuOiJEkaS4toHnT/A7ioZvsZYN7B\nviBJkgCmjMM55gMba7Y3AQvG4bySJI1L0AEU6rbL43ReSdIk1zkO59gMLKzZXlB5b5jjjjuuvGHD\nhnG4HElSJDYAx2c9yHiU6O4GPl5ZPwPoI/TSHGbDhg2Uy+URX6EgWKa7u8xrr42872R5XX311RN+\nDTG+vG/eN+9b/l/AcWMRQmNRovtH4H3AkYS2uKuBqZXPVhJ6XC4Fngd2AZ/McrKzzoI77oCurixH\nkSRNFmMRdBe3sM/lWU9SrrTqXXutISdJat14dUbJrL9/+FJBsVic6EuIkvctHe9bOt63iVXfG3Ii\nlcvl5p0xd+6EOXPg/vvh7LPH8aokSROiUCjAGORUNCW6N94YvpQkqRXRBN2ePWG5d+/EXockKS7R\nBJ0lOklSGtEE3VVXheW110Jf38ReiyQpHtEE3a9/HZaPPw7Ll0/stUiS4hFN0E2fHpbTpsFrr1mq\nkyS1Jpqg+9rXwnLv3jDEwFKdJKkV0QTdrFlD693dsGrVxF2LJCke0QTdvn1hefLJ0NvrNGCSpNZE\nE3SDg2H5sY8ZcpKk1kUTdNUSXXUpSVIrogm6aomuupQkqRXRBJ0lOklSGtEEnSU6SVIa0QSdJTpJ\nUhrRBJ0lOklSGtEEXbUkZ9BJkpKIJuiqAWfVpSQpiWiCzhKdJCmN6ILOEp0kKYlogs7OKJKkNKIJ\nOkt0kqQ0ogk6S3SSpDSiCTpLdJKkNMYi6M4FngGeA65s8PmRwA+BR4AngEvSnMQSnSQpjaxB1wHc\nQAi7twMXAyfV7XM5sA44FSgC1wGdSU+0bx90dFiikyQlkzXoFgPPAy8C/cCtwHl1+/w/4NDK+qHA\nNmAg6YkGB2HqVEt0kqRkEpes6swHNtZsbwJOr9vnO8BPgJeANwEXpjnRvn3Q2WmJTpKUTNYSXbmF\nff4zoX3uGEL15d8QAi+Rffss0UmSkstaotsMLKzZXkgo1dVaAlxbWd8A/Br4V8BD9Qfr6enZv14s\nFikWi/u3rbqUpPZWKpUolUpjftxCxu93As8CZxOqJh8kdEh5umaf/wbsAL4KzAMeBk4Gttcdq1wu\nNy8gfuMbcP31sGQJfO97Ga9akpR7hUIBsudU5hLdAKFX5X2EHpg3EULussrnK4H/AtwMPEqoKv0i\nB4bcqCzRSZLSyBp0AGsqr1ora9ZfBT6U9STVNjo7o0iSkohmZhRLdJKkNDLXfY6hEdvo3v1uePZZ\nmDUL1q+Hrq5xvDJJ0rgbqza6aEp0r74Ku3bB1q2wfPlEX40kKRbRBN3UqWE5ZQq89hr09U3s9UiS\n4hBN0J1/fpgZZXAQ7r/fUp0kqTXRBN20afCmynwq3d2watXEXo8kKQ7RBF25DO98J8ydC729dkaR\nJLUmmqAbHISZM+Ed7zDkJEmtiyroOjocRydJSiaaoCuXhzqjSJLUqmiCzhKdJCmN6IJuhMlTJEk6\nQDRBZ9WlJCmNaIJucNCgkyQlF1XQ2UYnSUoqmqArl22jkyQlF03QWXUpSUojqqCz6lKSlFQ0QVet\nujToJElJRBN0jqOTJKURVdDZRidJSiqaoHPAuCQpjWiCzs4okqQ0ogk6O6NIktKIJuiqbXR2RpEk\nJRFV0FmikyQlNRZBdy7wDPAccGWTfYrAOuAJoJTmJHZGkSSl0Znx+x3ADcA5wGbgV8DdwNM1+3QB\nfwP8IbAJODLNiSzRSZLSyFqiWww8D7wI9AO3AufV7fNR4A5CyAG8muZEDhiXJKWRNejmAxtrtjdV\n3qt1AnA48FPgIeDfpzmRVZeSpDSyVl22Ur6aCrwLOBuYBfwSeIDQptcyZ0aRJKWRNeg2Awtrthcy\nVEVZtZFQXfl65fUz4BQaBF1PT8/+9WKxSLFY3L9tG50ktbdSqUSpVBrz4xYyfr8TeJZQWnsJeBC4\nmOGdUU4kdFj5Q2A6sBa4CHiq7ljl8ggNcMuWwXnnwVVXwdatGa9akpR7hUIBsudU5hLdAHA5cB+h\nB+ZNhJC7rPL5SsLQgx8CjwGDwHc4MORGZdWlJCmNrEEHsKbyqrWybvu/Vl6pWXUpSUojmplR7HUp\nSUojmqCzRCdJSiO6oHPAuCQpiWiCzqpLSVIa0QSdvS4lSWlEFXS20UmSkoom6KpPGLeNTpKURDRB\nZ9WlJCmNqILOqktJUlLRBF216rK6LklSKzJPljmGRpzU+aijYMECWLcOXn0VjjhiHK9MkjTuxmpS\n52hKdK+/HkIO4E//dGKvRZIUj2hKdHPmwM6dYX3LllDCkyS1r0lXolu0CM45B6ZMgQ9/GJYuhb6+\nib4qSVLeRRN0U6bAihWhI8rPfw5r1sDy5RN9VZKkvIsm6MrlEHaFSiG2uxtWrZrYa5Ik5V80QTc4\nGEJu5ky44ALo7YWurom+KklS3o3FE8bHRbVE19EBN98cOqdIkjSaqEp01apLZ0eRJLUqqqArFELY\nOTOKJKlV0QRdtepyyhRLdJKk1kUTdNWqS4NOkpREVEFXrbo06CRJrYom6GrH0Rl0kqRWRRN0tVWX\ndkaRJLUqqqCz6lKSlNRYBN25wDPAc8CVI+z3r4EB4N+lOYm9LiVJaWQNug7gBkLYvR24GDipyX4r\ngB+S8pELDhiXJKWRNegWA88DLwL9wK3AeQ32+zRwO7A17YkcMC5JSiNr0M0HNtZsb6q8V7/PecCN\nle1UMWXVpSQpjayTOrcSWv8duKqyb4ERqi57enr2rxeLRYrF4v5tB4xLUnsrlUqUSqUxP27WR5Sf\nAfQQ2ugAvgQMEtrjql6oOc+RwG7gUuDuumOVyyPUSc6dC08+CWeeCf/8z3DiiRmvXJKUa4XwANKs\nOZW5RPcQcAKwCHgJuIjQIaXWsTXrNwP3cGDIjaq26tI2OklSq7IG3QBwOXAfoWflTcDTwGWVz1dm\nPP5+Vl1KktIYiwevrqm8ajULuE+mPYkDxiVJaUQzM4q9LiVJaUQTdA4YlySlEVXQOWBckpRUNEFn\n1aUkKY1ogs5el5KkNKIKukLBNjpJUjLRBJ0DxiVJaUQTdFZdSpLSiCroHDAuSUoqiqCrVlXaRidJ\nSiq6oLONTpKURDRBN6VypVZdSpKSiCLoqh1RwKCTJCUTTdAVKo/eM+gkSUlEEXS1VZd2RpEkJRFF\n0NVXXdoZRZLUqmiCzqpLSVIaUQSdvS4lSWlFEXS1VZe20UmSkogm6GqrLm2jkyS1Koqgs+pSkpRW\nFEHngHFJUlrRBF216tI2OklSElEEXX3VpW10kqRWRRF0Vl1KktKKIui+8AXYtg2WLoWBAYNOktS6\nsQi6c4FngOeAKxt8/sfAo8BjwP8GTk56ghdegP5+WLMGHnzQoJMkta4z4/c7gBuAc4DNwK+Au4Gn\na/Z5Afg3wA5CKK4CzkhykhkzwrK7G447zqCTJLWukPH77wGuJgQYwFWV5deb7H8Y8DiwoMFn5XKT\nXiaPPQbveQ9s3gyLF4d2umOPhVtuga6uDFcvScqtQuhunzWnMlddzgc21mxvqrzXzH8A7k16ktmz\nYd68EGo7d8Kzz4ZqzOXLkx5JkjTZZK26TNLR//eBPwF+r9kOPT09+9eLxSLFYhEY3uty6tSw7O6G\nVasSXaskKcdKpRKlUmnMj5u1SHgG0MNQ1eWXgEFgRd1+JwN3VvZ7vsmxmlZdrl8PH/xgWH784/DU\nU3D//VZbSlI7y0vV5UPACcAiYBpwEaEzSq23EELuYzQPuRHVluhmzoRLLzXkJEmtyVp1OQBcDtxH\n6IF5E6HH5WWVz1cCXyF0Qrmx8l4/sDjJSXzwqiQpraxBB7Cm8qq1smb9U5VXaj69QJKUVhQzo/jg\nVUlSWtEFnSU6SVISUQadTy+QJLUqyqCzRCdJapVBJ0lqa9EFnZ1RJElJRBd0ttFJkpKIJugcMC5J\nSiOKoHPAuCQprSiCzjY6SVJa0QWdbXSSpCSiDDpLdJKkVhl0kqS2ZtBJktpadEFnZxRJUhLRBZ2d\nUSRJSUQTdA4YlySlEUXQOWBckpRWFEFnG50kKa3ogs42OklSElEGnSU6SVKrDDpJUluLLuhso5Mk\nJRFd0NlGJ0lKIpqgcxydJCmNaILONjpJUhpjEXTnAs8AzwFXNtnn+srnjwKnJT2BA8YlSWllDboO\n4AZC2L0duBg4qW6fpcDxwAnAcuDGpCexM4okKa2sQbcYeB54EegHbgXOq9tnGbC6sr4W6ALmJTmJ\nnVEkSWl1Zvz+fGBjzfYm4PQW9lkAbKk/WLXDCUBnJzz8MJx88uhtdCeeCM8+m/afIEkHX+3vNI2v\nrEHXatmqULfd5Hs9+9cGBoqcckqRF18cOeguvdSQk5R/AwNw+unw+usTfSX5VSqVKJVKY37crEG3\nGVhYs72QUGIbaZ8Flfca6DngnSVL4C/+onkb3fe/n/CKJWkCFAqwdu1EX0W+FYtFisXi/u2vfvWr\nY3LcrEH3EKGTySLgJeAiQoeUWncDlxPa784A+mhQbQmhaA/hL5+qU08duY1u9+4Dj9PREX6oymXY\nty9sw9B6/Wfum3zfvF9fbPvm/fpi2zdv19fZCQ8+aLXlRMkadAOEELuP0APzJuBp4LLK5yuBewk9\nL58HdgGfbHaw/v6wnDsXXn0VZs+Gb38b7rln9AHj1n9LkhrJGnQAayqvWivrti9PcsDjjw9B97vf\nwRVXwHvf27yNrrYUeM018E//lOjaJUltLpczoxx2WFiedhqsWjV8wHh9G121mrO7O+wrSVKtXAbd\nLbeEklq5DB/9KOza1biNrlyGPXvgiCOgq2virleSlF/13f4nUrlc08ukqwt27Ajrhx8OM2eG9reL\nLw49Le+8M3TTPeSQoeD7yEesupSkdlEInTMy59RYtNEdFNOnh2V3N2zdCr/5DWzeHMJv7tzw2c6d\noeTX32/VpSSpsVxWXQKcfXYIr97e4aE3fTr8/OewdGkIvvnzQ0mut9fqS0nSgXJbops3LwRbVxdc\ndFGokuztDT0wt2+HNWtCSa6ry+pKSVJzuS3R/exn8K1vhZLb4CD80R+FUJsxI3ze3R0GYr7wQtin\nr29ir1eSlE+5Dbq+Pvjtb0PJ7a67hgaMX3llKO319sKmTaGdbs0aWL58Yq9XkpRPua26nDkzLLu7\n4ZxzhoYXHHoonHJKKN1Vp9mxI4okqZnclug+8xl4y1uGOqM0GjD+iU/AscfaEUWS1Fxug+6oo8KE\nzl1dzSd13rcPLrjAkJMkNZfboDvkkDAjCjR/Ht3OnaEqU5KkZnLbRtdK0N11V3hMzwMPhGnDLNlJ\nkurltkR33XXw2GNh6MDu3Y3b6LZuDcML7HUpSWomt0G3cWMIuDVrwqtRG131PXtdSpKayW3V5axZ\nYdndDaef3rjqcvFi2LIlBKHVlpKkRnJbolu5Mgwr6O2FadMaP2G8XA4DyA05SVIzuQ26o44K0311\ndQ1/8Gpt0O3dOzThsyRJjeQ26KZPDw9VheG9Lms7o7zxRijtSZLUTK6Dbu/eUJprNmDcEp0kaTS5\nDbqOjvDq728+js4SnSRpNLkNOgiltTfeaB50lugkSaOJLuhq2+h+8xv41Kd8Hp0kqblcB93rr8OH\nPgT33jvUMaW2jW7PHnjoIWdGkSQ1l+ugGxwM81hu2gTf/W54r7bqssqZUSRJzWQNusOBXmA98COg\n0dDthcBPgSeBJ4DPtHrwzsq8LUccAZdcEtZrg27WLFi2zOfRSZKayxp0VxGC7m3Ajyvb9fqB/wS8\nAzgD+HPgpFYO/ta3wvvfD3/wBzB7dnivto1uYCCU9Aw5SVIzWYNuGbC6sr4aOL/BPi8Dj1TWfwc8\nDRzTysFnzYKvfS2U7BxeIElKI2vQzQO2VNa3VLZHsgg4DVjbysGbDS8ol8Nr716DTpI0slaeXtAL\nHN3g/S/XbZcrr2ZmA7cDnyWU7EY10ji6/v7hJT1JkhppJejeP8JnWwgh+DLwZuCVJvtNBe4A/gG4\nq9nBenp69q8Xi0VmzCiyZ0/joHOwuCS1l1KpRKlUGvPjFjJ+/6+AbcAKQkeULg7skFIgtN9tI3RK\naaZcLg8vEH74w3DRRXDbbWH9wgvDIPEzz4R16+CEE2D79oz/AklSLhXC89my5lTmNrqvE0p864Gz\nKtsQOpv8oLL+e8DHgN8H1lVe57Zy8JHa6CzRSZJakfUJ49uBcxq8/xLwgcr6L0gZqCO10dnjUpLU\nilx35agNuvonjH/xi/DKK85zKUkaWa6D7he/gOuug1/+EnbvDu9VB4y/8EKY69J5LiVJI8l10PX1\nhUB75RW48cbwXrWNrto+5zyXkqSR5DroqmE2Zw58+tNhvVp1+ZWvwOGHO8+lJGlkuQ66Cy+Ek06C\nQw6Ba64J7XH/8i8h6KZPh3e+05CTJI0s10E3Z054Ht3u3fDEE6E97nOfG+p16fACSdJoch1006aF\n8XLVoQXd3TBjRijVffnLQz0xJUlqJus4uoOqGnQnnQRTp8L3vx9KeIOD8PDDMH/+RF+hJCnvogg6\nCI/r6eoKj+6B8Ky6d71r4q5NkhSHKKou9+2Djo7w3ne/G5af/zy86U0Td22SpDhEE3SdlbLnkUeG\nZUeHnVEkSaOLIugGBoZKdNWOKXv2ONelJGl0UQRdbdUlhNLd7t2W6CRJo4uiM0pt1SWE0LvtNti2\nDZ58Em65xYHjkqTGcl2i++u/hrVr4cUXYdeuofc7O2HrVti40UmdJUkjy3XQbd4MO3aEkOvpGXq/\no2P4IHIndZYkNZPrqsuZM8Ny2jT4y78cer+zE5YsgWeecVJnSdLIcl2iW7EizHc5dy4cdtjQ+x0d\noSfmpZcacpKkkeU66I48Eo45Jjx/zl6XkqQ0ch10I/W63L3bcXSSpNFFEXS1A8bBEp0kqXVRBF39\ngPFqic6gkySNJoqgGxgYXnVZLdFZdSlJGk2uhxdUg25w0BKdJCmdKIKuUGjcRmeJTpI0mlwHXWdn\nqLasrldZopMktSpLG93hQC+wHvgRMNLQ7Q5gHXBPkhMUCiHgGo2j27cP/uzPYOlS6OtLfO2SpEki\nS9BdRQi6twE/rmw381ngKaCc9CSDg2H5gQ8MBVo19B5+2EmdJUkjyxJ0y4DVlfXVwPlN9lsALAX+\nFigkPcm+fWG5Zg1ccklYry3dOamzJGkkWYJuHrClsr6lst3IN4ErgMEM5wJCVSbAhg1hOXcu3H67\n811KkpobLeh6gccbvJbV7VemcbXkB4FXCO1ziUtzMFR6O+00uPnmsP7GG2G5dStccUWao0qSJovR\nel2+f4TPtgBHAy8DbyYEWr0lhFBcCswADgX+Hvh4owP21Dx0rlgsUiwWOeQQ2LkTfvKToZLb1Klh\n+e53W20pSe2iVCpRKpXG/LipSlkVfwVsA1YQOqJ0MXKHlPcBXwA+1OTzcrl8YKGwqys8fLX2owsu\ngLvugu3bhz++R5LUPgqhvSpLTgHZ2ui+TijxrQfOqmwDHAP8oMl3Eve6LDT4J86YEZaGnCRpNFkG\njG8Hzmnw/kvABxq8/78qr0QaBZ0kSa3K9aTOkiRllfug27UrLGtnQFm79sD3JElqJPdB12jA+LZt\nQ+85K4okaSS5D7raiZur7XXViZ7nzIFvfGP8r0mSFI/cB92SJWFZO2D81FPDcscOB4xLkkaW+6C7\n7Tb4yEeGDxifMycsnedSkjSaPHXebzhgvJG+vtA2t2qV81xKUrsaqwHjUQadJKn95WFmFEmScs+g\nkyS1NYNOktTWDDpJUlsz6CRJbc2gkyS1NYNOktTWDDpJUlsz6CRJbc2gkyS1NYNOktTWDDpJUlsz\n6CRJbc2gkyS1NYNOktTWDDpJUlsz6CRJbc2gkyS1tSxBdzjQC6wHfgR0NdmvC7gdeBp4Cjgjwzkl\nSUokS9BdRQi6twE/rmw38i3gXuAk4GRC4GmMlEqlib6EKHnf0vG+peN9m1hZgm4ZsLqyvho4v8E+\nc4Azgb+rbA8AOzKcU3X8D5SO9y0d71s63reJlSXo5gFbKutbKtv13gpsBW4G/i/wHWBWhnNKkpTI\naEHXCzze4LWsbr9y5VWvE3gX8O3KchfNqzglSRpzhQzffQYoAi8DbwZ+CpxYt8/RwC8JJTuA9xKC\n7oMNjvc8cFyG65EktZcNwPFZD9KZ4bt3A58AVlSWdzXY52VgI6HDynrgHODJJsfL/I+RJGksHQ7c\nz4HDC44BflCz3ynAr4BHgTsJHVQkSZIkSTE4l9Ce9xxwZZN9rq98/ihwWsLvtrO0924hoc30SeAJ\n4DMH9zJzJ8vPHEAHsA6452BdYE5luW+TeWKILPftS4T/p48DtwDTD95l5s5o9+1EQh+PPcDnE353\nXHUQOpksAqYCjxAGjtdaShhQDnA68ECC77azLPfuaODUyvps4NkG321XWe5b1eeA/0loh54sst63\n1cCfVNY7mTxNFFnu2yLgBYbC7XuE/g6TQSv3bS7QDVzD8KBLnA0He67LxZULehHoB24Fzqvbp3bg\n+VrCX4ZHt/jddpb23s0jdAJ6pPL+7wh/ZR9zcC83N7LcN4AFhF9Mf0u2XsmxyXLfJvPEEFnu287K\nd2YR/jiYBWw+6FecD63ct63AQ5XPk353mIMddPMJvS6rNlXea2WfY1r4bjtLe+8W1O2ziFBVsnaM\nry+vsvzMAXwTuAIYPFgXmFNZft4m88QQWX7etgPXAb8FXgL6CB38JoNW7tuYffdgB12jQeSNTKa/\nnFuV9t7Vfm82od3ks4SS3WSQ9r4VCOM7XyG0z022n8ksP2+TeWKILL/jjgP+I+GP0WMI/1//eGwu\nK/davW9j8t2DHXSbCR0jqhYS0nekfRZU9mnlu+0s7b2rVn1MBe4A/oHGYxzbVZb7toRQzfRr4B+B\ns4C/P2hXmi9Z7tumyutXlfdvJwTeZJDlvnUD/wfYRqjuvZPwMzgZZPn9nrts6CSMbF8ETGP0htoz\nGGqobeW77SzLvSsQfkF/86BfZf5kuW+13sfk6nWZ9b79jDAxBEAPYSKJySDLfTuV0Ct6JuH/7Grg\nzw/u5eZGkt/vPQzvjJLLbPi3hF5/zxO60gJcVnlV3VD5/FGG/yXY6LuTSdp7915CG9MjhGq4dYTu\nuJNFlp+5qvcxuXpdQrb7Npknhshy377I0PCC1YSamMlitPt2NKEtbgfwGqEtc/YI35UkSZIkSZIk\nSZIkSZIkSZIkSZIkSZIktZv/D1KP9ntDNK4lAAAAAElFTkSuQmCC\n", + "text": [ + "" + ] + } + ], + "prompt_number": 14 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "txind = Utils.closestPoints(mesh, [0., 30.], gridLoc='CC')\n", + "q = Utils.sdiag(1/mesh.vol)*np.zeros(mesh.nC)\n", + "q[txind] = 1." + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 15 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Ainv = SolverLU(An)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 16 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "p = np.zeros((mesh.nC, time.size))" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 17 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Proj = mesh.getInterpolationMat(np.r_[3., 30], 'CC')" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 18 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%%time\n", + "b0 = np.zeros(mesh.nC+mesh.nF)\n", + "for i in range(time.size-1):\n", + " s1 = ricker(400, time[i+1], 0.0025)\n", + " s0 = ricker(400, time[i], 0.0025)\n", + " s = np.r_[q*(s1-s0)*1/dt, np.zeros(mesh.nF)]\n", + " bn = Ainv*(s-Bn*b0) \n", + " p[:,i+1] = bn[0:mesh.nC]\n", + " b0 = bn.copy()" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Wall time: 3min 22s\n" + ] + } + ], + "prompt_number": 42 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "data = Proj*p" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 43 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import JSAnimation" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 44 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import urllib" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 45 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min(), mesh.vectorCCy.max()]\n", + "extent[:2]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 46, + "text": [ + "[-124.75, 124.75]" + ] + } + ], + "prompt_number": 46 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from JSAnimation import IPython_display\n", + "from matplotlib import animation\n", + "fig, ax = subplots(1,2, figsize = (16, 8))\n", + "ax[0].set_xlabel('Easting (m)')\n", + "ax[0].set_ylabel('Depth (m)')\n", + "extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min(), mesh.vectorCCy.max()]\n", + "ax[0].set_xlim(extent[:2])\n", + "ax[0].set_ylim(extent[2:])\n", + "ax[1].plot(Utils.mkvc(data), time, 'k--',)\n", + "ax[1].set_xlim(data.min(), data.max())\n", + "ax[1].set_ylim(time.min(), 0.04)\n", + "ax[1].invert_yaxis()\n", + "\n", + "line, = ax[1].plot([], [], color=\"black\", lw=2)\n", + "nskip = 20\n", + "def animate(i_id):\n", + " icount = i_id*nskip\n", + " frame = ax[0].imshow(np.flipud(p[:,icount].reshape((500, 500), order = 'F').T), cmap = 'binary', extent=extent) \n", + " tx = ax[0].plot(mesh.gridCC[txind,0], mesh.gridCC[txind,1], 'k.', ms = 10)\n", + " rx = ax[0].plot(10, 30, 'r.', ms = 10)\n", + " text_tx = ax[0].text(mesh.gridCC[txind,0]-15., mesh.gridCC[txind,1], 'Tx', fontsize = 18)\n", + " text_tx = ax[0].text(10+5., 30, 'Rx', fontsize = 18, color=\"red\")\n", + " ax[0].plot(np.r_[-12.5, 12.5, 12.5, -12.5, -12.5], np.r_[-12.5, -12.5, 12.5, 12.5, -12.5], 'w-', lw=2)\n", + " line.set_data([Utils.mkvc(data)[:icount]], [time[:icount]])\n", + " return frame, line\n", + "animation.FuncAnimation(fig, animate, frames=40, interval=40, blit=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Once \n", + " Loop \n", + " Reflect \n", + "
\n", + "
\n", + "\n", + "\n", + "\n" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 48, + "text": [ + "" + ] + } + ], + "prompt_number": 48 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "!ipython nbconvert --to html SeismicEx.ipynb" + ], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/notebooks/SeismicFWD.html b/notebooks/SeismicFWD.html new file mode 100644 index 0000000..224507d --- /dev/null +++ b/notebooks/SeismicFWD.html @@ -0,0 +1,22715 @@ + + + + + +SeimsicFWD + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+
+
+In [23]: +
+
+
+
from simpegseis import *
+%pylab inline
+from SimPEG.Utils import mkvc
+
+ +
+
+
+ +
+
+ + +
+
+
+Populating the interactive namespace from numpy and matplotlib
+
+
+
+
+ +
+
+ +
+
+
+
+In [32]: +
+
+
+
# Step1: set time
+time = np.linspace(0, 0.04, 2**9)
+dt = time[1]-time[0]
+# Step2: set Tx and Rx
+options={'tlag':0.0025, 'fmain':400} # You need to set waveform to set Tx
+rx = AcousticRx(np.vstack((np.r_[0, 1], np.r_[0, 1])))
+tx = AcousticTx(np.r_[0, 1], time, [rx], **options)
+# Step3: set survey (pass txlist)
+survey = SurveyAcoustic([tx])
+wave = tx.RickerWavelet()
+# Step4: set mesh
+cs = 0.5
+hx = np.ones(150)*cs
+hy = np.ones(150)*cs
+mesh = Mesh.TensorMesh([hx, hy], 'CC')
+# Step5: set problem (pass mesh) and pair with survey
+prob = AcousticProblemPML(mesh)
+prob.pair(survey)
+# Step6: set boundary
+prob.setPMLBC(20, dt, bcflag='all', const=0.)
+prob.storefield = True
+# Step7: set velocity model and check stability
+v = np.ones(mesh.nC)*2000.
+prob.stabilitycheck(v, time, 100.)
+
+ +
+
+
+ +
+
+ + +
+
+
+You are good to go:)
+>> Stability information
+   dt: 7.83e-05 s
+   Optimal dt: 1.25e-04 s
+   Cell per wavelength (G): 4.00e+01
+   Optimal G: 1.60e+01
+
+
+
+
+ +
+
+ +
+
+
+
+In [25]: +
+
+
+
# Step8: run forward
+U = prob.fields(v)
+
+ +
+
+
+ +
+
+ + +
+
+
+>> Start Computing Acoustic Wave
+>> dt: 7.83e-05 s
+>> Optimal dt: 1.25e-04 s
+>> Main frequency, fmain: 1.00e+02 Hz
+>> Cell per wavelength (G): 4.00e+01
+  Tx at (   0.00,    0.00):    1/   1
+>>Elapsed time: 2.55e+00 s
+
+
+
+
+ +
+
+ +
+
+
+
+In [26]: +
+
+
+
# Step9: project data
+data = survey.projectFields(U)
+
+ +
+
+
+ +
+
+
+
+In [27]: +
+
+
+
print data[0].shape
+
+ +
+
+
+ +
+
+ + +
+
+
+(2L, 512L)
+
+
+
+
+ +
+
+ +
+
+
+
+In [28]: +
+
+
+
from JSAnimation import IPython_display
+from matplotlib import animation
+fig, ax = plt.subplots(1,1, figsize = (8, 8))
+ax.set_xlabel('Easting (m)', fontsize = 16)
+ax.set_ylabel('Depth (m)', fontsize = 16)
+extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]
+ax.set_xlim(extent[:2])
+ax.set_ylim(extent[2:])
+
+nskip = 20
+def animate(i_id):
+    icount = i_id*nskip
+    frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent)
+    return frame
+animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)
+
+ +
+
+
+ +
+
+ + +
+ Out[28]:
+ +
+ + + +
+ +
+ +
+ + + + + + + + + +
+ Once + Loop + Reflect +
+
+ + + + +
+ +
+ +
+
+ +
+
+
+
+In [29]: +
+
+
+
prob.setPMLBC(30, dt, bcflag='all', const=1.)
+
+ +
+
+
+ +
+
+
+
+In [30]: +
+
+
+
U = prob.fields(v)
+
+ +
+
+
+ +
+
+ + +
+
+
+>> Start Computing Acoustic Wave
+>> dt: 7.83e-05 s
+>> Optimal dt: 1.25e-04 s
+>> Main frequency, fmain: 1.00e+02 Hz
+>> Cell per wavelength (G): 4.00e+01
+  Tx at (   0.00,    0.00):    1/   1
+>>Elapsed time: 2.75e+00 s
+
+
+
+
+ +
+
+ +
+
+
+
+In [31]: +
+
+
+
from JSAnimation import IPython_display
+from matplotlib import animation
+fig, ax = plt.subplots(1,1, figsize = (8, 8))
+ax.set_xlabel('Easting (m)', fontsize = 16)
+ax.set_ylabel('Depth (m)', fontsize = 16)
+extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]
+ax.set_xlim(extent[:2])
+ax.set_ylim(extent[2:])
+
+nskip = 20
+def animate(i_id):
+    icount = i_id*nskip
+    frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent)
+    return frame
+animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)
+
+ +
+
+
+ +
+
+ + +
+ Out[31]:
+ +
+ + + +
+ +
+ +
+ + + + + + + + + +
+ Once + Loop + Reflect +
+
+ + + + +
+ +
+ +
+
+ +
+
+
+ + diff --git a/notebooks/SeismicFWD.ipynb b/notebooks/SeismicFWD.ipynb new file mode 100644 index 0000000..21cc7bf --- /dev/null +++ b/notebooks/SeismicFWD.ipynb @@ -0,0 +1,20989 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:43cef8feaa23d5020ffbd5874e54f80865dd6632477d24adbab1e45f5867176c" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from simpegseis import *\n", + "%pylab inline\n", + "from SimPEG.Utils import mkvc" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "prompt_number": 23 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Step1: set time\n", + "time = np.linspace(0, 0.04, 2**9)\n", + "dt = time[1]-time[0]\n", + "# Step2: set Tx and Rx\n", + "options={'tlag':0.0025, 'fmain':400} # You need to set waveform to set Tx\n", + "rx = AcousticRx(np.vstack((np.r_[0, 1], np.r_[0, 1])))\n", + "tx = AcousticTx(np.r_[0, 1], time, [rx], **options)\n", + "# Step3: set survey (pass txlist)\n", + "survey = SurveyAcoustic([tx])\n", + "wave = tx.RickerWavelet()\n", + "# Step4: set mesh\n", + "cs = 0.5\n", + "hx = np.ones(150)*cs\n", + "hy = np.ones(150)*cs\n", + "mesh = Mesh.TensorMesh([hx, hy], 'CC')\n", + "# Step5: set problem (pass mesh) and pair with survey\n", + "prob = AcousticProblemPML(mesh)\n", + "prob.pair(survey)\n", + "# Step6: set boundary\n", + "prob.setPMLBC(20, dt, bcflag='all', const=0.)\n", + "prob.storefield = True\n", + "# Step7: set velocity model and check stability\n", + "v = np.ones(mesh.nC)*2000.\n", + "prob.stabilitycheck(v, time, 100.)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "You are good to go:)\n", + ">> Stability information\n", + " dt: 7.83e-05 s\n", + " Optimal dt: 1.25e-04 s\n", + " Cell per wavelength (G): 4.00e+01\n", + " Optimal G: 1.60e+01\n" + ] + } + ], + "prompt_number": 32 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Step8: run forward\n", + "U = prob.fields(v)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + ">> Start Computing Acoustic Wave\n", + ">> dt: 7.83e-05 s\n", + ">> Optimal dt: 1.25e-04 s\n", + ">> Main frequency, fmain: 1.00e+02 Hz\n", + ">> Cell per wavelength (G): 4.00e+01\n", + " Tx at ( 0.00, 0.00): 1/ 1\n", + ">>Elapsed time: 2.55e+00 s" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 25 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Step9: project data\n", + "data = survey.projectFields(U)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 26 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print data[0].shape" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "(2L, 512L)\n" + ] + } + ], + "prompt_number": 27 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from JSAnimation import IPython_display\n", + "from matplotlib import animation\n", + "fig, ax = plt.subplots(1,1, figsize = (8, 8))\n", + "ax.set_xlabel('Easting (m)', fontsize = 16)\n", + "ax.set_ylabel('Depth (m)', fontsize = 16)\n", + "extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]\n", + "ax.set_xlim(extent[:2])\n", + "ax.set_ylim(extent[2:])\n", + "\n", + "nskip = 20\n", + "def animate(i_id):\n", + " icount = i_id*nskip\n", + " frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent) \n", + " return frame\n", + "animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Once \n", + " Loop \n", + " Reflect \n", + "
\n", + "
\n", + "\n", + "\n", + "\n" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 28, + "text": [ + "" + ] + } + ], + "prompt_number": 28 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "prob.setPMLBC(30, dt, bcflag='all', const=1.)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 29 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "U = prob.fields(v)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + ">> Start Computing Acoustic Wave\n", + ">> dt: 7.83e-05 s\n", + ">> Optimal dt: 1.25e-04 s\n", + ">> Main frequency, fmain: 1.00e+02 Hz\n", + ">> Cell per wavelength (G): 4.00e+01\n", + " Tx at ( 0.00, 0.00): 1/ 1\n", + ">>Elapsed time: 2.75e+00 s" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 30 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from JSAnimation import IPython_display\n", + "from matplotlib import animation\n", + "fig, ax = plt.subplots(1,1, figsize = (8, 8))\n", + "ax.set_xlabel('Easting (m)', fontsize = 16)\n", + "ax.set_ylabel('Depth (m)', fontsize = 16)\n", + "extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]\n", + "ax.set_xlim(extent[:2])\n", + "ax.set_ylim(extent[2:])\n", + "\n", + "nskip = 20\n", + "def animate(i_id):\n", + " icount = i_id*nskip\n", + " frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent) \n", + " return frame\n", + "animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Once \n", + " Loop \n", + " Reflect \n", + "
\n", + "
\n", + "\n", + "\n", + "\n" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 31, + "text": [ + "" + ] + } + ], + "prompt_number": 31 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# !ipython nbconvert --to html SeimsicFWD.ipynb" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 34 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..6b7df82 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,4 @@ +numpy +scipy +ipython +matplotlib diff --git a/simpegseis/SeisAcousticTime.py b/simpegseis/SeisAcousticTime.py index 5b5f608..5637a79 100644 --- a/simpegseis/SeisAcousticTime.py +++ b/simpegseis/SeisAcousticTime.py @@ -109,7 +109,7 @@ class AcousticProblemSponge(Problem.BaseProblem): self.mesh.setCellGradBC('dirichlet') Utils.setKwargs(self, **kwargs) - def setSpongeBC(self, npad, dt, bcflag="all"): + def setSpongeBC(self, npad, dt, bcflag="all", const=1.): #TODO: performance of abosrbing self.bcflag = bcflag ax = self.mesh.vectorCCx[-npad] @@ -140,7 +140,7 @@ class AcousticProblemSponge(Problem.BaseProblem): temp[temp>1.] = 1. f = 1.- temp*0.1 - self.sig = (1.-f)/f*2./dt + self.sig = (1.-f)/f*2./dt*const def stabilitycheck(self, v, time, fmain): @@ -274,7 +274,7 @@ class AcousticProblemPML(Problem.BaseProblem): self.mesh.setCellGradBC('dirichlet') Utils.setKwargs(self, **kwargs) - def setPMLBC(self, npad, dt, bcflag="all"): + def setPMLBC(self, npad, dt, bcflag="all", const=1.): #TODO: performance of abosrbing self.bcflag = bcflag ax = self.mesh.vectorCCx[-npad] @@ -304,8 +304,8 @@ class AcousticProblemPML(Problem.BaseProblem): fx = 1.-tempx*0.1 fy = 1.-tempy*0.1 - self.sigx = (1-fx)/fx*2./dt - self.sigy = (1-fy)/fy*2./dt + self.sigx = (1-fx)/fx*2./dt*const + self.sigy = (1-fy)/fy*2./dt*const def stabilitycheck(self, v, time, fmain):