diff --git a/docs/src/tutorials.md b/docs/src/tutorials.md index 17077a5b..9eab0292 100644 --- a/docs/src/tutorials.md +++ b/docs/src/tutorials.md @@ -1,5 +1,4 @@ # [Tutorial notebooks](@id tutorials) Here are a few tutorials in the form of Jupyter Notebooks to help you get started using Korg. Each of the following bullets link to a folder on GitHub containing the associated Julia and Python notebooks (the `.ipynb` files) and any supporting data. You can download everything if you want to run the notebook yourself, or just view the notebook on GitHub. - [The Basics](https://github.com/ajwheeler/Korg.jl/tree/main/misc/Tutorial%20notebooks/basics) demonstrates how to read linelists, read or interpolate model atmospheres, and synthesize spectra. It also goes into some of the other data returned when you synthesize a spectrum. -- [Fitting](https://github.com/ajwheeler/Korg.jl/tree/main/misc/Tutorial%20notebooks/fitting) demonstrates the use of [`Korg.Fit.fit_spectrum`](@ref) to fit observational data. -- [EW fitting](https://github.com/ajwheeler/Korg.jl/tree/main/misc/Tutorial%20notebooks/EW%20fitting) demonstrates the use of [`Korg.Fit.ews_to_abundances`](@ref) to calculate abundances from equivalent widths. \ No newline at end of file +- [Fitting](https://github.com/ajwheeler/Korg.jl/tree/main/misc/Tutorial%20notebooks/fitting) demonstrates the use of [`Korg.Fit.fit_spectrum`](@ref) to fit observational data. \ No newline at end of file diff --git a/misc/Tutorial notebooks/EW fitting/EW fitting in python.ipynb b/misc/Tutorial notebooks/EW fitting/EW fitting in python.ipynb deleted file mode 100644 index 4c043d12..00000000 --- a/misc/Tutorial notebooks/EW fitting/EW fitting in python.ipynb +++ /dev/null @@ -1,648 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "b1c9787e-781e-439f-9599-ea5c0d2f644b", - "metadata": {}, - "source": [ - "# Equivalent width fitting\n", - "This notebook is an example of how to go from measured equivalent widths to abundances. We'll use the data from [Melendez 2014](https://ui.adsabs.harvard.edu/abs/2014ApJ...791...14M/abstract) and redo their differential analysis of iron lines 18 Sco. \n", - "\n", - "A massive thanks to Emily Griffith for helping translate this to python." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "25fa263f-d156-4a4a-ab1e-02bea37f9d9b", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:22.850731Z", - "iopub.status.busy": "2024-06-26T19:09:22.850385Z", - "iopub.status.idle": "2024-06-26T19:09:26.233538Z", - "shell.execute_reply": "2024-06-26T19:09:26.233141Z", - "shell.execute_reply.started": "2024-06-26T19:09:22.850700Z" - } - }, - "outputs": [], - "source": [ - "from matplotlib import pyplot as plt\n", - "from juliacall import Main as jl\n", - "import time\n", - "import numpy as np\n", - "import pandas as pd\n", - "\n", - "jl.seval(\"using Korg\")\n", - "Korg = jl.Korg" - ] - }, - { - "cell_type": "markdown", - "id": "988c6101-cff0-4320-a444-954379815fdc", - "metadata": {}, - "source": [ - "## parsing the linelist" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "b9e28b1e-d8e8-4749-b5e0-cc8decba9284", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:26.234343Z", - "iopub.status.busy": "2024-06-26T19:09:26.234216Z", - "iopub.status.idle": "2024-06-26T19:09:26.245501Z", - "shell.execute_reply": "2024-06-26T19:09:26.245240Z", - "shell.execute_reply.started": "2024-06-26T19:09:26.234335Z" - } - }, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
wlspeciesExPotlog_gfC6EW_18ScoEW_Sun
05044.21126.02.8512-2.0582.710000e-3174.874.3
15054.64226.03.6400-1.9214.680000e-3240.940.5
25127.35926.00.9150-3.3071.840000e-3297.596.1
35127.67926.00.0520-6.1251.200000e-3218.919.1
45198.71126.02.2230-2.1354.610000e-3299.398.1
........................
3183694.81066.10.1030-0.1102.800000e+0017.313.9
3194077.97066.10.1030-0.0402.800000e+0014.810.7
3204103.31066.10.1030-0.3802.800000e+0017.815.0
3214449.70066.10.0000-1.0302.800000e+004.33.3
3223694.19070.10.0000-0.3002.800000e+0081.273.9
\n", - "

323 rows × 7 columns

\n", - "
" - ], - "text/plain": [ - " wl species ExPot log_gf C6 EW_18Sco EW_Sun\n", - "0 5044.211 26.0 2.8512 -2.058 2.710000e-31 74.8 74.3\n", - "1 5054.642 26.0 3.6400 -1.921 4.680000e-32 40.9 40.5\n", - "2 5127.359 26.0 0.9150 -3.307 1.840000e-32 97.5 96.1\n", - "3 5127.679 26.0 0.0520 -6.125 1.200000e-32 18.9 19.1\n", - "4 5198.711 26.0 2.2230 -2.135 4.610000e-32 99.3 98.1\n", - ".. ... ... ... ... ... ... ...\n", - "318 3694.810 66.1 0.1030 -0.110 2.800000e+00 17.3 13.9\n", - "319 4077.970 66.1 0.1030 -0.040 2.800000e+00 14.8 10.7\n", - "320 4103.310 66.1 0.1030 -0.380 2.800000e+00 17.8 15.0\n", - "321 4449.700 66.1 0.0000 -1.030 2.800000e+00 4.3 3.3\n", - "322 3694.190 70.1 0.0000 -0.300 2.800000e+00 81.2 73.9\n", - "\n", - "[323 rows x 7 columns]" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# use pandas to read the contents of Table 1, which were downloaded from the journal website\n", - "lines = pd.read_csv(\"Table 1.dat\", delim_whitespace=True, skiprows=24,\n", - " names=[\"wl\", \"species\", \"ExPot\", \"log_gf\", \"C6\", \"EW_18Sco\", \"EW_Sun\"])\n", - "lines" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "7e400290-7608-473c-af8e-9686e8fb9a61", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:26.246109Z", - "iopub.status.busy": "2024-06-26T19:09:26.246021Z", - "iopub.status.idle": "2024-06-26T19:09:27.536472Z", - "shell.execute_reply": "2024-06-26T19:09:27.536042Z", - "shell.execute_reply.started": "2024-06-26T19:09:26.246100Z" - } - }, - "outputs": [], - "source": [ - "# convert the numbers in the species column to Korg.Species objects\n", - "# jl.broadcast is how you use Julia's broadcasting (\".\" syntax) from python\n", - "# this called Korg.Species on each element in the line.species column\n", - "lines.species = jl.broadcast(Korg.Species, lines.species)\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "41af8329-f8cc-47bc-a55d-2ef9266bc52b", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:27.538462Z", - "iopub.status.busy": "2024-06-26T19:09:27.538349Z", - "iopub.status.idle": "2024-06-26T19:09:28.064649Z", - "shell.execute_reply": "2024-06-26T19:09:28.064350Z", - "shell.execute_reply.started": "2024-06-26T19:09:27.538453Z" - } - }, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
wlspeciesExPotlog_gfC6EW_18ScoEW_Sun
05044.211Fe I2.8512-2.0582.710000e-3174.874.3
15054.642Fe I3.6400-1.9214.680000e-3240.940.5
25127.359Fe I0.9150-3.3071.840000e-3297.596.1
35127.679Fe I0.0520-6.1251.200000e-3218.919.1
865197.577Fe II3.2306-2.2208.690000e-3384.080.7
........................
816810.263Fe I4.6070-0.9864.500000e-3153.150.6
826837.006Fe I4.5930-1.6872.460000e-3220.018.9
836839.830Fe I2.5590-3.3503.950000e-3232.531.6
846843.656Fe I4.5480-0.8602.940000e-3164.163.3
856858.150Fe I4.6070-0.9303.240000e-3153.652.4
\n", - "

98 rows × 7 columns

