{ "cells": [ { "cell_type": "markdown", "id": "3bb2cd70-e98d-4e8d-b195-64d7e5de6e13", "metadata": {}, "source": [ "# Visualize Caustics\n", "\n", "Here we will demonstrate how to collect caustic lines using caustic! Since caustic (the code) uses autodiff and can get exact derivatives, it is actually very acurate at computing caustics. \n", "\n", "Conceptually a caustic occurs where the magnification of a lens diverges to infinity. A convenient way to measure the magnification in the image plane is by taking the determinant ($det$) of the jacobian of the lens equation ($A$), its reciprocal is the magnification. This means that anywhere that $det(A) = 0$ is a critical line in the image plane (magnification goes to infinity). If we take this line and raytrace it back to the source plane we can see the caustics which define boundaries for lensing phenomena." ] }, { "cell_type": "code", "execution_count": 1, "id": "f716feef", "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import torch\n", "from torch.nn.functional import avg_pool2d\n", "import matplotlib.pyplot as plt\n", "from ipywidgets import interact\n", "from astropy.io import fits\n", "import numpy as np\n", "from time import process_time as time\n", "\n", "import caustic" ] }, { "cell_type": "code", "execution_count": 2, "id": "bdede2df", "metadata": {}, "outputs": [], "source": [ "# initialization stuff for an SIE lens\n", "\n", "cosmology = caustic.FlatLambdaCDM(name = \"cosmo\")\n", "cosmology.to(dtype=torch.float32)\n", "sie = caustic.SIE(cosmology, name = \"sie\")\n", "n_pix = 100\n", "res = 0.05\n", "upsample_factor = 2\n", "fov = res * n_pix\n", "thx, thy = caustic.get_meshgrid(res/upsample_factor, upsample_factor*n_pix, upsample_factor*n_pix, dtype=torch.float32)\n", "z_l = torch.tensor(0.5, dtype=torch.float32)\n", "z_s = torch.tensor(1.5, dtype=torch.float32)\n", "x = torch.tensor([\n", " z_l.item(), # sie z_l\n", " 0., # sie x0\n", " 0., # sie y0\n", " 0.4, # sie q\n", " np.pi/5, # sie phi\n", " 1., # sie b\n", "])\n", "packparams = sie.pack(x)" ] }, { "cell_type": "markdown", "id": "38dd09e3", "metadata": {}, "source": [ "## Critical Lines\n", "\n", "Before we can see the caustics, we need to find the critical lines. The critical lines are the locus of points in the lens plane (the plane of the mass causing the lensing) at which the magnification of the source becomes theoretically infinite for a point source. In simpler terms, it is where the lensing effect becomes so strong that it can create highly magnified and distorted images of the source. The shape and size of the critical curve depend on the distribution of mass in the lensing object. These lines can be found using the Jacobian of the lensing deflection, specifically $A = \\mathbb{I} - J$. When ${\\rm det}(A) = 0$, that point is on the critical line. Interestingly, $\\frac{1}{{\\rm det}(A)}$ is the magnification, which is why ${\\rm det}(A) = 0$ defines the points of infinite magnification." ] }, { "cell_type": "code", "execution_count": 3, "id": "487b3030", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAGfCAYAAACtEPAuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABdV0lEQVR4nO3deXgTVdsG8DtJmxToCt1AylJAkB2K8BUFwRdZRBEXUJGlgIBYUAGF8oJsCpVFRAFZlE0FwQ3xRUEQUFBBFC07yGproWVvSoEuyXx/pAlN6ZLJZJLJ5P55zSVNZuYcAvTuc86ZGY0gCAKIiIhIFbSe7gARERG5DoOdiIhIRRjsREREKsJgJyIiUhEGOxERkYow2ImIiFSEwU5ERKQiDHYiIiIVYbATERGpCIOdiIhIRfw83YGymM1mnDt3DkFBQdBoNJ7uDhERiSQIArKzs1GtWjVotfLVkrdu3UJeXp7k8+j1egQEBIg6ZuHChZg9ezYyMjLQrFkzzJ8/H61bt5bcF6cJCpaWliYA4MaNGzduXr6lpaXJlhU3b94UoiN1LulndHS0cPPmTYfbXrt2raDX64Xly5cLhw8fFoYMGSKEhoYKmZmZsv1+y6MRBOU+BCYrKwuhoaH4589aCA7krAERkbcxXjejZsuzuHbtGkJCQuRpw2hESEgIzuyrieAg57PCmG1G7bh/kJWVheDgYIeOadOmDe69914sWLAAgGWkOSYmBiNHjkRSUpLTfZFC0UPx1uH34ECtpD8sIiLyLHdMpwYHuSYrjEaj3dcGgwEGg+GO/fLy8rBv3z6MHz/e9ppWq0WnTp2we/duyf1wFtOSiIhUwSSYJW8AEBMTg5CQENuWnJxcYnuXLl2CyWRCVFSU3etRUVHIyMiQ/fdbGkVX7ERERI4yQ4AZzs8uW49NS0uzG4ovqVpXMgY7ERGpghlmmCUeDwDBwcEOzbGHh4dDp9MhMzPT7vXMzExER0dL6Ik0HIonIiJygl6vR1xcHLZt22Z7zWw2Y9u2bYiPj/dYv1ixExGRKpgEASYJF3o5c+zo0aMxYMAAtGrVCq1bt8a8efOQk5ODgQMHOt0PqRjsRESkCq6aYxfj6aefxsWLFzFp0iRkZGSgefPm2Lx58x0L6tyJwU5ERCTBiBEjMGLECE93w4bBTkREqmCGAJObK3YlYrATEZEqeGIoXom4Kp6IiEhFWLETEZEqeGJVvBIx2ImISBXMhZuU49WAQ/FEREQqwoqdiIhUwSRxVbyUY5WEwU5ERKpgEiyblOPVgMFORESqwDl2C86xExERqQgrdiIiUgUzNDBBI+l4NWCwExGRKpgFyybleDXgUDwREZGKsGInIiJVMEkcipdyrJIw2ImISBUY7BYciiciIlIRVuxERKQKZkEDsyBhVbyEY5WEwU5ERKrAoXgLDsUTERGpCCt2IiJSBRO0MEmoV00u7IsnMdiJiNzMJMhzV3KdxrcHYQWJc+wC59iJiKgkcgW3lHZ9IfQ5x24h6590cnIy7r33XgQFBSEyMhI9e/bE8ePH5WySiEh2JsFc5qZE3tJPkk7WYP/pp5+QmJiIPXv2YOvWrcjPz0fnzp2Rk5MjZ7NERC6l1kBU2+/JJGglb2og61D85s2b7b5euXIlIiMjsW/fPrRv317OpomIRFNLwDmj+O/dG4fuzdDALKFeNUMdT4Fx6xx7VlYWAKBy5colvp+bm4vc3Fzb10aj0S39IiLf5cthXhbr5+KNAe/r3PYnZjab8corr+C+++5D48aNS9wnOTkZISEhti0mJsZd3SMiH8B5ZvG86XOyLp6TsqmB24I9MTERhw4dwtq1a0vdZ/z48cjKyrJtaWlp7uoeEamYN4WTUnnDZ8g5dgu3DMWPGDECGzduxM6dO1G9evVS9zMYDDAYDO7oEhGplNLDx9txiF75ZA12QRAwcuRIrF+/Hj/++CNq164tZ3NE5KMY5u5nEsyKC3fL4jkJD4FRyVC8rMGemJiINWvWYMOGDQgKCkJGRgYAICQkBBUqVJCzaSJSMQa5MiitejdLvKWsWlbFy/qnsWjRImRlZaFDhw6oWrWqbVu3bp2czRKRCnHBm3Lxz0RZZB+KJyKSgqHhHZQwNC91AZxJJZmljPETIqJCvBzNe3n6z8wMreRNLtOnT0fbtm1RsWJFhIaGlrhPamoqunfvjooVKyIyMhKvvfYaCgoKRLfFh8AQkSIwxEkqk6CBScIT2qQcW568vDz06tUL8fHxWLZs2Z1tm0zo3r07oqOj8euvv+L8+fPo378//P39MWPGDFFtsWInIo9hZa5O/DO909SpUzFq1Cg0adKkxPe3bNmCI0eO4JNPPkHz5s3RrVs3vPHGG1i4cCHy8vJEtcVgJyK34zd+koOpcFW8lA2w3M686Fb0Vudy2b17N5o0aYKoqCjba126dIHRaMThw4dFnYvBTkSy47y5b3L3n7dZ0EreACAmJsbu9ubJycmy9z0jI8Mu1AHYvrZeKu4oBjsRyYpBTt4mLS3N7vbm48ePL3G/pKQkaDSaMrdjx465ufdcPEdELsYgJ08pOpzu3PGWy92Cg4MRHBxc7v5jxoxBQkJCmfvExsY61HZ0dDT27t1r91pmZqbtPTEY7ETkEgx08jQzpK1sF/s3OCIiAhEREU63V1R8fDymT5+OCxcuIDIyEgCwdetWBAcHo2HDhqLOxWAnIkkY6ETlS01NxZUrV5CamgqTyYSUlBQAQN26dREYGIjOnTujYcOG6NevH2bNmoWMjAxMnDgRiYmJoh+OxmAnIqcw0ElppN5kRs4b1EyaNAmrVq2yfd2iRQsAwI4dO9ChQwfodDps3LgRw4cPR3x8PCpVqoQBAwZg2rRpottisBORwxjmpGTSbykrX7CvXLkSK1euLHOfmjVr4rvvvpPcFoOdiMrFQCfyHgx2IioVA528CZ/HbsFgJyI7DHPyVkoeincnBjsR2TDUyZtJv46dwU5EKsAwJ1IXBjuRj2Kgk9qYBQ3MUm5QI+NjW92JwU7kYxjopFZmiUPxcl7H7k4MdiIfwUAn8g0MdiKVY6CTryj66FVnj1cDBjuRSjHQydeYoIFJwrXoUo5VEgY7kcow0Il8G4OdSCUY6OTrOBRvwWAn8nIMdCILE6QNp5tc1xWPYrATeSkGOhGVhMFO5GUY6EQl41C8BYOdyEsw0InKxofAWDDYibwAQ52ofILEx7YKvNyNiOTEMCciZzDYiRSIoU4kHofiLRjsRArBMCeShk93s1DHjydEXo6hTkSuwoqdyIMY6ESuY5L42FYpxyoJg53IQxjqRK7FoXgLBjuRmzHQiUhODHYiN2GgE8nLDC3MEobTpRyrJAx2Ipkx0IncwyRoYJIwnC7lWCVRx48nRArFUCcid2PFTiQDBjqR+3HxnAWDnciFGOhEniNIfLqbwDvPEVFRDHUizzJBA5OEB7lIOVZJGOxEEjDMiUhpGOxERKQKZkHaPLlZcGFnPIjBTiQSq3QiZTJLnGOXcqySqON3QeQmDHUiEuvs2bMYPHgwateujQoVKqBOnTqYPHky8vLy7PY7cOAA2rVrh4CAAMTExGDWrFlOtceKncgBDHQi5TNDA7OEBXBSji3LsWPHYDabsWTJEtStWxeHDh3CkCFDkJOTgzlz5gAAjEYjOnfujE6dOmHx4sU4ePAgBg0ahNDQUAwdOlRUewx2onIw1Im8g1LvPNe1a1d07drV9nVsbCyOHz+ORYsW2YJ99erVyMvLw/Lly6HX69GoUSOkpKRg7ty5ooOdQ/FEpTAJZoa6QpkhiNqIxDAajXZbbm6uy9vIyspC5cqVbV/v3r0b7du3h16vt73WpUsXHD9+HFevXhV1bgY7UQkY6J7l6qBm4PsG6+I5KRsAxMTEICQkxLYlJye7tJ8nT57E/PnzMWzYMNtrGRkZiIqKstvP+nVGRoao83MonqgIBrp7eTpci7evVckNSnyVGRJvKVv455+Wlobg4GDb6waDocT9k5KSMHPmzDLPefToUTRo0MD2dXp6Orp27YpevXphyJAhTve1LAx2okIMdfl5OsjLU7R/DHnfFRwcbBfspRkzZgwSEhLK3Cc2Ntb263PnzqFjx45o27Ytli5dardfdHQ0MjMz7V6zfh0dHe1gzy0Y7OTzGOjyUXqQl8Xadwa89xAkrooXRB4bERGBiIgIh/ZNT09Hx44dERcXhxUrVkCrtZ8Jj4+Px4QJE5Cfnw9/f38AwNatW1G/fn2EhYWJ6hfn2MlncXGc66lxDltNvxe1sz7dTcomh/T0dHTo0AE1atTAnDlzcPHiRWRkZNjNnffp0wd6vR6DBw/G4cOHsW7dOrz77rsYPXq06PZYsZNPYqC7jq+EnhkCq3eFU+qd57Zu3YqTJ0/i5MmTqF69ut17gmD59xMSEoItW7YgMTERcXFxCA8Px6RJk0Rf6gYw2MnHMNCl85UgLwmH58kZCQkJ5c7FA0DTpk2xa9cuye0x2InIIb4c6OQdpA6nyzUU724MdvIJrNSdwzAvGSt3ZVLqLWXdjcFOqsZAdw4Dnch7MdhJlRjo4jDIncPKXVk4FG/BYCfVYag7joFOasJgt+B17KQqDHXH8Nps1+JnSUrCip1UgYFePqWFj9msweWrYbh0uTKuXA2FMTsIOTcqIDdPj4ICP2g0Avz8ClAhIBeBlXJQpfJVRFS5jKpRF+DnZ/J09+/A69w9jxW7BYOdvB5DvWyeDPT8fD+cOF0bR/+ui+Mn6+DU2Zo4mxqDtPRqOJcRhYICf9Hn1GpNqF7tPBrUO4nGDY6jVYsDiL/3D1SNuijD74C8CYPdQtZg37lzJ2bPno19+/bh/PnzWL9+PXr27Clnk+RDGOhlc3eg5+f74dDR+vgjpRn+PNAYfx1ojCN/10NeXslPxrIKC72GymHXEByYjUqVbiDAkAs/PxMEQYOCAh1u3KyA6zmVcOVqKC5cqoL8fD1S/62O1H+rY8uODrbzNKh3Ag932o7HH9mEVs0PQKOO79FEoska7Dk5OWjWrBkGDRqEJ554Qs6miKiQuwL91i09dv8Rh52/tsHPv7XGH381w81bFe7YLygwG/fcfRJ31z2FurXPonaNNNSM+RfVq51HVMQl+PsXONym2axBxoUInP6nJo7+XRf7DzXC7382w4Ej9+DYiXo4dqIe5i4ahnvu/hvP9/0UA579HIGVbrjyt112/zgc71ECpF2LrqzJKufJGuzdunVDt27d5GyCfBAr9Tu5I8wFAThw+B788FM7bN91H37Zey9u3Qqw2yc0JAutmu9Hi6aH0LLpITRvcgi1Yv51WfWs1QqoFn0B1aIv4P42v9tev3I1BNt33YcN33XBt1v/g6N/340xkyZj+jsv4eVhy/DSkGWoUCHXNZ0gxeJQvIWi5thzc3ORm3v7H5/RaPRgb0iJGOp3kjPUb93S48df2uJ/33fC5m0dkX6+qt371aIz0L7tHtz/f3txX+s/UL/uKWi17q97Kodl4ake3+GpHt8hyxiEtV89hgUfJuDE6VhMfutVrFr7FOZNn4LOHXfK3hdW7Z7DYLdQVLAnJydj6tSpnu4GKRAD3Z6cYX49pyI2b+uADd91webtHZB9Pcj2XsUKN9Dhvt34T/uf8WD7X9Cg3knFzWWHBGdjWMIneL7fGny24VFMnD4Wp8/WQo/nVmL4wFWYOXkG9Pp8T3eTSDaKCvbx48fbPXvWaDQiJibGgz0iJWCo25Mj1HNz9fh++wNY93UPfLf1Qbu58ruqnscjnX/AI11+QLv/+w0BAXkub18OOp0Zzz6xAY90/gHTZo/C/A8GYdGKAdh/uCG+WDEUlcOyPN1FcjFW7BaKCnaDwQCDoewVtOQ7GOi3yRHmggD8urcVVn/xONZ/2w1Xr4Xa3outdRZPdN+MHt2+R6vmBzwyvO4qQYE5mD31TTzY7hckjHgHv+69Fw89uRbff94H4VWuerp75EIMdgtFBTuRFUP9NleH+rmMSHzy+ZNYtbYXTp2pZXu9WnQGej/2P/R+/H9o0eSQ4obYperWaQe2f90bjzy7CoeP1UfPfsux5YtnUbHiLU93jcilZA3269ev4+TJk7avz5w5g5SUFFSuXBk1atSQs2nyUgz021wZ6CaTFpu3dcDy1c9g07aOMJt1AIDAStfxxCOb8OyTX6N9/G/Q6dT9+Tdq8Dc2fdYX/+m5Dn+kNMOIpDex/L1XXd4OF9DdptO4787lgqCBIKHqlnKsksga7H/88Qc6duxo+9o6fz5gwACsXLlSzqbJCzHULVwZ6BcvV8bKNb2x9KPnkJZ+l+31+Hv/wMA+6/Dko9+hUsWbLmvPGzSodwprlw1Hl6fWYM0XT6Drgz+id8+Nnu6WKllC3X3/rvk8dgtZg71Dhw4QBO+dmyP3Yai7NtD37W+ChcsG4Itvutvu/FY57Cr69/4CA/usQ/16p13Wljdq93+/I+nlhZjxzksYPXEyuv7nRwQHXfd0t4hcgnPs5FEMdNcFusmkxXdbH8S8Jc/jl99a216Pa3YAwxI+Rq8eG2W9SYvZRZWZ1k0PnRz/ygJ88U13/H2qDt5Z9Dwmj53nlnZJPlw8Z8HHthJ5kCtC/eZNA5as7Ium7X9Ar0FL8ctvreHvn4dnn1yPXd/2xC+beqL/01+6NNTNJfwn17nl4u9fgGnj5wAAFnw4ENnXK8nWli9y59y6lXWOXcqmBqzYySN8vVJ3RaBfywrC0o/6YsEHA3HhUjgAywNVBvf9FC8OWoVq0RcktyFnsIrtgxyV/GPdvsfddU7h71N18OmXPTF0wGqXt+GLPBHqdBuDndyOoS4t1C9erox3Fz+PJaues90VrmZMGl4eugwDnv3cJYvhlBDoxckR8BoNMKT/Grw2+XWs+fIxBruX41C8BYOd3MaXA90VFXrmxXDMW/Q8lqzqixs3KwIAGjU4jtEvLkHvxzaKekranf3znj8bVwd8z4c347XJr+O3fS1x8XJlRFS5IvmcvnypmyerdV7uZsFgJ7dgqDvvwqUqeHvhMCxd9ZztVq9xzQ5g/Cvz8fBD2yXdFc6bAr04M8wuCfeYu86jScOjOHjkHuza3QZPPLLJBb0jTxAkVuwMdiIH+WqoSw30K1dDMG/J81j4YQJyblgWdrVu+RfGv7IAXf+zw6k7w3lzkJfEVeH+f63+xMEj9+D3P5sx2CXg3LoyMNhJVgx18XJuVMCCDwdi7vtDkWUMBmCp0CePnYuHOuxkoBfjinBv1vgwAODoiXqu6BJ5iADLMxCkHK8GDHaSBQNdvPx8P6xY8zRmvDMSGRciAQBNGh7FpFffwSNdfhAd6GoO8+KkhntszVQAwOmzvNW1s5RQrZuhgYZ3nmOwk+v5YqhLCXRBAL7Z3BkTpo/DydO1AQC1aqRi6ri30euxjaLn0H0p0F2lWnQmAODi5Soe7gmRdAx2cimGuji//9UUY6dMxO7fWwEAIsMv4b+j5mPQc2uh1+eL6IPvfe7FSanaQ4KNAIBrWcEwmzVe/ZhaT1BCtQ5wVbwVg51cxtdCXUqgp5+PwoTp47D2q54AgAoBN/HSsGV4NXEJggJzRPTBtz5zuRj0eQAAQdDCbNZCqzU5fS5fvtTN08yCBhpex85gJ+l8LdAB50P95k0D5i4aijkLXsDNWxWg0ZjRt9dXmDLubdxVNVNE+773mctJLZWaJyilWqfbGOwkia+FurOBbp1HHztlAv5JiwFgeXTq229MQ8umh0S071uft1jODsffuGm5P4Benws/P+erdfIsQZC4Kl4lMzD8UYucxlB3zPETsej+zEd4evBi/JMWg7uqnsPHi0Zi+9e9HQ51uR+I4uuslxUGB0p7dKuvDcMrrVpX8kNgevTogRo1aiAgIABVq1ZFv379cO7cObt9Dhw4gHbt2iEgIAAxMTGYNWuWU20p60+FvIYvhboZglOhfuNGACbPHI1Wnb7D9l33Q6/PxdiXFuLArofQ67FvHbp8zROBbhIEl23e4nym5fLC6KiLHu6J91BaqCtdx44d8dlnn+H48eP48ssvcerUKTz11FO2941GIzp37oyaNWti3759mD17NqZMmYKlS5eKbotD8SSar4W6Mzb90BGvTJhiG3bv+p8dePuNqahTK1VE2+77nOUKYet5dc7cVcdJzgzH//PvXQCAu6pmyNElchMlr4ofNWqU7dc1a9ZEUlISevbsifz8fPj7+2P16tXIy8vD8uXLodfr0ahRI6SkpGDu3LkYOnSoqLYY7OQwBnr5zmVE4tVJk/DVxocBANWrncPbb0xDj65bHK7Q5eaJSrpom+4MeUedOBULAKgXe8bDPSEpXLUq3mg02r1uMBhgMBgk9a2oK1euYPXq1Wjbti38/f0BALt370b79u2h1+tt+3Xp0gUzZ87E1atXERYW5vD5OZZCDmGol3OMWYMPPuqD5g9sxVcbH4ZOV4BRLyxFyk+d8Vg3ZYS6UobHldCH4g4fuxsAUL/eSafP4Uvz60odhrcunpOyAUBMTAxCQkJsW3Jyskv6N27cOFSqVAlVqlRBamoqNmzYYHsvIyMDUVFRdvtbv87IEDeSpMw/HVIUXwl1Z+fS/z5ZG52fWoORSW/CmB2Ee1ukYPfmHkie9BYCK91woF155tGVPN+tpP4IApBysBEAoEWTwx7ujfIpNdRdKS0tDVlZWbZt/PjxJe6XlJQEjUZT5nbs2DHb/q+99hr++usvbNmyBTqdDv3794cgw78FDsVTmXwp1MUqKNBh7qKhmD73JeTmGlCxwg1MGz8Hwwd+BJ2u/M9NrgpdSaFZFpMgKGJY/tSZWrh0pQoMhlw0bnCs/ANK4EvVupJZqm4pc+yW/wcHByM4OLjc/ceMGYOEhIQy94mNjbX9Ojw8HOHh4bj77rtxzz33ICYmBnv27EF8fDyio6ORmWl/Lwvr19HR0aJ+Hwx2KhEDvWyHjtbH0FGz8OeBJgCAhzr8hPkzJ6JWTLoDbcpTnXsjJYT7zt2tAQD3tkhBQECeR/uidEqv1t29eC4iIgIRERFOtWU2W74P5ObmAgDi4+MxYcIE22I6ANi6dSvq168van4d4FA8+TBnq/RZ84cjvusG/HmgCUJDsvDhu2PwzeqBHgl1JQ6ze5sdu+4DALSP/83DPVE2pYe6kv32229YsGABUlJS8M8//2D79u149tlnUadOHcTHxwMA+vTpA71ej8GDB+Pw4cNYt24d3n33XYwePVp0e6zYyQ4r9dIdO1EHQ16Zjd//ag4A6P7QD1gwawKqOnDtsxyBrhaurtrFXOpmMmmxvTDYOz2wy8n2OAyvFAKkPVNdrn9VFStWxFdffYXJkycjJycHVatWRdeuXTFx4kTbavuQkBBs2bIFiYmJiIuLQ3h4OCZNmiT6UjeAwU5F+EKoO7vi/f3lAzBh+jjk5hoQEmzEnGnT0LfXV+WudndloMsR5s70T8pzz5Xm19/jcPlqZYSFXsO9LfZ7ujuK5S3VulKvY2/SpAm2b99e7n5NmzbFrl3O/YBZFIOdAPhGqDvj33PRGPLKbOz42VLVPdThJyyaMx7Vq5V/+YkSQ90VfSp6Dm8P+Q3fdQUAPPzQNvj7F4g+ntU6KRGDnXwi1J2p1L/a2A2JY6fj6rVQVAi4ibcmJWPogE/cVqVLDXN33OzG2obUgPfEIjqTSYsvvukOAHii+ya3tu1NvKVaB6DcsXg3Y7D7OIb6nbKvV8KY1yfho3W9AABxzQ5g5YJRqFen/LuSKSHUPfGwGGefquZqYvqwfVdbZFyIROWwq3iog/ThTzXyqlAHAKkPclHJ43sZ7D5M7aHuTJX+54HG6Df8PZw6UwsajRmvjliMSa/OK3eY1hVh6m1hXloflBDwjrD+4Nb7sf9Br88XfTyH4ZWHj221YLD7KIZ6sf3NGry3dBBeT34N+fl6VK92DisWjEK7//vdgbY8F+pKCHRvdOlyGL7Z3BkA0P+ZL0Qf7wuh7nXVOtkw2H2M2gMdEB/ql6+E4vlX5mDTDw8CAHo+vBnvzx6PymFZ5bQj7bN0JsxdNtRfymekc0FgKWVYviwfrXsKubkGtGx6EC2aHPJ0dxTHW0Ndqavi3Y3B7kMY6nfa/XtL9H3hPaSfrwaDIRdzpr6B5/utcetlbI6S0mZpQV7efq4Iendx9IeJggIdFq3oDwAY0n+1Qw/oIS8haKTNkzPYyZuoPdTFBrogAO8tHYwJ08eioMAf9WJPY/WSEWjaqOx7hbu7ShfbnqMBLvZ8YgNeyVX71991QVr6XYiocgnPPL6h/AOKUfswvLdW63Qbg90HMNTtXcsKwtBRs/DN5i4AgF49/of35/wXQYE55bSj3FB3daCXdH5vqt5LIwjA3PeHAQCGDliNChVyRR3PUFc2Lp6zYLCrnJpD3ZlV7wePNMAzQ97HqTO1oNfnYs7UN8sdjlVqoMsd5qW1J0fAS72G3dHRgR9+uh9/HmiCCgE3MXzgR5LaJAXidewAGOyqpuZQd8bar3pg+KvJuHmrAmLuSsfaD19EXLODZR7jzlB3pC0pYW4uoS9aH5pgFgRg+tsvAwAG912L8CpXRR3Pap28BYNdpdQc6mIr9fx8P4x/YzwWfDgQgOW2sCsXjEKVytfKaUfCYjUPBnpJAe7ovo4EvdKG5R2v1tthz744BATcwpjExTL3yruoJdS5Kt6Cwa5CDPXbLl0Ow3MvzMdPv7QFAIx7eQEmvToPOl3pn5G3BbqYIHf0XGqr5AUBmDJzDABgaP/VDj2Rryi1V+uqopLhdCkY7CrDUL/t4JEGeDJhKVL/rY7AStex7L1X8Vi3LeW0oYxQLy/QXRnmzlBa1V6e9d92w779TVGpYg5eHcFqvSi1VOt0G4NdJdQc6ID4UN+wqTMGjXwbOTcqIbbWWXyxYhga1j9RThvyh7qUQHckzB0dsi8vlN1ZuUtZOOfIMHxenj8mTB8LAHh52DJEhl8W2Yb3/AAjltpCnUPxFur6UyVVEhPqggC89e6LeHrwYuTcqIQH2/2Mn799XLZQNwmCx0PdBMG2OUrs/krk6Nz6ohX9ceafmoiOvIDRLy6VuVfkUYILNhVgxe7l1Fypi63Sb93SY9iYmVi3/jEAwIuDVmLWlOnw8zOV0Ya8VbocYe7KQJbzEjZHyP2o1ouXK2PGOyMBAJNem4vASjccPpaVujfSFG5Sjvd+DHZSJLGhfuFSFfQeuBh79sXBzy8f82ZMxvN915bThnvm00s9Rym/R2cDvbzfjTu/lSvlrnMTp49FljEYzRsfwgAnHvZC5I0Y7F6Klfptx07UQc9+y3A2tQZCQ7Lw6QcvouP9u8tpw/mhd0eUdX4xgV7qvg71ovRjisduaZW7WRAUuULekR8cftvXHKvW9gYAvDN9SplXQtx5fuX9nl1FvdU6eIOaQgx2UhSxof7Tr23w9ODFuJYVgto1/8GGjwfh7rpnymlDvlB3ZuhdTIVeXs9NguDQ8LYZ8lXwjlbrzg7DO3L+/Hw/JI6dAQDo//TniL/3T6faUhtVhzrAYC/EYPcyrNRvW7f+UQwZNQt5eQb8X6s/8PmKFxBR5UoZ51dWle5IhV7SGcvrT2nvFw9SKeEudU5e7rn1eYufx6GjDVAl7ApmTHxL1LFqrdZVH+pkwz9pUgSxK9/fWfQ8BiS+i7w8A5545DtsWtdPkaFe2upzZ0K9tBX4pmKbKyltGN6Rav3k6VqY/s5LAICZU6aLvnUseTHrY1ulbCrAit2LqLVaFxPqZrMG46ZOwPwPBgEARjy/ArOmvAmttozV5R4aencm0EsKc/t9y1d8H12xc8ldLcvFkVA3mzUY/toM3LoVgP+034Xnnlov4vze+bk4wleqdT7dzYLB7gXUGuiAuFDPy/PHkFGzbJezvTVpOl55YVk551duqJdVoRdtv3hQe3r1e1nD8I6Er5w/WLy/fAB27f4/VKqYg/kzJ5b51D5f4SuhTrcx2MljxIR6zo0KeHbIQmzZ0QF+fvn4YN5YPPvEhnLO7/5QlyvQS55rv/1rXZEAs+5r/XZuwu2qXU6eDvUTp2rj9RmvAQBmTJyJ2JppsrVFCsXFcwAY7IrGSt3i6rVgPN7/Q+z5oxUqVriBtR++iM4dd5ZzfvGfnasDHbAP9dIC3ZEwNxU7/R3hXULIF10c565wl4Ojq+AHjpyLm7cqoOP9v2BI/9Uizq/Ost4nK3Wp8+ScYydyjphQz7wYjkeeXYWDR+5BaEgWvv54EP6v1V/lnN+9oV7ivmVU6aWdxVTC+9bALu2Ykq5NNwklh7srKPXBL9PfGYk/UpohNCQLH8x7rcw1F0Rqx2BXILVW6mIvZ0tLr4qHn/4YJ07HIjryAjZ+OgCN7zlexvnlGXoXM5fuaKCXVKXfruCLfV0kTE1FKgqdpni/Cx/cAvtwd1Tx0BazIl6uYXhHzvvLb60w670XAQDvJb+O6tUyRJxfmT+oSOWT1ToAjWDZpByvBgx2UqQzqdXRtddq/JMWg5i70rH5s76oU/sfl7cj5dawYu7ZXlKoOxroplKGB62v3xnw0pUU6qVV656cW79yNQQDEufBbNbhuV5fonfPjQ4fy1BXIc6xA2CwKwordYuTp2uhS69PkH6+GurUPotN6/qiRvVzZZzf8/PppVXqjga6GXcGuRkau2q9OJ21H9amNUVec0DRb//lDbFLGYKX6w5zggC8MOYt/HuuGurGnsG86VOcakdNfDrUAc6xF2Kwk6IcPxGLrr1X43xmFBrUO4FNn/VF1aiLnu5Wmcpb9Q6UHOolzZ+bBI1doJtL+Uaj1QgwQQMdBJihgbbYkL+7vr3L9bAXR847/4NB+GZzF+j1ufj4/ZcQFJgj4vzq+AZOVBIGu0KosVoXW6kfPxGLLr3WIONCJBo1OI5Nn/VFZPjlMs6vrEq9pPn08gLdBI1dmJuF26Fu+do+4LQay5l1gmCpzkuo0kuKROuKeGv1XFa1XnwY3t1D8I6cd/fvLfHfN8cBAGZNnoEWTQ+LOL86Q93nq3WAQ/GFGOwKoMZQF+v4iVh0fupTZF6MQNOGR/Ddun4euRWoO0K96Bx60VDPLwzxfEF3+/3CkNNZjyrcR6dx7OaxYr/Vu/IWsnLNq2dcCMdzwxagoMAfvXr8D8MSPnb4WIa6yjHYAfBe8R6nxlA3QxBVrZ84VRtde69G5sUINGl4tNxQNxf+J0Zp91kvft4Sjy1h5buYUDfDPtRN0CBP0CJf0CIfWtwSdIWbH3IE/e3NbMAts79lEyxbHnSFgV/O/HuRt8RU63ecR0K17ozyzpuf74e+L8zHuYxoNKh3Au/P+a/Dd5djqJMS5Obmonnz5tBoNEhJSbF778CBA2jXrh0CAgIQExODWbNmOdUG/0Z4kBpDXazT/8Sga+9PcD4zCo3vOYZN6/qWG+pyEBPqJb1/O7wtP0BYH8hiDfSioV5SoN8S/JFjNuCG2YCcwu2GYLAFer7gB5Mg7p+rXNW6J4fgx039L37e0wZBgdn4bNkLoubV1YihXozggk1mY8eORbVq1e543Wg0onPnzqhZsyb27duH2bNnY8qUKVi6dKnoNjgU7yFqDXUxlXpaelV07bUa6eer4p67/3aoUhfLlXPqjqx8t7wO2+tFA9067J4vaAsDXoc86HDL7A8ztMgxGwpf9yvsuwZ6jck2r67XAGZBC53mdovW+XWdxjLnroVz1bqj8+qOkCvUV619Cu8vTwAArJg/GnfXPSPi/Oqr1hnqJVD4qvhNmzZhy5Yt+PLLL7Fp0ya791avXo28vDwsX74cer0ejRo1QkpKCubOnYuhQ4eKaofB7gEMdcsd5R5++mOk/lsddWPP4Lt1/Vy+UM4RzoS62Pl0a6jnQwuzoMEtwQ8maJBjNsAMLYzmAOQLfsg2BcAELfIFHXQQ4K8pALS50AOw5pJt8RwEW6hrS/jctZA2BF+a8gJYrnn13b+3xIhxbwIAJox+F4902ebwsWoMdZKX0Wi0+9pgMMBgMEg6Z2ZmJoYMGYKvv/4aFStWvOP93bt3o3379tDr9bbXunTpgpkzZ+Lq1asICwtzuC3+yOdmDHXLvd8feXYVTpyORY3q/2LTurIvaXP3HeXKq9RvH1/8uJLn062hni9okQetpUovHGI3mgNwzVQJVwoCcakgCJfyg3C1oBKyTQG4YTbglqC3LaCzsoa6ViNAC+GOar28f9RyrYKX63r1f/6thqcHL0Z+vh6Pd9+ECaPfE3FudYY6q/WSWe88J2UDgJiYGISEhNi25ORkSf0SBAEJCQl44YUX0KpVqxL3ycjIQFRUlN1r1q8zMhy/myLgpmBfuHAhatWqhYCAALRp0wZ79+51R7OKw1C3PKXt8f4f4uCRexAdeQGb1vVDzF3nyzi3c8PvcoS6CcIdw+9F59SLzqfnF1npnivokC9ocUvwww2zAddMFXHZFIi0/Co4kxuJlOs18LuxJvZerYXD2VVxMS8I102W6kAHM/w1BdBrTJb/wzI0ry0Mc2eG4Ity5aVtzijvvMbsQDzRfxkuXApH04ZHRN0HnqHug1w0x56WloasrCzbNn78+BKbS0pKgkajKXM7duwY5s+fj+zs7FLP42qyD8WvW7cOo0ePxuLFi9GmTRvMmzcPXbp0wfHjxxEZGSl384rBULesaO4zdCH2/NEKoSFZ2PjpAFluE1seV1XqxVe+3/n+7dXr+YIO+YWL5KwL5CzVeUUcz4pEboEfCsxaVNLnoZJfHirpLNGs05gLg7vw/4VfW14reQjeSuy8uhRyzKsXFOjQ94X3cPhYfURHXsCXq4YgsNINB8/NUCfnBQcHIzg4uNz9xowZg4SEhDL3iY2Nxfbt27F79+47hvNbtWqF5557DqtWrUJ0dDQyMzPt3rd+HR0dLar/sgf73LlzMWTIEAwcOBAAsHjxYnz77bdYvnw5kpKS5G6eZCQm1M1mDYaOnonvt3dAhYCbWP/RYJc/0EWuhXKOXM5mvelMfmGg5Au3h97zBZ1tlXuO2YB/8yrjUn4gfjpfFxcvBKPWOi0QqMXlJ3LhV9mIEP+bCPG7iRDdDQRpb6GSNhcB2nzoYYIe5nKH4HUoOWjLC/WSqnVPrYAXBOCVCVOxZYfl78sXK4eWObLjCxjqyhMREYGIiIhy93vvvffw5ptv2r4+d+4cunTpgnXr1qFNmzYAgPj4eEyYMAH5+fnw9/cHAGzduhX169cXNb8OyDwUn5eXh3379qFTp063G9Rq0alTJ+zevfuO/XNzc2E0Gu02NVBjtS72rnITp4/Fp18+Dp2uAJ9+kIj4e/+UqWfiSanUS37/zko9H5ZwzzYF4FJ+IM7fCsHF8yEwpOlRITULAZfyodEI8NOaUVGbh4raPARo8+GvKYC/pgA6mMscgi8a6kU5+g/c2cV0ci2WmzX/RXz4cR9oNGZ89P7LaNX8gMPHqrFaZ6g7RgOJc+wy9atGjRpo3Lixbbv77rsBAHXq1EH16tUBAH369IFer8fgwYNx+PBhrFu3Du+++y5Gjx4tuj1ZK/ZLly7BZDKVuCDg2LFjd+yfnJyMqVOnytklt2OoA+8v64+5i4YBABa/nYSu//mxjHMrZ6GcM6vfrZey3RJ0yBd0tuvTs80VkJkfgkv5gfj9Uk2kXwpFnU/M0F+8gtPPVEFutXzcV/M0ogxG1DBcRpDuFkJ1OaikyYO/pgABmgL4a8zwh7ncUC8+BO/MYjlPXqv+yeePY/JbrwIA3p42DY92/UHE+RnqPk3hl7uVJSQkBFu2bEFiYiLi4uIQHh6OSZMmib7UDVDY5W7jx4+3++nEaDQiJibGgz2ShqEObPz+PxgzaRIAYGrSHPTr/ZXL++Tso1cdfeyqmOvUrfd7NwvawkvX/Cxz62Z/3DDrccOsR77ZMjueW9kfZn0QcqvnIyIqC1UDshDufx1BOsvwuz9MhRW7yTavbn1Ea/Fv9a5eLOcpm7d1wLDRMwEAo4cvwYuDP3L4WIY6eYtatWpBKOH7VtOmTbFr1y7J55c12MPDw6HT6UpcEFDSYgBXXCuoFGoMdbH27W+CfsPfgyBoMbjvGowd+X6Z+8sxr+6Ke78XbcvRm8/cEvyRVzi3nm0OQJapIrIKKiArvwIC/AoQHnod5qFGCDoTHg07hyr+Oaiuv4JK2lyE6m4gQJOPAE0+/DWmUit1oOzr1ZV0WZsjlfqeP1rg2SELYTL54dkn1+PNCY7fTlONoU5OkHr3ODfcec4dZP1xUK/XIy4uDtu23b6ZhNlsxrZt2xAfHy9n0x6l1lAXU63/ey4aTyUsxc1bFdC54494d8bkMu/pLddiOUc5OgRf9DW744vcv91cGO5maAtD3/LPzF9jgl5bgCoBOYisdB11Qi7hntBM3GW4hnD/bATpbqKiNrcw0AvsVsCXVKmXdxOaojx9WVt5Dh5pgJ79ltv+viydO46XtbFaF89Fl7t5O9mH4kePHo0BAwagVatWaN26NebNm4ecnBzbKnm1YahbrlV/csAHOJ8ZhUYNjuOTxS/Bz6/0p5Ep6a5y9scX7l/CYrmyhuDzYFkwl1e4cM56g5lAXS78tSZUM2RBpzEjRHcTBm0+QnU5CNDkW4bfNSYEaPIL7zxXfqVetHpWYqXuyHlPnamJR/usxLWsEPxfqz/w6QeJ8PcvcPDcDHWi4mQP9qeffhoXL17EpEmTkJGRgebNm2Pz5s13LKhTA4a65bK2Ia/Mxv7DjRAZfglfrnoewUHXyzi3chbLlfT+necteQjeLJT8xDUdzNBrClCxMLQrQmO7XWzFwsvYrPPp1kC3XtLmyEI5K7krdblCPfXfauja+xNkXIhEk4ZHsf6j51Gp4k0Hz81QJ3tF7x7n7PFq4JbFcyNGjMCIESPc0RS5mNjFcsnzRuCrjQ/D3z8P65a9gFox6S7vk7OhXu55Hbhe/c62bge6NdTNRZ7CptWYEYB8aLVm6DS3AKDw8jUBlbS50MJsu0bdukjOX2MWdZ263JW6s8o77/nMCDz8zMdIS78L9WJPY+OnAxAW6tglrgx1KhHn2AEobFW8N1Njte7MCvg35owCAMx/6/Vyr1V357y6mMVyYobgS6LVmKEVzNBb39YAAdBCC+tT2ky2QLfcMvZ2oBf9vy2wHajUvW34/cKlKujW+xOcPF0bNWPS8N26foiKuOTguRnqRGVhsLuAGkNdrBOnamPQS3MBAC8kfISEZz8vc39331nujv1E/oAgtreWW8GabWFufQ1AYZDbBzoAp1a+A95XqV+8XBkPP/0xjp2oh7uqnsPmz5/jXeUY6q7Bih0Ag10ytYa6mGr9ek5F9B68CMbsILRt/TtmT32zzP09vViurPfLq9ZLY32EqhkCdDDZQtxU+A1bV3ivd8AS7Jb/Wyvz28PuABy+8YzlPfG3ibWcwzOXtF26HIaHn/4Yh442QNWoTHz/eV/UrvGvg+dnpU5l4xy7BYNdAoa65Z7eL742A0f/vhtVozKxZqnjK5rFkOvSNrvXS9tfxLmtz0uHYA10k+11a3BbH7laPNTLGnp3Z6jL5eLlyujW+xMcOtoA0ZEX8P3nfVA39qxH+qIUDHWSA4PdSQx1iw8+eg6ffd0DOl0BPlk8EtGRZc+TKunhLiW9X7S9ki7QK6kVHSw3mdZCgM46766x31OrKRLqhW3qC/eREuiWc7tu6F2uSt06p259Utv3n/fB3XXPOHh+VurkIC++pawrMdjJafsP3YNXJ08EAEyfMBP3tfmjzP3dPQTvDloIMBdewmaCxhbgxRUNdF2RfYrPpVv2Lfy/A1W4Kx+9KpdzGZHo1vsTHD9Zt3D43fFQJxKFc+wAGOxOYbUOZF+vhL4vzEdengGPdN6Kl4ctK+fc8lyvXuaxTsytl0dbuK+1Sgdg+2ZQ0vPRgdtBXnQO3fI67L+2vS7+UraS9rl9Ds9V6tbr1E+frYXq1c5h82d9HR5+Z6VOYnGO3YLBLhJD3WLM65Nw4nQs7qp6DkvmjivzdrHOcuUQPCB+JbyjrMFd9PK3olW5o4Fuee/OUHeoDxJCUK5Hr544VRvdnv4Y/56rhlo1UrH58+dkua8BEdljsIvAULf4amM3fLSuF7RaE1YuHIUqla+Vc37PXtoGiL/DXEl0GsvK+NuBKxR7/85zFg9zu9ds7xWrxoue08n5dMt5PFepHzzSAI88uwqZFyNwd51T+G5dP1SvluHg+Vmpk5M4FA+AwU4ipZ+PQuLY6QCAV0csRrv/+93DPbqTM6FdFuvwO3A73K2vW5lx52vW/VHCe8VvNFPiseWsei+LJ0P9171xeLz/MmQZg9G04RFsXDsAkeGXHTy/OkOd3ETiUDyD3YewUrcQBGDoqFm4ei0ULZsexOtj3i3n/Or53IqHu1XRkNeVkEmlhXlZQ+5KnE939Nzfb38Azzz/Pm7eqoD4e//AV6ue521iWamTm/FvHDls2epnsG1nOwQE3MKK+aNkuV4dcP0wvLN0sA9jLUquyK1b8f3s7u1u219TaqjroHFo6F2pob76i554svBRvV0e/BHfftrf50Od3IyPbQXAYC+XWqt1sf75txqSpv4XADAtaQ7q1zstSztSQ93Vw/Al0ZazAbfDvKxAL7rivaRAV8J8uiPnFgRg7vtDMPiluSgo8MczT3yNL1YMRcWKtxw8vzpDXafRslp3NwY7AA7Fl0mtoe7MEPzIcW/iek4g2rb+HSOeX+FAG9772VkD0PpDRvH58OI3rilpvrzoeazKG3IHSp9Ll3oXOblC3WTSYuyUCVi4bCAA4OVhHyL59WRotY79HVNzqBN5CoOdyrV2/WPYsqMDDIZcLJ6TVO43baWFuqsvcystyG3vlxCizoa61Cq9tP44orzz37xpwMCR7+Dr77oCAGZOnl7u/Qzsz6/OUCfP4XXsFgz2UrBat7h8JRSvTbLcXW78K/NlvWOY1PvBSxmGL7o4zqqkQLRV8eWEZWmR6IpAt5zfs/PpFy9XxlMJS/HbvpbQ63Ox7N1X0euxb0W0oc5QZ6VOSsBgL4FaQ90ZE2eMxaUrVdCw/nGMHv5BufvLfYc5OUcDSgr34uQO9NL2vX1+z1bpAHD8RCx69l+GM//URFjoNXy2fJioyx7VGupESsFgL0bNoS62Wv/9r6ZY+WlvAMD8t16HXp9fzvmV+dlpNZo7huOt4Vm8yi/p2vRSz1vO+6XOiysw0B1tY9vO+9Bn6EJkGYNRu+Y/2PDxIFGjOGoNdVbqCsEb1ADgqngqhdmswagJUyEIWvR56qtyH/Dilj65+QcHR1a/l6S0S9JKWule9Bip5A71paueQ4/nViDLGIz4e//Azo1PMtRJUaxz7FI2NWDFXoiVur3VXzyOP1KaISgwGzMmviVDr25z5bPWxSqtcnfmHCUp645xrphHt51LxqH3/Hw/vDrpdSxZ1Q8A8OyT67F4zngYDHkOtqHeQGelrkAqCWcpGOx0hxs3AjBl5hgAQNLLC8t9xjqg3GF4R4kNeEcqbHeEutxV+uUroegzbAF++qUtNBozpiXNwasjFjv80B81hzqRUjHYwWq9uHeXDkb6+aqoUf1fjHh+pQNtOP/5uatat4ZseZe+SR0SL++e7kqo0h1t59DR+ug1aDHO/FMTgZWuY8X80Xi06w8i2lBvqLNSVyjOsQNgsKs61J1x+Uoo5r4/FADwxvjZDg+3KoEOGrfcea4kUgMdUFaor/+2K55/eTZyblRC7Zr/4MuVQ9Gw/gkRbTDUyf14HbuFzwe7mjlTrb+9cBiyrwehWaPD6PXYRhl65VmOVu6OnMMR3hboJpMWk996FXMWvgAA6Hj/L/hk8chyH81r3w5DnciTfDrYWa3bu3CpChavtCyQmjJurkO3BVXa3LqjVXvxcC4t6MU+LrVoP8rtg5sC3dG2Ll6ujITEd7BtZzsAltvDTp8wE35+xW+iW1ob6g10gKHuFTgUD8CHg13toe5Mtf7O+0Nw42ZFtGq+H13/s0OGXimXswFelKPz82ICHXBPqO/9sxmeHbIQ6eeroWKFG1j8dhJ691TfiI2zGOregUPxFj4b7GTv6rVgfPBxHwDAhNHvObTqWWq1LtfCOVdcwia2LUcoMdAFAVi4LAHj30hCfr4edWPPYN2Hw9Gowd8i2mGlTqQkPhnsrNbvtPSjvrieE4gmDY+qplqXczGd2NXzSgz1LGMQhr+ajK82PgwAeOKR77D47SQEB10X0Q5DnRSEQ/EAfPDOcwz1O+Xl+WPxCsvc+isvfODwNcruJjYcgdLvAufMOYpujtAW+c+hdgqf2S51gZwj7aUcbIi23Tbgq40Pw98/D29Pm4rVS0Y4HOpaaBjqpDwKfh57rVq1oNFo7La33rK/+deBAwfQrl07BAQEICYmBrNmzXKqLZ+s2Mne5990x/nMKFSLzkCvHo49oUtpi+bKUzyMS6vkXXFrV8C5H0Lc0aYgAItX9MO4af9FXp4BMXelY/WSEWjdcr+IdtQd6ABDneQxbdo0DBkyxPZ1UFCQ7ddGoxGdO3dGp06dsHjxYhw8eBCDBg1CaGgohg4dKqodnwp2tVfrzlq6qi8AYEj/1eU+6MXTtNC65IcKVwW4lbNBLnXIXUzbl6+EYtjomdi45SEAwKNdtmDJ3HGoHJYloi11hzoD3bspffFcUFAQoqOjS3xv9erVyMvLw/Lly6HX69GoUSOkpKRg7ty5ooPdZ/4W+0KoOzMMf+hoffy2ryX8/PIxsM86GXqlbmKG2otzZ6jv+DkerR/6Fhu3PAS9PhdvT5uKz5a/wFAvgqGuAi4aijcajXZbbm6uS7r31ltvoUqVKmjRogVmz56NgoIC23u7d+9G+/btodfrba916dIFx48fx9WrV0W141MVu5o5E+oA8PFnTwIAHum8zaF7wiuBNcw8NR0gdZjdnYGem6vH1Fmj8M7iIRAELe6ucwofL3oJzRofFdGWugMdYKirhosWz8XExNi9PHnyZEyZMkXCiYGXXnoJLVu2ROXKlfHrr79i/PjxOH/+PObOnQsAyMjIQO3ate2OiYqKsr0XFhbmcFs+Eey+UK07o6BAh0+/7AkA6Nf7C4ePU8r8uquG5cW0J4UrAl1MP46dqIOExHeQcqgxAOD5fmswc/J0VKp4U0RbDHXyPWlpaQgODrZ9bTAYStwvKSkJM2fOLPNcR48eRYMGDTB69Gjba02bNoVer8ewYcOQnJxc6vmdpfpgZ6iXbsfP8bhwKRxVwq6gc8edbm9fp9FIvpZdrurdVYvfXBXmgON9Mps1WLS8PybMGIdbtwJQJewKFs0Zjx7dtopoS/2BDjDU1cZVc+zBwcF2wV6aMWPGICEhocx9YmNjS3y9TZs2KCgowNmzZ1G/fn1ER0cjMzPTbh/r16XNy5dG9cHuC5wdhl//bTcAwOOPbIa/f0E5eytb0dBzJuTlWMXuiVD/599qGD7mLWzfdT8AoNMDO7H0nbGoFn3BZX1RC4a6Crn5OvaIiAhEREQ41VRKSgq0Wi0iIyMBAPHx8ZgwYQLy8/Ph7+8PANi6dSvq168vahgeUPniOVbrpTObNfh2SycAQM+HN3u4N66ldeI/V3HFtehFOdo/QQBWrOmNVg9uwvZd96NCwE3Mmz4Z/1uTICrUfeH6dIChTu61e/duzJs3D/v378fp06exevVqjBo1Cn379rWFdp8+faDX6zF48GAcPnwY69atw7vvvms3hO8oVuxeztlq/a+DjZF5MQJBgdloH/+biPZc+8OSK4bjlcCV1TkgbgQh/XwUEsfOwOZtHQEA/9fqD3w4byzqxp4V2SYDnbybUi93MxgMWLt2LaZMmYLc3FzUrl0bo0aNsgvtkJAQbNmyBYmJiYiLi0N4eDgmTZok+lI3QMXBzmq9bD/8ZHmCV4f7div+2nUlc3WgA46HuiAAq9b2wtgpE2HMDoLBkIspY9/GS0OXQ6cT9/ffF0KdfIBCbynbsmVL7Nmzp9z9mjZtil27dkluT5XBzlAv385f2wAAHmz3i4d74n1VuxxhDoir0k+drYERY6djx8/3AQBat/wLS+aOwz13nxTZpm8EOit18iWqDHZf4ewwfEGBDnv+aAkAaBe/15Vdcpo3hLsSAr2gQIf5HwzEtNmjcPNWBQQE3MLk1+aySi8DQ92HKLRidzfVBTur9fIdP1kHOTcqISgwGw3rO/54TrkpLdzlCnIrsYv2/jzQGCPGTsefB5oAADrc/ysWzvov6tRKFdmubwQ6wFD3NZrCTcrxaqCqYPelUHe2WgeAvw42AgA0a3wEWq1yghRQRrjLHeiAuFA3Zgdi2uxReH95f5jNOoQEGzFz8nQMeOZz0U/iY6gTqZ+qgp0cc+xEXQBAIwVV60UVDVa5Q94dIW4ltkIXBODL/z2MsVMm4lyG5QYVvXt+g1lT3hR9+18GOvkEDsUDYLD7pFNnagEA6tU5I/pYd9/G1Rq8rg54dwY6ID7UT5yqjVcmTMG2nZarF+rUPot50yfjoQ7iV8wy1MlXKPVyN3dTTbD70jC8VOnnLdVfzF3nPNwTx7k7iF1FbKBnX6+E5HdGYP6HA5Gfr4fBkIvXRizCq4mLERCQJ7Jt7/zMnMVQJ1bsFqoIdl8LdSnz6wBwPtNyC8PoSN5mVC5iA91s1mDt+h6Y8GYSzmdanujUrdN2zJn6BurU/seJ9hnqRL5KFcFO4mRfDwQAhIU4/izuotw9HO9NnLk97d4/m2HM65Px+1/NAViG3edMfQPdOu1won3fCnSAoU7FqKTqlsLrg93XqnVXuHEzAABQocItD/dEPZwJ9H/+rYbXZ4zFZ1/3AAAEVrqOsS8twstDl8FgEDfsbumDb4U6A52K4xy7hdcHu6+ROgxflEYtf4s9yJlAzzIGYdZ7L2LBsgTk5hqg0ZjRt9dXmDZ+NqpGXXSiD74V6ABDnagsXh3srNado/fPR16eAbm5BqfP4cvD8c4+De7WLT2WrOqHme+9iCtXLU90euC+X/HW68lo0fSwk31hqBPZcPEcAC8PdnJOSLAR13MCcc0Y7OmueBVnA91k0mL1F49j2uxR+PdcNQBA/bon8dakZHT9zw7RN5mx9MX3Ah1gqFPZOBRv4bXBzmrdedGRl5B+vhrOnY9CXLODTp/HF6p2Kc9qN5s1WP9tV0ybPQrHT1puCnRX1fOYOOZd9Ov9Jfz8TE70h4FORGXz2mAn59WukYp9+5vi1Nmaks+l1nCXEuiCAHy75T+YOns0Dh65BwAQFnoNr41YhOEDP0KFCrlO9omhTlQmDsUD8NJgZ7UuTcMGfwP/A1IONXLJ+dQQ7lKC3Mps1mDj950w452RSDnUGAAQHJSNl4Yuw8ghKxASnO1k33wz0AGGOonDoXgLrwx2kubeFvsBAHv+aAlBgFNzvMV5a7i7ItBNJi02bOqC5HkjbBV6pYo5GD7wY4x+cQkqhzl3vwBL/xjqRCSO1wW7L1frrrrULf7efdDrc3E2tQaOnaiLe+4+6ZLzWkNSyQHviiC3ys3VY82XPTH3/aE4cToWABAUmI3hgz7CS0OWI7zKVSf76LthDjDQSQIOxQOAC7/LkdcIrHQDHe/fDQD44pvuLj+/tvA/JXFln7KMQXhn0fO4J/5HDH/1LZw4HYuw0Gv476j3cPy39piW9LbToe7rGOokieCCTQW8qmL35Wrd1fo8uR7fb++AFWueRtLLC+HvX+DyNjxRwcv5A8XZtLvw/rIErFjTG9nXgwAA1aIz8NLQZRjcdy2CAnOcPjerdAY6Scc5dguvCXaGumv1fPh7RIZfwrmMaHzy+RMY2Ocz2doqGrauDnm5RwYEAfhl771Y+GECNmzqDLNZBwC45+6/8dLQ5ejz5NdO3f61KF8PdSJyLa8JdnItgyEPYxKXYNzUCXjz7ZfRu+f/UKniTdnbVdoQfWmu51TEp189hiUr++HQ0Qa21x9s9zNeHrYMD3XYCa1W2o/3DHRW6uRinGMHwGD3acMGfIyFywYg9d/qmPDmOMybMcXTXfK4/YfuwYo1T2PNl4/DmG0Zbq8QcBPPPvk1EgevQqMGf0tug4FuwVAnV9MIAjSC8+ks5Vgl8YpgtwzD85uAqwUE5GHhrP/i0T4fYfHK/niw3S/o0W2rp7vldlnGIHz29aNYseZp/Hmgie31OrXPYmj/1ej/9OcICzVKboeBfhtDnUg+sgX79OnT8e233yIlJQV6vR7Xrl2TqymfoYXGpU93A4CHOvyMkUOWY/4Hg5Aw4h1s/vw5tG6536VtKFF+vh9++Kkd1nzZE//7/iHcumV5lK2/fx4e67oFCX0+w4PtfpE83A4w0ItioJOsOBQPQMZgz8vLQ69evRAfH49ly5bJ1Qy5wIyJb+HI8XrYtrMdejy3El9/PAj/1+ovT3fL5cxmDX7b1wLr1vfAl/97GBcvh9veu+fuvzGwz2fo8+R6l16qxlC/jaFOcuOqeAvZgn3q1KkAgJUrV8rVBLmIv38B1i0bjkf7rMTu31uha+/VWDBzAvr2Wu/prklmMmnxy95W+GZTF6z/tgvSz1ezvRcZfglP9diI53qtR8umB11yBz6AYV4cA53Ivbxijp3kF1jpBjZ+OgD9h7+Lb7d2wvMvv41tO+/HzMkzEBl+2dPdEyX7eiXs2NUWG7d0wqYfOtpV5kGB2Xi061Y8+8QGdLz/V6eesFYWhro9hjq5FYfiASgs2HNzc5Gbe/vJV0aj9AVL5LhKFW/i8xXDMOOdkZjxzkh8+uXj2LytI14fMw+DnluLgABp12vLxWTSYv/hhti+8z78sLMdft0bh7w8g+39sNBr6N75B/Ts9j06PbDL5b8PhvmdGOjkCRyKtxAV7ElJSZg5c2aZ+xw9ehQNGjQoc5/SJCcn24bwqWRyLKCzO79WwMQx76Fzx58wctyb2H+4EUa/PgWzFwzHsISPMeCZz1E16qJs7TuioECHA0fuwa97W2HX7jbYtac1rlwNs9snttZZdPvPDjzadSvua/2HLHfWAxjqJWGoE3mWRhAcv3Dv4sWLuHy57GHZ2NhY6PV629crV67EK6+84tCq+JIq9piYGFw6XgvBQfxmYSVnsBdVUKDD8tXPYOZ7w21z035++eh4/6/o2X0zOj2wCzWrn5O1D/n5fvj7VCxSDjVEysHG2Le/Cf460Bg3b1Ww2y8oMBvt439Dpwd2odMDP6Nu7BmXzZkXxzAvGQOdSmLMNiPs7tPIyspCcHCwPG0YjQgJCUHLZ6ZDpw9w+jymvFv4c+0EWfvqDqIq9oiICERERMjVFxgMBhgMhvJ39HFyV+1Wfn4mDB2wGgnPfobPNjyCDz9+Fnv+aIWtPz6ArT8+AACoVSMVrVumoEWTQ2hQ7xTq1j6L6tXOoUKF3HLObiEIQPb1QGRkRiDtXDWk/nsXzqTG4NSZWjh2oi5OnK5lN6xuFRqShTZxf+G+Nr+jffwexDU7KFtVbsVAJ1I2DsVbyDbHnpqaiitXriA1NRUmkwkpKSkAgLp16yIwMFCuZkkGen0++vZaj7691uPEqdr4cmM3bP6hI35PaYazqTVwNrUGPvu6h90xwUHZqBx2FcGB11Gxwk3o/EzQaAQU5PshN0+PnBsVkX09EFeuhSI3t+wf5oICs9H4nuNo0eQQWjQ9hNYtUlCvzhmXXGPuCAZ62Vipk2Jw8RwAGYN90qRJWLVqle3rFi1aAAB27NiBDh06yNUsyaxenTNIevl9JL38PrKMQdj7ZzPs298Uh440wLGTdXE2tTqu5wTCmB1kuyWrI4ICs3FX1QzE3HUOsTVTEVsrFfXrnkL9eidRs3q620K8KAZ6+RjqRI779ttvMW3aNBw4cAABAQF44IEH8PXXX9veT01NxfDhw7Fjxw4EBgZiwIABSE5Ohp+fuKiWLdhXrlzJa9hl5K7h+LKEBGfjoQ4/46EOP9teEwTgWlYwLl2ugktXwnA9pyJu3qwAk1kHs1kDPz8TDPo8VKqUg6BKOagcdg3hla+gYsVbHvyd3MYwdwwDnZRKqcPpX375JYYMGYIZM2bgwQcfREFBAQ4dOmR732QyoXv37oiOjsavv/6K8+fPo3///vD398eMGTNEtaWoy91IHCWEe3EaDRAWakRYqBH16pzxdHccxkB3HEOdFEsQLJuU42VQUFCAl19+GbNnz8bgwYNtrzds2ND26y1btuDIkSP44YcfEBUVhebNm+ONN97AuHHjMGXKFLtF6eXhv1DyaVpoGOoiMNTJFxiNRrut6NVazvjzzz+Rnp4OrVaLFi1aoGrVqujWrZtdxb579240adIEUVFRtte6dOkCo9GIw4cPi2qP/0q9HENJHGuQM9DF0Wm0DHVSPOuqeCkbAMTExCAkJMS2JScnS+rX6dOnAQBTpkzBxIkTsXHjRoSFhaFDhw64cuUKACAjI8Mu1AHYvs7IyBDVHv+lqgADqnwMcucx0MlrCC7YAKSlpSErK8u2jR8/vsTmkpKSoNFoytyOHTsGs9kMAJgwYQKefPJJxMXFYcWKFdBoNPj8889d/jFwjl0llDjf7mkMcucxzMmXBQcHO3SDmjFjxiAhIaHMfWJjY3H+/HkA9nPqBoMBsbGxSE1NBQBER0dj7969dsdmZmba3hODwa4iDHcLBjqRb9KYLZuU48Vw9KZtcXFxMBgMOH78OO6//34AQH5+Ps6ePYuaNWsCAOLj4zF9+nRcuHABkZGRAICtW7ciODjY7gcCRzDYVcYaar4U8Axy12GlTl5NoTeoCQ4OxgsvvIDJkycjJiYGNWvWxOzZswEAvXr1AgB07twZDRs2RL9+/TBr1ixkZGRg4sSJSExMFH1HVga7SvlC9c5Ady2GOpF8Zs+eDT8/P/Tr1w83b95EmzZtsH37doSFWR5gpdPpsHHjRgwfPhzx8fGoVKkSBgwYgGnTpolui8GuYmqq3hni8mGgk1oo+V7x/v7+mDNnDubMmVPqPjVr1sR3330nuS0Guw/wxoBnkLsHQ51URaE3qHE3BrsPKRqWSgt5Brl7MdBJjZRcsbsTg91HFQ9SdwU9A9yzGOhE6sdgJwClB67YwGdwKxdDnVRPoavi3Y3BTmViUHs/Bjr5Cg7FW/BfPJGKMdSJfA8rdiIVYqCTT+KqeAAMdiLVYaiTr+JQvAWDnUgFGOZEZMVgJyIideCqeAAMdiKvxSqdyB6H4i34nYHICzHUiag0rNiJvAgDnagMZsGySTleBRjsRF6AgU7kAM6xA2CwEykeQ53IMRpInGN3WU88i8FOpFAMdCJyBoOdSGEY6ERO4p3nADDYiRSDgU4kDS93s+B3EiIFYKgTkauwYifyIAY6kQtxVTwABjuRxzDUiVxLIwjQSJgnl3KskjDYidyIYU5EcmOwE7kJQ51IZubCTcrxKsBgJ5IRw5zIfTgUb8HvOkQyYagTkSewYidyMQY6kYdwVTwABjuRyzDQiTyMd54DwGAnkoyBTqQMvPOcBYOdyEkMdCJSIgY7kUgMdCKF4lA8AAY7kcMY6ETKpjFbNinHqwGDnagcDHQi8iYMdqJSMNCJvAyH4gHwBjVEd9BptAx1Im8kuGCTwY8//giNRlPi9vvvv9v2O3DgANq1a4eAgADExMRg1qxZTrXHip2oEMOciOTQtm1bnD9/3u61119/Hdu2bUOrVq0AAEajEZ07d0anTp2wePFiHDx4EIMGDUJoaCiGDh0qqj0GO/k8BjqROij1XvF6vR7R0dG2r/Pz87FhwwaMHDkSGo0GALB69Wrk5eVh+fLl0Ov1aNSoEVJSUjB37lzRwc7vaOSTrMPtDHUiFbHOsUvZYKmei265ubku7eY333yDy5cvY+DAgbbXdu/ejfbt20Ov19te69KlC44fP46rV6+KOj+/q5HPYZgTUVliYmIQEhJi25KTk116/mXLlqFLly6oXr267bWMjAxERUXZ7Wf9OiMjQ9T5ORRPPoFhTuQDBEh7pnrhSHxaWhqCg4NtLxsMhhJ3T0pKwsyZM8s85dGjR9GgQQPb1//++y++//57fPbZZxI6WjYGO6kaA53Id7hqjj04ONgu2EszZswYJCQklLlPbGys3dcrVqxAlSpV0KNHD7vXo6OjkZmZafea9eui8/OOYLCTKjHQiXyQAInXsYvbPSIiAhEREY6fXhCwYsUK9O/fH/7+/nbvxcfHY8KECcjPz7e9t3XrVtSvXx9hYWGi+sXvfqQaXBBHREq2fft2nDlzBs8///wd7/Xp0wd6vR6DBw/G4cOHsW7dOrz77rsYPXq06HZYsZPXY5ATEQDF33lu2bJlaNu2rd2cu1VISAi2bNmCxMRExMXFITw8HJMmTRJ9qRvAYCcvxkAnIjtmABqJx8tozZo1Zb7ftGlT7Nq1S3I7DHbyOgx0IqLSMdhJ8RjkROQIpd55zt0Y7KRoDHUicpjC59jdxSuCnd/cfQf/rImIpPGKYCf1Y6ATkWSs2AF4UbBbv/GbBJmXLZLbMMyJyKUY7AC8KNhJPRjoRETy8bpgZ+XufRjkROQWCr+O3V28LtitdBotw13hGOhE5E683M3Ca4MdYLgrEcOciDyGc+wAvDzYAQ7NexqDnIhIWWT7rnz27FkMHjwYtWvXRoUKFVCnTh1MnjwZeXl5srTHgHEfPkWNiBTJLEjfVEC2iv3YsWMwm81YsmQJ6tati0OHDmHIkCHIycnBnDlzZGmT1bvrMbyJyGtwKB6AjMHetWtXdO3a1fZ1bGwsjh8/jkWLFskW7FYMeOkY6ERE3smtc+xZWVmoXLlyqe/n5uYiNzfX9rXRaJTUHgPeMQxxIlIHiRU71FGxu+07+smTJzF//nwMGzas1H2Sk5MREhJi22JiYlzSNueD78R5ciJSHetQvJRNBUR/V09KSoJGoylzO3bsmN0x6enp6Nq1K3r16oUhQ4aUeu7x48cjKyvLtqWlpYn/HZWhaJj5SqAV/z370u+diMgXiR6KHzNmDBISEsrcJzY21vbrc+fOoWPHjmjbti2WLl1a5nEGgwEGg0Fsl5xWNODUMlzP0CYin2UWIGk43VdXxUdERCAiIsKhfdPT09GxY0fExcVhxYoV0GqVGzolBaISw57BTURUCsFs2aQcrwKyLZ5LT09Hhw4dULNmTcyZMwcXL160vRcdHS1Xsy5VWojKGfgMbiIikkK2YN+6dStOnjyJkydPonr16nbvCV6+QIHhS0SkQLyOHYCMq+ITEhIgCEKJGxERkcvxznMAVHCveCIiIgCs2AtxTJmIiEhFWLETEZE6CJBYsbusJx7FYCciInXgUDwADsUTERGpCit2IiJSB7MZgIT7jJh5gxoiIiLl4FA8AA7FExERqQordiIiUgdW7AAY7EREpBZ8uhsADsUTERGpCoOdiIhUQRDMkje5/P3333jssccQHh6O4OBg3H///dixY4fdPqmpqejevTsqVqyIyMhIvPbaaygoKBDdFoOdiIjUQZD4ABgZ59gfeeQRFBQUYPv27di3bx+aNWuGRx55BBkZGQAAk8mE7t27Iy8vD7/++itWrVqFlStXYtKkSaLbYrATEZE6WBfPSdlkcOnSJZw4cQJJSUlo2rQp6tWrh7feegs3btzAoUOHAABbtmzBkSNH8Mknn6B58+bo1q0b3njjDSxcuBB5eXmi2mOwExERFWE0Gu223NxcSeerUqUK6tevj48++gg5OTkoKCjAkiVLEBkZibi4OADA7t270aRJE0RFRdmO69KlC4xGIw4fPiyqPQY7ERGpg9ksfQMQExODkJAQ25acnCypWxqNBj/88AP++usvBAUFISAgAHPnzsXmzZsRFhYGAMjIyLALdQC2r63D9Y5isBMRkTq4aCg+LS0NWVlZtm38+PElNpeUlASNRlPmduzYMQiCgMTERERGRmLXrl3Yu3cvevbsiUcffRTnz593+cfA69iJiIiKCA4ORnBwcLn7jRkzBgkJCWXuExsbi+3bt2Pjxo24evWq7bzvv/8+tm7dilWrViEpKQnR0dHYu3ev3bGZmZkAgOjoaFH9Z7ATEZEqCGYzBI3zl6yJvdwtIiICERER5e5348YNAIBWaz9IrtVqYS4c/o+Pj8f06dNx4cIFREZGAgC2bt2K4OBgNGzYUFS/OBRPRETqoNBV8fHx8QgLC8OAAQOwf/9+/P3333jttddw5swZdO/eHQDQuXNnNGzYEP369cP+/fvx/fffY+LEiUhMTITBYBDVHoOdiIhIRuHh4di8eTOuX7+OBx98EK1atcLPP/+MDRs2oFmzZgAAnU6HjRs3QqfTIT4+Hn379kX//v0xbdo00e1xKJ6IiNTBLAAaZT4EplWrVvj+++/L3KdmzZr47rvvJLfFYCciInUQBAASbgurkqe7cSieiIhIRVixExGRKghmAYKEoXhBJRU7g52IiNRBMEPaULx8T3dzJwY7ERGpAit2C86xExERqYiiK3brT0/G6+oYHiEi8jXW79/uqIYLhFxJw+kFyHdhbzxH0cGenZ0NAKjZ8qxnO0JERJJkZ2cjJCRElnPr9XpER0fj5wzp14BHR0dDr9e7oFeeoxEUPKlgNptx7tw5BAUFQaPReLo7MBqNiImJQVpamkMPCPBl/Kwcw8/JMfycHKe0z0oQBGRnZ6NatWp33CvdlW7duoW8vDzJ59Hr9QgICHBBjzxH0RW7VqtF9erVPd2NOzj65B/iZ+Uofk6O4efkOCV9VnJV6kUFBAR4fSC7ChfPERERqQiDnYiISEUY7CIYDAZMnjxZ9CP0fBE/K8fwc3IMPyfH8bMiRS+eIyIiInFYsRMREakIg52IiEhFGOxEREQqwmAnIiJSEQa7k86ePYvBgwejdu3aqFChAurUqYPJkye75M5HajN9+nS0bdsWFStWRGhoqKe7oygLFy5ErVq1EBAQgDZt2mDv3r2e7pKi7Ny5E48++iiqVasGjUaDr7/+2tNdUqTk5GTce++9CAoKQmRkJHr27Injx497ulvkIQx2Jx07dgxmsxlLlizB4cOH8c4772Dx4sX473//6+muKU5eXh569eqF4cOHe7orirJu3TqMHj0akydPxp9//olmzZqhS5cuuHDhgqe7phg5OTlo1qwZFi5c6OmuKNpPP/2ExMRE7NmzB1u3bkV+fj46d+6MnJwcT3eNPICXu7nQ7NmzsWjRIpw+fdrTXVGklStX4pVXXsG1a9c83RVFaNOmDe69914sWLAAgOXZCDExMRg5ciSSkpI83Dvl0Wg0WL9+PXr27OnprijexYsXERkZiZ9++gnt27f3dHfIzVixu1BWVhYqV67s6W6QF8jLy8O+ffvQqVMn22tarRadOnXC7t27PdgzUoOsrCwA4PcjH8Vgd5GTJ09i/vz5GDZsmKe7Ql7g0qVLMJlMiIqKsns9KioKGRkZHuoVqYHZbMYrr7yC++67D40bN/Z0d8gDGOzFJCUlQaPRlLkdO3bM7pj09HR07doVvXr1wpAhQzzUc/dy5nMiIvklJibi0KFDWLt2rae7Qh6i6Me2esKYMWOQkJBQ5j6xsbG2X587dw4dO3ZE27ZtsXTpUpl7pxxiPyeyFx4eDp1Oh8zMTLvXMzMzER0d7aFekbcbMWIENm7ciJ07dyrykdfkHgz2YiIiIhAREeHQvunp6ejYsSPi4uKwYsUKaLW+MwAi5nOiO+n1esTFxWHbtm22xWBmsxnbtm3DiBEjPNs58jqCIGDkyJFYv349fvzxR9SuXdvTXSIPYrA7KT09HR06dEDNmjUxZ84cXLx40fYeKy57qampuHLlClJTU2EymZCSkgIAqFu3LgIDAz3bOQ8aPXo0BgwYgFatWqF169aYN28ecnJyMHDgQE93TTGuX7+OkydP2r4+c+YMUlJSULlyZdSoUcODPVOWxMRErFmzBhs2bEBQUJBtnUZISAgqVKjg4d6R2wnklBUrVggAStzI3oABA0r8nHbs2OHprnnc/PnzhRo1agh6vV5o3bq1sGfPHk93SVF27NhR4t+dAQMGeLprilLa96IVK1Z4umvkAbyOnYiISEV8Z1KYiIjIBzDYiYiIVITBTkREpCIMdiIiIhVhsBMREakIg52IiEhFGOxEREQqwmAnIiJSEQY7ERGRijDYiYiIVITBTkREpCIMdiIiIhX5f3UTZ1WRRjGeAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Conveniently caustic has a function to compute the jacobian of the lens equation\n", "A = sie.jacobian_lens_equation(thx, thy, z_s, packparams)\n", "# Note that if this is too slow you can set `method = \"finitediff\"` to run a faster version. You will also need to provide `pixelscale` then\n", "\n", "# Here we compute A's determinant at every point\n", "detA = torch.linalg.det(A)\n", "\n", "# Plot the critical line\n", "im = plt.imshow(detA, extent = (thx[0][0], thx[0][-1], thy[0][0], thy[-1][0]), origin = \"lower\")\n", "plt.colorbar(im)\n", "CS = plt.contour(thx, thy, detA, levels = [0.], colors = \"b\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3b099cfd", "metadata": {}, "source": [ "## Caustics\n", "\n", "Critical lines show us where the magnification approaches infinity, they are important structures in understanding a lensing system. These lines are also very useful when mapped into the source plane. When the critical lines are raytraced back to the source plane they are called caustics (see what we did there?). In the source plane these lines deliniate when a source will be multiply imaged. " ] }, { "cell_type": "code", "execution_count": 4, "id": "0c1f1177", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABeHklEQVR4nO3dd3zTdf4H8FdGk3Sle09ogTJbaKGCAiKoOBBO3ANBRO9OvIHendx56p13h9458E5+eg7cioqLU0TZs6yyoYPuvdskTZr9/f2RtlpldCT9Junr+Xjk0fbbb9o3X9Lklc+UCIIggIiIiMhDSMUugIiIiKgvGF6IiIjIozC8EBERkUdheCEiIiKPwvBCREREHoXhhYiIiDwKwwsRERF5FIYXIiIi8ihysQtwNrvdjpqaGgQGBkIikYhdDhEREfWCIAjQ6XSIjY2FVHrhthWvCy81NTVISEgQuwwiIiLqh8rKSsTHx1/wHK8LL4GBgQAc/3i1Wi1yNURERNQbWq0WCQkJ3a/jF+J14aWrq0itVjO8EBEReZjeDPnggF0iIiLyKAwvRERE5FEYXoiIiMijMLwQERGRR2F4ISIiIo/C8EJEREQeheGFiIiIPArDCxEREXkUhhciIiLyKAwvRERE5FEYXoiIiMijMLwQERGRR2F4ISIiol4pb9bjue8KUNSgE7UOr9tVmoiIiJxv7Z5S/PWrMwCAWo0Rz96cLlotDC9ERER0XlqjBVc+vxP1WlP3sVuyEkSsiOGFiIiIzmN7fgOWvHWox7E9f5iF+BA/kSpyYHghIiKiHtoMZvz5y9P43/Ga7mMjIgPw9a+mQyEXf7gswwsREREBAARBwOdHq/GHT0/AYhO6jz923WjcN324iJX1xPBCREREKGrQ4bEvTmF/SUuP4+8unYLpIyJEqurcGF6IiIiGsA6zDS9tP4tXd5X0aG0J81fg3aXZGBOrFrG6c2N4ISIiGqK25dfjiQ2nUdnS0eP4qKhAvLlkMmKDfUWq7MIYXoiIiIaYogYdnvoqDzsLG3/yvWkpYXj5rkwE+fqIUFnvMLwQERENERqDBS9uPYt3cspgtQvwkUkgk0pgtNgBAAsyYvHPm9LdYkbRhTC8EBEReTmbXcCHByvw3HcFaDVYAABzRkdCrfLBZ0erAQCLpyXj8evHQCqViFlqrzC8EBERebG9RU146qszyK9z7Ec0IjIAf75+DHYUNGLt3lIAwG/mjMCvZ4+AROL+wQVgeCEiIvJKZ2q0eHpTPnZ1jmsJ8vXBiitH4tbJCXjsi1NYn1sFAHj8+jG497JhYpbaZwwvREREXqSq1YDnNxfi86PVEATARybBndlJ+PXsEfBXyvHQh0fw7el6yKQS/HPhBCzMjBe75D5jeCEiIvICbQYz1mwvwts55TBbHQNwr58Qg99dPQpJYf4wWW34xXu52JrfAIVMiv/cMRFXj40Wuer+YXghIiLyYEaLDW/vK8Oa7UXQGq0AgKnDw/DoNWlITwjuPufn7+ViR0EjlHIpXluUhRkj3WvV3L5geCEiIvJAJqsNHx2qxJrtRajXmgAAadGB+MM1abh8ZET34NsOsw33v3sYu882QeUjxRv3TMalqeFilj5gDC9EREQexGKzY31uFV7aVoTqNsfKuHHBvvjtlSPxs4lxkP1gqrPBbMV9bx/GvuJm+ClkWLt4Mi4ZHiZW6U4zKKvQrFmzBsnJyVCpVMjOzsbBgwd7db9169ZBIpFgwYIFri2QiIjIzVk7Q8vs53Zi5WcnUd3WgSi1Ek/NH4ttj8zETZnxPYKL3mTF4jcPYV9xM/wVMrx97xSvCC7AILS8fPTRR1ixYgVeeeUVZGdnY/Xq1bj66qtRUFCAyMjI896vrKwMjzzyCKZPn+7qEomIiNyWzS7gqxM1eHHLWZQ06QEA4QEK/OLyVNyZnQiVj+wn92k3WbF47UEcLm9FgFKOt++dgsykkMEu3WUkgiAIFz+t/7KzszF58mS89NJLAAC73Y6EhAQ89NBDePTRR895H5vNhhkzZuDee+/F7t270dbWhi+++KJXv0+r1SIoKAgajQZqtfvthElERNQbXaHlpW1FONvQDgAI8fPBz2em4O6pSfBTnLv9QW+y4p7O4BKokuPdpdnI6By468768vrt0pYXs9mM3NxcrFy5svuYVCrFnDlzkJOTc977/fWvf0VkZCSWLl2K3bt3X/B3mEwmmEym7q+1Wu3ACyciIhKJ2WrH50er8PKOYpQ1GwAAapUc988YjsWXDkOA8vwv3UaLDfe9fRiHy1uhVsnx3n3ZmBAfPEiVDx6XhpempibYbDZERUX1OB4VFYX8/Pxz3mfPnj144403cOzYsV79jlWrVuEvf/nLQEslIiISldHimD30353FqNEYAThaWu69dBgWTUu+6C7PJqsN97+bi5yS5u6uIm8MLoCbzTbS6XS4++678dprryE8vHfTuFauXIkVK1Z0f63VapGQkOCqEomIiJyq3WTFe/vL8fruUjS1O3oSIgKVuH/6cNyRnQj/C7S0dDFb7Xjw/SPYVdgIXx/HrKKJid4zxuXHXBpewsPDIZPJUF9f3+N4fX09oqN/uqpfcXExysrKMG/evO5jdrtjlUC5XI6CggKkpKT0uI9SqYRSqXRB9URERK7T3G7COznleGtfGTQdjp2e44J98fPLU3BzZvw5B+Kei9Vmx28+OooteQ1QyKV4/Z4sTBkW6srSRefS8KJQKJCZmYmtW7d2T3e22+3YunUrli9f/pPz09LScPLkyR7HHnvsMeh0Orz44otsUSEiIo9X0tiON/aUYn1uFUydy/gPD/fHLy5PwYKJcfCR9X4VE5tdwCOfHMfGk3XwkUnw37szPX4But5webfRihUrcM899yArKwtTpkzB6tWrodfrsWTJEgDAokWLEBcXh1WrVkGlUmHcuHE97h8cHAwAPzlORETkSQ6XteDVXSXYnFePrnm+E+KDcP+M4bhmXEyPNVp6w24X8MfPTuKLYzWQSSV46Y5JmDXq/EuQeBOXh5dbb70VjY2NePzxx1FXV4eMjAxs2rSpexBvRUUFpNJBWSuPiIhoUNnsAjafqcOru0pwpKKt+/jstEgsmzEc2cNCu5fx7wtBEPDU12fw0eFKSCXA6lszPHaTxf5w+Tovg43rvBARkdjaTVZ8mluFtXtLUd453Vkhk+JnE+OwbMYwpEYGDujnv7TtLJ79rhAA8OzN6bgpM37ANYvNbdZ5ISIiGkpKm/R4e18Z1udWod3k2OE5yNcHd1+ShEXTkhAZqBrw73j/QHl3cHn8+jFeEVz6iuGFiIhoAOx2ATvPNuKtvWXYWdjYfXx4hD/umZqMmzLjezXduTe+PlGLx744BQBYPisV9142zCk/19MwvBAREfWD1mjB+sNVeHd/OUo79xySSIBZoyJxz7RkTE8Nh7SPg3AvZPfZRvzmo6MQBOCO7EQ8fNVIp/1sT8PwQkRE1AeF9Tq8t78cn+ZWQW+2AQAClXLcnJWARVOTkBzu7/TfeayyDQ+8mwuLTcB142Pw1Pxx/Rro6y0YXoiIiC7CaLFh48lafHCgAofLW7uPp0YG4J5pybhxYpzTuoZ+rKhBhyVvHoTBbMP0EeF4/tb0Pk+r9jYML0T9UNTQjha92etXsSQa6s7W6/D+gQp8dqQKWqNjAK5MKsHstEgsmpqMS1PDXNoCUtPWgbvfOIhWgwXpCcF45a5MKOW9W3nXmzG8EPXDvW8dQkWLAfMzYvHibRPFLoeInKirleXDgxU4VPZ9K0tcsC9un5KAm7MSEKUe+Kyhi9F0WLD4zYOo1RiRGhmAtxZPdlnrjqfhVSDqh4oWx7oNXx6rwbwJsZgzJuoi9yAid3emRouPD1fi86PV3XsNdbWy3JGdiOkjIgatu8ZkteGBdw+jsL4dUWol3r53CkL8FYPyuz0BF6kj6odWvRkTn9rc/fUrd03C3HExIlZERP3R3G7Cl8dqsD63Cmdqtd3HB7uV5YfsdgG//fgYvjxWgwClHB8/MBVjYr3/9YyL1BG5WIi/Ag/MHI7/7iwBADz4wVGsutGKW7K4eSiRu7PY7Nie34D1uVXYlt8Aq93xHl4hk+LKMVG4KSseMwaxleXH/vVdAb48VgO5VIKX75o0JIJLXzG8EPXTAzNS8MbuUljtAmx2Ab9ffwLFDe34/dy0IT8TgMjdCIKAvFod1udW4ctj1WjWm7u/lx4fhJsy4zEvPRbBfuJ2zby7vxwv7ygGADy9cAKmj4gQtR53xfBC1E+h/grckB6Lz45Wdx/7764SFNTr8OJtExHk6yNidUQEAJUtBmw4XoMvj1WjsL69+3hEoBI3TozDwsx4jIwa2D5DzrL5TD2e+NKxeu6KK0cOyWX/e4vhhWgA7pqahM+OVkMpl+LJG8biL/87jR0Fjbj2xd148bYMZCVzKjXRYGvUmfD1iRp8ebwGR3+wk7NCJsWcMZG4OTMB00eEQy6TilfkjxyrbMNDHx6BXQBum5yAh65IFbskt8bwQjQAExOCkRLhj+JGPZRyKdb/fBp++f4RVLQYcMt/c/Cr2SOwfFaqWz1JEnkjndGCb0/X48tj1dhb1ITOYSyQSoBpKeG4ISMWV4+NdssW0eq2Dtz39mEYLXbMHBmBpxYM7dVze4OzjYgGaPWWQqzechazRkXgzSVToDNa8MSXp7u7k9ITgvHMwvFIi+bjkciZtEYLtuU14JtTtdhe0Aiz1d79vfSEYMxPj8X1E2IQOcizhfqi3WTFTS/vQ36dDmnRgVj/i2kIGKJrufTl9ZvhhWiAihraMef5nZBLJch97EoE+Tne2X1xtBp//uIUdCYr5FIJfj4zBcuvSIXKh6tjEvWXxmDB5rx6fHOyFrvPNsFs+z6wDI/wx4KMONyQHuuS/YWczWYXcP87h7E1vwERgUp8+eCliA32Fbss0XCqNNEgSo0MwIjIAJxtaMfuokZcPyEWALBgYhwuGR6Gx788he/O1OOl7UXYeLIWT94wFjNGcgYBUW81t5uw+Uw9Np6qw76ipu6pzYAjsFw7LgZzx0VjbKzao7pbVm3Mw9b8BijlUry2KGtIB5e+YnghcoLLR0XgbEM7tud/H14AIDpIhVcXZWHTqVo8/uVplDTpsWjtQVw+KgJ/unY0RrjJLAcid1PVasDWvAZ8e7oO+0ua8YO8glFRgbhmfDSuHR+DEZEBHhVYunxwoAKv7ykFADx3SzoyEoLFLcjDMLwQOcHloyLx2u5S7CxshN0uQPqjdV7mjovBtNRwvLjlLN7JKcOOgkbsPtuEO6Yk4lezRyAiUClS5UTuwW4XcLyqDVvzGrAlrx75dboe3x8Xp8Y1nS0sKREBIlXpHPuKmvD4D6ZE//AND/UOx7wQOYHJakPGXzajw2LDt7+ZgVHR529RKW3SY9XGPHx3ph4AoPKRYtHUZNw/YzjCAxhiaOjoMNuwp6gJW87UY2t+A5raTd3fk0qArKRQzBkTibljY5AY5idipc5T0tiOBWv2Qmu0Yn5GLFbfmuGRLUeuwDEvRINMKZchIyEYOSXNOFzecsHwMizcH68uykJOcTOe2ZSPY5VteHVXCd7NKceiaUm4f/pwhDHEkJcqa9Jj19lG7ChoxN6iJph+MEMoQCnHzJERmDMmEpePjPS6jQg1BguWvn0YWqMVkxKD8czCCQwu/cTwQuQkWckhyClpRm5ZK+7MTrro+VNTwvD5L6dhR2EjVm8uxPEqDf67swRv7S3DjZPisfSyYUiN9OzmcSK9yYqc4mbsLGzErrONKG829Ph+fIgv5oyOwpzRUZgyLBQKuXeuiWSzC/jVuqMobdIjLtgXry7K4szDAWB4IXKSzKQQAMDh8tZe30cikWDWqEhcPjIC2wsa8OLWIhyvbMOHByvw4cEKXJEWifumD8PU4WF8h0YewW4XkFenxa7CJuwqbMTh8hZYbN+PTvCRSZCZFIIZIyNwRVokRkUFDonH9r++LcDOwkaofKR4dVEmu4gHiOGFyEkmJYVAIgEqWgxo1Jn6NAhXIpHgirQozBoViUNlrXhtdwm25NVjW34DtuU3IDUyALdPScTCSXGibxxH9EOCIKC4UY+c4ibsK27G/pJmtBosPc5JDPXDzJERmDEyAlNTwobcImwbjtfglZ2OzRb/dVM6xsYGiVyR5xtajyAiF1KrfDAs3B8ljXrk1WoREdj3tVwkEgmmDAvFlGGhKG3SY+2eUqzPrUJRQzue+uoMntmUj2vHReP2KYmYMix0SLxjJfdT2WLAvuIm5BQ3Y19xMxp0ph7f91fIkD08DDNHRmDmyAiPWDDOVU5Va/D79ccBAD+fmYJ56ZxZ5AwML0ROlBYdiJJGPQrqdANeiG5YuD+eWjAOv587Cl8eq8EHBypwplaLL47V4ItjNUgI9cUN6bGYnxHnNrvikvex2wUUNuhwuKwVh8tacKisFdVtHT3OUcilyEoKwbSUMExNCceE+CD4cD8vNLeb8MC7uTBa7Lh8VAR+d/UosUvyGgwvRE6UFq3GxpN1yKvTOu1nBqp8cNclSbgzOxEnqzX48GAFNhyrQWVLB9ZsL8aa7cVIiw7EDRmxmDchFgmh3jGllMRhtNhwokqDQ2UtOFzWgtzyVmiN1h7nyKUSZCQEd4eViYnBHHz6IxabHb98/wiq2zowLNwfL942ETIpW0qdheGFyIm6pkgX/GiBLWeQSCSYEB+MCfHBePz6sdiSV48Nx2uwo6AB+XU65G8qwD83FWB0jBpXjo7ElWOiMS7Os5ZLp8ElCAIqWzpwvKoNJ6racKSiDSerND32CwIAP4UMExODkZUUisnJoZiYGAz/ITZupa/+9tUZHChtQYBSjtcWZbrlbtaejI8+IidK6wwvZxvaz7nSrrP4KmSYlx6Leemx0Bgs+OZULb48VoMDpc3Iq9Uir1aLf28rQrRahdmjIzF7dCSyh4XxBWeIa9SZcKKqDcerNDhe6QgsPx5cCwARgUpMTg7pDiujYwIhZzdQr32aW4W3c8oBAM/fko7USHbrOhufyYicKC7YF3KpBGarHfU6I2KCXL/RWpCfD26bkojbpiSiRW/G9nzH8uo7CxtRpzXi/QMVeP9ABeRSCSYmBmNaSjguTQ1HRkKw166pMdQJgoB6rQlnajXIq9XhVLUGJ6o0PxmrAjimLo+OUSM9PhjpCcGYnByCxFA/ttj1U16tFn/64iQA4NezR+CqsdEiV+SdGF6InEgukyIuxBflzQZUNBsGJbz8UKi/Agsz47EwMx5Giw05Jc3YcqYeu842orKlA4fKWnGorBUvbj0LP4UMU4Z93w2QHs+uAE9kttpR1NDe3eJ2pvPjuVpUJBIgNSIAE+KDkZ4QhPT4YKTFBEIp53gVZ9B0WPDz9xwDdGeOjMCvZ48QuySvxWcqIidLDPVDebMB5S0GZA8PE60OlY8Ms0ZFYtaoSABARbMBe4ubsLfIMcW1WW/GjgLHMu2AYy+ZtGg1JiUFY1JiCCYmhiAp1M9lXV/UN0aLDSWNepxt0KG4oR1Fje0oamhHaZO+xyJwXWRSCVIi/DE6Ro0xMWpMiA/GuDg1AlUce+EKgiDgkU+Oo7zZgLhgX6y+NYN/Oy7E8ELkZImds30qWwwXOXNwJYb5ITEsEbdPSYTdLiC/ToeckmYcqWjF0fJW1GiMONP5zv29/RUAHOt1jIoOdLwAxqoxOkaNtOhA+Cn41OEKNruAmrYOVLQYUN5sQFmzHkUNjpBS2WrA+bbRDVTJu0PKmBjH/9OIqADOABpE/91Vgs1n6qGQSfHyXZO8bl8md8NnICIniwlSAQAatKaLnCkeqVSCMbGOQLIUwwAAdRojjlS04kh5K3IrWnG6Rgu92YYjFY5ZKF0kEiAp1A8pEQEYFu6PYRH+GBbuj+HhAYhSKzlW4gIEQUCL3oxajbE7pHQFlYoWA6paDedsRekS7OeD1IgApEb2vMUF+/K6iyinuBn/3JQPAHjyhrGYEB8sbkFDAMMLkZN1bQvQ2O6+4eVcooNUuHZ8DK4dHwMAsNrsKG3Sd7fG5NXqcKZGi6Z2E8qaDShr/mnLkq+PDMPC/REf4ovYYF/EBKkQE+yL2CAVYoN9ERmo9NpZK0aLDc16M5rbTWjQmlCr6UCNxoi6zqBSpzWiVmOE2Wq/4M9RyKSID/VFcpg/EkP9kBIZgBGdISXMX8GQ4mbqtUY89OER2AVg4aR43D4lQeyShgSGFyIn69pwrVHnWeHlx+QyKUZEBWJEVCDmZ8R1H2/UmVBQp0Npsx6ljXqUNjnGXVS2dqDDYusOO+cilQCRgSqEBSgQ6v+Dm58CoQGOjyH+CgSq5AhQyuGvdHxUyqWD8qItCAIMZhu0Rgt0Rit0Rgu0Rmv35zqjFdoOC1r0ZjS1m9GsN6G53YwWvRntJuvFf0Gn8AAlYoJUSAj1RVKYP5JC/ZAY5oekMH9Eq1VczMxDWGx2PPj+ETS1m5EWHYi/LRjHcDlIGF6InKy75cXDw8v5RAQqERGoxGUjwnsct9jsqGwxoLRJj+q2DtS0GVGr6UBtmxE1mg7UaYyw2gXUaY2o0xr79DvlUgn8lXL4+sigkEvhI5PARyaFUi6Fj8xxU8il6HrdEARAgCOMfP+1AJtdgNlqh6nzZrTYHJ93fbxIq8jF+MgkCPN3XJ+YIFV3y1NMkArR6s7WJ7WSs3u8xL++LcDh8lYEKuV45a5M+Cr4/zpYGF6InKwrvDS1m1y6UJ278ZFJMTwiAMMjAs75fbtdQFO7CbUaI1r0ZjTrzWj90ccWvQmtBgvaTVboTVYYzDYAgNUuQNNhgabjp9N/XUEulSBQJUegyqfz4/efq1U+CPNXICxAiVB/BcIDHJ+HBSgQqJTznfcQsS2/Hq/uKgEA/Ovm9CG9+aQYGF6InCy0c5aB1S5AZ7JyWfBOUqkEkWoVItWqXt/HZhegNzuCjN5khdHiaB2x2Bw3c+fnJqvjcwFAV3SQSCSQwDHAWCIBJJBAKpVAKXe02Kh8ZJ2fy6D0kXZ/HqCUQ+UzON1U5JlqNR14+GPHTtGLpyVj7jguRDfYGF6InEwpl0Ehk8Jss0PP8DIgMqkEapUP1FybhNyE1WbHrz88hlaDBePi1Fh5bZrYJQ1J3jnsn0hkfkpH37fB3PtBnETk/v699SwOljk2XHzp9kkcvyQShhciF/DvXMSt3WQTuRIicpa9RU34z/YiAMA/bhzPcS4iYnghcgH/rpaXPkyfJSL31agz4TcfHYMgALdNTsAN6bFilzSkMbwQuUDXBod9WfuDiNyT3S5gxcfH0KgzYWRUAJ6YN1bskoY8hhciF/Dt3FOmw8JuIyJP98quYuw+2wSVjxRr7pjE9VzcAMMLkQt0rZB6vo30iMgznKhqw/PfFQIA/nrDOIyIChS5IgIYXohcomuNEDvTC5HHMpit+M26Y7DaBVw7Pho3Z8WLXRJ1YnghcoGuRXVtdoYXIk/11Fd5KGnSI1qtwj9+Np4LF7oRhhciF5BK2G1E5Mm+O12HDw9WQCIBnr8lHcF+CrFLoh9geCFyga6WF3YbEXmeBp0Rj352EgCwbPpwTEsNv8g9aLAxvBC5AHuLiDyTIAj43Scn0KI3Y0yMGg9fNVLskugcGF6IXMBstQMAlD78EyPyJG/vK8POwkYo5VK8eFsGl/93U3xmJXKBrvCikPGJj8hTFNbr8I9v8gEAf7x2NKdFuzGGFyIXMNk6W17k/BMj8gRmqx2/WXcMZqsdl4+KwKKpSWKXRBfAZ1YiFzB1rqyrYHgh8ggvbS/CmVotQvx88M+bJnBatJvjMyuRC5jZ8kLkMU5WabCmc7fopxaMQ2SgSuSK6GL4zErkAiZL55gXhhcit2ay2vDwJ8dgswu4bkIMrp/A3aI9AZ9ZiVxAa7QAANS+PiJXQkQX8sLmsyisb0d4gAJPzR8ndjnUSwwvRE5mswvQGa0AALWK4YXIXR2paMWru4oBAP/42XiE+nMVXU/B8ELkZO2dwQUAgtjyQuSWjBYbHvnkOOwC8LOJcbhqbLTYJVEfMLwQOZmmw9Fl5Osj45gXIjf1r28LUNKoR5RaiSfnjRW7HOojPrMSOVlXeGGrC5F7OljagrV7SwEATy+cgCA//q16mkEJL2vWrEFycjJUKhWys7Nx8ODB85772muvYfr06QgJCUFISAjmzJlzwfOJ3A3DC5H7Mlps+P364xAE4NasBMwaFSl2SdQPLg8vH330EVasWIEnnngCR44cQXp6Oq6++mo0NDSc8/wdO3bg9ttvx/bt25GTk4OEhARcddVVqK6udnWpRE7RrDcBAEL8GV6I3M3qLWdR1mxAtFqFP10/WuxyqJ9cHl6ef/55LFu2DEuWLMGYMWPwyiuvwM/PD2vXrj3n+e+//z5++ctfIiMjA2lpaXj99ddht9uxdetWV5dK5BQNWkd4iVJzoSsid3KqWoPXdpcAAP62YBxnA3owl4YXs9mM3NxczJkz5/tfKJVizpw5yMnJ6dXPMBgMsFgsCA0NPef3TSYTtFptjxuRmOq1RgAML0TuxGqz4w+fnoDNLuD6CTGYMyZK7JJoAFwaXpqammCz2RAV1fNBEhUVhbq6ul79jD/84Q+IjY3tEYB+aNWqVQgKCuq+JSQkDLhuooGo1zlaXiIDlSJXQkRdXt9TitM1WgT5+uAJzi7yeG492+jpp5/GunXr8Pnnn0OlOve72JUrV0Kj0XTfKisrB7lKop4aOlteItnyQuQWSpv0eGFzIQDgz9ePQQTfWHg8uSt/eHh4OGQyGerr63scr6+vR3T0hRcEevbZZ/H0009jy5YtmDBhwnnPUyqVUCr5QCT30dDZ8hLFJ0gi0QmCgJWfnYDJasf0EeFYOClO7JLICVza8qJQKJCZmdljsG3X4NupU6ee937//Oc/8dRTT2HTpk3IyspyZYlETiUIAse8ELmRjw5VYn9JC3x9ZPjHz8ZDIpGIXRI5gUtbXgBgxYoVuOeee5CVlYUpU6Zg9erV0Ov1WLJkCQBg0aJFiIuLw6pVqwAAzzzzDB5//HF88MEHSE5O7h4bExAQgICAAFeXSzQgbQYLDGYbACA6iOGFSEz1WiP+vjEPAPDwVSOREOonckXkLC4PL7feeisaGxvx+OOPo66uDhkZGdi0aVP3IN6KigpIpd83AL388sswm8246aabevycJ554Ak8++aSryyUakPIWAwAgSq2EykcmcjVEQ9tf/ncaOqMV6fFBWHLpMLHLISdyeXgBgOXLl2P58uXn/N6OHTt6fF1WVub6gohcpLxZDwBICvUXuRKioW1HQQM2nqyDTCrBqhsnQCZld5E3cevZRkSepqLZ0fKSGMbmaSKxGC02PLHhNABg8bRkjIlVi1wRORvDC5ETdXUbJbJvnUg0/7ejGOXNBkSplfjtlSPFLodcgOGFyIkqOsNLElteiERR2qTHKzuKAQCPXz8WAcpBGR1Bg4zhhciJuruN2PJCNOgEQcDjX56C2WbHjJERuHb8hdcTI8/F8ELkJAazFXWda7wkh3HALtFg+/pkLXafbYJCLsVfbxjLNV28GMMLkZOcrW8HAIQHKBHirxC5GqKhRWe04K//OwMA+OXlKUgO5xsIb8bwQuQkhfU6AMDIKC6mSDTYXth8Fg06E5LD/PDzmSlil0MuxvBC5CRnGxwtLyOjAkWuhGhoKazX4e2cMgDAX+eP4wKRQwDDC5GTdLW8jGDLC9GgEQQBf/3fGdjsAq4aE4UZIyPELokGAcMLkZMU1nV1G7HlhWiwfHemHnuKHIN0H7tujNjl0CBheCFyAp3RghqNY6bRyEiGF6LBYLTY8LevHYN0l00fxpWthxCGFyInKOhsdYlWqxDk5yNyNURDwxt7SlHZ0oEotRK/vDxV7HJoEDG8EDnBiSoNAGBcHPdQIRoMdRoj1mwvAgCsvGY0/LmS7pDC8ELkBKeqHeFlfFywuIUQDRFPf5MHg9mGzKQQzM+IFbscGmQML0ROcLIrvMSz5YXI1XLLW/DFsRpIJMCT87iS7lDE8EI0QHqTFcWNjjVexsUFiVwNkXez2wU8ucExSPeWzASMj+ff3FDE8EI0QGdqtbALQJRaichAldjlEHm1L45V42S1BoFKOX43d5TY5ZBIGF6IBuhkVdd4F74DJHIlo8WGZ78tAAD8YlYKwgOUIldEYmF4IRqgrvEu7DIicq0395ahRmNEbJAK9146TOxySEQML0QDdLi8BQAwKTFE5EqIvFdzuwn/1zk1+pGrR3H/oiGO4YVoABq0RlS2dEAqASYmBotdDpHX+s+2IuhMVoyNVWNBRpzY5ZDIGF6IBuBweSsAYFS0GoEqrqxL5AqlTXq8t78cAPDHa0dDKuXU6KGO4YVoAA6XOcJLVhK7jIhc5Zlv8mG1C7h8VAQuTQ0XuxxyAwwvRAOQ2zneJSuZ4YXIFQ6XtWDT6TpIJY5tAIgAhheifusw23C6RgsAyGTLC5HTCYKAp7/JBwDckpWAUdHcsZ0cGF6I+uloZSusdgHRahXign3FLofI6+woaMTh8lYo5VL8Zs5IscshN8LwQtRP+4ubAQDZw0O5twqRk9ntAv7ZuSDd4mnJiA7i6tX0PYYXon7a2xlepqWEiVwJkff5+mQt8mq1CFDK8fOZKWKXQ26G4YWoH9pNVhyvbAMATEvh7AciZ7La7Hh+cyEAYNn04QjxV4hcEbkbhheifjhU2gKrXUBiqB8SQv3ELofIq3x6pAqlTXqE+iuwdDq3AaCfYngh6oe9RU0A2GVE5GxGiw0vbjkLAPjl5SkIUMpFrojcEcMLUT/s6xrvwgWziJzqgwMVqNEYEa1W4a5LksQuh9wUwwtRH7XozThT61jfZepwtrwQOYveZMWazs0XfzV7BDdfpPNieCHqo33Fji6jUVGBiAhUilwNkfd4b385mvVmJIX54easeLHLITfG8ELUR9vyGwAAM0dFiFwJkfcwmK14dVcJAODBWanwkfHlic6Pjw6iPrDbBewsaAQAzBoVKXI1RN7jgwMVaNabkRDqi59NjBO7HHJzDC9EfXCiWoNmvRmBSjk3YyRykg6zDa/sdLS6LGerC/UCHyFEfdDVZTR9ZDifYImc5IODFWhqNyEu2Bc/m8ixLnRxfPYl6oMdBY7wcjm7jIicwmix4ZWdxQAcY10Ucr4s0cXxUULUSw06I05UaQAAl3OwLpFTrDtYgUadCbFBKtyUyVYX6h2GF6Je6hqoOz4uCJGB3OGWaKCMFhte7mx1+QVbXagP+Egh6qXtnV1Gs9LYZUTkDOtzq1CvNSEmSIVbuK4L9QHDC1EvGC027OhseZnN8EI0YFabvXtdl/tnDIdSztV0qfcYXoh6YVdhIwxmG+KCfTEhPkjscog83sZTdahoMSDEzwe3Tk4QuxzyMAwvRL2w6VQdAODqsdGQSCQiV0Pk2QRBwCs7HGNdFk8bBj8Fd46mvmF4IboIs9WOLXn1AIC546JFrobI8+0624QztVr4+siwaCp3jqa+Y3ghuoj9Jc3QGq0ID1AiM4mr6hINVFery+1TEhHirxC5GvJEDC9EF/FNZ5fRVWOjIJOyy4hoII5VtiGnpBlyqQT3TR8mdjnkoRheiC7AZhew+YwjvFzDLiOiAetqdZmfEYfYYF+RqyFPxfBCdAGHy1rQ1G6GWiXHJcPDxC6HyKOVNunxbeebgZ/PHC5yNeTJGF6ILmDD8RoAwFVjo7kRI9EAvbm3FIIAXJEWiRFRgWKXQx6Mz8ZE52G22vH1yVoAwPyMWJGrIfJsGoMFnxyuAgAsvYxjXWhgGF6IzmNPUSPaDBaEBygxlV1GRAOy7lAFOiw2pEUHYloK/55oYBheiM5jwzFHl9H1E2IgZ5cRUb9ZbXa8va8MAHDvpcO40CMNGJ+Ric7BYLbiuzOOheluYJcR0YBsOl2HGo0RYf4K/j2RUzC8EJ3DlrwGGMw2JIT6YmJCsNjlEHm0tXtKAQB3XpIElQ83YKSBY3ghOoeuLqP56XFs4iYagKMVrThS0QaFTIq7LkkUuxzyEgwvRD/S3G7CzsIGAOwyIhqoN/eWAQDmpcciMlAlbjHkNRheiH7ki2M1sNgEjI8LwkiuRUHUb406E7455VhuYMmlyeIWQ16F4YXoBwRBwCeHKwEAN2fFi1wNkWf7+HAlLDYB6QnBGBcXJHY55EUYXoh+4HSNFvl1OijkUtyQzi4jov6y2QV8cKACAHBXNse6kHMxvBD9QFery1VjohDspxC5GiLPtbOwAdVtHQjy9cE8vhEgJxuU8LJmzRokJydDpVIhOzsbBw8evOD5n3zyCdLS0qBSqTB+/Hhs3LhxMMqkIc5kteHLzr2Mbs5KELkaIs/23n5Hq8tNmfGcHk1O5/Lw8tFHH2HFihV44okncOTIEaSnp+Pqq69GQ0PDOc/ft28fbr/9dixduhRHjx7FggULsGDBApw6dcrVpdIQt+VMA9oMFsQEqXBZarjY5RB5rKpWA7YXOJ7j72SXEbmAy8PL888/j2XLlmHJkiUYM2YMXnnlFfj5+WHt2rXnPP/FF1/E3Llz8bvf/Q6jR4/GU089hUmTJuGll15ydak0xH2S6+gyunFSHGRSru1C1F8fHqyAIACXpoZheESA2OWQF3JpeDGbzcjNzcWcOXO+/4VSKebMmYOcnJxz3icnJ6fH+QBw9dVXn/d8k8kErVbb40bUV3UaI3YVNgIAbspklxFRf5mtdnx0yPFG4K7sJJGrIW/l0vDS1NQEm82GqKioHsejoqJQV1d3zvvU1dX16fxVq1YhKCio+5aQwBce6ruPDlXCLgCTk0MwLNxf7HKIPNbmM/VoajcjMlCJOWOiLn4Hon7w+NlGK1euhEaj6b5VVlaKXRJ5GKvNjg8POgYX3sl3ikQD8nHnjL1bshLgw93YyUXkrvzh4eHhkMlkqK+v73G8vr4e0dHR57xPdHR0n85XKpVQKpXOKZiGpG35DajTGhHqr8A148/9OCOii6vTGLH7bFf3Kxd5JNdxaSxWKBTIzMzE1q1bu4/Z7XZs3boVU6dOPed9pk6d2uN8ANi8efN5zycaqPc6F9K6JSsBSjmndBL116dHqmAXgCnJoUhm9yu5kEtbXgBgxYoVuOeee5CVlYUpU6Zg9erV0Ov1WLJkCQBg0aJFiIuLw6pVqwAAv/71rzFz5kw899xzuO6667Bu3TocPnwYr776qqtLpSGovFmPXYWNkEg4pZNoIARBwPrcKgDATdxag1zM5eHl1ltvRWNjIx5//HHU1dUhIyMDmzZt6h6UW1FRAan0+wagadOm4YMPPsBjjz2GP/7xjxgxYgS++OILjBs3ztWl0hDUtXz5zJERSAj1E7kaIs91uLwVpU16+ClkuG58jNjlkJeTCIIgiF2EM2m1WgQFBUGj0UCtVotdDrkxo8WGqau2otVgweuLsjgzgmgAfr/+OD4+XIWbMuPx7M3pYpdDHqgvr98cCk5D1jenatFqsCAu2Bez0iLFLofIYxnMVnx9ohYAcDMH6tIgYHihIatr75XbpyRwRV2iAdh4sg56sw3JYX6YMixU7HJoCGB4oSHpeGUbcstb4SOT4JbJXNiQaCA+O9I5UDczHhIJ3wiQ6zG80JC0dm8pAGBeeiwiA1UiV0Pkueq1RuSUNAMA5mfEiVwNDRUMLzTk1GmM3f3z9146TORqiDzb/47XQBCAzKQQztijQcPwQkPOOzllsNoFZA8Lxbi4ILHLIfJoG47XAADmZ8SKXAkNJQwvNKR0mG34oHMfo3svY6sL0UCUNulxokoDmVSCa7m2Cw0ihhcaUj49UoU2gwWJoX6YM5rruhANxIZjjlaXS1PDER7APeZo8DC80JBhtwt4s3Og7uJpyZweTTQAgiDgy+PVAID56ewyosHF8EJDxs6zjShu1CNQKef0aKIBOl2jRUmjHkq5FFeNZSsmDS6GFxoyXt1ZAgC4ZXICApQu39aLyKttPOmYsXdFWiQCVT4iV0NDDcMLDQnHKtuQU9IMuVSCpRyoSzQggiBg06k6AMDccdEiV0NDEcMLDQmv7CgGANyQEYvYYF+RqyHybGcb2lHSpIdCJsUV3BeMRMDwQl6vuLEd355xvEv8+cwUkash8nxdrS6XjQhnlxGJguGFvN5ru0ogCMCc0ZEYGRUodjlEHu8bdhmRyBheyKvVa4347IhjOucvLmerC9FAlTfrkVerhUwqwZVcK4lEwvBCXm3tnlKYbXZMTg5BZlKo2OUQebxvTztaXS4ZHooQf4XI1dBQxfBCXkvTYcH7BxxbAXCsC5FzdHcZjWWXEYmH4YW81lt7y9BusmJUVCBmjeKMCKKBamo34VhlGwDgKoYXEhHDC3klrdGCN/Y4FqVbfkUqpNwKgGjAdhY0QhCAsbFqRKlVYpdDQxjDC3mlt/eWQWu0IjUygLvdEjnJtoIGAODaLiQ6hhfyOjqjBa/vcWzA+NAVqdyAkcgJrDY7dhU2AgAuZzcsiYzhhbzOOznl0HRYMDzCH9dP4G63RM6QW94KndGKED8fZCQEi10ODXEML+RV9CYrXt/tGOvCVhci59le4Gh1mTkygn9XJDqGF/Iq7+SUo9VgwbBwf8xjqwuR02zPd4x3mcXxLuQGGF7IaxjMVrzW2eqyfFYq5DI+vImcobqtAwX1OkgljpYXIrHx2Z28xjs55WjRm5EU5of5GWx1IXKWroG6ExNDEOzHVXVJfAwv5BW0Rgte2VkMAHjoihFsdSFyor1FTQCAy1LDRa6EyIHP8OQVXt9dijaDBamRAfjZxDixyyHyGna7gJziZgDApQwv5CYYXsjjNbeb8EbnWJeHrxzJmRBETlTYoEOz3gxfHxmnSJPbYHghj/fyjmLozTaMjwvC3HHcb4XImfYWOVpdJg8LhULOlwxyD3wkkker1XTgnf3lAIBHrh4FiYStLkTOtK9zvMulKWEiV0L0PYYX8mj/3loEs9WOKcNCMWME++OJnMlqs+NAaQsAjnch98LwQh6rrEmPjw9XAgB+x1YXIqc7Ua1Bu8mKIF8fjI5Ri10OUTeGF/JYL2wphM0uYNaoCExODhW7HCKvc6DE0epyyfBQDoQnt8LwQh7pVLUGXx6rAQA8fNUokash8k655Y7wwjcH5G4YXsjjCIKAv3+dBwCYnxGLcXFBIldE5H0EQUBueSsAIDMpRORqiHpieCGPs72gATklzVDIpXiErS5ELlHcqEerwQKlXIqxsXyDQO6F4YU8itVmx6qN+QCAJdOSkRDqJ3JFRN6pq8soPT6Y67uQ2+EjkjzKx4ercLahHcF+PvjlrFSxyyHyWofLOruMktllRO6H4YU8ht5kxfObCwEAv7piBIJ8fUSuiMh7dY13yeJ4F3JDDC/kMf67qwRN7SYkhfnhrkuSxC6HyGtpDBaUNOkBAJMSGV7I/TC8kEeo1xrx2i7H5ot/mJvGPngiFzpZrQEAJIb6IcRfIXI1RD/FVwDyCM9+W4AOiw2TEoNxDTdfJHKprvAynssQkJtieCG3d7yyDZ/kVgEAHrt+DLcBIHKxU53hhWsokbtieCG3ZrcLePJ/pwEAN06MY/870SBgywu5O4YXcmtfHKvG0Yo2+Clk+MM1aWKXQ+T1NAYLKloMAIBxcdyMkdwTwwu5Lb3Jiqe/cSxI9+CsVESpVSJXROT9TtV8P1g32I+Ddck9MbyQ21qzvQgNOhMSQ/2w9LJhYpdDNCScYpcReQCGF3JL5c16vL67FADw2HWjofKRiVwR0dBQUK8DAKRFB4pcCdH5MbyQW/rb13kw2+yYPiIcV46JErscoiHjbH07AGBEFMMLuS+GF3I7u882YvOZesikEjzOqdFEg8ZuF3C2wdHyMjIqQORqiM6P4YXcitlqx5MbHFOjF01N4rs/okFU2WqA0WKHQi5FUpi/2OUQnRfDC7mV1/eUoLhRjzB/BX4ze6TY5RANKYWdXUapEQGQSdniSe6L4YXcRmWLAf/eehYA8KfrRiPIj7tGEw2mwnp2GZFnYHght/GX/52G0WJH9rBQ/GxinNjlEA05RQ0crEuegeGF3MLmM/XYktcAuVSCvy0Yx0G6RCIoa9YDAIaHc7wLuTeGFxKdwWztHqS7bMZwvusjEklFs2NbgIRQP5ErIbowhhcS3b+3FqG6rQNxwb741RUjxC6HaEhqN1nRrDcDAJLCGF7IvTG8kKjO1uvw+u4SAMCTN4yFr4Ir6RKJobyzyyjUX4FAFQfLk3tjeCHRCIKAx744BatdwJzRUVxJl0hEXV1GiewyIg/A8EKi+ehQJQ6UtsDXR4YnbxgjdjlEQ1p5C8MLeQ6XhZeWlhbceeedUKvVCA4OxtKlS9He3n7B8x966CGMGjUKvr6+SExMxK9+9StoNBpXlUgiatAa8feNeQCAh68aifgQPmESiamiM7xwvAt5ApeFlzvvvBOnT5/G5s2b8dVXX2HXrl24//77z3t+TU0Nampq8Oyzz+LUqVN46623sGnTJixdutRVJZKIHv/yNHRGK9Ljg7Dk0mFil0M05NVpjACA2GBfkSshuji5K35oXl4eNm3ahEOHDiErKwsA8J///AfXXnstnn32WcTGxv7kPuPGjcOnn37a/XVKSgr+/ve/46677oLVaoVc7pJSSQSbTtVi0+k6yKUSPL1wApchJ3ID9VpHeIlWq0SuhOjiXNLykpOTg+Dg4O7gAgBz5syBVCrFgQMHev1zNBoN1Go1g4sX0XRY8PiXjjVdfj4zBaNj1CJXREQAUK81AQAi1UqRKyG6OJekgrq6OkRGRvb8RXI5QkNDUVdX16uf0dTUhKeeeuqCXU0AYDKZYDKZur/WarV9L5gGzaqNeWjQmTA8wh/Lr0gVuxwiAmC12dGsdzyPRrHlhTxAn1peHn30UUgkkgve8vPzB1yUVqvFddddhzFjxuDJJ5+84LmrVq1CUFBQ9y0hIWHAv59cY19xE9YdqgQAPH3jBKh8uKYLkTtoajdDEAC5VIJQP4XY5RBdVJ9aXh5++GEsXrz4gucMHz4c0dHRaGho6HHcarWipaUF0dHRF7y/TqfD3LlzERgYiM8//xw+PhdeLGnlypVYsWJF99darZYBxg0ZLTb88bOTAIA7sxMxZVioyBURUZeu8S4RgUpIOQaNPECfwktERAQiIiIuet7UqVPR1taG3NxcZGZmAgC2bdsGu92O7Ozs895Pq9Xi6quvhlKpxIYNG6BSXbz5UqlUQqlkH627e2FLIcqaDYhWq/CHa9LELoeIfqArvEQG8rmUPINLBuyOHj0ac+fOxbJly3Dw4EHs3bsXy5cvx2233dY906i6uhppaWk4ePAgAEdwueqqq6DX6/HGG29Aq9Wirq4OdXV1sNlsriiTBsnRila8tsuxBcBTC8ZBzaXHidxKq8Gxp1FYAMMLeQaXTeN5//33sXz5csyePRtSqRQLFy7Ev//97+7vWywWFBQUwGBwLIx05MiR7plIqak9B3KWlpYiOTnZVaWSCxktNjz8yXHYBWBBRiy3ACByQ5oOCwAgyJdvLMgzuCy8hIaG4oMPPjjv95OTkyEIQvfXl19+eY+vyTs8+20BShr1iAxU4skbxopdDhGdQ1d4Uau4LAV5Bu5tRC5zsLQFb+wtBQA8vXA8gjmLgcgtaTusANjyQp6D4YVcwmC24nfrj0MQgJsz43FFGruLiNxVd8sLwwt5CIYXcomnv8lHebMBsUEq/Hked4wmcmcc80KehuGFnG5fURPeySkHADxz0wTOLiJyczqjI7wE8m+VPATDCzmVzmjB79afAOBYjG76iIuvC0RE4jJZ7QAApQ9fEsgz8JFKTvWPjXmobutAfIgv/njtaLHLIaJeMHeFFzlfEsgz8JFKTvPd6Tp8eLASEgnwr5vS4a/ktEsiT2C2MbyQZ+EjlZyiQWvEHz51dBctmz4cU1PCRK6IiHrLZOkKL9wslTwDwwsNmN0u4JH1J9BqsGBMjBoPXzVS7JKIqA/Y8kKeho9UGrC3c8qwq7ARSrkUL96WwXdvRB7GZHHsH6dgeCEPwUcqDUhBnQ6rvskHAPzputEYERUockVE1F8SSMQugahXGF6o30xWG3697ijMVjtmjYrA3ZckiV0SEfWDVOIILXbuL0ceguGF+u1fmwqQX6dDmL8C/7wpHRIJ37UReaKuP10bwwt5CIYX6pc9Z5vw+h7Hpov/vGkCIgKVIldERP0lkzrSi8DwQh6C4YX6rLndhBUfHwPgWEV39mhuukjkyb7vNhK5EKJeYnihPrHbBTz8yXE06ExIifDHY9dx00UiT9fV5Wu1Mb2QZ2B4oT55bXcJdhQ4pkWvuXMSfBWcFk3k6fw6/447OqdME7k7hhfqtSMVrfjXtwUAgCfmjUVatFrkiojIGbrCi8FsFbkSot5heKFe0RgseOiDo7DaBVw/IQa3T0kQuyQicpKAzn3I9CaGF/IMDC90UYIg4PefHkd1WwcSQ/2w6sbxnBZN5EX8OsNLu4ndRuQZGF7oot7dX45vT9fDRybBS3dMRKDKR+ySiMiJApTsNiLPIhe7AHJvp6o1+NtXeQCAldeMxoT4YHELIrdlswtoN1nRbrKiw2yDyWqDyWqHyWL//nOrHSaLDRabAAECBAEQAEAQIDg+dK814iOXwkcmhUImhaLzcx+ZBAq5FEq5FH4KOQKUjpu/Us59eQbAX+F4KdAZGV7IMzC80HnpjBY89OFRmG12zBkdhSWXJotdEg0CQRCg6bCgWW9Gc7sZze0mNOnNaGk3o1lvQqvBAp3RAp3R+oOPjtAiJoVMCn+lDEG+Pgj1V3TfQvwVCPNXIMRPgbAABSICVIgJViHMX8Huz06h/goAQIveLHIlRL3D8ELnJAgCfr/+BEqb9IgNUuFfN03gE70XMFvtqNcaUasxolbTgVqNEXUaI2raOlCndXzeojfDOoDVyhQyKXwVMijlUih9pFDKHZ+rfDqPyaWQy6SQwLEsvQQSx8fOz7v2BrTa7DBb7bDYBJi7P3d8NFntMJgdgclosTv+bTY7zAY7Wg0WlDUbLl6nXIqYIBViglSIDfJFTLAKscG+SAr1x7AIf8SoVZBKh8ZjvmuF7EadSeRKiHqH4YXO6Y09pfjmVB18ZBKsuXMSQjrfmZF7s9sFNOhMKG/Wo7zFgIpmQ+dHParbjGhq7/2LU6BKjvAAJcI6WzDCApQID1Ag2E+BQJUcapUcgSofBPb4KIdSPrhr/1htduhNNrSbrdCbrGgzWNCiN6NFb0arwdF61Gowo1lvRovehAatCY3tJpitdpQ3G1B+nqCjlEsxLNy/xy0lMgBp0YHwU3jXU2d4gCO89OXxQSQm7/oLJKc4WNqCVd/kAwAev34MJiaGiFwR/ZjGYEFRow5FDe04W9+OsmY9ypsNqGgxwGS1X/C+P2xxiAny7fF5dJAK4QFKhPj7DHoI6S+5TIogPymC/Ho/kLyrBaqmzdH6VKPpQG2b4+uyZn33dcyv0yG/TtfjvhIJMCzMH6Nj1BgTq8bomECMjlEjWq3y2NZJtryQp2F4oR4adEY8+MER2OwCFmTE4q5LksQuaUhr1JlQWO8IKUUN7TjboENRg/6C75BlUgnign2RFOaHxFC/zo/+iA9xBJVQjvWAQi5FQqgfEkL9zvl9q82O6rYOlDTpUdqoR2mTHiVN7Sisb0ejzoSSJj1KmvT4+mRt931C/RWYmBCMSUkhmJgYjPT4YPgrPeMptju8sOWFPIRn/GXRoLDa7Fj+wVE06kwYGRWAf3A9l0FjtdlR0qRHXq0WZ2q1OFOjRV6t7oIhJSZIhdTIAKREBCAlwh+JYf5ICvVDXIgvfGSceTMQcpkUSWH+SArzx6xRPb/XqDMhr1bb/X+VV6tFcaMeLXoztuY3YGt+AwBHiBwVFYhJScGYnByKaSnhbrv7elSgCgDQZrDAYLZ6XbcYeR8+Qqnbv74twMHSFgQo5Xj5rkw+gbmIyWpDfq0Ox6vacKpag7xaHQrqdTCfo7tHIgGSQv2QGhmI1MiA7ltKhD/X2xFJRKASEYERmDEyovuY0WLDmVotjpS34mhFG45UtKJWY3QE0Vot3ttfAQBIiw7EtJRwXJoahuzhYd0r24otyM8HapUcWqMVlS0dGBUdKHZJRBfkHn85JLpNp2rx310lAIB/3TQBKREBIlfkHWx2ASWN7ThW2YYTVRqcqGrDmVotLOfYvddfIUNazPdjKMbEqDHKCweHeiOVjwyTEkMw6Qfjw2o1HThS3obc8lbklDQjr1bbPYZm7d5SyKUSZCQEY/boKFw5JhIpEQGitnQmhvnhVLUWFS0Ghhdye3xWJJQ0tuORT04AAJZNH4ZrxseIXJHnam434XB5K46Ut+J4VRtOVmmgN/90yfUQPx+kJwRjfFwQxnQO/EwI8RsyU3OHgpggX1w3wRfXTXD8PTW3m5BT0oy9Rc3YW9SEihYDDpe34nB5K57ZlI/kMD/MGR2FK8dEITMpBPJB7vpLCvXHqWotypv1g/p7ifqD4WWI0xktWPbOYbSbrJiSHIrfz00TuySPIQgCypoNOFTWgsNlLThc3oqSxp8+8fv6yDA+LgjpCUGYEB+MjIRgxIf4cjzREBMWoMT1E2Jx/YRYAEBliwE7Chux5Uw9coqbUdZswOt7SvH6nlKE+PngmvExmJ8ei8nJoYMSahPDHIOXS5sYXsj9MbwMYXa7gN9+dAzFjXpEq1V46c6JHOh5AVabHadrtJ1hpRWHy1vQ1P7TFUlHRgUgMykUExOCMSEhCKkRAYP+LprcX0KoH+6+JAl3X5KEdpMVuwsbsTmvHtvyG9BqsOCDAxX44EAFYoJUmJceixvSYzE2Vu2y0DsyytFVXFivu8iZROJjeBnCVm8pxJa8BijkUry6KBORnTMOyMFuF5Bfp8O+4ibkFDfjQGnLT5bAV8ikSE8IQlZyKLKSQpCZFIJgPy7oR30ToJTjmvExuGZ8DKw2O/aXtODLY9XYdKoOtRojXt1Vgld3lWBkVABun5KIGyfG92ldm95Ii1YDAPLrdBAEgS2D5NYkQtcuaF5Cq9UiKCgIGo0GarVa7HLc1qZTtfj5e0cAAM/fko4bJ8WLXJH4BEFAcaMeOcVN2FfcjP0lzWg1WHqco1bJMTk5FFnJoZicHIJxcUFQ+XjGYm7keYwWG3YUNGLD8WpsyWvonpGmlEtx3fgY3J6diKykEKcEDbPVjjGPb4LVLmDvo1cgLth3wD+TqC/68vrNlpchKL9OixUfHwcALL1s2JAOLq16M3YXNWFnQSP2FDWiXttzXRV/hQxThoViakoYpqWEY3SMGjIOqqVBovKRYe64aMwdFw1NhwUbjlXj/QMVyK/T4bOj1fjsaDVGRgVg6WXDMD8jbkBBWiGXIiUiAAX1OuTXahleyK2x5WWIaTOYccNLe1HRYsBlqeF4a8nkITUew2qz43hVG3YWNGLn2SacqGrDD/8CFHIpspJCMC0lDFNTwjEhPojjgMitCIKA41UafHigAhuO16DD4pjNFh6gxD1Tk3DXJUn93ovstx8dw+dHq/Hr2SPw2ytHOrNsootiywudU9cKuhUtBiSE+uI/t08cEsGlQWfE9vwG7CxsxJ6zTdAae45bSYsOxMyRjkXHMpNC2A1Ebk0icawPk5EQjD9dPxrrDlbgzb1lqNUY8dzmQqzZUYTbJifil5enIFLdt3Fsk5JC8PnRauSWt7qoeiLnYHgZQv6xMR97iprg6yPDq3dnee1O0YLgGGi75Uw9tuQ34HhlW4/vB/n6YPqIcMwYGYEZIyIQHcSByuSZ1Cof3D8jBUsuHYaNJ2vx6q4SnK7R4q19ZfjwYAXuzE7Czy8f3uvB+FlJjkX2jlS0wmqzD4k3N+SZGF6GiPcPlGPt3lIAwHO3pGN0jHd1qZmsNhwoacHWvHpsyWtAdVtHj++nxwfh8lGRuHxUBCbEB3PcCnkVH5kU8zPicEN6LPYUNWH1lrPILW/F2r2l+OBgOe6ZmoxfzkpFkO+FZyiNjApEoFIOncmK/DodxsUFDdK/gKhvGF6GgD1nm/D4l6cBAA9fORLXeskKulqjBVvz6rH5TD12FjT2WMlWKZdi+ohwzB4dhdlpkX1uPifyRBKJBNNHROCy1HDsOtuEFzYX4lhlG/67qwQfH67Eb68ciTumJJ63RUUmlWBiUgh2FTbiQGkLwwu5LYYXL1fU0I5fvJ8Lm13AgoxYLL8iVeySBqTNYMZ3Z+rxzcla7Clq6rFHUESgEnNGR2J2WhQuTQ2Hr4JjV2hokkgkjnFcI8KxvaAB/9iYj6KGdjz+5Wm8k1OOx68f02NjyR+6LDUMuwobsauwEUsvGzbIlRP1DmcbebFWvRkL/m8vypsNyEwKwfv3ZXvkYNSmdhO+O12Pb07VIqe4GVb79w/Z1MgAzB0bjSvHRGF8XBD3BiI6B6vNjg8PVuD5zYXdaxfNz4jFn68fg/AAZY9zC+t1uOqFXVDKpTj2+FV8E0CDhrONCGarHQ+8l4vyZgPiQ3zx37szPSq4NOiM2HSqDhtP1uJgaQt+kFeQFh2Ia8fH4Jpx0RgRxd1viS5GLpPi7qnJuCEjDi9uOYu39pXiy2M12FnYiD9dOxo3ZcZ3L3Q3IjIAsUEq1GiM2F/SjFlpkSJXT/RTDC9eSBAE/OnzkzhY2oIApRxrF0/+ybsrd6TpsODb03XYcKwG+4qbegSWCfFBmDsuGteMi8GwcH/xiiTyYEG+Pnh83hgsmBiLRz89iTO1Wvxu/Ql8fbIW/7opHRGBSkgkEsxKi8T7Byqw6VQdwwu5JXYbeaFXdhbj6W/yIZUAaxdPxuWj3PfJx2ixYXt+A744Vo3t+Y0w2+zd38tICMZ142Mwd1w0EkL9RKySyPtYbHas3VOK5zcXwmS1IzxAgWdvTsfloyKRU9yM21/bD7VKjsOPXQmFnFOmyfXYbTSEfXWiBk9/kw8AeGLeWLcMLlabHTklzfjyWA2+PVUH3Q82OxwRGYAFE+Mwb0IsEsMYWIhcxUcmxQMzUzArLRK/+vAo8ut0WPzmIdx32TD8fm4aIgOVaNCZsPtsI2aPjhK7XKIe2PLiRQ6UNOPuNw7CbLNj8bRkPHnDWLFL6qGgTof1uZX4/GgNmtq/30MoLtgX89JjMT8jFmnRgdzNlmiQGS02rNqYh7dzygEAl6WGI1KtxGdHqjF3bDReuTtT5AppKGDLyxBU1KDDsncOw2yz4+qxUfjz9WPELgmAY8bThuM1WJ9bhZPVmu7jIX4+uG5CDOZnxCEzMYSzhIhEpPKR4S/zx2FqShh++9Fx7Clq6v7e5rx61Go6EBPEjRrJfTC8eIEGrRH3rD0ErdGKSYnBePG2iaKuIGu12bGzsBHrc6uwJa++ey0WH5kEs9OicFNmPGaOiuCGh0RuZu64GCSG+mPZO4e7V6m22QV8eKACK64aJXJ1RN9jt5GH05usuPXVHJyq1mJYuD8+/cU0hIq0Z1FRQzs+PlyJz45U9+gWGhurxk2Z8ZifESdabUTUe83tJtz71iEcr/q+tfTkk1chUHXh7QWIBoLdRkOExWbHL98/glPVWoT5K/DWksmDHg6MFhu+PV2HDw5U4EBpS/fxMH8FFkyMw8JJ8RgT6/0hksibhAUo8d592bj3rUM4VObYYXrlZyfx0h2TRK6MyIHhxUMJgoA/f3EKOwsbofKR4o3Fk5EUNnjrnxQ1tGPdwQp8eqSqe8VOqQS4Ii0St2QlYFZaJLuFiDxYoMoH7y7NRtqfNwEAvjpRi0evMSA+hLMASXwMLx7q+c2FWHeoElIJ8NLtk5CREOzy33m+VpaYIBVunZyAWycncFAfkRdR+chw9M9XYuJTmwEAFS0ML+QeGF480No9pfjPtiIAwFMLxmHOGNeuwVDRbMC7+8uwPvenrSx3ZCdi5shIUQcIE5HrhPgrcPLJq3Cssg1Th4eJXQ4RAIYXj/P50Sr89aszAIBHrhqJO7OTXPJ7BEHAnqImvL2vDFvzG9A1rJutLERDT6DKB9NHnHsXaiIxMLx4kG359XjkkxMAgHsvHYYHZ6U6/XfoTVZ8dqQKb+eUo6ihvfv4jJERWHRJEmalsZWFiIjExfDiIQ6VteAX7x2BzS7gxolxeOy60U5dibasSY+3c8qw/nBV93L9/goZbsqMx6JpyUiJCHDa7yIiIhoIhhcPkFerxb1vHYLJascVaZF45qYJTlmRtqtraO2eUuwobOzuGhoe7o9FU5OwMDOe6zoQEZHbYXhxcxXNBixaexA6oxWTk0Ow5o5JA56CbLba8dWJGry6qwT5dbru47NGReCeacmYMSKCy/UTEZHbYnhxY/VaI+564wAadSakRQfi9Xsmw1ch6/fP0xotWHewAmv3lKFOawQA+ClkuCUrAfdMS8aw8MFbJ4aIiKi/GF7cVFO7CXe8th8VLQYkhvrhnaVTEOTbvy6cmrYOvLm3FB8erER753iWiEAlFk9Lxl3ZSQjyY9cQERF5DoYXN9SqN+Ou1w+guFGP2CAV3r8vG5GBqj7/nPw6Lf67swT/O14Dq90xoGVEZACWzRiO+RmxUMr734pDREQkFpet397S0oI777wTarUawcHBWLp0Kdrb2y9+RzgGkl5zzTWQSCT44osvXFWiW9J0WHD32gPIr9MhMlCJD5ZdgoTQvq1oeayyDfe9fRhzV+/G50erYbULmDo8DG8unoxvfzMDt2QlMLgQEZHHclnLy5133ona2lps3rwZFosFS5Yswf33348PPvjgovddvXq1U6cBe4p2kxWL3zzYvdHiB8uykdyHcSgHSprx0vYi7D7bBACQSIBrx8XggZnDMSE+2EVVExERDS6XhJe8vDxs2rQJhw4dQlZWFgDgP//5D6699lo8++yziI2NPe99jx07hueeew6HDx9GTEyMK8pzSx1mG+596xCOVrQh2M8H792XjdTIwIveTxAE7CxsxJrtRd27v8qkEizIiMMvLk9BaiTXZyEiIu/ikvCSk5OD4ODg7uACAHPmzIFUKsWBAwfws5/97Jz3MxgMuOOOO7BmzRpER0f36neZTCaYTKbur7Va7cCKF4HRYsOydw7jYGkLApVyvHtvNkbHqC94H7tdwHdn6rFmexFOVmsAAAq5FLdkxeOBGSl97moiIiLyFC4JL3V1dYiMjOz5i+RyhIaGoq6u7rz3++1vf4tp06Zh/vz5vf5dq1atwl/+8pd+1yo2o8WGX7yXiz1FTfBTyPDWvVMwPj7ovOcLgoAteQ14YXMhztQ6gpqvjwx3XZKIZdOHI1Ld94G9REREnqRP4eXRRx/FM888c8Fz8vLy+lXIhg0bsG3bNhw9erRP91u5ciVWrFjR/bVWq0VCQkK/ahhsRosND7ybi52FjVD5SLF28WRkJoWc89yu7qEXNhfieJWjpSVAKcfiacm497JhCPVXDGbpREREoulTeHn44YexePHiC54zfPhwREdHo6Ghocdxq9WKlpaW83YHbdu2DcXFxQgODu5xfOHChZg+fTp27NhxzvsplUoolcre/hPcRofZ0VW0p6gJvj4yvLE4C5ecZ7v5fUVNeG5zIXLLHWNafH1kWHxpMu6fPhwhDC1ERDTE9Cm8REREICLi4tuiT506FW1tbcjNzUVmZiYARzix2+3Izs4+530effRR3HfffT2OjR8/Hi+88ALmzZvXlzLdnsFsxdK3DiOnpBl+ChneXDwZ2ecILgdLW/D85gLsL2kBACjlUtx9SRJ+fnkKwgM8L7ARERE5g0vGvIwePRpz587FsmXL8Morr8BisWD58uW47bbbumcaVVdXY/bs2XjnnXcwZcoUREdHn7NVJjExEcOGDXNFmaLQm6xY8tYhHCxtQYBSjreWTEZWcmiPc07XaPDMpgLsKmwEAChkUtw+JQG/nJWKKI5pISKiIc5l67y8//77WL58OWbPng2pVIqFCxfi3//+d/f3LRYLCgoKYDAYXFWC22k3WbF47UEcLm9FoFKOt5dOwaTE78e4VLYY8PzmQnx+tBoAIJdKcMvkBCyflYrYYF+xyiYiInIrEkEQBLGLcCatVougoCBoNBqo1ReebjyYqts6cOnT2wAAapUc7y7NRnpCMADHdgBrthfhnZxymG12AMC89Fg8ctVIJIVxs0QiIvJ+fXn95t5Gg6Cp3dQdXADg7XunID0hGEaLDW/uLcP/7SiCzujYMHFaShgevSaNK+ISERGdB8OLi1W3deDu1w90f/3szelIjw/Gp7lVePa7AtRqjACAtOhAPHpNGmaOjBiSWyMQERH1FsOLCxU3tuPu1w+gRmNEXLAv3rsvG60GM3728j4cr2wDAMQGqfDwVaOwYGIcZFKGFiIiootheHGRU9UaLFp7EC16M1Ii/PHszen499az3YNx/RUyLL9iBJZcmgyVD3d4JiIi6i2GFxc4WNqCpW8dgs5kxYjIAFyaGo47XjuADosNEglwc2Y8Hrl6FCIDOe2ZiIiorxhenOy703V46MOjMFkds4YadCa8ta8MAJCVFIIn5o294N5FREREdGEML0709r4yPPm/0/jh5HNNhwUxQSqsvHY05k2I4WBcIiKiAWJ4cQK7XcDTm/Lx6q6SHsd9ZBI8MCMFD85Kha+C41qIiIicgeFlgIwWGx755Di+OlHb4/jU4WF4asE4pEYGiFQZERGRd2J4GYA2gxm3/nc/Cup13cfCAxR47LoxmJ8Ryy4iIiIiF2B46afixnbMfm5nj2N3XZKI312VhiA/H5GqIiIi8n4ML/3waW4VHv7keM9jv5iGzKSQ89yDiIiInIXhpQ8EQcDP/m8fjnWujgsAt2TF46/zx3GhOSIiokHC8NJLdruA4X/c2OMYW1uIiIgGH8NLL/1wUG6gSo7cx66EQi4VsSIiIqKhieGll0ZGBeKfN01AmL8Cs0dHiV0OERHRkMXw0ksyqQS3ZCWIXQYREdGQx34PIiIi8igML0RERORRGF6IiIjIozC8EBERkUdheCEiIiKPwvBCREREHoXhhYiIiDwKwwsRERF5FIYXIiIi8igML0RERORRGF6IiIjIozC8EBERkUdheCEiIiKP4nW7SguCAADQarUiV0JERES91fW63fU6fiFeF150Oh0AICEhQeRKiIiIqK90Oh2CgoIueI5E6E3E8SB2ux01NTUIDAyERCIRrQ6tVouEhARUVlZCrVaLVoc74rW5MF6f8+O1uTBen/PjtTk/d7k2giBAp9MhNjYWUumFR7V4XcuLVCpFfHy82GV0U6vV/EM5D16bC+P1OT9emwvj9Tk/Xpvzc4drc7EWly4csEtEREQeheGFiIiIPArDi4solUo88cQTUCqVYpfidnhtLozX5/x4bS6M1+f8eG3OzxOvjdcN2CUiIiLvxpYXIiIi8igML0RERORRGF6IiIjIozC8EBERkUdheHGilpYW3HnnnVCr1QgODsbSpUvR3t5+0fvl5OTgiiuugL+/P9RqNWbMmIGOjo5BqHjw9PfaAI5VF6+55hpIJBJ88cUXri1UJH29Pi0tLXjooYcwatQo+Pr6IjExEb/61a+g0WgGsWrXWLNmDZKTk6FSqZCdnY2DBw9e8PxPPvkEaWlpUKlUGD9+PDZu3DhIlYqjL9fntddew/Tp0xESEoKQkBDMmTPnotfTk/X1sdNl3bp1kEgkWLBggWsLFFFfr01bWxsefPBBxMTEQKlUYuTIke71tyWQ08ydO1dIT08X9u/fL+zevVtITU0Vbr/99gveZ9++fYJarRZWrVolnDp1SsjPzxc++ugjwWg0DlLVg6M/16bL888/L1xzzTUCAOHzzz93baEi6ev1OXnypHDjjTcKGzZsEIqKioStW7cKI0aMEBYuXDiIVTvfunXrBIVCIaxdu1Y4ffq0sGzZMiE4OFior68/5/l79+4VZDKZ8M9//lM4c+aM8Nhjjwk+Pj7CyZMnB7nywdHX63PHHXcIa9asEY4ePSrk5eUJixcvFoKCgoSqqqpBrtz1+nptupSWlgpxcXHC9OnThfnz5w9OsYOsr9fGZDIJWVlZwrXXXivs2bNHKC0tFXbs2CEcO3ZskCs/P4YXJzlz5owAQDh06FD3sW+++UaQSCRCdXX1ee+XnZ0tPPbYY4NRomj6e20EQRCOHj0qxMXFCbW1tV4bXgZyfX7o448/FhQKhWCxWFxR5qCYMmWK8OCDD3Z/bbPZhNjYWGHVqlXnPP+WW24Rrrvuuh7HsrOzhQceeMCldYqlr9fnx6xWqxAYGCi8/fbbripRNP25NlarVZg2bZrw+uuvC/fcc4/Xhpe+XpuXX35ZGD58uGA2mwerxD5jt5GT5OTkIDg4GFlZWd3H5syZA6lUigMHDpzzPg0NDThw4AAiIyMxbdo0REVFYebMmdizZ89glT0o+nNtAMBgMOCOO+7AmjVrEB0dPRiliqK/1+fHNBoN1Go15HLP3LLMbDYjNzcXc+bM6T4mlUoxZ84c5OTknPM+OTk5Pc4HgKuvvvq853uy/lyfHzMYDLBYLAgNDXVVmaLo77X561//isjISCxdunQwyhRFf67Nhg0bMHXqVDz44IOIiorCuHHj8I9//AM2m22wyr4ohhcnqaurQ2RkZI9jcrkcoaGhqKurO+d9SkpKAABPPvkkli1bhk2bNmHSpEmYPXs2zp496/KaB0t/rg0A/Pa3v8W0adMwf/58V5coqv5enx9qamrCU089hfvvv98VJQ6KpqYm2Gw2REVF9TgeFRV13utQV1fXp/M9WX+uz4/94Q9/QGxs7E8Cn6frz7XZs2cP3njjDbz22muDUaJo+nNtSkpKsH79ethsNmzcuBF//vOf8dxzz+Fvf/vbYJTcKwwvF/Hoo49CIpFc8Jafn9+vn2232wEADzzwAJYsWYKJEyfihRdewKhRo7B27Vpn/jNcwpXXZsOGDdi2bRtWr17t3KIHkSuvzw9ptVpcd911GDNmDJ588smBF05e6emnn8a6devw+eefQ6VSiV2OqHQ6He6++2689tprCA8PF7sct2O32xEZGYlXX30VmZmZuPXWW/GnP/0Jr7zyitildfPM9uVB9PDDD2Px4sUXPGf48OGIjo5GQ0NDj+NWqxUtLS3n7fKIiYkBAIwZM6bH8dGjR6OioqL/RQ8SV16bbdu2obi4GMHBwT2OL1y4ENOnT8eOHTsGUPngcOX16aLT6TB37lwEBgbi888/h4+Pz0DLFk14eDhkMhnq6+t7HK+vrz/vdYiOju7T+Z6sP9eny7PPPounn34aW7ZswYQJE1xZpij6em2Ki4tRVlaGefPmdR/rejMpl8tRUFCAlJQU1xY9SPrzuImJiYGPjw9kMln3sdGjR6Ourg5msxkKhcKlNfeK2INuvEXXoMvDhw93H/v2228vOOjSbrcLsbGxPxmwm5GRIaxcudKl9Q6m/lyb2tpa4eTJkz1uAIQXX3xRKCkpGazSB0V/ro8gCIJGoxEuueQSYebMmYJerx+MUl1uypQpwvLly7u/ttlsQlxc3AUH7F5//fU9jk2dOtWrB+z25foIgiA888wzglqtFnJycgajRNH05dp0dHT85Pll/vz5whVXXCGcPHlSMJlMg1m6y/X1cbNy5UohKSlJsNls3cdWr14txMTEuLzW3mJ4caK5c+cKEydOFA4cOCDs2bNHGDFiRI/prlVVVcKoUaOEAwcOdB974YUXBLVaLXzyySfC2bNnhccee0xQqVRCUVGRGP8El+nPtfkxeOlsI0Ho+/XRaDRCdna2MH78eKGoqEiora3tvlmtVrH+GQO2bt06QalUCm+99ZZw5swZ4f777xeCg4OFuro6QRAE4e677xYeffTR7vP37t0ryOVy4dlnnxXy8vKEJ554wuunSvfl+jz99NOCQqEQ1q9f3+MxotPpxPonuExfr82PefNso75em4qKCiEwMFBYvny5UFBQIHz11VdCZGSk8Le//U2sf8JPMLw4UXNzs3D77bcLAQEBglqtFpYsWdLjSaK0tFQAIGzfvr3H/VatWiXEx8cLfn5+wtSpU4Xdu3cPcuWu199r80PeHF76en22b98uADjnrbS0VJx/hJP85z//ERITEwWFQiFMmTJF2L9/f/f3Zs6cKdxzzz09zv/444+FkSNHCgqFQhg7dqzw9ddfD3LFg6sv1ycpKemcj5Ennnhi8AsfBH197PyQN4cXQej7tdm3b5+QnZ0tKJVKYfjw4cLf//53t3pjJBEEQRjcjioiIiKi/uNsIyIiIvIoDC9ERETkURheiIiIyKMwvBAREZFHYXghIiIij8LwQkRERB6F4YWIiIg8CsMLEREReRSGFyIiIvIoDC9ERETkURheiIiIyKMwvBAREZFH+X+NH3w+3BGTDQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get the path from the matplotlib contour plot of the critical line\n", "paths = CS.collections[0].get_paths()\n", "caustic_paths = []\n", "for path in paths:\n", " # Collect the path into a descrete set of points\n", " vertices = path.interpolated(5).vertices\n", " x1 = torch.tensor(list(float(vs[0]) for vs in vertices))\n", " x2 = torch.tensor(list(float(vs[1]) for vs in vertices))\n", " # raytrace the points to the source plane\n", " y1,y2 = sie.raytrace(x1, x2, z_s, packparams)\n", "\n", " # Plot the caustic\n", " plt.plot(y1,y2)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "1e77da2e", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "PY39", "language": "python", "name": "py39" }, "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.9.5" } }, "nbformat": 4, "nbformat_minor": 5 }