Fix import bug of sys

This commit is contained in:
GudniRos
2015-04-03 10:44:56 -07:00
parent d7aa612a23
commit 6b3f9b9478
10 changed files with 6988 additions and 1779 deletions
+95
View File
@@ -0,0 +1,95 @@
{
"metadata": {
"name": "",
"signature": "sha256:74c8f271a46170528d2e89be943c0f766fa137ce8e32b8bc8ad86623db691b9d"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Maxwell's equations are\n",
"\n",
"\\begin{align}\n",
"\\frac{\\partial \\vec{B}}{\\partial t} &= - \\nabla \\times \\vec{E} \\hspace{2cm} \\text{Faraday equations}\\\\\n",
"\\frac{\\partial \\vec{D}}{\\partial t} &= - \\nabla \\times \\vec{H} - \\vec{J} - \\vec{s} \\hspace{2cm} \\text{Ampere's law.}\n",
"\\end{align}\n",
"\n",
"where $\\vec{E}$ is the electric field, $\\vec{H}$ is the magnetic field, $\\vec{J}$ is the electric flux, $\\vec{D}$ is the displacement flux, $\\vec{B}$ is the magnetic flux and $\\vec{s}$ is the source."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The constitutive realtations between the fields and flux are as follows.\n",
"\n",
"\\begin{align}\n",
"\\vec{J} &= \\sigma \\vec{E}\\\\\n",
"\\vec{B} &= \\mu \\vec{H}\\\\\n",
"\\vec{D} &= \\epsilon \\vec{E}\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the contitutive relations, the Maxwell's equations can written in terms of the electric field and the magnetic flux \n",
"\\begin{align}\n",
"\\mu \\frac{\\partial \\vec{B}}{\\partial t} &= - \\nabla \\times \\vec{E} \\\\\n",
"\\epsilon \\frac{\\partial \\vec{E}}{\\partial t} &= - \\nabla \\times \\frac{1}{\\mu} \\vec{B} - \\sigma \\vec{E} - \\vec{s} \n",
"\\end{align}\n",
"\n",
"Futher by droping the $\\epsilon$ term and using a $e^{-i\\omega t}$ Fourier time relation, Maxwell's equations can be written in the frequency domain as\n",
"\n",
"\\begin{align}\n",
"-i \\omega \\mu \\hat{\\vec{B}} &= - \\nabla \\times \\hat{\\vec{E}} \\\\\n",
"- \\vec{s} &= - \\nabla \\times \\frac{1}{\\mu} \\hat{\\vec{B}} - \\sigma \\hat{\\vec{E}}\n",
"\\end{align}\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"The weak form is \n",
"\\begin{align}\n",
"-i \\omega \\mu (\\hat{\\vec{B}},\\hat{\\vec{F}}) + (\\nabla \\times \\hat{\\vec{E}},\\hat{\\vec{F}}) &= 0 \\\\\n",
"(\\nabla \\times \\mu^{-1} \\hat{\\vec{B}},\\hat{\\vec{W}}) + (\\sigma \\hat{\\vec{E}},\\hat{\\vec{W}}) &= (\\hat{\\vec{s}},\\hat{\\vec{W}})\n",
"\\end{align}\n",
"\n",
"\n",
"$ (\\mu^{-1} B, \\nabla \\times W) + (\\sigma E, W) = (\\mu^{-1}B W )|_0^{end} $\n",
"\n",
"$ (i \\omega B,F) + (\\nabla \\times E , f) = 0 $"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
+95
View File
@@ -0,0 +1,95 @@
{
"metadata": {
"name": "",
"signature": "sha256:74c8f271a46170528d2e89be943c0f766fa137ce8e32b8bc8ad86623db691b9d"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Maxwell's equations are\n",
"\n",
"\\begin{align}\n",
"\\frac{\\partial \\vec{B}}{\\partial t} &= - \\nabla \\times \\vec{E} \\hspace{2cm} \\text{Faraday equations}\\\\\n",
"\\frac{\\partial \\vec{D}}{\\partial t} &= - \\nabla \\times \\vec{H} - \\vec{J} - \\vec{s} \\hspace{2cm} \\text{Ampere's law.}\n",
"\\end{align}\n",
"\n",
"where $\\vec{E}$ is the electric field, $\\vec{H}$ is the magnetic field, $\\vec{J}$ is the electric flux, $\\vec{D}$ is the displacement flux, $\\vec{B}$ is the magnetic flux and $\\vec{s}$ is the source."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The constitutive realtations between the fields and flux are as follows.\n",
"\n",
"\\begin{align}\n",
"\\vec{J} &= \\sigma \\vec{E}\\\\\n",
"\\vec{B} &= \\mu \\vec{H}\\\\\n",
"\\vec{D} &= \\epsilon \\vec{E}\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the contitutive relations, the Maxwell's equations can written in terms of the electric field and the magnetic flux \n",
"\\begin{align}\n",
"\\mu \\frac{\\partial \\vec{B}}{\\partial t} &= - \\nabla \\times \\vec{E} \\\\\n",
"\\epsilon \\frac{\\partial \\vec{E}}{\\partial t} &= - \\nabla \\times \\frac{1}{\\mu} \\vec{B} - \\sigma \\vec{E} - \\vec{s} \n",
"\\end{align}\n",
"\n",
"Futher by droping the $\\epsilon$ term and using a $e^{-i\\omega t}$ Fourier time relation, Maxwell's equations can be written in the frequency domain as\n",
"\n",
"\\begin{align}\n",
"-i \\omega \\mu \\hat{\\vec{B}} &= - \\nabla \\times \\hat{\\vec{E}} \\\\\n",
"- \\vec{s} &= - \\nabla \\times \\frac{1}{\\mu} \\hat{\\vec{B}} - \\sigma \\hat{\\vec{E}}\n",
"\\end{align}\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"The weak form is \n",
"\\begin{align}\n",
"-i \\omega \\mu (\\hat{\\vec{B}},\\hat{\\vec{F}}) + (\\nabla \\times \\hat{\\vec{E}},\\hat{\\vec{F}}) &= 0 \\\\\n",
"(\\nabla \\times \\mu^{-1} \\hat{\\vec{B}},\\hat{\\vec{W}}) + (\\sigma \\hat{\\vec{E}},\\hat{\\vec{W}}) &= (\\hat{\\vec{s}},\\hat{\\vec{W}})\n",
"\\end{align}\n",
"\n",
"\n",
"$ (\\mu^{-1} B, \\nabla \\times W) + (\\sigma E, W) = (\\mu^{-1}B W )|_0^{end} $\n",
"\n",
"$ (i \\omega B,F) + (\\nabla \\times E , f) = 0 $"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
+1880 -526
View File
File diff suppressed because one or more lines are too long
File diff suppressed because it is too large Load Diff
File diff suppressed because one or more lines are too long
+458 -440
View File
@@ -1,447 +1,465 @@
{
"metadata": {
"name": "",
"signature": "sha256:74ba4ac0804b189d65dce7f7af6f88a605b83052cf38c5acb492600297a13398"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
"cells": [
{
"cells": [
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"%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": "markdown",
"metadata": {},
"source": [
"Simple notebook on performing a 1D MT problem.\n"
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#from IPython.display import Latex"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Maxwell's equations in 1D are as follows\n",
"\n",
"\n",
"$i \\omega b = - \\partial_z e $ \n",
"\n",
"$ s = \\partial_z(\\mu^{-1} b) - \\sigma(z) e$\n",
"\n",
"$b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0$\n",
"\n",
"\n",
"\n",
"where $e = \\widehat{\\overrightarrow{E}}_x $ and $b = \\widehat{\\overrightarrow{B}}_y$\n",
"\n",
"In weak form the equations become\n",
"\n",
"$i \\omega(b,f) = - (\\partial_ze,f)$ \n",
"\n",
"$(s,w) = -(\\mu^{-1} b, \\partial_z w) - (\\sigma(z) e, w)$\n",
"\n",
"\n",
"where f and w are abritrary functions living in same discritizational space as b and e, respectivily.\n",
"\n",
"We consider e on nodes and b on cell centers. This way the derivative of any nodal function becomes\n",
"\n",
"$ (\\partial_z u)_k \\approx $\n",
"\n",
"$h_{k}^{-1} ( u_{k+\\frac{1}{2}} - u_{k-\\frac{1}{2}}) + O(h^2) $\n",
"\n",
"Matrix form\n",
"\n",
"$ e_z \\approx \\textbf{L}^{-1} \\textbf{G} e = \\begin{bmatrix} h_1^{-1} & & & \\\\\\\\ & h_2^{-1} & & \\\\\\\\ & & \\ddots & \\\\\\\\ & & & h_n^{-1} \\end{bmatrix}^{(n,n)}\n",
"\\begin{bmatrix} -1 & 1 & & & & \\\\\\\\ & -1 & 1 & & & \\\\\\\\ & & \\ddots & \\ddots & \\\\\\\\ & & & -1 & 1 \\end{bmatrix}^{(n,n+1)} \n",
"\\begin{bmatrix} e_1 \\\\\\\\ \\\\\\\\ \\vdots \\\\\\\\ \\\\\\\\ e_{n+1} \\end{bmatrix}^{(n+1,1)} $\n",
"\n",
"where $ \\textbf{L} = diag(h) $ is the cell size and $ \\textbf{G}$ is the gradient operator with -1,1 representing the topology of the mesh, taking the difference between adjoint cells.\n",
"\n",
"We need to compute 2 inner products, on cell centers and from nodes to cell centers.\n",
"\n",
"Cell centers inner product is\n",
"\n",
"$ (b,f) \\approx \\sum\\limits_k h_k \\textbf{b}_k \\textbf{f}_k + O(h^2) $\n",
"\n",
"and in matrix from\n",
"\n",
"$ (b,) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f}$ and $ (\\mu^{-1} b,f) \\approx \\textbf{b}^T \\textbf{M}_{\\mu}^f \\textbf{f}$ \n",
"\n",
"where $ \\textbf{M}_{\\mu}^f = diag(\\textbf{h} \\odot \\mu^{-1}) $ and $ \\textbf{M}^f = diag(\\textbf{h}) $ are the matrices.\n",
"Nodes to cell centers inner product is\n",
"\n",
"$ (\\sigma e, w) \\approx \\sum\\limits_k \\frac{h_k \\sigma_k}{4} ( e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} + e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} )$\n",
"\n",
"and in matrix from\n",
"\n",
"$ (\\sigma e, w ) \\approx (\\textbf{h} \\odot \\sigma )^T ( \\textbf{A}_v (\\textbf{e} \\odot \\textbf{w} = \\textbf{w}^T diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) \\textbf{e} $\n",
"\n",
"Here $\\odot$ is a point wise Hadamard product and $\\textbf{A}_v$ is the averaging operator/matrix from nodes to cell centers\n",
"\n",
"$ \\textbf{A}_v = \\begin{bmatrix} \\frac{1}{2} & \\frac{1}{2} & & \\\\\\\\ & \\ddots & \\ddots & \\\\\\\\ & & \\frac{1}{2} & \\frac{1}{2} \\end{bmatrix} ^{(n+1 , n)} $\n",
"\n",
"The sigma mass matrix is defined as \n",
"\n",
"$ \\textbf{M}_{\\sigma}^{e} = diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) $\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the MT problem there is no source in the domain, so $ s = 0 $. How ever the boundary conditions provide the right hand side where \n",
"\n",
"$ (\\partial_z \\mu^{-1} b, w ) = - (\\mu^{-1} b, \\partial_z w ) + (\\mu^{-1} b w )|_0^{end} $\n",
"\n",
"where \n",
"\n",
"$ (\\mu^{-1} b w )|_0^{end} = \\textbf{bc}^T (\\textbf{BC w}) $\n",
"\n",
"here $\\textbf{BC}$ is an matrix operator that extracts the boundary elements from $ \\textbf{w}$ and $\\textbf{bc} $ are the known boundary condintions. For the 1D case with homogenous boundary conditions we have \n",
"\n",
"$ \\textbf{B} = \\begin{bmatrix} -1 & 0 \\\\\\\\ \\vdots & \\vdots \\\\\\\\ 0 & 1 \\end{bmatrix} $ \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The weak form is \n",
"\n",
"$ (\\mu^{-1} b, \\partial_z w) + (\\sigma e, w) = (\\mu^{-1} b w )|_0^{end} $\n",
"\n",
"$ (i \\omega b,f) + (\\partial_z e , f) = 0 $\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the above matrix represntation we get the Maxwells equations in following form\n",
"\n",
"$ i \\omega \\textbf{f}^T \\textbf{M}^f \\textbf{b} + \\textbf{f}^T \\textbf{M}^f \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
"\n",
"$ \\textbf{w}^T \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} \\textbf{b} + \\textbf{w}^T \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{w}^T \\textbf{bc}^T \\textbf{BC} $ \n",
"\n",
"Here we use that \n",
"\n",
"$ (\\textbf{b},\\textbf{f}) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f} = \\textbf{f}^T \\textbf{M}^f \\textbf{b} $ \n",
"\n",
"since $\\textbf{M}^f$ is a symmetric diaoganal matrix of size n by n and $\\textbf{b}$ and $\\textbf{f}$ are vectors of length n."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We eliminate testing vectors and get system of equations to solve\n",
"\n",
"$ i \\omega \\textbf{b} + \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
"\n",
"$ \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} \\textbf{b} + \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{bc}^T \\textbf{BC} $ \n",
"\n",
"and as $ \\textbf{A} \\textbf{x} = \\textbf{bc} $ system that we will solve, where \n",
"\n",
"$ \\textbf{A} = \\begin{bmatrix} \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} & \\textbf{M}_{\\sigma}^e \\\\\\\\ i \\omega & \\textbf{L}^{-1} \\textbf{G} \\end{bmatrix} $\n",
"\n",
"$ \\textbf{x} = \\begin{bmatrix} \\textbf{b} \\\\\\\\ \\textbf{e} \\end{bmatrix} $\n",
"\n",
"$ \\textbf{bc} = \\begin{bmatrix} \\textbf{bc}^T \\textbf{BC} \\\\\\\\ \\textbf{0} \\end{bmatrix} $ "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import sys\n",
"sys.path.append('C:/GudniWork/Codes/python/simpeg')\n",
"import SimPEG as simpeg, numpy as np, scipy, scipy.sparse as sp"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Efficiency Warning: Interpolation will be slow, use setup.py!\n",
"\n",
" python setup.py build_ext --inplace\n",
" \n"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have\n",
"\n",
"$ i \\omega b = - \\partial_z e \\hspace{1cm} ; \\hspace{1cm} \\partial_z(\\mu^{-1} b) - \\sigma(z) e \\hspace{1cm} ;\\hspace{1cm} b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0 $\n",
"\n",
" \n",
"To deal with boundary: we assume that below depth L both $ \\sigma $ and $ \\mu $ are constants ($ z < - L $). At the boundary we have that \n",
"\n",
"$ e = c \\exp(ikz) \\hspace{0.2cm} where \\hspace{0.2cm} k = \\sqrt{i\\omega\\mu\\sigma} $. \n",
"\n",
"Therefore for $ z < - L $ we have that\n",
"\n",
"$\\omega b - k e = 0 $.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We discretize the e field on the nodes and b field at the cell centers. The system we want to solve is \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\begin{bmatrix} i \\omega & \\frac{\\partial}{\\partial z} \\\\\\\\ \\frac{1}{\\mu} \\frac{\\partial}{\\partial z} & -\\sigma \\end{bmatrix} \n",
"\\begin{bmatrix} b \\\\\\\\ e \\end{bmatrix} = \\begin{bmatrix} s1 \\\\\\\\ s2 \\end{bmatrix}$\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
" # Set up the problem\n",
"mu = 4*np.pi*1e-7\n",
"eps0 = 8.85e-12\n",
"# Frequency\n",
"fr = np.array([1e1]) #np.logspace(0,5,200) #np.array([2000]) #np.logspace(-4,5,82)\n",
"omega = 2*np.pi*fr\n",
"# Mesh\n",
"sig0 = 1e-2\n",
"#L = 3*np.sqrt(2/(mu*omega[0]*sig0))\n",
"#nn=np.ceil(np.log(0.3*L + 1)/np.log(1.3))\n",
"#h = 5*(1.3**(np.arange(nn+1)))\n",
"\n",
"h = np.ones(18)\n",
"x0 = np.array([0])\n",
"#sig = sig0*np.ones((len(h),1)) \n",
"#sig[0:50] = 0.1\n",
"#sig[50:100] = 1\n",
"# Make the mesh\n",
"mesh = simpeg.Mesh.TensorMesh([h],x0)\n",
"sig = np.zeros(mesh.nC) + 1e-8\n",
"sig[mesh.vectorCCx<=0] = 1e-2"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fr,omega"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 4,
"text": [
"(array([ 10.]), array([ 62.83185307]))"
]
}
],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Make the operators\n",
"G = mesh.nodalGrad\n",
"Av = mesh.aveN2CC\n",
"Li = scipy.sparse.spdiags(1/mesh.hx,0,mesh.nNx,mesh.nNx)\n",
"Mmu = scipy.sparse.spdiags(mesh.hx/mu,0,mesh.nCx,mesh.nCx)\n",
"Msig = scipy.sparse.spdiags(Av.T.dot(mesh.hx*sig.ravel()),0,mesh.nNx,mesh.nNx)\n",
"# The boundaries\n",
"bc_b = np.zeros((mesh.nCx,1))\n",
"bc_b[0] = -1 # Set the top b field to 1\n",
"bc_e = np.zeros((mesh.nNx,1))\n",
"# Make the sparse matrix\n",
"bc = sp.vstack((bc_b,bc_e))\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"b = np.empty((mesh.nCx,len(omega)),dtype=np.complex64)\n",
"e = np.empty((mesh.nNx,len(omega)),dtype=np.complex64)\n",
"# Loop all the frequencies\n",
"for nrOm, om in enumerate(omega):\n",
" # Left hand side\n",
" A = sp.vstack((sp.hstack(( -G.conj().T.dot(Mmu), - Msig)), sp.hstack((1j*om*scipy.sparse.identity(mesh.nCx) , G))))\n",
" #A = A.tocsr\n",
" # Solve the system\n",
" bef = scipy.sparse.linalg.spsolve(A,bc)\n",
" # Sort the output\n",
" b[:,nrOm] = bef[0:mesh.nCx]\n",
" e[:,nrOm] = bef[mesh.nCx::]\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stderr",
"text": [
"/home/Gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:90: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format\n",
" SparseEfficiencyWarning)\n"
]
}
],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import matplotlib.pyplot as plt\n",
"# Plot the solution\n",
"z=e[0,:]/(b[0,:]/mu)\n",
"app_res = ((1./(8e-7*np.pi**2))/fr)*np.abs(z)**2\n",
"app_phs = np.arctan(z.imag/z.real)*(180/np.pi)\n",
"ax_res = plt.subplot(2,1,1)\n",
"ax_res.loglog(fr,app_res)\n",
"ax_phs = plt.subplot(2,1,2)\n",
"ax_phs.semilogx(fr,app_phs)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 7,
"text": [
"[<matplotlib.lines.Line2D at 0x7fa8671e1fd0>]"
]
}
],
"prompt_number": 7
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plt.show()"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Calculate the impedance\n",
"z = e[0,:]/(b[0,:]/mu)\n",
"z"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 9,
"text": [
"array([-5714285.5+0.00050081j], dtype=complex64)"
]
}
],
"prompt_number": 9
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"app_res"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 10,
"text": [
"array([ 4.13555827e+17])"
]
}
],
"prompt_number": 10
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 10
}
],
"metadata": {}
"source": [
"%pylab inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Simple notebook on performing a 1D MT problem.\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#from IPython.display import Latex"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Maxwell's equations in 1D are as follows\n",
"\n",
"\n",
"$i \\omega b = - \\partial_z e $ \n",
"\n",
"$ s = \\partial_z(\\mu^{-1} b) - \\sigma(z) e$\n",
"\n",
"$b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0$\n",
"\n",
"\n",
"\n",
"where $e = \\widehat{\\overrightarrow{E}}_x $ and $b = \\widehat{\\overrightarrow{B}}_y$\n",
"\n",
"In weak form the equations become\n",
"\n",
"$i \\omega(b,f) = - (\\partial_ze,f)$ \n",
"\n",
"$(s,w) = -(\\mu^{-1} b, \\partial_z w) - (\\sigma(z) e, w)$\n",
"\n",
"\n",
"where f and w are abritrary functions living in same discritizational space as b and e, respectivily.\n",
"\n",
"We consider e on nodes and b on cell centers. This way the derivative of any nodal function becomes\n",
"\n",
"$ (\\partial_z u)_k \\approx $\n",
"\n",
"$h_{k}^{-1} ( u_{k+\\frac{1}{2}} - u_{k-\\frac{1}{2}}) + O(h^2) $\n",
"\n",
"Matrix form\n",
"\n",
"$ e_z \\approx \\textbf{L}^{-1} \\textbf{G} e = \\begin{bmatrix} h_1^{-1} & & & \\\\\\\\ & h_2^{-1} & & \\\\\\\\ & & \\ddots & \\\\\\\\ & & & h_n^{-1} \\end{bmatrix}^{(n,n)}\n",
"\\begin{bmatrix} -1 & 1 & & & & \\\\\\\\ & -1 & 1 & & & \\\\\\\\ & & \\ddots & \\ddots & \\\\\\\\ & & & -1 & 1 \\end{bmatrix}^{(n,n+1)} \n",
"\\begin{bmatrix} e_1 \\\\\\\\ \\\\\\\\ \\vdots \\\\\\\\ \\\\\\\\ e_{n+1} \\end{bmatrix}^{(n+1,1)} $\n",
"\n",
"where $ \\textbf{L} = diag(h) $ is the cell size and $ \\textbf{G}$ is the gradient operator with -1,1 representing the topology of the mesh, taking the difference between adjoint cells.\n",
"\n",
"We need to compute 2 inner products, on cell centers and from nodes to cell centers.\n",
"\n",
"Cell centers inner product is\n",
"\n",
"$ (b,f) \\approx \\sum\\limits_k h_k \\textbf{b}_k \\textbf{f}_k + O(h^2) $\n",
"\n",
"and in matrix from\n",
"\n",
"$ (b,) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f}$ and $ (\\mu^{-1} b,f) \\approx \\textbf{b}^T \\textbf{M}_{\\mu}^f \\textbf{f}$ \n",
"\n",
"where $ \\textbf{M}_{\\mu}^f = diag(\\textbf{h} \\odot \\mu^{-1}) $ and $ \\textbf{M}^f = diag(\\textbf{h}) $ are the matrices.\n",
"Nodes to cell centers inner product is\n",
"\n",
"$ (\\sigma e, w) \\approx \\sum\\limits_k \\frac{h_k \\sigma_k}{4} ( e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} + e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} )$\n",
"\n",
"and in matrix from\n",
"\n",
"$ (\\sigma e, w ) \\approx (\\textbf{h} \\odot \\sigma )^T ( \\textbf{A}_v (\\textbf{e} \\odot \\textbf{w} = \\textbf{w}^T diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) \\textbf{e} $\n",
"\n",
"Here $\\odot$ is a point wise Hadamard product and $\\textbf{A}_v$ is the averaging operator/matrix from nodes to cell centers\n",
"\n",
"$ \\textbf{A}_v = \\begin{bmatrix} \\frac{1}{2} & \\frac{1}{2} & & \\\\\\\\ & \\ddots & \\ddots & \\\\\\\\ & & \\frac{1}{2} & \\frac{1}{2} \\end{bmatrix} ^{(n+1 , n)} $\n",
"\n",
"The sigma mass matrix is defined as \n",
"\n",
"$ \\textbf{M}_{\\sigma}^{e} = diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) $\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the MT problem there is no source in the domain, so $ s = 0 $. How ever the boundary conditions provide the right hand side where \n",
"\n",
"$ (\\partial_z \\mu^{-1} b, w ) = - (\\mu^{-1} b, \\partial_z w ) + (\\mu^{-1} b w )|_0^{end} $\n",
"\n",
"where \n",
"\n",
"$ (\\mu^{-1} b w )|_0^{end} = \\textbf{bc}^T (\\textbf{BC w}) $\n",
"\n",
"here $\\textbf{BC}$ is an matrix operator that extracts the boundary elements from $ \\textbf{w}$ and $\\textbf{bc} $ are the known boundary condintions. For the 1D case with homogenous boundary conditions we have \n",
"\n",
"$ \\textbf{B} = \\begin{bmatrix} -1 & 0 \\\\\\\\ \\vdots & \\vdots \\\\\\\\ 0 & 1 \\end{bmatrix} $ \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The weak form is \n",
"\n",
"$ (\\mu^{-1} b, \\partial_z w) + (\\sigma e, w) = (\\mu^{-1} b w )|_0^{end} $\n",
"\n",
"$ (i \\omega b,f) + (\\partial_z e , f) = 0 $\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the above matrix represntation we get the Maxwells equations in following form\n",
"\n",
"$ i \\omega \\textbf{f}^T \\textbf{M}^f \\textbf{b} + \\textbf{f}^T \\textbf{M}^f \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
"\n",
"$ \\textbf{w}^T \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} \\textbf{b} + \\textbf{w}^T \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{w}^T \\textbf{bc}^T \\textbf{BC} $ \n",
"\n",
"Here we use that \n",
"\n",
"$ (\\textbf{b},\\textbf{f}) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f} = \\textbf{f}^T \\textbf{M}^f \\textbf{b} $ \n",
"\n",
"since $\\textbf{M}^f$ is a symmetric diaoganal matrix of size n by n and $\\textbf{b}$ and $\\textbf{f}$ are vectors of length n."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We eliminate testing vectors and get system of equations to solve\n",
"\n",
"$ i \\omega \\textbf{b} + \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
"\n",
"$ \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} \\textbf{b} + \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{bc}^T \\textbf{BC} $ \n",
"\n",
"and as $ \\textbf{A} \\textbf{x} = \\textbf{bc} $ system that we will solve, where \n",
"\n",
"$ \\textbf{A} = \\begin{bmatrix} \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} & \\textbf{M}_{\\sigma}^e \\\\\\\\ i \\omega & \\textbf{L}^{-1} \\textbf{G} \\end{bmatrix} $\n",
"\n",
"$ \\textbf{x} = \\begin{bmatrix} \\textbf{b} \\\\\\\\ \\textbf{e} \\end{bmatrix} $\n",
"\n",
"$ \\textbf{bc} = \\begin{bmatrix} \\textbf{bc}^T \\textbf{BC} \\\\\\\\ \\textbf{0} \\end{bmatrix} $ "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Efficiency Warning: Interpolation will be slow, use setup.py!\n",
"\n",
" python setup.py build_ext --inplace\n",
" \n"
]
}
],
"source": [
"import sys\n",
"sys.path.append('C:/GudniWork/Codes/python/simpeg')\n",
"import SimPEG as simpeg, numpy as np, scipy, scipy.sparse as sp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have\n",
"\n",
"$ i \\omega b = - \\partial_z e \\hspace{1cm} ; \\hspace{1cm} \\partial_z(\\mu^{-1} b) - \\sigma(z) e \\hspace{1cm} ;\\hspace{1cm} b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0 $\n",
"\n",
" \n",
"To deal with boundary: we assume that below depth L both $ \\sigma $ and $ \\mu $ are constants ($ z < - L $). At the boundary we have that \n",
"\n",
"$ e = c \\exp(ikz) \\hspace{0.2cm} where \\hspace{0.2cm} k = \\sqrt{i\\omega\\mu\\sigma} $. \n",
"\n",
"Therefore for $ z < - L $ we have that\n",
"\n",
"$\\omega b - k e = 0 $.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We discretize the e field on the nodes and b field at the cell centers. The system we want to solve is \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\begin{bmatrix} i \\omega & \\frac{\\partial}{\\partial z} \\\\\\\\ \\frac{1}{\\mu} \\frac{\\partial}{\\partial z} & -\\sigma \\end{bmatrix} \n",
"\\begin{bmatrix} b \\\\\\\\ e \\end{bmatrix} = \\begin{bmatrix} s1 \\\\\\\\ s2 \\end{bmatrix}$\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
" # Set up the problem\n",
"mu = 4*np.pi*1e-7\n",
"eps0 = 8.85e-12\n",
"# Frequency\n",
"fr = np.array([1e1]) #np.logspace(0,5,200) #np.array([2000]) #np.logspace(-4,5,82)\n",
"omega = 2*np.pi*fr\n",
"# Mesh\n",
"sig0 = 1e-2\n",
"#L = 3*np.sqrt(2/(mu*omega[0]*sig0))\n",
"#nn=np.ceil(np.log(0.3*L + 1)/np.log(1.3))\n",
"#h = 5*(1.3**(np.arange(nn+1)))\n",
"\n",
"h = np.ones(18)\n",
"x0 = np.array([0])\n",
"#sig = sig0*np.ones((len(h),1)) \n",
"#sig[0:50] = 0.1\n",
"#sig[50:100] = 1\n",
"# Make the mesh\n",
"mesh = simpeg.Mesh.TensorMesh([h],x0)\n",
"sig = np.zeros(mesh.nC) + 1e-8\n",
"sig[mesh.vectorCCx<=0] = 1e-2"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 10.]), array([ 62.83185307]))"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fr,omega"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Make the operators\n",
"G = mesh.nodalGrad\n",
"Av = mesh.aveN2CC\n",
"Li = scipy.sparse.spdiags(1/mesh.hx,0,mesh.nNx,mesh.nNx)\n",
"Mmu = scipy.sparse.spdiags(mesh.hx/mu,0,mesh.nCx,mesh.nCx)\n",
"Msig = scipy.sparse.spdiags(Av.T.dot(mesh.hx*sig.ravel()),0,mesh.nNx,mesh.nNx)\n",
"# The boundaries\n",
"bc_b = np.zeros((mesh.nCx,1))\n",
"bc_b[0] = -1 # Set the top b field to 1\n",
"bc_e = np.zeros((mesh.nNx,1))\n",
"# Make the sparse matrix\n",
"bc = sp.vstack((bc_b,bc_e))\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/Gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:90: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format\n",
" SparseEfficiencyWarning)\n"
]
}
],
"source": [
"b = np.empty((mesh.nCx,len(omega)),dtype=np.complex64)\n",
"e = np.empty((mesh.nNx,len(omega)),dtype=np.complex64)\n",
"# Loop all the frequencies\n",
"for nrOm, om in enumerate(omega):\n",
" # Left hand side\n",
" A = sp.vstack((sp.hstack(( -G.conj().T.dot(Mmu), - Msig)), sp.hstack((1j*om*scipy.sparse.identity(mesh.nCx) , G))))\n",
" #A = A.tocsr\n",
" # Solve the system\n",
" bef = scipy.sparse.linalg.spsolve(A,bc)\n",
" # Sort the output\n",
" b[:,nrOm] = bef[0:mesh.nCx]\n",
" e[:,nrOm] = bef[mesh.nCx::]\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fa8671e1fd0>]"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"# Plot the solution\n",
"z=e[0,:]/(b[0,:]/mu)\n",
"app_res = ((1./(8e-7*np.pi**2))/fr)*np.abs(z)**2\n",
"app_phs = np.arctan(z.imag/z.real)*(180/np.pi)\n",
"ax_res = plt.subplot(2,1,1)\n",
"ax_res.loglog(fr,app_res)\n",
"ax_phs = plt.subplot(2,1,2)\n",
"ax_phs.semilogx(fr,app_phs)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([-5714285.5+0.00050081j], dtype=complex64)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Calculate the impedance\n",
"z = e[0,:]/(b[0,:]/mu)\n",
"z"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 4.13555827e+17])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"app_res"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
}
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
+925
View File
@@ -0,0 +1,925 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import SimPEG as simpeg\n",
"import simpegMT as simpegmt\n",
"from scipy.constants import mu_0\n",
"\n",
"def omega(freq):\n",
" \"\"\"Change frequency to angular frequency, omega\"\"\"\n",
" return 2.*np.pi*freq"
]
},
{
"cell_type": "code",
"execution_count": 115,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Define the mesh\n",
"z = 100.\n",
"hz = [(z,5,-1.5),(z,10),(z,5,1.5)]\n",
"M = simpeg.Mesh.TensorMesh([hz],'C')\n",
"# sig = np.zeros(M.nC) + 1e-8\n",
"conds = [1,1e-2]\n",
"elev = 300\n",
"sig = np.zeros(M.nC) + conds[1]\n",
"sig[np.logical_and(M.gridCC>-200,M.gridCC<0)] = conds[0]\n",
"sig[M.gridCC>elev] = 1e-8\n",
"sig[M.gridCC<-500] = 1e-1\n",
"sig[M.gridCC<-900] = conds[1]"
]
},
{
"cell_type": "code",
"execution_count": 116,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([-2478.125, -1718.75 , -1212.5 , -875. , -650. , -500. ,\n",
" -400. , -300. , -200. , -100. , 0. , 100. ,\n",
" 200. , 300. , 400. , 500. , 650. , 875. ,\n",
" 1212.5 , 1718.75 , 2478.125])"
]
},
"execution_count": 116,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.vectorNx"
]
},
{
"cell_type": "code",
"execution_count": 141,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"\n",
"# Calculate the analytic fields\n",
"freqs = np.logspace(4,-4,33)\n",
"Zana = []\n",
"for freq in freqs:\n",
" Ed, Eu, Hd, Hu = simpegmt.Utils.getEHfields(M,sig,freq,np.array([elev]))\n",
" Zana.append((Ed + Eu)/(Hd + Hu))\n",
"ZanaArr = np.concatenate(Zana)"
]
},
{
"cell_type": "code",
"execution_count": 169,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Calculate the synthetic solution\n",
"Zsyn = [] \n",
"Qex = M.getInterpolationMat(np.array([elev]),'Ex')\n",
"Qfx = M.getInterpolationMat(np.array([elev]),'Fx')\n",
"for freq in freqs:\n",
" e = simpegmt.Utils.get1DEfields(M,sig,freq,sourceAmp=None)\n",
" h = -(M.nodalGrad*e)/(1j*omega(freq)*mu_0)\n",
" Zsyn.append((Qfx*e).conj()/(Qex*h).conj())\n",
"ZsynArr = np.concatenate(Zsyn)"
]
},
{
"cell_type": "code",
"execution_count": 170,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(array([12]),)"
]
},
"execution_count": 170,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.where(freqs==10)"
]
},
{
"cell_type": "code",
"execution_count": 183,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(0.000197894493906+0.000181285727378j)\n",
"[ 1.00551946+0.00563967j] [ 2776.87342066-2515.31936068j]\n"
]
}
],
"source": [
"print ZsynArr[-1]\n",
"print (Qfx*e).conj(), (Qex*h).conj()"
]
},
{
"cell_type": "code",
"execution_count": 171,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def appResPhs(freq,z):\n",
" app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2\n",
" app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)\n",
" return app_res, app_phs\n",
"appAna_r, appAna_p = appResPhs(freqs,ZanaArr)"
]
},
{
"cell_type": "code",
"execution_count": 172,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 99.9982606 100.03796243 99.63440411 99.33759941 106.29388449\n",
" 120.23667968 122.60558712 101.23963102 70.5016044 45.28589223\n",
" 28.43953028 17.6198611 10.84432222 6.83742416 4.62137881\n",
" 3.61487549 3.53761463 4.30735297 5.99370277 8.77556001\n",
" 12.87536967 18.46051811 25.52731592 33.81862185 42.83565417\n",
" 51.95659295 60.60094607 68.35309357 75.00153541 80.50894752\n",
" 84.9530758 88.47009307 91.21391269]\n",
"[ 44.9980325 45.00483104 45.04645859 44.4151918 43.72907936\n",
" 46.72033941 54.60169587 63.94626768 71.19189243 75.42083811\n",
" 77.37907592 77.73640994 76.058035 71.73541117 64.19535734\n",
" 53.36978096 41.16658711 30.81763511 24.08982259 20.77904552\n",
" 19.99288987 20.90312313 22.8822938 25.46014498 28.27778379\n",
" 31.06840674 33.65137133 35.92447275 37.84909319 39.43106457\n",
" 40.70229309 41.70642543 42.48936455]\n"
]
}
],
"source": [
"print appAna_r\n",
"print appAna_p"
]
},
{
"cell_type": "code",
"execution_count": 173,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"appSyn_r, appSyn_p = appResPhs(freqs,ZsynArr)"
]
},
{
"cell_type": "code",
"execution_count": 174,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 148.47846913 175.57907495 170.44187937 154.82139762 151.64437874\n",
" 157.04097616 142.11024881 103.42128503 66.07607566 41.13098488\n",
" 26.10536117 16.64437152 10.50362449 6.72485409 4.58294174\n",
" 3.60236081 3.53787296 4.31623649 6.00810219 8.793518\n",
" 12.8970497 18.48663881 25.55782824 33.85235462 42.87059902\n",
" 51.99052864 60.6320452 68.38024378 75.02434148 80.52754625\n",
" 84.96791113 88.48173456 91.22293906]\n",
"[ 24.90447028 35.50740917 43.45028334 46.88513715 48.17221695\n",
" 52.65990261 61.40691455 69.96355574 75.14671736 77.09349558\n",
" 77.37322373 76.95784679 75.19022669 71.03384758 63.68398281\n",
" 53.00055411 40.91266512 30.67116253 24.02349762 20.75523371\n",
" 19.98819197 20.90730695 22.89080611 25.4705669 28.28863141\n",
" 31.07873966 33.66064252 35.93243359 37.85570553 39.43641995\n",
" 40.70654829 41.709758 42.49194635]\n"
]
}
],
"source": [
"print appSyn_r\n",
"print appSyn_p"
]
},
{
"cell_type": "code",
"execution_count": 175,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAYEAAAEHCAYAAABIsPrhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
"AAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYVOX5//H3TZUOSy8rxQ4aBUFAkGpZJYqxo5iIgvhL\n",
"sKDGEsHZ9YsRUMSvkqJRsHyjIIkaWygqKBG7QREBQUHKwgpSpNf798cMuOKyzOzOzpnyeV0X18Wc\n",
"Oec5H8fduXme55znmLsjIiKZqVzQAUREJDgqAiIiGUxFQEQkg6kIiIhkMBUBEZEMpiIgIpLBVARE\n",
"RDKYioCISAYr0yJgZi3N7HEzm1yW5xERkZIp0yLg7kvcfWBZnkNEREou5iJgZuPNrMDM5h6wPcfM\n",
"FpjZIjO7PX4RRUSkrJSkJzAByCm8wczKA+Mi21sD/czsuNLHExGRshRzEXD3WcD6AzafAix296Xu\n",
"vguYCPQ1sywz+ytwknoHIiLJp0Kc2mkKLC/0egXQ0d3XAdcVd6CZaRlTEZEScHcrbRvxmhgu1Re5\n",
"u1vhP0BeNNuieS9exxy4/8GOT4Wch2qztOdUztTLmYifzVTJmUK/63ERr57ASiC70Otswr2BqJhZ\n",
"LjDT3WdGNs0sYreitkXzXjyPieb44totyTlLcsyhjj9Um7GeM9b9o23jUO3Get5Y94+2jUO1G+t5\n",
"Y90/2jaKa7ck5yzJMdG0UVy7JTlnSY6J5vji2i3JOaM6xsx6AD1K0H7R3D3mP0ALYG6h1xWAryPb\n",
"KwFzgOOibMtLkiHRf4DcoDMop3IqpzIWyunxaKckl4g+B8wGjjaz5WY2wN13A0OAqcCXwCR3nx9D\n",
"m7mR6pbMZgYdIEozgw4QpZlBB4jSzKADRGlm0AGiNDPoAFGYGXSA4phZj8joSXzai1SUwJiZexzH\n",
"t0REMkG8vju1dpCISAaL18RwqRQxMSwiIkWI98SwhoNERFKQhoNERKTUNBwkIpJCNBwkIiIaDhIR\n",
"kdLTcJCISArRcJCIiGg4SERESk9FQEQkg6kIiIhkME0Mi4ikEE0Mi4hI3L47k6MnkGcTgAnALA8F\n",
"XJVERDJIsswJzAP+AiyyPBtmeXZ40IFERDJB0gwHWZ4Z0AG4CrgU+JRw7+BFD/m2ACOKiCSdeA0H\n",
"JUURAPIoNDFseXYYcD4wgHBheJ5wQfhQw0UikskKTQyH0qgI+CRguDuLfvZ+nmUDVxIuCN8BN3vI\n",
"P0hwTBGRpJJmPQEfBgwF/gHc407+z/bLs/LAr4H/Ad4B7vSQf5vQsCIiSSKtioC7mxl1gTuBq4FH\n",
"gVHubPjZ/nlWDfg9cD3wGHCfh/yHRGYWEQla2hWBH1+TDeQC5wH3A4+487OJYcuzpsAIICey/xMe\n",
"8t2JyCwiErS0LQI/buc44F7CE8N5wJPu/OxL3vKsLTAGaAjc4iGfUsaRRUQCl/ZF4Mf36QiMBBoD\n",
"dwEvuPOT0JHLS88l3HNYCtzqIZ9bZqFFRAKWMUUgvA8GnEm4GOwE7nBnxs/2y7OKwGDgbuBpINdD\n",
"vjn+qUVEgpVWzxMws9zIta9FcsfdmQqcDDwEPGHGv8046Sf7hXyXh3wccDzQAJhneXZuGUYXEUko\n",
"M+sRWXQzPu2lQk/g58dQCbgWGAa8QfgegyU/2y/PehFejmIecIOHfEUcIouIBC6tegKxcmenO+OA\n",
"o4BFwMdmPGxGg5/sF/K3gBOBz4D/Wp7dZHmWFIvmiYgkg5TsCfy8DRoQ7hVcATwMPOjOpp/sk2fH\n",
"EO4V1AIGe8g/Ls05RUSClFETw9G3RSvCdxT3Inx56WPu7Nz/fvgqov7AaMLrEQ3XjWYikopUBIpt\n",
"k7bAfYSHi4YBk9zZu//9PMsCRhG+0ex2YKKHfG9RbYmIJCMVgajapifhL/sKhJekmFb4HgPLsy7A\n",
"WMAI31vwdlnkEBGJNxWBqNvHgAsJDw+tIHyPwUf738+zcsAlhHsOnwO3e8gXlFUeEZF4UBGI+TxU\n",
"JLwcdQiYDdzlzlf738+zysAQ4A7Cq5nmesgLyjqXiEhJpEQRMLNqwJ+BHYQfGvNsEfsk9EHzZlQF\n",
"bgRuIfxlP8Kd/fcPROYLhhFetvoh4EEP+dZE5RMRiUaqFIErgXXu/pqZTXT3y4rYJ6FF4MfzUpfw\n",
"pPBA4EVgtDsL97+fZ0cAfwS6AMOBpz3kexKdU0SkKIHdLGZm482swMzmHrA9x8wWmNkiM7s9srkp\n",
"sDzy96T6AnXne3duI3wF0bfALDP+YUZ7AA/51x7ySwnPJ1wNzLE8uzAyhyAikhZK8oU2gfCllfuZ\n",
"WXlgXGR7a6CfmR1HeCI2uxTnKnORYnAP0BKYBbxgxnQzepthkUdZdiN8ddEdhO88vkDFQETSQYmG\n",
"g8ysBfCKu58Qed2Z8EOPcyKv74js+jDh4rAdmOXuzxXR1r4Hze+z/4HzQYisS3Q54aGiTYRXLn3J\n",
"nb2Rm836EH6ITUXCuV/SPQYiUtYKPWB+n+AeNF9EEbgIOMvdB0Ve9wc6uvv1UbQVyJzAoZhRDuhL\n",
"uAdQk/CzCp51Z9sBxaAC4WLwLxUDEUmUZFtArlSzy4daSjoI7ux150WgI/BbwnMDy8wYTa638JC/\n",
"SvipZ8MJX030qeXZrzRMJCJlKd5LScfrC2slP479E/l71Ms2u3tukENAxYk8y+Atd84BOhH+zD4y\n",
"41VyPYdcfw1oz0+LQa8AI4tIGnP3me6eG6/24jUcVAFYCPQG8oEPgX7uPj+KtpJyOKg4kXsNLiN8\n",
"c1kNwquTTiDXNgDnE76/YCbhpSjWBJVTRNJXkJeIPkf4jtujzWy5mQ1w992EvxCnAl8Ck6IpAIXa\n",
"TLrhoOK4s9Wd8YSfdPZroB3wDbn+GLm+BGgDrAW+sDy7KjKHICJSanqyWJIyoyHhG8+uI/wks+Hk\n",
"2h7gUWAL4WcYLCymCRGRqKXEHcNRBfjxEtFALw2Nl8glplcTnh/4mMobcrmzTjfgbsKXy97nId8R\n",
"ZEYRSV2FLhUN7hLReEqXnsCBzKgCDCZ8g9lMOo39Czk330h4qGiwh1K/4IlIcNKqJ5CORWAfM6oT\n",
"ni+5GXidq7vM5vDZw4A3gd97yNcGGlBEUlKy3SdQKqk2MRwLdza7M5LwGkVLGf/uHxldMJ0tdfcQ\n",
"Xo/olIAjikgK0cRwiousXvp7YBCdH3iHM3/fFeNGD/18mW0RkYPRcFCKi1xN9AhNPmzL1adVocLO\n",
"pwg/+F5LT4jIIakIpIHIoy9/Q7WCBxjUcQO1vv0Co7+HfHPQ2UQkuWlOIA1ElqR4ki0NO/HIgu9Z\n",
"cP7J7Kn4geVZ86CziUhy0pxAmgo/A3nvcE594EZ6hvZQcft5HvL/BJ1LRJKThoPSlBldOOblf3JB\n",
"/5p4uZv9vg1/DTqTiCQfFYE0ZkYtms1+mosuP5u95SeR9c1Ver6xiBSmOYE05s5GX35qX559+Xf8\n",
"kH0pa4+ZZ0Ob1wo6l4gET3MCGcayFh9J72H/odF/K7Gm9fE+8cX8oDOJSPDSqicgB+frjlzMtPtb\n",
"kN++gAZzF9nZN7QOOpOIpA8VgRTgG7O3s6BvG1Z2/JhjXp1j3e7tFnQmEUkPGg5KMXbxJRPJfv9i\n",
"3r3tcv9gyKSg84hIMDQclKF88vOXsfqkhznt3mft5MdvCTqPiKS2pCgCujooNv7sy0NZ0/pOev9h\n",
"tJ34zP9Glp8QkQygq4NkPxt88gBqL32M1x95jS8uv8id3UFnEpHE0M1iAoDd1Oo8qqz7B6/9aQ5z\n",
"r+jljhafE8kAKgKyn91R+zTK7Z3G64+s4rPfdHVH9xKIpDlNDMt+PnLDLCps68I5Q+rR7m9zzDgy\n",
"6EwikhpUBNKE37PrUypv7kjO0Ip0+POHZhwfdCYRSX4qAmnEQz6fSltO4Yzb9tLpoXfN6BB0JhFJ\n",
"bioCacZDvohKWzrQ665tdBn1lhndg84kIskrKYqA7hOILw/5EiptPYXu92yk2z2vmXFO0JlEJD50\n",
"n4BEzfKsCbuqzOa9m+syI+8a31v++aAziUh86OogOSQPeT4Vt3Wk84OrOeP2J6zcnquDziQiyUVF\n",
"IM15yAuouK0zHf60nLNveMjK7bkp6Ewikjw0HJQhLM+y2F15Bp9fcTiv/vUh9la8x51g/+eLSInp\n",
"jmGJmeVZLXZXeoP5F7bihaefxCvcqkIgkppUBKRELM9qsKfiFL765dFMnvhP9lb6rTt7g84lIrHR\n",
"xLCUiId8E+V3ncnRr8zl8vMuoPz28WaUDzqXiASjTIuAmbU0s8fNbHJZnkdi4yHfQvnd59DyzY+5\n",
"8sxzqLj5OTMqBp1LRBKvTIuAuy9x94FleQ4pGQ/5dsrv7kv2e29zVc+eVN7wTzMqB51LRBIrqiJg\n",
"ZuPNrMDM5h6wPcfMFpjZIjO7vWwiSlnxkO+i/O7LaPzf17mm66lUXfOqGVWCziUiiRNtT2ACkFN4\n",
"g5mVB8ZFtrcG+pnZcWZ2pZmNNbMm8Y0qZcFDvodyewZQb/4kruncjuqrpppRPehcIpIYURUBd58F\n",
"rD9g8ynAYndf6u67gIlAX3d/xt2Hunu+mWWZ2V+Bk9RTSF4e8r2U2zuEOkueYFDH46n17Ztm1Ao6\n",
"l4iUvQqlOLYpsLzQ6xVAx8I7uPs64LpDNXTAYkgz3X1mKXJJCXjI3fLsdmqs3MygjjfzxH/eNjuy\n",
"lzvrgs4mIuGF44Ae8W63NEUgbjcYuHtuvNqSkvOQO3CP5ZbbzKBOIca/8x+z1t3dWRN0NpFMF/nH\n",
"8cx9r80sFI92S3N10Eogu9DrbMK9gZhpKenk4rl7H6TKutsYeGozGv33PTMaB51JRMICW0razFoA\n",
"r7j7CZHXFYCFQG8gH/gQ6Ofu82MKoDuGk5blWX92VP8zT0/fwMpOXdx/MvwnIgFK6B3DZvYcMBs4\n",
"2syWm9kAd98NDAGmAl8Ck2ItAIXaV08gCXnI/4/Km6/iN73r0Hzme2a0CDqTSKbTQ2Uk4SzPcthV\n",
"ZTKT/rmFxWef6s43QWcSyXRaO0gSxkM+hYrbfsmlF1bj2BfeN+OooDOJSHwkRRHQcFDy85C/TcVt\n",
"p3Nh/8r84pn3zTgu6EwimUjDQRIoy7NfsLPq20wdA59cd5o7XwSdSSQTaThIAuEh/5xKWztx1q17\n",
"6DzmXTNOCjqTiJRcUhQBDQelFg/5QiptaU+P3G10z5tlxslBZxLJFBoOkqRhedaE7TXf55Nrs3jz\n",
"j719T8UPgs4kkik0HCSB85Dnc9gP7Wj3+GrOvmGGVdzWNehMIhKbpCgCGg5KXR7ytVTZ0J42k5dw\n",
"7uDpVmlzz6AziaQzDQdJUrI8q86W+u+wtFtrXn30XN9ad3rQmUTSmYaDJKl4yDdTbU0Xmv/nv/zq\n",
"169ajVXnBJ1JRA5NRUDixkO+jeoF3Wk0ZzYXXPmi1VnSN+hMIlI8FQGJKw/5Tmrmn0HdhTP41a8n\n",
"W70FFwWdSUQOLimKgCaG04uHfDe1VvSh1rLXueDKZ63h55cHnUkkXWhiWFKG5ZmxvuWz7KhxES//\n",
"baCvPOWpoDOJpAtNDEvS85A7dZZcTsVtT/Or3zxuzWcNCjqTiPyUegKSEDak9TjK77iOVx+9yb8+\n",
"fVzQeURSnXoCklJ83JdD2FX1Qc67+iE7+tWbg84jImFJUQQ0MZwZ/M9zb2N71r2cO3i0tf7HH4LO\n",
"I5KKNDEsKc8Gt7+NmsvvY+qDo/zzK1QMREpAw0GSsvzRj0ezsfn15Nx0m7V74pGg84hkMvUEJDA2\n",
"qOOVZC2ewIx7/o+PfjfAnWB/GEVSSLy+O1UEJFA2qNP5ZC36B7P+8BLv3XKxCoFIdFQEJG3YNV1O\n",
"p97Cf/Pe0DeYdVcfd/YGnUkk2akISFqxa7p2pu5XM/n42g+YMaKnO3uCziSSzFQEJO3YoE4nUuvb\n",
"95l7+ZdMG9PRnd1BZxJJVml1dZDuExAA/9v7n7H+iF/QZvKx9Lnuc6u7qFLQmUSSje4TkLRngzpl\n",
"U3XNFyzr+j1v3NfGNzXZFnQmkWSj4SBJa3ZthwZU3vQlBSfsYNYfWvuqthuDziSSTFQEJO3Z//tF\n",
"LYx5bG5Yjdm3tPXFOUuDziSSLNJqTkCkKP6XzzdScWsrKm9aRq/h8+0Xf+8QdCaRdKMiIEnN/3fx\n",
"Tpp90JYdNWbSI/c96/Dn84POJJJONBwkKcOu+OWjNP5kIO8NvdHfvU3PJJCMpuEgyTj+91cHs/zU\n",
"YXS5/2HrmXt/0HlE0oF6ApJy7OJLL6flW8/w2a9fYNqYS7TekGSilLk6yMz6An2AmsAT7j79gPdV\n",
"BCRmdlG/bjR/5w0Wn/0Jn/fv6kt6aJkJySgpUwT2n8isNvCAuw88YLuKgJSIXXzZMTT+9GMKTviO\n",
"T649wReftTXoTCKJkvA5ATMbb2YFZjb3gO05ZrbAzBaZ2e3FNDEM0GSexI1PnriQguNbUmdJNbqO\n",
"WmZtJzQJOpNIqollYngCkFN4g5mVJ/zFngO0BvqZ2XFmdqWZjTWzJhY2Cvi3u8+JW3IRwCe9sJYd\n",
"NZpjuws4dcwS6zL6vKAziaSSqIuAu88C1h+w+RRgsbsvdfddwESgr7s/4+5D3T0fuB7oDVxkZoPj\n",
"FVxkH5/w9g5azDqetcc8Q5f7X7Kcmx81Q0OMIlGoUMrjmwLLC71eAXQsvIO7Pww8XFwjB6yIN9Pd\n",
"Z5Yyl2QYD7kDA+38Aa/T+h/PUWtpdzvpvFN9zlXrgs4mEg+RlZZ7xL3dWCaGzawF8Iq7nxB5fSGQ\n",
"4+6DIq/7Ax3d/foY2tTEsMSV/eo3jan71ftU2tKAuf3O81l3Tj/0USKpJVluFlsJZBd6nU24NxAT\n",
"PU9A4slffGoV2e+3YO0xL9F57FTLuXlM0JlE4iXQ5wkU0ROoACwkPOafD3wI9HP3+TG0qZ6AlBk7\n",
"b+AVHPXvJ1l+6pcsPquLfzpwc9CZROIhiEtEnwNmA0eb2XIzG+Duu4EhwFTgS2BSLAWgUNvqCUiZ\n",
"8Jcf/ztfnXsktZbX55Q/rbZuI7oGnUmkNPRkMZESsEsuLof5Cxw+61zm9hvpUx66K+hMIqWRcncM\n",
"HzSAioAkkJ07+FqOen0cG1puYGm3Qf7WiH8FnUmkJNKqCAB56NJQSRA7f0B1Km59mlZvns+KTgtZ\n",
"2aGfzwzpRkZJCYUuFQ2lTRFQT0CCYGff1Jw6Xz9P9uwOfH3GOxScdKnPuqMg6Fwi0UirnoCKgATJ\n",
"coZ2odGc/6PegsNZfPazbMwe6DPydgSdS6Q4aVUE0HCQJAE754b+ZL/7CJU3VeHrM0aytX6ez8jV\n",
"swokqWg4SKQMWc9co/rqXFpNv4NtWVtYecoDbK03Rj0DSTZp1RNQEZBkYz1Dlam17FEaf3IhNVZX\n",
"ZWWHeaxvNY7qqx/35yfvDTqfSFoVATQcJEnMTr+jF3WW3knjT7pRcVt5Vnb4kI2HP0Cnh1+MLFwn\n",
"kjAaDhIJiPXMNSpsvYisb26h6Qft2VN5N6vavs0Pze71KWPfCTqfZJa06gmoCEiqsTbPl6PZB4Op\n",
"u+i3NHu/DVvrbmPdkQvY0uAtdlX5u78+TvcdSJlSERBJEtb+0SpkLRpI9dV9qLX8JOp/2YA9lfaw\n",
"7qgVbGryIdvqvEClzS/5i09pclniJq2KAJoTkDRibZ4vR4MvzqTG6oupsbILWYtbUCO/Mt8fvY6N\n",
"h89jW9Ysdh82hWrfvatJZomV5gREUpB1Gd2KmvlXUH3VmdTIP5raS+tSZX151rfcxKYmy9ha73N2\n",
"1HiHcntf9Zcfi/mZHJJ50qonoCIgmci6jG5Fte/6UPX7blT77gRqrmxGna+rsaPmbjY2X8vmhovY\n",
"XudjdlabBrzlrz+yM+jMkjxUBETSkLV5viIN5vWg6vdnUnVtR6qtPpLay+pTfXUFNjTfyqamK9lS\n",
"bx7ba7/L3kqv+usPLwg6swRDRUAkg9ipYxpRfXUfqqzrQdW1J1Ij/3Cyvq7J3vLOhhbr2dRkMVvr\n",
"vcfO6i+xvc47Wu4i/akIiGQ465lr7KnQnmrf9aHa2q7UWNmGOl/Xp/Kmcqw7aj0/NF0UKQwvUyP/\n",
"bU1Cp5e0KgLo6iCRuLHu/9OGamsuoOqa7tTIb03W1w2otKkc647cwA/NFrKl4TR2Vnvap4z9Ouis\n",
"EjtdHSQiMdtfGKoV9Kb2tydQb34WW+vvYN0R37CpyTtsr/0sO2rN0jBS6kirnoCKgEhi2aljDqNa\n",
"wYVUL7iAmis6UP/LppTbA2uOy+eHZh+yre4LVNg+yV9+bHfQWaVoKgIiEjfWM9ewPadRraAfNfO7\n",
"U/erVlRdU4mCX6xgQ4u32Jb1qE8Z+17QOeVHKgIiUqas190nU33VddReegaNPjucndV2sab1An5o\n",
"9i92Vv+LTx2zKuiMmUxFQEQSxnrmlqfiloupkf9rshZ1ov78Oqw7aiPrjviIzY2fZFvWs5pPSCwV\n",
"AREJjHUfUY9qBf+PWst/RcPPj6fC9nKsavcFG5o/xZ5Kf/EpY7cHnTHdqQiISFIIP2dh+7nU+nYI\n",
"Db7oQs2VVcg/+RvWt5rM9tr3+/RR64LOmI7Sqgig+wRE0ob1vqsjNVfcSr35vam/oA6rT1zFuqNe\n",
"YUuD0T59pO5NKCXdJyAiKcN6ho6gRv5t1P3qXBrNacya1mtY03oyW+vl+fRR3wWdL5WlVU9ARUAk\n",
"/Vn3expRc8VwGsy7mAbz6rOq3VLWHvMkWxqM9hl524LOl2pUBEQkZVmvu4+j1rJcGn52NrWWV2fF\n",
"KfNYf8QjbKv7N11lFB0VARFJC9b7ru7U+WYYTT7uRoUd5chv/wEbWozwKQ9OCTpbMlMREJG0Yj1z\n",
"jYpbLiVr8S1kz27HpiabWdV2Elvr/8Gnj1obdL5koyIgImnLeoaqUHXN3TT6/Grqz6vPis5f8v1R\n",
"f6TTw895KOAvrSShIiAiGcF6DW9P3UUjyZ7dnZ3VdpHf/mU2Nr/d3xzxbdDZgqQiICIZxXrmlqfK\n",
"97fQYN4QGn+azcoO37D22LFsq/unTJxMTokiYGbHAjcCdYGp7v5EEfuoCIhITKzX3cdQe8lImr1/\n",
"DpizvPM/2dxoaCbde5ASRWD/SczKARPd/ZIi3lMREJESsZ65xmHrbqbRZ7fQcG4jlnWdw9qjb/dp\n",
"D0wPOltZS2gRMLPxQB/gO3c/odD2HOAhoDzwuLuPKuLYc4HfAn9z9xeKeF9FQERKzXr/oTP1Fo6h\n",
"+TudWH/E96xq+yfK7RmRrg/GSXQROA3YDDy9rwiYWXlgIXA6sBL4COgHtAfaAfe7e36hNv7l7n3L\n",
"6j9ERATAet1dixr599P0/Ss47IdKLOv6b9a3vMnfvPeboLPFU8KHg8ysBfBKoSLQmfACRjmR13cA\n",
"uPvIQsd0By4ADgPmu/tDRbS7bwG5fbSQnIiUmhnGmTdfRcO5w2j6UUuWd17Amta3+tQxrwedrSQK\n",
"LRy3T2IXkCuiCFwEnOXugyKv+wMd3f36mAKoJyAiZcx6DW9D3YXjaDmzG+uOXEt+u1Fsqzc2la8q\n",
"itd3Z4VSHBu3D8/MclEPQETKiL/1P/OAntbr7jrUXPEQR069j70VRljf5U/xQ9NbfcY9W4LOGK0i\n",
"egSla68UPYFOQG6h4aA7gb1FTQ4fol31BEQkocL3HKy9i2Yf3kTN5bVY2n0aG1v81qePXBJ0tmgl\n",
"w5xABcITw72BfOBDoJ+7z48pgIqAiATIzrz1Ehp9fh9NPmrJsi5zWXvcUJ92/1tB5zqUeH13lovy\n",
"ZM8Bs4GjzWy5mQ1w993AEGAq8CUwKdYCUKj93EgXR0QkoXzaA8/709OO4KPftWdP5S20Hf+GXdVj\n",
"meUMvSLobEUxsx6RIfT4tKdlI0REfmS9hjemzpK/0vKtX7Ixez0rO/wxGSeRU+qO4WID6BnDIpKE\n",
"rOvI6tRd9L+0mNmfXVV3sbzLn7C9dwV985meMSwikkDWM7c81VffS/bsIVTaXJGl3Z9la93rfdoD\n",
"mwPNlU49ARUBEUl24XWK1t9I0w+HUfvb2izt8TobD7/Wp49cHUiedCoCaDhIRFKInXXLZTT+72ga\n",
"ft6MJb1msL7V1T59ZEKeb6DhIBGRJGFn3H4WjeY8QtOPjmRpj/f5/qhrfPqoEl0lGfO506knoCIg\n",
"IqnMzrijKw2++CvZ77bm225zWNN6kL9x3ydlek4VARGR5GK9h7WjwReP0fztdqzoPJ/v2gz2aff/\n",
"p0zOlU5FAM0JiEgasd53HUu9BY/Tcsap5Lf/htUnDvFp90+JS9uaExARSQ3We1gLshY9Tqs3e7H6\n",
"xGWsavc7n3b/a3FpO516AioCIpLOrNewpmR9PZ5Wb5zBdycsJ7/d9T7tgZdL1aaKgIhIarFewxtH\n",
"isGZfNcmn1Vtb/CpD75YorbSqQigOQERySDWa3hD6nwznlbTc1jbehX5J9/oU8f8M6pjNScgIpIe\n",
"rNfw+tT55glavdGH748pIP/koT5l7KSojk2nnoCKgIhkMus1vC61lz7OEdPP4/ujCshvf5NPGft8\n",
"sceoCIiIpBfrNiKLeguf4Ihp57H22FXkd7jBpz7wQpH7qgiIiKQn6z6iHvUWjOeIaX34rk0++Sf/\n",
"7sCriVQERETSnPUaXp+sr5+k1Rs5FBy/glUn/3bffQZpVQTQ1UEiIgcVvrR08QRavXEmnzZbzUcF\n",
"b/LDqv5pUwTUExAROTTrNawpdRdPoNUbpzP6e1MREBHJQNZrWDYz7l2mIiAikqHi9d1ZLh5hREQk\n",
"NakIiIhkMBUBEZEMpiIgIpLBKgQdAMDMctF9AiIih1RoFdH4tKerg0REUo+uDhIRkVJTERARyWAq\n",
"AiIiGUxFQEQkg6kIiIhkMBUBEZEMVuZFwMyqmdlHZtanrM8lIiKxSURP4DZgUgLOU6YiN2gkPeWM\n",
"L+WMr1TImQoZ4ymqImBm482swMzmHrA9x8wWmNkiM7u9iOPOAL4E1sQnbqB6BB0gSj2CDhClHkEH\n",
"iFKPoANEqUfQAaLUI+gAUegRdIBEirYnMAHIKbzBzMoD4yLbWwP9zOw4M7vSzMaaWROgO9AJuBwY\n",
"ZGZR3d1WVCUurjqXpHKXttof7PhUyHmoNmM9Zzz+5aScJd8/2jaS7WfzYG0kW85U+V0vqaiKgLvP\n",
"AtYfsPm/p0S+AAADtklEQVQUYLG7L3X3XcBEoK+7P+PuQ909392HuftQ4FngMY9+jYoeUW6L5r14\n",
"HhPN8cW1W5JzluSYQx1/qDZjPWes+0fbxqHajfW8se4fbRuHajfW88a6f7RtFNduSc5ZkmOiaaO4\n",
"dktyzpIcE83xxbVbknOW5JhSi3rtIDNrAbzi7idEXl8EnOXugyKv+wMd3f36mAKEHzQvIiIxisfa\n",
"QaVZRTQuX95aPE5EJDiluTpoJZBd6HU2sKJ0cUREJJFKUwQ+Bo4ysxZmVgm4FHg5PrFERCQRor1E\n",
"9DlgNnC0mS03swHuvhsYAkwlfBnoJHefX3ZRRUQk3gJ/qIyIiARHaweJiGSwpCwCqbDekJkda2Z/\n",
"MbPnzeyaoPMcjJn1NbPHzGxi5A7upGRmLc3scTObHHSW4kR+Np+KfKaXB52nOCn0mabKz2hK/M5D\n",
"bN+hSTkcZGZ5wCZgvru/FnSe4phZOWCiu18SdJbimFlt4AF3Hxh0luKY2WR3vzjoHAdjZlcC69z9\n",
"NTOb6O6XBZ3pUJL9M90nhX5Gk/53Ppbv0IT0BGJZeyjI9YZiXSPJzM4FXiN8t3TS5owYRniZj4Qp\n",
"6ZpTiRZjzqbA8sjf9yRhvsDEIWdCfkZLkzORv/MlzRnzd6i7l/kf4DSgLTC30LbywGKgBVARmAMc\n",
"B4wAxhK+6uglIr2VZMt5wHH/SlTGEnyeBowCeicyY0k/T2ByMucE+gN9Ivs8l4T5roz8/jRJ9Gda\n",
"0pyJ/hkt7ecZ2b/Mf+dL8XnG9B1amjuGo+busyLLThS2f+0hADPbt/bQsMjr3wBrPPJfnmw5zawB\n",
"cAFwGDAjURkhtpzA6UBvoKaZHenujyZjTjMrAP4InGRmt7v7qGTMCTwMjIuMtSbkvpgYf39GAs9E\n",
"tmWRwM+0FDlvIIE/o6XI2Z0E/s6XNCfhHlXU36EJKQIHUbhbDeG7jTvue+HuTyU8UdGKzOnubwNv\n",
"BxOpSAfLeT3wSDCRinSwnOuA64KJVKSD5dwKXB1MpJ8o9vcHIEk+02hyPky4uAYpmpzJ8Dt/yJz7\n",
"RPsdGuTVQck3I1005Ywv5YyPZM+3j3LGV9xzBlkEUmXtIeWML+WMj2TPt49yxlfccwZZBFJl7SHl\n",
"jC/ljI9kz7ePcsZX/HMmaDb+OSAf2EF4PGtAZPvZwELCs913JiKLcipnquVM9nzKmdo5k/JmMRER\n",
"SYykXDZCREQSQ0VARCSDqQiIiGQwFQERkQymIiAiksFUBEREMpiKgIhIBlMREBHJYP8f64VSM4TA\n",
"JcIAAAAASUVORK5CYII=\n"
],
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc5b9632550>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"loglog(freqs,np.abs(ZanaArr),freqs,np.abs(ZsynArr))\n",
"gca().invert_xaxis()"
]
},
{
"cell_type": "code",
"execution_count": 176,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
"AAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FPW5x/HPw82gQCIi8agoBLmooGgVELRGrCKgpbXF\n",
"a7Eq4pF6beudc7rZnqL12J56q9hDxapHUTx6qhUBqxiqUpRqsYpXJFSwioqEgiEI+Jw/ZgMbSMhu\n",
"djezm/m+X695yc7O/ObJmjzz29/85hlzd0REJDrahB2AiIi0LCV+EZGIUeIXEYkYJX4RkYhR4hcR\n",
"iRglfhGRiFHiFxGJGCV+EZGIyXriN7P+ZjbVzGaa2YRsty8iIpmxXN25a2ZtgIfc/bScHEBERJol\n",
"pR6/mU03s1Vm9vp2608ys7fN7D0zuyZp/SnALOCh7IYrIiKZSqnHb2bHAOuB+9x9YGJdW+Ad4BvA\n",
"h8Ai4Ex3fytpv8fdfWwuAhcRkeZpl8pG7v68mfXcbvVgYKm7Lwcws4eAsWbWHTgVKAKey1qkIiKS\n",
"FSkl/kbsA6xIer0SGOLu84H5Te1sZioLKiLSDO5umeyfyayejBO3u1vdAsSTXze2Lp33M92+oX3S\n",
"jTMbx2wtcaYaU77F2dj+hRBna/gbUpw77JOxTBL/h0CPpNc9CHr9zVWZ4rp03s90+1Tb2Fm7zTlm\n",
"c/ZJpY2dtducYzZnn1Ta2Fm7zTlmc/ZJZf+dtducYzZnn6b2b6rNdI+Z7vapttFUu+keN93tU22j\n",
"qXbTPW6622eHu6e0AD2B15NetwPeT6zvACwGDkyjPQcqgPJU9wljASrCjqG1xFkIMSpOxZmvC1Ce\n",
"yJmeaVupTuecASwA+prZCjM7z903A5cAc4E3gYc9aUZPiiedCnevTGefEFSGHUCKKsMOIAWVYQeQ\n",
"osqwA0hRZdgBpKgy7ABSVBl2ADvj7pXuXpGNtixxJmlxZuaepfEqEZGoyEbuVK0eEZGIyWQ6Z8bM\n",
"rAKoLIDhHhGRUJlZOcE4f+ZtaahHRKRwaKhHRETSpsQvIhIxGuMXESkAGuMXEYkojfGLiEjalPhF\n",
"RCJGiV9EJGJ0cVdEpADo4q6ISETp4q6IiKRNiV9EJGKU+EVEIkaJX0QkYjSrR0SkAGhWj4hIRGlW\n",
"j4iIpE2JX0QkYpT4RUQiRolfRCRilPhFRCJG0zlFRAqApnOKiESUpnOKiEjalPhFRCJGiV9EJGKU\n",
"+EVEIkaJX0QkYpT4RUQiRolfRCRilPhFRCJGd+6KiBQA3bkrIhJRunNXRETSpsQvIhIxSvwiIhGj\n",
"xC8iEjFK/CIiEaPELyISMUr8IiIRo8QvIhIxSvwiIhGjxC8iEjFK/CIiEaMibSIiBaDVF2mzuI0B\n",
"XvSYVyetKwGGe8xntVSMIiL5pjUXaXsRmJJI9nVJf0pifT0WtzF12yWtK0mcPEREZDt5mfgTPf3J\n",
"BMl/NHA/MNljXm1x62JxK07aPOWThIiI5Gnit7h1TiT/m4FZwJNJwz7fAq5L2vx4YDlB8u8J3EDi\n",
"JNFiAYuIFJC8S/wWt87AKxa37sBVQC/gkLoevcf8Po/5tUm7vAQ8TnCSqAK2AOMbaDfrQ0IaZhKR\n",
"QhR64k8kyrMsbnsCeMzXAUcDMYKe+3K2DfuUbL+/x3wl8AnbThJtgUeT2r/T4vZ1cnPdQMNMIlJw\n",
"Qk38SYnyUOpPUzqSpOGapDH/4Ttpo+4kcT0wOSlx3wQsTmrjTxa3IXX7AAMsbkVJ7bVlJwnd4vat\n",
"xDYAa4GDgBsTw0xTgKcS60VE8lKo0zmp4E4yHI9Pd+qnxW0C8Fugl8d8ucVtJjDJY7468f6HwNeA\n",
"WoJEfjzwAnBl4uLyfcDFiW8mWNyOBj4ClgKDgArgVI+5W9w6AEUe839qiqqIZENrmM55c6YXYT3m\n",
"s7Zvw2Ne3UjSLwEOJxgSusriVuIxP60u6SfsC6xKurjcD5iS9O3jnLqkn/AG8KNEmxcC53ls69n0\n",
"cODBxL81LCQieSHsxH9VQ+P2udDAkFCD1w085p7orZew7brBlQ3F2VSbHvOFwCmJf1cD7wNPJw0L\n",
"afaRiLS4sBN/oxdtc2A4zb9u0FicTbaZ1PsHuJ3gW0EVWfi2IyLSHKGXbMjHce5cjccnnVBuBuLA\n",
"vR7zeZnGKyLRkY0x/tATfygHzyKzAaOhx2XQuQjW1cKK29zfeGqH7eJWwqpd7ufeYR2p6daOw97s\n",
"wvFvt6PTlq+r5y8iqcpG7gy1OmehC5L+UbfCtAO2rZ3Y22wAOyT/RbtfzrwxB7Hh/jIA/gq8PX4Z\n",
"I2ZdTtD7FxFpEerxZ8Bs1ByYPXLHd0bPcX9qVLANuwFPwpWHwy+6NLQtFbN/RzCTqDKX8YpI4WsN\n",
"0zkLlhntYN99Gn63U8ekFxuA/4DP3tvJtquAz7IboYhIw3Iy1GNmY4ExQBfgbnf/Yy6Ok0s7G7s3\n",
"4yrgR9ClqOG912+o+5c7XwHzzFY1ktjXb1BPX0RaUk56/O7+uLtfCFwEnJ6LY+TStrH72SNh5rHB\n",
"f4+6NVgPBDdtHQdzz4aJS+vvfcH78MHtO7a64rZUtrW4PWlxO2G7dSr8JiJZk/IYv5lNJ+jFf+Lu\n",
"A5PWnwTcQlAc7bfuflPSe78A/sfdFzfQnsNJcxudBZPibJl0tk19u6bH7uu3ud+lwZDN+g3wwe07\n",
"j7Pxbc3YhyN+fSmDf3wIvxvWgZpu7dj1s818f8EGSjeO1+wfEcnK9VF3T2kBjgEOA15PWteWoEZN\n",
"T6A9sBg4EDCC4mjH76Q9B3e44D04eHT99w4eHax337bsuF062+5sO/BdwIeD/xB8Bly/of52dcu4\n",
"ylQ/r+Ys4AfAk29QdFkto3/glLzvjP6B0/F77zf0s2vRoiV6S5C2M2sj5TF+d3/ezHput3owsNTd\n",
"lwOY2UPAWOAbBMXNupjZAe7+m8ZbnnYAnHujGbsDX7rzSNArT54iWbfdOVPM2IWg5v5moLbxbU+f\n",
"bMZHwFfB9v2vb3i70ZcSlFL4FbAImAPv9KCBO3qTx+5zwZ2lZnespHb2wSx8F67oDVNfhQ2HlSXi\n",
"bPCbhIhIOjK9uLsPsCLp9UpgiLtfSlCeoAkVif++1RsePhdOXwY8EgzFNKR0X+D7BN802gI1jW9b\n",
"dggwPbFdG+jXq+HtOnV05x2CkxgAZm9/ChO3m5/f2Nh9tnUuoqgaht4K0/4MX/stzJsCtfVmColI\n",
"RJhZOT2YyMdF7dnU/fBstJlp4s/wJoCKxH9fft799KSx83W1DW+/5C/ufCt5jdm6OQ1v+9oL7oza\n",
"tt2rc4AGxu137MW7v/GU2QCCXnbTY/dZtetnmymfnEj2JbC6P4yYDJWrN+X82CLSYixuY5hatiur\n",
"+k7Yet2x9N27mbSsJrk0jLtX2l6992R0u2nMfamY2t0zP3iaY0s9qT/GPxSYk/T6OuCaVMepgnHz\n",
"CUtTG4/fcbt0tk2nzVDH78bsHgvG9JPiPObwTzmv471hx6ZFi5amFyoYQ2nZODhpDoyrhJPmUFo2\n",
"jgrG1NuutGwcY/tWU7Qm+DsvWuOM7VtNadm4rdsc9HAb8G4w5o+ULHNO/5bTkmP8jfgL0Ccx9v8P\n",
"gqmbZ6a+++g5DfWk0+lxp7ptqL34NPiTn8fNBiyqH+fKuey/4dGm9xaR0E0t25Wh7aYxd0YxtSUE\n",
"Q7dDhjK1bCKxpO1W976Al35ezIjJsOAqOLYCvvh2Mas2/tqMayh9bV++dUN33jxtDey1EgzWNDJi\n",
"na6Uz2IwgyC5byQY1z8vsX4U8A7B7J7r0mjPCcZ6ysM+Q2vRokXLzpZUe/HuDpw0h6I1zpiLnJ7P\n",
"BDPzitY4u577d04fuwz8FfDv0nXUnznneKekyqnA6f6aM/ABhwteBz+STv/Yn/J/77y1zfZ/cA44\n",
"Iis9/vA+yCwEH7WFCvZk3Hd+CL532LFo0RKlZWfDMuCdwfeggl2o4GfBicGdPV8PEnpJVbBP24ve\n",
"Z/QPbgY/MtgncYIY/YNgm7oTBKNm7/z4eKY/j2r1FJYuGCcBs80oDjsYkUJncRtje/UeZzZqjtlp\n",
"lWaj5thevcftcKf8qr4TmPtSMSdfBCVVwYSLuX8u5rQ1D9OuZhUwHvgSqKbN2lqKquHIqXBLFQy7\n",
"ORju2fL3d33Wr69yZ5E76yh9925GDlnLvClQ3TOY0DFy8FpK35m+Q6CTltWwcPNEas9qZDJLmkI7\n",
"g2qop5mfmxv4HeDzwHcJOx4tWgp5aeoCKxWMooLOW3vx5w+r34vffeyr9HyubdNt9ql30TbRdmL4\n",
"aNTsoP1RsxsbPnJ3gPJEzvRMf+5Q6/G7e0WYxy9E7rgZl1O6+ElWDbrX7JD7YZ9LUyltISLbWdV3\n",
"AnNnbLvAespEeHpmMaum/RJ4hCDZvgPrgl78x4PgsQeCXvy8KbDmy1X+efmWem1OWlbD1LKJ1J51\n",
"PnTqSO36DSzcMp1Jy2qSN/OYz0pc7H0klVDdvRKoNLNYU9s2RfX4C5DFrR1ftZnPr6Z3Zd1re8F/\n",
"JT0LeOJS+PPlSv4SVanOjwewQ4e9xWe39aemG1zRC+beDO+eAqunvOx+35Ct2+3Ve1wwU+elbTN1\n",
"Rg5ey8ItE/3j91NK3Fn7+fQErmjymG+2uB3NukefgSf61393axkKJX6Jpp1MpzRsINDZY74AAPti\n",
"Ax3WwaB7to3Hf1EKfPZ5vTZT7MUXjNDG1jSrJ/PPcODgv1G02rHN2272KlrjDBz8t7Bj06IlrKXe\n",
"bJm9XnHGfScxW2bsPCo4gQpO37ptiuPx+bRkI3eG2uM3swqg0l0PImmW93f7hPEnwmvj4eXLg57N\n",
"iMlQ2WlV2KGJhKZtp07UlgRj9lf0goWXBOVPKGrvse0eClVAPXkzKye45pB5W4kzSIvTGH/mzAaM\n",
"puuBv+aobj158Zrga+pza6vYsPgSjfFLa5PK2L3FrYTqTp/ym6p2HBcLkn/dhdjas3Z4nkYhykbu\n",
"VOIvcGYDRrPn7ldz8QvH8uuj5/Np92nujz4Qdlwi2dbwBdYha9nzgznsW3uFx/xjAOu733gO6nh7\n",
"PlyIzQU9bF2gYskCLn5hCTCaSQt2pWjaf5oxLOywRLKu7iaqE38M/7IocRPVS8Us7t2P5ErBZ6/4\n",
"fNvNTqfNp/asOSzcMjEfh2/Coh5/AbO4lQBTgMkETz07jlUDJ/C75wayYY9D3VkTboQi2WN29gJ4\n",
"4ChO/BEM+1UwC6e6J3DafPeZ5SGH12IKvsdvZhWJCxbSPMOByR7zao/5Go/5Y5S+fjYDZ/wVmGaG\n",
"TqxS8Cxu3W3kjx+H/Y6gqBrabaxfCiHHT8bLF2ZWnpgQk3lb6vG3Lha3tmwoGcNNa+LAXe7s5LGX\n",
"IuFK+YLth0fezP91eYnhK37RWsfuU6UbuKQhHelY/V26v34Onwycaca97jTyRDORkDV2s9Xq9q9a\n",
"3Ko85m96zKuBibaPjSmUqZf5Tj3+VsyM9u7okY2St8xGzaFoxkhGXA+LLoYj7wymXu41cgHnvny8\n",
"x1ydlu0U/Bi/5FiFdbS46eQqeaxzEbUl8OlBcPGAYN59bQks33+Tkn7uKPG3bv8LDAo7CJGG2BF3\n",
"dYHd96WoGrovgdveidwF27CoZEPr9s26XpPZgNHQ4zKVb5a8UR5fyOqL1nHI0HU8vbAztSXbHkay\n",
"cMuODyOJOJVskLQESf+oW4PKnXVUvllansWtrcd8C4BdV1xC0T+HB7N6+gUXbFm/gdJ3pjdUQlkC\n",
"KtkgKbGj+71OzTUD+Ov5270zulXULpH80+A0zRMWzWf46rEe86Fhx1fINJ1TUuO7dGDFUfXXFVVD\n",
"n9X7hBOQtHpbp2k+mJimuRZ2HTyUx3b5MRk/P0oypYu7UbC42woG35G4aMa28s3vq3yz5EhdXZ3z\n",
"joFD7k/U1Xm5mL8d+t2wQxP1+KOh5pP/4rl9ejHyh2WsPApKX4N5/k9qV/0q7NCktfqXUmpL4Mm7\n",
"YMIxQYmF2hKCcXwJm3r8EeD+xlNs+OulvPW3lzjlX+HlJS9QO+IVeGNe2LFJ62Jxa2unf3cG7boe\n",
"SFE1DJwRybo6+U7TOaOiYskC4BXgDC6ZfxXMn6wbZCQHvmLvv1TR86OrOXjIT7fW1dE0zYxpOqek\n",
"Jbl8s8e8OvH6DuCSRB0UkYxY3Hp4zFckvR6jaZq5oemckhKL2xjgxbokb3FrDzwP3OIxfyjU4KSg\n",
"NDhN88jFjzPm44uBQzzmX4UdY2unxC/NZnEzj4X0P18KVqOPP1y06UL/cNnMsOOLAhVpk2arS/oW\n",
"tzbW9f0uZtxkplle0oS6aZrfHg+jJ217/OE/+p8XdmiSOiV++SWX9jsZOBy4NuxgJN8lqmk+OwUG\n",
"37WtmqamaRYUJX65gTZbHgbOAy4z44iwA5L8ZKMv6U7xbgdTVA1H/EbTNAuYxvhlKzPOAGLA19zR\n",
"U42kHru+8/d5a8CP8M/3j/rjD8Oki7uSNRa3w4GTqfD+MGM3uG8XlXCW7WmaZvgKvkibbuDKK1XA\n",
"C3DS29D/Tpi9x7a3JvY2G4CSf/RY3EYAe3rMHwbwmM9KFFlT776F6QYuyRmzUXNg9sgd31EJ56gx\n",
"oy0VdhBQ7DF/Iex4JFDwPX7JQ4es3od20+DVC4DE75ZKOLdqO9yUZf+spUt5NV/bvavH/MSw45Ps\n",
"06weqW/5rp8w+HbovDJ4rRLOrV9d7fyiGSNh5rEcNWIkp906jr8/9r9hhya5ocQv9f3zs19y76HL\n",
"OObnULI8SPrPra2iRiWcW626m7JGTA7+n3d9Hx58tQ2ftf122KFJbijxSz1bSzgv+tt8rugFr769\n",
"iA1j3oA3Zocdm+RK5yJsCyz5LlzRC164Dr74F3RTVuulxC87qliygItfWAL8DxfO/4yOJ/YBzgg7\n",
"LMmVdbX0eQq+cb1uyooIzeqRepJLOANrgWLW9PxvfvNKObVdD3Xno3AjlGxruPCabsrKV7qBS7Ju\n",
"+xLOiXUlPP2fd7Hgqt2Ab7qjqp6tgMWtKzAC2KCbsgqHEr+0CIvbCL5qO4Gfbj4YuMWd34Udk2TO\n",
"4tYbOM1jfmPYsUjqNI9fWsqfaLPlDWAvoHvYwUhmzCgGfgZ+jTtK+hGkHr9IhNgJ1+5Hx88WMvv2\n",
"x9nc8Qcatis8GuqRFmVxKyJ4Vu/lHvMvwo5H0mNGMW03zmX4Teto89WJ/lyFkn4BKvihHhVpKywe\n",
"81qL25NAbdixSNPql2Io2Q32HkhR34WUx2/3+BYl/QKjIm2SF8wGjIYel6l8c37aOk3z6ReL+fa5\n",
"MPeXcPQ317Jws6ZpFrCC7/FL4bKR3a6j5OQfUf27btvWqnxzXlnVdwJzZwSlGBZdBENuDZ6PW3vW\n",
"+aiscqTpzl1pns7tL6Lr2d3qrSu6+QAOafvTkCKSHXQKno+74Co4+xRYcLWejyuAEr801+zBK+j/\n",
"+8St/Wyr4rm0p8b/88WpL/ZkwANBCQaVYpAkSvzSPDVfrmfeFDj+eujxYpD0502Bmk3rwg5NwIwO\n",
"/J0b6VOxlnlToLpn8P9n5OC1lL4zPez4JFwa45dmWnEbtVf15u8jDmDC0UGPsvbaD+CD28OOLOqs\n",
"7aZLof3RnPLRfUwtmxiM6XfqSO36DSzcMp1Jy2rCjlHCpVk90my2V+9xHFPzCyoP+Yhj3ivj+QMm\n",
"+6dPTws7riizET+J0fW9q3l0xgB3qsKOR7IvG7lTQz3SLBa3EiYtK2fAx4f6p3OHcmhVXy7546BE\n",
"dU8JgRnHUvmTi6k6/lQlfdkZJX5pruHA5KQqnpuAYwiqPUoLs45rDgUewdud6a9cMDfseCS/aahH\n",
"ssbitq/HfGXYcUSNxa0ja3tUce+z1/nqPveEHY/klmr1iERM/TIMiTumS9+9m0nLijzm94cdn+Se\n",
"xvglL1lFm6ttvxf/aEa/sGNpdaaW7crQdtMomjGSNg8eS9GMkQxtN42pZbp/QlKm6ZySfeavsaVD\n",
"O+AeM45xZ0vYIbUayWUYdv0MOqyHx1SGQdKjHr9kncd8Lv848ufARuCHYcfTunTeVoZhwEx46jaV\n",
"YZC0KfFLTrjzFcP+8zr2WfhvZvQPO57WY10tRauTyjD8l8owSNqU+CV3TrymmMF3PE0w5NM27HBa\n",
"hcH2MRf2RmUYJBOa1SM5ZUYb+OWN8Owg6LSL6vY3nxkD6Pd/L/DVA3fxXs2hwfDO+g2UvjOdSctq\n",
"POazwo5Rci8vp3OaWS9gMlDs7uN2sp0SfwTYyV1jzBtzDtxSxoY9gpUdxy9jxKz7/MnP4+FGVzis\n",
"wxelbNrtJWCyOw+EHY+EJy+nc7p7lbtfkO12pUA9N+gYTmxXxvfLwbYE49HHdSnjucOODju0QmFx\n",
"My4aVMl+z/9BSV+yIaXEb2bTzWyVmb2+3fqTzOxtM3vPzK7JTYhS0Gq6tWPur2DFcChekVS+eY/2\n",
"YYdWKDzmzh5Lv0nZs5eFHYu0Dqn2+O8BTkpeYWZtgTsS6w8CzjSzA7MbnhS+dbXUlsCL18IVvWDB\n",
"lYnph5qFkg6P+Xv+XIUekC5ZkVLid/fngTXbrR4MLHX35e6+CXgIGGtmXc3sLmCQvgUIrLiNjuOX\n",
"bZ1+OP4E2O+MjVCj8s1NsLidYXG7Luw4pPXJ5M7dfYAVSa9XAkPc/XPgolQaMLOKpJeV7l6ZQTyS\n",
"jyqWLGDV0je5d/gyaq5uz5Nd2zN4xX58/sgY4LGww8tXZhgVzAa6hh2LhMvMyoHybLaZSeLP+Gun\n",
"u1dk2obkveGUbhzvXzxbV74ZO+Ha/eixcJHZqSe782SYweWL+sXX9iiG0v7cOegSfrC4GlRbP8oS\n",
"HeLKutdmFsu0zUxm9XwI9Eh63YOg1y+ylcd8VlLN/mDdH3/+AftX/geH3nu/Gd3Dii2v1BVf6zxt\n",
"JKfXDKXTj0o4qubXTC3bNezQpPXJJPH/BehjZj3NrANwOvBEOg2YWUXia4xEzdDb59Fm81PAXWGH\n",
"khdW9Z3A3JeKOeYGeP0M+PoNMPelYlb1Oz/s0CQ/mFn5dsPjzW8rlRu4zGwGcCywB/AJ8BN3v8fM\n",
"RgG3AG2Bu939xpQPrBu4Is+MDsBAd14JO5awmZ1WCTOPpWR5MPvplqqgJAOnzXefWR5qcJJX8vLO\n",
"3ZQPrMQvgMWtHXCixzzSJRzsuIOr+PinPSmbF1TeHHZzcL9D7Vlz3J8aFXZ8kj/y8s5dkTR1AsZZ\n",
"3KJ9Q1ebVXcy4Op1Kr4mLSHUB7Ekxqs0jTPCEhd+zws7jtB9ffWbTC2eEDxQpVNHatdvYOGW6Uxa\n",
"VhN2aJIfsjmtU0M9kjcsbt2569WN/tFha8OOpSVY3PYDLgOu8lhIf4hScLKRO/XoRckfGzvPpssT\n",
"+5v9ZDF0bJf0IPHWWnJ4NTBfSV9amnr8kjds37LTOazTfTzzpw7UlgSVPEcOWcvCzRP94/dbzfNk\n",
"g2cUcAtwhzvvhh2PFJaCv7irefxSz4f9zuOZP3VgxGQoqQoqebbGuezfHv8M3d4cBnwQdihSOLI5\n",
"jz/UxO/uFbqwK9skHiT+0WFwRVkwrbGVPUjcjHOoGnEwvZ8+w53asOORwuHuldkqc6Mxfskj62op\n",
"qoa9/wK/eTlpLnthlnCuX3+ncxG07UjRGf05eOY1/j+zl4Ydn0SXEr/kj9J372bokKHMfamY2hJY\n",
"0yeYy75wS2HOZa+rvzN3RjHFH8C+C2Dv677gmY2rww5Nok2JX/LHpGU1TC2bWG8ue6+qLhy2uTTs\n",
"0JplVd8JzJ1RzIjJ8Nap0O8JeGzBbsHPR6u5WC2FRzdwSd7wmM8iKDi7NSla3PoA74cVU2YS1ywW\n",
"XLWt/k4ru2YhLUc3cIkUAOt83LN8vf8IDHjxGtXfkawo+OmcIqmyuPWziYOXWLuNhfNc565Vv6X7\n",
"72t49kbV35G8ojF+KRRLeeVfZ7JllyfMGOLO52EH1BgzdgNmcE7ZH5j993NVf0fyjYZ6pKCY8Uvg\n",
"UGCUO5vCjmd7tt8LnRjxkyr+cNezfN73LHe+CjsmaV1Uj18ix4y2nHjlElZ9Uc1ry6uD+fHramHF\n",
"be5vhFrT34wi4An6PLmZTZ1O8aryLWHGI61TwRdp06weSZc7W+yM+17g5W+Ph9kdtr7RcXwfO7nr\n",
"kf7k5/Ew4rJ+T3SEbz4KrOa9k8e7o6QvWZXNWT0q2SCF5w8DejKsXQeKEs9wL6qG47qU8dxhR4cW\n",
"06jLf8+AGUXAeHc2hxaHtFoq2SDRVtOtHfOmwAlXQcc1ULMHPHMT1F7YIk/x2rEUw7paBn35v3zn\n",
"7I/99TOV9CXvKfFLAVpXS20JPH99UMxt641RLVTTp64Uwx9/V8yXxdCuFvYfMpS7ek2kokUiEMmI\n",
"Er8UoBW30XF8H4Z1KeOWquDGqOfWVnHEkkrba3FP/3jQ8pwevq4Uw5mnwIph0GF9UD5apRikQCjx\n",
"S+GpWLKAVUvf5N7hy6i5uj2Vqzdx7ou1tNvnRN7qcKUZ33Nnbi4ObQfPbA9lvaktgSemwaUHqhSD\n",
"FBwlfilEwyndON6/eLa6boXFrQSWDeezg9YDD1rbL+/hqw6xbM6usbabhjHh5qf5aOAmNlTDkNvZ\n",
"+o2jgMtHS/SEOo8fiKPpnJJlZpQy9ryX+ehr63n5kuNhwOHQ47Lmzvk3oxvwc2AUZU/fwBeTPgnK\n",
"LSfKRxdV15WPblWPiJT8kjSdM6YbuEQaYJccuDvTn7+Svv9+PO+u35MN95dtfbPj+GWMmHVf8pz/\n",
"BmfqlL57Nxct25u3Tp3MzEcfBGLurN22bb+gFAPrN1D6zvRW/FB4ySO6c1ekCbbbiGcY3f14lp4I\n",
"i88PeucjJkPl28/4F8+esHW7vXqP27EXP2Qtr22cxMm7vO23v/PXMH8OkToFf+euSM7VdGvHwkvh\n",
"2J/B8hHBePz6vaD3ymPMeBFYw5BbO7Pbwf2Ye1/w0JTNu8DuVfD4S8XUnnWO364SytK6KPFLK7eu\n",
"lpXDYdbUbQ9DcQP+vBi4Ftid9jVHsbFDr3oPTZn6qmbqSKulxC+tXCNz/jd88FN3ng+2ue4Js1GH\n",
"UVTdg2E3a6aOtHoa45dWzeJWwqpd7ufe4UXU7NGeXVdv4vsv1lK6cbzHfNt00AbH+DVTR/KPxvhF\n",
"mtbInH+uQsQuAAAFZklEQVSGA9tm4DT0oHc9NEVaKfX4RUQKSMH3+FWPX0QkNdmsx68ev4hIAclG\n",
"7gz1QSwiItLylPhFRCJGiV9EJGKU+EVEIkaJX0QkYpT4RUQiRolfRCRilPhFRCJGiV9EJGKU+EVE\n",
"IkaJX0QkYpT4RUQiRtU5RUQKgKpziohElKpziohI2pT4RUQiRolfRCRilPhFRCJGiV9EJGKU+EVE\n",
"IkaJX0QkYpT4RUQiRolfRCRilPhFRCJGiV9EJGKU+EVEIkaJX0QkYpT4RUQiRolfRCRisv4gFjPb\n",
"DbgT2EjwkJUHs30MERFpvlz0+E8FZrr7hcA3c9B+i0o89SbvFUKchRAjKM5sU5z5J6XEb2bTzWyV\n",
"mb2+3fqTzOxtM3vPzK5JrN4HWJH495YsxhqW8rADSFF52AGkoDzsAFJUHnYAKSoPO4AUlYcdQIrK\n",
"ww6gpaTa478HOCl5hZm1Be5IrD8IONPMDgRWAj3SbL/Bs21TZ+B0z9DZOKOnG2dzjqk48yvOxvYv\n",
"hDhbw99Qc44bpTibI6XE7O7PA2u2Wz0YWOruy919E/AQMBZ4DPiOmd0JPJFGLOUprkvn/Uy3T7WN\n",
"nbXbnGM2Z59U2thZu805ZnP2SaWNnbXbnGM2Z59U9t9Zu805ZnP2aWr/ptpM95jpbp9qG021m+5x\n",
"090+1Taaajfd46a7fVak/LB1M+sJ/MHdByZefxcY6e4TE6+/Bwxx90tTbC+cp7yLiBS4TB+2nsms\n",
"nowSd6aBi4hI82Qyq+dDto3lk/j3yszCERGRXMsk8f8F6GNmPc2sA3A66Y3pi4hICFKdzjkDWAD0\n",
"NbMVZnaeu28GLgHmAm8CD7v7W7kLVUREsiHli7siItI6qFaPiEjE5FXiN7PdzGyRmY0JO5aGmFl/\n",
"M5tqZjPNbELY8TTGzMaa2X+b2UNmdkLY8TTGzHqZ2W/N7JGwY9mZxO/lvYnP9Kyw42lKAX2uef97\n",
"Wih/85Be/syroR4ziwPrgLfcfVbY8TTGzNoAD7n7aWHHsjNmVgL8wt0vCDuWnTGzR9x9XNhxNMbM\n",
"xgOfu/ssM3vI3c8IO6ZU5PvnWqcQfk8L4W8+nfyZ0x5/OjV+Emf8N4FPcxlTJjEm1p8CzCK4Uzlv\n",
"40z4N4KyGi2mmXG2uEKoP9VKP8uG5Pz3NJMYW/Jvvrlxpp0/3T1nC3AMcBjwetK6tsBSoCfQHlgM\n",
"HAj8DPgVwSyh35P4NpLrJZ0Yt9vv8ZaIr5mfpQE3Ace3ZIzN/TyBR/I5TuB7wJjENjPyNMbxib+f\n",
"vVv6c21unC35e5rpZ5nYPud/8xl8lmnlz6zX40/m7s8nSj0k21rjB8DMHgLGuvu/JV5/H/jUEz9x\n",
"rqUTo5l1Jyg7XQQ81xLx1UknTuAbwPFAFzM7wN1/k49xmtkq4AZgkJld4+435WOcwG3AHYmx0xa7\n",
"VyXNv5+fA/cn1nWlBT/XDOK8jBb6Pc0gxmNpwb/55sZJ8K0p5fyZ08TfiOSvzRDc7Tuk7oW739vi\n",
"Ee2owRjdfT4wP5yQGtRYnJcCt4cTUoMai/Nz4KJwQmpQY3HWAOeHE9IOdvr3A5Ann2sqcd5GcFIN\n",
"Syox5sPffJNx1kk1f4Yxqyd/riY3rhBiBMWZbYUQZyHECIURZyHECDmIM4zEXwg1fgohRlCc2VYI\n",
"cRZCjFAYcRZCjJCDOMNI/IVQ46cQYgTFmW2FEGchxAiFEWchxAi5iDPHV6hnAP8gePD6CuC8xPpR\n",
"wDsEV6qvy/WV8kKPUXFGM85CiLFQ4iyEGFsyzry6gUtERHIvr0o2iIhI7inxi4hEjBK/iEjEKPGL\n",
"iESMEr+ISMQo8YuIRIwSv4hIxCjxi4hEzP8DY8K8YajT1GUAAAAASUVORK5CYII=\n"
],
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc5b8713f90>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"loglog(freqs,appAna_r,'bo--',freqs,appSyn_r,'gx:')\n",
"gca().invert_xaxis()"
]
},
{
"cell_type": "code",
"execution_count": 177,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAXcAAAEFCAYAAAAYKqc0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
"AAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4VNW5x/Hvj2u4B7wErMhNIyr1rqAUDQgiYK2tgtZT\n",
"tFVpSy1orVgtpx2m52C12IvoKadSbamn1aptrVUEUYhWEFSq1BsgN8ULQYVgEIKI7/ljT0gCk2Qy\n",
"mZk9M3k/z8OTzJ59+RGSN4u1115LZoZzzrn80iLsAM4551LPi7tzzuUhL+7OOZeHvLg751we8uLu\n",
"nHN5yIu7c87loQaLu6QbJb0q6WVJf5LUVlI3SQskrZb0uKTCTIR1zjmXmHqLu6TewATgRDP7PNAS\n",
"uBi4AVhgZsXAk7HXzjnnskRDLfePgN1Ae0mtgPbAu8B5wJzYPnOA89OW0DnnXKPVW9zNbAvwc+At\n",
"gqJebmYLgCIzK4vtVgYUpTWlc865RmmoW6YfcA3QGzgE6CjpazX3sWD+Ap/DwDnnskirBt4/GVhi\n",
"Zh8CSPorcBqwSVJ3M9skqQewOd7BkrzoO+dcEsxMTTm+oT73lcAgSe0kCRgOvAb8A7gsts9lwEP1\n",
"Baz6A0Rrvq5rWyLvpfOYbMyZaKZsy1nX8bmQs6FzpuLf0HPmVs4M1qQmq7flbmYrJP0BeAH4DPgX\n",
"cCfQCbhf0hXABmBcgtcrTXBbIu+l8phEzlHfeZO5ZjLHJHKO+s6bzDWTOSaR4+s7bzLXTOaYho5v\n",
"6JyNvWZj90/0HA2dt7HXbez+iZ6jofM29rqN3T/Rc9R33mSumcwxTSZL45S/kixVv4XSSdI0M5sW\n",
"do765EJG8Jyp5jlTK4dyNrl2NtTn3lyUhh0gAaVhXFQaMBp6ToZOBVBRCRtnmr0yt55DSjOVrYlK\n",
"ww6QoNKwAySoNOwACSoNO0CmeMu9mVFUY5jVtz1lxVfsLdhFq+9i4rodFrFHa+2rAaPhtNtg9uHV\n",
"WyesgWevbqDAO+eaIBW104t7nlBUY4DFFrHyGtsKgcE1i7a69xvLoFazmb+sC5WFUFAOIwduY+mn\n",
"E2zT2gdqnfO4Y5ez+ukTqawxu0RBORSfudxWrDg5/X8r55qnVNROnzgsfywGpscKelVhnx7bTmxb\n",
"Bz7sdyXzl3Vh2FT4/D1w9rUwf1kXys6fJfHdWmdc02cnw6YGBR2Cj8OmwprelZn6SznnkuPFPU/E\n",
"WuxTgZsUVW+Cwr4L6Fdjt/vou+NgKgthyRS44FJ45asELfPyd7ix8whFdezevXu/1YbF1wcFvXBD\n",
"8HHhdNjR70iJQRn7yznnGs2Le36ZTDBNxHpgBvBHgqGqAFjEvsjqTmUUlMPpM+BX66H/Q7GW+Xvv\n",
"0rbiWmDN3rOd8fo7tLh+A0umwDV9YOX5UDllLfS9F7hZomXVrtKA0dKoedK40uDjgNGZ+Ss75+Lx\n",
"Pvccp6hkkeAfUVH1Aa4HbgGmAFNr9sFDXX3up25j6Z79+tz37j9kx60s67+Zry7+PPd+7uv21vr7\n",
"au3jN16dSym/odrMKar2wDPAEKA1ZW3vYc7p7dhxYCvaf/Aply2pZOn/TLF/XfF6jWNio2WOvBw6\n",
"toPtOyladXfc0TLV/fZTLWLlNV8D3YGPLWIbpVHz4LGR+yccPc9s7qi0fQGcy1Ne3PNcIiNgFFUv\n",
"i9ibOrdbhIVjLmXnPX33nqBgciVH6A37923H7nfyJl6fYCbQTyxi/6djB/6bN+Z/fr9RNUeMfNn+\n",
"vSypazvXnHlxz3N1tJzvB/5lEau1QEo9ref5ZnPPSWvODkOf4IKWZ/GXP8GOg6tH1ZSufMI+fnJE\n",
"Oq/tXD7yoZB5rsYImOk1RsB8G7h1/707FcQ/S8c6tqdQ23f+h4o3tzE0UmNUjX3EjrJfpv3azrm4\n",
"vLhnuViBf4bYCBiL2DqL2Af771lRx9jz7TvTGA8A27L6bzzU9hKee+UprukDL7z8DJXnV8Ar3dJ9\n",
"bedcfF7cs1ysK2Y8cBEwpeohpf1tnBmMUKnpyrXw1u1pjhiY9uoSrnrmVeAovvPP/py4/jLg5xIH\n",
"ZuT6zrlavM89i9U3WmXfIY5QNSTxsEl7R8Hw1u2ZGIoYJ+eRwGQe/s3ttvybK9N9fefyjd9QzXOK\n",
"6qvAYw3NFxO2ROe1cc4lxot7nlNUfwR+ZxF7IuwsyVBU3wKwiP0m7CzO5RKfzz3/XQbsqetNiSnA\n",
"KjMezlykRpkH+CRjzoXAW+45SmIEMAc42Yx3w86TCInTgVfM+CjsLM5lM2+55yFFJYJ1aqdbxDbE\n",
"3UccBtwDXJwLhV1RdQam03p7K3Z3OEDiIjPS16pwzvlQyGwTmwTsPuDteO9LtAUeAH5uljNLhlUA\n",
"z3Dg69cC/WDW//gMks6ll3fL5BiJXwC9gQtysfUrXXY5HHIn/LRl9VafQdK5mny0TJ5RVCOBJyxi\n",
"e2+i7r9A9aF/hdl/NmNbiFGTplG9VrPy90ewYWj1Rl+6z7lavM89jyiqdgSjY54BPoY650nvBwPe\n",
"hhxt5X5wwG5OvR02ncDe+eSHTYVSX7rPuVRqsM9d0pGSXqzxZ5ukyZK6SVogabWkx6W6Hot3ibCI\n",
"7bSIXWIR+7h6a8/JtQs7BK8Pm5TZdCn0QtFGHr47ztJ9uyvCjuZcPmmwuJvZKjM7wcxOAE4CdgB/\n",
"A24AFphZMfBk7LVLgqJqHf+dOmd6bJe+NOm2cSaVU9awbHKwdN+SKVB5XebmwHGumWjsaJnhwBoz\n",
"2wicRzDOmtjH81MZLJ8pqjFVE4Apqs8ByxRV19hj/DWEN9Njupi9MpeihT/ki6dt5vUemxh+2lsU\n",
"LbrRb6Y6l1qNLe4XA/fGPi8ys7LY52UEK/O4xCwmmKO90CL2DnAB8N+x7TVsnAnffa/2tgzO9JgG\n",
"iqqQietK6L31SLvv3R4M2HQcE9eV1D3bpXMuGQmPlpHUBngHONrM3pe01cy61nh/i5l12+cYHy1T\n",
"hxozKc6gjsWsAaRH/g0PGuzYmsmZHtMl7iRj3/9cVzq/e7pPMuZcINOjZUYBy83s/djrMkndzWyT\n",
"pB7A5jpCTqvxstTMSpNKmn+2ETyotB7oE7+wczKcWwjn9jNjd8YTpsF+i3BHWk3gvmev491TziYS\n",
"VirnwiWpBChJ5TkbU9y/SnWXDMDDBEP3bol9fCjeQWY2Ldlwea4IGAv0IViEI17L/fvAr/KlsMfV\n",
"Ys8mWu5aDFwFXB92HOfCEGv0lla9ltTkpk5C3TKSOgBvAn3MrCK2rRvBYs2HARuAcWa1i5N3y8SX\n",
"yCIcEr2B5UCffJ9oS6Iv8BzQy4yPG9rfuXyXsQWyzexjMzuwqrDHtm0xs+FmVmxmZ+9b2F29BlOj\n",
"kNdYCHtwjX36AT/L98IOYMY6tOefwKVhZ3EuX/j0AyFQVF8GDrSIzQ47SzZQVJ35pMO/ueWDnewp\n",
"OMaMz8LO5FyYMtZydyn3L2Bp2CGyhUXsI/TZIPYULAFfUNu5VPCWu3POZRlvueeguqcacIrqAEXV\n",
"P+wczuUDL+4ZFFtlabmiOizu+0JSs/43GQqcG3YI5/KBd8tkmKLqbBGLOwJGYiTwTTMuyHAs51wW\n",
"8cU68ozEAuAeM/4QdpawSbTwUTOuufI+9xyiqAoVVXGd74vjgaMI1k9t1jSl+610Xfua5IvJOJcs\n",
"L+6ZM4DgEfu6fB+43YxPMpQne3Use4Y9rbfh00g7lzTvlskCEj2BFUBfM/xJX0BiLDDZjCFhZ3Eu\n",
"07xbJn/0AP7LC3stf6PF7l4SJ4UdxLlc5C33DFBU4wnmMF8XdpZcoah6s+WgF5g5bhdsfiNYlWrj\n",
"zFyey965RGV6PneXvALg05obpAGjgwWwOxV44Ypj2jFH0+6k7XBHL+CQYOOEftIA/OvkXMO85R6C\n",
"oLCfdhvMPrx664Q18OzVXrgCOu7Y5ax++kQqa6y+V1AOxWcutxUrTg4vmXPp533uOavn5NqFHYLX\n",
"h00KJ08WWtNnJ8OmQtFL0G1NUNiHTYU1vetYNNw5V5MX9zRSVL0U1Z/3f6dTQfwjOrZLb6IcsuOT\n",
"7SycDud8Dw5/LCjsC6fDjt0VDR/snPM+9/QqA362/+aKOlqf23emNU1O2TiTyin9+PvvDueaPvCr\n",
"9VB53Vp46/awkzmXC7zPPQRBn/tZc+C2GnOXX7kWlk72Pvdq6t5vLEN23Erpce9RsqIH/2x/nW1a\n",
"+0DYuZxLN59bJospqg7ADovE/wJLj62EP1bAJx8HLfa3bvfCXq3murLAATw19R4G3rGSgm3XxllI\n",
"3Lm84kMhs9tNwIvA7/d9Q+JwGNUVRg0wqz1E0u21d51ZRVXInratufO55Uw+cjDwaNjhnMt23nJP\n",
"IUU1huBhpfLY3O2tgA7AYIvY3oIk8SPgYDN8dEyCJL4GXGDGl8PO4ly6+VDI7LMYmK6oCmPdMR0I\n",
"uhYW77Nfd+CPmQ6X4x4DzpLwEUXOJcBb7ikW6yueA3yPYKbHqd5H3HSK6iBevXAFDzwwwcy7ZVx+\n",
"y1jLXVKhpAclvS7pNUkDJXWTtEDSakmPSyps+EzNwqdAG2AtMMMLe8p8yPtHzYfPxoQdxLlckGi3\n",
"zG3AXDM7CjgWWAncACwws2LgydhrF/SzrwP6AFNiLXnXRBaxzyj9ySRocW3YWZzLBQ12y0jqArxo\n",
"Zn332b4SONPMyiR1B0rNrP8++zSrbpmaw/f2jvKo8TrcdM65XJGpbpk+wPuSfifpX5JmS+oAFJlZ\n",
"WWyfMqCoKUHyxARgSVUhj32cSjCszzWRomqhqJ5TVN3CzuJctkukuLcCTgR+bWYnAh+zTxeMBc3/\n",
"9N2ZzR0LgU01N1jEyi1ij0q0kfi1rwuaPIvYZ8DFwNawsziX7RIpNG8Db5vZ87HXDwI3ApskdTez\n",
"TZJ6AJvjHSxpWo2XpWZW2oS8Wc0itryet88GjvWHlprGFzxx+UhSCVCS0nMmMhRS0tPAlWa2Olas\n",
"28fe+tDMbpF0A1BoZjfsc1yz6nOvj8SfgGfM+HXYWXKdjnykiHdP7GoVh6wMO4tz6ZCxuWUkHQf8\n",
"luohft8AWgL3A4cBG4BxZrVvGjan4q6oHgKujdeylOgAvAMcYcb7GQ+XZ3TtYS/yyK8rbfW5p4Wd\n",
"xbl0yNjcMma2AjglzlvDm3LxPPNj4M063jsPeNYLe4r8/a6zWDfiTYlOZvj87s7F4U+oZoDEQ8Bf\n",
"zLgn7Cz5QmI+cKcZfwk7i3Op5nPLZAFF1So2vW99vk5wI9qlyiHPP0XXtWPDjuFctvLi3nQn0UDh\n",
"NqPcDF9lKZUuuqArhy4dLdEy7CjOZSPvlkkBRdXKIuZDHDNM4ifAL8183LvLL74Sk3PO5SHvcw+Z\n",
"ojpZUXUPO0dzpaiKFZVP7eBcHF7cm+YsoH9db0oM9z7htPoccGTYIZzLRt4tkyYSRwFPAIeZsSfs\n",
"PM653OHdMtntq8B9Xtidc2Hw4p4ERdVaUf1MUcV9wldCwCXAnzKbrPnRj9r21VlTV0h0DTuLc9nE\n",
"i3ty2gJvxBv+KA0YDRctgamHwKjpwWuXNq0+qaCy6x5gVNhRnMsm3ueeQkEhP+02mH149dYJa+DZ\n",
"q81emRtesvwmcSUw3IyLw87iXCr4OPcsI42aB4+N3P+d0fPM5nrLMk0kugOvA0VmfBJ2Hueaym+o\n",
"hkBRfUtRXR3/3U4F8bd3bJe+RI5p2sH44YY+HRJ2FOeyhS/51nj3AZ3iv1VRGX/7dp9XJr0qeKnP\n",
"8/Dd30pb3gz+HTbO9K4w15x5cW8ki9g2YFv8dzfOhAn9ave5X7kW3ro9I+Gaq+e7/pg1u4qxP/QG\n",
"egPQbvwROrfbKfbIlmiY0ZwLixf3BiiqMcBii1i5oupmEduiqAqBwRaxR2vua/bKXGkAMHpS0BWz\n",
"fSe8dbu3INNs0fFDGNqpNwu3QmVXKCiHoZ37suiEL4QdzbmweHFv2GJguqKaCjyhqC4Dvg1MrbmT\n",
"RDugVayQezHPpB0HtmLhdJh0BNz7dzj2j7BwOlR+s3XY0ZwLi99QbYBFrJygkE8HLiRW2GPba7oM\n",
"uDPD8RwAFZVUFsL/PQZXDoYlU6CyEL/X4Zozb7knINYlMwNYD/TZt7DHnki9CrgmjHxu40zajT+C\n",
"Ezr35Vfr4fQZsGjbenb6vQ7XfPk49wYoKgF3AB2BCDCFfVruEkOA2cBRZqTvC+riUlSFlLW9hzmD\n",
"C9hzQCdaf7Sb8aXlFO0aH+d/WM5lPR/nnhldgD7AdRaxDcS6aGI3VatcBfzaC3toBlO0a7x9/OQI\n",
"vv18Tw7+/qsU7RoP+FzvrtnylnsDao6WqbFt72gZiR7Aa0Bvs7qGSLpMUZvtw9jd8WYzTg07i3PJ\n",
"ytj0A5I2AB8Be4DdZnaqpG7An4FewAZgnNm+fdG5W9wVVUdgtEXs/nr3E/2AYWbMzkwyVx+JtsD7\n",
"QC9fW9Xlqkx2yxhQYmYnmFlVi+gGYIGZFQNPxl7nk4OoZ5WlKmas9cKePczYRdFLL9Lx3XPCzuJc\n",
"mBJtua8HTjazD2tsWwmcaWZlkroDpWbWf5/jcrbl7nKXvnXScp6cvs7WnDM27CzOJSOT3TLrCB65\n",
"3wP8xsxmS9pqZl1j7wvYUvU6lQEzTVEdCmy1iH0cdhaXHIkjga+Y8dOwsziXjFTUzkTHuQ82s/ck\n",
"HQQsiLXa9zIzkxT3t4SkaTVelppZaVJJM+cy4F3gd2EHcckxYxV4YXe5Q1IJUJLSczZ2tIykCLAd\n",
"mEDQD79JUg9gUb50yygqWaT+L4xEezN2ZCqTaxxFdQbwukXs/bCzONdYGbmhKqm9pE6xzzsAZwMv\n",
"Aw8TtHKJfXyoKUHCpqjaVH2eQGHvCKzzdTuz2lDg0LBDOBeWBlvukvoAf4u9bAX80cx+GhsKeT9w\n",
"GDk4FHKf2R5bAy8AXwaO2ne2x/2OFd8CRprxlQxEdc41M77MXhPEHkSaTmwqAUXVC7ie+JOCVR8X\n",
"zCOzArjWjCcyk9Y515z49ANNECvg/wncoqh6k0BhjxkMtCEY2++ymL5y6U909F+uDTuHc2FotsU9\n",
"ZhKwmWC2xxkJTjJ1FTDL55HJAS0/6UDFIReHHcO5MORVcVdUY/aZ0AtFVRjrX0dRHaGovlvj7f8F\n",
"uhFMDDZl32P3O3/QJfMeMCe1yV1aPHjfNN4+7WiJDmFHcS7T8qq4U71qUiHs7Vf/RWw7BA9iba3x\n",
"XoSgK2YD8Wd7BEAaMFoaNQ/GLYJRR8OA09P/V3FNZUYF8C/gjLCzOJdpeXdDtcaN0hnAD4CzgJMs\n",
"YhX77FfvbI97t2nAaDjtttqLXk9YA89e7WujZj+NuvohNp6+01656KthZ3EuUT5aJt41g2GNxwPP\n",
"EXS3vNnQuPV6z6dR8+Cxkfu/M3qe2dxRyZ7XZYbO++Z3WDp5sm0e0OAkcM5lCx8tE9/5BFMH9CFY\n",
"NalL007XqSD+9o7tmnZelxH/uPNO3h/gN1Vds5NXxT3WrVICDGmoHz2h8wlBlx7x3/XFl3OBGZ+a\n",
"8VLYOZzLtLQXd2nUvKDfOt57wY1KaVxpffs1tK+ikqI6gWAM+lSL2FbYO5Z9KkkstybRAvgFXNgG\n",
"vr2h9rtXrgVffDlXKKrJiuqCsHM4l0mJzgrZBI+NhAn9pAHUvAFZx43K/fZLcN/uwH8D51nE9tQ8\n",
"Nlbg651OYF8SbYC7gV4w8gT4/ukwelLQFbN9J7x1u99MzSmPgy+B6JqXtN9QZe+zPhf+Ex78JrDT\n",
"jDfrvlF58bNw3/eBlgS/fD6CUTclelMz+EXQc3LQV15RCRtnNqYQx8ZEPwjsBi4yw7tfnHMZlcn5\n",
"3FPgyJMJZo5cDvxH3Tcq+wwAfgl8SrA4yIo69z2y1clqWzGBTzq9ALwKA4Yn+r+BevQnmAhtkhmf\n",
"JniMy3KxB9A6mLE97CzOZUIGi/uLT5lRo5VdURl/vxWLa+8HUsW8/XbTZ3DMa614/72RbOl0NdAX\n",
"ztkNt3auvePsw4MuFfYr7vW08pc39m/nstylw3/HO6ecAj89JuwozmVChkbLxLsBuXFm8DBQQ/sB\n",
"Y5Y9S7vx62ptK7h0Pbu2/Mo+LL7QjAHAAfDhmv2OBaBjO4n+EsdKwS+06n78x0bC/WcGH0+7rb6b\n",
"ui6HvT3wJp768aESrcOO4lwmZKDlPnpevBuQZq/MlQaQ0I3KU7bexmEPnMycd9fR9+MiKuwjRq/Y\n",
"StGu26rPx05pUx2r7mzfSTBiZgpwqMSLMORQmNW79n51t/JdbrMnp6+WWAMMAv4Zdh7n0i1nnlCt\n",
"Ma3ABuBk4Fv7zuIYf1TNlWth6eSqXxoShcBJ8N3ZcEef/a807imz+0tSkdllF+mzmyko/8x2dvth\n",
"2Fmcq0+O3VBtmtiCGjMIpuftE2963kT+N2BGOfCktHY1wVOs+/CHk/LW5UMOYOPpI2GGF3eXdRTV\n",
"GGb1bU9Z8RUpOV8Otdz7AdcSTAg2hcQW1qgnW8OtfJdfdOJvO/LilY8BQ30klMs26t5vLINazWb+\n",
"si5UdqVZTBwW65J5FvixReyBfZfISz7fgNFwmD+c5JwLnTRqHgX3jmTYVJj762ZT3McQzMn+kUXs\n",
"s9i2/abnda4hiqoT0Nki9k7YWVz+q93VEhtyXbT6Liau21FravFRk4+GxS8wb3k7CjdAeZ8mF/ec\n",
"6HOPV8CTmVbAOWAccCBwS9hBXDMwq2/7oKvl3i5UFkJBOQwaOIg/9LpK6C/AhRYxo/2WDaxvu5mC\n",
"8l6cPiMl4/WyvuWuqDoDIy1iD6QolmumEm1FOZcq1V0tP4Sil+GDI2HBrVB5yTymPXYbsKBqPqxU\n",
"97nnQsv9YKA47BAuD9TViprVdwKRsMO5XJJIQ0EXXfA4B/Y6kw8KYcn1cE0f+PtdUFkIdGxnEav9\n",
"5P3EdTuY1XcClZdcDpzT5JBm1uAfgkm8XgT+EXvdDVgArCaYca+wjuMskfP7H/+TiT9wzjwKthrn\n",
"X2oc8Ygx+jtGwVaDUY+Fnc3/5NYfivqO5UvF5cH3jwXfR187ZDvDD/jB3n3GfHsKXUYto2Br8L1W\n",
"uD7h77lU1M5Epx+4GniN6ikebwAWmFkx8GTstXNZrlMBlYWwbhj8x7mwZMreVlTYyVyOKSu+gqfn\n",
"dmHMRCjcAMOmwkszOvBm7/OrdrFHZs2gYNWtjBy4jYXTobw3LJwOI0/dRtGqu9MdscHiLulQYDTw\n",
"W6CqD+g8YE7s8zkES9ulnKK6RVEdl45zu+aoopKCcjj0OfjVejh9RtA14w+uuQRJtJAYAcecyEGv\n",
"Q/lhQXfLkinwyiXwRt9dtQ6YuG4HSz+dQOUl82DcU1ReMo+leyYwcd2OdGdNpM/9lwQPDdWcbbHI\n",
"zMpin5cBRakOFvMPgukGnGu6otV3MWjgoOCGVWF1K2rpnrS3olz2a6gfXedOPIZvvPICv3vqddjy\n",
"AW994SAOf6y6obBwOlTWbihYxB6N3c/J+ICQeou7pHOBzWb2oqSSePuYmQWLctR5jmk1XpaaWWmi\n",
"4SxizyS6r3MNqnnDqkubgzm88liWtvw2E99IeyvK5YB4N9xHHTWE3/cKbrif8r+v8cyUcWYt/qHu\n",
"T41l0MDZqWooxOprSSr/OvUOhZR0EzCeYOGMAoLW+1+BU4ASM9skqQewyMz6xzneLInhPIqqJdDS\n",
"IvZJY491LhGKqiulkX9TOm28GaVh53Hhq/WE6JIpQWv80zawfHmpffD00Fr77m3lH3n53ifci1bd\n",
"naphtSkZRl5fcd/nYmcC15nZFyX9DPjQzG6RdAPBaJn9bqo2obifAVxnETuvscc6lyiJgcCbZmwK\n",
"O4sLn3Tpcww86RQKymHotKC7pbw3YcwUm4ri3tjFOqp+E9wMjJC0GhgWe50yFrGngYtSeU7n9mXG\n",
"Mi/sTsc8MELir3DIcaw9Gzq/kxc33LP+CVXn0iXW/fd/wNctYrsa2t/lpvpulPL8xE30XLKE3yy/\n",
"noNO/YDTtv/P3n70gvKqfvQJtmltRm+I5uV87opqCLDCIvZR2FlcfrOI7VFUcxre0+W0fW+UHvUg\n",
"HHXDaczqe6Vt+vVyDZ3W3j5ruUfRf42pfkK0Yzsqt+9k6Z67MzFsMR2yruWuqGYBP7OIrU9TLOdc\n",
"M7LfjdILL4J5v4S3/3ue2dxRYeeLJy9b7haxiWFncM2LRAugjRmVYWdxaXBQ6x4c8nBQ2K/pU+NG\n",
"aX4/mdzYG6rO5Z/LB7/CybNua3hHl0skOkvcyq6e/SnYEtwgzYMbpYnKmuKuqE5QVJPCzuGaocXX\n",
"R1k+4fiwY7jU0Zcv+186b1wNFNJm5VV0nxXK/C5hypriDmwBXg87hGuGVn3pr1irwyU+F3YUlxhF\n",
"NUbd+42VRs2TxpVKo+ape7+xsVXbALXmmPu/Y8aVfHfhe2HN7xKmrLuh6lwY1GrXPXTY/KJt6/mL\n",
"sLO4htVe2KIQuq6B87+wg7kdvp7pYYvpEMZDTCmjqMbE1kFFUbWIfSys/s3rXAZdOqI1J935rbBj\n",
"uASVFV/B/GVdGDY1mHJ3yE/hzcvbU1Z8edjRskWY3TKLgemxAv+sojoWmB7b7lxmPTLrChb913qJ\n",
"lmFHcYnoVEDnt+GVi4MRME//CBbeBHTK6xEwjRFacY8tcD2VoKB/B/gWMDW23bmMss3HfGzGOWbs\n",
"CTuLq5+EoOPB9J0PX7i5WY2AaYxQb6jGCvkM4AVghhd2F6ZYt6AvDpPFdPJvjgbm0W5kKw7+7Xb+\n",
"+sdmNQKmMUIt7rEumWlAH2BKVR+8cyE5GviPsEO4+DR0mjj1jmX0XPwq137j+yz75PLmNgKmMUIb\n",
"LRMr5LcBR1vETom9no53zTjn6qAv3NzRnrlhe9g50i2nR8sAgwkW3j4VavXBDw4xk3MuRPuNXz9o\n",
"SKm+3+r5qlF0zaGwp4qPc3euBl3TpzOvjX2QBT871wxfCSzD9hu/XrAVzj9+O4taXZ4P49cTldMt\n",
"d0V1sKI6JazrOxdX1w2VtKk4hpaVQxve2aVcWfEVvDCnC18eH4xfH/af8NCKjsFydq4xwpwVsi9w\n",
"NvB8iBmcq8Ui9onEbcD5wPyw8zQ/PQ/FWsK64dUzOFYWku8zOKZDmOPcl1rEfhLW9Z2rx9+A82NT\n",
"AbsM0Og+QQYEAAAO10lEQVRJbdRm+03QtR9bjoADVvv49Sbyb17n9mHGG5x3eRv6PHl22Fmajc5v\n",
"P8ppv7iQA1+YxMiBzW4Gx3QIpVtGUZ0ByCL2VBjXd65BHxy1gA/6jwLmhR2lWfiw+Kvs6vIR3104\n",
"Ip+WugtTKKNlFNVZBMX9ibRd3LkmkDgYkBllYWfJV4rqG8DTFrG1YWfJNjm7zJ5F7MkwrutcoszY\n",
"HHaGfKGoxjCrb3vKiq+ATgVQUUnR6ruYyFHAM2Hny1c+zt25OBTVGN5rewUvHN6H5Udvqy5I63ZY\n",
"xB4NO18uqTV23QQyGDlwG0s/ndCcxq43RtrHuUsqkLRM0kuSXpP009j2bpIWSFot6XEp8TlhFNV1\n",
"impgU0I7l3az+rZnyyHDWfX48XD/mRTcO5JBrWYzq2/7sKPlnL1zr/8QLhsKo6+C+cu6+Nj19Kq3\n",
"uJtZJTDUzI4HjgWGSvoCcAOwwMyKgSdjrxO1FNiYZF7nMqOs+Ar+8a9OnDE99jDNVC9ISetUQGUh\n",
"LLkeDnkxGAHjY9fTrsGhkGZWdZe6DdAS2AqcB8yJbZ9D8MBHQixiz1jE3m1kTucyrKogXRc8TLNk\n",
"ihekJGjc2Bac9tbRdHgvGLPuY9czpsHiLqmFpJeAMmCRmb0KFJlZ1SiCMqAojRmdC0FFJQXlMPxG\n",
"ePMLXpCS9f4xxsEtNjHyjO0+dj2zGhwtY2afAcdL6gLMlzR0n/dNUp13ZSVN2/viYr5If75iEXsz\n",
"+cjOZUDR6rsYNHAQjyzrwq7O0PajoCAt3eMFKQGKqoVF7DNbNM0Ujd7IrL7tfex63SSVACUpPWdj\n",
"RstI+hGwE7gSKDGzTZJ6ELTo+8fZv9YdX0V1BLDWIvZZ06M7lz7Vw/eOvBy6dII+p3LAy9cw6ZE3\n",
"fbRM/RRVa4K1kEdbxD4IO08uSsVomXqLu6QDgU/NrFxSO4KJlKLASOBDM7tF0g1AoZntd1PVh0K6\n",
"fKHP/+kPtNt6gD131Ziws2QziVZmfKqoevn/0JOXiSl/ewALY33uy4B/mNmTwM3ACEmrgWGx1/WH\n",
"jcpvRLncdfKd/0d5rzMkPhd2lGyx38Ia/Uaupd2fXlSk5Rgv7OHLyENMiqoT8CrQ1yL2adou6Fwa\n",
"SZwKvGjG7rCzZINaDyd92hbGXgiVb2xnyZ5mtbBGOuTMYh0WsQrgcC/sLpeZ8RzTtEdR+WyqUOPh\n",
"pKnQsQy29YbHnvOFNbJExr5JLWK+ZJnLB3cCXww7RFZoV9CJUZNg2SR/FiALpb24K6p2iqo43ddx\n",
"LkN+ADwcdoissHNnBSu/BANn+sNJWSgTLfdi4McZuI5zaWcR+9AiabxRlQMU1YEAFL1xF8VTt7Hw\n",
"Jn84KQv5rJDONZKiasHiKTNZNulHtq3n1rDzZJKi6gYsAk4BRlQ/C9CxHWzfSdGqu33mzKbL2fnc\n",
"nctxxoErR9D2o+00btK8nBVbT7bIzN5TVCdbxHYDjxIBwEfGZKG0t9yZxukWsWfTdhHnQiBxJPBP\n",
"oJ8ZFWHnSSeNuepAdnabz6LoSrMW/xF2nuYgV4ZC7jctgXO5zoxVBNNdTww7Syrt92CSvrSQxSct\n",
"4dClXTnmwW+Gnc8lLlMPMRUCg70fzuUTnfK/wzhi7lwevrObbe+eF5Ng1XowqbILFGyD4UN28fyO\n",
"8f5gUubkRMs9VtinE0wk5Fz+OHfiIp7/9rN83H1s2FFSpurBpHOuhitPhbN+CE/8s60/mJR7MnFD\n",
"dTow1SJWnoFrOZcxFjHTNL4E+dTnHlukpDQaPJj04AP+YFKOykSf+wwv7C5vTdMQeve+WAcNKQ36\n",
"qEfNU/d+YxVVzs0eqagO5POrOlJQ7qsm5YFMFPcpsa4Z5/LPrL7tOfOT33LKcWfmwULaBzJw9RpG\n",
"DtzmqyblvkwU96nAdC/wLi+VFV/Bn19tjyxnF9JWm+1fkzjZIraSQyvvYemnE6i8ZB6Me4rKS+ax\n",
"dM8EXzUp9/hoGeeaQBpXCvefSeGGoI965irYUgyMe8rs/pJw09VPojUl0x6h/QenMveOM8x4OexM\n",
"LpAzT6jG+ty9sLs8FFtI+/QZ8Jvn4YrBcMfrsDO7+qirlw0svgI6FcAnRruLiviox8fsaX2sGRvD\n",
"zuhSy6cfcK4pqhbSnr+sC5WFMHsZDD3bWHzm8rCj1TKrb/tg/Pq9Xej+EuxpBcdN3MXzlZfapje8\n",
"sOchnzjMuSaovZB2bPKswp7PUHTOFKzVRFv1xXvDzgggjZpHwb0jGTYVNh8N/Z6Av/8OKi+ZZzZ3\n",
"VNj5XG050y3jXL6yiMWdPEvDbyyj6OXfS19834wnQglXU58dB7G+S7CgxjV9gmGOPn49r3lxdy4N\n",
"7Imf/lZtK94A2oSdRVG14Myuh/LRKhh4e/X49YXToTK77g241PFuGefSTFEdDBxkEXs1Y9c853sF\n",
"rD73h6w/630zbq89Z0xh8GDSyFO3sXTPBJ8zJvt4t4xzuWEgcCSQ0uK+/wiYikqKVt/F4NN6seHC\n",
"iZw0uy3rzxoBwMR1O5jVdwKVlwT3Biq372Tpnrt9/Hr+8pa7cxlSZzFOcuWi/Vrj3VbD6Rfu5ono\n",
"R1R++SpKpt1vi6Y16yUBc1VGZoWU1FPSIkmvSnpF0uTY9m6SFkhaLelxyZ9Ada5es/q2Z8Tue+g5\n",
"dWRKpiqomsFx2NTg6divnQOv/bw1lXe/aMafvbA3b4lMP7Ab+J6ZHQMMAq6SdBTB8mILzKyYYNGC\n",
"ZrHcmHNJKyu+gpf/qy2DflU9VcGzD3Rhc/8ra+62/4IZdUxGNmDrIRz11+oRMH94AtaNANq1zuDf\n",
"ymWpBou7mW0ys5din28HXgc+B5wHzIntNgc4P10hncsPnQpYcRksuDUoxkuugy9/HQqOHi5xEICi\n",
"GsidfTowqNVsCu6tbuEPbH0Xf/7GZRr24xel2M/dBy3K2dKvxgyOP/cZHN1ejbqhKqk3cAKwjGCx\n",
"3LLYW2VAUUqTOZd3akxVUDUccc5CqBy/0Iz3FVV74Ke8V/wJ8+/rwlk3wPYe0GEzPL60E22vGcJh\n",
"FfcBNwNgZbdz/DeP3dvnXjWD49I9PoOjS7y4S+oI/AW42swqpOq+/mB2MMXt35M0rcbLUjMrTS6q\n",
"czlu36kKqovxnQAWsR3AME0bV0plIbz0DZgwqPqBo8odq2ze77+393w+AiZvSCoBSlJ5zoSKu6TW\n",
"BIX9HjN7KLa5TFJ3M9skqQewOd6xZjYtJUmdy3UJF+NYC/+4P9T7wFFdT8e63BNr9JZWvZYUaeo5\n",
"GxwKqaCJPgf40My+V2P7z2LbbpF0A1BoZjfsc6wPhXSukfyBI5eK2plIcf8C8DTwb6Bq5xuB54D7\n",
"gcOADcA4s9rL6Xlxd67x4k5GVrTq7mTHw7vck5Hi3qSTe3F3zrlGy8hDTM4553KPF3fnnMtDXtyd\n",
"cy4PeXF3zrk85MXdOefykBd355zLQ17cnXMuD3lxd865POTF3Tnn8pAXd+ecy0Ne3J1zLg95cXfO\n",
"uTzkxd055/KQF3fnnMtDXtydcy4PeXF3zrk85MXdOefykBd355zLQ17cnXMuD3lxd865POTF3Tnn\n",
"8pAXd+ecy0MNFndJd0sqk/RyjW3dJC2QtFrS45IK0xvTOedcYyTScv8dcM4+224AFphZMfBk7HXO\n",
"klQSdoaG5EJG8Jyp5jlTK1dypkKDxd3M/gls3WfzecCc2OdzgPNTnCvTSsIOkICSsAMkqCTsAAkq\n",
"CTtAgkrCDpCgkrADJKgk7ACZkmyfe5GZlcU+LwOKEjko3m/N+n6TJvNbNhW/mT1n045J5PhcyNnQ\n",
"ORt7zTD+zZO5rufMru/NZDX5hqqZGWAJ7l6S4LZE3kvlMYmco77zJnPNZI5J5Bz1nTeZayZzTCLH\n",
"13feZK6ZzDENHd/QORt7zcbun+g5GjpvY6/b2P0TPUdD523sdRu7f6LnqO+8yVwzmWOaTEFtbmAn\n",
"qTfwDzP7fOz1SqDEzDZJ6gEsMrP+cY5LtOg755yrwczUlONbJXncw8BlwC2xjw/F26mp4ZxzziWn\n",
"wZa7pHuBM4EDCfrXfwz8HbgfOAzYAIwzs/K0JnXOOZewhLplnHPO5RZ/QtU55/KQF3fnnMtDGS3u\n",
"kjpIel7SmExetzEk9Zc0S9L9kq4IO09dJH1J0p2S7pM0Iuw8dZHUR9JvJT0Qdpb6xL4358S+ppeE\n",
"nac+OfQ1zZXv0Zz4mYfG1dCM9rlLigIVwOtm9mjGLpwESS2A+8xsXNhZ6hOb1+dWM7sy7Cz1kfSA\n",
"mY0NO0ddJI0HtpjZo5LuM7OLw87UkGz/mlbJoe/RrP+Zb0wNbVLLPd6kYrHt50haKekNST+IbRsB\n",
"vAa835RrpjtnbPsXgUeB+7I5Z8x/AndkLmXSOTOukTk/B2yMfb4nC/OFJgU5M/I92pScmfyZTzZn\n",
"o2uomSX9BxgCnAC8XGNbS2AN0BtoDbwEHAX8N/BLYD7BuHg15drpyrnPcX/PVMYkvp4ieM7grExm\n",
"TPbrCTyQzTmBrwFjYvvcm4X5xsd+fg7J9Nc02ZyZ/h5t6tcztn/af+ab8PVsVA1N9iEmIJhULPb0\n",
"ak2nAmvMbAOApPuAL5nZf8ZeXwa8b7G/USY0Jqekg4GvAAXAokxlhMblBIYDZwGdJR1uZr/JxpyS\n",
"yoCbgOMl/cDMbsnGnMBM4I5YX+bD2ZbPzG4G7olt60YGv6ZNyDmZDH6PNiHnmWTwZz7ZnAT/A0q4\n",
"hjapuNeh5n9vAd4GBla9MLM5+x0Rjrg5zewp4KlwIsVVV85JwO3hRIqrrpxbgG+HEymuunLuAC4P\n",
"J1It9f78AGTJ1zSRnDMJfmmGKZGc2fAz32DOKonW0HSMlsmVp6I8Z2p5ztTI9nxVPGdqpTxnOor7\n",
"O0DPGq97EvwWyjaeM7U8Z2pke74qnjO1Up4zHcX9BeAISb0ltQEuIkN9mI3kOVPLc6ZGtuer4jlT\n",
"K/U5m3jX917gXWAXQX/RN2LbRwGrCO7+3piJO+We03PmWs5sz+c5czunTxzmnHN5yOeWcc65POTF\n",
"3Tnn8pAXd+ecy0Ne3J1zLg95cXfOuTzkxd055/KQF3fnnMtDXtydcy4PeXF3zrk89P+s0mbHL+3g\n",
"UgAAAABJRU5ErkJggg==\n"
],
"text/plain": [
"<matplotlib.figure.Figure at 0x7fc5b8c06050>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"semilogx(freqs,appAna_p,'bo--',freqs,appSyn_p,'gx:')\n",
"gca().invert_xaxis()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
+1 -1
View File
@@ -1,7 +1,7 @@
from SimPEG import Survey, Problem, Utils, Models, np, sp, SolverLU as SimpegSolver
from scipy.constants import mu_0
from SurveyMT import SurveyMT, FieldsMT
import multiprocessing
import multiprocessing, sys
def omega(freq):
"""Change frequency to angular frequency, omega"""
@@ -0,0 +1,48 @@
import unittest
from SimPEG import *
import simpegMT as MT
TOL = 1e-6
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
return app_res, app_phs
def appResNorm(sigmaHalf):
nFreq = 26
m1d = Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C'])
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC[:]>200] = 1e-8
# Calculate the analytic fields
freqs = np.logspace(4,-4,nFreq)
Z = []
for freq in freqs:
Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
Z.append((Ed + Eu)/(Hd + Hu))
Zarr = np.concatenate(Z)
app_r, app_p = appResPhs(freqs,Zarr)
return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf)
class TestAnalytics(unittest.TestCase):
def setUp(self):
pass
def test_appRes2en1(self):self.assertLess(appResNorm(2e-1), TOL)
def test_appRes2en2(self):self.assertLess(appResNorm(2e-2), TOL)
def test_appRes2en3(self):self.assertLess(appResNorm(2e-3), TOL)
def test_appRes2en4(self):self.assertLess(appResNorm(2e-4), TOL)
def test_appRes2en5(self):self.assertLess(appResNorm(2e-5), TOL)
def test_appRes2en6(self):self.assertLess(appResNorm(2e-6), TOL)
if __name__ == '__main__':
unittest.main()
+1 -1
View File
@@ -21,7 +21,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
# Set the boundary conditions
Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx)
Etot = Ed + Eu
Etot = (Ed + Eu).conj()
if sourceAmp is not None:
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: need to use conjugate of the analytic solution. It is derived with e^iwt