Invert Laplacian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalised inverse Laplacian\n",
    "\n",
    "This demonstrates a manual method to compute the inverse Laplacian of a scalar field using bits and pieces we can find lying around in windspharm. This is obviously a massive hack but it demonstrates that it is possible!\n",
    "You can do all this just using `pyspharm`, which may be a better choice, but here we'll let windspharm do the work of setting up the required pyspharm objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data loading and reference results to show correctness\n",
    "\n",
    "I'm going to use u and v wind as one normally would with windspharm, but these fields can be anything you like (omega even); in general you just need some fields with the right shape to initialise everything. I will use Iris to load, you can use any interface you like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import iris\n",
    "\n",
    "u = iris.load_cube(\"uwnd_mean.nc\")\n",
    "v = iris.load_cube(\"vwnd_mean.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a VectorWind instance from the u and v wind components, then compute the vorticity and streamfunction.\n",
    "I'm going to demonstrate that the manual method of computing the inverse Laplacian gives the same result as\n",
    "the built-in version, which can be verified if I use vorticity as the scalar field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from windspharm.iris import VectorWind\n",
    "\n",
    "vw = VectorWind(u, v)\n",
    "vrt = vw.vorticity()\n",
    "psi = vw.streamfunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute and force result back into a cube\n",
    "\n",
    "We're going to need two things to compute the inverse Laplacian, a `Spharmt` instance for spectral transforms and the (private) `invlap` function from pyspharm.\n",
    "A `windspharm.standard.VectorWind` instance has an attribute `s` that holds the `Spharmt` object we need.\n",
    "The metadata interfaces are privately holding on to a windspharm.standard.VectorWind instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the Spharmt instance\n",
    "S = vw._api.s  # (just vw.s if you used windspharm.standard.VectorWind)\n",
    "\n",
    "# Transform your scalar field to spectral domain, you can't use metadata here\n",
    "from windspharm.tools import prep_data, recover_data\n",
    "scalar_field, shape_info = prep_data(vrt.data, \"tyx\")\n",
    "scalar_field_spec = S.grdtospec(scalar_field)\n",
    "\n",
    "# Compute the inverse Laplacian\n",
    "import spharm.spharm\n",
    "ilap_spec = spharm.spharm._spherepack.invlap(scalar_field_spec, S.rsphere)\n",
    "ilap = recover_data(S.spectogrd(ilap_spec), shape_info)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inverse_laplacian / (m/s)           (time: 12; latitude: 73; longitude: 144)\n",
      "     Dimension coordinates:\n",
      "          time                           x             -              -\n",
      "          latitude                       -             x              -\n",
      "          longitude                      -             -              x\n",
      "     Scalar coordinates:\n",
      "          air_pressure: 200.0 hPa\n",
      "     Attributes:\n",
      "          Conventions: CF-1.5\n",
      "          history: Wed May  9 11:39:18 2018: ncatted -a actual_range,time,d,, uwnd_mean.nc\n",
      "Wed...\n",
      "          parent_stat: Mean\n",
      "          statistic: Long Term Mean\n"
     ]
    }
   ],
   "source": [
    "ilap_cube = u.copy()\n",
    "ilap_cube.data = ilap\n",
    "ilap_cube.rename(\"inverse_laplacian\")\n",
    "print(ilap_cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check the results are OK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f487d90e860>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAAEDCAYAAADtKBX8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJztnXnQNVdd57+/e983ewgkkSUQyMgSAhQgICCLZCC4gAR0BmRRQhVIMVNTU1CgUMCwCSNSDmMpggwgQWeEAVGWsIOTAsQARgkYJSSQQEIWDBBCIGR5nt/80af7dp8+a9/uvud57vdT9b7PvX3WPuf06e/5nV/3FVUFIYSQMlhsugKEEEJWcFImhJCC4KRMCCEFwUmZEEIKgpMyIYQUBCdlQggpCE7KewgRuURETttQ2WeKyKtHzO8jInLGWPmNgYioiNxl0/VYFxF5qIhcKCLXicgTZiz34SJywVzl7Ve2elIWkbNF5Psicqh1fNQJiPRR1V9W1XdMlb+InCoil02Vf6ucZ4jIZ6cuJ5NXAXiDqh6lqu+bqhD7Jqaqn1HVk6cqb1vY2klZRE4C8HAACuD0jVZmZkTkwAbLFhHZqnEnIsuZi7wTgPNnLpOMhapu5T8ALwPwdwBeD+Cs1vFnA7gJwI0ArgPwQXP8EgC/DeDLAH4E4G0AbgPgIwB+COCTAG7Vyud0VBfGNQDOBnBKK+yFAL5t0l0A4FHm+CsA/BWA/2vC/hHAfVrpLgHwAlOHH5h4h7XCfwXAl0yZnwNwbyvtC03aGwAcAHACgPcC+DcAFwP4r4H2OhPAnwD4kKnb5wHcuRX+EABfNPX6IoCHtMLOBvAa097XA7iLOfYsE36eaev6nwI4NaEdne0B4EhTzm4rzxMAPBDA35u8rgDwBgCHtPJTAHfxnP8zAHzDnPvFAJ4G4BQAPwGwY8q4ptVWbwLwYVRj5TQAhwL4AwDfAnAVgD8FcLiJfysAZ5l++L75fAer/V5t+vQ6AB8EcByA/wPgWtPeJ5m4Xzfnfb2Je6hpp9Na+b0CwP82n08y532GqdvVAF7SirsE8GKT7w8BnAvgRACfNul+ZMr5dQCnArislfYUU/drTB+enjqetvnfxiuwsRMHLgLwnwHcH9UkfJtW2JkAXm3FvwTAOagm4tsD+A6qSfNnzMD/WwAvN3HvZgbrowEcBPA7prxDAJwM4FIAJ5i4J9WD0VwsNwH4jybdC1BNAAdbdfgCqgnmWAD/CuA5Jux+pk4PMhfSGSb+oa20XzIX1OGoVknnoro5HQLgp1FNOr/oaa8zAXwP1cR2ANWE8C4TdiyqyeQ3TdhTzPfjTPjZ5oK/pwk/iNakbJXzbABfBXCLUDsmtMepaE0Q5tj9ATzY1OEkE/+5rXDnpIxqkr8WwMnm++0A3NN8fgaAzzra6gcAHmra+TAAfwjgA6aeR6OaWH/PxD8OwH8AcIQJew+A97XyO9uc950BHAPgXwB8DdVkfwDAnwN4uzVWTwt8fwX6k/JbzLi4D6qb9ikm/LcBfAXVuBUTfpyrvdptbvrrIlQT+iEAHolq8j251UbO8bTt/7ZqGVkjIg9DtcR7t6qei0oFPDUh6R+r6lWq+m0AnwHweVX9J1W9AcDfoJqggUo1fEhVP6GqN6FSSIejUpM7qCbxe4jIQVW9RFW/3irjXFX9K5Pu9agu6Ae3wv9IVS9X1e+hurDva47/FoA3q+rnVXVHK3vtDY60l6rq9QB+FsBPqeqrVPVGVf0GqgvzyYHz/2tV/YKq3ozqIqrLfiyAC1X1L1T1ZlV9J6qJ9XGttGeq6vkm/CZX5qZfXo1KUV0bacdYe/RQ1XNV9RxTh0sAvBnAIwLn22YXwL1E5HBVvUJVY+aB96vq36nqLqp++C0Az1PV76nqDwH8d5i2VtXvqup7VfXHJuw1jnq9XVW/rqo/QLU6+7qqftL0xXuwGntDeaWqXq+q56FaudzHHH8WgJeq6gVacZ6qfjchvwcDOArAa834+ltUK4CntOL4xtNWs5WTMioV+XFVvdp8/0tzLMZVrc/XO74fZT6fAOCbdYC5MC8FcHtVvQjAc1Gple+IyLtE5IRWPpda6S4z+dVc2fr841aZdwLwfBG5pv6HShU78zbxT7DivxjVSsCHr+zO+Rq+iWpF4Sq7h4icCODdAM5Q1a+58m23Y0KdXGXcTUTOEpErReRaVBPj8aF6mXJ/hOoG8RwAV4jIh0Tk7pFk7fP9KVQq+NxWW3/UHIeIHCEibxaRb5p6fRrALS1bdOrYG4qvHU9EJVpyOQHApabPauwxkdx328TWTcoicjiAJwF4hLk4rwTwPAD3EZFaHaz76rzLUU16dZmCanB/GwBU9S9VtVbrCuD3W2lPbKVbALiDyS/GpQBeo6q3bP07wqjWGrXiX2zFP1pVH5N3qgCs8zXcEeZ8HWV3MH3yPgB/qKof8eVrt2MEV3lvQqXg76qqt0B1E5KEvKCqH1PVR6MyXXwV1arCV459/GpUE+c9W219jKrWk9DzUZkHHmTq9fPmeFLdEvgRqptCzW0z0l6KymySy+UATrQ2de0xQRxs3aQM4AmoTAj3QLVcui+qDYnPAHi6iXMVKhvrUN4N4LEi8igROYjqorsBwOdE5GQReaRxw/sJqot1p5X2/iLya8ZD4rkm3TkJZb4FwHNE5EHGw+FIEXmsiBztif8FANeKyAtF5HARWYrIvUTkZwec74cB3E1EnioiB0Tk11G171mJ6f8MwFdV9XXWcW87JuR5FYDjROSY1rGjUdmGrzNK9z+lVE5EbiMip4vIkab867Dqs6sA3EFEDvGlN2rxLQD+p4jc2uR5exH5xVa9rgdwjYgcC+DlKfXK4EsAniwiB0XkAaj2LFJ5K4DfFZG7mnF1bxE5zoSFrpPPo7oZ/I4p91RU5qx3DTuF7WEbJ+UzUNnnvqWqV9b/UO3EP81Mhm9DZfO9RkSy/TxV9QIAvwHgj1GppMcBeJyq3ojKnvxac/xKALdGpdhq3o9qqVxvnP2azwZrlfkPqOyWbzBpL0K1CeWLv2PqdV9Um4lXo7oAj/GlCeT1XVSeH88H8F1UG3K/0jIPxXgygF81DzvU/x4eacdYnb4K4J0AvmH68QRUG6dPRbXh9BZU3hopLMy5XY5qc+oRqDaJgWqD93wAV4pI6HxfiKpPzjEmik+iUsdAtQl4uDnHc1CZNsbkv6FSu98H8EpU5rpUXo/q5vhxVDe0t6GqK1CZ4N5h2vdJ7USmj04H8MuozuuNAJ5u+oUEEFW+5L4UROQVqHazf2PTdSGEbIZtVMqEEFIsnJQJIaQgaL4ghJCCoFImhJCC4KRMCCEFMfnbwg6XI/Un+PHUxRBCyFR8U1VPmquwyW3KIqKnLZ4Uj1g6nadFyaxs15s+SWF8cvfdUNWxnq6MMst7dWUx2/k40d0RbjyxiYGTdhhOrCTCpucJLzNf2rxSCCGkIDgpE0JIQcxjvljm/xrOKCaHpvxe5tl5ROtjL8+3wZyxhklikqVqKSaSwvp+zGsphVnMEKX09QTs3zMjhJA9yDw/oDngrjbqT032lMvSCo4riaY+ERXU5BU658KUVI/E/spSRFPkmcOsG7Xxwbu2es2o7zptWtezVPU7qF655UTf0TguVMqEEFIQ8yjlpXvul0X8nqC7IyiYiIuhiEe1ONXI0gS508jCXd9OfN+dekoFnaEOouojkJc3baz8oUpMxlJw9p7AiHZYx1jxrgQTx4DuZuipNcbVpm3/64xFAGnjarQxNA5UyoQQUhCb8b5IUMjetDY5StqjflQ8ebgUdm1j89qY3Urap6A7xeWon0ySFI9HdWSpX29c9/GU1VIxSmagevaOL2fenvFuj6d2tOg+xzyq2ssY+wm+PDLGW3SsFTLOqJQJIaQgZrIp20p5xDuSba9O2dVWW3WY+tmq26GMeqrHVtNeJd1XQEPU9GBy7MB23HXUiB3Hp1ZyVMoUds4h3hCJyllSzs23iqvH5NLKoxPfGls+VZ2ggqdcrdkM2n/ojdXu9+j46+Q1wlicACplQggpiM0+0TeF4nEVFVFBWiuIWnW74mutgCOq2gQ3CifLLp3vP+1jLTtwrvpoh/fCrLxybYMu1lEyPnXr27oY0AdqK1KX95Ev39RVXCBN8h5JyCsk1R8/g0GeFJFV2qCVWWSMdqPOr1uplAkhpCA26qfcMLUNxy7eViN1BEvtdjDKIKqqa0Vtn1O7zBQ1jYD/dA5D7ME+teFRGE41UZdrp0m1MUfqOpiYHdhWpKHqefKSFJ9ne4x5lGdvvIXip67m6q8uRe3zAply/yNHjcYUcopNOTRuW+GbgkqZEEIKYqZ3XwxQR760U+BTT22lUYtFjaggW1GHFE1ITSPRvzVC0m60V3V4FIWtgkN5+pSyr1+nVikxW25obKb6KQ/xZ/bYjKVZVrlXZN08rPFj24VttdvOw6Oqe94fiUp6LVxjwDdeYgrZkVczniN7KcFyJ4RKmZBtZ4xXGZDRmEcp+1RTTc7daM4713LpVz318Z5SrgZ4o6gDiiZLybTyTiLmg5miIHwq19WftgL2qZAUH+h1+zikVHt7AD5PhQzlN4ZKXHoUsY0r3DcGbZXtWdVVQe6VXTMGm/0PS427/KjHfG+IzxMoqpyldSiy0vOUVcWl90VZTPyjsvuGMW+UhTzquk30xME+YhMubeuy92pMCCH7mHnNFzbBRyBjjubjv2A9iK0mfEtB+1xdS0fbxGHHjblHZZC8dGsfs5eIPtOEMw+PCcS3wZey2RIjx1yxsM0XGa/s9JkpppA2sb4ObfT5TCH2WPSNw04cz8a1b7PaUe46r99NfjgktjndCvOPxYQHTWaASpkQQgpiQy5xGepohFf2DYrTI9E1qVES1sZIW2XZSiayaVjjVDQ+vC5ngU06+3tPWZi0S4+SdqTRJi9P/cbc1HOqR6tYX/+lPPY8q4TxPPcdUvKxsWjH23WMu5iajiloRz1irp9BYmo21X2zHTfmRscXEhFCCKnZzGPWqa+GBNLcVnzH11HTMXyqZGmFO23KHrcsW8GkukmlEHCF6ylhn+3NUshOFeyxIffi5vRjqjo0zdhRw5YNWXv9IZ20TTX2kueN9+Enzwqsp4bbKtcaizu7njQB109bRa+jlGtSH/EP2Y1jyjh0jcwIlTIhhBTEPEr5gFVMqvoFktWuOpVyOE1S+Taxh0lqbNUWVMox23JIKdoSL2KDD3ldWEpBa2Vhq1tzPEUpxxSys988ONsSWLX1olZoATuzr18W3fiKTMXeqegApTXF60gbVWsVFVTKdrvUSjlRQQNRO/SqqLhXRvSlQVM/4LQB1UylTAghBTGTUja20SnVbUBp9MpJLcOZmc9+Z2WplmoLPR7rsoW60gx5nNejBjptYqtYCStlr4Ju5yt2mn7cZCz12nlPk2pP5brszDU9Ne1R0k38AfbPJPU/VA7l2NxTx2TILn/Ao4x9Chro27J99mjbw6QOD66cEz0mXGo4de+E3heEEEJqZlHKeojbphxUFCPYg6MKOSGPfqZpqsRrH24f86VJ/R4ispJIUspL63vP5tzPWy2FXCvmHNuxzcoGav60PCpEpdWedWEB3+KearPK8nlrJFXUfY5Z5x6TSSl5edR/42+SoJRtVYsDi0641MddHh69sMSnYUPEnm2Iqd/259jq0bVHMiNUyoQQUhCzKuXo012utF6F7EkQymvEW1Dv/fP28/7Wzd9pm+y9I8BT1ho+s97VQoYdWHvqw0rXalc7L88vXSUhjQgT6/iqPXS3bSs14YGny3pqOvJTYcGwiIoKKuTE8dsbs2us6qQx/3ZXHkBfPavPS2UnbCeuPnvSNuEDVoAeP2XvCtC1h+GxGSelmREqZUK2nD31oMwWMItS3jncFNO7QwUSpd6hUkxsI9zt7IHbG8eegS19UdJnyosipc0jClgt5RzMs/6c69/psCs2yjcgrGSxUtKNYt5dle2zR4vdECEPjwjeFVhO2/eO12ni+yK+8d3/kQTruGvYmbB6tSG2jbnpE79t3uvlsYZnS7IHVWi/yvYisleIdl4Yd3WdCpUyIYQUxKxKOcu+uGFfwQbrHQm2avNJ4EYtpBQxxbl6suyqWo+Ss9Rur99C9uLUc7FV0rKtbus4nT+r4y3Fpdq2PZt4bRu3xx7dexDSXgl5v7QTuQ97907g6GtbsNveK73jrYOpq5GeDddk1T5s25/VjttVu061Xatszypynd8BjtnW+2O0bReu47jVdqxP5oZKOcR+trWVctOz6F3Qm2bE+ox68x3z8d8xxnmp10qZwzwIJ2VCCCmIWcwXNx3pfmH3KKpohDu0tx6O9XkvrtfxPXCLXkSWXglZ2CSbhhxqzWee8IevIT8i5iBXd67ME+aDZepYbQgac0s78a69zPaZodYYR9GltPnb6vdkk1CzKeXJM1CuWOe8imiHt8PCZozecdf495r03MdDRM87di1JIKyJU5acplImhJCCmEUp33yEZ+53qJPkjbRm0yBwl4vcmWNqya0kerlECnFkkbOZg0y3nHVE7IiCwd+PXZcz7YX380iJ4wv3KuMmTVe5h/JPVVp2P/r62xW3Od7b2IrUwZN/lSjsrukuyDq6xiZdLs7xPnD1NuaYngsqZUIIKYh5bMpHRGxaIQXks3/5bFyBON64HrWUk+cqkRXdoZL8Ssr67snTGSeTQfb8AW0etUk26VIqFFG9IbxpBjRkomqL9msoLFZGQr36hSTkkYtzbKa5mA3aS8nNy5WFd9WUeGxiqJQJIaQgZrIpmw899du/g3sVcuS4+4U/3a8xO6dNiiKL7qK7HNc9cbzHfWVlEDyX1HZKWLX4HuENrYqq8ASvkBgJ8aN5Znm8iDONr1+D9cjt64lspVH1GjtXBC6r6AojoYJW2ug1GtuHaMepx3ck/dRQKRNCSEHMo5QPq/6uXm5ivje+jqu4zV3Kp4h9eYR8inNsSDGSlbHnO9B6vNOKEymjjU/RRO1lLuVgxUlfnfTz7L/ox19ulNx2GMF2Gz2ekEeSg85Q2+gICjnswdH924vrG7uOfkhR1aMTWvWljufOy5XGrV4KVMqEEFIQ87yQ6NDqr3jUbkdpeZSwL41X5Vn5Or/HGKC4YooZaPlhRmzLTfwhDgI5duFMBTFbm0ds7E32HvXmPJZr18xRzAlpommHMND2HlpR1DRxYgq502/qD3N9t8sKEH2OIbSCtuaS2CrcWd4MUCkTQkhBzKKUdw8xH2J24pQ4KXn48kokySYZUwwBW1tUKafUJ0LMtpzU5rZCjtma22GReiS9+yFJla3aM2VVEit/kDIuhUj9gqsCj7rttZ944jvy7Zfnc73x1dhBTxlLN+ve+GsdaH4KzITZc4oV3okzI1TKc1H6BU22l+Lel7rdzKSUa99V951qkH0zpqDbjGBLjtrSUuyfvjS20vPVY4jjtKddXEo52csi4H2RWtWkdkrd6Q+1YxNmbUIU4BkwORnjqTdsFlYcn23e1dG5StlXP8DfZo1CNordtiHXU047vTXf2D83pi6bMpXy/qX4pW8hsJ3mh21eFvPalOsfY6zvTM1/cKgvcR4Pem4k7vz33uEbGpQxO2+KDTlmjwvlEalHc7w5J/vkHEmiyteokIgd35V/skJuEvQ/R32+F9pNE7Ih+2ylKSoudbKKqDkXc1gMss7Rp4B7fWApZ6D/gjYTpzkcU+wp9K5t90WsDhvyatzWqtpUYLdz2Ol9sX1P9AWWv7lMMsjXURCJm3iT18POasQ2H3XAhm5Iifg29QaxDepxhLFZ2PvhV9gT8h6C5gtCCCmIWcwXujS3q3ojxhx3OXirtdSwjfIxM0ebtVVhaBPKjhMxSTjzyHEz8uQVJW7NaLVx3ZbaiZzkCtcr16qsJ3KovWLt1MTxLbGBlezwbUzFFHrCW+79ZjO/nazXGhOsOvrHPSaJQJzmnUsOs4W/LHcb+5DAQGrGpG2ms20j9i9pW0O5k7Y+lXpsGFXdmMU6v5Y+/1KASpkQQgpiHqV8oKu8YCuz1p3VVmX2PTRZ1cEvWJJVScYGYPSR0oD7kT8vjzpJrWMnU/93u52S+6BJl+PLZGUSyiLWLnUcn/uWK40t4G0lGCS8Q9x7e2zGuJvC9hnd6LPUcChOf5ylq+G0tk1Lq71BYfWrdOeFrrttdzPQ/m2LlcpuD6D5jdJUyoQQUhCzKGUsfAqjL8V66tZOajmN2wlDNtN++SOQoviAiG3ZvtvH044qraz+sO33q3juw90OjBkQPfVOUMpOu7B4FLSVRnxtHSo/hu+Hdi076EqROepXt+mUpkvP6mB1wN9u6yjkJkV9za4xZrNtuwFbt5hAtY67Vrt8IREhhGw589iUranf+dNNvUT2d4+a88VvHbPvsj6VPQopdlA7booidMVzlTMIh/2tH5yTlZ8c6eE5p+4DOpZ9GSPZSFMq1N8S6ZbZKMR+Tp7f6vUTql5iHoN8in3eDk14en961a59OOt69A2SeF6N+dlO097/2YDDOpUyIYQUxDw25VS7awifKTKkmD226/gu+Rp3xyE208RdchfRqmYo015eto9oWi45kdMInUPb/pdTz5gCTEmbyagb+QOGqF1+014p3jNeZT+CkhyjXbyrPH/9ev3hWvlswKhMpUwIIQUxk1KOBCfdbC213fOVjSbxP+HXeyow4MqRS46P8YC7smW+7JP4ZF0w897x9CyyGWDabZLWjhXt8L78H1CpxCxifstjVmOM03C6KtWB1j6Mz36+IQatOgqpewwqZUIIKYh5lLJFlp+kL4/6Q9Ae3PXU6PtAh32e/fmORKpqTYgWq52G7IdDXj7eDh7BvplUVmxH3fYLduZheURMpThdee+FFxd7zqU5nLPS2gvn62LDr5ajUiaEkIKYRyn7dm6zniTyEPI3tR/esp+4ij7B1lfdyTfRqW+2np99arAqar89q5vU5z9qLTFiT4ZlICmqtia1LXfTKzRI+SVnbvvFp1Rkr7BH1a8Lr4vzZh/po1ImhJCC4KRMCCEFMe9GX281XD+cMGAN63UBC8TxrJn9LyxqPW6Z4oLXTpvjiTbA8T267O25wlnpWhXqvRjc3hQcw40uddMnYZMutchoflXq4NdBDOnPOSjFVLJOPdZpQp93pG3q60ThY9aEELLVzLTRN0Bx+fDtSblecJ764IC9Y+XYffK/rMd+YMDtXue64/Y231Jdq1LielUBPBFWkYJP3yLRpbGXKPaUS4DYxllog3bokJtSMefGmZq5N6XXIXUOCT05Hn1UXL0hc0ClTAghBTGLUpaeWtTW/y3bcudoExjOO3Dn9IbFhLHzxSS+CnWNzT77dEdBW+0xSDHHlLvHXuZM3zs3SzH7bP8pP1cVLHjV9ymucc5TUfTP2VVcrhqcQiJtWhWPqFg39nyFvTKN7S0NqmhgPpoBKmVCCCmIWR8e6SkvSzFXcTxy0fac8N0hB4gBr6J2Hq4fOPEZXK0PGT/gKBkPP9h5pYenSNL8agxPW6tyaX3rBOXjWlFEPDeSmmldebhp7wvDJCp3bkHpWbx5Xz/bPZhf0IxQKRNCSEHM633hs1W271xTvPdnqDQI1EV8irlWXj07dcctpDq021WHuugeD2wXOwpM9QNO8IQZxVsm0dPFit/eX1DrU88GPwDvNkOKe3Xua1ALUcaTsiGFvJ+hUiaEkIKY9+egZr7LpXpfrPNinZ4kbk61q6R115F2aeLWYXVcS/3KboZibr5GVJ0rq/qYz8vC/t7ED9SjV07XW2WVR/9py9WPElg2+JbtvbvWiLdPqtp2itxcib7hV0D6SLO7pjHI8WaCrZM0X/k1w2eCSpkQQgpi1p+DUluJpdyZmjgjeF3UplJvHvG7bU99+27h2lV+7XQ99WzdGu0XsDtexdFT1/2n2iLnIo7PVr+s+suqVy8PR1m91Ug3Tv8JR1ceHg8W2zy/W6tfS4UD/faxPF6ihFYUqaSUlZrnGGquV5+ECvps7jmLgdS9itA5Tql2B/zYxhRQKRNCSEHMopR1aT3FlaKYPQrZfsdFT7mGPCZ6Sq+vYoPxnXHcaWv7Z2fzvvn5KfN9YdlIaxUnXXUn9a2zZZfW+lijDi1baUzBOM7Nt5KRhc+G3O+/dDu+2xbfU8HtYy4howAW3Tbo2O+1Fa+dxxA7asZY65SVVYaVRU41U1dxa2S5ytv+nlJWTCGnK9VJHFsK2QOgUp4J78MmhBDSYp53Xxzo3VZNgCOuRz37j9sZ+O92tjqz0zpV3sC5tJ9MWz7MtTeB+b7oGnHVso3WSrqJh75ngq2Ym1JTNqU9q5DoagV2vHSlHMTyYGkO2+fqOS5Lh9z2eHKMoWaTaZflySM6Jge5O6xBRBGvusJR1kC7vXOFaq/eInkMqocL2pQLY8wOKWNlVD5sJzIme3A8cVImhJCCmMV8sTiwU31IXLJVB92bcL7Nt6xHqRPzcFZr4GZA16asrf8dcXxmjt1FKy6cce3vIa+1hoipKNl0BEy6WeLKWrVl5rDbonVsFd/aEO3trKVUJCFONmEzxXpmjBFJ2eDztmFgx7Yd3D4UM2cG0qbXqyyolAkhpCBmUcoHDt0JhgdfVD8gTSpDXOEWIyrB3cim3Eopm78tn7ioqu69KCnFvy+8OhllEy/GgOa1z7XzewKedrBXGt74CYwxFleZ2XnXH+w+iZe51usDPPjGaPfg0N3x/jlJMybdx0t5NHpMqJQJIaQgZlHKhx1y0+C0o6qQ5DLjcaZUzDUhVVJ/9qltW2Xn0FfKwx+uses1hNRzcMXz2u2tCL4y1ql3TU47xWylrr4YczXZexx+jbZvwlLLDtTHt3rzXYebmDfGgkqZEEIKYh6lfPDm0fLKUag+BTpG3lOSqpwBYBe2Uu5+j9mt2/gU3RyrAhe5KjWo1jLbJUUh5q5CctRbzOvI1Sf+Fc36K50cfP0Way/nKiBy3nWahUePT3WOU0KlTAghBTGLUj7i4I1zFLNvsNVvMG5P8bnT5q4agGlWDrF65KrP1PPyqbfclUYorzFJUcaruN0wWzXm5BXL20dKv6X2Vcq5+hRzTl5NnMIcmKmUCSGkIGZRyocfHO59sQ64ChzWAAAPKElEQVSl3QFrcpRwJ12GiozFdYWnKqgx7fop55Sr8pNUG2I25n4eqauSWLo2sbbMUYA9xRxR0LEyxmadNyWmtkOSMl5DVc8BlTIhhBTEPN4Xy8r7Yo470RB1PGW9hthym7QjPK60Tvk1pbXPmGns4yEF3Ys7ov3eR46qi6nG0LUxVD2OudKp2z7nGo7Vbx3lvCmolAnZckqblLadebwvDqzvfZF69wwNsIXsesMAYDnx4NzJ9iyY7p7ZVi+5F2WsHcdmnXbwqlnPKqSO7yozR1WnHE8hSyl7rpG+gu733xyeNrF+nMtDKDZ+N32TolKeidwJmZC52PQkRLrM432x7HpfrOMV4bvLhQbW0qsg8vMak2Rb24SKOYUUZVz6he1Xsd223UlQv3Uan9q2b8Bj9F+oD+wVnl8x71rf433mu3bsdnIRU8opeyZTipmpV8ZDoVIma1P6hEzIXoKTMiGEFMQs5osjlzc4jw9RWL4lxwLpmxdLpC3jljNvaO1kLHPHdLuqSe2PlHZJNVGV0sb2UtoVz27zHUvT9DcC4/3pW57Hltau8d6EJY77WPwcQuMxt51CJospxn1pUCkTQkhBzKKUb3HgJ2vn4VNUtiILKS97o8PexEjZ0PJtfIxByuZJzRSbf1GXwcwN0yrN+O5HsTxD2KqtxqfAXPHjm4Ph8HVIGX++/lgnrY/QOExtl97xQJ7+zdX9oy/3z5kQQsg+YBalfPQyXSnPodZsdeZTXkPsnTnufus8Rj2FMoidr+/cgqsTb9umuGOltX+OuoutMHxKGnC4uvVsoWn26jEIjbPUcRuyS+fisp/b7RGz27tXJWE7/iqe+/iQVcrGXVA3WjohhJAOsyjlY5Y/zk5TgmobYrtc5xHkMe/QtqIY81xS8srxkqniD7dLr/LOt0/HVGxYOdtKMKLCR/AcyHngIVUJD1kRpqzUYl4Vfa8Mh9qO5OFr8yF26ZS0c0ClTAghBTGLUr7lAKVsk6pAp1Bx7jzm9a9NxXeXX6e+Xv/WAXZN755AUr8N229IIcf22LSxlSTFLzm13DG8fMa0GfcINJdXvSb4glfp+5n3FbMvbbrPc2gV1OS3AdVMpUwIIQUxk1L+0RzFABhmVxyiIqf0V16LETf6Y0orbAce9vLx0EpnaJ7r4PRf9rTxmP7IYzDnGE0591wf8VC+qfb8UL3m9JbJgUqZEEIKYh4/5cX6T/TNQbHqdwLWUZVJdvvEF67H4qflGU06nIy8dwsbPpO2i2HIOffUq0TCneXanhzpT2TG8to0VMqEEFIQMynl8M9BlXan2lZy1fM66nZVZn7aHCWxnHNobeMwDpzzjrfrIz/GmlDsjoQbO6Tgx1DqU0KlTAghBTHP+5Rlx3m8uSPtYbtdCpu+87YZxf8143TWUbs+lbtMaM8S1MYyoub2OzsT7NHEVHRTZlDBr6/Up6SEsbuvKWlC3hY4qMlehuOXEEIKYhbzxWGRZVzOMsdtCAmzyeXIQXNum1bMY7r7jbnRtq4pItVEsJig/VPqvs0cHLF5ck0hu6H4pl4xM8amoFImhJCCmEUpHyr2I5HdO9RB1wtIPHcx++6bo4Kn2HjwlmUVNeTx7zGYwiVsiEJMvfunKN+Y6h1Wv+ENtZ829BYj6LTdCdamBz3Hh6jdev7xKfk55wkXVMqEEFIQsyjlg74XkQTuSPVdLHbXCtqOLBYz3ABrjTDrQwsJTHn3XUcp5ijUVAWclWdG3ddVkVPYtUtkgeXkZdTXfY57Zq3gY+NosWFbM5UyIYQUxCxK2adGXHcs20ZkqwtbGaeop53mrjqeUvEp9Ok1wvjMoeDW8VTIrd9U6je/HmVpnjHsxWOyju3Zvs52NJ6XreB91/AcL3IKUVYvEULIljOLUs66Q0d/8sf+SZm4/SemcHLs0r56bDNTKe1cW/UQJTik7mMo4NJU6ybIaYOYqvb1SUhBp/T9kLlhXTgyCCGkIGZSymG7cDdu5n1ihB8wHUP1lvp0UE2pvrRTKMYxlPuY9mCq4uGs6/Ps6scU+3PNJjxmOFoIIaQgZlHKNuvcfWyVXYoK2fSO7V6jNJ/dKT0lpnjCzcUU18JcdZ+DHIW8ScqY0QghhADYkFJeh9IUFtkf7BUVFWJnH6nabWY+pSxHxf+mxGFe5dWn1Lza/0qozzbkVVp91s1rA4hO7DUgssZv2RNCyOb5pqqeNFdhk0/KhBBC0uFGHyGEFAQnZUIIKQhOyoQQUhCclAkhpCA4KRNCSEFwUiaEkILgpEwIIQXBSZkQQgqCkzIhhBQEJ2VCCCkITsqEEFIQnJQJIaQgOCkTQkhBcFImhJCC4KRMCCEFwUmZEEIKgpMyIYQUBCdlQggpCE7KhBBSEJyUCSGkIDgpE0JIQXBSJoSQguCkTAghBcFJmRBCCoKTMiGEFAQnZUIIKQhOyoQQUhCclAkhpCA4KRNCSEFwUiaEkILgpEwIIQVxICXS8XJbvRE3rg6IQNoRpPmvEwfROL5jrS9WsPtYan6AhtK7ynfm0T8WKrM6JlW6aDxHWaH4I6Xppc89P+t40rkmhOuQegTDde1z652db0hJ+2OvRXqXiBmdgXxXeYgnb/sS6IW1g0XtajbxxTpfV5x2vH45VhqJhMOui/bOMRo/UFY7rp22m68VX4D2bOdOY4evjpz75Rs+pqq/hAySJuUbcSMetHg0ZFH3wKL6LEZoL6SaqBeLVS/Vn00akUUTzxxYxam/1/FlVU4vTStPiEBFVnq/Dq+PA1WYrPJo4re/S+v7os4XrXLRKadO00wYTXw07VF/98Wp81+F12U7wjvH0E3TCbfr5YrTz6OXJ1bx2+lDaUL16n1v1SNcLw2eSy8O7Dy0U0Yd3xVHfGlQh2uThzjSiGhrmHbTiDnWfBdg0aTpHwOAhZmUFk24dsNFTRxt4izM50Urj4UrHO3vu85jy04eu00eS/u7KBaovi+x2+SxlN0mz6XsVnk2Zew2+VR5rtIvzbEF6jx2m3KqeKs0ddl1GUvsNumaurbyXJr2qOu5+t4+N6zqCWApwBLSTDFLESwgWJoBsoCYY4vW90XzfXm7C49HJjRfEEJIQXBSJoSQguCkTAghBcFJmRBCCoKTMiGEFAQnZUIIKQhOyoQQUhCclAkhpCA4KRNCSEFwUiaEkILgpEwIIQXBSZkQQgqCkzIhhBQEJ2VCCCkIUe2/57UXSeSjALJfQUcIIVvO1bnvU06alAkhhMwDzReEEFIQnJQJIaQgtm5SFpEnisj5IrIrIg8IxPslEblARC4SkRfNWcd1EJFjReQTInKh+XsrT7xLROQrIvIlEfmHueuZQ6wvpOKPTPiXReR+m6jnUBLO71QR+YHpqy+JyMs2Uc8hiMifich3ROSfPeF7tu8Szm1Yv6nqVv0DcAqAkwGcDeABnjhLAF8H8NMADgFwHoB7bLruief3OgAvMp9fBOD3PfEuAXD8puubcD7RvgDwGAAfQfVzeA8G8PlN13vk8zsVwFmbruvA8/t5APcD8M+e8L3cd7FzG9RvW6eUVfVfVfWCSLQHArhIVb+hqjcCeBeAx09fu1F4PIB3mM/vAPCEDdZlDFL64vEA/lwrzgFwSxG53dwVHcheHmtRVPXTAL4XiLJn+y7h3AaxdZNyIrcHcGnr+2Xm2F7gNqp6BQCYv7f2xFMAHxeRc0Xk2bPVLp+UvtjL/ZVa958TkfNE5CMics95qjYLe7nvUsjutwNT12gTiMgnAdzWEfQSVX1/ShaOY8X4DobOLyObh6rq5SJyawCfEJGvmjt/aaT0RdH9FSGl7v8I4E6qep2IPAbA+wDcdfKazcNe7rsYg/ptX07KqnramllcBuDE1vc7ALh8zTxHI3R+InKViNxOVa8wy8DvePK43Pz9joj8DapldImTckpfFN1fEaJ1V9VrW58/LCJvFJHjVfXqmeo4JXu574IM7TeaL9x8EcBdReTficghAJ4M4AMbrlMqHwBwhvl8BoDeykBEjhSRo+vPAH4BgHMHuQBS+uIDAJ5udvIfDOAHtQlnDxA9PxG5rYiI+fxAVNftd2ev6TTs5b4LMrjfNr2DuYEd019FdXe+AcBVAD5mjp8A4MOteI8B8DVUO+Mv2XS9M87vOACfAnCh+XusfX6odvrPM//OL/38XH0B4DkAnmM+C4A/MeFfgcerptR/Cef3X0w/nQfgHAAP2XSdM87tnQCuAHCTue6euV/6LuHcBvUbH7MmhJCCoPmCEEIKgpMyIYQUBCdlQggpCE7KhBBSEJyUCSHEIvayISvuHUXk/4nIP5mXKj1mnbI5KZPiEZGnmcH+ZRH5nIjcZ428XiMil4rIdWPWkew7zgSQ+oshLwXwblX9GVR+5m9cp2BOymQvcDGAR6jqvQH8LoD/tUZeH0T19CIhXtTxsiERubOIfNS8L+YzInL3OjqAW5jPx2DNJxLpp0w2ioicBOCjAD6L6tWN5wF4O4BXonqZ0tNU9Qut+LdC9arE3ktrROS1AE4HcDOAj6vqCwLlXqeqR413JmS/YcbmWap6L/P9U6geDLlQRB4E4PdU9ZHmdQYfB3ArAEcCOE1Vzx1a7r589wXZc9wFwBMBPBvVY8dPBfAwVBPsi9F9/egzUb1/t4OIHIvqac27q6qKyC2nrjTZHkTkKAAPAfAe8+Q0ABxq/j4FwJmq+j9E5OcA/IWI3EtVd4eUxUmZlMDFqvoVABCR8wF8ykysXwFwUh1JRP49qkn5YY48rgXwEwBvFZEPAThr8lqTbWIB4BpVva8j7Jkw9mdV/XsROQzA8fC8DCylIEI2zQ2tz7ut77swwkFE7g3grQAer6q9l7qo6s2obMXvRaWsPyoiy9ZP8bxqyhMg+xut3vh2sYg8EWh+xqrecP4WgEeZ46cAOAzAvw0ti0qZFI+I3BHAXwP4TVX9mifOUQCO0OoVieeg+jWPHQAuZUNIEBF5J6qfczpeRC4D8HIATwPwJhF5KYCDqH4l5jwAzwfwFhF5HqpNv2foGpt1nJTJXuBlqN5+90Zjz7tZVR8AACLyYQDPQnUxvN8sHQXA81wZicjrUNmsjzAX21tV9RWTnwHZU6jqUzxBPTc5Vf0XAA8dq2x6XxBCSEHQpkwIIQXBSZkQQgqCkzIhhBQEJ2VCCCkITsqEEFIQnJQJIaQgOCkTQkhBcFImhJCC+P9QsAOfR/FQ2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import iris.quickplot\n",
    "\n",
    "iris.quickplot.pcolormesh(psi[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f487d816f98>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAAEDCAYAAADtKBX8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJztnXvQfVdZ37/POb9cucU03MItCimgFikqWLWWDqkNGR3aP6RgGYGmWjpTnVroyGCr1mnLoExtmZapFRUEwaK0knIHW5S2pFRaAqGCXBIkhktCDAEjufzep3+cvffZe631rMu+nfW+7/cz8/u9Z6+9bnutddb+rmc9ex9RVRBCCKmDzaErQAghZA8nZUIIqQhOyoQQUhGclAkhpCI4KRNCSEVwUiaEkIrgpEyOJSLyahH55zPkc5mIqIicmZjP20XkuVPrQwgn5VOKiNwoIlccuh4nBVV9uqq+5tD1IMcfTspkNaaqUUJOA5yUCUTkeSLy30Xk5SLyJyJyg4g8vTn3LBH5fSf+j4vINc3n85p0fyQiXxCRfy8iFzTnnioiN4nIT4jI5wH8qohcIiJvEZHbReQ2EXmfiGya+JeKyJtE5JamDj+WWf+va/K8pan/W0Tk4b3z7xWRl4rIB0TkyyLyZhG52Mjr+SLyByLyFRH5tIj8Pef8M0TkQyJyh4h8SkSu7JXxd5vPjxaR/yoiXxKRW0Xk10Xkol4eN4rIi0Tkw019/qOInJ9zreTkw0mZtDwFwMcBXALg5wD8sogIgGsAPFZELu/F/UEAr28+vwzAnwfwRACPAfAwAD/Vi/sQABcDeBSAHwHwQgA3AXgggAcDeAkAbSbm/wLguiaPpwH4hyLy1zPqvgHwq00ZjwTwZwD+rRPnhwD8HQCXArgXwCuMvL4I4PsA3B/A8wH8gog8CQBE5MkAfg3APwZwEYDvAXBjIA8B8NKmrMcDeASAn3HiPBPAlQC+HsATADwvfZnkNMBJmbR8RlV/SVXPAngNgIcCeLCq3gngzQCeDQDN5Pw4ANc0k/YPA/hxVb1NVb8C4F8CeFYv3yMAP62qd6nqnwG4p8n7Uap6j6q+T3cvYPl2AA9U1Z9V1btV9dMAfsnJK4iqfklV36SqdzZ1+BcA/ooT7bWqer2q/imAfwrgmSKyDeT1VlX9lO74XQDvAvCXm9NXA/gVVX23qh6p6h+r6scCeXyyiXOXqt4C4F8F6vMKVb1ZVW/D7mb0xNR1ktMBJ2XS8vn2QzMRA8B9m7+vRzMpY6eSf7uJ80AAFwL4YGOOuB3AO5rwlltU9Wu9458H8EkA72rMAy9uwh8F4NI2nyavl2CnpqOIyIUi8osi8hkRuQPA7wG4yJl0P9v7/BkA52C3KnDzerqIXNuYVm4HcFUv3iMAfCqjPg8Skd8QkT9u6vO6QFmf732+E/u2JqccTsokh3cBuEREnojd5NyaLm7FzlTwTap6UfPvAaran2AGryFU1a+o6gtV9RsAfD+AfyQiT8Nu0ryhl89Fqno/Vb0qo34vBPBYAE9R1ftjZ1YAdmaElkf0Pj8SO8V+az8TETkPwJsAvBy7VcJFAN7Wy+ezAB6dUZ+XYnfdT2jq8xynLoSYcFImSVT1XgC/hZ3KvRjAu5vwI+xMDL8gIg8CABF5WMwOLCLfJyKPaUwfdwA42/z7AIA7mk3BC0RkKyLfLCLfnlHF+2F3c7i92cD76UCc54jIN4rIhQB+FsBvNaaaPucCOA/ALQDubTY7v7d3/pcBPF9EniYim+ZaH2fU56tNfR6GnQ2akCw4KZNcXg/gCgC/2UzSLT+BnTni2map/h7sVKvF5U2crwJ4P4BXqup7mwny+7Gzrd6AnYp9FYAHZNTtXwO4oElzLXYmFJfXAng1dmaD8wF4nh2NPfrHALwRwJ9gZ6q5pnf+A2g2/wB8GcDvYmd2cflnAJ7UxHkrgP+UcQ2EAACEL7knJx0ReS+A16nqqw5dF0JSUCkTQkhFcFImhJCKoPmCEEIqgkqZEEIqgpMyIYRUxOJv7bpA7qNfw53piIQQUiefUdXL1ipscZuyiOgVm2cuWsYq6NGha3B6ES7oyOF4z9EboaqrPZG5yvttZXPYJ0z1aIYbT2pi4KQdhxMrSXDoecJk5a82vymEEFIRnJQJIaQi1jFfbL3X1iaZxeTQle9lXpxHsj7u8vw0mDMmmCQWWarWYiKprO/n/C7lsIoZopa+XoCTe2WEEHIMWeeHLEfc1fzfhJiAp1y2zum0kujqk1BBXV6xa65MSXlk9leRIloizxJW3ahND97J6rWgvlPatK1nrep3VL1Ky7mnvIgpUCkTQkhFrKOUt+G5Xzbpe4IezaBgEi6GIoZqCaqRbXMqnEY24foO4lt36iUVdIE6SKqPSF5m2lT5Y5WYzKXg3D2BGe2wgbFirgQzx4AeFeipCePq0Lb/KWMRQN64mm0MzQOVMiGEVMRhvC8yFLKZ1qVESRvqR8XII6SwWxubaWMOK2lLQQ+KK1E/hWQpHkN1FKlfM244PGe1VI2SGamezfEVzNsY7+546kdL7nOso6pN5thPsPIoGG/JsVbJOKNSJoSQiljJpuwq5RnvSK69OmdXW13V0dTPVd0BZeSpHldNm0raV0Bj1PRoSuzAbtwpasSNY6mVEpWyhJ1zjDdEpnKWnGuzVnHtmNw6eQziO2PLUtUZKnjJ1ZrLqP0Hb6wOj5Pjb5DXDGNxAaiUCSGkIg77RN8SiidUVEIFaasgWtUdiq+tAk6o6uZ0p3CK7NLl/tMWk+zApeqjf9475+RVahsMMUXJWOrW2roY0QfqKtKQ95GVb+4qLpIme48k5hWS649fwChPisQqbdTKLDFGh1HX161UyoQQUhEH9VPuWNqG4xbvqpE2gqN2BzTKIKmqW0XtXlO/zBw1jYj/dAlj7MGW2jAURlBNtOW6aXJtzIm6jiZlB3YVaax6Rl6S4/PsjjFDeXrjLRY/dzXXHoYUteUFsuT+R4kaTSnkHJtybNz2zh8KKmVCCKmIld59MUIdWWmXwFJPfaXRikVNqCBXUccUTUxNI9O/NUHWbrSpOgxF4argWJ6WUrb6dWmVkrLlxsZmrp/yGH9mw2Ys3bIqvCIb5uGMH9cu7Krdfh6Gqva8PzKV9CRCY8AaLymFHMirG8+JvZRouQtCpUzIaWeOVxmQ2VhHKVuqqaXkbrTmnWu7tVVPG+4p5d0A7xR1RNEUKZle3lmkfDBzFISlckP96SpgS4Xk+EBP7eOYUvX2ACxPhQLlN4dK3BqK2CV03hqDrso2VnW7U+GVXTcGu/0PR42H/KjnfG+I5QmUVM7SC0qs9IyydnHpfVEXC/+o7IlhzhtlJY+6niY8cXCCOIRL21SOX40JIeQEs675wiX6CGTK0Xz+F6xHcdWEtRR0rzW0dHRNHG7clHtUAdlLt36Yu0S0TBPBPAwTiLXBl7PZkqLEXLFxzRcFr+y0zBRLSJtUX8c2+ixTiDsWrXE4iGNsXFub1YFyp7x+N/vhkNTmdO+cPRYzHjRZASplQgipiAO5xBWooxle2Tcqjkema1KnJJyNkb7KcpVMYtOwJahoLEyXs8gmnXvsKYsm7dZQ0oE02uVl1G/OTb2genSKtfov57HnVSWM8dx3TMmnxqIb7ygw7lJqOqWgA/VIuX5GSanZXPfNftyUGx1fSEQIIaTlMI9Z574aEshzW7HCp6jpFJYq2TrngzZlwy3LVTC5blI5RFzhPCVs2d4chRxUwYYN2Ytb0o+56rBpxoEadmzI6vWHDNJ21ThOnjfmw0/GCsxTw32V64zFs0dGmojrp6uipyjlltxH/GN245Qyjn1HVoRKmRBCKmIdpXzGKSZX/QLZaleDSjmeJqt8l9TDJC2uaosq5ZRtOaYUXYmXsMHHvC4cpaCtsnDVbROeo5RTCjnYbwbBtgT2bb1pFVrEzmz1y2YYX1Go2AcVHaG0lngdaadqnaKiStltl1YpZypoIGmH3heV9spIvjRo6QecDqCaqZQJIaQiVlLKjW10SXUbURpeObllBDOz7HdOluqottjjsSFbaCjNmMd5DTUwaBNXxUpcKZsKup+vuGn8uNk46nXwniZVT+WG7Mwtnpo2lHQXf4T9M0v9j5VDJTb33DEZs8ufMZSxpaAB35Zt2aNdD5P2fHTlnOkxEVLDuXsn9L4ghBDSsopS1nPDNuWoopjBHpxUyBl5+JnmqRLTPtwPs9LkHsdIrCSylPLWOfZszn7e6ijkVjGX2I5d9jbQ5k/Po0JUeu3ZFhbxLfZUm1OW5a2RVdHwNRZde0om5eRlqP/O3yRDKbuqFmc2g/PShoc8PLxzmU/Dxkg925BSv/3PqdVjaI9kRaiUCSGkIlZVysmnu0JpTYVsJIjlNeMtyHv/vPu8v3PzD9omvXcEGGVN8Jk1VwsFdmD11IeTrteubl7GL11lIZ0IEyd83x561LeVNucjT5d5ajrxU2HRcwkVFVXImePXG7MTVnXSmX+HKw/AV89qeamcjduJd5+NtN35EStAw0/ZXAGG9jAMm3FWmhWhUibklHOsHpQ5BayilM9e0BTj3aEiiXLvUDkmthnudu7A9caxMbDFFyU+S34pcto8oYDVUc7RPNvPpf6dAbtip3wjwko2eyXdKeajfdmWPVrchoh5eCQwV2Albe+Ft2nS+yLW+PZ/JMEJDw275ly72hDXxtz1iW2bN708Jni2ZHtQxfarXC8id4Xo5oV5V9e5UCkTQkhFrKqUi+yLB/YV7HDekeCqNksCd2ohp4glrtXIcqhqDSXnqF2v32L24txrcVXStq9u2ziDP/vwnuJS7duem3h9G7dhj/YehHRXQuZBP1E42Nw7QaCvXcHueq944b3A3NWIZ8NtsuoHu/ZndeMO1W5Qbbcq21hFTvkd4JRt3R+jfbtwGyestlN9sjZUyjFOsq2tlpueg/eFPjQz1mfWm++cj//OMc5r/a7UOcyjcFImhJCKWMV8cc99wi/snkUVzXCHNusRWJ97cU3H98gtepNYemVk4ZJtGgqoNcs8YZ+fID8S5qBQd+7NE80Hx9Sx3xBszC39xEfuMtsyQ00YR8mldPO31+/ZJqFuU8rIM1KuONe8j+ie75+LmzG88ND4N0164fAYyetOfZckcq6LU5ecplImhJCKWEUp33uhMfcH1En2Rlq3aRC5yyXuzCm1FFYSXi6JQgJZlGzmoNAtZ4qInVEw2P04dDlT77yfR04c67ypjLs0Q+Ueyz9Xabn9aPV3KG4X7m1sJepg5L9LFHfXDBfkhE7YpCslON5Hrt7mHNNrQaVMCCEVsY5N+cKETSumgCz7l2XjisQx4xpqqSTPfSInekAl2UrKOTbyDMYpZJQ9f0SbJ22SXbqcCiVUbwwzzYiGzFRtyX6NnUuVkVEvv5CMPEoJjs08F7NReymleYWyMFdNmWELQ6VMCCEVsZJNufngqV//Dm4q5ER4+IU/w8OUndMlR5Eld9FDjutGHDPcKquA6LXktlPGqsV6hDe2Ktqdz/AKSZERP5lnkceLBNNY/RqtR2lfL2QrTarX1LUi8rVKrjAyKuikTX5HU/sQ/Tjt+E6kXxoqZUIIqYh1lPL5u7/7l5s0x52v4z5ud5eyFLGVR8ynuMSGlCJbGRvHQO/xTidOoow+lqJJ2stCysGJk7868fP0X/Rjl5uktB1msN0mwzPyyHLQGWsbnUEhxz04hn+9uNbYDfRDjqqendiqL3c8D16uNG/1cqBSJoSQiljnhUTn7f6KoXYHSstQwlYaU+U5+QaPU4xQXCnFDPT8MBO25S7+GAeBErtwoYJYrc0TNvYue0O9BcNK7ZolijkjTTLtGEba3mMripYuTkohD/pN7XOhY7esCMnnGGIraGcuSa3Cg+WtAJUyIYRUxCpK+ejc5kPKTpwTJycPK69MsmySKcUQsbUllXJOfRKkbMtZbe4q5JStuX8uUY+sdz9kqbJ9e+asSlLlj1LGtZCoX3RVYKhbr/3EiB/I1y/Pcr2xahzAU8YyzNobf72A7qfAmnPunOKcH8RZESrltaj9C01OL9W9L/V0s5JSbn1Xw3eqUfbNlILuM4MtOWlLy7F/WmlcpWfVY4zjtNEuIaWc7WUR8b7IrWpWO+Xu9MfasTvnbEJU4BmwOAXjyRs2GyeOZZsPdXSpUrbqB9ht1inkRrG7NuR2yumnd+Yb9+fGNGRTplI+uVS/9K0EttP6sM3rYl2bcvtjjO2dqfsPAfUlwfCo50bmzr/3Dt/YoEzZeXNsyCl7XCyPRD268O6a3IsLJEkq30aFJOz4ofyzFXKXwP+c9Pne6DBNzIZs2UpzVFzuZJVQcyHWsBgUXaOlgL0+cJQz4L+grYnTBacUew7edzv8JdaADXk/bltV3VTgaBAc9L44fU/0RZa/pSwyyKcoiMxNvMXr4WY1Y5vPOmBjN6RMrE29UZwG9TjD2Kzs/fB73An5GEHzBSGEVMQq5gvdNrerdiOmCQ85eKuz1HCN8ikzR5/JqjC2CeXGSZgkgnmUuBkZeSVJWzN6bdy2pQ4iZ7nCeeU6lTUix9or1U5dHGuJDexlh7UxlVLoGW+5t81mtp3Ma40FVh1+uGGSiMTp3rkUMFvYZYXb2EIiA6kbk66ZzrWNuL+k7QzlQdr2Utqx0ajqziw2+LX09ZcCVMqEEFIR6yjlM0PlBVeZ9e6sripz76HZqg62YMlWJQUbgMlHSiPuR3ZehjrJreMgU/vYbafsPujSlfgyOZnEski1SxvHct8KpXEFvKsEo8R3iL23xxaMuyVsn8mNPkcNx+L44yxfDee1bV5a9QaF068ynBeG7rbDzUD3ty32Krs/gNY3SlMpE0JIRayilLGxFIYvxTx16yZ1nMbdhDGbqV/+DOQoPiBhW3bv9um0s0orpz9c+/0+Xjh42IEpA6JR7wylHLQLi6GgnTRitXWs/BTWD+06dtC9IgvUr23TJU2XxupgH2C32xSF3KVov7MTxmyxbTdi65bmpDrhodUuX0hECCGnnHVsys7UH/zpJi+Re2yoOSt+L8y9y1oqexZy7KBu3BxFGIoXKmcUAfubf7okK5sS6WFc0/ABHce+jJlspDkV8rdEhmV2CtHPyfitXptY9TLzGOVTbHk7dOfz+9NUu25w0ffRGiTpvDrzs5umv/9zAId1KmVCCKmIdWzKuXbXGJYpMqaYDdt1epd8wt1xjM00c5c8RLKqBcrUy8v1Ec3LpSRyHrFr6Nv/SuqZUoA5aQuZdSN/xBB1y+/aK8d7xlT2MyjJOdrFXOXZ9fP6I7TyOYBRmUqZEEIqYiWlnDiddbN11LbnK5tMYj/h5z0VGHHlKKXEx3jEXdkxX/pkPlkXzdwLz8+imBGm3S5p61jRP+/L/xGVyswi5bc8ZzXmuIygq1J70tmHseznB2LUqqOSuqegUiaEkIpYRyk7FPlJWnm0H6L24KGnhu8DHfd5tvOdiVzVmhEtVTuN2Q/HvHy8f3oG+2ZWWakdddcvOJiH4xGxlOIM5X0cXlxsXEsXXLLSOg7XG+LAr5ajUiaEkIpYRylbO7dFTxIZxPxN3Ye33Ceukk+w+ao7+ya69M3W+NmnDqei7tuzhkkt/1FniZF6MqwAyVG1LblteZRfoVHKLztz1y8+pyLHhWOqfkOYLs6HfaSPSpkQQiqCkzIhhFTEuht93mq4fThhxBrWdAGLxDHWzPYLi3qPW+a44PXTlniijXB8Ty57PVc4J12vQt6Lwd1NwTnc6HI3fTI26XKLTOa3Sx09HMWY/lyDWkwlU+oxpQkt70jX1DeIwsesCSHkVLPSRt8IxWVh7UmFXnCe++CAu2MV2H2yX9bjPjAQdq8L3XG9zbdc16qcuKYqgBFhHyn69C0yXRq9RKmnXCKkNs5iG7Rjh9ySirk0ztKsvSk9hdw5JPbkePJRcTXPrAGVMiGEVMQqSlk8tai9/3u25UFodzKed+TOaZ5LCePgi0msCg2NzZZ9eqCgnfYYpZhTyt2wlwXTe9fmKGbL9p/zc1XRgvd9n+MaF7wUhX/NoeJK1eASEunQqnhGxXqw5yvclWlqb2lURSPz0QpQKRNCSEWs+vCIp7wcxbyLY8hF13PCukOOEAOmog4Gtw+cWAZX50PBDzhKwcMPbl7553MkaXk1xqdtVbn0jganygmtKBKeG1nNNFUeHtr7omERlbu2oDQWb+brZ4eB5QWtCJUyIYRUxLreF5atsn/nWuK9P2OlQaQuYinmVnl5duqBW8gu6GioDnUzDI9sFwcKzPUDzvCEmcVbJtPTxYnf319Q55Nngx+Buc2Q415d+hrUSpTxohxIIZ9kqJQJIaQi1v05qJXvcrneF1NerONJ4u5Sh0pajwJpt03c9lwb11G/clSgmLvDhKoLZdWGWV4W7nEXP1IPr5yht8o+D/9py/2PEjg2+J7tfbjWSLdPrtoOitxSiX7gV0Ba5Nld8xjleLPA1kmer/zE8ytBpUwIIRWx6s9BqavEcu5MXZwZvC5aU6mZR/pu66lv6xauQ+XXT+epZ+fW6L6APfAqDk9d+0+1Ja5FAp+dftn3l1MvL49AWd5qZBjHf8IxlIfhweKa549a9euocMBvH8fjJUlsRZFLTlm5ec6h5rz6ZFTQsrmXLAZy9ypi17ik2h3xYxtLQKVMCCEVsYpS1q3zFFeOYjYUsvuOC0+5xjwmPKXnq9ho/GCccNrW/jnYvO9+fqo53jg20lbFyVDdSXvr7NmltQ3r1KFjK00pmMC1WSsZ2Vg2ZL//8u34YVu8p4L7YSEhowA2wzYY2O+1F6+fxxg7asFYG5RVVIaTRUk1c1dxE7Lc5+0e55SVUsj5SnURx5ZK9gColFfCfNiEEEJ6rPPuizPebbU5EYhrqGc73M3Avtu56sxNG1R5I+dSP5n2fJhbb4LmeDM04qpjG22VdBcPvmeCq5i7UnM2pY1VSHK1AjdevlKO4niwdMHutRrhsg3IbcOTYw41m02/LCOP5Jgc5e4wgYQi3ndFoKyRdvvgCtVdvSXyGFWPELQpV8acHVLHyqh+2E5kTo7heOKkTAghFbGK+WJz5uzuQ+aSbRcY3oSzNt+KHqXOzCNYrZGbAUObsvb+D8SxzBxHm15cBOO6xzGvtY6EqSjbdAQsulkSylq1Z+Zw26IXto/vbIh6O2s5FcmIU0zcTDHNjDEjORt8ZhtGdmz7p/tBKXNmJG1+veqCSpkQQipiFaV85ryz0fPRF9WPSJPLGFe4zYxK8CixKbdXys3fnk9cUlV7L0rK8e+Lr05m2cRLMaJ53Wsd/J6A0Q7uSsOMn8EcY3GfmZt3+8Htk3SZk14fYGCN0WHg2N1x/5qkG5Ph8FoejZ4TKmVCCKmIVZTy+efeMzrtrCoku8x0nCUVc0tMlbSfLbXtquwSfKU8/uEat15jyL2GUDzTbu9EsMqYUu+WknZK2UpDfTHnatJ7HH5C23fncsuO1MdavVnfw0PMG3NBpUwIIRWxjlI+597Z8ipRqJYCnSPvJclVzgBwBFcpD49Tdus+lqJbY1UQolSlRtVaYbvkKMTSVUiJekt5HYX6xF7RTF/plGD1W6q9gquAxHW3aTaGHl/qGpeESpkQQipiFaV84Tl3r1HMicFVv9G4nuILpy1dNQDLrBxS9ShVn7nXZam30pVGLK85yVHG+7jDc65qLMkrlbdFTr/l9lXOtVqKuSSvLk5lDsxUyoQQUhGrKOULzhnvfTGF2u6ALSVKeJCuQEWm4obO5yqoOe36OddUqvKzVBtSNmY/j9xVSSpdn1RblihATzEnFHSqjLmZ8qbE3HbIUsYTVPUaUCkTQkhFrON9sd15X6xxJxqjjpes1xhbbpd2hseVppTfUlv7zJnGDY8paC/ujPZ7ixJVl1KNse/GWPU450qnbfuS73CqflOU86GgUibklFPbpHTaWcf74sx074vcu2dsgG3kyDwHANuFB+fZYs+C5e6ZffVS+qVMtePcTGkHU80aq5A2fqjMElWdE55DkVI2viO+gvb7bw1Pm1Q/ruUhlBq/h75JUSmvROmETMhaHHoSIkPW8b7YDr0vpnhFWHe52MDamgqiPK85yba1LaiYc8hRxrV/sW0VO2zbsxnqt01jqW33BjxH/8X6wF3h2Yr5yDlO95n13XHbKURKKefsmSwpZpZeGY+FSplMpvYJmZDjBCdlQgipiFXMF/fZ3hUMH6OwrCXHBvmbF1vkLeO2K29onS1Y5s7pdtWS2x857ZJroqqljd2ldCie2+ZnHU3jbwSm+9NanqeW1qHx3p3LHPep+CXExmNpO8VMFkuM+9qgUiaEkIpYRSnf/8zXJudhKSpXkcWUl7vR4W5i5GxoWRsfc5CzedKyxOZf0mWwcMN0l2Z+96NUnjFc1dZiKbBQ/PTmYPz8FHLGn9UfU9JaxMZhbrt44ZE87c3Vk6MvT86VEELICWAVpXy/bb5SXkOtuerMUl5j7J0l7n5THqNeQhmkrte6tujqxGzbHHesvPYvUXepFYalpIGAq5tnC82zV89BbJzljtuYXbqUkP3cbY+U3T68Konb8ffxwuFjVikHd0E9aOmEEEIGrKKUH7C9szhNDaptjO1yyiPIc96hXUUx57Xk5FXiJbOLP94uvc+73D6dUrFx5ewqwYQKn8FzoOSBh1wlPGZFmLNSS3lV+F4ZAbWdyMNq8zF26Zy0a0ClTAghFbGKUr5ohFJ2yVWgS6i4cB7r+tfmYt3lp9TX9G8dYdc09wSy+m3cfkMOJbbHro2dJDl+ybnlzuHlM6fN2CPSXKZ6zfAF36X3M/cVs5U23+c5tgrq8juAaqZSJoSQilhJKf/pGsUAGGdXHKMil/RXnsSMG/0ppRW3A497+XhspTM2zykE/ZeNNp7TH3kO1hyjOdde6iMeyzfXnh+r15reMiVQKRNCSEWs46e8mf5E3xpUq34XYIqqzLLbZ75wPRU/L89k0vEU5H1U2fBZtF0axlyzp14lcT5YruvJkf9EZiqvQ0OlTAghFbGSUo7/HFRtd6rTSql6nqJu92WWpy1REts1h9ZpHMaRaz5rdn3ix1gzij0r8caOKfg5lPqSUCkTQkhFrPM+ZTkbDO/uSMfYbpfDoe/iIlobAAAKvUlEQVS8fWbxfy24nClq11K524z2rEFtbBNq7qRzdoE9mpSK7sqMKvjpSn1Jahi7J5qaJuTTAgc1Oc5w/BJCSEWsYr44P7GMK1nmhA0hcQ65HDmnubZDK+Y53f3m3GibaorINRFsFmj/nLqfZs6ZsXlKTSFHsfhNvVJmjENBpUwIIRWxilI+T9xHIod3qHNCLyAx7mLu3bdEBS+x8WCW5RQ15vHvOVjCJWyMQsy9++co35TqHVe/8Q11kjb0NjPotKMF1qbnGOFj1G47/1hKfs15IgSVMiGEVMQqSvkc60UkkTtSexdL3bWitiOHzQo3wFYjrPrQQgZL3n2nKMUShZqrgIvyLKj7VBW5hF27RjbYLl5G+70vcc9sFXxqHG0ObGumUiaEkIpYRSlbaiR0x3JtRK66cJVxjno6291V51MqlkJfXiPMzxoKboqnQmn9llK/5fWoS/PMYS+ekym2Z/d7dlbTebkK3voOr/Eipxh19RIhhJxyVlHKRXfo5E/+uD8pk7b/pBROiV3aqsdpZimlXWqrHqMEx9R9DgVcm2o9BCVtkFLVVp/EFHRO34+ZG6bCkUEIIRWxklKO24WHcQvvEzP8gOkcqrfWp4NaavWlXUIxzqHc57QHUxWPZ6rPc6gfc+zPLYfwmOFoIYSQilhFKbtMufu4KrsWFXLoHdvjRm0+u0t6SizxhFuIJb4La9V9DUoU8iGpY0YjhBAC4EBKeQq1KSxyMjguKirG2ROkak8z6ylluW/6b04c5lVffWrNq/+vhvqchrxqq8/UvA6A6MJeAyITfsueEEIOz2dU9bK1Clt8UiaEEJIPN/oIIaQiOCkTQkhFcFImhJCK4KRMCCEVwUmZEEIqgpMyIYRUBCdlQgipCE7KhBBSEZyUCSGkIjgpE0JIRXBSJoSQiuCkTAghFcFJmRBCKoKTMiGEVAQnZUIIqQhOyoQQUhGclAkhpCI4KRNCSEVwUiaEkIrgpEwIIRXBSZkQQiqCkzIhhFQEJ2VCCKkITsqEEFIRnJQJIaQiOCkTQkhFcFImhJCK4KRMCCEVwUmZEEIqgpMyIYRUBCdlQgipiDM5kS6Rh+jduHsfIALpR5Duv0EcJONYYb0D53Q4LDc/QGPpQ+UH8/DDYmXuwmSXLhkvUFYs/kxpvPSl1+eEZ11rxnkdU4/oeZ18bd7VWUNK+h+9FvG+Is3ojOS7z0OMvN2vgHeuf1rUrWYXX5zrDcXpx/PLcdJI4jzcuqh3jcn4kbL6cd20w3yd+AL0Z7twGvf8PuSDH77rnap6JQrImpTvxt14yuavQTZtD2x2n6UR2hvZTdSbzb6X2s9NGpFNF68J2Mdpj9v4si/HS9PLEyJQkb3eb8+34cDunOzz6OL3j6V3vGnzRa9cDMpp03QTRhcfXXu0x1acNv/9+bbswPlBGIZpBufdeoXi+Hl4eWIfv58+liZWL++4V494vTR6LV4cuHnooIw2fiiOWGnQntcuDwmkEdHeMB2mkSasOxZg06XxwwBg00xKm+68Ds+LNnG0i7NpPm96eWxC59E/PgqGbQd5HHV5bN1jUWywO97iqMtjK0ddnls52uXZlXHU5bPLc59+24Rt0OZx1JWzi7dP05bdlrHFUZeuq2svz23THm0998f9a8O+ngC2Amwh3RSzFcEGgm0zQDaQJmzTO950x9uHfuISFELzBSGEVAQnZUIIqQhOyoQQUhGclAkhpCI4KRNCSEVwUiaEkIrgpEwIIRXBSZkQQiqCkzIhhFQEJ2VCCKkITsqEEFIRnJQJIaQiOCkTQkhFcFImhJCKEFX/Pa9eJJF3ACh+BR0hhJxybi19n3LWpEwIIWQdaL4ghJCK4KRMCCEVceomZRH5ARH5qIgcici3ReJdKSIfF5FPisiL16zjFETkYhF5t4h8ovn7dUa8G0XkIyLyIRH5/bXrWUKqL2THK5rzHxaRJx2inmPJuL6nisiXm776kIj81CHqOQYR+RUR+aKIXG+cP7Z9l3Ft4/pNVU/VPwCPB/BYAO8F8G1GnC2ATwH4BgDnArgOwDceuu6Z1/dzAF7cfH4xgJcZ8W4EcMmh65txPcm+AHAVgLdj93N43wHgfx263jNf31MBvOXQdR15fd8D4EkArjfOH+e+S13bqH47dUpZVf9AVT+eiPZkAJ9U1U+r6t0AfgPAM5av3Sw8A8Brms+vAfA3DliXOcjpi2cA+DXdcS2Ai0TkoWtXdCTHeawlUdXfA3BbJMqx7buMaxvFqZuUM3kYgM/2jm9qwo4DD1bVzwFA8/dBRjwF8C4R+aCI/MhqtSsnpy+Oc3/l1v0vich1IvJ2Efmmdaq2Cse573Io7rczS9foEIjIewA8JHDqJ1X1zTlZBMKq8R2MXV9BNt+lqjeLyIMAvFtEPtbc+Wsjpy+q7q8EOXX/PwAepapfFZGrAPw2gMsXr9k6HOe+SzGq307kpKyqV0zM4iYAj+gdPxzAzRPznI3Y9YnIF0Tkoar6uWYZ+EUjj5ubv18Ukf+M3TK6xkk5py+q7q8Eybqr6h29z28TkVeKyCWqeutKdVyS49x3Ucb2G80XYf43gMtF5OtF5FwAzwJwzYHrlMs1AJ7bfH4uAG9lICL3EZH7tZ8BfC+A4A5yBeT0xTUAfqjZyf8OAF9uTTjHgOT1ichDRESaz0/G7nv7pdVrugzHue+ijO63Q+9gHmDH9G9id3e+C8AXALyzCb8UwNt68a4C8IfY7Yz/5KHrXXB9fw7A7wD4RPP3Yvf6sNvpv67599Hary/UFwBeAOAFzWcB8O+a8x+B4VVT67+M6/sHTT9dB+BaAN956DoXXNsbAHwOwD3N9+7qk9J3Gdc2qt/4mDUhhFQEzReEEFIRnJQJIaQiOCkTQkhFcFImhJCK4KRMCCEOqZcNOXEfKSL/TUT+b/NSpaumlM1JmZwIROQcEfngoetBTgyvBpD7iyH/BMAbVfUvYudn/sopBXNSJieF7wbwPw9dCXIy0MDLhkTk0SLyjuZ9Me8Tkce10QHcv/n8AEx8IpGTMqkaEblMRD4mIq8SketF5NdF5AoR+R/NO6Of3ES9EsDbm6cV39q8BOZ6Eflbh6w/OVH8BwA/qqrfCuBF2CvinwHwHBG5CcDbAPzolEI4KZPjwGMA/BsATwDwOAA/iJ0yfhGAlzRx/ip278i+EsDNqvotqvrNAN6xem3JiUNE7gvgOwH8poh8CMAvAmhfMfpsAK9W1Ydj93Tma0Vk9Nx6Il9IRE4cN6jqRwBARD4K4HdUVUXkIwAuE5FLAdymqnc2YS8XkZdh94Lx9x2w3uTksAFwu6o+MXDuajT2Z1V9v4icD+ASGC8DyymIkNq5q/f5qHd8hJ2weDqAdwKAqv4hgG/F7j0KLz1OP51E6kV3b3y7QUR+AOh+xupbmtN/BOBpTfjjAZwP4JaxZXFSJieBK7H7SSE0qvlOVX0dgJdj93M9hBQhIm8A8H4AjxWRm0TkagB/G8DVItK+yKv9hZgXAvjhJvwNAJ6nE14qRPMFOe5sAVyuqh9rjv8CgJ8XkSPs3t719w9WM3JsUdVnG6c8NzlV/X8AvmuusvmWOHKsEZHvBvAcVX3BoetCyBxwUiaEkIqgTZkQQiqCkzIhhFQEJ2VCCKkITsqEEFIRnJQJIaQiOCkTQkhFcFImhJCK+P8ADdrofEN5UgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "iris.quickplot.pcolormesh(ilap_cube[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum absolute difference: 96.0\n",
      "Mean values: 3684024.8 3684025.0\n",
      "Standard deviations: 92327940.0 92327944.0\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "diff = psi.data - ilap_cube.data\n",
    "print(\"Maximum absolute difference:\", np.max(np.abs(diff)))\n",
    "print(\"Mean values:\", np.mean(psi.data), np.mean(ilap_cube.data))\n",
    "print(\"Standard deviations:\", np.std(psi.data), np.std(ilap_cube.data))"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}