{ "cells": [ { "cell_type": "markdown", "id": "41ae9ef8-a389-4246-a783-9eb9f103cfe9", "metadata": {}, "source": [ "# Example: Fit Helium" ] }, { "cell_type": "code", "execution_count": null, "id": "2c2fa30e-6bfd-47bc-9672-3245cb376f59", "metadata": {}, "outputs": [], "source": [ "import sunbather\n", "import sunbather.tools\n", "\n", "sunbather.tools.get_sunbather_project_path()\n", "sunbather.firstrun(quiet=True)" ] }, { "cell_type": "markdown", "id": "d7370e97", "metadata": {}, "source": [ "## License\n", "The code in this notebook is free to be used, edited and redistributed by anyone free of charge. Please cite Linssen et al. (2024) when making use of _sunbather_ and/or the code in this notebook." ] }, { "cell_type": "markdown", "id": "90d0f8d2", "metadata": { "jp-MarkdownHeadingCollapsed": true }, "source": [ "## Scientific goal\n", "In this notebook, we demonstrate how sunbather can be used to obtain mass-loss rate estimates from spectral observations. Fitting Parker wind models to helium data is a commonly applied method, but it can result in a degeneracy between the free parameters of the Parker wind model: the temperature and the mass-loss rate. Constraints on the mass-loss rate are not very stringent in that case. In this notebook, we follow the approach described in Linssen et al. (2022), to mitigate this. We fit the helium line and combine this with additional constraints on the Parker wind temperature parameter, in order to break the degeneracy and get better constraints on the mass-loss rate. We first do this for a simple case of a spectrum that we generate from random data, which will reproduce the results of Fig. 6 from Linssen et al. (2024). Then, we move on to fitting the observed He 10830 Å signal of TOI-2134 b from Zhang et al. (2023). We will perform fits with both an atmospheric composition that is a hydrogen/helium mixture in the solar ratio, and a 100x solar metallicity atmosphere. This reproduces the results of Fig. 7 from Linssen et al. (2024)." ] }, { "cell_type": "markdown", "id": "9c3b972d", "metadata": {}, "source": [ "## Methodology\n", "\n", "We need to run a grid of Parker wind models with different temperatures ($T_0$) and mass-loss rates ($\\dot{M}$), which can be computationally expensive depending on the size of the grid. Here, we include the commands to run your own grid (as commented out Python lines), but by default this notebook uses the model runs that we have pre-run to save computation time.\n", "\n", "One analysis typically consists of five main steps, which are the same as those of *sunbather*'s \"fit_helium\" example notebook:\n", "\n", "- Step 1. Generate the $T_0$-$\\dot{M}$ grid of isothermal (!) Parker wind models.\n", "- Step 2. Run the Parker wind models through Cloudy to find a non-isothermal temperature profile and the metastable helium density.\n", "- Step 3. Make synthetic helium transit spectra for each model and fit those to the observed value (to obtain the _likelihood_).\n", "- Step 4. Constrain the temperature parameter-space by assessing the model self-consistency (to obtain the _prior_).\n", "- Step 5. Combine the data fit with the model self-consistency (to obtain the _posterior_).\n", "\n", "Steps 4 and 5 can optionally also be skipped in the case of spectrally resolved observations, since then the $T_0-\\dot{M}$-degeneracy is usually less strong, and no prior is needed to obtain good constraints.\n", "\n", "Some of the steps we need to take to go through this analysis have to be taken outside of *sunbather* (e.g. storing the planet parameters in the planets.txt file and making the stellar SEDs available to *Cloudy*). Some steps are part of the \"core functionality\" of the *sunbather* package, which means they can be executed with a single command-line call to a Python script (e.g. generating the Parker wind profiles, running them through the Cloudy algorithm). And some of the steps are performed inside a user-made Python file/Jupyter notebook (e.g. saving the helium spectrum of each Parker wind profile, doing the Bayesian analysis to extract the mass-loss rate). In this notebook, we will provide all neccessary code." ] }, { "cell_type": "markdown", "id": "3e620f6b-9d14-4796-bef8-adf0ced1ca22", "metadata": {}, "source": [ "## Preparation\n", "\n", "To reproduce these results, we assume you have all codes set-up. That is; you have downloaded *sunbather* and installed its dependencies (the Python packages, including *p-winds*). You have installed *Cloudy* and have the path to it stored as the *\\$CLOUDY_PATH* environmental variable. You have created your \"project\" folder, and have the path to it stored as the *\\$SUNBATHER_PROJECT_PATH* environmental variable. These steps are described in more detail in the \"installation\" section of the _sunbather_ wiki.\n", "\n", "Before *sunbather* can create Parker wind profiles, we need to make sure the parameters of the system are available to the code. The parameters are stored in the *$SUNBATHER_PROJECT_PATH/planets.txt* file. To run the models for the generic hot Neptune mock retrieval, as well as the TOI-2134 b models, you will need to add their parameters.\n", "\n", "**You can manually add the parameters or execute the following cell:**" ] }, { "cell_type": "code", "execution_count": null, "id": "29bd070a-1891-4325-a53a-7d4467175935", "metadata": {}, "outputs": [], "source": [ "# Adds two lines to _$SUNBATHER_PROJECT_PATH/planets.txt_:\n", "with open(\n", " f\"{sunbather.tools.get_sunbather_project_path()}/planets.txt\", \"r+\", encoding=\"utf-8\"\n", ") as planetfile:\n", " planets = planetfile.read()\n", " if \"TOI2134b\" not in planets:\n", " planetfile.write(\n", " \"\\n\"\n", " \"TOI2134b,TOI-2134 b,0.240,0.709,0.078,0.0287,0.744,0.20,TOI2134.spec \"\n", " \"#parameters from Zhang et al. (2023)\\n\"\n", " )\n", " if \"hotNeptune\" not in planets:\n", " planetfile.write(\n", " \"hotNeptune,generic hot Neptune,0.5,1.0,0.05,0.1,1.0,0.,solar.spec\\n\"\n", " )" ] }, { "cell_type": "markdown", "id": "0aa89533-96e2-4b80-b884-a64f215815d0", "metadata": {}, "source": [ "The last column of the _planets.txt_ file specifies the **name** of the stellar SED that we want to use. The SED with exactly this name must be available to *Cloudy*, so it must be placed in its source folder, specifically: _\\$CLOUDY_PATH/data/SED/_. In the */sunbather/stellar_SEDs/* folder, the solar.spec file is provided, which is combined from Woods et al. (2005) and Rottman (2005). The TOI2134.spec file is also provided, which is downloaded from the online material of Zhang et al. (2023) ( https://iopscience.iop.org/article/10.3847/2041-8213/aced51 ) and put in the format accepted by *sunbather* (see wiki).\n", "\n", "**You can manually copy the SEDs or execute the following cell:**" ] }, { "cell_type": "code", "execution_count": null, "id": "ede7a4e4-8bb3-47b4-b18e-76f97acd3064", "metadata": {}, "outputs": [], "source": [ "import os\n", "cloudy_sed_folder = f\"{sunbather.tools.get_cloudy_path()}/data/SED\"\n", "local_sed_folder = f\"{sunbather.tools.get_sunbather_project_path}/stellar_SEDs/\"\n", "\n", "for file in [\"TOI2134.spec\", \"solar.spec\"]:\n", " if not os.path.isfile(f\"{cloudy_sed_folder}/{file}\"):\n", " shutil.copyfile(\n", " f\"{local_sed_folder}/{file}\",\n", " f\"{cloudy_sed_folder}/{file}\"\n", " )" ] }, { "cell_type": "markdown", "id": "84e276fe", "metadata": {}, "source": [ "## Function definitions\n", "Since this notebook contains quite a lot of functions for performing the various calculations, we have grouped them all together here, so that the rest of the notebook has improved readability. Once the functions are defined, most parts of the analysis are really only one line of code!" ] }, { "cell_type": "code", "execution_count": null, "id": "4d7362e7", "metadata": {}, "outputs": [], "source": [ "import os.path\n", "import shutil\n", "\n", "import traceback\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.integrate import (trapezoid, cumulative_trapezoid)\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "import scipy.stats as sps\n", "from scipy.optimize import curve_fit\n", "import spectres # to resample synthetic spectra to the wavelength grid of observations (pip install spectres)\n", "\n", "# import sunbather modules\n", "from sunbather import construct_parker\n", "from sunbather import tools\n", "from sunbather import RT\n", "\n", "# for interactive matplotlib plots\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "id": "1084514e-cde2-4120-85d8-5075df89e2f1", "metadata": {}, "outputs": [], "source": [ "projectpath = tools.get_sunbather_project_path()" ] }, { "cell_type": "code", "execution_count": null, "id": "7dd9f6cd", "metadata": {}, "outputs": [], "source": [ "# Useful functions\n", "\n", "\n", "def get_xyz(df):\n", " \"\"\"\n", " Takes as input a Pandas DataFrame which has the x-values\n", " as the index, and the y-values as the column names,\n", " and the z-values as the dataframe contents.\n", " Returns three 1D arrays with the x, y, z values\n", " that can e.g. be plotted with matplotlib tricontour.\n", " \"\"\"\n", "\n", " stacked_data = df.stack()\n", " stacked_data = stacked_data.dropna()\n", " x_values = stacked_data.index.get_level_values(0).tolist() # x-values\n", " y_values = stacked_data.index.get_level_values(1).tolist() # y-values\n", " z_values = stacked_data.tolist() # z-values\n", "\n", " return x_values, y_values, z_values\n", "\n", "\n", "def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=1000):\n", " \"\"\"\n", " Cuts a colormap to a smaller range.\n", " from https://stackoverflow.com/questions/18926031/how-to-extract-a-subset-of-a-colormap-as-a-new-colormap-in-matplotlib\n", " \"\"\"\n", "\n", " new_cmap = matplotlib.colors.LinearSegmentedColormap.from_list(\n", " \"trunc({n},{a:.2f},{b:.2f})\".format(n=cmap.name, a=minval, b=maxval),\n", " cmap(np.linspace(minval, maxval, n)),\n", " )\n", " return new_cmap\n", "\n", "\n", "def get_array_from_contour(cn, levelnum):\n", " \"\"\"\n", " Takes an ax.contour object and returns the 'levelnum'-th contour as an\n", " array of x and y values. You can then extract the min/max values from\n", " such a contour/array.\n", " \"\"\"\n", "\n", " lines = []\n", " for line in cn.collections[levelnum].get_paths():\n", " lines.append(line.vertices)\n", "\n", " if lines:\n", " x = np.concatenate(lines, axis=0)[:, 0]\n", " y = np.concatenate(lines, axis=0)[:, 1]\n", " else:\n", " print(\"No 1 sigma contour could be found\")\n", " x, y = None, None\n", "\n", " return x, y\n", "\n", "\n", "# Functions that calculate statistics\n", "\n", "\n", "def calc_chisqs_He10830(\n", " Tstrucpath,\n", " observed_wavs,\n", " observed_ea,\n", " observed_sig_ea,\n", " instrument_R=None,\n", " T0low=2000,\n", " T0up=12000,\n", " T0step=100,\n", " Mdotlow=8,\n", " Mdotup=12,\n", " Mdotstep=0.05,\n", " **kwargs,\n", "):\n", " \"\"\"\n", " Reads in all Parker wind profiles present in the Tstrucpath folder,\n", " calculates the metastable helium line, compares it to the observed\n", " profile and calculates the chi-squared value.\n", " Expects excess absorption and error in units of %.\n", " \"\"\"\n", "\n", " Mdots = [\n", " \"%.3f\" % Mdot for Mdot in np.arange(Mdotlow, Mdotup + 1e-5, Mdotstep)\n", " ] # set up a grid of Mdot\n", " T0s = [\n", " \"%i\" % T0 for T0 in np.arange(T0low, T0up + 1e-5, T0step)\n", " ] # set up a grid of T_0\n", "\n", " chisqs = pd.DataFrame(columns=Mdots, index=T0s, dtype=float)\n", "\n", " for Mdot in Mdots:\n", " for T0 in T0s:\n", " try:\n", " # read in the converged simulation for this combination of T0 and Mdot by specifying the path\n", " sim = tools.Sim(Tstrucpath + \"parker_\" + T0 + \"_\" + Mdot + \"/converged\")\n", " # we calculate the model spectrum on a high-resolution grid\n", " highres_wavs = np.linspace(10830, 10836, 100)\n", " FinFout, lines_found, lines_not_found = RT.FinFout(\n", " sim, highres_wavs, \"He\", **kwargs\n", " )\n", " # convert to excess absorption in units of %\n", " model_excess_absorption = (np.max(FinFout) - FinFout) * 100\n", " # convolve with instrumental resolution, if any was passed\n", " if instrument_R is not None:\n", " model_excess_absorption = RT.convolve_spectrum_R(\n", " highres_wavs, model_excess_absorption, instrument_R\n", " ) # careful - it gets the same variable name\n", " # resample model to the observed wav points\n", " model_excess_absorption_ondata = spectres.spectres(\n", " observed_wavs,\n", " highres_wavs,\n", " model_excess_absorption,\n", " fill=0.0,\n", " verbose=False,\n", " )\n", " # calculate chi squared\n", " chisq = np.sum(\n", " ((model_excess_absorption_ondata - observed_ea) / observed_sig_ea)\n", " ** 2\n", " )\n", " # save chisq value in the pandas dataframe\n", " chisqs.loc[Mdot][T0] = chisq\n", "\n", " except (\n", " FileNotFoundError\n", " ) as e: # then this Parker wind model was not calculated\n", " pass\n", "\n", " except Exception as e: # if something else went wrong\n", " traceback.print_exc()\n", "\n", " return chisqs\n", "\n", "\n", "def calc_EWs_He10830(\n", " Tstrucpath,\n", " T0low=2000,\n", " T0up=12000,\n", " T0step=100,\n", " Mdotlow=8,\n", " Mdotup=12,\n", " Mdotstep=0.05,\n", " **kwargs,\n", "):\n", " \"\"\"\n", " Reads in all Parker wind profiles present in the Tstrucpath folder,\n", " calculates the metastable helium line and integrates it to get the EW.\n", " \"\"\"\n", "\n", " Mdots = [\n", " \"%.3f\" % Mdot for Mdot in np.arange(Mdotlow, Mdotup + 1e-5, Mdotstep)\n", " ] # set up a grid of Mdot\n", " T0s = [\n", " \"%i\" % T0 for T0 in np.arange(T0low, T0up + 1e-5, T0step)\n", " ] # set up a grid of T_0\n", "\n", " model_EWs = pd.DataFrame(columns=Mdots, index=T0s, dtype=float)\n", "\n", " for Mdot in Mdots:\n", " for T0 in T0s:\n", " try:\n", " # read in the converged simulation for this combination of T0 and Mdot by specifying the path\n", " sim = tools.Sim(Tstrucpath + \"parker_\" + T0 + \"_\" + Mdot + \"/converged\")\n", " # set up the wavelength array in vacuum angstrom units\n", " wavs = np.logspace(np.log10(10831), np.log10(10835), num=100)\n", " # run the radiative transfer (to check if all lines were calculated, print lines_found)\n", " FinFout, lines_found, lines_not_found = RT.FinFout(\n", " sim, wavs, \"He\", **kwargs\n", " )\n", " # convert from Fin/Fout to excess absorption - this assumes we reach the continuum somewhere\n", " absorption = np.max(FinFout) - FinFout\n", " # integrate to get EW\n", " EW = trapezoid(absorption, x=wavs)\n", " # save EW value in the pandas dataframe\n", " model_EWs.loc[Mdot][T0] = EW\n", "\n", " except (\n", " FileNotFoundError\n", " ) as e: # then this Parker wind model was not calculated\n", " pass\n", "\n", " except Exception as e: # if something else went wrong\n", " traceback.print_exc()\n", "\n", " return model_EWs\n", "\n", "\n", "def metaHe_weighted_T(sim):\n", " \"\"\"\n", " Calculates the mean temperature and its standard devation of the atmosphere weighted by\n", " the metastable helium number density (Eq. 4 & 5 in Linssen et al. 2022).\n", " \"\"\"\n", "\n", " # in Cloudy, not all bins have the same thickness, so we must compensate for that to do a fair weighting\n", " bin_thickness = np.diff(np.insert(sim.den.depth.values, 0, 0.0))\n", " # the metastable helium state is the second lowest energy level and thus stored in the He[2] column of the density file\n", " T_He = np.sum(sim.ovr.Te.values * sim.den[\"He[2]\"].values * bin_thickness) / np.sum(\n", " sim.den[\"He[2]\"].values * bin_thickness\n", " )\n", " sigmaT = np.sqrt(\n", " np.sum(\n", " sim.den[\"He[2]\"].values * bin_thickness * (sim.ovr.Te.values - T_He) ** 2\n", " )\n", " / np.sum(sim.den[\"He[2]\"].values * bin_thickness)\n", " )\n", "\n", " return T_He, sigmaT\n", "\n", "\n", "def calc_dT_helium(\n", " Tstrucpath, T0low=2000, T0up=12000, T0step=100, Mdotlow=8, Mdotup=12, Mdotstep=0.05\n", "):\n", " \"\"\"\n", " Reads in all Parker wind profiles present in the Tstrucpath folder,\n", " calculates the mean temperature weighted by the metastable helium density\n", " and its standard deviation. Calculates the temperature difference between\n", " the isothermal value and the He-weighted value, which can be used as\n", " a measure of model self-consistency.\n", " \"\"\"\n", "\n", " Mdots = [\n", " \"%.3f\" % Mdot for Mdot in np.arange(Mdotlow, Mdotup + 1e-5, Mdotstep)\n", " ] # set up a grid of Mdot\n", " T0s = [\n", " \"%i\" % T0 for T0 in np.arange(T0low, T0up + 1e-5, T0step)\n", " ] # set up a grid of T_0\n", "\n", " T_He = pd.DataFrame(columns=Mdots, index=T0s, dtype=float) # stores the weighted T\n", " sigmaT = T_He.copy() # stores standard deviation of T(r) around T_He\n", " dT = T_He.copy() # stores T_He - isothermal T\n", "\n", " for Mdot in Mdots:\n", " for T0 in T0s:\n", " try:\n", " sim = tools.Sim(\n", " Tstrucpath + \"/parker_\" + T0 + \"_\" + Mdot + \"/converged\"\n", " ) # load the converged simulation\n", " ClT, ClsigmaT = metaHe_weighted_T(\n", " sim\n", " ) # find the T_He and sigT_He for this simulation\n", " T_He.loc[Mdot][T0] = ClT\n", " sigmaT.loc[Mdot][T0] = ClsigmaT\n", " dT.loc[Mdot][T0] = ClT - int(T0)\n", "\n", " except FileNotFoundError: # then this profile was not calculated\n", " pass\n", "\n", " return dT, sigmaT\n", "\n", "\n", "### Functions that calculate the Bayesian terms ###\n", "\n", "\n", "def calc_likelihood_resolved(chisq_fit):\n", " \"\"\"\n", " Calculates the likelihood of a spectrally resolved fit to data,\n", " based on the chi-squared value.\n", "\n", " See e.g. https://philuttley.github.io/statistical-inference/10-mle_model_fitting/index.html\n", " for how the chi-squared statistic relates to the likelihood.\n", " \"\"\"\n", "\n", " likelihood = np.exp(-chisq_fit.values / 2.0)\n", "\n", " likelihood = pd.DataFrame(\n", " columns=chisq_fit.columns.values,\n", " index=chisq_fit.index.values,\n", " data=likelihood,\n", " dtype=float,\n", " ) # turn into pd.DataFrame\n", "\n", " return likelihood\n", "\n", "\n", "def calc_likelihood_unresolved(nsig_fit):\n", " \"\"\"\n", " Calculates the likelihood of spectrally unresolved data,\n", " expects a pandas dataframe that stores the number of\n", " errorbars difference between model and data.\n", " \"\"\"\n", "\n", " likelihood = sps.norm.pdf(nsig_fit.abs())\n", "\n", " likelihood = pd.DataFrame(\n", " columns=nsig_fit.columns.values,\n", " index=nsig_fit.index.values,\n", " data=likelihood,\n", " dtype=float,\n", " ) # turn into pd.DataFrame\n", "\n", " return likelihood\n", "\n", "\n", "def calc_prior(dT, sigmaT):\n", " \"\"\"\n", " Calculates a prior based on the model temperature self-consistency.\n", " See Linssen et al. (2022) for details on the choice of prior evaluation.\n", " \"\"\"\n", "\n", " # check that dT and sigmaT are defined on the same T-Mdot grid\n", " assert np.array_equal(\n", " dT.index.values, sigmaT.index.values\n", " ), \"Different temperature grids.\"\n", " assert np.array_equal(\n", " dT.columns.values, sigmaT.columns.values\n", " ), \"Different mass-loss rate grids.\"\n", "\n", " T0grid = dT.index.values.astype(float)\n", " Mdotgrid = dT.columns.values.astype(float)\n", "\n", " # assume prior is a normal distribution around dT=0 with std dev sigmaT\n", " number_sigma = dT / sigmaT\n", " prior = sps.norm.pdf(number_sigma)\n", " # now normalize the prior\n", " prior[np.isnan(prior)] = 0.0\n", " prior_sum = trapezoid(trapezoid(prior, axis=1, x=Mdotgrid), x=T0grid)\n", " prior = prior / prior_sum\n", "\n", " prior = pd.DataFrame(\n", " columns=dT.columns.values, index=dT.index.values, data=prior, dtype=float\n", " ) # turn into pd.DataFrame\n", "\n", " return prior\n", "\n", "\n", "def calc_posterior(prior, likelihood):\n", " \"\"\"\n", " Combines the prior from the T0-T_He analysis with the likelihoods\n", " from the data fit with Bayes' theorem to calculate the posterior.\n", " \"\"\"\n", "\n", " if isinstance(prior, pd.DataFrame):\n", " # check that the prior and likelihoods are defined on the same T-Mdot grid\n", " assert np.array_equal(\n", " prior.index.values, likelihood.index.values\n", " ), \"Different temperature grids.\"\n", " assert np.array_equal(\n", " prior.columns.values, likelihood.columns.values\n", " ), \"Different mass-loss rate grids.\"\n", "\n", " T0grid = likelihood.index.values.astype(float)\n", " Mdotgrid = likelihood.columns.values.astype(float)\n", "\n", " integrand = prior * likelihood\n", " integrand[np.isnan(integrand)] = 0.0\n", " evidence = trapezoid(trapezoid(integrand, axis=1, x=Mdotgrid), x=T0grid)\n", " posterior = integrand / evidence\n", "\n", " posterior = pd.DataFrame(\n", " columns=likelihood.columns.values,\n", " index=likelihood.index.values,\n", " data=posterior,\n", " dtype=float,\n", " ) # turn into pd.DataFrame\n", "\n", " # calculate marginalized posteriors\n", " post_Mdot = trapezoid(posterior.values, axis=0, x=T0grid)\n", " post_T0 = trapezoid(posterior.values, axis=1, x=Mdotgrid)\n", "\n", " post_Mdot_cum = cumulative_trapezoid(post_Mdot, x=Mdotgrid)\n", " lowsig_Mdot = Mdotgrid[np.argmax(post_Mdot_cum > 0.16)]\n", " mid_Mdot = Mdotgrid[np.argmin(np.abs(post_Mdot_cum - 0.5)) + 1]\n", " upsig_Mdot = Mdotgrid[np.argmax(post_Mdot_cum > 0.84) + 1]\n", " print(\n", " \"Constraints from marginalized posteriors (not necessarily normally distributed!):\"\n", " )\n", " print(f\"log10(Mdot) = {mid_Mdot} + {upsig_Mdot-mid_Mdot} - {mid_Mdot-lowsig_Mdot}\")\n", "\n", " post_T0_cum = cumulative_trapezoid(post_T0, x=T0grid)\n", " lowsig_T0 = T0grid[np.argmax(post_T0_cum > 0.16)]\n", " mid_T0 = T0grid[np.argmin(np.abs(post_T0_cum - 0.5)) + 1]\n", " upsig_T0 = T0grid[np.argmax(post_T0_cum > 0.84) + 1]\n", " print(f\"T0 = {mid_T0} + {upsig_T0-mid_T0} - {mid_T0-lowsig_T0}\")\n", "\n", " return posterior\n", "\n", "\n", "# Functions that plot stuff\n", "\n", "\n", "def plot_chisq_fit(\n", " chisq_fit,\n", " dof=None,\n", " bounds_T0=None,\n", " bounds_Mdot=None,\n", " title=None,\n", " save=None,\n", " fig=None,\n", " ax=None,\n", "):\n", " \"\"\"\n", " Makes a standard plot of the Parker wind models fitted to the data.\n", " In this case, we fit the resolved line shape, and\n", " the colormap shows the chi squared statistic.\n", " \"\"\"\n", "\n", " pixT = float(chisq_fit.index[1]) - float(chisq_fit.index[0]) # T step size\n", " pixM = float(chisq_fit.columns[1]) - float(chisq_fit.columns[0]) # Mdot step size\n", "\n", " if dof != None: # then we assume you want to plot the reduced chi squared\n", " plotted_vals = chisq_fit / dof\n", " cbar_label = r\"Reduced $\\chi ^2$\"\n", " else:\n", " plotted_vals = chisq_fit\n", " cbar_label = r\"$\\chi ^2$\"\n", "\n", " if fig == None and ax == None:\n", " _showplot = True\n", " fig, ax = plt.subplots(1)\n", " else:\n", " _showplot = False\n", "\n", " im = ax.imshow(\n", " plotted_vals.T,\n", " cmap=\"Blues_r\",\n", " norm=matplotlib.colors.LogNorm(),\n", " origin=\"lower\",\n", " aspect=\"auto\",\n", " interpolation=\"none\",\n", " extent=[\n", " float(chisq_fit.index[0]) - 0.5 * pixT,\n", " float(chisq_fit.index[-1]) + 0.5 * pixT,\n", " float(chisq_fit.columns[0]) - 0.5 * pixM,\n", " float(chisq_fit.columns[-1]) + 0.5 * pixM,\n", " ],\n", " )\n", " ax.set_facecolor(\"grey\")\n", " ax.set_xlabel(r\"$T_0$ [K]\")\n", " ax.set_ylabel(r\"log$_{10}$($\\dot{M}$ [g s$^{-1}$])\")\n", " fig.colorbar(im, label=cbar_label)\n", " if title != None:\n", " ax.set_title(title)\n", " if bounds_T0 != None:\n", " ax.set_xlim(*bounds_T0)\n", " if bounds_Mdot != None:\n", " ax.set_ylim(*bounds_Mdot)\n", " if save != None:\n", " plt.savefig(save, bbox_inches=\"tight\", dpi=300)\n", " if _showplot:\n", " plt.show()\n", "\n", "\n", "def plot_nsig_fit(\n", " nsig_fit, bounds_T0=None, bounds_Mdot=None, title=None, save=None, fig=None, ax=None\n", "):\n", " \"\"\"\n", " Makes a standard plot of the Parker wind models fitted to the data.\n", " In this case, we fit the EW of the model to the EW of the data, and\n", " the colormap simply indicates how many data errorbars difference there is.\n", " \"\"\"\n", "\n", " cmap = plt.get_cmap(\"Blues_r\")\n", " normalize = matplotlib.colors.Normalize(vmin=0, vmax=5)\n", "\n", " pixT = float(nsig_fit.index[1]) - float(nsig_fit.index[0]) # T step size\n", " pixM = float(nsig_fit.columns[1]) - float(nsig_fit.columns[0]) # Mdot step size\n", "\n", " if fig == None and ax == None:\n", " _showplot = True\n", " fig, ax = plt.subplots(1)\n", " else:\n", " _showplot = False\n", "\n", " im = ax.imshow(\n", " nsig_fit.T.abs(),\n", " cmap=cmap,\n", " norm=normalize,\n", " origin=\"lower\",\n", " aspect=\"auto\",\n", " interpolation=\"none\",\n", " extent=[\n", " float(nsig_fit.index[0]) - 0.5 * pixT,\n", " float(nsig_fit.index[-1]) + 0.5 * pixT,\n", " float(nsig_fit.columns[0]) - 0.5 * pixM,\n", " float(nsig_fit.columns[-1]) + 0.5 * pixM,\n", " ],\n", " )\n", " ax.set_facecolor(\"grey\")\n", " ax.set_xlabel(r\"$T_0$ [K]\")\n", " ax.set_ylabel(r\"log$_{10}$($\\dot{M}$ [g s$^{-1}$])\")\n", " fig.colorbar(im, label=r\"$\\Delta\\sigma$\")\n", " if title != None:\n", " ax.set_title(title)\n", " if bounds_T0 != None:\n", " ax.set_xlim(*bounds_T0)\n", " if bounds_Mdot != None:\n", " ax.set_ylim(*bounds_Mdot)\n", " if save != None:\n", " plt.savefig(save, bbox_inches=\"tight\", dpi=300)\n", " if _showplot:\n", " plt.show()\n", "\n", "\n", "def plot_selfcons(\n", " dT, sigmaT, bounds_T0=None, bounds_Mdot=None, title=None, fig=None, ax=None\n", "):\n", " \"\"\"\n", " Makes a standard plot of the self-consistency of the Parker wind parameter space.\n", " Self-consistent models are white, while blue/red models indicate models for which\n", " Cloudy indicates that the temperature in the helium line-forming region is\n", " cooler/hotter than that assumed in creating the isothermal profiles, respectively.\n", " The dotted lines indicate 1 standard deviation discrepancy.\n", " \"\"\"\n", "\n", " cmap = truncate_colormap(plt.get_cmap(\"seismic\"), 0.2, 0.8)\n", " normalize = matplotlib.colors.Normalize(vmin=-4000, vmax=4000)\n", "\n", " pixT = float(dT.index[1]) - float(dT.index[0]) # T step size\n", " pixM = float(dT.columns[1]) - float(dT.columns[0]) # Mdot step size\n", "\n", " if fig == None and ax == None:\n", " _showplot = True\n", " fig, ax = plt.subplots(1)\n", " else:\n", " _showplot = False\n", "\n", " im = ax.imshow(\n", " dT.T,\n", " cmap=cmap,\n", " norm=normalize,\n", " origin=\"lower\",\n", " aspect=\"auto\",\n", " interpolation=\"none\",\n", " extent=[\n", " float(dT.index[0]) - 0.5 * pixT,\n", " float(dT.index[-1]) + 0.5 * pixT,\n", " float(dT.columns[0]) - 0.5 * pixM,\n", " float(dT.columns[-1]) + 0.5 * pixM,\n", " ],\n", " )\n", " # plot the 1 sigma line (i.e. dT = sigmaT)\n", " sig1lines = ax.contour(\n", " dT.index.astype(float),\n", " dT.columns.astype(float),\n", " (dT / sigmaT).T,\n", " levels=[-1, 1],\n", " zorder=1,\n", " colors=\"black\",\n", " linestyles=\"dotted\",\n", " )\n", " ax.set_facecolor(\"grey\")\n", " ax.set_xlabel(r\"$T_0$ [K]\")\n", " ax.set_ylabel(r\"log$_{10}$($\\dot{M}$ [g s$^{-1}$])\")\n", " fig.colorbar(im, label=r\"$T_{He} - T_0$ [K]\")\n", " if title != None:\n", " ax.set_title(title)\n", " if bounds_T0 != None:\n", " ax.set_xlim(*bounds_T0)\n", " if bounds_Mdot != None:\n", " ax.set_ylim(*bounds_Mdot)\n", " if _showplot:\n", " plt.show()\n", "\n", "\n", "def plot_posterior(\n", " posterior, bounds_T0=None, bounds_Mdot=None, title=None, save=None, log=False\n", "):\n", " \"\"\"\n", " Plots the posterior distribution and the marginalized distributions.\n", " Prints the 1 sigma confidence intervals, at the precision of the T0/Mdot-grid.\n", " \"\"\"\n", "\n", " T0grid = posterior.index.values.astype(float)\n", " Mdotgrid = posterior.columns.values.astype(float)\n", " post = posterior.values\n", "\n", " pixT = np.diff(T0grid)[0] # T step size\n", " pixM = np.diff(Mdotgrid)[0] # Mdot step size\n", "\n", " # calculate marginalized posteriors\n", " post_Mdot = trapezoid(post, axis=0, x=T0grid)\n", " post_T0 = trapezoid(post, axis=1, x=Mdotgrid)\n", "\n", " post_Mdot_cum = cumulative_trapezoid(post_Mdot, x=Mdotgrid)\n", " lowsig_Mdot = Mdotgrid[np.argmax(post_Mdot_cum > 0.16)]\n", " mid_Mdot = Mdotgrid[np.argmin(np.abs(post_Mdot_cum - 0.5)) + 1]\n", " upsig_Mdot = Mdotgrid[np.argmax(post_Mdot_cum > 0.84) + 1]\n", "\n", " post_T0_cum = cumulative_trapezoid(post_T0, x=T0grid)\n", " lowsig_T0 = T0grid[np.argmax(post_T0_cum > 0.16)]\n", " mid_T0 = T0grid[np.argmin(np.abs(post_T0_cum - 0.5)) + 1]\n", " upsig_T0 = T0grid[np.argmax(post_T0_cum > 0.84) + 1]\n", "\n", " fig = plt.figure(figsize=(5, 5))\n", " gs = fig.add_gridspec(\n", " 2,\n", " 2,\n", " width_ratios=(3, 1),\n", " height_ratios=(1, 3),\n", " left=0.1,\n", " right=0.9,\n", " bottom=0.1,\n", " top=0.9,\n", " wspace=0.05,\n", " hspace=0.05,\n", " )\n", " ax = fig.add_subplot(gs[1, 0])\n", " ax_T0 = fig.add_subplot(gs[0, 0])\n", " ax_Mdot = fig.add_subplot(gs[1, 1])\n", "\n", " im = ax.imshow(\n", " post.T,\n", " cmap=plt.get_cmap(\"Greys\"),\n", " origin=\"lower\",\n", " aspect=\"auto\",\n", " interpolation=\"none\",\n", " extent=[\n", " T0grid[0] - 0.5 * pixT,\n", " T0grid[-1] + 0.5 * pixT,\n", " Mdotgrid[0] - 0.5 * pixM,\n", " Mdotgrid[-1] + 0.5 * pixM,\n", " ],\n", " )\n", " if log:\n", " im.set_norm(matplotlib.colors.LogNorm())\n", " ax_T0.plot(T0grid, post_T0, color=\"k\")\n", " ax_Mdot.plot(post_Mdot, Mdotgrid, color=\"k\")\n", "\n", " ax_T0.axvline(mid_T0, color=\"blue\", linewidth=0.7)\n", " ax_T0.axvline(lowsig_T0, color=\"blue\", linewidth=0.7, linestyle=\"dotted\")\n", " ax_T0.axvline(upsig_T0, color=\"blue\", linewidth=0.7, linestyle=\"dotted\")\n", " ax_Mdot.axhline(mid_Mdot, color=\"blue\", linewidth=0.7)\n", " ax_Mdot.axhline(lowsig_Mdot, color=\"blue\", linewidth=0.7, linestyle=\"dotted\")\n", " ax_Mdot.axhline(upsig_Mdot, color=\"blue\", linewidth=0.7, linestyle=\"dotted\")\n", "\n", " ax_T0.set_xticks([])\n", " ax_T0.set_yticks([])\n", " ax_Mdot.set_xticks([])\n", " ax_Mdot.set_yticks([])\n", "\n", " ax.set_xlim(T0grid[0] - 0.5 * pixT, T0grid[-1] + 0.5 * pixT)\n", " ax_T0.set_xlim(T0grid[0] - 0.5 * pixT, T0grid[-1] + 0.5 * pixT)\n", " ax.set_ylim(Mdotgrid[0] - 0.5 * pixM, Mdotgrid[-1] + 0.5 * pixM)\n", " ax_Mdot.set_ylim(Mdotgrid[0] - 0.5 * pixM, Mdotgrid[-1] + 0.5 * pixM)\n", " if bounds_T0 != None:\n", " ax_T0.set_xlim(*bounds_T0)\n", " ax.set_xlim(*bounds_T0)\n", " if bounds_Mdot != None:\n", " ax_Mdot.set_ylim(*bounds_Mdot)\n", " ax.set_ylim(*bounds_Mdot)\n", "\n", " ax.set_xlabel(r\"$T_0$ [K]\")\n", " ax.set_ylabel(r\"log$_{10}$($\\dot{M}$ [g s$^{-1}$])\")\n", " if title != None:\n", " ax_T0.set_title(title)\n", "\n", " if save != None:\n", " plt.savefig(save, bbox_inches=\"tight\", dpi=300)\n", " plt.show()\n", "\n", "\n", "def plot_joint_constraint_resolved(\n", " dT,\n", " sigmaT,\n", " chisq_fit,\n", " posterior,\n", " bounds_T0=None,\n", " bounds_Mdot=None,\n", " title=None,\n", " save=None,\n", " fig=None,\n", " ax=None,\n", " post_uptosigma=1,\n", " cmap_T=None,\n", " cmap_fit=None,\n", " show_colorbar=True,\n", "):\n", " \"\"\"\n", " Makes a standard plot of the posterior distribution,\n", " on top of the prior (self-consistency) and likelihood (data fit) constraints.\n", " Also prints out the 1sigma bounds on T and Mdot.\n", "\n", " chisq_fit must be given, but dT, sigmaT and posterior can be passed as None\n", " if they should not be plotted.\n", " \"\"\"\n", "\n", " assert post_uptosigma <= 3, \"Maximum value of post_uptosigma is 3\"\n", "\n", " if cmap_T is None:\n", " cmap_T = truncate_colormap(plt.get_cmap(\"autumn\"), 0.35, 0.65)\n", " if cmap_fit is None:\n", " cmap_fit = truncate_colormap(plt.get_cmap(\"winter_r\"), 0.35, 0.65)\n", "\n", " if fig == None and ax == None:\n", " showplot = True\n", " fig, ax = plt.subplots(1)\n", " else:\n", " showplot = False\n", "\n", " # Plot the model temperature discrepancy (i.e. prior)\n", " if dT is not None and sigmaT is not None:\n", " nsig_T = dT / sigmaT\n", "\n", " x, y, z = get_xyz(nsig_T)\n", "\n", " ax.tricontour(\n", " x, y, z, cmap=cmap_T, levels=[0], zorder=0\n", " ) # plots the orange line\n", " prior_contour = ax.tricontourf(\n", " x, y, np.abs(z), cmap=cmap_T, levels=[0, 1, 2], alpha=0.5, zorder=0\n", " ) # plots the orange contours\n", "\n", " # Plot the fit constraints (i.e. likelihood):\n", " likelihood = calc_likelihood_resolved(chisq_fit)\n", " likelihood = likelihood.fillna(0.0)\n", " highest_likelihood = likelihood.max().max()\n", " sigma1_likelihood = highest_likelihood * np.exp(-1 / 2)\n", " sigma2_likelihood = highest_likelihood * np.exp(-4 / 2)\n", " sigma3_likelihood = highest_likelihood * np.exp(-9 / 2)\n", "\n", " x, y, z = get_xyz(likelihood)\n", "\n", " likelihood_contour = ax.tricontourf(\n", " x,\n", " y,\n", " z,\n", " cmap=cmap_fit,\n", " levels=[sigma2_likelihood, sigma1_likelihood, highest_likelihood],\n", " alpha=0.5,\n", " zorder=1,\n", " ) # plots the blue contours\n", "\n", " # Plot the posterior\n", " if posterior is not None and post_uptosigma > 0:\n", " # calculate the posterior values that corresponds to the 1,2,3 sigma contours\n", " posterior_sorted = np.sort(posterior.values.flatten())[::-1]\n", " posterior_sum = np.sum(posterior_sorted)\n", " posterior_cumsum = np.cumsum(posterior_sorted)\n", " p1sigma = posterior_sorted[posterior_cumsum > 0.3935 * posterior_sum][0]\n", " p2sigma = posterior_sorted[posterior_cumsum > 0.8647 * posterior_sum][0]\n", " p3sigma = posterior_sorted[posterior_cumsum > 0.9889 * posterior_sum][0]\n", " post_levels = [p3sigma, p2sigma, p1sigma][3 - post_uptosigma :]\n", " post_linestyles = [\"dotted\", \"dashed\", \"solid\"][3 - post_uptosigma :]\n", "\n", " x, y, z = get_xyz(posterior)\n", "\n", " # plot the joint posterior 1/2/3 sigma lines:\n", " cn_post = ax.tricontour(\n", " x,\n", " y,\n", " z,\n", " levels=post_levels,\n", " linestyles=post_linestyles,\n", " colors=\"black\",\n", " zorder=2,\n", " )\n", "\n", " # print the posterior 1sigma credible intervals:\n", " pT0, pMdot = get_array_from_contour(\n", " cn_post, -1\n", " ) # get last level (=1 sigma) from posterior\n", " if pT0 is not None and pMdot is not None:\n", " bestT0_index, bestMdot_index = np.unravel_index(\n", " np.nanargmax(posterior), posterior.shape\n", " )\n", " bestT0, bestMdot = float(posterior.index[bestT0_index]), float(\n", " posterior.columns[bestMdot_index]\n", " ) # T, Mdot of the maximum of the posterior\n", " max_pMdot, min_pMdot, max_pT0, min_pT0 = (\n", " np.max(pMdot),\n", " np.min(pMdot),\n", " np.max(pT0),\n", " np.min(pT0),\n", " ) #\n", " print(\n", " f\"1 sigma constraints:\\nlog10(Mdot) = {bestMdot} + {max_pMdot-bestMdot} - {bestMdot-min_pMdot}\"\n", " )\n", " print(f\"T0 = {bestT0} + {max_pT0-bestT0} - {bestT0-min_pT0}\")\n", "\n", " if show_colorbar:\n", " divider = make_axes_locatable(ax)\n", "\n", " if dT is not None and sigmaT is not None:\n", " # plot the prior colorbar\n", " cax1 = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", " cbar1 = fig.colorbar(prior_contour, cax=cax1)\n", " cax1.yaxis.set_major_locator(\n", " plt.NullLocator()\n", " ) # turn the ticks for cbar1 off\n", "\n", " # plot the likelihood colorbar\n", " cax2 = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", " cbar2 = fig.colorbar(likelihood_contour, cax=cax2)\n", " cax2.set_yticklabels([r\"2$\\sigma$\", r\"1$\\sigma$\", r\"0$\\sigma$\"])\n", " cax2.invert_yaxis() # because the p-values and sigma values are in opposite direction\n", "\n", " ax.set_xlabel(r\"$T_0$ [K]\")\n", " ax.set_ylabel(r\"log$_{10}$($\\dot{M}$ [g s$^{-1}$])\")\n", " if bounds_T0:\n", " ax.set_xlim(*bounds_T0)\n", " if bounds_Mdot:\n", " ax.set_ylim(*bounds_Mdot)\n", " ax.set_yticks(np.arange(bounds_Mdot[0], bounds_Mdot[1] + 0.5, 0.5))\n", " if title != None:\n", " ax.set_title(title)\n", " if save != None:\n", " plt.savefig(save, bbox_inches=\"tight\", dpi=300)\n", " if showplot:\n", " plt.show()\n", "\n", "\n", "def plot_joint_constraint_unresolved(\n", " dT,\n", " sigmaT,\n", " nsig_fit,\n", " posterior,\n", " bounds_T0=None,\n", " bounds_Mdot=None,\n", " title=None,\n", " save=None,\n", " fig=None,\n", " ax=None,\n", " cmap_T=None,\n", " cmap_fit=None,\n", " post_uptosigma=1,\n", " show_colorbar=True,\n", "):\n", " \"\"\"\n", " Makes a plot of the 1sigma posterior contour\n", " on top of the prior (self-consistency) and likelihood (data fit) constraints.\n", " Also prints out the 1sigma bounds on T and Mdot, which are calculated\n", " based on the matplotlib contour, and thus involve interpolation to\n", " higher precision than the data itself.\n", " \"\"\"\n", "\n", " assert post_uptosigma <= 3, \"Maximum value of post_uptosigma is 3\"\n", "\n", " if cmap_T is None:\n", " cmap_T = truncate_colormap(plt.get_cmap(\"autumn\"), 0.35, 0.65)\n", " if cmap_fit is None:\n", " cmap_fit = truncate_colormap(plt.get_cmap(\"winter\"), 0.35, 0.65)\n", "\n", " if fig == None and ax == None:\n", " _showplot = True\n", " fig, ax = plt.subplots(1)\n", " else:\n", " _showplot = False\n", "\n", " # Plot the model temperature discrepancy (i.e. prior)\n", " if dT is not None and sigmaT is not None:\n", " nsig_T = dT / sigmaT\n", "\n", " x, y, z = get_xyz(nsig_T)\n", "\n", " ax.tricontour(\n", " x, y, z, cmap=cmap_T, levels=[0], zorder=0\n", " ) # plots the orange line\n", " prior_contour = ax.tricontourf(\n", " x, y, np.abs(z), cmap=cmap_T, levels=[0, 1, 2], alpha=0.5, zorder=0\n", " ) # plots the orange contours\n", "\n", " # Plot the fit constraints (i.e. likelihood):\n", " x, y, z = get_xyz(nsig_fit)\n", " ax.tricontour(x, y, z, cmap=cmap_fit, levels=[0], zorder=1) # plots the blue line\n", " likelihood_contour = ax.tricontourf(\n", " x, y, np.abs(z), cmap=cmap_fit, levels=[0, 1, 2], alpha=0.5, zorder=1\n", " ) # plots the blue contours\n", "\n", " # Plot the posterior\n", " if posterior is not None and post_uptosigma > 0:\n", " # calculate the posterior values that corresponds to the 1,2,3 sigma contours\n", " posterior_sorted = np.sort(posterior.values.flatten())[::-1]\n", " posterior_sum = np.sum(posterior_sorted)\n", " posterior_cumsum = np.cumsum(posterior_sorted)\n", " p1sigma = posterior_sorted[posterior_cumsum > 0.3935 * posterior_sum][0]\n", " p2sigma = posterior_sorted[posterior_cumsum > 0.8647 * posterior_sum][0]\n", " p3sigma = posterior_sorted[posterior_cumsum > 0.9889 * posterior_sum][0]\n", " post_levels = [p3sigma, p2sigma, p1sigma][3 - post_uptosigma :]\n", " post_linestyles = [\"dotted\", \"dashed\", \"solid\"][3 - post_uptosigma :]\n", "\n", " x, y, z = get_xyz(posterior)\n", "\n", " # plot the joint posterior 1/2/3 sigma lines:\n", " cn_post = ax.tricontour(\n", " x,\n", " y,\n", " z,\n", " levels=post_levels,\n", " linestyles=post_linestyles,\n", " colors=\"black\",\n", " zorder=2,\n", " )\n", "\n", " # print the posterior 1sigma credible intervals:\n", " pT0, pMdot = get_array_from_contour(\n", " cn_post, -1\n", " ) # get last level (=1 sigma) from posterior\n", " if pT0 is not None and pMdot is not None:\n", " bestT0_index, bestMdot_index = np.unravel_index(\n", " np.nanargmax(posterior), posterior.shape\n", " )\n", " bestT0, bestMdot = float(posterior.index[bestT0_index]), float(\n", " posterior.columns[bestMdot_index]\n", " ) # T, Mdot of the maximum of the posterior\n", " max_pMdot, min_pMdot, max_pT0, min_pT0 = (\n", " np.max(pMdot),\n", " np.min(pMdot),\n", " np.max(pT0),\n", " np.min(pT0),\n", " ) #\n", " print(\n", " f\"1 sigma constraints:\\nlog10(Mdot) = {bestMdot} + {max_pMdot-bestMdot} - {bestMdot-min_pMdot}\"\n", " )\n", " print(f\"T0 = {bestT0} + {max_pT0-bestT0} - {bestT0-min_pT0}\")\n", "\n", " if show_colorbar:\n", " divider = make_axes_locatable(ax)\n", "\n", " if dT is not None and sigmaT is not None:\n", " # plot the prior colorbar\n", " cax1 = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", " cbar1 = fig.colorbar(prior_contour, cax=cax1)\n", " cax1.yaxis.set_major_locator(\n", " plt.NullLocator()\n", " ) # turn the ticks for cbar1 off\n", "\n", " # plot the likelihood colorbar\n", " cax2 = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", " cbar2 = fig.colorbar(likelihood_contour, cax=cax2)\n", " cax2.set_yticklabels([r\"0$\\sigma$\", r\"1$\\sigma$\", r\"2$\\sigma$\"])\n", "\n", " ax.set_xlabel(r\"$T_0$ [K]\")\n", " ax.set_ylabel(r\"log$_{10}$($\\dot{M}$ [g s$^{-1}$])\")\n", " if bounds_T0:\n", " ax.set_xlim(*bounds_T0)\n", " if bounds_Mdot:\n", " ax.set_ylim(*bounds_Mdot)\n", " ax.set_yticks(np.arange(bounds_Mdot[0], bounds_Mdot[1] + 0.5, 0.5))\n", " if title != None:\n", " ax.set_title(title)\n", " if save != None:\n", " plt.savefig(save, bbox_inches=\"tight\", dpi=300)\n", " if _showplot:\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "aa8609df", "metadata": {}, "source": [ "## Retrieval of a mock spectrum of a hot Neptune" ] }, { "cell_type": "markdown", "id": "dea2e2a7-78c6-4974-b3e9-4028834d87fd", "metadata": {}, "source": [ "### Step 1: Create a grid of Parker wind profiles with p-winds\n", "\n", "This step can take a long time and is therefore skipped by default. We provide pre-run models that are used in this notebook. Uncomment the cell below if you want to run your own models.\n", "\n", "This step can be done with the `construct_parker.run()` function, with the proper arguments. Running `help(construct_parker.run)` will give an overview of the available arguments. In our case, we will make the Parker profiles with a solar composition (so including metals at the solar abundances). This is done by passing `z=1`. We will run a grid of models that span temperatures from 4000 to 7000 K in steps of 100 K, and log10-mass-loss rates from 10.5 to 11.5 in steps of 0.05 dex. To run a grids of these parameters, we simply supply a list or array. Since in the future, we might want to explore different compositions, `construct_parker.run()` always expects you to give a folder name `pdir` where we want to store our Parker profiles. We recommend using a descriptive name, so in this case we will go with `\"z_1\"` and the path where our profiles will be saved is then \\$SUNBATHER_PROJECT\\_PATH/parker\\_profiles/hotNeptune/z\\_1/. Then, tidal gravity is included by default, but we will turn it off here in order to reproduce the results from Linssen et al. (2024), by passing `no_tidal=True`. Finally, with the `cores` argument you can specify the number of parallel CPUs used, we will use 8 but feel free to change it. \n", "\n", "If you run this step, for this set of parameters, it can take a long time (hours). In the \\$SUNBATHER\\_PROJECT_PATH/parker\\_profiles/hotNeptune/z\\_1/ folder, there should now be many different .txt files with the isothermal Parker wind structures, feel free to inspect them!" ] }, { "cell_type": "code", "execution_count": null, "id": "4c61a5f8-b946-4220-859d-3bee7000a0cf", "metadata": {}, "outputs": [], "source": [ "temp_min = 4000\n", "temp_max = 7000\n", "temp_step = 100\n", "temp = np.linspace(\n", " temp_min,\n", " temp_max,\n", " 1 + (temp_max - temp_min) // temp_step\n", ")\n", "mdot_min = 10.5\n", "mdot_max = 11.5\n", "mdot_step = 0.05\n", "mdot = np.linspace(\n", " mdot_min,\n", " mdot_max,\n", " 1 + int((mdot_max - mdot_min) / mdot_step)\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "2c63dd23", "metadata": {}, "outputs": [], "source": [ "# The full command to create a grid of solar composition Parker wind models is:\n", "# (only uncomment if you want to run these models which takes a long time)\n", "\n", "# import numpy as np\n", "# from sunbather import construct_parker\n", "#\n", "# construct_parker.run(\n", "# plname=\"hotNeptune\",\n", "# pdir=\"z_1\",\n", "# temp=temp,\n", "# mdot=mdot,\n", "# z=1,\n", "# no_tidal=True,\n", "# cores=8,\n", "# )" ] }, { "cell_type": "markdown", "id": "4d05d58c-6821-498e-9f63-df7a6cdd0523", "metadata": {}, "source": [ "### Step 2: Run the Parker wind profiles through Cloudy\n", "\n", "This step can take a long time and is therefore skipped by default. We provide pre-run models that are used in this notebook. Uncomment the cell below if you want to run your own models.\n", "\n", "This step can be done with the `convergeT_parker.run()` function, with the proper arguments (`help(convergeT_parker.run)` will explain these). The temperature and mass-loss rate grid-commands are the same as in Step 1. We also again need to specify a folder name where we want to save our *Cloudy* simulations (which is a different directory from where the isothermal Parker wind profiles were stored). In this case we will indicate that we used _Cloudy_ with a solar composition/metallicity, i.e. a metallicity of z=1: `workingdir=\"z_1\"` and so our simulations will be saved in \\$SUNBATHER\\_PROJECT\\_PATH/sims/1D/hotNeptune/z\\_1/. We now again need to specify the folder where we want to read the Parker wind profiles from: `pdir=\"z_1\"`. Since the *Cloudy* simulations are more time-consuming, you can easily run the different simulations in parallel by using the `cores` argument. We then need to choose the atmospheric composition in _Cloudy_. We go with the default solar composition so we do not have to provide an argument for it. The last thing we need to think about, is for which atomic/ionic species we want to save *Cloudy's* output. The default behavior is to save everything that's available, but this generally results in large file sizes (~5 MB per model, if metals are included). In our particular case, we want to fit helium observations, so we are fine with just saving the densities of the different atomic helium energy levels: `save_sp=[\"He\",]`. \n", "\n", "If you run this step, for this set of parameters, it should take many hours (roughly 10 minutes per model, and we are running roughy 30x20 models), and it is adviced to run it on a compute cluster. In the \\$SUNBATHER\\_PROJECT\\_PATH/sims/1D/hotNeptune/fH\\_0.9/ folder, there should now be many different sub-folders, with the output of the _Cloudy_ simulations, feel free to inspect the files! The _converged.png_ file shows the converged temperature structure, and the other _converged.*_ files are the _Cloudy_ output files." ] }, { "cell_type": "code", "execution_count": null, "id": "bbdc7986-6d55-4e35-aa04-5948e0d717f0", "metadata": {}, "outputs": [], "source": [ "# The full command to run our grid of hydrogen/helium Parker wind models through Cloudy is:\n", "# (only uncomment if you want to run these models which takes a long time)\n", "\n", "\n", "# import sunbather.convergeT_parker\n", "#\n", "# sunbather.convergeT_parker.run(\n", "# plname=\"hotNeptune\",\n", "# workingdir=\"z_1\",\n", "# pdir=\"z_1\",\n", "# temp=temp,\n", "# mdot=mdot,\n", "# save_sp=[\"He\",],\n", "# )" ] }, { "cell_type": "markdown", "id": "cc14c0e4", "metadata": {}, "source": [ "### Step 3: Make transit spectra of the metastable helium triplet and fit the data (to obtain _likelihood_)\n", "\n", "Since we first have to make our mock spectrum, and also because we perform two analyses (one fitting the spectrally resolved line shape, and one fitting the equivalent width), step 3 is broken up into three substeps here.\n", "\n", "To make transit spectra, we can make use of the `FinFout()` module in `RT.py`. In our case, we will first calculate and plot the helium spectrum of the chosen model, which will end up being our mock spectrum. Then, we run all models and calculate their spectra, which we can compare/fit to the mock spectrum and calculate the likelihood of each model.\n", "\n", "The `FinFout()` function takes a `Sim` object, which is a class defined in the `tools` module. For more information on how to use this class, we refer to the wiki. The `FinFout()` function takes three required arguments: the `Sim` object, a wavelength array in angstrom and vacuum, and the atomic/ionic species (or a list thereof) which to include in the calculations. In our case, we only need atomic helium, i.e. \"He\". There are a few optional arguments (see wiki), for example to specify a limb-darkening law or transit phase, but we will not use those here." ] }, { "cell_type": "markdown", "id": "146b8339", "metadata": {}, "source": [ "### Step 3a: Creating a mock spectrum\n", "\n", "First, we create our mock spectrum by selecting one of the models that we ran, calculating its spectrum, and adding noise to it. If you ran your own grid of models in steps 1 and 2, you will already have the model and **should not execute the cell below**. If instead, you want to use our pre-run models, the cell below copies the necessary files, from the /materials folder to your projectpath." ] }, { "cell_type": "code", "execution_count": null, "id": "1392a8c0", "metadata": {}, "outputs": [], "source": [ "materials_folder = os.path.join(os.getcwd(), \"materials\")\n", "dirs_to_copy = {\n", " \"parker_profiles/hotNeptune/\": f\"{projectpath}/parker_profiles/hotNeptune\",\n", " \"sims/1D/hotNeptune\": f\"{projectpath}/sims/1D/hotNeptune\"\"\n", "}\n", "\n", "def copy_if_not_exists(src, dst):\n", " \"\"\"Recursively copy files from src to dst without overwriting existing files.\"\"\"\n", " for root, dirs, files in os.walk(src):\n", " rel_path = os.path.relpath(root, src) # Get relative path from source root\n", " dest_path = os.path.join(dst, rel_path) # Construct destination path\n", " \n", " os.makedirs(dest_path, exist_ok=True) # Ensure the folder exists\n", " \n", " for file in files:\n", " src_file = os.path.join(root, file)\n", " dest_file = os.path.join(dest_path, file)\n", " \n", " if not os.path.exists(dest_file): # Only copy if the file doesn't exist\n", " shutil.copy2(src_file, dest_file)\n", "\n", "for folder, destination in dirs_to_copy.items():\n", " source = os.path.join(materials_folder, folder)\n", "\n", " if os.path.exists(source):\n", " copy_if_not_exists(source, destination)\n", " else:\n", " print(f\"Warning: Source folder {source} does not exist.\")" ] }, { "cell_type": "code", "execution_count": null, "id": "5c08a50d", "metadata": {}, "outputs": [], "source": [ "# set random seed for consistent results\n", "np.random.seed(0)\n", "from sunbather import tools\n", "\"\"\"\n", "if due to a different numpy version, your generated mock spectrum is different,\n", "our generated mock spectrum can be found in 'materials/hotNeptune_helium_mock_spectrum.txt',\n", "which has stored np.column_stack((wavs, observed_excess_absorption, observed_excess_absorption_errorbar))\n", "\"\"\"\n", "\n", "# Let's make the spectrum of the generic hot Neptune\n", "hotNep = tools.Sim(\n", " projectpath + \"/sims/1D/hotNeptune/z_1/parker_5100_10.950/converged\"\n", ")\n", "\n", "wavs = RT.constantR_wavs(10830, 10836, 80000)\n", "true_spec, _, _ = RT.FinFout(hotNep, wavs, \"He\")\n", "\n", "# draw errors from normal distribution\n", "observed_errorbar = np.ones_like(wavs) * 0.0025 # just pick a number\n", "observed_spec = np.random.normal(loc=true_spec, scale=observed_errorbar)\n", "\n", "# convert into excess absorption in % which is what the calc_chisqs_He10830() function expects\n", "# we pretend we don't know the continuum level and take it as the average away from the line (last 15 points)\n", "observed_excess_absorption = 100 * (np.mean(observed_spec[-15:]) - observed_spec)\n", "observed_excess_absorption_errorbar = 100 * observed_errorbar\n", "\n", "highres_wavs = np.linspace(10830, 10836, 1000) # for smooth plotting\n", "highres_true_spec, _, _ = RT.FinFout(hotNep, highres_wavs, \"He\")\n", "\n", "\n", "fig, ax = plt.subplots(1)\n", "ax.plot(highres_wavs, highres_true_spec)\n", "ax.errorbar(wavs, observed_spec, yerr=observed_errorbar, fmt=\"o\")\n", "ax.set_xlabel(\"Wavelength [Å]\")\n", "ax.set_ylabel(r\"$F_{in}$/$F_{out}$\")\n", "plt.show()\n", "\n", "\n", "# integrate the difference with the continuum to obtain the equivalent width\n", "true_EW = trapezoid(np.max(true_spec) - true_spec, x=wavs) # units: Å\n", "# let's choose the errorbar as 1/8 of the EW so that we will get around a 8 sigma detection\n", "observed_EW_errorbar = true_EW / 8.0\n", "observed_EW = np.random.normal(loc=true_EW, scale=observed_EW_errorbar)\n", "print(\n", " f\"The true EW is {true_EW}. We have drawn an observed EW = {observed_EW} +- {observed_EW_errorbar} \"\n", " f\"(which is a {observed_EW/observed_EW_errorbar} sigma detection)\"\n", ")" ] }, { "cell_type": "markdown", "id": "42e2baf3", "metadata": {}, "source": [ "### Step 3b: Fitting the spectrally resolved line\n", "\n", "With our mock spectrum and its EW calculated, we can now proceed to fitting it with a grid of Parker wind models." ] }, { "cell_type": "markdown", "id": "0084a6f0", "metadata": {}, "source": [ "Uncomment the next cell if you have run your own sunbather simulations and want to calculate the chi-squared values here." ] }, { "cell_type": "code", "execution_count": null, "id": "afb89ddc", "metadata": {}, "outputs": [], "source": [ "# chisqs = calc_chisqs_He10830(tools.get_sunbather_project_path()+'/sims/1D/hotNeptune/z_1/',\n", "# wavs, observed_excess_absorption, observed_excess_absorption_errorbar)\n", "# chisqs.to_csv('materials/hotNeptune_chisqs.csv') #this overwrites the supplied file" ] }, { "cell_type": "code", "execution_count": null, "id": "98436f6e", "metadata": {}, "outputs": [], "source": [ "# read in the chi-squared values from file\n", "chisqs = pd.read_csv(\"materials/hotNeptune_chisqs.csv\", index_col=0, dtype=float)\n", "\n", "# calculate the likelihood (an actual p-value, contrary to the chi-squared statistic)\n", "likelihood_resolved = calc_likelihood_resolved(chisqs)" ] }, { "cell_type": "markdown", "id": "3d86307a", "metadata": {}, "source": [ "We can already plot the chi-squared values of the resolved line fit to see what they look like:" ] }, { "cell_type": "code", "execution_count": null, "id": "133461ec", "metadata": {}, "outputs": [], "source": [ "plot_chisq_fit(\n", " chisqs, bounds_T0=(temp_min, temp_max), bounds_Mdot=(mdot_min, mdot_max), title=\"Resolved line fit\"\n", ")" ] }, { "cell_type": "markdown", "id": "62b139bf", "metadata": {}, "source": [ "### Step 3c: Fitting the line equivalent width (EW)" ] }, { "cell_type": "markdown", "id": "4dd90a98", "metadata": {}, "source": [ "Uncomment the next cell if you have run your own sunbather simulations and want to calculate the model temperature self-consistency here." ] }, { "cell_type": "code", "execution_count": null, "id": "c86ecca2", "metadata": {}, "outputs": [], "source": [ "# EWs = calc_EWs_He10830(tools.get_sunbather_project_path()+'/sims/1D/hotNeptune/z_1/')\n", "# EWs.to_csv('materials/hotNeptune_EWs.csv') #this overwrites the supplied file" ] }, { "cell_type": "code", "execution_count": null, "id": "bc6f2960", "metadata": {}, "outputs": [], "source": [ "# read in the EW values from file\n", "EWs = pd.read_csv(\"materials/hotNeptune_EWs.csv\", index_col=0, dtype=float)\n", "\n", "# calculate the number of errorbars discrepancy between model and data\n", "# this is what the calc_likelihood_unresolved() function expects\n", "nsig_fit = (EWs - observed_EW) / observed_EW_errorbar\n", "\n", "likelihood_unresolved = calc_likelihood_unresolved(nsig_fit)" ] }, { "cell_type": "markdown", "id": "aa848b68", "metadata": {}, "source": [ "We can already plot the number of errorbars discrepancy of the EW fit to see what they look like:" ] }, { "cell_type": "code", "execution_count": null, "id": "0b5ded03", "metadata": {}, "outputs": [], "source": [ "plot_nsig_fit(\n", " nsig_fit,\n", " bounds_T0=(temp_min, temp_max),\n", " bounds_Mdot=(mdot_min, mdot_max),\n", " title=\"Equivalent width fit\",\n", ")" ] }, { "cell_type": "markdown", "id": "9a2de0c2", "metadata": {}, "source": [ "### Step 4: Constraining the parameter space by assessing temperature self-consistency (to obtain _prior_)\n", "\n", "We now assess the self-consistency of the Parker wind profiles in the helium line-forming region. To do that, we compare Cloudy's converged nonisothermal temperature profile to the initial isothermal value to identify models that are not self-consistent. As discussed in Linssen et al. (2022) and Linssen&Oklopcic (2023), there is not one single best way to do this comparison and it depends on the aim of the modeling effort. Since we are fitting helium observations here, we care most about the self-consistency of the Parker wind models in the region where the helium line forms. Thus, we calculate the mean of the nonisothermal temperature profiles weighted by the number density of metastable helium (this will be done by the `metaHe_weighted_T()` function and stored in the `T_He` variable and its standard deviation in the `sigmaT` variable). We then calculate the difference between this weighted temperature and the isothermal value $T_0$, which is the `dT` variable. This `dT` will be used for our prior on our Bayesian analysis. We are not confident in regions of the parameter space where |dT| >> 0 so we use a low prior there, while we are more confident in Parker wind models for which dT ~ 0 so we use a high prior there. This prior will only be used for the EW analysis. For the resolved fit analysis we use a flat prior." ] }, { "cell_type": "markdown", "id": "ce96d12c", "metadata": {}, "source": [ "Uncomment the next cell if you have run your own sunbather simulations and want to calculate the model equivalent width values here." ] }, { "cell_type": "code", "execution_count": null, "id": "a6d078e4", "metadata": {}, "outputs": [], "source": [ "# dT, sigmaT = calc_dT_helium(tools.get_sunbather_project_path()+'/sims/1D/hotNeptune/z_1/')\n", "# dT.to_csv('materials/hotNeptune_dT.csv', float_format='%.3e') #these overwrite the supplied files\n", "# sigmaT.to_csv('materials/hotNeptune_sigmaT.csv', float_format='%.3e')" ] }, { "cell_type": "code", "execution_count": null, "id": "14a91e2d", "metadata": {}, "outputs": [], "source": [ "# read in the dT and sigmaT values from file\n", "dT = pd.read_csv(\"materials/hotNeptune_dT.csv\", index_col=0, dtype=float)\n", "sigmaT = pd.read_csv(\"materials/hotNeptune_sigmaT.csv\", index_col=0, dtype=float)\n", "\n", "# now we can calculate the prior based on the model self-consistency\n", "prior = calc_prior(dT, sigmaT)" ] }, { "cell_type": "markdown", "id": "e79ea998", "metadata": {}, "source": [ "We can plot the model self-consistencies to see what they look like (the dotted line in the figure indicates where dT = sigmaT, i.e. the 1 $\\sigma$ region):" ] }, { "cell_type": "code", "execution_count": null, "id": "67f04b54", "metadata": {}, "outputs": [], "source": [ "plot_selfcons(\n", " dT,\n", " sigmaT,\n", " bounds_T0=(temp_min, temp_max),\n", " bounds_Mdot=(mdot_min, mdot_max),\n", " title=\"Model self-consistency\",\n", ")" ] }, { "cell_type": "markdown", "id": "996eda04", "metadata": {}, "source": [ "### Step 5: Combine likelihood and prior (to obtain _posterior_)\n", "\n", "Now that we have both the fit of the models to the equivalent width, and the constraints on the self-consistency of the parameter space, we can combine the two in a Bayesian framework to get a posterior distribution of the temperature and atmospheric mass-loss rate of the planet. As a reminder, Bayes' theorem reads as follows:\n", "\n", "$P(T_0,\\dot{M} | data) = \\dfrac{P(data | T_0, \\dot{M}) \\times P(T_0, \\dot{M})}{P(data)} = \\dfrac{P(data | T_0, \\dot{M}) \\times P(T_0, \\dot{M})}{\\int P(data | T_0, \\dot{M}) \\times P(T_0, \\dot{M}) \\; dT_0 \\; d\\dot{M}}$\n", "\n", "Here, $P$ denotes any probability quantity. $P(T_0,\\dot{M} | data)$ is the *posterior* distribution of $T_0$ and $\\dot{M}$ (given the data, i.e. our observed helium line). This is what we want to know. $P(data | T_0, \\dot{M})$ is the *likelihood* of the data, so this is what we get from fitting our data with a model as done in Step 3. $P(T_0, \\dot{M})$ is the *prior* expectation of how likely each model is (i.e. the self-consistency of the Parker wind models based on the temperature and not the data, which we calculated in Step 4). $P(data)$ is what is called the *(model) evidence* and is basically a normalization factor. Here we write it as the integral over all considered models, which assumes that the true values of $T_0$ and $\\dot{M}$ are within our considered parameter space. In our example, this assumption is probably fine, but if you run a scenario with a very limited extent of the $T_0-\\dot{M}$ parameter space, the assumption might break and the posterior that you get from this code will not be robust.\n", "\n", "The `calc_posterior()` function below evaluates Bayes' theorem and gives us the posterior distribution. It also calculates the marginalized posteriors, and prints the 1 $\\sigma$ constraints on $T_0$ and $\\dot{M}$. The `plot_posterior()` function plots the 2D posterior. Remember that for a 1D Gaussian, the 1, 2 and 3 $\\sigma$ contours are defined by the area that comprises 68%, 95% and 99.7% of the probability density function. However, for a 2D Gaussian, these values are 39%, 86% and 89.9% (see e.g. https://corner.readthedocs.io/en/latest/pages/sigmas/), and that is what the plotting function uses to draw the contours." ] }, { "cell_type": "code", "execution_count": null, "id": "d2f8579b", "metadata": {}, "outputs": [], "source": [ "# for the resolved line fit, we calculate the posterior with a flat prior\n", "# we do this only so that calc_posterior() prints the constrained T0 and Mdot based on\n", "# the likelihood alone, which we can then quote\n", "posterior_resolved = calc_posterior(1.0, likelihood_resolved)\n", "\n", "# for the EW fit, we calculate the posterior with the prior based on model self-consistency\n", "posterior_unresolved = calc_posterior(prior, likelihood_unresolved)" ] }, { "cell_type": "markdown", "id": "dcf1c0bc", "metadata": {}, "source": [ "We can plot the posteriors to see what they look like" ] }, { "cell_type": "code", "execution_count": null, "id": "44e9ab5c", "metadata": {}, "outputs": [], "source": [ "plot_posterior(\n", " posterior_resolved,\n", " bounds_T0=(temp_min, temp_max),\n", " bounds_Mdot=(mdot_min, mdot_max),\n", " title=\"Resolved fit posterior\",\n", ")\n", "\n", "plot_posterior(\n", " posterior_unresolved,\n", " bounds_T0=(temp_min, temp_max),\n", " bounds_Mdot=(mdot_min, mdot_max),\n", " title=\"Equivalent width fit posterior\",\n", ")" ] }, { "cell_type": "markdown", "id": "b1691601", "metadata": {}, "source": [ "### Results\n", "\n", "Finally, we plot the results in a way that reproduces Fig. 6 of Linssen et al. (2024)." ] }, { "cell_type": "code", "execution_count": null, "id": "950d3c60", "metadata": {}, "outputs": [], "source": [ "with plt.rc_context(\n", " {\"font.family\": \"serif\", \"mathtext.fontset\": \"dejavuserif\"}\n", "): # identical to paper\n", "\n", " fig, axes = plt.subplots(ncols=1, nrows=3, figsize=(4.5, 12))\n", "\n", " #######\n", "\n", " # re-calculate the best-fit spectrum\n", " bestfitsim = tools.Sim(\n", " tools.get_sunbather_project_path() + \"/sims/1D/hotNeptune/z_1/parker_4800_10.900/converged\"\n", " )\n", " bestspec, _, _ = RT.FinFout(bestfitsim, highres_wavs, \"He\")\n", "\n", " axes[0].plot(\n", " highres_wavs,\n", " highres_true_spec,\n", " color=\"blueviolet\",\n", " label=\"True model\",\n", " lw=2,\n", " zorder=-20,\n", " )\n", " axes[0].errorbar(\n", " wavs,\n", " observed_spec,\n", " yerr=observed_errorbar,\n", " fmt=\"o\",\n", " color=\"black\",\n", " label=\"Mock observation\",\n", " lw=2,\n", " zorder=-10,\n", " )\n", " axes[0].plot(\n", " highres_wavs,\n", " bestspec,\n", " color=\"dodgerblue\",\n", " label=\"Retrieved model\",\n", " linestyle=\"dashed\",\n", " lw=2,\n", " zorder=-5,\n", " )\n", " axes[0].text(\n", " 0.96, 0.04, r\"$\\bf{a)}$\", ha=\"right\", va=\"bottom\", transform=axes[0].transAxes\n", " )\n", " axes[0].set_xlabel(\"Wavelength [Å]\")\n", " axes[0].set_ylabel(r\"$F_{in}$/$F_{out}$\")\n", " axes[0].legend(loc=\"lower left\", edgecolor=\"none\", facecolor=\"none\", framealpha=0)\n", " axes[0].xaxis.set_major_formatter(\n", " matplotlib.ticker.FuncFormatter(lambda x, _: \"{:g}\".format(x))\n", " )\n", "\n", " #######\n", "\n", " plot_joint_constraint_resolved(\n", " None,\n", " None,\n", " chisqs,\n", " None,\n", " bounds_T0=(4000, 7000),\n", " bounds_Mdot=(10.5, 11.5),\n", " fig=fig,\n", " ax=axes[1],\n", " show_colorbar=False,\n", " )\n", " axes[1].scatter(\n", " 5100, 10.95, marker=\"x\", color=\"blueviolet\", s=100, zorder=100\n", " ) # the true model\n", " axes[1].text(4870, 11, \"True model\", color=\"blueviolet\", ha=\"left\", va=\"bottom\")\n", " axes[1].text(\n", " 4400,\n", " 10.75,\n", " \"Fit likelihood\",\n", " ha=\"left\",\n", " va=\"top\",\n", " color=plt.get_cmap(\"winter_r\")(0.65),\n", " )\n", " axes[1].text(\n", " 0.96,\n", " 0.96,\n", " r\"$\\bf{b)}$ Resolved line fit\",\n", " ha=\"right\",\n", " va=\"top\",\n", " transform=axes[1].transAxes,\n", " )\n", " axes[1].text(\n", " 4900,\n", " 10.887,\n", " r\"1$\\sigma$\",\n", " va=\"center\",\n", " ha=\"center\",\n", " color=plt.get_cmap(\"winter_r\")(0.55),\n", " fontsize=9,\n", " )\n", " axes[1].text(\n", " 4120,\n", " 10.765,\n", " r\"2$\\sigma$\",\n", " va=\"center\",\n", " ha=\"center\",\n", " color=plt.get_cmap(\"winter_r\")(0.45),\n", " fontsize=9,\n", " )\n", " axes[1].set_yticks(np.arange(10.6, 11.41, 0.1), minor=True)\n", "\n", " #######\n", "\n", " plot_joint_constraint_unresolved(\n", " dT,\n", " sigmaT,\n", " nsig_fit,\n", " posterior_unresolved,\n", " bounds_T0=(4000, 7000),\n", " bounds_Mdot=(10.5, 11.5),\n", " fig=fig,\n", " ax=axes[2],\n", " post_uptosigma=2,\n", " show_colorbar=False,\n", " )\n", " axes[2].scatter(\n", " 5100, 10.95, marker=\"x\", color=\"blueviolet\", s=100, zorder=100\n", " ) # the true model\n", " axes[2].text(4870, 11.05, \"True model\", color=\"blueviolet\", ha=\"left\", va=\"bottom\")\n", " axes[2].text(\n", " 4550, 11.4, \"Prior\", ha=\"left\", va=\"top\", color=plt.get_cmap(\"autumn\")(0.35)\n", " )\n", " axes[2].text(\n", " 6100,\n", " 10.86,\n", " \"Fit likelihood\",\n", " ha=\"left\",\n", " va=\"top\",\n", " color=plt.get_cmap(\"winter_r\")(0.65),\n", " )\n", " axes[2].text(4500, 10.73, r\"Posterior\", ha=\"left\", va=\"top\", color=\"black\")\n", " axes[2].text(\n", " 6300,\n", " 10.972,\n", " r\"1$\\sigma$\",\n", " va=\"center\",\n", " ha=\"center\",\n", " color=plt.get_cmap(\"winter_r\")(0.55),\n", " fontsize=9,\n", " )\n", " axes[2].text(\n", " 6300,\n", " 10.913,\n", " r\"2$\\sigma$\",\n", " va=\"center\",\n", " ha=\"center\",\n", " color=plt.get_cmap(\"winter_r\")(0.45),\n", " fontsize=9,\n", " )\n", " axes[2].text(\n", " 4680,\n", " 11.3,\n", " r\"1$\\sigma$\",\n", " va=\"center\",\n", " ha=\"center\",\n", " color=plt.get_cmap(\"autumn\")(0.45),\n", " fontsize=9,\n", " )\n", " axes[2].text(\n", " 5000,\n", " 11.3,\n", " r\"2$\\sigma$\",\n", " va=\"center\",\n", " ha=\"center\",\n", " color=plt.get_cmap(\"autumn\")(0.55),\n", " fontsize=9,\n", " )\n", " axes[2].text(\n", " 4850, 10.866, r\"1$\\sigma$\", va=\"center\", ha=\"center\", color=\"black\", fontsize=9\n", " )\n", " axes[2].text(\n", " 4500, 10.8, r\"2$\\sigma$\", va=\"center\", ha=\"center\", color=\"black\", fontsize=9\n", " )\n", " axes[2].text(\n", " 0.96,\n", " 0.96,\n", " r\"$\\bf{c)}$ Equivalent width fit\",\n", " ha=\"right\",\n", " va=\"top\",\n", " transform=axes[2].transAxes,\n", " )\n", " axes[2].set_yticks(np.arange(10.6, 11.41, 0.1), minor=True)\n", "\n", " #######\n", "\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "a1c75069", "metadata": {}, "source": [ "In the above plots, the orange contours represent the prior: the self-consistency of the Parker wind models. It presents the same quantity ($dT = T_{He} - T_0$) as the red-white-blue plots from Step 4. The orange line are models with $dT=0$, and the dark and light shaded areas represent the 1 and 2 $\\sigma_T$ contours (Eq. 5 of Linssen et al. 2022). The blue contours represent the likelihood: the fits of the model to the observed helium line. In panel b), the blue contours presents the likelihood converted into a posterior assuming a flat prior, showing the 1 and 2 $\\sigma$ regions. In panel c), the blue contours presents the same quantity ($\\Delta \\sigma$) as the plots from Step 3. The blue line are the models for which the EW is equal to that of the data, and the dark and light shaded areas represent 1 and 2 errorbars offset in EW. The black line represents the 1 $\\sigma$ contour of the posterior and can be understood as the final constraints on the temperature and mass-loss rate of the planet from the EW fit analysis.\n", "\n", "The plotting function also prints the constraints. You may note that there are small differences between these constraints and those printed between the `calc_posterior()` function. This is because the `calc_posterior()` function quotes the median and 1 $\\sigma$ bounds of the marginalized $T_0$ and $\\dot{M}$ posteriors, at the numerical precision of the parameter grid. The `plot_joint_constraint_unresolved()` function instead quotes the maximum and $\\sigma$ bounds of the 2D posterior surface, at the (interpolated and hence overestimated) numerical precision of the matplotlib contour plot." ] }, { "cell_type": "markdown", "id": "a10c1b00", "metadata": {}, "source": [ "## Retrieval of the observed helium spectrum of TOI-2134 b\n", "\n", "Having seen how a typical analysis is performed, we now move on to interpreting the helium spectrum of a real planet: TOI-2134 b. We perform the analysis for two different atmospheric compositions. The steps in principle are the same as before, so we do not repeat much of the explanations. As we are working with spectrally resolved data, we use a flat prior here. However, we also calculate the model self-consistency like we did before, simply to demonstrate that there is interesting physics going on here, since the self-consistent temperature is very different from that indicated by the data fit. So in this case we do not convert that self-consistency measure into a prior." ] }, { "cell_type": "markdown", "id": "db4cf422-94cd-4454-adfd-90810a8bf728", "metadata": {}, "source": [ "### Step 1: Create a grid of Parker wind profiles with p-winds\n", "\n", "\n", "This step can take a long time and is therefore skipped by default. We provide pre-run models that are used in this notebook. Uncomment the cell below if you want to run your own models.\n", "\n", "Here, when calculating the pure hydrogen/helium models, instead of passing the `fraction_hydrogen=0.91` argument which would result in using `p-winds` standalone, we use the `z=0` argument, so that _Cloudy_ is still involved to calculate a slightly more accurate mean particle mass than `p-winds` would calculate. The composition is the same, however." ] }, { "cell_type": "code", "execution_count": null, "id": "f20ef53e-c45a-45c8-9b32-d560ef4b014a", "metadata": {}, "outputs": [], "source": [ "# H/He profiles:\n", "t_min = 2000\n", "t_max = 8000\n", "t_step = 100\n", "temp = np.linspace(\n", " t_min,\n", " t_max,\n", " 1 + (t_max - t_min) // t_step\n", ")\n", "mdot_min = 8.5\n", "mdot_max = 10.0\n", "mdot_step = 0.02\n", "mdot = np.linspace(\n", " mdot_min,\n", " mdot_max,\n", " 1 + int((mdot_max - mdot_min) / mdot_step)\n", ")\n", "\n", "# construct_parker.run(\n", "# plname=\"TOI2134b\",\n", "# pdir=\"z_0\",\n", "# temp=temp,\n", "# mdot=mdot,\n", "# z=0,\n", "# no_tidal=True,\n", "# cores=8,\n", "# )\n", "\n", "# 100x solar metallicity profiles:\n", "\n", "# construct_parker.run(\n", "# plname=\"TOI2134b\",\n", "# pdir=\"z_100\",\n", "# temp=temp,\n", "# mdot=mdot,\n", "# z=100,\n", "# no_tidal=True,\n", "# cores=8,\n", "# )" ] }, { "cell_type": "markdown", "id": "5cc2ee1b", "metadata": {}, "source": [ "### Step 2: Run the Parker wind profiles through Cloudy\n", "\n", "This step can take a long time and is therefore skipped by default. We provide pre-run models that are used in this notebook. Uncomment the cell below if you want to run your own models." ] }, { "cell_type": "code", "execution_count": null, "id": "0fd7a53f", "metadata": {}, "outputs": [], "source": [ "# sunbather.convergeT_parker.run(\n", "# plname=\"TOI2134b\",\n", "# workingdir=\"z_0\",\n", "# pdir=\"z_0\",\n", "# temp=temp,\n", "# mdot=mdot,\n", "# z=0,\n", "# save_sp=[\"He\",],\n", "# cores=8\n", "# )\n", "\n", "# sunbather.convergeT_parker.run(\n", "# plname=\"TOI2134b\",\n", "# workingdir=\"z_100\",\n", "# pdir=\"z_100\",\n", "# temp=temp,\n", "# mdot=mdot,\n", "# z=100,\n", "# save_sp=[\"He\",],\n", "# cores=8\n", "# )" ] }, { "cell_type": "markdown", "id": "a4271b4d", "metadata": {}, "source": [ "### Step 3: Make transit spectra of the metastable helium triplet and fit the data (to obtain _likelihood_)\n", "\n", "We will again use a $\\chi^2$-statistic, which assumes that the errorbars on the data are normally distributed, and that the individual datapoints are independent. In the mock spectrum analysis, this was true by definition, since we independently drew the datapoints from a normal distribution, however for real data it is worth thinking about carefully!\n", "\n", "First, we plot the observed spectrum, which was downloaded from https://iopscience.iop.org/article/10.3847/2041-8213/aced51 . We then fit the data with a Gaussian to see if there is any offset from the rest-frame wavelength. Such offsets can definitely have a physical origin worth investigating, but as our model cannot reproduce offsets, here we shift the observed feature into the rest-frame wavelength before fitting it." ] }, { "cell_type": "code", "execution_count": null, "id": "4937b2d6", "metadata": {}, "outputs": [], "source": [ "def gaussian(x, amplitude, mean, std_dev):\n", " return amplitude * np.exp(-(((x - mean) / std_dev) ** 2) / 2)" ] }, { "cell_type": "code", "execution_count": null, "id": "4459f3c1", "metadata": {}, "outputs": [], "source": [ "observed_spec_toi2134 = pd.read_table(\n", " \"materials/dbf2.txt\",\n", " skiprows=13,\n", " names=[\"wav\", \"ea\", \"sig_ea\"],\n", " delim_whitespace=True,\n", ")\n", "observed_spec_toi2134.dropna(inplace=True)\n", "\n", "# fit a Gaussian to the peak of the data to obtain the line center\n", "params, covariance = curve_fit(\n", " gaussian,\n", " observed_spec_toi2134.wav[\n", " (observed_spec_toi2134.wav > 10832.5) & (observed_spec_toi2134.wav < 10833.8)\n", " ],\n", " observed_spec_toi2134.ea[\n", " (observed_spec_toi2134.wav > 10832.5) & (observed_spec_toi2134.wav < 10833.8)\n", " ],\n", " p0=[0.4, 10833.0, 0.5],\n", ")\n", "offset = params[1] - 10833.25\n", "print(\n", " \"There is a\",\n", " offset,\n", " \"angstrom shift of the line from rest-frame, which we now shift back.\",\n", ")\n", "observed_spec_toi2134.wav = observed_spec_toi2134.wav - offset\n", "\n", "\n", "# we will not fit the whole dataset, but only the part between 10831 and 10836 angstroms\n", "obs_wav_toi2134 = observed_spec_toi2134.wav[\n", " (observed_spec_toi2134.wav > 10831) & (observed_spec_toi2134.wav < 10836)\n", "].values\n", "obs_ea_toi2134 = observed_spec_toi2134.ea[\n", " (observed_spec_toi2134.wav > 10831) & (observed_spec_toi2134.wav < 10836)\n", "].values\n", "obs_sig_ea_toi2134 = observed_spec_toi2134.sig_ea[\n", " (observed_spec_toi2134.wav > 10831) & (observed_spec_toi2134.wav < 10836)\n", "].values\n", "\n", "\n", "fig, ax = plt.subplots(1)\n", "ax.errorbar(obs_wav_toi2134, obs_ea_toi2134, yerr=obs_sig_ea_toi2134)\n", "ax.axvline(10832.057472, color=\"k\", linestyle=\"dotted\")\n", "ax.axvline(10833.216751, color=\"k\", linestyle=\"dotted\")\n", "ax.axvline(10833.306444, color=\"k\", linestyle=\"dotted\")\n", "ax.set_xlabel(\"Wavelength [Å]\")\n", "ax.set_ylabel(\"Excess absorption [%]\")\n", "ax.set_title(\"TOI-2134 b observations by Zhang et al. 2023 (shifted)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "212c51d2", "metadata": {}, "source": [ "Uncomment the lines of the next cell and run them, if you ran your own models and want to calculate the EWs yourself, instead of using our pre-ran values." ] }, { "cell_type": "code", "execution_count": null, "id": "2a3b9cf2", "metadata": {}, "outputs": [], "source": [ "# chisqs_z0 = calc_chisqs_He10830(tools.get_sunbather_project_path()+'/sims/1D/TOI2134b/z_0/',\n", "# obs_wav_toi2134, obs_ea_toi2134, obs_sig_ea_toi2134,\n", "# instrument_R=32000, Mdotstep=0.02)\n", "# chisqs_z0.to_csv('materials/TOI2134b_chisqs_fit_z0.csv') #overwrites the supplied file\n", "#\n", "# chisqs_z100 = calc_chisqs_He10830(tools.get_sunbather_project_path()+'/sims/1D/TOI2134b/z_100/',\n", "# obs_wav_toi2134, obs_ea_toi2134, obs_sig_ea_toi2134,\n", "# instrument_R=32000, Mdotstep=0.02)\n", "# chisqs_z100.to_csv('materials/TOI2134b_chisqs_fit_z100.csv') #overwrites the supplied file" ] }, { "cell_type": "code", "execution_count": null, "id": "0ed5d4f3", "metadata": {}, "outputs": [], "source": [ "# read the chi-squared values from file\n", "chisqs_z0 = pd.read_csv(\n", " \"materials/TOI2134b_chisqs_fit_z0.csv\", index_col=0, dtype=float\n", ")\n", "\n", "chisqs_z100 = pd.read_csv(\n", " \"materials/TOI2134b_chisqs_fit_z100.csv\", index_col=0, dtype=float\n", ")\n", "\n", "# from chi-squared, we can calculate the likelihood\n", "likelihood_z0 = calc_likelihood_resolved(chisqs_z0)\n", "\n", "likelihood_z100 = calc_likelihood_resolved(chisqs_z100)" ] }, { "cell_type": "markdown", "id": "6acdabb3", "metadata": {}, "source": [ "Plot the chi-squared values of the pure H/He fit:" ] }, { "cell_type": "code", "execution_count": null, "id": "73648232", "metadata": {}, "outputs": [], "source": [ "plot_chisq_fit(\n", " chisqs_z0,\n", " bounds_T0=(2000, 8000),\n", " bounds_Mdot=(8.5, 10),\n", " title=\"TOI-2134 b, pure H/He\",\n", ")" ] }, { "cell_type": "markdown", "id": "0dd58ee4", "metadata": {}, "source": [ "The above figure may look weird. This is because we first ran a rough grid of models, using $\\dot{M}$ steps of 0.1 dex. This allowed us to see which region of parameter space fitted best. We then ran a finer grid around this region, using steps of 0.02 dex, resulting in the \"striped\" pattern. \n", "\n", "However, notice that on top of the striped pattern, there are some \"holes\" in the figure. These are failed models. The problem can be in _p-winds_, which sometimes fails to find a solution for specific combinations of $T_0$ and $\\dot{M}$. If this is the case, the isothermal Parker wind profile will not be present in the \\$SUNBATHER\\_PROJECT\\_PATH/parker\\_profiles/TOI2134b/z\\_0/ folder. The problem can also be a failed _Cloudy_ simulation, which can happen for example at high densities. If this is the case, you can go to the \\$CLOUDY\\_PATH/project/sims/1D/TOI2134b/z\\_0/parker\\_T0\\_Mdot/ folder to find out why _Cloudy_ failed. Finally, in some cases the `convergeT_parker.run()` function may fail to find a solution for the nonisothermal temperature structure. If this is the case, in the \\$SUNBATHER\\_PROJECT\\_PATH/sims/1D/TOI2134b/z\\_0/parker\\_T0\\_Mdot/ folder, you would find \"iteration\" files but no \"converged\" files. The algorithm may converge if you try the model again using a different starting criterion with the `startT` argument, or if you let it run for longer by specifying a higher maximum number of iterations with the `maxit` argument. In the case of this notebook, the missing models are due to failed convergence of the temperature structure. We are continuing to improve the *sunbather* convergence algorithm, aiming to be able to converge any model in the future." ] }, { "cell_type": "markdown", "id": "e983074d", "metadata": {}, "source": [ "### Step 4: Constraining the parameter space by assessing temperature self-consistency\n", "\n", "Although we will nog actually use this as a prior, we still calculate the model self-consistency for plotting and showing that it is different from the temperature indicated by the line width." ] }, { "cell_type": "markdown", "id": "ca1630ce", "metadata": {}, "source": [ "Uncomment the lines of the next cell and run them, if you ran your own models and want to calculate the model self-consistency yourself, instead of using our pre-ran values." ] }, { "cell_type": "code", "execution_count": null, "id": "60e3c65c", "metadata": {}, "outputs": [], "source": [ "# dT_z0, sigmaT_z0 = calc_dT_helium(tools.get_sunbather_project_path()+'/sims/1D/TOI2134b/z_0', Mdotstep=0.02)\n", "# dT_z0.to_csv('materials/TOI2134b_dT_z0.csv', float_format='%.3e') #overwrites the supplied files\n", "# sigmaT_z0.to_csv('materials/TOI2134b_sigmaT_z0.csv', float_format='%.3e')\n", "\n", "# dT_z100, sigmaT_z100 = calc_dT_helium(tools.get_sunbather_project_path()+'/sims/1D/TOI2134b/z_100', Mdotstep=0.02)\n", "# dT_z100.to_csv('materials/TOI2134b_dT_z100.csv', float_format='%.3e') #overwrites the supplied files\n", "# sigmaT_z100.to_csv('materials/TOI2134b_sigmaT_z100.csv', float_format='%.3e')" ] }, { "cell_type": "code", "execution_count": null, "id": "1087d9c1", "metadata": {}, "outputs": [], "source": [ "# load the dT and sigmaT values from file\n", "dT_z0 = pd.read_csv(\"materials/TOI2134b_dT_z0.csv\", index_col=0, dtype=float)\n", "sigmaT_z0 = pd.read_csv(\"materials/TOI2134b_sigmaT_z0.csv\", index_col=0, dtype=float)\n", "\n", "dT_z100 = pd.read_csv(\"materials/TOI2134b_dT_z100.csv\", index_col=0, dtype=float)\n", "sigmaT_z100 = pd.read_csv(\n", " \"materials/TOI2134b_sigmaT_z100.csv\", index_col=0, dtype=float\n", ")" ] }, { "cell_type": "markdown", "id": "7e5c9c04", "metadata": {}, "source": [ "### Results" ] }, { "cell_type": "code", "execution_count": null, "id": "18bfd116", "metadata": {}, "outputs": [], "source": [ "# we calculate the posterior with a flat prior\n", "# we do this only so that calc_posterior() prints the constrained T0 and Mdot based on\n", "# the likelihood alone, which we can then quote\n", "posterior_z0 = calc_posterior(1.0, likelihood_z0)\n", "\n", "posterior_z100 = calc_posterior(1.0, likelihood_z100)" ] }, { "cell_type": "markdown", "id": "a123fcdf", "metadata": {}, "source": [ "Finally, we plot the results in a way that reproduces Fig. 7 of Linssen et al. (2024). The plot also shows four spectra overplotted on the observations. So, we first need to calculate the spectra of these four different models again. If you ran your own grid of TOI-2134 b models in steps 1 and 2, you will already have these models and can readily calculate their spectra. If instead, you made use of our pre-calculated grids, the four specific models that we are plotting must be copied to the right folders in SUNBATHER\\_PROJECT\\_PATH. If you did not do this yet in the \"hotNeptune\" section above: We have supplied all necessary files in the /sunbather/examples/materials/ folder, you just need to copy it to the correct folder structure within your project path. That is; the /sunbather/examples/materials/parker\\_profiles folder must be copied/merged with your SUNBATHER\\_PROJECT\\_PATH/parker\\_profiles folder, and the /sunbather/examples/materials/sims folder must be copied/merged with your SUNBATHER\\_PROJECT\\_PATH/sims folder (if you run your own models, *sunbather* automatically creates these folder structures)." ] }, { "cell_type": "code", "execution_count": null, "id": "230fec93", "metadata": {}, "outputs": [], "source": [ "sim_z0_5000 = tools.Sim(\n", " tools.get_sunbather_project_path() + \"/sims/1D/TOI2134b/z_0/parker_5000_9.280/converged\"\n", ") # best-fit\n", "sim_z0_2600 = tools.Sim(\n", " tools.get_sunbather_project_path() + \"/sims/1D/TOI2134b/z_0/parker_2600_8.800/converged\"\n", ") # self-consistent T\n", "sim_z100_9700 = tools.Sim(\n", " tools.get_sunbather_project_path() + \"/sims/1D/TOI2134b/z_100/parker_9700_9.680/converged\"\n", ") # best-fit\n", "sim_z100_4000 = tools.Sim(\n", " tools.get_sunbather_project_path() + \"/sims/1D/TOI2134b/z_100/parker_4000_9.100/converged\"\n", ") # self-consistent T\n", "\n", "# calculate spectra\n", "wavsHe = np.linspace(10831, 10836, 200)\n", "spec_z0_5000, _, _ = RT.FinFout(sim_z0_5000, wavsHe, \"He\")\n", "spec_z0_2600, _, _ = RT.FinFout(sim_z0_2600, wavsHe, \"He\")\n", "spec_z100_9700, _, _ = RT.FinFout(sim_z100_9700, wavsHe, \"He\")\n", "spec_z100_4000, _, _ = RT.FinFout(sim_z100_4000, wavsHe, \"He\")\n", "\n", "# convert Fin/Fout to excess absorption in %\n", "ea_z0_5000 = (np.max(spec_z0_5000) - spec_z0_5000) * 100\n", "ea_z0_2600 = (np.max(spec_z0_2600) - spec_z0_2600) * 100\n", "ea_z100_9700 = (np.max(spec_z100_9700) - spec_z100_9700) * 100\n", "ea_z100_4000 = (np.max(spec_z100_4000) - spec_z100_4000) * 100\n", "\n", "# convolve to instrument resolution\n", "ea_z0_5000 = RT.convolve_spectrum_R(wavsHe, ea_z0_5000, 32000)\n", "ea_z0_2600 = RT.convolve_spectrum_R(wavsHe, ea_z0_2600, 32000)\n", "ea_z100_9700 = RT.convolve_spectrum_R(wavsHe, ea_z100_9700, 32000)\n", "ea_z100_4000 = RT.convolve_spectrum_R(wavsHe, ea_z100_4000, 32000)" ] }, { "cell_type": "code", "execution_count": null, "id": "edb4de21", "metadata": {}, "outputs": [], "source": [ "with plt.rc_context({\"font.family\": \"serif\", \"mathtext.fontset\": \"dejavuserif\"}):\n", " fig, axes = plt.subplots(ncols=1, nrows=3, figsize=(4.5, 12))\n", "\n", " ##########\n", "\n", " plot_joint_constraint_resolved(\n", " dT_z0,\n", " sigmaT_z0,\n", " chisqs_z0,\n", " None,\n", " bounds_T0=(2000, 12000),\n", " bounds_Mdot=(8.5, 10),\n", " fig=fig,\n", " ax=axes[0],\n", " show_colorbar=False,\n", " )\n", "\n", " axes[0].text(\n", " 0.96,\n", " 0.96,\n", " r\"$\\bf{a)}$\" + \" 91% H, 9% He\",\n", " ha=\"right\",\n", " va=\"top\",\n", " transform=axes[0].transAxes,\n", " )\n", " axes[0].text(\n", " 3600,\n", " 9.7,\n", " \"Temperature self-consistency\",\n", " ha=\"left\",\n", " va=\"top\",\n", " color=plt.get_cmap(\"autumn\")(0.35),\n", " )\n", " axes[0].text(\n", " 5400,\n", " 9.24,\n", " \"Fit likelihood\",\n", " ha=\"left\",\n", " va=\"top\",\n", " color=plt.get_cmap(\"winter_r\")(0.65),\n", " )\n", "\n", " ##########\n", "\n", " plot_joint_constraint_resolved(\n", " dT_z100,\n", " sigmaT_z100,\n", " chisqs_z100,\n", " None,\n", " bounds_T0=(2000, 12000),\n", " bounds_Mdot=(8.5, 10),\n", " fig=fig,\n", " ax=axes[1],\n", " show_colorbar=False,\n", " )\n", "\n", " axes[1].text(\n", " 0.96,\n", " 0.96,\n", " r\"$\\bf{b)}$\" + \" 100x solar metallicity\",\n", " ha=\"right\",\n", " va=\"top\",\n", " transform=axes[1].transAxes,\n", " )\n", " axes[1].text(\n", " 5100,\n", " 8.7,\n", " \"Temperature self-consistency\",\n", " ha=\"left\",\n", " va=\"top\",\n", " color=plt.get_cmap(\"autumn\")(0.35),\n", " )\n", " axes[1].text(\n", " 9000,\n", " 9.58,\n", " \"Fit likelihood\",\n", " ha=\"left\",\n", " va=\"top\",\n", " color=plt.get_cmap(\"winter_r\")(0.65),\n", " )\n", "\n", " ##########\n", "\n", " axes[2].errorbar(\n", " obs_wav_toi2134,\n", " obs_ea_toi2134,\n", " yerr=obs_sig_ea_toi2134,\n", " fmt=\"o\",\n", " color=\"black\",\n", " zorder=-100,\n", " lw=2,\n", " )\n", " axes[2].plot(\n", " wavsHe, ea_z0_5000, label=r\"H/He, $T_0$=5000\", color=\"dodgerblue\", lw=2\n", " )\n", " axes[2].plot(wavsHe, ea_z100_9700, label=r\"100x, $T_0$=9700\", color=\"crimson\", lw=2)\n", " axes[2].plot(\n", " wavsHe,\n", " ea_z0_2600,\n", " label=r\"H/He, $T_0$=2600\",\n", " color=\"dodgerblue\",\n", " linestyle=\"dashed\",\n", " lw=2,\n", " )\n", " axes[2].plot(\n", " wavsHe,\n", " ea_z100_4000,\n", " label=r\"100x, $T_0$=4000\",\n", " color=\"crimson\",\n", " linestyle=\"dashed\",\n", " lw=2,\n", " )\n", " axes[2].text(\n", " 0.96, 0.04, r\"$\\bf{c)}$\", ha=\"right\", va=\"bottom\", transform=axes[2].transAxes\n", " )\n", " axes[2].xaxis.set_major_formatter(\n", " matplotlib.ticker.FuncFormatter(lambda x, _: \"{:g}\".format(x))\n", " )\n", " axes[2].set_xlabel(\"Wavelength [Å]\")\n", " axes[2].set_ylabel(\"Excess absorption [%]\")\n", " axes[2].legend(\n", " loc=\"upper right\", edgecolor=\"none\", facecolor=\"none\", framealpha=0, ncols=2\n", " )\n", " axes[2].set_ylim(-0.15, 0.55)\n", "\n", " ##########\n", "\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "d44d5c35", "metadata": {}, "source": [ "We see that the temperature obtained from the line fit is very different from the self-consistency constraint. In Linssen et al. (2024), we discuss the temperature discrepancy between the line fit and the self-consistency in more detail." ] } ], "metadata": { "jupytext": { "formats": "ipynb,md:myst" }, "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.12.9" } }, "nbformat": 4, "nbformat_minor": 5 }