10um linescan with 100x objective (modified from https://github.com/omsai/microscope-optical-input)

raster.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Laser Beam at Back Focal Plane-\n",
      "Measured Power:   0 uW\n",
      "\n",
      "Laser Beam at Objective Specimen Plane-\n",
      "Calculated Power:   0 uW\n",
      "Measured Width: 0.8 um, in FWHM\n",
      "\n",
      "Beam Profile Pixelation at Back Focal Plane-\n",
      "Calibration: 0.150 um/pixel\n",
      "Size: 14 x 14 pixels, covering 3 standard deviations\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Calculate power distribution of laser scanning photo-activation device\n",
    "\n",
    "Laser beam symbols used are from the beam power formula:\n",
    "[1] en.wikipedia.org/wiki/Gaussian_beam#Power_through_an_aperture\n",
    "\n",
    "References:\n",
    "[2] radiantzemax.com/kb-en/Goto50125.aspx\n",
    "[4] en.wikipedia.org/wiki/Beam_diameter#D4.CF.83_or_second_moment_width\n",
    "[5] en.wikipedia.org/wiki/Normal_distribution\n",
    "\n",
    "Theory contributors in alphabetical order:\n",
    "Browne, M. (Andor Technology)\n",
    "Karunarathne, A. (Washington University in St. Louis)\n",
    "Magidson, V. (Wadsworth Institute, NYSDOH)\n",
    "\n",
    "Code Author:\n",
    "Nanda, P. (Andor Technology)\n",
    "\"\"\"\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "\n",
    "## Create the round laser beam profile at the specimen plane\n",
    "\n",
    "# Gaussian laser beam power (micro-Watts) at back focal plane, Symbol: Ps\n",
    "# Measured by user's PD-300 laser meter using single pixel FRAP\n",
    "Ps = 4.8 / 40\n",
    "\n",
    "To = 0.84 # Percent transmission through above objective at 445nm \n",
    "Tm = 0.8 # Percent transmission through media, etc (estimated)\n",
    "# Gaussian laser beam power (micro-Watts), Symbol: Po\n",
    "Po = Ps * To * Tm\n",
    "\n",
    "# FWHM Beam diameter / waist size (microns) at objective, Symbol: Wo\n",
    "# Measured at objective specimen plane using FRAPPA slide\n",
    "Wo_FWHM = 0.8\n",
    "\n",
    "# Convert FWHM to 1/e^2 [2]\n",
    "Wo = Wo_FWHM * 0.8493218 * 2\n",
    "\n",
    "camera_pixels = 512 # e2v CCD97 pixel\n",
    "fov_width = 76.6 # microns\n",
    "# Pixel calibration (microns / pixel)\n",
    "cal = fov_width / camera_pixels\n",
    "\n",
    "# Beam axis pixel length on camera (pixels), Symbol: px_len\n",
    "# Wo for a single mode beam is 4 sigma wide [4]\n",
    "px_len_2sig = Wo / cal\n",
    "# Extend the beam diameter out to 6 sigma to contain 99% of Po\n",
    "px_len = (6 / 4.) * px_len_2sig\n",
    "\n",
    "px_fit = int(px_len) / px_len\n",
    "px_to_edge = int(3 / px_fit)\n",
    "edge = px_fit * px_to_edge\n",
    "steps = int(px_len / px_fit)\n",
    "#x = linspace(-edge, edge, steps)  # keep for debugging 1D construction\n",
    "\n",
    "# Map 1D arrays radially\n",
    "x, y = np.ogrid[-edge:edge:steps*1j,\n",
    "                -edge:edge:steps*1j]\n",
    "z = np.hypot(x,y)\n",
    "\n",
    "# Create 2D gaussian beam profile\n",
    "def norm(x, mean=0, sigma=1):\n",
    "    \"\"\"\n",
    "    Returns single value from Gaussian / Normal distribution curve [5]\n",
    "    \"\"\"\n",
    "    return (1/(sigma * np.sqrt(2 * np.pi))) * np.e ** (-0.5 * ((x - mean) / sigma) ** 2)\n",
    "\n",
    "beam_profile = norm(z)\n",
    "\n",
    "# norm() gives a 1D probability distribution, the area of which is 1.  The 2D\n",
    "# beam profile thus needs to be re-normalized back to 1, for it to be\n",
    "# multiplied by the Po scalar, obj_mag^2, and transmission losses\n",
    "beam_profile /= beam_profile.sum()\n",
    "beam_profile *= Po\n",
    "\n",
    "## Setup ROI\n",
    "\n",
    "w_microns = 10\n",
    "w = int(np.rint(w_microns / cal)) # pixels\n",
    "h = 1 # pixels\n",
    "\n",
    "length = beam_profile.shape[0]  # length in pixels of beam_profile square\n",
    "x = w + length - length % 2\n",
    "y = h + length - length % 2\n",
    "roi = np.ones((x,y))\n",
    "\n",
    "## Laser beam tiles ROI with overlap\n",
    "\n",
    "roi = signal.fftconvolve(roi, beam_profile, 'full')\n",
    "\n",
    "\n",
    "## Show results\n",
    "\n",
    "# Area calculation\n",
    "w_um = w * cal\n",
    "h_um = h * cal\n",
    "x_um = x * cal\n",
    "y_um = y * cal\n",
    "actual_area = x_um * y_um\n",
    "drawn_area = w_um * h_um\n",
    "\n",
    "# Power calculation\n",
    "peak_pixel_power = roi.max() # uW\n",
    "actual_power = roi.sum() # uW\n",
    "drawn_power = roi[length/2:w+length/2,length/2:h+length/2].sum()\n",
    "density_drawn_power = drawn_power / (drawn_area / 1000**2) # uW/mm^2\n",
    "\n",
    "# Energy calculation\n",
    "dwell_time = 80  # micro-seconds\n",
    "repeats = 1\n",
    "peak_pixel_energy = roi.max() * dwell_time * repeats / 1000 # uJ\n",
    "actual_energy = actual_power * dwell_time * repeats / 1000 # uJ\n",
    "drawn_energy = drawn_power * dwell_time * repeats / 1000 # uJ\n",
    "density_drawn_energy = drawn_energy / (drawn_area / 1000**2) # uJ/mm^2\n",
    "\n",
    "print('Laser Beam at Back Focal Plane-')\n",
    "print('Measured Power: {0:3.0f} uW'.format(Ps))\n",
    "\n",
    "print('')\n",
    "print('Laser Beam at Objective Specimen Plane-')\n",
    "print('Calculated Power: {0:3.0f} uW'.format(Po))\n",
    "print('Measured Width: {0:1.1f} um, '.format(Wo_FWHM) + 'in FWHM')\n",
    "\n",
    "print('')\n",
    "print('Beam Profile Pixelation at Back Focal Plane-')\n",
    "print('Calibration: {0:1.3f} um/pixel'.format(cal))\n",
    "print('Size: {0} x {1} pixels, '.format(beam_profile.shape[0],\n",
    "                                        beam_profile.shape[1]) +\n",
    "      'covering 3 standard deviations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb28f161210>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn0AAAESCAYAAACIOWPfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXmYXFW1t98fAcIUUEaRIQgBBJTpKuaKShSVQQVUVBKZ\nxOvlu4BwHRCZBBQU8SqKiuIAgjIJMolRECWIIBIIYUqAMIUpIPMUhAy/74+9Kzmp1Kmu6tPd1dW9\n3uc5T9fZZw/rnK5evc7ae68l2wRBEARBEARDmyU6LUAQBEEQBEHQ/4TRFwRBEARBMAwIoy8IgiAI\ngmAYEEZfEARBEATBMCCMviAIgiAIgmFAGH1BEARBEATDgDD6AgAkfVTSQ5JekLRFg+vzJa3fCdkG\nEknvlHRPfg67SJooaa98bR9J13ZaxiAIgiDoDUt2WoChgqQHgdWBucDLwJ+AA23P7qRcbfAd4ADb\nl5dcHy4BHb8OnGL7R/n8srrrw+U5BEEQBEOM8PT1HQY+ZHtFYGvgbcBRnRBE0oheNBsNTGvWbS/F\nGTS0+Fx6eg5BEARB0JWE0de3CMD2LOCPwFsAJK0p6VJJT+epw//K5SMlzZa0cj4/UtIcSSvk869L\n+l7+vLSk/5M0U9IsSadKGpmvbSfpYUlfkTQLOH0xwRJHSXpQ0uOSfiVpVO73RdJ34TZJM5rc34ck\n3SfpX5JOqut/P0nT8j3+UdK6hWvfz1PHz0uaLOldhWvHSPqtpF/nKdVbJW0o6auSnsjyvr/0gUsP\n5Lp35rF/KWnpZs9F0uckzZD0lKRLJL0hl98LvAm4PMuylKSrJe1XMvabJV2Zx50u6RNNnl0QBEEQ\ndJQw+voBSesAOwNTctH5wEPAG4BPAN+UNM72q8CNwHa53nuAB4Ft8/l2wKT8+dvAGGDz/HMt4GuF\nYd8AvA5YF/jvBmJ9Btg797k+MAr4se3XbI8iGaxvtb1hk1vbjeTF3BrYtWYMSdoV+Gq+vhpwLXBu\nod2NWe7XA+cAF9QMs8yHgTOz/FOBP2d53ggcD/ysiUwAE4APABsAG7Ooh3WR5yLpfcA3gd2BNUm/\nl/MBbI8BHiZ7bG3PKRtQ0nLAlcBvgFWBPYAfS3pzD7IGQRAEQUcIo69vuUTSM8DfgKuBb0laG/hP\n4DDbc2zfCvyCZICR626Xpx43B07J5yOBt+frAJ8DvmD7edsvAycC4wtjzwOOyWO82kC2CcD3bM/M\n6wwPB/aQVPwO9DSFe2Ie/xHg+4Xx9we+Zfse2/OzbFtm4xfb59h+zvZ82ycDI0nGWY1rbV+V214A\nrJzHmgecB4yWtGITuX5o+zHbzwEn9PBcJgC/tH1rNuoOB/6z6Jls4TlAMlQfsH2WE7cCF5GM+iAI\ngiAYdITR17fsantl22+y/flsZLwReKZuQ8dMkqcO4BrgvSTv2W0kL9c4YCwww/ZzklYDlgNulvRM\nNiz/CKxS6PPJZp6pLMfMOhmWBNZo4/4eqWv/xvx5NPCDgmxPk9Y4rgUg6ct56vdZSc8CK5K8YzWe\nKHx+BXjKtgvnAlbohVyw+HNZ5DlkA/ppFv4+WmU0MLZ2z/m+JpA8i0EQBEEw6Ijdu31LIw/RY8DK\nkpbPBgakqcZH8+frSV6vjwLX2L4re512JhmEAE8Bs4HN8nrBRvS0q/QxkqFSYzQwh0UNrp5YB5he\naP9Y/vwwcLztc+sb5PV7hwLvtT0tlz1D324MWafwuSgXLP5cFnkOkpYnGc+P0B4PA5Ns79BmuyAI\ngiDoCOHp62fyVOj1pKnekZI2Bz4L/DpffwW4GTiQhUbe9cD/q51nr9fPge9nrx+S1pL0wTZEORf4\ngqT18kaRE4Dz8pRqqxwq6XV52vZg0tQrwE+BIyRtmmVbSdLu+dooknH5dN408rVc1pccmJ/HysAR\nBbkacS7wGUmb5yn0bwI32H64zTEvBzaStKekJfOmj7fFmr4gCIKgVSTtKOkupU2eh5XUOSVvPpwq\nactC+eF5E+Ntks6uWyvfkDD6+o5mnrbxpF2hjwG/A462fXXh+jXACNKGh9r5CixczwdwGHAvcIOk\n50ibCDZqQ77TSYbm34D7SJ7Dg1uUv3b9UpKBOgX4fe4T25eQ1vGdl2W7Ddgxt7siH/cAD+Rx2zWw\nepLtHNLzuBeYQTJoG3dk/wU4mrT+7lHS72WPJmM1HNv2S8AHc9vH8nEi0OMfXRAEQRDkNfU/AnYA\nNgPG1zsOJO0EbJA3We5PcrIgaTRprf9WtjcnzdwW/5c1HnPh0qkg6D4kPQB81vZfOy1LEARBELSK\npLGkjYY75fOvkib3vl2o81Pgatvn5/PppHX/rwH/IG0UfRG4GPiB7auajRmeviAIWuZ1ktXe8WCn\nZQ6CIOgP+kAfrsWiM1+PsPimwvo6jwJr2X4W+C4p7NijwHM9GXwQGzmC7idc1QPI88CxbdQ/dtHN\nQ0EQBEOGTupDSesDXyD1+TxwoaQJts9p1i6MvqCrsb1+p2UYbizbaQGCIAgGCc304X3A/c2bP0qK\n5lFjbRZG9ijWWadBne2A62w/AyDpIuCdpDXupcT0bhAEbbFkG0cQBMFQppn+2xjYqXA0YDIwRtLo\nvPN2D+CyujqXkZM55DWAz9l+AribFCt2GUkCtmdhSLWm8g5aJMXUXRD0I7bbjpe4VH8IEgRB0IVU\n0Ye250k6iBR9YglStqjpkvZPl/0z2xMl7Zxzw79MSqmK7VslnUWKqDEPuIWeU5YO7t27yeg7pkmN\nSaRNLJ2g2dhV/y02s8X/TEozW0Z/jv1HSt9X+oS5Ta5dQdrVXkazZCRVx/4L6SWqv8ZuxiT67zt+\nXNtGnySf2kb9A+idYRkEQTDY6UZ9OKg9fUEQDD7C0xcEQZDoNn3YsTV9rUShDoJg8BFr+oIgCBLd\npg87YvS1EoW6NdbrU7m6Y+xOblYd08GxN+jg2G/q4NjrdXDsxizVxtGIimmHGraVdJKk6bn+7ySt\nmMvfL+kmSbdKmizpvYU2V+e+bpE0RdKqlR5MEATDjqr6cKDplKdvG2CG7Zm255Bype7afjfr9a1U\nXTF2J42fDTs4dicNzk4a2ut1cOzGLNvGUU/FtEPN2l4JbGZ7S1IqvsNz+ZPAh21vAexLznldYLzt\nrWxvbfup9p9GEATDmSr6sBN0yuhrJQp1EASDkIpvtq288O0KnAVg+5/ASpLWaNbW9lW25+f2N5Bi\nWWH7VtuP5893AstIKooWYauCIOg14ekLgmBIU3ENS2/SDtXqtPqyuB9pq/kiSNodmJINxhq/ylO7\nRzUWNwiCoJxuW9PXKTlaiUKdmVT4vB6DcborCLqDB/NRjWZvrLfmo49pOcSBpCOBOfWpiCRtBnyL\nReMdTbA9S9LywEWS9rT9mz6ROAiCYcFg8eC1SqeMvgVRqIFZpCjU4xtXHTdgQgXB0GY9Fn1puqZX\nvTRTGv+Rjxr1C+iolnZo6WZtJe0L7Ay8r9iZpLWBi4C9bD9YK7c9K/98WdI5pOnjMPqCIGiZweLB\na5WOTO/angfUolDfCZxnu8f0IUEQdJ6Ka1iqpB0qbStpR+BQYBfbr9Y6krQScDlwmO0bCuUjJK2S\nPy8FfBi4o3dPJAiC4Uq3renrmJFq+0+k1HRBEHQRVZRGxbRDDdvmrn9I8gT+OaWh5AbbB5BeLjcA\nvibpGMDAB4HZwBWSlgRGAFcBP69wa0EQDEO6zdPX5WnYBiv9mQptMI9dlWap0HqiP9Ow9ffYnaJ3\nadjaWbO3BZ1POxQEQdAfdKM+7DYjNQiCDjNYpimCIAg6TbfpwzD6giBoi1AaQRAEiW7Th90m7wBS\nxX4fVXHsFSu0Xbni2J2MG/5KhbbPVBz7hQptX6w4dndND3fbm20QBEF/0W36MIy+IAjaotuUXBAE\nQX/RbfowjL4gCNoilEYQBEGi2/Rht8kbBEGHWXaZNir/u9/ECIIg6Djdpg/D6AuCoC2WGtFpCYIg\nCAYH3aYPw+gLgqAtlgytEQRBAHSfPuxIGrYgCLqXpZZs/QiCIBjKVNWHknaUdJekeyQdVlLnFEkz\nJE2VtGUu20jSLZKm5J/PSzq4J3lDLQdB0B5dNp0RBEHQb1TQh5KWAH4EbA88BkyWdKntuwp1dgI2\nsL2hpHcAPwXG2r4H2KrQzyPAxT2NGUZfEATtEVojCIIgUU0fbgPMsD0TQNJ5wK7AXYU6uwJnAdj+\np6SVJK1h+4lCnfcD99l+uH/FDYJg+BFaIwiCIFFNH64FFA21R0iGYLM6j+ayotH3KeDcVgYM9R0E\nQXuM7LQAQRAEg4QO60NJSwG7AF9tpX4YfUEQtEdojSAIgkQTfTjpZZg0u2nrR4F1C+dr57L6Ous0\nqbMTcLPtJ3sWNnbvBkHQLku2cTSgt7vVmrWVdJKk6bn+7yStWLh2eO5ruqQPFsq3lnRb7uv7vX8g\nQatIulrSfn3Y308kHdlX/QV9h6QHJL2v03L0O03037iV4Ng1Fx4NmAyMkTRa0tLAHsBldXUuA/YG\nkDQWeK5uPd94WpzahTD6giBolxFtHHUUdqvtAGwGjJf05ro6C3arAfuTdqv11PZKYDPbWwIzgMNz\nm02BTwKbkN6IT5Wk3OYnwGdtbwRsJGmHKo8lSEh6UNJsSS9ImiXpDEnLtdnHaEnz8++8VraPpGuL\n9Wz/j+0T+kr2IGibCvrQ9jzgIJL+uhM4z/Z0SftL+u9cZyLwgKR7gdOAA2rt89/V+4GLWhU3JmqC\nIGiPDu1WA95U1tb2VYX2NwAfz593ISnSucCDkmYA20iaCYyyPTnXOwvYDbii0t0FAAY+ZPtqSWuS\n/qEdBRzRRh/K/ahB2bBD0ohsIFTpYwnb8/tKpiBT0Yqy/Sdg47qy0+rODyppOxtYrZ3xhrDRt1TF\n9itXaLtWxbE3qdB2w2pDv65a80o8V6XxjIqDT6/Qtn4JRrs8U6HtnIpj94KB3632SC5rpS3Afiyc\n7lgL+EfhWm3n29zcvn6MoG8QgO1Zkv4IvGWxCsnjeiTwX8AywJ+Az9t+EbgmV3tOkoEPkjyzS0p6\nEZhje2VJZwAP2/6apO2A3wAnA4eRfsdH2v5VHm9l4EzgPaSXjCuBcbbf3UC20cADJE/zsbn4e7a/\nm68vDZwEfIJkiF4AfMX2HEmTgB/YvljStsC1JCP4j3m687u2a/HV9gO+DKwB3Ajsb/uhfG0+yQv0\nvyQ/0QYN5Pwt8O78/G4FDrA9LV87A3gFGJ3veVdJN5G85TsCLwO/qHlKJR0DjLG9V90zWNL2fElX\n53t5H7A5cD0wwfYzuf5ewDeA5fPvYHjQZVZUTO8GQdAeFaYzeol6rpIrpvVdc2y3vMYl6D8krQPs\nDExpcPkzpLVK2wHrA6OAH+dr78k/V7S9ou0bgP8H/MP2KNtlb+VvyP28kWRM/ljSSvnaqcCLwOrA\nvsA+9Ow5HEcytnYADiusUTuK9MKxObBF/nxUvnZNble7j/sK97MdMAlA0q6kHZe7kbw117L42qxd\ngbcDm5bINzHLtzrpGZ9dd3088A3bo4DrSAbfKGC9LOPekj5TqF//POrPx5Oe22qkfatfzveyKen5\nfpr07FdhuLxEDbw+rEQYfUEQtMcy5cekV+DYxxYeDaiyW61pW0n7kgyMCS321WxHXFCNSyQ9A/wN\nuBr4VoM6E0jes5l5mupwYI+8jq9m6Lds8GdeIxk582z/EXgJ2Dj3+THga7ZftT2d5PXriWNt/9v2\nHcAZJKOnJvtxtp+2/TRwHLBXvnYNybiDZOx9q3C+HQu9mPsD37J9T552PRHYMhvKNb5p+3nbrzYS\nzvavbM+2PQf4OrCFpFGFKpdmgxnStMCngK/mNjOB7xbkboUzbN+X5fktUNtk9XHg97avy7IczXCZ\nim+iDxc7BgFh9AVB0B7NdqutBsdutPBoQJXdaqVtJe0IHArsUvcP8jKSIbG0pDcBY4AbbT8OPC9p\nmzzNuDdwabUHExTY1fbKtt9k+/MlRssbgZmF85mkb9Ia9N5geLpu3dpsYAWSZ2oEi07p95S9wHX1\nZ2aZyT8fKrn2D9LGoNVJXsCzgHUkrULyCNaMvtHADyQ9kw3kp/OYRQ9ZcfxFkLSEpBMl3SvpOdJU\nrIFVS+5xVdLzrZe7HY/c44XPtWcL6d4XjJWN+Kfb6Ld7qRjNYKDpiNEnaW1Jf5V0p6Tb1UKS4CAI\nBgkd2q1W1jZ3/UPSP6A/KyUgPzW3mUbySEwjTYUdYLtmUBwI/BK4h7RB5E998XgCoDUP3WMkw6fG\naJI36gkaG31VPEdPktb4rV0oW6ekbg3V1VmXJDM0lv0xANuvADcDhwB35E1E/wC+CNxr+9nc5iHS\nGr6V8/F62ysUPHPQ/J4nAB8B3mf7daQpW7Hosy+2f4r0fOvlrnm4XwaKu6wbBxlpzCwKzyrvKl2l\njfbdS5dN73bK9pwLfNH2VEkrADdLurKYZDgIgkFKZ3erLdY2l5fuYLL9LRpML9q+GXhra1IH/cC5\nwFck/YlkkJxAMuTnS3oSmE9ar1bbpfUEsLakpfIUYsvkPi8CjpX0OZKxszeLehobcXR+GVmftAax\ntnTgXOCovDEC0nTmrwvt/kZ6QTkpn08ifQfPKtQ5DfiGpFttT8trDz9g+8IWb2sU8CrwrKTlc/+l\nRmJ+Br8FTpC0D8ko+0JBxqmk38c6wAu0mOEhcyFwg6R3kjzyX6f9qfnuZJB48FqlI54+24/bnpo/\nv0TaOjk8Fn0GQbfTZdMZwYDTzDtVvHY6yVD6G2mzw2zgYFjgLTsBuC5Pf24D/JXk4X1c0r96Icvn\nSfEJZpHW851DMpqacQ1wL/Bn4CTbf8nlxwM3AbeRds3elOUttlsh31vxvDa1i+1LSOv4zsvTs7eR\ndtU2kr0RZ5G8hY8Cd5B20/bEwaTnfH+W7Te2z8jyXAWcn+WYDPy+rm0zg3IayXN+Lsnj+TRNpqaH\nFF2mD7VwpqNDAkjrkd6C3pINwOI1wzG97DlCtvSKCNnSC7o1ZMtx2G7rbVyS/Yk26l9A22MEwUAg\n6URgDdufaXBtNMkwWipi2wVldKM+7Kjtmad2LwQOqTf4FjKp8Hm9fARB0D4P5qMig+SNNQjaQdLG\nwNK2b8+ew8+SYjqWNhkYyYKupsv0YcfElbQkyeD7te0mu+bGDZBEQTDUWY9FX5quaVytJ0ZWlyQI\nOsAo4NycJeQJ4Du266cwiwyPkCNBNbpMH3bSRj0dmGb7Bx2UIQiCdumyN9sgALB9Ey2uf8kx7AbJ\nfstgUNNl+rAj4ubUNJ8Gbpd0C+mN6ogImRAEXUCXKbkgCIJ+o8v0YUfEtX0d8RYVBN1Jlym5IBhI\nJE0A9ra9Y4+Vm/ezSC7cYJDSZfowMnIEQdAeXRaMNBhaSNpD0jRJL0makWeOkDRB0ouSXsjHy5Lm\nS9qqhT43k/QnSU9KmldFPtvnVDX4it01KpS0naR5+T6fl3R3Lbh5Xb1DJd2Tn8WDkr6Zs9nUrp8h\n6et9JOvwpMv0YRh9QRC0R5fFpQqGDpI+QApCvI/tFUi5be+HBcbWKNsr2l6RlMnlPtu3tND1HFKM\numa7eQcbj+Z7XQn4X+BUSZvWLkr6IfBfwJ6kTSw7AduTMtQEfUWX6cMw+oIgaI8uU3LBkOJY4Ou2\nJwPYnmV7VkndfcgZMCQtJekWSQfl8yUk/V3SUbmfe3KQ4mmtCJE9iJ+XdJ+kf0k6qXBtH0nX5s//\nmb2Ha+XzLXKw6Y3y+ZqSLsx93Cfp8714Jtj+Iykg8qa53w2B/wEm2L7R9vycsvDjwI6SxvVmnKAB\nXaYPB4kYzehtkOVRFcetEmB522pDb1khZeHu1YZmbMX2Vbih5yqlXFgxKPXUKsG4r6s2dqUAyy9W\nHLsXdFmIgmBoIGkJ4G3AZZJmkL6JlwJftv1qXd3RwLtJqdOwPUfSnsDfJF1FMn6WYNEsGu2yG7A1\n6Z/NXyTdZfv0fM153H9I+ilwpqQPkzKQHGn7HkkiZb24GPgUKXftVbmfP7cqRO7nI8BKQM2r+T7g\n4ZxqcAG2H5F0A/ABFg2CG/SWLtOH4ekLgqA9uuzNNhgyrEHyAnyc9Ga9JbAVcFSDunsD1+bQKwDY\nvpOUPu0S4IvAnq6WkupE28/bfgT4PjC+pN5xpFxHN5IMsZ/k8rcDq9o+wfY82w8CvwD2aHH8tSQ9\nA7wCXATsZfu+fG1VUrq5RszK14O+oMv0YRh9QRC0R5cpuWDI8Er+eYrtf9l+BvgesHODunsBv2pQ\nfhYwGpho+/6K8hRzy84E3tioku25WZbNSPLWGE023PLxLHA4sHqL4z9qe2WSp/EHwBHZ6wfwFLBm\nSbs18/WgL+gyfRhGXxAE7dFlu9WCoYHt51jU0IIGu1vzbt41gd816OZU0pTqDpLeWVGkdQqf1wUe\na1Qpr+c7BjgD+J6k2pqlh4H7ba+cj9fbXsn2R9oRwvYc4Kuk6d1aeJe/AutIeludLOuQFvFc1c4Y\nQRMq6kNJO0q6K++yPqykzil5p/pUSVsWyleSdIGk6ZLulPSOnsQNoy8IgvbosjfbYEhxBvB5SatJ\nej3wBZIRV2Qf4He2Xy4WStqLtAZvX+AQ4CxJyxWujySt0JKkkcXQJiUcKul12ZA6BDivicw/t/1f\nJMPw+Fx+I/CipK9IWkbSiBw65m0l/ZSSDb/vAofl8xnAacDZkt6RN65sRkp9eqXtq9sdIyihgj7M\n61R/BOxA8gSPl/Tmujo7ARvY3hDYH/hp4fIPSF7rTYAtgOk9iRtGXxAE7VHR6Kv4ZtuwraTdJd2R\nY5dtXSifkHdtTsk/50naPF+blPuqXY91ToOfbwA3AfcAdwI3A9+sXcyG2+7UTe1mw+x7pHVvs22f\nC0wGTs7XR5Omj28neQ9fAe7qQZZL8/hTSIbn6fUVJB0MrAZ8LRftB+wraVvb84EPk9YmPgD8C/g5\nsGLPj6EhpwOrSdoFwPaBpDWCvyHt+JpI8gAWt/tFfuGqVNOH2wAzbM/Mhvt5wK51dXYl70K3/U9g\nJUlrSFoReHfedY7tubZfaEXcIAiC1qmwW63wZrs9yesxWdKltu8q1FnwZpunK34KjO2h7e3AR0ne\njQXYPgc4J/f7FuBi27fVLgPjW4zjFgwC8vq4A/PR6PqrwGJb8W0/TDK+imXjC59n0r4TZKLtHzUY\n60zgzPz5FOCUwrVZpA0ptfPHgQkl93Jc2cC2ryFNKRfLXqFuPaDt7wDfadJPN8UlHJxU2727Fmma\nv8YjJEOwWZ1Hc9k84ClJZ5C8fDcBh+TvQSlh9AVB0B7VtMaCN1sASbU326JXZZE327xuZQ3gTWVt\nbd+dy0Q541l8Ci5mO4Ig6D1N9OGku2DS3f068tbAgbZvkvR90trOY3pqFARB0DrVNmj05s32kVzW\nSttmfArYpa7sV5LmABfZPr5BmyBoREyLBokm+nDcZumocVz96tPktSt6bNfOZfV11imp87Dtm/Ln\nC8lrOpsRRl8QBO3R7M12ejr6mGbeu9Y6kLYBXrZdzLgwwfYsScsDF0na0/Zvqo4VDH1sx970IFHN\nipoMjMlrSmeRYjTWx3u8jLSc4XxJY4HnbD8BIOlhSRvZvoe07KXHjDJh9AVB0B5NtMa4t6ajxnEX\nL1alypvt0i20LWMP4NxiQS19l+2XJZ1D8hqG0RcEQetUsKJsz1NKDXglaanJL21Pl7R/uuyf2Z4o\naWdJ9wIvk7PMZA4m7dBeipSD+jP1Y/ShuEEQDEs69GYr6akW2kKdZzCv8/sk8K5C2Qjgdbafzgrz\nw0DLqa+CwYuk+cCYPgi+POSQtA/wX7bf3WlZhgwVrSjbfwI2riur35B2UEnbW0mZXVomjL4gCNqj\nwsRWlTfbsrYAknYDfkhKL3W5pKm2d8rDvgd4KKe5qjESuELSkvmOriKFywi6n1hv15x4Pn1Jl030\nh9EXBEF7LFOtecU328Xa5vJLSDlVG7W5BnhnXdlsoO0guEFXUHkNaFuDSSNsz+u2voM+oqI+HGgi\nXEEQBO0RadiCLiF7jKdIel7STEnHFK6NlPRrSU9JelbSPyWtlq+tKOkXkh7Li+W/UQsHJGkfSX+X\n9L285GCxEBmS3i5pch53lqT/y+WjJc2X9DlJj+bjS4V2xyil1fq1pOeAfZT4qqR7JT0p6TylbCS1\nNr/NYzybA45vWri2sqTLshw3ABv0w2Me3nSZPuwCT19vRextUPMam/S+6ZarVBv6/3rfdPz2iwWF\nb4vdubBS+ypcuP3uPVcq4dyxFWOMfrnC72xqhe8KAM9UaNs0Dmf/0AVaIwgyL5GycNyZg3P/WdIt\nti8jpWtbkRQK6DVSZozaH9SZpHWj6wMrAJcDD7FwCcA7SEG/VwdquXSL/AD4vu2zlVK9vaXu+jiS\nATYG+GuW6a/52i7A7rb3yhlGDs5l7waeIgV7/jELgzpPJKWWmwN8Gzgb2CpfOxWYTQoIvQFwBWnB\nf9BXdJk+DE9fEATtEbl3gy7B9t9s35k/30EKzr1dvjwHWAXYyIlbbL8kaXVgJ+ALtv9t+yng+yy6\naehR26fanp+zgNTzGmnT0So57duNddePzX3fQcrNW+z7H7Z/n2V+lZRv9Ujbs3Kqrq8Du+cMNdj+\nVR6jdm0LSaPy9Y8BR+ex7iRnCgn6kC7Th2H0BUHQHl02nREMXyS9Q9JfJf0rT5fuT9rsA/Brkufr\nPEmPSDox7+oeTfLezZL0jKRnSakAi7mZi0HCG/FZ0trTu/K08YcK10wKLF5jJvDGJn2PBi7OsjxD\nisU2B1hD0hJZ7nvz/T2Q+1+VlHZuRIOxgr6ky/ThILE9gyDoGkJrBN3D2aTp0B1sz5F0Msm7V8vj\n+w3gG5LWBf4I3J1//htYxXbZTtemO2Bt30eefpX0ceBCSbWcwCLFobwnn69LyiVd1vdDwH62/1E/\njqQ9gY8A77P9kKSVgGfzGE8CcxuMFfQlXaYPO+rpy28pUyRd1kk5giBogy6bzgiGNSsAz2aDbxsW\nroND0jhJb8nToC+RvGfzbD9OCgt0cp4mlaT1Jb2n1UElfVpSzTP4PMmQm1+ocrSkZSVtRgpJVJ8T\nushpwDeknB9FAAAgAElEQVSzYYqk1STV0gmOAl4FnlXKLPOtPBa25wMXAcfmsTYlrWMM+pIu04ed\nnt49hBbShgRBMIgY2cYRBANP0VN2AMmT9zxwFHB+4dobSPlKnwfuBK5mYUaWvUkZYKaRdlpdkOu3\nyo7AnZJeAE4GPlW39u8a4F5SQPCTbP+lSV8/AC4Frsz3cT0Lc06fRfIEPgrcka8V+TzJMJwFnJ6P\noC/pMn3YMdtT0trAzsAJwBc7JUcQBG0ySN5Yg6ARxby4ti8iebsa1TuPEg+b7RdJBuMBDa6dSQ8b\nImzv1ewycLrtXzRod1yDMpM2kny/wbWXgd3qin9TuP4Uafq3yGIhZoIKdJk+7KS4JwOHAit1UIYg\nCNqly5RcEAwyBjR4dNDPdJk+7Ii4eSfTE7anShpH0z+CYjrM9YnYkkHQW+6jT0J0DZJdaEHQpUQa\ntKFEl+nDTtmo2wK7SNoZWBYYJeks23svXvUDAyxaEAxVNmDRl6Zmy4ia0GVvtkEwWLA9k64zE4Km\ndJk+7Ii4to8AjgCQtB3wpcYGXxAEg44uU3JBEAT9Rpfpw07v3g2CoNvoshAFQVCPpBclrddpOZqR\nc/xe26Gxd5DUcANMH/R9kKQT+6PvjtBl+rDjRp/ta2zv0nPNIAgGA3NHtn4EgaRjJJ3VaTmK2B5l\n+8FOy9ECnVr/dzwp5l9TJI3IBvTbC2WfljS/Qdn0fPpzoBjHsKvpNn3YcaMvCILuYt6SrR+NkLSj\npLsk3SPpsJI6p0iaIWmqpC17aitpd0l3SJonaetC+WhJs3MQ+CmSTi1c21rSbbmvxcJhBMFwRNLb\ngBVtT+6pru15pNiAxcDV7wamNyi7Jrd5FZhIioXY9VTVhwNNGH1BELTF3BFLtHzUk7Mf/AjYAdgM\nGC/pzXV1dgI2sL0hKVfqT1toezvwUfI/ljrutb11Popx134CfNb2RsBGknbo9UMZ5kj6vqSHJD0v\nabKkd+XyHUjrtz+VPUK3lLQ/LOe/fUHSdEnvzeXHSLpA0nn52k2SNi+0W1PShTm37n2SPl+4toSk\nI3Je2ppca+Vr8yWtnz+fIenHkiZmGf8m6Q35np6VNE3SFi2OeYyk8yWdmeW9vfgS0uC+50v6fO7n\nX5JOavcZtzJuM5kbsBOFv6P84jQ///3Vyq6WtF8+vZbFDbxvA9vVlf2tcH4NUMxH3LVU0YedYHBI\nEQRB1zBvySVbPhqwDTDD9kzbc0jBcXetq7MrKdMAtv8JrCRpjWZtbd9tewaNwz8tVibpDcCogjfj\nLBYPchu0zo3A5sDrgXOACyQtbfsK4JvA+XlKdav6hpI2Ag4E/sP2iiSj/sFClV1ImTReD5wLXJKn\nFQX8HrgFWBPYHjhEUi3kw5eATwE72l4J2A+Yna/VT5t+gmScrkJKx3YDMBlYGfgdKa4sLYwJKRjy\nOaQYtL8HftzDs9sN2DofuxaMqXoaPuOexm1R5iJvJeUgLtJsmvlvpIgcKE3ZLgf8lpw1JJdtwqJG\n33RgC4YAFfVhb2Y+tiqUPyjpVkm3SLqxFXkHicOxGUv1st3KPVdpyoa9b7p7tZHHb9/7TDnnzPhs\ntcEbfuUGho99+4+9b7x9tbHP3b1Mz7bA1ArfFSDpv97yTMWx22feiEoRJ9YCHi6cP8LClFLN6qzV\nYttGrCdpCind1tG2/577eqTBGEEvsH1O4fRkSUcDG5M8sD0xj5Ty7C2Snrb9UN31m21fDCDpe6QM\nTmNJxtmqtk/I9R6U9AtgD1KA188CX7Z9b5axKEv9i8DFtqfmMS4G/sf22fn8fJJRCun71mxMgL9n\nYxdJvyalG23GibafB57PywzG0yBdWgvPuGzcVmQu8jrgxR5kLvJPYDlJbyXFhPq77X9Lur9Q9oDt\n4t/biwyRxAxV9KEWzl5sDzwGTJZ0qe27CnUWzHxIegdphmJsvjwfGGf72VbH7AKjLwiCwcS8gQ8z\nViWDwWPAurafzdNdlyglng/6EElfJnnS1sxFo4CWFurbvk/S/wLHAptKugL4ou3Hc5WHC3Ut6VHg\njbloLUm1Nx+RZq9qHqV1aD0a+ROFz680OF8hf163hzEBHi98ng0sI2kJ2/NLxi4aQzNZeG+L0MIz\nbjhuizIXeTb33RK2X81epu1IGRRqO46vK5TVjzWK9BLW9VTUhwtmLwAk1WYv7irUWWTmQ9JKktaw\n/QQLf5ctE0ZfEARtMbeJkrt+0lyunzS3WfNHSf+Eaqydy+rrrNOgztIttF2EPA38bP48RdJ9wEZN\nxgjaJK8tOxR4r+1puewZFhrrPe5AreXBlbQC8DPSmrB98uUFv6c8Vbk2yZifB9xve+OSbh8ieZmm\ntXtPTXi4hzF7wzosdPevS7q3RZD0bpo/42a0K/NtpL+RGi/nn8sBL+XPb6hrU1vXtx5pd26tbE/g\nTcCpdfU3AW5tUZ5BTTN92AK9mfl4NJc9Qfrb+rOkecDPbP+cHog1fUEQtMVrjCw93jZueQ4+dqUF\nRwMmA2Py4vClSVNMl9XVuYy8s0/SWOC5/FbbSlso/COUtGptAbrSwv0xpH+Aj5Om07bJhsTewKW9\nfyrDmlGkqdanJS0t6Wss6il6gjTF3tBAkbSRpPfm3+lrJM9a0Sv2H5J2kzQC+ALwb9KauxuBFyV9\nRdIyeZ3fZkq7TwF+CXxD0pg8zlslvb6X91iTvacxm7Ut41BJr5O0DmlK9rwGdVag+TPuS5knAuNq\nJ7afIhkaeyptjtmPxfOh/g14L7CO7ZoBe13uZwsW9/RtB1RYzzN4aKYP649+YFvbWwM7AweqsLmn\njDD6giBoi3mMaPmoJ4d4OAi4ErgTOM/2dEn7S/rvXGci8ICke4HTgAOatQXIRsHDpLUul0uq/UN5\nD3BbXtP3W2B/28/laweSDIN7SFMsf+rrZzVMuCIf9wAPkKYWi56JC0gGyNOSbmrQfiRwIvAkycu1\nGnB44fqlpA0ZzwKfBj5qe16eLv0wsGUe918kL9OKud33SL/zKyU9D/yClPYT2o9/Z4AWxixt24RL\ngZuBKaQNF40Wdff0jPtMZtu3AM+pEGcP+BzwFeApkpfuurpm1+f+bij08zTpd/qE7ftq5ZKWIRkp\nZ/Ygf1fQTP9dN2ku3z129oKjAVVmPrA9K/98EriYFtY4yx68uZ8lGb7by9ZjKo5eIV708dVGHn/k\n8NzIwbd733TChr+sNPS5J1TYyHFUpaFp7KxqlXsrtP0StttaLyfJM7x2y/U31CNtjzHckDSS5AlZ\nOh+X5lSVwx5Jx5AWsQ+JmG71SJoPjLHd6trDASHv7P0f2x/rh74PAta2/dW+7nugqaoPs/f6btJG\njlkkr+z4grcUSTsDB9r+UJ75+L7tsZKWA5aw/ZKk5Ukvw8fZvrKZDLGmLwiCtqi4hiWoIy+Ef6/t\n2fmfwHWStrVd700JggHB9p9pvLO3L/r+UX/02ymq6EPb87IRfCVp5vWXtZmPdNk/sz1R0s555uNl\n4DO5+RrAxck5xpLA2T0ZfBBGXxAEbTIv1EafY7s29zOSpPxbDsEQdDWDd6otaImq+jAvK9m4ruy0\nuvODGrR7gDRl3xahvYMgaIsOhGwZ8uTNJjeTFsj/tLZDc7hj+7hOy9Cf2I4/pi6n2/RhGH1BELRF\ntym5biAvtt9K0oqkjQfb2V4kpVyexgmCoJ/ozfrjbtOHsXs3CIK2eJWRLR9Be9h+AfgD0DCchu1B\ncxxzzDEdlyHkCXn66ugt3aYPw9MXBEFbdNub7WBHKTfpHNvPS1oW+AAwpKc1g2Co0G36MIy+IAja\notuUXBewJWkXXk0fX2H7L50UKAiC1ug2fRhGXxAEbdFtSq4/yJkd1rF9Wx90dwfwbttTcxqymyW9\n2YWk6/3B00/P5vTTb+l1+8cfX4XvfGfwRJUJeZoz1OTZb7+tWGWV5fpQot7RbfowjL4gCNpiuMbp\nkzSJFLV9SdJO239Jus72F6v065QS7vH8+SVJ00m5NfvV6HvyyZf5yleqhmJrpf1SJeXLlpQ3+7dU\n1mYpUlayGS32VSZTX7I6KZVqfzOnpLw+B3bt+ZTVf6XJGGX5tMvalI1RT++/fx/5yEaDwujrNn3Y\nBUZfb0UsUw4t8roKbcdWG3p3Lux944oZNY69uFr7SmNXaLv7RRWeGXDu2AoZOap8VwCeq/JdHfg/\n4WEcp28l2y9I+i/gLNvHSOoLT98CJK1Hmu79Z1/22z7tGmoAK5eUl2VHekfj4jc3McjKdGtZtLL1\nSsqb/c0uU1Je9rXvzZ9DmQ1VVv7vJn09V1L+YEn51JLyG0rKAe4qM+LKvqZlmYKeaTJIVQOyM3Sb\nPuwuaYMg6DivsXSnRegUS0paE/gkcGRfd56ndi8EDrH9UqM6xx577ILP48aNY9y4cX0tRqbsX0OZ\nYQewdePiN2zeuLwkCdcWh5RbH1/g5Ibl+/zrt40bTCzp6PrSIcrtkjLDa16Tvsoocw6VGZzNHvsm\nJeXjGxefecgnG5afzBdKh7j1ByXW9onvalz+eFkq4imlY8ATJeX9Y/RNmjSJSZMmVe6n2/RhGH1B\nELRFt01n9CFfJyW9/7vtyZLWZ/H5xF6RN3FcCPza9qVl9YpGX//Sh56+MkdfiR2xA1eUjlBq3JWY\n4LPPbVx+zculQ5ROyL5Y3qTPGFVS3iy763bLNy5frsTo2+eExs9w2uqblo5x69iSX1bZ7/bxMiu1\n2fen7DvXbNq599S/NB13XO82zHebPuyY0SdpJeAXwFuA+cB+tjs8pREEQU9023RGX2H7AuCCwvn9\nwMf7qPvTgWm2f9BH/XWANtfJlXi1RjRznb1aUl5ixL1YUv5C+Qilxl2zNv1NU3lL7nG5MsO25Bk2\nfe5lHshSBmLN5OCg2/RhS9JK+gTwJ9svSjqK5Mc/3nYzX21P/ACYaPsT+S238ysygyDokW7brdZX\nSFoN+BxppdgC3Wm7woJQkLQtsCcwT9L/ANOBI5xycgZBMIjpNn3Yqol6tO0LJL0LeD/wHeAnlK7E\nbU5ONfRu2/sC2J5LZ1+kgiBokW5Tcn3IpcC1wFX0biVXQ2xfJ+k9wEukDSIli+OCIBhsdJs+bDUN\nW03BfQj4me0/QKXVi28CnpJ0hqQpkn6WI9EHQTDImcuIlo9GSNpR0l2S7pHUcL+5pFMkzZA0VdKW\nPbWVtLukOyTNk7R1ofz9km6SdKukyZLeW7h2de7rlqyHVu3h1pezfZjt39r+Xe1o9bk1w/bfgWf7\noq8gCAaOqvpwoGnV6HtU0mnAp4CJkka20bYRS5KmiH+c32pnU7qXKwiCwcQ8lmz5qEfSEsCPgB2A\nzYDxkt5cV2cnYAPbGwL7Az9toe3twEeBa+qGfBL4sO0tgH2BX9ddH297K9tb236qh1u/XNLOPdQJ\ngmAYUUUfdoJWpfgksCPwf7afy2ELDq0w7iPAw7ZvyucXUhph7o+Fz2OADSsMGwTDmRmUx89qnYoh\nCrYBZtieCSDpPGBXFg1EvCtwFoDtf0paSdIapBmChm1t353LVBzM9q2Fz3dKWkbSUrZrcSDaeXk9\nBDhC0mssjCNh22XxKYIgGOIMqZAtkor7ricVyl4FbmrUphVsPyHpYUkb2b4H2B6Y1rj2Tr0dJgiC\nRdiQRV+aerdPoOIalrWAhwvnj5AMwZ7qrNVi21Ik7Q5MKRh8AL+SNAe4yPbxzdrbLouoMWAMXJy+\nIBja9FWcvm5b09eTp+9mwIAaXDOwfoWxDwbOlrQUcD/wmQp9BUEwQDRbmzJj0mPcO+mxvh6ykf5p\nrwNpM+BbwAcKxRNsz5K0PHCRpD1t/6aHfnYB3pNPJ9m+vKpsxe7p4V4HLk5fEAxtIk5fA2y/qb8G\nztMub++v/oMg6B+arU1Zf9y6rD9u3QXnVxy3WFSnR4F1C+dr57L6Ous0qLN0C20XQ9LawEXAXrYf\nrJXbnpV/vizpHJLXsNTok3QiSWednYsOkbSt7cN7kqEFGa8G3gWMkPQs8EXbZ1TtNwiC/mWwrNVr\nlZbWsyixp6Sj8/m6klqeVgmCYOgwjxEtHw2YDIyRNFrS0sAewGV1dS4D9gaQNBZ4zvYTLbaFgrcs\nB4G/HDjM9g2F8hGSVsmflwI+DNzRw63vDHzA9um2Tyetc/5QD216JG9QWYe0aHkkMBP4R9V+gyDo\nfyrqw0rRDPK1JXL0gUa6cDFaXcR8KvCfwIR8/iLw4xbbBkEwhKii5GzPAw4CrgTuBM6zPV3S/pL+\nO9eZCDwg6V7gNOCAZm0BJO0m6WFScq/LJdV2gB0EbAB8rS40y0jgCklTSQlBHwF+3sLtv67weaXW\nn1pTFmxuyesNaxtUgiAY5FTRh1WiGRQ4hNI9EYvTql/yHba3lnQLgO1n85t2EATDjKprWHKmiY3r\nyk6rOz+o1ba5/BLgkgblJwAnlIjythZFrvEt4JY8FSvS2r6+CDVVaYNKAB+65nImPlbZ6Tq0qMs7\nvPMb/8AftvtwZ2QZwlTUh72OZpA3xK5NmoE4AfhiKwO2avTNkTSCtHmjlo5ofottgyAYQrzGyE6L\nMODkUDB/J3kSa2uRD7P9eOekCmqEwdcz8Yz6h4r6sDfRDB7NZU8AJ5PC57U869Cq0XcKcDGwuqQT\ngN2Bo1odJAiCoUO3hSjoC2xb0kTbb6XxOsIqtLK5BYiQLUHv2fmNf+i0CIOKbg/ZIulDwBO2p0oa\nR4tRDloy+myfLelmUjw9AbvV1tIEQTC8GI5GX2aKpLfbntzH/S7YoALMIm1QGd+oYoRsaYzHN/5/\n98S5DYtT0NkSyraDD0Ry+LIo32s1aTOupHyNht+goEZfhWxppg8fmvQAD016sFnzKtEMdgd2yVmC\nlgVGSTrL9t7NBmzJ6JP0fttXUZhnlrSP7TNbaR8EwdCh2+JS9SHvAD4taSbwMukF2LY3722HOWD0\nsaRsI9eQMn38Ml6qg6A7aKYP3zhuDG8cN2bB+XXH1WeJbOmF7zLgQOD8umgGR+QDSdsBX+rJ4IPW\np3e/JunjwJeBFYBfkLJyhNEXBMOMbotL1Yfs0A991nIGnwZ82fZigQ2DIBi8VNGHtudJqkUkWIL8\nwidp/3TZP7M9UdLOOZrBy1RMZNGqtNsBXwKm5vOv2S5xnAdBMJQZxtO7awJ32n4RQNKKwCakuHq9\noixncBAE3UFVfVglmkHh+jWkmYIeadXoez1pR8l9pPnk0ZJk2y22D4JgiDCMjb6fAFsXzl9qUBYE\nwTCi2/Rhq0bfDcCJtk+XtCzwbeA64J39JlkQBIOSV4dhyJbMIi+6tudL6lGHSvozsEaxiBT+6kjb\nv+97MYMgGCi6TR+2avS93/ZDALZfAQ6W9J4e2gRBMATptjfbPuR+SQeTvHuQMoXc31Mj2x/oKwEi\nZEsQ9A3dHrKltzQ1+iS92fZdwKo5dVGRl/pPrCAIBivdpuT6kP9Hill6FMlT9xfgv/uw/x7X9UXI\nliDoGwYiZMtgpCdP3xdJSu275Gwcmdr0xPv6Sa4gCAYp3abk+grb/yKFVOgzJO0G/BBYlZQzeKrt\nnfpyjCAI+o9u04dNjT7btbfYnUlTGe8iGXvXsnCKIwiCYcRwi9Mn6Su2T5L0QxZ9+QXA9sEVun8n\nadbkadJGuUrhGIIgGFi6TR+2uqbvTFJA8lPy+QRSAuBP9odQQRAMXoZhnL5aoOSbaGD0VeRK4Kt5\nU8iJwOH56ELmtFf9342Lm3pOytbML9+4eFRJ+Yovlw8xEJk3yhhVUl6WqQPK77HsmZQ9w6bPveR3\nVU6b34Uuptv0YavSvsX2poXzqyVN6w+BgiAY3HTbdEZVCjtsp5Ei4K/HQt1p0gtwb/u+qnB6A/Dx\n3vbVt5T9036lSZtnGhffu3bj8hsaF1/xjvIY2Juu3vjfzj4n/LZh+XLbNu5np2b5Tkpuo9Twmdek\nrzLK/oSWKSlfuUlfm5SU79y4+MzVG/tqrmgWe7zkd8W9ZQ3KHmKz7093Gordpg9bNfqmSBpr+wYA\nSe8gvfUGQTDMeJWlK7WXtCPwfRZGoP92gzqnADuRItDva3tqs7aFdGabAG8vZraQdDiwHzAXOMT2\nlbl8a+BXpH+1E23/bw+i/wY4lJRFY35v7r0H9gPO64d+e8HckvKyf+YAJclEHi/xnf3vOxoW3/rT\nsaUj7Dv2/MblWzYuZ72SjkoMIqDc8Cr7b9kbR0/Z4y0rb+Zpe66kvCx9wtSS8jLDDuCuMoPsnyXl\nZdZgs+9P2c0Pbqrqw4Gm1a/rfwDXS3oon68L3C3pdirmngyCoLuoMp0haQngR8D2wGPAZEmX5igB\ntTo7ARvY3jC/YP4UGNtD22I6s+J4m5CWoWxCCix/laQNc7y9nwCftT1Z0kRJO9i+oon4T9m+rBf3\n3GOcPklHAnNsn9Osr4EL2VL2T76ZN6ZsYrQ+f3yNGxsX39Xk+3XXsiUXliopL+urrH43UvY7KTOi\neuPFLeurrM3g99r1XciWoTm9u2O/StGU3lr/zb7ALVD29tQKzd6YWuDC7XfvdduPffuPlcY+tlLr\niizm72mdC+n9MwOq/c6qfFeAat/VgX87rjidsQ0ww/ZMAEnnAbsCdxXq7EqeMrX9T0krSVoDeFNZ\n2ybpzHYFzrM9F3hQ0gxgG0kzgVG2J+d6ZwG7Ac2MvmMl/RK4ipR7nCzjRc1uuKc4fZL2JfmeeoyG\n0JchW1ZbbXlOOqnPQggGwYCy2mplixZbI0K2NKGmZIMgCCoqubWAhwvnj5AMwZ7qrNVi20bj/aNw\n/mgum5vb14/RjH1JOTKXZOH0roGmRl8z8nT1ocB7bL/aU/2+ZJVVluPQQ0sWvQVB0BLdZvQt0WkB\ngiDoLuYyouWjj+gxaPEA8Tbbb7O9j+3P5GO/in2eDWwIPCzpBUln9oGcA0JfTI31JSFPc0Ke/qED\n+rASHTP6JB0u6U5Jt0k6W1J3rYYMgmHKPJYsPZ6ZdAcPHvubBUcDHiWtCa6xNosv+noUWKdBnVba\nNhqvrK9G5c24XtKmPdRpl9G2l7a9HHAk3bAYKjPY/mmHPM0JefqHZvqw/hgMdEQKSaOBzwFvtv2a\npPNJke57HfogCIKB4bUmu9WWHbcNy45bOOP6yHGL/UlPBsZkHTCL9Hc/vq7OZcCBwPmSxgLP2X5C\n0lMttIVFPYOXAWdLOpk0fTsGuNG2JT0vaZss094sjENaxlhgqqQHSGv6RMWNbLaL6SyXB57qbV9B\nEAw8zfThYKRTpucLwGvA8pLmA8uRduMFQTDIqbKGxfY8SQeRghLXwq5Ml7R/uuyf2Z4oaWdJ95JC\ntnymWVsoT2dme5qk35Ji7M0BDsg7dyEZlr9iYciWP/Ugfr9saJN0PMnonA00jmMSBMGgpNvW9Gmh\n/hvggaXPAd8jKborbe/VoI5T2t/eMKaKeMAuvW96fLWRxx95eq/bnjPjs9UGP6xa80pU2L07YcNf\nVhr63BMqLM06qtLQJGdUbymNjtoCX8J2W+vlJPktLgm10YA7tE3bYww1WgnZkusdRpr9aJiKLenD\nIAj6i+GgDzs1vbs+8AVgNPA8cKGkCT3FqAqCoPMMlrUp3UJPIVsKnANMbNLPsDaeg2Aw0m36sFPS\nvg24zvYzAJIuIiUeb2D0FcNmbUB1D14QDFfuBe6r3Eu3TWcMZiSNsV1z1+5Geb6EIAgGId2mDztl\n9N0NHC1pGdKC6O1Ji6kb0CQfYBAEbTCGRV+aruxVL92m5AY5J0raiJTB9X7gfzosTxAEbVBVH/Y2\nLaWkkcDfgKXzcantI3oaryNGn+1bJZ0F3ExSdrcAP+uELEEQtMdgiTc1FLBdMZVMEASdpIo+rJKW\n0varkt5re7akEcB1kra1fV2zMTsWp8/2d2xvZnvzHOy0a+JTBcFw5jVGtnwE1ZC0u6Q7JM2TtHWh\nfLSk2ZKm5OPUTsqTrx0uaYak6ZI+OBDyNJDvGEmPFJ7LgKcQlbSjpLsk3ZM353QcSQ9KulXSLZJa\n33nQd+P/UtITkm4rlL1e0pWS7pZ0haSVOixPr747FfXhgrSU2QaqpZYsskhaSqCWlhLbs3OdkSR7\n7tme5I2MHEEQtMU8RrR8BJW5HfgocE2Da/fa3jofB3RSHkmbAJ8ENiFNQ53aIA/yQPG9wnPpKQxP\nn1Lw3OwAbAaMl/TmgZShhPnAONtb2e4pdWF/cAaLr9X6KnCV7Y2BvwKHd1ge6MV3p6I+LEs52axO\nLZUkkpaQdAvwODDJ9rSe5O2ubSdBEHScefPDmBsobN8NUGJADbhR1USeXYHzbM8FHpQ0g+TF+OcA\niwidTdu3wHMDIKnmubmraav+R3R2Zu/vOah6kV2B7fLnM4FJJEOwU/JAL747zfThnEnXM+eaf5Re\nr4rt+cBWklYErpS0ne1GL4gLCKMvCIK2mDs3jL5BwnqSppDCXh1t++8dlGUtoPjfbYE3ogMcJGkv\n4CbgS7afH8CxG3luOuFZq8fAnyXNA35m++edFghY3fYTALYfl7R6pwWiF9+dZvpQ73o3S7/r3QvO\nX/n6yfVVqqSlXIDtFyT9gRQZJYy+IAj6jnlzQ230Ja0Gb67jMWBd28/mtXWXSNq0Lq3bQMozYDST\nDzgV+HpOs3c8KQFAxYj1Q4Jtbc+StBrJ+Jve4ZeERnQ6+HivvjsV9WGVtJSrAnNsPy9pWeADwHE9\nDdgF2ru3+zueqTjujN43vXDDSiOfO7ZCdojtKw3N7hddWK2DClxI7zcynvuXCs8sDV6BCt8VoNp3\ndeD3P80LT1+f0kbw5mKbOeRF27anSLoP2AiY0gl5aMEb0Ve0Id/PgYE2Ulvx3Aw4tmfln09Kupjk\nfey00feEpDWyAfMG4F+dFMb2k4XTlr87VfRhlbSUwJrAmXmpxRLAr23/pacxu8DoC4JgMBFGX8dY\nsN4ov+U/Y3t+znA0hhTnryPykLwRZ0s6mTTFOQboxC7RN9h+PJ9+DLhjgEVoxXMzoEhaDljC9kuS\nlp2xD+wAABVRSURBVAc+SAseof4QhcW/M/uSEnDuA1zaSXl6+92pqg/zhpGN68pOqzs/qEG724Gt\n68t7Ioy+IAja4rV/RyiWgULSbsAPgVWByyVNtb0T8B7g65JeI+3M3N/2c52Sx/Y0Sb8FppHczwe4\nM4ndT5K0JemZPAjsP5CDl3luBlKGBqwBXKyUu3lJ4GzbvYvM3ksknQOMA1aR9BBwDHAicIGk/YCZ\npN3fnZTnvb357nSbPgyjLwiC9ghP34Bh+xLgkgblFwEXDRZ58rVvAd8aWIkWk2HvTo6fZVjMc9NJ\nbD8AbNlhGSaUXHr/gAqSKZHnjF511mX6MIy+IAjao8uUXBAEQb/RZfowjL4gCNpjbifDoAVBEAwi\nukwfRkaOIAjaY24bRwNaSVMl6ZSc0mtqXmfTtG1ZSidJE3LqqSn55zxJm+drk3JfteurVn84QRAM\nKyrqw4EmjL4gCNqjgpJrJU2VCgnGSYupf9pC24YpnWyfk1NPbQ3sBdxvu5Zv08D42nXbT1V8MkEQ\nDDfC6AuCYEhTTclVSTDerO2upFRO5J+7NRh7fG5TJHRgEAS9p8uMvljTFwRBe/y7UutW0lSVJSFv\n1naNFlI6fQrYpa7sV5LmABfZPr6dGwmCIKioDwecMPqCIGiPZm+st0yCqZP6esTerJReJEacpG2A\nl21PKxRPyKmplgcukrSn7d9UETQY+kj6CLCJ7ZM6LUswCBgkHrxWCaMvCIL2aKbk3jouHTV+tVjg\n/yoJxpdu0vbxHlI67QGcWywopKZ6OQdr3QYIoy9oSs4/3HJ6N0nqUKDqYCDoMqMv1rMEQdAe1daw\nLEhTJWlpkjF2WV2dy4C9AYoJxntoW0vpBHUpnXJuyk9SWM8naYSkVfLnpYAPM/Apu4JBRv5uTZd0\nRt4JfrakD0i6Lp+/XdI+kn6Y668u6aK8y/wWSWNzH3dJOlPS7cDaksZLui0fJxbGe1HS8bn99ZJW\ny+WfkHR77nNSZ55G0BJdtqYvjL4gCNqjgpKzPQ+opam6EzivlmBc0n/nOhOBB3KC8dOAA5q1zV1/\nG/iApLuB7Ukpnmq8B3jI9oOFspHAFZKmAlNI6wN/3utnEgwlNgC+k3eCbwzsYXtb4FDgCNLSgZrn\n7hRgku0tSXlQ78zlY4Af2X4r6S/hRFLary2Bt0uqrS1dHrg+t78W+FwuPxr4oO2tWHwdajCY6DKj\nL6Z3gyBoj4rKq7cJxsva5vJnKEnpZPsa4J11ZbOBt7UleDBceKCw9vPO/9/e3QfbVZV3HP/+AmRE\nkQhYoE1MQuRNM1LMaECliKA1Usd0ppby0irGaZkCg60vw4sdqB1mLDidVocWpcYY5F1gSoYixACp\nxZlgJAlQCBCEBIgk2AKCiBDD0z/2OuHk5J5zz747++y77vl9Zva49757nWed683i2XutvRawLO3f\nD8zsuPZYiqmASF24L0raG9gQESvTNe8F7kx/o0i6kuJGZAnwarrJAbiH1/+G7wIWp/WMB77cnpUw\nTpK5fmWQ9I31N/pCxbgV1shes3e10F/cZ8xFr/7kgkqhrz6yWvlKVlQoe33F2Gv+r0LhquupV/lb\nbaDFyayRMyvplbb919qOX2PH/2Z2G6v3Usdxt5eRtrTtb219fkScLum9FMMO7pE0JyKeG63i1oDM\n2kN375pZOb8psZnlp8zb4reThh9ImiRpzxE+4yfA0ZL2lrQLxXyRy3tWQJoVESsj4gKKl5Le1ut6\na1Bm7WEGT/rMbFzJ7M7WrKTosj/S8d8Al0n6LMW/jL8GNrVfl+aNPIfXE73/jIibu3xey9ckHZT2\nl7WtImPjTWbtYa1Jn6SFFI+nN0dEa73LvYBrgRnAeuCEiPhlnfUws50os0bOrF8RsQE4rO14QZef\ntVaMeYaRV385rP0gIq6l+O9eZ7w92/ZvAG5I+38y5i9hg5VZe1h39+4iinUy2424RqaZZSKzt9XM\nzGqTWXtYa9IXEXcBnYNP+1kj08zGq8waOTOz2lRsDyXNS/M6PiLp7C7XfEPSujSf4+Hp3DRJd0h6\nIM3peFY/1W1iTN++fayRaWbjlZM5M7NChfZQ0iTgEoq5RX8OrJR0U0Q81HbNx4C3R8RBko4Avgkc\nmSJ/PiLWSNqD4i3vpe1lRzIeXuTw8jRmOXm56QqYmY0T1drDucC6NF4USddQ9Ia2J27zeX0M6d2S\npqQlJzdRvDRERPxK0lpgakfZHTSR9G0eZY3MDre37R8AzKqxamYT2WPA49U/Zmv1jzAzmxCqtYdT\ngSfbjp+iSAR7XbMxndvcOiFpJsVqL3ePFnAQSZ/Yfs6i1hqZF9GxRubIjqupWmbDZhbb3zTdObaP\ncfeumVmh4fYwde1eD3wuIn412vV1T9lyFcV6g/tIegK4gGINwu9LWgBsoFgI3cxy4aTPzKzQqz1c\nvxw2LO9VeiMwve14WjrXec3bRrpG0q4UCd/3ImKUB2iFWpO+iDi5y49GXCPTzDLgpM/MrNCrPZx2\nTLG1/OgrnVesBA6UNAN4GjiRYsWWdkuAM4BrJR0JPN96GRb4DvBgRHy93+qOhxc5zCwnW0a/xMxs\nKFRoDyNiq6QzgaUUU+gtjIi1kk4rfhyXRcQtko6X9CjFms6nAkj6AHAKcL+k1RQvxZ4XEbf2iumk\nz8zK8YscZmaFiu1hStIO6Tj3rY7jM0co92Ngl7Lx6l6Rw8wmmooLjI91MtJeZSXtJWmppIcl3SZp\nSjo/Q9KvJa1K27+1lZkj6b70Wf9S7ZdiZkOpYns4aE76zKycLSW2Dm2TkX4UmA2cJOnQjmu2TUYK\nnEYxGeloZXst7/hoRMxJ2+lt5y8FPhsRBwMHS+pcMtLMrLcK7WETMujeHetv6sWKcTtfoCnjx9VC\nr3lHhbIHVYv9lmrFK3m+SuF1FYOvrVC2yt8KVPtbbaAlqdadMebJSCkm6uxWdj7wwVR+MbCcIhGE\n7aeMIpXdH3hzRKxMpy6nWBLytkrfzsyGS2bDXfykz8zKqbbW5EiTkU7t85peZfdrX94RaF/ecWbq\n2r1T0lFtMZ4apR5mZr1lthZ5Bk/6zGxcGXzjtcOTuj60lnd8GpgeEc9JmgP8h6R37ryqmdlQGyfJ\nXL+c9JlZOb16lJ9ZDr9Y3qt0lclIJ/cou2mk5R0j4lXg1bS/StLPgIN7xDAz6984GavXL3fvmlk5\nW3ts+xwDh/7969uOtk1GKmkyxWSkSzquWQJ8CqBjMtJeZVvLO0Lb8o6S3ppeAEHSLOBA4LHUBfxL\nSXMlKcXra0Z7M7NterWHnds44Cd9ZlZOhakHxjgZ6Wd6lU0ffRFw3QjLOx4N/IOkV4HXgNMiovXK\n0BnAd4E3ALeMNqmpmdkOxslULP1y0mdm5VTszhjrZKTdyqbzzzLC8o4RcSNwY5fPugd4V98VNzPr\nlFn3rpM+MytnnHRTmJk1LrP20EmfmZWT2dtqZma1yaw9dNJnZuVk1siZmdUms/bQSZ+ZlZPZGBYz\ns9pk1h466TOzcjIbw2JmVpvM2kMnfWZWTmZTFJiZ1Saz9tBJn5mVk1l3hplZbTJrD530mVk5mXVn\nmJnVJrP2cAInfVXT72czjb129Et6eX73auUreblC2Sq/M4AXKpR9sWLszG4VM3tbzcysNpm1hxM4\n6TOzWmTWyJmZ1Saz9nBS0xUws8xsKbGZmU1kFdtDSfMkPSTpEUlnd7nmG5LWSVoj6d1t5xdK2izp\nvn6rW2vSN1KFJF0saW2q/A2S9qyzDma2k20tsZmZTWQV2kNJk4BLgI8Cs4GTJB3acc3HgLdHxEHA\nacClbT9elMr2re4nfSNVaCkwOyIOB9YB59ZcBzPbmX5TYjMzm8iqtYdzgXURsSEitgDXAPM7rpkP\nXA4QEXcDUyTtl47vAp4rU91ak76RKhQRyyLitXS4AphWZx3MbCdz966ZWaFaezgVeLLt+Kl0rtc1\nG0e4pm9Nj+lbAPyg4TqYWRkVu3fHMIbl8NHKStpL0lJJD0u6TdKUdP7Dkn4q6V5JKyV9qK3Mnemz\nVktaJemt1X4xZjZ0Mhvu0tjbu5K+DGyJiKuaqoOZjUGMvWjbGJbjgJ8DKyXdFBEPtV2zbQyLpCOA\nbwJHjlL2HGBZRFycksFz07lfAB+PiE2SZgO3sX3vwkkRsXrs38jMhlrP9nB52rraCExvO56WznVe\n87ZRrulbI0mfpFOB44FjR796edv+zLSZWXnr09aobWNYACS1xrA81HbNdmNYJLXGsBzQo+x84IOp\n/GKKhuOciLi39aER8YCkN0jaLY2fgeZ7O8xswjombS1f6bxgJXCgpBnA08CJwEkd1ywBzgCulXQk\n8HxEbG77udLWl0EkfdtVSNI84EvA0RHxyujFj6mrXmZDZibb3zT9VxOVGGkMy9w+rpk6Stn9Wg1h\neqq3b2dgSZ8EVrUlfADflbQFuDEiLhzD9zEzG5OI2CrpTIoXXCcBCyNiraTTih/HZRFxi6TjJT0K\nvAR8plVe0lUUSdI+kp4ALoiIRb1i1pr0jVQh4DxgMvBDSQArIuL0OuthZjvTwN/Q6Psuts12nS6p\na/erwEfaTp8cEU9LehNwo6Q/j4grKtTTzIZOtfYwIm4FDuk4962O4zO7lD25bLxak74uFeqZhZrZ\neNdrCvofpa2rKmNYJvcou0nSfhGxWdL+wDOtiyRNA24E/iIi1rfOR8TT6X9fSjeocwEnfWZWQl5L\ncng8i5mV1GtOgvcBZ7dtO9g2hkXSZIoxLEs6rlkCfAqgYwxLr7JLgFPT/qeBm1L5twA3A2dHxIpW\nAEm7SNon7e8GfBz4nzH8MsxsqOU1h5XX3jWzksZ+Z1tlDEu3sumjLwKuk7QA2ACckM6fAbwdOF/S\nBRTdvn8I/Bq4TdKuwC7AMuDfx/zFzGxI5fWkTxEV5l+omaQohgE2YbcKZd9cMXaVlen2rhh794rl\nq3i5QtlnK8Z+oULZFyvGbuoO8CtERKnxcsW/yU0lSuxfOoaZWQ5ybA/9pM/MShof3RRmZs3Lqz10\n0tdVlf8jqz75afKJV5N/ElUek1f9h9dk7Nzk1Z1hZlafvNpDJ31mVlKVmxIzs4kkr/bQSZ+ZlTRs\nTzbNzLrJqz100mdmJeXVnWFmVp+82kMnfWZWUl53tmZm9cmrPXTSZ2Yl5XVna2ZWn7zaQyd9ZlZS\nXne2Zmb1yas9dNJnZiXldWdrZlafvNpDJ31mVlJeUxSYmdUnr/bQSZ+ZlZRXd4aZWX3yag+d9JlZ\nSXk1cmZm9cmrPZzUdAWqWT+EsX/WUFyAdQ3GfrTB2I81GHt9g7G7+W2JzcxsIsurPXTSl13sJhOQ\nJhOvJpPdxxuMvb7B2N1sKbGZmU1kebWH7t41s5LGxx2rmVnz8moPnfSZWUnj447VzKx5ebWHioim\n69CVpPFbObMJICJU5vri3+Q/lSjxhdIxzMxykGN7OK6f9DX9yzGzkeTVnWFmVp+82sNxnfSZ2XiU\nV3eGmVl98moPx3X3rpmNL5LWAzNKFNkQETPrqY2ZWXNybA+znLJF0jxJD0l6RNLZA4w7TdIdkh6Q\ndL+kswYVu60OkyStkrRkwHGnSPq+pLXp+x8xwNjnppj3SbpS0uQaYy2UtFnSfW3n9pK0VNLDkm6T\nNGWAsS9Ov/M1km6QtGcdsfsVETMjQiW2mU3W18ysLjm2h9klfZImAZcAHwVmAydJOnRA4X8LfD4i\nZgPvA84YYOyWzwEPDjgmwNeBWyLiHcDvA2sHEVTSDOAvgXdHxGEUQxJOrDHkIoq/rXbnAMsi4hDg\nDuDcAcZeCsyOiMMpZseuK7aZmU1w2SV9wFxgXURsiIgtwDXA/EEEjohNEbEm7f+KIvGZOojYUDxp\nBI4Hvj2omCnunsAfRMQigIj4bUS8MKDwLwCvAm+StCvwRuDndQWLiLuA5zpOzwcWp/3FwB8PKnZE\nLIuI19LhCmBaHbHNzGziyzHpmwo82Xb8FANMvFokzQQOB+4eYNh/Br4EDHog5gHA/0palLqWL5O0\n+yACR8RzFO/EPwFsBJ6PiGWDiN1m34jYnOqzCdh3wPFbFgA/aCi2mZllLsekr3GS9gCuBz6XnvgN\nIuYfAZvTk0albVB2BeYA/xoRc4BfU3R51k7SLOBvKQbL/h6wh6STBxG7h4G//STpy8CWiLhq0LHN\nzGxiyDHp2whMbzuels4NROpivB74XkTcNKi4wAeAT0h6DLga+JCkywcU+yngyYj4aTq+niIJHIT3\nAD+OiGcjYitwI/D+AcVu2SxpPwBJ+wPPDDK4pFMpuvWbTnbNzCxjOSZ9K4EDJc1Ib3GeCAzyTdbv\nAA9GxNcHGJOIOC8ipkfELIrvfEdEfGpAsTcDT0o6OJ06jsG9TPIwcKSkN0hSil33SySdT1KXAKem\n/U8DdSb728WWNI+iS/8TEfFKjXHNzGyCy25y5ojYKulMircaJwELI2JQb5J+ADgFuF/SaopuvvMi\n4tZBxG/YWcCVknYDHgM+M4igEXFveqJ5D7AVWA1cVlc8SVcBxwD7SHoCuAD4R+D7khYAG4ATBhj7\nPGAy8MMi52VFRJxeR3wzM5vYPDmzmZmZ2RDIsXvXzMzMzEpy0mdmZmY2BJz0mZmZmQ0BJ31mZmZm\nQ8BJn5mZmdkQcNJnZmZmNgSc9Fnf0pq7h46x7OOS9t7ZdTIzM7P+ZDc5szUnIv6qSvGdVhEzMzMr\nzU/6bAdpibu1kq6Q9KCk6yTtLulOSXMkTZf0iKS9VfiRpA+nsqdIulvSKkmXpqXTIC0tJumNkm6W\ntFrSfZL+tLEvamZmNkSc9Fk3hwCXRMQ7gReA00lP6yLiCYqlyb4JfAF4ICKWpa7fPwPeHxFzgNco\nlq1rNw/YGBHvjojDgGFYws7MzKxxTvqsmyciYkXavxI4qv2HEfEdYE/gNOCL6fRxwBxgZVqb+Fjg\ngI7PvR/4iKSvSjoqIl6s6wuYmZnZ6zymz/q13Zg8SbsD09LhHsBLFF24iyPiy10/JGKdpDnA8cCF\nkpZFxIU11dnMzMwSP+mzbqZLOiLtnwz8N2lcXnIRcAVwPvDtdO524JOSfgdA0l6Sprd/qKTfBV6O\niKuAr1E8GTQzM7OaOemzbh4GzpD0IDAFuJT0tE/S0cB7gIsi4mrgFUmfjoi1wN8BSyXdCywF9k+f\n13pS+C7gJ6n793zAT/nMzMwGQBGeScO2J2kGcHNEvKvpupiZmdnO4Sd91o3vBszMzCYQP+kzMzMz\nGwJ+0mdmZmY2BJz0mZmZmQ0BJ31mZmZmQ8BJn5mZmdkQcNJnZmZmNgSc9JmZmZkNgf8HwCkY0aFt\nMQEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb2b0fede10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "ax1.set_title('Power of beam profile')\n",
    "im1 = ax1.imshow(beam_profile, interpolation='nearest')\n",
    "ax1.set_xlabel('pixels')\n",
    "ax1.set_ylabel('pixels')\n",
    "fig.colorbar(im1, ax=ax1)\n",
    "fig.subplots_adjust(right=1.5)\n",
    "ax2.set_title('Plotting power around\\n{0}x{1} pixel ROI\\n laser spread\\n at specimen plane (uW)'.format(w, h))\n",
    "im2 = ax2.imshow(roi,\n",
    "                 extent=[-x_um, x_um, -y_um, y_um],\n",
    "                 interpolation='nearest')\n",
    "ax2.set_xlabel('microns')\n",
    "ax2.set_ylabel('microns')\n",
    "fig.colorbar(im2, ax=ax2)\n",
    "r_x = np.array([-w_um/2, +w_um/2, +w_um/2, -w_um/2, -w_um/2])\n",
    "r_y = np.array([-h_um/2, -h_um/2, +h_um/2, +h_um/2, -h_um/2])\n",
    "ax2.plot(r_x, r_y, 'b-', label='ROI', linewidth=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb28efb9a10>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEoCAYAAACtnQ32AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeclPW59/HPtRRlwaUtZQVcml2xxAJRYNETFX0S4pOm\nyYntsSQRWx4bxhzxmORokseoMccSsZ0cYyzJiYk1UdaCBTWIRFBBEZAqRQFRQLieP373rrO70xZ2\n5r5n5/t+vea1M3e9ZmZnrvmV+/czd0dERCSdirgDEBGR5FKSEBGRjJQkREQkIyUJERHJSElCREQy\nUpIQEZGMlCREyoiZDTKztWZmRT5vXzN7xsw+MrNfFPPcsn2UJKTNmdm3zexlM1tnZovN7GEzO6wN\njnuFmd3dFjGWCzObb2ZHNDx290XuXuXFv0DqTGCFu3d394uKfG7ZDkoS0qbM7IfAtcBPgL7ALsBv\ngC/HGReAmXWIY18BoBaYHXcQsg3cXTfd2uQGVAHrgP+dZRsDLgXmAR8A9wI9onW1wFbgJGABsAK4\nLFp3NLAxuq0DZqSc8zZgCbAIuAqwaN3JwHOEpLUS+Pc08VQAl0XxfAS8DAyI1m0FfgC8DbwTLfsi\nMB1YA7wEjEo51inAO8Da6O+J0fJhQD3wYfScfp/htdkB+K8o1obj98n1PKP1ZxC+hNcC/wT2B+4G\ntgAfR8svTHmNK6L9aoA/A6ui53l6yjGvAP4A3BXtPws4MMt7m/a1Ae4ANkXv3VrgiLj/V3Vrxec6\n7gB0az+36It8U8MXUIZtzgOej76cOgE3AfdE6xq+wG4BOgMjgE+B3aP1VwB3Nzven4D/BHYEqoEX\ngTOidScDm6Mv+gpghzTxXATMBIZHj/cFekb3twKPA92jL/CewGrg29HxToge9wQqCUmm4Tj9gD2j\n+/cAk6L7nYEvZnhtzoy+sHcgJNMDgG55PM9vRInjwOjxUGBQdH8+MC7lHLVR4mhIEs8Av47ei/0I\nSawu5fXeEL2vBvwMeCFD7Blfm2j9HaRJ0rol/xZ7ALq1n1v0BbEkxzazm31p1TQklpQvsJqU9S8B\n34zuN0kShOqsT1O//KMvp6ei+ycD7+WI503gf2VYtxUYm/L4X4EXm23zPKHkUxl9KR4P7Nhsm7uA\nm4lKKFliOZVQ8tm32fJMz/PJ6P5jwDkZjjmflF/uqUkCGERIopUp638G3J7yej+Rsm5P4OMM58n4\n2kT3lSRK9KY2CWlLq4BqM8v2f1UL/MnMVpvZakLS2Ez45d1gecr9DUC3LMfqBCyNjreG8GVcnbLN\nohwxDwLezbL+/ZT7OxOqwVItIHz5bwC+BXw/iucvZrZ7tM1FhC/l6WY2y8xOzXCuuwkll3vN7H0z\nuzpqC8n0PPukPId3cjzPdGqA1VHsTZ5PyuNlKfc3ADtmeH8zvjbbEJckiJKEtKUXCPXOX82yzUJg\nvLv3im493b2ruy/N4/jNe+QsIvzC7p1yrB7uPiLLPuniGZbnOZcAg5ut3wVYDODuf3P3o4D+wFvA\nb6PlK9z9THcfAHwP+E8zG9riRO5b3P0qd9+bUL//ZUIpJdfzXJTlOWR7/kuAXmbWNd3zaaWsr42U\nLiUJaTPuvpZQRfEbM5tgZl3MrKOZHWNmV0eb3QL8zMx2ATCzPmb2lZTDZOu/vxwY3NDH392XAU8A\nvzKznSwYamZjWhH2FOAqMxsexbOvmfXMsO0jwK5mdoKZdTCzbxGqYP4aXQfwFTOrJJSM1hOqdTCz\nr5tZwy/qDwnVWFubH9zM6sxsn+iX+vroOFvyeJ63ARea2YHRcYaZ2aCU16x5Qmp4/d4nVAn9h5nt\nYGYjgP9DaDzPJNP7k/G1yXIsKQFKEtKm3P1a4IfA5YRG0IXA2cD/RJtcT2icfcLMPiJ8SR2Seojm\nh0y5fz/hS2qVmb0SLTuZ0Bg8m9AmcD/hl3y+rgXuS4nnNqBLuljcfTXwvwi9hFZGf4+LlldEz3tx\ntG4MoeoJ4GDgJTNbS3gdznX399LE0h94gNAA/gYwFfhdtO6kTM/T3R8AfgrcE53jT0CvaL//AH4c\nVVP9MM3zOhEYQigJPAj82N2nZnm90pZMcrw2GfeT5GvoKigiItKCShIiIpKRkoSIiGSkJCEiIhkp\nSYiISEZKEmWi+Wig8jkzG2tmuS66KxnR6LuD444jGzM72cyejTsOyU1JQiTYpm5+SRy+3N13ytDF\nNmnUtbIEKElIQRRyaG0N2y1SPEoSZcjMDjaz581sTTQp0K/NrGPK+l+Z2fJoFrGZZrZXtLyzmf3S\nzBaY2VIz+08z2yFaN9bMFpnZxWa2FLg9zXmHmVm9mX1oZivM7Pcp67aa2Tlm9k607ucp6042s+fM\n7FozW0m4qhszO83MZpvZKjN7tOEq7mjddWa2MHoOL5vZ4SnrdjSzO6MLzP5JuNgt2+uV9lhmdjRh\nmPFvRVU8MzLsf4mFsZjWmtkcMxsXLb/CzO43s3ujda9EVz037FdjZg9Er8c7ZnZOyroKM7vMzOal\nxDUg5bUcGt2/w8x+Y2aPRDE+Y2b9o+e0Jnr99svznFeY2R/M7K4o3lkWXeWd4XlnfE/zfY3zOW+2\nmKUNxD3CoG7FuZEyGihwIOEqZyOMr/MG4SpggKMIcyrsFD3eHegX3f8V4Yrh7kBXwpXTP43WjSUM\nI/EzwmB06YblzjhkNmGYiiejYw8kjH10WrSuxZDfwATC/Ae78fmcENNSjvdtoEe07gJgKdA5Wnc1\n8HR0rgGEeRIWZnntsh2rxfDlzfbdjXDVecNruAswJGXfjYSRYzsA/5cw2GCH6L15BfhR9HgwYc6L\nL0X7ZhvifAswNLp/B+HK9/2j1/xJ4D3gO9E5ruLzUXNznTPvocPzfE+facVrnPa8uWLWrQ2+O+IO\nQLcivdHNhoxutu484MHo/jjC8NmHkjKpTbRufcMXXPR4FPBudH8sYRC6TlliyDhkdvSF8qWUx98H\n/hbdbzHkN2GsoFNTHlcQJtcZlOHcq4mG4CaMmJp6rjPIkiRyHCtXkhhGGEn1SKBjs3VXAM+nPDbC\nsB6HEZJ48+d8KTAlup9riPPUJHFLyrqJwBspj/chjARL9J5nO2feQ4fn+Z4+k2Xf5q9x2vPmilm3\n7b+puqkMmdmuFoayXmpmHxLG/akG8DBuz42EKUeXm9nNZtbNzPoQ5kx41T4f5vtRoHfKoT9w981Z\nTp1ryOzUYbkXEIafbtC891EtcH1KLKsIDaENVS4XRlUpaywMrV3F50OI75zmXBnlOFZW7v4OcD4w\nmfB63mNmqWNLLUrZ1glJYufo+Q1oeH7ReScR5paA3EOcp0odev2TNI8bhmLfJcc5If+hwxtke08b\n5fEaZzpvPjHLdlCSKE83AXOAYe7eg1BUbxzd091vdPeDgL0I1U0XEQZt2wDs7Z8P893D3bunHDdr\nbxXPPWT2oJT7uxAGnct07IXAWd50yPFu7v5iVJ99EfD1aHlPwrSZDc9xabNz1WaKOY9j5eyh4+73\nuvvolPNck+45m5kRqmUapih9t9nz6+7uDXOF5xrifFvkOue2yPaeAmBmo8n+Ghc7ZkmhJFGedgLW\nuvsGM9uDz0crxcwOMrNDLDRkf0KoQtoa/cr9LXBdVKrAzAaY2VH5ntRyD5l9kZn1sDDM9XmE+a8z\nuQW4zD5vVO9uZl9PeX6bCaPFdjazf4uWNbgPmBSdayChCiaTXMdqMnx5mue8m5mNM7POhBn4Pmn2\nnL9gZl+10GPrAsLr/SJhruh1FjoC7Ghh+O29zeygaL/WDHGeS0Psuc6Zbd9M8nlPu5H9NW7rmKUV\nlCTKR+ov3guB71gYVvoWmn5wqwjJYDWhHWMl8Ito3SWERsEXo2qqJwgNs/nKNWT2n4FXgX8AfyFN\nD6nGJ+P+P4QG6HujWF4HjolWPx7d3o6ewwaaVlddSfglPp8w9We26xxyHSvd8OWpdoji/IDwK7oP\noTok9Tl/C1hDaEw+3sPkQ1sJQ2/vH513BeF9qYr2y3uI8zw4QB7nzLhvFvm8p7le47aOWVqh4EOF\nm9kxwHWEhDTF3a9Js80NwHhCw+Mp7v5atLw74Z9/H8Kvr9Pc/aWCBiyxMLOthJ46+dazlzwzu4JQ\n5XdS3LEUQjm+p+1RQUsSUcPSjYSua3sDJ0bVG6nbjCd8UHYFziL0fmlwPfCIu+8J7EeoRxcRkSIp\ndHXTIcBcd18Q9Xq5l9C/PdUEouJ+VErobmb9zKwKGO3ud0TrPvMwPaa0Txqiof3Re9oOdMy9yXYZ\nQNO6xfdpOlVlum0WR8u2ACvN7A5CKeIV4Dx3/6Rw4Upc3L3shtpw9yvjjqGQyvE9bY+S3HDdkXBl\n8G/c/UBCY9al8YYkIlJeCl2SWEzoG91gYLSs+TaDMmyzyN0beow8QOhd04KZqVgrItJK7p7zWpRC\nlyReBoabWW3UT/wE4KFm2zwEnARgZiOBD919ubsvBxaZWUMXyyOB2ZlOFPel67luV1xxRZsfE0oj\nzlJ5PRVn/HHEEWe5fo7yVdCShLtvMbOJhP70DV1g55jZWWG13+ruj5jZsWY2j9AFNnWohnOB/zaz\nToQhCJoP4yAiIgVU6Oom3P0xwtAOqctuafY47RWv7j6THMM4i4hI4SS54bpdqaurizuEvCjOtqU4\n25biLL6CX3FdDGbm7eF5tJYZlOHTFmlT5fo5MjM8AQ3XIiJSwpQkREQkIyUJERHJSElCREQyUpIQ\nEZGMlCRERCQjJQkREclISUJERDJSkhARkYyUJEREJCMlCRERyUhJQkREMlKSEBGRjJQkREQko4JP\nOiQJsnEjzJ4Nn3wSdyQiba9rV9hrL+jUKe5I2hUlifZuyxa4/nq47z6YNQuGDoWddoo7KpG299FH\nsGAB7LcffPvb8P3vQ4UqS7aXJh0qYTknS3n3XTjpJOjcGa64Ag46KPzaEmmv1q2D6dPhxz+Gykq4\n4w4YNCjrLpp0KDul2fZq2jQ49FD42tfg73+HsWOVIKT922knOPJIeOYZOOII+MIXYMaMuKMqaSpJ\nlLCMv4C2boVDDoELL4QTTih6XCKJceutcM89MHVq+MCkoZJEdipJtEd/+EOoi/3Wt+KORCRep50G\nH3wADz8cdyQlSyWJEpb2F9DGjbDHHnDnnaGKSaTc/fWvcMklMHMmdGzZV0cliexUkmhvfvMb2Hdf\nJQiRBscdB336hEZsaTWVJEpYi19An3wCtbVQXx/6i4tI8PLL8NWvhi6yzUoTKklkp5JEe/LCCzB8\nuBKESHMHHww9e8Jrr8UdSclRkmhPpk6FcePijkIkmcaNC58RaRUlifZk6lSoq4s7CpFkqqtTktgG\nBW+TMLNjgOsICWmKu1+TZpsbgPHAx8Cp7j4jWv4e8BGwFdjs7odkOIfaJD7+GPr1g+XLddGcSDor\nV8KwYbBqVZN2CbVJZFfQkoSZVQA3AkcDewMnmtkezbYZDwxz912Bs4CbUlZvBerc/YBMCUIizz8P\n+++vBCGSSXV16Njx6qtxR1JSCl3ddAgw190XuPtm4F5gQrNtJgB3A7j7S0B3M+sXrbMixNg+1Ner\nPUIkl3HjwmdF8lboL+ABwKKUx+9Hy7JtszhlGwf+ZmYvm9kZBYuyPVB7hEhuapdotaT/Sj/M3Q8E\njgXONrPD4w4okdavh9dfh1Gj4o5EJNnGjg1Vs5s3xx1JySj0fBKLgV1SHg+MljXfZlC6bdx9afT3\nAzP7E6H66rl0J5o8eXLj/bq6OurK6Vf1tGlw4IFhaGQRyaxXrzCnyiuvlN2Pqvr6euq3oaqtoL2b\nzKwD8BZwJLAUmA6c6O5zUrY5Fjjb3Y8zs5HAde4+0swqgQp3X29mXYEngCvd/Yk05ynv3k2TJoU5\nI668Mu6QRJLvggugb9/wuUG9m3JtV9DqJnffAkwkfMG/Adzr7nPM7CwzOzPa5hFgvpnNA24BfhDt\n3g94zsxmAC8Cf0mXIAR49lkYMybuKERKw5gx4TMjedHYTSWs8RdQnz6hTaKmJu6QRJLvrbfg2GPh\nnXcAlSRybtcevlzLOkms+RAGDgzTNmaYVEVEUmzaFGawW7cOOndWksgh6b2bJJd33gmD+ilBiOSn\nc2cYMCCMCCs5KUmUunnzQpIQkfwNHx4+O5KTkkSpU5IQaT0libwpSZQ6JQmR1lOSyJuSRKlTkhBp\nPSWJvClJlDolCZHWU5LIm7rAljAz8C6VYeymCuV7kbx9+in06AHr12OdOqoLbBb6Zil1Q4cqQYi0\n1o47hqE5Fi3KvW2Z07dLqVNVk8i2UZVTXpQkSp2ShMi2UZLIi5JEqdt117gjEClNu+6qJJEHJYlS\np5KEyLZRSSIvShKlTklCZNsoSeRFXWBL1SefYJVd8M+2QIcOcUcjUno+/hiqq7FPP1EX2CxUkihV\n774b/ipBiGybrl3DdKaSlZJEqVIxWWT7qbo2JyWJUqUkIbL9lCRyUpIoVbpSVGT77bJL3BEknpJE\nqVq6NO4IREpf//5xR5B4ShKlatmyuCMQKX01NXFHkHhKEqVKJQmR7aeSRE5KEqVKJQmR7aeSRE5K\nEqXo449h8+a4oxApff36hb9bt8YbR4IpSZSiZcv0C0ikLXTuHP6uWhVvHAmmJFGKli1TXapIW1L1\nbUZKEqVo6VKVJETakjqCZKQkUYpUkhBpWypJZKQkUYpUkhBpWypJZFTwJGFmx5jZm2b2tpldkmGb\nG8xsrpm9Zmb7N1tXYWb/MLOHCh1ryVBJQqRtqSSRUUGThJlVADcCRwN7Ayea2R7NthkPDHP3XYGz\ngJubHeY8YHYh4yw5KkmItC2VJDIqdEniEGCuuy9w983AvcCEZttMAO4GcPeXgO5m1g/AzAYCxwK3\nFTjO0qKShEjbUkkio0IniQFA6nCl70fLsm2zOGWbXwEXAWU4b1QWS5cqSYi0JZUkMuoYdwCZmNlx\nwHJ3f83M6oCs0+xNnjy58X5dXR11dXWFDC8+W7bAypXQt2/ckYi0H2VQkqivr6e+vr7V+xV0jmsz\nGwlMdvdjoseXAu7u16RsczMw1d3/ED1+ExhLaIv4V+AzoAuwE/BHdz8pzXnKZ47rZctgxAhYsQIz\nynJuXpG2ZAa+Y5dw1XVlZdzhFE1S5rh+GRhuZrVm1hk4AWjeS+kh4CRoTCofuvtyd7/M3Xdx96HR\nfk+lSxBlR0NyiLS9/v3LojSxLQpa3eTuW8xsIvAEISFNcfc5ZnZWWO23uvsjZnasmc0DPgZOLWRM\nJU/tESJtr3//8NkaOjTuSBKn4G0S7v4YsHuzZbc0ezwxxzGeBp5u++hKkEoSIm2vpkYliQx0xXWp\nUUlCpO01lCSkBSWJUqOShEjbU0kiIyWJUqOShEjbU0kiIyWJUqOShEjbU0kio1YnCTPraWYjChGM\n5EElCZG2p5JERnklCTOrN7MqM+sF/AP4rZldW9jQJC2VJETankoSGeVbkuju7muB/w3c7e6HAv9S\nuLAkrXXrwiXW3brFHYlI+9K3L3zwQRj2RprIN0l0NLMa4JvAXwsYj2TTUIqwnFfSi0hrdOoEPXuG\ncdGkiXyTxJXA48A8d3/ZzIYCcwsXlqS1bBn06xd3FCLtU79+qnJKI98rrpe6e2Njtbu/qzaJGKxc\nCX36xB2FSPvUp49KEmnkW5L4dZ7LpJCUJEQKR0kirawlCTMbBXwR6GNmP0xZVQV0KGRgksYHH0B1\nddxRiLRP1dXhMyZN5Kpu6gx0i7bbKWX5WuDrhQpKMli5EgY0n9hPRNpEdbVKEmlkTRINo6+a2Sfu\n/vPUdWb2DdR4XVwrV8J++8UdhUj7VF0Nb70VdxSJk2+bxAlplk1qy0AkDytXqrpJpFBUkkgrV5vE\neOBYYICZ3ZCyqoowragUk5KESOEoSaSVq01iCfAK8BXg1ZTl64ALChWUZKDeTSKFo95NaeVqk5gJ\nzDSze9x9c5FikkzUu0mkcNS7Ka18L6YbbGb/AewF7Niw0N01IWyxfPopbNwIO+2Ue1sRab3evUNJ\nwl1D36TIt+H6DuAmQjvEOOBu4HeFCkrSWLUq/NLRP69IYXTpEsZwWr8+7kgSJd8k0cXdnwTM3Re4\n+2TguMKFJS2o0Vqk8NR43UK+SWKjmVUAc81sopkdT7jITopFSUKk8JQkWsg3SZwHVALnAl8A/hU4\nuVBBSRrq2SRSeOrh1EKu6yT2B2a6+8vRovXAqQWPSlpSzyaRwlMPpxZy9W66DRhqZq8CzwPTgBfc\nfV3BI5OmVN0kUniqbmoha3WTux8EDAR+CmwkVDfNM7OZZvafRYhPGihJiBSekkQLOdsk3H2Du9cD\n1wO/An4DdAWOKWxo0oSShEjhKUm0kKtN4tuE+ST2J5QkXgZeAg53d83zV0xKEiKFpyTRQq42iVuA\nt4CbgWfc/e3WnsDMjgGuI5Raprj7NWm2uQEYD3wMnOLur5nZDsAzhDktOgN/dvfLWnv+dkO9m0QK\nT72bWsiVJHoA+xFKE5PNbHdgKfACoQH7qWw7R9dW3AgcSRgs8GUz+7O7v5myzXhgmLvvamaHEhLS\nSHffaGbj3H2DmXUAppnZYe4+bRufa2lT7yaRwlPvphZyNVxvcfd/uPuN7v5twrDhjxG6wf4tj+Mf\nAsyNrtLeDNwLTGi2zQTCMB+4+0tAdzPrFz3eEG2zQxTrmvyeVjvjHn7d9O4ddyQi7Zuqm1rI1SYx\nglCKaLh1JnSF/TWhO2wuA4BFKY/fJySObNssjpYtj0oirwLDgJvdfXYe52x/1q8PY8p06RJ3JCLt\nW69esGYNbNkCHTrEHU0i5KpuuhN4DngUuNzdFxY8ohTuvhU4wMyqgCfMbGw0pWoLkydPbrxfV1dH\nXV1dUWIsCjVaixRHx45QVQUfftjuSu719fXU19e3er9c80kc2HyZmR3o7v/I8/iLgV1SHg+MljXf\nZlC2bdx9rZk9DBwE5EwS7Y6ShEjxNFQ5tbMk0fzH85VXXpnXfvmO3ZTqtlZs+zIw3MxqzawzYa7s\nh5pt8xBwEoCZjQQ+dPflZlZtZt2j5V2ALwGvbUO8pU89m0SKRz2cmsh30qFUeU9o4O5bzGwi8ASf\nd4GdY2ZnhdV+q7s/YmbHmtk8QhfYhrGhaoC7zMyiff8rGq68/KgkIVI8arxuYluSRH5llIi7Pwbs\n3mzZLc0eT0yz3yygRXVXWVL3V5HiUTfYJvJOEmb2FWBMdH+Lu/+lYFFJUypJiBSPShJN5NUmEc1v\nfR4wO7qda2Y/K2RgkkJJQqR4lCSayLckcRywf9QlFTO7C5gBlO8wGcWkJCFSPNXV8MYbcUeRGK3p\n3dQj5X73tg5EslDvJpHiUe+mJvItSfwHMMPMphJ6N40BLi1YVNKUShIixaPqpiZyJomoC+pzwEjg\n4GjxJRoqvIjUu0mkeNS7qQlz99wbmc1y932LEM82MTPP53mUpC1bYIcd4NNPw5ABKczC2H8isu1a\nfI4+/BBqa+Gjj2KLqRjMDHfPed1bvm0S/zCzg3NvJm3uww/DWDIdt+WSFhFpte7dYcMG2LQp7kgS\nId9vnkOB75jZAsJV0Ua4YnpEwSKTQO0RIsVlFsZtWrUKamrijiZ2+SaJowsahWSmJCFSfA2N10oS\neVc31QCro8mDFhAm/+lfuLCkkbq/ihSfusE2yjdJ3ASsT3m8PlomhaaeTSLFpx5OjfJNEk26D0VX\nXqsltRhU3SRSfLpWolG+SeJdMzvXzDpFt/OAdwsZmESUJESKT0miUb5J4nuEOa4XE+apPhQ4s1BB\nSQolCZHiU5JolFeVkbuvIMwqJ8WmJCFSfNXV8NJLcUeRCFmThJld7O4/N7NfAy2u7XX3cwsWmQTq\n3SRSfOrd1ChXSWJO9PcV0iQJKQL1bhIpPvVuapQ1SaTMPjebMHfE4JR9HLi7YJFJoOomkeJTm0Sj\nfAf4ewu4CJgFbG1YHl1YF7t2O8Dfpk3QtWv4ay3H4dIAfyLbL+3naMOGMDTHhg1pP3vtQb4D/OV7\nrcNKd39oO2OS1lq1KvyjttN/UpHEqqwMn7sNG8IPtTKWb5KYbGZTgL8DGxsWuvsfCxKVBKpqEolP\nQ5WTkkReTgF2j7ZvqG5yQEmikNSzSSQ+DT2camvjjiRW+SaJg9x994JGIi2pJCESHzVeA/lfcf28\nme1V0EikJXV/FYmPusEC+ZckRgKvmdl8QpuEJh0qBpUkROKjkgSQf5I4pqBRSHorV8KwYXFHIVKe\nlCSA/MduSsT1EGVn5Uo49NC4oxApT9XVMHNm3FHELt82iW1mZseY2Ztm9raZXZJhmxvMbK6ZvWZm\n+0fLBprZU2b2hpnNMrPyGydK1U0i8VFJAihwkjCzCuBGwhzZewMnmtkezbYZDwxz912Bs4Cbo1Wf\nAT90972BUcDZzfdt99QFViQ+GuQPKHxJ4hBgbjQ39mbgXmBCs20mEI0B5e4vAd3NrJ+7L3P316Ll\n6wmDDQ4ocLzJot5NIvFR7yag8EliALAo5fH7tPyib77N4ubbmNlgYH+gfAZ4d1d1k0icVN0ElMA8\n1WbWDXgAOC8qUaQ1efLkxvt1dXXU1dUVPLaCahhYrLIy7khEylPv3rB6NWzdChUFb74tuPr6eurr\n61u9X16jwG4rMxsJTHb3Y6LHlxKur7gmZZubganu/ofo8ZvAWHdfbmYdgb8Cj7r79VnO0/5GgV2w\nAEaPhoULM26iUWBFtl/Wz1GPHjB/PvTsWdSYiiHfUWALnR5fBoabWa2ZdSZMgdp8NNmHgJOgMal8\n6O7Lo3W3A7OzJYh2S1VNIvFTlVNhk4S7bwEmAk8AbwD3uvscMzvLzM6MtnkEmG9m84BbgO8DmNlh\nwHeAI8xshpn9w8zK56I+9WwSiZ96OBW+TcLdHyOMIJu67JZmjyem2W8a0KGw0SWYejaJxE89nAp/\nMZ1sI1U3icRP1U1KEomlJCESPyUJJYnEUpIQiZ+ShJJEYqlNQiR+apNQkkis5cuhf/+4oxApb/37\nh89iGVOSSKqlS5UkROLWv3/4LJYxJYkkcodly6CmJu5IRMpbTU34LJYxJYkkWrcujBXQrVvckYiU\nt759Q8OCljvPAAAUGklEQVT1li1xRxIbJYkkUilCJBk6doRevcq68VpJIonUHiGSHGXeLqEkkUQq\nSYgkR5m3SyhJJNGyZSpJiCRF//5KEpIwqm4SSQ5VN0niqLpJJDlU3SSJo5KESHKoJCGJo5KESHKo\nJCGJo5KESHKUeUnCPOMM4KXDzLw9PA8ANm+Gykr49FPokH1ivqwTuItIXnJ+jtatC6WJ9euLFlMx\nmBnubrm2U0kiaVasCMMT50gQIlIk3bqFLLJuXdyRxEJJImnUHiGSLGZl3S6hJJE0ao8QSZ4ybpdQ\nkkgalSREkkclCUkMDckhkjxlPDSHkkTSLF2qkoRI0tTUqLpJEkIlCZHkUUlCEkMN1yLJo4ZrSQw1\nXIskjxquJRHcVZIQSSKVJArHzI4xszfN7G0zuyTDNjeY2Vwze83MDkhZPsXMlpvZ64WOMxHWroVO\nnaBr17gjEZFUffvC6tXw2WdxR1J0BU0SZlYB3AgcDewNnGhmezTbZjwwzN13Bc4CbkpZfUe0b3lQ\nKUIkmTp0gN69w7A5ZabQJYlDgLnuvsDdNwP3AhOabTMBuBvA3V8CuptZv+jxc8CaAseYHGqPEEmu\nMm2XKHSSGAAsSnn8frQs2zaL02xTHlSSEEmuMm2XUMN1krz/PgwcGHcUIpLOwIGwaFHu7dqZjgU+\n/mJgl5THA6NlzbcZlGObnCZPntx4v66ujrq6utYeIn7z5sH++8cdhYikM3w4vPNO3FFss/r6eurr\n61u9X0EnHTKzDsBbwJHAUmA6cKK7z0nZ5ljgbHc/zsxGAte5+8iU9YOBv7j7vlnO0z4mHTrySLj0\nUvjSl/LaXJMOiWy/vD9HDz4Iv/sd/OlPBY+pGBIx6ZC7bwEmAk8AbwD3uvscMzvLzM6MtnkEmG9m\n84BbgB807G9m9wDPA7uZ2UIzO7WQ8cZu3rzwa0VEkmf48PAZLTOavjQpPv0UevQIUyR2zK8WUCUJ\nke2X9+do/fpwvcT69VBR+s25iShJSCu8+y7U1uadIESkyLp1g+7dYcmSuCMpKiWJpFBVk0jy7bpr\n2VU5KUkkhZKESPKVYbuEkkRSKEmIJF8ZJglVgCfFvHnw5S/HHUUTgwcPZsGCBXGHIWWutraW9957\nL+4wguHD4b774o6iqJQkkiKBJYkFCxZQ8r3GpOSZ5eyAUzxlWJJQF9gk2LQJqqpg3bowVHieCt0F\nNuoiV7gTiOSh0P+HrfocffQRDBgQPqtJSl7bQF1gS8l774VxYVqRIEQkBt27Q2UlLF8edyRFoySR\nBAmsahKRDMqsyklJIgmUJERKh5KEFN3cuUoSbei5555jzz333KZ9Fy1aRFVVVVHaYk499VR69erF\nyJEjc267YMECKioq2Lp1a8HjkhyGDw+f2TKhJJEEKklskyFDhvDUU0+1WH744YczZ86cNHvkNmjQ\nINauXdvYo2bcuHHcfvvt2xVnOs899xxPPvkkS5Ys4cUXX8xrn0T18okMHjyYyspKqqqqqKmp4bvf\n/S7r1q1rss3zzz/PkUceSVVVFT179mTChAlN3p+nn36aQYMGNT90cpXZVddKEkmgJFF23nvvPQYP\nHsyOO+7YJseLqxeamfHwww+zdu1aZs6cyaxZs/jJT37SuP6FF17g6KOP5vjjj2fp0qXMnz+fESNG\ncNhhhzW59iGJCTAjVTdJUW3eDAsXwpAhcUfSbjT/ZTpkyBB++ctfMmLECKqqqjj99NNZsWIFxx57\nLN27d+eoo47io48+AppW61x++eU8++yzTJw4kaqqKs4991wALrjgAvr160f37t3Zb7/9mD17dto4\nli5dyoQJE+jduze77bYbt912GwC33347Z5xxBi+88AJVVVVceeWVLfbdunUrF154IX369GH48OE8\n/PDDTdaPGzeOyy+/nMMPP5yuXbsyf/587rzzTvbaay+qqqoYPnw4t956a+P2dXV1/CmaB2HatGlU\nVFTw6KOPAvDUU09xwAEHAHDXXXcxevRoLrroInr16sWwYcN47LHHsr7eDQmqb9++HH300bzxxhuN\n6y655BJOOeUUJk6cSNeuXenRowdXXXUVI0eObDJRWElpSBLl0j3c3Uv+Fp5GiZozx33w4G3atdBP\nO+mv6+DBg/3JJ59ssby+vt4HDRrUZLtRo0b5Bx984EuWLPG+ffv6gQce6DNnzvSNGzf6EUcc4f/+\n7//u7u7vvfeeV1RU+JYtW9zdva6uzqdMmdJ4rMcff9wPOuggX7t2rbu7v/nmm75s2bK08Y0ePdon\nTpzomzZt8tdee8379OnjU6dOdXf3O++800ePHp3xud10002+5557+uLFi33NmjU+bty4FnHV1tb6\nnDlzfMuWLb5582Z/5JFHfP78+e7u/swzz3hlZaXPmDHD3d3/7d/+zc8991x3d//Zz37mw4cP90sv\nvbRx3fnnn98YV+fOnX3KlCm+detWv+mmm3znnXfO6z1YtGiR77vvvo2v5YYNG7xDhw5eX1/fYr87\n7rij8bjN36/mCv1/uE2Hr652X7y4zWMppuh1zfn9qpJE3J55Bg4/PO4otplZ29wK7ZxzzqG6upqa\nmhpGjx7NyJEjGTFiBJ07d+b4449nxowZeR2nU6dOrFu3jtmzZ+Pu7L777vTr16/Fdu+//z4vvPAC\n11xzDZ06dWK//fbj9NNP5+67787rPPfffz/nn38+O++8Mz169GDSpEkttjnllFPYY489qKiooGPH\njowfP57BgwcDMHr0aI466iieffZZAMaOHcvTTz8NwDPPPMOkSZMaHz/99NOMHTu28bi1tbWcdtpp\nmBknn3wyy5YtY8WKFRlj/epXv0pVVRW77LILw4YN40c/+hEAq1evZuvWrdTU1LTYp6amhpUrV+b1\nWiTSF78YPrtlQEkibvX1UIrzcUfC77DtvxVa6hd5ly5dWjxev359XscZN24cEydO5Oyzz6Zfv358\n73vfS7vvkiVL6NWrF5WVlY3LamtrWbw4v+nblyxZ0qTKrLa2tsU2zRt7H330UUaNGkXv3r3p2bMn\njz76aOMX8ahRo3j77bdZsWIFM2fO5KSTTmLRokWsWrWK6dOnM2bMmMbj9O/fv/F+ly5dcPesr8+f\n//xn1q5dS319PVOnTuXVV18FoGfPnlRUVLB06dIW+yxdupTq6uq8XotEqqsLn90yoCQRJ3eYOhXG\njYs7EskgXYPqxIkTeeWVV5g9ezZvvfUWv/jFL1pss/POO7N69Wo+/vjjxmULFy5kwIABeZ23pqaG\nRYsWNT5ON9BiamybNm3i61//OhdffDEffPABa9asYfz48Y3tBV26dOELX/gC119/Pfvssw8dO3Zk\n1KhRXHvttQwfPpxevXrlFVc6DecYM2YMEydO5OKLLwagsrKSUaNGcf/997fY57777uNf/uVftvmc\nsRs3Lnx2y4CSRJzefhs6d1aj9XbYtGkTGzdubLxt2bJlu4/pKUWbfv368e677zY+fuWVV5g+fTqf\nffYZXbp0Yccdd6QizVSWAwcO5Itf/CKTJk1i48aNvP7660yZMoXvfve7ecXwzW9+kxtuuIHFixez\nZs0arrnmmqzbb9q0iU2bNlFdXd3YKP3EE0802WbMmDHceOONjVVLdXV1TR63hfPPP5/p06czffp0\nAK6++mruuusubrzxRtavX8+aNWu4/PLLefHFF7niiiva7LxFN2IErFxZFrPUKUnEaerUUGwtpe5/\nCXPcccdRWVlJly5dqKysTNtTqHlpIFd3y9T15513Hvfffz+9e/fm/PPPZ+3atZxxxhn06tWLIUOG\nUF1dzUUXXZT2OL///e+ZP38+O++8M1/72te46qqrGJdnqfGMM87g6KOPZr/99uOggw7ia1/7Wtbn\n0K1bN2644Qa+8Y1v0KtXL+69914mTJjQZJuxY8eyfv36xqqlhse5kkS216v5uurqak455RSuvvpq\nAA477DAef/xxHnzwQWpqahgyZAgzZ85k2rRpDB06NPuLkGQVFTBmTFlUOWkU2DidcAIccwyccso2\n7a5RYKUcJGoU2FQ33AD//CekdDUuJRoFNuncS77RWqSs1dWVRbuEkkRc3nwTunSBqMuiiJSYffaB\nNWvg/ffjjqSglCTi0tAeISKlqaICxo5t9+0SShJxqa9X11eRUjduXLtPEmq4jsPKlbDHHjBzZpgK\ncRup4VrKQWIbriEMGX744eFvVVWbxlVoarhOsquughNP3K4EISIJsOuuoYdimgsq2wuVJIpt3jwY\nORJmz4a+fbfrUIUuSQwePDjtlb4ixVRbW9tkWPG2tt2fo4UL4YADYNYs2HnnNour0PItSRQ8SZjZ\nMcB1hFLLFHdvcemomd0AjAc+Bk5x99fy3TfarnSSxLe+Fa7WjAZB2x6FThIi5aBNPkeXXAKrV8Nv\nf9smMRVDIqqbzKwCuBE4GtgbONHM9mi2zXhgmLvvCpwF3JzvvqWkvr4eXnoJpk2DCy6IO5yM6kuk\nEU5xti3FuZ0mTYKHHgoX15HgOLdBodskDgHmuvsCd98M3AtMaLbNBOBuAHd/CehuZv3y3Ldk1N90\nE0yYEOouU0YGTZpS+edWnG1LcW6nHj3gpz8N7ROPPZbcOLdBoZPEAGBRyuP3o2X5bJPPvsm3di2c\neio88QQ8+GBosBaR9uf00+Guu+DMM+HhhyFlBOBSlsTeTe1rtLt//jOM9Pr978Nhh8UdjYgU0pFH\nwuuvh2mJ586NO5o2UdCGazMbCUx292Oix5cSpsy7JmWbm4Gp7v6H6PGbwFhgSK59U46h5lsRkVbK\np+G6Y4FjeBkYbma1wFLgBKB5fctDwNnAH6Kk8qG7LzezlXnsC+T3REVEpPUKmiTcfYuZTQSe4PNu\nrHPM7Kyw2m9190fM7Fgzm0foAntqtn0LGa+IiDTVLi6mExGRwkhiw3WrmdnBZjbdzGZEfw+KO6ZM\nzOwcM5tjZrPM7Oq448nGzP6vmW01s22fALmAzOzn0Wv5mpk9aGaJGTzHzI4xszfN7G0zuyTueNIx\ns4Fm9pSZvRH9P54bd0zZmFmFmf3DzB6KO5ZMzKy7md0f/V++YWaHxh1TOmY2KYrvdTP7bzPrnGnb\ndpEkgJ8Dl7v7AcAVQCIHUjGzOuDLwL7uvi/wy3gjyszMBgJfApI8LscTwN7uvj8wF5gUczxASV0I\n+hnwQ3ffGxgFnJ3QOBucB8yOO4gcrgcecfc9gf2AxFWRR+28ZwAHuPsIQrPDCZm2by9JYinQPbrf\nA1gcYyzZfB+42t0/A3D3lTHHk82vgPSTNyeEu//d3bdGD18EBsYZT4qSuBDU3Zc1DIHj7usJX2iJ\nvBYp+tFyLHBb3LFkEpVkR7v7HQDu/pm7r405rHTWApuArmbWEagElmTauL0kiUuBa81sIaFUkYhf\nlGnsBowxsxfNbGpSq8XM7CvAInefFXcsrXAa8GjcQURK7kJQMxsM7A+8FG8kGTX8aElyI+oQYKWZ\n3RFVi91qZl3iDqo5d18D/D9gIeEH9Yfu/vdM2xe6C2ybMbO/Af1SFxH+YS4HzgHOcff/MbOvA7cT\nqkqKLkecHYGe7j7SzA4G7gOGFj/KnHFeRtPXL7Yuxlni/JG7/yXa5kfAZne/J4YQS56ZdQMeAM6L\nShSJYmbHAcvd/bWoyjapXd47AgcCZ7v7K2Z2HeEH7BXxhtWUmQ0FLgBqgY+AB8zs25k+PyWTJNw9\n45e+mf2uYb27P2BmU4oXWVM54vwe8Mdou5ejRuHe7r6qaAFGMsVpZvsAg4GZZmaEKpxXzewQd19R\nxBCB7K8ngJmdQqiGOKIoAeVnMbBLyuOBJLQKNKpueAD4L3f/c9zxZHAY8BUzOxboAuxkZne7+0kx\nx9Xc+4QS+CvR4weAJHZaOAiY5u6rAczsj8AXgbRJor1UN801s7EAZnYk8HbM8WTyP0RfZma2G9Ap\njgSRjbv/0937u/tQdx9C+Mc/II4EkUs0lPxFwFfcfWPc8aRovIg06jVyAuGi0SS6HZjt7tfHHUgm\n7n6Zu+/i7kMJr+VTCUwQuPtyYFH02QY4kmQ2tL8FjDSzHaMfgkeSpYG9ZEoSOZwF/Cb6QH4KnBlz\nPJncAdxuZrOAjUDi/tHTcJJbvP810Bn4W/hf50V3/0G8IZXOhaBmdhjwHWCWmc0gvNeXuftj8UZW\n0s4F/tvMOgHvEl0cnCTuPtPM7gZeBbYAM4BbM22vi+lERCSj9lLdJCIiBaAkISIiGSlJiIhIRkoS\nIiKSkZKEiIhkpCQhIiIZKUmIbAMz+7KZXRx3HCKFpuskRIrAzMz1YZMSpJKESDPRcBpzotE834om\nZfmSmU2LHh9sZieb2a+j7fua2R+jyY9mmNnI6Bhvmtld0RX2A83sxGiSl9dTJ5wys3Vm9pNo/+fN\nrE+0/BvRZEAzzKw+nldDyp2ShEh6w4BfuPvuwO7ACe5+GGGsqMsIQ1g0lAxuAOqjyY8OBN6Ilg8H\nbowmmPoMuBqoIwzJfXA0JDtAV+D5aP9nCRPCAPwYOCqaTKthW5GiUpIQSW++uzcMzvYG0DDe/izC\nKLmpjgBuAvBgXbR8gbu/HN0/GJjq7qujiZL+GxgTrdvk7o9E919NOf5zwF1mdjrtZ5w1KTFKEiLp\npY4quzXl8VZafmFnamv4uNnjTAMlbk65v6Xh+NFghT8CBhGGa++ZI2aRNqckIZJea0a+fRL4AYT5\nraNpLJsfYzphVsJeZtYBOBGozxqA2VB3f9ndrwBWEJKFSFEpSYik5xnup3t8PjDOzF4HXgH2bL6d\nuy8jzFJWTxia+RV3/2uG4zX4RUNDN2GSmNdb/SxEtpO6wIqISEYqSYiISEZKEiIikpGShIiIZKQk\nISIiGSlJiIhIRkoSIiKSkZKEiIhkpCQhIiIZ/X+4N3Sdfi0qfAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb28f31dd10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Center cross section of\\nlaser spread at specimen plane')\n",
    "center_line = roi[:, y/2]\n",
    "fig = plt.plot(np.linspace(-x_um/2, x_um/2, center_line.shape[0]), center_line, 'r-')\n",
    "plt.xlabel('microns')\n",
    "plt.ylabel('micro-Watts')\n",
    "center_line = roi[:, y/2]\n",
    "# Draw ROI boundaries\n",
    "plt.axvline(x=-w_um/2, label='Limits of drawn ROI')\n",
    "plt.axvline(x=w_um/2)\n",
    "plt.legend(loc='lower center')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

raster.R

## Calculate power distribution of laser scanning photo-activation device
##
## Laser beam symbols used are from the beam power formula:
## [1] en.wikipedia.org/wiki/Gaussian_beam#Power_through_an_aperture
##
## References:
## [2] radiantzemax.com/kb-en/Goto50125.aspx
## [4] en.wikipedia.org/wiki/Beam_diameter#D4.CF.83_or_second_moment_width
## [5] en.wikipedia.org/wiki/Normal_distribution
##
## Theory contributors in alphabetical order:
## Browne, M. (Andor Technology)
## Karunarathne, A. (Washington University in St. Louis)
## Magidson, V. (Wadsworth Institute, NYSDOH)
##
## Code Author:
## Pariksheet Nanda (Andor Technology)

library(spatstat) # calulate image convolution
library(ggplot2)  # display complex plots
library(reshape2) # reshaping matrix as data.frame for ggplot

## Create the round laser beam profile at the specimen plane

# Gaussian laser beam power (micro-Watts) at back focal plane, Symbol: Ps
# Measured by user's PD-300 laser meter using single pixel FRAP
Ps <- 4.8 / 40

To <- 0.84 # Percent transmission through above objective at 445nm 
Tm <- 0.8 # Percent transmission through media, etc (estimated)
# Gaussian laser beam power (micro-Watts), Symbol: Po
Po <- Ps * To * Tm

# FWHM Beam diameter / waist size (microns) at objective, Symbol: Wo
# Measured at objective specimen plane using FRAPPA slide
Wo_FWHM <- 0.8

# Convert FWHM to 1/e^2 [2]
Wo <- Wo_FWHM * 0.8493218 * 2

camera_pixels <- 512 # e2v CCD97 pixel
fov_width <- 76.6 # microns
# Pixel calibration (microns / pixel)
cal <- fov_width / camera_pixels

# Beam axis pixel length on camera (pixels), Symbol: px_len
# Wo for a single mode beam is 4 sigma wide [4]
px_len_2sig <- Wo / cal
# Extend the beam diameter out to 6 sigma to contain 99% of Po
px_len <- (6 / 4.) * px_len_2sig

px_fit <- floor(px_len) / px_len
px_to_edge <- floor(3 / px_fit)
edge <- px_fit * px_to_edge
steps <- floor(px_len / px_fit)
#x <- seq(-edge, edge, length=steps)  # keep for debugging 1D construction

# Map 1D arrays radially
x <- seq(-edge, edge, length=steps)
mx <- matrix(rep(x, length(x)), length(x), length(x)) # square matrix of x
my <- t(mx)
z <- sqrt(mx^2 + my^2)

# Create 2D gaussian beam profile
norm <- function(x, mean=0, sigma=1) {
    # Returns single value from Gaussian / Normal distribution curve [5]
    (1/(sigma * sqrt(2 * pi))) * exp(-0.5 * ((x - mean) / sigma) ^ 2)
}

beam_profile <- norm(z)

# norm() gives a 1D probability distribution, the area of which is 1.  The 2D
# beam profile thus needs to be re-normalized back to 1, for it to be
# multiplied by the Po scalar, obj_mag^2, and transmission losses
beam_profile <- beam_profile / sum(beam_profile)
beam_profile <- beam_profile * Po

## Setup ROI

w_microns <- 10
w <- round(w_microns / cal) # pixels
h <- 1 # pixels

length <- ncol(beam_profile)  # length in pixels of beam_profile square
x <- w + length - length %% 2
y <- h + length - length %% 2
roi <- matrix(1, x, y)

## Apply the laser beam over the ROI

roi <- convolve.im(as.im(roi), as.im(beam_profile))

## Show results

# Area calculation
w_um <- w * cal
h_um <- h * cal
x_um <- x * cal
y_um <- y * cal
actual_area <- x_um * y_um
drawn_area <- w_um * h_um

# Power calculation
peak_pixel_power <- max(roi) # uW
actual_power <- sum(roi) # uW
drawn_power <- sum(roi[length/2:w+length/2,length/2:h+length/2])
density_drawn_power <- drawn_power / (drawn_area / 1000**2) # uW/mm^2

# Energy calculation
dwell_time <- 80  # micro-seconds
repeats <- 1
peak_pixel_energy <- max(roi) * dwell_time * repeats / 1000 # uJ
actual_energy <- actual_power * dwell_time * repeats / 1000 # uJ
drawn_energy <- drawn_power * dwell_time * repeats / 1000 # uJ
density_drawn_energy <- drawn_energy / (drawn_area / 1000**2) # uJ/mm^2

sprintf('Laser Beam at Back Focal Plane-')
sprintf('Measured Power: %3.0f uW', Ps)
sprintf('')
sprintf('Laser Beam at Objective Specimen Plane-')
sprintf('Calculated Power: %3.0f uW', Po)
sprintf('Measured Width: %1.1f um, in FWHM', Wo_FWHM)
sprintf('')
sprintf('Beam Profile Pixelation at Back Focal Plane-')
sprintf('Calibration: %1.3f um/pixel', cal)
sprintf('Size: %d x %d pixels, covering 3 standard deviations',
        nrow(beam_profile), ncol(beam_profile))

m <- as.matrix(beam_profile)
mt <- melt(t(m))
names(mt)[names(mt)=="Var1"] <- "y"
names(mt)[names(mt)=="Var2"] <- "x"
mt$x <- (mt$x - mean(mt$x)) * cal
mt$y <- (mt$y - mean(mt$y)) * cal
ggplot(mt, aes(x=x, y=y, fill=value)) +
  scale_fill_distiller(palette="Spectral") +
  geom_raster(hjust=1, vjust=1) +
  coord_fixed(expand=FALSE) +
  labs(list(
    title=expression(paste('Beam profile power of laser (', mu, 'W)')),
    x=expression(paste(mu, 'm')),
    y=expression(paste(mu, 'm'))
  )) +
  geom_blank()

# 2D plot of laser at specimen plane
m <- as.matrix(roi)
mt <- melt(t(m))
names(mt)[names(mt)=="Var1"] <- "y"
names(mt)[names(mt)=="Var2"] <- "x"
mt$x <- (mt$x - mean(mt$x)) * cal
mt$y <- (mt$y - mean(mt$y)) * cal
ggplot(mt, aes(x=x, y=y, fill=value)) +
  scale_fill_distiller(palette="Spectral") +
  geom_raster(hjust=1,vjust=1) +
  coord_fixed(expand=FALSE) +
  labs(list(
    title=expression(paste('Laser power at specimen plane (', mu, 'W)')),
    x=expression(paste('Specimen plane (', mu, 'm)')),
    y=expression(paste('Specimen plane (', mu, 'm)'))
    )) +
  geom_blank() +
  annotate("rect", color="black", alpha=0,
           xmin=-w_um/2, xmax=w_um/2,
           ymin=-h_um/2, ymax=h_um/2)

# 1D cross section of laser at specimen plane
m <- t(as.matrix(roi))
y <- m[nrow(m)/2,]
x <- 1:length(y)
x <- (x - mean(x)) * cal
line <- data.frame(x, y)
ggplot(line, aes(x=x, y=y)) +
  geom_line(color="red") +
  geom_vline(xintercept=c(-w_um/2,+w_um/2)) +
  labs(list(
    title="Cross section of laser spread at specimen plane",
    x=expression(paste('ROI length (', mu, 'm)')),
    y=expression(paste('Laser power (', mu, 'W)'))
  ))