make_Lyotstop_obs0.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import skimage.transform\n",
    "import astropy.io.fits as fits\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import PIL.ImageDraw\n",
    "import PIL.Image\n",
    "%pylab inline --no-import-all\n",
    "matplotlib.rcParams['image.origin'] = 'lower'\n",
    "matplotlib.rcParams['image.interpolation'] = 'nearest'\n",
    "matplotlib.rcParams['image.cmap'] = 'gray'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dictionary to map file keys from STScI to Krist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prim_key_map = {'hex1':'hex1', 'hex2':'hex2', 'hex3':'hex3', 'hex4':'hex4',\n",
    "                 'key24':'keystone24', 'pie08':'piewedge8', 'pie12':'piewedge12'}\n",
    "telap_dir = os.path.normpath('/astro/opticslab1/SCDA/Apertures/JPL/offset_masks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set basic parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "overwrite = False\n",
    "N = 250 # pupil array quadrant width after binning\n",
    "iD = 32\n",
    "oD = 90 # _LS_OD_map = {'hex1':76, 'hex2':82, 'hex3':81, 'hex4':82, 'pie08':90,\n",
    "aligntol = 2 # units of thousandths of pupil diameter\n",
    "L = 1024 # starting array width\n",
    "D = 2*N * (L//(2*N)) # starting pupil diameter (forced to be <= L)\n",
    "symm = 'half' # set to either 'quart' or 'half'\n",
    "shape = 'ann' # set to 'ann' or 'hex' or a primary mirror key\n",
    "#shape = 'hex4' # set to 'ann' or 'hex' or telap key\n",
    "edge = 'round' # set to 'gray' or 'round'\n",
    "write_LDZ = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Destination /astro/opticslab1/SCDA/Apertures/InputMasks/LS already exists\n"
     ]
    }
   ],
   "source": [
    "LS_dir = os.path.normpath(\"/astro/opticslab1/SCDA/Apertures/InputMasks/LS\")\n",
    "#LS_dir = os.path.expanduser(\"~/Box Sync/scda/InputMasks_wfirst/LS\")\n",
    "#LS_dir = os.path.normpath(\"/astro/opticslab1/SCDA/Apertures/InputMasks_WFIRST/LS\")\n",
    "\n",
    "if not os.path.exists(LS_dir):\n",
    "    os.mkdir(LS_dir)\n",
    "    print(\"created {:s} for LS arrays\".format(LS_dir))\n",
    "else:\n",
    "    print(\"Destination {:s} already exists\".format(LS_dir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define annular stop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inD = float(iD)/100\n",
    "outD = float(oD)/100\n",
    "xs = np.linspace(-L/2 + 0.5, L/2 - 0.5, L)\n",
    "[Xs, Ys] = np.meshgrid(xs, xs)\n",
    "\n",
    "if shape is 'hex':\n",
    "    xycent = (L/2-0.5,L/2-0.5)\n",
    "    hexagon_obj = matplotlib.patches.RegularPolygon(xycent,6,radius=outD*D/2/np.cos(np.pi/6))\n",
    "    hexagon_verts = hexagon_obj.get_verts()\n",
    "    hexagon_img = PIL.Image.new('L', (L,L), 0)\n",
    "    PIL.ImageDraw.Draw(hexagon_img).polygon(list(hexagon_verts.ravel()), outline=1, fill=1)\n",
    "    hexagon_mask = np.array(hexagon_img)\n",
    "    Lyot_stop = hexagon_mask\n",
    "    inside_ann = np.less_equal(Xs**2 + Ys**2, (inD*D/2)**2)\n",
    "    Lyot_stop[inside_ann] = 0.\n",
    "elif shape in prim_key_map: # Lyot stop with outer edge tracing the aperture perimeter \n",
    "    # Numerically erase the segment gaps\n",
    "    prim_fname = os.path.join(telap_dir, \"{0:s}_{1:04d}pix_offset_no_central_obsc.fits\".format(\n",
    "                                          prim_key_map[shape], D))\n",
    "    prim = fits.getdata(prim_fname)\n",
    "    max_shift = 1e-6\n",
    "    shift_range = np.linspace(-max_shift,max_shift,3)\n",
    "    [Xshifts, Yshifts] = np.meshgrid(shift_range, shift_range)\n",
    "    allowed_shifts = np.less_equal(Xshifts**2 + Yshifts**2, max_shift**2)\n",
    "    XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])\n",
    "    nogap_prim = np.zeros(prim.shape)\n",
    "    for (xshift,yshift) in XYshifts_allowed:\n",
    "        shift = skimage.transform.SimilarityTransform(translation=(xshift,yshift))\n",
    "        prim_shifted = skimage.transform.warp(prim, shift, order=1)\n",
    "        nogap_prim = np.logical_or(nogap_prim, prim_shifted)\n",
    "    perim = np.ndarray.astype(nogap_prim, np.float)\n",
    "    # Define the undersized stop outer edge based on a shift-and-multiply padding loop\n",
    "    w = D*(1. - outD)/2\n",
    "    XYshifts_allowed = [(-w,0), (+w,0),\n",
    "                        (+w*np.sin(np.pi/6),+w*np.cos(np.pi/6)), (-w*np.sin(np.pi/6),+w*np.cos(np.pi/6)),\n",
    "                        (+w*np.sin(np.pi/6),-w*np.cos(np.pi/6)), (-w*np.sin(np.pi/6),-w*np.cos(np.pi/6))]\n",
    "    padded_perim = np.ones(perim.shape)\n",
    "    for (xshift,yshift) in XYshifts_allowed:\n",
    "        shift = skimage.transform.SimilarityTransform(translation=(xshift,yshift))\n",
    "        perim_shifted = skimage.transform.warp(perim, shift, order=1)\n",
    "        padded_perim *= perim_shifted\n",
    "    Lyot_stop = padded_perim.copy()\n",
    "    # define inner, circular edge\n",
    "    inside_ann = np.less_equal(Xs**2 + Ys**2, (inD*D/2)**2)\n",
    "    Lyot_stop[inside_ann] = 0.\n",
    "    plt.figure(figsize=(8,8))\n",
    "    plt.imshow(Lyot_stop - perim)\n",
    "    plt.figure(figsize=(8,8))\n",
    "    plt.imshow(Lyot_stop)\n",
    "else:\n",
    "    inside_ann = np.less_equal(Xs**2 + Ys**2, (inD*D/2)**2)\n",
    "    outside_ann = np.greater_equal(Xs**2 + Ys**2, (outD*D/2)**2)\n",
    "    Lyot_stop = np.ones((L,L))\n",
    "    Lyot_stop[inside_ann] = 0.\n",
    "    Lyot_stop[outside_ann] = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bin to final array resolution, crop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#scalefac = float(N)/(D/2)\n",
    "#Lyot_stop_binned = skimage.transform.rescale(Lyot_stop, scalefac, order=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binning the original LS array 2x\n"
     ]
    }
   ],
   "source": [
    "N_orig = D/2\n",
    "Lyot_stop_crop = Lyot_stop[L/2-N_orig:L/2+N_orig,L/2-N_orig:L/2+N_orig]\n",
    "scalefac = int(N_orig/N)\n",
    "print(\"Binning the original LS array {0:d}x\".format(scalefac))\n",
    "Lyot_stop_binned = np.reshape(Lyot_stop_crop, (Lyot_stop_crop.shape[0]/scalefac, scalefac, \n",
    "                                               Lyot_stop_crop.shape[1]/scalefac, scalefac)).mean(1).mean(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 1000)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lyot_stop_crop.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 500)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lyot_stop_binned.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(500, 250)\n"
     ]
    }
   ],
   "source": [
    "L_bin = Lyot_stop_binned.shape[0]\n",
    "if symm is 'half':\n",
    "    Lyot_stop_binned_crop = Lyot_stop_binned[L_bin/2-N:L_bin/2+N,L_bin/2:L_bin/2+N]\n",
    "    Lyot_stop_binned_crop_binary = np.round(Lyot_stop_binned_crop).astype(int)\n",
    "    print Lyot_stop_binned_crop.shape\n",
    "else:\n",
    "    Lyot_stop_binned_crop = Lyot_stop_binned[L_bin/2:L_bin/2+N,L_bin/2:L_bin/2+N]\n",
    "    Lyot_stop_binned_crop_binary = np.round(Lyot_stop_binned_crop).astype(int)\n",
    "    print Lyot_stop_binned_crop.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAMuCAYAAACaYU1JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X+MrXl9H/b3x8Gsg9MFNQm7dVNkIhrcqJHTuwSKHJiZ\nYsXFSMSRK5eVJWKjyrICkbVSGxoVaWagqlpHgS0JWPxB7VhJWtF1EY4VQwiZObUxBulC7LrGtMRQ\nQshug03ZFdheA9/+cc65nJ29d+78OOc8z/k+r5c0Yu85z979ztHc++bzeT7f71OttQBAT75l6AUA\nwLoJNwC6I9wA6I5wA6A7wg2A7gg3ALoj3ADozjOGXsBFVNUfT/J9ST6b5PeHXQ0AA/q2JN+Z5AOt\ntd+500U7EW6ZB9s/GHoRAIzGDyf5h3d6c1fakp8degEAjMpnz3tzV8JNKxKAVefmwq6EGwBcmHAD\noDvCDYDuCDcAuiPcAOiOcAOgO8INgO4INwC6I9wA6I5wA6A7wg2A7gg3ALoj3ADojnADoDvCDYDu\nCDcAuiPcAOiOcAOgO8INgO4INwC6I9wA6I5wA6A7wg2A7gg3ALoj3ADojnADoDvCDYDuCDcAuiPc\nAOiOcAOgO8INgO4INwC6I9wA6I5wA6A7wg2A7jxj6AVcxs2bN3Pjxo0r/buz2Synp6cXvn557WX+\nHQDGoVprQ6/hrqrqRpKb1wm3dbtbWB4dHW1tLQAT9EBr7eN3elNbEoDuqNy2ZDabJXl6m1OFB3Al\n51Zuwm2Ejo+Pn/aaEAR4Cm1JAKZF5baDbjfMorIDJkblBsC0qNw6tFrZnZ6e2qsH9MhACXNnQ2/1\nfwF2jLYkANOicuMpe/BUdMCO0JbkalZDzzQmMDLakgBMi8qNKzk+PtbCBIakLcnmLY8Ms/UA2BJt\nSQCmReXGxqzuqzOQAqyZtiTj4V4dsCbakgBMi8qNQRlEAa5I5QbAtKjcGJ3j42MDKMDdqNzYLYeH\nh2mtpbV2q125v78/9LKAHSLcAOiOcGPU9vb2sre3l5OTk1vVnJYlcDfCjZ1ztm25v7+vbQk8hXAD\noDvPGHoBcB3LluXSct+c1iVMm8qNrhweHj6lbSnkYJqEGwDdEW50zfAJTJN7bkzG6v252Wx2q2Xp\nTEvoj8oNgO6o3JgkVRz0TbgxeWeDLplvJRB0sLu0JQHojnCDFatnWXoiAewu4QZAd4Qb3IEqDnaX\ncIMLEHSwW4QbAN2xFQAuaW9vL0lycnJya+uAKg7GReUG17BsV3oKAYyLcAOgO8IN1siz5GAchBts\nyOHhoclKGIhwA6A7piVhg0xWwjCEG2zJMuhaazk+Pk4S9+VgQ7QlAeiOyg0GcHh4mOSbLUrPj4P1\nUrkB0B2VGwzIwAlshnCDkVgdODk4OEgSrUq4Im1JALqjcoMROjk5SZLMZjNtSrgClRuM2OoTB4Qc\nXJxwA6A72pKwI0xUwsWp3GCHnH04qpCD2xNuAHRHWxJ2mFYl3J5wgx1n8zc8nbYkAN1RuUFHbP6G\nOZUbAN0RbtChvb29nJ6e5vT0NPv7+6o4JkdbEjq1+jidJDk4ODBowmSo3ADojsoNJsKeOKZE5QYT\nsjy+y704eifcAOiOcIMJ2tvby8nJSU5OTlRvdMk9N5i4k5OTHB8fJ0mOjo6GXQysicoNgO6o3IAc\nHh4miSETuqFyA25xsgm9EG4AdEdbEngKx3bRA5UbAN1RuQHnslWAXSTcgLsyTcmu0ZYEoDvCDbiw\n1a0CMGbCDbgUTxZgFwg3ALoj3IAr8WQBxky4Ade2DDghx1gINwC6Y58bsBarx3UlMVHJoFRuwFq5\nD8cYCDcAuqMtCWyEMykZknADNsaZlAxFWxKA7gg3YOOcScm2CTcAuiPcgK1YPXAZNk24AVulRck2\nCDcAuiPcgK3TomTThBswGC1KNkW4AdAdJ5QAg9rb20vyzacIOMmEdVC5AaPgPhzrJNwA6I5wA0bF\nkAnrINyA0dGi5LqEGwDdEW7AaO3t7eXo6MjDTrm0a4VbVf3XVfWNqnrrmdffXFVfqKqvVtUHq+oF\nZ96/p6reUVVfrKonquqRqnruddYCAEtXDreq+gtJfizJr515/Y1J3rB478VJvpLkA1X1zJXLHk7y\nqiQ/mOTlSb4jyc9ddS1Avw4PD3N4eGj/G5dypXCrqj+W5O8n+S+S/H9n3v6JJG9prf1Ca+03krw2\n8/D6gcW/e2+S1yV5qLU2a619IsmPJvmeqnrx1b4NoHcnJyfZ398XclzIVSu3dyT5R621f7b6YlU9\nP8n9ST60fK219niSjyZ56eKlF2V+MsrqNZ9K8rmVawDgyi59/FZVvSbJn888pM66P0lL8tiZ1x9b\nvJck9yV5chF6d7oG4GlOTk6SJAcHB7YJcK5LhVtV/anM75d9b2vtDzezJIDznZyc5ODgIEmEHLd1\n2crtgSR/MsnHq6oWr/2RJC+vqjck+a4klXl1tlq93ZfkE4t/fjTJM6vq3jPV232L9+7ooYceyrOf\n/eynvPbggw/mwQcfvOS3AUDPLhtu/zTJnzvz2s8k+WSS/7619ttV9WiSVyT59eTWAMlLMr9PlyQ3\nk3xtcc17F9e8MMnzknzkvP/42972tty4ceOSSwZ6tNqiTFRwPNWlwq219pUkv7n6WlV9JcnvtNY+\nuXjp4SRvqqpPJ/lskrck+XyS9y1+j8er6t1J3lpVX0ryRJK3J/lwa+1j1/hegAlabvA2RcmqdTzP\nrT3lF639ZFU9K8m7kjwnyS8leWVr7cmVyx5K8vUkjyS5J8n7k7x+DWsBgFRr7e5XDayqbiS5efPm\nTW1J4LZms5nqbVoeaK19/E5vOlsS6IJH5bBKuAHQHeEGdGP5HDhPEUC4AdAd4QZ0Z/kUAQMm0yXc\ngC6dnJzcepIA0yPcAOiOcAO6ZrhkmoQb0DX736ZJuAHQHeEGdM/+t+kRbsBk2CIwHcINgO6s45E3\nADvDQ06nQbgBk+Qhp33TlgSgO8INmKTlBKW2ZJ+EGzBpe3t7Jig7JNwA6I5wAybPEwT6I9wA6I5w\nA1hwPFc/hBvAgicI9EO4AdAd4QawYrn/zXDJbhNuALexPIOS3STcAOiOcAO4A8Mlu0u4AdyB+2+7\nS7gB0B3hBnAXhkt2j3ADuAD333aLcAOgO88YegEAu2Bvby9Jbg2XqOLGTeUGQHdUbgCXsBwuqaqB\nV8J5VG4AV6AtOW7CDYDuCDeAK1ieXOL0knFyzw3gitx/Gy+VGwDdEW4A12S4ZHyEG8A1uf82PsIN\ngO4YKAFYg6OjoyRRvY2Eyg1gDZYPNl2GHMMSbgB0R7gBrNHh4eHQSyDCDYAOCTeANfPU7uEJN4A1\nWw6XmJwcjnADoDvCDWBDbAsYjnAD2JDlvjcht33CDYDuVGtt6DXcVVXdSHLz5s2buXHjxtDLAbg0\nz3xbuwdaax+/05sqN4At0J7cLuEGQHe0JQG2SHtybbQlAcbCySXbIdwA6I5wA9ii5bFcjubaLE/i\nBtgyT+3ePJUbAN0RbgBb5qkBmyfcAOiOcAMYiBNLNke4AQzE5OTmCDcAumMrAMCAbAvYDJUbwIBM\nTm6GcAOgO8INYARMTq6XcAMYAZOT6yXcAOiOaUmAkTA5uT4qN4CRMDm5PsINgO4IN4CRMTl5fcIN\ngO4IN4CRsS3g+kxLAoyQycnrUbkB0B2VG8AI7e3tDb2EnaZyAxixo6Mj05NXINwA6E611oZew11V\n1Y0kN2/evJkbN24MvRyArauqoZcwNg+01j5+pzdVbgA7QGvycoQbAN0RbgA74PDwcOgl7BThBrAj\nnFpyccINgO7YxA2wIxzJdXEqNwC6o3ID2BGO5Lo4lRvAjrHn7e6EGwDdEW4AO8aet7sTbgA7yJ63\n8wk3ALpjWhJgB9nzdj7hBrCDbAs4n7YkAN3ZqcrtgQceOPf9ZXl+pzLdhBHQm6OjI/vebmOnnsS9\nid979Ydif39fqQ/snIk+pduTuAGYlslXbnezWtlpawJjdHBwkCQ5PT0ddiHbdW7lJtwuafW+nrAD\nxuD4+DjJ5M6c1JYEYFpUbmuwegyOag4YysQGS1RuAEyLym1DVHPAth0cHExpqMRAyRgsb/QKOmBT\njo+PpzRUoi0JwLSo3AawbFceHR05EQVYqwkNlWhLjtlq0CVO+gauZ0J73rQlAZiWnXoqQI+Wk02r\nJ5+o4oCr8vDSOW3JETs6OjJdCVzJBO69aUsCMC3CbcSOjo5SVamqzGazoZcD7JDVgySmSLjtiP39\n/RwcHOTg4EDQAXcl3ACgMwZKdpTjvICL6HiwxEAJANOicttxjvICztPxkwJUbj07PT3N6elp9vf3\nM5vNDJsATzHVoRLhBkB3tCU7tL+/n5OTk6GXAYxEp0Ml2pJTc3p6emtPHMAUCTcAuuOpAJ1aTkdV\n1a1/Nk0J07TcFzuBZ7zdItwmYDktdXp6KuBggqY4MaktCUB3TEtOjGO7YLo6m5o8d1pSuE3Usk1h\nywBMR2enldgKAMC0CLeJWh7bZS8cTMeUBkuEGwDdEW4T5zQTmI4pVW4GSrjFmZTQv44mJg2UADAt\nwo1btCihf/v7+5NoTwo3nmI5RXl8fDz0UoANEG63UVU/XlW/VlVfXnz9SlX9p2eueXNVfaGqvlpV\nH6yqF5x5/56qekdVfbGqnqiqR6rquev4ZgAgueRASVW9KsnXk/zfSSrJjyT5r5L8+dbaJ6vqjUne\nmOS1ST6b5L9N8ueS/AettScXv8dPJXllkr+a5PEk70jy9dbay8757xooGYCnCUB/ZrNZki4mJzd7\n/FZV/U6S/7K19tNV9YUkf6u19rbFe/cmeSzJX22tvWfx63+T5DWttfcurnlhkk8m+Y9bax+7w39D\nuA3I0wSgPx1MTW5mWrKqvqWqXpPkWUl+paqen+T+JB9aXtNaezzJR5O8dPHSizJ/zM7qNZ9K8rmV\nawDgWi4dblX1H1bVE0n+IMk7k/yVRUDdn6RlXqmtemzxXpLcl+TJRejd6RpGpoP2BXBG73+ur1K5\n/VaS707y4iQ/leRnq+q71roqRscWAeiLcDujtfa11tpvt9Y+0Vr7b5L8WpKfSPJo5kMm9535V+5b\nvJfF/z5zce/tTtcAwLU8Yw2/x7ckuae19pmqejTJK5L8enJroOQlmU9EJvOhkK8trlkdKHleko+s\nYS1syHJyclm9OaYLdlvvldulwq2q/rskv5j5AMi/leSHk+wl+UuLSx5O8qaq+nTmWwHekuTzSd6X\nzAdMqurdSd5aVV9K8kSStyf58J0mJQHgsi5buT03yd9L8u8k+XLmFdpfaq39syRprf1kVT0rybuS\nPCfJLyV55XKP28JDme+VeyTJPUnen+T11/km2J5lBTebzWwPgB3W+59fTwXgymzyht2243vdPBUA\ngGlZx0AJE3V0dJTEcAnsquWf4eX/9kTlxpV5ggAwVsINgO5oS3JtR0dHt/bMGC6B3dHzXjfTkqzF\n8g+J+2+we3Z0atK0JADToi3JWtjcDYyJyo216rmHD72yFQAAdoBwY+2Oj4/tfQMG5Z4ba7dscezv\n77v/BgxC5QZAd+xzY2P29/fte4MdsYN73exzA2BahBsbc3p6mtlsltlsNvRSgIkxUMJGeSwOMASV\nGwDdUbmxUavHciWeGgBj1duDS4UbW6E9CWyTtiQA3VG5sRXL9iTANqjc2CpnTsI47e/vd/VUD+EG\nQHccv8XW7cLPHEzVDh3D5fgtxsUjcYBNE24AdEe4sXWnp6emJ4GNEm4AdMc+N7ZO1Qbj1csxXCo3\nBmOoBNgU4QZAd4Qbg9GeBDZFuDEY4QZsinADoDumJRnUcqjk8PBw4JUASbo5PNnZkgxq+QfJQ0xh\nXHbgjElnSwIwLdqSDMpQCbAJKjdGwYZuYJ2EGwDdEW4AdEe4AdAd4cYo7PoJ5NCbXd/vJtwAeBrh\nBgAjI9wYjdlsltlsNvQygA7YxM1oLDd07+3tDbsQYOep3ADojnADoDvakozGcjuAx98A16VyA6A7\nwg2A7gg3ALoj3ADojnBjdDzbDYbn+C0AurPrhykINwC6I9wA6I5wA6A7wg2A7gg3Rmf5dACAqxJu\nAHRHuDE6KjfguoQbAN0RbgB0R7gB0B3hBkB3hBsA3RFuAHRHuAHQHeEGQHeEGwDdEW4AdEe4AdAd\n4QZAd4QbAN0RbgB0R7gB0B3hBkB3hBsA3RFujM7+/v7QSwB2nHADoDvCjdFRuQHXJdwA6I5wA6A7\nwg2A7gg3ALoj3AB4mtlsNvQSrkW4MTqHh4dDLwEm7/T0dOglXItwA6A7wg2A7gg3ALoj3ADozjOG\nXgAsHR0dDb0EoBMqNwC6I9wA6I62JKPhaQDAugg3RmNvb2/oJQCd0JYEoDvCDYCncfwWrIFtADAu\nwg0ARka4AdAd4QZAd6q1NvQa7qqqbiS5OfQ62Jxd+DmEKamqoZdwNw+01j5+pzftc2NQNm4Dm6At\nCUB3VG4MSuUG4zKbzYZewloINwZ1eHg49BKAFbu+v21JWxKA7gg3BqMlCWyKcGMwwg3YFOEGQHeE\nG4MxTAJsimlJtk47Esarlyd0qNwA6I7Kja1TuQGbpnIDoDueCsDW7cLPHEzVDjwNYMlTARiPXm5W\nQ296OVNySVsSgO6o3Ngqe9tgnHo5MHlJuLEVJiSBbdKWBKA7Kje2wiAJsE22ArBRy3bkycnJsAsB\nzrVDWwCWzt0KoC0JQHe0Jdko7UhgCMKNjdnf38/e3t7QywAmSFsSgO4INzZGSxJ2w/Hx8dBLWDvh\nBkB33HNj7ZYVm/ttwFCEG2vn/EhgaNqSAHRHuLFWvZ0sDlPQ4/CXtiRrsTxmy302YAxUbgB0R+XG\nWvTY1oDezWazoZewMcKNazs6OtKOhB3U8z1ybUkAuqNy48qWQyT2tQFj42GlXNmypaElCbtpBx9Q\nump9Dyutqr9ZVR+rqser6rGqem9V/ZnbXPfmqvpCVX21qj5YVS848/49VfWOqvpiVT1RVY9U1XMv\nsxYAuJPL3nN7WZK/k+QlSb43ybcm+SdV9UeXF1TVG5O8IcmPJXlxkq8k+UBVPXPl93k4yauS/GCS\nlyf5jiQ/d8XvgQGcnp5mb29P1QaM0qXuubXWvn/111X1I0n+3yQPJPnlxcs/keQtrbVfWFzz2iSP\nJfmBJO+pqnuTvC7Ja1prs8U1P5rkk1X14tbax67+7bBpNmtDH3reBpBcf1ryOUlakt9Nkqp6fpL7\nk3xoeUFr7fEkH03y0sVLL8o8VFev+VSSz61cAwBXduVpyZrfiXw4yS+31n5z8fL9mYfdY2cuf2zx\nXpLcl+TJRejd6RpGZlmxnZycDLsQYC163uOWXG8rwDuT/Nkk37OmtQDAWlwp3Krq7yb5/iQva639\n65W3Hk1SmVdnq9XbfUk+sXLNM6vq3jPV232L9xghFRv0pffK7dL33BbB9peTHLTWPrf6XmvtM5kH\n1CtWrr838+nKX1m8dDPJ185c88Ikz0vykcuuh83r/Q8BTFHvf64vVblV1TuTPJjk1Um+UlX3Ld76\ncmvt9xf//HCSN1XVp5N8Nslbknw+yfuS+YBJVb07yVur6ktJnkjy9iQfNikJwDpcti3545kPjJye\nef1Hk/xskrTWfrKqnpXkXZlPU/5Skle21p5cuf6hJF9P8kiSe5K8P8nrL7t4NssJJNCf3rcALF12\nn9uF2pittaMkR+e8/wdJ/vriixFy0j/0qfd25JKnAgDQHU8F4Cmc9A99m0rlJty4ZX9/38g/dG4q\n4aYtCUB3VG44WgsmYiqTkolwmzytSJiOqbQkE21JADqkcpsorUiYHpUbAOwwldvEHB0dJbGPDaZo\nSpWbcJuQ09NTR2rBBE1pSnJJWxKA7qjcJsDp/jBtU2pHLlVrbeg13FVV3cj8IadckGlIYKmqhl7C\nJjzQWvv4nd7UlgSgO9qSHXLqCDB1wq0j7q0BZx0fHw+9hEFoSwLQHZXbjlsOjhwdHanYgKeZ4qRk\nItx2lpNGgIuYarhpSwLQHZXbDtGCBC5qqoMkSyo3ALqjctsRDj0GLmOq99qWhNuIHR0dGRgBrmTq\n4aYtCUB3VG4js7+/f2vMXxsSuKwpPrvtdoTbwFYnIBOBBlzP1NuRS9qSAHTH89wGYL8asCmdPrvt\nds59npu25JY4LgvYtKlv3F6lLQlAd1RuG7K/v3+r/ahaA7bBMMk3Cbc1EGTAGAi3b9KWBKA7KrdL\nWlZo+/v7qjRgFAySPJ3KDYDuqNzuYjnCn7ifBoyTe21PN/lwWw2v/f19m6qBnSPcnk5bEoDudFW5\nrQ573I62ItAbwyS3t1NnS968eTM3btwYejkAozGhsyTPOvdsSW1JALrTVVsSYCo8lPR8wg1gB61O\nevN02pIAdEflBrCD7G07n8oNYMcY/7874QZAd4QbwI4xTHJ37rkB7Ajj/xencgOgOyo3gB2hHXlx\nKjcAuqNyA9gR9rZdnMoNYAfY23Y5wg2A7gg3gB1gmORy3HMDGDHtyKtRuQHQHZUbwIhpR16NcAMY\nIUdtXY+2JADdUbkBjJB25PUIN4CRmc1mTiO5Jm1JALoj3ABGRkvy+oQbAN1xzw1gJJbj/+63XZ9w\nAxgJ7cj10ZYEoDsqN4ARMP6/Xio3gBHQklwv4QZAd7QlAQZkQnIzhBvAgLQjN0NbEoDuqNwABmJC\ncnNUbgAD0ZLcHOEGQHe0JQG2zITk5qncAOiOyg1gy9xr2zzhBrBFJiS3Q1sSgO4IN4At2t/fH3oJ\nk6AtCbAFx8fHQy9hUlRuAHRH5QawBSYkt0u4AWyQduQwtCUB6I7KDWBDZrOZduRAVG4AGyLYhiPc\nAOiOtiTAmjn1f3jCDWDNnEIyPG1JALoj3ADWyL62cRBuAHTHPTeANVgOkRj/HwfhBrAGQm1ctCUB\n6I7KDeCaZrOZPW0jo3IDuCb72sZHuAHQHW1JgCs6ODgYegncgXADuAL32cZNWxKA7gg3gCswRDJu\n2pIAl+A+225QuQHQHZUbwAV4AOluUbkB0B2VG8AFGCDZLSo3gLswRLJ7hBsA3dGWBLgDQyS7S7gB\n3IH7bLtLWxKA7gg3gNswRLLbtCUBVrjP1geVGwDdUbkBLMxmM0MknVC5ASwcHR0NvQTWRLgB0B1t\nSWDylpORhkj6oXIDoDsqN2DSZrOZiq1Dwg2YpOV+NtORfdKWBKA7Kjdgkoz99024AZNiMnIatCUB\n6I7KDZiM4+NjFdtECDege8vJSPfZpkNbEoDuqNyArjnpf5pUbkDXtCKnSbgB0B1tSaBL9rNNm3AD\numPkH21JALqjcgO6YT8bSyo3ALpz6XCrqpdV1c9X1b+qqm9U1atvc82bq+oLVfXVqvpgVb3gzPv3\nVNU7quqLVfVEVT1SVc+9zjcCTNtyP5s9bSRXq9y+Pck/T/LXkrSzb1bVG5O8IcmPJXlxkq8k+UBV\nPXPlsoeTvCrJDyZ5eZLvSPJzV1gLMHGz2cxGbZ7m0vfcWmvvT/L+JKmqus0lP5HkLa21X1hc89ok\njyX5gSTvqap7k7wuyWtaa7PFNT+a5JNV9eLW2seu9J0AwMJaB0qq6vlJ7k/yoeVrrbXHq+qjSV6a\n5D1JXrT4765e86mq+tziGuEGXJjhEW5n3dOS92feqnzszOuPLd5LkvuSPNlae/ycawDOZZM25zEt\nCUB31l25PZqkMq/OVqu3+5J8YuWaZ1bVvWeqt/sW793RQw89lGc/+9lPee3BBx/Mgw8+eN11Azvk\n4OBAxca51hpurbXPVNWjSV6R5NeTZDFA8pIk71hcdjPJ1xbXvHdxzQuTPC/JR877/d/2trflxo0b\n61wysEO0IrmoS4dbVX17khdkXqElyZ+uqu9O8ruttX+Z+Zj/m6rq00k+m+QtST6f5H3JrQGTdyd5\na1V9KckTSd6e5MMmJQFYh6tUbi9KcpL54EhL8rcXr/+9JK9rrf1kVT0rybuSPCfJLyV5ZWvtyZXf\n46EkX0/ySJJ7Mt9a8PorfQfAJGhFchlX2ec2y10GUVprR0mOznn/D5L89cUXwB0dHx8n0YrkckxL\nAtAdTwUARms2m9mkzZWo3ADojsoNGJ3lc9kchsxVCTdgVJzwzzpoSwLQHZUbMApakayTcAMGJdTY\nBG1JALqjcgMGY3iETRFuwNZpRbJp2pIAdEflBmyVViTbINyArdCKZJu0JQHojsoN2DitSLZN5QZA\nd1RuwMYsn6LtmWxsm3ADNuLg4CCnp6dDL4OJ0pYEoDsqN2CtDg4OkkTVxqCEG7AWQo0x0ZYEoDsq\nN+DaDI8wNsINuJLZbHZrxF+wMTbakgB0R+UGXIoDkNkFwg24MGdEsiu0JQHojsoNuCtnRLJrhBtw\nLmP+7CJtSQC6o3IDnmI5DWkPG7tM5QZAd1RuwC1G/emFcANMQ9IdbUkAuqNyg4kz6k+PhBtMkBP9\n6Z22JADdUbnBhDjRn6kQbjAR7q0xJdqSAHRH5QadcowWUybcoENOGmHqtCUB6I7KDTpycHCQRAsS\nVG4AdEflBjvO3jV4OuEGO8zeNbg9bUkAuqNygx2iBQkXI9xgR2hBwsVpSwLQHZUbjJiTRuBqhBuM\nkM3YcD3akgB0R+UGI2ESEtZHuMGAVh9LowUJ66MtCUB3VG4wgOPj4yTffJAosF4qNwC6o3KDLTEw\nAtsj3GCDDIzAMLQlAeiOyg025Pj42MAIDES4wRoJNBgHbUkAuqNyg2swAQnjJNzgkkxAwvhpSwLQ\nHZUbXIBqDXaLcIM7EGiwu7QlAeiOyg1WqNagDyo3ALqjcmPyZrPZrVNFVGvQB+HGJAk06Ju2JADd\nUbkxGao1mA7hRtec0g/TpC0JQHdUbnTl+Pg4SVRrMHHCjZ3mPhpwO9qSAHRH5cbOMSQC3I1wY9Sc\n9QhchbYkAN1RuTE62o7AdQk3BrUc3T89PdV2BNZGWxKA7qjc2Krj4+NbFZpKDdgUlRsA3VG5sTGz\n2exWdWZABNgm4cZaGAwBxkRbEoDuqNy4EoMhwJgJN+5oefTV6empe2bATtGWBKA7KjeeUqFpNQI9\nEG4TsjoJ1/WDAAAGN0lEQVSaL8SAnmlLAtAdlVuHzlZoqjNgaoTbDloNryXTjADfpC0JQHdUbiO0\nPMpqlcoM4OJUbgB0R+W2Jat7yVapyADWT7hd0e2GOlYJLYDhaEsC0J3JVG53q7TOcoIHwO7aqXB7\n4IEHhl4CADtAWxKA7gg3ALoj3ADojnADoDvCDYDuCDcAuiPcAOiOcAOgO8INgO4INwC6I9wA6I5w\nA6A7wg2A7gg3ALoj3ADojnADoDvCDYDuCDcAuiPcAOiOcAOgO8INgO4INwC6I9wA6I5wA6A7wg2A\n7gg3ALoj3ADojnADoDvCDYDuCDcAuiPcAOiOcAOgO4OGW1W9vqo+U1W/V1W/WlV/Ycj1ANCHwcKt\nqv7zJH87yWGS/yjJryX5QFX9iaHWBEAfhqzcHkryrtbaz7bWfivJjyf5apLXDbgmADowSLhV1bcm\neSDJh5avtdZakn+a5KVDrAmAfgxVuf2JJH8kyWNnXn8syf3bXw4APdmVaclvG3oBAIzKubkwVLh9\nMcnXk9x35vX7kjx6m+u/c9MLAmCnfOd5bz5jS4t4itbaH1bVzSSvSPLzSVJVtfj122/zr3wgyQ8n\n+WyS39/SMgEYn2/LPNg+cN5FNZ/j2L6q+qEkP5P5lOTHMp+e/M+SfFdr7d8MsigAujBI5ZYkrbX3\nLPa0vTnzduQ/T/J9gg2A6xqscgOATdmVaUkAuDDhBkB3hBsA3dmJcPP0gPNV1WFVfePM12+euebN\nVfWFqvpqVX2wql4w1HqHUFUvq6qfr6p/tfh8Xn2ba879jKrqnqp6R1V9saqeqKpHquq52/sutu9u\nn1tV/fRtfvb+8ZlrJvW5VdXfrKqPVdXjVfVYVb23qv7Mba7z87ZBow83Tw+4sN/IfOr0/sXXX1y+\nUVVvTPKGJD+W5MVJvpL5Z/jMAdY5lG/PfCL3ryV52hTVBT+jh5O8KskPJnl5ku9I8nObXfbgzv3c\nFn4xT/3Ze/DM+1P73F6W5O8keUmS703yrUn+SVX90eUFft62oLU26q8kv5rkf1z5dSX5fJK/MfTa\nxvKVefB//Jz3v5DkoZVf35vk95L80NBrH+jz+kaSV1/mM1r8+g+S/JWVa164+L1ePPT3NODn9tNJ\n/rdz/h2f2/ws3W8k+Ysrr/l52/DXqCs3Tw+4lH9/0Tr6F1X196vq30uSqnp+5v9vevUzfDzJR+Mz\nTHLhz+hFme8LXb3mU0k+F5/j/qL99ltV9c6q+rdX3nsgPrfnZF71/m7i521bRh1u8fSAi/rVJD+S\n5PsyP/Hl+Un+96r69sw/pxaf4Xku8hndl+TJxV9Cd7pmin4xyWuT/CdJ/kaSvST/eHGcXjL/bCb7\nuS0+h4eT/HJrbXkf3M/bFgx2Qgnr01pbPWPtN6rqY0n+nyQ/lOS3hlkVU9Bae8/KL//Pqvo/kvyL\nJPtJTgZZ1Li8M8mfTfI9Qy9kasZeuV326QEkaa19Ocn/leQFmX9OFZ/heS7yGT2a5JlVde8510xe\na+0zmf+5XU7+TfZzq6q/m+T7k+y31v71ylt+3rZg1OHWWvvDJMunByR5ytMDfmWodY1dVf2xzP9y\n+cLiL5tH89TP8N7MJ7l8hrn1F/LdPqObSb525poXJnleko9sbbEjV1V/KskfT7L8y3ySn9si2P5y\nkoPW2udW3/Pzth270JZ8a5KfWTwiZ/n0gGdl/kQBklTV30ryjzJvRf67SY6T/GGS/2VxycNJ3lRV\nn878sUFvyXzi9H1bX+xAFvcfX5D5/2NOkj9dVd+d5Hdba/8yd/mMWmuPV9W7k7y1qr6U5InMH8/0\n4dbax7b6zWzReZ/b4usw8/H0RxfX/Q+Zdw0+kEzzc6uqd2a+HeLVSb5SVcsK7cutteUju/y8bdrQ\n45oX+cp8j81nMx+V/UiSFw29pjF9JfmfM/+D8XuZT1P9wyTPP3PNUebjx1/N/C+eFwy97i1/RnuZ\nj1F//czX/3TRzyjJPZnvX/pi5n/Z/K9Jnjv09zbU55b5c7Xen3mw/X6S307yU0n+5JQ/tzt8Xl9P\n8toz1/l52+CXpwIA0J1R33MDgKsQbgB0R7gB0B3hBkB3hBsA3RFuAHRHuAHQHeEGQHeEGwDdEW4A\ndEe4AdCd/x97BDnxx7IuAAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110a04790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if symm is 'half':\n",
    "    plt.figure(figsize=(20,10))\n",
    "    #plt.imshow(Lyot_stop_binned_crop_binary)\n",
    "    plt.imshow(Lyot_stop_binned_crop)\n",
    "else:\n",
    "    plt.figure(figsize=(10,10))\n",
    "    plt.imshow(Lyot_stop_binned_crop)\n",
    "    #plt.imshow(Lyot_stop_binned_crop_binary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write the binned, cropped Lyot stop to a dat (ASCII) file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote binned, cropped LS array to /astro/opticslab1/SCDA/Apertures/InputMasks/LS/LS_half_ann32D90_clear_N0250.dat\n"
     ]
    }
   ],
   "source": [
    "if symm is 'half':\n",
    "    if shape in ['ann', 'hex']:\n",
    "        LS_dat_fname_tail = \"LS_half_{0:s}{1:02d}D{2:02d}_clear_N{3:04d}.dat\".format(\n",
    "                             shape, int(round(100*inD)), int(round(100*outD)), N)\n",
    "    else:\n",
    "        LS_dat_fname_tail = \"LS_half_{0:s}P{1:02d}D{2:02d}_clear_N{3:04d}.dat\".format(\n",
    "                             shape, int(round(100*inD)), int(round(100*outD)), N)\n",
    "else:\n",
    "    if shape in ['ann', 'hex']:\n",
    "        LS_dat_fname_tail = \"LS_quart_{0:s}{1:02d}D{2:02d}_clear_N{3:04d}.dat\".format(\n",
    "                             shape, int(round(100*inD)), int(round(100*outD)), N)\n",
    "    else:\n",
    "        LS_dat_fname_tail = \"LS_quart_{0:s}P{1:02d}D{2:02d}_clear_N{3:04d}.dat\".format(\n",
    "                             shape, int(round(100*inD)), int(round(100*outD)), N)\n",
    "        \n",
    "LS_dat_fname = os.path.join(LS_dir, LS_dat_fname_tail)\n",
    "\n",
    "if not os.path.exists(LS_dat_fname) or overwrite==True:\n",
    "    if edge is 'gray':\n",
    "        np.savetxt(LS_dat_fname, Lyot_stop_binned_crop, fmt='%.6f', delimiter=\" \")\n",
    "    else: # store version rounded to binary integers\n",
    "        np.savetxt(LS_dat_fname, Lyot_stop_binned_crop_binary, fmt='%d', delimiter=\" \")\n",
    "    print(\"Wrote binned, cropped LS array to {0:s}\".format(LS_dat_fname))\n",
    "else:\n",
    "    print(\"LS array {0:s} already exists, will not overwrite\".format(LS_dat_fname))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#if os.path.exists(LS_dat_fname): os.remove(LS_dat_fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Lyot plane dark zone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lyot plane dark zone width in binned array: +/- 1 pixels\n"
     ]
    }
   ],
   "source": [
    "orig_LS = np.round(Lyot_stop_binned[L_bin/2-N:L_bin/2+N,L_bin/2-N:L_bin/2+N]).astype(int)\n",
    "dz_width_fac = float(aligntol)/1000              # dark zone is +/- this fraction of pupil diameter\n",
    "dz_width = np.ceil(2*N*dz_width_fac).astype(int) # dark zone is +/- this number of pixels in binned pupil array\n",
    "print(\"Lyot plane dark zone width in binned array: +/- {:d} pixels\".format(dz_width))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "max_shift = dz_width\n",
    "shift_range = range(-max_shift,max_shift+1,1)\n",
    "[Xshifts, Yshifts] = np.meshgrid(shift_range, shift_range)\n",
    "allowed_shifts = np.less_equal(Xshifts**2 + Yshifts**2, max_shift**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True, False],\n",
       "       [ True,  True,  True],\n",
       "       [False,  True, False]], dtype=bool)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "allowed_shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])\n",
    "len(XYshifts_allowed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fat_LS = np.ones(orig_LS.shape)\n",
    "for (xshift,yshift) in XYshifts_allowed:\n",
    "    LS_shifted = np.roll(np.roll(orig_LS, yshift, 0), xshift, 1)\n",
    "    fat_LS *= LS_shifted\n",
    "\n",
    "inv_thin_LS = np.ones(orig_LS.shape)\n",
    "for (xshift,yshift) in XYshifts_allowed:\n",
    "    inv_LS_shifted = 1-np.roll(np.roll(orig_LS, yshift, 0), xshift, 1)\n",
    "    inv_thin_LS *= inv_LS_shifted\n",
    "thin_LS = 1-inv_thin_LS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Lyot_dz = np.logical_xor(thin_LS, fat_LS)\n",
    "L = Lyot_dz.shape[0]\n",
    "\n",
    "if symm is 'half':\n",
    "    Lyot_dz_crop = Lyot_dz[:,N:]\n",
    "else:\n",
    "    Lyot_dz_crop = Lyot_dz[N:,N:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1173c9a90>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQwAAAT/CAYAAABAa78BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3V+orOtdH/DfI5scCVaF0iSVIgpCKLSk5NjYIMQLi614\n02IvAkJLS5GCQggUUtArA72wkAYhF16IWNpelEgJFWqwBkSMKE0kpdQGxB6shIQarAkROYl9erHW\n+z1PZs+atWatWfP++3wgnMmsOfvMXnuv553nO9/fPK33XgAAAAAAVVXfMPcTAAAAAACWQ2AIAAAA\nAITAEAAAAAAIgSEAAAAAEAJDAAAAACAEhgAAAABACAwBAAAAgHgx9xN4iNbaX6yqv1NVr1XVn837\nbAAAAABgdb6xqr6jqj7ee//iqQeuIjCsm7Dw3839JAAAAABg5X6kqv79qQesZST5tbmfAAAAAABs\nwGv3PWAtgaExZAAAAAB4untztrUEhgAAAADAFQgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgXsz9BAAAeD6997mfwrNprc39FAAANknDEAAAAAAIDUMAgGe0\n5Ybf3Ob+3mo4AgBbpWEIAAAAAITAEAAAAAAII8kAAIO5x1wvbctjs3P/WV36v7/lPysAYF00DAEA\nAACA0DAEADZn7ubZSGvs+cz9vb3037On/Hpzfy8AgG3RMAQAAAAAQmAIAAAAAISRZABgkeYaKzba\nyUNd+u/KU/7OG2cGAC5JwxAAAAAACA1DAOAqrtkY1JhijZ7y9/aa7UQ/XwCwfRqGAAAAAEBoGAIA\nT/JczUEtJni4a7YTH/p4P8MAsF4ahgAAAABACAwBAAAAgDCSDAAcdelRY+OJsEzn/mw+dG0wugwA\n66VhCAAAAACEhiEA7IzmIPAUD/2Z10QEgPXSMAQAAAAAQmAIAAAAAISRZADYoEuMHRv/A57imqPL\n1isAuCwNQwAAAAAgNAwBYIWe0iDUxAGW5BJNxPvWROseAJxHwxAAAAAACIEhAAAAABBGkgFgQS5x\nWEmV8Ttge06ta/etnQ9dW62dAHBDwxAAAAAACA1DAJiJg0sALuO+NfGh6+2xx1lvAdgjDUMAAAAA\nIASGAAAAAEAYSQaAZ3Tu2LHRN4DLe+jaemzNPrWOW7MB2CoNQwAAAAAgNAwB4AkcXAKwHcfW5VPr\n/H3XAOs8AGulYQgAAAAAhMAQAAAAAAgjyQBwD2PHAPt1ah2/7/rgwBQA1krDEAAAAAAIDUMAGJzb\nJtQQAdiv+64B5x6Y4poCwFJoGAIAAAAAITAEAAAAAMJIMgC7Y+wYgGs498AUh6QAsBQahgAAAABA\naBgCsGnahAAs0bHrjUNSAFgKDUMAAAAAIASGAAAAAEAYSQZgMx46fmyEC4AlusQhKa5xAFyChiEA\nAAAAEBqGAKzGuQeYVGlaALANDz0k5b5rpesiAA+hYQgAAAAAhMAQAAAAAAgjyQAskvFjADjtoWPK\n933d9ROAQxqGAAAAAEBoGAIwu4e2CTUgAOC0+66VDz0oxTUXYN80DAEAAACAEBgCAAAAAGEkGYCr\nMn4MAPN56EEpxpQB9k3DEAAAAAAIDUMAno02IQAsn9YhAIc0DAEAAACAEBgCAAAAAGEkGYCLMH4M\nANthTBlg3zQMAQAAAIDQMATgbNqEALA/WocA+6FhCAAAAACEwBAAAAAACCPJAJxk/BgAuIsxZYBt\n0jAEAAAAAELDEICjTjULNQIAgLuc2zr0ugJgeTQMAQAAAIAQGAIAAAAAYSQZAOPHAMCzOjWm7EAU\ngOXRMAQAAAAAQsMQYGdOtQlH3tkHAJ7T9Frj1IEoxx4PwPPTMAQAAAAAQmAIAAAAAISRZIANM34M\nACzdqQNR7rvPaxiA56FhCAAAAACEhiHABp1qFnonHgBYunNbh17fAFyWhiEAAAAAEBqGACunTQgA\n7MGp1qHPNwS4LA1DAAAAACAEhgAAAABAGEkGWJFT48cjIzgAwB5Mr3lOHYhy7PEAnKZhCAAAAACE\nhiHAwt3XKvROOQCwd6cORLnrPq+hAO6mYQgAAAAAhMAQAAAAAAgjyQALYvwYAOAy7htTdigKwN00\nDAEAAACA0DAEWIBTzULvdAMAXMb4uurUoShefwF7p2EIAAAAAITAEAAAAAAII8kAV2b8GABgfqcO\nRXEgCrB3GoYAAAAAQGgYAlzBqVZhlXesAQCWYHpNdupAlPFxAFulYQgAAAAAhMAQAAAAAAgjyQAX\nZvwYAGDdTh2Icnj71L8DsFYahgAAAABAaBgCXMipZqF3nAEA1m18PXfqUBSv+4At0DAEAAAAAEJg\nCAAAAACEkWSARzB+DACwX6cORXEgCrAFGoYAAAAAQGgYAjzQqVZhlXeOAQD2bHoteOpAlPFxAEum\nYQgAAAAAhMAQAAAAAAgjyQD38MHVAAA81KkDUcbbXk8CS6ZhCAAAAACEhiHA4NTBJt4FBgDgMcbX\nkdPrTVMswJJpGAIAAAAAITAEAAAAAMJIMrB7p8aQq4yGAABwOdNry2OvQcf7vAYF5qRhCAAAAACE\nhiGwWz5oGgCAuRx73Tm+Pp1ue30KzEHDEAAAAAAIgSEAAAAAEEaSgV0xhgwAwFKNr0un161Gk4E5\naBgCAAAAAKFhCGzasUbhxLu0AAAs1fRa9bBpeOwxAJemYQgAAAAAhMAQAAAAAAgjycDmGEMGAGAr\nDkeTR+N9XucCl6RhCAAAAACEhiGwGT4IGgCArTr2unZ8/Tvd9voXuAQNQwAAAAAgBIYAAAAAQBhJ\nBlbNGDIAAHs1vu6dXhc7CAW4BA1DAAAAACA0DIHVONYmHHkHFQCAvZpeCx87COXY4wBO0TAEAAAA\nAEJgCAAAAACEkWRg8U6NIhupAACANxw7CGU03ed1NHCKhiEAAAAAEBqGwCL5gGYAAHiaY6+fp9fZ\n4+ttr7OBQxqGAAAAAEAIDAEAAACAMJIMLIpRZAAAeD7Ta+vxdbeDUIBDGoYAAAAAQGgYArPTKgQA\ngOsaX287CAU4pGEIAAAAAITAEAAAAAAII8nAbIwiAwDA/ByEAhzSMAQAAAAAQsMQuCqtQgAAWCYH\noQATDUMAAAAAIDQMgavQLAQAgPXwuYawbxqGAAAAAEAIDAEAAACAMJIMXNyx8eOJ8QUAAFiP+w5C\nOfY4YP00DAEAAACA0DAELuJUq7DKO44AALB2xw5CmYz3ee0P66dhCAAAAACEwBAAAAAACCPJwJP4\nwGMAANiXY6/3x33BdNu+ANZLwxAAAAAACA1D4GxahQAAwGjcD0z7BQehwHppGAIAAAAAITAEAAAA\nAMJIMvBgRpEBAID7THsEB6HAemkYAgAAAAChYQicpFUIAAA8hoNQYL00DAEAAACAEBgCAAAAAGEk\nGTjKKDIAAHApDkKBddEwBAAAAABCwxAIrUIAAOA5OQgF1kHDEAAAAAAIgSEAAAAAEEaSAaPIAADA\n1TkIBZZLwxAAAAAACA1D2CmtQgAAYAkchALLo2EIAAAAAITAEAAAAAAII8mwM0aRAQCApXIQCiyD\nhiEAAAAAEAJDAAAAACAEhgAAAABACAwBAAAAgHDoCeyAg04AAIA1Gfcr035m3NfYz8Dz0jAEAAAA\nAOJJgWFr7V+01v5fa+1DB/f/VGvtc621P22t/Upr7bsOvv5Ka+0jrbU/aq19ubX20dbaW57yXICX\n9d5fahe21rwbBwAArMaxPcyxvQ5wOY8ODFtrf7OqfrSqPnNw/weq6sdvv/auqvpKVX28tfam4WEf\nrqofqqofrqr3VNW3VdUvPva5AAAAAACX8ajAsLX2TVX1b6vqn1bV/z348vuq6oO991/qvf/3qvqH\ndRMI/r3bf/ebq+qfVNX7e++/1nv/nar6x1X1va21dz3utwEAAAAAXMJjG4Yfqar/1Hv/xHhna+07\nq+ptVfWr03299y9V1W9V1btv7/ruujlsZXzMZ6vqD4bHAE9wWM+fKvxGkQEAgLU6tq8xmgzP4+xT\nkltr762qv1E3wd+ht1VVr6ovHNz/hduvVVW9tapevw0S73oMAAAAADCDswLD1tpfqZvPH/zbvfev\nPs9TAh7j2LtqGoUAAMAWTXudaR90OGEFPM25I8mvVtVfqqpPt9a+2lr7alV9X1W9r7X2et20BFvd\ntAhHb62qz9/e/nxVven2swzvegwAAAAAMINzA8P/UlV/vW5Gkt9x+7//WjcHoLyj9/77dRP6ff/0\nL9wGg99TVZ+8vetTVfW1g8e8vaq+vap+81G/CwAAAADgIs4aSe69f6Wq/sd4X2vtK1X1xd77797e\n9eGq+snW2u9V1WtV9cGq+sOq+tjtr/Gl1trPVdWHWmt/XFVfrqqfqarf6L3/9hN+L7BLRpEBAIC9\nOhxNHm/bF8HjnX3oyRFfl1b03n+6tfbmqvrZqvrWqvr1qvrB3vvrw8PeX1V/XlUfrapXquqXq+rH\nLvBcAAAAAIAnaGs4fry19s66GWUGSqsQAADgLvZLcK9Xe++fPvWAcz/DEAAAAADYMIEhAAAAABCX\n+AxD4EpU6wEAAE5zEAo8nYYhAAAAABAahrBwWoUAAADnG/dN075q3F/ZV8HdNAwBAAAAgBAYAgAA\nAAAhMAQAAAAAQmAIAAAAAIRDT2ChHHYCAABwGdNeatxnTbfts+BlGoYAAAAAQAgMAQAAAIAwkgwL\nYgwZAADg+Yz7q2n/Ne7D7L/ghoYhAAAAABAahrAAmoUAAADX5SAUuJuGIQAAAAAQAkMAAAAAIIwk\nw0yMIQMAAMzPQSjwMg1DAAAAACAEhgAAAABACAwBAAAAgBAYAgAAAADh0BO4MoedAAAALNO0Nxv3\nbdNt+zb2RMMQAAAAAAgNQ7gCrUIAAID1GPdr035u3NfZz7F1GoYAAAAAQGgYwjPSLAQAAFg3n2vI\nHmkYAgAAAAAhMAQAAAAAwkgyXJgxZAAAgO1xEAp7omEIAAAAAITAEAAAAAAIgSEAAAAAEAJDAAAA\nACAcegIX4rATAACAfZj2euM+cLptH8gWaBgCAAAAAKFhCE+gVQgAALBf4/5v2h+O+0T7Q9ZKwxAA\nAAAACIEhAAAAABBGkuERjCIDAAAwchAKW6JhCAAAAACEwBAAAAAACIEhAAAAABACQwAAAAAgHHoC\nD+SgEwAAAO4z7hOnfeS4n7SPZA00DAEAAACA0DCEe2gWAgAA8BjT3nHcV0637StZMg1DAAAAACAE\nhgAAAABACAwBAAAAgBAYAgAAAADh0BM4wkEnAAAAXMq4n5z2m+O+036TpdEwBAAAAABCYAgAAAAA\nhJFkGBhFBgAA4DlNe8xx/zndtv9kKTQMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEQ0/YPQed\nAAAAcG3jvnPal477U/tS5qRhCAAAAACEhiG7pVkIAADAEkx70XGfOt22T2UOGoYAAAAAQAgMAQAA\nAIAQGAIAAAAAITAEAAAAAMKhJ+yKg04AAABYqnF/Ou1fx32s/SvXomEIAAAAAITAEAAAAAAII8ns\nglFkAAAA1mTas4772em2/SzPTcMQAAAAAAiBIQAAAAAQAkMAAAAAIASGAAAAAEA49ITNctAJAAAA\nazfuY6d97rjftc/lOWgYAgAAAAAhMAQAAAAAQmAIAAAAAITAEAAAAAAIh56wOQ47AQAAYIumve24\n751u2/dySRqGAAAAAEBoGLIZh81C764AAACwReN+d9oLaxpySRqGAAAAAEAIDAEAAACAEBgCAAAA\nACEwBAAAAADCoSes2uFBJ1U+4BUAAID9mPbAh4efjF+Dc2kYAgAAAAAhMAQAAAAAQmAIAAAAAITA\nEAAAAAAIh56wSg47AQAAgDccHn4y3rZf5lwahgAAAABACAwBAAAAgBAYAgAAAAAhMAQAAAAAwqEn\nrIaDTgAAAOC0cZ887aPH/bR9NA+hYQgAAAAAhIYhi6dZCAAAAOeb9s7jvnq6bV/NKRqGAAAAAEAI\nDAEAAACAEBgCAAAAACEwBAAAAADCoScskoNOAAAA4DLG/fS03x733fbbHNIwBAAAAABCYAgAAAAA\nhMAQAAAAAAiBIQAAAAAQDj1hURx2AgAAAM9n2mOP++/ptv03Ew1DAAAAACAEhgAAAABACAwBAAAA\ngBAYAgAAAADh0BNm56ATAAAAuK5x3z3ty8f9uX35vmkYAgAAAAAhMAQAAAAAQmAIAAAAAITAEAAA\nAAAIh54wG4edAAAAwPymvfi4T59u26fvk4YhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIBw6AlX\n5aATAAAAWKZxfz7t38d9vP37fmgYAgAAAAAhMAQAAAAAQmAIAAAAAITAEAAAAAAIh55wFQ47AQAA\ngPWY9uzjfn66bT+/fRqGAAAAAEAIDAEAAACAEBgCAAAAACEwBAAAAADCoSc8GwedAAAAwLqN+/hp\nnz/u9+3zt0nDEAAAAAAIgSEAAAAAEAJDAAAAACAEhgAAAABACAwBAAAAgBAYAgAAAAAhMAQAAAAA\n4sXcT4Dt6b2/dF9rbYZnAgAAAFzKtLcf9/3Tbfv+bdEwBAAAAABCYAgAAAAAhMAQAAAAAAiBIQAA\nAAAQDj3hIhx0AgAAAPsw7venPGDMBeQB66dhCAAAAACEwBAAAAAACIEhAAAAABACQwAAAAAgHHrC\nkzjsBAAAAPZrygDGfGC6LR9YLw1DAAAAACAEhgAAAABACAwBAAAAgBAYAgAAAAAhMAQAAAAAQmAI\nAAAAAMSLuZ8A6zMelT5xVDoAAADs15gLTLnBmB/IDdZFwxAAAAAACIEhAAAAABACQwAAAAAgBIYA\nAAAAQDj0hAdz2AkAAABwnykrGHOE6bYcYR00DAEAAACA0DAEANiwYxMCW6GhAADwPDQMAQAAAIAQ\nGAIAAAAAYSQZAGBFtjxifK5zvxdGmAEAHkbDEAAAAAAIDUNOOvbOvXfnAeAyrtEW3PJ1+9zv32O+\n31v+/gHAcxuvo9N1eLweu84ul4YhAAAAABACQwAAAAAgjCQDADyjS48dG915w7nfi8f8WZz6d/xZ\nAABbpWEIAAAAAISGIQDAEzylQaihdl2P+X6f+vO978/eny8AsFYahgAAAABACAwBAAAAgDCSzFHH\nRmyM1QCwV8aO9+vUn999fy8cmAIAb5iufeP1cbrturg8GoYAAAAAQGgYAgAMzm0Tekd8v+77sz/3\nwBR/lwCApdAwBAAAAABCYAgAAAAAhJFkAGC3Hjp+bFSUxzj3wBRjygDAUmgYAgAAAAChYUh4VxuA\nLdMmZEmO/T3TOgRgD8br2HSdG693rnPLoGEIAAAAAITAEAAAAAAII8kAwOYYP2aNjCkDAEuhYQgA\nAAAAhIYh3qUGYNW0CdkyrUMAtmy6Vo3Xsem269i8NAwBAAAAgBAYAgAAAABhJBkAWKVTo8hGWNiy\nc8eU/TwAAOfSMAQAAAAAQsMQAFg8bUI47VTr0IEoAMC5NAwBAAAAgBAYAgAAAABhJBkAWJRT48cj\nI5Vw2vQzcupAlGOPBwDQMAQAAAAAQsNwp7yrDMCS3NcqdI2Cxzt1IMpd9/mZA+CaxuvOsUO7XJeu\nT8MQAAAAAAiBIQAAAAAQRpIBgNn4iAyYx31jytNtP48AsE8ahgAAAABAaBgCAFdx6mATLSaY330f\nOH/scQDANmkYAgAAAAAhMAQAAAAAwkjyzhgrAeCaTo0hV7kGwVJNP5vHfobH+/wMA3Bpx65BDuO6\nPg1DAAAAACA0DAGAi9Noh2049nOr8QEA26dhCAAAAACEwBAAAAAACCPJAMBFGEOGfRh/rqefeweh\nAMC2aBgCAAAAAKFhuAMaHwBc2rFry8h1BvZh+lk/dhDKsccBwENptM9LwxAAAAAACIEhAAAAABBG\nkgGABzs1imwsBPbr2NjYaLrPOgEA66BhCAAAAACEhiEAcK/DxpCWEHCXY+vD4YfVW0MAYNk0DAEA\nAACAEBgCAAAAAGEkGQA46tjBBcYIgceY1o7D0eTxawDAcmgYAgAAAAChYbhhmiEAPIbrB/BcDpuG\n423rDAB3cf24Pg1DAAAAACAEhgAAAABAGEkGAIwhA1c1ri8OQgGA5dEwBAAAAABCwxAAdkyzEJib\nD7IHgOXRMAQAAAAAQsMQAHZGqxBYIp9rCADLoWEIAAAAAITAEAAAAAAII8kbY8wMAAAA2CIfX3E9\nGoYAAAAAQJwVGLbW/llr7TOttT+5/d8nW2t/9+AxP9Va+1xr7U9ba7/SWvuug6+/0lr7SGvtj1pr\nX26tfbS19pZL/GYAgLv13l9qorfWvBMLLM6xtenYGgYAPI9zG4b/u6o+UFXvrKpXq+oTVfWx1tpf\nrapqrX2gqn68qn60qt5VVV+pqo+31t40/Bofrqofqqofrqr3VNW3VdUvPuH3AAAAAABcSHvqu3St\ntS9W1T/vvf98a+1zVfWveu//+vZr31xVX6iqf9R7/w+3////VNV7e+//8fYxb6+q362qv9V7/+07\n/hvvrKpPPemJ7oTPMATgLq4RwNpYtwC4i2vEk7zae//0qQc8+jMMW2vf0Fp7b1W9uao+2Vr7zqp6\nW1X96vSY3vuXquq3qurdt3d9d90ctDI+5rNV9QfDYwCAC5lG+A4/DNooMrAGx9arY+saAHBZZ5+S\n3Fr7a1X1m1X1jVX15ar6+733z7bW3l1VvW4ahaMv1E2QWFX11qp6/TZIvOsxAAAAAMBMzg4Mq+p/\nVtU7qupbquofVNW/aa2956LPCgB4EiMawNZMa9i4vk23rW8AcFlnjyT33r/We//93vvv9N5/oqo+\nU1Xvq6rPV1Wrmxbh6K23X6vbf77p9rMM73oMAAAAADCTR3+G4cGv8Urv/X/VTej3/dMXboPB76mq\nT97e9amq+trBY95eVd9eN2POAAAAAMCMzhpJbq39y6r6z3VzSMlfqKofqarvq6ofuH3Ih6vqJ1tr\nv1dVr1XVB6vqD6vqY1U3h6C01n6uqj7UWvvjuvkMxJ+pqt+464RkHsboGQCuBcAeHB6AMv7z8OsA\nbJePqnhe536G4Vuq6heq6i9X1Z9U1X+rqh/ovX+iqqr3/tOttTdX1c9W1bdW1a9X1Q/23l8ffo33\nV9WfV9VHq+qVqvrlqvqxp/wmAAAAAIDLaMfaCEvTWntn3YwzcwetEgBcC4C9se4B4FrwKK/23j99\n6gGX+AxDAAAAAGAjBIYAAAAAQJz7GYYAwMIYwwD2ygfeA8Dz0DAEAAAAAELDEABWSKsQ4A3j+jet\nj+M6aX0EgPNoGAIAAAAAITAEAAAAAMJIMgCsiFFkgNMchAIAT6dhCAAAAACEhuGKaZkAAAAAe3bq\n4CsZyeNpGAIAAAAAITAEAAAAAMJIMgAsnI+gADjfqRG1w68DAF9PwxAAAAAACIEhAAAAABACQwAA\nAAAgBIYAAAAAQDj0BAAWymEnAJcxrZ3jujrdtq4CwMs0DAEAAACA0DAEgAXRKgR4PuN6Oq2347pr\nvQWAGxqGAAAAAEAIDAEAAACAEBgCAAAAACEwBAAAAADCoScr5APxAbbH2g5wXdMaO66/023rL8A6\nHa7tDrZ6PA1DAAAAACAEhgAAAABACAwBAAAAgBAYAgAAAADh0BMAmImDTgDmN667PiQfAG5oGAIA\nAAAAITAEAAAAAEJgCAAAAACEwBAAAAAACIeeAMCVOewEYJmmtXhcp6fb1mkA9kTDEAAAAAAIgSEA\nAAAAEAJDAAAAACAEhgAAAABAOPRkRQ4/JN8HLwOsh4NOANZjXJ+n9Xtcx63fAMvmEKun0zAEAAAA\nAEJgCAAAAACEwBAAAAAACIEhAAAAABAOPQGAZ+SwE4B188H5AOyRhiEAAAAAEAJDAAAAACAEhgAA\nAABACAwBAAAAgBAYAgAAAAAhMAQAAAAA4sXcTwAAtqb3/tJ9rbUZngkAlzKu49M6P6731nkAtkTD\nEAAAAAAIgSEAAAAAEAJDAAAAACAEhgAAAABACAwBAAAAgBAYAgAAAAAhMAQAAAAA4sXcT4DTeu8v\n3ddam+GZAHAfazbAPkxr+7juT7et+wDLMa7J0zo9rt3W7LtpGAIAAAAAITAEAAAAAEJgCAAAAACE\nwBAAAAAACIEhAAAAABACQwAAAAAgXsz9BABgzXrvL93XWpvhmQBwbeN6P10PxuuC6wEAa6VhCAAA\nAACEwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAA\nAACEwBAAAAAACIEhAAAAABACQwAAAAAgXsz9BABgjXrvL93XWpvhmQCwBNM1YLw+TLddHwBYGw1D\nAAAAACAEhgAAAABACAwBAAAAgBAYAgAAAAAhMAQAAAAAQmAIAAAAAITAEAAAAAAIgSEAAAAAEAJD\nAAAAACAEhgAAAABACAwBAAAAgHgx9xPguN77S/e11mZ4JgBMrM0A3Ge8LkzXjfH64boBMI9p/R3X\n5Om2tfllGoYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAAAAAACIEhAAAAABAC\nQwAAAAAPU0E7AAAgAElEQVQgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAA\nAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAA\nAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAA\nAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAA\nAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAA\nAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAA\nAAAAiBdzPwGOa61VVVXvPfdNt6evAXBd4/o7rcnjOm19BmC8LkxcHwDmZ30+j4YhAAAAABACQwAA\nAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAEAAAAAEJgCAAAAACEwBAAAAAACIEhAAAAABACQwAA\nAAAgBIYAAAAAQLyY+wkAwBq11qqqqvee+6bb09cA2I/xejBxPQBgrTQMAQAAAIAQGAIAAAAAITAE\nAAAAAEJgCAAAAACEwBAAAAAACIEhAAAAABACQwAAAAAgBIYAAAAAQAgMAQAAAIAQGAIAAAAAITAE\nAAAAAEJgCAAAAADEi7mfAACsWWstt3vvX/fPw68DsC3jej+x7gOwBRqGAAAAAEAIDAEAAACAEBgC\nAAAAACEwBAAAAABCYAgAAAAAhMAQAAAAAIgXcz8BTmut5Xbv/ev+OX4NgPlN6/K0To+3rdkA2zGu\n8xPrPMDyWK8fT8MQAAAAAAiBIQAAAAAQAkMAAAAAIASGAAAAAEAIDAEAAACAEBgCAAAAACEwBOD/\nt3d/odbl5X3Av08yzhhtrRfWmUqRGCTThhYbx2gkOE1QkpqASTE0hkCDXoiNglgKNtQSUOiFQTOY\nKHghNdIrUURbqEZNUvw/dIz2QiNIFDVmJhpFhzEy/vn1Yu/1uDyz333+7XP2v88Hhlln7fWed70z\n79p7nWd9n98DAAAA7aZtnwAAHJqq6u0xxg/9++TrAOyH+fv4xPs5AIdKwhAAAAAAaAqGAAAAAEBT\nMAQAAAAAmoIhAAAAANAUDAEAAACApmAIAAAAALSbtn0CAHDIqipJMsbofdP29BoAu2v+/j3x/g3A\noZMwBAAAAACagiEAAAAA0BQMAQAAAICmYAgAAAAANENP9sjJhfPnCzBbeBlgt83fp72PA+w2g04A\n9pv38cuTMAQAAAAAmoIhAAAAANAUDAEAAACApmAIAAAAADRDTwDgmp0cYjXfthgzwPZYJB8AFiQM\nAQAAAICmYAgAAAAANAVDAAAAAKApGAIAAAAAzdATANiS+UL600L78wX3LbQPcPUMOgGAh5IwBAAA\nAACagiEAAAAA0BQMAQAAAICmYAgAAAAANENP9tC0CPN8geZp2wLNAPvJezvA9TLsBODwnHxv975+\ncRKGAAAAAEBTMAQAAAAAmoIhAAAAANAUDAEAAACAZugJAOyQ+cLM06LN88WbLdwMcHEGnQDA2UgY\nAgAAAABNwhAAdtSUepknYqZtiRiAs5MsBIDzkTAEAAAAAJqCIQAAAADQFAwBAAAAgKZgCAAAAAA0\nQ08AYMfNF+afFu6fL+Bv4X6AhzLoBAAuTsIQAAAAAGgKhgAAAABA05K8x7SoAQAAAMfMEhRXQ8IQ\nAAAAAGgShgCwR6anpfMnqdO2J6kAkiYAsAkShgAAAABAUzAEAAAAAJqWZADYQwZfAfyANmQA2CwJ\nQwAAAACgSRgCwJ4zCAU4VpKFAHA1JAwBAAAAgHaugmFV/W5V3V1V36yq+6rqHVX1kyuOe2VVfbmq\nvlVV762qJ554/Zaqen1VfbWq7q+qt1XVYy/7hwEAAAAALue8CcNnJPnDJE9L8qwkD0vyJ1X1Y9MB\nVfXyJC9J8sIkT03yQJL3VNXNs+9zV5JfSfLcJHcmeVySt1/wzwAAAAAAbEitWvfjzL+46jFJ/jbJ\nnWOMDy73fTnJ748x/mD59aOS3Jfkt8cYb11+/ZUkzxtjvGN5zO1JPp3kZ8cYd6/4fZ6c5J4Ln+gR\nsH4LAD4LgGPjfQ8AnwUXcscY4+PrDrjsGoaPTjKSfC1JquoJSW5L8v7pgDHGN5N8LMnTl7ueksWw\nlfkxn0nyhdkxnFNVPeSCGGOsvHAAOEzTZ8H882D6LPB5AByKVe9rq97/ADhsq+5xfRZszoULhrX4\nP3BXkg+OMT613H1bFgXE+04cft/ytSS5NcmDy0LijY4BAAAAALbgpkv82jck+akkP7ehcwEAAAAA\ntuxCCcOq+qMkv5zk58cYfzN76d4klUWKcO7W5WvTMTcv1zK80TEAwCVYqgI4NFrPAOD6nLtguCwW\n/mqSXxhjfGH+2hjjc1kU/Z45O/5RWUxV/vBy1z1JvnvimNuTPD7JR857PgAAAADA5pyrJbmq3pDk\nN5M8J8kDVTUlCb8xxvj2cvuuJK+oqs8m+XySVyX5UpJ3JoshKFX1piSvraqvJ7k/yeuSfGjVhGQA\n4OJODkCZ//vk6wC7xuRLANiO865h+KIshpr8+Yn9z0/yliQZY7y6qh6R5I1ZTFH+QJJnjzEenB3/\nsiTfS/K2JLckeXeSF5/35AEAAACAzap9WMuoqp6cRSszp/AUFoAb8RkB7BvvWwDciM+IS7ljjPHx\ndQdcZkoyALBHphuo+c3VtO3mCtglfggEgO260JRkAAAAAOAwSRgeGIvbAwAAAIdIAv36SBgCAAAA\nAE3BEAAAAABoWpIB4MhYvgLYRdrMAGB3SBgCAAAAAE3CEACO2JTemSd7pm3JHuA6SBYCwO6RMAQA\nAAAAmoQhAGBdQ+BaSRUCwG6TMAQAAAAAmoIhAAAAANC0JB8wC9kDcBE+P4CrohUZgIvw+XH9JAwB\nAAAAgCZhCACstG4Qiie6wHmcTIZ4DwGA3SZhCAAAAAA0BUMAAAAAoGlJBgBOdXIQioWngRtZ9f4w\n8T4BAPtBwhAAAAAAaBKGAMCZnUwazs33SRHBcVmXKky8JwDAvpEwBAAAAACagiEAAAAA0LQkH4F5\nC8iqxeq1iABwXqs+O+afLdO2zxg4bAYgAXBVfMZsl4QhAAAAANAkDAGAjZBoh+Mg8QEAh0/CEAAA\nAABoCoYAAAAAQNOSDABs3NSeuGoQyqrjgN206rqdcw0DwGGSMAQAAAAAmoThkVmX+PCEGIBNWzUI\nZc5nEOymdclC1ysAV0lXym6QMAQAAAAAmoIhAAAAANC0JAMA12JVK8nUcjJvPdFyAtuhBQwAmEgY\nAgAAAABNwhAA2Jp1w7hWHQdsxrqhJolrDgCOnYQhAAAAANAUDAEAAACApiX5SM3bTCw4D8C2rfpc\nmtOmDBd3WvvxxDUFwLa419s9EoYAAAAAQJMwBAB2yqqnyetSh54+w2rrkoWuGwBgHQlDAAAAAKAp\nGAIAAAAATUsyALDz1rUpWyQbtB8DAJslYQgAAAAANAlDAGAvTampdQNRVh0P+25dmnDO33kA4KIk\nDAEAAACApmAIAAAAADQtyaxs6Zq2tbIAsOvWDUQ5bZ/POXad9mMADpn7s90lYQgAAAAANAlDAODg\nSB2yj6QJAYBdIWEIAAAAADQFQwAAAACgaUmmzdtbppaYeWuM9hcA9pk2ZXaJ9mMAjpV7rf0gYQgA\nAAAANAlDAOBoXSZ1uO57QHL2FOHE3yUAYFdIGAIAAAAATcEQAAAAAGhakgEAZs7apnyW1270/TgM\n5205nvP3AgDYZRKGAAAAAECTMGSl6an3/Mn5tO2JOADHZt1n32kpMwNT9psUIQBsxqrPVJ+Vu0vC\nEAAAAABoCoYAAAAAQNOSDABwCae10lxmYMpFfj9u7DLtxav4fwEAHCoJQwAAAACgSRgCAFyhywxM\n2cSvOeQU3KYTg6sc8n8/AIAbkTAEAAAAAJqCIQAAAADQtCSz1rwNZ2r7mbf/aNMBgIu7yOfoedtw\nr6Ntd1+4bwGA67XqPsTn8X6QMAQAAAAAmoQhAMAeOe9T+UNOGEooAABcDQlDAAAAAKApGAIAAAAA\nTUsyAMAB07YLAMB5SRgCAAAAAE3CkDObEgrzxdOnbekFAAAAIFk9dE3dYL9IGAIAAAAATcEQAAAA\nAGgKhgAAAABAUzAEAAAAAJqhJ5zbfKHSaSHT+YKmFjIFAACA42LQyWGRMAQAAAAAmoIhAAAAANAU\nDAEAAACApmAIAAAAADQFQwAAAACgKRgCAAAAAO2mbZ8A+20akT4fnz5tG58OAAAAh21eD5ioB+w/\nCUMAAAAAoCkYAgAAAABNwRAAAAAAaAqGAAAAAEAz9ISNmC9oOi14Ol/41IKnAAAAcBgMOjl8EoYA\nAAAAQFMwBAAAAACagiEAAAAA0BQMAQAAAIBm6AkbNy10Ol8Eddq2CCoAAADsJ8NOjoeEIQAAAADQ\nFAwBAAAAgKZgCAAAAAA0BUMAAAAAoCkYAgAAAABNwRAAAAAAaAqGAAAAAEC7adsnwOGqqt4eY/zQ\nv0++DgAAAOye+c/xEz/PHz4JQwAAAACgKRgCAAAAAE3BEAAAAABoCoYAAAAAQDP0hGsxLYg6Xyx1\n2rZYKgAAAOwWw06Om4QhAAAAANAUDAEAAACApmAIAAAAADQFQwAAAACgGXrCtZovkDotoDpfSNUC\nqgAAALAdBp0wkTAEAAAAAJqCIQAAAADQFAwBAAAAgKZgCAAAAAA0Q0/Ymmnh1PmiqtO2RVUBAADg\nehh2wkkShgAAAABAUzAEAAAAAJqCIQAAAADQFAwBAAAAgGboCVs3X0h1Wmh1vuCqhVYBAABgsww6\nYR0JQwAAAACgKRgCAAAAAE3BEAAAAABoCoYAAAAAQDP0hJ0yLbA6X3x12rb4KgAAAFyOYSechYQh\nAAAAANAUDAEAAACApmAIAAAAADQFQwAAAACgGXrCTpovuDotyDpfmNWCrAAAAHA2Bp1wXhKGAAAA\nAEBTMAQAAAAAmoIhAAAAANAUDAEAAACAZugJO29aiHW+SOu0bZFWAAAAWM2wEy5KwhAAAAAAaBKG\n7I35U5DpKcn8aYmnJAAAABw7qUI2QcIQAAAAAGgKhgAAAABAUzAEAAAAAJqCIQAAAADQDD1hL00L\nts4Xc522LeYKAADAsTHshE2SMAQAAAAAmoIhAAAAANAUDAEAAACApmAIAAAAADRDT9hr8wVcpwVe\nDT8BAADgWJwcduJnYTZBwhAAAAAAaAqGAAAAAEBTMAQAAAAAmoIhAAAAANAMPeFgTAu7nhx+Mn8N\nAAAA9t3JQSeJn3vZLAlDAAAAAKBJGHJwTiYN59ueuAAAALCvJAu5LhKGAAAAAEBTMAQAAAAAmoIh\nAAAAANAUDAEAAACAZugJB2u+8Ou0MOx8gVgLwwIAALDrDDphGyQMAQAAAICmYAgAAAAANAVDAAAA\nAKApGAIAAAAAzdATjsK0IOx8sdhp22KxAAAA7BrDTtgmCUMAAAAAoCkYAgAAAABNSzJHZR7fnuLd\n85i3eDcAAADbog2ZXSFhCAAAAAA0BUMAAAAAoCkYAgAAAABNwRAAAAAAaIaecLSmhWPni8pO2xaV\nBQAA4LoYdsKukTAEAAAAAJqEIUdv/tRmeqozf7rjqQ4AAACbJlXILpMwBAAAAACagiEAAAAA0BQM\nAQAAAICmYAgAAAAANENPYGZaYHa++Oy0bfFZAAAALsuwE/aBhCEAAAAA0BQMAQAAAICmJRlWmMfB\np7j4PDYuLg4AAMBZaUNm30gYAgAAAABNwRAAAAAAaAqGAAAAAEBTMAQAAAAAmqEncIppIdr5IrXT\ntkVqAQAAuBHDTthXEoYAAAAAQJMwhDOaPwWanhLNnxZ5SgQAAIBUIYdAwhAAAAAAaAqGAAAAAEBT\nMAQAAAAAmoIhAAAAANAMPYELmBasnS9mO21bzBYAAOD4GHbCIZEwBAAAAACagiEAAAAA0LQkwyXM\n4+VT/HweQxc/BwAAOFzakDlUEoYAAAAAQJMwhA0xCAUAAOA4SBZy6CQMAQAAAICmYAgAAAAANAVD\nAAAAAKApGAIAAAAAzdAT2LD5QrfTQrjzBXEthAsAALB/DDrhmEgYAgAAAABNwRAAAAAAaFqS4QpN\n8fR5dH3aFl0HAADYfVqROUYShgAAAABAkzCEa2AQCgAAwP6QKuTYSRgCAAAAAE3BEAAAAABoWpLh\nmhmEAgAAsJu0IsOChCEAAAAA0BQMAQAAAICmYAgAAAAANGsYwpbM18GY1smYr5dhnQwAAICrZ91C\neCgJQwAAAACgKRgCAAAAAE1LMuyAKe4+j8JP26LwAAAAm6cVGW5MwhAAAAAAaBKGsEMMQgEAALg6\nUoVwNhKGAAAAAEBTMAQAAAAAmpZk2FEGoQAAAGyGVmQ4HwlDAAAAAKApGAIAAAAATcEQAAAAAGgK\nhgAAAABAM/QEdtx8Id5pod75gr0W6gUAAHgog07g4iQMAQAAAIAmYQh7ZHoaNn9SNm17UgYAACBZ\nCJsgYQgAAAAAtHMXDKvqGVX1rqr666r6flU9Z8Uxr6yqL1fVt6rqvVX1xBOv31JVr6+qr1bV/VX1\ntqp67GX+IAAAAADA5V0kYfjIJJ9I8jtJHpLzraqXJ3lJkhcmeWqSB5K8p6punh12V5JfSfLcJHcm\neVySt1/gXOAoVVX/Mxlj9D8AAADHZNXPQ6t+bgLO5txrGI4x3p3k3UlSq6+6lyZ51Rjjfy2P+fdJ\n7kvya0neWlWPSvKCJM8bY/yf5THPT/LpqnrqGOPuC/1JAAAAAIBL2+gahlX1hCS3JXn/tG+M8c0k\nH0vy9OWup2RRqJwf85kkX5gdA5zRqidmkoYAAMAxWPWzj1QhXN6mh57clkWb8n0n9t+3fC1Jbk3y\n4LKQeKNjAAAAAIAtMCUZAAAAAGibLhjem6SySBHO3bp8bTrm5uVahjc6BjindYNQAAAADokBJ3C1\nNlowHGN8Loui3zOnfcvC4NOSfHi5654k3z1xzO1JHp/kI5s8HwAAAADgfM49JbmqHpnkiVkkCZPk\nJ6rqSUm+Nsb4YpK7kryiqj6b5PNJXpXkS0nemSyGoFTVm5K8tqq+nuT+JK9L8iETkmEzpqdq0xO3\nk0/eAAAA9s2q7ik/38DVOHfBMIspx3+WxXCTkeQ1y/1/nOQFY4xXV9UjkrwxyaOTfCDJs8cYD86+\nx8uSfC/J25LckuTdSV58oT8BAAAAALAxtQ/rm1XVk7NoZQbOwRM4AADgUPj5BjbmjjHGx9cdcJGE\nIbAnTrYmz7d9sAIAAPtAoRCu36anJAMAAAAAe0zCEI7A/OmbQSgAAMCukyqE7ZIwBAAAAACagiEA\nAAAA0BQMAQAAAICmYAgAAAAANENP4MhMCwXPFxGeti0iDAAAbJNhJ7AbJAwBAAAAgKZgCAAAAAA0\nLclwpOax/in2P4//i/0DAADXQRsy7B4JQwAAAACgSRgCBqEAAADXTrIQdpeEIQAAAADQFAwBAAAA\ngKYlGWgGoQAAAFdJGzLsBwlDAAAAAKBJGAIrGYQCAABsimQh7BcJQwAAAACgKRgCAAAAAE1LMrCW\nQSgAAMBFaEOG/SVhCAAAAAA0CUPgzAxCAQAATiNZCPtPwhAAAAAAaAqGAAAAAEDTkgycm0EoAADA\nnDZkOCwShgAAAABAkzAELmXdIJRVxwEAAPtt1f3+nHt/2H8ShgAAAABAUzAEAAAAAJqWZGAjVg1C\nmZv2aU8AAID9tK4V2X0+HBYJQwAAAACgSRgCG7fq6eL0NHL+VNJTSAAA2G0GGsJxkjAEAAAAAJqC\nIQAAAADQtCQD12JqW5i3NBiEAgAAu0krMhw3CUMAAAAAoEkYAtdq/lTSIBQAANgdUoXARMIQAAAA\nAGgShsDWWNcQAAC2T7IQOEnCEAAAAABoCoYAAAAAQNOSDGydQSgAAHC9tCED60gYAgAAAABNwhDY\nKQahAADA1ZEsBM5CwhAAAAAAaAqGAAAAAEDTkgzspNMGoaw6DgAA+IFV988T99HAOhKGAAAAAECT\nMAR23qpBKJP5Pk9JAQA4dutShYl7ZuBsJAwBAAAAgKZgCAAAAAA0LcnA3ljVPjFvuZi2tVkAAHBs\nDAgENknCEAAAAABoEobAXps/NZ2eqhqEAgDAMZAqBK6KhCEAAAAA0BQMAQAAAICmJRk4GFP7xapB\nKCePAQCAfbKq/XjOfS6wSRKGAAAAAECTMAQOzqpBKKu+9hQWAIBdty5Z6H4WuCoShgAAAABAUzAE\nAAAAAJqWZOCgnWzTWDUQRSsHAAC7xvA+YJskDAEAAACAJmEIHJVVA1EMQgEAYBesGnDi/hTYBglD\nAAAAAKApGAIAAAAATUsycLSm9o5Vg1BWHQcAAJuy6r5zzj0osE0ShgAAAABAkzAEjt6qQShz0z5P\neQEAuKx1yUL3m8CukDAEAAAAAJqCIQAAAADQtCQDzKxqA5naRubtI9pFAAA4K4P1gH0jYQgAAAAA\nNAlDgFNMT3/nT4Y9JQYAYJV1Q00S94zAfpAwBAAAAACagiEAAAAA0LQkA5zRvH1kVavJtE+bCQDA\n8VnXiuz+ENg3EoYAAAAAQJMwBLiAVU+Jp6fKBqIAABw2aULg0EkYAgAAAABNwRAAAAAAaFqSATZk\naj9ZNxBlfhwAAPtjXRty4h4POCwShgAAAABAkzAE2LB1A1FObq/7NQAAbIc0IXDsJAwBAAAAgKZg\nCAAAAAA0LckA12DetrJuKIr2FgCA7VnXiuw+DTgmEoYAAAAAQJMwBLhm64aiGIgCAHA9pAkBbkzC\nEAAAAABoCoYAAAAAQNOSDLADpraXdQNR5scBAHB+69qQE/daABMJQwAAAACgSRgC7JB1A1FObq/7\nNQAAx06aEODiJAwBAAAAgKZgCAAAAAA0LckAO27eLnPaUJRVvwYA4NCd1n48cY8EcDYShgAAAABA\nkzAE2COnDUU5uc9TdADgkK1LFroPArg4CUMAAAAAoCkYAgAAAABNSzLAnlvXpmwgCgBwKLQfA1wf\nCUMAAAAAoEkYAhyg6Sn7uoEoq44HANgF69KEc+5hAK6GhCEAAAAA0CQMAQ7YuvUNT9vniT0AcB2k\nCQF2j4QhAAAAANAUDAEAAACApiUZ4Mict01Z+w8AcBXWtSK7/wDYLglDAAAAAKBJGAKwNnVoIAoA\ncFnShAD7RcIQAAAAAGgKhgAAAABA05IMwEpTe9C6gSirjgcAjtu69uM59w4Au0vCEAAAAABoEoYA\nrLVuIMpp+yQHAOCwSRMCHCYJQwAAAACgKRgCAAAAAE1LMgDnpk0ZAI6P9mOA4yFhCAAAAAA0CUMA\nNkLqEAAOhzQhwHGTMAQAAAAAmoIhAAAAANC0JANwZbQpA8Du034MwEkShgAAAABAkzAE4FpJHQLA\n9kgTAnAWEoYAAAAAQFMwBAAAAACalmQAtu4ybcqnfR8AOCZnbTme8/kJwEkShgAAAABAkzAEYCed\nNXV42utSEwAcKmlCAK6KhCEAAAAA0BQMAQAAAICmJRmAvXFaG9VZB6VoxwJg35y1/dhnHACbIGEI\nAAAAADQJQwAOxlkHpaxLaUhmALAt5x1i4jMLgKsiYQgAAAAANAVDAAAAAKBpSQbgoJ21TXnda1q+\nANg07ccA7DIJQwAAAACgSRgCcHTWpTQMSQFgU6QIAdhXEoYAAAAAQFMwBAAAAACalmQAmNnEkJTT\nvh8Ah+G8LcdzPh8A2GUShgAAAABAkzAEgFOcd0jKWV+XLgHYfVKEABwjCUMAAAAAoCkYAgAAAABN\nSzIAXMJp7WbnHZiifQ1ge87bfuw9G4BDJWEIAAAAADQJQwC4QucdmHLWdItUC8DZGVwCAOcjYQgA\nAAAANAVDAAAAAKBpSQaALVnV5nbWtrnTjtNCBxyTy7Qcz3nvBIAFCUMAAAAAoEkYAsAOOWu65bQ0\nzbrXJWiAfWRwCQBcHwlDAAAAAKApGAIAAAAATUsyAOyh09rr1rXunbWtTwsfcB02MbDE+xUAbJaE\nIQAAAADQJAwB4ACtS9ucNc0jiQhcxiaSg3PeawDg+kgYAgAAAABNwRAAAAAAaFqSAeDInLWtT+sy\nsIpWYwA4fBKGAAAAAECTMAQAVtpWEvG8vz9wOZtODE5cwwCwvyQMAQAAAICmYAgAAAAANC3JAMCl\nbLp1+aLHz2mF5NhcVVvxKq4vADh8EoYAAAAAQJMwBACuxXlTSZdJTEknso+uMyU45+88AHCShCEA\nAAAA0BQMAQAAAICmJRkA2EmXaZPcVjvzKto9D9e2WohX8fcMANgkCUMAAAAAoEkYAgAHZ1vpxOv4\nfud1yMmzbf+33bRD/n8FAOwXCUMAAAAAoEkYAgDMbDrlte0U3LZ//0MmEQgAHCoJQwAAAACgKRgC\nAAAAAE1LMgDAFdp22+ohtyRv+78tAMChkjAEAAAAAJqEIQDAAZPCAwDgvCQMAQAAAICmYAgAAAAA\nNAVDAAAAAKApGAIAAAAATcEQAAAAAGgKhgAAAABAUzAEAAAAAJqCIQAAAADQFAwBAAAAgKZgCAAA\nAAA0BUMAAAAAoCkYAgAAAABNwRAAAAAAaAqGAAAAAEBTMAQAAAAAmoIhAAAAANAUDAEAAACApmAI\nAAAAADQFQwAAAACgKRgCAAAAAE3BEAAAAABoCoYAAAAAQFMwBAAAAACagiEAAAAA0BQMAQAAAICm\nYAgAAAAANAVDAAAAAKApGAIAAAAATcEQAAAAAGgKhgAAAABAUzAEAAAAAJqCIQAAAADQFAwBAAAA\ngKZgCAAAAAA0BUMAAAAAoCkYAgAAAABNwRAAAAAAaAqGAAAAAEBTMAQAAAAAmoIhAAAAANAUDAEA\nAD+/YeAAAAaOSURBVACApmAIAAAAADQFQwAAAACgKRgCAAAAAE3BEAAAAABoCoYAAAAAQFMwBAAA\nAACagiEAAAAA0BQMAQAAAICmYAgAAAAANAVDAAAAAKApGAIAAAAATcEQAAAAAGhbLRhW1Yur6nNV\n9fdV9dGq+pltng8AAAAAHLutFQyr6jeSvCbJ7yX56SSfTPKeqnrMts4JAAAAAI7dNhOGL0vyxjHG\nW8YYf5nkRUm+leQFWzwnAAAAADhqWykYVtXDktyR5P3TvjHGSPK+JE/fxjkBAAAAANtLGD4myY8m\nue/E/vuS3Hb9pwMAAAAAJPszJfnh2z4BAAAAADgAp9bZtlUw/GqS7yW59cT+W5Pcu+L4H7/qEwIA\nAACAI/Djpx1w0zWcxEOMMb5TVfckeWaSdyVJVdXy69et+CXvSfJbST6f5NvXdJoAAAAAcCgenkWx\n8D2nHViLWSPXr6r+XZI3ZzEd+e4spib/epJ/Nsb4ylZOCgAAAACO3FYShkkyxnhrVT0mySuzaEX+\nRJJfUiwEAAAAgO3ZWsIQAAAAANg9+zIlGQAAAAC4BgqGAAAAAEBTMAQAAAAA2t4UDKvqxVX1uar6\n+6r6aFX9zLbPCVivqp5RVe+qqr+uqu9X1XNWHPPKqvpyVX2rqt5bVU888fotVfX6qvpqVd1fVW+r\nqsde358COKmqfreq7q6qb1bVfVX1jqr6yRXHub5hz1TVi6rqk1X1jeU/H66qf3PiGNc2HICq+s/L\ne/TXntjvGgf2o2BYVb+R5DVJfi/JTyf5ZJL3LKcsA7vrkVlMQP+dJA+ZsFRVL0/ykiQvTPLUJA9k\ncW3fPDvsriS/kuS5Se5M8rgkb7/a0wZO8Ywkf5jkaUmeleRhSf6kqn5sOsD1DXvri0lenuTJSe5I\n8qdJ3llV/zxxbcOhWAZwXpjFz9bz/a5xIMmeTEmuqo8m+dgY46XLryuLm5nXjTFevdWTA86kqr6f\n5NfGGO+a7ftykt8fY/zB8utHJbkvyW+PMd66/PorSZ43xnjH8pjbk3w6yc+OMe6+7j8H8FDLB3h/\nm+TOMcYHl/tc33AgqurvkvynMcZ/d23D/quqf5DkniT/Icl/TfIXY4z/uHzNNQ4k2YOEYVU9LIun\nm++f9o1FlfN9SZ6+rfMCLqeqnpDktvzwtf3NJB/LD67tpyS56cQxn0nyhbj+YZc8OosU8dcS1zcc\niqr6kap6XpJHJPmwaxsOxuuT/M8xxp/Od7rGgbmbtn0CZ/CYJD+axVONufuS3H79pwNsyG1ZFBhW\nXdu3LbdvTfLg8kblRscAW7RM/d+V5INjjE8td7u+YY9V1b9I8pEkD09yf5J/O8b4TFU9Pa5t2GvL\nhwD/KovC30k+v4G2DwVDAGB3vSHJTyX5uW2fCLAxf5nkSUn+UZJfT/KWqrpzu6cEXFZV/dMsHvI9\na4zxnW2fD7Dbdr4lOclXk3wviycZc7cmuff6TwfYkHuTVNZf2/cmuXm5VsqNjgG2pKr+KMkvJ/n5\nMcbfzF5yfcMeG2N8d4zxV2OMvxhj/JcshiK8NK5t2Hd3JPnHST5eVd+pqu8k+ddJXlpVD2aREnSN\nA0n2oGC4fPJxT5JnTvuW7U/PTPLhbZ0XcDljjM9lcVMxv7YflcXU1enavifJd08cc3uSx2fRKgVs\nybJY+KtJfmGM8YX5a65vODg/kuQW1zbsvfcl+ZdZtCQ/afnP/03yP5I8aYzxV3GNA0v70pL82iRv\nrqp7ktyd5GVZLL785m2eFLBeVT0yyROzeFKZJD9RVU9K8rUxxhezaIl4RVV9Nsnnk7wqyZeSvDNZ\nLLJcVW9K8tqq+noW6yi9LsmHTGCD7amqNyT5zSTPSfJAVU1JhG+MMb693HZ9wx6qqv+W5H9nMcDg\nHyb5rSwSSL+4PMS1DXtqjPFAkk/N91XVA0n+bozx6eUu1ziQZE8Khsvx7Y9J8sosos6fSPJLY4yv\nbPfMgFM8JcmfZbF48kjymuX+P07ygjHGq6vqEUnemMWU1Q8kefYY48HZ93hZFssSvC3JLUneneTF\n13P6wA28KItr+s9P7H9+krckiesb9tZjs/ic/idJvpHk/yX5xWmaqmsbDs74oS9c48BSjTFOPwoA\nAAAAOAo7v4YhAAAAAHB9FAwBAAAAgKZgCAAAAAA0BUMAAAAAoCkYAgAAAABNwRAAAAAAaAqGAAAA\nAEBTMAQAAAAAmoIhAAAAANAUDAEAAACApmAIAAAAALT/DyG3wC+08B2gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110a61650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(16,16))\n",
    "plt.imshow(Lyot_dz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alignment tolerance check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The LDZ accomomdates a translation +/-0.2% of D=500 pixels = +/-1.00 pixels, up to +/-1 whole pixels\n",
      "Testing an (x,y) translation of (1, 0) pixels. Within the design tolerance? True\n",
      "LDZ encompasses the LS transmission error region and does not overlap with tel pupil obscurations? True\n",
      "Total unconstrained \"leak\" area after translation = 0 pixels\n"
     ]
    }
   ],
   "source": [
    "max_shift_tol = 2*N*dz_width_fac\n",
    "max_shift_tol_int = int(np.round(max_shift_tol))\n",
    "test_shift = (max_shift_tol_int,0)\n",
    "print(\"The LDZ accomomdates a translation +/-{0:.1f}% of D={1:d} pixels = +/-{2:.2f} pixels, up to +/-{3:d} whole pixels\".format(\n",
    "      float(aligntol)/10, 2*N, max_shift_tol, max_shift_tol_int))\n",
    "print(\"Testing an (x,y) translation of {0:} pixels. Within the design tolerance? {1:}\".format(\n",
    "      test_shift, test_shift[0]**2 + test_shift[1]**2 <= max_shift_tol))\n",
    "shift_LS = np.roll(np.roll(orig_LS, test_shift[0], axis=1), test_shift[1], axis=0)\n",
    "LS_err_mask = np.ceil(np.abs(shift_LS - orig_LS)).astype(bool)\n",
    "\n",
    "LDZ_valid = ~np.any(np.logical_and(LS_err_mask, ~Lyot_dz))\n",
    "print(\"LDZ encompasses the LS transmission error region and does not overlap with tel pupil obscurations? {0:}\".format(LDZ_valid))\n",
    "\n",
    "print(\"Total unconstrained \\\"leak\\\" area after translation = {0:d} pixels\".format(\n",
    "      int(np.sum(np.logical_and(LS_err_mask, ~Lyot_dz)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write the Lyot dark zone file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if aligntol > 0 and LDZ_valid and write_LDZ:\n",
    "    if symm is 'half':\n",
    "        if shape in ['ann', 'hex']:\n",
    "            LDZ_dat_fname_tail = \"LDZ_half_{0:s}{1:02d}D{2:02d}_clear_Tol{3:02d}_N{4:04d}.dat\".format(\n",
    "                                  shape, int(round(100*inD)), int(round(100*outD)), aligntol, N)\n",
    "        else:\n",
    "            LDZ_dat_fname_tail = \"LDZ_half_{0:s}P{1:02d}D{2:02d}_clear_Tol{3:02d}_N{4:04d}.dat\".format(\n",
    "                                  shape, int(round(100*inD)), int(round(100*outD)), aligntol, N)            \n",
    "    else:\n",
    "        if shape in ['ann', 'hex']:   \n",
    "            LDZ_dat_fname_tail = \"LDZ_quart_{0:s}{1:02d}D{2:02d}_clear_Tol{3:02d}_N{4:04d}.dat\".format(\n",
    "                                  shape, int(round(100*inD)), int(round(100*outD)), aligntol, N)\n",
    "        else:\n",
    "            LDZ_dat_fname_tail = \"LDZ_quart_{0:s}P{1:02d}D{2:02d}_clear_Tol{3:02d}_N{4:04d}.dat\".format(\n",
    "                                  shape, int(round(100*inD)), int(round(100*outD)), aligntol, N)            \n",
    "\n",
    "    LDZ_dat_fname = os.path.join(LS_dir, LDZ_dat_fname_tail)\n",
    "    \n",
    "    if not os.path.exists(LDZ_dat_fname) or overwrite==True:\n",
    "        np.savetxt(LDZ_dat_fname, Lyot_dz_crop, fmt='%d', delimiter=\" \")\n",
    "        print(\"Wrote binned, cropped LDZ array to {0:s}\".format(LDZ_dat_fname))\n",
    "    else:\n",
    "        print(\"LDZ array {0:s} already exists, will not overwrite\".format(LDZ_dat_fname))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:py27]",
   "language": "python",
   "name": "conda-env-py27-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}