make_luvoir_Lyotstop.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 142,
   "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": 143,
   "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": "code",
   "execution_count": 144,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prim_key_map = {'luvoir15m':'LUVOIR15m'}\n",
    "secobs_key_map = {'Cross':'cross', 'X':'x'}\n",
    "telap_dir = os.path.expanduser(\"~/SCDA/Apertures/GSFC\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set basic parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "overwrite = False\n",
    "N = 256 # pupil array quadrant width after binning\n",
    "iD = 3\n",
    "oD = 95 # _LS_OD_map = {'hex1':76, 'hex2':82, 'hex3':81, 'hex4':82, 'pie08':90,\n",
    "aligntol = 2 # units of thousandths of pupil diameter\n",
    "L = 4096 # starting array width\n",
    "D = 2*N * (L//(2*N)) # starting pupil diameter (forced to be <= L)\n",
    "symm = 'quart' # set to either 'quart' or 'half'\n",
    "#shape = 'ann' \n",
    "shape = \"luvoir15m\" # 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": 146,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Destination /Users/neil/SCDA/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",
    "#LS_dir = os.path.expanduser(\"~/SCDA/InputMasks_wfirst/LS\")\n",
    "LS_dir = os.path.expanduser(\"~/SCDA/InputMasks/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": 147,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAHVCAYAAADPSuPPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+sXOV95/HPt8DSuAkC2xHyD7R2IjcS0NQJVyxurFKS\nG/+gViDSKnKlApFCbAmw03ZXXVCl1fQPpGxVmspmg+L8kI3bhkVqKywrToypm8qSCb2kBDCExdSu\nsONg2d6WRFbZjfPdP+ZxGC537pyZe+Y8P877JV3duWfOzDzPec55Puc595wz5u4CAADx/VLsAgAA\ngC5CGQCARBDKAAAkglAGACARhDIAAIkglAEASAShDABAIghlAAASQSgDAJCIS2MXYJB58+b5lVde\nGbsYQJZOnTpV6/stWrSo1vcD2uLUqVNn3P39g+ZLPpSvvPJKbd68OXYxgOx0Op3a35NtERhNp9P5\nlyrzcfgaKNA4Anmc7wugi1AGACARhDJQmHGPZhktA+NDKAMAkAhCGShIU6NYRsvAeFQOZTO7xMz+\nycz2hr/nm9mTZvZq+H1Vz7wPmNlRM3vFzNb2TL/BzF4Iz20zM6u3OkB7NR2UBDNQv2FGyl+Q9HLP\n3/dLesrdV0h6KvwtM7tW0kZJ10laJ+nLZnZJeM0jkj4vaUX4WTen0gMAUJBK1ymb2VJJvy3pQUl/\nECbfJum3wuNdkv5e0n8L0x9z97ckHTOzo5JuNLPjkq5w96fDez4q6XZJ++qoCNBmw4xaq8xb9f06\nnQ4jZqBGVUfKfy7pDyX9vGfa1e5+8XZBP5Z0dXi8RNLrPfOdCNOWhMfTp7+LmW0ysykzmzp//nzF\nIgIAkLeBoWxmGySddvdn+83j7i7J6yqUu+9w9wl3n5g3b15dbwsUqe5R8rjeE8BgVUbKH5P0qXD4\n+TFJHzezv5D0hpktkqTw+3SY/6Ska3pevzRMOxkeT58OoAGrV68ean7CFmjewFB29wfcfam7L1P3\nBK6/c/fflbRH0l1htrskPREe75G00cwuN7Pl6p7Q9Uw41P2mmd0Uzrq+s+c1AEYwTHBOTk4mUQ4A\n/c3lOuUvSvqkmb0qaTL8LXc/IulxSS9J+rake939QnjNPZK+JumopNfESV7AyJo4xMxhbKBZQ31L\nlLv/vbpnWcvdz0r6RJ/5HlT3TO3p06ckXT9sIQEAaAPu6AVkqMkRLKNloDnJf58yUILcw2rTpk3a\nsWNHpXnrrmvuyw4YBiNloGB1BdrixYtreR8AsyOUgTGLNdLbtGlTre8Xqx6MlNEmhDJQKEa3QH4I\nZWCMShtdllYfIDWEMoChLFu2LMrnEsxoA0IZGJNSR5Wf/exnx/r+QJsRysAYxArkpkaxW7dubeRz\npmO0jNIRykBBmhrFzp8/v5HPAdqGm4cANWvLHbBi3SSk0+lkvdyA2TBSBiJZsGBB7CIkZZhD4oQy\nSkUoAzUaJiy2bNkyvoJkiEPiAKEM1KYth63HiWWItiOUAQBIBKEM1IARXn1YlmgzQhlAclavXh27\nCEAUXBKFVir1blulmJyc1KFDhyrNy/c3oySMlNE6sTrdWHfByhU7TmgjQhloCJf8ABiEUEarMPrK\nC+2FtiGUASSNr4pEmxDKaA1GXXniqyLRJoQyWoFgzBujZbQFoQyMEZ16PRgtoy0IZRSPS6DKwL8f\n0AaEMjAmXAJVDoIZTeGOXiga91EuC3fvQukYKaNYBDIGYR1BaghlAAASQSijSIyAUBXrClJCKANo\nvQULFsQuAiCJUEaBGPlgWFu2bKk8L+sMxolQRmtxHTF6sTOHFHBJFKKL1cFxHTHmos71lpDHRYyU\nERV3aUJKWB8RG6GMaOiIgLexPUAilNFCdH6YDesHYiKUEUWsji/WVwAiL7FOAmSHAIQyWoWvAEQV\nnASIWAhlNI6vUkQOOOkLMRDKaI2SRz/Hjh2LXYRGtKWeaC9CGY1i9DEey5cvj12ERjRdT9ZXNI2b\nh6Ax3DEJOfqN3/gNrVmzprb3q7pudzodtoMWGjhSNrNfNrNnzOwHZnbEzP44TO+Y2Ukzey783Nrz\nmgfM7KiZvWJma3um32BmL4TntpmZjadaAFCPOgMZGKTK4eu3JH3c3X9d0kpJ68zspvDcl9x9Zfj5\nliSZ2bWSNkq6TtI6SV82s0vC/I9I+rykFeFnXX1VQcoYJQNdbAuYzcBQ9q6fhj8vCz8+y0tuk/SY\nu7/l7sckHZV0o5ktknSFuz/t7i7pUUm3z634AJCf1atXxy4CElXpRC8zu8TMnpN0WtKT7v698NQW\nM3vezL5hZleFaUskvd7z8hNh2pLwePr0mT5vk5lNmdnU+fPnh6gOUsTIAHinycnJyvOyTbRLpVB2\n9wvuvlLSUnVHvdereyj6A+oe0j4l6aG6CuXuO9x9wt0n5s2bV9fbInGMHtAm7KxiJkNdEuXu/yrp\noKR17v5GCOufS/qqpBvDbCclXdPzsqVh2snwePp0FGyYzmSY0QMAlKjK2dfvN7Mrw+P3SPqkpB+G\n/xFf9GlJL4bHeyRtNLPLzWy5uid0PePupyS9aWY3hbOu75T0RI11QWIYCQCzYxvBdFWuU14kaVc4\ng/qXJD3u7nvNbLeZrVT3pK/jkjZLkrsfMbPHJb0k6WeS7nX3C+G97pG0U9J7JO0LP0gIGz6Qrrq3\nT7b39AwMZXd/XtJHZph+xyyveVDSgzNMn5J0/ZBlREO4exHQPG4Sgl7cZhOSCEYgpk2bNkX5XLb7\n9BDKiIpOAZAWL14cuwhIBKEMDlsDCWA7hEQoA7Xat68d5y7u3r07dhEa0ZZ6Ih2Ecsuxd16v9evX\nxy5CI+64o+95nkVpup5sjyCUAQAEcyII5RZjrxxID9tHu1W5eQgKxJ2EgHTFukkI10zHx0gZAIBE\nEMotxCgZaBe2+XwQyujrlltuiV0EADUhbPNAKLfMMBvmzTffPL6CAEgWAR4PodwiHMIC2o0+IH2E\nMgAAiSCUW4I95HQdPnw4dhEaQT3TQF+QNq5TThQbQ3usWrUqdhEaQT3TsWzZMh0/frzSvLGumW4r\nRsp4BzYYoHyf/exnYxcBfRDKCYoVjFu3bo3yuQCax21200Qo4xfmz58fuwgA0GqEcmLYewXQFPqb\n9BDKCWFFBdAW9HczI5TBxgG0GNt/WgjlRMTaMJYtWxblcwGkI9ZJnuwQvBuh3HJcGgGAkzzTwc1D\nEsAddgDEFusmIZ1Oh36tByPljCxYsCB2EQAU6MyZM7W/5zCHxAnltxHKkQ2zMm7ZsmV8BQHQWgsX\nLqz9PTkkPhpCOSIOWwMoGX3c8AhlAAASQShHwh4kgDagrxsOoQwAGKvVq1fHLkI2uCSqIu4RCwCj\nmZyc1KFDhyrN2/bvb2akXAF32wKAueGuYdUQygnjblsAShHzEqmcgplQHoDD1gBQD/q1wQhlAEDx\nctkhIJRnwSgZAOpF/zY7QrkPVhwAGI9NmzZF+dwc+nVCOTE5rDSj2r9/f+wiNIJ6lqUt9WzS4sWL\nYxchWYTyDDhsPR5r1qyJXYRGUM+ytKWeTaOfnRmh3FLHjh2LXQQAiCLlYOaOXtO05T6ty5cvj10E\nAC3X9rt3zYSR8ohKaHwgJ/v27YtdhEZQz9GVMKgilHuk2kgApPXr18cuQiOoZ7sNDGUz+2Uze8bM\nfmBmR8zsj8P0+Wb2pJm9Gn5f1fOaB8zsqJm9YmZre6bfYGYvhOe2mZmNp1rDK2EPCwDaLve+vMpI\n+S1JH3f3X5e0UtI6M7tJ0v2SnnL3FZKeCn/LzK6VtFHSdZLWSfqymV0S3usRSZ+XtCL8rKuxLgAA\naMGCBbGLMLKBoexdPw1/XhZ+XNJtknaF6bsk3R4e3ybpMXd/y92PSToq6UYzWyTpCnd/2t1d0qM9\nr4kq9z0rAMDbtmzZUnne1Pr0Sv9TNrNLzOw5SaclPenu35N0tbufCrP8WNLV4fESSa/3vPxEmLYk\nPJ4+fabP22RmU2Y2df78+cqVGTe+qBsA8pBa2FZVKZTd/YK7r5S0VN1R7/XTnnd1R8+1cPcd7j7h\n7hPz5s2r623n7KMf/WjsIgAAapZSgA919rW7/6ukg+r+L/iNcEha4ffpMNtJSdf0vGxpmHYyPJ4+\nPRvbtm2LXQQAQAW5/luyytnX7zezK8Pj90j6pKQfStoj6a4w212SngiP90jaaGaXm9lydU/oeiYc\n6n7TzG4KZ13f2fOaqHJtPADAu23fvr3yvKn16VXu6LVI0q5wBvUvSXrc3fea2WFJj5vZ5yT9i6TP\nSJK7HzGzxyW9JOlnku519wvhve6RtFPSeyTtCz8AANTm7NmzsYswMuv+Ozhdixcv9s2bNzfyWYyY\nASBvqfbjnU7nWXefGDQfd/QCACARhHKPVPewAACDldCHE8rTDNNQU1NT4ysIAKB1COU52Lt3b+wi\nAABUxihZIpRn1HTj7t+/f87vkYO21BNAs4Y5aplyIEuEchLWrFkTuwiNaEs9ATSrpKOWhHIfpRwK\nAYCSldZXE8qzyKEBAQDlIJRrQoADQLNKGyVLhPJAXCIFAGgKoVyjkk42AICUlThKlgjlSkptfADI\nUUmXQE1HKNcstxUAAHIzMTHwex1+Ibc+mVCuILdGBQDkiVAegEPXAJCeUvtmQhkAkKUNGzbELkLt\nCOVZlLonBgAlKPF/y4RyDeayt3bmzJkaS5KuttQTQLNKGzwRyn0M03jD7K1Nt3DhwpFfm5O21BMA\n5oJQnkFpe14AULKS+mxCGQCARBDK05S0xwUAbVFK300oAwCKsGDBgthFmDNCuUcpe1oA0EZbtmyp\nPG+qfTihPIKtW7fGLgIAYAa5D64I5WCYxpk/f/74CgIAaC1CWfnvWQEA3pZzn04oAwCQiNaHco57\nVIcPH45dhEZQz7JQTzQpx75dIpSzbLhVq1bFLkIjqGdZqCealGPfLhHKAAAkg1BWvntUAIB3y7lP\nJ5SDYRrm3Llz4ysIAKC1COURbNu2LXYRAAAzyHmULBHK75B7YwJAm23fvr3yvKn24YQyAKAIZ8+e\njV2EOSOUp2G0DAD5KaXvJpQBAEgEoTyDUva4AKANSuqzCeU+hmm4qampkT/nzJkzI782J9SzLG2p\nJ9C0S2MXoAR79+7VxMTESK9duHBhzaVJE/VsRpOjgNRHHGiHkkbJEiPlWZXW2Cjb7t27G/081nnE\nNsxRylzWV0IZKMRrr70WuwhAo/bu3Ru7CLUjlAdgtIwcxFr3WOcRS6l988BQNrNrzOygmb1kZkfM\n7AthesfMTprZc+Hn1p7XPGBmR83sFTNb2zP9BjN7ITy3zcxsPNUCACA/5u6zz2C2SNIid/++mb1P\n0rOSbpf0GUk/dfc/nTb/tZK+KelGSYslHZD0q+5+wcyekbRV0vckfUvSNnffN9vnL1682Ddv3jxS\n5epU6l4Z8jfb+lbnutjU5wCDTE1NVT50ncq62el0nnX3gWcEDxwpu/spd/9+ePwTSS9LWjLLS26T\n9Ji7v+XuxyQdlXRjCPcr3P1p7+4JPKpuuCePQEaO6l4XWbeRimGudsltvR3qf8pmtkzSR9Qd6UrS\nFjN73sy+YWZXhWlLJL3e87ITYdqS8Hj69Jk+Z5OZTZnZ1Pnz54cpYlQbNmyIXQS0TNMdTr/Py63j\nQ/5KHSxVDmUze6+kv5b0e+7+pqRHJH1A0kpJpyQ9VFeh3H2Hu0+4+8S8efPqetuRDNOYo16rDIzi\n4Ycf7vtcjE4op44PSFWlUDazy9QN5L9097+RJHd/w90vuPvPJX1V3f8hS9JJSdf0vHxpmHYyPJ4+\nvQh0SGhav7tqjXtdZF1HKkocLVc5+9okfV3Sy+7+Zz3TF/XM9mlJL4bHeyRtNLPLzWy5pBWSnnH3\nU5LeNLObwnveKemJmuoxFrk0Iton1XUz1XIBuagyUv6YpDskfXza5U9/Ei5vel7SLZJ+X5Lc/Yik\nxyW9JOnbku519wvhve6R9DV1T/56TdKsZ17HVOIeGMrX1LrIOo9UlNZXVzn7+pC7m7t/2N1Xhp9v\nufsd7v5rYfqnwkj44msedPcPuvuHei95cvcpd78+PHefD7oeqyWOHTsWuwjISOodS+rlQ3lKOsmW\nO3rNoOk9r+XLl8/5PdAOs307UypnYkvSl770peYKgtYr6RIpQnkOSto7Qx5mO+M6hn4d3L/92781\nWxC0XimHsQnlabgECqmabfSZYieTYpmA1A28zWZsTd9ms2pHQoeDQc6cOdPIyDaFdbGJMixdulR3\n33332D8HeUt1xFzbbTbbJIXODeVI7VBz7k6cODF4JiBzhHKQ6t4V8tS2S5PaVl+kK/e+nFAGABRl\n9erVsYswMkJZ+e9ZIS1tHTW2td5Iz+TkZOV5U1ufCOUhbN26NXYRAEnpdSQXNVWu2a7XBqR8B1uX\nxi5ATl544QXdfPPNsYuBhM22cae04Y9TnfXs914PP/xwa5Ynxi+ldYmR8hAOHjwYuwhI2IEDB/o+\nl9JGnxN2cjCqXNcPQln5HuZAWg4dOjTjdNaZuWH5YVg59+mEMlCD1DbstmC5ozSEcpDznhXSxbpS\nD5Yjqsq9LyeUe6TYQEgf601cLP887du3b/BMLUQoj4iOAIOwjtRrtuU520l2SNP69etrf8/cR8kS\nofwuJTQqmtNvHWDdGI9+y7XfSXZoj6mpqcrzprx9EsrAiFLesNuI9mi3vXv3xi5CLQjlGTBaxlyw\nTowXyxfTldRnE8rACFLfsNuKdkHuCOU+StrzQnNYF5rBcsZFpfXVhPIscmhANI/1Im20D3JGKNeE\njgCsA82abXkfPny4uYIgmtJGyRKhPFCJjY7RcQlUWvot9+985zvNFgSNK+USqOkIZaCinDZs0F6l\nK+USqOkI5QoYLWM2tHlcLP/2KXWULBHKlRDKmA1tHhfLv30mJiYqz5vb+kEoAxXltnG3He2FHBHK\nAzBKRhW0fRws9/YqtW8mlIEh5LRxtxnt1A4bNmyIXYTaEcqzKHVPDHOzdu3aGaezDjSL5Y0S/7dM\nKNegxL019Ldq1aq+z+Wy4ZeszjY4c+ZMbe+VspzrWdrgiVDuY5jGG2ZvDWXIYeMuWb/lv3r16lo/\nZ+HChbW+X6raUs8cEMozKG3PC81inRiv2Zbv5ORkcwVBMkrqswllYESpb9xtQ3ugBITyNCXtcSEe\n1o3xYLmin1L6bkIZmIOUN+42oR0g1X9OQQyEco9S9rSQBtaRerE8Mcgw5xSkuj4RyiNItTERx2zr\nA+tKPVjGqCr3wRWhHKTYOMgH608cJRyuBHoRysp/zwppY52ZGy6BwrBy7tMJZaAmqW3cpWN5o0St\nD+Wc96iQD9ad0eS43A4fPhy7CFC+fTuhnGnDIS+sO6PJcbnNdm90NCfXvn1gKJvZNWZ20MxeMrMj\nZvaFMH2+mT1pZq+G31f1vOYBMztqZq+Y2dqe6TeY2QvhuW1mZuOpFtC8lDbsNmB5o0RVRso/k/Rf\n3P1aSTdJutfMrpV0v6Sn3H2FpKfC3wrPbZR0naR1kr5sZpeE93pE0uclrQg/62qsy8hy3aNCHlhn\n5ma25ffwww83VxBkI+c+fWAou/spd/9+ePwTSS9LWiLpNkm7wmy7JN0eHt8m6TF3f8vdj0k6KulG\nM1sk6Qp3f9rdXdKjPa+JbpiG+e53vzu+giBLqW3YbZHzVw4ivhS326H+p2xmyyR9RNL3JF3t7qfC\nUz+WdHV4vETS6z0vOxGmLQmPp0+f6XM2mdmUmU2dP39+mCI24uDBg7GLgIRwY4vxYxmjqtzXh8qh\nbGbvlfTXkn7P3d/sfS6MfL2uQrn7DnefcPeJefPm1fW2A+V8yAPpYR2pF8sTg5TQh1cKZTO7TN1A\n/kt3/5sw+Y1wSFrh9+kw/aSka3pevjRMOxkeT58OZCvVDbttaAeUosrZ1ybp65Jedvc/63lqj6S7\nwuO7JD3RM32jmV1uZsvVPaHrmXCo+00zuym85509r0lGCXtaiI91YzxYruinlL67ykj5Y5LukPRx\nM3su/Nwq6YuSPmlmr0qaDH/L3Y9IelzSS5K+Leled78Q3useSV9T9+Sv1yTtq7MydUm5wZAO1pO0\n0B4oQZWzrw+5u7n7h919Zfj5lrufdfdPuPsKd59093M9r3nQ3T/o7h9y930906fc/frw3H3hf9FZ\noyPAdKwT48UlUpiulFGyxB29+iqpkVNQ2qUrtHlc/ZZ/3evZ7t27a32/VOVcz6mpqcrz5rDdEspo\nxMKFC2MXoTZcnpO2OtvgjjvuqO29UpZzPffu3Ru7CLUilGfBaBnDYB1oFssbJfbRhDIwhFw27Laj\nnZArQnmAEvfEUD/aPg6We3uV2jcTyhXk1KAYH9aDNPVrF9oLOSKUa0ZH0D60OdCs0s647kUoV1Dq\nYRIMZ+3atTNOp83jYvm3z8TEROV5c1s/COUabdiwIXYRMEarVq3q+1xuG34b0CZlK3WwRCgPMExj\nDrP3hjzltHG3Af9PRmkI5VmUuieG8WAdaBbLGyX20YQyMKRcNu62on2QM0K5jxL3wDB+rAvNYDnj\notL6akIZGEEOG3cb0S7tVNJJtoTyDErb88J49PuSDdaJ8WL5YrqSLpEilOeg33WrOTh27FjsImTv\nvvvu6/tc6ht+iVjm7VbKYIpQnmaYxprtutXULV++PHYRipDyxl0iLoFC6QjlHqXsaSENrCP1Ynli\nkBL6cEIZmKNUN+62YPmj14IFC2IXYU4I5aCEPSykh3WlHixHVLVly5bK86a4XhHKQA1S3LjbgOWO\nmaxevTp2EUZGKItRMurBSUjj0W/59bskDZicnKw8b2rbJ6E8hNQaD/lg3RnNbMtttkvSgFwHW5fG\nLkBODhw4MNQeGNqn0+k0NmJOqSPp1US5Uq078pTS+sRIeQiHDh2KXQQgaSl1bmi3XNdFQln5HuZA\nmppaR1JbF9tab6Rn+/btledNbX0ilIExWLp0aewiAK119uzZ2EUYGaEcMFpGne6+++5GPieVdZFR\nMlKRe19OKPdIsYHGYd++fbGL0IjY9WzL+tSUD37wg7GLAIyduXvsMsxq8eLFvnnz5sY+L/e9LJTr\nwIEDfU82jLkuzvbZbCNoUsr9d6fTedbdB37HJCPlaVJuVLTbbJfjxVoXCWSkYmpqqvK8Ka+bhDKQ\nkZQ7k1453+YQedq7d2/sItSCUJ5B06Pl/fv3z/k9gAiH4/o+x0120KRSRskSoZyENWvWxC4CMpJ6\np5J6+VCeUkbJEqHcF/9bRo64NAltU1pfTSjPIocGRDulum6mWi4gF4RyTeiM0LR+J1ONe11kXUcq\nShslS4TyQCU2OsoQ4xIpLoFCKko6uasXoQxkLJXOJpVyoD1KOrmrF6FcAaNl5KjudZF1G6kouU/m\nNpsVlLwCoAyp3mYTGJfc+mVuswm0CF8VCZSBUB4gt70xtFNTXxU5Hes8Yim1byaUgULcd999jX4e\no3PEtmHDhthFqN3AUDazb5jZaTN7sWdax8xOmtlz4efWnuceMLOjZvaKma3tmX6Dmb0QnttmZlZ/\ndepV6p4YyrRw4cJGPy/W6By4aGJi4L9ofyGXPvrSCvPslPSwpEenTf+Su/9p7wQzu1bSRknXSVos\n6YCZ/aq7X5D0iKTPS/qepG9JWicp7rfQ16TEvTXkKZeOB6hLp9OpvN4PM28sA0fK7v4Pks5VfL/b\nJD3m7m+5+zFJRyXdaGaLJF3h7k9793TvRyXdPmqhmzBMww2ztwYAQD9z+Z/yFjN7PhzevipMWyLp\n9Z55ToRpS8Lj6dNnZGabzGzKzKbOnz8/hyKOhsPWAJCPkvrsUUP5EUkfkLRS0ilJD9VWIknuvsPd\nJ9x9Yt68eXW+NQAAyRoplN39DXe/4O4/l/RVSTeGp05KuqZn1qVh2snwePr05JS0xwUAbVFK3z1S\nKIf/EV/0aUkXz8zeI2mjmV1uZsslrZD0jLufkvSmmd0Uzrq+U9ITcyg3AADvcMstt8QuwpxVuSTq\nm5IOS/qQmZ0ws89J+pNwedPzkm6R9PuS5O5HJD0u6SVJ35Z0bzjzWpLukfQ1dU/+ek0Jnnldyp4W\nALTRzTffXHneVPvwKmdf/467L3L3y9x9qbt/3d3vcPdfc/cPu/unwkj44vwPuvsH3f1D7r6vZ/qU\nu18fnrvPU7/p9ixSbUwAaLvcB1fc0StIsXEAAO1CKCv/PSsAwNty7tMJZQBAcRYsWBC7CCNpfSjn\nuEd1+PDh2EVoBPUsS1vqiTRs2bKl8ryp9O0SoZxlKK9atSp2ERpBPcvSlnoiDTn27RKhDABAMghl\n5btHBQB4t5z7dEI5GKZhDhw4ML6CAAAakVogS4TySA4dOhS7CACAGaQYtMMglHvkfMgDANquhD6c\nUAYAIBGE8jQl7GkBQNuU0ncTyjNIucEAAOUilOeIAAeAuEoZJUuEcl9NNfK+fcl9rfRYtKWeAJo1\nNTVVed7UA1kilKNbv3597CI0oi31BNCsvXv3xi5CrQjlWZR0SAQASlNiH00oAwCQCEJ5gBL3xAAg\nd6X2zYRyBTk1KAAgX4RyzQhwABiv0s647kUoV1DqYRIAyNHExETleXPrkwnlGm3YsCF2EQCgFUod\nLBHKAwzTmMPsvQEAMB2hPItS98QAoAQl9tGEMgAAiSCU+yhxDwwASlNaX00oJ+ArX/lK7CI0oi31\nRFn2798fuwgYoKSTbAnlGTS957V58+Y5v0cO2lJPlGXNmjWxi4ABSrpEilCeg7Vr18YuAgBA5RzG\nJpSnGaaxVq1aNb6CAABah1DuUcqeFgC0UQl9OKEMAEAiCOWghD0sAGi73PtyQhkAUJRbbrkldhFG\nRigr/z0rAMDbbr755srzptanE8pDSK3xAAAzy3WwRSgPYfv27bGLAACo4Ny5c5XnJZQzdfbs2dhF\nAABUsG3btthFGAmhrHwPcwAA3i3nPp1QBgAgEYRykPOeFQCgK/e+nFDuMUwDDXMSQWr27dsXuwiN\noJ5l2b17d+wiNIJ6tpu5e+wyzGrx4sXe5Ff+5b6XBQBtlXL/3el0nnX3gd8xeemgGczsG5I2SDrt\n7teHafNEIFAgAAANbElEQVQl/S9JyyQdl/QZd/8/4bkHJH1O0gVJW939O2H6DZJ2SnqPpG9J+oIn\nuEfQ6XQqN1bdjUrIA2iTWH1eyn1tlcPXOyWtmzbtfklPufsKSU+Fv2Vm10raKOm68Jovm9kl4TWP\nSPq8pBXhZ/p7tl7Oh8QBAHM3MJTd/R8kTU+L2yTtCo93Sbq9Z/pj7v6Wux+TdFTSjWa2SNIV7v50\nGB0/2vOa5MTai8r1ujoAGBaj5JmNeqLX1e5+Kjz+saSrw+Mlkl7vme9EmLYkPJ4+HdOkvsIAwFzt\n3LkzdhGSNeezr8PIt9b/DZvZJjObMrOp8+fP1/nWlRGOADAex48fj/K5OfTro4byG+GQtMLv02H6\nSUnX9My3NEw7GR5Pnz4jd9/h7hPuPjFv3rwRizh3HF4BgHrRr85u1FDeI+mu8PguSU/0TN9oZpeb\n2XJ1T+h6JhzqftPMbjIzk3Rnz2sAAIAqhLKZfVPSYUkfMrMTZvY5SV+U9Ekze1XSZPhb7n5E0uOS\nXpL0bUn3uvuF8Fb3SPqauid/vSYpizsesFcHAPWgPx1s4HXK7v47fZ76RJ/5H5T04AzTpyRdP1Tp\nWm6Ya6YBIGU/+tGPYhchC9xmswKCEQDmZseOHVE+N7f+e+BIGV2x7t7FaBlA7lK+/WVqGCkDAJAI\nQjkS9hwBtAF93XAI5YhYAQEAvQjlTBDgAHLDKHl4hHJkrLQASnTgwIHK89K3vY1QBgDU7tChQ7GL\nkCUuiUrAMJc9xbo0C0DZuNtWGhgptxwbBADutpUOQjkRhCOAWLjbVjoI5YRw+AhA09j+00IoAwAa\nx87AzAjlxDBaBtAU+pv0EMr4BU72AIC4COUExdqLjHWyB4DmMUpOE6GMd2CDAcrHdp4ubh6SKG4S\nAiAF9B3NYqTcEtxjG4BEX5A6QhkzOnfuXOwiAEDrEMotMsxe77Zt28ZXEABRMEpOH6HcMsNsaNu3\nbx9fQQA0apijXwRyPIQy+jp79mzsIgCoCUe/8kAotxCHsIB2YZvPB6EMAEAiuE65pTqdTuU9Yq6Z\nBprF3bbai5Fyi7HhA+lh+2g3QhkAwM5AIgjllmO0jFHs27cvdhEa0XQ92R5BKCMa7hqWr/Xr18cu\nQiPaUk+kg1BGtL1krpsE3sYoGRKhjMjoEABp586dsYuARBDKkEQ4AjEdP348yuey3aeHUMYvcPgM\naB7rP3px8xC8AzcKAdLF9lQ+RspIAp0N2oh7UmM6QhljRacDzOy73/1u5XnZNtqDUAaACA4ePBi7\nCEgQoYyxY7QMvBPbBPohlAEASAShjEYwMgC62BYwGy6JQmOG+Q7nhx56aOA8S5YsqfzZGzdurDwv\n0OuLX/yili1bVmnekydPjrcwKB6hjCT95Cc/GTjPD3/4wwZKgrb793//9yjrGqPkduLwNRrFXcPG\n4ytf+UrsIjRi//79jX4e6yuaNqdQNrPjZvaCmT1nZlNh2nwze9LMXg2/r+qZ/wEzO2pmr5jZ2rkW\nHhjGj370o9hFGJvNmzfHLkIj1qxZE7sIwFjVMVK+xd1XuvtE+Pt+SU+5+wpJT4W/ZWbXStoo6TpJ\n6yR92cwuqeHzkZlYo4AdO3ZE+VzkiVEyYhjH4evbJO0Kj3dJur1n+mPu/pa7H5N0VNKNY/h8oC86\nPFRR8lEVpG2uoeySDpjZs2a2KUy72t1Phcc/lnR1eLxE0us9rz0RpqGFCEekLNZRFbYLzDWUV7v7\nSknrJd1rZr/Z+6S7u7rBPRQz22RmU2Y2df78+TkWEani8CBSxPqBmOYUyu5+Mvw+Lelv1T0c/YaZ\nLZKk8Pt0mP2kpGt6Xr40TJvpfXe4+4S7T8ybN28uRQSALLAzAGkOoWxmv2Jm77v4WNIaSS9K2iPp\nrjDbXZKeCI/3SNpoZpeb2XJJKyQ9M+rnowyMlpES1kfENpebh1wt6W/N7OL7/JW7f9vM/lHS42b2\nOUn/IukzkuTuR8zscUkvSfqZpHvd/cKcSo8i1N0hVX2/AwcOaHJystbPRjsQohiXkUPZ3f9Z0q/P\nMP2spE/0ec2Dkh4c9TOBOh06dIhQxi8QtEgBd/RCcbjhP4bFOoNUEMoAACSCUEaRGPmgKtYVpIRQ\nBgAgEYQyisUICIOwjiA1fJ8yitbpdCp3prEuzUJ1XEeM0jFSBsaELzUAMCxCGcXjqyLLwCgZbUAo\nA2NEh14PliPaglBGK9CpYxSsN2gaoYzW4PBnnlh+aBNCGQBmwM4AYiCU0SqMlvNCe6FtCGWgIVwi\nBWAQQhmtwyVS0uHDh2MXYSBGyWgjc/fYZZjV4sWLffPmzbGLAcyK2zXWa/v27Tp79myleVmeyEGn\n03nW3ScGzcdIGUByqgYyUBpCGagBI+X6sCzRZoQyUBMCAsBcEcpABAT4zBglo+0IZaBGhMroDhw4\nUHlelh1KRSgDSMKhQ4diFwGI7tLYBQBK0+l0Ko/k6h7xNTmC5DpioH6MlIGCNBVYO3fubORzgLYh\nlIExKH00d/z48SifW/pyBQhlYExKPbxLMALjQygDyAI7A2gDQhkYo9JGy6XVB0gNoQwUqu6Tsc6d\nO1fr+wF4N0IZGLNYo7y6T8batm1bre9XFaNktAmhDBSsrkDjEiigGdw8BGhAzjcJkYYbdTOyBUbH\nSBnIUJP32OZ+3kBzCGUAABJBKAOZamIEyygZaBahDGRsmCDcvn37UO89zCVQBDJQD0IZaImzZ88O\nNX+sS6CANiOUgcyN4xAzh62BOAhlAAASwXXKQAE6nc5YRsFNvhcARspAMQhIIH+EMoCRsBMA1I9Q\nBgrSVFASyMB4EMoAACSCUAYKM+5RLKNkYHwaD2UzW2dmr5jZUTO7v+nPBwAgVY2GspldIul/Slov\n6VpJv2Nm1zZZBqANxjWaZZQMjFfTI+UbJR1193929/8r6TFJtzVcBqAVcv8OZ6CNzN2b+zCz/yxp\nnbvfHf6+Q9J/cvf7ps23SdKm8Of1kl5srJDNWCjpTOxC1Iw65YE65YE65WGYOv1Hd3//oJmSvKOX\nu++QtEOSzGzK3SciF6lW1CkP1CkP1CkP1Kmapg9fn5R0Tc/fS8M0AABar+lQ/kdJK8xsuZn9B0kb\nJe1puAwAACSp0cPX7v4zM7tP0nckXSLpG+5+ZMDLdoy/ZI2jTnmgTnmgTnmgThU0eqIXAADojzt6\nAQCQCEIZAIBEJBvKOd+O08yOm9kLZvacmU2FafPN7EkzezX8vqpn/gdCPV8xs7XxSv42M/uGmZ02\nsxd7pg1dBzO7ISyLo2a2zcys6br0lGWmOnXM7GRoq+fM7Nae53Ko0zVmdtDMXjKzI2b2hTA927aa\npU7ZtpWZ/bKZPWNmPwh1+uMwPed26lenbNuppzyXmNk/mdne8Hdz7eTuyf2oexLYa5I+IOk/SPqB\npGtjl2uI8h+XtHDatD+RdH94fL+k/xEeXxvqd7mk5aHelyRQh9+U9FFJL86lDpKekXSTJJO0T9L6\nxOrUkfRfZ5g3lzotkvTR8Ph9kv53KHu2bTVLnbJtq/D57w2PL5P0vVCunNupX52ybaeesv6BpL+S\ntDf83Vg7pTpSLvF2nLdJ2hUe75J0e8/0x9z9LXc/JumouvWPyt3/QdK5aZOHqoOZLZJ0hbs/7d21\n9NGe1zSuT536yaVOp9z9++HxTyS9LGmJMm6rWerUTw51cnf/afjzsvDjyrud+tWpn+TrJElmtlTS\nb0v6Ws/kxtop1VBeIun1nr9PaPaNMjUu6YCZPWvdW4ZK0tXufio8/rGkq8PjnOo6bB2WhMfTp6dm\ni5k9Hw5vXzwslV2dzGyZpI+oO2Ipoq2m1UnKuK3CIdHnJJ2W9KS7Z99OfeokZdxOkv5c0h9K+nnP\ntMbaKdVQzt1qd1+p7rdh3Wtmv9n7ZNhzyvpatBLqEDyi7r9JVko6JemhuMUZjZm9V9JfS/o9d3+z\n97lc22qGOmXdVu5+IfQLS9UdTV0/7fns2qlPnbJtJzPbIOm0uz/bb55xt1OqoZz17Tjd/WT4fVrS\n36p7OPqNcEhD4ffpMHtOdR22DifD4+nTk+Hub4SO5eeSvqq3/3WQTZ3M7DJ1w+sv3f1vwuSs22qm\nOpXQVpLk7v8q6aCkdcq8nS7qrVPm7fQxSZ8ys+Pq/tv042b2F2qwnVIN5Wxvx2lmv2Jm77v4WNIa\ndb/lao+ku8Jsd0l6IjzeI2mjmV1uZsslrVD3BIEUDVWHcLjnTTO7KZx5eGfPa5JwcUMLPq23v5Es\nizqFMnxd0svu/mc9T2XbVv3qlHNbmdn7zezK8Pg9kj4p6YfKu51mrFPO7eTuD7j7Undfpm7u/J27\n/66abKcqZ4PF+JF0q7pnXb4m6Y9il2eIcn9A3bPxfiDpyMWyS1og6SlJr0o6IGl+z2v+KNTzFUU+\n67CnTN9U99DT/1P3/yGfG6UOkibU3Shfk/Swwl3kEqrTbkkvSHo+bGCLMqvTanUPpT0v6bnwc2vO\nbTVLnbJtK0kflvRPoewvSvrvYXrO7dSvTtm207T6/ZbePvu6sXbiNpsAACQi1cPXAAC0DqEMAEAi\nCGUAABJBKAMAkAhCGQCARBDKAAAkglAGACAR/x+wGSNT2Ip30gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11ccdfad0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAHVCAYAAADPSuPPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+MXeV95/HPtzYlqAmKabKWY1vBtpwZTLTr1JbjKFEV\nW5sypVFN/olcqQVpCYMEmw7Zxl1mIraJQJBNmthBu0GeNAjTJmEthQqLjRsBO1JUacE7TvhlsI1t\njGBk7I2yEck/3sV894/7uFyG+XF/nHueH+f9kq7mzJl7Z77P+fF8z/e5z7lj7i4AABDf78QOAAAA\ntJCUAQBIBEkZAIBEkJQBAEgESRkAgESQlAEASARJGQCARJCUAQBIBEkZAIBELI0dwGLMjI8cAwDk\n7pfu/sHFnkSlDADA4L3SyZNIygAAJIKkDABAIkjKAAAkgqQMAEAiSMoAACSi46RsZkvM7Bdm9mj4\n/goze8zMXgpfl7U9d9zMTpjZMTO7pm39JjN7LvzsXjOzapsDAEC+uqmUxyS92Pb97ZKecPf1kp4I\n38vMNkjaKelqSSOSvmtmS8Jr7pN0k6T14THSV/QAABSko6RsZqsk/Ymkv2tbvUPSvrC8T9J1besf\ncvfz7v6ypBOStpjZCkmXu/uT7u6SHmx7DYCauPuiDwBxdFop75H015Lealu33N3PhOXXJS0Pyysl\nvdr2vNfCupVhefb6dzGzUTObNrPpDuMDACB7iyZlM/uspHPufni+54TKt7LLa3efdPfN7r65qt8J\nQB1XwVTLQBydfPb1JyX9qZldK+k9ki43s3+QdNbMVrj7mTA0fS48f0bS6rbXrwrrZsLy7PUAajA+\nPt7V891dzMUEatbJ+0tt7zN9WtKjYfmbkm4Py7dL+kZYvlrSM5IulbRG0ilJS8LPDknaKskkHZR0\nbQd/03nw4NH/oxexY+bBo6DHdCd5tp//EvV1SfvN7Ea1Pmj785Lk7kfMbL+kFyS9KelWd78QXnOL\npAckXaZWUj7Yx98H0KFeh6OploF6WervHfGvG4H+9XOek5SBShzuZJ4Un+gFFK7fC+/UL9yBkpCU\nASxqepq7E4E6MHwNFKzK85thbKAvDF8DTVZ1dZv6BTxQAiploFCDOLeploGeUSkDTTWoi+3UL+KB\n3JGUAXRl27ZtsUMAisXwNVCYOs5phrGBrjF8DTRNXVXsyZMna/k7QNNQKQMFqfN8ploGutJRpdzP\nZ18D6FPqF8ULqTp2kjzA8DUQzdDQUOwQksKQOMDwNRBN6udeDFTLKBgTvYBUkZDnxnZB05GUAQBI\nBEkZqBnV4MLYPmgykjKA5IyPj8cOAYiCiV5AjVI/31LCpC8UholeQEq45ac7XMCgiUjKQE3Wrl0b\nOwQAiSMpAzWg6usN2w1NQ1IGkDT+VSSahIlewIClfo7lgElfKAATvQCUgWoZTUGlDAxQ6udXTqiW\nkTkqZSAmboGqFhc4aAKSMjAg3AIFoFtLYwcApIJKLH1V7yOGxJEaKmVAJOSmYr8jNSRlAAASQVJG\n41EtNRv7HykhKQNovKGhodghAJK4TxkNl/rxj/ow6QsDxn3KwEK4jxjtuEBDCkjKaCzuIwaQGpIy\nGomqCHPhuEBsJGUAABJBUkbjUA1hIRwfiImkjEbhXwCiE0wCRCzcEoVGSf14Rzq4RQoV45YooB3V\nD7rBBRxiICmjMUq+BWpqaip2CLVoSjvRXAxfoxFSP86RLoaxURGGr5Evd6/0AfTqy1/+MsciarNo\nUjaz95jZITN7xsyOmNnXwvqvmtmMmT0dHte2vWbczE6Y2TEzu6Zt/SYzey787F7jEhRA4r75zW/G\nDgENsrSD55yXtN3df2tml0j6ZzM7GH62293/tv3JZrZB0k5JV0v6kKTHzewj7n5B0n2SbpL0lKSf\nSBqRdFBAG6oJlMzdGRLHvBatlL3lt+HbS8JjoV5zh6SH3P28u78s6YSkLWa2QtLl7v6kt3rdByVd\n11/4AJCf8fHx2CEgUR29p2xmS8zsaUnnJD3m7k+FH33RzJ41s/vNbFlYt1LSq20vfy2sWxmWZ6+f\n6++Nmtm0mU130RYUgCoZTXD33XfHDgGJ6igpu/sFd98oaZVaVe9H1RqKXitpo6Qzkr5VVVDuPunu\nmzuZqYZyUD2gSbgAxVy6mn3t7r+WNCVpxN3PhmT9lqTvSdoSnjYjaXXby1aFdTNhefZ6QBLVAwB0\nMvv6g2b2/rB8maTPSDoa3iO+6HOSng/LByTtNLNLzWyNpPWSDrn7GUlvmNnWMOv6ekmPVNgWZIyq\nAU3EcY/ZOpl9vULSPjNbolYS3+/uj5rZ35vZRrUmfZ2WdLMkufsRM9sv6QVJb0q6Ncy8lqRbJD0g\n6TK1Zl0z8xoAgIBP9EJ0qR+DwKBxi1Qj8IleAJCD6WluNEELlTKiSv34A+pCtVw8KmWkjYQMvI3z\nARJJGajU2NhY7BBqMTIyEjuEWjSlnUgHw9eIIvXjDoiFYexiMXwNAEBOSMqoHVUyMD/Oj2br5MND\n0HB0EkC9qj7nGBLPB5UyAACJICljQVTJQP44j/NBUsa87rzzztghAKgIiTkP3BKFeaV+bADoDu8t\nR8UtUegdCRkoD+d1+kjKAAAkgqSMd+Fqul579uyJHUItaGcaOL/TxnvKeJfUjwkA/dm+fbumpqZi\nh9E0Hb2nTFLGO6R+PACoBpO+asdEL3Tn5MmTsUMAUBMuwNNEUsa/WLt2bewQAKDRSMqQxFUz0ESc\n9+khKQMAkAiSMrhaBhqM8z8tJOWG27ZtW+wQAETGJM90cEtUw6W+/wHUg1ukBq6jW6KW1hEJqkMS\nBTAIVfctJPneMHydkaGhodghACjQ8ePHK/+dDIn3huHrjKS+rwCgHdXyO/CJXiUhIQPIDf1W90jK\nAAAkgqScAa42AeSK/qs7JGUAwECNj4/HDiEbTPRKXOr7BwA6waQvJnplj0/bAlAKbpHqDJVywlLf\nNwDQjYZXy1TKOSMhAygN/driSMoAACSCpJwgriYBlIr+bWEkZQBAraanp2OHkCyScmJKvorctWtX\n7BBqQTvL0pR21mnTpk2xQ0gWs68Tkvq+AIAqNWw2NrOvMb+pqanYIQAAZqFS7kPq2w4AmiTxyptK\neZBIyEC9xsbGYodQC9rZuxL6ZSrlHqW+3QCgiRKulquplM3sPWZ2yMyeMbMjZva1sP4KM3vMzF4K\nX5e1vWbczE6Y2TEzu6Zt/SYzey787F5LeOsthIQMAGnKvX/uZPj6vKTt7v5vJG2UNGJmWyXdLukJ\nd18v6Ynwvcxsg6Sdkq6WNCLpu2a2JPyu+yTdJGl9eIxU2BYAADQ0NBQ7hJ4tmpS95bfh20vCwyXt\nkLQvrN8n6bqwvEPSQ+5+3t1flnRC0hYzWyHpcnd/0luXMg+2vSYbuV+FAUDpjh49GjuEnnU00cvM\nlpjZ05LOSXrM3Z+StNzdz4SnvC5peVheKenVtpe/FtatDMuz18/190bNbNrMkvrYF/5RNwDkIdcC\nqqOk7O4X3H2jpFVqVb0fnfVzV6t6roS7T7r75k7eFK/TF77whdghAAAK1tUtUe7+a0lTar0XfDYM\nSSt8PReeNiNpddvLVoV1M2F59vpsrFu3LnYIAIAOZDqPuKPZ1x80s/eH5cskfUbSUUkHJN0QnnaD\npEfC8gFJO83sUjNbo9aErkNhqPsNM9saZl1f3/aabOS6owGgKYaHh2OH0LOlHTxnhaR9YQb170ja\n7+6Pmtn/lLTfzG6U9Iqkz0uSux8xs/2SXpD0pqRb3f1C+F23SHpA0mWSDoYHAACVOXbsWOwQesaH\nh/Qo9e0GAE2U8GgmH7MJAEBOSMo9SvhqDAAaqYR+maTch8nJydghAAAKwnvKfUp9+wFAE2RQJfOe\nch2qOBB27dpVQSTpa0o7AdSrpFFLKuUKpL4NAaBkGVTJEpVyfTI5IACgOKX1vyRlAAASQVKuSGlX\nawCQuhL7XZJyhUqabAAAqB8TvSqW+vYEgBJkWCUz0SuGDA8UAMhKyaOSVMoVS317AkAJMiyAqJQB\nAMgJSblCVMkAUI9S+1uSMgAgS3v37o0dQuV4T7kiqW9HAChRRu8t855yXfq5Wjt+/HiFkaSrKe0E\nUK/SCiIq5Qqkvg0BoGSZVMtUynUgIQNAXCX1wyRlAAASQVLuQ0lXZwCQs1L6Y5IyAKAIQ0NDsUPo\nGxO9epT6dgOAJkp40hcTvQbl5MmTsUMAAMwh94KJpNyDtWvXxg4BAFAgknKXcr8KA4DS5dxPk5QB\nAEgESbkLqVx97dmzJ3YItaCdZaGdqFMq/XW3mH3dpdS3FwAgyVnYzL4GACAnJOUuJXj1BQBok3M/\nTVLuwalTp2KHAAAoEEm5B+vWrYsdAgBgDjlXyRJJuWe573gAKM3w8HDsEPpGUgYAFOHYsWOxQ+gb\nSbkPVMsAkIZS+mOSMgAAiSAp96mUqzMAyFVJ/TBJuQKTk5M9v/b48eMVRpIu2lmWprQTqBsfs1mR\n1LcjmqHOioFjHinIqErmYzbrlNGBgUKNjIzU+vc45hFbP6OUqaJSrlDq2xJli5EkOeYRU2YXhlTK\ndcvsAEFBYh17HPOIpdRjb9GkbGarzWzKzF4wsyNmNhbWf9XMZszs6fC4tu0142Z2wsyOmdk1bes3\nmdlz4Wf3WqlbFQCAHizt4DlvSvord/+5mb1P0mEzeyz8bLe7/237k81sg6Sdkq6W9CFJj5vZR9z9\ngqT7JN0k6SlJP5E0IulgNU1Jg5kxpIdaxb625ZhH3Up8L/miRStldz/j7j8Py7+R9KKklQu8ZIek\nh9z9vLu/LOmEpC1mtkLS5e7+pLfO4AclXdd3CxJD5wQAgzU6Oho7hIHp6j1lM7tS0sfUqnQl6Ytm\n9qyZ3W9my8K6lZJebXvZa2HdyrA8e/1cf2fUzKbNbLqb+GLbu3dv7BDQMLGr5ItSiQPNUWoB1HFS\nNrP3SvqxpNvc/Q21hqLXStoo6Yykb1UVlLtPuvvmTmaqpaTkqzek56qrroodwjuQmIH+dZSUzewS\ntRLyD9z9YUly97PufsHd35L0PUlbwtNnJK1ue/mqsG4mLM9eX4RSr9qQrqNHj8YOAYiqxH63k9nX\nJun7kl5092+3rV/R9rTPSXo+LB+QtNPMLjWzNZLWSzrk7mckvWFmW8PvvF7SIxW1A2iUVKvSVOMC\nctHJ7OtPSvoLSc+Z2dNh3YSkPzOzjZJc0mlJN0uSux8xs/2SXlBr5vatYea1JN0i6QFJl6k167qI\nmdclXq0BQA7cvaiLQT7RqwL9bsOpqSlt27atomhQuhw6oNT7FZRlcnJSN998c+wwFtPRJ3qRlPuU\n+vZDWY4fP66hoaHYYSzqwx/+sE6fPh07DDRIBherfMzmoHELFOqWQ0KWpFdeeSV2CGiYUgokKuU+\npL7tUJYrr7wyu2THOYI6JV4tUykPEp0NFnP8+HGZWWWP3BKypErb/4lPfCJ2c5C4EvplKuUepb7d\nEF/iV+1Z4rzDYhI+76iUB4WOAYtJuGPIGtsVi8m9fyYpAwCKMj4+HjuEnjF83aXUtxfio5obPM5D\nLCbB85Dh66qdPHkydggA1JpEBywk1ws3knIXfvjDH8YOAYlL8Oq8SLncrw10i6TchTvuuCN2CEjY\nxMRE7BAahQsgLCTX44P3lHuQ+jZDHLl2AjnjXMRcEj0XeU8ZqEuinUDx2O4oDUm5B3QEAJCm3Ptn\nkjLQp9w7gdyx/fM0NjYWO4QkkZR7REcApINJdvn5zne+U/nvLKFfJin3oYQDAP3hGEjDPffcEzsE\nRDY5ORk7hEow+7pPqW8/DA4JOT2cj82VwfnI7Os6ZHAgAEDRSuqHScpAD0rqBErCfkHuSMoVoCMA\ngDhK639JykCXSusESsP+Qc5IyhWhIwDSsWfPntghoAYl9rsk5QqVeIDgndjHefjSl74UOwQMWCm3\nQM3GLVEVS317onck5PxwPpYrw/ORW6JiyPBAAYCslFolS1TKlUt9e6J3XHDlh/OxXBmej1TKQJXo\n4PPC/kKOSMoVohMAgHqU2t+SlIEulNoRlIb91Ax79+6NHULlSMoVoRNojt27d8cOAYCk0dHR2CFU\njqRcgRKv1jC/2267LXYIWECVF8jHjx+v7HelLOd2llYQMfu6AqlvQwxGhrM/izc+Pq677747dhio\nWSbnYkezr0nKfUp9+2FwMukIGoXzsbkyOB+5JQoYJBJAWtgfKAFJuQ90AgCQhlL6Y5Iy0IdSOoLc\nsR8gteYU5I73lHuU+nZDfTJ4L6t4nI+4KOHzkfeUB4UOAO04HuJi+6Nd7scDSRmoQO4dQa5KGK4E\n2jF83aXUtxfiSXjYrFicj5hPgucjw9dAnUgQ9WJ7o0Qk5S7QCQBYzJ49e2KHAOXbXzN83aXUtxfi\nSnDIrFici1hIgudiNcPXZrbazKbM7AUzO2JmY2H9FWb2mJm9FL4ua3vNuJmdMLNjZnZN2/pNZvZc\n+Nm9luBWA3rF4VwvtjdK1Mnw9ZuS/srdN0jaKulWM9sg6XZJT7j7eklPhO8VfrZT0tWSRiR918yW\nhN91n6SbJK0Pj5EK21ILOgIgHVdddVXsEJCgnPvpRZOyu59x95+H5d9IelHSSkk7JO0LT9sn6bqw\nvEPSQ+5+3t1flnRC0hYzWyHpcnd/0lvjTg+2vSYrd911V+wQkJicO4GcHT16NHYIQKW6muhlZldK\n+pikpyQtd/cz4UevS1oelldKerXtZa+FdSvD8uz1c/2dUTObNrPpbuKryx133BE7BCSEhBwX2x/t\ncj8eOk7KZvZeST+WdJu7v9H+s1D5Vjbrwt0n3X1zJ2+Kx5L7jgeA0pTQL3eUlM3sErUS8g/c/eGw\n+mwYklb4ei6sn5G0uu3lq8K6mbA8ez2QrRI6gRKwH1CKTmZfm6TvS3rR3b/d9qMDkm4IyzdIeqRt\n/U4zu9TM1qg1oetQGOp+w8y2ht95fdtrskRHAABpKKU/XtrBcz4p6S8kPWdmT4d1E5K+Lmm/md0o\n6RVJn5ckdz9iZvslvaDWzO1b3f1CeN0tkh6QdJmkg+EBZKmUTqAUZsa9y8geHx5SgdS3IQaDpJye\n4eFhvfjii7HDQM0yORf57Ou6ZHJARHX8+PHYIVSKfZ6mqm+RGhnJ7qMUepJzOycnJ2OHUCkq5Yqk\nvh1RHRJy+jgfmyOj85FKuU4ZHRgAUIQS+12SMtCFEjuBErGfkCuScoXoCACgHqX2tyRloEOldgKl\nYn8hRyTlitERAMBglTbjuh1JuWLM+izX7t27Y4cAQNLo6GjsEAaGpFyhvXv3xg4BA3TbbbfFDgFd\n4AK5bKXuX+5TrlDq2xLV4C2K9HEuNkNm5yL3KdeJTgAA6lViv0tSBrpUYkdQEvYPckZSrgCdAADE\nUVr/S1IGelBaR1AK9kszlTTJlolefUp9+2Fwrrrqqsr/KxH6w/nYXBlM+mKi16DlfN/q1NRU7BCy\nx//tTQsJudlK2f9Uyn1IfduhHhlcoRePcxFS8ucilfIg0QkAQFpK6JdJykCfSugIcsb2R7uhoaHY\nIfSFpNwDOgEASFPuky9JykAFuFCLg+2OuYyPj8cOoWdM9OpS6tsLcSU+0aQow8PDzIDHvBI8Fzua\n6EVS7kLq2wrxJdgRFIvzEYtJ7Hxk9nXVJiYmYoeAxJEo6sF2RqmolLuU+vZCfIldnReJ8xCLSfA8\npFIehAR3NBJDwhgsti8WMzw8HDuEnpGUgQHYunVr7BCAxjp27FjsEHrG8HWPUt9uiI9Rlepx3mEx\nCZ93DF9jbmNjY7FDqEXsdpJAqnXNNdfEDgEYOCrlPqS+7VCWiYkJ3XPPPbHD6ArnCOqUcJUsUSkP\nXuIHAApz9913xw6hKyRk1GlycjJ2CJUgKQMZySXR5fwxh8jTzTffHDuESpCU+1RFtbxr164KIgHS\nkVtVj7yVUiVLvKdcidS3IcqT8lsnnA+oW8rnQxveU65LJgcEABSntP6XpAxkKNVqNNW4gFyQlCtS\n2tUa0sdkKjRdif0uSblCJR4gSFdqk6moklGnkiZ3tWOiV8VS354oTwoXgxz3qFsKx32XmOgVQ4YH\nCgBkpeR+lqRcMSoG1C32MRf776N5Sj7mSMpAAfhXkUAZeE+5QqlvS5QtxpAexzxiymwYm/eUgSap\n+x+7U50jtr1798YOoXKLVspmdr+kz0o65+4fDeu+KukmSf87PG3C3X8SfjYu6UZJFyT9pbv/NKzf\nJOkBSZdJ+omkMe/gMjuXSpmKASmos3LgmEcKMqqWK6uUH5A0Msf63e6+MTwuJuQNknZKujq85rtm\ntiQ8/z61Evn68Jjrd2apxKs15Mnda3sAKSjtWFw0Kbv7zyT9qsPft0PSQ+5+3t1flnRC0hYzWyHp\ncnd/MlTHD0q6rtegUzM6Oho7BABAAfp5T/mLZvasmd1vZsvCupWSXm17zmth3cqwPHv9nMxs1Mym\nzWy6j/hqUdpVGgDkpqR+uNekfJ+ktZI2Sjoj6VuVRSTJ3SfdfXMn4+8AAJSip6Ts7mfd/YK7vyXp\ne5K2hB/NSFrd9tRVYd1MWJ69PmslXZ0BQM5K6Y97SsrhPeKLPifp+bB8QNJOM7vUzNaoNaHrkLuf\nkfSGmW211lS56yU90kfcAAC8w5133hk7hL51ckvUjyR9WtIHJJ2V9Dfh+42SXNJpSTeHxCsz+4qk\nfyfpTUm3ufvBsH6z3r4l6qCkL+Z8S1QpV2UAUJKEb5Hq6JYoPtGrB6lvMwBoskQTM5/oBQBATkjK\nXaJKBoC05dxPk5QBAMUZGhqKHUJPSMpdSOXqa8+ePbFDqAXtLEtT2ok0HD16NHYIPWGiV5dS314A\ngCQnezHRCwCAnJCUu5Tg1RcAoE3O/TRJuQcTExOxQwAAFIj3lHuU+nYDgCZKuErmPeVBSnjHA0Aj\nldAvk5QBAEgESbkPJVyVAUAJSumPScoAACSCpNynUq7OACBXJfXDJOUK9HNAjI2NVRhJuprSTgD1\nmpycjB1CpbglqiKpb0cAKFFGVTK3RNUpowMDAIpQYr9LUgYAIBEk5QqVeNUGACkqtb8lKQMAkAiS\ncsVKvXoDgFSUNuO6HbOvK5b69gSAEmRYADH7um579+6NHQIANEKpBRCVcoVS35YAUJLMqmUq5TqR\nkAGgXiX2uyRlAAASQVKuQIlXawCQg9L6X5JyAjZvXvRthiI0pZ0oy65du2KHgEWUNMmWiV59Sn37\nAUATZDDpi4leg7Z79+7YIQAAVE6BRKXch9S3HQA0SeLVMpXyIJGQASAtJfTLJGUAABJBUu5BCVdj\nAFCi3PtnkjIAoCh33nln7BB6xkSvLqW+vQAASU76YqJX1UjIAJCHXPtrknIXhoeHY4cAAOjAqVOn\nYofQE4avu5T69gIAMHzdGAnuaABAm5z7aZIyAACJICn3IOerMAAoWe79M0m5R7lOIpCksbGx2CHU\ngnaWZWRkJHYItaCdzcZErz6kvu0AoEkSr5KrmehlZveb2Tkze75t3RVm9piZvRS+Lmv72biZnTCz\nY2Z2Tdv6TWb2XPjZvZb41uuEmVX6AIAmoQ99t06Grx+QNHuc4XZJT7j7eklPhO9lZhsk7ZR0dXjN\nd81sSXjNfZJukrQ+PBi7mCXnIXEAQP8WTcru/jNJv5q1eoekfWF5n6Tr2tY/5O7n3f1lSSckbTGz\nFZIud/cnvTXm+2DbaxCsW7cudggAUItSKtuq9TrRa7m7nwnLr0taHpZXSnq17XmvhXUrw/Ls9ZiF\nAxVA6bZv3x47hGQt7fcXuLtXPRnLzEYljVb5OwEAaZiamoodQrJ6rZTPhiFpha/nwvoZSavbnrcq\nrJsJy7PXz8ndJ919cycz1UpEtQygVPRvC+s1KR+QdENYvkHSI23rd5rZpWa2Rq0JXYfCUPcbZrY1\nzLq+vu01AABAHQxfm9mPJH1a0gfM7DVJfyPp65L2m9mNkl6R9HlJcvcjZrZf0guS3pR0q7tfCL/q\nFrVmcl8m6WB4YB5mxn3QAIpClbw4PjwkcanvHwDoxOHDh7V5cyPfkbyoow8PISknLvX9AwCdoErm\nXzcWgQMZQO7oxzpHUgYAIBEk5QxwlQkgV/Rf3SEpAwCQCJJyJrjaBJAb+q3ukZQzwgEOIBcTExOx\nQ8gSt0RlJvX9BQASRcQcOrolqu9/SIF6VX2gk+QBSCTRVDB83XCciAAOHz4cOwQEDF+DahloOC7O\na8EneqEznJBAc3H+p4WkDABAIkjKkMTVMtBEnPfpISnjXzDZAwDiYqIX3iH14wFANaiSa8dEL3SP\nExUoH+d5ukjKAAAkgqSMd+EqGigX53faSMqY06lTp2KHAACNQ1LGnNatWxc7BAAVo0pOH0kZ8xoe\nHo4dAoCKMPqVB26JwoJSPz4AdIYqOTpuiUL/OJGB/HEe54OkDABAIpbGDgDpq/oqmyFxYGFUts1F\npYza0eEA8+P8aDaSMgAAiSApIwqqgbyNjY3FDqEWdbeT8wLcEoVoTp48qbVr18YOA0gGSblo3BKF\ntPGpYcDbSMiQSMqIjI4IkLZv3x47BCSC4WtEl/oxCAwaF6eNwPA18kCHhCbj+Ec7kjIAAIkgKSMJ\nVAtoIo57zEZSRjLooNAkd911V+wQkCAmeiEpqR+PQFW4CG0cJnohP3RUaAKOc8yHpAwAQCJIykgO\nVQRKxvGNhfD/lJGklStXLvqcj3/84x3/vocffrifcNBgy5Yt07Zt2zp67lNPPTXgaFA6JnqhEVI/\nzpEuKltUhIlewEWld6ybNy96rhdh165dtf690o8bpKevStnMTkv6jaQLkt50981mdoWk/ybpSkmn\nJX3e3f9PeP64pBvD8//S3X/awd+gxEElpqentWnTpthhICMkZVSoo0q5iqS82d1/2bbuG5J+5e5f\nN7PbJS1z9/9oZhsk/UjSFkkfkvS4pI+4+4VF/gZJGZVhGBudIiGjYtGGr3dI2heW90m6rm39Q+5+\n3t1flnRCrQQN1IaOFp04fPhw7BDQUP0mZZf0uJkdNrPRsG65u58Jy69LWh6WV0p6te21r4V1AJCU\nprxHj/QMbVQRAAAHvklEQVT0e0vUp9x9xsz+laTHzOxo+w/d3XsZfg4JfnTRJwI9MDOGsTEvRlMQ\nU1+VsrvPhK/nJP2jWsPRZ81shSSFr+fC02ckrW57+aqwbq7fO+numzsZfwcAoBQ9J2Uz+z0ze9/F\nZUl/JOl5SQck3RCedoOkR8LyAUk7zexSM1sjab2kQ73+faAfVEOYC8cFYuunUl4u6Z/N7Bm1kut/\nd/d/kvR1SZ8xs5ck/dvwvdz9iKT9kl6Q9E+Sbl1s5jUwSBMTE7FDAIB34BO90GipH/+oD1UyBoxP\n9AIWQ0cMieMA6SApAwCQCJIyGo8qqdnY/0gJSRkAgESQlAFRLTUV+x2p6fcTvYBiVN1BM7O7eiRR\nlI5KGRgQ/qkBgG5xnzIwQKmfXzmhSkbmuE8ZiI1EUg22I5qCpAwAQCJIysCAUeX1h+2HJiEpAwCQ\nCJIyUAOqvd6w3dA0JGWgJtwiBWAxJGWgJps3L3o3RG327NkTO4RFUSWjibhPGahZ6udcCoaHh3Xs\n2LHYYQBV4j5lAHkiIaOpSMpAzRiWXRjbB01GUgYAIBEkZSACqsG5sV3QdCRlIBIS0DtNTEzEDgGI\njtnXQESpn3914iIFheto9vXSOiIBMLeqE1GdSZ4kClSP4WugIHUlyu3bt9fyd4CmYfgaKEwd5zRV\nMtA1PjwEaKJBJ0wSMjA4JGUAABJBUgYKNKhqlioZGCySMlCoqidjnTp1qtLfB+DdmOgFFKzK85sq\nGegLE72ApqsqkXILFFAPKmWgcFWc41TJQN+olAH0n1BJyEB9SMoAACSCpAw0QK/VLlUyUC+SMtAQ\nw8PDXT2fW6CA+jHRC2iQbs53qmSgUkz0AvBOnSZaEjIQB0kZAIBELI0dAIB6UQUD6aJSBgAgESRl\nAAASQVIGACARJGUAABJBUgYAIBG1J2UzGzGzY2Z2wsxur/vvAwCQqlqTspktkfRfJf2xpA2S/szM\nNtQZAwAAqaq7Ut4i6YS7n3L3/yvpIUk7ao4BAIAk1f3hISslvdr2/WuSPj77SWY2Kmk0fHte0vOD\nD61WH5D0y9hBVIw25YE25YE25aGbNn24kycl+Yle7j4paVKSzGy6kw/xzgltygNtygNtygNt6kzd\nw9czkla3fb8qrAMAoPHqTsr/S9J6M1tjZr8raaekAzXHAABAkmodvnb3N83s30v6qaQlku539yOL\nvGxy8JHVjjblgTblgTblgTZ1wLr5p+cAAGBw+EQvAAASQVIGACARySblnD+O08xOm9lzZva0mU2H\ndVeY2WNm9lL4uqzt+eOhncfM7Jp4kb/NzO43s3Nm9nzbuq7bYGabwrY4YWb3mpnV3Za2WOZq01fN\nbCbsq6fN7Nq2n+XQptVmNmVmL5jZETMbC+uz3VcLtCnbfWVm7zGzQ2b2TGjT18L6nPfTfG3Kdj+1\nxbPEzH5hZo+G7+vbT+6e3EOtSWAnJa2V9LuSnpG0IXZcXcR/WtIHZq37hqTbw/Ltkv5zWN4Q2nep\npDWh3UsSaMMfSvoDSc/30wZJhyRtlWSSDkr648Ta9FVJX57jubm0aYWkPwjL75N0PMSe7b5aoE3Z\n7qvw998bli+R9FSIK+f9NF+bst1PbbH+B0k/lPRo+L62/ZRqpVzix3HukLQvLO+TdF3b+ofc/by7\nvyzphFrtj8rdfybpV7NWd9UGM1sh6XJ3f9JbR+mDba+p3Txtmk8ubTrj7j8Py7+R9KJan5yX7b5a\noE3zyaFN7u6/Dd9eEh6uvPfTfG2aT/JtkiQzWyXpTyT9Xdvq2vZTqkl5ro/jXOikTI1LetzMDlvr\nI0Mlabm7nwnLr0taHpZzamu3bVgZlmevT80XzezZMLx9cVgquzaZ2ZWSPqZWxVLEvprVJinjfRWG\nRJ+WdE7SY+6e/X6ap01SxvtJ0h5Jfy3prbZ1te2nVJNy7j7l7hvV+m9Yt5rZH7b/MFw5ZX0vWglt\nCO5T622SjZLOSPpW3HB6Y2bvlfRjSbe5+xvtP8t1X83Rpqz3lbtfCP3CKrWqqY/O+nl2+2meNmW7\nn8zss5LOufvh+Z4z6P2UalLO+uM43X0mfD0n6R/VGo4+G4Y0FL6eC0/Pqa3dtmEmLM9enwx3Pxs6\nlrckfU9vv3WQTZvM7BK1ktcP3P3hsDrrfTVXm0rYV5Lk7r+WNCVpRJnvp4va25T5fvqkpD81s9Nq\nvW263cz+QTXup1STcrYfx2lmv2dm77u4LOmP1PovVwck3RCedoOkR8LyAUk7zexSM1sjab1aEwRS\n1FUbwnDPG2a2Ncw8vL7tNUm4eKIFn9Pb/5EsizaFGL4v6UV3/3bbj7LdV/O1Ked9ZWYfNLP3h+XL\nJH1G0lHlvZ/mbFPO+8ndx919lbtfqVbe+R/u/ueqcz91MhssxkPStWrNujwp6Sux4+ki7rVqzcZ7\nRtKRi7FL+n1JT0h6SdLjkq5oe81XQjuPKfKsw7aYfqTW0NP/U+v9kBt7aYOkzWqdlCcl/ReFT5FL\nqE1/L+k5Sc+GE2xFZm36lFpDac9Kejo8rs15Xy3Qpmz3laR/LekXIfbnJf2nsD7n/TRfm7LdT7Pa\n92m9Pfu6tv3Ex2wCAJCIVIevAQBoHJIyAACJICkDAJAIkjIAAIkgKQMAkAiSMgAAiSApAwCQiP8P\nIjC+mrGYbvIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x127a538d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "    if shape == 'luvoir15m':\n",
    "        prim_fname = os.path.join(telap_dir, \"{0:s}_{1:04d}pix.fits\".format(prim_key_map[shape], D))\n",
    "    else:\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",
    "    # Reset inner region to 1 in order to offer a \"blank slate\" for the inner padding array\n",
    "    clear_inside = np.less_equal(Xs**2 + Ys**2, (0.5*D/2)**2)\n",
    "    padded_perim[clear_inside] = 1\n",
    "    # inner aperture padding, computed separately\n",
    "    w = D*inD/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_inner = 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_inner *= perim_shifted\n",
    "    \n",
    "    Lyot_stop = padded_perim.copy()\n",
    "    Lyot_stop *= padded_inner\n",
    "    \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": 148,
   "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": 149,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binning the original LS array 8x\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": 150,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4096, 4096)"
      ]
     },
     "execution_count": 150,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lyot_stop_crop.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(512, 512)"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lyot_stop_binned.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(256, 256)\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": 153,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAJCCAYAAADQsoPKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2QXFWd+P/3ySMUWSWL4Un4LS4JIeEpaABdZIWFlYdd\nebJgAwJZXIECBQywBihAFHkKJIEsgQJZJV/KBbMCSxCQBZVlqVp0EypCeI4QNiAxEUGcFJlkMuf3\nx/SdmYyZOw+3u8/tnverKtV37nTf8+m+pyenP/fT54QYI5IkSdq8YakDkCRJKjMHS5IkSTkcLEmS\nJOVwsCRJkpTDwZIkSVIOB0uSJEk5HCxJkiTlcLAkSZKUw8GSJElSjhGpAwAIITiNuCRJqrffxRjH\n9XUnM0uSJGmoerM/d3KwJEmSlMPBkiRJUg4HS5IkSTkcLEmSJOVwsCRJkpTDwZIkSVKOPgdLIYSd\nQwg/DyG8GEJ4IYRwfmX/lSGEt0MISyv/jur2mEtCCMtDCK+EEA6v5ROQJEmqpf5MStkGXBhjfDaE\n8GfAkhDC45XfzY0x3tj9ziGEycA0YA9gR+CJEMJuMcaN1QxckiSpHvrMLMUY34kxPlvZ/iPwEvDx\nnIccA9wbY2yNMb4BLAf2r0awkiRJ9TagmqUQwi7AvsAvKrvODSE8F0L4XghhbGXfx4GV3R72FvmD\nK0mSpNLq92AphDAGuA/4eozxA+A24C+BKcA7wOyBNBxCODOEsDiEsHggj5MkSaqnfg2WQggj6Rgo\n/SDGeD9AjPG3McaNMcZ24Lt0XWp7G9i528N3quzbRIzxjhjj1Bjj1CJPQJIkqZb68224APwr8FKM\ncU63/Tt0u9txwLLK9iJgWghhdAjhE8AE4JfVC1mSJKl++vNtuAOBU4HnQwhLK/suBU4KIUwBIrAC\nOAsgxvhCCGEh8CId36T7qt+EkyRJjSrEGFPHQAghfRCSJGmoWdKfciBn8JYkScrhYEmSJCmHgyVJ\nkqQcDpYkSZJyOFiSJEnK4WBJkiQpR3/mWaq5vfbai0ceeaRf950xYwYAP/rRj2oZkiRJEmBmSZIk\nKZeDJUmSpBylmMF76tSpcfHixf2678qVKwGYOHEiAB9++GHN4pIkSU3NGbwlSZKKKkWB90DsvPPO\nAMycOROAK6+8MmE0kiSp2ZlZkiRJytFwNUuZlpYWACZPngx01TJJkiT1kzVLkiRJRTVczVJmzJgx\nAFx//fUAnHzyySnDkSRJTcrMkiRJUo6GzSxlpk2bBsCtt94KwNNPP50yHEmS1GTMLEmSJOVo+MxS\nCAGAuXPnAnDAAQcA0N7eniwmSZLUPMwsSZIk5Wj4zFJm6tSOaRKmT58OwPe///2U4UiSpCZhZkmS\nJClH02SWMtdeey0ADzzwAADvv/9+ynAkSVKDM7MkSZKUw8GSJElSjqa7DLfddtsB8O1vfxuA8847\nL2U4kiSpwZlZkiRJyhFijKljYOrUqXHx4sVVPWZbWxsA++67LwDLli2r6vElSVLDWxJjnNrXncws\nSZIk5Wi6mqXMiBEdT23OnDkAfP7zn08ZjiRJalBmliRJknI0bWYp87d/+7cA/P3f/33nvh//+Mep\nwpEkSQ3GzJIkSVKOps8sZebOndu5/cQTTwCwbt26VOFIkqQGYWZJkiQpx5DJLI0fP75zO5vVe9as\nWanCkSRJDcLMkiRJUo6mncE7z/vvvw/ApEmTAFi1alXd2pYkSaXhDN6SJElFDZmape623nprAK66\n6ioAzjjjjJThSJKkEjOzJEmSlMPBkiRJUo6GKPDuHuOGDRsAGDVqVOF2s+Puv//+ANSzyFySJCVn\ngbckSVJRDVHgHULo3K5GRqnncefMmQPA5z73OWDTTJYkSRrazCxJkiTlaIjMUq0ddNBBAFx33XUA\nrFixoi7tZgv5DhtWvTHrQw89BMDq1aurdkxJkoYyM0uSJEk5zCx1841vfCN1CIXdfffdAJx22mmJ\nI5EkqTmYWZIkScrR8Jml1tZWAEaPHp30GGVxyimnAHDrrbcC8Mwzz6QMR5KkhmdmSZIkKUfDZ5aq\nkQ1qhoxSJps7at68eQAccMABgHNHSZI0WGaWJEmScjR8Zkmbt99++wEwffp0AO66666E0UiS1LjM\nLEmSJOUIZahlmTp1aly8eHHqMJrSb37zGwAmTpwIQEtLS8pwJEkqkyUxxql93cnMkiRJUg4HS5Ik\nSTks8G5yO+64IwCXXnrpJreSJKl/zCxJkiTlsMB7iFi3bh0Ae+yxBwCvv/56ynAkSSoDC7wlSZKK\nsmZpiNhiiy0AuOGGGwD44he/mDIcSZIahpklSZKkHKXKLLW2tgJdC9v2/DmVgcRRlph7c/zxxwPw\n+c9/vnPff/7nf6YKR5Kk0jOzJEmSlKNUmaWe2ZiyZGcGEkdZYu7LLbfc0rmdfUNuw4YNqcKRJKm0\nzCxJkiTlKFVmSfUzYcKEzu2vfe1rAMydOzdVOJIklZaZJUmSpBzO4C3ee+89ACZOnAjAmjVrUoYj\nSVK9OIO3JElSUdYsibFjxwJw9dVXA3DmmWemDEeSpFIxsyRJkpTDwZIkSVIOC7zVqa2tDYD99tsP\ngKVLl6YMR5KkWrPAW5IkqaiGL/CuxsK1ZV/8tl5GjOjoDvPmzQPgc5/7HABlyD5KkpSKmSVJkqQc\nDZ9ZqkY2aKhnlHo66KCDADjhhBMAWLhwYcpwJElKysySJElSjobPLKl2rr/+egAWLVoEwLp161KG\nI0lSEmaWJEmScphZUq922WUXAC688EKgazkUSZKGEjNLkiRJOZzBW3368MMPAdhtt90AeOutt1KG\nI0lStTiDtyRJUlENUbPUPfu1YcMGAEaNGlWz9pzRe1NbbrklAD/+8Y8BePXVV+vSbnau29vbq3bM\nH/zgBwD85Cc/qdoxJUnNzcySJElSDmuWNKSsXLkSgEmTJgGwdu3alOFIktKyZkmSJKkoB0uSJEk5\nGqLAW6qWnXfeGYCZM2cCcMUVV6QMR5LUAMwsSZIk5bDAW0NStijwxIkTAfi///u/lOFIktKwwFuS\nJKmoPgdLIYSdQwg/DyG8GEJ4IYRwfmX/n4cQHg8hvFa5HdvtMZeEEJaHEF4JIRxeyycgDcYWW2zB\nFltswQ033MANN9yQOhxJUon1J7PUBlwYY5wMfBr4aghhMnAx8NMY4wTgp5WfqfxuGrAHcARwawhh\neC2ClyRJqrU+vw0XY3wHeKey/ccQwkvAx4FjgIMrd1sAPAnMrOy/N8bYCrwRQlgO7A/8T7WDl4o6\n8cQTAZg/f37nvqeeeipVOJKkEhpQzVIIYRdgX+AXwHaVgRTAKmC7yvbHgZXdHvZWZV/PY50ZQlgc\nQli8Zs2aAYYtSZJUH/2eZymEMAa4D/h6jPGDEELn72KMMYQwoK/VxRjvAO6Ajm/DDeSxUrXdfPPN\nnduf+tSngOou4CtJalz9yiyFEEbSMVD6QYzx/sru34YQdqj8fgdgdWX/28DO3R6+U2WfJElSw+kz\nsxQ6Ukj/CrwUY5zT7VeLgOnAdZXbB7vt/7cQwhxgR2AC8MtqBi1V25QpUzq3v/zlLwNw5513pgpH\nklQi/bkMdyBwKvB8CGFpZd+ldAySFoYQ/gl4EzgRIMb4QghhIfAiHd+k+2qMcWPVI5ckSaoDZ/CW\neli9uuOK8m677QbAH/7wh5ThSJJqxxm8JUmSiur3t+GkoWLbbbcF4NprrwXgnHPOSRmOJCkxM0uS\nJEk5HCxJkiTlsMBb6kVbWxsA++67LwDLli1LGY4kqfos8JYkSSrKAm+pFyNGdLw95s2bB8Chhx4K\nQBmysZKk+jGzJEmSlMPMktSHQw45BIBjjz0WgAceeCBlOJKkOjOzJEmSlMPMktRPs2bNAuDRRx8F\nYN26dSnDkSTViZklSZKkHGaWpH4aP348ADNmzAC6lkORJDU3M0uSJEk5nMFbGqCWlhYAJkyYAMCq\nVatShiNJGjxn8JYkSSrKmiVpgMaMGQN01SydfvrpKcORJNWYmSVJkqQcDZFZ6l5XtWHDBgBGjRpV\ns/ZaW1sBGD16dM3aaEQbN24EYPjw4Q3bXjXP7fTp0wG4/fbbAXjmmWcKH1OSVD5mliRJknI4WJIk\nScrREJfhQgid27W8/Jbx8tvm1evyWy3bq+a5zfrlnDlzADjwwAM7f1eGKTkkSdVhZkmSJClHQ2SW\npDL7zGc+A8App5zSue++++6rebvZlx3a29urdszsWGbGJKmLmSVJkqQcLnciVUn3DM+wYbX/HJK1\nV8227rnnHgBOPvnkqh1TkkrM5U4kSZKKaviapba2NgBGjBj8U6nGMdSYGrn/1CJ7NW3aNABuvfXW\nzn1PP/101duRpEZiZkmSJClHw6dSus/BlPIYQ91glybJHjeYx1aD/WdT2XOZO3du574DDjgAqO63\n7iSpkZhZkiRJytHwmaVqZCNSZDSazWBfw9Svvf1n86ZO7fpyyGmnnQbAXXfdlSgaSUrLzJIkSVIO\n51mSlGvVqlUATJgwAYCWlpaU4UhSNTnPkiRJUlENX7Mkqba23357AC6//HIAZs6cmTIcSao7M0uS\nJEk5HCxJkiTlsMBbUr+sX78egEmTJgHw+uuvpwxHkqrBAm9JkqSiLPCW1C+jRo0CupZCOeaYY1KG\nI0l1Y2ZJkiQpR6kyS62trQCMHj16sz+nMpA4yhJzXwaygO1gF8mtlp7t573GPe+b/ZzV5o0YkabL\n97dfZPfrz31TOfroowE47LDDAHjiiSdShiNJNWdmSZIkKUepMks9P0mX5ZP1QOIoS8x9GUiWKPVC\nsT3bz3uNe943deyZ/vaLRuk/ALfddhsAe+21FwDr1q1LGY4k1YyZJUmSpBylyixJahzjx48H4Lzz\nzgNg1qxZKcORpJoxsyRJkpTDGbwlFfL+++8DXTN7r1q1KmU4kjQQzuAtSZJUlDVLkgrZeuutAbjq\nqqsAOOOMM1KGI0lVZ2ZJkiQphzVLkqoi+1uy//77A+B7WlIDsGZJkiSpKAdLkiRJOSzwllQVIQQA\n5syZA8DnPve5zt+V4XK/JA2WmSVJkqQcDZ9Zam1tBYotQFqNYwx1GzduBAa+cG32uME8thrsP9V3\n0EEHAXDCCSd07lu4cGGqcCSpMDNLkiRJOZw6QFJNvPXWW53bEyZMAGDdunWpwpGkzXHqAEmSpKIa\nvmapmgZbd1Om9urxHOr9OjWKetdflf087LTTTp3bM2fOBOBb3/pWqnAkadDMLEmSJOWwZklSzX34\n4YcATJ48GYAVK1YkjEaSOlmzJEmSVFRD1Cx1z35t2LABgFGjRtWsPefN2bxmqOmq97m1L3XYcsst\nAZg9ezYAd955Z83bXL9+fU2O++abbwKwfPnymhxfUvmYWZIkScphzZIkDcBvfvMbACZOnAhAS0tL\nynAkFWPNkiRJUlENUbOUSnt7O9C1mnp2OxCp58KpRvupn0MzKPIa+vqXy4477gh0zR11+eWXpwxH\nUh2YWZIkScrhYEmSJCmHBd6SNAjZosCTJk0CnGhTalAWeEuSJBVlgbckDcIWW2wBdE20+cUvfjFl\nOJJqyMySJElSDjNLklTA8ccfD8DBBx8MwJNPPpkuGEk1YWZJkiQph5klSaqCm2++GYB9990X6JrU\nVlLjM7MkSZKUw8ySJFXB3nvvDcBZZ50FwG233ZYyHElVZGZJkiQphzN4S1IV/e53vwNgwoQJALz/\n/vspw5GUzxm8JUmSimq6mqUsU5bdDhvW93hw48aNAAwfPrx2gQ2g/d7i2dz+1LEPNI6yxNsfjfKc\nBtJ+6liHgo997GMAXHPNNQCcc845KcORVAVmliRJknJYsyRJNdDW1gbAfvvtB8DSpUtThiNp86xZ\nkiRJKsrBkiRJUo6mK/CWpDIYMaLjz+ucOXMAOPTQQzt/V4byB0n9Z2ZJkiQph5klSaqhQw45BIAT\nTjihc9/ChQtThSNpEMwsSZIk5TCzJEl1cP3113duL1q0CIB169alCkfSAJhZkiRJymFmSZLqYJdd\nduncnjFjBgDXXnttomgkDYSZJUmSpBx9DpZCCN8LIawOISzrtu/KEMLbIYSllX9HdfvdJSGE5SGE\nV0IIh9cq8Hpob2+nvb290DE2btzYuXhpCtVoP/VzaAZFXkNf/+Zz6aWXcumll7LTTjux0047pQ5H\nUh/6k1m6CzhiM/vnxhinVP49AhBCmAxMA/aoPObWEILLm0uSpIbVZ81SjPGpEMIu/TzeMcC9McZW\n4I0QwnJgf+B/Bh1hQsOGFb9KOXx42rFiNdpP/RyaQZHX0Ne/+YwZMwboqlk69dRTU4YjqQ9FRgPn\nhhCeq1ymG1vZ93FgZbf7vFXZ9ydCCGeGEBaHEBavWbOmQBiSJEm1M9jB0m3AXwJTgHeA2QM9QIzx\njhjj1Bjj1HHjxvX7cdWoIypDG42o3q9LLdprhufQyOr5etSqrWrWkH3pS1/iS1/6Ep/+9Kf59Kc/\nXZVjSqq+QQ2WYoy/jTFujDG2A9+l41IbwNvAzt3uulNlnyRJUkMa1DxLIYQdYozvVH48Dsi+KbcI\n+LcQwhxgR2AC8MvCUW7adjUPl6yNRlTv16UW7TXDc2hk9Xw9atVWNWvIshjnzZsHwAEHHABAjLFq\nbUgqrs/BUgjhHuBg4GMhhLeAbwIHhxCmABFYAZwFEGN8IYSwEHgRaAO+GmP0O8+SJKlh9efbcCdt\nZve/5tz/auDqIkFJkiSVRcMtd+JluHSa4RJWMzyHRtYMl+FqYb/99gPgtNNOA2DBggUpw5HUg8ud\nSJIk5Wi4zFItZV8HrtckgLVorx7Pod6vU6Po/nXyerw2nofmc9111wGwxx571KW9rA+1tbVV7Zgr\nV3ZMtXfHHXdU7ZhSamaWJEmScphZ6qben9Br0V49noOZjM1rhv6jtLbffnsA/vmf/zlxJMW9/PLL\nADz11FOJI5GKM7MkSZKUo+EzS9Wo27D2o7jBvob1rvPprX37T1pFXsNsAsfsthoLYA9GkefQjH3o\n5ptvBuBTn/oUgMv+qKGZWZIkScrR8JmlanwSa6ZPc6kM9jVM/drbf8qhyGuYzaeUel6lIs+hGfvQ\nlClTAPjyl78MwJ133pkyHKkQM0uSJEk5QhkWbJw6dWpcvHhx6jAkSVW2evVqACZMmADABx98kDIc\nqaclMcapfd3JzJIkSVKOhq9ZkiSV17bbbgvAN7/5TQAuvPDClOFIg2JmSZIkKYc1S5Kkmlu/fj0A\ne+21FwCvvvpqynCkjDVLkiRJRTlYkiRJymGBtySp5kaNGgXAjTfeCMDRRx+dMhxpQMwsSZIk5TCz\nJGnI6u8Cto204HPZF+X9whe+AMDhhx/eue+xxx5LFY7UL2aWJEmScpQqs9TzE1FZPiE106e6zEA+\nKad+TgPpFz1/l02Nkd0OG5bm80GjZDAGolHeF3mvac+f29vbgT9dnDf1uRhI+6lj7a958+Z1bu+z\nzz4ArFu3LlU4Ui4zS5IkSTlKlVnq61NfKs34qa6RntNA+kXP3/XMEKTS39cw9Ws9EI3ShwbSdqrM\n41C02267dW6fd955AMyaNStVOFIu/zJIkiTlcLkTSVJSLS0tAEyYMAGAVatWpQxHQ4vLnUiSJBVV\nqpolSdLQM2bMGACuvfZaAE4//fSU4Uh/wsySJElSDjNLkqRSmD59OgC33347AM8880zKcKROZpYk\nSZJymFmSJJVCNh/anDlzADjwwAOBrhn4pVTMLEmSJOVwsCRJkpSj4S/DVWOBztQLxTaDwb6GqReO\ntf+UQ5HXsNEWS672Y5vRZz7zGQCmTZsGwD333JMyHMnMkiRJUh6XO5EkldJbb70FwKRJk4CuZVGk\nKnK5E0mSpKIavmapmupdN1CL9lpbWwEYPXp01Y7Z0/r16wEYOXJk577sK79DWb3rrxq1/6Rsr+z8\nG7SpnXbaCYCZM2cCcPnllxc+pjQYZpYkSZJyWLMkSSq1Dz/8EIDJkycDsGLFioTRqMlYsyRJklSU\nNUuSNABZNr4edXr1bKvMttxySwCuvfZaAE466aSU4WgIMrMkSZKUw8zSZjib7uY1wzd1muE5NLJ6\nvh61asvMUjr/8A//AMDBBx8M1L4fbdiwAYD29vaqHfP3v/890LVIsHNHNQYzS5IkSTnMLG2GWYDN\nq/frUov2muE5NLJ6vh61aquea8+lWueurLIM2/bbb584ksHL5o7K5ozK5pBSuflOlCRJyuE8Szmy\n69RFPt2lrlmpxmy6zrJcXCOvSJ/VbWTtD+b9YB8qphp9IHU/KtIHmrH/ZCshZOvevf766ynDGcqc\nZ0mSJKkoB0uSJEk5LPDOUY3iytQFvtVIWzdT6juVIv0gdR/qvmDyYNmHiqlGH0jdj4r0gWbsP6NG\njQJg9uzZABx33HEpw1EfzCxJkiTlMLMkSVIixx57LACHHXZY574nnngiVTjqhZklSZKkHGaWJElK\nbO7cuZ3b++yzD1DdZVZUjJklSZKkHGaWJElKbM899+zcPvvsswGYP39+qnDUg5klSZKkHC53IklS\nibz77rsATJw4cZOfVRMudyJJklSUNUukX2CySPtlWWCyv3GUJd7+6O95KVv/yX4GyDLHI0aM2Ox9\nM6nPy0DaTx1rfw2kX6TuQ9m3rkIIm9z29lp372NZzKnPSzP9Ddpmm20AuOiiiwC45JJLUoYjzCxJ\nkiTlMrNE+jWTirRflk9H/Y2jLPH2R3/PS9n6T148vf0u9XkZSPupY+2vgfSL1H2ot3Uwe3utNxdv\n6vPSjH+DfvrTn6YOQRVmliRJknKYWZIkqUQWLVoEuEZcmZhZkiRJyuFgSZIkKYeX4SRJSmz9+vWd\n2xdeeGHCSLQ5ZpYkSZJymFmSJCmxm266qXN7+fLlCSPR5phZkiRJymFmSZKkRFavXg3AVVddlTgS\n5TGzJEmSlMPMkiRJiVx88cUAtLS0JI5EecwsSZIk5TCzlKO9vR2AEMImtwOxceNGIN1Cma2trUCx\nxSOrcYyhrkg/SN2HNmzYsEn7vS26msc+VEw1+kDqflSkDzRj/1m8eDEACxYsSByJ+sPMkiRJUg4z\nSzkG8wm6p1Sf4jLV+CTWTJ/mUinSD1L3oZEjRxY+hn2omGr0gdT9qEgfaMb+M2PGDKDrCobKzcyS\nJElSDjNLkuqmnnUzqWt0VDv1Ore1aOdHP/oRAE8//XTVjqnaM7MkSZKUw8ySpLqpZ5bHjFLzqte5\nrWY769atA7rmVVJjMbMkSZKUw8GSJElSDi/DSVWSFYNC4xYw13vyv2acbLCIehelN3ofytqqV3tF\n3HLLLQD8+te/ThyJBsPMkiRJUg4zS1KV1LuguBbt1fvTedmzAfVmHypvW4O1atUqAK688sq0gagQ\nM0uSJEk5zCxJSq5I7UyMcZPbaixTpMZT1sWGs6kC1q5dW7Vjqv78qyJJkpTDzJKk5Ip8kg8hbHKr\noalsiw0vWbIEgLvvvrtqx1Q6ZpYkSZJymFmSJKlKstq5r33tawC0t7enDEdVYmZJkiQph5klSZKq\n5Ic//CEAzzzzTOJIVE1mliRJknL0mVkKIXwP+HtgdYxxz8q+Pwd+COwCrABOjDG+V/ndJcA/ARuB\n82KMj9UkckmSSmDdunWd2xdddFHCSFQr/cks3QUc0WPfxcBPY4wTgJ9WfiaEMBmYBuxRecytIYT6\nzt8vSZJURX1mlmKMT4UQdumx+xjg4Mr2AuBJYGZl/70xxlbgjRDCcmB/4H+qE64kSeXyne98p3P7\n7bffThiJamWwNUvbxRjfqWyvArarbH8cWNntfm9V9v2JEMKZIYTFIYTFa9asGWQYkiRJtVW4wDt2\nTCoRB/G4O2KMU2OMU8eNG1c0DEmSpJoY7NQBvw0h7BBjfCeEsAOwurL/bWDnbvfbqbJPkqSmsmLF\nCgBuuummtIGo5gabWVoETK9sTwce7LZ/WghhdAjhE8AE4JfFQpQkSUqnP1MH3ENHMffHQghvAd8E\nrgMWhhD+CXgTOBEgxvhCCGEh8CLQBnw1xrixRrFLSmjjxo63dn8WHx3Ifasta7s/7WdLU7g4b+0N\n5Lyk7D95cdx7770ArF27NllMqo/+fBvupF5+dWgv978auLpIUJIkSWURskX/Upo6dWpcvHhx6jAk\nSeq31as7ynV33333zn3vvfdeqnA0OEtijFP7upPLnUiSJOVwIV1JkgZh2223BeCKK67o3DdjxoxU\n4aiGzCxJkiTlsGZJkqQCun+zb++99wbgxRdfTBWOBsaaJUmSpKKsWZIkqYDu8z/Nnj0bgCOPPDJV\nOKoBM0uSJEk5zCxJklQlRxxxBAB/93d/B8DDDz+cMhxViZklSZKkHGaWJEmqsrlz5wLw2GOPAdDW\n1pYyHBVkZkmSJCmHgyVJkqQcXoaTlFw2qV/3r2D3VzaxbnY7bJifAYeiIn2omsfITJgwAYDzzz8f\n6JpSQI3JvyqSJEk5XO5EUqf29nagftmZeren2qvnOc3aqld7g/HBBx8AMH78eADWrFmTMhz9KZc7\nkSRJKsqaJalBVbO+oucx6/Upvd7tlV0tzmm926vnOe2+gG1Z+9BHPvIRAK655hoAzjjjjJThaJDK\n2bskSZJKwpolSZJqLMuC7b///gA8++yzKcNRF2uWJEmSirJmSVLd1LMmp971P6qfep3baraTHWPW\nrFkAHHbYYYWPqfoxsyRJkpTDzJKkuqlnlseMUvOq17mtRTuHHnooAF/4whcAeOihh6rehqrPzJIk\nSVIOM0s5NmzYAHR9uhjMPB6tra0AjB49unqB1bn91M+hGRSpfUhde9MM74NGV7Z1zwajSB9oxv5z\n4403AvDYY48BsH79+pThqA9mliRJknKYWcoxcuTIwsdI/UmoGu2nfg7NoMin+dS1N83wPmh01fw2\nVipF+kAz9p/ddtsNgPPOOw/oyjSpnMwsSZIk5XCwJEmSlMPlTiRJSqSlpQWAXXfdtXPf6tWrU4Uz\nFLnciSRJUlEWeEuSlMiYMWMAuOaaazr3feUrX0kVjnphZkmSJCmHNUuSJCXW3t7eub3ffvsB8Oyz\nz6YKZyhMg2ZMAAAUD0lEQVSxZkmSJKkoa5YkSUqs+zJCN998MwAHHXRQqnDUg5klSZKkHGaWJEkq\nkc9+9rMAnHjiiQAsXLgwZTjCzJIkSVIuM0tAa2srkG6xxiLtp469bHGksHHjRiDdQqU9289+Bsi+\n7TpixIjN3jeT+vwNpP3UsfbXQPpF6j6UfRMrhLDJbW+vdfc+lsWc+rz0t/3UcQ7ExIkTU4egCjNL\nkiRJOZxnSZKkElm5ciUAkyZNAmDt2rUpw2l2zrMkSZJUlDVLkiSVyDe+8Q3AjFKZmFmSJEnK4WBJ\nkiQph5fhJElK7Omnn+7cvvfeexNGos0xsyRJkpTDzJIkSYlkE4LOmDEjcSTKY2ZJkiQph5klSZIS\nWbBgAQBOzFxuZpYkSZJymFnKsWHDBqBrochhwwY+tky9aGM12k/9HJpBkYVSUy+y2gzvg0ZXjT6Q\nuh81w4Lh1dTS0gLAxRdfnDgS9YeZJUmSpBxmlnKMHDmy8DFSfxKqRvupn0MzKPJpPlUmINMM74NG\nV40+kLofFekDzdh/rrrqKgBWr16dOBL1h5klSZKkHGaWJNVNPetmatVWMzyHRlev16UW7SxfvhyA\nm266qWrHVO2ZWZIkScphZklS3dQzQ1KrtprhOTS6er0utWjnsssuA2D9+vVVP7Zqx8ySJElSDjNL\nkiTV2H/9138B8O///u+JI9FgmFmSJEnK4WBJkiQph5fhJCXTjMtYFFHvqQJq0V49z2nWVr3aG4zs\nNb7gggsAaG9vTxmOBsnMkiRJUg4zS5KSKWs2IJV6TxVQi/bqeU4bof98//vfB+DZZ59NHImKMLMk\nSZKUI8QYU8fAVlttFSdPnjyox2bXf0eNGjXo9k899VQAzjnnnEEfY6gbbO1D9rjBPFbNo0jtTPY3\nLLsdNizNZ8Aiz8FlTYqrxmtYzfPwwQcfALD77rsD8M477xQ+pmpiSYxxal93MrMkSZKUoxSZpRBC\n0iA+8pGPAF0LHI4bNy5lOJKkBnfRRRcBMHv27MSRqA9mliRJkooys9TNV77yFQC++93vJo5EktSI\nXnvtNQCyOty2traU4ahvZpYkSZKKMrPUTfYtmiVLlgAwZcqUlOFIkhrMscceC8CDDz6YOBL1k5kl\nSZKkoswsbcZnP/tZAP77v/87cSSSpLL7yU9+0rl95JFHJoxEg2BmSZIkqSgzSzl++MMfAnDiiScm\njkSSVDbZjN977713574XX3wxVTgaHDNLkiRJRTlYkiRJyuFluBx/8Rd/AcCyZcsAGDNmTMpwkkm9\nyGfq9lVcynOYt1hzz7iyhblDCJvcpl7weSCvX6O8XwbympblOfWM46abbgJgxowZyWJSYV6GkyRJ\nKmpE6gDK7M033wTgO9/5DgDXXXddynCSSf1pLnX7Ki7lOcxru+fvsolpB3KMehhI+6lj7a9GfE49\n43Di4qHDzJIkSVIOa5b6YfTo0QC88MILnft23XXXVOFIkkpk2rRpndvZlDNqGNYsSZIkFWVmaQCO\nO+64zu37778/YSSSpLJYsWJF5/aee+4JwNq1axNFowEysyRJklSUmaVBeuKJJwA49NBDE0ciSSqL\nyy67DICrr746cSTqJzNLkiRJRZlZGqTsuvTSpUuB8swDIklKp6WlBYDdd98dgLfffjtlOOqbmSVJ\nkqSizCwVNH/+fADOOeecxJFIksri3nvvBeCkk05KHIn6YGZJkiSpKDNLBW2zzTYAvPbaawCMHTs2\nZTjJDHZV8NSruasciqwqn/0Ny257W9+t1oo8hyKPVYdqvIbVPA9Zf/yrv/orAJ555pnCx1RNmFmS\nJEkqakSRB4cQVgB/BDYCbTHGqSGEPwd+COwCrABOjDG+VyxMSZKkNApdhqsMlqbGGH/Xbd8s4Pcx\nxutCCBcDY2OMM/s4TsNehsuce+65AMybNy9xJJKksvjFL34BdF2Oa29vTxmO/lSyy3DHAAsq2wuA\nY2vQhiRJUl0UzSy9AfyBjstwt8cY7wghvB9j3Lry+wC8l/2cc5yGzyyNGNFxRfO5554DYNKkSSnD\nkdSA6l3oXYv26vkcGukLIv/4j/8IwIIFC/LvqHrrV2apUM0S8NkY49shhG2Bx0MIL3f/ZYwx9jYQ\nCiGcCZxZsH1JkqSaqtrUASGEK4EW4Azg4BjjOyGEHYAnY4wT+3hsw2eWMocddhgAjz/+eOJIJEll\nsWrVKgDGjx8PwNq1a1OGoy61rVkKIWwVQvizbBv4PLAMWARMr9xtOvDgYNuQJElKbdCZpRDCXwIP\nVH4cAfxbjPHqEMI2wELg/wPepGPqgN/3caymySxlHnywY4x49NFHJ45EKo8U9SzVbiv7m9lRkllb\n9WyrkdSrH9Wineuuuw6ASy65pGrHVCG1rVmKMb4O7LOZ/e8Chw72uJIkSWXicic1kl2XfuGFFwAY\nNWpUynAkVUk2T049llWpZ1uqjz/+8Y8A7LHHHgCsXLkyZThyuRNJkqTizCzV2OWXXw7AiSeemKT9\n7BNpNg/UYIwbNw4YuosEV4OLrKqIsi0SW+/2U8deC/fccw8AJ598cuJIhjwzS5IkSUWZWVKfDjnk\nEAB+9rOfJY5EkppD9n/vX//1XwPw9NNPpwxnKDOzJEmSVJSZJfXbfffdB8Dxxx+fOBJJag6LFy8G\n4IADDujcl30LUnVhZkmSJKkoB0uSJEk5vAynftt1112Brok2R48enTIcSWoap59+euf2XXfdlS6Q\nocfLcJIkSUWZWdKAXXXVVQBcdtlliSORpOawatWqzu0JEyYA0NLSkiqcocTMkiRJUlFmljRgW221\nFQCvvfZa574ddtghVTiS1FRmzZoFwMyZMxNHMiSYWZIkSSrKzJIG7ZRTTuncvvvuuxNGIknNY/36\n9QBMmjQJgNdffz1lOM3OzJIkSVJRZpY0aMOGdY21f/7znwNdi0IOJRs3bgRg+PDhQ7J9bd5Azkvq\nc5gtrxFC2OS2t7iy/d1/l/o59Lf91HEOxLx58wA4//zzE0fS1MwsSZIkFWVmSVWxzz77APDss88C\nm2adJEn919bWBsCUKVOArlUTVBNmliRJkooakToANYdf/epXANx5550AnHnmmSnDkaSGddtttwFm\nlMrEzJIkSVIOa5ZUVePGjQPglVdeAWDs2LEpw5GkhvDuu+92bk+cOPFP9qlmrFmSJEkqysGSJElS\nDgu8VVVr1qwB4KqrrgJgzpw5KcORpIbwrW99q3Pby2/lY2ZJkiQphwXeqonRo0cDXVMKZAWLkqQu\ny5YtA7omoIRNl5NRzVngLUmSVJQ1S6qJ1tZWAL7+9a8D8Oijj6YMR5JKacaMGYDZpLIzsyRJkpTD\nmiXVxcMPPwzAUUcdlTiSNLJPjcOHD6/rY9UcqtEHUvcj3wObWrRoEQDHHHNM4kiGPGuWJEmSijKz\npLrYfffdAXjuuecAGDlyZMpwJCmJ9evXA7DHHnsAsHz58pThyMySJElScX4bTnXx8ssvA/Av//Iv\nAFxwwQUpw9EQUKs6l+ybntlcYrVUz7ZUH7fddhtgRqnRmFmSJEnKYc2S6uqjH/0oAL/+9a8B2Gab\nbVKGI0l1sXr1aqCrfvO9995LGY66WLMkSZJUlDVLqqs//OEPAFx66aUA3H777SnDkUql3vMJ1aK9\nej6H7rNel30Opm9/+9uAGaVGZWZJkiQph4MlSZKkHBZ4K4lhwzrG6YsXL+7ct++++6YKR5Jq4vnn\nnwfgk5/8JABtbW0pw9GfssBbkiSpKAu8lUR7ezsAF154Yee+n/3sZ3WPoxqT/jlxYHFFioKz7Hh2\nm2Ut661IP7APNZ+sP2YT8JpRamxmliRJknJYs6TSePjhhwE46qijEkciScX8x3/8BwDHHXdc4kjU\nB2uWJEmSijKzpNLIlgFYunQpYP2GpMazbt06APbcc0+ga2knlZaZJUmSpKL8NpxK4+WXXwZg/vz5\nQNe3SCSpUdxyyy2AGaVmY2ZJkiQphzVLKp2tt94a6Mo0bbfddinDkaRcq1at6tweP348AGvXrk0V\njgbGmiVJkqSirFlS6bz//vsAXH755QDccccdKcORpFwXX3xx57YZpeZkZkmSJCmHNUsqrWyNryVL\nlgAwZcqUlOFI0iayv037779/575s3Us1DGuWJEmSinKwJEmSlMMCb5VWls5+6qmngMFdhtu4cSMA\n2eXmESPSdPnW1lZg6C7hkp2H4cOHJ2t7c+33PC9ZnwshbHKb3a/7fetpIP1nqPe1enr88ccBL70N\nBWaWJEmScljgrdLKJndbtmwZ4CdlSeXSc9FccJmTBmSBtyRJUlHWLKm0Zs+eDZhRklROW2yxBQA3\n3nhj577jjjsuVTiqITNLkiRJOaxZUukceeSRADzyyCOJI5GkvnX/f/Rv/uZvAHjyyScTRaMBsmZJ\nkiSpKDNLKo1sDqRf/epXAEyePDllOJI0YM8//zwAn/zkJwFoa2tLGY76ZmZJkiSpKL8Np9I4++yz\nATNKkhrXXnvtBcBZZ50FwPz581OGoyoxsyRJkpTDmiUlNXbs2M7tN954A4CPfvSjdWu/GutouRZX\ncUXWjsv+hmW3w4al+QxYpB/Yh5rP6tWrAdh9990BeO+991KGo95ZsyRJklSUNUtK6rLLLuvcrmdG\nKVONT/JmA4obTEYpE0LY5DaVIv3APtR8tt12WwCuuOIKAGbMmJEyHBVkZkmSJCmHgyVJkqQcFngr\nifHjxwPw4osvdu4bOXJkqnCkUihS6N4I7Q1F2Wu89957A5v+zVMpWOAtSZJUlAXeSmLu3LmA2SSp\nu3pneMwo1V72Gs+ePRvoWihcjcXMkiRJUg5rllRXRxxxBACPPvpo4kgkqf6OPfZYAB588MHEkajC\nmiVJkqSizCypLkaM6CiPe+655wCYNGlSynA0BNTqm171XJrEZVCaT/Y3cN999wWgvb09ZTgysyRJ\nklSc34ZTXZx77rmAGSXVT62+6VXPLI8ZpeaTzbd09tlnAzB//vyU4aifzCxJkiTlsGZJNTVu3DgA\nXnnlFQDGjh2bMpxkitTPOMuyqtEH7Efl8u677wIwceLETX5W3VmzJEmSVJQ1S6qpK6+8Ehi6GaVM\nkU/zZgJUjT5gPyqXbbbZBoArrrgCgPPPPz9lOOqDmSVJkqQc1iypJvbZZx8AlixZAvipVpI2Z8OG\nDUDXvEsAL7zwQqpwhiJrliRJkopysCRJkpTDAm/VxE033QR4+U2S8owcORKAOXPmdO47/PDDU4Wj\nXphZkiRJymGBt6rquOOOA+D+++9PHIkkNaajjz4agIceeihxJENC2gLvEMIRIYRXQgjLQwgX16od\nSZKkWqpJZimEMBx4Ffhb4C3gf4GTYowv9nJ/M0sNLlvw86WXXgLgE5/4RMpwJKlhvfrqqwDstdde\nAKxfvz5lOM0uaWZpf2B5jPH1GON64F7gmBq1JUmSVDO1+jbcx4GV3X5+CzigRm2pBC644AJgaGaU\nUi9Qmrp9bd5Azkvqc9je3g5ACGGT295k8YL9rhZ22203AM4991wAZs+enTIckXDqgBDCmcCZqdqX\nJEnqj1rVLH0GuDLGeHjl50sAYozX9nL/NcBa4HdVD0a18DE8V43A89QYPE+NwfPUGAZ6nv4ixjiu\nrzvVarA0go4C70OBt+ko8D45xtjrgjchhMX9KbJSep6rxuB5agyep8bgeWoMtTpPNbkMF2NsCyF8\nDXgMGA58L2+gJEmSVFY1q1mKMT4CPFKr40uSJNVDmZY7uSN1AOo3z1Vj8Dw1Bs9TY/A8NYaanKdS\nLHciSZJUVmXKLEmSJJVOKQZLriNXXiGEFSGE50MIS0MIiyv7/jyE8HgI4bXK7djUcQ41IYTvhRBW\nhxCWddvX63kJIVxSeX+9EkI4PE3UQ08v5+nKEMLblffU0hDCUd1+53lKIISwcwjh5yGEF0MIL4QQ\nzq/s9z1VIjnnqebvqeSX4Qa6jpzqK4SwApgaY/xdt32zgN/HGK+rDG7HxhhnpopxKAoh/DXQAvy/\nGOOelX2bPS8hhMnAPXQsQ7Qj8ASwW4xxYy+HV5X0cp6uBFpijDf2uK/nKZEQwg7ADjHGZ0MIfwYs\nAY4F/hHfU6WRc55OpMbvqTJkllxHrvEcAyyobC+go7OqjmKMTwG/77G7t/NyDHBvjLE1xvgGsJyO\n951qrJfz1BvPUyIxxndijM9Wtv8IvETHsl2+p0ok5zz1pmrnqQyDpc2tI5f35FVfEXgihLCkskQN\nwHYxxncq26uA7dKEph56Oy++x8rn3BDCc5XLdNmlHc9TCYQQdgH2BX6B76nS6nGeoMbvqTIMllRu\nn40xTgGOBL5auazQKXZcx/UrlSXjeSm124C/BKYA7wCukloSIYQxwH3A12OMH3T/ne+p8tjMear5\ne6oMg6W3gZ27/bxTZZ9KIMb4duV2NfAAHSnM31auHWfXkFeni1Dd9HZefI+VSIzxtzHGjTHGduC7\ndF0W8DwlFEIYScd/wD+IMd5f2e17qmQ2d57q8Z4qw2Dpf4EJIYRPhBBGAdOARYljEhBC2KpSREcI\nYSvg88AyOs7P9MrdpgMPpolQPfR2XhYB00IIo0MInwAmAL9MEJ/o/E83cxwd7ynwPCUTQgjAvwIv\nxRjndPuV76kS6e081eM9VbPlTvrLdeRKbTvggY7+yQjg32KMPwkh/C+wMITwT8CbdHwTQXUUQrgH\nOBj4WAjhLeCbwHVs5rzEGF8IISwEXgTagK/6rZ366OU8HRxCmELHJZ0VwFngeUrsQOBU4PkQwtLK\nvkvxPVU2vZ2nk2r9nko+dYAkSVKZleEynCRJUmk5WJIkScrhYEmSJCmHgyVJkqQcDpYkSZJyOFiS\nJEnK4WBJkiQph4MlSZKkHP8/H7QVSF6BrX8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11ab18dd0>"
      ]
     },
     "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": 154,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote binned, cropped LS array to /Users/neil/SCDA/InputMasks/LS/LS_quart_luvoir15mP03D95_clear_N0256.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": 155,
   "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": 156,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lyot plane dark zone width in binned array: +/- 2 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": 157,
   "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": 158,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False, False,  True, False, False],\n",
       "       [False,  True,  True,  True, False],\n",
       "       [ True,  True,  True,  True,  True],\n",
       "       [False,  True,  True,  True, False],\n",
       "       [False, False,  True, False, False]], dtype=bool)"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "allowed_shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])\n",
    "len(XYshifts_allowed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "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": 161,
   "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": 162,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x127e5bf50>"
      ]
     },
     "execution_count": 162,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5EAAAOICAYAAABPC3XsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3U2odVeaF/Dn6bxttWiLVXQbQhLsGmSSEmwhFA3tQAWt\n+IGpURFByaAgkxJKECRxIj3rkTixBkEbA34UAW0qNI0SY4MT7XRiV9OddIcKVheVkKrYiPgxKKli\nOXjPMTfJOes8597zsdbevx+E9959z33vfm/22Xv/9/OstbK1FgAAAFDxY9feAQAAAOYhRAIAAFAm\nRAIAAFAmRAIAAFAmRAIAAFAmRAIAAFAmRAIAAFAmRAIAAFAmRAIAAFB279o7EBGRme3a+wAAALBy\nf9Ba++lDL1KJBAAAICLiO5UXCZEAAACUCZEAAACUCZEAAACUCZEAAACUCZEAAACUCZEAAACUCZEA\nAACUCZEAAACUCZEAAACUCZEAAACUCZEAAACUCZEAAACUCZEAAACU3bv2DgDADFprd/r+zDzRngDA\ndZUqkZn5+5n525n5zcx8fbPtM5n5SmZ+a/Pnp2+8/vnMfCcz387ML5xr5wEAALisY9pZ/3xr7Wdb\na09sPn8uIl5trT0WEa9uPo/MfDwino6Iz0XEkxHxtcx84IT7DAAX01q7cxXylH8PAFzbXcZEPhUR\nL24+fjEivnhj+9dbaz9orX07It6JiM/f4ecAAAAwiGqIbBHx7zPzjcx8drPtwdba+5uPvxcRD24+\nfjgivnvje9/dbAOAaZyrcqgiCcDsqhPr/NnW2nuZ+Sci4pXM/L2bX2yttcw86oq4CaPPHnwhAFzQ\npQLe9ueYcAeA2ZQqka219zZ/fhARvxz321O/n5kPRURs/vxg8/L3IuLRG9/+yGbbx//OF1prT9wY\nYwkAAMDgDobIzPwjmfmT248j4i9FxO9ExMsR8czmZc9ExDc2H78cEU9n5qcy87MR8VhEvHbqHQeA\nS8vMW/+3j9ZWAGZTaWd9MCJ+eXMBvBcR/7K19m8z8zci4qXM/HJEfCcivhQR0Vp7MzNfioi3IuKH\nEfGV1tqPzrL3AAAAXFSO8AT02PGUAHAOvWviXccunvPvBoATeaMy3LA6sQ4ALNYlAt7279n1s0yy\nA8BM7rJOJAAAACsjRAIAAFAmRAIAAFBmTCQAq3WNyW6MjQRgdiqRAAAAlAmRAAAAlAmRAKzSvlbW\nzLxIO2nv54ywhjMA7CNEAgAAUCZEAsAV9aqRKpIAjEiIBAAAoMwSHwCsyjWW9QCAJREiAViFkcOj\ntSMBmIl2VgAAAMqESAAAAMqESAAAAMqESAAWb994yMwcaqxhb38s9wHAKIRIAAAAyoRIAAAAyoRI\nABartdZtZR1Vr6VVWysA1yZEAgAAUCZEAsCADk2yoyIJwLUIkQAAAJTdu/YOAMCp9ap0I4+FBIAZ\nCJEA3Mps7ZSzhsftfu/6fY/8/2DW3zcAh2lnBQAAoEyIBOAoJnW5jtkqe44TgOUSIgEAACgTIgEo\nm7WyNFsVb5/esh+jmvWYAWA/E+sAsFizBa6q3mQ7AHBuKpEAAACUqUQCcJB1F8c04u++txTJiPsL\nwPFUIgEAACgTIgEAACjTzgrAXtpYOVZv0h9trQDLoBIJAABAmRAJAJxcb01LS5MAzE2IBAAAoMyY\nSAB22lctMp4NANZNJRIAAIAyIRKAj2itqUJyMr1xkcZGAsxJiAQAAKBMiAQAAKDMxDoARER/2QVt\nrNzF9vjZdYxttznGAOahEgkAAECZEAkAAECZdlaAldPGyqVoawVYBpVIAAAAyoRIAOCietVGa0cC\njE+IBAAAoEyIBFgx4yEBgGOZWAdghYRHrs0kOwDzUokEAACgTIgEAACgTIgEAACgzJhIgBUxFpLR\nGBsJMB+VSAAAAMpUIgHOxKLpcBojvpdUR4E1EyIBTmzEG94eN8OMoNfWOiKttsCaaWcFAACgTIgE\nOKFZqihbqiiMZrZjcrb3PMApCJEAAACUGRMJcAKzVSNmq/awLsZHAoxNiAQ4MzeWcDsjv3esawms\nmXZWAAAAylQiAe6g126nGgEALJFKJAAAAGUqkQC3oAIJ69ab/MfYSGDpVCIBAAAoEyIBAAAoEyIB\njqSVFdjqvednWecS4FhCJAAAAGVCJAAAAGVCJAAAAGWW+AAoMhYS2MVyH8DaqEQCAJxAZu4Ni601\nE+0AiyFEAgAAUKadFeAAbawAAB9SiQQAAKBMiATo2FeF7I19Atbt0NhIgNkJkQAAAJQJkQAAAJQJ\nkQA79Kbj18YKVFjuA1gqIRIAAIAyIRIAAIAy60QC3GBNSOCUtueNXeeW7TbnFmA2KpEAAACUCZEA\nAGdm7UhgSYRIAAAAyoyJBNiwpAcAwGEqkQAAAJSpRAIXM+O4H1VI4JQysztT66icC4GbVCKBixj9\nBgmA/ZzDgZuESAAAAMq0swJnNevTa61bwLlszy+znR+3++v8CKhEAgAAUKYSCXCDJ+zApcxakQQQ\nIoGz6N0UCWoAHxr5nNibSXbk/QbOSzsrAAAAZUIkAAAAZUIkAAAAZUIkcHL7xkNmpjE0ABPpnbdN\nCATrJUQCANDVC5LCJKyPEAkAAECZJT6Ak7GsBwDA8qlEAgAAUKYSCdyZCiTA8m3P57vO+dttzvmw\nDiqRAAAAlAmRAAAAlAmRwJ1YExJgXawdCQiRAAAAlAmRAAAAlAmRAAAAlAmRwK201rrjIQFYtt64\nSGMjYdmESAAAbuXQJDvCJCyTEAkAAEDZvWvvADCX3lNlbawAAMunEgkAAECZSiRQogIJwD7b68Cu\na8V2m2sFLIdKJAAAAGVCJAAAAGVCJHCQVlYAKnrXBMt9wHIIkQAAAJQJkQAAAJQJkQAAAJRZ4gMG\nMdtYEWMhAdilstzHiFzXoE6IhCsb+YK6i4ssABW9MDki61lCnXZWAAAAylQi4UpmeTK75cksALeh\nIgnLoxIJAABAmRAJHORpLAB3Ndu1ZJbKKVyDdla4gn0XptkusABwjFGvcwIjHEclEgAAgDIhEi6o\ntaYKCQCD2XcN7l23Yc2ESAAAAMqESAAAAMpMrAMX0GuF0cYKANfXW8/S2pHwUSqRAAAAlAmRAAAA\nlAmRAAAAlAmRcGa9JT2MrQCAsfSuz5b7gPuESAAAAMqESAAA+JheNVJFkrWzxAeciWU9AABYIpVI\nAAAAylQi4cRUIAFgGbbX7V3X9u0213bWSCUSAACAMiESAACAMu2scCLaWAFgmbS1wkepRAIAAFAm\nRAIAAFAmRAIAAFAmRMIJGA8JAMvXu6b37gVgaUysA3cgPALAuphkB1QiAQAAOIIQCQAAQJkQCQAA\nQJkxkUxjlgHrxkEAwPJVxkaOxP0Jp6QSCQAAQJlKJMMb8WkeAMBMzBzLKQmRDGu28OikDADr02tr\nHZEwySmU21kz84HM/M3M/JXN55/JzFcy81ubPz9947XPZ+Y7mfl2Zn7hHDsOAADA5R0zJvKrEfG7\nNz5/LiJeba09FhGvbj6PzHw8Ip6OiM9FxJMR8bXMfOA0u8tazPI0b8vTPABYt9nuBWa712IspRCZ\nmY9ExF+NiH9yY/NTEfHi5uMXI+KLN7Z/vbX2g9batyPinYj4/Gl2FwAAgGuqjon8RxHx9yLiJ29s\ne7C19v7m4+9FxIObjx+OiP9843XvbrZ9RGY+GxHPHrW3LF7vqdhsT/gAgHUZ+V6ltxTJyPvNmA5W\nIjPzr0XEB621N/a9pt0/Ao+qibfWXmitPdFae+KY7wMAAOB6KpXIn4+Iv56ZfyUifiIi/lhm/vOI\n+H5mPtRaez8zH4qIDzavfy8iHr3x/Y9stgEAADC5g5XI1trzrbVHWms/E/cnzPkPrbW/GREvR8Qz\nm5c9ExHf2Hz8ckQ8nZmfyszPRsRjEfHayfecRWmt7W1lzUxtFgAAd9C7n+rdh8Eud1kn8hcj4qXM\n/HJEfCcivhQR0Vp7MzNfioi3IuKHEfGV1tqP7rynAAAAXF2O8NQhM6+/E1yVCXUAAC7HvRd7vFGZ\ns+YulUi4MycwAACYS2mdSAAAAIgQIrkiVUgAgOvo3WuNMNyNsQmRAAAAlAmRAAAAlAmRAAAAlJmd\nlYszFhIA4Pq291277s2229ybsYtKJAAAAGVCJAAAAGXaWbkYbawAAOPR1sqxVCIBAAAoEyIBAIDI\nzL0Vx15HGeujnZWL2Hfi0RoBAABzUYkEAACgTIjkrFprqpAAABPptbRqayVCiAQAAOAIQiQAAABl\nJtbhLKwJCQAwL2tH0qMSCQAAQJkQCQAA7GTtSHYRIgEAACgzJnJisz390TcPALAso96Puu88L5VI\nAAAAyoTICc240KunQQAA85rtXm7G++WZCJGTme3N0BuMDQDAPGa8r5vt3nkWQiQAAABlJtbhLGZ7\nSgUAQM32Pm+WKt92P92fno5KJAAAAGUqkZPoPenxVAUAgEsb8R50luro7ITIwQmPAABQ02u11dZ6\nOtpZAQAAKBMiAQAAKBMiAQAAKBMiB7ZvPOSMC70CAMCl9O6XTb5zd0IkAAAAZUIkAAAAZULkgFpr\n3VZWAADgsF5Lq7bW2xMiAQAAKBMiAQCAxTo0yY6K5PHuXXsHAAAAzqUXEg0Vux2VSAAAAMqESAAA\nAMqESAAAAMqEyAEdGvwLAAAc1ls2z3jI2xMiAQAAKBMiB2YqYgAAYDSW+AAAABbFsh7npRIJAABA\nmUrk4LZPSnY9Tdlu8zQFAABUIC9FJRIAAIAyIRIAAIAy7ayT0NYKAAC7aWO9LJVIAAAAyoRIAAAA\nyoRIAAAAyoTIyfR6unu94AAAsETGQ16eEDmhzNz7hmitCZMAAMDZCJEAAACUWeIDAACYjjbW61GJ\nBAAAoEwlcmLbJyy7nsJst3kKAwDAkqhAXp9KJAAAAGVCJAAAAGVC5AIcWvIDAACWYN+9be9+mNMT\nIgEAACgTIgEAACgTIgEAACgTIhekNy7S2EgAAGbVu581FvLyrBO5MNaOBABgKawJOSaVSAAAAMqE\nSAAAAMqESAAAAMqEyIXqLbhqkh0AAOC2TKyzUAYhAwAA56ASCQAAQJkQCQAADOnQEC3DtK5DiAQA\nAKDMmMiFMRYSAAA4J5VIAAAAyoTIBdlXhez1kgMAwOgsXzcWIRIAAIAyIRIAAIAyIXIBetMba2MF\nAGApLPcxBiESAACAMiESAACAMutETsyakAAArM32PnfXvfB2m3vh81KJBAAAoEyIBAAApnNo7UgT\n7ZyPEAkAAECZMZETMhYSAAC4FpXIyQiQAADwod49sJbW8xAiAQAAKBMiAQAAKBMiAQAAKDOxziSM\nhQQAgN2298O77pm329wzn45KJAAAAGVCJAAAAGXaWQenjRUAAGq0tV6GSiQAAABlQiQAALAombm3\n4tjr9KNGiAQAAKDMmMiB7XtKoo8bAAC4FpXIAbXWBEgAALijXkurttbbEyIBAAAoEyIBAAAoEyIB\nAAAoEyIBAAAoEyIBAAAoEyIBAAAoEyIHlJnd6YgBAIDDesvmWTrv9oRIAAAAyoRIAAAAyoRIAAAA\nyoTIgfXGRRobCQAAu/Xul42FvDshcnCHJtkRJgEAgEsSIgEAACi7d+0dAAAAOIVel5421tNRiQQA\nAKBMJXIS2ycnu56ubLd5ugIAwBqpQF6WSiQAAABlQiQAAABlQuRkeuV4y30AALA2WlkvT4gEAACg\nTIgEAACgzOysEzJTK8BxltLu79wO8CFtrNejEgkAAECZSuTEVCQB+pZSgdxybgdgBCqRAAAAlKlE\nArA4S6tAfpyKJLBmxkJen0okAAAAZSqRC3BobKQnMsBaLL0C+XEqksDa7DvPOw9elhAJwCos4QZj\n382TB4YAXJJ2VgAAAMqEyAXZ9xS6tba6Fi9gfUy0ALBcvftZ5/jLEyIBAAAoMyYSgKmtqQJ5aCK1\nm68BgHMRIhfGDQYAAEuxpgeFM9HOCgAAQJlK5EJZOxJYujU/ndZ1AsA1qUQCAABQJkQulGmQAZYv\nM7vLOwHAOWhnBWA6HpIBwPWoRAIAAFCmErkwa55oAlg+57jdMtMkO8AimUhsTCqRAAAAlAmRAAAA\nlAmRAAAAlBkTuSBmKwSWyljIw4wbApbs0DnO+e2yVCIXwJqQAACsQW9tXOvjXo4QCQAAQJkQCcDQ\nep0Wui0+qfd78ZQegFMQIgEAACgzsc7ETDQBwD6ZaZIdYJFMJHZ9QiQAQ/KgDADGpJ0VAACAMpXI\nCXk6DyyZc9zpaPkClsw57npUIgEAACgTIgEAgGn1qo2WNjoP7ayT0eYFnMJsF1Xnt7uptHyNyP93\ngDGpRAIAAFCmEjkJFUjgFEauOnF++9aOHJWJMYAqk+xc1sFKZGb+RGa+lpm/lZlvZuYvbLZ/JjNf\nycxvbf789I3veT4z38nMtzPzC+f8BwAAAHA5lXbWH0TEX2it/emI+NmIeDIzfy4inouIV1trj0XE\nq5vPIzMfj4inI+JzEfFkRHwtMx84x84DUDdTBeomT45PKzOn+53OeuwCLNXBENnu+9+bT39881+L\niKci4sXN9hcj4oubj5+KiK+31n7QWvt2RLwTEZ8/6V4DAABwFaWJdTLzgcz8ZkR8EBGvtNZ+PSIe\nbK29v3nJ9yLiwc3HD0fEd298+7ubbR//O5/NzNcz8/Vb7/0KtNb2PoGd8WkyQJVz3Hn5/QJL1Du3\n9e6rOU5pYp3W2o8i4mcz849HxC9n5p/62NdbZh71f6S19kJEvBARcez3AlBnYi56RjwGTIwBMLaj\nlvhorf2PiPi1uD/W8fuZ+VBExObPDzYvey8iHr3xbY9stgEAADC5yuysP72pQEZm/uGI+IsR8XsR\n8XJEPLN52TMR8Y3Nxy9HxNOZ+anM/GxEPBYRr516x9dAGysAABzvUFsrd1NpZ30oIl7czLD6YxHx\nUmvtVzLzP0XES5n55Yj4TkR8KSKitfZmZr4UEW9FxA8j4iubdlgAAAAmlyMkcWMid+tVIgEOMRaS\n2TmGgVNwLjnKG621Jw69qDSxDpflQAfuykMolmB7vO6baMfxDHAdR02sAwAAwLoJkQAAwGJZO/L0\nhEgAAADKjIkEWBBjqlmizNw7LnL7dYB9XBtPTyVyQNa1AQAARiVEAgAAUKadFWABtOqwdIeW+7j5\nGoAtS16dh0okAAAAZULkwExFDAAfZap+oKJ3PlCFvDshEgAAgDJjIgEmZiwkAHBpKpGD07YD7CNA\nsma9Y9y1EdbrUBur6+NpCJEAAACUCZEAAACUCZEAAACUmVhnEhZZBraMhYT7XBuBLdfGy1KJnIyJ\nBAAAgGsSIgEAACjTzgowCa06sJu2Vlg318fLU4kEAACgTCVyQp64AsAnHbo+ujbCsqhAXo8QCTCB\nfRdKF0kA4NK0swIAAFCmEjkxba2wfFp14HiZ6doIC+baeH0qkQAAAJQJkQAAAJQJkQAAAJQZE7kA\npjSH5THeA+7GvAGwTGYrH4NK5ILse/O01ro3pAAAAFVCJAAAAGXaWQEGo1UHTseQD1gGwzzGohIJ\nAABAmUrkwphIAAA+KTNdG2FCKpBjEiIBBuFCCQDMQDsrAAAAZULkQmVmd8kPAACA2xAiAQAAKDMm\ncqGMrQIAAM5BiFwY4RHmZT07OB/rr8KcrDwwJu2sAAAAlAmRAAAAlAmRAAAAlAmRC9Ib76FXHObR\nW57HEj1wnN77xrUR5mH5urEIkQAAAJQJkQADOvTE1VNXAOBaLPGxAJb1AICPcm2EZcpMy30MQCUS\nAACAMiESYGAmEoDjmWgOls2Qj+sTIgEAACgzJnJixnsAAACXJkROSHiE9TGRABzm+gjrsn1fuz5e\nnnZWAAAAyoRIAAAAyoRIAAAAyoTIyRjvAetlSnPYrXf8W9YDlq/3HndtPA8hEmAywiQAcE1CJAAA\nAGWW+JiENlYA+CjXRmDLch+XpRIJAABAmRAJMCkTCbBmqpDALuYNuAwhEgAAgDJjIgfnSSsAADAS\nIXJgvTWvACJMJMD6eLgKVBy6Pjpf3I12VgAAAMqESAAAAMqESAAAAMqMiRyQ8R7AsYyNZOlcG4Hb\nyEzXxjNQiQRYkEPrYwEA3JUQCQAAQJkQCQAAQJkQCQAAQJmJdQAWxOQjAMC5CZEDOjTLohtB4OOE\nR5bODMTAbey7Pjpf3I12VgAAAMqESAAAAMqESAAAAMqMiRxYZhr7AXQZC8naGBsJVLg+npdK5OAy\nc++B3lrrvkEAAABOTYgEAACgTIgEmJRWHdasd4zr0oH16nXq9Tr8OI4QCQAAQJmJdSZhIgFgSwUS\n7nNtBLZcGy9LJRIAAIAyIRIAAIAyIXIyJhKA9TJZAOxmOSxYN62slydEAgAAUCZEAgAAUCZEAgAA\nUGaJjwmZ0hzWpzcWErjv0PXR+wWWxVjI6xEiJyZMAsAnZaZrIyyY8Hh92lkBAAAoU4kEGJinrQDA\naFQiAQAAKFOJXAATCcDyqEDC3Zg3AJbJRHNjUIkEAACgTIgEAACgTIhckH1l/NZatzUOGEuvVUe7\nDhyn975xbYR59O5nXRsvT4gEAACgTIhcmENPXD11BQAA7kKIBAAAoMwSHwCDsKwHnE9mWu4DJuTa\nOCaVyIUykQAAADMz0dy4hEgAAADKhMiFMg0yAAAzM1nkuIRIAAAAyoRIgEEYywwAzMDsrAtjBisA\n+CRDPGBO2/eo2ZXHohIJAABAmRAJMBgTCcDpmGgOlsGQj7EIkQAAAJQZE7kgnrQCAADnJkQugMl0\nYHlMJAB349oIy5SZro0D0M4KAABAmRAJAABAmRAJAABAmTGREzPeA5bv0NhI73X4JBPNwbKZN+D6\nVCIBAAAoEyIBJtBbYNkiywDAJWlnnZA2VgD4JNdHWBdtrdejEgkAAECZEDkZT1lhvTJTWyvs0Dv+\ne+8bYBl673HXxvMQIgEAACgTIgEAACgzsc4ktLECWyYSgPtcG4Et18bLUokEAACgTIgEAACgTIgE\nAACgTIgcnGnLgX1Mac6aGQ8J7GI5rMswsQ7AxEwkwNoIj0DFoeuj88XdqEQCAABQphI5sF4bKwAA\nwDWoRAIAAFCmEjkg4z2AYxkbydK5NgK3kZmujWegEgkAAECZEAkAAECZEAmwIIfWx4IZWS8ZuAtr\nR56eEAkAAECZiXUAFshEAgBwn4m5Tk8lEgAAgDKVyAEdmqrfExMA1kD1ALir3phqbk+IBFgoa0cy\nK+ERYGzaWQEAACgTIgdmKmIAADhe735ZR8PdCZEAAACUGRMJsHCVsZEj8qT4vEb+f7+L4wFgHELk\n4EyMAayVc9x5zBYeAY5hYq7L0M4KAABAmRA5icw00Q5wJ7M+gXV+O51Zf5ezHrsASyVEAgAAUGZM\nJMCK9MZZw2hUIIEqYyEvSyUSAACAMpXIyWTm3ictrTVPWoCSkc8VZqM+PU/ogSVzjrs8lUgAAADK\nhEgAAADKtLNOqDcxhpYvYHbOcaejxQtYMue461GJBAAAoEyIBGBImbn3SbIlSg7b9zvq/V4BoEI7\n68S0fAEAsDbaWK9PJRIAAIAyIXIBtHwBS9Y7vznHfVLv9+IJPQCncDBEZuajmflrmflWZr6ZmV/d\nbP9MZr6Smd/a/PnpG9/zfGa+k5lvZ+YXzvkPAAAA4HIqlcgfRsTfba09HhE/FxFfyczHI+K5iHi1\ntfZYRLy6+Tw2X3s6Ij4XEU9GxNcy84Fz7DwAALAeJg0bw8EQ2Vp7v7X2XzYf/6+I+N2IeDginoqI\nFzcvezEivrj5+KmI+Hpr7QettW9HxDsR8flT7zifpOULWKpDbfvOcYfbWN1cATPTqj+Wo8ZEZubP\nRMSfiYhfj4gHW2vvb770vYh4cPPxwxHx3Rvf9u5m28f/rmcz8/XMfP3IfQYAAOBKyiEyM/9oRPzr\niPg7rbX/efNr7f5jgaMeA7fWXmitPdFae+KY7wNgvUwkBgDXVwqRmfnjcT9A/ovW2r/ZbP5+Zj60\n+fpDEfHBZvt7EfHojW9/ZLMNAACAyVVmZ82I+KcR8buttX9440svR8Qzm4+fiYhv3Nj+dGZ+KjM/\nGxGPRcRrp9tleowbAlgfYyGBpTLee0z3Cq/5+Yj4WxHx25n5zc22vx8RvxgRL2XmlyPiOxHxpYiI\n1tqbmflSRLwV92d2/Upr7Ucn33MAAAAuLkeoTGXm9XdiYXr/Xz2xAZbAec7vAFg+57mLe6MyZ81R\ns7MCAACwbkIkAAAAZZUxkQAwnG0b065Wp+22pbY6ae8C4JpUIgEAAChTiVwYT6cBAFiKNXedjEyI\nBGBqa7rB8KAQgBFoZwUAAKBMiFyQfU+oM9MTamDxeue5EdZEBuD2evezzvGXJ0QCAABQZkwkAKuw\n5CfVuk2AtcjMVYyBH50QuQAmWgC4rzfJzhI5xwNwDdpZAQAAKFOJnJgKJMBuS69IOscDa7ampZ1G\npRIJAABAmUokAIu1tIqkJ+sAjEAlEgAAgDKVyAkZCwlwHOdGgOUxNvJ6VCIBAAAoEyIBAAAoEyIn\no5UVAAA+1LsHXsrEaqMRIgEAACgTIgEAgKll5t6KZGtNRfLEzM46CW2sAADACFQiAQAAKFOJHJwK\nJAAA1Fg78jJUIgEAACgTIgEAACgTIge2r5W1N/sUAACs3aHZWrkbIRIAAIAyIRIAAIAyIRIAAIAy\nIXJArbXueEgAAOCw3rhIYyNvT4gEAACgTIgEAAAW69BMrSqSx7t37R0AAAA4l15INFTsdlQiAQAA\nKBMiAQAAKBMiAQAAKBMiB3Ro8C8AAHBYb9k84yFvT4gcmFmkAACA0QiRAAAAlFniAwAAWBTLepyX\nSiQAAABlKpGD2z4p2fU0ZbvN0xQAAFCBvBSVSAAAAMqESAAAAMq0s05CWysAAOymjfWyVCIBAAAo\nEyIBAAAoEyIBAAAoEyIn0+vp7vWCAwDAEhkPeXlC5IQyc+8borUmTAIAAGcjRAIAAFBmiQ8AAGA6\n2livRyUSAACAMpXIiW2fsOx6CrPd5ikMAABLogJ5fSqRAAAAlAmRAAAAlAmRC3BoyQ8AAFiCffe2\nvfthTk+kGempAAARM0lEQVSIBAAAoEyIBAAAoEyIXJBeS6u2VgAAZtW7n9XGenlCJAAAAGXWiVwY\na0cCALAU1oQck0okAAAAZUIkAAAAZUIkAAAAZULkQvUWXDVTKwAAcFsm1lkog5ABAIBzUIkEAACg\nTIgEAACGdGiIlmFa1yFEAgAAUGZM5MIYCwkAAJyTSuSC7AuQvTYAAAAYnZUHxiJEAgAAUCZEAgAA\nUCZEAgAAUCZELkBvemNjIQEAWArLfYxBiAQAAKBMiAQAAKDMOpETsyYkAABrs73P3XUvvN3mXvi8\nVCIBAAAoEyIBAIDpZKaJdq5EiAQAAKDMmMgJGQsJAABci0rkZARIAAD4UO8eWEvreQiRAAAAlAmR\nAAAAlAmRAAAAlJlYZxLGQgIAwG7b++Fd98zbbe6ZT0clEgAAgDIhEgAAgDLtrIPTxgoAADXaWi9D\nJRIAAIAyIRIAAFiUzNxbcex1+lGjnXVg+w5wJXgAAOBaVCIBAAAoEyIH1FpThQQAgDvqtbRqa709\nIRIAAIAyIRIAAIAyIRIAAIAyIRIAAIAyIRIAAIAyIRIAAIAyIXJAmdmdjhgAADist2yepfNuT4gE\nAACgTIgEAACgTIgcWK+lVVsrAADs1rtf1sZ6d0IkAAAAZULk4A5NsqMiCQAAXJIQCQAAQNm9a+8A\nAADAKfS69IyFPB0hchLbg37XG2PkllZvVgCAZRr5HvQm96Onp50VAACAMpVIzmr7hMoTIACAZZil\nAsn5qEQCAABQJkROZtaKnidWAABzm3F5uVnvnUcnRAIAAFBmTOSEejO1AgDAqc1236kCeV5C5MRG\nfnP0liIZeb8BAPiQdRfZRTsrAAAAZUIkAAAAZUIkAAAAZcZEcha9yX+MjQQAGJuxkPSoRAIAAFAm\nRAIAAFAmRHJWmbm35WG29YYAANZg3z1a776OdREiAQAAKBMiuYheNVJFEgAA5mF2VgAAwIyslKlE\nAgAAUKYSycVYOxIAYDwqkBxLJRIAAIAyIRIAAIAy7axcnLZWAIDr08bKbalEAgAAUCZEAgAAUCZE\nAgAAUCZEcjW9Xvtejz4AAHdjPCR3IUQCAABQZnZWrspMrQAAl6MCySmoRAIAAFAmRAIAAFCmnZUh\naGsFADgfbayckkokAAAAZSqRTGPEZT88uQMAPm7EexY4JSGSofTaWkek1RYA2Jrl/mXL/Qu3pZ0V\nAACAMiGSIc32ZGy2J48AwGnNdi8w270WYxEiAQAAKDMmkmEZHwkAjG6W+5Qt9ymcghDJ8EY82c12\nwQAALm/Eexg4Be2sAAAAlKlEwi30Wm21tQLA8vW6ktwDsHQqkQAAAJQJkXAHmbn3aWNrzdhJAAAW\nR4gEAACgzJhIAAAoMhYSVCLhJHoXDS2tALAMAiTcJ0QCAABQJkQCAABQdjBEZuYvZeYHmfk7N7Z9\nJjNfycxvbf789I2vPZ+Z72Tm25n5hXPtOAAAAJdXqUT+s4h48mPbnouIV1trj0XEq5vPIzMfj4in\nI+Jzm+/5WmY+cLK9hYFZ7gMAlql3He9d/2GpDobI1tp/jIj//rHNT0XEi5uPX4yIL97Y/vXW2g9a\na9+OiHci4vMn2lcAAACu7LZjIh9srb2/+fh7EfHg5uOHI+K7N1737mbbJ2Tms5n5ema+fst9AAAA\n4MLuvE5ka61l5tF9eq21FyLihYiI23w/jGrb0rKr7WW7TdsLAIzPkh6w220rkd/PzIciIjZ/frDZ\n/l5EPHrjdY9stgEAALAAtw2RL0fEM5uPn4mIb9zY/nRmfiozPxsRj0XEa3fbRZjToYl2AABgRgfb\nWTPzX0XEn4uIn8rMdyPiH0TEL0bES5n55Yj4TkR8KSKitfZmZr4UEW9FxA8j4iuttR+dad8BAOAs\nerOxwtrlCBURYyJZMhchAJiP6zcr9UZr7YlDL7ptOytQZO1IAJjHoTUhASESAACAIwiRAAAAlAmR\nAAAAlB2cnRW4u+0Yil1jLLbbjLMAgOvpzVPgGg0fpRIJF2TtSAAAZidEAgAAUKadFQCAVbOkBxxH\nJRIAAIAylUi4gszsTrIzIk9jAbirka9zH+e6B/upRAIlrbWpLv4AjMU1BJZDiAQAAKBMOytcSW/t\nyJFZ1xKAY8x2nXN9g8NUIgEAAChTiYQrm7UiCQA9s13XVCChToiEQYx88erNJDvyfgNweb3w6JoB\ny6CdFQAAgDIhEgAAgDIhEgAAgDIhEjioN4ZltokTADgf4yFhHYRIAAAAyoRIAAAAyizxAZT01rO0\n3AfAumljhXVRiQQAAKBMiASOkpl7nyq31ky0AwCwcEIkAAAAZcZEAgBwK8ZCwjqpRAK3cqitFYBl\n23eu710fgGUQIgEAACgTIgEAACgTIgEAACgTIoE7sdwHwLr0zu/GQsI6CJEAAACUCZEAAACUWScS\nuLNt+9Ku9qbtNi1OAHOzJiSwpRIJAABAmRAJnExvgWmT7ADMqzeRjiokrI8QCQAAQJkQCQAAQJkQ\nCZyctSMBlsGakMAuQiQAAABlQiQAAABlQiQAAABl9669A8AybcfK7BpLM/K4SGN8gEsb+Zy4i/Mk\nIEQCZ9ULkyPa7qebJODcZjkvAnycdlYAAADKVCKBi8jMqZ66q0gC5zLTufAm50NgSyUSAACAMiES\nuJjM9CQbYELO3cBN2lmBixv1ZmTfTLKj7i8wn32trM4zwExUIgEAACgTIgE29lUCWmvTToQBAHBq\nQiQAAABlxkQCAJxZr5vBeEhgNiqRAAAAlKlEAtywrQjsm6n15msADlGBBJZIJRIAAIAyIRIAAIAy\nIRJgh8zsLvkBcMi+c0Xv/AIwAyESAACAMiESAACAMiESoKPX0qqtFdild37QxgosgRAJAABAmRAJ\ncMChSXZUJAGANREiAQAAKLt37R0AAFiCXleCsZDAkgiRAEXbm8BdN4rbbW4UYX2ER2BttLMCAABQ\nJkQCAABQJkQCAABQJkQCHKk3xslyH7AuxkMCayREAgAAUCZEAtxCZu6tMrTWVCQBgMWyxAcAwJG0\nsQJrphIJAABAmUokwB1sKw67qhLbbaoSsBwqkAAqkQAAABxBJRLgzEacZEfFhBmM+N4BQIgEOIle\nW+uItNoyslneR1veR8DaaGcFAACgTIgEOKHZKhKzVXxYvtmOydne8wCnIEQCAABQZkwkwIkZHwnH\nm+X9suX9AqyZEAlwJiPeZM52ow4RY76XANZMOysAAABlKpEAK9JrtdXWyjX0quOORYAxqUQCAABQ\nJkQCAABQJkQCAABQZkwkwAoZG8m1GQsJMC+VSIAV692sWw4EANhFiAQAAKBMiAQALkorK8DchEgA\nAADKTKwDsHIm2eFSVCABlkElEgAAgDIhEgAAgDLtrABEhLZWzkcbK8CyqEQCAABQJkQC8BGZubc6\n1KsoAQDroJ0VADibfQ8etLECzEslEgAAgDIhEoCdei2t2lo5pHecqEICzE2IBAAAoEyIBAAAoEyI\nBAAAoMzsrADstR27tmts23ab8W3c1Bsv61gBWAaVSAAAAMqESAAAAMq0swJwUKWtdURLbZ8c+Xe+\ny1L/PwCslUokAAAAZUIkAGWzVZR6C97PaGn/HgDmJEQCcJTMnDJMzm7Wf8NsxwoAhwmRAAAAlJlY\nB4Bb6U22AyqQAMulEgkAAECZSiQAdzJixam3FMmI+9vTq/TO9m8BYBlUIgFYnN7kP2Y4BYC7ESIB\nAAAo084KAAPSxgrAqFQiAQAAKBMiAVisQ2MjR7Vv33r/HgC4FCESAACAMiESAACAMiESgMWbZbmP\n3v5oYwVgFEIkAAAAZUIkAAAAZUIkAAAAZfeuvQMAcAnbMYW7xhxut11r3GFvXKaxkACMRiUSgFUZ\nbe1Ia0ICMBshEgAAgDIhEgAAgDIhEgAAgDIhEoBV6o2LvMTYyN7PMRYSgJEJkQAAAJQJkQAAAJRZ\nJxKA1brG2pHWhARgdiqRAAAAlAmRAKxeZl59oh0AmIUQCQAAQJkxkQBwAcZCArAUKpEAsNELc3dp\naRUgAVgSIRIAAIAy7awAUGSCHQBQiQQAAOAIKpEAcMN2jOK5q47GQgIwKyESAHY4V5gUHgGYnXZW\nAAAAylQiAaDjVBVJFUgAlkIlEgAAgDKVSAAoUEkEgPtUIgEAACgTIgEAACgTIgEAACgTIgEAACgT\nIgEAACgTIgEAACgTIgEAACgTIgEAACgTIgEAACgTIgEAACg7W4jMzCcz8+3MfCcznzvXzwEAAOBy\nzhIiM/OBiPjHEfGXI+LxiPgbmfn4OX4WAAAAl3OuSuTnI+Kd1tp/ba3934j4ekQ8daafBQAAwIWc\nK0Q+HBHfvfH5u5tt/19mPpuZr2fm62faBwAAAE7s3rV+cGvthYh4ISIiM/9bRPyfiPiDa+0PFP1U\nOE4Zn+OUGThOmYHjlBmc8jj9k5UXnStEvhcRj974/JHNtp1aaz+dma+31p440/7ASThOmYHjlBk4\nTpmB45QZXOM4PVc7629ExGOZ+dnM/EMR8XREvHymnwUAAMCFnKUS2Vr7YWb+7Yj4dxHxQET8Umvt\nzXP8LAAAAC7nbGMiW2u/GhG/esS3vHCufYETcpwyA8cpM3CcMgPHKTO4+HGarbVL/0wAAAAmda4x\nkQAAACzQECEyM5/MzLcz853MfO7a+8N6ZeYvZeYHmfk7N7Z9JjNfycxvbf789I2vPb85bt/OzC9c\nZ69Zk8x8NDN/LTPfysw3M/Orm+2OU4aRmT+Rma9l5m9tjtNf2Gx3nDKczHwgM38zM39l87njlKFk\n5u9n5m9n5jcz8/XNtqsep1cPkZn5QMT/a+/+Qa4s4zCOf6/CIipoqEQ00MHFgmwJwQYTCqPIpjAo\nHIIWg4IgqqXVKVpqqkjojwglSUNhFrREiRGU6SD9IcV6h4hqMayr4TzFQRrO8vI89n4/cDj3/Xve\nF+7h4nB+PPf9HF4A7gI2AQ8k2TTuqrSCvQrsuKD2FHCk7UbgyDBnyOku4Mbhf14c8iwtp/PAE203\nAVuAPUMWzamm5Bywve3NwGZgR5ItmFNN02PAibm5OdUU3d5289xPeYya09GbSOBW4FTbb9r+AewH\ndo68Jq1QbT8Gfr6gvBPYN4z3AffN1fe3Pdf2W+AUszxLy6bt2bafD+PfmH3xWYs51YR05vdhump4\nFXOqiUmyDrgbeGmubE51MRg1p1NoItcCP8zNTw81aSpWtz07jH8EVg9js6tRJVkP3AJ8ijnVxAxb\nBL8AloDDbc2ppuh54Engr7maOdXUFPggybEkjwy1UXO6bD/xIf0ftW0SH2ms0SW5CngLeLztr0n+\nvWZONQVt/wQ2J7kGOJjkpguum1ONKsk9wFLbY0m2/dffmFNNxG1tzyS5Hjic5OT8xTFyOoU7kWeA\nG+bm64aaNBU/JVkDMLwvDXWzq1EkWcWsgXy97dtD2Zxqktr+AnzE7GyOOdWUbAXuTfIds+NU25O8\nhjnVxLQ9M7wvAQeZbU8dNadTaCKPAhuTbEhyGbODoIdGXpM07xCwexjvBt6Zq+9KcnmSDcBG4LMR\n1qcVJLNbji8DJ9o+N3fJnGoyklw33IEkyRXAHcBJzKkmpO3Tbde1Xc/s++eHbR/EnGpCklyZ5Op/\nxsCdwFeMnNPRt7O2PZ/kUeB94FLglbbHR16WVqgkbwLbgGuTnAaeBfYCB5I8DHwP3A/Q9niSA8DX\nzJ6YuWfYviUtp63AQ8CXw3kzgGcwp5qWNcC+4YmAlwAH2r6b5BPMqabPz1NNyWpmRwJg1ru90fa9\nJEcZMadp3eYtSZIkSVrMFLazSpIkSZIuEjaRkiRJkqSF2URKkiRJkhZmEylJkiRJWphNpCRJkiRp\nYTaRkiRJkqSF2URKkiRJkhZmEylJkiRJWtjff/nd9olHr6oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11ac8bdd0>"
      ]
     },
     "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": 163,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The LDZ accomomdates a translation +/-0.2% of D=512 pixels = +/-1.02 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": 164,
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}