LSST simulation generator-variable PSF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LSST Simulation generation\n",
    "This notebook is a wrapper for StarFast, which generates simulated images with realistic stellar distributions and spectra, psf and filters, and detailed atmospheric differential chromatic refraction (DCR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import imp\n",
    "import numpy as np\n",
    "import pickle\n",
    "import galsim\n",
    "from lsst.afw.geom import Angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "# from matplotlib import pylab as plt\n",
    "import lsst.afw.display.rgb as rgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ramp_distribution(rng, x0, x1, x2, invert=True):\n",
    "    # A distribution with a linear increase in probability up to a point, and flat after that\n",
    "    # The probability of measuring x' < x is P(x'<x) = Integral{x0,x: c1*x dx}\n",
    "    # if x is greater than x1, the probability is P(x' < x) = P(x' < x1) + Integral(x1,x: c2 dx)\n",
    "    # The probability of measuring less than x2 is 1.0, and c1*x1 = c2, which allows us to solve for c1:\n",
    "    c1 = 1./(x1*x2 - (x0**2. + x1**2.)/2.)\n",
    "    p1 = c1*(x1**2. - x0**2.)/2.\n",
    "    if invert:\n",
    "        # Calculate x from P(x)\n",
    "        p = rng.rand()\n",
    "        if p <= p1:\n",
    "            x = np.sqrt(2.*p/c1 + x0**2.)\n",
    "        else:\n",
    "            x = (p - p1)/(c1*x1) + x1\n",
    "        return((p, x))\n",
    "    else:\n",
    "        # Calculate P(x) from x\n",
    "        x = rng.rand()*(x2 - x0) + x0\n",
    "        if x <= x1:\n",
    "            p = c1*(x**2. - x0**2.)/2.\n",
    "        else:\n",
    "            p = p1 + c1*x1*(x - x1)\n",
    "        return((x, p)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'elevation')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwAAAAIPCAYAAADEoLnvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYLVdZL+DfB5FMkIQEEAQZJSQXcCCIECSEQQWUQQgaGQS5oCg4Ac4gqOi9KkocQdGAwNUgIFExDEqAIIhKEOMQAgJBRiOEJEASIOG7f1S1NM3uc/qcrh72qfd9nv3U2bXWrlq7T50+9atVa1V1dwAAgHm4xk43AAAA2D4CAAAAzIgAAAAAMyIAAADAjAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMCAAAAzIgAAAAAMyIAAADAjBy00w1YdlX1viRHJLlwh5sCAMCB7eZJLuvuW2xmIwLA5h1x6KGHHn388ccfvdMNAQDgwHX++efniiuu2PR2BIDNu/D4448/+txzz93pdgAAcAA74YQT8va3v/3CzW7HGAAAAJgRAQAAAGZEAAAAgBkRAAAAYEYEAAAAmBEBAAAAZkQAAACAGREAAABgRgQAAACYEQEAAABmRAAAAIAZEQAAAGCjqobXEhMAAABgRgQAAACYEQEAAABmRAAAAIAZOWinGwAAALvSngb7Lirr3rq2TEgPAAAAzIgeAAAAWGTRFf2VK/9LcrV/ET0AAAAwIwIAAADMiAAAAAAzIgAAAMCMLGUAqKpvrarXVtUHq+qKqnpvVb20qu6yTv0Tq+qsqrq4qi6vqvOq6keq6prb3XYAAJZY91IPAE6WMABU1S8neWWSOyR5dZLfSPL2JA9M8uaqesSa+g9Mck6Sk5K8IsnvJLlWkmcnOWP7Wg4AADtvqaYBraobJnlKkv9K8tXdfdGqsnskOTvJzyd58bjuiCTPS3J1kpO7+23j+qeNdU+pqlO7WxAAAGAWlq0H4GYZ2vz3q0/+k6S7X5/kk0muv2r1KeP7M1ZO/se6VyZ56vj2+7e0xQAAsIssWwB4d5LPJrlTVV1vdUFVnZTkOkn+ZtXqe47LVy/Y1jlJLk9yYlUdvAVtBQCAXWepbgHq7our6ieS/HqSf6+qM5N8PMmtkjwgyV8n+b5VH7nNuHzXgm1dVVXvS3LbJLdMcv6e9l1V565TdNw+fQkAANhBSxUAkqS7T6uqC5OcnuRxq4r+I8kL1twadOS4vHSdza2sP2rSRgIAwC61bLcApap+PMnLkrwgw5X/w5OckOS9Sf5fVf3KvmxuXO51LqfuPmHRK8k79+kLAACwcVXDi8ksVQCoqpOT/HKSv+juJ3X3e7v78u5+e5JvT/KhJE+uqluOH1m5wn/kl24tSXLEmnoAAHBAW6oAkOTbxuXr1xZ09+VJ/iHDd/q6cfUF4/LYtfWr6qAkt0hyVYbeAwAAOOAtWwBYma3n+uuUr6z/7Lg8e1zeZ0Hdk5IcluQt3f2ZaZoHAAC727IFgDeNy++tqhuvLqiq+ya5a5Irk7xlXP2yJB9LcmpV3XFV3UOSPHN8+5wtbTEAAOwiyzYL0MsyzPN/7yTnV9Urknw0yfEZbg+qJD/Z3R9Pku6+rKoeN37uDVV1RpKLM0wZeptx/Uu2/VsAAPCl9jTYd1FZ73UeFxZYqgDQ3Z+vqvsleUKSUzMM/D0sw0n9WUl+s7tfu+YzZ1bV3ZP8TJKHJDkkw5ShTxrrO3IAALbLyom8U7Ads1QBIEm6+3NJThtfG/3Mm5Pcb8saBQDA5i0KBQLD5JZtDAAAALAJAgAAAMyIAAAAADOydGMAAABYEmb12ZUEAAAAdi+hYHICAAAAW8OsPruSMQAAADAjAgAAAMyIAAAAADMiAAAAwIwYBAwAwPYx+HfH6QEAAIAZEQAAAGBGBAAAgDmr2vMTezngCAAAADAjAgAAAMyIAAAAADMiAAAAwIx4DgAAwFzsabDvojJz9h+Q9AAAAMCM6AEAAJiLRVf0V678u9o/G3oAAABgRgQAAACYEQEAAABmRAAAAIAZMQgYAGDODP6dHT0AAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMCAAAAzIgAAAAAMyIAAADAjAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAASVI1vOAAJwAAAMCMCAAAADAjAgAAAMyIAAAAADNy0E43AABg2+1psO+isu6tawtsMz0AAAAwI3oAAID5WXRFf+XKv6v9HOD0AAAAwIwIAAAAMCMCAAAAzIgAAAAAM2IQMABAYvAvs6EHAAAAZkQAAACAGREAAABgRgQAAACYEQEAAABmRAAAAIAZEQAAAGBGBAAAAJgRAQAAAGZkqQJAVT26qnovr6sXfO7Eqjqrqi6uqsur6ryq+pGquuZOfA8AANgpB+10A/bRO5L83Dpld0tyzySvWr2yqh6Y5OVJrkzykiQXJ7l/kmcnuWuSh25VYwEAYLdZqgDQ3e/IEAK+RFX93fjH31+17ogkz0tydZKTu/tt4/qnJTk7ySlVdWp3n7GlDQcAgF1iqW4BWk9V3S7JnZN8KMlfrSo6Jcn1k5yxcvKfJN19ZZKnjm+/f7vaCQAAO+2ACABJvm9c/mF3rx4DcM9x+eoFnzknyeVJTqyqg7eycQAAsFss1S1Ai1TVoUkekeTzSf5gTfFtxuW71n6uu6+qqvcluW2SWyY5fy/7OXedouP2qcEAALCDDoQegO9IclSSV3X3B9aUHTkuL13nsyvrj9qKhgEAwG6z9D0ASb53XP7efny2xmXvrWJ3n7BwA0PPwB32Y98AALDtlroHoKr+V5ITk3wwyVkLqqxc4T9yQVmSHLGmHgAAHNCWOgBk/cG/Ky4Yl8euLaiqg5LcIslVSd67Nc0DAIDdZWkDQFUdkuSRGQb//uE61c4el/dZUHZSksOSvKW7PzN9CwEAYPdZ2gCQ4Qm+101y1oLBvyteluRjSU6tqjuurBzDwzPHt8/Z0lYCAMAussyDgFcG//7+ehW6+7KqelyGIPCGqjojycVJHpBhitCXJXnJVjcUAAB2i6XsAaiq45N8Y9Yf/Ps/uvvMJHfP8OCvhyT5wSSfS/KkJKd2915nAAIAgAPFUvYAdPf5+cIUnhup/+Yk99u6FgEAwHJYyh4AAABg/wgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMCAAAAzIgAAAAAMyIAAADAjAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMCAAAAzIgAAAAAMyIAAADAjAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAABQNbwAZkAAAACAGREAAABgRgQAAACYEQEAAABm5KCdbgAAbKs9DfZdVNa9dW0B2AF6AAAAYEb0AAAwL4uu6K9c+Xe1H5gBPQAAADAjAgAAAMyIAAAAADMiAAAAwIwYBAwABv8CM6IHAAAAZkQAAACAGREAAABgRgQAAACYEQEAAABmRAAAAIAZ2bJpQKvq8CRHJbnmovLu/s+t2jcAALDY5AGgqh6Z5CeSHL+Har0V+wYAAPZs0pPwqnp0ktOTXJ3kTUk+kOSqKfcBAADsv6mvwj8lySeSfGN3nz/xtgEAgE2aehDwVyV5mZN/AADYnaYOABcnuXLibQIAABOZOgC8MsnJVVUTbxcAAJjA1AHgp5IcnOS5VXXtibcNAABs0tSDgF+a5PIkj03ysKp6d5JLFtTr7r7XxPsGAAD2YuoAcPKqPx+e5GvXqdcT7xcAANiASQNAd099SxEAADAhJ+wAADAjAgAAAMzIlgSAqjq1qv6mqj5eVVdV1cVV9ddVdepW7A8AANiYSccAjPP/vzDJw5JUkquT/HeS6yW5V5J7VtUDuvthU+4XAADYmKl7AL4vycOTvD3JvZMc0t03SnLI+P7cJN9ZVY+feL8AAMAGTB0AHpPkwiQndffZ3X11knT31d19dpK7j+X/e7M7qqq7VdXLq+ojVfWZcfnaqrrfgronVtVZ461Il1fVeVX1I1V1zc22AwAAlsnUAeB/JXlFd1+xqHBcf2aS4zezk6p6apJzkpyU5NVJfi3JXya5br74WQSpqgeuqvuKJL+T5FpJnp3kjM20AwAAls3UDwLrDPf+78neyvf84aqHJvmFJH+T5MHd/ck15V+26s9HJHlehrEIJ3f328b1T0tydpJTqurU7hYEAACYhal7AM5P8uCqOnRR4bj+QUn+fX82XlXXSPLLSS5P8rC1J/9J0t2fW/X2lCTXT3LGysn/WOfKJE8d337//rQFAACW0dQB4PQkN01yTlXdq6oOSpKqumZV3SPJ65PcbKy3P05McoskZyX5RFV9a1X9RFX9cFXdZUH9e47LVy8oOydDkDixqg7ez/YAAMBSmfoWoN9Lcrck35XktUk+X1UXJzk6Q9ioJH/a3c/dz+1//bj8rwwzDd1+dWFVnZPklO7+73HVbcblu9ZuqLuvqqr3Jbltkltm6L1YV1Wdu07RcRtrOgAA7LxJewB68PAMU4GeneTSDCf/l47vH97dm3kY2A3G5eOTHJphatHrJLldktdkGOj70lX1jxyXl66zvZX1R22iTQAAsDSm7gFIknT3nyT5ky3Y9Mq0nZXhSv8/j+//raq+PcOV/rtX1V26++82sL2VAcm9t4rdfcLCDQw9A3fYwL4AAGDHTT0GYKt9Yly+d9XJf5L/mWL0NePbO43LlSv8R2axI9bUAwCAA9qyBYALxuUl65SvBISVWYhW6h+7tuI4QPkWSa5K8t6pGggAALvZpgJAVX2+qq6qqmNXvb96A6+r9nOX52Q4Yb91VV1rQfntxuWF4/LscXmfBXVPSnJYkrd092f2sz0AALBUNjsG4JwM989fvub9lujuj1XVSzIMMv7ZfGEu/1TVNyX5lgy386xM+/myDM8NOLWqfmvVg8AOSfLMsc5ztqq9AACw22wqAHT3yXt6v0WelOQbkvxMVZ2U5B8yPFvg2zM88fdx3X3J2J7LqupxGYLAG6rqjCQXJ3lAhilCX5bkJdvQZgAA2BWWbQxAuvuiDAHg2Um+MskPZXjg118luVt3v3RN/TOT3D1D78RDkvxgks9lCBKndveW9VgAAMBuM+k0oFX13iSndfdv7qHOE5I8ubtvub/76e6LM5zAP2mD9d+c5H77uz8AADhQTN0DcPPs/aFaR2W4ZQcAANhmO3EL0LWTfHYH9gsAALO36VuAquqma1YdtWBdMjzF96ZJTol59wEAYEdMMQbgwnzx1J8/PL7WU9ngvfsAAMC0pggAL8wQACrJdyc5L8k7FtS7OsnHk7yuu187wX4BAIB9tOkA0N2PXvlzVX13kld0989vdrsAAMD0Jp0GtLuX7rkCAAAwJ07YAQBgRibtAVhRVV+f5FuS3DjJwQuqdHf/763YNwAAsL6pnwRcSV6Q5BEZBgWvDA5e0avWCwAAALDNpr4F6IlJHpnkRUnumOFk/7QkJyb56SSfTHJGkltOvF8AAGADpr4F6FFJLliZGWjoEMgl3f3WJG+tqtckeWuSv07y/In3DQAA7MXUPQC3SXL2mnX/EzK6+5+SvDLJD0y8XwAAYAOmDgCV5NJV7z+d5Og1dd6d5LiJ9wsAAGzA1AHgQxlm/lnx3iQnrKlz6wzBAAAA2GZTB4B/yBef8L8qyZ2q6mlVdduqekKSB2YYBwAAAGyzqQPAy5Ncs6puMb7/lSTvT/JzSc5L8ltJLknykxPvFwAA2IBJZwHq7jOTnLnq/cVV9XVJHpfkVkkuTPLC7v7IlPsFAAA2ZkueBLxad1+a5FlbvR8AAGDvJr0FqKruW1VT31YEAABMZOqT9b9K8oGq+pWqut3E2wYAADZp6gDw+0kOSfKUJP9cVf9YVU+sqmMm3g/AzqkaXgCwhCYNAN39+CQ3SvKdGaYA/Zokv5HkQ1X1Z1X1gKra8nEHAADAYpPfr9/dn+3ul3b3tyW5SZIfS3JBkgcleUWSD1fVaVPvFwAA2LstHbDb3Rd1969399ck+bokv5nkyCQ/uJX7BQAAFtuWGXuq6tgk35HkwUm+bDv2CQAAfKktux+/qo5KcmqSRyW5U5JKclmSP0zygq3aL8Ck9jTYd1FZ99a1BQAmMGkAGJ8BcN8MJ/33T3KtJJ3kdRlO+v+su6+ccp8AAMDGTd0D8OEk189wtf9dSf4oyQu7+0MT7wdgeyy6or9y5d/VfgCW0NQB4JAkz0vygu5+68TbBgAANmnqAPDl3f2ZibcJAABMZOoHgX3RyX9VXbeqvnLKfQAAAPtv8mlAq+raVfVrVfXRJB9L8r5VZd9QVWdV1R2m3i8AALB3kwaAqjoyyd8l+dEMA4LPzzAgeMW/JLlbku+acr8A26rbAGAAltbUPQA/k+S2SR7d3XdI8tLVhd19eZI3JrnXxPsFAAA2YOoA8OAkr+nuF+6hzvuT3Hji/QIAABswdQC4SZLz9lLnU0mOnHi/AADABkwdAD6Z5AZ7qXOLDIODAQCAbTZ1APjHJN9WVddZVFhVN0pyvyR/O/F+AQCADZg6APxGkmOSnFVVx68uGN+/NMPTgn9z4v0CAAAbMOmTgLv7NVX1jCTPSPKvST6XJFX1sSTXzTAl6E9091um3C8AALAxkz8IrLt/PsM0n3+R5BNJrk7SSc5Kcu/u/tWp9wkAAGzMpD0AK7r79UlevxXbBgAA9t/kPQAAAMDuJQAAAMCMbOoWoKr6fIb7+/dVd/eW3H4EAACsb7Mn4edk/wIAAACwAzYVALr75InaAfAFVV/4c7vGAABT2rLbcKrq8CTHJrl2d79pq/YDAABs3OSDgKvqJlX18gzPAHhbVk0HWlXfWFX/XlUnT71fAABg7yYNAFV1oyR/n+SBSV6Z5O8yPP13xd8nuUGS75xyvwAAwMZM3QPw9Awn+Pfu7gcn+evVhd39uSRvSnLXifcLAABswNQB4H5J/qK737CHOv+Z5Csm3i+wrKq+9LUv5QDAPpk6AHx5knfvpc7nkhw+8X4BAIANmHoWoIuTfOVe6hyb5KMT7xdYVoum+TQNKABsmal7AN6c5AFVdcNFhVV16yT3yaqZgQAAgO0zdQD41SSHJHljVd03yWHJ8EyA8f1fJvl8kl+beL8AAMAGTHoLUHf/fVV9b5LnZpgGdMVl4/KqJI/p7n+bcr8AAMDGTP4k4O5+flX9bZIfSHLnJMckuTTJW5P8dndfMPU+AQCAjZk8ACRJd787yY9uxbaBGTDwFwC2zNRjAAAAgF1MAAAAgBkRAIDN84ReAFgaAgAAAMyIAAAAADMiAAAAwIwIAAAAMCNb8hwA4ACzMsC3e8+DfReVmdMfAHaVpesBqKoLq6rXeX10nc+cWFVnVdXFVXV5VZ1XVT9SVdfc7vYDAMBOWtYegEuTnLZg/afWrqiqByZ5eZIrk7wkycVJ7p/k2UnumuShW9dMOAAtuqK/uocAANjVljUAXNLdz9hbpao6Isnzklyd5OTuftu4/mlJzk5ySlWd2t1nbGVjAQBgt1i6W4D20SlJrp/kjJWT/yTp7iuTPHV8+/070TAAANgJy9oDcHBVPSLJTZN8Osl5Sc7p7qvX1LvnuHz1gm2ck+TyJCdW1cHd/Zktay0sE4N8AeCAtqwB4IZJXrRm3fuq6nu6+42r1t1mXL5r7Qa6+6qqel+S2ya5ZZLz97TDqjp3naLjNtZkAADYect4C9Dzk9wrQwg4PMntk/xekpsneVVVfc2qukeOy0vX2dbK+qOmbyYsgaovvarf/aWvvZWtrQcA7FpL1wPQ3T+3ZtW/Jnl8VX0qyZOTPCPJt29wcytnPns9c+nuExZuYOgZuMMG9wcAADtqGXsA1vPccXnSqnUrV/iPzGJHrKkHAAAHtAMpAFw0Lg9fte6CcXns2spVdVCSWyS5Ksl7t7ZpAACwOxxIAeAu43L1yfzZ4/I+C+qflOSwJG8xAxAAAHOxVAGgqm5bVUcvWH+zJL89vn3xqqKXJflYklOr6o6r6h+S5Jnj2+dsUXNhd1kZ8Lv6tZGyxCBfADiALNsg4Icm+cmqen2S9yX5ZJJbJfnWJIckOSvJs1Yqd/dlVfW4DEHgDVV1RpKLkzwgwxShL0vykm39BgAAsIOWLQC8PsOJ+9dluOXn8CSXJPnbDM8FeFH3F1+m7O4zq+ruSX4myUMyBIX/SPKkJL+5tj4csBYd6itX+v0zAIDZWKoAMD7k6417rfiln3tzkvtN3yIAAFguSzUGAAAA2BwBAAAAZkQAAACAGVmqMQDAxAz+BYDZ0QMAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMCAAAAzIgAALtN1fACANgCAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMH7XQDYNb2NNh3UVn31rUFAJgFAQB2gll+AIAdIgDATlp0RX8lHLjaDwBsAWMAAABgRgQAAACYEQEAAABmxBgA2Gr7OtMPAMAWEgBgtzH4FwDYQgIAbDUz/QAAu4gxAAAAMCMCAAAAzIgAAAAAMyIAAADAjBgEDPtrMwN5Df4FAHaIHgAAAJgRAQAAAGZEAAAAgBkRAAAAYEYMAoaNWBnwu9Eyg3wBgF1KDwCsp2rPJ/4AAEtIDwBsxKIr+puZBhQAYIfoAQAAgBkRAAAAYEYEAAAAmBFjACAxyw8AMBsCAOwvIQAAWEICACRm+QEAZsMYAAAAmBEBAAAAZkQAAACAGREAAABgRgwChvUY/AsAHID0AAAAwIwIAAAAMCMCAAAAzIgAAAAAMyIAAADAjAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMCAAAAzIgAAAAAMyIAAADAjAgAAAAwIwIAAADMyNIHgKp6ZFX1+HrsOnW+rareUFWXVtWnqurvq+pR291WAADYaUsdAKrqK5P8VpJP7aHOE5P8ZZLbJXlxkucl+YokL6iqZ21HOwEAYLdY2gBQVZXk+Uk+nuS569S5eZJnJbk4yR27+wnd/aNJvjrJe5I8uarusi0NZn1VwwsAgC23tAEgyQ8luWeS70ny6XXqPCbJwUl+u7svXFnZ3Z9I8kvj28dvYRsBAGBXWcoAUFXHJ/m/SX6ju8/ZQ9V7jstXLyh71Zo6AABwwDtopxuwr6rqoCQvSvKfSX56L9VvMy7ftbaguz9SVZ9OcpOqOqy7L9/Lfs9dp+i4vbQBAAB2jaULAEl+NsnXJfnG7r5iL3WPHJeXrlN+aZLDx3p7DAAAAHAgWKoAUFV3ynDV/9e6+++m2OS47L1V7O4T1mnTuUnuMEFb5mFPg30XlfVe/2oAANgHSzMGYNWtP+9K8rQNfmzlyv+R65QfMS4v20TTAABgaSxNAEhy7STHJjk+yZWrHv7VSZ4+1nneuO608f0F4/LYtRurqhtluP3ng3u7/58JdX/payNlAABMYpluAfpMkj9cp+wOGcYF/G2Gk/6V24POTnLXJPdZtW7FfVfVAQCAWViaADAO+H3sorKqekaGAPBH3f0Hq4qen+THkzyxqp6/8iyAqrpuvjCD0MKHiAEAwIFoaQLA/uju91XVjyX5zSRvq6qXJPlsklOS3CTTDSYGAIClcEAHgCTp7t+qqguTPCXJd2cY9/DvSZ7a3X+0k20DAIDtdkAEgO5+RpJn7KH8L5P85Xa1h31ksC8AwLZZplmAAACATRIAAABgRgQAAACYEQEAAABmRAAAAIAZEQAAAGBGBAAAAJgRAQAAAGZEAAAAgBkRAAAAYEYEAAAAmBEBAAAAZkQAAACAGREAAABgRgQAAACYEQEAAABmRAAAAIAZEQAAAGBGBAAAAJgRAQAAAGZEAAAAgBkRAAAAYEYEAAAAmBEBAAAAZkQAAACAGREAAABgRgSAOagaXgAAzJ4AAAAAMyIAAADAjAgAAAAwIwIAAADMyEE73QAmtqfBvovKureuLQAA7Dp6AAAAYEb0AByoVl/ZX7ny72o/AMDs6QEAAIAZEQAAAGBG3AK0zPZ1wC8AALOnB2BZOcEHAGA/6AFYdmsH9hrwCwDAHugBAACAGREAAABgRgQAAACYEWMAloHZfgAAmIgegAORAcAAAKxDD8AyWHRCb7YfAAD2gx4AAACYEQEAAABmRAAAAIAZEQAAAGBGDAJeVgb/AgCwH/QAAADAjAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAAAAMCMCAAAAzIgAAAAAM1LdvdNtWGpV9fFDDz306OOPP36nmwIAwAHs/PPPzxVXXHFxdx+zme0IAJtUVe9LckSSC5McN65+5441iGXhWGEjHCdslGOFjXKsLLebJ7msu2+xmY0IABOqqnOTpLtP2Om2sLs5VtgIxwkb5VhhoxwrJMYAAADArAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMmAUIAABmRA8AAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAAAAADMiAAAAwIwIAHtRVTepqtOr6sNV9ZmqurCqTquq6+7jdo4eP3fhuJ0Pj9u9yVa1ne2z2eOkqg6vqodX1R9X1Tur6tNV9cmqeltVPbmqrrXV34HtMdXvlDXbPKmqrq6qrqpnTtledsaUx0lV3b6qXlhVHxi3dVFVvbGqvnsr2s72mvA85Rur6s/Hz19ZVf9ZVWdV1X22qu3sHA8C24OqulWStyS5QZI/T/LOJHdKco8kFyS5a3d/fAPbOWbczrFJzk7yj0mOS/LAJBcluUt3v3crvgNbb4rjZPwF+6okFyd5fZL/SHJ0kvsnueG4/Xt195Vb9DXYBlP9TlmzzeskOS/J9ZJcO8kvdvdTp2w322vK46SqHp3kD5JcnuSVSS5MclSS2yX5cHf2FOkSAAAOMklEQVSfOnHz2UYTnqd8f5LfTfLpJK9I8sEkN0ny4CSHJXlqd//iVnwHdkh3e63zSvKaJJ3kB9es//Vx/XM3uJ3fG+v/+pr1PzSuf/VOf1evnT1Oknxtkocnudaa9ddJcu64nSfv9Hf12vljZcE2T88QHH963MYzd/p7eu2O4yTJnZNcleQdSW64oPzLdvq7eu38sZLky5JckuSKJLdZU3Z8kiszBMiDd/r7ek330gOwjqq6ZZL3ZLhacqvu/vyqsusk+UiSSnKD7v70HrZzeJL/TvL5JDfq7k+uKrvGuI+bj/vQC7BkpjpO9rKPhyX5f0le2d3333Sj2RFbcaxU1QOTnJnkkUkOSvL86AFYalMeJ1V1TpK7Jbl9d//rljWaHTHhecqXJ/lokvO6+2sWlJ+X5PZJrtf72EPJ7mUMwPruOS5fu/ofVZKMJ/FvztAtdue9bOcuSQ5N8ubVJ//jdj6f5LXj23tsusXshKmOkz353Li8ahPbYOdNeqxU1Q2SPC/Jmd394ikbyo6a5DgZx5fdLcnbkvxbVd2jqp4yjim613gBiuU21e+UizJcqDy2qm69uqCqjk1y6yTvcPJ/YPELYH23GZfvWqf83ePy2G3aDrvTdvz9PmZcvnoT22DnTX2s/H6G3+GP30yj2HWmOk6+flX9s8fXryZ5VpK/SfKOqvqqTbSTnTfJsdLDrSBPyPD75Nyq+qOq+j9V9cIMt6D+W5KHTtBedpGDdroBu9iR4/LSdcpX1h+1Tdthd9rSv9+qemKS+2S4h/f0/dkGu8Zkx0pVPSbDJALf2d3/NUHb2D2mOk5uMC6/I8nHMgzmfF2S6yd5eobbxv6qqm7f3Z/d/+aygyb7ndLdL62qDyf5kySrZ4f6rwy3FrpF+QCjB2D/1bjc7CCKqbbD7rTff79V9eAkp2W4N/Mh3f25vXyE5bahY6Wqbp7huHhpd//pFreJ3Wejv1OuuWr52O5+RXdf1t3vSfKoDLcGHZvkIVvTTHaBDf//U1WPyNAz9KYMA38PG5evS/LbSc7YojayQwSA9a0k5yPXKT9iTb2t3g6705b8/VbVgzL8wr0oyckGiB8QpjpWTs8wW8cPTNEodp2pjpNPjMvPJDlrdcF4y8efj2/vtK8NZNeY5FgZ7/M/PcOtPo/s7nd29xXd/c4MPUXnJnloVZ28+SazWwgA67tgXK5379zKQJn17r2bejvsTpP//VbVQ5O8NEPX6927+4K9fITlMNWxcocMt3f89/jgr66qztBNnyQ/M647c3PNZYdM/X/PJ9cOEB2tBIRD96Ft7C5THSvfnGEq0DcuGEz8+STnjG9P2J9GsjsZA7C+14/Lb66qayyYXuuuGa7CvXUv23nrWO+uVXWdBdOAfvOa/bFcpjpOVj7zsCQvTPKhJPdw5f+AMtWx8sIM3fNr3TrJSRnGi5yb5J823WJ2wlTHyXkZ7v2/XlV9+YKxIrcblxduvsnskKmOlYPH5fXXKV9Zb6zIAUQPwDrG+yRfm2GO/iesKf65JIcneeHquXWr6riqOm7Ndj6V5EVj/Wes2c4Tx+2/xonecprqOBnXPyrDsfKfSU5yTBxYJvyd8kPd/di1r3yhB+CvxnW/s2Vfhi0z4XFyVYaHUCbJr6ye9rOqbp/k0RmmFn7ZxF+BbTLh/z9vGpenVNVXry6oqq9NckqGcQRnT9d6dpoHge3Bgkdsn5/kGzLM2f+uJCeunhd37IZPd9ea7RwzbufYDP+A/iHD4JoHZrjH+8TxHzJLaIrjpKrukWEA1jUy3Iv5gQW7uqS7T9uir8E2mOp3yjrbfnQ8COyAMOH/PYdlGMR55ww9Qm/IcDX3IRlu/Xlyd//6Fn8dttCEx8rpSb4nw1X+VyR5f4Zg8aAk10pyWnf/6BZ/HbaRALAXVfWVSX4+w1SMx2R4st6ZSX6uuy9eU3fd/6yr6ugMU689KMmNknw8yauS/Gx3f3ArvwNbb7PHyaqTtz15f3fffLpWsxOm+p2yYLuPjgBwwJjw/57Dkvx4klOT3CLJlUn+McmvdfertvI7sD2mOFaqqjLMDvXoJF+T5DpJLssQHJ/X3WYBOsAIAAAAMCPGAAAAwIwIAAAAMCMCAAAAzIgAAAAAMyIAAADAjAgAAAAwIwIAAADMiAAAAAAzIgAAAMCMCAAAADAjAgAAAMyIAABwgKiqF1bVRVV1+Ca3c2FVXThRs3bMTn2PqnpyVX2uqo7b7n0DbIQAALDLVFWveV1dVR+rqrOr6uHrfOaOSR6R5P9296e3t8U7o6reUFW90+1Y4HeTXJTkWTvdEIBFDtrpBgCwrp8bl1+W5DZJHpTkHlV1Qnc/aU3dX0pyWZLnbGP7drt77cROu/uKqvqNJL9cVSd291t2oh0A66nu3XjxBGC+Vq5qd3etWX+vJH89vr1ld184rj82yTuT/EF3f+8E+79w3P/NN7utrVRVb0hy97U/p92gqr4iyX8mOaO7H7HT7QFYzS1AAEuiu1+X4US/knz9qqLHjOtest5nq+pbquqs8Vaiz1TVe6rqV6vqqH1pQ1V9V1W9vqo+UVVXVtX5VfXUqjp4VZ0bj7ctvX0P23n1eHvT7Vate3RVvbyq3ltVV1TVZVX15qp6xJrP3nwMSXcf36++XeoNq+otHANQVQdX1U9W1XlVdfm4nzdV1XcsqHvzcbsvGP98xvgzvLKq3lZV37bo+3X3h5O8KckpVXXEHn6kANvOLUAAy2Xlavfq7tt7J7k6yVsXfqDqZzPcTnRxkldmuD/9q5M8Jcn9quou3X3ZXndc9YcZwsYHk/xZkkuS3DnJLyS5V1V9U3df1d0fqqq/SfLNVXX77v6XNdu50djmc7v7X1cVPSfJvyc5J8lHkhyT5H5JXlRVt+nup431Lhm/z6OT3CxfuFUqSS7cy3e4VpLXZAgP70zyO0kOS3JKkpdU1dd2908v+OjNkvxDkvcmeVGSo5N8Z5I/r6p7d/frF3zmzUlOTnJShp87wO7Q3V5eXl5eu+iV4eS+F6y/d5LPj6+bjesOT3JVkn9ZZ1v3GLf3liRHrSl79Fj27DXrL0xy4Tp1/yzJoWvKnjGW/fCqdd81rnvWgjb92Fj2g2vW32pB3WsleV2SzyW58ZqyNyz6Oe3le/zUuO+zkhy0av0Nxvqd5MRV62++8veR5OlrtvUtK9taZ/8PHMt/ZaePKS8vL6/VL7cAAexSVfWM8fWLVfWyJK/O0ANwWne/f6x24yTXzHDFfJEfGpeP6+5LVhd09wuSvCPJwpmF1vjhDEHjMd19xZqyX0jy8TXbOTPJpUkeXlXXXFP/URlO6P9kTXves3an3f3ZDFfpD8o0g3ofk+Gk/EndfdWq/VyU4XskyWMXfO79SZ65pm2vyXCf/53W2ddHx+VNN9NggKm5BQhg93r6uOwMt728KckfdveLV9U5Zlx+Yp1t3CXDyfZDq+qhC8qvleT6VXVMd3980Qaq6rAkX5PkY0l+pGrhmNvPJDl+5U0PM+H8aZLHZbhSfta4rROS3DbJK7r7Y2v2c9MkP5HhRP+mSQ5ds48br/MdN6SqrpPkq5J8qLvfuaDK2ePy6xaUvaO7r16w/gMZfsaLXDwur7dPDQXYYgIAwC7VG5vdZuVq/CHrlB+T4Xf909cpX3HtDFfxF7luhp6H629gO6u9IEMAeFTGADD+OUn+aHXFqrplhnvsr5sh6Lw2Qw/C1Rluw3lUkoOzOUeOy/V6S1bWLxoYfcmCdcnQK7Jeb/pKgFnbYwKwowQAgOV20bg8Zp3yS5Nco7uP3sQ+Lh2X/9Tdd9joh7r7LVX17iQPHGcb+nSGsQEfyxcCwYonZfgO3zPemvQ/quq78oXgsBkr3+OG65TfaE29zVr5O7loj7UAtpkxAADL7SNJ/jvDg8IWeWuS61bVbfd3B939qST/luS2VbWvQeKPMly5/84k35rhdpg/7u7Pran3VePy5Qu2cfd1tn11kiwYY7BQd38yyXuS3Liqbr2gyj3G5brTl+6j48blOybaHsAkBACAJdbdnWHazOtV1VctqPLscfm88eFUX6SqDq+qO29gV7+eYbzA6YueHVBV162qRb0DL8wwa9F3j69kuDVorQvH5clrtvstWTwoN/nCLUv7Msj29Ay3M/3q6uBQVddL8rRVdaaw8nNdNEUowI5xCxDA8nt5kodkGGz7H6sLuvt1VfWTSf5PkndX1VlJ3pfhnv+bZbi6/rdJ7rOnHXT36eMA3h9I8p6qWpkB5+gkt8gw1/3zkzx+zec+UFWvzzCwd2W60n9asIvfTfI9SV5aVS9P8qEktxvb9acZehDWel2Shyb5s/F7XZHk/d39oj18lWcluW+GKTr/efzcYeN2bpBhys6/3dPPYiOq6hoZvvMF/cXPOgDYcQIAwPJ7eZL/ynCF/XfWFnb3L1fVmzNMCfqNGU5+L81wkv37Sf54Izvp7idU1asynOTfO8Ng2YszBIFfTfLidT76ggwnwwdlzeDfVds+r6rukWGqzfuNdf85yYMzDMBdFAD+IEOIOTXJj4+feWOGB3Wt9x0+W1XflGHMwcOS/GCGYPLPSX6ku/9kvc/uo3tnmLXoRyfaHsBkaug9BmCZVdVPJfmlJHdY5wo722jsxbh7hoebTTWoGGASAgDAAaCqDklyQZLzuvv+O92eOauqr80wkPiHuvu3d7o9AGsZBAxwAOjuK5M8MsnbqurwnW7PzN0ow4Di5+50QwAW0QMAAAAzogcAAABmRAAAAIAZEQAAAGBGBAAAAJgRAQAAAGZEAAAAgBkRAAAAYEYEAAAAmBEBAAAAZkQAAACAGREAAABgRgQAAACYEQEAAABmRAAAAIAZ+f8uoR9bZs6zQAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 263,
       "width": 384
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng1 = np.random.RandomState(15)\n",
    "el1 = []\n",
    "prob1 = []\n",
    "air1 = []\n",
    "for _i in range(20):\n",
    "    vals = ramp_distribution(rng1, 30., 70., 85.)\n",
    "    _elevation = vals[1]\n",
    "    _airmass = 1./np.cos((90. - _elevation)*np.pi/180.)\n",
    "    prob1.append(vals[0])\n",
    "    air1.append(_airmass)\n",
    "    el1.append(_elevation)\n",
    "plt.plot(prob1, el1,'r+')\n",
    "plt.xlabel('P(elevation)')\n",
    "plt.ylabel('elevation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quad_distribution(rng, x0, x1, x2, invert=True):\n",
    "    # As the ramp_distribution, but now the probability distribution increases quadratically from 0 between x0 and x1\n",
    "    c1 = 1./((x1**3. - x0**3.)/3. + x1**2.*(x2 - x1))\n",
    "    p1 = c1*(x1**3. - x0**3.)/3.\n",
    "    if invert:\n",
    "        # Calculate x from P(x)\n",
    "        p = rng.rand()\n",
    "        if p <= p1:\n",
    "            x = (3.*p/c1 + x0**3.)**(1./3.)\n",
    "        else:\n",
    "            x = (p - p1)/(c1*x1**2.) + x1\n",
    "        return((p, x))\n",
    "    else:\n",
    "        # Calculate P(x) from x\n",
    "        x = rng.rand()*(x2 - x0) + x0\n",
    "        if x <= x1:\n",
    "            p = c1*(x**3. - x0**3.)/3.\n",
    "        else:\n",
    "            p = p1 + c1*x1*2.*(x - x1)\n",
    "        return((x, p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'elevation')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwAAAAIPCAYAAADEoLnvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu8LWVdP/DPV0huCghqmqaoySW1EswUlYtaKaWYYpKXJH9almamdlVTyy7eyS5aKir6SwwM+mWIpoAopiWmlCGacLwb4RFQucjl+f0xs2O7Xfvsfc6etddee97v12tew5rnWTPPOgyH+cwzzzPVWgsAADAON5l1AwAAgPUjAAAAwIgIAAAAMCICAAAAjIgAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAAAjIgAAAMCICAAAADAiAgAAAIyIAAAAACMiAAAAwIgIAAAAMCICAAAAjMjOs27AvKuqi5PsmWTLjJsCAMDmtl+SK1prd1rLTgSAtdtzt9122+eggw7aZ9YNAQBg87rgggty1VVXrXk/AsDabTnooIP2Oe+882bdDgAANrFDDjkkH/vYx7asdT/GAAAAwIgIAAAAMCICAAAAjIgAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAAAjIgAAAMCICAAAADAiAgAAAIyIAAAAAKtV1S1zTAAAAIAREQAAAGBEBAAAABgRAQAAAEZk51k3AAAANqRtDfadVNba9NoyID0AAAAwInoAAABgkkl39Bfu/M/J3f5J9AAAAMCICAAAADAiAgAAAIyIAAAAACNiEDAAAKzWHA/+XaAHAAAARkQAAACAEREAAABgRAQAAAAYkbkMAFX1U1X1nqr6YlVdVVUXVdXJVXXfZeofWlWnV9XWqrqyqs6vqmdW1U7r3XYAALZD1Y1v32UQcxcAquolSd6Z5OAkZyT50yQfS3J0knOr6vFL6h+d5JwkhyU5NclfJLlpklclOWn9Wg4AALM3V9OAVtVtkjwnyX8n+aHW2iWLyo5McmaS30/y1n7bnklel+T6JEe01j7ab39+X/eYqjq2tSYIAAAwCvPWA3DHdG3+yOKL/yRprZ2V5BtJbrVo8zH955MWLv77ulcneV7/8Zen2mIAANhA5i0AfCbJt5Pcu6puubigqg5LcvMk7120+YH9+owJ+zonyZVJDq2qXabQVgAA2HDm6hGg1trWqvqtJK9M8p9VdVqSryW5S5KHJ/mnJL+06CsH9OtPT9jXdVV1cZK7Jblzkgum2XYAAHLjgN5Jb9Td1mDfSWWb4K28szBXASBJWmvHV9WWJCckecqiov9K8qYljwbt1a8vX2Z3C9v3Xum4VXXeMkUHrvRdAADYKObtEaBU1W8mOSXJm9Ld+d8jySFJLkryf6vqpduzu34tPgIAzFpr372spoztMlc9AFV1RJKXJDm1tfasRUUfq6qfSfeoz7Or6rWttYty4x3+vTLZnv16uR6C/9VaO2SZNp2XbkpSAADY8OatB+Cn+/VZSwtaa1cm+Zd0v+me/eYL+/X+S+tX1c5J7pTkunS9BwAAsOnNWwBYmK3nVsuUL2z/dr8+s18/ZELdw5LsnuRDrbVrhmkeAAD/a+EtvouX1ZQxVfMWAD7Qr3+xqm63uKCqHprkfkmuTvKhfvMpSS5NcmxV3WtR3V2TvLj/+JqpthgAADaQuRoDkO6C/r1JHpzkgqo6NclXkxyU7vGgSvLbrbWvJUlr7Yqqekr/vbOr6qQkW9NNGXpAv/3t6/4rAADGYFtTfa52EK/BvoObqwDQWruhqo5K8rQkxyb5mXSP8WxNcnqSV7fW3rPkO6dV1eFJnpvkUUl2TTdl6LP6+s4qAABGY64CQJK01q5Ncny/rPY75yY5amqNAgCAOTFvYwAAAIA1EAAAAGBE5u4RIAAA5pjhlzOnBwAAYMzMwT86AgAAAIyIAAAAACMiAAAAwIgIAAAAMCJmAQIAGIttDfadVGbGnk1JDwAAAIyIHgAAgLGYdEd/4c6/u/2joQcAAABGRAAAAIAREQAAAGBEBAAAABgRg4ABAMbM4N/R0QMAAAAjIgAAAMCICAAAADAiAgAAAIyIAAAAACMiAAAAwIgIAAAAMCICAAAAjIgAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAAAjIgAAAMCICAAAADAiAgAAAIyIAAAAkCRV3QKbnAAAAAAjIgAAAMCICAAAADAiAgAAAIzIzrNuAADAutvWYN9JZa1Nry2wzvQAAADAiOgBAADGZ9Id/YU7/+72s8npAQAAgBERAAAAYEQEAAAAGBEBAAAARsQgYACAxOBfRkMPAAAAjIgAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAAAjIgAAAMCICAAAADAiAgAAAIyIAAAAACMiAAAAwIgIAAAAMCICAAAAjIgAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAAAjIgAAAMCICAAAADAiAgAAAIyIAAAAACMiAAAAwIgIAAAAMCICAAAAjIgAAAAAIyIAAADAiMxVAKiq46qqrbBcP+F7h1bV6VW1taqurKrzq+qZVbXTLH4HAADMys6zbsB2+niSFy1T9oAkD0zyrsUbq+roJO9IcnWStyfZmuRhSV6V5H5JHj2txgIAwEYzVwGgtfbxdCHgu1TVP/f/+NeLtu2Z5HVJrk9yRGvto/325yc5M8kxVXVsa+2kqTYcAAA2iLl6BGg5VXX3JPdJ8qUk/7io6Jgkt0py0sLFf5K01q5O8rz+4y+vVzsBAGDWNkUASPJL/foNrbXFYwAe2K/PmPCdc5JcmeTQqtplmo0DAICNYq4eAZqkqnZL8vgkNyR5/ZLiA/r1p5d+r7V2XVVdnORuSe6c5IIVjnPeMkUHbleDAQBghjZDD8DPJtk7ybtaa19YUrZXv758me8ubN97Gg0DAICNZu57AJL8Yr/+qx34bvXrtlLF1tohE3fQ9QwcvAPHBgCAdTfXPQBV9YNJDk3yxSSnT6iycId/rwllSbLnknoAALCpzXUAyPKDfxdc2K/3X1pQVTsnuVOS65JcNJ3mAQDAxjK3AaCqdk3yhHSDf9+wTLUz+/VDJpQdlmT3JB9qrV0zfAsBAGDjmdsAkO4NvrdIcvqEwb8LTklyaZJjq+peCxv78PDi/uNrptpKAADYQOZ5EPDC4N+/Xq5Ca+2KqnpKuiBwdlWdlGRrkoenmyL0lCRvn3ZDAQBgo5jLHoCqOijJ/bP84N//1Vo7Lcnh6V789agkv5rk2iTPSnJsa23FGYAAAGCzmMsegNbaBblxCs/V1D83yVHTaxEAAMyHuewBAAAAdowAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAFXdAjACAgAAAIyIAAAAACMiAAAAwIgIAAAAMCI7z7oBALCutjXYd1JZa9NrC8AM6AEAAIAR0QMAwLhMuqO/cOff3X5gBPQAAADAiAgAAAAwIgIAAACMiAAAAAAjYhAwABj8C4yIHgAAABgRAQAAAEZEAAAAgBERAAAAYEQEAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAAAAgBERAAAAYEQEAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABiRnae146raI8neSXaaVN5a+/y0jg0AAEw2eACoqick+a0kB22jWpvGsQEAgG0b9CK8qo5LckKS65N8IMkXklw35DEAAIAdN/Rd+Ock+XqS+7fWLhh43wAAwBoNPQj4B5Kc4uIfAAA2pqEDwNYkVw+8TwAAYCBDB4B3Jjmiqmrg/QIAAAMYOgD8TpJdkry2qm428L4BAIA1GnoQ8MlJrkzy5CSPrarPJLlsQr3WWnvQwMcGAABWMHQAOGLRP++R5EeWqdcGPi4AALAKgwaA1trQjxQBAAADcsEOAAAjIgAAAMCITCUAVNWxVfXeqvpaVV1XVVur6p+q6thpHA8AAFidQccA9PP/n5jksUkqyfVJ/ifJLZM8KMkDq+rhrbXHDnlcAABgdYbuAfilJI9L8rEkD06ya2vttkl27T+fl+QxVfXUgY8LAACswtAB4ElJtiQ5rLV2Zmvt+iRprV3fWjszyeF9+f8Z+LgAAMAqDB0AfjDJqa21qyYV9ttPS3LQwMcFAABWYegA0NI9+78tK5UDAABTMnQAuCDJI6tqt0mF/fZHJPnPgY8LAACswtAB4IQkd0hyTlU9qKp2TpKq2qmqjkxyVpI79vUAAIB1Nug0oEn+KskDkvxckvckuaGqtibZJ13YqCR/21p77cDHBQAAVmHQHoDWeVy6qUDPTHJ5uov/y/vPj2uteRkYAADMyNA9AEmS1trbkrxtGvsGAAB23NBjAAAAgA1MAADYXlXdAgBzaE2PAFXVDUluSPKDrbVP95/bKr7aWmtTefwIAABY3lovws9Jd8F/5ZLPAADABrSmANBaO2JbnwEAgI1lbscAVNUDquodVfWVqrqmX7+nqo6aUPfQqjq9qrZW1ZVVdX5VPbOqdppF2wEAYFYGDQBVdVFVPWOFOk+rqovWeJznpXvc6LAkZyR5RZJ/SHKLJEcsqXv0orqnJvmLJDdN8qokJ62lHcAILAz4XbyspgwANqihB+Lul2TvFersneSOO3qAqnp0kj9I8t4kj2ytfWNJ+fcs+uc9k7wuyfVJjmitfbTf/vx0LyY7pqqOba0JAgAAjMIsHgG6WZJv78gXq+omSV6SbtDxY5de/CdJa+3aRR+PSXKrJCctXPz3da5O8rz+4y/vSFuAkWjtu5fVlAHABrXmHoCqusOSTXtP2JYkOyW5Q7qL8h19BOjQJHdKckqSr1fVTyW5e5Krk/xLa+2fl9R/YL8+Y8K+zkkXJA6tql1aa9fsYJsAAGBuDPEI0JZ859Sfv9Yvy6kkz9rBY/1ov/7vJB9Lco/v2HHVOUmOaa39T7/pgH796aU7aq1dV1UXJ7lbkjsnuWBbB66q85YpOnB1TQcAgNkbIgCcmC4AVJKfT3J+ko9PqHd9kq8leV9r7T07eKxb9+unJrk4yYOTfCTdmIJXJPnJJCfnxoHAe/Xry5fZ38L2lcYtAADAprDmANBaO27hn6vq55Oc2lr7/bXudxkL03ZWujv9n+g/f7Kqfibdnf7Dq+q+Ex4HmmRhyo4VH9xtrR0ycQddz8DBqzgWAADM3KCDgFtrN5nixX+SfL1fX7To4n/h2FcleXf/8d79euEO/16ZbM8l9QBWZsAvAHNs3l4EdmG/vmyZ8oWAsNuS+vsvrVhVO6cbUHxddnxQMgAAzJWh3wOQJKmqH033PP7tkuwyoUprrf2fHdj1Oeku2O9aVTdtrS2dTvTu/XpLvz4zyeOSPCTJ25bUPSzJ7knOMQMQAABjMWgAqKpK8qYkj0/3fP3C4OAFbdH27Q4ArbVLq+rt6S7qfy83zuWfqvrxdKHj8tw47ecp6d4bcGxV/dmiF4HtmuTFfZ3XbG87AABgXg39CNDTkzwhyVuS3Cvdxf7x6ebv/90k30hyUrppN3fUs5L8V5LnVtU5VfXyqjo5ybvSzTT0lNbaZUnSWrsiyVPSDR4+u6peX1UvTTdL0X3TBYS3r6EtAAAwV4Z+BOiJSS5cmBmo6xDIZa21Dyf5cFW9O8mHk/xTkjfuyAFaa5dU1Y+lu/v/M0nuky5Y/GOSP+6Ptbj+aVV1eJLnJnlUkl3TBYhnJXl1a0byAQAwHkMHgAPSvRdg4jFaa/9WVe9M8ivZwQDQ72drugv4Vb1QrLV2bpKjdvR4AACwWQz9CFDlO6fU/FaSfZbU+Uy8PRcAAGZi6ADwpXQz/yy4KMnSF2jdNV0wAAAA1tnQAeBf8p0X/O9Kcu+qen5V3a2qnpbk6HTjAAAAgHU2dAB4R5KdqupO/eeXJvlckhclOT/Jn6V7iddvD3xcAABgFQYdBNxaOy3JaYs+b62qe6abivMu6V7QdWJr7StDHhcAAFidqbwJeLHW2uVJXj7t4wAAACsb9BGgqnpoVQ39WBEwRlXdAgAMauiL9X9M8oWqemlV3X3gfQNj4cIfAKZm6ADw1+netPucJJ+oqn+tqqdX1b4DHwcAANgBgwaA1tpTk9w2yWPSTQH6w0n+NMmXqurvqurhVTX1cQcAAMBkgz+v31r7dmvt5NbaTye5fZLfSHJhkkckOTXJl6vq+KGPCwAArGyqA3Zba5e01l7ZWvvhJPdM8uokeyX51WkeF5gzCwN+Jw38XVpmfAAArMm6zNhTVfsn+dkkj0zyPetxTGBOuKAHgHU1tefxq2rvJMcmeWKSeyepJFckeUOSN03ruMAm19qsWwAAc23QANC/A+Ch6S76H5bkpklakvelu+j/u9ba1UMeE5hzky7oF/cKuOAHgEEN3QPw5SS3Sne3/9NJ3pzkxNbalwY+DgAAsAOGDgC7Jnldkje11j488L4BAIA1GjoAfG9r7ZqB9wkAAAxk6BeBfcfFf1Xdoqq+f8hjACPi+X8AGNzg04BW1c2q6hVV9dUklya5eFHZj1XV6VV18NDHBTaR1lz8A8CUDBoAqmqvJP+c5NfTDQi+IN2A4AX/nuQBSX5uyOMCAACrM3QPwHOT3C3Jca21g5OcvLiwtXZlkvcnedDAxwVmyRt6AWBuDB0AHpnk3a21E7dR53NJbjfwcQEAgFUYOgDcPsn5K9T5ZpK9Bj4uAACwCkMHgG8kufUKde6UbnAwAACwzoYOAP+a5Ker6uaTCqvqtkmOSvLBgY8LAACswtAB4E+T7Jvk9Ko6aHFB//nkdG8LfvXAxwWmafEg34V/XrwsrTepDADYEAZ9E3Br7d1V9cIkL0zyH0muTZKqujTJLdJNCfpbrbUPDXlcAABgdQZ/EVhr7ffTTfP5/5J8Pcn1SVqS05M8uLX2sqGPCayjhZd0LV5WUwYAbAiD9gAsaK2dleSsaewbAADYcYP3AAAAABvXVHoAgDm2rYG7k8o85gMAc2VNAaCqbkj3fP/2aq014QMAANbZWi/Cz8mOBQBgo5p0R3/hzv9yd/v1AgDA3FhTAGitHTFQO4BZWOnCHgDYdKb2GE5V7ZFk/yQ3a619YFrHAQAAVm/wWYCq6vZV9Y507wD4aBZNB1pV96+q/6yqI4Y+LgAAsLJBA0BV3TbJR5IcneSdSf453dt/F3wkya2TPGbI4wIAAKszdA/AC9Jd4D+4tfbIJP+0uLC1dm2SDyS538DHBabJW30BYNMYOgAcleT/tdbO3kadzyf5voGPC6yk6ruX1ZQBAJvK0AHge5N8ZoU61ybZY+DjAgAAqzD0LEBbk3z/CnX2T/LVgY8LrGRH5vcHADadoXsAzk3y8Kq6zaTCqrprkodk0cxAAADA+hk6ALwsya5J3l9VD02ye9K9E6D//A9JbkjyioGPCwAArMKgjwC11j5SVb+Y5LXppgFdcEW/vi7Jk1prnxzyuAAAwOoM/ibg1tobq+qDSX4lyX2S7Jvk8iQfTvLnrbULhz4mAACwOoMHgCRprX0mya9PY9/AgAz+BYDRGXoMAAAAsIEJAAAAMCICAAAAjIgAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAAAjIgDARlPVLQAAUyAAAADAiAgAAAAwIgIAAACMiAAAAAAjsvOsGwCjtq3BvpPKWpteWwCAURAAYBbM8gMAzIgAALM06Y7+Qjhwtx8AmAJjAAAAYEQEAAAAGBEBAAAARsQYAJi27Z3pBwBgigQA2GgM/gUApkgAgGkz0w8AsIEYAwAAACMiAAAAwIjMXQCoqi1V1ZZZvrrMdw6tqtOramtVXVlV51fVM6tqp/VuPwAAzNK8jgG4PMnxE7Z/c+mGqjo6yTuSXJ3k7Um2JnlYklcluV+SR0+vmQAAsLHMawC4rLX2wpUqVdWeSV6X5PokR7TWPtpvf36SM5McU1XHttZOmmZj4bsY/AsAzMjcPQK0nY5JcqskJy1c/CdJa+3qJM/rP/7yLBrGJlBlHn8AYO7Maw/ALlX1+CR3SPKtJOcnOae1dv2Seg/s12dM2Mc5Sa5McmhV7dJau2ZqrQUAgA1iXgPAbZK8Zcm2i6vqF1pr71+07YB+/emlO2itXVdVFye5W5I7J7lgKi0FAIANZB4DwBuTfCDJJ5N8I93F+9OT/GKSd1XVfVtrn+jr7tWvL19mXwvb917poFV13jJFB66m0QAAsBHMXQBorb1oyab/SPLUqvpmkmcneWGSn1nl7hYe4DYiEwCAUZi7ALANr00XAA5btG3hDv9e3109SbLnknrLaq0dMml73zNw8CrbyLza1mDfSWVm+QEANqjNNAvQJf16j0XbLuzX+y+tXFU7J7lTkuuSXDTdpgEAwMawmQLAffv14ov5M/v1QybUPyzJ7kk+ZAYgVtTady+rKQMA2GDmKgBU1d2qap8J2++Y5M/7j29dVHRKkkuTHFtV91pUf9ckL+4/vmZKzWXemecfANiE5m0MwKOT/HZVnZXk4nSzAN0lyU8l2TXJ6UlevlC5tXZFVT0lXRA4u6pOSrI1ycPTTRF6SpK3r+svAACAGZq3AHBWugv3e6Z75GePJJcl+WC69wK8pbXvfP6itXZaVR2e5LlJHpUuKPxXkmclefXS+gAAsJnNVQDoX/L1/hUrfvf3zk1y1PAtAgCA+TJXAQA2FJ1HAMAcEgAgMc8/ADAaczULEAAAsDZ6ACCZfEd/4c6/u/0AwCaiBwAAAEZEAAAAgBERAAAAYEQEAAAAGBGDgGE5Bv8CAJuQHgAAABgRAQAAAEZEAAAAgBERAAAAYEQEAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAAAAgBERAAAAYEQEAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAGD2qroFAICpEwAAAGBEBAAAABgRAQAAAEZEAAAAgBHZedYNYGS2Ndh3Ullr02sLAMAI6QEAAIAR0QPA+pp0R3/hzr+7/QAAU6cHAAAARkQAAACAEREAAABgRAQAAAAYEYOAmT2DfwEA1o0eAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAAAAgBERAAAAYEQEAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAAAAgBERAAAAYETmPgBU1ROqqvXLk5ep89NVdXZVXV5V36yqj1TVE9e7rQAAMGtzHQCq6vuT/FmSb26jztOT/EOSuyd5a5LXJfm+JG+qqpevRzsBAGCjmNsAUFWV5I1JvpbktcvU2S/Jy5NsTXKv1trTWmu/nuSHknw2ybOr6r7r0mAAANgA5jYAJHlGkgcm+YUk31qmzpOS7JLkz1trWxY2tta+nuSP+o9PnWIbAQBgQ5nLAFBVByX5kyR/2lo7ZxtVH9ivz5hQ9q4ldQAAYNPbedYN2F5VtXOStyT5fJLfXaH6Af3600sLWmtfqapvJbl9Ve3eWrtyheOet0zRgSu0AQAANoy5CwBJfi/JPZPcv7V21Qp19+rXly9TfnmSPfp62wwAAACwGcxVAKiqe6e76/+K1to/D7HLft1WqthaO2SZNp2X5OAB2gIAAFM3N2MAFj368+kkz1/l1xbu/O+1TPme/fqKNTQNAADmxtwEgCQ3S7J/koOSXL3o5V8tyQv6Oq/rtx3ff76wX++/dGdVddt0j/98caXn/wEAYLOYp0eArknyhmXKDk43LuCD6S76Fx4POjPJ/ZI8ZNG2BQ9dVAcAAEZhbgJAP+D3yZPKquqF6QLAm1trr19U9MYkv5nk6VX1xoV3AVTVLXLjDEITXyIGAACb0dwEgB3RWru4qn4jyauTfLSq3p7k20mOSXL7DDeYGAAA5sKmDgBJ0lr7s6rakuQ5SX4+3biH/0zyvNbam2fZNgAAWG+bIgC01l6Y5IXbKP+HJP+wXu0BAICNap5mAWJHVXULAACjJwAAAMCICAAAADAiAgAAAIyIAAAAACOyKWYBYpFtDfadVNba9NoCAMCGowcAAABGRA/AZjPpjv7CnX93+wEARk8PAAAAjIgAsNl46RcAANsgAAAAwIgIAAAAMCIGAY+Bwb8AAPT0AAAAwIjoAZhnXvoFAMB20gMAAAAjogdgnnnpFwAA20kPwLwy3z8AADtAAAAAgBERAAAAYEQEAAAAGBGDgDcbg38BANgGAWAemO8fAICBeAQIAABGRA/APDDfPwAAA9EDAAAAIyIAAADAiAgAAAAwIgIAAACMiEHA88rgXwAAdoAeAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAAAAgBERAAAAYESqtTbrNsy1qvrabrvtts9BBx0066YAALCJXXDBBbnqqqu2ttb2Xct+BIA1qqqLk+yZZEuSA/vNn5pZg5gXzhVWw3nCajlXWC3nynzbL8kVrbU7rWUnAsCAquq8JGmtHTLrtrCxOVdYDecJq+VcYbWcKyTGAAAAwKgIAAAAMCICAAAAjIgAAAAAIyIAAADAiJgFCAAARkQPAAAAjIgAAAAAIyIAAADAiAgAAAAwIgIAAACMiAAAAAAjIgAAAMCICAArqKrbV9UJVfXlqrqmqrZU1fFVdYvt3M8+/fe29Pv5cr/f20+r7ayftZ4nVbVHVT2uqv6mqj5VVd+qqm9U1Uer6tlVddNp/wbWx1B/pyzZ52FVdX1Vtap68ZDtZTaGPE+q6h5VdWJVfaHf1yVV9f6q+vlptJ31NeB1yv2r6u/7719dVZ+vqtOr6iHTajuz40Vg21BVd0nyoSS3TvL3ST6V5N5JjkxyYZL7tda+tor97NvvZ/8kZyb51yQHJjk6ySVJ7ttau2gav4HpG+I86f+CfVeSrUnOSvJfSfZJ8rAkt+n3/6DW2tVT+hmsg6H+Tlmyz5tytKY4AAAOXklEQVQnOT/JLZPcLMkfttaeN2S7WV9DnidVdVyS1ye5Msk7k2xJsneSuyf5cmvt2IGbzzoa8Drll5P8ZZJvJTk1yReT3D7JI5PsnuR5rbU/nMZvYEZaa5ZlliTvTtKS/OqS7a/st792lfv5q77+K5dsf0a//YxZ/1bLbM+TJD+S5HFJbrpk+82TnNfv59mz/q2W2Z8rE/Z5Qrrg+Lv9Pl48699p2RjnSZL7JLkuyceT3GZC+ffM+rdaZn+uJPmeJJcluSrJAUvKDkpydboAucusf69luEUPwDKq6s5JPpvubsldWms3LCq7eZKvJKkkt26tfWsb+9kjyf8kuSHJbVtr31hUdpP+GPv1x9ALMGeGOk9WOMZjk/zfJO9srT1szY1mJqZxrlTV0UlOS/KEJDsneWP0AMy1Ic+TqjonyQOS3KO19h9TazQzMeB1yvcm+WqS81trPzyh/Pwk90hyy7adPZRsXMYALO+B/fo9i/+jSpL+Iv7cdN1i91lhP/dNsluScxdf/Pf7uSHJe/qPR665xczCUOfJtlzbr69bwz6YvUHPlaq6dZLXJTmttfbWIRvKTA1ynvTjyx6Q5KNJPllVR1bVc/oxRQ/qb0Ax34b6O+WSdDcq96+quy4uqKr9k9w1ycdd/G8u/gJY3gH9+tPLlH+mX++/TvthY1qPf79P6tdnrGEfzN7Q58pfp/s7/KlraRQbzlDnyY8uqn9mv7wsycuTvDfJx6vqB9bQTmZvkHOldY+CPC3d3yfnVdWbq+qPq+rEdI+gfjLJowdoLxvIzrNuwAa2V7++fJnyhe17r9N+2Jim+u+3qp6e5CHpnuE9YUf2wYYx2LlSVU9KN4nAY1pr/z1A29g4hjpPbt2vfzbJpekGc74vya2SvCDdY2P/WFX3aK19e8ebywwN9ndKa+3kqvpykrclWTw71H+ne7TQI8qbjB6AHVf9eq2DKIbaDxvTDv/7rapHJjk+3bOZj2qtXbvCV5hvqzpXqmq/dOfFya21v51ym9h4Vvt3yk6L1k9urZ3aWruitfbZJE9M92jQ/kkeNZ1msgGs+v8/VfX4dD1DH0g38Hf3fv2+JH+e5KQptZEZEQCWt5Cc91qmfM8l9aa9Hzamqfz7rapHpPsL95IkRxggvikMda6ckG62jl8ZolFsOEOdJ1/v19ckOX1xQf/Ix9/3H++9vQ1kwxjkXOmf8z8h3aM+T2itfaq1dlVr7VPpeorOS/Loqjpi7U1moxAAlndhv17u2bmFgTLLPXs39H7YmAb/91tVj05ycrqu18Nbaxeu8BXmw1DnysHpHu/4n/7FX62qWrpu+iR5br/ttLU1lxkZ+v8931g6QLS3EBB22462sbEMda78RLqpQN8/YTDxDUnO6T8esiONZGMyBmB5Z/Xrn6iqm0yYXut+6e7CfXiF/Xy4r3e/qrr5hGlAf2LJ8ZgvQ50nC995bJITk3wpyZHu/G8qQ50rJ6brnl/qrkkOSzde5Lwk/7bmFjMLQ50n56d79v+WVfW9E8aK3L1fb1l7k5mRoc6VXfr1rZYpX9hurMgmogdgGf1zku9JN0f/05YUvyjJHklOXDy3blUdWFUHLtnPN5O8pa//wiX7eXq//3e70JtPQ50n/fYnpjtXPp/kMOfE5jLg3ynPaK09eemSG3sA/rHf9hdT+zFMzYDnyXXpXkKZJC9dPO1nVd0jyXHpphY+ZeCfwDoZ8P8/H+jXx1TVDy0uqKofSXJMunEEZw7XembNi8C2YcIrti9I8mPp5uz/dJJDF8+L23fDp7VWS/azb7+f/dP9B/Qv6QbXHJ3uGe9D+/+QmUNDnCdVdWS6AVg3Sfcs5hcmHOqy1trxU/oZrIOh/k5ZZt/HxYvANoUB/9+ze7pBnPdJ1yN0drq7uY9K9+jPs1trr5zyz2GKBjxXTkjyC+nu8p+a5HPpgsUjktw0yfGttV+f8s9hHQkAK6iq70/y++mmYtw33Zv1Tkvyotba1iV1l/2fdVXtk27qtUckuW2SryV5V5Lfa619cZq/gelb63my6OJtWz7XWttvuFYzC0P9nTJhv8dFANg0Bvx/z+5JfjPJsUnulOTqJP+a5BWttXdN8zewPoY4V6qq0s0OdVySH05y8yRXpAuOr2utmQVokxEAAABgRIwBAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAAAAgBERAAAAYEQEAAAAGBEBAGCTqKoTq+qSqtpjjfvZUlVbBmrWzMzqd1TVs6vq2qo6cL2PDbAaAgDABlNVbclyfVVdWlVnVtXjlvnOvZI8PsmftNa+tb4tno2qOruq2qzbMcFfJrkkyctn3RCASXaedQMAWNaL+vX3JDkgySOSHFlVh7TWnrWk7h8luSLJa9axfRvdg2Zx0NbaVVX1p0leUlWHttY+NIt2ACynWtuIN08AxmvhrnZrrZZsf1CSf+o/3rm1tqXfvn+STyV5fWvtFwc4/pb++PutdV/TVFVnJzl86Z/TRlBV35fk80lOaq09ftbtAVjMI0AAc6K19r50F/qV5EcXFT2p3/b25b5bVT9ZVaf3jxJdU1WfraqXVdXe29OGqvq5qjqrqr5eVVdX1QVV9byq2mVRndv1jy19bBv7OaN/vOnui7YdV1XvqKqLquqqqrqiqs6tqscv+e5+fUg6vP+8+HGpsxfVmzgGoKp2qarfrqrzq+rK/jgfqKqfnVB3v36/b+r/+aT+z/DqqvpoVf30pN/XWvtykg8kOaaq9tzGHynAuvMIEMB8Wbjbvbj79sFJrk/y4YlfqPq9dI8TbU3yznTPp/9QkuckOaqq7ttau2LFA1e9IV3Y+GKSv0tyWZL7JPmDJA+qqh9vrV3XWvtSVb03yU9U1T1aa/++ZD+37dt8XmvtPxYVvSbJfyY5J8lXkuyb5Kgkb6mqA1prz+/rXdb/nuOS3DE3PiqVJFtW+A03TfLudOHhU0n+IsnuSY5J8vaq+pHW2u9O+Oodk/xLkouSvCXJPkkek+Tvq+rBrbWzJnzn3CRHJDks3Z87wMbQWrNYLBbLBlrSXdy3CdsfnOSGfrljv22PJNcl+fdl9nVkv78PJdl7SdlxfdmrlmzfkmTLMnX/LsluS8pe2Jf92qJtP9dve/mENv1GX/arS7bfZULdmyZ5X5Jrk9xuSdnZk/6cVvgdv9Mf+/QkOy/afuu+fkty6KLt+y38+0jygiX7+smFfS1z/KP78pfO+pyyWCyWxYtHgAA2qKp6Yb/8YVWdkuSMdD0Ax7fWPtdXu12SndLdMZ/kGf36Ka21yxYXtNbelOTjSSbOLLTEr6ULGk9qrV21pOwPknxtyX5OS3J5ksdV1U5L6j8x3QX925a057NLD9pa+3a6u/Q7Z5hBvU9Kd1H+rNbadYuOc0m635EkT57wvc8lefGStr073XP+917mWF/t13dYS4MBhuYRIICN6wX9uqV77OUDSd7QWnvrojr79uuvL7OP+6a72H50VT16QvlNk9yqqvZtrX1t0g6qavckP5zk0iTPrJo45vaaJActfGjdTDh/m+Qp6e6Un97v65Akd0tyamvt0iXHuUOS30p3oX+HJLstOcbtlvmNq1JVN0/yA0m+1Fr71IQqZ/bre04o+3hr7foJ27+Q7s94kq39+pbb1VCAKRMAADaotrrZbRbuxu+6TPm+6f6uf8Ey5Qtulu4u/iS3SNfzcKtV7GexN6ULAE9MHwD6f06SNy+uWFV3TveM/S3SBZ33pOtBuD7dYzhPTLJL1mavfr1cb8nC9kkDoy+bsC3pekWW601fCDBLe0wAZkoAAJhvl/TrfZcpvzzJTVpr+6zhGJf3639rrR282i+11j5UVZ9JcnQ/29C30o0NuDQ3BoIFz0r3G36hfzTpf1XVz+XG4LAWC7/jNsuU33ZJvbVa+HdyyTZrAawzYwAA5ttXkvxPuheFTfLhJLeoqrvt6AFaa99M8skkd6uq7Q0Sb0535/4xSX4q3eMwf9Nau3ZJvR/o1++YsI/Dl9n39UkyYYzBRK21byT5bJLbVdVdJ1Q5sl8vO33pdjqwX398oP0BDEIAAJhjrbWWbtrMW1bVD0yo8qp+/br+5VTfoar2qKr7rOJQr0w3XuCESe8OqKpbVNWk3oET081a9PP9knSPBi21pV8fsWS/P5nJg3KTGx9Z2p5Btieke5zpZYuDQ1XdMsnzF9UZwsKf66QpQgFmxiNAAPPvHUkelW6w7X8tLmitva+qfjvJHyf5TFWdnuTidM/83zHd3fUPJnnItg7QWjuhH8D7K0k+W1ULM+Dsk+RO6ea6f2OSpy753heq6qx0A3sXpiv9twmH+Mskv5Dk5Kp6R5IvJbl7366/TdeDsNT7kjw6yd/1v+uqJJ9rrb1lGz/l5Ukemm6Kzk/039u938+t003Z+cFt/VmsRlXdJN1vvrB957sOAGZOAACYf+9I8t/p7rD/xdLC1tpLqurcdFOC3j/dxe/l6S6y/zrJ36zmIK21p1XVu9Jd5D843WDZremCwMuSvHWZr74p3cXwzlky+HfRvs+vqiPTTbV5VF/3E0kemW4A7qQA8Pp0IebYJL/Zf+f96V7Utdxv+HZV/Xi6MQePTfKr6YLJJ5I8s7X2tuW+u50enG7Wol8faH8Ag6mu9xiAeVZVv5Pkj5IcvMwddtZR34txeLqXmw01qBhgEAIAwCZQVbsmuTDJ+a21h826PWNWVT+SbiDxM1prfz7r9gAsZRAwwCbQWrs6yROSfLSq9ph1e0butukGFL921g0BmEQPAAAAjIgeAAAAGBEBAAAARkQAAACAEREAAABgRAQAAAAYEQEAAABGRAAAAIAREQAAAGBEBAAAABgRAQAAAEZEAAAAgBERAAAAYEQEAAAAGBEBAAAARuT/A9nJOwxf0cLkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 263,
       "width": 384
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng2 = np.random.RandomState(15)\n",
    "el2 = []\n",
    "prob2 = []\n",
    "air2 = []\n",
    "for _i in range(20):\n",
    "    vals = quad_distribution(rng2, 30., 70., 85.)\n",
    "    _elevation = vals[1]\n",
    "    _airmass = 1./np.cos((90. - _elevation)*np.pi/180.)\n",
    "    prob2.append(vals[0])\n",
    "    air2.append(_airmass)\n",
    "    el2.append(_elevation)\n",
    "plt.plot(prob2, el2,'r+')\n",
    "plt.xlabel('P(elevation)')\n",
    "plt.ylabel('elevation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x8486df710>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAusAAAH0CAYAAACEkWPuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X2UbWddJ/jvTxPpgLwI3CI96BiTIeTlMomkW4kocJWhbR3bqCnGpdJBB7vFW8M0ktUtmm6D3S60X2iUlOMrpgHnhXIWOKOANHqV+DZ2X0TW5YZ3E0ENVQExgCHy8swf+xxSVOpU1a3adc4+pz6ftfba5+y9z3Ofyt2p+62nnue3q7UWAABgeD5v1h0AAAC2J6wDAMBACesAADBQwjoAAAyUsA4AAAMlrAMAwEAJ6wAAMFDCOgAADJSwDgAAAyWsAwDAQAnrAAAwUMI6AAAMlLAOAAADJawDAMBACesAADBQwjoAAAzUebPuwDRV1Z8meViSO2bcFQAAFttFSe5prX3ZQRo5UmE9ycMuuOCCR15++eWPnHVHAABYXLfffnvuvffeA7dz1ML6HZdffvkjT58+Pet+AACwwK655pq85S1vueOg7ZizDgAAAyWsAwDAQAnrAAAwUMI6AAAMlLAOAAADJawDAMBACesAADBQwjoAAAyUsA4AAAMlrAMAwEAJ64dgdTVZX9/9uvX17loAANhOL2G9qq6vqpdV1W1VdU9Vtap61T7b+saqemNVfaCq7q2q91XVWlVd20dfD9vqarKykpw4sXNgX1/vrllZEdgBANheXyPrNyVZSXJ1kj/fbyNV9RNJfi3JE5O8IclPJnlLkm9O8ntV9V0H7+rhWl5OrrgiOXt2cmAfB/WzZ7trl5en308AAIavr7D+/CSXJnlYkufup4GqujDJjUk+mOSK1tpzWms/2Fq7Psk/SFJJfrSn/h6apaXk1KnJgX1rUD91qvsMAABs1UtYb62daq29u7XWDtDMl4768/+11j5nPLq1dirJR5McO0D7UzMpsAvqAACci/Nm3YFN3p3kb5N8RVU9urV29/hEVT0lyUOTvHZWnTtX48A+DufHj3fHNzYEdQAA9mYwYb219uGq+hdJXpLkbFW9NsmHklyS5B8l+c9J/ule2qqq0xNOXdZHX/dqHNiPH+9CepIcOyaoAwCwN4MJ60nSWntpVd2R5OVJvnfTqfckuXXr9BgAAFhkg6qzXlX/PMmvJLk13Yj6Q5Jck+R9SX65qv7tXtpprV2z3ZbkHYfU9W2N56hvbHQj6seOda93K+sIAADJgMJ6VT0tyU8k+X9aaz/QWntfa+1vWmtvSfIt6UpCvqCqLp5lP/dq62LSM2e6bbeyjgAAMDaYsJ7kfxztT2090Vr7myR/lK6/Xz7NTu3HpKovu5V1BACAzYYU1h802k8qzzg+/rdT6Mu+7VaeUWAHAGCvph7Wq+r8qrqsqi7Zcuq20f6fVNVjt3zmHyZ5cpJPJPn9KXRz39bWdq+jvjWwr61Nv58AAAxfL9Vgquq6JNeN3l442l9bVbeOXt/dWrtx9PqxSW5PcmeSizY18ytJ3pTk6Ulur6rXJLkryeXppshUkh9srX2ojz4flpMnu/3y8s7lGceBfW3t/s8AAMBmfZVuvDrJDVuOXTzaki6Y35gdtNY+U1XfkORkkm9Pt6j0wUk+nOR1SX6qtfbGnvp7qPYavpeWBHUAACbrJay31m5OcvMer70j3Sj5duc+meSlow0AAI60IS0wXWirq3tbSLq+3l0LAADC+hSsriYrK7tXfhlXkllZEdgBABDWp2J5efdSjVtLPi4vT7+fAAAMi7A+BbvVVt+tNjsAAEeTsD4lkwK7oA4AwCR9lW5kD8aBfRzOjx/vjm9sCOoAADyQkfUpGwf2Y8e6kL6x0b0W1AEA2EpYBwCAgRLWp2w8R308oj4eYd+trCMAAEePsD5FWxeTnjnTbbuVdQQA4GgS1qdkUtWX3co6AgBwdAnrU7BbeUaBHQCA7QjrU7C2tnsd9a2BfW1t+v0EAGBY1FmfgpMnu/3y8s7lGceBfW3t/s8AAHB0GVmfkpMnHxjUV1cfON1laemBQX19vbsWAICjRVifkdXVZGVl9/np4/nuKysCOwDAUSOsz8jy8u4LSrcuTF1enn4/AQCYHWF9RnarALNbBRkAABafsD5DkwK7oA4AQKIazMyNA/s4nB8/3h3f2BDUAQCOOiPrAzAO7MeOdSF9Y6N7LagDABxtwjoAAAyUsD4A4znq4xH18Qj7bmUdAQBYbML6jG1dTHrmTLftVtYRAIDFJ6zP0KSqL7uVdQQA4GgQ1mdkt/KMAjsAAML6jKyt7V5HfWtgX1ubfj8BAJgdddZn5OTJbr+8vHN5xnFgX1u7/zMAABwNwvoM7TV8Ly0J6gAAR5FpMAAAMFDCOgAADJSwDgAAAyWsAwDAQAnrAAAwUMI6AAAMVC9hvaqur6qXVdVtVXVPVbWqetU5tvHs0ed22j7dR38BAGAe9FVn/aYkVyX5WJIPJLlsH228NcmLJpz7miRfm+T1++odAADMob7C+vPThfT3JHlqklPn2kBr7a3pAvsDVNUfjF7+3H47CAAA86aXsN5a+2w4r6o+mvysqjqe5ElJ/jzJr/faOAAADNg8LDD9p6P9L7bWzFkHAODI6GsazKGoqguSfFeSzyT5hXP43OkJp/Yzlx4AAGZi6CPrz0zyiCSvb629f9adAQCAaRr0yHqSfzLa/+y5fKi1ds12x0cj7k88aKf2ZTyXv7WZ/PEAAMyfwY6sV9UVSb4qXZWZ1824OwAAMHWDDeuxsBQAgCNukGG9qv5OkmelW1j6izPuDgAAzMTU56xX1flJLknyydbaeydctpzki5L8moWlAAAcVb2E9aq6Lsl1o7cXjvbXVtWto9d3t9ZuHL1+bJLbk9yZ5KIJTY4Xls7nE0t3ejDUducsOgUAYBt9jaxfneSGLccuHm1JF8xvzB5U1eVJvjoWlgIAcMT1EtZbazcnuXmP196RZOLQc2vt9p3Oz4XtRsqVbgQA4BwNcoEpAAAgrAMAwGAJ6wAAMFDCOgAADNTU66wfWRaWAgBwjoysAwDAQAnrAAAwUMI6AAAMlLAOAAADJawDAMBACesAADBQwvosVXUbAABsQ1gHAICBEtYBAGCghHUAABgoYR0AAAbqvFl34MjYaSHpdudaO7y+AAAwF4ysAwDAQBlZn5btRsrHI+pG0QEA2IaRdQAAGChhHQAABkpYBwCAgRLWAQBgoCwwnSULSwEA2IGRdQAAGChhHQAABkpYBwCAgRLWAQBgoIT1oam6/8mmAAAcacI6AAAMlLAOAAADJawDAMBACesAADBQnmA6SzstJN3unCeeAgAcKUbWAQBgoHoJ61V1fVW9rKpuq6p7qqpV1asO0N7XVNX/XVV/WVX3jfZvrKpv6KO/g9HaA7e9nAMA4EjoaxrMTUmuSvKxJB9Ictl+G6qqm5L86yR3J/m1JH+Z5NFJvjzJ05K87oB9BQCAudBXWH9+upD+niRPTXJqP41U1XK6oP6mJN/aWvvolvPnH7CfAAAwN3oJ6621z4bz2ufTN6vq85L8RJK/SfIdW4P66M/55H77CAAA82ZI1WC+KsmXJfmVJH9VVd+Y5HiSTyT5o9baH8yycwAAMG1DCut/f7T/YJK3JHnC5pNV9eYk17fWNnZrqKpOTzi177n0U2MhKQAAI0Mq3bg02n9fkguSPD3JQ9ONrv9GkqckWZtN1wAAYPqGNLL++aN9pRtB/5PR+7dX1bckeVeSp1bVtbtNiWmtXbPd8dGI+xP76jAAABymIY2s/9Vo/75NQT1J0lq7N93oepJ8xVR7BQAAMzKksP7O0f4jE86Pw/wFU+gLAADM3JDC+puTfCrJ46rqC7Y5f3y0v2NqPQIAgBmaelivqvOr6rKqumTz8dba3Un+ryQPT/Kvtnzmf0jyD5L8dZI3TKuvAAAwS70sMK2q65JcN3p74Wh/bVXdOnp9d2vtxtHrxya5PcmdSS7a0tQPJPnKJD9cVU9J8kdJvjTJtyT5dJLvba1NmiYDAAALpa9qMFcnuWHLsYtHW9IF8xuzi9baelV9ZZKb0gX0JyX5aJJfT/Li1tof9tRfAAAYvF7Cemvt5iQ37/HaO9KVZ5x0/sPpRth/oIeuAQDA3BrSAlMAAGATYR0AAAZKWAcAgIES1plodTVZX9/9uvX17loAAPolrLOt1dVkZSU5cWLnwL6+3l2zsiKwAwD0TVhnW8vLyRVXJGfPTg7s46B+9mx37fLy9PsJALDIhHW2tbSUnDo1ObBvDeqnTnWfAQCgP8I6E00K7II6AMB09PUEUxbUOLCPw/nx493xjQ1BHQDgsBlZZ1fjwH7sWBfSNza614I6AMDhEtYBAGCghHV2NZ6jPh5RH4+w71bWEQCAgxHW2dHWxaRnznTbbmUdAQA4OGGdiSZVfdmtrCMAAP0Q1tnWbuUZBXYAgMMnrLOttbXd66hvDexra9PvJwDAIlNnnW2dPNntl5d3Ls84Duxra/d/BgCAfgjrTLTX8L20JKgDABwG02AAAGCghHUAABgoYR0AAAZKWAcAgIES1gEAYKCEdQAAGChhHQAABkpYBwCAgRLWAQBgoIR1AAAYKGEdAAAGSlgHAICBEtYBAGCghHUAABgoYR0AAAZKWAcAgIES1gEAYKB6CetVdX1Vvayqbquqe6qqVdWr9tHOHaPPbrfd1UdfAQBgXpzXUzs3JbkqyceSfCDJZQdo66+TvHSb4x87QJsAADB3+grrz08X0t+T5KlJTh2grY+01m7uo1MAADDPegnrrbXPhvOq6qNJAAA48voaWe/Tg6rqu5L8t0k+nuRtSd7cWvv0bLsFAADTNcSwfmGSV2459qdV9d2ttd/ZSwNVdXrCqYPMpQcAgKkaWunGX0rydekC+0OSPCHJzya5KMnrq+qq2XUNAACma1Aj6621F205dCbJ91XVx5K8IMnNSb5lD+1cs93x0Yj7Ew/YTQAAmIqhjaxP8jOj/VNm2gsAAJiieQnr66P9Q2baCwAAmKJ5CevXjvbvm2kvAABgiqYe1qvq/Kq6rKou2XL8yqp65DbXf2mSW0ZvXzWNPgIAwBD0ssC0qq5Lct3o7YWj/bVVdevo9d2ttRtHrx+b5PYkd6ar8jK2nOQHq+pUkj9N8tEklyT5xiR/J8nrkvz7PvoLAADzoK9qMFcnuWHLsYtHW9IF8xuzs1NJHp/ky9NNe3lIko8k+d10dddf2VprPfUXAAAGr5ew3lq7OV1Zxb1ce0eS2ub47yTZ00OPAADgKJiXBaYAAHDkCOsAADBQwjoAAAyUsA4AAAMlrAMAwEAJ6wAAMFDCOgAADJSwDgAAAyWsAwDAQAnrAAAwUMI6AAAMlLAOAAADJawzd1ZXk/X13a9bX++uBQCYV8I6c2V1NVlZSU6c2Dmwr69316ysCOwAwPwS1pkry8vJFVckZ89ODuzjoH72bHft8vL0+wkA0AdhnbmytJScOjU5sG8N6qdOdZ8BAJhHwjpzZ1JgF9QBgEVz3qw7APsxDuzjcH78eHd8Y0NQBwAWh5F15tY4sB871oX0jY3utaAOACwKYZ1zU9VtAAAcOmGduTWeoz4eUR+PsO9W1hEAYF4I68ylrYtJz5zptt3KOgIAzBNhnbkzqerLbmUdAQDmjbDOXNmtPOPSUvKsZyWPf/zugX193dNNAYBhE9aZbLyYdPO2l3OHaG1t5zrqq6vJC1+YtHZ/YF9be2A749C/siKwAwDDJawzV06eTG65ZXJ5xuXlLsi/611dYH/xi7vPbLZ1dH55eTp9BwA4V8I6k7X2wG0v5w7ZyZOT66hvnrf+rnclr3zl506D8ZRTAGCeCOssnEkLTQV1AGDenDfrDsBhGAf2cTg/frw7vrEhqAMA88PIOgtrHNjHD0saPzxJUAcA5oWwDgAAAyWsc26mvJj0IMZz1Mcj6uMRdg9LAgDmhbDOQtq6mPTMmW7zdFMAYJ4I6yycSVVfJlWJAQAYKmGdhbJbeUaBHQCYJ72E9aq6vqpeVlW3VdU9VdWq6lU9tPusUVutqp7TR19ZbGtru9dR3xrY19am308AgL3oq876TUmuSvKxJB9IctlBG6yqL0nyslGbX3jQ9jgaTp7s9svLO5dnHAf2tbX7PwMAMDR9TYN5fpJLkzwsyXMP2lhVVZJfSvKhJD9z0PY4Wk6e3Fsd9aUlQR0AGLZeRtZba6fGr7ucfWDPS/K1SZ422gMAwJEzuAWmVXV5kh9P8pOttTfPuj8AADArgwrrVXVeklcm+bMkP3SAdk5vt6WHufQsttXVvVWHWV/vrgUAOEyDCutJ/lWSL0/y7NbavbPuDEfL6mqysrJ7OcdxeciVFYEdADhcgwnrVfUV6UbT/0Nr7Q8O0lZr7ZrttiTv6KWzLKTl5d3rr2+t4768PP1+AgBHxyDC+qbpL+9K8i9n3B2OqN0emLTbA5cAAPo2iLCero76pUkuT/KJTQ9Cakl+ZHTNz4+OvXRmvWThTQrsgjoAMAt9PRTpoO5L8osTzj0x3Tz2303yziQHmiIDuxkH9nE4P368O76xIagDANM19bBeVecnuSTJJ1tr702S0WLS50y4/uZ0Yf0/tdZ+YVr95GgbB/bjx7uQniTHjgnqAMB09RLWq+q6JNeN3l442l9bVbeOXt/dWrtx9PqxSW5PcmeSi/r48wEAYBH1NbJ+dZIbthy7eLQlXTC/MTAnxnPUNza6EfWke/2UpyRvfvPOo+vr68naWnLy5HT6CgAsrl4WmLbWbm6t1Q7bRZuuvWPrsT22bQoMU7F1MemZM8kP/EB37p3v7AL7pDrsarADAH0aSjUYGIRJVV++53uSSy/trpkU2NVgBwD6JqzDyE7lGZeWkttumxzYlXYEAA6DsA4ja2s7h+3tAvvLXy6oAwCHR1iHkZMnk1tu2Tlsbw3sL3lJV95RUAcADsNQHooEg7CXCi7jwK4GOwBw2IysAwDAQAnrcI621mA/dqx7feLE5JKOAAD7IazDOdiuBvuZM93rs2cFdgCgX8I67NGkqi9LS91rgR0A6JuwDnuwW3lGgR0AOAzCOuzBbjXYkwcG9rW16fcTAFgsSjfCHoxLOi4v71yecRzY19b2VgYSAGAnwjrs0V7D99KSoA4A9MM0GAAAGChhHQAABkpYBwCAgRLWAQBgoIR1AAAYKGEdAAAGSlgHAICBEtYBAGCghHUAABgoYZ3FUNVtAAALRFgHAICBEtZhoFZXk/X13a9bX++uBQAWj7AOA7S6mqysJCdO7BzY19e7a1ZWBHYAWETCOgzQ8nJyxRXJ2bOTA/s4qJ892127vDz9fgIAh0tYZ/6MF5Nu3vZybo4sLSWnTk0O7FuD+qlT3WcAgMUirMNATQrsgjoAHB3nzboDcM5ae+Cx8Qj6dufm2Diwj8P58ePd8Y0NQR0AjgIj6xwNCzAl5tixLqRvbHSvBXUAWHzCOgAADJSwDgM3nqM+HlEfj7DvVtYRAJh/wjoM2NbFpGfOdNt40elVVyVvf/ve2lGHHQDmjwWmLIYFW1ia7Fz15dSpLqjfdVdyzTXJ6dPJlVfu3k6SnDw5nf4DAAfXy8h6VV1fVS+rqtuq6p6qalX1qn208xNV9ZtV9f6qureqPlxVf1xVP1JVj+qjrxwBC1CHfbfyjEtLyZvelDzoQcl993WBfbsRdg9OAoD51tc0mJuSrCS5OsmfH6Cd5yd5SJL/nOQnk/xykk8luTnJ26rqSw7WTZgPa2u711G/8spuRH0c2J/+dA9OAoBF09c0mOcn+UCS9yR5apJT+2znYa21T2w9WFU/luSHkrwwyffvt5McEQtQh308VWV5eeeAPQ7sT396NyXmxIkulCeCOgAsgl7Cemvts+G8DjCtYLugPvLqdGH9cftuHObMXueWX3ll8id/4sFJALCI5qUazDeN9m+baS9goCY9OOlZz9rb51WLAYBhGmQ1mKq6MckXJnl4kr+X5KvTBfUf3+PnT084dVkvHYQ5cO+9yQtfmLzylTuPrqsWAwDDNciwnuTGJI/Z9P4NSZ7dWtuYUX9g0LY+OCnpXj/oQV0IH89l3xrYVYsBgGEb5DSY1tqFrbVKcmGSb01ycZI/rqon7vHz12y3JXnHIXabIWttbhaXnqudHpx0332fG9hViwGA+TLIsD7WWvtga+01SZ6R5FFJXjHjLsGgTArc4znskwK7oA4A82HQYX2stXZnkrNJrqyqR8+6PxwxA3140l4enLRdYD9+vNsEdQAYvrkI6yP/zWj/6Zn2AgZiLw9O2hrYv/ALP7dajKAOAMM29bBeVedX1WVVdcmW45dV1YXbXP95o4ciLSX5/dbaX02rrzBkJ08mt9yye+AeB/YXvzi54ILp9Q8AOLheqsFU1XVJrhu9HQfua6vq1tHru1trN45ePzbJ7UnuTHLRpma+Psm/q6o3J3lvkg+lqwjz1HQLTO9K8r199BcWxbmUWXzlKx9YLWZSlRgAYBj6Kt14dZIbthy7eLQlXTC/MTt7U5KfS/LkJFcleUSSjyd5V5JXJvmp1tqHe+ovHBnbzW1P7j8msAPAcPUS1ltrNye5eY/X3pHkAav1WmtnkngcC7O100LS7c4NvBzkTotQT50S2AFg6OZpgSkMz0ArxSTnVi1muzrsAMDsCeuw2fjhSZu3vZwboHOtFnP2bPcZAGA4+pqzDgzMePHp8vLeqsWsrZ3bglUA4PAJ67DA9hq+l5YEdQAYItNggD1bXd3bvPb19e5aAOBghHXYq/Fi0s3bXs4tiNXVZGUlueqq5O1vn3zdeGHrykryzGdOr38AsIiEddjNHCwmnYbl5eTCC5O77kquuWb7wL65As2DHtTNgzfCDgD7J6zDXi1QpZj9WFpK3vSmLoTfd98DA/vWoH7ffV2VmeXl2fUZAOadsA7s2ZVXJqdPPzCwTwrqHrQEAAejGgxwTsaB/ZprulB+1VXJIx6RfOhDgjoA9M3IOnDOxoH98z8/+fSnu6BeJagDQN+EdWBfjh3rRtTHWkse9ShBHQD6JKzDQSzoYtLdjOeoj0fUxz7ykWRjY3b9AoBFI6wD52TrYtLxiPp4Ssykso4AwLkT1oE9m1T15ezZ5E/+ZHJZRwBgf4R1YE92K884qawjALB/wjqwJ2tru9dR3xrYn/70LuQDAPsjrAN7cvJk9zTS3cozjgP7hRcmd93VhXwAYH88FAnYs1e/Olld7UL7TuUZr7yym8O+ttaFfABgf4R14JzsNXwvLQnqAHBQpsEAAMBACeswT6o+9ylEAMBCE9YBAGCghHUAABgoYR0AAAZKWAcAgIFSuhGGaqeFpNuda+3w+gIAzISRdVg0KsYAwMIwsg5Dtd1I+TiEG0UHgCPByDqwMFZXk/X13a9bX++uBYChE9aBhbC6mqysJCdO7BzY19e7a1ZWBHYAhk9YBxbC8nJyxRXJ2bOTA/s4qJ892127vDz9fgLAuRDWYZ6NF5Nu3nY7t6ALUJeWklOnJgf2rUH91KnuMwAwZMI6zJPWLC7dwaTALqgDMK9Ug4F5pmLMA4wD+zicHz/eHd/YENQBmD+9jKxX1fVV9bKquq2q7qmqVlWvOsc2HlVVz6mq11TVe6rq3qr666r63ar6n6vKbwGAPRkH9mPHupC+sdG9FtQBmDd9jazflOSqJB9L8oEkl+2jjeUk/1uSv0xyKsmfJXlMkm9N8gtJ/mFVLbd2RIcLAQA4cvoarX5+kkuTPCzJc/fZxruS/KMkX9xa+87W2gtba9+TLvi/P8m3pQvuwF7tZwHqAhjPUR+PqI9H2Hcr6wgAQ9NLWG+tnWqtvfsgo96ttd9qrf2/rbXPbDl+V5KfGb192gG6CRwBWxeTnjnTbbuVdQSAIZqXeeCfHO0/NdNewDzYXDFm/HrztvW67c7NqUlVX3Yr6wgAQzX4sF5V5yX5x6O3b5hlX4Dh2q08o8AOwDyah9KNP57keJLXtdZ+Yy8fqKrTE07tZ+ErMAfW1navo761rOPaWnLy5PT7CgB7NeiwXlXPS/KCJO9I8qwZdwcYsHHoXl7euTzjOLAL6gDMg8GG9ao6meQnk5xN8nWttQ/v9bOttWsmtHk6yRP76SGQ5HOryMx43vtew/fSkqAOwHwY5Jz1qvpnSW5JcibJiVFFGKAPC7KYFACOgsGF9ar6F0n+Y5K3pgvqloABAHAkTT2sV9X5VXVZVV2yzbl/mW5B6el0U1/unnb/AABgKHqZs15V1yW5bvT2wtH+2qq6dfT67tbajaPXj01ye5I7k1y0qY0bkvxokk8nuS3J8+qBT1S8o7V269aDAACwiPpaYHp1khu2HLt4tCVdML8xO/uy0f7zk/yzCdf8TpJb99E/oA8P/AF69/PmxwPAvvUyDaa1dnNrrXbYLtp07R1bj+2xjWqtPa2P/gIAwDwYbOlGYIC2GyUfUOlGAFg0g6sGAwAAdIR1AAAYKGEdAAAGSlgHAICBssAUOBiLSgHg0BhZBwCAgRLWAQBgoIR1AAAYKGEdAAAGSlgHAICBEtaBxVbVbQAwh4R1AAAYKGEdAAAGSlgHAICBEtYBAGCgzpt1BwB6s9NC0u3OtXZ4fQGAHhhZBwCAgTKyDiyO7UbKxyPqRtEBmENG1gEAYKCEdQAAGChhHQAABkpYBwCAgbLAFFhsFpYCMMeMrAMAwEAJ6wAAMFDCOgAADJSwDgAAAyWsAwDAQAnrAAAwUMI6AIO1upqsr+9+3fp6dy3AohHWARik1dVkZSU5cWLnwL6+3l2zsiKwA4tHWAdgkJaXkyuuSM6enRzYx0H97Nnu2uXl6fcT4DAJ6wAM0tJScurU5MC+NaifOtV9BmCRCOsADNakwC6oA0dFL2G9qq6vqpdV1W1VdU9Vtap61azaAWBxbA3sx493m6AOHAXn9dTOTUmuSvKxJB9IctmM2wFggYwD+/HjycZGd+zYMUEdWHx9TYN5fpJLkzwsyXMH0A4AAMy9XsJ6a+1Ua+3drbU2hHYAWCzjOeobG92I+rFj3evdyjoCzDsLTAEYtK2LSc+c6bbdyjoCLAJhHYDBmlT1ZbfEOfSFAAAO/klEQVSyjgCLoq8FpoNSVacnnLJgFWBO7FaecRzYx9ecOGHBKbB4jKwDMEhra7uXZ9w6wr62Nv1+AhymhRxZb61ds93x0Yj7E6fcHQD24eTJbr+8vPNo+Tiwr63d/xmARbGQYR2AxbDX8L20JKgDi8k0GAAAGKipj6xX1flJLknyydbae6f95wMw56q6vUdyAEdAL2G9qq5Lct3o7YWj/bVVdevo9d2ttRtHrx+b5PYkdya56ADtAADAQutrZP3qJDdsOXbxaEu6YL6XkN1XOwAAMPd6mbPeWru5tVY7bBdtuvaOrcf20w4AACw6C0wBAGCglG4EYLjGi0n3es6iU2DBGFkHAICBMrIOwHBtN1KudCNwhBhZBwCAgRLWAQBgoIR1AAAYKGEdAAAGygJTAOaLhaXAEWJkHQAABkpYBwCAgRLWAYBBW11N1td3v259vbsWFomwDgAM1upqsrKSnDixc2BfX++uWVkR2FkswjoAMFjLy8kVVyRnz04O7OOgfvZsd+3y8vT7CYdFWAcABmtpKTl1anJg3xrUT53qPgOLQlgHAAZtUmAX1DkK1FkHAAZvHNjH4fz48e74xoagzmIzsg4AzIVxYD92LFnfqKxvVI4dE9RZbMI6AAAMlLAOAMyF8Rz1jY37j21s7F7WEeaZsA4ADN7WxaRju5V1hHknrAMAg7Zd1Zexnco6zjtPbiVRDQYAGLD19WTpMZW3jw+cTfKY+89vd279g23uF5yOn9z60z+98wLazT/IJMnJk9PrI9NhZB0AGKy1tel8Zmg8uZUxYR0AGKyTJ5PVW1rWP9iStmkb23Rs/YMtq7e0hRhdnvjk1qqk6tAeCGXqzfAI6wDAoJ08ubcgurS0WNNAtgvsY4cV1FdWdp/7P/5BYWVFYJ8GYR0AYKC2BvaxvoN6cu5Tby68MHna03Zv1yj8wQjrAAADtvnJrWOH8eTWc5l6c+GFyV13Jc98plH4wyasAwCQZO9Tb970pgmj8KNgn1gA2xdhHQCYP1sXmi6icfCtytJjKusb9dlT6xvdsc3X9GUvU2+uvHLnGveHtQB2J4u6OFZYBwDgc+xl6s20F8DuZJEXxwrrAAADtP7BliuvaKl0+/UPtj2dm6ZpLoDdycTFsQswLUdYBwAYmN2mkUxcDHpQ+5h6M60FsDvZ7b/HLKbl9EVYBwAYmLW13YPl1oA6kye3zmBO/SRDmpbTJ2EdAGBgTp5Mbrll92A5Dqi33NLPA6HOZerNEA1lWk6fegnrVXV9Vb2sqm6rqnuqqlXVq/bZ1hdX1cur6i+q6r6quqOqXlpVX9RHXwEA5sG0n9x6rlNvrryi5e1ntg/wDwj8U6zcM4RpOX3qa2T9piQrSa5O8uf7baSqLklyOsl3J/mjJP8xyfuS/K9J/qCqHnXwrgIAzKlDLFm5n6k3T3/6535m7FDm0+9mRqUuD1tfYf35SS5N8rAkzz1AOz+dZCnJ81pr17XWfrC19rXpQvvjk/zYgXsKAMADnMvUm1e/+v6nmE51AewR1EtYb62daq29u7X9/6hXVRcneUaSO5JsrXz5I0k+nuRZVfWQfXcUAICJ9jr15rd/e3JQH5v6AtjWJs6rn+W0nIMa0gLTrx3t39ha+8zmE621jyb5vSQPTvKkaXcMAODI2mbqzawWwO5kuzn3Y/M8yj+ksP740f5dE86/e7S/dAp9AQBgB9uOwm8T7PtaALuTmdWln4LzZt2BTR4+2v/1hPPj44/YraGqOj3h1GXn2ikAAIbtXBbHjkP92trh/xDRhyGF9d2Ml+3OzyQjAAAO3Th0Ly/vbVrOvAT1ZFhhfTxy/vAJ5x+25bqJWmvXbHd8NOL+xHPvGgAAQ7Zt+N5mIek0puX0aUhz1t852k+ak/640X7SnHYAAFgoQwrr4zW7z6iqz+lXVT00yZOT3JvkD6fdMQAAmIWph/WqOr+qLhs9rfSzWmvvTfLGJBcl2frLiRcleUiSV7TWPj6VjgIAwIz1Mme9qq5Lct3o7YWj/bVVdevo9d2ttRtHrx+b5PYkd6YL5pt9f5LfT/JTVfV1o+u+MsmJdNNffriP/gIAwDzoa4Hp1Ulu2HLs4tGWdMH8xuyitfbeqvp7SX40ydcn+YYkf5nkp5K8qLX24Z76CwAAg9dLWG+t3Zzk5j1ee0fuL8O43fn3J/nuPvoFAADzbEgLTAEAgE2EdQAAGChhHQAABkpYBwCAgRLWAQBgoIR1AAAYqGqtzboPU1NVH7rgggseefnll8+6KwAALLDbb789995774dba486SDtHLaz/aZKHJbkjyWWjw++YWYeYF+4V9sq9wl64T9gr98p8uyjJPa21LztII0cqrG9WVaeTpLV2zaz7wrC5V9gr9wp74T5hr9wrJOasAwDAYAnrAAAwUMI6AAAMlLAOAAADJawDAMBAHdlqMAAAMHRG1gEAYKCEdQAAGChhHQAABkpYBwCAgRLWAQBgoIR1AAAYKGEdAAAGaqHCelV9cVW9vKr+oqruq6o7quqlVfVF59jOI0efu2PUzl+M2v3iw+o703PQ+6SqHlJV31lV/3tVvaOqPl5VH62q/1pVL6iqLzjsr4Hp6Ot7ypY2n1JVn66qVlX/ps/+Mjt93itV9YSqekVVvX/U1npV/U5V/ePD6DvT02NO+eqq+tXR5z9RVX9WVa+rqq8/rL4zOwvzUKSquiTJ7ydZSvKrSd6R5CuSnEjyziRPbq19aA/tPGrUzqVJfivJf0lyWZJvTrKe5NrW2vsO42vg8PVxn4y+Gb4+yYeTnEryniSPTPJNSS4ctf91rbVPHNKXwRT09T1lS5sPTfK2JI9O8oVJfqy1dlOf/Wb6+rxXqurZSX4hyd8k+bUkdyR5RJLjSf6itfbtPXefKekxpzw3yU8n+XiS1yT5QJIvTvKtSR6c5KbW2o8dxtfAjLTWFmJL8htJWpL/Zcvxl4yO/8we2/nZ0fUv2XL8eaPjb5j112qb7X2S5Ook35nkC7Ycf2iS06N2XjDrr9U2+3tlmzZfnu6HvB8atfFvZv112oZzryR5UpJPJXlrkgu3OX/+rL9W22zvkyTnJ/lIknuTPH7LucuTfCLdD3oPmvXXa+tvW4iR9aq6OMl7041AXNJa+8ymcw9N8pdJKslSa+3jO7TzkCQbST6T5O+21j666dznjf6Mi0Z/htH1OdPXfbLLn/EdSX45ya+11r7pwJ1mJg7jXqmqb07y2iTPSnJekl+KkfW51+e9UlVvTvI1SZ7QWjtzaJ1m6nrMKY9JcleSt7XWrtrm/NuSPCHJo9s5/uaP4VqUOetfO9q/cfP/AEkyCty/l+5XQ0/apZ1rk1yQ5Pc2B/VRO59J8sbR2xMH7jGz0Nd9spNPjvafOkAbzF6v90pVLSX5+SSvba29qs+OMnO93CujNVFfk+S/Jnl7VZ2oqhtH62C+bjRgxPzq63vKerpBxUur6nGbT1TVpUkel+StgvpiWZT/+R8/2r9rwvl3j/aXTqkdhmkaf7/fM9q/4QBtMHt93ys/l+777fcdpFMMUl/3yt/fdP1vjbZ/l+TfJ3lTkrdW1X93gH4yW73cJ62bDnEy3feT01X1n6rqxVX1inTTMN+eZLmH/jIg5826Az15+Gj/1xPOj48/YkrtMEyH+vdbVStJvj7dfNOX76cNBqO3e6WqvifdAvX/qbX2wR76xrD0da8sjfbPTHJ3usWCv5nkWJIfSTd96ter6gmttb/df3eZkd6+p7TW1qrqL5L8H0k2Vwj6YLrpdabpLphFGVnfTY32B52g31c7DNO+/36r6luTvDTdXMJva619cpePMN/2dK9U1UXp7ou11tqrD7lPDNNev698/qb9c1prr2mt3dNae2+SG9JNj7k0ybcdTjeZsT3/+1NV35Xuty23pVtU+uDR/jeT3JLk/zykPjIjixLWxz+RPnzC+Ydtue6w22GYDuXvt6quS/fNcT3J0yw+Xgh93SsvT1e14fv76BSD1Ne98lej/X1JXrf5xGjqw6+O3n7FuXaQQejlPhnNS395uukuz2qtvaO1dm9r7R3pfvtyOslyVT3t4F1mKBYlrL9ztJ8012u8CGPSXLG+22GYev/7rarlJGvpfv341NbaO3f5CPOhr3vliemmN2yMHoLUqqql+1V1kvzw6NhrD9ZdZqjvf38+unUB4sg4zF9wDn1jOPq6T56Rrnzj72yzUPUzSd48envNfjrJMC3KnPVTo/0zqurztimJ9OR0o1t/uEs7fzi67slV9dBtSjc+Y8ufx3zp6z4Zf+Y7krwiyZ8nOWFEfaH0da+8It2vqLd6XJKnpFvfcDrJHx+4x8xKX/fK29LNVX90VT1mm/UNx0f7Ow7eZWagr/vkQaP9sQnnx8eta1ggCzGyPprT98Z0NdBPbjn9oiQPSfKKzbVLq+qyqrpsSzsfS/LK0fU3b2lnZdT+bwhl86mv+2R0/IZ098qfJXmKe2Kx9Pg95Xmtteds3XL/yPqvj46tHtoXw6Hq8V75VLqH8iXJv91cqrGqnpDk2elKwv5Kz18CU9Djvz+3jfbXV9V/v/lEVV2d5Pp0895/q7/eM2sL8VCkZNvH+N6e5CvT1UR/V5Kv2lx3dPSr6LTWaks7jxq1c2m6m/2P0i3c+OZ0c5K/avQ/HXOoj/ukqk6kW9zzeenmDr5/mz/qI621lx7Sl8EU9PU9ZULbz46HIi2MHv/9eXC6RYJPSvfblt9ON1L6bemmv7ygtfaSQ/5yOCQ93icvT/Ld6UbPX5PkznQ/BFyX5AuSvLS19vxD/nKYooUJ60lSVV+S5EfTlc97VLongr02yYtaax/ecu3Ef1ir6pHpSmVdl+TvJvlQktcn+VettQ8c5tfA4TvofbIpaO3kztbaRf31mlno63vKNu0+O8L6Qunx358HJ/nnSb49yZele3z8f0nyH1prrz/Mr4HD18d9UlWVrkLQs5NcleShSe5J9wPez7fWVINZMAsV1gEAYJEsxJx1AABYRMI6AAAMlLAOAAADJawDAMBACesAADBQwjoAAAyUsA4AAAMlrAMAwEAJ6wAAMFDCOgAADJSwDgAAAyWsAwDAQAnrAAAwUMI6AAAMlLAOAAADJawDAMBACesAADBQ/z/a2UfvlvkLdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 250,
       "width": 373
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(prob1, air1,'bx')\n",
    "plt.plot(prob2, air2,'r+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(15)\n",
    "rand_test = []\n",
    "for _i in range(100):\n",
    "    rand_test.append(np.abs(0.5*rng.randn()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x848a77208>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAusAAAH0CAYAAACEkWPuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu0tFddJ/jvj4t2CEm4aGSWcRkSc1sLR0zoQEgLufTESLeaQXBm9RAgimsQ7AANS12CDXEapZ0WQri03dIhBpxmGl3KtNKSaQjXSNO+yHh7cyH4culwkcSGEEIEs+ePeg4cTk6d63OqdlV9PmvVet5T+6l99nueqqe+tWvv/VRrLQAAQH/uN+8GAAAAmxPWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6NQD5t2AWaqqv0pybJIjc24KAADL7cQkX2ytPXI/lYwS1qvqKUmemOTRSb4vyTFJfqu19rR91ntpkmuHH3+qtfaGfTU0Ofaoo4562BlnnPGwfdYDAABTHT58OHffffe+6xmrZ/0lmYT0LyX5VJLT91thVX1XktcMdT54v/UNjpxxxhkPO3To0EjVAQDAfZ111ln58Ic/fGS/9Yw1Zv0FSU7NZIjJT++3sqqqJG9McnuSX99vfQAAsIhG6VlvrV2/9u9Jzt63y5NckOS8YQsAACunu9VgquqMJK9I8urW2nvn3R4AAJiXrsJ6VT0gyZuSfCLJL8y5OQAAMFe9Ld34z5N8f5J/0Frb8/TZqpo2g3TfE18BAGBWuulZr6qzM+lN/7XW2h/Nuz0AADBvXfSsrxv+cnOSX9xvfa21s6b8nkNJztxv/QAAMAu99Kw/OJOlH89I8pWqamu3JC8d9vmN4b4r59ZKAACYoS561pPck+TfTSk7M5Nx7O9PclMSQ2QAAFgJMw/rVfXAJCcn+Wpr7dYkGSaTPmvK/i/LJKz/ZmvtDbNqJwAAzNsoYb2qLklyyfDjI4btOVV1zfDvz7fWXjT8+zuTHE7y8SQnjvH7AQBgGY3Vs/7oJM/YcN9Jwy2ZBPMXBQBYbGtXKm9tvu2AFTHKBNPW2staa7XF7cR1+x7ZeN8O6zYEBgCAldLLajAAAMAGwjoAAHRKWAcAgE71ss46ANCbtcmkOy0z6RRGp2cdAAA6pWcdANjcZj3llm6EmdKzDgAAnRLWAQCgU8I6AAB0SlgHAIBOmWAKAOyciaUwU3rWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAsMqqJjegS8I6AAB0SlgHAIBOCesAANApYR0AADr1gHk3AACYka0mkm5W1trBtQXYkVF61qvqKVX1mqp6X1V9sapaVb15l3U8vKqeVVW/W1Ufraq7q+oLVfX+qvrJqvItAAAAK2WsnvWXJPm+JF9K8qkkp++hjqcm+ddJPp3k+iSfSPIdSZ6c5A1Jfqiqntqaj/kAsCebvYWu9ah7e4UujRXWX5BJSP9okidmErZ36+YkP5LkD1pr967dWVW/kORDSX4sk+D+O/tuLQAALIBRhpa01q5vrd2yn17v1tq7Wmv/cX1QH+7/TJJfH348bx/NBACAhbIo48C/Omy/NtdWAADADHW/GkxVPSDJ04cf/3CHjzk0pWgvY+kBAGAuug/rSV6R5FFJ3t5ae8e8GwMAS8XEUuha12G9qi5P8sIkNya5dKePa62dNaW+Q0nOHKd1AABwsLods15Vz03y6iR/meT81todc24SAADMVJdhvaqen+S1Sf48k6D+mTk3CQAAZq67sF5VP5fkVUk+kklQ/9ycmwQAAHMx87BeVQ+sqtOr6uRNyn4xkwmlh5Jc2Fr7/KzbBwAAvRhlgmlVXZLkkuHHRwzbc6rqmuHfn2+tvWj493cmOZzk40lOXFfHM5L8UpK/S/K+JJfX2iWQv+FIa+2ajXcCAMAyGms1mEcnecaG+04abskkmL8oW3vksL1/kudP2ec9Sa7ZQ/sAAGDhjDIMprX2stZabXE7cd2+Rzbet8M6qrV23hjtBQCARdDdBFMAAGBCWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsA2ynanIDgBkT1gEAoFPCOgAAdEpYBwCATgnrAADQqQfMuwEAXdlqIulmZa0dXFsAWHl61gEAoFOjhPWqekpVvaaq3ldVX6yqVlVv3mNdJ1TV1VV1W1XdU1VHqurKqnroGG0F2FJr973tpAwADsBYw2BekuT7knwpyaeSnL6XSqrq5CQ3JDk+yduS3Jjk7CTPS3JxVZ3bWrt9lBYDAEDnxhoG84IkpyY5NslP76Oe12cS1C9vrV3SWvv51toFSV6V5LQkL993SwEAYEGMEtZba9e31m5pbe/fCVfVSUkuSnIkyes2FL80yV1JLq2qo/fcUAD2x9VcAWaqpwmmFwzb61pr964vaK3dmeQDSR6U5HGzbhgAAMxDT0s3njZsb55SfksmPe+nJnnnVhVV1aEpRXsaSw+sOBNJAZiTnnrWjxu2X5hSvnb/Q2bQFgAAmLueeta3szZIctsurtbaWZtWMOlxP3PMRgEAwEHpKayv9ZwfN6X82A37AXCQXM0VYO56GgZz07A9dUr5KcN22ph2AABYKj31rF8/bC+qqvutXxGmqo5Jcm6Su5N8cB6NA1g5m/WUr/Wo60WHcXltMcXMe9ar6oFVdfpwtdKva63dmuS6JCcmee6Gh12R5Ogk17bW7ppJQwEAYM5G6VmvqkuSXDL8+Ihhe05VXTP8+/OttRcN//7OJIeTfDyTYL7ec5LckOSqqrpw2O+xSc7PZPjLi8doLwAALIKxhsE8OskzNtx30nBLJsH8RdlGa+3WqnpMkl9KcnGSJyX5dJKrklzRWrtjpPYCAED3RgnrrbWXJXnZDvc9km8sw7hZ+SeTXDZGuwAAYJH1NMEUgN6Z/Ab7Z1lUdqGnpRsBWEZVW4cTAKbSsw4AMEuWRWUX9KwDAECnhHUAAOiUsA4AAJ0yZh2WiTGPzJtVLgBGJawDAMybD65MIawDMB6rXACMyph1AADolLAOAACdMgwGFpWJfACw9PSsAwBAp/Ssw6LqaSKfCYRsxfMCYM/0rAMAQKeEdQAA6JSwDgAAnRLWAQCgUyaYwjKZxUQ+S0YCwMzoWQcAgE7pWV8FltVjTD0tGQkAS07POgAAdEpYBwCATgnrAADQKWEdAHpWtfUqTMBSM8F02VhWj3nwPFo+Jg0DdEHPOgAAdErP+rKxrB4AwNLQsw4AAJ3Ssw4AvTDviEXn2/zRCetAX5zoZ09ABOiWsA4AvTDvCNhAWF8FTvDAVgREgG6ZYAoAAJ3Ssw4AwO6Z7zITwjowP070ALAlYR0AeuZDKr0y32UmhHVgfpzo++XvD9AFE0wBAKBTwjoAAHRKWAcAgE4Zsw4AwDjMdxmdsA70xYkeAL5utGEwVXVCVV1dVbdV1T1VdaSqrqyqh+6ynn9QVW8bHv+VqvpEVb29qi4eq60AALAIRgnrVXVykkNJLkvyoSSvSvKxJM9L8kdV9fAd1vPTSd6X5MJh+6ok70nyxCT/qapePEZ74T6qtr5ADwDAHIw1DOb1SY5Pcnlr7TVrd1bVK5O8IMnLkzx7qwqq6oFJfiXJV5Kc1Vq7aV3ZLyf5kyQvrqp/1Vq7Z6R2AwBAt/bds15VJyW5KMmRJK/bUPzSJHclubSqjt6mqoclOS7JzeuDepK01g4nuTnJUUkevN82AwDAIhhjGMwFw/a61tq96wtaa3cm+UCSByV53Db1fC7JXyc5tapOWV9QVacmOSXJR1prt4/QZgAA6N4YYf20YXvzlPJbhu2pW1XSWmtJnju06VBV/WZV/UpVXZvJePi/SPLUEdoLAAALYYwx68cN2y9MKV+7/yHbVdRae2tV3Zbk3yd5+rqizyZ5YyaTVrdVVYemFJ2+k8ez5LaaSLpZmaUEAYA5mcUVTNfSz7aJp6qeluQ/Z7ISzBmZDJ85I8k7k7w2yVsOqI0AANCdMXrW13rOj5tSfuyG/TY1jEu/OsmfJrl03fj3G6vq0kyG2zy1qs5rrb17q7paa2dN+R2Hkpy51WNZAZv1lK/1qOtFB4Dted+cmTF61tdWbpk2Jn1tsui0Me1rLkrywCTv2WSi6r1J3jv8uGkQBwCAZTNGWL9+2F5UVd9UX1Udk+TcJHcn+eA29XzrsP32KeVr9//tXhoJAACLZt9hvbV2a5LrkpyYyWou612R5Ogk17bW7lq7s6pOr6qNkz3fN2yfUlX/4/qCqnp0kqdkMu79XfttMwAALIKxrmD6nCQ3JLmqqi5McjjJY5Ocn8nwlxdv2P/wsP360huttQ9V1RuTXJbkv1bV7yb5eCYfAi5J8i1Jrmyt/cVIbQYAgK6NEtZba7dW1WOS/FKSi5M8Kcmnk1yV5IrW2h07rOonMxmb/swkP5jkmCRfTPL+JL/RWrMaDAfDBBkA2Jwlj+dqrJ71tNY+mUmv+E723fSoDxdGuma4AQDAShstrAPch6W9ABafJY/nahYXRQIWTdXWX3sCADMhrAMAQKeEdQDYK99CAQdMWAcAgE6ZYAqMw9JeAKvDOXxmhHVYdUI2AHRLWAfGsdulvSz7BQDbEtZh1QnSsDO+hQLmwARTAADolJ51ANgJ30IBc6BnHVaZNaIBoGt61oGDs9bbaKzv+PToAqwEYR24LwEQALogrAMHz1hfANgTYR0A9sqHTeCACeuwKowbB4CFI6wD9M4HLYCVJazDqjBuHAAWjrAOzIcPCDu3rB+0luH/AHDAXBQJAAA6JawDAECnhHUAAOiUMeuwyowVBoCuCesAi2jRPmhZfpLdMPkYvs4wGAAA6JSedQAO3rIuPwlwwPSsAwBAp4R1AADolGEwAMD8mHwMW9KzDgAAndKzDuycCYGMabfPI8+/5WTyMWxJzzoAAHRKWAcAgE4J6wAA0ClhHQAAOmWCKbA5y6kxT55/q83xhK/Tsw4AAJ3Ssw5sznJqzJPnH0ASPesAANAtYR0AADolrAMAQKeEdYBZqdp6lROAnjhndcEEU2DnTOxjnjz/gBU0Ws96VZ1QVVdX1W1VdU9VHamqK6vqoXuo63ur6tqq+uRQ1+eq6j1V9fSx2gsAAL0bpWe9qk5OckOS45O8LcmNSc5O8rwkF1fVua2123dY1zOTvCHJl5P8fpIjSR6S5FFJnpTk2jHaDAAAvRtrGMzrMwnql7fWXrN2Z1W9MskLkrw8ybO3q6SqHpdJUP/zJBe31j6zofyBI7UXAJaT9ehhqex7GExVnZTkokx6wF+3ofilSe5KcmlVHb2D6n41yf2TPG1jUE+S1tpX99dagBlZm5i1/raTMoB5cM7q1hg96xcM2+taa/euL2it3VlVH8gkzD8uyTunVVJVJyT5gSR/nOQvqur8JGclaUk+kuT6jfUDAMAyGyOsnzZsb55SfksmYf3UbBHWk/z9dfu/K8l5G8r/rKqe3Fr76HYNqqpDU4pO3+6xsJR8LT57m/2tHQdgFvZyrnHO6tYYq8EcN2y/MKV87f6HbFPP8cP2x5OckeTJQ93fk+RNSb43yR9U1bfsvakAALA4ZrHO+trApu0+lt1/3fZZrbXfH37+YlU9I5MA/5gkP5bk329VUWvtrE0bMulxP3MnjQaA7m01dnizMj2ksHDG6Flf6zk/bkr5sRv2m+Zvhu09Sd6+vqC11jJZEjKZLAkJAABLb4ye9ZuG7alTyk8ZttPGtG+s584pE0nXwvxRu2gbLC5jBYHtGGcMS2+MsH79sL2oqu63PmhX1TFJzk1yd5IPblPPnyb5fJJvq6rvaK19dkP5o4btkf03GZaYr8X75W8NjO0gz/nOWV3Y9zCY1tqtSa5LcmKS524oviLJ0Umuba3dtXZnVZ1eVd+0Mktr7WtJ/s3w469W1f3W7f+9SZ6Z5GtJfnu/bQYAgEUw1gTT5yS5IclVVXVhksNJHpvk/EyGv7x4w/6Hh+3Gj3y/nOTCJE9P8r1V9e4k357JpNK/l+SFO1m6EVaar8UBVodz/tIbY4LpWu/6Y5Jck0lIf2GSk5NcleSc1trtO6zny5mE9SuSPCiTnvofyeSDwJNaa68co70AALAIRlu6sbX2ySSX7XDfqQOshsD+suEGq8E4c2Aszg+wVEbpWWdBVW0dEgEAmKtZXBSJRWOs2+wZcwgAbEJYh1Ug8AOsDuf8pWIYDAAAdEpYBwCAThkGsyp2u9oIAABzJ6wznSUD58vfGgBWnrC+KrZabQQAWDzLvGrYMv/fdklY574vBC8QANgZ75kcMBNMAQCgU8I6AAB0SlgHAIBOGbO+yoyvA4DFsNslmBfpPX6Z/28jENa5rxV7EQDAjgiVzIGwDjA2q0MAY9tqCeZFP9cs8/9tBMI6ADCOZQ9YQiVzYIIpAAB0SliHWaly1VgAYFeEdQAA6JQx6wD7YXUIYF6W+XyyzP+3XRLWAYDd80F1Yln/X3RDWAfYD6tDAHCAhHU4CHqcgGXngyrMhAmmAADQKT3rcBD0OAEAI9CzDgAAndKzDjA2354AMBJhHQAYhw+qMDrDYAAAoFN61mFW9DgBALukZx0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAPAIqqa3IClJqwDAECnhHUAAOiUsA4AAJ0aLaxX1QlVdXVV3VZV91TVkaq6sqoeuo86n1BVf1dVrar+xVhtBQCARfCAMSqpqpOT3JDk+CRvS3JjkrOTPC/JxVV1bmvt9l3WeUyS30zy5SQPHqOdAFOtTdRrbb7tgM1sNZF0szLPY1gaY/Wsvz6ToH55a+2S1trPt9YuSPKqJKclefke6nx1kuOS/MpIbQQAgIWy77BeVScluSjJkSSv21D80iR3Jbm0qo7eRZ0/muSyJJcnuW2/bQRYGZbzW06t3fe2kzJg4Y3Rs37BsL2utXbv+oLW2p1JPpDkQUket5PKqur4JL+R5Pdaa28eoX0AALCQxhizftqwvXlK+S2Z9LyfmuSdO6jv32byIeLZe21QVR2aUnT6XusEAIBZGyOsHzdsvzClfO3+h2xXUVX9RJIfTfK/tNY+O0LbAO7LZD0AFsQoq8FsY+2db8t3u6o6McmVSd7aWvsP+/mFrbWzpvyOQ0nO3E/dAAAwK2OE9bWe8+OmlB+7Yb9prk5yd5LnjNAmgOk26ylfxKUbfUOw2hxPWAljTDC9adieOqX8lGE7bUz7mjMzWf7xr4eLILWqakneOJS/eLjv9/bXXAAAWAxj9KxfP2wvqqr7rV8RZriw0bmZ9Jh/cJt6rs1k1ZiNTknyhCQfSXIoyZ/su8UAy2BZviFYJY4P8+T5t5D2HdZba7dW1XWZrPjy3CSvWVd8RZKjk/yb1tpda3dW1enDY29cV8/lm9VfVc/MJKz/QWvtJfttLwAALIqxJpg+J8kNSa6qqguTHE7y2CTnZzL85cUb9j88bF25AwAAphhjzHpaa7cmeUySazIJ6S9McnKSq5Kc01q7fYzfA3BgXPkRgA6NtnRja+2TSS7b4b477lFvrV2TyYcAAABYKbNYZx2AWfHtQD8srck8ef4tjVGGwcBoqrY+wQAArBA96wBwECytyTx5/i0NPeuwHb39sDq83oHOCOsAANApYR1g3vTmAjCFMevMj5nqAABbEtYBYFZ0OjBPnn8LSVhnfnqcqa63H1aH1zuwAIxZBwCATulZh/V67O1nuejN7YfXO8yO19ae6VkHAIBO6VkHmCW9uQDsgrBOX4QVAICvMwwGYK9czAiWj9c1ndGzznIac1iB3n5YHV7vi8lQsr6YSD8qPesAMG96c4Ep9KwDzJteJWCZmEg/Kj3rAADQKT3rADthDCYsH69rFoCwzmLbboynky3AchGwx2eISteEdYCdMAaTsQib/fC6ZgEI6yyHjSdVJ1uA5SRgLybHZs+EdQCYJWET2AWrwQAAQKf0rLM4djvOEwD4ZuZMLBxhHWCvvInB8vG6pjPCOovDOE8ANuM9YOe8ly4cYR0A5k1IAqYwwXTRVBmfDQCLyvv4wVjiv6uwDiy23Z6gl/iEDsDyEdYBAKBTxqyz2IzzBID98V7aNT3r7J9hBQAAB0LPes9cuAAAFpf38YOxYn9XYR1YHGNexXYJT+gALB9hvWcuXAAAi8v7+MFYsb+rsA4sjjFO0Et8Qgdg+Qjr7M6KjRMDAJgnq8EAAECn9KyzOys2TgwAYJ6E9UUjEAPA4vI+fjCW+O8qrLMY9N4zzW6fE55DACyQ0casV9UJVXV1Vd1WVfdU1ZGqurKqHrrDxx9dVf9bVf1fVXVjVd1VVXdW1R9X1Qur6lvGaisAACtiwa+0PkrPelWdnOSGJMcneVuSG5OcneR5SS6uqnNba7dvU80PJHlzkjuSXJ/k95I8LMkPJ/lXSZ5cVRe21r4yRpsBAKB3Yw2DeX0mQf3y1tpr1u6sqlcmeUGSlyd59jZ1fCbJ05K8tbX2t+vqOCbJu5M8Pslzk/zaSG1mLIYVAAAciH0Pg6mqk5JclORIktdtKH5pkruSXFpVR29VT2vtI62131of1If778w3Avp5+20vABy4Bf/aHejHGGPWLxi217XW7l1fMATtDyR5UJLH7eN3fHXYfm0fdbAo1t7k1t92UgYAsGTGCOunDdubp5TfMmxP3cfv+Ilh+4f7qAMAgGW2hB1+Y4xZP27YfmFK+dr9D9lL5VX1M0kuTvKRJFfv8DGHphSdvpc2MGMuvAQAkGQ266yvfWzZdcqqqicnuTKTyac/1lr76jYPYREI3gCz1+O5t8c2sdiWsMNvjLC+1nN+3JTyYzfstyNVdUmStyT5XJLzW2sf2+ljW2tnTanzUJIzd9MO6M6Cn3Rg6Wz1VfpmZV67wC6MMWb9pmE7bUz6KcN22pj2+6iqpyZ5a5LPJnlia+2mbR6yuBZozBTAQnF+xXOAJTBGz/r1w/aiqrrf+hVhhjXSz01yd5IP7qSyqvonSa5N8t+yyx51AJi5JfzaHejHvnvWW2u3JrkuyYmZXLRovSuSHJ3k2tbaXWt3VtXpVXWfyZ5V9Ywkb0ryiSRPENT5uta86QEAK2esCabPSXJDkquq6sIkh5M8Nsn5mQx/efGG/Q8P269/N1VV52ey2sv9Mumtv6zu+9XVf2+tXTlSm5kFYzkBZq/Hc2+PbWI1LPhzaZSw3lq7taoek+SXMllm8UlJPp3kqiRXtNbu2EE1351v9PT/xJR9Pp7J6jCLx1ei7IU3NwBYaaMt3dha+2SSy3a4731SRmvtmiTXjNWebk0LX4sQvPbygcNYTujbsrwefbD9Zj2eew+6TZ4DLKlZrLMO7FWPb7jA9rw+gZEI67O28QQueAHsnw+2eA6wpMZYZx0AgPWs8b4z/k7b0rN+EHY7bg4AADYhrLO5g5yo4+tI2Jmxv8I3AW+19Xg8e2wTdEZYPwi7GTenp53d8uYGACtDWJ+3XoOXiTqwfFbtdb2M/yd2x3OAJSCsA6yCZQ7ls+Dvx1YMMdsZf6c9sRoMAIvNahLAEtOzDtADPU6wuFZtiNle+TvtibA+Kwf9JPRkBwBYOsI6O+eDAPO07B9I59XjtKx/z60s+3MJWCrCOsCymceQmmUKwIYkAR0R1gFYHII0sGKEdYBlYxLX/vj7MQbPlZ3xd9qWsL6IFqVnyZsbMDZBGlgxwjrQn0X5QHrQlvX/tZ0xw7fnUh/GOqY+mLGChPVFpGcJmCcBGJiXFcw7wjrQHx9IGYvnErDghHWAVTBmMF3FALys/y92btmf43RLWGcci/K1uJMtLB+vZ2CJCeu9EiqBVbAoH/TZubGOqecGJBHWl8e8T1Kr+LU4AHBwfGBLIqwDi2JJT8Irbx4f9D2XDtZYx1QnECQR1gEYg/DEMtGj2wcf2JII6ywzJ1sAYMEJ6z0QKlfbCvYSAHRNjy4dEdY5OPM+oTnZAgALTljvgVAJ8A3Oe8tnrGPqucEKEtYBAFgMK/iB7X7zbgAAALA5Peswtq2GMJlMDLCYnI+ZE2Gd1eJkC0BibhgLQ1jvlZPHcjKZGADYBWPWAQCgU8I6AAB0yjAY2A8TRgGAAySsAwDLbdE7VsxtWmnCOuzHGBNGnXwBgCmEdQDgvpapN9dKXCwwE0wBAKBTwjoAAHTKMBjYK1+hAjC2RZ8My+iEdRibEycAMJLRhsFU1QlVdXVV3VZV91TVkaq6sqoeust6HjY87shQz21DvSeM1VZgBVRt3UMFfMPa62X9bSdli6y1PjtX1tq1/raTMpbWKD3rVXVykhuSHJ/kbUluTHJ2kuclubiqzm2t3b6Deh4+1HNqkncleUuS05NcluQfVdU5rbWPjdFmAADo3Vg966/PJKhf3lq7pLX28621C5K8KslpSV6+w3p+OZOg/qrW2oVDPZdkEvqPH34PADAmvbnQrX2H9ao6KclFSY4ked2G4pcmuSvJpVV19Db1HJ3k0mH/l24ofu1Q/w8Ovw9maxW/IgYA5m6MnvULhu11rbV71xe01u5M8oEkD0ryuG3qOSfJUUk+MDxufT33Jrlu+PH8fbcYAAAWwBhj1k8btjdPKb8lk573U5O8c5/1ZKhnS1V1aErR6ds9Fjbl6nd9s9QZsMycs1baGD3rxw3bL0wpX7sVfNaqAAAKt0lEQVT/ITOqBwAAlsIs1llf69ba78fCHdfTWjtr0womPe5n7rMdQG988wHj89qBLozRs77W433clPJjN+x30PUAAMBSGCOs3zRsp40lP2XYThuLPnY9AACwFMYI69cP24uq6pvqq6pjkpyb5O4kH9ymng8O+507PG59PffLZJLq+t8H82W9YQDggO07rLfWbs1kWcUTkzx3Q/EVSY5Ocm1r7a61O6vq9Kr6ppVZWmtfSvKmYf+XbajnZ4b63+EKpgAArIqxJpg+J8kNSa6qqguTHE7y2EzWRL85yYs37H942G5cU+0XkpyX5J9V1aOTfCjJGUl+NMnnct8PAwCb860HAEtgjGEwa73rj0lyTSYh/YVJTk5yVZJzWmu377Ce2zO5ONJVSb5nqOexSd6Y5Kzh9wAAwEoYbenG1tonk1y2w32nXsGktXZHkucNNwAAWFmj9KwDAADjE9YBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBT1VqbdxtmpqpuP+qoox52xhlnzLspAAAsscOHD+fuu+++o7X28P3Us2ph/a+SHJvkyBx+/enD9sY5/G5my7FeHY716nCsV4djvToO+lifmOSLrbVH7qeSlQrr81RVh5KktXbWvNvCwXKsV4djvToc69XhWK+ORTnWxqwDAECnhHUAAOiUsA4AAJ0S1gEAoFPCOgAAdMpqMAAA0Ck96wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWD9gVXVCVV1dVbdV1T1VdaSqrqyqh867bexcVT28qp5VVb9bVR+tqrur6gtV9f6q+smq2vS1VFWPr6q3V9UdVfXlqvrTqnp+Vd1/1v8H9qeqLq2qNtyeNWWff1xV7x6eG1+qqv9SVc+YdVvZvar6gar6nar69HCu/nRVXVdVT9pkX6/rBVVV/2g4rp8azuMfq6q3VtU5U/Z3rDtWVU+pqtdU1fuq6ovD+fnN2zxm18d03ud2F0U6QFV1cpIbkhyf5G1JbkxydpLzk9yU5NzW2u3zayE7VVXPTvKvk3w6yfVJPpHkO5I8OclxSX4nyVPbuhdUVf3ocP9XkvzfSe5I8sNJTkvy2621p87y/8DeVdV3JfmzJPdP8uAkP9Vae8OGfX4myWuS3J7J8f7bJE9JckKSX2utvWimjWbHquolSf6PJJ9P8vuZvM6/Lcn3J7m+tfaz6/b1ul5QVfUvk/xsJq/R38vkeH9Pkh9J8oAkT2+tvXnd/o5156rqI0m+L8mXknwqyelJfqu19rQp++/6mHZxbm+tuR3QLck7krQk/3TD/a8c7v/1ebfRbcfH8oJMXtD323D/IzIJ7i3Jj627/9gkn0tyT5LHrLv/72XyAa4l+V/n/f9y29GxryT/OcmtSf7P4dg9a8M+J2Zy8r89yYnr7n9oko8Ojzln3v8Xt02P71OH4/P/Jjlmk/IHrvu31/WC3oZz9d8l+UyS4zeUnT8cu4851ot1G47dKcN5+rzhuLx5yr67Pqa9nNsNgzkgVXVSkouSHEnyug3FL01yV5JLq+roGTeNPWitvau19h9ba/duuP8zSX59+PG8dUVPSfLtSd7SWvvjdft/JclLhh9/+uBazIguz+TD2mWZvG438xNJvjXJa1trR9bubK39TZJfHn589gG2kT0Yhq/9yyRfTvJPWmt3btyntfbVdT96XS+u785k6O9/aa19bn1Ba+36JHdmcmzXONYLoLV2fWvtljYk6G3s5Zh2cW4X1g/OBcP2uk0C3p1JPpDkQUkeN+uGMbq1N/Ovrbtv7fj/4Sb7vzeTcPD4qvrWg2wY+1NVZyR5RZJXt9beu8WuWx3v/7RhH/rx+CSPTPL2JH8zjGf+uap63pQxzF7Xi+uWTIYvnF1V37a+oKqekOSYTL5BW+NYL5+9HNMuzu3C+sE5bdjePKX8lmF76gzawgGpqgckefrw4/oX89Tj31r7WpK/ymSM5EkH2kD2bDi2b8pkmNMvbLP7Vsf705n0yJ9QVQ8atZHs198ftp9N8uFMxqu/IsmVSW6oqvdU1freVq/rBdVauyPJz2Uy1+gvq+rfVtWvVNV/SHJdJsOg/vd1D3Gsl89ejmkX53Zh/eAcN2y/MKV87f6HzKAtHJxXJHlUkre31t6x7n7Hf/H980wmGD6ztXb3Nvvu9HgfN6Wc+Th+2D47yVFJ/mEmPayPymTO0ROSvHXd/l7XC6y1dmUmiwI8IMlPJfn5TOYsfDLJNRuGxzjWy2cvx7SLc7uwPj81bC3Hs6Cq6vIkL8xklZ9Ld/vwYev4d6iqzs6kN/3XWmt/NEaVw9bx7svaUm2V5CmttXe21r7UWvuLJP9zJqtLPHHasn6bcJw7VlU/m+S3k1yT5OQkRyc5K8nHkvxWVf3qbqobto718tjLMZ3J80BYPzjbfdo6dsN+LJCqem6SVyf5yyTnD1+xruf4L6h1w19uTvKLO3zYTo/3F/fRNMb3N8P2Y621/299wfBtytq3ZWcPW6/rBVVV52Uymfj/aa39s9bax1prX26tfTiTD2b/LckLh8UhEsd6Ge3lmHZxbhfWD85Nw3bamPRThu20Me10qqqen+S1Sf48k6D+mU12m3r8hzD4yEwmpH7soNrJnj04k+N2RpKvrLsQUstkJack+Y3hviuHn7c63v9DJj14n2qtffmA287urB23/z6lfC3MH7Vhf6/rxfOPh+31GwuG1+WHMslE3z/c7Vgvn70c0y7O7cL6wVk7IVy08eqWVXVMknOT3J3kg7NuGHtXVT+X5FVJPpJJUP/clF3fNWwv3qTsCZmsBHRDa+2e8VvJPt2T5N9Nuf3JsM/7h5/Xhshsdbx/aMM+9OO9mbw5n1JV37JJ+aOG7ZFh63W9uNZW+Pj2KeVr9//tsHWsl89ejmkf5/Z5L2i/zLe4KNJS3TIZEtGS/HGSh22z77FJ/jouqLFUtyQvy+YXRXpkOrhwhtuejumbh+PzLzbc/z8luTeTXveHDPd5XS/oLcmPD8fnM0m+c0PZDw3H+u4kD3esF/OWnV0UaVfHtJdzew2/lANQVSdn8gQ4PsnbkhxO8thMrrh1c5LHt9Zun18L2amqekYmk5L+LpPLDm82TvFIa+2adY+5JJPJTF9J8pZMLmv8Ixkua5zkx5sX4EKpqpdlMhTmp1prb9hQ9k+TXJV5XpKaXauq4zO57sX3JHlfJsMhvjuTccwtk4slvXXd/l7XC2j4hvsdmaz4c2eS380kuJ+RyRCZSvL81tqr1z3Gse7ccIwuGX58RJIfzGQYy/uG+z6//ty7l2Paxbl93p+Elv2W5LuSvDHJp4cD/PFMJiZu2TPr1tct3+hR3er27k0ed26GC65k0mvzZ0lekOT+8/4/ue3refCsKeU/nOQ9mYSBu5L81yTPmHe73bY9rg/L5BvPvxrO07dn0sHyuCn7e10v4C3JA5M8P5Php1/MZAjU5zJZX/8ix3rxbjt4bz4yxjGd97ldzzoAAHTKBFMAAOiUsA4AAJ0S1gEAoFPCOgAAdEpYBwCATgnrAADQKWEdAAA6JawDAECnhHUAAOiUsA4AAJ0S1gEAoFPCOgAAdEpYBwCATgnrAADQKWEdAAA6JawDAECnhHUAAOjU/w9EFL3aoyYOYgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 250,
       "width": 373
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(rand_test, 'r+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_exp = 24\n",
    "expId0 = 100\n",
    "elevation_min = 30.\n",
    "elevation_flat = 70.\n",
    "elevation_max = 85.\n",
    "rng = np.random.RandomState(15)\n",
    "az = []\n",
    "el = []\n",
    "airmasses = []\n",
    "for expId in range(expId0, expId0 + n_exp):\n",
    "    azimuth = rng.rand()*360.\n",
    "    elevation = quad_distribution(rng, elevation_min, elevation_flat, elevation_max)[1]\n",
    "    az.append(azimuth)\n",
    "    el.append(elevation)\n",
    "    _airmass = 1./np.cos((90 - elevation)*np.pi/180.)\n",
    "    airmasses.append(_airmass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from lsst.sims.photUtils import matchStar\n",
    "\n",
    "imp.load_source('calc_refractive_index', '/Users/sullivan/LSST/code/StarFast/calc_refractive_index.py')\n",
    "imp.load_source('fast_dft', '/Users/sullivan/LSST/code/StarFast/fast_dft.py')\n",
    "imp.load_source('StarFast', '/Users/sullivan/LSST/code/StarFast/StarFast.py')\n",
    "from StarFast import StarSim\n",
    "\n",
    "lsst_lat = Angle(np.radians(-30.244639))\n",
    "lsst_lon = Angle(np.radians(-70.749417))\n",
    "lsst_alt = 2663."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the general parameters to use for all of the simulations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 8  # Seed for the random number generator. Simulations with the same seed are reproduceable\n",
    "dimension = 1024  # Number of pixels on a side for the simulated image\n",
    "n_star = 10000  # Number of stars to model in the simulated catalog.\n",
    "# The catalog covers an area ~4x larger than the area, to allow simulated rotations and dithering observations\n",
    "pixel_scale = 0.25  # plate scale, in arcseconds/pixel\n",
    "hottest_star = 'B'  # Hottest star to include (types are 'OBAFGKMR')\n",
    "coolest_star = 'M'  # Coolest star to include\n",
    "wavelength_step = 10  # Wavelength resolution of the spectra and calculation of filter and DCR effects. In nm.\n",
    "ra_offset = Angle(0)  # Additional offset in RA from the field center, for dithering. In radians as an LSST Angle object\n",
    "dec_offset = Angle(0)  # Additional offset in Dec from the field center, for dithering. In radians as an LSST Angle object\n",
    "sky_rotation = 0.0  # Sky rotation angle, in Degrees. I realize this is different than RA and Dec\n",
    "exposureId = 0  # Unique exposure identification number. Also used as the \"OBSID\"\n",
    "instrument_noise = 0.  # Adds noise akin to instrumental noise (post-PSF). Set to 1.0 for default value, can be scaled up or down\n",
    "photon_noise = 1./15.  # Adds poisson noise akin to photon shot noise. Set to 1.0 for default value, can be scaled up or down\n",
    "sky_noise = 0  # Adds noise prior to convolving with the PSF.\n",
    "output_directory = \"/Users/sullivan/LSST/simulations/test%1i_variablePSF/\" % seed  # Path to the directory to save output\n",
    "band_dict = {'u': 0, 'g': 1, 'r': 2, 'i': 3, 'z': 4, 'y': 5}  # LSST filter numbers used by the butler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The \"base\" PSF is used, e.g.,  for setting the buffer size to avoided sources right at the edge\n",
    "# It should be larger than most PSFs used\n",
    "psf_fwhm_min = 0.25 - 0.0625\n",
    "psf_fwhm_base = 0.75  # FWHM of the PSF, in arcseconds\n",
    "psf_fwhm_variance = 0.125"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "attenuation = 20. # attenuation factor that was used in the simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra = lsst_lon + ra_offset\n",
    "dec = lsst_lat + dec_offset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load simulated Kurucz stellar SEDs. Cache the list for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "pickle_file = \"sed_list.pickle\"\n",
    "if os.path.exists(pickle_file):\n",
    "    with open(pickle_file, 'rb') as dumpfile:\n",
    "        sed_list = pickle.load(dumpfile)\n",
    "else:\n",
    "    matchStarObj = matchStar()\n",
    "    sed_list = matchStarObj.loadKuruczSEDs()\n",
    "    with open(pickle_file, 'wb') as dumpfile:\n",
    "        pickle.dump(sed_list, dumpfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the PSF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "gsp = galsim.GSParams(folding_threshold=1.0 / (dimension), maximum_fft_size=12288)\n",
    "psf_base = galsim.Kolmogorov(fwhm=psf_fwhm_base / pixel_scale, flux=1, gsparams=gsp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# g-band simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "band_name = 'g'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_bp = True\n",
    "sim = StarSim(psf=psf_base, pixel_scale=pixel_scale, x_size=dimension, y_size=dimension,\n",
    "              band_name=band_name, wavelength_step=wavelength_step,\n",
    "              sed_list=sed_list, ra=ra, dec=dec, sky_rotation=sky_rotation,\n",
    "              use_mirror=use_bp, use_lens=use_bp, use_atmos=use_bp,\n",
    "              use_filter=use_bp, use_detector=use_bp, attenuation=attenuation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate a catalog of stars, with fluxes and SEDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number and flux contribution of stars of each type:\n",
      " [M 7745| 1.21%] [K 1164| 1.99%] [G 716| 2.45%] [F 300| 3.59%] [A 62| 9.34%] [B 13| 81.42%] [O 0| 0.00%]\n"
     ]
    }
   ],
   "source": [
    "sim.load_catalog(n_star=n_star, hottest_star=hottest_star, coolest_star=coolest_star, seed=seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the raw simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing 695 stars brighter than 16.0 mag to reference catalog in 1 bands\n",
      "Min/max magnitude:  5.724159303122955 22.49404801300464\n"
     ]
    }
   ],
   "source": [
    "sim.make_reference_catalog(output_directory=output_directory + \"input_data/\", filter_list=['g'], magnitude_limit=16.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulating 2428 stars within observable region\n",
      "Time to model 2421 stars: [39.672s | 0.01639s per star]\n",
      "Time to model 7 bright stars: [4.945s | 0.70650s per star]\n"
     ]
    }
   ],
   "source": [
    "sim.simulate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_exp = 24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FFT timing for 15 DCR planes: [2.194s | 0.146s per plane]\n",
      "FFT timing for 15 DCR planes: [8.169s | 0.545s per plane]\n",
      "FFT timing for 15 DCR planes: [1.741s | 0.116s per plane]\n",
      "FFT timing for 15 DCR planes: [7.027s | 0.468s per plane]\n",
      "FFT timing for 15 DCR planes: [1.853s | 0.124s per plane]\n",
      "FFT timing for 15 DCR planes: [8.211s | 0.547s per plane]\n",
      "FFT timing for 15 DCR planes: [1.984s | 0.132s per plane]\n",
      "FFT timing for 15 DCR planes: [8.567s | 0.571s per plane]\n",
      "FFT timing for 15 DCR planes: [1.923s | 0.128s per plane]\n",
      "FFT timing for 15 DCR planes: [7.553s | 0.504s per plane]\n",
      "FFT timing for 15 DCR planes: [1.810s | 0.121s per plane]\n",
      "FFT timing for 15 DCR planes: [7.372s | 0.491s per plane]\n",
      "FFT timing for 15 DCR planes: [2.163s | 0.144s per plane]\n",
      "FFT timing for 15 DCR planes: [7.692s | 0.513s per plane]\n",
      "FFT timing for 15 DCR planes: [1.960s | 0.131s per plane]\n",
      "FFT timing for 15 DCR planes: [8.591s | 0.573s per plane]\n",
      "FFT timing for 15 DCR planes: [2.530s | 0.169s per plane]\n",
      "FFT timing for 15 DCR planes: [7.617s | 0.508s per plane]\n",
      "FFT timing for 15 DCR planes: [1.727s | 0.115s per plane]\n",
      "FFT timing for 15 DCR planes: [7.262s | 0.484s per plane]\n",
      "FFT timing for 15 DCR planes: [1.752s | 0.117s per plane]\n",
      "FFT timing for 15 DCR planes: [7.861s | 0.524s per plane]\n",
      "FFT timing for 15 DCR planes: [1.754s | 0.117s per plane]\n",
      "FFT timing for 15 DCR planes: [7.497s | 0.500s per plane]\n",
      "FFT timing for 15 DCR planes: [1.880s | 0.125s per plane]\n",
      "FFT timing for 15 DCR planes: [8.006s | 0.534s per plane]\n",
      "FFT timing for 15 DCR planes: [1.886s | 0.126s per plane]\n",
      "FFT timing for 15 DCR planes: [7.320s | 0.488s per plane]\n",
      "FFT timing for 15 DCR planes: [2.129s | 0.142s per plane]\n",
      "FFT timing for 15 DCR planes: [8.893s | 0.593s per plane]\n",
      "FFT timing for 15 DCR planes: [2.175s | 0.145s per plane]\n",
      "FFT timing for 15 DCR planes: [8.730s | 0.582s per plane]\n",
      "FFT timing for 15 DCR planes: [1.782s | 0.119s per plane]\n",
      "FFT timing for 15 DCR planes: [8.789s | 0.586s per plane]\n",
      "FFT timing for 15 DCR planes: [1.820s | 0.121s per plane]\n",
      "FFT timing for 15 DCR planes: [7.155s | 0.477s per plane]\n",
      "FFT timing for 15 DCR planes: [1.738s | 0.116s per plane]\n",
      "FFT timing for 15 DCR planes: [7.031s | 0.469s per plane]\n",
      "FFT timing for 15 DCR planes: [1.877s | 0.125s per plane]\n",
      "FFT timing for 15 DCR planes: [7.003s | 0.467s per plane]\n",
      "FFT timing for 15 DCR planes: [1.720s | 0.115s per plane]\n",
      "FFT timing for 15 DCR planes: [8.488s | 0.566s per plane]\n",
      "FFT timing for 15 DCR planes: [1.726s | 0.115s per plane]\n",
      "FFT timing for 15 DCR planes: [7.636s | 0.509s per plane]\n",
      "FFT timing for 15 DCR planes: [1.918s | 0.128s per plane]\n",
      "FFT timing for 15 DCR planes: [7.096s | 0.473s per plane]\n",
      "FFT timing for 15 DCR planes: [1.783s | 0.119s per plane]\n",
      "FFT timing for 15 DCR planes: [7.304s | 0.487s per plane]\n"
     ]
    }
   ],
   "source": [
    "expId0 = exposureId + 100*band_dict[band_name]\n",
    "elevation_min = 30.\n",
    "elevation_flat = 70.\n",
    "elevation_max = 85.\n",
    "rng = np.random.RandomState(15)\n",
    "for expId in range(expId0, expId0 + n_exp):\n",
    "    azimuth = 0 if rng.rand() < 0.5 else 180.\n",
    "    elevation = quad_distribution(rng, elevation_min, elevation_flat, elevation_max)[1]\n",
    "    psf_fwhm_use = psf_fwhm_min + np.abs(psf_fwhm_variance*rng.randn())\n",
    "    psf = galsim.Kolmogorov(fwhm=psf_fwhm_use/pixel_scale, flux=1, gsparams=gsp)\n",
    "    exposure = sim.convolve(elevation=elevation, azimuth=azimuth,\n",
    "                            instrument_noise=instrument_noise, sky_noise=sky_noise,\n",
    "                            photon_noise=photon_noise, exposureId=expId, obsid=expId,\n",
    "                            psf=psf)\n",
    "    filename = \"lsst_e_%3.3i_f%i_R22_S11_E000.fits\" % (expId, band_dict[band_name])\n",
    "    exposure.writeFits(output_directory + \"images/\" + filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}