{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physical Units Demo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates the abilites of the `PhysicalFresnelWavefront` class.\n",
    "\n",
    "The most important change is that the `wavefront` array now contains values that have units of an actual electric field, i.e. V/m. This is of special importance once non-linear effects are of interest. For example, the thermal blooming effect results from a high-energy laser which heats the air and induces refractive-index changes. For this effect to be implemented, one needs to have the intensity distribution in physical units at disposal. The class `PhysicalFresnelWavefront` is designed to provide exactly this.\n",
    "\n",
    "An additional functionality is a rudimentary ability to emulate a beam quality factor. This has been done by scaling the wavelength accordingly. The capability to measure the beam quality factor is also implemented. This is useful in cases of arbitrary initial field distributions, where it might not be known."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demo 1: Beam Power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The beam power is very easy to apply. However, be careful when to apply a power, since it may result in unexpected behaviour. When the intensity is scaled before applying an aperture for example, the power will change. Let us begin by importing relevant packages and defining some constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import poppy\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.units as u\n",
    "\n",
    "w0 = 10e-2              # beam radius (m)\n",
    "P0 = 10e3               # beam power (W)\n",
    "w_extend = 6            # weight of the spatial extend\n",
    "wavelength = 1064e-9    # wavelength in vacuum (m)\n",
    "R0 = 100                # initial radius of curvature (m)\n",
    "npix = 256              # spatial resolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets define a wavefunction and apply an aperture, lens, and power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P = 10.00 kW\n"
     ]
    }
   ],
   "source": [
    "wf = poppy.PhysicalFresnelWavefront(beam_radius=w_extend*w0*u.m,\n",
    "                                    wavelength=wavelength,\n",
    "                                    npix=npix,\n",
    "                                    oversample=2)\n",
    "\n",
    "wf *= poppy.GaussianAperture(w=w0*u.m)\n",
    "wf *= poppy.QuadraticLens(f_lens=R0*u.m)\n",
    "\n",
    "wf.scale_power(P0)\n",
    "\n",
    "print(\"P = {0:.2f} kW\".format(wf.power*1e-3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let us see what happens when we apply the power before the aperture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P = 0.11 kW\n"
     ]
    }
   ],
   "source": [
    "wf = poppy.PhysicalFresnelWavefront(beam_radius=w_extend*w0*u.m,\n",
    "                                    wavelength=wavelength,\n",
    "                                    npix=npix,\n",
    "                                    oversample=2)\n",
    "\n",
    "wf.scale_power(P0)\n",
    "wf *= poppy.GaussianAperture(w=w0*u.m)\n",
    "wf *= poppy.QuadraticLens(f_lens=R0*u.m)\n",
    "\n",
    "print(\"P = {0:.2f} kW\".format(wf.power*1e-3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power values differ from the first result because the power is applied and calculated using the current intensity distribution, which differ as well. So always take care of when to apply the power."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demo 2: Beam Quality Factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can assign a beam quality factor to the beam which only scaled the wavelength appropriately. However, from the analytical equation for the beam radius one can see that this should work as expected. Let us test it by using the implemented beam quality measurement calculation based on ISO Standard 11146 (see https://www.rp-photonics.com/beam_quality.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M2 = 3.83\n"
     ]
    }
   ],
   "source": [
    "wf = poppy.PhysicalFresnelWavefront(beam_radius=w_extend*w0*u.m,\n",
    "                                    wavelength=wavelength,\n",
    "                                    npix=npix,\n",
    "                                    oversample=2,\n",
    "                                    M2=3.827)\n",
    "wf *= poppy.GaussianAperture(w=w0*u.m)\n",
    "wf *= poppy.QuadraticLens(f_lens=30*u.cm)\n",
    "wf.scale_power(P0)\n",
    "\n",
    "M2, z, caustic, z_fine, w_fit, rayleigh_length = wf.M2()\n",
    "\n",
    "print(\"M2 = {0:.2f}\".format(M2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I did not test the algorithm to very large extent, but it worked for my purposes. The results of the $M^2$ fitting routine is shown along with the values used to obtain the fit (red dots) in the following Figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x24ec58864a8>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XlclXX+///HG1kE1MrcNTjVT0TQxCXNLIfSscYaNZ2b02SNaUbLTElapuXSpuUKIaLilgmmljVTqKnlNq4hKiooruCWqTUfNcH1vL9/gP0cAj1Hzjnvs7zut5s3D9e5ONeTC67Xuc7rel/XpbTWCCGE8Bx+pgMIIYSwjxRuIYTwMFK4hRDCw0jhFkIIDyOFWwghPIwUbiGE8DBSuIUQwsNI4RZCCA8jhVsIITyMvzNetEaNGtpisTjjpYUQwitlZWWd0lrXtGVepxRui8XC5s2bnfHSQgjhlZRSBbbOK60SIYTwMFK4hRDCw0jhFkIIDyOFWwghPIwUbiGE8DDuU7jT08FiAT+/4v/T000nEkII27i4fjllOKDd0tMhLg4KC4u/Ligo/hqgVy9zuYQQ4kYM1C/ljFuXtWrVSts1jttiKf5hSwsPh/x8R8USQgjHc1D9Ukplaa1b2TKve7RKDh2yb7oQQrgLA/XLPQp3WJh904UQwl0YqF82FW6lVH+l1E6lVI5SKt7hKUaOhJCQ/5mkg4OLpwshhBs7PWgQhaUnhoQ4tX7dsHArpZoAzwOtgWbA40qphg5N0asXpKZCeDhaKfKBuQ89JAcmhRBub0BWFi/4+XGpXj1Qqri3nZrq1Pplyx53Y2Cj1rpQa30ZWA084fAkvXpBfj7KauWD557j2WXL2L17t8MXI4QQjpKZmcnMmTOpM2AAAUePgtVafEDSyTudthTunUB7pdTtSqkQoDNwhzNDjRo1ipCQEOLj43HGqBchhKgoq9XKq6++Su3atRk2bJhLl33Dwq213gWMBpYD3wLZwOXS8yml4pRSm5VSm0+ePFmhULVq1eKdd95h6dKlZGRkVOi1hBDCGdLT09m4cSMfffQR1apVc+my7R7HrZQaBRzRWqeUN4/d47jLcOnSJZo1a8bFixfJyckhKCioQq8nhBCOcvbsWSIiIggLC2PDhg34+VV8gJ7Dx3ErpWqV/B8GdAc+u/l4tgkICCAxMZH9+/eTkJDg7MUJIYTNRo4cyfHjx0lKSnJI0baXrUtcqJTKBb4B/qG1/q8TM/2mU6dOdO3alQ8++ICjR4+6YpFCCHFde/fuJSEhgd69e9OmTRsjGWwq3FrrB7XWUVrrZlrr750d6loTJkzg8uXLDB482JWLFUKIMg0YMIDAwEA+/PBDYxnc48zJ67jrrrsYOHAgaWlprFu3znQcIYQPW7JkCRkZGQwbNoy6desay+EeF5m6gV9//ZXIyEhq167NDz/8QKVKlRz22kIIYYuLFy/StGlTtNbs2LHD4QMmPO8iUzdQpUoVxo4dy5YtW5g1a5bpOEIIHzRx4kT27NlDQkKC8VFuHrHHDaC1pn379uzevZu9e/dy6623OvT1hRCiPMePHyciIoIHHniAxYsXO2UZXrfHDaCUIikpiZ9//pl33nnHdBwhhA956623OH/+vNsMTfaYwg3QvHlz4uLiSE5OJjc313QcIYQPyMzMZNasWfTv359GjRqZjgN4UKvkqlOnTtGwYUNatmzJ8uXLUUo5ZTlCCGG1WmnXrh0HDx5kz549Tj213StbJVfVqFGD9957j++//55//etfpuMIIbxYWloaGzdu5MMPP3T59Uiux+P2uAEuX75MTEwM586dIzc3l+DgYKctSwjhm65ej+SOO+5g48aNTj+13av3uAH8/f1JSkoiPz+f8ePHm44jhPBCV69HMnHiRCPXI7ke90pjh4cffpgePXowatQoDh8+bDqOEMKL7NmzhwkTJvD3v//d2PVIrsdjCzfAuHHj0FozaNAg01GEEF5Ca018fDyVK1dm9OjRpuOUyaMLt8Vi4c0332TevHmsWbPGdBwhhBfIyMhgyZIljBgxgjp16piOUyaPPDh5rcLCQiIjI6levTpZWVlyHRMhxE07f/480dHRBAUFkZ2dTUBAgMuW7fUHJ68VEhLC+PHjyc7OZtq0aabjCCE82Pjx4zlw4ABJSUkuLdr28vg9bijuST388MNs376dvXv3Ur16dZctWwjhHQ4dOkRkZCR/+tOfWLhwocuX71N73FB8HZOPP/6Y//u//2P48OGm4wghPNDrr7+O1tojhhh7ReEGuOeee3jppZeYPHkyO3bsMB1HCOFBVqxYweeff87gwYOxWCym49yQV7RKrvrll19o2LAhTZs2ZeXKlXIdEyHEDV26dInmzZsbPxPb51olV1WvXp2RI0eyevVqFixYYDqOEMIDpKSkkJOTQ0JCgsdcPsOr9rgBrly5wr333suJEyfYvXs3VapUMZJDCOH+Tpw4QUREBG3atOHbb781+indZ/e4ASpVqkRycjJHjx5l5MiRpuMIIdzYkCFDOHfuHB9//LFHtVa9rnAD3H///fTu3ZsjY8ZwqX598PMDiwXS001HE0KYlp4OFgvaz49hM2cy+5FHiIyMNJ3KLjYVbqXUa0qpHKXUTqXUZ0qpys4OVlGJrVsz1Wol4Ngx0BoKCiAuToq3EL4sPb24DhQUoLTGAvxtxQqPqws37HErpeoDa4EorXWRUmoBsFhr/Ul532Oyx/0bi6W4WJcWHg75+a5OI4RwB25cF5zR4/YHgpVS/kAIcOxmw7nMoUP2TRdCeD8vqQs3LNxa66PAOOAQ8CNwWmu9zNnBKiwszL7pQgjv5yV14YaFWyl1G9AVuBOoB4QqpZ4uY744pdRmpdTmkydPOj6pvUaOhJCQ/5lkDQ4uni6E8EmHX3qJc6UnhoR4XF2wpVXSETiotT6ptb4EfAncX3omrXWq1rqV1rpVzZo1HZ3Tfr16QWoqhIejleKQUiQ0blw8XQjhc7TW/P3bb3ktNJQrDRqAUsW97dRUj6sLthTuQ8B9SqkQVTzQsQOwy7mxHKRXL8jPR1mtfPbhh7y+ZQtLliwxnUoIYcDnn3/OqlWraDFuHJUOHwartfiApIcVbbDxzEml1LvAX4HLwFagn9b6Qnnzu8WoklIuXLjAPffcg9aaHTt2EBQUZDqSEMJFzp07R2RkJDVr1iQzM9Mtb7ji8FElWusRWutIrXUTrfUz1yva7iooKIikpCT27t1LQkKC6ThCCBf68MMPOXLkCBMnTnTLom0vrzxzsjyPPPII3bp14/333+fIkSOm4wghXGDfvn2MHTuWZ555hnbt2pmO4xA+VbgBEhISsFqtvP7666ajCCFc4LXXXiMoKMht79h+M3yucFssFoYMGcL8+fNZuXKl6ThCCCdavHgxGRkZDB8+nLp165qO4zBed1lXWxQVFREdHU1ISAhbt25165uCCiFuzvnz52nSpAn+/v5s376dwMBA05Guy6cv62qL4OBgEhMTycnJYdKkSabjCCGcYNy4cezfv5/k5GS3L9r28sk9bigejP/YY4+xbt068vLyqFOnjulIQggHyc/Pp3HjxnTp0oX58+ebjmMT2eO2wdU7w58/f57BgwebjiOEcKD4+HgqVarkEXdsvxk+W7gBGjZsyMCBA5k9ezbr1683HUcI4QCLFy/m3//+N8OHD6dBgwam4ziFz7ZKrvKEM6qEELa5ekAyICCA7Oxsj+ptS6vEDqGhoYwfP56tW7eSmppqOo4QogLGjh3rtQckr+Xze9xQfKCyY8eObN26lT179lCjRg3TkYQQdjp48CBRUVF07dqVefPmmY5jN9njtpNSiokTJ3L27Fnefvtt03GEEDfh6gHJcePGmY7idFK4S0RFRfHqq68ybdo0POnTghACFi1axNdff82IESO89oDktaRVco0zZ84QERGBxWJh/fr1+PnJ+5oQ7u78+fNER0cTFBTEtm3bPLa3La2Sm1StWjXGjh3Lpk2bmDlzpuk4QggbjBkzhgMHDnj9AclryR53KVprYmNj2blzJ3l5eXKgUgg3duDAAaKjo+nWrRufffaZ6TgVInvcFaCUYtKkSZw+fZohQ4aYjiOEuI74+Hj8/f194oDktaRwl6FJkybEx8czffp0Nm7caDqOEKIMGRkZfPPNN4wYMYL69eubjuNS0iopx9mzZ4mMjKR27dpyRqUQbubqpZmDg4PZtm2bV1yaWVolDlC1alUSEhLYunUrkydPNh1HCHGNMWPGcPDgQZKTk72iaNtL9rivQ2tNp06dyMzMJC8vj9q1a5uOJITPO3DgAFFRUXTv3p25c+eajuMwssftIEopkpOTKSws5I033jAdRwgB9O/fn4CAAJ87IHktKdw30KhRI9544w3mzJnDmjVrTMcRwqd98803ZGRk8M4771CvXj3TcYy5YatEKdUIuPYWEncBw7XWieV9j7e0Sq4qLCwkKiqKKlWqyD0qhTCkqKiIqKgoQkJCvOaA5LUc2irRWudprWO01jFAS6AQ+KqCGT1KSEgIH3/8MTk5OSQlJZmOI4RPGj16NPn5+UyaNMnrira97G2VdAD2a60LnBHGnXXp0oXHH3+cnW+9xeUGDcDPDywWSE83HU0I75WeDhYL2s+PPu++S3LbtsTGxppOZZy9hftJwLPPK71JSilmdOhA8sWL+B89ClpDQQHExUnxFsIZ0tOLt6+CApTWhAMvbdsm2xt2DAdUSgUCx4BorfVPZTwfB8QBhIWFtSwo8MKdcouluFiXFh4O+fmuTiOEd/Ox7c1ZwwH/BGwpq2gDaK1TtdattNatatasacfLepBDh+ybLoS4ebK9lcuewv03fLRN8puwMPumCyFunmxv5bKpcCulQoA/Al86N46bGzkSQkL+Z5I1OLh4uhDCoQ698ALnSk8MCZHtDRsLt9a6UGt9u9b6tLMDubVevSA1FcLD0UpxSCkSGjcuni6EcBir1crfMjIYUKUKVxo0AKWKe9upqbK9IWdO2q9XL8jPR1mtzPvoI17fsoWMjAzTqYTwKrNmzWL9+vW0nTiRSocPg9VafEBSijYgF5mqkIsXLxITE0NRURE5OTmElGqjCCHsd+rUKRo1akR0dDSrV69GKWU6kkvIRaZcJDAwkJSUFPLz8xk1apTpOEJ4hTfffJMzZ84wefJknyna9pLCXUGxsbE8/fTTjBkzhl27dpmOI4RHW7duHTNnzmTAgAFER0ebjuO2pFXiACdOnCAyMpJ77rmHlStXyl6CEDfh0qVLtGjRgtOnT5Obm0uVKlVMR3IpaZW4WK1atRg9ejSrV69m9uzZpuMI4ZGSkpLYuXMnSUlJPle07SV73A5itVp58MEHycvLIy8vj9tvv910JCE8xuHDh2ncuDEPPfQQX3/9tU9+apU9bgP8/PyYMmUKp0+fZtCgQabjCOFR4uPjsVqtTJw40SeLtr2kcDtQ06ZNGTBgADNnzuQ///mP6ThCeITFixfz5ZdfMmzYMCwWi+k4HkFaJQ527tw5oqOjCQ0NZevWrQQGBpqOJITbKiwspEmTJgQFBZGdne3T24u0SgwKDQ0lOTmZ3NxcJkyYYDqOEG5t1KhRHDx4kJSUFJ8u2vaSwu0Ejz/+OE888QTvvfceBw8eNB1HCLe0e/duxowZw9NPP81DDz1kOo5HkcLtJElJSVSqVIl//OMfOKMdJYQn01rz4osvEhoayrhx40zH8ThSuJ2kQYMGvP/++yxZsoSFCxeajiOEW/nkk09YvXo1Y8aMoXbt2qbjeBw5OOlEly9f5t577+XEiRPs2rWLatWqmY4khHEnT54kMjKSxo0bs2bNGvz8ZP8R5OCk2/D392fq1Kn8+OOPDBs2zHQcIdzCwIEDOXv2LKmpqVK0b5KsNSdr3bo1L730EsnJyWRlZZmOI4RR33//PXPmzGHQoEFERUWZjuOxpFXiAqdPnyYyMpL69euzadMmKlWqZDqSEC53/vx5mjZtCsD27dsJDg42nMi9SKvEzdxyyy0kJiaSlZVFSkqK6ThCGDFy5Ej27dvHlClTpGhXkOxxu4jWmkcffZQNGzawe/du6tWrZzqSEC6Tm5tLTEwMTz75JJ9++qnpOG5J9rjdkFKKlJQULl26RHx8vOk4QriM1WrlhRdeoGrVqowfP950HK8ghduF7r77boYOHcrnn3/OokWLTMcRwiVmzpzJ2rVrGTduHDVr1jQdxytIq8TFLl68SIsWLThz5gw5OTlUrVrVdCQhnOann34iMjKSZs2ayd2hbsDhrRKl1K1KqS+UUruVUruUUm0rFtF3BQYGMm3aNI4cOcLQoUNNxxHCqQYMGEBhYSFTpkyRou1AtrZKPga+1VpHAs0AuStuBbRt25aXX36ZiRMnsmnTJtNxhHCKZcuWMXfuXIYMGUJkZKTpOF7lhq0SpVQ1IBu4S9vYV5FWyY2dOXOGqKgoqlevTlZWFgEBAaYjCeEwhYWFNG3aFH9/f7Zv305QUJDpSG7P0a2Su4CTwCyl1Fal1HSlVGiFEgqqVatGSkoKO3bsYOzYsabjCOFQH3zwAQcOHGDq1KlStJ3AlsLtD7QAJmutmwPngMGlZ1JKxSmlNiulNp88edLBMb1Tly5d+Mtf/sJ7773H3r17TccRwiF27tzJ2LFj6dOnD7GxsabjeCVbWiV1gI1aa0vJ1w8Cg7XWj5X3PdIqsd2PP/5I48aNadGiBd9//70cwBEe7cqVK7Rr1479+/eza9cuatSoYTqSx3Boq0RrfRw4rJRqVDKpA5BbgXziGnXr1mXMmDGsXLmSWbNmmY4jRIUkJyezadMmkpKSpGg7kU3juJVSMcB0IBA4APTRWv+3vPllj9s+VquV2NhYdu7cya5du+TC8sIj5efnEx0dTWxsLBkZGfLp0U4OH8ettd6mtW6ltb5Ha93tekVb2M/Pz4/U1FTOnTtH//79TccRwm5aa1544QX8/PyYPHmyFG0nk1Pe3URkZCRDhw5l/vz5cjq88Dhz5sxh2bJlfPTRR4SFhZmO4/XklHc3IqfDC0904sQJGjduTGRkJP/5z3/krjY3Sa4O6KHkdHjhifr378+vv/7K9OnTpWi7iKxlNyOnwwtPkpGRwbx58xg6dCiNGzc2HcdnSKvEDV17OvzmzZsJDAw0HUmI3zlz5gzR0dHceuutZGVlyd9pBUmrxMNdezr86NGjTccRokyDBw/m6NGjzJgxQ4q2i0nhdlNdunThySef5P3332fnzp2m4wjxP9auXcvkyZPp378/rVu3Nh3H50irxI2dPHmSqKgo7rzzTtavX4+/v7/pSEJw/vx5YmJiuHDhAjt37iQ0VK455wjSKvESNWvWJDk5mczMTBISEkzHEQIovvJfXl4eU6dOlaJtiBRuN9ezZ0+6devGsGHDyMvLMx1H+Ljt27czevRoevfuTadOnUzH8VlSuN3c1bvDBwcH89xzz2G1Wk1HEj7q0qVL9O3bl9tuu03u1m6YFG4PULduXRITE1m3bh2TJk0yHUf4qLFjx5KVlcXkyZO5/fbbTcfxaXJw0kNorencuTNr1qxh586d3HnnnaYjCR+Sk5NDixYt6NatG/PnzzcdxyvJwUkvpJRi6tSpVKpUiX79+uGMN1whynL58mX69OlDtWrVSE5ONh1HIIXbo4SFhTF27FhWrFjB9OnTTccRPmL8+PFkZmaSkpJCzZo1TccRSKvE41itVjp27MjmzZvJycnhjjvuMB1JeLHc3FyaN2/On//8Z7744gvTcbyatEq8mJ+fH9OmTePKlSu8+OKL0jIRTnO1RVK1alU5KO5mpHB7oLvvvptRo0axePFi0tLSTMcRXiohIYEffviB5ORkuZ2em5FWiYe6cuUK7du3Jzc3l5ycHOrVq2c6kvAiu3fvJiYmhs6dO7Nw4UK5FZkLSKvEB1SqVIlZs2Zx4cIFnn/+eWmZCIe5cuUKffr0ITQ0lJSUFCnabkgKtweLiIjgo48+YvHixcycOdN0HOElEhMT2bhxIxMnTqROnTqm44gySKvEw1mtVjp06EBWVhY7duwgPDzcdCThwfLy8oiJieGRRx7hq6++kr1tF5JWiQ/x8/Nj1qxZaK3p27evXMtE3LTLly/Tu3dvgoODmTx5shRtN2ZT4VZK5SuldiiltimlZFfazVgsFiZMmMCKFStISUkxHUd4qNGjR7Np0yZSUlKoW7eu6TjiOmxqlSil8oFWWutTtryotEpc79prmWzbto2GDRuajiQ8yNatW2ndujU9evRg3rx5puP4JGmV+CClFNOnTycwMJBnn32WK1eumI4kPMT58+d55plnqFmzppxo4yFsLdwaWKaUylJKxTkzkLh59evXZ+LEiaxfv17umCNsNmzYMHJycpgxY4ZcrtVD2Noqqae1PqaUqgUsB17RWq8pNU8cEAcQFhbWsqCgwBl5xQ1orenRoweLFi1iy5YtREdHm44k3NiaNWuIjY0lLi6OKVOmmI7j0+xpldg9HFAp9Q7wq9Z6XHnzSI/brBMnThAdHU14eDgbNmwgICDAdCThhs6ePUuzZs1QSpGdnU2VKlVMR/JpDu1xK6VClVJVrz4GOgE7KxZROFOtWrWYOnUqWVlZvPvuu6bjCDc1cOBA8vPz+fTTT6Voexhbety1gbVKqWzgB2CR1vpb58YSFdW9e3f69OnDhx9+yNq1a03HEW5m0aJFTJs2jUGDBtGuXTvTcYSd5MxJL3b27FliYmKwWq1kZ2dTrVo105GEGzh16hRNmjShVq1aZGZmEhQUZDqSQIYDihJVq1YlLS2NQ4cO8corr5iOI9yA1poXX3yRX375hTlz5kjR9lBSuL1c27ZtGTp0KJ9++ikLFiwwHUcYNnPmTBYuXMgHH3xAs2bNTMcRN0laJT7g0qVLPPDAA+zdu5ft27fToEED05GEAXv27KF58+bcd999LF++HD8/2W9zJ9IqEf8jICCAtLQ0Ll68SO/eveVCVD7o4sWLPPXUUwQFBfHpp59K0fZw8tvzEQ0bNiQxMZEVK1aQmJhoOo5wsREjRpCVlcX06dOpX7++6TiigqRw+5DnnnuObt26MWTIELKzs03HES6ycuVKRo8eTb9+/ejevbvpOMIBpMftY06dOkXTpk2pXr06mZmZhISEmI4knOiXX37hnnvuITQ0lC1bthAaGmo6kiiH9LhFuWrUqMHs2bPJzc0lPj7edBzhRFpr4uLiOHHiBHPnzpWi7UWkcPugTp068eabbzJt2jTmz59vOo5wkmuH/rVs2dJ0HOFA0irxUZcuXeIPf/gDOTk5bN26lbvuust0JOFAubm53HvvvTL0z4NIq0TcUEBAAHPnzsXPz4+//vWvXLx40XQk4SCFhYX07NmTKlWqkJaWJkXbC8lv1IdZLBZmzJjB5s2bGTJkiOk4wkFeeeUVcnNzSUtLk3tHeikp3D6ue/fuvPzyy0yYMIFFixaZjiMqKC0tjZkzZ/LWW2/xxz/+0XQc4STS4xacP3+eNm3acPToUbKzs+UEDQ+1e/duWrVqRYsWLVixYgX+/v6mIwk7SI9b2KVy5crMnz+foqIievXqJTca9kBFRUX07NmT4OBgPvvsMynaXk4KtwAgMjKSlJQUVq9ezYgRI0zHEXZ67bXX2LFjB3PmzJFPTD5ACrf4Te/evenbty8jR44kIyPDdBxho3nz5jF16lTefPNNHn30UdNxhAtIj1v8j6KiItq1a8fBgwfJysqS8d1uLicnhzZt2hATE8PKlSvlxtAeTHrc4qYFBwfzxRdfANCjRw+KiooMJxLlOX36NN27d6dq1aosWLBAirYPkcItfueuu+4iLS2Nbdu28fLLL+OMT2WiYqxWK88++yz79+9nwYIF1KtXz3Qk4UJSuEWZHnvsMYYNG8Ynn3zC9OnTTccRpYwePZp//etfjBs3jgcffNB0HOFi0uMW5bpy5QqdO3dm1apVrFu3jlatbGq/CSdbvnw5jz76KD179mTu3LkopUxHEg5gT49bCre4rlOnTtGyZUuUUmRmZlKzZk3TkXxaQUEBLVu2pG7dumzcuFEu1epFnHJwUilVSSm1VSkl48R8SI0aNVi4cCHHjx9n8oMPosPDwc8PLBZITzcdzzekp4PFgvbzIzAigq7nzvHll19K0fZh9vS4+wO7nBVEuK9WrVrxXd++DMzLQx06BFpDQQHExUnxdrb09OL1XFCA0pq6Fy8yVWsa/vCD6WTCIJsKt1KqAfAYIEepfNQDixfzu/27wkJ4+20TcXzH228Xr+dr+F+4IOvdx9m6x50IDAKsTswi3NmhQ/ZNF44h612U4YaFWyn1OHBCa511g/nilFKblVKbT5486bCAwk2Ehdk3XTjExTp1yn5C1rtPs2WPux3QRSmVD8wDHlZKpZWeSWudqrVupbVuJSMPvNDIkVDqjvBFSnFOPrI7zU8//cTACxcoLD3cLySk+PchfNYNC7fWeojWuoHW2gI8CazQWj/t9GTCvfTqBampEB4OSlFUqxZxSvHnzz6T2545QWFhId26dWNGURE/vf/+b+ud8PDi30OvXqYjCoPkzElhu169ID8frFaCf/qJTp98wsqVK+nXr5+cFu9AV65c4amnnmLTpk2kp6dz59tv/7beyc+Xoi2w62rrWutVwCqnJBEe55lnniE/P5/hw4dz55138u6775qO5PG01sTHx/Pvf/+bpKQknnjiCdORhBuS22SIChk6dCj5+fm89957WCwW+vTpYzqSRxs/fjzJyckMHDiQV155xXQc4aakcIsKUUoxZcoUDh8+TFxcHPXr16dTp06mY3mkefPm8cYbb9CzZ0/GjBljOo5wY9LjFhUWEBDAF198QVRUFE888QQbNmwwHcnjfPfdd/Tu3Zv27dsze/Zs/Pxk0xTlk78O4RDVqlVj6dKl1KtXj86dO7N9+3bTkTzGhg0b6Nq1K40aNeKrr76icuXKpiMJNyeFWzhMnTp1WL58OaGhoXTq1Il9+/aZjuT2tm3bRufOnalXrx7Lli2jevXqpiMJDyCFWziUxWJh+fLlXL58mY4dO3LkyBHTkdxWXl4enTp1omrVqnz33XfUKe8sSSFKkcItHK4g6/1lAAAJ7ElEQVRx48YsXbqUX375hQ4dOnDs2DHTkdzOwYMH6dixI0opvvvuO8LDw01HEh5ECrdwipYtW7JkyRKOHTtGbGwsR48eNR3Jbezbt4/27dtz7tw5li1bRkREhOlIwsNI4RZO065dO5YuXcrx48eJjY2VtgnF7ZE//OEPFBUVsXLlSpo1a2Y6kvBAUriFU91///0sXbqUn376idjYWA758OVIc3NziY2N5dKlS1K0RYVI4RZO17ZtW5YtW8bJkydp164du3b53o2Utm3bxkMPPYTWmlWrVtG0aVPTkYQHk8ItXOK+++5j1apVXLp0iQceeICNGzeajuQy33//Pe3btycoKIjVq1cTFRVlOpLwcFK4hcs0b96cdevWcdttt9GhQweWLFliOpLTzZs3jz/96U+Eh4ezfv16GjVqZDqS8AJSuIVL3X333axbt46IiAi6dOnCjBkzTEdyCq01CQkJ/O1vf+O+++5jzZo1NGjQwHQs4SWkcAuXq127NqtXr+bhhx+mX79+vPbaa1y+fNl0LIe5cOEC/fr1Y8CAAfTo0YNly5Zx2223mY4lvIgUbmFEtWrVWLRoEf379ycxMZFxLVpw5Y47wM8PLBZITzcd0T7p6WCxoP38+KVaNYpmzmT48OEsWLBArj0iHE4u6yqM8ff3JzExkS6//kqbGTOodPWJggKIiyt+7Al3e0lPL85bWIgC6l68yKeBgfhHRBS/EQnhYMoZt5xq1aqV3rx5s8NfV3gpi6W4WJcWHl58qy43p8PDUWWNT/eQ/MI9KKWytNatbJlXdgeEeeWclKM94GSdQ4cOlZ/TA/ILzySFW5gXFlbm5MPArFmz3PJGxFarlUmTJhEdHc0RpcqeqZyfS4iKksItzBs5EkJC/meStXJlPomIoG/fvjzyyCPk5eUZCvd727dvp3379vzzn/+kbdu2VJ4w4Xf5CQkp/rmEcAIp3MK8Xr0gNbW4J6wUhIfjN306Q3NzmTRpEhs3bqRJkyb079+fn3/+2VjM48eP8/zzz9O8eXN27drFJ598wtKlS6kVH/+7/KSmesaBVeGZtNYO/9eyZUsthKMcP35cv/DCC9rPz0/fcssteujQofrkyZMuW/6xY8f066+/rkNDQ3VAQIB+7bXX9M8//+yy5QvfAGzWNtbYG88AlYEfgGwgB3j3Rt8jhVs4w44dO3T37t01oENCQvQ///lPnZ2d7ZRlWa1WnZmZqZ9//nkdFBSk/fz89FNPPaX37NnjlOUJ4ejCrYAqJY8DgE3Afdf7HincwplycnL0M888owMDAzWgW7durRMTE/X+/fvL/6a0NK3Dw7VWqvj/tLQyZ8vLy9MfffSRjomJ+e0NIi4uTu/bt88pP4sQV9lTuO0ax62UCgHWAi9prTeVN5+M4xau8PPPP5OWlsaMGTPYsWMHUHzbtHbt2tGmTRuaNWuGxWKhxrJlqJITZK7SISGcGTeOnGbNyMnJYf369axdu/a3Gxy3atWKvn378tRTT3HLLbcY+fmEb7FnHLdNhVspVQnIAv4/YJLW+s3rzS+FW7ja/v37+eabb/j222/54Ycf+O9///vbcwVKEVbG33k+cGfJ4xo1anD//ffTsWNHunbtSpgM5RMu5vDCfc0L3wp8Bbyitd5Z6rk4IA4gLCysZUFZZ8IJ4QJaa/bu3Utubi4FBQW8Gh9PWSOtNbBk0SIiIiK4++67UeWNxxbCBZxWuEtefARwTms9rrx5ZI9buBUPP6Ve+AaHnvKulKpZsqeNUioY6AjsrlhEIVyojBN85AQZ4clsOQGnLrBSKbUdyASWa60znBtLCAcq4wQfOUFGeDK5OqAQQrgBuTqgEEJ4MSncQgjhYaRwCyGEh5HCLYQQHkYKtxBCeBinjCpRSp0EbvbUyRrAKQfGcRTJZR/JZR/JZR9vzBWuta5py4xOKdwVoZTabOuQGFeSXPaRXPaRXPbx9VzSKhFCCA8jhVsIITyMOxbuVNMByiG57CO57CO57OPTudyuxy2EEOL63HGPWwghxHUYK9xKqUeVUnlKqX1KqcFlPB+klJpf8vwmpZTFTXI9q5Q6qZTaVvKvnwsyzVRKnVBK7SzneaWUSirJvF0p1cLZmWzMFauUOn3Nuhruolx3KKVWKqV2KaVylFL9y5jH5evMxlwuX2dKqcpKqR+UUtklud4tYx6Xb4825nL59njNsisppbYqpX53tVSnry9bb07pyH9AJWA/cBcQSPEd5KNKzfMyMKXk8ZPAfDfJ9SyQ7OL11R5oAews5/nOwBKKb+x8H7DJTXLFAhkG/r7qAi1KHlcF9pTxe3T5OrMxl8vXGTbcENzQ9mhLLpdvj9csewAwt6zfl7PXl6k97tbAPq31Aa31RWAe0LXUPF2B2SWPvwA6KOffW8qWXC6ntV4D/HKdWboCn+piG4FblVJ13SCXEVrrH7XWW0oenwV2AfVLzebydWZjLpcrWQe/lnwZUPKv9MEvl2+PNuYyQinVAHgMmF7OLE5dX6YKd33g8DVfH+H3f8C/zaO1vgycBm53g1wAPUo+Xn+hlLrDyZlsYWtuE9qWfNRdopSKdvXCSz6iNqd4b+1aRtfZdXKBgXVW8rF/G3CC4pullLu+XLg92pILzGyPicAgwFrO805dX6YKd3n3brV3HkezZZnfABat9T3Ad/z/76ommVhXtthC8Wm8zYCJwL9cuXClVBVgIRCvtT5T+ukyvsUl6+wGuYysM631Fa11DNAAaK2UalJqFiPry4ZcLt8elVKPAye01lnXm62MaQ5bX6YK9xHg2nfGBsCx8uZRSvkDt+D8j+U3zKW1/llrfaHky2lASydnsoUt69PltNZnrn7U1VovBgKUUjVcsWylVADFxTFda/1lGbMYWWc3ymVynZUs8/+AVcCjpZ4ysT3eMJeh7bEd0EUplU9xO/VhpVRaqXmcur5MFe5MoKFS6k6lVCDFzfuvS83zNdC75PFfgBW6pNNvMlepPmgXivuUpn0N/L1kpMR9wGmt9Y+mQyml6lzt6ymlWlP89/azC5argBnALq31hHJmc/k6syWXiXWmbLshuMu3R1tymdgetdZDtNYNtNYWimvECq3106Vmc+r68nfUC9lDa31ZKfVPYCnFIzlmaq1zlFLvAZu11l9T/Ac+Rym1j+J3qifdJNerSqkuwOWSXM86O5dS6jOKRxvUUEodAUZQfKAGrfUUYDHFoyT2AYVAH2dnsjHXX4CXlFKXgSLgSRe8+ULxHtEzwI6S/ijAW0DYNdlMrDNbcplYZ3WB2UqpShS/USzQWmeY3h5tzOXy7bE8rlxfcuakEEJ4GDlzUgghPIwUbiGE8DBSuIUQwsNI4RZCCA8jhVsIITyMFG4hhPAwUriFEMLDSOEWQggP8/8A94RBmT6d2x4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot((z_fine-z_fine.min())*1e5, w_fit*1e6, '-k')\n",
    "plt.plot((z-z_fine.min())*1e5, caustic*1e6, 'or')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demo 3: Normal propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, lets have a quik look on how to propagage by examining the field at the focus and printing the power to the console."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(0m) = 10.00 kW\n",
      "r(0m) = 10.00 cm\n",
      "P(100m) = 9.05 kW\n",
      "r(100m) = 23.33 cm\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARAAAAElCAYAAAAsvb0LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH/9JREFUeJzt3Xm0ZWV95vHvc++5A0WhFIVDMUOEOKRdGJEMtlEjEWzTQLdEMTEpbAwtabtXL4clSjqdEOnEDM3Sxo6iUcohgmDbllMTBnGlDSBlFhEBoQpQqa5i0KLQqjvVvffXf+x3H9+7a597T+177lT3+ax11tlnT+c9p+556t3vu/e7FRGYmTXRt9QFMLOVywFiZo05QMysMQeImTXmADGzxhwgZtaYA8R6TtKHJf2XRXqvl0raKmmPpHMX4z3tZxwgK5yk70s6o8t1b5X0loUuU0S8NSL+NL3nKyRtX8C3uwy4MiLWRsT/XqzPaAUHiK10xwP39Gpnkvp7ta/VwAFyEJF0gaT/K+mvJD0p6WFJr0nLLgdeBlyZqvtXpvnPlXSjpF2S7pf0+mx/V0v6kKSvSPqppDsk/VxaJklXSHpc0lOSviPpF7Lt3ifpUOBrwFHpPfdIOkrSiKT12fu8WNITkgZqPtPpkm6TtFvSTklXShpMyx4ETgK+lPb9Zw0/499I+qqkvcAre/qPcrCLCD9W8AP4PnBGmr4A2Af8PtAPXAzsAJSW3wq8Jdv2UOAR4M1AC/hF4EfAC9Lyq4FdwOlp+WeAa9KyM4FvA4cDAp4HbMi2e1+afgWwvVLmrwIXZ6+vAP5Hh8/3YuCX0/ufANwH/Oe6zz+Pz/gU8FKK/1CHl/rfdCU9XAM5+PwgIj4aEVPAJmAD8KwO6/4m8P2I+ERETEbEPwGfB87L1vlfEfGtiJikCJBT0/x9wGHAcykC6r6I2NllGTcBb4L2IcMbgU/VrRgR346I21P5vg98BHh5l+/T7Wf8YkR8MyKmI2LsAPa96rWWugDWc4+WExExIglgbYd1jwd+SdLubF6LmT/mR7PpkXJfEXFLOkT4EHCcpC8A74yIn3RRxi8CH5Z0EnAK8FREfKtuRUmnAP8dOA1Yk8r37S7eo9TNZ3zkAPZnGddAVpfqpdePAN+IiMOzx9qIuLirnUV8MCJeDLyAIgje1cV7kv6X/xzwO8Dv0qH2kfwN8D3g5Ih4GvBeikOmjsWqvO7mM/qS9IYcIKvLYxSNjqUvA6dI+l1JA+nxEknPm2tHab1fSg2fe4ExYKrDe66X9PTK/E9StNmcDXx6lrc6DPgJsEfScynadWbTs89oc3OArC4fAM5LPTQfjIifAq8GzqdobH0UeD8w1MW+ngZ8FHgS+AHwY+CvqitFxPeAzwIPpZ6Uo9L8bwLTwD+lto1O3gn8NvDT9H7XLuJntDmUrfNmi07SLcDfRcTHlros1owDxJaEpJcANwLHplqCrUA+hLFFJ2kTcBPF+RwOjxXMNRAza8w1EDNrzAFiZo05QKxrC3Fp/mKOHWK95wBZwdIVse9KA+qMSvqhpP9WXq263JRXC+fzIhs75AD31fU4KLZwHCAr2weBi4Dfozhj8zXAGcA1i10QSb6uahVygKxQkk4G/gD4nYi4LV1peg/wOuC1kl6e1psxQle1FiDpA5IekfQTSd+W9LJs2SFpvIwnJd0LvKRShu9Lerek7wB7JbUkXSLpwTR+yL2S/k1a93nAh4FfSWN17E7zr5b0vmyf50i6K5XnQUlnNfhufjPtY7ekf5T0wkqZ35nGL3lK0rWShtOyIyV9OW23S9I/SPJvZBb+clauV1GMszHjKtaIeAS4neL07W7cSXGJ/hHA3wHXlT8o4L8CP5ceZwIba7Z/I/Ba4PB0yf+DFIP6PB34E+DTkjZExH3AW4Hb0sVsh1d3JOl0imtk3kUxzsivUYz30TVJvwh8HPj3wHqKy/83S8pPXX89cBZwIvBCimtyAN4BbAeeQTEEwnvxhXazcoCsXEcCncbf2EnxI5hTRHw6In6cajB/TXGNyM+nxa8HLo+IXSmYPliziw9GxCMRMZr2d11E7Ehja1wLbKUYkKgbFwIfj4gb0/b/L11LcyB+H/hIRNwREVMRsQkYpxiUKC/zjojYBXyJmWOcbACOj4h9EfEP4ROlZuUAWbl+RPHHXmcD8EQ3O5H0Dkn3per8boqaw5Fp8VHMHCvjBzW7mDGWhqTfyw4fdgO/kO1vLsdS1GDm43jgHeX7pzIcS/FZSrVjnAB/CWwD/l7SQ5IumWdZDnoOkJXrFuDYVO1vk3Qsxf+230iz9lIMxFN6drbuy4B3U9Q01qXDiqf42XgbOyl+fKXjasrR/h9a0vEUV8y+DVif9vfdbH9z/W/+CMXh0nw8QlFrysf/WBMRn51rw4j4aUS8IyJOAv418HZJr5pneQ5qDpAVKiIeoGiU/IykX5bUL+kFFMP1/SPFtSYAdwH/VtIaSc+hOEwoHQZMUtRWWpL+iOIy/dLngPdIWifpGOA/zlGsQylC4gkASW+mqIGUHgOOmaWb+W+BN0t6laQ+SUenMUA6GZA0nD1aFAH2VhVjlUjSoZJeK+mwOcpeNr4+R5IoxiCZon6ME0scICvb24CPUQzIM0Lxv/0PgHMjYjqtcwUwQfHj3UQxrmnpBopR0x9I240x85DkT9L8h4G/Z/aRw4iIe4G/Bm5L7/cvgG9mq9xCcQuGRyX9qGb7b1EMfnwFRU3oGxSHJJ18FRjNHn8cEVso2kGupBirZBs/aySdy8kUwbsnfYb/GRG3drntquSL6Q4iki4DzgV+LSJ2z7W+2Xw5QA4ykt4GbIuI/7PUZbGDnwPEzBpzG4iZNbYir1+Q5GqT2QKLiNlunwG4BmJm8+AAMbPGHCBm1pgDxMwac4CYWWMOEDNrzAFiZo05QMysMQeImTXmADGzxhwgZtaYA8TMGnOAmFljDhAza6wnASLpLEn3S9pWNxR+uhvaE2m4/7sqd0rbqOLerlsl1d24yMyWq4iY1wPop7iXx0nAIPDPwPMr61wAXFmz7RHAQ+l5XZpe18V7hh9++LGwj25+/72ogZxOMQbnQxExQXFj53O63PZM4MZ057MngRspbjloZitALwLkaGbeCmB7mlf1unRD4+vTzY8OZFskXSRpi6QtPSizmfVALwKkbtizqLz+EnBCRLyQ4r4bmw5g22JmxFURcVpEnNa4pGbWU70IkO3MvP3hMcCOfIV08+bx9PKjwIu73dbMlq9eBMidwMmSTky3LDwf2JyvICm/CfTZwH1p+gbg1enWieuAV6d5ZrYCzHtU9oiYTDczuoGiR+bjEXFPukvalojYDPwnSWdT3Id1F+lWgxGxS9KfUoQQwGURsWu+ZTKzxbEibyzl2zqYLTzf1sHMFpQDxMwac4CYWWMOEDNrzAFiZo05QMysMQeImTXmADGzxhwgZtaYA8TMGnOAmFljDhAza8wBYmaNOUDMrDEHiJk15gAxs8YcIGbWmAPEzBpzgJhZYw4QM2vMAWJmjTlAzKwxB4iZNeYAMbPGHCBm1pgDxMwac4CYWWMOEDNrzAFiZo05QMysMQeImTXmADGzxhwgZtaYA8TMGnOAmFljDhAza6wnASLpLEn3S9om6ZKa5W+XdK+k70i6WdLx2bIpSXelx+ZelMfMFociYn47kPqBB4DfALYDdwJvjIh7s3VeCdwRESOSLgZeERFvSMv2RMTaA3zP+RXazOYUEZprnV7UQE4HtkXEQxExAVwDnFMpyNcjYiS9vB04pgfva2ZLrBcBcjTwSPZ6e5rXyYXA17LXw5K2SLpd0rmdNpJ0UVpvy/yKa2a90urBPuqqObWHGJLeBJwGvDybfVxE7JB0EnCLpLsj4sH9dhhxFXBV2o8PYcyWgV7UQLYDx2avjwF2VFeSdAZwKXB2RIyX8yNiR3p+CLgVeFEPymRmi6AXAXIncLKkEyUNAucDM3pTJL0I+AhFeDyezV8naShNHwm8FLgXM1sR5n0IExGTkt4G3AD0Ax+PiHskXQZsiYjNwF8Ca4HrJAH8MCLOBp4HfETSNEWY/Xnee2Nmy9u8u3GXgttAzBbeYnXjmtkq5QAxs8YcIGbWmAPEzBpzgJhZYw4QM2vMAWJmjTlAzKwxB4iZNdaLq3HNZkiXK3RlJZ4JbT/jALF5O5DAmG1bh8nK4wCxxuYKjk7LOwWFw2TlcYDYAesUDN3WROrWqwZGuY6DZHlzgNgBqfvx5/OaHM5ERMfAkOQQWcYcIDanTgFRTs+1vKoaCOXruiDJ59Vta0vLAWKz6qbGMVuQzHa4Uj5XaxmzhYZrJMuLA8RqVX/41UDIn+tC5EACpAyMiNgvSPJ9dAocWzoOEJtTNSyq0319fR3XyZ+hPjzK6enp6fb6dYc5+XyHyPLgALH9zHYoUhciZYDU1UY61UCqAVIXKJKYnp6eER4OkeXFAWIzdAqPukOW8tHf3187v1OAlDWNPDzqAqXaFlIud4gsHw4Qa6trHK0++vr62jWOcrq/v3+/+eVzdb/VgChrGNPT00xOTrany/lTU1P7bTNXI6wtHgeIAd2FRx4a+aPVau03L6+ZVOUBUU6XtZLp6en2e+XrVLkmsjw4QGzWw5ZqO0de2yinBwYGZswr18sDJP+Bl7WK/DkPkOnpaaampmZsUw0St4ksDw6QVe5AwiMPjfy51WrRarXaoZHXTKrtIHlg5GExNTUF0J6WtF+w9PX1zeipKffn4Fg6DhCbYbYQyYMjD4yhoaH263xZGSBVZVvH1NQUk5OT7RApl/X19bVf57WO8vCmWl6Hx9JxgBjQuQ0kb9coAyOvdbRarXaAlG0h5fTAwEDHGkgZHmWAlKHS19fH5ORke5u69o88SKqHLw6UxeUAWcU6ddNWG0Tz4BgaGmJgYICBgQFarRaDg4OsWbOGgYEBBgcH91uWnyNS9qDs27evHSD79u1j3759TExMMDAw0J4ul8HMWkneNlLXwOoQWVwOkFWuruaRT+chUtYqyoAow+SQQw5hcHCw/ciX570y5Q9+37597dpHGRitVmtGD0ypDBGY2QWc1zjKZbb4HCCrVF3jaTndqd0jP3Qpw2N4eJhDDjmEoaGhGSEyODjI8PBwex+l6enpdg1jYmKiHR6tVmvGoUup1Wq1g6e/v39GgHQ6jKlO28JxgBjQ/WFMfqgyPDzcDpDh4WGGhobaj7JmUoZPaXp6mvHx8fahy/j4OOPj47RaLSYmJvZrL5mYmJhx3khfX9+MAHEtZGk5QKxtthPH8trHwMDAjLA45JBD2iFSBknZNlIGSH4IMz4+3q59lIc7/f39jI+PAzMPVcpDm/JRduVWT5V3zWNpOEBWodnOOu0UGIODgwwNDbFmzZr2c/lYv359e7oMk6GhIdauXdvumckDZHR0tF3zGB0dZXR0lJGREQCGhoYYGRlpb7dv377aruCymxfqe2XKz+YwWVgOkFVsrnaQukOYsnelfJQ1kDxQhoeHWbNmTTtAyjNSyy7c/v5+BgcHZ9ROAIaHh2d0705OTrbbRvIT1PKenbzMroUsPgeI1Z48Vj37ND9tvayR5AFShsihhx7aDpHDDjusY4CMjY219wtFzWF4eLjdM1M+VwOs7MatPqrtIrY4HCCrTN2ZnNXX1RpIOV2GQR4ieTtIWfMoD2U61UDya2qA9inrQ0ND7W7dffv2tdtH8h6aag2k7vCm+nkcKAvHAWLA/oP/5D/QvBs3bxspe2PKEClrHvkhTX6hXbU3BZhxLUy5v8HBwXbPTH6KfH9/fztI6srqoFh8Pbk3rqSzJN0vaZukS2qWD0m6Ni2/Q9IJ2bL3pPn3SzqzF+Wx5uoaVfPn/DqXagNree5H3r1bhkn+nPfYDA8Pz9g+r3XkV/RWy5GX1ZbOvGsgkvqBDwG/AWwH7pS0OSLuzVa7EHgyIp4j6Xzg/cAbJD0fOB94AXAUcJOkUyJiCltQnWob1dd1NY/qIUweFGU7yKGHHsratWvbjaX5IQyw37ypqal2TWZ8fLz9HuX75u0wZbnyK3Nna1S1hdOLGsjpwLaIeCgiJoBrgHMq65wDbErT1wOvUvGvfQ5wTUSMR8TDwLa0P1sis/XM5NPVHpHqNTN119FUl9et36mnpZty2eLrRYAcDTySvd6e5tWuExGTwFPA+i63tSXQ6Uda14VarlNXc6nrMckPRaqP6v46lcmWh140otb9q1brjp3W6WbbYgfSRcBFB1Y0a6ru/ixAuxG07rYL1bFOp6am2teyVNfJ91Md6zRft1OZbHnoRQ1kO3Bs9voYYEendSS1gKcDu7rcFoCIuCoiTouI03pQZuugU3BUp/PTy/PAKB/V0caqj07rV/dbHUB5tmlbfL0IkDuBkyWdKGmQolF0c2WdzcDGNH0ecEsU//KbgfNTL82JwMnAt3pQJjtAs/1YqzWH6nCE+Zmj+WX65Tkd5cVz5TUw5bLqNnVhMtuI7HXTtrjmfQgTEZOS3gbcAPQDH4+IeyRdBmyJiM3A3wKfkrSNouZxftr2HkmfA+4FJoH/4B6YxVG9ijWfzoOiOnpYeRHc+Pg4Y2NjtFot9u7dO+N6l6mpqfY4HnXngezZs4exsTFGRkbYu3cve/fuZc+ePYyMjDA6OsrY2NiMsKkLmE6HPXXTtnB6ciJZRHwV+Gpl3h9l02PAb3XY9nLg8l6Uw5rr1AZRBkKnMClrFmV3bX6GaX6lbR4uIyMjjI2NtcNiZGRkRmCU44XkwVG+f7ftJbY4fCbqKpPXNjq9rtZCyunyx1w9PCnPGi17XoD2AECzBUgZHmWYlPuqO8yp3kemLkg6fV5bOA4QmzFEYLWXJK95lMMRloP/tFqt9iX5+Rmi+dilnQJkYmKiHR51hy5liFQPpbptF7HF4QBZxbptB6kevuTDELZaLUZHR2dsmx9uHMh4INUAKQ9lZmtczT9L3bQtLAeItc3W41KOyTExMUFfXx8TExPt9o5yMKD8UKccirB6b5iIYHR0dMbhTxkiszWglocxnXpnbGk4QFaham0jV15fko/4la83MDDQDomyHQRgZGRkXmOijo+Ps3v37nbvztjYGBMTE+3nTr0xs4WJg2XhOUAM2H9M0epNrqemptpDDJbrlQYHB/cLlcHBwfb9XOYalb18lLWPsbGxGb0xcx3CuBaydBwgq9Rs7R/lcxki+b1q6269MDAwMKN7Nz9prNv7wuQBUm5f1407Wy+Max+LzwGyytWFR94jk4cIzLzRU7nN6OjojEAoz/8ohyTMB/yJ6HxnumrNo3rY4trH8uMAWcW66YWBmVfaTk5OztgeaJ9lmvfSlAMnl4cv+fgd1baMMijKWkc+Jmo318nkZalO28JygNh+P7gyTKo1kbIWUgZAeeZpeW5I2WCaj1pWdwl+eSPtMnDKfZcBkncZ14VIfp5KWd582haPA8RmyE8qy+flvTPlsvIHXc4vA6VsPC0DpNqNWw2EMjAmJiZqr9bttt3DFp8DZJWrO4yp/kjzO8Ll16OU96ote2ryQYTyAZBh5hCD5T7y57JxNQ+VTrWOujNQbWk4QKxjiOTyxtT8B5y3cfT19bW7e8tDnLpDmGoo5Cer1bV3zHX6ev45bHE5QAyY/eSyXHnIUv0hV8dHzXtuOl2sl9cqygDJQ6VTrcPhsXw4QKyjuh9l3rWb/5jLeXm3bX4CWV6r6RQiZe2mU2h0OuPUlo4DxNqq54Hk8+u6e8ugqLvdZL687n3qwqSuobTa25I/V6dt8TlAbD+d2kTyWkfdyWflc/VWDNVDmOpzp7CoO1xxeCwvDhCbU6eG1bwHplNodKqBdHquq51Ut6tO29JxgFitas0i/8HOFiZ1NZJO+64GQqdu2ble29LRSvzHkLTyCr2CdbrJVF1IdJrOdQqEuhqGw2PpRMScd/JygFhXZrtL3HzuIDfbYYnDY2k5QGxBdAqI+dx6stPf4Ur8+zxYdBMgbgOxA1Zt56jOL3V7CNN0HVt6DhBrrFOQVJcf6P5s5XCA2Lx1W/PoZltbWRwg1nMOhdWjFzfXNrNVygFiZo05QMysMQeImTXmADGzxhwgZtaYA8TMGnOAmFljDhAza8wBYmaNOUDMrLF5BYikIyTdKGlrel5Xs86pkm6TdI+k70h6Q7bsakkPS7orPU6dT3nMbHHNa0AhSX8B7IqIP5d0CbAuIt5dWecUICJiq6SjgG8Dz4uI3ZKuBr4cEdcf4Pv6ai2zBdbNgELzPYQ5B9iUpjcB59YU4oGI2JqmdwCPA8+Y5/ua2TIw3wB5VkTsBEjPz5xtZUmnA4PAg9nsy9OhzRWShmbZ9iJJWyRtmWeZzaxH5jyEkXQT8OyaRZcCmyLi8GzdJyNiv3aQtGwDcCuwMSJuz+Y9ShEqVwEPRsRlcxbahzBmC64nY6JGxBmdlkl6TNKGiNiZwuDxDus9DfgK8IdleKR970yT45I+AbxzrvKY2fIx30OYzcDGNL0R+GJ1BUmDwBeAT0bEdZVlG9KzKNpPvjvP8pjZIppvL8x64HPAccAPgd+KiF2STgPeGhFvkfQm4BPAPdmmF0TEXZJuoWhQFXBX2mZPF+/rQxizBeb7wphZY4vRjWtmq5gDxMwac4CYWWMOEDNrzAFiZo05QMysMQeImTXmADGzxhwgZtaYA8TMGnOAmFljDhAza8wBYmaNOUDMrDEHiJk15gAxs8YcIGbWmAPEzBpzgJhZYw4QM2vMAWJmjTlAzKwxB4iZNeYAMbPGHCBm1pgDxMwac4CYWWMOEDNrzAFiZo05QMysMQeImTXmADGzxhwgZtaYA8TMGnOAmFljDhAza2xeASLpCEk3Stqantd1WG9K0l3psTmbf6KkO9L210oanE95zGxxzbcGcglwc0ScDNycXtcZjYhT0+PsbP77gSvS9k8CF86zPGa2iBQRzTeW7gdeERE7JW0Abo2In69Zb09ErK3ME/AE8OyImJT0K8AfR8SZXbxv80KbWVciQnOtM98ayLMiYmd6s53AMzusNyxpi6TbJZ2b5q0HdkfEZHq9HTi60xtJuijtY8s8y2xmPdKaawVJNwHPrll06QG8z3ERsUPSScAtku4GflKzXseaRURcBVyVyuQaiNkyMGeARMQZnZZJekzShuwQ5vEO+9iRnh+SdCvwIuDzwOGSWqkWcgywo8FnMLMlMt9DmM3AxjS9EfhidQVJ6yQNpekjgZcC90bR+PJ14LzZtjezZSwiGj8o2jFuBram5yPS/NOAj6XpXwXuBv45PV+YbX8S8C1gG3AdMNTl+4YffvixsI9ufovz6oVZKm4DMVt4i9ELY2armAPEzBpzgJhZYw4QM2vMAWJmjTlAzKwxB4iZNeYAMbPGHCBm1pgDxMwac4CYWWMOEDNrzAFiZo3NOaDQMvUj4Ac92M+RaV8rzUost8u8eHpR7uO7WWlFXs7fK5K2RMRpS12OA7USy+0yL57FLLcPYcysMQeImTW22gPkqqUuQEMrsdwu8+JZtHKv6jYQM5uf1V4DMbN5cICYWWMHfYBIOkLSjZK2pud1HdabknRXemzO5p8o6Y60/bWSBpdDmSWdKuk2SfdI+o6kN2TLrpb0cPZ5Tl3g8p4l6X5J2yTtd4N1SUPpu9uWvssTsmXvSfPvlzTnfZEXscxvl3Rv+m5vlnR8tqz2b2UZlPkCSU9kZXtLtmxj+nvaKmljzwo1n/vCrIQH8BfAJWn6EuD9Hdbb02H+54Dz0/SHgYuXQ5mBU4CT0/RRwE7g8PT6auC8Rfp++4EHKe7xM0hx/5/nV9b5A+DDafp84No0/fy0/hBwYtpP/zIp8yuBNWn64rLMs/2tLIMyXwBcWbPtEcBD6Xldml7Xi3Id9DUQ4BxgU5reBJw7y7ozSBLw68D1TbafhznLHBEPRMTWNL2D4raiz1iEslWdDmyLiIciYgK4hqL8ufzzXA+8Kn235wDXRMR4RDxMcYOx05dDmSPi6xExkl7eTnHr1aXUzffcyZnAjRGxKyKeBG4EzupFoVZDgDwrInYCpOdndlhvWNIWSbdLKn+w64HdUdy7F2A7cPTCFhfovswASDqd4n+lB7PZl6fq9xXlrUUXyNHAI9nruu+ovU76Lp+i+G672XYhHOj7Xgh8LXtd97ey0Lot8+vSv/v1ko49wG0P2Eq9FmYGSTcBz65ZdOkB7Oa4iNgh6STgFkl3Az+pWa8n/d49KjPppuafAjZGxHSa/R7gUYpQuQp4N3BZ89LOXoSaedXvqNM63Wy7ELp+X0lvorhV68uz2fv9rUTEg3Xb91A3Zf4S8NmIGJf0Vopa3693uW0jB0WARMQZnZZJekzShojYmX5sj3fYx470/JCkW4EXAZ8HDpfUSv9zHgPsWC5llvQ04CvAH0bE7dm+d6bJcUmfAN7ZizJ3sB04Nntd9x2V62yX1AKeDuzqctuF0NX7SjqDItBfHhHj5fwOfysLHSBzljkifpy9/Cjw/mzbV1S2vbUXhVoNhzCbgbLVeSPwxeoKktaV1XxJRwIvBe6NogXq68B5s22/ALop8yDwBeCTEXFdZdmG9CyK9pPvLmBZ7wROTr1VgxSNpNWeifzznAfckr7bzcD5qZfmROBkiputL7Q5yyzpRcBHgLMj4vFsfu3fyjIp84bs5dnAfWn6BuDVqezrgFenefO32K3JS9B6vR64Gdiano9I808DPpamfxW4m6Jl+27gwmz7kyj+qLcB1wFDy6TMbwL2AXdlj1PTslvS5/gu8Glg7QKX918BD1D8L3xpmncZxY8PYDh9d9vSd3lStu2labv7gdcs4t/FXGW+CXgs+243z/W3sgzK/GfAPalsXweem23779L3vw14c6/K5FPZzayx1XAIY2YLxAFiZo05QMysMQeImTXmADGzxhwgZtaYA8TMGnOA2IKQdIKk70n6mKTvSvqMpDMkfTONSbEYV93aAnOA2EJ6DvAB4IXAc4HfBv4lxbU5713CclmPOEBsIT0cEXdHcZXwPcDNUZz6fDdwwpKWzHrCAWILaTybns5eT3OQXAm+2jlAzKwxB4iZNearcc2sMddAzKwxB4iZNeYAMbPGHCBm1pgDxMwac4CYWWMOEDNr7P8Dc42VUDHl1zUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAE0CAYAAADUu30UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG6NJREFUeJzt3XuYZVV95vHvW/eu7qavSNOKkAiOiqOMKN4NDkSQqCijiHgflciMmXGijtFM1DBRY55kHBEniEpAFDBiUIx4QXkEVESBByKIBhBItw1K36q6urruv/ljrXPYdayqdbrrXvV+nqeeOmfvffZe+1zes/bea62jiMDMbCot810AM1v4HBRmVuSgMLMiB4WZFTkozKzIQWFmRQ4KO2CSzpf0F3O0redKultSn6SXz8U27REOikVC0v2STmxy2e9Leutslyki3h4R/ztv83hJW2dxc+cA50XEqoj46lztoyUOClssDgfunKmVSWqdqXUtBw6KRUjSmyT9QNLfStol6T5JL87zPgw8HzgvV9PPy9OfIOkaSTsl/VLS6ZX1XSTpU5K+IWmPpJskPS7Pk6SPS/qtpB5J/yLpyZXH/ZWklcA3gc15m32SNkvql7Shsp1jJT0sqX2CfTpO0o2Sdkt6UNJ5kjryvHuB3we+ntf90QPcx7+XdLWkvcALZ/RFWeoiwn+L4A+4Hzgx334TMAy8DWgFzga2Acrzvw+8tfLYlcAW4M1AG/A0YDtwdJ5/EbATOC7P/yJweZ53EnALsBYQ8ETg0Mrj/irfPh7Y2lDmq4GzK/c/Dnxykv07FnhW3v4RwF3AOyfa/2nsYw/wXNIXZNd8v6aL6c81isXrgYj4TESMAhcDhwKHTLLsS4D7I+IfImIkIm4FvgK8srLMP0XETyJihBQUx+Tpw8Bq4AmkILorIh5ssowXA6+DelX/NcAlEy0YEbdExI9z+e4HPg38QZPbaXYfvxYRP4yIsYgY2I91L3tt810AO2AP1W5ERL8kgFWTLHs48ExJuyvT2hj/oX2ocru/tq6IuDZX7T8FPFbSlcC7I6K3iTJ+DThf0u8Djwd6IuInEy0o6fHA/wGeDnTn8t3SxDZqmtnHLfuxPqtwjWJpauwSvAW4LiLWVv5WRcTZTa0s4tyIOBY4mvSBf08T2yR/a/8j8Frg9UxSm8j+HvgFcFREHAS8n3SoM2mxGu43s4/uKn2AHBRL029IJ/9q/hl4vKTXS2rPf8+Q9MTSivJyz8wnIPcCA8DoJNvcIGlNw/TPk86pvAz4whSbWg30An2SnkA67zKVGdtHK3NQLE2fAF6Zr4icGxF7gBcBZ5BOej4EfAzobGJdBwGfAXYBDwA7gL9tXCgifgFcBvwqX7nYnKf/EBgDbs3nHibzbuBMYE/e3pfmcB+toHaW3GzWSLoWuDQiPjvfZbED46CwWSXpGcA1wGH5W98WIR962KyRdDHwXVJ7CIfEIuYahZkVuUZhZkUOCjMrclDMoNnoaj2XYz7MBo8jsTQs6aDIPR/fk9+o+yT9m6SP1HolLjS1XqHVaVEZ82E/1vM+SVc3TLt7kmlnHHiJmzJuHImZXnkz43RIOkjS/82vf5+ke/L9jTNdnpJaj9u53u50LemgAM4FzgLeQGr592LgRODyuS6IpLnsV3M98NzcEQtJm4B24GkN047My86mSceRyEE+q+/B/KXwPVLz85NJDcieQ2o4dtwBrG9e+0fN2/bnu/vqbP0BR5GaGh/XMP0wYBD4g3z/+4zvrvwm4AeV+58g9SPoJXVSen5l3gpS9+VdwM9JfSC2VubfD7wX+Je8zTbgz4B7SS0Qfw68Ii/7RB5pHt0H7M7TLyJ35c73TwVuy+W5Fzh5gn3vIHXsOjbfPx34B+C6hmn3lPYT2AzsA9ZXlv0PpC7c7fn+fyZ1C98FfBs4PE+/l9Qqc1/ep878fH8Y+GGefmTexlWkru73AG+rbOtDpP4in8/P2Z3A0/O8SxrW/z8neC7eSmruvWqK90oAR1bu159zcvf5/Do+lLd5F/CSyvJt+fl4Wr7/5bxsDymIa13dzyL1xh3K5f36gWw/T39Jfh/sBn4EPGU2P09LuUZxAulDO663YkRsAX5Mau7bjJ+SulyvBy4FviypK8/7IPC4/HcS8MYJHv8a4I+AtZG6cN9LGnRlDfCXwBckHRoRdwFvB26MVE1f27giSceRPjDvIY0P8QJSGI0TEUPATXk++f8NwA8aplVrExPuZ0RsA24E/lNl2TOBKyJiOJ93eD9wGnBw3s5luRyPA/4NeGnep8H8+NeTPjSrSc3CLyN9GDaTuoV/RNIJle29jFQLXEsKlPPy+l/fsP6/aXwuSDXIb0VE3wTzmrWJ9Lwcnst9Gel1rTkJ2B6pazukQXyOAh4F3Erqtk9EXJBv/00u70sPZPuSngZcCPwxsIHUJf8qSbPWXH0pB8VGYLJxEx4kvamLIuILEbEj0hgHf0f6Vvx3efbpwIcjYmcOoHMnWMW5EbElIvbl9X05IrZFGhPhS8DdNF8FfgtwYURckx//60h9LCZyHY+EwvNJH+AbGqZd1+R+Xkr+YCj1Zz8jT4P0Zv1opHEqRoCPAMdIOnyK/bgoIu7My28Cnge8NyIGIuI24LOkMKn5QURcHWnsjUuAp06x7kYbmPx90Kwx4IMRMZhfx0uBl0nqzvPP5JHng4i4MCL25GD8EPDUCTrLTWf7bwM+HRE3RcRoRFxMqrE+axrbmNJSDortpMFcJnIo8HAzK5H0Lkl3KQ0Dt5tUE6idBNvM+DEOHphgFePGQJD0Bkm35Y5Tu4EnV9ZXchipRtKM64HnSVoHHBwRd5OqqM/J055MpUZR2M8rgGfnjl4vIFWVb8jzDgc+UdmfnaTu4Y+eomzV52QzsDPGt9x8oOHxjWNldO3HsfoOJn8fNOvhqAx0ExH3kA4/XprD4mXkoJDUKumvJd0rqZdHanzTOXE6bvuk5/xdtec8P++HkZ7LWbGUg+Ja4LBcXa+TdBgpeWvfpntJA6XUbKos+3zSseHpwLp8ONDDI+MkPEh6gWoeO0E56k1f87fsZ4B3ABvy+u6orK/UTHYL6TCnGTeSPuxnkc4HEGmwmW152raIuC+Xa8r9jIjdwHfy/DOByyIfKOcy/XGMHwdiRUT8aIqyVfdzG7Be0urKtMcCv25yP0vP2XeBk5TG9ZxMP5O8B6bYRu3w41Tg5zk8ID0/p5IOedaQhvWDqV/j/d3+FlJNtvqcd0fEZROse0Ys2aCIiH8Fzge+KOlZOemPJg2P9iPSGwjSCaHTJHVLOpJUva9ZDYyQah9tkj5AOmte84/A+yStk/QY4E8KxVpJetEfBpD0ZtI3e81vgMdMcfn2c8CbJZ0gqUXSo/PYDRPt/z7gZuBPeeTbH9J5ij9l/PmJ0n5C+sZ8A+lcxaWV6eeTnoOj8z6tkfSqSco/UTm3kF6Pj0rqkvQU0mvwxSZX0TguRaNLSB+srygNvtsiaYOk90s6JS9zG3Bmfo+cTHND8F1OOs91NuOfj9Wkw4AdpA//R5oo7/5u/zPA25XGCZGklZL+qCFsZ9SSDYrsHaTj3S+QUvsOUrX25RExlpf5OOks9G9IYzxW36DfJp2Y+tf8uAHGV5v/Mk+/j/SNO9UITkTEz4G/I33b/wb49+Rv++xa0ln9hyRtn+DxPyENHvtx0jf+daRq6GSuI51Qq7bNuCFPqwZFaT8hnUQ8CvhNRNxeKdOVpHEfLs9V7TtIl6H3x2tI37zbgCtJx+PXNPnYjwL/K1fB3904M58nOJE0etY1pKs6PyEdCtyUF/vvwEtJVxBeCxTbe0QaN/RG0qXW6tgZnyc9h78mXdX6ccNDPwc8KZe3tp392n5E3Ew6T3Ee6UrTPaSrdbNmWXUKk3QO8HLgBbk6bWZNWFZBASDpHaT2A9+a77KYLRbLLijMbP8t9XMUZjYDHBRmVuSgWEZy79SfKf0m6ENKv8X5O03F56gskS9H2yLgoFgmJL2LdBnzPaSGQM8iXVq9Zop2Gwe6rVnt4TjfPTiXIwfFMiDpIFKbjz+JiG9FxHCk39g4nRQWr5P0IUlXSPqS0i+a3yrpqZV1bJb0FaVfI79P0n+rzKs99gu5LcWbNPWvk9facNyuND7Eq/P0tymNFbFT0lW5yXhtGyHpv0q6m9Q/xuaQg2J5eA7QBfxTdWLuUflN4A/zpFNJXaRrPUi/qvSLWy3A14HbSX0wTgDeKemkyupOJfUJWUtqtDYK/A9Sw6Zn58f8l7zdWse0p+ZelF+S9B9JjadOJ/XNeIDfHTfk5cAzgScd8DNhB8RBsTxsJHWDHplg3oM80mHploi4IiKGST8Y3EU6RHkGqWPZORExFBG/IjUjro6OdWNEfDX3at0X+//r5K8l9Yy9NbemfB+pI9oRlWU+mnvq7tvfJ8Cmx8d6y8N2YKOktgnC4tA8HyrNtiNiTGn8z82k/imbNf6XwlsZ34eksZfs/v46+WbS2A217fdJ2kGqwdw/0TZs7rhGsTzcSOqodFp1Yu5R+WLSUHFQ6QmbDzceQ+p/sQW4r6G34uqIOKWyusaWe/v76+TbqPRbyWXbwPhepG4dOE8cFMtARPSQTmZ+UtLJ+bzDEaTzEVt5pDPbsZJOy1cV3kkKlx+TOlH1SnqvpBW5l+OTlX4ucDKlXydv7EV5Kaln7DFKIzV9BLgppv5hY5sjDoplIg8T937SL5H3knpObgFOqAxR9zXg1aQeia8HTstXSEZJvRuPIfWU3U7qlTvVqE2lXyf/EHBxvipyekR8D/gL0jAAD5LG3ZjtEcKtSe7rYUC6xEka4PV1810WW3hcozCzIgeFmRX50MPMilyjMLOiBd3gSpKrO2azLCKmat8CuEZhZk1wUJhZkYPCzIocFGZW5KAwsyIHhZkVOSjMrMhBYWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysyEFhZkUOCjMrclCYWZGDwsyKHBRmVuSgMLMiB4WZFTkozKzIQWFmRQ4KMytyUJhZkYPCzIocFGZW5KAwsyIHhZkVOSjMrMhBYWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysyEFhZkUOCjMrclCYWZGDwsyKHBRmVuSgMLMiB4WZFc1IUEi6UNJvJd0xyfzjJfVIui3/fWAmtmtmc6NthtZzEXAe8PkplrkhIl4yQ9szszk0IzWKiLge2DkT6zKzhWcuz1E8W9Ltkr4p6eg53K6ZTdNMHXqU3AocHhF9kk4BvgocNdGCks4CzpqjcplZExQRM7Mi6QjgnyPiyU0sez/w9IjYXlhuZgpnZpOKCJWWmZNDD0mbJCnfPi5vd8dcbNvMpm9GDj0kXQYcD2yUtBX4INAOEBHnA68EzpY0AuwDzoiZqsqY2aybsUOP2eBDD7PZt2AOPcxscXNQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysyEFhZkUOCjMrclCYWZGDwsyKHBRmVuSgMLMiB4WZFTkozKzIQWFmRQ4KMytyUJhZkYPCzIocFGZW5KAwsyIHhZkVOSjMrMhBYWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysyEFhZkUOCjMrclCYWZGDwsyKHBRmVuSgMLMiB4WZFbXNdwFscZNUXCYi5qAkNpscFLbfauHQ+L/RRAHh0FicHBTWtGowTHa7JiLqoVC7HRFIclgsQg4Ka0pjMEiipaVl3P3aMrUgGBsbqwfE2NjYuPU5LBYXB4UV1QKgFgwtLS20trbW/7e2ttanwyM1iNHRUcbGxhgdHa3froaGw2LxcFBYU2pBUPtrb2+ntbWVtrY22tvbx9UuamEwMjLCyMgIo6OjDA8PjwuLGofF4uCgsElVaxItLS20tbXR2tpKe3s7nZ2ddHR00NHRQVdXV32epHpIDA4OMjg4yPDwMAMDAwwPDzMyMgKkwxLXLBYPB4VNqXo+ohYSHR0ddHd309XVRVdXF93d3XR2do4LiuHhYfr7+xkYGGBgYKB+IhP4nVqFT3AufA4Km1D15GW1RlGrTXR3d7Ny5Uq6u7tZu3YtnZ2dtLe3I4nR0VGGhobo7e1l7969tLW1MTo6ClA/d+FzFYuLg8ImVa1NtLW10dnZyYoVK1i1ahUbN25kzZo1rF27loMPPrheq2hpaWF4eJh9+/axY8cOenp66OnpAajXOKpB0Vi7sIXJQWGTmuiwo7Ozk66uLg466CDWrVvHhg0b2LRpE6tXr2bFihVIYnh4mL1799aviAD09fUxMjLC0NBQffrIyMiEl1Vt4XFQ2IQmO/SoHXasXr2adevWsXHjRg455BDWrl1Ld3c3ra2tDA4O0tfXVz+0GB0dZffu3QwNDTE4OMi+ffvqIWGLw4wEhaSTgU8ArcBnI+KvG+Z3Ap8HjgV2AK+OiPtnYts28xobVzUeeqxcuZL169dz8MEHc+ihh3LEEUewfv16Vq1aRUtLC4ODg/XDjZaWFkZHR9m1axeDg4MMDAzUr5C0tLQwNjbmk5mLwLSDQlIr8CngD4GtwE8lXRURP68s9hZgV0QcKekM4GPAq6e7bZsb1cOPtra2+iXR7u5uVq1aVT9XUQ0KgFWrVrFy5UpWrFhBR0cH7e3ttLW11dtc1NZtC99MdDM/DrgnIn4VEUPA5cCpDcucClycb18BnCC/QxaFifp1VK+AdHR01M9brFixon5Ss/bX3t5eb5xVO4RpbPJtC99MBMWjgS2V+1vztAmXiYgRoAfYMNHKJJ0l6WZJN89A2WyW1K5YVJtoV5tq16bXrm5UO4b5MGPxmYlzFBN9LTS+E5pZJk2MuAC4AECS31HzrPrhHh0drTfLrp6Y3Lt3L11dXUC6BDowMEBfXx/9/f3s27ev3jqzFiSN4WEL30wExVbgsMr9xwDbJllmq6Q2YA2wcwa2bbOs2vuzFhKDg4P09/fT19dHb28vu3btAmBgYKB+1aOnp4fdu3ezZ88e+vv76825R0ZGHBKL0EwExU+BoyT9HvBr4AzgzIZlrgLeCNwIvBK4NvwuWbAax5Go1SSGhobYt28ffX197Ny5s35FA2DHjh2sWLGClpYWhoaG6Ovr48EHH2T79u1s376d3bt315t012oXPhRZPKYdFBExIukdwLdJl0cvjIg7JZ0D3BwRVwGfAy6RdA+pJnHGdLdrs69am6j1AK3VJnp7e+uNqVpbW9mzZ8+4lpn9/f08/PDD7Nq1a1xIDA0N1XuUNg5uYwuXFvIL5HMU86d6ZaKtrY22trZ6R7Du7m7WrVvHmjVrOOigg1i/fj3d3d10dHTQ0tLCyMgIAwMD7Ny5kz179tDb28v27dvp7++v/w0PDzM8PDzufIXNj4goXn5yy0ybUPWDW6tRDA4O1g9FIJ2T2LNnDz09PXR2dtLW1lZffmhoiL1799ZrEr29vQwODjI0NPQ7hx628DkobEq1b3xJ9bEkAPr7++uHI0NDQ/VaR+0xIyMj9SsetRaZtVpEtTOYw2JxcFDYpKqD4dZqFTW14KhdJq01y649rlYDGRkZqYdJ7dyEr3osPg4Km1J1oNza/Vpo1HqBNjbLri1fuxRau2rS2BCrun5b2BwU1pSJWlhWW2U2dhevjjdR/fOVjsXJQWFFjYPhVs9bTDVcf2OzbZ+TWLwcFNa0xh/wqd0v/VKYw2Hxc1DYfmn8wFcHzS0t77BYvBwUdkCa+dA7GJYOB4VNmwNh6ZuJ8SjMbIlzUJhZkYPCzIocFGZW5KAwsyIHhZkVOSjMrMhBYWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysyEFhZkUOCjMrclCYWZGDwsyKHBRmVuSgMLMiB4WZFTkozKzIQWFmRQ4KMytyUJhZkYPCzIocFGZW5KAwsyIHhZkVOSjMrMhBYWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysaFpBIelVku6UNCbp6VMsd7+kn0m6TdLN09mmmc29tmk+/g7gNODTTSz7wojYPs3tmdk8mFZQRMRdAJJmpjRmtiDN1TmKAL4j6RZJZ021oKSzJN3sQxSzhaNYo5D0XWDTBLP+PCK+1uR2nhsR2yQ9CrhG0i8i4vqJFoyIC4AL8rajyfWb2SwqBkVEnDjdjUTEtvz/t5KuBI4DJgwKM1t4Zv3QQ9JKSatrt4EXkU6CmtkiMd3Lo6+QtBV4NvANSd/O0zdLujovdgjwA0m3Az8BvhER35rOds1sbili4Z4G8DkKs9kXEcXLlm6ZaWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysyEFhZkUOCjMrclCYWZGDwsyKHBRmVuSgMLMiB4WZFTkozKzIQWFmRQ4KMytyUJhZkYPCzIocFGZW5KAwsyIHhZkVOSjMrMhBYWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KMysyEFhZkUOCjMrclCYWZGDwsyKHBRmVuSgMLMiB4WZFTkozKzIQWFmRQ4KMytyUJhZkYPCzIocFGZW1DbfBSjYDjww34WYARtJ+7IcLJd9XSr7eXgzCykiZrsgy56kmyPi6fNdjrmwXPZ1uexnjQ89zKzIQWFmRQ6KuXHBfBdgDi2XfV0u+wn4HIWZNcE1CjMrclCYWZGDYhZIepWkOyWNSZr0Epqk+yX9TNJtkm6eyzJOh6STJf1S0j2S/myC+Z2SvpTn3yTpiLkv5fRJulDSbyXdMcn84yX15NfvNkkfmOsyzhUHxey4AzgNuL6JZV8YEccslmvyklqBTwEvBp4EvEbSkxoWewuwKyKOBD4OfGxuSzljLgJOLixzQ379jomIc+agTPPCQTELIuKuiPjlfJdjlhwH3BMRv4qIIeBy4NSGZU4FLs63rwBOkKQ5LOOMiIjrgZ3zXY6FwEExvwL4jqRbJJ0134Vp0qOBLZX7W/O0CZeJiBGgB9gwJ6Wbe8+WdLukb0o6er4LM1sWel+PBUvSd4FNE8z684j4WpOreW5EbJP0KOAaSb/I32IL2UQ1g8Zr7M0ssxTcChweEX2STgG+Chw1z2WaFQ6KAxQRJ87AOrbl/7+VdCWpWr/Qg2IrcFjl/mOAbZMss1VSG7CGJViFj4jeyu2rJf0/SRsjYil0FhvHhx7zRNJKSatrt4EXkU6CLnQ/BY6S9HuSOoAzgKsalrkKeGO+/Urg2liCLfskbaqde5F0HOnztGN+SzU7XKOYBZJeAXwSOBj4hqTbIuIkSZuBz0bEKcAhwJX5fdYGXBoR35q3QjcpIkYkvQP4NtAKXBgRd0o6B7g5Iq4CPgdcIukeUk3ijPkr8YGTdBlwPLBR0lbgg0A7QEScTwrBsyWNAPuAM5ZiIIKbcJtZE3zoYWZFDgozK3JQmFmRg8LMihwUZlbkoDCzIgeFmRU5KOyASTpC0i8kfVbSHZK+KOlEST+UdHdurWhLgIPCputI4BPAU4AnAGcCzwPeDbx/HstlM8hBYdN1X0T8LCLGgDuB7+VmzD8DjpjXktmMcVDYdA1Wbo9V7o/hvkRLhoPCzIocFGZW5N6jZlbkGoWZFTkozKzIQWFmRQ4KMytyUJhZkYPCzIocFGZW9P8BlagRuTCePBEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "wf = poppy.PhysicalFresnelWavefront(beam_radius=w_extend*w0*u.m,\n",
    "                                    wavelength=wavelength,\n",
    "                                    npix=npix,\n",
    "                                    oversample=2,\n",
    "                                    M2=3.827)\n",
    "wf *= poppy.GaussianAperture(w=w0*u.m)\n",
    "wf *= poppy.QuadraticLens(f_lens=30*u.m)\n",
    "wf.scale_power(P0)\n",
    "\n",
    "print(\"P(0m) = {0:.2f} kW\".format(wf.power*1e-3))\n",
    "print(\"r(0m) = {0:.2f} cm\".format(wf.radius[0]*1e2))\n",
    "\n",
    "plt.figure(1)\n",
    "wf.display()\n",
    "\n",
    "wf.propagate_fresnel(100*u.m, attenuation_coeff=1e-3)\n",
    "\n",
    "plt.figure(2)\n",
    "wf.display()\n",
    "\n",
    "print(\"P(100m) = {0:.2f} kW\".format(wf.power*1e-3))\n",
    "print(\"r(100m) = {0:.2f} cm\".format(wf.radius[0]*1e2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The class `PhysicalFresnelWavefront` implemented an alternative units convention for the electric field and intensity. At the moment, this doesn't change any propagation results. However, if accepted to be included into poppy, I can provide a thermal blooming phase screen (implemented as a `OpticalElement`) which relies on physical units of the intensity because it is a non-linear effect. In its current state, the class does not follow with poppy's conventions, e.g. I did not craft it to use astropy's units package whenever possible. My intent is basically to find out if the community and poppy developers are interested in the capabilites this class provides and features that can build upon it (for example a thermal blooming phase screen)."
   ]
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}