{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notebook to solve idealized wave problem. Barotropic current interacting with a step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving this equation:\n",
    "\n",
    " $-\\sigma^2\\nabla^2\\psi + N^2\\frac{\\partial^2}{\\partial x^2}\\psi = 0 $ or \n",
    "\n",
    "$\\frac{\\partial^2}{\\partial z^2}\\psi - \\frac{N^2-\\sigma^2}{\\sigma^2}\\frac{\\partial^2}{\\partial x^2}\\psi = 0 $\n",
    "\n",
    "Boundary conditiions are no flow through sloid boundaries so $\\psi=0$.\n",
    "\n",
    "Assume $\\psi = \\psi_b + \\psi'$ barotopic and baroclinic parts.\n",
    "\n",
    "Set $\\psi_b = zc\\sin(kx)$\n",
    "\n",
    "where $c = \\sqrt{gh}$ and $k=\\sigma/c$. Then the differentatial equation simplifies to:\n",
    "\n",
    "$\\frac{\\partial^2}{\\partial z^2}\\psi' -\\frac{N^2 -\\sigma^2}{\\sigma^2}\\frac{\\partial^2}{\\partial x^2}\\psi' = \n",
    "\\frac{N^2 -\\sigma^2}{\\sigma^2}\\frac{\\partial^2}{\\partial x^2}\\psi_b$.\n",
    "\n",
    "We can subsitute these expressions for the derivaives of $\\psi_b$\n",
    "\n",
    "$\\frac{\\partial^2}{\\partial x^2}\\psi_b = - zck^2\\sin(kx)$ \n",
    "\n",
    "\n",
    "So, further simplication gives\n",
    "$\\frac{\\partial^2}{\\partial z^2}\\psi' -\\frac{N^2 -\\sigma^2}{\\sigma^2}\\frac{\\partial^2}{\\partial x^2}\\psi' = \n",
    "-k^2\\frac{N^2 -\\sigma^2}{\\sigma^2}\\psi_b$.\n",
    "\n",
    "Note that the domain depth $h$ is piecewise constant and so $\\psi_b, c, k$ are piecewise too.\n",
    "\n",
    "Setting up constants:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "L = 20e3 #units m\n",
    "H = 50 # units m\n",
    "h0 = 25 #units m\n",
    "N2 = 0.001 #units 1/s^2\n",
    "g = 10 #units m/s^2\n",
    "\n",
    "Nx = 100\n",
    "Nz = 50\n",
    "dx = L/Nx\n",
    "dz = H/Nz\n",
    "\n",
    "sigma  = 2*np.pi/(12*3600) #units 1/s\n",
    "\n",
    "x = np.arange(0, L,dx)\n",
    "z =  np.arange(0, H,dz)\n",
    "\n",
    "#meshgrid\n",
    "xx ,zz = np.meshgrid(x,z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up topography and currents..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x0 = L/2\n",
    "h = h0*(xx<x0) + H*(xx>=x0)\n",
    "\n",
    "c= np.sqrt(g*h)\n",
    "k = sigma/c\n",
    "\n",
    "psi_b = zz*c*np.sin(k*xx)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize $\\psi'$ = - $\\psi_b$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 20000.0, 50, 0]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAEACAYAAABmohcVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXuwXVWd5z8/khAQH5jCITzSBpyggtoC3dEZxQ4zxMae\nFqS6BWmHYdTpssRXzUwpD5026bHxUTO2TXVpO612IzWgmUEQBrUJNpfRKjGtREVClKCZMkECKigI\n5Obe/OaPs0/uvidrn7Wf9+y97/dTtSv7rLXXY+978rv7fn/r91vm7gghhJgsh0x6AkIIIWSMhRCi\nFcgYCyFEC5AxFkKIFiBjLIQQLUDGWAghWkAlY2xmZ5vZdjO7z8wurWtSQggxCSZp06zsOmMzWwL8\nEDgL2A38E3Chu99b3/SEEGJhmLRNq/JmvBbY4e473X0f8Hng3HqmJYQQC85EbVoVY3wc8NPU511J\nmRBCdJGJ2rQqxlhx1EKIPjFRm7a0QtvdwKrU51UMfpMcwMxksIUQuXF3K9u2jL0ZGS9q05qkijH+\nNrDGzFYDDwAXABeOXrShwgB94nbgzElPoiXoWcyx2J/FB1ILCMxK2+EDfLDAte8/uCiXTWuK0sbY\n3WfM7B3APwBLgM9oJYUQYpIsq9B20jatypsx7v4V4Cs1zUUIISpRyaAxWZtWde6Lgjoe0vNq6qcM\nVd4WmuD5wOEl23b9Czv6s3gJsGISE6lAm38GZb9XbaDNz7VXnDjpCbSI5016Ai1izaQn0DPa9uJR\nBBljIURv6LJB6/LchRBiHnoznjBN30TTP2DNPz9dv5eFMBZ9uIeydNmgdXnuQggxjzb/ooghYyyE\n6A0yxkII0QK0tG0BqDrRqr8xNf5kx08z6bl0/VlOevwm6YxBC9CJuU/6y6Px28Okn0WfnmUfafMv\nihj6bgghekOXDVqX5y6EEPPQm7EQQrSALhu01s29zISK/Dbsev9lxmpr/2XaFXm+Xe+/7Fh9+o63\ndZwmaJ0xFkKIsmhpmxBCtIAuvxlX2ZBUCCFaxdICRxHM7PVmdo+ZzZrZaany1Wb2pJltTY5PVJn7\nxKlTp2tC54z12ZTelrffOufX9udS188/T19F+p3EvTb2vYtcvDRnZ8uWFBi0JmJzn8dMoa7vBs4D\nPhWo2+HupxbqLUArjLEQiwX9h2uWvL8ogELG2N23Qz2bpmah74YQojdM4m0cOMHMtgK/At7v7t8o\n04mMsRCiN4x7M75jZnBkYWabgZWBqivc/eaMZg8Aq9z9kURLvtHMTnH3x3JPOkHGWAjRG5Ytz647\nazmclfr8wYfm17v7+qLjufs0MJ2c32Vm9zPY2vCuon11whhXddr0vX3b51dXH3X0M/FnneN/XEz3\njP0pHm0fm0PsT/1Y+0lalYUZ+4BwbGZHAY+4+6yZncjAEP+4TKda2iaE6A8NrW0zs/PM7KfAy4Fb\nzOwrSdXvAd9LNOP/BbzV3R8tO3UhhOgHDVk0d78BuCFQfj1wfR1jyBgLIfrDZFZT1MLEjHFZ7a4J\nzbCsDtiEPtlEn3141rH2peZSMbghpN1W1mshbFDK6LRNacNNa85V6PDrZYenLoQQI4xZTdF2ZIyF\nEP2hwxatw1MXQogROmzROjx1IYQYQQ68+qjLAVWn062Ig6tM0EDTc51U/7F2pX7WGReEHGe1Otia\ndqqVcXo1PdesOTXhTKyL1lm0/HR46kIIMUKHLVqHpy6EECN02KJ1eOpCCDGClrYtPFV13CI6apFr\nmxhr0vdaJNAic6xARV6dN9Ymq/9C2mlenbOIjtpE+6w+irRv4lmV0amboLMWLUeiIDP7rJntMbO7\nU2UrzGyzmf3IzG41syObnaYQQuRgSYGjZeTJ2vZ3wNkjZZcBm939JOBryWchhJgsDWVtWwiixtjd\nvw48MlJ8DnB1cn418Lqa5yWEEMXpsDEuO6Wj3X1Pcr4HOLqm+QghRHlaKD/kpfLvB3d3M/Os+ttT\n56sZpMEfpU6n00I5xdLlbXLQLagDMqnIdKoVCLoo5YBbSKfXkkh90w7EhXSwDcurOhsjz3JqF0xt\n2JBxUUkaeuM1s//KQBFw4BfAv3f3nyZ1lwNvBmaBd7n7raXGcM+0o+mJrAZudvcXJ5+3A+vc/UEz\nOwa43d1fEGjnG0bKml4BIGOcr5+sa2WMI+1ljGszxgB8fM7+mBnubgc3yIeZub+jwPV/Te7xzOwZ\nw01GzeydwG+7+38ws5OBa4HfBY4DbgNOcvf9Redfdtulm4CLk/OLgRtL9iOEEPXR0GqKkd2enw78\nPDk/F7jO3fe5+05gB7C2zNSjL/Vmdh2DfZ6OSvaA+jPgw8AmM3sLsBM4v8zgQghRKw065szsL4CL\ngCeZM7jHAnemLtvF4A25MNGpu/uFGVVnlRkwxkL9aV72T/cyMsOkZYhC46cKyyTfKZ1op0zQxELK\nDGX+NC8rLZSRVMqMVefPqkybJqgwjpltBlYGqq5w95vd/X3A+8zsMuDjwJsyuoprvwFauMBDCCFK\nMkZ+mPp/gyMLd1+fc5RrgS8n57uBVam645OywsgYCyH6wxiLtu55g2PIxm/k79bM1rj7fcnHc4Gt\nyflNwLVm9jEG8sQaYEuBGR9AxlgI0R+as2gfMrPnM1i+dj/wNgB332Zmm4BtwAxwiedZohZAxlgI\n0R8aytrm7n88pu5K4MqqY3TOGMcmHHOKxeqrOtAm4WyMOhAzOhg63rIypR2or8MpVpeDrEibqk7D\nqo6sMv3nadf159IknbNoc3R46kIIMUKHLVqHpy6EECMs5twUQgjRGjps0SY29TK7KGfV16kTV+2r\naX05NFatmnCdeuJCarJ1jRXTcYvo51U146bHqvNZF/m5N2l1ZIyFEKIFSKYQQogWcNikJ1AeGWMh\nRH/osEXr8NSFEGIEyRTVKBLIkdexVyRQoqxTLm+7Ilndok691IdYVrXKDrq8jpo82b+adkrV1X9V\nB2EsACbm6EqXV3WMVp1L1v3X5QxtglZYtHJ0eOpCCDFChy1ah6cuhBAjSKYQQogWoNUU9RHThIeU\n1XFDmnKR4ItQeUyfzuorrz48byeNIjttTFqHLaNJlg0aiems48bMGj/Uruyzit1fGf26yFzr1JTL\n+BIWytLozVgIIVpAhy1ah6cuhBAjdNiiHTLpCQghRG0sLXCUwMz+s5ntN7MVyefVZvakmW1Njk9U\nmboQQvSDBjVjM1sFrAdGtzXd4e6nVu2/tca4SKBHiFhWtiKBGDFnXtUMcVmBHENnXSi4Y167mKOl\nzuCBMo6g9HkRp1RVB14TTqumn2WevvryLJugWYv2MeC9wJea6FwyhRCiPywvcBTAzM4Fdrn79wPV\nJyQSxZSZvbLs1Fv7ZiyEEIUZY9Gmvjk4sjCzzcDKQNX7gMuBV6cvT/59AFjl7o+Y2WnAjWZ2irs/\nVmziMsZCiD4xxqKtO2NwDNn48fn17r4+1M7MXgScAHzPzACOB75jZmvd/SFgOml/l5ndD6wB7qpx\n6u0kr6Yb03FDfabblQ36iAaYZOjDB+qLJPopo+1NSkdtQueseq9V5xeaV9WfVdm+0jRxr7GfewpP\n2s2m6tPnM0vm1NEjwl2UpwGL5u4/AI4efjaznwCnu/svzewo4BF3nzWzExkY4h+XGadzxlgIIbLw\nrF+WNQ+TOn8V8Odmtg/YD7zV3R8t06mMsRCiN8wugEVz9xNT518EvlhHvzLGQojesBDGuCk6PPVy\nFFknnFc/rnMdcaFEP3Vpf3nWgy7UWGV10MX+rKpqzqnzvJrvbOqLO7tkboDZZLAZ0mXhH2zdmvHe\n5YcWuHq65tGrseiM8ZDYjZcJNOktRTygix09q4mS/qXQNfTVEEL0htnM5S7tR8ZYCNEbZmSMhRBi\n8mRp010gOvMkU9HngH/GYH3d/3D3q5IUcl8AngvsBM7PWl9X1SlWOdCixv5D5Zn1yYesnToOBHWk\n4+SbdtQ0HTTS9FiHReqbupcuOgADTjmYc8btXX5IquxgZ1z6T/6QM26Wg5126fMsB1762uOoly7L\nFHkSBe0D/qO7nwK8HHi7mb0QuAzY7O4nAV9LPgshxMSYZUnuo21E34zd/UHgweT8cTO7l8EvtHOA\n30suuxqYQgZZCDFB9lJkaVu7KCSwmNlq4FTgW8DR7r4nqdpDKnZbCCEmQa814yFm9nTgeuDd7v5Y\nkr0IAHd3M/PMxgtIE5pwTJ8+fNJBHVmac9M6aGjcpjXbsvda9VksD5Q18ayz+kpp5U0HZQzrp1Nv\nmUU04+G1WTpxkxJBG+WHvOQyxma2jIEhvsbdb0yK95jZSnd/0MyOAR4Ktb2dOWH6BOD5FSe8UCjo\nowDdfRkRE2Tr1K/5x6kNtfbZa2Nsg1fgzwDb3D2dAfQm4GLgI8m/Nwaac2aeQYQQi45T1z2TV63b\ncODzxo0bK/fZ93XGrwD+LfB9M9ualF0OfBjYZGZvIVna1sgMhRAiJ73WjN39G2QvgTur3ukIIUR5\nei1TtIEimdJi1BVUAnOOu2hQRywQ4rBIfZFAiphTrUj/hwXKygSlpOfVJmdl6LmXddDV5eADZpK+\nsnbHGDrmQk45mPtTfTr1sMoEZaTb7A301UYH3nSDS9vM7J3AJcAscIu7X5qUXw68OSl/l7vfWqb/\nThhjIYTIQ1OasZmdySC24iXuvs/MnpOUnwxcAJzMIP7iNjM7yd33Fx1DxlgI0Rsa1IzfBnzI3fcB\nuPvDSfm5wHVJ+U4z2wGsBe4sOkCecGghhOgEDYZDrwFeZWZ3mtmUmf1OUn4ssCt13S5Kptxo3Ztx\nSLONUUTzHdcmXR7TiSGlCQd0YkhpxbFEQE3ppGV0zDKBHDHNO6uv0Fhln0UTQRtlNOWMREaxQI30\nDhV5E/Wk9dEimu0wZDim6WaNP9SiiyQS2jtvrs2ZnSp6tJltBlYGqt7H4Kf5bHd/uZn9LrAJODFw\nLczfsDQ3rTPGQghRlnGa8Q+mfsE9U7/MrHf39Vl1ZvY2ko1H3f2fzGy/mR0F7AZWpS49PikrjIyx\nEKI3TM/702k+J607lpPWHXvg86aNO4p0fSPwr4A7zOwk4FB3/7mZ3QRca2YfYyBPrAG2lJi6jLEQ\noj80uGzus8BnzexuBjuZ/jsAd99mZpuAbcAMcIm7S6YQQixumlralqyWuCij7krgyqpjtMIYF5lE\naHeNcdeN9p/XmZfVvlBWtrxOoyJBH0WCOvIGbRRxihVxanVxrjU66EI7aWQFahQJqtgbcKBVDfoI\nXZsVFDJ0HFbN+tYEvQ6HFkKIrqBwaCGEaAEyxkII0QJkjBsmNsksfXdcfZbmHAwaCejEkJEIKKaT\nFgl0iOmYZepDmmpWfV6dNlbfprnEAmhq1IRDgQ6z83Ta/DtplOkrpPNm1RcJColpxnk15SbYO2Zp\nW9vphDEWQog86M1YCCFagIyxEEK0gL5vuySEEJ1A64xLUCRoI0TWxGNZ32JZ3Q44+FIXLM04z52V\nrYjTqMxOIFn9N+EAjNWXdJDlfhZZ/ceCRiL34oHdNUJZ1ao66GJOtfQ1saxsWUEdsaCQeP14B15d\nDr4mkEwhhBAtQMZYCCFawN4G37qbRsZYCNEbpBnXSJEAj3FtigSCzDtPLs7a8Tm403PVQIes5DZ5\nk99k6ayxQIqFGr+OvmL6ec578VT9dKp+uPvy9PKwDhvXaQ/WUYtoxlU13zL6dda9LOS91o1kCiGE\naAEyxkII0QK0zlgIIVpAU5qxmX0eeH7y8UjgUXc/1cxWA/cC25O6b7r7JWXGkDEWQvSGpmQKd3/D\n8NzM/hvwaKp6h7ufWnWMThjjkNMuKyikTNa2ef63YSa2QHY2IL4TR6i+iNOqTFBFnUEfVfsv62Cr\nGhQScFbOpOY9DOAYBm9AOIAj5nQK7ciRpz4UCFEkaCPkoKvTAVjEwVe1fZNSQpMBJQBmZsD5wJl1\n931I/BIhhOgGMyzJfZTkDGCPu9+fKjvBzLaa2ZSZvbJsx514MxZCiDyM04wfntrGw1P3Ztab2WZg\nZaDqCne/OTm/ELg2VfcAsMrdHzGz04AbzewUd3+s6NxljIUQvWGcZrxi3YtZse7FBz7fu/GGefXu\nvn5c32a2FDgPOC3VZhqYTs7vMrP7gTXAXUXnvuDGOKbpjruujvZB/Tlw8bygjyLJaaoGZVTVbGPt\nQ5psnf2XvbbEs/SgJnzw7hsA00vyabYhbXZw7cK0H/QxXpNtQlNuun3WfddNw+uMzwLudfcHhgVm\ndhTwiLvPmtmJDAzxj8t0rjdjIURvaHid8QXAdSNlrwL+3Mz2AfuBt7r7owe1zIGMsRCiNzSZm8Ld\n3xQo+yLwxTr6lzEWQvSGppe2NYmMsRCiN/Q2HNrMDgPuYOB2ORT4krtfbmYrgC8AzwV2AueX1Uny\nUNVBF9rVIxToAXPBHvOcelnnoaxtRYI+8jrA6nQg1hnUUWSnj2F5xaxr6UCOkLNu6KiDsNMo5lQK\nOarS5SFHWrq+rAOvzqCLuhx0TT+rJuhyCs2xQR/u/hRwpru/FHgJcGayqPkyYLO7nwR8LfkshBAT\nZZYluY+2EY3Ac/cnktNDGby3PAKcA1ydlF8NvK6R2QkhRAG6bIyj7/RmdgiDBczPAz7p7veY2dHu\nvie5ZA9wdINzFEKIXLTRyOYlaozdfT/wUjN7FvAPZnbmSL2bmdc1oTK7Rock3aydPA6UZezksTSk\n88b046YS9eRNxFNEE25CU86z00fsXksk+gnpw0V02CZ01JhmnKWzxuZaRlOe9L3kCXCpm73zvlTd\nIrfa7e6/MrNbgNOBPWa20t0fNLNjgIey2t3OnBZyAnBKldkKIXrD9NQ32TD12Vr77O2bcRLqN+Pu\nj5rZ4cB6YCNwE3Ax8JHk3xuz+jgzNogQYlFy6Lp/wYZ1Fxz4vHHjxsp99tYYA8cAVye68SHANe7+\nNTPbCmwys7eQLG1rdppCCBGnt+uM3f1uUhmKUuW/ZJA0ozT7gMPH1M9UrM/qf3jD+2bg8ESnnLfO\neHg+S3x351By+dB5bMdjyL+muKp+W2f7UF8l9evhrs3pHZtjieDr1Ifj9ePXFMeSBoXrw/OPJbqv\nus64yFxi9UF9e3+qfiapfyr1s3gyFSVXs+u/y+uMJzbzcYZ0QeoPi1ywkH6A2FhV59p0+wX8FqWN\nSojYm1HsP2u8PtZ/tfqFpOpcY896aIgXkjY936J099eIEEKMIGMshBAtYO+0EgUJIcTEmZ3prknr\nxMxjk4wGeKTPk86WZf01kzf5T9a1IadVkUCIqk67mNMsEmgRvLaqUzB1bWh3DphL+jO9vJzTLZa8\nZtguFshQNhAjr9Msz+4XcQfewQ7MWFBHkUQ+0aRGiYNub8opN5vSh4cOuv17Uz/stH78lM2d1+3A\na0inNrPfBv4GOILBCrI3Dve5M7PLgTczcPu/y91vLTOGdocWQvSG2ZkluY+CfBp4r7u/BLgBeA+A\nmZ3MYAeQk4GzgU8kS4ELI2MshOgNM/uW5D4Kssbdv56c3wb8UXJ+LnCdu+9z953ADmBtmbnLGAsh\nesP+2aW5j4LcY2bnJuevB1Yl58cCu1LX7QKOKzP3VmjGWZLpuGtjbZZl1A8DPIKBHumLZ8i/e3Ms\nkCMW9BFL5BNLrlMkOX2Z5D2x+eW4dpj0Z75OfHDSn7LBB2WS4xTSSWvSlJsaK66fH1wfDQRJrUwI\nBW1ENeGnSJWlzmdojgqasZltBlYGqq5goAlfZWb/hUE6iOkxXZVKnNYKY9xKYoEOYg49q9x0eY+2\nTvDUGJP2rSnYMpVZ7e7rI73/PoCZnQT8m6RsN3NvyQDHJ2WFkTEWQvSHcW/dp68bHEP+On9iIjN7\njrs/nDjn3g98Mqm6CbjWzD7GQJ5YA2wpNOcEGWMhRH9oTgK50Mzenpxf7+5/D+Du28xsE7AtGf0S\nd5dMIYRY5DRkjN39KuCqjLorgSurjtFaYxzb8SOL4E4fsbuM7egcy8oWqy8S9FHAKRbN9BYLCqka\nVBKoLxvUEQqKKJIpLeb0KhOUEQ4aCQda5B0/5mCss6/ovaQcdOnItemnkqCSlNOOp1Jad8hBty91\nvpeD69PnT9Ic++KXtJXWGmMhhCjM7KQnUB4ZYyFEf2hy2VzDyBgLIfrDU/FL2konjHEoEVBWcqDo\n7tKJtBYM9Eif50kUlDeoI9ZXTBOOJQIqG7SRNygkEuiR1omzduqIBXWEkt+USRSUpSnHdgdpQlMu\nc3919HVgrqkdN0JJfdI7bswL4Biu1c3SeWOa8EykvkmDqTdjIYRoATLGQgjRAmSMhRCiBWhpmxBC\ntAAtbWsfoUxt886L7N4Rc8DFgjqKZE0rEoBSxgEX28mjgANwmIltGNAB5YM6mgiqiO30UcbBV2f/\nsQCW2E4god03YM5ZN50K1AgGcKR33Ag52GIOuLL1jWZta7DvhumtMRZCLEK0tE0IIVqA3oyFEKIF\nyBiPp0zQRtncPAf6T+/kkdJcl8U6qJooKBb0ETpfSM23TFBIQCeG/Ml/BudlNN/iQR1P8rRS7Rdq\nfkV2Asmcy/TBiXyCARzpROuhAI69hOufGrku69qyQSEK+giiN2MhRH/Q0jYhhGgBHV7apt2hhRD9\n4akCRwHM7PVmdo+ZzZrZaSN1l5vZfWa23cxenSqfSsq2JsdR48bQm7EQoj80pxnfDZwHfCpdaGYn\nAxcAJzPYA+82M1uTbL3kwJ+4+115Blh8xjgU1FEm01pWX7H6UF9NO+BKOgj9iMG/VTOxDcrrcaA9\nweGV2rd9fvOuTe3EEXLWBTOtQX4H3OOBsvR5kaCO0FhZ43dwpw933w5gZqNV5wLXufs+YKeZ7QBe\nBtyZ1B/UIAvJFEKI/jBb4KiHY4Fdqc+7krIhVycSxftjHS2+N2MhRH8ZJ1PsmoLdU5nVZrYZWBmo\nusLdby4xmze6+wNm9nTgejO7yN2vybpYxlgI0R/GGeOV6wbHkC0b51W7+/oSI+4GVqU+H5+U4e4P\nJP8+bmbXAmuB/hjjWNDIAck3lBwoiyI6cJEAkbyabhHNuYmgkCPmikK7doR0YojrsPMDMNod1PFE\n0m/5nTwOvr8nUnPNG8gBc/pwdCeOOjXf0HnT/TfBwqwzTuvANwHXmtnHGDjw1gBbzGwJ8Gx3/7mZ\nLQNeC9w6rtNcxjjp+NvALnd/rZmtAL4APBfYCZzv7o8WvCEhhKiXvfFLymBm5wFXAUcBt5jZVnd/\njbtvM7NNwDYG7+WXuLub2WHAVxNDvATYDPztuDHyvhm/OxnsGcnny4DN7v5RM7s0+XxZwfsTQoh6\naWhpm7vfANyQUXclcOVI2W+A3ykyRnQ1hZkdD/wB8GnmXs/PAa5Ozq8GXldkUCGEaIR9BY6WkefN\n+C+B9wDPTJUd7e57kvM9wNF1T0wIIQrT4XDoscbYzP4QeMjdt5rZutA1iT7iWX3czpyD7UTgRTkn\nluWgi7ZL7iiYqQ3CQR0E6mOBHunyLAdZ3gCSIkEdsaxvRYJCEsddyGkHc467kNMO4MkkwKFIIEPI\nqZW+tt72BzsTYwEqWTtx1NY+EsgBGVnXQs60Ig60UPsiQSFlgjp+E5nfb6bYsGGKWulx1rZ/CZxj\nZn/A4L/xM83sGmCPma109wfN7BjgoawOzoRUTJIQQiQcsY4NG9Yd+Lhx48bsa/PSYWM8VjN29yvc\nfZW7nwC8AfhHd7+IwXKOi5PLLgZubHaaQgiRg55rxmmGcsSHgU1m9haSpW11TkoIIUrR0NK2hSC3\nMXb3O4A7kvNfAmc1MaG8+nBod49Miuy4nHd3jnQfRRIJFQn6yBvUUTYRUHIeSwRURKcte23enTCK\nJOIJXZulSQ+vje3eXGj+qR2bn3g80YxjgRxQLugipimn9d3fBMpCmm+6rIhmPbx2EomCOixTdC4C\nTwghMmmh/JAXGWMhRH/o69I2IYToFJIphBCiBcgYV6NIUrXotUvn/5vZWVYgRmigWFa2mLOvbNBG\n3qCOkNMvXR/YvQPiWdlimcxCTq+yDrK5oIpqWduyHYgHO9hiWdWeTDkAY0EpT+5PAmBSgRzTT6Xm\nlzjwgo66wWBz/CZQH3KQZdWHgi5iQRuhAI0sB2HIwVdk/t3P2tYIrTDGQghRC4thaZsQQrQeyRRC\nCNECJFPUP3CsPrS7x7z6mNBcZHeNLP04tvtzKKgjVh8LCokFdUTqY4mAYjpqLJAiFAiRvjam6YZ2\nB0lfGwsqCem8We2LaNZB/TyV9OfJRBNOJ/8hpRnzuA0HmqNIUEUs0U9IHy7S/2OBeVXVlDPfUtMd\n15y5pqGlbWb2emAD8AJgrbt/JylfD3wIOBSYBt7j7rcndacDf8/gf+CX3f3d48bQ7tBCiP4wU+Ao\nxt3AecD/ZS4tBMDDwB+6+0sY5OlJ73H3SeAt7r4GWGNmZ48bQDKFEKI/NLfTx3YAMxst/27q4zbg\n8GSrpaOAZ7j7lqTucww24fhq1hgyxkKI/jBZzfiPgO+4+z4zOw7YlarbzWDD0kxkjIUQ/aHCm7GZ\nbQZWBqqucPebI21PYZDNcn3Z8VtnjIs47g6UpRotyxvAEcvaFgv0SJfH6rOyqpUJCknPK7STR6B+\nJlUfCvCI7dQRy5QWc4ql+4pleIs56LLmkj+oJH9Qx7y5JBnYhtnXICMD2+MZQR2xQI0yDrZQGcR3\n2hi2Kxu0EcrKFjSC6Qt+nVG+OtSwIaaSI4y7lzKkyT6hXwQucvefJMW7geNTlx2flGXSOmMshBDN\nsC45hpTeWeSAcGxmRwK3AJe6+zeH5e7+MzP7tZm9DNgCXARcNa5TraYQQogIZnaemf0UeDlwi5l9\nJal6B/A84ANmtjU5jkrqLgE+DdwH7HD3TOcd6M1YCNErmvHgufsNwA2B8g8CH8xo8x3gxXnH6IQx\nDurEGdcuDWm+8y4IlIX04VigR/qaLM021lcoqKPITh5LxtcPteK9y+f+AArt9JwVKBFLBPREIFAi\nK2gjtJN0SP/N0oRjiYJCmnDl/mNBHY+nAzyGZRxclj7P0nlD7bL6iu2kEdKnQ0EZZTXloD4c0oSz\ndOImNeMruGyRAAAHgUlEQVTuxkN3whgLIUQ+uhsPLWMshOgRTW6w1ywyxkKIHqE340oDFtKEQ2VF\nEgENz2cZvzvzXuCIMfXp89juzkXWIceS14f04VSfnqofasXTy9Pa6PgdlaczEuXkT7Qzfh1vkXXC\nRdYhx/Xr8YmIQsnhn3hsbq651xGny2Oab9X2sUQ+MX26dKKfX49cmC5Llz8WKBu9tm6kGXeP5ZH6\nIyL1i4i00QuRNnCLntguFo9H6kVF9GYshBAtQG/GQgjRAvRmLIQQLUCrKRacrIkvC1XEkgeVSSRU\n5NpYop8sB92SketGz5Nr00679E4eQ8fd/N0rQjsm50/eU8ZpN+jj8HltyvaVdgCGnHFFAliemE45\nExPH3ZOpRECFgjrKONWydm8OJeoJOeNicyky19yBHOnydNljkfqsvupGMoUQQrQAyRRCCNEC9GYs\nhBAtQG/GtREK9ogmnK+qCYfapYM+YppwluYbSz4fSvQTujaiOc+m2qSTx89pwvkTrocSuscSAcV0\n2vljhfsK7b4c0qpDbdLtnojVP5EqS+3ePD3Uioskh3+8QH1Icw4FWkBcMy6jKafrcyeCz9J5H4vU\nxzRlacYhWmeMW4OCPoToIHozFkKIFtDdpW3a6UMI0SP2FTjyY2avN7N7zGzWzE5Lla9N7fDxfTO7\nIFU3ZWbbAzuABNGbsRCiRzSmGd8NnAd8KlB+urvvN7OVwA/M7H+7+yzgwJ+4+115BshljM1sJwNV\nfhbY5+5rzWwF8AXgucBO4Hx3fzRPf0II0QyNbbu0HcDMRsvTusjhwK8SQzxkfoMx5H0zdmCdu/8y\nVXYZsNndP2pmlyafLxtt+AH3vHPpNVNTU6xbt66RvtM/xGdmnLeJJp9FlKdlnK9Y6IkMmOizmMfh\ngfOjJzGRiiz8agozWwv8HXACcOFI9dVmtg+4PtkvL5MimvGohT8HuHo4IPC6An0tOqampiY9hdag\nZzGHnkXdjNOI7wX+T+qYj5ltNrO7A8drx43o7lvc/RTgNOCvzOxZSdUb3f1FwBnAGWZ20bh+irwZ\n32Zms8Cn3P1vgaPdfU9Sv4du/hoVQvSKcW/Gq5m/Aeqt82rdfX2Vkd19u5ndD/xz4Dvu/kBS/riZ\nXQusBa7Jap/XGL/C3X9mZs8BNpvZ9pFJuJlJjxBCTJgFWdp2QCUws9XALnefMbPnAmuA+8xsCfBs\nd/+5mS0DXsuo9R/t1Atqumb2AQaxPX/KQEd+0MyOAW539xeMXCsDLYTIjbvndniNUsbe5B3PzM4D\nrgKOAn4FbHX31yTSw6XMaSF/5u5fNbMjgDsYBBUvATYD/8nHGNyoMTazpwFL3P2xZIBbgY3AWcAv\n3P0jZnYZcKS7H+TAE0IIESePMT4BuCH5uBT4n+7+oWRp2ybgt9DSNiGEqERhmUIIIUT9NBIObWZn\nJ2GA9yVrkHuJme1MQiC3mtmWpGxFskTmR2Z2q5kdmbr+8uSZbDezV6fKT0+W0NxnZn81iXspipl9\n1sz2mNndqbLa7t3MlpvZF5LyOxPnSCvJeBYbzGxXKhT2Nam6Xj4LM1tlZrcnYcM/MLN3JeWL8ntR\nGHev9WAgVu9gsIZkGfBd4IV1j9OGA/gJsGKk7KPAe5PzS4EPJ+cnJ89iWfJsdjD3l8kWYG1y/mXg\n7EnfW457PwM4Fbi7iXsHLgE+kZxfAHx+0vdc8Fl8gIHDZvTa3j4LYCXw0uT86cAPgRcu1u9F0aOJ\nN+O1wA533+nu+4DPA+c2ME5byBsMcy5wnbvvc/edDL54L0tWojzD3bck132ODgTQuPvXgUdGiuu8\n93Rf1wP/uvabqImMZwHhUNjePgt3f9Ddv5ucP84gyuI4Fun3oihNGOPjgJ+mPu9KyvrIMBjm22b2\np0lZVjDMsQyexZDhcxkt3013n1ed937ge+TuM8CvEqdxl3inmX3PzD6T+tN8UTyLZP3tqcC30Pci\nF00Y48XkEXyFu58KvAZ4u5mdka70wd9Si+l5HGAx33vCJxnkKngp8DPgv092OguHmT2dwVvru909\nvcWHvhdjaMIY7wZWpT6vYv5vud7g7j9L/n2YwfK/tcAeG6TSI/lz66Hk8tHncjyD57I7OU+X7252\n5o1Rx73vSrX5raSvpcCzfH6iqlbj7g95AvBpBt8N6PmzSKLNrgeucfcbk2J9L3LQhDH+NrDGzFab\n2aEMRPabGhhnopjZ08zsGcn5EcCrGeQ2vQm4OLnsYmD4hbwJeIOZHZqs3V4DbHH3B4Ffm9nLzMyA\ni1JtukYd9/6lQF9/DHxtIW6gLhKjM+Q8Bt8N6PGzSOb9GWCbu388VaXvRR6a8Aoy+LP9hwwE+csn\n7aVs6B5PYOAJ/i7wg+F9MkjGeBvwIwbRikem2lyRPJPtwO+nyk9n8J91B3DVpO8t5/1fBzwATDPQ\n8N5U570z2IZ1E3AfcCewetL3XOBZvJmB0+n7wPcYGJ+j+/4sgFcC+5P/E1uT4+zF+r0oeijoQwgh\nWoD2wBNCiBYgYyyEEC1AxlgIIVqAjLEQQrQAGWMhhGgBMsZCCNECZIyFEKIFyBgLIUQL+P/9Cfbr\n4s7frgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3d7353ff28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "psi_c = -psi_b\n",
    "plt.pcolormesh(xx,zz,psi_c)\n",
    "plt.colorbar()\n",
    "plt.axis([0,L, H,0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerical integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def integrate(psi):\n",
    "    \n",
    "    for i in np.arange(2, Nx-1):\n",
    "        for j in np.arange(1,Nz-1):\n",
    "            if zz[j,i+1] > h[j,i+1]:\n",
    "                psi[j,i+1] = -psi_b[j,i+1]\n",
    "            else: \n",
    "                psi[j, i+1] = 2*psi[j,i] - psi[j,i-1] + \\\n",
    "                          dx**2/dz**2*sigma**2/(N2-sigma**2)*(psi[j-1,i] - 2*psi[j,i] + psi[j+1, i]) + \\\n",
    "                          k[j,i]**2*dx**2*psi_b[j,i]\n",
    "\n",
    "    \n",
    "    return psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 20000.0, 50, 0]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAEACAYAAABmohcVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+wXVWV5z/LJATEH0jhEH6kDThBBbUFuqPVih16iI09\nLUh1A9IOw6jTZYm/amZK+aHTJj0O/qgZx6a6dJxWu5Ea0MwgCIPaBJtHa42YVqIiIUrQTJkgARUU\nhOTlvaz5456bd97NPneffX68e+5930/VqZy79z5773PezXrnfddea5u7I4QQYrQ8bdQTEEIIIWMs\nhBCdQMZYCCE6gIyxEEJ0ABljIYToADLGQgjRAWoZYzM728y2mdn9ZnZZU5MSQohRMEqbZlXXGZvZ\nEuCHwFnALuCfgIvc/b7mpieEEAvDqG1anTfjNcB2d9/h7vuAzwPnNjMtIYRYcEZq0+oY4+OAn+Y+\n78zKhBBiHBmpTatjjBVHLYSYJEZq05bWuHYXsDL3eSW93yQHMDMZbCFEadzdql5bxd4MjBe1aW1S\nxxh/G1htZquAB4ELgYsGG62vMcAkcQdw5qgn0RH0LOZY7M/iA7kFBGaV7fABPpjQ9v0HF5WyaW1R\n2Ri7+4yZvQP4e2AJ8BmtpBBCjJJlNa4dtU2r82aMu38F+EpDcxFCiFrUMmiM1qbVnfuioImH9PyG\n+qlCnbeFNngBcFjFa8f9Czv4s3gpcOQoJlKDLv8Mqn6vukCXn+tEceKoJ9Ahnj/qCXSI1aOewITR\ntRePFGSMhRATwzgbtHGeuxBCzENvxiOm7Zto+wes+Zdn3O9lIYzFJNxDVcbZoI3z3IUQYh5d/kUR\nQ8ZYCDExyBgLIUQH0NK2BaDuROv+xtT4ox0/z6jnMu7PctTjt8nYGLQAYzH3UX95NH53GPWzmKRn\nOYl0+RdFDH03hBATwzgbtHGeuxBCzENvxkII0QHG2aB1bu5VJpTy23Dc+68yVlf7r3JdyvMd9/6r\njjVJ3/GujtMGnTPGQghRFS1tE0KIDjDOb8Z1NiQVQohOsTThSMHMzjeze81s1sxOy5WvMrOnzGxL\ndnyiztxHTpM6XRs6Z6zPtvS2sv02Ob+uP5emfv5l+krpdxT32tr3LtJ4acnOli1JGLQhYnOfx0xS\n1/cA5wGfCtRtd/dTk3oL0AljLMRiQf/h2qXsLwogyRi7+zZoZtPUIvTdEEJMDKN4GwdOMLMtwK+A\n97v7N6p0ImMshJgYhr0Z3znTO4ows03AikDVle5+S8FlDwIr3f3RTEu+ycxOcffHS086Q8ZYCDEx\nLFteXHfWcjgr9/mDD8+vd/d1qeO5+zQwnZ3fbWYP0Nva8O7UvsbCGNd12kz69V2fX5k+FmouI3/W\nJR5ETPeM/SkevT42h9if+mV/mEXt2rQ6C2PRDgjHZnYU8Ki7z5rZifQM8Y+rdKqlbUKIyaGltW1m\ndp6Z/RR4BXCrmX0lq/p94HuZZvy/gLe6+2NVpy6EEJNBSxbN3W8EbgyU3wDc0MQYMsZCiMlhNKsp\nGmFkxriqdteGZlhVB2xDn2yjz0l41rHrK82lZnBDSLutrddC2KDErqtSX1UbrqI5L5SlGePXyzGe\nuhBCDDBkNUXXkTEWQkwOY2zRxnjqQggxwBhbtDGeuhBCDCAHXnM05YBq0umW4uCqEjTQ9lxH1X/s\nuko/64IGIcdZow62ppxqKf2Puq+ifpq876bpnEUrzxhPXQghBhhjizbGUxdCiAHG2KKN8dSFEGIA\nLW1beOrquCk6akrbNsYa9b2mBFoUjhWoaDKoIqj1pui8beiobVxf1EfK9W08q7rXN8XYWrQSiYLM\n7LNmttvM7smVHWlmm8zsR2Z2m5kd0e40hRCiBEsSjo5RJmvb3wJnD5RdDmxy95OAr2WfhRBitLSU\ntW0hiBpjd/868OhA8TnANdn5NcDrG56XEEKkM8bGuOqUjnb33dn5buDohuYjhBDV6aD8UJbavx/c\n3c3Mi+rvyJ2vopcGf5AmnU4L5RTLl3fJQbegDsisotCp1rYDbiGdXksi9QvlQBzVXJZE6ivMf2on\nTK1fX9CoIi298ZrZf6KnCDjwC+DfuPtPs7orgDcDs8C73P22SmO4F9rR/ERWAbe4+0uyz9uAte7+\nkJkdA9zh7i8MXOfrB8raXgEgY1yun6K2MsaR62WMGzPGAHx8zv6YGe5uB19QDjNzf0dC+7+m9Hhm\n9sz+JqNm9k7gt93935rZycB1wO8CxwG3Aye5+/7U+Vfddulm4JLs/BLgpor9CCFEc7S0mmJgt+dn\nAD/Pzs8Frnf3fe6+A9gOrKky9ehLvZldT2+fp6OyPaD+AvgwsNHM3gLsAC6oMrgQQjRKi445M/vP\nwMXAU8wZ3GOBu3LNdtJ7Q04mOnV3v6ig6qwqA8ZYqD/Nq/7pXkVmGLUMkTR+rrBK8p3KiXaqBE0s\n5J/2Vf40ryoNVJFUqozV5M+qyjVtUGMcM9sErAhUXenut7j7+4D3mdnlwMeBNxV0Fdd+A3RwgYcQ\nQlRkiPww9f96RxHuvq7kKNcBX87OdwErc3XHZ2XJyBgLISaHIRZt7fN7R58N3yjfrZmtdvf7s4/n\nAluy85uB68zsY/TkidXA5oQZH0DGWAgxObRn0T5kZi+gt3ztAeBtAO6+1cw2AluBGeBSL7NELYCM\nsRBicmgpa5u7/+mQuquAq+qOMXbGODbhmFMsVl/XgTYKZ2PUgVjQQd/xFnLKzatvwinWlNMnxalW\n12nY9vyqXjeJz6Upxs6izTHGUxdCiAHG2KKN8dSFEGKAxZybQgghOsMYW7SRTb3KLspF9U3qxHX7\naltfDo3VqCbcpJ64kNpjU2PFduKoq5/H5rSQYzX5rOtq3k0hYyyEEB1AMoUQQnSAQ0c9gerIGAsh\nJocxtmhjPHUhhBhAMkU9UgI5yjr2UgIlqjrlyl6XktUt6tTLfYhlVavtoCvrqEkJZEhxSrXhLGzb\nQRhzqpUJfqjrgKvirExxujXlDG2DTli0aozx1IUQYoAxtmhjPHUhhBhAMoUQQnQAraZojpgm3Keq\njhvSlFOCL0LlMX26qK+y+vC8nTRSdtpoSoetqvNW0SSrBo3U1VmrBH00eX9V9OuUuTapb1fxJZQJ\ndmkCvRkLIUQHGGOLNsZTF0KIAcbYoj1t1BMQQojGWJpwVMDM/oOZ7TezI7PPq8zsKTPbkh2fqDN1\nIYSYDFrUjM1sJbAOGNzWdLu7n1q3/84a45RAjxCxrGwpgRgxZ17dDHFFgRx9Z10ouGPedTFHS1uB\nElV2+khxStV14LXhtKobiLEQzs6mnJEL6QBsinYt2seA9wJfaqNzyRRCiMlhecKRgJmdC+x09+8H\nqk/IJIopM3tV1al39s1YCCGSGWLRpr7ZO4ows03AikDV+4ArgNfkm2f/PgisdPdHzew04CYzO8Xd\nH0+buIyxEGKSGGLR1p7RO/ps+Pj8endfF7rOzF4MnAB8z8wAjge+Y2Zr3P1hYDq7/m4zewBYDdzd\n4NS7SVlNN6bjhvrMX1c16CMaYFKgDx+oT0n0U0XbW0gdtW3NuO69thGUUfdnVaav2PyauteEn7sH\n5je7NHw+s2ROHT384Mvq0YJFc/cfAEf3P5vZT4DT3f2XZnYU8Ki7z5rZifQM8Y+rjDN2xlgIIYoI\n/VJoY5jc+auBvzSzfcB+4K3u/liVTmWMhRATw+wCWDR3PzF3/kXgi030K2MshJgYFsIYt8UYT70a\nKeuEy+rHTa4jTkr005T2V0YnbXus2PgLpQmPw7Oqok8n6Nf9P/Vjmu9s7os7u2RugNlssBnyZUsP\nqofmNeO9yw9JaD3d8Oj1WHTGuE/sxqsEmkwsKR7QxY6e1UjJ/1IYN/TVEEJMDLPBZSfjgYyxEGJi\nmJExFkKI0TM7xiYtOvMsU9HngH9Gb33d/3D3q7MUcl8AngfsAC4oWl9X1ylWO9Ciwf5D5YX12Yei\nnToOBHXk4+TbdtS0HTQyqrGWR+oXowOwICij74zbu/xpubKDnXH5P/lDzrhZDnba5c/LOPB+i2YZ\nZ5miTKKgfcC/c/dTgFcAbzezFwGXA5vc/STga9lnIYQYGbMsKX10jeibsbs/BDyUnT9hZvcBxwHn\nAL+fNbsGmEIGWQgxQvaSsrStWyQJLGa2CjgV+BZwtLvvzqp2k4vdFkKIUTDRmnEfM3sGcAPwbnd/\nPMteBIC7u5l54cULSBuacEyfPmzUQR1FmnPbOmho3LY126r3Oi7PukRfns2h7aCM6dxbZopm3G9b\npBO3KRF0UX4oSyljbGbL6Bnia939pqx4t5mtcPeHzOwY4OHQtXcwJ0yfALyg5oQXCgV9JDC+LyNi\nhGyZ+jX/MLW+0T4n2hhb7xX4M8BWd89nAL0ZuAT4SPbvTYHLObPMIEKIRcepa5/Fq9euP/B5w4YN\ntfuc9HXGrwT+FfB9M9uSlV0BfBjYaGZvIVva1soMhRCiJBOtGbv7NyheAndWs9MRQojqTLRM0QVS\nMqXFaCqoBOYcd9GgjlggxKGR+pRAiphTLaX/QwNlVYJS8vPqkgMt9NxHFXQSCNCYzs0vvztG3zEX\ncsrB3J/q07mHVSUoI3/N3kBfVR14bS4/m26xbzN7J3ApMAvc6u6XZeVXAG/Oyt/l7rdV6X8sjLEQ\nQpShLc3YzM6kF1vxUnffZ2bPzcpPBi4ETqYXf3G7mZ3k7vtTx5AxFkJMDC1qxm8DPuTu+wDc/ZGs\n/Fzg+qx8h5ltB9YAd6UOUCYcWgghxoIWw6FXA682s7vMbMrMficrPxbYmWu3k94bcjKdezMOabYx\nUjTfYdfky2M6MeQ04YBODDmtOJYIqC2dtIqOWSWQI6Z5F/UVGqvqsygboJLyLKpoyoeG68vungEw\nvbw32Viinrw+WkWzjQVlpGjGsaCQvfPm2p7ZqePAM7NNwIpA1fvo/TSf4+6vMLPfBTYCJwbawvwN\nS0vTOWMshBBVGaYZ/2DqF9w79cvCendfV1RnZm8j23jU3f/JzPab2VHALmBlrunxWVkyMsZCiIlh\net6fTvM5ae2xnLT22AOfN27YntL1TcAfAHea2UnAIe7+czO7GbjOzD5GT55YDWyuMHUZYyHE5NDi\nOuPPAp81s3vo7WT6rwHcfauZbQS2AjPApe4umUIIsbhpa2lbtlri4oK6q4Cr6o7RCWOcMonQ7hrD\n2g32X9aZV3R9Ula2sk6jlKCPlKCOskEbKU6xFKfWOM61QQddaCeNokCNlKCKfttY0EVK0EeobVFQ\nSN9xWDfrWxtMdDi0EEKMCwqHFkKIDiBjLIQQHUDGuGVikyzSd4fVF2nOwaCRgE4MBYmAYjppSqBD\nTMesUh/SVIvqy+q0sfoyc+mXx/TzunNJCaDJzaWKJhwKdCjWgYe3jdWH9NmQzlu1PjRWVU25zWQ+\ne4csbes6Y2GMhRCiDHozFkKIDiBjLIQQHWDSt10SQoixQOuMK5AStBGiaOKxrG+xrG4HHHy5BksL\nzktnZYs5jWJOq6pOrTYcgLH6lHut8iyqPuvIvXhWX5RVbe/yLNNZJGgj5qCLOc3ybWJZ2YqCOmL1\nVeaaUh9y4BXNpWkkUwghRAeQMRZCiA7Q5v56bSNjLISYGKQZN0hKgMewa1ICQeadZ42LdnwO7vRc\nN+iiKLlN2eQ3RTprLKhjocZvoq+Yfh5LNJTVe64+tPtyf5cNSNVh03XWIk04rukeXF8lKKRJTTkl\nqESacZjOGWMhhKiKjLEQQnQArTMWQogO0JZmbGafB16QfTwCeMzdTzWzVcB9wLas7pvufmmVMWSM\nhRATQ1syhbu/oX9uZv8FeCxXvd3dT607xlgY45DTrigopErWtnn+t34mtkB2NiC+E0eoPsVpVSWo\nosmgj7r9pzjYYjtx1HRWzuTm3Q/g6AdvQDjDWszpFAtkCO3Yke+rTCBEWQddkw7Aug6+Kte3QZsZ\n4QDMzIALgDOb7vtp8SZCCDEezLCk9FGRM4Dd7v5AruwEM9tiZlNm9qqqHY/Fm7EQQpRhmGb8yNRW\nHpm6r7DezDYBKwJVV7r7Ldn5RcB1uboHgZXu/qiZnQbcZGanuPvjqXOXMRZCTAzDNOMj176EI9e+\n5MDn+zbcOK/e3dcN69vMlgLnAaflrpkGprPzu83sAWA1cHfq3BfcGMc03WHtmrg+qD8HGs8L+ijS\nTMvqoClBGXU129j1IX26yf6rtq3yLAP6cCi5D8D0knKabUib7bVdmOt7faRrsnU15bavL7rvpml5\nnfFZwH3u/mC/wMyOAh5191kzO5GeIf5xlc71ZiyEmBhaXmd8IXD9QNmrgb80s33AfuCt7v7YQVeW\nQMZYCDExtJmbwt3fFCj7IvDFJvqXMRZCTAxtL21rExljIcTEMLHh0GZ2KHAnPbfLIcCX3P0KMzsS\n+ALwPGAHcEFVnaQMdR10oV09QoEeMBfsMc+pV3QeClRICfoo6wBr0unVZFBHyk4f/fKqQSFZX/lA\njr3L55bJz2Y/xL6jDsJOoyKnU7j+YAdbzClV1YHXZNBFUw662L2mPKuiYJimGecUmkODPtx9D3Cm\nu78MeClwZrao+XJgk7ufBHwt+yyEECNlliWlj64RjcBz9yez00Povbc8CpwDXJOVXwO8vpXZCSFE\nAuNsjKPv9Gb2NHoLmJ8PfNLd7zWzo919d9ZkN3B0i3MUQohSdNHIliVqjN19P/AyM3s28PdmduZA\nvZuZNzWhKrtGhyTdop08DpQV7OSxNKTzxvTjthL1lE3Ek6IJt6Epl9npo6y+nN99I5LoJ6QP19Vh\n6+qoMc24SGeNzXUc76VQH9+fW/HQcHacNpMQtU1ptdvdf2VmtwKnA7vNbIW7P2RmxwAPF113B3PP\n+wTglDqzFUJMDPum/i/r//HeRvuc2DfjLNRvxt0fM7PDgHXABuBm4BLgI9m/NxX1cWZsECHEomTZ\n2t9j/R+cf+Dzhg0bavc5scYYOAa4JtONnwZc6+5fM7MtwEYzewvZ0rZ2pymEEHEmdp2xu99DLkNR\nrvyX9JJmVGYfcNiQ+pma9UX992943wwclumT89YZ989nie/uHEouHzqP7XgM5dcU19Vvm7w+ovmm\n6Nf9XZvzOzbHEsHHdmpO0VRj63BD62RD65TT6sOacCzRfd11xilzidUH9e2cJjwzk9XvySU1eiqn\nGTfs+h/ndcYjm/kwQ7og9YdGGiykHyA2Vt25tn39An6L8kalWv3wm439Z479GVy3fiGpO9fYW2jf\nEC8kXXq+qYzvrxEhhBhAxlgIITrA3mklChJCiJEzOzO+Jm0sZh6bZDTAI3+edbas6K+Zssl/itqG\nnFYpgRB1nXYxp1nIgRYLOqnrFMy19UAgB8wl/ZleXn73ipREP/3rYoEMRQltQk7BKk6zMrtfxB14\nVRyU5RP5lA3a2Jtzys3m9OG+g27/3twPO68f77G586YdeC3p1Gb228B/Bw6nt4Lsjf197szsCuDN\n9Nz+73L326qMod2hhRATw+zMktJHIp8G3uvuLwVuBN4DYGYn09sB5GTgbOAT2VLgZGSMhRATw8y+\nJaWPRFa7+9ez89uBP8nOzwWud/d97r4D2A6sqTJ3GWMhxMSwf3Zp6SORe83s3Oz8fGBldn4ssDPX\nbidwXJW5d0IzLpJMh7WNXbOsoL4f4BEM9Mg3nqH87s2xQI5Y0EcskU9Eh01KTp+SvCfUf4p+HEj6\nE9u9uWrwQZXkOEk6aUBTTknuU3esmH4d188Pro8GguRWJoSCNqKa8B5yZQXnJ9MsNTRjM9sErAhU\nXUlPE77azP4jvXQQ00O6qpQ4rRPGuJPEAh3EHHpWpRnnPdrGgj1DTNq3pmDzVGG1u6+L9P6HAGZ2\nEvAvs7JdzL0lAxyflSUjYyyEmBxmhtSdvrZ39Pnr8omJzOy57v5I5px7P/DJrOpm4Doz+xg9eWI1\nsDlpzhkyxkKIyWGYMa7HRWb29uz8Bnf/OwB332pmG4Gt2eiXurtkCiHEIqclY+zuVwNXF9RdBVxV\nd4zOGuPYjh9FBHf6iN1lbEfnWFa2WH1K0EeCUyya6S0WFFI3qCRQXzWoIxQUkZIprYqDLWX8ug6+\nsg7GlL5iWdti95J30OUj16b3ZA7CnNOOPTmtO+Sg25c738vB9fnzp2iPffEmXaWzxlgIIZKZHfUE\nqiNjLISYHNrTjFtHxlgIMTnsiTfpKmNhjEOJgIqSA0V3l85ktmCgR/68TKKgskEdsb5imnAsEVDV\noI2yQSGRQI+8Tly0U0csqKPu7hyxQIZYop+mNOW691e1r+BccztuhJL65HfcmBfA0V+rW6TzxjTh\nmUh9mwZTb8ZCCNEBZIyFEKIDyBgLIUQH0NI2IYToAFra1j1Cmdrmnafs3hFzwMWCOlKypqUEoFRx\nwMV28khwAPYzsfUDOmB0QR2lAx0qOvja6D92r7GdQEK7b8Ccs246F6gRDODI77gRcrDFHHB169tA\nMoUQQnQALW0TQogOoDdjIYToADLGw6kStFE1N8+B/vM7eeQ012WxDuomCooFfYTOF1LzrRIUEtCJ\noXzyn955Fc23uaCOmOb7FId1Yn758sK5TB+cyCcYwJFPtB4K4NhLuH7PQLuitlWDQqQZB9GbsRBi\nctDSNiGE6ABjvLRNu0MLISaHPQlHAmZ2vpnda2azZnbaQN0VZna/mW0zs9fkyqeysi3ZcdSwMfRm\nLISYHNrTjO8BzgM+lS80s5OBC+ntc30ccLuZrc62XnLgz9z97jIDLD5jHArqqJJpraivWH2or7Yd\ncBUdhH5479+6mdh65cMdXHWDOtp2wFWZ35M8Pfn6eW1zO3GEnHXBTGtQ3gH3RKAsf54S1BEaq2j8\nMdzpw923AZjZYNW5wPXuvg/YYWbbgZcDd2X1B11QhGQKIcTkMJtwNMOxwM7c551ZWZ9rMoni/bGO\nFt+bsRBichkmU+ycgl1ThdVmtglYEai60t1vqTCbN7r7g2b2DOAGM7vY3a8taixjLISYHIYZ4xVr\ne0efzRvmVbv7ugoj7gJW5j4fn5Xh7g9m/z5hZtcBa4DJMcaxoJEDkm8oOVARKTpwSoBIWU03RXNu\nIyjk8Lmi0K4dIZ0Y4jrsUznNtI2giScznbjM9THNtq/vVt/Jo9yziAVywJw+HN2Jo0nNN3TeVv/5\nfptmYdYZ53Xgm4HrzOxj9Bx4q4HNZrYEeI67/9zMlgGvA24b1mkpY5x1/G1gp7u/zsyOBL4APA/Y\nAVzg7o8l3pAQQjTL3niTKpjZecDVwFHArWa2xd1f6+5bzWwjsJXee/ml7u5mdijw1cwQLwE2AX8z\nbIyyb8bvzgZ7Zvb5cmCTu3/UzC7LPl+eeH9CCNEsLS1tc/cbgRsL6q4Crhoo+w3wOyljRFdTmNnx\nwB8Bn2bu9fwc4Jrs/Brg9SmDCiFEK+xLODpGmTfj/wa8B3hWruxod9+dne8Gjm56YkIIkcwYh0MP\nNcZm9sfAw+6+xczWhtpk+ogX9XEHcw62E4EXl5xYkYMuel12R8FMbRAO6iBQHwv0yJcXOcjKBpCk\nBHXEsr6lBIVkjruQ0w7mHHchpx2UD7TIl6c5+A6+PhRIUXx9ggMtshNHY9dHAjmgIOtayJmW4kAL\nXZ8SFFIlqOM3kfn9Zor166dolAnO2vZ7wDlm9kf0/hs/y8yuBXab2Qp3f8jMjgEeLurgTMj5vIUQ\nIuPwtaxfv/bAxw0bNhS3LcsYG+OhmrG7X+nuK939BOANwD+4+8X0lnNckjW7BLip3WkKIUQJJlwz\nztOXIz4MbDSzt5AtbWtyUkIIUYmWlrYtBKWNsbvfCdyZnf8SOKuNCZXVh0O7exSSsuNy2d058n2k\nJBJKCfooG9RRNRFQdh5LBJSi06ZpugcHcMR2wggFeuTbhvrMt00ZP7a79ZPzNO/hOzY/+USmGccC\nOaBa0EVMU87ru78JlIU033xZimbdbzuKREFjLFOMXQSeEEIU0kH5oSwyxkKIyWFSl7YJIcRYIZlC\nCCE6gIxxPVKSqkXbLp3/b2FnRYEYoYFiWdlizr6qQRtlgzpCTr98fWD3DohnZYtlMivrdCtqG3KQ\nFTno+m2fCjjN5l9f5EAMZX0bvhPHU7m5RHca2Z8FwOQCOab35OaXOfCCjrreYHP8JlAfcpAV1YeC\nLmJBG6EAjSIHYcjBlzL/MdzpYyHohDEWQohGWAxL24QQovNIphBCiA4gmaL5gWP1od095tXHhOaU\n3TWK9OPY7s+hoI5YfSwoJBbUEamPJQKK6aixQIpQIES+bUzTDSUPyreNBZWEdN6i61M066B+nkv6\n81SmCeeT/5DTjHkiyz5bFAgRC6qIJfoJ6cMp/T8emFddTbnwLTXfccOZa1pa2mZm5wPrgRcCa9z9\nO1n5OuBDwCHANPAed78jqzsd+Dt6/wO/7O7vHjaGdocWQkwOMwlHGvcA5wH/yFxaCIBHgD9295fS\ny9OT3+Puk8Bb3H01sNrMzh42gGQKIcTk0N5OH9sAzGyw/Lu5j1uBw7Ktlo4Cnunum7O6z9HbhOOr\nRWPIGAshJofRasZ/AnzH3feZ2XHAzlzdLnoblhYiYyyEmBxqvBmb2SZgRaDqSne/JXLtKfSyWa6r\nOn7njHGK4+5AWe6iZWUDOGJZ22KBHvnyWH1RVrUqQSH5eYV28gjUz+TqQwEesZ06YpnSipxisZ08\nQs64mIOu2ME2PKgktFNHLKhj3lyyDGz97GtQkIHtiYKgjligRhUHW6gM4jtt9K+rGrQRysoWNIL5\nBr8uKF8VurAlprIjjLtXMqTZPqFfBC52959kxbuA43PNjs/KCumcMRZCiHZYmx19Ku8sckA4NrMj\ngFuBy9z9m/1yd/+Zmf3azF4ObAYuBq4e1qlWUwghRAQzO8/Mfgq8ArjVzL6SVb0DeD7wATPbkh1H\nZXWXAp8G7ge2u3uh8w70ZiyEmCja8eC5+43AjYHyDwIfLLjmO8BLyo4xFsY4qBMXtF0a0nznNQiU\nhfThWKBHvk2RZhvrKxTUkbKTx5Lh9X2teO/yuT+AQjs9FwVKxBIBPXlAhx2uE+fHqBtAUpQoKKQJ\nh/pPSkSBFZe7AAAHrElEQVQUC+p4Ih/g0S/j4LL8eZHOG7quqK/YThohfToUlFFVUw7qwyFNuEgn\nzj+EVaHOajC+8dBjYYyFEKIc4xsPLWMshJgg2szP2S4yxkKICUJvxrUGTNKEQ2UpiYD657MM3515\nL3D4kPr8eWx355R1yLHk9SF9ONen5+r7WvH08rw2OnxH5elA8px8eWxNcWwdb8o64ZR1yHH9engi\nolBy+Ccfn5tr6XXE+fKY5lv3+lgin5g+XTnRz68HGubL8uWPB8oG2zaNNOPxY3mk/vBI/SIib/RC\n5A3comdPpP6JSL2oid6MhRCiA+jNWAghOoDejIUQogNoNcWCUzTxZaGKWPKgKomEUtrGEv0UOeiW\nDLQbPM/a5p12+Z08+o67+TtahHZMLp+8p4rTrtfHYfOuqdpX3gEYcsalBLA8OZ1zJmaOu6dyiYCS\ngjqqONWKdm8OJeoJOeNic0mZa+lAjnx5USBHWQdfG0imEEKIDiCZQgghOoDejIUQogPozbgxQsEe\n0YTzdTXh0HX5oI+YJlyk+caSz4cS/YTaRjTn2dw1+eTxc5pw+YTroYTusURAMZ12/ljhvkK7L8cS\nBVXRt598MndNbvfm6b5WnJIc/omE+pDmHAq0gLhmXEVTzteXTgQfS/RTVVOWZhyic8a4MyjoQ4gx\nRG/GQgjRAcZ3aZt2+hBCTBD7Eo7ymNn5Znavmc2a2Wm58jW5HT6+b2YX5uqmzGxbYAeQIHozFkJM\nEK1pxvcA5wGfCpSf7u77zWwF8AMz+9/uPgs48GfufneZAUoZYzPbQU+VnwX2ufsaMzsS+ALwPGAH\ncIG7P1amPyGEaIfWtl3aBmBmg+V5XeQw4FeZIe4z/4IhlH0zdmCtu/8yV3Y5sMndP2pml2WfLx+8\n8APuZecy0UxNTbF27dpW+s7/EJ9VcN4l2nwWUZ5ecH7kQk+kx0ifxTwOC5wfPYqJ1GThV1OY2Rrg\nb4ETgIsGqq8xs33ADdl+eYWkaMaDFv4c4Jr+gMDrE/padExNTY16Cp1Bz2IOPYumGaYR3wf8n9wx\nHzPbZGb3BI7XDRvR3Te7+ynAacBfmdmzs6o3uvuLgTOAM8zs4mH9pLwZ325ms8Cn3P1vgKPdfXdW\nv5vx/DUqhJgohr0Zr2L+Bqi3zat193V1Rnb3bWb2APDPge+4+4NZ+RNmdh2wBri26PqyxviV7v4z\nM3susMnMtg1Mws1MeoQQYsQsyNK2AyqBma0Cdrr7jJk9D1gN3G9mS4DnuPvPzWwZ8DoGrf9gp56o\n6ZrZB+jF9vw5PR35ITM7BrjD3V840FYGWghRGncv7fAapIq9KTuemZ0HXA0cBfwK2OLur82kh8uY\n00L+wt2/amaHA3fSCypeAmwC/r0PMbhRY2xmTweWuPvj2QC3ARuAs4BfuPtHzOxy4Ah3P8iBJ4QQ\nIk4ZY3wCcGP2cSnwP939Q9nSto3Ab6GlbUIIUYtkmUIIIUTztBIObWZnZ2GA92drkCcSM9uRhUBu\nMbPNWdmR2RKZH5nZbWZ2RK79Fdkz2WZmr8mVn54tobnfzP5qFPeSipl91sx2m9k9ubLG7t3MlpvZ\nF7LyuzLnSCcpeBbrzWxnLhT2tbm6iXwWZrbSzO7IwoZ/YGbvysoX5fciGXdv9KAnVm+nt4ZkGfBd\n4EVNj9OFA/gJcORA2UeB92bnlwEfzs5Pzp7FsuzZbGfuL5PNwJrs/MvA2aO+txL3fgZwKnBPG/cO\nXAp8Iju/EPj8qO858Vl8gJ7DZrDtxD4LYAXwsuz8GcAPgRct1u9F6tHGm/EaYLu773D3fcDngXNb\nGKcrlA2GORe43t33ufsOel+8l2crUZ7p7puzdp9jDAJo3P3rwKMDxU3ee76vG4B/0fhNNETBs4Bw\nKOzEPgt3f8jdv5udP0EvyuI4Fun3IpU2jPFxwE9zn3dmZZNIPxjm22b251lZUTDMsfSeRZ/+cxks\n38X4Pq8m7/3A98jdZ4BfZU7jceKdZvY9M/tM7k/zRfEssvW3pwLfQt+LUrRhjBeTR/CV7n4q8Frg\n7WZ2Rr7Se39LLabncYDFfO8Zn6SXq+BlwM+A/zra6SwcZvYMem+t73b3/BYf+l4MoQ1jvAtYmfu8\nkvm/5SYGd/9Z9u8j9Jb/rQF2Wy+VHtmfWw9nzQefy/H0nsuu7DxfvqvdmbdGE/e+M3fNb2V9LQWe\n7fMTVXUad3/YM4BP0/tuwIQ/iyza7AbgWne/KSvW96IEbRjjbwOrzWyVmR1CT2S/uYVxRoqZPd3M\nnpmdHw68hl5u05uBS7JmlwD9L+TNwBvM7JBs7fZqYLO7PwT82sxebmYGXJy7Ztxo4t6/FOjrT4Gv\nLcQNNEVmdPqcR++7ARP8LLJ5fwbY6u4fz1Xpe1GGNryC9P5s/yE9Qf6KUXspW7rHE+h5gr8L/KB/\nn/SSMd4O/IhetOIRuWuuzJ7JNuAPc+Wn0/vPuh24etT3VvL+rwceBKbpaXhvavLe6W3DuhG4H7gL\nWDXqe054Fm+m53T6PvA9esbn6El/FsCrgP3Z/4kt2XH2Yv1epB4K+hBCiA6gPfCEEKIDyBgLIUQH\nkDEWQogOIGMshBAdQMZYCCE6gIyxEEJ0ABljIYToADLGQgjRAf4/0dX4TvoTrZoAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3d713faf28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "psi_int = integrate(psi_c)\n",
    "plt.pcolormesh(xx,zz,psi_int)\n",
    "plt.colorbar()\n",
    "plt.axis([0,L, H,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}