{
 "metadata": {
  "name": "",
  "signature": "sha256:17c0a1c35987faeca5ef26a566c86ff0c186213b20581a5b474efebecffdc9c5"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "ATI 2014 Practicum 1"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Matthew Kenworthy** November 2014"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This practicum looks at Fourier optics, and how the geometry of the telescope pupil relates to the point spread function (PSF) of the resultant image in the telescope focal plane.\n",
      "\n",
      "The pupil plane and the focal plane are related by a Fourier transform, which is a continuous function.\n",
      "\n",
      "In computers, these functions are approximated to Discrete Fourier Transforms, and these are typically implemented with the Fast Fourier Transform (FFT) algorithm. Nearly every computer language has a version of these.\n",
      "\n",
      "Part 1\n",
      "\n",
      "  * Making amplitude and phase masks\n",
      "  * Putting them into a complex number\n",
      "  * FFT to the focal plane complex amplitude\n",
      "  * Making an intensity image\n",
      "  * Comparing PSFs"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is a Python notebook - you can run Python commands interactively in cells - the cells have a green or a grey box around them."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from numpy.fft import fft2,fftshift,ifft2,ifftshift\n",
      "%pylab inline\n",
      "\n",
      "rcParams['image.cmap'] = 'gray'\n",
      "rcParams['image.interpolation'] = 'nearest'\n",
      "rcParams['image.origin'] = 'lower'\n",
      "\n",
      "# python subroutines needed for this practicum - you don't need to understand them.\n",
      "\n",
      "def complex_cmap():\n",
      "        cdict = {'red': ((0.0, 1.0, 1.0),(0.25, 0.0, 0.0),(0.5, 0.0, 0.0),(0.75, 1.0, 1.0),(1.00, 1.0, 1.0)),\n",
      "               'green': ((0.0, 0.0, 0.0),(0.25, 1.0, 1.0),(0.5, 0.0, 0.0),(0.75, 1.0, 1.0),(1.0, 0.0, 0.0)),\n",
      "                'blue': ((0.0, 0.0, 0.0),(0.25, 0.0, 0.0),(0.5, 1.0, 1.0),(0.75, 0.0, 0.0),(1.0, 0.0, 0.0))}\n",
      "        my_cmap = mpl.colors.LinearSegmentedColormap('my_colormap',cdict,256)\n",
      "        return my_cmap\n",
      "\n",
      "def displC(c,trim=0,minlog=0):\n",
      "    c2 = np.copy(c)\n",
      "    if (trim>0):\n",
      "        (nx,ny) = c.shape\n",
      "        dx = (nx-trim) / 2 + 1\n",
      "        dy = (nx-trim) / 2 + 1\n",
      "        c2 = c[dx:dx+trim,dy:dy+trim]\n",
      "        \n",
      "    my_cmap=complex_cmap()\n",
      "\n",
      "\n",
      "        \n",
      "    fig=plt.figure(figsize=(10,8))\n",
      "    axre = fig.add_subplot(221)\n",
      "    axim = fig.add_subplot(222)\n",
      "    axamp = fig.add_subplot(223)\n",
      "    axpha = fig.add_subplot(224)\n",
      "    \n",
      "    im = axre.imshow(c2.real)\n",
      "    im = axim.imshow(c2.imag)\n",
      "    im = axpha.imshow(np.angle(c2),cmap=my_cmap)\n",
      "\n",
      "    if (minlog<0):\n",
      "        intens = c2*np.conj(c2)   # calculate intensity as modulus squares\n",
      "        intensn = intens.real / np.max(intens.real)  # normalise to peak value\n",
      "        im = axamp.imshow(np.log10(intensn),vmax=0,vmin=minlog)\n",
      "        axamp.set_title('log10(intensity)')\n",
      "    else:\n",
      "        axamp.set_title('Amplitude')\n",
      "        im = axamp.imshow(np.abs(c2))\n",
      "    \n",
      "    axre.set_title('Real part')\n",
      "    axim.set_title('Imaginary part')\n",
      "    axpha.set_title('Phase')\n",
      "    \n",
      "    show()\n",
      "    return(c2)\n",
      "    \n",
      "def make_complex_amplitude(ampli, phase, f=5):\n",
      "    # pad it out by this factor\n",
      "    (nx, ny) = ampli.shape\n",
      "    bignx = nx*f + 1\n",
      "    bigny = ny*f + 1\n",
      "    big_ampli = np.zeros((bignx,bigny))\n",
      "    big_phase = np.zeros_like(big_ampli)\n",
      "\n",
      "    dx = (nx * (f-1)) / 2 + 1\n",
      "    dy = (ny * (f-1)) / 2 + 1\n",
      "\n",
      "    big_ampli[dx:dx+nx,dy:dy+ny] = ampli\n",
      "    big_phase[dx:dx+nx,dy:dy+ny] = phase\n",
      "\n",
      "    # make this a complex number\n",
      "    Cp = big_ampli * np.exp(1j*big_phase)\n",
      "\n",
      "    return(Cp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# introduction to python numbers and arrays\n",
      "data1=np.zeros((5,8))\n",
      "print data1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.]]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# we can see how big this array is:\n",
      "print 'the number of numbers in data1 is:', data1.size\n",
      "print 'the shape of the data1 array is:', data1.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "the number of numbers in data1 is: 40\n",
        "the shape of the data1 array is: (5, 8)\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# we can now change the array with simple commands:\n",
      "print data1 + 5."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# and we can assign it to a new variable:\n",
      "data2 = data1 + 5.\n",
      "print data2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# to change just a few values in the array, you can either select single pixels:\n",
      "data2[2,3] = -7\n",
      "print data2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5. -7.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]\n",
        " [ 5.  5.  5.  5.  5.  5.  5.  5.]]\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# and you can specify a region with a range like this - NOTICE!!! the range of numbers includes the \n",
      "# starting value up to 1 less than the ending value (ie, up to but not including the ending value) \n",
      "\n",
      "data2[1:4,1:6] = 10.\n",
      "print data2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  5.   5.   5.   5.   5.   5.   5.   5.]\n",
        " [  5.  10.  10.  10.  10.  10.   5.   5.]\n",
        " [  5.  10.  10.  10.  10.  10.   5.   5.]\n",
        " [  5.  10.  10.  10.  10.  10.   5.   5.]\n",
        " [  5.   5.   5.   5.   5.   5.   5.   5.]]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# if you now make two arrays that contain \n",
      "# the x-value of the array element they are in, and another array containing\n",
      "# the y-value of the array they contain, you can now make telescope pupils\n",
      "nx, ny = data2.shape\n",
      "\n",
      "(x,y) = np.mgrid[0:ny,0:nx]\n",
      "\n",
      "imshow(x)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "<matplotlib.colorbar.Colorbar instance at 0x4a0b5a8>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAANEAAAD7CAYAAAD0DXG/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFpFJREFUeJzt3W1QVNf9B/DvMtAU1CBPIrI7jQIprAi7hmTVahBTGmKC\nJoFMMDEyGEnGpkN18qKT9EXavCA2D9OqvEmaatpYk6aZzmgi0pbUjREjqKBU7SAYHAEdklVBkSAP\ne/8vGPfPAveBPbvLvfD9zOwMZM+ec8zw5dx79uwPkyRJEojIZyGTPQEio2OIiAQxRESCGCIiQQwR\nkSCGiEiUJCg7O1sCwAcfmh7Z2dmiP3JSVFSU5vGioqKEx1MjHCLAty5ee+010aE5ngHH8/XnZXQf\nbrdb08Mf46kJBZEBSTo6I8AQkSExRABWrlzJ8Tiez9xu96SOP5JJEoy0yWTS1W8F0jd//LyYTCbc\nvn1bU9u77ror4D+fvJwjQ9LTL26GiAyJISISpKcQqZ5YaGpqgt1u9zwiIyOxY8eOYMyNSJY0/B6n\n6iMYJrSx4Ha7kZiYiLq6OlgsluEOuLFAE+CvjYWenh5NbWfOnKmvjYXq6mokJSV5AkQ0WfS0xT2h\nEH388cd45plnAjUXIs30dPWj+XKuv78fiYmJOHfuHOLi4v6/A17O0QT463Kuq6tLU9vZs2fr53Lu\n4MGDuO+++7wCdMdvfvMbz9crV66c9HezST+cTiecTqff+xUJRldXFzZt2oSzZ8/CZDJh165dWLJk\nyZh2x48fx9KlS/HJJ5/gySeflO1P80pUVFSERx55BMXFxd4dcCWiCfDXSnTt2jVNbaOjo8eMV1xc\njOzsbGzcuBGDg4O4desWIiMjvdoMDQ0hNzcXERERKCkpQUFBgfx8tITo1q1b+NGPfoTW1lbMmjVr\nzD+IISKt/BWiq1evamobExPjNV53dzfsdju++eYbxdf94Q9/wA9+8AMcP34cjz32mGKINH2ydcaM\nGXC5XGMCRDRZfH2fqLW1FXFxcSgpKcHixYtRWlqK3t5erzYdHR3Yt28fNm/eDGA4tEr48XAyJLfb\nrekx2uDgIOrr6/Hzn/8c9fX1mDFjBrZt2+bVZsuWLdi2bZtn1VRbOXnshwxJ7ge7pqYGR48elX2d\n2WyG2WzG/fffDwAoLCwcE6KTJ0+iqKgIAOByuXDw4EGEhYVhzZo14/bJEJEhyYVo2bJlWLZsmef7\nt99+2+v5uXPnwmKx4Pz587j33ntRXV2NhQsXerUZeb9UUlKC/Px82QABDBEZlMjmxM6dO/Hss8+i\nv78fSUlJ2LVrF959910AwIsvvjjh/vihPAoqf+3OXb58WVPbefPm6efNViI90dMvboaIDIkhIhJk\n2FPcRHox5VailpYWf3RDpNmUCxFRsDFERIIYIiJBDBGRIO7OEQniSkQkiCEiEsQQEQliiIgE6SlE\nqh8P7+rqQmFhIdLS0mC1WnHs2LFgzItIkZ5qcauuRL/85S+xevVqfPrpp57yQkSTzTBb3N3d3fjq\nq6/w5z//ebhxaOiY+lxEk8Ewl3NaygsRTQaRy7l77rkHGRkZsNvteOCBB8Zt43Q6YbfbkZ6erlrR\nV3ElulNeqKKiAvfff7+nlNDrr7/u1W7k3ytyOBxwOByKg9L0UVtbi9raWr/3K7ISmUwmOJ1OREdH\nj/t8V1cXXnrpJfzzn/+E2WyGy+VS7E8xRFrKCwFAWVmZ1vnTNDP6l+rOnTv90q/o5ZzS6/fu3YuC\nggKYzWYAQGxsrGJfipdzI8sLARi3vBDRZBC5nDOZTPjpT3+KrKws/PGPfxzzfHNzM65du4acnBxk\nZWXhww8/VJyL6u7c6PJCu3fv1vjPJAocuYCcOHECJ0+eVHxtTU0NEhIS8N133yE3NxepqalYsWKF\n5/mBgQHU19fjiy++QG9vL5YuXYolS5YgJSVl3P5UQ5SZmYnjx4+rNSMKKrkt7sWLF2Px4sWe7997\n770xbRISEgAAcXFxeOKJJ1BXV+cVIovFgtjYWISHhyM8PBwPPvggTp8+LRsi1uImQ/L1cq63txc3\nb94EMPzXTv71r39h0aJFXm3Wrl2LI0eOYGhoCL29vaitrYXVapWdC4/9kCH5urHQ2dmJJ554AsDw\n7vOzzz6Ln/3sZ14VUFNTU5GXl4eMjAyEhISgtLRUMUR+qYDa3Nws0gVNIykpKX6pgKr1+NmSJUtY\nAZVoPHo6scAQkSExRESCGCIiQYY5xU2kV1yJiARNuRBxi5uCbcqFiCjYGCIiQQwRkSCGiEgQt7iJ\nBHElIhLEEBEJYoiIBDFERIL0FCJNHw/XUuyOKJjcbremh5yhoSHY7Xbk5+ePec7lciEvLw82mw3p\n6en44IMPFOeiaSVSK3ZHFGyiK9H27dthtVo99RZGqqiogN1uxxtvvAGXy4Uf//jHWL9+PUJDx4+L\n5kIlelo+iUTqzrW3t6OyshKbNm0at01CQgJu3LgBALhx4wZiYmJkAwRoDJFasTuiYBMJ0datW/HW\nW28hJGT8H//S0lKcPXsW8+bNQ2ZmJrZv3644F02Xc2rF7vbs2eP5OiMjAxkZGVq6pWmgsbERjY2N\nfu/X1yujzz//HHPmzIHdbofT6Ry3TXl5OWw2G5xOJy5cuIDc3FycPn0as2bNGre9phCpFbtbv379\nBP8pNF2M/qW6d+9ev/QrF6IzZ87gzJkzsq87evQo9u/fj8rKSvT19eHGjRvYsGED/vKXv3i1+fWv\nfw0ASEpKwvz589HU1ISsrKxx+1S9nNNS7I4o2OQu3xYuXIinn37a8xitvLwcbW1taG1txccff4xV\nq1Z5BQgAUlNTUV1dDWC4Tl1TUxMWLFggOxfVlUiu2B3RZPLXAVSTyQQAXsUbX331VZSUlCAzMxNu\ntxtvvvmm4s60X4o3VlZWinRB08jq1av9Urzx73//u6a2Tz31FIs3Eo1HT2+5MERkSAwRkSCGiEgQ\nQ0QkiCEiEsQaC0SCuBIRCWKIiARNuRC1tLT4oxsizaZciIiCjSEiEsQQEQniFjeRIK5ERIIYIiJB\nDBGRID2FSHPdOSI9ESmZBShXQAWAsrIypKSkIDMzEw0NDYpzYYjIkERDdKcC6p0aCyNVVlaipaUF\nzc3NeO+997B582bFuWgKkVpqiYJNpBa3WgXU/fv3o7i4GADgcDjQ1dWFzs5O2bloCpFSaokmQyAr\noHZ0dMBisXi+N5vNaG9vl52LaojUUks0GeRCc/78eRw4cMDzGG1kBVSln+fRzyktIKq7c3dSe6fA\nN5EeyAUgOTkZycnJnu+rqqq8ntdSATUxMRFtbW2e79vb25GYmCg7F8UQaalbDMCr7lxKSgpSUlKU\nuqVppLm5Gc3NzX7v19erovLycpSXlwMAvvzyS7z99ttjKqCuWbMGFRUVKCoqwrFjxzB79mzEx8fL\n9qkYIi2pBYYL8hGNZ/Qv1dErg6/8dWsxXgXU1atXo7KyEsnJyZgxYwZ2796t3IfWCqh3UvvZZ5+N\nmcSOHTt8mT9NQ2VlZX6pgPrOO+9oavvyyy/rqwIqd+dIL/S0yaU5RNnZ2cjOzg7kXIg0M2SIiPSE\nISISxBARCWKIiAQxRESCWGOBSBBXIiJBDBGRIIaISNCUC1EgTukSKZlyISIKNoaISBC3uIkEcSUi\nEqSnELHuHBmSr9V++vr64HA4YLPZYLVa8corr4xp89e//hWZmZnIyMjAT37yEzQ2NirOhSsRGZKv\nK9EPf/hDHDp0CBERERgcHMTy5ctx5MgRLF++3NNmwYIFOHz4MCIjI1FVVYUXXngBx44dk+2TISJD\nErmci4iIAAD09/djaGgI0dHRXs8vXbrU87XD4VCsOQfwco4MSqR4o9vths1mQ3x8PHJycmC1WmXH\n+dOf/qRaiIcrERmSyBZ3SEgITp06he7ubjz88MNwOp1YuXLlmHaHDh3Crl27UFNTo9ifaoj6+vqQ\nnZ2N27dvo7+/H2vXrsUbb7zh8z+AyB/kVplLly55FV5UEhkZiUcffRQnTpwYE6LGxkaUlpaiqqoK\nUVFRiv2ohkjLjRhRsMmFyGKxeNXR/vrrr72ed7lcCA0NxezZs/H999/j3//+N1577TWvNpcuXcKT\nTz6JPXv2eFVTlaPpck7tRowo2HzdWLhy5QqKi4s9fzXiueeew0MPPeRVvPH111/H9evXPX9SJSws\nDHV1dbJ9agqR2+3G4sWLceHCBWzevFnxRowoGHwN0aJFi1BfXz/mv7/44ouer99//328//77mvvU\nFCK1G7Ha2lrP14mJiTCbzZonQFNbe3s7Ojo6/N6vnk4sTGh3Tu5GzOFw+HteNEWYzWavX6pKl0UT\noacQqb5P5HK50NXVBQCeGzG73R7wiREpEflLef6muhLJ3YgRTSY9rUSqIZK7ESOaTIYKEZEeMURE\nghgiIkEMEZEg1lggEsSViEgQQ0QkiCEiEsQQEQmaciFqaWnxRzdEmk25EBEFG7e4iQTpaSViySwy\nJF9LZrW1tSEnJwcLFy5Eeno6duzYITvG8ePHERoain/84x+Kc+FKRIbk60oUFhaG3//+97DZbOjp\n6cF9992H3NxcpKWlebUbGhrCr371K+Tl5amOxZWIDMnXlWju3Lmw2WwAgJkzZyItLQ2XL18e027n\nzp0oLCxEXFyc6lwYIjIkkQqod1y8eBENDQ1jyht0dHRg3759nmo/JpNJsR9ezpEhyQWks7MT3377\nrerre3p6UFhYiO3bt2PmzJlez23ZsgXbtm2DyWTSFEaGiAxJbos7Li7O6xLs7NmzY9oMDAygoKAA\n69evx+OPPz7m+ZMnT6KoqAjAcI2RgwcPIiwsDGvWrBl3TIaIDMnXjQVJkvD888/DarViy5Yt47b5\n5ptvPF+XlJQgPz9fNkCAhhC1tbVhw4YN+Pbbb2EymfDCCy+grKzMh+kT+Y+vIaqpqcGePXuQkZHh\nqVpVXl6OS5cuAfAu4qiVaoi0bgkSBZOvIVq+fPmETjvs3r1btY3q7pzWLUGiYPLH7py/TOieSG5L\nkCjY9HTsR3OIlLYEm5ubPV9HR0cjJibGfzMkQ7t69SquXbvm934NFyK1LcGUlBS/T4ymhpiYGK9f\nqv762IyhTnFr2RIkCjY9rUSqGwt3tgQPHToEu90Ou92OqqqqYMyNSJahNhYmuiVIFAx6Wol4YoEM\niSEiEsQQEQliiIgE6ek+nSEiQ+JKRCSIISISxBARCWKIiARNuRCNPMVNFAy+7s5t3LgRBw4cwJw5\nc/Df//533DZOpxNbt27FwMAAYmNj4XQ6FftkySwyJF/PzpWUlCie/ezq6sJLL72Ezz77DGfOnMGn\nn36qOheGiAzJ1xCtWLECUVFRsv3u3bsXBQUFMJvNAIDY2FjVuTBEZEiBOsXd3NyMa9euIScnB1lZ\nWfjwww9VX8ONBTKkQG0sDAwMoL6+Hl988QV6e3uxdOlSLFmyRPGDpwwRGZJciLq7u9Hd3e1zvxaL\nBbGxsQgPD0d4eDgefPBBnD59WjFEvJwjQ5K7fLv77rthsVg8j4lau3Ytjhw5gqGhIfT29qK2thZW\nq1XxNVyJyJB83eJet24dvvzyS7hcLlgsFvz2t7/FwMAAgOHCjampqcjLy0NGRgZCQkJQWlqqGiKT\nJHhxaTKZkJycLNIFTSMtLS3C9zMmkwkPPPCAprZ1dXUBf2OWKxEZkp5OLKjeE23cuBHx8fFYtGhR\nMOZDpImeCpWohkjtHV6iyWCoEKm9w0s0GfQUIt4TkSHp6Z7ILyG6evWq5+vw8HBERET4o1uaAnp7\ne/H999/7vd8pV2OBBexJTkREhNcv1evXr/ul3ym3EhEFm55CpLqxsG7dOixbtgznz5+HxWLR9JfD\niALNUBsLH330UTDmQTQhelqJeDlHhsQQEQliiIgETbktbqJg40pEJIghIhLEEBEJ0lOIWGOBDEnk\nzdaqqiqkpqYiJSUFv/vd78Y873K5kJeXB5vNhvT0dHzwwQeKc2GIyJB8DdHQ0BB+8YtfoKqqCufO\nncNHH32E//3vf15tKioqYLfbcerUKTidTrz88ssYHByUnYtfLudaWlr80Q2RZr5ucdfV1SE5ORn3\n3HMPAKCoqAj79u1DWlqap01CQgIaGxsBADdu3EBMTAxCQ+WjwnsiMiRf74k6Ojq8SmmZzWbU1tZ6\ntSktLcWqVaswb9483Lx5E5988olinwwRGZJciPr6+nD79m3Z15lMJtW+y8vLYbPZ4HQ6ceHCBeTm\n5uL06dOYNWvWuO15T0SGJHcPdNddd+Huu+/2PEZLTExEW1ub5/u2tjZP8fo7jh49iqeeegoAkJSU\nhPnz56OpqUl2LgwRGZKvGwtZWVlobm7GxYsX0d/fj7/97W9Ys2aNV5vU1FRUV1cDADo7O9HU1IQF\nCxbIzoWXc2RIvt4ThYaGoqKiAg8//DCGhobw/PPPIy0tDe+++y6A4Sqor776KkpKSpCZmQm32403\n33wT0dHRsn36pQIq0UT4owJqQkKCprZXrlxhBVSi8fAUN5EgQx37UTsiQTQZ9FRjQTFEWo5IEE0G\nw4Ro5BGJsLAwzxEJoslmmBCNd0Sio6Mj4JMiUqOnECluLHD7mvTKMLtzWo5IEE0Gw+zOaTkiQTQZ\nDHM5J3dEgmiy6Wkl4rEfCjp/HPuZOXOmprY9PT089kM0Hj2tRAwRGRJDRCTIMFvcRHrFlYhIkJ5C\nxI+HkyEFsngjAJSVlSElJQWZmZloaGhQnYwQAHzwMaGHP37mQkNDNT1Gjzc4OCglJSVJra2tUn9/\nv5SZmSmdO3fOq82BAwekRx55RJIkSTp27JjkcDgU58OViAxJ8nEl0vLJhP3796O4uBgA4HA40NXV\nhc7OTtm5MERkSL6GSMsnE8Zr097eLjsXbiyQIfm6xa31hM3oACq9Tnglys7OFu2CppFg/7yMPh6k\n5ZMJo9u0t7cjMTFRfhCf7+6IDGhgYEBasGCB1NraKt2+fVt1Y+Hrr79W3Vjg5RxNK1qKN65evRqV\nlZVITk7GjBkzsHv3bsU+hU9xE013k7I7F8wyXBs3bkR8fDwWLVoU0HHuaGtrQ05ODhYuXIj09HTs\n2LEjoOP19fXB4XDAZrPBarXilVdeCeh4wHAVKLvdjvz8/ICPZQgBuvSUpeXNLn86fPiwVF9fL6Wn\npwdsjJGuXLkiNTQ0SJIkSTdv3pTuvffegP77JEmSbt26JUnS8PW+w+GQvvrqq4CO984770jPPPOM\nlJ+fH9BxjCLoK1Gwy3CtWLECUVFRAet/tLlz58JmswEY3hlKS0vD5cuXAzpmREQEAKC/vx9DQ0OK\nxddFtbe3o7KyEps2bdLV+bXJFPQQTacyXBcvXkRDQwMcDkdAx3G73bDZbIiPj0dOTg6sVmvAxtq6\ndSveeusthITwffo7gv5/Yrp8nLynpweFhYXYvn275o8y+yokJASnTp1Ce3s7Dh8+DKfTGZBxPv/8\nc8yZMwd2u52r0AhBD9F0KMM1MDCAgoICrF+/Ho8//njQxo2MjMSjjz6KEydOBKT/o0ePYv/+/Zg/\nfz7WrVuH//znP9iwYUNAxjKUYN+EaXmzy99aW1uDtrHgdrul5557TtqyZUtQxvvuu++k69evS5Ik\nSb29vdKKFSuk6urqgI/rdDqlxx57LODjGEHQV6KRb3ZZrVY8/fTTAS3DtW7dOixbtgznz5+HxWJR\nfeNMVE1NDfbs2YNDhw7BbrfDbrejqqoqYONduXIFq1atgs1mg8PhQH5+Ph566KGAjTfSdLk0V8M3\nW4kEcYuFSBBDRCSIISISxBARCWKIiAQxRESCGCIiQQwRkaD/A+Pusqo/fakyAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x42ea610>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size = 51\n",
      "(x,y) = np.mgrid[0:size,0:size]\n",
      "\n",
      "# you can make masks in the array and make telescope pupils\n",
      "\n",
      "# To make a circle centred at (20,15) with a radius of 6.5 pixels, first make an array where every \n",
      "# element is the radial distance (r) from the centre (20,15).  Then, a circle of radius 6.5 can be \n",
      "# defined by selecting elements of r which are less than 6.5 (as shown below)\n",
      "\n",
      "xc = 20\n",
      "yc = 15\n",
      "\n",
      "x1 = x - xc\n",
      "y1 = y - yc\n",
      "r = np.sqrt((x1*x1)+(y1*y1))\n",
      "\n",
      "imshow(r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "<matplotlib.image.AxesImage at 0x47e9890>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAD+CAYAAADVndu7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGlBJREFUeJztnV9sFNcVxr8FLEVVaEja7NgqkUxoHGJw8Aoa85DKUHfN\nQ4prhEqDKmeV5LWVyEtNVFWhfYBtUzUlVV7aUHVppbT0hbopRQZhWw3IOKJ2HRIkRyRUDrJXoa7T\nhgRh7OlDxcoez/qePdyZHXq/n2TJO3vvzNk/Pr7nm3POTfm+74MQ4iTLqm0AIaR60AEQ4jB0AIQ4\nDB0AIQ5DB0CIw9ABEOIwKySD6uvr8dnPfhbLly9HTU0NhoaGMDU1hW9+85v4xz/+gfr6ehw9ehSr\nVq2K2l5CiEVEK4BUKoX+/n4MDw9jaGgIAJDP55HNZjE2Noa2tjbk8/lIDSWE2EccAgTzhXp6epDL\n5QAAuVwOx44ds2sZISR6fAFr1qzxm5ub/U2bNvm/+MUvfN/3/VWrVpWen5ubW/D4Fk1NTT4A/vCH\nPzH/tLa2Sv60fZEGcObMGdTV1eHDDz9ENpvFunXrFjyfSqWQSqUWzXvrrbfQ29uLI0eO4KmnnsLs\n7OyiMXNzc0s+lowJmxNcsYRlPAfn3Rrzhz/8Ad/4xjcWjZee1wbB9zPs/QUW2ho2ZtmyhYs8yXmD\nc4KPJWPC5ixfvrz0XVhqTCWPbY1ZsWLxn8Ly5cvx85//HN/5zndCx5SbYxpjOo9mzq1j+/fvx/79\n+8t+X4KIQoC6ujoAwP3334+dO3diaGgInudhcnISADAxMYF0Oh0698iRI/j73/+OI0eOYHR0VGQU\nIaQyBgYGsH//fvT392P//v3iecYVwCeffILZ2VmsXLkS165dQ29vL1544QV0dHSgUCigu7sbhUIB\nnZ2dofP37NmDubk57NmzB7Ozs5iZmVnwfFwrAMl55x+fmZkJ/e9erdqpMI+eSqVKtpbD9J85zhXA\n3NzcAnvLjanksXSMZuUWtFdCXN+P4Of2+OOPY+vWrZibm8P3v/99/OAHPxCdx+gAisUidu7cCQC4\nefMmvvWtb6G9vR2bN2/G7t27cfjw4dJtwHI0NTWJjEkKjzzySLVNEHMn2Qrced+FL33pS9U2oSJa\nW1srGp+Kshw4lUrhz3/+c+lxNTWASlYA5c5R7lgclFsBmEjSCiAYH8epAQRjZkmsHrSvpqbGOCcu\nDSBoS3DMihUrRN9VZgIS4jB0AIQ4jOg24O0wX0CJMwQod4tvqWOmkCBsTrVvA84nbEltOo9meS8Z\nIxH4wpbqUd3yNY0JmxO2FDdh4/ugDfc0cAVAiMPQARDiMHQAhDhM4jQAjU4gSfPVxI62bgMG52ji\nOUlcKBkjidU1OoHkFp/kMwmex5aWYxoTdlvNBpp4XvL9CBsT9rdjgisAQhyGDoAQh6EDIMRh6AAI\ncZjEiYBhAk5wXpJqAaqZCFTNWgCTWKcVAU3JQrY+k7AkpCRj6/sQhCsAQhyGDoAQh6EDIMRhYtUA\nNPG9ZIy2GMg0RpMIpNUENMkhmv5+toqBTDG/pNAnTBMKFuDYKPSRjqkWks/x5s2bxjHUAAghFUEH\nQIjD0AEQ4jCRawDzY5ewmC8uDSCqYqAk5QHYKgayoQFICn0k9+I1hT5hSJq9RIGkaMfG51hungmu\nAAhxGDoAQhyGDoAQh6EDIMRhql4MVE0R0MVEoGqKgFF9JkFsdfjVoOnmHJXgK4ErAEIchg6AEIeh\nAyDEYf4vNICw82oagmg0gKgSTCSNO2xoAGFzNB1+NRqArR18TOewRVRFO5omLZIxErgCIMRh6AAI\ncRg6AEIchg6AEIdJXDWgRijUbjlmmnOnJQJJRCAbAl/YGFvde21UW2qEWU3VnuQ8WvE2CEVAQoh1\nRA5gdnYWmUwGO3bsAABMTU0hm82ioaEB7e3tmJ6ejtRIQkg0iBzAoUOH0NjYWFp25PN5ZLNZjI2N\noa2tDfl8PlIjCSHRYNQAPvjgAxw/fhzf+9738NOf/hQA0NPTg4GBAQBALpfD1q1byzoBUyJQMHHC\nViKQrYKhIEnqCBREUtijie8lY4KfQVjxjY1EIAlh19Z04QliowuzZEzY+x20XzJGgnEF8Nxzz+HF\nF19ccMFisQjP8wAAnuehWCxWfGFCSPVZ0gG8/vrrSKfTyGQyZT1xKpVSlSESQqrPkiHA2bNn0dPT\ng+PHj+P69ev497//ja6uLnieh8nJSdTW1mJiYgLpdLrsOXp7e0u/r1mzBg8++KA96wkhAIBz587h\nzTffrHheyhcGWQMDA/jJT36CP/3pT/jud7+Lz33uc+ju7kY+n8f09HSoBpBKpXDgwIHSY4kGYCsP\nIHhMszOQZkfhqJDE95p4M6wzryRXwNThNywOl3QFDs4zPbY1JmxOTU3NkmOCz0vHaM4bPGYas3bt\nWpF+UlEewK0vz759+3Dy5Ek0NDTg9OnT2LdvXyWnIYQkBHEmYGtrK1pbWwEA9913H06dOhWZUYSQ\neGAmICEOQwdAiMPEWgwUFPwAXTGQJHnIRiKQrUSV4Jioiko0HX4liVdhAlmSt1WTYCOpJ6rPJOxv\nQNPNWQJXAIQ4DB0AIQ5DB0CIw8TaFdiWBmDSBABdUo8keSjJDUHCEms0XYGDMb8kicrGduy20DTh\n0MzRFu1oioE0mpYErgAIcRg6AEIchg6AEIepelPQ+RoBEB7HmAqGtE1B4yoG0uQBBJHEkpL7x5JC\nH8n7ENcOSRKiiueDaO7Fa+7pMw+AEBILdACEOAwdACEOQwdAiMMkrhgoKAoCZkFPu+OQphjIJHZF\nlQgUJvBoCk9MO/oAsoSd4HnCOtSYzivB1k47GhEwOEaT5KNJ6gn7rgY/JyYCEUJuGzoAQhyGDoAQ\nh6m6BiDpCmxKFrKlAWiKgWwVuNgoVpEUA2mKSsLeO1PML3kf4myMYhqjSerRagCmMSwGIoTEAh0A\nIQ5DB0CIw9ABEOIwkYuA84UTbSKQpiOQjUQgjdASZyKQpMOvRgQMvoYwcVGDjSpI7TbpNkRAyftt\noxu1pCI27DOhCEgIqQg6AEIchg6AEIeJtSuwJMlHM0abCBSMqzQdgWx1xtF0irXR4VeS7BRVt58w\ne03JNpIdfCQddWzE6raKdiS2mM4hnReEKwBCHIYOgBCHoQMgxGEi1wDmxyWSe/yaMVoNIHhMExcm\nSQPQdPi1FUtKCL6msPwN0+uWaAC2CnDiKtqxMafcMRNcARDiMHQAhDjMkg7g+vXraGlpQXNzMxob\nG/H8888DAKamppDNZtHQ0ID29nZMT0/HYiwhxC5LOoC77roLfX19GBkZwejoKPr6+vDGG28gn88j\nm81ibGwMbW1tyOfzcdlLCLGIUQT8zGc+AwC4ceMGZmdnce+996KnpwcDAwMAgFwuh61bt5Z1AvOF\ntjCRQiLWmUQ/SSdhScGQptNQkKiKgcKKPyTbfEUlWgaRiHXB904i6EW1dXaSxLrg+xJM3pLMKTfP\nhFEDmJubQ3NzMzzPw7Zt27B+/XoUi0V4ngcA8DwPxWKx4gsTQqqP0WUsW7YMIyMj+Oijj7B9+3b0\n9fUteD6VSi353+v8+fOl39PpNGpra2/DXEJIGMPDwxgdHa14nnjNcM899+CJJ57A+fPn4XkeJicn\nUVtbi4mJCaTT6bLzNm3aVPo9bKlOCLl9MpkMtmzZUnr86quviuYt6QCuXr2KFStWYNWqVfj0009x\n8uRJvPDCC+jo6EChUEB3dzcKhQI6OzvLnmN+rCIpyJHE3ZoCIon+ELRFkwgk2UVH0rAi+Fhif1gM\naIr5JTv6hGGK1TXxfdixqLbOtlG0I/msNQk72iQfjZ6zpAOYmJhALpfD3Nwc5ubm0NXVhba2NmQy\nGezevRuHDx9GfX09jh49WvGFCSHVZ0kH0NTUhL/97W+Ljt933304depUZEYRQuKBmYCEOAwdACEO\nE2tX4Kiq9jSdhAGzABmWPBQUfkyJQVqCIpWm20/YMUmHX9O22ED1qhe1XXhsVEXa6BglOY+2S5Mm\nCY0rAEIchg6AEIehAyDEYe5IDUCSsKPpNCRJStIkhwSPSRJgJMlDmg6/pusCi98HW0U7Eg3AFM9L\nkm9s7QKl0Q00sbqNOeWOmeAKgBCHoQMgxGHoAAhxmFi7Akvum2p2+tUWGQXHaObYygsw3feXFPpI\nkBTtSOJ505g4d87VFMrYKOzR3q/X5BdQAyCEWIcOgBCHoQMgxGHoAAhxmFhFQEmChkaMkQiHYWNM\nIqBGkNQmAgXPExQFbYk+ms49moQXSdFOmP2a90GToGNDVJPYIpkn+Rw1nackcAVAiMPQARDiMHQA\nhDhM5BrA/FjFVhdVU4MQQNdtWKIBmM5ra3twiTYSRJLUI0nyMcXhYfZF1dRCs5OR5NoaDcBG7B7l\nGA1cARDiMHQAhDgMHQAhDkMHQIjDxJoIFKcIKBGPNF2BTXPCxEZJIpAk4SWIJqlH8r5I3t+ots2K\nq1tOteZoz2s6hxauAAhxGDoAQhyGDoAQh4lcA5iPNrbRFGVodAJNh2LJVuUSDcCU4CKJ7yUJO5It\nxYPY+kw02IqHTe9v0pN6ojovVwCEOAwdACEOQwdAiMNUvRjIRhdVbXxkauYhuUcu2YFIsstPEM1O\numH6gynGj/Mziev+dtILcGxcl3kAhJDbxugAxsfHsW3bNqxfvx4bNmzAyy+/DACYmppCNptFQ0MD\n2tvbMT09HbmxhBC7GB1ATU0NXnrpJbz99tsYHBzEK6+8gosXLyKfzyObzWJsbAxtbW3I5/Nx2EsI\nsYjRAdTW1qK5uRkAcPfdd+ORRx7BlStX0NPTg1wuBwDI5XI4duxYtJYSQqxTkQh4+fJlDA8Po6Wl\nBcViEZ7nAQA8z0OxWIzEQAnaLjE2xkgEs6AwWFNTU/F5JbZoxmg7GNkgyUKchCTbJkUsAn788cfY\ntWsXDh06hJUrVy54LpVKhWaqEUKSjWgFMDMzg127dqGrqwudnZ0A/vdff3JyErW1tZiYmEA6nQ6d\n+84775R+X7VqFT7/+c9bMJsQMp+RkRFcuHCh4nlGB+D7Pp599lk0NjZi7969peMdHR0oFAro7u5G\noVAoOYYgjY2Npd9v3LhRsYGEEDPNzc147LHHSo9//etfi+YZHcCZM2fw29/+Fo8++igymQwA4ODB\ng9i3bx92796Nw4cPo76+HkePHtVZbgFJYo0kRNGMkXTZDcb8ku22o7I3iOS9iwpbr7FaJNk2KUYH\n8Pjjj5cVik6dOmXdIEJIfDATkBCHoQMgxGHoAAhxmMirAecLJZKtq8JEqeAxW4JZsINOsJouzJbg\nHE0SksSW4GOJcBjWEcj03sT5mcQl+sUloEZFnLZxBUCIw9ABEOIwdACEOEysXYG1sY0mjpXE6qYY\nWtJJOIikc49GAwiL74PHNJ2D4/xMNNiKh00JT3HqBpr3JqrzcgVAiMPQARDiMHQAhDhM5BrA/NhL\nEqNK7jlr4vuw8wY75mqackiuI8H0GsO6+0rGmO7XS3IdJJ+J6TqSOWHzJPZqchCqNUd7XtM5tHAF\nQIjD0AEQ4jB0AIQ4DB0AIQ5T9WIgGyKgZEssyVZVNkQ/yXUkwpAkESj4GiXJQsE5cYqANgQ9bfen\nuERAiS1RjdHAFQAhDkMHQIjD0AEQ4jCxJgKFxajBJB6NThCWAKOJ74NoYvUwPUKCqUhHmwhk0gkk\nyUOSWN2k04SNkRRFxZVgZKvgSaNR2OhGLb32ojkVzyCE/N9AB0CIw9ABEOIwdACEOEysIqAkoUQi\n6AWFQ1tJPkEkQkvQFon9miQZScKORiiUJA9JkpA0CV2aMZI5YURRvRhVR2VbFY8SuAIgxGHoAAhx\nGDoAQhwmcg1gfvwY1rlHEs+bxoSdNyxuDaJJDgleW5v4Y0LSvVcSz5vGSDolaXQCrbZgI1bXnFeT\n7KSZEzaGGgAhpCrQARDiMHQAhDjMHakBBB/X1NQY7dA0fAiLxUz3/SX224oLJTkUppg/7L2TxNQ2\nNICoioxsxfOa+/U2cgVsNVORwBUAIQ5jdADPPPMMPM9DU1NT6djU1BSy2SwaGhrQ3t6O6enpSI0k\nhESD0QE8/fTTOHHixIJj+Xwe2WwWY2NjaGtrQz6fj8xAQkh0GB3Al7/8Zdx7770LjvX09CCXywEA\ncrkcjh07Fo11hJBIUYmAxWIRnucBADzPQ7FYLDvWJAJKinbC5pkICiI3b940jgkKK2FJPqbEH03R\nUZgtQbSimqYrcFAY1HQbjqrISHte03sV5/ZnpuIlWwVPEm5bBEylUpG1LCaERItqBeB5HiYnJ1Fb\nW4uJiQmk0+myY8+fP1/6PZ1Oo7a2VnNJQsgSDA8PY3R0tOJ5KgfQ0dGBQqGA7u5uFAoFdHZ2lh27\nadOm0u8zMzOayxFCDGQyGWzZsqX0+NVXXxXNMzqAPXv2YGBgAFevXsUDDzyAH/7wh9i3bx92796N\nw4cPo76+HkePHi07f37sFRY/S7b61hC8lqa5h6SLsaQwSYKNpBOJBqCJl8OShUzn0dgSdkyiWcTV\naCSqpieSwi+JFqIJxY0O4LXXXgs9furUqYovRghJFswEJMRh6AAIcZhYm4JKinZu9xrljkmKJyTx\nfJI1AE3TEO19dVOugKTISFK8ZMveKM6bJD2i3DETXAEQ4jB0AIQ4DB0AIQ5DB0CIw0QuAkYh/Gk6\nqGgSgSQiYJCoioGiSjrRdBIOOyYpIAqOsZW4FFVXYNMYbXcijS2mc0jnLZpT8QxCyP8NdACEOAwd\nACEOE2tXYG18HEQT30s6/AYLiDS7DtvSADQ7wWjiY22jEVPML0ny0RQZhZ03rp2MqrkDsqYoSgJX\nAIQ4DB0AIQ5DB0CIw9ABEOIwkYuAYaLNfCRdTIJjgmKdVowxdQDSdCiOMxFIUzEY1TbekkQgzXbm\ntroNR5FgFJUtrAYkhMQCHQAhDkMHQIjDxKoBaOJ9YPGuPlElAmk6AkkInle1jbOyI5ApjtV27zUl\npkiSfDRJPba0hbg0gGoWGUngCoAQh6EDIMRh6AAIcZiq5wEECds9yHS/Oyz2Dcbdkl1+NMVAlT5f\nDpMuYKsYyEYn4bBjmlhdogFI5mgKkWxoC9rmJCadwFZ+gQSuAAhxGDoAQhyGDoAQh6EDIMRhYu0K\nrE3Y0WybpenwKykG0nQA0iQCaToCRbUNVVyiWtgYjXCo2XZcY2/St0CXwBUAIQ5DB0CIw9ABEOIw\nVS8GshHrhiUPSZJ64ioGsoGtYiBJ8pCNWFejG4TNs5E8FHbMRuMO7WtkIhAhJBHclgM4ceIE1q1b\nh4ceegg/+tGPbNlECIkJtQOYnZ3Ft7/9bZw4cQLvvPMOXnvtNVy8eNGmbYSQiFE7gKGhIXzxi19E\nfX09ampq8OSTT+KPf/yjTdsIIRGjFgGvXLmCBx54oPR49erVOHfu3OILrFiBS5cuYe3ataIEGInY\nFRTiJCKVRNC79Xh8fHzBa5uPra3ATEjEUWChrbYSgTQioDR56NZ3YakxS53HlggorTJ866230NTU\nFDrG1tbqNjs1Dw4OYsuWLYueK4d6BVBJW6tLly5pL1MVPvjgg2qbIOZOshW4874LFy5cqLYJFTE4\nOFjRePUK4Atf+ALGx8dLj8fHx7F69epF43p7e3Hp0iX09vZizZo1ePDBB7WXJISU4dy5c3jzzTcx\nODiIn/3sZ+J56hXA5s2b8e677+Ly5cu4ceMGfv/736Ojo2PRuPb2dqxduxbt7e384yckIlpaWrB3\n715s2bIFe/fulU/0b4Pjx4/7DQ0N/tq1a/0DBw4ser61tdUHwB/+8Cfmn9bWVtHfcMqPSsUihCQe\nZgIS4jB0AIQ4TKQOIOmpws888ww8zyvd5wWAqakpZLNZNDQ0oL29HdPT01W0cCHj4+PYtm0b1q9f\njw0bNuDll18GkFybr1+/jpaWFjQ3N6OxsRHPP/88gOTaC/wvwzWTyWDHjh0Akm1rfX09Hn30UWQy\nGTz22GMAKrc3MgdwJ6QKP/300zhx4sSCY/l8HtlsFmNjY2hra0M+n6+SdYupqanBSy+9hLfffhuD\ng4N45ZVXcPHixcTafNddd6Gvrw8jIyMYHR1FX18f3njjjcTaCwCHDh1CY2NjKc8lybamUin09/dj\neHgYQ0NDABT23s5dgKU4e/asv3379tLjgwcP+gcPHozqcmref/99f8OGDaXHDz/8sD85Oen7vu9P\nTEz4Dz/8cLVMM/L1r3/dP3ny5B1h87Vr1/zNmzf7Fy5cSKy94+Pjfltbm3/69Gn/a1/7mu/7yf4+\n1NfX+1evXl1wrFJ7I1sBhKUKX7lyJarLWaNYLMLzPACA53koFotVtiicy5cvY3h4GC0tLYm2eW5u\nDs3NzfA8rxS+JNXe5557Di+++OKCtOak2gr8bwXw1a9+FZs3b8Yvf/lLAJXbG1lDEM0OuEkjlUol\n8nV8/PHH2LVrFw4dOoSVK1cueC5pNi9btgwjIyP46KOPsH37dvT19S14Pin2vv7660in08hkMujv\n7w8dkxRbb3HmzBnU1dXhww8/RDabxbp16xY8L7E3shWANFU4aXieh8nJSQDAxMQE0ul0lS1ayMzM\nDHbt2oWuri50dnYCSL7NAHDPPffgiSeewPnz5xNp79mzZ9HT04M1a9Zgz549OH36NLq6uhJp6y3q\n6uoAAPfffz927tyJoaGhiu2NzAFIU4WTRkdHBwqFAgCgUCiU/siSgO/7ePbZZ9HY2Lgg3TOpNl+9\nerWkQn/66ac4efIkMplMIu09cOAAxsfH8f777+N3v/sdvvKVr+A3v/lNIm0FgE8++QT/+c9/AADX\nrl1Db28vmpqaKrc3KoHC982pwtXmySef9Ovq6vyamhp/9erV/q9+9Sv/n//8p9/W1uY/9NBDfjab\n9f/1r39V28wSf/3rX/1UKuVv3LjRb25u9pubm/2//OUvibV5dHTUz2Qy/saNG/2mpib/xz/+se/7\nfmLtvUV/f7+/Y8cO3/eTa+t7773nb9y40d+4caO/fv360t9XpfYyFZgQh2EmICEOQwdAiMPQARDi\nMHQAhDgMHQAhDkMHQIjD0AEQ4jB0AIQ4zH8BjRFwif3EYjMAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x480a910>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# now make a mask to select the points within a radius of 6.5 and make the amplitude 2.0\n",
      "mask1 = (r < 6.5) * 2.0\n",
      "imshow(mask1)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<matplotlib.colorbar.Colorbar instance at 0x456e7e8>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAASsAAAEACAYAAADrz1BBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGklJREFUeJzt3X9M1Pf9B/DnOXB1YhR/cCKQnQooCNydo560aVelh9RO\nRDFOspELOmqci6FpFmubZdE/8GazGS1dootYpi2WLJ0yh2yigFOKOASxlQC2MH4It17Pqz9oBeHz\n/cN43576uftwcOPzxucjuaR39773vdDy9P353OfeL40kSRKIiFRuwlgXQESkBMOKiITAsCIiITCs\niEgIDCsiEgLDioiEwLAiolHX2dmJZcuWYdGiRYiLi8P+/fufOG7btm2IioqCXq9HfX29xzkD/FEo\nET3dAgMDsXfvXhgMBty5cwc/+tGPYDabERMT4xpTWlqK69evo7W1FRcvXsSWLVtQU1MjO6eilZVO\np0NCQgKMRiOWLFkCAHA4HDCbzYiOjkZKSgqcTucIfzwiGi9mz54Ng8EAAAgKCkJMTAxu3LjhNqak\npAQWiwUAYDKZ4HQ6YbPZZOdUFFYajQaVlZWor69HbW0tAMBqtcJsNqOlpQXJycmwWq0+/VBENL61\nt7ejvr4eJpPJ7fHu7m5ERES47oeHh6Orq0t2HsXnrB79Vs53U9FiseD48eNKpyKip8SdO3ewbt06\n7Nu3D0FBQY89/2iuaDQa2bkUnbPSaDR4+eWX8b3vfQ+bN29GTk4ObDYbtFotAECr1T5x+ebpjYnI\nv0b6td/p06fj5s2bisYGBQXh9u3bbo8NDAwgIyMDP//5z5Genv7Ya8LCwtDZ2em639XVhbCwMPk3\nkRS4ceOGJEmS9N///lfS6/XSuXPnpGnTprmNCQ4Ofux1AHjjjbcxuo0UAGloaEjR7dH3GxoakrKy\nsqTc3FzZ+f/+979Lr7zyiiRJkvTJJ59IJpPJYz2KVlahoaEAgFmzZmHNmjWora2FVqtFb28vZs+e\njZ6eHoSEhCiZiogEIvm4Ortw4QKOHj3q+mAOAPLy8tDR0QEA2Lx5M1auXInS0lJERkZi8uTJOHz4\nsMc5NZKXavr6+jA4OIgpU6bg7t27SElJwW9/+1uUl5djxowZ2L59O6xWK5xO52Mn2XkYSDR2fA2a\nhzQaDe7fv69obEBAwIjfz2s93sKqra0Na9asAQDcv38fP/vZz7Bjxw44HA6sX78eHR0d0Ol0KC4u\nxrRp09wnZ1gRjZnRCKv+/n5FYydOnDj2YTWiyRlWRGNmNMLq3r17isZ+//vf93tY8Qp2IpLl7wAa\nDoYVEcliWBGREBhWRCQEhhURCYFhRURCGBoaGusSXBhWRCSLKysiEgLDioiEwLAiIiEwrIhICAwr\nIhKCmsKKrbiISNbQ0JCi26M2btwIrVaL+Pj4J85rt9uRmpoKg8GAuLg4vP/++15rYVgRkSxJkhTd\nHpWdnY2ysjLZefPz82E0GtHQ0IDKykq88cYbXvfOYlgRkSxfw+qFF15AcHCw7LyhoaG4desWAODW\nrVuYMWMGAgI8n5XiOSsikuWvc1Y5OTlYvnw55syZg9u3b6O4uNjraxhWRCRLLqyqq6tRXV3t87x5\neXkwGAyorKzE559/DrPZjCtXrmDKlCmyr2FYEZEsubBKSkpCUlKS6/4f/vCHYc1bXV2Nt99+GwAw\nf/58zJ07F83NzUhMTJR9Dc9ZEZEsX89ZebNw4UKUl5cDAGw2G5qbmzFv3jyPr+HKiohk+brrQmZm\nJqqqqmC32xEREYGdO3diYGAAwIM2XG+99Rays7Oh1+sxNDSEPXv2YPr06R7nZMMIonFqNBpG/Oc/\n/1E09oc//CEbRhDR2FHTFewMKyKSxbAiIiEwrIhICAwrIhIC92AnIiFwZUVEQmBYEZEQGFZEJASG\nFREJgWFFREJgWBGREHjpAhEJQU0rK+5nRUSyfN3Pylt3GwCorKyE0WhEXFwcXnrpJa+1KAqrwcFB\nGI1GrFq1CgDgcDhgNpsRHR2NlJQUOJ1OJdMQkWD81d3G6XRi69at+Nvf/oZPP/0Uf/nLX7zWoiis\n9u3bh9jYWNf+VFarFWazGS0tLUhOTobValUyDREJxl/dbT788ENkZGQgPDwcADBz5kyvtXgNq66u\nLpSWluIXv/iFq6iSkhJYLBYAgMViwfHjx72+ERGJx1/bGre2tsLhcGDZsmVITEzEkSNHvL7G6wn2\n119/He+8846rxxfwYM9krVYLANBqtbDZbMMulojUTy6I/v3vf6Ours7neQcGBnD58mWcOXMGfX19\nSEpKwtKlSxEVFSX7Go9hdfLkSYSEhMBoNKKysvKJYzQaDbcvJhqn5C5dWLx4MRYvXuy6f/DgwWHN\nGxERgZkzZ2LSpEmYNGkSXnzxRVy5csVjWHk8DKyurkZJSQnmzp2LzMxMnD17FllZWdBqtejt7QUA\n9PT0ICQkZFiFEpEY/HUYuHr1apw/fx6Dg4Po6+vDxYsXERsb6/E1HsMqLy8PnZ2daGtrw7Fjx7B8\n+XIcOXIEaWlpKCwsBAAUFhYiPT192MUSkfr5GlaZmZl47rnn0NzcjIiICBQUFODAgQM4cOAAgAet\nuFJTU5GQkACTyYScnByvYaW4u01VVRV+//vfo6SkBA6HA+vXr0dHRwd0Oh2Ki4sxbdq0xyfn4SHR\nmBmN7jY1NTWKxi5dutTvF5CyFRfRODUaYfXJJ58oGpuUlMRWXEQ0dtT0dRuGFRHJYlgRkRC46wIR\nCYErKyISAsOKiITAsCIiITCsiEgIDCsiEgLDioiEwEsXiEgIXFkRkRDUFFbsbkNEsvzZ3QYALl26\nhICAAHz88cdea2FYEZEsf3W3AR50zdq+fTtSU1MVreAYVkQky1/dbQDg3Xffxbp16zBr1ixFtfCc\nFRHJ8tengd3d3Thx4gTOnj2LS5cuKdr7jmFFRLL8dYI9NzcXVqsVGo1G8T7uDCsikiUXIlevXsXV\nq1d9nreurg4bNmwAANjtdpw6dQqBgYFIS0uTfQ3DiohkyYVVXFwc4uLiXPePHTs2rHm/+OIL139n\nZ2dj1apVHoMKYFgRkQe+HgZmZmaiqqoKdrsdERER2LlzJwYGBgAAmzdv9mlONowgGqdGo2HEX//6\nV0Vj16xZw4YRRDR21HQFO8OKiGTxi8xEJASurIhICAwrIhICw4qIhMCwIiIhMKyISAgMKyISAi9d\nICIhcGVFREJgWBGREBhWRCQENYUV92AnIln+6m7zwQcfQK/XIyEhAc8//zwaGxu91uIxrL799luY\nTCYYDAbExsZix44dAACHwwGz2Yzo6GikpKTA6XQq+bmJSDD+6m4zb948nDt3Do2NjfjNb36D1157\nzWstHsPqmWeeQUVFBRoaGtDY2IiKigqcP38eVqsVZrMZLS0tSE5OhtVqVfBjE5FohoaGFN0e5a27\nTVJSEqZOnQoAMJlM6Orq8lqL18PAH/zgBwCA/v5+DA4OIjg4GCUlJbBYLAAAi8WC48ePe30jIhKP\nryur4Th06BBWrlzpdZzXE+xDQ0NYvHgxPv/8c2zZsgWLFi2CzWaDVqsFAGi1WthsthEVS0Tq5O8T\n7BUVFSgoKMCFCxe8jvUaVhMmTEBDQwO+/vprrFixAhUVFW7PazQabl9MNE7JhVVzczNaWlpGNHdj\nYyNycnJQVlbmtSEqMIxLF6ZOnYpXX30VdXV10Gq16O3txezZs9HT04OQkJARFU1E6iQXVtHR0YiO\njnbdP3ny5LDm7ejowNq1a3H06FFERkYqeo3Hc1Z2u931Sd8333yD06dPw2g0Ii0tDYWFhQCAwsJC\npKenD6tQIhKDr+esMjMz8dxzz6G5uRkREREoKCjAgQMHcODAAQDArl27cPPmTWzZsgVGoxFLlizx\nWovH7jZXr16FxWJxnfHPysrCr3/9azgcDqxfvx4dHR3Q6XQoLi7GtGnTHp+ch4dEY2Y0utv88Y9/\nVDT2l7/85dh2t4mPj8fly5cfe3z69OkoLy/3W1FEpA7cdYGIhKCmr9swrIhIFsOKiITAsCIiITCs\niEgIDCsiEgI/DSQiIXBlRURCYFgRkRAYVkQkBIbVOOKvv0x+r5LUgGFFREJgWBGRENR06QJbcRGR\nLH+14gKAbdu2ISoqCnq9HvX19V5rYVgRkSx/teIqLS3F9evX0draioMHD2LLli1ea2FYeeDvrh5q\nfW+ih3wNK2+tuL7bIctkMsHpdHptPMOwIiJZ/mrF1d3djYiICNf98PBwr70DeYKdiGTJBVFbWxva\n29tHdW5vl+swrIhIllxY6XQ66HQ61/2qqqphzRsWFobOzk7X/a6uLoSFhXl8DQ8Dv0O0c0Ki1Uvi\n8bV9vDdpaWn485//DACoqanBtGnTXI2T5XBlRUSyfP1HMDMzE1VVVbDb7YiIiMDOnTsxMDAAANi8\neTNWrlyJ0tJSREZGYvLkyTh8+LDXORlWRCTL17AqKiryOiY/P39YczKsiEiWmk4vMKyISBbDioiE\nwLAiIiEwrIhICGradYFhRUSyuLJSCTX9RYyGJ/083HGURkJNvyNPdVgRkWcMKyISAsOKiITAsCIi\nITCsiEgIvHSBiISgppWV1/2sOjs7sWzZMixatAhxcXHYv38/AMDhcMBsNiM6OhopKSlwOp1+L5aI\n/rdGsq1xWVkZFi5ciKioKPzud7977Hm73Y7U1FQYDAbExcXh/fff91iL17AKDAzE3r178dlnn6Gm\npgbvvfcempqaYLVaYTab0dLSguTkZFitVmU/PREJw9ewGhwcxK9+9SuUlZXh2rVrKCoqQlNTk9uY\n/Px8GI1GNDQ0oLKyEm+88Qbu378vW4vXsJo9ezYMBgMAICgoCDExMeju7nbrTmGxWHD8+PFh/SEQ\nkfr5Gla1tbWIjIyETqdDYGAgNmzYgBMnTriNCQ0Nxa1btwAAt27dwowZMxAQIH9maljnrNrb21Ff\nXw+TyQSbzebahlSr1Xpto0NE4vH1nNWTutdcvHjRbUxOTg6WL1+OOXPm4Pbt2yguLvY4p+I92O/c\nuYOMjAzs27cPU6ZMcXtOo9Hwax1E45Cve7AryYO8vDwYDAbcuHEDDQ0N2Lp1K27fvi07XtHKamBg\nABkZGcjKykJ6ejqAB6up3t5ezJ49Gz09PQgJCVEyFREJRG5ldePGDfT09Mi+7tHuNZ2dnQgPD3cb\nU11djbfffhsAMH/+fMydOxfNzc1ITEx84pxeV1aSJGHTpk2IjY1Fbm6u6/G0tDQUFhYCAAoLC10h\nRkTjh9w5qtDQUCxevNh1e1RiYiJaW1vR3t6O/v5+fPTRR0hLS3Mbs3DhQpSXlwMAbDYbmpubMW/e\nPNlavK6sLly4gKNHjyIhIQFGoxEAsHv3brz55ptYv349Dh06BJ1O5/V4k4jE4+s5q4CAAOTn52PF\nihUYHBzEpk2bEBMTgwMHDgB40OHmrbfeQnZ2NvR6PYaGhrBnzx5Mnz5ddk6N5MervtR+HktNF7z5\ni9r/Dsh/Rvr/t0ajwcaNGxWNLSgo8PvvE69gJyJZavoHnWFFRLIYVkQkBH6RmYiEwJUVEQmBYUVE\nQmBYqcSjH+ur6S/GF7xMgUabmn4nnuqwIiLPGFZEJASGFREJgZcuEJEQuLIiIiEwrIhICGoKK8U7\nhRLR08ef3W0AoLKyEkajEXFxcXjppZc81sKVFRHJ8nVl9bC7TXl5OcLCwvDss88iLS0NMTExrjFO\npxNbt27FP/7xD4SHh8Nut3uckyur73i4l7woe8qLVi+Jx5/dbT788ENkZGS4tjueOXOmx1oYVkQk\ny9eGEU/qbtPd3e02prW1FQ6HA8uWLUNiYiKOHDnisRYeBhKRLF8PA5Ws9AcGBnD58mWcOXMGfX19\nSEpKwtKlSxEVFfXE8QwrIpIlF1Z2ux1fffWV7OuUdLeJiIjAzJkzMWnSJEyaNAkvvvgirly5IhtW\nPAz04NFzQk+6jcf3JnpI7hzVjBkzEB0d7bo9Skl3m9WrV+P8+fMYHBxEX18fLl68iNjYWNlauLIi\nIln+7G6zcOFCpKamIiEhARMmTEBOTo7HsHqqu9uMBn/98T0Nf3bkX6PR3eaVV15RNPbUqVPsbkNE\nY0dNV7AzrIhIFnddGEd4uEbjGVdWRCQEhhURCYFhRURCYFgRkRAYVkQkBH4aSERC4MqKiITAsCIi\nITCsiEgIagorr1vEbNy4EVqtFvHx8a7HHA4HzGYzoqOjkZKSAqfT6dciiWhsjKRhxGjzGlbZ2dko\nKytze8xqtcJsNqOlpQXJycmwWq1+K5CIxo5QYfXCCy8gODjY7bGSkhJYLBYAgMViwfHjx/1THRGN\nKV/3YAeUteICgEuXLiEgIAAff/yxx1p82inUZrNBq9UCALRaLWw2my/TEJHK+bqyetiKq6ysDNeu\nXUNRURGampqeOG779u1ITU31ukIb8bbG3GKXaPzyZysuAHj33Xexbt06zJo1y2stPoWVVqtFb28v\nAKCnpwchISG+TENEKudrWClpxdXd3Y0TJ05gy5YtALxvt+RTWKWlpaGwsBAAUFhYiPT0dF+mISKV\nkwunW7duobu723V7lJKjrdzcXFitVmg0GkUn6r1eZ5WZmYmqqirY7XZERERg165dePPNN7F+/Xoc\nOnQIOp0OxcXFXgsjIvHIBUhQUBCCgoJc93t6etyeV9KKq66uDhs2bADwoLXXqVOnEBgY+FgXnIfY\nMIJonBqNhhFGo1HR2Pr6erf3u3//PhYsWIAzZ85gzpw5WLJkCYqKihATE/PE12dnZ2PVqlVYu3at\n7HvwCnYikuXrrgtKWnENF1dWROPUaKysEhISFI1tbGxkKy4iGjtq+m4gw4qIZDGsiEgIDCsiEgLD\nioiEwLAiIiGwYQQRCYErKyISAsOKiITAsCIiITCsiEgIDCsiEgLDioiEoKZLF0a8BzsRjV8jacXl\nrbvNBx98AL1ej4SEBDz//PNobGz0WAtXVkQky9fDwIfdbcrLyxEWFoZnn30WaWlpbpvvzZs3D+fO\nncPUqVNRVlaG1157DTU1NbJzcmVFRLL82d0mKSkJU6dOBQCYTCZ0dXV5rIVhRUSy/Nnd5rsOHTqE\nlStXeqyFh4FEJEvuMPDevXu4d++e7OuGs0twRUUFCgoKcOHCBY/jGFZEJEsurCZOnIiJEye67t+5\nc8fteSXdbYAH2yHn5OSgrKwMwcHBHmvhYSARyRoaGlJ0e1RiYiJaW1vR3t6O/v5+fPTRR4+12Oro\n6MDatWtx9OhRREZGeq2FKysikuXrp4FKutvs2rULN2/edHVkDgwMRG1treyc7G5DNE6NRnebWbNm\nKRr75ZdfsrsNEY0dft2GiITAsCIiITCsiEgIavoiM8OKiGRxZUVEQmBYEZEQGFZEJASGFREJgWFF\nREJgWBGRENR06cKIdl3wtscyEYltJHuwjzafv8g8ODiIBQsWuO2xXFRU5LbHMr/ITDR2RuOLzN/d\ns8qT/v5+v4eWzysrJXssE5HY/NndBgC2bduGqKgo6PV61NfXe6zF57Aa7h7LRCQeX8PqYXebsrIy\nXLt2DUVFRWhqanIbU1paiuvXr6O1tRUHDx507Wslx+ew4iEe0fjnz+42JSUlsFgsAB50t3E6nbDZ\nbLK1+BxWSvdYJiJx+bO7zZPGeGrH5XNYKdljmYjE5use7EqPvB4NOk+v8zmsvrvHcmxsLH7605+6\nfRIIAD/+8Y99nZ6IRuB//bsXFBTkdl/JkdejY7q6uhAWFib7Hn7dg52Ink7379/HggULcObMGcyZ\nMwdLlix57NKm0tJS5Ofno7S0FDU1NcjNzfXYPp5XsBPRqFPS3WblypUoLS1FZGQkJk+ejMOHD3uc\nkysrIhKCX5ucqv3rOBs3boRWq0V8fLzrMYfDAbPZjOjoaKSkpMDpdI5hhe46OzuxbNkyLFq0CHFx\ncdi/fz8A9db87bffwmQywWAwIDY2Fjt27ACg3nqBB9cHGY1GrFq1CoC6a9XpdEhISIDRaMSSJUsA\nqLvekfJbWCm5KGysZWdno6yszO0xq9UKs9mMlpYWJCcnw2q1jlF1jwsMDMTevXvx2WefoaamBu+9\n9x6amppUW/MzzzyDiooKNDQ0oLGxERUVFTh//rxq6wWAffv2ITY21vWplJpr1Wg0qKysRH19vas5\nqJrrHTHJT6qrq6UVK1a47u/evVvavXu3v97OZ21tbVJcXJzr/oIFC6Te3l5JkiSpp6dHWrBgwViV\n5tXq1aul06dPC1Hz3bt3pcTEROnTTz9Vbb2dnZ1ScnKydPbsWeknP/mJJEnq/v9Bp9NJdrvd7TE1\n1ztSfltZifp1HJvNBq1WCwDQarUer6gdS+3t7aivr4fJZFJ1zUNDQzAYDNBqta5DWLXW+/rrr+Od\nd97BhAn//2uh1lqBByurl19+GYmJifjTn/4EQN31jpTfPg0cD1/H0Wg0qvw57ty5g4yMDOzbtw9T\npkxxe05tNU+YMAENDQ34+uuvsWLFClRUVLg9r5Z6T548iZCQEBiNRlRWVj5xjFpqfejChQsIDQ3F\nl19+CbPZjIULF7o9r7Z6R8pvKytRv46j1WrR29sLAOjp6UFISMgYV+RuYGAAGRkZyMrKQnp6OgD1\n1wwAU6dOxauvvoq6ujpV1ltdXY2SkhLMnTsXmZmZOHv2LLKyslRZ60OhoaEAgFmzZmHNmjWora1V\ndb0j5bewEvXrOGlpaSgsLAQAFBYWugJBDSRJwqZNmxAbG4vc3FzX42qt2W63uz6N+uabb3D69GkY\njUZV1puXl4fOzk60tbXh2LFjWL58OY4cOaLKWgGgr68Pt2/fBgDcvXsX//znPxEfH6/aekeFP0+I\nlZaWStHR0dL8+fOlvLw8f76VTzZs2CCFhoZKgYGBUnh4uFRQUCB99dVXUnJyshQVFSWZzWbp5s2b\nY12my7/+9S9Jo9FIer1eMhgMksFgkE6dOqXamhsbGyWj0Sjp9XopPj5e2rNnjyRJkmrrfaiyslJa\ntWqVJEnqrfWLL76Q9Hq9pNfrpUWLFrl+v9Ra72jgRaFEJAS/XhRKRDRaGFZEJASGFREJgWFFREJg\nWBGREBhWRCQEhhURCYFhRURC+D9BMQC0CcH9CAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x41c0350>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# you can also select a single line from the 2D array and plot the values of a single line or column:\n",
      "lin = r[:,12]\n",
      "xlin = x1[:,12]\n",
      "print lin.shape\n",
      "plot(xlin, lin)\n",
      "minlin = np.min(lin)\n",
      "print 'minimum value of lin is ', minlin"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(51,)\n",
        "minimum value of lin is "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 3.0\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEACAYAAACuzv3DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVWW+x/HP9lY5mtptY1HpoIQoAV20i9omA8oL3srR\nypik5mRjk6NjZlMj1WRYdlGrk6es6GZaKVoZgeb2Vq884yU9leVY5A0wUhzxhsA6fzyJkQqy2bD2\n2vv7fr14hXvD2j9W+mXtZz3P73FZlmUhIiKO1cjuAkREpG4U5CIiDqcgFxFxOAW5iIjDKchFRBxO\nQS4i4nDVBvnBgwfp1q0bcXFxREdHM2HCBADS09MJDw8nPj6e+Ph4srOzG6RYERE5lqumeeT79++n\nefPmlJWV0b17d6ZMmcLixYtp2bIlY8aMaag6RUTkBGocWmnevDkApaWllJeX06ZNGwC0jkhEJDDU\nGOQVFRXExcXhdrtJSEigc+fOAEyfPp3Y2FjS0tIoLi6u90JFROT4ahxaOWLPnj0kJyeTkZFBdHQ0\nZ599NgAPPfQQ+fn5zJw5s14LFRGR42tysl/YqlUr+vTpw7/+9S88Hk/l43fccQf9+vU75us7dOjA\n5s2b/VKkiEioiIiI4N///netvqfaoZWioqLKYZMDBw6Qm5tLfHw8BQUFlV8zb948YmJijvnezZs3\nY1mWPiyLiRMn2l5DoHzoXOhc6FxU/+HLBXC1V+T5+fmkpqZSUVFBRUUFw4cPp1evXtx2222sW7cO\nl8tF+/btmTFjRq1fWERE/KPaII+JiWHNmjXHPP7666/XW0EiIlI7WtnZAH59TyHU6VwcpXNxlM5F\n3Zz0rJVaH9jlop4OLSIStHzJTl2Ri4g4nIJcRMThFOQiIg6nIBcRcTgFuYiIwynIRUQcTkEuIuJw\nCnIREYdTkIuIOJyCXETE4RTkIiIOpyAXEXE4BbmIiMMpyEVEHE5BLiLicApyERGHU5CLiDicglxE\nxOEU5CIiDqcgFxFxOAW5iIjDVRvkBw8epFu3bsTFxREdHc2ECRMA2LVrF4mJiURGRpKUlERxcXGD\nFCsiIsdyWZZlVfcF+/fvp3nz5pSVldG9e3emTJnCggULOOuss7jvvvuYPHkyu3fvJiMjo+qBXS5q\nOLSIiPyGL9lZ49BK8+bNASgtLaW8vJw2bdqwYMECUlNTAUhNTSUrK8uHckVExB9qDPKKigri4uJw\nu90kJCTQuXNnCgsLcbvdALjdbgoLC+u9UBGRYPfss759X5OavqBRo0asW7eOPXv2kJyczJIlS6o8\n73K5cLlcx/3e9PT0ys89Hg8ej8e3KkVEgpTX62XJEi/LlsGXX/p2jBqD/IhWrVrRp08fVq9ejdvt\npqCggLCwMPLz8znnnHOO+z2/DnIRETnWNdd4+OgjD0VF8PXX0Lbtw7U+RrVDK0VFRZUzUg4cOEBu\nbi7x8fGkpKSQmZkJQGZmJgMGDPChfBGR0FZRASNHwtKl5iMszLfjVDtrZcOGDaSmplJRUUFFRQXD\nhw9n3Lhx7Nq1iyFDhrBlyxbatWvHnDlzaN26ddUDa9aKiMgJlZXBH/8IW7fCBx/A6aebx33Jzhqn\nH/pKQS4icnyHDsHQoea/770Hv0wOBOpp+qGIiPjPvn2QkgKNG0NWVtUQ95WCXESkgezZA8nJ0LYt\nvPMONGvmn+MqyEVEGkBREVx7LcTFwSuvQJOTnjNYMwW5iEg927EDrrnGXI1Pnw6N/Jy8CnIRkXqU\nlwc9e8Ktt8KkSXCC9ZN1oiAXEaknGzeaEB89Gn5pHlsv/DhKIyIiR3z5Jdxwg7kK/+Mf6/e1FOQi\nIn72+ecwYAA8/zzceGP9v56CXETEjxYvhmHDIDPTXJE3BI2Ri4j4yQcfmBB/992GC3FQkIuI+MU7\n78Cdd8JHH5mphg1JQS4iUkcvvwxjx0JuLlx+ecO/vsbIRUTq4JlnYOpU8HqhY0d7alCQi4j4wLLg\n0UfhzTdh2TK44AL7alGQi4jUkmXBuHGQkwPLl8MvWxjbRkEuIlIL5eVmV5/1681wyhln2F2RglxE\n5KQdPgypqZCfb25stmxpd0WGglxE5CQcPAh/+IO5Il+4EE47ze6KjtL0QxGRGpSUQN++cOqpMHdu\nYIU4KMhFRKpVXAxJSXDhhfD22/7b1cefFOQiIiewcyckJJhFPi+9ZPbZDEQKchGR49i2zfQS79cP\nnn3W/7v6+FMAlyYiYo/Nm6FHD0hLg0ceqZ9dffyp2iDfunUrCQkJdO7cmS5dujBt2jQA0tPTCQ8P\nJz4+nvj4eLKzsxukWBGR+vbVV6bp1f33m0U/TuCyLMs60ZMFBQUUFBQQFxdHSUkJl156KVlZWcyZ\nM4eWLVsyZsyYEx/Y5aKaQ4uIBJzVq6FPH3jqKbjlFntq8CU7q51HHhYWRlhYGAAtWrSgU6dObN++\nHUAhLSJBZflyGDzY3NTs39/uamrnpMfI8/LyWLt2LVdccQUA06dPJzY2lrS0NIqLi+utQBGR+vbJ\nJzBoELz1lvNCHGoYWjmipKQEj8fDgw8+yIABA9i5cydnn302AA899BD5+fnMnDmz6oFdLiZOnFj5\nZ4/Hg8fj8W/1IiJ1NHcu3HUXZGXBVVc1/Ot7vV68Xm/lnx9++OFaj3jUGOSHDx+mb9++3HDDDYwe\nPfqY5/Py8ujXrx8bNmyoemCNkYtIgHv9dRg/3iy5j4+3uxrDl+ysdmjFsizS0tKIjo6uEuL5+fmV\nn8+bN4+YmJhalioiYq8XXoC//x0+/TRwQtxX1V6Rr1ixgp49e3LxxRfj+mUi5aRJk5g1axbr1q3D\n5XLRvn17ZsyYgfs3DXl1RS4igSojw9zUXLQI2re3u5qqfMnOkxojb6hiRETqk2WZq/CsLNOG9rzz\n7K7oWH6ffigiEiwqKuDee+Gzz8zWbGedZXdF/qMgF5GgV1Zmltt//70ZE2/Vyu6K/EtBLiJB7dAh\nuPlm01P8k0+geXO7K/I/Nc0SkaC1f//RBT4LFgRniIOCXESC1J49kJxsdrifPRtOOcXuiuqPglxE\ngk5REVx7LcTGwquvQpMgH0RWkItIUNmxw7ShTU6G6dMDe0MIfwmBH1FEQsUPP5gNIYYPh0mTAn9D\nCH9RkItIUPjmG7M129ixZlOIUBLkI0ciEgrWrDEbQjzxhLkaDzUKchFxtBUrTC/xGTNg4EC7q7GH\nglxEHCsnx2zJ9tZbkJRkdzX20Ri5iDjSvHlmGCUrK7RDHBTkIuJAb7wBd98NH38MV19tdzX209CK\niDjKCy+YfuKffgqdOtldTWBQkIuIYzz+OMycadrQtmtndzWBQ0EuIgHPsuCBB0zjq2XL4Nxz7a4o\nsCjIRSSgVVTAPffAqlWwdGlwbQjhLwpyEQlYZWUwYgT8+CMsXgynn253RYFJQS4iAenQIRg6FA4e\nNLNTgrWXuD9o+qGIBJx9+6BvX2jaFObPV4jXREEuIgGluBgSE+H882HWLGjWzO6KAp+CXEQCxs6d\n4PFA167w8svQuLHdFTlDtUG+detWEhIS6Ny5M126dGHatGkA7Nq1i8TERCIjI0lKSqK4uLhBihWR\n4LV1q2lD278/PPNMaGwI4S8uy7KsEz1ZUFBAQUEBcXFxlJSUcOmll5KVlcWrr77KWWedxX333cfk\nyZPZvXs3GRkZVQ/sclHNoUVEKm3aZIZT7rnH9BMPZb5kZ7W/88LCwoiLiwOgRYsWdOrUie3bt7Ng\nwQJSU1MBSE1NJSsr67jfv3t3rWoRkRC0YYMZTnngAYW4r076zUteXh5r166lW7duFBYW4na7AXC7\n3RQWFh73ezweOMFTIiJ88QVcdx1MmQJ/+pPd1TjXSc0jLykpYfDgwUydOpWWLVtWec7lcuE6wcZ4\nZ5yRTnQ03HYb9O/vwePx1LlgEQkOS5bAkCFml/u+fe2uxj5erxev11unY1Q7Rg5w+PBh+vbtyw03\n3MDo0aMBiIqKwuv1EhYWRn5+PgkJCWzcuLHqgX8Z55k6FZ5+GnJzITKyTrWKSJD44AOzYnPOHEhI\nsLuawOL3MXLLskhLSyM6OroyxAFSUlLIzMwEIDMzkwEDBpzwGPfeCxMnmmGWL7+sVW0iEoRmzYI7\n7oCPPlKI+0u1V+QrVqygZ8+eXHzxxZXDJ48//jhdu3ZlyJAhbNmyhXbt2jFnzhxat25d9cC/+a3y\n7rswapTZzePKK+vppxGRgPY//wOPPALZ2dCli93VBCZfrshrHFrxZzHZ2Wa8/O23zQ0OEQkdU6aY\nTSFycyEiwu5qAlfABznA8uUweLD5zVzNiIyIBAnLgocegvffNyEeHm53RYHNlyBv8O6HPXqYK/M+\nfWDvXrN5qogEp4oKc59s5UqzIcTZZ9tdUXCypY3tJZeY/faSkkyY3323HVWISH0qK4O0NPj+ezPV\nsFUruysKXrb1I+/UyfyGTkyEPXvg/vvhBNPRRcRhjvQSP3AAPvlEbWjrm61tadq3N2H+1lsmyNWa\nRcT5jvQSb9xYvcQbiu39xc491+zDt2QJjBwJ5eV2VyQivjrSSzw8HN55B045xe6KQoPtQQ5w5plm\nP75vvzU3Pw8ftrsiEamtwsKjvcRnzoQm2kiywQREkAO0bAkLF5qbn4MGmbE1EXGGLVvMjLQBA9RL\n3A4BdbpPOw3mzoUWLaB3bxPqIhLYvv3WhPjIkZCerkkLdgioIAez2eqbb8JFF0GvXvDzz3ZXJCIn\nsm6d6ZcycSL89a92VxO6Ai7Iwdzt/u//Nn9BrrkGduywuyIR+a3PPoPkZJg61XQyFPsE7O0Ilwsm\nT4bWrc3btkWLzHRFEbFfbi7cfDO88QZcf73d1UjABvkREyaYFWE9e5qFBdHRdlckEtrmzoW77oJ5\n86B7d7urEXBAkINZwn/66XDttfDhh3DZZXZXJBKaXnvNXFxlZ5tWGxIYHBHkALfeaqYo9u5teptf\nc43dFYmElqlT4amnzOK9qCi7q5Ffa/A2tnX16aemh8Orr5oOiiJSvywLHn3UzCbLzYULL7S7ouDm\niH7k/vDFF9C/v1l4MGxYvbyEiGDa0I4day6gcnLA7ba7ouDniH7k/tCtm5nFkpxsOifedZfdFYkE\nn7Iy+NOfYONG8HqhTRu7K5ITcWSQg9nv79dtcMePt7sikeBx6BDccov5t5WTY1ZbS+AKyAVBJysi\nwmwd9/rr5k662uCK1N2+fZCSYoZVPvxQIe4Ejg5ygPPOM21wFy0y0xTVBlfEd7t3m527zj0X5sxR\nG1qncHyQA5x1lmmDu3Gj2uCK+OpIG9rLL1cbWqcJiiAHs2Bo4UIoKYGBA9UGV6Q2fvzRtMIYNEht\naJ2oxv9dI0aMwO12ExMTU/lYeno64eHhxMfHEx8fT3Z2dr0WebJOOw3ef9/0Z7n+evjPf+yuSCTw\nbdxoQvzPfzZdDNWG1nlqDPLbb7/9mKB2uVyMGTOGtWvXsnbtWq4PoK45TZuam58xMaZ74k8/2V2R\nSOBas8b8O3n0Ubj3XrurEV/VGOQ9evSgzXEmkNbXYh9/aNQIpk+HG24wzba2bbO7IpHAs2yZeef6\nwguQmmp3NVIXPo+ETZ8+ndjYWNLS0iguLvZnTX7hcsE//wl33GHeNm7aZHdFIoFj4UK48UZ4+21z\nT0mczaf70iNHjuQf//gHAA899BBjx45l5syZx3xdenp65ecejwePx+NTkXUxdqwZM/d4zF/e2NgG\nL0EkoMyaZXbz+eADs0pa7OX1evF6vXU6xkn1WsnLy6Nfv35s2LDhpJ+rz14rvnjvPXMzZ+5cuPpq\nu6sRsceLL5p3qtnZZnW0BB5fstOnoZX8/PzKz+fNm1dlRkuguvFGcxN04ECzQYVIKLEsePxxePJJ\nMzauEA8uNV6RDxs2jKVLl1JUVITb7ebhhx/G6/Wybt06XC4X7du3Z8aMGbh/0xYt0K7Ij/jsMxPm\nzz0HN91kdzUi9c+yTC+ijz82FzHnnmt3RVKdkGljW1fr15sZLenpcOeddlcjUn/Ky2HkSPjySxPk\nZ5xhd0VSk5BpY1tXF19s2nImJZneEvfdZ3dFIv536JBpWbFrl2lhoeZXwStkF+J27AgrVkBmJtx/\nvzonSnA50sGwvBw++kghHuxCNsjBdE5ctszsfnLXXeqcKMFh1y7Tp/+882D2bHUwDAUhHeQAZ55p\n3nb++99m27jSUrsrEvFdfr7ZmPzKK9XBMJSEfJADtGxp3n4ePmzeju7bZ3dFIrX3/ffQvbu5IJky\nRc2vQomC/Benngrvvgtt25q3pbt3212RyMn7v/8zfYX+9jd44AGFeKhRkP9Kkybm7egVV5i3p79a\n9yQSsD7/HHr1Mot9Ro60uxqxg4L8Nxo1gqeegj/8wbxN/f57uysSObFPPjHDga+9ZoZUJDTpVshx\nuFzw97+bG6E9e5qFFA7oQiAhZs4cuOceyMpS/6BQpyCvxl13QZs2cN11MG8eXHWV3RWJGDNmwCOP\nQG6uWeAmoS0kl+jXVnY23HababoVQJshSQiyLMjIgJdeMiEeEWF3ReJv6rVSj44025o2zYyfizQ0\ny4Jx48y4uJpfBS/1WqlHV10FixaZZlu7dml2gDSssjLT4G3jRli6VM2vpCoFeS3ExJgl/UlJUFQE\nDz6o+bpS/w4cgKFDzarjRYvgd7+zuyIJNJp+WEu//71ptvX++zB6NFRU2F2RBLM9e8y7wObNYf58\nhbgcn4LcB2Fhpg3umjVm9/HDh+2uSILRzp2QkGB283nrLWjWzO6KJFApyH3UurW54VRcbG6C7t9v\nd0USTPLyzIK0lBSYPt0sVBM5Ef31qIPmzc1mzmecYcbNi4vtrkiCwVdfQY8e8Je/mF2sdB9GaqIg\nr6OmTc3y6MsvN6tA1Z9F6uJI35TJk2HUKLurEadQkPtBo0bw9NOm18XVV5ve5iK19fHH0L+/uTC4\n+Wa7qxEn0fRDP3G5YMIEOPtsc2X+0UcQH293VeIUb70FY8eamSlXXml3NeI0CnI/u+MOM2aenGya\nGnk8dlckgW7qVNNxc/Fi6NzZ7mrEibREv54sWWKW8s+YYWa1iPyWZZlFZe+9Bzk5cOGFdlckgcCX\n7KxxjHzEiBG43W5iftXHddeuXSQmJhIZGUlSUhLFmq5xjIQEMz3xz3+Gl1+2uxoJNOXlprtmTo5Z\nYKYQl7qoMchvv/12srOzqzyWkZFBYmIi3333Hb169SIjI6PeCnSy+HizpH/SJPMRwm9Q5FcOHoSb\nboLNm+HTT819FZG6qDHIe/ToQZs2bao8tmDBAlJTUwFITU0lKyurfqoLAh06wMqVMHs23HuvlvSH\nuuJic/+kWTNzQ7xlS7srkmDg0/TDwsJC3G43AG63m8LCQr8WFWzatjVX5uvXm2llhw7ZXZHYIT/f\n7AUbGwtvvw2nnGJ3RRIs6jxrxeVy4TrB0rP09PTKzz0eD54QnsLRqpXZoOKWW6BvX7MiVFdjoWPT\nJnMlfscdZpqqVmvKEV6vF6/XW6djnNSslby8PPr168eGDRsAiIqKwuv1EhYWRn5+PgkJCWzcuLHq\ngUN81sqJlJfD3XfD6tWwcCGcc47dFUl9W73a/PJ+9FET5CLVqZdZK8eTkpJCZmYmAJmZmQwYMMCX\nw4Skxo3hxRehd2+zCvT77+2uSOrTkc1IXnxRIS71p8Yr8mHDhrF06VKKiopwu9088sgj9O/fnyFD\nhrBlyxbatWvHnDlzaN26ddUD64q8Rs8/b2azfPihVoEGo7ffNj3r33/fNMESORnas9OB3nvPDLXM\nmmWaJUlwePppePZZM3zWpYvd1YiTKMgdaulSM6942jSzpZc4V0UF3HefCfBPPoHzz7e7InEabb7s\nUNdcY/ps9O4NBQXm7bg4T2kp3H47/PijWa2pDZKloSjIA0RMjFk4lJwMO3ZARoZ2hXGSvXth0CCz\np2ZuLpx2mt0VSShRVASQCy4wV3IrV8Jtt5krPAl8hYWmy+Xvf2/ueSjEpaEpyAPMmWeaKWt795q5\nx//5j90VSXW++w6uusrsrfnii9BE73HFBgryAHTaaWbKWocOZpOKHTvsrkiO57PPzP+fv/8dJk7U\nak2xj4I8QDVpYuaZDx1qrvi++sruiuTX5s2DAQPMtmwjRthdjYQ6vREMYC4X3H8/hIfDtdeaHYeu\nucbuquTIQq6PP4ZLL7W7GhHNI3eMxYvN5s7Tp5udh6ThVVSYhlfz55sQb9/e7ookGGkeeRDr1ctM\na+vbF7ZtgzFjNCbbkA4dOjpHfOVKc1NaJFBojNxBYmPNDbZXX4V77oGyMrsrCg27dkFSkgnzRYsU\n4hJ4FOQOc/755orw22/NzbaSErsrCm6bN5ubzZdfDu++qzniEpgU5A7UqpXp5eF2a3piffr8c+je\nHf7yF5gyRSttJXDpr6ZDNW0KL78MN94IV14Jv+z5IX7y3ntmkc/LL5vulCKBTLNWgsCsWWZj5zff\nNGO54jvLMlffU6fCBx+oT7w0PLWxDWHLl5tWuI8+CnfeaXc1znT4sBlGWbnS7HCvFrRiBwV5iPvu\nOzM9sW9fePJJs62cnJzdu2HIELOidvZsOP10uyuSUNVge3ZKYIqMhC++gPXrzfiuGm6dnO++gyuu\nMDv5fPihQlycR0EeZNq0MasOL7jA3ATV5s7VW7zY7Kf5t7/BM8/oXYw4k4I8CDVtCi+8ACNHwtVX\nm/FzOdaLL8Itt5ihFN1XECfTGHmQy8mBW2+FyZPNEnMxK2LHjDEtDz74wLQLFgkUutkpx/XNN9Cv\nn/l44glzxR6qfv7ZNB9zucyVeOvWdlckUpVudspxdeoEq1aZQE9MhJ077a7IHmvXwmWXQVycmV6o\nEJdgUacr8nbt2nH66afTuHFjmjZtyqpVq44eWFfkAae83Oxk8/rrZuVi1652V9Rw3njDDKc8/7yZ\nZigSqBp8aKV9+/asXr2aM844wy/FSMOYNw/+67/g8cchLc3uaurX4cMmwLOzzc/dpYvdFYlUz5Z+\n5Apr5xk40Ay3DBwI//u/Zjn6KafYXZX/FRSY1a6tW5ufU0MpEqzqNEbucrm47rrruOyyy3jppZf8\nVZM0gKgoM27+009m+7gff7S7Iv9ascKMh/fqZXb0UYhLMKvTFfnKlStp27YtP/30E4mJiURFRdGj\nR4/K59PT0ys/93g8eDyeuryc+FnLlmas/KmnTL/t5583V7BOVlYG//wnzJgBM2dC7952VyRSPa/X\ni9frrdMx/Db98OGHH6ZFixaMHTvWHFhj5I7yr3+ZaXk9e8K0afC739ldUe39+KNZ4NO8OWRmQtu2\ndlckUnsNOv1w//797N27F4B9+/aRk5NDTEyMr4cTm112GaxZYzYYvuQS87mTzJ5t3lUMGGBubCrE\nJZT4PLRSWFjIwIEDASgrK+OWW24hSc2wHa1lS7Mf6DvvwPXXw/jx8Ne/BvbOOCUlpvXsihWmx8yl\nl9pdkUjD08pOOa68PLj5ZjPE8sIL0LGj3RUda8kSM42ye3czHNSihd0VidSdVnaK37RrB8uWQXKy\n6aI4cSIcOGB3VUZ+vvklc/vtpu/6K68oxCW0KcjlhJo0Me1d160zy/u7dDFL2+1SVmbmvF98sflF\n89VX0L+/ffWIBAoNrchJy8mBUaMgOtoE6oUXNtxrf/aZ2QT5zDPNNMmoqIZ7bZGGpKEVqVdJSbBh\ng5nhcsklcP/9sGVL/b7m6tWmDe9NN5nXW7RIIS7yWwpyqZVTToEHHzQBe+iQ2WV+0CD49FOzA70/\nHDoEb71lxuYHD4aYGDO0M3SoaT8rIlVpaEXqpKTEhO706WYO+qhRMHy4mcpYW9u2mRWZL71kwnvU\nKLORtLZfk1CijSXENpYFS5fCc8+ZfTBjY80QyEUXmf9GRZl9RBs3NlfcmzfDxo3m49tvzX83bTIr\nM+++2zT1EglFCnIJCDt3wvr1RwP6SFgXFcHZZ0NhoblR+uugv+giE/6aRiihTkEuAW3fPtNa9vzz\noVkzu6sRCUwKchERh9P0QxGREKQgFxFxOAW5iIjDKchFRBxOQS4i4nAKchERh1OQi4g4nIJcRMTh\nFOQiIg6nIBcRcTgFuYiIwynIRUQczucgz87OJioqio4dOzJ58mR/1iQiIrXgU5CXl5czatQosrOz\n+frrr5k1axbffPONv2sLGl6v1+4SAobOxVE6F0fpXNSNT0G+atUqOnToQLt27WjatClDhw5l/vz5\n/q4taOgv6VE6F0fpXBylc1E3PgX59u3bOf/88yv/HB4ezvbt2/1WlIiInDyfgtylrcxFRAKH5YPP\nP//cSk5OrvzzpEmTrIyMjCpfExERYQH60Ic+9KGPWnxERETUOpN92uqtrKyMiy66iMWLF3PuuefS\ntWtXZs2aRSdtfS4i0uCa+PRNTZrw3HPPkZycTHl5OWlpaQpxERGb1NvmyyIi0jD8vrJz3LhxdOrU\nidjYWAYNGsSePXsqn3v88cfp2LEjUVFR5OTk+PulA867775L586dady4MWvWrKnyXKidCwjtRWQj\nRozA7XYTExNT+diuXbtITEwkMjKSpKQkiouLbayw4WzdupWEhAQ6d+5Mly5dmDZtGhCa5+PgwYN0\n69aNuLg4oqOjmTBhAuDDufDlZmd1cnJyrPLycsuyLGv8+PHW+PHjLcuyrK+++sqKjY21SktLrR9+\n+MGKiIio/Lpg9c0331jffvut5fF4rNWrV1c+HornoqyszIqIiLB++OEHq7S01IqNjbW+/vpru8tq\nMMuWLbPWrFljdenSpfKxcePGWZMnT7Ysy7IyMjIq/60Eu/z8fGvt2rWWZVnW3r17rcjISOvrr78O\n2fOxb98+y7Is6/Dhw1a3bt2s5cuX1/pc+P2KPDExkUaNzGG7devGtm3bAJg/fz7Dhg2jadOmtGvX\njg4dOrBq1Sp/v3xAiYqKIjIy8pjHQ/FchPoish49etCmTZsqjy1YsIDU1FQAUlNTycrKsqO0BhcW\nFkZcXBwALVq0oFOnTmzfvj1kz0fz5s0BKC0tpby8nDZt2tT6XNRr06xXXnmF3r17A7Bjxw7Cw8Mr\nnwvlRUQgXjXNAAACfUlEQVSheC60iOxYhYWFuN1uANxuN4WFhTZX1PDy8vJYu3Yt3bp1C9nzUVFR\nQVxcHG63u3LIqbbnwqdZK4mJiRQUFBzz+KRJk+jXrx8Ajz32GM2aNePmm28+4XGCYWHRyZyLkxEM\n56I6wf7z1ZXL5Qq5c1RSUsLgwYOZOnUqLVu2rPJcKJ2PRo0asW7dOvbs2UNycjJLliyp8vzJnAuf\ngjw3N7fa51977TUWLlzI4sWLKx8777zz2Lp1a+Wft23bxnnnnefLyweUms7F8QTruajOb3/mrVu3\nVnlXEorcbjcFBQWEhYWRn5/POeecY3dJDebw4cMMHjyY4cOHM2DAACC0zwdAq1at6NOnD6tXr671\nufD70Ep2djZPPvkk8+fP59RTT618PCUlhXfeeYfS0lJ++OEHNm3aRNeuXf398gHL+tUsz1A8F5dd\ndhmbNm0iLy+P0tJSZs+eTUpKit1l2SolJYXMzEwAMjMzKwMt2FmWRVpaGtHR0YwePbry8VA8H0VF\nRZUzUg4cOEBubi7x8fG1Pxf+vgPboUMH64ILLrDi4uKsuLg4a+TIkZXPPfbYY1ZERIR10UUXWdnZ\n2f5+6YAzd+5cKzw83Dr11FMtt9ttXX/99ZXPhdq5sCzLWrhwoRUZGWlFRERYkyZNsrucBjV06FCr\nbdu2VtOmTa3w8HDrlVdesX7++WerV69eVseOHa3ExERr9+7ddpfZIJYvX265XC4rNja2Mic+/vjj\nkDwf69evt+Lj463Y2FgrJibGeuKJJyzLsmp9LrQgSETE4bTVm4iIwynIRUQcTkEuIuJwCnIREYdT\nkIuIOJyCXETE4RTkIiIOpyAXEXG4/wcZhK2HbFrjewAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x42ebbd0>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# you can make other rectangular shapes and combine them together\n",
      "mask2 = (x>10) * (x<12)\n",
      "mask3 = (y>20) * (y < 41)\n",
      "mask_combined = np.logical_or(mask2, mask3) * 3.2\n",
      "imshow(mask_combined)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 58,
       "text": [
        "<matplotlib.image.AxesImage at 0x53e4a10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAD+CAYAAADVndu7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADlRJREFUeJzt2l9sU3Ufx/HPGS55QoIMhB4Xppl/2HBjrI0TboyAs+MC\nN0aWIMYsC6Dxxgu4EOHKOyhCYkbilRHTeIFyNSfBhRk2FIlZQiCAYmYik4ltFbYhfwNsv+fCUJ+J\nD+3Gztbyfb+SXvTPTr+0v745PT2ec84JgEkFUz0AgKlDAADDCABgGAEADCMAgGEEADDsoWweVFpa\nqocffljTpk1TYWGhenp6NDAwoFdeeUW//PKLSktLtW/fPhUVFQU9L4AJlNUegOd56u7u1vHjx9XT\n0yNJisViikaj6u3tVW1trWKxWKCDAph4WX8F+Of5Qu3t7WppaZEktbS0qK2tbWInAxA4L5szAZ98\n8knNnDlT06ZN05tvvqk33nhDs2bN0uDgoKS/4jB79uz09fTGPS+YqQFklNVJvi4Lv/32m3POud9/\n/91VV1e7r7/+2hUVFY16zKxZs+76O0lcJvGSD959992pHmFM/nfeqX5/g1gPWX0FKC4uliTNnTtX\nq1evVk9Pj3zfVzKZlCQlEgmFQqFsNgUgh2QMwLVr13T58mVJ0tWrV3Xw4EFVVVWpoaFB8XhckhSP\nx9XY2BjspAAmXMafAVOplFavXi1Jun37tl577TXV1dWppqZGa9as0UcffZT+GRDIZNmyZVM9wpjk\n27xjldVBwHFvnIOAkyrAtxLKv/WczXrgTEDAMAIAGEYAAMMIAGAYAQAMIwCAYQQAMIwAAIYRAMAw\nAgAYRgAAwwgAYBgBAAwjAIBhBAAwjAAAhhEAwDACABhGAADDCABgGAEADCMAgGEEADCMAACGEQDA\nMAIAGEYAAMMIAGAYAQAMIwCAYQQAMIwAAIYRAMAwAgAYRgAAw7IKwPDwsCKRiOrr6yVJAwMDikaj\nKisrU11dnYaGhgIdEkAwsgpAa2urKioq5HmeJCkWiykajaq3t1e1tbWKxWKBDgkgGBkD8Ouvv+rA\ngQN6/fXX5ZyTJLW3t6ulpUWS1NLSora2tmCnBBCIjAHYtGmTdu7cqYKCvx+aSqXk+74kyfd9pVKp\n4CYEEJh7BmD//v0KhUKKRCLp//3/yfO89FcDAPnloXvdefToUbW3t+vAgQO6ceOG/vzzTzU3N8v3\nfSWTST366KNKJBIKhUKTNS+ACeS5//df+z8cPnxYu3bt0hdffKHNmzfrkUce0TvvvKNYLKahoaF/\nPRDInsHkyvKtxDjl23rOZj2M6TyAOy/Ali1b1NnZqbKyMh06dEhbtmwZ34QAplTWewDj2nieFTPf\nsQcQrHxbzxO+BwDgwUIAAMMIAGAYAQAMIwCAYQQAMIwAAIYRAMAwAgAYRgAAwwgAYBgBAAwjAIBh\nBAAwjAAAhhEAwDACABhGAADDCABgGAEADCMAgGEEADCMAACGEQDAMAIAGEYAAMMIAGAYAQAMIwCA\nYQQAMIwAAIYRAMAwAgAYRgAAwwgAYNg9A3Djxg0tWbJE4XBYFRUV2rp1qyRpYGBA0WhUZWVlqqur\n09DQ0KQMC2Biec45d68HXLt2TdOnT9ft27f1/PPPa9euXWpvb9ecOXO0efNm7dixQ4ODg4rFYndv\n3PMCGxx3y/BW4j7l23rOZj1k/Aowffp0SdLNmzc1PDysWbNmqb29XS0tLZKklpYWtbW13eeoAKZC\nxgCMjIwoHA7L930tX75clZWVSqVS8n1fkuT7vlKpVOCDAph4D2V6QEFBgU6cOKFLly5pxYoV6urq\nGnW/53l5t2sE4C9Z/wowc+ZMrVy5UseOHZPv+0omk5KkRCKhUCgU2IAAgnPPAFy4cCF9hP/69evq\n7OxUJBJRQ0OD4vG4JCkej6uxsTH4SQFMuHv+CnDq1Cm1tLRoZGREIyMjam5u1ttvv62BgQGtWbNG\n586dU2lpqfbt26eioqK7N85Xg0nFrwDByrf1nM16yPgz4P3Itxcs3xGAYOXbep6QnwEBPLgIAGAY\nAQAMIwCAYQQAMIwAAIYRAMAwAgAYRgAAwwgAYBgBAAwjAIBhBAAwjAAAhhEAwDACABhGAADDCABg\nGAEADCMAgGEEADCMAACGEQDAMAIAGEYAAMMIAGAYAQAMIwCAYQQAMIwAAIYRAMAwAgAYRgAAwwgA\nYFjGAPT392v58uWqrKzUwoULtXv3bknSwMCAotGoysrKVFdXp6GhocCHBTCxPOecu9cDksmkksmk\nwuGwrly5omeffVZtbW36+OOPNWfOHG3evFk7duzQ4OCgYrHY6I17XqDDY7QMbyXuU76t56zWgxuj\nVatWuc7OTldeXu6SyaRzzrlEIuHKy8vveqwkLpN4QbCm+v0NYj1k3AP4X319fVq6dKlOnz6txx9/\nXIODg9Jfz6TZs2enr9+Rb8XMd2N4KzEO+baes1kPWR8EvHLlipqamtTa2qoZM2aMus/zvLx7cQBk\nGYBbt26pqalJzc3NamxslCT5vq9kMilJSiQSCoVCwU0JIBAZA+Cc04YNG1RRUaGNGzemb29oaFA8\nHpckxePxdBgA5I+MxwCOHDmiF154QYsWLUrv5m/fvl2LFy/WmjVrdO7cOZWWlmrfvn0qKioavXG+\nFkwqjgEEK9/WczbrYUwHAccq316wfEcAgpVv63lCDwICePAQAMAwAgAYRgAAwwgAYBgBAAwjAIBh\nBAAwjAAAhhEAwDACABhGAADDCABgGAEADCMAgGEEADCMAACGEQDAMAIAGEYAAMMIAGAYAQAMIwCA\nYQQAMIwAAIYRAMAwAgAYRgAAwwgAYBgBAAwjAIBhBAAwjAAAhj0U9BM454J+CmBS5NNa9jwvq8ex\nBwAYljEA69evl+/7qqqqSt82MDCgaDSqsrIy1dXVaWhoKNAhAQQjYwDWrVunjo6OUbfFYjFFo1H1\n9vaqtrZWsVgssAEBBMdzWXyx6evrU319vU6dOiVJWrBggQ4fPizf95VMJrVs2TL9+OOPd2/c8/Lq\nexPwoMj2szeug4CpVEq+70uSfN9XKpW65yCYHMQ2WA/iWr7vg4Ce5z2QLwxgwbgCcGfXX5ISiYRC\nodCEDgVgcowrAA0NDYrH45KkeDyuxsbGCR0KwCRxGaxdu9YVFxe7wsJCV1JS4vbs2eMuXrzoamtr\n3fz58100GnWDg4P/+reSuEziBcGa6vc3iPWQ1a8A48WxgckV4FsJ5d96zmY9cCYgYBgBAAwjAIBh\nBAAwjAAAhhEAwDACABhGAADDCABgGAEADCMAgGEEADCMAACGEQDAMAIAGEYAAMMIAGAYAQAMIwCA\nYQQAMIwAAIYRAMAwAgAYRgAAwwgAYBgBAAwjAIBhBAAwjAAAhhEAwDACABhGAADDCABgGAEADCMA\ngGH3FYCOjg4tWLBA8+fP144dOyZqJgCTxHPOufH84fDwsMrLy/XVV19p3rx5eu6557R3714988wz\nf2/c8yZsUGQ2zrcSWcq39ZzNehj3HkBPT4+efvpplZaWqrCwUGvXrtXnn38+3s0BmALjDsD58+f1\n2GOPpa+XlJTo/PnzEzIUHlzd3d1TPcKY5Nu8YzXuAOTb7hByQ759oPJt3rEadwDmzZun/v7+9PX+\n/n6VlJRMyFAAJse4A1BTU6OffvpJfX19unnzpj777DM1NDRM5GwAgubuw4EDB1xZWZl76qmn3LZt\n2+66f+nSpU4SFy5cJvmydOnSrD7D4/4ZEED+40xAwDACABgWaABy/VTh9evXy/d9VVVVpW8bGBhQ\nNBpVWVmZ6urqNDQ0NIUTjtbf36/ly5ersrJSCxcu1O7duyXl7sw3btzQkiVLFA6HVVFRoa1bt0rK\n3Xmlv85wjUQiqq+vl5Tbs5aWlmrRokWKRCJavHixpLHPG1gAhoeH9dZbb6mjo0M//PCD9u7dqzNn\nzgT1dOOybt06dXR0jLotFospGo2qt7dXtbW1isViUzTd3QoLC/X+++/r+++/13fffacPPvhAZ86c\nydmZ//Of/6irq0snTpzQyZMn1dXVpSNHjuTsvJLU2tqqioqK9HkuuTyr53nq7u7W8ePH1dPTI2kc\n897PrwD3cvToUbdixYr09e3bt7vt27cH9XTjdvbsWbdw4cL09fLycpdMJp1zziUSCVdeXj5Vo2W0\natUq19nZmRczX7161dXU1LjTp0/n7Lz9/f2utrbWHTp0yL388svOudxeD6Wlpe7ChQujbhvrvIHt\nAeTrqcKpVEq+70uSfN9XKpWa4on+XV9fn44fP64lS5bk9MwjIyMKh8PyfT/99SVX5920aZN27typ\ngoK/Pxa5Oqv01x7ASy+9pJqaGn344YeSxj7vQ0EOl+88z8vJf8eVK1fU1NSk1tZWzZgxY9R9uTZz\nQUGBTpw4oUuXLmnFihXq6uoadX+uzLt//36FQiFFIpH/e/pvrsx6x7fffqvi4mL98ccfikajWrBg\nwaj7s5k3sD2AfD1V2Pd9JZNJSVIikVAoFJriiUa7deuWmpqa1NzcrMbGRkm5P7MkzZw5UytXrtSx\nY8dyct6jR4+qvb1dTzzxhF599VUdOnRIzc3NOTnrHcXFxZKkuXPnavXq1erp6RnzvIEFIF9PFW5o\naFA8HpckxePx9IcsFzjntGHDBlVUVGjjxo3p23N15gsXLqSPQl+/fl2dnZ2KRCI5Oe+2bdvU39+v\ns2fP6tNPP9WLL76oTz75JCdnlaRr167p8uXLkqSrV6/q4MGDqqqqGvu8QR2gcC7zqcJTbe3ata64\nuNgVFha6kpISt2fPHnfx4kVXW1vr5s+f76LRqBscHJzqMdO++eYb53meq66uduFw2IXDYffll1/m\n7MwnT550kUjEVVdXu6qqKvfee+8551zOzntHd3e3q6+vd87l7qw///yzq66udtXV1a6ysjL9+Rrr\nvJwKDBjGmYCAYQQAMIwAAIYRAMAwAgAYRgAAwwgAYBgBAAz7L+ZMyZQ+f/KQAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4e91590>"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# QUESTION make a circular pupil in a square array of 101 pixels that fills the array to the edge.\n",
      "# Call it pupil1 and show the pupil and add a colorbar\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# QUESTION make a circular pupil in a square array of 101 pixels filling the array with a central obscuration of 30% of the diameter\n",
      "# call it pupil2 and imshow it\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# QUESTION make a circular pupil in a square array of 101 pixels filling the array with a central obscuration of 20% of the diameter\n",
      "# and secondary support structures that cross through the centre of the secondary obscuration with a width of 5% of the diameter\n",
      "# imshow it and call it pupil3\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 61
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# to carry out a fourier transform, you need the complex amplitude of the wavefront at the pupil.\n",
      "# \n",
      "test_ampl = pupil1\n",
      "test_phase = np.ones_like(pupil1) * 0.0\n",
      "PupC = make_complex_amplitude(test_ampl,test_phase)\n",
      "test = displC(PupC)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# now we Fourier transform this and look at the resultant complex amplitude\n",
      "psf = fftshift(fft2(fftshift(PupC)))\n",
      "tmp = displC(psf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# it is hard to see the PSF core, so we can zoom in on a square box centered on the PSF\n",
      "(c2) = displC(psf, trim=51,minlog=-4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# cut out the central line of the intensity PSF\n",
      "modsq = c2*np.conj(c2)\n",
      "intensity = modsq.real\n",
      "intensity /= np.max(intensity)\n",
      "plot(intensity[:,25])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 67
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# you can't see much, so try the log10\n",
      "plot(np.log10(intensity[:,25]))\n",
      "# additionally you can add more line plots and they'll be displayed in additional colours\n",
      "plot(np.log10(intensity[:,32]))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# QUESTION: plot the cut across the PSF of pupil1 with that of pupil2, both normalised to their peak intensity and on a log10 scale.\n",
      "# What differences do you see?\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# QUESTION: image display the PSF for pupil1 and pupil3 with trim=101. \n",
      "# Describe the effect the secondary supports have on the PSF.\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 319
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}