Note

Phase Curves Part 1¶

Citation: Robbins-Blanch et al. 2022 “Cloudy and Cloud-free Thermal Phase Curves with PICASO: Applications to WASP-43b” ApJ

From the previous two tutorials you should understand how:

How to convert your GCM input to PICASO’s required xarray
How to post-process output to append to your GCM output

Here you will learn:

How to compute a phase curve
How to analyze the resulting output

[1]:

import pandas as pd
import numpy as np

from picaso import justdoit as jdi
from picaso import justplotit as jpi
jpi.output_notebook()

Loading BokehJS ...

Analyze Phase Curves¶

All original plotting tools still work by selecting an individual phase

[12]:

#same old same old
wno =[];fpfs=[];legend=[]
for iphase in allout.keys():
    w,f = jdi.mean_regrid(allout[iphase]['wavenumber'],
                               allout[iphase]['fpfs_thermal'],R=100)
    wno+=[w]
    fpfs+=[f*1e6]
    legend +=[str(int(iphase*180/np.pi))]
jpi.show(jpi.spectrum(wno, fpfs, plot_width=500,legend=legend,
                     palette=jpi.pals.viridis(n_phases)))

Phase Curve Plotting Tools: Phase Snaps¶

Phase snaps allows you to show snapshots of what is happening at each phase that was computed. x,y,z axes can be swapped out for flexibility.

Allowable x: 'longitude', 'latitude', 'pressure'
Allowable y: 'longitude', 'latitude', 'pressure'
Allowable z: 'temperature','taugas','taucld','tauray','w0','g0','opd' (the latter three are cloud properties)
Allowable collapse: Collapse tells picaso what to do with the extra axes. For instance, taugas is an array that is [npressure,nwavelength,nlongitude,nlatitude]. If we want pressure x latitude we have to collapse wavelength and the longitude axis. To do so we could either supply an integer value to select a single dimension. Or we could supply one of the following as str input: [np.mean, np.median, np.min, np.max]. If there are multiple axes you want to collapse, supply an ordered list.

[13]:

fig=jpi.phase_snaps(allout, x='longitude', y='pressure',z='temperature',
                y_log=True, z_log=False, collapse='np.mean')#this will collapse the latitude axis by taking a mean

../_images/notebooks_9e_PhaseCurves_22_0.png

[14]:

fig=jpi.phase_snaps(allout, x='longitude', y='pressure',z='taugas', palette='PuBu',
                y_log=True, z_log=True, collapse=['np.mean','np.mean'])
#this will collapse the latitude axis by taking a mean, and the wavelength axis also by taking the mean

../_images/notebooks_9e_PhaseCurves_23_0.png

[15]:

fig=jpi.phase_snaps(allout, x='longitude', y='latitude',z='taugas', palette='PuBu',
                y_log=False, z_log=True, collapse=[20,'np.mean'])
#this will collapse the pressure axis by taking the 20th pressure grid point
#and the wavelength axis by taking the mean
print('Pressure at', allout[iphase]['full_output']['layer']['pressure'][20,0,0])

Pressure at 0.0030016929695969906

../_images/notebooks_9e_PhaseCurves_24_1.png

Phase Curve Plotting Tools: Phase Curves¶

We can use the same collapse feature as in the phase_snaps tool to collapse the wavelength axis, when plotting phase curves. For example, we can select a single wavelength point at a given resolution, or we can average over all wavelengths.

Allowable collapse: - 'np.mean' or np.sum - float or list of float: wavelength in microns (will find the nearest value to this wavelength). Must be in wavenumber range.

For float option, user must specify resolution.

[16]:

to_plot = 'fpfs_thermal'#or thermal
collapse = [1.1,1.4,1.6]#micron
R=100
phases, all_curves, all_ws, fig=jpi.phase_curve(allout, to_plot, collapse=collapse, R=100)
jpi.show(fig)

Run Cloudy Thermal Phase Curve w/ 3D Input¶

Here we will create a hypothetical cloud path and add it to our fake H2O profile to demonstrate the basics of running a cloudy phase curve

[17]:

# create coords
lon = all_gcm.coords.get('lon').values
lat = all_gcm.coords.get('lat').values
pres = all_gcm.coords.get('pressure').values
pres_layer = np.sqrt(pres[0:-1]*pres[1:])
wno_grid = np.linspace(1e4/2,1e4/1,10)#cloud properties are defined on a wavenumber grid

#create box-band cloud model
fake_opd = np.zeros((len(lon), len(lat),len(pres_layer), len(wno_grid))) # create fake data
where_lat = np.where(((lat>-60) & (lat<60)))#creating a grey cloud path
where_lon = np.where(lon>0)#creating a grey cloud path
where_pres = np.where(((pres_layer<0.01) & (pres_layer>1e-4)))#creating a grey cloud band

for il in where_lat[0]:
    for ip in where_pres[0]:
        for ilo in where_lon[0]:
            fake_opd[ilo,il,ip,:]=10 #optical depth of 10 (>>1)

#make up asymmetry and single scattering properties
fake_asymmetry_g0 = 0.8+ np.zeros((len(lon), len(lat),len(pres_layer), len(wno_grid)))
fake_ssa_w0 = 0.9+ np.zeros((len(lon), len(lat),len(pres_layer), len(wno_grid)))