\n", - "
" - ], - "text/plain": [ - " wl species ExPot log_gf C6 EW_18Sco EW_Sun\n", - "0 5044.211 Fe I 2.8512 -2.058 2.710000e-31 74.8 74.3\n", - "1 5054.642 Fe I 3.6400 -1.921 4.680000e-32 40.9 40.5\n", - "2 5127.359 Fe I 0.9150 -3.307 1.840000e-32 97.5 96.1\n", - "3 5127.679 Fe I 0.0520 -6.125 1.200000e-32 18.9 19.1\n", - "86 5197.577 Fe II 3.2306 -2.220 8.690000e-33 84.0 80.7\n", - ".. ... ... ... ... ... ... ...\n", - "81 6810.263 Fe I 4.6070 -0.986 4.500000e-31 53.1 50.6\n", - "82 6837.006 Fe I 4.5930 -1.687 2.460000e-32 20.0 18.9\n", - "83 6839.830 Fe I 2.5590 -3.350 3.950000e-32 32.5 31.6\n", - "84 6843.656 Fe I 4.5480 -0.860 2.940000e-31 64.1 63.3\n", - "85 6858.150 Fe I 4.6070 -0.930 3.240000e-31 53.6 52.4\n", - "\n", - "[98 rows x 7 columns]" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# look at Fe lines only, and sort by wavelength\n", - "is_iron = [Korg.get_atoms(s) == jl.Vector([26]) for s in lines.species]\n", - "lines = lines.loc[is_iron]\n", - "lines = lines.sort_values(\"wl\")\n", - "lines" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "5b75be64-3518-45de-8f0d-9fd08ef6d2be", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:28.065169Z", - "iopub.status.busy": "2024-06-26T19:09:28.065088Z", - "iopub.status.idle": "2024-06-26T19:09:28.180243Z", - "shell.execute_reply": "2024-06-26T19:09:28.179886Z", - "shell.execute_reply.started": "2024-06-26T19:09:28.065162Z" - } - }, - "outputs": [], - "source": [ - "# to pass the lines to Korg, we need to turn each row of the \"lines\" dataframe into a Korg.Line object\n", - "linelist = jl.broadcast(Korg.Line, lines.wl, # can be in either cm or Å (like these), but NOT nm\n", - " lines.log_gf,lines.species, # needs to be a Korg.Species, which we handled in the cell above\n", - " lines.ExPot)" - ] - }, - { - "cell_type": "markdown", - "id": "9e0db71d-1010-4a39-918c-bd1ce71c2069", - "metadata": {}, - "source": [ - "## Set up params for each star" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "69f5b674-9538-4fd8-b95e-ae6e17a0f9d8", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:28.180826Z", - "iopub.status.busy": "2024-06-26T19:09:28.180727Z", - "iopub.status.idle": "2024-06-26T19:09:36.323403Z", - "shell.execute_reply": "2024-06-26T19:09:36.323100Z", - "shell.execute_reply.started": "2024-06-26T19:09:28.180819Z" - } - }, - "outputs": [], - "source": [ - "# solar params\n", - "sun_Teff, sun_logg, sun_Fe_H, sun_vmic = 5777, 4.44, 0.0, 1.0\n", - "\n", - "# vector of abundances for the sun\n", - "sun_A_X = Korg.format_A_X(sun_Fe_H)\n", - "\n", - "# interpolate a model atmosphere for the sun\n", - "sun_atm = Korg.interpolate_marcs(sun_Teff, sun_logg, sun_A_X) \n", - "\n", - "# and likewise for 18 Sco\n", - "sco_teff, sco_logg, sco_fe_h, sco_vmic = (5823, 4.45, 0.054, sun_vmic + 0.02)\n", - "sco_A_X = Korg.format_A_X(sco_fe_h)\n", - "sco_atm = Korg.interpolate_marcs(sco_teff, sco_logg, sco_A_X)" - ] - }, - { - "cell_type": "markdown", - "id": "e1e6a43b-2d28-4d03-b5d3-3a1e0ef6d135", - "metadata": {}, - "source": [ - "## Calculate the abundances for each star\n", - "Because Julia uses \"just in time\" compilation, the first call to `Korg.Fit.ews_to_abundances` will take several seconds. The second is much faster, because the code is already compiled. If you rerun the cell below, both calculations will be fast, even though no data is cached." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "fc38fcb2-2bae-4fe8-af6f-c3604cc05832", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:36.323891Z", - "iopub.status.busy": "2024-06-26T19:09:36.323822Z", - "iopub.status.idle": "2024-06-26T19:09:44.196202Z", - "shell.execute_reply": "2024-06-26T19:09:44.195932Z", - "shell.execute_reply.started": "2024-06-26T19:09:36.323884Z" - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "7.6071089999999995\n", - "0.276288000000001\n" - ] - } - ], - "source": [ - "# calculate abundances from the EWs for each star\n", - "\n", - "# in order to pass the EWs to Korg.Fit.ews_to_abundances, we convert them to numpy arrays, which \n", - "# are handled automatically, unlike pandas series. We also wrap the returned julia vector in np.array.\n", - "t1 = time.process_time()\n", - "lines[\"A_sun\"] = np.array(Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, np.array(lines.EW_Sun), vmic=sun_vmic))\n", - "t2 = time.process_time()\n", - "lines[\"A_18Sco\"] = np.array(Korg.Fit.ews_to_abundances(sco_atm, linelist, sco_A_X, np.array(lines.EW_18Sco), vmic=sco_vmic))\n", - "t3 = time.process_time()\n", - "print(t2-t1)\n", - "print(t3-t2)" - ] - }, - { - "cell_type": "markdown", - "id": "31560282-f1b1-4538-a3c0-cc1961b52ecb", - "metadata": {}, - "source": [ - "# Plot the results\n", - "Plotting as a function of wavelength and equivalenth width is left as an excercise for the reader." - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "id": "836607ca-df58-4513-816a-33e09aa87739", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:44.196685Z", - "iopub.status.busy": "2024-06-26T19:09:44.196604Z", - "iopub.status.idle": "2024-06-26T19:09:44.198455Z", - "shell.execute_reply": "2024-06-26T19:09:44.198204Z", - "shell.execute_reply.started": "2024-06-26T19:09:44.196678Z" - } - }, - "outputs": [], - "source": [ - "# get a bitmask for the lines of Fe I vs Fe II\n", - "neutrals = np.array([spec.charge == 0 for spec in lines.species])" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "e6dd0392-6ed9-4640-b22c-b67383a19910", - "metadata": { - "execution": { - "iopub.execute_input": "2024-06-26T19:09:44.198824Z", - "iopub.status.busy": "2024-06-26T19:09:44.198758Z", - "iopub.status.idle": "2024-06-26T19:09:44.279104Z", - "shell.execute_reply": "2024-06-26T19:09:44.278737Z", - "shell.execute_reply.started": "2024-06-26T19:09:44.198818Z" - } - }, - "outputs": [ - { - "data": { - "text/plain": [ - "''" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAADZCAYAAADCHOODAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABF0UlEQVR4nO3deVxU9foH8M+ALIKsKmuI4IIiKrmB5pJb4AZmedPrFTWv3hbLtLxqdVXSUjNNTdO6ueZVy7JE83KvornST0NxidTkkriAe4yggs7M749pyIGZ4ZyZMzNnhs/79eKlzJw5PANzZp7zPc/3+So0Go0GIuzbtw9dunSBm5ubmIcREREREdmFQmzC+6j79++joqJC7zZfX1+LgyIiIiIikoqL2AfcvXsXEyZMQFBQELy9vREQEKD3RUREREQkJ6IT3ilTpmDPnj1YsWIFPDw88NlnnyE9PR1hYWFYv369NWIkIiIiIjKb6JKGRo0aYf369XjyySfh6+uLY8eOoWnTpvj888+xadMm7Ny501qxEhERERGJJnqE99atW4iOjgagrde9desWAKBr167Yv3+/tNEREREREVlIdMIbHR2NgoICAECLFi3w5ZdfAgC2b98Of39/SYMjIiIiIrKU6JKGDz/8EK6urnj11Vexe/duDBo0CBqNBg8ePMCiRYswceJEa8VKRERERCSaRW3JAODChQvIyclB06ZN0aZNG6niIiIiIiKShMUJrzNSq9W4cuUKfHx8oFAo7B0OEREREVWh0Whw584dhIWFwcXFdJVuHXN+wNGjR7F3715cu3YNarVa775FixaZs0tZuXLlCiIiIuwdBhERERHV4OLFi3jsscdMbiM64X3vvffw9ttvIyYmBsHBwXojoM4yGurj4wNA+wvkynFERERE8qNUKhEREVGZt5kiOuFdsmQJVq9ejdGjR5sTm0PQJe6+vr5MeKmSSq3BkYJbuHbnPoJ8PNEpKhCuLs5xkkdEROSohAy4ik54XVxc8MQTT5gVEJGjyjxdhPTteSgquV95W6ifJ2YOikVyXKgdIyMiIqKaiO7DO2nSJCxfvtwasRDJUubpIry44ZhesgsAxSX38eKGY8g8XWSnyIiIiEgI0V0a1Go1BgwYgHPnziE2NhZubm5692/dulXSAO1BqVTCz88PJSUlLGmo5VRqDbrO31Mt2dVRAAjx88TBqb1Y3kBERGRDYvI10SUNr776Kvbu3YuePXuifv36TjNRjciQIwW3jCa7AKABUFRyH0cKbqFzk/q2C4yIiBySRqPBw4cPoVKp7B2K7Lm6uqJOnTqS5JqiE95169bh66+/xoABAyz+4URyd+2O8WTXnO2IiKj2qqioQFFREe7evWvvUByGl5cXQkND4e7ubtF+RCe8gYGBaNKkiUU/lMhRBPl4SrodERHVTmq1GgUFBXB1dUVYWBjc3d15ldwEjUaDiooKXL9+HQUFBWjWrFmNi0uYIjrhnTVrFmbOnIk1a9bAy8vL7B9M5Ag6RQUi1M8TxSX3YajYXVfD2ykq0NahERGRA6moqIBarUZERATzJ4Hq1q0LNzc3XLhwARUVFfD0NH9wSXTCu3TpUuTn5yM4OBiNGzeuNmnt2LFjZgdDJDeuLgrMHBSLFzccgwLQS3p15+UzB8U67YQ19h4mIpKWJaOUtZFUvy/RCe/gwYMl+cFEjiI5LhQr/tKuWh/eECfvw8vew0RE5CxEtyWrSUVFBW7cuIGwsDApd2tTbEtGhtSm0U5d7+Gqbw66Z7viL+2Y9BIRiXD//n0UFBQgKirKokvztY2p35uYfE2SceKbN29i3bp1eOaZZ9CgQQNs3LhRit0SyYqriwKdm9RHanw4Ojep77TJrkqtQfr2PIM1y7rb0rfnQaWW9FyZiIjIasxOeM+dO4cPPvgA3bp1Q3h4OFatWoXOnTvjxx9/xBtvvCFljERkQ2J6DxMRkfMbPXo0FApFta/z58+btb9Zs2YhPj5e2iBrILiGV6PR4PDhw8jIyMC2bdtw+fJl9OnTB2PHjsU333yDBg0aWDNOIrIR9h4mIpI3e5TYJScnY82aNXq3NWzY0Ko/U0qCE97g4GDUqVMHAwYMwAcffIA+ffqwBoXICbH3MBGRfNlrQrGHhwdCQkIM3rdt2zakp6cjLy8PYWFhGDVqFN566y3UqSO6N4LVCI4kIyMDCQkJgpskDxgwAJ999hlCQzmxhciRsPcwEZE8GZtQXFxyHy9uOGaXCcUHDhxAWloali5dim7duiE/Px/jx48HAMycOdOmsZgiuIY3MTFR1Iog+/fvx71798wKiojsR9d7GPijK4NObeg9TEQkR/aeULxjxw7Uq1ev8mvo0KHan5mejmnTpmHUqFGIjo5G3759MXv2bHzyySdWicNc8hlrJiLZqK29h4nI9mpTy0dLiJlQ3LlJfcl/fs+ePbFixYrK7729vQEAJ06cwKFDh/Duu+9W3qdSqXD//n3cvXtXNqvKMeElIoOS40LRNzaEH0REZDVc4EY4e08o9vb2RtOmTavdXlpaivT0dAwZMqTafXKa68WEl4iM0vUeJiKSmhzrUeVMrhOK27Vrh7NnzxpMhuWECS8RERHZVE31qApo61H7xobwqtLv5DqheMaMGRg4cCAaNWqEZ599Fi4uLjhx4gROnz6NOXPm2DQWUyRZaY2IiIhIKC5wI55cJxQnJSVhx44d+O9//4uOHTsiMTERH374ISIjI20aR02sNsL75ptvIjCQbYuIiIhIn73rUR2VvSYUr1271uT9SUlJSEpKEry/WbNmYdasWZYFJZLohLe8vBz/93//hwsXLuDu3bto2LAhHn/8cURFReltN336dMmCJCIiIuch13pUR8AJxeYRnPAeOnQIS5Yswfbt2/HgwQP4+fmhbt26uHXrFsrLyxEdHY3x48fjhRdegI+PjzVjJiIiIgcm13pUR8EJxeIJquFNSUnBc889h8aNG+O///0v7ty5g5s3b+LSpUu4e/cufvnlF7z99tvIyspC8+bNsWvXLmvHTURERA5KrvWo5LwEjfAOGDAAX3/9Ndzc3AzeHx0djejoaIwaNQp5eXkoKiqSNEgiIiJyLlzghmxJUML7t7/9TfAOY2NjERsba3ZAREREVDuwHpVsxawuDb/99hu++uor5OfnY8qUKQgMDMSxY8cQHByM8PBwqWMkIiIiJ8V6VLIF0QnvyZMn0adPH/j5+eHXX3/FuHHjEBgYiK1bt6KwsBDr16+3RpxERERERGYRvfDE5MmTMXr0aPzyyy96ayT3798f+/fvlzQ4IiIiIiJLiU54jx49arCmNzw8HMXFxWYFsXz5cjRu3Bienp5ISEjAkSNHTG6/ZcsWtGjRAp6enmjdujV27typd79CoTD4tWDBArPiIyIiIiLHJTrh9fDwgFKprHb7uXPn0LBhQ9EBfPHFF5g8eTJmzpyJY8eOoW3btkhKSsK1a9cMbn/48GEMHz4cY8eOxfHjxzF48GAMHjwYp0+frtymqKhI72v16tVQKBR45plnRMdHREREVJuNHj3a4EDi+fPnzdrfrFmzEB8fb/R7axCd8KakpOCdd97BgwcPAGhHUwsLCzF16lSzEspFixZh3LhxGDNmDGJjY7Fy5Up4eXlh9erVBrdfsmQJkpOTMWXKFLRs2RKzZ89Gu3btsGzZssptQkJC9L62bduGnj17Ijo6WnR8RERERLVdcnJytQHFqqvsypnohHfhwoUoLS1FUFAQ7t27hx49eqBp06bw8fHBu+++K2pfFRUVyMnJQZ8+ff4IyMUFffr0QXZ2tsHHZGdn620PaNdwNrb91atX8d1332Hs2LGiYiMiIiKSJbUKKDgAnPpK+69aZfUf6eHhUW1A0dXVFQCwbds2tGvXDp6enoiOjkZ6ejoePnxo9ZjEEN2lwc/PD7t27cKhQ4dw4sQJlJaWol27dtWSUCFu3LgBlUqF4OBgvduDg4Nx5swZg48pLi42uL2x+uF169bBx8cHQ4YMMRpHeXk5ysvLK783VLJBREREZHd5GUDmVEB55Y/bfMOA5PlAbIrNwzlw4ADS0tKwdOlSdOvWDfn5+Rg/fjwAYObMmTaPxxjBI7xpaWm4c+dO5ff16tXDuHHj8Pe//92sZNdWVq9ejREjRuh1lKhq7ty58PPzq/yKiIiwYYREREREAuRlAF+m6Se7AKAs0t6el2G1H71jxw7Uq1ev8mvo0KEAgPT0dEybNg2jRo1CdHQ0+vbti9mzZ+OTTz6xWizmEJzw/utf/8K9e/cqv+/WrRsuXrxo0Q9v0KABXF1dcfXqVb3br169ipCQEIOPCQkJEbz9gQMHcPbsWfz1r381Gcf06dNRUlJS+WXp8yIiIiKSlFqlHdmFxsCdv9+WOc1q5Q09e/ZEbm5u5dfSpUsBACdOnMA777yjlwyPGzcORUVFuHv3rlViMYfgkgaNRmPye3O4u7ujffv2yMrKwuDBgwEAarUaWVlZmDBhgsHHdO7cGVlZWXjttdcqb9u1axc6d+5cbdtVq1ahffv2aNu2rck4PDw84OHhYfbzICIiIrKqC4erj+zq0QDKy9rtorpJ/uO9vb3RtGnTareXlpYiPT3dYOmoqavrtmbW0sJSmjx5MkaNGoUOHTqgU6dOWLx4McrKyjBmzBgA2lKK8PBwzJ07FwAwceJE9OjRAwsXLsSAAQOwefNm/Pjjj/j000/19qtUKrFlyxYsXLjQ5s+JiIiISFKlV2veRsx2EmnXrh3Onj1rMBmWE1EJb15eXuXkMI1GgzNnzqC0tFRvmzZt2ogK4LnnnsP169cxY8YMFBcXIz4+HpmZmZUT0woLC+Hi8kflRZcuXbBx40a8/fbbePPNN9GsWTN8++23iIuL09vv5s2bodFoMHz4cFHxEBEREclOveCatxGznURmzJiBgQMHolGjRnj22Wfh4uKCEydO4PTp05gzZ45NYzFFVMLbu3dvvVKGgQMHAtD24tVoNFAoFFCpxNeOTJgwwWgJw/fff1/ttqFDh1YWSxszfvz4ylmCRERERA4tsou2G4OyCIbreBXa+yO72DSspKQk7NixA++88w7mz58PNzc3tGjRosb5U7am0Agsxr1w4YKgHUZGRloUkBwolUr4+fmhpKQEvr6+9g6HiIiIHNz9+/dRUFCAqKgo82tbdV0aAOgnvQrtP39ab5fWZNZk6vcmJl8TPMLrDIksERERkcOKTdEmtQb78M5zumRXSoIS3sLCQjRq1EjwTi9fvozw8HCzgyIiIiIiA2JTgBYDtN0YSq9qa3YjuwAurvaOTNYE9eHt2LEj/va3v+Ho0aNGtykpKcE///lPxMXF4euvv5YsQCIiIiJ6hIurtvVY62e1/zLZrZGgEd68vDy8++676Nu3Lzw9PdG+fXuEhYXB09MTt2/fRl5eHn766Se0a9cO77//Pvr372/tuImIiIiIBBE0wlu/fn0sWrQIRUVFWLZsGZo1a4YbN27gl19+AQCMGDECOTk5yM7OZrJLRERERLIiqi1Z3bp18eyzz+LZZ5+1VjxERERETkuKlWprE6l+X4JGeGui0Wjw73//m4kwERERkQFubm4AgLt379o5Esei+33pfn/msmhp4YKCAqxevRpr167F9evX0adPH4uCISIiInJGrq6u8Pf3x7Vr1wAAXl5eUCgUdo5KvjQaDe7evYtr167B398frq6WTcwTnfCWl5fjq6++wqpVq3Dw4EGoVCp88MEHGDt2LBdpICIiIjIiJCQEACqTXqqZv79/5e/NEoIT3pycHKxatQqbNm1C06ZNMXLkSGzatAmPPfYYkpKSmOwSERERmaBQKBAaGoqgoCA8ePDA3uHInpubm8UjuzqCE96EhAS88sor+OGHHxATEyPJDyciIiKqbVxdXSVL5EgYwQlv7969sWrVKly7dg0jR45EUlISa0+IiIiISPYEd2n4z3/+g59++gkxMTF48cUXERoaiokTJwIAE18iIiIiki1RbckiIiIwY8YMFBQU4PPPP8f169dRp04dpKam4s0338SxY8esFScRERERkVkUGgs7+t6+fRsbNmzA6tWrcfLkSahUKqlisxulUgk/Pz+UlJRwMh4RERGRDInJ1yxOeB917NgxtGvXTqrd2Q0TXiIiIiJ5E5OvSbLSmo4zJLtERERE5FwkTXiJiIiIiOSGCS8REREROTUmvERERETk1CxKeOfNm4fffvtNolCIiIiIiKRnUcL73nvv4datW1LFQkRERCRrKrUG2fk3sS33MrLzb0KllqzZFVmR4KWFDZGwoxkRERGRrGWeLkL69jwUldyvvC3UzxMzB8UiOS7UjpFRTVjDS0QkAY76EDm3zNNFeHHDMb1kFwCKS+7jxQ3HkHm6qMZ98H3Cfiwa4c3Ly0NYWJhUsRAROSSO+hA5N5Vag/TteTCUnmoAKACkb89D39gQuLooDO6D7xP2ZdEIb0REBFxdXaWKhYjI4Ugx6kNE5rHViOmRglvVjvFHaQAUldzHkQLD85r4PmF/Fo3wEhHVZlKM+hCReWw5YnrtjvFkt6bt+D4hD6zhJSIyk6WjPkRkHluPmAb5eJq9ndD3iR/yb5obHgnAhJeIyEyWjPoQkXlqGjEFtCOmUpY3dIoKRKifJ4yNvyqgHV3uFBVY7T6hx/+4z3/Ekt2/cCKblTDhJSIykyWjPkRkHntcWXF1UWDmoFgAqJb06r6fOSjWYEmC0OP/boUKH+4+h/ZzdrGm1wosSngvXboEtVotVSxERA7FklEfIjKPva6sJMeFYsVf2iHETz+BDfHzxIq/tDNaN1zT+0RVv919gBc4kU1yFk1aa9q0KY4fP46WLVtKFQ8RkcPQjfq8uOEYFIDeJdaaRn2IyDxirqyo1BocKbiFa3fuI8hHe/JpyfGYHBeKvrEhovb56PuEGJzIJi1RI7xlZWVQKpW4fv06li1bhocPH2LlypW4cuUKlEolysrKzApi+fLlaNy4MTw9PZGQkIAjR46Y3H7Lli1o0aIFPD090bp1a+zcubPaNj///DNSUlLg5+cHb29vdOzYEYWFhWbFR0RkjLmjPnLHBvkkV0KvrNwuq0DX+Xsw/J8/YOLmXAz/5w/oOn+PxSOnri4KdG5SH6nx4ejcpL6ghDQ5LhTju0cJHuUFOOFVagqNgPWBV61ahTfeeANKpbLyNhcXF3z00UeYPn263u3e3t6YOnUq3nrrLUEBfPHFF0hLS8PKlSuRkJCAxYsXY8uWLTh79iyCgoKqbX/48GF0794dc+fOxcCBA7Fx40bMnz8fx44dQ1xcHAAgPz8fnTp1wtixYzF8+HD4+vrip59+QmJiosF9VqVUKuHn54eSkhL4+voKeh5EVLtJPZJkT2yQT3Kn69IAGL6yMr57FD7dX1BtYpvuflufjOriFXvauGRYPFLjw60SkzMQk68JSnijo6PRr18/PP3003Bzc4ObmxuaNGmC4OBgKJVK5OXloby8HCqVCocOHcKcOXOgVCrh4eFRY7AJCQno2LEjli1bBgBQq9WIiIjAK6+8gmnTplXb/rnnnkNZWRl27NhReVtiYiLi4+OxcuVKAMCwYcPg5uaGzz//vMafbwgTXiJyRkKScmMfzPZKFIiMMXZi9o8BLTH7u5+NTmxTQHsF5uDUXjY5KVWpNeg6f4/JiXbGbBqXiM5N6lshKucgJl8TVMN78eJFTJs2DREREdXu8/X1RWJiYuX3vXr1wpw5c3D16lU0atTI5H4rKiqQk5OD6dOnV97m4uKCPn36IDs72+BjsrOzMXnyZL3bkpKS8O233wLQJszfffcd/v73vyMpKQnHjx9HVFQUpk+fjsGDBwt5ukRETkfIqC0b5JPcmDpJM1ZPK6aLgy2SyZriMUSXlHPCq3QEJby7d+9GaKjwM/qsrCwEBwfXuN2NGzegUqmqbRscHIwzZ84YfExxcbHB7YuLiwEA165dQ2lpKebNm4c5c+Zg/vz5yMzMxJAhQ7B371706NGj2j7Ly8tRXl5e+f2jJRpERI7O2Kitrkm/btRWbokC1W5CTtJ09bSPklt/bLE/hxNerUPQpLUePXqgTh3hDR26du0qqJzBGnRt0lJTUzFp0iTEx8dj2rRpGDhwYGXJQ1Vz586Fn59f5ZehkWwiInPYe/KXmCb9cksUqPayZCU1ufXHFvtzHH3Cq1xZ1Jbs0qVLAIDHHnvMrMc3aNAArq6uuHr1qt7tV69eRUhIiMHHhISEmNy+QYMGqFOnDmJjY/W2admyJQ4ePGhwn9OnT9crk1AqlUx6ichicpj8JWbUVm6Jgk2pVcCFw0DpVaBeMBDZBXBxtXdUtZKlpTW6Lg7FJfcN7sPW5QI1xQMAvp51MHNgLMICvBx6wquciV54Qq1W45133oGfnx8iIyMRGRkJf39/zJ49W/QiFO7u7mjfvj2ysrL09p+VlYXOnTsbfEznzp31tgeAXbt2VW7v7u6Ojh074uzZs3rbnDt3DpGRkQb36eHhAV9fX70vIiJjhIzaWjJCJSUxo7a1diGNvAxgcRywbiDw9Vjtv4vjtLeTzVm6kpolq6JZg6l4dJT3H+KDXedQcq+Cya6ViB7hfeutt7Bq1SrMmzcPTzzxBADg4MGDmDVrFu7fv493331X1P4mT56MUaNGoUOHDujUqRMWL16MsrIyjBkzBgCQlpaG8PBwzJ07FwAwceJE9OjRAwsXLsSAAQOwefNm/Pjjj/j0008r9zllyhQ899xz6N69O3r27InMzExs374d33//vdinS7WAM7WTIutztMlfYkZtbbKQhtxGUvMygC9HVr9deUV7+58+B2JTbB9XLSZFaY2uP3bVYzXETu31jMXzqKo19SQtQW3JHhUWFoaVK1ciJUX/DWDbtm146aWXcPnyZdFBLFu2DAsWLEBxcTHi4+OxdOlSJCQkAACefPJJNG7cGGvXrq3cfsuWLXj77bfx66+/olmzZnj//ffRv39/vX2uXr0ac+fOxaVLlxATE4P09HSkpqYKiodtyWoPOVxyJschtGVXdv5NDP/nDzXuzxYth3QtkWq6vPtoiyarHRd5GUDmVG0yqeMbBiTPt09SqVYBC5oC90w0968bCEw5z/IGkSwZSJDy+JHbgEbFQzUS5+7GrbIHBu+3dcs0Ryd5H95HeXp64uTJk2jevLne7WfPnkV8fDzu3bsnPmKZYcJbO5hKXjQAJvVphsYNvGXxJkn2V1MvzUc/qHacvIKJm3Nr3KetmsrX1KTf0IiS5IlCXgbwZVqVCB6J4k/rbZ/0/m8fsF7Az0zLAKKrd/ghwyw9YTLnJM1RyOlk2BmIyddE1/C2bdu2cpGIRy1btgxt27YVuzsiuxAyc/3D3b9IuhwlOTYxdYUN6gnrUiN0O0uZs/yxOcunGqVWaUd2TR1xmdO029lSwQFpt7MDe3cBqUqK2nW51eBKiZ1Q7Ed0De/777+PAQMGYPfu3ZUTxbKzs3Hx4kXs3LlT8gCJrEFsI3DWVpGYD6oG3gITWSvnJlVHafdN6YmcC7dtf3n3wmH9MoZqNIDysna7qG7Wj0dH6FOXaV4lp5IslVqDH/53E9O+PiW6dt3Q1QSpa3DlUtpQqzuh2JnohLdHjx44e/YsPv7448rFIYYMGYKXXnoJYWFhkgdIZA1iz565yhSJ+aAS+vq6UVZe80ZmMpUM2aKMQk/p1Zq3EbOdVCK7AlggcDt5EbqYiK1iMTUZS8fQwiU1Je2GVlIT+/4rpxMDubVMq03M6sMbHh4uuhsDkZyYc/bMVaZqNzEfVMbaJVVlrVEcOSVDALTdGKTcTipR3YC6AcC928a3qRto21FnAezVBcTQKOmuvGKDrzVTdCeEQl+nlrzfyu1YsEknFDJIdA3vmjVrsGXLlmq3b9myBevWrZMkKCJrq6nfqCmsraqdxNQV2rOfrZiV1Wwmsou2G4Op34hvuHY7W3JxBQYtNb3NoCWy69BgaZ9ac2SeLkLX+Xsw/J8/VM5teGLeHkzbariEwZQgH0+bvE5leSzAvJp6spzohHfu3Llo0KBBtduDgoLw3nvvSRIUkbUJaQRuDGurai+hH1T2nHRjj2SoRi6u2tZjAIz+RpLn2SexjE3R9tr1qZJk+ITJtgevrSc+GZ2IpryP3+4abq9lyKMnerZ4ncryWPhdclwoDk7thU3jErFkWDw2jUvEwam9mOxakeiShsLCQkRFRVW7PTIyEoWFhZIERWQLQhqBP4q1VQRAcF2hpZNuzJ1kIzTJOXT+um0n8MSmaFuPGezDO8++iWWLAYCnn7YbgwLamt2obrIb2dWx5cQnU6OkYlQ90bNF0i73jgi6TiiWksuEPLkTnfAGBQXh5MmTaNy4sd7tJ06cQP36rGu0FF+4tlU1efn1xl0s3n0OAGuryDihH1TmTrqxZJKN0CRn2d580fu2WGyKNrmU20pr1ZLwf9lvMQwBbDnxSWxHG2OqnujZImm3V0cEsZ/jlnzuy2lCntyJTniHDx+OV199FT4+PujevTsAYN++fZg4cSKGDRsmeYC1CV+49lE1eYkJqSeb5SjJcRj70BI7imPpJJtOUYHw93ITdanZphN4XFzlMwnM2GIYyiLt7fZYDEMAW058snT009/LDcuHt0NilV7Otkja7dERQeznuCWf+3KbkCd3oldaq6iowMiRI7FlyxbUqaPNl9VqNdLS0rBy5Uq4u7tbJVBbssdKa0KXLCXb4Eg7iZF5ugizMn5CsfKPNmMhvh6YldJK1HErZjU3Y69HlVqD9nN2iUp4he7bqahVwOI4E/2BFdpyi9dOyba0wVSyJEU7L0D4ymBVCfnsyjxdhBd+XwHQkJUSfO4JXWVQivd8sZ/jlnzuS/Fe4QzE5GuiR3jd3d3xxRdfYM6cOcjNzUXdunXRunVrREZGmh1wbSdkJum0r0/Bx9MNidEWrnhEgkhVW0XOz9iHdrGyHC9sOCbqQ1vMJBtjr88jBbdEJ7tC9+1U5LoYhgjGSmZ25RVXS4bMvVooZJTU38sNHnVc9E/4bHRVrKZEVUgtvRRXV8W2irO0tZwU7xW1jVl9eAGgWbNmaNasGR4+fIj799mmyRJCaqR+u/cAIz77P5Y4EMmISq3BtK2nTG4zbespwf1QpZhkY+kl6FrTdk+ui2GIVPXkXOrL3ELKJ+YOaS16RFmX8BkjpJew0ETVVC29VL8vsQmopQmr3CfkyZHgtmTbt2/H2rVr9W579913Ua9ePfj7++Opp57C7dsmmneTUWJekGLWIyci6/oh/2aNo6m/3X2AH/JvCtqfFJNsLJ2AU2va7sl1MQwLWKvvrJB2fLrEOzU+HJ2b1Hwl0tKWYUZbpRn5jDQUn5S/L6Gf4/8+XYTs/JsoLrknaHtj+5XjEsUqtQbZ+TexLfcysvNv2ry/cU0Ej/AuWrQIzz77bOX3hw8fxowZM/DOO++gZcuWeOuttzB79mwsWrTIKoE6MzEvSC5xSyQf2f+7IXi7J5pV719elRSTbGrahzG1ru2ebjEMZRGqTVoDUFnDa+vFMESoejlfrdFY7TK3VMv86lgyQinVSnNSlgUI/Rxfn30B67MvINBb2HwnY/uV2xLFjjDpXvAI708//YQuXf448L/66iv07dsXb731FoYMGYKFCxdi+/btVgnSmanUGqjVGvjXdRP8GHs2yyaiP0Yyzl0tFfgIYUmBFAtWmLOoSq1suyfnxTAEMLTy2cv/Mj4B7FHmXuYWO4priiUjlFItKCH091Bccq/GkUtddxShbpdVmLy/ptUY7bm4TVViR9vtRXDCe+fOHb0+uwcPHkTv3r0rv2/VqhWuXDE1AYCq0r1hjVj1f/jtnvhJJqzNIbK9RxON/+YJq+8UM5omxbKjxvYR6ueJv3WPQiiXNNXSLYbhW+V5+4bJtiUZYDzBEPo5IoeyFUuW35aqflXo72H2dz/rnVh0nb/H4iTO1NUXoQmrHJYoluvyzYYILmkIDw/Hzz//jEaNGqG0tBQnTpzAhx9+WHn/zZs34eXlZZUgnZGxQnkx5PCmRVSbmHPcBnhpu6uIIcXlY1P7+HtyS7bdA7StyeoGAH3SgbLrgHdD7RLD9l4MwwRLVj6TU9mKJb2EpapfFVr+c6vKaKyhCW3mdkcBAG8PV5SVqyq/F9PhQupSE7EcqVuE4IR36NCheO211/Dmm29i586dCAkJQWJiYuX9P/74I2JiYqwSpLORYqlGfy83WbxpEdUW5h63c4e0NuvDR4rWeGyvZ4LBFdbCtGUOMk12AfNXPpNj2Yq5y29LVb9qKuk2xVCdsCVXXMvKVZjUpzkaN/AymrCaar9mz+PckbpFCE54Z8yYgcuXL+PVV19FSEgINmzYAFfXP94UNm3ahEGDBlklSGcjxVKNY7pEyeZNi6g2EHvcym3Cho4jTC6xOgddYQ0Qnjj413XTK3GQ62qR5oxQSrnSnLGkO9DbDbfKjI/YVh25tPSK6+ajhUYXiZDzMSvHbhHGCE5469ati/Xr1xu9f+/evZIEVBtIsVTjhF5NJYqGiIQQc9zW93bHvik94V5H8DQJm+BSpNCWMWROheHxvN/H7jKnAS0GyHKkV2jisPzP7eDy+8ij3MtWzBmhNHd02Ni+qibdxcr7mPRFbo2P1b0vmNsdRcfYZX+5H7Ny6xZhitkLT5D5LD3TmWfmJVIiMp+Y4/ZmWQVyLtyWVTmBVK2cHJ6Dr7AmNMFItLCLgklqlfb3U3pV26fYTjXPUtavVk26swX2zr5xpxwqtcbs8ohH/fv3iXC65+AIx6yUo+3WJq/hh1pCyOxUfy83hPhWn2EtxdriRCSe7rgVSg41a4+SqpWTw3PwFdbs3o4qLwNYHAesGwh8PVb77+I47e12IGWrtEfV9DmtM/u7nyu7NhjrmlDPQ9jJwPrsC3pdIBzlmJVDtwghOMJrB0LOiOaZsVQjEVmP7rh9YYOwXqcNvD2sHJE4jjS5xKqcYIU1KS/ni+LAtc9iiRmxrVpeUPWzu31kALq/vwfFynJBP1u3v+efaCxoezkcs/buFiGEQqPR2L85mswolUr4+fmhpKQEvr6+Vvs5ci5EJyLDdp4swoRNx1BTW8kQX0/MSpHPsZydfxPD//lDjdttGpcoq1IMyalV2hFJoyusAVC4AM+sAeIG2zIy0UzN3Jdc5e/NcDmIBgpUeIWgzuTTcK3jPGNphj6nDdGVkpiaePbi7yfLQpIuBYCAGibO6Tj9MWuCmHyNCa8Btkp4ARu/YRGRJHaevIKXNh43uY3uKJbLJT2VWoOu8/fUWPtp7APbqRgbqdSjcKoRS4sVHNCWL9TgJbd3kJL6J0lf8/b+nFSpNVh7qACzv/u5xm1NJZ9Ck+dHBXq743ZZhcljdt+Unsi5cLtW5hFi8jVRp2E7duzAkSNHkJSUhCeeeAJ79uzBBx98ALVajSFDhmD8+PEWBV4bsU8mkePp3yYMK10UmJXxk9HLlHKZVKIj9eQSeychFolNAYauBb4aA2jUxreTcbcGmxNY01yn7Jqk3QPkcCXU1UWBBj7CSpRMlRc8etn/36eLsD77Qo37GxwfhjWHfjV6zKa0DUWPBXt5pVgAwZPWPvnkEzz99NPYuXMn+vfvjw0bNmDw4MEIDw9H48aN8dprr2HJkiXWjJWISDaS40Kx8E/xJreRy6QSHakmlzy6vLKQ5VZVag2y829iW+5lZOfflMUyo/CqbzrZfbRbAwmuab4GfwDSLCdrbAllXY2r7vVmi9eXVP1mdYNc/QQea31jQ4wes+O7R+HT/QU1/n5IS/AI79KlS/Hxxx9j3Lhx2Lt3L/r374+FCxfipZdeAgAkJibi/fffx8SJE60WLBGRnNwoFTYJRQ6TSnQsnVwiti+oHEboDLKwW4NDj3CbI7KLdiU6I7XPag1QjPo4om4hyXKyQltyqdXA7O+s//qSut+smP25uigMToTrsWCvrFuWyY3gEd6CggIkJSUBAHr27AmVSoXu3btX3v/kk0/iwoWah+eJiJyFI60y9ChzWznVlIQA+iN7Qkfo7MKCbg1iR7idgosrVEnzoEH1dFc3oJr+YCTUj6QVlpzoCW3J9dJG27y+pG4HJ3Z/VY/ZnAu37dKyTJZXawQSnPDWr1+/MqG9cuUKHj58iMLCwsr7L1y4gMBA+6+kQURkK0J6aofKZJUhKYjpCyo2ObY53Yilqb+eb7h2u0fIOom3oszTReiaUQ8vVExEkUb/9VyM+njxwWv4j7qT3u2WnOhZkixb6/Uldb9ZS/ZnjzaDjn6iJ7ikITU1FWPHjsWoUaOQkZGBtLQ0vP7663BxcYFCocCUKVPw1FNPWTNWIiJZcaRVhqQg5kNWTHJsl4m7Lq5A8vzfuzUY+eslz9ObsOYIK19Zw6NlLEXohF3lHdDJ5QyC8BuuwR9H1C30RnalWE7W0qsi1np9Sd1v1tz92frqktyXOBZCcMI7f/58VFRUYPPmzejSpQs++ugjLF26FKmpqXjw4AF69OiBuXPnWjNWIiLZsdsiAHYg5kN2d16xoG3tWt8cm6JtPZY5Vb+/rG+YNtmt0pJM9km8FRhK8tVwwQ/qWIPbS3WiV1ONq1BCX19iarKl7q5kzv6krik2xVlO9AQnvN7e3vj000/1bnvjjTcwYcIEPHjwAD4+PpIHR0TkCJLjQtGrRTA+z/4VF27dRWSgF0Z2bgz3Os61ervQD9n2kQF4eWOOoH3avb45NkXbeuzCYe0EtXrB2jIGA63IauNqdVWTfBeoTY7uSnWiV9PVE6FJsJDXl2wnVppQ08qPGkh3dclZTvQsfjf29PS0ONldvnw5GjduDE9PTyQkJODIkSMmt9+yZQtatGgBT09PtG7dGjt37tS7f/To0VAoFHpfycnJFsVIRGRM5uki9FiwF7O/+xnrsy9g9nc/o8eCvZLVtslloojQiTY5F24LWiGqvre7POqbXVyBqG5A62e1/xrpu+uokxQt8WjynuRyBAc9XsVm9zlY6r4Mm93n4KDHq0hyOYK0zpHYNC4RB6f2kixJNFbjGuzrgdd6N4V/XTejjxVaP2/rmmy5HMtiOMuJnmTr/128eBEzZ87E6tWrRT3uiy++wOTJk7Fy5UokJCRg8eLFSEpKwtmzZxEUFFRt+8OHD2P48OGYO3cuBg4ciI0bN2Lw4ME4duwY4uLiKrdLTk7GmjVrKr/38JDXuvZE5BysXdsmt9EnISUc23IvC9pXanyYrC+BVmXLy8hyoUvek1yOYIXb4mr3h+AWVrgtxi++zRHTJK7a/ZaqWuP664272HSkEIuzzht9jNCyCltfqpfyWNbFboyUsTvLiZ5k19tu3bqFdevWiX7cokWLMG7cOIwZMwaxsbFYuXIlvLy8jCbOS5YsQXJyMqZMmYKWLVti9uzZaNeuHZYtW6a3nYeHB0JCQiq/AgICzHpeRETGWLsTgVw7AiTHheLg1F7YNC4RS4bFVxvZE/rB1zc2xJphSk7q1lSOoFNUIMJ93TDTbT0AoOpTc1EAUADNj78LqFVWiUFX4+pRxwWLd59DsdL0SKLQrgliLtVbSupj2ZaxO0s3GsEjvBkZGSbv/9///if6h1dUVCAnJwfTp0+vvM3FxQV9+vRBdna2wcdkZ2dj8uTJerclJSXh22+/1bvt+++/R1BQEAICAtCrVy/MmTMH9esbri0pLy9HefkfDeSVSqXo50JEtY81a9vkPlHE1EQbIROOZPEBqVYJqt19VG2apAho/84fJt5F2H7jiZML8MeqdFHdrBKHqeNBx9/LDcuHt0OiwN7StrpUb41j2ZZlBs7SjUZwwjt48GAoFApoNMZfbgqFuCd748YNqFQqBAfrN/YODg7GmTNnDD6muLjY4PbFxX/MCE5OTsaQIUMQFRWF/Px8vPnmm+jXrx+ys7Ph6lr9zWzu3LlIT08XFTsRkTU/dBx5oohDfEDmZRjpzjC/WneGqqRuTSV3nRo+FLah0NXrzFDT8QAAv919ABcXheC/g60u1VvjWLZ1mYEznOgJTnhDQ0Px8ccfIzU11eD9ubm5aN++vWSBWWLYsGGV/2/dujXatGmDJk2a4Pvvv0fv3r2rbT99+nS9UWOlUomIiAibxEpEjsuaHzqOPlFE1h+QeRm/99+tMoCjLNLe/qf1NSa9UremkjULVqWTijWOB1vVZDty7I9y9BM9wQlv+/btkZOTYzThrWn015AGDRrA1dUVV6/qnxVevXoVISGGa7tCQkJEbQ8A0dHRaNCgAc6fP28w4fXw8OCkNqqVxPSepOqs+aHjDBNFZPkBqVZpR3ZNXWDOnKZtVVZDeUOtoVuVTlkEw783hfb+KqvSSckax4OtrkQ4cuyGfq6jnugJnrQ2ZcoUdOli/MXctGlT7N27V9QPd3d3R/v27ZGVlVV5m1qtRlZWFjp37mzwMZ07d9bbHgB27dpldHsAuHTpEm7evInQUPkPuRPZiqMvE2lPutZCO05ewbCOEZV1eI+y9ENHbhNFzG2npPuATI0PR2eBtZVWdeGwfhlDNZo/6lFJS7cqHQCjr/Qqq9JJzVrHg9TLBRvSKSoQ/l7GW6gBQICXmyxjdyaCR3i7dTNdiO7t7Y0ePXqIDmDy5MkYNWoUOnTogE6dOmHx4sUoKyvDmDFjAABpaWkIDw+vXMVt4sSJ6NGjBxYuXIgBAwZg8+bN+PHHHysXxSgtLUV6ejqeeeYZhISEID8/H3//+9/RtGlTJCUliY6PyBk5wzKR9mKotZDuw+y3u3/0nrX00r2c6mDl1hrNIkLrTK1Yj+qQRK5KJzVrHg9yuBJhbjdeOcTuKCTrw2uu5557DtevX8eMGTNQXFyM+Ph4ZGZmVk5MKywshIvLHwPRXbp0wcaNG/H222/jzTffRLNmzfDtt99W9uB1dXXFyZMnsW7dOvz2228ICwvDU089hdmzZ7NsgQjyn/0vZ8ZOFEp+T3Qn9WmOxg28JPvQSY4LxfI/P463t53WW8jBlnWwTndyJIN6VLFkU3okYlU6a7BmXbg1L9UfKbildzJsyG93H5g9AdWRywxsSaERUHi7f/9+dOnSBXXqCMuPDx06hA4dOjhsgqlUKuHn54eSkhL4+vraOxwiSWXn38Twf/5Q43abxiXyTfQRKrUGXefvMTrbWleve3BqL8mSEUMjq4He7piTGof+bayfZNrjOVudWgUsjqu5HvW1U7Ko4XWq0XWJyOYEQKBtuZcxcXNujdstGRaP1Phw6wfkRMTka4JqePv27avX9qsmvXv3rjaxjIjkwdFn/9uLLRu9A8Yb1d8uq8DLG22z6IStn7NNyKAeVShZLjyiVgEFB4BTX2n/tdJiE6bIri68Bs4wAdUZCBqyjYiIwIwZMzBw4EDUqVMHbm5uiImJQXR0NK5fv47jx4/j/v37UKvV+OGHH6DRaKr1yiWyJ0cbEbAmvvmax5YnCnIpO3HakyM716MKIZfXgB4LehfXZu0jA+CiAEzN83RRaLcj6xGU8Kanp2PatGlYt25dZesxhUKB+fPnY9asWbh79y4A7SppISEhmDt3rsOWM5Dz4SVBffbo3+gMbHmiIJdFJ5z65OjRetQ7RUDZdcC7IVA3QDtqaecRXrm8BipJ0Lu4tqg6wKJWa0wmu4A2Gc65cJtlZFYkKOEdMWIERowYUfl9eXk5/vWvf2H8+PF4/fXX8Y9//AP16tWzWpBE5nK6CTcSkNPsf0diyxMFuYysOv3JkYsrcO82sHum7EYt5fIaAMDexSIY7OJS13RLMh2Hu1LiYAT34X2Uh4cHnn/+ebi6umL06NFMdkmWarokCGgvCQrtJ+pM2L9RPN2JAiB9z92q5DKyapXnLIMa0Eq6UcuqfXl1o5Z5GfaJC8L/tg3qeZjVH1kU9i4WxFjN9W/3THdo0HHIKyUOxKK2ZOfPn0d4OGcUkjzJ7pKgzLB/o3i2Wi5XTiOrkj5nOdWAynzUUshrwN/LDa9/mYtiZXnl7VYp17Jh72JHnW9haoClJg5/pcRBWJTwRkRESBUHkeRkdUlQpti/UTxbnCiYKjvB798P69hIsp9XE0mes9xqQMWMWkaZXnjJGmoqPdIAuG2gt6tVyrVs1LvYkedb1DTAYozuCBrWsRF2nLziUEm+ozGrpIHIEcjlsjA5H1u0RTJWdqLz4e5zNl0K2qLnXONoKrSjqbYsb3CAFddMlR4ZW6rWKuVakV20I/GmFvb1DdduZyZZtmATQejASdV6Xn8vN/h5ueHD3ee4xLuViU54b968iZdffhmxsbFo0KABAgMD9b6I5MJaa68T2UpyXCgOTu2FSX2aGbzfUZIBWdaAOsiKa7rXwKZxiVgyLB6bxiXig2fbmly5S/L+yFbuXewM8y2EDpwsH9Gu8m85qU8z3L77oNrf0mGOawcjuqRh5MiROH/+PMaOHYvg4GAoFBx2J3liNwJyFpuPXjR4u8MsBS3H0VTdqGVNK65ZMGoplaqlR9tyLwt6nKTlWlbsXewM8y2E1t0nRmuvjuhWMTTEYY5rByM64T1w4AAOHjyItm3bWiMeIknZapIRkbU4QzIgy9FU3ajll2mAsVNimay4VpXdyrUe7V1celX794rsYvHvSK7zLcRMoBM7wOIUx7WDEZ3wtmjRAvfu3bNGLERWwW4E5MjkmgyIItfRVAdYcc0Qu3bxcHGVfBKfHOdbmDOBTswAi1Mc1w5GdML78ccfY9q0aZgxYwbi4uLg5qZfgO3r6ytZcERSYTcCclRyTAZEk/NoqpVGLa3J2cq15NSGD7BswSKhAyxOcVw7GNGT1vz9/aFUKtGrVy8EBQUhICAAAQEB8Pf3R0AA14EmIpKS00y+1I2m+lZJFHzD7L8srW7UsvWz2n9lnOzqONPiMbZc1KUmUkygE9LRxGmOawcieoR3xIgRcHNzw8aNGzlpjYjIypxqNM8BR1PlzJnKteQy38JWtbVOdVw7CIVGoxHV58PLywvHjx9HTEyMtWKyO6VSCT8/P5SUlLBEg4hkwZGb8hMJZe+V1rblXsbEzbk1brdkWDxS4y1faZbHtWXE5GuiR3g7dOiAixcvOnXCS0QkN840mkdkjL3nW9i6tpbHte2ITnhfeeUVTJw4EVOmTEHr1q2rTVpr06aNZMEREdEf7J0MEDk7e0yg43FtG6JLGlxcqs9zUygU0Gg0UCgUUKlsuDSklbCkgYiIqHbSdWkADNfWOtqkQGdm1ZKGgoICswMjIiIikjO5TKAjaYlKeB88eIBevXphx44daNmypbViIiIiIrIb1tY6H1EJr5ubG+7fd/5VP3RVHkql0s6REBERkb20auiGVg21c5XKSu/YORqqSpenCanOFV3D+9577+HcuXP47LPPUKeO6IoIh3Dp0iVERETYOwwiIiIiqsHFixfx2GOPmdxGdML79NNPIysrC/Xq1UPr1q3h7e2td//WrVvFRyozarUaV65cgY+Pj+QLayiVSkRERODixYucEOdg+LdzbPz7OTb+/RwX/3aOTc5/P41Ggzt37iAsLMxgU4VHiR6i9ff3xzPPPGN2cI7AxcWlxjMFS/n6+sruhUPC8G/n2Pj3c2z8+zku/u0cm1z/fn5+foK2E53wrlmzRnQwRERERET2Ynr8l4iIiIjIwYlOeK9evYqRI0ciLCwMderUgaurq94Xmebh4YGZM2fCw8PD3qGQSPzbOTb+/Rwb/36Oi387x+Ysfz/Rk9b69euHwsJCTJgwAaGhodUmdaWmpkoaIBERERGRJUQnvD4+Pjhw4ADi4+OtFBIRERERkXRElzREREQIavBLRERERCQHohPexYsXY9q0afj111+tEA4RERERkbQElTQEBATo1eqWlZXh4cOH8PLygpubm962t27dkj5KIiIiIiIzCerDu3jxYiuHUTssX74cCxYsQHFxMdq2bYuPPvoInTp1sndYVIP9+/djwYIFyMnJQVFREb755hsMHjzY3mGRQHPnzsXWrVtx5swZ1K1bF126dMH8+fMRExNj79CoBitWrMCKFSsqryi2atUKM2bMQL9+/ewbGJll3rx5mD59OiZOnMi8wgHMmjUL6enperfFxMTgzJkzdorIMoIS3lGjRlk7Dqf3xRdfYPLkyVi5ciUSEhKwePFiJCUl4ezZswgKCrJ3eGRCWVkZ2rZti+effx5Dhgyxdzgk0r59+/Dyyy+jY8eOePjwId5880089dRTyMvLq7Y0OsnLY489hnnz5qFZs2bQaDRYt24dUlNTcfz4cbRq1cre4ZEIR48exSeffII2bdrYOxQSoVWrVti9e3fl93XqiF6vTDZEd2lwdXVFUVFRtSTt5s2bCAoKgkqlkjRAZ5GQkICOHTti2bJlAAC1Wo2IiAi88sormDZtmp2jI6EUCgVHeB3c9evXERQUhH379qF79+72DodECgwMxIIFCzB27Fh7h0IClZaWol27dvj4448xZ84cxMfHc4TXAcyaNQvffvstcnNz7R2KJERPWjOWH5eXl8Pd3d3igJxRRUUFcnJy0KdPn8rbXFxc0KdPH2RnZ9sxMqLap6SkBIA2cSLHoVKpsHnzZpSVlaFz5872DodEePnllzFgwAC9z0ByDL/88gvCwsIQHR2NESNGoLCw0N4hmU3w2PTSpUsBaEe4PvvsM9SrV6/yPpVKhf3796NFixbSR+gEbty4AZVKheDgYL3bg4ODHbYWhsgRqdVqvPbaa3jiiScQFxdn73BIgFOnTqFz5864f/8+6tWrh2+++QaxsbH2DosE2rx5M44dO4ajR4/aOxQSKSEhAWvXrkVMTAyKioqQnp6Obt264fTp0/Dx8bF3eKIJTng//PBDANoR3pUrV+otI+zu7o7GjRtj5cqV0kdIRCSRl19+GadPn8bBgwftHQoJFBMTg9zcXJSUlOCrr77CqFGjsG/fPia9DuDixYuYOHEidu3aBU9PT3uHQyI9Ojm0TZs2SEhIQGRkJL788kuHLCkSnPAWFBQAAHr27ImtW7ciICDAakE5mwYNGsDV1RVXr17Vu/3q1asICQmxU1REtcuECROwY8cO7N+/H4899pi9wyGB3N3d0bRpUwBA+/btcfToUSxZsgSffPKJnSOjmuTk5ODatWto165d5W26K8LLli1DeXm53uAZyZu/vz+aN2+O8+fP2zsUs4iu4d27dy+TXZHc3d3Rvn17ZGVlVd6mVquRlZXFWjQiK9NoNJgwYQK++eYb7NmzB1FRUfYOiSygVqtRXl5u7zBIgN69e+PUqVPIzc2t/OrQoQNGjBiB3NxcJrsOprS0FPn5+QgNDbV3KGYxq7/EpUuXkJGRgcLCQlRUVOjdt2jRIkkCczaTJ0/GqFGj0KFDB3Tq1AmLFy9GWVkZxowZY+/QqAalpaV6Z7QFBQXIzc1FYGAgGjVqZMfISIiXX34ZGzduxLZt2+Dj44Pi4mIAgJ+fH+rWrWvn6MiU6dOno1+/fmjUqBHu3LmDjRs34vvvv8d//vMfe4dGAvj4+FSrlff29kb9+vVZQ+8A3njjDQwaNAiRkZG4cuUKZs6cCVdXVwwfPtzeoZlFdMKblZWFlJQUREdH48yZM4iLi8Ovv/4KjUajd9mC9D333HO4fv06ZsyYgeLiYsTHxyMzM7PaRDaSnx9//BE9e/as/H7y5MkAtP2p165da6eoSKgVK1YAAJ588km929esWYPRo0fbPiAS7Nq1a0hLS0NRURH8/PzQpk0b/Oc//0Hfvn3tHRqR07t06RKGDx+OmzdvomHDhujatSt++OEHNGzY0N6hmUV0H95OnTqhX79+SE9Ph4+PD06cOIGgoCCMGDECycnJePHFF60VKxERERGRaKITXh8fH+Tm5qJJkyYICAjAwYMH0apVK5w4cQKpqamVS0ASEREREcmB6Elr3t7elXW7oaGhyM/Pr7zvxo0b0kVGRERERCQB0TW8iYmJOHjwIFq2bIn+/fvj9ddfx6lTp7B161YkJiZaI0YiIiIiIrOJLmn43//+h9LSUrRp0wZlZWV4/fXXcfjwYTRr1gyLFi1CZGSktWIlIiIiIhJNdMJLRERERORIzOrDq1NaWgq1Wq13m6+vr0UBERERERFJSfSktYKCAgwYMADe3t7w8/NDQEAAAgIC4O/vzxXYiIiIiEh2RI/w/uUvf4FGo8Hq1asRHBwMhUJhjbiIiEiktWvXVq7eOHHiRCxevFjS/T/55JPYt28fAOD48eOIj4+XdP9ERNYiOuE9ceIEcnJyEBMTY414iIjIAr6+vjh79iy8vb0Fbf/111/jT3/6EwoLCxEeHl7t/mbNmmHQoEFYtGgRtm7divz8fHTq1EnqsImIrEp0SUPHjh1x8eJFa8RCREQWUigUCAkJgY+Pj6DtU1JSUL9+faxbt67affv378f58+cxduxYAEBgYKDDLitKRLWb6IT3s88+w/z587Fu3Trk5OTg5MmTel9ERGS5MWPGoGHDhnj//fcrbzt16hRcXV2xfft2UfsqLy/HG2+8gfDwcHh7eyMhIQHff/89AMDNzQ0jR47E2rVrqz1u9erVSEhIQKtWrSx5KkREdie6pOH69evIz8+vrBMDtCMKGo0GCoUCKpVK0gCJiGqjRYsWoUuXLnjxxRfxzDPPoEmTJvjHP/6Bjh07YtCgQaL2NWHCBOTl5WHz5s0ICwvDN998g+TkZJw6dQrNmjXD2LFjsWjRIuzfvx/du3cHoO3C89VXX+HDDz+0xtMjIrIp0Qnv888/j8cffxybNm3ipDUiIisJCAjAuHHj8MUXX2DDhg3o378/tm3bht27d4vaT2FhIdasWYPCwkKEhYUBAN544w1kZmZizZo1eO+99xAbG4vExESsXr26MuH98ssvodFoMGzYMMmfGxGRrYlOeC9cuICMjAw0bdrUGvEQEdEj0tLSMHv2bBw6dAg9e/ZE7969RT3+1KlTUKlUaN68ud7t5eXlqF+/fuX3zz//PCZNmoSPPvoIPj4+WL16NYYOHSq4FpiISM5EJ7y9evXCiRMnmPASEdnAkCFD8NJLL+H8+fM4fPiw6MeXlpbC1dUVOTk5cHV11buvXr16lf8fNmwYJk2ahC+//BLdu3fHoUOHMHfuXIvjJyKSA9EJ76BBgzBp0iScOnUKrVu3hpubm979KSkpkgVHRFTb1atXD02aNIGPjw86d+4s+vGPP/44VCoVrl27hm7duhndzsfHB0OHDsXq1auRn5+P5s2bm9yeiMiRiE54X3jhBQDAO++8U+0+TlojIpJWZmYmTp48icDAQFRUVMDd3V3U45s3b44RI0YgLS0NCxcuxOOPP47r168jKysLbdq0wYABAyq3HTt2LLp164aff/4ZU6dOlfqpEBHZjei2ZGq12ugXk10iIuloNBq8/fbbGDp0KNzd3fHdd9+ZtZ81a9YgLS0Nr7/+OmJiYjB48GAcPXoUjRo10tuua9euiImJgVKpRFpamhRPgYhIFkSP8Jpy6tQpZGRkICoqCn/+85+l3DURUa2zdetW5ObmYuPGjfjkk0/w+eef4+mnnxa9Hzc3N6SnpyM9Pb3Gbc+cOWNOqEREsiZ6hPdRKpUKe/bswWuvvYaoqCh0794deXl51WYDExGROGq1GjNmzEBaWhqaN2+OtLQ0ZGRkYN26dbh9+7bRx5WUlKBevXpWKUno168fF6EgIoek0Gg0GjEPUCqV+Pe//41t27YhMzMTfn5+GDRoEFJTU9GjRw/UqSPpoDERUa30+eef469//SvOnj2Lxo0bAwD+8Y9/YPny5Rg1apTBBSHu3LmDq1evAgD8/f3RoEEDSWO6fPky7t27BwBo1KiR6HpiIiJ7EZzwfvTRR8jIyMCBAwfQunVrpKSkICUlBW3btrV2jEREREREZhOc8CYnJyM1NRUpKSkIDw+vcftDhw6hQ4cO8PDwsDhIIiIiIiJziS5pEMrX1xe5ubmIjo62xu6JiIiIiASxaNKaKVbKo4mIiIiIRLFawktEREREJAdMeImIiIjIqTHhJSIiIiKnZrWEV6FQWGvXRERERESCcdIaERERETk1q7UlIyIiIiKSg/8HCxgtnxsGq14AAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure(figsize=(8, 2))\n", - "plt.scatter(lines.ExPot[neutrals], (lines.A_18Sco - lines.A_sun)[neutrals], label=\"Fe I\")\n", - "plt.scatter(lines.ExPot[~neutrals], (lines.A_18Sco - lines.A_sun)[~neutrals], label=\"Fe II\")\n", - "plt.ylabel(\"A(Fe)_\\mathrm{18 Sco} - A(Fe)_\\mathrm{sun}\")\n", - "plt.xlabel(\"χ [eV]\")\n", - "plt.legend()\n", - ";" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.11.5" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/misc/Tutorial notebooks/EW fitting/EW fitting.ipynb b/misc/Tutorial notebooks/EW fitting/EW fitting.ipynb deleted file mode 100644 index 03b626a9..00000000 --- a/misc/Tutorial notebooks/EW fitting/EW fitting.ipynb +++ /dev/null @@ -1,220 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "b1c9787e-781e-439f-9599-ea5c0d2f644b", - "metadata": {}, - "source": [ - "# Equivalent width fitting\n", - "This notebook is an example of how to go from measured equivalent widths to abundances. We'll use the data from [Melendez 2014](https://ui.adsabs.harvard.edu/abs/2014ApJ...791...14M/abstract) and redo their differential analysis of iron lines 18 Sco." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "25fa263f-d156-4a4a-ab1e-02bea37f9d9b", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mPrecompiling Korg [acafc109-a718-429c-b0e5-afd7f8c7ae46]\n" - ] - } - ], - "source": [ - "using Revise, Korg, DataFrames, CSV, PyPlot" - ] - }, - { - "cell_type": "markdown", - "id": "988c6101-cff0-4320-a444-954379815fdc", - "metadata": {}, - "source": [ - "## parsing the linelist" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "eeac654b-0dcd-420c-82b1-bdee393e2f44", - "metadata": {}, - "outputs": [], - "source": [ - "# use the CSV library to read the contents of Table 1, which were downloaded from the journal website\n", - "# the \"lines\" objects is a DataFrame, which we'll use to organize our data\n", - "lines = CSV.read(\"Table 1.dat\", DataFrame; skipto=25, delim=' ', ignorerepeated=true,\n", - " header=[\"wl\", \"species\", \"ExPot\", \"log_gf\", \"C6\", \"EW_18Sco\", \"EW_Sun\"])\n", - "\n", - "# convert the numbers in the species column to Korg.Species objects\n", - "lines.species = Korg.Species.(lines.species)\n", - "\n", - "# let's look at Fe lines only\n", - "filter!(lines) do row\n", - " Korg.get_atoms(row.species) == [26]\n", - "end\n", - "\n", - "#sort by wavelength\n", - "sort!(lines, :wl)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "87d2b027-79bc-4b86-9480-82da48dcae3a", - "metadata": {}, - "outputs": [], - "source": [ - "# to pass the lines to Korg, we need to turn each row of the \"lines\" DataFrame into a Korg.Line object\n", - "linelist = Korg.Line.(lines.wl, # can be in either cm or Å (like these), but NOT nm\n", - " lines.log_gf,\n", - " lines.species, # needs to be a Korg.Species, which we handled in the cell above\n", - " lines.ExPot) # excitation potential, i.e. lower level energy (must be in eV)\n", - ";" - ] - }, - { - "cell_type": "markdown", - "id": "9e0db71d-1010-4a39-918c-bd1ce71c2069", - "metadata": {}, - "source": [ - "## Set up params for each star" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "69f5b674-9538-4fd8-b95e-ae6e17a0f9d8", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mloading the model atmosphere grid into memory. This will take a few seconds, but will only happen once per julia session.\n" - ] - }, - { - "data": { - "text/plain": [ - "Korg.PlanarAtmosphere{Float64, Float64, Float64, Float64, Float64} with 56 layers" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# solar params\n", - "sun_Teff, sun_logg, sun_Fe_H, sun_vmic = 5777, 4.44, 0.0, 1.0\n", - "\n", - "# vector of abundances for the sun\n", - "sun_A_X = Korg.format_A_X(sun_Fe_H)\n", - "\n", - "# interpolate a model atmosphere for the sun\n", - "sun_atm = Korg.interpolate_marcs(sun_Teff, sun_logg, sun_A_X) \n", - "\n", - "# and likewise for 18 Sco\n", - "sco_teff, sco_logg, sco_fe_h, sco_vmic = (5823, 4.45, 0.054, sun_vmic + 0.02)\n", - "sco_A_X = Korg.format_A_X(sco_fe_h)\n", - "sco_atm = Korg.interpolate_marcs(sco_teff, sco_logg, sco_A_X)" - ] - }, - { - "cell_type": "markdown", - "id": "e1e6a43b-2d28-4d03-b5d3-3a1e0ef6d135", - "metadata": {}, - "source": [ - "## Calculate the abundances for each star\n", - "Because Julia uses \"just in time\" compilation, the first call to `Korg.Fit.ews_to_abundances` will take several seconds. The second is much faster, because the code is already compiled. If you rerun the cell below, both calculations will be fast, even though no data is cached." - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "fc38fcb2-2bae-4fe8-af6f-c3604cc05832", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " 18.033960 seconds (103.51 M allocations: 7.324 GiB, 5.10% gc time, 84.61% compilation time)\n", - " 2.616027 seconds (15.51 M allocations: 2.969 GiB, 8.50% gc time)\n" - ] - } - ], - "source": [ - "# calculate abundances from the EWs for each star\n", - "@time lines.A_sun = Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, lines.EW_Sun, vmic=sun_vmic)\n", - "@time lines.A_18Sco = Korg.Fit.ews_to_abundances(sco_atm, linelist, sco_A_X, lines.EW_18Sco, vmic=sco_vmic)\n", - ";" - ] - }, - { - "cell_type": "markdown", - "id": "31560282-f1b1-4538-a3c0-cc1961b52ecb", - "metadata": {}, - "source": [ - "# Plot the results\n", - "Plotting as a function of wavelength and equivalenth width is left as an excercise for the reader." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "e4acef29-6ad0-4ae3-8263-e72d5e2f4506", - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAADZCAYAAAApK1iAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDuUlEQVR4nO3deVxTV/o/8E+CLJXNBRWx4ooLoiIKiNW2KgpdtNQZq7YVa60dO9pRaa3Lr4qMdax20xmt2k5dWku1q7ba0lr81hUVoVAVq5ZRcRRwBwRBTPL7I5MIkpB7k5vcG/J5v155KTcnlwNJyHPPec5zVDqdTgciIiIiIheglrsDRERERESOwuCXiIiIiFwGg18iIiIichkMfomIiIjIZTD4JSIiIiKXweCXiIiIiFwGg18iIiIichkMfomIiIjIZTSSuwNKp9VqcfHiRfj6+kKlUsndHSIiIiK6h06nQ1lZGYKCgqBW1z+2y+DXgosXL6Jt27Zyd4OIiIiILDh//jzuv//+etsw+LXA19cXgP6X6efnJ3NviIiIiOhepaWlaNu2rTFuqw+DXwsMqQ5+fn4MfsnlabQ6HD5zDZfKKtHS1wtRHZrBTc10ICIiUgYhKaoMfolIkLRjhUj5Lg+FJZXGY639vZA8IhTxYa1l7BkREZFwrPZARBalHSvES5uyawW+AFBUUomXNmUj7VihTD0jIiISh8EvEdVLo9Uh5bs86EzcZziW8l0eNFpTLYiIiJSFaQ9EVK/DZ67VGfGtSQegsKQSh89cQ0yn5o7rGBGRk9PpdLhz5w40Go3cXXEK7u7ucHNzs/k8DH6JqF6XyswHvta0IyIi4Pbt2ygsLERFRYXcXXEaKpUK999/P3x8fGw6j2TB7/79+3H27NlaVy+JiYlSnZ6IZNLS10vSdkRErk6r1eLMmTNwc3NDUFAQPDw8uJGWBTqdDpcvX8Z///tfhISE2DQCLEnwO27cOBQVFaFPnz7GzvBJJGoYojo0Q2t/LxSVVJrM+1UBCPTXlz0jIiLLbt++Da1Wi7Zt26Jx48Zyd8dptGjRAmfPnkV1dbX8wW9ubi7y8vKkOBURKYybWoXkEaF4aVM2VECtANhwiZs8IpT1fonIZq5WS9zSNrxUm1QDq5IEv1FRUTh58iS6du0qxemISGHiw1pj9bMRder8BrLOLxFJhLXEyVEkCX5zcnLQu3dvdO3aFZ6entDpdFCpVDh8+LAUpyciBYgPa41hoYEuNSpDRI5hqCV+b2qVoZb46mcjGACTZCQJfrdt2ybFaYhI4dzUKpYzIyJJWaolroK+lviw0EBebJMkJEk2adeunckbERERUX3E1BIn+T333HNQqVR1bn/88YdV51u4cCHCw8Ol7aQFkoz8fvzxxyaPs9QZERER1Ye1xG0jxyLB+Ph4rF+/vtaxFi1a2PV7SkmS4Pfo0aPG/1dVVWHnzp3o1asXg18iIiKqF2uJW0+uRYKenp4IDAw0ed+2bduQkpKCvLw8BAUFYcKECfh//+//oVEj5eyrJklP3nrrrVpf37x5EwkJCVKcmoiIiBow1hK3jhIXCe7duxeJiYn45z//iUGDBiE/Px8vvvgiACA5OdmhfamPXQrMqVQqnDt3zh6nJiIiogbEUEscuFs73IC1xE2ztEgQ0C8S1GhNtbDd9u3b4ePjY7yNHj1a/z1TUjBnzhxMmDABHTt2xLBhw7Bo0SKsXbvWLv2wliQjv5GRkcbCwxqNBoWFhXjttdekODURERE1cKwlLo6YRYL2qNAzePBgrF692vi1t7c3AP2mZ/v378fixYuN92k0GlRWVqKiokIxu9lJEvx++eWXd0/YqBFatmwJd3d3KU5NRERELoC1xIWTe5Ggt7c3OnfuXOf4zZs3kZKSglGjRtW5z8tLOTnbkgS/p06dQkxMDHx8fPD+++8jJycHSUlJ6NatmxSnJyIiIhfAWuLCKHWRYEREBE6ePGkyMFYSSXJ+X331Vfj4+ODgwYP49NNPERsbi0mTJklxaiIiIiKqwbBI0NyYuAr6qg+OXiS4YMECfPzxx0hJScHx48dx4sQJbN68Ga+//rpD+2GJpAvetm7diilTpuCpp55CRUWFlKcmIiIiIih3kWBcXBy2b9+On376CZGRkejfvz/ee+89xW18JknwGxQUhPHjx+Ozzz7D448/jqqqKmg0GqvOtWrVKrRv3x5eXl6Ijo7G4cOH623/xRdfoFu3bvDy8kLPnj3x/fff17rf1C4kKpWqTnk2IiIiImdhWCQY6F87tSHQ38uuZc42bNiArVu3mr0/Li4O+/fvR0VFBUpKSnDo0CFMnjzZbPuFCxciJydH+o7WQ7IFbz/++CMWLFiApk2borCwEG+//bbo82zZsgVJSUlYs2YNoqOjsXz5csTFxeHkyZNo2bJlnfYHDhzAuHHjsGTJEjz++ONITU1FQkICsrOzERYWBgAoLCys9ZgffvgBkyZNwp/+9CfrflgiIiIiBeAiQeuodDqdVUXgqqurUVRUhIqKCrRo0QLNmtmeVxIdHY3IyEisXLkSAKDVatG2bVu8/PLLmDNnTp32Y8aMQXl5ObZv32481r9/f4SHh2PNmjUmv0dCQgLKysqQnp4uqE+lpaXw9/dHSUkJ/Pz8rPipiIiIiO6qrKzEmTNn0KFDB0VVQVC6+n5vYuI1UWkPZWVlWL16NR566CH4+fmhffv26N69O1q0aIF27dph8uTJyMzMFP/TALh9+zaysrIQGxt7t3NqNWJjY5GRkWHyMRkZGbXaA/rhdnPti4uLsWPHjnoX41VVVaG0tLTWjYiIiIgaBsHB77vvvov27dtj/fr1iI2NxdatW5GTk4NTp04hIyMDycnJuHPnDoYPH474+HicPn1aVEeuXLkCjUaDVq1a1TreqlUrFBUVmXxMUVGRqPYbN26Er6+vyfpzBkuWLIG/v7/x1rZtW1E/BxEREREpl+Cc38zMTOzZswc9evQweX9UVBSef/55rFmzBuvXr8fevXsREhIiWUelsG7dOjzzzDP1TjHMnTsXSUlJxq9LS0sZABMRERE1EIKD388++0xQO09PT0yZMkV0RwICAuDm5obi4uJax4uLixEYGGjyMYGBgYLb7927FydPnsSWLVvq7Yenpyc8PT1F9p6IiIiInIGkdX5t4eHhgb59+9ZaiKbVapGeno6YmBiTj4mJiamzcG3nzp0m23/00Ufo27cvevfuLW3HiYiIiMhpWBX8lpSU4MUXX0Tnzp3RvXv3OuXErJWUlIQPP/wQGzduxIkTJ/DSSy+hvLwcEydOBAAkJiZi7ty5xvbTp09HWloa3nnnHfz+++9YuHAhjhw5gmnTptU6b2lpKb744gu88MILkvSTiIiIiJyTVXV+p06divz8fCxbtgzPPvssbt26BQCYOXMmOnXqVCf4FGrMmDG4fPkyFixYgKKiIoSHhyMtLc24qK2goABq9d14fcCAAUhNTcXrr7+OefPmISQkBFu3bjXW+DXYvHkzdDodxo0bZ1W/iIiIiKhhsKrOb/PmzfHzzz+jT58+8PX1RW5uLjp27Ii0tDTMnz/f6nJnSsQ6v0RERCQlZ67z+9xzz2Hjxo11jp8+fRqdO3cWfb6FCxcaK4iZ+romqer8WjXyq9Pp4OvrW+d4SEiI6BJnREREROQ84uPjsX79+lrHWrRoIVNvxLMq5/eRRx7Bp59+Wud4eXk5VCpuqUdERETkEFoNcGYvcPRL/b9ajd2/paenJwIDA2vd3NzcAADbtm1DREQEvLy80LFjR6SkpODOnTt275MYVo38LlmyBP369QOgHwVWqVSorKzEokWLEBERIWkHiYiIiMiEvG+BtNlA6cW7x/yCgPilQOhIh3dn7969SExMxD//+U8MGjQI+fn5ePHFFwEAycnJDu+POVaN/AYHB+PAgQM4cOAAKioqEBUVhSZNmmD37t1YunSp1H0kIiIiopryvgU+T6wd+AJAaaH+eN63dvvW27dvh4+Pj/E2evRoAEBKSgrmzJmDCRMmoGPHjhg2bBgWLVqEtWvX2q0v1hA18ltRUYHGjRsDADp37owff/wRBQUFyM3Nhbu7O6Kjo9G0aVO7dJSIiIiIoE9tSJsNwFTNAh0AFZA2B+j2GKB2k/zbDx48GKtXrzZ+7e3tDQDIzc3F/v37sXjxYuN9Go0GlZWVtWJIuYkKfv39/VFYWIiAgADjseDgYAQHB0veMSIiIiIy4dyBuiO+teiA0gv6dh0GSf7tvb29TVZ2uHnzJlJSUjBq1Kg69ympqoWo4Fej0UCr1Rq/HjRoEL788ktjHV4iIiIisrObxdK2k0hERAROnjxpVckzR7JqwZtBTk4OysvLpeoLkVPRaHU4fOYaLpVVoqWvF6I6NIObmtVOiIjIznwEDjoKbSeRBQsW4PHHH0dwcDD+/Oc/Q61WIzc3F8eOHcMbb7zh0L7Ux6bgl8hVpR0rRMp3eSgsqTQea+3vheQRoYgPay1jz4iIqMFrN0Bf1aG0EKbzflX6+9sNcGi34uLisH37dvz973/H0qVL4e7ujm7duuGFF15waD8sEbXDm1qtxrvvvosHH3wQPXv2RLNmzYy7uzVU3OGN7pV2rBAvbcqu8+fGMOa7+tkIBsBERGSWJDu8Gao9AKgdAP/v0+ipj2Upd2ZPUu3wJqrU2aBBg5CcnIx+/frBx8cHFRUVSE5OxgcffIDMzExUVVWJ/0mInIhGq0PKd3lm19cCQMp3edBoRe8aTkREJFzoSH2A63fPYItfUIMMfKUkKu1h9+7dAPT7N2dlZSE7OxvZ2dmYM2cObty4gUaNGqF79+7Izc21S2eJ5Hb4zLVaqQ730gEoLKnE4TPXENOpueM6RkRErid0pL6c2bkD+sVtPq30qQ52KG/WkFiV8xsSEoKQkBCMHTvWeOzMmTM4cuQIfv31V8k6R6Q0l8rMB77WtCMiIrKJ2s0u5cwaMsHBb0FBQb31fDt06IAOHToYd/m4cOEC2rRpY3sPiRSkpa+w3Cyh7YiIiMixBOf8RkZG4i9/+QsyMzPNtikpKcGHH36IsLAwfPXVV5J0kEhJojo0Q2t/L5graKaCvupDVIdmjuwWERERCSR45DcvLw+LFy/GsGHD4OXlhb59+yIoKAheXl64fv068vLycPz4cURERGDZsmV49NFH7dlvIlm4qVVIHhGKlzZlQwWT62uRPCKU9X6JiMgiEQW3CNL9vkSVOgOAW7duYceOHdi3bx/OnTuHW7duISAgAH369EFcXBzCwsIk6ZhSsNQZmcI6v0REZC2NRoNTp06hZcuWaN6ci6OFKikpwcWLF9G5c2e4u7vXuk9MvCY6+HU1DH7JHO7wRkRE1iosLMSNGzfQsmVLNG7cGCoVPz/qo9VqcfHiRbi7uyM4OLjO70tMvMYd3ois5KZWsZwZEZELs2UQJDAwEABw6dIle3axQVGr1SYDX7EkC35LS0uRk5ODnJwc/O1vf5PqtERERESKY2v6m0qlQuvWrdGyZUtUV1fbs6sNhoeHB9RqUfuzmWRV2kNBQYEx0DXczp07B51OB29vb5SVldncMaVg2gMRERHVxG3ulcduaQ9DhgxBbm4url+/Dn9/f4SGhiIsLAwFBQX46KOPMHToULRt29amzhMROTvmgxM1XJa2uVdBv839sNDAet/3/DshH1HB7759+zBr1iz89a9/rbWBxbp16xAVFcXAl4hcHiuBEDVsUmxzz78T8hKVOHHo0CHs3bsXU6dOxalTp+zVJyIip2SYCr33g7GopBIvbcpG2rFCmXpG1PBptDpk5F/FtpwLyMi/Co3WPsWsbN3mnn8n5Cdq5LdPnz7Ys2cPUlNTERcXh0cffRTJycn26hsRkdOQaiqUiMRz5EiqLdvc8++EMli1ZO7pp5/G8ePH0bRpU/To0QNarRYajUbqvhEROQ0xU6FEJB1Hj6Tass290L8Tr3yeg/2nr9ht9NrVWV0vonHjxnjjjTdw6NAhPP744xg6dCjefvtt3Lp1S8r+ERE5BVunQolIPEsjqYB+JFXKINKwzT2AOgGwpW3uhb7/t+ZcxDMfHULfN3YyDcIObC6W1rFjR2zbtg2bNm3C+vXr0bFjRyn6RUTkVGyZCiUi68g14xIf1hqrn41AoH/t93Ogv1e9Zc7Evv9vVFRjCvOAJSfZJhfDhw9Hbm4uVq5cKdUpiYichmEqtKik0uQolAr6D0ZTU6FEZB0xMy5SlxaLD2uNYaGBos5p6e+EOQu/Pc48YAlJur1xo0aNMGPGDClPSUTkFAxToS9tyoYKqPXBZmkqlIisI3Qk9eyVCgxcukvyBXFit7k3/J2Ysilb1PcpKq2qt3QaiWP7HnESW7VqFdq3bw8vLy9ER0fj8OHD9bb/4osv0K1bN3h5eaFnz574/vvv67Q5ceIERo4cCX9/f3h7eyMyMhIFBQX2+hGIyEVZOxVKRNYRsvisSWN3LP/5lNOXFuN6AenYHPyeOnUKd+7ckaIv2LJlC5KSkpCcnIzs7Gz07t0bcXFxuHTpksn2Bw4cwLhx4zBp0iT8+uuvSEhIQEJCAo4dO2Zsk5+fj4EDB6Jbt2745Zdf8Ntvv2H+/Pnw8mLeHRFJLz6sNfbNHoLPJvfHirHh+Gxyf+ybPURxga+YmqiOqp9KJJalxWeGV6ojF8TVx7BAzxpcLyAdlU6ns+kZd3Nzw4kTJ9ClSxebOxMdHY3IyEhj3rBWq0Xbtm3x8ssvY86cOXXajxkzBuXl5di+fbvxWP/+/REeHo41a9YAAMaOHQt3d3d88sknVvVJzF7RRETOQExNVO5ERUpRX86uudfp2Mi2eO/n0xbP/dnk/g5JKcjIv4pxHx4U/bhAP0/snzOUaVP1EBOv2Zzza2PsbHT79m1kZWVh7ty5xmNqtRqxsbHIyMgw+ZiMjAwkJSXVOhYXF4etW7cC0AfPO3bswGuvvYa4uDj8+uuv6NChA+bOnYuEhAST56yqqkJVVZXx69LSUtt+MCIiBTHURL33L7dhCrhmeoaYtkT2ZOkizNzis+2/XRR0fkelFFj7fRaO7MHAV0KKyfm9cuUKNBoNWrVqVet4q1atUFRUZPIxRUVF9ba/dOkSbt68iTfffBPx8fH46aef8OSTT2LUqFHYvXu3yXMuWbIE/v7+xlvbtm0l+OmIiORPHxBTE1WO+qmKotUAZ/YCR7/U/6vlRk5yEbqJhWHx2RPhbRDTqTnc1CrFlSAU+32aNHbHGl5kSk7Sag9Ko9VqAQBPPPEEZs6cCQAIDw/HgQMHsGbNGjz00EN1HjN37txao8mlpaUMgInILKHlk5SQPiC2JqrQtg1uBXret0DabKC0xqihXxAQvxQIHSlfv1yQrdsBK60EoZBSZ56N1Ijv0Qp/jmiLASEBHPG1A8WM/AYEBMDNzQ3FxcW1jhcXFyMwMNDkYwIDA+ttHxAQgEaNGiE0NLRWm+7du5ut9uDp6Qk/P79aN3Idco/MkXNJO1aIgUt3YdyHBzF9cw7GfXgQA5fuqrN63NHbr5ojpiaqy+5Yl/ct8Hli7cAXAEoL9cfzvpWnXy7K1k0sbNmNzR7q649B1R0ttuUW4rWvf8POPNMz32QbxQS/Hh4e6Nu3L9LT043HtFot0tPTERMTY/IxMTExtdoDwM6dO43tPTw8EBkZiZMnT9Zqc+rUKbRr107in4CcndBAhggQHtAqKX1AzBSww6aLlZReoNXoR3zNPls6IG0OUyCsYO3AghQXYUorQWiuP/dytlJszkRRaQ9JSUmYMGEC+vXrh6ioKCxfvhzl5eWYOHEiACAxMRFt2rTBkiVLAADTp0/HQw89hHfeeQePPfYYNm/ejCNHjuCDDz4wnnPWrFkYM2YMHnzwQQwePBhpaWn47rvv8Msvv8jxI5JCcWEPiSFmKlbMyJW90wfETgHbfbpYaekF5w7UHfG9V+kFfbsOgxzTpwbAlpQfqS7CrNmNzZ4M/Tn4n6uY+mk2btyqrtNGSFoHWUcxI7+AvnTZ22+/jQULFiA8PBw5OTlIS0szLmorKChAYeHdK6ABAwYgNTUVH3zwAXr37o0vv/wSW7duRVhYmLHNk08+iTVr1mDZsmXo2bMn/v3vf+Orr77CwIEDHf7zkTIpaWSOnIOYgFZJ6QNipoDtPl2sxPSCMoEjbELbOZgS07ZsTfkRsolFa4EXYaYWxMnJTa2CWqUyGfgaWErrIOvYPPI7e/ZsNG8u3WjFtGnTMG3aNJP3mRqtHT16NEaPHl3vOZ9//nk8//zzUnSPGiAljcyRcxAT0Ab4eApqK7SdLTRaHfzv88DEB9pja85FXCu/bbwv0MRInGF69t5RO1NtRbGYXqDSpxd0ewxQu1n3PaxRflnadg6khAWVBoZFoEWllVi0/bjVi9WAhr9tuJIujl2JzcGvIQWByFnxjw+JJWYqVit09M3Og3SmgqNm3u54MrwNYkMDzU4B22W62GJ6gU6e9ALvFtK2cxAlpW2Zep2ZY2pgwVT1FKkvwoRWaHEEpZVicxWKyvklkgP/+JBYYnJnhRbZv1JeZbmRlcwFR9fLq7Fu/1lEWvjwN0wXS+ZmseU2YtpJxVdgECW0nQPYWgpMSuZeZ5YYBhYsjV5LcRGmpBFyQHml2FyFonJ+ieQgZU4ZuQYx+bByX1wpMqfdp5XlNmLaSaXdAP2Cu/r4tdG3UwhbS4FZ69784tt3tGZfZ5a09PUSlBtsa86uUkoO1qS0UmyugsEvuTz+8SFrCC2fJPfFlVzBUb2MQWY9vxU5gky1m77SRH39in/TsXnIFsiRtmWqLGT/JemCUh1qMrz2+7ZravcLNEVeBP6P0kqxuQLJ0h6KioqQmZkJAIiKiqqz7TCRktltYQ81aEKmYqVYsGNLjqLQoOeH/416OST/0RBkfp4ImPutyBVkho4EnvrYRAm2Nvo+KWyHN0fPLJhLbai5eFKImq/9rHPX7b7oWOkLm6VK61BSPrOSSRL8pqamYsGCBYiNjYVOp8Mrr7yCv//97xg7dqwUp3dpfCE7jtLqQJJzEJIPa8vFla05ikKDno8zzuHjjHOOy380G2QGyR9kdnsM8PLXb7qhAtBuoH7hnYJGfA0cmTNa3+ipWDVf+9tyLgh6jC2j13IubBb6Oe6mViGqQzNj28Nnron6DFJaPrOSSRL8Ll26FJmZmWjatCkA4Pr163j44YcZ/NqIL2THk3xhD7kkcyvWxV5cSbGK3xAcCZ2SdmiFgNCR+kDz3AH94jafVvpUBzmDTJMbb3wq38YbFjiyFJil0dP6qAA08/bA6491R6D/fbVe+44YvZYr917M57gtn/lKqvjhDCTJ+dVqtfDx8TF+7ePjA61WK8WpXZYSE/OJyLK0Y4V44M30WvmQD7yZLnrBjlQ5im5qFUb2Fv6h5/D8R7WbflS155/lH11V4sYbAljKGR0WGijJ5hfWjooaXuWLnwzDkxH313ntGy7Q6mNrXryY3HupNgsR8zluy2e+kvOZlUqSkd9nn30WAwYMwJ/+9CcAwNdff43x48dLcWqXJOSFPO+boxjSrRU8GnHNIpFSpB0rxJRN2XWOF5VWYcqmbKwRMfoiVY6iRqvDt7niLpblzn+UhVI33hDI3MzCzrwiDFy6S5IZRKGjos283XGt/O6uZZbSewwXaGv3nDF7zpG9W1scva4vvUDoCPnOvCJJZlzFlKDD//5vbbk6peczK5Ekwe/s2bMxdOhQ7N+/HwCwevVq9O3bV4pTuyQhU0vXyqvRf0k6/vFkGKcyiBRAo9VhztdH620z5+ujguutSpWjaMtUtUtt7KLUjTdEuDdtS+qpcKH5xbtnDUbWueuC03uEXKB9m1uI1+K7mz2PkJQBS7n3ACT7fYmtsmJL8MqNmsSTJPh97bXXMG/ePPTr1w+APud3zpw5ePPNN6U4vcsR+gK9Vn6buTxECnEw/ypuVFTX2+ZGRTUO5l/FAyEBFs8nVY6iLR94LrWxi1I33hCh5shngLcnFn4r7eYXQkdPPRqpRY0wCrlAqy/4ExPkmxshB4CBS3dJ9vsSU2XFz0tYKGbunHLXEjdF6Yv1JQl+d+7ciWXLlhm/btq0KX766ScGv1YS+wJ11O49RGRexn+uCG4nJPiVahW/NR94LrmrlFI33hBIzLbCgPVT4fYoC2nLyKU1O9yZWtickX9V0OjrwfyrUKtVFoM6MVVWhDJ3TqXtEucMi/UlCX61Wi3Kysrg6+sLACgtLUV1df0jIGSepRdyTczlIZKXYYTjVPFNgY+QdpTN0kWvmL8nYs/doBg23igthOm8X5X+fgXt7mZg7bbCgHUzA1KXhbRl5FKqfFehv4epqdm4cetufGMuqIvq0AxNGrtbnA0SwlLw6siKH5Y4S9UJSVZLTZ8+HQMHDsQ//vEPLF68GIMGDcLMmTOlOLVLqrnjmFDM5SFyvJo7Xf2UJ2w63JpRNlt2fqpvB0NTXHZXKePuboDZvR4VtrsbYHvtXWunwm3dargmW3ZBlCrfVejvoWbgC9i/ApPQ4FUJu8Q5U9UJSUZ+n3/+eURFReH//u//oFKpkJqaih49ekhxapcVH9YaM2K74L2fTwlq71K5eUQKYM1oW9PG7ujfUdwMjRSjbOamqlv7e2Heo91xqbQS565VoF2zxhgf0941q8hoNcB9TYH+LwG/bQEqrt69Twkbb5hh7YJGJaW22DJyKVW+q9gZEgNzqRWHz1yTZNTXVDqJuXxauTdqcqaqEzYFv+fOnUOTJk3g7++PsLAwFBcXY+vWraisrETnzp3h6ekpVT9dRs0XdbVGWK3kJo3dFfEHjMhVWDvatmRUT6s+iKTYfMXUB+P18ios2nGi1gfWv/edUVRunkOY2tiicQDQ6ymg66Pyb7xRD2tm/ZSY2mJtLrFU+a71BeCWmArqbJmNnTa4M0Ja+ZgMXi3l08q5UZMzVZ2wKfgdPXo0tm7dCn9/f2RlZeGpp57CvHnzcOzYMfzlL3/Bhg0bJOqmaxC7YMFg4oAOivkD1hApfdUqOZ7Y0TalLPao+cGYdqwQU1N/VXxunt0ZNra49zdRcRU4uBoIjlFs4AtYN+tny+I0e7Jm5FLKfFdzAXiT+9zrpDuYUjOos2U21t1NhSfC29Q5rvR8WiVWnTDHpuC3srISQUFBAIBPPvkEL774Il555RXodDr06tVLkg66CmsXLDRp7I5pQzrbpU/kHKtWyfHEjFw09/bA7lmDFZVKYM0K+QbJyTe2AISNfLby88Q7T4Xjys0qxV/AWzNyKWUFClMBuFanwzP/PmTxsTWDOmvTKABg44GzCG7WuNY20M7wnlVa1Yn62BT8arVaaLVaqNVq/Pzzz1ixYgUAQKVS5ptKqWxZsPCmldOoZJnSr7JJPmJGLq6W30bWueuy57jV5Ey5eXbVQDa2sDTyuXBkDzzQ2XJ5PWcmZb7rvQG4RquzGMg2aewOrVYHjVYHN7XKpjSKaxXVmPl5LoC7gy3+93ko/j2rpKoTltg0FDFmzBgMGzYMY8aMgVqtxuDBgwEA//nPf4xlz8gyaxYstPb3ErVVKonjTKtWyfEMIxxCFZXcsmNvxHOm3Dy7ErphRZl9VvJLRfaV/loNcGYvcPRL/b9ajX2/nxlSVqC497yWKqbcqKjGMx8dwsClu4yVH8w9L2IYBlt+zisS1F7u96zsr0WBbBr5nT9/PmJjY1FUVIRhw4ZBrdbH0nfu3MHKlSsl6aArEPpinTa4E0Ja+Sp+2qoh4MgY1cfwYThlU7ag9ot2nMB9Hm6K+cPvTLl5diV0w4q0OUAjL0VWezCQbaW/icWCt7xaoSA6GZ0ferrBfE6ZS624170zg/c+L1fK9ItMhTKkNHyTc0FQeyW8Z+WuOiGEzaXOYmJi6hzr0qWLrad1KUJfrA90bsFAy0E4MkaWxIe1xvtPR2DaZ9mwNAFwXWFbkTtTbp5dWdzY4n8qruoXxT31saIDYIev9DezWNDzVjFCfvkr5u0/i4cTnpf8NS/XImRDUHfwP1cx9dNsk4vgTOXf1nxeNFod/r3vjKhcYB2Aa+XVaObtgevlt53iPStn1QkhbF6BcejQIZw5cwYAsG/fPrz99tv47rvvbO6YK7GlwDfZB0fGSIhHe7XGynF9LLZTWqpMfdO4SsvNs6t6N7YwIW2ObFP6ilPPYkHDy+Zv1R9h6qYjkm4AUXNjmembczDuw4O1Ug3szU2tglqlqrf6Q82ZQVOPF7PpTE0J4UEmH1fzPQvot2relnMBGflXFfH3RolsCn5nzJiBV155BePGjcPrr7+O+fPnAwDWrl2LGTNmSNE/l8APIuXhBQkJ9WivIKx5NgLNvN3rbVffB6Ic5MrN02h1yvpwDh2pH9FtbGmUqsbiN7K4WFCtAoJUVxGp/l2yiz7DIuR70w5q7rLmiNeXrTOD1uYCDwsNrPc9C0DWCwNnYlPaw88//4yjR4+iqqoK999/Py5cuABPT0/MnDkT4eHhEnXRNUhZqoVs50yrVkl+8WGtcatai5lbciy2VVKqjKNz8xRbOjB0JHCnEvh6suW2JhbJuWQtcIGLBVviBg5KsD5CSKmvOV8fxcJv81BUat/XlxQzgzXfe0WllVi0/TiulZseTa6Z0uCmVpl8z+7MK2J1IhFszvm9c+cOqqqqUF1djcrKSnh6ekKj0UCj4dSQWM6QJO5KeEFCYgT6OWeqjBS5eUKCP8WXDvQV+L3vWSSn2IDezjTeLSGk8vElNNH/a+NFn5BFyPrthGsHkPZ4fUm5q5zhvXefuxov/W8BraXBFlOl2OSoAezMF302Bb+TJk1C9+7dodFosHjxYowZMwYhISE4cOAARo0aJVUfXYrSk8RdDS9ISChXXUQmJPjTaHWY8/VRRRfot7z4TaW/v90A4xHFB/R2knasEIu+rcQXumYIxDWYesq0OqAIzXFY2w2A7Rd91gbP9nh92WNm0JbBFjmqEzn7RZ9Kp9PZlBBz5coVAEBAQABu3LiBn3/+GW3btkV0dLQkHZRbaWkp/P39UVJSAj8/P7m7Q0QKZgiGANMfiA0tGDIX/N378674+RTe+/m0xfN9Nrm/vBf/xuoFgMlnsEa1B41Wh4FLd5kNOgwXO/tmD2lQF8s1n/M49WGsdl8OALUCYEOa7UvVM/CTNkqS30NG/lWM+/Cg9R2H9K8vewSA1oymbsu5gOmbcyyee8XYcJPbJosl9H3vaGLiNZvTHgIC7u4a06RJE/z5z38GoK8C0VACYCIiIQyjN/fmHTbEVBmhU61DurXC+v1nBZ1T9nxow+K3e+rWwi8IiH+zVpkzV6wFfu9z/qM2Ci9Vz0Cy+8cIwt2FnEVojpTq8fhJGwVAmvURtmwXbCD09SU0ALXHzKA1s7+OrE7kDNssC2Fz8GvO6NGjUVBQYK/TExEpWO2PBhsn2BRJaPD3ScbZestC1aSIfOjQkUC3x/TVDG4W63N82w3Ql0WrwRVrgd/7nKuhRQl8sLR6LJqrSnFV54diNMNhbTdooZZ0GtyW7YINhLy+xI7mKiFV0ZEpVw3los+m4Pepp54yeVyn0+HaNevK+axatQpvvfUWioqK0Lt3b/zrX/9CVFSU2fZffPEF5s+fj7NnzyIkJARLly7Fo48+arz/ueeew8aNG2s9Ji4uDmlpaVb1j4jIHHPTgcWlVQ0uB1RoUHfuWoWgdk0auysnH1rtBnQYVG8TV6wFXvM5j1Mf1o/4qu5+1l/UNUNKdSKejemAR8JaS74+wlxebNPG7qi6o0XFbdML7YUGf47O4ZZqwZgjqxM1lIs+m0udffLJJ/Dx8al1XKfTYc+ePaLPt2XLFiQlJWHNmjWIjo7G8uXLERcXh5MnT6Jly5Z12h84cADjxo3DkiVL8PjjjyM1NRUJCQnIzs5GWFiYsV18fDzWr19v/NrT01N034iI6tNQpgOFEhrUtWvWWFC7iQM6ONXvxRUXOBqe85q5vjUF4hpWuy/Hab8u6NoprM79UqiZarAzrwhbcy7iWvlts+2FBn+Ofv9KnS/sqOpEDeWiz6bg9+GHH4avry8efPDBOvf16tVL9PneffddTJ48GRMnTgQArFmzBjt27MC6deswZ86cOu1XrFiB+Ph4zJo1CwCwaNEi7Ny5EytXrsSaNWuM7Tw9PREYGCi6P0REQjWU6UChhAZ/42PaW9zOtUljd0wb0tmOvRVAq7GY6lCTK9YCj+rQDG383JFc9TEA1KnyoFYBWgBdfl0MPDS23t+fLdzUKpTcuo31+89aTH8QGvw58v1rrxFmR1QnaigXfTbt8Pb111+bDHwBYOfOnaLOdfv2bWRlZSE2NvZu59RqxMbGIiMjw+RjMjIyarUH9CkN97b/5Zdf0LJlS3Tt2hUvvfQSrl69arYfVVVVKC0trXUjIrKkoUwHCiV0Z0qPRmqL27m+OaqnvEFi3rfA8jBg4+PAV5P0/y4P0x+vh1y75MnFTa3Ce/0rEKQyXd4M0AcVKjvvhFffKK1Bk8bu+HRSNPbNHiLoeXDU+9fSCDNg2zbohhzkJ8LbIKZTc8nfVw1lR1qrg9/q6mqcP38eJ0+etDq/t6YrV65Ao9GgVavaBcRbtWqFoqIik48pKiqy2D4+Ph4ff/wx0tPTsXTpUuzevRuPPPKI2U04lixZAn9/f+Otbdu2Nv5kROQKGsp0oBhCgz9z7Vr7e2GN3EGiobzZvVv1lhbqjwsIgPfNHoLPJvfHirHh+Gxyf8EBlzOKanFHWEOBO8BZw9IoLaDf8EKtVgkOwhz1/hUzwqxUDeGiT1TaQ1lZGTZt2oTNmzfj8OHDuH37NnQ6HVQqFe6//34MHz4cL774IiIjI+3VX9HGjh1r/H/Pnj3Rq1cvdOrUCb/88guGDh1ap/3cuXORlJRk/Lq0tJQBMBFZ1FCmA8USOtWqyA1jtBp9WbP6Mj3T5uirP1hIgWgIqSyC3LPDnc3trGCPUVpHvX8bygyRIt/PIggOft99910sXrwYnTp1wogRIzBv3jwEBQXhvvvuw7Vr13Ds2DHs3bsXw4cPR3R0NP71r38hJCREcEcCAgLg5uaG4uLaV4vFxcVm83UDAwNFtQeAjh07IiAgAH/88YfJ4NfT05ML4ohINFfKATW1Sl1I8Ke4IPHcgbojvrXoAMMUvoXqDy7Dip3wpGaPUVpHvX8b0gyR4t7PIggOfjMzM7Fnzx706NHD5P1RUVF4/vnnsWbNGqxfvx579+4VFfx6eHigb9++SE9PR0JCAgBAq9UiPT0d06ZNM/mYmJgYpKenY8aMGcZjO3fuRExMjNnv89///hdXr15F69bKH5YnIuW7NxBc9XQfLNpxwq4rruXk7Nua1iJ0at6OU/hiSVUey2pqNyB+6f92wjMTJsa/abfFboD9RmkdUTGhb7um+oWB9aT0qlX6dmQ/goPfzz77TFA7T09PTJkyxarOJCUlYcKECejXrx+ioqKwfPlylJeXG6s/JCYmok2bNliyZAkAYPr06XjooYfwzjvv4LHHHsPmzZtx5MgRfPDBBwCAmzdvIiUlBX/6058QGBiI/Px8vPbaa+jcuTPi4uKs6iMRkYG5QHD+Y6Fo6u0heYAid+Dj6DqodqeAKXwxFHPhIWInPHuw5yitvafzs85drzfwBfSBcda56047quoM7LbDmzXGjBmDy5cvY8GCBSgqKkJ4eDjS0tKMi9oKCgqgVt9dozdgwACkpqbi9ddfx7x58xASEoKtW7caa/y6ubnht99+w8aNG3Hjxg0EBQVh+PDhWLRoEVMbiMgm9QWCU1P1geAT4W0k/X5yBj4Nso6xAqbwhVLchYfAnfDsxZ6jtPaczm8oOb/OTqWzYt/NkpISzJo1C7t27YK7uzt27drVYNMISktL4e/vj5KSEvj5+cndHZKI3CNo5Nw0Wh0GLt1ldtW2Ydp13+whkhXENxX4GM7siMAnI/8qxn140GK7zyb3d64RK0O1BwAmxxCf+tjuI5mWOPr1JpjI2sj24Gx/yxvs+0gBxMRrVo38Tp06Ffn5+Vi2bBmeffZZ3Lp1CwAwc+ZMdOrUyWyOLpESyD2CRs7PkQXxlTLi2mBHrMxN4TduDvR6CrivqT7Ic3BQV5MiN1DJ+9ZM2sNSh14sONuiK+b8KoNVdX5/+OEHvP/++xg1ahTc3O7+QYiLi8PGjRsl6xyR1AwjaPd+kBimDtOOFcrUM3ImjgwElVIXtCGtUq8jdCQw4xgwYTvQ/6/6wLfiCnDwfcEbXtiT4i48bKyN7Go0Wh0y8q9iW84FfJJxVnDOL9mPVSO/Op0Ovr6+dY6HhITg9OnTNneKyB6UMoKmZM42hSgXRwaCSgl8GnwdY7UbcOs6cHA16uT/GoI6mVIghL6OAnw8kZF/1b7vX4lqI7sKUzONQjjdDIqTsSr4feSRR/Dpp58iOTm51vHy8nKoVPygJGVS5NShgjAdRDhHBoJKGXFt8HWMFRzUCXm9NWnsjlc+z0FRaZXxuF3evw6sjezsF+PmcvWFcMoZFCdiVdrDkiVLsGrVKqSkpBh3eKusrMSiRYsQEREhdR+JJKGUETQlYjqIOI7c394Q+NR3piaN3aHV6qCxNJ9qo4awralZYoI6B7P0etMBuF5RXSvwBez0/nVQbeS0Y4UYuHQXxn14ENM352DchwcxcOkup/lbVN9MY31UAJp7e6Co5BYy8q/a/T3tqqwa+Q0ODsaBAwcwdepUVFRUICoqCmVlZfDz88P3338vdR+JJKGUETSlYTqIdRxREB+of8TV4EZFNZ756JBDRuolr4OqgIoBABS/4UV9r7db1RrcqKiu8xi7vH8dUBtZcWXdrGBpptEcHYCr5bcx8/NcAJx9sxdRwW9FRQUaN24MAOjcuTN+/PFHFBQUIDc3F+7u7oiOjkbTplyhSMrU4HMWrcR0EOs5an97c4HPvRwVHEi2wl4hFQMAOMWGF6Zeb1qtDs98dMjsYyR//9q5NnJDuRiXagbRmQJ+ZyIq7cHf3x9XrlypdSw4OBgjRoxAfHw8A19SNEdOVTsTpoPYxhAIPhHeBjGdmtvt9RMf1hr7Zg/Bpy9Eo8l97ibbGAKGlO/ylD9dqrSKAYagzmyCiQrwayP7hhf3vt6ulFdZfhAkfP8atjcGYPYvqQ3bGyuluomthM4gzn+sO94bE45m3g3gPe1ERAW/Go0GWq3W+PWgQYNQXKycPc+JLGnQOYtWYjqI83BTq6BWqXDjVt0pbgOnCA4sLi6DfnGZVuO4Ptk5qLMXWd6/htrIfvf8vfQLsrkihlIvxmuWKxOSi2spV18FfUrDcw90QKCfF66VO/l72snYtL1xTk4OysvLpeoLkUM4aqraWTAdxLkoNTgQxYEVA0Qxt+GFX5A+8JV5pzdTZHv/2ml7YyVejFtTCUdMdZQG8Z52MjYFv0TOytl2BbKnBl/CqoFRYnAgmpIXl9kpqLMXWd+/ajfJL06UdjFuy+I7oYtiG8R72smIDn5TU1Px4IMPomfPnvboDxHJwFGVC8h2SgsOrKL0xWV2COrsqSG9f5V0MS7F4jshM40N4j3tZFQ6nU5wBvVDDz2EnJwclJWVwd3dHXfu3MHTTz+NQYMGoU+fPujVqxc8PT3t2V+HKy0thb+/P0pKSuDn5yd3d4jsytmLyrsKw2gUYDo4UHz+ulaj3zLYUsWAGUcVO+KqRA3p/auETXcy8q9i3IcHLbb7bHJ/m2cSnf49rQBi4jVRwa/B6dOnkZWVhezsbOPtxo0baNSoEbp3747c3FyrO680DH6JSImUEBzYxFDtAYDJj3uZthIm5ZA7mN+WcwHTN+dYbLdibDieCG9j8/dz+ve0zMTEa1bl/IaEhCAkJARjx441Hjtz5gyOHDmCX3/91ZpTEhGRCE6/cNMJF5eRY8m9NsPRubhO/552IoJHfgsKChAcHCz4xBcuXECbNrZfCcmNI79ERHaklB3eiO6h0eowcOkui7m4+2YPYYCqAGLiNcF1fiMjI/GXv/wFmZmZZtuUlJTgww8/RFhYGL766ivhPSYiItdkWFzW88/6fxn4kkJwY6SGS3DaQ15eHhYvXoxhw4bBy8sLffv2RVBQELy8vHD9+nXk5eXh+PHjiIiIwLJly/Doo4/as99EREREdtWQKmnQXaIXvN26dQs7duzAvn37cO7cOdy6dQsBAQHo06cP4uLiEBYWZq++yoJpD0RERK5N7sV3QjlLP+3B7tUeXAmDXyIiIlI6V68WYZecXyIiIiJSHkOd4JqBL3B3J7q0Y4Uy9UyZJAt+S0tLsWfPHvzzn/+U6pREREREVA9LO9EB+p3oNFpO9BtYVee3oKAAOTk5tW7nzp2DTqeDt7c3/va3v0ndTyIiIiK6x+Ez1+qM+NakA1BYUonDZ67JWjdZSUQFv0OGDEFubi6uX78Of39/hIaGIiwsDAUFBfjoo48wdOhQtG3b1l59JSIiIqIaLpWZD3ytaecKRKU97Nu3D1OmTMH58+dx/fp17N+/H2vXroVKpUJUVBQDXyIiIiIHcvROdA2BqOD30KFD2Lt3L6ZOnYpTp07Zq09ERGSGRqtDRv5VbMu5gIz8q8zjI3JxUR2aobW/V52NOAxU0Fd9iOrQzJHdUjRRaQ99+vTBnj17kJqairi4ODz66KNITk62V9+IiKgGVy9lRER1GXaie2lTNlRArYVv3InONKuqPTz99NM4fvw4mjZtih49ekCr1UKj0UjdNyIi+h+WMiIicww70QX6105tCPT3wupnI3hxfA+bN7n4z3/+g5kzZyIjIwOvvfYapk6divvuu0+q/smOm1wQkdw0Wh0GLt1ldkW3CvoPuX2zh3B0h8iFcYc3B21y0bFjR2zbtg2bNm3C+vXr0bFjR5vOt2rVKrRv3x5eXl6Ijo7G4cOH623/xRdfoFu3bvDy8kLPnj3x/fffm207ZcoUqFQqLF++3KY+EhE5kphSRkTkutzUKsR0ao4nwtsgplNzlwl8xZJsk4vhw4cjNzcXs2fPtvocW7ZsQVJSEpKTk5GdnY3evXsjLi4Oly5dMtn+wIEDGDduHCZNmoRff/0VCQkJSEhIwLFjx+q0/eabb3Dw4EEEBQVZ3T8iIjmwlBERkXQk3d64UaNGmDFjhtWPf/fddzF58mRMnDgRoaGhWLNmDRo3box169aZbL9ixQrEx8dj1qxZ6N69OxYtWoSIiAisXLmyVrsLFy7g5Zdfxqeffgp3d3er+0dEJAeWMiIiko6kwa8tbt++jaysLMTGxhqPqdVqxMbGIiMjw+RjMjIyarUHgLi4uFrttVotxo8fj1mzZqFHjx726TwRkR2xlBERkXQUE/xeuXIFGo0GrVq1qnW8VatWKCoqMvmYoqIii+2XLl2KRo0aCd5yuaqqCqWlpbVuRERyMpQyAlAnAGYpIyIicRQT/NpDVlYWVqxYgQ0bNkClEvahsGTJEvj7+xtv3LWOiJSApYyIiKQhapMLewoICICbmxuKi4trHS8uLkZgYKDJxwQGBtbbfu/evbh06RKCg4ON92s0GrzyyitYvnw5zp49W+ecc+fORVJSkvHr0tJSBsBEpAjxYa0xLDTQZUsZERFJQTHBr4eHB/r27Yv09HQkJCQA0OfrpqenY9q0aSYfExMTg/T09FqL7Hbu3ImYmBgAwPjx403mBI8fPx4TJ040eU5PT094enoavzaUQWb6AxEpRY8W7ujRQr94t/xmmcy9ISKSnyFOE7R9hU5BNm/erPP09NRt2LBBl5eXp3vxxRd1TZo00RUVFel0Op1u/Pjxujlz5hjb79+/X9eoUSPd22+/rTtx4oQuOTlZ5+7urjt69KjZ79GuXTvde++9J7hP58+f10FfRpM33njjjTfeeOONNwXfzp8/bzG2U8zILwCMGTMGly9fxoIFC1BUVITw8HCkpaUZF7UVFBRArb6bpjxgwACkpqbi9ddfx7x58xASEoKtW7ciLCxMsj4FBQXh/Pnz8PX1FZw3LJQhpeL8+fPcPc7J8Llzbnz+nBufP+fF5865Kfn50+l0KCsrE7Sfg83bG5P1uHWy8+Jz59z4/Dk3Pn/Oi8+dc2soz1+DrvZARERERFQTg18iIiIichkMfmXk6emJ5OTkWtUlyDnwuXNufP6cG58/58Xnzrk1lOePOb9ERERE5DI48ktERERELoPBLxERERG5DAa/REREROQyGPwSERERkctg8CuTVatWoX379vDy8kJ0dDQOHz4sd5dIoD179mDEiBEICgqCSqXC1q1b5e4SCbRkyRJERkbC19cXLVu2REJCAk6ePCl3t0iA1atXo1evXvDz84Ofnx9iYmLwww8/yN0tstKbb74JlUqFGTNmyN0VEmDhwoVQqVS1bt26dZO7W1Zj8CuDLVu2ICkpCcnJycjOzkbv3r0RFxeHS5cuyd01EqC8vBy9e/fGqlWr5O4KibR7925MnToVBw8exM6dO1FdXY3hw4ejvLxc7q6RBffffz/efPNNZGVl4ciRIxgyZAieeOIJHD9+XO6ukUiZmZlYu3YtevXqJXdXSIQePXqgsLDQeNu3b5/cXbIaS53JIDo6GpGRkVi5ciUAQKvVom3btnj55ZcxZ84cmXtHYqhUKnzzzTdISEiQuytkhcuXL6Nly5bYvXs3HnzwQbm7QyI1a9YMb731FiZNmiR3V0igmzdvIiIiAu+//z7eeOMNhIeHY/ny5XJ3iyxYuHAhtm7dipycHLm7IgmO/DrY7du3kZWVhdjYWOMxtVqN2NhYZGRkyNgzItdTUlICQB9EkfPQaDTYvHkzysvLERMTI3d3SISpU6fiscceq/UZSM7h9OnTCAoKQseOHfHMM8+goKBA7i5ZrZHcHXA1V65cgUajQatWrWodb9WqFX7//XeZekXkerRaLWbMmIEHHngAYWFhcneHBDh69ChiYmJQWVkJHx8ffPPNNwgNDZW7WyTQ5s2bkZ2djczMTLm7QiJFR0djw4YN6Nq1KwoLC5GSkoJBgwbh2LFj8PX1lbt7ojH4JSKXNHXqVBw7dsyp89ZcTdeuXZGTk4OSkhJ8+eWXmDBhAnbv3s0A2AmcP38e06dPx86dO+Hl5SV3d0ikRx55xPj/Xr16ITo6Gu3atcPnn3/ulGlHDH4dLCAgAG5ubiguLq51vLi4GIGBgTL1isi1TJs2Ddu3b8eePXtw//33y90dEsjDwwOdO3cGAPTt2xeZmZlYsWIF1q5dK3PPyJKsrCxcunQJERERxmMajQZ79uzBypUrUVVVBTc3Nxl7SGI0adIEXbp0wR9//CF3V6zCnF8H8/DwQN++fZGenm48ptVqkZ6eztw1IjvT6XSYNm0avvnmG+zatQsdOnSQu0tkA61Wi6qqKrm7QQIMHToUR48eRU5OjvHWr18/PPPMM8jJyWHg62Ru3ryJ/Px8tG7dWu6uWIUjvzJISkrChAkT0K9fP0RFRWH58uUoLy/HxIkT5e4aCXDz5s1aV7tnzpxBTk4OmjVrhuDgYBl7RpZMnToVqamp2LZtG3x9fVFUVAQA8Pf3x3333Sdz76g+c+fOxSOPPILg4GCUlZUhNTUVv/zyC3788Ue5u0YC+Pr61smt9/b2RvPmzZlz7wReffVVjBgxAu3atcPFixeRnJwMNzc3jBs3Tu6uWYXBrwzGjBmDy5cvY8GCBSgqKkJ4eDjS0tLqLIIjZTpy5AgGDx5s/DopKQkAMGHCBGzYsEGmXpEQq1evBgA8/PDDtY6vX78ezz33nOM7RIJdunQJiYmJKCwshL+/P3r16oUff/wRw4YNk7trRA3ef//7X4wbNw5Xr15FixYtMHDgQBw8eBAtWrSQu2tWYZ1fIiIiInIZzPklIiIiIpfB4JeIiIiIXAaDXyIiIiJyGQx+iYiIiMhlMPglIiIiIpfB4JeIiIiIXAaDXyIiIiJyGQx+iYgaqA0bNkClUkGlUmHGjBmSn//hhx82nj8nJ0fy8xMR2QODXyKiBszPzw+FhYVYtGiRoPZfffUV3NzccOHCBZP3h4SEGHc1/Prrr3H48GHJ+kpE5AgMfomIGjCVSoXAwED4+voKaj9y5Eg0b94cGzdurHPfnj178Mcff2DSpEkAgGbNmjnt9qZE5LoY/BIROYGJEyeiRYsWWLZsmfHY0aNH4ebmhu+++07UuaqqqvDqq6+iTZs28Pb2RnR0NH755RcAgLu7O8aPH48NGzbUedy6desQHR2NHj162PKjEBHJisEvEZETePfdd/GPf/wD8+bNQ35+PgBg/vz5iIyMxIgRI0Sda9q0acjIyMDmzZvx22+/YfTo0YiPj8fp06cBAJMmTcLp06exZ88e42Nu3ryJL7/80jjqS0TkrBj8EhE5gaZNm2Ly5Ml4+OGHsWnTJmRmZmLbtm1YvHixqPMUFBRg/fr1+OKLLzBo0CB06tQJr776KgYOHIj169cDAEJDQ9G/f3+sW7fO+LjPP/8cOp0OY8eOlfTnIiJytEZyd4CIiIRLTEzEokWLsH//fgwePBhDhw4V9fijR49Co9GgS5cutY5XVVWhefPmxq+ff/55zJw5E//617/g6+uLdevWYfTo0YJzh4mIlIrBLxGRExk1ahT++te/4o8//sCBAwdEP/7mzZtwc3NDVlYW3Nzcat3n4+Nj/P/YsWMxc+ZMfP7553jwwQexf/9+LFmyxOb+ExHJjcEvEZET8fHxQadOneDr64uYmBjRj+/Tpw80Gg0uXbqEQYMGmW3n6+uL0aNHY926dcjPz0eXLl3qbU9E5CwY/BIROZG0tDT89ttvaNasGW7fvg0PDw9Rj+/SpQueeeYZJCYm4p133kGfPn1w+fJlpKeno1evXnjssceMbSdNmoRBgwbhxIkTmD17ttQ/ChGRLLjgjYjISeh0Orz++usYPXo0PDw8sGPHDqvOs379eiQmJuKVV15B165dkZCQgMzMTAQHB9dqN3DgQHTt2hWlpaVITEyU4kcgIpIdR36JiJzE119/jZycHKSmpmLt2rX45JNP8OSTT4o+j7u7O1JSUpCSkmKx7e+//25NV4mIFIsjv0RETkCr1WLBggVITExEly5dkJiYiG+//RYbN27E9evXzT6upKQEPj4+dklbeOSRR7jhBRE5HZVOp9PJ3QkiIqrfJ598ghdeeAEnT55E+/btAeg3uVi1ahUmTJiA9957r85jysrKUFxcDABo0qQJAgICJO3ThQsXcOvWLQBAcHCw6PxjIiI5MPglIiIiIpfBtAciIiIichkMfomIiIjIZTD4JSIiIiKXweCXiIiIiFwGg18iIiIichkMfomIiIjIZTD4JSIiIiKXweCXiIiIiFwGg18iIiIichn/H5NnM39r/73lAAAAAElFTkSuQmCC", - "text/plain": [ - "Figure(PyObject
)" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "# get a bitmask for the lines of Fe I vs Fe II\n", - "neutrals = [spec.charge == 0 for spec in lines.species]\n", - "\n", - "figure(figsize=(8, 2))\n", - "scatter(lines.ExPot[neutrals], (lines.A_18Sco - lines.A_sun)[neutrals], label=\"Fe I\")\n", - "scatter(lines.ExPot[.! neutrals], (lines.A_18Sco - lines.A_sun)[.! neutrals], label=\"Fe II\")\n", - "ylabel(L\"A(Fe)_\\mathrm{18 Sco} - A(Fe)_\\mathrm{sun}\")\n", - "xlabel(\"χ [eV]\")\n", - "legend()\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a1f62ac9-37f0-4b91-809f-4de341a4b250", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.8.5", - "language": "julia", - "name": "julia-1.8" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.8.5" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/misc/Tutorial notebooks/EW fitting/Table 1.dat b/misc/Tutorial notebooks/EW fitting/Table 1.dat deleted file mode 100644 index f2d00f01..00000000 --- a/misc/Tutorial notebooks/EW fitting/Table 1.dat +++ /dev/null @@ -1,347 +0,0 @@ -Title: 18 Sco: a solar twin rich in refractory and neutron-capture elements. - Implications for chemical tagging. -Authors: Melendez J., Ramirez I., Karakas A.I., Yong D., Monroe T.R., Bedell M., - Bergemann M., Asplund M., Maia M.T., Bean J., do Nascimento J.-D., - Bazot M., Alves-Brito A., Freitas F.C., Castro M. -Table: Adopted atomic data and equivalent widths -================================================================================ -Byte-by-byte Description of file: apj497906t1_mrt.txt --------------------------------------------------------------------------------- - Bytes Format Units Label Explanations --------------------------------------------------------------------------------- - 1- 9 F9.3 0.1nm Wave Wavelength in Angstroms - 11- 14 F4.1 --- Ion Ion numeric identifier (1) - 16- 21 F6.4 eV ExPot Excitation potential - 23- 28 F6.3 [-] log(gf) Log of the oscillator strength times - the statistical weight - 31- 38 E8.2 --- C6 Damping constant - 40- 44 F5.1 10-13m EW-18S The 18 Sco differential equivalent width (2) - 46- 51 F6.2 10-13m EW-Sun Solar differential equivalent width (2) --------------------------------------------------------------------------------- -Note (1): Where the number before the digit is the atomic number and the number - after the digit is the ionization level, e.g. 26.1 = Fe II. -Note (2): In units of milli-Angstroms. --------------------------------------------------------------------------------- - 5044.211 26.0 2.8512 -2.058 2.71e-31 74.8 74.3 - 5054.642 26.0 3.640 -1.921 4.68e-32 40.9 40.5 - 5127.359 26.0 0.915 -3.307 1.84e-32 97.5 96.1 - 5127.679 26.0 0.052 -6.125 1.20e-32 18.9 19.1 - 5198.711 26.0 2.223 -2.135 4.61e-32 99.3 98.1 - 5225.525 26.0 0.1101 -4.789 1.23e-32 72.5 72.1 - 5242.491 26.0 3.634 -0.967 4.95e-32 88.3 86.9 - 5247.050 26.0 0.0872 -4.946 1.22e-32 67.4 66.9 - 5250.208 26.0 0.1212 -4.938 1.23e-32 66.3 65.9 - 5295.312 26.0 4.415 -1.49 6.54e-31 31.0 30.3 - 5322.041 26.0 2.279 -2.80 4.29e-32 62.9 61.5 - 5373.709 26.0 4.473 -0.77 7.04e-31 65.2 63.9 - 5379.574 26.0 3.694 -1.514 5.02e-32 62.9 61.5 - 5386.334 26.0 4.154 -1.74 5.27e-31 34.8 33.6 - 5466.396 26.0 4.371 -0.565 4.40e-31 81.3 79.4 - 5466.987 26.0 3.573 -2.23 2.80e+00 35.8 35.2 - 5522.446 26.0 4.209 -1.31 3.02e-31 45.8 43.7 - 5546.506 26.0 4.371 -1.18 3.91e-31 52.6 51.4 - 5560.211 26.0 4.434 -1.16 4.79e-31 53.8 52.0 - 5577.02 26.0 5.0331 -1.455 2.80e+00 11.9 11.2 - 5618.633 26.0 4.209 -1.276 2.90e-31 51.8 50.2 - 5636.696 26.0 3.640 -2.56 5.19e-32 20.7 19.7 - 5638.262 26.0 4.220 -0.81 2.88e-31 78.8 77.6 - 5649.987 26.0 5.0995 -0.8 2.77e-31 37.8 35.9 - 5651.469 26.0 4.473 -1.75 4.83e-31 20.3 18.9 - 5661.348 26.0 4.2843 -1.756 3.24e-31 24.4 23.0 - 5679.023 26.0 4.652 -0.75 8.13e-31 61.7 59.6 - 5696.089 26.0 4.548 -1.720 5.78e-31 14.5 13.7 - 5701.544 26.0 2.559 -2.216 4.95e-32 86.4 84.7 - 5705.464 26.0 4.301 -1.355 3.02e-31 39.7 38.1 - 5855.076 26.0 4.6075 -1.478 5.74e-31 23.8 22.9 - 5905.672 26.0 4.652 -0.69 6.23e-31 61.7 60.1 - 5916.247 26.0 2.453 -2.936 4.29e-32 55.4 55.5 - 5927.789 26.0 4.652 -1.04 6.07e-31 46.0 43.6 - 5934.655 26.0 3.928 -1.07 5.69e-31 79.1 77.9 - 5956.694 26.0 0.8589 -4.605 1.55e-32 52.3 52.8 - 5987.065 26.0 4.795 -0.212 1.55e-32 70.3 68.2 - 6003.012 26.0 3.881 -1.060 4.83e-31 85.5 84.5 - 6005.541 26.0 2.588 -3.43 2.80e+00 23.9 23.0 - 6024.058 26.0 4.548 -0.02 3.88e-31 114.0 111.2 - 6027.050 26.0 4.0758 -1.09 2.80e+00 65.7 64.0 - 6056.005 26.0 4.733 -0.45 6.79e-31 74.2 72.3 - 6065.482 26.0 2.6085 -1.530 4.71e-32 118.8 117.3 - 6079.009 26.0 4.652 -1.10 5.13e-31 47.5 46.5 - 6082.711 26.0 2.223 -3.573 3.27e-32 35.3 34.0 - 6093.644 26.0 4.607 -1.30 4.41e-31 32.8 30.8 - 6096.665 26.0 3.9841 -1.81 5.75e-31 39.3 37.6 - 6151.618 26.0 2.1759 -3.299 2.55e-32 51.2 49.9 - 6157.728 26.0 4.076 -1.22 2.80e+00 64.7 62.7 - 6165.360 26.0 4.1426 -1.46 2.80e+00 47.4 45.4 - 6173.335 26.0 2.223 -2.880 2.65e-32 69.0 68.7 - 6187.990 26.0 3.943 -1.67 4.90e-31 49.4 47.6 - 6200.313 26.0 2.6085 -2.437 4.58e-32 74.6 73.6 - 6213.430 26.0 2.2227 -2.52 2.62e-32 83.9 82.7 - 6219.281 26.0 2.198 -2.433 2.58e-32 91.2 90.3 - 6226.736 26.0 3.883 -2.1 4.15e-31 31.5 30.1 - 6240.646 26.0 2.2227 -3.233 3.14e-32 49.3 48.7 - 6252.555 26.0 2.4040 -1.687 3.84e-32 123.1 121.6 - 6265.134 26.0 2.1759 -2.550 2.48e-32 86.9 85.9 - 6270.225 26.0 2.8580 -2.54 4.58e-32 52.9 52.2 - 6271.279 26.0 3.332 -2.703 2.78e-31 25.3 24.6 - 6380.743 26.0 4.186 -1.376 2.80e+00 54.3 53.2 - 6392.539 26.0 2.279 -4.03 3.38e-32 17.3 16.7 - 6430.846 26.0 2.1759 -2.006 2.42e-32 113.4 114.0 - 6498.939 26.0 0.9581 -4.699 1.53e-32 47.2 46.5 - 6593.871 26.0 2.4326 -2.422 3.69e-32 86.6 86.3 - 6597.561 26.0 4.795 -0.98 4.76e-31 46.2 44.7 - 6625.022 26.0 1.011 -5.336 2.80e+00 15.8 14.9 - 6677.987 26.0 2.692 -1.418 3.46e-32 132.8 129.9 - 6703.567 26.0 2.7585 -3.023 3.66e-32 38.3 37.6 - 6705.102 26.0 4.607 -0.98 2.80e+00 47.9 46.3 - 6710.319 26.0 1.485 -4.88 2.01e-32 15.7 15.7 - 6713.745 26.0 4.795 -1.40 4.30e-31 22.3 21.2 - 6725.357 26.0 4.103 -2.19 4.82e-31 19.2 18.6 - 6726.667 26.0 4.607 -1.03 4.82e-31 48.4 47.4 - 6733.151 26.0 4.638 -1.47 3.41e-31 28.7 27.0 - 6739.522 26.0 1.557 -4.79 2.10e-32 13.2 12.4 - 6750.152 26.0 2.4241 -2.621 4.11e-32 75.3 75.3 - 6752.707 26.0 4.638 -1.204 3.37e-31 39.0 36.8 - 6793.259 26.0 4.076 -2.326 2.80e+00 14.3 13.0 - 6806.845 26.0 2.727 -3.11 3.46e-32 35.8 35.4 - 6810.263 26.0 4.607 -0.986 4.50e-31 53.1 50.6 - 6837.006 26.0 4.593 -1.687 2.46e-32 20.0 18.9 - 6839.830 26.0 2.559 -3.35 3.95e-32 32.5 31.6 - 6843.656 26.0 4.548 -0.86 2.94e-31 64.1 63.3 - 6858.150 26.0 4.607 -0.930 3.24e-31 53.6 52.4 - 5197.577 26.1 3.2306 -2.22 8.69e-33 84.0 80.7 - 5234.625 26.1 3.2215 -2.18 8.69e-33 83.5 80.7 - 5264.812 26.1 3.2304 -3.13 9.43e-33 48.6 45.8 - 5325.553 26.1 3.2215 -3.16 8.57e-33 42.7 41.3 - 5414.073 26.1 3.2215 -3.58 9.30e-33 28.5 26.7 - 5425.257 26.1 3.1996 -3.22 8.45e-33 43.0 41.3 - 6084.111 26.1 3.1996 -3.79 7.87e-33 22.6 20.9 - 6247.557 26.1 3.8918 -2.30 9.43e-33 54.6 52.9 - 6369.462 26.1 2.8912 -4.11 7.42e-33 20.4 18.7 - 6416.919 26.1 3.8918 -2.64 9.30e-33 42.6 40.5 - 6432.680 26.1 2.8912 -3.57 7.42e-33 43.9 42.4 - 6456.383 26.1 3.9036 -2.05 9.30e-33 64.9 62.4 - 5052.167 06.0 7.685 -1.24 2.80e+00 34.0 34.5 - 5380.337 06.0 7.685 -1.57 2.80e+00 20.5 20.6 - 6587.61 06.0 8.537 -1.05 2.80e+00 15.6 15.2 - 7111.47 06.0 8.640 -1.07 2.91e-30 14.9 14.3 - 7113.179 06.0 8.647 -0.76 2.97e-30 24.9 23.6 - 7771.944 08.0 9.146 0.37 8.41e-32 73.0 71.1 - 7774.166 08.0 9.146 0.22 8.41e-32 65.6 63.7 - 7775.388 08.0 9.146 0.00 8.41e-32 51.7 51.1 - 4751.822 11.0 2.1044 -2.078 2.80e+00 11.9 11.6 - 5148.838 11.0 2.1023 -2.044 2.80e+00 11.8 11.7 - 6154.225 11.0 2.1023 -1.547 2.80e+00 36.9 36.9 - 6160.747 11.0 2.1044 -1.246 2.80e+00 54.0 54.3 - 4571.095 12.0 0.000 -5.623 2.80e+00 105.6 106.0 - 4730.040 12.0 4.340 -2.389 2.80e+00 69.3 69.0 - 5711.088 12.0 4.345 -1.729 2.80e+00 106.4 105.6 - 6319.236 12.0 5.108 -2.165 2.80e+00 27.6 25.6 - 6696.018 13.0 3.143 -1.481 2.80e+00 37.0 36.3 - 6698.667 13.0 3.143 -1.782 2.80e+00 21.7 21.2 - 7835.309 13.0 4.021 -0.68 2.80e+00 43.6 43.6 - 7836.134 13.0 4.021 -0.45 2.80e+00 60.9 57.7 - 8772.866 13.0 4.0215 -0.38 9.71e-30 75.0 73.9 - 8773.896 13.0 4.0216 -0.22 9.71e-30 95.3 92.3 - 5488.983 14.0 5.614 -1.69 2.80e+00 21.4 19.8 - 5517.540 14.0 5.080 -2.496 2.80e+00 15.3 13.9 - 5645.611 14.0 4.929 -2.04 2.80e+00 37.7 36.1 - 5665.554 14.0 4.920 -1.94 2.80e+00 43.2 41.1 - 5684.484 14.0 4.953 -1.55 2.80e+00 64.3 62.0 - 5690.425 14.0 4.929 -1.77 2.80e+00 51.2 49.3 - 5701.104 14.0 4.930 -1.95 2.80e+00 41.4 40.5 - 5793.073 14.0 4.929 -1.96 2.80e+00 45.8 44.0 - 6125.021 14.0 5.614 -1.50 2.80e+00 33.8 32.1 - 6145.015 14.0 5.616 -1.41 2.80e+00 40.8 39.1 - 6243.823 14.0 5.616 -1.27 2.80e+00 47.0 44.6 - 6244.476 14.0 5.616 -1.32 2.80e+00 48.3 45.8 - 6721.848 14.0 5.862 -1.12 2.80e+00 48.0 44.8 - 6741.63 14.0 5.984 -1.65 2.80e+00 17.3 16.4 - 6046.000 16.0 7.868 -0.15 2.80e+00 22.2 20.0 - 6052.656 16.0 7.870 -0.4 2.80e+00 12.8 12.3 - 6743.54 16.0 7.866 -0.6 2.80e+00 9.0 8.5 - 6757.153 16.0 7.870 -0.15 2.80e+00 20.1 19.5 - 8693.93 16.0 7.870 -0.44 1.51e-30 13.2 12.7 - 8694.62 16.0 7.870 0.1 1.51e-30 30.3 28.6 - 7698.974 19.0 0.000 -0.168 1.04e-31 157.3 157.0 - 5260.387 20.0 2.521 -1.719 7.27e-32 35.9 34.0 - 5512.980 20.0 2.933 -0.464 2.80e+00 88.1 85.7 - 5581.965 20.0 2.5229 -0.555 6.40e-32 96.5 94.8 - 5590.114 20.0 2.521 -0.571 6.36e-32 94.7 92.4 - 5867.562 20.0 2.933 -1.57 2.80e+00 25.1 24.0 - 6166.439 20.0 2.521 -1.142 5.95e-31 71.9 70.2 - 6169.042 20.0 2.523 -0.797 5.95e-31 95.4 92.7 - 6455.598 20.0 2.523 -1.34 5.09e-32 58.0 56.8 - 6471.662 20.0 2.525 -0.686 5.09e-32 94.4 93.3 - 6499.650 20.0 2.523 -0.818 5.05e-32 88.6 86.5 - 4743.821 21.0 1.4478 0.35 5.97e-32 9.3 9.1 - 5081.57 21.0 1.4478 0.30 2.80e+00 9.6 9.1 - 5520.497 21.0 1.8649 0.55 2.80e+00 7.4 7.4 - 5671.821 21.0 1.4478 0.55 2.80e+00 15.3 15.2 - 5526.820 21.1 1.770 0.140 2.80e+00 78.5 76.6 - 5657.87 21.1 1.507 -0.30 2.80e+00 70.3 68.6 - 5684.19 21.1 1.507 -0.95 2.80e+00 39.4 38.6 - 6245.63 21.1 1.507 -1.030 2.80e+00 36.4 35.3 - 6279.76 21.1 1.500 -1.2 2.80e+00 30.6 30.3 - 6320.843 21.1 1.500 -1.85 2.80e+00 9.7 9.0 - 6604.578 21.1 1.3569 -1.15 2.80e+00 39.1 37.1 - 4465.802 22.0 1.7393 -0.163 3.98e-32 40.6 40.4 - 4555.485 22.0 0.8484 -0.488 4.42e-32 66.6 66.1 - 4758.120 22.0 2.2492 0.425 3.84e-32 45.2 45.2 - 4759.272 22.0 2.2555 0.514 3.86e-32 47.1 46.9 - 4820.410 22.0 1.5024 -0.439 3.78e-32 45.9 44.3 - 4913.616 22.0 1.8731 0.161 3.86e-32 54.5 52.1 - 5022.871 22.0 0.8258 -0.434 3.58e-32 73.0 72.6 - 5113.448 22.0 1.443 -0.783 3.06e-32 28.6 27.5 - 5147.479 22.0 0.0000 -2.012 2.08e-32 37.6 37.5 - 5219.700 22.0 0.021 -2.292 2.08e-32 29.0 29.1 - 5295.774 22.0 1.0665 -1.633 2.58e-32 13.6 13.3 - 5490.150 22.0 1.460 -0.933 5.41e-32 22.9 22.1 - 5739.464 22.0 2.249 -0.60 3.86e-32 8.7 8.5 - 5866.452 22.0 1.066 -0.840 2.16e-32 48.6 48.0 - 6091.174 22.0 2.2673 -0.423 3.89e-32 16.3 15.8 - 6126.217 22.0 1.066 -1.424 2.06e-32 23.1 22.8 - 6258.104 22.0 1.443 -0.355 4.81e-32 52.6 52.3 - 6261.101 22.0 1.429 -0.479 4.68e-32 49.8 49.1 - 4470.857 22.1 1.1649 -2.06 2.80e+00 64.8 64.0 - 4544.028 22.1 1.2429 -2.53 2.80e+00 44.4 41.5 - 4583.408 22.1 1.165 -2.87 2.80e+00 33.8 32.2 - 4636.33 22.1 1.16 -3.152 2.80e+00 20.8 20.3 - 4657.212 22.1 1.243 -2.47 2.80e+00 47.8 46.4 - 4779.985 22.1 2.0477 -1.26 2.80e+00 65.8 64.9 - 4865.611 22.1 1.116 -2.81 2.80e+00 41.5 40.3 - 4874.014 22.1 3.095 -0.9 2.80e+00 37.9 36.7 - 4911.193 22.1 3.123 -0.537 2.80e+00 53.8 52.7 - 5211.54 22.1 2.59 -1.49 2.80e+00 35.1 33.5 - 5336.778 22.1 1.582 -1.630 2.80e+00 73.8 72.2 - 5381.015 22.1 1.565 -1.97 2.80e+00 62.2 60.1 - 5418.767 22.1 1.582 -2.11 2.80e+00 51.1 49.1 - 4875.486 23.0 0.040 -0.81 1.98e-32 46.1 46.6 - 5670.85 23.0 1.080 -0.42 3.58e-32 19.6 19.7 - 5727.046 23.0 1.081 -0.011 4.35e-32 40.2 40.1 - 6039.73 23.0 1.063 -0.65 3.98e-32 13.5 13.0 - 6081.44 23.0 1.051 -0.578 3.89e-32 14.2 14.4 - 6090.21 23.0 1.080 -0.062 3.98e-32 34.4 34.6 - 6119.528 23.0 1.064 -0.320 3.89e-32 21.9 21.8 - 6199.20 23.0 0.286 -1.28 1.96e-32 13.7 13.8 - 6251.82 23.0 0.286 -1.34 1.96e-32 14.4 14.9 - 4801.047 24.0 3.1216 -0.130 4.52e-32 51.0 49.3 - 4936.335 24.0 3.1128 -0.25 4.32e-32 47.1 45.4 - 5214.140 24.0 3.3694 -0.74 2.06e-32 18.2 17.6 - 5238.964 24.0 2.709 -1.27 5.19e-32 17.2 16.5 - 5247.566 24.0 0.960 -1.59 3.92e-32 83.7 82.4 - 5272.007 24.0 3.449 -0.42 3.15e-31 26.3 24.9 - 5287.20 24.0 3.438 -0.87 3.09e-31 11.7 10.8 - 5296.691 24.0 0.983 -1.36 3.92e-32 94.2 93.6 - 5300.744 24.0 0.982 -2.13 3.92e-32 61.0 60.4 - 5345.801 24.0 1.0036 -0.95 3.92e-32 115.2 113.3 - 5348.312 24.0 1.0036 -1.21 3.92e-32 101.1 100.4 - 5783.08 24.0 3.3230 -0.43 8.02e-31 33.4 32.1 - 5783.87 24.0 3.3223 -0.295 7.98e-31 46.3 44.9 - 6661.08 24.0 4.1926 -0.19 4.67e-31 13.9 13.4 - 4588.199 24.1 4.071 -0.594 2.80e+00 71.5 69.3 - 4592.049 24.1 4.073 -1.252 2.80e+00 49.0 46.6 - 5237.328 24.1 4.073 -1.087 2.80e+00 54.9 52.7 - 5246.767 24.1 3.714 -2.436 2.80e+00 16.8 15.4 - 5305.870 24.1 3.827 -1.97 2.80e+00 27.7 25.6 - 5308.41 24.1 4.0712 -1.846 2.80e+00 28.7 26.8 - 5502.067 24.1 4.168 -2.049 2.80e+00 20.1 18.3 - 5004.891 25.0 2.9197 -1.63 3.14e-32 14.3 13.9 - 5399.470 25.0 3.85 -0.104 2.80e+00 40.6 39.3 - 6013.49 25.0 3.073 -0.251 2.80e+00 87.1 86.0 - 6016.64 25.0 3.073 -0.084 2.80e+00 96.7 96.3 - 6021.79 25.0 3.076 +0.034 2.80e+00 90.4 89.6 - 5212.691 27.0 3.5144 -0.11 3.39e-31 20.6 20.9 - 5247.911 27.0 1.785 -2.08 3.27e-32 17.9 18.2 - 5301.039 27.0 1.710 -1.99 3.01e-32 19.8 20.4 - 5342.695 27.0 4.021 0.54 2.80e+00 33.9 33.7 - 5483.352 27.0 1.7104 -1.49 2.89e-32 51.2 51.5 - 5530.774 27.0 1.710 -2.23 2.26e-32 19.4 20.4 - 5647.23 27.0 2.280 -1.56 4.14e-32 14.2 14.3 - 6189.00 27.0 1.710 -2.46 2.06e-32 10.7 11.1 - 6454.995 27.0 3.6320 -0.25 3.78e-31 14.0 14.4 - 4953.208 28.0 3.7397 -0.66 3.25e-31 57.6 57.2 - 5010.938 28.0 3.635 -0.87 3.90e-31 50.7 50.1 - 5176.560 28.0 3.8982 -0.44 3.84e-31 59.9 59.0 - 5589.358 28.0 3.898 -1.14 3.98e-31 28.9 27.9 - 5643.078 28.0 4.164 -1.25 3.79e-31 15.7 15.1 - 5805.217 28.0 4.1672 -0.64 4.10e-31 44.3 42.6 - 6086.282 28.0 4.266 -0.51 4.06e-31 45.8 44.0 - 6108.116 28.0 1.676 -2.44 2.48e-32 65.7 65.3 - 6130.135 28.0 4.266 -0.96 3.91e-31 23.8 22.7 - 6176.811 28.0 4.088 -0.26 3.92e-31 65.0 64.4 - 6177.242 28.0 1.826 -3.51 2.80e+00 14.8 15.4 - 6204.604 28.0 4.088 -1.14 2.77e-31 22.7 22.5 - 6223.984 28.0 4.105 -0.98 3.93e-31 28.8 27.9 - 6378.25 28.0 4.1535 -0.90 3.91e-31 33.2 32.1 - 6643.630 28.0 1.676 -2.1 2.14e-32 94.1 93.9 - 6767.772 28.0 1.826 -2.17 2.80e+00 80.1 79.7 - 6772.315 28.0 3.657 -0.99 3.56e-31 52.0 51.2 - 7727.624 28.0 3.678 -0.4 3.43e-31 93.3 93.2 - 7797.586 28.0 3.89 -0.34 2.80e+00 80.3 78.6 - 5105.541 29.0 1.39 -1.516 2.80e+00 93.1 94.1 - 5218.197 29.0 3.816 0.476 2.80e+00 48.8 48.4 - 5220.066 29.0 3.816 -0.448 2.80e+00 17.5 17.0 - 7933.13 29.0 3.79 -0.368 2.80e+00 28.2 28.1 - 4722.159 30.0 4.03 -0.38 2.80e+00 69.9 70.0 - 4810.534 30.0 4.08 -0.16 2.80e+00 74.3 74.2 - 6362.35 30.0 5.79 0.14 2.80e+00 23.4 23.0 - 4607.338 38.0 0.00 0.283 6.55e-32 47.2 44.9 - 4854.867 39.1 0.9923 -0.38 2.80e+00 57.7 54.4 - 4883.685 39.1 1.0841 0.07 2.80e+00 62.0 59.3 - 4900.110 39.1 1.0326 -0.09 2.80e+00 58.7 55.9 - 5087.420 39.1 1.0841 -0.17 2.80e+00 51.1 47.7 - 5200.413 39.1 0.9923 -0.57 2.80e+00 42.0 38.4 - 4050.320 40.1 0.713 -1.06 2.80e+00 26.2 23.9 - 4208.980 40.1 0.713 -0.51 2.80e+00 45.5 42.1 - 4442.992 40.1 1.486 -0.42 2.80e+00 26.8 24.4 - 3158.16 42.0 0.000 -0.31 2.80e+00 13.3 11.6 - 3193.97 42.0 0.000 0.07 2.80e+00 16.2 15.0 - 3436.736 44.0 0.148 0.165 2.80e+00 8.8 7.0 - 3498.942 44.0 0.000 0.322 2.80e+00 27.4 25.4 - 3242.698 46.0 0.8138 0.07 2.80e+00 27.2 25.5 - 3404.576 46.0 0.8138 0.33 2.80e+00 32.9 30.2 - 3516.944 46.0 0.9615 -0.21 2.80e+00 16.2 14.1 - 3280.68 47.0 0.0000 -0.022 2.80e+00 40.6 38.5 - 3382.89 47.0 0.0000 -0.334 2.80e+00 29.3 26.9 - 5853.67 56.1 0.604 -0.91 5.30e-32 67.7 63.4 - 6141.71 56.1 0.704 -0.08 5.30e-32 121.2 114.9 - 6496.90 56.1 0.604 -0.38 5.30e-32 106.8 99.9 - 4662.50 57.1 0.0000 -1.24 2.80e+00 7.6 6.1 - 4748.73 57.1 0.9265 -0.54 2.80e+00 4.8 3.9 - 5303.53 57.1 0.3213 -1.35 2.80e+00 4.1 3.2 - 3942.151 58.1 0.000 -0.22 2.80e+00 13.4 11.9 - 3999.237 58.1 0.295 0.06 2.80e+00 20.4 17.2 - 4042.581 58.1 0.495 0.00 2.80e+00 12.1 10.0 - 4073.474 58.1 0.477 0.21 2.80e+00 16.4 14.4 - 4364.653 58.1 0.495 -0.17 2.80e+00 13.6 11.4 - 4523.075 58.1 0.516 -0.08 2.80e+00 15.0 12.9 - 4562.359 58.1 0.477 0.21 2.80e+00 27.4 24.7 - 5274.229 58.1 1.044 0.13 2.80e+00 9.8 8.5 - 5259.73 59.1 0.633 0.114 2.80e+00 2.9 2.3 - 4021.33 60.1 0.320 -0.10 2.80e+00 14.8 12.4 - 4059.95 60.1 0.204 -0.52 2.80e+00 7.5 5.4 - 4446.38 60.1 0.204 -0.35 2.80e+00 11.9 9.7 - 5234.19 60.1 0.550 -0.51 2.80e+00 6.7 5.2 - 5293.16 60.1 0.822 0.10 2.80e+00 12.5 10.0 - 5319.81 60.1 0.550 -0.14 2.80e+00 13.3 10.8 - 3760.70 62.1 0.18 -0.42 2.80e+00 9.9 7.1 - 4318.936 62.1 0.277 -0.25 2.80e+00 14.1 10.4 - 4467.341 62.1 0.659 0.15 2.80e+00 14.5 12.6 - 4519.630 62.1 0.543 -0.35 2.80e+00 6.7 5.3 - 4676.902 62.1 0.040 -0.87 2.80e+00 6.6 5.0 - 3819.67 63.1 0.000 0.51 2.80e+00 47.4 33.2 - 3907.11 63.1 0.207 0.17 2.80e+00 28.8 19.0 - 4129.72 63.1 0.000 0.22 2.80e+00 70.1 56.3 - 6645.10 63.1 1.379 0.12 2.80e+00 8.0 6.6 - 3331.387 64.1 0.000 -0.28 2.80e+00 10.2 7.8 - 3697.733 64.1 0.032 -0.34 2.80e+00 6.0 4.5 - 3712.704 64.1 0.382 0.04 2.80e+00 15.5 12.3 - 3768.396 64.1 0.078 0.21 2.80e+00 17.7 14.0 - 4251.731 64.1 0.382 -0.22 2.80e+00 13.3 10.4 - 3531.71 66.1 0.000 0.77 2.80e+00 40.3 36.1 - 3536.02 66.1 0.538 0.53 2.80e+00 24.4 19.8 - 3694.81 66.1 0.103 -0.11 2.80e+00 17.3 13.9 - 4077.97 66.1 0.103 -0.04 2.80e+00 14.8 10.7 - 4103.31 66.1 0.103 -0.38 2.80e+00 17.8 15.0 - 4449.70 66.1 0.000 -1.03 2.80e+00 4.3 3.3 - 3694.19 70.1 0.000 -0.30 2.80e+00 81.2 73.9 diff --git a/src/fit.jl b/src/fit.jl index b9804fc8..f2ed84d7 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -436,6 +436,11 @@ end """ ews_to_abundances(atm, linelist, A_X, measured_EWs; kwargs... ) +!!! danger + + Functions using equivalent widths have been dissabled while we + [fix some bugs](https://github.com/ajwheeler/Korg.jl/pull/331). + Compute per-line abundances on the linear part of the curve of growth given a model atmosphere and a list of lines with equivalent widths. @@ -466,6 +471,7 @@ A vector of abundances (`A(X) = log10(n_X/n_H) + 12` format) for each line in `l """ function ews_to_abundances(atm, linelist, A_X, measured_EWs; ew_window_size::Real=2.0, wl_step=0.01, blend_warn_threshold=0.01, synthesize_kwargs...) + throw("This function is currently disabled while we fix some bugs. See https://github.com/ajwheeler/Korg.jl/pull/331 for details.") synthesize_kwargs = Dict(synthesize_kwargs) if length(linelist) != length(measured_EWs) throw(ArgumentError("length of linelist does not match length of ews ($(length(linelist)) != $(length(measured_EWs)))")) @@ -537,6 +543,11 @@ end """ ews_to_stellar_parameters(linelist, measured_EWs, [measured_EW_err]; kwargs...) +!!! danger + + Functions using equivalent widths have been dissabled while we + [fix some bugs](https://github.com/ajwheeler/Korg.jl/pull/331). + Find stellar parameters from equivalent widths the "old fashioned" way. This function finds the values of ``T_\\mathrm{eff}``, ``\\log g``, ``v_{mic}``, and [m/H] which satisfy the following conditions (using a Newton-Raphson solver): @@ -607,6 +618,7 @@ function ews_to_stellar_parameters(linelist, measured_EWs, Korg._sdss_marcs_atmospheres[1][3][end])], fix_params=[false, false, false, false], callback=Returns(nothing), max_iterations=30, passed_kwargs...) + throw("This function is currently disabled while we fix some bugs. See https://github.com/ajwheeler/Korg.jl/pull/331 for details.") if :vmic in keys(passed_kwargs) throw(ArgumentError("vmic must not be specified, because it is a parameter fit by ews_to_stellar_parameters. Did you mean to specify vmic0, the starting value? See the documentation for ews_to_stellar_parameters if you would like to fix microturbulence to a given value.")) end diff --git a/test/fit.jl b/test/fit.jl index d34a2fae..678f5267 100644 --- a/test/fit.jl +++ b/test/fit.jl @@ -123,174 +123,179 @@ using Random 1:end-1]) end - @testset "ews_to_abundances" begin - @testset "require sorted linelists" begin - sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) - sun_A_X = Korg.format_A_X(sun_Fe_H) - sun_atm = Korg.read_model_atmosphere("data/sun.mod") - - linelist = [ - Korg.Line(5054.642 * 1e-8, -1.92100, Korg.Species("26.0"), 3.64, 4.68e-32), - Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, 2.71e-31) - ] - sun_ews = [74.3, 40.5] - @test_throws ArgumentError Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, - sun_ews, vmic=sun_vmic) - end - - @testset "length of linelist and ews should match" begin - sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) - sun_A_X = Korg.format_A_X(sun_Fe_H) - sun_atm = Korg.read_model_atmosphere("data/sun.mod") - linelist = [Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, - 2.71e-31)] - @test_throws ArgumentError Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, [], - vmic=sun_vmic) - end + if false # SKIP EWs tests until we fix them + @testset "ews_to_abundances" begin + @testset "require sorted linelists" begin + sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) + sun_A_X = Korg.format_A_X(sun_Fe_H) + sun_atm = Korg.read_model_atmosphere("data/sun.mod") + + linelist = [ + Korg.Line(5054.642 * 1e-8, -1.92100, Korg.Species("26.0"), 3.64, 4.68e-32), + Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, 2.71e-31) + ] + sun_ews = [74.3, 40.5] + @test_throws ArgumentError Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, + sun_ews, vmic=sun_vmic) + end - @testset "disallow molecules" begin - sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) - sun_A_X = Korg.format_A_X(sun_Fe_H) - sun_atm = Korg.read_model_atmosphere("data/sun.mod") - linelist = [ - Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, 2.71e-31), - Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("106.0"), 2.8512, 2.71e-31)] - sun_ews = [74.3, 40.5] - @test_throws ArgumentError Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, - sun_ews, vmic=sun_vmic) - end + @testset "length of linelist and ews should match" begin + sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) + sun_A_X = Korg.format_A_X(sun_Fe_H) + sun_atm = Korg.read_model_atmosphere("data/sun.mod") + linelist = [Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, + 2.71e-31)] + @test_throws ArgumentError Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, + [], + vmic=sun_vmic) + end - @testset "Melendez et al. (2014) sanity check" begin - sun_ews = [74.3, 40.5, 96.1, 19.1, 80.7] - sco_ews = [74.8, 40.9, 97.5, 18.9, 84.0] - linelist = [ - Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, 2.71e-31), - Korg.Line(5054.642 * 1e-8, -1.92100, Korg.Species("26.0"), 3.64, 4.68e-32), - Korg.Line(5127.359 * 1e-8, -3.30700, Korg.Species("26.0"), 0.915, 1.84e-32), - # don't convert these ones to cm. It should still work. - Korg.Line(5127.679, -6.12500, Korg.Species("26.0"), 0.052, 1.2e-32), - Korg.Line(5197.577, -2.22000, Korg.Species("26.1"), 3.2306, 8.69e-33) - ] - sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) - sun_A_X = Korg.format_A_X(sun_Fe_H) - sun_atm = Korg.read_model_atmosphere("data/sun.mod") - - sco_teff, sco_logg, sco_fe_h, sco_vmic = (5823, 4.45, 0.054, sun_vmic + 0.02) - sco_A_X = Korg.format_A_X(sco_fe_h) - # Note: NOT true, but just for testing the whole performance - sco_atm = sun_atm - - sun_abundances = Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, sun_ews; - vmic=sun_vmic) - sco_abundances = Korg.Fit.ews_to_abundances(sco_atm, linelist, sco_A_X, sco_ews; - vmic=sco_vmic) - diff_abundances = sco_abundances .- sun_abundances - - mean_diff_abundances = sum(diff_abundances) / length(diff_abundances) - @test abs(mean_diff_abundances) < 0.05 - # TODO: test for stddev? - end - end + @testset "disallow molecules" begin + sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) + sun_A_X = Korg.format_A_X(sun_Fe_H) + sun_atm = Korg.read_model_atmosphere("data/sun.mod") + linelist = [ + Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, 2.71e-31), + Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("106.0"), 2.8512, 2.71e-31)] + sun_ews = [74.3, 40.5] + @test_throws ArgumentError Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, + sun_ews, vmic=sun_vmic) + end - @testset "stellar parameters via EWs" begin - # used in both tests below - good_linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 4.1), - Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 4.2), - Korg.Line(7e-5, -1.92100, Korg.species"Fe I", 4.3), - Korg.Line(8e-5, -1.92100, Korg.species"Fe II", 4.4)] - - @testset "validate arguments" begin - # vmic can't be specified - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, - ones(length(good_linelist)); - vmic=1.0) - - # number of EWs, EW uncertainties, and lines should match - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, [1.0]) - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, - ones(length(good_linelist)), - [1.0]) - - # different elements - linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 0, 0), - Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 0, 0), - Korg.Line(7e-5, -1.92100, Korg.species"Mn I", 0, 0), - Korg.Line(8e-5, -1.92100, Korg.species"Fe II", 0, 0)] - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, - ones(length(linelist))) - - # can't use a molecule - linelist = [Korg.Line(5e-5, -2.05800, Korg.species"CO I", 0, 0), - Korg.Line(6e-5, -1.92100, Korg.species"CO I", 0, 0), - Korg.Line(7e-5, -1.92100, Korg.species"CO I", 0, 0), - Korg.Line(8e-5, -1.92100, Korg.species"CO II", 0, 0)] - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, - ones(length(linelist))) - - # no ions - linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 0, 0), - Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 0, 0), - Korg.Line(6.5e-5, -1.92100, Korg.species"Fe I", 0, 0), - Korg.Line(7e-5, -1.92100, Korg.species"Fe I", 0, 0)] - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, - ones(length(linelist))) - - # not enough lines - linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 0, 0), - Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 0, 0), - Korg.Line(7e-5, -1.92100, Korg.species"Fe II", 0, 0)] - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, - ones(length(linelist))) - - # parameter ranges must have positive measure - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, - ones(length(good_linelist)); - parameter_ranges=[(5777, - 5777), - (3.0, 4.0), - (1.0, 2.0), - (-0.5, 0.0)]) - - # vmic0 can't be 0 - @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, - ones(length(good_linelist)); - vmic0=0) + @testset "Melendez et al. (2014) sanity check" begin + sun_ews = [74.3, 40.5, 96.1, 19.1, 80.7] + sco_ews = [74.8, 40.9, 97.5, 18.9, 84.0] + linelist = [ + Korg.Line(5044.211 * 1e-8, -2.05800, Korg.Species("26.0"), 2.8512, 2.71e-31), + Korg.Line(5054.642 * 1e-8, -1.92100, Korg.Species("26.0"), 3.64, 4.68e-32), + Korg.Line(5127.359 * 1e-8, -3.30700, Korg.Species("26.0"), 0.915, 1.84e-32), + # don't convert these ones to cm. It should still work. + Korg.Line(5127.679, -6.12500, Korg.Species("26.0"), 0.052, 1.2e-32), + Korg.Line(5197.577, -2.22000, Korg.Species("26.1"), 3.2306, 8.69e-33) + ] + sun_Teff, sun_logg, sun_Fe_H, sun_vmic = (5777, 4.44, 0.0, 1.0) + sun_A_X = Korg.format_A_X(sun_Fe_H) + sun_atm = Korg.read_model_atmosphere("data/sun.mod") + + sco_teff, sco_logg, sco_fe_h, sco_vmic = (5823, 4.45, 0.054, sun_vmic + 0.02) + sco_A_X = Korg.format_A_X(sco_fe_h) + # Note: NOT true, but just for testing the whole performance + sco_atm = sun_atm + + sun_abundances = Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, sun_ews; + vmic=sun_vmic) + sco_abundances = Korg.Fit.ews_to_abundances(sco_atm, linelist, sco_A_X, sco_ews; + vmic=sco_vmic) + diff_abundances = sco_abundances .- sun_abundances + + mean_diff_abundances = sum(diff_abundances) / length(diff_abundances) + @test abs(mean_diff_abundances) < 0.05 + # TODO: test for stddev? + end end - @testset "basic fit" begin - # 2 Å wide window around each line - synth_wls = map(good_linelist) do line - wl = line.wl * 1e8 - wl-1.0:0.01:wl+1.0 - end - sol = synthesize(interpolate_marcs(5777.0, 4.44), good_linelist, format_A_X(), - synth_wls; hydrogen_lines=false) - - # EWs you get for the fake linelist with solar params - # the real implementation uses the trapezoid rule, but this is close enough - EWs = [sum((1 .- sol.flux./sol.cntm)[r]) for r in sol.subspectra] * 10 #0.01 Å -> mÅ - EW_err = ones(length(EWs)) * 0.5 - best_fit_params, stat_err, sys_err = Korg.Fit.ews_to_stellar_parameters(good_linelist, - EWs, EW_err) - - @test sys_err == [0.0, 0.0, 0.0, 0.0] - for (i, p) in enumerate([5777.0, 4.44, 1.0, 0.0]) - @test best_fit_params[i]≈p atol=stat_err[i] + @testset "stellar parameters via EWs" begin + # used in both tests below + good_linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 4.1), + Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 4.2), + Korg.Line(7e-5, -1.92100, Korg.species"Fe I", 4.3), + Korg.Line(8e-5, -1.92100, Korg.species"Fe II", 4.4)] + + @testset "validate arguments" begin + # vmic can't be specified + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, + ones(length(good_linelist)); + vmic=1.0) + + # number of EWs, EW uncertainties, and lines should match + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, [1.0]) + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, + ones(length(good_linelist)), + [1.0]) + + # different elements + linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 0, 0), + Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 0, 0), + Korg.Line(7e-5, -1.92100, Korg.species"Mn I", 0, 0), + Korg.Line(8e-5, -1.92100, Korg.species"Fe II", 0, 0)] + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, + ones(length(linelist))) + + # can't use a molecule + linelist = [Korg.Line(5e-5, -2.05800, Korg.species"CO I", 0, 0), + Korg.Line(6e-5, -1.92100, Korg.species"CO I", 0, 0), + Korg.Line(7e-5, -1.92100, Korg.species"CO I", 0, 0), + Korg.Line(8e-5, -1.92100, Korg.species"CO II", 0, 0)] + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, + ones(length(linelist))) + + # no ions + linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 0, 0), + Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 0, 0), + Korg.Line(6.5e-5, -1.92100, Korg.species"Fe I", 0, 0), + Korg.Line(7e-5, -1.92100, Korg.species"Fe I", 0, 0)] + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, + ones(length(linelist))) + + # not enough lines + linelist = [Korg.Line(5e-5, -2.05800, Korg.species"Fe I", 0, 0), + Korg.Line(6e-5, -1.92100, Korg.species"Fe I", 0, 0), + Korg.Line(7e-5, -1.92100, Korg.species"Fe II", 0, 0)] + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(linelist, + ones(length(linelist))) + + # parameter ranges must have positive measure + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, + ones(length(good_linelist)); + parameter_ranges=[(5777, + 5777), + (3.0, 4.0), + (1.0, 2.0), + (-0.5, 0.0)]) + + # vmic0 can't be 0 + @test_throws ArgumentError Korg.Fit.ews_to_stellar_parameters(good_linelist, + ones(length(good_linelist)); + vmic0=0) end - # check that fixing parameters works - for (i, param) in enumerate([:Teff0, :logg0, :vmic0, :m_H0]) - fixed_params = zeros(Bool, 4) - fixed_params[i] = true - kwargs = Dict(param => best_fit_params[i]) - # start teff and logg close to right answer to make it faster. - bestfit_fixed, _, _ = Korg.Fit.ews_to_stellar_parameters(good_linelist, EWs, EW_err; - Teff0=5740.0, logg0=4.4, - fix_params=fixed_params, - kwargs...) - - for i in 1:4 - @test bestfit_fixed[i]≈best_fit_params[i] atol=stat_err[i]/10 + @testset "basic EW fit" begin + # 2 Å wide window around each line + synth_wls = map(good_linelist) do line + wl = line.wl * 1e8 + wl-1.0:0.01:wl+1.0 + end + sol = synthesize(interpolate_marcs(5777.0, 4.44), good_linelist, format_A_X(), + synth_wls; hydrogen_lines=false) + + # EWs you get for the fake linelist with solar params + # the real implementation uses the trapezoid rule, but this is close enough + EWs = [sum((1 .- sol.flux./sol.cntm)[r]) for r in sol.subspectra] * 10 #0.01 Å -> mÅ + EW_err = ones(length(EWs)) * 0.5 + best_fit_params, stat_err, sys_err = Korg.Fit.ews_to_stellar_parameters(good_linelist, + EWs, EW_err) + + @test sys_err == [0.0, 0.0, 0.0, 0.0] + for (i, p) in enumerate([5777.0, 4.44, 1.0, 0.0]) + @test best_fit_params[i]≈p atol=stat_err[i] + end + + # check that fixing parameters works + for (i, param) in enumerate([:Teff0, :logg0, :vmic0, :m_H0]) + fixed_params = zeros(Bool, 4) + fixed_params[i] = true + kwargs = Dict(param => best_fit_params[i]) + # start teff and logg close to right answer to make it faster. + bestfit_fixed, _, _ = Korg.Fit.ews_to_stellar_parameters(good_linelist, EWs, + EW_err; + Teff0=5740.0, + logg0=4.4, + fix_params=fixed_params, + kwargs...) + + for i in 1:4 + @test bestfit_fixed[i]≈best_fit_params[i] atol=stat_err[i]/10 + end end end end