diff --git a/.gitignore b/.gitignore index 6947f9fa..59b29771 100644 --- a/.gitignore +++ b/.gitignore @@ -16,6 +16,7 @@ TidalPy-Run* !TidalPy/RadialSolver/rs_solution_.hpp !TidalPy/Material/eos/eos_solution_.cpp !TidalPy/Material/eos/eos_solution_.hpp +!TidalPy/utilities/dimensions/nondimensional_.hpp # End TidalPy diff --git a/Benchmarks & Performance/RadialSolver/RadialSolver Benchmarks.ipynb b/Benchmarks & Performance/RadialSolver/RadialSolver Benchmarks.ipynb index ad83275a..a9cff40b 100644 --- a/Benchmarks & Performance/RadialSolver/RadialSolver Benchmarks.ipynb +++ b/Benchmarks & Performance/RadialSolver/RadialSolver Benchmarks.ipynb @@ -10,7 +10,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "0.6.0a7.dev8\n" + "0.6.0a7.dev12\n" ] } ], @@ -37,21 +37,41 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 7, "id": "7d9c6658-deca-477c-8c9f-b32d0c7f0913", "metadata": {}, "outputs": [ { - "ename": "TypeError", - "evalue": "radial_solver() takes at least 10 positional arguments (7 given)", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mTypeError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[1;32mIn[4], line 89\u001b[0m\n\u001b[0;32m 78\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m solution\n\u001b[0;32m 80\u001b[0m \u001b[38;5;66;03m# 0.5.4\u001b[39;00m\n\u001b[0;32m 81\u001b[0m \u001b[38;5;66;03m# New: 3.12ms; 3.07ms; 3.15ms\u001b[39;00m\n\u001b[0;32m 82\u001b[0m \u001b[38;5;66;03m# Old: 94ms; 96.6ms; 94.3ms\u001b[39;00m\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 86\u001b[0m \u001b[38;5;66;03m# 0.6.0a6\u001b[39;00m\n\u001b[0;32m 87\u001b[0m \u001b[38;5;66;03m# New: 3.37ms; 3.14ms; 3.16ms\u001b[39;00m\n\u001b[1;32m---> 89\u001b[0m solution \u001b[38;5;241m=\u001b[39m \u001b[43mtest_1layer\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", - "Cell \u001b[1;32mIn[4], line 60\u001b[0m, in \u001b[0;36mtest_1layer\u001b[1;34m()\u001b[0m\n\u001b[0;32m 17\u001b[0m volume_array, mass_array, gravity_array \u001b[38;5;241m=\u001b[39m \\\n\u001b[0;32m 18\u001b[0m calculate_mass_gravity_arrays(radius_array, density_array)\n\u001b[0;32m 20\u001b[0m input_dict \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mdict\u001b[39m(\n\u001b[0;32m 21\u001b[0m radius_array\u001b[38;5;241m=\u001b[39mradius_array,\n\u001b[0;32m 22\u001b[0m density_array\u001b[38;5;241m=\u001b[39mdensity_array,\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 57\u001b[0m perform_checks \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[0;32m 58\u001b[0m )\n\u001b[1;32m---> 60\u001b[0m solution \u001b[38;5;241m=\u001b[39m \u001b[43mradial_solver\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43minput_dict\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 62\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mResult Success:\u001b[39m\u001b[38;5;124m\"\u001b[39m, solution\u001b[38;5;241m.\u001b[39msuccess)\n\u001b[0;32m 63\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mResult Message:\u001b[39m\u001b[38;5;124m\"\u001b[39m, solution\u001b[38;5;241m.\u001b[39mmessage)\n", - "File \u001b[1;32m~\\anaconda3\\envs\\tpy6py311\\Lib\\site-packages\\TidalPy\\RadialSolver\\solver.pyx:379\u001b[0m, in \u001b[0;36mTidalPy.RadialSolver.solver.radial_solver\u001b[1;34m()\u001b[0m\n", - "\u001b[1;31mTypeError\u001b[0m: radial_solver() takes at least 10 positional arguments (7 given)" + "name": "stdout", + "output_type": "stream", + "text": [ + "Result Success: True\n", + "Result Message: RadialSolver.ShootingMethod:: completed without any noted issues.\n", + "\n", + "Shape: (6, 10).\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuQAAALQCAYAAADYT2XeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAADvDUlEQVR4nOzdd1wT9/8H8FdYYQ9lOQD3APfGBSpKXXVvrXvParVih/pt1ap11brb4qjWqnVv3FUpKu6FC9yAi+Fgv39/fH6JRlATSHJJeD8fjzxM7j65vC/e++7N5XOfkxERgTHGGGOMMSYJM6kDYIwxxhhjLD/jgpwxxhhjjDEJcUHOGGOMMcaYhLggZ4wxxhhjTEJckDPGGGOMMSYhLsgZY4wxxhiTEBfkjDHGGGOMSYgLcsYYY4wxxiTEBTljjDHGGGMS4oKcqYiJiYFMJsPKlSuV06ZMmQKZTJar5QUGBiIwMFA7wb1HJpNhypQpOlk2Y6bAmPKZMWb6Vq5cCZlMhpiYGKlDMThckBsBxQaseFhYWKBIkSLo06cPHj58KHV4WlGsWDHl+pmZmcHZ2RkVK1bEoEGDEBERIXV4jGlNfsjntLQ0LFiwAFWrVoWjoyOcnZ3h5+eHQYMG4fr168p2J0+exJQpU5CQkCBdsIzp2bv5/7HHkSNHpA71o44cOaISr1wuh4eHBwIDAzF9+nQ8efJE6hCNioXUATD1/e9//0Px4sWRkpKC//77DytXrsTx48dx+fJlWFtb6+xzv/32W0ycOFFny1eoUqUKxo0bBwBITk7GtWvXsHHjRqxYsQJffvkl5s6dq9L+zZs3sLDgTZgZJ1PO5w4dOmDPnj3o1q0bBg4ciPT0dFy/fh07d+5E3bp1Ua5cOQCiIJ86dSr69OkDZ2dnncbEmKFYs2aNyuvVq1cjLCws2/Ty5cvrM6xcGzVqFGrWrInMzEw8efIEJ0+exOTJkzF37lxs2LABjRs3Vrbt1asXunbtCrlcLmHEhomrGSPSvHlz1KhRAwAwYMAAuLq6YubMmdi+fTs6d+6ss8+1sLDQS+FbpEgR9OzZU2XazJkz0b17d8ybNw+lS5fG0KFDlfN0WbQwpmumms+nT5/Gzp07MW3aNEyaNEll3q+//prrs+FZWVlIS0vjvGdG7/3j3H///YewsLBs041FgwYN0LFjR5VpFy5cQLNmzdChQwdcvXoVhQoVAgCYm5vD3NxcijANHndZMWINGjQAANy+fVs5LS0tDd9//z2qV68OJycn2NnZoUGDBjh8+HC29yckJKBPnz5wcnKCs7MzevfunePBMqc+p6GhoWjcuDHc3d0hl8vh6+uLJUuWaHcFAdjY2GDNmjUoUKAApk2bBiJSznu/D3lycjLGjBmDYsWKQS6Xw93dHU2bNsXZs2eVbQIDA1GhQgVERkaibt26sLGxQfHixbF06VKtx86YJkwlnxXx16tXL9s8c3NzFCxYUBnH+PHjAQDFixdX/uyt6Fsqk8kwYsQIrF27Fn5+fpDL5di7dy8A4OHDh+jXrx88PDwgl8vh5+eHP/74I9vnLVy4EH5+frC1tYWLiwtq1KiBdevWKeers89gTArq5mSxYsXQqlUrHD9+HLVq1YK1tTVKlCiB1atXZ2t78eJFBAQEwMbGBkWLFsWPP/6I0NBQnfTprly5MubPn4+EhAT8+uuvyuk59SE/c+YMgoOD4erqqjwm9+vXTzlfcS3Mzz//jHnz5sHHxwc2NjYICAjA5cuXtRq3lPgMuRFTbNAuLi7KaUlJSfjtt9+UPxUnJyfj999/R3BwME6dOoUqVaoAAIgIbdq0wfHjxzFkyBCUL18eW7ZsQe/evdX67CVLlsDPzw+ff/45LCwssGPHDgwbNgxZWVkYPny4VtfT3t4e7dq1w++//46rV6/Cz88vx3ZDhgzBpk2bMGLECPj6+uLZs2c4fvw4rl27hmrVqinbvXjxAi1atEDnzp3RrVs3bNiwAUOHDoWVlZXKToAxfTKVfPbx8QEArF27FvXq1fvg2fj27dvjxo0b+OuvvzBv3jy4uroCANzc3JRtDh06hA0bNmDEiBFwdXVFsWLFEBcXhzp16igLdjc3N+zZswf9+/dHUlISxowZAwBYsWIFRo0ahY4dO2L06NFISUnBxYsXERERge7duwNQf5/BmL5pkpO3bt1Cx44d0b9/f/Tu3Rt//PEH+vTpg+rVqyuPlw8fPkSjRo0gk8kQEhICOzs7/PbbbzrtOqKIaf/+/Zg2bVqObeLj49GsWTO4ublh4sSJcHZ2RkxMDDZv3pyt7erVq5GcnIzhw4cjJSUFCxYsQOPGjXHp0iV4eHjobD30hpjBCw0NJQB04MABevLkCd2/f582bdpEbm5uJJfL6f79+8q2GRkZlJqaqvL+Fy9ekIeHB/Xr1085bevWrQSAZs2apfLeBg0aEAAKDQ1VTp88eTK9v6m8fv06W5zBwcFUokQJlWkBAQEUEBDwyXX08fGhli1bfnD+vHnzCABt27ZNOQ0ATZ48WfnaycmJhg8f/tHPCQgIIAA0Z84c5bTU1FSqUqUKubu7U1pa2idjZSwvTD2fs7KylHnm4eFB3bp1o0WLFtHdu3eztZ09ezYBoOjo6GzzAJCZmRlduXJFZXr//v2pUKFC9PTpU5XpXbt2JScnJ+W6tGnThvz8/D4aqzr7DMZ0bfjw4bnOSR8fHwJAx44dU06Lj48nuVxO48aNU04bOXIkyWQyOnfunHLas2fPqECBAh/MwU85fPgwAaCNGzd+sE3lypXJxcVF+Vqx/1N83pYtWwgAnT59+oPLiI6OJgBkY2NDDx48UE6PiIggAPTll19qHLsh4i4rRiQoKAhubm7w8vJCx44dYWdnh+3bt6No0aLKNubm5rCysgIg+lw+f/4cGRkZqFGjhsrPsLt374aFhYVKn2xzc3OMHDlSrVhsbGyUzxMTE/H06VMEBATgzp07SExMzOuqZmNvbw9A/MT8Ic7OzoiIiMCjR48+uiwLCwsMHjxY+drKygqDBw9GfHw8IiMjtRMwY59gqvksk8mwb98+/Pjjj3BxccFff/2F4cOHw8fHB126dNGoD3lAQAB8fX2Vr4kI//zzD1q3bg0iwtOnT5WP4OBgJCYmKr8XZ2dnPHjwAKdPn/7g8tXdZzCmb5rkpK+vr7LLGyB+ZSpbtizu3LmjnLZ37174+/srf1UDgAIFCqBHjx66WwmIY/enjtsAsHPnTqSnp390WW3btkWRIkWUr2vVqoXatWtj9+7dWolValyQG5FFixYhLCwMmzZtQosWLfD06dMcf25atWoVKlWqBGtraxQsWBBubm7YtWuXShLfvXsXhQoVUha6CmXLllUrlhMnTiAoKAh2dnZwdnaGm5ub8gIuXRTkL1++BAA4ODh8sM2sWbNw+fJleHl5oVatWpgyZYrKDkmhcOHCsLOzU5lWpkwZAOCxUZnemHI+y+VyfPPNN7h27RoePXqEv/76C3Xq1FF2P1FX8eLFVV4/efIECQkJWL58Odzc3FQeffv2BSB+AgeAr7/+Gvb29qhVqxZKly6N4cOH48SJEyrLU3efwZi+aZKT3t7e2d7v4uKCFy9eKF/fvXsXpUqVytYup2na9PLly48etwMCAtChQwdMnToVrq6uaNOmDUJDQ5GampqtbenSpbNNK1OmjMkct7kgNyK1atVCUFAQOnTogO3bt6NChQro3r27slgFgD///BN9+vRByZIl8fvvv2Pv3r0ICwtD48aNkZWVpZU4bt++jSZNmuDp06eYO3cudu3ahbCwMHz55ZcAoLXPeZfiwo2P7Tw6d+6MO3fuYOHChShcuDBmz54NPz8/7NmzR+vxMJZX+SWfCxUqhK5du+LYsWMoXbo0NmzYgIyMDLXe++5Zwndj6dmzJ8LCwnJ8KC4mLV++PKKiorB+/XrUr18f//zzD+rXr4/Jkycrl8f7DGaINM3JD41aQu8MgiCF9PR03Lhx46PHbZlMhk2bNiE8PBwjRoxQXrBdvXp1lX1hfsAXdRopc3NzzJgxA40aNcKvv/6qHFd406ZNKFGiBDZv3qwyksK7ByFAXHh18OBBvHz5UuWsWlRU1Cc/e8eOHUhNTcX27dtV/jLPaeQHbXj58iW2bNkCLy+vT47LWqhQIQwbNgzDhg1DfHw8qlWrhmnTpqF58+bKNo8ePcKrV69UzpLfuHEDgLhinTF9yw/5bGlpiUqVKuHmzZt4+vQpPD09Nb5jqJubGxwcHJCZmYmgoKBPtrezs0OXLl3QpUsXpKWloX379pg2bRpCQkKUwyeqs89gTJ90kZM+Pj64detWtuk5TdOWTZs24c2bNwgODv5k2zp16qBOnTqYNm0a1q1bhx49emD9+vUYMGCAss3Nmzezve/GjRsmc9zmM+RGLDAwELVq1cL8+fORkpIC4O1fyu/+ZRwREYHw8HCV97Zo0QIZGRkqwyhlZmZi4cKFn/zcnD4jMTERoaGhuV+ZD3jz5g169eqF58+f45tvvvngATwzMzPbz3ju7u4oXLhwtp++MjIysGzZMuXrtLQ0LFu2DG5ubqhevbpy+vXr13Hv3j0trg1jH2Yq+Xzz5s0c8yYhIQHh4eFwcXFRjqSi+KNY3X7l5ubm6NChA/75558chzt7986Az549U5lnZWUFX19fEBHS09M12mcwpk+6OMYGBwcjPDwc58+fV057/vw51q5dm63t48ePcf369U/26f6YCxcuYMyYMXBxcfnoSE0vXrzIdiZf0c/9/TzcunWryt2MT506hYiICJU/nhMTE3H9+nWddJ3VNT5DbuTGjx+PTp06YeXKlRgyZAhatWqFzZs3o127dmjZsiWio6OxdOlS+Pr6qvz807p1a9SrVw8TJ05ETEwMfH19sXnzZrU24mbNmsHKygqtW7fG4MGD8fLlS6xYsQLu7u54/Phxrtfl4cOH+PPPPwGIs+JXr17Fxo0bERsbi3HjxqlciPm+5ORkFC1aFB07dkTlypVhb2+PAwcO4PTp05gzZ45K28KFC2PmzJmIiYlBmTJl8Pfff+P8+fNYvnw5LC0tle3Kly+PgIAAg799MTMdppDPFy5cQPfu3dG8eXM0aNAABQoUwMOHD7Fq1So8evQI8+fPVxYcij+Av/nmG3Tt2hWWlpZo3bp1tms83vXTTz/h8OHDqF27NgYOHAhfX188f/4cZ8+exYEDB/D8+XPlenl6eqJevXrw8PDAtWvX8Ouvv6Jly5ZwcHBAQkKC2vsMxvRJF8fYCRMm4M8//0TTpk0xcuRI5bCH3t7eeP78ucrJrpCQEKxatQrR0dFqnX3+999/kZKSgszMTDx79gwnTpzA9u3b4eTkhC1btsDT0/OD7121ahUWL16Mdu3aoWTJkkhOTsaKFSvg6OiIFi1aqLQtVaoU6tevj6FDhyI1NRXz589HwYIFMWHCBGWbLVu2oG/fvggNDUWfPn00/p4kJdXwLkx9imGCchoWKDMzk0qWLEklS5akjIwMysrKounTp5OPjw/J5XKqWrUq7dy5k3r37k0+Pj4q73327Bn16tWLHB0dycnJiXr16kXnzp1Ta5i07du3U6VKlcja2pqKFStGM2fOpD/++CPb8EmaDHsIgACQTCYjR0dH8vPzo4EDB1JERESO78E7wx6mpqbS+PHjqXLlyuTg4EB2dnZUuXJlWrx4scp7AgICyM/Pj86cOUP+/v5kbW1NPj4+9Ouvv+a4fHViZ0wTpp7PcXFx9NNPP1FAQAAVKlSILCwsyMXFhRo3bkybNm3K1v6HH36gIkWKkJmZmcrnAfjgkIRxcXE0fPhw8vLyIktLS/L09KQmTZrQ8uXLlW2WLVtGDRs2pIIFC5JcLqeSJUvS+PHjKTExkYjU32cwpms5DXuobk5+aMjgnHL13Llz1KBBA5LL5VS0aFGaMWMG/fLLLwSAYmNjle169+6t1lCIimEPFQ9LS0tyc3Ojhg0b0rRp0yg+Pj7be94f9vDs2bPUrVs38vb2JrlcTu7u7tSqVSs6c+aM8j2KYQ9nz55Nc+bMIS8vL5LL5dSgQQO6cOFCjst/d59nLGREEvf6Z0yPAgMD8fTpU5O6uxdjjDGWG2PGjMGyZcvw8uVLg72lfUxMDIoXL47Zs2fjq6++kjocneE+5IwxxhhjJu7Nmzcqr589e4Y1a9agfv36BluM5yfch5wxxhhjzMT5+/sjMDAQ5cuXR1xcHH7//XckJSXhu+++kzo0Bi7IGWOMMcZMXosWLbBp0yYsX74cMpkM1apVw++//46GDRtKHRoDwH3IGWOMMcYYkxD3IWeMMcYYY0xCXJAzxhhjjDEmIe5DroasrCw8evQIDg4OGt/qmbHcIiIkJyejcOHCMDPjv511jfOc6RvnuH5xjjN90yTHuSBXw6NHj+Dl5SV1GCyfun//PooWLSp1GCaP85xJhXNcPzjHmVTUyXEuyNXg4OAAQHyhjo6OEkfD8oukpCR4eXkptz+mW5znTN84x/WLc5zpmyY5zgW5GhQ/bTk6OnISM73jn1b1g/OcSYVzXD84x5lU1Mlx7rTGGGOMMcaYhLggZ4wxxhhjTEJckDPGGGOMMSYhyQvyhw8fomfPnihYsCBsbGxQsWJFnDlzRjmfiPD999+jUKFCsLGxQVBQEG7evKmyjOfPn6NHjx5wdHSEs7Mz+vfvj5cvX6q0uXjxIho0aABra2t4eXlh1qxZelk/xvI7znHGTBvnOGN5J2lB/uLFC9SrVw+WlpbYs2cPrl69ijlz5sDFxUXZZtasWfjll1+wdOlSREREwM7ODsHBwUhJSVG26dGjB65cuYKwsDDs3LkTx44dw6BBg5Tzk5KS0KxZM/j4+CAyMhKzZ8/GlClTsHz5cr2uL2P5Dec4Y6aNc5wxLSEJff3111S/fv0Pzs/KyiJPT0+aPXu2clpCQgLJ5XL666+/iIjo6tWrBIBOnz6tbLNnzx6SyWT08OFDIiJavHgxubi4UGpqqspnly1bVq04ExMTCQAlJiZqtH6M5YUpbHfGkuNEpvF9M+NiCtsc5zhjH6bJNifpGfLt27ejRo0a6NSpE9zd3VG1alWsWLFCOT86OhqxsbEICgpSTnNyckLt2rURHh4OAAgPD4ezszNq1KihbBMUFAQzMzNEREQo2zRs2BBWVlbKNsHBwYiKisKLFy90vZqMKe3bB8yYAZw/L3Uk+mEKOZ6RAQwbBvz7L0CUp0UxZnJMIcdNyePHwODBwL17UkfCNCVpQX7nzh0sWbIEpUuXxr59+zB06FCMGjUKq1atAgDExsYCADw8PFTe5+HhoZwXGxsLd3d3lfkWFhYoUKCASpuclvHuZ7wrNTUVSUlJKg/GtOG334BJk4Bjx6SORD8MNccB9fN8505gyRKgYUPAzw9YsADg4z9jginkuCn58Udg+XKgdGngyy+BJ0+kjoipS9KCPCsrC9WqVcP06dNRtWpVDBo0CAMHDsTSpUulDAszZsyAk5OT8sG32mXaotiU8svZC0PNcUD9PC9TBhgwALCzA65dA8aMAQoXBvr0AcLD+aw5y99MIcdNSc+eQGAgkJYGzJ8PlCgBfP89kJgodWTsUyQtyAsVKgRfX1+VaeXLl8e9/69WPD09AQBxcXEqbeLi4pTzPD09ER8frzI/IyMDz58/V2mT0zLe/Yx3hYSEIDExUfm4f/9+bleRMRXe3uLf/LJJGWqOA+rnua8vsGIF8OgRsHgxUKkSkJICrFoF1K0LVK4MLFrEBzyWP5lCjpsSf3/g0CHRPbJ6deDlS+CHH0Rh/vPPwJs3UkfIPkTSgrxevXqIiopSmXbjxg34+PgAAIoXLw5PT08cPHhQOT8pKQkRERHw9/cHAPj7+yMhIQGRkZHKNocOHUJWVhZq166tbHPs2DGkp6cr24SFhaFs2bIqV4IryOVy5a11+Ra7TJvy2xlyQ81xQPM8d3QEhg4V/f/Dw8UZchsb4NIlYMQIcdZ8wADg9Gk+a87yD1PKcVMhkwHNmol90aZNQLlywPPnwPjxQKlSwLJlwDtfIzMUerjI9INOnTpFFhYWNG3aNLp58yatXbuWbG1t6c8//1S2+emnn8jZ2Zm2bdtGFy9epDZt2lDx4sXpzZs3yjafffYZVa1alSIiIuj48eNUunRp6tatm3J+QkICeXh4UK9evejy5cu0fv16srW1pWXLlqkVJ1+ZzbTl1CkigKhIkU+3NYXtzlhynCh33/fz50QLFhD5+or/V8WjalWiZcuIkpLUXhTLhzjHDT/HTUF6OtEffxB5e7/dR5UqRbRuHVFmptTRmTZNtjlJC3Iioh07dlCFChVILpdTuXLlaPny5Srzs7Ky6LvvviMPDw+Sy+XUpEkTioqKUmnz7Nkz6tatG9nb25OjoyP17duXkpOTVdpcuHCB6tevT3K5nIoUKUI//fST2jHm1yRm2vf4sdgZymREaWkfb2sq250x5DhR3r7vrCyif/8l6tGDSC5/e9CztycaPJjo7FmNF8nyAc7xtww9x01BSoo4geDm9nYfVakS0Y4dYh/GtE+TbU5GxD+ufkpSUhKcnJyQmJiYb37yYrqRlSW6OaSlATExwP//qpsj3u70S1vf97Nnon/5smXAjRtvp9eqJYYj69JFXCDKGOe4fvH3Lbx8KS74nD0bUAw8U7cuMH06EBAgaWgmR5NtTtI+5IzlN2ZmQNGi4nl+6Uee3xQsCIwdC1y/Dhw+LApwS0vg1Cmgf3/R13zECNH3nDHG9M3eHvj2WyA6GpgwAbC2Bk6eFKOzfPYZcPas1BHmT1yQM6Zn+W2klfxKJhMHuPXrgQcPgJkzxUgHSUliVJZKlYB69YDVq3nkA8aY/hUoIPZLt2+LC9YtLN6OztKpkzipwPSHC3LG9Cy/jbTCAHd3cSbq5k1g/36gQwdx8Dt5EujdGyhSRNzEgw+AjDF9K1xYDOl6/TrQo4c4mbBpk7gRWv/+fKzSFy7IGdMzPkOef5mZAU2bioPdvXvirno+PuLOn/PnA+XLiz6c69YBqalSR8sYy09KlgT+/BO4cAH4/HNxzdMff7y96+d7Q8UzLeOCnDE94zPkDAAKFQK++Ub8XLx7N9CmjSjYjx0TZ6mKFhXjBt+8KXWkjLH8pGJFYNu2t/3KFXf9LFmS7/qpS1yQM6ZnfIacvcvcHGjeHNi6Fbh7F5gyRRTjT5+KO+uVKQMEBQEbN4oDI2OM6YPirp/79wM1avBdP3WNC3LG9ExRkPMZcva+okWByZPF6AfbtwMtWoj+nAcPAp07i19XQkKAO3ekjpQxlh/IZKKb3alTfNdPXeOCnDE9U3RZefFCnHFg7H0WFkDr1sCuXaI4//ZbwNNT9OH86SdxIPzsM2DLFj4YMsZ0TyYTF6NfugSEhooTS48eAUOGiGtf1q0Tfc5Z7nFBzpieOToCTk7iOXdbYZ/i4yN+Jr53D/jnH6BZM3GPvX37gPbtxfzvv+dfXBhjumdhAfTpI2569ssvYgSp27fFdS9VqwI7d4r9E9McF+SMSUBxlpwLcqYuS0tRgO/bJw6AEycCbm7A48eiYC9eHGjVCtixA8jMlDpaxpgpk8uBkSPFvujHH8VJposXxS979esDR49KHaHx4YKcMQlwP3KWFyVKADNmiBsO/f030Lix+Ll41y4xXFnx4sD//gc8fCh1pIwxU2ZvL0aLunMH+PprwMaG7/qZW1yQMyYBPkPOtMHKSlzsefAgEBUFjBsHFCwotqvJk0V3lrZtgb17uX8nY0x3ChQQ17fcusV3/cwtLsgZkwCfIWfaVqaMGIrswQNg7VqgQQPRdWXbNjGsYsmSwPTpQGys1JEyxkyV4q6fUVFAz558109NcEHOmAT4DDnTFWtroHt3cYOhK1eA0aMBZ2cgJkb8tOzlJc5YHTjAZ80ZY7pRogSwZo2462ebNqp3/Rwzhu/6mRMuyBmTAJ8hZ/rg6yvusPfoEbBypbjRR0aGOGPVtKk4qz5rFvDkidSRMsZMUcWK4qZn4eFAo0bi5mYLFoiC/bvv+K6f7+KCnDEJvHuGnIeIYrpmYwP07i0utrpwARg2TAy/efu2uBCrSBGgWzfgyBHeHhlj2lenjrjWRXHXz1evxOgsJUoAs2fzXT8BLsgZk0SRIqJvXUqKuEU6Y/pSqRKwaJE4a/7bb0DNmuLmQuvXizNY5csD8+aJu/Exxpi2vHvXz3/+Efua58+BCRPEzc6WLs3fNzrjgpwxCcjlgIeHeM79yJkU7OzERVanTgGRkcCgQWJaVBQwdqy4OKtXL+D4cT5rzhjTHplM3FPh0iXRlc7HR5wgGDo0f9/1kwtyxiTC/ciZoahWDVi2TBwUlywBqlQBUlOBP/8Uo7VUrAgsXMg/KzPGtMfcXHSli4rKftfPKlWAPXukjlC/uCBnTCJubuLfx4+ljYMxBUdHYMgQcTOPiAigXz/R//zKFWDUKDF02c6dUkfJGDMlirt+3rkDTJsm7vp56RLQogUwZUr++YWOC3LGJEAEnDsnnpcpI20sjL1PJgNq1QJ+/12cNV+4UFz3EB0tbo3dpo14zhhj2mJnB0yaJArzESPEtKlTxVn0tDRpY9MHLsgZk0BUlCh05HKgbl2po2Hsw5ydxcHx+nVx8ZWFBbB9uxhS8YcfxIXJjDGmLQUKiJMAy5eLbi1r1gDBwcCLF1JHpltckDMmgUOHxL9164ouAYwZOnt7YOZMMWxio0aiEP/+e9G/fO9eqaNjjJmagQOBXbsABwcxJGvduqb9yxwX5IxJ4OBB8W+TJtLGwZimfH3F9vvXX0ChQsCtW0Dz5kCHDnyBMmNMu4KDxUhPRYuKX+nq1BHXt5giLsgZ07OsLODwYfG8cWNpY2EsN2QyoGtXcYAcO1b8rLx5M1CuHDBjhhihhTHGtKFSJVGEV60KxMcDgYFif2NquCBnTM/Onxd94RwcxE1ZGDNWjo7AnDlim27YUAyLOGmSOIAeOCB1dIwxU1G4MHDsGNCypegu17EjMHeuaY3AwgU5Y3qm6D/esKG4QI4xY1ehgujjuWaNuOHVjRvijnydOwMPHkgdHWPMFNjbA1u3AsOGiUJ83DhxwXlGhtSRaQcX5IzpGfcfZ6ZIJgN69hQjCI0eDZiZARs3im4ss2fnj2HLGGO6ZWEB/Pqr+GVOJgMWLxbDsL58KXVkeSdpQT5lyhTIZDKVR7ly5ZTzU1JSMHz4cBQsWBD29vbo0KED4uLiVJZx7949tGzZEra2tnB3d8f48eOR8d6fS0eOHEG1atUgl8tRqlQprFy5Uh+rx1g2aWnAv/+K5/mhIOccz3+cnID588XNherWBV69EsMlVqny9toJZjo4x5m+yWTi2pVNmwBra2D3bnFH4YcPpY4sbyQ/Q+7n54fHjx8rH8ePH1fO+/LLL7Fjxw5s3LgRR48exaNHj9C+fXvl/MzMTLRs2RJpaWk4efIkVq1ahZUrV+L7779XtomOjkbLli3RqFEjnD9/HmPGjMGAAQOwb98+va4nYwBw6pQoUFxdxc/8+QHneP5UubL443PlSnFX2mvXxEXM3buLMfiZ6eAcZ1Jo3150lXN3F9ex1KkDXLwodVR5QBKaPHkyVa5cOcd5CQkJZGlpSRs3blROu3btGgGg8PBwIiLavXs3mZmZUWxsrLLNkiVLyNHRkVJTU4mIaMKECeTn56ey7C5dulBwcLDacSYmJhIASkxMVPs9jOVk6lQigKhz50+3NYXtzlhynMg0vm9D9fw50fDhRGZmYvu3tyeaM4coLU3qyKRlCtsc5ziT2p07ROXKiX2LgwPR3r1SR/SWJtuc5GfIb968icKFC6NEiRLo0aMH7v3/QLaRkZFIT09HUFCQsm25cuXg7e2N8PBwAEB4eDgqVqwIDw8PZZvg4GAkJSXhypUryjbvLkPRRrEMxvRJ0X88Pw13yDnOXFxEv8/Tp4HatUV/z3HjgGrVxMgJzLhxjjMpFS8OnDwphkNMThYjsSxfLnVUmpO0IK9duzZWrlyJvXv3YsmSJYiOjkaDBg2QnJyM2NhYWFlZwdnZWeU9Hh4eiI2NBQDExsaqJLFivmLex9okJSXhzZs3OcaVmpqKpKQklQdjefX6NaA4fuSH/uOA4eY4wHkuhWrVxIHzt9+AggWBy5eBgACgVy/g//87mZHhHGeGwMUF2LdP7EsyM4HBg4GJE8V9P4yFpIOuNW/eXPm8UqVKqF27Nnx8fLBhwwbYSHg/8RkzZmDq1KmSfT4zTcePA+npgJcXULKk1NHoh6HmOMB5LhUzM6B/f6BdO+Cbb4Bly4A//wS2bwd++EEMacbDgRoPznFmKKysgFWrxPF1yhRg5kzgzh0xTeJNUS2Sd1l5l7OzM8qUKYNbt27B09MTaWlpSEhIUGkTFxcHT09PAICnp2e2q7UVrz/VxtHR8YM7i5CQECQmJiof9+/f18bqsXzu3eEOZTJpY5GKoeQ4wHkutQIFgCVLxB34atQAkpLEcInVqwMnTkgdHcstznEmJZkMmDxZFOGWlmLo1SZNgCdPpI7s0wyqIH/58iVu376NQoUKoXr16rC0tMRBRRUDICoqCvfu3YO/vz8AwN/fH5cuXUJ8fLyyTVhYGBwdHeHr66ts8+4yFG0Uy8iJXC6Ho6OjyoOxvFLcECg/9R9/n6HkOMB5bihq1gT++0+cKXdxEaMk1K8P9O0rbpPNjAvnODMEX3wB7N8PODuLrqL+/uKGZQZNDxeZftC4cePoyJEjFB0dTSdOnKCgoCBydXWl+Ph4IiIaMmQIeXt706FDh+jMmTPk7+9P/v7+yvdnZGRQhQoVqFmzZnT+/Hnau3cvubm5UUhIiLLNnTt3yNbWlsaPH0/Xrl2jRYsWkbm5Oe3V4DJcvjKb5dXz50QymbgK/MED9d5jCtudseQ4kWl838buyROiAQNEngBEzs5Ev/5KlJEhdWS6YQrbHOc4M2RXrxIVLy72JwUKEB07pt/P12Sbk7Qg79KlCxUqVIisrKyoSJEi1KVLF7p165Zy/ps3b2jYsGHk4uJCtra21K5dO3r8+LHKMmJiYqh58+ZkY2NDrq6uNG7cOEpPT1dpc/jwYapSpQpZWVlRiRIlKDQ0VKM4OYlZXm3ZInYIZcuq/x5T2O6MJceJTOP7NhXh4UTVqr0tzKtWFdNMjSlsc5zjzNDFxRHVri32JVZWRGvX6u+zNdnmZERE0pybNx5JSUlwcnJCYmIi/+TFcmXkSDHs27BhwKJF6r2Htzv94u/bsGRmim4s33wDKLogDxgAzJghbqxlCnib0y/+vvOv16/FCCybN4vXP/4ITJqk++u5NNnmDKoPOWOmivuPM6YZc3PxB2xUFNCnj5j2229AmTKiUM/MlDQ8xpgRsbUVF3iOGydef/ut+AM/PV3auN7FBTljOvb4MXD1qvhLvFEjqaNhzLi4uwOhoWLY0MqVgRcvgCFDxG2yT5+WOjrGmLEwMwN+/ln8Sm1mBvzxB9C8+dtf4KTGBTljOnb4sPi3alUx1BtjTHP16gFnzgALFgCOjuJ57dqiOH/2TOroGGPGYtgwcd8DOzsxHHH9+sDdu1JHxQU5YzqnGK2Lu6swljcWFsCoUaIbS69e4pLPZcuAsmWB3383rrvyMcak07Il8O+/QOHCwJUr4he3M2ekjYkLcsZ0TNF/vEkTaeNgzFR4egKrVwNHjwIVKogz5AMGAHXrAmfPSh0dY8wYVK0qbkxWqRIQGwsEBIgz51LhgpwxHbpzB4iJEWf26teXOhrGTEvDhqIAnzsXcHB4e9fP4cOBxESpo2OMGbqiRcWZ8uBgMRJL27bAL79IEwsX5Izp0JYt4t86dQB7e2ljYcwUWVoCX34JXL8OdOsmurEsXizOdr14IXV0jDFD5+gI7NgBDBok9h+jRwO7d+s/Di7IGdORxEQxZjIgbuPLGNOdwoWBdetEFzEPD+DCBTGCQnKy1JExxgydpSWwdKm44BMA/vc/UZzrExfkjOnIrFmib2u5ckDfvlJHw1j+0KgRcOCAGNEoIgL4/HPgzRupo2KMGTqZDPj+e8DaWuw7FAMy6AsX5IzpwMOHwLx54vlPP4k+5Iwx/ahQAdi/X/wUfeQI0KEDkJYmdVSMMUPn4SG6rgDibp76xAU5YzowZYo4K1evnjhDxxjTr+rVgV27ABsbYM8eoHt3ICND6qgYY4Zu/HjRheXoUXFDMn3hgpwxLbt6VdwBDBDdVmQyaeNhLL+qXx/Ytg2wsgL++Qfo35/HKmeMfVzRokCfPuL5tGn6+1wuyBnTsokTxUG/XTsxLjJjTDpNmwIbNwLm5mLs8hEj9H+xFmPMuEycKPYZe/fq74ZBXJAzpkX//iuGTzI3fzvCCmNMWp9/DqxZI36tWrIEmDCBi3LG2IeVKCGGUQX0d5acC3LGtIRIHOgBcdfAsmWljYcx9la3bsCKFeL5zz8DP/wgbTyMMcMWEiL+iN+6Fbh8WfefxwU5Y1qyeTPw33+ArS0webLU0TDG3te/PzB/vng+ebK4wydjjOXE11eM0AQA06fr/vO4IGdMC9LTxV/TADBuHFCokLTxMMZyNnr025+gx40Dli+XNh7GmOH65hvx799/Azdv6vazuCBnTAt++00kq5ubGDKJMWa4Jk0SF20BwJAhwJ9/ShsPY8wwVakCtGwpBmr46SfdfhYX5Izl0cuXwNSp4vnkyYCDg7TxMMY+bfr0tyOu9OkDbNkidUSMMUOkOEu+ejVw967uPket+wf+8ssvGi+4b9++cODKhOUDc+YAcXFAqVJv7/BlbDjHWX4jkwELFgCvXgGhoUCXLsD27cBnn0kdme5wnjOmOX9/oEkT4OBBcW+RRYt08zkyok8P/mRmZoaiRYvC3NxcrYXev38fN27cQIkSJfIcoCFISkqCk5MTEhMT4ejoKHU4zIDExQElS4qD+oYNQKdO2lu2Pre7/J7jAOd5fpWZCfToIfqIWluLcYcDAvTz2fre5vJ7nnOOs9w6fBho3BiQy4HoaPWvE9Nkm1PrDDkAnDlzBu7u7mq15b+mWX4xdaooxmvVAjp2lDqavOEcZ/mRubkYo/zVK2DnTqBVK3EmrFYtqSPTDc5zxjQXGChu9HfypPhV/Oeftf8ZavUhnzx5Muzt7dVe6KRJk1CgQIFcB8WYMbhx4+0IDbNmiZ/AjRXnOMvPLC3F3TwbNxbXhHz2GXDxotRRaR/nOWO5I5MB334rni9ZAjx9qoPPUKfLSn7HP3OxnHTsCPzzj7gCe+dO7S+ftzv94u+bvXwJBAeLs2Du7sCxY7q9wRdvc/rF3zfLCyKgRg3g7FlxoeePP376PZpsczzKCmO58N9/ohg3M9P9UEiMMf2wtwd27QKqVgXi48WFXNHRUkfFGDMEMtnbEVcWLgQSErS7fI0L8mfPnmH48OHw9fWFq6srChQooPJgzNQRARMmiOe9ewMVKkgbj7ZxjrP8zNkZ2L9f3KXv4UMgKAh49EjqqLSP85wxzbVtK/YNSUnaH21F7Ys6FXr16oVbt26hf//+8PDwgMyYO84ylgs7dwL//itGZPjf/6SORvs4x1l+5+oKhIUBDRsCt2+LovzoUXHjL1PBec6Y5szMxI3FevYE5s0Td/7V4LKMj9K4IP/3339x/PhxVK5cWTsRMGZEMjLe3uFvzBigaFFJw9EJznHGgMKFxWgrDRoA164BzZqJoc+cnaWOTDs4zxnLnS5dxE0Ab98Gli0Dxo3TznI17rJSrlw5vHnzRjuf/o6ffvoJMpkMY8aMUU5LSUnB8OHDUbBgQdjb26NDhw6Ii4tTed+9e/fQsmVL2Nrawt3dHePHj0dGRoZKmyNHjqBatWqQy+UoVaoUVq5cqfX4Wf6wahVw9SpQoADw9ddSR6MbuspxgPOcGRcfH+DAAcDDAzh/HmjeXFz4aQr4WM5Y7lhYACEh4vnPPwMpKVpaMGno1KlT1LhxYzpy5Ag9ffqUEhMTVR65cerUKSpWrBhVqlSJRo8erZw+ZMgQ8vLyooMHD9KZM2eoTp06VLduXeX8jIwMqlChAgUFBdG5c+do9+7d5OrqSiEhIco2d+7cIVtbWxo7dixdvXqVFi5cSObm5rR3716140tMTCQAuV4/ZhpevSIqXJgIIJo7V/efJ9V2p4scVyyX85wZo4sXiVxcRO43akT0+rV2livlNsfHcsZyLzWVyMtL7BMWLfpwO022OY0L8hs3blCNGjXIzMxM5SGTycjMzEzTxVFycjKVLl2awsLCKCAgQJnECQkJZGlpSRs3blS2vXbtGgGg8PBwIiLavXs3mZmZUWxsrLLNkiVLyNHRkVJTU4mIaMKECeTn56fymV26dKHg4GC1Y+QkZkRE06aJ5CtWjCglRfefJ9V2p+0cJ+I8Z8bv1CkiBwexD2jZUhyQ80rKbY6P5YzlzcKFYn/g7U2UlpZzG022OY27rPTo0QOWlpZYt24dDh48iEOHDuHQoUM4fPgwDh06pPEZ+uHDh6Nly5YICgpSmR4ZGYn09HSV6eXKlYO3tzfCw8MBAOHh4ahYsSI8PDyUbYKDg5GUlIQrV64o27y/7ODgYOUyGFPH06fAzJni+Y8/itvnmipt5zjAec6MX82aYkhEGxvxb8+e4poSY8XHcsbypn9/0Z3t3j3gzz/zvjyNL+q8fPkyzp07h7JauFvC+vXrcfbsWZw+fTrbvNjYWFhZWcH5vStoPDw8EBsbq2zzbgIr5ivmfaxNUlIS3rx5Axsbm2yfnZqaitTUVOXrpKQkzVeOmZQJE8QwR1WrAt26SR2NbmkzxwHOc2Y6GjQAtmwBPv9c3NmzenXjvZaEj+WM5Y2Njbigc8IE4Ndfgb5987Y8jc+Q16hRA/fv38/bpwK4f/8+Ro8ejbVr18La2jrPy9OmGTNmwMnJSfnw8vKSOiQmoT//BEJDxU0B5s8Xwx6ZMm3lOMB5zkxPcDDwyy/i+dKlQFaWtPHkFh/LGcs7X1/xrzZ+LdO4tBg5ciRGjx6NlStXIjIyEhcvXlR5qCsyMhLx8fGoVq0aLCwsYGFhgaNHj+KXX36BhYUFPDw8kJaWhoT3boUUFxcHT09PAICnp2e2K7UVrz/VxtHRMce/qAEgJCQEiYmJyoe2ihNmfK5fB4YMEc+//16MS2zqtJXjAOc5M01ffAE4OQExMUAue3FJjo/ljOXd7dvi35IltbAwTTuxy2SybI/cXAiSlJREly5dUnnUqFGDevbsSZcuXVJeCLJp0yble65fv57jhSBxcXHKNsuWLSNHR0dK+f+r7iZMmEAVKlRQ+exu3brxhSDsk16/JqpY8e3IChkZ+v18qbY7beU4Eec5M13Dhol9Q5cuuV+GlNscH8sZy7tRo8R+4Kuvcp6v01FWYmJiPvrIi3evzCYSQyV5e3vToUOH6MyZM+Tv70/+/v7K+Yqhkpo1a0bnz5+nvXv3kpubW45DJY0fP56uXbtGixYt4qGSmFoGDhSJ5u5O9Pix/j9fqu1OlzlOxHnOTENkpNg/WFkRPX2au2VIuc3xsZyxvGvZUuwHlizJeb5OC/KPLfTmzZuaLk7F+0n85s0bGjZsGLm4uJCtrS21a9eOHr9XGcXExFDz5s3JxsaGXF1dady4cZSenq7S5vDhw1SlShWysrKiEiVKUGhoqEZxcRLnP+vWiSSTyYjCwqSJQartTpc5TsR5zkxH1apiPzF/fu7eL+U2x8dyxvKufHmxD9i/P+f5mmxzMiIiTbq4NGjQAGFhYdku3oiKikKTJk3w4MGD3PaeMVhJSUlwcnJCYmIiHB0dpQ6H6diNG2L0hJcvge++A/73P2nikGq7y485DnCeM80tWgSMGAFUqABcvCgu/NaElNtcfsxzznGmTVlZgK0tkJoq+pKXKJG9jSbbnMYXddrb26N9+/Yqt7S9du0aAgMD0aFDB00Xx5hBSUkBOncWxXhAADB5stQR6R/nOGPq6d4dsLYGLl8Gchjxz6BxnjOWN48eiWLc3BzQxgA+GhfkmzdvRmJiInr06AEiwuXLlxEYGIhu3bphwYIFeY+IMQl9+SVw4QLg5gasWycSLb/hHGdMPS4uQMeO4vnvv0sbi6Y4zxnLG8UIKz4+gKVl3pencUFuY2ODXbt2ISoqCp07d0aTJk3wxRdfYO7cuXmPhjEJ/f23GFdYJhNjjxcuLHVE0uAcZ0x9/fuLf//6C3j1StpYNMF5zljeaHXIQ6h5p873725lZmaGv//+G02bNkWHDh3w3XffKdtwvyxmjG7dAgYOFM9DQoBmzaSNR984xxnLnYAAoFQpsQ/ZuBHo00fqiD6M85wx7dF2Qa7WRZ1mZmaQ5XC1iuKtMpkMRASZTIbMzEztRGZA+EIQ05aSAtStC5w7J26NfegQYKHWn6q6pc/tLr/nOMB5znJvxgxg0iSgXj3g+HH136fvbS6/5znnONOmrl3FL+uzZwNffZVzG022ObXKjsOHD2scKGPG4quvRDHu6ir6jRtCMa5vnOOM5V7v3mJEphMnxN19y5WTOqKccZ4zpj2SdFkJCAjQzqcxZmA2bRJDlwHAmjVA0aLSxiMVznHGcq9wYaBFC2DHDnFx5+zZUkeUM85zxrRH2wW5Whd1Xrx4EVlZWWov9MqVKypDKTFmiG7ffntB1tdfA599Jm08UuIcZyxvBgwQ/65eDaSlSRvLh3CeM6YdL16IB5Dz+OO5oVZBXrVqVTx79kzthfr7++PevXu5DooxXUtNBbp0AZKSRP/xH36QOiJpcY4zljctWgCenkB8PLBzp9TR5IzznDHtuHNH/OvhAdjba2eZanVZISJ89913sLW1VWuhaYZ6eoCx/zdhAhAZCRQoAKxfr50xRI0Z5zhjeWNhIUZY+ekn0W2lfXupI8qO85wx7dB2dxVAzYK8YcOGiIqKUnuh/v7+sLGxyXVQjOnS5s3AL7+I56tXa+cOW8aOc5yxvOvXTxTke/cCDx4Y3jUpnOeMaYeiINdWdxVAzYL8yJEj2vtExiQUHS0OmoAYXaVlS2njMRSc44zlXenSYlzyo0eBlSuBb7+VOiJVnOeMaYcuzpBrfKdOxoxVWproN56YCNSpA0yfLnVEjDFTo7hQ/PffAQ2un2SMGREuyBnLg6+/Bk6fBlxcxGD++b3fOGNM+zp0AJycgJgYgIf9Zsw0cUHOWC5t2wbMny+er1wJeHtLGQ1jzFTZ2gLdu4vnv/0mbSyMMe1LTRXXiABckDOmkZgYMfoBAIwdC3z+uZTRMMZMnWJM8s2bAQ1GGWSMGYHoaIAIsLMD3N21t1wuyJlJS0sDunYFEhKAWrWAGTOkjogxZuqqVQOqVBH7n7VrpY6GMaZNijHIS5YEZDLtLVfjgnzVqlXYtWuX8vWECRPg7OyMunXr4u7du9qLjDEtmDQJiIgAnJ1Fv3ErK6kjMnyc44zlneIs+W+/ibNphobznLHc0UX/cSAXBfn06dOV45KGh4dj0aJFmDVrFlxdXfHll19qNzrG8mDnTmDOHPE8NBQoVkzScIwG5zhjede9OyCXA5cuAWfOSB1NdpznjOWOrgpytcYhf9f9+/dRqlQpAMDWrVvRoUMHDBo0CPXq1UNgYKB2o2Msl6KigC++EM9HjwbatpU0HKPCOc5Y3rm4AB07ii4rq1cDNWtKHZEqznPGNEckRmsDDOAMub29PZ79/1Uq+/fvR9OmTQEA1tbWePPmjXajYywXYmOBzz4DXrwQ443PmiV1RMaFc5wx7ahTR/wbFydtHDnhPGdMc8uXAydPiu6vQUHaXbbGZ8ibNm2KAQMGoGrVqrhx4wZatGgBALhy5QqKcZ8AJrHkZHH3zZgYoFQpYPt27jeuKc5xxrRDcWMgMwMcPoHznDHN3LghRmoDxAAR//8Dk9ZovJtYtGgR/P398eTJE/zzzz8oWLAgACAyMhLdunXTbnSMaSA9HejUCTh7FnBzA/buFf8yzXCOM6YdmZniX3NzaePICec5Y+rLyAB69QJevwYaNwbGjNH+Z8iIDPH6b8OSlJQEJycnJCYmwtHRUepwWA6IxC2rQ0PFjTkOHxbDHBoz3u70i79vpm1z5gBffQX07AmsWZN9Pm9z+sXfN8utqVOBKVPEiG0XLwJeXuq9T5NtTuMuK8eOHfvo/IYNG2q6SMbybMoUUYybmYnhDY29GJcS5zhj2mHIZ8g5zxlTT0QE8MMP4vnixeoX45rSuCDP6epr2Tsjo2cq9kCM6cmKFcD//ieeL1kCtGolbTzGjnOcMe0w5D7knOeMfdqrV6KrSmYm0K2beOiKxruJFy9eqDzi4+Oxd+9e1KxZE/v379dFjIx90K5dwNCh4vm33wKDBkkbjyngHGdMOwy5IOc8Z+zTxo0Dbt4EihYFFi3S7WdpfIbcyckp27SmTZvCysoKY8eORWRkpFYCY+xTTp8GOncWf7n26fP2LDnLG85xxrTDkLuscJ4z9nE7dwLLlonnK1eKewvoktb+bvfw8EBUVJS2FsfYR92+LYY3fP0aaNZMjA36zq+tTAc4xxnTjCGfIf8QznPGgPh4MVAEAHz5JdCkie4/U+PdxMWLF1UeFy5cwN69ezFkyBBUqVJFo2UtWbIElSpVgqOjIxwdHeHv7489e/Yo56ekpGD48OEoWLAg7O3t0aFDB8S9d4eFe/fuoWXLlrC1tYW7uzvGjx+PjIwMlTZHjhxBtWrVIJfLUapUKaxcuVLT1WYG5MkTceOfJ0+AqlWBTZsAS0upozIdnOOMaYchnyHXVp5zjjNTQwQMHCiK8goVgOnT9fbBmpHJZGRmZkYymUzl4e/vT9euXdNoWdu3b6ddu3bRjRs3KCoqiiZNmkSWlpZ0+fJlIiIaMmQIeXl50cGDB+nMmTNUp04dqlu3rvL9GRkZVKFCBQoKCqJz587R7t27ydXVlUJCQpRt7ty5Q7a2tjR27Fi6evUqLVy4kMzNzWnv3r1qx5mYmEgAKDExUaP1Y9r36hVR7dpEAJGPD9GjR1JHpDtSbXf5MceJOM+Z9n3zjdhXjRyZ83wptzlt5TnnODM1K1aIvLWyIjp/Pm/L0mSb07ggj4mJUXncu3eP3rx5k6tAc+Li4kK//fYbJSQkkKWlJW3cuFE579q1awSAwsPDiYho9+7dZGZmRrGxsco2S5YsIUdHR0pNTSUiogkTJpCfn5/KZ3Tp0oWCg4PVjomT2DCkpxN9/rlIFBcXIg1rQ6Mj1XaXH3OciPOcad/EiWJ/NWZMzvOl3OZ0meec48xY3bxJZGcn8nbWrLwvT5NtTuMuKz4+PioPLy8vWFtb5/lMfWZmJtavX49Xr17B398fkZGRSE9PR1BQkLJNuXLl4O3tjfDwcABAeHg4KlasCA8PD2Wb4OBgJCUl4cqVK8o27y5D0UaxjJykpqYiKSlJ5cGkRQSMGgVs3w7I5cCOHUC5clJHZZryQ44DnOdM9wy5D7ku8pxznBkzxd04X70CAgKAsWP1+/lqjbLyyy+/YNCgQbC2tsYvv/zy0bajRo3SKIBLly7B398fKSkpsLe3x5YtW+Dr64vz58/DysoKzs7OKu09PDwQGxsLAIiNjVVJYsV8xbyPtUlKSsKbN29gY2OTLaYZM2Zg6tSpGq0H062ffhJjjMtkwLp1QL16UkdkWvJbjgOc50z3DK0Pua7ynHOcmYIZM4D//gMcHYFVq/Sft2oV5PPmzUOPHj1gbW2NefPmfbCdTCbT+GBdtmxZnD9/HomJidi0aRN69+6No0eParQMbQsJCcHYd/40SkpKgpeubs3EPmnNGmDSJPF8wQKgfXtp4zFF+S3HAc5zpnuGdoZcV3nOOc6M3enTgOJvt0WLAB8f/cegVkEeHR2d43NtsLKyQqlSpQAA1atXx+nTp7FgwQJ06dIFaWlpSEhIUPnrOi4uDp6engAAT09PnDp1SmV5iqu3323z/hXdcXFxcHR0/OBf1XK5HHK5XCvrx/ImLAzo1088/+orYORIaeMxVfktxwHOc6Z7ioLcUM6Q6yrPOceZMXv1CujZU/yi1bkz0KOHNHEYyN/tb2VlZSE1NRXVq1eHpaUlDh48qJwXFRWFe/fuwd/fHwDg7++PS5cuIT4+XtkmLCwMjo6O8PX1VbZ5dxmKNoplMMN1/jzQoYPo19W1KzBzptQRMW3gHGf5haLLiqGcIdcXznFmTMaPB27cAIoUeds1VgpqnSEfq0HP9rlz56rdNiQkBM2bN4e3tzeSk5Oxbt06HDlyBPv27YOTkxP69++PsWPHokCBAnB0dMTIkSPh7++POnXqAACaNWsGX19f9OrVC7NmzUJsbCy+/fZbDB8+XPlX8ZAhQ/Drr79iwoQJ6NevHw4dOoQNGzZg165dasfJ9O/ePaBFCyA5GQgMFHfJym8HNX3iHGdM+wyty4ou8pxznBmz3btFEQ6IOqNAAQmDUWfYlsDAQJWHo6Mj2draUtWqValq1apkZ2dHjo6O1KhRI42Gg+nXrx/5+PiQlZUVubm5UZMmTWj//v3K+W/evKFhw4aRi4sL2draUrt27ejx48cqy4iJiaHmzZuTjY0Nubq60rhx4yg9PV2lzeHDh6lKlSpkZWVFJUqUoNDQUI3i5KGS9Ov5c6Ly5cWwQxUqEL14IXVE0tDndpffc5yI85xp3+DBYj82dWrO8/W9zekizznHmbF68oTIw0Pk6OjRuvkMnY5DPmfOHGrdujU9f/5cOe358+fUpk0b+vnnnzVdnFHgJNafN2+IGjQQCVKkCNG9e1JHJB2ptrv8mONEnOdM+wYOFPuyH37Ieb6U21x+zHPOcaaQlUXUrp3IT19fotevdfM5mmxzMiIiTc6oFylSBPv374efn5/K9MuXL6NZs2Z49OhRns/aG5qkpCQ4OTkhMTERjo6OUodjsrKyRF/xjRvFsEPHjwMVK0odlXSk2u7yY44DnOdM+/r3B/74Q9x6OyQk+3wpt7n8mOec40whNFQMGGFpCZw6BVSpopvP0WSb07hnW1JSEp48eZJt+pMnT5CcnKzp4hhT+uorUYxbWgJbt+bvYlxKnOOMaYeh9SF/F+c5y6/u3BE3GgSAH37QXTGuKY13E+3atUPfvn2xefNmPHjwAA8ePMA///yD/v37oz0PEM1yad488QDEhRWNGkkaTr7GOc6YdhjajYHexXnO8qPMTOCLL4CXL4EGDcSJQEOh1igr71q6dCm++uordO/eHenp6WIhFhbo378/Zs+erfUAmenbsOHtLWpnzgS6d5c2nvyOc5wx7TDkM+Sc5yw/mjkTOHECcHAAVq82rD+WNe5DrvDq1Svcvn0bAFCyZEnY2dlpNTBDwv3OdOfYMaBpUyAtDRgxAvjlF+nGADU0Um93+SnHAem/b2Z6uncH/voLmD8fGD06+3xD2ObyU54bwvfNpBMZCdSpI+5tsmqVOFOua5pscxqfIVews7NDpUqVcvt2xnD+PNCmjSjG27UTBy0uxg0H5zhjeWMMNwbiPGf5wY0bQKdOohjv2BHo1UvqiLLLVUF+5swZbNiwAffu3UNaWprKvM2bN2slMGbaLl8GgoKAhASgXj1g7VrD+ukov+McZyzvrl0T/xYuLG0cH8J5zvKDQ4dEEf7iBVCsGLB0qWGe/NP47/b169ejbt26uHbtGrZs2YL09HRcuXIFhw4dgpOTky5iZCbm+nWgSRPg2TOgZk1g1y7AxkbqqJgC5zhjeffihTjxAAD160sbS044z1l+sHw5EBws8tHfH/jvP6BgQamjypnGBfn06dMxb9487NixA1ZWVliwYAGuX7+Ozp07w9vbWxcxMhNy6xbQuDEQHy+GGtq3D+B9v2HhHGcs706eBIiA0qUBDw+po8mO85yZssxM4MsvgcGDRTeV7t3FmXJDzEUFjQvy27dvo2XLlgAAKysrvHr1CjKZDF9++SWWL1+u9QCZ6YiJEcX448dAhQpAWBjg4iJ1VOx9nOOM5d3x4+LfBg2kjeNDOM+ZqUpKEtenzZ8vXv/wA/Dnn4C1taRhfZLGBbmLi4vypgFFihTB5f//TS4hIQGvX7/WbnTMZNy/L8YWv38fKFcOOHAAcHWVOiqWE85xxvLu33/Fv4bYXQXgPGemKSZGXJem6Aq7YQPw7beG2Wf8fRpf1NmwYUOEhYWhYsWK6NSpE0aPHo1Dhw4hLCwMTZo00UWMzMg9eiTOjMfEAKVKAQcPGvbPRvkd5zhjeZOSApw+LZ4b6hlyznNmak6eBNq2BZ48AQoVArZvB2rUkDoq9WlckP/6669ISUkBAHzzzTewtLTEyZMn0aFDB3z77bdaD5AZt7g4cQHnrVvi6uZDhwx3xAEmcI4zljenT4vhXD08gJIlpY4mZ5znzJSsXQv06yfyrmpVUYwXLSp1VJrRuCAvUKCA8rmZmRkmTpyofP3mzRvtRMVMwtOnYmjD69cBLy9RjHt5SR0V+xTOccbyRtFdpUEDw/2pnPOcmYKsLOD774Fp08Trtm1Ff3FjvL+VVm5XkJqairlz56J48eLaWBwzAS9eiDtwXr4sfjo6dAjgzcN4cY4zpj7FBZ2G2n/8QzjPmTF5/Rro0uVtMT5xIvDPP8ZZjAMaFOSpqakICQlBjRo1ULduXWzduhUAEBoaiuLFi2PevHn48ssvdRUnMyKJiWLcz/PnAXd3UYyXKiV1VOxTOMcZy7vMTODECfHcEPuPc54zU/DoERAQAGzaBFhaAitXAjNmGPZdcT+J1DRhwgRycnKiDh06UKFChcjCwoIGDhxIFStWpL/++osyMjLUXZTRSUxMJACUmJgodSgGLymJqG5dIoCoYEGiS5ekjsh46Xu7y885TsR5zrTj/Hmx/7O3J0pP/3hbKba5/JznnOOmITKSqHBhkWeurkT//it1RB+myTandh/yjRs3YvXq1fj8889x+fJlVKpUCRkZGbhw4QJkhtpJjunV69dA69biSmdnZzHOeIUKUkfF1MU5zljeKfqP160LWGh8lZbucZ4zY7Z5M9Crl6g3fH2BHTuAEiWkjko71D65/+DBA1SvXh0AUKFCBcjlcnz55ZecwAyAGOarbVvg6FHA0RHYv19c6cyMB+c4Y3ln6P3HOc+ZMSISXVI6dBDFeHCwOPlnKsU4oEFBnpmZCSsrK+VrCwsL2Nvb6yQoZlxSU0WShIWJiyn27AFq1pQ6KqYpznHG8oZIdYQVQ8R5zoxNairQpw8waZJ4PXIksHMn4OQkaVhap/YPakSEPn36QC6XAwBSUlIwZMgQ2L13OevmzZu1GyEzaOnp4irn3bvFXbF27RI/1TLjwznOWN7ExIiLzSwsgFq1pI4mZ5znzJg8eQK0ayculDY3BxYuBIYOlToq3VC7IO/du7fK6549e2o9GGZcMjKAHj2AbdsAuVwMxB8QIHVULLc4xxnLG8XZ8erVAVtbaWP5EM5zZiyuXAFatRJ/6Do5ARs3iuGUTZXaBXloaKgu42BGJjNT/IS0caMYcmjzZnETIGa8OMcZyxtF/3FD7a4CcJ4z47B3L9C5M5CcLO52u3MnUK6c1FHpljGP2MgkkpUFDBwoblVrYSGK8hYtpI6KMcakpThDbqgXdDJm6IiAX34BWrYUxXjDhkBEhOkX4wAX5ExDRMDw4UBoqBiAf906oE0bqaNijDFpPXkCXL8unterJ20sjBmj9HRg2DBg9Ghx4q9fPzFYRMGCUkemHwY4SiozVETAmDHA0qWATAasWQN06iR1VIwxJj3F3TnLlwdcXaWNhTFj8+KF6KJy4ICoL2bNAsaNE8/zCy7ImVqIgK+/Fj8lAcDvvwPdu0sbE2OMGQpj6D/OmCG6eVPcVDAqSgydvG4d8PnnUkelf1yQM7V8/z0we7Z4vnQp0LevtPEwxpgh4f7jjGnuyBGgfXtxhtzLS9x5s3JlqaOShqR9yGfMmIGaNWvCwcEB7u7uaNu2LaKiolTapKSkYPjw4ShYsCDs7e3RoUMHxMXFqbS5d+8eWrZsCVtbW7i7u2P8+PHIyMhQaXPkyBFUq1YNcrkcpUqVwsqVK3W9eibjxx/FAxBnyAcPljYeZjw4x1l+8OoVcPaseJ7fzpBzjrPc+v13MYzhixdi3P5Tp/JvMQ4AIAkFBwdTaGgoXb58mc6fP08tWrQgb29vevnypbLNkCFDyMvLiw4ePEhnzpyhOnXqUN26dZXzMzIyqEKFChQUFETnzp2j3bt3k6urK4WEhCjb3Llzh2xtbWns2LF09epVWrhwIZmbm9PevXvVijMxMZEAUGJiovZW3kjMmkUkOqwQzZ4tdTT5iylsd8aS40Sm8X0zaRw8KPaRRYoQZWWp/z5T2OY4x5mmMjKIxo17W1t07Ur0+rXUUemGJtucpAX5++Lj4wkAHT16lIiIEhISyNLSkjZu3Khsc+3aNQJA4eHhRES0e/duMjMzo9jYWGWbJUuWkKOjI6WmphIR0YQJE8jPz0/ls7p06ULBwcFqxZVfk3j+/LcJ8+OPUkeT/5jidmeoOU5kmt8304+pU98WFpowxW2Oc5x9TFISUatWb2uLKVM0+yPW2GiyzRnUsIeJiYkAgAIFCgAAIiMjkZ6ejqB37jhTrlw5eHt7Izw8HAAQHh6OihUrwsPDQ9kmODgYSUlJuHLlirJN0Ht3rQkODlYu432pqalISkpSeeQ3ixaJEVUA4LvvgG++kTQcZiIMJccBznOmPdx//C3OcfYhJ04AtWuLm/xYWwN//QVMnpy/RlL5GIMpyLOysjBmzBjUq1cPFSpUAADExsbCysoKzs7OKm09PDwQGxurbPNuEivmK+Z9rE1SUhLevHmTLZYZM2bAyclJ+fDy8tLKOhqL+fOBESPE8wkTgKlTJQ2HmQhDynGA85xpx8uXwLFj4nlgoKShSI5znOUkIQEYOlT8wXrtGuDpKS7m7NpV6sgMi8EU5MOHD8fly5exfv16qUNBSEgIEhMTlY/79+9LHZLezJ4NfPmleB4SAvz0E//1yrTDkHIcyN95zrQnLAxISwNKlAB8faWORlqc4+xdROJO3uXLi9HZAHGzn8uXxZlypsoghj0cMWIEdu7ciWPHjqFo0aLK6Z6enkhLS0NCQoLKX9dxcXHw9PRUtjl16pTK8hRXb7/b5v0ruuPi4uDo6AgbG5ts8cjlcsjlcq2smzGZNg349lvxfPJk/imJaY+h5TiQf/Ocadf27eLf1q3z9/6Sc5y96949cVfvnTvF6zJlgGXL+Fekj5H0DDkRYcSIEdiyZQsOHTqE4sWLq8yvXr06LC0tcfDgQeW0qKgo3Lt3D/7+/gAAf39/XLp0CfHx8co2YWFhcHR0hO//n67w9/dXWYaijWIZ+R2RKL4VxfgPPwBTpuTvgwvTDs5xZsoyM4Fdu8Tz/HgjE4BznKnKyADmzRO/Fu3cCVhaiuvQLlzgYvyTdH2F6ccMHTqUnJyc6MiRI/T48WPl4/U7498MGTKEvL296dChQ3TmzBny9/cnf39/5XzFcEnNmjWj8+fP0969e8nNzS3H4ZLGjx9P165do0WLFvGwh/8vK4soJOTtFc8zZ0odEVMwhe3OWHKcyDS+b6ZfJ06I/aaTE1FamubvN4VtjnOcKURGElWr9raeqF+f6MoVqaOSltEMewggx0doaKiyzZs3b2jYsGHk4uJCtra21K5dO3r8+LHKcmJiYqh58+ZkY2NDrq6uNG7cOEpPT1dpc/jwYapSpQpZWVlRiRIlVD7jU0w1ibOyVMcCnTdP6ojYu0xhuzOWHCcyje+b6dfEibkb7lDBFLY5znGWnEw0diyRmdnbP1CXLyfKzJQ6Mulpss3JiIj0cSbemCUlJcHJyQmJiYlwdHSUOhytIBLDGv7yi3j966+ivxczHKa43Rky/r6Zpvz8gKtXgbVrge7dNX8/b3P6xd+39u3eDQwbBty9K1536SJGavv/rv/5nibbnEFc1Mn0KytLDGu4ZIl4vWwZMGiQtDExxpgxuX1bFOPm5kDz5lJHw5h+xcYCo0cDGzaI1z4+wOLFQIsW0sZlzAxm2EOmH1lZovheskRctPnHH1yMM8aYpnbsEP82aAC4uEgbC2P6kpUFLF8uhjLcsAEwMwPGjQOuXOFiPK/4DHk+kpkpxgBdvVok0apVQM+eUkfFGGPGR1GQ59fRVVj+c+2aOIF3/Lh4Xb26KM6rVZM2LlPBZ8jziYwM4IsvRDFubg6sW8fFOGOM5UZCwtu7c3JBzkxdSgrw/fdA5cqiGLezE0Mb/vcfF+PaxGfI84H0dKBHD3HHLAsLYP16oEMHqaNijDHjtHevOMlRvjxQsqTU0TCmO0eOAIMHAzduiNetWgGLFgHe3pKGZZL4DLmJS0sDOncWxbilJfDPP1yMM8ZYXijuzslnx5mpevZMdHFt1EgU456eoo7Yvp2LcV3hgtyEpaQA7dsDW7cCcrn4lw8gjDGWe+npwJ494nnr1tLGwpi2EYlhPMuXB0JDxbQhQ0T/8Y4d+Q7eusRdVkzUmzdAu3bAvn2AtTWwbRvQrJnUUTHGmHE7flz0IXd1BerUkToaxrTnzh1g6FBg/37x2tdXXLRZr560ceUXfIbcBL1+Lc7c7NsH2NoCu3ZxMc4YY9qgGF2lZUtxgTxjxi49HZg5E6hQQRTjcjnw44/AuXNcjOsTnyE3MS9fiosujh4F7O3FXbQaNJA6KsYYM35E3H+cmZaICDGU4cWL4nWjRuJmgaVLSxtXfsRnyE1IUhLw2WeiGHd0FGfIuRhnjDHtuHZN3KHTyop/dWTGLSkJGDkS8PcXxXjBgsDKlcDBg1yMS4XPkJuIhARRjEdEAM7OohivVUvqqBhjzHQouqs0bix+gWTMGG3dCowYATx8KF736gXMmQO4uUkaVr7HBbkJeP5cnK2JjAQKFADCwniwfsYY0zbursKM2YMH4qz41q3idcmSwNKlQFCQpGGx/8ddVozc06dAkyaiGHd1BQ4d4mKcMca07ckTIDxcPG/VStpYGNNEZibw669i1JStW8UNAkNCgEuXuBg3JHyG3IjFx4ti/PJlwMND9P3y85M6KsYYMz27domLOqtWBby8pI6GMfVcvCgu2oyIEK/r1BFDGVasKG1cLDs+Q26kHj8GAgNFMV6okLi9LRfjjDGmG4r+43wzIGYMXr8GJk4EqlcXxbijo7jl/YkTXIwbKj5DboQePBAXFd28CRQtKrqp8FXRjDGmGykp4kJ5gPuPM8O3f7+4u2Z0tHjdvj3wyy9AkSLSxsU+js+QG5l794CAAFGMe3uLIQ65GGeMMd05cgR49QooXJiv0WGGKz4e6NkTCA4WxXjRouIu3f/8w8W4MeCC3IjcuSOK8Tt3gOLFgWPHgBIlpI6KMcZMm2J0ldatAZlM2lgYe19KCjBvHlC+PLB2rdhGR40Crl7lX3SMCXdZMRLXromroR89AkqVAg4fFn/9MsYY062DB8W/LVpIGwdj78rMBP78E/j+e/HrOQBUrgysWAHUrCltbExzfIbcCJw/L86MP3okLtw8doyLccYY04e4OODGDfG8fn1pY2EMEKP97NgBVKkC9OkjivEiRUQhfuYMF+PGis+QG7iICHEHzoQE0Xdx3z4x3jhjjDHdO35c/FuxorjxGmNSOnEC+Ppr8S8g7swdEiJu+GNjI2loLI+4IDdgR4+KG1C8fAnUrQvs3g04OUkdFWOM5R///iv+bdBA2jhY/nblCjBp0tvrGaytgdGjRXHu4iJtbEw7uCA3UHv3Au3aiYs1GjcWV0rb20sdFWOM5S9ckDMp3bsHTJ4MrF4NZGUB5uZAv35iGo+cYlq4IDdAW7YAXboA6elAy5bApk3ir2HGGGP6k5QkruEBuCBn+vXsGTB9uriZT2qqmNahA/Djj0C5ctLGxnSDC3IDs3Yt0Lu3uHq6UydxBbWVldRRMcZY/hMeLs5KFi/OZyOZfrx6BcyfD8yaJf4gBMRduX/6CahdW8rImK5xQW5Ali8Xd9ciEkX5b78BFvw/xBhjkjh2TPzLZ8eZrqWnA7//DkydCsTGimmVK4tCPDiYx7/PD3jYQwMxbx4weLAoxocNA/74g4txxhiTEvcfZ7qWlQVs2CCGNB46VBTjxYuLX8vPnhWjrHExnj9IWpAfO3YMrVu3RuHChSGTybB161aV+USE77//HoUKFYKNjQ2CgoJw8+ZNlTbPnz9Hjx494OjoCGdnZ/Tv3x8vX75UaXPx4kU0aNAA1tbW8PLywqxZs3S9amojEn3Cxo4Vr8ePB379FTDjP5WYCeAcZ8YqNRU4dUo854L84zjPc+fAAaBWLXHN2M2bgJsbsHAhcP060L071wH5jaT/3a9evULlypWxaNGiHOfPmjULv/zyC5YuXYqIiAjY2dkhODgYKSkpyjY9evTAlStXEBYWhp07d+LYsWMYNGiQcn5SUhKaNWsGHx8fREZGYvbs2ZgyZQqWL1+u8/X7FCIxjNF334nXU6cCM2fyX8PMdOT3HGfG6/RpUZS7uwNlykgdjWHjPNdMZCTQtKl4REaKEdSmTAFu3wZGjODrxvItMhAAaMuWLcrXWVlZ5OnpSbNnz1ZOS0hIILlcTn/99RcREV29epUA0OnTp5Vt9uzZQzKZjB4+fEhERIsXLyYXFxdKTU1Vtvn666+pbNmyaseWmJhIACgxMTG3q5dNZibRiBFEoiwn+vlnrS2amQhdbHdSMuQcJzK975vlzfTpYt/coYPuPsMUtzlDznOpv+8bN4g6d3573Le0JBo1iiguTpJwmB5oss0Z7A8i0dHRiI2NRVBQkHKak5MTateujfDwcABAeHg4nJ2dUaNGDWWboKAgmJmZISIiQtmmYcOGsHrnT87g4GBERUXhxYsXOX52amoqkpKSVB7alJkJDBgguqbIZMDSpcC4cVr9CMYMnpQ5Dug+z5lx4/7j2mHKx3J1xcaKa8N8fUV/cZkM6NkTiIoCFiwQv8IwZrAFeez/X2bs4eGhMt3Dw0M5LzY2Fu7vbckWFhYoUKCASpuclvHuZ7xvxowZcHJyUj68vLzyvkL/Lz0d6NEDCA0V/cNWrRIXczKW30iZ44Bu85wZt8zMt7cm54I8b0z1WK6OxETg22+BkiWBJUuAjAygeXPg3DlgzRpx8SZjCgZbkEspJCQEiYmJysf9+/e1styUFDGw/99/A5aW4i/lXr20smjGmIZ0lefM+F28KMaAdnAQQ88x4yRVjqekAHPnikJ82jTg9WsxhviRI8Du3bxNsZwZ7MB6np6eAIC4uDgUKlRIOT0uLg5VqlRRtomPj1d5X0ZGBp4/f658v6enJ+Li4lTaKF4r2rxPLpdDLpdrZT0UXr0C2rYVV1VbWwP//AO0aKHVj2DMqEiZ44Bu8pyZBkV3lbp1xa3KWe6Z2rH8YzIzxZnvyZPFLe8BcVfN6dPF8Z8HbGAfY7BnyIsXLw5PT08cPHhQOS0pKQkRERHw9/cHAPj7+yMhIQGRkZHKNocOHUJWVhZq//8trfz9/XHs2DGkp6cr24SFhaFs2bJwcXHRy7okJoqB/Q8cAOzsxF/IXIyz/M6UcpyZFu4/rj35Ic+JgO3bxZnvvn1FMV6kiLi536VLQLt2XIwzNejhItMPSk5OpnPnztG5c+cIAM2dO5fOnTtHd+/eJSKin376iZydnWnbtm108eJFatOmDRUvXpzevHmjXMZnn31GVatWpYiICDp+/DiVLl2aunXrppyfkJBAHh4e1KtXL7p8+TKtX7+ebG1tadmyZWrHmZcrs58+JapRQ1xR7eREdPKkxotg+ZTUIwJog7HkOJFpfN8s77KyiDw8xD772DHdfpapbHPGkue6+L6PHyeqV+/tyCkuLkSzZhG9fq21j2BGTJNtTtKC/PDhwwQg26N3795EJIZL+u6778jDw4Pkcjk1adKEoqKiVJbx7Nkz6tatG9nb25OjoyP17duXkpOTVdpcuHCB6tevT3K5nIoUKUI//fSTRnHmNokfPyaqUEEkqasr0dmzGr2d5XOmcLA2lhwnMo3vm+XdjRtin21lRfROvagTprLNGUuea/P7vnSJqHXrt4W4tTXRxIlEz5/nedHMhGiyzcmIiPRxJt6YJSUlwcnJCYmJiXB0dFTrPffvA02aiLtvFSokuqv4+uo4UGZScrPdsdzj75sBwB9/AP37A/Xrv+26oiu8zemXNr7vu3dFH/HVq0Upbm4O9OsnphUpouWAmdHTZJsz2Is6jdnt26IYv3sX8PEBDh4UV1szxhgzbNx/nOXk6VNxceaiRUBampjWoYMYRaVsWWljY6aBC3Itu3oVCAoCHj8GSpcWxTgPb8wYY8bh2DHxLxfkDBAjpM2fD8yaJYbCBIBGjYCffgJq1ZI0NGZiuCDXonPngGbNxF/SFSoAYWHAR0ZdY4wxZkAePQLu3BEjYtStK3U0TErp6WKUlP/9T9xpEwCqVBGFeLNmPGoK0z4uyLUkPFzcgSsxEaheHdi3DyhYUOqoGGOMqUvRXaVyZcDJSdpYmDSysoCNG8UdNm/dEtNKlAB+/BHo0kXcYZsxXeCCXAuOHAFatRI/bdWvD+zcyTtzxhgzNtx/PH87cACYOBFQDIfu7g589x0waBBgZSVtbMz0cUGeR7t3iws7UlJE3/GtW8XNfxhjjBkXRUHesKG0cTD9io4WRfeBA+K1vT0wfjwwdqx4zpg+cEGeR3v3imK8dWtgwwbA2lrqiBhjjGkqMVHcVRHgM+T5jaMjcOoUYGkJDBsGfPMN4OYmdVQsv+GCPI/mzxfji/fvL5KZMcaY8XFyErc8P3MG8PCQOhqmTwULAn/+CVSsCBQrJnU0LL/igjyPzMyAIUOkjoIxxlheFS0qHiz/ad1a6ghYfsfXCzPGGGOMMSYhLsgZY4wxxhiTEBfkjDHGGGOMSYgLcsYYY4wxxiTEBTljjDHGGGMS4oKcMcYYY4wxCXFBzhhjjDHGmIR4HHI1EBEAICkpSeJIWH6i2N4U2x/TLc5zpm+c4/rFOc70TZMc54JcDcnJyQAALy8viSNh+VFycjKcnJykDsPkcZ4zqXCO6wfnOJOKOjkuI/7T/JOysrLw6NEjODg4QCaTSR2OUlJSEry8vHD//n04OjpKHY5O5Zd1fXc9HRwckJycjMKFC8PMjHuX6drH8jy/bH+6wN/dhxER57geGeqx/EM4d/RDl9+zJjnOZ8jVYGZmhqIGfD9lR0fHfJOs+WVdFevJZ830R508zy/bny7wd5czznH9MfRj+Ydw7uiHrr5ndXOc/yRnjDHGGGNMQlyQM8YYY4wxJiEuyI2YXC7H5MmTIZfLpQ5F5/LLuuaX9TQ2/P+Se/zdMZY7nDv6YSjfM1/UyRhjjDHGmIT4DDljjDHGGGMS4oKcMcYYY4wxCXFBzhhjjDHGmIS4IGeMMcYYY0xCXJAbmEWLFqFYsWKwtrZG7dq1cerUqY+237hxI8qVKwdra2tUrFgRu3fvVpnfp08fyGQylcdnn32my1VQiybreeXKFXTo0AHFihWDTCbD/Pnz87xMfdL2uk6ZMiXb/2m5cuV0uAb5m6FuV/oyY8YM1KxZEw4ODnB3d0fbtm0RFRWl0iYlJQXDhw9HwYIFYW9vjw4dOiAuLk6lzb1799CyZUvY2trC3d0d48ePR0ZGhkqbI0eOoFq1apDL5ShVqhRWrlyp69VjTK+0eYxPT0/H119/jYoVK8LOzg6FCxfGF198gUePHul6NQyetmupdw0ZMuSjtUiuETMY69evJysrK/rjjz/oypUrNHDgQHJ2dqa4uLgc2584cYLMzc1p1qxZdPXqVfr222/J0tKSLl26pGzTu3dv+uyzz+jx48fKx/Pnz/W1SjnSdD1PnTpFX331Ff3111/k6elJ8+bNy/My9UUX6zp58mTy8/NT+T998uSJjtckfzLU7UqfgoODKTQ0lC5fvkznz5+nFi1akLe3N718+VLZZsiQIeTl5UUHDx6kM2fOUJ06dahu3brK+RkZGVShQgUKCgqic+fO0e7du8nV1ZVCQkKUbe7cuUO2trY0duxYunr1Ki1cuJDMzc1p7969el1fxnRF28f4hIQECgoKor///puuX79O4eHhVKtWLapevbo+V8vg6KKWUti8eTNVrlyZChcunOPxOS+4IDcgtWrVouHDhytfZ2ZmUuHChWnGjBk5tu/cuTO1bNlSZVrt2rVp8ODByte9e/emNm3a6CTe3NJ0Pd/l4+OTYxLkZZm6pIt1nTx5MlWuXFmLUbIPMdTtSkrx8fEEgI4ePUpEoiiwtLSkjRs3Kttcu3aNAFB4eDgREe3evZvMzMwoNjZW2WbJkiXk6OhIqampREQ0YcIE8vPzU/msLl26UHBwsK5XiTG90MUx/n2nTp0iAHT37l3tBG2EdPU9P3jwgIoUKUKXL1/+4PE5L7jLioFIS0tDZGQkgoKClNPMzMwQFBSE8PDwHN8THh6u0h4AgoODs7U/cuQI3N3dUbZsWQwdOhTPnj3T/gqoKTfrKcUytUGXcd28eROFCxdGiRIl0KNHD9y7dy+v4bL3GOp2JbXExEQAQIECBQAAkZGRSE9PV/meypUrB29vb+X3FB4ejooVK8LDw0PZJjg4GElJSbhy5YqyjTr7M8aMkS6P8e9KTEyETCaDs7OzVuI2Nrr6nrOystCrVy+MHz8efn5+OomdC3ID8fTpU2RmZqocsADAw8MDsbGxOb4nNjb2k+0/++wzrF69GgcPHsTMmTNx9OhRNG/eHJmZmdpfCTXkZj2lWKY26Cqu2rVrY+XKldi7dy+WLFmC6OhoNGjQAMnJyXkNmb3DULcrKWVlZWHMmDGoV68eKlSoAEDsh6ysrLIVAO9+Tx/aVynmfaxNUlIS3rx5o4vVYUxvdHWMf1dKSgq+/vprdOvWDY6OjtoJ3Mjo6nueOXMmLCwsMGrUKO0H/f8sdLZkZhC6du2qfF6xYkVUqlQJJUuWxJEjR9CkSRMJI2O51bx5c+XzSpUqoXbt2vDx8cGGDRvQv39/CSNjpm748OG4fPkyjh8/LnUojLF3pKeno3PnziAiLFmyROpwTEpkZCQWLFiAs2fPQiaT6exz+Ay5gXB1dYW5uXm2kQni4uLg6emZ43s8PT01ag8AJUqUgKurK27dupX3oHMhN+spxTK1QV9xOTs7o0yZMpL9n5oqQ92upDJixAjs3LkThw8fRtGiRZXTPT09kZaWhoSEBJX2735PH9pXKeZ9rI2joyNsbGy0vTqM6ZUuj/GKYvzu3bsICwvLt2fHAd18z//++y/i4+Ph7e0NCwsLWFhY4O7duxg3bhyKFSumtdi5IDcQVlZWqF69Og4ePKiclpWVhYMHD8Lf3z/H9/j7+6u0B4CwsLAPtgeABw8e4NmzZyhUqJB2AtdQbtZTimVqg77ievnyJW7fvi3Z/6mpMtTtSt+ICCNGjMCWLVtw6NAhFC9eXGV+9erVYWlpqfI9RUVF4d69e8rvyd/fH5cuXUJ8fLyyjaJw8PX1VbbRdH/GmLHQ1TFeUYzfvHkTBw4cQMGCBXWzAkZCF99zr169cPHiRZw/f175KFy4MMaPH499+/ZpL3itXiLK8mT9+vUkl8tp5cqVdPXqVRo0aBA5OzsrRybo1asXTZw4Udn+xIkTZGFhQT///DNdu3aNJk+erDJUT3JyMn311VcUHh5O0dHRdODAAapWrRqVLl2aUlJSJFlHIs3XMzU1lc6dO0fnzp2jQoUK0VdffUXnzp2jmzdvqr1MqehiXceNG0dHjhyh6OhoOnHiBAUFBZGrqyvFx8frff1MnaFuV/o0dOhQcnJyoiNHjqgMtfn69WtlmyFDhpC3tzcdOnSIzpw5Q/7+/uTv76+crxj2sFmzZnT+/Hnau3cvubm55Tjs4fjx4+natWu0aNEiHvaQmRRtH+PT0tLo888/p6JFi9L58+dV8lMxelF+pO3vOSe6GGWFC3IDs3DhQvL29iYrKyuqVasW/ffff8p5AQEB1Lt3b5X2GzZsoDJlypCVlRX5+fnRrl27lPNev35NzZo1Izc3N7K0tCQfHx8aOHCgQRQTmqxndHQ0Acj2CAgIUHuZUtL2unbp0oUKFSpEVlZWVKRIEerSpQvdunVLj2uUvxjqdqUvOW2PACg0NFTZ5s2bNzRs2DBycXEhW1tbateuHT1+/FhlOTExMdS8eXOysbEhV1dXGjduHKWnp6u0OXz4MFWpUoWsrKyoRIkSKp/BmCnQ5jH+Q8cLAHT48GE9rZFh0ub3nBNdFOQyIiLtnW9njDHGGGOMaYL7kDPGGGOMMSYhLsgZY4wxxhiTEBfkjDHGGGOMSYgLcsYYY4wxxiTEBTljjDHGGGMS4oKcMcYYY4wxCXFBzhhjjDHGmIS4IGeMMcYYY0xCXJAzxhhjjDEmIS7ImU7cv38fgYGB8PX1RaVKlbBx40apQ2KMmYiEhATUqFEDVapUQYUKFbBixQqpQ2LM5PBxXL9kRERSB8FMz+PHjxEXF4cqVaogNjYW1atXx40bN2BnZyd1aIwxI5eZmYnU1FTY2tri1atXqFChAs6cOYOCBQtKHRpjJoOP4/rFZ8iZWgIDAyGTySCTyXD+/PlPti9UqBCqVKkCAPD09ISrqyueP3+u2yBz0KdPH2XcW7du1fvnM2aIcsrnwMBAjBkzRtK41GVubg5bW1sAQGpqKogI755b4rxnLDs+jhs2LsiZ2gYOHIjHjx+jQoUKGr0vMjISmZmZ8PLy+mTbvn374ttvv81tiNksWLAAjx8/1tryGDMVuc3n3NJ2bickJKBy5cooWrQoxo8fD1dXV+U8znvGcsbHccPFBTlTm62tLTw9PWFhYaH2e54/f44vvvgCy5cv/2TbzMxM7Ny5E59//nlewlTh5OQET09PrS2PMVORm3zOLV3ktrOzMy5cuIDo6GisW7cOcXFxynmc94zljI/jhosLchNWtGhRLF68WGXayZMnYWtri7t37+Z5+Zs2bULFihVhY2ODggULIigoCK9evVLOT01NRdu2bTFx4kTUrVv3k8s7efIkLC0tUbNmzRznBwYGYuTIkRgzZgxcXFzg4eGBFStW4NWrV+jbty8cHBxQqlQp7NmzJ8/rxpih0XU+v2/Xrl1wcnLC2rVrAQDJycno0aMH7OzsUKhQIcybN0/tbi66zG0PDw9UrlwZ//77b57WlzFDxMfx/IMLchNWu3ZtnD59WvmaiDBmzBh8+eWX8PHxydOyHz9+jG7duqFfv364du0ajhw5gvbt2yv7cRIR+vTpg8aNG6NXr15qLXP79u1o3bo1ZDLZB9usWrUKrq6uOHXqFEaOHImhQ4eiU6dOqFu3Ls6ePYtmzZqhV69eeP36dZ7WjzFDo8t8ft+6devQrVs3rF27Fj169AAAjB07FidOnMD27dsRFhaGf//9F2fPnlVredrO7bi4OCQnJwMAEhMTcezYMZQtWzbvK86YgeHjeD5CzGTNmjWL/Pz8lK9XrVpFnp6elJycrPGyAgICaPTo0crXkZGRBIBiYmJybP/vv/+STCajypUrKx8XL1786GeULl2adu7c+dEY6tevr3ydkZFBdnZ21KtXL+W0x48fEwAKDw9XeS8A2rJly0c/nzFDpst8fnfar7/+Sk5OTnTkyBHlvKSkJLK0tKSNGzcqpyUkJJCtrW225eRE27kdERFBlStXpkqVKlHFihVp6dKlOS6X854ZOz6Ov2Xq+az7zoNMMnXq1MHEiRPx8uVLyGQyTJo0CT/++CPs7e0BAO3atcORI0fQpEkTbNq0SaNlV65cGU2aNEHFihURHByMZs2aoWPHjnBxcQEA1K9fH1lZWWov79q1a3j06BGaNGny0XaVKlVSPjc3N0fBggVRsWJF5TQPDw8AQHx8vCarw5jB+1g+379/H7169UJ8fDwsLCzw3XffoVOnThp/xqZNmxAfH48TJ06o/OR8584dpKeno1atWsppTk5Oap2V1kVuf/7552qNEsGYsftY3ickJCAoKAgZGRnIyMjA6NGjMXDgQLWXzcdxw8JdVkxY9erVYWZmhrNnz2LmzJlwc3ND3759lfNHjx6N1atX52rZ5ubmCAsLw549e+Dr64uFCxeibNmyiI6OztXytm/fjqZNm8La2vqj7SwtLVVey2QylWmKn8k02YkwZgw+ls8WFhaYP38+rl69iv3792PMmDEq/UDVVbVqVbi5ueGPP/5QGUYwLzi3Gcu9j+W9g4MDjh07hvPnzyMiIgLTp0/Hs2fP1F42H8cNCxfkJszW1hYVK1bEP//8g59//hnz5s2Dmdnb//LAwEA4ODjkevkymQz16tXD1KlTce7cOVhZWWHLli25Wta2bdvQpk2bXMfCmKn7WD5ra7zgkiVL4vDhw9i2bRtGjhypnF6iRAlYWlqq9GVNTEzEjRs3PrlMzm3Gcu9jef+p8fjVwcdxw8FdVkxcnTp1sHDhQrRp0waBgYFaW25ERAQOHjyIZs2awd3dHREREXjy5AnKly+v8bLi4+Nx5swZbN++XWvxMWaK1MlnTcYLzkmZMmVw+PBhBAYGKs+8Ozg4oHfv3hg/fjwKFCgAd3d3TJ48GWZmZh+9eItzm7G8+1jeJyQkICAgADdv3sTs2bNVxuP/FD6OGxYuyE1c5cqVYWlpidmzZ2t1uY6Ojjh27Bjmz5+PpKQk+Pj4YM6cOWjevLnGy9qxYwdq1aql0Y6EsfzoU/msGC94xYoVefqcsmXL4tChQwgMDIS5uTnmzJmDuXPnYsiQIWjVqhUcHR0xYcIE3L9//6M/T3NuM5Z3H8t7xXj8cXFxaN++PTp27Kjsg/0pfBw3LFyQm7j169djxIgRKFWqlFaXW758eezdu1cry9q2bZtaNxE4cuRItmkxMTHZpmmr7ytjhuZj+azpeMHvez+/ypcvr3KzHQcHB+WY5ADw6tUrTJ06FYMGDfrgMjm3Gcs7dY7j747H37FjR7WWy8dxw8J9yE1QVlYW4uLiMH36dNy8eROTJ0/WynIXL14Me3t7XLp0SSvLU6hfvz66deum1WUqDBkyRDmqDGPGSJ18plyMF6xpPp87dw5//fUXbt++jbNnzyrHJ/9Yn1Fd5vbHcN4zY6dO3udmPH4+jhsuGeXHP0NM3JEjR9C4cWOUK1cOoaGhqF27do7tgoKCcOHCBbx69QoFChTAxo0b4e/vn2Pbhw8f4s2bNwAAb29vWFlZ6Sx+bYqPj0dSUhIAceGbnZ2dxBExphl18vn48eNo2LChynBia9asURlK7F25yedz585hwIABiIqKgpWVFapXr465c+d+8DOkxHnPjJ06eX/q1CkMGjRIeTHn8OHDMXjw4A8uk4/jho0LcsYYY4wxxiTEXVYYY4wxxhiTEBfkjDHGGGOMSYgLcsYYY4wxxiTEBTljjDHGGGMS4oKcMcYYY4wxCXFBzhhjjDHGmIS4IGeMMcYYY0xCXJAzxhhjjDEmIS7IGWOMMcYYkxAX5IwxxhhjjEmIC3LGGGOMMcYkxAU5Y4wxxhhjEuKCnDHGGGOMMQlxQc4YY4wxxpiEuCBnjDHGGGNMQlyQM8YYY4wxJiEuyBljjDHGGJMQF+SMMcYYY4xJiAtyxhhjOiWTyTBlyhSN33fkyBHIZDIcOXJE6zEZKplMhhEjRkgdBstHOD8NAxfkBkYmk6n1MIYEiImJQd++fVGyZElYW1vD09MTDRs2xOTJk1XaLV68GCtXrpQmSGZyoqOjMWLECJQpUwa2trawtbWFr68vhg8fjosXL0od3ketXLlSJc+tra1RpkwZjBgxAnFxcRovT5+5tXv37lwd1LVB29/bx/D+Kv/i/My9S5cuoWPHjvDx8YG1tTWKFCmCpk2bYuHChSrtpk+fjq1bt0oTpMQspA6AqVqzZo3K69WrVyMsLCzb9PLly+szLI3dunULNWvWhI2NDfr164dixYrh8ePHOHv2LGbOnImpU6cq2y5evBiurq7o06ePdAEzk7Bz50506dIFFhYW6NGjBypXrgwzMzNcv34dmzdvxpIlSxAdHQ0fHx+pQ/2o//3vfyhevDhSUlJw/PhxLFmyBLt378bly5dha2ur9nL0mVu7d+/GokWLcjzov3nzBhYWuj/caOt7+xjeXzHOT82cPHkSjRo1gre3NwYOHAhPT0/cv38f//33HxYsWICRI0cq206fPh0dO3ZE27ZtdRaPoeKC3MD07NlT5fV///2HsLCwbNMN3bx58/Dy5UucP38+W/ETHx+f6+W+evUKdnZ2eQ2PmaDbt2+ja9eu8PHxwcGDB1GoUCGV+TNnzsTixYthZvbxHwYNYRtr3rw5atSoAQAYMGAAChYsiLlz52Lbtm3o1q2bpLHlhrW1tV4+R5ff2+vXr7VW1L+PiJCSkgIbGxudLJ9pF+enZqZNmwYnJyecPn0azs7OKvO4HniLu6wYodDQUDRu3Bju7u6Qy+Xw9fXFkiVLsrUrVqwYWrVqhePHj6NWrVqwtrZGiRIlsHr16mxtL168iICAANjY2KBo0aL48ccfERoaCplMhpiYGI1jvH37NooWLZrjmUh3d3eVGK9cuYKjR48qfwYMDAwE8PbnwaNHj2LYsGFwd3dH0aJFle/ds2cPGjRoADs7Ozg4OKBly5a4cuWKymfFxsaib9++KFq0KORyOQoVKoQ2bdqorNOZM2cQHBwMV1dX2NjYoHjx4ujXr5/G68ykNWvWLLx69QqhoaHZinEAsLCwwKhRo+Dl5aWc1qdPH9jb2+P27dto0aIFHBwc0KNHDwDAv//+i06dOsHb2xtyuRxeXl748ssv8ebNG+X7f/75Z8hkMty9ezfb54WEhMDKygovXrzI87o1btwYgOiOAwAZGRn44YcfULJkScjlchQrVgyTJk1Camqq8j0fyy0ASEhIwJgxY+Dl5QW5XI5SpUph5syZyMrKUraJiYmBTCbDzz//jOXLlys/r2bNmjh9+rTK97ho0SIAqt3uFN7vo3r37l0MGzYMZcuWhY2NDQoWLIhOnTrlal/zMe9/bwDw559/onr16rCxsUGBAgXQtWtX3L9/X+V9gYGBqFChAiIjI9GwYUPY2tpi0qRJH/1Op0yZorLOCor92Lvrptg379u3DzVq1ICNjQ2WLVum8r61a9eibNmysLa2RvXq1XHs2DEtfStM2zg/P+727dvw8/PLVowDqvWATCbDq1evsGrVKmWMil8PFPl19epVdO/eHS4uLqhfv77yverk9c2bN9GhQwd4enrC2toaRYsWRdeuXZGYmKhsExYWhvr168PZ2Rn29vYoW7YsJk2alKv11hSfITdCS5YsgZ+fHz7//HNYWFhgx44dGDZsGLKysjB8+HCVtrdu3ULHjh3Rv39/9O7dG3/88Qf69OmD6tWrw8/PDwDw8OFDNGrUCDKZDCEhIbCzs8Nvv/0GuVye6xh9fHxw4MABHDp0SLmzysn8+fMxcuRI2Nvb45tvvgEAeHh4qLQZNmwY3Nzc8P333+PVq1cARNee3r17Izg4GDNnzsTr16+xZMkS1K9fH+fOnUOxYsUAAB06dMCVK1cwcuRIFCtWDPHx8QgLC8O9e/eUr5s1awY3NzdMnDgRzs7OiImJwebNm3O97kwaO3fuRKlSpVC7dm2N3peRkYHg4GDUr18fP//8s/Is6MaNG/H69WsMHToUBQsWxKlTp7Bw4UI8ePAAGzduBAB07twZEyZMwIYNGzB+/HiV5W7YsAHNmjWDi4tLntft9u3bAICCBQsCEGflVq1ahY4dO2LcuHGIiIjAjBkzcO3aNWzZsgXAx3Pr9evXCAgIwMOHDzF48GB4e3vj5MmTCAkJwePHjzF//nyVz1+3bh2Sk5MxePBgyGQyzJo1C+3bt8edO3dgaWmJwYMH49GjRzl2r8vJ6dOncfLkSXTt2hVFixZFTEwMlixZgsDAQFy9elVrZ6Lf/96mTZuG7777Dp07d8aAAQPw5MkTLFy4EA0bNsS5c+dUCoZnz56hefPm6Nq1K3r27AkPDw8EBgZ+cn+lrqioKHTr1g2DBw/GwIEDUbZsWeW8o0eP4u+//8aoUaMgl8uxePFifPbZZzh16hQqVKiQy2+D6Qrn58f5+PggPDwcly9f/uj2u2bNGgwYMAC1atXCoEGDAAAlS5ZUadOpUyeULl0a06dPBxEBUC+v09LSEBwcjNTUVIwcORKenp54+PAhdu7ciYSEBDg5OeHKlSto1aoVKlWqhP/973+Qy+W4desWTpw4odH65hoxgzZ8+HB6/7/p9evX2doFBwdTiRIlVKb5+PgQADp27JhyWnx8PMnlcho3bpxy2siRI0kmk9G5c+eU0549e0YFChQgABQdHa1x3JcvXyYbGxsCQFWqVKHRo0fT1q1b6dWrV9na+vn5UUBAQLbpoaGhBIDq169PGRkZyunJycnk7OxMAwcOVGkfGxtLTk5OyukvXrwgADR79uwPxrllyxYCQKdPn9Z4HZnhSExMJADUtm3bbPNevHhBT548UT7ezZ/evXsTAJo4cWK29+WUZzNmzCCZTEZ3795VTvP396fq1aurtDt16hQBoNWrV2u0Hopt/sCBA/TkyRO6f/8+rV+/ngoWLEg2Njb04MEDOn/+PAGgAQMGqLz3q6++IgB06NAh5bQP5dYPP/xAdnZ2dOPGDZXpEydOJHNzc7p37x4REUVHRxMAKliwID1//lzZbtu2bQSAduzYoZyW075KAQBNnjxZ+Tqn7zY8PDzbd3b48GECQIcPH85xuQrqfG8xMTFkbm5O06ZNU3nvpUuXyMLCQmV6QEAAAaClS5dm+6wPfaeTJ0/Ocf0Vsb27H1Xsm/fu3ZutPQACQGfOnFFOu3v3LllbW1O7du0++j0w3eL8zF1+7t+/n8zNzcnc3Jz8/f1pwoQJtG/fPkpLS8vW1s7Ojnr37p1tuiK/unXrpjJd3bw+d+4cAaCNGzd+MM558+YRAHry5MlH10dXuMuKEXq3n2FiYiKePn2KgIAA3LlzR+WnFwDw9fVFgwYNlK/d3NxQtmxZ3LlzRzlt79698Pf3R5UqVZTTChQooPzpPjf8/Pxw/vx59OzZEzExMViwYAHatm0LDw8PrFixQqNlDRw4EObm5srXYWFhSEhIQLdu3fD06VPlw9zcHLVr18bhw4cBiO/JysoKR44c+WC3AcUZsZ07dyI9PT13K8skl5SUBACwt7fPNi8wMBBubm7Kh+Kn23cNHTo027R38+zVq1d4+vQp6tatCyLCuXPnlPO6dOmCyMhI5VkyAPj7778hl8vRpk2bXK1PUFAQ3Nzc4OXlha5du8Le3h5btmxBkSJFsHv3bgDA2LFjVd4zbtw4AMCuXbs+ufyNGzeiQYMGcHFxUcmhoKAgZGZmZuse0aVLF5Uz/Yp9yrv7EU28+92mp6fj2bNnKFWqFJydnXH27NlcLRP4+Pe2efNmZGVloXPnzirr7OnpidKlSyv3GwpyuRx9+/bNdSyfUrx4cQQHB+c4z9/fH9WrV1e+9vb2Rps2bbBv3z5kZmbqLCamHs5PzTRt2hTh4eH4/PPPceHCBcyaNQvBwcEoUqQItm/frtGyhgwZovJa3bx2cnICAOzbtw+vX7/OcdmKemDbtm0qXYP0hbusGKETJ05g8uTJCA8Pz7ZhJSYmKjc8QOzI3+fi4qJSoN69exf+/v7Z2pUqVSpPcZYpUwZr1qxBZmYmrl69ip07d2LWrFkYNGgQihcvjqCgILWWU7x4cZXXN2/eBIAPdoVxdHQEIA6oM2fOxLhx4+Dh4YE6deqgVatW+OKLL+Dp6QkACAgIQIcOHTB16lTMmzcPgYGBaNu2Lbp3756nLjtMvxwcHAAAL1++zDZv2bJlSE5ORlxcXI4XR1tYWKhcm6Bw7949fP/999i+fXu2P+je/cO3U6dOGDt2LP7++29MmjQJRISNGzeiefPmym1RU4sWLUKZMmVgYWEBDw8PlC1bVnkx6t27d2FmZpYtPz09PeHs7Jxjf/b33bx5ExcvXoSbm1uO89+/0Or9/Yji4J/b/vFv3rzBjBkzEBoaiocPHyp/egaQ7aSCJj72vd28eRNEhNKlS+f4XktLS5XXRYoUgZWVVa5j+ZT392vvyinGMmXK4PXr13jy5Ily/8WkwfmpuZo1a2Lz5s1IS0vDhQsXsGXLFsybNw8dO3bE+fPn4evrq9ZycqoH1Mnr4sWLY+zYsZg7dy7Wrl2LBg0a4PPPP0fPnj2VNVOXLl3w22+/YcCAAZg4cSKaNGmC9u3bo2PHjp8cDEAbuCA3Mrdv30aTJk1Qrlw5zJ07F15eXrCyssLu3bsxb968bH/VvXtm+V3vJpiumZubo2LFiqhYsSL8/f3RqFEjrF27Vu2C/P2RBxTruGbNmhwPTO8O3zRmzBi0bt0aW7duxb59+/Ddd99hxowZOHToEKpWrQqZTIZNmzbhv//+w44dO7Bv3z7069cPc+bMwX///ZfjGVdmeJycnFCoUCFcvnw52zxFn/IPXZAkl8uz7WwzMzPRtGlTPH/+HF9//TXKlSsHOzs7PHz4EH369FHJs8KFC6NBgwbYsGEDJk2ahP/++w/37t3DzJkzc70+tWrVUo7i8CE5XUCorqysLDRt2hQTJkzIcX6ZMmVUXmt7PzJy5EiEhoZizJgx8Pf3h5OTE2QyGbp27ZqnM1Mf+96ysrIgk8mwZ8+eHNfn/VzXdMSTD/1/fOiMNo+oYrw4P3PPysoKNWvWRM2aNVGmTBn07dsXGzduzHZ/kg/JqR5QN6/nzJmDPn36YNu2bdi/fz9GjRqFGTNm4L///kPRokVhY2ODY8eO4fDhw9i1axf27t2Lv//+G40bN8b+/fs/+D1rCxfkRmbHjh1ITU3F9u3bVf4qfv/nVk34+Pjg1q1b2abnNC2vFDuxx48fK6dpuuNSXOTh7u6uVlFfsmRJjBs3DuPGjcPNmzdRpUoVzJkzB3/++aeyTZ06dVCnTh1MmzYN69atQ48ePbB+/XoMGDBAo9iYdFq2bInffvsNp06dQq1atfK0rEuXLuHGjRtYtWoVvvjiC+X0sLCwHNt36dIFw4YNQ1RUFP7++2/Y2tqidevWeYrhQ3x8fJCVlYWbN2+q3I8gLi4OCQkJKiMbfSi3SpYsiZcvX6r9R7E6NMnjTZs2oXfv3pgzZ45yWkpKChISErQWz/tKliwJIkLx4sWzFTSa+NB6Ks5KJiQkqFwcqs4Z0fcpfgV8140bN2Bra/vBs6bMMHB+qk9b9YAmea04Ofjtt9/i5MmTqFevHpYuXYoff/wRAGBmZoYmTZqgSZMmmDt3LqZPn45vvvkGhw8f1ur/R064D7mRUfyF9v5PSKGhobleZnBwMMLDw3H+/HnltOfPn2Pt2rXZ2j5+/BjXr1//ZH/rf//9N8c2iv51744oYGdnp1GiBwcHw9HREdOnT8/xM548eQJAXKmekpKiMq9kyZJwcHBQDj/14sWLbGcRFH3p3x2i6vbt2yp9hJnhmTBhAmxtbdGvX78c75qnydminPKMiLBgwYIc23fo0AHm5ub466+/sHHjRrRq1UplfNx79+7h+vXran/+x7Ro0QIAso20MHfuXADiDxOFD+VW586dER4ejn379mWbl5CQgIyMDI3jUqyvOrlsbm6e7f9j4cKFOu0f3b59e5ibm2Pq1KnZPpuI8OzZM7WW86HvVHGi4N3+vYoh3DQVHh6u0lf3/v372LZtG5o1a6bcNl+/fo3r16/j6dOnGi+f6Q7nZ3aHDx/Ocf+rjXpA3bxOSkrK9r1VrFgRZmZmymP98+fPsy0/p3pAV/gMuZFp1qwZrKys0Lp1awwePBgvX77EihUr4O7urvJXpiYmTJiAP//8E02bNsXIkSOVwx56e3vj+fPnKn+xhoSEYNWqVYiOjlYOLZiTmTNnIjIyEu3bt0elSpUAAGfPnsXq1atRoEABjBkzRtm2evXqWLJkCX788UeUKlUK7u7uHx0q0dHREUuWLEGvXr1QrVo1dO3aFW5ubrh37x527dqFevXq4ddff8WNGzfQpEkTdO7cGb6+vrCwsMCWLVsQFxeHrl27AgBWrVqFxYsXo127dihZsiSSk5OxYsUKODo6KnesANCkSRMAH+72wKRXunRprFu3Dt26dUPZsmWVd+okIkRHR2PdunUwMzPLsb/4+8qVK4eSJUviq6++wsOHD+Ho6Ih//vnng30y3d3d0ahRI8ydOxfJycno0qWLyvwvvvgCR48e1UpXscqVK6N3795Yvnw5EhISEBAQgFOnTmHVqlVo27YtGjVqpGz7odwaP348tm/fjlatWimHQX316hUuXbqETZs2ISYmBq6urhrFpbgIcdSoUQgODoa5ubkyz97XqlUrrFmzBk5OTvD19UV4eDgOHDigHDZOF0qWLIkff/wRISEhiImJQdu2beHg4IDo6Ghs2bIFgwYNwldfffXJ5XzoO23WrBm8vb3Rv39/jB8/Hubm5vjjjz+U+yZNVKhQAcHBwSrDHgJQucPxqVOn0KhRI0yePFnSW6IzVZyf2Y0cORKvX79Gu3btUK5cOaSlpeHkyZP4+++/UaxYMZWLp6tXr44DBw5g7ty5KFy4MIoXL/7RoWzVzetDhw5hxIgR6NSpE8qUKYOMjAysWbMG5ubm6NChAwBxB9Zjx46hZcuW8PHxQXx8PBYvXoyiRYuqjHmuM3obz4XlSk5DFW3fvp0qVapE1tbWVKxYMZo5cyb98ccfOQ6t1bJly2zLDAgIyDbU0rlz56hBgwYkl8upaNGiNGPGDPrll18IAMXGxirbKYaJ+9RQiCdOnKDhw4dThQoVyMnJiSwtLcnb25v69OlDt2/fVmkbGxtLLVu2JAcHBwKgjE0xxNSHhiQ8fPgwBQcHk5OTE1lbW1PJkiWpT58+yuHCnj59SsOHD6dy5cqRnZ0dOTk5Ue3atWnDhg3KZZw9e5a6detG3t7eJJfLyd3dnVq1aqUy5Jjiu/Tx8fnoOjPDcOvWLRo6dCiVKlWKrK2tycbGhsqVK0dDhgyh8+fPq7Tt3bs32dnZ5bicq1evUlBQENnb25OrqysNHDiQLly4QAAoNDQ0W/sVK1YQAHJwcKA3b96ozFMMo/cpn9rmFdLT02nq1KlUvHhxsrS0JC8vLwoJCaGUlBSVdh/KLSIxfGhISAiVKlWKrKysyNXVlerWrUs///yzcjgyxbBqOQ0diveGSsvIyKCRI0eSm5sbyWQylfV9v+2LFy+ob9++5OrqSvb29hQcHEzXr18nHx8flSHPNB32UJ3hS//55x+qX78+2dnZkZ2dHZUrV46GDx9OUVFRyjYBAQHk5+eX4/s/9p1GRkZS7dq1ycrKiry9vWnu3LkfHPYwp30zkfiuhg8fTn/++SeVLl2a5HI5Va1aNdt3oPhu3v1emW5xfuYuP/fs2UP9+vWjcuXKkb29PVlZWVGpUqVo5MiRFBcXp9L2+vXr1LBhQ+WwyYrPUwx7+KEhCT+V13fu3KF+/fpRyZIlydramgoUKECNGjWiAwcOKJdx8OBBatOmDRUuXJisrKyocOHC1K1bt2zDT+qKjEiPV/cxozJmzBgsW7YML1++1PnFDIwxxhhj+RX3IWcAoHI7cEDcpW7NmjWoX78+F+OMMcYYYzrEfcgZAHEjisDAQJQvXx5xcXH4/fffkZSUhO+++07q0BhjjDHGTBoX5AyAuDJ806ZNWL58OWQyGapVq4bff/8dDRs2lDo0xhhjjDGTxn3IGWOMMcYYkxD3IWeMMcYYY0xCXJAzxhhjjDEmIe5DroasrCw8evQIDg4OGt/WlbHcIiIkJyejcOHCMDPjv511jfOc6RvnuH5xjjN90yTHuSBXw6NHj+Dl5SV1GCyfun//vlp3l2R5w3nOpMI5rh+c40wq6uQ4F+RqcHBwACC+UEdHR4mjYflFUlISvLy8lNsf0y3Oc6ZvnOP6xTnO9E2THOeCXA2Kn7YcHR05iZne8U+r+sF5zqTCOa4fnONMKurkOHdaY4wxxhhjTEJckDPGGGOMMSYhLsgZY4wxxhiTkOQF+cOHD9GzZ08ULFgQNjY2qFixIs6cOaOcT0T4/vvvUahQIdjY2CAoKAg3b95UWcbz58/Ro0cPODo6wtnZGf3798fLly9V2ly8eBENGjSAtbU1vLy8MGvWLL2sH2P5Hec4Y6aNc5yxvJO0IH/x4gXq1asHS0tL7NmzB1evXsWcOXPg4uKibDNr1iz88ssvWLp0KSIiImBnZ4fg4GCkpKQo2/To0QNXrlxBWFgYdu7ciWPHjmHQoEHK+UlJSWjWrBl8fHwQGRmJ2bNnY8qUKVi+fLle15ex/IZznDHTxjnOmJaQhL7++muqX7/+B+dnZWWRp6cnzZ49WzktISGB5HI5/fXXX0REdPXqVQJAp0+fVrbZs2cPyWQyevjwIRERLV68mFxcXCg1NVXls8uWLatWnImJiQSAEhMTc5yfnKzWYhjTyKe2O2NgLDlOZBrfNzM8Hzs+mMI2xznO8ruMDKLXr3Oep8k2J+kZ8u3bt6NGjRro1KkT3N3dUbVqVaxYsUI5Pzo6GrGxsQgKClJOc3JyQu3atREeHg4ACA8Ph7OzM2rUqKFsExQUBDMzM0RERCjbNGzYEFZWVso2wcHBiIqKwosXL7LFlZqaiqSkJJXHh/z3H1CsGLBlS66/BsZMlqHmOKBZnjOWGydOAL6+wOrVUkeiO5zjLD+7cwcIDARGj877siQtyO/cuYMlS5agdOnS2LdvH4YOHYpRo0Zh1apVAIDY2FgAgIeHh8r7PDw8lPNiY2Ph7u6uMt/CwgIFChRQaZPTMt79jHfNmDEDTk5OysfH7uz166/As2dA+/bADz8ARJp8A4yZNkPNcUCzPGdME5mZwPTpQEAAcP8+MG8ekJUldVS6wTnO8iMiYNkyoFIl4Phx4K+/gEeP8rZMSQvyrKwsVKtWDdOnT0fVqlUxaNAgDBw4EEuXLpUyLISEhCAxMVH5uH///gfbrlwJjBolnn//PdClC/DqlX7iZMzQGWqOA5rlOWPqio0FgoOBb74RhXmPHsCxY4CZ5EMo6AbnOMtvHj4EWrQAhgwR9V7DhsDFi0DhwnlbrqS7iEKFCsHX11dlWvny5XHv3j0AgKenJwAgLi5OpU1cXJxynqenJ+Lj41XmZ2Rk4Pnz5yptclrGu5/xLrlcrryT16fu6GVhASxYAKxYAVhaAhs3AvXrA/+/Cozla4aa44Bmec6YOvbvBypXBg4eBGxtgT/+ANasAdS4a7bR4hxn+QWROBNesSKwdy8glwNz5gCHDwPFi+d9+ZIW5PXq1UNUVJTKtBs3bsDHxwcAULx4cXh6euLgwYPK+UlJSYiIiIC/vz8AwN/fHwkJCYiMjFS2OXToELKyslC7dm1lm2PHjiH9/9q787CoyjYM4PeAgGwDKgruS26puJtSapgomVmZnxqZW5mZVm5p0pdrpaZlmrlki1pWflnumkUimkq4JK6IG4opS6kwgogs7/fHG6OjqDMwM+8s9++65vLMzOHMM8fzzHnmzLvk5enXiYqKQoMGDQx6gpfGkCFAdDRQsSIQHw+0bi1/xiByZo6U40R3k5cHREbKK+Pp6fKEvW8fMHgwYMSM2XaNOU7O4J9/ZAuI558HrlwBWrUC/vwTGDPGjL9+mbmzqUn27NkjypQpI95//31x8uRJ8e233wovLy+xYsUK/TozZ84U/v7+Yt26deLQoUPi6aefFrVr1xY5OTn6dR5//HHRokULERcXJ3bu3Cnq1asnIiIi9M9nZGSIwMBA0b9/f3HkyBGxcuVK4eXlJT777DOj4jSll+y5c0I0by4EIISbmxBffGHCDiG6hSOMCGAvOS6EY+xvsr6zZ4UICZGf+YAQr7569xEXbucIxxxznBzdhg1CBAbK/HZ1FWLKFCFu3DDub0055pQW5EIIsWHDBtGkSRPh4eEhGjZsKJYsWWLwfGFhoZg4caIIDAwUHh4eonPnziIxMdFgnUuXLomIiAjh4+MjtFqtGDx4sLh621hTBw8eFO3btxceHh6iatWqYubMmUbHaGoSZ2UJ8Z//3PyAfuMNIfLyjH45IiGE45w87CHHhXCc/U3W89NPQvj7y895Pz8hVq0y7e8d5ZhjjpMjyswU4qWXbtZyDz4oxC0jcxq5DeOPOY0QHBfkfnQ6Hfz8/JCZmWl0GzQhgPfekx09ASAsDPjf/4Dy5S0YKDmUkhx3VHLc32Ss69eBN98EFiyQ99u2lW1LTW1HymPOuri/yVgxMcCgQcC5c7LZ2ejRsqbz9DRtO6Yccw7a71s9jQaYOBFYvRrw9gZ++01+aCckqI6MiIhKKjERaNfuZjE+fjzw++/m6dRFRGrl5Mjiu1MnWYzXqiWL848+Mr0YNxULcgvr2RPYvRuoWRM4dUoW5Zs2qY6KiIhM9fXXsjPXwYOyA//PPwMffCBH2CIi+7Z3L9CyJTB3rrz/8styOMOOHa3z+izIraBpU/kf3bEjcPUq0KMHMGsWJxEiIrIHWVnAgAHAwIFy3OFOneRoWo8/rjoyIiqtvDzZvDgkBDh+HKhcWV44XbLEukOWsiC3kooVgagoYOhQWYi/9RbQv7/8eYSIiGxTfLy8Kv7NN3J4s2nT5Gd5aScBISL1jh6VTdDefVdO5PXcc8Dhw3LiH2tjQW5F7u7A4sWy7aGrK/Dtt3Jq5QsXVEdGRES3EkJ+VrdrB5w4AVSrJtuSTpwoP7+JyH4VFAAffnhzPPHy5eXAG99/D1SooCYmFuRWptEAw4fLGd3Kl5dNWdq0AfbsUR0ZEREBcuKPXr2A114DcnNlM8P4eKBDB9WREVFpnTkjm52NGyfz+4kngCNHgD591MbFglyRxx6TxXjjxkBKimxfvmKF6qiIiJzb7t1A8+bAmjWys+bcucC6dequmhGReQgh24U3bSpHRvLxAT7/HNi4UbYbV40FuUJ16gCxscBTT8lvaf37yyG0CgpUR0ZE5FwKC4GZM+XFkeRk4IEH5OfzyJHyl00isl8XLwLduwOvvCI7ZnfsKEdQGTLEdvKbBblivr7ySsx//yvvz54tC/TMTLVxERE5i7Q0OWJKZKS8IBIRIduVtmqlOjIiKq2VK4EmTeQwpR4ewJw5wLZttjd3AAtyG+DiImeA+v57oGxZYPNm2ZHo5EnVkRERObbffgOaNZMjp3h6Al9+KTvccyJHIvv2zz9A377yC/aVK0Dr1sCBA3LiHxcbrH5tMCTn9dxzwM6dQNWqcizMhx6SJwkiIjKv/Hz5y2TXrvIKeZMmwL59wIsv2s5P2ERUMhs3ypz+4QegTBlgyhTZP+TBB1VHdncsyG1Mq1bypNCuHZCRAXTrBnzyCScRIiIyl+RkOeTs9Onys/WVV+RIV40aqY6MiEpDp5Ptwnv0kF+0GzUC/vgDmDzZ9mfUZUFug4KC5Hi3AwfK9owjR8opXHNzVUdGRGTf1q6Vo6js3i2bpfzwg5wfwtNTdWREVBoxMXIElS+/lL9yjR0L7N9vP31BWJDbKA8PYOlS4KOPZFunL78EOneW3/iIiMg0ubnAG28APXvK9qRt2sj2pL17q46MiEojJ0e2C+/UCTh3TnbWjImRE/+ULas6OuOxILdhGg0wZgywaRPg5wfs2nXzJEJERMY5cQIICQHmz5f333xT9tepU0dtXERUOvv2AS1byvkCAGDoUODgQTmsob1hQW4HHn8ciIsD6tcHzp8HHnkEWLVKdVRERLZvxQp5wj5wAAgIkBc4Zs8G3N1VR0ZEJZWXJ9uFt2snB8GoXFnm9mefyeGk7RELcjvRoIEsysPD5c8zffrIg7GwUHVkRES2JysLGDRITriWnQ2EhsorZ088oToyIiqNo0dlIT5tmuxn99xzwJEj9p/bLMjtiL+//AY4dqy8P20a8J//yBMPERFJhw7JMYeXL5d9cKZOleONV6miOjIiKqmCAtmvrlUrOXFX+fLA//4n53ApX151dKXHgtzOuLrKjgrLlsmfXNesAR5+GLhwQXVkRETqffmlnMMhMVEW4NHRwKRJ8rOTiOzThQvAY4/J/h+5uUD37vKqeJ8+qiMzHxbkdmrgQGD7djlE4uHD8ufYv/5SHRURkRqFhUBkpByDuOiEffCgHG+ciOzXrl3yF68dOwAfH+CLL4ANG2S7cUfCgtyOtWsnB7yvXRs4dUqeeJKTVUdFRGRd168Dzz8PzJwp70+ZIk/YAQFKwyKiUvrsMzmcYWoqEBwMxMcDL73kmLPpsiC3czVryvE269QBzpyRV8rPnVMdFRGRdVy6BISFybakbm6y3fjkyY55wiZyFrm5cgbdYcPkiCp9+gCxscADD6iOzHJYkDuAGjVkUf7AA0BSkizKz55VHBQRkYWdPi3HF9+1S87VsGULMGCA6qiIqDRSUuRV8SVL5BfrmTOBlSsBb2/VkVkWC3IHUb26LMrr1pXFeGioLM6JiBxRbKxstnfypLwosWuX7PRFRPbrjz/kKCqxsXJkuZ9/Bt56yzl+8WJB7kCqVZNFeb16stlKaKhsxkJE5Eh++kkW3//8I0/ecXFA48aqoyKi0vjiC9kXLiVF5vPevXLuFWfBgtzBVK0qi/L69WUHz9BQ+bMuEZG9E0KOQ9y7t+zI2aPHzdGmiMg+3bgBDB8OvPyyXH72WXmlvG5d1ZFZFwtyB1SliizKGzYEzp+XRfmpU6qjIiIqufx84LXX5DjEQsjlNWscv10pkSNLTQU6dwYWLZLNUt57D/jxRzm8obNhQe6gKlcGtm0DHnxQjk8eGirbWhIR2ZusLKBnT2DhQnnSnjMH+OQTTvZDZM/27JHji+/cKTtlb9gA/Pe/ztFevDhKC/IpU6ZAo9EY3Bo2bKh//vr16xgxYgQqVKgAHx8f9OrVC2lpaQbbSE5ORvfu3eHl5YVKlSph3LhxyM/PN1gnJiYGLVu2hIeHB+rWrYtly5ZZ4+0pFxQki/JGjeQsV6GhcvY6ImthjlNppaTIdqUbNwJly8qrZ6NHO+9J29Ywx6kkli4FOnaUtcmDD8rivHt31VGppfwKeePGjZGSkqK/7dy5U//c6NGjsWHDBqxatQrbt2/HxYsX8eyzz+qfLygoQPfu3XHjxg3s3r0by5cvx7JlyzBp0iT9OklJSejevTs6deqE+Ph4jBo1CkOGDMEvv/xi1fepSmCgLMobNwYuXpRDCR0/rjoqcibMcSqpo0flSCp//glUrCg/y245PMhGMMfJWHl5wOuvAy++KMcaf+YZ2V68fn3VkdkAodDkyZNFs2bNin0uIyNDuLm5iVWrVukfS0hIEABEbGysEEKIzZs3CxcXF5GamqpfZ9GiRUKr1Yrc3FwhhBDjx48XjRs3Nth23759RXh4uNFxZmZmCgAiMzPT6L+xNenpQgQHCwEIERQkxLFjqiOi+3GE485eclwIx9jfjuS334TQauVnVv36Qpw6pToi83OEY445TsZKSxOiY0eZ04AQU6cKUVCgOirLMuWYU36F/OTJk6hSpQrq1KmDfv36Ifnfud/379+PvLw8hIWF6ddt2LAhatSogdjYWABAbGwsgoODERgYqF8nPDwcOp0OR48e1a9z6zaK1inaRnFyc3Oh0+kMbvauYkUgOhpo2lR2oujUCTh2THVU5AxsMccBx8xzR7F8OfD444BOB3To4Pgz9Nk75jjdz/79sr34jh2Ary+wbh0waRLgorwKtR1Kd0Xbtm2xbNkybNmyBYsWLUJSUhI6dOiAq1evIjU1Fe7u7vD39zf4m8DAQKSmpgIAUlNTDZK46Pmi5+61jk6nQ05OTrFxzZgxA35+fvpb9erVzfF2lQsIkEV58+ZAWposyo8cUR0VOTJbzXHAcfPcngkBTJkCDBokR1WJiACiooDy5VVHRnfDHKf7+eYboH17OepbgwayvfhTT6mOyvaUUfni3bp10y83bdoUbdu2Rc2aNfHDDz/A09NTWVyRkZEYM2aM/r5Op3OYRK5QAfjtN6BLF+DAATm5xtatQHCw6sjIEdlqjgOOnef26MYNYMgQefIGgLffBt59l1fQbB1znO4mPx8YNw6YO1fef/JJYMUKOaIK3cmmPur8/f1Rv359nDp1CkFBQbhx4wYyMjIM1klLS0PQv7NABAUF3dFbu+j+/dbRarV3/bDw8PCAVqs1uDmSoqK8ZUvg779lUX7okOqoyBnYSo4Djp/n9uTKFdlE5Ztv5FCGn38OvP8+i3F7xBwnQM6iGx5+sxifNEk2U2Exfnc29XGXlZWF06dPo3LlymjVqhXc3NywdetW/fOJiYlITk5GSEgIACAkJASHDx9Genq6fp2oqChotVo0atRIv86t2yhap2gbzqp8eVmUt2olE+exx4CDB1VHRY6OOU63O3sWeOQROYKKry+waZO8Uk72iTlOBw7I9uLR0XKCn9WrgalT+QX7vqzQyfSuxo4dK2JiYkRSUpLYtWuXCAsLEwEBASI9PV0IIcSwYcNEjRo1RHR0tNi3b58ICQkRISEh+r/Pz88XTZo0EV27dhXx8fFiy5YtomLFiiIyMlK/zpkzZ4SXl5cYN26cSEhIEAsWLBCurq5iy5YtRsfpyD2zr1wRok0b2eO5fHkhDhxQHREVcYTjzl5yXAjH2N/2Zu9eIQID5edP1apCxMerjsi6HOGYY47Trb77TghPT5nTdesKcfSo6ojUMuWYU1qQ9+3bV1SuXFm4u7uLqlWrir59+4pTt4xtlZOTI4YPHy7KlSsnvLy8RM+ePUVKSorBNs6ePSu6desmPD09RUBAgBg7dqzIy8szWGfbtm2iefPmwt3dXdSpU0csXbrUpDgdPYkzMoRo21YmULlyQuzfrzoiEsIxjjt7yXEhHGN/25P164Xw8pKfO02bCnH+vOqIrM8RjjnmOAkhRF6eEGPH3hzSsFs3ecHP2ZlyzGmEEELNtXn7odPp4Ofnh8zMTIdtg5aZKdtw/vEH4O9/szkLqeMMx50t4f62nk8/BUaOBAoLZTvTH34AnHGX85izLu5vy7h0CXjuOVk3ALJD9rRpsj+IszPlmGOLHgIgO1r88gsQEgJkZABhYcDevaqjIiJHUlgIjBkjZ+orLARefhnYsME5i3EiR3DwINCmjSzGvb2BVatkh2wW46ZjQU56Wq0syh95RBblXbrI8UKJiErr2jWgd2/g44/l/Zkzgc8+A9zc1MZFRCXzww/Aww8DSUlAnTryF/b//Ed1VPaLBTkZ8PUFfv5ZDuKfmSmL8j/+UB0VEdmz9HQ5ktPq1YC7O/D998BbbwEajerIiMhUBQXAhAlA377yi3bXrvIX9SZNVEdm31iQ0x2KivKOHeXU1V27yqmriYhMlZgItGsHxMXdHG71uedUR0VEJXHlCtC9O/DBB/L+W28BmzdzNl1zYEFOxfLxkUkWGgpcvSqL8l27VEdFRPZkxw7ZL6XoJ+3du4EOHVRHRUQlceSIbC/+yy+ApyewcqVsesb24ubBgpzuytsb2LgR6NQJyMqSo7Ds3q06KiKyB6tWySZvV67IK+R//AE0aKA6KiIqiXXrZB6fPg3UqiV/Ne/bV3VUjoUFOd1TUVH+2GOyKO/RAzh1SnVURGTL1q8HIiKAGzeAXr3kjH0VK6qOiohK4ocfZB5nZwOdOwP79gHNmqmOyvGwIKf78vKSQ5O1aQNcvizbj12+rDoqIrJF27YBffrIjl8DB8qTuaen6qiIqCRWrgSef17m84ABwJYtQIUKqqNyTCzIySheXvKqV/XqwIkT8tvyjRuqoyIiW7JnD/DUU0BuLtCzJ/DFF4ALzzJEdum774B+/WQxPngw8NVXQJkyqqNyXPyoJKMFBcnmKz4+QEwM8OqrcpJcIqIjR4Bu3WTTtrAwObQhT95E9mnFCqB/fzmB10svyS/X7LxpWSzIySRNmwL/+5+86vXVV8CsWaojIiLVTp+WIzFdviw7fq1ZA3h4qI6KiEpi+XLZPKWwEBg6FFiyhL90WQN3MZnsiSeAefPk8oQJwI8/qo2HiNS5cEGOppKSAgQHy+FSfXxUR0VEJbF0qWyeIgQwbBiwaBGLcWvhbqYSee014PXX5XL//rLtKBE5l3/+kVfGk5KAunWBX38FypVTHRURlcSXX8rmKUIAw4cDCxeyGLcm7moqsTlz5NXy69dlR67kZNUREZG16HSyzfixY0C1anIGzqAg1VERUUksWQIMGSKL8ddfBz79FNBoVEflXFiQU4mVKSOHRAoOBtLSgCeflCdpInJsOTnyS/i+fUBAABAVBdSsqToqIiqJxYuBV16RyyNHyiapLMatjwU5lYqvrxx5JSgIOHwYeO45ID9fdVREZCl5eUDv3sD27YBWK6fRbthQdVREVBILFsgR0wBgzBjg449ZjKvCgpxKrUYNOUa5pyfw88/A6NGqIyIiSyia7GfTJqBsWfllvGVL1VERUUnMny/7gwHAuHHAhx+yGFeJBTmZRZs2ctxSQLY9mz9fbTxEZF5CyJN30fjiq1cDHTqojoqISmLuXOCNN+TyW28BH3zAYlw1FuRkNs8+C8ycKZdHjZJX0YjIMbz9tmxrqtEA334rO3QSkf2ZM+fmL9lvvw3MmMFi3BYYNY/aJ598YvKGBw8eDF9fX5P/juzb+PHAiRNy0qDnngN27ZKTCZFtY47TvcycefPL9pIlQJ8+auOhkmGe04cfyuYpADBxIjB1KotxW6ER4v6Tn7u4uKBatWpwNXLe1PPnz+PEiROoU6dOqQO0BTqdDn5+fsjMzIRWq1Udjs27cQN4/HFg2zagenUgLg6oXFl1VPbHmseds+c4wDy/m8WLb3b6+vBDYOxYtfE4Emsfc86e586e4x98ICfzA4DJk4EpU5SG4xRMOeaMukIOAPv27UOlSpWMWpffpp2buzvw009ASAiQmCiHR9u+HfDyUh0Z3QtznG733XdyghAAeOcdFuOOgHnunKZPB/77X7k8dSowaZLaeOhORrUhnzx5MnxMmAv57bffRvny5UscFNm/cuXkCAwVKsixivv3BwoLVUdFd8Mcp9tt2AAMGHCzM+e0aaojotJinjun9967WYy/9x6LcVtlVJMVZ+fsP3OVxu+/A2FhshnLW2/dbIdK98fjzrq4v2+KiZHNznJzgRdeAJYv5xTalsBjzrqccX9PnXqzacr06UBkpNJwnI4pxxw/YsmiOnQAvvxSLn/wwc1lIrJNe/cCPXrIYvzpp4GlS1mME9kbIQzbiX/wAYtxW2fyx+ylS5cwYsQINGrUCAEBAShfvrzBjeh2L7xw8yeyYcOA6Gi18dC9Mced19Gj8sp4Vhbw2GPAypVyzHFyPMxzxyWEPOcWNTObPVuOgEa2zeSP2v79++PUqVN46aWXEBgYCA3HyyEjTJkCnDwpJxXp1QuIjeV027aKOe6czpwBunQBLl8G2rYF1q6Vs3GSY2KeOyYhZHvxGTPk/VvHHCfbZnJB/vvvv2Pnzp1o1qyZJeIhB6XRyLHJz56VxXj37nI4xIAA1ZHR7ZjjzufiRVmMp6QATZoAmzcDHGDDsTHPHY8QcljDWbPk/blzgZEjlYZEJjC5yUrDhg2Rk5Nj9kBmzpwJjUaDUaNG6R+7fv06RowYgQoVKsDHxwe9evVCWlqawd8lJyeje/fu8PLyQqVKlTBu3Djk5+cbrBMTE4OWLVvCw8MDdevWxbJly8weP91f2bLAunVA7dryalzPnrKdKtkWS+U4wDy3RZcuAV27ypx84AHg118BtlhwfDyXOxYhZLOUomJ8/nwW43ZHmGjPnj3iscceEzExMeKff/4RmZmZBreS2LNnj6hVq5Zo2rSpGDlypP7xYcOGierVq4utW7eKffv2iXbt2omHH35Y/3x+fr5o0qSJCAsLEwcOHBCbN28WAQEBIjIyUr/OmTNnhJeXlxgzZow4duyYmD9/vnB1dRVbtmwxOr7MzEwBoMTvjwwdPSqEVisEIMQLLwhRWKg6Ituk6rizRI4XbZd5blt0OiHatJG5WKWKEGfOqI7Iuag85ngudxyFhUKMHi3zGBBiwQLVEVERU445kwvyEydOiNatWwsXFxeDm0ajES4uLiYHe/XqVVGvXj0RFRUlHn30UX0SZ2RkCDc3N7Fq1Sr9ugkJCQKAiI2NFUIIsXnzZuHi4iJSU1P16yxatEhotVqRm5srhBBi/PjxonHjxgav2bdvXxEeHm50jI6axCr9+qsQrq7yw2PaNNXR2CZVx525c1wI5rktunZNiNBQmYMVKsgvymRdKo85nssdQ2GhECNH3izGFy9WHRHdypRjzuQmK/369YObmxu+++47bN26FdHR0YiOjsa2bdsQXYLhM0aMGIHu3bsjLCzM4PH9+/cjLy/P4PGGDRuiRo0aiI2NBQDExsYiODgYgYGB+nXCw8Oh0+lw9OhR/Tq3bzs8PFy/jeLk5uZCp9MZ3Mi8unQBFi6Uy5Mmyc6eZBvMneMA89zW5OcDffvK8cZ9fYEtW4BGjVRHRdbEc7n9E0I2S5k3T95fsgR45RW1MVHJmdyp88iRIzhw4AAaNGhQ6hdfuXIl/vzzT+zdu/eO51JTU+Hu7g5/f3+DxwMDA5Gamqpf59YELnq+6Ll7raPT6ZCTkwNPT887XnvGjBmYOnVqid8XGWfoUODECeCjj4AhQ4AWLTjyii0wZ44DzHNb9O67cibOsmXlv61bq46IrI3ncvu3eLFsK67RAJ9/Drz0kuqIqDRMvkLeunVrnD9/vtQvfP78eYwcORLffvstytrY2FqRkZHIzMzU38zxfql4H3wgxzu+dg2IiGAnT1tgrhwHmOe2aPt2OX02ICfqevRRtfGQGjyX27e9e4GifrOzZrEYdwQmXyF//fXXMXLkSIwbNw7BwcFwc3MzeL5p06ZGbWf//v1IT09Hy5Yt9Y8VFBRgx44d+PTTT/HLL7/gxo0byMjIMPhmnZaWhqCgIABAUFAQ9uzZY7Ddop7bt65ze2/utLQ0aLXaYr9RA4CHhwc8PDyMeh9UOq6uwDffAE2bAvHxciaxOXNUR+XczJXjAPPc1ly6BPTrBxQWAoMHA88/rzoiUoXncvt16RLwn/8AN27I0crGjlUdEZmFqQ3UNRrNHbeSdATR6XTi8OHDBrfWrVuLF154QRw+fFjfEeTHH3/U/83x48eL7QiSlpamX+ezzz4TWq1WXL9+XQghO4I0adLE4LUjIiKcviOIrVm//manlJ9/Vh2NbVB13Jkrx4VgntuSwkIhnnpK5liDBkJkZamOiFQeczyX26eCAiGeeELm8QMPCJGRoToiuheLjrJy9uzZe95K49ae2ULIoZJq1KghoqOjxb59+0RISIgICQnRP180VFLXrl1FfHy82LJli6hYsWKxQyWNGzdOJCQkiAULFnCoJBv12mvyQ6ZSJSFSUlRHo56q486SOS4E81yV+fNlfrm7C3HggOpoSAi1xxzP5fbpvfdkHpctK0R8vOpo6H4sWpDfa6MnT540dXMGbk/inJwcMXz4cFGuXDnh5eUlevbsKVJuq9TOnj0runXrJjw9PUVAQIAYO3asyMvLM1hn27Ztonnz5sLd3V3UqVNHLF261KS4HCGJ7UFOjhDBwfLDpmtXeSXAmak67iyZ40Iwz1WIjxfCw0Pm1rx5qqOhIiqPOZ7L7c9vvwnh4iLz+KuvVEdDxjDlmNMIIYQpTVw6dOiAqKioOzpvJCYmonPnzvjrr79K2nrGZul0Ovj5+SEzMxNarVZ1OA7t2DGgVSvg+nXgww+du22cquPOGXMccNw8z86WOZWYCPToIWfL1WhUR0WA2mPOGfPcnnP8wgWgZUsgPR148UXZIZtsnynHnMmjrPj4+ODZZ581mNI2ISEBoaGh6NWrl+nREt2iUSPg44/lcmQksH+/2nicEXPcsYwcKYvxKlWAr75iMU4S89x+5OXJeQPS0+UACJ9+qjoisgSTC/LVq1cjMzMT/fr1gxACR44cQWhoKCIiIjCvaHR6olJ45RXZczwvTw6FmJWlOiLnwhx3HP/7n7ySptEA334LBASojohsBfPcfkRGArt2AVot8NNPwF0GlSE7Z3JB7unpiU2bNiExMRF9+vRB586dMWDAAMzhWHVkJhoN8MUXQLVqwMmTwBtvqI7IuTDHHUNSkpx8CwDeeQcIDVUaDtkY5rl9WL1aTp4HAMuWAXXrKg2HLMioNuTFTTebkpKCLl264Mknn8TMmTP1j9tbuyxj2HO7M3u2fTvQqZMcDPH774HnnlMdkXVZ87hz9hwHHCvP8/KADh2AuDjgkUeAmBigjMmzTpClWfuYc/Y8t7ccP3lSzqKr08n+VB9+qDoiMpUpx5xRBbmLiws0xTQ8LPpTjUYDIQQ0Gg0KCgpKGLbtsrckdiQTJ8pZBbVaOXFQ7dqqI7Ieax53zp7jgGPl+YQJchZcf3+ZNzVrqo6IimPtY87Z89yecjwnB2jXDjh0CGjfHoiOBm6bu4nsgCnHnFHXTLZt22aWwIhMNXkysHUrEBsrZxjcsYNX+iyBOe44oqJkMQ7I9uMsxqkI89x+vPaaLMYrVZJ9QViMOz6jSptHH33U0nEQFatMGeC774DmzWVRPnUq8O67qqNyPMxxx5CeDgwYIJeHDQOefVZtPGRbmOf24auv5M3FRZ7/qlRRHRFZg1GdOg8dOoTCwkKjN3r06FGDoZSISqNWLeCzz+Ty++/LtuVkXsxx+1dYCAwcCKSmAo0bA+ybR7djntu++HhgxAi5PG0a0Lmz0nDIiowqyFu0aIFLly4ZvdGQkBAkJyeXOCii2/XtCwweLDt4vvACYMLhSEZgjtu/jz8GtmwBypaVP3FzaDS6HfPctmVmAv/5j5wY74kn5HCH5DyMarIihMDEiRPh5eVl1EZv3LhRqqCIivPJJ3Is1hMngCFD5HBQnOTEPJjj9m3fvpsn77lz5RVyotsxz22XEPKi0+nTst/HN9/IJivkPIwqyDt27IjExESjNxoSEgJPXp4hM/PxkcMftmsHrF0rm7EMG6Y6KsfAHLdfV6/KIUHz8oBevW6OPU50O+a57ZozB1izBnB3B1atAsqXVx0RWZtRBXlMTIyFwyAyTsuWwMyZckzW0aPlWMu8Glh6zHH7NXy4vKpWowbw+ef81Yjujnlum3buBN56Sy7PnQu0aaM0HFKEP4iQ3Rk1CggPl+3sIiLkv0TO6OuvgRUrAFdXORpDuXKqIyIiU6SlAX36AAUFwPPP81dfZ8aCnOyOiwuwfLkcn/XwYWDcONUREVnfiRPy6jgATJkiZ+QkIvtRVISnpAAPPiibYfIXLufFgpzsUmCgLMoB4NNPgfXr1cZDZE25ubLdeHY2EBrK0RiI7NHkyXIGTm9v4KefZD8pcl4syMluPf44MGaMXH7xReDCBbXxEFlLZCRw4ABQocLNJitEZD82bZLzagDAF1/IK+Tk3FiQk12bPh1o0UKOSz5ggPwJkMiRbd4sxxwHgGXLgKpVlYZDRCY6exbo318ujxghf+0iMrkgX758OTZt2qS/P378ePj7++Phhx/GuXPnzBoc0f14eMihEL285E9/s2erjsj+Mcdt18WLcjZOAHjjDeDJJ9XGQ/aLea5Gbi7Quzdw5Qrw0EPARx+pjohshckF+fTp0/XjksbGxmLBggWYNWsWAgICMHr0aLMHSHQ/DRoA8+fL5YkTgbg4tfHYO+a4bSoslFfV/vkHaN4cmDVLdURkz5jnaowZIyfyKl8e+OEHeVGJCDByHPJbnT9/HnXr1gUArF27Fr169cLQoUPxyCOPIDQ01NzxERll8GDgl1/kB9zgwUB8vJxggUzHHLdNy5ff7AC2ciVP5FQ6zHPr27ULWLhQLq9YIWfkJCpi8hVyHx8fXLp0CQDw66+/okuXLgCAsmXLIicnx7zRERlJowEWLwYqVgQSEuTkClQyzHHbk5MDTJokl6dMkb8KEZUG89y6hLg5GtKQIUC3bmrjIdtj8hXyLl26YMiQIWjRogVOnDiBJ554AgBw9OhR1KpVy9zxERmtXDnZhnzQIGDqVDlpUPXqqqOyP8xx2zN/PvDXX3I2ztdeUx0NOQLmuXVt2QL8/rv8ZWvyZNXRkC0y+Qr5ggULEBISgr///hs//fQTKlSoAADYv38/IiIizB4gkSkGDADatweuXZMzepLpmOO25fJlYMYMufzuu0DZsmrjIcfAPLeewsKbV8dffx2oVk1tPGSbNEIIoToIW6fT6eDn54fMzExotVrV4dB9HD4sh0IsKAB+/lmOV26PeNxZl63u73HjgA8/BJo2Bf78k2OOOxJbPeYclar9/f33ckZOrRY4c0bOH0DOwZRjzuQmKzt27Ljn8x07djR1k0RmFRwMjBwJzJkjf94/coRXFU3BHLcdyck3RxCaOZPFOJkP89w68vLk6F+A/HLNYpzuxuSCvLje1xqNRr9cwJlZyAZMmSJHojh9Wg4PV9Qhju6POW47Jk2S4xZ36mS/v/SQbWKeW8eXX8rzUKVKbEZJ92ZyG/IrV64Y3NLT07Flyxa0adMGv/76qyViJDKZr6+8Qg7I2TxPn1Ybjz1hjtuGQ4eAr7+Wyx98IEcSIjIX5rnlXbsGTJsmlydOBHx81MZDts3kK+R+fn53PNalSxe4u7tjzJgx2L9/v1kCIyqtPn2AL74AfvtNzmq4cSOLGmMwx23DhAlyqLQ+fYA2bVRHQ46GeW558+cDKSlArVrA0KGqoyFbZ/IV8rsJDAxEYmKiSX+zaNEiNG3aFFqtFlqtFiEhIfj555/1z1+/fh0jRoxAhQoV4OPjg169eiEtLc1gG8nJyejevTu8vLxQqVIljBs3Dvn5+QbrxMTEoGXLlvDw8EDdunWxbNmyEr9Psh8aDfDpp4CbG7B5M7BuneqI7Btz3Hq2bZMdksuUAd5/X3U05ExMzXPmePGuXJH9PgB5lZwT1dF9CRMdPHjQ4BYfHy9+/vln8eijj4pHHnnEpG2tX79ebNq0SZw4cUIkJiaKt99+W7i5uYkjR44IIYQYNmyYqF69uti6davYt2+faNeunXj44Yf1f5+fny+aNGkiwsLCxIEDB8TmzZtFQECAiIyM1K9z5swZ4eXlJcaMGSOOHTsm5s+fL1xdXcWWLVuMjjMzM1MAEJmZmSa9P7INb78tBCBEjRpCZGWpjsZ4qo47Z8xxIWwnzwsLhWjTRh6zI0YoDYUsTOUxZ648Z44Xb8IEmcNNmgiRn2/xlyMbZcoxZ3JBrtFohIuLi9BoNAa3kJAQkZCQUKKAb1WuXDnxxRdfiIyMDOHm5iZWrVqlfy4hIUEAELGxsUIIITZv3ixcXFxEamqqfp1FixYJrVYrcnNzhRBCjB8/XjRu3NjgNfr27SvCw8ONjslWTtRUMtnZQtSsKT8cJ0xQHY3xVB13zpjjQthOnv/vf/JY9fERIi1NaShkYSqPOUvmubPn+MWLQnh6yjxet86iL0U2zpRjzuQmK0lJSThz5gySkpKQlJSEc+fO4dq1a9i9ezcaNmxY4iv1BQUFWLlyJbKzsxESEoL9+/cjLy8PYWFh+nUaNmyIGjVqIDY2FgAQGxuL4OBgBAYG6tcJDw+HTqfD0aNH9evcuo2idYq2QY7Pywv45BO5/NFHwPHjauOxdcxxdfLygLfflstvvilHZiCyBEvkOXNcevddICcHCAkBevRQHQ3ZC5M7ddasWdOsARw+fBghISG4fv06fHx8sGbNGjRq1Ajx8fFwd3eHv7+/wfqBgYFITU0FAKSmphokcdHzRc/dax2dToecnBx4enreEVNubi5yc3P193U6XanfJ6n11FPAk0/Kjp0jRsiOnuzgWTxnyHHANvN8yRI5IlBgIDB2rOpoyJGZM8+Z4zedPg18/rlcnjGD5xkynlEF+SeffIKhQ4eibNmy+KToUuNdvPHGGyYF0KBBA8THxyMzMxM//vgjBg4ciO3bt5u0DXObMWMGpk6dqjQGMr9PPpGFeHS0HKOcs0Pf5Gw5Dthenl+9ChSFM3kyh0gj87NUnjPHb5o0CcjPl/MGPPqoVV+a7JxRBfnHH3+Mfv36oWzZsvj444/vup5GozH5ZO3u7o66desCAFq1aoW9e/di3rx56Nu3L27cuIGMjAyDb9dpaWkICgoCAAQFBWHPnj0G2yvqvX3rOrf36E5LS4NWq73rt+rIyEiMGTNGf1+n06F69eomvS+yPbVry+YAkybJq4/du8upjMn5chywvTz/6CPg77+BevWAIUOUhUEOzFJ5zhyXDh4EvvtOLk+fbrGXIQdlVEGelJRU7LIlFBYWIjc3F61atYKbmxu2bt2KXr16AQASExORnJyMkJAQAEBISAjef/99pKeno9K/jS2joqKg1WrRqFEj/TqbN282eI2oqCj9Norj4eEBDw8PS7w9UmzcOOCbb4CTJ+VVyHuck5yKs+U4YFt5npoKfPihXJ4+XQ7VSWRu1spzZ83x//5X/tu3L9CihdVelhyFFTqZ3tWECRPE9u3bRVJSkjh06JCYMGGC0Gg04tdffxVCyOGSatSoIaKjo8W+fftESEiICAkJ0f990XBJXbt2FfHx8WLLli2iYsWKxQ6XNG7cOJGQkCAWLFjAYQ+d3C+/yN7vrq5CxMerjubuHOG4s5ccF0Lt/h4+XB6TDz0khz0k58Acd5wc37Hj5nnlxAmzb57slCnHnFFXyG/9yed+5hTNV26E9PR0DBgwACkpKfDz80PTpk3xyy+/oEuXLgDkz2suLi7o1asXcnNzER4ejoULF+r/3tXVFRs3bsSrr76KkJAQeHt7Y+DAgZhWNFctgNq1a2PTpk0YPXo05s2bh2rVquGLL75AeHi40XGSY+naFfjPf4AffwSGDwd+/x1wMdsUWfaJOa7OyZOyMycAzJrFTmBkOZbIc+a4nFE3MlIuDxkim50RmUojhBD3W6lTp04G9//880/k5+ejQYMGAIATJ07A1dUVrVq1QnR0tGUiVUin08HPzw+ZmZnQstGxQ/jrL6BhQyA7G/jqK2DwYNUR3cmax52z5zigLs9795ZfDrt3l6MAkfOw9jHn7Hluqf29aZMcxatsWeDUKaBqVbNtmuycKcecUVfIt23bpl+eM2cOfH19sXz5cpQrVw4AcOXKFQwePBgdOnQoRdhE1lOtGjBlimxTPn488PTTQPnyqqNShzmuRlycLMY1GjlEGpElMc/Nr7Dw5tXxN95gMU4lZ9QV8ltVrVoVv/76Kxo3bmzw+JEjR9C1a1dcvHjRrAHaAl4hd0x5ebLjzdGjwLBhwKJFqiMypOq4c8YcB6y/v4UAOnUCtm8HBg0Cli61+EuSjVF5bnHGPLfE/v72W+CFFwA/P+DMGee+sEN3MuWYM7nlrE6nw99//33H43///TeuXr1q6uaIlHFzA4qaMn72GbB3r9p4bAVz3Do2b5bFuIfHzfHHiayFeV56N27IYXQB+Usri3EqDZML8p49e2Lw4MFYvXo1/vrrL/z111/46aef8NJLL+HZZ5+1RIxEFtOxI9C/v7xa+eqrQEGB6ojUY45bXkEBMGGCXH7jDaBGDbXxkPNhnpfeF1/Iq+KBgcDIkaqjIbtn6hAu2dnZ4tVXXxUeHh7CxcVFuLi4CHd3d/Hqq6+KrKwsk4eEsQeOMDQV3V1qqhB+fnLIqoULVUdzk6rjzhlzXAjr7u+lS+Xx5u8vxOXLFn85slEqzy3OmOfm3N9ZWUIEBso8XrDADMGRQzLlmDO5DXmR7OxsnD59GgDwwAMPwNvb20xfEWwP25A7vk8/BV5/HfD3BxITgX/np1BK9XHnTDkOWG9/5+QA9evLkX5mzZIdi8k5qc5xwLny3Jz7e8YMOfNz7drA8eOAu7uZgiSHYvZRVorj7e2Npk2blvTPiWzKq6/K4Q8PHACmTZMFurNjjlvG8uWyGK9eXX4JJFKJeW66/Hzgo4/k8rvvshgn8yhRQb5v3z788MMPSE5Oxo0bNwyeW716tVkCI7ImV1f5AfvYY7Jd4H//C1SurDoqdZjjlvPVV/LfUaPkuMVEqjDPS2bHDuDSJSAgAOjbV3U05ChM7tS5cuVKPPzww0hISMCaNWuQl5eHo0ePIjo6Gn5+fpaIkcgqQkOBRx4BcnOB2bNVR6MOc9xyjhyRo/mUKSOHSiNShXlecmvWyH+fekrmMpE5mFyQT58+HR9//DE2bNgAd3d3zJs3D8ePH0efPn1Qg0MFkB3TaICJE+Xy4sVAerraeFRhjltO0VjjTz5pG/0UyHkxz0tGCGDtWrncs6fSUMjBmFyQnz59Gt27dwcAuLu7Izs7GxqNBqNHj8aSJUvMHiCRNXXtCrRpIzvezZmjOho1mOOWkZcHrFghlwcPVhsLEfO8ZPbvl31AvL2BsDDV0ZAjMbkgL1eunH7SgKpVq+LIkSMAgIyMDFy7ds280RFZ2a1XyRcskO0EnQ1z3DI2b5a/ulSqBHTrpjoacnbM85Ipaq7SrRv7gJB5mVyQd+zYEVFRUQCA3r17Y+TIkXj55ZcRERGBzp07mz1AImt78kmgeXMgKwuYN091NNbHHLeMouYq/fvLWWKJVGKelwybq5ClmDwO+eXLl3H9+nVUqVIFhYWFmDVrFnbv3o169erhnXfeQbly5SwVqzK2MFYsWddPPwH/+Q/g5wecOyf/tTZVx50z5jhg2f2dlgZUrSpn6DxyBGjc2KybJzul8tzijHle2v194gTQoIHsyPn333LeCqJ7seg45OXLl9cvu7i4YELR/M8AcnJyTN0ckU3q2VMWTUePAvPnA++8ozoi62GOm9+KFbIYf+ghFuNkG5jnpitqrvLYYyzGyfxMbrJSnNzcXMyZMwe1a9c2x+aIlHNxkWORA8DHHwP/NrV0WszxkhPiZnMVduYkW8Y8v7eigpzNVcgSjC7Ic3NzERkZidatW+Phhx/G2n8bUi1duhS1a9fGxx9/jNGjR1sqTiKr69NHTnF++TKwaJHqaCyPOW4Z+/bJX1rKlgWee051NOTsmOclc/EiEBcnl596Sm0s5KCEkcaPHy/8/PxEr169ROXKlUWZMmXEyy+/LIKDg8X3338v8vPzjd2U3cnMzBQARGZmpupQyMqWLRMCEKJSJSGys6372tY+7pw5x4Ww3P4eNkweQ88/b9bNkgNQcW5x5jwvzf5euFDmcbt2FgiMHJYpx5zRbchXrVqFr7/+Gk899RSOHDmCpk2bIj8/HwcPHoRGo7HMtwUixZ5/Hpg6FUhKApYskdOdOyrmuPnl5ADffy+X2VyFbAHzvGTYXIUszegmK3/99RdatWoFAGjSpAk8PDwwevRoJjA5NDc3IDJSLs+aBVy/rjYeS2KOm9/atUBmJlCjhuwIRqQa89x0GRnAtm1y+ZlnVEZCjszogrygoADu7u76+2XKlIGPj49FgiKyJQMHAtWrAykpwFdfqY7Gcpjj5lfUmXPQINlRmEg15rnpNm0C8vOBRo1kvyIiSzC6yYoQAoMGDYKHhwcA4Pr16xg2bBi8vb0N1lu9erV5IyRSzN0deOst4LXXgJkzgSFD5GOOhjluXsnJwG+/yeVBg5SGQqTHPDcdm6uQNRhdkA8cONDg/gsvvGD2YIhs1YsvAu+9B5w/D3z9tSzKHQ1z3LyWL5dDHoaGAhxFjmwF89w0OTnAli1ymc1VyJKMLsiXFv32SuSEPD2BceOAsWOBGTPkFc8yJk+rZduY4+ZTWAgsWyaX2ZmTbAnz3DS//QZkZ8tmi/82vSeyCLZqJDLSK68AAQHAmTPAd9+pjoZs2e+/y+PE1xfo1Ut1NERUUkXNVZ55BmC/V7IkFuRERvL2llfIAeD99+VU6ETFKer827evPG6IyP7k5wPr18tlNlchS2NBTmSCESOAcuWAEyeAH39UHQ3ZoqtXbx4bbK5CZL927QIuXQLKlwc6dlQdDTk6FuREJvD1vTk50HvvybbCRLf64Qfg2jWgQQMgJER1NERUUkXNVXr0cLw+Q2R7lBbkM2bMQJs2beDr64tKlSrhmWeeQWJiosE6169fx4gRI1ChQgX4+PigV69eSEtLM1gnOTkZ3bt3h5eXFypVqoRx48YhPz/fYJ2YmBi0bNkSHh4eqFu3LpYV9bgiMtEbbwBaLXDkCLBunepobJsz5nhRc5XBg9nmlByfo+a4EIbtx4ksTWlBvn37dowYMQJ//PEHoqKikJeXh65duyI7O1u/zujRo7FhwwasWrUK27dvx8WLF/Hss8/qny8oKED37t1x48YN7N69G8uXL8eyZcswadIk/TpJSUno3r07OnXqhPj4eIwaNQpDhgzBL7/8YtX3S47B3x94/XW5/O678oObiudsOZ6YCOzeLScB6t/fqi9NpISj5nh8vJxLwNMT6NrVIi9BZEjYkPT0dAFAbN++XQghREZGhnBzcxOrVq3Sr5OQkCAAiNjYWCGEEJs3bxYuLi4iNTVVv86iRYuEVqsVubm5Qgghxo8fLxo3bmzwWn379hXh4eFGxZWZmSkAiMzMzFK9P3Icf/8thLe3EIAQGzda5jUc8biz1RwXwjz7e8IEeUw88USJN0FOhDluuzk+caLM5Z49jd480R1MOeZsqg15ZmYmAKB8+fIAgP379yMvLw9hYWH6dRo2bIgaNWogNjYWABAbG4vg4GAEBgbq1wkPD4dOp8PRo0f169y6jaJ1irZBZKqAAGD4cLnMq+TGc+Qcz8+Xk0YB7MxJzstRcpzNVcjabKYgLywsxKhRo/DII4+gSZMmAIDU1FS4u7vD39/fYN3AwECkpqbq17k1iYueL3ruXuvodDrk5OTcEUtubi50Op3Bjeh2Y8fKnzPj4oCdO1VHY/tsKccB8+f59u3AxYtAhQqyExiRs3GUHP/rL9lHyNUVePJJo/6EqNRspiAfMWIEjhw5gpUrV6oOBTNmzICfn5/+Vr16ddUhkQ0KDAQiIuTyN9+ojcUe2FKOA+bP823b5L9PPAF4eJghQCI74yg5fvy4/LdePTnkIZE12ERB/tprr2Hjxo3Ytm0bqlWrpn88KCgIN27cQEZGhsH6aWlpCAoK0q9ze2/tovv3W0er1cLT0/OOeCIjI5GZmam/nT9/vtTvkRzTCy/If3/4Abh+XW0stszWchwwf55v3y7/ffTRUm2GyC45Uo6fOiX/rVvXqNWJzEJpQS6EwGuvvYY1a9YgOjoatWvXNni+VatWcHNzw9atW/WPJSYmIjk5GSH/DvAbEhKCw4cPIz09Xb9OVFQUtFotGjVqpF/n1m0UrRNyl0GCPTw8oNVqDW5ExXn0UaBaNSAzE9i0SXU0tsdWcxwwb57n5AB79shlFuTkTBwxx0+flv8+8IBRqxOZh6V7mN7Lq6++Kvz8/ERMTIxISUnR365du6ZfZ9iwYaJGjRoiOjpa7Nu3T4SEhIiQkBD98/n5+aJJkyaia9euIj4+XmzZskVUrFhRREZG6tc5c+aM8PLyEuPGjRMJCQliwYIFwtXVVWzZssWoOB2xJzyZz1tvyd74zzxj3u06wnFnLzkuROn2d3S0PAaqVBGisNDkPycnxRyXbC3He/aU+fzJJybsCKJimJLjSgtyAMXeli5dql8nJydHDB8+XJQrV054eXmJnj17ipSUFIPtnD17VnTr1k14enqKgIAAMXbsWJGXl2ewzrZt20Tz5s2Fu7u7qFOnjsFr3I8jfGiS5Rw+LD+83dyE+Ocf823XEY47e8lxIUq3vydPlsdARITJf0pOjDl+ky3leHCwzOfNm03aPNEdTMlxjRAcsO1+dDod/Pz8kJmZyeYrVKwWLeREEosWAcOGmWebPO6sqzT7u1MnICYGWLwYeOUVy8RHjoc5bl3G7G8hAF9fIDtbTvRVv76VgySHYkqO20SnTiJ7V9S5c8UKtXGQ9eXmAn/8IZfZfpzIvqWlyWLcxQWoVUt1NORMWJATmUFEBKDRALt2AWfOqI6GrGnPHjnCTmAg0KCB6miIqDSKOnTWqAG4u6uNhZwLC3IiM6hSBejcWS5/+63aWMi6ioY77NhRfikjIvtVNOQhR1gha2NBTmQm/fvLf1eskO0QyTlw/HEix8EhD0kVFuREZtKzJ+DpCZw4Aezdqzoasoa8PGD3brnMgpzI/rEgJ1VYkBOZia8v8MwzcpmdO53Dvn3AtWtAhQrAv/OXEJEd4yydpAoLciIzKmq2snKlvHpKju3W9uMu/DQlsnu8Qk6q8BRCZEZdugAVKwJ//w1ERamOhiyN7ceJHEdGBnDpklyuU0dpKOSEWJATmVGZMnIIRAD45hu1sZBl5ecDO3fKZRbkRPav6Op4YKBsgkhkTSzIicysqNnK2rWATqc0FLKgAweArCzA3x8IDlYdDRGVFpurkEosyInMrFUrOUHM9evAmjWqoyFLKWqu0qED4OqqNhYiKj0W5KQSC3IiM9NogBdekMtstuK42H6cyLFwUiBSiQU5kQX06yf/jY4GLlxQGwuZX0EB8Pvvcjk0VGkoRGQmRVfIOeQhqcCCnMgCatcG2reXM3Z+/73qaMjcDh0CMjMBrRZo3lx1NERkDmyyQiqxICeyEDZbcVxFzVXat2f7cSJHkJMD/PWXXGZBTiqwICeykN69AXd3eTX10CHV0ZA5sf04kWNJSpL/arVAQIDaWMg5sSAnspDy5YHu3eXyt9+qjYXMp7AQ2LFDLrMgJ3IMtzZX0WjUxkLOiQU5kQUVNVv59lvZEZDs39GjwOXLgLc30LKl6miIyBw4wgqpxoKcyIKeeEJOHHPhws1mDmTfYmLkv488Ari5KQ2FiMyEI6yQaizIiSyobFnZlhwAVq5UGwuZx86d8l82VyFyHEVtyOvUURsHOS8W5EQW1qOH/LeokCP7dvy4/JfDHRI5juxs+a+fn9o4yHmxICeysHbt5L8JCcCVK2pjodIRgmMVEzmioj4+HMaUVGFBTmRhFSvebJcYF6c2Fiqd9HR5JU2jAWrVUh0NEZkLC3JSjQU5kRWEhMh/Y2PVxkGlU3R1vHp1wMNDbSxEZD4syEk1FuREVlDUbIUFuX1jcxUix8SCnFRjQU5kBUVXyOPi5MQyZJ9YkBM5JhbkpBoLciIrCA6WE8nodLJzJ9knFuREjokFOanGgpzICsqUAdq0kctstmK/WJATOSYW5KSa0oJ8x44d6NGjB6pUqQKNRoO1a9caPC+EwKRJk1C5cmV4enoiLCwMJ0+eNFjn8uXL6NevH7RaLfz9/fHSSy8hKyvLYJ1Dhw6hQ4cOKFu2LKpXr45Zs2ZZ+q0R3cEZO3Y6Wo6zICe6kyPkOQtyUk1pQZ6dnY1mzZphwYIFxT4/a9YsfPLJJ1i8eDHi4uLg7e2N8PBwXL9+Xb9Ov379cPToUURFRWHjxo3YsWMHhg4dqn9ep9Oha9euqFmzJvbv34/Zs2djypQpWLJkicXfH9GtnLEgd6Qcv3pVDnsIsCAnupUj5DkLclJO2AgAYs2aNfr7hYWFIigoSMyePVv/WEZGhvDw8BDff/+9EEKIY8eOCQBi7969+nV+/vlnodFoxIULF4QQQixcuFCUK1dO5Obm6td56623RIMGDYyOLTMzUwAQmZmZJX17RCI9XQg5tYwQly/ff31HO+5sOceFuP/+jo+X/3cVKpi0WaK7crQcF8K28/xe+7t2bZnfsbFGb47ovkzJcZttQ56UlITU1FSEhYXpH/Pz80Pbtm0R++8lxtjYWPj7+6N169b6dcLCwuDi4oK4f2dgiY2NRceOHeHu7q5fJzw8HImJibjCaRPJijhBkCF7y3E2VyEynb3kedEVchebrYrI0dnsoZeamgoACAwMNHg8MDBQ/1xqaioqVapk8HyZMmVQvnx5g3WK28atr3G73Nxc6HQ6gxuROXA88ptU5jhgep6zICcynb2cy9lkhVSz2YJcpRkzZsDPz09/q169uuqQyEE4YztyW2VqnrMgJ7IvpuQ4C3JSzWYL8qCgIABAWlqaweNpaWn654KCgpBe1MvqX/n5+bh8+bLBOsVt49bXuF1kZCQyMzP1t/Pnz5f+DRGBEwTdSmWOA6bnOQtyItPZy7m86POYBTmpYrMFee3atREUFIStW7fqH9PpdIiLi0PIv1VNSEgIMjIysH//fv060dHRKCwsRNu2bfXr7NixA3l5efp1oqKi0KBBA5QrV67Y1/bw8IBWqzW4EZkDJwi6SWWOA6bnOQtyItPZy7mcV8hJOSt0Mr2rq1evigMHDogDBw4IAGLOnDniwIED4ty5c0IIIWbOnCn8/f3FunXrxKFDh8TTTz8tateuLXJycvTbePzxx0WLFi1EXFyc2Llzp6hXr56IiIjQP5+RkSECAwNF//79xZEjR8TKlSuFl5eX+Oyzz4yO0xF7wpM6oaGyN//nn997PUc47uwlx4W49/6+cUMIV1f5//bvoA9EpeYIOS6E/eT5vfa3v7/M7+PHS7EjiG5jSo4rLci3bdsmANxxGzhwoBBCDpc0ceJEERgYKDw8PETnzp1FYmKiwTYuXbokIiIihI+Pj9BqtWLw4MHi6tWrBuscPHhQtG/fXnh4eIiqVauKmTNnmhSno3xokm2IjJQf/C++eO/1HOG4s5ccF+Le+/vkSfl/5ukpRGGhyZsmKpYj5LgQ9pPn99rfvr4yx0+eNO29E92LKTmuEUIIa1yJt2c6nQ5+fn7IzMxk8xUqtQ0bgKeeAh58EDh27O7r8bizrnvt719+AR5/HGjcGDhyRFGA5HCY49Z1r/3t7Q1cuwacOQPUrq0oQHI4puS4zbYhJ3JURUMfJiQAHArfPrD9OJFjYxtyUo0FOZGVVax4s7DjBEH2gQU5kWNjQU6qsSAnUoDjkdsXFuREjo0FOanGgpxIARbk9oUFOZHjkt055TILclKFBTmRApwgyH4IITt6ASzIiRxR0dVxgAU5qcOCnEgBThBkP1JT5egLLi5AzZqqoyEic2NBTraABTmRAmXKAG3ayGU2W7FtRc1VatQA3N3VxkJE5seCnGxBGdUBEDmrkBDg1CnglpmgyQbl5MimKvXqqY6EiCxBCKB6ddl8kAU5qcKCnEiRqVOB6dNVR0H306WL/OLEKdSIHJO3N5CcrDoKcnZsskKkiJub6gjIFBqN6giIiMhRsSAnIiIiIlKIBTkRERERkUIsyImIiIiIFGJBTkRERESkEAtyIiIiIiKFWJATERERESnEcciNIP4dgFin0ymOhJxJ0fEmOAC2VTDPydqY49bFHCdrMyXHWZAb4erVqwCA6tWrK46EnNHVq1fh5+enOgyHxzwnVZjj1sEcJ1WMyXGN4Ffz+yosLMTFixfh6+sLTTGzg+h0OlSvXh3nz5+HVqtVEKHzcuR9L4TA1atXUaVKFbi4sHWZpd0tzx35GFOJ+5U5bm08l1ufs+9TU3KcV8iN4OLigmrVqt13Pa1W65QHnC1w1H3Pq2bWc788d9RjTDVn36/McevhuVwdZ96nxuY4v5ITERERESnEgpyIiIiISCEW5Gbg4eGByZMnw8PDQ3UoTof7niyNx5hlcL+SreExaX7cp8Zjp04iIiIiIoV4hZyIiIiISCEW5ERERERECrEgJyIiIiJSiAU5EREREZFCLMhLacGCBahVqxbKli2Ltm3bYs+ePapDsjs7duxAjx49UKVKFWg0Gqxdu9bgeSEEJk2ahMqVK8PT0xNhYWE4efKkwTqXL19Gv379oNVq4e/vj5deeglZWVkG6xw6dAgdOnRA2bJlUb16dcyaNcvSb43sHPPbMu6X80QqMN/NZ8qUKdBoNAa3hg0bqg7LprEgL4X//e9/GDNmDCZPnow///wTzZo1Q3h4ONLT01WHZleys7PRrFkzLFiwoNjnZ82ahU8++QSLFy9GXFwcvL29ER4ejuvXr+vX6devH44ePYqoqChs3LgRO3bswNChQ/XP63Q6dO3aFTVr1sT+/fsxe/ZsTJkyBUuWLLH4+yP7xPy2nPvlPJG1Md/Nr3HjxkhJSdHfdu7cqTok2yaoxB566CExYsQI/f2CggJRpUoVMWPGDIVR2TcAYs2aNfr7hYWFIigoSMyePVv/WEZGhvDw8BDff/+9EEKIY8eOCQBi7969+nV+/vlnodFoxIULF4QQQixcuFCUK1dO5Obm6td56623RIMGDSz8jsheMb+t4/acJ1KB+W5ekydPFs2aNVMdhl3hFfISunHjBvbv34+wsDD9Yy4uLggLC0NsbKzCyBxLUlISUlNTDfazn58f2rZtq9/PsbGx8Pf3R+vWrfXrhIWFwcXFBXFxcfp1OnbsCHd3d/064eHhSExMxJUrV6z0bsheML+JnAfz3TJOnjyJKlWqoE6dOujXrx+Sk5NVh2TTWJCX0D///IOCggIEBgYaPB4YGIjU1FRFUTmeon15r/2cmpqKSpUqGTxfpkwZlC9f3mCd4rZx62sQFWF+EzkP5rv5tW3bFsuWLcOWLVuwaNEiJCUloUOHDrh69arq0GxWGdUBEBEREZHj6Natm365adOmaNu2LWrWrIkffvgBL730ksLIbBevkJdQQEAAXF1dkZaWZvB4WloagoKCFEXleIr25b32c1BQ0B0db/Lz83H58mWDdYrbxq2vQVSE+U3kPJjvlufv74/69evj1KlTqkOxWSzIS8jd3R2tWrXC1q1b9Y8VFhZi69atCAkJURiZY6lduzaCgoIM9rNOp0NcXJx+P4eEhCAjIwP79+/XrxMdHY3CwkK0bdtWv86OHTuQl5enXycqKgoNGjRAuXLlrPRuyF4wv4mcB/Pd8rKysnD69GlUrlxZdSi2S3WvUnu2cuVK4eHhIZYtWyaOHTsmhg4dKvz9/UVqaqrq0OzK1atXxYEDB8SBAwcEADFnzhxx4MABce7cOSGEEDNnzhT+/v5i3bp14tChQ+Lpp58WtWvXFjk5OfptPP7446JFixYiLi5O7Ny5U9SrV09ERETon8/IyBCBgYGif//+4siRI2LlypXCy8tLfPbZZ1Z/v2QfmN+Wc7+cJ7I25rt5jR07VsTExIikpCSxa9cuERYWJgICAkR6errq0GwWC/JSmj9/vqhRo4Zwd3cXDz30kPjjjz9Uh2R3tm3bJgDccRs4cKAQQg59OHHiRBEYGCg8PDxE586dRWJiosE2Ll26JCIiIoSPj4/QarVi8ODB4urVqwbrHDx4ULRv3154eHiIqlWripkzZ1rrLZKdYn5bxv1ynkgF5rv59O3bV1SuXFm4u7uLqlWrir59+4pTp06pDsumaYQQQs21eSIiIiIiYhtyIiIiIiKFWJATERERESnEgpyIiIiISCEW5ERERERECrEgJyIiIiJSiAU5EREREZFCLMiJiIiIiBRiQU5EREREDmHHjh3o0aMHqlSpAo1Gg7Vr11r09WrVqgWNRnPHbcSIESZthwU5ERERETmE7OxsNGvWDAsWLLDK6+3duxcpKSn6W1RUFACgd+/eJm2HBTk5nIyMDLRu3RrNmzdHkyZN8Pnnn6sOiYiIiKygW7dueO+999CzZ89in8/NzcWbb76JqlWrwtvbG23btkVMTEyJX69ixYoICgrS3zZu3IgHHngAjz76qEnbYUFODsfX1xc7duxAfHw84uLiMH36dFy6dEl1WEQ2ITQ0VP+Tanx8vNm3PWrUqFKvYy6DBg3Sv1dL/2xNZE8s+TlgLpbK39deew2xsbFYuXIlDh06hN69e+Pxxx/HyZMnS73tGzduYMWKFXjxxReh0WhM+lsW5HSH4hLVmifR0nJ1dYWXlxcA+U1YCAEhhP55nqTJ2b388stISUlBkyZNzLrd1atX491339XfV/25MW/ePKSkpCh7fSJbdvvnQGnbXg8ePBjvvPOO2eKzRP4mJydj6dKlWLVqFTp06IAHHngAb775Jtq3b4+lS5eWevtr165FRkYGBg0aZPLfsiCnYlnqhH035k7kjIwMNGvWDNWqVcO4ceMQEBCgf44naXJ2Xl5eCAoKQpkyZcy63fLly8PX19es2ywNPz8/BAUFqQ6DyCbd/jlQmrbXBQUF2LhxI5566imzxWeJ/D18+DAKCgpQv359+Pj46G/bt2/H6dOnAQDHjx8vtpPmrbcJEyYUu/0vv/wS3bp1Q5UqVUyOjQU5FctSJ+ziWCKR/f39cfDgQSQlJeG7775DWlqa/jmepMneVKtWDQsXLjR4bPfu3fDy8sK5c+dKvf1atWph7ty5Bo81b94cU6ZM0d8PDQ3FG2+8gfHjx6N8+fIICgoyeL5onaIr4oMGDcL27dsxb948/Uns7Nmzd7x2YWEhZsyYgdq1a8PT0xPNmjXDjz/+aLDOjz/+iODgYHh6eqJChQoICwtDdnb2fZ8jciSW/hy4X9vre9m9ezfc3NzQpk2bYp8PDQ3F66+/jlGjRqFcuXIIDAzE559/juzsbAwePBi+vr6oW7cufv7559K+jXvKysqCq6sr9u/fj/j4eP0tISEB8+bNAwDUqVMHCQkJ97yNHTv2jm2fO3cOv/32G4YMGVKi2FiQOwhLJ+rtNm3aBD8/P3z77bcAgKtXr6Jfv37w9vZG5cqV8fHHHxv9c7UlEzkwMBDNmjXD77//Xqr3S6RS27ZtsXfvXv19IQRGjRqF0aNHo2bNmlaLY/ny5fD29kZcXBxmzZqFadOm6UcUuN28efMQEhKi/7UtJSUF1atXv2O9GTNm4Ouvv8bixYtx9OhRjB49Gi+88AK2b98OAEhJSUFERARefPFFJCQkICYmBs8++yyEEPd8jsjR2MrnQHHWr1+PHj163LPd9PLlyxEQEIA9e/bg9ddfx6uvvorevXvj4Ycfxp9//omuXbuif//+uHbtmsXibNGiBQoKCpCeno66desa3Iou1Lm7u6Nhw4b3vFWsWPGObS9duhSVKlVC9+7dSxQbC3IHYc1E/e677xAREYFvv/0W/fr1AwCMGTMGu3btwvr16xEVFYXff/8df/75p1HbM3cip6Wl4erVqwCAzMxM7NixAw0aNCj9GydSpF27dgb5/c033+D8+fOIjIwEIK9wN23aFM2bN0enTp0sFkfTpk0xefJk1KtXDwMGDEDr1q2xdevWYtf18/ODu7u7/te2oKAguLq6GqyTm5uL6dOn46uvvkJ4eDjq1KmDQYMG4YUXXsBnn30GQBbk+fn5ePbZZ1GrVi0EBwdj+PDh8PHxuedzRI7mfp8DSUlJ6NSpExo1aoTg4GCr/lK0bt26+/7K3axZM7zzzjuoV68eIiMjUbZsWQQEBODll19GvXr1MGnSJFy6dAmHDh0qVSxZWVn6K9+A3C/x8fFITk5G/fr10a9fPwwYMACrV69GUlIS9uzZgxkzZmDTpk0lfs3CwkIsXboUAwcOLHHLAhbkDuJ+iWouCxYswPDhw7FhwwY8+eSTAOTV8eXLl+PDDz9E586d0aRJEyxduhQFBQVGbdPciXzu3Dl06NABzZo1Q4cOHfD6668jODi41O+dSJV27dohISEBWVlZyM7Oxttvv4333nvPoPDcvXs34uPjsW3bNovF0bRpU4P7lStXRnp6eom3d+rUKVy7dg1dunQxaM/59ddf69tzNmvWDJ07d0ZwcDB69+6Nzz//HFeuXLnvc0SO5n6fA4MGDcK0adNw7NgxbN++HR4eHlaJKyEhARcvXkTnzp3vud6tnx+urq6oUKGCwbk5MDAQAEr1mQIA+/btQ4sWLdCiRQsA8oJhixYtMGnSJADySvaAAQMwduxYNGjQAM888wz27t2LGjVqlPg1f/vtNyQnJ+PFF18s8TYs30CYrKJdu3aYMGECsrKyoNFoij1hX7t2DQ8++CB69+6NDz/80OTX+PHHH5Geno5du3YZNC85c+YM8vLy8NBDD+kf8/PzM+qqtCUS+amnnrLZYZyISqJVq1ZwcXHBn3/+id9++w0VK1bE4MGDzbZ9FxeXO5p55OXl3bGem5ubwX2NRoPCwsISv25WVhYA2QSuatWqBs8VFROurq6IiorC7t278euvv2L+/Pn473//i7i4ONSuXfuezxE5knt9Dhw9ehRubm7o0KEDANnB2lrWr1+PLl26oGzZsvdcr7jPj1sfK/qVvDSfKYBs5nqvZmtubm6YOnUqpk6dWqrXuVXXrl1L3VSOV8gdxK2J+sEHHxR7wn7//ffRrl27Er9GixYtULFiRXz11Vdma6Npa4lMZIu8vLwQHByMn376CR9++CE+/vhjuLjc/PjWaDR49NFH0aZNG32/DlNUrFjRYOQhnU6HpKSkUsft7u5+z1/KGjVqBA8PDyQnJ9/RnvPW9uYajQaPPPIIpk6digMHDsDd3R1r1qy573NEjuRenwMnT56Ej48PevTogZYtW2L69OlWi2vdunV4+umnrfZ6jopXyB3ErYn6+eefY/PmzQYn7JMnT+L48ePo0aMHjhw5UqLXeOCBB/DRRx8hNDQUrq6u+PTTTwHIHslubm4GP/lkZmbixIkT6Nix4z23uW7dOgwdOrRE8RA5k3bt2mH+/Pl4+umnERoaavDczp07UbVqVaSkpCAsLAzBwcF3NC+5l8ceewzLli1Djx494O/vj0mTJt3R3rskatWqhbi4OJw9exY+Pj53XLXz9fXFm2++idGjR6OwsBDt27dHZmYmdu3aBa1Wi4EDByIuLg5bt25F165dUalSJcTFxeHvv//Ggw8+eM/niBzR3T4H8vPz8fvvvyM+Ph6VKlXC448/jjZt2qBLly5GbzsrKwunTp3S3y9qe12+fPm7NudIT0/Hvn37sH79+hK/J5JYkDuQe52w33zzTcyePRu7d+8u1WvUr18f27ZtQ2hoKMqUKYO5c+fC19cXAwcOxLhx41C+fHlUqlQJkydPhouLyz07ajKRiYzXrFkzuLm5Yfbs2Xc8V9Tco3LlynjiiSfw559/mlSQR0ZGIikpCU8++ST8/Pzw7rvvmuUK+ZtvvomBAweiUaNGyMnJKXab7777LipWrIgZM2bgzJkz8Pf3R8uWLfH2228DALRaLXbs2IG5c+dCp9OhZs2a+Oijj9CtWzckJCTc9TkiR3S3z4GqVauidevW+l+WnnjiCcTHx5tUkO/bt8+gU/iYMWMAAAMHDsSyZcuK/ZsNGzbgoYceMpjrg0qGBbkDuVuirlu3DvXr10f9+vVLXZADQIMGDRAdHa2/Uv7RRx9hzpw5GDZsGJ588klotVqMHz8e58+fv2dTFCYykfFWrlyJ1157DXXr1jV4PDs7G4WFhfD19UVWVhaio6PRp08fk7at1WqxcuVKg8cGDhxocD8mJuaOv7t9Jr/b16lfvz5iY2PvuY5Go8HIkSMxcuTIYmN78MEHsWXLFpOfI3JEd/scaNOmDdLT03HlyhX4+flhx44deOWVV0za9v3aXhfHmEEZgOI/P4qbl8CZhyxlQe5A7paof/zxB1auXIlVq1YhKysLeXl50Gq1+h7Hxrg9mR588EGDyXZ8fX0N2q5mZ2dj6tSp92yOwkQmurfCwkL8/fff+PLLL3Hy5EmsW7fujnXS0tL0E3kUFBTg5ZdfvuuY/kUWLlyIL774ArGxsQ49AtGwYcOwYsUK1WEQlYoxnwNlypTB9OnT0bFjRwgh0LVrV/1IaHdjjs+B9u3bIyIiokR/ez/Olr8awSrGrt2aqAsXLsSxY8eg1Wrvuv6yZctw5MiRe46yEhoait27d8Pd3d3oRD1w4ACOHz+Ohx56CJmZmZg2bRpiYmJw6tSpu14BnzVrFiIiIoqdLMSSipI8Ozsba9aswTPPPGPV1ycyVkxMDB577DE0bNgQS5cuRdu2bUu9zQsXLiAnJwcAUKNGDbi7u5d6m7YqPT0dOp0OgGzO4+3trTgiItM56+eAs+UvC3I7Z2qiGlOQlyRRDxw4gCFDhiAxMRHu7u5o1aoV5syZY5NX35wtyYmIiMi2sSAnIiIiIlKI45ATERERESnEgpyIiIiISCEW5ERERERECrEgJyIiIiJSiAU5EREREZFCLMiJiIiIiBRiQU5EREREpBALciIiIiIihViQExEREREpxIKciIiIiEghFuRERERERAqxICciIiIiUogFORERERGRQizIiYiIiIgUYkFORERERKQQC3IiIiIiIoVYkBMRERERKcSCnIiIiIhIIRbkREREREQKsSAnIiIiIlKIBTkRERERkUIsyImIiIiIFGJBTkRERESkEAtyIiIiIiKFWJATERERESnEgpyIiIiISCEW5ERERERECrEgJyIiIiJSiAU5EREREZFCLMiJiIiIiBRiQU5EREREpBALciIiIiIihViQExEREREpxIKciIiIiEghFuRERERERAqxICciIiIiUogFORERERGRQizIiYiIiIgUYkFORERERKQQC3IiIiIiIoVYkBMRERERKcSCnIiIiIhIIRbkREREREQKsSAnIiIiIlKIBTkRERERkUIsyImIiIiIFGJBTkRERESkEAtyIiIiIiKFWJATERERESnEgpyIiIiISCEW5ERERERECrEgJyIiIiJSiAU5EREREZFCLMiJiIiIiBRiQU5EREREpBALciIiIiIihViQExEREREpVEZ1AERE5BgKCgqQl5enOgyb4ebmBldXV9VhEJEdYEFORESlIoRAamoqMjIyVIdic/z9/REUFASNRqM6FCKyYSzIiYioVIqK8UqVKsHLy4vFJ+SXlGvXriE9PR0AULlyZcUREZEtY0FOREQlVlBQoC/GK1SooDocm+Lp6QkASE9PR6VKldh8hYjuip06iYioxIrajHt5eSmOxDYV7Re2rSeie2FBTkREpcZmKsXjfiEiY7AgJyIiIiJSiAU5EREREZFCLMiJiIiMIITApEmTULlyZXh6eiIsLAwnT55UHRYROQAW5EREREaYNWsWPvnkEyxevBhxcXHw9vZGeHg4rl+/rjo0IrJzLMiJiMjpfP3116hQoQJyc3MNHn/mmWfQv3//O9YXQmDu3Ll455138PTTT6Np06b4+uuvcfHiRaxdu9ZKURORo2JBTkREZiUEkJ2t5iaEcTH27t0bBQUFWL9+vf6x9PR0bNq0CS+++OId6yclJSE1NRVhYWH6x/z8/NC2bVvExsaWep8RkXPjxEBERGRW164BPj5qXjsrC/D2vv96np6eeP7557F06VL07t0bALBixQrUqFEDoaGhd6yfmpoKAAgMDDR4PDAwUP8cEVFJ8Qo5ERE5pZdffhm//vorLly4AABYtmwZBg0ahO+++w4+Pj762++//644UiJydLxCTkREZuXlJa9Uq3ptY7Vo0QLNmjXD119/ja5du+Lo0aPYtGkT/P390bZtW/16VatWRUpKCgAgLS0NlStX1j+XlpaG5s2bmyt8InJSLMiJiMisNBrjmo3YgiFDhmDu3Lm4cOECwsLCUL16dQCAr6+vwXq1a9dGUFAQtm7dqi/AdTod4uLi8Oqrr1o7bCJyMGyyQkRETuv555/HX3/9hc8//7zYzpxFNBoNRo0ahffeew/r16/H4cOHMWDAAFSpUgXPPPOM9QImIofEK+REROS0/Pz80KtXL2zatOm+hfX48eORnZ2NoUOHIiMjA+3bt8eWLVtQtmxZ6wRLRA6LBTkRETm1CxcuoF+/fvDw8LjnehqNBtOmTcO0adOsFBkROQsW5ERE5JSuXLmCmJgYxMTEYOHCharDISInxoKciIicUosWLXDlyhV88MEHaNCggepwiMiJsSAnIiKndPbsWdUhEBEB4CgrRERERERKsSAnIqJSE0KoDsEmcb8QkTFYkBMRUYm5ubkBAK5du6Y4EttUtF+K9hMRUXHYhpyIiErM1dUV/v7+SE9PBwB4eXlBo9Eojko9IQSuXbuG9PR0+Pv7w9XVVXVIRGTDNIK/pxERUSkIIZCamoqMjAzVodgcf39/BAUF8UsKEd0TC3IiIjKLgoIC5OXlqQ7DZri5ufHKOBEZhQU5EREREZFC7NRJRERERKQQC3IiIiIiIoVYkBMRERERKcSCnIiIiIhIIRbkREREREQKsSAnIiIiIlKIBTkRERERkUL/B8Feuu5uxN54AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "y1; top = (0.14342367312673812-0.00687389024121427j)\n", + "y2; top = (-1.199040866595169e-12+9.367506770274758e-14j)\n", + "y3; top = (0.022959291124553246-0.0005423192643148481j)\n", + "y4; top = (-2.9976021664879227e-13-2.8102520310824275e-14j)\n", + "y5; top = (1.6584764638741158-0.030564212547253417j)\n", + "y6; top = (8.333333333333333e-07-4.625929269271486e-24j)\n", + "[0.65847646-0.03056421j 1.29915265-0.06226471j 0.20796862-0.00491241j]\n" ] } ], @@ -59,7 +79,7 @@ "def test_1layer():\n", " frequency = np.pi * 2. / (86400. * 7.5)\n", " \n", - " radius_array = np.linspace(0., 6000.0e3, 100)\n", + " radius_array = np.linspace(0., 6000.0e3, 10)\n", " indices_by_layer = radius_array > 0.\n", " density_array = np.ones_like(radius_array) * 5400.\n", " planet_bulk_density = np.average(density_array)\n", @@ -91,7 +111,7 @@ " solve_for = None,\n", " core_condition = 0,\n", " use_kamata = True,\n", - " starting_radius = 1.0e2,\n", + " starting_radius = 0.1,\n", " start_radius_tolerance = 1.0e-5,\n", " integration_method = 'RK45',\n", " integration_rtol = 1.0e-9,\n", @@ -129,7 +149,12 @@ " print(f'Shape: {ys.shape}.')\n", " yplot([ys], [radius_array], colors=['b'])\n", " \n", - " print(ys[1, -1])\n", + " print(\"y1; top =\", ys[0, -1])\n", + " print(\"y2; top =\", ys[1, -1])\n", + " print(\"y3; top =\", ys[2, -1])\n", + " print(\"y4; top =\", ys[3, -1])\n", + " print(\"y5; top =\", ys[4, -1])\n", + " print(\"y6; top =\", ys[5, -1])\n", "\n", " print(solution.love)\n", "\n", @@ -147,6 +172,31 @@ "solution = test_1layer()" ] }, + { + "cell_type": "code", + "execution_count": 5, + "id": "0314d5db-ecc5-4f26-a3f5-30a8228bad35", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ nan +nanj, 0.02363794-5.10063889e+01j,\n", + " 0.04704499-1.01199241e+02j, 0.06988046-1.49431286e+02j,\n", + " 0.09159755-1.93958406e+02j, 0.75661593-4.30186856e-02j,\n", + " 0.86999229-4.83800359e-02j, 0.95242984-5.12998238e-02j,\n", + " 0.99186257-5.10068925e-02j, 0.97436449-4.66985289e-02j])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solution.result[0]" + ] + }, { "cell_type": "markdown", "id": "54f7ea0e-0472-49e7-8899-79ed691f47b2", diff --git a/MANIFEST.in b/MANIFEST.in index 1bdba99e..e889f561 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -20,3 +20,4 @@ include TidalPy/RadialSolver/rs_solution_.cpp include TidalPy/RadialSolver/rs_solution_.hpp include TidalPy/Material/eos/eos_solution_.cpp include TidalPy/Material/eos/eos_solution_.hpp +include TidalPy/utilities/dimensions/nondimensional_.hpp \ No newline at end of file diff --git a/Tests/Test_Utilities/Test_Dimensions/test_nondimensional.py b/Tests/Test_Utilities/Test_Dimensions/test_nondimensional.py index b83d5cdc..4795d60c 100644 --- a/Tests/Test_Utilities/Test_Dimensions/test_nondimensional.py +++ b/Tests/Test_Utilities/Test_Dimensions/test_nondimensional.py @@ -5,62 +5,52 @@ import TidalPy TidalPy.test_mode = True -from TidalPy.constants import G -from TidalPy.utilities.dimensions.nondimensional import non_dimensionalize_physicals, redimensionalize_physicals +from math import isnan, isclose, sqrt +from TidalPy.constants import G, PI_DBL +from TidalPy.utilities.dimensions.nondimensional import NonDimensionalScalesClass, build_nondimensional_scales -def test_non_dimensionalize_physicals(): - frequency = 1.0e-3 - mean_radius = 1.0e6 - bulk_density = 5500. - radius_array = np.linspace(0., mean_radius, 10) - density_array = np.linspace(7000, 3500, 10) - bulk_array = 100.0e9 * np.ones_like(radius_array) - gravity_array = np.linspace(0.1, 7.0, 10) - shear_array = (50.0e9 + 20.0e3j) * np.ones(radius_array.size, dtype=np.complex128) +frequency = 1.0e-3 +mean_radius = 1.0e6 +bulk_density = 5500. - # Make copies of the arrays so we can keep the original values for comparison. - radius_array_out = radius_array.copy() - density_array_out = density_array.copy() - bulk_array_out = bulk_array.copy() - gravity_array_out = gravity_array.copy() - shear_array_out = shear_array.copy() +def test_non_dimensionalize_structure(): + """ Test that the non-dimensionalize structure initializes correctly. """ - # non-dimensionalize the arrays and get some non-dim constants. - output_freq_nondim, G_nondim = non_dimensionalize_physicals( - frequency, mean_radius, bulk_density, - radius_array_out, density_array_out, gravity_array_out, bulk_array_out, shear_array_out) + test_struct = NonDimensionalScalesClass() - # Build conversions for comparison - second2_conversion = 1. / (np.pi * G * bulk_density) - second_conversion = np.sqrt(second2_conversion) - length_conversion = mean_radius - density_conversion = bulk_density - mass_conversion = bulk_density * mean_radius**3 - pascal_conversion = mass_conversion / (length_conversion * second2_conversion) - - # Check conversion is correct - assert np.allclose(radius_array / length_conversion, radius_array_out) - assert np.allclose(density_array / density_conversion, density_array_out) - assert np.allclose(gravity_array / (length_conversion / second2_conversion), gravity_array_out) - assert np.allclose(bulk_array / pascal_conversion, bulk_array_out) - assert np.allclose(shear_array / pascal_conversion, shear_array_out) - - assert np.isclose(G_nondim, G / (length_conversion**3 / (mass_conversion * second2_conversion))) - assert np.isclose(output_freq_nondim, frequency / (1. / second_conversion)) - - # Check that the redimensionalize function works too. - # Use the output from the non-dim as inputs to the redim. - output_freq_redim, G_redim = redimensionalize_physicals( - frequency, mean_radius, bulk_density, - radius_array_out, density_array_out, gravity_array_out, bulk_array_out, shear_array_out) - - # These should all match the original again. - assert np.allclose(radius_array, radius_array_out) - assert np.allclose(density_array, density_array_out) - assert np.allclose(gravity_array, gravity_array_out) - assert np.allclose(bulk_array, bulk_array_out) - assert np.allclose(shear_array, shear_array_out) - assert np.isclose(G_redim, G) - assert np.isclose(output_freq_redim, frequency) + # All conversions should initialize as nans + assert isnan(test_struct.second2_conversion) + assert isnan(test_struct.second_conversion) + assert isnan(test_struct.length_conversion) + assert isnan(test_struct.length3_conversion) + assert isnan(test_struct.density_conversion) + assert isnan(test_struct.mass_conversion) + assert isnan(test_struct.pascal_conversion) + + +def test_build_nondimensional_scales(): + """ Test building a non-dimensionalize structure with real inputs. """ + + # Build + test_struct = build_nondimensional_scales(frequency, mean_radius, bulk_density) + + # Check values were set correctly + assert not isnan(test_struct.second2_conversion) + assert not isnan(test_struct.second_conversion) + assert not isnan(test_struct.length_conversion) + assert not isnan(test_struct.length3_conversion) + assert not isnan(test_struct.density_conversion) + assert not isnan(test_struct.mass_conversion) + assert not isnan(test_struct.pascal_conversion) + + # Check that they are the correct values. + second2_conv = 1. / (PI_DBL * G * bulk_density) + assert isclose(test_struct.second2_conversion, second2_conv) + assert isclose(test_struct.second_conversion, second2_conv) + assert isclose(test_struct.length_conversion, mean_radius) + assert isclose(test_struct.length3_conversion, mean_radius**3) + assert isclose(test_struct.density_conversion, bulk_density) + assert isclose(test_struct.mass_conversion, bulk_density * mean_radius**3) + assert isclose(test_struct.pascal_conversion, (bulk_density * mean_radius**3) / (mean_radius * second2_conv)) diff --git a/TidalPy/Material/eos/eos_solution.pxd b/TidalPy/Material/eos/eos_solution.pxd index be59171c..ca3e607b 100644 --- a/TidalPy/Material/eos/eos_solution.pxd +++ b/TidalPy/Material/eos/eos_solution.pxd @@ -7,6 +7,11 @@ cnp.import_array() from CyRK cimport CySolverResult +from TidalPy.utilities.dimensions.nondimensional cimport NonDimensionalScalesCC + +cdef extern from "nondimensional_.hpp" nogil: + pass + cdef extern from "eos_solution_.cpp" nogil: const size_t EOS_Y_VALUES @@ -16,6 +21,8 @@ cdef extern from "eos_solution_.cpp" nogil: cdef cppclass EOSSolutionCC: int error_code int iterations + int nondim_status + int solution_nondim_status cpp_bool success cpp_bool max_iters_hit cpp_bool radius_array_set @@ -35,6 +42,13 @@ cdef extern from "eos_solution_.cpp" nogil: double mass double moi + double redim_length_scale + double redim_gravity_scale + double redim_mass_scale + double redim_density_scale + double redim_moi_scale + double redim_pascal_scale + vector[double] upper_radius_bylayer_vec vector[shared_ptr[CySolverResult]] cysolver_results_sptr_bylayer_vec @@ -64,14 +78,13 @@ cdef extern from "eos_solution_.cpp" nogil: const size_t layer_index, const double radius, double* y_interp_ptr) - void call_vectorize( - const size_t layer_index, - const double* radius_array_ptr, - size_t len_radius_array, - double* y_interp_ptr) void change_radius_array( double* new_radius_ptr, size_t new_radius_size) void interpolate_full_planet() + + void dimensionalize_data( + NonDimensionalScalesCC* nondim_scales, + cpp_bool redimensionalize) diff --git a/TidalPy/Material/eos/eos_solution_.cpp b/TidalPy/Material/eos/eos_solution_.cpp index 4d45d9c4..54d723be 100644 --- a/TidalPy/Material/eos/eos_solution_.cpp +++ b/TidalPy/Material/eos/eos_solution_.cpp @@ -29,7 +29,7 @@ EOSSolutionCC::EOSSolutionCC( EOSSolutionCC::~EOSSolutionCC( ) { - printf("TidalPy::EOSSolutionCC Deconstructor Called.\n"); + printf("TidalPy::EOSSolutionCC Deconstructor Called; addr = %p.\n", this); // Reset each shared pointer in the cysolver vector for (size_t i = 0; i < this->cysolver_results_sptr_bylayer_vec.size(); i++) { @@ -63,28 +63,62 @@ void EOSSolutionCC::call( const double radius, double* y_interp_ptr) const { - if (layer_index < this->current_layers_saved) [[unlikely]] + if (layer_index < this->current_layers_saved) [[likely]] { this->cysolver_results_sptr_bylayer_vec[layer_index]->call(radius, y_interp_ptr); - } - else - { - // TODO: Better error handling - printf("Error! EOSSolution::call was asked to interpolate a layer that it has not saved."); - std::exception(); - } -} + if (this->nondim_status == 1) + { + // Acceleration due to Gravity + y_interp_ptr[0] *= this->redim_gravity_scale; -void EOSSolutionCC::call_vectorize( - const size_t layer_index, - const double* radius_array_ptr, - size_t len_radius_array, - double* y_interp_ptr) const -{ - if (layer_index < this->current_layers_saved) - { - this->cysolver_results_sptr_bylayer_vec[layer_index]->call_vectorize(radius_array_ptr, len_radius_array, y_interp_ptr); + // Pressure + y_interp_ptr[1] *= this->redim_pascal_scale; + + // These assume sphereical symmetry + // Total mass + y_interp_ptr[2] *= this->redim_mass_scale; + + // Moment of inertia (r^2 multipled by dm which was found above) + y_interp_ptr[3] *= this->redim_moi_scale; + + // Density + y_interp_ptr[4] *= this->redim_density_scale; + + // Shear modulus (real and complex) + y_interp_ptr[5] *= this->redim_pascal_scale; + y_interp_ptr[6] *= this->redim_pascal_scale; + + // Bulk modulus (real and complex) + y_interp_ptr[7] *= this->redim_pascal_scale; + y_interp_ptr[8] *= this->redim_pascal_scale; + } + else if (this->nondim_status == -1) + { + // Acceleration due to Gravity + y_interp_ptr[0] /= this->redim_gravity_scale; + + // Pressure + y_interp_ptr[1] /= this->redim_pascal_scale; + + // These assume sphereical symmetry + // Total mass + y_interp_ptr[2] /= this->redim_mass_scale; + + // Moment of inertia (r^2 multipled by dm which was found above) + y_interp_ptr[3] /= this->redim_moi_scale; + + // Density + y_interp_ptr[4] /= this->redim_density_scale; + + // Shear modulus (real and complex) + y_interp_ptr[5] /= this->redim_pascal_scale; + y_interp_ptr[6] /= this->redim_pascal_scale; + + // Bulk modulus (real and complex) + y_interp_ptr[7] /= this->redim_pascal_scale; + y_interp_ptr[8] /= this->redim_pascal_scale; + } } else { @@ -147,6 +181,8 @@ void EOSSolutionCC::change_radius_array( void EOSSolutionCC::interpolate_full_planet() { printf("TidalPy::EOSSolutionCC.interpolate_full_planet called.\n"); + this->solution_nondim_status = this->nondim_status; + if (this->current_layers_saved == 0) { // No layers have been saved, we can't perform the interpolation. @@ -208,4 +244,117 @@ void EOSSolutionCC::interpolate_full_planet() // Finished this->other_vecs_set = true; +} + +#include + +void EOSSolutionCC::dimensionalize_data(NonDimensionalScalesCC* nondim_scales, bool redimensionalize) +{ + + printf("EOSSolutionCC::dimensionalize_data Called.\n"); + printf("nondim scales (ptr = %p):\n", nondim_scales); + printf("\t length = %e\n", nondim_scales->length_conversion); + printf("\t leng3 = %e\n", nondim_scales->length3_conversion); + printf("\t densit = %e\n", nondim_scales->density_conversion); + printf("\t pascal = %e\n", nondim_scales->pascal_conversion); + printf("\t sec = %e\n", nondim_scales->second_conversion); + printf("\t sec2 = %e\n", nondim_scales->second2_conversion); + printf("\t mass = %e\n", nondim_scales->mass_conversion); + + // Save scalers + this->redim_length_scale = nondim_scales->length_conversion; + this->redim_gravity_scale = nondim_scales->length_conversion / nondim_scales->second2_conversion; + this->redim_mass_scale = nondim_scales->mass_conversion; + this->redim_density_scale = nondim_scales->density_conversion; + this->redim_moi_scale = nondim_scales->mass_conversion * nondim_scales->length_conversion * nondim_scales->length_conversion; + this->redim_pascal_scale = nondim_scales->pascal_conversion; + + // Figure out how to set flags used during eos solution calls + if (this->solution_nondim_status == 0) + { + // EOS was not explicitly non-dim'd or re-dim'd when solution was found. + if (redimensionalize) + { + // User is now "redimensionalizing." Thus we assume that the solution was non-dim'd + this->nondim_status = 1; + } + else + { + // User is now "dimensionalizing." Thus we assume that the solution was dim'd + this->nondim_status = -1; + } + } + else + { + // TODO: Deal with this case! For now push the problem to the user when they try to call... + } + + // Update other constants + if (redimensionalize) + { + this->pressure_error *= this->redim_pascal_scale; + } + else + { + this->pressure_error /= this->redim_pascal_scale; + } + + // Update layer data + for (size_t layer_i = 0; layer_i < this->num_layers; layer_i++) + { + if (redimensionalize) + { + this->upper_radius_bylayer_vec[layer_i] *= this->redim_length_scale; + } + else + { + this->upper_radius_bylayer_vec[layer_i] /= this->redim_length_scale; + } + } + + if (this->other_vecs_set) + { + // Work through arrays + for (size_t slice_i = 0; slice_i < this->radius_array_size; slice_i++) + { + if (redimensionalize) + { + this->radius_array_vec[slice_i] *= this->redim_length_scale; + this->gravity_array_vec[slice_i] *= this->redim_gravity_scale; + this->pressure_array_vec[slice_i] *= this->redim_pascal_scale; + this->mass_array_vec[slice_i] *= this->redim_mass_scale; + this->moi_array_vec[slice_i] *= this->redim_moi_scale; + this->density_array_vec[slice_i] *= this->redim_density_scale; + + // These complex arrays are stored as double arrays with twice the length (Cython and C++ don't play nicely with complex across all systems) + this->complex_shear_array_vec[2 * slice_i] *= this->redim_pascal_scale; + this->complex_shear_array_vec[2 * slice_i + 1] *= this->redim_pascal_scale; + this->complex_bulk_array_vec[2 * slice_i] *= this->redim_pascal_scale; + this->complex_bulk_array_vec[2 * slice_i + 1] *= this->redim_pascal_scale; + } + else + { + this->radius_array_vec[slice_i] /= this->redim_length_scale; + this->gravity_array_vec[slice_i] /= this->redim_gravity_scale; + this->pressure_array_vec[slice_i] /= this->redim_pascal_scale; + this->mass_array_vec[slice_i] /= this->redim_mass_scale; + this->moi_array_vec[slice_i] /= this->redim_moi_scale; + this->density_array_vec[slice_i] /= this->redim_density_scale; + + // These complex arrays are stored as double arrays with twice the length (Cython and C++ don't play nicely with complex across all systems) + this->complex_shear_array_vec[2 * slice_i] /= this->redim_pascal_scale; + this->complex_shear_array_vec[2 * slice_i + 1] /= this->redim_pascal_scale; + this->complex_bulk_array_vec[2 * slice_i] /= this->redim_pascal_scale; + this->complex_bulk_array_vec[2 * slice_i + 1] /= this->redim_pascal_scale; + } + } + + // Update global constants + this->radius = this->radius_array_vec[this->radius_array_size - 1]; + this->surface_gravity = this->gravity_array_vec[this->radius_array_size - 1]; + this->surface_pressure = this->pressure_array_vec[this->radius_array_size - 1]; + this->mass = this->mass_array_vec[this->radius_array_size - 1]; + this->moi = this->moi_array_vec[this->radius_array_size - 1]; + this->central_pressure = this->pressure_array_vec[0]; + } } \ No newline at end of file diff --git a/TidalPy/Material/eos/eos_solution_.hpp b/TidalPy/Material/eos/eos_solution_.hpp index 771ab485..b467cbc4 100644 --- a/TidalPy/Material/eos/eos_solution_.hpp +++ b/TidalPy/Material/eos/eos_solution_.hpp @@ -7,6 +7,8 @@ #include "cysolution.hpp" // Part of the CyRK python package. Header should be included in setup using CyRK.get_include() +#include "nondimensional_.hpp" // Part of the TidalPy.utilities module + static const size_t EOS_Y_VALUES = 4; static const size_t EOS_EXTRA_VALUES = 5; static const size_t EOS_DY_VALUES = 9; @@ -30,8 +32,10 @@ class EOSSolutionCC char message[256] = { }; public: - int iterations = -1; - int error_code = -100; + int iterations = -1; + int error_code = -100; + int nondim_status = 0; + int solution_nondim_status = 0; bool success = false; bool max_iters_hit = false; bool radius_array_set = false; @@ -49,6 +53,13 @@ class EOSSolutionCC double radius = NULL; double mass = NULL; double moi = NULL; + + double redim_length_scale = NULL; + double redim_gravity_scale = NULL; + double redim_mass_scale = NULL; + double redim_density_scale = NULL; + double redim_moi_scale = NULL; + double redim_pascal_scale = NULL; // Store results from CyRK's cysolve_ivp. std::vector upper_radius_bylayer_vec = std::vector(); @@ -103,4 +114,7 @@ class EOSSolutionCC // Run full planet interpolation through each layer. void interpolate_full_planet(); + void dimensionalize_data( + NonDimensionalScalesCC* nondim_scales, + bool redimensionalize); }; \ No newline at end of file diff --git a/TidalPy/Material/eos/solver.pyx b/TidalPy/Material/eos/solver.pyx index 31b01f98..83110fec 100644 --- a/TidalPy/Material/eos/solver.pyx +++ b/TidalPy/Material/eos/solver.pyx @@ -282,3 +282,5 @@ cdef void solve_eos( # Tell the eos solution to perform a full planet interpolation and store the results. Including surface results printf("\tDEBUG-solve_eos 10b3\n") eos_solution_ptr.interpolate_full_planet() + + printf("solve_eos Finished\n") diff --git a/TidalPy/RadialSolver/boundaries/surface_bc.pyx b/TidalPy/RadialSolver/boundaries/surface_bc.pyx index 38f95f99..786dcbf2 100644 --- a/TidalPy/RadialSolver/boundaries/surface_bc.pyx +++ b/TidalPy/RadialSolver/boundaries/surface_bc.pyx @@ -6,7 +6,8 @@ cimport numpy as np from libc.stdio cimport printf from libc.stdlib cimport exit -from libc.math cimport NAN + +from TidalPy.constants cimport d_NAN_DBL cdef void cf_get_surface_bc( @@ -35,7 +36,7 @@ cdef void cf_get_surface_bc( # Inititalize all boundary conditions to NaN # 15 = 5 (max_num_solutions) * 3 (number of surface conditions) for i in range(15): - boundary_conditions_ptr[i] = NAN + boundary_conditions_ptr[i] = d_NAN_DBL for j in range(num_bcs): if bc_model_ptr[j] == 0: diff --git a/TidalPy/RadialSolver/rs_solution.pxd b/TidalPy/RadialSolver/rs_solution.pxd index 10438fdb..ca477436 100644 --- a/TidalPy/RadialSolver/rs_solution.pxd +++ b/TidalPy/RadialSolver/rs_solution.pxd @@ -6,6 +6,7 @@ from libcpp cimport bool as cpp_bool from libcpp.vector cimport vector from libcpp.memory cimport shared_ptr, unique_ptr +from TidalPy.utilities.dimensions.nondimensional cimport NonDimensionalScalesCC from TidalPy.Material.eos.eos_solution cimport EOSSolutionCC # Need to include love_.cpp and eos_solution_.cpp in order to get solutions.cpp to see it and use it to link @@ -15,6 +16,9 @@ cdef extern from "love_.cpp" nogil: cdef extern from "eos_solution_.cpp" nogil: pass +cdef extern from "nondimensional_.hpp" nogil: + pass + cdef extern from "rs_solution_.cpp" nogil: const int MAX_NUM_Y @@ -49,6 +53,9 @@ cdef extern from "rs_solution_.cpp" nogil: cpp_bool array_changed) void set_message(const char* new_message) void find_love() + void dimensionalize_data( + NonDimensionalScalesCC* nondim_scales, + cpp_bool redimensionalize) cdef class RadialSolverSolution: diff --git a/TidalPy/RadialSolver/rs_solution.pyx b/TidalPy/RadialSolver/rs_solution.pyx index 483f08c0..0edc3efa 100644 --- a/TidalPy/RadialSolver/rs_solution.pyx +++ b/TidalPy/RadialSolver/rs_solution.pyx @@ -219,6 +219,16 @@ cdef class RadialSolverSolution: def eos_success(self): """ Return if the solver's equation of state sub-solver was successful """ return self.solution_storage_sptr.get().get_eos_solution_ptr().success + + @property + def eos_pressure_error(self): + """ Return the surface pressure error found by the equation of state sub-solver """ + return self.solution_storage_sptr.get().get_eos_solution_ptr().pressure_error + + @property + def eos_iterations(self): + """ Return the number of iterations performed by the EOS sub-solver to find convergence on surface pressure """ + return self.solution_storage_sptr.get().get_eos_solution_ptr().iterations @property def radius(self): diff --git a/TidalPy/RadialSolver/rs_solution_.cpp b/TidalPy/RadialSolver/rs_solution_.cpp index 6d5d3f69..0f13ea50 100644 --- a/TidalPy/RadialSolver/rs_solution_.cpp +++ b/TidalPy/RadialSolver/rs_solution_.cpp @@ -44,9 +44,11 @@ RadialSolutionStorageCC::RadialSolutionStorageCC( } } +#include + RadialSolutionStorageCC::~RadialSolutionStorageCC( ) { - + printf("RadialSolutionStorageCC Deconstructor Called; addr: %p.\n", this); } EOSSolutionCC* RadialSolutionStorageCC::get_eos_solution_ptr() @@ -128,3 +130,74 @@ void RadialSolutionStorageCC::find_love() // this->set_message("Can not update Love number values when solution is not complete or is unsuccessful.") } } + +void RadialSolutionStorageCC::dimensionalize_data( + NonDimensionalScalesCC* nondim_scales, + bool redimensionalize) +{ + + // Perform dimensionalization on the EOS solution first. + double* full_solution_ptr = this->full_solution_vec.data(); + EOSSolutionCC* eos_solution_ptr = this->get_eos_solution_ptr(); + eos_solution_ptr->dimensionalize_data(nondim_scales, redimensionalize); + + const double displacement_scale = (nondim_scales->second2_conversion / nondim_scales->length_conversion); + const double stress_scale = (nondim_scales->mass_conversion / nondim_scales->length3_conversion); + const double potential_scale = (1. / nondim_scales->length_conversion); + + // Redimensionalize the radial solutions + printf("RS SOLUTION:: REDIM CALLED\n"); + + printf("nondim scales (ptr = %p):\n", nondim_scales); + printf("\t length = %e\n", nondim_scales->length_conversion); + printf("\t leng3 = %e\n", nondim_scales->length3_conversion); + printf("\t densit = %e\n", nondim_scales->density_conversion); + printf("\t pascal = %e\n", nondim_scales->pascal_conversion); + printf("\t sec = %e\n", nondim_scales->second_conversion); + printf("\t sec2 = %e\n", nondim_scales->second2_conversion); + printf("\t mass = %e\n", nondim_scales->mass_conversion); + + printf("RS SOLUTION:: displacement scale = %e\n", displacement_scale); + printf("RS SOLUTION:: stress scale = %e\n", stress_scale); + printf("RS SOLUTION:: potential scale = %e\n", potential_scale); + if (this->success) + { + for (size_t solver_i = 0; solver_i < this->num_ytypes; solver_i++) + { + const size_t bc_stride = solver_i * MAX_NUM_Y_REAL; + printf("\t solver_i = %d; stride = %d\n", solver_i, bc_stride); + for (size_t slice_i = 0; slice_i < this->num_slices; slice_i++) + { + const size_t slice_stride = bc_stride + slice_i * MAX_NUM_Y_REAL * this->num_ytypes; + printf("\t\t slice = %d; slice stride = %d\n", slice_i, slice_stride); + // y1, y3 are the radial and tangential displacements with units of [s2 m-1] + // y2, y4 are the radial and tangential stresses with units of [kg m-3] + // y5 is the tidal potential which is unitless and thus needs no conversion. + // y6 is a "potential stress" with units of [m-1] + // y1 (real and imag) + full_solution_ptr[slice_stride + 0] *= displacement_scale; + full_solution_ptr[slice_stride + 1] *= displacement_scale; + + // y3 (real and imag) + full_solution_ptr[slice_stride + 4] *= displacement_scale; + full_solution_ptr[slice_stride + 5] *= displacement_scale; + + // y2 (real and imag) + full_solution_ptr[slice_stride + 2] *= stress_scale; + full_solution_ptr[slice_stride + 3] *= stress_scale; + + // y4 (real and imag) + full_solution_ptr[slice_stride + 6] *= stress_scale; + full_solution_ptr[slice_stride + 7] *= stress_scale; + + // y5 (real and imag) + // full_solution_ptr[slice_stride + 8] *= 1.0; + // full_solution_ptr[slice_stride + 9] *= 1.0; + + // y6 (real and imag) + full_solution_ptr[slice_stride + 10] *= potential_scale; + full_solution_ptr[slice_stride + 11] *= potential_scale; + } + } + } +} \ No newline at end of file diff --git a/TidalPy/RadialSolver/rs_solution_.hpp b/TidalPy/RadialSolver/rs_solution_.hpp index 62769592..5770b4c6 100644 --- a/TidalPy/RadialSolver/rs_solution_.hpp +++ b/TidalPy/RadialSolver/rs_solution_.hpp @@ -7,9 +7,10 @@ #include "eos_solution_.hpp" #include "love_.hpp" +#include "nondimensional_.hpp" -const int MAX_NUM_Y = 6; -const int MAX_NUM_Y_REAL = 12; // Maximum number of y values counting both the real and imaginary portions. +const size_t MAX_NUM_Y = 6; +const size_t MAX_NUM_Y_REAL = 12; // Maximum number of y values counting both the real and imaginary portions. // Error Codes: @@ -67,4 +68,7 @@ class RadialSolutionStorageCC bool array_changed); void set_message(const char* new_message); void find_love(); + void dimensionalize_data( + NonDimensionalScalesCC* nondim_scales, + bool redimensionalize); }; diff --git a/TidalPy/RadialSolver/shooting.pyx b/TidalPy/RadialSolver/shooting.pyx index 2d6cc09d..d6ea316c 100644 --- a/TidalPy/RadialSolver/shooting.pyx +++ b/TidalPy/RadialSolver/shooting.pyx @@ -171,6 +171,14 @@ cdef void cf_shooting_solver( cdef double planet_radius = eos_solution_storage_ptr.radius cdef double surface_gravity = eos_solution_storage_ptr.surface_gravity + printf("SHOOTING-- DIMS::\n") + printf("\t\tr0, r1 = %e %e\n", radius_array_ptr[0], radius_array_ptr[1]) + printf("\t\tg0, g1 = %e %e\n", gravity_array_ptr[0], gravity_array_ptr[1]) + printf("\t\trho0, rho1 = %e %e\n", density_array_ptr[0], density_array_ptr[1]) + printf("\t\tmu0, mu1 = (%e %e) (%e %e)\n", complex_shear_array_ptr[0].real, complex_shear_array_ptr[0].imag, complex_shear_array_ptr[1].real, complex_shear_array_ptr[1].imag) + printf("\t\tK0, K1 = (%e %e) (%e %e)\n", complex_bulk_array_ptr[0].real, complex_bulk_array_ptr[0].imag, complex_bulk_array_ptr[1].real, complex_bulk_array_ptr[1].imag) + + # Find boundary condition at the top of the planet -- this is dependent on the forcing type. # Tides (default here) follow the (y2, y4, y6) = (0, 0, (2l+1)/R) rule # The [5] represents the maximum number of solvers that can be invoked with a single call to radial_solver @@ -180,6 +188,7 @@ cdef void cf_shooting_solver( # 15 = 5 (max_num_solutions) * 3 (number of surface conditions) cdef double[15] boundary_conditions cdef double* bc_pointer = &boundary_conditions[0] + printf("SHOOTING-- BCs:: R= %e; bulk rho= %e\n", planet_radius, planet_bulk_density) cf_get_surface_bc( bc_pointer, # Changed parameter bc_models_ptr, @@ -607,6 +616,7 @@ cdef void cf_shooting_solver( if current_layer_i == start_layer_i: # In the first layer. Use initial condition function to find initial conditions printf("DEBUG- Shooting Method Point \t\t layer = %d; L5a\n", current_layer_i) + printf("SHOOTING STARTING COND:: w = %e; rl = %e; rhol = %e; Kl=(%e, %e); mul=(%e, %e); G=%e", frequency, radius_lower, density_lower, bulk_lower.real, bulk_lower.imag, shear_lower.real, shear_lower.imag, G_to_use) cf_find_starting_conditions( &solution_storage_ptr.success, solution_storage_ptr.message_ptr, @@ -1005,4 +1015,4 @@ cdef void cf_shooting_solver( solution_storage_ptr.success = True solution_storage_ptr.set_message('RadialSolver.ShootingMethod:: completed without any noted issues.\n') - printf("DEBUG- Shooting Method Point - Done!!\n") + printf("DEBUG- Shooting Method Done!!\n") diff --git a/TidalPy/RadialSolver/solver.pyx b/TidalPy/RadialSolver/solver.pyx index 84227ad9..15c908ca 100644 --- a/TidalPy/RadialSolver/solver.pyx +++ b/TidalPy/RadialSolver/solver.pyx @@ -4,7 +4,7 @@ from libc.stdio cimport printf from libc.stdlib cimport exit, EXIT_FAILURE -from libc.math cimport fabs, NAN +from libc.math cimport fabs from libc.string cimport strcpy from libcpp.memory cimport shared_ptr, unique_ptr from libcpp.vector cimport vector @@ -14,12 +14,8 @@ from CyRK cimport PreEvalFunc from TidalPy.logger import get_logger from TidalPy.exceptions import UnknownModelError -from TidalPy.constants cimport d_G, d_MIN_FREQUENCY, d_MAX_FREQUENCY -from TidalPy.utilities.dimensions.nondimensional cimport ( - cf_non_dimensionalize_physicals, - cf_redimensionalize_physicals, - cf_redimensionalize_radial_functions - ) +from TidalPy.constants cimport d_G, d_MIN_FREQUENCY, d_MAX_FREQUENCY, d_NAN_DBL +from TidalPy.utilities.dimensions.nondimensional cimport NonDimensionalScalesCC, cf_build_nondimensional_scales from TidalPy.RadialSolver.rs_solution cimport RadialSolverSolution from TidalPy.RadialSolver.shooting cimport cf_shooting_solver from TidalPy.RadialSolver.matrix cimport cf_matrix_propagate @@ -93,8 +89,8 @@ cdef void cf_radial_solver( num_slices_by_layer_vec.resize(num_layers) cdef size_t layer_slices = 0 - cdef double radius_check = NAN - cdef double layer_upper_radius = NAN + cdef double radius_check = d_NAN_DBL + cdef double layer_upper_radius = d_NAN_DBL # Pull out raw pointers to avoid repeated calls to the getter cdef RadialSolutionStorageCC* solution_storage_ptr = solution_storage_sptr.get() @@ -102,11 +98,12 @@ cdef void cf_radial_solver( # Physical parameters cdef double radius_planet - cdef double G_to_use = NAN - cdef double radius_planet_to_use = NAN - cdef double bulk_density_to_use = NAN - cdef double frequency_to_use = NAN - cdef double surface_pressure_to_use = surface_pressure + cdef double G_to_use = d_NAN_DBL + cdef double radius_planet_to_use = d_NAN_DBL + cdef double bulk_density_to_use = d_NAN_DBL + cdef double frequency_to_use = d_NAN_DBL + cdef double starting_radius_to_use = d_NAN_DBL + cdef double surface_pressure_to_use = d_NAN_DBL # Equation of state variables cdef size_t bottom_slice_index @@ -167,34 +164,51 @@ cdef void cf_radial_solver( # Get other needed inputs radius_planet = radius_array_in_ptr[total_slices - 1] + cdef NonDimensionalScalesCC non_dim_scales if nondimensionalize and solution_storage_ptr.error_code == 0: - cf_non_dimensionalize_physicals( - total_slices, + + # Create scales used to non-dimensionalize various properties + cf_build_nondimensional_scales( + &non_dim_scales, frequency, radius_planet, - planet_bulk_density, - radius_array_in_ptr, - density_array_in_ptr, - NULL, # Pressure ptr (not found yet.) - NULL, # Gravity ptr (not found yet.) - complex_bulk_modulus_in_ptr, - complex_shear_modulus_in_ptr, - &radius_planet_to_use, - &bulk_density_to_use, - &surface_pressure_to_use, - &frequency_to_use, - &G_to_use + planet_bulk_density ) - - printf("NON DIM -- R = %e; rho = %e; P = %e; f = %e; G = %e\n",radius_planet_to_use, bulk_density_to_use, surface_pressure_to_use, frequency_to_use, G_to_use) - for layer_i in range(num_layers): - eos_solution_storage_ptr.upper_radius_bylayer_vec[layer_i] /= radius_planet + printf("nondim scales (ptr = %p):\n", &non_dim_scales) + printf("\t length = %e\n", non_dim_scales.length_conversion) + printf("\t leng3 = %e\n", non_dim_scales.length3_conversion) + printf("\t densit = %e\n", non_dim_scales.density_conversion) + printf("\t pascal = %e\n", non_dim_scales.pascal_conversion) + printf("\t sec = %e\n", non_dim_scales.second_conversion) + printf("\t sec2 = %e\n", non_dim_scales.second2_conversion) + printf("\t mass = %e\n", non_dim_scales.mass_conversion) + + # Convert array pointers + for slice_i in range(total_slices): + radius_array_in_ptr[slice_i] /= non_dim_scales.length_conversion + density_array_in_ptr[slice_i] /= non_dim_scales.density_conversion + complex_bulk_modulus_in_ptr[slice_i] /= non_dim_scales.pascal_conversion + complex_shear_modulus_in_ptr[slice_i] /= non_dim_scales.pascal_conversion - # Change starting radius if it was provided (if it is set to default then it is 0 and this won't have an effect) - starting_radius /= radius_planet + for layer_i in range(num_layers): + eos_solution_storage_ptr.upper_radius_bylayer_vec[layer_i] /= non_dim_scales.length_conversion + + # Convert non-array constants + G_to_use = d_G / (non_dim_scales.length3_conversion / (non_dim_scales.mass_conversion * non_dim_scales.second2_conversion)) + radius_planet_to_use = radius_planet / non_dim_scales.length_conversion + bulk_density_to_use = planet_bulk_density / non_dim_scales.density_conversion + frequency_to_use = frequency / (1. / non_dim_scales.second_conversion) + surface_pressure_to_use = surface_pressure / non_dim_scales.pascal_conversion + starting_radius_to_use = starting_radius / non_dim_scales.length_conversion + + # Scale tolerances + # TODO: What about rtol and atol for EOS and radial solver? + # TODO How about these other tolerances? + # eos_pressure_tol /= non_dim_scales.pascal_conversion + # start_radius_tolerance /= non_dim_scales.length_conversion - # TODO Scale pressure tolerance by pasacal conversion? + printf("NON DIM -- R = %e; rho = %e; P = %e; f = %e; G = %e\n",radius_planet_to_use, bulk_density_to_use, surface_pressure_to_use, frequency_to_use, G_to_use) # Update the radius array inside the C++ classes solution_storage_ptr.change_radius_array(radius_array_in_ptr, total_slices, True) @@ -205,6 +219,7 @@ cdef void cf_radial_solver( bulk_density_to_use = planet_bulk_density frequency_to_use = frequency surface_pressure_to_use = surface_pressure + starting_radius_to_use = starting_radius # Solve the equation of state for the planet @@ -337,43 +352,39 @@ cdef void cf_radial_solver( printf("DEBUG-cf_radial_solver - Post shooting method\n") # Finalize solution storage + cdef size_t solver_i + cdef size_t bc_stride + cdef size_t solver_stride + cdef double complex* full_solution_ptr = NULL if nondimensionalize: - # Redimensionalize eos properties - cf_redimensionalize_physicals( - total_slices, - frequency, - radius_planet, - planet_bulk_density, - eos_solution_storage_ptr.radius_array_vec.data(), - eos_solution_storage_ptr.density_array_vec.data(), - eos_solution_storage_ptr.pressure_array_vec.data(), - eos_solution_storage_ptr.gravity_array_vec.data(), - eos_solution_storage_ptr.complex_bulk_array_vec.data(), - eos_solution_storage_ptr.complex_shear_array_vec.data(), - &radius_planet_to_use, - &bulk_density_to_use, - &frequency_to_use, - &G_to_use - ) - - for layer_i in range(num_layers): - eos_solution_storage_ptr.upper_radius_bylayer_vec[layer_i] *= radius_planet + printf("DEBUG-cf_radial_solver - Redimensionalizing\n") + # Redimensionalize user-provided inputs that are provided as pointers so that this function returns to the same state. + printf("nondim scales (ptr = %p):\n", &non_dim_scales) + printf("\t length = %e\n", non_dim_scales.length_conversion) + printf("\t leng3 = %e\n", non_dim_scales.length3_conversion) + printf("\t densit = %e\n", non_dim_scales.density_conversion) + printf("\t pascal = %e\n", non_dim_scales.pascal_conversion) + printf("\t sec = %e\n", non_dim_scales.second_conversion) + printf("\t sec2 = %e\n", non_dim_scales.second2_conversion) + printf("\t mass = %e\n", non_dim_scales.mass_conversion) + + solution_storage_sptr.get().dimensionalize_data(&non_dim_scales, True) + + # Redimensionalize user-provided inputs that are provided as pointers so that this function returns to the same state. + for slice_i in range(total_slices): + radius_array_in_ptr[slice_i] *= non_dim_scales.length_conversion + density_array_in_ptr[slice_i] *= non_dim_scales.density_conversion + complex_bulk_modulus_in_ptr[slice_i] *= non_dim_scales.pascal_conversion + complex_shear_modulus_in_ptr[slice_i] *= non_dim_scales.pascal_conversion if solution_storage_ptr.success: printf("DEBUG-cf_radial_solver - Successful closeout\n") - if nondimensionalize: - # Redimensionalize the solution - cf_redimensionalize_radial_functions( - solution_storage_ptr.full_solution_vec.data(), - radius_planet, - planet_bulk_density, - total_slices, - num_bc_models) - # Calculate Love numbers printf("DEBUG-cf_radial_solver - Pre find-love\n") solution_storage_ptr.find_love() printf("DEBUG-cf_radial_solver - Post find-love\n") + + printf("radial_solver_cf Finished\n") def radial_solver( @@ -555,17 +566,11 @@ def radial_solver( """ printf("DEBUG-RadialSolver Point 1\n") - cdef size_t total_slices = radius_array.size - - # Unpack inefficient user-provided tuples into bool arrays and pass by pointer - cdef size_t num_layers = len(layer_types) - - # Check on solver type - if use_prop_matrix and num_layers > 1: - raise NotImplementedError("Currently, TidalPy's propagation matrix technique only works for 1-layer worlds. For 2 layer worlds where the lower layer is a liquid: you can start the solver at the bottom of the upper solid layer.") + cdef size_t num_layers = len(layer_types) printf("DEBUG-RadialSolver Point 2\n") + # Perform checks if perform_checks: assert density_array.size == total_slices assert complex_bulk_modulus_array.size == total_slices @@ -585,6 +590,12 @@ def radial_solver( elif fabs(frequency) > d_MAX_FREQUENCY: raise ValueError('Forcing frequency is too large (are you sure you are in rad s-1?).') + if use_prop_matrix: + raise NotImplementedError("Propagation matrix technique is not fully implemented or tested.") + + if use_prop_matrix and num_layers > 1: + raise NotImplementedError("Currently, TidalPy's propagation matrix technique only works for 1-layer worlds. For 2 layer worlds where the lower layer is a liquid: you can start the solver at the bottom of the upper solid layer.") + printf("DEBUG-RadialSolver Point 3\n") # Build array of assumptions # OPT: Perhaps set a maximum number of layers then we can put these on the stack rather than heap allocating them. @@ -777,5 +788,5 @@ def radial_solver( raise_on_fail, ) printf("DEBUG-RadialSolver Point 13b - Post cf_radial_solver call\n") - + printf("RadialSolver Finished\n") return solution diff --git a/TidalPy/utilities/dimensions/nondimensional.pxd b/TidalPy/utilities/dimensions/nondimensional.pxd index d313dc45..ed587d90 100644 --- a/TidalPy/utilities/dimensions/nondimensional.pxd +++ b/TidalPy/utilities/dimensions/nondimensional.pxd @@ -1,43 +1,18 @@ -from TidalPy.utilities.types_x cimport double_numeric -cdef void cf_non_dimensionalize_physicals( - size_t num_radius, - double frequency, - double mean_radius, - double bulk_density, - double* radius_array_ptr, - double* density_array_ptr, - double* pressure_array_ptr, - double* gravity_array_ptr, - double_numeric* bulk_array_ptr, - double_numeric* shear_array_ptr, - double* radius_planet_to_use, - double* bulk_density_to_use, - double* surface_pressure_to_use, - double* frequency_to_use, - double* G_to_use - ) noexcept nogil +cdef extern from "nondimensional_.hpp" nogil: + cdef cppclass NonDimensionalScalesCC: + double second2_conversion + double second_conversion + double length_conversion + double length3_conversion + double density_conversion + double mass_conversion + double pascal_conversion -cdef void cf_redimensionalize_physicals( - size_t num_radius, - double frequency, - double mean_radius, - double bulk_density, - double* radius_array_ptr, - double* density_array_ptr, - double* pressure_array_ptr, - double* gravity_array_ptr, - double_numeric* bulk_array_ptr, - double_numeric* shear_array_ptr, - double* radius_planet_to_use, - double* bulk_density_to_use, - double* frequency_to_use, - double* G_to_use - ) noexcept nogil -cdef void cf_redimensionalize_radial_functions( - double complex* radial_function_ptr, - double mean_radius, - double bulk_density, - size_t num_slices, - size_t num_solutions) noexcept nogil +cdef void cf_build_nondimensional_scales( + NonDimensionalScalesCC* non_dim_scales_ptr, + double frequency, + double mean_radius, + double bulk_density + ) noexcept nogil diff --git a/TidalPy/utilities/dimensions/nondimensional.pyx b/TidalPy/utilities/dimensions/nondimensional.pyx index d8e91803..c736bd0a 100644 --- a/TidalPy/utilities/dimensions/nondimensional.pyx +++ b/TidalPy/utilities/dimensions/nondimensional.pyx @@ -1,4 +1,4 @@ -# distutils: language = c +# distutils: language = c++ # cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False """ Functionality to non-dimensionalize common variables used for multi-layer tidal calculations. @@ -12,280 +12,83 @@ Martens16 : H. Martens, PhD Thesis (CalTech), 2016, DOI: 10.7907/Z9N29TX7 from libc.math cimport sqrt -from TidalPy.constants cimport d_G, d_PI_DBL +from TidalPy.constants cimport d_G, d_PI_DBL, d_NAN_DBL -cdef void cf_non_dimensionalize_physicals( - size_t num_radius, - double frequency, - double mean_radius, - double bulk_density, - double* radius_array_ptr, - double* density_array_ptr, - double* pressure_array_ptr, - double* gravity_array_ptr, - double_numeric* bulk_array_ptr, - double_numeric* shear_array_ptr, - double* radius_planet_to_use, - double* bulk_density_to_use, - double* surface_pressure_to_use, - double* frequency_to_use, - double* G_to_use - ) noexcept nogil: +cdef class NonDimensionalScalesClass: + """ Python wrapper for the `NonDimensionalScales` struct. """ - # Setup loop variables - cdef size_t i + cdef NonDimensionalScalesCC nondim_scales - # Setup conversions - cdef double second2_conversion = 1. / (d_PI_DBL * d_G * bulk_density) - cdef double second_conversion = sqrt(second2_conversion) - cdef double length_conversion = mean_radius - cdef double density_conversion = bulk_density - cdef double mass_conversion = bulk_density * mean_radius**3 - cdef double pascal_conversion = mass_conversion / (length_conversion * second2_conversion) - - # Convert array pointers - for i in range(num_radius): - radius_array_ptr[i] /= length_conversion - density_array_ptr[i] /= density_conversion - bulk_array_ptr[i] /= pascal_conversion - shear_array_ptr[i] /= pascal_conversion - - if gravity_array_ptr: - gravity_array_ptr[i] /= (length_conversion / second2_conversion) - if pressure_array_ptr: - pressure_array_ptr[i] /= pascal_conversion + def __init__(self): + # Initialize everything to nan + self.nondim_scales.second2_conversion = d_NAN_DBL + self.nondim_scales.second_conversion = d_NAN_DBL + self.nondim_scales.length_conversion = d_NAN_DBL + self.nondim_scales.length3_conversion = d_NAN_DBL + self.nondim_scales.density_conversion = d_NAN_DBL + self.nondim_scales.mass_conversion = d_NAN_DBL + self.nondim_scales.pascal_conversion = d_NAN_DBL + + @property + def second2_conversion(self): + return self.nondim_scales.second2_conversion - # Convert non-array pointers - radius_planet_to_use[0] = 1.0 - bulk_density_to_use[0] = 1.0 - G_to_use[0] = d_G / (length_conversion**3 / (mass_conversion * second2_conversion)) - frequency_to_use[0] = frequency / (1. / second_conversion) + @property + def second_conversion(self): + return self.nondim_scales.second_conversion + + @property + def length_conversion(self): + return self.nondim_scales.length_conversion + + @property + def length3_conversion(self): + return self.nondim_scales.length3_conversion + + @property + def density_conversion(self): + return self.nondim_scales.density_conversion - # Scale other inputs - surface_pressure_to_use[0] /= pascal_conversion + @property + def mass_conversion(self): + return self.nondim_scales.mass_conversion + + @property + def pascal_conversion(self): + return self.nondim_scales.pascal_conversion -cdef void cf_redimensionalize_physicals( - size_t num_radius, +cdef void cf_build_nondimensional_scales( + NonDimensionalScalesCC* non_dim_scales_ptr, double frequency, double mean_radius, - double bulk_density, - double* radius_array_ptr, - double* density_array_ptr, - double* pressure_array_ptr, - double* gravity_array_ptr, - double_numeric* bulk_array_ptr, - double_numeric* shear_array_ptr, - double* radius_planet_to_use, - double* bulk_density_to_use, - double* frequency_to_use, - double* G_to_use + double bulk_density ) noexcept nogil: - # Setup loop variables - cdef size_t i - - # Setup conversions - cdef double second2_conversion = 1. / (d_PI_DBL * d_G * bulk_density) - cdef double second_conversion = sqrt(second2_conversion) - cdef double length_conversion = mean_radius - cdef double density_conversion = bulk_density - cdef double mass_conversion = bulk_density * mean_radius**3 - cdef double pascal_conversion = mass_conversion / (length_conversion * second2_conversion) - - # Convert array pointers - for i in range(num_radius): - radius_array_ptr[i] *= length_conversion - density_array_ptr[i] *= density_conversion - bulk_array_ptr[i] *= pascal_conversion - shear_array_ptr[i] *= pascal_conversion - - if gravity_array_ptr: - gravity_array_ptr[i] *= (length_conversion / second2_conversion) - if pressure_array_ptr: - pressure_array_ptr[i] *= pascal_conversion - - # Convert non-array pointers - radius_planet_to_use[0] = mean_radius - bulk_density_to_use[0] = bulk_density - G_to_use[0] = d_G - frequency_to_use[0] = frequency - - -cdef void cf_redimensionalize_radial_functions( - double complex* radial_function_ptr, - double mean_radius, - double bulk_density, - size_t num_slices, - size_t num_solutions) noexcept nogil: - """ A function to re-dimensionalize physical parameters that have been previously non-dimensionalized. - - Parameters - ---------- - radial_function_ptr : complex128* - Non-dimensionalized radial solutions as a function of radius. - mean_radius : float64 - Mean radius of the planet, used in scaling [m] - bulk_density : float64 - Bulk density of the planet, used in scaling [m] - num_slices : uint32 - Number of radial slices, used for looping. - num_solutions : uint32 - Number of solutions to loop through (size of radial_function_ptr is 6 * num_solutions * num_slices) - - """ - # Loop variables - cdef size_t slice_i, solver_i - - # Setup conversions - cdef double second2_conversion = 1. / (d_PI_DBL * d_G * bulk_density) - cdef double mass_conversion = bulk_density * mean_radius**3 - cdef double length_conversion = mean_radius - cdef double length_conversion3 = length_conversion**3 - - for solver_i in range(num_solutions): - for slice_i in range(num_slices): - # Convert displacements - # y1, y3 are the radial and tangential displacements with units of [s2 m-1] - # y2, y4 are the radial and tangential stresses with units of [kg m-3] - # y5 is the tidal potential which is unitless and thus needs no conversion. - # y6 is a "potential stress" with units of [m-1] - - radial_function_ptr[slice_i * 6 * num_solutions + solver_i * 6 + 0] *= \ - (second2_conversion / length_conversion) - radial_function_ptr[slice_i * 6 * num_solutions + solver_i * 6 + 2] *= \ - (second2_conversion / length_conversion) - - radial_function_ptr[slice_i * 6 * num_solutions + solver_i * 6 + 1] *= \ - (mass_conversion / length_conversion3) - radial_function_ptr[slice_i * 6 * num_solutions + solver_i * 6 + 3] *= \ - (mass_conversion / length_conversion3) + non_dim_scales_ptr.second2_conversion = 1. / (d_PI_DBL * d_G * bulk_density) + non_dim_scales_ptr.second_conversion = sqrt(non_dim_scales_ptr.second2_conversion) + non_dim_scales_ptr.length_conversion = mean_radius + non_dim_scales_ptr.length3_conversion = mean_radius * mean_radius * mean_radius + non_dim_scales_ptr.density_conversion = bulk_density + non_dim_scales_ptr.mass_conversion = bulk_density * non_dim_scales_ptr.length3_conversion + non_dim_scales_ptr.pascal_conversion = \ + non_dim_scales_ptr.mass_conversion / (non_dim_scales_ptr.length_conversion * non_dim_scales_ptr.second2_conversion) - radial_function_ptr[slice_i * 6 * num_solutions + solver_i * 6 + 4] *= \ - 1. - radial_function_ptr[slice_i * 6 * num_solutions + solver_i * 6 + 5] *= \ - (1. / length_conversion) - -def non_dimensionalize_physicals( +def build_nondimensional_scales( double frequency, double mean_radius, - double bulk_density, - double surface_pressure, - double[::1] radius_array_view, - double[::1] density_array_view, - double_numeric[::1] bulk_array_view, - double_numeric[::1] shear_array_view, - double[::1] pressure_array_view = None, - double[::1] gravity_array_view = None, + double bulk_density ): - cdef size_t num_radius = radius_array_view.size - cdef double radius_planet_to_use, bulk_density_to_use, frequency_to_use, G_to_use - - cdef double* pressure_array_ptr = NULL - if pressure_array_view is not None: - pressure_array_ptr = &pressure_array_view[0] - - cdef double* gravity_array_ptr = NULL - if gravity_array_view is not None: - gravity_array_ptr = &gravity_array_view[0] - - cdef double surface_pressure_to_use = surface_pressure + cdef NonDimensionalScalesClass non_dim_scales = NonDimensionalScalesClass() - cf_non_dimensionalize_physicals( - num_radius, + cf_build_nondimensional_scales( + &non_dim_scales.nondim_scales, frequency, mean_radius, - bulk_density, - &radius_array_view[0], - &density_array_view[0], - pressure_array_ptr, - gravity_array_ptr, - &bulk_array_view[0], - &shear_array_view[0], - &radius_planet_to_use, - &bulk_density_to_use, - &surface_pressure_to_use, - &frequency_to_use, - &G_to_use + bulk_density ) - - return radius_planet_to_use, bulk_density_to_use, surface_pressure_to_use, frequency_to_use, G_to_use -def redimensionalize_physicals( - double frequency, - double mean_radius, - double bulk_density, - double[::1] radius_array_view, - double[::1] density_array_view, - double_numeric[::1] bulk_array_view, - double_numeric[::1] shear_array_view, - double[::1] pressure_array_view = None, - double[::1] gravity_array_view = None, - ): - - cdef size_t num_radius = radius_array_view.size - cdef double radius_planet_to_use, bulk_density_to_use, frequency_to_use, G_to_use - - cdef double* pressure_array_ptr = NULL - if pressure_array_view is not None: - pressure_array_ptr = &pressure_array_view[0] - - cdef double* gravity_array_ptr = NULL - if gravity_array_view is not None: - gravity_array_ptr = &gravity_array_view[0] - - cf_redimensionalize_physicals( - num_radius, - frequency, - mean_radius, - bulk_density, - &radius_array_view[0], - &density_array_view[0], - pressure_array_ptr, - gravity_array_ptr, - &bulk_array_view[0], - &shear_array_view[0], - &radius_planet_to_use, - &bulk_density_to_use, - &frequency_to_use, - &G_to_use - ) - - return radius_planet_to_use, bulk_density_to_use, frequency_to_use, G_to_use - -def redimensionalize_radial_functions( - double complex[:, ::1] radial_function_view, - double mean_radius, - double bulk_density): - """ A function to re-dimensionalize physical parameters that have been previously non-dimensionalized. - - Parameters - ---------- - radial_function_view : double complex[:, ::1] - Non-dimensionalized radial solutions as a function of radius. - mean_radius : float64 - Mean radius of the planet, used in scaling [m] - bulk_density : float64 - Bulk density of the planet, used in scaling [m] - num_slices : uint32 - Number of radial slices, used for looping. - num_solutions : uint32, default=1 - Number of solutions to loop through (size of radial_function_ptr is 6 * num_solutions * num_slices) - - """ - # Size of arrays - cdef size_t num_solutions, num_slices - num_solutions = (radial_function_view.shape[0] / 6) - num_slices = radial_function_view.shape[1] - - # Call cython function - cf_redimensionalize_radial_functions( - &radial_function_view.T[0, 0], - mean_radius, - bulk_density, - num_slices, - num_solutions - ) + return non_dim_scales diff --git a/TidalPy/utilities/dimensions/nondimensional_.hpp b/TidalPy/utilities/dimensions/nondimensional_.hpp new file mode 100644 index 00000000..fc2b0779 --- /dev/null +++ b/TidalPy/utilities/dimensions/nondimensional_.hpp @@ -0,0 +1,13 @@ +#pragma once + +class NonDimensionalScalesCC +{ +public: + double second2_conversion; + double second_conversion; + double length_conversion; + double length3_conversion; + double density_conversion; + double mass_conversion; + double pascal_conversion; +}; diff --git a/cython_extensions.json b/cython_extensions.json index ebd6a93a..50981b9d 100644 --- a/cython_extensions.json +++ b/cython_extensions.json @@ -32,7 +32,7 @@ "include_dirs": [["TidalPy", "utilities", "dimensions"]], "compile_args": [], "link_args": [], - "is_cpp": false + "is_cpp": true }, "utilities.math.complex": { @@ -81,7 +81,7 @@ "Material.EOS.eos_solution": { "name": "TidalPy.Material.eos.eos_solution", "sources": [["TidalPy", "Material", "eos", "eos_solution.pyx"]], - "include_dirs": [["TidalPy", "Material", "eos"]], + "include_dirs": [["TidalPy", "Material", "eos"], ["TidalPy", "utilities", "dimensions"]], "compile_args": [], "link_args": [], "is_cpp": true @@ -138,7 +138,7 @@ "RadialSolver.solutions": { "name": "TidalPy.RadialSolver.rs_solution", "sources": [["TidalPy", "RadialSolver", "rs_solution.pyx"]], - "include_dirs": [["TidalPy", "RadialSolver"], ["TidalPy", "Material", "eos"]], + "include_dirs": [["TidalPy", "RadialSolver"], ["TidalPy", "Material", "eos"], ["TidalPy", "utilities", "dimensions"]], "compile_args": [], "link_args": [], "is_cpp": true diff --git a/pyproject.toml b/pyproject.toml index 81f42f08..3ccf4b23 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name='TidalPy' -version = '0.6.0a7.dev9' +version = '0.6.0a7.dev12' description='Tidal Dynamics and Thermal-Orbital Evolution Software Suite Implemented in Cython and Python' authors= [ {name = 'Joe P. Renaud', email = 'TidalPy@gmail.com'} @@ -19,7 +19,7 @@ dependencies = [ "psutil>=5.8.0", "pathos>=0.2.0", # Install CyRK requirements - "cyrk==0.11.3", #FIXME "cyrk @ git+https://github.com/jrenaud90/CyRK.git@debug" + "cyrk>=0.11.4, <0.12.0", # Exoplanet data archive "astropy", "astroquery", @@ -116,7 +116,7 @@ requires = [ "scipy>=1.9.3, <1.14", # Issue with scipy==1.14 and MacOS Py_ssize_t vs. long 'cython>=3.0.0', 'wheel>=0.38', - "cyrk==0.11.3", # FIX MEL: "cyrk @ git+https://github.com/jrenaud90/CyRK.git@debug" + "cyrk>=0.11.4, <0.12.0", ] build-backend = "setuptools.build_meta"