# put data into a dataset
ds_cld= jdi.xr.Dataset(
    data_vars=dict(
        opd=(["lon", "lat","pressure","wno"], fake_opd,{'units': 'depth per layer'}),
        g0=(["lon", "lat","pressure","wno"], fake_asymmetry_g0,{'units': 'none'}),
        w0=(["lon", "lat","pressure","wno"], fake_ssa_w0,{'units': 'none'}),
    ),
    coords=dict(
        lon=(["lon"], lon,{'units': 'degrees'}),#required
        lat=(["lat"], lat,{'units': 'degrees'}),#required
        pressure=(["pressure"], pres_layer,{'units': 'bar'}),#required
        wno=(["wno"], wno_grid,{'units': 'cm^(-1)'})#required for clouds
    ),
    attrs=dict(description="coords with vectors"),
)
ds_cld['opd'].isel(pressure=30,wno=0).plot(x='lon',y='lat')

[17]:

<matplotlib.collections.QuadMesh at 0x7f8047995e90>

../_images/notebooks_9e_PhaseCurves_28_1.png

[18]:

case_cld = jdi.inputs()
#setup geometry of run
n_phases = 4
min_phase = 0
max_phase = 2*np.pi
phase_grid = np.linspace(min_phase,max_phase,n_phases)#between 0 - 2pi
#send params to phase angle routine
case_cld.phase_angle(phase_grid=phase_grid,
                    num_gangle=6, num_tangle=6,calculation='thermal')

#setup 4d atmosphere
case_cld.atmosphere_4d(all_gcm, shift = np.zeros(n_phases), zero_point='night_transit',
                                     plot=True,verbose=False)
#no need to input shift here, it will take it from atmosphere_4d
#in fact, clouds_4d must always be run AFTER atmosphere_4d
case_cld.clouds_4d(ds_cld, iz_plot=30,
                                     plot=True,verbose=False)
#same old planet and star properties
case_cld.gravity(radius=1,radius_unit=jdi.u.Unit('R_jup'),
                mass=1, mass_unit=jdi.u.Unit('M_jup')) #any astropy units available
case_cld.star(opacity,5000,0,4.0, radius=1, radius_unit=jdi.u.Unit('R_sun'))

/home/nbatalh1/anaconda3/envs/pic311/lib/python3.11/site-packages/xesmf/smm.py:131: UserWarning: Input array is not C_CONTIGUOUS. Will affect performance.
  warnings.warn('Input array is not C_CONTIGUOUS. ' 'Will affect performance.')

../_images/notebooks_9e_PhaseCurves_29_1.png

../_images/notebooks_9e_PhaseCurves_29_2.png

Run cloudy phase curve!

[19]:

cldout = case_cld.phase_curve(opacity, n_cpu = 3,#jdi.cpu_count(),
                             full_output=True)

[20]:

#same old same old
wno =[];fpfs=[];legend=[]
for iphase in cldout.keys():
    w,f = jdi.mean_regrid(cldout[iphase]['wavenumber'],
                               cldout[iphase]['fpfs_thermal'],R=100)
    wno+=[w]
    fpfs+=[f*1e6]
    legend +=[str(int(iphase*180/np.pi))]
jpi.show(jpi.spectrum(wno, fpfs, plot_width=500,legend=legend,
                     palette=jpi.pals.viridis(n_phases)))

[21]:

to_plot = 'fpfs_thermal'#or thermal
collapse = [1.1,1.4,1.6]#micron
R=100
phases, all_curves, all_ws, fig=jpi.phase_curve(cldout, to_plot, collapse=collapse, R=100)
jpi.show(fig)

[22]:

fig=jpi.phase_snaps(cldout, x='longitude', y='pressure',z='taucld', palette='PuBu',
                y_log=True, z_log=False, collapse=['np.max','np.mean'])
#this will collapse the latitude axis by taking a mean, and the wavelength axis also by taking the mean

../_images/notebooks_9e_PhaseCurves_34_0.png

[23]:

fig, ax = jpi.plt.subplots(figsize=(6, 4))
case_cld.inputs['clouds']['profile']['opd'].isel(phase=0,
    lat=3,wno=0).plot(
    x='lon',y='pressure',ax=ax)
ax.set_ylim([3e2,1e-6])
ax.set_yscale('log')

../_images/notebooks_9e_PhaseCurves_35_0.png

[ ]:

Phase Curves Part 1¶

Run Thermal Phase Curve w/ 3D Input¶

Add Chemistry (randomized for purposes of the example)¶

Determining the `zero_point` and the `shift` parameter¶

Analyze Phase Curves¶

Phase Curve Plotting Tools: Phase Snaps¶

Phase Curve Plotting Tools: Phase Curves¶

Run Cloudy Thermal Phase Curve w/ 3D Input¶

Phase Curves Part 1¶

Run Thermal Phase Curve w/ 3D Input¶

Add Chemistry (randomized for purposes of the example)¶

Determining the zero_point and the shift parameter¶

Analyze Phase Curves¶

Phase Curve Plotting Tools: Phase Snaps¶

Phase Curve Plotting Tools: Phase Curves¶

Run Cloudy Thermal Phase Curve w/ 3D Input¶

Determining the `zero_point` and the `shift` parameter¶