sp_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#  ------------------------------------------------------------------------------------\n",
    "#  ------------------------------------------------------------------------------------\n",
    "#\n",
    "#      ,o888888o.        ,o888888o.       d888888o.   8 888888888o  `8.`8888.      ,8' \n",
    "#     8888     `88.   . 8888     `88.   .`8888:' `88. 8 8888    `88. `8.`8888.    ,8'  \n",
    "#  ,8 8888       `8. ,8 8888       `8b  8.`8888.   Y8 8 8888     `88  `8.`8888.  ,8'   \n",
    "#  88 8888           88 8888        `8b `8.`8888.     8 8888     ,88   `8.`8888.,8'    \n",
    "#  88 8888           88 8888         88  `8.`8888.    8 8888.   ,88'    `8.`88888'     \n",
    "#  88 8888           88 8888         88   `8.`8888.   8 888888888P'      `8. 8888      \n",
    "#  88 8888           88 8888        ,8P    `8.`8888.  8 8888              `8 8888      \n",
    "#  `8 8888       .8' `8 8888       ,8P 8b   `8.`8888. 8 8888               8 8888      \n",
    "#     8888     ,88'   ` 8888     ,88'  `8b.  ;8.`8888 8 8888               8 8888      \n",
    "#      `8888888P'        `8888888P'     `Y8888P ,88P' 8 8888               8 8888      \n",
    "#\n",
    "#  ------------------------------------------------------------------------------------\n",
    "#  ------------------------------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# import packages \n",
    "# --------------------\n",
    "\n",
    "import h5py\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "\n",
    "#import matplotlib.cbook as cbook\n",
    "\n",
    "from numpy import cos, sin, sqrt, exp\n",
    "from scipy.optimize import curve_fit, leastsq\n",
    "\n",
    "from cospy.interface import * \n",
    "%matplotlib inline  \n",
    "\n",
    "from matplotlib import cm\n",
    "\n",
    "import matplotlib.pyplot as plt \n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from sympy.interactive import printing\n",
    "printing.init_printing(use_latex=True)\n",
    "\n",
    "from IPython.display import display, HTML\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "from IPython.html import widgets\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "\n",
    "time_step = '1'\n",
    "\n",
    "# --------------------\n",
    "# exact solution\n",
    "# --------------------\n",
    "\n",
    "zi = 1\n",
    "w  = 0.0005\n",
    "\n",
    "def solution(x, y, z, t):\n",
    "\n",
    "    h      = 0.38\n",
    "    hh     = 0.125657*1.75\n",
    "    sigma2 = 1**2\n",
    "    r1     = sqrt(x**2 + y**2 + z**2) \n",
    "    r2     = x**2 + y**2 + z**2\n",
    "\n",
    "    return h*exp(-hh*r2)*cos(w*t)\n",
    "\n",
    "# --------------------\n",
    "# data\n",
    "# --------------------\n",
    "\n",
    "file_list = list()\n",
    "\n",
    "file_list.append('/Users/fuentes/workspace/sp_x0_32_t1/output.h5')\n",
    "file_list.append('/Users/fuentes/workspace/sp_x0_64_t1/output.h5')\n",
    "file_list.append('/Users/fuentes/workspace/sp_x0_128_t1/output.h5')\n",
    "\n",
    "# --------------------\n",
    "# computing abs-error\n",
    "# --------------------\n",
    "\n",
    "E = list()\n",
    "\n",
    "for fl in file_list:\n",
    "    \n",
    "    h5file = h5py.File(fl)\n",
    "    \n",
    "    u = h5file['psiR/' + time_step]\n",
    "    t = u.attrs['time']\n",
    "    \n",
    "    xx = np.array(h5file['coordinates/x'])\n",
    "    yy = np.array(h5file['coordinates/y'])\n",
    "    zz = np.array(h5file['coordinates/z'])\n",
    "    \n",
    "    Udiff = np.abs(u - solution(xx,yy,zz,1))\n",
    "    E.append(np.linalg.norm(Udiff.flatten()))\n",
    "    h5file.close()\n",
    "    \n",
    "# --------------------\n",
    "# least squares\n",
    "# --------------------\n",
    "\n",
    "y = E\n",
    "x = np.array([1/float(2**n) for n in range(7,4,-1)])\n",
    "\n",
    "funcLine = lambda tpl, x: tpl[0]*x    + tpl[1]\n",
    "funcQuad = lambda tpl, x: tpl[0]*x**2 + tpl[1]*x + tpl[2]\n",
    "\n",
    "func = funcLine\n",
    "ErrorFunc = lambda tpl, x, y: func(tpl, x) - y\n",
    "tplInitial1 = (1.0, 2.0) # arbitrary\n",
    "\n",
    "tplFinal1,success = leastsq(ErrorFunc, tplInitial1[:], args = (x, y))\n",
    "\n",
    "# --------------------\n",
    "# linear fit\n",
    "# --------------------\n",
    "\n",
    "x1 = np.linspace(x.min(), x.max(), 50)\n",
    "y1 = func(tplFinal1, x1)\n",
    "\n",
    "# --------------------\n",
    "# quadratic fit\n",
    "# --------------------\n",
    "\n",
    "func = funcQuad\n",
    "tplInitial2 = (1.0, 2.0, 3.0) # arbitrary\n",
    "\n",
    "tplFinal2, success = leastsq(ErrorFunc, tplInitial2[:], args = (x, y))\n",
    "\n",
    "x2 = x1\n",
    "y2 = func(tplFinal2, x2)\n",
    "\n",
    "# --------------------\n",
    "# maps\n",
    "# --------------------\n",
    "\n",
    "fig = plt.figure()\n",
    "map = fig.add_subplot(111)\n",
    "\n",
    "map.set_title('Convergence test')\n",
    "map.set_xlabel('$1/N$')\n",
    "map.set_ylabel('$|E|$')\n",
    "map.grid(True)\n",
    "\n",
    "map.loglog(x,   y, 'o' )\n",
    "map.loglog(x2, y2, '--')\n",
    "\n",
    "plt.xlim([1/512,1])\n",
    "plt.ylim([-100,1000])\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# maps wavefunction\n",
    "# --------------------\n",
    "\n",
    "file_sp ='/Users/fuentes/workspace/sp_x4_64g_t12/output.h5'\n",
    "data_sp = h5py.File(file_sp)\n",
    "\n",
    "# --------------------\n",
    "\n",
    "nx = data_sp['coordinates/x'][:,0,0].shape[0]\n",
    "x_slider = widgets.IntSlider(min=0,max=nx-1,value=25,description='X index')\n",
    "\n",
    "ny = data_sp['coordinates/y'][0,:,0].shape[0]\n",
    "y_slider = widgets.IntSlider(min=0,max=ny-1,value=55,description='Y index')\n",
    "\n",
    "nz = data_sp['coordinates/z'][0,0,:].shape[0]\n",
    "z_slider = widgets.IntSlider(min=0,max=nz-1,value=nz/8,description='Z index')\n",
    "\n",
    "# --------------------\n",
    "\n",
    "def plot_timestep(xi, yi, zi):\n",
    "\n",
    "\n",
    "    fig = plt.figure()\n",
    "    fig.clf()\n",
    "    ax = Axes3D(fig)\n",
    "\n",
    "    # iter_time = str(t)\n",
    "\n",
    "    X = data_sp['coordinates/x'][:,30:100,zi]\n",
    "    Y = data_sp['coordinates/y'][:,30:100,zi]\n",
    "    U = data_sp['psiR/1'][:,30:100,zi]\n",
    "    \n",
    "    ax.plot_wireframe(X, Y, U, color = 'Black',linewidth=.125)\n",
    "    ax.w_xaxis.set_pane_color((1.0,1.0,1.0))\n",
    "    ax.w_yaxis.set_pane_color((1.0,1.0,1.0))\n",
    "    ax.view_init(elev=xi,azim=yi)  \n",
    "    ax.grid(b=False)\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$y$')\n",
    "    ax.set_zlabel('Re $\\Psi$')\n",
    "    \n",
    "    #ax.set_zlim3d([-.75,.75])\n",
    "    \n",
    "    fig.savefig('ho_0.png', dpi=200)\n",
    "\n",
    "interact(plot_timestep, xi = x_slider, yi = y_slider, zi = z_slider)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAESCAYAAADXMlMiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9+P/XeyaZTFiCCrGsQgRkCYtFpcYiDIsCKmJo\nrYr9aivykw9g7fIBoQgMIOBWrQpCKVBFWtQPYgti6wKfFIFRoBaRLWxBCQmLiiyWZLKc3x+ZzCcJ\niSSZ5c7yfj4e95HMnTP3vu9k5r5zz7nnHDHGoJRSStWXzeoAlFJKRTdNJEoppQKiiUQppVRANJEo\npZQKiCYSpZRSAdFEopRSKiCaSCKUiLitjiHc9Jjjgx5z7BHtRxKZRMQYY8TqOMJJjzk+6DHHHr0i\nUUopFRBNJEoppQIS81VbIhLbB6iUUiFS2+q4hFAHEgliPVkqpVSwidS+SUertpRSSgVEE4lSSqmA\naCJRSikVEE0kSimlAqKJRCkVNQoLC5kyZQrNmzfHZrMhIrrUc7Hb7Vx33XXk5uYG/HeJi7u2lFKx\nYfjw4TidTjZv3swVV1xBQoKewurL6/XyzDPPkJmZydatWwPaVlz0I4n1Y1TRyePxkJWVhcvlIiMj\nw+pwokJiYiJnzpwhOTnZ6lBigtfrJTk5mZKSkgueE5Fa9yPRRKKUBTweDwMGDKCgoACn08n69es1\nmdSC7+RmdRgxpab3tC6JRNtIlAojj8fD3LlzWbZsGYWFhQAUFBTgdrvxeDwWR6dU/egViVJh4vF4\nGDhwIOfPn6+xzObNm/XK5DvoFUnw6RWJUlEkKysLr9frf1x1CAq73U5WVlaYo1LB1K5dO773ve/x\n7bff+tctXrwYl8tVp+3079+f1NRUUlJS6NmzJ3/729/8z61du5Y+ffpwySWX0Lx5c0aPHs3Zs2eD\ndQj1oolEqTBxuVw4HA7/bavGmEp3HTkcjjqfcFTdud3ukG6/uLiY559/PqBtPP/88+Tn53PmzBkW\nLVrET3/6U/Lz8wE4ffo0jz32GHl5eezZs4fc3FwmTJgQjNDrTROJUmGSkZHB6tWrGTRoEDfeeCNj\nx45l27ZtDB8+nNTUVNatW6fVWmEwY8aMkG5/woQJPPPMM3zzzTf13kaPHj38/2SICEVFRRw5cgSA\nkSNHMmTIEBo0aMCll17K6NGj2bRpU1Biry+9CVupMCi/1ffLL79k+fLlpKam+p/761//yhtvvEFK\nSoreEhwDrr32WlwuF8888wyPP/54ped69OjBF198Ue3rRo4cyUsvveR/fNttt/HBBx9QWFjI4MGD\nufbaa6t93YYNG0hPTw/eAdSHMSaml7JDVMo6mzdvNk6n0wAmISHBbN68+YIyxcXFZvTo0f5yTqez\n2nLxrr7f5+nTpxvggmX69OlBja9t27bm/fffN5999plJSUkxJ06cMH/84x9Nv3796rU9r9dr3nnn\nHfPss89W+/x7771nLrnkEpOdnV3vmGt6T33ra3We1aotpUIsKyvLf6tvcXFxpQb18vp6u93ON998\nQ0FBAVA2FIg2vAeP2+2u+M+l//dQtZd069aN2267jSeeeCKg7SQmJjJ06FDeffddVq9eXem5jz76\niJEjR7Jy5UquuuqqgPYTsNpmnGhd0CsSZaHyK4yqi9PpNMb833+DFyunygTj+xzKc0L5FYkxxuzf\nv980btzYuN1u/xVJ165dTcOGDatdHnrooRq3O3DgwEpXJZ988olJTU01q1evDjjmmt4P9IpEqchw\n6NAhRo4c6W84TUpK4t577yUnJ6facna7/TvLqcBNnz49LPvp0KEDd911Fy+88IJ/3a5duzh37ly1\ny8KFCwHYu3cvf//73zl//jxFRUUsX76cDRs20K9fPwB27tzJkCFDePHFFxk2bFhYjuViwp5IRGSI\niGSLyAERmVTN831F5BMRKRaRH1dYf7WIeERkl4jsEJG7whu5UnXXokULUlJSKCkpISkpiaKiIrKz\ns2nRooW/H4mI0LJlS7KzsyktLSUxMZGioiJSUlJo3ry5xUcQe0J9+29F06ZNq9SnpDaMr8rt8ssv\nJzU1leeff57XX3+dXr16AfC73/2OkydPMmrUKBo1akSjRo0sb2wPa892EbED+4CbgFxgK3CPMWZ3\nhTLtgBTgv4HVxpiVvvVXUXaptV9EWgL/AroYY77zHjvt2a6slpmZybFjx1i4cCGLFi0iPz+fVatW\nAZV7FY8YMYLmzZuTmJhIcXFxpXKqjPZsD75g9GwPdyLJANzGmMG+x5MBjDFzqyn7MvB2eSKp5vlP\ngR8bY/ZfZJ+aSJSlPv30U/Ly8hg6dOgFz1X3JX7iiSd49NFHL+j5rjSRhEI0DpHSCjhS4XGub12d\niEhvwAEcDFJcSoWEx+NhypQpNQ57Xl19fffu3Vm+fDlz587VgRxVVAh3h8Tqslud/r0QkRbAq8D9\nxpjSap53A+FpTVPqO5QP0lhQUMD69eur7bleXX19w4YNyczMpLS0FIfDoT3elWVEpOL5eYYxxl1d\nuXBfkeQCbSo8bg3k1fbFIpICrAUeM8Z8VF0ZY4zbGCPlS0DRKhWA8v4jxhi8Xm+t+4V4PB6Ki4sp\nKSmp0+uUCraK59KakgiEP5FsBTqKSJqIOIC7gdUXeQ0AvvJvAcuMMf8TwhiVCgqXy4Xdbsdms9Vp\nQEaXy0VCQkKdX6eUVcI+H4mI3AL8HrADS40xs0VkJrDNGLNaRK6jLGFcChQAx4wx6SLyU+BPwK4K\nm/uZMWb7Rfanje3KMuPGjaN169Z1Hjtr5cqVrFq1iocfflirtSrQxvbgi7q7tqygiURZ6amnnmLi\nxIl1fp0xhqeeeopHH300BFFFL00kwReNd20pFTdOnDjB5ZdfXq/X6q2/KppoIlEqRDweT0DVUk6n\n0z+Io1KRTBOJUiGyb9++gEZlveaaa/jXv/4VxIhUqAVrql0omyUxLS2Nhg0b0qVLF/bt23dBmZ//\n/OeICAcOHAgk7IBpIlEqRIwxAVVRXXPNNWzbti2IEan8/Hz69evHsWPHQraPYEy1u3jxYpYsWcLa\ntWs5d+4cb7/9Ns2aNatUZuPGjRw8GBl9sjWRKBUCRUVFleZjr4/k5GSt2gqyWbNmsXHjRmbOnBmy\nfQQ61W5paSkzZszgueeeo2vXrogI7du357LLLvOXKS4u5uGHH2bevHnBCjsgmkiUCoHt27dz9dVX\nB2VbepdS4JKTkxERFixYQGlpKQsWLEBEahy6JhAVp9qtqkePHlxyySXVLmPHjgUgNzeX3Nxcdu7c\nSZs2bUhLS2P69OmUlv7fQB7PPfccffv2pUePHkGPvz40kSgVZB6Ph1mzZgUlARQVFTFx4kQdcytA\n5fO9NGjQAIAGDRqEdL6XmTNn8uKLL3Ly5MlK63fs2ME333xT7VI+X3tubi4A7733Hp999hn/+7//\ny4oVK1iyZAkAR44c4Q9/+ENIr6rqShOJUkHk8XgYMGAAa9as4bbbbgsoAXg8Hh5//HGeeeYZBgwY\noMkkAOXzwhQUFPjvhgvlfC+BTLVbfpU0ceJELrnkEtq1a8dDDz3EO++8A8Avf/lLpk2bRpMmTYIa\ncyA0kSgVRBXnZw903vWsrCy8Xm9QtqXg+PHjjBkzho8++ogxY8aEtMEdYMaMGfzxj3/k6NGj/nXp\n6en+yaiqLmPGjAGgU6dOOByOGm/UWLduHRMmTKB58+b+RJiRkcFf/vKXkB7Pd9Ge7UoFSU2N406n\nk/Pnz1u2rVgS6T3b27Vrx+LFixk0aBAAo0ePZtWqVXTv3r1O/wjcd999fP3116xYsYLTp08zaNAg\nJkyYwKhRozhx4kSl9pIWLVrg8Xjo2bNnvdp8tGe7UhGkvB4+MTERCGze9fJtJSUlAeBwOHQO9yhU\nn6l2AebNm0ejRo1o2bIlGRkZjBw5kgceeACAyy+/3H81Un5F0qxZs5DcOFBb4Z6PRKmYVV4PX1xc\nDBDQvOvl2yoqKgp4Wyp8Dh8+XOlxmzZt6nULd0pKCq+99lqtykbCFZpekSgVRMePH+f6668HCLge\nvrxOH+Dmm28OeZ2+UvVmjInppewQlQq96dOnG8pm/Ky0TJ8+3dJtxRL9PgdfTe+pb32tzrPa2K5U\nEHm9XhYsWMAvf/nLoFU5iAhPPfUUEyZMCMr2olmkN7ZHI21sVyrC7N69my5dulgdhlJhpYlEqSDa\nvn07PXv2ZPr06UHb5vTp07nsssv46quvgrZNpYJJq7aUCqL6zoh4Mf/61784ffo0AwYMCPq2o4lW\nbQWfVm0pFWFCNbNheno6u3btCsm2lQqUJhKlgmTz5s2sX78+JGNiOZ1OsrOzmTt3ro65pSKOVm0p\nFQQej4f+/fvj9XpxOp2sW7cuoGl2q9t+3759McbgcDiCvv1ooVVbwadVW0pFiKysLIqKijDG4PV6\ngz7AYlZWFiUlJZSUlIRk+yo4gjXV7tSpU+nevTsJCQm43e5Kz61du5Y+ffpwySWX0Lx5c0aPHs3Z\ns2f9z3/99dfcddddNGvWjGbNmnHvvfdy5syZQA7rosKeSERkiIhki8gBEZlUzfN9ReQTESkWkR9X\nee5+EdnvW+4PX9RKfTeXy4XNZsNms+FwOOo1R/fFtp+YmBiy7avgCcZUux06dOCpp57i1ltvveC5\n06dP89hjj5GXl8eePXvIzc2t1Mfoscce49SpUxw6dIiDBw9y/PjxC5JRsIU1kYiIHZgPDAW6AveI\nSNcqxb4Afgb8pcprLwOmAz8AegPTReTSUMesVG1kZGTw0EMP8fjjj4ek2ikjI4MVK1Zw1113xW21\nVjB4PJ6QtzMFOtUuwP3338/QoUNp3LjxBc+NHDmSIUOG0KBBAy699FJGjx7Npk2b/M/n5ORwxx13\nkJKSQpMmTcjMzAz5jRrhviLpDRwwxhwyxniB14DhFQsYYw4bY3YApVVeOxh43xjztTHmFPA+MCQc\nQStVG23btmXy5MkhO8nfcccd9OrVS5NIPXk8HgYOHMjUqVMZOHBgyJJJoFPt1tWGDRtIT0/3Px43\nbhxvv/02p06d4tSpU7z55psMHTq03sdTG+Ee/bcVcKTC41zKrjDq+9pWQYpLqYCUlpaG7Nbfcjab\nTRuaA1A+UVjFdqZQJeWZM2fywx/+kEceeaTS+h07dgR1P++//z6vvPIKH3/8sX9dr1698Hq9NG3a\nFICBAwfWO0nVVrivSKr7ptX2mxHIa5UKqcOHD5OWlmZ1GOo7uFwuHA4Hdrs95O1MgUy1W1sfffQR\nI0eOZOXKlVx11VX+9XfeeSdXXXUVZ8+e5cyZM7Rv356f/vSnIYsDwp9IcoE2FR63BvKC+VoRcYuI\nKV/qHalSdbBr1y66dq3a3Bd8iYmJ/jlKVN1kZGSwbt06Zs2aFZZ2pvpOtVsb//73v7n99ttZunQp\nAwcOrPTcp59+ykMPPUTDhg392y2f772uKp5LRcRdU7lwV21tBTqKSBpwFLgbGFnL174LzKnQwH4z\nMLlqIWOMG3CXP9ZkosJh//79DBkS+ia7Dh06sH///rAkrViUkZERtjamDh06cNddd/HCCy/QvXt3\ngFo3ehcVFVFSUkJpaSnFxcUUFBSQmJiI3W5n586dDBkyhBdffJFhw4Zd8NrrrruOxYsX89RTTwGw\naNEievbsWa9jiMh+JMaYYmA8ZUlhD/CGMWaXiMwUkdsBROQ6EckF7gT+ICK7fK/9GphFWTLaCsz0\nrVPKcsXFxf4pdkMpPT2d3bt3h3w/KjjqO9Xu6NGjSU5OZsWKFcyePZvk5GReffVVAH73u99x8uRJ\nRo0a5b+aqdjYvnTpUg4fPkzr1q1p1aoVhw4d4uWXXw7WIVVLe7YrFQShGqyxqtLSUp599ln++7//\nO+T7ikTasz34tGe7UhGgtLQUmy08XyW9c0tFIk0kSgXo8OHDtGvXzuowlLKMJhKlArR79+5KddSh\nlpCQoHduqYiiiUSpAL3//vu8+eabYRvevbi4mIkTJ+pw8ipiaGO7UgEoH969uLgYp9PJ+vXrQ3p7\nqcfjweVy+YerD/X+Io02tgefNrYrZbGsrCyKi4sBKCwsDPnw7uXDfIRrf0rVhl6RKFVPycnJFBQU\nXLDe6XRy/vz5qN9fJNIrkuDTKxKlLHTo0CEyMzNJSCgbICIpKYl7772XnJyckO1v5MiRJCUlhWV/\nStWWJhKl6qlFixb+WQudTidFRUWkpKTQvHnzkO0vJSWFoqIiEhIS8Hq9Id2fUrUV7rG2lIopX3zx\nBQ888AAPP/wwixYtIj8/P6T7O378OGPGjKFVq1bs27ePY8eOhXR/KnK9/PLLLF68mI0bN1odiiYS\npQIxcuRI/zSn8+fPD/n+Vq1aBcC2bdvIyMigf//+Id+nUhejVVtKRaGrrrqK7Oxsq8NQVeTl5fGj\nH/2I1NRU0tLSeOGFFwC45ZZb+M1vfuMvd9ddd/HAAw8AcPDgQQYMGEDTpk1p1qwZ9957b6Vpeo8c\nOcKIESNITU2ladOmjB8/nj179jBmzBg8Hg+NGjXikksuAWDt2rV8//vfJyUlhTZt2oR8rvZyekWi\nVBRKSUnh7NmzVodhuTNnzrBz586Q7qNbt26kpKRctFxpaSnDhg1j+PDhrFixgtzcXAYNGkSnTp1Y\nunQpPXr04NZbbyU/P5+tW7fy6aefAmCMYfLkyfTt25czZ87wox/9CLfbze9//3tKSkq47bbbGDBg\nAK+++ip2u51t27bRpUsXFi5ceEHVVsOGDVm2bBnp6ens3LmTm266iauvvpo77rgjZO8PaCJRqt4K\nCwtxOBxWh6EixNatWzl58iTTpk0D4Morr2T06NG89tprDB48mIULF3L//fdz/vx5/vrXv9K4cWOg\nbN6SDh06AJCamsqvf/1rZsyYAcCWLVvIy8vj6aef9t8d2KdPnxpjqDjrY48ePbjnnnv45z//qYlE\nqUh18OBB/wlAWSMlJYUbbrjB6jAA+Pzzz8nLy/NXMwGUlJRw4403AnDbbbcxfvx4OnXqVCkZnDhx\ngl/84hd8+OGHnD17ltLSUi69tGz+viNHjtC2bVt/ErmYjz/+mEmTJrFz5068Xi+FhYXceeedQTzK\n6mkbiVL1lJ2dTadOnSzbf3Jyctx0RIwGbdq0IS0tjW+++ca/nD171j/N7ZQpU+jSpQv5+fmsWLHC\n/7rJkycjIuzYsYMzZ86wfPlyfwfBNm3a8MUXX/hHT6hI5MK+giNHjuT222/nyJEjnD59mjFjxoSl\nA6cmEqXqKScnx9Lh4zt27Mj+/fst27+qrHfv3qSkpPDkk09y/vx5SkpK2LlzJ1u3bmXDhg386U9/\nYtmyZSxbtoyHH37YP5f72bNn/Q3mR48e5emnn660zRYtWjBp0iS+/fZbCgoK2LRpEwDf+973yM3N\n9Q+ZU76tyy67DKfTyZYtW/jLX/4SlmPXRKJUPZWUlNS6yiEUOnfuzN69ey3bv6rMbrezZs0atm/f\nTlpaGs2aNePBBx8kPz+f++67j3nz5tGqVSv69OnDqFGj+PnPf44xhunTp/PJJ5/QpEkTbr31VkaM\nGHHBNg8cOMAVV1xB69atef311wEYMGAA6enpNG/enGbNmgHw0ksvMW3aNBo3bszMmTP5yU9+EpZj\n17G2lKoHj8fDzJkzmTZtmmWj727atInZs2czderUuBkBWMfaCr5gjLWliUSpOvJ4PAwYMIDCwkKc\nTifr1q0L+4nc4/EwcOBACgoKLIvBCppIgk8HbVTKAuVDuRtj8Hq9lgzlHgkxKFVOE4lSdeRyuUhI\nSMBms+FwOCrdux/OGBwOByJCYmKiJTEoVU6rtpSqh8mTJ5OYmMjQoUMtq1LyeDw8//zz3HPPPQwf\nPtySGMJNq7aCLxhVW2G/5UREhgDPA3ZgsTHmiSrPJwHLgGuAr4C7jDGHRSQRWAz0oizuZcaYuWEN\nXimfyy67zD9Yo1UyMjJISkri66+/tjSOcEpISOD8+fMkJydbHUpM8Hq92GyBV0yFtWpLROzAfGAo\n0BW4R0S6Vik2CjhljOkAPAc86Vt/J5BkjOlOWZJ5SETahSNupSLVVVddxb59+6wOI2wGDhzI3Xff\nzcGDB6vtpKdqz+v18swzz9CrV6+AtxXuNpLewAFjzCFjjBd4Dah6TT4ceMX3+0pgoJR14TRAQxFJ\nAJIBL3AmPGErFZkaNWrEt99+a3UYYfO3v/2Nbt260adPH5KSkhARXeq5JCcn89Zbb/HWW28F/HcJ\nd9VWK+BIhce5wA9qKmOMKRaR00BTypLKcCAfaAD8yhgTP9f0KmJ4vV4SExOtDiMuJSUlMXv2bGbP\nnm11KKqCcF+RVNdwU7WVp6YyvYESoCWQBvxGRK68YAcibhEx5UugAStV1cGDB2nfvr3VYSgVchXP\npSLirqlcuBNJLtCmwuPWQF5NZXzVWE2Ar4GRwD+MMUXGmBPAJuDaqjswxriNMVK+hOAYVJyzerDG\nqho0aBBX1VsqfCqeS40x7prKhTuRbAU6ikiaiDiAu4HVVcqsBu73/f5jYL3v/t0vgAFSpiFwPaAD\nDamwy8nJIS0tzeow/Dp16hRXDe4q8oQ1kRhjioHxwLvAHuANY8wuEZkpIrf7ii0BmorIAeDXwCTf\n+vlAI2AnZQnpT8aYHeGMXymA4uLiiGoj6dSpkw7eqCylHRKVqqObbrqJ999/3+ow/EpLS3n22Wc5\nd+5c2OboVrFPx9pSKoQ++OADq0OoxGazYYzxT8+qVLhpIlGqDr788kurQ1Aq4mgiUaoW3G43IkJq\naiqAv1OX1VVJ5XFNnDgxouJS8UXbSJSqJY/HQ58+fSgtLSUhIYENGzZExBwgkRqXim51aSPRRKJU\nLSQnJ1NQUHDBeqfTyfnz5y2IqEykxqWinza2KxVkhw4d4qabbqq07uabbyYnJ8eiiMpEalwqvmgi\nUaoWWrRo4R8WJSEhARGhffv2NG/ePCLiEhFsNlvExKXiS9jnI1EqWuXl5fHDH/6Q+fPns2jRIvLz\n860OCYDjx4/zX//1XzgcDrxeb8TEpeKHtpEoVUu7d+8mJyeHW2+91epQqvX0009bPtmWih0haSMR\nkRUi8qiIDBWR1PqHp1R02rt3L507d7Y6jBo1atSIs2fPWh2GikN1aSNZCPwH+BHwvoi8JSKNQxOW\nUpEnJyeHdu3aWR1GjTp16kR2drbVYag4VJdEcgXwT2CMMeZqyiaa0jEZVNwoLS3FbrdbHUaNOnfu\nrIlEWaIuje1XUTZvelcROQvsAPqJyFpghzHmZCgCVErVTosWLcjLqzq9j1KhV5dE8kJ5shCRSymb\nVGoocC/wBHBd8MNTKjIYYxCJ7HnSIj0+Fbvqkkje8yWQPUA20BlYYYx5JCSRKRVBjh07pn0zlKpB\nrROJMeb7ImKnLIF0Aj4A/hGqwJSKJJF+x1Y5m81GSUlJRLflqNhTp57txpgSY8wuY8wqY8zbvhkP\nlYp5e/fujah52mvStm1bvvjiC6vDUHFGh0hR6iI8Hg+rV69m586dVodyUQUFBcyYMQOPx2N1KCqO\naM92pb6Dx+Nh4MCBFBQU4HQ6WbduXcQO0e7xeBgwYACFhYURH6uKfDr6r1JBkpWVhdfrxRiD1+sl\nKyvL6pBqlJWVRVFRUVTEqmJLXYZIERH5qYhM8z2+QkR6hy40pazncrlITExERHA4HLhcLqtDqpHL\n5cLhcERFrCq21LpqS0QWAKXAAGNMF9+twO8ZYyK6/4hWbalAvfzyy3z44Yc8+OCDEV9V5PF4mDlz\nJtOmTYv4WFVkC1XV1g+MMeOAAgBjzCnAUY/ghohItogcEJFJ1TyfJCKv+57/WETaVXiuh4h4RGSX\niHwmIs667l+pukpKSmLGjBlRcWLOyMjgjjvuoGvXrlaHouJIXRJJka8fiQHwjQBcWped+V4/n7Ie\n8V2Be0Sk6id+FHDKGNMBeA540vfaBGA5ZWN9pQMuoKgu+1eqPnJzc2nVqpXVYdSajrmlwq0uieQF\n4C3gchGZDWwE5tRxf72BA8aYQ8YYL/AaMLxKmeHAK77fVwIDpWzsh5spG9PrUwBjzFfGmJI67l+p\neomm4Uc6d+7M3r17rQ5DxZG69Gz/s4j8CxjoW3Un0L2O+2sFHKnwOBf4QU1ljDHFInIaaErZoJFG\nRN4FUoHXjDFP1XH/SsW8yy+/nOPHj1sdhoojF70iEZEUEZksIvMoG0r+Jd/r1gA/qeP+qvu3rmpL\neE1lEoA+lA0S2QfIFJGBVQuKiFtETPlSx/iUqqS0tDSqrkYguq6eVGSreC4VEXdN5WpzRfIqcArw\nAA8Cv6XsZD/cGLO9jnHlAm0qPG4NVB33urxMrq9dpAnwtW/9P40xXwKIyDtAL2BdxRcbY9yAu/yx\nJhMViC+++IK2bdtaHYZSlgjmXVtXGmN+Zoz5A3APZY3kg+uRRAC2Ah1FJE1EHMDdwOoqZVYD9/t+\n/zGw3nf/7rtADxFp4Esw/YDd9YhBqVqLlsEaq7Lb7RQX61B4Kjxqk0j8d0b5GrdzjTEF9dmZb5DH\n8ZQlhT3AG8aYXSIyU0Ru9xVbAjQVkQPAr4FJvteeAp6lLBltBz4xxqytTxxK1da+ffvo2LGj1WHU\nWVpaGocPH7Y6DBUnalO11VNEzvh+FyDZ91gAY4xJqcsOjTHvAO9UWTetwu8FlDXkV/fa5ZTdAqxU\nWJSPWxVtyudv79Chg9WhqDhw0SsSY4zdGJPiWxobYxIq/F6nJKJUNMnPz2fBggUcO3bM6lDqrH37\n9nzyySf069cvKuNX0UUHbVSqBuPGjSMnJ4exY8daHUqdJSUl8frrr/Phhx9GZfwquugw8kpVkZyc\nTEHBhc2ATqeT8+fPWxBR3UR7/Coy6DDySgXg0KFD9OzZs9K6nj17kpOTY1FEdRPt8avoo4lEqSpa\ntGhBWlpapXVXXnklzZs3tyiiuon2+FX0qfUQKUrFExGhW7du9O/fn6NHjxJt1aMiQv/+/WnSpAk2\nmy3q4lfRRdtIlKrB008/zYQJE6wOo97+85//sHTpUsaPH291KCoKaRuJUooGDRpo47oKC00kSlXj\nq6++ommagVaSAAAWzUlEQVTTplaHoVRU0ESiVDV27dpFenq61WEEhVbtqlDTRKJUNXbt2hUT09W2\nbt2ao0ePWh2GinGaSJSqxrlz52jcuLHVYQQsPT2dXbt2WR2GinGaSJSKYeWDNyoVSppIlKoiltoU\nkpKSKCwstDoMFeM0kShVxdq1a9m6dSsej8fqUILi888/Z+7cuTFzPCryaIdEpSrweDy4XC6Ki4tJ\nSkpi3bp1ZGRkWB1WvXk8Hvr27UtpaWlMHI8KH+2QqFQ9ZWVlUVRURGlpKV6vl6ysLKtDCkhWVhal\npaUxczwqMmkiUaoCl8uF3W7HbrfjcDhwuVxWhxQQl8tFYmIiNpstJo5HRSat2lKqivHjx9OqVStc\nLldMVANt3LiROXPmMHXq1Jg4HhUedana0tF/laqgqKiI9u3b86tf/crqUIKmT58+9O/fX5OIChmt\n2lLKx+12k52dTefOna0OJWTcbrfVIagYpIlEKZ8ZM2awfft2rr76aqtDCbomTZpw+vRpZsyYYXUo\nKgaFPZGIyBARyRaRAyIyqZrnk0Tkdd/zH4tIuyrPXyEi50Tkv8MVs4of+fn5MTmTYI8ePdixY4fV\nYagYFdZEIiJ2YD4wFOgK3CMiVUfGGwWcMsZ0AJ4Dnqzy/HPA30Mdq4oPbrcbEUGkrE1x4sSJ2Gy2\nmKoCcrvdZGRk0LdvXwD/8cbSMSprhfWuLRHJANzGmMG+x5MBjDFzK5R511fGIyIJwDEg1RhjROQO\n4IfAt8A5Y8wztdin3rWlLio/P5+WLVsydepUZs6caXU4ITF16lQef/zxmL3qUsEVyR0SWwFHKjzO\n9a2rtowxphg4DTQVkYbAo4BW8qqgmzVrFgCbN2+2OJLQ+eCDDwBiNlEq64T79t/qslvVy4WayswA\nnjPGnCuvhlAqUMnJyRQUFPgfr1u3DhHB6XTGzDS1VY9xwYIFLFiwIKaOUVkr3FckuUCbCo9bA3k1\nlfFVbTUBvgZ+ADwlIoeBXwK/FZHxVXcgIm4RMeVL8A9BxZJDhw5x0003VVp38803k5OTY1FEwRcP\nx6hCo+K5VETcNZUL9xXJVqCjiKQBR4G7gZFVyqwG7gc8wI+B9b5GjhvLC/gO6JwxZl7VHRhj3IC7\nQllNJqpGLVq04PTp05XWnT59OqbaEOLhGFVoRGTPdmNMse8q4l3ADiw1xuwSkZnANmPMamAJ8KqI\nHKDsSuTucMao4o/T6SQhIYGSkhLsdjtJSUlWhxR0TqcTu91OaWlpzB6jso6OtaUU8N5777FgwQIm\nTpwYs0OJeDweZs6cybRp02L2GFXw6FhbStWRzWZj6tSp9OrVy+pQQiYjI4NBgwbRu3dvq0NRMUaH\nSFEK2LFjB927d7c6jJDr0qULe/bssToMFWM0kSgFFBcXk5iYaHUYIXfdddexZcsWq8NQMUYTiVJA\nvPRNSk1N5auvvrI6DBVjNJGouHf06FFatmxpdRhhozefqGDTRKLi3pYtW+KqAVp7tKtg00Si4prH\n42HRokWcPHnS6lDC5vvf/z6vvvoqc+fOxePxWB2OigHaj0TFLY/Hw8CBAykoKMDpdLJu3bq46F+x\nfv16hgwZQmlpKQ6HI26OW9VNJI/+q1TEyMrKorCwEGMMXq+XrKwsq0MKi48//pji4mJKSkri6rhV\n6GgiUXHL5XKRmJiIzWbD4XDgcrmsDiksXC4Xdrsdu90eV8etQkertlRce+yxxygtLWXYsGFxVb0z\nZcoUAG677ba4Om5VezpEilK1lJKSwsSJE60OI+x+9rOfsXPnTk0iKii0akupONShQwf2799vdRgq\nRmgiUXEpPz+f66+/ngYNGlgdiiXKe/Ln5+fTr18/jh07ZnFEKpppIlFxadasWWzZsoWPP/7Y6lAs\n07hxY6ZOncrGjRt1HncVEG1sV3Gl6vzl5eKtt7e+D+pitB+JUjXQ+cvLHDp0iEGDBlVaF4/vgwoO\nTSQqruj85WVatGjBmTNnKq2Lx/dBBYfe/qvijtPpxGYr+x/KZrPF7fzl+j6oYNE2EhWXxo0bR6tW\nrejfv39c96VYsmQJHo+HUaNGxfX7oC6kHRKVuojWrVszefJkq8Ow3P3338+pU6c0iaiAaBuJiju5\nubm0adPG6jAiQkJCAsXFxVaHoaKcXpGouOHxeMjKyuLcuXM88sgjVocTMZo3b86xY8fIyckhKysL\nl8ulVyiqTsLeRiIiQ4DnATuw2BjzRJXnk4BlwDXAV8BdxpjDInIT8ATgALzABGPM+lrsT9tIlH/u\nEa/Xi4iwYcMGPVn65OXlsXDhQp555hm8Xq/OUaKACO5HIiJ2YD4wFOgK3CMiXasUGwWcMsZ0AJ4D\nnvSt/xIYZozpDtwPvBqeqFUseOyxx/B6vZSUlFBSUqJzcFTQsmVLtmzZ4n9/vF4vjz32mNVhqSgS\n7jaS3sABY8whY4wXeA0YXqXMcOAV3+8rgYFSdlnxb2NMnm/9LsDpu3pR6qLWr1+Pw+HAZrORmJio\nc3BU0aFDBxITE/1zlKxff9GLfaX8wt1G0go4UuFxLvCDmsoYY4pF5DTQlLIrknI/Av5tjCkMYawq\nRuTn5wPwxhtvMH/+fKZMmaLVNlX85Cc/oUePHnz11Vd069aN22+/nWPHjmkHRVUr4b4iqa6+rWoD\nxneWEZF0yqq7Hqp2ByJuETHlS70jVVHP7XYjIrRs2RKAYcOG8Y9//IMPPvjA4sgiT0ZGBitXruS3\nv/0tt99+O1DW+11EcLvd1ganLFPxXCoi7hrLhbMhWkQyALcxZrDv8WQAY8zcCmXe9ZXxiEgCcAxI\nNcYYEWkNrAd+bozZVMt9amN7nNKBCesmMTGx2luB9f2KTxHb2A5sBTqKSJqIOIC7gdVVyqymrDEd\n4MfAel8SuQRYC0yubRJR8a26ARoHDBigAxPWYMGCBdxwww2V1ulAjqo2wppIjDHFwHjgXWAP8IYx\nZpeIzBSR233FlgBNReQA8Gtgkm/9eKADMFVEtvuWy8MZv4ou1Q3Q+O2332q9fw3uvfdecnNzK63T\ngRxVbYS9Q6Ix5h3gnSrrplX4vQC4s5rXPQ48HvIAVUxxOp0kJCRQWloKoAMTfofk5GRsNht2ux0o\nq9rQ90vVhg7aqGKex+Nh9uzZ/PrXv2bAgAFWhxPR3n77bU6dOkVubq72cI9zkdxGolRYeDwe5s6d\ni8fjISMjgxtuuEGTSC3cfPPNHD9+3D+gZfl7qNR30bG2VMypOByKw+Fg0aJFpKenWx1WVHA4HHi9\nXjZv3sygQYN0yBRVK3pFomJOVlZWpeE+/vznPzN48GCrw4oaGRkZvPLKK5XeQx1SRn0XTSQqpuTn\n5/PGG2/4h/tITEykffv2OJ1Oq0OLGi6XC2MMDofDP2RKt27d6NevH8eOHbM6PBWBtGpLxZRZs2ax\nY8cOhg8fznXXXcf58+d54IEHrA4rqogI1157LcOGDWPnzp24XC5effVVNm7cyMyZM3nppZesDlFF\nGL1rS8WEmnqxJyQkUFRUZEFE0e3s2bMsWbKEyZMn6+gAcUrv2lJxp7pe7L1792b58uUWRRTdGjdu\nzPnz58nOzr7gfdXe7qoqTSQqJlTXiz0vL48777ygb6uqpczMTD766KML3lft7a6q0kSiYoLH4+HE\niRPY7Xb/0rBhQ2w2/YjXV+fOndm9ezdJSUkkJCT439fjx49r3xJVibaRqKjmdrsZPHiwv9+I3W7n\ngQceoKioiIULF5KQoPeTBGLXrl189tlntG3blmXLlvGnP/2J4uJiHA4H99xzD0uWLLE6RBUi2kai\n4kJ+fj4zZsxgzZo1labRtdvt3HLLLZpEgiA9PZ19+/ZxzTXXcMUVV1BcXOzvW7J06VK9HVgBmkhU\nFBs3bhwAmzZtqtTn4T//+Q+ZmZkWRxc77rvvPpYvX47L5fK/z+XGjh1rYWQqUui/bCrqVJ2AacOG\nDQDYbDbmzJlDjx49EKnVFbmqhXbt2vHll1/SqVMn/5VfubfeegsR0dus45xekaios2rVqgvuGurZ\nsye7d+/m7NmzOjhjCIwdO5YXXniBI0eO0LNnz0rPdejQgQkTJmgDfBzTxnYVVcoHZKzaGS4zM5Or\nr76aX/3qVzRu3Nii6GLbhx9+yNdff83LL7/MX//6V//68qouHdwxtmhju4pJa9as4ZZbbqGwsNC/\nrnXr1owYMYJjx47Rp08fTSIhdOONN7J37168Xi8jRoygd+/epKen+29yKCws5JZbbuHtt9+2OlQV\nZnpFoqKCx+Ohb9++/rYRm81GUlIS69ato3Hjxqxbt45HHnnE4ihjX0FBAdOnT8ftdpOcnOy/Qiws\nLPTPQpmQkMCGDRv0yiTK6RWJinr5+fm0a9eOt99+m8TERG644YZKDeylpaWUlJRwxRVX8Oabb/KL\nX/zCwmjjh9PpZOLEicyaNYuSkhIyMjIoKSnxJxGA4uJibrjhBhITE1mzZg3t2rXT24RjnCYSFZHG\njRvH559/zh133IEx5oIe6jfffDPbtm1j3rx5/Pa3v9W7tMKoadOmjB49mscff5ySkhIOHz58wXhc\nNpsNYwyZmZl8/vnneptwjNNEoiyXn59Pv379/FcfIsJbb70F4K9/r1htIiJceumlrFy5ErfbTWJi\nopXhx6W0tDTuu+8+pk6dSuPGjWnfvr3/NmD4vyvG8luFy28TTkpK8v+99SolhhhjYnopO0QVKfLy\n8kzfvn1Nfn6+f11mZqYBjN1uNyJigAsWETEjRoww27dvN0OGDDE9evQwpaWlFh6JMsaYU6dOmUmT\nJplBgwaZsWPHmu3bt5sRI0bU+HdMSEgwffv2NSJiMjMzjTHVfyaU9XznztqdZ2tbMFgLMATIBg4A\nk6p5Pgl43ff8x0C7Cs9N9q3PBgbXcn/BfG9VHVU9SZQnjczMTONwOKo92VRdunbtajZv3myOHDli\npkyZYjZt2mTxUamKSkpKzJIlS8xzzz1nCgoKjDHGbN682XTp0qVWf9/ypFOeWIzR5BIJIjaRAHbg\nIHAl4AA+BbpWKTMWWOj7/W7gdd/vXX3lk4A033bstdhnkN/e2JKXl2fatm3r/8JW9wWuWuZi5das\nWWPmzJljNm/ebDIzM42I1PgfanWLzWYz/fr1M4AZMWKEuemmm8ycOXPMggULTGFhYfjeHFUnR48e\nNbNnzzbz5883p0+fNpmZmWbEiBH+q82KSeO7lqpXLZs3bzZz5swxq1evvujnsK6fVU1UNYvkRJIB\nvFvh8WRgcpUy7wIZvt8TgC8BqVq2YrmL7LPeb2QoTrKRsK2qJ3rA9O3b18yZM6faL3Dfvn39/zFW\nXCci/tdV3Fb5SaM+S9euXc2mTZvMoUOHzK233mrmzJljli9frgkkihw/ftz8/ve/N3PnzjX/8z//\nYx599NELPkv1Wco/V1U/cxU/h9/1Wa3uM11xW8FMSqHcVnXlQpEsIzmR/BhYXOHx/wPmVSmzE2hd\n4fFBoBkwD/hphfVLgB/XYp/1ehONMSE5yVq9rUBP9KFYbDab6d27t+nWrZtJT083TzzxhHnttdfM\nV199Ve+/nYoM2dnZZuHCheaJJ54wTz75pOnevbu5/vrrzYABAyz/3FVdqn7XAklKwUxwVbdVn3NA\nfdQlkYS1Q6KI3ElZ28aDvsf/D+htjHm4QpldvjK5vscHgd7ATMBjjFnuW78EeMcY8+ZF9mnqeow1\nzf+tAufr5ET37t0pKCjg+uuv58Ybb+SWW26hZcuWehtvjCsqKuLw4cPk5eWxceNG5s2bR4MGDTh0\n6BBXXnklOTk5hPOcFE+cTucFQwt9l7p0SAz36L+5QJsKj1sDeTWUyRWRBKAJ8HUtX4uIuIHpgQQ5\nduxYnn322UA2ERVsNluljmRV2e32SiO91rSuum117dqVZs2a0axZM1q1asWOHTs4efIkixcv5oYb\nbmDHjh3BOxAVNRITE+nYsSMdO3akX79+TJkyBSg7aR08eBCPx8ODDz5IamoqPXr04OjRo3z55Zec\nPHmSPXv2VNpWMD+/5f/gxLLyaRfqwtemVW6GMcZdbbkwX5EkAPuAgcBRYCsw0hizq0KZcUB3Y8wY\nEbkbGGGM+YmIpAN/oezqpCWwDuhojLnwU1F5n3W+IoGyQQArDkxXVSAn2UjYVsUT/apVq+jXrx+f\nfPIJvXr1qvQF/vLLL0lNTeWf//wnI0aMqLSuR48e/gRRcVsjRozAGMOqVauqjdHtduN2u2s8BhV/\nLvaZGDFiBCLi/6xW/MyV/6NS/vmt6bNa3We64rbKpyOoTiR93+uz/czMzBq/jzWpyxVJ2MfaEpFb\ngN9TdgfXUmPMbBGZCWwzxqwWESfwKvB9yq5E7jbGHPK9dgrwAFAM/NIY8/da7K9eiaTqBzdYJ9lI\n2NbRo0crnehrc2Kv7clfk4QKpVB9Vsu/78FKSsFMcFWv6utzDviuf+xqUpdEEtbGdisWAmhsN8aY\n6dOnB6VMJG9LKVW9SP1uB3uf1SFSG9utoKP/KqVU3enov0oppcJGE4lSSqmAaCJRSikVEE0kSiml\nAqKJRCmlVEDC3bPdEjrshlJKhU7M3/4brXy3LcdVBtRjjg96zLFHq7aUUkoFRBOJUkqpgGgiiVwz\nrA7AAnrM8UGPOcZoG4lSSqmA6BWJUkqpgGgiUUopFRBNJBFGRO4UkV0iUioi11Z5brKIHBCRbBEZ\nbFWMoSQibhE5KiLbfcstVscUCiIyxPd3PCAik6yOJxxE5LCIfOb7u26zOp5QEJGlInJCRHZWWHeZ\niLwvIvt9Py+1MsZQ0EQSeXYCI4BK07WJSFfgbiAdGAK8JCL28IcXFs8ZY672Le9YHUyw+f5u84Gh\nQFfgHt/fNx709/1dr7140aj0MmXfz4omAeuMMR0pm9k15v5x0EQSYYwxe4wx2dU8NRx4zRhTaIzJ\nAQ5QNu2wij69gQPGmEPGGC/wGmV/XxXljDEbKJvZtaLhwCu+318B7ghrUGGgiSR6tAKOVHic61sX\ni8aLyA5fNUHMVQMQX3/Ligzwnoj8S0T+P6uDCaPvGWPyAXw/L7c4nqCLi7G2Io2IfAA0r+apKcaY\nv9X0smrWReW92991/MACYBZlxzYL+B3wQPiiC4uY+VvW0Q+NMXkicjnwvojs9f0Hr6KcJhILGGMG\n1eNluUCbCo9bA3nBiSi8anv8IvJH4O0Qh2OFmPlb1oUxJs/384SIvEVZFV88JJLjItLCGJMvIi2A\nE1YHFGxatRU9VgN3i0iSiKQBHYEtFscUdL4vWrlMym4+iDVbgY4ikiYiDspuolhtcUwhJSINRaRx\n+e/AzcTm37Y6q4H7fb/fD9RU6xC19IokwohIJvAikAqsFZHtxpjBxphdIvIGsBsoBsYZY0qsjDVE\nnhKRqymr6jkMPGRtOMFnjCkWkfHAu4AdWGqM2WVxWKH2PeAt35QOCcBfjDH/sDak4BORFYALaCYi\nucB04AngDREZBXwB3GldhKGhQ6QopZQKiFZtKaWUCogmEqWUUgHRRKKUUiogmkiUUkoFRBOJUkqp\ngGgiUUopFRBNJEoppQKiiUQppVRANJEoZQER+V8Rucn3++Mi8oLVMSlVXzpEilLWmA7M9I2E+33g\ndovjUaredIgUpSwiIv8EGgEuY8xZq+NRqr60akspC4hId6AFUKhJREU7TSRKhZlvqPw/UzYF67ci\nMtjikJQKiCYSpcJIRBoAq4DfGGP2UDYLpNvSoJQKkLaRKKWUCohekSillAqIJhKllFIB0USilFIq\nIJpIlFJKBUQTiVJKqYBoIlFKKRUQTSRKKaUCoolEKaVUQP5/3GDvrbJkVSoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109756190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# --------------------\n",
    "# maps wf profile\n",
    "# --------------------\n",
    "\n",
    "file_sp32 ='/Users/fuentes/workspace/sp_x0_32_t1/output.h5'\n",
    "data_sp32 = h5py.File(file_sp32)\n",
    "\n",
    "file_sp64 ='/Users/fuentes/workspace/sp_x0_64_t1/output.h5'\n",
    "data_sp64 = h5py.File(file_sp64)\n",
    "\n",
    "file_sp128 ='/Users/fuentes/workspace/sp_x0_128_t1/output.h5'\n",
    "data_sp128 = h5py.File(file_sp128)\n",
    "\n",
    "# --------------------\n",
    "\n",
    "zi = 0\n",
    "w  = 0.0005\n",
    "\n",
    "def solution(x, y, z, t):\n",
    "\n",
    "    h      = 0.12758447\n",
    "    hh     = 0.22\n",
    "    sigma2 = 1**2\n",
    "    r1     = sqrt(x**2 + y**2 + z**2) \n",
    "    r2     = x**2 + y**2 + z**2\n",
    "\n",
    "    return h*exp(-hh*r2)*cos(w*t)\n",
    "\n",
    "\n",
    "X32 = data_sp32['coordinates/x'][:,zi,zi]\n",
    "U32 = data_sp32['psiR/1'][:,zi,zi]\n",
    "\n",
    "X64 = data_sp64['coordinates/x'][:,zi,zi]\n",
    "U64 = data_sp64['psiR/1'][:,zi+1,zi+1]\n",
    "\n",
    "X128 = data_sp128['coordinates/x'][:,zi,zi]\n",
    "U128 = data_sp128['psiR/1'][:,zi+3,zi+3]\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.clf()\n",
    "\n",
    "map = fig.add_subplot(111)\n",
    "\n",
    "map.plot(X32,U32,'k+',linewidth=.15,label='N=32')\n",
    "map.plot(X64,U64,'k*',linewidth=.15,label='N=64')\n",
    "map.plot(X128,U128,'k.',linewidth=.15,label='N=128')\n",
    "map.plot(X128,solution(X128,0,0,zi),'k-',linewidth=.45,label='exacta')\n",
    "\n",
    "map.set_xlabel('$x$')\n",
    "map.set_ylabel('Re $\\psi$')\n",
    "\n",
    "legend = map.legend(loc='upper right', shadow=False)\n",
    "\n",
    "plt.show()\n",
    "fig.savefig('sp_x4_pf_conv.png', dpi=300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# maps potential\n",
    "# --------------------\n",
    "\n",
    "file_sp ='/Users/fuentes/workspace/sp_x1_64_t1/output.h5'\n",
    "data_sp = h5py.File(file_sp)\n",
    "\n",
    "# --------------------\n",
    "\n",
    "nx = data_sp['coordinates/x'][:,0,0].shape[0]\n",
    "x_slider = widgets.IntSlider(min=0,max=nx-1,value=nx/2,description='X index')\n",
    "\n",
    "ny = data_sp['coordinates/y'][0,:,0].shape[0]\n",
    "y_slider = widgets.IntSlider(min=0,max=ny-1,value=ny/2,description='Y index')\n",
    "\n",
    "nz = data_sp['coordinates/z'][0,0,:].shape[0]\n",
    "z_slider = widgets.IntSlider(min=0,max=nz-1,value=nz/8,description='Z index')\n",
    "\n",
    "#frames = len(data_sp['U/1'.keys())-1\n",
    "#time_slider = widgets.IntSlider(min=1,max=frames,value=1,description='Time')\n",
    "\n",
    "# --------------------\n",
    "\n",
    "def plot_timestep(xi, yi, zi):\n",
    "\n",
    "    fig = plt.figure()\n",
    "    fig.clf()\n",
    "    ax = Axes3D(fig)\n",
    "\n",
    "   # iter_time = str(t)\n",
    "\n",
    "    X = data_sp['coordinates/x'][:,:,zi]\n",
    "    Y = data_sp['coordinates/y'][:,:,zi]\n",
    "    U = data_sp['UI/1'][:,:,zi]\n",
    "\n",
    "    ax.plot_wireframe(X, Y, U, color = 'Black')\n",
    "    fig.savefig('ipotential.png', dpi=200)\n",
    "\n",
    "    \n",
    "interact(plot_timestep, xi = x_slider, yi = y_slider, zi = z_slider)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# --------------------\n",
    "# clear workspace\n",
    "# --------------------\n",
    "\n",
    "%reset"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}