diff --git a/dev/notebooks/spergel_hlr_flux_radius_approx.ipynb b/dev/notebooks/spergel_hlr_flux_radius_approx.ipynb
new file mode 100644
index 00000000..e5bee6ab
--- /dev/null
+++ b/dev/notebooks/spergel_hlr_flux_radius_approx.ipynb
@@ -0,0 +1,805 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "d7b8bc37-8799-433c-9399-de95a21a1727",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import galsim\n",
+ "import numpy as np\n",
+ "\n",
+ "import ultraplot as uplot"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1cc6d839-6f5d-48aa-ba7a-dfd1a3371d33",
+ "metadata": {},
+ "source": [
+ "## Spergel HLR"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "850a3ad0-74c4-417b-b79d-2c5590f70027",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nu_vals = np.linspace(-0.85, 4.0, 500)\n",
+ "NU_FRAC = 0.5\n",
+ "\n",
+ "frac_vals = []\n",
+ "for nu in nu_vals:\n",
+ " prof = galsim.Spergel(nu, scale_radius=1)\n",
+ " frac_vals.append(prof.calculateFluxRadius(NU_FRAC))\n",
+ "\n",
+ "frac_vals = np.array(frac_vals)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "c781af78-6de8-474b-ab56-551455d981b4",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAI1CAYAAADcqV49AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAewgAAHsIBbtB1PgAAXndJREFUeJzt3XlcVOX+B/DPDAyDiKggA+KCgIK4U7jlbmpqehM3LDPLpMXyZpZRai65cK1upWWaaJlmVmiaW2n+LPGqqeASigqCCoowIIIsMsDM/P7wwuU4gyIMc85wPu/Xy9crv8/h8JVj8PGc5zyPIj8/3wgiIiIiAgAoxW6AiIiISEoYjoiIiIgqYDgiIiIiqoDhiIiIiKgChiMiIiKiChiOiIiIiCpgOCIiIiKqgOGIiIiIqAKGIyIiIqIKGI6IiIiIKrAXuwEx+fr6orCwEM2bNxe7FSIiIqqma9euwcnJCcnJyRY5n6zDUWFhIUpKSqBU8gZabdDr9QAAOzs7kTuhe/HaSBevjXTx2khXcXGxRc8n63DUvHlzKJVKnDt3TuxW6iStVgsA0Gg0IndC9+K1kS5eG+nitZGutm3bQqFQWOx8vGVCREREVAHDEREREVEFDEdEREREFTAcEREREVXAcERERERUgWzeVtNqtcjKyhLUdDodHBwcROqIiIiIpEg24SgyMhIREREmdXd3dxG6ISIiIqmSTTgKCwtDSEiIoBYaGso7R0RERCQgm3Ck0WhMFu5Sq9VcHZuIiIgEmAyIiIiIKmA4IiIiIqrAIuEoNjYWkydPhr+/Pxo3boxmzZph8ODBWLduXflGfURERES2oMbhKCoqCo8//ji2bt2KtLQ0lJSUIDc3F0ePHsUbb7yBMWPGoKSkpMrnmzdvHpydne/7a/LkyTVtm4iIiMisGoWjrKwsTJ8+HaWlpejatSv27duHtLQ0nD9/HgsXLoSDgwP279+PDz/8sMrnTExMrElLRERERDVSo7fVtmzZgvz8fHh5eWH37t1wcnICALi4uOCtt95CcXExlixZgvXr12POnDlVOuelS5cAAHv37kWvXr1q0h4RERHRQ6vRnaPY2FgAwMiRI8uDUUWjR48GANy4cQPZ2dkPPJ/BYMDly5cBAIGBgTVpjYiIiKhaahSOtFotAKBVq1Zmxxs0aFD+30aj8YHnS01NRVFRETw9PeHq6lqT1oiIiIiqpUaP1X755Zf7jh8+fBgA4OnpCTc3tweer+yRWkBAACIjI7FhwwacP38eDg4OaN++PSZPnoxnnnmGCzcSERFRrbH4CtlFRUVIT09HdHQ05s6dCwCYNWtWlT62bDL2oUOHcPDgQcE5jx49iqNHj2LXrl3YuHEjVCqVpVsnIiIiG3OnRI/sOyVwc7LcdmAWDUerV6/G22+/Xf57JycnfPLJJ3jppZeq9PFld44MBgOef/55vPHGG/D29sb169exdu1arFixArt27cLixYuxcOHCSs8THBxcpc+XnJwMHx+f8seDZFmZmZlit0CV4LWRLl4b6eK1kZZivQGbzmTis6PXkVlg2XBUq8+nCgsL8fvvv+PGjRtVOl6v1yMwMBDh4eH44osv0KZNGzg4OMDHxwdLliwpf+Nt5cqVyMrKqs3WiYiISIJKDUZ8/7cWj605g3d/v4L0/KqvpVhVivz8/AfPlH5I6enpOHbsGBYsWIDExES0adMGf/31F9RqdY3OW1RUhNatWyMnJwfr16/H2LFja3S+4OBgKJVKnDt3rkbnIfPK7sjdu+EviY/XRrp4baSL10ZceoMRP5y6jgX7EnApq0A4+O1raOteHzExMRb5XLVy58jT0xNPPfUUdu3ahYYNGyIxMRHbtm2r8XkdHR3LH5mVPYIjIiKiustgMGLr32no9PGfePb7U6bBqBbU6mO1Zs2aoW/fvgCAuLg4i5yz7K23oqIii5yPiIiIpMdoNGJXfAYe/TQaY7+NRXxGvtnj+vi6ooVLzZ5M3avaE7LT0tIQFBQE4O5ikM2bNzd7nLu7OwAgLy/vvufLzMzE6dOnoVKp0L9//0qPu337NgBUaWkAIiIisi1GoxH7E7Lw/m8XcCwlp9LjurVshEVDAzDY3x2Bn9tZtIdqhyMPDw8YjUYUFhbi0qVLlYajstfzmzVrdt/z5eXlISQkBMDdsBUQEGByTGlpafmq3GXBjIiIiOqG/yTfxJxfLyA6ufJdNTp7uWDR0ACMaOcBhUJRK31U+7GanZ1d+SOzb7/91uwxJ06cKF8I8vHHH7/v+Xx9fdGpUycAwNKlS80es3btWmi1WrRs2RI9e/asbutEREQkIaeu5WJ45DH0WXmk0mDUVuOMn557FCff7IuR7T1rLRgBNZxzNG3aNABAVFQUXnjhBZw5cwYFBQVITU3FN998gzFjxsBgMGDo0KGCtYeCgoIQFBSE+fPnC85Xtljk1q1bERYWhvPnz0On0yElJQVLlixBeHg4AGDRokWws7PsLTQiIiKyrovafIRuiMUjn0bj1wvm1xz0c3PCxmeCcHZWf4zr7AWlsvZCUZkaLQI5cOBAvP322/j4448RFRWFqKgok2O6deuGyMhIQa3sUVt6erqgHhISgtdeew0rV67E5s2bsXnzZpPzzZ8/H2PGjKlJ20RERCSilFuFWLgvAetPpMJQyYJCLRo5Yt5gf0zu2gIqO+tuG1bjFbIXLFiA3r17Y82aNYiJiUF2djacnZ3Rrl07hIaG4rnnnnuorT6WLVuG/v37IzIyErGxscjNzYWrqysee+wxvP766+jRo0dNWyYiIiIRaPN0WPp/iVh15CqK9Qazx2icHTB3kD9e6tkSantxnhJZZPuQQYMGYdCgQVU+Pj/f/Ot4ZYYNG4Zhw4bVtC0iIiKSgJw7Jfj4zyR8Fp2MgmK92WMa1VPhnQF++GdvH9RXW3zr14ci7mcnIiKiOqtAV4ovDl/BsgOXcOuO+W0+nBzs8EYfH8zq74fGFtwfrSZkE460Wq3Jfmw6nQ4ODtK4EERERHVFcakBkX9dxeL9iUjP05k9RmWnwCs9W2H2463h6eJo5Q7vTzbhKDIyEhERESb1skUqiYiIqGb0BiO+i72GBfsu4kr2HbPHKBXA5OAWmDfEH61cnazcYdXIJhyFhYWVLzJZJjQ0lHeOiIiIashoNGJbXDrm/nYB5yvZ5gMAxnZqikVDA9DWo4EVu3t4sglHGo3GZCdltVoNpdK6rwcSERHVFUajEb8nZGL2nguIvZZb6XFD27pj8dC2eLRFI+s1VwOyCUdERERkOX9dvYX3dp/Hn0k3Kz2mV6vGWDo8EH39bGs/VIYjIiIiqrLzGXmY8+sFbItLr/SYLl4uWDq8LYa21dTqNh+1heGIiIiIHuhazh0s2JuAb06kVLqqtb97fSwa2hZjOzW1yjYftYXhiIiIiCp1q7AY/zpwCSsOXUZRqflVrVs0csSCIQF4Lrg57K281UdtYDgiIiIiE3dK9Pj80GVEHLiEnEoWcHRzUmHOoDZ49bFWcFTVnQ3hGY6IiIioXKnegPUnUrFgXwKu5xaZPcbJwQ4z+/ri7f5+aFiv6vun2gqGIyIiIoLRaMT2s+mYvecCLmjNr1Vkr1QgrEdLvD/YH00ltqq1JTEcERERydzBpCy8u/sC/rp6q9JjQrt4YfGwtmjdpL4VOxMHwxEREZFM/Z12G+/tOY8957WVHjOoTRP868lAm1nA0RIYjoiIiGTmSnYh3v/tAjadvA5jJa/lP9K8If41PBCDA+S3B6lswpFWq0VWVpagptPpuLcaERHJRma+Dkv2J2LVkaso1pt/Ld/PzQlLhrXFuM5eNr1WUU3IJhxFRkYiIiLCpO7uLr9ETERE8lJYXIpPo5Ox7EAS8nSlZo/xaKDGvMH+mNq9JRzsbX+topqQTTgKCwtDSEiIoBYaGso7R0REVGfpDUZ8eyIV7/92EWm3zb+W30Btj3cG+GFGX184q2UTC+5LNl8FjUYDjUYjqKnVaiiV8k7HRERU9xiNRvx2QYt3dp3H2fQ8s8c42CkxrZc3Zj/eBu7Oait3KG2yCUdERERycPJaDmbtPI8Dl7LMjisUwMRHmmHR0LZo5epk5e5sA8MRERFRHXA1uxBzf7uA72KvV3rMYP8m+HBEO3Rp1tCKndkehiMiIiIbdquwGBH/dwkr/nMZuko2hu3U1AUfjQzEkACN2XESYjgiIiKyQbpSPb48fAWL9yciu9D8xrDNGzpi8bC2ePbR5rCT6Wv51cFwREREZEOMRiN+PJ2G2Xsu4HJ2odljXBzt8d7A1nijry/qqeys3KHtYzgiIiKyEQeTsjBr53mcSM0xO26vVGBar1aYO4hvoNUEwxEREZHExafn4d3d57EzPqPSY8Z2aoqIJwNlsTFsbWM4IiIikqj020WYv/ci1h5LgaGSPdB6tWqMj//RHj28G1u3uTqM4YiIiEhiCnSl+PfBZHz4xyUUFOvNHuPvXh/LngzEUx08oVBwsrUlySYcceNZIiKSOoPBiO9OXsPsPRdwPdf8dh8aZwcseCIAU7u3hMqOuzzUBtmEI248S0REUnYwKQszd8Tj5LVcs+P1VEq81c8P7wxojQaOsvnxLQrZfHW58SwREUnRpawCvLMrHtvi0s2OKxXAC11b4oOhAfBq6Gjl7uRJNuGIG88SEZGU3CosxqLfE/HF4cso0ZufbT2oTRP8+x/t0cnLxcrdyZtswhEREZEUlOgNWHXkChbuS6h0Zeu2Gmf8+x/tMKythpOtRcBwREREZAVGoxE7z2Vg1q54JGQWmD2mSX0HLHwiAGE9ONlaTAxHREREtezUtVy8tfMc/rh00+y4g50Sb/TxwexBbdConsrK3dG9GI6IiIhqSVpuEeb+egHrY1JhrGQRx3Gdm+JfTwbC140rW0sFwxEREZGFlS3iuOyPSyisZBHHri0a4dOn2qOXj6uVu6MHYTgiIiKykKos4tiikSP+9WQgJnRpBqWSk62liOGIiIjIAo5czsYbv5xFTKr5RRyd1XZ4b2AbvNnPF/VUdlbujh4GwxEREVENpN66g/Dd57H51HWz40oF8GL3lvjgiQB4unARR1vAcERERFQNhcWl+OiPJCz74xLulBjMHsNFHG0TwxEREdFDMBqN+PF0Gt7ZFY/UHPPzigLc6+Pf/2iP4YFcxNEWySYcabVaZGVlCWo6nY57qxERUZWdSS/Awh8TcPjKLbPjjeqpsGCIP6b1asVFHG2YbMJRZGQkIiIiTOru7u4idENERLYk/XYR3tyTjB/jMmFuuSKlAni5pzc+eCIATZzVVu+PLEs24SgsLAwhISGCWmhoKO8cERFRpXSlenwWfRmL9ycgX2d+vaKBrZvgs1Ht0bEp5xXVFbIJRxqNBhqNRlBTq9VQKnnbk4iIhIxGI345m463dsYj+Wah2WN83Zzw75Ht8FQHT84rqmNkE46IiIiqIu7GbczYfg4HLmWZHXdW22HuIH/M6OsDtT3XK6qLJBmOYmNjsWLFChw9ehSZmZlwcnJCu3btMGHCBDz//POws+NfRiIisqysfB3m7b2Ir45ehcHMxCIFgNCO7vhkdBc05XpFdZrkwlFUVBTCwsJQWlpaXsvNzcXRo0dx9OhR7Ny5E1FRUVCpuGsxERHVXKnegFVHrmLe3ovIuVNi9pherRpjXl8vdGnqDA2DUZ0nqQk3WVlZmD59OkpLS9G1a1fs27cPaWlpOH/+PBYuXAgHBwfs378fH374oditEhFRHfDnpSwEfRKNf24/azYYNW/oiM3PPoJDr/dCl6bOInRIYpDUnaMtW7YgPz8fXl5e2L17N5ycnAAALi4ueOutt1BcXIwlS5Zg/fr1mDNnjsjdEhGRrUq9dQdv74zHT2fSzI7XUykRPqA1Zg3wg5ODpH5UkhVI6s5RbGwsAGDkyJHlwaii0aNHAwBu3LiB7Oxsq/ZGRES2r6hEjyX7E9D2wz8qDUZPBzXDxfCBmP9EAIORTEnqqmu1WgBAq1atzI43aNCg/L+NRnPLcBEREZm3Kz4DM7afRVIlr+Z38XLB5yEd0NvXzcqdkdRIKhz98ssv9x0/fPgwAMDT0xNubvzLS0RED5aYmY8Zv5zDnvNas+ON66mwZHhbvNTDG3ZKrldEEgtH5hQVFSE9PR3R0dGYO3cuAGDWrFkid0VERFKXryvFkv2J+ORgMor1BpNxhQJ4qYc3Fg/llh8kJOlwtHr1arz99tvlv3dycsInn3yCl1566b4fFxwcXKXzJycnw8fHp/xxHllWZmam2C1QJXhtpIvXpuaMRiO2nb+JD/5IwY1886/md23mjKWDWqGTZ30YCnOhNf+kTYDXRrpKS0stusSPpMPRvQoLC/H7779j5MiRaNq0qdjtEBGRxJzTFmLO/is4mppndlxTX4V5/VtgbPsm3PKDKqXIz8+X/Mzm9PR0HDt2DAsWLEBiYiLatGmDv/76C2p1zW6DBgcHQ6lU4ty5cxbqlCoquyN37552JD5eG+nitame7MJizPvtIlYduWJ2dWt7pQIz+vri/cFt4OJYvTsMvDbS1bZtWygUCsTExFjkfDZx58jT0xNPPfUUgoOD0a1bNyQmJmLbtm2YMGGC2K0REZGI9AYj1h1Lwew953Gz0PwjtCH+7lg+qj3aejQwO050L0mtc/QgzZo1Q9++fQEAcXFxIndDRERiOpGSgx4rDuHlLX+bDUatXOth2/PB+O2l7gxG9FAkc+coLS0NQUFBAO4uBtm8eXOzx7m7uwMA8vLMP08mIqK6LbuwGLP3XMCav67C3JJ3jvZKvPd4G8wa4Id6Km5UTg9PMuHIw8MDRqMRhYWFuHTpUqXhKDExEcDdu0hERCQfBoMR60+kInz3eWQVFJs9Zkynpvj3yHbwdjXdZYGoqiTzWM3Ozq78kdm3335r9pgTJ06ULwT5+OOPW603IiIS1+nruej9xWG8+NMZs8GorcYZv7/cA1smBzMYUY1J5s4RAEybNg2//fYboqKiAAAzZsxA69atkZ2djf3792P+/PkwGAwYOnRoldcyIiIi25V7pwTv/3YRKw9fNvsWmpODHeYN9sebfX3hYC+Zf++TjZNUOBo4cCDefvttfPzxx4iKiioPSRV169YNkZGRInRHRETWYjQasenkdby9Mx4ZeTqzx4zu6IlPn2qPlo15p4gsS1LhCAAWLFiA3r17Y82aNYiJiUF2djacnZ3Rrl07hIaG4rnnnrPoKphERCQtZ2/cxms/xyE6OdvsuJ+bE74Y3RFD23K9IaodkgtHADBo0CAMGjRI7DaIiMiK8opKsXDfRXx26DL0Zp6hOdorMXtQG8zq7wdHvoVGtUiS4YiIiOTDaDQi6swNvPnLOaTdLjJ7zIh2Hlg+qj183epbuTuSI9mEI61Wi6ysLEFNp9PBwcFBpI6IiOiiNh+v/xyH/YlZZsdbudbDilEdMLK9p5U7IzmTTTiKjIxERESESb1sUUkiIrKeOyV6LNmfiA//uIQSvekjNAc7Jd4Z4If3Hm8NJwfZ/KgiiZDN37iwsDCEhIQIaqGhobxzRERkZb9d0OK1n+OQfLPQ7PgQf3d8ProD/N2drdwZ0V2yCUcajcZkJ2W1Wg2lkutiEBFZw43bRZix/Rx+OpNmdrx5Q0d8Nqo9RndsCoVCYeXuiP5HNuGIiIjEoTcYsfrIFcz+9QJuF5WajNsrFZjZzxfvD/aHs5o/lkh8/FtIRES15tS1XLy85W+cSM0xO96rVWOsHtsJHZq6WLcxovtgOCIiIovLKyrFvL0XsOKQ+W0/GtdT4aOR7fBC1xZQKvkIjaSF4YiIiCzGaDRi+9l0/HPbWVzLNb9m0XPBzfHxyHZwd1ZbuTuiqmE4IiIii7iaXYjp285iZ3yG2fEA9/pYNbYTBrRuYuXOiB4OwxEREdVIid6Az6KTsWBfAgqL9Sbjansl5gxqg3cG+EFtz20/SPoYjoiIqNqOXsnGy1v+RtyNPLPjg9o0wZdjOqIN1ywiG8JwRERED+1WYTHe23MBXx29anZc4+yAz57qgAlBXlyziGwOwxEREVVZ2Sax/9x+Fhl5OpNxhQJ4uYc3lg5vi8ZO3IGAbBPDERERVUnqrTuY9nMcdlUy4bpTUxd8Na4Teng3tnJnRJYlm3Ck1WqRlSXc9Vmn03FvNSKiB9AbjPjy8BXM/vU88nWmE66dHOzwwRMBeKOPD+ztuCUT2T7ZhKPIyEhERESY1N3d3UXohojINsTduI2wn87gWEqO2fEnAzX4ckxHtGzsZN3GiGqRbMJRWFgYQkJCBLXQ0FDeOSIiMqOoRI9Fvyfgwz+SUGpmiWuPBmqsGNUB4zpzk1iqe2QTjjQaDTQajaCmVquhVPIWMBFRRX9eysJLUX8jMavA7PjU7i3x4YhATrimOks24YiIiO7vVmExZu08j3XHU8yOt2lSH2vGdUJ/rnBNdRzDERGRzD3o9Xx7pQLhA1tj7qA2cFRxhWuq+xiOiIhkLOVWIV77+Wylr+d3b9kIkeM7o2NTFyt3RiQehiMiIhnSG4xYefgy5vx6wezr+c5qOywdFohpvVrBTskJ1yQvDEdERDITn56HF386g7+u3jI7PqKdB74c3REtGtezcmdE0sBwREQkEyV6A5YduIRFvyeiWG8wGfdooMbnIR0wthNfzyd5YzgiIpKBU9dyMeXH0ziddtvsOF/PJ/ofhiMiojqsqESPD/67mKPezGKOrZvURyRfzycSYDgiIqqjjlzOxos/ncEFbb7JmFIBvNXPDwuHBqAeX88nEpBNOOLGs0QkFwW6Usz+9QI+/89lGE1vFqGDZwN8HdoFXVs2snpvRLZANuGIG88SkRz8X0ImwqL+xuXsQpMxlZ0Ccx5vg/cebwMHe26dRFQZ2YQjbjxLRHVZ7p0SvL0zHmuPmd/6I7hFQ3wd2oWLORJVgWzCETeeJaK6aue5dLyyJQ5pt4tMxhztlfhgaADe7OsLezt+vyOqCtmEIyKiuiYrX4c3tp/D96eumx3v4+uKteM7w9/d2cqdEdk2hiMiIhv00+k0vL4tDpn5xSZjzmo7LHuyHV7p6Q0lt/4gemgMR0RENiQzX4dpW+Ow5e8bZsefCHDHV2M7wdvVycqdEdUdDEdERDZiy5k0TPvZ/N2iRvVU+Oyp9nguuDm3/iCqIYYjIiKJy8rX4fVtZ/Hj6TSz4yEdPbFydEc0dXG0cmdEdRPDERGRhG2Lu4FXtvwNrZm7RW5OKqwc3RHju3jxbhGRBTEcERFJ0M2CYvxz29lK30QL6eiJVWM6waOB2sqdEdV9DEdERBLzW+IthP9+Gul5OpOxxvVU+GJ0Bzwd1Ix3i4hqCcMREZFE3Cosxuu7khB1Lsvs+D/ae2D12E6cW0RUyxiOiIgkYFd8Bl6KOoMbt03vFjWqp8LnIR0w8RHeLSKyBtmEI61Wi6ws4b/GdDod91YjIlHl3CnBjO1n8W3MNbPjI9p54KuxneDVkHeLiKxFNuEoMjISERERJnV3d3cRuiEiAn67oMXUn87geq7pnmgNHe2xfFQHrltEJALZhKOwsDCEhIQIaqGhobxzRERWl68rxayd8Vh99KrZ8YE+DfHts13RvFE9K3dGRICMwpFGo4FGoxHU1Go1lEruUk1E1nPkcjae23wKSTcLTcZcHO2xcEALPN3RHR4MRkSikU04IiISU3GpAfP3XsSHf1yCwWg6PsTfHWvHd4a6JM/6zRGRAMMREVEti7txG5O+P4UzabdNxpwc7PDvke3wck9vKBQKaLUMR0Ris8gzpYyMDMydOxfBwcHlj6+6deuGRYsWIScn56HONW/ePDg7O9/31+TJky3RNhFRrdIbjPjwwCUEf3rIbDB6rFVjnHmrH155rBUnXRNJSI3vHMXHx2PkyJHIyMgwqcfHx+OHH37AL7/8gtatW1fpfImJiTVtiYhIdMk3CzB582n853K2yZjKToEPngjArAGtYadkKCKSmhqFI6PRiBdeeAEZGRnw8/PDRx99hN69eyMvLw8HDhzAnDlzcPXqVYwfPx7Hjh2DSqV64DkvXboEANi7dy969epVk/aIiKzOaDRi3bEUvLnjHPJ1epPxjk0bYOMzQejs1VCE7oioKmoUjg4cOIBz585BpVJh27Zt8PX1BQA4OTnh6aefRvfu3dGjRw8kJCRg+/btGDdu3H3PZzAYcPnyZQBAYGBgTVojIrK69NtFmPrTGew+rzUZUyiAWf398MHQAKjt7UTojoiqqkZzjv744w8AQP/+/cuDUUW+vr7lawsdPnz4gedLTU1FUVERPD094erqWpPWiIisasuZNHT46E+zwcjH1QnR0x7DshHtGIyIbECNwlFSUhKA+9/lKVtbqKCg4IHnK3ukFhAQgMjISPTp0wdNmjSBl5cXBg8ejO+++w4Gg6EmLRMRWVTOnRJM+v4kxm2Ixc3CEpPxsB4tceatfujt6yZCd0RUHTV6rPbqq69i9OjRaNeuXaXHnD59GgDg7e39wPOVTcY+dOgQDh48WF4vKirC0aNHcfToUezatQsbN26s0vwlIqLa9OelLDy3+RRSc0y3//BooMa68Z3xZDsPETojopqoUTjq27fvfcf3799f/uhtxIgRDzxf2Z0jg8GA559/Hm+88Qa8vb1x/fp1rF27FitWrMCuXbuwePFiLFy4sNLzBAcHV6n/5ORk+Pj4QKs1vQ1ONZeZmSl2C1QJXpuaKdYbsOzQNaw8dgNm1nPEyABXLBvSCm5Oiof+/sJrI128NtJVWlpq0ZsmtbZ3xqZNmzBx4kQAwJgxY9ClS5cHfoxer0dgYCDCw8PxxRdfoE2bNnBwcICPjw+WLFmCOXPmAABWrlyJrKys2mqdiKhSiTfv4MmN5/CFmWDkorbDyhF+iHyqNdyceHebyFYp8vPzzf3Dp9ri4uIQHh6O6OhoAEDv3r2xbds21KtX832CioqK0Lp1a+Tk5GD9+vUYO3Zsjc4XHBwMpVKJc+fO1bg3MlX2L+Z797Qj8fHaPDyj0YjVR6/irR3ncKfEdO7jwNZNsH5CF7RoXLPvdbw20sVrI11t27aFQqFATEyMRc5nsTtHubm5ePPNN9GrVy9ER0dDpVJhzpw52LVrl0WCEQA4OjqWPzIrewRHRFTbtHk6/OPrE5i2Nc4kGDnYKfHxyHb4/eUeNQ5GRCQNFtlbLSYmBpMmTUJqaioAYOTIkVi8eDH8/PwscXoBN7e7b3wUFZlOgCQisrQ95zPwwg+noc0vNhlr5+GM7599hAs6EtUxNQ5HBw8exLhx41BYWAhvb2+sWrXqgRO1zcnMzMTp06ehUqnQv3//So+7ffvu/kRlIYmIqDbcKdFj1s54rDx8xez4671a4cOR7VBPxXWLiOqaGoWjmzdvYtKkSSgsLES/fv2wefNmuLi4VOtceXl55QtGxsbGIiAgwOSY0tJSxMbGAgCCgoKq3zgR0X2cvp6LZzadxPmMfJMxjwZqfBPaGcMC+Yo+UV1VozlHa9asQXZ2Nry9vREVFVXtYATcXU27U6dOAIClS5eaPWbt2rXQarVo2bIlevbsWe3PRURkjsFgxMd/JKHb8kNmg9HIdh74+61+DEZEdVyNwtHu3bsBAFOmTIGTk1OVPy4oKAhBQUGYP3++oD5r1iwAwNatWxEWFobz589Dp9MhJSUFS5YsQXh4OABg0aJFsLPjrWwispxrOXcw+Ku/MGtXPEr0wpd466mUWD22I36Z0hWaBmqROiQia6n2Y7WSkpLyV+Dnz59vEnTuNW3aNHz44YcA/rcSdnp6uuCYkJAQvPbaa1i5ciU2b96MzZs3m5xn/vz5GDNmTHXbJiIysS3uBl788Qxu3THd/uOR5g2x6ZkgtPVoIEJnRCSGaocjrVaLkhLTbyQ1tWzZMvTv3x+RkZGIjY1Fbm4uXF1d8dhjj+H1119Hjx49LP45iUie7pTo8daOc1h15KrJmEIBhA9ojYVPBMDBvtbWyyUiCap2OGrWrBny802fyVfFgz5u2LBhGDZsWLXOTURUFefS8zBhYyzOpueZjLVo5IiNzwShn18TETojIrFZZJ0jIiJbYTQaseavq5ix/RyKSk1Xuh7f2QtfjeuERvW4/QeRXDEcEZFs3CosRljU39j69w2TMScHO6wY1QFTurWAQqEQoTsikgrZhCOtVmuyWa1Op4ODg4NIHRGRNf0n+Sae2XQSqTmmq+t3auqCHyY9gkBOuiYiyCgcRUZGIiIiwqTu7u4uQjdEZC16gxFL9idi4b6LMJjZZnt6bx98OCIQjlzpmoj+SzbhKCwsrHwF7jKhoaG8c0RUh13LuYOJm04iOjnbZMzNSYVvJnTByPaeInRGRFImm3Ck0Wig0WgENbVaDaWSr+gS1UXb427gxZ/OILvQdMmR/n5u+G5iEJo1rCdCZ0QkdbIJR0QkD3dK9Hh7Rzy+PHLFZMxOqcDCJ/zx7sA2sFNy0jURmcdwRER1xoWMPIzfGIu4G6ZrF7VsXA+bJz6Cx3xcReiMiGwJwxER1QmbYq/h5S1/o6BYbzI2tlNTrBnXCY2dOMeQiB6M4YiIbFphcSn+ue0c1h1PMRmrp1Ji+agOmNq9JdcuIqIqYzgiIpt1ISMP4zaY3wKkg2cD/DjpUbTz5NpFRPRwGI6IyCZtjEnFq1vjzD5Gm9q9JZaPag8nB36LI6KHx+8cRGRTCotLMX3bWXx9PNVkrL6DHVaP7YRnH20uQmdEVFcwHBGRzYhPv/s22jkzj9E6Nm2AnyY9irbcAoSIaojhiIhswob/PkYrNPMYLaxHSywf1QH1uAUIEVmAbMIRN54lsk0FulK8vu0s1p8wfYzmrLbDV2M74ZlH+BiNiCxHNuGIG88S2Z749DyM2xCD+Ix8k7FOTV3w03OPIkDjLEJnRFSXySYcceNZItuyKfYaXtryt9nHaC/39ManT7XnYzQiqhWyCUfceJbINuhK9Xjzl3NYdeSqyZiz2g5rxnbG0480E6EzIpIL2YQjIpK+q9mFGLchFidSc0zGOnvdfYzm787HaERUuxiOiEgSfrugxcRNJ5FdWGIyFtajJVaM6gBHPkYjIitgOCIiUekNRnywLwGL9ifAaBSO1VMpsWpMJ0zu2kKc5ohIlhiOiEg0Wfk6TNx0CvsSMk3GWjepj62Tg9HJy0WEzohIzhiOiEgUx67ewrgNMUjNKTIZC+noiW9Cu6BhPZUInRGR3DEcEZFVGY1GrDx8BTN3nEOJXvgczU6pwLInAzGzny8UCoVIHRKR3DEcEZHV5OtKEfbTGfxwOs1krKmLGj9OehR9fN1E6IyI6H8YjojIKs5n5GHMtzE4b2a1635+bvjh2Ufg6eIoQmdEREIMR0RU67b+nYbnfziNfJ3patfhA1pj8bAA2NtxQVYikgaGIyKqNXqDEXN/vYB/HbhkMtbQ0R4bng7CPzp4itAZEVHlZBOOtFotsrKyBDWdTse91Yhqyc2CYjzz3Umzr+l38XLBlsnB8GtSX4TOiIjuTzbhKDIyEhERESZ1d3d3EbohqttOX89FyPoTuJJ9x2RscnBzrBrbiZvGEpFkySYchYWFISQkRFALDQ3lnSMiC9sUew1hUWdwp8QgqNsrFVg+qgNefcybr+kTkaTJJhxpNBpoNBpBTa1WQ6nkJFAiSyjRGzBrZzyWH7psMubZQI0tk4PRy8dVhM6IiB6ObMIREdWejDwdQjfG4mDSTZOxnt6NsWVyMLwa8jV9IrINDEdEVCPHU25h9PoYXM813QbklZ7eWD6qAxzseYeWiGwHwxERVdu6YymYtjUOxXrh/CK1vRJfju6IKd1bitQZEVH1MRwR0UMrLjXgn9vP4qujV03GWjRyxNbJXdG1ZSPrN0ZEZAEMR0T0ULR5Ooz5Ngb/uZxtMtbfzw0/TnoUmgZqETojIrIMhiMiqrJT13Lx1DfHkZpjOr/ozb6++HBEILcBISKbx3BERFXy0+k0PP/DKZP1i+qplFg7vjOeeaS5SJ0REVkWwxER3ZfBYMS8vRexZH+iyViLRo745YVuCGreUITOiIhqB8MREVXqdlEJJn1/CjvOZZiM9fZxxZbJwfDg/CIiqmNkE4648SzRw0nKKsA/vj6O+Ix8k7GwHi3xRUhHrl9ERHWSbMIRN54lqrr/S8jEuA2xuHWnRFC3Uyqw/Kn2mNarFfdHI6I6SzbhiBvPEj2Y0WjE5/+5jJk74qE3GAVjbk4qRE0OxoDWTUTqjojIOmQTjrjxLNH96Ur1mLY1Dl8fTzUZ6+DZADumdIOPm5MInRERWZdswhERVU6bp0PI+hM4cuWWydioDp7Y8HQQGjjy2wURyQO/2xHJ3NkbtzFi3XFcvXXHZGzeYH/MH+IPpZLzi4hIPhiOiGRsz/kMTNh4Enm6UkHdycEOG57ugjGdvETqjIhIPJIMRxkZGfj888/x22+/ISUlBQDQqlUrjBw5EtOnT0ejRo3EbZDIxhmNRqyJSceCP1Jwz7xrtGxcDzumdEVnLy7sSETyJLlwFB8fj5EjRyIjI8OkHh8fjx9++AG//PILWrduLVKHRLatRG/AO/uuYMNprclYD+/G2P5CVy7sSESyJqlXtYxGI1544QVkZGTAz88PP//8M7RaLZKSkhAZGQmNRoOrV69i/PjxKCkpefAJiUjgVmExhkUeMxuMng5qhj9e7clgRESyJ6lwdODAAZw7dw4qlQrbtm3DkCFD4OTkBA8PDzz99NPYv38/nJyckJCQgO3bt4vdLpFNSczMR48V/8H/JWaZjH0wNACbJgbBUWUnQmdERNIiqXD0xx9/AAD69+8PX19fk3FfX9/yhRwPHz5s1d6IbNmfl7LQffl/kJBZIKg72ivx46RH8f5gf654TUT0X5Kac5SUlAQACAwMrPSYsoUcCwoKKj2GiP5n3bEUvLLlb5TeM/NaU1+FnVO7o1vLxiJ1RkQkTZIKR6+++ipGjx6Ndu3aVXrM6dOnAQDe3t6VHhMcHFylz5ecnAwfHx9otabzL6jmMjMzxW5B1gxGIxb9mYovj98wGQtwdcAXgzzRyrGEf/8lhv/fSBevjXSVlpZCpVJZ7HySCkd9+/a97/j+/fvLH72NGDHCGi0R2aTCEj1e35WE3QmmK14/0boRFvVsBCeVpJ6qExFJhqTC0f1s2rQJM2fOBACMGTMGXbp0qfTYmJiYKp0zODgYSqXSZM81six+fa1Lm6fDhM3HcSwlx2RsVn8/RDwZiJtZd/8FzGsjXbw20sVrIz329paNM5IPR3FxcQgPD0d0dDQAoHfv3li9erXIXRFJ04WMPAxfexyXswsFdXulAl+N7YQp3VuK1BkRke2QbDjKzc3FggULsG7dOhgMBqhUKrzzzjuYNWuWxRMiUV0QnXQTo745gVt3hGuANXS0x8/Pd8XANk1E6oyIyLZIMmXExMRg0qRJSE1NBQCMHDkSixcvhp+fn8idEUnT9yev4YUfzqBYbxDUvRvXw56p3dHOs4FInRER2R7JhaODBw9i3LhxKCwshLe3N1atWvXAidpEcmU0GrH0/xIx99eLJmPBLRpi55Ru8HRxFKEzIiLbJalwdPPmTUyaNAmFhYXo168fNm/eDBcXF7HbIpKkEr0Br2z5G18fTzUZ+0d7D3w/8RHUV0vqf3EiIpsgqe+ca9asQXZ2Nry9vREVFQUnJyexWyKSpNw7JRi3IQa/J5huBTK9tw8+fao97JRc8ZqIqDokFY52794NAJgyZQqDEVElUm/dwZPrjiHuRp6grlAAn/yjPWb0Nd16h4iIqk4yq8CVlJTg3LlzAID58+fD2dn5vr/eeecdkTsmsr4zabnovuKQSTCqp1Ji6+RgBiMiIguQzJ0jrVaLkpKSBx9IJFMHErMQsv4EbheVCuruzg7YOaUbuntzjzQiIkuQTDhq1qwZ8vPzxW6DSJJ+OHUdz20+hRK9cPPYthpn7J7aDb5u9UXqjIio7pFMOCIi8z49mISZO+JN6n18XfHLC13R2MlBhK6IiOou2YQjrVaLrCzhmz06nQ4ODvzBQtJkMBgRvvs8Pv4zyWRsTKem+O6ZIDiq7ETojIiobpNNOIqMjERERIRJ3d3dXYRuiO6vuNSAKT+exqaT103GXuvVCstHdeCr+kREtUQ24SgsLAwhISGCWmhoKO8ckeTkFZVi9PoT2J9ouobR0uFt8e7A1lAoGIyIiGqLbMKRRqOBRqMR1NRqNZRKyaxmQIT020UYvvYYTl2/LajbKRVYN74zJndtIVJnRETyIZtwRCR1CZn5GLrmGC5nFwrq9R3ssGVyMIa21VTykUREZEkMR0QScOzqLYxYdxxZBcWCuruzA3a/2B1dWzYSpzEiIhliOCIS2d4LWoz+NgaFxXpB3dfNCXtf6oHWTbiGERGRNTEcEYnox1PXMcnM4o6PNG+IPVO7w6OBWqTOiIjki+GISCSrj1zBtJ/jYBTmIgzxd8eWycFo4Mj/PYmIxMDvvkRWZjQasfT/EjH314smY88ENcM3E7rAwZ5vURIRiYXhiMiKDAYj3t4Zj0+jk03GXv/v4o5KLu5IRCQqhiMiKynVGzD1pzP4Nuaaydj8If6YP8SfizsSEUkAwxGRFRSV6DFhYyx+OZdhMrZ8VHv8s4+vCF0REZE5sglH3HiWxHK7qARPfX0CfybdFNTtlAqsn9AFzz7aXKTOiIjIHNmEI248S2LIzNdhWOQxxF7LFdQd7ZWImhyMEe08ROqMiIgqI5twxI1nydqu5dzBoNVHcTGzQFB3cbTHzind0NfPTaTOiIjofmQTjrjxLFlT8s0CPL76KK5k3xHUNc4O+C2sB4KaNxSpMyIiehDZhCMia7mQkYfHV/+FtNtFgrp343r4/eUeaOPuLFJnRERUFQxHRBZ0Ji0Xg7/6C5n5wg1k22qcsf+VHmjWsJ5InRERUVUxHBFZyPGUW3hizTHk3CkR1Lt4uWDfyz3g7sx90oiIbAHDEZEFRCfdxIh1x5GnKxXUu7dshF/DuqOxEyf+ExHZCoYjohrad1GLUd+cwJ0Sg6De19cVu17szg1kiYhsDL9rE9XAL2fTMX5DLIr1wmD0RIA7fn4+GE4O/F+MiMjW8Ds3UTX9cOo6nv3+FPQGo6A+qoMnfpj0CNT2diJ1RkRENcFwRFQNXx9LwdSoMzAKcxGeDmqGb5/uApUd188iIrJVDEdED2n1kSt4dWucSf3Fbi3x1bhOsFMqROiKiIgshf+8JXoIK/9z2Www+mcfH6xhMCIiqhNkc+dIq9UiKytLUNPpdNxbjapsxaFkvLH9nEn9vcdbY8mwtlAoGIyIiOoC2YSjyMhIREREmNTd3d1F6IZszacHkzBzR7xJfcEQf8x/IkCEjoiIqLbIJhyFhYUhJCREUAsNDeWdI3qgj/9IwqxdpsFo0dAAzB3sL0JHRERUm2QTjjQaDTQajaCmVquhVHLaFVVu2YFLeHf3eZN6xPC2ePfxNiJ0REREtU024YjoYS3dn4g5v14wqX84IhCzBrQWoSMiIrIGhiMiMxb9noB5v100qf/7H+0ws5+fCB0REZG1MBwR3WPh3otYsC/BpP7pU+0xo6+vCB0REZE1MRwRVfDBvgSzwWjFqA6Y3sdHhI6IiMjaGI6I/ivi/xIxf6/po7SVoztiWq9W1m+IiIhEwXBEhLuv68/eYzr5evXYjni5ZyvrN0RERKLhe+wke8ujk82uY7RqDIMREZEcMRyRrH15+Apm/GK6JciKUR3wymOtrN8QERGJjuGIZCvyr6t47WfTTWT//Y92nHxNRCRjsplzxI1nqaL1x1Px8pa/Ter/ejKQ6xgREcmcbMIRN56lMt/FXsOUn07DaBTWPxgagPCBXPmaiEjuZBOOuPEsAcCPp65j8uZTJsHo/cFt8D43kSUiIsgoHHHjWdoedwMTvz8Fwz3BKHxAayx8IkCcpoiISHKYDEgWfr+YidCNJ6G/JxnN7OeLiCfbQqFQiNQZERFJDcMR1XmHL2dj1PoTKNYbBPXXe7XCxyPbMRgREZEAwxHVaaeu5WL42mMoLNYL6lO7t8SKkA4MRkREZMKi4chgMMDPzw8TJ06s1sfPmzcPzs7O9/01efJkS7ZMddj5jDwMWfMXbheVCuoTunhh9dhODEZERGSWRcPR3r17kZGRUe2PT0xMtGA3JGeXbxZi8Fd/IaugWFAf0c4DG54Jgp2SwYiIiMyz2NtqSUlJCA8Pr9E5Ll26BOBuyOrVq5cl2iIZSsstwqCvjuJ6bpGgPqC1G3567lGo7Pg0mYiIKlejcHT69Gls2rQJsbGxiImJgcFgePAHVcJgMODy5csAgMDAwJq0RTKWla/D4K+OIvlmoaDevWUj/PJCN9RT2YnUGRER2YoahaPDhw9j1apVFmkkNTUVRUVF8PT0hKurq0XOSfJyu6gEQyOPIT4jX1Dv1NQFe8K6o4GjbJb1IiKiGqjR84UJEybg2LFj5b+mTp1a7XOVPVILCAhAZGQk+vTpgyZNmsDLywuDBw/Gd999V6M7U1S3FRaXYsS644i9liuot2lSH/te7gFXJ66ETkREVVOjf0q7ubnBzc2t/Pc12aesbDL2oUOHcPDgwfJ6UVERjh49iqNHj2LXrl3YuHEjVCpV9ZumOqdEb8D4DbE4lJwtqLdo5Ij9r/SARwO1SJ0REZEtksxzhrI7RwaDAc8//zzeeOMNeHt74/r161i7di1WrFiBXbt2YfHixVi4cOF9zxUcHFylz5mcnAwfHx9otdoa90+mMjMza/1zGI1GvLEnGbvPZwnq7vVV+GlcABxL8qHV5lfy0fJljWtD1cNrI128NtJVWlpq0RsnknltR6/XIzAwEOHh4fjiiy/Qpk0bODg4wMfHB0uWLMGcOXMAACtXrkRWVtYDzkZysejPVPx4Vvj3oaHaDj+ObwtfV0eRuiIiIlsmmTtHn3766X3HZ8yYgS+++AI5OTn4888/MXbs2EqPjYmJqdLnDA4OhlKpNNmQliyrtr6+//4zCSuP3xDU6qmU2B3WA718OKm/Kvh3X7p4baSL10Z67O0tG2ckc+foQRwdHcsfl5U9giP52hCTird3xgtqdkoFop4LZjAiIqIasZlwBKB88ndRUdEDjqS6bHd8Bqb8eMakvm58ZzzZzkOEjoiIqC6RxGO1zMxMnD59GiqVCv3796/0uNu3bwOA4A05kpcjl7MxbkMM9AajoP7RiHaY3LWFSF0REVFdIolwlJeXh5CQEABAbGwsAgICTI4pLS1FbGwsACAoKMiq/ZE0nEvPw4h1x3GnRLje1dv9/fD2AD+RuiIiorpGEo/VfH190alTJwDA0qVLzR6zdu1aaLVatGzZEj179rRmeyQBKbcK8cSav3DrTomg/lxwcyx7ktvNEBGR5Vg9HAUFBSEoKAjz588X1GfNmgUA2Lp1K8LCwnD+/HnodDqkpKRgyZIl5ZvaLlq0CHZ23B9LTnLulGD42uMmG8k+GajB2vGdoVQqROqMiIjqIqs/VitbCTs9PV1QDwkJwWuvvYaVK1di8+bN2Lx5s8nHzp8/H2PGjLFKnyQNulI9Qr45gXPpeYL6Y60a46fnHoXKThI3P4mIqA6R1E+WZcuWISoqCkOGDIGbmxvs7e2h0WgwatQo7N+/v/zuEsmDwWDElB/O4M+km4J6oIczdr7YDU4OkpgyR0REdYxFf7rMmTOnfCXryuTn338rh2HDhmHYsGGWbIts1JxfL+D7U9cFNc8Gavw6tTs3kiUiolojqTtHRGVWH7mCfx0QLvbprLbDnqnd4e3qJFJXREQkB7J5LqHVak32ZNPpdHBw4B0Iqdl5Lh2v/RwnqNkpFdjyXDCCmjcUqSsiIpIL2YSjyMhIREREmNTd3d1F6IYqcyIlBxO+O4l71njEmrGd8ERb7mdERES1TzbhKCwsrHyhyTKhoaG8cyQhyTcLMGLdMRQW6wX1eYP9MaV7S5G6IiIiuZFNONJoNCY7KavVaiiVnHYlBTcLijEs8hi0+cWC+vNdW2DBE/4idUVERHLEZECi05XqEbL+BBIyCwT1wf5NsGZcJygUXOSRiIish+GIRGU0GvFy1N84lJwtqHf2csGWycFc5JGIiKyOP3lIVP86cAnfxlwT1Jo1dMTuqd3g4qgSqSsiIpIzhiMSzda/0zB7zwVBrb6DHXa92A3NGtYTqSsiIpI7hiMSRUxqDiZ9f0pQUyiA7yc+gi7NuJYRERGJh+GIrC711h2MXHccd0oMgvpHI9rhHx08ReqKiIjoLoYjsqp8XSlGfn0c6Xk6QX1q95aY2c9XpK6IiIj+h+GIrEZvMOKZ707iTNptQX1g6yb4ckxHvrJPRESSwHBEVvPOrnjsjM8Q1Pzd62PL5Ef5yj4REUmGbFbI5saz4lr711V8cjBZUHN1UmHXi93Q2InXgIiIpEM24Ygbz4rn8OVsTPs5TlBT2Snw8/PBaOPuLFJXRERE5skmHHHjWXFcv63DmO9Oo0RvFNS/GtsJ/fyaiNQVERFR5WQTjrjxrPXdKTHghW2JyLjnzbS3+vnihW4tReqKiIjo/pgMqFYYjUbM/C0ZZ9KFm8kO8XfHshHtROqKiIjowRiOqFZ8/GcSfo6/Kai1aVIfP0x6BHZKvrJPRETSxXBEFvfbBS3Cd58X1Bqo7fHLlK58M42IiCSP4Ygs6qI2HxM2xsJYYf61QgFsmhiEQI8G4jVGRERURQxHZDG5d0rw1NfHkVtUKqgvHtoWI9tzzzQiIrINDEdkEQaDEZO+P4WLmcIJ2E+1dcV7j7cWqSsiIqKHx3BEFhFxINFka5AOGid8OsyXe6YREZFNYTiiGtt3UYv3f7soqDWp74D1o/1R38FOpK6IiIiqRzaLQFLtuJJdiKe/OymYgK1UAD88+whaNDRW/oFEREQSxTtHVG1FJXqM+TYG2YUlgvrS4YF43J971hERkW2SzZ0jrVaLrKwsQU2n03FvtWoyGo147ec4nLyWK6iHdPTEOwP8ROqKiIio5mQTjiIjIxEREWFSd3fnHY7qWHssBV8fTxXU/N3rY/2ELpyATURENk024SgsLAwhISGCWmhoKO8cVcPxlFt4/eezglp9Bztse74rXBxVInVFRERkGbIJRxqNBhqNRlBTq9VQKjnt6mFk5usw9tsYFOsNgvrXoV3QzpMrYBMRke1jMqAqMxiMeHbTKaTmFAnqM/v5YnwXL5G6IiIisiyGI6qyiAOJ2JeQKaj183PDsicDReqIiIjI8hiOqEoOJmVh3j0LPXo2UOOHZx+BvR3/GhERUd3Bn2r0QNo8HZ7+7iQM9y70OOkReLo4itcYERFRLWA4ovsyGIx49vuTuHFbJ6h/MDQA/fyaiNQVERFR7WE4ovta+n+J+D1BuHjmYP8meG9gG5E6IiIiql0MR1SpPy9lYf5e4Tyjpi5qfPfMI1AqudAjERHVTQxHZFZGJfOMNj/7CDQN1OI1RkREVMsYjsiE3mDEs5tOIj2P84yIiEh+ZLNCNjeerbp/HUjE/kTh12qIvzvnGRERkSzIJhxx49mqOXolG/P3JghqXi6O2PhMEOcZERGRLMgmHHHj2QfLvVOCZzadhL7CRCPOMyIiIrmRTTjixrP3ZzQa8erWOFzJviOozxvsj75+biJ1RUREZH1MBgQA2Bh7DZtPXRfUerVqjDmDOM+IiIjkheGIcCmrAK/9HCeoNXS0x6aJ3DeNiIjkhz/5ZK641IBnvjuJfJ1eUF8zrjO8XZ1E6oqIiEg8kg9HBoMBfn5+mDhxotit1Enz917EidQcQW1KtxYY38VLnIaIiIhEJvlwtHfvXmRkZIjdRp10IDELy/64JKj5u9fH8lEdROqIiIhIfJIOR0lJSQgPDxe7jTop504JJm8+BWOF7UFUdgp8P/EROKtl8xIjERGRCcn9FDx9+jQ2bdqE2NhYxMTEwGAwiN1SnfT6z3G4llskqC0dFohHWzQSpyEiIiKJkFw4Onz4MFatWiV2G3Va1Jk0bDopfG1/YOsmmNnPV6SOiIiIpENyj9UmTJiAY8eOlf+aOnWq2C3VKTduF+GVLX8Lag0d7bF+QhduD0JERAQJ3jlyc3ODm9v/VmTm3meWYzQa8eKPZ5BdWCKofzG6I1o0ridSV0RERNIiuXBkCcHBwVU6Ljk5GT4+PtBqtbXckTR8eyoDv14Q/llHBrhicHNVrXwNMjMzLX5OsgxeG+nitZEuXhvpKi0thUqlstj5JPdYjWpHcnYR5v+RIqhp6qvw4ROtoFDwcRoREVGZOnnnKCYmpkrHBQcHQ6lUmmxIW9eU6g146ocjuFMifPNv/dNBaOvtUeufv65/fW0Zr4108dpIF6+N9NjbWzbO8M6RDHz4RxL+unpLUHu5pzeGBdZ+MCIiIrI1DEd13Ln0PCzclyCo+bk54eOR7UTqiIiISNoYjuqwUr0BL/xwGsX6/z1OUyqAjc8EcRVsIiKiSjAc1WH/PphssqnsW/380LOVqzgNERER2QCGozrqQkYe5u+9KKgFuNfHwqEBInVERERkGxiO6iC9wYgXfjwDXen/HqcpFMDXoV1QT2UnYmdERETSx3BUB30WnWzydtqbfX3xmA8fpxERET0Iw1Edk5CZj7m/XhDUWjepj0V8nEZERFQlkg9Hc+bMQX5+PjZt2iR2K5KnNxgx5YfTKDJ5nNYZTg58O42IiKgqJB+OqOpWHr6Mw1eEj9Om9/ZBH1+3Sj6CiIiI7iWb2wlarRZZWVmCmk6ng4ODg0gdWVbKrULM3iN8nObr5oSlw9qK1BEREZFtkk04ioyMREREhEnd3d1dhG4sy2g04rWfz6KgWC+orxvfGfW52CMREdFDkc1PzrCwMISEhAhqoaGhdeLO0da/b2BXfIagFtajJfq3biJSR0RERLZLNuFIo9GY7KSsVquhVNr2tKucOyWYvu2soObRQI1lTwaK1BEREZFts+1kQHhv93mk5+kEteVPtUdjJ9u/I0ZERCQGhiMbdvhyNlYfvSqoDQ/UYHwXL5E6IiIisn0MRzaquNSAl6LOCGpODnZYObojFAqFSF0RERHZPoYjG/XhH5cQn5EvqC0aGoBWrk4idURERFQ3MBzZoITMfCzenyioPdK8If7Z20ekjoiIiOoOhiMbYzQaMW1rHHQVtghRKoDIcZ1gb8fLSUREVFP8aWpjos7cwP8lClf6ntHXF480byROQ0RERHUMw5ENySsqxZu/nBPUWjRyxMInAkTqiIiIqO5hOLIhH/yegLTbRYLap0+1hzO3CCEiIrIYhiMbcS49D59FJwtqTwS4Y3THpiJ1REREVDfJ5paDVqtFVpZwro5Op7OJvdWMRiNe/zkOpQZjec3BTonPQzpwTSMiIiILk004ioyMREREhEnd3d1dhG4ezg+n0vBn0k1BbdYAP7RxdxapIyIiorpLNuEoLCwMISEhglpoaKjk7xzdLirBWzuFk7BbNq6H2Y+3FqkjIiKiuk024Uij0UCj0QhqarUaSqW0p10t3JeAG7dNN5Z1cpDNpSMiIrIqaScDmbuozceKQ5cFteGBGjzVwVOkjoiIiOo+hiMJe2vHOZNJ2MtHcRI2ERFRbWI4kqjfLmix+7xWUJvZzxetm9QXqSMiIiJ5YDiSoBK9wWQlbM8Gasx+vI1IHREREckHw5EErTpyBRe0+YLa0uFt0cCRk7CJiIhqG8ORxGTl6zB/b4Kg9mjzhpgc3EKkjoiIiOSF4Uhi5u9NQM6dEkFt+agOUCo5CZuIiMgaGI4kJO7Gbaw+ekVQm9DFC718XMVpiIiISIYYjiTk7R3xqPDmPuqplFg2IlC8hoiIiGRINjN8pb7x7L6LWuxLyBTU3hnQGi0bO4nUERERkTzJJhxJeeNZvcGId3adF9S8XBwxq7+fSB0RERHJl2zCkZQ3nt108hrOpN0W1BYNDUB9tWwuDxERkWTI5qevVDeevVOix9xfLwhq7T0bYHJXvrpPREQkBk7IFtmKQ5eRmlMkqH04IhB2fHWfiIhIFAxHIsrK12Hp/yUKagNbN8GwtppKPoKIiIhqG8ORiBbvT8TtolJB7cMRgVAoeNeIiIhILAxHIknKKsCXR64Ias8ENcOjLRqJ0g8RERHdxXAkkrm/XkCJ/n8rPjrYKbFkeFsROyIiIiKA4UgUZ9Jy8cPpNEFteu9WaOXKBR+JiIjExnAkgvd/vSj4fUNHe8we1EakboiIiKgihiMr++vqLeyMzxDU3u7vB1cn8RejJCIiIoYjq5uzR7jgY5P6Dnijj69I3RAREdG9GI6s6EBiFg5cEm5+O/vx1mjgKJuFyomIiCRPNj+VtVotsrKEwUSn01ltbzWj0Yg592wT0qyhI159rJVVPj8RERFVjWzCUWRkJCIiIkzq7u7uVvn8u89r8dfVW4La+4PbwFFlZ5XPT0RERFUjm3AUFhaGkJAQQS00NNQqd44MBqPJXCNfNydM6day1j83ERERPRzZhCONRgONRrhnmVqthlJZ+9Ouos6k4e8btwW1hU8EQGXHKV9ERERSw5/OtUxvMGLBvgRBrZ2HM54OaiZSR0RERHQ/NQ5HRqMRkZGR6NOnDzw9PdGiRQsMHz4ce/bssUR/Nm/LmTRc0OYLaouGtoWdkpvLEhERSVGNwpHRaMTEiRPx5ptv4tSpU8jPz8etW7cQHR2N8ePH41//+tdDnW/evHlwdna+76/JkyfXpGWrMhiMWLQ/UVDr7OWCkI6eInVERERED1KjcLRy5Urs2LEDarUay5cvR1paGpKSkjBjxgwAwJIlS3Do0KEqny8xMfHBB9mQbWdv4Fx6nqA2b7A/FAreNSIiIpKqak/I1ul0+OSTTwDcDUEvvvgiAMDFxQWLFy9GVlYWvvvuO3z00Ufo06dPlc556dIlAMDevXvRq1ev6rYmCQaDER/sE4a9Dp4NMKoD7xoRERFJWbXvHB0+fBharRaurq6YMmWKyfjMmTMBAAcPHkROTs4Dz2cwGHD58mUAQGBgYHXbkoyd8Rkmb6jNG+IPJecaERERSVq1w1F0dDQAoG/fvmbXCvL394e3tzf0ej2OHDnywPOlpqaiqKgInp6ecHV1rW5bkmA0GvHB76ZvqI3p2FSkjoiIiKiqqh2OyuYHdezYsdJjysaqMpeo7JFaQEBA+dtvTZo0gZeXFwYPHozvvvsOBoOhuu1a1Z7zWpy8liuozR3Eu0ZERES2oNpzjq5duwYAaNas8vV6vLy8AAApKSkPPF9ZgDp06BAOHjxYXi8qKsLRo0dx9OhR7Nq1Cxs3boRKpapu27XO3F0jf/f6GN/FS6SOiIiI6GFUOxzl599du8fZ2bnSY8rGCgoKHni+sjtHBoMBzz//PN544w14e3vj+vXrWLt2LVasWIFdu3Zh8eLFWLhw4X3PFRwcXKU/Q3JyMnx8fKDVaqt0fFUcSM7B8ZQcQe2f3TxwMyvTYp/DVmRmyu/PbCt4baSL10a6eG2kq7S01KI3Tqr9WK24uBgA7rs3WVmjhYWFDzyfXq9HYGAgwsPD8cUXX6BNmzZwcHCAj48PlixZgjlz5gC4u3xAVlZWdduudZ8euS74fatGaoS0ayJSN0RERPSwqn3nqCwU6XS6So8pG1Or1Q8836effnrf8RkzZuCLL75ATk4O/vzzT4wdO7bSY2NiYh74+YC7d5iUSqXJnmvV9Z/kmzh+Xbga9rwn2sLL08Mi57dVlvr6kuXx2kgXr4108dpIj729ZbeKrfado7JHZmWP18zJy7u7AGL9+vWr+2nKOTo6lj8uK3sEJzXL/kgS/L5FI0c8+2hzkbohIiKi6qh2OGre/O4P/evXr1d6zI0bNwTH1pSbmxuAu5O0pebsjdvYFZ8hqL3Vzw8qO+7tS0REZEuqfR/K398fABAXF1fpMWfPnhUcW5nMzEycPn0aKpUK/fv3r/S427fvLqpYFpKk5KM/hXeNXJ1UmNq9pUjdEBERUXVV+7ZG3759AdxdDLJscnZFCQkJSElJgZ2dHXr37n3fc+Xl5SEkJAQjRozAxYsXzR5TWlqK2NhYAEBQUFB1264VKbcK8f1J4R2013v5oL7ass9AiYiIqPZVOxz16tULHh4euHXrFtavX28yvmLFCgDAgAEDHrjita+vLzp16gQAWLp0qdlj1q5dC61Wi5YtW6Jnz57VbbtWfHIwGaUGY/nv66mUmN67lXgNERERUbVVOxw5ODiU7582e/ZsbNiwAXl5ecjIyMD8+fOxfv16KJVKvPvuu4KPCwoKQlBQEObPny+oz5o1CwCwdetWhIWF4fz589DpdEhJScGSJUsQHh4OAFi0aBHs7Oyq27bF3SwoRuQx4SKXL3ZriSbOD35Dj4iIiKSnRs99pk2bhqNHj2L79u2YNm0apk2bJhj/4IMP0KNHD0GtbCXs9PR0QT0kJASvvfYaVq5cic2bN2Pz5s0mn2/+/PkYM2ZMTVq2uJWHr6CwWF/+ezulAm/19xOxIyIiIqqJGr1KpVAosHHjRnz22WcICgpC/fr10ahRI/Tr1w9bt27FjBkzHup8y5YtQ1RUFIYMGQI3NzfY29tDo9Fg1KhR2L9/f/ndJakoLC7F5/+5LKiFdvZCK1cnkToiIiKimqrxjGGFQoGpU6di6tSpVTr+fusiAcCwYcMwbNiwmrZlFV8fT0VWgXAyevjA1iJ1Q0RERJbARXiqSW8w4pODyYLasLYadPJyEakjIiIisgTZvGuu1WpN9mTT6XT33RvufnacS8flbOGeceEDOdeIiIjI1skmHEVGRiIiIsKk7u7uXq3zfRotvGv0aPOG6OsrvcUpiYiI6OHIJhyFhYUhJCREUAsNDa3WnaPY1BwcSs4W1Gb09YVCoahRj0RERCQ+2YQjjUZjspOyWq2GUvnw064+OyS8a9TURY3xnb1q1B8RERFJAydkP6S03CL8cCpNUHu9lw8c7PmlJCIiqgv4E/0hfXnkimCrEEd7JV7qwQ1miYiI6gqGo4dwp0SP1UeuCGrPBTfnViFERER1CMPRQ/gu9hpuFpYIam/08RWpGyIiIqoNDEdVZDQa8dk9r+8/EeCOdp4NROqIiIiIagPDURXtu5iJ+Azh1idv9uVdIyIiorqG4aiKlh8SbjAb6OGMIQHVW0CSiIiIpIvhqAouZRXg1wtaQW1GHy76SEREVBcxHFXBqnveUGtcT4VnH20mTjNERERUqxiOHqCwuBRfH08V1KZ0awEnB9ksLk5ERCQrsvkJr9VqkZWVJajpdLoH7q32w6k05NwRvr7/ymOtLN0eERERSYRswlFkZCQiIiJM6u7ulU+qNhqNWHnPI7Whbd3Rukl9S7dHREREEiGbcBQWFoaQkBBBLTQ09L53jo6l5ODktVxB7bVePrXSHxEREUmDbMKRRqOBRqMR1NRqNZTKyqddrTwsfH2/lWs9DGurqeRoIiIiqgs4IbsSmfk6/HT6hqD2Ss9WsFPy9X0iIqK6jOGoEuuOpaBYbyj/vdpeiRe7tRCxIyIiIrIGhiMz9AYjVh+9KqiFdvFCE2e1SB0RERGRtTAcmbE7PgNXb90R1Kbx9X0iIiJZYDgy48t7Xt9/tHlDdGvZSJReiIiIyLoYju6RmJmPvRczBbXXerXiPmpEREQywXB0j1VHhHONGtdTYUIQ91EjIiKSC4ajCgqLS/HNCdN91Oqp7ETqiIiIiKyN4aiCzffso6ZQAK9yIjYREZGsyGaF7AdtPGs0Gk1WxB4aoIEf91EjIiKSFdmEowdtPPvX1Vs4df22YGxar1bWaI2IiIgkRDbh6EEbz648fEUwxn3UiIiI5Ek24eh+G89q83SIOiPcR+1V7qNGREQkS5yQDWCtmX3UpnAfNSIiIlliOAKw6p4VsSdwHzUiIiLZkn04ul1Uimu5RYLa6719ROqGiIiIxCb7cHSrwrpGANC9ZSMEt2gkTjNEREQkOtmHo4LiUsHvX+Pr+0RERLIm+3BUkbuzA8Z19hK7DSIiIhIRw1EFU7u3hCP3USMiIpI1hqP/UiqAV3p6i90GERERiYzh6L+e6uCJlo2dxG6DiIiIRMZw9F+vPdZK7BaIiIhIAmSzfYhWq0VWVpagptPpAABtNc4Y2KaJGG0RERGRxMgmHEVGRiIiIsKk7tSkKWY/3hoKBfdRIyIiIhmFo7CwMISEhAhqoaGhcHBwwKRg7qNGREREd8kmHGk0Gmg0GkFNrVZDqeS0KyIiIvofJgMiIiKiCiQZjoxGIyIjI9GnTx94enqiRYsWGD58OPbs2SN2a0RERFTHSe6xmtFoxMSJE7Fjxw5BPTo6GtHR0Zg7dy7effddkbojIiKiuk5yd45WrlyJHTt2QK1WY/ny5UhLS0NSUhJmzJgBAFiyZAkOHTokbpNERERUZ0kqHOl0OnzyyScA7oagF198ES4uLvDw8MDixYvx7LPPwmg04qOPPhK5UyIiIqqrJBWODh8+DK1WC1dXV0yZMsVkfObMmQCAgwcPIicnx8rdERERkRxIKhxFR0cDAPr27QsHBweTcX9/f3h7e0Ov1+PIkSPWbo+IiIhkQFLhKDExEQDQsWPHSo8pGys7loiIiMiSJBWOrl27BgBo1qxZpcd4eXkBAFJSUqzSExEREcmLpF7lz8/PBwA4OztXekzZWEFBQaXHBAcHV+nzJScnw8fHB1qt9iG6pKrKzMwUuwWqBK+NdPHaSBevjXSVlpZCpVJZ7HySunNUXFwMAGbnG5Up+8MXFhZapSciIiKSF0ndOSoLRTqdrtJjysbUanWlx8TExFTp8wUHB0OpVJrsuUaWxa+vdPHaSBevjXTx2kiPvb1l44yk7hyVPTIre7xmTl5eHgCgfv36VumJiIiI5EVS4ah58+YAgOvXr1d6zI0bNwTHEhEREVmSpMKRv78/ACAuLq7SY86ePSs4loiIiMiSJBWO+vbtC+DuYpBlk7MrSkhIQEpKCuzs7NC7d29rt0dEREQyIKlw1KtXL3h4eODWrVtYv369yfiKFSsAAAMGDICrq6uVuyMiIiI5kNzbajNnzkR4eDhmz54NR0dHhISEoLCwEF9++SXWr18PpVKJd9991yKf79q1aygpKUH79u0tcj4S0uv1AAA7OzuRO6F78dpIF6+NdPHaSNfly5fvuwzQw1Lk5+cbLXY2CzAajZg0aRK2b99udnzx4sWYMWOGRT6Xr68vCgsLBZO79Xo9bt26hcaNG1vsfwBbOGdt9JicnAzg7tfZEmzh61gb57SFawPYxp/bFnrktZFuj7w20j1nQkICFAoFcnNzLdCdBMMRcDcgrVu3Dt9++y0SEhKgUqnQuXNn/POf/8QTTzxRq587Pj4e3bp1w/Hjx9GuXTvZnLM2eixbqbyq6049iC18HWvjnLZwbQDb+HPbQo+8NtLtkddGuue09LWR1GO1MgqFAlOnTsXUqVPFboWIiIhkRlITsomIiIjExnBEREREVAHDEREREVEFDEf3aNKkCd577z00adJEVuesjR4tzRa+jrVxTlu4NoBt/LltocfaYAt/blvosTbYwp/bVs5pSZJ8W43qhtp4s4Msg9dGunhtpIvXRrosfW1454iIiIioAoYjIiIiogoYjoiIiIgq4JwjIiIiogp454iIiIioAoYjIiIiogoYjoiIiIgqYDgiIiIiqoDhiIiIiKgChiOyKKPRiMjISPTp0weenp5o0aIFhg8fjj179ojdGhFRtcTGxmLy5Mnw9/dH48aN0axZMwwePBjr1q2DXq8Xuz2qBXyVnyzGaDRi4sSJ2LFjh9nxuXPn4t1337VyV1QZg8GANm3aoEePHti0aZPY7cheRkYGPv/8c/z2229ISUkBALRq1QojR47E9OnT0ahRI3EblKmoqCiEhYWhtLTU7PigQYMQFRUFlUpl5c6oMu+//z4+/fRTfPzxx3jllVeqdQ7eOSKLWblyJXbs2AG1Wo3ly5cjLS0NSUlJmDFjBgBgyZIlOHTokLhNUrm9e/ciIyND7DYIQHx8PB577DF89tlnuHDhAgoLC1FYWIj4+HgsW7YMvXr1wqVLl8RuU3aysrIwffp0lJaWomvXrti3bx/S0tJw/vx5LFy4EA4ODti/fz8+/PBDsVul/zpy5AiWL19e4/MwHJFF6HQ6fPLJJwDuhqAXX3wRLi4u8PDwwOLFi/Hss8/CaDTio48+ErlTAoCkpCSEh4eL3Qbh7h3XF154ARkZGfDz88PPP/8MrVaLpKQkREZGQqPR4OrVqxg/fjxKSkrEbldWtmzZgvz8fHh5eWH37t147LHH4OLighYtWuCtt97CrFmzAADr168Xt1ECAOTl5eGll16CwWCo8bkYjsgiDh8+DK1WC1dXV0yZMsVkfObMmQCAgwcPIicnx8rdEQCcPn0as2bNwsCBAxEUFITk5GSxWyIABw4cwLlz56BSqbBt2zYMGTIETk5O8PDwwNNPP439+/fDyckJCQkJ2L59u9jtykpsbCwAYOTIkXBycjIZHz16NADgxo0byM7OtmpvZCo8PBxXrlyxyLkYjsgioqOjAQB9+/aFg4ODybi/vz+8vb2h1+tx5MgRa7dHuBtgV61ahePHj1vkX1ZkGX/88QcAoH///vD19TUZ9/X1RUhICIC715CsR6vVArg798ucBg0alP+30cjpu2LavXs3NmzYgA4dOqB79+41Ph/DEVlEYmIiAKBjx46VHlM2VnYsWdeECRNw7Nix8l9Tp04VuyXC3UecABAYGFjpMRqNBgBQUFBglZ7orl9++QX5+fmYPn262fGysOrp6Qk3NzdrtkYVZGZm4vXXX4eDgwPWrl1rkcnx9hboiwjXrl0DADRr1qzSY7y8vACg/E0csi43NzfBN3B3d3cRu6Eyr776KkaPHo127dpVeszp06cBAN7e3lbqiipTVFSE9PR0REdHY+7cuQBQPveIxDF9+nRkZmbigw8+QIcOHSxyToYjsoj8/HwAgLOzc6XHlI3xX79E/9O3b9/7ju/fv7/80duIESOs0RJVYvXq1Xj77bfLf+/k5IRPPvkEL730kohdyduGDRuwa9cu9OzZs/zNaEvgYzWyiOLiYgAwO9+oTNmtzsLCQqv0RGTrNm3ahIkTJwIAxowZgy5duojbEAkUFhbi999/x40bN8RuRZauXr2K8PBw1K9fH1999RWUSstFGoYjsoiyUKTT6So9pmxMrVZbpSciWxUXF4fhw4fj5ZdfRkFBAXr37o3Vq1eL3ZbsvfLKK8jPz8elS5ewadMmtGnTBr/++iuGDx9+3+99ZHkGgwFhYWHIy8vD0qVLzb7MUBMMR2QRZY/Myh6vmZOXlwcAqF+/vlV6IrI1ubm5ePPNN9GrVy9ER0dDpVJhzpw52LVrF+rVqyd2e/Rfnp6eeOqpp7Br1y40bNgQiYmJ2LZtm9htycry5ctx5MgRDB48GC+++KLFz89wRBbRvHlzAMD169crPabs1nPZsUT0PzExMejRowciIyNhMBgwcuRIxMTE4L333oO9PaeHSlGzZs3K54zFxcWJ3I18XL58GYsWLYKrqytWrVpVK5+D4Ygswt/fH8D9v0GcPXtWcCwR3XXw4EEMHz4cqamp8Pb2xp49e7B582b4+fmJ3ZpspaWlwcPDAx4eHuVv45pT9tZn2Z1xqn2pqakoLi5GdnY2WrduDWdnZ8Gv//znPwCAt99+u7z2sIsPMxyRRZT96yk6Orp8cnZFCQkJSElJgZ2dHXr37m3t9ogk6+bNm5g0aRIKCwvRr18/HD169IFvsFHt8/DwgNFoREFBwX33tStbt+1+y5iQ7WE4Iovo1asXPDw8cOvWLbP7DK1YsQIAMGDAALi6ulq5OyLpWrNmDbKzs+Ht7Y2oqCi4uLiI3RIBsLOzKw+p3377rdljTpw4Ub4Q5OOPP2613uSub9++yM/Pr/RX2T/AP/744/Jao0aNHupzMByRRTg4OJTvnzZ79mxs2LABeXl5yMjIwPz587F+/XoolUq8++67IndKJC27d+8GAEyZMsXs/l0knmnTpgEAoqKi8MILL+DMmTMoKChAamoqvvnmG4wZMwYGgwFDhw5FcHCwyN2SJTEckcVMmzYNo0aNQlFREaZNm4amTZvCz88P//73vwEAH3zwAXr06CFyl0TSUVJSgnPnzgEA5s+fbzJ34t5f77zzjsgdy8vAgQPLF32Miooqv0MeGBiI6dOnIzs7G926dUNkZKTInZKlMRyRxSgUCmzcuBGfffYZgoKCUL9+fTRq1Aj9+vXD1q1bLbp6KVFdoNVqUVJSInYbdB8LFizA9u3bMXz4cGg0Gtjb26NRo0Z47LHHsHz5cuzduxeNGzcWu02yMEV+fj63EiYiIiL6L945IiIiIqqA4YiIiIioAoYjIiIiogoYjoiIiIgqYDgiIiIiqoDhiIiIiKgChiMiIiKiChiOiIiIiCpgOCIiIiKqgOGIiIiIqAKGIyIiIqIKGI6IiIiIKmA4IiIiIqqA4YiIiIioAoYjIiIiogoYjoiIiIgq+H8WgdBK8Q1EpgAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "Figure(nrows=1, ncols=1, refwidth=2.5)"
+ ]
+ },
+ "metadata": {
+ "image/png": {
+ "height": 282,
+ "width": 291
+ }
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig, axs = uplot.subplots()\n",
+ "\n",
+ "axs.plot(nu_vals, frac_vals)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "e786f3d4-557a-4f76-8836-fc8c61763b7c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import scipy.optimize\n",
+ "import numpy.polynomial\n",
+ "from numpy.polynomial import Polynomial\n",
+ "\n",
+ "numpy.polynomial.set_default_printstyle(\"ascii\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "00c73375-0648-46e5-bae6-4048a0a867fd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RATNL_ORDER = 9\n",
+ "\n",
+ "\n",
+ "def get_ratnl_func_polys(coeff):\n",
+ " p_coeff = coeff[0 : RATNL_ORDER + 1]\n",
+ " q_coeff = np.concatenate([[1], coeff[RATNL_ORDER + 1 :]])\n",
+ " pm = Polynomial(p_coeff)\n",
+ " qm = Polynomial(q_coeff)\n",
+ " return pm, qm\n",
+ "\n",
+ "\n",
+ "def ratnl_func(x, *coeff):\n",
+ " pm, qm = get_ratnl_func_polys(coeff)\n",
+ " return pm(x) / qm(x)\n",
+ "\n",
+ "\n",
+ "res = scipy.optimize.curve_fit(\n",
+ " ratnl_func, nu_vals, frac_vals, p0=np.ones(2 * RATNL_ORDER + 1), full_output=True\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "b2d7adf1-ba49-4ac3-9924-39a7938e0b71",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "coeff = res[0]\n",
+ "\n",
+ "pm, qm = get_ratnl_func_polys(coeff)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "0d14d88c-cb83-4c21-a11b-0021614c0a1f",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "def _spergel_hlr(x):\n",
+ " pm = 1.2571513771129166 + x * (\n",
+ " 3.7059053890269102 + x * (\n",
+ " 2.8577090425861944 + x * (\n",
+ " -0.30570486567039273 + x * (\n",
+ " 0.6589831675940833 + x * (\n",
+ " 3.375577680133867 + x * (\n",
+ " 2.8143565844741403 + x * (\n",
+ " 0.9292378858457211 + x * (\n",
+ " 0.12096941981286179 + x * (\n",
+ " 0.004206502758293099\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " qm = 1.0 + x * (\n",
+ " 2.1939178810491837 + x * (\n",
+ " 0.8281034080784796 + x * (\n",
+ " -0.5163329765186994 + x * (\n",
+ " 0.9164871490929886 + x * (\n",
+ " 1.8988551389326231 + x * (\n",
+ " 1.042688817291684 + x * (\n",
+ " 0.22580140592548198 + x * (\n",
+ " 0.01681923980317362 + x * (\n",
+ " 0.00018168506955933716\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " return pm / qm\n"
+ ]
+ }
+ ],
+ "source": [
+ "def make_poly_code(pm, head=\"\", base_indent=0):\n",
+ " res = \"\"\n",
+ " indent = base_indent\n",
+ " for c in pm.coef:\n",
+ " if c == pm.coef[-1]:\n",
+ " end = \"\"\n",
+ " else:\n",
+ " end = \" + x * (\"\n",
+ "\n",
+ " if c == pm.coef[0]:\n",
+ " _hd = head\n",
+ " else:\n",
+ " _hd = \"\"\n",
+ " res += \" \" * 4 * indent + f\"{_hd}{c}{end}\\n\"\n",
+ " if c != pm.coef[-1]:\n",
+ " indent += 1\n",
+ "\n",
+ " for _ in pm.coef[:-1]:\n",
+ " indent -= 1\n",
+ " res += \" \" * 4 * indent + \")\\n\"\n",
+ "\n",
+ " return res\n",
+ "\n",
+ "\n",
+ "code = \"def _spergel_hlr(x):\\n\"\n",
+ "code += make_poly_code(pm, head=\"pm = \", base_indent=1)\n",
+ "code += make_poly_code(qm, head=\"qm = \", base_indent=1)\n",
+ "code += \" return pm / qm\"\n",
+ "\n",
+ "print(code)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "8d5d7152-ea9e-4569-b5b7-9a10037ecdfa",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "9.072944411797624e-09 6.270160715637907e-08\n",
+ "9.072944411797624e-09 6.270160715637907e-08\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAI1CAYAAADcqV49AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAewgAAHsIBbtB1PgAAX1dJREFUeJzt3XlYVGXDBvCbdRAQF2QgNwQUBMWi0Mw91zRNEZXKzFyw0sxdTC01RTJNc3s1R8tSX33DPdxSyyU1C94sBBcUFVFhQERZZGAYvj984eM4gyIMc85w7t91eV0fz3M83HneT27POfM8FtnZ2UUgIiIiIgCApdgBiIiIiKSE5YiIiIioFJYjIiIiolJYjoiIiIhKYTkiIiIiKoXliIiIiKgUliMiIiKiUliOiIiIiEphOSIiIiIqheWIiIiIqBRrsQOIydPTE7m5uWjYsKHYUYiIiKiCkpOTYW9vj8TERKOcT9blKDc3FwUFBbC05A20qlBYWAgAsLKyEjkJPY7XRrp4baSL10a68vPzjXo+WZejhg0bwtLSEnFxcWJHqZbUajUAQKlUipyEHsdrI128NtLFayNdzZs3h4WFhdHOx1smRERERKWwHBERERGVwnJEREREVArLEREREVEpLEdEREREpcjm02pqtRrp6emCMY1GA1tbW5ESERERkRTJphypVCpERETojbu4uIiQhoiIiKRKNuUoNDQUQUFBgrGQkBDeOSIiIiIB2ZQjpVKpt3CXQqHg6thEREQkwGZAREREVArLEREREVEpRilHMTExGD58OLy9vVGnTh00aNAAPXr0wIYNG0o26iMiIiIyB5UuR5GRkejWrRt27NiB27dvo6CgAPfv38eZM2cwYcIEBAcHo6CgoNzn++yzz+Do6PjEX8OHD69sbCIiIiKDKlWO0tPTMX78eGi1WrRu3Ro///wzbt++jQsXLmDevHmwtbXFkSNH8OWXX5b7nAkJCZWJRERERFQplfq02vbt25GdnY369etj3759sLe3BwA4OTlhypQpyM/PR3h4ODZu3IhZs2aV65xXrlwBABw6dAjt27evTDwiIiKiZ1apO0cxMTEAgH79+pUUo9IGDhwIALhz5w4yMjKeej6dTodr164BAHx9fSsTjYiIiKhCKlWO1Go1AKBJkyYG52vWrFnyfxcVFT31fDdv3kReXh7c3NxQt27dykQjIiIiqpBKPVbbs2fPE+dPnToFAHBzc4Ozs/NTz1f8SM3HxwcqlQo//PADLly4AFtbW7Ro0QLDhw/H22+/zYUbiYiIqMoYfYXsvLw8pKSk4MSJE5g9ezYAYNq0aeX6vcUvY588eRLHjx8XnPPMmTM4c+YMoqKisGnTJtjY2Bg7OhEREZmZnJxsZGemo2Yd4+2VatRytHbtWkydOrXka3t7eyxduhRjxowp1+8vvnOk0+nw3nvvYcKECXB3d8etW7ewfv16rFixAlFRUViwYAHmzZtX5nkCAwPL9f0SExPh4eFR8niQjCstLU3sCFQGXhvp4rWRLl4badFo8nBq9zeo/9c3sMu7C8B45ahKn0/l5ubi8OHDuHPnTrmOLywshK+vL8LCwrBq1So0a9YMtra28PDwQHh4eMkn3lavXo309PSqjE5EREQSVFBQgGM7v8GVua3RKnoh6hXeNfr3sMjOzn76m9LPKCUlBWfPnsXcuXORkJCAZs2a4ffff4dCoajUefPy8tC0aVNkZmZi48aNGDRoUKXOFxgYCEtLS8TFxVXqPGRY8R25xzf8JfHx2kgXr4108dqIS6vV4uf/rIHtr4tQv+CWYO6NPYDNc80RHR1tlO9VJXeO3Nzc0L9/f0RFRaFWrVpISEjArl27Kn1eOzu7kkdmxY/giIiIqPoq1Bbi5+0qHB3bFE1+/livGFWFKn2s1qBBA3Tq1AkAEBsba5RzFn/qLS8vzyjnIyIiIunR6XT4Ze8mHBrXHA1/GoNGmhsGj0uo+TwKnBoY9XtXuBzdvn0brq6ucHV1RXJycpnHubg8ekEqKyvriedLS0vD4cOHcezYsSce9+DBAwAo19IAREREZF50Oh1OHIzEvnH+cNvxLprkGX5SdNXeF2lv/gf9VvwXNewdjZqhwp9Wc3V1RVFREXJzc3HlyhU0bNjQ4HHFH89v0ODJrS4rKwtBQUEAHq287ePjo3eMVqstWZU7ICCgotGJiIhIgn4/tg8p22fDO+sc6pVxzHU7Lyj6fIbX+71TZeseVvisVlZWJY/Mvv/+e4PH/PnnnyULQXbr1u2J5/P09ESrVq0AAAsXLjR4zPr166FWq9G4cWO88sorFY1OREREEnLujxPYMbEdnL7rC++scwaPSbZtjJt916LX6kvo1v/dKl0QulJnHjt2LAAgMjISI0aMwN9//42cnBzcvHkT3333HYKDg6HT6fDaa68J1h4KCAhAQEAA5syZIzhf8WKRO3bsQGhoKC5cuACNRoOkpCSEh4cjLCwMADB//nxYWVlVJjoRERGJ7OL5/+LHaT1hu7ozfO+dMXjMHZv6uN5zObquuYpeg9+HlXXV//yv1CKQXbt2xdSpU7FkyRJERkYiMjJS75g2bdpApVIJxooftaWkpAjGg4KCMG7cOKxevRpbt27F1q1b9c43Z84cBAcHVyY2ERERieh64iWc3RAGv+Sf0BI6g8eorV3woP0UvDZ0ImwruRTQs6r0Ctlz585Fhw4dsG7dOkRHRyMjIwOOjo7w8/NDSEgI3n333Wfa6mPRokXo0qULVCoVYmJicP/+fdStWxft2rXDRx99hLZt21Y2MhEREYkg5fZNHF8/Ez5Xt8EfWoPHZFjWRtrLE9B7+HTY1bA3ccJHjLJ9SPfu3dG9e/dyH5+dnf3E+d69e6N3796VjUVEREQSkHE3DUfWfwrPC9/Dv8jwUjwPLByR/MIY9Br5GWo61TJxQiGjbzxLREREBABZD+7j0HcL0PDcWrTUGb4x8tBCgavN30W3UfPR1sXVxAkNk005UqvVevuxaTQa2NraipSIiIioetLk5eHA94vgfHYF/AozDB5TAGtc9ByCjqPC8VLDJqYN+BSyKUcqlQoRERF648WLVBIREVHlaLVaHPz3CtgfXwxvbYrBYwphifgG/dB6ZARCmvqaOGH5yKYchYaGliwyWSwkJIR3joiIiCpJp9PhyK7vUHRoPjzL2OYDAOJcuqHl8EUI8X/JhOmenWzKkVKp1NtJWaFQVOkiUkRERNWZTqfDyUPbkb3nM3g8vFTmcRfqtIXHmxEY3LaL6cJVgmzKERERERnPn7/9jFvbZsA76y+U9YLKFUd/1AsOR3DXfibNVlksR0RERFRu8f9EI+77aWiRfgzeZRxz3a4pavSbh7593jTLJzQsR0RERPRUSdcTcEYVBr/kPWhRxqrWt2wborD7LPQKDjXJNh9VheWIiIiIynQ3LRVHVDPR7NJm+CPf4DFqaxdkdZyG196eAJtq8EEnliMiIiLSk5OTjYMbPkeD/66Bf5HhBRwzLZ2QEvgRXhsxE/b2DiZOWHVYjoiIiKhEQX4+Dmxehlq/LYFvYbrBYx5aKHDVdwR6hM5Hu7r1TJyw6rEcEREREXQ6HY7u3oiiA/PQND/J4DEFsMLFJsHoMPoLvNTIw8QJTYfliIiISOZOH92DuztmwSsnrsxjzrv2QMDIJQhp3sqEycTBckRERCRT/8ScRsKmafC9dxq1yzjmUq3WcH/7SwwxkwUcjYHliIiISGYSr1xA9Pqp8LtzAL4oMnjMtRreqDlgAYJeG2zidOKTTTlSq9VITxe+WKbRaLi3GhERyUbqnWQcU30Cn6vb0BJag8fcsamPgu6z8dqgMWa9VlFlyKYcqVQqRERE6I27uJS16DkREVH1kJ31AAfWz0GTv7+Bf9FDg8dkWNVBettJ6P3uNCjs7EycUFpkU45CQ0MRFBQkGAsJCeGdIyIiqra0Wu2jj+WfWIQWhXcNHpNjUQPXW41Br1Fz4VSrtmkDSpRsypFSqYRSqRSMKRQKs9zzhYiI6El0Oh2O798Gzd7Z8NJcM3hMPqxxqemb6DI6Aq2fa2jihNImm3JEREQkB3+dPY7rm6fA50GMwXkdLBD/XG8Ejl6CkKa+Jk5nHliOiIiIqoFrVy7iz/VT0PLOfviUcczFWq3hNewrDGnd0aTZzA3LERERkRm7m5aKI+s+gfflzWiJAoPH3LDzgv2AhRjYe4iJ05knliMiIiIzlPcwF/u/Dcdz0Svhr8syeEyaVT1kdfkEr739Mayt+SO/vPgnRUREZEZ0Oh0O/bgWiiPz0bwgxeAx2Rb2SHrhQ/QOnQsHB0cTJzR/LEdERERm4vTRPcjYPgOeuRcNzhfAChe93kSX0C/Qhp9AqzCWIyIiIomL//tPxG+cAr+Mk2XugRbn0g0vjF4qi41hqxrLERERkUTdTr6Ok99MhV/SLvhBZ/CYK47+eO6tJRjcoaeJ01VfLEdEREQSk/XgPg6qPoVHrAr+RXkGj7ll2xDoPQ99B7zHBY2NTDbliBvPEhGR1BVqC3Fw6wo4/roQLQrTDR6TYVkb6e2moM+702CrUJg4oTzIphxx41kiIpKy00f34F7kdHg8vGxwPg8KXPEbgV7vL0St2nVMnE5eZFOOuPEsERFJUcLFf3BONQEt0o8ZfNm6EJaIb9gf7d9fghcbe5o6nizJphxx41kiIpKSu2mpOLJ2Oppf+TdaQGvwmEu1WsNr+NcIeamdidPJm2zKERERkRTkazTY92043P74usyVrZNtG8Om/0L07/MW/xEvApYjIiIiE9DpdPhl7yYU7ZsNn/xkg8dkWtaC+pUp6DN8Ol+2FhHLERERURU798cJXNs0CT4P/mtwPh/WuOQ9DD0+WIR2zvygkNhYjoiIiKpIclIiTq2djBa39sIHRQaPiVN2R0DoMoR4tzRxOioLyxEREZGRFS/i6Bm7Dv5FGoPHJNo3R703l2Jw594mTkdPw3JERERkJOVZxFFt7YLcbp+i95CxsLK2MnFCKg+WIyIiIiM4e/wg0rZNhGfuJYPzORY1cOP5D9B7zOdwcHA0cTp6FixHRERElXAjMQFnvxmPlimHUNPAfCEsEd94ADq+/xVaN2xi6nhUASxHREREFZCd9QAHvpkFr1gVWsLwe0VcxNE8sRwRERE9A51Oh0M/rkWNw3PRQptm8Jhk20awfmMh+r/+NhdxNEOyKUdqtRrp6cKX4zQaDfdWIyKicov/+3c82D0bTbNjDc4/sHDE7TYT0XfUbC7iaMZkU45UKhUiIiL0xl1cuNgWERE92e3k6zi28iO0StkPpYH1igphifgmg/Dqh0vR1q2BCAnJmGRTjkJDQxEUFCQYCwkJ4Z0jIiIqU97DXEStmwP3v1bjhaKHBo+55PQSmo5ciZCAV0ycjqqKbMqRUqmEUqkUjCkUCj4LJiIiPTqdDkd3b4Tl/tnwK7hj8Jg7Ns9B12cB+g94jz9LqhnZlCMiIqLyiP3rDK58Ox4+D2IMzudY1MCNgHHoO2Ye7GrYmzgdmYIkq25MTAyGDx8Ob29v1KlTBw0aNECPHj2wYcMGFBYWih2PiIiqIXXKLfxnTgjwdQeDxUgHC5xzex3158dh0ITFLEbVmOTuHEVGRiI0NBRarbZk7P79+zhz5gzOnDmDn376CZGRkbCxsRExJRERVRcF+fmI+nYhnvv9K/gXZRs85oqjP2q+8Tm6B7TTe0WDqh9J3TlKT0/H+PHjodVq0bp1a/z888+4ffs2Lly4gHnz5sHW1hZHjhzBl19+KXZUIiKqBk4d2YOjHzWHz5l5cDJQjNKs6uFG79Xou/IcWgRwIUe5kNSdo+3btyM7Oxv169fHvn37YG//6Jalk5MTpkyZgvz8fISHh2Pjxo2YNWuWyGmJiMhc3UhMwNk149BSfRh1DMznQYGr/qHo/X44HGs6mTwfiUtSd45iYh494+3Xr19JMSpt4MCBAIA7d+4gIyPDpNmIiMj85ebmIHLZJKR/7o+W6sMGjznv1gvOc2IxeOpKFiOZktSdI7VaDQBo0qSJwfmaNf9/S7+iIv1FuIiIiMryy95NKNo7Ay0Kbhucv27XFHXf+hpDurxu4mQkNZIqR3v27Hni/KlTpwAAbm5ucHZ2NkUkIiIyc5fjzyF23Tj43jttcP6BpSNS2oWh74gZsLaW1I9FEonk/1eQl5eHlJQUnDhxArNnzwYATJs2TeRUREQkdQ/uZ+Lgv6bB++JG+EKrN6+DBeLcg/Hq2K+55QcJSLocrV27FlOnTi352t7eHkuXLsWYMWOe+PsCAwPLdf7ExER4eHiUPM4j40pLM7xbNYmP10a6eG0qT6fT4dT+Taj32yK0LLxr8JgE+xaoPTACrz7/MgCU6+cAr410abVaoy7xI+ly9Ljc3FwcPnwY/fr1w3PPPSd2HCIikphL52OQvmMGfLP/MTh/16oOUl+Zhg79RnDLDyqTRXZ2tuTfbE5JScHZs2cxd+5cJCQkoFmzZvj999+hUCgqdd7AwEBYWloiLi7OSEmptOJ/iXHBNOnhtZEuXpuKSU9LwdHVk+B37UdYQac3XwArXPR+F73GLkbtOhV7Z5XXRrqaN28OCwsLREdHG+V8ZlGb3dzc0L9/f0RFRaFWrVpISEjArl27xI5FREQi02q12LMhApen+8D/2jaDxehi7ZdhO/UsQmZ9W+FiRPJiFuWoWIMGDdCpUycAQGxsrMhpiIhITDGnf8HB8S3R7MRM1NY90JtPsXbDrf4bMGDZafj6vyRCQjJXknnn6Pbt2wgICADwaDHIhg0bGjzOxcUFAJCVlWWybEREJB3paSk4unI8WtzYAU/ovxmSB1tcbfU++nywEA4OjiIkJHMnmXLk6uqKoqIi5Obm4sqVK2WWo4SEBACP7iIREZF8FGoLsX/zUjgfD4e/7r7BY+JcuiLwg9UY3LS5idNRdSKZx2pWVlYlj8y+//57g8f8+eefJQtBduvWzWTZiIhIXH//eRL7Pn4BXr9OR20DxSjZtjHS347E4CVH4cFiRJUkmXIEAGPHjgUAREZGYsSIEfj777+Rk5ODmzdv4rvvvkNwcDB0Oh1ee+21cq9lRERE5uteRjr+M28oLFd1QdOc83rzDy0UiH9xKjqvvIROvQaJkJCqI8k8VgOArl27YurUqViyZAkiIyMRGRmpd0ybNm2gUqlESEdERKai0+lwcOsqOB2dC//CewaPiXN5Fa3HrsFLnj4mTkfVnaTKEQDMnTsXHTp0wLp16xAdHY2MjAw4OjrCz88PISEhePfdd426CiYREUnL+b9+x5UNH8I765zB+Ts29WEV9BUGv/6maYORbEiuHAFA9+7d0b17d7FjEBGRCd3PvIeDqybBN2ETvA2sV/ToU2gf4PUPF8Le3kGEhCQXkixHREQkHzqdDj9HroP9odnwL2MvtPi6HfDCh2sw2LulidORHMmmHKnVaqSnpwvGNBoNbG1tRUpEREQXz/8XF9Z9AJ/7fxqcT7F2Q9EbizCo/7smTkZyJptypFKpEBERoTdevKgkERGZTk5ONvatmgqf+A3wgVZvPh/WuOw3Gn3GLoJjTScREpKcyaYchYaGIigoSDAWEhLCO0dERCZ2bN82FO6ajJYFdwzOX6z9Mlq8vwZD/AJMnIzoEdmUI6VSqbeTskKhgKWlpJZ6IiKqtm7dvIZTK95HS/Vhg/NpVvWg6bMQAwaO4t/NJCrZlCMiIhKHVqtF1IZw1D/zJVoW5erNF8AKF32Go/e4r+BUq7bpAxI9huWIiIiqzLk/TiD5u/fhnXvR4PwVR394jl6HkIC2Jk5GVDaWIyIiMrr7mfdwcOXH8Lvyb3gaWLPogaUj0jp9iteHTYGVtZUICYnKxnJERERGo9PpcHT3RthGhcG/MN3gMbH1+6LLx2vQ9rmGJk5HVD4sR0REZBTXrlxEzOox8Ms4aXA+2bYR7ENWIKT7ANMGI3pGLEdERFQp+RoNflr3GZrELIdfkUZvXgMbXGn1IfqNjYBdDXsREhI9G5YjIiKqsD9OHsLdzWPhm5docP5SrdZo8cE6DPZ7wbTBiCqB5YiIiJ7Z3bRUHFn5EfxvbIejgfkMy9rI6rkA/UM+5JpFZHZYjoiIqNyKN4l1PDQT/oX39OdhgTj3Qeg+fiWcXVxFSEhUeSxHRERULjcSE/DnypHwy/jN8LydF+oN+xdCOvQ0cTIi45JNOVKr1UhPF36sVKPRcG81IqKn0Gq1+Gn9AjT8/Uv4FT3Um39oocD1wEl4Y8w82PDvVKoGZFOOVCoVIiIi9MZdXFxESENEZB5i/zqD6+tGwyc33uD8hbrt8dL4DXjJ08fEyYiqjmzKUWhoKIKCggRjISEhvHNERGRAbm4OolZMhs+FDfBCod58hlUdZPdaiKDBY/jCNVU7silHSqUSSqVSMKZQKPj/1EREjzl1ZA9yt41Dy4JbBudjGwWh+4Q1fOGaqi3ZlCMiInqyu2mpOPL1h/BP3oU6BuZv2TSA/ZurEdK9v8mzEZkSyxERkcw97eP5BbDCJb/R6Dv+K9jbO4iQkMi0WI6IiGTseuIlRK8cXebH86/a+6HJmPUYEvCKiZMRiYfliIhIhrRaLX5SfY6GZ5cY/Hh+jkUNJLedjn6jZ8Pamj8qSF74v3giIpmJ//tPXP1mBHxy4gzP1+2A1uO/RWvPZiZORiQNLEdERDKRr9Fgz+owNPt7Nbyg1ZvPsKqDnN4RGBgcyk/ykqyxHBERycC5P07g9oZRaJF3xeA8P55P9P9YjoiIqrFHizlOQvMLG9AEOr352zYNUIMfzycSYDkiIqqmzh4/iPub30fL/CS9uUJYIt57OPpOXAEHB0cR0hFJl2zKETeeJSK5yHpwH/u/HocWV/+NmijSm09SeMDlvfUIaddVhHRE0iebcsSNZ4lIDk4e2oGCyI/gX5CiN1cAa1z2/wBvfLQYCjs7EdIRmQfZlCNuPEtE1dm9jHT8vHQM/G/uMjifaO8D9zHfYTAXcyR6KtmUI248S0TV1dE9P8Bmz2T4F97Vm8uDLa4FTsQb78+HDf8xSFQusilHRETVjTrlFo4tG42WKQcNzifUfB7Nx36HYL8AEycjMm8sR0REZujgf9ai5sEZaKm7rzeXY1EDt9p9gn4jZ8LK2kqEdETmjeWIiMiMpN5Jxomv3kOLtKMG5y/WfhkB4zeiddPmJk5GVH2wHBERmYlDkevguH86Whi4W/TAwhHpXedhwDsT+S4lUSWxHBERSZw65RaOLR2BlqmHDc7H1euCdhO/RdtGHiZORlQ9sRwREUnY4Z3fosZPU9BSl6k3l2nphAe9vkDwkPd5t4jIiFiOiIgkKC31Dn5dOrLMT6LF1euCDpO/x3MNGps4GVH1x3JERCQxZw7vQN2js9GyMENv7oGlIzJ6LETwm+N4t4ioirAcERFJxN20VBxd8h6eTzV8tyjeuSNemfQ93y0iqmIsR0REEvDL3k2w3j0Rzxu6W2ThiIwe4Rj41ke8W0RkArIpR2q1Gunp6YIxjUbDvdWISFQZd9NweMko+N/+yeB8fN0OaDvpe7Rt7GniZETyJZtypFKpEBERoTfu4uIiQhoiIuDYvm2w2DEe/oXpenNZFg5I6/o5BnLdIiKTk005Cg0NRVBQkGAsJCSEd46IyOQe3M/Ega/GwP9GpMH5+Fpt0GbyZrzcpJmJkxERIKNypFQqoVQqBWMKhYL/IiMikzp7/CCyN42Cf8FtvblsC3tcbzsDnQaEws3NTYR0RATIqBwREYlJk5eHPcvGw/fit6gJnd78xdov46UJm9DEsZYI6YioNJYjIqIqFvvXGdxcOwwt867qzT20UCC5w6cYMPITWFpaQq1Wi5CQiEozyjOl1NRUzJ49G4GBgSWPr9q0aYP58+cjMzPzmc712WefwdHR8Ym/hg8fbozYRERVSqvVYseK6Sj8uhOaGChGVxxawmnGH+g/ehYf8RNJSKXvHMXHx6Nfv35ITU3VG4+Pj8e2bduwZ88eNG3atFznS0hIqGwkIiLRXbl8HnHLh8I3+x+9uQJYIyFgAgZ89AWsrXkDn0hqKvVPlaKiIowYMQKpqanw8vLCzp07oVarcfXqVahUKiiVSty4cQNDhgxBQUFBuc555coVAMChQ4eQnZ1t8Nf3339fmdhERFVGp9Nh77eLcG9hGzQzUIxu2HkCH/+KQROXsBgRSVSlytEvv/yCuLg42NjYYNeuXejZsyfs7e3h6uqKt956C0eOHIG9vT0uX76M3bt3P/V8Op0O165dAwD4+vpWJhoRkcndTr6OXVM6oenxGXAoeiiY08ECsT4j8erXsXj+pQ4iJSSi8qhUOfr1118BAF26dIGnp/7qrZ6eniVrC506deqp57t58yby8vLg5uaGunXrViYaEZFJHYpch+ufPg/fDP2/61Js3JA94ieEzNwAuxr2IqQjomdRqXJ09eqjFwyfdJeneG2hnJycp56v+JGaj48PVCoVOnbsiHr16qF+/fro0aMHNm/eDJ1O/yOwRERiybibhh9nvI5GUe+jtu6B3nxs44F4aUk82nZ5XYR0RFQRlXrg/eGHH2LgwIHw8/Mr85hz584BANzd3Z96vuKXsU+ePInjx4+XjOfl5eHMmTM4c+YMoqKisGnTJtjY2FQmOhFRpZ06sgcFW0PRUpumN5dhVQf5QSsQ0u8dEZIRUWVUqhx16tTpifNHjhwpefTWt2/fp56v+M6RTqfDe++9hwkTJsDd3R23bt3C+vXrsWLFCkRFRWHBggWYN29emecJDAwsV/7ExER4eHhwXZEqkpam/wODpIHXpnI0mjyc+PZTPJ+4CZYo0pv/x7kLWo1agXours/89wuvjXTx2kiXVqs16k2TKltYY8uWLRg6dCgAIDg4GC+88MJTf09hYSF8fX0RFhaGVatWoVmzZrC1tYWHhwfCw8Mxa9YsAMDq1auRnq6/USMRUVW7ejkWfy/ojoDEH/SKUZaFA+I7LkKX6VtRz8VVpIREVFkW2dnZ+v/sqYTY2FiEhYXhxIkTAIAOHTpg165dqFGjRqXPnZeXh6ZNmyIzMxMbN27EoEGDKnW+wMBAWFpaIi4urtLZSF/xv5gf39OOxMdr8+x0Oh32bliIxr8tgB00evOXnF7Ci5O2wt2zcpvF8tpIF6+NdDVv3hwWFhaIjo42yvmMdufo/v37mDRpEtq3b48TJ07AxsYGs2bNQlRUlFGKEQDY2dmVPDIrfgRHRFTVUm7fxM4pneH926d6xSgf1rjQegbeWHa20sWIiKTBKCuQRUdHY9iwYbh58yYAoF+/fliwYAG8vLyMcXoBZ2dnAI/uIhERVbVfo7bAZsdH8NNl6s3dVLjD7f3NCOa6RUTVSqXL0fHjxzF48GDk5ubC3d0da9aseeqL2oakpaXh3LlzsLGxQZcuXco87sGDRx+VLS5JRERVIScnG1GLQ+F/bZvB+ViPt9B32jo4ODiaOBkRVbVKlaO7d+9i2LBhyM3NRefOnbF161Y4OTlV6FxZWVklC0bGxMTAx8dH7xitVouYmBgAQEBAQMWDExE9wd9/nkSKahj8NTf05jKs6qBg4CqE9H1bhGREZAqVeudo3bp1yMjIgLu7OyIjIytcjIBHq2m3atUKALBw4UKDx6xfvx5qtRqNGzfGK6+8UuHvRURkSKG2EDtWfgKs6opGBopRfN2OaLrgb7zKYkRUrVWqHO3btw8AMHLkSNjbl39J/ICAAAQEBGDOnDmC8WnTpgEAduzYgdDQUFy4cAEajQZJSUkIDw9HWFgYAGD+/PmwsrKqTHQiIoGk6wnYO+ll+EZ/ARtoBXN5UCChUzgGfnUMbvUbiZSQiEylwuWooKCg5CPwc+bMgaOj4xN/TZ8+veT3JiQkICEhASkpKYJzBgUFYdy4cQCArVu3onXr1nB2doafnx8iIiJQWFiIOXPmIDg4uKKxiYj0HN75LW7PexE+D2L05q7V8EaNqafQf9RMWFpW2dJwRCQhFX7nSK1Wo6CgwJhZAACLFi1Cly5doFKpEBMTg/v376Nu3bpo164dPvroI7Rt29bo35OI5CknJxtRX46C//Uf9eZ0sEB881HoP2klFHZ2IqQjIrFUuBw1aNAA2dnZFfq9T/t9vXv3Ru/evSt0biKi8og79wdu/utN+Guu6c2prV1g+7YKQ7r1FyEZEYnNKOscERGZC51Oh5+++wKNTsxDY+TrzZ9X9kD3GVtQ19lFhHREJAUsR0QkG3fTUvHLorfRIu0XvbmHFgrc7vw5Bg2fyneLiGRONuVIrVbrbVar0Whga2srUiIiMqXfj+1D3qYRaKHV31n9hp0X3MdtQ79WgSIkIyKpkU05UqlUiIiI0Bt3ceGtc6LqTKvVYtfyyWj+z2o4Qac3H+v5NvpNWwd7ewcR0hGRFMmmHIWGhpaswF0sJCSEd46IqrGk6wn475IhaJF1Tm8u09IJDwesREj/d00fjIgkTTblSKlUQqlUCsYUCgXfLSCqpg7v/BaOP02Ety5Lb+5yzQC8NO1HNHJvKkIyIpI62ZQjIpKHR2sXjYb/9f/ozWlhiUsvfIyg8Ythbc2//ojIMP7tQETVxoXYGFxfNQT+eYl6c6nWSti/+z0Gd35NhGREZE5YjoioWti/ZSWUh6fDvShPby7OpRu6hm2Bs4urCMmIyNywHBGRWcvOeoB9X7wH/+RdenN5UCC581wEvzed7xcSUbmxHBGR2boQG4MbKwcb3AIkSeGBxuP+gzeeby1CMiIyZyxHRGSW9m3+Gq5HPkFjA4/RYhsF4fWwjXCs6SRCMiIydyxHRGRWHj1Gexf+yXv05nIt7KDusQghQz8WIRkRVRcsR0RkNuL//hNJq4fAX3Ndb+6GnSeafPQj+vi/ZPpgRFStsBwRkVnYt2kZ3I5+gsZFGr252MYD0XfG93BwcBQhGRFVN7IpR9x4lsg8ZT24j/1fvAv/W3v15nIsaiCt5yKEvD1ehGREVF3Jphxx41ki8xP/95+4uXow/DU39OZu2HnBY/yPaN3yRRGSEVF1JptyxI1niczL/i0r4Xp4GhoZeozmPgh9w77jYzQiqhKyKUfceJbIPOQ9zMWeL0bA//qPenM5FjWQ3msxQt4aJ0IyIpIL2ZQjIpK+a1cuIu6rIPjnXtSbu27nBa8JkWjtFyBCMiKSE5YjIpKEY/u2QbF9DDx1WXpzsY0Hot8nP8De3kGEZEQkNyxHRCQqrVaLncsmwO/8GliiSDCXBwVud1uIkHcni5SOiOSI5YiIRKNOuYXfIoLRMvOs3txtmwZQfvgj+r7UToRkRCRnLEdEJIo/fzuMnO+Gork2TW8url4XdJ8ViTp164mQjIjkjuWIiExKp9Nhj2o+PE8vgAO0gjktLJHQejqCx4bzk6REJBqWIyIymQf3M3EwPAQtU3/Wm0u3qgubYd8j+NW+IiQjIvp/LEdEZBLx/0Tj5qpBaGlgtevLNQPQZsZO1G/YxPTBiIgew3JERFXu5+0q1ImagEZFD/XmzvuMRtDU1bDhavVEJBEsR0RUZbRaLXYu/hAtL67Xm8uycEDWgNUYMmC4CMmIiMomm3KkVquRnp4uGNNoNNxbjaiKpKXewcmFQQY/pn/drim8J+3Ey839RUhGRPRksilHKpUKEREReuMuLi4ipCGq3v7+8yTurh2C5toUvbnY+v3Qd/a/uWksEUmWbMpRaGgogoKCBGMhISG8c0RkZPu3rITbz9PgBo1gvABWuNZ+DgaPnsWP6RORpMmmHCmVSiiVSsGYQqHgX9JERpKv0WDXotHwv7pZb+6uVV3YDN+MAZ17i5CMiOjZyKYcEVHVuXMrCWcjBsA/6y+9uSsOLfHijD1o2NhThGRERM+O5YiIKiX69BHkrH8L3oXpenOx7oMxYOYPUNjZiZCMiKhiWI6IqMJ++u5LNDo2Cy6PbQOigQ2SXw1HyHvTREpGRFRxLEdE9Mw0eXnYvXAY/G9s15tTW7ug5qht6NeuqwjJiIgqj+WIiJ5Jyu2bOBveF/7Z/+jNXa4ZgLYz98CtfiMRkhERGQfLERGV27k/TiDzm0Fopk3Tm4tt+i4Ghqm4DQgRmT2WIyIql4P/WQvl/olQPrZ+UR4USOm1GCFvjxcpGRGRcbEcEdETFWoLseOrcWgZ/43enNraBbXf344+bTqJkIyIqGqwHBFRmTLv3cWR+UFoefek3lyCYyu0mfkTnmvQWIRkRERVRzbliBvPEj2bKxdjcXVpP/hpbujNxTYeiAGztnD9IiKqlmRTjrjxLFH5nTy0Azbb3kMjXbZgXAtLJLb7DINDP+XWO0RUbcmmHHHjWaKn0+l02PPNXHj9Hg5r6ARzmZZO0A39HgO6DxAnHBGRicimHHHjWaIny3uYiz0L3oZ/8h69uSSFB5pNjYKXt58IyYiITEs25YiIypZy+yb+WNAH/jnn9ebi63VGj093oVbtOiIkIyIyPZYjIpk7/9fvSFs1AE21qfpzLccieNIKWFlbiZCMiEgcLEdEMvZr1BY4bg+Fa9FDwfhDCwXu9l2JIYNCRUpGRCQeSb5wk5qaitmzZyMwMLDkXaE2bdpg/vz5yMzMFDsekdnT6XQ4umUx6kW+C4fHilGqtRLW44+gJ4sREcmU5O4cxcfHo1+/fkhNTdUbj4+Px7Zt27Bnzx40bdpUpIRE5i1fo8HRFWMRcGuX3txVhxYInLWfCzsSkaxJ6s5RUVERRowYgdTUVHh5eWHnzp1Qq9W4evUqVCoVlEolbty4gSFDhqCgoEDsuERm525aKvZN62CwGJ1364VuS86yGBGR7EmqHP3yyy+Ii4uDjY0Ndu3ahZ49e8Le3h6urq546623cOTIEdjb2+Py5cvYvXu32HGJzMrl+HP478xA+NyP1puLa/UxBkXsh729gwjJiIikRVLl6NdffwUAdOnSBZ6ennrznp6eJQs5njp1yqTZiMzZqSN7kLG4IxrkJwvG82CLpD5rMHjKcq75RUT0P5J65+jq1asAAF9f3zKPKV7IMScnxySZiMzdT999icbHZsIGhYLxu1Z1UGP0j3itXXeRkhERSZOkytGHH36IgQMHws+v7FV4z507BwBwd3cv85jAwMByfb/ExER4eHhArVY/U04qn7S0NLEjyFphYSF++WY6Aq5t1pu7pvBEjSH/gk/TVvzfv8Tw/2+ki9dGurRaLWxsbIx2PkmVo06dOj1x/siRIyWP3vr27WuKSERmKScnC3+sGIWAjON6c+frtIfH21/B3sFRhGRERNInqXL0JFu2bMHkyZMBAMHBwXjhhRfKPDY6Wv+FU0MCAwNhaWmpt+caGRf/fE0r5fZNnP/yDfjnxuvNxXqPQHDYOmRkZADgtZEyXhvp4rWRHmtr49YZyZej2NhYhIWF4cSJEwCADh06YO3atSKnIpKmC7ExuLO8L7wKUgTjBbBC0qsRCHlvmkjJiIjMh2TL0f379zF37lxs2LABOp0ONjY2mD59OqZNm2b0hkhUHZz+5SdYbnobbrpswXiWhQMK3t6Efj2DREpGRGReJNkyoqOjMWzYMNy8eRMA0K9fPyxYsABeXl4iJyOSpv3/Xon6hybDFlrBeKq1K1w//gl+z7cWKRkRkfmRXDk6fvw4Bg8ejNzcXLi7u2PNmjVPfVGbSK50Oh12LJ+CFue+1ptLtPfBi7MOon7DJibPRURkziRVju7evYthw4YhNzcXnTt3xtatW+Hk5CR2LCJJytdosOvzEPgn79Gbi3fuiF5zf0JNp1oiJCMiMm+SKkfr1j36FI27uzsiIyNhb28vdiQiSbqXkY5f5/aB//0/9eZiPd9G8Kzv+W4eEVEFSepvz3379gEARo4cyWJEVIYbiQmIX/QamuclCsZ1sMDll2chZOx8kZIREVUPktlMqaCgAHFxcQCAOXPmwNHR8Ym/pk+fLnJiItP7O+Y33FjQFu6PFaM8KHCn/3oMZDEiIqo0ydw5UqvVKCgoEDsGkWT99vMu2P77HdQryhWM37OsBdtRkejRoYdIyYiIqhfJlKMGDRogOzv76QcSydCBbf9C/QMTYPPYR/WTbRvDc9o+NPVuKVIyIqLqRzLliIgM27l6Fpr/sVBvPKHm8+gw5xCcXVxFSEVEVH3Jphyp1Wqkp6cLxjQaDWxtbUVKRPRkhdpCbF8UCv/L3+nNxbl0xeuf74W9vYMIyYiIqjfZlCOVSoWIiAi9cRcXFxHSED2ZJi8Pe+YOhP+dA3pzsR5vInj2Jn5Un4ioisjmb9fQ0FAEBQn3lgoJCeGdI5Kc+5n38MtnvdDSwBpG8QGTMfjjxbC0lMwHTYmIqh3ZlCOlUgmlUikYUygU/CFDknI7+TrOLegJn4cJgnEtLJHUbTEGvTtZpGRERPIhm3JEJHWX4v/CraV90KQgRTCea2GH7MHfoe/rb4qUjIhIXliOiCTgz98Oo2DDYLjp7gvG71nWgl3oTnRp11WkZERE8sNyRCSy4/v/A8cfh6N2kUYwfsfmOTSeehDNmrcSKRkRkTyxHBGJ6OC2NXjuwMd6izteq+GNFz89jOcaNBYpGRGRfPFtZCKR7FkfjoYHxukVo4u1X0aHL35nMSIiEgnvHBGZmE6nw47lU9Di3Nd6c+fdXkP/ebugsLMzfTAiIgLAckRkUoXaQmyPGAn/Kz/ozcV6vIVBszfBytpKhGRERFSM5YjIRAry87FzziD43/5Jby6u5TgMnrKC624REUkAyxGRCeTm5mD/7D7wv3tCb+5S208x+MPPRUhFRESGyKYcceNZEkvmvbs49mkP+GX9JRjXwhLJPZchaOjHIiUjIiJDZFOOuPEsiSH1TjKi53WH98NLgvE82CIzeD36vDFMpGRERFQW2ZQjbjxLppZ0PQEXw7vBI/+mYDzbwh7ad7eha9d+IiUjIqInkU054sazZEpXLp9H0qIeaKgV7pOWYVkbjh/uQZs2nURKRkRETyObckRkKhdiY5C2rBfcCu8KxlOtXdFgykF4+70gTjAiIioXliMiI/o75jdkr+qLeo9tIJts2xi+s4+ikXtTkZIREVF5sRwRGUn06SPQrgtCnaJswfh1u6Z4ae6vcH2uoUjJiIjoWfCFGyIjOP3LTyha9wacHitGV+398PKC31iMiIjMCO8cEVXS8QM/oua2d2EHjWD8cs0X8OqCX1Crdh2RkhERUUXwzhFRJRzZtRG1tg3VK0YXa7+M7hHHWYyIiMwQyxFRBR3Y9i8od4+CLbSC8fh6ndHni1/gWNNJpGRERFQZLEdEFfDTxsVodOAjWEMnGD/v1gtvLDwEuxr2IiUjIqLKYjkiekZ71ofD69fpsESRYDy2YRAGhkfBVqEQKRkRERkDyxHRM9j9zTw0OzlbbzzWaygGzYuEtTU/40BEZO5k8ze5Wq1Genq6YEyj0XBvNSq3XWs+g8/v8/XGz/uGYvD0tdyKhoiompBNOVKpVIiIiNAbd3FxESENmZudq2eh+R8L9cbj/D/CkKkrRUhERERVRTblKDQ0FEFBQYKxkJAQ3jmip9qx8hP4Rn+hNx73/AQMnvy16QMREVGVkk05UiqVUCqVgjGFQsFHIfRE25dPg99/l+iNxwdMweCJ+uNERGT+ZFOOiJ7V9mWT4Xdumd74hcDpGDR+kQiJiIjIFFiOiAyI/GoCWvyzQm/8YptPEDxO/90jIiKqPliOiB4TuWQ8WsSu0hu/+PJsDByr/2k1IiKqXvjCDVEpkV99bLAYXXplDosREZFM8M4R0f9s/3oKWvyj/7H8y+0/R9CYT0VIREREYmA5IsKjj+v7/bVUbzyhUzgGjJopQiIiIhILH6uR7O3816cG1zFK6LgA/VmMiIhkh+WIZG33uvlofnaB3vilV+ag/+hZIiQiIiKxsRyRbO399gt4n/pMb/xim08Q9MFc0wciIiJJkM07R9x4lkqL+n4pPI/rPzKLf3EqBnEdIyIiWZNNOeLGs1Rs/5YVaPLLVFiiSDAe1+pjDJ6wWKRUREQkFbIpR9x4lgDg4LY1aPTzJL1idL7FhxgyZblIqYiISEpkU4648Swd3vkt6h/4CFbQCcbP+4zGoKn6Cz8SEZE8sRmQLJw4GAmXPe/D+rFiFNtsOAbN+IYlmYiISvAnAlV7vx8/AIdtw2ALrWA81uMtDJ75LYsREREJ8KcCVWvn/jgBy+8GoUaRRjAe2ygIgz/bzGJERER6jPqTQafTwcvLC0OHDq3Q7//ss8/g6Oj4xF/Dhw83ZmSqxuL/iUbumjfgWJQrGD/v2hPBc39kMSIiIoOM+tPh0KFDSE1NrfDvT0hIMGIakrOrl+OR9nVv1NbdF4zH1+2AAfP3wtpaNp9FICKiZ2S0nxBXr15FWFhYpc5x5coVAI9KVvv27Y0Ri2QoOSkR17/sjucKhYt+XnJ6Eb0XHICtQiFSMiIiMgeVKkfnzp3Dli1bEBMTg+joaOh0uqf/pjLodDpcu3YNAODr61uZWCRj6pRbiFvQFY0K7gjGr9r7oev8I3BwcBQpGRERmYtKlaNTp05hzZo1Rgly8+ZN5OXlwc3NDXXr1jXKOUleMu/dxZ9zu8FDc0MwfsPOC698fhS1atcRKRkREZmTSr1z9Oabb+Ls2bMlv0aPHl3hcxU/UvPx8YFKpULHjh1Rr1491K9fHz169MDmzZsrdWeKqrfsrAc4/mk3eDy8JBi/ZdMAL3z2C+q5uImUjIiIzE2l7hw5OzvD2dm55OvK7FNW/DL2yZMncfz48ZLxvLw8nDlzBmfOnEFUVBQ2bdoEGxubioemaidfo8Ghz/rAN+tvwbja2gXNZh7Fcw0ai5SMiIjMkWQ+slN850in0+G9997DhAkT4O7ujlu3bmH9+vVYsWIFoqKisGDBAsybN++J5woMDCzX90xMTISHhwfUanWl85O+tLS0Kv8eOp0OvywdhRcyTgnGMyzrwH5UJOwd6/D6GmCKa0MVw2sjXbw20qXVao1640QyC70UFhbC19cXYWFhWLVqFZo1awZbW1t4eHggPDwcs2bNAgCsXr0a6enpTzkbycXRb6bjhdT9grEHFg4oGrYFTZryxX4iInp2krlztGzZsifOT5w4EatWrUJmZiaOHTuGQYMGlXlsdHR0ub5nYGAgLC0t9TakJeOqqj/fHas+QUDiJsFYHhTAiEh07Ny7Sr5ndcP/7UsXr4108dpIj7HXrpPMnaOnsbOzK3lcVvwIjuRr36Zl8P3zC8GYFpbIHLQebVmMiIioEsymHAEoefk7Ly9P5CQkpl9+2oxGR6bqjSd1W4yu/d4RIREREVUnknislpaWhnPnzsHGxgZdunQp87gHDx4AgOATciQvZ48fRO3to2EN4bIOFwJnIPjdySKlIiKi6kQS5SgrKwtBQUEAgJiYGPj4+Ogdo9VqERMTAwAICAgwaT6Shrhzf6Bo42DYQSMYj/UegZDxESKlIiKi6kYSj9U8PT3RqlUrAMDChQsNHrN+/Xqo1Wo0btwYr7zyiinjkQRcT7yE9BWvw0mXLRiPrd8Xg8JUIqUiIqLqyOTlKCAgAAEBAZgzZ45gfNq0aQCAHTt2IDQ0FBcuXIBGo0FSUhLCw8NLNrWdP38+rKysTB2bRJRxNw0Xv3gNLo9tJHuhbnsEzd0OK2v+74GIiIzH5I/VilfCTklJEYwHBQVh3LhxWL16NbZu3YqtW7fq/d45c+YgODjYJDlJGvIe5uLEnF7w1lwXjF9xaIlen++HrUIhTjAiIqq2JPFYrdiiRYsQGRmJnj17wtnZGdbW1lAqlRgwYACOHDlScneJ5KFQW4i9c4PhnfWXYPymwh3t5h2GY00nkZIREVF1ZpGdnV0kdgixFC8CGRcXJ3aUaql4246KLpj2Y0QoWl5cLxi7a1UXDWeegkfT5pXOJ2eVvTZUdXhtpIvXRrqaN28OCwuLci8C/TSSunNEVGzP+nC9YpRjUQM1P9jFYkRERFVKEh/lNwW1Wq23J5tGo4Gtra1IiagsR/f8AM+TnwnGtLBE7pDv0LlNJ5FSERGRXMimHKlUKkRE6K+F4+LiIkIaKkvM6V9Qe9cYWD22yOONV79Avz4hIqUiIiI5kU05Cg0NLVloslhISAjvHEnIlcvn8VAVjLpFwkUez7cciyHv8WV8IiIyDdmUI6VSqfcSnUKhgKUlX7uSgrTUO7i2uDca6DIF47EN3sDgKSvFCUVERLLEZkCiy3uYi9Pze6NBfrJg/GKt1gia8yMLLBERmRR/6pCodDod9nz+Jppl/S0Yv27nha7zDnCRRyIiMjmWIxLVzhXT4H/7J8FYmlU9+M88iNp1nEVKRUREcsZyRKL5ebsKfn8tFYzlWtih9oc70ci9qUipiIhI7liOSBT//f1X1PtpvGBMBwtkDliL51t3FCkVERERyxGJ4EZiAnLWDYIdhB/Zv9Q6DN0HDBcpFRER0SMsR2RSD+5nIm5RbzgXZgjGYxsFIWhsuEipiIiI/h/LEZmMVqvFz3PfQJO8q4LxS04vIejTrfzIPhERSQJ/GpHJ7PhiNPwyTgrGbtk2RKfP9vEj+0REJBmyWSGbG8+Ka++3i+Cf8L1g7L5lTbhP+gnOLq4ipSIiItInm3LEjWfF8/vxA3A/PlswVgBrFL2zCd5+L4gTioiIqAyyKUfceFYct25eg/b7obCBVjCe1DUC/br1FykVERFR2WRTjrjxrOnl5ubg1rp34Vl4TzAe22w4QoZPFSkVERHRk7EZUJXQ6XQ4s+p9eOZdFoxfrP0ygmesFykVERHR07EcUZXYtfoTPJ92WDB2y6YBOs/eA2tr2dywJCIiM8RyREZ3bN82+EQvFozlWNRAg49385NpREQkeSxHZFQXz/8XNSJHwxJFJWM6WOBB0Fr4tQoUMRkREVH5sByR0dzLSEfS8gGoWZQjGL/4wkR06/+uSKmIiIieDcsRGUWhthBH5w9Ew/ybgvF/6nXDwAlLREpFRET07FiOyCh2rpyqtzXINTsvtBn3DZdLICIis8KfWlRpxw/8CN9zywVjmZa18NzoTXB0rClSKiIioophOaJKSbxyAbY/Cl/ALoQlCt/cgIbuXiImIyIiqhiWI6qw3NwcXFgyALV0WYLxSy9ORsdewSKlIiIiqhzZrManVquRnp4uGNNoNNxbrYJ0Oh1+Ch8K/4fCFbDj6nVB8PhFIqUiIiKqPNmUI5VKhYiICL1xFxcXEdKYv6iNX8I/eY9g7JZtQ/SYvZ0vYBMRkVmTTTkKDQ1FUFCQYCwkJIR3jiog+vQRND7+qWAs18IODT7aidp1nEVKRUREZByyKUdKpRJKpVIwplAoeJfjGaXeSUb2hrehhFYwnt57GV57vrVIqYiIiIyHzYDKrVBbiFNfDIJSmyYYj202HK+FfCBSKiIiIuNiOaJy27lyKppnnhWMXa4ZgIHT14mUiIiIyPhYjqhcTh/dg+bnVgjG7lrVRZsZO2HD97aIiKgaYTmip0q5fRPaLSNhBV3JWCEsYTX0W9Rv2ES8YERERFWA5YieqFBbiNNfBKNeYYZg/OILH6Ndt/4ipSIiIqo6LEf0RDtXTEHz+38Kxi7Wao2B45eIlIiIiKhqsRxRmU4d2YPmf68UjKVb1UW7GTtgZW0lUioiIqKqxXJEBt25lQTdv/XfM7Ie+i3c6jcSMRkREVHVYjkiPVqtFr8vGgRnvmdEREQyJJsVsrnxbPntWjkNLR5/z6j2y3zPiIiIZEE25Ygbz5bPHycPweex9YzSrZzRLiyS7xkREZEsyKYccePZp7uXkY7sje/B8fH3jN75ju8ZERGRbMimHHHj2SfT6XQ4vOgdtNSmCMYv+o/F4K79REpFRERkemwGBAA4sGU5WqYcEoxdcfTHgI+/EikRERGROFiOCAkX/4Hy6CzBWJaFA/wn/ch904iISHZYjmROk5eHS8tD4FD0UDB+r/cSeDRtLlIqIiIi8Ui+HOl0Onh5eWHo0KFiR6mW9iwbD8/ci4Kx2Ib98VrIByIlIiIiEpfky9GhQ4eQmpoqdoxq6befd8Hv4gbB2C3bhug9faM4gYiIiCRA0uXo6tWrCAsLEztGtZRxNw2F28bAEkUlYwWwhvOoTXCqVVu8YERERCKT3Ef5z507hy1btiAmJgbR0dHQ6XRP/030zI4sHo6WhcIVw68ETkVw2y7iBCIiIpIIyZWjU6dOYc2aNWLHqNYORX6DlncOCMYuOb2EAR8uECkRERGRdEjusdqbb76Js2fPlvwaPXq02JGqlVs3r8Fp/zTBWJaFA16ctJXbgxAREUGCd46cnZ3h7Oxc8jX3PjMenU6H378aCl9dlmD8bo+FeNmzmUipiIiIpEVy5cgYAgMDy3VcYmIiPDw8oFarqziRNPz640r43zsjGPvHuQu6dBtSJX8GaWlpRj8nGQevjXTx2kgXr410abVa2NjYGO18knusRlXj+pULaBq9WDB216oOng9dyf3liIiISqmWd46io6PLdVxgYCAsLS31NqStbgry83F28zh4FWkE49qBq+Dr27LKv391//M1Z7w20sVrI128NtJjbW3cOsNbBjKwZ/UMeOXECcZi3Qfh1b5vi5SIiIhIuliOqrm4c3+g6bmVgrE7NvXRe4pKpERERETSxnJUjRXk5+PaN8NhC23JWCEsUXP4t1wFm4iIqAwsR9XYnjWz9DaVjfcejjYde4mUiIiISPpYjqqpC7Ex8PrvcsFYsm0j9J24QqRERERE5oHlqBrSarW4smY4FCgoGdPBArXeWQcHB0cRkxEREUkfy1E1tGftp3qfTotr9i5e7vyaSImIiIjMB8tRNXMp/i94/LlUMHbbpgFen8DHaUREROUh+XI0a9YsZGdnY8uWLWJHkTytVotLq4fDDvklYzpYwGHoWjjWdBIxGRERkfmQfDmi8vtJ9TmaZscKxuK83sYrr/YVKREREZH5qZbbhxiiVquRnp4uGNNoNLC1tRUpkXFdT7yERmeFe6fdsXkOfSauFikRERGReZJNOVKpVIiIiNAbd3FxESGNcel0OkSvHA2/ojzBuN2ba1DTqZZIqYiIiMyTbMpRaGgogoKCBGMhISHV4s7R4R3r4Zfxm2AstvFAhHTvL1IiIiIi8yWbcqRUKvV2UlYoFLC0NO/XrjLupsH+wCfCMas66P7xv0RKREREZN7MuxkQDi8fB+fCDMFYds9wOLu4ipSIiIjIvLEcmbHfjx+A/41IwdiFOu3Qc8j7IiUiIiIyfyxHZkqTl4f7Wz4UjD20UCDgo/Vm/6iQiIhITPwpaqb2/msGGmluCMZutJ4Mz6a+IiUiIiKqHliOzNCl+L/Q9G/hC9fXanijX+hccQIRERFVIyxHZkan0yF+7RgoUFAyVghL1B+xDjbVYFkCIiIisbEcmZmfI9fB5360YCy+2TAEvNxZpERERETVC8uRGbmfeQ/2h2YLxtTWLug9/mtxAhEREVVDLEdm5ODqKahXeFcwlt87HE61aosTiIiIqBpiOTITcef+QPPLPwjGLtZ+Gd0HjhIpERERUfUkm+1D1Go10tPTBWMajcYs9lbT6XRIWP8BvFFYMpYPa7R8fy3XNCIiIjIy2ZQjlUqFiIgIvXEXFxcR0jybQ/9ZA++svwRjl/1GYYjfC+IEIiIiqsZkU45CQ0MRFBQkGAsJCZH8naPMe3fhcPgzwViqtRJ9xn4pUiIiIqLqTTblSKlUQqlUCsYUCoXkH0sdWj0F/o9tLKt9PQKONZ1ESkRERFS9SbsZyNzF8/9F84TNgrELddqh24D3xAlEREQkAyxHEha34WPYPPYSdqsP1kj+bhcREZE5409ZiTq2bxt8M04Jxi43fw/NmrcSKREREZE8sBxJUL5Gg/zdYYKxu1Z10etDvoRNRERU1ViOJGjft+FomJ8kGMvs/Alq1a4jUiIiIiL5YDmSGHXKLTx3dplg7FoNH/QZOkmkRERERPLCciQxv66ZAqeibMGY81vLYGVtJVIiIiIieWE5kpDYv87A73qkYOy8a0+07dxbpERERETyw3IkIQkbJ8EKupKv86BAmw9WiJiIiIhIfmSzQrbUN549fuBHNM88Kxi72nI0XvT0ESkRERGRPMmmHEl541mtVouc3TNROkm6lTNeez9ctExERERyJZtyJOWNZw9tXQmPvKuCsfudwlDTqZZIiYiIiORLNuVIqhvP5uRkw+GY8I5WkqIJer/Dj+4TERGJgS9ki+zA+nlQatMEY4o3FsDaWja9lYiISFJYjkSkTrmFxn/9SzB2yekldO7zlkiJiIiIiOVIRL+umwHHolzBmPvQJaI/6iMiIpIz/hQWyZWLsWh+dZtg7Lzba3ixbRdxAhEREREAliPR/PfbqbCBtuTrfFgjMHSpiImIiIgIYDkSxd8xv6Fl6s+CsUvNhsKzqa9IiYiIiKgYy5EIEjbPFHydZeGAbu9/IVIaIiIiKo3lyMT+/O1n+GWcFIwltRqNei5uIiUiIiKi0liOTCz5P8K7RpmWtdBz5ByR0hAREdHjWI5M6Lefd8HnQYxgLCVwHGrVriNSIiIiInqcbJZhVqvVSE9PF4xpNBqT7a2m0+lwb/enqFtqLM2qHl4bMbPM30NERESmJ5typFKpEBERoTfu4uJiku9/LGoLvHLiBGOZ7SfB3t7BJN+fiIiIykc25Sg0NBRBQUGCsZCQEJPcOSrUFkKzb55g7I7Nc+j9zpQq/95ERET0bGRTjpRKJZRKpWBMoVCYZKuOw9vXwT3vqmBM8+oM2CoUVf69iYiI6NnwhewqptVqYXlUuIbRTYU7eoaMFSkRERERPUmly1FRURFUKhU6duwINzc3NGrUCH369MH+/fuNkc/sHdmuQsP8JMGYRa9PYW0tm5t2REREZqVS5aioqAhDhw7FpEmT8NdffyE7Oxv37t3DiRMnMGTIEHzxxbOt+vzZZ5/B0dHxib+GDx9emcgmVagtBH75UjB23c4L3YNGiJSIiIiInqZS5Wj16tXYu3cvFAoFli9fjtu3b+Pq1auYOHEiACA8PBwnT5588klKSUhIqEwcyTm6+1s01lwXjFn3/MQk7zkRERFRxVT42Y5Go8HSpY92kQ8PD8eoUaMAAE5OTliwYAHS09OxefNmLF68GB07dizXOa9cuQIAOHToENq3b1/RaJJQqC1E4eFFgrEkhQd69H9PnEBERERULhW+hXHq1Cmo1WrUrVsXI0eO1JufPHkyAOD48ePIzMx86vl0Oh2uXbsGAPD1Nf/d6X/9aZPeJ9QseoTBytpKpERERERUHhUuRydOnAAAdOrUyeBaQd7e3nB3d0dhYSFOnz791PPdvHkTeXl5cHNzQ926dZ96vJTpdDpoDgkXnLypcEf3oNEiJSIiIqLyqnA5Kn4/yN/fv8xjiufK8y5R8SM1Hx+fkk+/1atXD/Xr10ePHj2wefNm6HS6isY1qWP7/g2Ph5cFY7pXp/GuERERkRmo8DtHycnJAIAGDRqUeUz9+vUBAElJSWUeU6y4QJ08eRLHjx8vGc/Ly8OZM2dw5swZREVFYdOmTbCxsalo7Cqn0+mQs3+hYOyWbUP0GPy+SImIiIjoWVS4HGVnZwMAHB0dyzymeC4nJ+ep5yu+c6TT6fDee+9hwoQJcHd3x61bt7B+/XqsWLECUVFRWLBgAebNm/fEcwUGBpbrvyExMREeHh5Qq9XlOr48/ji2F165FwRj99uMRUZGhtG+h7lIS0sTOwKVgddGunhtpIvXRrq0Wq1Rb5xU+LFafn4+ADxxb7LioLm5uU89X2FhIXx9fREWFoZVq1ahWbNmsLW1hYeHB8LDwzFr1iwAj5YPSE9Pr2jsKqf9dbng69s2z+GVPu+KlIaIiIieVYXvHBWXIo1GU+YxxXOKcuwhtmzZsifOT5w4EatWrUJmZiaOHTuGQYMGlXlsdHT0U78f8OgOk6Wlpd6eaxX1+7F98M49Lxh72HHyEx89yoGx/nzJ+HhtpIvXRrp4baTH2LtOVPjOUfEjs+LHa4ZkZWUBABwcHCr6bUrY2dmVPC4rfgQnNcl7hCuCq61d0OvN8SKlISIiooqocDlq2LAhAODWrVtlHnPnzh3BsZXl7OwM4NFL2lJz/q/f4Zfxm2As46X3YVuOu2ZEREQkHRW+D+Xt7Q0AiI2NLfOY8+fPC44tS1paGs6dOwcbGxt06dKlzOMePHgA4P9LkpTE/RiO0osa3LesiR7vTBUtDxEREVVMhe8cderUCcCjxSCLX84u7fLly0hKSoKVlRU6dOjwxHNlZWUhKCgIffv2xaVLlwweo9VqERMTAwAICAioaOwqcT3xEprfPiAYu+k3DDWdaomUiIiIiCqqwuWoffv2cHV1xb1797Bx40a9+RUrVgAAXn311aeueO3p6YlWrVoBABYuXGjwmPXr10OtVqNx48Z45ZVXKhq7SpzdvAA2KCz5Og8KdBk2U8REREREVFEVLke2trYl+6fNnDkTP/zwA7KyspCamoo5c+Zg48aNsLS0xIwZMwS/LyAgAAEBAZgzZ45gfNq0aQCAHTt2IDQ0FBcuXIBGo0FSUhLCw8MRFhYGAJg/fz6srKSz0nRa6h00TYwUjCV4DITSTd6fUCMiIjJXlfrs29ixY3HmzBns3r0bY8eOxdixYwXzn3/+Odq2bSsYK14JOyUlRTAeFBSEcePGYfXq1di6dSu2bt2q9/3mzJmD4ODgykQ2umObI9Ci6P+XM9DCEq3f+VTERERERFQZFb5zBAAWFhbYtGkTvv76awQEBMDBwQG1a9dG586dsWPHDkycOPGZzrdo0SJERkaiZ8+ecHZ2hrW1NZRKJQYMGIAjR46U3F2SiuysB2hw/gfB2EW3nvBs6itSIiIiIqqsSq+aZGFhgdGjR2P06PLtOP+kdZEAoHfv3ujdu3dlY5nE4X8vg4/uvmCs+RDeNSIiIjJnlbpzJGdarRa1/lwrGLtQ5xW0eqmdSImIiIjIGIy73raEqdVqvT3ZNBrNE/eGe5Jf9/6ABgXC96ae6xdW4XxEREQkDbIpRyqVChEREXrjLi4uFTpf7i8rBF9fq+GD3q/2q9C5iIiISDpkU45CQ0MRFBQkGAsJCanQnaP//n4MzbL+FoxZdBwLS0s+pSQiIjJ3silHSqVSbydlhUJRoUJzZfditCz1dbpVXXQLHlPJhERERCQFvNXxjJKTEuFz55BgTO0/HAo7O5ESERERkTGxHD2j01u/fGyrEFt0emuKiImIiIjImFiOnkFOTjYaXhSu3J3QuB+3CiEiIqpGWI6ewdH/rEJt3QPBmP8gfnyfiIioOmE5KiedTgfF78JFHy/Wfhl+z7cWKRERERFVBZajcjp5MBKNNDcEY869JomUhoiIiKoKy1E5pf+8XPD1TYU7Or42WKQ0REREVFVYjsoh4eI/8L13RjCW9/L7XPSRiIioGuJP93L47/algq8fWDqi25BxIqUhIiKiqsRy9BTZWQ/gfmWHYOyGVzAcazqJlIiIiIiqkmy2D1Gr1UhPTxeMaTSap+6t9uuOb+BVlC0Ye3HQZKPnIyIiImmQTTlSqVSIiIjQG3dxcSnz9+h0Olj+vl4wdqFOWwQ3b2X0fERERCQNsilHoaGhCAoKEoyFhIQ88c5RzOkj8Hh4WTBWp+vYKslHRERE0iCbcqRUKqFUKgVjCoXiiZ84uxa1HC1LfZ1i7YZOfd6qooREREQkBXwhuwypd5LhfednwVhGq3dgbS2bPklERCRLLEdlOPmfr2ELbcnXGtigU8hE8QIRERGRSbAcGaDVauEcu1kwdrl+LyjdGoiUiIiIiEyF5ciAY1Fb4KpNFYx5vv6xSGmIiIjIlFiODHjw6xrB19dq+OCldt1ESkNERESmxHL0mMvx59A886xgTPfKaO6jRkREJBP8if+Yv3YsE3z9wNIRXQd9IFIaIiIiMjWWo1Kysx7A/epOwdgNr2A4ODiKlIiIiIhMjeWolF+3r4VTqX3UdLDAS4OmiJiIiIiITE02Kxo+beNZnU4Hq7PCfdQu1WmL4Ob+JstIRERE4pNNOXraxrPRpw6jycMEwVydrh+aJBsRERFJh2zK0dM2nr2+j/uoERERkYzK0ZM2nk25fRPedw4L5u49P4z7qBEREckQX8gG8Nt/lunto9ZxyAQRExEREZFYZF+OioqK4BK7STB2uf5r3EeNiIhIpmRfjnKy7sOlUPgptmb9J4mUhoiIiMQm+3JUlJsh+PqqvR9ebPuqSGmIiIhIbLIvR7aFeYKvLduHipSEiIiIpED25ai0e5a18OpAliMiIiI5Yzkq5ZbPENjbO4gdg4iIiETEcvQ/hbBEmyHcR42IiEjuWI7+52K9Tmji6SN2DCIiIhIZy9H/uPX8SOwIREREJAGy2R9DrVYjPV24npFGo4ElgGTbxujeI8jwbyQiIiJZkU05UqlUiIiI0Bt3reOIgs6TYGnJm2hEREQko3IUGhqKoCDh3aGQkBDY2tri9XcmihOKiIiIJEc25UipVEKpVArGFAoF7xgRERGRAJsBERERUSmSLEdFRUVQqVTo2LEj3Nzc0KhRI/Tp0wf79+8XOxoRERFVc5J7rFZUVIShQ4di7969gvETJ07gxIkTmD17NmbMmCFSOiIiIqruJHfnaPXq1di7dy8UCgWWL1+O27dv4+rVq5g4cSIAIDw8HCdPnhQ3JBEREVVbkipHGo0GS5cuBfCoBI0aNQpOTk5wdXXFggUL8M4776CoqAiLFy8WOSkRERFVV5IqR6dOnYJarUbdunUxcuRIvfnJkycDAI4fP47MzEwTpyMiIiI5kFQ5OnHiBACgU6dOsLW11Zv39vaGu7s7CgsLcfr0aVPHIyIiIhmQVDlKSEgAAPj7+5d5TPFc8bFERERExiSpcpScnAwAaNCgQZnH1K9fHwCQlJRkkkxEREQkL5L6KH92djYAwNHRscxjiudycnLKPCYwMLBc3y8xMREeHh5Qq9XPkJLKKy0tTewIVAZeG+nitZEuXhvp0mq1sLGxMdr5JHXnKD8/HwAMvm9UrPg/Pjc31ySZiIiISF4kdeeouBRpNJoyjymeUygUZR4THR1dru8XGBgIS0tLvT3XyLj45ytdvDbSxWsjXbw20mNtbdw6I6k7R8WPzIofrxmSlZUFAHBwcDBJJiIiIpIXSZWjhg0bAgBu3bpV5jF37twRHEtERERkTJIqR97e3gCA2NjYMo85f/684FgiIiIiY5JUOerUqROAR4tBFr+cXdrly5eRlJQEKysrdOjQwdTxiIiISAYkVY7at28PV1dX3Lt3Dxs3btSbX7FiBQDg1VdfRd26dU2cjoiIiORAcp9Wmzx5MsLCwjBz5kzY2dkhKCgIubm5+Ne//oWNGzfC0tISM2bMMMr3S05ORkFBAVq0aGGU85FQYWEhAMDKykrkJPQ4Xhvp4rWRLl4b6bp27doTlwF6VhbZ2dlFRjubERQVFWHYsGHYvXu3wfkFCxZg4sSJRvlenp6eyM3NFbzcXVhYiHv37qFOnTpG+38AczhnVWRMTEwE8OjP2RjM4c+xKs5pDtcGMI//bnPIyGsj3Yy8NtI95+XLl2FhYYH79+8bIZ0EyxHwqCBt2LAB33//PS5fvgwbGxs8//zz+Pjjj9GrV68q/d7x8fFo06YN/vjjD/j5+cnmnFWRsXil8vKuO/U05vDnWBXnNIdrA5jHf7c5ZOS1kW5GXhvpntPY10ZSj9WKWVhYYPTo0Rg9erTYUYiIiEhmJPVCNhEREZHYWI6IiIiISmE5IiIiIiqF5egx9erVwyeffIJ69erJ6pxVkdHYzOHPsSrOaQ7XBjCP/25zyFgVzOG/2xwyVgVz+O82l3MakyQ/rUbVQ1V8soOMg9dGunhtpIvXRrqMfW1454iIiIioFJYjIiIiolJYjoiIiIhK4TtHRERERKXwzhERERFRKSxHRERERKWwHBERERGVwnJEREREVArLEREREVEpLEdkVEVFRVCpVOjYsSPc3NzQqFEj9OnTB/v37xc7GhFRhcTExGD48OHw9vZGnTp10KBBA/To0QMbNmxAYWGh2PGoCvCj/GQ0RUVFGDp0KPbu3Wtwfvbs2ZgxY4aJU1FZdDodmjVrhrZt22LLli1ix5G91NRUrFy5EgcPHkRSUhIAoEmTJujXrx/Gjx+P2rVrixtQpiIjIxEaGgqtVmtwvnv37oiMjISNjY2Jk1FZPv30UyxbtgxLlizBBx98UKFz8M4RGc3q1auxd+9eKBQKLF++HLdv38bVq1cxceJEAEB4eDhOnjwpbkgqcejQIaSmpoodgwDEx8ejXbt2+Prrr3Hx4kXk5uYiNzcX8fHxWLRoEdq3b48rV66IHVN20tPTMX78eGi1WrRu3Ro///wzbt++jQsXLmDevHmwtbXFkSNH8OWXX4odlf7n9OnTWL58eaXPw3JERqHRaLB06VIAj0rQqFGj4OTkBFdXVyxYsADvvPMOioqKsHjxYpGTEgBcvXoVYWFhYscgPLrjOmLECKSmpsLLyws7d+6EWq3G1atXoVKpoFQqcePGDQwZMgQFBQVix5WV7du3Izs7G/Xr18e+ffvQrl07ODk5oVGjRpgyZQqmTZsGANi4caO4QQkAkJWVhTFjxkCn01X6XCxHZBSnTp2CWq1G3bp1MXLkSL35yZMnAwCOHz+OzMxME6cjADh37hymTZuGrl27IiAgAImJiWJHIgC//PIL4uLiYGNjg127dqFnz56wt7eHq6sr3nrrLRw5cgT29va4fPkydu/eLXZcWYmJiQEA9OvXD/b29nrzAwcOBADcuXMHGRkZJs1G+sLCwnD9+nWjnIvliIzixIkTAIBOnTrB1tZWb97b2xvu7u4oLCzE6dOnTR2P8KjArlmzBn/88YdR/mVFxvHrr78CALp06QJPT0+9eU9PTwQFBQF4dA3JdNRqNYBH734ZUrNmzZL/u6iIr++Kad++ffjhhx/QsmVLvPzyy5U+H8sRGUVCQgIAwN/fv8xjiueKjyXTevPNN3H27NmSX6NHjxY7EuHRI04A8PX1LfMYpVIJAMjJyTFJJnpkz549yM7Oxvjx4w3OF5dVNzc3ODs7mzIalZKWloaPPvoItra2WL9+vVFejrc2Qi4iJCcnAwAaNGhQ5jH169cHgJJP4pBpOTs7C/4Cd3FxETENFfvwww8xcOBA+Pn5lXnMuXPnAADu7u4mSkVlycvLQ0pKCk6cOIHZs2cDQMm7RySO8ePHIy0tDZ9//jlatmxplHOyHJFRZGdnAwAcHR3LPKZ4jv/6Jfp/nTp1euL8kSNHSh699e3b1xSRqAxr167F1KlTS762t7fH0qVLMWbMGBFTydsPP/yAqKgovPLKKyWfjDYGPlYjo8jPzwcAg+8bFSu+1Zmbm2uSTETmbsuWLRg6dCgAIDg4GC+88IK4gUggNzcXhw8fxp07d8SOIks3btxAWFgYHBwc8M0338DS0niVhuWIjKK4FGk0mjKPKZ5TKBQmyURkrmJjY9GnTx+8//77yMnJQYcOHbB27VqxY8neBx98gOzsbFy5cgVbtmxBs2bNcODAAfTp0+eJf/eR8el0OoSGhiIrKwsLFy40+GGGymA5IqMofmRW/HjNkKysLACAg4ODSTIRmZv79+9j0qRJaN++PU6cOAEbGxvMmjULUVFRqFGjhtjx6H/c3NzQv39/REVFoVatWkhISMCuXbvEjiUry5cvx+nTp9GjRw+MGjXK6OdnOSKjaNiwIQDg1q1bZR5TfOu5+Fgi+n/R0dFo27YtVCoVdDod+vXrh+joaHzyySewtubroVLUoEGDknfGYmNjRU4jH9euXcP8+fNRt25drFmzpkq+B8sRGYW3tzeAJ/8Fcf78ecGxRPTI8ePH0adPH9y8eRPu7u7Yv38/tm7dCi8vL7Gjydbt27fh6uoKV1fXkk/jGlL8qc/iO+NU9W7evIn8/HxkZGSgadOmcHR0FPz67bffAABTp04tGXvWxYdZjsgoiv/1dOLEiZKXs0u7fPkykpKSYGVlhQ4dOpg6HpFk3b17F8OGDUNubi46d+6MM2fOPPUTbFT1XF1dUVRUhJycnCfua1e8btuTljEh88NyREbRvn17uLq64t69ewb3GVqxYgUA4NVXX0XdunVNnI5IutatW4eMjAy4u7sjMjISTk5OYkciAFZWViUl9fvvvzd4zJ9//lmyEGS3bt1Mlk3uOnXqhOzs7DJ/Ff8DfMmSJSVjtWvXfqbvwXJERmFra1uyf9rMmTPxww8/ICsrC6mpqZgzZw42btwIS0tLzJgxQ+SkRNKyb98+AMDIkSMN7t9F4hk7diwAIDIyEiNGjMDff/+NnJwc3Lx5E9999x2Cg4Oh0+nw2muvITAwUOS0ZEwsR2Q0Y8eOxYABA5CXl4exY8fiueeeg5eXF7766isAwOeff462bduKnJJIOgoKChAXFwcAmDNnjt67E4//mj59usiJ5aVr164liz5GRkaW3CH39fXF+PHjkZGRgTZt2kClUomclIyN5YiMxsLCAps2bcLXX3+NgIAAODg4oHbt2ujcuTN27Nhh1NVLiaoDtVqNgoICsWPQE8ydOxe7d+9Gnz59oFQqYW1tjdq1a6Ndu3ZYvnw5Dh06hDp16ogdk4zMIjs7m1sJExEREf0P7xwRERERlcJyRERERFQKyxERERFRKSxHRERERKWwHBERERGVwnJEREREVArLEREREVEpLEdEREREpbAcEREREZXCckRERERUCssRERERUSksR0RERESlsBwRERERlcJyRERERFQKyxERERFRKSxHRERERKX8H2OT5QZhZUi0AAAAAElFTkSuQmCC",
+ "text/plain": [
+ "Figure(nrows=1, ncols=1, refwidth=2.5)"
+ ]
+ },
+ "metadata": {
+ "image/png": {
+ "height": 282,
+ "width": 291
+ }
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "exec(code)\n",
+ "\n",
+ "preds = _spergel_hlr(nu_vals)\n",
+ "diff = np.abs(frac_vals - preds)\n",
+ "print(np.mean(diff), np.max(diff))\n",
+ "\n",
+ "preds = ratnl_func(nu_vals, *coeff)\n",
+ "diff = np.abs(frac_vals - preds)\n",
+ "print(np.mean(diff), np.max(diff))\n",
+ "\n",
+ "\n",
+ "fig, axs = uplot.subplots()\n",
+ "\n",
+ "# axs.plot(nu_vals, hlr_vals)\n",
+ "axs.plot(nu_vals, frac_vals)\n",
+ "axs.plot(nu_vals, preds)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "a929b493-9997-4c04-b663-edbe2483e187",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "9.079122786914695e-09 6.302919142164853e-08 -0.8412399010431769\n"
+ ]
+ }
+ ],
+ "source": [
+ "rng = np.random.RandomState()\n",
+ "\n",
+ "diffs = []\n",
+ "nus = []\n",
+ "for _ in range(10000):\n",
+ " nu = rng.uniform(low=-0.85, high=4.0)\n",
+ " re = 1.0\n",
+ "\n",
+ " try:\n",
+ " ival = _spergel_hlr(nu)\n",
+ " tval = galsim.Spergel(nu, scale_radius=re).calculateFluxRadius(NU_FRAC)\n",
+ " except Exception:\n",
+ " pass\n",
+ " else:\n",
+ " diffs.append(ival - tval)\n",
+ " nus.append(nu)\n",
+ "\n",
+ "print(np.median(np.abs(diffs)), np.max(np.abs(diffs)), nus[np.argmax(np.abs(diffs))])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7bbac0ff-e7bd-4a99-bcf6-2478a69434da",
+ "metadata": {},
+ "source": [
+ "# CODE NOT USED BELOW HERE"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9f5a7d62-b320-469f-a875-2c8a796b94c3",
+ "metadata": {},
+ "source": [
+ "## Try rad at fixed nu"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "51614604-c11d-4cbb-bb28-c10a5fe3c37d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "NU = 0.5\n",
+ "\n",
+ "n_pts = 500\n",
+ "n_log = 200\n",
+ "n_mid = n_pts - 2 * n_log\n",
+ "frac_cut = 0.4\n",
+ "frac_min = 1e-8\n",
+ "eps = 1e-12\n",
+ "logfracs = np.logspace(np.log10(frac_min), np.log10(frac_cut), n_log)\n",
+ "frac_vals = np.concatenate(\n",
+ " [\n",
+ " logfracs,\n",
+ " np.linspace(frac_cut + eps, 0.5 + np.abs(0.5 - frac_cut) - 2 * eps, n_mid),\n",
+ " 1.0 - logfracs[::-1] - eps,\n",
+ " ]\n",
+ ")\n",
+ "\n",
+ "prof = galsim.Spergel(NU, scale_radius=1)\n",
+ "rad_vals = []\n",
+ "for frac in frac_vals:\n",
+ " rad_vals.append(prof.calculateFluxRadius(frac))\n",
+ "\n",
+ "rad_vals = np.array(rad_vals)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "db923966-1218-4978-9472-b662a05d52b7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, axs = uplot.subplots()\n",
+ "\n",
+ "axs.plot(np.log(frac_vals), np.log(rad_vals))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0022f769-4787-4f86-ad6e-e701a836e29d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import scipy.optimize\n",
+ "import numpy.polynomial\n",
+ "from numpy.polynomial import Polynomial\n",
+ "\n",
+ "numpy.polynomial.set_default_printstyle(\"ascii\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2eb066bf-ffdd-45a6-a17f-541ff744d71f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RATNL_ORDER = 7\n",
+ "\n",
+ "\n",
+ "def get_ratnl_func_polys(coeff):\n",
+ " p_coeff = coeff[0 : RATNL_ORDER + 1]\n",
+ " q_coeff = np.concatenate([[1], coeff[RATNL_ORDER + 1 :]])\n",
+ " pm = Polynomial(p_coeff)\n",
+ " qm = Polynomial(q_coeff)\n",
+ " return pm, qm\n",
+ "\n",
+ "\n",
+ "def ratnl_func(x, *coeff):\n",
+ " pm, qm = get_ratnl_func_polys(coeff)\n",
+ " return pm(x) / qm(x)\n",
+ "\n",
+ "\n",
+ "res = scipy.optimize.curve_fit(\n",
+ " ratnl_func,\n",
+ " np.log(frac_vals),\n",
+ " np.log(rad_vals),\n",
+ " p0=np.ones(2 * RATNL_ORDER + 1),\n",
+ " full_output=True,\n",
+ " maxfev=100000,\n",
+ ")\n",
+ "print(res[-2])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "99e31d6c-4c99-41da-9d84-068a946c471b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "coeff = res[0]\n",
+ "\n",
+ "pm, qm = get_ratnl_func_polys(coeff)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d49e0d06-58d9-49fe-9667-030aeb588b19",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rng = np.random.RandomState()\n",
+ "\n",
+ "diffs = []\n",
+ "fracs = []\n",
+ "for _ in range(10000):\n",
+ " frac = rng.uniform(low=frac_min, high=1.0 - frac_min - eps)\n",
+ " re = 1.0\n",
+ "\n",
+ " try:\n",
+ " ival = ratnl_func(np.log(frac), *coeff)\n",
+ " tval = np.log(galsim.Spergel(NU, scale_radius=re).calculateFluxRadius(frac))\n",
+ " except Exception:\n",
+ " pass\n",
+ " else:\n",
+ " diffs.append(ival - tval)\n",
+ " fracs.append(frac)\n",
+ "\n",
+ "print(\n",
+ " np.median(np.abs(diffs)), np.max(np.abs(diffs)), (fracs[np.argmax(np.abs(diffs))])\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "de92808e-05ba-4bc1-b736-b3437956e4d7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, axs = uplot.subplots()\n",
+ "\n",
+ "axs.plot(np.log(frac_vals), ratnl_func(np.log(frac_vals), *coeff) - np.log(rad_vals))\n",
+ "# axs.plot(np.log(frac_vals), )"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "71967d9d-4e2c-4edf-9a9f-07e8a61bd3a1",
+ "metadata": {},
+ "source": [
+ "## Try Full Flux Radius Computation w/ Polynomial Product"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4baaabaf-ddcf-459d-8d3b-5c9f397ff476",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RATNL_ORDER = 5\n",
+ "n_params = 2 * (RATNL_ORDER + 1) + 2 * RATNL_ORDER\n",
+ "\n",
+ "\n",
+ "def get_ratnl_func_polys(coeff):\n",
+ " n_p = RATNL_ORDER + 1\n",
+ " n_q = RATNL_ORDER\n",
+ " p_coeff = coeff[0:n_p]\n",
+ " q_coeff = np.concatenate([[1], coeff[n_p : n_p + n_q]])\n",
+ " pm = Polynomial(p_coeff)\n",
+ " qm = Polynomial(q_coeff)\n",
+ "\n",
+ " n_p = RATNL_ORDER + 1\n",
+ " n_q = RATNL_ORDER\n",
+ " p_coeff = coeff[n_p + n_q : n_p + n_q + n_p]\n",
+ " q_coeff = np.concatenate([[1], coeff[n_p + n_q + n_p :]])\n",
+ " pm2 = Polynomial(p_coeff)\n",
+ " qm2 = Polynomial(q_coeff)\n",
+ "\n",
+ " return pm, qm, pm2, qm2\n",
+ "\n",
+ "\n",
+ "def ratnl_func(x, *coeff):\n",
+ " pm, qm, pm2, qm2 = get_ratnl_func_polys(coeff)\n",
+ " return pm(x[0, :]) / qm(x[0, :]) * pm2(x[1, :]) / qm2(x[1, :])\n",
+ "\n",
+ "\n",
+ "n_pts = 500\n",
+ "n_log = 200\n",
+ "n_mid = n_pts - 2 * n_log\n",
+ "frac_cut = 0.4\n",
+ "frac_min = 1e-8\n",
+ "eps = 1e-12\n",
+ "logfracs = np.logspace(np.log10(frac_min), np.log10(frac_cut), n_log)\n",
+ "fracs = np.concatenate(\n",
+ " [\n",
+ " logfracs,\n",
+ " np.linspace(frac_cut + eps, 0.5 + np.abs(0.5 - frac_cut) - 2 * eps, n_mid),\n",
+ " 1.0 - logfracs[::-1] - eps,\n",
+ " ]\n",
+ ")\n",
+ "nu_vals = np.linspace(-0.0, 0.01, n_pts)\n",
+ "\n",
+ "rad_vals = []\n",
+ "x_vals = []\n",
+ "for nu in nu_vals:\n",
+ " for frac in fracs:\n",
+ " try:\n",
+ " rad_vals.append(\n",
+ " galsim.Spergel(nu, scale_radius=1).calculateFluxRadius(frac)\n",
+ " )\n",
+ " except Exception:\n",
+ " rad_vals.append(np.nan)\n",
+ " x_vals.append([nu, np.log(frac)])\n",
+ "rad_vals = np.array(rad_vals)\n",
+ "x_vals = np.array(x_vals).T\n",
+ "\n",
+ "msk = np.isfinite(rad_vals)\n",
+ "rad_vals = np.log(rad_vals[msk])\n",
+ "x_vals = x_vals[:, msk]\n",
+ "\n",
+ "print(x_vals.shape, x_vals[0, 0:10], x_vals[1, 0:10])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7fb46e58-d6ff-4fc3-9eed-5d40b51ade7a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "eps = 1e-6\n",
+ "res = scipy.optimize.curve_fit(\n",
+ " ratnl_func,\n",
+ " x_vals,\n",
+ " rad_vals,\n",
+ " p0=np.ones(n_params),\n",
+ " full_output=True,\n",
+ " maxfev=1000000,\n",
+ " ftol=eps,\n",
+ " xtol=eps,\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cf64003d-3928-4cb9-b712-fbd009cbaa20",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "coeff = res[0]\n",
+ "\n",
+ "pred = ratnl_func(x_vals, *coeff)\n",
+ "diff = np.abs(rad_vals - pred)\n",
+ "\n",
+ "print(np.median(diff), np.max(diff))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "249e67c8-4385-46f6-a089-bf3cb95e3ba5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nu_vals = np.linspace(-0.85, 4.0, 500)\n",
+ "NU_FRAC = 0.005\n",
+ "\n",
+ "frac_vals = []\n",
+ "ifrac_vals = []\n",
+ "for nu in nu_vals:\n",
+ " prof = galsim.Spergel(nu, scale_radius=1)\n",
+ " frac_vals.append(prof.calculateFluxRadius(NU_FRAC))\n",
+ " xvals = np.array([[nu, np.log(NU_FRAC)]]).T\n",
+ " ifrac_vals.append(np.exp(ratnl_func(xvals, *coeff)))\n",
+ "\n",
+ "frac_vals = np.array(frac_vals)\n",
+ "ifrac_vals = np.array(ifrac_vals)\n",
+ "\n",
+ "fig, axs = uplot.subplots()\n",
+ "\n",
+ "axs.plot(nu_vals, frac_vals)\n",
+ "axs.plot(nu_vals, ifrac_vals)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "08d5aade-ab62-40f1-ab87-1026c7e06f6d",
+ "metadata": {},
+ "source": [
+ "## Try full Flux Radius computation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7987d13c-ddd2-40e2-9cc9-96bef3ecaf01",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from numba import njit\n",
+ "\n",
+ "NU_FRAC_ORDER = 3\n",
+ "NU_FRAC_ORDER_PLUS_ONE = NU_FRAC_ORDER + 1\n",
+ "n_params = 2 * (NU_FRAC_ORDER_PLUS_ONE * (NU_FRAC_ORDER_PLUS_ONE + 1) // 2) - 1\n",
+ "\n",
+ "\n",
+ "@njit\n",
+ "def compute_poly_matrix_from_coeff(coeff):\n",
+ " n_per = NU_FRAC_ORDER_PLUS_ONE * (NU_FRAC_ORDER_PLUS_ONE + 1) // 2\n",
+ " p_coeff = coeff[:n_per]\n",
+ " q_coeff = np.zeros(n_per) # p_coeff.copy()\n",
+ " q_coeff[0] = 1.0\n",
+ " q_coeff[1:] = coeff[n_per:]\n",
+ "\n",
+ " loc = 0\n",
+ " c_p = np.zeros((NU_FRAC_ORDER_PLUS_ONE, NU_FRAC_ORDER_PLUS_ONE))\n",
+ " c_q = np.zeros((NU_FRAC_ORDER_PLUS_ONE, NU_FRAC_ORDER_PLUS_ONE))\n",
+ " for i in range(NU_FRAC_ORDER_PLUS_ONE):\n",
+ " for j in range(NU_FRAC_ORDER_PLUS_ONE):\n",
+ " if i + j <= NU_FRAC_ORDER:\n",
+ " c_p[i, j] = p_coeff[loc]\n",
+ " c_q[i, j] = q_coeff[loc]\n",
+ " loc += 1\n",
+ " return c_p, c_q\n",
+ "\n",
+ "\n",
+ "@njit\n",
+ "def ratnl_func(nu_frac, *coeff):\n",
+ " nu, frac = nu_frac[0, :], nu_frac[1, :]\n",
+ " c_p, c_q = compute_poly_matrix_from_coeff(np.array(coeff))\n",
+ " return numpy.polynomial.polynomial.polyval2d(\n",
+ " nu, frac, c_p\n",
+ " ) / numpy.polynomial.polynomial.polyval2d(nu, frac, c_q)\n",
+ "\n",
+ "\n",
+ "n_pts = 500\n",
+ "n_log = 150\n",
+ "n_mid = n_pts - 2 * n_log\n",
+ "frac_cut = 0.4\n",
+ "frac_min = 1e-8\n",
+ "eps = 1e-12\n",
+ "logfracs = np.logspace(np.log10(frac_min), np.log10(frac_cut), n_log)\n",
+ "fracs = np.concatenate(\n",
+ " [\n",
+ " logfracs,\n",
+ " np.linspace(frac_cut + eps, 0.5 + np.abs(0.5 - frac_cut) - 2 * eps, n_mid),\n",
+ " 1.0 - logfracs[::-1] - eps,\n",
+ " ]\n",
+ ")\n",
+ "nu_vals = np.linspace(-0.85, 4.0, n_pts)\n",
+ "\n",
+ "rad_vals = []\n",
+ "x_vals = []\n",
+ "for nu in nu_vals:\n",
+ " for frac in fracs:\n",
+ " try:\n",
+ " rad_vals.append(\n",
+ " galsim.Spergel(nu, scale_radius=1).calculateFluxRadius(frac)\n",
+ " )\n",
+ " except Exception:\n",
+ " rad_vals.append(np.nan)\n",
+ " x_vals.append([nu, np.log(frac)])\n",
+ "rad_vals = np.array(rad_vals)\n",
+ "x_vals = np.array(x_vals).T\n",
+ "\n",
+ "\n",
+ "msk = np.isfinite(rad_vals)\n",
+ "rad_vals = rad_vals[msk]\n",
+ "x_vals = x_vals[:, msk]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "049901f2-3cda-46ab-8733-ac9faafa4cbe",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tol = 1e-2\n",
+ "res = scipy.optimize.curve_fit(\n",
+ " ratnl_func,\n",
+ " x_vals,\n",
+ " np.log(rad_vals),\n",
+ " p0=np.ones(n_params),\n",
+ " full_output=True,\n",
+ " xtol=tol,\n",
+ " ftol=tol,\n",
+ " maxfev=100000,\n",
+ ")\n",
+ "res[-2]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "05c74099-8c84-43bc-8a45-3a05e11c2741",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "coeff = res[0]\n",
+ "\n",
+ "pred = ratnl_func(x_vals, *coeff)\n",
+ "diff = np.abs(np.log(rad_vals) - pred)\n",
+ "\n",
+ "print(np.median(diff), np.max(diff))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "40cc1dcc-3a46-49e5-be97-eb1056d2e489",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nu_vals = np.linspace(-0.85, 4.0, 500)\n",
+ "NU_FRAC = 0.9\n",
+ "\n",
+ "frac_vals = []\n",
+ "ifrac_vals = []\n",
+ "for nu in nu_vals:\n",
+ " prof = galsim.Spergel(nu, scale_radius=1)\n",
+ " frac_vals.append(prof.calculateFluxRadius(NU_FRAC))\n",
+ " xvals = np.array([[nu, np.log(NU_FRAC)]]).T\n",
+ " ifrac_vals.append(np.exp(ratnl_func(xvals, *coeff)))\n",
+ "\n",
+ "frac_vals = np.array(frac_vals)\n",
+ "ifrac_vals = np.array(ifrac_vals)\n",
+ "\n",
+ "fig, axs = uplot.subplots()\n",
+ "\n",
+ "axs.plot(nu_vals, frac_vals)\n",
+ "axs.plot(nu_vals, ifrac_vals)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "856a6679-15ba-46b7-b8fa-c4914e5fb5e0",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0c0756b9-e21c-49ba-81eb-c48fd14d0edd",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python [conda env:jax-galsim]",
+ "language": "python",
+ "name": "conda-env-jax-galsim-py"
+ },
+ "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.10.15"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/jax_galsim/spergel.py b/jax_galsim/spergel.py
index a5a7389b..ba6493f5 100644
--- a/jax_galsim/spergel.py
+++ b/jax_galsim/spergel.py
@@ -190,6 +190,54 @@ def calculateFluxRadius(alpha, nu, zmin=0.0, zmax=30.0):
)
+def _spergel_hlr_pade(x):
+ """A Pseudo-Pade approximation for the HLR of the Spergel profile as a function of nu.
+
+ See dev/notebooks/spergel_hlr_flux_radius_approx.ipynb for code to generate this routine.
+ """
+ # fmt: off
+ pm = 1.2571513771129166 + x * (
+ 3.7059053890269102 + x * (
+ 2.8577090425861944 + x * (
+ -0.30570486567039273 + x * (
+ 0.6589831675940833 + x * (
+ 3.375577680133867 + x * (
+ 2.8143565844741403 + x * (
+ 0.9292378858457211 + x * (
+ 0.12096941981286179 + x * (
+ 0.004206502758293099
+ )
+ )
+ )
+ )
+ )
+ )
+ )
+ )
+ )
+ qm = 1.0 + x * (
+ 2.1939178810491837 + x * (
+ 0.8281034080784796 + x * (
+ -0.5163329765186994 + x * (
+ 0.9164871490929886 + x * (
+ 1.8988551389326231 + x * (
+ 1.042688817291684 + x * (
+ 0.22580140592548198 + x * (
+ 0.01681923980317362 + x * (
+ 0.00018168506955933716
+ )
+ )
+ )
+ )
+ )
+ )
+ )
+ )
+ )
+ # fmt: on
+ return pm / qm
+
+
@implements(
_galsim.Spergel,
lax_description=r"""The fully normalized Spergel profile (used in both standard GalSim and JAX-GalSim) is
@@ -235,7 +283,7 @@ def __init__(
else:
super().__init__(
nu=nu,
- scale_radius=half_light_radius / calculateFluxRadius(0.5, nu),
+ scale_radius=half_light_radius / _spergel_hlr_pade(nu),
flux=flux,
gsparams=gsparams,
)
@@ -282,7 +330,7 @@ def _inv_r0_sq(self):
@lazy_property
@implements(_galsim.spergel.Spergel.half_light_radius)
def half_light_radius(self):
- return self._r0 * calculateFluxRadius(0.5, self.nu)
+ return self._r0 * _spergel_hlr_pade(self.nu)
@lazy_property
def _shootxnorm(self):
diff --git a/tests/jax/test_spergel_hlr_pade.py b/tests/jax/test_spergel_hlr_pade.py
new file mode 100644
index 00000000..964d8acd
--- /dev/null
+++ b/tests/jax/test_spergel_hlr_pade.py
@@ -0,0 +1,31 @@
+import numpy as np
+import pytest
+
+from jax_galsim.core.testing import time_code_block
+from jax_galsim.spergel import _spergel_hlr_pade, calculateFluxRadius
+
+
+@pytest.mark.parametrize(
+ "nu", [-0.85, -0.6, -0.5, 0.0, 0.1, 0.5, 1.0, 1.1, 1.5, 2, 2.7, 4.0]
+)
+def test_spergel_hlr_pade_correct(nu):
+ np.testing.assert_allclose(
+ _spergel_hlr_pade(nu), calculateFluxRadius(0.5, nu), rtol=0, atol=1e-7
+ )
+
+
+@pytest.mark.parametrize(
+ "nu", [-0.85, -0.6, -0.5, 0.0, 0.1, 0.5, 1.0, 1.1, 1.5, 2, 2.7, 4.0]
+)
+def test_spergel_hlr_pade_time(nu):
+ print("\n", end="", flush=True)
+
+ with time_code_block("root finding") as tr_rf:
+ for _ in range(100):
+ calculateFluxRadius(0.5, nu)
+
+ with time_code_block("pade") as tr_pade:
+ for _ in range(100):
+ _spergel_hlr_pade(nu)
+
+ assert tr_rf.dt > tr_pade.dt