make_Lyotstop_obs1.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",
    "%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": [
    "telap_key_map = {'hex1':'hex1', 'hex2':'hex2', 'hex3':'hex3', 'hex4':'hex4',\n",
    "                 'key24':'keystone24', 'pie08':'piewedge8', 'pie12':'piewedge12'}\n",
    "secobs_key_map = {'Cross':'cross', 'X':'x'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set basic parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 336,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "overwrite = False\n",
    "N = 250 # pupil array quadrant width after binning\n",
    "iD = 30\n",
    "oD = 76\n",
    "aligntol = 2 # units of thousandths of pupil diameter\n",
    "pad = 0 # units of thousandths of pupil diameter\n",
    "L = 1000 # starting array width\n",
    "D = 1000 # starting pupil diameter\n",
    "symm = 'quart' # set to either 'quart' or 'half'\n",
    "shape = 'ann' # set to 'ann' or 'hex'\n",
    "edge = 'round' # set to 'gray' or 'round'\n",
    "secobs_key = \"X\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 337,
   "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.abspath(\"../InputMasks/LS\")\n",
    "LS_dir = os.path.abspath(\"/astro/opticslab1/SCDA/Apertures/InputMasks/LS\")\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": [
    "## Load secondary obscuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 338,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#telap_dir = os.path.abspath('../Apertures/JPL/offset_masks')\n",
    "telap_dir = os.path.abspath('/astro/opticslab1/SCDA/Apertures/JPL/offset_masks')\n",
    "\n",
    "secobs_fname = os.path.join(telap_dir, \"{0:s}_spiders_{1:04d}pix_2.5cm_offset.fits\".format(\n",
    "                            secobs_key_map[secobs_key], D))\n",
    "\n",
    "secobs_hdulist = fits.open(secobs_fname, \"readonly\")\n",
    "secobs = secobs_hdulist[0].data\n",
    "secobs_hdulist.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x145e6c510>"
      ]
     },
     "execution_count": 339,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJKCAYAAADweDtMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl0lPXd///Xe5ZMyAoIGJawKC4o3CKoqAckylKkBEEF\npRUUt6/2/LTV3kXooty9T3tbW9ujd4+9v60tLrd1FwWrVbCNwmlFxS8WqmhVUEBZFCWsySRz/f4g\nsRQhkGRmPtfyfJzDOWEyybySK5O85/X5zDXmeZ4AAACQPTHXAQAAAMKGAQsAACDLGLAAAACyjAEL\nAAAgyxiwAAAAsowBCwAAIMsSrgPsy8w4ZwQAAAgMz/PsQJf7rsHyPI9/Afp3yy23OM/AP45Z2P9x\nzIL3j2MWvH9tOWYt8d2ABQAAEHQMWAAAAFnGgIV2qaqqch0BrcQxCx6OWfBwzIIn28fMDrWGmE9m\n5vkpDwAAwMGYmbygbHIHAAAIOgYsAACALGtxwDKz35nZJjNbuc9lnc1skZm9Y2bPm1nHfd43x8z+\nYWarzWzsPpcPNbOVTe+7IzdfCgAAgD8cqsGaJ2ncfpfNlrTI87xjJb3Q9H+Z2QmSLpJ0QtPH3GVm\nzeuSv5J0hed5x0g6xsz2/5wAAACh0eKA5XneEkmf7XfxREn3Nr19r6RJTW+fJ+lBz/PSnuetlfSu\npGFm1l1Sqed5rzRd7759PgYAACB02rIH60jP8zY1vb1J0pFNb/eQtH6f662X1PMAl29ouhwAACCU\n2rXJvemcCpxXAQAAYB9tebHnTWZW4Xnexqblv81Nl2+QVLnP9Xppb3O1oentfS/fcLBPPnfu3C/e\nrqqq4mRtAADAF2pqalRTU3NY1z3kiUbNrK+khZ7nDWr6/22SPvU87ydmNltSR8/zZjdtcv+9pNO0\ndwlwsaT+nud5ZrZM0vWSXpH0B0l3ep73xwPcFicaBQAAgdDSiUZbbLDM7EFJIyV1MbN1km6WdKuk\nR8zsCklrJU2VJM/z3jSzRyS9KalB0jf2mZa+IekeSR0kPXOg4QoAACAseKkcAACANuClcgAAAPKI\nAQsAACDLGLAAAACyjAELAAAgyxiwAAAAsowBCwAAIMsYsAAAALKMAQsAACDLGLAAAACyjAELAAAg\nyxiwAAAAsowBCwAAIMsYsAAAALKMAQuBsGLFCpWUlKixsdF1FAAR0djYqJKSEq1YscJ1FASQeZ7n\nOsMXzMzzUx74i5kpFosxZAHIi3g8rkwmI/4u4WDMTJ7n2YHeR4OFwPjBD36gTCajmpoa11EAhFxN\nTY0ymYx+8IMfuI6CgKLBQqA8/fTTqq6u5hElgJwyMy1cuFATJkxwHQU+RoOF0JgwYYLMTBdffLHr\nKABC6uKLL5aZMVyhXRiwEDgXXXSRHnnkEdcxAITUI488oosuush1DAQcS4QIJDPTggULVF1d7ToK\ngBBZuHChJk6cyDYEHJaWlggZsBBIL774oqqqqvglCCCrzEw1NTUaOXKk6ygIAPZgIXRGjhypWCzG\nHgkAWTNhwgTFYjGGK2QFDRYCK5PJKB6P02IByAozU2Njo2IxugccHhoshFIsFlPfvn1VXl7OkAWg\nzTzPU3l5ufr27ctwhazhJwmB9thjj6m2tlbf+ta3XEcBEFDf+ta3VFtbq8cee8x1FIQIS4QIvGQy\nqYaGBlosAG1iZkokEkqn066jIGBYIkSopdNplZeXq0+fPq6jAAiYPn36qLy8nOEKWUeDhdAwM6XT\naSUSCddRAARAQ0ODkskk7TfajAYLkXDeeecplUq5jgEgIFKplM477zzXMRBSNFgIFTPT4sWLNWrU\nKNdRAPjYCy+8oNGjR9NeoV04kzsio1u3borFYtq4caPrKAB8rKKiQplMRps3b3YdBQHGgIVIMTPt\n3r1bhYWFrqMA8KE9e/aoQ4cOtFdoN/ZgIVK+/vWvq6ioyHUMAD5VVFSkr3/9665jIORosBBKZqb5\n8+dr0qRJrqMA8JEnn3xSkydPpr1CVtBgIXKWLFmiyZMna86cOa6jAPCJOXPmaPLkyVqyZInrKIgA\nGiyEVjweVyaT4ZEqAEl724ZYLKbGxkbXURASNFiIpKuvvlqS9MYbbzhOAsC15t8Dzb8XgFyjwUKo\n3XPPPZo5cyYtFhBxZqZ58+bpsssucx0FIUKDhchq/mX6jW98w20QAM403/8ZrpBPNFgIvdraWpWX\nl9NiARFlZtq2bZvKyspcR0HI0GAh0srKytSpUydVVla6jgIgzyorK9WpUyeGK+QdAxYi4e6779b6\n9etdxwCQZ+vXr9fdd9/tOgYiiAELkXD++efLzDRjxgzXUQDkyYwZM2RmOv/8811HQQSxBwuRweuP\nAdHC65Ii13ixZ6BJRUWFPM/Tpk2bXEcBkENHHnmkzEwbN250HQUhxoAF7MPMVFBQoLq6OtdRAORA\nKpVSfX09bTVyjmcRAvu49NJLVV9fr7ffftt1FABZ9vbbb6u+vl6XXnqp6yiIOBosRE5DQ4PKysrU\nr18//f3vf3cdB0AWnXjiiVqzZo1qa2uVSCRcx0HIsUQIHEDTHcN1DABZxP0a+cQSIXAAI0eOVDKZ\ndB0DQJYkk0mNHDnSdQxAEg0WIs7MtHz5cg0ZMsR1FADt8Prrr2vo0KG0V8grGizgIO6//34NHTrU\ndQwA7TR06FDdf//9rmMAX2DAQqRdcsklkqT/83/+j+MkANqq+f7bfH8G/IAlQkTe+PHj9dxzz6mx\nsdF1FABtEI/H9ZWvfEXPPPOM6yiIGJ5FCByCmemll17SiBEjXEcB0ApLlizRWWedxd4rOMEeLOAQ\nnnrqKZ111ln6xS9+4ToKgMP0i1/8QmeddZaeeuop11GAL6HBAprEYjF5nscjYSAgzExmpkwm4zoK\nIooGCzgMu3fvViqV0sknn+w6CoBDOPnkk5VKpbR7927XUYADosEC9sOZoAH/434KP6DBAlpp7ty5\nriMAOAjunwgCGixgPx9//LF69OjBo2PAp8xMH330kbp37+46CiKOBgtohe7du6uoqEgDBgxwHQXA\nfgYMGKCioiKGK/geDRZwEOzxAPyH+yX8hAYLaIOzzz5biUTCdQwAkj777DMlEgmdffbZrqMAh4UG\nC2iBmemVV17Rqaee6joKEGlme0sC/kbAT3ipHKCNysvL1blzZ61Zs8Z1FCDSzExlZWXatm2b6yjA\nFxiwgHYwMxUVFWnnzp2uowCRVFxcrMrKSq1evdp1FOBfsAcLaIfq6mrt2rVLH3/8sesoQOR8/PHH\n2rVrF8MVAocBCziEBQsW6Oabb1aPHj1cRwEip0ePHrr55ptdxwBajSVC4DCZmW677TZ95zvfcR0F\niISf/vSnmjVrFhvb4VvswQKyYO3aterXrx+/7IE8MTOtWbNGffv2dR0FOCD2YAFZ0PxL/v3333cb\nBIiA5vsZwxWCigELaIXbb79dRx99tOsYQOgdffTRuv32213HANqMAQtohRtvvFGS9P3vf99xEiC8\nmu9fzfc3IIh4HRCglTZu3KiKigq9/vrreuaZZ1zHAUJl/PjxevbZZ7Vx40bXUYB2YZM70AYlJSXa\nuXMnG96BLDMzFRcXa8eOHa6jAIfEJncgy3bs2KG+ffuqrKzMdRQgNMrKytS3b1+GK4QCAxbQRo8/\n/ri2b9/uOgYQGtu3b9fjjz/uOgaQFQxYQBsNGTJEiURCVVVVrqMAgVdVVaVEIqEhQ4a4jgJkBXuw\ngHZqWoN3HQMINO5HCCL2YAE5dMIJJ6ioqMh1DCCwioqKdMIJJ7iOAWQVDRaQBWamdevWqVevXq6j\nAIGyfv16VVZW0l4hkGiwgBz74Q9/qMrKStcxgMCprKzUD3/4Q9cxgKyjwQKyhD0kQOtxv0GQ0WAB\neTB06FClUint2bPHdRTA9/bs2aNUKqWhQ4e6jgLkBA0WkEVmJjNTJpNxHQXwtVgsJs/zaK8QaDRY\nQJ7ccccd8jxP8+fPdx0F8K358+fL8zzdcccdrqMAOUODBWTZ0qVLNWLECB6ZAwdhZlqyZImGDx/u\nOgrQLjRYQB4NHz5c8Xhc48aNcx0F8J1x48YpHo8zXCH0GLCAHLj66qv13HPPuY4B+M5zzz2nq6++\n2nUMIOdYIgRyxMw0b948XXbZZa6jAL5wzz33aObMmSyfIzRaWiJkwAJy5I033tDgwYP5YwI0MTOt\nWLFCJ510kusoQFawBwtw4KSTTlIikdCIESNcRwGcGzFihBKJBMMVIoMGC8gxzlQNcD9AONFgAQ4N\nHDhQRUVFSqfTrqMAeZdOp1VUVKSBAwe6jgLkFQMWkGPz58/X7t27dfnll7uOAuTd5Zdfrt27d3Py\nXUQOS4RAHhQWFqquro4lEkSOmfEanQgtnkUI+EBFRYUymYw2b97sOgqQF926dVMsFtPGjRtdRwFy\nggEL8Akz086dO1VUVOQ6CpBTu3btUnFxMa0tQo1N7oBPXHbZZSopKXEdA8i5kpISTrKLSGvzgGVm\nc8zs72a20sx+b2YpM+tsZovM7B0ze97MOu53/X+Y2WozG5ud+ECwzJs3T57n6dFHH3UdBciZRx99\nVJ7nad68ea6jAM60aYnQzPpK+pOkAZ7n1ZnZw5KekXSipE88z7vNzG6S1MnzvNlmdoKk30s6VVJP\nSYslHet5Xma/z8sSIUKvT58+qq2t1WeffeY6CpATnTp1UllZmT744APXUYCcysUSYa2ktKQiM0tI\nKpL0kaSJku5tus69kiY1vX2epAc9z0t7nrdW0ruSTmvjbQOB9sEHH+jzzz/X559/7joKkHXNP9sM\nV4i6Ng1YnudtlXS7pA+1d7D63PO8RZKO9DxvU9PVNkk6suntHpLW7/Mp1mtvkwVE0nXXXadOnTq5\njgFkXadOnXTddde5jgE416YBy8yOlvQtSX21d3gqMbNL9r1O01pfS+t9rAUisu68805J0m9/+1vH\nSYDsaf55bv75BqIs0caPO0XSXzzP+1SSzOwJSWdI2mhmFZ7nbTSz7pKaT/izQVLlPh/fq+myL5k7\nd+4Xb1dVVamqqqqNEQF/W7lypQYNGqRly5bp17/+tes4QLtcffXV+s1vfqOVK1e6jgLkTE1NjWpq\nag7rum3d5H6SpAe0d9P6Hkn3SHpFUh9Jn3qe9xMzmy2p436b3E/TPze5999/Rzub3BE1BQUFSqfT\nnCsIgWdmSiaTqq+vdx0FyJusb3L3PO8NSfdJek3S35ou/rWkWyWNMbN3JJ3T9H95nvempEckvSnp\nWUnfYJICpFmzZkmSXnrpJcdJgLZr/vlt/nkGwJncAecWLlyoiRMn0mIhsMxMCxYsUHV1tesoQF5x\nJnfAx6qrq2VmmjZtmusoQKtNmzZNZsZwBeyHBgvwgbq6OhUWFtJiIXDMTHv27FEqlXIdBcg7GizA\n51KplLp3766uXbu6jgIctq5du6p79+4MV8AB0GABPmJm2r59Oy8IDd/bsWOHSktLaV0RaTRYQEDE\nYjH2YiEQpk2bpliMPyHAwdBgAT7S2NioRCJBKwDfMzM1NDQoHo+7jgI4Q4MFBEQ8Hlffvn1VXl7O\nkAVf8jxP5eXl6tu3L8MV0AIaLMCHzEzxeFwNDQ2uowD/IpFIqLGxkQcAgGiwgMC54YYb1NjYqFde\necV1FOALr7zyihobG3XDDTe4jgL4Hg0W4FMPPfSQpk2bRlMA3zAzPfjgg7r44otdRwF8oaUGiwEL\n8DE2EsMveAIG8GUsEQIBNXHiRBUUFLiOAaigoEATJ050HQMIDBoswOfMTC+88ILOOecc11EQUX/6\n0580atQo2itgPzRYQID94Q9/0KhRo1zHQISNGjVKf/jDH1zHAAKFAQvwufHjx8vMNHXqVNdREEFT\np06VmWn8+PGuowCBwhIhEAD19fVKpVIs0SDvzEx1dXXsBQQOgCVCIOCa/7g99dRTjpMgSpp/3hiu\ngNZjwAIC4qWXXtKkSZP03e9+13UURMB3v/tdTZo0SS+99JLrKEAgsUQIBEg8Hlcmk2GpEDlnZorF\nYmpsbHQdBfAtlgiBkEin0yopKdExxxzjOgpC7JhjjlFJSYnS6bTrKEBg0WABAWRmymQyMjvgAyeg\nzTzPUywWoyUFDgMNFhAyY8aMUTKZdB0DIZRMJjVmzBjXMYDAo8ECAmjbtm3q2LEjLQOyzsz0+eef\nq7y83HUUwPdosICQKS8vV8eOHdW7d2/XURAivXv3VseOHRmugCygwQICzMyUTqeVSCRcR0HANTQ0\nKJlM0ooCrUCDBYTU5MmTlUqlXMdACKRSKU2ePNl1DCA0aLCAgDMzPf/882xMRpstWrRIY8eOpb0C\nWokGCwix5557TmPHjtWPf/xj11EQQD/+8Y81duxYPffcc66jAKFCgwWEgJmpoKBAdXV1rqMgYFKp\nlOrr62mvgDagwQJCbvr06aqvr9c777zjOgoC5J133lF9fb2mT5/uOgoQOjz1CAiB++67T6eeeqqO\nO+44mggctuOOO0533nmnrrvuOtdRgNChwQJCovmP5E033eQ4CYKg+eeE4QrIDfZgASHyySefqGvX\nrrRYOCQz05YtW9SlSxfXUYDAYg8WEBFdunRRaWmp+vfv7zoKfKx///4qLS1luAJyiAELCJlf/vKX\neu+991zHgI+99957+uUvf+k6BhBqDFhAyMyYMUOSdO211zpOAj9q/rlo/jkBkBvswQJCaPv27Sor\nK2MvFr7EzFRbW6vS0lLXUYDAYw8WEDGlpaXq1KmTevXq5ToKfKRXr17q1KkTwxWQBzRYQIiZmQoL\nC7V7927XUeBYhw4dtGfPHlpNIItaarAYsIAQi8Vi8jyPP6qQmcnMlMlkXEcBQoMlQiCidu3apVQq\npaFDh7qOAoeGDh2qVCqlXbt2uY4CRAYNFhABTY+yXMeAIxx/IDdosICIO+2005RKpVzHgAOpVEqn\nnXaa6xhA5NBgARFhZnrrrbd0/PHHu46CPFm9erUGDBhAewXkCA0WAN11110aMGCA6xjIowEDBuiu\nu+5yHQOIJBosIELYixMtHG8gt2iwAEiSqqqqlEwmXcdAHiSTSVVVVbmOAUQWDRYQMWam1157jVM3\nhNjy5ct1yimn0F4BOUaDBeALDzzwgE455RTdfffdrqMgB+6++26dcsopeuCBB1xHASKNBguIILO9\nD7i4v4UPxxbIHxosAP9i+/btSiaTOuOMM1xHQRadccYZSiaT2r59u+soQOQxYAERVFJSoptuukkv\nv/yy6yjIopdfflk33XSTSkpKXEcBIo8lQiDCzEw/+9nP9O1vf9t1FLTT7bffrn//939naRDIo5aW\nCBmwgAhbs2aNjjrqKP4oh4CZ6f3331e/fv1cRwEigz1YAA6oX79+SqVSGjx4sOsoaIfBgwcrlUox\nXAE+QoMFgDN+BxzHD3CDBgtAi04//XQVFBS4joE2KCgo0Omnn+46BoD90GAB0MaNG9W9e3dakAAy\nM3388ceqqKhwHQWIHBosAC2qqKhQcXGxjjvuONdR0ArHHXeciouLGa4AH6LBAvAFM1NZWZm2bdvm\nOgoOoby8XLW1tbSOgEM0WAAOy6hRo1RbW6utW7e6joIWbN26VbW1tRo1apTrKAAOggELwBcWL16s\nG264QUcccYTrKGjBEUccoRtuuEGLFy92HQXAQbBECOBLzEy/+tWvdM0117iOgv38z//8j6699lqW\nBgEf4EzuAFplwIAB+vDDD7Vz507XUbCf4uJi9e7dW2+99ZbrKEDkMWABaDUz04YNG9SjRw/XUdDk\no48+Us+ePWmvAJ9gkzuAVps7d6569uzpOgb20bNnT82dO9d1DACHgQYLwEGZmW699VbddNNNrqNE\n3k9+8hPNnj2b9grwEZYIAbTJBx98oL59+/JH3QfMTGvXrlWfPn1cRwHQhCVCAG3Sp08fFRYWatCg\nQa6jRNqgQYNUWFjIcAUECAMWgBbdeuutWrVqlZ588knXUSLpySef1KpVq3Trrbe6jgKgFVgiBHBI\nZnsbcO6f+cf3HvAvlggBtMumTZsUi8V07rnnuo4SKeeee65isZg2bdrkOgqAVqLBAnBYMpmM4vE4\nTUoemZkaGxsVi/FYGPAjGiwA7RaLxXT00UerrKzMdZRIKCsr09FHH81wBQQUDRaAVjEzbdmyRV26\ndHEdJbQ++eQTde3albYQ8DkaLABZk0gkdMEFF7iOEWoXXHCBEomE6xgA2oEGC0CrNT1qcx0jtPj+\nAsFAgwUgq0488UQVFRW5jhFKRUVFOvHEE13HANBONFgA2sTM9OGHH6qystJ1lNBYt26devfuTXsF\nBAQNFoCs+9GPfqTevXvr+eefdx0lFJ5//nn17t1bP/rRj1xHAZAFNFgA2oyzjGcP30sgeGiwAOTE\nKaecIknavXu34yTB1vz9a/5+Agg+BiwAbfbqq69qxowZKi4udh0l0IqLizVjxgy9+uqrrqMAyBKW\nCAG0m5np0Ucf1YUXXug6SuA89thjmjJlCkuDQAC1tETIgAWg3f7617/qzDPPZEhoAzPTX/7yF51x\nxhmuowBoJfZgAcipM844Q/F4XGPHjnUdJVDGjh2reDzOcAWEEA0WgKzwPE+xWIwWqxXMTJlM5otn\nEAIIFhosADnXPCTMmzfPcZJgaP4+MVwB4USDBSBrVq5cqX/7t3+jxToMZqa//e1vGjRokOsoANqI\nBgtAXgwaNEiJRELDhw93HcXXhg8frkQiwXAFhBgNFoCsMzP17NlT69evdx3Fd3r16qUNGzbQ8gEh\nQIMFIK8GDhyoDRs2KJ1Ou47iK+l0Whs2bNDAgQNdRwGQYwxYALJu5cqVmjJlilKplOsovpJKpTRl\nyhStXLnSdRQAOcYSIYCcqKurU2FhIUth+zAz7dmzh8ETCAmWCAHkXSqVUkVFhbp27eo6ii907dpV\nFRUVDFdARNBgAcgpM9P27dtVUlLiOoozO3bsUGlpKW0eEDI0WACcmTlzpkpLS13HcKq0tFQzZ850\nHQNAHtFgAcg5M9NDDz2kiy66yHWUvHv44Yd18cUX014BIdRSg8WABSDnXnnlFQ0bNiySQ4aZadmy\nZTrttNNcRwGQZQxYAJwzM23dulWdOnVyHSVvPvvsM3Xu3DmSgyUQBTnZg2VmHc3sMTN7y8zeNLNh\nZtbZzBaZ2Ttm9ryZddzn+nPM7B9mttrMxrb1dgEE0ze/+U117txZr732musoefHaa6+pc+fO+uY3\nv+k6CgAH2txgmdm9kl70PO93ZpaQVCzpe5I+8TzvNjO7SVInz/Nmm9kJkn4v6VRJPSUtlnSs53mZ\n/T4nDRYQYmZ7H+hF4X4epa8ViKqsN1hmVi5phOd5v5Mkz/MaPM/bJmmipHubrnavpElNb58n6UHP\n89Ke562V9K4kNiQAEdN8BvMrr7zScZLcav76OGM7EF1tarDMbLCk/yvpTUknSVou6VuS1nue16np\nOiZpq+d5nczsvyW97HneA03vu1vSs57nPb7f56XBAkIuCueE4txfQDTkYg9WQtIQSXd5njdE0k5J\ns/e9QtOk1NJv0PD+dgVwUCUlJTriiCPUo0cP11FyokePHjriiCMYroCIS7Tx49Zrb1v1atP/H5M0\nR9JGM6vwPG+jmXWXtLnp/RskVe7z8b2aLvuSuXPnfvF2VVWVqqqq2hgRgF/dc889qq6udh0jJz7+\n+GMtXLjQdQwAOVBTU6OamprDum57Nrm/JOlKz/PeMbO5koqa3vWp53k/MbPZkjrut8n9NP1zk3v/\n/dcDWSIEoiMWi2nq1Kl66KGHXEfJmosvvliPPPKIMpnMoa8MIPBych4sMztJ0t2SCiS9J2mmpLik\nRyT1lrRW0lTP8z5vuv53JV0uqUHSNz3Pe+4An5MBC4iI+vp6pVKpUO3FMjPV1dWpoKDAdRQAecCJ\nRgH4Us+ePVVXV6dPPvnEdZR269Kli1KplDZsOODuBwAhxIAFwLfC8Iy7KDwzEsCX5eRM7gCQDVdd\ndZVKS0u1atUq11HaZNWqVSotLdVVV13lOgoAH6HBAuBUJpNRPB7Xscceq7ffftt1nFY77rjj9M47\n76ixsVGxGI9ZgShhiRCAr/Xr109r165VJpP54iVmgsDzPMViMfXt21dr1qxxHQdAnrFECMDX1qxZ\no9GjRyuZTLqO0irJZFKjR49muALwJTRYAHzDzPTyyy9r2LBhrqMc0rJly3T66aezsR2IMBosAIHw\n8MMP6/TTT3cd47Ccfvrpevjhh13HAOBTDFgAfGPq1KkyM82cOdN1lBbNnDlTZqapU6e6jgLAp1gi\nBOArkyZN0sKFC9XY2Og6ykHF43FVV1frySefdB0FgEM8ixBAoJiZFi1apNGjR7uO8iWLFy/WmDFj\n2HsFgD1YAILl2Wef1ZgxY1zHOKAxY8bo2WefdR0DgM8xYAHwnXHjxsnMdOGFF7qO8i8uvPBCmZnG\njRvnOgoAn0u4DgAAB1JXV6eCggINGjRIK1eudB1HgwYN0qpVq1RfX+86CoAAYA8WAN+qrKzU+vXr\nfbHfyczUq1cvrVu3znUUAD7BHiwAgfTggw9KkubMmeM0R/PtN+cBgEOhwQLga1u2bFG3bt2ctlhm\nps2bN6tr167OMgDwHxosAIHVtWtXlZSU6JhjjnFy+8ccc4xKSkoYrgC0Cg0WgEAwM2UyGZkd8MFi\nTniep1gs5os9YAD8hwYLQOB95StfUTKZzOttJpNJfeUrX8nrbQIIBxosAIFhZlqyZImGDx+e89ta\nunSpRowYQXsF4KB4qRwAodC5c2cVFxfn5VQJlZWV2rlzp7Zu3Zrz2wIQTCwRAgiFrVu3av369Uqn\n0wd8/6xZs3TyySfLzA757+STT9asWbMO+HnS6bTWr1/PcAWgzRiwAATK5MmTlUqlvmixfvrTn34x\nNP30pz/VihUrDuvzrFix4ksfK0nr1q1TKpXS5MmTc/Y1AAg/lggBBE4+nknI7yIAh8ISIYDQ+NrX\nvhaq2wEQTjRYAAIjkUiosbExb7cXj8fV0NCQt9sDECw0WAACb86cOXkdriSpsbHR+esgAggmGiwA\nvlddXa2hUJd/AAAgAElEQVSnn37a2e1PmDBBCxcudHb7APyJBgtAYM2ZM8fpcCVJTz/9NE0WgFah\nwQLgW/nec3Uo7MkCsC8aLACB87Wvfc1Xw5W0d08Wzy4EcDhosAD4Uj7OddVW/J4CINFgAQiYW2+9\n1XWEFvk9HwD3aLAA+I6f26tm/K4CQIMFIDCaXxPQ74KSE4AbNFgAfCUI7VUzfl8B0UaDBQAAkEcM\nWAB8Y9asWa4jtErQ8gLIHwYsAL6xaNEi1xFaJWh5AeQPe7AA+EaQ9l8143cWEF3swQIAAMgjBiwA\nAIAsY8ACAADIMgYsAACALGPAAgAAyDIGLAAAgCxjwAIAAMgyBiwAAIAsY8AC4BuDBw92HaFVgpYX\nQP4wYAHwjTFjxriO0CpBywsgf3ipHAC+EqSXy+H3FRBtvFQOAABAHjFgAfCV2267zXWEwxKUnADc\nYIkQgO8EYZmQ31UAWCIEECj/9V//5TpCi/yeD4B7NFgAfMnPLRa/pwBINFgAAmjatGmuIxyQX3MB\n8BcaLAC+lUgk1NjY6DrGF+LxuBoaGlzHAOATNFgAAqmhoUGzZ892HUOSNHv2bIYrAIeNBguA71VX\nV+vpp592dvsTJkzQwoULnd0+AH+iwQIQaAsXLnTWZM2ePZvhCkCr0WABCIx878lizxWAltBgAQiF\nhoaGvD2Lb9q0aQxXANqMBgtA4OTjHFn8LgJwKDRYAELj/PPPl5npww8/lOd5WXtNwNtuu02e5+nD\nDz+Umen888/PyucFEE00WAACxcxUX1+vZDL5pffNmjVLixYt0ooVKw75eQYPHqwxY8YccEBLp9Mq\nKCigxQLQopYaLAYsAIHRuXNnFRcXa926dTm/rcrKSu3cuVNbt27N+W0BCCaWCAGEwmeffaYHH3ww\nL7f14IMP6rPPPsvLbQEIHwYsAIEwbtw4xeNxDR8+PC+3N3z4cMXjcY0bNy4vtwcgXFgiBBAIZqZM\nJpOXZxA28zxPsViMvVgADoglQgCBVlpaqv79++d1uJL2/vLs37+/SktL83q7AIKPBguAr23ZskXd\nunVz2iKZmTZv3qyuXbs6ywDAf2iwAATS0qVL1a1bN2evQ9hs9uzZ6tatm5YuXeo0B4DgoMEC4FuV\nlZVav369L/ZAmZl69eqVl1NEAAgGGiwAgZNOp7V+/XoNHDjQdRRJ0sCBA7V+/Xql02nXUQAEAAMW\nAF9KpVK64IILtHLlStdRJEkrV67UBRdcoFQq5ToKgABgiRCA7/zxj3/Uueee64ulwf2ZmZ599lnO\njwWgxSXCRL7DAMChnHvuuVq0aJHrGAe0aNEijRkzxpfDHwD/YIkQgK9MmjRJsVhMo0ePdh3lgEaP\nHq1YLKZJkya5jgLAxxiwAPjKggULNGPGDNcxWjRjxgwtWLDAdQwAPsYeLAC+8cgjj+iiiy4KxPKb\nmenhhx/W1KlTXUcB4EhLe7AYsAD4hpnp5Zdf1rBhw1xHOaRly5bp9NNPD8QwCCA3OA8WAN8bM2aM\n4vF4IIYrSRo2bJji8bjGjBnjOgoAH+JZhACc69evn9auXatMJuM6Squk02nFYjH169dPa9ascR0H\ngI/QYAFwKpPJaO3atTr22GNldsCm3bfMTMcee2wgh0MAucWABcCpa665RpL0+OOPO07SNs25m78O\nAJDY5A7AMTPT9u3bVVJS4jpKm+3YsUOlpaVseAcihmcRAvClnj17qq6uTp988onrKO3WpUsXpVIp\nbdiwwXUUAHnCgAXAd+rr65VKpULV+piZ6urqVFBQ4DoKgDzgNA0AfKewsFAXXXSR6xhZddFFF6mw\nsNB1DAA+QIMFIO+efvppVVdXh6q9amZmWrhwoSZMmOA6CoAcY4kQgK906dJFBQUF+uijj1xHyboe\nPXqovr4+FPvKALSMAQuAb0ThGXdheGYkgENjDxYAX1i1apVKS0t1xRVXuI6SU1dccYVKS0u1atUq\n11EAOEKDBSBvms/UHoX7eZS+ViCqaLAAOPetb31LkvTqq686TpIfzV9n89cNIFposADkhZlp69at\n6tSpk+soefPZZ5+pc+fOtFhASNFgAXDqlVdekaRIDVfSP7/e5q8fQHQwYAHIuWHDhumhhx5yHcOJ\nhx56SMOGDXMdA0CeMWAByKnLL79ckkJ31vbD1fx1N38fAERDu/ZgmVlc0muS1nueV21mnSU9LKmP\npLWSpnqe93nTdedIulxSo6TrPc97/gCfjz1YQMhwTqhonPsLiKJc7sH6pqQ3JTX/1pgtaZHnecdK\neqHp/zKzEyRdJOkESeMk3WVmtGdAyHXv3l1dunSJ9HAlSSUlJerSpYu6d+/uOgqAPGnzkGNmvSSN\nl3S3pObpbaKke5vevlfSpKa3z5P0oOd5ac/z1kp6V9Jpbb1tAP5XV1enjRs3asuWLa6j+MKWLVu0\nceNG1dXVuY4CIA/a0yL9QtJ3JGX2uexIz/M2Nb29SdKRTW/3kLR+n+utl9SzHbcNwOemT5/+xck2\nsZeZafr06a5jAMiDRFs+yMwmSNrsed7/M7OqA13H8zzPzFracMBmBCCkBg0apFWrVqm+vt51FF+p\nq6tTQUGBBg0apJUrV7qOAyCH2jRgSTpT0kQzGy+pUFKZmd0vaZOZVXiet9HMukva3HT9DZIq9/n4\nXk2XfcncuXO/eLuqqkpVVVVtjAjAlVWrVqlnz55KJpOuo/hKMplUz549eY1CIKBqampUU1NzWNdt\n95nczWykpH9vehbhbZI+9TzvJ2Y2W1JHz/NmN21y/7327rvqKWmxpP77P2WQZxECwZdMJjVs2DAt\nXbrUdRTfGj58uJYtW6Z0Ou06CoB2yMeZ3JunolsljTGzdySd0/R/eZ73pqRHtPcZh89K+gaTFBA+\nK1euVENDA8PVISxdulQNDQ0sEwIhxmsRAsgaM9Pvfvc7zZw503UU35s3b54uv/xyzo0FBFhLDRYD\nFoCs8DxPsViMgaEVzEyZTIZnWwIBxYs9A8i5ZDKpMWPGuI4RKGPGjOGJAEBI0WABaLe//vWvOvPM\nM2mv2sDM9Je//EVnnHGG6ygAWokGC0BOnXnmmXr00UddxwikRx99VGeeeabrGACyjAELQLtceuml\nMjNdeOGFrqME0oUXXigz06WXXuo6CoAsauuJRgFAp556ql577TXt2rXLdZRA27lzp4qKivTmm2/q\n1VdfdR0HQBawBwtAmzU/+437bfvxvQSChz1YALLuxz/+sSTpueeec5wkHJq/j83fVwDBRoMFoE3M\nTB9++KEqKysPfWUclnXr1ql37960WEBA0GAByKqBAweqQ4cODFdZVllZqQ4dOmjgwIGuowBoJxos\nAK3W9KjNdYzQ4vsLBAMNFoCsSSaTOuuss1zHCLWzzjqLM7wDAUeDBaBVzExbtmxRly5dXEcJrU8+\n+URdu3alxQJ8jgYLQFb0799fpaWlDFc51qVLF5WWlqp///6uowBoIxosAIclk8koHo/TquSRmamx\nsVGxGI+FAT+iwQLQLps3b1YymdS4ceNcR4mUcePGKZlMavPmza6jAGglGiwAh8RZxt3hew/4Fw0W\ngDa74447JEnz5893nCSamr/vzccBQDDQYAFoUYcOHdS/f3+tXLnSdZTIGjRokN59913t3r3bdRQA\n+2ipwWLAAnBQH3zwgfr27cvylA+YmdauXas+ffq4jgKgCUuEANqkb9++uvXWW13HgKRbb71Vffv2\ndR0DwGFiwAJwQP/xH/8hSbrpppscJ4H0z+PQfFwA+BtLhAAOyMy0YcMG9ejRw3UUNPnoo4/Us2dP\nlmwBn2CJEECrDBgwQEVFRQxXPtOjRw8VFRVpwIABrqMAOAQGLABfsnr1at1+++2uY+AAbr/9dq1e\nvdp1DACHwIAF4F/ceOONkqRrrrnGcRIcSPNxaT5OAPwp4ToAAP8YPXq0XnjhBX366aeuo6AFn376\nqY444gj97W9/0+LFi13HAXAAbHIH8AUzU1lZmbZt2+Y6Cg6hvLxctbW1bHgHHGKTO4BDKikp0bHH\nHstwFRDbtm3Tscceq5KSEtdRABwADRYAbdy4Ud27d6cNCSAz08cff6yKigrXUYDIocEC0KLJkycr\nmUy6joE2SCaTmjx5susYAPZDgwWg+VGY6xhoI44f4AYNFoCDKiws1EknneQ6BtrhpJNOUmFhoesY\nAPZBgwVE2Jo1a3TUUUfRfoSAmen9999Xv379XEcBIoMGC8ABHXXUUfrZz37mOgay4Gc/+5mOOuoo\n1zEANGHAAiLqBz/4gSTp29/+tuMkyIbm49h8XAG4xRIhEEE7duxQ586dNXToUP31r391HQdZcsYZ\nZ2j58uXaunUr58cC8qClJUIGLCCCzPb+PuD+Fj4cWyB/2IMF4Au///3vJUm/+c1vHCdBLjQf1+bj\nDMANGiwgYsxMr732moYOHeo6CnJk+fLlOuWUU2ixgByjwQIgSTr77LOVSCQYrkJu6NChSiQSOvvs\ns11HASKLBguIEM74HS0cbyC3aLAA6Fe/+pXrCHCA4w64QYMFRISZ6a233tLxxx/vOgryZPXq1Row\nYAAtFpAjNFhAxA0bNkwFBQUMVxFz/PHHq6CgQMOGDXMdBYgcGiwgAtiLE20cfyA3aLCAiNqzZ48K\nCws1ZMgQ11Hg0JAhQ1RYWKg9e/a4jgJEBgMWEGLTp09XXV2d5s+f7zoKHJo/f77q6uo0ffp011GA\nyGCJEAgxM1NhYaF2797tOgoc69Chg/bs2cNSIZBFLBECEdS5c2f17NmT4QqSpN27d6tnz57q3Lmz\n6yhAJNBgASG0fft2lZWV0VbgS8xMtbW1Ki0tdR0FCDwaLCBiysrKdM0117iOAR+65pprVFZW5joG\nEHo0WEDI3Hfffbr00ktpr3BQZqZ7771XM2bMcB0FCLSWGiwGLCBkysrK1K1bN7377ruuo8Cn+vfv\nr82bN6u2ttZ1FCDQGLCAiPjkk0/UtWtX2isckplpy5Yt6tKli+soQGCxBwuIiK5du2rWrFmuYyAA\nZs2apa5du7qOAYQWDRYQEv/93/+t66+/nvYKh83MdOedd+q6665zHQUIpJYarES+wwDIvhkzZuj+\n++/X22+/7ToKAuTtt9/Wcccdp1dffVX33Xef6zhAqNBgASFgZiooKFBdXZ3rKAiYVCql+vp6mk+g\nDdiDBYTY888/L0m65ZZbHCdBEDX/3DT/HAHIDhosIODMTM8//7zGjBnjOgoCatGiRRo7diwtFtBK\nNFhASJ1//vmKxWIMV2iXMWPGKBaL6fzzz3cdBQgNGiwgwMxM6XRaiQTPV0H7NDQ0KJlM0mIBrUCD\nBYRQp06dVFlZyXCFrEgkEqqsrFSnTp1cRwFCgQYLCKBt27apY8eOtA3IOjPT559/rvLyctdRAN+j\nwQJCZsqUKYrH465jIITi8bimTJniOgYQeDRYQACZmTKZjMwO+MAJaDPP8xSLxWhHgcNAgwWERCaT\nUWlpqfr3789whZwwM/Xv31+lpaXKZDKu4wCBRYMFBEg8Hlcmk6FdQM6ZmWKxmBobG11HAXyLBgsI\ngSVLliiTyWjOnDmuoyAC5syZo0wmoyVLlriOAgQSDRYQEGamJ598Uuedd57rKIiIp556SpMmTaIx\nBQ6ipQaLAQsIgPr6eqVSKf7QIe/MTHV1dSooKHAdBfAdlgiBgCssLOSp83BiypQpKiwsdB0DCBwa\nLMDnnnnmGX31q1+lvYIzZqY//OEPGj9+vOsogK+wRAgEmJnphRde0DnnnOM6CiLqT3/6k0aNGsWQ\nD+yHJUIgoM477zzFYjGGKzh1zjnnKBaL8QQLoBVosAAfMzM1NDTwsjhwrrGxUYlEghYL2AcNFhBA\nDz30kCQxXMEXmn8Om38uAbSMAQvwoRtvvFHTpk3TsmXLXEcBvrBs2TJNmzZNN954o+sogO+xRAj4\nkJkpHo+roaHBdRTgXyQSCTU2NrJUCIglQiBQ+vXrp7KyMqXTaddRgC9Jp9MqKytTv379XEcBfI0G\nC/ARNhIjKHgCBkCDBQRGQUGBJkyY4DoGcEgTJkzg5XOAFtBgAT5iZtq+fbtKSkpcRwFatGPHDpWW\nltK2ItI4kzsQAD169FA6ndaWLVtcRwEOS9euXZVMJvXRRx+5jgI4wYAF+FxdXZ0KCwtpAxA4ZqY9\ne/YolUq5jgLkHXuwAJ/r0KGDLr74YtcxgFa7+OKL1aFDB9cxAN+hwQIcW7hwoSZOnEh7hcAyMy1Y\nsEDV1dWuowB5RYMF+NT3v/99TZw4US+++KLrKECbvfjii5o4caK+//3vu44C+AYNFuBQQUGB0uk0\n7RUCz8yUTCZVX1/vOgqQNzRYgA+tWrVK6XRaV111lesoQLtdddVVSqfTWrVqlesogC/QYAGOmJnu\nvvtuXXHFFa6jAFnx29/+VldeeSWNLCKDBgvwmeuvv16SGK4QKs0/z80/30CU0WABDpiZPvvsM3Xs\n2NF1FCCrPv/8c3Xq1IkWC5FAgwX4SJ8+fdSxY0eGK4RS8892nz59XEcBnGLAAvLsww8/1K9//WvX\nMYCc+fWvf60PP/zQdQzAKQYsII9mzpwpM9OUKVNcRwFyZsqUKTIzzZw503UUwBn2YAF5ZGbauXOn\nioqKXEcBcmrXrl0qLi5mLxZCjRd7BnygoqJCmUxGmzdvdh0FyItu3bopFotp48aNrqMAOcGABThW\nWFiouro6Hs0jcsxMqVRKe/bscR0FyLqsP4vQzCrN7M9m9nczW2Vm1zdd3tnMFpnZO2b2vJl13Odj\n5pjZP8xstZmNbduXAgTPu+++q7q6Ol1yySWuowB5d8kll6iurk7vvvuu6yhAXrWpwTKzCkkVnuet\nMLMSScslTZI0U9InnufdZmY3Serked5sMztB0u8lnSqpp6TFko71PC+z3+elwULoDBo0SO+99562\nbdumZDLpOg6QV+l0WuXl5Tr66KO1cuVK13GArMr5EqGZPSnpl03/Rnqet6lpCKvxPO94M5sjKeN5\n3k+arv9HSXM9z3t5v8/DgIXQaboDuo4BOMX9AGGU0xONmllfSSdLWibpSM/zNjW9a5OkI5ve7iFp\n/T4ftl57mywg1JLJpIYPH+46BuDc8OHDaXARKe0asJqWBx+X9E3P87bv+76mKqqlhys8lEGovfHG\nG2poaNCSJUtcRwGcW7JkiRoaGvTGG2+4jgLkRaKtH2hmSe0dru73PO/Jpos3mVmF53kbzay7pObn\no2+QVLnPh/dquuxL5s6d+8XbVVVVqqqqamtEwKnBgwdr3rx5rmMAvjFv3jwNHjyYpUIEVk1NjWpq\nag7rum3d5G6S7pX0qed5N+xz+W1Nl/3EzGZL6rjfJvfT9M9N7v3333DFHiyExTe+8Q396le/4g8J\nsB8z07XXXqu77rrLdRSg3bK+yd3Mhkt6SdLf9M+lvjmSXpH0iKTektZKmup53udNH/NdSZdLatDe\nJcXnDvB5GbAQColEQqNHj9Yf//hH11EAXxk3bpwWL16shoYG11GAduNEo0AeLV26VCNGjKC9Ag7C\nzLRkyRKeAILAy+mzCAH805133qkRI0boiSeecB0F8K0nnnhCI0aM0J133uk6CpAzNFhAFpmZzEyZ\nTObQVwYiLBaLyfM8ml4EGg0WkAennHKKCgoKtGvXLtdRAN/btWuXCgoKdMopp7iOAuQEDRaQJZyp\nGmg97jcIMhosIMf+8z//03UEILC4/yCMaLCALDAzrVu3Tr169XIdBQiU9evXq7KykhYLgUSDBeTQ\niSeeqA4dOjBcAW3Qq1cvdejQQSeeeKLrKEBW0WAB7cQeEqD9uB8hiGiwgBxJJpMaOXKk6xhA4I0c\nOVLJZNJ1DCBraLCANnr99dc1dOhQHnUDWWJmWr58uYYMGeI6CnBYeKkcIAf69eunTz/9VLW1ta6j\nAKFQVlamI444QmvWrHEdBTgsDFhAlpWUlGjnzp20V0CWmZmKi4u1Y8cO11GAQ2IPFpBFmzZt0s6d\nO3Xuuee6jgKEzrnnnqudO3dq06ZNrqMA7cKABbRSRUWFvve97+mZZ55xHQUInWeeeUbf+973VFFR\n4ToK0C4sEQKt8POf/1zf/va3WRoEcszMdPvtt+vGG290HQU4KPZgAVliZnrvvfd01FFHuY4ChNr7\n77+vo48+mgcz8DX2YAFZsHbtWkliuALyoPl+1ny/A4KGAQs4TP369dNtt93mOgYQGbfddpv69evn\nOgbQJgxYwGG45ZZbJEnf+c53HCcBoqP5/tZ8/wOChD1YwCFMnDhRCxcu1EcffaTu3bu7jgNEyscf\nf6wePXqwFwu+xCZ3oB3MTEVFRdq5c6frKEAkFRcXq7KyUqtXr3YdBfgXDFhAG5WXl6tz5868dAfg\nmJmprKxM27Ztcx0F+EJLA1Yi32GAIKmtrdXixYtdxwAg8bqfCBQ2uQMHcc455ygej+vUU091HQWI\nvK1btyoej+ucc85xHQU4LCwRAgfRVP26jgFgH9wv4SecaBRopeLiYh1//PGuYwDYz/HHH6/i4mLX\nMYBDosEC9sPTwgF/MzNOmwJfoMECWqFHjx6c2BDwsVtuuUU9evRwHQNoEQ0WsB/2eAD+x/0UfkCD\nBRyGuro6FRYWavDgwa6jADiEwYMHq7CwUHV1da6jAAdEgwU0icVi8jyPR8VAQJiZzEyZTMZ1FEQU\nDRZwCAsWLJDnefr5z3/uOgqAw/Tzn/9cnudpwYIFrqMAX0KDBWjvo5CXXnpJI0aMcB0FQCssWbJE\nZ511Fs0znKDBAlowfvx4xWIxhisggEaMGKFYLKbx48e7jgL8CwYsRN6zzz6rK6+80nUMAG105ZVX\n6tlnn3UdA/gXLBEi0v73f/9X06dPZ3kBCDgz0/33369LLrnEdRRESEtLhAxYiDQz0/LlyzVkyBDX\nUQC0w+uvv66hQ4fyYAl5xR4s4ACqqqqUSCQYroAQGDJkiBKJhKqqqlxHASTRYCHCOBM0ED7cr5FP\nNFjAPhoaGlRUVKQTTjjBdRQAWXbCCSeoqKhIDQ0NrqMg4hiwEDlXXnmldu/erSeeeMJ1FABZ9sQT\nT2j37t08MxjOsUSIyDEzFRQU8BpmQEilUinV19ezVIic41mEQJOKigp5nqdNmza5jgIgh4488kiZ\nmTZu3Og6CkKMAQuQtGfPHnXo0IFHtUBEmJl2796twsJC11EQUmxyByQVFRVp+vTprmMAyJPp06er\nqKjIdQxEFA0WIuGJJ57QBRdcQHsFRIyZ6fHHH9f555/vOgpCiCVCRF7nzp1VXFysdevWuY4CII8q\nKyu1c+dObd261XUUhBADFiKttrZW5eXltFdARJmZtm3bprKyMtdREDLswUKklZeX69prr3UdA4Aj\n1157rcrLy13HQMTQYCHU7rnnHs2cOZP2Cog4M9O8efN02WWXuY6CEKHBQiRde+21mjlzplasWOE6\nCgDHVqxYoZkzZ9JmI29osBBa8XhcmUyG9gqApL1tQywWU2Njo+soCAkaLETO0qVLlclkNHv2bNdR\nAPjE7NmzlclktHTpUtdREAE0WAglM9P8+fM1adIk11EA+MiTTz6pyZMn02wjK2iwECmXXHKJzIzh\nCsCXTJo0SWamSy65xHUUhBwNFkKH1x8D0BJelxTZwolGERndunVTLBbTxo0bXUcB4GMVFRXKZDLa\nvHmz6ygIMJYIERlbtmzRAw884DoGAJ974IEHtGXLFtcxEGIMWAiNSZMmKRaLadSoUa6jAPC5UaNG\nKRaLsVcTOcMSIULDzJROp5VIJFxHARAADQ0NSiaT7MVCm7FEiNDr2LGjevfuzXAF4LAlEgn17t1b\nHTt2dB0FIUSDhcBLJpNqaGjgUSiANjEzJRIJpdNp11EQMDRYCK3ly5eroaFB119/vesoAALq+uuv\nV0NDg5YvX+46CkKEBguB1q9fP23dulWff/65zA74IAIAWuR5njp27KjOnTtrzZo1ruMgQDgPFkIp\nk8koHo+zNAggK8xMjY2NisVY3MHhYYkQoZRMJvXVr37VdQwAIfHVr35VyWTSdQyEBA0WAunFF19U\nVVUV7RWArDIz1dTUaOTIka6jIABosBA6VVVVWrBggesYAEJmwYIFqqqqch0DIcCAhcCZNm2azEzV\n1dWuowAImerqapmZpk2b5joKAo4BC4Hz8MMPa+rUqa5jAAipqVOn6uGHH3YdAwHHHiwEytNPP63q\n6mr2XgHIKTPTwoULNWHCBNdR4GPswUIo3Hzzzaqurtaf//xn11EAhNyf//xnVVdX6+abb3YdBQFF\ng4XAMDPFYjE1Nja6jgIgAuLxuDKZDI05DooGC4G3YsUKFRcXq76+3nUUABFRX1+v4uJirVixwnUU\nBBANFgAAQBvQYAEAAOQRAxYAAECWMWABAABkGQMWAABAljFgAQAAZBkDFgAAQJYxYAEAAGQZAxYA\nAECWMWABAABkGQMWAABAljFgAQAAZBkDFgAAQJYxYAEAAGQZAxYAAECWMWABAABkGQMWAABAljFg\nAQAAZBkDFgAAQJYxYAEAAGQZAxYAAECWMWABAABkGQMWAABAljFgAQAAZBkDFgAAQJYxYAEAAGRZ\nXgcsMxtnZqvN7B9mdlM+bxsAACBf8jZgmVlc0i8ljZN0gqRpZjYgX7eP3KipqXEdAa3EMQsejlnw\ncMyCJ9vHLJ8N1mmS3vU8b63neWlJD0k6L4+3jxzgl0jwcMyCh2MWPByz4AnygNVT0rp9/r++6TIA\nAIBQyeeA5eXxtgAAAJwxz8vP3GNmp0ua63neuKb/z5GU8TzvJ/tchyEMAAAEhud5dqDL8zlgJSS9\nLWmUpI8kvSJpmud5b+UlAAAAQJ4k8nVDnuc1mNn/J+k5SXFJv2W4AgAAYZS3BgsAACAqfHMmd05C\n6j9mVmlmfzazv5vZKjO7vunyzma2yMzeMbPnzazjPh8zp+kYrjazse7SR5uZxc3s/5nZwqb/c8x8\nzDCVYYcAAAOkSURBVMw6mtljZvaWmb1pZsM4Zv7WdAz+bmYrzez3ZpbimPmLmf3OzDaZ2cp9Lmv1\nMTKzoU3H+R9mdsfh3r4vBixOQupbaUk3eN7/3879hFhZxWEc/z5Qgqa0CcyYgRlihIqChiBLJAgX\nEWGtyoUiIW2MsE1QLVq3iWpjkJmpxIBUiEEEUpsg6A9aQZOLIskRZoz+Wa0SnxbnXHy9xDC3LneO\n+XxW7/29750Z7jNz5vfee87xLcAG4PGay9PAMdvrgQ/qYyTdDDxCyfA+YI+kJn7HrkC7gVkurt5N\nZm17GXjP9k3AbcBJklmzJE0AjwHTtm+lTHvZSjJrzX7K6901SEa9yeuvADttTwFTkvq/5j9qJeBs\nQtog2/O2v6jHfwDfUPYu2wIcqJcdAB6qxw8CM7b/sn0K+JaSbYyQpDHgfuA1oDdAJLNGSboW2GT7\ndSjzVW3/RjJr2TnKDeiquoBrFWXxVjJriO2PgF/6yoNkdKekdcAa25/W6w52nrOoVhqsbELauHrH\ndjvwCbDW9kI9tQCsrcc3ULLrSY7L40XgKeBCp5bM2jUJ/Chpv6TjkvZKuoZk1izbPwMvAD9QGqtf\nbR8jmV0OBs2ov36GJWbXSoOVmfYNk7QaeBvYbfv37jmXVRKL5ZdsR0jSA8BZ2ye4+O7VJZJZc64C\npoE9tqeBP6kfW/Qks7ZIuhF4Epig/ANeLWlb95pk1r4lZPSftNJgnQHGO4/HubRjjGUi6WpKc3XI\n9pFaXpB0fT2/Djhb6/05jtVajM7dwBZJ3wMzwL2SDpHMWjYHzNn+rD5+i9JwzSezZt0BfGz7J9vn\ngXeAu0hml4NBxsK5Wh/rqy8pu1YarM8pE8cmJK2gTDQ7usw/0xWvTvDbB8zafqlz6iiwox7vAI50\n6lslrZA0CUxRNpSNEbH9rO1x25OUSbcf2t5OMmuW7XngtKT1tbQZ+Bp4l2TWqpPABkkr6zi5mbKo\nJJm1b6CxsP59nqsrewVs7zxnUSPbaHQx2YS0WRuBbcBXkk7U2jPA88BhSTuBU8DDALZnJR2mDDTn\ngV3ORmvLrff6J7O2PQG8WW8wvwMepYyFyaxBtr+UdJDy5sAF4DjwKrCGZNYMSTPAPcB1kk4Dz/Hv\nxsJdwBvASspq3/eX9P2TcURERMRwtfIRYURERMT/RhqsiIiIiCFLgxURERExZGmwIiIiIoYsDVZE\nRETEkKXBioiIiBiyNFgRERERQ5YGKyIiImLI/gaoC5QxrG0fpwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x145647b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(secobs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pad the secondary obscuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 340,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_shift = int(round(D*float(pad)/1000))\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)\n",
    "XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])\n",
    "len(XYshifts_allowed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "padded_secobs_accum = np.ones(secobs.shape)\n",
    "for (xshift,yshift) in XYshifts_allowed:\n",
    "    secobs_shifted = np.roll(np.roll(secobs, yshift, 0), xshift, 1) \n",
    "    padded_secobs_accum *= secobs_shifted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#plt.figure(figsize=(10,10))\n",
    "#plt.imshow(padded_secobs_accum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine with annular stop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 343,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "L = secobs.shape[0]\n",
    "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",
    "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.round(padded_secobs_accum)\n",
    "Lyot_stop[inside_ann] = 0.\n",
    "Lyot_stop[outside_ann] = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x14636e390>"
      ]
     },
     "execution_count": 344,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAK5CAYAAABT3oOfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3V+o7edd5/HP18Y4re0Yi5C0TYYGSaCREazaqjMizJRS\nRNJetREswQlz0xmtcyE2XngtgmhhqOCf1rTYzIQqIYWijVUYQTB1WrU2zbSVCfZEcuL/qhdDQp+5\n2M9u9tnd/9b/5/f7vV5X+6y99znrrN+znue9n/1bv1WttQAAAMnXHfoOAADAKMQxAAB04hgAADpx\nDAAAnTgGAIDupkPfgZOqyqUzAADYudZanXW7nWMAAOjEMQAAdOIYAAA6cQwAAJ04BgCAThwDAEAn\njgEAoBPHAADQiWMAAOjEMQAAdOIYAAC6C+O4qt5fVder6jMnbntlVT1eVZ+vqo9X1S0nPvdgVX2h\nqp6qqjefuP07q+oz/XPv3c1/BQAANnPZzvEHkrzl1G3vSfJ4a+3uJJ/of05V3ZPkHUnu6d/zvqqq\n/j2/lOSB1tpdSe6qqtN/JwAAHNyFcdxa+4Mkf3/q5nuTPNQ/fijJ2/rHb03ycGvt+dba00m+mOSN\nVfWqJK9orT3Rv+6DJ74HAACGsc45x7e21q73j68nubV//Ook10583bUkrznj9mf67QAAMJSNXpDX\nWmtJ2pbuCwAAHNQ6cXy9qm5Lkn7KxHP99meS3HHi627P0Y7xM/3jk7c/s8a/CwAAO7VOHD+W5P7+\n8f1JHj1x+31VdXNV3ZnkriRPtNaeTfLlqnpjf4HeO098DwAADOOmiz5ZVQ8n+YEk31JVX0ryM0l+\nNskjVfVAkqeTvD1JWmtPVtUjSZ5M8kKSd/XTLpLkXUl+PclLk3ystfbb2/+vAADAZurFfj28qhrn\nzgAAMFuttTrrdu+QBwAAnTgGAIBOHAMAQCeOAQCgE8cAANCJYwAA6MQxAAB04hgAADpxDAAAnTgG\nAIBOHAMAQCeOAQCgE8cAANCJYwAA6MQxAAB04hgAADpxDAAAnTgGAIBOHAMAQCeOAQCgE8cAANCJ\nYwAA6MQxAAB04hgAADpxDAAAnTgGAIBOHAMAQCeOAQCgE8cAANCJYwAA6MQxAAB04hgAADpxDAAA\nnTgGAIBOHAMAQCeOAQCgE8cAANCJYwAA6MQxAAB04hgAADpxDAAAnTgGAIBOHAMAQCeOAQCgE8cA\nANCJYwAA6MQxAAB04hgAADpxDAAAnTgGAIBOHAMAQCeOAQCgE8cAANCJYwAA6MQxAAB04hgAADpx\nDAAAnTgGAIBOHAMAQCeOAQCgE8cAANCJYwAA6MQxAAB04hgAALqbDn0HALhRa22nf39V7fTvB5gy\ncQywJ7uO3qu66v0Q0cASiWOADYwSvLuwyv9NSANzIY4BrmjOIbypsx4bwQxMkTgG6MTvdl32eIpn\nYETiGFg0QXw4Jx97oQyMQhwDiyCCx3be8RHNwL6JY2C2BPH02V0G9k0cA7MghOfPi/6AfRDHwKSJ\n4mU7Pv4iGdgWcQxMhhDmPHaVgW0Rx8DQBDHrcr4ysA5xDAxFDLMLp8eVWAbOI46BIYhi9sm5ysB5\nxDFwMIKYQ3PqBXCaOAb2ThQzIrvJQCKOgT0Qw0yJ85Nh2cQxsDOimDmwowzLIo6BrRLEzJXzk2EZ\nvu7QdwCYD2HMUhjrMF92joGNiASWyk4yzJM4BlYihuFreREfzIfTKoArE8ZwNZ4rMF12joELWeRh\nPU67gGkSx8CZRDFsj8vBwXSIY+AGohh2RyTD+MQxkEQUwz6JZBiXOIaFE8VwOCIZxiOOYaFEMYxD\nJMM4XMoNFkgYw5g8N+Hw7BzDglh4YXx2keGwxDEsgCiG6RHJcBjiGGZMFMP0iWTYL3EMMySKYX5E\nMuyHF+TBDFk8YZ48t2H37BzDTJzcLa6qVJUdZJiR4zA+/VwHtsvOMcyACIZl8tyH7RPHMGGttTMX\nR+cmwryctWt87Lx5AFiP0ypggiyEwGl+KIbtsHMME3PVMLZQwjxctGt8Fj88w2bEMUzEOr86Fcgw\nbauG8TGnWsD6nFYBg7PAAevyAzKszs4xDGpbOz8WR5imdXeNz2InGa5OHAMAQCeOYTC72OGxewzT\nss1d45PsIMPlxDEMYteLlkCGadhVGJ8kkuF84hgGYJECDsHcA19LHMOBWZyAQzIHwY1cyg0O5BAL\n0slTKyyIMJ5DPTeddgUvsnMMByBMgRGZm0Acw16N9CIYO0QwllGekyPNU3AI4hgWyK9QYSz7uEIF\ncDXOOYY9sOABU+OHaJbKzjHs2KhhPOr9gqUa9Tk56v2CXRHHsEOjLyqttVSVnSE4kOPn3xTmClgK\np1XADlhIgLlxmgVLYecYtmxqYTy1+wtzM7Xn4NTuL6xq7Tiuqger6rNV9Zmq+nBVfUNVvbKqHq+q\nz1fVx6vqllNf/4Wqeqqq3ryduw9jmfKiYTcI9mvKz7kpz3VwmVpngFfVa5P8XpLXtdb+X1X9zyQf\nS/JtSf6mtfZzVfVTSb65tfaeqronyYeTfHeS1yT53SR3t9a+curv9WxjkuawULiUFOzPnJ5vU458\nlq21dubgXXfn+MtJnk/ysqq6KcnLkvxVknuTPNS/5qEkb+sfvzXJw62151trTyf5YpI3rPlvw1Dm\nsLgBrMscyNysFcettb9L8vNJ/jJHUfwPrbXHk9zaWrvev+x6klv7x69Ocu3EX3EtRzvIMGlzWhS8\n2Ab2Y067xsfm9H+BteK4qr41yU8keW2OwvflVfUjJ7+mHT1TLnq2eCYxWXN9e1WBDLs1xzA+Ntd5\nkeVZ97SK70ryh621v22tvZDkt5J8b5Jnq+q2JKmqVyV5rn/9M0nuOPH9t/fbAABgGOvG8VNJvqeq\nXlpHPwa/KcmTST6a5P7+NfcnebR//FiS+6rq5qq6M8ldSZ5Y/27DYSxhZ8TuMezGnHeNT1rCPMm8\nrfUmIK21P62qDyb54yRfSfKpJL+c5BVJHqmqB5I8neTt/eufrKpHchTQLyR5V/PMAQBgMGtdym1X\nXMqNUY30PNmXpexywT4s9fnkN1CMbNuXcgMAuNDSfhhgHsQxXGDJ58459xi2Y6m7xseWPI8yTeIY\nzmEyF8iwqaWH8UkeA6ZCHMMZTOIA22duZQrEMZxi8r6R3WNYj13js3k8GJ04hhNM2mcTyLAaYXwx\njwsjE8fQmawB9secy6jEMcQkfRUeI1iN58zlPEaMSByzeCbn1Ti1Ai7mObIaczCjEccsmkl5Nc49\nhos513g9Hi9GIo5ZLJMxwDjMyYxCHLNIJuH12T2Gs9k13pzHjhGIYxbH5Ls5gQw3Esbb4zHk0MQx\ni2LSBRifuZpDEscshsl2u+wewxG7xrvh8eRQxDEAAHTimNlrrdmB2BGPKxzxXNgN8zeHII6BjTm1\ngqUy9mF+xDGzZsdh95x7zFI513h/PMbskzhmtkymAPNhTmdfxDGzZBLdL7vHLI1d48PweLMP4hjY\nCoHMUghjmLebDn0HYJssVgDz5gdxds3OMQAAdHaOmQU7xmM4uaPjmDBHxvY47CCzK3aOgZ2wYDE3\nxjQsgzhm8uziACyXNYBtE8dMmklxTH7dydy4QsXYHBe2SRwDOyGQmQthDMsijpksCxUAx6wJbIs4\nZpJMgtNg95ips2s8LY4T2yCOmRyT37QIZKZKGE+T48WmxDGTYtID4DLWCjYhjpkMk910OXZMlbE7\nXY4d6xLHwN44tYKpMFZhucQxk2AHYPqce8xUONd4PhxD1iGOGZ7JDYB1WUNYlThmaCa1ebF7zOjs\nGs+T48kqxDGwVwKZUQljIBHHDMwCBcC2WFO4KnEMAACdOGZIfsKft9ZaWmtOrWAYVfXVccl8Ob5c\nhThmOCavZRHIHJoxuCzWGC4jjgEAoBPHDMVP9MviyhUcmitULJPjzUXEMcMwWS2T486hGYPL5Lhz\nHnEMAACdOGYIfoJftuMrVzi9gn05Hm/mnmVz/DmLOObgTE4AHIo1iNPEMTAECxT7ZswBZxHHHIwL\n7nOaNwZhH5xOwWnWI04Sx8BwBDK7YmwBlxHHAADQiWMOwq+vOI83BmFXvOEHlzE2SMQxB2Dy4TIC\nmW0TxlyVMYI4BgCAThyzV34i56rsHrMtdo1ZlbGybOKYvTHZADAV1qzlEsfAsOwesym7xsCqxDEA\nAHTimJ3zzkNswu4x67JrzKasX8skjoHhCWRWJYyBdYljAADoxDE7ZdeGbTGWWJUxw7YYS8sijtkZ\nkwnb1lpLVTm9gnMdjw/zD9tmTC2HOAYAgE4cA5Ni94bLGCPAJsQxO2FxYtecWsFpxgS7Zm1bBnEM\nTI5Lu3GaS7cB2yKO2TqLEwBzZY2bP3HMVpk02Be7xxyza8y+GWvzJo6ByRLICGNg28QxAAB04pit\nsXPDIdg9Xi67xhyScTdf4pitMEkAsDTWvnkSx8DkWaCWy7EHtk0cA7Ph1IrlcKyBXRHHbMzODSNw\n7vFyONeYkRiH8yOOAQCgE8dsxE/MjMTu8fzZNWZExuO8iGNgVixS8+cYA7skjlmbBYqR2T2eH8eU\nkVkT50McAwBAJ45Zi5+QGVlrLa01O40zUlVfPa4wKuNzHsQxMGsCefocQ2Cfbjr0HQDgYpvGod0s\ngKurkSbNqhrnznCukcYMXGZKl/7a9w7plB6TKdxXOOa3HdPQWjvzQIljVjLSeIGrGimwprJojvRY\njXBfYFVTea4v2Xlx7LQKgB2Z8uJ41n0XqcASiGNg9o6vXHF8xYNdmnIQX+bk/21fj6MgB/bNaRVc\n2UhjBdaxi+CacwyvYhePqTmHqTM/jM1pFQBbZNG7kaAF5kIccyUWPOZg0zcGEcSX29apF+Yc5sCb\nEU2TNwEBFmedxcoCtzqPMzBFdo6BRVnlxXlCbXOr7CQ7NQMYgZ1jgDMI4+3zmAJT4GoVXGqkMQLb\nctYupXjbr7Mee/MNc2RuGdN5V6uwc8yFLFTM1emxbfHav9OPufmGuTK2p8U5x8CiieLD8vgDo7Fz\nDAAAnZ1jYFHsVI5rn29PDXAeO8ecy+LE3Ajj6XCsmBtr6nTYOQZmT2hNkytYAIdg55gzWYyYC2E8\nfY4hc2FtnQZxDMyWqJoPxxLYF6dVALMjpObJaRbAPtg5BmZFGM+fYwzskjjma9iVYYqqSjQtiOPN\nVFljxyeOAQCgE8fcwE+0TI0dxGVz/Jkia+3YxDEAAHSuVgFMkt1CTnIlC2Bb7BwDkyOMOY+xAWxK\nHPNVdlyYAvHDZYwRpsCaOy5xDEyG6OGqjBVgXeIYmASxw6qMGWAd4hgYnshhXcYOsCpxDAxN3LAp\nYwhYhTgmiRcGAMC+WXvH5DrHwJDs9rFNroMMXJWdYywWAHAg1uDxiGNgOHaN2RVjC7jM2nFcVbdU\n1Ueq6nNV9WRVvbGqXllVj1fV56vq41V1y4mvf7CqvlBVT1XVm7dz94G5ES/smjEGXGSTneP3JvlY\na+11Sb49yVNJ3pPk8dba3Uk+0f+cqronyTuS3JPkLUneV1V2rYEbiBb2xVgDzrNWoFbVNyX5/tba\n+5OktfZCa+0fk9yb5KH+ZQ8leVv/+K1JHm6tPd9aezrJF5O8YZM7DsyLWGHfjDngLOvu3t6Z5K+r\n6gNV9amq+pWq+sYkt7bWrvevuZ7k1v7xq5NcO/H915K8Zs1/my3yQgAAOCxr8VjWjeObkrw+yfta\na69P8i/pp1Aca0dH+qKjbSQASezgcTjGHnDaunF8Lcm11ton+58/kqNYfraqbkuSqnpVkuf6559J\ncseJ77+93wYAAMNYK45ba88m+VJV3d1velOSzyb5aJL7+233J3m0f/xYkvuq6uaqujPJXUmeWPte\nsxV+jcOhVZWdOw7OOGQE1uRxbPIOeT+W5Deq6uYkf5HkR5O8JMkjVfVAkqeTvD1JWmtPVtUjSZ5M\n8kKSdzWjAACAwdRIjVpV49yZBRjp2LM8duoYkXmRQzIv7ldr7cwH3LWGAQCgE8fA3tkdYVTGJiCO\nAQCgE8cL5bw6DsXOHKMzRjkUa/MYxDGwN6KDqTBWYbnEMQAAdOIYAAA6cbxAzmniEPyamqkxZjkE\na/ThiWNg50QGU2XswvKIYwAA6MQxAAB04hjYKb+WZuqMYVgWcQzsjKhgLoxlWA5xDAAAnTheGJeI\nYV/stDE3xjT7Yq0+LHEMAACdOAa2zg4bc2Vsw/yJYwAA6MQxAAB04nhBnOAPANNgzT4ccQxslXMy\nmTtjHOZNHAMAQCeOga2xo8ZSGOswX+IYAAA6cQwAAJ04BgCAThwDW+EcTJbGmId5EscAANCJYwAA\n6MQxsDG/XmapjH2YH3G8EN6GEgCmxdp9GOIYAAA6cQwAAJ04BgCAThwDAEAnjoGNeLU+S+c5APMi\njgEAoBPHAADQiWMAAOjEMQAAdOJ4AbzDDgBMkzV8/8QxAAB04hgAADpxDAAAnTgGAIBOHANr885g\ncMRzAeZDHAMAQCeOAQCgE8cAANCJYwAA6MQxAAB04hgAADpxDAAAnTgGAIBOHAMAQCeOAQCgE8cA\nANCJYwAA6MQxAAB04hgAADpxDKyttXbouwBD8FyA+RDHAADQiWMAAOjEMQAAdOIYAAA6cbwAVXXo\nuwAArMEavn/iGAAAOnEMAACdOAYAgE4cAwBAJ46BjXhnMJbOcwDmRRwDAEAnjgEAoBPHAADQiWMA\nAOjE8UJ4hx0AmBZr92GIY2BjXq3PUhn7MD/iGAAAOnEMAACdOAa2wq+XWRpjHuZJHAMAQCeOAQCg\nE8cAANCJY2BrnIPJUhjrMF/iGAAAOnEMbJUdNebOGId5E8cL4m0oAWAarNmHI44BAKATxwAA0Ilj\nYOuck8lcGdswf+IYAAA6cbwwTvBnX+ywMTfGNPtirT4scQwAAJ04BnbGThtzYSzDcohjYKdEBVNn\nDMOyiGMAAOjEMQAAdOIY2Dm/lmaqjF1YHnG8QC4RwyGIDKbGmOUQrNGHJ44BAKATxwAA0IljYG/8\nmpqpMFZhucTxQjmniUMRHYzOGOVQrM1jEMcAANCJY2Dv7MwxKmMTEMcAANCJY+Ag7NAxGmMSSMQx\nAAB8lTheMK+K5dBaa3brODjjkBFYk8chjgEAoBPHwMHZteNQjD3gNHG8cH6NAwCHZS0eizgGhmAH\nj30z5oCziGNgGGKFfTHWgPOIY2AoooVdM8aAi4hjYDjihV0xtoDLiGO8EAAADsQaPJ6bDn0HAM5y\nvMNn4WAb7BgDV2XnmCQCBAD2zdo7JnEMDM2OH5syhoBViGNgeOKGdRk7wKrEMTAJIodVGTPAOjaK\n46p6SVV9uqo+2v/8yqp6vKo+X1Ufr6pbTnztg1X1hap6qqrevOkdB5ZH7HBVxgqwrk13jt+d5Mkk\nx7PQe5I83lq7O8kn+p9TVfckeUeSe5K8Jcn7qsqu9WC8MIApED1cxhhhCqy541o7UKvq9iQ/mORX\nkxwf4XuTPNQ/fijJ2/rHb03ycGvt+dba00m+mOQN6/7bwLKJH85jbACb2mT39heS/GSSr5y47dbW\n2vX+8fUkt/aPX53k2omvu5bkNRv828DCtdaEEF9lPADbslYcV9UPJXmutfbpvLhrfIN2NEtdNFOZ\nxQAAGMq675D3fUnuraofTPKvkvzrqvpQkutVdVtr7dmqelWS5/rXP5PkjhPff3u/jcFUld0XJsU7\n6S2b+YopMl+Nba2d49baT7fW7mit3ZnkviS/11p7Z5LHktzfv+z+JI/2jx9Lcl9V3VxVdya5K8kT\nm911AADYrnV3jk87/tH9Z5M8UlUPJHk6yduTpLX2ZFU9kqMrW7yQ5F3Nj/vDsnvMFNlBXhZzFFNl\njhpfjTTBVNU4d2bhRhoXsCqLz7yZn5gy89M4WmtnHoxt7RwDDMMu8jyJYmAfvBEHMFtiaj4cS2Bf\nxDFnsuPGXIiq6XMMmQtr6zQ4rQKYPadZTJMoBg7BzjHnEhLMjdiaDseKubGmToedY2BRTkaXxWos\nghgYgZ1jAADoxDGwaK01O5YH5PEHRiOOuZBfOzNXp8e2QNu/04+5+Ya5MranxTvkcSUjjRPY1PFC\nddG4tpjtxlUec/MNc2IuGdd575Bn5xjgDAJt+zymwBS4WgWwKKvsTrqyxeZWCeLWWqoqVSWkgYOx\ncwwszjrhJdZW53EGpsjOMVdiJ4c52HT3107y5bY1T5hzmAPzxDSJY4A1eEvqGwlZYC5crYKVjDRe\nYBX7uhLCUmJ5X4+jOYepWspcMGXnXa3CzjEwe/sMrTmferHPUPXiPOBQxDHAjpwVdVMJZkEKLJXT\nKljZSGMGLjOlX8/vO5yn9JhM4b7Csan8ELx0551WIY5Zy0jjBs4zl7Da5lU2pmgux5FlEMbT4Zxj\ngIkShQD7401AgFkTltPnGAL7JI5Zi18bMTJXOZifk1evgFEZn/MgjgEAoBPHrM1PyIzMrvH8OKaM\nzJo4H+IYmBUL1Pw5xsAuiWM2YpFiJC75NX/Hx9bcw0iMx3kRxwAA0IljNuYnZkZg13g57B4zEuNw\nfsQxMBvCeDkca2BXxDEweXZulsuxB7ZNHLMVFigAlsbaN0/imK0xSXAIzjVeLucec0jG3XyJYwAA\n6MQxMFl2jbF7DGybOGarLFDsizDmmEBm34y1eRPHbJ1JA4C5ssbNnzgGJseuMafZPQa2RRyzExYo\ndk0Yc5oxwa5Z25ZBHAOTYnHiMsYIsAlxDAAAnThmZ+zesG1VldaaX59zruPxYf5h24yp5RDH7JTJ\nhG0xlliVMcO2GEvLIo4BAKATx8DwXLqNVbm0G7AucczOVZUFirUJY9YlkNmU9WuZxDEAAHTiGBiW\nXWM2ZfcYWJU4Zm8sTgBMhTVrucQxe2Wy4arsGrMtdo9ZlbGybOIYAAA6ccze+Ymcy9g1ZtvsHnNV\nxgjimIMw+XAeYcyuCGQuY2yQiGMAAPgqcQwMx64xu2JsAZcRxxyMdx7itKoSL+xca83cww2sR5wk\njoEhWJjYN2MOOIs45uAsUAAcijWI08QxQzA5Ldvx6RROqWBfjsebuWfZHH/OIo4BAKATxwzDT/DL\n5LhzaMbgMjnunEccMxST1bJ4ww8OzRuDLJPjzUXEMQAAdOKY4fiJflnsGnNoxuCyWGO4jDhmSCav\neTu+4L4oYRTHV64w98yb48tViGMAAOjEMcPyEz4A22JN4arEMbBXrlDBqFy5AkjEMYOzSM2LMGZ0\nAnmeHE9WIY4ZnkkNgHVZQ1iVOGYSTG7TZ9eYqbB7PB+OIesQx8DeCGOmwliF5RLHTIYdgOly7Jgq\nY3e6HDvWJY6ZFJMdAJexVrAJcczkmPSmxbnGTJVzj6fJ8WJT4phJMvlNgzBm6gTytDhObIM4ZrJM\nggAcsyawLeIY2Am7xsyF3WNYFnHMpFmsxiSMmRuBPDbHhW0Sx0yeSRFguawBbJs4BnbCrjFzY0zD\nMtx06DsA2+DX+GNwHJi71ppxPgg7xuyKnWMAAOjsHDMrdnQA5s2OMbtm5xjYCj+YsBSuXAHzJo6Z\nJYvWfgljlkYgH4bHm30Qx8yWSRRgPszp7Is4ZtZMprtn15ilsnu8Px5j9kkcAxsTxiyVsQ/zI46Z\nvaqy67AjHlc44rmwG+ZvDkEcAwBAJ45ZDLsP2+VcYzji3OPd8HhyKOKYRTHZAozPXM0hiWMWx6S7\nObvGcCO7x9vjMeTQxDGLZPJdnzCGswnkzXnsGIE4ZrFMwgDjMCczCnHMopmMV2PXGC5m93g9Hi9G\nIo5ZPJPyaoQxXMxzZDXmYEYjjiEm56vwGMFqPGcu5zFiROIYOpM0wP6YcxmVOIYTTNZnc64xrMa5\nxxfzuDAycQynmLRvJIxhPQL5bB4PRieO4Qwmb4DtM7cyBeIYzmESt2sMm7J7/CKPAVMhjuECVbXY\nCV0Yw3YsPZCXPI8yTeIYANgJUcwUiWO4gqXtfNg1hu1a4u7xkv6vzIs4BgCAThzDCpawg2zXGHZj\nKbvHS5gnmTdxDAAAnTiGNcx1Z8SuMezWnHeP5zovsjziGDYwp4VAGMN+zDGQ5/R/AXEMG7IoAEtm\nDmRuxDFswdR/nWjXGPZrDrvHU5/34DziGLZoyguFMIb9mvJzbspzHVxGHMOWTW3RmNr9hbmZ2nNw\navcXVnXToe8AzJHTFIC5EcUshZ1j2KHRF5OqSmtNxMOBHD//pjBXwFKIY9ixUReVUe8XLNWoz8lR\n7xfsitMqYA+cZgFMjShmqewcwwKJdRjLHC7tBnMhjmGPRrouqDCGsYzynBxpnoJDWCuOq+qOqvr9\nqvpsVf15Vf14v/2VVfV4VX2+qj5eVbec+J4Hq+oLVfVUVb15W/8BmCILDzAicxMktc5PqlV1W5Lb\nWmt/UlUvT/K/k7wtyY8m+ZvW2s9V1U8l+ebW2nuq6p4kH07y3Ulek+R3k9zdWvvKqb93jB+bYY/2\nuVvkdAoY3yGep6KYJWqtnTnw19o5bq0921r7k/7xPyf5XI6i994kD/UveyhHwZwkb03ycGvt+dba\n00m+mOQN6/zbMDcWJeCQzEFwo43POa6q1yb5jiR/lOTW1tr1/qnrSW7tH786ybUT33YtRzENxOIE\nHIa5B77WRnHcT6n4zSTvbq3908nPtaPfB130OyG/14UTdv0iGKdUwDTs48oVXnQH51s7jqvq63MU\nxh9qrT3ab77ez0dOVb0qyXP99meS3HHi22/vtwGn7GLREsYwLbsKZFEMl1v3ahWV5NeSPNla+8UT\nn3osyf394/uTPHri9vuq6uaqujPJXUmeWO8uAwDAbqx7tYp/n+R/JfmzvHh6xIM5Ct5HkvybJE8n\neXtr7R/69/x0kv+U5IUcnYbxO2f8vba14JRNd3vtGsN0bev5a7cYvtZ5V6tYK453RRzD+dZ9ropj\nmK5Nn7+iGM631Uu5Afu3zrmCwhimbd1zj51bDOsTxzAxV13whDHMw6qBLIphMzcd+g4AqxO+wGmi\nGLbDzjFM2Hm/OhXPMC8X7R47hQK2SxzDDFgYYZk892H7nFYBM3F6kbRrDPPSWrNLDHtg5xhmSBjD\nPHluw+7zcdpLAAAHyElEQVTZOYYZcs4xzI8dY9gPcQwzJpJh+kQx7Jc4hgUQyTA9ohgOQxzDgohk\nGJ8ohsPygjxYIIsvjMlzEw7PzjEslF1kGIcohnGIY1g4kQyHI4phPOIYSCKSYZ9EMYxLHAM3EMmw\nO6IYxieOgTOJZNgeUQzTIY6BC51c1IUyXJ0ghmlyKTfgyiz2cDWeKzBddo6BlZxe9O0mgxiGObFz\nDGxEFLB0ngMwL3aOgY05L5mlEcQwX+IY2CqhzFwJYlgGcQzsjMvBMQeiGJZFHAM750V8TIkYhmUT\nx8De2VFmRKIYSMQxcEDOT+bQBDFwmjgGhmA3mX0SxcB5xDEwFOcnswtiGLgqcQwMzakXrEsQA+sQ\nx8BknBU7gplECAPbI46BSXOu8rKJYmDbxDEwC3aV508IA/sgjoHZcr7y9AliYN/EMbAI50WWaB6D\nCAZGIY6BRbO7fDiCGBiROAboLos18bwa8QtMkTgGuCIv+jufEAbmQhwDbGCVKJxaSAteYInEMcCe\nXDU2dx3RohfgfOIYYDDiFeBwvu7QdwAAAEYhjgEAoBPHAADQiWMAAOjEMQAAdOIYAAA6cQwAAJ04\nBgCAThwDAEAnjgEAoBPHAADQiWMAAOjEMQAAdOIYAAA6cQwAAJ04BgCAThwDAEAnjgEAoBPHAADQ\niWMAAOjEMQAAdOIYAAA6cQwAAJ04BgCAThwDAEAnjgEAoBPHAADQiWMAAOjEMQAAdOIYAAA6cQwA\nAJ04BgCAThwDAEAnjgEAoBPHAADQiWMAAOjEMQAAdOIYAAA6cQwAAJ04BgCAThwDAEAnjgEAoBPH\nAADQiWMAAOjEMQAAdOIYAAA6cQwAAJ04BgCAThwDAEAnjgEAoBPHAADQiWMAAOjEMQAAdOIYAAA6\ncQwAAJ04BgCAThwDAEAnjgEAoBPHAADQiWMAAOjEMQAAdOIYAAA6cQwAAJ04BgCAThwDAEAnjgEA\noBPHAADQiWMAAOjEMQAAdOIYAAA6cQwAAJ04BgCAThwDAEAnjgEAoBPHAADQiWMAAOjEMQAAdOIY\nAAA6cQwAAJ04BgCAbq9xXFVvqaqnquoLVfVT+/y3AQDgMtVa288/VPWSJP8nyZuSPJPkk0l+uLX2\nuRNfs587AwDAorXW6qzb97lz/IYkX2ytPd1aez7J/0jy1j3++wAAcKF9xvFrknzpxJ+v9dsAAGAI\n+4xjp0wAADC0fcbxM0nuOPHnO3K0ewwAAEPY5wvybsrRC/L+Y5K/SvJETr0gDwAADummff1DrbUX\nquq/JvmdJC9J8mvCGACAkext5xgAAEbnHfIAAKAbJo69e95yVNUdVfX7VfXZqvrzqvrxfvsrq+rx\nqvp8VX28qm458T0P9rHxVFW9+XD3nl2oqpdU1aer6qP9z8bCAlXVLVX1kar6XFU9WVVvNBaWpx/X\nz1bVZ6rqw1X1DcbBMlTV+6vqelV95sRtKx/7qvrOPn6+UFXvXfV+DBHH/d3z/nuStyS5J8kPV9Xr\nDnuv2KHnk/y31tq3JfmeJP+lH+/3JHm8tXZ3kk/0P6eq7knyjhyNjbckeV9VDTF22Zp3J3kyL17y\n0VhYpvcm+Vhr7XVJvj3JUzEWFqWqXpvkPyd5fWvt3+boNUr3xThYig/k6DietMqxP37Hu19K8kBr\n7a4kd1XV6b/zQqMMIO+etyCttWdba3/SP/7nJJ/L0RvC3Jvkof5lDyV5W//4rUkebq0931p7OskX\nczRmmIGquj3JDyb51STHE5uxsDBV9U1Jvr+19v7k6EXcrbV/jLGwNF/O0QbKy/pVrl6WoytcGQcL\n0Fr7gyR/f+rmVY79G6vqVUle0Vp7on/dB098z5WMEsfePW+h+i7BdyT5oyS3ttau909dT3Jr//jV\nufGa2MbHvPxCkp9M8pUTtxkLy3Nnkr+uqg9U1aeq6leq6htjLCxKa+3vkvx8kr/MURT/Q2vt8RgH\nS7bqsT99+zNZcUyMEscumbFAVfXyJL+Z5N2ttX86+bl2dBmVi8aFMTMDVfVDSZ5rrX06L+4a38BY\nWIybkrw+yftaa69P8i/pvz49ZizMX1V9a5KfSPLaHEXOy6vqR05+jXGwXFc49lsxShx797yFqaqv\nz1EYf6i19mi/+XpV3dY//6okz/XbT4+P2/ttTN/3Jbm3qv5vkoeT/Ieq+lCMhSW6luRaa+2T/c8f\nyVEsP2ssLMp3JfnD1trfttZeSPJbSb43xsGSrbIeXOu3337q9pXGxChx/Mc5OmH6tVV1c45OsH7s\nwPeJHeknzP9akidba7944lOPJbm/f3x/kkdP3H5fVd1cVXcmuStH77DIxLXWfrq1dkdr7c4cvejm\n91pr74yxsDittWeTfKmq7u43vSnJZ5N8NMbCkjyV5Huq6qV9rXhTjl6saxws10rrQZ9LvtyvdlNJ\n3nnie65kb++QdxHvnrc4/y7JjyT5s6r6dL/twSQ/m+SRqnogydNJ3p4krbUnq+qRHE2QLyR5V/Pu\nNXN1fFyNhWX6sSS/0TdJ/iLJj+ZoTTAWFqK19qdV9cEcbZp9JcmnkvxyklfEOJi9qno4yQ8k+Zaq\n+lKSn8l668G7kvx6kpfm6Ao4v73S/TCGAADgyCinVQAAwMGJYwAA6MQxAAB04hgAADpxDAAAnTgG\nAIBOHAMAQPf/AQDQpQFrj2eSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x145f0c690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12,12))\n",
    "plt.imshow(Lyot_stop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bin to final array resolution, crop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "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": 346,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binning the original aperture array 2x\n"
     ]
    }
   ],
   "source": [
    "N_orig = D/2\n",
    "scalefac = int(N_orig/N)\n",
    "print(\"Binning the original aperture array {0:d}x\".format(scalefac))\n",
    "Lyot_stop_binned = np.reshape(Lyot_stop, (Lyot_stop.shape[0]/scalefac, scalefac, \n",
    "                                          Lyot_stop.shape[1]/scalefac, scalefac)).mean(1).mean(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(250, 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": 348,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1463a3a90>"
      ]
     },
     "execution_count": 348,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAJKCAYAAADnWquFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF8FJREFUeJzt3U/IbHd5B/DvU6OL2oIVy03UWAM10JSCWXgpWPFCS9CN\nsRurpRBaEcG2iquaLjS7SsHSRbGbxpKWNm1ADHFRahSlrrQW/0RjqgEveNN4Y8FS3SX118U7Sea+\nef/dZ2beOWfm84GQuWfmvu9Jzj3v/c7398w5NcYIAADX72e2vQMAAHMlSAEANAlSAABNghQAQJMg\nBQDQJEgBADTdsI1vWlWuuQAAzMYYo47arpECAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCA\nJkEKAKBJkAIAaDoxSFXVzVX1+ar6VlV9s6rev9j+8qp6uKq+U1WfqaqXLf2eu6vqu1X1WFXdsen/\nAACAbakxxvFPVt2Y5MYxxteq6ueS/EeStyf5/ST/Pcb486r6kyS/MMb4UFXdluQfk7whyauSfDbJ\nrWOMnx76usd/UwCAiRlj1FHbT2ykxhg/GGN8bfH4J0m+nYOA9LYk9y1edl8OwlWS3Jnk/jHG02OM\ny0keT3Jx5b0HAJigM89IVdVrk9ye5EtJLowxri6euprkwuLxK5NcWfptV3IQvAAAds6ZgtRiWe+T\nST4wxvjx8nPjYG3wpKU6y3gAwE46NUhV1YtzEKL+fozx4GLz1cX8VKrqpiRPLbY/keTmpd/+6sU2\nAICdc9qn9irJvUkeHWP85dJTDyW5a/H4riQPLm1/Z1W9pKpuSfK6JF9e7y4DAEzDaZ/a+40k/5bk\nG3l+ie7uHISjB5K8JsnlJO8YY/zP4vf8aZI/SPJMDpYC//WIr2u5DwCYjeM+tXdikNoUQQoAmJPW\n5Q8AADieIAUA0CRIAQA0CVIAAE2CFABAkyAFANAkSAEANAlSAABNghQAQJMgBQDQJEgBADQJUgAA\nTYIUAECTIAUA0CRIAQA0CVIAAE2CFABAkyAFANB0w7a+8Rhj7V+zqtb+NQEAjqORAgBo2lojtQmb\naLmOo/0CADRSAABNghQAQNNOLe2dp1WXES0NAsD8aaQAAJo0UlvSabS0WAAwLRopAIAmjdSMXE+L\npb0CgM3TSAEANGmkdtT1zmBpsADg+mmkAACaBCkAgCZLeyQ521Kg5T8AuJZGCgCgSSPFmR3XWmmq\nANhXGikAgCaNFCszXwXAvtJIAQA0aaQ4Fye1VtoqAOZKIwUA0CRIAQA0Wdpj6yz7ATBXGikAgCaN\nFJPmIqAATJlGCgCgSSPFLB3VVGmpADhvGikAgCaNFDvDp/8AOG8aKQCAJo0Ue+FwW6WhAmAdNFIA\nAE2CFABAk6U99pLBdADWQSMFANCkkYJDDKYDcFYaKQCAJo0UnMLtaAA4jkYKAKBJIwUNx33qT1MF\nsF80UgAATYIUAECTpT1YI4PpAPtFIwUA0KSRgg1zgU+A3aWRAgBo0kjBOTNHBbA7NFIAAE2CFEzA\nGOPYi3wCMF2CFABAkyAFANBk2BwmxKUSAOZFIwUA0KSRgglzqQSAadNIAQA0CVIwMy6VADAdghQA\nQJMZKZg5c1QA26ORAgBoEqQAAJos7cHMnGXZ7tnlPkt8AJulkQIAaNJIwUwcbpfOcgkEg+gAm6WR\nAgBo0kjBzKx6MU43RgZYH40UAECTRgom7DzaIp/wA+jTSAEANGmkYCY2faNin/ADuH4aKQCAJkEK\nAKDJ0h5M0FSW1AyiA5xMIwUA0KSRggnb9ID5WbmIJ8DRNFIAAE0aKZiQuTQ9yw3VXPYZYBM0UgAA\nTRopmIDDrc5UZqPOwif7gH2mkQIAaBKkAACaLO3BhMxpSe8wl0gA9pFGCgCgSSMFW7TLrY1LJAD7\n4NRGqqo+UVVXq+qRpW33VNWVqvrq4p+3Lj13d1V9t6oeq6o7NrXjAADbVqfNZFTVm5L8JMnfjTF+\nbbHtI0l+PMb4i0OvvS3JPyZ5Q5JXJflsklvHGD899Lox51kQWJdnm5p9Ox80VMDcjDGO/MF1aiM1\nxvhikh8d8dRRX/DOJPePMZ4eY1xO8niSi9exnwAAs7HKsPkfV9XXq+reqnrZYtsrk1xZes2VHDRT\nwEJVPffPvhpj7F0LB+ymbpD66yS3JHl9kieTfOyE1/ppCQDspNan9sYYTz37uKr+JsmnF798IsnN\nSy999WLbC9xzzz3PPb506VIuXbrU2RUAgK05ddg8SarqtUk+vTRsftMY48nF4w8mecMY43eXhs0v\n5vlh818+PFlu2Jx9tryk5zx43j4vdQLTd9yw+amNVFXdn+TNSV5RVd9P8pEkl6rq9TlYtvtekvcu\nvsmjVfVAkkeTPJPkfRITALCrztRIrf2baqTYQ0c1Ls6DF9JMAVPUvvwBAABHc4sY2LDDDYsW6mRu\nfgzMiUYKAKBJIwXnRBPV4+bHwJRppAAAmgQpAIAmS3uwIZah1u/ZZT7/b4Gp0EgBADRppGDDDJmv\nn2YKmAqNFABAk0YK1khDcr40U8C2aaQAAJo0UrAGbki8XUf9v9ZSAedBIwUA0CRIATtpjKEVBDZO\nkAIAaBKkAACaDJvDCg4PNFtKmh6XSAA2SSMFANCkkYI10ERNn2YK2ASNFABAk0YKGrQa87XcHjqO\nwKo0UgAATRopWIHZqHkzNwWsSiMFANAkSAEANFnagzOy/LO7LPEBXRopAIAmjRSc4qiWwpD5btJM\nAddLIwUA0KSRgjPSQu0PF+0EzkojBQDQpJGCY2giSMxNASfTSAEANAlSAABNlvbgFIbMSSzxAUfT\nSAEANGmk4BCNAyfRTAHLNFIAAE0aKTiG2ShO4qKdQKKRAgBo00hBNAqsxtwU7C+NFABAkyAFANBk\naY+9dtRSjCFzuizxwf7RSAEANGmkIFoo1kszBftDIwUA0KSRYi9pCjgPLtoJu08jBQDQpJFir5mN\n4ryYm4LdpJECAGgSpAAAmiztsVcsqwCwThopAIAmjRR74XATZcicbTF0DrtFIwUA0CRIAWzBGEMz\nCjtAkAIAaDIjxc46agZFA8DUmJmCedNIAQA0aaTYC5oops4NjmGeNFIAAE2CFABAk6U9do5lEebO\nADrMh0YKAKBJI8XOMmAOwKZppAAAmjRS7AzzJOwas1IwfRopAIAmjRSzd/jdutkodo1mCqZLIwUA\n0CRIAQA0WdpjZ1jSY9e5Hx9Mj0YKAKBJI8UseTfOvjOADtOgkQIAaNJIMXtmowDYFo0UAECTRopZ\nMQ8C1zIrBdulkQIAaBKkAACaLO0xSwbM4VqW+GA7NFIAAE0aKWbBu2wApkgjBQDQpJFi0g43UWaj\n4GRubAznSyMFANCkkWIWNFFw/XySDzZPIwUA0CRIAQA0WdpjcixDwHpZ4oPN0UgBADRppJg0Q+YA\nTJlGCgCgSSPFZJjfgM0yKwXrp5ECAGjSSDE55qJgs9xGBtZHIwUA0KSRYuu8IwZgrjRSAABNghQA\nQJOlPbbm8JKeIXM4fy6JAKs5tZGqqk9U1dWqemRp28ur6uGq+k5VfaaqXrb03N1V9d2qeqyq7tjU\njgMAbNtZlvb+NslbDm37UJKHxxi3Jvnc4tepqtuS/E6S2xa/5+NVZfmQE40xtFGwZc5D6Dk15Iwx\nvpjkR4c2vy3JfYvH9yV5++LxnUnuH2M8Pca4nOTxJBfXs6sAANPSbYsujDGuLh5fTXJh8fiVSa4s\nve5Kklc1vwc7qKqe+wcA5m7lYfMxxqiqk/rgI5+75557nnt86dKlXLp0adVdAQA4V90gdbWqbhxj\n/KCqbkry1GL7E0luXnrdqxfbXmA5SLGfzGPA9PgUH1yf7tLeQ0nuWjy+K8mDS9vfWVUvqapbkrwu\nyZdX20UAgGk6tZGqqvuTvDnJK6rq+0k+nOSjSR6oqncnuZzkHUkyxni0qh5I8miSZ5K8b6gdAIAd\nVdvIOVUlX+2Zo5YJ/BmA6bK0B9caYxx5UrjGEwBAk1vEsFFuAwPzZOgczkYjBQDQpJEC4FjLLbJ2\nCl5IIwUA0KSRYiPMRgGwDzRSAABNghQAQJOlPTbKkh7sDpdEgBfSSAEANGmkWCvvVAHYJxopAIAm\njRQbYTYKdpdZKXieRgoAoEkjxcq8KwVgX2mkAACaBCkAgCZLe7QdtaRnyBz2h6Fz0EgBALRppFiZ\nFgqAfaWRAgBo0khx3cxDAMuWW2k/H9g3GikAgCaNFG1mo4DDfJKPfaORAgBo0khxZt5hAsC1NFIA\nAE2CFABAk6U9rpshc+A0hs7ZFxopAIAmjRQn8m4SAI6nkQIAaNJIcaSjmiizUcD1MivFrtNIAQA0\naaQ4kRYKAI6nkQIAaBKkAACaLO1xDQOhwCYsjwn4OcMu0UgBADRppDiSIXMAOJ1GCgCgSSNFEjML\nANChkQIAaNJIcQ2zUcCmuW0Mu0QjBQDQJEgBADRZ2ttjanUAWI1GCgCgSSO1h45qogyZA+fN0Dm7\nQCMFANCkkdpjWigAWI1GCgCgaWuN1CbWxDUsJzOHAEyRWSnmTCMFANAkSAEANO3UsHm3Ft63JcF9\n++8FgE3RSAEANO1UI9V1liZrzi2OAU4A2AyNFABAk0bqjI5rdabaVLkNDDA3yz+jNOnMhUYKAKBJ\nI7Wik941aYAAYLdppAAAmjRSG7SNtspsFLAL3DaGudBIAQA0CVIAAE2W9rbkcF297uU3y3kAsHka\nKQCAJo3URKw6JG4gEwDOn0YKAKBJIzVhnTkqs1HALnEZBKZOIwUA0KSRmhHvyABgWjRSAABNghQA\nQJOlvR1hyBzYZYbOmSqNFABAk0ZqR2z6ljMAwAtppAAAmjRSO2rVW84AAKfTSAEANAlSe6SqfOIF\nmLUxhnadSRGkAACaBCkAgCZBag9Z4gOA9RCkAACaXP5gj7lEAgCsRiMFANCkkeIaz7ZUmilgypZ/\nRpn5ZJs0UgAATYIUR/LJPgA4nSAFANAkSAEANBk250TLy3sG0AHgWhopAIAmQYozM4AOTNEYQ2PO\n1ghSAABNZqS4bodbKe8EAdhXGikAgCaNFCvzyT4A9pVGCgCgSZBirXyyD4B9stLSXlVdTvK/Sf4v\nydNjjItV9fIk/5zkl5JcTvKOMcb/rLifAACTs2ojNZJcGmPcPsa4uNj2oSQPjzFuTfK5xa8BAHbO\nOpb2Dq/jvC3JfYvH9yV5+xq+BzNjiQ84by7MyTaso5H6bFV9pares9h2YYxxdfH4apILK34PAIBJ\nWvXyB28cYzxZVb+Y5OGqemz5yTHGqCpvDwCAnbRSIzXGeHLx7x8m+VSSi0muVtWNSVJVNyV5atWd\nBACYonaQqqqfraqfXzx+aZI7kjyS5KEkdy1edleSB1fdSebLrBQAu6y6g3lVdUsOWqjkYInwH8YY\nf7a4/MEDSV6TYy5/YLlv/xgABc6LN29swhjjyD9Y7SC1CkFqfwlUwHkRqFin44KUK5sDADQJUgAA\nTate/gCuy3LVbpkPgLnTSAEANAlSbI1LIwAwd4IUAECTIMXWaaYAmCtBCgCgyaf2ANhJy58M1nqz\nKRopAIAmQQoAoEmQYjIMnQMwN4IUAECTYXMmx21kAJgLjRQAQJMgxaSZmwJgygQpAIAmQQoAoEmQ\nAgBoEqQAAJoEKWbB0DmwijGGy6mwEYIUAECTIAUA0CRIAQA0CVLMilkpAKZEkAIAaHLTYmbJjY0B\nmAKNFABAkyAFANAkSAEANAlSAABNghSz55IIwFm5VQzrJkgBADQJUgAATYIUAECTIMXOMCsFwHkT\npAAAmgQpAIAmQQoAoEmQAgBoEqTYOYbOATgvghQAQJMgBQDQJEgBADTdsO0dgE1ZnpNyk1Jg2bM/\nE8xTsiqNFABAkyAFANAkSAEANAlSAABNghQAQJMgBQDQJEixF9w2BoBNEKQAAJoEKQCAJkEKAKBJ\nkAIAaBKkAACaBCkAgCZBir3iMggArJMgBQDQdMO2dwAAtmWM8dxjbTUdGikAgCZBCgCgSZACAGgS\npNhLPr0HwDoIUgAATYIUAECTIAUA0CRIAQA0CVIAAE2CFABAkyAFANAkSAEANAlS7DUX5gRgFYIU\nAECTIAUA0CRIAQA0CVIAAE2CFABAkyAFANAkSAEANAlSAABNghQAQJMgBQDQJEgBADQJUgAATYIU\nAEDTDdveAZiCqnru8Rhji3sCbMuz5/7yzwM4jUYKAKBJkAIAaBKkAACaBCkAgCZBCgCgSZACAGgS\npAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCAJkEKAKBJkAIAaBKkAACaNhKkquotVfVYVX23qv5k\nE98DAGDbaoyx3i9Y9aIk/5nkt5I8keTfk7xrjPHtpdes95vCGq37nADmpaq2vQtM0BjjyD8Ym2ik\nLiZ5fIxxeYzxdJJ/SnLnBr4P0PCFL3xh27vAChw/mJZNBKlXJfn+0q+vLLYBE+Av4nlz/GBaNhGk\nrIsAAHthE0HqiSQ3L/365hy0UgAAO2UTw+Y35GDY/DeT/FeSL8ewOQAwY8cNm9+wgW/0TFX9UZJ/\nTfKiJPcuh6iTdgYAYE7W3kgBAOwLVzYHAGg69yDlqufzU1WXq+obVfXVqvryYtvLq+rhqvpOVX2m\nql627f3kQFV9oqquVtUjS9uOPV5VdffifHysqu7Yzl6THHvs7qmqK4vz76tV9dal5xy7Camqm6vq\n81X1rar6ZlW9f7Hd+bfDzjVILa56/ldJ3pLktiTvqqpfOc99oGUkuTTGuH2McXGx7UNJHh5j3Jrk\nc4tfMw1/m4NzbNmRx6uqbkvyOzk4H9+S5ONVpanenqOO3UjyF4vz7/Yxxr8kjt1EPZ3kg2OMX03y\n60n+cPF3nPNvh533AXPV8/k6/AGBtyW5b/H4viRvP9/d4ThjjC8m+dGhzccdrzuT3D/GeHqMcTnJ\n4zk4T9mCY45d8sLzL3HsJmeM8YMxxtcWj3+S5Ns5uCC182+HnXeQctXzeRpJPltVX6mq9yy2XRhj\nXF08vprkwnZ2jTM67ni9Mtde5805OU1/XFVfr6p7l5aFHLsJq6rXJrk9yZfi/Ntp5x2kfERwnt44\nxrg9yVtzUFW/afnJcfDRT8d2Js5wvBzLafnrJLckeX2SJ5N87ITXOnYTUFU/l+STST4wxvjx8nPO\nv91z3kHKVc9naIzx5OLfP0zyqRxUz1er6sYkqaqbkjy1vT3kDI47XofPyVcvtjERY4ynxkKSv8nz\nSz+O3QRV1YtzEKL+fozx4GKz82+HnXeQ+kqS11XVa6vqJTkYsnvonPeB61BVP1tVP794/NIkdyR5\nJAfH7a7Fy+5K8uDRX4GJOO54PZTknVX1kqq6JcnrcnA3AiZi8Rfvs347B+df4thNTlVVknuTPDrG\n+Mulp5x/O2ztVzY/yVmues7kXEjyqYOfD7khyT+MMT5TVV9J8kBVvTvJ5STv2N4usqyq7k/y5iSv\nqKrvJ/lwko/miOM1xni0qh5I8miSZ5K8b7hK79Yccew+kuRSVb0+B0s+30vy3sSxm6g3Jvm9JN+o\nqq8utt0d599Oc2VzAIAm16sAAGgSpAAAmgQpAIAmQQoAoEmQAgBoEqQAAJoEKQCAJkEKAKDp/wGJ\nUawEj4uI+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x145f02950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(Lyot_stop_binned_crop_binary)\n",
    "#plt.imshow(np.floor(Lyot_stop_binned_quad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write the binned, cropped Lyot stop to a dat (ASCII) file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 349,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote binned, cropped LS array to /astro/opticslab1/SCDA/Apertures/InputMasks/LS/LS_quart_ann30D76_X025Pad00_N0250.dat\n"
     ]
    }
   ],
   "source": [
    "if symm is 'half':\n",
    "    LS_quart_dat_fname_tail = \"LS_half_ann{0:02d}D{1:02d}_{2:s}025Pad{3:02d}_N{4:04d}.dat\".format(\n",
    "                               int(round(100*inD)), int(round(100*outD)), secobs_key, pad, N)\n",
    "else:\n",
    "    LS_quart_dat_fname_tail = \"LS_quart_ann{0:02d}D{1:02d}_{2:s}025Pad{3:02d}_N{4:04d}.dat\".format(\n",
    "                               int(round(100*inD)), int(round(100*outD)), secobs_key, pad, N)\n",
    "    \n",
    "LS_quart_dat_fname = os.path.join(LS_dir, LS_quart_dat_fname_tail)\n",
    "#np.savetxt(LS_quart_dat_fname, Lyot_stop_binned_quad, fmt='%.6f', delimiter=\" \")\n",
    "np.savetxt(LS_quart_dat_fname, Lyot_stop_binned_crop_binary, fmt='%d', delimiter=\" \")\n",
    "print(\"Wrote binned, cropped LS array to {0:s}\".format(LS_quart_dat_fname))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 350,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['LS_quart_ann20D99_clear_N0250.dat',\n",
       " 'LS_quart_ann25D90_clear_N0250.dat',\n",
       " 'LS_quart_ann15D90_clear_N0250.dat',\n",
       " 'LS_quart_ann25D99_clear_N0250.dat',\n",
       " 'LS_quart_ann15D99_clear_N0250.dat',\n",
       " 'LS_half_ann25D81_clear_N0250.dat',\n",
       " 'LDZ_half_ann25D81_clear_Tol02_N0250.dat',\n",
       " 'LDZ_quart_ann25D81_clear_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D81_X025Pad00_N0250.dat',\n",
       " 'LS_quart_ann25D81_X025Pad05_N0250.dat',\n",
       " 'LDZ_quart_ann25D81_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D82_X025Pad00_N0250.dat',\n",
       " 'LS_quart_ann25D82_clear_N0250.dat',\n",
       " 'LS_quart_ann25D98_clear_N0250.dat',\n",
       " 'LS_quart_ann20D90_clear_N0250.dat',\n",
       " 'LS_quart_ann20D81_clear_N0250.dat',\n",
       " 'LS_quart_ann20D82_clear_N0250.dat',\n",
       " 'LS_quart_ann25D81_clear_N0250.dat',\n",
       " 'LS_quart_ann20D98_clear_N0250.dat',\n",
       " 'LS_quart_ann20D76_clear_N0250.dat',\n",
       " 'LS_quart_ann15D98_clear_N0250.dat',\n",
       " 'LS_quart_ann15D82_clear_N0250.dat',\n",
       " 'LS_quart_ann15D81_clear_N0250.dat',\n",
       " 'LS_quart_ann15D76_clear_N0250.dat',\n",
       " 'LS_quart_ann25D76_clear_N0250.dat',\n",
       " 'LDZ_quart_ann25D82_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D76_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann25D76_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D90_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann25D90_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D90_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D90_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D90_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann30D90_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D81_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann30D81_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D81_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D81_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D76_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D76_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D82_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D82_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D82_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann30D82_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D76_X025Pad00_N0250.dat']"
      ]
     },
     "execution_count": 350,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.listdir(LS_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Lyot plane dark zone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 351,
   "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": 352,
   "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": 353,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True, False],\n",
       "       [ True,  True,  True],\n",
       "       [False,  True, False]], dtype=bool)"
      ]
     },
     "execution_count": 353,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "allowed_shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 354,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 354,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])\n",
    "len(XYshifts_allowed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 355,
   "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": 356,
   "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": 357,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x14725fb50>"
      ]
     },
     "execution_count": 357,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6EAAAOZCAYAAADvcibqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3U2oLtuZF/DnSY6hFdF0aLj5FDNIBhkp0hmIIIKG0EiS\nUToO5KKZZWBPlE564NSmJzrqkR9cAh25IISIg046NqiTjkiC2rdjksEFb6BPaGjFpicJWQ723jd1\n6tT7/dZTtap+PwhnZ5997n73fmutqmf910e21gIAAAAqvGPpFwAAAMB+KEIBAAAoowgFAACgjCIU\nAACAMopQAAAAyihCAQAAKPNsiW+amc6FAQAA2LDWWk59XhIKAABAGUUoAAAAZRShAAAAlFGEAgAA\nUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEKAABA\nGUUoAAAAZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABl\nFKEAAACUUYQCAABQRhEKAABAGUUoAAAAZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRR\nhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEKAABAGUUoAAAAZRShAAAAlFGEAgAAUEYR\nCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEKAABAGUUo\nAAAAZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEA\nAACUUYQCAABQRhEKAABAGUUoAAAAZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIA\nAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEKAABAGUUoAAAAZRShAAAAlFGEAgAAUEYRCgAA\nQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEKAABAGUUoAAAA\nZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACU\nUYQCAABQRhEKAABAGUUoAAAAZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJR5tvQLAIA1\na60t/RJekJlLvwQAuIkkFAAAgDKSUAA2Y22p5Rzm+BmlqwBUkoQCAABQRhIKwCqtJdVcW0o4x+/l\nmv/m2n4vAPRDEgoAAEAZSSgAZeZON/eQzs3xM17zvlzyb/bwvgBwPkkoAAAAZSShANzFvVNO6Vmd\na37Xl7zf53yt9xtgPxShAJzlHkWmQmM7Lnkvz7l2Tn2NawdgO0zHBQAAoIwkFICIuD3plFRxyDnX\nxqnrT1IKsB2SUAAAAMpIQgF26JrUU9LEnE5dX9ckpa5ZgHWShAIAAFBGEgqwUZemnVIj1uyapPRY\nG3C9AyxHEgoAAEAZSSjABkg92bupa/pYu7CGFGA5ilCAjlxSbHqgZu+OtYFLpu9qSwD3ZTouAAAA\nZSShACsl9YT5XDJ919RdgPuShAIAAFBGEgqwEuckn9IXmM+h9nXO+lFtE+B8klAAAADKSEIBFiD1\nhH6cs37UulGA80lCAQAAKCMJBSgg+YRtGbdX60YBzqcIBbgzBSfsjym7AOczHRcAAIAyklCAG0k+\ngSmm7AJMk4QCAABQRhIKcIVj6ackA5hy6bpRfQmwVZJQAAAAykhCAY6w3hOY07F1o4f6H30O0DtF\nKMCIwhNYyrBvOdQXmbIL9M50XAAAAMpIQoHds8kQsEaH+p9jU3b1WUAPJKEAAACUkYQCuyT9BHp1\nbN2o9aJADyShAAAAlJGEArsylYBKC4BenXPEiz4OWBtJKAAAAGUkocCmOewd2JOp9aJ20AXWRhEK\nbI5NhwB+1t/ZvAhYG9NxAQAAKCMJBTbDpkMAL7N5EbA2klAAAADKSEKBLln3CXCdczYvGn8dwD1J\nQgEAACgjCQW6ciwBBeB+rBcF5iIJBQAAoIwkFFi1Q8lnZh78u9aakXuAI07NKhn2seOv1b8Ct8ol\nprZlpvl0wEGXbjpkkyKA8+ljgSqttclOwnRcAAAAypiOC6zGtUcEPH2NTYsArnOsrx3/3bCvtXkR\ncA1JKAAAAGUkocBijm06dO/vYZQe2Lt7zRYZ9qc2LwKuIQkFAACgjCQUKDdXAmptKMBl5up3zUIB\njlGEAiWu3XTo3t/fAxGwNxUDc4c2Lxp+b/0v8MR0XAAAAMpIQoFZLZGAZqYpuQBHVPTDEY5zAaZJ\nQgEAACgjCQVmsdbt+o3EA3uxhhkhx45z0Q/DfklCAQAAKJNLjJJl5vJDc8DdLb0D7pRjfdzSrw1g\nLmvu+9Z4rwDm0VqbbNySUAAAAMpYEwrcbM2j2ocOUgfYmzX3y9aJwr4oQoGLHSroent48NADbE1P\nA27HNi2a+hpgO86ajpuZ78zMb2fmv3/8/+/JzG9k5vcy8+uZ+e7B134pM7+fmd/NzE/M9cIBAADo\nz7lrQn8lIt6IiKfhqS9GxDdaax+NiG8+/v/IzI9FxC9HxMci4pMR8ZuZad0pbMihqbdrH63u4TUC\nzKGH/u/Qa+wp2QXOd7JAzMwPRsQvRcS/jIin3uFTEfHa48evRcRnHj/+dER8pbX249bamxHxg4j4\n+D1fMAAAAP06J6X85xHxTyLip4PPvdJae/748fOIeOXx4/dHxFuDr3srIj5w64sEltdae2EN5fB/\nvRv+bAA92ko/NnV/2crPBvzM0SI0M/9uRPyotfbt+FkK+oL20Csc6xn0GgAAAETE6d1x/3pEfCoz\nfykifi4i/kJmfjkinmfme1trf5iZ74uIHz1+/Q8j4kODf//Bx88BndraToWObAH2ovf+OuLhZxjv\nnLuFnwv27mgS2lr7tdbah1prH46Iz0XEf2yt/f2I+FpEvPr4Za9GxFcfP/5aRHwuM9+VmR+OiI9E\nxLfmeekAAAD05tJzQp+ig1+PiNcz8/MR8WZEfDYiorX2Rma+Hg876f4kIr7QxA3Qla2cAXoNo+xA\nb/bwmDWewbK1GTqwR7lE55WZ2+8xoUOHjl/ZqlP935Z/dqBve+6/9navgp611iYbpzM8AQAAKHPp\ndFxgg4ajynsaTbZJEbA1e+jDxz/j+AgxYP0koQAAAJRRhMKOjQ8AN4L8IgkpsEb6phcN713j+xqw\nTopQAAAAylgTCjtke/sXWRsK9E4/Pn2My95/L7BWklAAAADKSEJhJ5yrdh2j6cBamK1x2qFEdPh3\nwPJyiQ4tM/WiUETxeblj/aLfHVBNn3Q990BYVmttssGZjgsAAEAZ03Fho2w+BMDejY9vGf7pvgjL\nkYQCAABQxppQ2BgJ6P1YhwUsSR80D/dJqGNNKAAAAItThMKGjLeiN7p7G78/YI30TbcZ3x8dfQP1\nbEwEG+ActHo2tgDmoiiqMTxTVJ8OtSShAAAAlJGEQsdsrjC/4Ug5wJL08fPITMe3QDFJKAAAAGUk\nodAhCeh6tNb8/oG7MetiGeNZLxJRmJckFAAAgDK5xIhbZhrmgytIQJfn8HhgDvqWdXG/hftorU02\nHkkoAAAAZawJhU44C3Qd7JYLVNHXL2dqjaj3A+5HEQorp/jsh40sgEsZ0Fq3YTGqj4f7MR0XAACA\nMpJQWCmbIqybabnAXPT365OZjm+BO5KEAgAAUEYSCisjAe2fUXLgFLMo+jO1WdHw88D5JKEAAACU\nkYTCSkhA+zRcJwRwK33/+klE4XaSUAAAAMooQmEFxmeBGk3dBgkpMEXfsA3j+7X3Fc6XSzSYzNRK\nIV4uPunbsf7U+wvoI7bNPR1e1lqbbAySUAAAAMrYmAgWYBOibRpvVgFwiv5/O4Yb1dmsCI6ThAIA\nAFBGEgqFJKD7ZVQc9svsiP1wfAucRxIKAABAGUkoFLFr3n5YGwqc4j6wbVOJqPccfkYSCgAAQBlJ\nKEAh64NgP8yGAJiWS3SQmalXZjdMw923U32sawK2R7tnzLMAe9Vam7zgTccFAACgjCIUZtJaM/IJ\nALzwDDB+PoA9UoQCAABQxppQuLNxm5KAEnF8jZhrBLZDW+cUzwnsiTWhAAAALM4RLXAnRjY5Znxw\nObAf7gcMje8Hju5ijxShcAc2IOIWHkCgfwaYuNRUMeo+wF6YjgsAAEAZSShAIdNyYT+kWgDTJKEA\nAACUkYTCDawF5Z6sDYX+mNXArYYzZNwH2AtJKAAAAGVyiRG8zDRsSNccx8K9ONge+qTtMgfPF2xN\na23yIpaEAgAAUMaaULiAEUruLTOtKYMNcV/gFlNnhw4/D1uhCIUzKD5ZgocPWB+DRlRQjLJ1puMC\nAABQRhIKsLDxiDfQHwkVwPkkoQAAAJRRhMIJw3QqM412U05CCuugLVJt/NzhGmQrFKEAAACUsSYU\nDhgnoDA3a0OhP+4PVBjeH+yUyxZIQgEAACgjCQXogJFvWI7ZCQD3lUt0rJmpN2e1xm3CQz9LOdY/\nuy5hftoga+Q5hZ601iYvUNNxAQAAKGM6Ljwyssja2KgI1se9gaWN7w2Wa9AjSSgAAABlJKHsngSU\nHhn5hvmYfUAPphJR9wR6IQkFAACgjCIUYOWMbMM6aIsA92E6Lrs2nHLl4YIemZYL92MaLj0aTst1\nT6AXklAAAADKSELZJQkovXFcCyzHfYIeZKZjW+iGJBQAAIAyilCADZCQwm20IYA6ilAAAADK5BIj\nf5lpuJFFjK93ayXo1bG+23UN59OW2CLPO6xFa23y4pOEAgAAUMbuuOyCEUG2xm65MB/3CHo3vkfY\nLZe1UYQCbIyHDTjNAA7AckzHBQAAoIwklM0bjnZLhtga03Lhftwj2Jqpabmuc9ZAEgoAAEAZSSjA\nRlkbCi8zawBgeZJQAAAAyuQSI4KZaRiS2VkLyt4c68+1AdBGwLMR1VprkxeaJBQAAIAyilA2p7Vm\npA8AYGT4TDR+XoJKNiYC2Ihjx7XYpIg9Mw0XYF0koQAAAJSRhLIZ45Fuo9vs1bFEFHjgHsFeje8R\nZsqwBEkoAAAAZSShADtixJs9MRsAYJ0koQAAAJSRhNI9a0FhmrWh8DL3CHhgbShLkoQCAABQRhIK\nsENGvNky6T/AuuUSHXVmujtwF8Pr18M0HHaqr9d+2ALXOVzPMxVzaK1NXkym4wIAAFDGdFyAHbBJ\nEXsm1QFYF0koAAAAZSShdMm6Bbiv1pq2RNek/HCb4YwZm9cxN0koAAAAZSShdGU80m2EDi5jbSh7\n4h4Bl8vMt+8RElHmoggFICI8bNAnAyoA/TEdFwAAgDKSULpgGi7cl2m5bJl7BNxmfI8wU4Z7k4QC\nAABQRhIKwAuMeNMDKT5AvyShAAAAlMklRhIz0/AlZ7EWFGocuxdod6yJaxXqeR7jWq21yYtFEgoA\nAEAZRSgARrXpnmsYoB82JmK1hlM/PFzAcmxUxBrYiAiWM3Vki3sCt5CEAgAAUEYSCkBEvDzSDT2Q\nxgD0RxIKAABAGUUoAGeRkLIU1x7AtihCAQAAKGNNKKtjV1xYlrWh9MD9AeoN7w92TucWklAAAADK\nSEIBOJuRbypJ4wG2KZfo4DPTXYWXjK9FD7mwDsfuE9opc3DNwfp5buMcrbXJC8N0XAAAAMqYjgvA\nUXvfqGhtP/de04a9/twAWyQJBQAAoIwklMVZUwD96nWjorWlm5e45LV7X4C5jGfJ9Ho/YBmSUAAA\nAMpIQgE4S2Z2lVTd67WubVT/kp9rK7vM9vRaAThNEQrAzVprixYKlxacPRc1l7z2Y7+Xqb/r6T0E\noF+m4wIAAFAmlxh5zEzDnUTEiyPfPScTsDdLTvM8976lTzluLb/HrUwZhr3zTMeU1trkxSAJBQAA\noIw1oQCs0paPIlmDqd/Z1O/80Pvgdw7AtSShAAAAlLEmlEVYNwDbcO/1fOfck/QZteZ6T6wFhW3y\njMeQNaEAAAAszppQAK72NMp97awaadj6HXofhu/d+H289r3zngPsgyIUgFk8FSbDwuJUsaoI6cex\n9/VYUbrEMiAA1sV0XAAAAMrYmIhS95qyBazPNfcTfcA2uRZg3zzv8cTGRAAAACzOmlAAyhgN34fx\n+2wdKABDklAAAADKSEIpYW0A7Jf2jt1xYV/Gx3dN7ZbOvklCAQAAKCMJBeBqx1ItI95MGSckT4b/\n37UDsG2KUAAuNlV8Khy4xLHNi0zdA9g203EBAAAok0tsEJCZdiXYCRsSwbZo01RwncG2aNP71Vqb\nfLMloQAAAJSxJhSAo6z/pNqh4x2GfwdAvyShAAAAlJGEAjBJAsrSpo5zsXMuQP8koQAAAJSxOy6z\nsYYH+mQXQ9bM9Qn98my4P4d2xzUdF4CI8HBPHw5tWuR6BeiH6bgAAACUUYQC8NIUKakSaze+TpdY\nXgTAdRShAAAAlLExEXdn0Tmsn+NX2BrXNPTDs+J+HNqYSBIKAABAGUUowI601iRGbNLUNXzoegdg\nWY5oAdgBx6+wB+Pr2jEuAOskCQUAAKDM0SI0M38uM38vM7+TmW9k5j97/Px7MvMbmfm9zPx6Zr57\n8G++lJnfz8zvZuYn5v4BADjO8SvslWNcANbp5O64mfnnWmt/mpnPIuK/RMQ/johPRcQftdZ+IzN/\nNSJ+vrX2xcz8WET8VkT8YkR8ICJ+JyI+2lr76ei/6S6wYXY8g3XRJtk7bQDWRZvcj6t3x22t/enj\nh++KiHdGxB/HQxH62uPnX4uIzzx+/OmI+Epr7cettTcj4gcR8fHrXzY9GW8AoVOBZQ3bpASUPRte\n/zYrguWNZyhok/tzsgjNzHdk5nci4nlE/G5r7fcj4pXW2vPHL3keEa88fvz+iHhr8M/fiodEFAAA\nAE7vjvs4lfavZOZfjIjfzsy/Nfr7dmJ6raENAAAAIuKC3XFba/83Iv5DRPy1iHieme+NiMjM90XE\njx6/7IcR8aHBP/vg4+cAAADg5O64v/C0821m/tmI+DsR8e2I+FpEvPr4Za9GxFcfP/5aRHwuM9+V\nmR+OiI9ExLfmeOEATLM+G6ZZhwawDqem474vIl7LzHfEQ8H65dbaNzPz2xHxemZ+PiLejIjPRkS0\n1t7IzNcj4o2I+ElEfKHp4Tdv/BZ74IVlaItw2nCDouGf2gvU0hb37eQRLbN8U0e0bIoHX1gHbRHO\np73AOmiL23b1ES0ArN94+q2bOBw3bicmbgHUUYQCAABQ5uQRLQCslw2I4DbDdWnWpAHUkIQCAABQ\nRhEK0CHHsMB9Ob4FoI4iFAAAgDLWhHI1W2oDAHAL54XukyIUoCMGf2A+HoYBapiOCwAAQBlJKEAH\nJKBQRyIKMC9JKAAAAGUUoQArNz6KRRoDNcbtzbEtAPehCAUAAKCMIhQAAIAyucTUksw0n6Vz4+mB\nwP1pZ7Au2iTMTzvbltba5JsoCQUAAKCMIhQAAIAyilAAAADKKEIBVqa1Zk0MrND4uBZHtgBcRxEK\nAABAmWdLvwAAHoxTFQkorM9Tu3xqr09/aq8A55OEAgAAUEYRCgAAQJlcYlF9ZlrJ3ymbpcD9mYYL\n/dJ+YR6eObehtTb55klCAQAAKKMIBQAAoIwiFAAAgDKKUIAFjde8WPcCfRm32yX22gDojSIUAACA\nMopQAAAAyihCAQAAKPNs6RcAsEfOP4NteWrHrbW327e2DTBNEgoAAEAZRSgAAABlcomtxDPT/uWd\nGV8nphjBdbQl2D7tHO5DW+pfa23yTZOEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFDm2dIv\nAGAP7PAH+/HUvp/a/dOf2j3AA0koAAAAZRShAAAAlFGEAgAAUEYRCgAAQBkbEwHMyIZEsF82KAKY\nJgkFAACgjCSUo6Q4AAAswWyC7ZKEAgAAUEYRCgAAQBlFKMBMhtPZM9P0IdipcfsfL3UB2BtFKAAA\nAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUebb0CwDYmvHRLAARP+sPWmtv9xP6CGCP\nJKEAAACUUYQCAABQRhEKAABAGUUoAAAAZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRR\nhAIAAFAmW2v13zSz/ptykfF1kZkLvRLoh3YDnEt/AZfTbvrTWpt8kyShAAAAlFGEAgAAUEYRCgAA\nQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEKAABAGUUoAAAA\nZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACU\nUYQCAABQRhEKAABAGUUoAAAAZRShAAAAlFGEAgAAUEYRCgAAQBlFKAAAAGUUoQAAAJRRhAIAAFBG\nEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEKAABAmWyt1X/TzPpvylXG10dmLvRKoB/DdqPNAGP6\nCLiM59F+tdYm3yxJKAAAAGUUoQAAAJRRhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCAABQRhEK\nAABAGUUoAAAAZRShAAAAlHm29AsA2JrMjIiI1lq01l74HLBfT/1BhD4B2DdJKAAAAGUUoQAAAJRR\nhAIAAFBGEQoAAEAZRSgAAABlFKEAAACUUYQCzGR4BMPwuBZgX8bt3/EswN4pQgEAACijCAUAAKDM\ns6VfAOv2NGXoaRrR05+mEgEAMKfxMhbPn9shCQUAAKCMJBRgRmYTwH5JcQCmSUIBAAAoowgFAACg\njCIUAACAMopQAAAAytiYCKCADYpgP2xIBHCcJBQAAIAyilAAAADKKEIBAAAoowgFAACgjCIUAACA\nMjnewa3km2bWf1PuYni92O0PrqctwTZp23Af2tI2tNYm3zxJKAAAAGUUoQAAAJR5tvQLANijp6lF\nrbW3pxyZbgT9MnUQ4HySUAAAAMooQgEAACijCAUAAKCMIhRgQcO1Y8P1oUAfxu3WelCA0xShAAAA\nlFGEAgAAUEYRCgAAQJlcYv1RZlr01DnrX+D+xv2xtgXrpb3CPDxjbktrbfJNlIQCAABQRhEKAABA\nmWdLvwAAHjxNO3qaivT0p+lIsB6m4QLcThIKAABAGUUowMpk5gvpyhIbyAEvG2+YIgUFuI4iFAAA\ngDKKUAAAAMooQgEAACiTS6w1ykwLnDbADoFQw8HdsCxtEObnuXKbWmuTb6QkFAAAgDKKUAAAAMoo\nQgFWbnxciyNboMa4vZkeCHAfilAAAADKPFv6BQBw2lMC85TKPP0pmYH7s0EKwLwkoQAAAJSRhAJ0\nRCIK85GAAtRQhHI1D8MAANzC4M8+mY4LAABAGUUoQIcy86WjW4DrjY9ikcYAzEcRCgAAQJmTRWhm\nfigzfzczfz8z/2dm/qPHz78nM7+Rmd/LzK9n5rsH/+ZLmfn9zPxuZn5izh8AYM+GiU1rTSIKFxq2\nGwkoQI089cCSme+NiPe21r6TmX8+Iv5bRHwmIv5BRPxRa+03MvNXI+LnW2tfzMyPRcRvRcQvRsQH\nIuJ3IuKjrbWfDv6bnpI2xIJyWN54KiFwHm0HluU5cttaa5Nv6MkktLX2h6217zx+/CcR8QfxUFx+\nKiJee/yy1+KhMI2I+HREfKW19uPW2psR8YOI+PhNrx6Ao8brQyWicNy4nXjwBahz0ZrQzPzLEfFX\nI+L3IuKV1trzx796HhGvPH78/oh4a/DP3oqHohUAAICdO/uc0MepuP8uIn6ltfb/RqPu7cQUW0Py\nG+a8UFgHbRFOM/UP1kFb3LezktDM/DPxUIB+ubX21cdPP39cLxqZ+b6I+NHj538YER8a/PMPPn4O\ngAKOb4FpjmEBWIdzdsfNiPhXEfFGa+1fDP7qaxHx6uPHr0bEVwef/1xmviszPxwRH4mIb93vJQMA\nANCrc3bH/RsR8Z8i4r/Hz6bVfikeCsvXI+IvRcSbEfHZ1tr/efw3vxYR/zAifhIP03d/e/TfNCy/\nQaZVwLrYdAVepE3Aenhu3IdDu+OeLELnoAjdJp0JrIsHbniRNgHr4blxHxShlHGTh3XRJtk7bQDW\nRZvcj6vPCQUAAIB7UYQCbNx4p1y75bIX4+td4gKwDopQAAAAyjxb+gUAML+nBOgpFbIhBFt0KOV3\nfQOsiyQUYEcyc/KB3BRdejd1DR+63gFYliIUAACAMo5oYTY2g4B+aK/0yrUL/dBe98cRLQAAACxO\nEQqAY1zojuNXAPqlCAUAAKCMI1oAiIjDx7hImFgTxwsB9M/GRMzKwwL069CRF1DNtQh98zy4XzYm\nAgAAYHGm4wIwaTw9d/ixUWwqSEABtkkSCgAAQBlJKABHjY9vGf4plWIO1o8BbJskFAAAgDJ2x6WE\nUW3YFmv1uDfXFGyP5z/sjgsAAMDirAkF4GLHds4dfw1MOTYTy7UDsG2KUACuNrVp0fj/KygYOlR8\nuk4A9sN0XAAAAMpIQikxnronIYH9GCZf2vw+LbEJIrAcyzM4RRIKAABAGUkoAGWMju+D5BOAYySh\nAAAAlMklRisz0xDpzlkjBtty7nEbp+45+oM+XfK+OpoFts0zHkOttcmLQBIKAABAGWtCAZjF1Aj4\n+HOHzhY99O9Zj2sTzfFu6QDsjyIUgKvdWkgcK0oP/bcVp7XOeY/v9Z44vgtgH0zHBQAAoIyNiViU\nqXfQp4rNZS65P+k/7mOJ37mNimAbPNMxxcZEAAAALM6aUABW6dBI+lRyNvU5I/HHnZt6+j0CcG+S\nUAAAAMpYE8rixtegUXdYr7Wu37v0XraXfqan38tary3gOM9xHGNNKAAAAIuzJhSAmy098j31/Y8l\naz3vvHuvGUxr+7ky824/GwDrpggF4Cy9FQjHiqxLfpbefu6htRWa12qtbeZnAcB0XAAAAApJQlnc\n0+j2U9rw9KdRb1i/XtvpJa97bUlor7/zc4zvB8B62ZCIW0hCAQAAKCMJBeCovadSRvfXwSwZgO2Q\nhAIAAFBWregsAAAThklEQVQmlxjhzsx9D6tz1PCaNOINyzl2f9A2qeAahHXyrMa5WmuTF4jpuACc\nzcMGlWxUBLBNpuMCAABQRhLK6gxHvm1EAfWkTvTA/QHqmYbLvUhCAQAAKKMIBeAsRr1ZimsPYFsU\noQAAAJSxJhSAiLAWlD5ZGwrQH0koAAAAZSShrFZmvj3CbaQblqPdsQbODIXljNud+wK3koQC4MGe\n7rmGAfqhCAUAAKBMLjFymJmGK7mIaSAwj2P3AO2MNXPtwvw8f3Gr1trkRSMJBQAAoIyNiQB4gZFu\nemCjIoB+SUIBAAAoIwmlC+MRb0e2wG2kR2yZewTcxlpQ5iYJBQAAoIwkFICIMNJNn6wNBeiPIpSu\nTE3L9eAM5/Ogzp6YlguXG94ntB3mYjouAAAAZSShdGmYiBrphttpP/QuMyX9cAMJKJUkoQAAAJSR\nhALsgISIPTNjBmBdJKEAAACUySVGxzPTkDx340BlOO5YP6+9sEWueTif5yjm1FqbvKBMxwXYIQ8Z\nbJmzQwHWzXRcAAAAykhC6d54xNsGFPBACgQvc4+AB6bhsiRJKAAAAGUkoQA7YqSbPbE2FGCdJKEA\nAACUkYSyGdaGwgOpD5zmHsFeWQvKGkhCAQAAKCMJBdiIYwmokW727NjaUIkoQD1JKJuTmS88TJia\nCADw4jPR+HkJKilCAQAAKJNLpESZKZqizHjUD7bIVFw4j7bCHnkWYimttckLThIKAABAGRsTAWyU\n0W542bFNigCoIQkFAACgjCSUzcvMt0e8bcXP1khz4H7cI9ia8T3Ctc1aSEIBAAAoIwkF2Bgj3XCa\ntaEAy1GEsgvjhw1TruidB2eYj3sEvTMNl7UzHRcAAIAyucRoemYawmdRDm2mV8f6bNcyXE/bYgsk\noKxNa23yIpSEAgAAUEYRCrABRrvhNtoQQB1FKAAAAGXsjssuDXfLtQsiPbAbLizHfYIe2O+CnkhC\nAQAAKCMJZdcy09mhdM31CvczPlMaeiABpUeSUICV80AM66AtAtyHIhQAAIAypuOye+PpV6bl0gPX\nJ8zHtFx6ML4+3RfoiSQUAACAMpJQeCQRZW2kMLA+7g0sTQLKFkhCAQAAKJNLjLRnpuF9Vs+W5yzl\nWL/sWoR62iRr4dmE3rTWJi9U03EBOuBhA5ZjoyKA+zIdFwAAgDKSUDhgOPJtIwoqSFmgP+4PVDAN\nl62RhAIAAFBGEQonDEcch6koVDHqDeugLVJt/NzhGmQrFKEAAACUsSYUYGHSdeiftaEA55OEAgAA\nUEYSCmcYnxFnxJsKri9YH2eGUmF8fbkfsDWKULiAYpR78yAL29Jac0/gaopP9sJ0XAAAAMrkEqPw\nmWnon02wbTrXOtb3upagH9oy9+KZgi1qrU1ezJJQAAAAylgTCjcYrhG1PpRbuXagPzYq4lYSUPZI\nEgoAAEAZSShAIWkJ7IcZMgDTJKEAAACUkYTCHWSms0O5mmsF+mdtKJdyJih7pgiFOxk/gChGGfJg\nCvvlfsCQ4hNMxwUAAKBQLjE6n5kiATbPSCcRpxNQ1wVsj3bPFM8F7FFrbfJCl4QCAABQRhEKM8nM\nF0Y5rQkEgH0aPgOMnw9gjxShAAAAlLEmFIqMR0HZvmP9q2sAtk8fgHs/e2dNKAAAAItzTihAISPh\nsB/j86MBeKAIhSKZ+faDiIPLt80DJ3CK+8C2OY4FjjMdFwAAgDKSUCg0npplJHw/vMewX6bl7ocE\nFM4jCQUAAKCMJBQWMJWIGi3tn5QDuJQZMdvhOBY4nyQUAACAMrnEyH1migtgwOhp3xxID1xCn7Et\n7uFwWGttslFIQmEFhjet1pppnRvhYQSYom/YhvH92vsK51OEAgAAUMbGRLASjm/pk9QauCcb1a2f\nY1jgdpJQAAAAykhCYWUkov3zXgGnjPt61k8CCvcjCQUAAKCMJBRWaioRNeq6HtILYC5mwKyPXXDh\nviShAAAAlJGEwsoNE1Gj4+vmfQEuZW3ouklAYR6KUOhEZtqsaAU8KAJV9PXLsQkRzMt0XAAAAMrk\nEqP6mSlKgBsYoa13rK/0+wfuRV+zLPdXuK/W2mQjkoQCAABQxppQ6NDU8S3Dz1PH7xy4p+H6f+pI\nQKGWJBQAAIAyklDo2FQiavT2viQSwFqY9TIPx7BAPUkoAAAAZSShsAHDRNRIeQ2/X2Au41kuzEMC\nCsuRhMKGDG+iw4KU6/j9AWukb7rN+P6oAIV6ilAAAADK5BKjaZlpCA9mZrv56zksHlgTfdJ9uC9C\nvdbaZEOThAIAAFDGxkSwUVPHt4z/DgC2bCpFdg+E5UlCAQAAKGNNKOyE0eDTrLsC1kwfdT73PFiH\nQ2tCTceFnRgf3zL80435ML8bYC2cH3qazYegD6bjAgAAUEYSCjt0aNOivY4YSxWA3unHJaDQE0ko\nAAAAZRShsGOZOblWlAdG0oE10je9aHwEmd8PrJ8iFAAAgDIni9DM/NeZ+Twz/8fgc+/JzG9k5vcy\n8+uZ+e7B330pM7+fmd/NzE/M9cKB+xmOHLfWXvjflu3hZwT2ZQ/92tR9SgIKfTknCf03EfHJ0ee+\nGBHfaK19NCK++fj/IzM/FhG/HBEfe/w3v5mZ0lYAAAAi4owitLX2nyPij0ef/lREvPb48WsR8ZnH\njz8dEV9prf24tfZmRPwgIj5+n5cKzG1qFHmro+qnDn03og6s2al+aqv99tTPpb+G/lybUr7SWnv+\n+PHziHjl8eP3R8Rbg697KyI+cOX3ABbw9GAztWnRFh9qhhSfQG/20G+N7z+H7lNAP26eKtseeoVj\nT6bbfmoFAADgbNcWoc8z870REZn5voj40ePnfxgRHxp83QcfPwd0amvHuOwh0QWI2EZ/5/gV2KZr\ni9CvRcSrjx+/GhFfHXz+c5n5rsz8cER8JCK+ddtLBAAAYCuenfqCzPxKRPzNiPiFzPzfEfFPI+LX\nI+L1zPx8RLwZEZ+NiGitvZGZr0fEGxHxk4j4Qut9CA6IiHjpCJepv+tV768fYNhH98zGQ7APuURn\nlZl995CwY70/IPT++gGO6b2P6/31Ay9qrU024JNJKMBQZr70kDA8LHytek8HAK7Vcx+95tcMXO/m\n3XEBAADgXJJQ4GLjkemnEezxLoZr18NrBLhET2tDTb2F/VKEAjebeuhZy/SvHh7EACqsuV9e+jUB\ntUzHBQAAoIzdcYFZLH2My6m+zag7sHXH+sEl+sCl7wtAvUO740pCAQAAKGNNKDCL8TrRtWxaZOQd\n2Is1bFJk/ScwRRIKAABAGUkoMKslds61Iy7AcUv0wxJQ4IkkFAAAgDKSUKDEcAR8ap3o+Gvm/P4A\ne1KxNvTQf1vfC0xRhALlDj0Q3To9zDRcgMvM1e8qPoFjTMcFAACgjCQUWMx4pHyO41yMxgM8uNe0\nXJsOAbeShAIAAFBGEgqsxrHjXMZfM2QtKMBtjq0NPdbHSkCBa0hCAQAAKJNLJAiZKbYAznJs50Wj\n8wDXOfX8d6yP1b8C52qtTXYYilCgK+f0WR6QAM6jTwXmdKgINR0XAACAMjYmArpyryMGADhOAgrM\nRRIKAABAGUko0KVjxwgMU1Ij+QAvmppJoq8EKklCAQAAKCMJBTZjar3oeMTfaD+wN46zAtZGEgoA\nAEAZSSiwOcOR/XEC8PT/jf4DW3coAdX/AUtThAKbNn7YsnkRsGU2HQJ6YDouAAAAZSShwK7YvAjY\nEpsOAT2ShAIAAFBGEgrs0jmbF42/DmANpJ9A7yShAAAAlJGEArt3aAfd8cfH/g3AvR1LPJ/oi4Ae\nSUIBAAAoIwkFGDm2XnTq85II4F6kn8AeKEIBjjj0sOeIF+AebDIE7JHpuAAAAJSRhAJc4dwjXqa+\nHtgnU20BHkhCAQAAKCMJBbjRsSNeDn1O2gHbJ/kEmCYJBQAAoIwkFODOppIN60Zh26SeAOdThAIU\nMGUXtkXRCXA903EBAAAoIwkFWIApu9APqSfAfUlCAQAAKCMJBViJa9aNHvq3wOXOSTyfaHMA15OE\nAgAAUEYSCrBS56wbPfZ5SQ0cJvUEWI4kFAAAgDKSUICOHEpkLlk/euy/A1tySdoZoV0AVFGEAmzA\nJVN3D/2dB3B6puAE6IfpuAAAAJSRhAJs1LGk59Lpu6f+ezC3S5POCNcswFpJQgEAACgjCQXYoUvX\nkJ7z91InbnFN0jnk+gPohyQUAACAMpJQACLidJJ0a1J6zvdgm25NOSNcOwBbIgkFAACgjCQUgLPc\nmpSe+zXnfj+WdY90c8j7DbAfilAA7uKcIuKSwsWRHHXuXVCOeV8AGDIdFwAAgDKSUADKXJKIXZPO\nzZHorS3Fmzu1PNfafi8A9EMSCgAAQBlJKACrdE3SNkdKuJbkcU5STQAqSUIBAAAoIwkFYDPmSPTW\nloRKLQHonSQUAACAMpJQADhC8ggA9yUJBQAAoIwiFAAAgDKKUAAAAMooQgEAACijCAUAAKCMIhQA\nAIAyilAAAADKKEIBAAAoowgFAACgjCIUAACAMopQAAAAyihCAQAAKKMIBQAAoIwiFAAAgDKKUAAA\nAMooQgEAACijCAUAAKCMIhQAAIAyilAAAADKKEIBAAAoowgFAACgjCIUAACAMopQAAAAyihCAQAA\nKKMIBQAAoIwiFAAAgDKKUAAAAMooQgEAACijCAUAAKCMIhQAAIAyilAAAADKKEIBAAAoowgFAACg\njCIUAACAMopQAAAAyihCAQAAKKMIBQAAoIwiFAAAgDKKUAAAAMooQgEAACijCAUAAKCMIhQAAIAy\nilAAAADKKEIBAAAoowgFAACgjCIUAACAMopQAAAAyihCAQAAKKMIBQAAoIwiFAAAgDKKUAAAAMoo\nQgEAACijCAUAAKCMIhQAAIAyilAAAADKKEIBAAAoowgFAACgjCIUAACAMopQAAAAyihCAQAAKKMI\nBQAAoIwiFAAAgDKKUAAAAMooQgEAACijCAUAAKCMIhQAAIAyilAAAADKKEIBAAAoowgFAACgjCIU\nAACAMopQAAAAyihCAQAAKKMIBQAAoIwiFAAAgDKKUAAAAMooQgEAACijCAUAAKCMIhQAAIAyilAA\nAADKKEIBAAAoowgFAACgjCIUAACAMopQAAAAyihCAQAAKKMIBQAAoIwiFAAAgDKKUAAAAMooQgEA\nACijCAUAAKCMIhQAAIAyilAAAADKKEIBAAAoowgFAACgjCIUAACAMopQAAAAyihCAQAAKKMIBQAA\noIwiFAAAgDKzFKGZ+cnM/G5mfj8zf3WO7wEAAEB/srV23/9g5jsj4n9FxN+OiB9GxH+NiL/XWvuD\nwdfc95sCAACwKq21nPr8HEnoxyPiB621N1trP46IfxsRn57h+wAA/7+9OwixqgzDOP5/HBsKWkQI\nZjqgC4MmiCKIqCgmKCRCW5VB4aJdSdIizBZta1eLaJWFuDCEQIwWJuailSGMVI5SQgMqNG4kighG\nfFvcI94RqyGcc+Z4/7/V971nZu63eLjwzj33PZIk9cxSNKFrgbND+3NNTZIkSZI04paiCfVWW0mS\nJEnSdS1FE3oemBjaTzD4NFSSJEmSNOKWogk9DmxMsj7JOPAicHAJXkeSJEmS1DMrb/QfrKpLSbYD\nh4AxYPfwZNzmZ647JUmSJEmSdHO74Y9okSRJkiTpnyzF7biSJEmSJF1XJ01okk1JTif5OcnOLs4g\nLUaST5PMJflhqHZnksNJfkrydZI7hq7tanJ9Oskz3ZxaWijJRJKjSU4m+THJG03dLKtXktya5FiS\nE0lmkrzX1M2yeifJWJLpJF82e3OskdF6E5pkDPgI2ARMAi8lubftc0iL9BmDrA57GzhcVfcAR5o9\nSSYZDOKabH7n4yTebaDlYB54s6ruAx4BXm/ed82yeqWq/gKmquoB4H5gKsnjmGX10w5ghquPNzTH\nGhldBPhh4ExVzVbVPPA5sKWDc0j/qaq+BS5eU94M7GnWe4Dnm/UWYF9VzVfVLHCGQd6lTlXVr1V1\noln/AZwC1mKW1UNV9WezHGcwAPEiZlk9k2Qd8CzwCXBlYKc51sjoogldC5wd2p9ralJfrK6quWY9\nB6xu1nez8Jm4ZlvLTpL1wIPAMcyyeijJiiQnGGT2aFWdxCyrfz4A3gIuD9XMsUZGF02o43h106jB\neOl/y7R517KR5HbgC2BHVf0+fM0sqy+q6nJzO+464IkkU9dcN8ta1pI8B1yoqmmufgq6gDnWza6L\nJvQ8MDG0n2Dhf3ek5W4uyV0ASdYAF5r6tdle19SkziW5hUEDureqDjRls6zeqqrfgK+AhzDL6pdH\ngc1JfgH2AU8l2Ys51gjpogk9DmxMsj7JOIMvWh/s4BzS/3UQ2NastwEHhupbk4wn2QBsBL7r4HzS\nAkkC7AZmqurDoUtmWb2SZNWViaFJbgOeBqYxy+qRqnqnqiaqagOwFfimql7BHGuErGz7BavqUpLt\nwCEGAwV2V9Wpts8hLUaSfcCTwKokZ4F3gfeB/UleBWaBFwCqaibJfgaT7i4BrzW300hdewx4Gfg+\nyXRT24VZVv+sAfY0k0FXMPhk/0iTa7OsvrqSSd+TNTJihiVJkiRJbfEZQ5IkSZKk1tiESpIkSZJa\nYxMqSZIkSWqNTagkSZIkqTU2oZIkSZKk1tiESpIkSZJaYxMqSZIkSWqNTagkSZIkqTV/AywXmc9Q\ncgTJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x14723a410>"
      ]
     },
     "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": 358,
   "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": 359,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote binned, cropped LDZ array to /astro/opticslab1/SCDA/Apertures/InputMasks/LS/LDZ_quart_ann30D76_X025Pad00_Tol02_N0250.dat\n"
     ]
    }
   ],
   "source": [
    "if aligntol > 0 and LDZ_valid:\n",
    "    if symm is 'half':\n",
    "        LDZ_dat_fname_tail = \"LDZ_half_{0:s}{1:02d}D{2:02d}_{3:s}025Pad{4:02d}_Tol{5:02d}_N{6:04d}.dat\".format(\n",
    "                              shape, int(round(100*inD)), int(round(100*outD)), secobs_key, pad, aligntol, N)\n",
    "    else:\n",
    "        LDZ_dat_fname_tail = \"LDZ_quart_{0:s}{1:02d}D{2:02d}_{3:s}025Pad{4:02d}_Tol{5:02d}_N{6:04d}.dat\".format(\n",
    "                              shape, int(round(100*inD)), int(round(100*outD)), secobs_key, pad, 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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 360,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['LS_quart_ann20D99_clear_N0250.dat',\n",
       " 'LS_quart_ann25D90_clear_N0250.dat',\n",
       " 'LS_quart_ann15D90_clear_N0250.dat',\n",
       " 'LS_quart_ann25D99_clear_N0250.dat',\n",
       " 'LS_quart_ann15D99_clear_N0250.dat',\n",
       " 'LS_half_ann25D81_clear_N0250.dat',\n",
       " 'LDZ_half_ann25D81_clear_Tol02_N0250.dat',\n",
       " 'LDZ_quart_ann25D81_clear_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D81_X025Pad00_N0250.dat',\n",
       " 'LS_quart_ann25D81_X025Pad05_N0250.dat',\n",
       " 'LDZ_quart_ann25D81_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D82_X025Pad00_N0250.dat',\n",
       " 'LS_quart_ann25D82_clear_N0250.dat',\n",
       " 'LS_quart_ann25D98_clear_N0250.dat',\n",
       " 'LS_quart_ann20D90_clear_N0250.dat',\n",
       " 'LS_quart_ann20D81_clear_N0250.dat',\n",
       " 'LS_quart_ann20D82_clear_N0250.dat',\n",
       " 'LS_quart_ann25D81_clear_N0250.dat',\n",
       " 'LS_quart_ann20D98_clear_N0250.dat',\n",
       " 'LS_quart_ann20D76_clear_N0250.dat',\n",
       " 'LS_quart_ann15D98_clear_N0250.dat',\n",
       " 'LS_quart_ann15D82_clear_N0250.dat',\n",
       " 'LS_quart_ann15D81_clear_N0250.dat',\n",
       " 'LS_quart_ann15D76_clear_N0250.dat',\n",
       " 'LS_quart_ann25D76_clear_N0250.dat',\n",
       " 'LDZ_quart_ann25D82_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D76_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann25D76_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann25D90_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann25D90_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D90_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D90_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D90_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann30D90_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D81_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann30D81_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D81_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D81_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D76_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D76_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann20D82_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann20D82_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D82_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann30D82_X025Pad00_Tol02_N0250.dat',\n",
       " 'LS_quart_ann30D76_X025Pad00_N0250.dat',\n",
       " 'LDZ_quart_ann30D76_X025Pad00_Tol02_N0250.dat']"
      ]
     },
     "execution_count": 360,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.listdir(LS_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}