From 046182ba124137ad04681a09d61686a2c4dc8a83 Mon Sep 17 00:00:00 2001 From: Nick Murphy Date: Tue, 23 Jul 2024 14:52:10 -0400 Subject: [PATCH] Update astropy.units notebook --- notebooks/astropy-units-completed.ipynb | 248 ++++++++++++++++++++---- 1 file changed, 205 insertions(+), 43 deletions(-) diff --git a/notebooks/astropy-units-completed.ipynb b/notebooks/astropy-units-completed.ipynb index 3dbdf88..2520872 100644 --- a/notebooks/astropy-units-completed.ipynb +++ b/notebooks/astropy-units-completed.ipynb @@ -46,7 +46,9 @@ " !pip install plasmapy==2024.5.0\n", "\n", "import numpy as np\n", - "import astropy.units as u" + "import astropy.units as u\n", + "import matplotlib.pyplot as plt\n", + "from astropy import visualization" ] }, { @@ -92,31 +94,12 @@ "This notebook introduces [astropy.units] with an emphasis on the functionality needed to work with [plasmapy.particles] and [plasmapy.formulary]. We typically import this subpackage as `u`." ] }, - { - "cell_type": "markdown", - "id": "ad28288e", - "metadata": {}, - "source": [ - "## Contents\n", - "\n", - "1. [Unit basics](#Unit-basics)\n", - "2. [Unit operations](#Unit-operations) \n", - "3. [Unit conversions](#Unit-conversions)\n", - "4. [Detaching units and values](#Detaching-units-and-values)\n", - "5. [Equivalencies](#Equivalencies)\n", - "6. [Physical constants](#Physical-constants)\n", - "7. [Units in PlasmaPy](#Units-in-PlasmaPy)\n", - "8. [Optimizing unit operations](#Optimizing-unit-operations)\n", - "9. [Physical Types](#Physical-types)" - ] - }, { "cell_type": "markdown", "id": "a9fad673", "metadata": {}, "source": [ - "## Unit basics\n", - "\n" + "## Unit essentials" ] }, { @@ -124,7 +107,7 @@ "id": "e9ccbb06", "metadata": {}, "source": [ - "We can create a physical quantity by multiplying or dividing a number or array with a unit." + "We can create a _physical quantity_ by multiplying or dividing a number or array with a unit." ] }, { @@ -145,7 +128,7 @@ "source": [ "[`Quantity`]: https://docs.astropy.org/en/stable/api/astropy.units.Quantity.html#astropy.units.Quantity\n", "\n", - "This operation creates a [`Quantity`]: a number, sequence, or array that has been assigned a physical unit." + "This operation creates a [`Quantity`] object: a number, sequence, or array that has been assigned a physical unit." ] }, { @@ -165,7 +148,7 @@ "source": [ "[`Quantity`]: https://docs.astropy.org/en/stable/api/astropy.units.Quantity.html#astropy.units.Quantity\n", "\n", - "We can also create an object by using the [`Quantity`] class itself." + "We can create an object by using the [`Quantity`] class itself." ] }, { @@ -205,7 +188,7 @@ "source": [ "[`Quantity`]: https://docs.astropy.org/en/stable/api/astropy.units.Quantity.html#astropy.units.Quantity\n", "\n", - "We can even create [`Quantity`] objects that are explicitly dimensionless." + "We can even create [`Quantity`] objects that are explicitly dimensionless. " ] }, { @@ -571,19 +554,31 @@ "source": [ "[`astropy.constants`]: https://docs.astropy.org/en/stable/constants/index.html\n", "\n", - "We can use [`astropy.constants`] to access the most commonly needed physical constants." + "fWe can use [`astropy.constants`] to access the most commonly needed physical constants." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "746a79a7", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Name = Speed of light in vacuum\n", + " Value = 299792458.0\n", + " Uncertainty = 0.0\n", + " Unit = m / s\n", + " Reference = CODATA 2018\n", + "The history saving thread hit an unexpected error (OperationalError('attempt to write a readonly database')).History will not be written to the database.\n" + ] + } + ], "source": [ - "from astropy.constants import c, e, k_B\n", - "\n", - "print(c)" + "from astropy import constants\n", + "print(constants.c)" ] }, { @@ -600,13 +595,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "d2d59d08", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "6.962710872930049 MK\n" + ] + } + ], "source": [ "thermal_energy_per_particle = 0.6 * u.keV\n", - "temperature = thermal_energy_per_particle / k_B\n", + "temperature = thermal_energy_per_particle / constants.k_B\n", "print(temperature.to(\"MK\"))" ] }, @@ -615,31 +618,58 @@ "id": "7c145497", "metadata": {}, "source": [ - "Electromagnetic constants often need the unit system to be specified, or will result in a `NameError`." + "Electromagnetic constants often need the unit system to be specified, or will result in an exception." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "015de7fc", "metadata": { "tags": [ "raises-exception" ] }, - "outputs": [], + "outputs": [ + { + "ename": "TypeError", + "evalue": "Constant 'e' does not have physically compatible units across all systems of units and cannot be combined with other values without specifying a system (eg. e.emu)", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[5], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;241;43m2\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mconstants\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43me\u001b[49m\n", + "File \u001b[0;32m/run/media/namurphy/d423d80e-c227-4a33-b50f-545b44160ce3/namurphy/Applications/miniconda3/envs/pldev/lib/python3.12/site-packages/astropy/constants/constant.py:49\u001b[0m, in \u001b[0;36mConstantMeta.__new__..wrap..wrapper\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 47\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msystem \u001b[38;5;129;01mand\u001b[39;00m name_lower \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_has_incompatible_units:\n\u001b[1;32m 48\u001b[0m systems \u001b[38;5;241m=\u001b[39m \u001b[38;5;28msorted\u001b[39m(x \u001b[38;5;28;01mfor\u001b[39;00m x \u001b[38;5;129;01min\u001b[39;00m instances \u001b[38;5;28;01mif\u001b[39;00m x)\n\u001b[0;32m---> 49\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(\n\u001b[1;32m 50\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mConstant \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mabbrev\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;124m does not have physically compatible \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 51\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124munits across all systems of units and cannot be \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 52\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcombined with other values without specifying a \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 53\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124msystem (eg. \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mabbrev\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m.\u001b[39m\u001b[38;5;132;01m{\u001b[39;00msystems[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m)\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 54\u001b[0m )\n\u001b[1;32m 56\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m meth(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n", + "\u001b[0;31mTypeError\u001b[0m: Constant 'e' does not have physically compatible units across all systems of units and cannot be combined with other values without specifying a system (eg. e.emu)" + ] + } + ], "source": [ - "2 * e" + "2 * constants.e" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "1a90c979", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/latex": [ + "$3.2043533 \\times 10^{-19} \\; \\mathrm{C}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "2 * e.si" + "2 * constants.e.si" ] }, { @@ -650,6 +680,108 @@ "Code within PlasmaPy generally uses SI units." ] }, + { + "cell_type": "markdown", + "id": "cd539f42-ce70-4ecf-9bdc-e1f56a86dd07", + "metadata": {}, + "source": [ + "\n", + "## Plotting quantities" + ] + }, + { + "cell_type": "markdown", + "id": "90ced3d3-aa0d-47b7-b526-359e466e5cb5", + "metadata": {}, + "source": [ + "[`Quantity`]: https://docs.astropy.org/en/stable/api/astropy.units.Quantity.html#astropy.units.Quantity\n", + "\n", + "Astropy has built-in support for plotting [`Quantity`] objects. Let's plot the number density of electrons in the solar wind using an empirical formula given by [Kruparova et al. (2023)](https://iopscience.iop.org/article/10.3847/1538-4357/acf572), which has a range of validity from 13 to 50 solar radii." + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "id": "2fce5587-16cf-40ea-baa5-dbd9023281c3", + "metadata": {}, + "outputs": [], + "source": [ + "radii = np.linspace(13, 50, num=50) * constants.R_sun" + ] + }, + { + "cell_type": "markdown", + "id": "d73d6343-d563-4f3d-ae41-a1eefaed4618", + "metadata": {}, + "source": [ + "Next we can apply the formula to get the electron density:\n", + "\n", + "$$ n_e(R) = \\left( 343466\\ \\mbox{cm}^{-3} \\right) × \\left( \\frac{R}{R_☉} \\right)^{-1.87} $$" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "id": "7a734940-a970-4991-9ada-f0844ae80567", + "metadata": {}, + "outputs": [], + "source": [ + "n_e = 343_466 * u.cm**-3 * (radii / constants.R_sun) ** -1.87" + ] + }, + { + "cell_type": "markdown", + "id": "1ed451a9-909f-40ae-bec7-c67d8767545c", + "metadata": {}, + "source": [ + "First let's do some imports:\n", + "\n", + "We can use the [`astropy.visualization.quantity_support`] to help with plotting `Quantity` objects against each other. First let's do some imports." + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "60a5352b-b0ff-46fd-91c5-7acdab398412", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "from astropy.visualization import quantity_support" + ] + }, + { + "cell_type": "markdown", + "id": "07547260-da42-40b7-98de-693b60dbb788", + "metadata": {}, + "source": [ + "Will make use make use of [`astropy.visualization.quantity_support`](https://docs.astropy.org/en/stable/api/astropy.visualization.quantity_support.html). This is a [_context manager_](https://realpython.com/python-with-statement/), which means that we use the `with` statement." + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "a4563556-64d3-4b7e-80f6-7df22ccf4488", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAGyCAYAAADwPVBzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGmUlEQVR4nO3deXxU1f3/8fdkm6yTELJDCAn7riyFqKAIBRQXFK18tYJV8SsFW9RaqvXn0n790tpqXapYaxVbd/sVW3HByKoQFpEIBIgsgQDJJJCQmezb3N8fIVMiYQkJcyeT1/PxuA+SuSeTz/E+cN7cc+45FsMwDAEAAKAZP7MLAAAA8EaEJAAAgBYQkgAAAFpASAIAAGgBIQkAAKAFhCQAAIAWEJIAAABaEGB2AR2Vy+VSfn6+IiIiZLFYzC4HAACcBcMwVFZWpqSkJPn5nf5eESHpHOXn5ys5OdnsMgAAwDk4ePCgunfvfto2hKRzFBERIanxP7LNZjO5GgAAcDacTqeSk5Pdn+OnQ0g6R01DbDabjZAEAEAHczZTZZi4DQAA0AJCEgAAQAsISQAAAC0gJAEAALSAkAQAANACQhIAAEALCEkAAAAtICQBAAC0gJAEAADQAkISAABACwhJAAAALSAkAQAAtICQ5KUMwzC7BAAAOjVCkpf5R+Z+Xf7UKr26dr/ZpQAA0KkRkrxMZW2D9h2pUObeYrNLAQCgUyMkeZkxaV0lSRtzi9XgYsgNAACzEJK8zKAkmyKsAXJW12tngdPscgAA6LQISV4mwN9Po1KjJUnr9zHkBgCAWQhJXmhMGiEJAACzEZK8UHpajCRpQ24J85IAADAJIckLDTw+L6msul478pmXBACAGQhJXsjfz6IfMC8JAABTEZK8VHqvxqUAMglJAACYgpDkpZrWS9qUW6L6BpfJ1QAA0PkQkrzUgESbbMEBKqup1w7WSwIAwOMISV6qcV5S490k5iUBAOB5hCQv1rReEvu4AQDgeYQkL+ael7T/GPOSAADwMEKSFxuYaFNkSKDKa+qVzXpJAAB4FCHJi/mdsF4SSwEAAOBZhCQv1zTkxuRtAAA8i5Dk5dJZLwkAAFMQkrxc/4QIRYYEqqK2QdsOO8wuBwCAToOQ5OX8/Cwa7d7HrcTkagAA6DwISR1A0z5uzEsCAMBzCEkdwH/WSypRHfOSAADwCEJSB9AvPkJdQgNVybwkAAA8hpDUATTOS2LIDQAATyIkdRDs4wYAgGcRkjqIMccnb3+9/xjzkgAA8ABCUgfRNy5C0WFBqqpr0NZDzEsCAOB8IyR1EM3XS2LIDQCA842Q1IGwjxsAAJ5DSOpA0k+Yl1Rbz7wkAADOJ0JSB9InLvyEeUmlZpcDAIBPIyR1IBaLxb0UAENuAACcX4SkDibdPS+JzW4BADifCEkdTNPk7a8PlKimvsHkagAA8F2EpA6md1y4YsKDVF3nYr0kAADOI0JSB2OxWDS6aciNLUoAADhvCEkdUNOQ2zpCEgAA5w0hqQO6pHeMpMZ5SWXVdSZXAwCAbyIkdUCpMWFKiwlTXYOhr3YfNbscAAB8EiGpgxrfP06StGJXkcmVAADgm7wqJC1cuFCjRo1SRESE4uLiNG3aNOXk5DRrc9lll8lisTQ77r777mZt8vLyNHXqVIWGhiouLk4PPPCA6uvrm7VZtWqVhg8fLqvVqt69e2vx4sXnu3vt6vLjIWllzhG5XIbJ1QAA4Hu8KiStXr1ac+fO1fr165WRkaG6ujpNmjRJFRUVzdrNnj1bBQUF7uPJJ590n2toaNDUqVNVW1urdevW6fXXX9fixYv1yCOPuNvk5uZq6tSpGj9+vLKysjR//nzdeeedWrZsmcf62lajekYr3Bqgo+U12naYpQAAAGhvFsMwvPY2xJEjRxQXF6fVq1dr3LhxkhrvJF1wwQV65plnWvyZTz/9VFdddZXy8/MVHx8vSXrppZe0YMECHTlyREFBQVqwYIE+/vhjbd++3f1zM2bMUGlpqT777LOzqs3pdCoyMlIOh0M2m61tHT1Hc97YrE+32/XzCX107w/7mlIDAAAdSWs+v73qTtL3ORyNd0iio6Obvf7mm28qJiZGgwcP1oMPPqjKykr3uczMTA0ZMsQdkCRp8uTJcjqdys7OdreZOHFis/ecPHmyMjMzT1lLTU2NnE5ns8NszEsCAOD8CTC7gFNxuVyaP3++Lr74Yg0ePNj9+s0336yUlBQlJSVp69atWrBggXJycvTBBx9Ikux2e7OAJMn9vd1uP20bp9OpqqoqhYSEnFTPwoUL9fjjj7drH9tqfL/GkLTtsENFzmrF2YJNrggAAN/htSFp7ty52r59u7766qtmr991113ur4cMGaLExERNmDBBe/fuVa9evc5bPQ8++KDuu+8+9/dOp1PJycnn7fedjdgIq4Z1j9S3hxxalXNEPxplbj0AAPgSrxxumzdvnpYuXaqVK1eqe/fup207evRoSdKePXskSQkJCSosLGzWpun7hISE07ax2Wwt3kWSJKvVKpvN1uzwBk1Dbst3FZ6hJQAAaA2vCkmGYWjevHlasmSJVqxYodTU1DP+TFZWliQpMTFRkpSenq5t27apqOg/83QyMjJks9k0cOBAd5vly5c3e5+MjAylp6e3U088Z0L/xmHDr3YfVU19g8nVAADgO7wqJM2dO1dvvPGG3nrrLUVERMhut8tut6uqqkqStHfvXv32t7/V5s2btX//fv373//WzJkzNW7cOA0dOlSSNGnSJA0cOFC33nqrvv32Wy1btkwPP/yw5s6dK6vVKkm6++67tW/fPv3yl7/Url279OKLL+q9997Tvffea1rfz9WgJJtiI6yqqG3QptxjZpcDAIDP8KqQtGjRIjkcDl122WVKTEx0H++++64kKSgoSF988YUmTZqk/v376/7779f06dP10Ucfud/D399fS5culb+/v9LT0/XjH/9YM2fO1G9+8xt3m9TUVH388cfKyMjQsGHD9NRTT+mVV17R5MmTPd7ntvLzs2h8v1hJDLkBANCevHqdJG/mDeskNflsu113v7FZPbuGatUD402tBQAAb+Yz6yTh7FzSJ0aB/hbtL67UviPlZpcDAIBPICT5gHBrgEandpXEwpIAALQXQpKPuJzVtwEAaFeEJB/RFJI25paorLrO5GoAAOj4CEk+omdMmNJiwlTvMvTl7qNmlwMAQIdHSPIhDLkBANB+CEk+pCkkrcopksvFyg4AALQFIcmHjOwZrXBrgI6W12rrYYfZ5QAA0KERknxIUICfxvWNkcSQGwAAbUVI8jHj+zUOua0kJAEA0CaEJB9z2fGQtO2wQ0XOapOrAQCg4yIk+ZjYCKuGJUdJklbmcDcJAIBzRUjyQZf3YykAAADaipDkg5qWAvhy91HV1DeYXA0AAB0TIckHDUqyKS7CqsraBm3MLTG7HAAAOiRCkg/y87O4n3JjyA0AgHNDSPJR40/YosQwWH0bAIDWIiT5qLF9YmQN8NOB4krtKHCaXQ4AAB0OIclHhVkD3ENuS7cWmFwNAAAdDyHJh109LEmS9NG3+Qy5AQDQSoQkH3Z5/ziFBvnr0LEqfXuIDW8BAGgNQpIPCwny18QB8ZIa7yYBAICzR0jycVcNTZQkfby1QC4XQ24AAJwtQpKPu7RfrCKsAbI7q7U575jZ5QAA0GEQknycNcBfkwYlSGLIDQCA1iAkdQJXDWsccvtkW4HqG1wmVwMAQMdASOoELukdo6jQQB0tr9UG9nIDAOCsEJI6gUB/P10xuHHIbelWhtwAADgbhKRO4uqhjQtLfrrdrjqG3AAAOCNCUicxOq2rYsKtKq2s01d7jppdDgAAXo+Q1En4+1k0dQhPuQEAcLYISZ3IVcf3csvILlR1XYPJ1QAA4N0ISZ3IiB5dlBgZrLKaeq3+7ojZ5QAA4NUISZ2In59FU4c0rpm0dGuBydUAAODdCEmdzNXHh9y+2FGoytp6k6sBAMB7EZI6maHdI9UjOlRVdQ1asavI7HIAAPBahKROxmKxaOrQxiE3nnIDAODUCEmdUNPCkitzjqisus7kagAA8E6EpE5oQGKE0mLDVFvvUsaOQrPLAQDAKxGSOiGLxeK+m8RTbgAAtIyQ1EldPaxxXtKa746otLLW5GoAAPA+hKROqndchPonRKjeZWhZtt3scgAA8DqEpE6sac2kj75lyA0AgO8jJHViVx1fCmDd3qM6UlZjcjUAAHgXQlInltI1TBckR8llSB98c8jscgAA8CqEpE7uplHJkqR3Nx2UYRgmVwMAgPcgJHVyVw9LUmiQv/YdrdDG3BKzywEAwGsQkjq5cGuArjk+gfudTQdNrgYAAO9BSIJm/KCHJOmTbQVyVLJNCQAAEiEJkoZ1j1T/hAjV1Lu0ZAsTuAEAkAhJUOM2Jf91/G7SO0zgBgBAEiEJx027oJusAX7aZS/Tt4ccZpcDAIDpCEmQJEWGBurKIY2LS76zMc/kagAAMJ9XhaSFCxdq1KhRioiIUFxcnKZNm6acnJxmbaqrqzV37lx17dpV4eHhmj59ugoLC5u1ycvL09SpUxUaGqq4uDg98MADqq+vb9Zm1apVGj58uKxWq3r37q3Fixef7+55vRnH10z697f5Kq+pP0NrAAB8m1eFpNWrV2vu3Llav369MjIyVFdXp0mTJqmiosLd5t5779VHH32k999/X6tXr1Z+fr6uv/569/mGhgZNnTpVtbW1WrdunV5//XUtXrxYjzzyiLtNbm6upk6dqvHjxysrK0vz58/XnXfeqWXLlnm0v97mB6nRSosNU2Vtgz76Nt/scgAAMJXF8OJZukeOHFFcXJxWr16tcePGyeFwKDY2Vm+99ZZuuOEGSdKuXbs0YMAAZWZmasyYMfr000911VVXKT8/X/Hx8ZKkl156SQsWLNCRI0cUFBSkBQsW6OOPP9b27dvdv2vGjBkqLS3VZ599dla1OZ1ORUZGyuFwyGaztX/nTfLymr363092aVj3SP1r3iVmlwMAQLtqzee3V91J+j6Ho3ECcXR0tCRp8+bNqqur08SJE91t+vfvrx49eigzM1OSlJmZqSFDhrgDkiRNnjxZTqdT2dnZ7jYnvkdTm6b3aElNTY2cTmezwxdNH95dgf4WfXvIoR35vtlHAADOhteGJJfLpfnz5+viiy/W4MGDJUl2u11BQUGKiopq1jY+Pl52u93d5sSA1HS+6dzp2jidTlVVVbVYz8KFCxUZGek+kpOT29xHb9Q13KpJAxMkSe9uYgI3AKDz8tqQNHfuXG3fvl3vvPOO2aVIkh588EE5HA73cfCg727h0bTp7ZIth1Vd12ByNQAAmMMrQ9K8efO0dOlSrVy5Ut27d3e/npCQoNraWpWWljZrX1hYqISEBHeb7z/t1vT9mdrYbDaFhIS0WJPVapXNZmt2+KpLeseoe5cQOavr9cm2ArPLAQDAFF4VkgzD0Lx587RkyRKtWLFCqampzc6PGDFCgYGBWr58ufu1nJwc5eXlKT09XZKUnp6ubdu2qaioyN0mIyNDNptNAwcOdLc58T2a2jS9R2fn52fRTSMb7ya9s9F375gBAHA6XhWS5s6dqzfeeENvvfWWIiIiZLfbZbfb3fOEIiMjdccdd+i+++7TypUrtXnzZv3kJz9Renq6xowZI0maNGmSBg4cqFtvvVXffvutli1bpocfflhz586V1WqVJN19993at2+ffvnLX2rXrl168cUX9d577+nee+81re/e5saRyfKzSBv3l2hPUbnZ5QAA4HFeFZIWLVokh8Ohyy67TImJie7j3Xffdbf505/+pKuuukrTp0/XuHHjlJCQoA8++MB93t/fX0uXLpW/v7/S09P14x//WDNnztRvfvMbd5vU1FR9/PHHysjI0LBhw/TUU0/plVde0eTJkz3aX2+WEBmsy/vHSWICNwCgc/LqdZK8ma+uk3SiL3YU6s6/f63osCCtf3CCggK8KlMDANBqPrNOEsx1Wb9YxdusKqmoVcaOwjP/AAAAPoSQhFMK8PfTj5omcDPkBgDoZAhJOK2mkPTl7qM6WFJpcjUAAHgOIQmnlRwdqrF9YiRJ/1h/wORqAADwHEISzugnF/eUJL29IU9l1XXmFgMAgIcQknBGl/WNU++4cJXV1OvdTSwuCQDoHAhJOCM/P4vuvKRx9fPX1u5XXYPL5IoAADj/CEk4K9Mu7KaY8CAdLq1iPzcAQKdASMJZCQ7018z0npKkv365T6xBCgDwdYQknLUfj0lRcKCfth92av2+ErPLAQDgvCIk4axFhwXphhHdJUmvfLnP5GoAADi/CElolTsuSZPFIi3fVaQ9RWVmlwMAwHlDSEKrpMaE6YcD4iVJf/sq1+RqAAA4fwhJaLXZ49IkSf/3zWEdKasxuRoAAM4PQhJabWRKF12QHKXaehdblQAAfBYhCa1msVg0e2zj3aR/ZO5XVW2DyRUBAND+CEk4J5MHxSs5OkTHKuv0f98cMrscAADaHSEJ5yTA30+3X9y4VcnfvsqVy8XikgAA30JIwjn70chk2YIDlHu0Ql/sLDS7HAAA2hUhCecszBqgW8akSGrcqgQAAF9CSEKb3HZRTwX6W7Rp/zFtyTtmdjkAALQbQhLaJN4WrGuGdZMkvfIli0sCAHwHIQltdufYxgncn24v0MGSSpOrAQCgfRCS0GYDEm0a2ydGLkN6eQ1zkwAAvoGQhHYx57JekqR3Nx3U4dIqk6sBAKDtCEloFxf1itGYtGjVNrj05xW7zS4HAIA2IySh3dw/qZ8k6f2vD+lAcYXJ1QAA0DaEJLSbUT2jNa5vrOpdhp5dzt0kAEDHRkhCu7r/h30lSR9uOaw9ReUmVwMAwLkjJKFdDUuO0sQB8XIZ4m4SAKBDIySh3d13/G7SR9/ma5fdaXI1AACcG0IS2t3AJJumDkmUJP0p4zuTqwEA4NwQknBezJ/YRxaLtCy7UNsOOcwuBwCAViMk4bzoEx+haRc07un2dEaOydUAANB6hCScNz+f0Ef+fhatzDmizQeOmV0OAACtQkjCedMzJkzTh3M3CQDQMRGScF7dc3kfBfpbtHZPsdbvKza7HAAAzlqrQ9KxY8dUUlIiSTpy5Ig++OADZWdnt3th8A3J0aG6aVSyJOnpz7+TYRgmVwQAwNlpVUh65ZVXNGLECI0cOVKLFi3Sddddp+XLl2vGjBl65ZVXzleN6ODmje+joAA/bdxfoi93HzW7HAAAzkpAaxo/99xzys7OVlVVlXr06KHc3FzFxsbK4XDo0ksv1Z133nm+6kQHlhAZrB+PTtGra3P1VMZ3GtsnRhaLxeyyAAA4rVbdSQoICFBISIiio6PVu3dvxcbGSpIiIyP50MNpzbmsl0IC/fXtwVKt2FVkdjkAAJxRq0KSv7+/qqurJUmrV692v15ezkamOL3YCKtmXdRTkvSHZTmqb3CZWxAAAGfQqpD0xRdfyGq1Smq8e9SksrJSL7/8cvtWBp9z96VpigwJ1C57md7edNDscgAAOK1WhaTvD6vZ7XZJUlxcnEaNGtW+lcHnRIUGuTe/ferzHJVW1ppcEQAAp9amdZImTZrUXnWgk7hldA/1i49QaWWdnvlit9nlAABwSm0KSax5g9YK8PfTI1cPlCT9Y/0B5djLTK4IAICWtSkk8UQbzsXFvWM0eVC8GlyGfrM0m7ANAPBKbEsCUzw8daCCAvy0dk+xlmUXml0OAAAnISTBFMnRobprbJok6YlPdqi6rsHkigAAaK5NIcnf37+96kAn9NPxvZRgC9bBkir97atcs8sBAKCZNoWkLVu2tFcd6IRCgwL0qyv6S5JeWLlHdke1yRUBAPAfDLfBVNdekKQRKV1UWdug33+2y+xyAABwa9UGt6dSXV2trVu3qqioSC5X8+0mrrnmmvb4FfBRFotFj109SNe88JWWbDmsH49J0YiULmaXBQBA2+8kffbZZ+rRo4fGjBmja665RtOmTXMf1113Xavfb82aNbr66quVlJQki8WiDz/8sNn52267TRaLpdkxZcqUZm1KSkp0yy23yGazKSoqSnfcccdJ+8tt3bpVY8eOVXBwsJKTk/Xkk0+2ula0jyHdI3XjiO6SpMc/ypbLxZIAAADztTkk3XPPPbrxxhtVUFAgl8vV7GhoaP0TSxUVFRo2bJheeOGFU7aZMmWKCgoK3Mfbb7/d7Pwtt9yi7OxsZWRkaOnSpVqzZo3uuusu93mn06lJkyYpJSVFmzdv1h/+8Ac99thj7D9nogcm91eENUBbDzn0z28OmV0OAABtH24rLCzUfffdp/j4+PaoR1dccYWuuOKK07axWq1KSEho8dzOnTv12WefadOmTRo5cqQk6fnnn9eVV16pP/7xj0pKStKbb76p2tpavfrqqwoKCtKgQYOUlZWlp59+ulmYgufERlj1swl99MQnO/XkZzm6YnCCIoIDzS4LANCJtflO0g033KBVq1a1Qylnb9WqVYqLi1O/fv00Z84cFRcXu89lZmYqKirKHZAkaeLEifLz89OGDRvcbcaNG6egoCB3m8mTJysnJ0fHjh1r8XfW1NTI6XQ2O9C+Zl3UU2kxYTpaXqPnV+wxuxwAQCfX5jtJf/7zn3XjjTfqyy+/1JAhQxQY2Pxf/z/72c/a+iuamTJliq6//nqlpqZq7969euihh3TFFVcoMzNT/v7+stvtiouLa/YzAQEBio6Olt1ulyTZ7XalpqY2a9N0J8xut6tLl5MnDi9cuFCPP/54u/YFzQUF+On/XTVQP1m8Sa9+latpF3TTwCSb2WUBADqpNoekt99+W59//rmCg4O1atWqZvu5WSyWdg9JM2bMcH89ZMgQDR06VL169dKqVas0YcKEdv1dJ3rwwQd13333ub93Op1KTk4+b7+vsxrfP06TB8VrWXahFvzfVi356UUK8GelCgCA57X50+fXv/61Hn/8cTkcDu3fv1+5ubnuY9++fe1R42mlpaUpJiZGe/Y0Ds8kJCSoqKioWZv6+nqVlJS45zElJCSosLD5fmFN359qrpPVapXNZmt24Pz47bWDZQsO0LbDDv31S1biBgCYo80hqba2VjfddJP8/Mz51/6hQ4dUXFysxMRESVJ6erpKS0u1efNmd5sVK1bI5XJp9OjR7jZr1qxRXV2du01GRob69evX4lAbPCvOFqxHrh4kSfrTF99p75HyM/wEAADtr83JZtasWXr33XfboxZJUnl5ubKyspSVlSVJys3NVVZWlvLy8lReXq4HHnhA69ev1/79+7V8+XJde+216t27tyZPnixJGjBggKZMmaLZs2dr48aNWrt2rebNm6cZM2YoKSlJknTzzTcrKChId9xxh7Kzs/Xuu+/q2WefbTacBnNNH95Nl/aNVW29Swv+uZW1kwAAHmcxDKNNnz4/+9nP9Pe//13Dhg3T0KFDT5q4/fTTT7fq/VatWqXx48ef9PqsWbO0aNEiTZs2TVu2bFFpaamSkpI0adIk/fa3v222BEFJSYnmzZunjz76SH5+fpo+fbqee+45hYeHu9ts3bpVc+fO1aZNmxQTE6N77rlHCxYsOOs6nU6nIiMj5XA4GHo7Tw6XVmnS06tVUdugx64eqNsuTj3zDwEAcBqt+fxuc0hqKdCcaOXKlW15e69FSPKMf6w/oP/34XaFBPrr83vHKTk61OySAAAdmEdDUmdFSPIMl8vQjL+u18bcEl3cu6veuGN0sycoAQBojdZ8frd5TtLChQv16quvnvT6q6++qt///vdtfXt0cn5+Fj05faiCA/20dk+x3t100OySAACdRJtD0l/+8hf179//pNcHDRqkl156qa1vD6hnTJju/2E/SdITH++U3VFtckUAgM6gzSHJbre7H78/UWxsrAoKCtr69oAk6fZLUjUsOUplNfX69ZJtYpQYAHC+tTkkJScna+3atSe9vnbtWvcj90Bb+ftZ9IcbhirQ36Llu4r072/zzS4JAODj2hySZs+erfnz5+u1117TgQMHdODAAb366qu69957NXv27PaoEZAk9Y2P0D2X95EkPfbvbB0trzG5IgCAL2vz3m0PPPCAiouL9dOf/lS1tbWSpODgYC1YsEAPPvhgmwsETjTnsl76dLtdOwucevTf2Xrh5uFmlwQA8FHttgRAeXm5du7cqZCQEPXp00dWq7U93tZrsQSAebYfdujaF9aqwWXo2RkX6NoLupldEgCgg/DoEgBNwsPDNWrUKA0ePNjnAxLMNbhbpOaO7y1J+vWS7TpQXGFyRQAAX2TOrrRAG/3s8t76Qc9oldfU6563t6i23mV2SQAAH0NIQocU4O+nZ2ZcoMiQQG095NAflu0yuyQAgI8hJKHDSooK0R9uGCpJ+uuXuVqZU2RyRQAAX0JIQoc2aVCCZqWnSJLuf+9bFTpZjRsA0D4ISejwHrxygAYk2lRSUat7381Sg4vVuAEAbUdIQocXHOivP998oUKD/LVub7FeWr3X7JIAAD6AkASf0Cs2XI9fM0iS9HTGd9p8oMTkigAAHR0hCT7jhhHdNe2CJDW4DP3s7Sw5KuvMLgkA0IERkuAzLBaL/ue6IerZNVSHS6u04P+2qp0WlAcAdEKEJPiUcGuAnv+v4Qr0t+izbLve2JBndkkAgA6KkASfM6R7pBZM6S9J+u3SHcrOd5hcEQCgIyIkwSfdcUmqLu8fp9p6l+76+2YVl9eYXRIAoIMhJMEnWSwW/elHF7jnJ8158xv2dwMAtAohCT4rMjRQr8waqXBrgDbmlujxj7LNLgkA0IEQkuDTesdF6Ln/ukAWi/Tmhjz9Y/0Bs0sCAHQQhCT4vMv7x+uByf0kSY//O1vr9xWbXBEAoCMgJKFTmHNpL10zLEn1LkM/ffMbHSypNLskAICXIyShU7BYLHryhqEa0i1SJRW1mv33r1VRU292WQAAL0ZIQqcRHOivl2eOUEy4VbvsZbr/vW/lcrEiNwCgZYQkdCqJkSH6y60jFOTvp8+y7XpuxW6zSwIAeClCEjqdESld9D/XDZYkPfPFbn22vcDkigAA3oiQhE7pRyOTdfvFqZKke9/9VjvynSZXBADwNoQkdFoPXdlfY/vEqKquQT9ZvFGHjvHEGwDgPwhJ6LQC/P305/8arn7xESp01mjm3zayxxsAwI2QhE4tMjRQr9/+A3WLCtG+oxW6ffEmlgYAAEgiJAFKiAzW67f/QF1CA/XtIYfufmMzm+ECAAhJgCT1jgvXq7eNUkigv77cfVS/eJ81lACgsyMkAcdd2KOLXrp1hAL8LPr3t/n67cc7ZBgEJQDorAhJwAku7RurP944TJL02tr9WrR6r8kVAQDMQkgCvmfahd308NQBkqQnP8vRe5sOmlwRAMAMhCSgBXeOTdPdl/aSJP3qg63K2FFockUAAE8jJAGnsGBKP90wortchjTvrW+0MbfE7JIAAB5ESAJOwWKx6HfXD9GE/nGqqXfp9sWbtPkAQQkAOgtCEnAaAf5++vPNwzUmLVrlNfWa+beN+no/QQkAOgNCEnAGIUH+eu22Hyg9rasqahs069WN2kRQAgCfR0gCzkJIkL9evW2ULur1n6DEHCUA8G2EJOAshQT562+zRumS3jGqrG3Qba9t1IZ9xWaXBQA4TwhJQCuEBPnrlVkjNbZPY1D6yeJNWk9QAgCfREgCWik40F9/nXlCUHptkzL3EpQAwNcQkoBz0BSUxvWNVVVdg25fvEnr9h41uywAQDsiJAHnKDjQXy/fOkKXnhiU9hCUAMBXEJKANggO9Ndfbh2hy/rFqrrOpZ8s3qRl2XazywIAtANCEtBGTUFp4oB41dS7NOeNzXprQ57ZZQEA2oiQBLQDa4C/XvrxcN00MlkuQ3poyTY988V3MgzD7NIAAOeIkAS0kwB/P/1u+hDdc3lvSdIzX+zWQ0u2q8FFUAKAjsjrQtKaNWt09dVXKykpSRaLRR9++GGz84Zh6JFHHlFiYqJCQkI0ceJE7d69u1mbkpIS3XLLLbLZbIqKitIdd9yh8vLyZm22bt2qsWPHKjg4WMnJyXryySfPd9fQCVgsFt0/qZ9+e+0gWSzS2xvzNOeNzaquazC7NABAK3ldSKqoqNCwYcP0wgsvtHj+ySef1HPPPaeXXnpJGzZsUFhYmCZPnqzq6mp3m1tuuUXZ2dnKyMjQ0qVLtWbNGt11113u806nU5MmTVJKSoo2b96sP/zhD3rsscf08ssvn/f+oXO4Nb2nFt0yXEEBfvp8R6Fu/dsGOSrrzC4LANAKFsOLJ01YLBYtWbJE06ZNk9R4FykpKUn333+/fvGLX0iSHA6H4uPjtXjxYs2YMUM7d+7UwIEDtWnTJo0cOVKS9Nlnn+nKK6/UoUOHlJSUpEWLFunXv/617Ha7goKCJEm/+tWv9OGHH2rXrl1nVZvT6VRkZKQcDodsNlv7dx4+Yf2+Ys3++9cqq65X3/hwvX77D5QYGWJ2WQDQabXm89vr7iSdTm5urux2uyZOnOh+LTIyUqNHj1ZmZqYkKTMzU1FRUe6AJEkTJ06Un5+fNmzY4G4zbtw4d0CSpMmTJysnJ0fHjh1r8XfX1NTI6XQ2O4AzGZPWVe/fna54m1XfFZZr+ovrtLuwzOyyAABnoUOFJLu9cf2Z+Pj4Zq/Hx8e7z9ntdsXFxTU7HxAQoOjo6GZtWnqPE3/H9y1cuFCRkZHuIzk5ue0dQqfQP8Gm/5tzkXrFhinfUa0bXsrUWhadBACv16FCkpkefPBBORwO93Hw4EGzS0IH0r1LqP5590Ua3iNKjqo6zXx1oxavzWWJAADwYh0qJCUkJEiSCgsLm71eWFjoPpeQkKCioqJm5+vr61VSUtKsTUvvceLv+D6r1SqbzdbsAFqjS1iQ3po9Rtdf2E0NLkOPfbRDD36wTbX1LrNLAwC0oEOFpNTUVCUkJGj58uXu15xOpzZs2KD09HRJUnp6ukpLS7V582Z3mxUrVsjlcmn06NHuNmvWrFFd3X+eNsrIyFC/fv3UpUsXD/UGnVFwoL+e+tEw/frKAfKzSO9sOqhbXlmvo+U1ZpcGAPgerwtJ5eXlysrKUlZWlqTGydpZWVnKy8uTxWLR/Pnz9T//8z/697//rW3btmnmzJlKSkpyPwE3YMAATZkyRbNnz9bGjRu1du1azZs3TzNmzFBSUpIk6eabb1ZQUJDuuOMOZWdn691339Wzzz6r++67z6ReozOxWCyaPS5Nf7ttlCKCA7Rp/zFd8/xXys53mF0aAOAEXrcEwKpVqzR+/PiTXp81a5YWL14swzD06KOP6uWXX1ZpaakuueQSvfjii+rbt6+7bUlJiebNm6ePPvpIfn5+mj59up577jmFh4e722zdulVz587Vpk2bFBMTo3vuuUcLFiw46zpZAgDtYe+Rcs1+/WvtO1qh4EA/PXXjBZo6NNHssgDAZ7Xm89vrQlJHQUhCe3FU1emet7dozXdHJEn3XN5b907sKz8/i8mVAYDv8dl1kgBfFBkSqNduG6XZY1MlSc+v2KP/fmOzyqpZoRsAzERIAryAv59Fv546UE/dOExB/n7K2FGoq5//StsPM08JAMxCSAK8yPQR3fXuf49RUmSw9hdX6voX1+kf6w+wnhIAmICQBHiZC3t00Sc/H6uJA+JU2+DS//twu+a9tUVOht8AwKMISYAXigoN0l9njtTDUwcowM+ij7cV6Ornv9K2Qwy/AYCnEJIAL2WxWHTn2DS9f3e6ukWF6EBxpaYvWqfX1+1n+A0APICQBHi5C3t00Sc/G6sfDoxXbYNLj/47Wz998xs5qhh+A4DziZAEdACRoYF6+dYReuSqgQr0t+jT7XZd9fyX2pJ3zOzSAMBnEZKADsJisej2S1L1z7svUvcuITpYUqXpi9bpj8ty2CQXAM4DQhLQwQxLjtLHPxuray9IksuQ/rxyj6a9sFa77E6zSwMAn0JIAjqgyJBAPTvjQr14y3B1CQ3UjgKnrnl+rRat2qsGF5O6AaA9EJKADuzKIYladu8495pKv/9sl370l0ztP1phdmkA0OERkoAOLi4iWH+dOVJ/uGGowq0B2nzgmK549kv9I5OlAgCgLQhJgA+wWCy6cWSyPps/VulpXVVV16D/969szXx1o/JLq8wuDwA6JEIS4EO6dwnVm3eO1mNXD1RwoJ++3H1UP3x6tV5bm8tcJQBoJUIS4GP8/Cy67eJUffKzsRqR0kUVtQ16/KMduu7Ftdp+mG1NAOBsEZIAH5UWG673/ztdT1w3WBHBAdp6yKFrX1irJz7eocraerPLAwCvR0gCfJifn0W3jE7R8vsu1dShiWpwGfrrl7n64dNrtGJXodnlAYBXIyQBnUCcLVgv3Dxcr902St2iQnS4tEq3L/5ac9/8RkXOarPLAwCvREgCOpHx/eOUcd84/fe4NPn7WfTxtgJNeGq1/p65X/UNbG0CACeyGCykck6cTqciIyPlcDhks9nMLgdotex8hx76YJu+PdQ4mbt/QoQeuWqgLuodY3JlAHD+tObzm5B0jghJ8AUNLkNvbjigpz7/To6qOknS5EHx+vWVA9Wja6jJ1QFA+yMkeQAhCb7kWEWtnvniO72xIU8NLkNB/n66Y2yq5o7vrXBrgNnlAUC7ISR5ACEJvui7wjL9dukOfbn7qCQpNsKqX07up+nDu8vPz2JydQDQdoQkDyAkwVcZhqEvdhbpiY93aH9xpSRpaPdIPXr1QI1IiTa5OgBoG0KSBxCS4Otq6hu0eO1+Pb9ij8prGhefnDwoXr+Y1E994iNMrg4Azg0hyQMISegsjpTV6I/LcvT+5oNyGZKfRZo+vLvm/7CvukWFmF0eALQKIckDCEnobHYXlumPn+doWXbjSt1B/n66NT1FP72sl7qGW02uDgDODiHJAwhJ6Ky25B3T7z/bpfX7SiRJ4dYAzR6bpjvGpvIkHACvR0jyAEISOjPDMPTl7qN6ctkubT/slCR1DQvS3PG9dfPoHgoO9De5QgBoGSHJAwhJgORyGfp0u11//DxHuUcrJElxEVb996W9dPMPeigkiLAEwLsQkjyAkAT8R12DS+9/fUh/XrFb+Y7GDXNjwoN017g03TI6RWEMwwHwEoQkDyAkASerrXfp/745pBdX7dHBkipJUpfQQN05Nk0z01MUERxocoUAOjtCkgcQkoBTq2tw6cMth/XCyj3uBSkjQwJ1+8Wpuu3inooMISwBMAchyQMIScCZ1Te4tHRrgZ5fsVt7jzTOWYqwBujH6Sm67aKeircFm1whgM6GkOQBhCTg7DW4DH2yrUB/XrFHOYVlkqRAf4umXdBNs8elqS8reAPwEEKSBxCSgNZzuQx9sbNQf/1ynzbtP+Z+fXy/WM0el6b0tK6yWNhIF8D5Q0jyAEIS0Dbf5B3TX9fs02fZdjX9X2hIt0jNHpemKwcnKMDfz9wCAfgkQpIHEJKA9rH/aIVeXZur974+qOo6lySpW1SIfnJxT904IlmRoUzyBtB+CEkeQEgC2ldJRa3eWH9Ar6/br+KKWklSSKC/pl3YTTPTUzQgkb9nANqOkOQBhCTg/Kiua9CSLYf1+rr92mUvc7/+g57RmnlRiiYPSlAgQ3EAzhEhyQMIScD5ZRiGNu0/ptcz92vZdrvqXY3/q4qLsOqW0Sn6r9HJiotgCQEArUNI8gBCEuA5dke13tqYp7c25OloeY2kxiUEJg9K0H/9oIfS07rKz4+n4gCcGSHJAwhJgOfV1rv06fYC/SPzgL4+8J8lBJKjQ3TTyGTdMCJZCZHcXQJwaoQkDyAkAebaftihdzcd1IdbDquspl6S5GeRLu8fp5tG9dD4frEsIwDgJIQkDyAkAd6hqrZBn2wr0LubDmrj/hL363ERVt04srt+NDJZKV3DTKwQgDchJHkAIQnwPnuKyvXe1wf1f5sPuZcRkKRRPbvougu7a+qQRNZdAjo5QpIHEJIA71Vb79LynYV6e9NBfbn7iHtF76AAP00cEKfrL+yuS/vFspQA0AkRkjyAkAR0DHZHtf6VdVgffHPYvbmuJEWHBemaYUm67sJuGto9kj3jgE6CkOQBhCSgYzEMQzsKnFryzWF9mJXvXkpAknrFhunqYUm6amiSeseFm1glgPONkOQBhCSg46pvcOmrPUf1wTeH9fkOu3vPOEnqnxBxPDAlMuEb8EGEJA8gJAG+oay6Tp9nF2rp1nx9ufuoe2VvSRrSLVJXDU3U1KGJ6t4l1MQqAbQXQpIHEJIA31NaWatl2XYt3VqgdXuL1XBCYLqwR5SuHJyoyYMS1KMrgQnoqFrz+d3hHu147LHHZLFYmh39+/d3n6+urtbcuXPVtWtXhYeHa/r06SosLGz2Hnl5eZo6dapCQ0MVFxenBx54QPX19Z7uCgAvExUapJtG9dA/7hitjQ9N0P9MG6wxadGyWKQteaV64pOdGveHlbri2S/1zBffaZfdKf6dCfiuALMLOBeDBg3SF1984f4+IOA/3bj33nv18ccf6/3331dkZKTmzZun66+/XmvXrpUkNTQ0aOrUqUpISNC6detUUFCgmTNnKjAwUP/7v//r8b4A8E5dw6368ZgU/XhMioqc1fp0u13Lsu3akFuinQVO7Sxw6pkvdiula6imDErQpEEJujA5ij3kAB/S4YbbHnvsMX344YfKyso66ZzD4VBsbKzeeust3XDDDZKkXbt2acCAAcrMzNSYMWP06aef6qqrrlJ+fr7i4+MlSS+99JIWLFigI0eOKCgo6KzqYLgN6JyOVdTqi52FWpZt15rdR1Vb/59J33ERVk0cGK8J/eN0Ua8YhQT5m1gpgJa05vO7Q95J2r17t5KSkhQcHKz09HQtXLhQPXr00ObNm1VXV6eJEye62/bv3189evRwh6TMzEwNGTLEHZAkafLkyZozZ46ys7N14YUXtvg7a2pqVFPzn0eGnU7n+esgAK/VJSxIN45M1o0jk1VRU69VOUe0LNuuFbuKVFRWo7c25OmtDXmyBvjp4t4xmjAgTpf3j1NiZIjZpQNopQ4XkkaPHq3FixerX79+Kigo0OOPP66xY8dq+/btstvtCgoKUlRUVLOfiY+Pl91ulyTZ7fZmAanpfNO5U1m4cKEef/zx9u0MgA4tzBqgqceffqupb9C6vcVasbNIK3YV6XBplVbsavxakgYm2tyBaVh3huWAjqDDhaQrrrjC/fXQoUM1evRopaSk6L333lNIyPn7l9qDDz6o++67z/290+lUcnLyeft9ADoWa4C/xveL0/h+cfqNYSinsEzLdxZp+c5CbTlYqh0FTu0ocOr5FXvUNSxIl/SJ0aV9YzW2T6xiI6xmlw+gBR0uJH1fVFSU+vbtqz179uiHP/yhamtrVVpa2uxuUmFhoRISEiRJCQkJ2rhxY7P3aHr6ralNS6xWq6xW/kcG4MwsFov6J9jUP8GmueN7q7i8RqtyjmjFriKt+e6Iiitq9a+sfP0rK1+SNCjJpnF9YzWuT6xGpHRRUECHe/AY8EkdPiSVl5dr7969uvXWWzVixAgFBgZq+fLlmj59uiQpJydHeXl5Sk9PlySlp6friSeeUFFRkeLi4iRJGRkZstlsGjhwoGn9AOC7uoZbNX1Ed00f0V219S59k3dMa747ojW7j2j7Yaey8xuPRav2KizIX+m9YnRp3xhd1DtGaTFh7CsHmKTDPd32i1/8QldffbVSUlKUn5+vRx99VFlZWdqxY4diY2M1Z84cffLJJ1q8eLFsNpvuueceSdK6deskNS4BcMEFFygpKUlPPvmk7Ha7br31Vt15552tWgKAp9sAtIcjZTX6as8RrfnuqL7cfURHy2ubnU+MDFZ6r666uFeMLu4do4TIYJMqBXyDT6+4PWPGDK1Zs0bFxcWKjY3VJZdcoieeeEK9evWS1LiY5P3336+3335bNTU1mjx5sl588cVmQ2kHDhzQnDlztGrVKoWFhWnWrFn63e9+12y9pTMhJAFoby5X4ya8q787oq92H9XmvGPNlhiQpLTYsOOBqavGpHVVVOjZLVsCoJFPhyRvQUgCcL5V1zXo6/3HtHbvUa3bc1TbDjt0wk4pslik/gk2jU6N1pi0aP0gtauiwwhNwOkQkjyAkATA0xxVdVq/r1jr9hzV2r3F2lNUflKbvvHhGp3aVaPTojU6tStPzgHfQ0jyAEISALMVOau1IbdEG3KLtTG3RN8Vnhya0mLDNColWiN7dtHIntHq2TWUieDo1AhJHkBIAuBtistrtGl/idbvK9GG3JLjG/A2bxMTHqThPbpoZM8uGpESrSHdIllyAJ0KIckDCEkAvJ2jsk4b95fo6wMl2rz/mLYecqi2oflEcGuAn4Z1j9LwlC66sEeULkyOUpyNJ+jguwhJHkBIAtDR1NQ3aPthhzbtP6av9x/T5gMlOlZZd1K7blEhuuB4YLqwRxcNSrIpOJDNeuEbCEkeQEgC0NEZhqF9Ryv09f4SZR0s1Za8UuUUlp00RBfob9HARJsuSI7S0O5RGpYcqdSYcPmz/xw6IEKSBxCSAPiisuo6bTvk0JbjoSnr4LGTFriUpLAgfw3uFqlhyVEa0i1Sw7pHKTk6hEnh8HqEJA8gJAHoDAzD0KFjVfomr3FO09ZDpdp+2KmquoaT2kaFBmpIt0gN7hapwUmRGpRkU4/oUPlxxwlehJDkAYQkAJ1VfYNLe49U6NtDpdp6qFRbDzm0s8CpuoaTP04irAEamGTToKRIDe5m0+BukUqLCVOAP0/UwRyEJA8gJAHAf9TUNyjHXqZthx2NG/YedminveykbVUkKTjQT/3iIzQg0eY++idGyBYcaELl6GwISR5ASAKA06trcGlPUbmy853aftih7HyHduQ7VVF78lCdJHXvEuIOTQMTI9QvoXG4jgniaE+EJA8gJAFA67lchvYXV2hnQZl2FjjdR76jusX2wYF+6hMXoX4JEeoXH6G+CRHqnxChuAgrk8RxTghJHkBIAoD2U1pZ2zw42Z3aXViumhaG66TGSeJ94yPUNz5cfeMj1DsuXH3iIhQTHkR4wmkRkjyAkAQA51eDy1BeSaVy7E7l2MuVU+hUjr1MuUcr5DrFJ1dUaKD6xkWod3y4+hwPTr3jwhVv484TGhGSPICQBADmqK5r0N4j5cqxl2l3Ubl2F5ZrT1GZDpRUnrQQZpNwa4B6xYYpLTZcvWLD1Cs2XL3iwpXSNVTWAFYT70wISR5ASAIA79IUnvYcD067i8q0u7BcB0oq1XCKW09+FqlHdKjSYsOVFhOm1NgwpcY0Hgm2YO4++SBCkgcQkgCgY6itdymvpEJ7iiq090j58aNC+4rKVVZTf8qfCwn0V8+YsMbwdPzoGROqnl3DFB3G3KeOqjWf3wEeqgkAAFMEBfipd1yEesdFNHvdMAwdKavRnuOhaf/RCuUeP/JKKlVV1+CeSP59EdYApcSEKqVrmHp2bfqz8etYnrzzGdxJOkfcSQIA31XX4NLBkkp3aNp3tEK5RxrDU76j6pRzn6TGO1DJ0SHqER2mHtGh6hEdopSuYUqODlX3LiEKDmQOlJm4kwQAQBsE+vs1zlOKDT/pXHVdgw4dq1Tu0UodKK7Q/uIKHSiu1P7iCh0+VqWqugZ9V1iu7wrLW3zvBFuwkqNDlNwlVN2jQ5XcJUTdu4QqOTpEiZEhLJ7pRbiTdI64kwQA+L7aepfyS6uUV1KpAyWVOlhSqbzixq/ziitOudp4kwA/i5KiQtS9S+PRLSpU3dxfhyghMliB7HvXJtxJAgDABEEBfuoZE6aeMWEnnTMMQ8cq63SguEIHj1XpYEmlDh2r0qFjjWHqcGmV6hoa14bKK6ls8f39LI13orodD03duoQoKSpESZHH/4wKVgR74LUbQhIAAB5gsVgUHRak6LAgXdijy0nnXS5DhWXVOlhS5Q5Nh49VNf55/OvaBpfyHdXKd1Rrk461+HsirAFKigpRYlTw8QAVrMTIECVGBisxKkQJtmCFBDEv6mww3HaOGG4DAHiSy2XoaHmNDp0Yno5VqcBRpfzSauU7qlRaWXdW7xUVGvif4HT8iLcFKyEyWAm2YMVHBivCGuCTT+kx3AYAgI/x87MozhasOFuwhrdwJ0qSKmvrGwNTaWN4Onz8a7ujWgWOKhU4qlVZ26DSyjqVVta1uLxBk9Ag/8bAdDw8xduCFW+zuv+MiwhWnM3q0yuWE5IAAPARoUEB6h0Xrt5xJz+VJzXOi3JW18vuaLzz1BieqlVQWiW7s1qFzmrZHdVyVtersrZB+44vf3A6XUIDFX88vMVFWBVvsyo23Or+PjaiMVB1xCE+QhIAAJ2ExWJRZEigIkMC1S8h4pTtKmvrVeiskd1xPDgdD09FZdUqdNao0FmtImeNahtcOlZZp2OVddplLzvt746wBii2KTTZghUTHtT4fbhVMcf/jIuwKjosSAFe8gQfIQkAADQTGhSg1JgApbbwlF4TwzBUWlmnorLG0NR0HCmrUdHxo/HralXXuVRWU6+ymvoz3pmyWKTo0MYAdfPoHpqZ3rOde3f2CEkAAKDVLBaLuoQFqUtY0GnvShmGobKa+sbA5GwMTUfKanS0vFZHymp0pLxGR4//WVxeI5chFVfUqriiVmXVp95bzxMISQAA4LyxWCyyBQfKFhyoXi2sYH6iBpehY5W1x0NUjZK7hHqoypYRkgAAgFfw97MoJtyqmHCr2aVIkrxjZhQAAICXISQBAAC0gJAEAADQAkISAABACwhJAAAALSAkAQAAtICQBAAA0AJCEgAAQAsISQAAAC0gJAEAALSAkAQAANACQhIAAEALCEkAAAAtCDC7gI7KMAxJktPpNLkSAABwtpo+t5s+x0+HkHSOysrKJEnJyckmVwIAAFqrrKxMkZGRp21jMc4mSuEkLpdL+fn5ioiIkMViMbucVnE6nUpOTtbBgwdls9nMLscj6DN99lX0mT77ovPZX8MwVFZWpqSkJPn5nX7WEXeSzpGfn5+6d+9udhltYrPZOsVfthPR586BPncO9Nn3na/+nukOUhMmbgMAALSAkAQAANACQlInZLVa9eijj8pqtZpdisfQ586BPncO9Nn3eUt/mbgNAADQAu4kAQAAtICQBAAA0AJCEgAAQAsISQAAAC0gJPmwNWvW6Oqrr1ZSUpIsFos+/PDDZudvu+02WSyWZseUKVPMKbYdLFy4UKNGjVJERITi4uI0bdo05eTkNGtTXV2tuXPnqmvXrgoPD9f06dNVWFhoUsVtdzZ9vuyyy066znfffbdJFbfdokWLNHToUPcic+np6fr000/d533tGktn7rOvXeOW/O53v5PFYtH8+fPdr/nitW7SUn998To/9thjJ/Wpf//+7vNmX2NCkg+rqKjQsGHD9MILL5yyzZQpU1RQUOA+3n77bQ9W2L5Wr16tuXPnav369crIyFBdXZ0mTZqkiooKd5t7771XH330kd5//32tXr1a+fn5uv76602sum3Ops+SNHv27GbX+cknnzSp4rbr3r27fve732nz5s36+uuvdfnll+vaa69Vdna2JN+7xtKZ+yz51jX+vk2bNukvf/mLhg4d2ux1X7zW0qn7K/nmdR40aFCzPn311Vfuc6ZfYwOdgiRjyZIlzV6bNWuWce2115pSjycUFRUZkozVq1cbhmEYpaWlRmBgoPH++++72+zcudOQZGRmZppVZrv6fp8NwzAuvfRS4+c//7l5RXlAly5djFdeeaVTXOMmTX02DN++xmVlZUafPn2MjIyMZv301Wt9qv4ahm9e50cffdQYNmxYi+e84RpzJ6mTW7VqleLi4tSvXz/NmTNHxcXFZpfUbhwOhyQpOjpakrR582bV1dVp4sSJ7jb9+/dXjx49lJmZaUqN7e37fW7y5ptvKiYmRoMHD9aDDz6oyspKM8prdw0NDXrnnXdUUVGh9PT0TnGNv9/nJr56jefOnaupU6c2u6aS7/59PlV/m/jidd69e7eSkpKUlpamW265RXl5eZK84xqzwW0nNmXKFF1//fVKTU3V3r179dBDD+mKK65QZmam/P39zS6vTVwul+bPn6+LL75YgwcPliTZ7XYFBQUpKiqqWdv4+HjZ7XYTqmxfLfVZkm6++WalpKQoKSlJW7du1YIFC5STk6MPPvjAxGrbZtu2bUpPT1d1dbXCw8O1ZMkSDRw4UFlZWT57jU/VZ8k3r7EkvfPOO/rmm2+0adOmk8754t/n0/VX8s3rPHr0aC1evFj9+vVTQUGBHn/8cY0dO1bbt2/3imtMSOrEZsyY4f56yJAhGjp0qHr16qVVq1ZpwoQJJlbWdnPnztX27dubjW37ulP1+a677nJ/PWTIECUmJmrChAnau3evevXq5eky20W/fv2UlZUlh8Ohf/7zn5o1a5ZWr15tdlnn1an6PHDgQJ+8xgcPHtTPf/5zZWRkKDg42Oxyzruz6a8vXucrrrjC/fXQoUM1evRopaSk6L333lNISIiJlTViuA1uaWlpiomJ0Z49e8wupU3mzZunpUuXauXKlerevbv79YSEBNXW1qq0tLRZ+8LCQiUkJHi4yvZ1qj63ZPTo0ZLUoa9zUFCQevfurREjRmjhwoUaNmyYnn32WZ++xqfqc0t84Rpv3rxZRUVFGj58uAICAhQQEKDVq1frueeeU0BAgOLj433qWp+pvw0NDSf9jC9c5++LiopS3759tWfPHq/4+0xIgtuhQ4dUXFysxMREs0s5J4ZhaN68eVqyZIlWrFih1NTUZudHjBihwMBALV++3P1aTk6O8vLyms3t6EjO1OeWZGVlSVKHvc4tcblcqqmp8clrfCpNfW6JL1zjCRMmaNu2bcrKynIfI0eO1C233OL+2peu9Zn629IUCF+4zt9XXl6uvXv3KjEx0Tv+PntkejhMUVZWZmzZssXYsmWLIcl4+umnjS1bthgHDhwwysrKjF/84hdGZmamkZuba3zxxRfG8OHDjT59+hjV1dVml35O5syZY0RGRhqrVq0yCgoK3EdlZaW7zd1332306NHDWLFihfH1118b6enpRnp6uolVt82Z+rxnzx7jN7/5jfH1118bubm5xr/+9S8jLS3NGDdunMmVn7tf/epXxurVq43c3Fxj69atxq9+9SvDYrEYn3/+uWEYvneNDeP0ffbFa3wq33+6yxev9YlO7K+vXuf777/fWLVqlZGbm2usXbvWmDhxohETE2MUFRUZhmH+NSYk+bCVK1cakk46Zs2aZVRWVhqTJk0yYmNjjcDAQCMlJcWYPXu2YbfbzS77nLXUV0nGa6+95m5TVVVl/PSnPzW6dOlihIaGGtddd51RUFBgXtFtdKY+5+XlGePGjTOio6MNq9Vq9O7d23jggQcMh8NhbuFtcPvttxspKSlGUFCQERsba0yYMMEdkAzD966xYZy+z754jU/l+yHJF6/1iU7sr69e55tuuslITEw0goKCjG7duhk33XSTsWfPHvd5s6+xxTAMwzP3rAAAADoO5iQBAAC0gJAEAADQAkISAABACwhJAAAALSAkAQAAtICQBAAA0AJCEgAAQAsISQAAAC0gJAEAALSAkAQAANACQhKATum2226TxWKRxWJRYGCgUlNT9ctf/lLV1dVmlwbASwSYXQAAmGXKlCl67bXXVFdXp82bN2vWrFmyWCz6/e9/b3ZpALwAd5IAdFpWq1UJCQlKTk7WtGnTNHHiRGVkZJhdFgAvQUgCAEnbt2/XunXrFBQUZHYpALwEw20AOq2lS5cqPDxc9fX1qqmpkZ+fn/785z+bXRYAL0FIAtBpjR8/XosWLVJFRYX+9Kc/KSAgQNOnT2/W5sMPP9Qrr7yi2tpa3XTTTbrjjjtMqhaApzHcBqDTCgsLU+/evTVs2DC9+uqr2rBhg/72t7+5z7/55pt67733tGjRIr3xxhvasWOHnnjiCRMrBuBJhCQAkOTn56eHHnpIDz/8sKqqqiRJL7/8sl5//XUlJycrLi5OTz31lFatWqWysjKTqwXgCYQkADjuxhtvlL+/v1544QUVFxerR48eCgwM1CuvvKLbbrtNkjR69Gh999135hYKwCMISQBwXEBAgObNm6cnn3xSwcHBKigokNQYnhYuXChJysnJUXJyspllAvAQi2EYhtlFAIA3euSRRxQVFaX77rtPkvTee+/pk08+0eLFi80tDIBHEJIA4BTq6+v18MMP69NPP5XFYtGoUaP0zDPPKCwszOzSAHgAIQkAAKAFzEkCAABoASEJAACgBYQkAACAFhCSAAAAWkBIAgAAaAEhCQAAoAWEJAAAgBYQkgAAAFpASAIAAGgBIQkAAKAF/x+wENi6jRw7JAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "with quantity_support():\n", + " plt.figure()\n", + " plt.plot(radii.to(u.R_sun), n_e)\n", + " plt.draw()" + ] + }, { "cell_type": "markdown", "id": "f4da2ab0", @@ -665,12 +797,42 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 63, "id": "e3bc2348", "metadata": {}, "outputs": [], "source": [ - "volume = 0.62 * (u.barn * u.Mpc)" + "volume = 1.6 * (u.barn * u.Mpc)" + ] + }, + { + "cell_type": "markdown", + "id": "6ac2b3db-2c6d-486f-8f1c-baaa563c710f", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": 64, + "id": "2a109221-baf0-495d-b914-e40dfa5ed75b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$1.001656 \\; \\mathrm{tsp}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 64, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "volume.to(u.imperial.tsp)" ] }, {