Skip to content

Instantly share code, notes, and snippets.

@jobovy
Last active November 21, 2022 17:31
Show Gist options
  • Save jobovy/326b8ef3b6ef6b0ea3c9ac9d8e1ff980 to your computer and use it in GitHub Desktop.
Save jobovy/326b8ef3b6ef6b0ea3c9ac9d8e1ff980 to your computer and use it in GitHub Desktop.
Frequency vs. action for different vertical potentials
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/Users/bovy/miniforge3/envs/py39/lib/python3.9/site-packages/jax/_src/lib/__init__.py:34: UserWarning: JAX on Mac ARM machines is experimental and minimally tested. Please see https://github.com/google/jax/issues/5501 in the event of problems.\n",
" warnings.warn(\"JAX on Mac ARM machines is experimental and minimally tested. \"\n",
"\n",
"%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"import numpy\n",
"import astropy.units as u\n",
"from astropy.constants import G\n",
"import galpy.util.conversion\n",
"galpy.util.conversion._APY_UNITS=True\n",
"from galpy.potential import (MWPotential2014, toVerticalPotential,\n",
" IsothermalDiskPotential,turn_physical_on,\n",
" evaluatelinearPotentials)\n",
"from galpy.potential.mwpotentials import McMillan17\n",
"from galpy.orbit import Orbit\n",
"from galpy.actionAngle import actionAngleVertical\n",
"from galpy.util.conversion import get_physical, freq_in_Gyr\n",
"%pylab inline\n",
"scott_ro, scott_vo= 0.3, 30.\n",
"from galpy.util import plot as galpy_plot\n",
"galpy_plot.start_print(axes_labelsize=17.,text_fontsize=12.,\n",
" xtick_labelsize=15.,ytick_labelsize=15.)\n",
"from matplotlib.ticker import FuncFormatter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Computing frequency vs. action for different vertical potentials"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Potential models"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we define three potentials:\n",
"* Scott's isothermal potential\n",
"* `MWPotential2014` as a function of $z$ at the solar position (8 kpc)\n",
"* The McMillan (2017) as a function of $z$ as its solar position"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Scott's isothermal\n",
"iso= IsothermalDiskPotential(amp=1.,sigma=30./220.)\n",
"iso._H= (300./8000.)\n",
"turn_physical_on(iso,ro=8.,vo=220.)\n",
"# MWPotential2014\n",
"mwp14_v= toVerticalPotential(MWPotential2014,1.)\n",
"turn_physical_on(mwp14_v,ro=8.,vo=220.)\n",
"# McMillan (2017)\n",
"mcm_v= toVerticalPotential(McMillan17,1.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $H(J)$ vs. $J$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we compute the Hamiltonian as a function of action for the three potentials:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1062d3f40>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAFCCAYAAACKKGi2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABQtElEQVR4nO3deVxU1f/H8dcVBXdZ3HfBfRdGLbM0xbIyMxX12/ItMzHbS3MrSzMzl9JyBbdKMxdM+7ZpYFn+XBIYFdyVURQVRXFwZZ3z+2OAUFG2WeHzfDx8JMPcued6zTfn3HPOR1NKIYQQQhR3pezdACGEEMIWJPCEEEKUCBJ4QgghSgQJPCGEECWCBJ4QQogSQQJPCCFEiVDa3g0orKpVq6qGDRsW+XMSEhKoVq1a0RvkoOT6nJtcn3OT67O9yMjIi0qpXBvltIHXsGFDIiIiivw5Op3OIp/jqOT6nJtcn3OT67M9TdNi7/Y9GdIUQghRIkjgCSGEKBFKfOAFBgbauwlWJdfn3OT6nJtcn2PRnHUvTZ1Opxxt7FgIIYR9aZoWqZTS5fa9Et/DE0IIUTLYJfA0TfPWNG2dpmn+t73uq2namMxfd3xfCCGEKCybL0vIEWLeuXzbXyk1I/N97sAJTdN6KqX0tmqfEEKI4snmPTylVJhSKgxIzPm6pmm+wPgc7zMCEYD08oQQQhSZwzzDy+zFBdz2sjdgtH1rrEOv1xMQEICfnx9hYWEFPt5oNN7ydVhYGB4eHhZqXcGFhYXh4+Njt/MLIURBOEzggbn3l/V7TdO8AU9gbW7vTUhIQKfTZf8KDg62VTMLzdfXl169eqHT6fD3L3jHde3aW/8o/P390elynYxkE/7+/nh75zYyLYRwRoaEaxy/cM3ezSiQ4ODgW7IAqHq39zry1mJBQM/Moc07VKtWzeG2tLG2oKAgp1v3IoRwfJevp/LllmOs3BVLt6bVWPpiR3s3Kd8CAwNv+XdR07SLd3uvQ/XwsmiaNgaYXtwnqxgMBsLCwtDr9YwdOzb79ZCQEMLCwrL/C+bhQ6PRSEhICHr9rX8ser0evV7PiBEjbnk9ODiYsLAwZsyYkX2sj48PYWFhBAQEsHHjRjw8PNDr9YSFhTFixAgMBkP2ZxkMhuzPCgsLIywsjLFjx94xtCqEcE6p6SaWbDPQbeaffLvzJIM61uOzAW3t3SyrcbgenqZpA4GwrLDTNM1bKWXI47C7mvzTAQ6evWKx9uWmZe3KfPRkqwIfFxQUxPjx43F3d88OF71ej8FgYMyYMQCMGDEiewjU3d2dgQMH5vpZvr6+uLu7ExYWhr+/P3q9npiYGAIDA9HpdIwdO5agoCCCgoLw9PRk+vTpeHt7M3/+fMA8PBkaGkpISAhjxowhICCAoKAgpk+fDsC6desICgoCYNq0admvCyGcj1KKzQfOM+23Q8ReukG3ptWY8HgLmtWsZO+mWZVD9fAylywYc4SdO+Br10ZZ0eDBg/Hz82PEiBHZz8KCgoLw9f33krN6ZPeSdayXl1f2a2vWrMHLyys7QLMCNTExEV9f31ueveU8Puv3np6et5wjK+zgzskzQgjnER2XxODgXbyyMhJXl1J8PbQj37zUqdiHHdhnHZ4v5qUGOmBsZg8uOHOSSmjme3Ie4leU8xWm52Ur3t7eREZGZg8xxsTE5Pq+xMRbVnBgMBhuCSx3d/c7jjEajfTq1Ss7PENDQ4E7g+z243P7LICxY8fSq1evXI8XQji++KRkZm4+wg974vAs78rUp1szWFeP0i4O1e+xKnusw9MrpWYopTyUUr2UUsGZrxuUUlouv4rtc7zg4ODsYcqsZ2MBAQG3PKMLDw/PntGZFTY5n63dTUBAQHbIAXc89ytoO728vG6ZWVqUzxNC2M7N1AzmhB3l4Vlb+WnfWUY85MOf73Xn2c4NSlTYgYMNaRZ3er2e0NBQIiIisocpsyaheHt74+7unv2sLmvSSs7hzhEjRtyy/CJruHLGjBkYDAZCQ0NZt24dRqMRf3//7OHQnBNfIiIisj/jXscHBQURFhaGwWDA39+fmJiYW4ZWExMTs493hiUhQpQ0JpPiB30cD8/aypywY/RoXp0to7ox7rHmVC5bxt7NswupliCEEMVMxMlEpvx8kH1xSbStW4WJfVrSsWHJeBxxr2oJDjdLUwghROGcTrzBZ5sO80vUOWpWLssXg9rRr30dSpXS8j64BJDAE0IIJ3ctJZ2FW4+zeNsJSmnwVs8mjOjmTXlX+Sc+J/nTEEIIJ2UyKUL0cczcfISEqyk83aEOY3o3o1aVcvZumkOSwBNCCCe0+0QiH/98gP1nrtChvjvBz/vRob79NpN3BhJ4QgjhRE4n3uCz3w7zS/Q5alUpy5dD2tO3Xe3b1y+LXEjgCSGEE7ieks7CrTEEbzNQSoO3/Zsw4iEfyrm62LtpTkMCTwghHJjJpNiw5wwzNh/m/JUU+rWvzdjHmstzukKQhec2ZDAYGDFiBH5+d+6WZjQa8fDwYMaMGYB5QbqHh0d2FQWDwZBdPDbnRtN+fn6EhITc8v2QkBCCg4MZO3ZsnjuiFHVfzMIcL4Vjhcgf/anLPL1wB6PW7aNm5bKsH9mFOUM6SNgVkvTwbMjb25uAgAASExPv2A/TYDDg6emZXSVh4MCBhIaGZgeDt7c348ePZ/jw4dnH+fr6Mn369Owtv3r16kVkZOQtFRU8PDw4ceLEXffIXLt2bZFq7OV2vI+Pz133BYU7C8eGhISQmJhIZGQkAQEB2dcTEhKCu7s7RqMRb2/v7H1BjUZj9rZsubU967jCFNkVwhHEJyUzfdNhNuw5Q/VKbnwe0I6nO8h6uqKSHp4dDB48+JbqA3cTEBDAunXrsr/OCsqcvaqcmznntrGzp6fnPffezE877iW34yMjI/N9fNa2aoGBgQQFBREQEACYQy00NBR/f//svUaz3GuHnaxt0YRwRslpGcz74xg9Pt/KL9HneO1hH/4c3Z0BfnUl7CxAAs8OBg4ceMu+lLf39rL4+/vf8Y/7oEGDWLt2LWAOi5ylhHKTmJiIt7d3gYrK3l44Niws7I5CsXc7PmuYNUtehWMTExNvCShPT0/0ej1r1669ZdjT3d09+xxZ+43mJiIigl69et3zz0QIR6OUYtP+c/h/8Rezfj/KQ02qseXdbrz3aHMquMlAnKUU/z/J38ZBfLR1z1GzDTz2WYEO8fb2zg6srCG73Oh0uuyirp6engQEBDB9+nQCAwPvKBsE/1ZRNxqNhIeHs2XLlux6ePkpKnu3wrE6nXlruqxCsVltuv3422vt5VU41t/f/5ahx6x6fWFhYbeEmqenZ67Xm5Ner0en00klB+FUDsdfYfL/DrLTcInmNSux6uXOdGlc1d7NKpaKf+A5qBEjRrBmzZo8e2hZw5qenp7Z780a9suNt7d3doBkBdGIESNuOSarikJu1dNzFo6FW0sR5VZoNi8FKRw7duzYW4Zwb5efCTJ36/kJ4WguX0/li9CjfPdPLJXLlWFKv9b8p2PJqk9na8U/8ArY87IVf39/RowYwYgRI+4ZeoMGDWLs2LG3BJZOp2PGjBlFmmxyt6KydyscC/cOk7sNy+a3cGxISAiDBw/OPm/WZJWc7b1bLxjMw7BZzyvDw8O5dOkS3t7e9zxGCHtIzzDx/e5TfB56lCs303j+vga806sp7uVd7d20Yk9+lLAjf39/QkJC7vked3f37CDKEhAQwJo1a/LdmylIUdmCFo69V1Ha/BaOzRra9fX1zR5+HTRo0C0zPY1G4z1/MAgMDGTgwIEMHDgQb29vevXqJWEnHM7OmEv0mft/TPzxAC1qVubXtx5k8lOtJexspPj38ByIwWBg+vTp6PV6xowZc0uvLWtqfnBw8B09t8GDB98SGoMGDbpjeM9gMLBmzZrsZ3g53+/v73/Ls73cispmfZ3zvfDvJJKsQrFZyyViYmLQ6XR3HJ+zKKy/vz/Tp08nLCwsOxhvLxyr0+no2bPnLd+/fPkyYF5mkdWOnLM0w8LCCA0NzX72eXug6vX67O9J6AlHcMZ4k09/PcQvUeeo416OBc/68ljrmrIdmI1JAVghhLCS5LQMgv82sGDrcZSCkd19ZDswK5MCsEIIYUNKKX4/eJ4pPx8k7vJNHm9TkwmPt6CuR3l7N61Ek8ATQggLOn7hGpN/OsC2YxdpVqMSq4Z3pouPLDNwBBJ4QghhAVeT0/hqyzGWbz9JeVcXJj3ZkufuayDLDByIBJ4QQhRBVjWDzzYd5uK1FAbr6vHeo83wquhm76aJ20jgCSFEIe0/k8SHP+5Hf8pI+3ruLPmvjnb13O3dLHEXEnhCCFFAl6+nMuv3I6zafQrP8q7MGNiWgb6ywbOjk8ATQoh8yjApVoefYubmI1xNTufFLg15278pVcqVsXfTRD7I01QbKkgBWEsen7PgqhRfFaJw9Kcu02/+dt7fsJ9mNSrxy5td+ejJVhJ2TkQCz4ayCsB6e3vfsRXX7QVgLXF8VrDlLLh6e/HVwjAajcyYMYPg4OBbXsu5G4oQxcWlaymMCdlH/wU7uHA1mS+HtGd14H00r1nZ3k0TBSSBZwf5LQBb1OMLUoi1IHLb4SZrqzAfHx98fHzy7K0K4egyTIpvd57k4Vlb+UF/hhEPebNlVHeeal9HtgRzUvIMzw4GDhzItGnTsr++W6WBrI2ls76XtXlyfo7X6/UEBATcsgFzbrL2qgwNDWX8+PG4u7sTFhZGQEAAW7ZsITEx8ZaadvDvfpu3y9oDM6vtuZUfEsIZRMYmMnHjAQ6eu8IDjb2Y3LcVjatXsnezRBHZJfA0TfMGpgNBSqmwHK+7A4GAAfAGwpRSRarmOX33dA4nHi7KR+SpuWdzxnYq2HBeXgVg9Xo9oaGhBAUFZQ8X5gydvI6/vRDr3eRWoNXf3/+uBV/vJmclg5CQkHu+VwhHdfFaCtN/O8y6yDhqVSnL/Gd8ebyNbPJcXNg88DRNy/qXMLd/jdcBI5RShsz3hmqaFqCUMtqqfbaSVwHYNWvW0KtXL8BcIuj2Icz8FpDNy70KtBam4KvRaMRgMEghVuFUMkyKVf/EMnPzEW6kZvBKNx/e6NGYCm4yCFac2PxuZvXoNE27pQJpZu/OOyvsMhkAf+DeRePuoaA9L1vJbwFYax2f5V4FWgsTWtOmTWPw4MGFbo8QtqY/dZkPf9zP/jNX6OLjxcdPyfClPaRlpFHGxbozXh1p0ooOMN72mhHoZfOW2Mi9CsAOHjw4z0Ks+Skgey/5LdBaEGFhYdK7E04h8XoqY0Oi6L9gBwlXU5j3TAe+e7mzhJ2N7UvYx5t/vMnILSOtfi5HCjx3IPG21y4Bd3Y9nFRWAdis2YsBAQHZvbOcBWDB/EzMz8+PkJCQ7BApyPE5i6ze7ff+/v7ExMRkT1yBWwu0zpgxA4PBQGhoKOvWrcse8swqwLpu3bpbjs0iRVeFIzOZFKv+OUWPz7eyXh/H8AcbsWVUd/q0rS3P6mxEKcW2uG28uOlFnvv1OfQX9PhW9yXDlGHV89qtAKymaaHA9BxDnAOB8UopvxzvGQN0VEoF3H58gwYNVLVq1bK/DgwMvKNSuBBC5LT/TBIfbNzP3tNGOjXy5JN+rWlaQ3p0tpJuSuf3k7+zdP9Sjl4+Ss0KNXmh5Qv0b9Kf8mUKVyswODj4ljXBkZGRsUqphrm915GeyBox9/Jy8uLOXh8A1apVy3U9mBBC3C7pZhpf/H6EFbti8azgxuzB7egn6+lsJjk9mR+P/8jyA8s5c+0M3lW8+eSBT3i80eNFfm53e2dH07SLd3uvIwVeBHcOX7oDoXe+VQgh8qaU4se9Z/nkl0MkXk/h+fsa8O4jzWQ7MBu5knqFtUfWsuLgChKTE2lbrS1jO46lW71ulNJs/0TNYQJPKWXUNC1C07ScMzV1gGNOsxRCOLTjF67ywcb97DIk0q6eO18P7UjrOlXs3awS4eLNi6w4uIK1R9ZyLe0aD9R+gGFthqGrobNrr9oe6/B8MS810AFjMwMuawA2AAjUNM2Aubc3vDiuwRNCWM/N1Azm/nGMxdsMlCvjwif9WvOfTvVxkdI9Vnf6ymmWH1jOj8d/JF2l80iDRxjWZhjNPZvbu2mAfdbh6QE9cMdGi5nhJhswCiEKZcuh83z0vwPEXb7JAN+6jH+8OVWl8rjVHUk8wtLopWyO3YyL5sJTjZ9iaKuh1K9c395Nu4XDDGkKIURhnTXeZPJPB9h84DxNqldkTeB9dPbO/w5BonD05/UsiV7CtjPbKF+6PC+0fIHnWz5PtfLV8j7YDiTwhBBOKy3DxNfbTzI77CgmpRjbuznDujbCtbQjLTEuXpRSbDuzjSXRS9hzYQ8ebh683v51hjQfQhU3x35GKoEnhHBKkbGXeX9DNIfjr9KzeXUm9W1FPc/CreUSecswZfB77O8sjV7KkctHqFmhJuM6jaN/k/6UK13O3s3LFwk8IYRTMd5IZfqmI3y/+xS1q5Ql6Hk/HmlZo/itqTsTCf83G/rOhXIedmtGakYq/4v5H8v2L+P01dM0qtKIKQ9M4YlGT1h970tLk8ATQjgFpRQb9pxh6i+HMN5M4+WujXinV9PiV9HgwiH44xM4/DOU94ILh6HB/TZvxo20G4QcDeGbA99w4eYFWnm1Ynb32fSo38Mua+gsoZj9TRFCFEeGhGt8sHE/O2Iu0b6eO98+3ZpWtR37eVGBJRpg62cQtRbcKsHD78N9I82/t6GklCRWHV7Fd4e+IykliU41OzGl6xTur3W/0/eiJfCEEA4rJT2DhVtjWPBnDG5lSjGlX2ue7VSfUsVpTd2Vs/D3TNB/C6XKwANvwgNvQ3nb7pufcCOBFQdXsObIGm6k36B73e4MazOM9tXb27Qd1iSBJ4RwSDtiLvLBhv0YLl7nyXa1mdinBdUrlbV3syzn+iXYPht2LwZTBvgNhYdGQ6WaNm1G3NU4vj7wNRuObSBdpdO7YW+GtRlGU4+mNm2HLUjgCSEcSuL1VKb+coj1+jjqeZbjm5c60a2pY67rKpSUq7BzAeyYC2nXoe1g6D4OPBratBkGo4El0Uv49cSvaJrGUz5P8VLrlxxusbglSeAJIRyCUor1+jNM/eUgV5PTebW7D2/0aEI5Vxd7N80y0pIhYils+xxuXIIWT8LDH0B12267deDSAZZELWHLqS2ULV2WZ1o8wwstX6BGhRo2bYc9SOAJIezuxMXrvL8hmh0xl/Br4MGnT7ehWc1iUqcuIx32fgd/TYcrZ8C7O/T8EOr45XmoJUWej2Rx9GK2n9lOpTKVGN52OM+1eA6PsvZb8mBrEnhCCLtJTTcR9FcMc/88jlvpUkx9ujX/6VhMJqWYTHBwI/w5FS4dhzo6eHoRNHrIZk1QSrHj7A6Co4LRX9DjWdaTt3zfYkizIVR0rWizdjgKCTwhhF1EnExk/A/RHLtwjSfa1uKjPi2pXrkYTEpRCo5vgS2TIT4KqrWAIaug2eNgo2n9JmXiz1N/sjh6MQcuHaBG+RpOtyuKNUjgCSFs6kpyGtN/O8x3/5yijns5lr6go2eLYvL86PRuCJsMsf8H7vXh6SBoEwClbPMcMt2UzuaTm1kSvYTjxuPUq1SPSfdPoq9PX6fbFcUaJPCEEDahlGLzgXg+/PEAF6+lMKxrI94tLjulnD8If0yBI79ChWrw2AzwexFK26Y0UVpGGj8ZfmJJ9BJOXz2NTxUfpj04jd4Ne1O6VDH487UQ+ZMQQljduaSbfPjjAUIPnqdlrcoseUFH27ru9m5W0V2Oha3TYN9q844oPT6AziPBzTbPx5LTk9lwfAPL9i8j/no8LTxbMKf7HB6u/7DTbv9lTRJ4QgirMZkUK/+JZcamI6SbTIx/zFy+p7SLk/9jfC3BvLwgYilopaDL69D1XZvtjnIj7QZrj6zl6wNfcyn5Eh2qd+DD+z6ka52uTr/9lzVJ4AkhrOLo+auMWx+F/pSRro2r8unTbajv5eTle5KvwM75sHMepN2EDs9Bt7FQpY5NTn8l9QrfH/qelYdWYkwxcl+t+5jZdia6GjoJunyQwBNCWFRKegYL/oxhwdbjVHArzecB7ejvW8e5/0FOT4GIZeY9L29cghZ9ocdEqGab7bcuJ19m5aGVrDq0imtp13io7kMEtg2kXbV2Njl/cSGBJ4SwmMjYRMauj+b4hWv0a1+biX1a4lXRNhM3rMKUAdHr4I+pkHTKvIbOf5LNFo1fvHmRbw58w5oja0hOT8a/gT+BbQNp7mnb3VmKCwk8IUSRXUtJZ+amw3y7K5baVcqxfGhHHm5W3d7NKjyl4FgohE2CCwegVjt4cg749LDJWrr46/Es37+c9cfWk2ZKo3fD3gS2DcTH3cfq5y7OJPCEEEXy5+ELvL8hmnNXknnh/oaMfrQZFZ15qcHpcAj7CGK3g0cjGLgMWj4Npaw/0SbuahxL9y9l4/GNoOBJnycZ1mYYDSo3sPq5SwIn/lsphLCnxOupfPzTATbuPUuT6hUJeaULfg2ceF/GhKPm3VEO/wwVqsPjs8xr6WywYDv2SiyLoxbzs+FnSmml6N+4Py+1eYk6FW0zGaakkMATQhSIUor/7TvL5J8OcjU5jbd6NuHVh31wK+2kVQ2unDVXGt+zEsqUy6w0/qpN1tIZjAaCo4P57cRvlClVhiHNh/BiqxepWcG2NfFKCgk8IUS+nUu6ycSN+wk7dIF29dyZMaCt81Y1uGmE7V/CroVgSodOw+Gh96BCVauf+ujlowRHBfP7yd8pW7os/235X15o9QJVy1n/3CWZBJ4QIk9KKb7ffZppvx4izWTigydaMPSBRrg4Y1WD9BRzlfFts+DmZWgzCHq8b5MCrIcuHSIoKogtp7ZQoUwFXm7zMs+3fL5EleixJwk8IcQ9xV66zrj10ew0XKKLjxef9W/rnAvITSaIXvvvEgOfnuD/kXkGppXtv7ifoH1BbI3bSqUylRjZbiTPtniWKm5VrH5u8S8JPCFErjJMiuXbTzDr9yOUKVWKaf3bMKRjPedbQK4UxGyB0ElwPhpqtYen5poLsVrZ3gt7WRS1iO1ntlPZtTKvt3+dZ1o8QyVXJx0GdnISeEKIOxy/cJX3QqLYc8pIz+bV+eTp1tSq4oR11M7ugdAP4cTf5iHLAUuhVX+rLzHQn9ezcN9Cdp3bhbube4kuuupIJPCEENnSMkwE/23gy7BjlHdz4csh7enbrrbz9eoun4QtU2B/CJT3gt7TQfcSlHa16mnD48NZtG8Ru+N341nWk1F+oxjUbBDlyzjhEHAx5HCBp2maOzAox0sGpVSYnZojRIlx8OwVxqzfx/4zV3i8TU0m921NtUpOti3YjUTzfpe7F0Op0vDgKHjgbShb2WqnVEqxO343C/ctJPJ8JFXLVWVMxzEMbDqwRFcXd0QOF3hAoFJqRtYXmqZN1zQtQilltGObhCi2UtNNzP/zOPP/PI57+TIsfNaXx9rUsnezCibtpnl5wf/NhtRr5ioG3cdD5dpWO6VSil3ndrFo3yL0F/RUL1edcZ3GMaDJAMqWLmu184rCc8TAGwzMyPH1JcAb0NunOUIUX/vPJDF63T4Ox1+lX/vafPRkKzwqWHfYz6JMGebiq39OhStnoOlj5pmX1VtY7ZRKKXae3cnCfQvZm7CX6uWrM6HzBPo36Y+bi5P1iEsYRww8g6ZpkUBA5tdeSikJOyEsKCU9g7lbjrPwrxi8Kriy+L86erWsYe9m5V/2zMuP4Px+qO0L/YOhYVcrnlKx4+wOFuxbQFRCFDXK1+CDzh/wdJOncXVxoh8SSjCHCzylVICmaeuAGCBMKdXL3m0SojiJijPy3roojpy/ygDfunzYpyVVylt/v0iLORcFoRPBsNU883LgMvPMSytNrFFKsf3sdhbuXUjUxShqVqjJxPsm0q9xPwk6J6MppezdhltomjYQ8AQMQBBgBHre/gyvQYMGqlq1atlfBwYGEhgYaLuGCuFkUtIz+GrLMRb9ZaBqRVem9W9Dj+ZO1KsznoY/PoGoNVDO3VxpXPcSlLbOMOLtQVerQi2Gtx1OP59+lLHBhtIif4KDgwkODs7+OjIyMlYp1TC39zpU4Gma5g2MUEqNzfHaOswzNcfmfK9Op1MRERG2bqIQTik6zvys7sj5qwz0q8vEJ5yoV5ecBNu+ME9KAbjvFej6rjn0rOD2oKtdoTbD2w7nKZ+nJOicgKZpkUopXW7fc7QhTV8g/LbXhgPT7dAWIZxezmd1VSu6suxFnfP06jLSIGKZuZLBzURoOxh6fADu9a1yutyC7qP7P5KgK0YcLfDCMIdbSI7XdMA6+zRHCOeVcwamUz2rU8pcky70I0iMgUYPQa8pULu9lU4nQVdSOFTgKaWMmqYFaZo2BvOzO4BEpVTIPQ4TQuSQlmFeVzfvj+N4VHBl6Qs6erZwkl5dXAT8/gGc2gnVmsMz66BJL6tMSMlaXjB/33yiEszP6D68/0N5RleMOVTgAWQuQZBlCEIUwpH4q4xat5f9Z67wVPvaTO7bCvfyTjCT8PJJ2PIx7F9vrjbeZw50eB5cLP9PlFKKned2snCveR2dBF3J4XCBJ4QouPQME8HbDMwJPUalsqVZ9JwfvVs7QdXsm5dh2+fwTxBoLvDQGHjgTXCzfDUBpRT/xP/Dgr0L2HNhDzXK15DlBSWMBJ4QTi4m4Rqj1u5j72kjj7epyZSnWuNV0cF3/EhPNU9I+eszc+Xx9s+ai7BaaSuw8Phw5u+dT+T5SFkwXoJJ4AnhpEwmxdc7TjJj82HcSjtJZQOl4PAv5pI9iTHmmnSPfAI121jldJHnI1mwdwG743dTvVx1xncaz8CmAyXoSigJPCGc0OnEG7wXso9dhkR6NK/OZ/3bUL2yg29YfHYPbH4fYrdD1WZWnZCy98Je5u+dz65zu6harirjOo1jYNOBstdlCSeBJ4QTUUqxNuI0U34+BMD0AW0YpHPwKuRJceYJKVFroHxVeOIL8H3BKhNSohOimb9vPtvPbMezrCejdaMZ1GyQlOkRgASeEE7jwtVkxq+PZsvhC9zn7cnMge2o5+nAhUVTrsL2L2HHXPNQZtd3zDukWKE23cFLB1mwdwF/xf2Fu5s77/q9y+Bmg6XwqriFBJ4QTuC36HNM2BDNjdQMPuzTkhe7NKRUKQft1ZkyYM9K876X1y9A64Hmkj1W2CHlSOIRFuxdwB+n/6Cya2Xe7PAmz7R4hgplKlj8XML5SeAJ4cCSbqYx6X8H2LDnDG3rVuGLQe1oXN3yU/YtxrDV/Jzu/H6o2wmGrIJ6HS1+mhhjDPP3zic0NpRKZSrxavtXea7Fc1RydeA/G2F3EnhCOKjtxy/y3rp9nL+awtv+TXjt4caUcSll72bl7uIx8w4pRzeZe3IDl0Orpy0+IeVk0kkW7lvIbyd+o1zpcgS2DeS/Lf9LFbcqFj2PKJ4k8IRwMMlpGUzfdJjl20/iXa0CP4zsQrt67vZuVu5uJMJf0yF8CZQuB/6ToPNIKGPZGaNxV+NYtG8RPxl+ws3FjaGth/JiqxfxKOth0fOI4k0CTwgHEh2XxDtr93L8wjVe7NKQsb2bU87Vxd7NulNGmjnktn4GKVfMsy4ffh8qVsv72AKIvx5PcFQwG45toJRWimeaP8OwNsOoWq6qRc8jSgYJPCEcQHqGiUV/xTAn7BheFV1ZMawTDzaxbHhYhFJwdDP8/j5cOg7eD8Ojn0KNlhY9zcWbF1kSvYS1R9aiUAxoOoDhbYZTo4KTbIItHJIEnhB2FnvpOu+u3Udk7GX6tK3FJ/1aO+aGz+cPwuYJYPgTvJrAM2uhySMWfU53Ofkyy/cv5/vD35NmSqOvT19eafcKtStaZ8sxUbJI4AlhJ0op1kXEMfmnA5QqpTFncHueau+AW4Ndvwh/fgqRy8GtMvSeDh2HgQUrC1xNvcq3B79lxcEV3Ei7wRPeTzCy3UjqV7ZOsVdRMkngCWEHiddTmfBDNJsOxNO5kSdfDG5PHXcH2w0kPRV2B8NfMyD1GnQcDt3HQXlPi53iRtoNVh1exfL9y7mSeoVeDXrxartXaezR2GLnECKLBJ4QNvb30QRGr9vH5RupTHi8OS939XasReRKmZcXbH7fvMFz417w6FSo1sxip0jJSGHtkbUsiV5CYnIiD9V9iNfav0ZLL8s+CxQiJwk8IWwkOS2Dz347zNc7TtKkekWWD+1Iq9oOtn7swiHYNN78nK5qU3g2xLzBs4WkmdLYcGwDQVFBXLhxgc41O/N6h9dpX729xc4hxN1I4AlhA4fjr/DW93s5cv4qL3ZpyLjHmlO2jAMtN7iRaH5OF7EM3Cpa/DldhimDX0/8yoK9C4i7Fke7au2Y1nUanWp1ssjnC5EfBQ48TdMaAt6Ae+ZLRsCglDppqUYJUVyYTIrlO04y/bfDVC5Xhq+HdqR7s+r2bta/MtIgfClsnWbe7Fn3Ejw8wWLP6ZRSbDm1hXl75hGTFENzz+bM7zmfB+s86HiTc0Sxl6/A0zStAzAC8AAMmb8SM7/tAzyiaZo7cBlYo5Taa/GWCuFkLlxJZnRIFH8fTcC/RXWmD2jrWJXIj28xD19ePGJeT9d7GlRvYZGPVkqx4+wO5u6Zy4FLB2hYuSEzu83kkQaPUEpz0O3RRLGXZ+BpmvYZcBwYq5RKyuO9VYBBmqYNVkqNt1AbhXA6YQfPM2Z9FDdS0/mkX2ue7VzfcXo0l2LME1KO/gYejWDI99DsMYutp9tzYQ9f6r8k8nwktSvUZsoDU+jj3YfSpeQJirCve/4N1DTtZaXUuPx+WGYgLs5x7JIitk8Ip5KclsHUXw6xYlcsLWtV5qv/tHec6gbJV+DvmbBrIZQuC70+hs6vQGnL9DoPJx7mK/1XbDuzDa+yXkzoPIEBTQbg6uKAi+hFiXTPwCtKYEnYiZLm0LkrvLV6D0fPX2P4g40Y/Wgz3Eo7wMQUkwn2rYKwyeb6dO2fhZ4fQSXLbNN1Mukk8/fOZ9PJTVR2rczbvm/zn+b/keKrwuHk1cNrKJNRhLg3pRTf7oxl6q+HqFy2DN++1ImHmjrIPpinw+G3MXBWD3U7wjOroY6fRT46/no8i/YtYuPxjbi6uBLYNpAXWr1AZVfLVzQXwhLyGlSfDgy2RUOEcEaXrqUwJiSKLYcv0KN5dWYMbEtVR5iYcuUchE2CqNVQsSY8HQxtAqBU0SeMJCYnsiR6CWsOr0Gh+E/z/0gFA+EU8gq8AE3TBgJhQCgQJjMwhTDbfvwi76zZi/FGGh892ZIXuzS0/8SU9BTzM7q/Z0JGKnR9Bx4cBW5Ff454LfUa3x78lm8OfENyRjJ9ffoyst1I2dhZOI28Ai8EWAN0BIYAMzRNU0gAihIsLcPEF6FHWfRXDN5VK/D10E60rO0Aw3hHf4dN48zbgTXtbS7b4+VT5I/N2gZscdRiLqdcpleDXrze/nW83b0t0GghbCevwAtXSq0H1kP2soOOgD//BuBlYIxSaqlVWyqEAzideIM3vt/D3tNGhnSsx4dPtqS8q52n21+KMa+nO7bZXLbn2fXQxL/IH5tuSuenmJ9YsG8B8dfjub/W/bzl+xatqrayQKOFsL28ZmnOvO3rJMy9u7Cs1zKHPMdpmnZCKfWHVVophAP4ad9ZJvwQDRrMe6YDfdraeSgv5RpsmwU754OLG/SakrnMoGjLAJRS/HHqD77a8xWGJAOtvVoz5YEp3FfrPgs1XAj7KPKPpkqpECBE07RpgASeKHZupmYw+acDrA4/jW99d74c0oF6nnaccq8U7F8Pv0+Eq2eh3TPg/xFUqlnkj959bjdz9HOIvhhNoyqNmN19Nj3r97T/s0khLMCSYzERFvwsIRzCkfirvL5Kz/GEa7za3Yd3ejWljIsdt8aK329eZhC7HWq1h0HfQL2ib8B88NJBvtR/yY6zO6hRvgYfd/mYJ32elN1RRLFS5L/NmqYdA8YCqujNyf5MX8zPCQ2Ap1Iq2FKfLUR+KKVYHX6aSf87QKXMtXUPNrHj2rqbl83VDMKXQFl3ePJL6PA8lCrawvZTV04xd89cNp3cRBW3KozWjWZI8yG4uTjA0gohLMwSP77NBB4BPrPAZ2WF3XilVEDm15GapkUopfSW+Hwh8nI1OY3xP0Tzc9Q5HmxSlS8GtadaJTsFgMkEe1bAlsnm0NMNs0g1g4s3L7Jo3yLWH11PGZcyDG8znKGth1LJ1UG2QRPCCizxDM/Sva/FQECOr3sqpYwWPocQuYqOS+K1VXrOGG/y3qPNGNnNx37VyM9Ewi+jzbuk1L8fHp8JNdsU6SOvpV5j+YHlrDi4grSMNAY0HcAr7V6RReOiRHCoAfrMEkPeSilDZk/PqJQy2LlZogRQSvHNjpN8+uthqlZ0ZU3gfegaWqYmXIFdvwRbJoF+BVSsbt4lpe2gIlUzSM1IZc2RNQRHBWNMMfJYw8d4vcPr1K9c33LtFsLB3TXwMtfc9VRK/VCYD9Y0rT/mhelXCnCYDkjMsbuLTtO0sUqpEYVpgxD5kXQzjTEh+9h84Dz+Laozc2A7PCrYYYd/UwZELoctU8zFWO9/DbqNhbKFX9RuUiZ+MfzCvD3zOHv9LPfVuo+3/d6mlZespRMlj6bU3eeaZIbeeOD3/K6x0zStJ+YJJ9MKGHZZa/rWAR5Zw5iapoUCQZnLH7I1aNBAVav27ySCwMBAAgMDC3I6Idh32shrq/TEJyUz7rHmDOvayD5T8E+Hw6+j4Nw+aPigefiyCMVYlVJsP7udOZFzOHL5CC08W/C239t0qd3Fgo0Wwv6Cg4MJDv73yVpkZGSsUqphbu+9Z+Blv8kcYgGYZ2JGYp49mZM35t6ZB+ZwKtR6PE3T/IF1SimPHK8FAdzey9PpdCoiQlZCiMJRSvH1jpN8+ushqlcqy7xnOtChvkfeB1ra9YsQ9hHsWQmVasGjU6FV/yINXx64eIDZkbP5J/4f6lSsw5sd3qR3o95SaVyUCJqmRSqldLl9L1/P8JRSW4AtmR/WE3OwZW2kZ8QcgOvyqoieD3dLMGMRP1eIbFeS0xizLopNB+Lxb1GDWQFtcS9v4yHM7OHLjyH1OnR5wzx8WYRNnk9fOc1Xe75i08lNeLh5MK7TOAY1HUQZlzIWbLgQzqvAk1Yyw88qlFJGTdPCNE3zzjFZRQcMt9Y5Rcmy/0wSr36n56zxJu8/3oKXH7TDEGZcBPzyrnn4stFD8NhMqN680B+XmJxI0L4g1h5dS5lSZQhsG8jQVkOp6FrRgo0WwvnlVQC2vR2qIQwHxmuadgnwAsbKGjxRVEopvvvnFB//dBCviq6sGXEffg1sPAvzRqK5Rp3+W/M2YAOWQusBhR6+vJF2g5WHVrJs/zKS05Pp36Q/I9uNpFp5Byk+K4SDyauHNxjYa4N2ZMucrDLWlucUxdv1lHQmbIjmx71n6da0GrMHt8fTlrMwsxaPh02C5CTz7Mvu4wo9fJluSufH4z8yf+98Em4m0KNeD97yewvvKlKuR4h7ySvwRmTO1AwFthR01qUQ9nbs/FVGfqfHkHCN0Y805dXujW27kPxclHn4Mi4c6neBJ2ZBjcItCVBK8Xfc38yOnE1MUgztqrXj8+6f06F6Bws3WojiKa/AMwCNgVcApWmakVuLv57MeqOmaT2kPJBwJBv3nGH8D9FUcHNh5bDOdGlsw91EkpPMe1/uDoZyntBvEbQbUujhywMXDzArYhYR5yNoULmBVDEQohDyCrw1WTXxNE1rBPQCAjP/WyVHAIZjLgwrgSfsLiU9gyk/H2TlrlN0bOjBvGd8qVG5rG1OnlW6Z/MEuHYBdC9Bz4lQrnBLHuKuxvGV/it+O/kbnmU9eb/z+wxoOoAypWTmpRAFle8CsEqpE0CwpmlVlFIzNU3zxrzA3B9zD7CRVVsqRD7EXb7Ba9/p2ReXxIiHvBn9aDPblfO5eAx+GQUn/jKX7vnPaqjjW6iPSkpJIjgqmO8Pf4+L5iIzL4WwgMLspakAMpcNBGf+QtO0RRZslxAF9tfRBN5avYeMDMWi5/zo3broBVHzJe0mbPsctn8JpcvB47PMPbtClO5JzUjl+8PfExQVxLXUa/Rr3I/X2r9GjQo1rNBwIUqWvJYljFZKzbr95bu8PdQyTRKiYEwmxdw/jjNny1Ga1ajEwuf8aFS1gm1OfiwUfh0Nl09C28HQawpUKng4KaXYdHITX+q/5My1MzxQ5wHe9XuXph5NLd9mIUqovHp4vTRNW6eUis3xWq57kSml1luuWULkj/FGKu+s2cufRxLo36EOU59uQznXohVFzZcrZ2HTODj4I3g1gRd+Mi8iLwT9eT2zImYRfTGaZh7NCOoVJHteCmEFeQYeYLhtcopPbm/UNK2yLFsQtrT/TBIjv4skPimZKf1a81zn+taftZiRDuGL4Y9PwJQOPT6ALm9C6YIXiI29EsvsyNlsObWF6uWr88kDn9DHuw8uRaxiLoTIXV6BNwMIwhx8vpgnp3hrmhbIv8sTQpVS+zDP3rx9+FMIqwiJjOP9DdF4VnBl7Yj7bbPx8xk9/Py2eUuwxv7migaeBV/sbUw2sihqEWsOr8HVxZU3OrzB8y2fp1zpcpZvsxAiW16BF5Q1OzPrhcyF6L34d3bmDE3TFOYNniXwhFWlppv4+OcDrNx1ivu9vZj7TAeqVix476pAkq+Ye3Thi6FCdQj4Glr2K/CautSMVFYdWkVwVDDX068zoMkAXm3/qlQbF8JG8lqWcCKX15KAkMxfOQNwnDUaKESW81eSGbkyEv0pIyMe8ua9R5tR2ppLDpQyP6PbNA6uxkOn4eYhzLJVCvgxit9jf2d25GzOXDtD1zpdGeU3isYeja3UcCFEbgqzLOEWWQGYuTBdCKsIP5nIyJV6bqSmM/8ZX55oW8u6J7wcC7++B8c2Q802MPg7qOtX4I+JSohiZvhM9ibspYlHE5mQIoQdFTnwcgjO+y1CFIxSihW7Yvn4p4PU8yzPquGdaVqj8DXj8pSRBrsWwNbPAA0emQqdXwGXgv2vcu7aOebo5/DriV/xKuvFpPsn0a9xP5mQIoQdWSzwLFD8VYhbJKdl8MHG/YRExtGzeXW+GNyeKuWsuKVWXCT89Cac3w/NnoDHpoN7vQJ9xPW06yyNXsq3B78FILBtIC+1fokKZWy0LlAIcVeW7OEJYTFnjTd5ZWUkUXFJvNmzCW/3bGK9KgfJV+CPKbB7MVSqBYNXQosnC/QRGaYMNh7fyNw9c7mUfIk+3n14y/ctalaw0W4vQog8SeAJh/OP4RKvrdKTnGYi+Hk/HmllxdA49LP5Wd3Vc5mTUiZC2coF+oh/zv3DjPAZHL18lA7VOzC3x1zaVGtjpQYLIQpLAk84DKUUK3fFMvmng9T3LM/qQD8aV7fS87orZ81Bd/hnqNEaBq+AuroCfUTslVhmRcxi6+mt1KlYh5ndZvJog0elZI8QDkoCTziElPQMPvrxAKvDT9OjeXXmDGlP5bJWeF5nMkHEUgibDKY08J9srkDukv9zJaUkERQVxPeHv8fNxY23fN/i+ZbP4+Zi5fWAQogikcATdnfhajIjV+qJjL3M6w835p1eTXGxxvO6C4fNk1JO/wPe3aHP7ALtlJJuSifkaAjz984nKSWJ/k3683qH12XhuBBOQgJP2FVUnJERKyIx3kiz3vq69BRz+Z5tX4BbpUJVH99xdgczw2dy3HicjjU7MqbjGJp7Nrd8W4UQViOBJ+zmx71nGBMSRdWKboSMvJ9WtQu2g0m+xO409+ouHoU2g6D3NKiQ/x7ZyaSTzIqYxV9xf1G3Yl3mPDyHHvV6yHM6IZyQBJ6wuQyTYtbvR1i4NYZOjTxZ+KwvXpbeDzP5CoRNMj+vq1Ifnltv3vA5n66mXiVoXxDfHf4ONxc33vF7h+daPIeri6tl2ymEsBkJPGFTV5PTeHv1XrYcvsAznesz6clWuJa28H6YR36Dn9+Fa/Fw32vQ431wzd/C7wxTBhuOb2DunrlcTr5Mv8b9eNP3TXlOJ0QxIIEnbObUpRsM+yYcw8XrTHmqFc/d18CyQ4PXLsBvY+DABqje0ryAvAD7X0aej+Sz3Z9xOPEwvtV9Wei/kJZeLS3XPiGEXUngCZvYGXOJV7+LxKTg25c68UBjC/aYlIJ9q2HzeEi9Dg9/AA+8BaXzN/x47to5Po/8nM0nN1OzQk1ZTydEMSWBJ6zu+92nmLhxPw28yrP0hY40rGrBfSWNp+CntyFmC9TtBE/Ng2rN8nXozfSbLN+/nGX7l6Gh8Wq7V3mx9YtSiFWIYkoCT1hNhkkx9ZdDLNt+goeaVmPeMx0st5g8ewH5JHMP77GZ0PFlKJX388Cs+nSfR3zOuevneLTho4zyG0WtilYuOSSEsCsJPGEVV5PTeOP7PWw9ksDQBxry/uMtLFes9eJx+N8bcGoH+PSAPnPAo0G+Dj2SeITPdn9GxPkImnk049Oun6KrWbAtxYQQzkkCT1jc6UTz5JSYhOtMfbo1z3bOXxjlKSMdds6DrdOgtBs8tQDaP5OvBeRJKUnM2zOPtUfXUtm1MhPvm8iAJgOkPp0QJYgEnrCoyNjLjFgRQWq6iW+GdqJrEwtNTjl/EH58Fc7ugeZ94InPoVLeVRQyTBmsP7aeuXvmciX1CoObDea19q9Rxc0Ki9yFEA7NoQNP0zR/wF0pFWLvtoi8/W/fWUav20etKmVZHdiRxtUrFv1DM9Lg/2bDXzPMZXsGLodWT+erV7fnwh6m/TONQ4mH0NXQMa7TOJp55m9CixCi+HHYwNM0zR0IAqbbuSkiD0op5v1xnM9Dj9KpkSdBz/nhUcECO5Kc2wcbX4Pz0dB6ADw2I1/bgiXcSOCLyC/42fAzNcrXYOZDM3m0oSwzEKKkc9jAAwYBYfZuhLi3lPQMxq+P5oc9Z+jvW4dp/dvgVrqIz8XSU+HvmfB/X0A5Txj8HbTok+dhaRlpfHfoOxbuW0iaKY3hbYbzcpuXKV+mfNHaI4QoFhwy8DKHMsOA/G+TIWzOeCOVwBWR7D6RyKheTXm9R+Oi96LO7jH36i4cgLZDzJs9l/fM87CdZ3cybfc0TiSdoFvdbozpOIb6lesXrS1CiGLF4QIvcyjTXSllkCEox3Xy4nVe+jqcuMs3+XJIe55qX6doH5ieCn/PMJfwqVgd/rMGmvXO87D46/HMCJ9BaGwodSvWZV6PeXSr161obRFCFEsOF3iAf34mqSQkJKDT/bt+KjAwkMDAQKs2TJhFxiYy/NtIlFJ8N7wzHRvm3QO7p7N7YeOr5l5du2eg96dQzuOeh6RmpPLtwW8JjgrGpEy81v41hrYeKlXHhShhgoODCQ4OzvnSXR/0a0op67conzRN8wWMSilD5tdBQKRSKvj29+p0OhUREWHrJpZ4v0Sd4521e6ldpSzLh3aiUVG2CUtPhW2z4O9ZUKEa9P0Kmj6a52E7zuzg092fEnsllh71ejCm0xjqVCxiD1MIUSxomhaplMp1NwlH6+F5ArocQ5n+gKemaeQWesJ2lFIs2XaCqb8ewq+BB4v/q8OzKDMx46Nhw0jzDMy2Q+Cxz/Ls1Z27do6ZETMJjQ2lfqX6LPRfSNc6XQvfBiFEieJQgaeUumVWpqZpvYBQCTv7yjApPv7pAN/sjOWJNrX4fFA7ypYp5EzMjPTMdXXTzQE3ZBU0f+Keh6RlpPHtwW8JigpCKcUbHd7gxVYvSjFWIUSBOFTg5aRpWiDmHp67pmmJsvjcPm6mZvDm6j2EHjzP8AcbMf6xFpQqVcjJRAlHYMMrcFZvXlf3+Kw8Z2D+c+4fpv4zlRNJJ2T4UghRJA4beJm9OunZ2dGlaykM+yaCfXFGJj3ZkhcfaFS4DzJlwK4FsGWKufJ4wNfm3VLuIeFGAjMjZvLbid+oU7EO83vO56G6DxXu/EIIgQMHnrCv2EvXeWHZbs4lJbPwWT96t85738pcJZ4wz8A8tQOaPQFPzjEvO7iLdFM6a46sYe6euaRmpPJKu1cY1noYZUuXLdz5hRAikwSeuENUnJGhy8PJUIpVwzvj16AQyw6UgsivYfP7UMoF+i2CdkPuuQdmVEIUU3ZN4XDiYbrU7sKEzhNoUNlClRaEECWeBJ64xdYjFxi5Uo9XRVe+eakTPtUKsQH01Xj48XU4HgqNusFT88G93l3fnpSSxJf6Lwk5GkK1ctWY1W0WjzR4RPa+FEJYlASeyLY+Mo6x66NoWqMSXw/tSPXKhRhG3P8D/PIupCXnWYVcKcXPhp+ZFTGLpJQknmv5HK+1f40KZYqwtk8IIe5CAk+glCLobwOf/XaYBxp7seg5PyqVLVOwD7l5GX59D6LXQR0/eDoIqja569sNSQY+2fUJ4fHhtK3WluBewVK6RwhhVRJ4JZzJpPjkl0Ms236CPm3Na+wKXO3AsNU8MeVqPHSfAA+OApfc/2olpyezOHoxy/Yvo1zpcnx4/4cMaDKAUlruvUAhhLAUCbwSLDXdxHsh+/hx71le7NKQD/u0LNgau7Rk2PIx7JoPXk3g5VBz7+4udpzdwSe7PuH01dP08e7DaN1ovMp5WeBKhBAibxJ4JdSN1HReWann76MJvPdoM17t7lOwSSLx0bB+OCQcgk6B4D8ZXHOvO3fx5kVmhs/k1xO/0qByA5Y8soTOtTpb6EqEECJ/JPBKoMvXUxn6dThRcUamD2jD4I4FqBtnyoCd88yLyMt7wXProbF/7m9VJn449gNfRH5Bcnoyr7R7hZfbvCwVDYQQdiGBV8KcS7rJ80t3cyrxBguf8+PRVgVYUJ4UZ94a7OQ2aPEkPPnVXbcGizHG8PHOj9Ff0KOroWPi/RPxruJtoasQQoiCk8ArQQwJ13h+6W6SbqbxzdBO3O9TgOdn+9fDz++Ye3h950GH53JdRJ6SkUJwVDDL9i+jQpkKfNzlY/o17idr6oQQdieBV0LsP5PEC8t2A7A68D5a16mSvwNTrsKvY2DfKqijgwGLwTP3nlp4fDiTd04m9kosT3o/yeiOo/EsW8TisEIIYSESeCXA7hOJDPs6nMrlyrBiWCe887t7yulw+OFlMJ6Ch8ZAtzHgcuf6vKSUJL6I/IIfjv1A3Yp1CeoVRJfaXSx8FUIIUTQSeMXc1iMXeGVlJLXdy7FyWGdqu5fL+yBTBmz7ArZOg8p1YOhvUP++O96mlGJz7Gam/TONpJQkhrYeysh2IylXOh/nEEIIG5PAK8Z+iTrH22v20KxmJb4Z2gmvivmYHWk8DT8EmqsbtB4Ifb6AsncOf8Zfj2fqrqlsjdtKS6+WLPJfRAuvFla4CiGEsAwJvGJqbcRpxq2Pwq+BB0tf7Ejl/GwVdvBH+N8b5h7eXaobmJSJtUfWMkc/hwxTBqN1o3m2xbOULiV/lYQQjk3+lSqGvt5+gkk/HeTBJlUJet6P8q553ObUG7B5AkQuh9odYMBS8PK5422GJAOTd0xGf0HP/bXuZ+L9E6lX6e5VEIQQwpFI4BUz8/88zszNR3i0VQ2++k+HvPfFPH8QQoZCwmF44C14+AMo7XrLW9JMaXy9/2sW7ltIudLl+OSBT+jr01eWGgghnIoEXjGhlGLW70eY/2cM/drXZlZAO0q73GNDZqUgYpm5Z+dWGZ7fAD497njboUuH+HDHhxxOPMwjDR5hfOfxVC1X1YpXIoQQ1iGBVwwoZa54sPT/TjCkYz2mPt0Gl3ttAn3TCD+9aX5m59MTnl4EFavf8paUjBQW7VvE8v3L8SjrwZzuc+jZoKd1L0QIIaxIAs/JmUyKD/+3n5W7TvFil4Z89GTLew81xkWYhzCvnDVv+NzlzTsKtO69sJeJ2ydy8spJ+jXux2jdaKq45XOhuhBCOCgJPCeWYVKM/yGKtRFxjOjmzbjeze8ediaTuYxP2CSoXBuGboJ6HW95y830m8zdM5eVB1dSs0JNgvyD6FJHFpALIYoHCTwnlZ5h4r2QKDbsOcObPZvwjn+Tu4fdjUTzps/HNps3fe47F8p53PKWiPgIPtzxIaevnmZws8G84/cOFcpUsMGVCCGEbUjgOaH0DBPvrN3HT/vOMvqRprzeo8nd33xqF4S8BNcT4LGZ0Gn4LWvrbqTdYI5+Dt8f/p66Feuy7NFldKzZ8e6fJ4QQTkoCz8mkZZh4a/Uefo2OZ9xjzXml253r5QDzLMwdc81DmO71YNjv5jV2OYTHhzNx+0TOXDvDsy2e5c0Ob1K+TO5FXIUQwtlJ4DmR1HQTb36/h00H4vngiRa8/OBd6svdSISNr8LR36BFX3hq3i3bg91Iu8FXe77iu0PfUa9SPZY/uhxdTZ2NrkIIIexDAs9JpKabeH2Vnt8PnufDPi15qWuj3N94JhLWvghXz0Hv6dB5xC1DmPrzeiZun8ipq6ekVyeEKFEk8JxAzrCb9GRLXnwgl7BTCsKXwKbxUKkWvLQZ6vplfzs5PZm5e+ay4uAKalesLc/qhBAljgSeg0vLMPHG9+awm9y3FS90aXjnm1Kuwc9vQ/Q6aPIIPB0E5f8tvLr/4n4m/N8ETiSdYFDTQYzSjZJenRCixJHAc2BpGeZndpsPmHt2uYZdwlFY+zxcPAo9PoCuo7IXkqdlpLEoahFLo5dStVxVWVcnhCjRHC7wNE3zBfwzv+wIBCmlwuzYJLtIzzDx9uq9/LY/ng/73GUY88BG+PE1KO0Gz/0APg9nf+v45eNM+L8JHEo8RF+fvoztNJbKrpVtdwFCCOFgHC7wAH+l1AwATdPcgROapvVUSunt2yzbyTAp3l27j1+iz/HBEy3unKCSkQ5bJpmXHdTtCAFfQ5W6gLle3YqDK/hK/xUVXSvKHphCCJHJoQIvs3c3HpgBoJQyapoWgbnHVyICz2RSvBeyj//tO8vY3s3vXHpwLcG8F+bJbdDxZXh0WnY5n7PXzvLB9g8Ijw/n4XoP89H9H+FVzssOVyGEEI7HoQJPKaXXNC3gtpe9AaMdmmNzJpNiwoZoftCfYVSvpozsftui8jORsOZ5uHEJ+i2E9s8A5moJPxt+5tN/PkWh+LjLx/Rr3E/q1QkhRA4OFXgAOZ/XaZrmDXgCa+3XIttQSjHppwOsDj/NGz0a80bP27YL06+AX0ZBxRrmJQe12wNgTDby8a6PCY0Nxbe6L1O7TqVupbq2vwAhhHBwDhd4twkCeiqljLd/IyEhAZ3u391BAgMDCQwMtGHTLEcpxbTfDvPtzlgCH/Lm3V5N//1meipsHm9eY+fdHQYsgwrmYcqdZ3fywf99QGJKIm/5vsXQVkNxKZVHhXMhhChGgoODCQ4OzvnSXStUa0op67eoEDRNGwPo7zZDU6fTqYiICBu3yjpmhx7lyy3H+O/9DZjct9W/Q5HXLsDaF+DUDujyBvScBC6lSc1IZY5+DisOrqBRlUZ89uBntPRqaddrEEIIR6BpWqRSKte9Eh2yh6dp2kAgLGtmpqZp3kopg52bZRXBf8fw5ZZjBPjVZdKTOcLu7F5Y/az5ed2ApdBmIAAxxhjG/D2Go5ePMrjZYEbpRlGudDn7XYAQQjgJhws8TdP8AWOOsHMHfIFiF3grd8Xy6a+HeaJtLT4b0JZSpTLDLjrEvL6ufFUYthlqtUMpxdoja5kZMZPypcszr8c8utXrZt8LEEIIJ+JQgZc5SSU08/c5v+WX6wFObOOeM0z8cT89mldn9qD2uJTSzFXJ//wEtn0O9bvAoG+hYjUuJ1/mox0f8efpP3mg9gN80vUTqpa76zC1EEKIXDhU4GUOWxb7ufRhB88zat0+OjfyZMGzvriWLgUpV+GHEXDkF/B9AR6fBaVd2X1uN+O3jScxJZH3dO/xXMvnKKWVsvclCCGE03GowCsJdsZc4tVVelrXrsySFzpStowLGE/B9/+BC4fgsRnQKZB0lcHCPXNZHLWYBpUbMK/nPFp4tbB384UQwmlJ4NlQdFwSw7+NoIFneb4e2omKbqXh9G5Y/Yx5+cGz66BxT85dO8fYbWPZc2EP/Rr3Y3yn8VLdQAghikgCz0ZiEq7xwvLduJcvw4phnfGo4GqenLLxVahcG15cC9Wa8uepP/lg+wdkqAw+e/AznvB+wt5NF0KIYkECzwbOJd3k+SX/UEqDFcM6U7OyG2ydDls/NU9OGbyS1LKVmL17OisPraSlV0tmPjST+pXr27vpQghRbEjgWdnl66k8v3Q3V5PT+T7wPhq5l4YNr0DUamj3H3jyS+JuJjD6t9c4cOkAz7Z4lnf93sXVxdXeTRdCiGJFAs+KbqZmMOybcE4l3uCboZ1o7WmCFf0h9v/g4Q/godFsOf0HE/9vIoCU8hFCCCuSwLOStAwTr63Ss/e0kQXP+nK/5zVYGgCXT0D/JaS17seciFl8e/BbWnm1Yla3WbLpsxBCWJEEnhUopZjwQzR/HL7A1Kdb09szHpYMgowUeH4D56s14b3Nw9hzYQ9Dmg3hvY7vyRCmEEJYmQSeFXwRepR1kXG81bMJz3oeg+X/hfJe8OLP/JNuZMzPg7iZfpMZD83gsUaP2bu5QghRIsiWHRa2clcsc/84zpCO9Xi7ajh8Pxi8vFHDfmdJ/DYCQwNxd3Nn9ROrJeyEEMKGpIdnQaEHz/Phj/vp0awaU6tvQftxEnh359rTi/ggYjpbTm3h0YaPMrnLZCqUqWDv5gohRIkigWche05d5o3v9bStXYmgGj/gsmUhtB6I4eExvBUWyOmrp3lP9x7Pt3z+9o2xhRBC2IAEngXEXrrOy99EULOiC6trrKDM7rXQ+RW2tOjFhE3/pWzpsix+ZDEda3a0d1OFEKLEksArosvXUxm6PBwXUwq/1PiGsgd/x9R9AouqVGThX+/Q2qs1sx+eTc0KNe3dVCGEKNEk8IogJT2DESsiSbx8mb/rB1Ph5A6uP/opE24e5Y+oP+jr05cP7/8QNxc3ezdVCCFKPAm8QlJKMTYkisMnT7O11jwqx+8n7vHpvHFuM4YkA2M7juXZFs/K8zohhHAQEniFNCfsGH/tPcyWqrPxTDIQ0XsS7xi+I0NlsNB/IV1qd7F3E4UQQuQg6/AKYcOeOFZuieDXKtOpmhzLhh7vMPzoctzd3Fn1+CoJOyGEcEDSwyugiJOJzAjZxsaK06iekcDs+59hecxq7qt1H7O6zaKKWxV7N1EIIUQuJPAK4HTiDcZ/G8Zq10+oqiXyXvtHCI0LY1DTQYzrPI4ypcrYu4lCCCHuQgIvn64mpzFqeSgLTZOp4JrIy807En1xL6N1o/lvy//K5BQhhHBwEnj5kGFSTPjuL6YkvY9Lucv817sFF26cZXb32VK/TgghnIQEXj7M+TmcwNhR3CyfyPB6DUCls/TRpbSr1s7eTRNCCJFPEnh52PjPUbpFvEZixQTeq1ULr7LuLOq1iAaVG9i7aUIIIQpAliXcw76TF/D85WXiK53h7RrVaODuzYrHV0jYCSGEE5Ie3l1cSLpB/DcvYqwSw9SqXvjV8OWrHl9RybWSvZsmhBCiECTwcpGSls4/i17hQuUo5nh60q1uN2Z1m0XZ0mXt3TQhhBCFJEOat1FKEbp0InFufzLH04PHGvZm9sOzJeyEEMLJSeDd5u8fl2BI/o4FHu709e7LtAc/kwXlQghRDEjg5XAocit7T05lsUcVnvZ+iildp+BSysXezRJCCGEB8gwv08VzJ9n0dyDLPCvRp25vJnX9mFKa/DwghBDFhcMFnqZp7kAgYAC8gTCllN6a50xLucnKNf1Z5ulGd8/7+OThzyTshBCimHG4wAPWASOUUgYATdNCNU0LUEoZrXXC4OUBLPVIo4NrY2Y/sVCGMYUQohhyqG5MZu/OOyvsMhkAf2ud8/sf3mdJ2ZM0yajC4kGrKV3KEX8GEEIIUVQOFXiADjDe9poR6GWNk+3a+xvzjBupnu7C4kH/w83FzRqnEUII4QAcLfDcgcTbXrsEeFr6RBkZ6UwLHwPAtG4L8apo8VMIIYRwII44fpev5ElISECn02V/HRgYSGBgYL5P4uJSmhGtxnDzZhK+TboUvJVCCCHsLjg4mODg4JwvVb3bezWllPVblE+apvkDQUopnxyvTQfclVIjcr5Xp9OpiIgIWzdRCCGEA9M0LVIppcvte442pBnBnT08dyDU9k0RQghRnDhU4GUuPYjQNM07x8s6IMw+LRJCCFFcOOIzvAAgUNM0A+be3nBrrsETQghRMjhc4GWG2wx7t0MIIUTx4lBDmkIIIYS1SOAJIYQoEUp84N22fqPYketzbnJ9zk2uz7FI4DnZDSsouT7nJtfn3OT6HEuJDzwhhBAlgwSeEEKIEsGhthYrCE3TEoBYC3xUVeCiBT7HUcn1OTe5Pucm12d7DZRS1XL7htMGnhBCCFEQMqQphBCiRJDAE0IIUSI43NZitqBpmjsQCBgAbyBMKaW3a6MsSNO0MYAXsAbzfqQBt5dXcjaZG4pPx1w+KizH6+4Ug3t5j+srFvdS0zRfwD/zy47kuM7icA/zuD6nv4eZ1+eJuXqNN4BSakbm99xxlvunlCpxvzCXG/K+7Wt3e7fLgtc3Bric+Wuds18b5n9I/IFIwL+43cs8rq9Y3EtgTI7fu2dej28xuof3uj6nv4c5ryfza+WM96/EDWlm/jTirZQy5HjZwL8/nRUHRqWUR+avAOXk1SaUUmHK/NNyYs7Xi8u9vNv1ZXL6e5nZOxif9XXmNUQA/sXhHt7r+jJfcvp7CPRUmb22zHsGYHS2+1fiAg9zfT3jba8ZgV42b4mVaZrme1ttweJG7qUTyPyHMuC2l70x3yunv4d5XF+2YnAPswwCQjJDzqnuX0kMPHfu/En6EndWWndqmqYNxPyTlq+madPt3R4rcUfupVNQtz6X9MZ8j9ZSTO7hPa4v6zWnv4eapnlrmhYI9FJKZQW8O050/0pi4IGD3gxLUUoFK6VClFJGpVQIMFDTNIccYrAAuZfOJwjzEJkx8+vidg9vub7icg+VUgalVDAQqmnauhzfcpr7VxIDz4j5p5KcvMj9+YlTynymkJMeBx1iKCIjci+dSuaMxek5hsiMFKN7mMv1Fbt7mBl6/pnXasSJ7l9JDLwI7vyJxB3zzCKnl/k/15bbXnYHYmzfGquTe+lEMof1siboZA39FZt7mNv1FYd7mPns8fb2GgAfnOz+lbjAy5pBddvDYx0QlvsRziXzJ8uxt73sTY7nCcWF3EvnkTmEZ7xtpp9vcbmH97i+4nAPjdx5P7yBUGe7fyVyL83bFkp6AhHKURdKFkKORbBGzD+FrXHm68txPeMx/0S5LnNYpVjcyzyuz+nvZeY/hrn1aPyUUnpnv4f5uL7icA/9yVxwDvgBkc74/2CJDDwhhBAlT4kb0hRCCFEySeAJIYQoESTwhBBClAgSeEIIIUoECTwhhBAlggSeEEKIEkECT4hiQNO0MZqmrdM0TWX+1yl35RfCmmQdnhDFRObWVuuUUpq92yKEI5IenhDFR0ccdEsnIRyBBJ4QxYc/DrpprxCOQAJPiGIga7NipIcnxF1J4AlRPOggu8KCECIXEnhCFA+9kN6dEPckgSdE8SDP74TIgyxLEKIY0DRNkVl/zd5tEcJRSQ9PCCejaVqgpmkxOb4eSI5q20KI3EngCeF8/ICxOb4ef9vXQohcyJCmEE4mc9sw38wvOwKhSimZsCJEHiTwhBBClAgypCmEEKJEkMATQghRIkjgCSGEKBEk8IQQQpQIEnhCCCFKBAk8IYQQJYIEnhBCiBJBAk8IIUSJIIEnhBCiRPh/0Botn3QwQ5AAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 504x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"figsize(7,5)\n",
"for pot,label in zip([iso,mwp14_v,mcm_v],\n",
" ['Isothermal','MWPotential2014','McMillan17']):\n",
" ro, vo= get_physical(pot)['ro'], get_physical(pot)['vo']\n",
" # The following Emax is in (km/s)^2\n",
" Emax= evaluatelinearPotentials(pot,4.*u.kpc)\n",
" Es= numpy.linspace(0.*u.km**2./u.s**2.,Emax,101)\n",
" # Setup Orbits to easily pass the required (z,vz) to the actionAngleVertical calculation\n",
" orbs= Orbit([[0.,0.,0.,0.,numpy.sqrt(2.*(E-evaluatelinearPotentials(pot,x=0.)))\\\n",
" .to_value(u.km/u.s)/220.]\n",
" for E in Es])\n",
" # The following action J is in internal units of ro*vo [different for McMillan17!] \n",
" J= numpy.array([actionAngleVertical.actionAngleVertical(orb,pot=pot).Jz() for orb in orbs])\n",
" # We now convert the action to Scott's units of 30 km/s x 0.3 kpc\n",
" plot(J*ro*vo/(scott_ro*scott_vo),\n",
" Es/(scott_vo)**2.,\n",
" label=rf'$\\mathrm{{{label}}}$')\n",
"xlabel(r'$J$')\n",
"ylabel(r'$H(J)$')\n",
"legend(frameon=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $\\Omega(J)$ vs. $J$\n",
"\n",
"Next, we compute frequency vs. action"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/Users/bovy/Repos/galpy/galpy/actionAngle/actionAngleVertical.py:194: RuntimeWarning: divide by zero encountered in double_scalars\n",
" return 1./_JzIntegrand(z,Ez,pot)\n",
"\n",
"/Users/bovy/Repos/galpy/galpy/actionAngle/actionAngleVertical.py:86: IntegrationWarning: Extremely bad integrand behavior occurs at some points of the\n",
" integration interval.\n",
" self._Tz= 4.*integrate.quad(_TzIntegrand,0.,zmax,\n",
"\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1690f32e0>"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcAAAAFCCAYAAABmXQ2gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABWcklEQVR4nO3dd1RUx9vA8e+lCVgo9oIF7F3AGhOjQuwdsHeFGDWmWhJ900xB84vGxCRgb1EUE3sDjVFjA9Zeo1jArrioYEG47x8LBJCmlN2F53MOx7A7s/ch5x4eZu7MPIqqqgghhBCFjYm+AxBCCCH0QRKgEEKIQkkSoBBCiEJJEqAQQohCSRKgEEKIQkkSoBBCiELJTN8B5JZSpUqpRYsWpXTp0voO5aVER0djY2NjNNfJyee8bN/sts+qXU7ev3PnjlHdU/l1P+XmtfLrnsqt+yk7beSe0u+1Un5OWFjYXVVV0/8frqpqgfhycXFRXVxcVGMzevRoo7pOTj7nZftmt31W7XLyvrHdU/l1P+XmtfLrnsqt+yk7beSe0u+1Un4OEKpmkDdkClTPunXrZlTXycnnvGzf7LbPql1O3zcm+fmzGNs9lVv3U3bayD2l32tl93MUtYCcBOPq6qoChIaG6jsUUYC4urrKPSVyldxT+UtRlDBVVV3Te69AjQC9vb31HYIoYOSeErlN7inDUaBGgPJXlRBCiJQKzQhQCCGEyC5JgEIIIQolSYBCCCEKJUmAQgghCiVJgEIIkUMajQZPT09cXFwIDg5+6f5arTbV98HBwdjZ2eVSdC8vODgYJycnvV0/v+glASqKYqsoykRFUTwS/3XOpO0aRVEc8yWwmyfg/pV8uZQQouBwdnbG3d0dV1dX3NzcXrr/6tWrU33v5uaGq2u6CxfzhZubG46O+fNrV5/0NQJcAwSqqhqoquoMwFdRFNsM2noAFxVFUVN8XcyTqLZNgTmNYfVQiJQtFUKI/OHn56fvEAqlfE+AiYnOUVXV8BQvhwMv/NmU2NZTVVUl6QtwBzzzJLheftByHFz8C+a3hwUd4PQGSIjPk8sJIQqm8PBwgoOD0Wg0TJo0Kfn1wMBAgoODk/8F3XSjVqslMDAQjUaT6nM0Gg0ajQYfH59Ur/v7+xMcHMyMGTOS+zo5OREcHIynpyfr1q3Dzs4OjUZDcHAwPj4+hIeHJ39WePh/v36Dg4MJDg5m0qRJL0zFFnT6qAbhCmjTvKZFl9gCU76oqqo25WuJCdFWVdWXn2TPDpuK8NZX0GYiHFkOB3+B1YPBrhq0eAeaDASLonlyaSHEy/li4ylOX3+Qp9eoW6EEn3Wr99L9/Pz8mDJlCra2tsnJRqPREB4ezsSJEwHw8fFJnjK1tbXFw8Mj3c9ydnbG1taW4OBg3Nzc0Gg0XLx4EW9vb1xdXZk0aRJ+fn74+flhb2+Pr68vjo6OzJ07F9BNZwYFBREYGMjEiRPx9PTEz88PX19fANasWZM8Av3222+TXy8M9DEFagtEpXntHmCfjb5TVFUNzLpZDhUpDi3GwPgj4LkYrEvC1o/hh7oQ/AU8uJHnIQghjFffvn1xcXHBx8cn+Vman58fzs7/LXdIGrFlJqlvyZIlk18LCAigZMmSyQk1KcFGRUXh7Oyc6tldyv5J/21vn/pXbcrpVxkB5o/sJLtUFEXxAEIyev/OnTupHhp7e3vn/Mw9UzOo10v3dfUQHPgJ9s2C/T9BA09oORbK1U/VJSrmGUXMTChapMCUWhTCIL3KyCy/ODo6EhYWljwlefFi+ssWoqJSjwXCw8NTJTBbW9sX+mi1Wtzd3ZOTaVBQEPBiYkvbP73PApg0aRLu7u7p9jdG/v7++Pv7p3ypVEZt9fFbWotuFJhSSV4cFaY1BWif0ZulS5fO2xPWKzfXfUWFw8HfdFOkx34Hxzeh5Xio3h4UhW+2nCEwLJJSxSyoZGdNZXvdl4O9FQ6J/13exgpTEyXvYhVC6JW/vz8TJ07Ew8ODqKgotFotnp6eaDSa5FWiISEhydONScknbQJMj6enJ0FBQcmfo9FoUo0sXzbOkiVLJk+t5vTzDEHawY+iKHczaquPBBjKiyNAWyAoow6Jz/6cE58J6pe9I3SeAW2nQOgiOOwPK/pA6drQcix9GranWqmiRETFEnE/liMR99l84gbxCf8dOm5molDRzorK9tYvJMnK9tbYWJmjKJIghTAWGo2GoKCg5MUvoFvw4ujoiKOjI7a2tri5uSW/r9VqU02P+vj44O/vn/x90vTmjBkz8PDwICgoiIsXLyY/M0x5HXt7e4KDgwkNDcXf3x9vb+9M+/v5+REaGkp4eDhubm74+voSHBycnISjoqKS+yd9XkGll2oQiqIEAT5JK0EVRQkD2quqqk3a85dylaiiKG7AGlVVM9wZqrdqEM+fwcm1cOBnuHUSipaGZt7gOhKK6ubtn8cncCP6CVejYrkaFUtEin8j7j8mKuZZqo8sbmmGQ1JiLGmNg91/o8eKdlYUMTPN/59TCCGMUGbVIPSVAG0Bb3TbH+zRlazXJL7ni26lp0+K9m7AJFVV3TP6TL2XQ1JVuPQ37P8ZLgSBmSU06q97TliqRqZdHz6JIyLqMRH3UyfHq4kJ8tnzhOS2igLlSljiYG+dIklaJf936eJFZPQohBCJDC4B5oX8ToCPnz/mbuxdHEo4vPjm7bNwcC4cC4D4p1CzE7QaB1Ve02Wwl5CQoHLn0VPd6PGeblo1efQY9ZibD56kam9pbpI8rZpy5OiQ+FVMFucIIQoRSYB54Ldjv+F/3J9BdQbh3dCbYhbFXmz06DaEzNd9xd6D8o11G+3r9QRT81yJ40lcPJH3/xs9Jo0cr0Y9JiIqlkdPn6dqX7KoBZXsk547/jdyrFKqKBVsLGX0KIQoUCQB5oE7sXeYc2QO6y+sx87SjvFNxtOrei9MTdJ5Phf3GI6thAO/wL1/oUQlaO4DLkPB0ibPYlRVFW1sXOJUauqR49WoWK5pH6danFOyqAVNKtvSpLIdTRxsaehgKyNGIYRRkwSYh07dO8WMwzPQ3NZQy64Wk5pNomm5puk3TkiAf3foFsxc3gsWxcB5CDR/G+yq5G/g/Lc4JyIqlot3HnE0IpojEfcJvxMD6GZra5YpnpgUdYmxeulimMgWDiGEkZAEmMdUVWX7le3MCp3F9ZjrtK/cng9dP8SheDrPB5NcPwoH5sKpP0BNgLo9dPsJK7nkW9wZ0cY+42iElqMRWo5c1XLk6n0ePNFNpRYvYkZDBxuaONjRpLItjR1sKVmsiJ4jFkKI9EkCzCdPnj9h6emlzD8xn+cJzxlUdxDeDTJ4PpgkOhIO+UHYEngaDZVb6laO1uoM6U2n6kFCgsqlezEcuarlaMR9jlzVcvbmw+Tp0yolrWnsYEsTB90osU75EliYSalJIYT+SQLMZ7djb/Oj5kc2XNxAScuSvOv8Lj2ceqT/fDDJ04egWQYHf4Xoq7oN9y3egcYDDPIA7thnzzkRGc2RCN0I8chVLbcfPgXAwsyEBhVtdEkxcepUFtiIgiw8PBxfX19CQ0MJCwtL9Z5Wq6VatWpMmTKFiRMnEhgYyOjRo/H29sbX15fw8HAmTZpEeHg4a9aswdHREY1Gw+jRo5kyZQrOzs7J70+ZMoWoqCguXrxI3759Mz2xRavVZnj8WXa8Sv+kyhMZHf2mD5IA9eTU3VP4hvhy5PYRatvXZmLTiRk/H0wS/xzObtTtJ7wWCpa24DpCt2imeLl8iftVqKrKjegnyVOmRyK0nLgWnbyHsUzxIolTpnY0q2ZHYwc7OQ5OFCjBwcHJVRZSHmeWVC0+ZVLw8fHBxcUl+ZSVpISXMnkmVX8A3ZFlYWFhqQ6utrOz49KlSxkmqZye4pJefycnpyyTm7u7e/L5pIGBgURFRREWFoanp2fyzxMYGIitrS1arRZHR8fkRK7VavH398fW1jbd2JP6vUzR4cwSoMxT5aF6peqxpOMSZr4xk+in0YzYPoIPdn9A5MPIjDslHcA9KhhGbIdqr+sO4J5VH/4cAzdP5t8P8BIURaGCrRVdGpZnate6rB3TipOfd2DDuNf4ons9WjmV5OzNh/huO0ufXw/gOj2IDwKOsun4dR48idN3+ELkir59+2aruK2npydr1qxJ/j4qKorw8PBU1RhSHk6d3kHV9vb2qer6pZXTIrvp9U87us2MRqPB0dERb29v/Pz88PTUlXHVarXJZ5l6eHikqpeY2SBGq9XmeuFgSYB5TFEUOlbryIaeGxjXeBz7ru2j+7ruzA6bTUxcTGYdoXIL6Lsc3tWA63A4vQ5+ew2W9oR/g3WnzxgwCzMTGlayZWirqszu14S/P25L2FQ3furfhDdrlWHXuduM+/0Izl8GMWDeQebvDefS3Uz+nwhh4Dw8PFKVOMrocGs3N7cXftl7eXmxevVqIHsHUkdFReHo6PhSRXbTFtINDg5+oXBuRv01Gg0uLv8t0suqkG5UVFSqhGVvb49Go2H16tU4OTklv25ra5t8jaTaiOkJDQ3F3T3Dw8BeiWzyyieWZpb4NPKhZ/WezDkyhwUnF7DuwjomOE+gu1P3zJ8P2jtC55nw5hQIW6xbNLOiD5Suo1sw09ALzIxjJWbJYkXo1qgC3RpV4Hl8AkcitOw8c5tdZ28xffMZpm8+g2PporSvXYZ2tcviWtUOc1P5O02kY+tkuHkib69RrgF0+u6luiQ9w3N2dk6e4kuPq6tr8jSnvb09np6e+Pr64u3t/UKZJCDVQdohISHs3LkzuR5gdorsZlRIN6mMXFLh3KSY0vZPW2swq0K6bm5uqaYqk+oVBgcHp0py9vb26f68KWk0GlxdXVMl89wgCTCflS1alq9bf02/Wv3wDfHl//b/HyvPrmRSs0m4lM1iC4S1Pbz+ge40maQDuDeMg51fQrPRuueEebixPreZmZrQtKo9TavaM7lTbSKiYtl55hY7z95myf4rzNt7iRKWZrSpVYb2tcvQpmZp7Ipa6DtsITLl4+NDQEBAliO4pGlQe3v75LZJ04TpcXR0TE4oSYnJx8cnVZ+kIrvpVZdPWUgXSDV9ml7h3ay8TCHdSZMmpZryTSs7hXhzsqAnI5IA9aRB6QYs67SMrZe2Mkszi2HbhvFWlbd43+V9KhWvlHlnMwto3B8a9fvvAO6/vtYdudbxO90zRCNccelgb82w16ox7LVqPHr6nH3/3mXnmVv8de42G49dx0QB1yr2tKujS4jVyxSTlaWF2UuOzPKLm5sbPj4++Pj4ZJoEvby8mDRpUqoE5urqyowZM3K0eCWjIrsZFdKFzJNLRtO42S2kGxgYmGrFatLil5TxZlYD0d/fP/l5Z0hICPfu3UsuM5VTMrekR4qi0NmxMxt6buCdxu+w99peeqzrwY+aHzN/PvjfB+gK8g4KBO/dULw8BA6HFZ5w/3IeR5+3ihUxo2P9csz0bMThT9xYN/Y1xrWtTsyz53y39Szus/bQZuZuvtp0mpPXoikoq5lFweDm5kZgYGCmbWxtbZMTUxJPT08CAgKyPdpJKrKbJCQkJHmUmLLIblLblEkvq+nEtP1TSllIN7PPS5oKdnZ2Tp6u9fLySrWSVKvVZvqHgre3Nx4eHnh4eODo6Ii7u3uuJD+QEaBBsDKzYkyjMfSq3osfNT8y/8T8VM8HTZRs/J1SoQmM3gWH58Gur2BuC3hzkm66NJcO3tYXExOFxg66U2c+eKsWN6Ifs+vsbXaeuc2yA1dYsO8StcsVx8OlEj0aV6R0ceN4HioKjqR9gBqNhokTJ6Ya1SVtBUhvW0Hfvn1TJREvL68XpgPDw8MJCAhIfgaYsv3LFNlNr5BuZoVzMyrS6+/vn61Cuq6urrRv3z7V+/fv3wd0WyWS4ki5CjQ4OJigoKDkZ6dpE6xGo0l+LzeSoOwDNEDH7xzHN8SX43eOU7dkXSY1nYRz2cyfJ6QSfQ22ToSzm6BMPeg2Gxya5Vm8+qSNfcbG4zcIDIvkWIQWUxOFN2uWpo9LJdrXKSPFg4Uo5GQjvBFSVZUtl7YwK2wWt2Jv0aFqB953eZ+KxSpm/0POboEtH8ODa7ptFO0/AyvbPItZ3y7cfshazTX+0ERy68FTbKzM6daoPB4uDjSqZCPPC4UohCQBGrHYuFiWnFrCwpMLSVATGFpvKKMajMLa3Dp7H/D0Eez+Fg7+AtaloOO3UL+PUS6Sya74BJV/LtxlrSaSbSdv8vR5Ak6li9LHpRK9m1SinI2lvkMUQuQTSYAFwM2Ym8zWzGZz+GZKW5VmgvMEujl1y97zQdBVn9j0Hlw/Ak7tocv/wL5aXoZsEB48iWPL8Rus1UQScvk+igKtq5fCw6USb9Uth5WFTJEKUZBJAixAjt05hu9hX07cPUG9kvWY1GwSTco0yV7nhHjdVomdX0FCHLSZBK3GG/0imey6fDeGPzSRrNVc45r2MTZW5vRvVpmhrapQ3sZK3+EJIfKAJMACJkFNYHP4ZmZrZnM79jadqnbifZf3KV+sfPY+IPoabJsEZzbqTpPpNlt37FohkZCgcvDSPZYduML2UzcxURS6NCzPyNbVaFjJVt/hCSFykSTAAio2LpZFpxax6OQiAIbVG8aI+iOy/3zw3FbdIpnoCHAZBm6fg5VdnsVriCKiYlm8/zIBIRE8evqcplXtGNm6Gu51y0m1CiEKAEmABdyNRzeYpZnF1ktbKWNVhgkuE+jq2DV7zweTF8n8qjtqreN3BX6RTHoePokjICSCxfsvE3n/MQ72VgxrVQ0v10oUtywcU8RCFERSDqmAK1+sPDPemMGyTssoY12GT/d9yqAtgzh6+2jWnYsUgw5fg/dfYOMAa0fC8t4QlXGZlYKouKU5o153ZPdHb/LrQGfKFrfkq02nafXtLr7adJqIqFh9hygMWHh4eHKNv7S0Wi12dnbMmDEj1/sHBwcnV1ZI+d8ieyQBFiCNyzRmRZcVfNP6G27F3GLw1sFM3DORmzE3s+5cvpGuBmGnmRARAr+0hD3fw/NneR+4ATEzNaFTg/IEjmnF+rGv0bZ2GRbvv0ybmX8xdoWGMzce6DtEYYAcHR3x9PTE0dHxhaPDwsPDsbe3T67YkBv9kxKdm5tbqpNecno6ilarZcaMGfj7+6d6LeVpLQWJJMACxkQxoZtTNzb22ohPQx92Xd1Ftz+7MffoXGLjshjFmJhCc28YdxhqvKU7Us3vdbhyIH+CNzCNHGyZ078J+ya1xfsNJ/acv0OnH/fyzoowzt18qO/whAHKbkHcnPZ/mcK0LyO9x0hJR5s5OTnh5OSU5WjWmEgCLKCsza0Z12QcG3puoK1DW3479hvd1nVj48WNJKgJmXcuUQH6LoP+q+BZDCzqCBvGQ2zmNbsKqvI2VkzuVJt9k9rxbrvq7Dl/l44/7mH8yiNcuC2JUPwnuwVxAwMDk4vNpjxEOjv90xamzUh6BWszKoCbJKOCtPfv3+fixYtcvHiRefPmZTqaNSZyGHYBV6FYBWa0mUH/Ov3xPezLJ/s+YdXZVUxsNpFGpRtl3rlWJ6j6Ovz9HRz4RXe0WsdvoYFnoVskA2Bjbc4Hb9Vi+GvVmLc3nMX7L7Pp+HV6NKrAu+1r4Fi6mL5DLFR8D/tyNupsnl6jtn1tJjV7uem/rAriajQagoKC8PPzS55eTDnqy6p/2sK0GUmvYK2bm1uGBXAzkrJSQ2BgYKZtjY2MAAuJJmWa8HuX35n+2nRuxNxg0JZBfLrvU6KfRmfesUgxeGu6rtySXRX4YzQs6wn3LmberwCzK2rBxI612TuxLd5vOLL91C3cfvibD1Yf5fLdbJSxEgVaUkHcjAQEBODu7g7oSiKlnfLMqn92ZVaw9lUK4Gq1WsLDw/OkMK2+yAiwEDFRTOhRvQfuVdyZf2I+i04u4sD1A3ze6nPeqPRG5p3LN4SRQRC6EIK/0C2Scf8Cmr9dKEeDACWLFWFKpzqMau2I398XWXbwCuuPXqd3k4pMcKtBJbts7scUr+RlR2b5JbsFcfOqf5LMCta+ShL79ttv6du37yvHY4hkBFgIWZtb867zu6zosgKbIjaM3TmWaf9M4+GzLJ5nmZhCs9EwLgSc2sK2ybCyH8Tcy5/ADVTp4kWY2rUueye2ZUjLKqw/dp12//ubGdvO8ujpc32HJ/Qgs4K4ffv2zbIwbXYK6mYmuwVrX0ZwcHCBGv2BJMBCrW7JugR0DWBUg1FsuLiB3ht6s//6/qw7liivWyDT0Rcu7oLfXoNLe/M+YANXpoQln3Wrx+6P3qRLg/L8svsib87czarDV4lPKBgHToj0JRXETVod6enpmTx6S1kQF3TP1FxcXAgMDExOKi/TP2XR2Yz+283NjYsXL6ZaUJOyYO2MGTMIDw8nKCiINWvWpFokk/Rayr5JcqsSu6GQk2AEoCvC++m+T7n84DJeNb340PXD7B2pduMYrBmu2zjfZiK8MRFMZWYd4GiElq82nSbsyn3qlC/BtC51aFW9lL7DEqJQkaPQRLY8ef6En4/8zNLTS6lQrAJfvfYVTcs1zbrj04e6M0WPrYTKraDPPLCplPcBGwFVVdl0/AbfbT3LNe1j3OqU5ZPOtWXFqBD5RBKgeCmaWxqm/jOViIcRDKwzkAnOE7Ayy0a5oGOrYNMHYGYBPeZC7S55H6yReBIXz4J9l/jlrws8fZ7AkJZVmdC+BjbWcs6oEHlJEqB4abFxsczWzGbl2ZVUKVGF6a9Np3GZxll3vHsBAofDzePQzBvcvwJzqcCe5PbDJ/yw4zwBoRHYW1swtWsdejauiFJIV9IKkdckAYpXdvjGYab9M42bsTcZWncoY5uMpYhpkcw7PX8KwZ/DwV+gbAPwXASlauRLvMbi5LVoPl13kmMRWl6rXpLpPRtQrVRRfYclRIEjCVDkSExcDN+Hfk/g+UAcbRz5uvXX1C9VP+uO57bBujG6hNjle2jUv9DuGUxPfILK74euMGPbOZ7GJzCubXV82jhSxMxU36EJUWBIAhS54p9r//B/+/+Pe4/vMaL+CMY0GoO5aRbPsB5ch7Wj4co+aOAFXX+AIsXzJ2AjcevBE77cdJrNx2/gWLoo3/RqQAvH7J/QIYTImCRAkWsePHvAjMMzWH9xPTXtavJ166+pbV87804J8brSSn9/B3ZVwWMhVGiSL/Eak93nbjNt/Ukioh7j4VKJTzrXwb6ohb7DEsKoSQIUuW53xG6+OPAF2idafBr5MLLBSMxNshgNXv5Hd5boo9u6Y9RavCNTomk8fhbPnF3/Mm9PODZW5nzdqz4d65fXd1hCGC1JgCJPaJ9o+ebwN2y9tJW6Jevy9WtfU92ueuadYqNg/Tg4t1lXc7Dnr1BUNoendfbmAz5ac4yT1x7Qo3EFvuheD1trGQ0K8bIkAYo8FXQliK8OfMWjuEeMbTyWYfWGYWqSyUIOVYXD82DHp2Blr9s4Xy2Lw7gLobj4BH756yI/7foXu6IWfNe7Ae3rlNV3WEIYlcwSoJwFKnLMvYo7f/b4kzcd3mS2ZjZDtg3hUvSljDsoiq7y/KidugUxS7rDrq8hXg6OTsnc1IQJbjVYN/Y1Sha1YOSSUD5cfYzox3H6Dk2IAkFGgCLXqKrK1ktb+frQ1zyNf8q7Td5lUN1BmCiZ/J319BFsnQhHV+iK73oulinRdDx7nsCcnf/y698XKVO8CN/1aUibmqX1HZYQBk+mQEW+uhN7hy8OfMHfkX/jXMaZ6a9Nx6GEQ+adjv4OG9+DYmWg73Ko0Dg/QjU6xyK0fLjmGBduP2L4a1WZ3Km27BsUIhMyBSryVWnr0vzU7ie+eu0rzt8/T5+NfVh1dhUJakLGnRoPgBHbQE2AhR3g+Or8C9iINHKwZdP41gxrVZVF/1ym59z9XLidRR1HIUS6JAGKPKEoCj2r9+TPHn/SpEwTvj70Nd5B3tyMuZlxp4rO4P03VHTRbZfY9ok8F0yHpbkpn3evx4Khrtx68ISuP+1j5eGrFJTZHCHyi16mQBVFsQW8gXDAEQhWVTXDcsWKojgDbont7VVV9U/bRqZADZeqqgT+G8j3Id9jZmLGF62+wK2KW8Yd4uNg+6dw2E+3OtRjMRSVk1HSc/vBEz5YfYx9F+7SqX45vuvdUCpMCJGCwT0DVBQlCPBRVTU8xfeeqqpq02nrDExRVdUz8fswYHTahCkJ0PBdfXCViXsmcureKTxqejCx6cTMyywdWQGb3odiZaHfcijfKP+CNSIJCSrz9oYzc/s5yhQvwux+TWhWzV7fYQlhEAzqGWDi6M8xKfklCkc3wkvPPGBSiu/bZzZaFIarconKLOu0jOH1hxN4PpB+m/pxLupcxh2aDIQRW0GNhwUd4ERg/gVrRExMFHzaOPHHO62wMDOh/7yD+O+5KFOiQmRBH88AXQFtmte0gHvahimTpaIozoqiOKY3ShTGw9zUnA9cPsDP3Y8Hzx4wYPMAVpxZkfEv64ou4L1bd3bo2pG6qVF5LpiuhpVs2Ti+NR3qleWbLWfxWRYmewaFyIQ+EqAtEJXmtXtAenM2rkCUoigeJD4vVBTFL2/DE/mhVYVWrO2+lhYVWvDd4e8Yv2s8UU/S3haJipWBoRt0BXYP/AzLe+uOVBMvKG5pztwBzkzrWpddZ2/T7ad9nLwWre+whDBI+loFmt0HFLb8t0hGq6pqMLok6JG24Z07d3B1dU3+8vd/YZ2MMDD2lvb83O5nJjebzP7r+/HY4MHBGwfTb2xqDp1nQo+5cPUA+LeBmyfyN2AjoSgKI1tXI8CnBc+eJ9D71/2sklWiopDw9/dPlQuADE/WyPdFMIqiuAF+qqo6pXjNF7BVVdUnnbZrVFW1S/GaH0DatrIIxridizrHx3s+5nL0ZYbXH864JuMyri4RGQYBg+DxfejxMzR44e8hkejeo6e8F3CUvf/excOlEtN71sfSXDbOi8LDoBbBAKG8OAK0BYIyaJsebS7GIwxALftaBHQNoE/NPiw8uZAhW4YQ8SAi/caVkp4LNtY9F9wxTVdzULygZLEiLB7ejHfbVScwLJK+/ge5Gf1E32EJYRDyPQEmLmIJVRTFMcXLrkAwgKIojknvJbYNTqdtQP5EK/KTlZkVn7X8jB/e/IErD6/gsdGDjRc3pt+4eFkYsgGajoL9c2B5H3kumAFTE4UP3qqF32AXLtx6SNef9hF6Wf5fCWEIG+HtgdCkrQ1pp0MT205Bt1CmJBCU+CwwFZkCLVhuPLrB5L2T0dzW0MWxC1ObT6WYRbH0G2uWwuYPwcYBBqyGUlnUJCzEzt96iPfSUK5pH/NF9/oMaF5Z3yEJkacMbiN8XpAEWPDEJ8Tjf8Kf3479RoWiFfB9w5eGpRum3/jqQVg1EBKeQ99lUl8wE9GxcYxfdYQ95+8wsHllPutWDwszORVRFEyG9gxQiGwxNTFlTKMxLO64mHg1nqFbhzL/xPz0D9Wu3AJG74Ti5WBZLwhbkv8BGwkba3MWDWuKTxtHVhy6ypCFh9DGPtN3WELkO0mAwuA1KdOEwO6BtK/Snh81P+K9w5vbsbdfbGhXFUbugGptYOO7sGOqLI7JgKmJwpROdZjVtxGaK1p6/bKf8DuP9B2WEPlKEqAwCiUsSjDzjZl82epLjt89Tp8NfdgdsfvFhpY2uueAzbxh/0+67RJP5Rd7Rno1qcTvo5vz4HEcvX7Zz/4Ld/UdkhD5RhKgMBqKotCrRi9WdV1FuaLlGL9rPN8c+oan8U9TNzQ1022a7zQTzm+DhR0hOlI/QRsB16r2rBv7GmWKF2HIwsOsOnxV3yEJkS8kAQqj42jjyIrOKxhcdzArz65k8JbBXH2Qzi/t5t4wYA1or8C8dnAtLP+DNRIO9tasfacVraqXYvIfJ/hmyxkSEgrGAjkhMiIJUBglC1MLJjadyJy2c7j26Bpem7zYfnn7iw1ruOmeC5oVgUWd4dS6fI/VWJSwNGfhUFcGt6iC/55w3l11hCdx8gxVFFySAIVRa1u5LWu6rcHJ1omP/v6I6QenvzglWqYOjNqlqye4Zijs+R4KyPaf3GZmasKXPeoxpVNtNh2/wZCFh2WFqCiwJAEKo1ehWAUWd1zMsHrDCDgXkP6UaLHSupNjGnjBrq/gz7fh+dP0P7CQUxRdfcE5/Ztw9KqWPr/uJyIqVt9hCZHrJAGKAsHcxJwPXT/k53Y/cz3mOl6bvNh2aVuaRpbQ2x/aToXjq2BpD4iRVY8Z6d6oAktHNuPOw6f0/nW/lFUSBY4kQFGgtHFow5qua6huW52P93z84pSookCbj8FjEVw/AvPd4O4F/QVs4Fo4lmTtmFZYmJrQz/+gbJMQBYokQFHglC9WnkUdFzG83nACzgUwaMsgrjy4krpR/d4wdBM8fQAL3HRHqYl01ShbnLVjWlHR1ophi0LYfPyGvkMSIldIAhQFkrmJOR+4fsDc9nO5EXMDr43pTIk6NIVRwWBlD0u6w8k/9BOsEShnY8lqn5Y0rGTDuJUalh24rO+QhMgxSYCiQHuj0hsEdgukpl1NPt7zMV8d+Cr1lKi9oy4JVmgCgcNh32xZIZoBG2tzlo9qTvvaZZi2/hQ/7DgnVeaFUZMEKAq8ckXLsbDjQobXH87q86sZuHkgl6Mv/9fA2h6GrId6vSH4M9j8AcQ/11u8hszS3JTfBrng5VqJObsu8NmGU7JhXhgtSYCiUDA3MecDF92U6K3YW/Td1Jct4VtSNLCEPgvgtfcgdCGsGiBniGbAzNQE3z4N8X7DkaUHrvDhmmM8j0+nQocQBk4SoChU3qj0Bmu6raGWfS0m7Z3EFwe++G9K1MQE3L+ArrPgQhAs7gwPb+o3YAOlKApTOtXmo7dq8ueRa4xZoZFTY4TRkQQoCp1yRcuxoMMCRtYfSeD5QAZvGUzEg4j/GriOgP4Buu0R893g9hn9BWvAFEVhXLsafN6tLkGnbzFicQgxT2XqWBgPSYCiUDI3Mec9l/eY234u1x5do++mvuy8uvO/BjXfguFbID4OFrwF4bv1FquhG/ZaNf7n2YiD4fcYvOAQD57E6TskIbJFEqAo1N6o9Aaru62mSokqvPfXe3wf8j1xCYm/wCs01q0QLVERlveBoyv1Gqsh6+NSibkDnDlxLZqB8w5xP0bODxWGTxKgKPQqFqvIkk5L6FerH0tOL2Hk9pHcirmle9PWAUZuhyqvwbq3Yc9M2SaRgU4NyuM32IVztx7Sf95B7jyUs1aFYZMEKAS68kqftviUmW/M5FzUObw2eXHg+gHdm5Y2MDAQGvaFXdNh0/uyTSID7WqXZeHQply+F0M//wPcjH6i75CEyJAkQCFS6FitIyu7rsTe0h6fIB9+PfYrCWoCmFlALz9o/T6ELYKAQfAsRt/hGqTWNUqxZHgzbkY/oa//Aa5rH+s7JCHSJQlQiDSSKs53dezKL0d/YUzwGKKeROkO0nb7HDp/D/9uhyXdpJpEBpo7lmTZqOZEPXpGP/+DXJMkKAyQJEAh0mFtbs3Xrb/ms5afEXozFM+Nnhy9fVT3ZrPR0Hc53DoFC9whKlyvsRoq58p2LB3ZjPsxz+jnf4DI+1JTUBgWSYBCZEBRFDxqerC883IsTCwYvm04S08t1Z1/WbsLDN0Ij7Uw3x0iw/QdrkFqUtmO5aOaEx0bR1+/g1JYVxgUSYBCZKFOyToEdAvgjUpvMDN0Jh/s/oCHzx6CQzMYGQQWRWFJVzi3LesPK4QaOdjy++gWPHr6nH7+kgSF4ZAEKEQ2lLAowey2s/nI9SP+iviLfpv6cS7qHJSqrtsrWKomrOoPYYv1HapBql/RhhWjmvPwSRz95x2U6VBhECQBCpFNiqIwtN5QFnZYyOPnjxm0ZRDrL6yHYmVg2GZwag8bJ8Cur2WvYDp0SbAFDx7rkqAsjBH6JglQiJfkXNaZ1d1W07B0Q6b+M1V3oLaZOfRfCU0GwZ4ZsGGc7BVMR4NKNiwb2RxtbBz9/Q9yI1qSoNAfSYBCvIJSVqXwc/dLdaB2ZOwt6P4ztJkER5brpkRlr+ALGjnYsmxkc+7HPGPAvEPcfiCb5YV+SAIU4hWZmZjxnst7zGk7h8iHkfTd1Jc91/ZC208SSyoFw+KuslcwHY0dbFk8oim3HjxhwPxD3H0kx6aJ/CcJUIgcalu5LQFdA6hQrAJjd45ljmYO8c5DdXsFb5+WvYIZcKliz6JhTYm8H8ug+XKAtsh/kgCFyAUOJRxY1mkZvWv0Zt6JefgE+xBVtWXiXsH7upJK14/oO0yD09yxJPOHNCX8bgyDFhwi+rGUUhL5RxKgELnE0sySL1p9wZetvuTo7aO602OKWOj2CppZwaIuumlRkUrrGqXwG+zC+VsPGbbosBTVFflGEqAQuaxXjV4s67Qs+fSY5bcPoY7cASUd4fe+UlcwHW1rleGn/k04HhnNqCWhPImL13dIohCQBChEHkg6PaZ1pdb4hvgy8cgPxA4M/K+u4N7/yV7BNDrWL6+rLH/pHm8vD+PZ8wR9hyQKOEmAQuSREhYl+LHtj0xwnsCOKzvov9OH8K4zoYEn7PwStnwMCTLSSalnk4p826sBu8/dYcKqIzyPlyQo8o4kQCHykIliwqgGo/B390f7VEv/bUPY1qQPtBoPIfNgzTCIk31wKfVrVplpXeuy9eRNJv9xgoQEGSmLvCEJUIh80Lx8c1Z3XU0Nuxp8vHcivrbFiXvrKzizAZb31q0UFclGtq7Ge241CAyL5MtNp3UVOITIZZIAhcgnZYuWZVGHRQysM5DlZ5YzMjqU291mQ8RhWNgJoq/pO0SDMqF9DUa2rsbi/Zf5Iei8vsMRBdBLJ0BFUaoqitJOUZTeiV/tFEWpmgexCVHgmJuaM7nZZGa8MYOzUWfxvLCYkK7fQnSkbsP87TP6DtFgKIrC1C516OvqwE+7LjB/rxwmIHJXthKgoihNFEX5TVGUAMAHcAKUxC8n4O3E979VFKVxnkUrRAHRqVonVnZZiU0RG0afnMuitu+gJsTDwg5wZb++wzMYiqLwTe8GdG5Qjumbz7A6JELfIYkCRMlqbl1RlO+AC8AaVVWjs2hrA3gBjqqqTsm1KLPB1dVVDQ0Nzc9LCpFjMXEx/N8//8eOKztoX64lX50Po7g2AvrMh7rd9R2ewXj6PJ5RS0L558Jd5g5wplOD8voOSRgJRVHCVFV1Tfe9zBKgoiijVFWd/4oXfeW+r0ISoDBWqqqy7PQyZoXNomLRcsy6/5gakUeh80xoNlrf4RmM2GfPGTT/ECevPWDx8Ka0ql5K3yEJI5BZAsx0CjQnCSw/k58QxkxRFIbUG8KCDguIjX/KQMsYNlVvCVs+guAvZMN8ImsLMxYOa0q1UkUZvTSU45FafYckjFymCVAWtwiRf5IK7dYtVY8p8ZFMr92SZ/t+gHXvQLwcEg1ga23B0pHNsCtqwbBFIVy4/UjfIQkjltUiGN98iUIIAegK7c57ax7D6g0j4Ok1htdy5ubJAN0Zok/llz1A2RKWLBvZHBMFhiw4xHWtVJUXryarBOipKEq8oijbFUX5SFZ4CpH3zE3M+dD1Q3548wcuqo/xqladA9cPwOIu8Oi2vsMzCNVKFWXx8GY8ePKcoQsPo42VWoLi5WWVAAPRreo8AvQDNLmREBVFsVUUZaKiKB6J/zpn0naioii+iqI4K4ripiiK36tcUwhj417FnZVdVmJfrAI+5Urj/ySChAXucO+ivkMzCPUr2uA/xIUr92IZsTiE2GdSRkm8nKwSYIiqqmtVVZ2cuIrGDuhA6oR4T1GUkS953TVAoKqqgaqqzgB8FUWxzaS9N7AT3R7ESS95LSGMVjWbavze5Xc6VuvIT7ZFmWD5jOiFb0FkmL5DMwitnEoxp39jjkZoeWeFhjg5PFu8hKxWgc5M8320qqrBSQlRVVUTdElpjKIo7bJzwcRE56iqaspjHcIBtwy6aFVVtUv88lRVVZud6whRUFibW+P7ui+Tm01mn1UR+pW05uzvPeD8dn2HZhA61i/P9J66ChKT1h6Xc0NFtuX4LNDEUZwr4J7NLq6ANs1r2qz6J06BOr50gEIUAIqiMLDOQBZ1XMyzoqUYVMaOPzeOBM1SfYdmEAY0r8z7bjX5Q3ON77ad1Xc4wkjk5mHY2d2FbgtEpXntHmCfUQdFUTzQjRKdFUVJd2XqnTt3cHV1Tf7y9/fPZjhCGI/GZRqzunsgjcu68H+l7Pj8n894+tfXslcQeLd9dQa1qIzf3+Fybmgh5u/vnyoXABmemJDlUWhZURTlXxKfy6mq+kc22nsAU1RVdUnx2kSgqaqqntnofxHwUVU1OOXrchKMKEziE+L5WfMj808tos7TZ/xQ+nUqdf8NTM30HZpexSeojPtdw9aTN5nTvwndG1XQd0hCz175JJhsmgm8BWiy2V6LbhSYUkleHBUCuqnPNC9pyP50qxAFkqmJKRNcP2BO2x+JtCxK3/v72ft7d3gWo+/Q9MrURGFW38Y0q2bPh6uP8s+Fu/oOSRiw3HgG6K+q6tuqql7OZpdQXpzutAWC0jZMTH4702kr68CFANpWbkdAz/WUty7L2PirzF3WlviHN/Udll5Zmpsyb4grjqWK4bMsjJPXMj3DXxRi+V4QN3EVZ2iaBS2uQDCAoiiOSe+pqqrhxW0PjsDqfAhVCKPgUMKBZR5b6F7ald/MHvNOgDv3bxzRd1h6ZWNlzpIRzShhacbwxSFERMXqOyRhgDJMgIqi2CiK0vtVPzixWG6JDN72BDwSN8J7A6NTbG9Iu9cvNHEzvHfiAhjZCiFEGlZmVnzVeRGf1RpKqJmK19ZBnDhVuP9OLGdjyZIRzXgaF8/QhYeJipHTYkRqWZVDsgGmADtUVd2VrQ9UlPbo9vR9q6rqg1yJMhtkEYwQOqcubufDvz/ilonK5Gq98HrjSxRF0XdYehNyOYqB8w9Rr0IJfh/VAisLU32HJPLRK9cDTPEB7dGN2lQgDN2WhJQc0U1j2gF+2U2WuUkSoBD/iY66yJR1nuw1jaOrTW2mdVmCtbm1vsPSm20nbzBmhYb2tcvy2yBnzEzz/emP0JMcJ8A0H9Ye3UKUpGd4WnQJMTSrivF5SRKgEKklPH3IvNU9mRt/CydzG2Z1WUZV28J7lsSS/Zf5bMMpBjavzPSe9Qv1qLgwySwBvvSmIVVV067KFEIYIJMixfEZsJ0G60cwWRtKv/W9+PL1b3jLsYu+Q9OLoa2qcj36MX5/h1PexpJx7WroOyShZ9meB1AUpbGiKO3SLmxRFKVa4utVcz06IUTOmJrRqtcSVlfrh9OTx3y4dzIzDnxFXELhLLA7qUNtejauwPc7zhMYFqnvcISeZZkAExPfBXQb0IOB+4qihCSVQlJV9VLiMz85nl4IQ6QolGv7fyxu/iUDHsSw7PxqRm4exK2YW/qOLN+ZmCjM8GhEK6eSTF57nL3/3tF3SEKPMk2AiatA56M7eswksfqDPfAd8ImiKAEpRn738zRSIUSOmDfux5SuS5gRFcvZe6fw2tCbQzcO6TusfGdhZsJvg12oXqYYY5ZrOHVdNsoXVlmNAEcD7VM+90ssibRWVVUvdHX6XBRF6UP2j0ITQuhLtdfpNHATq6LBNiYK7x2jmXd8Hglq4aqjV8LSnMXDm1Hc0ozhi0K4pn2s75CEHmSVAI9ktrIzRTJMSohCCENXpg6Ow4NYGV+KDo9imHNkDuN3jSf6aeEaCZWzsWTx8GY8jotn2MLDRMcWzueihVlWCVBqrAhREJUoj/WwrfjaNOHTu1Hsj9yL10YvTt09pe/I8lWtcsXxG+zC5XsxeC8L5enzeH2HJPJRVgnQNj+CEELoQZFiKP1X0a+WF0uv3UCNvcvgrYMJOBtQqKqqt3IqxfeejTh0KYqP1hwnIaHw/OyFXVYJ0CmT8zxRFKWqoigfKYrSR1GUkFyOTQiR10zNoMsPNGgzldWXw2n+3ITph6Yzee9kYuMKzwHSPRpXZFLH2mw8dh3f7VJRvrDINAGqqjoT0CiKMjIx2ZVI3BbxkaIooehWh36vqupaZLQohHFSFHhtArZ9FjI38irjn5iw7dI2+m/uz0Vt4ak89nYbx+SK8ksPXNZ3OCIfZGcjvCswBt1xZ/fR1edzQrc6dAokb5comVdBCiHyQb1emAzbhHf0Q/yiHqGNvUP/zf3ZFL5J35HlC0VR+KJ7fdzqlOHzDacIOl349kkWNlkmQFVVtYnnqNkB1VVVLamq6pik1aGKolQDqqmqmrbIrRDC2Dg0g1HBtDCzY83li9SxLM2UvVP48sCXPI1/qu/o8pypicKc/k1oUNGG8Ss1HI3Q6jskkYeyfRRa4paHS+m8fklV1aO5GpUQQn/sHWHkDsqUd2HBiX0Mt6nHmvNrGLxlMBEPI/QdXZ6ztjBjwbCmlCluycjFIVy9V3iehRY2UhNECPEia3sY/CdmDfvywdGtzLGqQ+SjSPpu7MvOqwX/PPxSxYqwaHhT4lWVYYsOc1+K6RZIkgCFEOkzKwK9/ODNT2h7ejurnxTHoVgF3vvrPWaGzCzwB2o7lS7GvCGuRGofM3ppKE/iZI9gQSMJUAiRMUWBNydB7/lUitSw7PJF+lXpzNLTSxm+bTg3Y27qO8I81bSqPbO8GhN29T4frD4qewQLGEmAQoisNfSEoRuxeBLNp4cCmFlnFBe0F/Dc6MneyL36ji5PdWlYnk8712HLiZt8veWMvsMRuUgSoBAieyq3gFHBYF2Kjtu+YpXTQMpYl+Gdne/wo+ZHnic813eEeWZk62oMa1WVBfsuseifF9YCCiMlCVAIkX32jjAqCKq0ourWqaywqkuf6r2Zf2I+I7ePLLA1BhVFYVrXurxVtyxfbjrN9lMFe+q3sJAEKIR4OVZ2MGgtuAzD8p85fB55kW9afMaZqDN4bfJi/7X9+o4wT5iaKPzYrwmNKtny7sojaK5KCVRjJwlQCPHyTM2h62zo8A2c2US33T+xqs0c7C3teTv4bX468hPxCQVv1aSVhSkLhrpSzsaSUUtCuXw3Rt8hiRyQBCiEeDWKAi3HQv+VcO8CjgHD+N15Ej2r98T/uD+jdoziduxtfUeZ60oWK8Li4c1QE/cI3ntU8E/IKagkAQohcqZWJxixHUzMsFrSky9tnfm69decuncKz42eBXJKtFqposwf2pQb0U8YtTSUx88K3mi3MJAEKITIuXL1YfQuKNcA1gyle+Q5Vnb+PXlKdI5mToFbJepSxY4f+zXmaISWCauOEC97BI2OJEAhRO4oVgaGboSGfeGv6Tjt/Jbf31pIrxq9mHdiXoFcJdqxfnmmdanLjtO3+GrT6UJVSLggkAQohMg95pa649PafwYn12K1rDdf1H+bb1p/w5moM3hu9GRP5B59R5mrRrSuxojXqrF4/2UW7JM9gsZEEqAQIncpCrz+AfRbAXfOwby2dLOsQEDXAEpbl2bszrH8EPpDgTpLdGqXOnSqX46vt5xhy4kb+g5HZJMkQCFE3qjdBUbuABNzWNSZalfDWNF5BX1r9WXRqUUM2zaMa4+u6TvKXGFiojCrb2NcKtvxXsBRQi5H6TskkQ2SAIUQeadcffD+Cyo0gbUjsfx7BlObfcL3bb4nXBuO50ZPgq8E6zvKXGFpbsq8Ia5UsrNi1JJQLtx+pO+QRBYkAQoh8lbRUjBkAzgPgb3/g1UD6FC+Fau7raZK8Sq8v/t9ph+cXiAqztsVtWDJ8GaYmyoMW3SY2w+f6DskkQlJgEKIvGdmAd3mQKeZ8O8OmO+OQ1wcSzstZWjdoQScC2DA5gGER4frO9Icc7C3ZuGwptx79IwRi0OIeVqwtn8UJJIAhRD5Q1GguTcM/gMe3QT/tphf3stHTT9ibvu53I69Tb9N/Vh3YZ3RbydoWMmWuQObcPr6A8b+ruF5fIK+QxLpkAQohMhfjm/C6L+gRAVY3gcOzOWNiq8T2C2Q+qXqM+2faUzZN4WYOOM+Z7Nd7bJ83asBu8/dYeq6k0af1AsiSYBCiPxnXw1GBkGtzrD9E/jzbcpalGCe+zzGNR7H1ktb8dzoyam7p/QdaY70b1aZd9tVZ1VIBHN2XtB3OCINSYBCCP0oUgy8lsGbn8DxVbCoE6YPb+DTyIdFHRYRlxDHoK2DWHJqCQmq8U4hvu9ekz7OlZgVfJ7VIRH6DkekIAlQCKE/Jibw5iTotxLuXgD/N+HKfpzLOhPYLZA2ldrwfej3vLPzHe4+vqvvaF+Joih816cBr9coxZQ/T/DXuYJXIcNYSQIUQuhf7c4weicUKQFLusHhedhYlGDWm7OY1mIaoTdD8djgYbSVJcxNTfh1kAu1yxVn7AoNxyO1+g5JIAlQCGEoStfSVZRwag9bPoL141CeP8Wrlhcru6zEtogtPsE+/C/0f8TFG98xasWKmLFoWFPsrC0YsTiEq/di9R1SoScJUAhhOKxsof8qaDMJji6HRR0hOpIadjVY1XUVfWv1ZfGpxQzaOogrD67oO9qXVqaEJUtGNON5gsqQhYekmK6eSQIUQhgWExNo+wn0+133XNCvDVzag6WZJVNbTGV229lce3QNz42eRrlnsHqZYiwY6sqN6CeMWBJK7DPZKK8vkgCFEIapdhfdlKi1PSztCft/BlWlfeX2BHYLpEGpBkz7Zxof7/mYB88e6Dval+JSxZ45/ZtwIlLL+N+PyEZ5PZEEKIQwXKVr6pJg7c6w41MIHAFPH1GuaDn83f2Z4DyBnVd24rHBg7BbYfqO9qV0qFeOL3vUZ+fZ23z6p2yU1wdJgEIIw1akuG6/oNvncHodzHeDuxcwNTFlVINRLOu8DDMTM0ZsH8FPR34yqjqDg1pUYXy76gSERjAr6Ly+wyl0JAEKIQyfokDr92HQHxBzW7df8MwmAOqXqs+abmvo5tgN/+P+DNs6jIgHxrPh/AP3mni5VmLOrgssP2h8C3uMmSRAIYTxcGoL3n9DqRoQMBCCPoP45xQ1L8r01tOZ2WYmlx5cwmOjh9EskFEUhW96NaBd7TJMW3+SbSelonx+0UsCVBTFVlGUiYqieCT+65zNfm6KonjkdXxCCANm6wAjtoHLcPhnNizrCY90p6t0rNqRtd3WUrdkXab9M42P/v6I6KfReg03O8xMTZg7wJnGDra8u+ooh8Lv6TukQkFfI8A1QKCqqoGqqs4AfBVFsc2sQ+L7foB93ocnhDBoZkWg22zo+StEhoDfG3D1IADli5Vn/lvzmeA8gV1Xd9F7Q28O3Tik33izwcrClIVDm+JgZ8WopaGcvWlcK1uNUb4nwMRE5qiqasrKl+GAWxZdvYDgvIpLCGGEGg+AUcFgZgmLu8CBuaCqyQtklndejrWZNaN2jOL7kO95Fv9M3xFnyq6oBUtGNMPawpShCw8TeV9Oi8lL+hgBugLaNK9pAfeMOiiK4oYkPyFEeso1AJ+/oWZHXWml1UPgiW70VK9UPQK6BuBV04slp5fQf3N//r3/r54DzlwlO2uWjmjO42fxDFl4mKgYw07axkwfCdAWiErz2j0ymNpMHDHaphkxCiHEfyxtoO9ycP8Kzm4G/zZw8wQA1ubWTGs5jZ/b/czdx3fpt6kfy04vM+gSS7XKFWf+0KZcu/+Y4YtDiHkqp8XkBX09A3yZ53huqqoGZtXozp07uLq6Jn/5+/vnIDwhhNFRFHjtXRi2CZ7F6vYLapYlv93GoQ1ru6+lZYWWzAiZgU+QDzdjbuox4Mw1q2bPzwOcORGp5e3lYTx7brgJ25D4+/unygVAqYzaKvm9TDhxOtNPVVWnFK/5ohvl+aRp6wxok0Z/iqL4AWGqqr6Q3VxdXdXQ0NC8DV4IYRwe3Ya1I+HSHmg0ALp8DxZFAVBVlcB/A5kZMhNzE3OmtZhGx2od9RxwxgJCrjJp7Qm6N6rA7L6NMTFR9B2SUVEUJUxVVdf03tPHCDCUF0eAtkBQOm3tATdFUbwVRfFGt1DGPfG/hRAifcXKwOB10GYyHFsJ89rDnXOAbt+dZ01P1nRbQ9USVfl4z8dM3jvZYM8T7du0Mh93qMWGY9f5ctNpo9jbaCzyPQGqqqoFQhVFcUzxsiuJi1wURXFMek9V1WBVVf2TvgANEJTeCFAIIVIxMYW2U2DwHxBzR3d6zLFVyW9XKVGFJZ2W8E7jd9h2aRu91/fm4I2D+os3E++86cSI16qxeP9l5v51Qd/hFBj6egboCXgkboT3BkYnJkYAH2BS2g4pRoCeshleCJFtTu3g7X1QoQn86QPrx+qeEQJmJmaMaTSG5Z2XY2Vmxegdo/E97MuT50/0HHRqiqIwtUsdejepyPc7zrPikByZlhvy/RlgXpFngEKITMU/h93fwt7/Qena4LkYytROfvvx88fMCpvFyrMrcbRx5JvXv6FeyXr6izcdcfEJ+CwL469zt/mpfxO6Nqyg75AMnqE9AxRCiPxnagbtp/03JTqvLRxZDomDACszKz5p/gl+7n48invEoM2D+PXYrwZVXcLc1IRfBjrTtIo97wccZc/5O/oOyahJAhRCFC5O7WDMP1DRRTcd+oc3PH2Y/HarCq34o/sfdKjWgV+O/sKQLUMIjzacbciW5qbMG+pK9TLF8VkWRtiV+/oOyWhJAhRCFD7Fy8GQ9fDmJ3AyEPzawI1jyW/bFLHhu9e/439t/kfko0i8Nnqx/PRyg9k8b2NlztIRzShboggjFofIuaGvSBKgEKJwMjGFNyfB0E0Q91i3cf6QX/KUKMBbVd/izx5/0qJ8C3xDfBm1YxTXHl3TY9D/KV28CMtGNsfS3ITBCw5z5V6MvkMyOpIAhRCFW9XXdKtEndrB1omwagDE/ndaYymrUvzU7ie+bPUlp++dpvf63qw9v9Yg9uM52FuzfGRznscnMGjBIW49MKzVq4ZOEqAQQhQtCf1XQcfv4EIw/PoaXN6X/LaiKPSq0Yu13ddSv1R9Pj/wOe/sfIfbsbf1GLROjbLFWTy8GVGPnjFo/iHuy+HZ2SYJUAghQHeWaIsxMDIIzK1gcVfY9bVu+0SiisUqMu+teUxuNpnQm6H0XN+TjRc36n002MjBlvlDm3IlKpahiw7z8InhrFw1ZJIAhRAipQqNwWePrtbgnhmwuDPc/2/juYliwsA6AwnsHoijjSOf7PuE9/56j7uP7+ovZqClU0l+HejM6esPGLkklMfP4vUajzGQBCiEEGkVKQY9f4He8+H2GfjtdTi5NlWTKiWqsKTjEj50+ZB91/bRa30vtl3epqeAddrXKcv/vBoRcjmKMSukgkRWJAEKIURGGnrC23uhdC0IHAF/jkm1Z9DUxJRh9YexuttqKhWrxMd/f8yHuz8k6knakqf5p0fjinzTqwG7z93hvYAjPI+XJJgRSYBCCJEZu6owfCu8MRGOr9KNBiPDUjVxsnViWedlTHCewK6IXfRa34ugK+kVuMkf/ZtVZmqXOmw5cZNJa0+QkKD/FauGSBKgEEJkxdQM2n0KwzZDwnNY4A57ZkLCf8/ZzEzMGNVgFKu7rqZc0XJ8sPsDPvr7I72NBke97sgH7jVZq4nk/zac1PtCHUMkCVAIIbKrSivdnsF6PWHXdN1KUe3VVE1q2NVgeefljG8ynp1Xd9JrfS92XN6hl3DHt6uOTxtHlh+8yjdbzkgSTEMSoBBCvAwrW+izAHr5w80Tuj2Dx1enamJuYo53Q28CugZQrmg5Pvz7Qz7c/SH3Ht/L11AVRWFyx9oMbVmFeXsvMSv433y9vqGTBCiEEC9LUaBRXxizD8rUhT9G6xbJPE59MHVNu5os77ycd5u8y18Rf+lWil7alq8jMUVR+KxbPbxcKzFn57/8uvtivl3b0EkCFEKIV2VXVfdcsO1UOL0efm0Nl/akamJuYs7ohqNZ3XU1FYtV5OM9H/P+7vfzdd+giYnCt70b0qNxBXy3nWXBvkv5dm1DJglQCCFywtQM2nwMI3eAuSUs6Q7bP4XnT1M1q25XnWWdl/G+y/vsjdxLj3U92HBxQ76NBk1NFP7n2YiO9crx1abTLDsoVeUlAQohRG6o6KI7QcZ1BBz4Gfzbws2TqZqYmZgxov6I5FNkPt33KWN3juVmzM18CdHM1IQ5/ZvQvnYZpq07SUDI1aw7FWCSAIUQIrdYFIWuP8CANf9Vnf/nx1TbJQCq2VRjccfFujNFb4XSa30v1pxfky+jQQszE+YOdOaNmqWZ/McJ/tBE5vk1DZUkQCGEyG0134J3DkLNDhD0f7rtEvcvp2piamLKwDoDWdt9LfVK1uPLA18yascoIh5E5Hl4luam+A92oZVTST5ac4z1Rw2jxmF+kwQohBB5oWhJ8FoGPX+DWyd12yU0S1MV3AVwKO7AvLfm8VnLz3T1Bjf0ZsmpJcQn5O1h1pbmpswf0pRm1ex5P+Aom45fz9PrGSJJgEIIkVcUBRr3hzH/QIUmsGE8rOwHD2+laabgUdODdT3W0bx8c74P/Z7BWwfz7/283bdnZWHKgqFNcalix4RVR9ly4kaeXs/QSAIUQoi8ZlsZhmyADt9C+G74pQWcWvdCs7JFy/JTu5/wfd2XyIeReG3yYu7RuTyLz7sit0WLmLFoeDOaONjy7sojbDtZeJKgJEAhhMgPJibQ8h3dSlG7qrBmKASOhNjUZ4UqikJnx86s77meDlU78Nux3/Da6MXR20fzLLRiRcxYNLwpDSvZMO73wpMEJQEKIUR+Kl1LV3W+7adweh380hLOb3+hmZ2lHd+9/h2/tP+FmOcxDNk6hG8PfUtMXEyehFXc0pwlI5rRIDkJ5s/WDH2SBCiEEPnN1AzaTITRu8DaHn73gvVj4cmDF5q+Xul11vVYR//a/Vl5diU91/dkT+SedD4054pbmrM0OQlqCvxIUBKgEELoS/lG4L0bWr8PR3/XjQYv/vVCs6LmRZnSfApLOy2lqFlRxu4cy8Q9E/PkcO2kJJg0HVqQF8ZIAhRCCH0yKwJun8OIHWBuBct6wqYP4OmjF5o2LtOY1d1W806jdwi6EkSP9T1Yd2Fdrm+gT5oObeRgy/iVRwrsFglJgEIIYQgcmsLbe6HlOAhdCL+2gkt7X2hmYWrBmMZjCOwWSLUS1Zj2zzRGB43O9Q30SUnQubItE1YdLZCb5SUBCiGEoTC3gg5fw/CtYGIKS7rC5o/SHQ062TqxpNMSprWYxqm7p+i1oRfzT8wnLiEu18IpVsSMxcOb4VrFjvcDjha4Y9MkAQohhKGp0hLe/geavw0h83Sjwcv7XmhmopjgVcuLdT3W0bpia37U/Ei/Tf04fud4roVSNHGLRAvHkny45hirQ/L+qLb8IglQCCEMkYU1dPKFYVtAMYHFXTIcDZYtWpbZbWczu+1stE+1DNoyiG8OfcOjZy+2fRXWFmYsHNaU12uUZuLa4ywvIKWUJAEKIYQhq/oajNkPzcdAyHzdaDD873Sbtq/cnvU91tOvdj9WnV1Fj/U92Hl1Z66EkXSAtludMkxdd5L5e8Nz5XP1SRKgEEIYOgtr6PRd4rNBM1jaHTa9n+6+wWIWxfik+Scs77wc2yK2vPfXe0zYNSFXag5ampvyy0AXOjcox/TNZ5j714Ucf6Y+SQIUQghjUaUlvL0vcaXoIt1o8EJwuk0blm7Iqq6reN/lffZf30+PdT1Yfnp5jqtMWJiZMKdfE3o2rsDM7ef4fvu5fKtqn9skAQohhDGxsNatFB0ZpFs1uryP7hSZx9oXmpqbmDOi/gj+7PEnTco2wTfElwFbBnDq3qkchWBmasL/vBrTr6kDP/91ga82nTHKJCgJUAghjJFDU/DZm3iKzEpdhYmzW9JtWql4JX5t/ysz35jJrZhbDNg8AN/Dvjk6V9TUROGbXg0Y1qoqC/+5xCd/niA+wbiSoCRAIYQwVuaWulNkRu8E65Kwqr+uwkTMi0ekKYpCx2od2dBrA541PVlxZgXd13Vn55Wdrzx6MzFR+KxbXca1rc7KwxG8H3CUuPiEHP5Q+UcSoBBCGLsKTWD0X/DmJ3B6PcxtCicCX6g+D1DCogRTW0z9b5HM7vcYv2s81x+92nFniqLwUYdaTOxYiw3HrjNmeRhP4vK2mn1uUYxx3jY9rq6uamhoqL7DEEII/bp1WvdM8LoGanWGLv+DEhXSbRqXEMeK0yv45dgvALzd6G0G1x2MuYn5K1162YHLTFt/ipaOJZk31JViRcxe+cfILYqihKmq6preezICFEKIgqRsXd0CGfev4OIumNsCwpakOxo0NzFnWP1hrO+xnublmzMrbFaOiu8OblmVWX0bcfhyFAPnHeR+TN5Vss8NMgIUQoiC6t5F2PAuXNkHVV+H7nPA3jHD5ruu7uLbw99yM+YmfWr04X2X97EpYvPSlw06fYuxv2uoYm/NspHNKWdjmZOfIkdkBCiEEIVRSScYuhG6zoYbx+CXVrD/J8hgL2C7yu1Y32M9Q+sOZd2FdXT7sxvrL6x/6UUy7nXLsmR4M65rH+Px234u3c2bKvY5JSNAIYQoDKKvweYP4fxWqOAMPX6GsvUybH4u6hxfHfyKY3eO4VrWlaktpuJk6/RSlzweqWXYohBMFFg8vBn1K778aDKnZAQohBCFnU1F6L8SPBaC9ir4vQG7psPzp+k2r2Vfi6WdlvJZy884f/88Hhs8mBU2i9i42GxfsmElW1b7tMTC1IR+/gc5cDH3K9jnhCRAIYQoLBQF6veBcSFQ3wP2zITfWsPVg+k2N1FM8KjpwcZeG+ns2JmFJxfSa30v/o5I/zDu9FQvU4y177SinI0lQxcdZtvJnJ9JmlskAQohRGFjbQ+9/WDQWoh7Ags76KZH0zlcG8De0p6vW3/Nog6LsDKzYtyucUzYNYEbj25k63LlbaxY49OSuuVL8M6KMH4/dDU3f5pXppdngIqi2ALeQDjgCASrqqrJoK0zYA/YJrZFVdUZadvJM0AhhHgFTx/ppkIP/abbL9jlB6jVMcPmcfFxLD29FL/jfsDL7R2MffacsSs0/HXuDu+71eTd9tVRFCXXfpT0ZPYMUF8JMAjwUVU1PMX3nqqqatNpex9on5QgFUVRAZe0CVMSoBBC5EBkKGwYD7dPQ73eumK8xcpk2Pz6o+t8e/hbdkfsprptdaa2mIpLWZcsLxMXn8DktSdYq4lkQPPKfNm9HmameTcZaVCLYBJHf45JyS9ROOCWQZeUyc828TVtXsUnhBCFUiVX8P4b2k6Fs5vg56ZwZHm6G+gBKhSrwE/tfmJO2znExMUwbNswpu6bStSTqEwvY25qwveeDRnzphO/H7rKmBUavR2dpo9ngK68mMC0gHt6jdOM9LyAwDTJUwghRG4ws4A2H8Pb/0CZuroj1Zb2gKiMf+W2rdyWdT3WMbL+SDaHb6bbn91Yc34NCWrGh2IrisKkjrX5ons9Dly8p7d9gvpIgLZA2j8R7qF7zpcuRVEcFUXxBtxVVfXMw9iEEEKUrgnDNuueB14/Ar+0hH2zIf55us2tza15z+U9ArsHUtOuJl8e+JLBWwZz5t6ZTC8ztFVV/v74TeqUL5EHP0TW8v0ZoKIoHsAUVVVdUrw2EWiaVXLLLAlWqVJFLV26dPL33t7eeHt7517gQghRGD24Dls+1k2LlmsA3X/SVZ/IgKqqbArfxPeh36N9qmVA7QGMbTyWYhbF8iVcf39//P39k78PCwu7oqpq1fTa6iMBugF+qqo6pXjNF7BVVdUnG/3vA9+mXQkqi2CEECIPnd6gS4Qxt6HFO9D2E7AommHz6KfRzNHMYc35NZSyKsXEphPpULVDnq/6TMugFsEAobw43WkLBKVtqCiKs6IoF9O8HA683Hk8QgghcqZudxh7CJyHwIGfdRXoL+zMsLlNERumtZzGis4rKGVVio/3fIx3kDdXHlzJx6Azl+8JMHGrQ6iiKCmPJHcFgiH5eV/Se9qk11NwJJ1kKYQQIo9Z2UK3H2HYFjAtAst7w9rREHM3wy4NSjdgZZeVTGk2hZN3T9JrfS9+PvIzT54/yb+4M2AIG+HtgdAUWx1STYcmTpkmJUQXIExVVf+0nylToEIIkY/insDe/8G+WVCkOHT4Bhr10x23loG7j+8yM2QmWy5twaG4A580/4TWFVvnaZgGtxE+L0gCFEIIPbh9RldzMPIwOLaFrrPAvlqmXQ7dOMT0g9O5/OAy7lXcmdh0IuWKlsuT8AztGaAQQoiCokwdGLEdOn+vO03ml5bwz48ZbpkAaF6+OWu7r2V8k/HsidxDj3U9WHJqCc8TMu6TF2QEKIQQIndEX4MtH8G5LVCuoa4CfSZbJgAiHkbw7aFv2XttLzXtajKtxTQal2mcayHJCFAIIUTes6kI/X4Hr6Xw6BbMawfbP4VnGZ/04lDcgbnt5zLrzVlEP41m8NbBfL7/c7RPtHkerowAhRBC5L7HWgj+DMIWg21l3bPB6hkd+awTGxfLL0d/YfmZ5ZS0LMmm3puwMrPKURiyCEYIIYR+XP4HNk6Ae//CgDVQ860su5y/f56Td0/Su0bvHF8+swRoluNPF0IIITJS9TV4ex8cWQbV22erS027mtS0q5nHgUkCFEIIkdfMLaHZaH1H8QJZBCOEEKJQkgQohBCiUJIEKIQQolCSBCiEEKJQkgQohBCiUJIEKIQQolAqUAnQ3/+FKklC5IjcUyK3yT1lOCQBCpEJuadEbpN7ynAUqARojDZu3GhU18nJ57xs3+y2z6pdTt83Jvn5sxjbPZVb91N22sg9pd9rZfdzJAHqmSTAnLeXBPifgvzLKqd9JQG+moJ8TxWYw7AVRbkDxAB39R3LS7IBoo3oOjn5nJftm932WbXLyfulMK57Kr/up9y8Vn7dU7l1P2WnjdxT+r1Wys+poqpq6fQaFZgEKIQQQrwMmQIVQghRKEkCFEIIUSgZTQJUFMVWUZSJiqJ4JP7rnEV7R0VR1iiKknkJYiHSkHtH5FRG99DL/h4TecuY6gGuAXxUVQ0HUBQlSFEUT1VVtWkbprjpHPMxPlEAyL0jciqLeyjbv8dE3jOKBKgoii3gmHTTJAoH3IDAtO1VVQ1O7BeVLwGKAkPuHZFTGd1DL/t7TOQ9Y5kCdQW0aV7TAu75HokQQrwa+T1mYIwlAdoCaf8ivwfY538oQgjxSmyR32MGxVgSIMhNIoQwfvJ7zIAYSwLUovvrKaWSvPjXlBBCGCot8nvMoBhLAgzlxb+cbIGg/A9FCCFeifweMzBGkQATlwiHKoqSclmxK5C02soxzXtCCGFQsvo9JvKfUWyDSOQJeCuKEo7ur6jRKfbO+KD7S8oHIHFzqRu6m2uSoiiOqqpKES6RJbl3RE5lcQ9l9ntM5DM5DFsIIUShZBRToEIIIURukwQohBCiUJIEKIQQolCSBCiEEKJQkgQohBCiUJIEKIQQolCSBChEAZRYbHWNoihq4r9yUIQQacg+QCEKKEVRPIA1qqoq+o5FCEMkI0AhCq6myDFbQmRIEqAQBZcbctCyEBmSBChEAaQoii3gjIwAhciQJEAhCiZXAFVVNfoORAhDJQlQiILJHRn9CZEpSYBCFEzy/E+ILMg2CCEKIEVRVMBFpkCFyJiMAIUwcoqieCuKcjHF9x6AVpKfEJmTBCiE8XMBJqX4fkqa74UQ6ZApUCGMXOIxZ86J3zYFglRVlQUwQmRBEqAQQohCSaZAhRBCFEqSAIUQQhRKkgCFEEIUSpIAhRBCFEqSAIUQQhRKkgCFEEIUSpIAhRBCFEqSAIUQQhRKkgCFEEIUSv8PJK/nEBA+2zoAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 504x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"figsize(7,5)\n",
"for pot,label in zip([iso,mwp14_v,mcm_v],\n",
" ['Isothermal','MWPotential2014','McMillan17']):\n",
" ro, vo= get_physical(pot)['ro'], get_physical(pot)['vo']\n",
" # The following Emax is in (km/s)^2\n",
" Emax= evaluatelinearPotentials(pot,4.*u.kpc)\n",
" Es= numpy.linspace(0.*u.km**2./u.s**2.,Emax,101)\n",
" # Setup Orbits to easily pass the required (z,vz) to the actionAngleVertical calculation\n",
" orbs= Orbit([[0.,0.,0.,0.,numpy.sqrt(2.*(E-evaluatelinearPotentials(pot,x=0.)))\\\n",
" .to_value(u.km/u.s)/220.]\n",
" for E in Es])\n",
" # The following action J and frequency OJ are in internal units of ro*vo \n",
" # [different for McMillan17!] \n",
" J= numpy.array([actionAngleVertical.actionAngleVertical(orb,pot=pot).Jz() for orb in orbs])\n",
" OJ= numpy.array([2.*numpy.pi/actionAngleVertical.actionAngleVertical(orb,pot=pot).Tz() for orb in orbs])\n",
" # We now convert the action and frequency to Scott's units of 30 km/s x 0.3 kpc\n",
" # and the frequency to the implied scale in this unit system\n",
" semilogx(J*ro*vo/(scott_ro*scott_vo),\n",
" OJ*freq_in_Gyr(vo,ro)/freq_in_Gyr(scott_vo,scott_ro),\n",
" label=rf'$\\mathrm{{{label}}}$')\n",
"xlabel(r'$J$')\n",
"ylabel(r'$\\Omega(J)$')\n",
"gca().xaxis.set_major_formatter(FuncFormatter(\n",
" lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\\\n",
" numpy.maximum(-numpy.log10(y),0)))).format(y)))\n",
"legend(frameon=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $\\mathrm{d}\\ln \\Omega(J)/ \\mathrm{d} \\ln J$ vs. $J$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we compute the logarithmic derivative of the frequency with respect to action:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f at J=1.79 for Isothermal is 0.1249207757247324\n",
"f at J=1.79 for MWPotential2014 is 0.1965060677124203\n",
"f at J=1.79 for McMillan17 is 0.17178291476723556\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1691e2340>"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAFCCAYAAACEgRbZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABPdklEQVR4nO3deVhU1/3H8fdBXECBAcQN1wH3HcZEY9IsYvbFRNEkzdKkEdI2abqk0qRbmraxmqZN01+aSpqtaZKqaEw0UQPGbGqiMO47DCiKCwgDIso25/fHHQiOIIswC3xfz8MDzNx756AjH8+555yv0lojhBBCiG/5eboBQgghhLeRcBRCCCFcSDgKIYQQLiQchRBCCBcSjkIIIYQLCUchhBDChb+nG+AuPXv21IMHD27z1ykuLiYkJMTrrtmSazT3nKYe35TjGjsmPz+fiIiIJrfN27XF+8ZTr+up92tzz5P3a8u1l/drRkZGgda6/r8YrXWH+IiNjdXuMHfuXK+8Zkuu0dxzmnp8U45r7Bh3/X26S1u8bzz1up56vzb3PHm/tlx7eb8C6bqBzJBh1VZ22223eeU1W3KN5p7T1OObclxb/Dl6M0/9vO3p/drc8+T92nLt6f3aEKU7yA45FotFp6ene7oZopVYLBbk71P4Cnm/eielVIbW2lLfc9JzFD4pISHB000Qosnk/ep7pOcohBCiQ5KeoxBCCNEMEo5CCCGECwlHIYQQwoWEoxBCCOFCwtFLWa1W4uPjiY2NJS0trdnn2+32875PS0sjNDS0lVrXfGlpaURFRXns9YUQojkkHL1UTEwM06dPx2KxEBcX1+zzlyxZct73cXFxWCz1Tspyi7i4OMxms8deXwjhG7TW5J8ux3q4iGPFZz3Wjg6zt2pHs2jRIllbJYTwSlprTp4ux5Z/huyCM+ScOkNOwRkOnSojt6iMsopqAH59y0geucoz/6mWcPQRNpsNm81GWFgYixcvZsGCBQCkpKRgMpmw2+2YTCbi4uJIS0vDbreTkpKC2WwmJiam9jpWqxUwwnPRokW1jycnJ2M2m7FarSQkJJCWlkZSUlLtcd/97nd56KGHWLduHYWFhSxdupSkpCTsdjuLFi0iKSmptmdYMwycmprKU089hclkctOfkhDCmzgcmqP2sxw4cZr9J06TebKUrJOlZOWfobS8qva4Lv5+DAwLZHB4IFOjezIwLIABYYGM7uf+zc1rSDjW8fuVu9mTV9KmrzGqXzC/u210s89btGhRbdDYbDbACDqbzca8efMASExMrB2GNZlMzJo1q95rxcTEYDKZSEtLIy4uDqvVSlZWFgkJCVgslvNCMSwsjAULFmA2m3n55ZcBY4g0NTWVlJQU5s2bR3x8PIsWLaoN7KVLl9YG7/z582sfF0K0X+cqq9l7rIQ9x0rYk1fCvuOn2X/89Hkh2Du4K0N7BTEzJpKoXj0Y0rM7Q3p2p19IAH5+yoOtv5CEo4+YM2cOsbGxxMXFkZiYCBiBGR8fX3tMVFQUaWlpDYYiUNu7Cw8Pr31s8eLFhIeH1/Yqa8K3sLDwvF6n6/k1X4eFhZ13TN0eqevEICGE76uqdrD/xGm2Hraz44idHUeKOXiylGqHseNaUDd/RvYJ5q6YSEb0CWZ4nx5E9woiJKCzh1vedBKOdbSkR+cuZrOZjIwM0tLSiI+PJysrq97jCgsLz/veZrOdNxGmviFOu93O9OnTa4MwNTUVuDD0XM9vaLg0KSmJ6dOn13u+EML3FJ+txHqoiPRDhWQcKmJ7bjFnK437gqGBnRnb38S0kb0YGxnC6H4h9A8NQCnv6gk2l4Sjj0hOTmbevHnMmjWLwsJC7HY78fHxWK3W2tmsW7ZsqR3CrAkm13CsT3x8PKmpqbXXsVqtF/QYm9PO8PDw2uHaS72eEML9Ss5V8o2tkE1Zp/jadoq9x0vQGvz9FKP6BTNn0gAmDjQxcUAoA8J8PwjrI+HopaxWK6mpqdhsttoJLjUTbMxmc+3km5rn7XY7iYmJtUGYmJhYO8mm5no2m42FCxcya9YsUlNTycrKqr1HWfd1wsLCSEtLIz09neTkZBISEi56/qJFi0hPT8dmsxEXF8eCBQtIS0urDejCwsLa82uuJ4TwHlXVDrYfsfP5gQK+PJjP9lw7Dg1d/f2IGRjKE9OGctngMCYMNBHYpWPEhlTlEEKIDsheVsFn+/NZt+8kXxzIp/hsJX4KxvU3cWV0T6ZG92TiQBPdOnfydFPbzMWqcnSM/wIIIYTgWPFZ1u46ztrdJ9icU0i1Q9OzRxemj+rNNcMjuCo6gpBA35k005YkHIUQoh07WXKOVTuOsWpHHtbDdgCG9urBo1ebiRvZm/H9TV63jMIbSDgKIUQ7U1pexeqdx1ix7Sgbs06hNYzsG8yT1w/jprF9iYro4ekmej0JRyGEaAe01nyTXciSLbms3nWcs5XVDAoP5PHrhnL7+H5E95JAbA4JRyGE8GGFZypIycjlvc25ZBecIairPzMmRjIrNpKYgaHtcpmFO0g4CiGED9p5pJg3NmazascxKqocWAaF8ti10dw8ti8BXdrvDFN3kXAUQggfUe3QpO45zr+/zCb9UBGBXTox29Kf+ycPZnifIE83r12Reo5eymazkZiYSGxs7AXP2e12QkNDWbhwIWBsDhAaGkpSUlLtuTWFkutuUh4bG0tKSsp5z6ekpJCcnExSUlLtjjYNudR9UltyvhRJFsLY1Pu/Xx/iuhc+49H/Wjlx+hy/uXUUXz89jT/OGCvB2Aak5+ilzGYz8fHxFBYWXrAFXE3pqppqHDU71tSEiNls5qmnnmLu3Lm158XExLBgwYLaLeKmT59ORkbGeZuUh4aGkp2d3eCeqUuWLLmk3W3qOz8qKqrBfWLhwiLJKSkpFBYWkpGRQXx8fO3PU7d0V90yXXa7neTkZEwmU71trzmvJQWlhWhrZ8qreOebQ7z6ZTb5p8sZP8BE0o0juGF0HzrJ8os25ZGeo1LKpJSap5Sa5fzc4MabSqkYpVRcnWPnteQ6vmrOnDnnVbloSHx8PEuXLq39viZU6/bW6m4EXt+m4GFhYbU9zfo0pR0XU9/5GRkZTT7farViNptJSEg4ryKJ3W6v3Rt21qxZtT1ogIvtilRTi1IIb1NWUcW/Ps/iqoXree7jfQzr3YN3H7mcFT+8gpvH9pVgdANPDasuBVK01ila64XAAqWUqYFj1wGFLsfGtOA6PmnWrFm1e55CwxuJx8XFXRAEs2fPZsmSJUDTNv8uLCzEbDaTkpJCWlpa7WfgvALKdYdfk5OTSUtLY+HChdjtdtLS0ggNDcVqtZKWllZbXqu+82uGemukpaXVFlmubwi2sLDwvDALCwvDarWyZMmS84ZeTSZT7WvU1LasT3p6OtOnT7/on4kQ7lReVc3rX2XznYXr+fPqfYyJDGHZD67gnUcmc0V0T5l56kZuH1Z1hpdZa123i2ID4oCUek6ZprW21jkXwN6C6zRu9S/h+M4WndpkfcbCTX9u1ilms7k23GqGDetjsVhqCxiHhYURHx/PggULSEhIuKCUFXDepuVbtmxh3bp12Gy2JhdQbqhIssVibFVYUxS5pk2u58fExJz3szRWJDkuLu684c+aepNpaWnnBWBYWFi9P29dVqsVi8XS6H1WIdzB4dB8sP0of1l7gKP2s0wxh7Po/mHEDpKyb57iiXuOFsDu8pgdmE49oVYTjE6zMXqKNqVUXHOu48sSExNZvHhxoz2/mqHVsLCw2mPrFkN2ZTaba8OmJrQSExObXEC5oSLJNdeG84sqN6Y5RZKTkpLOG0Z21ZTJPw31KIVwp01Zp/jTx3vYdbSEMZHB/HnmWK4aGuHpZnV4nghHE+D63/pTQINFB5VSZowe4XStdc1v7mZfp1HN7NG5S1xcHImJiSQmJl40IGfPnk1SUtJ54WaxWFi4cOElTaRpqIByQ0WS4eLB09DQcFOLJKekpDBnzpza162ZiFO3vRerYZmcnFx7f3XLli2cOnWqthSYEO6SW1jGcx/vZfWu40SaAvj73RO4bVw/2efUS3jqnmOzxgq01jatdTKQqpSq211o8nXy8/OxWCy1H8nJyc1pgsfFxcWRknLxDrHJZKoNrRrx8fEsXry4yb2kmgLKNbZs2VLbu6xbQLnm2LqB2NgQpev5ddUtknyx69UML8fExNQOAc+ePfu8Ga92u/2i/4lISEhg1qxZzJo1C7PZzPTp0yUYhducq6zmb6kHmPbXz/lsfz5PXj+MdT+/mjsmREowtrHk5OTzcgDo2dCxnug52jF6fXWFc2Ev8AJa62Sl1ALnjFVrc64TERFx0ZmL3sZms7FgwQKsVivz5s07rzdYs5yhvsLBc+bMOS9gZs+efcEQo81mY/HixbX3HOse35wCyvUVSb5YUeSGCjAnJyc3qUiyxWJh2rRp5z1fVFQEGEtTatpRd7ZqWloaqamptfdqXcPXarXWPicBKdraZ/tP8tsPdnO4sIzbxvfjVzePpE9IN083q8NISEg473emUqqgoWPdXuzYOZEmW2sdWuexRUCq1jrF5dgYYKnWOqrOYxlAOpDU1OuAFDsWQnjOydPn+P3KPXy04xhREd35w4wxXBHVYKdFuIlXFTvWWtuVUulKqbozTS0YYVdzfxHnc3YgzeUSZmB+Y9cRQghP01qTknGEP360l7MV1fw0bhiPXmOmq7/sfertPLVDTjyQoJSyYdw3nKu1tjufS8QYLk10zkpdqpSq6QfHAkl1eoYXu44QQnjMseKzPLV8J5/tz2fS4FDm3zVOykb5ELcPq3qKDKsKIdxBa81y61GeWbmbqmpN0o3DeWDKYJls44W8alhVCCHaq8IzFfzq/Z2s3nWcSYND+Uv8eAaFd/d0s0QLSDgKIUQr2JBZwE8Xb6OorIJf3jSCuVeZZQ9UHybhKIQQl6Cy2sELnxxg0RdZmHt25/XvTWJMZIinmyUukYSjEEK00FH7WR5/14r1sJ17LhvAb24dRWAX+bXaHkixYy/VnGLHrXl+3eLCUmhYiIat33eSW176kgMnSvnHPROZf9c4CcZ2RMLRS9UUOzabzRdst+Za7Lg1zq8JwbrFhV0LDbeE3W5n4cKF523XZ7fbz9vFRghfUu3QvPDJfh56cwt9QwJY9fiV3Da+n6ebJVqZhKOXa2qx40s9vzlFh5ujvuUzNdvBRUVFERUV1WgvWAhvYS+r4KE3t/CPTzOZbenP+z+8gsE9ZTZqeyTh6OWaWuw4JSWltpBw3Q27m3K+a9HhhtRXjLih4sY1Gio2XFRURFZWFllZWbz66qsX7QUL4Q32HS/h9v/bwKasAp67cywLZ42nW2fZ6aa9kgHyOhZsXsC+wn1t+hojwkaQdFnzhhQbK3ZstVpJTU1l0aJFtUOWdXuLjZ3vWnS4IfUVI46Li2uwuHFD6lbMSElJueixQniDNbuO87Ml2+jR1Z//JUwhdlBo4ycJnyY9Rx9QU+y4IYsXL2b69OmAUbbKdRi1sfOb6mLFiFtS3Nhut2Oz2aTosPBaWmv+se4gj/43g6G9g1j5+JUSjB2E9BzraG6Pzl2aWuy4rc6vcbFixC0JuPnz5zNnzpwWt0eItlReVc0vl+3k/a1HmTGhH3+eOU6GUTsQ6Tn6iIsVO54zZ06jRYebUiz5YppajLg50tLSpNcovFLhmQru+/c3vL/1KE9eP4y/zZkgwdjBSDh6qZpixzWzOOPj42t7fXWLHYNxDy82NpaUlJTawGnO+XULCjf0dVxcHFlZWedN7qlbjHjhwoXYbDZSU1NZunTpeRN2ah6re24NKTAsvM2hU2eY+cpGth8p5h/3TOSx64ailGwD19FIVQ4hhHDalmvn+29uoVprXnvQQuygC28hiPZDqnIIIUQj1u87yQ/fsdIzqAtvPnQZURFSe7Ejk3AUQnR4y61H+EXKDkb2DeKN711GRFBXTzdJeJiEoxCiQ/v3lzb++NFerogKZ9H9sQR16+zpJgkvIOEohOiQtNb8NfUA//g0k5vH9uFvcybQ1V9mpAqDhKMQosNxODTPrtrDmxtzuHvSAP5051gpTCzOI+EohOhQqh2aXy7bwdKMIzxy5RB+dctIWaohLiDhKIToMKqqHfxsyXY+3J7HE9OG8pM4WcMo6ifhKIToECqqHDzxv62s3nWcpBtH8INrpJC3aJiEoxCi3auocvCjd62k7jnBb28dxcNXDvF0k4SXk3AUQrRrdYPx2TtG88CUwZ5ukvABEo5CiHarosrBD9+xkrZXglE0j2w8LoRolyqrHfz4va0SjKJFJByFEO1OzazUNbuP89tbR0kwimaTcBRCtCsOh2Zeyg5Wbs/jqZtGyOQb0SISjkKIdkNrzW8+2MXyrUf5+fRhJF4tyzVEy0g4CiHaBa01f169j3e+OcyjV0fx2HXRnm6S8GESjkKIduH/Ps1k0Rc2HpgyiKQbh8vON+KSeGQph1LKBCQANsAMpGmtrQ0cGwPEOb+dBCzSWqc5n5sHhAOLgTAgXmud2LatF0J4m7c25vBC6gHuionkmdtGSzCKS+apdY5LgUSttQ1AKZWqlIrXWtvrOTZOa73QeZwJyFZKTasTpgnOjzRgbpu3XAjhVT7YdpTffbib6aN6s3DmOPykuoZoBW4fVnUGnLkmGJ1sfNs7rHtsDPBUzffO8Eyvc6xdax3q/GgoXIUQ7dT6fSf5+ZLtXD4kjH/cMxH/TnKnSLQOT7yTLIDd5TE7MN31QGfvMN7lYbPr+UqpGKWUudVaKITwehmHivjBOxmM6BvEvx+00K2zFCoWrccT4WgCCl0eO4Vxz/ACNfcXAZwBGAYsqfPYLIyeZ4xSakFDL5qfn4/FYqn9SE5ObvlPIITwqMyTp/n+W1voE9yNNx+6jKBunT3dJNESWkN1ldteLjk5+bwcAHo2dKzSWrutYVAbZk9prWPrPDYPmKS1du0lup6bCiRdZPJOFsa9zDTX5ywWi05PT7+0xgshPO5Y8Vlm/nMjFdWa5T+4goHhgZ5ukmiJskJY9VPo3hNuecEjTVBKZWitLfU954meox2j91hXOBf2Js/jDNAFdYPReU+yLiv1DM8KIdqH4rJKHnx9MyXnqnjzoUkSjL7qYBq8cgXs+whC+hs9SC/jidmq6Vw4hGoCUhs6wdnbrF3u4RxeNQHrgFCX62S1XlOFEN6ivKqahLfTyS44w5sPXcaYyBBPN0k017kS+OTXYH0LIkbAvYuh73hPt6pebu851sw4dZlAY8FYioFSylz3OaVUHMas1JpgNAExzu+TXC5vps79SCFE++BwaJ5cuoNvsgv5S/x4pkY3eKtIeKuDqfDPKbD1bZj6BCR87rXBCJ5b5xgPJCilbBi9yLl1lmEkYvQAE50hmQq4LuqtuV+Z7hxutQNRGJsA2BFCtCsL1uxj5fY8km4cwR0TIj3dHNEcZwpgzVOwcwn0HA7fT4X+9d7m8yoeCUdngC1s4LmkOl/bgAZX9Dp7j/VOzhFCtA//2ZTDoi9s3D95EI9eLSu2fIbDYfQSU38LFWfg6iS46ufg39XTLWuSBsNRKTUEWAC05E6pcp6XpLXOaVnThBAd3bq9J3jmw93EjezFM7fLtnA+I28bfPwkHNkCg6bCrX+DiOGeblWzNBiOWutsYLYb2yKEELV2HS3m8fe2MqpfMC/dM5FOsi2c9yvNh0//ANb/GEs0ZrwC4+8BH/xPjafuOQohRIOOFZ/l+29twRTQmdcfnERgF/lV5dUqz8HmRfDFX6CyDCb/wBhGDTC1+ktVOapYZVvF7VG346fabk7pRd9xSqnBMiwqhHCn0vIqHnpjC2fKq0n5wRR6BXfzdJNEQxzVsHMpfPpHKM6FoTfA9X+EiGFt8nKHSw7z9FdPsz1/O6auJq4ZcE2bvA403nNcAMxps1cXQog6qh2aH7+3lYMnS3nje5MY0SfY000S9dEa9n9shOLJPcaSjDteBvPVbfJyDu1g8f7F/C3jb/j7+bPgqgVtGozQeDjG1yzAx1hSkaa13tamLRJCdFh/+mgvn+47yR9mjOE7wyI83RzhSmtjveJnz0HeVgiLglmvw6g7wa9thjhzT+fyzMZn2Hx8M1P7TeX3V/ye3t17t8lr1dVYOKZgFBKeBNwNLFRKaSQshRCt7L9fH+L1Ddl874rB3D95kKebI+pyOGD/R8Y9xWPbwDQQbv8HjL8XOrXN/eBqRzXv7XuPl7a+hJ/y43dTfsfMoTPdNmO5sZ9qi9Z6GbAMQCkVghGUcXwblkXAPK31a23aUiFEu/XVwQJ+9+Furh0ewW9uHeXp5ogaVeXGPcUNL0HBfggd7AzFe6BT21VCySzK5JlNz7A9fztTI6fyu8m/o2+Pvm32evW5aDhqrZ93+b4Yo9dYt4zULOCXSqlsrfWnbdJKIUS7lZVfyg/eySA6oocs2fAWZYWQ8QZ8kwylx6H3GLjr3zD6zjbrKQKcqzpH8o5k3tj9Bj069+C5K5/jVvOtHlnfesk/pdY6BUhRSs0HJByFEE1WXFbJ3LfS6dzJj38/aJG6jJ52bDtsftXoLVadA/O1MONliJrW5msVvzjyBc998xxHS49ye9TtPGl5ktBuoY2f2EZa878AUixRCNFkVdUOfvSuldyiMt55ZDIDwqT8lEdUnIHd70P6G3A0HToHwvi74bJE6N32Q9y5p3N5fsvzrM9dz5CQIbx2/Wtc1veyNn/dxlxyOCqlDmJUx/C+glxCCK/1h1V7+CqzgIUzx3HZENcqdqJNaQ25m2HbO7BrOVScNjYFv2E+TLgHAtq+x1ZWWca/d/6bt3a/RSe/Tvwk5ic8MOoBOrfhvczmaI2e4/PA9cCfW+FaQogO4N1vDvPWpkM8cuUQZk8a4OnmdBynsowh0x2LodBm9BJHzYCY+2HgFLds81btqObDrA95aetLFJwt4BbzLfw05qduWZ7RHK1xzzG5NRoihOgYvrGd4rcf7OLqYRE8dfNITzen/Su0wZ4PjKHTY9sBBYOvhKuehFG3Q9cgtzRDa82XR7/kReuLHCw6yLiIcbx47YuMj/DOmo6yYaEQwm2OFJXxg3esDAwLlJmpbUVrIwT3r4Z9q+DELuPxyFi44Tmjpxji3pqY1hNWXtr6EhknMhgQNIDnr36eGwbd4NVVViQchRBuUVZRxdz/ZFBZ7eDVBy2EBHjHvaV24VwJZH9u7F5zMBVO5wEKBk42AnHkbcbCfTfbnr+dV7a9woa8DfQM6MmvLv8VM4fO9Jr7ihcj4SiEaHNaa55cup39x0t47XuTiIro4ekm+baqcjiaAbbPwfaZUTdRV0OXIIi6FobdCEOvhx7u34JPa03GiQySdySz6dgmQruG8rPYn3H3iLsJ8A9we3taqtXCUSn1pNb6L611PSFE+/Hy+kw+3nmcp28ewbXDe3m6Ob6n/LQRgIe/hkMbja+rzgEK+k2EqU9AdBwMuKxNd665GId28FnuZ7yx6w225W8jrFsYP4v9GXOGzyGws+8t02lWOCqlpmFU6hji+hQQAkg4CiHOs27vCV5IPcCMCf2Ye5XZ083xftWVcHKvsbH30Qzj4+Qe0A5QftBnLMQ+ZEyqGTzVLcsuLqassowPsj7gnb3vcKjkEJE9Inn68qe5M/pOuvn7brmx5vYcpwPTnNvInUcp9YvWaZIQor3IPFnKE//bxuh+wfx55jivnoDhEWeLjCA8sRuO74DjO+HEHqguN57vFmJMpBlxq9Er7D8JunlHGS9bsY2l+5eyInMFpZWljO05lue/8zxxg+Lw9/P9O3bN/QlS6wtGuHAfViFEx1ZyrpKE/6TT1d+PRfdb6Na5k6eb5BkOB5QchVOZxkfBQWMT75P7jH1LawSEGr3Cy+YaQ6X9JkKY2S1rD5vqbNVZ1h1ex7IDy0g/kY6/nz/TB03n3hH3MqHXBE83r1U1Nxwb3AVHKXWdbDwuhACjaPFP/reNw4VlvDt3MpEm35mI0Wxaw5kCKM41Puy5YD8ERYegKMf4qOkJAnTpAT2HGhNneo2EiJHQezQE9/OqIKzh0A62ntzKyqyVrM1ZS2llKQOCBvBEzBPcGX0n4QHhnm5im2huOF6vlEoC7ECh8zNAKDANGNpqLRNC+KwX0w4YRYvvGO27W8NVV0HZKTiTD2dOQulJKD0Bp08YPb6SY8aSiZJj54cfGLNGQwcZITjsBqMHGB4F4dEQ1NcrQ7AurTV7C/eyJmcNa7PXkncmjwD/AKYPms6M6BnE9o7FT7VNcWNv0dxwjAPm820o1jABcqddCMGaXcf4x6eZzLEM4D5vKFpceQ7OFUN5ifH5nB3O2o37fTVflxXC2cJvP58pMJ6rT+dACOpjhFykBUb2hZABEBwJpgHG1wGhXh+Arqod1ews2Mmnhz8l9VAqR0qP4K/8mdxvMo9NfIxpA6f55KzTlmpuOM7VWm+t7wmllK0V2iOE8GEHT5zm50u2M36Aid/fMfrSJuBUnjOWMJSXOD/X/ajzWEWp8flcifPxkvO/r664+Ot07m6EWWAoBIRByFgI7Ande0JgOPToBd0joEdv4+suPXwu+BpyuuI0Xx/7mi+OfMEXR76g8Fwh/n7+XN73ch4Z+wjTBk7D1M3k6WZ6RLPCsaFgdPLsfGIhhEcVn60k4e0MArr4s+i+WLp1UkbvrKzQ+Hy26PweW00v7lyx86Pk2x5e+enGQw3Az9/YG7RLkDGLs2uQEWLhQ53fB9f5HALdTMbnAJPxdYAJ/Lu24Z+Kd6lyVLHn1B425W1iY95Gtudvp1pXE9QliCsjr+Sa/tdwVf+rCOrinv1WvdlFw9FZwLgpFMY9x0mX3CIhhPeqroTTx7+9/1ZqfK1LT3Jgz0H+XJrPuNBKAhbZjeFJ7Wj4Wp27OwPLGVw9ehn35GrCrGtQnXALqvNR83wP8O/WbnpxbaHaUc2+on2kH083Pk6kU1pZikIxImwED495mKmRUxkfMb5dLL9oTY39acQCi7jwHmN95J6jEL7M4TACrzgX7Iedsy+PQPHRbyeenMmnvknr5/yDCa3oQVB4HwL6RBvDkoHhEBhmfA4IPf+jW4jHdnJpz8oqy9h9ajdbT27FetLK9pPbKa0sBWBg0EBuGHwDk/tO5vK+lxPaTQb7LqaxcExqZCi1ltxzFMIHOBxQcgQKDsApGxRlGyWNCrMvXHIAxtBjSH9j8knf8RDUD4L7Qo8+Rk+vR29SDzuY+84OZlv6s2DmOOnJuYnWmkMlh9hZsJMd+TvYnr+dA0UHqNbVAESborl5yM3E9I4htncsfbr38XCLfctFw7Gpweg8NvvSmyOEaBUOB9hzjJ1XTuw2wrDgABRkQtXZb4/rHAihQ5xLDq4H0yCjekPIACMUG9mNJSu/lJ+mbGB8/xCevWOM7IDTRrTWHD9znN2ndhsfBbvZdWoXpytOAxDoH8iYnmN4eMzDTOg1gfER4wnpGuLhVvs2jwwyK6VMQAJgwxiOTdNaWxs4NgZjCQkY9zQXaa3TmnsdIdqtqnI4vgvyrEbtvhO7jS3IKs84D1BG4PUcBoO/YwRhz2HGursevVvc0ystryLx7Qy6+vvxyn2xHXcHnFbm0A6Onj7KnsI97D21l72Fe9l7ai9F5UUA+Ct/okOjuWHwDYwJH8PYiLFEhUTRyU/+/FuTp+7ALgUStdY2AKVUqlIqXmttr+fYOK31QudxJiBbKTXNGYLNuY4Qvs9RbezFeTTD2Jg6z2oEoaPSeD4gFHqPgZj7jV1Xeo+GiBHQpXurNkNrzZNLtpNdcIb/fv9y+rXnHXDaUHl1OZn2TA4UHmBf4T72Fe5jf9F+zjj/Y1MThNcMuIZR4aMYFT6KYaHDfHpDb1/h9nB0Bpy5JtCcbBi9wxSXY2OAp4CFAFpru1IqHYhz3uNs0nWE8FkVZUYQ5n5tlCvK3WwsdQBjUku/iXDFY9/uxRkywC33/P75WRZrdh/n17eMZEpU+9w+rDVprTlRdoIDRQe+/Sg8QE5JTu09wgD/AIaFDuNW862MDBvJyPCRRJui6dKpi4db3zG1Zj3H+Vrrp5pwqIULZ7/aMSp+nBdqWmurUire5Viz8/gmX0cIn1F5Dg5vMgrY5nwFx7aBo8p4rtcoGDsLBkyG/haPbUr9+YF8/vLJfm4f34/vX+lavU6UVJSQZc/iYNFBDhQd4GDRQQ7aD9beHwTo270vw0OHc93A6xgRNoLhYcMZEDSg3W/J5ksaW+eYyUU2G68jHKOeY1PC0YSxL2tdp2hgKUjN/UVne8xAGLAEo4fY5OsI4ZUcDqNUke0zsK03eodV54zF7ZEWuOJxGDjFKFUU6Pk9SnMLy/jxe1sZ3juIP88c26En4JRVlpFlzyLTnml8Ls4ksyiTE2Unao/p0bkH0aZobhp8E0NDhzIsdBjRodEEd/GOslOiYY31HNO01o/WfKOUmggXzmJ1FkF2DaqLaem/8kUY9STtzn+UTb5Ofn4+Foul9vuEhAQSEhJa2AwhLsG5EshMg/2rIWudsbk1GNUZYh8yqjUMmmoscvciZyuqSXw7A601i+6PJbBLx1g0frriNLZiGza7jSx7FlnFWdjsNvLO5NUe08WvC1GmKCb1mUS0KZqhoUOJNkXTt3vfDv0fCG+TnJxMcnJy3Yd6NnRsY0s5HnV5KLS+slRa63VKqeua2D47Ru+xrnAaCVel1DxgQZ3ZqM26TkREBOnp6U1sohCtzJ4LB9bA/o8h+0tjAk1AGAydDlHXwZCrjfWDXkprzdPv72Tv8RJe/94kBoW37gQfT9Nak382n+zi7NogrPk6/2x+7XFd/LowJGQIE3pN4K6Qu4g2RRNliqJ/UH/ZYcYHuHaKlFIFDR3rib/NdC7s8ZmA1IZOUErNos4yDefwarOvI4RbndwHu983AvH4DuOx8GiY/AMYfrNR2d1Hpt+/tTGH97ce5efTh3Ht8F6ebk6LlVeXc6jkEDnFOeSU5JBdnE1OcQ7ZJdm1M0QBunfujjnEzJR+UzCHmDGHmIk2RdOvRz9ZMtFBNDccY5VShVrrbXUfVEoNxthqrtFixzUzTpVSdWeaWoAk57XMzuNqlmfEAfY6wWgCYrTWKRe7jhAeUXQIdi0zPk7sAhQMnAzTnzUCsafvlTzdnF3IHz/aS9zI3vzo2mhPN6dRDu3gxJkT5JQYAZhTnGMEYkkOeaV56DrTKHoH9mZwyGBuj7qdISFDGBIyBHOImYiACBkO7eCU1k2Zb1PnBKWWAEMwlk2AMQHGprWe04xrmPh28X4YkF4n/BYAJq11ojMos+q5RKxzJmuD13FlsVi0DKuKNlF60ugh7kyBI5uNx/pPgjGzYPSdENTbs+27BCdKznHLS18R3M2fFY9NJbib9+yHWlxeTE5Jznk9wUMlhzhccphz1edqjwv0D2RQ8CAGhwxmSPCQ2q8HBw/uUPUJxYWUUhlaa0u9zzU3HJ0XHALEOL+1+sLWcRKOolVVnoU9H8L2dyH7C6P6RK/RMHYmjJkJoYM93cJLVlHl4O7kTew7fpoPfjSVob3dX8aovLqcwyWHa4MvuzibQyWHOFRyCHu5vfY4f+VP/6D+DAweyODgwUYAOj/3CuwlvUBRr4uFY4PDqkqpCa7DpzWcYej1gShEqzuxGzLegh3/M2oPmgbBlT8z1h/2Gunp1rWqZ1ftxnrYzj+/G9Omwai15mTZyW/vATqHQusbBu0V0ItBIYOYPmj6eQEYGRRJZz/v6dUK33exe45RSqmnMdYOLq1vlqoQHULFGdi1HKxvwZEt0KkLjLwNYh6EwVeBX/tbuL0kPZf/fn2YxKvN3Dy2dWbRVjoqyT2de95M0OzibLKLsymrKqs9LsA/gMHBgxkXMY47ou6oHQIdFDxIhkGF2zQYjlrrZcAypVQIMFsp9SjG/b/FDfUohWhX8rZBxpvGvcSK08Zm3Tc8B+Puhu7td8u0nUeK+fWKXUyNDucX1w9v9vmVjkoOlxwm055Zu0DeZrdx6PQhqmp2+wH6dO/DkOAhzIiewZCQIbX3BGUYVHiDRmeraq2LgVeBV51BmeDsUWZhVMjIadsmCuFG1VWw90PY9DIcTTcqzY++0+glDpzc7msVFp6p4NH/ZhDRoysv3T0R/04N94q11hw7c+yCbdJySnJqQ1ChGBA0ALPJzDUDriHKFIU5xMyQkCHSCxRerVlLOZxB+TzUTsqJV0pZgC1AigSl8FnnSmDr2/D1v6D4sLFv6U0LYdxso9JFB1BV7eCxd63kl5aT8ugUwnt0rX3uXNU5suxZtZUjajbPrqkyDxDZI5JoUzRX97+aKFMU0aZohoQMkQoSwie1eBMA56ScukH5qPPzFiBZa13SOk0Uog0VH4Fv/mVMsikvgYFXwI3zYfhNPrNAv7UsXLufjVmn+OOd0VR2zuLtPd/WE8wuzq6tHtG9c3eGhw7nFvMtDAsdZuwXaoqmRxfv2u5OiEvRoqUcF72gsf/qHIy1kEne0puUpRziPHlbYeP/GesTAUbdAVMeg/6xnm2Xm52rOse+wn38b8cGVuz5mtDQk5xxHKudIRoREMHI8JGMCBthfISOIDIoUqpHiHahRUs5Wsq5KfnWRg8Uwt20hoOfwIa/w6EN0CXI2Mrt8kQwDfR069qcQzvIKc5he/52dhbsZFfBLg4WHaRKG/cHA4JMxPabyOieMxgdPpqRYSOJCIzwcKuF8AzZKVe0f1obVTDWPwd5VqMg8PV/gpgHoFv7LR1UWlHKzoKdbM/fzrb8bezI31FbUzCocxCje47mnuEP8P7XflSVRfLRj26lV7DcHxQCWhCOzn1UzXxbEcOOsX1cTms1SohWobVREmr9fGPmachAuO0lmHAvdGp/C8YLzhaQcSID6wkr1pNWDhQdwKEdKBRRpiiuH3Q94yPGMz5iPINDBoNWfP+tLRScLOB/CZMlGIWoo0nh6LyPmAiEYuxjauPb0lBRwPXOfU6LkHWQwtO0NooHfzYfcr8xeoq3vggTvgv+XTzdulahteZI6RHSj6djPWnFesLK4dOHAWMR/biIcSSOS2RCxATGRowlqMuFO9y8kLqf9fvz+cOMMcQO8nwhZSG8SaPhqJT6M5CJMbmmuJFjazYMmKO1fqqV2ihE02V/YQyfHt4EwZFwy19h4n3g37Xxc71cXmkeW45vYfPxzWw5voVjZ44BYOpqIqZXDLOHzyamVwwjwkc0upXaml3H+cenmcyxDOC+y9v//VYhmuui4aiUekRr/cumXqzOhgE15/77EtsnRNPkfGUMnx76CoL6ws1/Me4p+nAoFpwt4OtjX7P52GY2H9/M0dKjAIR2DcXSx8LDYx5mUp9JmEPMzdpR5uCJ0/x8yTbGDzDx+ztGy240QtTjouF4KeEmwSjcIm8bpP7G6DH26G0s3I95EDr73v2zs1VnyTiRwaa8TWw6tomDRQcBCO4SjKW3hftH3c+kPpOINkW3eClF8dlKEt7OIKBLJ/51XwzdOnestZxCNFVjPcfBMtFGeKWSPFj3B9j+HgSGwQ3zwfIQdA7wdMuazKEd7D21l03HNrEpbxNbT26l0lFJF78uTOw9kZ/E/IQp/aYwImxEq6wrdDg0P128jdzCMt6dO5m+Ib7zZyWEuzV2z3EBxoJ+IbxDxRnY+A9jraKjCqb+GK76OXQL8XTLmsR+zs6GvA18dfQrNuZtpPCcMa9teOhwvjvyu0zpO4WJvScS4N/6wfVi2gE+3XeSZ+8YzWVDZAKOEBfTWDjGK6VmAWlAKpAmM1GFRzgcsGMxrHsWTufBqBkQ9wyEDfF0yy7KoR3sObWHL49+yVdHv2JXwS4c2kFo11CuiLyCqf2mMqXfFHoG9GzTdqzdfZyXPs0kPrY/908e1KavJUR70Fg4pgCLgUnA3cBCpZRGwlK406GNsPZpY8u3fhNh1uswaIqnW9Wg0xWn2XB0A18c+YINeRsoPFeIQjGm5xgSxyVyZeSVjA4fTSc37d164MRpfrZ4G+P7h/CHGWNkAo4QTdBYOG6pqesItUs1JgFxfBuWRcA8rfVrbdpS0fEUZkPqb40SUsGRcGcyjI33yuLCuadz+Tz3cz478hkZxzOo0lWYupq4ot8VXBl5JVMjpxLWzf1DmfayCub+J53Arv4sut8iE3CEaKLGZqs+7/J9MUavMa3mMeew6y+VUtla60/bpJWiYzlXDF88D98sAj9/uPZXxqbgXbyn/l+1o5qdBTv5LPczPj/yOZn2TADMIWbuH30/1/S/hvER493WO6xPVbWDx9/bSp79LP9LmEyfEN+bwSuEp1zy3qpa6xQgRSk1H5BwFC3nqIaMN4xF/GWFxo421/0agvt6umWAsdRiY95GPsv9jC+OfEHhuUI6qU7E9o7lF5ZfcM2AaxgY7D0L6heu3c+XBwtYMHOs7IAjRDO15sbjUg9KtNyx7bDyCeO+4uCr4IY/Qd/xnm4VJRUlfJ77OZ8e/pQNeRs4W3WWoM5BXNn/Sq7pfw1TI6cS0tX7Zsqu2HqU5C9sPDBlEHMmeU9gC+ErLjkclVIHgSSgdQtDio6hvNToKX7zCgT2hJmvwZiZ4MFJI/ll+azPXc+6w+vYfGwzVbqKiIAIbo+6nWkDp2HpY2l0ezZP2nmkmKRlO7h8SBi/uXWUp5sjhE9qjZ7j88D1wJ9b4VqiI9n3MXz8Cyg5ApaHYdrvIMDkkabkluSy7vA61h1ex/b87Wg0A4MGcv/o+5k2cBpje471iQK/+afLSXg7nZ49uvLP78bQuZP3t1kIb9Qa9xyTW6MhogMpPgqr58G+VdBrFMz6BAZe7vZmZBdn80nOJ6QeSmV/0X4ARoaN5IcTfsi0gdOINkX71LKH8qpqEt9Ox15WScoPphDew3f3lRXC06TYsXAfRzVsfhU+/YPxddwzxixUN9ZWtBXb+CTnEz459Ent3qUTe03kF5ZfcN3A6+gf1N9tbWlNWmt+9f4urIftvPLdGEb38777oEL4kgbD0bmmcZrWenlLLqyUugtjk4CSljZOtCN522DVT4wJN9FxRtUMN+1uY7PbWHtoLZ/kfFK75CKmVwy/vOyXxA2Mo3f33m5pR1v695fZpGQc4YlpQ7lprHfM7hXClzUYjlrrYqXUOmc9x0+auoZRKTUNY5OA+RKM4oIJN7Neh9F3tfmEmyx7Vm0PMdOeiUIxsdfEdhWINdbvP8n81Xu5aUwfnpg21NPNEaJdaGwTgGKMBf7TlFL/wpiRmgHYXA41AxYgFFgkhY4FAPs+go/nuW3CTXZxNmuy19QbiNMHTadXYK82e21PyTx5mh+/u5URfYJ5YfZ4/Px85x6pEN6sSfcctdbrgHVQ2zMMxQhEADtGWC51hqno6IqPwOok54Sb0RD/Bgy4rE1e6ljpMdbkrGF19mr2Fu5FoYjpHcNTlz1F3KC4dhmINexlFTzyVjpdO/vx6oMWArvIFAIhWkuz/zU5g1KIC2kN6a8b+6E6qiHu9zDlR60+4ebU2VN8cugTVmevZuvJrQCM7TmWeZPmccPgG9p1INaorHbw2LtbOWo/y3tzJxNpktqMQrSmJoejUmoCEAak172XqJQaAgwBbE0tjKyUMgEJGD1OM8bEHetFjjdj1JZcpLWuu6/rPCAco3JIGBCvtU5s6s8kWlHxEfjwccj6FMzXwG1/h9DBrXb5kooS1h1ax5qcNXxz7BuqdTXRpmgen/g4Nw2+iQHBA1rttbyd1ppnPtzNV5kFLJw1Dstg2RpOiNbWaDg6QzGFb4dRtVLKCszVWm/TWmcD2UqpUxhB1RRLgUSttc35GqlKqXittb2e149zfml2fc4pwfmRBsxt4uuL1qI1bP+fMYzqqIJb/mrcX2yFCTdnq87y+ZHPWW1bzZdHv6TSUUlkj0geHvMwNw65kWGhw1rhB/A9b2zI4Z1vDpP4HTOzLR3nPwVCuNNFw9G5nOPfGEG2rs5jccDTztqOSc4eY1FTXtDZazTXBKOTzXnNFNfja3qKSqnCei5n11qHNuV1RRsoPQkrfwL7P4KBV8CMlyGsof/DNE2lo5KNRzfycfbHrM9dz9mqs/QM6Mmc4XO4achNjO051qcW5re2T/ed4I8f7eH6Ub1JunGEp5sjRLvVWM9xLsZax9qJNs6vlwHLaoJSKRULNDgs6sKCMYmnLjswnXrCsSmUUjEYQek6i1a0ld3vw6qfQcUZuP5PMPkH0MLyTFprdhTsYFXWKtbkrMFebie4SzA3D7mZm4fcTGzvWI+WfvIW+46X8Pi7WxnZN5gX754gM1OFaEONhePWi81ArROU1PncGBPg2gs8RcPDphflrCeZhhHSiVrrpJZcRzRRWSF8/CTsWgb9YuDOf0HE8BZdKvd0Lqtsq/jI9hGHSg7RtVNXrh1wLbeYb2Fqv6l0duPOOd4u/3Q5338znR7d/HntwUkyM1WINtbYv7C2qrTRKjMIXPZ1TVFKLVBKpdadtFMjPz8fi8VS+31CQgIJCQmt0YyOY/9qo6xUWaFRZ3HqT6FT835JF5cXszZnLSuzVrItfxsKxaQ+k/j+mO8TNyiOoC5BbdR433Wuspq5/0nn1JlyliZeIUWLhWih5ORkkpPP2w68Z0PHNvabzdQaDXJhr+e64VzYm2yUUirGZZarFWN49oJwjIiIID1dSk62yLliWPM0bPsv9B4D9y2DPmObfHpFdQVfHPmClVkr+eLoF1Q5qogKieKJmCe41Xwrfbr3acPG+zatNU8u3c62XDv/ui+Gsf1lz1QhWsq1U6SUKmjo2MbCMUopFdzQNnBKqcHALCAb+KXWelIT2pfOhT1HE5DahHPrvnYMxsYEdSfkmICs5lxHNCJrPXzwGJzOg6uehKuTwL9Lo6dprdl6cisrbStZm7OW0xWnCe8Wzj0j7uE2822MCBvRoSfWNNXfUg+wascxkm4cwY1jZM9UIdylse3jnldKZSql5mMEUSHGvcE44G4gtWarOOcerI3SWtuVUulKqbozVi0YBZNr1jTS2OQarbVVKeV6f9EMLGlKO0QjKs4Yi/m3/BvCh8L3U6G/pdHTcopzWGlbyUe2jzhaepQA/wCuG3gdt5lv4/K+l+PvJ/fKmmpJei4vfZpJfGx/Hr360mYBCyGapym/qSwYw5SvYtyDtGMEUO0sVues1aaucQSIBxKUUjaMXuTcOmscEzF6gInOa8dghLEFSHKGas2gcbpzIwA7EIWxCUDNdURLHdoEK34ARTkw+Ucw7TfQueEdWOzn7Hyc/TGrbKvYWbATP+XH5X0u50cTfsS0gdMI7Bzovra3E18ezOfp5Tu5amhPnrurYy9fEcITlNZNm3PjDMAw56L/uo8PAUK01ttav3mtx2KxaLnn2IjKc0atxU0vg2kgzHgFBk+t99AqRxUb8zayInMF63PXU+WoYnjocG6Luo2bhtzUIbZwayt7j5UQ/69N9A8NYOmjUwjqJrN2hWgLSqkMrXW9Q2JNHuNy9hIvWNbhGpbCRx21wvuPQsF+Y4eb6X+Arj0uOMxWbGNF5gpWZa0i/2w+oV1DuXv43cyInsHwsJYt6RDfOlZ8lofe2EKPrv688dAkCUYhPERuAHV0jmrY8HdY/yfo0RvuWw7R08475HTFaVZnr+aDrA/Ykb+DTqoTV0VexYzoGXyn/3dkPWIrOX2ukofe2EJpeRVLEqfQN0Q2ExfCUyQcO7LiI0ZvMedLGH0n3Po3CDAm/zq0g2+OfcOKzBWsO7yO8upyok3RPGl5klvMt9AzoMHlQaIFKqsd/PAdKwdPlvLG9yYxql+wp5skRIcm4dhR7V5hLOivroQ7/gkT7gWlyC3J5YOsD/gw60OOnTlGUJcgZkTPYEb0DEaHj5aJIW1Aa82v39/FlwcLWDhzHN8ZFuHpJgnR4bVaOCqlntRa/6W1rifaSHmpUUFj238hMhbuepWy4L58kvUBKzJXkHEiA4ViSr8p/DT2p1w38Dq6durq6Va3a//4NJPF6bn8+LpoZk+SKhtCeINmhaNSahpGXcUhrk8BIYCEozc7mgHLHoHCbPSVP8c6+gZW7HuTtTlrOVt1loFBA/nxxB9zW9RtsmuNm7y3+TB/TT3AXTGR/HR6xyzBJYQ3am7PcTouVTpqKKV+0TpNEq3OUQ0bXoT1z3E8pC8fXvsjPjj1DYc/WUqgfyA3Dr6RGdEzmNhrogybutGaXcf51fs7uXZ4BAtmjpM/eyG8SHPDMbWhKh1a6+dboT2itRUfoXx5Ap8WbOX9ISP42nEanfMhlt4WEsYlMH3QdFmk7wFf207x4/9tZfwAEy9/N4bOnfw83SQhRB3NDccGdwxQSl2ntf70EtsjWtH+La/w/pa/sTKgMyW9etI3oDuJ0fdye9TtDAiSe1uesievhLlvpTMwLJDXpfyUEF6puf8qr3fuZ2rH2GfV7nw8FJgGDG21lokWKa0oZfXB91m+9RV2VZ+mc2BXpkVeyZ2j72dy38n4KemheFJuYRkPvrGZHt38+c/DlxHavfFN3IUQ7tfccIwD5vNtKNYw0cJixeLSaa3Zlr+N5QeXs9b2MWcdFURXVDKv1yRujXuB0O6ylZs3KCgt5/7XvqGiysG7j06hn0kW+QvhrZobjnO11lvre8K5ibhwo1NnT7HKtoplB5eRXZxNoPLn5pJi7nIEMvb211AN7Isq3K+0vIqH3tjC8ZJzvPPIZIb2lqLOQnizZoVjQ8HoFHqR50QrqXZUs+nYJpYfXF674ff40JE86wjlhsM7CBx9F9zyVwgwebqpwqm8qprEt9PZc6yEVx+IJXaQ/FMRwttdNByddRybQmHcc2xKsWPRAnmleazIXMH7me9z/MxxTF1N3DviXu4iiKh1zxnLNe54BcbfDbIkwGtUVTt44r1tbMg8xV/ix3PdiN6ebpIQogka6znGAou48B5jfeSeYyurqK5gfe56lh9czqa8TQBM6TeFJy1Pcm0vC13W/hq2vwv9J8FdyRAmfwXepNqheXLpdtbsPs7vbhvFrNj+nm6SEKKJGgvHpEaGUmvJPcfWk1mUyfLM5azKWkVReRF9uvfh0fGPMiN6Bv169IMj6fDqNLAfgu/Mg6vngVTG8Cpaa369YicrtuXxixuG89BU102lhBDe7KLh2NRgdB4rdR0vQVllGWtz1rLs4DK252/H38+fawdcy11D72JK3yl08utkDJ1+8Tysnw/B/eB7H8OgKZ5uunChtebZVXt4b3Muj10bzY+ujfZ0k4QQzSSrjz1Ia83Ogp0sP7ic1dmrKasqY0jIEJ60PMmt5lsJDwj/9uCiQ0Z5qcMbYcwsuOUFmXTjpf7yyX7e2JDDQ1MH8/PrZb9UIXyRhKMHFJ0rYpVtFcsPLifTnkmAfwA3DL6BmUNnMj5i/Pl7bGoNO5bAx08aX9+5CMbNkUk3Xurl9Zm8vD6Ley4bwG9vHSX7pQrho1qzZNV8rfVTrXW99qamePDyg8tZd3gdlY5KxvYcy2+n/JabBt9Ejy49LjzpbBGs+hnsXg4DpxjBGDrI/Y0XTfL6V9k8v3Y/Myb0448zxkowCuHDGlvKkclF9lOtIxyjZJWEo4uSihLe3fsuKzJXcLT0KMFdgpk9fDZ3Rt/J8LDhDZ9o+xxW/ABKT8C038LUn4BfJ7e1WzTP/zYf5tlVe7hxdB/+Ej+eTn4SjEL4ssZ6jmla60drvlFKTYQLJ+o46zwWtn7zfJ8ffryx6w3GRozliZgnGi8eXFUO656FTf8H4UPhkTToN9F9DRbNtmLrUZ56fyfXDI/gpXsm4i8VNoTweY3NVn3U5aHQ+ipvaK3XKaWua9WWtRM9uvTgk1mfENI1pPGDT+yB5XPhxC6Y9AhM/wN0kXJS3my59QhPLt3O5CHh/Ou+WLr4SzAK0R7IhBw3aDQYHQ745l+Q9gx0C4Z7l8CwG9zSNtFyS7bkkrR8B1dEhfPvBybRrbMMewvRXjQ3HGOVUoVa6211H1RKDcbYTUfqOTZXSZ5xb9H2GQy7CW7/B/SI8HSrRCPe/eYwT7+/k6uG9uTVBywSjEK0M83dePx5pdQSpdQQoGZHHDNg01rPafXWtXe7V8DKJ6C6Am59EWK/J0s0fMDbm3L4zQe7uXZ4BK/cFyvBKEQ71OxhVa31bGc4xjgfssruOM10rgTW/BK2vQP9YuCuV6Gn7KLiC97YkM3vV+4hbmRvXv7uRLr6SzAK0R616J6jMwwlEFvi8NewPAGKc2VfVB/z6hc2/vTxXm4c3YeX7pkok2+EaMdkQo67VFfC5wvgyxcgZAA8tAYGXu7pVokm+udnmSxcs59bxvblxbsn0FmWawjRrjUYjs6h0wU0bROAC053npektc5pWdPakYJMWP4I5G2FCffBjfONWanCJ/xj3UFeSD3AHRP68UL8eFnHKEQH0GA4OodOZ7uxLe2P1pDxBqz9Ffh3hfi3YPQMT7dKNJHWmhfTDvL3dQe5KyaS52fJzjdCdBQeGVZVSpmABIwZr2aMnXisFznejNGLXaS1TmvpddyqNB8+fAwOrAHztTDjn0aZKeETHA6j7NSbG3OYbenP/LvGSTAK0YF46p7jUiBRa20DUEqlKqXitdZ21wOVUnHOL+src9/k67jV/jVGMJ4rgRv/DJclgp8MxfmKiioHP1+6nZXb83jkyiE8ffNI/CQYhehQ3B6Ozt6euSbQnGxAHJDienxNT1Epdd7erc29jltUnIFPfg3pr0PvMfDAh9B7lEeaIlqmtLyKR9/O4KvMAn550wgSv2OW6hpCdECe6DlaALvLY3ZgOs0Ltda6Tus4ajX2RT2VBVc8Dtf9xrjPKHxGQWk5D7+5hd15JTw/axzxlgGebpIQwkM8EY4mLqzgcYr6h01b7Tr5+flYLJba7xMSEkhISGjmS9bDUQ1f/RU++zP06A0PfADmqy/9usKtcgvLeOD1zRwrPkvy/bFMG9nb000SQrSy5ORkkpOT6z7Us6FjPXXPMczd14mIiCA9Pb2VXtapKAeWJ0Lu1zD6Lrj1rxAQ2rqvIdrcnrwSHnxjMxVVDt555HJiB7XW21MI4U1cO0VKqYKGjvVEONoxen11hdP8epCtdZ3m0xq2vwcfzzP2Qr3rVRgbL/ui+qCvbaeY+1Y6Pbr58+6jUxjaO8jTTRJCeAFPhGM6F/b4TECqh67TPGWFsOqnsGcFDJoKd/4LTAPb9CVF21iz6xg//t82BoYF8p+HL6OfKcDTTRJCeAm3ry9wLrNId65drGEBamalml2ea9F12kTWenjlCtj3EcQ9Aw+ulGD0Ue98c4gfvmNldL9gliZOkWAUQpzHU/cc44EEpZQNo/c3t87axESMHmAigFIqBmN5hgVIUkqZtdbJTbhO69r2Hqx4FHoOg3v+B/0mtMnLiLblcGheTDvAS59mcu3wCF7+bgyBXWSLYSHE+ZTWLdk61fdYLBZ9SRNyygph40tGJY0uga3XMOE2ZRVV/HzJdlbvOk58bH+eu2usbCAuRAemlMrQWlvqe07+y9xUgWHGUKrwSUeKypj7nwz2Hy/h17eM5PtXDpHF/UKIBkk4inYvPaeQxLczqKhy8Nr3JnHt8F6ebpIQwstJOIp2bcmWXH61Yif9QwN59QEL0b16eLpJQggfIOEo2qWqagfPfbyP1zdkc9XQnvzfPTGEBHb2dLOEED5CwlG0O8VllTz2npUvDxbw0NTB/OrmkVKgWAjRLBKOol3Jyi9l7lvp5BaVsWDmWOZMknWoQojmk3AU7cbnB/J57F0rXTr58e7cyUwaLHukCiFaRsJR+DytNa99lc1zH+9leJ9gXn0glv6hshZVCNFyEo7Cp50pr+I3K3axfOtRbhzdhxdmj6d7V3lbCyEujfwWET5r19FiHn9vK4dOneEncUP58XVD8fOThf1CiEsn4Sh8Ts0w6oI1+wjv3pV3505msjnc080SQrQjEo7CpxSUlvOLpdtZvz+f6aN6s3DmOEK7d/F0s4QQ7YyEo/AZGzIL+MnibRSfreTZO0Zz/+RBsj+qEKJNSDgKr1dZ7eCvqQf41+dZREX04D8PX8bIvsGebpYQoh2TcBReLbewjMff28q2XDv3XDaA39w6SuovCiHanPyWEV7rw+15/Gr5TlDwf/dO5NZx/TzdJCFEByHhKLxOWUUVz3y4myXpR4gZaOLvd09kQJgs6hdCuI+Eo/Aqu/OMtYvZBWd47NponogbSmfZNFwI4WYSjsIrVDs0b27MYcHqfZgCO/PO9y/niuienm6WEKKDknAUHrfraDG/en8n248UM21EL56PH0+YrF0UQniQhKPwmDPlVfw19QBvbMgmrHsX/n73BG4f30/WLgohPE7CUXjE2t3HeebD3RwrPse9lw8k6YYRhAR29nSzhBACkHAUbpZnP8vvPtxN6p4TjOgTxP/dG0PsoFBPN0sIIc4j4SjcoqrawZsbc/hr6gG0hqduGsHDVw6RmahCCK8k4Sja3LZcO08v38meYyVcN6IXv799tKxbFEJ4NQlH0WZKzlXy/Jr9/PebQ/QK6sor343hxjF9ZMKNEMLrSTiKVqe15qOdx/j9yj2cKi3nwSmD+fn1wwjqJhNuhBC+QcJRtKrDp8r4zQe7+PxAPmMig3ntQQvj+ps83SwhhGgWCUfRKiqqHLz6pY2X1h3E30/x21tH8cCUQfjLhBshhA+ScBSXRGtN6p4TPL92PwdPlnLj6D787vZR9A0J8HTThBCixTwSjkopE5AA2AAzkKa1tjb3WKXUPCAcWAyEAfFa68S2br8wQjFt70leTDvA7rwSBocH8tqDFqaN7O3ppgkhxCXzVM9xKZCotbYBKKVSlVLxWmt7C45NcH6kAXPbvOUdnNaaT/ed5MW0g+w8Wsyg8ED+Ej+eGRP6yRCqEKLdcHs4OnuC5pqwc7IBcUBKM4+1a61lexU30Fqzfr8RijuOFDMwLJDnZ43jzomREopCiHbHEz1HC2B3ecwOTMclHJt6rFIqBiMobYhWpbXmswP5vJh2kO25dgaEBbBw5jjujImU3W2EEO2WJ8LRBBS6PHYK435is49VSs3CGFKNU0olaq2TWq2lHZjWms+dobgt107/0AAWzBzLXTH9JRSFEO2ep+45hrXGsVrr5DrfpiilFiilUrXWaa7H5ufnY7FYar9PSEggISGhGc3oGLTWfHmwgL+lHWDrYTuRpgDm3zWWmTH96eIvoSiE8F3JyckkJ9eNDRqsqO6JcLRj9AjrCufCHmKjxyqlYlxmuVoxhlwvCMeIiAjS09Nb0t4OQWvNhsxT/C3tABmHiog0BfDcnWOZFSuhKIRoH1w7RUqpgoaO9UQ4pnNhb9AEpDbnWOd9xnVAqMtzWa3RyI5Ca83GrFO8mHaALTlF9A3pxh9njCHe0p+u/p083TwhhPAIt4ej1tqulEpXStWdhWoBkgCUUmbncbaLHet8zvX+ohlY4o6foz3YmFXAi6kH2ZxTSN+QbvxhxhhmSygKIYTH7jnGAwlKKRtGz3BunXWLiRg9wMQmHJvu3AjADkRhbAJQ85xowCZnT/Gb7EJ6B3fl2TtGM2fSAAlFIYRwUlprT7fBLSwWi+7I9xzPVVazbu9J3v46h69thfQK6soPr4ni7ssG0q2zhKIQouNRSmVorS31PSd7q7ZjWmu25dpZZj3Cyu3HKD5bSd+QbvzutlHcI6EohBANknBsh44Vn2W59SjLrEew5Z+hW2c/bhjdh1mx/bkiqied/KTYsBBCXIyEYztxtqKatbuPs8x6hK8yC9AaLhscRuJ3zNw8tq8UGhZCiGaQcPRhWms2ZxeyzHqEj3cep7S8iv6hATx+3VBmxkQyKLy7p5sohBA+ScLRB+UWlrHMeoTl1qMcLiwjsEsnbh7bl5kx/bl8SBh+MmwqhBCXRMLRR5SWV/HxzmMsyzjCN9mFKAVXRIXzk7ih3DimD4Fd5K9SCCFai/xG9WIOh2aT7RQpGUdYs+s4ZyurGdKzO09eP4w7Y/oTaQrwdBOFEKJdknD0Qrb8UpZZj/C+9Sh5xecI6ubPnTGRzIzpT8xAE0rJsKkQQrQlCUcvUVxWycodeSyzHmHrYTt+Cr4zLIKnbh7J9FG9ZU2iEEK4kYSjB1VVO/jyYAEp1iOk7jlBRZWDYb178NRNI5gxMZLewd083UQhhOiQJBw9YN/xEpZlHGHFtjzyT5cTGtiZey8byMyY/oyJDJZhUyGE8DAJRzdLzylk1r824e+nuHZEL2bG9Oe6Eb2kZqIQQngRCUc3mzgwlD/OGMNNY/oQ3qOrp5sjhBCiHhKObtbJT3Hf5EGeboYQQoiLkLE8IYQQwoWEoxBCCOFCwlEIIYRwIeEohBBCuJBwFEIIIVxIOAohhBAuJByFT0pOTvZ0E4RoMnm/+h4JR+GT5JeN8CXyfvU9Eo6tbOXKlV55zZZco7nnNPX4phzXFn+O3sxTP297er829zx5v7Zce3q/NkTCsZW1p182Eo7u055+2Ug4tn/t6f3aEKW1dtuLeZJSKh845IaXCgGKvfCaLblGc89p6vFNOa6xY3oCBU1sly9oi/eNp17XU+/X5p4n79eWay/v10Fa64j6nugw4SiEEEI0lQyrCiGEEC4kHIUQQggXEo7CZyilzEqppUqpOE+3RYj6NPQeVUqZlFLzlFKznJ9jPNVG0TRSz1H4hDq/bMwebYgQDWjkPboUSNRa25zHpiql4rXWdne1TzSPhKPwCVrrNAClVKGn2yJEfRp6jyqlTIC5JhidbEAckOK2BopmkWFVIYRoWxbA7vKYHZju9paIJpNwFEKItmUCXEc8TgFh7m+KaCoJRyGEaHsShD5GwlEIIdqWHaP3WFc4F/YmhReRcBRCiLaVzoU9RxOQ6v6miKaScBRCiDbkXK6RrpSqu8TDAqR5pkWiKWQph/AJzkXTcRi/VJKUUmattRTJE16jkfdoPJCglLJh9CLnyhpH7yYbjwshhBAuZFhVCCGEcCHhKIQQQriQcBRCCCFcSDgKIYQQLiQchRBCCBcSjkIIIYQLCUchOhhnsd2lSint/Cw1MoVwIescheiAlFKzgKVaa+XptgjhjaTnKETHNAnZvkyIBkk4CtExxSEbXwvRIAlHIToYpZQJiEF6jkI0SMJRiI7HAqC1tnq6IUJ4KwlHITqe6UivUYiLknAUouOR+41CNEKWcgjRwSilNBArw6pCNEx6jkK0Y0qpBKVUVp3vZwF2CUYhLk7CUYj2LRZIqvP9Uy7fCyHqIcOqQrRjzq3hYpzfTgJStdYyGUeIRkg4CiGEEC5kWFUIIYRwIeEohBBCuJBwFEIIIVxIOAohhBAuJByFEEIIFxKOQgghhAsJRyGEEMKFhKMQQgjhQsJRCCGEcPH/jawIxJ4MGXEAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 504x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"figsize(7,5)\n",
"for pot,label in zip([iso,mwp14_v,mcm_v],\n",
" ['Isothermal','MWPotential2014','McMillan17']):\n",
" ro, vo= get_physical(pot)['ro'], get_physical(pot)['vo']\n",
" # The following Emax is in (km/s)^2\n",
" Emax= evaluatelinearPotentials(pot,4.*u.kpc)\n",
" Es= numpy.linspace(0.*u.km**2./u.s**2.,Emax,101)\n",
" # Setup Orbits to easily pass the required (z,vz) to the actionAngleVertical calculation\n",
" orbs= Orbit([[0.,0.,0.,0.,numpy.sqrt(2.*(E-evaluatelinearPotentials(pot,x=0.)))\\\n",
" .to_value(u.km/u.s)/220.]\n",
" for E in Es])\n",
" # The following action J and frequency OJ are in internal units of ro*vo \n",
" # [different for McMillan17!] \n",
" J= numpy.array([actionAngleVertical.actionAngleVertical(orb,pot=pot).Jz() for orb in orbs])\n",
" OJ= numpy.array([2.*numpy.pi/actionAngleVertical.actionAngleVertical(orb,pot=pot).Tz() for orb in orbs])\n",
" # Get the derivative simply using finite differencing\n",
" dOJdJ= numpy.gradient(OJ,J)\n",
" J= J[:-2]\n",
" OJ= OJ[:-2]\n",
" dOJdJ= dOJdJ[:-2]\n",
" # We now convert the action to Scott's units of 30 km/s x 0.3 kpc\n",
" semilogx(J*ro*vo/(scott_ro*scott_vo),\n",
" numpy.fabs(dOJdJ*J/OJ),\n",
" label=rf'$\\mathrm{{{label}}}$')\n",
" from scipy.interpolate import InterpolatedUnivariateSpline\n",
" print(\"f at J=1.79 for {} is {}\"\\\n",
" .format(label,InterpolatedUnivariateSpline(J[2:]*ro*vo/(scott_ro*scott_vo),numpy.fabs(dOJdJ*J/OJ)[2:])(1.79)))\n",
" # Save the data for the paper\n",
" numpy.savetxt(f\"dlnOdlnJ-{label}.dat\",(J*ro*vo,numpy.fabs(dOJdJ*J/OJ)),delimiter=',')\n",
"xlabel(r'$J$')\n",
"ylabel(r'$|\\mathrm{d} \\ln \\Omega(J) / \\mathrm{d} \\ln J|$')\n",
"gca().xaxis.set_major_formatter(FuncFormatter(\n",
" lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\\\n",
" numpy.maximum(-numpy.log10(y),0)))).format(y)))\n",
"legend(frameon=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Scott's time-to-spiral-dissolution calculation:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.8549879733383469 Gyr\n"
]
}
],
"source": [
"f= 0.2\n",
"P= 100*u.Myr\n",
"to= 10.*u.Gyr\n",
"print(((2.*f)**(-2./3.)*P**(2./3.)*to**(1./3.)).to(u.Gyr))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.9.13 ('py39')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.13"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "a82f9eb2f595400cd0f024a877274681853baaeedd921355b08c40b0424842e5"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment