{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Setup `PICASO`\n", "\n", "In this tutorial you will learn: \n", "\n", "1. What is a thermal emission spectrum\n", "2. How various atmospheric factors (temperature, abundances, clouds) influence an emission spectrum\n", "3. Given a spectrum, how do we analyze it's components \n", "\n", "What you should have already done:\n", "\n", "1. Complete all [Installation instructions](https://natashabatalha.github.io/picaso/installation.html) \n", " - This involves downloading two files, one of which is large (6 Gig). So plan accordingly! \n", "\n", "**Questions?** [Submit an Issue to PICASO Github](https://github.com/natashabatalha/picaso/issues) with any issues you are experiencing. Don't be shy! Others are likely experiencing similar problems\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings('ignore')\n", "import numpy as np\n", "import pandas as pd\n", "import os\n", "import astropy.units as u\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are the two main PICASO functions you will be exploring:\n", "\n", "`justdoit` contains all the spectroscopic modeling functionality you will need in these exercises.\n", "\n", "`justplotit` contains all the of the plotting functionality you will need in these exercises.\n", "\n", "Tips if you are not familiar with Python or `jupyter notebooks`:\n", "\n", "- Run a cell by clicking shift-enter. You can always go back and edit cells. But, make sure to rerun them if you edit it. You can check the order in which you have run your cells by looking at the bracket numbers (e.g. [1]) next to each cell.\n", "\n", "- In any cell you can write `help(INSERT_FUNCTION)` and it will give you documentation on the input/output\n", "\n", "- If you type `jdi.` followed by \"tab\" a box will pop up with all the available functions in `jdi`. This applies to any python function (e.g. `numpy`, `pandas`)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import picaso.justdoit as jdi\n", "import picaso.justplotit as jpi\n", "import picaso.opacity_factory as op\n", "jpi.output_notebook() #will force all our plots to appear in the notebook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you get an error regarding environment variables you can directly add them in the notebook. **You just need to make sure you run this line of code BEFORE the import the `picaso` functions.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#os.environ['picaso_refdata']='your_path'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Basic Inputs\n", "\n", "## Cross Section Connection\n", "\n", "As you will continue seeing throughout the ERS training sessions, all rapid radiative transfer codes rely on a database of pre-computed cross sections. Cross sections are computed by using line lists in combination with critical molecular pressure broadening parameters. Both can either be derived from theoretical first principles (e.g. [UCL ExoMol's line lists](https://www.exomol.com/)), measured in a lab, and/or some combination thereof (e.g. [HITRAN/HITEMP line lists](https://hitran.org/)). \n", "\n", "When cross sections are initially computed, a resolution ($\\lambda/\\Delta \\lambda$) is assumed. Cross sections are computed on a line-by-line nature and therefore usually computed for R~1e6. For JWST we are often interested in large bands (e.g. 1-14 $\\mu$m). Therefore we need creative ways to speed up these runs. You will usually find one of two methods: correlated-k tables, and resampled cross sections. [Garland et al. 2019](https://arxiv.org/pdf/1903.03997.pdf) is a good resource on the differences between these two. \n", "\n", "For this demonstration we will use the resampled cross section method. **The major thing to note about using resampled cross sections** is that you have to compute your model at ~100x higher resolution that your data. You will note that the opacity file you downloaded is resampled at R=10,000. Therefore you will note that **in this tutorial we will always bin it down to R=100**. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "opa = jdi.opannection(wave_range=[1,5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set Basic Planet and Stellar Inputs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Second step is to set basic planet parameters. To compute raw thermal flux, you only need gravity. However, if you want contrast units (relative flux of planet, to flux of star) you also need planet mass and radius, and steller radius. Below, we specify the planet's mass and radius. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "case1 = jdi.inputs()\n", "\n", "#here we are going to have to specify gravity through R and M since we need it in the Flux calc\n", "case1.gravity(mass=1, mass_unit=u.Unit('M_jup'), radius=1.2, radius_unit=u.Unit('R_jup'))\n", "\n", "#here we are going to have to specify R as well\n", "case1.star(opa, 4000,0.0122,4.437,radius=0.7, radius_unit = u.Unit('R_sun') )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Toy Models\n", "\n", "## How does climate structure influence resultant thermal flux?\n", "\n", "Before we load in full atmospheric profiles. Let's first gain intuition for how temperature-pressure profile influences your spectrum. \n", "\n", "Building a toy model: \n", "1. Select a pressure scale\n", " - More layers will increase the runtime of your code. The standard is ~50-90 levels.\n", " - It is critical that the range of pressure cover optically thick to optically thin regions. If you are concerned your range is too narrow, increase and test determine if it affects your spectrum. If it does, it is too narrow.\n", " - The standard is a pressure range is from 100 bar -1 microbar to be safe. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nlevels = 50\n", "pressure = np.logspace(-6,2,nlevels) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Specify vertical temperature profile that follows this scaling. We will try a few toy model examples: \n", " - isothermal \n", " - linearly increasing w/ pressure\n", " - linearly decreasing w/ pressure" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cases = ['iso','inc','dec']#will turn into the keys of our output\n", "min_t = 500\n", "max_t = 2000\n", "t_iso = np.zeros(nlevels) + 1300 #temperature in kelvin (totally arbitrary!)\n", "t_inc = np.linspace(min_t,max_t,nlevels)\n", "t_dec = np.linspace(max_t,min_t,nlevels)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Specify your atmospheric abundance breakdown. We will try these toy models: \n", " - \"well-mixed\" : equal abundance at all pressure levels \n", " \n", "For now, we will keep molecular abundances fixed, and vary the pressure-temperature profile. Let's put these all together." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "atmo_iso = {\n", " 'pressure':pressure,\n", " 'temperature':t_iso,\n", " 'H2': 0.80,\n", " 'He': 0.19,\n", " 'H2O':1e-3\n", "}\n", "atmo_inc = {\n", " 'pressure':pressure,\n", " 'temperature':t_inc,\n", " 'H2': 0.80,\n", " 'He': 0.19,\n", " 'H2O':1e-3\n", "}\n", "atmo_dec = {\n", " 'pressure':pressure,\n", " 'temperature':t_dec,\n", " 'H2': 0.80,\n", " 'He': 0.19,\n", " 'H2O':1e-3\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to loop through these three different cases. So we need to repeat the steps above three times, while changing the atmospheric input each time. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#let's store our three cases in toy_models\n", "toy_models={} \n", "\n", "#same as above but in a loop! three at once! \n", "for case, atmo in zip(cases,[atmo_iso,atmo_inc,atmo_dec ]):\n", " toy_models[case] = jdi.inputs()\n", " toy_models[case].gravity(mass=1, mass_unit=u.Unit('M_jup'), radius=1.2, radius_unit=u.Unit('R_jup'))\n", " toy_models[case].star(opa, 4000,0.0122,4.437,radius=0.7, radius_unit = u.Unit('R_sun') )\n", "\n", " #NOW we vary this to run our three different toy models\n", " toy_models[case].atmosphere(df = pd.DataFrame(atmo))\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally let's create three different toy model spectra " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "toy_out={}\n", "for case in toy_models.keys():\n", " #picaso also has transmission, and reflected light options but let's specify \n", " #thermal for this tutorial\n", " toy_out[case] = toy_models[case].spectrum(opa, calculation='thermal',full_output=True)\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally plot the sequence! " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wno,spec=[],[]\n", "fig = jpi.figure(height=400,width=500, y_axis_type='log',\n", " x_axis_label='Wavelength(um)',y_axis_label='Flux (erg/s/cm2/cm)')\n", "#for reference, let's plot our three pt profiles next door\n", "pt_fig = jpi.figure(height=400,width=300, y_axis_type='log',y_range=[1e2,1e-6],\n", " x_axis_label='Temperature',y_axis_label='Pressure(bar)')\n", "for i,case in enumerate(toy_out.keys()):\n", " x,y = jdi.mean_regrid(toy_out[case]['wavenumber'],\n", " toy_out[case]['thermal'], R=150)\n", "\n", " fig.line(1e4/x,y,color=jpi.Colorblind8[i],line_width=3,\n", " legend_label=case)\n", " \n", " pt_fig.line(toy_out[case]['full_output']['level']['temperature'],\n", " pressure,color=jpi.Colorblind8[i],line_width=3)\n", "fig.legend.location='bottom_right'\n", "\n", "jpi.show(jpi.row([fig,pt_fig]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What are the characteristic features of each spectrum? \n", "\n", "The first obvious feature is that our isothermal case follows a perfect blackbody. The second two might not be immediately noticeable if you are not familiar with the absorption cross section of H2O. Let's take a look." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How does absorption influence resultant thermal flux?\n", "\n", "### Molecular Absorption\n", "\n", "Our toy model included water absorption, along with H2/He. In order to see the cross section of H2O we can use `PICASO`'s opacity factory. This will show you the specific absorptive power of H2O at a specific pressure and temperature. Let's pick a single P and T that is representative of our toy models in order to gain an intuition for where water is absorbing" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#database should be the same for everyone\n", "db_filename = os.path.join(os.getenv('picaso_refdata'), 'opacities','opacities.db')\n", "species_to_get = ['H2O']\n", "t_to_get = [1500]#kelvin\n", "p_to_get = [1] #in bars\n", "data = op.get_molecular(db_filename, species_to_get, t_to_get,p_to_get)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x,y = jdi.mean_regrid(data['wavenumber'],data['H2O'][1500][1.0], R=150)\n", "#plot data\n", "h2o_fig = jpi.figure(height=300,y_axis_type='log',x_range=[1,5],y_range=[1e-24,5e-20]\n", " , x_axis_label='Micron',y_axis_label='Cross Section (cm2/species)')\n", "h2o_fig.line(1e4/x, y,line_width=4)\n", "jpi.show(h2o_fig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise:** This plot shows where the absorptive strength of water peaks. Find the peaks of water in the cross section plot and correlate them with the peaks and troughs of your toy model. \n", "\n", "In the case where temperature is increasing with increasing pressure, does the flux or your spectrum increase toward the peak of an absorption feature? And vice versa? What does this tell you about where the flux is emanating from at band center, band trough? \n", "\n", "### Continuum Absorption\n", "\n", "A second, more subtle contributor to the spectrum is the continuum absorption from H2, He, which we included. At first glance it might look like H2O contribution is the only contributor to our toy models. Let's take a closer look at one case by using the `get_contribution` function. \n", "\n", "### Get contribution function\n", "\n", "This output consists of three important items: \n", "\n", "`taus_per_layer`\n", "- Each dictionary entry is a nlayer x nwave that represents the per layer optical depth for that molecule. \n", "\n", "`cumsum_taus`\n", "- Each dictionary entry is a nlevel x nwave that represents the cumulative summed opacity for that molecule. \n", "\n", "`tau_p_surface` \n", "- Each dictionary entry is a nwave array that represents the pressure level where the cumulative opacity reaches the value specified by the user through `at_tau`.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "contribution = jdi.get_contribution(toy_models['inc'], opa, at_tau=1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#explore the output\n", "contribution['tau_p_surface'].keys()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wno=[]\n", "spec=[]\n", "labels=[]\n", "for j in contribution['tau_p_surface'].keys():\n", " x,y = jdi.mean_regrid(opa.wno, contribution['tau_p_surface'][j],R=100)\n", " if np.min(y)<5: # Bars\n", " wno+=[x]\n", " spec+=[y]\n", " labels +=[j]\n", "fig = jpi.spectrum(wno,spec,plot_width=600,plot_height=350,\n", " y_axis_label='Tau~1 Pressure (bars)',\n", " y_axis_type='log',x_range=[1,5],\n", " y_range=[1e2,1e-4],legend=labels)\n", "jpi.show(fig)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#this plot does the same thing and we will use it from now on\n", "jpi.show(jpi.molecule_contribution(contribution, opa,\n", " min_pressure=4.5,\n", " R=100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This plot shows where the optical depth=1 (i.e. tau=1) surface is for the different absorbing components of your model. In this case, the tau=1 surface of H2O is far above that of the continuum. Therefore, for the toy model we have specified, the continuum is not contributing to the overall spectrum. \n", "\n", "**Exercise:** In this toy model example, at what abundance of H2O does the continuum become important? In those cases, what is the result on the final spectrum? Run the abundance cases below to find out." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#let's store our three cases in toy_models\n", "vary_h2o_out={} \n", "\n", "#same as above but in a loop over a h2o abundance multilier \n", "h2o_factors = [100, 10, 1, 1e-1, 1e-2]\n", "for ifact in h2o_factors:\n", " pln = jdi.inputs()\n", " pln.gravity(mass=1, mass_unit=u.Unit('M_jup'), radius=1.2, radius_unit=u.Unit('R_jup'))\n", " pln.star(opa, 4000,0.0122,4.437,radius=0.7, radius_unit = u.Unit('R_sun') )\n", "\n", " #NOW we vary this to run our three different toy models\n", " pln.atmosphere(df = pd.DataFrame(atmo_inc))\n", " \n", " #to figure out the problem above, we can artificially decrease/increase H2O\n", " pln.inputs['atmosphere']['profile']['H2O'] = (ifact * \n", " pln.inputs['atmosphere']['profile']['H2O'])\n", " #and rerun the spectrum with those different abundances \n", " vary_h2o_out[ifact] = pln.spectrum(opa, calculation='thermal', full_output=True)\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally plot!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wno,spec=[],[]\n", "fig = jpi.figure(height=400,width=500, y_axis_type='log',\n", " x_axis_label='Wavelength(um)',y_axis_label='Flux (erg/s/cm2/cm)',\n", " title='Toy model with varied H2O abundance')\n", "\n", "#can create a little color scale for each of our spectra\n", "colors = jpi.pals.viridis(len(h2o_factors))\n", "for i,case in enumerate(vary_h2o_out.keys()):\n", " x,y = jdi.mean_regrid(vary_h2o_out[case]['wavenumber'],\n", " vary_h2o_out[case]['thermal'], R=150)\n", "\n", " fig.line(1e4/x,y,color=colors[i],line_width=3,\n", " legend_label=f'{case}xH2O')\n", "\n", "fig.legend.location='bottom_right'\n", "\n", "jpi.show(fig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two factors influencing the spectrum. The first is the increased opacity of the water itself. Why does flux increase with decreasing water abundance? \n", "\n", "The second, less noticeable is the new contribution from the continuum opacity. Can you eyeball where the continuum kicks in? One trick is to normalize each spectrum by its mean value to compare the relative size of the absorption features. This isolates the size of the absorption features while removing the increased flux from decreased water abundance. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wno,spec=[],[]\n", "fig = jpi.figure(height=400,width=500, #y_axis_type='log',\n", " x_axis_label='Wavelength(um)',y_axis_label='Flux (erg/s/cm2/cm)',\n", " title='Normalized toy model with varied H2O abundance')\n", "\n", "#can create a little color scale for each of our spectra\n", "colors = jpi.pals.viridis(len(h2o_factors))\n", "for i,case in enumerate(vary_h2o_out.keys()):\n", " x,y = jdi.mean_regrid(vary_h2o_out[case]['wavenumber'],\n", " vary_h2o_out[case]['thermal'], R=150)\n", " spec += [y]\n", " fig.line(1e4/x,y/np.mean(y),color=colors[i],line_width=3,\n", " legend_label=str(case))\n", "\n", "fig.legend.location='bottom_right'\n", "\n", "jpi.show(fig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now can you tell? When water is fully dominant (dark purple) the delta peak-trough flux is the highest. When the continuum is added, the H2H2/H2He opacity interferes with the H2O absorption. The window regions, which would otherwise be optically thin and sensitive to the highest pressures, are now blocked by the continuum opacity. This is an important effect as it will come back to haunt us in the cloud exercises. Let's take a look at the contribution plot for our 0.01xH2O case, to see how H2H2 and H2He now play a role." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#PLN is from the last run in the cell above\n", "contribution = jdi.get_contribution(pln, opa, at_tau=1)\n", "\n", "jpi.show(jpi.molecule_contribution(contribution, opa,\n", " min_pressure=4.5,\n", " R=100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ah ha! We see in this case the continuum and the molecular absorption from H2H2 and H2He now compete across wavelength space. \n", "\n", "**Confirm understanding:** Does this make sense with what you noted from the raw and normalized spectra? What is happening to the raw spectra in the regions most dominated by continuum opacity?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to determine what pressures you are sensitive to?\n", "\n", "One major aspect of understanding thermal emission is understanding where (in pressure) your flux is emanating from. The tau=1 surface plots should give you a rough idea of this. You can see that in our 0.01xH2O case above, the flux is coming from roughly 1 bars. \n", "\n", "Another way to visualize this, and relate it back to your spectrum is by comparing your raw spectra against blackbody curves. What is most helpful is to pick temperatures at known pressures along your climate profile.\n", "\n", "**Exercise:** Revisit the Pressure-temperature profile for the `inc` case. Even though we are increasing water abundance, all those cases still are computed from the same pressure-temperature profile. Pick three pressures along this curve. Determine what the temperature is at those pressures. Use the ``blackbody`` function to compute three blackbodies and compare against your thermal flux spectra. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wno,spec=[],[]\n", "fig = jpi.figure(height=400,width=500, y_axis_type='log',\n", " x_axis_label='Wavelength(um)',y_axis_label='Flux (erg/s/cm2/cm)',\n", " title='Toy model with varied H2O abundance')\n", "\n", "#for reference, let's plot our pt profile next door\n", "pt_fig = jpi.figure(height=400,width=300, y_axis_type='log',y_range=[1e2,1e-6],\n", " x_axis_label='Temperature',y_axis_label='Pressure(bar)')\n", "pt_fig.line(toy_out['inc']['full_output']['level']['temperature'],\n", " pressure,color='black',line_width=3)\n", "\n", "#same exact code as before\n", "colors = jpi.pals.viridis(len(h2o_factors))\n", "for i,case in enumerate(vary_h2o_out.keys()):\n", " x,y = jdi.mean_regrid(vary_h2o_out[case]['wavenumber'],\n", " vary_h2o_out[case]['thermal'], R=150)\n", "\n", " fig.line(1e4/x,y,color=colors[i],line_width=3,\n", " legend_label=f'{case}xH2O')\n", "#show for reference first\n", "jpi.show(jpi.row([fig, pt_fig])) \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#STEPS FOR EXERCISE\n", "#Step 1: pick a pressure\n", "at_pressure=1#bar\n", "#Step 2: what is the temperature of your planet at that pressure\n", "corresponding_t = 1630 #Kelvin, I have just eyeballed this from the plot \n", "#Step 3: use the plack function to compute the blackbody flux\n", "corr_intensity = jpi.blackbody(corresponding_t, 1/opa.wno)[0]\n", "corr_flux = np.pi * corr_intensity\n", "#Step4: add to your plots\n", "pt_fig.circle(corresponding_t, at_pressure, size=10, color='black')\n", "fig.line(1e4/opa.wno, corr_flux, color='black',line_width=4)\n", "jpi.show(jpi.row([fig, pt_fig])) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Does this checkout out with your tau=1 pressure curves for the 0.01xH2O case you made above? What about the increased water abundance plots? \n", "\n", "**Exercise:** Using this methodology, for each of these 1-5 micron spectra, determine the range of pressures your spectrum is sensitive to? \n", "\n", "Now that you have gone through the exercise, you can use the prebuilt `PICASO` function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f=jpi.flux_at_top(vary_h2o_out[1], pressures=[1,0.1,0.01],R=150)\n", "f.legend.location='bottom_right'\n", "jpi.show(f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Parameterized Pressure-Temperature Profiles\n", "\n", "In lecture you learned about the double gray model, which is an exoplanet-relevant analytic solution of the two-stream equation [Guillot et al. 2010](#References). You also encountered this in last week's transmission tutorial. The basic equation is: \n", "\n", "$$ T = \\left(\\frac{F_\\odot}{2 \\sigma} \\left[(\\frac{1}{\\gamma} - \\gamma)\\exp^{-\\tau/\\gamma} + 1 + \\gamma \\right] \\right)^{1/4} $$\n", "\n", "Where the limits are: \n", "\n", "- $\\gamma>>1$ : \"greenhouse\" limit with a hot deep atmosphere\n", "- $\\gamma=1$ : isothermal atmosphere with T $= \\frac{F_\\odot}{\\sigma}^{1/4}$\n", "- $\\gamma<<1$ : \"anti-greenhouse\" limit, thermal inversion" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pt_fig = jpi.figure(height=400,width=300, y_axis_type='log',y_range=[1e2,1e-6],\n", " x_axis_label='Temperature',y_axis_label='Pressure(bar)')\n", "#first, tau\n", "#let's assume the optical depth logarithmically increases from 1e-5-1000\n", "tau =np.logspace(-5,3,50)\n", "#this will roughly follow our pressure scale \n", "pressure = np.logspace(-6,2,50)\n", "\n", "F_sig = 1500**4 #we can play around with this insolation scaling\n", "\n", "#gamma \n", "gamma_gt_1 = 10\n", "gamma_lt_1 = 0.1\n", "gamma_1 = 1\n", "for i,ig in enumerate(zip([gamma_1,gamma_lt_1, gamma_gt_1],['g=1','g<1','g>1'])):\n", " g = ig[0]\n", " legend=ig[1]\n", " \n", " temperature = (F_sig * ((1/g -g)*np.exp(-tau/g) +1 + g)) **(0.25)\n", " pt_fig.line(temperature,\n", " pressure,color=jpi.Colorblind8[i],line_width=3,legend_label=legend)\n", "\n", "jpi.show(pt_fig)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Make sure you have an intuition for how the parameters gamma, tau, and Fsig affect the resulting PT profile. \n", "\n", "We've already learned that the isothermal profile will return back the blackbody. Let's use one of the \"greenhouse\" limit pressure-temperature profiles to proceed with creating a full thermal emission spectrum. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Combing parameterized climate with chemistry" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "case1 = jdi.inputs()\n", "case1.gravity(mass=1, mass_unit=u.Unit('M_jup'), radius=1.2, radius_unit=u.Unit('R_jup'))\n", "case1.star(opa, 4000,0.0122,4.437,radius=0.7, radius_unit = u.Unit('R_sun') )\n", "\n", "#NOW let's add our parameterized profile\n", "case1.atmosphere(df = pd.DataFrame({\n", " 'pressure':pressure,\n", " 'temperature':temperature}),verbose=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need to add the chemistry! `PICASO` has a prebuilt chemistry table that was computed by Channon Visscher. You can use it by adding it to your `case1`. Two more chemistry parameters are now going to be introduced: \n", "\n", "1. C/O ratio: Elemental carbon to oxygen ratio \n", "2. M/H: Atmospheric metallicity \n", "\n", "Let's choose Solar values. Feel free to explore the effect of this after the first completed spectrum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "log_mh = 0 #log relative to solar\n", "c_o = 1 #relative to solar\n", "case1.chemeq_visscher( c_o, log_mh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now you can check out what has been added to your `case1` bundle " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "case1.inputs['atmosphere']['profile'].head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to compute the spectrum" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "out_solar = case1.spectrum(opa, calculation='thermal',full_output=True)\n", "#get the contribution as well now that we have all the chemistry!\n", "contribution = jdi.get_contribution(case1, opa, at_tau=1)\n", "#regrid\n", "wno, fp = jdi.mean_regrid(out_solar['wavenumber'], out_solar['thermal'], R=150)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jpi.show(jpi.spectrum(wno, fp, y_axis_type='log',\n", " y_axis_label='Flux (erg/cm2/s/cm)',plot_width=400))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When comparing to transit observations we will mostly be looking at contrast units, relative to the stellar flux. Note the 1e6 multiplier gets us to PPM units. This allows you to orient your brain to an observers. Remember the hypothesized JWST noise from is 20 & 50 ppm for near-IR and mid-IR modes, respectively. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wno, fpfs = jdi.mean_regrid(out_solar['wavenumber'], out_solar['fpfs_thermal'], R=150)\n", "jpi.show(jpi.spectrum(wno, fpfs*1e6, \n", " y_axis_label='Relative Flux (ppm)',plot_width=400))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can revisit our contribution plot, which will look more complicated now." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jpi.show(jpi.molecule_contribution(contribution, opa,\n", " min_pressure=4.5,\n", " R=100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: \n", "\n", "1. Cross compare this contribution plot with your resultant spectrum. Can you properly account for all the molecules that you should be dominant. \n", "2. What molecules are barely visible. What is the approximate signal size of those features?\n", "3. Look through the JWST modes from the figure on the [PandExo webiste](https://exoctk.stsci.edu/pandexo/calculation/new). For the transit time series modes you see, what molecules are observable with what JWST models. \n", "4. Repeat the exercise such that your upper atmospheric temperature is ~600 K. \n", " - What major differences do you notice in your contribution plot? \n", " - What are the dominant carbon-bearing species?\n", " - What are the dominant nitrogen-bearing species? \n", " - What are the dominant continuum species?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Combing parameterized climate with chemistry AND clouds\n", "\n", "In this last module we will think about how clouds affect your thermal emission spectrum. We will use the same general procedure outlined above but add one additional step to add a box model cloud. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cld1 = jdi.inputs()\n", "cld1.gravity(mass=1, mass_unit=u.Unit('M_jup'), radius=1.2, radius_unit=u.Unit('R_jup'))\n", "cld1.star(opa, 4000,0.0122,4.437,radius=0.7, radius_unit = u.Unit('R_sun') )\n", "\n", "#NOW let's add our parameterized profile\n", "cld1.atmosphere(df = pd.DataFrame({\n", " 'pressure':pressure,\n", " 'temperature':temperature}),verbose=False)\n", "log_mh = 0 #log relative to solar\n", "c_o = 1 #relative to solar\n", "cld1.chemeq_visscher( c_o, log_mh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding a box model cloud \n", "\n", "Later in the ERS series, you will learn how to compute full cloud models. For now, we will use a simple box model cloud to understand the effect of adding a cloud. \n", "\n", "The `PICASO` box model is specified by a cloud layer with an asymmetry parameter (g0), a single scattering abledo (w0), an optical depth (opd) and a vertical pressure location (p,the pressure level in log10 bars) and finally the vertical cloud thickness (dp, the cloud thickness also in log10 bars). Such that:\n", "\n", "cloud_base(bars)=$10^p$\n", "\n", "cloud_top(bars)=$10^{p−dp}$\n", "\n", "The single scattering albedo controls how scattering the cloud is. The asymmetry controls the degree of forward scattering. Checkout the `PICASO` [radiative transfer tutorial](https://natashabatalha.github.io/picaso_dev#slide02) to see a visual of the asymmetry phase function. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# highly forward scattering cloud from 1.0 bar up to 0.1 bar\n", "cld1.clouds( g0=[0.9], w0=[0.8], opd=[0.5], p = [0.0], dp=[1.0]) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use the cloud input function to visualize what we just added to our code" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nwno = 196 #this is just the default number for the simple case above\n", "nlayer = cld1.nlevel-1 #one less than the number of PT points in your input\n", "jpi.show(jpi.plot_cld_input(nwno, nlayer,df=cld1.inputs['clouds']['profile']))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's similarly compute the spectrum and compare to our cloud free case" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "out_cld = cld1.spectrum(opa, calculation='thermal',full_output=True)\n", "#get the contribution as well now that we have all the chemistry!\n", "contribution_cld = jdi.get_contribution(cld1, opa, at_tau=1)\n", "#regrid\n", "wno, fp_cld = jdi.mean_regrid(out_cld['wavenumber'], out_cld['thermal'], R=150)\n", "wno, fpfs_cld = jdi.mean_regrid(out_cld['wavenumber'], out_cld['fpfs_thermal'], R=150)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jpi.show(jpi.spectrum([wno,wno], [fpfs*1e6, fpfs_cld*1e6], legend=['Cloud free','Cloudy'],\n", " y_axis_label='Relative Flux (ppm)',plot_width=600))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looks relatively minor! Why is this? Let's see the contribution plot with the cloud to find out" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jpi.show(jpi.molecule_contribution(contribution_cld, opa,\n", " min_pressure=4.5,\n", " R=100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do the minor modulations that you see in your cloudy thermal emission spectrum make sense with there the cloud tau=1 surface is? \n", "\n", "**Final Exercise:** Return to where we defined the box model. Increase the cloud thickness until you can see it in the contribution plot. \n", "\n", "1. What does your spectrum approach in the 100% cloud coverage? \n", "2. What spectral features are first made undetectable because of clouds? \n", "3. What spectral features are least inhibited by cloud coverage?\n", "4. What JWST spectral models in 1-5 micron region are most susceptible to cloud coverage? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# References \n", "\n", "[Guillot, Tristan. \"On the radiative equilibrium of irradiated planetary atmospheres.\" Astronomy & Astrophysics 520 (2010): A27.](https://arxiv.org/abs/1006.4702)" ] } ], "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" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }