larrybradley/ellipse_hl.ipynb

## ellipse_hl.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e258ce00-2ec0-4fc4-9f00-6aa2e32967c3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-01T02:32:29.744224Z",
     "iopub.status.busy": "2023-06-01T02:32:29.743777Z",
     "iopub.status.idle": "2023-06-01T02:32:29.754055Z",
     "shell.execute_reply": "2023-06-01T02:32:29.753304Z",
     "shell.execute_reply.started": "2023-06-01T02:32:29.744193Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from astropy.modeling.models import Gaussian2D\n",
    "from photutils.aperture import EllipticalAperture\n",
    "from scipy.optimize import root_scalar\n",
    "\n",
    "# create a simulated source\n",
    "xypos = (50, 50)\n",
    "theta = np.pi / 3\n",
    "m = Gaussian2D(100, xypos[0], xypos[1], 10, 5, theta)\n",
    "yy, xx = np.mgrid[0:101, 0:101]\n",
    "data = m(xx, yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7890b56a-d668-4e1d-8eaa-74003e3900c4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-01T02:32:30.426054Z",
     "iopub.status.busy": "2023-06-01T02:32:30.425718Z",
     "iopub.status.idle": "2023-06-01T02:32:30.431468Z",
     "shell.execute_reply": "2023-06-01T02:32:30.430678Z",
     "shell.execute_reply.started": "2023-06-01T02:32:30.426026Z"
    }
   },
   "outputs": [],
   "source": [
    "# define the objective function to minimize\n",
    "def ellip_opt_fcn(radius, data, aperture, ba_ratio, normflux):\n",
    "    aperture.a = radius\n",
    "    aperture.b = radius * ba_ratio\n",
    "    flux, _ = aperture.do_photometry(data)\n",
    "    return flux[0] - normflux\n",
    "\n",
    "\n",
    "# define the optimizer\n",
    "def ellip_opt(data, ellip_aper, normflux, opt_fcn, \n",
    "             min_radius=0.1, max_radius=100):\n",
    "    ba_ratio = ellip_aper.b / ellip_aper.a\n",
    "    args = (data, ellip_aper, ba_ratio, normflux)\n",
    "    result = root_scalar(opt_fcn, args=args, \n",
    "                         bracket=[min_radius, max_radius],\n",
    "                         method=None)\n",
    "    return result.root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5b48d85f-761c-4187-be6c-707605685514",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-01T02:32:30.976472Z",
     "iopub.status.busy": "2023-06-01T02:32:30.976109Z",
     "iopub.status.idle": "2023-06-01T02:32:30.998415Z",
     "shell.execute_reply": "2023-06-01T02:32:30.997571Z",
     "shell.execute_reply.started": "2023-06-01T02:32:30.976442Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11.785082078757133, 5.892541039378567)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define a starting ellipse (e.g., based on calculated \n",
    "# centroid and shape parameters);\n",
    "# xypos and theta will be fixed; a and b vary with fixed ratio\n",
    "aper = EllipticalAperture(xypos, a=10, b=5, theta=theta)\n",
    "ba_ratio = aper.b / aper.a\n",
    "\n",
    "target_flux1 = np.sum(data) / 2.0  # e.g., half the flux\n",
    "a_hl = ellip_opt(data, aper, target_flux1, ellip_opt_fcn)\n",
    "b_hl = a_hl * ba_ratio\n",
    "a_hl, b_hl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1271fb36-e3ce-48b2-8b25-149a220fb16e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-01T02:32:31.481122Z",
     "iopub.status.busy": "2023-06-01T02:32:31.480777Z",
     "iopub.status.idle": "2023-06-01T02:32:31.495481Z",
     "shell.execute_reply": "2023-06-01T02:32:31.494585Z",
     "shell.execute_reply.started": "2023-06-01T02:32:31.481093Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7.594084743091976, 3.797042371545988)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target_flux2 = np.sum(data) / 4.0  # e.g., quarter flux\n",
    "a_ql = ellip_opt(data, aper, target_flux2, ellip_opt_fcn)\n",
    "b_ql = a_ql * ba_ratio\n",
    "a_ql, b_ql"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6c88ab92-eb3f-4d38-9a09-9edf1e8cb8b9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-01T02:32:31.980797Z",
     "iopub.status.busy": "2023-06-01T02:32:31.980456Z",
     "iopub.status.idle": "2023-06-01T02:32:31.986237Z",
     "shell.execute_reply": "2023-06-01T02:32:31.985494Z",
     "shell.execute_reply.started": "2023-06-01T02:32:31.980770Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15707.962941461983 15707.962941461981 -1.8189894035458565e-12\n"
     ]
    }
   ],
   "source": [
    "# check the result\n",
    "aper_hl = EllipticalAperture(xypos, a=a_hl, b=b_hl, theta=theta)\n",
    "flux, _ = aper_hl.do_photometry(data)\n",
    "print(flux[0], target_flux1, target_flux1 - flux[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f61497eb-9e37-4d33-b106-2abd97eb4135",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-01T02:32:32.453254Z",
     "iopub.status.busy": "2023-06-01T02:32:32.452911Z",
     "iopub.status.idle": "2023-06-01T02:32:32.458644Z",
     "shell.execute_reply": "2023-06-01T02:32:32.457748Z",
     "shell.execute_reply.started": "2023-06-01T02:32:32.453226Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7853.981470730979 7853.981470730991 1.1823431123048067e-11\n"
     ]
    }
   ],
   "source": [
    "aper_ql = EllipticalAperture(xypos, a=a_ql, b=b_ql, theta=theta)\n",
    "flux, _ = aper_ql.do_photometry(data)\n",
    "print(flux[0], target_flux2, target_flux2 - flux[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5b1718fe-8bd5-48d6-b3c0-da267bb47caf",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-01T02:32:33.241002Z",
     "iopub.status.busy": "2023-06-01T02:32:33.240661Z",
     "iopub.status.idle": "2023-06-01T02:32:33.443287Z",
     "shell.execute_reply": "2023-06-01T02:32:33.442754Z",
     "shell.execute_reply.started": "2023-06-01T02:32:33.240974Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11fe14050>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAGgCAYAAAAtsfn1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABTRklEQVR4nO3dd3gUVdsG8HtbNn1DAqRAQiK9i4AYQQGJBEEEQRQFpAmKQUpsIMKLCILohwoqqK8CvoIoKqgoINIUDaE36T2UhBKSTd063x/B3ZnZZJPApk3u33XlujI7Z2bODmyenfOcohIEQQAREZECqCu6AkRERJ7CoEZERIrBoEZERIrBoEZERIrBoEZERIrBoEZERIrBoEZERIrBoEZERIrBoEZERIrBoEZERIpR6qD2xx9/oHfv3oiIiIBKpcLq1asl+wVBwLRp0xAeHg4fHx/ExcXhxIkTkjLp6ekYNGgQAgMDERQUhJEjRyI7O/u23ggREZG2tAfk5OSgdevWGDFiBPr16+eyf+7cuZg/fz6WLl2KmJgYTJ06FfHx8Th8+DC8vb0BAIMGDcLly5exYcMGWCwWDB8+HKNHj8by5ctLVAe73Y5Lly4hICAAKpWqtG+BiIgqmCAIyMrKQkREBNRqDzYaCrcBgLBq1SrHtt1uF8LCwoR33nnH8VpGRoag1+uFr7/+WhAEQTh8+LAAQNi5c6ejzNq1awWVSiVcvHixRNdNSUkRAPCHP/zhD3+q+E9KSsrthCEXpX5Sc+fMmTNITU1FXFyc4zWDwYAOHTogKSkJAwcORFJSEoKCgtCuXTtHmbi4OKjVaiQnJ+PRRx91Oa/JZILJZHJsCzcXFuiEntBC58m3QERE5cAKC7bhVwQEBHj0vB4NaqmpqQCA0NBQyeuhoaGOfampqahdu7a0ElotgoODHWXkZs+ejTfeeMPldS100KoY1IiIqpyCZxOPp5CqRO/HyZMnIzMz0/GTkpJS0VUiIqJKyKNBLSwsDACQlpYmeT0tLc2xLywsDFeuXJHst1qtSE9Pd5SR0+v1CAwMlPwQERHJeTSoxcTEICwsDBs3bnS8ZjQakZycjNjYWABAbGwsMjIysHv3bkeZTZs2wW63o0OHDp6sDhERVTOlzqllZ2fj5MmTju0zZ85g3759CA4ORlRUFCZMmICZM2eiYcOGji79ERER6Nu3LwCgadOm6NGjB0aNGoVFixbBYrFg7NixGDhwICIiIjz2xogqM5UK8A3yhU+AnsNSSHEEQUBelgm5Gbm42a+v3JQ6qO3atQtdu3Z1bCcmJgIAhg4diiVLluCVV15BTk4ORo8ejYyMDHTq1Anr1q1zjFEDgGXLlmHs2LHo1q0b1Go1+vfvj/nz53vg7RBVfoG1/NH9ufvRoH00NBoNwJhGSiMANpsNJ3eexW+L/oDxavlNrqEShPKOo7fPaDTCYDCgC/qw9yNVKRqtGmM+fxoRMaHw8/aHqmr01SIqNQF25ORn49KZNCwc+SVsVrtkv1WwYAt+RGZmpkf7SXi0Sz8RuRcUbkBgTX/4ewdAw48fKZoG/t4BCKyZg6CwQFy/kFEuV+XXRKJypFarbubQ2OZI1UHB/3e1pvxCDYMaEREpBts/iCoJrZcWGm35fM+0We2wmq3lci2i8sSgRlQJaL20iG4eWW7NNHabHWf/SamwwDbk2afQpFEzTHnxdQBAXn4eXpn2Ev7asQ05OTnYuWkPAgMC3R7zwCOd8fTAYRj21PASXfPCpQvo1qcLVn/1E5o2bubZN1SIxu0b4KN3FiKuy4Mu107evR1PPze40PdJt4dBjagS0GjVUGvUuHw6DeZ8c5ley8vbC+F3hEKjVcNaiktNmv4KjNlGfPzuIsnrnvgDvWrND9i1bydW/HclagTVQIB/8ZPcfrf0B/j4+N7S9Yryw8/f4615M7Fr816PnleuTau7sG1tUoneJ5UOgxpRJWLON8OUW7ZBrTJKuXAe9aPro1GDRiU+JrhGSBnWqGx56bxQq2atiq6GIrGjCBF51I2MG0icMgH39eyI1p1aoPfAnliz/uciyw959il8sexz7Ny7E43bN8CQZ58q0XUeeKQzlixf7Ng+dfYUnnzmCbTs2Aw9H4/H38l/oXH7Bvh9ywbJcSkXUzDkuUFo3akFHnnqYew9sAdAwRPn5BmvIis7C43bN0Dj9g2w4NMPirz+71s34NHBj6Blx2bo1qcrPvxsPqzWkjXnJu/ejsbtG8CYZQRQ8ITYrmsb/L5lA7r364aWHZth5AvDcDn1kuOYo8ePYMhzg9Cmc2vc1aU1+g3pg4OHDzr279q3C0+NGohWnZqjc69OmPnuDOTm5ZaoPkrCoEZEHmU2m9C8SQt8+t5nWLPiVzz+6EC88p+XcOCf/YWWXzD3Yzze9wm0adkG29YmYcHcj0t9TZvNhoSXxsDH2xsrF3+PGa/NwnsL5xVa9r2F8zBy8DNYvexnREdF48XXJ8JqtaJNq7vwWuLr8Pfzx7a1Sdi2NgkjBj9T6Dl27d2JV//zMp4eOAy/frMOM157Ez+s+QGLFpe+7v/Kz8/Hwi8+xtvT38HX//0WxqwsTJwywbH/pamJCKsdhu+W/oAfvvwRo4Y+C522oLHt/IVzGDVuBLp3jcdPy3/Be2/Nx+59u/DmXNclu5SOzY9EVGJbtm1Gm/tbSV6z2W2S7dDaYRg5xBkMhjzxNLZt/xNrN/yKVs1bu5wzyBAEb28f6HS6W26S+yt5G1IunMf/Fi1znGPimEQMHzvUpeyIwSPRpVPBVH/jRo9HrycewrkL51A/uj4C/AOgUqmKrceHny3A6KHP4tGH+wEAIutGYfyzE/DOgrkYO2rcLb0Hi9WCaa/8B61b3AkAmDN9LnoOiMeBf/ajVfPWuJR2CSOHjEL96PoAgOioaMexnyxZhN49HnF0momOisaUl6ZhyLNPYfqkGdDr9bdUp6qIQY2ISqxD23swfZL02//+Q/vx8rQXHds2mw2LFi/Eut9/RdrVNFgsFpjNZsn8r+7s2rsTo8aPdGy/MflNPPJQH7fHnDl3BmGh4ZJg1Kp5q0LLNm7QxPF7rZoFCxanp193BIuSOHriCPYc2I1Fixc6XrPZbTCZTMjLz4OPt0+Jz/UvrUaLls2cda4fXR+BAYE4deYUWjVvjeFPjcDrM1/Dj7+uxr1334secQ8hqm69gvocP4pjJ4/i53U/OY4XBAF2ux0XLqWgfkyDUtenqmJQI6IS8/HxQb3IaMlrqVekK9Z//r/P8OWKpXgtcQoaN2gMHx8fvDVvFiwWS4mu0aJpS6xe5vzjHBJc87brLfZvkx3gXHXZLtiLKl6o3LxcvDB6PLp37e6yT+9VNk9FL4wej4fjH8HWvzbjj7+3Yv6nH+C9WR/gwa7dkZuXi4H9nsSQJ552OS48rHqtfsKgRkQetWf/bnTr3A19evYFANjtdpw9f6bETwve3t4ugbM4MfVikJp2GdeuX0PNkIIgKO5EUVI6nQ42e/EBrlnj5jhz7nSp6+mO1WbFoSMHHU20p8+ehjHLiPoxzifImHoxiKkXg2FPjUDilAn4/ufv8GDX7mjWuDlOnj7p0fpUVQxqRJWIl7dXlb9GvahorN+4Dnv274EhMBCLl32Ba9evlWkTWMcOnRBZNwqvTn8ZL497FTm5OXh/0c2OIqWYZrNOeF3k5uYgacffaNyoCXy8fQptSkx4ZiyemzgaEWERiH+gB9RqNY6eOIrjp45j4pjEW3oPOq0Ob74zA6+/NBUajRZvvjMdd7a8E62at0Z+fj7mzp+D+G49UDciEqlXUnHw8AF0f6AHAGDU0NF4YvhjmDF3Ogb0eRw+Pr44eeYE/k7+C9NemX5L9amqGNSIKgGb1Q67zY7wO0LL5Xp2m91lKRBPGTMiASkXUzBy3HD4eHvj8b4DEdflQWRlZ5XJ9QBAo9Hgo3cX4vWZr+GxoY8isk4UXhn3Kp5LHF2q5sC7Wt+Fgf2fwoTXxiMj8wbGjnoBL4we71Luvtj7sei9T/HRfz/EZ0s/hVarxR3R9TGgz4Bbfg/e3t4Y9fRovPj6RKRdTUO7O9tj1tTZAAC1Ro2MzAy8+p+XcS39GmoEBaN71+4Yd7NuTRo2wf8+WY73F/4fnhr9JCAIiKwbhZ4P9rzl+lRVXE+NqBzVqheM0QsHIbRmGDTQSPZx7kfP2r1/N5565glsWLXR0aGisiqvmUzKmw02pF1LxadjluHquXTJPq6nRqRwVrO1VNNWkdSGzb/B19cX9SKjcT7lHGb935u4q3XbSh/QyLMY1IhIEXJyc/Duh3NxKfUSagTVwL3tO+LVCZMrulpUzhjUiEgR+vZ6FH17PVrR1bgl/Xr3R7/e/Su6GorAabKIiEgxGNSIiEgxGNSIiEgxGNSIiEgxGNSIiEgx2PuRqJLQpl6CJiO9+IIeYAsKhrWaTXRL1QODGlEloE29hOgB8VDn55XL9ezePji7cj0DWyEat2+Aj95ZiLguD+LCpQvo1qcLVn/1E5o2blai48tzdpDk3dvx9HODsXPTHgQGBLpce8GnH+D3Lb/jx+VFrzyuNAxqRJWAJiMd6vw8XJ7xfzDHlHxdr1vhdeYUwqe9CE1GeqmD2uXUS5j/6Xz8mfQHMjJuoFbNWujW+UEkPDMWNYJqlFGNpYY8+xSaNGqGKS++XubXCg8Nx7a1SR5/b5OmvwJjthEfv7vIo+eVGzH4GQx+3HU5GiVjUCOqRMwx9WFq0qKiq1GolAvn8cTIAYiOisG8me+hbp1InDh1Au/Mn4M/k7bimy++Q5AhqMyub7aY4aXz3AoDJTmfRqO55dW4KwM/Xz/4+fpVdDXKFTuKEFGJvDF3OnRaHb5YsAR3t+2AiLAIdO7YGYs/+hJpV9Lw3sJ5jrKN2zfA71s2SI5v17UNfvj5e8f2OwvmIr5/HFp3aoFufbri/YXvwWJ1LiS64NMP0Oep3li5+hs80KcLWnVsjknTX8GOPTvw5YolaNy+ARq3b4ALly4AAI6fPI5nxo1Am/tb4d74Dnh52otIF+Uohzz7FGbMnY5Z/zcTHeLaY+QLw4t9zxcuXUDj9g1w5Nhhx2sbt/6O7v26oWXHZhjy3CCsWvMDGrdvAGOWUXLsn0l/4KEB8WhzfyuMfGE4rly74nhfq375ARu3/u54D8m7txd6fbvdjk8WLyx4/52a45GnHsa6jWuLrbf8Hv5r0vRX8PxLz+HDz+bjngfb464urTFt9lSYLc5JR9dtXIveA3uiVafm6BDXDsOefxq5ebmO/StXf4OHBsSjZcdm6PFYdyxb+VWJ61Me+KRGRMXKyMzAtu1/YuKYRHh7e0v21apZC717PIK1G37B9FffcKwmXRw/Xz/MnjYXtWvVxvGTxzB11hT4+flh1NOjHWXOXziH9ZvW48O5H0Gt1qBOeB2cPX8GDes3wrhnJwAAgmsEw5hlxNDnB2NAn8cxOXEKTKZ8vLtgLiZMHocvFzr/6K76ZRWe7P8Uvv7vN7d0H1IupmD8pBcwZOBQDOjzOI4cP4y3P5jtUi4/Px9ffPU55r7xLtRqFV6e9iLefn8O/m/mPIwY/AxOnTmF7JxszJ72NgDAYDAUer1PlizCT2t/xBuTZiA6Mho79+7Ey9NeRHBQMO5u2+GW3kPSziTovfT436JluHj5IibPeBU1DEGY+PyLuHLtCl6cMhEvj3sFcV26Iyc3B7v27sS/i7n8tPZHfPDJB5j28n/QtHEzHDl2GFPfmgJfH188+nC/W6qPpzGoEVGxzqWchSAIRS70WT+mPjKNmUi/kY6Q4JASnfP5kQmO3+tG1MWZc2fwy4Y1kqBmsVgw9413EFzDeU6dTgdvbx9Js+BX3/4PzRo3Q2LCS47X3po6B50fvg9nzp1BTL0YAEB0ZD28Mu7Vkr3pQnzzw9eIqReDV8dPAgDcEX0Hjp86jkVffCwpZ7Fa8MbkGY4VAgYNGIKP//shgIJg7q33htlidtu0aTab8MnihVj80VK0aXUXACCybhR279+Fb1atuOWg5qXT4a1pc+Dj7eP4cjB3/tsY/9xEXL12BVabFQ92jUed8DoAgMYNGjuOXfDpB5g0YTK6PxBfUJ86kTh55iS++eFrBjUiqnqKW35Rpyv5+oa//vYLvvxmKVIunEduXi6sNiv8/fwlZSLCIyQBrShHTxxB8q5ktLm/lcu+8xfOO4Jac1G+ctHij/HJYmdHjV++XYeIYjrOnDl/Bi2aSa/RqpnrNX28fSRL3tSuWRvXb1wv9n2InUs5h7z8PIwYO0zyusViKXFPzMI0bthEspp3m5ZtkJubg8tpl9GkYVPEtr8XvZ/siU733IdOHTohvttDMAQakJuXi/MXzmPKm5MxddYUx/FWmxUB/gG3XB9PY1AjomJF1a0HlUqFU2dP4kF0d9l/6swpBNcIRmBAwWKPKpXKJQBarc5FSfce2IOXpiXihdHj0eme+xDgH4BffluDxcs+lxzj4+1bovrl5uai630P4KUXXnbZV6tmbef5fJznG9jvKTwU51wZurao3O3SaqV/Wgu7H8X5N4/1yXufIbS2dEV0T3aYEdNoNFj80VLsObAHf23/E//79n94b+E8fLv4e0cgfHPKLLRu0VpynFqtKex0FYJBjYiKVSOoBjp26Ijl3y3DsCdHSPJqV69dxc/rfsJTAwY7XguuEezoGAEAZ8+fRZ5oDN7eA3sQERaBMSOed7x2KfViieqi0+lgt9skrzVv0hzrN61HnfC6LgGlKEGGoFL31oyJisHWv7dKXjt4+GCpzgHcfA82m9sy9WMawMvLC5fSLt1yU2Nhjp04ivz8fMe/4b5D++Dr64fw0HAABQG4beu2aNu6LRKeeQFdH7kfv2/5DcMHjUTtWqFIuZiCRx7q47H6eBqDGlEl4nXmVKW9xtSXp2PgyAEYOW44Jjw3EXUj6uLE6RN4Z/7biI6KRsIzYx1l72kXi2Ur/4c2rdrAZrPj3QVzodM6mybrRUbjcupl/PLbGrRs1hJbtm1x6S1ZlDrhdbH/0H5cuHQBvr6+CAoMwlMDBuPb1d8g8fUJeGbIaAQZDDiXcg6//vYLZr7+FjQazzxJPNHvSSxZvhjvLJiLxx4ZgCPHD2PVmoIenSXtIAMAdSLqYNv2P3H67GkEBQUhwD9Acn8AwN/PHyMGP4PZ896CYBfQ9s62yMrOwp79e+Dv53/LOSyzxYIpMydjzIjncfHyRSz45AMMHjAYarUa+w/tQ9LOv9Gxw30ICQ7B/kP7kH4jHXdEF+RSx40eh5nvvokA/wDcF3s/zBYzDh0+CGNWJoYPGnlL9fE0BjWiSsAWFAy7tw/Cp71YLteze/vAFhRcqmOio6Lx3dJV+PDT+ZgweRyu37gOQRDQvWs85s54V5KneXXCZLw241UMGvUkateqjdcSp+Kfo4cc+7t1jsPQp4Zjxtw3YLaY0aVjF4wZkYAPP5tfbD1GDH4Gk954Gb0e74F8Uz42/rgFdSPq4uv/fot3F8zFyBeGwWw2IyK8Du6LvQ9qtedGLkXWicQHcxbg7Q9m48sVS3BnyzZ4bsTzmD5nWqmaBB/v+wR27E5G/6GPIjc3B18u+god2t7jUm7CcxMRHBSMT5YswoWLKQgICECzxs3x3PAxt/weYtvHol5kPQwa/RTMFjMe7v4wXhg9DkBBIN25ZyeWfr0E2TnZiAirg0kTJqNzx84AgAF9n4C3tw8+/99nmDt/Dnx9fNGofiMMfbL44RHlRSWUtqG3EjAajTAYDOiCPtCqSp6YJqpoteoFY/TCQQitGQYNpE8PVXHux/mfvI/Fyxdj8YcFf+Cro4VffIwV3y/H1l+2VXRVilVeM5n8ywYb0q6l4tMxy3D1nPT/tlWwYAt+RGZmJgIDAz12TT6pEVUS1rCIKjcX47hnJ6BOeF3sO7gPrZq39uhTUWW1bOVXaNmsFWoYgrD7wG58/r/PMOjxIRVdLbqJQY2Ibkv/Rx6r6CqUq3MpZ7Hwi4+RacxARFgEhg8aiWeHPVfR1aKbGNSIiErhtcTX8Vpi2U+mXBbmTJ9b0VUoc8pvKyAiomqDQY2oHAmCUOpBuERVWXn/n2dQIypHWddzYDXbYIf7gbdESmCHDVazDcZrOeV2TebUiMqRKceMnT/vx30D9QgOCoYalWd6ISJPssOG9Ix07Px5P8y55uIP8BAGNaJytmXx3wCA9r1bQ+ulKdVMFERVgSAIsJpt2Pnzfsf/9/LCoEZUzgQB2PzF3/hrxS4E1vRjUCPFEQQBxms55fqE9i8GNaIKYs4149r58v/QEykZO4oQEZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFiMKgREZFieDyo2Ww2TJ06FTExMfDx8UH9+vXx5ptvQhAERxlBEDBt2jSEh4fDx8cHcXFxOHHihKerQkRE1YzHg9rbb7+NhQsX4sMPP8SRI0fw9ttvY+7cuViwYIGjzNy5czF//nwsWrQIycnJ8PPzQ3x8PPLz8z1dHSIqayqV9IeoAnl85eu///4bffr0Qa9evQAA0dHR+Prrr7Fjxw4ABU9p77//Pl5//XX06dMHAPDll18iNDQUq1evxsCBAz1dJSIiqiY8/qR27733YuPGjTh+/DgAYP/+/di2bRseeughAMCZM2eQmpqKuLg4xzEGgwEdOnRAUlJSoec0mUwwGo2SHyIiIjmPP6lNmjQJRqMRTZo0gUajgc1mw6xZszBo0CAAQGpqKgAgNDRUclxoaKhjn9zs2bPxxhtveLqqRESkMB5/Uvv222+xbNkyLF++HHv27MHSpUvx7rvvYunSpbd8zsmTJyMzM9Pxk5KS4sEaE1VT8lzYrf546rxEHuDxJ7WXX34ZkyZNcuTGWrZsiXPnzmH27NkYOnQowsLCAABpaWkIDw93HJeWloY777yz0HPq9Xro9XpPV5WIiBTG409qubm5UKulp9VoNLDb7QCAmJgYhIWFYePGjY79RqMRycnJiI2N9XR1iIioGvH4k1rv3r0xa9YsREVFoXnz5ti7dy/mzZuHESNGAABUKhUmTJiAmTNnomHDhoiJicHUqVMRERGBvn37ero6RERUjXg8qC1YsABTp07F888/jytXriAiIgLPPvsspk2b5ijzyiuvICcnB6NHj0ZGRgY6deqEdevWwdvb29PVIaq+PJmnUpWiUUew3+I1ZPUVTdhAVFIqQah6/3OMRiMMBgO6oA+0Kl1FV4eocqpqQc3lPFXuTxOVglWwYAt+RGZmJgIDAz12Xs79SEREiuHx5kciKke38zRWmqevsjqvu6c6NkfSLeCTGhERKQaDGhERKQaDGhERKQZzakRVSXE5tFLks1Tq8p+aSrDL8mLi+hbXa1L83plfoyLwSY2IiBSDQY2IiBSDQY2IiBSDOTWiyqY0Y89kOTS3eTJ3+TZP5tfkeTNJFYrOmwl2tfyFoq/BMWxUBD6pERGRYjCoERGRYrD5kaii3UY3fZfmRnHZYpoUVeLrqm/j+61d1kyocZ7XZb50N02M8vdSquZIopv4pEZERIrBoEZERIrBoEZERIrBnBpRRXCXR3PXTV+eX5PlodzmyeTXvJ08mrvzivJoKlm+TZBf0l6aabLclGUXf7qJT2pERKQYDGpERKQYDGpERKQYzKkRlYdbzaHJ9qs08pxa0XkzlUYjK+tmTJucvKybqa/k+S3BJtqWvW+VzSYtK66CXX4fZOd1Vweim/ikRkREisGgRkREisHmR6KyonLTFV9crJhmQkk3fVlZl2PF2/KmymKGA5SYvBnQLm1SVImvI29ulJ1K3Bzp0t3fBqJS45MaEREpBoMaEREpBoMaEREpBnNqRJ4i77Yv7opfiiVi3OXJVFrZR1aeNxOXlXf3l5+3NCtsi6edkufJ7LLz2Nwkw+TTV7k5r+tQB2cXfy5LQ0XhkxoRESkGgxoRESkGgxoRESkGc2pEt8pNDq1gs+hxapLproobe6YTfUxlOTWXHJtWdKwspya4jFtzk1OT5b5U4rFpNtl7sRW9DIwAq3SXbDCaIFqaRiWrj+Ayqo2oeHxSIyIixWBQIyIixWDzI1FpuJn6qlSz64u73utkH0Odl/Q04iZFL51knyA/Vtx0qZU2YwqaUnThl/eQF3e3t8ret8WKIsm78MubKsVNpO66+wOcNotKhE9qRESkGAxqRESkGAxqRESkGMypEbnjptu726mvIMuj6aS5MElXfFmeTCUrK94vz6EJetm2zplHc+nCL6uvIF4lW96F3yrLfVnFU37JuuVLS0J8FZeVrjXy4QDi7v8eJP53k+fmSNH4pEZERIrBoEZERIrBoEZERIrBnBqRWHHLsYjzZu5yaIAkj+YynZUoT6byko5LE7xl23ov0e/SsWd2L1lOzctZB7vLtFgokmQaLAAqq3RbbXHmxuTTWallOStBtK2yyqb8ssoGm0mmEivFveegNSoCn9SIiEgxGNSIiEgx2PxI5K7Zy6WJUTS9lZuprwBpk6NKL21ShKRJUdbc6CPdtvs4mypt3tKPrE0vrYNd59wWZJ9uQfY+VZJmQmlZtUXapV9jch4ra1CUzLRfcC5R06C8C79sW7AW/b36tmbtZzf+aotPakREpBgMakREpBgMakREpBjMqRGJFbecjHhbnkOTdc2X5NG89ZJ9go9z2y7rwm/zk06TZfN1fkytPtL62byk9bPrnNv2YpaaEXfj11ik+zRmN8fZZO9bNqWWyiLabylmtW3x/ZRNvwUbu+1T6fFJjYiIFINBjYiIFINBjYiIFIM5Nap+5HkdlXhZFdk+d2PPdEVPfQVAOhbNR5pTs/s5t+U5NIuf9LxWX2f9LPKcmvS0sItOJcjfi2zollo0Ns0uy6EJ8vyW6Fi1yxRasrFnomPlY82Iyhqf1IiISDEY1IiISDHY/EjKV0wTmKTJ0c20WC7bOlkXfnm3fV9vx+/i5kYAsPo7j7UESD+GZn9pHSx+zvpZfaTvxeYt2YTdzSdaJW9+NIu6/+cXfRwAqK2isjpZc6PLdGGi1axlTaDyVbNLRbAXX4aqPT6pERGRYjCoERGRYjCoERGRYjCnRtWPSj51k6gLunypFHcrVsuXjJFt232d2+IcGgCYDc7zmgOk1zQHSPNQFn9xTk1aHbtemiizi1J+8kyifHkZTb44lyjdp5bNUGUTLT2jla094zJ0QFWK1ayJPIxPakREpBhlEtQuXryIwYMHIyQkBD4+PmjZsiV27drl2C8IAqZNm4bw8HD4+PggLi4OJ06cKIuqEFU7/l5eCA/wh4ZPSVQNebz58caNG+jYsSO6du2KtWvXolatWjhx4gRq1KjhKDN37lzMnz8fS5cuRUxMDKZOnYr4+HgcPnwY3t7ebs5ORGJhAf7oGF0PDWqFoFGtmmhYMwThgQEAAKvdjivZ2UjJyMTm02ew9tAxXDJmVXCNicqWShA8u+75pEmT8Ndff+HPP/8sdL8gCIiIiMCLL76Il156CQCQmZmJ0NBQLFmyBAMHDiz2GkajEQaDAV3QB1qVrtjyVA2p3Iw9czMVljyHJh97phJ96RKPQwMAu7804WULdObRTEHS/6cmg0b0uyyHFgDZtvMjavMp+L19RB0Ma90G3es3hFatxrnMDBy/cQ3HblzDiRvXYDSZEOYXgAj/ADSqUROdI6PhpdHip2NHMPuvP3AlJweaPOd1dVnSOuhksU+fKYh+lybc9BnSdWs0RuecW+rsPMk+Va50QJyQny/63STdZ5UlAUVL0Qh22Z8t+Rg2z/5ZozJgFSzYgh+RmZmJwMBAj53X409qP/30E+Lj4zFgwABs3boVderUwfPPP49Ro0YBAM6cOYPU1FTExcU5jjEYDOjQoQOSkpIKDWomkwkmk/M/vNFo9HS1iSo9vVaLXs0aY1jrNmhWqzZOpadjxt+b8MOJf5BlvhlIxLFJ9HfdT6dDv5gWmNAhFhtihmHMLz8j+URKudafqDx4PKd2+vRpLFy4EA0bNsT69esxZswYjBs3DkuXLgUApKamAgBCQ0Mlx4WGhjr2yc2ePRsGg8HxExkZ6elqE1VqDzdtjM3PjsCcbt1xOTsLT6/+Dg9+tRhL/9nrDGhu5FgsWHZwP7r9bzH2pl7G4j790LNJo3KoOVH58viTmt1uR7t27fDWW28BANq0aYNDhw5h0aJFGDp06C2dc/LkyUhMTHRsG41GBjaSKk2nCDdTYcln3lfp3M28737FanOAc9sUKO0HbwpSiX6XXsJicDalRRkMmNG9G+6vG4NfTh/FOwc341z2jYKdIYBOK2t2E82FZbNKr2nXapGBXIzY+D3+7/6emNsrHsezruPY9WuSKbMAQCttNYQgvmXyyf9l914lbvqTNwPaZfWVNyOWFKfMoiJ4/EktPDwczZo1k7zWtGlTnD9/HgAQFhYGAEhLS5OUSUtLc+yT0+v1CAwMlPwQKd2TLVvh1yFDEWMIxvB13+H5jT85A9ptsAp2vLJtHc5mZODDng/DRz4Wj6gK83hQ69ixI44dOyZ57fjx46hXrx4AICYmBmFhYdi4caNjv9FoRHJyMmJjYz1dHaIqJ8THB5/16YtZcQ9i9ZHD6P7dF9iUctqj1zDZrEj4dQ0iAgIxqdP9Hj03UUXyeFCbOHEitm/fjrfeegsnT57E8uXL8emnnyIhIQFAwaKBEyZMwMyZM/HTTz/h4MGDePrppxEREYG+fft6ujpEVcodNYOx6slBaB0ahmd+XIXXN/6OXKul+ANvwekb6fh09070a9ocPvIFT4mqKI//T27fvj1WrVqFyZMnY8aMGYiJicH777+PQYMGOcq88soryMnJwejRo5GRkYFOnTph3bp1HKNGniOZ+kq+fIybqbBky8lAPhWWt3PbJltOxhIg77bvvI44hwYAJuewTZhrFOSH2oZH4L+9+yItLxsD1n6Fy7lZQAjg4+/s+evnLe0U4qWVdnu32Z3XzDVL65Orkn6+bFYVvjt9CBPuuRcPNm+A1f8cceyza6T1FefU5Dk0d1QuXe/d5NBc8m/SbZdu/ESF8Pg4tfLAcWrkQv6HtjRBzUsUuPSycWmysWiCn3MsmtUgHZdmDpIGwPxgjeh36TVNwaLjatgRX78BPojvhb2plzH6z+9gtDgDmbenglqOLKgZC/Z/03sgbGY7hq743rFPny69n97pzj8T3jeknTT0N6RPkroM59gztTFXsg+50h4ogsks+l02Ts0ifW+CaJxasR1Fqt6ftWqnrMapce5Hogo2tFUbfNzzEfx2+iSG/vi9JKCVh++PH0JsvSiEB/iX63WJygKDGlEFeimuE6Z3eQCf792N8et+gdlmK/4gD/vl9DFYbDY82KhhuV+byNOYHaaqqZi8jmQqLDfTYgEARJ0kVF7SJju3y8n4ST8+ZvlYNFFOzRwkvaQp2I7EDvdi1N3tMWPXBiw+thMIKdgXZJA22dXyy3H8btBLm++0KmkzXK7VWb8MrWzaLpv0O2xuXkF9c2DCuYwM1AsOcubO3HzdVcma9lQ2N2PRZGVdsh2eGm/G5ka6iUGNqAIMatEK4+6OxVvbtmLxuZ2lP4EgIORKDiJTbiAnwAupkQbk+XkVf1wRLmVmISIwoPiCRJUcgxpROXug8R2Y0bkbFu/bg0/37oIquPhjAECfb0HcusPovO4Eok6lwy9H2mnkRogP1j/cHD8+fidyAvRFnKVwF41G3BkRXqpjiCojBjVSBrerWct6P8rHZIm79MubH91MhWUJlDc/Sps5zQbn76ab3fbbhIdjXv+eWJ9yDDMPr4cqWECtEOmU+HUDMiTbzTIvocuK44j98TR8ssw43ikUf45sgCv1A5EWHQifLDNCzucg6mA6+q3Yi/i1h/HFmx2xuWljyXmM+dJAl6txNtldyjKiV2DjopsdRa17Klnaz6XbvrhFsZhu+pJt+RRanAqLbgGDGlE5uaNGDXz+yKM4kJqGiX//BHsxeSC11Y6HVx7Ek5/vgE2nxt9962PvwCjcqOvnKGMXCgLphRbB2N8zEmsGtcLwaX9hzItbcPTTUFyuG1Siul00ZiHIxxt+Oh1yLGUz2JuoPLD3I1E5qOnriyV9++Nqbg5G/bwaZrv7Xo7h5zMwe8wqDFm4HX/3rY9pa/pg9YQ2koBWmBthfvh4XhdkBXtjypRfoc8rWYC6llPQGaWmn/vzE1V2DGpEZUylAt7r0RN6rQbDVv0Ao8n9OLTYzafwfyO+g1+WGa8t7IvvX26LfP+STzKQ7++FRe90RvjFTPT9dl+JjtHeHJButlmLKUlUubH5kaqO0qxmLd6Wd+GXr24tmlFEPA0WANh8pdsWUTd+s79slhDZCtbmoILmxdF3t8O9kVEYvGEFLnmlA8FAzeBsR7k7DNcdvz/0zUEMmb8dh+LrYPX0O+Hla0c93SXHfl+1NCDm2qV5Mp264AnQ2lSD5AdiEPfbEfwyoiWgUkGrkeaoxLdTryl4XyaLDRBc82Zq0bZrF355F39RYZv7PJm4i38VnNyIKiE+qRGVoZahoUjs1BGf/LMdf6eeK7Kcyi5g8PztGPrBdvw1rAG+m9MWFt/b+875R89GCE8xotGhK8WW9b4Z6E1WPqlR1cYnNaIy4qvT4f3ePXHkylXM2/dnkeVUNjsS3tyCezecwuLEe3F6WC2PXP+fuyKQXtMX7f44i+MtQ92WDfLxhtlmYycRqvL4pEZURqZ164pafn6YuOZXWOTd1f8lCBj6QRLu/f00PnizG9Y/1txj1xfUKlyIqYHQi1nFlq3p64frubnFliOq7PikRsogz7GJ8mhul5oBIOidnTDs8pyaj7Ssxd95XnOANIdmCXTmhHo2a4QBLVvgxS2/4iSuITBIGjCiDekAgD7L9qHHd4exZmpLZPXxQTNcRiPvy5KytTXOoKSWTYuVbpNOQizOsek1VtwI90XkiRvQa1ybFQW7s/4hPj5Iz8nFv6d3yalZne9NnlNTW2SFxXk0eTCX59jE+bfilpoRip5+i+hffFIj8rCwQH/M6BmHn04dwXcnDhVZ7r71xzHk42T88WxD7Hm8XpnU5VqYP0JSs4stV69GEC4ai3+iI6rsGNSIPOy1B7sgz2LFlG2/FVmmyf7LSJi1FZt6NcaWhEZlVpfsIG/4G90PIVABaBkWioOpqWVWD6LywuZHqrzczMTvtgs/IO3GL2tuhE425ks0NZZd1txo9ZMOBzD7O69jlq1raAkUcH9UPfRo2gjj1v2CfF8jvHwL9tUxZDrK+Wfl46UZG3ChVRC2zmiEFj4XJeepr7sq2TaonZ038gXp91CLIK2fTtRuaBNUUJvtsHhpYBNUsMpm6Ye14L3UC6yBQG9v/HMhDf9eSi3rL6K2OJv7NGbZ0ACrfNtZB0Hem1K+tE5RuUaiW8QnNSIP8dJoML3zA9h+IQU/Hz9aeCFBwAvzNkOfY8F3c9rCrivbj6A+zwKTt/vvrq1qhgEADl1OK9O6EJUHBjUiDxl+512INARh2paNRZZ5cN0RdN58Aj9NuxOZ4b5lXifDjXxkBXm7LdO6VhhSbmQgIy+/zOtDVNYY1Ig8oIavD55v1wHLDu7HifTrhZaJuJCB5z/YivUPNcOh+DrlUq+ga7nIDPZxW6ZlzTAcvMSnNFIG5tSo6hB323fThb9gW7T0jJsu/ABgF227dOH3leWw/Jw5Nau/M8/03IMdABXw3uE/YQ8syCMFBzpXrL7D7xpemL8JOSF6rJvcAq315x375Dm0ulppHkqncr63TNlEyBpZ93pxji3fpkPUsXTsvjcK+TYdTBbZ9GBmNdQqFVrUDMVH+7dDbXK+N41Jel6NKKemsshyaC5d+m2F/w7XqbAk21xqhjyAT2pEtykmuAaebNMKHx74GzdMeYWWabPxPJruSMW3L7eDybfkkxPfDj+jCeEXjTjVtOgZSlrUCoWfzgsHLrDnIykDgxrRbRrb6R5cyc7BkiO7C92vN1nQ/709OHBfHRzqVD7NjgBQ/2jBE+CpJrWLLNPzjka4npeLvSmXiixDVJUwqBHdhvDAAPRs2hifJ++CSd5d/aaBq3bBcC0P3ye2Lde61T96FTl+XrgcaSiyTM/6jbD+9AnY5LN3EFVRzKlR5eFmXBogG5smH5cmP9btODXptl3U5d3qI83NWXxlU2H5i38XMOTeNsgxm7Hi9CF415IOcg73MyIgIw/Dvk3C3icj4dvIhHooKBPtdc1RLkKWQzOopb0V7XDmmrIgDZz5grQpM8vmPPaBXZdxtGkYMq0FHUVMJmnZloGhiDIEYd2G36GRjc922c4XjVMzyYK3RTYWzepm6RmXcWqipWfcTYsFcGosKhE+qRHdogAvLwxs3hJfHdqP3CJmt39o1T/Q2AT8/Wz9cq2b7w0TWu+9gG33FX3dno0aIT0vD9svpJRjzYjKFoMa0S0a2LwV9Botlh7YW+h+ncmKnj8cxMaejZEX7FVombLSZFNBx4+/7mtQZJmHGjXChpMnYeWsHqQgbH6kyktV9Hculay5Ud6lX9LkKGtuFLzkzY/OY60+0vNafeXbBU1gWrUaw+9sg9WnDiMNmYAfYPCT9nzsu/kAAjPzkfx0fTTWSWfer6V2ztrvr5I2C4q78ANArqgbf5Zdek8ybNIB3FfNAQCAhuuvYG/LSJz1CQFujqm25Trfd7Pg2ogOqoE31m6CJk8FrazTpjZfGujETY4qs7S5USVrfpRMjVXctFjsxk8exic1olvQs2kjRPgH4rNDOwrdr7IL6LbsCA50icTVyIByrZtfhgmNd6ZiY6cmRZbpFdMEGXn52H6WTY+kLAxqRLdgxD1tseXCaRzPuFbo/nt2n0bYWSN+H9y0nGsGxP50Cna1Cps7Ni50v16jxZONWuPnf46y6ZEUh0GNqJQa1gpB87BQfHlkT5Fl4rccxqU7DDjduuiBz2VBY7Gh64qj2NkjGuk1/Aot069Bc9Tw9sGS5MLH1RFVZcypUdUhzrGpZd/H5NNkiXJqQjE5NZtelFPzlp7XKps20eYjoEfLRjDm5+PvGyeh8XXmjIK8CxJTOrMVnZNPYMfTMajlVbDwZpBGuvJ1gNr5hKSR5QctgjQPlS04e1ZetUsDVYolWLJ9xy9XUeNKHr4e0B6ZOdLKq3K0UAEY1exurD91ApcuGVHwCqDNlXaX1+bJVrfOLzqn5tql37ktFDNNlnQnnxrp9vFJjaiUejRsiN9Pn4LZXvhg67a7zsMvx4x/ukeUb8UEAX1X7MPuDpFIiQkutMgD0Xegfo1gfLZ3V/nWjaicMKgRlUKD4GA0qlkTvx4/XmSZ+7aexLl6wbhav3w7iNTbno6Yk9fx48A7iywz+q522H35EvakXi6yDFFVxqBGVAo9GzVGlsmEbefOFbpfbbMj9q9T+PP+oseHlQm7gM7vHcfxprVxoG3h80u2qh2KDnUi8eneneVbN6JyxJwaVSxxPkm+nIx8Kizxtnxcmkb2/Uwr2q8rOocGADZvlWgfZGWl2z0aN8Tv504hX2eFj146i0igVz7qnk6Hf44ZJ9rXRoDGueimt0paVpw9yhekOSmLLLd0yeqs/1lzTcm+c3kF2+3WnUH4P0ZMf+dhXM4tmOsxP0ta+dF3t8O5GxnYcuA0dIIKumznPp1LTk3atKrOF9XRTQ4NgHRsmjyHJpsKy2VqLMlOTotFpccnNaISiqlZA01CauHXk8eKLPPvzPinm5Rfr0evPCseXbAP+7rUxb6WkYXXyxCMHk0aYcnOPbAzWJCCMagRldADTeojx2LGHymFNz0CBUHtYlQQ8vzKb1qshz4/CL8ME76fWPQqAJPadsFlYxa+3Xew3OpFVBHY/EiVl7w5UtRUKZ8my7U50rkt6KT7BJ30vDaduPlRel673vlU0yoqDPuvXka+xgxoAL1O2uzmrbGg4bErONMsBN4aC3QqZzOcDdLz5oimu7LJuvBn2qXTZp2yOJ/6TplCpXU/ADzw1VF8PaQ9dhrqIf2Ks8u/2ljw8e5Qpy66RzXEiyt+gT3T7ujGr8txvjddrmxarDzpe1ObnM2nKpO0KVWQNT8Kopn55V36XWfeZzd+8iw+qRGVUIuIUBy45n6F6MhTN3C+YeHd6T1Nm2/DpDfX4UJkDawcWPhTmgrAlPu6YG/qJfx6sOgem0RKwSc1ohII9vVBHUMgDlwtOqjp8yzwzrPiRk3fIst4Uvx7/yDiYgbGLRoIi77wj/KAZi3QsnYoHlv5dbnUiaii8UmNqARahocBgNug5p9RsLJmVpB3kWU8pdHWVNzz9Rn8d0wnnLsjpNAyBr03Xu14P1YdPYzdly+VeZ2IKgM+qVHVIZ4ay820WAAgiLr0C1pZDs1Lum0X5dRk6SzHdvO6obiRl4eU/AzHp0arkeaDfLMLck0WXw00KgEWwVmHLLt0yqqronybXZDW54pNOmj7eH644/cjWWEIvpqDl6aux66OUVjU9T4g01l/c4azG//0th2hU6nxzto/oMtRwcsofW9e2c6cmjZH1oVfllODWZRHky+I6q5LP7vwUznjkxpRCbQMC8Wh1DS3ZbSWgj/mVlnHFE/yzTbhtZd+hcVLg48nd5GO8xNpVTMMT7Ztjfl//I1rObmFliFSIgY1ohJoERaKg5fdB7VsQ8FTUuCNfLflbpXaYsfLr/2GmmnZmPl/PWGs4VNoOV+tDvMf6I1/Lqfhq537yqQuRJUVmx+pfMmfLETd9t3OIAJImxzlM4i4zNrv3LbLu/RrZN32tSrRPulpBDWgU6sRGuCPs5kZUKmLbiJLqxEIm0YFw6VcmGxaZNqcQSfVYpCUNdqcebdcu3Tmj8uWIMn2iezaUNkFjH1jMxrvT0PCrIHYGxIFZAJZGdJOKdoMLWb0eBC1ffzw/Ceroc4E1De78HtlSeuuy3Y2E7p04c83S7bF3fjddeEv2BbPKMIu+1S++KRGVIwaPgXB6Xqu+2Y8u1aN9Fp+qJma7bZcaansAp5+Pwmxv5/C1Jd7Y2/LqCLL9mrcCANatsAbv2/C+WsZHq0HUVXAJzWiYtT0K3gaul6C3FRanUDUO5nusWtrzDaMnb4ZsRtPYfGLHbHxvqJX0q7jF4iZ3eOw5shRfP/PYXij8HwbkZLxSY2oGEHeBU2FGfnF58qSu8Sg1Y4LCEzPu+3rehsteHrMdrT/4yw+eLMbNvRrVmRZjUqF9+/rjSyTGa9v2Hjb1yaqqvikRpWWfCoslZsZ/eU5NkGcf5Pl5lxyahrxPnklAK+bs/zn26ySXudWm/Sa2VY9NnRtiqfnJ6HdurPYMyzasc9il55YrXKeyGiVjmu7lGdA6EUjJk1eh+CrORg7ayD2tYgEMoHr6f7S6qUXjDlI6HAP7qpZB0//dyVMaWZ4QQV9pjSH5pUtzW/pcpy5MU2uLIcmy6kJZtG2uy78gKQbv0sXfpdpstiNnzyLT2pExdBrbgY1+R/zQmQHeiO5Ywz6rtwPv5uDsUtFENBtzRHMG7YS+nwrXlvYtyCguXFfVD2M7xCLBTu2Y885DrKm6o1BjagYXjefAu32kvXk+zyhE/T5Fox4fRtUtpL3/qt91ogXxmxGwpyt+OuB+khc8hguRtdwe0yL2rWx8OFHsPXcWXyYvL3E1yJSKgY1omJcvdlBpJafXzElb5YPDcDc6fFouiMVI6b+BV+jmyc2QUDDnWkY9eKfeH3Arwi+lIMZ/9cLH0/ugnxf98vXRPoHYXGffjiZfh1jf/kZNjblETGnRpWIPE8mH3smzo3Jx7DJx7+JDhVccmry67qpkx24lFkwv1REQCBO5lx17DJZpB8fo8mZG/uzRUPMez0Oz83bihmP/IQtDzfC+QbBuBxpgKBWIeJ0BiJPp6PVzouodzodZ+sFY/74rtjQvRkuWg3ADed5MzOdY9GE9IIxbTW8fbCk50Bk55vx3NIfYc+1Qw81vDKdx3kZZePSjNLmU02OaDmZPGkOTTItFiDJowkW+Tg1N8vLcJwalTMGNaJiXDZmwS4IqBMYCOSU/Lg/4hrh4J110H/5HnTceAoPr3Au0GlXAWkRgTjdpBY+fL4r9rWp6wzMxaTuvLVafNG7LwL0ejz16Te4kXv7PS2JlIJBjagYFrsdV7JzUMcQCFwu3bE3avrhv+Puw/LEDtDnWRB60Qi1TcD5qBowexf0XMw0FT7dVWE0KhUWxD+MRsE1MXDVt0i5kVn8QUTVCIMala0iJty9tXOJptRymW5L3vzoLCvvwi+4q5OstUx1c/tSphERAQGAxXleU74055Uum7U/3+r8eGnUBfuOh9YGAJjztcDNYW85svPkZ0unzVJlFAQ/L40G73V5CF3qRWPM1z/i+Kmr0GdI66vPdNbByyhtFtTmyKbCEnXjV5nkXfhlzY/iJkc3XfgBzsRPFYtBjagEzt64gRahoRV2/UC9Hp8+3Aeta4dj/Mpf8OepcxVWF6LKjL0fiUpg7dHjaFy7JpqF1C73a4f7B2DlYwPRKLgmhn35HTYeO1XudSCqKhjUiErgzzPncC0nB/0bNS/X6zYJqoUfHn8SPlodHlv5NfZeKGVSj6iaYfMjVShVabrpu9kn77bv9pqyvI5KlApTy9JFakvBee0Q8PPBo+jbqhne3vonrHY7bLLVtnME6XRX+VrRMtqCtH42q2hpnFzpeTTZBWMOYiMjsSi+D86nZ2D0stW4lp3rmkPLkL4XcU5NlyXNi2lypOPlJFNhyXJqsMhybOIu/e6WmgHYjZ8qFJ/UiEpo1YHDqOnri67RMWV6HY1KhbEd7sGSR/tjz+VLGLJ4Ja5lc/VqopLgkxpRCR1Nu4YdFy9g6v1dsf3CBWTC84GmXmAQ3u/1MFqHheGjHcn4MHk7VObijyOiAmX+pDZnzhyoVCpMmDDB8Vp+fj4SEhIQEhICf39/9O/fH2lpaWVdFaLblrh+LQL1erwd192j59Wq1BjTugPW9x+OYF8fPP7tCryf9DesJZxvkogKlOmT2s6dO/HJJ5+gVatWktcnTpyIX375BStXroTBYMDYsWPRr18//PXXX2VZHaoM5FNhlYYk/yY7jzzHJtqW59Ag2xTn1FSy2Tw0smkbU69mYfLa3/Dxo49gwvVOmL8jybHPbpUtaaMS59Sk51GbnPXvGhqNyQ90Rv2QYCzZtQcL121HrtkC/c3vnLos53Hy5WT0srFoXqI8mngaLKCQqbBEeTTBIi3rdiqs0uTMOC6NylmZBbXs7GwMGjQIn332GWbOnOl4PTMzE59//jmWL1+OBx54AACwePFiNG3aFNu3b8c999zjci6TyQSTyfnXxWg0llW1iYr12/GTmPfHX0i8vyPC/P3xxtbNMNmKX5ZGzFenw6NNmmFo6zZoGByCfZcu49Gly3DkylXozFyxmuhWlVnzY0JCAnr16oW4uDjJ67t374bFYpG83qRJE0RFRSEpKUl+GgDA7NmzYTAYHD+Rke7XlyIqax8nJePlDevQr2kzrB88FH0bN4W3tvjviPWDgjE1tiuSRozGG50fwMn063hq+bd47H9f48iVq8UeT0TulcmT2ooVK7Bnzx7s3LnTZV9qaiq8vLwQFBQkeT00NBSpqamFnm/y5MlITEx0bBuNRgY2BXKZ+krexFiac4mbvWStZWqrtElMbXFua2QtdHZT0UMHVu88jAOnL2PSA53xXnxPzDTHYdvZczifkYnLxixcyclGqL8/GtWqiUYhIWhYMwT+ej3Sc/Pw9Y4D+HrXAVzOzII2D9CLlgrQZctWrM4SRL9L34xXpmzm/WznG1Dn5Ev2qfKkbamCqPUDFq5mTcrg8aCWkpKC8ePHY8OGDfD29i7+gBLQ6/XQ6/XFFyQqZ6ev38DolasRVSMIPZs3wr31otC9YQOEBfhDr9Ui32LFyevXceLqdWw4dhLHr13H9nMpsOTYij85EZWax4Pa7t27ceXKFdx1112O12w2G/744w98+OGHWL9+PcxmMzIyMiRPa2lpaQgLC/N0dYjKxfkbGVi4fQcWbt/heM3g7Q2jyQS7ILgO6na7iBsR3SqPB7Vu3brh4MGDkteGDx+OJk2a4NVXX0VkZCR0Oh02btyI/v37AwCOHTuG8+fPIzY21tPVIaowGfn5xRciIo/yeFALCAhAixYtJK/5+fkhJCTE8frIkSORmJiI4OBgBAYG4oUXXkBsbGyhPR+JSsSl275zWyXLAalkT00acU5NlkOzy1fJFp1KbZUPIyi6emr5UAFRvNPmyVaozpHl1LLFU19JT6TNliYB3S4nI58KS7S8jHgaLMB1KiyuZk1VRYXMKPLee+9BrVajf//+MJlMiI+Px8cff1wRVSEiIgUpl6C2ZcsWyba3tzc++ugjfPTRR+VxeSIiqiY4oTERESkGJzSminWr02bJ5kSU580EmyinJssPqS3ybWcyTJvvfoyVODdm18qXypFtiuqkls5CJZl+S5cnrY82V7qty3FeVD71lVq+nEyuM1kn5Mvm+JIvJyMai1bc8jEuY9MkOzkujSoPPqkREZFiMKgREZFisPmRlMFlNWtn85nKKmtuNEu3NfmidkOXJkTp9z5xM6JQzFdCtU3c/Citn8bs3NbKmh81ubKpr0RNjup8aROifOoryQrW5mJm3hd345c1P7qdCovNjVSJ8UmNiIgUg0GNiIgUg0GNiIgUgzk1qrxk3fbFy5+45HXkZa3OHJHKIs0XqbXS73LyqbHEXFbCFn1iBNlSOfIVtsXHyocRaEzObU2+9CLqvKJXrJZPfQU3U2EJZlnZ21lOhqiK4JMaEREpBoMaEREpBoMaEREpBnNq5HmqW1sAU5DnpFwKOPM88qVR5FNhiXNs8pyaSi3LqUnOI11rRtDKpotSu3lvsrSUWlQnlVk+Vs6ZR1OZZPWT5clU+W7yZPKxaO7GnsnukWRqrOJyaBybRlUEn9SIiEgxGNSIiEgx2PxIFUvS7CVbZtql2764S7+suUzW1KYSdemXd71Xy7bFzZ7yKbUEWfd/8bnkXfjhssK2myZQk6j5UTZ9lev0Vpai98lWrIaorNvVqwvbluxjcyNVTXxSIyIixWBQIyIixWBQIyIixWBOjTxPnI+5xe79hZ+26NWsBY2si794miz5eWTb4mVqYJXm9UqzMLfKJjuzqA4uOT9xHs0lTyabzkq0YrV8+Ri33fa5ejVVQ3xSIyIixWBQIyIixWBQIyIixWBOjSoPec5HJRtHJc7PFTNOzV0uz2V8mTgPpZblrNzlBOXTermZqguyPJlkfJnLWDPptmQ6q2KmvpJMJVbccjLMo5EC8UmNiIgUg0GNiIgUg82PVGm5zNov2naZpV/WVClA1ITn0kwo+y4nGg6gKs0QBJdmTFnToLj5Ud48Kqq/y1RXLk2Mou1imhQFd1OJsbmRqgE+qRERkWIwqBERkWIwqBERkWIwp0blyyUH5PxepVLLl0KRLUUjzlHJc1/ynJXkItJckqB2c97SKK7LvK3onJrbbvpulrQR3L1PeR2YQ6NqiE9qRESkGAxqRESkGAxqRESkGMypUeVll+fCnL+r5Hko+bHiZWo0su9u7pZgUctydW6Xa5FPiyWrr7iOLnky0Tg1N1NdFWy6GXvmUifm0ah645MaEREpBoMaEREpBpsfqdKQzyrv0sVf1P1fkH0dkzdHipvhBHmXfbWHvsvJz+uuidFNN/1iV6h21+TI5kYiCT6pERGRYjCoERGRYjCoERGRYjCnRmVLnvORT28lzhepivmOJS5rl5Z1ybG5m/rqdvJQ7vJkLmVLkTeT7nSzjzk0Inf4pEZERIrBoEZERIrBoEZERIrBnBpVHm6WpQEAlXgKK5cpqmQ5NvHEWfKpr4pbvsWdW8yFlWrsmcvBzKMRlRSf1IiISDEY1IiISDHY/Ejlq7gu/pKy7psjJadxWTVb5DZaG91x2y3ftbCbfWxeJPIUPqkREZFiMKgREZFiMKgREZFiMKdGFUucT3KXXwPcTqlVqvyWp5SmW77LscyjEZUFPqkREZFiMKgREZFiMKgREZFiMKdGlcdtjGFzq7glbW71vMWei3kzovLGJzUiIlIMBjUiIlIMNj9S5VWa5ki35/Fgk6LkvGxeJKps+KRGRESKwaBGRESK4fGgNnv2bLRv3x4BAQGoXbs2+vbti2PHjknK5OfnIyEhASEhIfD390f//v2Rlpbm6aoQEVE14/GgtnXrViQkJGD79u3YsGEDLBYLunfvjpycHEeZiRMn4ueff8bKlSuxdetWXLp0Cf369fN0VUhpBOHWfjx1nuLOS0QVTiUIZfvpvHr1KmrXro2tW7fi/vvvR2ZmJmrVqoXly5fjscceAwAcPXoUTZs2RVJSEu655x6Xc5hMJphMJse20WhEZGQkuqAPtCpdWVaflEDewYQBiajCWQULtuBHZGZmIjAw0GPnLfOcWmZmJgAgODgYALB7925YLBbExcU5yjRp0gRRUVFISkoq9ByzZ8+GwWBw/ERGRpZ1tYmIqAoq06Bmt9sxYcIEdOzYES1atAAApKamwsvLC0FBQZKyoaGhSE1NLfQ8kydPRmZmpuMnJSWlLKtNRERVVJmOU0tISMChQ4ewbdu22zqPXq+HXq/3UK2o2mFzI1G1UWZPamPHjsWaNWuwefNm1K1b1/F6WFgYzGYzMjIyJOXT0tIQFhZWVtUhIqJqwONBTRAEjB07FqtWrcKmTZsQExMj2d+2bVvodDps3LjR8dqxY8dw/vx5xMbGero6RERUjXi8+TEhIQHLly/Hjz/+iICAAEeezGAwwMfHBwaDASNHjkRiYiKCg4MRGBiIF154AbGxsYX2fCQiIiopjwe1hQsXAgC6dOkieX3x4sUYNmwYAOC9996DWq1G//79YTKZEB8fj48//tjTVSEiomqmzMeplQWj0QiDwcBxakREVVSVHadGRERUXhjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMRjUiIhIMSosqH300UeIjo6Gt7c3OnTogB07dlRUVYiISCEqJKh98803SExMxH/+8x/s2bMHrVu3Rnx8PK5cuVIR1SEiIoWokKA2b948jBo1CsOHD0ezZs2waNEi+Pr64osvvii0vMlkgtFolPwQERHJlXtQM5vN2L17N+Li4pyVUKsRFxeHpKSkQo+ZPXs2DAaD4ycyMrK8qktERFWItrwveO3aNdhsNoSGhkpeDw0NxdGjRws9ZvLkyUhMTHRsZ2ZmIioqClZYAKFMq0tERGXACgsAQBA8+0e83IPardDr9dDr9Y7tf5sft+HXiqoSERF5QFZWFgwGg8fOV+5BrWbNmtBoNEhLS5O8npaWhrCwsBKdIyIiAikpKRAEAVFRUUhJSUFgYGBZVLfKMxqNiIyM5D1yg/eoeLxHJcP7VLx/79H58+ehUqkQERHh0fOXe1Dz8vJC27ZtsXHjRvTt2xcAYLfbsXHjRowdO7ZE51Cr1ahbt67jiS0wMJD/gYrBe1Q83qPi8R6VDO9T8QwGQ5ncowppfkxMTMTQoUPRrl073H333Xj//feRk5OD4cOHV0R1iIhIISokqD3xxBO4evUqpk2bhtTUVNx5551Yt26dS+cRIiKi0qiwjiJjx44tcXNjUfR6Pf7zn/9IOpGQFO9R8XiPisd7VDK8T8Ur63ukEjzdn5KIiKiCcEJjIiJSDAY1IiJSDAY1IiJSDAY1IiJSDAY1IiJSjCob1LjIqNPs2bPRvn17BAQEoHbt2ujbty+OHTsmKZOfn4+EhASEhITA398f/fv3d5mqrDqZM2cOVCoVJkyY4HiN96jAxYsXMXjwYISEhMDHxwctW7bErl27HPsFQcC0adMQHh4OHx8fxMXF4cSJExVY4/Jls9kwdepUxMTEwMfHB/Xr18ebb74pmZi3ut2jP/74A71790ZERARUKhVWr14t2V+S+5Geno5BgwYhMDAQQUFBGDlyJLKzs0tfGaEKWrFiheDl5SV88cUXwj///COMGjVKCAoKEtLS0iq6ahUiPj5eWLx4sXDo0CFh3759Qs+ePYWoqCghOzvbUea5554TIiMjhY0bNwq7du0S7rnnHuHee++twFpXnB07dgjR0dFCq1athPHjxzte5z0ShPT0dKFevXrCsGHDhOTkZOH06dPC+vXrhZMnTzrKzJkzRzAYDMLq1auF/fv3C4888ogQExMj5OXlVWDNy8+sWbOEkJAQYc2aNcKZM2eElStXCv7+/sIHH3zgKFPd7tGvv/4qTJkyRfjhhx8EAMKqVask+0tyP3r06CG0bt1a2L59u/Dnn38KDRo0EJ588slS16VKBrW7775bSEhIcGzbbDYhIiJCmD17dgXWqvK4cuWKAEDYunWrIAiCkJGRIeh0OmHlypWOMkeOHBEACElJSRVVzQqRlZUlNGzYUNiwYYPQuXNnR1DjPSrw6quvCp06dSpyv91uF8LCwoR33nnH8VpGRoag1+uFr7/+ujyqWOF69eoljBgxQvJav379hEGDBgmCwHskD2oluR+HDx8WAAg7d+50lFm7dq2gUqmEixcvlur6Va758VYWGa1uMjMzAQDBwcEAgN27d8NisUjuWZMmTRAVFVXt7llCQgJ69eoluRcA79G/fvrpJ7Rr1w4DBgxA7dq10aZNG3z22WeO/WfOnEFqaqrkPhkMBnTo0KHa3Kd7770XGzduxPHjxwEA+/fvx7Zt2/DQQw8B4D2SK8n9SEpKQlBQENq1a+coExcXB7VajeTk5FJdr0qspyZ2K4uMVid2ux0TJkxAx44d0aJFCwBAamoqvLy8EBQUJCkbGhqK1NTUCqhlxVixYgX27NmDnTt3uuzjPSpw+vRpLFy4EImJiXjttdewc+dOjBs3Dl5eXhg6dKjjXhT2+asu92nSpEkwGo1o0qQJNBoNbDYbZs2ahUGDBgEA75FMSe5HamoqateuLdmv1WoRHBxc6ntW5YIauZeQkIBDhw5h27ZtFV2VSiUlJQXjx4/Hhg0b4O3tXdHVqbTsdjvatWuHt956CwDQpk0bHDp0CIsWLcLQoUMruHaVw7fffotly5Zh+fLlaN68Ofbt24cJEyYgIiKC96gSqHLNj55YZFSpxo4dizVr1mDz5s2oW7eu4/WwsDCYzWZkZGRIylene7Z7925cuXIFd911F7RaLbRaLbZu3Yr58+dDq9UiNDS02t8jAAgPD0ezZs0krzVt2hTnz58HAMe9qM6fv5dffhmTJk3CwIED0bJlSwwZMgQTJ07E7NmzAfAeyZXkfoSFheHKlSuS/VarFenp6aW+Z1UuqIkXGf3Xv4uMxsbGVmDNKo4gCBg7dixWrVqFTZs2ISYmRrK/bdu20Ol0knt27NgxnD9/vtrcs27duuHgwYPYt2+f46ddu3YYNGiQ4/fqfo8AoGPHji7DQY4fP4569eoBAGJiYhAWFia5T0ajEcnJydXmPuXm5kKtlv7p1Gg0sNvtAHiP5EpyP2JjY5GRkYHdu3c7ymzatAl2ux0dOnQo3QVvq5tLBVmxYoWg1+uFJUuWCIcPHxZGjx4tBAUFCampqRVdtQoxZswYwWAwCFu2bBEuX77s+MnNzXWUee6554SoqChh06ZNwq5du4TY2FghNja2Amtd8cS9HwWB90gQCoY7aLVaYdasWcKJEyeEZcuWCb6+vsJXX33lKDNnzhwhKChI+PHHH4UDBw4Iffr0UXR3dbmhQ4cKderUcXTp/+GHH4SaNWsKr7zyiqNMdbtHWVlZwt69e4W9e/cKAIR58+YJe/fuFc6dOycIQsnuR48ePYQ2bdoIycnJwrZt24SGDRtWny79giAICxYsEKKiogQvLy/h7rvvFrZv317RVaowAAr9Wbx4saNMXl6e8Pzzzws1atQQfH19hUcffVS4fPlyxVW6EpAHNd6jAj///LPQokULQa/XC02aNBE+/fRTyX673S5MnTpVCA0NFfR6vdCtWzfh2LFjFVTb8mc0GoXx48cLUVFRgre3t3DHHXcIU6ZMEUwmk6NMdbtHmzdvLvRv0NChQwVBKNn9uH79uvDkk08K/v7+QmBgoDB8+HAhKyur1HXhempERKQYVS6nRkREVBQGNSIiUgwGNSIiUgwGNSIiUgwGNSIiUgwGNSIiUgwGNSIiUgwGNSIiUgwGNSIiUgwGNSIiUgwGNSIiUoz/BxMaCPbvMmj7AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data)\n",
    "aper_hl.plot(color='w', label='Half-light ellipse')\n",
    "aper_ql.plot(color='r', label='Quarter-light ellipse')\n",
    "\n",
    "plt.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}