From c954e4febc52de09d1b872347a83fad35080ba82 Mon Sep 17 00:00:00 2001 From: Chris Papageorgakis Date: Fri, 14 Jul 2023 20:35:32 +0200 Subject: [PATCH] all changes for draft 1 --- annealing_depth.py | 87 +++++++++++++++ notebooks/oxygen_penetration_model.ipynb | 135 +++++++++++++++++++++++ pre_irr_index.py | 4 +- 3 files changed, 224 insertions(+), 2 deletions(-) create mode 100644 annealing_depth.py create mode 100644 notebooks/oxygen_penetration_model.ipynb diff --git a/annealing_depth.py b/annealing_depth.py new file mode 100644 index 0000000..cae3b4b --- /dev/null +++ b/annealing_depth.py @@ -0,0 +1,87 @@ +import matplotlib.pyplot as plt +import mplhep as hep +import numpy as np +from iminuit import Minuit +from iminuit.cost import LeastSquares +from jacobi import propagate + +# from scipy.special import erfc + +plt.style.use(hep.style.ROOT) + + +def sqrt_law(x, a): + return a * np.sqrt(x) + + +if __name__ == "__main__": + days = np.array([1, 5, 11, 18, 25]) + edays = np.ones_like(days) * 0.5 + + time_array = np.linspace(0, 32, 300) + + z_index = np.array([1.23, 1.23, 1.24, 1.16, 1.14]) + z_color = np.array([1.59, 2.09, 2.68, 3.32, 4.09]) + + ez_color = np.array([0.02, 0.04, 0.07, 0.11, 0.21]) + + z_index_mean = np.mean(z_index) + z_index_std = np.sqrt((np.std(z_index)) ** 2 + (0.05 * z_index_mean) ** 2) + + fig, ax = plt.subplots(figsize=(7, 6)) + plt.errorbar( + days, + z_color, + yerr=ez_color, + xerr=edays, + fmt="o", + color="black", + ecolor="black", + label="color depth", + capsize=3, + ) + ax.hlines( + z_index_mean, + -2, + 32, + color="black", + linestyles="dashed", + label="index boundary", + linewidth=2, + ) + ax.fill_between( + np.linspace(-2, 32, 100), + z_index_mean - z_index_std, + z_index_mean + z_index_std, + alpha=0.2, + color="black", + ) + plt.xlabel("Days after irradiation") + plt.ylabel("Annealing depth $z$ (mm)") + plt.xlim(-2, 32) + plt.ylim(0, 5) + + least_squares = LeastSquares(days, z_color - z_index, ez_color, sqrt_law) + m = Minuit(least_squares, a=0.1) + m.migrad() + m.hesse() + y_1, ycov = propagate(lambda p: sqrt_law(time_array, *p), m.values, m.covariance) + y_1 += z_index_mean + yerr_prop = np.diag(ycov) ** 0.5 + ax.plot(time_array, y_1, label="fit", color="red", lw=2) + ax.fill_between( + time_array, y_1 - yerr_prop, y_1 + yerr_prop, facecolor="red", alpha=0.2 + ) + + # Print the fit results + print("\nFit results:") + print(f" chi2 / ndof = {m.fval:.1f} / {m.ndof:.0f} = {m.fmin.reduced_chi2:.1f}") + for p, v, e in zip(m.parameters, m.values, m.errors): + print(f"\t{p} = {v:.4f} ± {e:.4f}") + print(f"\tD = {v**2:.4f} ± {2*v*e:.4f}") + print() + + plt.tight_layout() + plt.legend(loc="upper left") + plt.savefig("plots/annealing_depth.pdf") + # plt.show() diff --git a/notebooks/oxygen_penetration_model.ipynb b/notebooks/oxygen_penetration_model.ipynb new file mode 100644 index 0000000..9156f4f --- /dev/null +++ b/notebooks/oxygen_penetration_model.ipynb @@ -0,0 +1,135 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import mplhep as hep\n", + "from scipy.special import erfc" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "hep.style.use(hep.style.ROOT)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "def ox_diffusion(t, D, z_idx):\n", + " n_order = 10\n", + " u_critical = 0.05\n", + " l = 10 - 2 * z_idx\n", + " x = np.linspace(-l / 2, l / 2, 100)\n", + " u_n = lambda n, D, t: ((-1) ** n) * erfc(\n", + " ((2 * n + 1) * l / 2 - x) / (2 * np.sqrt(D * t))\n", + " ) + ((-1) ** n) * erfc(((2 * n + 1) * l / 2 + x) / (2 * np.sqrt(D * t)))\n", + " depths = np.empty_like(t)\n", + " for i, t_i in enumerate(t):\n", + " u = np.zeros_like(x)\n", + " for n in range(n_order):\n", + " u += u_n(n, D, t_i)\n", + " depths[i] = l / 2 - abs(x[np.argmin(u - u_critical)]) + z_idx\n", + " return depths" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "time_array = np.linspace(0, 32, 100)\n", + "z_index = np.array([1.23, 1.23, 1.24, 1.16, 1.14])\n", + "z_index_mean = np.mean(z_index)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/4f/rfs4xk8j61x737wyrrknwd980000gn/T/ipykernel_13536/1074346541.py:7: RuntimeWarning: divide by zero encountered in divide\n", + " ((2 * n + 1) * l / 2 - x) / (2 * np.sqrt(D * t))\n", + "/var/folders/4f/rfs4xk8j61x737wyrrknwd980000gn/T/ipykernel_13536/1074346541.py:7: RuntimeWarning: invalid value encountered in divide\n", + " ((2 * n + 1) * l / 2 - x) / (2 * np.sqrt(D * t))\n", + "/var/folders/4f/rfs4xk8j61x737wyrrknwd980000gn/T/ipykernel_13536/1074346541.py:8: RuntimeWarning: divide by zero encountered in divide\n", + " ) + ((-1) ** n) * erfc(((2 * n + 1) * l / 2 + x) / (2 * np.sqrt(D * t)))\n", + "/var/folders/4f/rfs4xk8j61x737wyrrknwd980000gn/T/ipykernel_13536/1074346541.py:8: RuntimeWarning: invalid value encountered in divide\n", + " ) + ((-1) ** n) * erfc(((2 * n + 1) * l / 2 + x) / (2 * np.sqrt(D * t)))\n" + ] + }, + { + "data": { + "text/plain": [ + "(0.0, 5.0)" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzAAAANFCAYAAACk5xlpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABoCUlEQVR4nO3dd3hUZcL+8XvSaaHXhN4DKF1RFFGkGRR3XfvaNcAICIqCDRAUbKDiIKC7+upa310XlaYIrIguoCggXTqZwNATCKTO+f3BjzG8JGQmmZkzZ+b7uS6u6xzPM09uOB7IndNshmEYAgAAAAALiDI7AAAAAAB4iwIDAAAAwDIoMAAAAAAsgwIDAAAAwDIoMAAAAAAsgwIDAAAAwDIoMAAAAAAsgwIDAAAAwDIoMAAAAAAsgwIDAAAAwDLKVGCuvfZa2Ww2r3699957fo4MAAAAIFKVqcBs27bN3zkAAAAAoFQxvn4gJydH+/btkyQNGzZMgwYNuuD4Dh06lC0ZAAAAAPwfPheYHTt2yDAMSVJqaqr69+/v91AAAAAAUByfLyH7/fffPcutWrXyaxgAAAAAuJAyF5jY2Fg1btzY74EAAAAAoCQ+X0J2tsA0b95cMTExKiws1L59+7R7925VqVJFLVu2VGJiot+DAgAAAIDPZ2DOPoGsfv36mjJliurUqaOmTZuqd+/e6tq1q2rWrKkBAwZo5cqVfg8LAAAAILLZjLN35HspKSlJGRkZpY6LiorSa6+9puHDh5c5HAAAAAAU5VOByc7OVuXKlT3rnTp10ksvvaRu3bqpoKBAv/32m2bNmqVPP/3UM2b+/PkaOHBgsfNVqlRJOTk5io6OVu3atcv+m7DZyvxZAAAAAP7h47mRcxw6dEiFhYVKSEhQdnZ2ieN8KjDr169X586dJUkDBw7URx99dE6hOWv69OkaPXq0JKlNmzbauHGjoqLOv1otOjpabrfb2y8PAAAAIMxFRUWpsLCwxO0+X0LmjcLCQl1yySVas2aNJGnNmjWe4lNUXFyc8vPzFRUVpXr16pX56xV3Bsblcqlu3bplnrMkhmEoIyNDDRo0CMiZn0Dlturc7MfwmJv9GB5zsx+DN7cVM7MfgzdvIOdmPwZv3kDObeZ+LE+1OHDggNxut2JjY5WXl1fyQCNAJk+ebEgyJBmzZs0qdkxSUpIhyUhKSvL712/btq3f5zQMw8jMzDQkGZmZmQGZP1C5rTo3+zE85mY/hsfc7MfgzW3FzOzH4M0byLnZj8GbN5BzW3U/etsNfH4KmbfatGnjWT5y5EigvgwAAACACOLze2C8VVBQ4FmuWrVqoL4MAABAUB05dURp89K0/ej2C47b2XenOs7qGJAMgZrb7XZLQ6Se/+hZ7P3L5WXFPxMrZrbqfjyYfdCrcT4VmDvuuEPr1q1T3bp1tXjx4gv+gWzevNmz3KpVK1++DAAAQMgasWiE/rX5X6UPrC6tc60LTIhAzl1P+u3Qb4GZ24p/JlbMLFlzP5Z83/45fCowHTt21EcffaSNGzfq22+/Vd++fYsdd/r0aX3wwQeSpMqVK6tnz56+fBkAAICQ9M2Ob/TRbx9JkmyyKT4mvsSxuTm5ik8oeXt5BGxuQ8rJyVFCQoIUgLdUWPHPxIqZrbofc2w5Xo3zqcDceuutGjt2rNxut+677z7997//VcOGDc8Zc/r0aQ0fPlw7d+6UJI0YMUIVKlTw5csAAACEnNP5pzV0/lDP+t+u/5vu7XRvieNTUlK0adOmgGQJ1NxZWVmqWrWqXJkuJSYm+n1+K/6ZWDGzVfdj8lvJcmY5Sx3nU4Fp2LChXn31VY0aNUpOp1Nt27ZVWlqaOnbsqJiYGG3evFkff/yxtm8/c01oly5d9PTTT5ftd1BOdrvdlK9bXoHMbcW52Y/hMTf7MTzmZj8Gb24rZg60UPgzef7757Xz2Jkf0F7Z+Erd0/Eev8xbFuzH4M1txcyBZnZun98DYxiGHn/8cU2bNu2CL6Hs16+fPvroI9WoUaPEMcnJyXI6nUpKSlJ6erovMUxzttFmZmYGpNEiONiP4YH9GB7Yj+Eh3PfjpkOb1HFWR+W78xUbFat1Q9apbe22Zsfyu3Dfj5HCqvvR227g82MJbDabXn75ZW3atEn333+/unTpolq1aikuLk6NGjXSrbfeqgULFmjhwoUXLC8AAABW4DbcSpuXpnx3viTpicufCMvyAlhFmR+j3Lp1a73zzjv+zAIAABBy3v31Xa3Yu0KS1Lx6cz15xZMmJwIiW8BeZAkAAGB1B7MPasziMZ71t657SxVieTgRYCYKDAAAQAke++YxHcs5Jkm6vcPturb5tSYnAkCBAQAAKMaSnUv0wfoz77WrllBN0/pOMzkRAIkCAwAAcJ6cgpxz3vnyYp8XVbdyXRMTATirzDfxAwAAhLKfM37WjZ/eqIwTGT5/1jAMGTrzponLGl6mBzo/4O94AMqIAuOj+Ph4jR8/XvHx8WZHQTmwH8MD+zE8sB/DQ6jtx9yCXP31339Velb53jMXExWj2amzFWWLjItWQm0/omzCfT/6/CJLf7LiiywBAEDom7x8sp5Z9owkqW6lumpUtZHPc8RFx2lI1yG686I7/R0PQDG87QacgQEAAGFl+9Htmrx8siQp2hatRXcuUsd6Hc0NBcBvIuN8KAAAiAiGYWjo/KHKLcyVJD1y6SOUFyDMUGAAAEDY+HjDx/p257eSpEZVG2nCVRPMDQTA7ygwAAAgLBw7fUyjvh7lWX9zwJuqHFfZxEQAAoECAwAAwsLYb8fqYPZBSdKf2v5Jg1oPMjkRgECgwAAAAMv7Ye8PmvPLHElSlbgqeqP/GyYnAhAoFBgAAGBpeYV5SpuX5lmffPVkJSUmmZgIQCBRYAAAgKVN++80bTy0UZLUpX4X2bvZTU4EIJB4DwwAADBdxokMHTt9zOfPHTl9RBO/myhJirJFac6gOYqOivZ3PAAhhAIDAABMNfvn2Royf0i55xnRfYQ61+/sh0QAQhmXkAEAANPsOb5Ho78ZXe55khOT9Vzv5/yQCECoC4kzMC6XSykpKcVus9vtstu5lhUAgHBjGIYeXviwTuWfkiT1bNRTrWu29nmeCjEVNKTrEFWJr+LviAD8zOFwyOFwFLvN5XJ5NYfNMAzDn6F8kZycLKfTqaSkJKWnp5sVAwAAmODzzZ/rz5/9WZJUv3J9bbZvVtWEqianAmAWb7sBl5ABAICgy8rN0vCFwz3rbwx4g/ICwCsUGAAAEHTPLH1GGScyJEkDWw7Un9v+2eREAKyCAgMAAILqJ+dPmrF6hqQz9684Bjpks9lMTgXAKigwAAAgaArcBUqblyZDZ27BnXDVBDWp1sTcUAAshQIDAACC5s3Vb+rXA79KkjrU6aBRl44yOREAq6HAAACAoNiXuU9PL31akmSTTbNTZys2OtbkVACshgIDAACCYsSiEcrOz5YkpXVJU4+GPUxOBMCKQuJFlgAAwPrWHlirLYe3FLtt57GdmrtlriSpbqW6mtJnShCTAQgnFBgAAFBuC35foNSPUj0351/Ia/1fU7WEaoEPBSAscQkZAAAol+y8bA2dP9Sr8jKgxQDd0u6WIKQCEK44AwMAAMplwn8maG/mXklSj+QeJRaUKvFVdHO7m3nnC4ByocAAAIAyW3tgraavnC5Jio+O1/s3vq8WNVqYnApAOOMSMgAAUCaF7kKlzUtToVEoSXr6yqcpLwACjgIDAADKZPaa2VrtXC1JalurrcZcNsbkRAAiAQUGAAD4LONEhsYtGedZn5U6S/Ex8SYmAhApKDAAAMBno74epazcLEnSfR3v05WNrzQ5EYBIQYEBAAA+Wfj7Qn228TNJUq2KtfTStS+ZnAhAJKHAAAAAr53KP6VhC4Z51l/t+6pqVqxpYiIAkYYCAwAAvPbcd89p9/HdkqTeTXrrrxf91dxAACIO74EBACCCFLoL9c4v72hP5h6fP1vgLvC88yUuOk5vXfcWL6UEEHQUGAAAIsjE7yZq0vJJ5Z7nyZ5PqnWt1n5IBAC+4RIyAAAixOZDmzV1xdRyz3NR3Ys0tudYPyQCAN9xBgYAgAjgNtxKm5emfHe+JGlY12G6KeUmn+eJiYpR5/qdeecLANNQYAAAiADvrX1P3+/9XpLUvHpzvdL3FVWIrWByKgDwHZeQAQAQ5g5lH9KYxWM86zOvm0l5AWBZFBgAAMLcY4sf09HTRyVJt7W/TX2b9zU5EQCUHQUGAIAwtmzXMr2/7n1JUrWEaprWb5rJiQCgfCgwAACEqdyCXA2ZP8SzPvWaqapXuZ6JiQCg/CgwAACEqakrpmrbkW2SpB7JPfRglwdNTgQA5RcSTyFzuVxKSUkpdpvdbpfdbg9yIgAArG3r4a16YcULks48+nh26mxF2fi5JQBzORwOORyOYre5XC6v5rAZhmH4M5QvkpOT5XQ6lZSUpPT0dLNiAAAQVgzD0DXvX6Nlu5dJkp64/AlN7VP+F1gCQCB52w1C4gwMAAA4l+ukS+OWjNOhU4d8/uzJvJP6z+7/SJKaVGuiZ3s96+d0AGAeCgwAACHGMAz99d9/1eKdi8s918yBM1UxtqIfUgFAaOBiWAAAQszHGz72S3kZ1nWYBrQc4IdEABA6OAMDAEAIOXb6mEZ9Pcqz/ulNn+rqplf7PE9cdJwS4xP9GQ0AQgIFBgCAEDL227E6mH1QknRjmxt1c7ubTU4EAKGFS8gAAAgRP+77UXN+mSNJqhxXWW8MeMPkRAAQeigwAACEgPzCfKXNS/OsT+49WcmJySYmAoDQRIEBACAETPvvNG04uEGS1KV+Fz3c/WGTEwFAaKLAAABgsl3HdmnidxMlSVG2KM1Ona3oqGiTUwFAaKLAAABgIsMwNGzBMJ0uOC1JGt59uLo06GJyKgAIXRQYAABM9L+b/leLti+SJCVVSdKk3pNMTgQAoY0CAwCASY7nHNfIRSM96zMGzFCV+ComJgKA0Md7YAAACJATuSd019y79Mv+X4rdfir/lA6fOixJur719RrcZnAQ0wGANVFgAAAIkKeWPqW5W+aWOq5SbCXNGDBDNpst8KEAwOIoMAAABMBPzp/05uo3JUkxUTGqXbF2seMqxlbU81c/r0ZVGwUzHgBYFgUGAAA/K3AXKG1emgwZkqQXrn5BYy4fY3IqAAgP3MQPAICfzVg1Q78e+FWSdFHdi/TIpY+YGwgAwggFBgAAP9qXuU/PLHtGkmSTTbNTZys2OtbkVAAQPigwAAD40fCFw5Wdny1JGtJ1iC5NvtTkRAAQXigwAAD4ydwtc/XF1i8kSfUq19ML17xgciIACD8UGAAA/OBE7gkNXzjcs/5av9dULaGaeYEAIExRYAAA8INnlz2r9Kx0SVL/Fv11c7ubTU4EAOGJAgMAQDn9sv8XvbH6DUlSQkyCHAMdvJQSAAKE98AAACCp0F3oufneF27DrbR5aXIbbknS+F7j1ax6M3/HAwD8fxQYAEDE23hwo655/xq5sl3lmqdd7XZ6tMejfkoFACgOl5ABACJaobtQ9395f7nLiyTe+QIAQcAZGABARJuzZo5WOVdJOvPo44vqXuTzHNG2aN3c7mZd3uhyf8cDAPwfFBgAQMTaf2K/xi4Z61n/5M+fqFeTXiYmAgCUhkvIAAARa9TXo5SVmyVJurfjvZQXALCAkDgD43K5lJKSUuw2u90uu90e5EQAgHC38PeF+nTjp5KkmhVq6qVrXzI5EQCEP4fDIYfDUew2l8u7exFthmEY/gzli+TkZDmdTiUlJSk9Pd2sGACACHMq/5TazWyn3cd3S5Leu+E93d3xbnNDAUCE87YbcAkZACDiTPpukqe8XNXkKt118V3mBgIAeI0CAwCIKBsObtAr/31FkhQXHae3rntLNpvN5FQAAG9RYAAAEcNtuJU2L00F7gJJ0tjLx6pNrTYmpwIA+IICAwCIGO/88o5+3PejJKlljZYad8U4kxMBAHwVEk8hAwDAHzJzMrXh4AYZOv/5NLkFuXri2yc867NSZykhJiGY8QAAfkCBAQCEhUPZh9Rpdic5TzhLHXvnRXfq6qZXByEVAMDfuIQMABAWHv3mUa/KS40KNfRq31eDkAgAEAicgQEAWN7SXUv1wfoPJEnVEqrpgU4PFDsuNjpWt7a/VXUq1QlmPACAH1FgAACWllOQoyHzhnjWX+rzkh7s8qCJiQAAgcQlZAAAS5u6Yqp+P/q7JOnyhpfr/s73m5wIABBIFBgAgGVtObxFU1ZMkSTFRMVodupsRdn4pw0Awhl/ywMALMkwDA2ZN0R5hXmSpDGXjVG7Ou1MTgUACDQKDADAkt5f976+2/OdJKlptaZ6+sqnTU4EAAgGCgwAwHIOnzqsR7951LM+87qZqhhb0cREAIBgocAAACzn8cWP68jpI5KkW9rdov4t+pucCAAQLBQYAIClfLf7O7279l1JUtX4qpreb7rJiQAAwcR7YAAAQfdzxs/6yflTmT77+qrXPctTrpmi+lXq+ysWAMACKDAAgKD6Ye8P6vVeLxUaheWa55KkS5TWNc1PqQAAVsElZACAoMkrzFPavLRyl5eEmATNGTSHd74AQATiDAwAIGhe/fFVbTy0UZLUsV5Hjbp0VJnm6Z7UXW1qtfFnNACARVBgAABBsfPYTj23/DlJUpQtSu8MekddGnQxORUAwGo49w4ACDjDMDRs/jDlFORIkkZ0H0F5AQCUCQUGABBwn238TF/v+FqSlJyYrOd6P2dyIgCAVVFgAAABdTznuEYuGulZnzFghqrEVzExEQDAyigwAICAenLJk3JluyRJ17e+XoPbDDY3EADA0igwAICAWZm+UrN+niVJqhRbSTMGzDA5EQDA6igwAICAyC/MV9q8NBkyJEnP9X5Ojao2MjkVAMDqKDAAgIB4beVrWu9aL+nMO19GXDLC5EQAgHDAe2AAAGXy474f9fFvH6vAXXDeNkOGPlj/gSTJJpvmpM5RTBT/5AAAyi8k/jVxuVxKSUkpdpvdbpfdbg9yIgDAhezN3Ku+H/RVdn52qWPt3ezqltQtCKkAAKHO4XDI4XAUu83lcnk1h80wDMOfoXyRnJwsp9OppKQkpaenmxUDAOADwzB0wyc36KttX5U6tmWNlvrpwZ9UNaFqEJIBAKzM224QEmdgAADWMXfLXE95qVe5nr689UvFRseeNy7KFqVWNVspISYh2BEBAGGMAgMA8NqJ3BMavnC4Z/2N/m9weRgAIKh4ChkAwGvPLHtGzhNOSdKAFgN0U8pNJicCAEQaCgwAwCtrMtZoxuozL6KsEFNBjoEO2Ww2k1MBACINBQYAUKoCd4EemveQ3IZbkjS+13g1rd7U5FQAgEhEgQEAlMqx2qFf9v8iSWpfp71G9xhtciIAQKSiwAAALig9K11PL3vasz4ndU6xTx0DACAYKDAAgAsasXCETuadlCSldUlTj4Y9TE4EAIhkFBgAQIm+3Pql/r3l35KkupXqaso1U0xOBACIdLwHBgAi2JqMNXp62dM6kXui2O2bDm3yLE/vN13VK1QPVjQAAIpFgQGACJWdl60/ffYn7c3cW+rYvs376tb2twYhFQAAF8YlZAAQoSZ+N9Gr8tKyRkvNTp3NO18AACGBMzAAEIHWHVinaf+dJkmKj47XuiHr1KJGi2LHRtmiKC8AgJBBgQGACOM23Eqbl6ZCo1CS9NQVT6l1rdYmpwIAwDtcQgYAEWb2z7O1yrlKktS6Zms9fvnjJicCAMB7FBgAiCD7T+zXuCXjPOuzU2crPibexEQAAPiGAgMAEWTU16OUmZspSbqn4z3q1aSXyYkAAPANBQYAIsSi7Yv06cZPJUk1K9TUy9e+bHIiAAB8R4EBgAhwKv+Uhs0f5ll/te+rqlWxlomJAAAoG78XmPz8fF1yySWy2WyaMGGCv6cHAJTBpO8madfxXZKkq5pcpbsuvsvkRAAAlI3fC8zTTz+t1atX+3taAEAZbTi4Qa/89xVJUlx0nN667i3e6wIAsCy/vgdm8eLFeumll/w5JQCgFNuPbtefPv2Tdh/fXez23MJcFbgLJEljLx+rNrXaBDEdAAD+5bcC43K59Ne//tVf0wEAvOA23Lrvi/v028HfSh3bskZLjbtiXKnjAAAIZX4pMG63W/fcc49cLpfq16+v/fv3+2NaAEAp3lv7nr7f+70kKTE+UcmJycWOq55QXTMGzFBCTEIw4wEA4Hd+KTDTp0/XokWLlJCQoLffflupqan+mBYAcAGHsg9pzOIxnvV//uWfurb5tSYmAgAg8Mp9E/+aNWs0btyZSxKmT5+udu3alTsUAKB0jy1+TEdPH5Uk3d7hdsoLACAilKvAnDhxQrfeeqvy8/M1ePBgpaWl+SsXAOAClu5aqvfXvS9JqpZQTdP6TjM5EQAAwVGuAmO327V9+3YlJSXpnXfe4bGcABAEOQU5Gjp/qGf9pT4vqW7luiYmAgAgeMp8D8wHH3ygDz74QDabTf/4xz9Us2bNMocwDENZWVll/nx8fLzi4+PL/HkAsJKpK6Zq25FtkqTLG16u+zvfb3IiAACk3Nxc5ebmlvnzhmF4Na5MBeb333/XsGHDJElPPvmkrrrqqrJM45GRkaGqVauW+fPjx4/XhAkTypUBAKxg6+GtmrJiiiQpJipGs1JnKcrm93cSAwDgsylTpmjixIkB/zo+F5i8vDzddtttOnnypC699FKNHz++3CEaNGigzZs3l/nznH0BEAkMw9CQ+UOUV5gnSXqsx2NqX6e9yakAADhj3LhxGj16dJk/37ZtW2VkZJQ6zucC8+STT2rNmjVKTEzURx99pNjY2DIFLMpmsykxMbHc8wBAOHt/3fv6z+7/SJKaVmuqZ3o9Y24gAACKKO9tHd7eT+9TgcnIyNCrr74qSbr55pu1detWbd269ZwxLpfLs7x9+3YtWrRIklSnTh117tzZly8HAGEnKzdLx04f8/lz2fnZevSbRz3rM6+bqYqxFf0ZDQAAS7AZ3t4tI2n37t1q2rRpmb7QDTfcoLlz557z35KTk+V0OpWUlKT09PQyzQsAVvHFli900//epAJ3QbnmuaXdLfrkpk/8lAoAgNDgbTfgzk8ACIKjp4/qoXkPlbu8JMYnanq/6X5KBQCA9fh0CVmTJk1KfbxZ0bM0PB0MAM4Y++1YHcw+KEnqUKeDUmqn+DxHfEy8Huz8oOpXqe/veAAAWEaZ3wMDAPDOD3t/0Nu/vC1JqhJXRQvvWKikxCSTUwEAYE1cQgYAAZRXmKe0eWme9clXT6a8AABQDhQYAAigV398VRsPbZQkdanfRfZudpMTAQBgbRQYAAiQncd26rnlz0mSomxRmjNojqKjok1OBQCAtfn9HhhvbvQHgHBnGIaGzR+mnIIcSdKI7iPUuT7vwgIAoLw4AwMAAfDZxs/09Y6vJUnJicl6rvdzJicCACA8UGAAwM+O5xzXyEUjPeszBsxQlfgqJiYCACB8UGAAwM+eXPKkXNkuSdINrW/Q4DaDzQ0EAEAY4T0wAOCjncd2auPBjcVuc2W7NOvnWZKkSrGVNGPAjGBGAwAg7FFgAMAHP2f8rJ5/76ncwtxSx07qPUkNqzYMQioAACIHl5ABgJcK3AV68KsHvSovnet31vBLhgchFQAAkYUzMADgpTdWvaG1B9ZKktrWaqs7L7qz2HEVYyvq9g63KyaKv2IBAPA3/nUFAC/szdyrZ5c9K0myyaZ3b3hXlyRfYnIqAAAiD5eQAUApDMPQwwseVnZ+tiRpaNehlBcAAExCgQGAUszdMldfbftKklSvcj29cM0LJicCACByUWAA4AJO5J7Q8IV/3Iz/ev/XVTWhqomJAACIbBQYALiAZ5Y9I+cJpyRpQIsB+kvKX0xOBABAZKPAAEAJ1mSs0YzVZ15EWSGmghwDHbLZbCanAgAgslFgAKAYhe5Cpc1Lk9twS5LG9xqvptWbmpwKAABQYACgGI6fHFqzf40kqX2d9hrdY7TJiQAAgMR7YABEqI0HN2retnkqNArP22YYhqb+MNWzPjt1tmKjY4MZDwAAlIACAyDiOLOcuvzvlyszN7PUsWld0nRZw8uCkAoAAHiDS8gARJyRi0Z6VV4aJjbUlGumBCERAADwVkicgXG5XEpJSSl2m91ul91uD3IiAOFq3rZ5+tfmf0mSalesrTmD5ijaFn3euChblC5NvlTVK1QPdkQAAMKWw+GQw+EodpvL5fJqDpthGIY/Q/kiOTlZTqdTSUlJSk9PNysGgAiRnZetlJkp2pu5V5L0wY0f6M6L7jQ5FQAAkLzvBlxCBiBiTPxuoqe8XNP0Gt3R4Q6TEwEAAF9RYABEhHUH1mnaf6dJkuKj4zXzupm8lBIAAAuiwAAIe27DrbR5aZ5HJj95xZNqVbOVyakAAEBZUGAAhL3ZP8/WKucqSVLrmq31xOVPmJwIAACUFQUGQFjbf2K/xi0Z51mflTpL8THxJiYCAADlQYEBENZGfT3K886Xezreo6uaXGVuIAAAUC4UGABha9H2Rfp046eSpJoVaurla182OREAACivkHiRJQCUxeZDm/Xqf1/VybyTxW5fvme5Z/mVvq+oVsVawYoGAAAChAIDwJJO55/WoI8HacexHaWO7dW4l+6++O4gpAIAAIHGJWQALOn575/3qrzUqlhLcwbN4Z0vAACECc7AALCcTYc26aUfXpIkxUbFaundS9UwsWGxY+tVrsdTxwAACCMUGACWcvallPnufEnSE5c/oZ6NepqcCgAABAuXkAGwlHd/fVcr9q6QJLWo0UJPXvGkyYkAAEAwUWAAWMbB7IMas3iMZ/2t695ShdgKJiYCAADBRoEBYBmPffOYjuUckyTd0eEO9WnWx+REAAAg2CgwACxh6a6l+mD9B5Kk6gnVNa3fNJMTAQAAM1BgAIS8nIIcDZk3xLP+Yp8XVadSHRMTAQAAs1BgAIS8Kd9P0e9Hf5ckXd7wct3f+X6TEwEAALNQYACEtC2Ht2jqD1MlSTFRMZqdOltRNv7qAgAgUvEeGACm2pu5Vw98+YCcJ5zFbj+UfUh5hXmSpDGXjVG7Ou2CGQ8AAIQYCgwA0xiGoXu/uFdLdy0tdWzTak319JVPByEVAAAIZRQYAKb5x/p/eMpLTFSMKsQU/06XWhVr6R9/+ocqxlYMZjwAABCCKDAATHHk1BGN/ma0Z/2LW7/QwJYDTUwEAACsgDthAZjiiW+f0OFThyVJN6XcRHkBAABeocAACLrv93yvv/36N0lSlbgqer3/6yYnAgAAVkGBARBUeYV5SpuX5ll/4ZoX1KBKAxMTAQAAK6HAAAiqV358RZsPb5YkdWvQTUO7DjU5EQAAsBIKDICg2XF0hyYtnyRJirZFa86gOYqOijY5FQAAsBIKDICgMAxDwxYMU05BjiRp5CUj1bFeR3NDAQAAywmJxyi7XC6lpKQUu81ut8tutwc5EQB/+2TDJ/pmxzeSpIaJDTWx90STEwEAgGBzOBxyOBzFbnO5XF7NYTMMw/BnKF8kJyfL6XQqKSlJ6enpZsUAEGDHTh9TW0dbubLP/MX0xa1f6PrW15ucCgAAhBJvuwGXkAEIuHFLxnnKy41tbqS8AACAMqPAAAio/+77r2avmS1JqhxXWW8MeMPkRAAAwMooMAACJr8wXw/Ne8izPrn3ZCUnJpuYCAAAWB0FBkDATF85XRsObpAkda7fWQ93f9jkRAAAwOooMAACYtexXZrwnwmSpChblOak8s4XAABQfhQYAH5nGIbsC+w6XXBakvRwt4fVpUEXk1MBAIBwQIEB4Hf/3PRPLdy+UJKUVCVJk6+ebHIiAAAQLigwAPwqMydTIxaN8KzPGDBDVeKrmJgIAACEEwoMAL96aulTOnDygCRpUKtBGtxmsLmBAABAWKHAAPCb1c7VmvnTTElSpdhKenPgm7LZbCanAgAA4YQCA8AvCtwFeuirh2TIkCRNvGqiGlVtZHIqAAAQbmLMDgDAOvZl7tOR00eK3fbvzf/WOtc6SdLFdS/WyEtHBjMaAACIEBQYAF6Zs2aO0uallTrOJpvmDJqjmCj+egEAAP7HJWQASrU3c69Gfz3aq7HDug1T96TuAU4EAAAiFT8iBVCqEQtHKDs/W5J0RaMr1K52u2LHJScm67HLHgtmNAAAEGEoMAAuaO6Wufpi6xeSpHqV6+nL275UtYRq5oYCAAARi0vIAJToRO4JDV843LP+Wr/XKC8AAMBUFBgAJXp22bNKz0qXJPVr3k83t7vZ5EQAACDSUWAAFOuX/b/ojdVvSJISYhI087qZvJQSAACYjgID4DyF7kKlzUuT23BLkp698lk1q97M5FQAAAAUGADFmPnTTP2c8bMkqV3tdnr0skdNTgQAAHAGBQbAOZxZTj219CnP+uzU2YqLjjMxEQAAwB8oMADOMXLRSJ3IOyFJerDzg7q80eUmJwIAAPgD74EBIszPGT9r06FNxW7bm7lX/9r8L0lS7Yq1NbXP1GBGAwAAKBUFBoggi7Yv0sAPB8qQUerY6f2mq0aFGkFIBQAA4D0uIQMiRHZetobMG+JVeenXvJ9u73B7EFIBAAD4hjMwQIR47rvntCdzjySpR3IP/fWivxY7rnJcZd3Y9kbe+QIAAEISBQaIAL+5ftO0ldMkSXHRcXpv8HtqVbOVyakAAAB8xyVkQJhzG249NO8hFbgLJElPXfEU5QUAAFhWSJyBcblcSklJKXab3W6X3W4PciIgfLy95m2tTF8pSWpds7WeuPwJkxMBAIBI5XA45HA4it3mcrm8msNmGEbpd/QGSHJyspxOp5KSkpSenm5WDCBsHTh5QG3ebKPM3ExJ0rK7l+mqJleZGwoAAKAY3nYDLiEDwtjor0d7ysvdF99NeQEAAJZHgQHC1Dc7vtHHGz6WJNWsUFOv9H3F5EQAAADlR4EBwtDp/NMaOn+oZ/2Vvq+oVsVaJiYCAADwDwoMEIYmL5+sncd2SpJ6Ne6luy++2+REAAAA/hESTyED4JsVe1foq61fyW24z9tWaBRqxuoZkqTYqFjNSp3FSykBAEDYoMAAFrP50GZd8/41yivMK3Xs2J5j1aZWmyCkAgAACA4uIQMsxDAMDZk/xKvyclHdi/TkFU8GIRUAAEDwcAYGsJD31r6n5XuWS5KaVW+m9254r9jLw6Jt0bq43sVKiEkIdkQAAICAosAAFnEo+5AeW/yYZ/2t697SFY2vMDERAABA8HEJGWARYxaP0dHTRyVJt7W/TX2b9zU5EQAAQPBRYAALWLZrmf5n3f9IkqrGV9W0ftNMTgQAAGAOCgwQ4nILcjVk/hDP+ot9XlS9yvVMTAQAAGAeCgwQ4qaumKptR7ZJknok99CDXR40OREAAIB5KDBACNt2ZJteWPGCJCkmKkazU2crysZhCwAAIhffCQEhyjAMDZ0/1PPOl0d7PKoOdTuYnAoAAMBcFBggRP1j/T+0dNdSSVKTak30bK9nTU4EAABgPgoMEIKOnDqi0d+M9qzPHDhTFWMrmpgIAAAgNFBggBD0xLdP6PCpw5Kkv6T8RQNaDjA5EQAAQGigwAAh5vs93+tvv/5NkpQYn6jX+r9mbiAAAIAQQoEBQkheYZ7S5qV51l+4+gU1qNLAxEQAAAChhQIDhJCXf3hZmw9vliR1T+quIV2HlPIJAACAyEKBAULE9qPbNfn7yZKkaFu0ZqfOVnRUtMmpAAAAQgsFBggBhmFo2PxhyinIkSQ9cukj6livo7mhAAAAQhAFBggBn2z4RIt3LpYkNaraSBOummBuIAAAgBBFgQFMduz0MT3y9SOe9TcHvKnKcZXNCwQAABDCKDCAycZ+O1YHsw9Kkv7U9k8a1HqQyYkAAABCFwUGMNGP+37UnF/mSJIqx1XW6/1fNzkRAABAaKPAACbJL8w/550vz1/9vJITk01MBAAAEPpizA4AhKv8wnzd/+X9WrJriQzDOH+7O1+HTx2WJHWp30X2bvZgRwQAALCckCgwLpdLKSkpxW6z2+2y2/nGDtYzfeV0fbD+g1LHRdmieOcLAACICA6HQw6Ho9htLpfLqzlsRnE/Gg6S5ORkOZ1OJSUlKT093awYgN/tPr5bKY4UnS44LZtsali1YbHj4qLj9Mglj8jenZIOAAAim7fdICTOwADhxDAM2RfYdbrgtCTp4e4P640Bb5icCgAAIDxwEz/gZ//a/C8t+H2BJKlBlQaafPVkkxMBAACEDwoM4EeZOZkasXCEZ/2N/m8oMT7RxEQAAADhhQID+NHTS5/W/pP7JUmprVL1p7Z/MjkRAABAeKHAAH6y2rlajp/OPFWjYmxFvTngTdlsNpNTAQAAhBcKDOAHBe4Cpc1Lk6EzD/WbeNVENa7W2ORUAAAA4YcCA/jBG6ve0NoDayVJF9W9SCMvGWluIAAAgDBFgQHKaW/mXj277FlJkk02zU6drdjoWJNTAQAAhCcKDFBOwxcOV3Z+tiRpSNchujT5UpMTAQAAhC8KDFAOc7fM1Zdbv5Qk1atcTy9c84LJiQAAAMIbBQYooxO5J/Twgoc966/1e03VEqqZFwgAACACUGCAMnp22bNynnBKkvo176eb291sciIAAIDwR4EBymBNxhq9sfoNSVJCTIJmXjeTd74AAAAEAQUG8FGhu1Bp89LkNtySpPG9xqtZ9WYmpwIAAIgMFBjAR46fHFqzf40kqX2d9nq0x6MmJwIAAIgcFBjAB+lZ6Xpq6VOedd75AgAAEFwUGMAHIxeN1Mm8k5Kkhzo/pMsaXmZyIgAAgMhCgQG89NXWr/T55s8lSXUq1dHUPlNNTgQAABB5Yso7wcmTJ7V9+3adOnVKLVu2VK1atXgaE8LOybyTenjhH+98md5vuqpXqG5iIgAAgMhUpjMwOTk5eu6559S4cWNVqVJFnTp10uWXX646deqoRo0aeuqpp5SVleXvrIBpJvxngvZm7pUkXdvsWt3W/jaTEwEAAEQmm2EYhi8fyM7O1tVXX63Vq1dfcFzt2rW1ZMkSdejQocQxycnJcjqdSkpKUnp6ui8xAL/af2K/thzeUuy2g9kHdcfnd6jQKFR8dLw2DNugFjVaBDkhAABAePO2G/h8CdmYMWM85aVTp06aMmWKLrroIiUkJGjDhg2aNGmSFi9erEOHDummm27Shg0bFBvLU5oQun5z/aYef+uh7PzsUsc+c+UzlBcAAAAT+XQGJjMzU3Xq1FFeXp6aN2+uTZs2KS4u7pwxbrdbvXv31vLlyyVJS5Ys0dVXX13sfJyBgdnchluX//1yrUxfWerYlNop+jXtV8VFx5U6FgAAAL4JyBmYzZs3Ky8vT5L02GOPnVdeJCkqKkqjRo3yFJh169aVWGAAs81ZM8dTXppVb6Zb2t1S7LiKsRV1T8d7KC8AAAAm86nA/P77757ldu3alTiuZcuWxX4GCCUHTh7Q2G/Hetb/dv3fdFWTq8wLBAAAgFL5VGBatGihKVOmSJJSUlJKHOd0Oj3L9erVK2M0ILBGfT1KmbmZkqR7Ot5DeQEAALAAnwpMjx491KNHjwuOyc/P19Spf7zg78YbbyxbMiCAvt7+tT7Z8IkkqWaFmnr52pdNTgQAAABvlPtFlhkZGVqxYoWOHj2qbdu26bPPPvOcgZk0adIFH6MMmOF0/mkNWzDMs/5K31dUq2ItExMBAADAW+UuMMuXL9dtt537Ur8aNWroww8/VP/+/b2awzCMcr34Mj4+XvHx8WX+PCLL5OWTtfPYTklSr8a9dPfFd5ucCAAAwPpyc3OVm5tb5s97+3DkcheY4hw9elSPPPKIpk+frgEDBpQ6PiMjQ1WrVi3z1xs/frwmTJhQ5s8jcmw8uFEv/fiSJCk2KlazUmfJZrOZnAoAAMD6pkyZookTJwb86/j0HpgLycrK0q5du/T5559r+vTpOnHihKKjo/XPf/5TgwcPLvYzZ5/13KBBA23evLnMX5szMPCG23Cr13u9tGLvCklnXkr5XO/nTE4FAAAQHsp7BqZt27bKyMgo9T0wfiswRa1YsUK9evWS2+1Wy5YttXnzZkVHR583jhdZIpje+eUdPfjVg5KkFjVa6LehvykhJsHkVAAAAJC87wZRgfjiPXv2VL9+/SSdeQ8M5QRmO5h9UI8vftyzPuu6WZQXAAAAC/KpwIwdO1Z33nmnnnnmmVLHtm3b1rOckZHhezLAjx795lEdyzkmSbrzojt1TbNrTE4EAACAsvCpwKxfv14ffvih3nnnnVLHFn2ZZbVq1XwOBvjLtzu/1T/W/0OSVD2hul7t+6rJiQAAAFBWPhWYVq1aSZIOHDigLVu2lDju5MmT+v777yWdKS8tWrQoR0Sg7HIKcjR0/lDP+kvXvqQ6leqYmAgAAADl4VOBueSSSzzLQ4cO1cmTJ88bk5ubq4cffthz2dhdd92l2NjYcsYEyuaF71/Q9qPbJUk9G/XUfZ3uMzkRAAAAysOnp5Dl5+erU6dO2rhxoyQpKSlJ999/v1q2bKno6Ght27ZNH3zwgXbs2CFJatGihVatWqUaNWoUOx9PIUMgbT60WRfPulj57nzFRMVobdpatavTzuxYAAAAKIa33cCnF1nGxsbqyy+/1HXXXactW7bI6XTqueeKf49Gjx499Mknn5RYXoBAMgxDQ+YPUb47X5L0+GWPU14AAADCgE8FRpKaNWumNWvW6P3339fXX3+tnTt3ateuXYqPj1fr1q3VqlUr3XDDDUpNTS323S9AMLy39j0t37NcktSsejM9feXTJicCAACAP/hcYCSpYsWKGjJkiIYMGeLvPEC5HT51WGMWj/Gsv3XdW6oQW8HERAAAAPCXgLzIEjDTY988piOnj0iSbmt/m/o272tyIgAAAPgLBQZhZdmuZfqfdf8jSaoaX1XT+k0zOREAAAD8iQKDsJFbkKsh8/+4rPHFPi+qXuV6JiYCAACAv1FgEDamrpiqbUe2SZJ6JPfQg10eNDkRAAAA/I0Cg7Cw9fBWvbDiBUlSTFSMZqfOVpSN/70BAADCDd/hwfIMw9DQ+UOVV5gnSRp96Wh1qNvB5FQAAAAIBAoMLO+D9R9o2e5lkqQm1Zro2V7PmpwIAAAAgUKBgaUdOXVEj37zqGfdMdChSnGVTEwEAACAQKLAwNIeX/y4Dp86LEn6S8pfNLDlQJMTAQAAIJAoMLCs5XuW6+9r/y5JSoxP1Gv9XzM3EAAAAAKOAgNLyivM05B5f7zz5YWrX1CDKg1MTAQAAIBgoMDAkl7+4WVtPrxZktStQTcN6TqklE8AAAAgHFBgYDnbj27XpOWTJEnRtmjNGTRH0VHRJqcCAABAMFBgYCmGYWjY/GHKLcyVJI28ZKQ61utobigAAAAEDQUGlvLxho+1eOdiSVLDxIaa2HuiyYkAAAAQTBQYWMax08c06utRnvU3B76pynGVTUwEAACAYKPAwDLGfjtWB7MPSpJubHOjrm99vcmJAAAAEGwUGFjCj/t+1Jxf5kiSKsdV1hsD3jA5EQAAAMwQY3YASXK5XEpJSSl2m91ul91uD3IihJL8wnylzUvzrE/uPVnJickmJgIAAEBZOBwOORyOYre5XC6v5rAZhmH4M5QvkpOT5XQ6lZSUpPT0dLNiIMS9uOJFjV0yVpLUpX4XrXpgFY9NBgAACDPedgMuIUNI23VslyZ+d+ZJY1G2KM1OnU15AQAAiGAUGIQswzA0bMEwnS44LUka3n24ujToYnIqAAAAmIkCg5D1v5v+V4u2L5IkJVVJ0qTek0xOBAAAALNRYBCSMnMyNXLRSM/6jAEzVCW+iomJAAAAEAooMAhJTy19SgdOHpAkDWo1SIPbDDY3EAAAAEICBQYhZ7VztWb+NFOSVCm2kt4c+KZsNpvJqQAAABAKKDAIKQXuAj301UMydObp3hOvmqhGVRuZnAoAAAChIiReZInIUeAu0B2f36Gvtn7lKSlFuQ238grzJEkX171YIy8ded4YAAAARC4KDIJqxqoZ+mzjZ6WOs8mmOYPmKCaK/0UBAADwB747RNDszdyrZ5Y9I+lMQelQt4NsOv/eltjoWD3Q6QF1T+oe7IgAAAAIcRQYBM2IhSOUnZ8tSRrSdYhmXjfT5EQAAACwGm7iR1DM3TJXX2z9QpJUr3I9vXDNCyYnAgAAgBVRYBBwJ3JPaPjC4Z711/q9pmoJ1cwLBAAAAMuiwCDgnl32rNKz0iVJ/Vv0183tbjY5EQAAAKyKAoOA+mX/L3pj9RuSpISYBDkGOngpJQAAAMqMAoOAKXQXKm1emtyGW5I0vtd4NavezORUAAAAsDIKDAJm5k8z9XPGz5KkdrXb6dEej5qcCAAAAFZHgUFAOLOcemrpU5712amzFRsda2IiAAAAhAMKDAJi5KKROpF3QpL0YOcHdXmjy01OBAAAgHBAgYHfzds2T//a/C9JUp1KdTS1z1STEwEAACBcUGDgV9l52bIvsHvWp/WdphoVapiYCAAAAOGEAgO/mvCfCdqbuVeS1KdZH93e4XaTEwEAACCcUGDgN+sOrNP0ldMlSfHR8Xrrurd45wsAAAD8igIDvzj7zpdCo1CS9PSVT6tFjRYmpwIAAEC4ocDAL2avma1VzlWSpLa12mrMZWNMTgQAAIBwRIFBue0/sV/jlozzrM9KnaX4mHgTEwEAACBcUWBQbo98/YiycrMkSfd1vE9XNr7S5EQAAAAIVxQYlMvC3xfqs42fSZJqVayll659yeREAAAACGcUGJTZqfxTGrZgmGf91b6vqmbFmiYmAgAAQLiLMTuAJLlcLqWkpBS7zW63y263F7sN5pr03STtPr5bktS7SW/99aK/mhsIAAAAIc3hcMjhcBS7zeVyeTWHzTAMw5+hfJGcnCyn06mkpCSlp6ebFQNl8JvrN3We01kF7gLFRcdp/ZD1al2rtdmxAAAAYFHedgMuIYPP3IZbQ+YPUYG7QJL0ZM8nKS8AAAAICgoMfPbOL+/ox30/SpJa1WylsT3HmpwIAAAAkYICA5+4Trr0xLdPeNZnXcc7XwAAABA8FBj4ZPQ3o3U857gk6a6L71Lvpr3NDQQAAICIQoGB177Z8Y0++u0jSVKNCjX0yrWvmJwIAAAAkYYCA6+czj+tYfP/eOfLy9e+rNqVapuYCAAAAJGIAgOvPP/989pxbIck6crGV+rejveanAgAAACRiAKDUm06tEkv/fCSJCk2Klazrpslm81mcioAAABEIgoMLshtuJU2L0357nxJ0hOXP6G2tduanAoAAACRigKDC3r313e1Yu8KSVLz6s315BVPmpwIAAAAkYwCgxIdzD6oMYvHeNbfuu4tVYitYGIiAAAARDoKDEr02DeP6VjOMUnS7R1u17XNrzU5EQAAACIdBQbFWrprqT5Y/4EkqVpCNU3rO83kRAAAAAAFBsXIKcjRkHlDPOsv9nlRdSvXNTERAAAAcAYFBueZumKqfj/6uyTpsoaX6YHOD5icCAAAADiDAoNzbDm8RVNWTJEkxUTFaHbqbEXZ+N8EAAAAoYHvTOFhGIaGzBuivMI8SdJjPR5T+zrtTU4FAAAA/IECA4/3172v7/Z8J0lqWq2pnun1jMmJAAAAgHNRYCBJOnzqsB795lHP+szrZqpibEUTEwEAAADno8BAkvT44sd15PQRSdIt7W5R/xb9TU4EAAAAnI8CA323+zu9u/ZdSVJifKKm95tuciIAAACgeBSYCJdbkKsh8/9458vUa6aqfpX6JiYCAAAASkaBiXAv/fCSthzeIkm6JOkSpXVNMzkRAAAAUDIKTAT7/cjvev775yVJ0bZo3vkCAACAkMd3qxHKMAwNnT9UuYW5kqRRl47SxfUuNjkVAAAAcGEUmAj14W8fasmuJZKkRlUbacJVE8wNBAAAAHiBAhOBjp4+qtFfj/asOwY6VCmukomJAAAAAO9QYCLQE4uf0KFThyRJf277Z6W2SjU5EQAAAOAdCkyEWbF3hd759R1JUpW4Knq9/+smJwIAAAC8F2N2AElyuVxKSUkpdpvdbpfdbg9yovCUV5intHl/PCb5+aufV1JikomJAAAAEEkcDoccDkex21wul1dz2AzDMPwZyhfJyclyOp1KSkpSenq6WTEixgvfv6Cnlj4lSeraoKtW3r9S0VHRJqcCAAAAvO8GXEIWIXYc3aFJyydJkqJsUZqdOpvyAgAAAMuhwEQAwzA0bMEw5RTkSJJGXjJSnet3NjkVAAAA4DsKTAT4dOOn+mbHN5Kk5MRkPdf7OZMTAQAAAGVDgQlzx3OO65FFj3jW3xzwpirHVTYvEAAAAFAOFJgw9+SSJ+XKPvNEhxta36Ab2txgciIAAACg7CgwYex4znG988uZd75Uiq2kGQNmmJwIAAAAKB8KTBhb8PsC5bvzJUkPdH5ADas2NDkRAAAAUD4UmDD2xdYvPMs3trnRxCQAAACAf1BgwlRuQa4W/L5AklSzQk1d3uhykxMBAAAA5UeBCVNLdy3VybyTkqTUVqmKiYoxOREAAABQfhSYMDV3y1zP8uA2g03LAQAAAPgTBSYMuQ23vtz2pSSpQkwF9W3e1+REAAAAgH9QYMLQaudqHTh5QJLUt3lfVYytaHIiAAAAwD8oMGGIy8cAAAAQrigwYehsgYmyRSm1Vaq5YQAAAAA/osCEmS2Ht2jrka2SpJ6NeqpWxVomJwIAAAD8hwITZs65fKz1YNNyAAAAAIFAgQkzX2z9wrN8Q5sbTEwCAAAA+B8FJozsP7FfK9NXSpI61OmgZtWbmZwIAAAA8C8KTBj5cuuXnmWePgYAAIBwRIEJI3O3zvUsU2AAAAAQjigwYSIrN0tLdi6RJDVMbKhO9TqZnAgAAADwPwpMmFi0fZHy3fmSzpx9sdlsJicCAAAA/I8CEyaKPj75htY8fQwAAADhiQITBvIK8zT/9/mSpGoJ1XRl4ytNTgQAAAAEBgUmDPxn93+UlZslSUptlarY6FiTEwEAAACBQYEJA0UvHxvcerBpOQAAAIBAo8BYnNtwe97/Eh8dr34t+pmcCAAAAAgcCozF/bjvRzlPOCVJfZr1UeW4yiYnAgAAAAInxuwAkuRyuZSSklLsNrvdLrvdHuRE1lDoLtSor0d51v/U9k8mpgEAAAAuzOFwyOFwFLvN5XJ5NYfNMAzDn6F8kZycLKfTqaSkJKWnp5sVw7JmrJqhEYtGSJJSaqfo17RfFRcdZ3IqAAAAwHfedgMuIbMoZ5ZTTy19yrM+O3U25QUAAABhjwJjUSMXjdSJvBOSpAc6PaCejXqanAgAAAAIPAqMBc3bNk//2vwvSVLtirX14rUvmpwIAAAACA4KjMVk52XLvuCPhxpM6zdNNSrUMDERAAAAEDwUGIuZ+N1E7c3cK0m6puk1uqPDHSYnAgAAAIKHAmMh6w6s07T/TpN05qWVb133lmw2m8mpAAAAgOChwFhEobtQafPSVGgUSpKeuuIptazZ0uRUAAAAQHBRYCxizpo5WuVcJUlqXbO1Hr/8cZMTAQAAAMFHgbGA/Sf2a+ySsZ712amzFR8Tb2IiAAAAwBwUGAsY/5/xysrNkiTd2/Fe9WrSy+REAAAAgDkoMCEuvzBf/9z0T0lS5bjKeunal0xOBAAAAJiHAhPiVuxdoWM5xyRJ17W8TrUq1jI5EQAAAGAeCkyIm7tlrmf5htY3mBcEAAAACAEUmBBmGIbmbp0rSYqNitXAlgPNDQQAAACYjAITwta51mlv5l5JUu+mvVU1oarJiQAAAABzUWBCWNHLxwa3HmxaDgAAACBUUGBCWNECc33r680LAgAAAIQICkyI2nVsl9a51kmSujXopqTEJJMTAQAAAOajwISoL7d+6Vke3GaweUEAAACAEBJT1g+6XC7NmjVLP//8s7Zt26Z9+/apadOmateunfr27at7771X0dHR/swaUc4+fUyiwAAAAABnlanALFy4ULfddpsyMzPP+e+bNm3Spk2b9L//+79666239Pbbb6tz585+CRpJjpw6ouV7lkuSWtZoqba12pqcCAAAAAgNPheYTZs26c9//rNOnz4tSRo0aJD69Omj+vXra+vWrfr444+1adMm/fLLL7r66qu1atUqtW7d2u/Bw9m8bfPkNtySzry80mazmZwIAAAACA0+F5iJEyd6ysubb76pYcOGnfMN9uOPP65Ro0Zp5syZyszM1MiRI7Vo0SL/JY4AX2z9wrPM5WMAAADAH3y6iT8nJ0f//ve/JUmdO3c+r7xIUlxcnKZPn6727dtLkr7++ms5nU4/xQ1/p/JPadH2M4WvTqU6ujT5UpMTAQAAAKHDpwKzbt065efnS5JuuummEi9tiouLU2pqqmf9119/LUfEyPLtzm91uuDMGa7rW12v6CgehAAAAACc5VOBcblcnuXGjRtfcGz9+vU9y2cvOUPpir68ksvHAAAAgHP5dA9M+/bt9e6770qSevXqdcGxP/30k2e5VatWZYgWeQrdhfpq21eSpEqxlXRNs2tMTgQAAACEFp8KTLNmzdSsWbNSx61fv16ffPKJJCk5OVlt2/IYYG/8uO9HHT51WJLUv0V/JcQkmJwIAAAACC1lfpFlSdasWaNBgwapoKBAkjRu3DjFxcVd8DOGYSgrK6vMXzM+Pl7x8fFl/nyo4PIxAAAAWFVubq5yc3PL/HnDMLwa57cCk5WVpSlTpuiVV17xlJfbb79dQ4YMKfWzGRkZqlq1apm/9vjx4zVhwoQyfz4UGIahuVvnSpKibdG6ruV15gYCAAAAfDBlyhRNnDgx4F+n3AUmPz9fb7/9tiZMmKBDhw55/vsDDzygmTNnKiqq9OcENGjQQJs3by5zhnA4+7Lx0EbtPLZTktSrSS9Vr1Dd5EQAAACA98aNG6fRo0eX+fNt27ZVRkZGqePKVWDWrl2r++6775zHJDdq1EizZs3SgAEDvJ7HZrMpMTGxPFEs75zLx1oPNi0HAAAAUBblva2jpFe0/F8+PUb5rMLCQj3zzDPq1q2bp7zUqFFDr776qrZu3epTecEZRQvMDW1uMC8IAAAAEMJ8PgNjGIbsdrtmz54t6UxTGj58uCZOnKhq1ar5O19E+M/u/2jN/jWSpM71O6tR1UYmJwIAAABCk88FZvLkyZ7yUrduXX3++ee67LLL/B4sUuQW5GrIvD8edDCkS+kPPQAAAAAilU8FJisrSy+++KIkqXr16vrhhx/UvHnzgASLFC/+8KK2HtkqSeqR3EP3d77f5EQAAABA6PLpHpgPP/xQ2dnZkqTnn3+e8lJO245s0wvfvyDpzKOTZ6XOUpStTLclAQAAABHBpzMwy5cv9yxXqlRJixYt8upz3bt3V40aNXxLFuYMw9DQ+UOVW3jmZT+P9nhUF9W9yORUAAAAQGjzqcC4XC7P8t133+3155YtW6arrrrKly8V9j787UMt3bVUktS4amM92+tZkxMBAAAAoc+n65WKFhiU3ZFTRzTq61Ge9ZnXzVSluEomJgIAAACswaczMBs3bgxUjojyxLdP6PCpw5Kkm1Ju0sCWA01OBAAAAFgDd4wH2fd7vtfffv2bJKlKXBW93v91kxMBAAAA1kGBCaK8wjylzUvzrL9wzQtqUKWBiYkAAAAAa6HABNErP76izYc3S5K6NeimoV2HmpwIAAAAsBYKTJDsOrZLk5ZPkiRF2aI0O3W2oqOiTU4FAAAAWAsFJkg+WP+BcgpyJEkjLxmpTvU7mZwIAAAAsB4KTJCsdq72LA/rNszEJAAAAIB1UWCCwDAM/ZTxkySpekJ1Na/e3OREAAAAgDVRYIJgX9Y+Hcw+KEnq2qCrbDabyYkAAAAAa6LABMFPzp88y92TupuYBAAAALA2CkwQFL3/pVuDbiYmAQAAAKyNAhMEZ+9/kaRuSRQYAAAAoKwoMAHmNtxas3+NJKlBlQZqUKWByYkAAAAA66LABNi2I9uUlZslicvHAAAAgPKiwAQYN/ADAAAA/kOBCbBz7n/hDAwAAABQLhSYACv6BLKuDbqamAQAAACwPgpMAOUV5mntgbWSpBY1Wqh6hermBgIAAAAsjgITQBsOblBuYa4kLh8DAAAA/CHG7ACS5HK5lJKSUuw2u90uu90e5ET+wQ38AAAAwB8cDoccDkex21wul1dzhESBqVu3rjZt2mR2DL/jBn4AAADgDxc6OZGcnCyn01nqHFxCFkBnb+CPtkWrU/1OJqcBAAAArI8CEyDZednaeGijJKldnXaqGFvR5EQAAACA9VFgAuTXA7/KbbglcfkYAAAA4C8UmADhBn4AAADA/ygwAcIN/AAAAID/UWAC5GyBSYhJUPs67U1OAwAAAIQHCkwAHD19VNuPbpckdazXUbHRsSYnAgAAAMIDBSYAfs742bPM5WMAAACA/1BgAoAb+AEAAIDAoMAEADfwAwAAAIFBgQmAswUmMT5RLWu2NDkNAAAAED4oMH7mzHIq40SGJKlrg66KsvFHDAAAAPgL3137GZePAQAAAIFDgfGzojfwU2AAAAAA/6LA+FnRMzA8gQwAAADwLwqMHxmG4XkHTN1KdZWcmGxyIgAAACC8UGD8aMexHTqWc0yS1C2pm2w2m8mJAAAAgPBCgfEj7n8BAAAAAosC40ernas9yxQYAAAAwP8oMH60Zv8az3K3JAoMAAAA4G8UGD/aeWynpDM38NeqWMvkNAAAAED4ocD4SaG7UAdOHpAkJSUmmZwGAAAACE8UGD85mH1QhUahJCmpCgUGAAAACAQKjJ9knMjwLDeo0sDEJAAAAED4osD4ifOE07PMGRgAAAAgMCgwfuLMKlJguAcGAAAACAgKjJ9wCRkAAAAQeBQYP+ESMgAAACDwKDB+whkYAAAAIPBizA4gSS6XSykpKcVus9vtstvtQU7ku7NnYOKj41WjQg2T0wAAAAChx+FwyOFwFLvN5XJ5NYfNMAzDn6F8kZycLKfTqaSkJKWnp5sVwy9qvFhDx3KOqWm1pto5cqfZcQAAAABL8bYbcAmZH5zOP61jOcck8QQyAAAAIJAoMH5Q9P4XbuAHAAAAAocC4wfcwA8AAAAEBwXGD3iEMgAAABAcFBg/cGb9UWA4AwMAAAAEDgXGD865B4ab+AEAAICAocD4AZeQAQAAAMFBgfEDbuIHAAAAgoMC4wdnz8BUT6iuCrEVTE4DAAAAhC8KTDkZhuG5iZ+zLwAAAEBgUWDK6VjOMeUW5kriBn4AAAAg0Cgw5VT0EcrcwA8AAAAEFgWmnLiBHwAAAAgeCkw58QhlAAAAIHgoMOXEGRgAAAAgeCgw5XTOPTDcxA8AAAAEFAWmnLiEDAAAAAgeCkw5nb2ELNoWrTqV6picBgAAAAhvFJhyOnsGpl7leoqOijY5DQAAABDeKDDlUOAukOukSxI38AMAAADBQIEphwMnD8iQIYkb+AEAAIBgoMCUQ9EnkDWozBkYAAAAINAoMOVQ9B0wnIEBAAAAAo8CUw48QhkAAAAILgpMORQ9A8NN/AAAAEDgUWDK4ZwzMFxCBgAAAAQcBaYczrmJnzMwAAAAQMBRYMrh7CVkFWMrqmp8VZPTAAAAAOGPAlMOZy8hS6qSJJvNZnIaAAAAIPzFmB1Aklwul1JSUordZrfbZbfbg5yodCfzTiorN0sSl48BAAAA3nA4HHI4HMVuc7lcXs0REgWmbt262rRpk9kxfMI7YAAAAADfXOjkRHJyspxOZ7HbiuISsjI65wb+ypyBAQAAAIKBAlNGnIEBAAAAgo8CU0bnvAOmCgUGAAAACAYKTBkVPQPDTfwAAABAcFBgyuicMzBcQgYAAAAEBQWmjIqegalfub6JSQAAAIDIQYEpo7NPIatVsZbiY+JNTgMAAABEBgpMGbgNt+cMDDfwAwAAAMFDgSmDI6eOKN+dL4kb+AEAAIBgosCUAY9QBgAAAMxBgSkDHqEMAAAAmIMCUwZnb+CXeIQyAAAAEEwUmDIoegkZZ2AAAACA4KHAlEHRS8i4BwYAAAAIHgpMGZxzEz+XkAEAAABBQ4Epg7NnYGKjYlWrYi2T0wAAAACRgwJTBmdv4q9fpb6ibPwRAgAAAMHCd98+yi3I1aFThyRxAz8AAAAQbBQYHx04ecCzzA38AAAAQHBRYHx0zg38FBgAAAAgqCgwPir6CGUuIQMAAACCiwLjo7M38Es8QhkAAAAINgqMj4peQsYZGAAAACC4KDA+KnoJGffAAAAAAMFFgfERZ2AAAAAA81BgfHT2DEyVuCqqEl/F5DQAAABAZKHA+MAwDM9N/NzADwAAAAQfBcYHWblZys7PlsTlYwAAAIAZKDA+4AZ+AAAAwFwxZgeQJJfLpZSUlGK32e122e32ICcqXtEb+CkwAAAAgG8cDoccDkex21wul1dz2AzDMPwZyhfJyclyOp1KSkpSenq6WTG8dir/lLYe3qqMExlqXK2x2tdpb3YkAAAAICx42w1C4gyMVVSMrahO9TupU/1OZkcBAAAAIhL3wAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMugwAAAAACwDAoMAAAAAMvwW4FxOByy2WyaMGGCv6YEAAAAgHP4rcD84x//8NdUAAAAAFAsvxSYd999VytXrvTHVAAAAABQopiyfjAzM1Pr16/Xu+++y9kXAAAAAEFRpgLTvXt3/fTTT/7OAgAAAAAXVKZLyA4ePOjvHAAAAABQqjIVmK1bt+r06dOeX1u2bPF3LgAAAAA4T5kuIYuPj7/gOgAAAAAEAi+yBAAAAGAZZX4KmT8ZhqGsrKwyfz4+Pp6zQAAAAICJcnNzlZubW+bPG4bh1biQKDAZGRmqWrVqmT8/fvx4TZgwwX+BAAAAAPhkypQpmjhxYsC/TkgUmAYNGmjz5s1l/jxnXwAAAABzjRs3TqNHjy7z59u2bauMjIxSx4VEgbHZbEpMTDQ7BgAAAIAyKu9tHTabzatx3MQPAAAAwDIoMAAAAAAsgwIDAAAAwDIoMAAAAAAsgwIDAAAAwDIoMAAAAAAsgwIDAAAAwDIoMAAAAAAswy8vsmzSpIkMw/DHVAAAAABQIs7AAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy4gxO4AkuVwupaSkFLvNbrfLbrcHOREAAAAAf3M4HHI4HMVuc7lcXs1hMwzD8GcoXyQnJ8vpdCopKUnp6elmxQAAAABgMm+7AZeQAQAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALCMGLMDSJLL5VJKSkqx2+x2u+x2e5ATAQAAAPA3h8Mhh8NR7DaXy+XVHDbDMAx/hvJFcnKynE6nkpKSlJ6eblYMAAAAACbzthtwCRkAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy4gxO4AkuVwupaSkFLvNbrfLbrcHOREAAAAAf3M4HHI4HMVuc7lcXs1hMwzD8GcoXyQnJ8vpdCopKUnp6elmxQAAAABgMm+7AZeQAQAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALCMGLMDSJLL5VJKSkqx2+x2u+x2e5ATAQAAAPA3h8Mhh8NR7DaXy+XVHDbDMAx/hvJFcnKynE6nkpKSlJ6eblYMAAAAACbzthtwCRkAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCoyPcnNzNWHCBOXm5podBeXAfgwP7MfwwH4MD+zH8MB+DA/hvh9thmEYZn3x5ORkOZ1OJSUlKT093awYPsnKylLVqlWVmZmpxMREs+OgjNiP4YH9GB7Yj+GB/Rge2I/hwar70dtuwBkYAAAAAJZBgQEAAABgGRQYAAAAAJZBgQEAAABgGTHlneDYsWPauXOnTpw4oQYNGqhFixaKiqIXAQAAAPC/MjeNbdu2adCgQapdu7a6du2q3r17q3Xr1mrSpIleffVVFRYW+jMnAAAAAJStwHz//ffq1KmT5s2bd15R2bdvnx577DH96U9/osQAAAAA8CufC8zhw4d144036tSpU4qKitJzzz2nffv26eTJk1q6dKk6deokSfryyy/13HPP+T2wtxwOh2lfuzwCmduKc7Mfw2Nu9mN4zM1+DN7cVswcaFb8M7Fi5kCz4p+JFTMHmum5DR+NGTPGkGRIMt58883ztmdmZhpNmjQxJBmVKlUyDh48WOJcSUlJhiQjKSnJ1xilatu2rd/nNIwzvz9JRmZmZkDmD1Ruq87NfgyPudmP4TE3+zF4c1sxM/sxePMGcm72Y/DmDeTcVt2P3nYDn87AFBYW6u9//7skqU6dOhoyZMh5YxITE/XYY49JkrKzs/Xpp5+WsVoBAAAAwLl8KjArV67UkSNHJEmDBg1SdHR0seOuv/56z/L8+fPLEQ8AAAAA/uBTgdm6datneeDAgSWOa9iwoS666CJJ0i+//FLGaAAAAABwLp8KzP79+z3LjRs3vuDYhg0bSpIOHjyo48eP+54MAAAAAP4PnwrMgQMHPMs1atS44NiaNWt6losWHwAAAAAoqxhfBhctIkULSnGKbs/Ozi52zMGDBz3zJiUl+RLlHDab7bz/5nK5lJycXOY5S2IYhiSpbdu2xX7d8gpUbqvOzX4Mj7nZj+ExN/sxeHNbMTP7MXjzBnJu9mPw5g3k3Gbux7NfuyzOniw52xFKYjN8+Cp9+/bV4sWLJUl5eXmKjY0tcey4ceM0depUSdLy5ct1xRVXnDcmOjpabrfb2y8PAAAAIMxFRUWpsLCwxO0+nYGpXbu2Z/n48ePnrP9fRe97qVChQrFjEhISlJOTo+jo6AvOVZpANEsAAAAAvinPGZhDhw6psLBQCQkJFxznU4GpX7++Z/no0aMXLB1Hjx71LFeuXLnYMSVdWgYAAAAAxfHpJv569ep5losWlOIcO3bMs1ye+1sAAAAA4CyfCkzRMzDr1q0rcZzb7dZvv/0mSWrUqJGqVKlSxngAAAAA8AefCkzXrl09y1999VWJ49asWeN5isBll11WxmgAAAAAcC6fCkzr1q3VunVrSdKSJUvOuUysqM8//9yzfOONN5YjHgAAAAD8wacCI0mjR4+WJOXm5mr48OHnPQb5119/1WuvvSZJatq0qQYPHlzukAAAAAAg+fgeGEnKz89Xz549tXr1aknSlVdeqXvuuUeJiYlavXq1Zs6cqZMnT8pms+mLL77QoEGDAhLcDMeOHdPOnTt14sQJNWjQQC1atFBUlM8dEEAQOJ1O7dmzR5LUuHFjHiYCmIRjEfCfkydPavv27Tp16pRatmypWrVq+fQ6kbA5Ho0y2L9/v9GpUydDUrG/YmNjjbfeeqssU4ekrVu3GqmpqUZ0dPQ5v8+GDRsar7zyilFQUGB2RJSiT58+Jf7/+n9/vfvuu2bHRRFvvvmmIckYP368V+MXLFhgdO3a9bz92qVLF2P+/PmBDYtiebsPCwoKjLi4OK+P1V27dgUlf6Q6cOCAMWHCBCM1NdVo1aqVUaFCBSMlJcX4y1/+Yrz99tul/tvHsRgayrofOR5Dx+nTp42JEycajRo1Ou/PvVq1asaTTz5pZGZmXnCOcDsey1RgDMMwcnNzjRkzZhg9evQwatSoYcTFxRlNmjQxHnjgAWP9+vX+zGiq5cuXGxUrVrzgQXv99ddTYkJccQc9BcYaLr30Uq8LzCuvvFLq/n355ZcDHxrn8HYf7ty50+vjlG+YAmvBggVG1apVL/jn37lzZ2PNmjXFfp5jMTSUZz9yPIaGkydPGt27dy/1z7927dolfv8djsejz5eQRZLDhw+rTZs2OnLkiKKiojRhwgTde++9ql69ulavXq1HH31Uv/76qyTp2Wef1cSJE01OjOLk5OSoYsWKMgxDw4YNK/Wyxg4dOlj3lGqYeffdd3XfffdJksaPH68JEyaUOHbp0qXq06ePDMNQrVq19Prrr2vAgAGSpAULFmjkyJE6cuSIJOnbb7/VNddcE/D88G0ffvPNN+rXr58kadasWWrcuPEF5+7Vq5cqVKjgt6w4Y9OmTeratatOnz4tSRo0aJD69Omj+vXra+vWrfr444+1adMmSVLVqlW1atUqzwN+JI7FUFHe/cjxGBqGDRumt956S5LUqVMnTZkyRRdddJESEhK0YcMGTZo0SYsXL5YktWrVShs2bFBsbKzn82F7PJrZnkLdmDFjPM30zTffPG97Zmam0aRJE0OSUalSJePgwYMmpERpNmzY4NmPCxYsMDsOSnH8+HFj+fLlxr333mvExsZ69t2Ffnrvdrs9P6GKiYkxfvnll/PGrFmzxoiJiTEkGd27dzfcbncAfxeRrSz70DD+uNRMknHy5MnghMV5br755nP+7fu/x0pubq4xbNgwz5h+/fp5tnEsho7y7EfD4HgMBcePH/dcxte8eXMjNzf3vDGFhYXGlVde6dlXS5Ys8WwL5+ORAlOCgoICo2bNmoYko06dOiVeIlb0AJ8xY0aQU8Ib//73vz37aPv27WbHwQV069atxNPbF/rmd+3atZ5xN910U4njbrrpJs+43377LQC/A5R1HxqGYYwcOdKQZCQnJwcnLM5z+vRpT+ns3Llzid/M5ObmGu3bt/fs2/T0dMMwOBZDRXn3o2FwPIaC//73v559c6F7y4t+nzNt2jTPfw/n45FHaJVg5cqVnlNqgwYNUnR0dLHjrr/+es/y/Pnzg5INvvn9998lSbGxsaWeAoe5Dh48WKbPzZs3z7N8ww03lDiu6DaO18Ao6z6U/jhWW7Vq5a848NG6deuUn58vSbrppptKfLpRXFycUlNTPetnL6fmWAwN5d2PEsdjKDi7DySpXbt2JY5r2bJlsZ8J5+MxxuwAoWrr1q2e5YEDB5Y4rmHDhrrooou0fv16/fLLL8GIBh+dPZibN2+umJgYFRYWat++fdq9e7eqVKmili1bKjEx0eSUkM4cd0aR2/L27NmjNm3aePW5sy50vBbdxvEaGGXdh9L53zDl5ORo165dOnDggGrVqqU2bdqcc203/M/lcnmWS/uBT/369T3LZ++z4FgMDeXdjxLHYyho0aKFpkyZIklKSUkpcZzT6fQs16tXz7MczscjBaYE+/fv9yyXdvA3bNhQ69ev18GDB3X8+HFVq1YtwOngi23btkk685f0lClT9Morr+jo0aOe7TExMerTp4/Gjx+vSy+91KyYkBQfH3/B9ZKcPV4rV66sGjVqlDiuRo0aqlSpkrKzs8/5ix3+U9Z9mJ+fr507d3o+k5aWpr///e8qKCjwjKlSpYpuueUWPfvss2rYsKH/QsOjffv2evfddyWduSn7Qn766SfP8tlvcjkWQ0N59yPHY2jo0aOHevToccEx+fn5mjp1qmf9xhtv9CyH8/FIgSnBgQMHPMsX2umSVLNmTc/y/v37KTAh5uxPkZYtW6Zly5adt72goECLFi3SN998o9dee03Dhw8PdkSU09njtbRjVTpzvGZnZ5/zQwqYb/fu3SosLJQkvf7668WOOXHihN555x19/vnnmjt3rq644opgRowIzZo1U7NmzUodt379en3yySeSpOTkZLVt21YSx2KoKO9+5HgMXRkZGVqxYoWOHj2qbdu26bPPPvOcgZk0aZI6dOjgGRvOxyP3wJSg6A4sWlCKU3R7dnZ2wDLBd9nZ2crIyPCsd+rUSYsXL9bx48d1+PBhLVu2TLfccoskye12a8SIEVqwYIFZcVFGZ4/X0o7VomM4VkNL0eu2JenWW2/Vr7/+quzsbO3evVtz585Vp06dJElHjx7V4MGDlZ6ebkbUiLdmzRr179/f89P4cePGKS4uThLHopVcaD9yPIau5cuX65ZbbtHQoUM1ffp0OZ1O1ahRQwsXLtTTTz99zthwPh4pMCXIysryLJf2XPOil0gUvX4U5tuxY4eio6MVHR2tQYMGafny5erTp4+qVq2qmjVr6qqrrtInn3yiadOmeT7z6KOPyu12m5gavjp7vHrzDoKzxyvHamgpeqxOnDhRH330kTp27KiKFSuqcePGuuGGG7Rq1Sr9+c9/lnTmm6ZJkyaZnDqyZGVlady4cbr00ks93xjdfvvtGjJkyDljJI7FUObNfuR4tJajR4/qkUce0cKFC8/572F9PJr9GLRQdfvtt3seKVfa+12GDBniGfvTTz8FKSH8qaCgwOjSpYtnP5b0dmkE165du7x6BG+DBg0MSUZKSkqpc7Zp08bz7iYEnrf70Fvp6elGpUqVDElGtWrVjPz8/PKHxAXl5eUZDofDqF279jmPxX7ggQeMvLy8c8ZyLIYuX/ajtzgezZOZmWmsXbvWePbZZ40qVaoYkozo6Gjj3//+t2dMOB+PnIEpQdGnchS94bs4RbdXrlw5YJkQONHR0efc+Fb0pkaEvrPHa2nHatExHKvWlJSU5Lmp9fjx49qxY4fJicLb2rVrdckll8hut+vQoUOSpEaNGmnBggV6++23z3sKFcdiaPJ1P3qL49E8iYmJuvjiizVx4kQtWLBAUVFRKiws1OOPP+65fymcj0cKTAmKPoautB1/7Ngxz3JSUlLAMiGwij7q9ew7gGANZ4/Xo0ePnvMI3//LMAzP8cpTc6yLYzXwCgsL9cwzz6hbt26ed4PUqFFDr776qrZu3aoBAwYU+zmOxdBS1v3oC45H8/Xs2VP9+vWTdOb+pbP3I4Xz8UiBKUHRMzDr1q0rcZzb7dZvv/0m6cxPM6pUqRLwbAiMoo+HrFq1qolJ4Kuzx2teXt4FHwG5detWz8vdLvRMfYQ2jtXAMgxDdrtdkydPVkFBgWw2m0aMGKEdO3Zo9OjRSkhIKPGzHIuhozz70Rccj4EzduxY3XnnnXrmmWdKHXv2CXKSPA8vCufjkQJTgq5du3qWv/rqqxLHrVmzxvOYussuuyzgueCbO+64Q+3bt9c111xT6o35mzdv9izz5mFr8fZ4LbqN4zV0ZGdnq2PHjmrfvr0ef/zxUsefPVZjY2NLfU8XfDd58mTNnj1bklS3bl2tWLFCr7/+ulevCOBYDB1l3Y8cj6Fj/fr1+vDDD/XOO++UOrboyyzP7uOwPh7Nu/0m9LVu3dqQZMTHxxtHjx4tdszYsWM9N8J9+umnQU6I0rz00kue/fP111+XOO7UqVNGs2bNDElG5cqVjVOnTgUxJUri7Q3gTqfTsNlshiTjkksuMdxu93lj3G63cemllxqSDJvNZmRkZAQwOc7ydh8OHDjQcwNpSX/fGoZhbNiwwYiJiTEkGX369AlA4siWmZnpuSm7evXqxvbt2336PMdiaCjvfuR4DA0jR470/P25efPmEsedOHHCc8N+tWrVPA9lCOfjkTMwFzB69GhJUm5uroYPH37eT/B//fVXvfbaa5Kkpk2bavDgwUFOiNLceuutioo687/5fffdp3379p035vTp0xo+fLjnrcMjRozw6pGDCB0NGjTQrbfeKklatWqV5syZc96Y2bNna+XKlZLOnJkrepkozHfnnXdKOvPT31tvvVW5ubnnjXE6nbr33ns9l6x489Nh+ObDDz/0vAfi+eefV/PmzX36PMdiaCjvfuR4DA2XXHKJZ3no0KE6efLkeWNyc3P18MMPey4bu+uuuzwPZQjn49FmGBe4qyfC5efnq2fPnlq9erUk6corr9Q999yjxMRErV69WjNnztTJkydls9n0xRdfaNCgQSYnRnFee+01jRo1SpJUqVIlpaWlqWPHjoqJidHmzZv18ccfa/v27ZKkLl266Pvvv6fAhIjdu3eradOmkqTx48drwoQJFxzbpUsXz0M37r77bg0YMEBut1vz58/Xhx9+KOnMy7rWrFnDpQ5B4u0+dLvd+vOf/6y5c+dKOnNP4UMPPaRWrVopKytLGzZs0N/+9jedOHFCkjRs2DA5HI5g/BYiym233eZ5M/v//M//qE6dOl59rnv37p63fXMsmq+8+5HjMTTk5+erU6dO2rhxo6QzD4q6//771bJlS0VHR2vbtm364IMPPE9/a9GihVatWuU5FqUwPh7NPgUU6vbv32906tTpnGemF/0VGxtrvPXWW2bHxAW43W7jscceM6Kiokrcj5KMfv36GUeOHDE7Lorw9R0i33///XnvOCj6q06dOsaPP/4Y+ODw8GUfnjhxwrj66qsveJzabDZjzJgxRkFBQXB+AxGmd+/eF/zzL+nXsmXLzpmHY9Fc/tiPHI+hYceOHZ53tFzoV48ePYw9e/YUO0c4Ho8xPvadiFOvXj2tXLlSc+bM0UcffaStW7fq5MmTatCggfr06aMRI0aoQ4cOZsfEBdhsNr388st64IEH9PLLL2vt2rXas2ePsrKyVK9ePV122WW666671L9/f9lsNrPjohx69uyp3377TW+88Ybmzp2rPXv2yGazqXHjxho8eLBGjhyp2rVrmx0TJahcubK+/fZbLVmyRA6HQ7///rv27NkjwzDUsGFDXX311Ro2bJjatWtndtSw5XK5/DIPx6K5/LEfOR5DQ7NmzbRmzRq9//77+vrrr7Vz507t2rVL8fHxat26tVq1aqUbbrhBqampio6OLnaOcDweuYQMAAAAgGVwEz8AAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy6DAAAAAALAMCgwAAAAAy/h/VSAQcbwQKdIAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "y_2 = ox_diffusion(time_array, D=0.004, z_idx=z_index_mean)\n", + "plt.plot(time_array, y_2, label=\"full equation\", color=\"green\", lw=2)\n", + "plt.xlim(-2, 32)\n", + "plt.ylim(0, 5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/pre_irr_index.py b/pre_irr_index.py index 9c399fa..e4ab43d 100644 --- a/pre_irr_index.py +++ b/pre_irr_index.py @@ -19,7 +19,7 @@ def sellmeier_eq(wl, B1, C1): wavelengths = np.array([470, 527, 635]) wl_space = np.linspace(400, 700, 1000) -fig, ax = plt.subplots(figsize=(7, 6)) +fig, ax = plt.subplots(figsize=(8, 6)) for material in input.keys(): index = np.array([]) @@ -62,7 +62,7 @@ def sellmeier_eq(wl, B1, C1): print(f"\nFit results for {material}:") print(f" chi2 / ndof = {m.fval:.1f} / {m.ndof:.0f} = {m.fmin.reduced_chi2:.1f}") for p, v, e in zip(m.parameters, m.values, m.errors): - print(f"\t{p} = {v:.3f} ± {e:.3f}") + print(f"\t{p} = {v:.4f} ± {e:.4f}") print() ax.set_xlabel(r"Wavelength $\lambda$ (nm)")