jmeyers314/Hankel transform demo.ipynb

## Hankel transform demo.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:32.009471Z",
     "start_time": "2018-04-11T17:46:31.176698Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import quad\n",
    "from scipy.special import j0\n",
    "from scipy.interpolate import interp1d\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:32.037047Z",
     "start_time": "2018-04-11T17:46:32.011502Z"
    }
   },
   "outputs": [],
   "source": [
    "# The optical transfer function of a Kolmogorov profile \n",
    "# (with lam/r0=1)\n",
    "\n",
    "def kolmogorovKValue(k):\n",
    "    return np.exp(-k**(5./3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:32.070972Z",
     "start_time": "2018-04-11T17:46:32.039267Z"
    }
   },
   "outputs": [],
   "source": [
    "# This function shows how to turn the Kolmogorov OTF into an image\n",
    "# using 2D Fourier transforms.\n",
    "\n",
    "maxk = 12\n",
    "nk = 1024\n",
    "subregion = 32 # Return only innermost (subregion x subregion) square\n",
    "\n",
    "def makeImage1():\n",
    "    kx = np.linspace(-maxk, maxk, nk, endpoint=False)\n",
    "    kx, ky = np.meshgrid(kx, kx)\n",
    "    kVals = kolmogorovKValue(np.hypot(kx, ky))\n",
    "    xVals = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(kVals)))\n",
    "    slc = slice((nk-subregion)//2, (nk+subregion)//2)\n",
    "    xVals = xVals[slc, slc]\n",
    "    return xVals.real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:32.592289Z",
     "start_time": "2018-04-11T17:46:32.072757Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1052e08d0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAD8CAYAAADnhGhBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+sJWWd5/H35/7oH4AIDEhYwAV322yIGZH00iTjbnCNTkMm0zqTITBZQcdMyy4YTdyMDJOs7hgT1hl1dZdA2rEzkDgiWWXsuL3DMMSNa7Joo0Hkp9MhrXSnpUHkR0P/uD+++8ep3j1973m+t27dc849p+7nlVT6nKp66jynTvVzq+r51vdRRGBm1mYTq10BM7NBc0NnZq3nhs7MWs8NnZm1nhs6M2s9N3Rm1npu6MxsoCRtlfS0pL2SbumxXJK+XC1/VNJldcpK+qikpyQ9LulzWR2m+vd1zMxOJmkSuB14D7Af2CNpV0Q80bXaVcCmatoC3AFsycpKehewDXh7RByT9KasHj6jM7NBuhzYGxHPRMRx4B46DVS3bcDd0fEQcIak85Yo+++A2yLiGEBEHMoqMdQzunVaHxs4dZgfOdo09IJjrOETPGvswZ+jvMbxOLaiA+S333Vq/OrFuVrr/ujRY/dHxNZklfOBZ7ve76dz1rbUOucvUfatwL+S9FngKPAfImJPqRJDbeg2cCpb9O7hfJgG0BiovyfAmmhYxyb1aPpZ/TbftMGab1iszy1dw3rk2+xfHX8QD654G796cY4f3v/mWutOnveP/0LSw12zdkTEjhVXYmlTwFnAFcC/BO6V9JYoPNPqe3RmdpIA5qndoL8QEZuT5QeAC7veX1DNq7POdFJ2P/CtqmH7oaR54Gzg+V6V8D06MztJEMzEXK2phj3AJkkXS1oHXAvsWrDOLuD6qvf1CuDliDi4RNm/Bd4FIOmtwDrghVIlfEZnZoss44wuFRGzkm4G7gcmgZ0R8bikG6vldwK7gauBvcDrwIeystWmdwI7JT0GHAduKF22ghs6M1sgCOb6eN8wInbTacy6593Z9TqAm+qWreYfB/5t3Tq4oTOzReZb1l29ZEMnaQPwPWB9tf5/j4hPSToL+AZwEbAPuCYifj24qvas3AC2ufzblmnvadOe2mSbavK9J0bkdqwa9p5G8p2TnlVNljZYrkfaU5v9nk17ZEu/5yolxQ1grmUNXZ2j/xjwbyLi7cClwNbqhuEtwIMRsQl4sHpvZi0wT9SaxsWSDV0VrXy4ejtdTUEnQvmuav5dwPsGUkMzG6oAZiJqTeOi1vWMpElJjwCHgAci4gfAuVUXMMAvgXMHVEczG6IgmKs5jYtanRERMQdcKukM4D5Jb1uwPCT1/NaStgPbATZwygqra2YDFzA3Pm1YLcu6Qx0RLwHfBbYCz1UP3lL92/Oh2ojYERGbI2LzNOtXWl8zG7DOkxH1pnGxZEMn6ZzqTA5JG+mkTHmKToTyDdVqNwDfHlQlzWyYxFzNaVzUuXQ9D7iryg01AdwbEd+R9H/oPEj7YeDnwDUDq2W/w0iSEIFGD9pnIQdNw0SycJAG5RqFpECeDKDBA/PpOMLz5XMEZeWSkJXi581nx0DTBAJ9Dj3JfrMBdgR0OiPGpxGrY8mGLiIeBd7RY/6vgCGlIjGzYenE0a2xhs7M1p75tXZGZ2Zri8/ozKz1AjHXsgxubujMbBFfug7KqPSsNuhBTXs0J0tPlZP2rGoy63VNlpXKDSC5QJNeV2W9j3PJsuwh/KScSj25c+WkkY0TCCS9tUPrke1DZ2wgjkdy3I6h0WnozGwkdAKGfelqZi3nzggza7UIMRc+ozOzlpv3GZ2ZtVmnM6JdTUO7vo2ZrZg7I0bNEENIAFQKFUnCS4plAKaT3T9RLpeGnkwVtpkmAsj2VcMHy0thGFmZ2dnyR2UhJJPJ+KIzvbeZRWFkCQSC5LMaJgoohp40HYOiD+YcR2dmbeYnI8xsTZh3r6uZtVnnoX43dGbWYoGY8SNgZtZmEThg2MzaTg4YXrEmWUoK4SBDDSGBYiaSRuEegKanG5XLwlJiqlD/rI5JFpXIQmeajP+QhYnMlkM3VAgTAfKwlEL9NVv+XmkoS7kWfQ896XvGk5qC/p7RSdoKfAmYBP4qIm5bsFzV8quB14EPRsSPs7KSPg38MfB8tZlbI2J3qQ4+ozOzRfrVGVENqnU7ndED9wN7JO2KiCe6VrsK2FRNW4A7gC01yn4xIv6yTj3adSFuZisWiPmoN9VwObA3Ip6JiOPAPcC2BetsA+6OjoeAM6qxouuUrcUNnZmdpDPc4VStCThb0sNd0/YFmzsfeLbr/f5qXp11lir7UUmPStop6czsO/nS1cwWWNbg1C9ExOZB1qbgDuAzdNrlzwCfB/6otLIbOjM7SdDXJyMOABd2vb+gmldnnelS2Yh47sRMSV8BvpNVwpeuZrbIXHVWt9RUwx5gk6SLJa0DrgV2LVhnF3C9Oq4AXo6Ig1nZ6h7eCe8HHssqseQZnaQLgbuBc+k09jsi4kvL7d5dUtPBW0r6HEICoFLIRxZCsn5d+bOyMJF15dCT2JAsm+5d//nCfCDdV9n9ZmUpQArZSyZmshCSZNnRmfKy48mywu8Zx46XyygJVykuWUHoSbKo/GGDO0eJUN/O6CJiVtLNwP10QkR2RsTjkm6slt8J7KYTWrKXTnjJh7Ky1aY/J+lSOj/JPuAjWT3qXLrOAp+IiB9LegPwI0kPVMtqd++a2XjodEb07xGw6gRo94J5d3a9DuCmumWr+R9YTh2WbOiqU8iD1etXJT3J4l4TM2uN9o0ZsaxvI+ki4B3AD6pZtbt3zWw8dDoj+hZHNxJqN3SSTgO+CXw8Il6h0737FuBSOmd8ny+U234ixmaGY32ospkN2hwTtaZxUSu8RNI0nUbuaxHxLajfvRsRO4AdAKfrrD6MI25mg3TiyYg2qdPrKuCrwJMR8YWu+edV9++gRveumY2PtTg4zm8BHwB+KumRat6twHXL6d5diWKWkmxwnCxLSjqYzfIzkWQhJLEhW7a+vGxjOYRkbkP5Z5vbWAovKX/nbFnT7CUTM6XwknKZySPlOIvJdUkIz5Fk2dHet0vyUJAyJVlDopSxhfx4jGJ2nuyzBndxFAEzWeaUMVSn1/X79D4umsfMmdnI6ly6rrGGzszWnmU86zoW3NCZ2UlOhJe0iRs6M1vAl65mtgZ4zAgza7VOr6uHO1yZJlkXSmWywXGyLCRZ9pIGmUjSEJJTyiEkc6eWl82eWq7H7Gnl+s+c0ntfza7PwkuKi8iuYJSMzzJRSCgydawcFjH9evnDpg4ny6bKy7JooZI09CQJL1ES8pEGg5TCdNIQjyYpT+pZkwHDZrb2+NLVzFrNva5mtia419XMWi1CzLqhM7O286XrgBQf3M/KZA/uTyR/kZKxGjSddEEWymUP52c9qzOnlz/r+OnlntVjp5e/28wbeu+T2Y3FIsyXq0hMlPsLNV/e/xOF1INTR8plpl8tL1ufJR6YXP6xM5n0GGsu6VmdS8brSJIcFHtWKR/HkfTVFv+/JN+rLt+jM7M1wQ2dmbWa4+jMbE1wHJ2ZtVoEzK61xJtmtvb40tXMWs336FZLqSs9CSFJx36YaPZQf6zrHVqQje+QPZyfhZAcPbNc/2NnlQ/C42/sHZIw+4Zk/IENyQPik0nIxFwynsTR3t9t6tXy95rbkIWQNLuU0lxh7IrZZH/MzJa3N5vsq2SZJsvLYq7w3bJwlcE90199tBs6M2u5tnVGtOuOo5mtWETnHl2dqQ5JWyU9LWmvpFt6LJekL1fLH5V02TLKfkJSSDo7q4PP6MxsATHXp15XSZPA7cB7gP3AHkm7IuKJrtWuAjZV0xbgDmDLUmUlXQi8F/jFUvXwGZ2ZLRKhWlMNlwN7I+KZiDgO3ANsW7DONuDu6HgIOEPSeTXKfhH4E5bIawpu6MxsgRPPuta8dD1b0sNd0/YFmzsfeLbr/f5qXp11imUlbQMORMRP6nwnX7qa2cki7/Bd4IWI2DzA2iwi6RTgVjqXrbUs2dBV18F3A+fSaex3RMSXJJ0FfAO4CNgHXBMRv843lmRdSMaSKGYpybKXZNvLQhWSzCaxoXcYydyGZuM7ZFlIshCSo2eXQyPmfqP3YA2nnXGkWOaNG48Wl62fKodaHJstf++Xj2zoOf/wS+U0KvNZ5pjk4kOl8Axgcqb3/p88lmSwOV6uh7LQk5nCQBlAzGaDb/T+rbPsPFE8vvvTW9rHXtcDwIVd7y+o5tVZZ7ow/58BFwM/qfbRBcCPJV0eEb/sVYk6l66zwCci4hLgCuAmSZcAtwAPRsQm4MHqvZmNuag6I+pMNewBNkm6WNI64Fpg14J1dgHXV72vVwAvR8TBUtmI+GlEvCkiLoqIi+hc0l5WauSgxhld9YEHq9evSnqSznXyNuDKarW7gP8FfLLONzez0baMS9clthOzkm4G7gcmgZ0R8bikG6vldwK7gauBvcDrwIeysk3qsax7dJIuAt4B/AA4t2oEAX5J59LWzFqgn09GRMRuOo1Z97w7u14HcFPdsj3WuWipOtRu6CSdBnwT+HhEvNJ9/yAiQlLPvwFVL8x2gA2cUvfjzGyVRLTvEbBaF9mSpuk0cl+LiG9Vs5+rYl2o/j3Uq2xE7IiIzRGxeVpJzm4zGxn9fDJiFCzZ0Klz6vZV4MmI+ELXol3ADdXrG4Bv9796ZrYaIupN46LOpetvAR8AfirpkWrercBtwL2SPgz8HLhmMFVMZAPgZCEkWYaSqXI4SEz3Xja3sVxm5pTlD2QD5SwkUA4hATjnTa/0nP/WM58vlnnzxheLy06bLIxyAxyeK5+h/+LIWT3n/2zdOcUyz3N6cdnxmXXFZZNHy/tx5kjv/T/1evk3mzySHAPJ8aHkuGK2HJZSPI7nBpyipCAQ82st8WZEfJ9ycM67+1sdMxsFY3SyVoufjDCzk7WwM8INnZkt1rJTOjd0ZraIz+jMrNUCmJ93Q2dmbRaAz+gGpJTVBIrd71l2hyx7SZr1JAlLmS+El8xPl7c3uz5ZVk7kkQ5mk2UiKYWRbH7jvmKZS9YvTCbx/50+Uc5s8sp87wwlAE+sW5hybGlHkqwhr72WZIh5pRzyUdr/2W9W+p0BJrKwpabZdArlIgufGnAQ2zjFyNUxOg2dmY0ON3Rm1m6106SPDTd0ZraYz+jMrNUCwr2uZtZ+bujGX4Me3qxc3oNX3tx8krUqNpQf6M7GeCg9oJ/1rF6yrjzUxxkT5UPkpfly72/JoY3lB/f3beydCADg8IZyLsP59eVe0tL+z36zfh8fY8mXrmbWem7ozKzVHDBsZmuBA4bNrP3c62pmbdd7qKvx5YbOzE4WuDNipKRhANkD1uVlkS5rUCZ7LnsiOZomy8vWT5XHHyiN8ZA9nJ+FkJw2UX5wH8rbLH1eNgZF9r2y/ZHtx9L+b/I7L1UuTTLR9FhdFXJnhJmtAT6jM7PWK2cJG0tu6MzsZC2Mo2vX4I1m1heKelOtbUlbJT0taa+kW3osl6QvV8sflXTZUmUlfaZa9xFJfy/pn2R1cENnZotFzWkJkiaB24GrgEuA6yRdsmC1q4BN1bQduKNG2b+IiN+MiEuB7wD/MauHGzozG6TLgb0R8UxEHAfuAbYtWGcbcHd0PAScIem8rGxEvNJV/lSWaHaXvEcnaSfwO8ChiHhbNe/TwB8DJwYouDUidi+1rb6bT75btix5vkXpsgZlkpu6yqLP58rLjs2Wf7bDc71TomTjO+RZSMohJC/Nl8NBXpl/Q8/5pfpB/r2y/ZHtRxWOgya/81Ll0uemmh6rq6SPAcPnA892vd8PbKmxzvlLlZX0WeB64GXgXVkl6pzR/TWwtcf8L0bEpdU0/EbOzAYj6DwCVmeCsyU93DVtH1o1I/4sIi4EvgbcnK275BldRHxP0kX9qZqZjYX6Z3QvRMTmZPkB4MKu9xdU8+qsM12jLHQaut3Ap0qVWMk9uo9WvR47JZ1ZWknS9hOt/UyUo+LNbHT0sdd1D7BJ0sWS1gHXArsWrLMLuL7qfb0CeDkiDmZlJW3qKr8NeCqrRNOG7g7gLcClwEHg86UVI2JHRGyOiM3TSlLqmtno6FOva0TM0rmsvB94Erg3Ih6XdKOkG6vVdgPPAHuBrwD/PitblblN0mOSHgXeC3wsq0ejgOGIeO7Ea0lfodO9a2Zt0cf+keoe/u4F8+7seh3ATXXLVvN/fzl1aHRGV3X9nvB+4LEm2zGz0VP3snWcUjnVCS/5OnAlnd6V/XRu+F0p6VI67f4+4CMDrGP/pV39STxIodzETHl7EzPlzU0ktyx1tDzgy8tHyqEivzjSe4CZJ9adX/6wRJb1pBRCAvDEsd6fV6of5N8r2x/Zfizt/+w36/fxMZbWWuLNiLiux+yvDqAuZjYixulsrQ4/1G9mi7mhM7NWG7P7b3W4oTOzxdzQmVnbZc9ojyNnLzGz1hudM7qsa77w5yXSDBTJn6Qsy8RcudzEzFxhfnl7U8eSZUfKXfhTr5b/Bh1+aWNx2c/WnVNcVnJo4+nFZdlgNlkmklIYyc9+Xa5f9r2mk/0xlSRfKe3/PCSo9+8MpMdHelwlx2PxOF7NUBZfuppZq7kzwszWBDd0ZtZ6bujMrM1E+3pd3dCZ2cl8j27EZL1SWe/YbHmsA82We9xU6I2bPFIuM/16ubdw+tVyr+vchvKy+enp4rLn6d2DeuR4ucy+jeUH7ddPlfdVNsZD6QH9rGd18lflOq57ubw/pl8t/6+cfr33cZD9ZqXfGfLjIzuu0uMxO45Xixs6M2s9N3Rm1na+dDWz9nNDZ2atFu51NbO1wGd0ZtZ2vke3EgFReBhZE9lDz71DC9T0Ieqkq18zSejJ0d4DEEyuK+/GqcPl8JL10+WQiZjMEsuUlx2fWddz/muvlet4eMMp5Y+azBIglOtfGuMhezg/CyFZ/2K5HutfKf+eU4cLIUFHl/87A5AcH9lxlR2PpWQAWdKK8vb61EK5oTOzVqs5Zus4cUNnZicRvnQ1szXADZ2ZtZ8bOjNrvZY1dB4zwsxOVmUvqTPVIWmrpKcl7ZV0S4/lkvTlavmjki5bqqykv5D0VLX+fZLOyOqw5BmdpJ3A7wCHIuJt1byzgG8AFwH7gGsi4tf1vnYDpfz4Sfh2GkIy2SwDhY73DjvQkSS8ZKr8tyQmy+EUGc2Vtzl5tPc2Z1/pHe4BML++vCwmknE55sv1nygMNZGN75BlIclCSNa9Uv49p17r/XvqSDmEpPQ7A3mGkvlyPdLQk1L2kkGPC5Hp00dLmgRuB94D7Af2SNoVEU90rXYVsKmatgB3AFuWKPsA8KcRMSvpPwN/CnyyVI86Z3R/DWxdMO8W4MGI2AQ8WL03s5bQfL2phsuBvRHxTEQcB+4Bti1YZxtwd3Q8BJwh6bysbET8fUSc+KvzEHBBVoklG7qI+B7wYo+K3VW9vgt431LbMbPx0cdL1/OBZ7ve76/m1VmnTlmAPwL+Z1aJpp0R50bEwer1L4FzSytK2g5sB9hAEoFvZqNheQHDZ0t6uOv9jojY0fc6FUj6M2AW+Fq23op7XSMipHLbXn3pHQCn66yW9eWYtVT9/6kvRMTmZPkB4MKu9xdU8+qsM52VlfRBOv0H7470ebnmva7PVdfQVP8eargdMxsxJ56M6NOl6x5gk6SLJa0DrgV2LVhnF3B91ft6BfBydcVYLCtpK/AnwO9GxOtLVaLpGd0u4AbgturfbzfcjpmNIPWpx7fqFb0ZuB+YBHZGxOOSbqyW3wnsBq4G9gKvAx/Kylab/m/AeuABSQAPRcSNpXrUCS/5OnAlnWvx/cCn6DRw90r6MPBz4Jrlff3FSllNAFSIfsjOVpUNOJJloFCSkWOyd0V0tBBLAaRJSBKaK3+3yZlyOMjMkd4fOLs+G2ynXI9I6p/9Z5goRGhMHVv+QDZQzkIC5RASgMnXev822W+WHh8zzTKbZAPgLHHV1btMad/3o33q80P9EbGbTmPWPe/OrtcB3FS3bDX/ny+nDks2dBFxXWHRu5fzQWY2Pvysq5m1nxs6M2s7n9GZWfu5oTOzVvMoYGbWds4w3A+lQT2UxDGUyswnZeaSTBLlUmg2GbDm2PHeZZLtZSaTv5oTs+WFk8eSbCmv9w49mU8G4smWpeE2SVjExEzvZaX5AJNHyr9ZOphNlomkEEaio71/Syj/zkCavSSSYy47HotZSrIBdQatQcjLKPMZnZkt4jM6M2s3jwJmZmuBOyPMrPXc0JlZuwXujDCz9nNnxCooZWrQRJYRollYRDqojnqHFqThKtmy5LMiyYSh4+V0I5NHSuEl5YwnTGThJeVi6X+Gwm82MVMOs1C27GjDwWwK+7F5CElyTZcdV9kZUiGMJMvoM3Bu6MyszRwwbGbtF9G3xJujwg2dmS3WrnbODZ2ZLeZLVzNrt6D8/O2YGp2GLnuAOXvgvyT5oYKkdy/ZZGmLSuoeyTLNlXtPNZv1Tia9glO9e1cnssErJsrLmj7UXxwjIevVTr5zOh5D1ktaGuOhac9qliyiyYP7TQ36gf92tXMj1NCZ2cjwpauZtZ57Xc2s3Zy9xMzarhMw3K6Wzg2dmS3m7CVm1nY+o+siaR/wKjAHzEbE5iULlXZgEsZQfui5HBaRPfCfjTXRJPQkSqEU5Dd10we90/CS5AH3qcJPmu3f5KF+pb9LFl5SGgch+84NQz7mlx+WkoaCpEkfGoaQZCFIozZmhO/R9fSuiHihD9sxs5HQvmddG0TimlnrRdSbapC0VdLTkvZKuqXHckn6crX8UUmXLVVW0h9IelzSvKQlryRX2tAF8A+SfiRpe68VJG2X9LCkh2foPfScmY2QagDrOtNSJE0CtwNXAZcA10m6ZMFqVwGbqmk7cEeNso8Bvwd8r85XWuml6zsj4oCkNwEPSHoqIk764IjYAewAOF1ntet82Kyt+tcZcTmwNyKeAZB0D7ANeKJrnW3A3dG5af2QpDMknQdcVCobEU9W82pVYkVndBFxoPr3EHBf9aXMbNxFzQnOPnHFVk0Lr+zOB57ter+/mldnnTpla2l8RifpVGAiIl6tXr8X+POm2zOz0aEkkmCBF2pFW6yylVy6ngvcV506TgF/ExF/15da1ZV22Q8v9CQ7fU4vAJLLA00mWTJmkxPxUohGkwwwkIaeNMrIkYVMpOMxJL91GnpSyqKS7N8mYTPQLIRkiXKrIuhnwPAB4MKu9xdU8+qsM12jbC2NG7rquvntTcub2WgS0c+A4T3AJkkX02mkrgX+cME6u4Cbq3twW4CXI+KgpOdrlK3FT0aY2WJ9augiYlbSzcD9wCSwMyIel3RjtfxOYDdwNbAXeB34UFYWQNL7gf8KnAP8D0mPRMRvl+rhhs7MFuvjI2ARsZtOY9Y9786u1wHcVLdsNf8+Oh2gtbihM7OT9fce3UhwQ2dmiyyj13UsuKEzswXqP941Lkanoct2bM3o55O3N4DQk0JEQmShG1kISRaWMpdsM81E0rtc3Qjyxdvrb3hJHrqRhZc0K1f8vAGEeww1hGSQDVEMePurYHQaOjMbHe26cnVDZ2aLOfGmmbWfGzoza7WI/HG8MeSGzswW8xmdmbWeG7pV0GRAnXR7zUJPSppnQ8nCEZqFpZSyckQh7GToGgaiNs4oUtxgw0wjDbfZbHur1NgEzfbpCBuPhs7MhihGL3XUCrmhM7OTBe6MMLM1wPfozKz13NCZWbv5of7R0u9EAJDfhC08vJ/30pXHJlD2wHy5WJ5EoFhoRA7cIfdoNu5B7XM98m2OyG9zQtC4d3xUjXdDZ2aDMWqN7wq5oTOzBfwImJm1XUA4js7MWs9PRphZ6/kenZm1WoR7XbtJ2gp8ic7gsn8VEbf1pVb90PQvUhaW0uS+RRIK0jT0IU0iUJKEq4wDh4kMWZu+Cyto6CRNArcD7wH2A3sk7YqIJ/pVOTNbDUEUsuCMq5Xk7rkc2BsRz0TEceAeYFt/qmVmq+ZEmqY605hYyaXr+cCzXe/3A1sWriRpO7AdYAOnrODjzGxoHF6yPBGxA9gBcLrOGp8/AWZrVDCAe6KrbCUN3QHgwq73F1TzzGychRNvdtsDbJJ0MZ0G7lrgD/tSKzNbVW3rjFCai3+pwtLVwH+hE16yMyI+u8T6zwM/r96eDbzQ+MP7x/U4metxsnGrxz+NiHNW8kGS/q76vDpeiIitK/m8YVhRQ7eiD5YejojNq/Lhrofr4XqsKSMyNJSZ2eC4oTOz1lvNhm7HKn52N9fjZK7HyVyPFli1e3RmZsPiS1cza71VaegkbZX0tKS9km5ZjTpU9dgn6aeSHpH08BA/d6ekQ5Ie65p3lqQHJP1j9e+Zq1SPT0s6UO2TR6oQokHX40JJ35X0hKTHJX2smj/UfZLUY6j7RNIGST+U9JOqHv+pmj/0Y6Qthn7pWmU9+RldWU+A61Yj64mkfcDmiBhqnJSkfw0cBu6OiLdV8z4HvBgRt1WN/5kR8clVqMengcMR8ZeD/OwF9TgPOC8ifizpDcCPgPcBH2SI+ySpxzUMcZ9IEnBqRByWNA18H/gY8HsM+Rhpi9U4o1vzWU8i4nvAiwtmbwPuql7fRec/2GrUY+gi4mBE/Lh6/SrwJJ2kEUPdJ0k9hio6Dldvp6spWIVjpC1Wo6HrlfVk6AdTJYB/kPSjKsvKajo3Ig5Wr38JnLuKdfmopEerS9uhXh5Jugh4B/ADVnGfLKgHDHmfSJqU9AhwCHggIlZ1f4y7td4Z8c6IuBS4CripupRbddG5n7Ba3eF3AG8BLgUOAp8f1gdLOg34JvDxiHile9kw90mPegx9n0TEXHVsXgBcLultC5av5jEydlajoRuZrCcRcaD69xBwH53L6tXyXHWP6MS9okOrUYmIeK76TzYPfIUh7ZPqXtQ3ga9FxLeq2UPfJ73qsVr7pPrsl4DvAlsZkWNkHK1GQ/f/sp5IWkcn68muYVdC0qnVDWcknQq8F3gsLzVQu4Abqtc3AN9ejUqc+I9UeT9D2CfVzfevAk9GxBe6Fg11n5TqMex9IukcSWdUrzfS6bh7ihE5RsZ1EPhUAAAAp0lEQVTRqgQMLzfryYDq8BY6Z3HQSVf1N8Oqh6SvA1fSyRDxHPAp4G+Be4E308nwck1EDLSjoFCPK+lcogWwD/hI132hQdXjncD/Bn4KnEiEdiud+2ND2ydJPa5jiPtE0m/S6WyYpHMycm9E/Lmk32DIx0hb+MkIM2u9td4ZYWZrgBs6M2s9N3Rm1npu6Mys9dzQmVnruaEzs9ZzQ2dmreeGzsxa7/8CWPT2TsLo8/sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e07a7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(makeImage1())\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:32.640826Z",
     "start_time": "2018-04-11T17:46:32.594276Z"
    }
   },
   "outputs": [],
   "source": [
    "# This is the code to compute the Hankel transform, alternatively\n",
    "# known as the Fourier-Bessel transform, which allows one to compute\n",
    "# the Kolmogorov surface brightness at a given radius by doing a 1D\n",
    "# integral.  We can tabulate such results over a range of radii, and\n",
    "# interpolate this onto our desired output grid.\n",
    "\n",
    "def KolmogorovIntegrand(k, x):\n",
    "    return j0(x*k)*kolmogorovKValue(k)*k\n",
    "\n",
    "def KolmogorovXValue(x):\n",
    "    return quad(KolmogorovIntegrand, 0, np.inf, x)[0]/(2*np.pi)\n",
    "\n",
    "def makeImage2():\n",
    "    dx = np.pi/maxk\n",
    "    xmaxTable = dx*subregion\n",
    "    xTable = np.arange(0, xmaxTable, dx)\n",
    "    xVals = [KolmogorovXValue(x) for x in xTable]\n",
    "    xInterp = interp1d(xTable, xVals, 'cubic', \n",
    "                       fill_value=(xVals[0], 0),\n",
    "                       bounds_error=False)    \n",
    "    \n",
    "    xmax = nk//2*dx\n",
    "    xs = np.linspace(-xmax, xmax, nk, endpoint=False)\n",
    "    slc = slice((nk-subregion)//2, (nk+subregion)//2)\n",
    "    xs = xs[slc]\n",
    "    xs, ys = np.meshgrid(xs, xs)\n",
    "    rs = np.hypot(xs, ys)\n",
    "    return xInterp(rs)*dx**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:33.056405Z",
     "start_time": "2018-04-11T17:46:32.642975Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x10ba7f518>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAD8CAYAAADnhGhBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+sJWWd5/H35/7oH4AIDEhYwAV322yIGZH00iTjbnCNTkMm0zqTITBZQcdMyy4YTdyMDJOs7hgT1hl1dZdA2rEzkDgiWWXsuL3DMMSNa7Joo0Hkp9MhrXSnpUHkR0P/uD+++8ep3j1973m+t27dc849p+7nlVT6nKp66jynTvVzq+r51vdRRGBm1mYTq10BM7NBc0NnZq3nhs7MWs8NnZm1nhs6M2s9N3Rm1npu6MxsoCRtlfS0pL2SbumxXJK+XC1/VNJldcpK+qikpyQ9LulzWR2m+vd1zMxOJmkSuB14D7Af2CNpV0Q80bXaVcCmatoC3AFsycpKehewDXh7RByT9KasHj6jM7NBuhzYGxHPRMRx4B46DVS3bcDd0fEQcIak85Yo+++A2yLiGEBEHMoqMdQzunVaHxs4dZgfOdo09IJjrOETPGvswZ+jvMbxOLaiA+S333Vq/OrFuVrr/ujRY/dHxNZklfOBZ7ve76dz1rbUOucvUfatwL+S9FngKPAfImJPqRJDbeg2cCpb9O7hfJgG0BiovyfAmmhYxyb1aPpZ/TbftMGab1iszy1dw3rk2+xfHX8QD654G796cY4f3v/mWutOnveP/0LSw12zdkTEjhVXYmlTwFnAFcC/BO6V9JYoPNPqe3RmdpIA5qndoL8QEZuT5QeAC7veX1DNq7POdFJ2P/CtqmH7oaR54Gzg+V6V8D06MztJEMzEXK2phj3AJkkXS1oHXAvsWrDOLuD6qvf1CuDliDi4RNm/Bd4FIOmtwDrghVIlfEZnZoss44wuFRGzkm4G7gcmgZ0R8bikG6vldwK7gauBvcDrwIeystWmdwI7JT0GHAduKF22ghs6M1sgCOb6eN8wInbTacy6593Z9TqAm+qWreYfB/5t3Tq4oTOzReZb1l29ZEMnaQPwPWB9tf5/j4hPSToL+AZwEbAPuCYifj24qvas3AC2ufzblmnvadOe2mSbavK9J0bkdqwa9p5G8p2TnlVNljZYrkfaU5v9nk17ZEu/5yolxQ1grmUNXZ2j/xjwbyLi7cClwNbqhuEtwIMRsQl4sHpvZi0wT9SaxsWSDV0VrXy4ejtdTUEnQvmuav5dwPsGUkMzG6oAZiJqTeOi1vWMpElJjwCHgAci4gfAuVUXMMAvgXMHVEczG6IgmKs5jYtanRERMQdcKukM4D5Jb1uwPCT1/NaStgPbATZwygqra2YDFzA3Pm1YLcu6Qx0RLwHfBbYCz1UP3lL92/Oh2ojYERGbI2LzNOtXWl8zG7DOkxH1pnGxZEMn6ZzqTA5JG+mkTHmKToTyDdVqNwDfHlQlzWyYxFzNaVzUuXQ9D7iryg01AdwbEd+R9H/oPEj7YeDnwDUDq2W/w0iSEIFGD9pnIQdNw0SycJAG5RqFpECeDKDBA/PpOMLz5XMEZeWSkJXi581nx0DTBAJ9Dj3JfrMBdgR0OiPGpxGrY8mGLiIeBd7RY/6vgCGlIjGzYenE0a2xhs7M1p75tXZGZ2Zri8/ozKz1AjHXsgxubujMbBFfug7KqPSsNuhBTXs0J0tPlZP2rGoy63VNlpXKDSC5QJNeV2W9j3PJsuwh/KScSj25c+WkkY0TCCS9tUPrke1DZ2wgjkdy3I6h0WnozGwkdAKGfelqZi3nzggza7UIMRc+ozOzlpv3GZ2ZtVmnM6JdTUO7vo2ZrZg7I0bNEENIAFQKFUnCS4plAKaT3T9RLpeGnkwVtpkmAsj2VcMHy0thGFmZ2dnyR2UhJJPJ+KIzvbeZRWFkCQSC5LMaJgoohp40HYOiD+YcR2dmbeYnI8xsTZh3r6uZtVnnoX43dGbWYoGY8SNgZtZmEThg2MzaTg4YXrEmWUoK4SBDDSGBYiaSRuEegKanG5XLwlJiqlD/rI5JFpXIQmeajP+QhYnMlkM3VAgTAfKwlEL9NVv+XmkoS7kWfQ896XvGk5qC/p7RSdoKfAmYBP4qIm5bsFzV8quB14EPRsSPs7KSPg38MfB8tZlbI2J3qQ4+ozOzRfrVGVENqnU7ndED9wN7JO2KiCe6VrsK2FRNW4A7gC01yn4xIv6yTj3adSFuZisWiPmoN9VwObA3Ip6JiOPAPcC2BetsA+6OjoeAM6qxouuUrcUNnZmdpDPc4VStCThb0sNd0/YFmzsfeLbr/f5qXp11lir7UUmPStop6czsO/nS1cwWWNbg1C9ExOZB1qbgDuAzdNrlzwCfB/6otLIbOjM7SdDXJyMOABd2vb+gmldnnelS2Yh47sRMSV8BvpNVwpeuZrbIXHVWt9RUwx5gk6SLJa0DrgV2LVhnF3C9Oq4AXo6Ig1nZ6h7eCe8HHssqseQZnaQLgbuBc+k09jsi4kvL7d5dUtPBW0r6HEICoFLIRxZCsn5d+bOyMJF15dCT2JAsm+5d//nCfCDdV9n9ZmUpQArZSyZmshCSZNnRmfKy48mywu8Zx46XyygJVykuWUHoSbKo/GGDO0eJUN/O6CJiVtLNwP10QkR2RsTjkm6slt8J7KYTWrKXTnjJh7Ky1aY/J+lSOj/JPuAjWT3qXLrOAp+IiB9LegPwI0kPVMtqd++a2XjodEb07xGw6gRo94J5d3a9DuCmumWr+R9YTh2WbOiqU8iD1etXJT3J4l4TM2uN9o0ZsaxvI+ki4B3AD6pZtbt3zWw8dDoj+hZHNxJqN3SSTgO+CXw8Il6h0737FuBSOmd8ny+U234ixmaGY32ospkN2hwTtaZxUSu8RNI0nUbuaxHxLajfvRsRO4AdAKfrrD6MI25mg3TiyYg2qdPrKuCrwJMR8YWu+edV9++gRveumY2PtTg4zm8BHwB+KumRat6twHXL6d5diWKWkmxwnCxLSjqYzfIzkWQhJLEhW7a+vGxjOYRkbkP5Z5vbWAovKX/nbFnT7CUTM6XwknKZySPlOIvJdUkIz5Fk2dHet0vyUJAyJVlDopSxhfx4jGJ2nuyzBndxFAEzWeaUMVSn1/X79D4umsfMmdnI6ly6rrGGzszWnmU86zoW3NCZ2UlOhJe0iRs6M1vAl65mtgZ4zAgza7VOr6uHO1yZJlkXSmWywXGyLCRZ9pIGmUjSEJJTyiEkc6eWl82eWq7H7Gnl+s+c0ntfza7PwkuKi8iuYJSMzzJRSCgydawcFjH9evnDpg4ny6bKy7JooZI09CQJL1ES8pEGg5TCdNIQjyYpT+pZkwHDZrb2+NLVzFrNva5mtia419XMWi1CzLqhM7O286XrgBQf3M/KZA/uTyR/kZKxGjSddEEWymUP52c9qzOnlz/r+OnlntVjp5e/28wbeu+T2Y3FIsyXq0hMlPsLNV/e/xOF1INTR8plpl8tL1ufJR6YXP6xM5n0GGsu6VmdS8brSJIcFHtWKR/HkfTVFv+/JN+rLt+jM7M1wQ2dmbWa4+jMbE1wHJ2ZtVoEzK61xJtmtvb40tXMWs336FZLqSs9CSFJx36YaPZQf6zrHVqQje+QPZyfhZAcPbNc/2NnlQ/C42/sHZIw+4Zk/IENyQPik0nIxFwynsTR3t9t6tXy95rbkIWQNLuU0lxh7IrZZH/MzJa3N5vsq2SZJsvLYq7w3bJwlcE90199tBs6M2u5tnVGtOuOo5mtWETnHl2dqQ5JWyU9LWmvpFt6LJekL1fLH5V02TLKfkJSSDo7q4PP6MxsATHXp15XSZPA7cB7gP3AHkm7IuKJrtWuAjZV0xbgDmDLUmUlXQi8F/jFUvXwGZ2ZLRKhWlMNlwN7I+KZiDgO3ANsW7DONuDu6HgIOEPSeTXKfhH4E5bIawpu6MxsgRPPuta8dD1b0sNd0/YFmzsfeLbr/f5qXp11imUlbQMORMRP6nwnX7qa2cki7/Bd4IWI2DzA2iwi6RTgVjqXrbUs2dBV18F3A+fSaex3RMSXJJ0FfAO4CNgHXBMRv843lmRdSMaSKGYpybKXZNvLQhWSzCaxoXcYydyGZuM7ZFlIshCSo2eXQyPmfqP3YA2nnXGkWOaNG48Wl62fKodaHJstf++Xj2zoOf/wS+U0KvNZ5pjk4kOl8Axgcqb3/p88lmSwOV6uh7LQk5nCQBlAzGaDb/T+rbPsPFE8vvvTW9rHXtcDwIVd7y+o5tVZZ7ow/58BFwM/qfbRBcCPJV0eEb/sVYk6l66zwCci4hLgCuAmSZcAtwAPRsQm4MHqvZmNuag6I+pMNewBNkm6WNI64Fpg14J1dgHXV72vVwAvR8TBUtmI+GlEvCkiLoqIi+hc0l5WauSgxhld9YEHq9evSnqSznXyNuDKarW7gP8FfLLONzez0baMS9clthOzkm4G7gcmgZ0R8bikG6vldwK7gauBvcDrwIeysk3qsax7dJIuAt4B/AA4t2oEAX5J59LWzFqgn09GRMRuOo1Z97w7u14HcFPdsj3WuWipOtRu6CSdBnwT+HhEvNJ9/yAiQlLPvwFVL8x2gA2cUvfjzGyVRLTvEbBaF9mSpuk0cl+LiG9Vs5+rYl2o/j3Uq2xE7IiIzRGxeVpJzm4zGxn9fDJiFCzZ0Klz6vZV4MmI+ELXol3ADdXrG4Bv9796ZrYaIupN46LOpetvAR8AfirpkWrercBtwL2SPgz8HLhmMFVMZAPgZCEkWYaSqXI4SEz3Xja3sVxm5pTlD2QD5SwkUA4hATjnTa/0nP/WM58vlnnzxheLy06bLIxyAxyeK5+h/+LIWT3n/2zdOcUyz3N6cdnxmXXFZZNHy/tx5kjv/T/1evk3mzySHAPJ8aHkuGK2HJZSPI7nBpyipCAQ82st8WZEfJ9ycM67+1sdMxsFY3SyVoufjDCzk7WwM8INnZkt1rJTOjd0ZraIz+jMrNUCmJ93Q2dmbRaAz+gGpJTVBIrd71l2hyx7SZr1JAlLmS+El8xPl7c3uz5ZVk7kkQ5mk2UiKYWRbH7jvmKZS9YvTCbx/50+Uc5s8sp87wwlAE+sW5hybGlHkqwhr72WZIh5pRzyUdr/2W9W+p0BJrKwpabZdArlIgufGnAQ2zjFyNUxOg2dmY0ON3Rm1m6106SPDTd0ZraYz+jMrNUCwr2uZtZ+bujGX4Me3qxc3oNX3tx8krUqNpQf6M7GeCg9oJ/1rF6yrjzUxxkT5UPkpfly72/JoY3lB/f3beydCADg8IZyLsP59eVe0tL+z36zfh8fY8mXrmbWem7ozKzVHDBsZmuBA4bNrP3c62pmbdd7qKvx5YbOzE4WuDNipKRhANkD1uVlkS5rUCZ7LnsiOZomy8vWT5XHHyiN8ZA9nJ+FkJw2UX5wH8rbLH1eNgZF9r2y/ZHtx9L+b/I7L1UuTTLR9FhdFXJnhJmtAT6jM7PWK2cJG0tu6MzsZC2Mo2vX4I1m1heKelOtbUlbJT0taa+kW3osl6QvV8sflXTZUmUlfaZa9xFJfy/pn2R1cENnZotFzWkJkiaB24GrgEuA6yRdsmC1q4BN1bQduKNG2b+IiN+MiEuB7wD/MauHGzozG6TLgb0R8UxEHAfuAbYtWGcbcHd0PAScIem8rGxEvNJV/lSWaHaXvEcnaSfwO8ChiHhbNe/TwB8DJwYouDUidi+1rb6bT75btix5vkXpsgZlkpu6yqLP58rLjs2Wf7bDc71TomTjO+RZSMohJC/Nl8NBXpl/Q8/5pfpB/r2y/ZHtRxWOgya/81Ll0uemmh6rq6SPAcPnA892vd8PbKmxzvlLlZX0WeB64GXgXVkl6pzR/TWwtcf8L0bEpdU0/EbOzAYj6DwCVmeCsyU93DVtH1o1I/4sIi4EvgbcnK275BldRHxP0kX9qZqZjYX6Z3QvRMTmZPkB4MKu9xdU8+qsM12jLHQaut3Ap0qVWMk9uo9WvR47JZ1ZWknS9hOt/UyUo+LNbHT0sdd1D7BJ0sWS1gHXArsWrLMLuL7qfb0CeDkiDmZlJW3qKr8NeCqrRNOG7g7gLcClwEHg86UVI2JHRGyOiM3TSlLqmtno6FOva0TM0rmsvB94Erg3Ih6XdKOkG6vVdgPPAHuBrwD/PitblblN0mOSHgXeC3wsq0ejgOGIeO7Ea0lfodO9a2Zt0cf+keoe/u4F8+7seh3ATXXLVvN/fzl1aHRGV3X9nvB+4LEm2zGz0VP3snWcUjnVCS/5OnAlnd6V/XRu+F0p6VI67f4+4CMDrGP/pV39STxIodzETHl7EzPlzU0ktyx1tDzgy8tHyqEivzjSe4CZJ9adX/6wRJb1pBRCAvDEsd6fV6of5N8r2x/Zfizt/+w36/fxMZbWWuLNiLiux+yvDqAuZjYixulsrQ4/1G9mi7mhM7NWG7P7b3W4oTOzxdzQmVnbZc9ojyNnLzGz1hudM7qsa77w5yXSDBTJn6Qsy8RcudzEzFxhfnl7U8eSZUfKXfhTr5b/Bh1+aWNx2c/WnVNcVnJo4+nFZdlgNlkmklIYyc9+Xa5f9r2mk/0xlSRfKe3/PCSo9+8MpMdHelwlx2PxOF7NUBZfuppZq7kzwszWBDd0ZtZ6bujMrM1E+3pd3dCZ2cl8j27EZL1SWe/YbHmsA82We9xU6I2bPFIuM/16ubdw+tVyr+vchvKy+enp4rLn6d2DeuR4ucy+jeUH7ddPlfdVNsZD6QH9rGd18lflOq57ubw/pl8t/6+cfr33cZD9ZqXfGfLjIzuu0uMxO45Xixs6M2s9N3Rm1na+dDWz9nNDZ2atFu51NbO1wGd0ZtZ2vke3EgFReBhZE9lDz71DC9T0Ieqkq18zSejJ0d4DEEyuK+/GqcPl8JL10+WQiZjMEsuUlx2fWddz/muvlet4eMMp5Y+azBIglOtfGuMhezg/CyFZ/2K5HutfKf+eU4cLIUFHl/87A5AcH9lxlR2PpWQAWdKK8vb61EK5oTOzVqs5Zus4cUNnZicRvnQ1szXADZ2ZtZ8bOjNrvZY1dB4zwsxOVmUvqTPVIWmrpKcl7ZV0S4/lkvTlavmjki5bqqykv5D0VLX+fZLOyOqw5BmdpJ3A7wCHIuJt1byzgG8AFwH7gGsi4tf1vnYDpfz4Sfh2GkIy2SwDhY73DjvQkSS8ZKr8tyQmy+EUGc2Vtzl5tPc2Z1/pHe4BML++vCwmknE55sv1nygMNZGN75BlIclCSNa9Uv49p17r/XvqSDmEpPQ7A3mGkvlyPdLQk1L2kkGPC5Hp00dLmgRuB94D7Af2SNoVEU90rXYVsKmatgB3AFuWKPsA8KcRMSvpPwN/CnyyVI86Z3R/DWxdMO8W4MGI2AQ8WL03s5bQfL2phsuBvRHxTEQcB+4Bti1YZxtwd3Q8BJwh6bysbET8fUSc+KvzEHBBVoklG7qI+B7wYo+K3VW9vgt431LbMbPx0cdL1/OBZ7ve76/m1VmnTlmAPwL+Z1aJpp0R50bEwer1L4FzSytK2g5sB9hAEoFvZqNheQHDZ0t6uOv9jojY0fc6FUj6M2AW+Fq23op7XSMipHLbXn3pHQCn66yW9eWYtVT9/6kvRMTmZPkB4MKu9xdU8+qsM52VlfRBOv0H7470ebnmva7PVdfQVP8eargdMxsxJ56M6NOl6x5gk6SLJa0DrgV2LVhnF3B91ft6BfBydcVYLCtpK/AnwO9GxOtLVaLpGd0u4AbgturfbzfcjpmNIPWpx7fqFb0ZuB+YBHZGxOOSbqyW3wnsBq4G9gKvAx/Kylab/m/AeuABSQAPRcSNpXrUCS/5OnAlnWvx/cCn6DRw90r6MPBz4Jrlff3FSllNAFSIfsjOVpUNOJJloFCSkWOyd0V0tBBLAaRJSBKaK3+3yZlyOMjMkd4fOLs+G2ynXI9I6p/9Z5goRGhMHVv+QDZQzkIC5RASgMnXev822W+WHh8zzTKbZAPgLHHV1btMad/3o33q80P9EbGbTmPWPe/OrtcB3FS3bDX/ny+nDks2dBFxXWHRu5fzQWY2Pvysq5m1nxs6M2s7n9GZWfu5oTOzVvMoYGbWds4w3A+lQT2UxDGUyswnZeaSTBLlUmg2GbDm2PHeZZLtZSaTv5oTs+WFk8eSbCmv9w49mU8G4smWpeE2SVjExEzvZaX5AJNHyr9ZOphNlomkEEaio71/Syj/zkCavSSSYy47HotZSrIBdQatQcjLKPMZnZkt4jM6M2s3jwJmZmuBOyPMrPXc0JlZuwXujDCz9nNnxCooZWrQRJYRollYRDqojnqHFqThKtmy5LMiyYSh4+V0I5NHSuEl5YwnTGThJeVi6X+Gwm82MVMOs1C27GjDwWwK+7F5CElyTZcdV9kZUiGMJMvoM3Bu6MyszRwwbGbtF9G3xJujwg2dmS3WrnbODZ2ZLeZLVzNrt6D8/O2YGp2GLnuAOXvgvyT5oYKkdy/ZZGmLSuoeyTLNlXtPNZv1Tia9glO9e1cnssErJsrLmj7UXxwjIevVTr5zOh5D1ktaGuOhac9qliyiyYP7TQ36gf92tXMj1NCZ2cjwpauZtZ57Xc2s3Zy9xMzarhMw3K6Wzg2dmS3m7CVm1nY+o+siaR/wKjAHzEbE5iULlXZgEsZQfui5HBaRPfCfjTXRJPQkSqEU5Dd10we90/CS5AH3qcJPmu3f5KF+pb9LFl5SGgch+84NQz7mlx+WkoaCpEkfGoaQZCFIozZmhO/R9fSuiHihD9sxs5HQvmddG0TimlnrRdSbapC0VdLTkvZKuqXHckn6crX8UUmXLVVW0h9IelzSvKQlryRX2tAF8A+SfiRpe68VJG2X9LCkh2foPfScmY2QagDrOtNSJE0CtwNXAZcA10m6ZMFqVwGbqmk7cEeNso8Bvwd8r85XWuml6zsj4oCkNwEPSHoqIk764IjYAewAOF1ntet82Kyt+tcZcTmwNyKeAZB0D7ANeKJrnW3A3dG5af2QpDMknQdcVCobEU9W82pVYkVndBFxoPr3EHBf9aXMbNxFzQnOPnHFVk0Lr+zOB57ter+/mldnnTpla2l8RifpVGAiIl6tXr8X+POm2zOz0aEkkmCBF2pFW6yylVy6ngvcV506TgF/ExF/15da1ZV22Q8v9CQ7fU4vAJLLA00mWTJmkxPxUohGkwwwkIaeNMrIkYVMpOMxJL91GnpSyqKS7N8mYTPQLIRkiXKrIuhnwPAB4MKu9xdU8+qsM12jbC2NG7rquvntTcub2WgS0c+A4T3AJkkX02mkrgX+cME6u4Cbq3twW4CXI+KgpOdrlK3FT0aY2WJ9augiYlbSzcD9wCSwMyIel3RjtfxOYDdwNbAXeB34UFYWQNL7gf8KnAP8D0mPRMRvl+rhhs7MFuvjI2ARsZtOY9Y9786u1wHcVLdsNf8+Oh2gtbihM7OT9fce3UhwQ2dmiyyj13UsuKEzswXqP941Lkanoct2bM3o55O3N4DQk0JEQmShG1kISRaWMpdsM81E0rtc3Qjyxdvrb3hJHrqRhZc0K1f8vAGEeww1hGSQDVEMePurYHQaOjMbHe26cnVDZ2aLOfGmmbWfGzoza7WI/HG8MeSGzswW8xmdmbWeG7pV0GRAnXR7zUJPSppnQ8nCEZqFpZSyckQh7GToGgaiNs4oUtxgw0wjDbfZbHur1NgEzfbpCBuPhs7MhihGL3XUCrmhM7OTBe6MMLM1wPfozKz13NCZWbv5of7R0u9EAJDfhC08vJ/30pXHJlD2wHy5WJ5EoFhoRA7cIfdoNu5B7XM98m2OyG9zQtC4d3xUjXdDZ2aDMWqN7wq5oTOzBfwImJm1XUA4js7MWs9PRphZ6/kenZm1WoR7XbtJ2gp8ic7gsn8VEbf1pVb90PQvUhaW0uS+RRIK0jT0IU0iUJKEq4wDh4kMWZu+Cyto6CRNArcD7wH2A3sk7YqIJ/pVOTNbDUEUsuCMq5Xk7rkc2BsRz0TEceAeYFt/qmVmq+ZEmqY605hYyaXr+cCzXe/3A1sWriRpO7AdYAOnrODjzGxoHF6yPBGxA9gBcLrOGp8/AWZrVDCAe6KrbCUN3QHgwq73F1TzzGychRNvdtsDbJJ0MZ0G7lrgD/tSKzNbVW3rjFCai3+pwtLVwH+hE16yMyI+u8T6zwM/r96eDbzQ+MP7x/U4metxsnGrxz+NiHNW8kGS/q76vDpeiIitK/m8YVhRQ7eiD5YejojNq/Lhrofr4XqsKSMyNJSZ2eC4oTOz1lvNhm7HKn52N9fjZK7HyVyPFli1e3RmZsPiS1cza71VaegkbZX0tKS9km5ZjTpU9dgn6aeSHpH08BA/d6ekQ5Ie65p3lqQHJP1j9e+Zq1SPT0s6UO2TR6oQokHX40JJ35X0hKTHJX2smj/UfZLUY6j7RNIGST+U9JOqHv+pmj/0Y6Qthn7pWmU9+RldWU+A61Yj64mkfcDmiBhqnJSkfw0cBu6OiLdV8z4HvBgRt1WN/5kR8clVqMengcMR8ZeD/OwF9TgPOC8ifizpDcCPgPcBH2SI+ySpxzUMcZ9IEnBqRByWNA18H/gY8HsM+Rhpi9U4o1vzWU8i4nvAiwtmbwPuql7fRec/2GrUY+gi4mBE/Lh6/SrwJJ2kEUPdJ0k9hio6Dldvp6spWIVjpC1Wo6HrlfVk6AdTJYB/kPSjKsvKajo3Ig5Wr38JnLuKdfmopEerS9uhXh5Jugh4B/ADVnGfLKgHDHmfSJqU9AhwCHggIlZ1f4y7td4Z8c6IuBS4CripupRbddG5n7Ba3eF3AG8BLgUOAp8f1gdLOg34JvDxiHile9kw90mPegx9n0TEXHVsXgBcLultC5av5jEydlajoRuZrCcRcaD69xBwH53L6tXyXHWP6MS9okOrUYmIeK76TzYPfIUh7ZPqXtQ3ga9FxLeq2UPfJ73qsVr7pPrsl4DvAlsZkWNkHK1GQ/f/sp5IWkcn68muYVdC0qnVDWcknQq8F3gsLzVQu4Abqtc3AN9ejUqc+I9UeT9D2CfVzfevAk9GxBe6Fg11n5TqMex9IukcSWdUrzfS6bh7ihE5RsZ1EPhUAAAAp0lEQVTRqgQMLzfryYDq8BY6Z3HQSVf1N8Oqh6SvA1fSyRDxHPAp4G+Be4E308nwck1EDLSjoFCPK+lcogWwD/hI132hQdXjncD/Bn4KnEiEdiud+2ND2ydJPa5jiPtE0m/S6WyYpHMycm9E/Lmk32DIx0hb+MkIM2u9td4ZYWZrgBs6M2s9N3Rm1npu6Mys9dzQmVnruaEzs9ZzQ2dmreeGzsxa7/8CWPT2TsLo8/sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1052ef710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(makeImage1())\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:33.353178Z",
     "start_time": "2018-04-11T17:46:33.058869Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x10c67a0b8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAD8CAYAAADnhGhBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+sJWWd5/H35/7oH4AIDEhYwAV322yIGZH00iTjbnCNTkMm0zqTITBZQcdMyy4YTdyMDJOs7hgT1hl1dZdA2rEzkDgiWWXsuL3DMMSNa7Joo0Hkp9MhrXSnpUHkR0P/uD+++8ep3j1973m+t27dc849p+7nlVT6nKp66jynTvVzq+r51vdRRGBm1mYTq10BM7NBc0NnZq3nhs7MWs8NnZm1nhs6M2s9N3Rm1npu6MxsoCRtlfS0pL2SbumxXJK+XC1/VNJldcpK+qikpyQ9LulzWR2m+vd1zMxOJmkSuB14D7Af2CNpV0Q80bXaVcCmatoC3AFsycpKehewDXh7RByT9KasHj6jM7NBuhzYGxHPRMRx4B46DVS3bcDd0fEQcIak85Yo+++A2yLiGEBEHMoqMdQzunVaHxs4dZgfOdo09IJjrOETPGvswZ+jvMbxOLaiA+S333Vq/OrFuVrr/ujRY/dHxNZklfOBZ7ve76dz1rbUOucvUfatwL+S9FngKPAfImJPqRJDbeg2cCpb9O7hfJgG0BiovyfAmmhYxyb1aPpZ/TbftMGab1iszy1dw3rk2+xfHX8QD654G796cY4f3v/mWutOnveP/0LSw12zdkTEjhVXYmlTwFnAFcC/BO6V9JYoPNPqe3RmdpIA5qndoL8QEZuT5QeAC7veX1DNq7POdFJ2P/CtqmH7oaR54Gzg+V6V8D06MztJEMzEXK2phj3AJkkXS1oHXAvsWrDOLuD6qvf1CuDliDi4RNm/Bd4FIOmtwDrghVIlfEZnZoss44wuFRGzkm4G7gcmgZ0R8bikG6vldwK7gauBvcDrwIeystWmdwI7JT0GHAduKF22ghs6M1sgCOb6eN8wInbTacy6593Z9TqAm+qWreYfB/5t3Tq4oTOzReZb1l29ZEMnaQPwPWB9tf5/j4hPSToL+AZwEbAPuCYifj24qvas3AC2ufzblmnvadOe2mSbavK9J0bkdqwa9p5G8p2TnlVNljZYrkfaU5v9nk17ZEu/5yolxQ1grmUNXZ2j/xjwbyLi7cClwNbqhuEtwIMRsQl4sHpvZi0wT9SaxsWSDV0VrXy4ejtdTUEnQvmuav5dwPsGUkMzG6oAZiJqTeOi1vWMpElJjwCHgAci4gfAuVUXMMAvgXMHVEczG6IgmKs5jYtanRERMQdcKukM4D5Jb1uwPCT1/NaStgPbATZwygqra2YDFzA3Pm1YLcu6Qx0RLwHfBbYCz1UP3lL92/Oh2ojYERGbI2LzNOtXWl8zG7DOkxH1pnGxZEMn6ZzqTA5JG+mkTHmKToTyDdVqNwDfHlQlzWyYxFzNaVzUuXQ9D7iryg01AdwbEd+R9H/oPEj7YeDnwDUDq2W/w0iSEIFGD9pnIQdNw0SycJAG5RqFpECeDKDBA/PpOMLz5XMEZeWSkJXi581nx0DTBAJ9Dj3JfrMBdgR0OiPGpxGrY8mGLiIeBd7RY/6vgCGlIjGzYenE0a2xhs7M1p75tXZGZ2Zri8/ozKz1AjHXsgxubujMbBFfug7KqPSsNuhBTXs0J0tPlZP2rGoy63VNlpXKDSC5QJNeV2W9j3PJsuwh/KScSj25c+WkkY0TCCS9tUPrke1DZ2wgjkdy3I6h0WnozGwkdAKGfelqZi3nzggza7UIMRc+ozOzlpv3GZ2ZtVmnM6JdTUO7vo2ZrZg7I0bNEENIAFQKFUnCS4plAKaT3T9RLpeGnkwVtpkmAsj2VcMHy0thGFmZ2dnyR2UhJJPJ+KIzvbeZRWFkCQSC5LMaJgoohp40HYOiD+YcR2dmbeYnI8xsTZh3r6uZtVnnoX43dGbWYoGY8SNgZtZmEThg2MzaTg4YXrEmWUoK4SBDDSGBYiaSRuEegKanG5XLwlJiqlD/rI5JFpXIQmeajP+QhYnMlkM3VAgTAfKwlEL9NVv+XmkoS7kWfQ896XvGk5qC/p7RSdoKfAmYBP4qIm5bsFzV8quB14EPRsSPs7KSPg38MfB8tZlbI2J3qQ4+ozOzRfrVGVENqnU7ndED9wN7JO2KiCe6VrsK2FRNW4A7gC01yn4xIv6yTj3adSFuZisWiPmoN9VwObA3Ip6JiOPAPcC2BetsA+6OjoeAM6qxouuUrcUNnZmdpDPc4VStCThb0sNd0/YFmzsfeLbr/f5qXp11lir7UUmPStop6czsO/nS1cwWWNbg1C9ExOZB1qbgDuAzdNrlzwCfB/6otLIbOjM7SdDXJyMOABd2vb+gmldnnelS2Yh47sRMSV8BvpNVwpeuZrbIXHVWt9RUwx5gk6SLJa0DrgV2LVhnF3C9Oq4AXo6Ig1nZ6h7eCe8HHssqseQZnaQLgbuBc+k09jsi4kvL7d5dUtPBW0r6HEICoFLIRxZCsn5d+bOyMJF15dCT2JAsm+5d//nCfCDdV9n9ZmUpQArZSyZmshCSZNnRmfKy48mywu8Zx46XyygJVykuWUHoSbKo/GGDO0eJUN/O6CJiVtLNwP10QkR2RsTjkm6slt8J7KYTWrKXTnjJh7Ky1aY/J+lSOj/JPuAjWT3qXLrOAp+IiB9LegPwI0kPVMtqd++a2XjodEb07xGw6gRo94J5d3a9DuCmumWr+R9YTh2WbOiqU8iD1etXJT3J4l4TM2uN9o0ZsaxvI+ki4B3AD6pZtbt3zWw8dDoj+hZHNxJqN3SSTgO+CXw8Il6h0737FuBSOmd8ny+U234ixmaGY32ospkN2hwTtaZxUSu8RNI0nUbuaxHxLajfvRsRO4AdAKfrrD6MI25mg3TiyYg2qdPrKuCrwJMR8YWu+edV9++gRveumY2PtTg4zm8BHwB+KumRat6twHXL6d5diWKWkmxwnCxLSjqYzfIzkWQhJLEhW7a+vGxjOYRkbkP5Z5vbWAovKX/nbFnT7CUTM6XwknKZySPlOIvJdUkIz5Fk2dHet0vyUJAyJVlDopSxhfx4jGJ2nuyzBndxFAEzWeaUMVSn1/X79D4umsfMmdnI6ly6rrGGzszWnmU86zoW3NCZ2UlOhJe0iRs6M1vAl65mtgZ4zAgza7VOr6uHO1yZJlkXSmWywXGyLCRZ9pIGmUjSEJJTyiEkc6eWl82eWq7H7Gnl+s+c0ntfza7PwkuKi8iuYJSMzzJRSCgydawcFjH9evnDpg4ny6bKy7JooZI09CQJL1ES8pEGg5TCdNIQjyYpT+pZkwHDZrb2+NLVzFrNva5mtia419XMWi1CzLqhM7O286XrgBQf3M/KZA/uTyR/kZKxGjSddEEWymUP52c9qzOnlz/r+OnlntVjp5e/28wbeu+T2Y3FIsyXq0hMlPsLNV/e/xOF1INTR8plpl8tL1ufJR6YXP6xM5n0GGsu6VmdS8brSJIcFHtWKR/HkfTVFv+/JN+rLt+jM7M1wQ2dmbWa4+jMbE1wHJ2ZtVoEzK61xJtmtvb40tXMWs336FZLqSs9CSFJx36YaPZQf6zrHVqQje+QPZyfhZAcPbNc/2NnlQ/C42/sHZIw+4Zk/IENyQPik0nIxFwynsTR3t9t6tXy95rbkIWQNLuU0lxh7IrZZH/MzJa3N5vsq2SZJsvLYq7w3bJwlcE90199tBs6M2u5tnVGtOuOo5mtWETnHl2dqQ5JWyU9LWmvpFt6LJekL1fLH5V02TLKfkJSSDo7q4PP6MxsATHXp15XSZPA7cB7gP3AHkm7IuKJrtWuAjZV0xbgDmDLUmUlXQi8F/jFUvXwGZ2ZLRKhWlMNlwN7I+KZiDgO3ANsW7DONuDu6HgIOEPSeTXKfhH4E5bIawpu6MxsgRPPuta8dD1b0sNd0/YFmzsfeLbr/f5qXp11imUlbQMORMRP6nwnX7qa2cki7/Bd4IWI2DzA2iwi6RTgVjqXrbUs2dBV18F3A+fSaex3RMSXJJ0FfAO4CNgHXBMRv843lmRdSMaSKGYpybKXZNvLQhWSzCaxoXcYydyGZuM7ZFlIshCSo2eXQyPmfqP3YA2nnXGkWOaNG48Wl62fKodaHJstf++Xj2zoOf/wS+U0KvNZ5pjk4kOl8Axgcqb3/p88lmSwOV6uh7LQk5nCQBlAzGaDb/T+rbPsPFE8vvvTW9rHXtcDwIVd7y+o5tVZZ7ow/58BFwM/qfbRBcCPJV0eEb/sVYk6l66zwCci4hLgCuAmSZcAtwAPRsQm4MHqvZmNuag6I+pMNewBNkm6WNI64Fpg14J1dgHXV72vVwAvR8TBUtmI+GlEvCkiLoqIi+hc0l5WauSgxhld9YEHq9evSnqSznXyNuDKarW7gP8FfLLONzez0baMS9clthOzkm4G7gcmgZ0R8bikG6vldwK7gauBvcDrwIeysk3qsax7dJIuAt4B/AA4t2oEAX5J59LWzFqgn09GRMRuOo1Z97w7u14HcFPdsj3WuWipOtRu6CSdBnwT+HhEvNJ9/yAiQlLPvwFVL8x2gA2cUvfjzGyVRLTvEbBaF9mSpuk0cl+LiG9Vs5+rYl2o/j3Uq2xE7IiIzRGxeVpJzm4zGxn9fDJiFCzZ0Klz6vZV4MmI+ELXol3ADdXrG4Bv9796ZrYaIupN46LOpetvAR8AfirpkWrercBtwL2SPgz8HLhmMFVMZAPgZCEkWYaSqXI4SEz3Xja3sVxm5pTlD2QD5SwkUA4hATjnTa/0nP/WM58vlnnzxheLy06bLIxyAxyeK5+h/+LIWT3n/2zdOcUyz3N6cdnxmXXFZZNHy/tx5kjv/T/1evk3mzySHAPJ8aHkuGK2HJZSPI7nBpyipCAQ82st8WZEfJ9ycM67+1sdMxsFY3SyVoufjDCzk7WwM8INnZkt1rJTOjd0ZraIz+jMrNUCmJ93Q2dmbRaAz+gGpJTVBIrd71l2hyx7SZr1JAlLmS+El8xPl7c3uz5ZVk7kkQ5mk2UiKYWRbH7jvmKZS9YvTCbx/50+Uc5s8sp87wwlAE+sW5hybGlHkqwhr72WZIh5pRzyUdr/2W9W+p0BJrKwpabZdArlIgufGnAQ2zjFyNUxOg2dmY0ON3Rm1m6106SPDTd0ZraYz+jMrNUCwr2uZtZ+bujGX4Me3qxc3oNX3tx8krUqNpQf6M7GeCg9oJ/1rF6yrjzUxxkT5UPkpfly72/JoY3lB/f3beydCADg8IZyLsP59eVe0tL+z36zfh8fY8mXrmbWem7ozKzVHDBsZmuBA4bNrP3c62pmbdd7qKvx5YbOzE4WuDNipKRhANkD1uVlkS5rUCZ7LnsiOZomy8vWT5XHHyiN8ZA9nJ+FkJw2UX5wH8rbLH1eNgZF9r2y/ZHtx9L+b/I7L1UuTTLR9FhdFXJnhJmtAT6jM7PWK2cJG0tu6MzsZC2Mo2vX4I1m1heKelOtbUlbJT0taa+kW3osl6QvV8sflXTZUmUlfaZa9xFJfy/pn2R1cENnZotFzWkJkiaB24GrgEuA6yRdsmC1q4BN1bQduKNG2b+IiN+MiEuB7wD/MauHGzozG6TLgb0R8UxEHAfuAbYtWGcbcHd0PAScIem8rGxEvNJV/lSWaHaXvEcnaSfwO8ChiHhbNe/TwB8DJwYouDUidi+1rb6bT75btix5vkXpsgZlkpu6yqLP58rLjs2Wf7bDc71TomTjO+RZSMohJC/Nl8NBXpl/Q8/5pfpB/r2y/ZHtRxWOgya/81Ll0uemmh6rq6SPAcPnA892vd8PbKmxzvlLlZX0WeB64GXgXVkl6pzR/TWwtcf8L0bEpdU0/EbOzAYj6DwCVmeCsyU93DVtH1o1I/4sIi4EvgbcnK275BldRHxP0kX9qZqZjYX6Z3QvRMTmZPkB4MKu9xdU8+qsM12jLHQaut3Ap0qVWMk9uo9WvR47JZ1ZWknS9hOt/UyUo+LNbHT0sdd1D7BJ0sWS1gHXArsWrLMLuL7qfb0CeDkiDmZlJW3qKr8NeCqrRNOG7g7gLcClwEHg86UVI2JHRGyOiM3TSlLqmtno6FOva0TM0rmsvB94Erg3Ih6XdKOkG6vVdgPPAHuBrwD/PitblblN0mOSHgXeC3wsq0ejgOGIeO7Ea0lfodO9a2Zt0cf+keoe/u4F8+7seh3ATXXLVvN/fzl1aHRGV3X9nvB+4LEm2zGz0VP3snWcUjnVCS/5OnAlnd6V/XRu+F0p6VI67f4+4CMDrGP/pV39STxIodzETHl7EzPlzU0ktyx1tDzgy8tHyqEivzjSe4CZJ9adX/6wRJb1pBRCAvDEsd6fV6of5N8r2x/Zfizt/+w36/fxMZbWWuLNiLiux+yvDqAuZjYixulsrQ4/1G9mi7mhM7NWG7P7b3W4oTOzxdzQmVnbZc9ojyNnLzGz1hudM7qsa77w5yXSDBTJn6Qsy8RcudzEzFxhfnl7U8eSZUfKXfhTr5b/Bh1+aWNx2c/WnVNcVnJo4+nFZdlgNlkmklIYyc9+Xa5f9r2mk/0xlSRfKe3/PCSo9+8MpMdHelwlx2PxOF7NUBZfuppZq7kzwszWBDd0ZtZ6bujMrM1E+3pd3dCZ2cl8j27EZL1SWe/YbHmsA82We9xU6I2bPFIuM/16ubdw+tVyr+vchvKy+enp4rLn6d2DeuR4ucy+jeUH7ddPlfdVNsZD6QH9rGd18lflOq57ubw/pl8t/6+cfr33cZD9ZqXfGfLjIzuu0uMxO45Xixs6M2s9N3Rm1na+dDWz9nNDZ2atFu51NbO1wGd0ZtZ2vke3EgFReBhZE9lDz71DC9T0Ieqkq18zSejJ0d4DEEyuK+/GqcPl8JL10+WQiZjMEsuUlx2fWddz/muvlet4eMMp5Y+azBIglOtfGuMhezg/CyFZ/2K5HutfKf+eU4cLIUFHl/87A5AcH9lxlR2PpWQAWdKK8vb61EK5oTOzVqs5Zus4cUNnZicRvnQ1szXADZ2ZtZ8bOjNrvZY1dB4zwsxOVmUvqTPVIWmrpKcl7ZV0S4/lkvTlavmjki5bqqykv5D0VLX+fZLOyOqw5BmdpJ3A7wCHIuJt1byzgG8AFwH7gGsi4tf1vnYDpfz4Sfh2GkIy2SwDhY73DjvQkSS8ZKr8tyQmy+EUGc2Vtzl5tPc2Z1/pHe4BML++vCwmknE55sv1nygMNZGN75BlIclCSNa9Uv49p17r/XvqSDmEpPQ7A3mGkvlyPdLQk1L2kkGPC5Hp00dLmgRuB94D7Af2SNoVEU90rXYVsKmatgB3AFuWKPsA8KcRMSvpPwN/CnyyVI86Z3R/DWxdMO8W4MGI2AQ8WL03s5bQfL2phsuBvRHxTEQcB+4Bti1YZxtwd3Q8BJwh6bysbET8fUSc+KvzEHBBVoklG7qI+B7wYo+K3VW9vgt431LbMbPx0cdL1/OBZ7ve76/m1VmnTlmAPwL+Z1aJpp0R50bEwer1L4FzSytK2g5sB9hAEoFvZqNheQHDZ0t6uOv9jojY0fc6FUj6M2AW+Fq23op7XSMipHLbXn3pHQCn66yW9eWYtVT9/6kvRMTmZPkB4MKu9xdU8+qsM52VlfRBOv0H7470ebnmva7PVdfQVP8eargdMxsxJ56M6NOl6x5gk6SLJa0DrgV2LVhnF3B91ft6BfBydcVYLCtpK/AnwO9GxOtLVaLpGd0u4AbgturfbzfcjpmNIPWpx7fqFb0ZuB+YBHZGxOOSbqyW3wnsBq4G9gKvAx/Kylab/m/AeuABSQAPRcSNpXrUCS/5OnAlnWvx/cCn6DRw90r6MPBz4Jrlff3FSllNAFSIfsjOVpUNOJJloFCSkWOyd0V0tBBLAaRJSBKaK3+3yZlyOMjMkd4fOLs+G2ynXI9I6p/9Z5goRGhMHVv+QDZQzkIC5RASgMnXev822W+WHh8zzTKbZAPgLHHV1btMad/3o33q80P9EbGbTmPWPe/OrtcB3FS3bDX/ny+nDks2dBFxXWHRu5fzQWY2Pvysq5m1nxs6M2s7n9GZWfu5oTOzVvMoYGbWds4w3A+lQT2UxDGUyswnZeaSTBLlUmg2GbDm2PHeZZLtZSaTv5oTs+WFk8eSbCmv9w49mU8G4smWpeE2SVjExEzvZaX5AJNHyr9ZOphNlomkEEaio71/Syj/zkCavSSSYy47HotZSrIBdQatQcjLKPMZnZkt4jM6M2s3jwJmZmuBOyPMrPXc0JlZuwXujDCz9nNnxCooZWrQRJYRollYRDqojnqHFqThKtmy5LMiyYSh4+V0I5NHSuEl5YwnTGThJeVi6X+Gwm82MVMOs1C27GjDwWwK+7F5CElyTZcdV9kZUiGMJMvoM3Bu6MyszRwwbGbtF9G3xJujwg2dmS3WrnbODZ2ZLeZLVzNrt6D8/O2YGp2GLnuAOXvgvyT5oYKkdy/ZZGmLSuoeyTLNlXtPNZv1Tia9glO9e1cnssErJsrLmj7UXxwjIevVTr5zOh5D1ktaGuOhac9qliyiyYP7TQ36gf92tXMj1NCZ2cjwpauZtZ57Xc2s3Zy9xMzarhMw3K6Wzg2dmS3m7CVm1nY+o+siaR/wKjAHzEbE5iULlXZgEsZQfui5HBaRPfCfjTXRJPQkSqEU5Dd10we90/CS5AH3qcJPmu3f5KF+pb9LFl5SGgch+84NQz7mlx+WkoaCpEkfGoaQZCFIozZmhO/R9fSuiHihD9sxs5HQvmddG0TimlnrRdSbapC0VdLTkvZKuqXHckn6crX8UUmXLVVW0h9IelzSvKQlryRX2tAF8A+SfiRpe68VJG2X9LCkh2foPfScmY2QagDrOtNSJE0CtwNXAZcA10m6ZMFqVwGbqmk7cEeNso8Bvwd8r85XWuml6zsj4oCkNwEPSHoqIk764IjYAewAOF1ntet82Kyt+tcZcTmwNyKeAZB0D7ANeKJrnW3A3dG5af2QpDMknQdcVCobEU9W82pVYkVndBFxoPr3EHBf9aXMbNxFzQnOPnHFVk0Lr+zOB57ter+/mldnnTpla2l8RifpVGAiIl6tXr8X+POm2zOz0aEkkmCBF2pFW6yylVy6ngvcV506TgF/ExF/15da1ZV22Q8v9CQ7fU4vAJLLA00mWTJmkxPxUohGkwwwkIaeNMrIkYVMpOMxJL91GnpSyqKS7N8mYTPQLIRkiXKrIuhnwPAB4MKu9xdU8+qsM12jbC2NG7rquvntTcub2WgS0c+A4T3AJkkX02mkrgX+cME6u4Cbq3twW4CXI+KgpOdrlK3FT0aY2WJ9augiYlbSzcD9wCSwMyIel3RjtfxOYDdwNbAXeB34UFYWQNL7gf8KnAP8D0mPRMRvl+rhhs7MFuvjI2ARsZtOY9Y9786u1wHcVLdsNf8+Oh2gtbihM7OT9fce3UhwQ2dmiyyj13UsuKEzswXqP941Lkanoct2bM3o55O3N4DQk0JEQmShG1kISRaWMpdsM81E0rtc3Qjyxdvrb3hJHrqRhZc0K1f8vAGEeww1hGSQDVEMePurYHQaOjMbHe26cnVDZ2aLOfGmmbWfGzoza7WI/HG8MeSGzswW8xmdmbWeG7pV0GRAnXR7zUJPSppnQ8nCEZqFpZSyckQh7GToGgaiNs4oUtxgw0wjDbfZbHur1NgEzfbpCBuPhs7MhihGL3XUCrmhM7OTBe6MMLM1wPfozKz13NCZWbv5of7R0u9EAJDfhC08vJ/30pXHJlD2wHy5WJ5EoFhoRA7cIfdoNu5B7XM98m2OyG9zQtC4d3xUjXdDZ2aDMWqN7wq5oTOzBfwImJm1XUA4js7MWs9PRphZ6/kenZm1WoR7XbtJ2gp8ic7gsn8VEbf1pVb90PQvUhaW0uS+RRIK0jT0IU0iUJKEq4wDh4kMWZu+Cyto6CRNArcD7wH2A3sk7YqIJ/pVOTNbDUEUsuCMq5Xk7rkc2BsRz0TEceAeYFt/qmVmq+ZEmqY605hYyaXr+cCzXe/3A1sWriRpO7AdYAOnrODjzGxoHF6yPBGxA9gBcLrOGp8/AWZrVDCAe6KrbCUN3QHgwq73F1TzzGychRNvdtsDbJJ0MZ0G7lrgD/tSKzNbVW3rjFCai3+pwtLVwH+hE16yMyI+u8T6zwM/r96eDbzQ+MP7x/U4metxsnGrxz+NiHNW8kGS/q76vDpeiIitK/m8YVhRQ7eiD5YejojNq/Lhrofr4XqsKSMyNJSZ2eC4oTOz1lvNhm7HKn52N9fjZK7HyVyPFli1e3RmZsPiS1cza71VaegkbZX0tKS9km5ZjTpU9dgn6aeSHpH08BA/d6ekQ5Ie65p3lqQHJP1j9e+Zq1SPT0s6UO2TR6oQokHX40JJ35X0hKTHJX2smj/UfZLUY6j7RNIGST+U9JOqHv+pmj/0Y6Qthn7pWmU9+RldWU+A61Yj64mkfcDmiBhqnJSkfw0cBu6OiLdV8z4HvBgRt1WN/5kR8clVqMengcMR8ZeD/OwF9TgPOC8ifizpDcCPgPcBH2SI+ySpxzUMcZ9IEnBqRByWNA18H/gY8HsM+Rhpi9U4o1vzWU8i4nvAiwtmbwPuql7fRec/2GrUY+gi4mBE/Lh6/SrwJJ2kEUPdJ0k9hio6Dldvp6spWIVjpC1Wo6HrlfVk6AdTJYB/kPSjKsvKajo3Ig5Wr38JnLuKdfmopEerS9uhXh5Jugh4B/ADVnGfLKgHDHmfSJqU9AhwCHggIlZ1f4y7td4Z8c6IuBS4CripupRbddG5n7Ba3eF3AG8BLgUOAp8f1gdLOg34JvDxiHile9kw90mPegx9n0TEXHVsXgBcLultC5av5jEydlajoRuZrCcRcaD69xBwH53L6tXyXHWP6MS9okOrUYmIeK76TzYPfIUh7ZPqXtQ3ga9FxLeq2UPfJ73qsVr7pPrsl4DvAlsZkWNkHK1GQ/f/sp5IWkcn68muYVdC0qnVDWcknQq8F3gsLzVQu4Abqtc3AN9ejUqc+I9UeT9D2CfVzfevAk9GxBe6Fg11n5TqMex9IukcSWdUrzfS6bh7ihE5RsZ1EPhUAAAAp0lEQVTRqgQMLzfryYDq8BY6Z3HQSVf1N8Oqh6SvA1fSyRDxHPAp4G+Be4E308nwck1EDLSjoFCPK+lcogWwD/hI132hQdXjncD/Bn4KnEiEdiud+2ND2ydJPa5jiPtE0m/S6WyYpHMycm9E/Lmk32DIx0hb+MkIM2u9td4ZYWZrgBs6M2s9N3Rm1npu6Mys9dzQmVnruaEzs9ZzQ2dmreeGzsxa7/8CWPT2TsLo8/sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ba9bb00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(makeImage2())\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:33.809123Z",
     "start_time": "2018-04-11T17:46:33.355067Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x11854f400>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAEDCAYAAAC/Cyi3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+MZWWd5/H3p6q7bH632A62CqJrR1dYbbUDxDETXMcROiboxCUwm5UxxnYMY5xkJ1njJuqYbOJuxpkdw4ymHYmQjI7s+ovMggaJWTAZVGARRXElA6y2SAvY3TSgTdX97h/3lNyqus/3njr33J/9eSU3XXXPved56tSpp8+5z/f5fhURmJnNk4VJd8DMrG0e2Mxs7nhgM7O544HNzOaOBzYzmzse2Mxs7nhgM5tDkq6WdFDSD1rY1xsk3dXz+LWkt7bRz1GR49jM5o+k3wOOAtdGxLkt7vd04D7ghRHxZFv7bZuv2MzmUETcAjzW+5ykfyXpa5LukHSrpJc32PXbgRuneVADD2xmx5P9wPsi4rXAnwN/12AflwGfb7VXI7Bl0h0ws9GTdDLwOuB/SFp9+lnVtj8EPtrnbQci4s09+9gJ/Bvg66Pt7fA8sJkdHxaAQxGxe/2GiPgS8KUa+7gU+HJEPN1259rmW1Gz40BEHAHul/TvANT1qk3u5nJm4DYUPLCZzSVJnwf+GXiZpJ9Jehfw74F3SfoecA9wySb2dzZwJvC/2+9t+xzuYWZzx1dsZjZ3xjp5sLTlxDhhafs4m2xVqP/zSi56S+8Z5n1NZG2NU9s/FzQ7jqP4nU2Dp44d4tjyk0Md5Te/4aR49LGVWq+94+7ffD0iLhqmvVEY68B2wtJ2Lnj5u8fZZKtC/c8XJbfzsVg+x7SSvK/QFgDJPin0JWtrnLLjQfYzJ/1vcvybHvusrWlw272fHnofjz62wne+flat1y7u/MmOoRscAYd7mNkaAXToTLobQ/HAZmZrBMHTUe9WdFp5YDOzDXzFZmZzJQhWpvyzxEE8sJnZBh3mfGCTtA24he6C2S3A/4yID1d5mb4AnA08AFwaEb8aXVfnTzrzmWn0vik5UbO+pz9Xwxlk27QAVqblfGmoToDub4B/GxGvAnYDF0m6APgAcHNE7AJurr43sznQIWo9ptXAgS26jlbfbq0eQXed2TXV89cAU50q2MzqCeDpiFqPaVVrSZWkRUl3AQeBmyLi28AZEfFQ9ZJfAGeMqI9mNkZBsFLzMa1qTR5ExAqwW9J24MuSzl23PaT+i00k7QP2AWzbetqQ3TWzkYt0ocdM2NQi+Ig4BHwTuAh4uMqouZpZ82DhPfsjYk9E7FnacuKw/TWzEeuuPKj3GETSNknfkfQ9SfdI+os+r7lQ0uGeKlgfGvZnGDiwSXpudaWGpBOANwH3AtcDV1QvuwL46rCdMbNpIFZqPmooTT6ud2tE7K4e/dKUb0qdW9GdwDWSFukOhNdFxD9J+mfguiqB3YN00wZPjaYLmZuEDjQON8gWyC+X/z9M7xIKfRn3wu5iew2PVfo725L8/1x4X9Pf2SjOq2lbWN+dPGgnhCa6CR/7TT6O1MCBLSLuBl7d5/lHgTeOolNmNjndOLbaA9sOSbf3fL8/Ivb3vqC6KLoDeCnwt9Xk43qvk3Q3cAD484i4Z/M9f4ZXHpjZBp36V2yPRMSe7AX9Jh8jordC/Z3AWRFxVNJe4CvArib9XuUMuma2xuoVW0ufsT2z37WTj73PH1mNlY2IG4CtkobK8+aBzczWCMQKC7UegySTj72veZ6qYqeSzqM7Lj06zM/gW1Ez22ATt6KDlCYf/wQgIj4FvB14r6Rl4CngshiyytRMD2wjmflM/hMqpZPOZuXS2c2tvmBuRXYeFI5/499ZMpMdnXbPx0nNlgbiWCy2s6/y5OOner6+CriqlQYrMz2wmVn7ugG6s/2frgc2M9tgsxMD08YDm5mtESFWwldsZjZnOr5iM7N50p08mO2hYbZ7b2at8+TBhDUO6ciqkTepEN54QXV5W7rHTpIwZqFwQmbn6ShKSDb5u8h+rkSj45j8ztJzJ+tiFgqSHONpWwQPsNJeHNtEzPTAZmbtW115MMs8sJnZBh3PiprZPOkugvfAZmZzJBBPt7SkalI8sJnZGhE4QNfM5o0coDtJTbN0ZCEd2dR7Z6n/5blWmmWDUCcJLWlaD2Fa/qMtdTEJ6VCyLc3G0eQ4Zm0lWVcWjpXjNrKQjux8LGUFmVx2D1+xmdkcmvXJg9nuvZm1LhCdqPcYpGZdUUn6hKT7JN0t6TXD/gy+YjOzNbrl91obGlbrih6VtBX4lqQbI+K2ntdcTLd4yy7gfOCT1b+NeWAzs3U2X6ilpGZd0UuAa6vX3iZpu6SdEfFQ03Z9K2pmawTdlQd1HlR1RXse+9bvT9KipLuAg8BNfeqKvgD4ac/3P6uea8xXbGa2wSau2NqoK9q6gQObpDOBa4Ez6A7m+yPibyR9BHg38MvqpR+sagI2koZuNFAqvDKorVJIB5TDCtK2tib7O7Zcfl8pSwek4SoLx471399S8qtu+dhDOXQj/ZkXk2j3rMhOJwnBKPzcejpNt1HclJ4fTyehLA3OxyZ/E20k5YjQSNaKRsQhSat1RXsHtgPAmT3fv7B6rrE6V2zLwH+MiDslnQLcIemmattfR8RfDtMBM5su3cmDdpZUSXou8HQ1qK3WFf2v6152PfCnkv6R7qTB4WE+X4MaA1vVwEPV149L+hFD3v+a2TRrteZBnbqiNwB7gfuAJ4F3Dtvopj5jk3Q23RqB3wZ+F3ifpHcAt9O9qvvVsB0ys8nqTh60Nitap65oAFe20mCl9rAs6WTgi8CfRcQRurEmLwF2072i+3jhfftWZ0yOLT/ZQpfNbNRWWKj1mFa1elYF1n0R+IeI+BJARDwcESsR0QE+DZzX770RsT8i9kTEnqUtJ7bVbzMbkTZXHkxKnVlRAZ8BfhQRf9XzfG8A3dtYO8thZjPseCjm8rvAfwC+XwXZAXwQuFzSbrq35A8A7xmmI00yGWTT4Vk2iCy8IcvUUZqyTzNPpGEF5X4s/Obp4ra0vUOF2/2TTijv74Sl4ram9Ov+/dcTT5X7cfopxW3p8UhCWUrHPz2GSfaUtBBQss/Izu/CedXkbyKJVKktAp7uzPnAFhHfon+xn8Yxa2Y2vbq3onM+sJnZ8aettaKT4oHNzNZoM9xjUjywmdk6vhU1sznkmgdmNle6s6Iuv1dbKAnRaPnKN52yT2RFVEqZOtKQjqRoSBp2cuSJcj92nFbeduTx/vtLQgdi29bitjTzR7JPHe0fdhKPH+37PADJz6XCzwUQzzm1vM9CH1VOMpJmZEmL0WTnQaJ0zkWDq6ZWsnsw3cG3dfiKzcw28K2omc0Vz4qa2VzyrKiZzZUIseyBzczmjW9FN0GRzM4lE0rFmdR0BjNZkJz8zkp1DSDJ15/NHGYL3R89Un7flvLM3MKRcl67TuH4xqHDxfdoe3nxeVqzIZl57pTaS2o5ZD9XejyS49jJZkwL0roMyfHoLCUL67PF6cW/ieRcLM32trEInvYGtlLNlHWvuRD4KnB/9dSXIuKjw7TrKzYz26DFK7a+NVMi4ofrXndrRLylrUY9sJnZGm3GsSU1U9YPbK2a7U8IzWwkOqjWgxoFk1etq5my3usk3S3pRknnDNt/X7GZ2RoRsFw/0eTAgsnQt2ZKrzuBsyLiqKS9wFeAXZvp83q+YjOzDdqsedCvZkqviDgSEUerr28AtkraMUz/fcVmZmu0+RlbqWbKutc8D3g4IkLSeXQvuB4dpt3pWQTfQNOF7lkP0kXwhVCFLB9/tpg9C2HonHZSedvd95Z3+aIz+z6//MD/K75nMVuoT7MF4Z0n+v/cW84+q/ie5fsfLG5beOXLy9sOl49xKRQkTk2O77OSpADJwv+FY83Ox2JTWU2Pwra2JjOjvVnRUs2Us7rtxKeAtwPvlbQMPAVcFmmRiMF8xWZmG7S1CD6pmdL7mquAq1ppsOKBzczWiPDKAzObO2Jl3svvmdnxp8XP2CbCA5uZreF8bGY2fyKd/J0JAwe20up8SacDXwDOBh4ALo2IXw1ssXTr3iBDRuNjn+Wtz0JICtkWYkuSSSTJ459m6UhCOrSlHI7Q2X5y4U1JfYUkXCWS0IfsfaX2iv0DdCD5uZLjsfDiF5Xft6PQXtPaBVk4UHIeZFlNGhnxyDPrqcHrHO3V1fmvAC4ArpT0CuADwM0RsQu4ufrezGZcVJMHdR7TamDPIuKhiLiz+vpxYHV1/iXANdXLrgHeOqpOmtl4RdR7TKtNfca2bnX+GVVKEoBf0L1VNbM5cNzMiq5fna+ez1CqNV59x+8qjck+gG1L5c+bzGw6dK/GZntgq3WTXFid/7CkndX2ncDBfu+NiP0RsSci9mzdcmIbfTazEWszu8ckDBzYktX51wNXVF9fQTdnuZnNgePhM7bS6vyPAddJehfwIHBprRZLs+wLDY5S00whydR7VnVs4dix/t04VA7biCOPF7eVCq9AOUsHDAiZ+Okv+j//2nJS0pUGRUMAVpIMJKX2dP/Pi+/hnJcWNy0eOlrc1vllOcONHnms//OnJgVsTinfWXS2LhW3jTWkY4SjSiA6UzzjWcfAgW3A6vw3ttsdM5sGU3wxVotXHpjZWsfL5IGZHWei5mMASWdK+qakH0q6R9L7+7xGkj4h6b6qoMtrhu2+r9jMbIMWr9jq1BW9mG7xll3A+cAnq38b8xWbma0RQKejWo+B+yqvXOp1CXBtdN0GbF8NJWvKA5uZrRWsFigZ/NiEpK7oC4Cf9nz/MzYOfpsy1ltRRRI+kCRUKN3MNy4M03A4j6XC4TrphOJ7snCJOHS4uC0rvpJm6iiEWSw8Vc7EsXLPj8v72/2K4rb4P/cUty2e87K+z3de/Pzy/u4o7y8LjVk4KQk72d5/tUskv7Pi73mQrChOgzov2blTfs/m2+lnE03vkHR7z/f7I2L/+hcNqCvaOn/GZmYb1R/YBhZMHlRXFDgA9AZuvrB6rjHfiprZOiKi3mPgnmrUFaW7iukd1ezoBcDhngQbjfiKzcw2ai9Ct05d0RuAvcB9wJPAO4dt1AObma0VEDVmPGvtql5d0QCubKXBigc2M+tjtlcezPTA1mTWCBgwA5s1WKi9cEJ5YXRsS2oGbC8vxF7McvIntQZKC9qzmc+v//yu4raL95ZnRbP3vbkw+alXlxfjb0lqF2S1F7LF551SHYKGM+rpOdf0vJpGM75YdKYHNjMbEQ9sZjZXVgN0Z5gHNjPbYJqTSNbhgc3MNmppVnRSPLCZ2QZtLc2aFA9sZrZWzVxr08wD26glYQWxdTF5Y3lbFvpQCkfIFrNnIR2Zi/f+UXHbwu7ChuTDm86prmI2HTafuWPaeGAzs418xWZmc6dBmqVp4oHNzNZyHJuZzSPPiprZ/Jnxgc2JJs1s7gy8YpN0NfAW4GBEnFs99xHg3cAvq5d9MCJuGFUnZ1oS3qDlLEd+w+wep/XP/5/VJ8ize5RDOm684XPFbW9+fv94jyy7x+LhJ4rbmmb3iJazexwvZv1WtM4V22eBi/o8/9cRsbt6eFAzmxdBd0lVnccAkq6WdFDSDwrbL5R0WNJd1eNDbfwIA6/YIuKWqmyWmR0v2rti+yxwFXBt8ppbI+ItrbXIcJ+xva8qR3+1pGeXXiRpn6TbJd1+bPnJIZozs3Hplsoc/BgkIm4BHht5h9dpOrB9EngJsBt4CPh46YURsT8i9kTEnqUtXjJjNhOi5qOqK9rz2NegtddVF0k3Sip/CLsJjcI9IuLh1a8lfRr4pzY6Y2ZTosW6ogPcCZwVEUcl7QW+AuwaYn9Awys2STt7vn0b0PeDQTObPXVvQ9uYOY2IIxFxtPr6BmCrpB3D7rdOuMfngQvpXnL+DPgwcKGk3XTH9QeA9wzbkSai6ZR9NpwnERgqhGDo1+XwCx0tf67YOXS4vO2JcuhDFqqg1/a/kl8852XF95QKr0CSpYNySEfWXmex3Pfl+x8sN5aEzSyc1D/EBWBh+2n9d3dy+WORrABPJKElTc+rksbFitowpkSTkp4HPBwRIek8ukfx0WH3W2dW9PI+T39m2IbNbHq1FcdWuDDaCr8tlvx24L2SloGngMuqOqND8ZIqM9uopYGtcGHUu/0quuEgrfLAZmZrtfT52SR5YDOzjTywmdm80YwnmnR2DzObO2O9YgslIRrZEFt6zygyNGRZNY4t93/+iaeK74nHj5bbSkIHtpx9VnFbZ/vJxW26/+f93/PickxHlnEjzU6SvK8U1rFQ6B9AvOpfF7ctHCofx86jvyrvs3D8lZ07C0kBnm1LyfsahoKU2somBwvbWkt861tRM5srnjwws7nkgc3M5o4HNjObJ2L2Z0U9sJnZWv6MrYHS7FA2S9X27Gc285lsi8XF/s+ffkq5rR39F2EDLBwpL5DPFoTrQJL//5yX9n067ijXPNjy4hcVt3VOLS8Wz2oUlPqfzXzGPfeV97dcTjTQpP+R/J6bJELo7rP8vnTGtNjYBOsyeGAzs7njgc3M5o1vRc1s/nhgM7O5Ep4VNbN5NONXbF4Eb2YbtFXzoEbBZEn6hKT7qkpVr2mj/+O/Yite4mZHqf+2UeSEjy2bX8i88Juk5sGRx8v729I/fARg4ZUvL27r3H1vcdtiYbF4JzlW8awkfCSRvq/QXraYPQvpyI4HSdjJwiP960rEqeU6CZ3s50pOD60k5+PKSnlbQeOaHm0YX8Hki+lWpdoFnE+3tOf5wzbqKzYzW6tuTdF2CiZfAlwbXbcB29dVwWvEn7GZ2RpiU+EeOyTd3vP9/ojYv4nmXgD8tOf7n1XPPbSJfWzggc3MNtjEwDZsweSR8MBmZhuNb1b0AHBmz/cvrJ4bij9jM7ONWvqMrYbrgXdUs6MXAIcjYqjbUPAVm5mt12J2jxoFk28A9gL3AU8C72yj3YEDm6SrgbcAByPi3Oq504EvAGcDDwCXRkQ58fzqvqJZiEZp2jsNzUimyrO88OokOf47/afsY6l8GOM5pxa3LTx6pLwtC2HIsln88tH+7zmpHN7QKPPEgPeV2svqE2RZOrKQDpbLoRSd0vFPzo8mGV4AOkvl45EOFKW/iSR8pPR31Noaz/EVTA7gynZae0adM/qzwEXrnvsAcHNE7AJurr43szmhTr3HtBo4sBXiUC4Brqm+vgZ4a8v9MrMJamvlwaQ0/YztjJ4P+H4BnFF6oaR9wD6AbVvLSRfNbEq0NzEwMUPPilb3yMXDEBH7I2JPROxZ2lLOxmpmU2R8s6Ij0XRge3h12UP178H2umRmk7S68mCWb0WbDmzXA1dUX18BfLWd7pjZNFAnaj2mVZ1wj35xKB8DrpP0LuBB4NI6jYWSjAVNhtgkdETLScGWbJeLSZhIIaxDTyeZG5I+FkMRyENBOjtOLm7TI/3XG2t7+fPNThY2k8jCbRYK7cXj5eweWeGYUpYOyI9j8fhn4UBJ+E5WCGjhWLNpwuI5l/xaonAWZ6FM9TvEVN9m1jFwYEviUN7Ycl/MbEpM821mHV55YGYbeWAzs3njKzYzmz8e2MxsrrhKlZnNm01m0J1KYx3Y0uwem691kRa7aJr5I5vOL4V1ZG1pudxUJis2kvbx1FP67++kE8r7a1o0JDv+J/cP3VD2niyrRnY8sv43yAyThe+k4UBby5k/srCfUhaPJplw2svuMdsjm6/YzGyDWb9icwZdM1urxSpVAJIukvTjqnbohhRnki6UdFjSXdXjQ8P+CL5iM7MN2po8kLQI/C3wJroVqL4r6fqI+OG6l94aEW9pp1VfsZlZHy0mmjwPuC8i/iUijgH/SDef40h5YDOztYLu5EGdx2CluqHrvU7S3ZJulHTOsD+Cb0XNbIMxFkwGuBM4KyKOStoLfAXYtcl9rDE1A1sWutFElt0jDRPZmoWJFKbls7aSEAAdK8eCdJ61tfy+LBzhlP5hFmnGihGIbYX+L2ShNsmm7HhkYSINMrJkoSCR9F9PJ2E4ydVNsVhRg7+JVrJ7wGZWHgwqmDywbmhEHOn5+gZJfydpR0Q8UrsX6/hW1MzWaDnR5HeBXZJeLGkJuIxuPsdn2pOepyrAUdJ5dMel/uXWapqaKzYzmxLRXhLJiFiW9KfA14FF4OqIuEfSn1TbPwW8HXivpGXgKeCyquRAYx7YzGyjFgN0I+IGuoWRe5/7VM/XVwFXtdeiBzYz62PWVx54YDOztQKY4noGdUzNwNZkwW86u5ksVs5m3xaOlWfLOkuFGc60rWYzpmmMUNJeZ+tS+X0FTY79ILHQf14qtpX7l81upvn/F5PjWNhnupg9mfnMzo+2z8fJLoJvaT8TMjUDm5lND9+KmtncmebSenV4YDOztY6H8ntmdnzpBujO9sjmgc3MNnLNAzObN8f1FZukB4DH6VYsWB6wGLZ16cLiTjK9nk31J7UXioucG+bP7yyV37dwLAkTyeo5FMIssrCTJvUmBip1sdQ/IJIuqlAXAJodxzQ5QbKYPV2YnoV0JP2fukHEn7EB8IZhVuGb2bRpb63opPhW1Mw2mraryE0aNm1RAN+QdIekff1eIGmfpNsl3X5s+ckhmzOzkYtWU4NPxLBXbK+PiAOSfge4SdK9EXFL7wuqbJr7AU478fmz/d+A2fHieL5ii4gD1b8HgS/TLdxgZrOuxfJ7k9B4YJN0kqRTVr8G/gD4QVsdM7PJUadT6zGthrkVPQP4cpXRdwvwuYj4Wiu9qimbek9DQbLwhqzkQWHKvmmyz8YLjZOQiaJxn4Ol9rKuZz/XSvmX1ug4Jr+z9NwZQUhH6TyeWBhI0Or5Iuki4G/oZtD9+4j42LrtqrbvBZ4E/jgi7hymzcYDW0T8C/CqYRo3s+kjorVBtWbB5IvpVqXaBZwPfLL6tzEXczGzjdqrK1qnYPIlwLXRdRuwXdLOYbrvgc3MNqo/sO1YDeeqHuvDvuoUTK5bVLk2B+ia2Vqb+4xtUF3RifDAZmYbtDjjObBgcs3XbIpvRc1snZq3ofU+YxtYMLn6/h3qugA4HBEPDfMTzPQVW5Mp9IHvS7KCFPeZTPOn4QEzHt09NZqEZ2S/syyDR5aBZATn40QErZ2bNQsm30A31OM+uuEe7xy23Zke2MxsRFqMY6tRMDmAK9tr0QObmfUxdVeRm+SBzcw28sBmZnMlAlamdx1oHR7YzGwjX7GZ2dzxwDadmn742XYWhqBZNog05CBT6Mu4PwwutZdmQkl+5vR4JMexFKk57tCMmfowPgDXPDCz+RJ5ybAZ4IHNzNYKPHlgZnNolm6d+/DAZmYbeWAzs/lSe4H71PLANkFNZ+bSk27aT8gR9K/xrLT1F8AUF2qpwwObmW007f9BDuCBzczW8ZIqM5s3AeE4NjObO155YGZzZwyfsUk6HfgCcDbwAHBpRPyqz+seAB4HVoDlOsVjXPPAzNaK6M6K1nkM5wPAzRGxC7i5+r7kDRGxu25FrKEGNkkXSfqxpPskZZ2aC4ro+0h1kkdTWXGNtttqW3Y8hisasrn2EqXf80wtZB9We8VcMpcA11RfXwO8ddgdrmo8sPWUrr8YeAVwuaRXtNUxM5uUIFZWaj0YXDA5c0ZPNapfAGcUOwTfkHRH3f0P8xnbb0vXA0haLV3/wyH2aWaTtrm0RWnBZEnfAJ7XZ9N/XtNkREgqNfr6iDgg6XeAmyTdGxG3ZJ0aZmDrV5b+/PUvqkbYfQDbtp42RHNmNjYthXtExO+Xtkl6WNLOiHhI0k7gYGEfB6p/D0r6Mt2LqnRgG/nkQUTsj4g9EbFnacuJo27OzIYUQHSi1mNI1wNXVF9fAXx1/QsknSTplNWvgT8AfjBox8MMbK2XpTezKRBVosk6j+F8DHiTpJ8Av199j6TnS1qtQ3oG8C1J3wO+A/yviPjaoB0Pcyv629L1dAe0y4A/GmJ/ZjYlqomB0bYR8Sjwxj7P/5xuZXiqz/Bftdl9K81BP+jN0l7gv/NM6fr/MuD1vwQerL7dATzSuPH2uB9ruR9rzVo/XhQRzx2mIUlfq9qr45GIuGiY9kZhqIFtqIal2+sG27kf7of7MR39mBVeeWBmc8cDm5nNnUkObPsn2HYv92Mt92Mt92MGTewzNjOzUfGtqJnNHQ9sZjZ3JjKwTUu6I0kPSPq+pLsk3T7Gdq+WdFDSD3qeO13STZJ+Uv377An14yOSDlTH5K4qVnHU/ThT0jcl/VDSPZLeXz0/1mOS9GOsx0TSNknfkfS9qh9/UT0/9nNkVo39M7Yq3dH/Bd5Ed+H8d4HLI2LsWUGqzJx7ImKsAZiSfg84ClwbEedWz/034LGI+Fg12D87Iv7TBPrxEeBoRPzlKNte14+dwM6IuLNaF3gH3dxcf8wYj0nSj0sZ4zGRJOCkiDgqaSvwLeD9wB8y5nNkVk3iiu236Y4i4hiwmu7ouFGlXHls3dMjS7q3yX6MXUQ8FBF3Vl8/DvyIbvaYsR6TpB9jFV1Hq2+3Vo9gAufIrJrEwNYv3dHYT57KphPYjVDdpHvj8D5Jd1e3qmO93ZF0NvBq4NtM8Jis6weM+ZhIWpR0F91UPjdFxESPx6w53icPXh8Ru+lmAb6yujWbuOh+PjCpOJxPAi8BdgMPAR8fV8OSTga+CPxZRBzp3TbOY9KnH2M/JhGxUp2bLwTOk3Tuuu2TPEem3iQGtqlJd9SbwA5YTWA3KQ9Xn/GsftbTN+neqEXEw9UfVQf4NGM6JtVnSV8E/iEivlQ9PfZj0q8fkzomVduHgG8CFzEl58gsmMTA9tt0R5KW6KY7un7cnWiawG6EBibdG4fVP5zK2xjDMak+LP8M8KOI+KueTWM9JqV+jPuYSHqupO3V1yfQnWi7lyk5R2bBRFYebDbd0Yj68BK6V2nQzUv3uXH1Q9LngQvppoZ5GPgw8BXgOuAsuqmdLo2IkX6wX+jHhXRvuYJurcf39HyuM6p+vB64Ffg+z9SQ+iDdz7fGdkySflzOGI+JpFfSnRxYpHvxcV1EfFTScxjzOTKrvKTKzObO8T55YGZzyAObmc0dD2xmNnc8sJnZ3PHAZmZzxwObmc1Wwbi4AAAADElEQVQdD2xmNnf+P7SHpNGEAE4NAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11855d7b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(makeImage1() - makeImage2())\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:34.004803Z",
     "start_time": "2018-04-11T17:46:33.810820Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.94166507995\n",
      "0.94166681392\n"
     ]
    }
   ],
   "source": [
    "print(makeImage1().sum())\n",
    "print(makeImage2().sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T17:46:34.031317Z",
     "start_time": "2018-04-11T17:46:34.006950Z"
    }
   },
   "outputs": [],
   "source": [
    "# The timing difference here isn't very large, but by implementing\n",
    "# this in c++ and by carefully choosing the resolution and extent of\n",
    "# the interpolating function, using the Hankel transform can be much\n",
    "# faster than doing a 2D FFT (even if the latter is also implemented\n",
    "# in c++)."
   ]
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}