{ "cells": [ { "cell_type": "markdown", "id": "ecd5af06", "metadata": {}, "source": [ "# Custom Particle Distribution and Cloud Pressure Height\n", "\n", "Following the custom haze layer tutorial, this tutorial you will learn: \n", "\n", "1. How to generate your own Mie coefficients for a given radius grid\n", "2. How to compute the optics of an aerosol for your custom radius grid with an arbitrary fsed, and pressure height\n", "3. How to inject this aerosol as a custom layer into an atmospheric model for PICASO\n", "\n", "You need to have downloaded PICASO, Virga. This is particularly useful for retrievals with PICASO." ] }, { "cell_type": "markdown", "id": "c8bbfed5", "metadata": {}, "source": [ "First, let's import all the necessary packages:" ] }, { "cell_type": "code", "execution_count": null, "id": "c264e8ce", "metadata": {}, "outputs": [], "source": [ "#standards\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "#main programs\n", "\n", "#radiative transfer and atmosphere code\n", "from picaso import justdoit as pdi\n", "\n", "#cloud code\n", "from virga import justdoit as vdi\n", "\n", "\n", "#plotting tools\n", "from picaso import justplotit as ppi\n", "from virga import justplotit as vpi\n", "ppi.output_notebook()\n", "\n", "import matplotlib.pyplot as plt " ] }, { "cell_type": "markdown", "id": "f61a3074", "metadata": {}, "source": [ "Okay, now let's grab those Mieff coefficients we just calculated." ] }, { "cell_type": "code", "execution_count": null, "id": "15478225", "metadata": {}, "outputs": [], "source": [ "mieff_dir = '/data/virga'\n", "qext, qscat, cos_qscat, nwave, radius, wave_in = vdi.get_mie('SiO2',directory=mieff_dir)" ] }, { "cell_type": "markdown", "id": "8ea7cfb8", "metadata": {}, "source": [ "## Pick the Cloud Distribution Function\n", "\n", "### Gaussian particle distribution " ] }, { "cell_type": "code", "execution_count": null, "id": "edfd8920", "metadata": {}, "outputs": [], "source": [ "sigma = 1 #width of the distribution\n", "mu = -7 # mean particle size \n", "logradius = np.log10(radius)\n", "dist = (1/(sigma * np.sqrt(2 * np.pi)) *\n", " np.exp( - (logradius - mu)**2 / (2 * sigma**2)))\n", "plt.plot(logradius+4, dist) \n" ] }, { "cell_type": "markdown", "id": "255bad1c", "metadata": {}, "source": [ "Pick an approximate particle density (if using a fitting code this will be a free parameter) " ] }, { "cell_type": "code", "execution_count": null, "id": "caa16cfe", "metadata": {}, "outputs": [], "source": [ "ndz = 5e6 #that's particles/cm^2 but over the entire region of atmosphere where we're sticking aerosol\n", "\n", "opd,w0,g0,wavenumber_grid=vdi.calc_optics_user_r_dist(wave_in, ndz ,radius, pdi.u.cm,dist,\n", " qext, qscat, cos_qscat)" ] }, { "cell_type": "markdown", "id": "9e11ac2c", "metadata": {}, "source": [ "## Pick $f_{sed}$ and cloud base pressure" ] }, { "cell_type": "code", "execution_count": null, "id": "f0d75db8", "metadata": {}, "outputs": [], "source": [ "nlayer = 30\n", "pressure = np.logspace(-6,2,nlayer)\n", "z = np.linspace(100,0,nlayer)#just toy model here\n", "scale_h = 10 #toy model (could grab from picaso calc)\n", "\n", "base_pressure = 1e-1 # let's start by setting the haze layer base at 0.1 bars\n", "fsed=1 #this is more just an exponential scaling to control the cloud drop off, and is not connected to particle size\n", "\n", "opd_h = pressure*0+10\n", "opd_h[base_pressure=pressure]=opd_h[base_pressure>=pressure]*np.exp(\n", " -fsed*z[base_pressure>=pressure]/scale_h)\n", "opd_h = opd_h/np.max(opd_h)\n", "plt.loglog(opd_h, pressure)\n", "plt.ylim(1e2,1e-6)\n", "plt.xlim(1e-2,10)" ] }, { "cell_type": "code", "execution_count": null, "id": "8ac93062", "metadata": {}, "outputs": [], "source": [ "#here's where we shove all the variables into a dataframe PICASO can read\n", "df_cld = vdi.picaso_format_slab(base_pressure,opd, w0, g0, wavenumber_grid, pressure, \n", " p_decay=opd_h)" ] }, { "cell_type": "markdown", "id": "92464e9a", "metadata": {}, "source": [ "Let's checkout our optical parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "62bef42a", "metadata": {}, "outputs": [], "source": [ "nwno = len(wavenumber_grid) \n", "ppi.show(ppi.plot_cld_input(nwno, nlayer,df=df_cld,pressure=pressure, wavelength=1e4/wavenumber_grid))" ] }, { "cell_type": "markdown", "id": "eb8730f4", "metadata": {}, "source": [ "## Now run `PICASO`" ] }, { "cell_type": "code", "execution_count": null, "id": "b9a3b03a", "metadata": {}, "outputs": [], "source": [ "opa = pdi.opannection(wave_range=[0.3,14])\n", "case1 = pdi.inputs()\n", "\n", "case1.phase_angle(0)\n", "\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=pdi.u.Unit('M_jup'),\n", " radius=1.2, radius_unit=pdi.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 = pdi.u.Unit('R_sun') )\n", "\n", "#atmo -- make sure your pressure grid matches the one you computed your haze on!\n", "case1.atmosphere( df = pdi.pd.DataFrame({'pressure':np.logspace(-6,2,31),\n", " 'temperature':np.logspace(-9,3,31)*0+600, #just an isothermal one for simplicity\n", " \"H2\":np.logspace(-9,3,31)*0+0.837,\n", " \"He\":np.logspace(-9,3,31)*0+0.163,\n", " \"H2O\":np.logspace(-9,3,31)*0+1e-4}))" ] }, { "cell_type": "markdown", "id": "414674fd", "metadata": {}, "source": [ "And compute and plot our clear transmission spectrum, without our haze:" ] }, { "cell_type": "code", "execution_count": null, "id": "01f59184", "metadata": {}, "outputs": [], "source": [ "\n", "df= case1.spectrum(opa, full_output=True,calculation='transmission')\n", "\n", "wno, rprs2 = df['wavenumber'] , df['transit_depth']\n", "wno, rprs2 = pdi.mean_regrid(wno, rprs2, R=300)\n", "full_output = df['full_output']\n", "\n", "ppi.show(ppi.spectrum(wno,rprs2*1e6,plot_width=500))\n" ] }, { "cell_type": "markdown", "id": "25c50292", "metadata": {}, "source": [ "Now, let's add in that hazy information. We need to make sure the pressure grids match (by reducing the layers by one because of PICASO's idiosyncrasies)." ] }, { "cell_type": "code", "execution_count": null, "id": "8b8381fa", "metadata": {}, "outputs": [], "source": [ "case1.clouds(df=df_cld.astype(float))\n", "hazy= case1.spectrum(opa, full_output=True,calculation='transmission')\n", "hazyx,hazyy =hazy['wavenumber'] , hazy['transit_depth']\n", "hazyx,hazyy = pdi.mean_regrid(hazyx,hazyy, R=300)" ] }, { "cell_type": "markdown", "id": "b09f7a85", "metadata": {}, "source": [ "Let's compare our hazy spectrum to our clear spectrum!" ] }, { "cell_type": "code", "execution_count": null, "id": "22253f07", "metadata": { "scrolled": true }, "outputs": [], "source": [ "ppi.show(ppi.spectrum([wno,hazyx],[(rprs2*1e6),(hazyy*1e6)],\n", " legend=['clear','hazy'],plot_width=900,plot_height=300))" ] }, { "cell_type": "markdown", "id": "f2010c3b", "metadata": {}, "source": [ "Very cloud! Verification we can see the SiO2 feature at 9um" ] }, { "cell_type": "code", "execution_count": null, "id": "1998f849", "metadata": {}, "outputs": [], "source": [ "ppi.heatmap_taus(hazy)" ] } ], "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 } }, "nbformat": 4, "nbformat_minor": 5 }