anim_interferometer_dance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Animation of the PSF of dancing interferometers\n",
    "\n",
    "_M Kenworthy_ Leiden Observatory kenworthy@strw.leidenuniv.nl\n",
    "\n",
    "Based on hcipy tutorials https://docs.hcipy.org/0.3.1/tutorials/index.html and examples, and http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pip install progressbar2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pip install ffpmeg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from hcipy import *\n",
    "from scenes import *\n",
    "\n",
    "from progressbar import progressbar \n",
    "\n",
    "anim_version='v1'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import animation, rc\n",
    "#from IPython.display import HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAF1CAYAAADoc51vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABQCklEQVR4nO3dd3gU1f7H8ffZdCAEAoEUegkl1BAQUIoiPQqIKIgCCooF9IoXlQvca2+oF3+IFCmCYEFEQEBBRWzUEIr0DoYkhBAINXXP748s3ICBhGRmZzf5vp4nT7K7s3M+s0n2uzNzzhmltUYIIYSwWR1ACCGEa5CCIIQQApCCIIQQwkEKghBCCEAKghBCCAcpCEIIIQApCMLNKKVqKKW0UsrT6ixCFDdSEIRllFLnc33ZlVKXct0eaHU+IUoa+ZQlLKO1LnP5Z6XUEWCY1vpH6xK5JqWUp9Y6K7/7hCgq2UMQLkcpZVNKvaiUOqiUOqWUWqCUCrzOsgFKqZlKqQSl1HGl1GtKKQ/HY3WUUr8opVKVUslKqS9zPS9CKfWDUipFKXVCKfUvx/0+SqmJSql4x9dEpZSP47GOSqk4pdS/HOs7kntPxvHcd5VSxxzrnKqU8rtO7tpKqdWO7UtWSs1XSpXL9fgRpdQLSqntwAXHtmil1FCl1DFgtWO5r5RSiY5t/FUpFeG4v6Ujg2eudfZVSm0t7O9FFH9SEIQrehroDXQAQoHTwOTrLDsHyALqAM2BLsAwx2OvAquA8kAVYBKAUsof+BH43rH+OsBPjueMBVoDzYCmQCtgXK72goGKQBgwGJiulKrneOxtINzx3DqOZf59ndwKeNPRfgOgKvDSNcsMAHoC5RzbCDmvSQOgq+P2d0BdoBIQC8wH0FpvAk4BnXOt70Hg0+vkEQK01vIlX5Z/AUeAOx0/7wY65XosBMgk5xBnDUA7fq4MpAN+uZYdAPzs+HkuMB2ock1bA4At18lxEOiR63ZX4Ijj547kvDGXzvX4AmA8OW/wF4DauR5rAxwu4Pb3zp3J8Xo8kuv25e2udYN1lHMsE+C4/QIw3/FzIHARCLH6dy1frvsl5xCEK6oOfKOUsue6L5ucAnDtcl5AglLq8n024C/Hz8+Ts5ewUSl1GnhPaz2LnE/jB6/TdihwNNfto477Ljuttb6Qx+NBQClgc64sCvDIqxGlVCXg/4B2gL8j9+lrFvvr2uflvs9xaOx1oJ+j/cuvV0UgFZgH7FZKlQHuA37TWifklUcIkENGwjX9BXTXWpfL9eWrtT6ex3LpQMVcy5XVWkcAaK0TtdaPaq1DgeHAR0qpOo7n1b5O2/HkFJrLqjnuu6y8Uqp0Ho8nA5eAiFxZAnSuE+fXeJOcT/NNtNZlyTmco65ZJq+piHPf9wDQC7gTCCBnL4LL63G8XuuAPsBDyOEikQ8pCMIVTQVeV0pVB1BKBSmlel27kOPT7irgPaVUWcfJ6NpKqQ6O5/VTSlVxLH6anDfTbGAZEKyU+ofjRLC/UuoWx3KfA+McbVYk5xzAvGuaflkp5a2UagdEA19pre3Ax8B/HZ/+UUqFKaW6kjd/4DxwRikVBoy+6VcpZx3p5JwrKAW8kccyc8nZU2oMfFOINkQJIgVBuKIPgKXAKqXUOWA9cMt1lh0EeAO7yHnTX0jOOQeAlsAGpdR5x/qe0Vof1lqfI+dk611AIrAfuN3xnNeAGGA78Cc5J2pfy9VeoqOdeHJO4D6utd7jeOwF4ACwXil1lpwT1/XI28tAJDmHdpYDi278kuRpLjmHrI47tn99Hst8g+MQ3DWHuoT4G6W1XCBHiIJQSnUE5mmtq+SzqEtRSh0EhmsZ4yHyIXsIQhRjSqm+5BwqW211FuH6DCkISqluSqm9SqkDSqkXr7NMR6XUVqXUTqXUL0a0K4S4PqXUGmAK8JTjHIcQN1TkQ0aOrm/7yDkmGwdsAgZorXflWqYcsBboprU+ppSqpLVOKlLDQgghDGXEHkIr4IDW+pDWOgP4gpyucLk9ACzSWh8DkGIghBCux4iCEMbVA2jiHPflFk5O/+01SqnNSqlBBrQrhBDCQEaMVL52MA38fUCNJ9AC6AT4AeuUUuu11vv+tjKlHgMeAyhdunSL+vXrGxBRCCFKhs2bNydrrYMK81wjCkIcOVMBXFaFq0d2Xl4m2dEP+oJS6ldyJg77W0HQWk8nZ/4ZoqKidExMjAERhRCiZFBKHc1/qbwZcchoE1BXKVVTKeUN9CdnEFBuS4B2SilPpVQpcgYZ7TagbSGEEAYp8h6C1jpLKTUCWEnORF6ztNY7lVKPOx6fqrXerZT6npzRn3ZghtZ6R1HbFkIIYRyXHqksh4yEEOLmKKU2a62jCvNcGakshBACkIIghBDCQQqCEEIIQAqCEEIIBykIQgghACkIQgghHKQgCCGEAKQgCCGEcJCCIIQQApCCIIQQwkEKghBCCEAKghBCCAcpCEIIIQApCEIIIRykIAghhACkIAghhHCQgiCEEAKQgiCEEMJBCoIQQghACoIQQggHKQhCCCEA8LQ6gBBW01qjtWbdunUkJCSQmJhIUFAQ999/Px9++CFbt24lKyuLrKws5s2bx8qVK/n444/x9PTE09OT5557jpo1a/Liiy/i6elJuXLl6NChA507d2bnzp2UKlWKkJAQfH19rd5UIW5ICoIo1pKTk9m+fTvx8fEkJCQwbNgwjh49yrPPPktCQgIJCQm89tprjBgxgjFjxlCxYkWCg4Np3bo1ALVr18bX1/fKm7/WmvDwcPr373+lSAQFBeHl5UWzZs3IzMzk9OnTnD17FoBXXnmF9evXk5iYSP369dm2bRuzZs3ip59+Ijg4mJCQEIYOHYqPjw8XLlwgKCjIypdLlHBKa211huuKiorSMTExVscQbiIuLo7NmzcTGxvL8ePHmTFjBhMmTGDZsmWEhYUREhLCCy+8gIeHB9u2bbvyhlyuXDmUUqZm01pz/vx5/P392blzJ7GxsSQmJpKQkMDYsWPZtWsXd911F2XLliUyMpLnn3+eNm3aEB8fT2hoqOn5RPGhlNqstY4q1HOlIAh3o7Xm0KFDxMbGEhsbS9myZRkzZgyDBg0iOTmZ5s2b06JFC/r06eNWb6R2u53Dhw+zZcsWmjZtSmBgIPXr18dmsxEZGcngwYPp378/SUlJBAUFudW2CecpSkGQQ0bC5WVnZ7NhwwaWLVtGs2bN6NevH/369aN69epERkbSpk0bAObOnWtx0qKx2WzUrl2b2rVrX7kvKSmJuLg4tmzZQmBgIAA9evTg4MGDNG/enOjoaEaNGmVVZFHMSEEQLunMmTNkZ2cDUL9+fUJDQ4mOjqZ+/foopYiNjbU4oXMopahatSpVq1a9cl9MTAxJSUnExsaSmpoKwMMPP4xSiujoaDp37oy/v79VkYUbk0NGwmWcPn2aWbNmsWzZMjZv3syUKVMYOHAgx48fJywszOp4Lu3QoUMsX76cZcuWcfLkSWJjY9m0aROBgYFX7XGI4q8oh4xkHIKwTEZGBj/++CPPPvss69atIysri3379vHcc8+RmJjIwIEDAaQYFECtWrUYOXIkK1euZNOmTQD8/PPP3HbbbTRo0IDRo0djt9stTilcnRQE4VRaa7Kzs9m4cSOVKlVi3LhxVKhQgdDQUIKCgpg2bRrR0dGUKlXK6qhuy8PDA4Dnn3+e48eP8+mnn1KrVi1sNhuvvfYaDz/8MOvXr8eVjw4Ia8ghI+EUKSkpzJ07l+nTp/P6668THR3N6dOnqVSpktXRSpSEhATmzZvH9OnTKV++POvWrbtSQETxIN1OhUu6vDdw+PBhWrZsSc+ePRk+fDjt2rWTLpMWs9vt7Nq1i0aNGjFixAjS09MZPnw4UVGFeh8RLkTOIQiXkpqayuTJk2natCmzZs2iTp06HDhwgPnz59O+fXspBi7AZrPRqFEjAMaNG0fNmjXp168fXbp0sTiZsJLsIQhDaK3JzMzk7NmzhIeHc+edd/L444/TsWNHbDb53OEO7HY7+/fvp169evTv35/y5cvz+OOP07RpU6ujiZsgewjCMna7na+//pqoqCjeeecdKlasyMGDB1mwYAF33HGHFAM3YrPZqFevHgDvvfceISEhREdHc99991mcTDiL7CGIQsvOzqZdu3akp6fz8ssv06NHDykAxUxWVhZxcXFUq1aNRx99lIEDB3LHHXdYHUvcgOwhCKf65ZdfeOedd/Dw8GDKlCnExMQQHR0txaAY8vT0pEaNGgDceeedDB8+nE6dOrF582ZrgwlTyH+wKLBNmzbRpUsXHnnkkSuDxZo2bSoniUsAm83GgAED2LVrFwMGDCApKYnMzEy2b99udTRhICkIIl9Hjx4FYO3atfTt25c9e/ZcGUUsShYvLy+GDRtG9+7d2bdvH127duX+++9n7969VkcTBpCCIK5r//79PPDAA9xyyy0kJyfzzDPPMHz4cLy8vKyOJlxAREQEBw4coFmzZtx2223s2bPH6kiiiKQgiDytWrWKNm3a0LBhQ/bv30/FihWtjiRcUOnSpRkzZgyHDx+mfv36fPDBB4wYMYKUlBSro4lCkIIgrtBaM2/ePH7++Wfat2/P/v37GTdunEylLPJVpkwZAB544AEAGjVqxMKFC62MJApBCoIA4NixY/Ts2ZMJEyZQrlw5fH19KV++vNWxhJsJCgriww8/ZOHChWzZsgXImdVWuAcpCAKtNYMHD+bWW28lJiaG5s2bWx1JuLm2bdvy+uuvk5qaSr169fjkk09kdlU3IAWhBNu3bx+DBg0iLS2NH374gbFjx8oJY2GogIAAvvnmGyZNmkTXrl2Ji4uzOpK4ASkIJVBWVhbvvPMOt956Ky1btsTb2xtPT7maqjBHs2bN2LBhA127dsXT05PTp09fuTyqcC2GFASlVDel1F6l1AGl1Is3WK6lUipbKXWvEe2KwtmyZQurV69m06ZNjBw5UubDF6bz9PTkueeeIzg4mLfeeov27duze/duq2OJaxS5ICilPIDJQHegITBAKdXwOsu9Dawsapvi5qWnpzN+/Hhef/11WrZsyffff39lSgIhnOnNN99k4MCBtG/fnokTJ1odR+RixB5CK+CA1vqQ1joD+ALolcdyI4GvgSQD2hQ34ciRI7Rp04YdO3bwyCOPWB1HlHA2m40nn3ySzZs306BBA7TWnD592upYAmMKQhjwV67bcY77rlBKhQF9gKn5rUwp9ZhSKkYpFXPy5EkD4okpU6YwaNAgFi1aREhIiNVxhACgWrVqdO3aldjYWBo1asTatWutjlTiGVEQ8prZ7Nr+ZROBF7TW+Z5J0lpP11pHaa2jgoKCDIhXMmmt+e9//8vWrVt5++23+cc//iGT0AmX1KJFC6ZPn07v3r2ZNm2adE+1kBEFIQ6omut2FSD+mmWigC+UUkeAe4GPlFK9DWhb5OHixYsMHDiQefPmyeAy4RZ69uzJH3/8wY4dO6QgWMiIgrAJqKuUqqmU8gb6A0tzL6C1rqm1rqG1rgEsBJ7UWi82oG2Rh379+uHp6cnvv/9O9erVrY4jRIHUrVuXSZMmkZyczF133cXx48etjlTiFLkgaK2zgBHk9B7aDSzQWu9USj2ulHq8qOsXBbd+/XoyMzOZMWMGc+bMwc/Pz+pIQty0oKCgK2NkfvvtN6vjlChyCc1iQGvNu+++y/vvv8+aNWuuXBdXCHf2/fffM2LECLZs2SITLN4EuYRmCZaVlcWAAQNYsGABGzdulGIgio1u3bqxa9cuypQpw4wZM0hLS7M6UrEnBcGNZWVl4enpyV133cVvv/1G1apV83+SEG7E29ubzMxMVq1aRfv27Tlx4oTVkYo1KQhuav/+/TRu3Jjjx48zcOBAfH19rY4khCm8vb358ssv6dGjB+3atUPGJ5lHZjRzQ9u3b6d79+68/PLLVy52L0RxppTipZdeonXr1lSoUIHs7GyZg8sEsofgZrTWPProo7z//vsMGzbM6jhCOFW3bt24cOECzZo1Y9u2bVbHKXakILiRTZs2kZ6ezm+//cb9999vdRwhLOHv78+///1vunTpwrp166yOU6xIQXATK1eupGfPnuzZswdvb2+r4whhqX79+vHJJ5/w8MMPk56ebnWcYkMKghtYuXIlDz30EEuWLKFZs2ZWxxHCJXTv3p2tW7fi5eXFnj17rI5TLEhBcAPbt29n8eLFtGnTxuooQrgUX19fdu3aRfv27fn999+tjuP2pJeRC/vxxx+x2+2MHj3a6ihCuKxGjRoxf/587rnnHr755htuvfVWqyO5LdlDcFE//fQTAwYMkPmIhCiAzp07M2/ePJYsWWJ1FLcmcxm5oM2bN9OtWze+/vpr2rdvb3UcIdxKbGwsvr6+NGz4tyv5lggyl1ExU6FCBebNmyfFQIhC2LVrF9HR0TKiuRCkILiQCxcu8NxzzxEcHEzXrl2tjiOEW3rwwQcZMGAAffr0kS6pN0kKgouw2+089NBDnDp1Ch8fH6vjCOHWXn31VZo2bcqxY8esjuJWpCC4iLFjx5KcnMy0adPk2sdCFJHNZmPy5MnUqFGDlStXWh3HbUhBcBG1atVi0aJFsncghIFSUlIYNmwYixYtsjqKW5BxCBb7/fffSU5O5tFHH7U6ihDFTuXKlVmyZAndunWjRo0aREZGWh3JpckegoUOHTpEv379ZKyBECaKjIxk6tSpLF++3OooLk/GIVgkNTWVtm3b8sQTTzBixAir4whRImzevJkGDRpQqlQpq6OYRsYhuKGLFy/yyCOPSDEQwommTp3KoEGDsNvtVkdxSVIQLDBjxgx8fX157rnnrI4iRIny4YcfkpiYyCuvvGJ1FJckBcHJfvzxR15++WWrYwhRIvn4+PD111+zatUqLl26ZHUclyMFwYnOnDnDI488wsyZMylfvrzVcYQokSpXrswff/yBh4eHFIVrSEFwog0bNtCvXz+6dOlidRQhSjSlFK+++qpMLX8N6WXkJElJSVSqVMnqGEIIhzNnztC4cWNmz57NnXfeaXUcw0gvIxeXlJRE06ZN2b9/v9VRhBAO5cqVY+bMmQwdOpS0tDSr47gEKQgm01rz2GOPMWTIEOrWrWt1HCFELl26dGHFihX4+vpaHcUlyNQVJluyZAmHDx/myy+/tDqKECIPERERzJo1i8DAQHr37m11HEvJHoLJevbsyYoVK2TSOiFcWL169XjiiSdISkqyOoqlpCCYxG63M2TIEI4dO0ZYWJjVcYQQN3DrrbcyePBghg8fjit3tDGbFASTfPTRR+zdu5fq1atbHUUIUQAvv/wy4eHhJfoqa9Lt1AQHDhygdevWrF27lvDwcKvjCCFuwsGDB/Hz8yM0NNTqKIUi3U5dzKVLl/jggw+kGAjhhmbPns3YsWOtjmEJ6WVksMOHD1OzZk0aN25sdRQhRCGMHj2a8PBwtm/fTpMmTayO41Syh2AgrTUPPvggixcvtjqKEKKQAgICGDduXImcEVX2EAy0ePFizp8/z4ABA6yOIoQoguHDh5fI/2PZQzCI3W7nX//6FxMmTMDDw8PqOEKIIvD29sbT05Onn366RF1MRwqCQWw2G0uXLpWZTIUoJgICAti4cSOfffaZ1VGcRgqCAc6dO8e///1v6tSpY3UUIYRBlFJMmDCBsWPHlpjJ76QgGGDChAkcPXoUpZTVUYQQBmrXrh1dunRh9+7dVkdxCjmpXETx8fFMnjyZLVu2WB1FCGGCjz/+GK016enpxX5OMtlDKKJff/2VJ554gmrVqlkdRQhhknfffZd//etfVscwnUxdUQTZ2dnSo0iIEiAxMZGIiAhiYmKoWbOm1XFuSKausMiAAQP44YcfrI4hhDBZcHAwTz/9NC+99JLVUUxlSEFQSnVTSu1VSh1QSr2Yx+MDlVLbHV9rlVJNjWjXSocPH2b16tW0bdvW6ihCCCcYOXIk3bt3tzqGqYpcEJRSHsBkoDvQEBiglGp4zWKHgQ5a6ybAq8D0orZrtY8++oghQ4ZQunRpq6MIIZwgMDCQvn37smfPHqujmMaIPYRWwAGt9SGtdQbwBdAr9wJa67Va69OOm+uBKga0a6lKlSrx5JNPWh1DCOFEe/bs4c477yQzM9PqKKYwoiCEAX/luh3nuO96hgLfXe9BpdRjSqkYpVTMyZMnDYhnvPPnzzN69Ghq1apldRQhhBM1btyY2rVr880331gdxRRGFIS8RmPl2XVJKXU7OQXhheutTGs9XWsdpbWOCgoKMiCesbTW3HbbbWzYsMHqKEIIC4wcOZJJkyZZHcMURgxMiwOq5rpdBYi/diGlVBNgBtBda33KgHYt8dtvv5GWlkbLli2tjiKEsEDv3r0pX7681TFMYcQewiagrlKqplLKG+gPLM29gFKqGrAIeEhrvc+ANi0zadIkRowYgc0mPXaFKIk8PT1p3bo1P/30k9VRDFfkdzWtdRYwAlgJ7AYWaK13KqUeV0o97ljs30AF4COl1FallOuONstHp06dGDRokNUxhBAWunjxIvfeey+nTrntwY48yUjlm3D48GFq1Kghk9gJIRgyZAgNGjTghReue0rUEjJS2QnS0tJo3bo1Bw4csDqKEMIFjBgxgo8++qhYXUBHCkIBLViwgGbNmlG3bl2rowghXEBUVBQ//PBDsTqfWHy2xGQffvghI0eOtDqGEMKFhIaGMmvWLKtjGEYKQgForRkzZkyxn8dECHFzvLy8GDVqFElJSVZHMYQUhAI4fPgwXbt2lamuhRBX8fHxoUuXLnz77bdWRzGEFIQCeOSRR1i9erXVMYQQLqhPnz4sXrzY6hiGkG6n+UhOTqZ27dokJibi5+dnaRYhhOu5ePEiZ8+eJTg42OoogHQ7NdWyZcvo3LmzFAMhRJ5KlSrFiRMn2LZtm9VRisyIuYyKtebNmxMREWF1DCGEC/vjjz9Yv349c+fOtTpKkcgewg2kpaUREhIiE9kJIW7o7rvvZvny5W5/nQQpCDewcuVKHnjgAatjCCFcXJUqVahduza//fab1VGKRA4Z3cDixYvp1atX/gsKIUq8efPmERZ2o2uDuT7ZQ7iOrKwsli1bJgVBCFEgderU4eeff8aVe27mRwrCdaSnp/PKK69QrVo1q6MIIdyAUopRo0axdetWq6MUmhSE68jIyOCJJ56wOoYQwk0opejVq5dbD1KTgnAdrVq1YteuXVbHEEK4kV69ern1NBZSEPJw6tQpkpKSqF+/vtVRhBBupGXLlkycONHqGIUmBSEPMTExREZGFqt5zoUQ5vPx8aF27dokJydbHaVQ5B0vD2XLlmXYsGFWxxBCuKFXX32Vzz77zOoYhSLjEPLQpk0b2rRpY3UMIYQbatmyJWvWrLE6RqHIHkIeIiMji80FL4QQzhUVFYXVszQXluwhXCM+Pp5jx44RFBRk2Drtds2h5AvsSTxL6qVM0jPtZNs1Xh4KHy8PKvn7EBEaQHCAr2FtuquUCxnsOJ5KQuol0jLtZGTZsdkUPp42yvh4El7Zn/DKZfD0KNmfZTKy7Ow7cY79See4kJ5NWmY2WoOPlw1fTw9Cy/nRKKws5Up5Wx3VcsfPXGLn8VROXcggPTObjGw7njYbvl4elCvlRYOQstSoUAqllCHtRURE8Oijj6K1NmydziIF4RoxMTFERUUV6ReZeimTH3ed4M/jqew4nsruhLNcyMjO93kVy3jTMDSARqFlaVG9PB3Cg4r1G5/WmvWHUth4OIUd8ansPJ5KfGpavs/z9rRRP9ifiNAAGoWVpXODylQqW7yLafyZS/y0J4kdcansiE9l/4nzZGTb831emKMwNAoNoHXtCrSsEeiEtNbJyLLz894kthw7w874VHbGnyXlQka+z/P38aRBaFkiQsvStEo57mxYmTI+hXt79PT05NlnnyUrKwsvL69CrcMqcoGca3z55ZckJibyzDPP3PRzY4+d5rMNx1i2PZ60zPz/WfMT5O/DfVFV6N+yGlUDSxV5fa7i1Pl0vtocxxcbj3Hk1MUir8/TpujUoBIP3FKd9nUrut2nsuvJtmvW7E1i/oZjrNmbhN2Af9U6lcowoFU17o2sQkAp93qzupEjyRf4fOMxFm6O41QBCkB+Snt7cHezMAbeUo1GYQE3/fwxY8YQGBjI6NGji5zlZhXlAjlSEIooK9vOlzF/MW/9MXYnnDWlDZuC9uFBDGlbg471KpnShjNs/esMM38/zModiQX6dFsYVQP9GNCqGoPb1KB0IT/hWS31Uiaf/HGELzcdK9AeU2H4eNro2TiEoe1qEhF68294rmLVzkTmrDvC2oOnMOutrEmVAB68pTp9W1TBw1awDxvz5s3j22+/5csvvzQn1A1IQTCI1poBAwYwa9YsSpXK/xP5juOpjF643bRCkJcejYN5+e5GBPn7OK3NorqQnsU73+/h0/VHDfmUWxBh5fx4457GdAg37lyQM6zcmcj4xTtIOpfulPY8bIph7Wry7J3h+Hp5OKVNIySkXmLcNzv4aY/zOn80rVqOCfc2Ibyyf77L7tmzhx49enDo0CEnJLuaFASDHDlyhLZt2xIfH3/D5dIys5n4435m/HaILGe9w+VSrpQX43s2pG+LKk5v+2b9uu8kYxb9yfEzlyxp/57IMP4d3dDlT64mn0/nP0t2svzPBEvar1mxNG/d05hbalWwpP2C0lrz2cZjvLViD+fSs5zevreHjSc61mbEHXXwusH5PbvdzpAhQ5g5c6bTzyNIQTDIwoULmTt3LkuXLr3uMlv/OsOoL7dyKPmC03JdT4fwIN7t19Ql9xYuZmTx7yU7Wbg5zuooVCzjw1v3NObOhpWtjpKnZdvjGb94B6cvWnu1LaVg4C3VGB/dEB9P19tbiD9ziVELtrL+UIrVUahX2Z//3t+MhqFlrY7yN0UpCMW3C0sh7Nmz54aXy1y95wT9p69ziWIA8Mu+k/Sdspajp1wjz2WnL2TwwMcbXKIYQM6n78c+jeHTdUesjvI3U385yIjPtlheDAC0hnnrjzF41kbOpVmfJ7f9J87Rd8palygGAHtPnOO+aetYe+D6U1TMmDGDjz76yImpik4KQi7jxo1j3LhxeT62ZOtxHpu72ZDeQ0Y6lnKRvlPWsTfxnNVRAEg6m0a/aevY+tcZq6Ncxa5h/JKdTPppv9VRrnjzu9289d0eq2P8zfpDKfSfvp7TBvTWMcL2uDP0m7aOBJNOsBfW+fQshnyyiZU7E/N8PCMjgz///NPJqYpGCkIun376KYmJf//lrvgzgVELtllyvqAgks+nM3DGBg6dPG9pjpQLGQycsYEDSdbmuJH3ftjH1F8OWh2Dd1fuZdovzj/hWFA7488yaNZGzlq8p7Ar/iwPzdzIGRfYg8pLRpadkZ9t4ee9fz+5HRoamu/5SFcjBSGX1157jTNnzlx139qDyTzzxRayXbQYXHa5KJx0Uu+Ua6VlZjNo1gb2u3AxuOyt7/awYNNflrU/6/fDfPjzAcvaL6g/j6cy9JNNZJnURTg/cacv8tDMDaRecs1icFlGtp3HP91M7LHTV91ftWpVixIVnhSEXOLj4wkNDb1yO/VSJs8t2EZmtmsXg8sSUtMYs2i7JW2/u3IvO447r/ttUb307U6OWHAuaE/iWZc8THQ9m46ctqR4aa3551fbDBlk5gzpWXZGfbmVixn/6/nUokULlixZYmGqmycFweHcuXPY7XbKlv1fr4GXl+50ueOW+flxdxILYpz76XfDoVPM+uOwU9ssqosZ2Tz31TbsTtzzy8iy8+yX20wblGeWD1cf4M+4VKe2OfuPIy5zArmgjpy6yJsr/lfstdaMHTuW7Oz8p61xFVIQHHx8fPj555+vTHvw/Y5EFm05bnGqwnn1213EnS76lBAFcSE9i38u3Oa0AWdG2nz0NNN+dd5x/A9+2ufUQYxGybJrRi3YSlqmc97YDp48zzsr3WcvKrd5G47y2/6TQM41lmfOnMmJEycsTlVwUhAcLly4QOXKOf3UT51PZ+w37tU7ILdz6VmM/mo7zhhj8tryXfyVYs2gMyP894d97Ek0/016y7HTTHXhk8j52Z90nvdW7TW9nWy7ZtSCbS7Xm6+gtIbnF26/cjLe3U4sS0FwWL58OWPGjAHg498Ou82xy+tZd+gUv+439zJ+R5Iv8IWFJ2eNkJFtZ+IP5ndFfW/VPpfvmJCfT9YeIemsuYdQv9+RyDYX67J8sxJS0/h03VFACoLbunxCOS0zmy83HbM6jiHmrD1i6vrnrjtq2oRizvTD7hOmTq1xIOkcv99gAJO7yMzWzNtg7v+G2X+zzjJv/VGysu1MmTKFO+64w+o4BSYFweFyQVi6Nd4lRo0aYc3eJNNGMV/MyOKrze69d3BZtl0zb/1R09Y/Z61563a2zzYcIyPLnMM5u+LPsvGIe51Ivp6E1DRW7TqBzWbj+HH3ORcpBcHh9ttvp0OHDsxxwekNCsuucz7Fm+Hr2OOcS3P+5GJm+WLjMVNOmp5Ly2RRrGtM4WGE5PPprDBpAr7isndw2Sdrj7BixQomTJhgdZQCk4Lg0KtXL6hYk53x7tcL5Ea+ivmLSwW4WtvNcsV5gYri9MVMvt1m/LHehZvjCnS1PHcy14TfferFTJZsc59P0gWx8XAK9rIhcg7BHbVo0YJvN1k/pYHRzqZlGT6vUGJqGvtOuP6I5Jtlxkn4X/edNHydVtvy1xnDp7TYeCTFbXsW3cjx7LJudcjIPS8pZYKjR4+y75T50z4E+HnxTKe6dImoTJC/D2cuZvLLvpP894d9pg2C2x53hja1jZvnfnvcGcPWlZfA0t6MvKMOzauWo0Fo2StTMf97yQ7TDoEB/GnCdv153LwBXXfUr0TPxiE0rVqOIH8fPG2KuNOXWL3nBFN/OWTalA9aw464VNrWqWjYOs38m2pSJYCRd9ShfnBZAkt74+Vh4/TFDLb+dYZpvxwk9ph5bZ/I9OHZZ581bf1Gk4LgkJWVxe4T5g7m8vfx5Osn2lKnUpkr91Uu68F9UVXpEB7EPR+tNaW3y3aD35S2mzxqNbisLw/fWtPUNvJy5NRFUi9lEuBnzAVNjp+5RPJ587ovD25TnQ7XXFK1XrA/9YL96dk4lJ7/95tpF5HZftzogmDe31R4ZX86Nwy+6r7KZX3pGhHMHfUrce+UtWwzqf3dJy4y7fkhpqzbDHLIyKHRLR1M71309J11rxSDqWsO0vTlVfxn6U4g5w90bM8GprRr9KcvowvMtc6mZTLjt0OM+CzW1N4/eTFyigYz9jhyy8i2M3fdEXr+32/UG/cdvSf/QbzjA0W1CqW4r6V5k6sZ/Tdl5p7U4eQLjF64jdveXk342O/o9N4vV8Y6eHnY6NUszLS2j6VcpHpdc/6vzWBIQVBKdVNK7VVKHVBKvZjH40op9X+Ox7crpSKNaNdIY9+danobfSNzLnl5MSOL937YS+qlTOasPXKla2jnhpUp62f8TttfKZcMndve7De6uNOXeG35bpZtTyD5vHNnb91m4LaZ9anzsn98sZV/L9nJzvizpGfZ2frXGWb/ceTK47UqljatbSM/0f+VcpEUEweCbj56mq9i4og7fYmMbDsHT56/alqaTJPnlkov7ZpX6stLkQuCUsoDmAx0BxoCA5RSDa9ZrDtQ1/H1GDClqO0a7b1pn5i6/qqBfgSWzrmu79FTF6+aQfXyCVovDxsRoQGmtG/UoagL6VnFZpxGXow8ZHf8tLlTeuTVe8nH63//0vEmTsxo5OsUZ/LrlJuHTVGnUhnuaZ6zV3AuLZOvTL6yn92vvKnrN5IRH0dbAQe01ocAlFJfAL2AXbmW6QXM1TmT66xXSpVTSoVora25ovg17HY7MbHbKNe+iWltBJX533WPr+2hkftyhRVLm3MxeKP62DtrgjOrGLl9zn6tgvx9GNymOpCzF2rm+Aetc7bP16vo115Oy3LO6/T7C7dTpXypK7dPnE3jsbkxpl/QqU69knXIKAzIPWQ1znHfzS5jmezsbGxe5rwRF8TlGVbNlG7Q6FKj1uOq0g3s+ujM1yokwJfPH21NkL8v2XbNcwu2mT51u1HbZ9bI5/xULuvL7IdbUT/Y39R2Hhoy1NT1G8mIgpDXu9m1M9wUZJmcBZV6TCkVo5SKOXnSOX24bTYb2m7uH+XJXMfCy/pe3YuljM//dtSSTTqW6mEzpuh4GrQeV+XpYdz2Oeu1qh1UmoWO3muZ2Xae/XIr3+3I+zq/RjJq+zyc8IEI4La3f6bOv1Zwx7trWL495+BEYGlvnusSbmq7C774zNT1G8mIghAH5O7OUAW4dmheQZYBQGs9XWsdpbWOCgoKMiBe/mw2G/96YbSpbfyVcunKibMaFUrjleuNJ7xyTs+jzGw7O+PNORFpxK49gI9B63FVvp7GbZ9Rr/mNNA4LYMHwNoSV8+NiRhaPzo1hqQkjrvNi1PY543W6LMuuOZR8gcm5rgJXs2KZGzyj6HZs22Lq+o1kREHYBNRVStVUSnkD/YGl1yyzFBjk6G3UGkh1lfMHkHPI5mS8+d0bv3Yc0/Xz9uC5zvUo6+fJ4LY1qF4hpzfID7tOcPaSOf3Gy/oa03uptLeH6Z98lYLypbwoX8rrqjeLUt4eV+43i5G9vMzoMZZbm9oV+OzRW6hQxoeUCxkM/HgDa/Y6Z6+6tLeHYXudZr9O46Mb0LlhZUICfPHyUFQN9OPxDrWuPG7WBJCX2bKtuc55YRT5N6G1zlJKjQBWAh7ALK31TqXU447HpwIrgB7AAeAi8HBR2zXa1DfHEXj/m6a28X8/7uf2epWoU6kMj3eszeMda195LOlsGq8v321Ku6W9PahRwZguiJ4eNsIr+7PLxCt/hZXz4/cX/j5l8IvdG/Bi95wTdDVeXG5K20b28moYGsDVp86M9Uynuvg7Dj8Glvbmm6duverx9YdO0X/6elPaNvJ1Cq/sj5eHMu3a5V0jghl6W608H7uQnsUHP5l7PQx1xn0mNzSkNGutV5Dzpp/7vqm5ftbAU0a0ZZrTcXjaFFkmXsTkXHoWfaes5R931qVzw8pU8vflzMUMftl3kvdNnLqiUVgANgM/1TepEmBqQbBS4yrGvdE1CTOnC7EraGLg6+Tr5UF4ZX/TJpb8bMMxbq9fiRoVShPg50W2XRN/5hLrD53i498OceSUeTMUlPHxZPc2c4qyGWTqCgdPG9QJKsWeE+buPqZeyuTlb3fx8re78l/YIEb+8+asr5ypV0qLO33JtD2AG/H39TR0MFeDkLJ4e9jIMGngk1mf/gvCyMIJOX+jZhWEj9Yc5KM11kxc2SC4NCtWrKBv376WtH+zZOoKh0mTJtHU4D9yV9GkSjmD11c8X6dGoQGGdgH29rRRz+QujVZpavjflLHrcxXVynDl0rzuQAqCw/3330+L6oFWxzBF82rlDF1fvWB/SnkXv95GRr9OAJEmrNNqgaW9qWHwtBiR1dxnNO/NCPZKIzg4OP8FXYQUBIfevXtji9+Ov0/xOorWtnaFq0ZnGsHLw0bv5i4zrtAQSkG/KOMngzNjnVbr16KK4eusF+xf7PY8K5T2pnJWEmFh7vO/IgXBISQkhFMn4ulrwh+7lQa1qWHKegebtF6rtK8bRE0TJoNrFBZAi+rF59OvTcGDraubsm6z/latMqBVNTq0u5V//vOfVkcpMCkIDqGhocTHxzOoTXWcNHDSdGHl/Ojc0JyZFusF+9O6VvE5xDakbQ3T1j3YxHU72x31K1M10Ng9zsvuahpCBZPm8nI2T5viwdbV8ff3p1GjRlbHKTApCA633HILNWvWpFZQGdrVdc4IabM92Lq6YYOH8mLmm6gz1ahQio71zPudd28UTCV/n/wXdANm/s59PD3o36p4HGLrGhFMcIAvTz75JIsWLbI6ToFJQXDo0aMHQ4YMAeCxdnkPYnEn/r6e9DfxAikAnRsGm3KYxdmGtqtl6gSDXh42S64AZ7QGIWW5ra5xV0nLy0Ota+Dn5tOjKAVD2+X8vuPj4wkNDbU4UcFJQXA4fPjwlYJwW92K3Ovm5xJeuiuC8ibvfnvYFBPubYI7z3fXqmYgA1tVM72dYe1q0iisrOntmMXLQ/FuP/Omh78sOMCXF7vXN70dMw1uU+NKrykpCG7Kz8+P77///srt/9zVkLByfhYmKryuEZWddnI8qkYgj7Z3zz2q0t4evNevqaGjuK/Hy8PG+/c1w9vTPf/lnulU17SLN11rUJvq3Gbg9ZqdqVZQ6asK2t133y0FwR0FBQWRkpJCZmbOxWr8fb2Y0K+J251grljGmzf6NHZqm6M6h1OvsvsNwBrbs6FpJ0jzEl7Zn3+aPNWyGZpVLccTHes4rT2lFBP6NcHfoAkZncXDpnivX9OrJmR8//33KV3afQ6rSkFw8PDwoE6dOiQnJ1+5r23tim7XvfL1Po2pUMa5JzB9PD14//6mV03p7eo61gvigVvMP1R0rWG31aJVDffpneXrZeO9+5qa2jkhLyEBfrx0V4RT2yyqJzrUpnmuAXZ79+7lnnvusTDRzZOCkMuuXbsICQm56r5/9Whgag8UI/2zSzhdI6wZFRkRGsC7/Zq6xfmEBiFl+eD+5pa0bbMpJg+MpFaQ639q9PJQfDggktpB5l4v4Hr6tqjC8A7ucTiyR+Ngnu189d7fsWPHSE015/omZpGCkMvSpUv5888/r7rP29PG1Adb0LZ2BYtSFcxTt9dmxB11Lc3Qq1kYb/V17cNsdSqVYd7QVgSYeE2F/AT5+/DZsNZUc+LhqpvlYVNMvL85d5o0jqWgxnRv4PLdmzvVr8QH/Zv/bS/K3U4ogxSEq3z//fesWbPmb/f7enkwa0hLOtWv5PxQBfBc53BGd3WNnhn3RVXlv/c1c8nDR43DAvjysdZOP6SWl+AAXxYMb3PlanmuxMfTxkcDI+nZJCT/hZ3gpbsjXHZPIbpJCFMebIGXx9/fSu12O3XqOO/cixGkIOQSERHB9u3b83zM18uDaQ+14L4o1+mO6uWheKNPY0Z2snbP4Fq9m4fx8aCoq64VbbV2dSvyuYsUg8uCA3z5anhblzqnEODnxZxHWll26PF6xnRvwLieDZx+LuNGhrStwf/1b37dnmMPP/ww//nPf5ycqmikIOQSFRVFTEzMdR/39LDxzr1NmfpgJEEWjzxtUiWAb0feZsmJ0YLoWK8Sq55tb/n5lzI+nrzSK4K5j7RyqQJ1WUApLz5/rDVjezTA18vaf8duEcH8MKo9rWu55uHRYe1q8c2Tbalv8ZTioQG+zH64JS/dHXHDLssvvPACZ8+614WkVM7FzFxTVFSUvtEbtNEyMjI4cuQI4eH5dw1MvZjJq8t3sXCzcy+P5+tlY1TncIbeVsulPi3dyDdb4njl212cvpjp1HY7hAfxxj2N3WY8ydFTF3jh6+2sP5Ti1HaD/H145e4Iujd2jUNE+cnMtjNlzUE+XH3AtIsP5UUpGHhLNV7s3iDfDxcpKSnUqFGD06dP4+Hh3JHXSqnNWuuoQj1XCsLVtm/fTsWKFQt8Mui3/Sd5aelODp4090prAO3Dg3jl7gjD56J3huTz6by2bBffbk8g28TLlAJULuvD813ru+XMtVprPtt4jP/+sI/k8xmmtuXloejTPIyxPRpaepK9sPafOMf4JTucUkDrB/vz8t0R3FLAvadVq1bxxhtv5HlO0mxSEAz02GOP0aRJE0aMGHFTz1t7MJnPNhxj1c4Thn5qKVfKi76RVXjglmqWdf8zUmJqGl9u+osvNx0j3sBrSNsUtKubM7agU/1KeOZxks+dZGbbWbXzBPM3HGXdoVMY+W9aNdCP/i2rcV9UVcsPfRphb+I5PttwlEVbjnMuLcuw9fp42ujZOIQHbqlG1E2e53nzzTc5deoU7777rmF5CkoKgoE+/vhjfv/9d+bMmVOo5586n85Xm+P4dls8+06cIzP75l/f0t4eNK1ajn5RVejeKOSqkY/FRbZd8/OeJL7Y9BcxR1M4U4jDSTYFNSuWplujYPq3rObUUcfOdDj5Ap9vPMYPu05w5NSFQhWHwNLetKoRSP9WVWlfN8gp03U426WMbL7dFs/C2Dj+jEvlUmb2Ta/D28NG/RB/7m4ayr0tqlCuVOHmA8vOzubixYv4+zv/fIcUBANt3bqVgQMHsnPnziKvKz0rm72J59gZf5Ydx1PZnXCWs2lZpGdlk5Wt8fa04eNpI8jfh4jQACJCy9IoLICaFUoXy3/YG/kr5SI741OvvFYJqWmkZWaTnmXHw6bw8bRRxseTupX9aeR4nRqGlqWUt+udKDbT+fQsdh53vE7xqRxIOs+F9CzSMu1orfHx8sDH00ZYOT8iwgKuvFahbnIexSh2u+bgyfPsiE9l5/Gz7Iw/y6kL6aRn2cnIsuPpofD19CDAz4sGIWVpFFaWiNAA6gX759mF9GbNmTOHBx54AC8v5x+Kk4JgoMzMTBYuXMiAAQOc2q4QonhITEwkIiKC5ORkU6dVv56iFAT3PtBqAi8vL3r37s25c+esjiKEcEMxMTFERUVZUgyKSgpCHl544QWmT59udQwhhBvatGkTUVGF+oBuuZJ1ALaAoqKiWLFihdUxhBBuaPDgwU4fe2AU2UPIQ34jloUQIi9aa86cOUO1aq45g0B+pCDkoV69evTp0we73XmjIIUQ7u/YsWP06NHD6hiFJgUhDx4eHkyYMMHqGEIIN/Pdd9/RqVMntzyhDFIQruvbb7/lwQcftDqGEMKNLF68mN69e1sdo9CkIFxH69atWbFiBenp6VZHEUK4iYcffphu3bpZHaPQpCBcR1BQEI0bN2b16tVWRxFCuIETJ07Qp08fS6arMIoUhBt44oknyM6++flQhBAlzz/+8Y9Cz4HmKmTqigLQWrvtSSIhhPnS09MJDg5m9+7dBAdbe7U5mbrCRH379mX9+vVWxxBCuLA1a9bQoEEDy4tBUUlByEfDhg1ZvHix1TGEEC4sNDTU7a6fnBcpCPno3bu3FAQhxHXZ7XaqVq1K165drY5SZFIQ8hEZGUnz5s1l9lMhRJ42bdpE+/btrY5hCJncLh9KKb744gurYwghXNTixYuJjo62OoYhZA+hAPbu3UufPn2sjiGEcEHuPjo5NykIBVCrVi1+/fVXjhw5YnUUIYQLyczM5IEHHnDb6x9cSwpCAXh5eTF48GCmTJlidRQhhAtJTU1l/Pjx2GzF4620eGyFEzz11FOcP3/e6hhCCBdx5swZwsPDSUlJsTqKYaQgFFDt2rWZPHmyXCNBCAHA7Nmz6dq1K4GBgVZHMYwUhJuwY8cObr31Vlx5ug8hhPnsdjuTJ09m5MiRVkcxlBSEm9CwYUNOnTrF77//bnUUIYSFMjMzGTduHG3atLE6iqGKVBCUUoFKqR+UUvsd38vnsUxVpdTPSqndSqmdSqlnitKmlWw2GyNGjODDDz+0OooQwkL79u1j8ODBxW7Sy6LuIbwI/KS1rgv85Lh9rSzgOa11A6A18JRSqmER27XMkCFDqFmzphw2EqKE2rdvH506dSqWF88qakHoBVyeAHwO0PvaBbTWCVrrWMfP54DdQFgR27VM2bJleeutt6THkRAl1OTJkxk2bBi+vr5WRzFcUaeuqKy1ToCcN36lVKUbLayUqgE0BzbcYJnHgMcAqlWrVsR45oiLi6NVq1YcPnwYHx8fq+MIIZzk3LlzfPrpp2zbts3qKKbIdw9BKfWjUmpHHl+9bqYhpVQZ4GvgH1rrs9dbTms9XWsdpbWOCgoKupkmnKZKlSpERETw1VdfWR1FCOFEPj4+LFq0iKpVq1odxRT5FgSt9Z1a60Z5fC0BTiilQgAc35PyWodSyoucYjBfa73IyA2wysiRI5k0aZLVMYQQTqK1ZuXKlcVmZtO8FPUcwlJgsOPnwcCSaxdQOafhZwK7tdbvF7E9l9GzZ08GDhwo11wWooRYtWoV//rXv4pdz6LciloQ3gI6K6X2A50dt1FKhSqlVjiWuRV4CLhDKbXV8dWjiO1azsPDg6effpqtW7dKjyMhijm73c7YsWMZN25csS4IRTqprLU+BXTK4/54oIfj59+BYvkK2u12HnroId55551iMx+6EOLvvv/+e2w2G/fdd5/VUUwlI5WLwGaz8fbbb/PCCy+QlZVldRwhhEm6d+/OihUrivXeAUhBKLLo6GgqVqzI559/bnUUIYQJZs6cycqVK6lYsaLVUUwnl9AsIqUUn3/+eYn4YxGipDl9+jRjxoxhzZo1VkdxCtlDMEBoaCjbtm1jxowZVkcRQhjozTffpE+fPjRs6Laz7dwUKQgGqVChAi+88AInTpywOooQwgBaa44dO8ZLL71kdRSnkYJgkFq1ajFo0CBeeeUVq6MIIQxw6dIlvvjiC0JCQqyO4jRSEAw0btw4du/eTWZmptVRhBBFsGXLFqKiokrcFRKlIBioQoUKrF69uth3TROiONNaM3r0aEaOHInNVrLeIkvW1jrJnXfeyR9//GF1DCFEIaxcuZK4uDiGDRtmdRSnk4JggqFDh/LMM8/IoSMh3FDdunX55JNP8PLysjqK00lBMMGDDz5IxYoVefPNN62OIoS4CfPmzSMgIIDWrVtbHcUSUhBMoJRi5syZdO7c2eooQogC+vXXX3n++edL9GSVUhBMEhYWRqtWrXj77bdJS0uzOo4Q4gbOnTvHkCFDmDZtGq56YS5nkIJgIpvNRkxMDOPHj7c6ihDiBjZs2EC3bt246667rI5iKeXKu0dRUVE6JibG6hhFkpycTJMmTfjiiy+K9ZWWhHBXycnJxWouMqXUZq11VGGeK3sIJqtYsSLTpk3j999/tzqKEOIaKSkpNGvWjN27d1sdxSXIHoIT/fXXX8X24txCuKMBAwZQuXJlJk6caHUUw8geghvIysrijjvuYMWKFfkvLIQw3YoVK9i6dat0D89FCoKTeHp6Mn36dB599FFOnTpldRwhSrzOnTuzcuVK/Pz8rI7iMqQgONHtt99Ov379+M9//mN1FCFKLK01w4cP58iRI1SrVs3qOC5FrpjmZG+++Sbp6elorWUSPCEs8OGHH7Jp0yaqV69udRSXI3sITubn50dAQADR0dH89ttvVscRokRZuXIlb7zxBgsXLsTb29vqOC5HCoIFlFI8/fTT9OvXj4MHD1odR4gSIzg4mK+//ppatWpZHcUlSUGwSNeuXRk/fjx33323zIoqhMmSk5MZM2YMjRs3pm3btlbHcVlSECz01FNPMW/evBI5za4QzpKenk6fPn1QSpW4C97cLHl1LNa8eXMWL17Ms88+a3UUIYqdyz2KKlWqxGuvvWZ1HJcnBcEF3H777axcuZKPPvrI6ihCFCtaa5o3b87cuXNl76AApNupCwgICODbb7/l1ltvpWHDhnTs2NHqSEK4vaVLlxIYGMgzzzxjdRS3ISXTRdSuXZvvvvuO5s2bWx1FCLcXGxvL0KFDZRTyTZKC4EKaN2+Ol5cXffr0ITk52eo4Qril+Ph4evXqxZQpU2jRooXVcdyKFAQXU6pUKerVq8c999zDpUuXrI4jhNs5f/48Y8eO5d5777U6ituRguCC3njjDapWrcq9995boq/vKsTNOH/+PP/5z3+oUaMGjz/+uNVx3JKcVHZBNpuNOXPmsGnTJpRS2O126SEhxA2cP3+eHj16UK9ePTw95W2tsORdxkV5enrSpk0bli9fTnR0NGlpaVZHEsIlpaenEx0dTd26dZk2bZp8eCoCeeVcXNeuXSlTpgx9+/YlPT3d6jhCuBxvb29GjhzJxx9/LMWgiOTVc3Genp7Mnz8fPz8/GWkpRC7nz5/nrrvuYufOnfTt21eKgQHkYJsb8PLy4vPPPyc9PZ3ExERKly6Nv7+/1bGEsMzp06fp0aMHDRs2pEGDBlbHKTakpLoJLy8vypQpw8cff0znzp1JSUmxOpIQlhk0aBBt2rRhxowZeHh4WB2n2JA9BDczbtw4UlNT6dixI6tWrSI4ONjqSEI4TXx8PBUqVGD27NlUqFBBrjpoMNlDcDNKKSZMmMB9993Hli1brI4jhNPs27ePtm3bsmrVKipWrCjFwARSENyQUopx48bRvXt3pkyZwoYNG6yOJISpfvnlFzp06MD48eO56667rI5TbElBcHNhYWHcddddzJw50+ooQphm48aNzJ07l6FDh1odpViTcwhu7u677+bXX3+ld+/eAPIPI4qNS5cu8eSTTzJ06FBGjx5tdZwSQQpCMVC/fn02btyIzWbj4MGDlCpVipCQEKtjCVFoR48e5Z577qFevXpERkZaHafEkENGxUTZsmUpU6YMP/30Ey1btmTdunVWRxKi0B5++GEGDhzI/PnzKVWqlNVxSowiFQSlVKBS6gel1H7H9/I3WNZDKbVFKbWsKG2KG3vssceYOnUqvXr1Ys6cOVbHEaLAtNbMmTOHCxcu8P333zNq1CjpSeRkRd1DeBH4SWtdF/jJcft6ngF2F7E9UQDR0dH8/vvv1KhRA7vdTkZGhtWRhLihixcvMmjQICZOnMiZM2fw9va2OlKJVNSC0Au4/DF0DtA7r4WUUlWAnsCMIrYnCig8PJwOHTqwYMECOnbsSHx8vNWRhMhTRkYG7du3R2vNH3/8QVhYmNWRSqyiFoTKWusEAMf3StdZbiLwPGDPb4VKqceUUjFKqZiTJ08WMZ6477776NGjBy1btmTNmjVWxxHiKvv27cPb25vJkyfz6aefyvkCi+VbEJRSPyqlduTx1asgDSilooEkrfXmgiyvtZ6utY7SWkcFBQUV5CniBmw2G+PGjWPGjBksXboUALs937oshKlSUlIYPHgwPXr04OLFi9xyyy1yvsAF5FsQtNZ3aq0b5fG1BDihlAoBcHxPymMVtwJ3K6WOAF8Adyil5hm4DaIAunfvzvvvv09ycjINGzZk+fLlVkcSJdSWLVto1KgR5cqVY+vWrbJX4EKKeshoKTDY8fNgYMm1C2itx2itq2itawD9gdVa6weL2K4opIoVK/LRRx/x9NNPM3DgQE6dOmV1JFFCJCQksGfPHsLDw1m0aBEffPABZcqUsTqWyKWoBeEtoLNSaj/Q2XEbpVSoUmpFUcMJc9xxxx1s376datWqkZWVxYULF9BaWx1LFFNaa2bPnk3Tpk1Zs2YNpUuXpnXr1lbHEnlQrvxGEBUVpWNiYqyOUew9+eSTxMXFMWXKFOnhIQz3yCOPsHXrVmbNmkWzZs2sjlPsKaU2a62jCvNcGaksmDhxIpGRkTRr1oy5c+daHUcUA9nZ2Xz66adkZWUxZswYNm7cKMXADUhBEHh7e/PSSy+xevVqypUrh9aa48ePWx1LuKndu3fTrl07pk+fTkpKCnXr1sXTU6ZNcwdSEMQVjRs35u6772bfvn00adKEZ555hhMnTlgdS7iRvXv30q5dOwYOHMgvv/xCpUrXG5okXJEUBPE39erVY9euXSilaNiwIX/++afVkYQLS0pKYtSoUUyaNInw8HAOHDjAU089hc0mby/uRn5jIk+VK1dm4sSJbNu2jYiICL788kteffVVzp07Z3U04SK01owfP54GDRqQmZnJvffei1KKcuXKWR1NFJIUBHFDVapUwWaz0aJFC/bs2UOdOnV4//33yc7OtjqasMiFCxf45ZdfUEoRHBzM5s2bmTRpklyDoxiQgiAKpE6dOsyfP58ff/yRlJQUbDYbf/75p8ykWoKkpaXxwQcfUKdOHT799FMAnnrqKWrUqGFtMGEYOfUvbkrjxo1p3LgxAO+88w5//PEHL730EgMHDsTDw8PidMIMWmuUUowdO5b9+/fz/fff07RpU6tjCRPIwDRRJL/++itjx46lWbNmTJo0ifT0dHx8fKyOJQxw8uRJZs+ezaxZs1i9ejWVKlWS7qNuQAamCcu0b9+eX3/9lTfeeIP09HRq1arFoEGDWLt2rUyH4YYu/86+/vpr6taty+7du/nkk08IDQ2VYlACSEEQRaaUwt/fHx8fH7Zt20aTJk0YMmQIr7/+OgCXLl2yOKHIT0pKCv/973+vdDPu1KkThw8fZvbs2TLvUAkiBUEYqmLFivzzn/9k7969PPvss6SkpFC1alWGDh3Kxo0bZa/BhWit0VqzevVqatWqxebNm/n444+vTE1dvvx1L5EuiikpCMIUSilKly5NYGAgu3btIjw8nP79+zNt2jQgp+uisMaZM2eYNGkSjRs3ZvXq1bRu3ZqDBw8yb948brvtNrlQTQkmJ5WF09jtdjIyMoiPj6dFixb07NmT6OhounXrJoOZTHbq1Cn8/PzYv38/HTp0oFu3bgwfPpyOHTtKAShm5KSycAs2mw1fX19q1arFrl27aNeuHfPnz2f58uVkZmby/vvvs2fPHjmsZJD4+HjeeustbrvtNmrVqsWmTZto0KAB+/bt44svvuD222+XYiCuIgVBWCIkJIThw4fz7bffMnDgQM6fP8+BAwfo0qULderU4fPPPweQEdE3IS0tje+++46nnnqKo0ePkpSUxPHjxxk/fjwnTpygQ4cOeHt7y4Rz4rrkkJFwKVpr/vzzT7y9valRowbVq1fn1ltvJTo6mh49ehAcHGx1RJdy/PhxQkNDWbJkCYMHD6Zp06ZER0fz8MMPExQUZHU8YYGiHDKSjsXCpSilaNKkyZXbO3bs4LvvvmPZsmUcO3aMl156iRdffJFKlSpduahPSTn/cP78ecqUKcOqVav45JNPiI2NJSkpia1bt9KxY0cOHz5MYGCg1TGFG5M9BOF25syZw6ZNm4iNjSUgIIDvvvuOqVOncubMGSIjI2nevLnbfzq+dOkSfn5+zJo1i+XLlxMbG0taWhrx8fH88ccfHDx4kMjISBo0aCADxsRVirKHIAVBFAvfffcdP/74I7Gxsezfv59jx46xatUqNmzYQJMmTQgLC6Nx48b4+flZHfWK9PR0EhMT8fHxITg4mFdffZWNGzcSGxtL06ZNWbFiBQsWLCA7O5vIyEjq1q0r1xgQ+ZKCIEQulydji4mJYdGiRezYsYP4+Hjmz5/P6dOniY6OJiQkhODgYJ588kn69OnD9OnTCQgIIDg4mOrVqxd6Bk+tNefOnePkyZPUrl2b7du389NPP5GQkEBCQgJjx44FoG3btpw/f57KlSszevRonn76aT744AOqVq1KZGQk1atXlx5AolDkHIIQuVx+I42KiiIq6ur/C7vdzq5du0hISCAxMZGaNWsCOdcBjouLIyEhgbp16zJ79mzuu+8+lixZgpeXF56enpw4cYIlS5YwatQoALKyssjMzKRnz57MnTuXHj16sHLlSkqVKkW1atXYsWMHiYmJHDlyhJCQECIiIggMDKR8+fLs3buXChUqXPWJ/5lnnnHSKyRE3mQPQYjryM7OJjMzk6ysLLKysggICODSpUucOnUKpRSenp54enri5+dH6dKlycjIwGazyTF9YSnZQxDCBB4eHn+7xkOpUqUoVapUnst7e3s7I5YQppEzVEIIIQApCEIIIRykIAghhACkIAghhHCQgiCEEAKQgiCEEMJBCoIQQghACoIQQggHKQhCCCEAKQhCCCEcpCAIIYQApCAIIYRwkIIghBACkIIghBDCQQqCEEIIQAqCEEIIBykIQgghACkIQgghHKQgCCGEAKQgCCGEcJCCIIQQAihiQVBKBSqlflBK7Xd8L3+d5coppRYqpfYopXYrpdoUpV0hhBDGK+oewovAT1rrusBPjtt5+QD4XmtdH2gK7C5iu0IIIQxW1ILQC5jj+HkO0PvaBZRSZYH2wEwArXWG1vpMEdsVQghhsKIWhMpa6wQAx/dKeSxTCzgJzFZKbVFKzVBKlS5iu0IIIQzmmd8CSqkfgeA8Hhp7E21EAiO11huUUh+Qc2hp/HXaewx4zHEzXSm1o4DtuJuKQLLVIUwk2+feZPvcV73CPlFprQvdqlJqL9BRa52glAoB1mit612zTDCwXmtdw3G7HfCi1rpnAdYfo7WOKnRAF1actw1k+9ydbJ/7Ksq2FfWQ0VJgsOPnwcCSaxfQWicCfymlLheKTsCuIrYrhBDCYEUtCG8BnZVS+4HOjtsopUKVUityLTcSmK+U2g40A94oYrtCCCEMlu85hBvRWp8i5xP/tffHAz1y3d4KFGYXZnqhw7m+4rxtINvn7mT73Feht61I5xCEEEIUHzJ1hRBCCMCFCkJxnwajoNvnWNbDMWZjmTMzFkVBtk8pVVUp9bPj97ZTKfWMFVlvhlKqm1Jqr1LqgFLqbyPxVY7/czy+XSkVaUXOwijAtg10bNN2pdRapVRTK3IWVn7bl2u5lkqpbKXUvc7MV1QF2T6lVEel1FbH/9sv+a5Ua+0SX8A75HRHhZxxCm9fZ7k5wDDHz95AOauzG7l9jsdHAZ8By6zObeT2ASFApONnf2Af0NDq7DfYJg/gIDmDK72BbdfmJedc2XeAAloDG6zObeC2tQXKO37u7i7bVtDty7XcamAFcK/VuQ3+/ZUjp0dnNcftSvmt12X2ECj+02Dku30ASqkqQE9ghnNiGSbf7dNaJ2itYx0/nyNnTqswZwUshFbAAa31Ia11BvAFOduZWy9grs6xHijnGJPj6vLdNq31Wq31acfN9UAVJ2csioL87iCnB+TXQJIzwxmgINv3ALBIa30MQGud7za6UkEo7tNgFGT7ACYCzwN2J+UySkG3DwClVA2gObDB/GiFFgb8let2HH8vYAVZxhXdbO6h5OwJuYt8t08pFQb0AaY6MZdRCvL7CwfKK6XWKKU2K6UG5bfSInU7vVnOngbD2Yq6fUqpaCBJa71ZKdXRwGiGMOD3d3k9Zcj5VPYPrfVZI7KZROVx37Xd8gqyjCsqcG6l1O3kFITbTE1krIJs30TgBa11tlJ5Le7SCrJ9nkALcoYG+AHrlFLrtdb7rrdSpxYErfWd13tMKXVCKRWi/zcNRl67N3FAnNb68qfKhVx/ym2nM2D7bgXuVkr1AHyBskqpeVrrB02KfFMM2D6UUl7kFIP5WutFJkU1ShxQNdftKkB8IZZxRQXKrZRqQs7hy+46Z9yRuyjI9kUBXziKQUWgh1IqS2u92CkJi6agf5vJWusLwAWl1K/kXH7gugXBlQ4ZFfdpMAqyfWO01lV0zrxP/YHVrlIMCiDf7VM5/3kzgd1a6/edmK2wNgF1lVI1lVLe5PxOll6zzFJgkKO3UWsg9fKhMxeX77YppaoBi4CHbvSp0kXlu31a65pa6xqO/7eFwJNuUgygYH+bS4B2SilPpVQp4BbyuxaN1WfLc50Rr0DORXb2O74HOu4PBVbkWq4ZEANsBxbj6AXh6l8F3b5cy3fEvXoZ5bt95Bxy0I7f3VbHVw+rs+ezXT3I+UR1EBjruO9x4HHHzwqY7Hj8TyDK6swGbtsM4HSu31WM1ZmN3L5rlv0EN+plVNDtA0aT86F5BzmHaG+4ThmpLIQQAnCtQ0ZCCCEsJAVBCCEEIAVBCCGEgxQEIYQQgBQEIYQQDlIQhBBCAFIQhBBCOEhBEEIIAcD/A5xL4dw8J6AjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.4, -0.13333333333333336, 0.1333333333333333, 0.4]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.NonUniformImage at 0x7fb9c1816f90>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAD4CAYAAADhGCPfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQQElEQVR4nO3de6wc5X3G8eeJL5AAFnYoYGMDhpq2oGBDzSWYNnYIAR+qmjSVQkkBqUiIEhCJeqOiQlGpqqaoSZoogByTlqRUpAoutsBAwIEgoHa4GYNxYgwEMHYgAZeL0wRsfv1j5yTL8bvn7O68MzvY34+0Ors778z7O7Pjxzu7e/bniBAAjPS+QRcAoJkIBwBJhAOAJMIBQBLhACBp/KALGM1E7xF7aq9BlwHssn6hbXorfunUskaHw57aSyf4lEGXAeyyVsfKjss4rQCQRDgASCIcACQRDgCSCAcASYQDgCTCAUAS4QAgiXAAkEQ4AEgiHAAkEQ4AkggHAEmEA4AkwgFAEuEAIIlwAJCUJRxsn277R7Y32r6sw5j5ttfYXmf7+znmBVCd0l8TZ3ucpK9JOlXSJkkP2l4eEU+2jdlX0tWSTo+I523vX3ZeANXK8czheEkbI+KZiHhL0o2SFo0Yc7akpRHxvCRFxMsZ5gVQoRzhcJCkF9pubyrua3eEpMm277H9sO1zO23M9gW2H7L90Nv6ZYbyAPQjx7dPp77WemR33vGSflfSKZLeL+l/bK+KiA07rRixWNJiSZrkKXT5BQYkRzhskjSj7fZ0SZsTY34WEdskbbN9r6TZknYKBwDNkOO04kFJs2zPtD1R0lmSlo8Ys0zS79keb/sDkk6QtD7D3AAqUvqZQ0Rst32xpDskjZP0jYhYZ/vCYvm1EbHe9u2S1kp6R9KSiHii7NwAquOI5p7WT/KUoOMVUJ3VsVKvx6vJdnh8QhJAEuEAIIlwAJBEOABIIhwAJBEOAJIIBwBJhAOAJMIBQBLhACCJcACQRDgASCIcACQRDgCSCAcASYQDgCTCAUAS4QAgiXAAkEQ4AEiqrZFuMe442zts/3GOeQFUp3Q4tDXSXSjpSEl/YvvIDuO+oNZX2ANouLoa6UrSJZJukkQTXeA9oJZGurYPkvQJSdeOtTEa6QLNkCMcummk+2VJfxMRO8baWEQsjoi5ETF3gvbIUB6AftTVSHeupBttS9J+koZsb4+ImzPMD6ACOcLhV410Jb2oViPds9sHRMTM4eu2/13SLQQD0Gy1NNItOweA+tFIF9iN0UgXQM8IBwBJhAOAJMIBQBLhACCJcACQRDgASCIcACQRDgCSCAcASYQDgCTCAUAS4QAgiXAAkEQ4AEgiHAAkEQ4AkggHAEmEA4AkwgFAUi2NdG1/2vba4vKA7dk55gVQnboa6T4r6SMRcbSkKyUtLjsvgGrV0kg3Ih6IiK3FzVVqdcUC0GC1NNId4XxJt3VaSCNdoBlytMPrppFua6C9QK1wOLnTxiJisYrTjkme0tyOO8Aurq5GurJ9tKQlkhZGxCsZ5gVQoRynFb9qpGt7olqNdJe3D7B9sKSlks6JiA0Z5gRQsboa6V4h6YOSrrYtSdsjYm7ZuQFUh0a6wG6MRroAekY4AEgiHAAkEQ4AkggHAEmEA4AkwgFAEuEAICnH31ZghPfNOVK3rfjPjss3bX9T5x/c8W/PKvXC352kJy+6uuPyM586Tf/3kZdqrOjX3v/9A3TzrDs6Lj/y6os04x8eqLGiX7vu+fs0ffzeHZcvHDpb76x5ssaKqscnJDO6Y/OanteZuewCHfHnP8hfzAj91HbG8Wdo+6YX8xfTZvz0g3TrD27teb3Tps3JX8wIG645Xs8u6v17ieqoLZfRPiFJOGTSzz++dlUdUD/53El67K86P1MYy4a3t+mSQ+ZlrOjXvvrc/Tpiwl59rz/7qot04JeqeSbR1MczN8KhQmUPonZDRy3Qjq1bxx7YpZy15T7Ym1rbuMmTtWLd3dm21/SQ4G8rKpLzAJekFevu1viZh2TZVu7acm6vqbWNn3lI1mCQ8v+udSIc+rTxSydWst1b719WehtVHZA5ttvk2nLs+5SqjpWqEQ598HEf0tOfuray7Zc50K97/r58hSSUqa3q/0XL/O5V1vb0p66Vj/tQZduvCuHQh9uXfavyOa545pG+1hvt7bZcxh2wfy3r9Krf373ffd2LOo6Z3AiHHo074vBa5pm3Z+8PzY0v1PMZgBWPfreWdfrRzz7oZ1/3o65jJxfCoUcr7rmptrnGTZ7c0/jJ4z5QUSXvHb3ug173cRl1Hjs5EA4NlvuV85zeurP7d1V6GVu3Ju/jQSMcdhGT759S63x3H9X9K/u9jM2h7n2xq6qrka5tf6VYvtb2sTnmrdubtx826BI6unHm9wZdQmM0eV80+Rgaqa5GugslzSouF0i6puy8g3D/0Utrn7POc+LdzSD27SCOoX7V0ki3uP3NaFklaV/bUzPMvct7Y/4Rgy5hl8W+HV1djXS7brZLI9132+ceGoRVhX07uhzh0E0j3a6b7UbE4oiYGxFzJ2iP0sXlNG/tH9U+Z84/xMK7DWLfDuIY6leOcOimkW5XzXabbu/Tnxl0CR2d9exHB11CYzR5XzT5GBqplka6xe1zi3ctTpT0WkRsyTA3ClvnvVrrfAvWjXxZKc/YHOreF7uq0uEQEdslDTfSXS/pv4Yb6Q4305W0QtIzkjZK+rqki8rOuzsYOmrBoEvoaOKpz1Uytm5N3seDluU7JCNihVoB0H7ftW3XQ9Jncsw1aEPzP1nbx2B7PSfeuuPnu/1HqLfu+HlP4+t83WFo/iclPV3bfGXxCcke7dhQz4N7/y/e6Xmds2acVEElOxs65uO1rNOPfvZBP/u6H3UdO7kQDn04fdE5lc/x94f19yHSTdvfzFzJzna89HIt6/Sq39+9333dizqOmdwIhz7Eg4/r8G9fOPbAPpX53sGqv/K+TG1Vf59imd+9ytoO//aFigcfr2z7VSEc+vSbn1tVyXbPmFf+lf2qDvQc221ybTn2fUpVx0rVCIcSch/oQ0ct0PZn87yyn7u2nNtram3bn30u+7sXTf/26dEQDiWdNm1Otv9Rc79yftq0OZp9Vbl3jTe8va2SA/y0aXO04e1tpbYx+6qLste2Y+vWbI/nezkYJPpWZEXHq97R8WqwaGpTM3pl9odemfUjHAAk0fEKQM8IBwBJhAOAJMIBQBLhACCJcACQRDgASCIcACQRDgCSCAcASYQDgKRS4WB7iu07bT9V/Nyp+aDtGbbvtr3e9jrbl5aZE0A9yj5zuEzSyoiYJWllcXuk7ZL+IiJ+R9KJkj6TaLQLoGHKhsMiSdcX16+XdObIARGxJSIeKa6/oVZvi2SfTADNUTYcDhjuXFX83H+0wbYPlXSMpNWjjKGRLtAAYza1sX2XpAMTiy7vZSLbe0u6SdJnI+L1TuMiYrGkxVLr+xx6mQNAPmOGQ0R8rNMy2y/ZnhoRW2xPlZRsTmB7glrBcENELO27WgC1KXtasVzSecX18yQtGznAtiVdJ2l9RHyx5HwAalI2HP5J0qm2n5J0anFbtqfZHu6dOU/SOZI+antNcRkqOS+AipVqpBsRr0ja6UseI2KzpKHi+n2Skt9RB6C5+IQkgCTCAUAS4QAgiXAAkEQ4AEgiHAAkEQ4AkggHAEmEA4AkwgFAEuEAIIlwAJBEOABIIhwAJBEOAJIIBwBJhAOAJMIBQBLhACCJcACQVHkj3bax42w/avuWMnMCqEcdjXSHXapWn0wA7wGVN9KVJNvTJZ0haUnJ+QDUpK5Gul+W9NeS3hlrgzTSBZqh8ka6tv9A0ssR8bDt+WONp5Eu0Ax1NNKdJ+kPixZ4e0qaZPs/IuJP+64aQOUqb6QbEX8bEdMj4lBJZ0n6HsEANF8djXQBvAc5ormn9ZM8JU7wTn16AWSyOlbq9Xg12eiaT0gCSCIcACQRDgCSCAcASYQDgCTCAUAS4QAgiXAAkEQ4AEgiHAAkEQ4AkggHAEmEA4AkwgFAEuEAIIlwAJBEOABIIhwAJBEOAJIIBwBJtTTStb2v7e/Y/qHt9bY/XGZeANWrq5Huv0q6PSJ+W9Js0VAXaLzKG+naniTp9yVdJ0kR8VZE/G/JeQFUrI5GuodJ+qmkf7P9qO0ltvfqtEEa6QLNMGY42L7L9hOJy6Iu5xgv6VhJ10TEMZK2qfPphyJicUTMjYi5E7RHl1MAyK2ORrqbJG2KiNXF7e9olHAA0Ax1NNL9iaQXbP9Wcdcpkp4sOS+AitXVSPcSSTfYXitpjqR/LDkvgIqNeVoxmoh4Ra1nAiPv3yxpqO32Gklzy8wFoF58QhJAEuEAIIlwAJBEOABIIhwAJBEOAJIIBwBJhAOAJMIBQBLhACCJcACQRDgASCIcACQRDgCSCAcASYQDgCRHxKBr6Mj2TyU9N+g62uwn6WeDLqIN9YytaTU1rZ5DIuI3UgsaHQ5NY/uhiGjMN1pRz9iaVlPT6hkNpxUAkggHAEmEQ28WD7qAEahnbE2rqWn1dMRrDgCSeOYAIIlwAJBEOLSxPcX2nbafKn5O7jDux7Yft73G9kO9rp+7JtszbN9te73tdbYvbVv2edsvFrWusT00cv0u6zjd9o9sb7S9U69Tt3ylWL7W9rHdrltRPZ8u6lhr+wHbs9uWJR+/iuuZb/u1tsfhim7XHZiI4FJcJP2zpMuK65dJ+kKHcT+WtF+/6+euSdJUSccW1/eRtEHSkcXtz0v6y5I1jJP0tKTDJE2U9Njw9tvGDEm6TZIlnShpdbfrVlTPSZImF9cXDtcz2uNXcT3zJd3Sz7qDuvDM4d0WSbq+uH69pDNrXr+vbUbEloh4pLj+hqT1kg7KMPew4yVtjIhnIuItSTcWdY2s85vRskrSvkXn9W7WzV5PRDwQEVuLm6skTS85Z6l6Klq3UoTDux0QEVuk1j84Sft3GBeSvmv7YdsX9LF+FTVJkmwfKukYSavb7r64eHr9jT5PdQ6S9ELb7U3aOXw6jelm3SrqaXe+Ws9qhnV6/Kqu58O2H7N9m+2jely3dqUa6b4X2b5L0oGJRZf3sJl5EbHZ9v6S7rT9w4i4d8A1yfbekm6S9NmIeL24+xpJV6r1D+JKSf8i6c96LTFx38j3wDuN6WbdXnW9TdsL1AqHk9vuzvr4dVnPI2r9HcObxes+N0ua1eW6A7HbhUNEfKzTMtsv2Z4aEVuKp8Qvd9jG5uLny7b/W62nhvdK6mr9KmqyPUGtYLghIpa2bfultjFfl3RLNzWNsEnSjLbb0yVt7nLMxC7WraIe2T5a0hJJC6PVEV7SqI9fZfW0hbUiYoXtq23v1+3vMgicVrzbcknnFdfPk7Rs5ADbe9neZ/i6pI9LeqLb9SuqyZKuk7Q+Ir44YtnUtpufaKu1Fw9KmmV7pu2Jks4q6hpZ57nFuxYnSnqtOA3qZt3s9dg+WNJSSedExIa2+0d7/Kqs58DicZLt49X6t/dKN+sOzKBfEW3SRdIHJa2U9FTxc0px/zRJK4rrh6n1ivJjktZJunys9Wuo6WS1noqulbSmuAwVy74l6fFi2XJJU/usY0itd0GeHv6dJV0o6cLiuiV9rVj+uKS5o62bYb+MVc8SSVvb9sdDYz1+FddzcTHfY2q9QHpSlfsnx4WPTwNI4rQCQBLhACCJcACQRDgASCIcACQRDgCSCAcASf8PWc9e2JgkCJAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# each Dish now has a Stage with each variable!\n",
    "\n",
    "class Dish:\n",
    "    def __init__(self, x=0.0, y=0.0, d=0.1, ampl=1.0):\n",
    "        'Dish - a single element in a telescope array at position (x,y) with diameter d and electric field amplitud ampl'\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.d = d\n",
    "        self.ampl = ampl\n",
    "        self.xsc = Stage()\n",
    "        self.ysc = Stage()\n",
    "        self.dsc = Stage()\n",
    "        self.asc = Stage()\n",
    "        \n",
    "        def set_polar(self, r, theta):\n",
    "            'set position of Dish with polar coordinates'\n",
    "            self.x = r*np.cos(theta)\n",
    "            self.y = r*np.sin(theta)\n",
    "            \n",
    "        def get_polar(self):\n",
    "            r = np.sqrt(self.x*self.x + self.y*self.y)\n",
    "            theta = np.atan2(self.y, self.x)\n",
    "            return (r, theta)\n",
    "            \n",
    "\n",
    "class TelArray:\n",
    "    def __init__(self):\n",
    "        self.dishes = []\n",
    "        \n",
    "    def show_array(self):\n",
    "        fig, ax = plt.subplots(figsize=(6,6))\n",
    "        ax.set_xlim(-0.6,0.6)\n",
    "        ax.set_ylim(-0.6,0.6)\n",
    "        c = plt.Circle((0.0, 0.0), 0.5, fc=None, fill=False, linestyle='--', color='black', zorder=-5)\n",
    "        ax.add_artist(c)\n",
    "        for i, (a) in enumerate(self.dishes):\n",
    "            draw_circle = plt.Circle((a.x,a.y), a.d/2)\n",
    "            ax.add_artist(draw_circle)\n",
    "            ax.text(a.x, a.y, i, fontsize=16, color='white', weight='bold', horizontalalignment='center', verticalalignment='center')\n",
    "            \n",
    "            \n",
    "        plt.title('Telescope array')\n",
    "        plt.show()\n",
    "            \n",
    "    def add(self, Dish):\n",
    "        'add a Dish to the Array'\n",
    "        self.dishes.append(Dish)\n",
    "        return len(self.dishes)\n",
    "        \n",
    "    def add_linear(self, n=5, length=1., diam=0.1, ampl=1.0):\n",
    "        'a linear array of n Dishes with centre of first dish to last dish of length, diameter diam, transmission ampl'\n",
    "        xpos = np.linspace(-length/2., length/2., n)\n",
    "        for x in xpos:\n",
    "            dish = Dish(x, 0., diam, ampl)\n",
    "            self.add(dish)\n",
    "            \n",
    "        return len(self.dishes)\n",
    "    \n",
    "    def add_middle(self, n=5, diam=0.1, ampl=1.0):\n",
    "        'add n Dishes at the origin diameter diam, transmission ampl'\n",
    "        for x in xpos:\n",
    "            dish = Dish(0, 0., diam, ampl)\n",
    "            self.add(dish)\n",
    "            \n",
    "        return len(self.dishes)\n",
    "    \n",
    "    def make_grid(self):\n",
    "        'returns a HCIpy Grid containing all the dishes in TelArray'\n",
    "        # this seems very daft - surely there's a pythonic way of getting a quantity of a list of Objects?\n",
    "        xd = []\n",
    "        yd = []\n",
    "        diam = []\n",
    "        ampl = []\n",
    "        for dish in self.dishes:\n",
    "            xd.append(dish.x)\n",
    "            yd.append(dish.y)\n",
    "            diam.append(dish.d)\n",
    "            ampl.append(dish.ampl)\n",
    "            \n",
    "        print(xd)\n",
    "        interf_array = UnstructuredCoords((xd,yd))\n",
    "\n",
    "        grid = Grid(interf_array)\n",
    "\n",
    "        # currently we keep all dishes the same diameter - the first Diah\n",
    "        dish = circular_aperture(diam[0]) # this is a Field generator\n",
    "\n",
    "        interferom = make_segmented_aperture(dish, grid, np.array(ampl))\n",
    "        \n",
    "        return interferom\n",
    "    \n",
    "    \n",
    "    \n",
    "    def longest_time(self):\n",
    "        '''calculates the longest time of all the Stages listed in the TelArray'''\n",
    "        longest_time = 0.\n",
    "        \n",
    "        # loop over all the Dishes in the TelArray\n",
    "        for di in self.dishes:    \n",
    "            # loop through the Stages in the Dish object\n",
    "            for st in (di.xsc, di.ysc, di.dsc, di.asc):\n",
    "                # find longest time of all the Stages\n",
    "                if (st.total_time() > longest_time):\n",
    "                    longest_time = st.total_time()\n",
    "\n",
    "        return longest_time\n",
    "\n",
    "    \n",
    "    def t(self, t):\n",
    "        'given a TelArray, makes an HCIpy Grid at SINGLE time t'\n",
    "        print(0)\n",
    "        # WARNING: should barf if number of Dishes is not the same as 4* number of Scenes in the Play    \n",
    "    \n",
    "        xd = []\n",
    "        yd = []\n",
    "        diam = []\n",
    "        ampl = []\n",
    "        for dish in self.dishes:\n",
    "            xd.append(dish.xsc.t(t))\n",
    "            yd.append(dish.ysc.t(t))\n",
    "            diam.append(dish.dsc.t(t))\n",
    "            ampl.append(dish.asc.t(t))\n",
    "            \n",
    "        print(xd)\n",
    "        interf_array = UnstructuredCoords((xd,yd))\n",
    "\n",
    "        grid = Grid(interf_array)\n",
    "\n",
    "        # currently we keep all dishes the same diameter - the first Diah\n",
    "        dish = circular_aperture(diam[0]) # this is a Field generator\n",
    "\n",
    "        interferom = make_segmented_aperture(dish, grid, np.array(ampl))\n",
    "        \n",
    "        return interferom\n",
    "    \n",
    "    def tidyup(self):\n",
    "        # loop through all Scenes in your Dishes in your TelArray, bring them up to the latest time value\n",
    "        return 0\n",
    "        \n",
    "    def tidyup(self):\n",
    "        '''adds constant value Acts to brings up all the other Stages to the longest time'''\n",
    "\n",
    "#        print('longest time is {}'.format(self.longest_time()))\n",
    "        longest_time = self.longest_time()\n",
    "    \n",
    "        # loop through all the Dishes\n",
    "        for di in self.dishes:\n",
    "            # loop over all the Stages in each Dish object\n",
    "            for s in (di.xsc, di.ysc, di.dsc, di.asc):\n",
    "                stage_time = s.total_time()\n",
    "                if (stage_time < longest_time):\n",
    "#                print(\"adding a time delta to this one\")\n",
    "                    dt = longest_time - stage_time\n",
    "                    last_value = s.t(stage_time)\n",
    "                    s.add(Act(last_value,last_value,dt))\n",
    "\n",
    "    \n",
    "    # you can now build up your animation by adding TelArray[0].xsc.Add(asdf)\n",
    "    # and then run tidyup on your Telarray\n",
    "\n",
    "arr = TelArray()\n",
    "\n",
    "print(arr.add_linear(4,1-0.2,0.2))\n",
    "arr.show_array()\n",
    "        \n",
    "# start with two dishes\n",
    "\n",
    "interfer = arr.make_grid()\n",
    "\n",
    "pupil_grid = make_pupil_grid(1024,1.5)\n",
    "inter_pupil = evaluate_supersampled(interfer, pupil_grid, 8)\n",
    "\n",
    "imshow_field(inter_pupil)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.NonUniformImage at 0x7fb9bb142350>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAD4CAYAAADhGCPfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOs0lEQVR4nO3de4xc9XnG8efBF2hNXOxQgo25mMahAYlbHUNi0jghbmwXFaJUCiEF2qJatCUiUquWCimKhFQ1lRqlkQLRYpJAS0Kq4IJl2aTGudCE2sWAMZiNL1wCG7s4AWqCpcQxfvvHnKXj9Tu7cznnzNh8P9Jo58z5nfN7Z8/ss3N2duZ1RAgAxjqm3wUAGEyEA4AU4QAgRTgASBEOAFKT+13AeKb62DhO0/pdBnDU+oX2aX/80tm6gQ6H4zRNF/nSfpcBHLU2xvqW6zitAJAiHACkCAcAKcIBQIpwAJAiHACkCAcAKcIBQIpwAJAiHACkCAcAKcIBQIpwAJAiHACkCAcAKcIBQIpwAJAqJRxsL7G9zfZO2ze1GLPI9mbbW21/v4x5AVSn54+Jsz1J0pckLZY0IukR26si4ummMSdIulXSkoh4wfZJvc4LoFplPHNYIGlnRDwbEfsl3SPp8jFjrpK0MiJekKSI2FPCvAAqVEY4nCLpxablkeK2Zu+SNMP292w/avuaVjuzvdz2JtubfqVfllAegG6U8enT2cdaj+3OO1nS70i6VNKvSfov2xsiYvthG0YMSRqSpOmeSZdfoE/KCIcRSac2Lc+RtCsZ87OI2Cdpn+2HJJ0n6bBwADAYyjiteETSPNtzbU+VdKWkVWPG3C/p/bYn2/51SRdJGi5hbgAV6fmZQ0QcsH2DpG9LmiTpKxGx1fb1xfovR8Sw7QckbZF0UNKKiHiq17kBVMcRg3taP90zg45XQHU2xnq9Fq+k7fD4D0kAKcIBQIpwAJAiHACkCAcAKcIBQIpwAJAiHACkCAcAKcIBQIpwAJAiHACkCAcAKcIBQIpwAJAiHACkCAcAKcIBQIpwAJAiHACkamukW4x7j+03bP9hGfMCqE7P4dDUSHeppLMlfcL22S3GfU6Nj7AHMODqaqQrSZ+SdK8kmugCR4BaGunaPkXSRyV9eaKd0UgXGAxlhEM7jXS/IOlvI+KNiXYWEUMRMT8i5k/RsSWUB6AbdTXSnS/pHtuSdKKkZbYPRMR9JcwPoAJlhMObjXQl/USNRrpXNQ+IiLmj121/TdJqggEYbLU00u11DgD1o5Eu8BZGI10AHSMcAKQIBwApwgFAinAAkCIcAKQIBwApwgFAinAAkCIcAKQIBwApwgFAinAAkCIcAKQIBwApwgFAinAAkCIcAKQIBwApwgFAqpZGurY/aXtLcXnY9nllzAugOnU10n1O0gci4lxJt0ga6nVeANWqpZFuRDwcEa8WixvU6IoFYIDV0kh3jOskrW21kka6wGAoox1eO410GwPtD6oRDpe02llEDKk47ZjumYPbcQc4ytXVSFe2z5W0QtLSiHi5hHkBVKiM04o3G+nanqpGI91VzQNsnyZppaSrI2J7CXMCqFhdjXQ/I+ntkm61LUkHImJ+r3MDqA6NdIG3MBrpAugY4QAgRTgASBEOAFKEA4AU4QAgRTgASBEOAFKEA4AU4QAgRTgASBEOAFKEA4AU4QAgRTgASJXxMXFvaR8b3qPlv/H/n4p35XMf0qsLX6m1hm/v2nzI8ln/eY3O+PiW2uY/5vyztXbN1w+5bemyq3Rw89O11fD8N8/VtvffdchtH5l9fm3zS9KMH87UPXO/8+by0N7ZuvfdJ9VaQ5n4sJcuHfzABVr3ja+2XD939Z/pXcsfqbSGO174geZMPr7l+jp+OMYG06DVMHLgdV13WsvPMy7F9qH36LnLbm+5fvEn/kTHfP/xSmvo1ngf9kI4dGGiH4hmVf1wtFtDVT8cnjJVD/z4v9sau+T0BYpf7S+9honCsVm/j0OVNfSCT4Iq0ceG93Q0vpMHTxX7nDP5eE0656zSa2g3GDod265J55zVdjBI/T8OUuePnX4jHDrU/PeFI8Wadd8sdX9/vmNnLduMp+z7VIcj7bFTVyNd2/5isX6L7QvLmLdux0yb1tV2z3z9/NJqqOI3YKeumPZ6LduUrczvXbfHtNvHUD/U1Uh3qaR5xWW5pNt6nbcf1u74YVfb7Vz0tXIL6cIvLlvQ7xJKMwj3pdtj2u1jqB9qaaRbLN8VDRsknWB7Vglzo00jlx49Z5BH030ZZHU10m272S6NdKsxZ/3BfpdQmqPpvgyyMsKhnUa6bTfbjYihiJgfEfOn6NieiyvT0nkLu9rund/743IL6cJxq8t/xaBfBuG+dHtMu30M9UMZ4dBOI922mu0OuoP79nW13W9dtbm0GgbhtfL79rX/EmIv25StzO9dt8e028dQP9TSSLdYvqZ41eJiSXsjYncJc9duaO/sfpfQsWWLP17q/m6b985athlP2fepDkfaY6fncIiIA5JGG+kOS/q30Ua6o810Ja2R9KyknZJul/QXvc7bL53+r3wVv+k72efIgdf1xtZtpdew5PT2XzHoZGy73ti6TSMH2n95tN/HQer8sdNv/Pt0l3hvRQPvreC9FX0xyOEwindl8q7MUUfiuzIJBwAp3ngFoGOEA4AU4QAgRTgASBEOAFKEA4AU4QAgRTgASBEOAFKEA4AU4QAgRTgASBEOAFKEA4AU4QAgRTgASBEOAFKEA4BUT+Fge6btdbZ3FF9nJGNOtf1d28O2t9q+sZc5AdSj12cON0laHxHzJK0vlsc6IOmvIuLdki6W9JdJo10AA6bXcLhc0p3F9TslXTF2QETsjojHius/V6O3RdonE8Dg6DUc3jHauar4Ou7ncNs+Q9IFkjaOM4ZGusAAmDzRANsPSjo5WXVzJxPZPl7SvZI+HRGvtRoXEUOShqTGR9N3MgeA8kwYDhHx4VbrbL9ke1ZE7LY9S9KeFuOmqBEMd0fEyq6rBVCbXk8rVkm6trh+raT7xw6wbUl3SBqOiM/3OB+AmvQaDv8gabHtHZIWF8uyPdv2mmLMQklXS/qQ7c3FZVmP8wKo2ISnFeOJiJclHdavLiJ2SVpWXP+BpLTdFoDBxX9IAkgRDgBShAOAFOEAIEU4AEgRDgBShAOAFOEAIEU4AEgRDgBShAOAFOEAIEU4AEgRDgBShAOAFOEAIEU4AEgRDgBShAOAFOEAIFV5I92msZNsP257dS9zAqhHHY10R92oRp9MAEeAyhvpSpLtOZJ+X9KKHucDUJO6Gul+QdLfSDo40Q5ppAsMhsob6dq+TNKeiHjU9qKJxtNIFxgMdTTSXSjpD4oWeMdJmm77XyPij7quGkDlKm+kGxF/FxFzIuIMSVdK+g7BAAy+OhrpAjgCOWJwT+une2Zc5MP69AIoycZYr9filbTRNf8hCSBFOABIEQ4AUoQDgBThACBFOABIEQ4AUoQDgBThACBFOABIEQ4AUoQDgBThACBFOABIEQ4AUoQDgBThACBFOABIEQ4AUoQDgFQtjXRtn2D7W7Z/ZHvY9nt7mRdA9epqpPvPkh6IiN+WdJ5oqAsMvMob6dqeLul3Jd0hSRGxPyL+t8d5AVSsjka6Z0r6qaSv2n7c9grb01rtkEa6wGCYMBxsP2j7qeRyeZtzTJZ0oaTbIuICSfvU+vRDETEUEfMjYv4UHdvmFADKVkcj3RFJIxGxsVj+lsYJBwCDoY5Guv8j6UXbZxU3XSrp6R7nBVCxuhrpfkrS3ba3SDpf0t/3OC+Aik14WjGeiHhZjWcCY2/fJWlZ0/JmSfN7mQtAvfgPSQApwgFAinAAkCIcAKQIBwApwgFAinAAkCIcAKQIBwApwgFAinAAkCIcAKQIBwApwgFAinAAkCIcAKQcEf2uoSXbP5X0437X0eREST/rdxFNqGdig1bToNVzekT8ZrZioMNh0NjeFBED84lW1DOxQatp0OoZD6cVAFKEA4AU4dCZoX4XMAb1TGzQahq0elribw4AUjxzAJAiHACkCIcmtmfaXmd7R/F1Rotxz9t+0vZm25s63b7smmyfavu7todtb7V9Y9O6z9r+SVHrZtvLxm7fZh1LbG+zvdP2Yb1O3fDFYv0W2xe2u21F9XyyqGOL7Ydtn9e0Lj1+FdezyPbepuPwmXa37ZuI4FJcJP2jpJuK6zdJ+lyLcc9LOrHb7cuuSdIsSRcW198mabuks4vlz0r66x5rmCTpGUlnSpoq6YnR/TeNWSZprSRLuljSxna3raie90maUVxfOlrPeMev4noWSVrdzbb9uvDM4VCXS7qzuH6npCtq3r6rfUbE7oh4rLj+c0nDkk4pYe5RCyTtjIhnI2K/pHuKusbWeVc0bJB0QtF5vZ1tS68nIh6OiFeLxQ2S5vQ4Z0/1VLRtpQiHQ70jInZLjR84SSe1GBeS/sP2o7aXd7F9FTVJkmyfIekCSRubbr6heHr9lS5PdU6R9GLT8ogOD59WY9rZtop6ml2nxrOaUa2OX9X1vNf2E7bX2j6nw21r11Mj3SOR7QclnZysurmD3SyMiF22T5K0zvaPIuKhPtck28dLulfSpyPiteLm2yTdosYPxC2S/knSn3ZaYnLb2NfAW41pZ9tOtb1P2x9UIxwuabq51OPXZj2PqfE+hteLv/vcJ2lem9v2xVsuHCLiw63W2X7J9qyI2F08Jd7TYh+7iq97bP+7Gk8NH5LU1vZV1GR7ihrBcHdErGza90tNY26XtLqdmsYYkXRq0/IcSbvaHDO1jW2rqEe2z5W0QtLSaHSElzTu8ausnqawVkSssX2r7RPbvS/9wGnFoVZJura4fq2k+8cOsD3N9ttGr0v6PUlPtbt9RTVZ0h2ShiPi82PWzWpa/GhTrZ14RNI823NtT5V0ZVHX2DqvKV61uFjS3uI0qJ1tS6/H9mmSVkq6OiK2N90+3vGrsp6Ti+Mk2wvU+Nl7uZ1t+6bffxEdpIukt0taL2lH8XVmcftsSWuK62eq8RflJyRtlXTzRNvXUNMlajwV3SJpc3FZVqz7F0lPFutWSZrVZR3L1HgV5JnR+yzpeknXF9ct6UvF+iclzR9v2xK+LxPVs0LSq03fj00THb+K67mhmO8JNf5A+r4qvz9lXPj3aQApTisApAgHACnCAUCKcACQIhwApAgHACnCAUDq/wAvvtIM/E97tgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ndish = 5\n",
    "\n",
    "ddish = .1\n",
    "\n",
    "dish_x = np.linspace(-0.5+(ddish/2.), 0.5-(ddish/2.), ndish)\n",
    "\n",
    "dish_y = np.zeros_like(dish_x)\n",
    "\n",
    "interf_array = UnstructuredCoords((dish_x,dish_y))\n",
    "\n",
    "grid = Grid(interf_array)\n",
    "\n",
    "dish = circular_aperture(ddish) # this is a Field generator\n",
    "\n",
    "interferom = make_segmented_aperture(dish, grid, np.array([1.0]))\n",
    "\n",
    "pupil_grid = make_pupil_grid(1024,1.5)\n",
    "inter_pupil = evaluate_supersampled(interferom, pupil_grid, 8)\n",
    "\n",
    "imshow_field(inter_pupil)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "focal_grid = make_focal_grid(8, 14)\n",
    "\n",
    "prop = FraunhoferPropagator(inter_pupil, focal_grid)\n",
    "\n",
    "wf = Wavefront(inter_pupil)\n",
    "focal_image = prop(wf).intensity\n",
    "imshow_field(np.log10(focal_image / focal_image.max()), vmin=-5, cmap='inferno')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HCIpy setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddish = Stage()\n",
    "ddish.add(Act(1.0, 0.05, 5,  'sig',10))\n",
    "ddish.add(Act(0.05, 0.1, 5,  'sig',10))\n",
    "ddish.add(Act(0.1, 0.2, 5,  'sig',10))\n",
    "\n",
    "print(ddish.total_time())\n",
    "\n",
    "t = np.arange(0,ddish.total_time(),0.05)\n",
    "\n",
    "plt.plot(t, ddish.t(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_time_animation = ddish.total_time() # seconds\n",
    "frame_rate = 10          # frames per second\n",
    "total_frames = total_time_animation * frame_rate # seconds\n",
    "supersampling = 4\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "\n",
    "anim = FFMpegWriter('interferometers_{}.mp4'.format(anim_version), framerate=frame_rate)\n",
    "\n",
    "times = np.arange(0,total_time_animation, 1./frame_rate)\n",
    "\n",
    "dd = ddish.t(times)\n",
    "\n",
    "ndish = 5\n",
    "pupil_grid = make_pupil_grid(1024,1.5)\n",
    "focal_grid = make_focal_grid(8, 14)\n",
    "\n",
    "for t in progressbar(range(np.int(total_frames))):\n",
    "\n",
    "    dish = circular_aperture(dd[t]) # this is a Field generator\n",
    "\n",
    "    dish_x = np.linspace(-0.5+(dd[t]/2.), 0.5-(dd[t]/2.), ndish)\n",
    "\n",
    "    dish_y = np.zeros_like(dish_x)\n",
    "\n",
    "    interf_array = UnstructuredCoords((dish_x,dish_y))\n",
    "\n",
    "    grid = Grid(interf_array)\n",
    "    \n",
    "    interferom = make_segmented_aperture(dish, grid, np.array([1.0]))\n",
    "\n",
    "    inter_pupil = evaluate_supersampled(interferom, pupil_grid, supersampling)\n",
    "\n",
    "    plt.clf()\n",
    "#    plt.suptitle('Timestep %d / %d' % (t, total_frames))\n",
    "    \n",
    "    plt.subplot(1,2,1)\n",
    "    imshow_field(inter_pupil, cmap='gray')\n",
    "    plt.title('Aperture')\n",
    "    \n",
    "    plt.text(0.05,0.05,anim_version, fontsize=12, color='white',\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    "            transform=plt.gca().transAxes)\n",
    "    \n",
    "    prop = FraunhoferPropagator(inter_pupil, focal_grid)\n",
    "\n",
    "    wf = Wavefront(inter_pupil)\n",
    "    focal_image = prop(wf).intensity\n",
    "    \n",
    "    plt.subplot(1,2,2)\n",
    "    imshow_field(np.log10(focal_image / focal_image.max()), vmin=-5, cmap='inferno')\n",
    "    #plt.title('frame {} at time {:.2f}'.format(timestep, timestep/frame_rate))\n",
    "    plt.title('PSF')\n",
    "\n",
    "    plt.text(0.95,0.05,\"@mattkenworthy\", fontsize=12, color='white',\n",
    "        horizontalalignment='right', verticalalignment='bottom',\n",
    "        transform=plt.gca().transAxes)\n",
    "    \n",
    "    anim.add_frame()\n",
    "\n",
    "plt.close()\n",
    "anim.close()\n",
    "\n",
    "# Show created animation\n",
    "anim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ffmpeg -i interferometers_v1.mp4 -vf \"fps=10,scale=400:-1:flags=lanczos\" -c:v pam -f image2pipe - | convert -delay 10 - -loop 0 -layers optimize interferometers_v1.gif"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}