# First steps with exo_k¶

Author: Jeremy Leconte (CNRS/LAB/Univ. Bordeaux)

The goal of exo_k is to provide a library to:

• Interpolate efficiently and easily in correlated-k and cross section tables.

• Convert easily correlated-k and cross section tables from one format to another (hdf5, LMDZ GCM, Exomol, Nemesis, PetitCode, TauREx, ExoREM, ARCIS, etc.).

• Adapt precomputed correlated-k tables to your needs by changing:

• the spectral and quadrature (g) grids,

• the pressure/temperature grid.

• Create tables for a mix of gases using tables for individual gases.

• Create your own tables from high-resolution spectra.

• Use your data in an integrated radiative transfer framework to simulate planetary atmospheres.

This tutorial will show you how to do all that.

If you are reading this online, the executable notebook can be found there: https://forge.oasu.u-bordeaux.fr/jleconte/exo_k-public/-/blob/public/tutorial-exo_k.ipynb

## Initialization¶

First, let’s make some common initializations.
We also import our library: exo_k.

For this to work, you need to have installed the library by either typing

pip install exo_k


anywhere or

pip install -e .


at the root of exo_k directory (the location where you first saw this tutorial).

[100]:

import exo_k as xk
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
import time,sys,os

[101]:

# Uncomment the line below if you want to enable interactive plots
#%matplotlib notebook
plt.rcParams["figure.figsize"] = (7,4)
from matplotlib import cycler
colors = cycler('color',[plt.cm.inferno(i) for i in np.linspace(0.1,1,5)])
plt.rc('axes', axisbelow=True, grid=True, labelcolor='dimgray', labelweight='bold', prop_cycle=colors)
plt.rc('grid', linestyle='solid')
plt.rc('xtick', direction='in', color='dimgray')
plt.rc('ytick', direction='in', color='dimgray')
plt.rc('lines', linewidth=1.5)


Important

Many global options can be changed after having imported the library using a syntax of the type xk.Settings().method_name(args). Many examples will be provided along the way, but you can have a look at exo_k.settings.Settings for all the possible options.

## Documentation¶

The documentation in html form can be found at http://perso.astrophy.u-bordeaux.fr/~jleconte/exo_k-doc/index.html

It is searchable, and has a general index to find any function.

To access it through a python interface, just run the following line for any class, function, etc.:

[8]:

help(xk.Ktable().bin_down)

Help on method bin_down in module exo_k.ktable:

bin_down(wnedges=None, weights=None, ggrid=None, remove_zeros=False, num=300, use_rebin=False, write=0) method of exo_k.ktable.Ktable instance
Method to bin down a kcoeff table to a new grid of wavenumbers (inplace).

Parameters
----------
wnedges: array
Edges of the new bins of wavenumbers (cm-1)
onto which the kcoeff should be binned down.
if you want Nwnew bin in the end, wnedges.size must be Nwnew+1
wnedges[0] should be greater than self.wnedges[0] (JL20 not sure anymore)
wnedges[-1] should be lower than self.wnedges[-1]
weights: array, optional
Desired weights for the resulting Ktable.
ggrid: array, optional
Desired g-points for the resulting Ktable.
Must be consistent with provided weights.
If not given, they are taken at the midpoints of the array
given by the cumulative sum of the weights



# Dealing with Ktable() objects¶

## Loading and saving a Ktable() object¶

For this tutorial to work, you should launch this notebook from a base_directory that contains a data/corrk/ directory where your correlated k files are storred. A directory with some sample files can be downloaded at the following url:

https://mycore.core-cloud.net/index.php/s/w2cHuigAiwcfBVW

Many thanks to K. Chubb and the Exomol project for providing these data.

### Loading a Ktable from a file¶

One of the main objects we will deal with is the Ktable() object. This contains a big matrix with k-coefficients for a species or a mix of species along with all the needed supporting information such as the Pressure, Temperature, wavenumber, and g grids onto which the k-coefficients have been computed.

To instantiate a Ktable() object from most currently supported formats, just run the following line. If you are specifying the file through the filename keyword, you have to give the full path to the file. The format should be recognized from the extension (see http://perso.astrophy.u-bordeaux.fr/~jleconte/exo_k-doc/units.html for supported formats).

exo_k.ktable.Ktable() for details on arguments and options.

[9]:

h2o_ktab=xk.Ktable(filename='data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5')
print(h2o_ktab)


file         : data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5
molecule     : H2O
p grid       : [1.00000000e-05 2.15443469e-05 4.64158883e-05 1.00000000e-04
2.15443469e-04 4.64158883e-04 1.00000000e-03 2.15443469e-03
4.64158883e-03 1.00000000e-02 2.15443469e-02 4.64158883e-02
1.00000000e-01 2.15443469e-01 4.64158883e-01 1.00000000e+00
2.15443469e+00 4.64158883e+00 1.00000000e+01 2.15443469e+01
4.64158883e+01 1.00000000e+02]
p unit       : bar
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : cm^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [  22   27 1538   20]



### Managing Units¶

#### Tracking units¶

Except for self-defining formats (such as hdf5), most formats handled by exo_k do not carry the information on units with them (Exo_transmit, HITRAN .cia, etc.).

However, previous experience has shown us that tracking units is essential to avoid errors and working with SI units is recommended when possible.

For this reason, the units for pressure, kdata, and wavenumbers are kept as attributes of any Ktable object. They can be seen by printing the object (as above). If the input format is self-defining (e.g. hdf5) the units are read in the file. For all the other formats, the units are assumed to always be the same and have been inferred from reading the codes using these formats and their documentation (see http://perso.astrophy.u-bordeaux.fr/~jleconte/exo_k-doc/units.html for some examples and a more general discussion).

#### Converting to new units¶

To control the unit we want to work with, we just need to specify what units we want with the p_unit and kdata_unit keywords.

These keywords accept string input with units recognized by the astropy.units library such as ‘Pa’, ‘mbar’, and ‘bar’ for pressure and e.g. ‘m^2/molecule’ and ‘cm^2/molecule’ for cross sections. Although ‘/molecule’ is not a unit recognized by astropy.units, it is automatically appended to the unit string provided by the user when necessary as a reminder that opacities per ‘moles’ or ‘kg’ are not supported.

To use SI units, just reload your ktable as shown below. Note that you can also force SI units to be used with the global option to be set once and for all.

xk.Settings().set_mks(True)

[10]:

h2o_ktab_SI=xk.Ktable(filename='data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5',
kdata_unit='m^2', p_unit='Pa')
# or simply
xk.Settings().set_mks(True)
h2o_ktab_SI=xk.Ktable(filename='data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5')
print(h2o_ktab_SI)


file         : data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5
molecule     : H2O
p grid       : [1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [  22   27 1538   20]



Now both the pressure unit and the pressure grid have changed!

#### Overriding default units in a file¶

If for some reason you know that the units used in a given input file are different from the default ones assumed for this format, you can always override the default units by explicitly stating what are the pressure (file_p_unit keyword) and k-coefficient (file_kdata_unit keyword) units.

Let’s say, for example, that you know the pressures in the above file where in fact given in mbar, you could specify it using

[11]:

xk.Ktable(filename='data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5', file_p_unit='mbar')

Be careful, you are assuming that p_unit is mbar
but the input file says that it is bar. The former will be used

[11]:


file         : data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5
molecule     : H2O
p grid       : [1.00000000e-03 2.15443469e-03 4.64158883e-03 1.00000000e-02
2.15443469e-02 4.64158883e-02 1.00000000e-01 2.15443469e-01
4.64158883e-01 1.00000000e+00 2.15443469e+00 4.64158883e+00
1.00000000e+01 2.15443469e+01 4.64158883e+01 1.00000000e+02
2.15443469e+02 4.64158883e+02 1.00000000e+03 2.15443469e+03
4.64158883e+03 1.00000000e+04]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [  22   27 1538   20]



Notice how the code pressure unit as changed, but not the actual values in the pressure grid. This is normal, you did not ask for a conversion, you just specified what units you were expecting from the file.

### Setting up the search path and searching files using regular expressions¶

At some point you might get tired of always typing the whole path to your files. To avoid that, we can tell exo_k where to look for files. By default, it searches the local directory where the library has been imported (“.”), as can be seen by running:

[12]:

xk.Settings().reset_search_path()
xk.Settings().search_path

[12]:

{'ktable': ['/Users/jleconte/atmosphere/RadiativeTransfer/exo_k'],


This also shows that exo_k keeps track of four different path types in this dictionary:

• 'ktable': for k-coefficient tables (i.e. Ktable objects)

• 'xtable': for cross-section tables (i.e. Xtable objects, see below)

• 'cia': for Collision Induced Absorption coefficients (i.e. Cia_table objects, see below)

• 'aerosol': for aerosol optical properties (i.e. Atable objects). Those are still experimental and will not be detailed here.

A new path can be added to the search path as follows

[102]:

xk.Settings().add_search_path('data/xsec', path_type='xtable')

[14]:

xk.Settings().search_path

[14]:

{'ktable': ['/Users/jleconte/atmosphere/RadiativeTransfer/exo_k',


If you do not want to add a path, but to reset the path (in order for exo_k not to see some files in the local directory for example), you can do it as follows:

[15]:

xk.Settings().set_search_path('data/corrk', path_type='ktable')
xk.Settings().search_path['ktable']

[15]:

['/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/corrk']


Notice that these methods can take as many directories as you want following one of the syntaxes below

[16]:

xk.Settings().set_search_path('.','data/corrk', path_type='ktable')
print(xk.Settings().search_path['ktable'])
xk.Settings().set_search_path(*['.','data/corrk'], path_type='ktable')
print(xk.Settings().search_path['ktable'])
xk.Settings().set_search_path(*('.','data/corrk'), path_type='ktable')
print(xk.Settings().search_path['ktable'])

['/Users/jleconte/atmosphere/RadiativeTransfer/exo_k', '/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/corrk']


Once we have done that, we can only give part of the name of the file, possibly in bits and pieces, as long as it is specific enouh to identify a unique file.

Warning

The pattern matching is done using regular expressions (python re module). Some special characters can have a special meaning. A dot, for example, stands for any single character. “.*” stands for any number of any characters. If you want to use these as usual characters, do not forget to put an escaping backslash in front of it.

[17]:

h2o_ktab_SI=xk.Ktable('H2O_R300_0.3-50mu.ktable.TauREx')
print('file found:',h2o_ktab_SI.filename)

h2o_ktab_SI=xk.Ktable('H2O_R300_0\.3\-50mu\.ktable\.TauREx')
print('file found:',h2o_ktab_SI.filename)

h2o_ktab_SI=xk.Ktable('H2O','R300.*tauREx')
print('file found:',h2o_ktab_SI.filename)

h2o_ktab_SI=xk.Ktable('H2','R300','tauREx')
print('file found:',h2o_ktab_SI.filename)

try:
h2o_ktab_SI=xk.Ktable('R300','tauREx',mol='H2')
except:
print("""There is no Ktable for H2.
Notice how the previous statement was looking for 'H2',
but settled for H2O because the string 'H2O' contains 'H2'.
If you want to avoid that, really specify a molecule with the mol keyword.""")


file found: /Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5
There is no Ktable for H2.
Notice how the previous statement was looking for 'H2',
but settled for H2O because the string 'H2O' contains 'H2'.
If you want to avoid that, really specify a molecule with the mol keyword.


Tip

At any moment, you can override the global path by specifying a local search path:  h2o_ktab_SI=xk.Ktable('H2','R300','tauREx', search_path='data/corrk/') 

This can be particularly usefull if you have different files with the same names in different directories.

By default, the filters you provide are compared to the file names in a case insensitive way to enable more flexibility. This can be changed with

xk.Settings().set_case_sensitive(True)


### Writing a Ktable into a file¶

To save any Ktable() object to a file, just use the relevant write_format('path_to_write') method. At the time of writing of this tutorial, the available writing methods are:

• write_hdf5(): General self defining format, which makes it the most optimal for exo_k (based on the Exomol Format). This method accepts the exomol_units=True option to force the output units to be the same as the native Exomol ones (i.e. directly compatible with the TauREX and petitRADTRANS codes).

• write_LMDZ(): Format used by the LMDZ Generic Global Climate Model (see LMDZ section below).

• write_nemesis(): Format used by the Nemesis retrieval code (binary; .kta).

• write_arcis(): Format used by the ARCIS forward model and retrieval code (fits; .fits).

• write_pickle(): General python/pickle format.

When needed, the data are converted to the right units before writing.

Notice that the full path (relative or absolute) must be provided, except for the file extension that will be added automatically.

[18]:

h2o_ktab_SI=xk.Ktable('H2O_R300_0.3-50mu.ktable.TauREx', kdata_unit='m^2', p_unit='Pa')
h2o_ktab_SI.write_hdf5('data/corrk/H2O_R300_0.3-50mu.ktable.SI')


This new file has saved the unit information. So now, it can be loaded with only:

[103]:

h2o_ktab_SI=xk.Ktable('H2O_R300_0.3-50mu.ktable.SI')

[20]:

h2o_ktab_SI

[20]:


molecule     : H2O
p grid       : [1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [  22   27 1538   20]



### What is inside a Ktable object¶

To see what is inside a Ktable() object, you can list the attributes using __dict__

Warning

Although you can access and see these attributes, they should NOT be changed manually as changing one of them usually implies that others need to be changed too. There is a method to change all the attributes that can possibly be changed.

The most useful attributes are the following:

• mol, Name of the mol described

• Ng, Number of g points

• filename, Name of the file used to load the Ktable

• ggrid, Abscissa of the g grid

• weights, quadrature weights corresponding to the g points

• kdata, Numpy array with the coeficients

• kdata_unit, Units of the kdata

• pgrid, pressure grid

• logpgrid, Log 10 of the pressure grid

• Np, Number of pressure points

• p_unit, Units of the pressure grid

• tgrid, Temperature grid

• Nt, Number of Temperature points

• wns, Central points of the wavenumber bins (in cm^-1)

• wnedges, Edges of the wavenumber bins (in cm^-1)

• wls, Wavelengths (10000./wns in microns)

• wledges, Wavelengths (10000./wnedges in microns)

• Nw, Number of wavenumber points

• shape, shape of the kdata array

[21]:

print(h2o_ktab_SI.mol)
print(h2o_ktab_SI.pgrid)

H2O
[1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]


the data themselves can also be accessed directly as follows (but they cannot be modified this way):

[22]:

print(h2o_ktab_SI.kdata[0,0,0,:])
print(h2o_ktab_SI[0,0,0,:])

[1.97818148e-30 1.99009429e-30 2.01152500e-30 2.04241642e-30
2.08267365e-30 2.13212071e-30 2.19045642e-30 2.25722151e-30
2.33171732e-30 2.41298132e-30 2.49969054e-30 2.59013679e-30
2.68222162e-30 2.77340390e-30 2.86082684e-30 2.94136768e-30
3.01184614e-30 3.06921780e-30 3.11082349e-30 3.13465184e-30]
[1.97818148e-30 1.99009429e-30 2.01152500e-30 2.04241642e-30
2.08267365e-30 2.13212071e-30 2.19045642e-30 2.25722151e-30
2.33171732e-30 2.41298132e-30 2.49969054e-30 2.59013679e-30
2.68222162e-30 2.77340390e-30 2.86082684e-30 2.94136768e-30
3.01184614e-30 3.06921780e-30 3.11082349e-30 3.13465184e-30]


### Handling zeros¶

At wavelengths where data are not available, the ktable can be filled with zeros that can become a problem for some codes (e.g. because of log interpolation).

If you want to avoid that, we have created the remove_zeros=True option that will replace zeros with a number 10 orders of magnitude smaller than the smallest number in the table.

h2o_ktab_SI=xk.Ktable(filename, remove_zeros=True, **kwargs)


You can also call manually the remove_zeros() method on a Ktable()object.

### Copying a Ktable object¶

If for some reason, as we will see below, you need to copy a Ktable object before modifying it inplace (to change resolution for example), you can use the copy method. This is similar to the deep copy in numpy.

[23]:

new_h2o_ktab=h2o_ktab_SI.copy()
new_h2o_ktab

[23]:


molecule     : H2O
p grid       : [1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [  22   27 1538   20]



Sometimes one needs to copy only the structure and metadata, but not the big self.kdata table. Then, you can use:

new_h2o_ktab=h2o_ktab_SI.copy(cp_kdata=False)


## Plotting opacities and k-distributions¶

plot_spectrum(self,ax,p=1.e-5,t=200.,g=0.,xscale=None,xscale=None,x=1.,**kwarg)


method is there to make it easy to plot at the requested (p,T,g) point. You can also normalize the cross section with the volume mixing ratio you want for the species with the x=vmr keyword.

[24]:

# %matplotlib notebook
p_plot=1.e5
t_plot=1000.
fig,axs=plt.subplots(2,2,sharex=False,sharey=False,figsize=(9,6.5))
h2o_ktab_SI.plot_spectrum(axs[0,0],p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='g=1')
h2o_ktab_SI.plot_spectrum(axs[0,0],p=p_plot,t=t_plot,g=0.,yscale='log',xscale='log',label='g=0')
axs[0, 0].set_title('Log-Log')
h2o_ktab_SI.plot_spectrum(axs[1,0],p=p_plot,t=t_plot,g=1.,xscale='log',label='g=1')
h2o_ktab_SI.plot_spectrum(axs[1,0],p=p_plot,t=t_plot,g=0.,xscale='log',label='g=0')
axs[1, 0].set_title('Lin-Log')
h2o_ktab_SI.plot_spectrum(axs[0,1],p=p_plot,t=t_plot,g=1.,yscale='log',label='g=1')
h2o_ktab_SI.plot_spectrum(axs[0,1],p=p_plot,t=t_plot,g=0.,yscale='log',label='g=0')
axs[0, 1].set_title('Log-Lin')
h2o_ktab_SI.plot_spectrum(axs[1,1],p=p_plot,t=t_plot,g=1.,label='g=1')
h2o_ktab_SI.plot_spectrum(axs[1,1],p=p_plot,t=t_plot,g=0.,label='g=0')
axs[1, 1].set_title('Lin-Lin')
for axs in axs:
for ax in axs:
ax.legend()
#ax.set_ylim(bottom=1.e-37)
fig.tight_layout()


The g gistribution is shown with plot_distrib. Notice that when you ask for a log scale in g, we switch to 1-g as it is the most relevant quantity to look at in log.

[25]:

p_plot=1.e5
t_plot=1000.
fig,axs=plt.subplots(2,2,sharex=False,sharey=False,figsize=(9,6.5))
h2o_ktab_SI.plot_distrib(axs[0,0],p=p_plot,t=t_plot,wl=1.0,yscale='log',xscale='log',label='$\lambda=1\mu$m')
h2o_ktab_SI.plot_distrib(axs[0,0],p=p_plot,t=t_plot,wl=10.,yscale='log',xscale='log',label='$\lambda=10\mu$m')
h2o_ktab_SI.plot_distrib(axs[1,0],p=p_plot,t=t_plot,wl=1.0,xscale='log',label='$\lambda=1\mu$m')
h2o_ktab_SI.plot_distrib(axs[1,0],p=p_plot,t=t_plot,wl=10.,xscale='log',label='$\lambda=10\mu$m')
h2o_ktab_SI.plot_distrib(axs[0,1],p=p_plot,t=t_plot,wl=1.0,yscale='log',label='$\lambda=1\mu$m')
h2o_ktab_SI.plot_distrib(axs[0,1],p=p_plot,t=t_plot,wl=10.,yscale='log',label='$\lambda=10\mu$m')
h2o_ktab_SI.plot_distrib(axs[1,1],p=p_plot,t=t_plot,wl=1.0,label='$\lambda=1\mu$m')
h2o_ktab_SI.plot_distrib(axs[1,1],p=p_plot,t=t_plot,wl=10.,label='$\lambda=10\mu$m')
for axs in axs:
for ax in axs:
ax.legend()
#ax.set_ylim(bottom=1.e-37)
fig.tight_layout()


## Interpolating and remapping tables¶

### Basic interpolation¶

To get the value of the cross section of the molecule or mix described by our Ktable() object, we can use:

self.interpolate_kdata(logp_array=None,t_array=None,log_interp=None,wngrid_limit=None)


This returns the cross section of the molecule along a LogP-T profile specified by logp_array and t_array. This returns an array of shape (logp_array.size, self.Nw, self.Ng), i.e (number of PT points, number of wavenumber points, number of g-points). logp_arrayand t_arraymust have the same size.

It also works if logp_array and/or t_array are floats

Important

Notice that the input is always in Log10(Pressure) because this is the quantity used for the interpolation

exo_k.data_table.Data_table.interpolate_kdata() for details on arguments and options.

[26]:

logp=3.
Temp=1000.
tmp_opa=h2o_ktab_SI.interpolate_kdata(logp,Temp)
print(tmp_opa.shape)
print(tmp_opa)

(1, 1538, 20)
[[[1.57577243e-26 1.70593813e-26 1.71642784e-26 ... 4.95008363e-22
8.13668554e-22 9.35286142e-22]
[2.90167053e-26 4.18911348e-26 4.71088691e-26 ... 8.99173698e-23
4.25319413e-22 6.23065168e-22]
[5.81193968e-27 5.82883586e-27 5.92632233e-27 ... 1.21096603e-25
3.52299381e-25 4.86654868e-25]
...
[2.43561538e-34 3.12899011e-34 3.55619421e-34 ... 2.18003828e-33
2.92876500e-33 3.94600721e-33]
[2.71942214e-34 3.29638704e-34 3.74713705e-34 ... 2.21029744e-33
3.13286094e-33 5.08112258e-33]
[2.86054904e-34 3.25393178e-34 3.70676091e-34 ... 2.45139953e-33
3.46404319e-33 7.47015700e-33]]]

[27]:

logp_array=np.linspace(1.,6.,6)
t_array=[1000. for logp in logp_array]
tmp_opa=h2o_ktab_SI.interpolate_kdata(logp_array,t_array)
print(tmp_opa.shape)

(6, 1538, 20)


### Interpolation options¶

The default interpolation scheme interpolates log(kdata). This can be changed directly using the log_interp=False keyword. But the default behavior can also be changed using Settings().set_log_interp. The keyword always supercedes the default behavior.

[28]:

logp=3.
Temp=1000.
tmp_opa_lin1=h2o_ktab_SI.interpolate_kdata(logp,Temp,log_interp=False)

[29]:

xk.Settings().set_log_interp(False)
logp=3.
Temp=1000.
tmp_opa_lin2=h2o_ktab_SI.interpolate_kdata(logp,Temp)
np.sum(tmp_opa_lin2-tmp_opa_lin1)
xk.Settings().set_log_interp(True) ## let's go back to default behavior


Linear interpolation is of course faster, but is believed to be less accurate.

[30]:

logp=3.
Temp=1000.
%timeit -n 100 -r 2 h2o_ktab_SI.interpolate_kdata(logp,Temp,log_interp=True)
%timeit -n 100 -r 2 h2o_ktab_SI.interpolate_kdata(logp,Temp,log_interp=False)

1.92 ms ± 13.5 µs per loop (mean ± std. dev. of 2 runs, 100 loops each)
289 µs ± 43.3 µs per loop (mean ± std. dev. of 2 runs, 100 loops each)

If we just want a part of the spectrum, we can specify a wavenumber range (in cm^-1) to cover as follows.
You can see that it goes much faster.
[31]:

logp=3.
Temp=1000.
wn_range=[6000.,5000.] # does not need to be sorted
%timeit -n 100 -r 2 h2o_ktab_SI.interpolate_kdata(logp,Temp,wngrid_limit=wn_range,log_interp=True)
%timeit -n 100 -r 2 h2o_ktab_SI.interpolate_kdata(logp,Temp,wngrid_limit=wn_range,log_interp=False)
tmp_opa=h2o_ktab_SI.interpolate_kdata(logp,Temp,wngrid_limit=wn_range)
print(tmp_opa.shape)
print(tmp_opa)

240 µs ± 23 µs per loop (mean ± std. dev. of 2 runs, 100 loops each)
222 µs ± 40.2 µs per loop (mean ± std. dev. of 2 runs, 100 loops each)
(1, 54, 20)
[[[4.07753036e-29 6.85197398e-29 1.02732246e-28 ... 2.41869028e-25
4.85016488e-25 1.14633633e-24]
[5.57875921e-29 7.97288442e-29 1.21104087e-28 ... 3.42969506e-25
1.22902497e-24 2.53733317e-24]
[6.17042832e-29 9.10883868e-29 1.31081205e-28 ... 3.53333315e-25
8.73232110e-25 3.80127738e-24]
...
[2.02444510e-30 3.93145431e-30 5.47901786e-30 ... 2.32477558e-27
7.86972907e-27 3.78365657e-26]
[2.33493343e-30 3.54354098e-30 4.84518937e-30 ... 2.32501304e-27
6.12582068e-27 2.35847913e-26]
[2.50675050e-30 3.68390571e-30 4.92245284e-30 ... 2.94029197e-27
1.10535798e-26 4.63705729e-26]]]


### (Log P, T) grid Remapping¶

If for computational efficiency reason we just want to use a subset of the (Log P, T) grid, we can remap the table. This will modify the Table in place. So we might want to make a copy before.

[32]:

Small_Grid_h2o_ktab=h2o_ktab_SI.copy()
logp_array=np.arange(6.)
t_array=np.linspace(200.,1000.,9)

Small_Grid_h2o_ktab.remap_logPT(logp_array=logp_array,t_array=t_array)
Small_Grid_h2o_ktab

[32]:


molecule     : H2O
p grid       : [1.e+00 1.e+01 1.e+02 1.e+03 1.e+04 1.e+05]
p unit       : Pa
t grid   (K) : [ 200.  300.  400.  500.  600.  700.  800.  900. 1000.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [   6    9 1538   20]



### Remapping in g-space¶

Sometimes, we also want to reduce the number of g-points, or combine two tables with different g-quadratures. For this we can use the following:

[33]:

# simple function to compute a set of gauss legendre points and weights
from exo_k.util.interp import gauss_legendre
weights, ggrid, gedges=gauss_legendre(8)
ggrid

[33]:

array([0.01985507, 0.10166676, 0.2372338 , 0.40828268, 0.59171732,
0.7627662 , 0.89833324, 0.98014493])

[34]:

h2o_ktab_8g=h2o_ktab_SI.copy()
h2o_ktab_8g.remap_g(ggrid=ggrid, weights=weights)
h2o_ktab_8g

[34]:


molecule     : H2O
p grid       : [1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.05061427 0.11119052 0.15685332 0.18134189 0.18134189 0.15685332
0.11119052 0.05061427]
data oredered following p, t, wn, g
shape        : [  22   27 1538    8]



The table now has only 8 g points. You can of course specify custom g-points and weights if you want.

### Extending the spectral range¶

As we will see afterwards, we need tables to have the same spectral grid to be able to combine them. If they cover the same spectral range, but with different resolution, the solution is to bin them down as explained a few sections down. If, however, we want to combine molecules that are absorbing in rather different part of the spectrum, we may often encounter tables that cover only the spectral region where each molecule has a significant absorption.

To circumvent this issue, exo_k allows you to extend the spectral range of your table adding a null opacity (or a very small one if one wants to use log interpolation afterwards) in the new spectral bins. This can be done as follows. It changes the table in place.

There are 3 options to specify the left and/or right extension of the spectral grid:

• The wngrid_left/right keyword allows you to specify the center of the new bins. The edges are computed automatically. This option is most useful with Xtables for which you most want to control the center of the bins.

• The wnedges_left/right keyword allows you to specify the the edges of the new spectral bin. Bin centers are computed automatically. This option is most useful with Ktables for which you most want to control the edges of the bins. Since we are re-using the last wnedges on the left and right of the current grid, the user only needs to specify N new values to create Nnew bins.

• You can specify both at the same time.

Warning

The user must make sure that the new left and right wavenumber grids do not overlap with the current one. If both wnedges and wngrids are specified, the user should also make sure that there is one (and only one) wavenumber point between each pair of wnegdes points.

[35]:

large_spectral_range_h2o=h2o_ktab_SI.copy()
new_wn_points_left=np.linspace(10,190,20)
new_wn_edges_left=np.linspace(5,185,20)
new_wn_points_right=np.linspace(33400,50000,20)
new_wn_edges_right=np.linspace(33450,50050,20)

# specifying only wngrid
#large_spectral_range_h2o.extend_spectral_range(wngrid_left=new_wn_points_left, wngrid_right=new_wn_points_right, remove_zeros=True)
large_spectral_range_h2o.extend_spectral_range(wngrid_left=new_wn_points_left, wnedges_left=new_wn_edges_left,
wngrid_right=new_wn_points_right, wnedges_right=new_wn_edges_right, remove_zeros=True)
#large_spectral_range_h2o.extend_spectral_range( wnedges_left=new_wn_edges_left,wnedges_right=new_wn_edges_right, remove_zeros=True)
print('Old wn grid: ',h2o_ktab_SI.wns)
print('New wn grid: ',large_spectral_range_h2o.wns)
print('Old wn edges: ',h2o_ktab_SI.wnedges)
print('New wn edges: ',large_spectral_range_h2o.wnedges)

p_plot=1.e5
t_plot=1000.
fig,ax=plt.subplots(1,1,sharex=False,sharey=False,figsize=(6,4.5))
large_spectral_range_h2o.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',color='black',linestyle='dashed',label='after extension')
h2o_ktab_SI.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,label='before extension')
ax.legend()
fig.tight_layout()

Old wn grid:  [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
New wn grid:  [1.00000000e+01 1.94736842e+01 2.89473684e+01 ... 4.82526316e+04
4.91263158e+04 5.00000000e+04]
Old wn edges:  [  199.57138937   200.23662733   200.90408276 ... 33112.21730444
33222.59136213 33333.33333333]
New wn edges:  [5.00000000e+00 1.44736842e+01 2.39473684e+01 ... 4.83026316e+04
4.91763158e+04 5.00500000e+04]


### Clipping data to a given spectral range¶

If for effciency reasons you want to only deal with a small part of the spectral domain, you can use the Ktable.clip_spectral_range(wl_range= , wn_range= ) method with either a wavelength range (wl_rangein microns) or a wavenumber range (wn_range in cm^-1). Data outside this range will be removed from the table.

## Binning down a Ktable¶

We now get to one of the core features of exo_k: we want to create a new Ktable() object with just the spectral resolution we need in the model we want to run (whether it is a Global Climate Model, a 1D climate model, or a simple forward model for a retrieval).

You first need to provide a grid, wnedges, containing the edges of the new wavenumber bins. So the final Ktable will have a spectral dimension with size wnedges.size-1.

You can of course provide any custom grid you want. But if you want a grid at fixed resolution you can use the following

xk.wavenumber_grid_R(min_wavenumber, max_wavenumber, R)


where $$R=\frac{\nu}{\delta \nu}$$.

[36]:

R=10
wavenumber_grid=xk.wavenumber_grid_R(200.,10000.,R)
wavenumber_grid

[36]:

array([  200.        ,   221.03418362,   244.28055163,   269.97176152,
298.36493953,   329.74425414,   364.42376008,   402.75054149,
445.1081857 ,   491.92062223,   543.65636569,   600.83320479,
664.02338455,   733.85933352,   811.03999337,   896.33781407,
990.60648488,  1094.78947835,  1209.92949288,  1337.17888846,
1477.81121979,  1633.23398251,  1805.00269989,  1994.83649096,
2204.63527613,  2436.49879214,  2692.747607  ,  2975.94634497,
3288.92935422,  3634.82907389,  4017.10738464,  4439.59025629,
4906.50603942,  5422.52778413,  5992.82000948,  6623.09039174,
7319.64688874,  8089.46087201,  8940.23689866,  9880.48982111,
10000.        ])


Then we just need to use the bin_down(wnedges) method. However, bin_down modifies the Ktable() in place. So again, if we want to keep the original one in memory, we just copy it or use bin_down_cp.

Optional keywords can be weights= and ggrid= where you can provide the weights you want (and the corresponding g points). In general, for an LMDZ user, just input the values inside g.dat (except for the final zero) in a w_array, and use weights=w_array. If you want gauss legendre weights, you can have a look at exo_k.gauss_legendre in the index.

exo_k.ktable.Ktable.bin_down() bin_down_cp() for details on arguments and options.

It should only take a few seconds. (However, the bigger the final resolution wanted, the longer it takes).

[37]:

start=time.time()
LowRes_h2o_ktab=h2o_ktab_SI.bin_down_cp(wavenumber_grid)
end=time.time()
print('computation time=',end-start,'s')
LowRes_h2o_ktab

computation time= 16.91690421104431 s

[37]:


molecule     : H2O
p grid       : [1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [ 210.51709181  232.65736762  257.12615657  284.16835052  314.05459683
347.08400711  383.58715079  423.9293636   468.51440396  517.78849396
572.24478524  632.42829467  698.94135904  772.44966345  853.68890372
943.47214947 1042.69798161 1152.35948561 1273.55419067 1407.49505412
1555.52260115 1719.1183412  1899.91959542 2099.73588355 2320.56703413
2564.62319957 2834.34697599 3132.4378496  3461.87921405 3825.96822926
4228.34882046 4673.04814786 5164.51691178 5707.67389681 6307.95520061
6971.36864024 7704.55388037 8514.84888534 9410.36335988 9940.24491055]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [22 27 40 20]


[38]:

p_plot=1.e0
t_plot=300.
fig,axs=plt.subplots(1,2,sharex=False,sharey=False)
h2o_ktab_SI.plot_spectrum(axs[0],p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='Original data')
LowRes_h2o_ktab.plot_spectrum(axs[0],p=p_plot,t=t_plot,g=1.,label='Binned data')
h2o_ktab_SI.plot_spectrum(axs[1],p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='Original data')
LowRes_h2o_ktab.plot_spectrum(axs[1],p=p_plot,t=t_plot,g=1.,label='Binned data')
axs[1].set_xlim(left=1.,right=6.)
axs[1].set_xscale('linear')
for ax in axs:
ax.legend()
fig.tight_layout()


## Rescaling and mixing Ktables¶

The Ktable we have seen up to now contain opacity data (in terms of cross section in surface per molecule) for pure species. To create a mix, you can dilute your tables with (multiply by) the Volume mixing ratio of the relevant species and add them.

Important

Although we are using the ‘+’ sign, the tables are actually combined using the Random overlap method.

Warning

Ktables should have the same dimensions and grids to be combined. For efficiency, the code only checks whether the dimensions are the same (and raises a TypeError if they are not). The user should make sure the grids are the same as well.

[39]:

h2o=xk.Ktable('H2O_R300_0.3-50mu.ktable.SI', remove_zeros=True)
co2=xk.Ktable('CO2_R300_0.3-50mu.ktable.SI', remove_zeros=True)
for ktab in [co2,h2o]:
ktab.remap_logPT(logp_array=np.arange(0,5),t_array=np.arange(200.,500.,50.))
start=time.time()
mix=0.01*h2o+co2*0.02
print('Mixing time=',time.time()-start,'s')

Mixing time= 1.5422921180725098 s

[40]:

p_plot=1.e0
t_plot=300.
fig,ax=plt.subplots(1,1,sharex=False,sharey=False)
h2o.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,x=0.01,yscale='log',xscale='log',label='H2O')
co2.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,x=0.01,label='CO2')
mix.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='mix')
ax.set_ylim(bottom=1.e-36)
ax.legend()
fig.tight_layout()


Important

For the moment, the Ktable for the mix inherits attributes (such as the name of the molecule) from the most leftward Ktable in the operation.

The renormalization operation also accepts an array of volume mixing ratio, but only if it has the dimensions (Number of P points, number of T points).

[41]:

x_h2o=np.array([0.01+0.02*logP+0.00001*T for logP in h2o.logpgrid for T in h2o.tgrid]).reshape((h2o.Np,h2o.Nt))
mix=x_h2o*h2o+co2*0.02


# Dealing with mixes: The Kdatabase() object¶

## Loading a Kdatabase() object¶

A Kdatabase() object mainly contains a dictionary of Ktable() objects identified by the name of their molecule (in the attribute self.ktables, but they can be accessed with self['mol_name']). In addition, most of the attributes of the Ktable() are reloaded as attributes of Kdatabase().

To load a database, just specify a list of molecules and a string filter that is common to all your filenames.

For details on the available arguments and options, check-out exo_k.kdatabase.Kdatabase().

[42]:

molecs=['H2O','CO2']
filt='R300_0.3-50mu.ktable.SI'
database=xk.Kdatabase(molecs, filt)
database

[42]:

The available molecules are:
All tables share a common spectral grid
All tables share a common logP-T grid
All tables share a common p unit
All tables share a common kdata unit


It will look for files starting with the molecule names, followed by _, ., or -, and containing the filters.

If you want the moleculer name to be followed by something else than these characters, then just change the setting by using the xk.Settings().set_delimiters method.

Of course, this is unpractical when you do not have a homogeneous set of Ktables to start with. Then you can just feed a dictionary of molecules with the associated files path.

[43]:

database=xk.Kdatabase({'H2O':'data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5',
'CO2':'data/corrk/CO2_R300_0.3-50mu.ktable.SI.h5'},
p_unit='Pa', kdata_unit='m^2/molecule')
database

[43]:

The available molecules are:
H2O->data/corrk/H2O_R300_0.3-50mu.ktable.TauREx.h5
CO2->data/corrk/CO2_R300_0.3-50mu.ktable.SI.h5
All tables share a common spectral grid
All tables share a common logP-T grid
All tables share a common p unit
All tables share a common kdata unit


Warning

The latest release of Exomol data is using a new convention to write molecule names in filenames to include isotopologues. For example, the name of a file describing the water molecule now starts with ‘1H2-16O’ instead of ‘H2O’. To use these, use the new name in the method call (database=xk.Kdatabase([‘1H2-16O’,’12C-16O2’], filt)). However, be aware that inside the database object, individual Ktables are recognized and called by the name of the molecule that is stored inside the file (‘H2O’ and ‘CO2’ in our case; see the accessing individual Ktables section below).

For FITS and Nemesis formats that do not yet include the molecule name inside the file, we recommend that you use the second loading method: database=xk.Kdatabase({‘H2O’: ‘data/corrk/1H2-16O__POKAZATEL__R1000_0.3-50mu.ktable.ARCiS.fits’, etc.}).

### Adding individual, already-loaded Ktables¶

Sometimes, you want to modify/remap/bin_down your Ktables before loading it into a Kdatabase. However, you do not necessarily want to have to store these intermediate steps into a file to read it into a database as described above.

To do this, you can add one or several Ktable using the add_ktables(*ktables) method. If you want to start from scratch, you’ll need to create an empty Kdatabase as follows.

[44]:

database=xk.Kdatabase(None) # create empty database
print(database)

h2o=xk.Ktable('H2O_R300_0.3-50mu.ktable.SI', remove_zeros=True)
co2=xk.Ktable('CO2_R300_0.3-50mu.ktable.SI', remove_zeros=True)
co2.remap_logPT(logp_array=np.arange(0,5),t_array=np.arange(200.,500.,50.))

print(database)

The available molecules are:
All tables share a common spectral grid
All tables share a common logP-T grid
All tables share a common p unit
All tables share a common kdata unit

Carefull, not all tables have the same PT grid.
You'll need to use remap_logPT
The available molecules are:
All tables share a common spectral grid
All tables do NOT have common logP-T grid
You will need to run remap_logPT to perform some operations
All tables share a common p unit
All tables share a common kdata unit



### Dealing with inhomogeneous databases¶

As we are loading Ktables from various possible sources, that we may have modified however we like, we are bound to end up with a inhomogeneous database because all the Ktable:

• do not have the same datatype.

• do not have the same weights and g grid,

• do not have the same logP-T grid,

• do not have the same spectral grid,

• do not have the same units

For the moment, a Kdatabase needs to be homogeneous to be used to create mixes or compute forward models (see below).

The 1st issue is due to the fact that exo_k can handle other types of data (e.g. Xtable that are regular cross section tables; see dedicated section below). Although these data types can be handled in a database, a Kdatabase can only deal with a single type at a time. Loading different types will result in an error.

The 2nd issue is delicate because there are no easy way to accurately re-interpolate a Ktable on a different g grid. So for the moment, there is no way to remedy this. The user needs to provide Ktable with the same g grid and weights. Loading tables with different g grids will result in an error.

The remaining issues are easily handled by following the suggestion in the warning message (see above):

• Using Kdatabase.remap_logPT(logp_array= , t_array= ) will homogenize the log P-T grid by interpolating (in place) all the Ktables to a common grid.

• Using Kdatabase.bin_down(wnedges= , **kwargs) will homogenize the spectral grid by binning down (in place) all the Ktables to a common grid.

• If you are using cross sections, the spectral grid can be homogenized using Kdatabase.sample(wngrid= ) instead of bin_down (see below).

• All the Ktables can be converted to the same units using Kdatabase.convert_p_unit() and Kdatabase.convert_kdata_unit(), or directly Kdatabase.convert_to_mks()

[45]:

print(database)
print('remapping...')
database.remap_logPT(logp_array=np.arange(0,5),t_array=np.arange(200.,500.,50.))
print(database)

The available molecules are:
All tables share a common spectral grid
All tables do NOT have common logP-T grid
You will need to run remap_logPT to perform some operations
All tables share a common p unit
All tables share a common kdata unit

remapping...
The available molecules are:
All tables share a common spectral grid
All tables share a common logP-T grid
All tables share a common p unit
All tables share a common kdata unit



Warning

A Ktable is not deep copied when loaded into a database! Modifying the database inplace thus modifies the original table as well.

[46]:

# See how the values in the P and T grids have changed
# compared to their values in the previous section.
h2o

[46]:


molecule     : H2O
p grid       : [1.e+00 1.e+01 1.e+02 1.e+03 1.e+04]
p unit       : Pa
t grid   (K) : [200. 250. 300. 350. 400. 450.]
wn grid      : [  199.90345491   200.56979976   201.23836576 ... 33057.12210094
33167.31250795 33277.87021631]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [   5    6 1538   20]



## Accessing individual Ktables¶

A Kdatabase is essentially a dictionary and individual Ktable objects can be accessed using the molecule name as a key. When adding Ktables one by one, the molecule name used as a key is inherited from the Ktable.mol attribute. If you are loading a new Ktable with a molecule name that is already in the database, the latter will be replaced.

[47]:

database['H2O'].filename

[47]:

'/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/corrk/H2O_R300_0.3-50mu.ktable.SI.h5'


Warning

You cannot load a new table this way!

database['SiO']=SiO_ktable
TypeError: 'Kdatabase' object does not support item assignment


## Create a Ktable() for a mix of gases¶

One of the main interests of databases is to actually be able to create new correlated-k tables for mixtures of gases (the other being running a full atmospheric model as we will see at the end).

Just create a dictionary with the names of the molecules and the corresponding volume mixing ratios. Notice that the sum of the mixing ratios does not need to be equal to one. In that case, it is assumed that the rest of the gas is composed of radiatively inactive species. Then, the opacity given is in m^2 per total number of molecules (couting the inactive ones).

If a molecule is not in the database, it will just be considered radiatively inactive.

[48]:

molecs=['H2O','CO2']
suffix='R300_0.3-50mu.ktable.SI'
database=xk.Kdatabase(molecs,suffix)

composition={'H2O':0.001,'CO2':400.e-6,'N2':'background'}

start=time.time()
mixH2O_CO2=database.create_mix_ktable(composition)

mixH2O_CO2.write_hdf5('data/corrk/mix_H2O_CO2')

end=time.time()
print("computation time ",end - start,'s')

Do not have all the molecules in my database
Molecules to be treated:  {'CO2', 'H2O'}
treating molecule  H2O
computation time  32.64994812011719 s


Here, the 'N2':'background'in the composition dictionary tells the code that N2 is the background ‘filler’ gas. So whatever the vmr given before, N2 will be adjusted to get a total vmr of 1. Because N2 is not in the database (because it does not have any strong lines), its opacity will however not be considered.

(Why doing this you might ask? Well, here including N2 does not change anything as it is considered transparent anyway. This will become important however when we will deal with atmospheric models where including N2 will allow us to have the right mean molecular weight for the atmosphere).

[49]:

ptest=1.e4
ttest=300.
fig,axs=plt.subplots(1,2,sharex=True,sharey=True)
for species,x in composition.items():
if x is not 'background':
database[species].plot_spectrum(axs[0],p=ptest,t=ttest,g=1.,x=x,linestyle='dashed',yscale='log')
database[species].plot_spectrum(axs[1],p=ptest,t=ttest,g=0.,x=x,linestyle='dashed',yscale='log')
mixH2O_CO2.plot_spectrum(axs[0],p=ptest,t=ttest,g=1.)
mixH2O_CO2.plot_spectrum(axs[1],p=ptest,t=ttest,g=0.)
axs[0].set_title('g=1')
axs[1].set_title('g=0')
for ax in axs :
ax.set_xscale('log')
ax.label_outer()
ax.set_ylim(bottom=1.e-40)
fig.tight_layout()


## Binning down a Kdatabase()¶

Just like Ktable() objects, you can change the spectral resolution of an entire database with the same bin_down(wnedges) method. It will be applied to all ktables in place (weights keyword available as well).

[50]:

R=10
wavenumber_grid=xk.wavenumber_grid_R(200.,10000.,R)
database.bin_down(wavenumber_grid)


Notice how the resolution has changed below when we run the exact same code.

[51]:

ptest=1.e4
ttest=300.

fig,axs=plt.subplots(1,2,sharex=True,sharey=True)
for species,x in composition.items():
if x is not 'background':
database[species].plot_spectrum(axs[0],p=ptest,t=ttest,g=1.,x=x,linestyle='dashed',label=species)
database[species].plot_spectrum(axs[1],p=ptest,t=ttest,g=0.,x=x,linestyle='dashed',label=species)
mixH2O_CO2.plot_spectrum(axs[0],p=ptest,t=ttest,g=1.,label='mix')
mixH2O_CO2.plot_spectrum(axs[1],p=ptest,t=ttest,g=0.,label='mix')
axs[0].set_title('g=1')
axs[1].set_title('g=0')
for ax in axs :
ax.set_xscale('log')
ax.set_yscale('log')
ax.label_outer()
ax.set_ylim(bottom=1.e-40)
ax.legend()
fig.tight_layout()


# Loading and creating correlated-k tables to use with LMDZ GCM¶

## Creating new corr-k for LMDZ-GCM¶

exo_k has also been designed to create corr-k tables to be used with the LMDZ gcm.

To produce a LMDZ usable directory with the corr-k and all the supporting files, nothing is simpler.

Important

This section only deals with corr-k tables without variable gases. For these, we need to deal with gas mixes and Ktable5D that are discussed a little further down.

First, let’s create a grid for our IR and VI channels and create our two ktables.

[52]:

R=10
IR_grid=xk.wavenumber_grid_R(220.,10000.,R)
VI_grid=xk.wavenumber_grid_R(5200.,33000.,R)

IR_h2o_ktab=h2o_ktab_SI.bin_down_cp(IR_grid)

VI_h2o_ktab=h2o_ktab_SI.bin_down_cp(VI_grid)


Once your low resolution corrk table has been produced, just run

[53]:

IR_h2o_ktab.write_LMDZ('data/corrk/lmdz/new_h2o_corrk',band='IR')
VI_h2o_ktab.write_LMDZ('data/corrk/lmdz/new_h2o_corrk',band='VI')

Directory was already there data/corrk/lmdz/new_h2o_corrk
Directory was already there: data/corrk/lmdz/new_h2o_corrk/IR39
Directory was already there data/corrk/lmdz/new_h2o_corrk
Directory was already there: data/corrk/lmdz/new_h2o_corrk/VI19


Should you care about units?

Actually… no! Because we have tracked units from the beginning, and we know LMDZ likes ktables in cm^2/molec and pressures in log mBar, the conversion has been done automatically for you.

The LMDZ aficionados may have noticed that in the corrk/new_h2o_corrk directory, two subdirectories have been created with the names IR self.Nw and VI self.Nw that contain the corr-k tables whereas we would like a single directory called Number_of_IR_bins x Number_of_VI_bins. Then just run:

[54]:

xk.finalize_LMDZ_dir('data/corrk/lmdz/new_h2o_corrk',IR_h2o_ktab.Nw,VI_h2o_ktab.Nw)

Everything went ok. Your ktable is in: data/corrk/lmdz/new_h2o_corrk/39x19
You'll probably need to add Q.dat before using it though!


You just have to create a Q.dat file in the directory and you are good to run LMDZ with new corrk. If you do not know what I am talking about, you have not been using the LMDZ GCM long enough!

To load a LMDZ corr-k table, the initialization method needs some more information:

• The directory to use (path=)

• The resolution (res=)

• Whether you want to load the ‘IR’ or ‘VI’ band (band=)

The mol keyword can be given to specify which species is described. This will be useful for mixes.

Here, if you do not specify units, you will keep the LMDZ ones (cm^2 and mBar). But we advise to use the SI units. In particular because our radiative transfer method use SI units.

[55]:

IR_h2o_ktab=xk.Ktable(path='data/corrk/lmdz/new_h2o_corrk',res='39x19',band='IR',mol='H2O',
kdata_unit='m^2',p_unit='Pa')
print(IR_h2o_ktab)


file         : data/corrk/lmdz/new_h2o_corrk
molecule     : H2O
p grid       : [1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [ 231.56880099  255.92310439  282.83877223  312.58518557  345.46005652
381.79240782  421.94586586  466.32229996  515.36584436  569.56734336
629.46926376  695.67112414  768.83549494  849.69462979  939.05779409
1037.81936442 1146.96777977 1267.59543418 1400.90960974 1548.24455953
1711.07486126 1891.03017532 2089.91155497 2309.7094719  2552.62373755
2821.08551953 3117.78167359 3445.68163456 3808.06713546 4208.56505219
4651.18370251 5140.35296264 5680.96860295 6278.44128649 6938.75072067
7668.50550426 8475.00926841 9366.33377387 9917.13029426]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
data oredered following p, t, wn, g
shape        : [22 27 39 20]



## Creating a LMDZ corr-k with a variable gas from a Kdatabase: the Ktable5d¶

Here, we will finally see one of the strengths of Kdatabase objects. Because we can create mix of gases and bin down, we have everything we need to create LMDZ like ktables for variable gases.

For this, use

Kdatabase.create_mix_ktable5d(vgas_comp=, bg_comp=, x_array=)


where:

• vgas_comp must contain a dictionary of the composition of the variable gas (it can be a single molecule).

• bg_comp does the same for the background gas.

• x_array is an array of the vol mix ratios of the variable gas that we want to compute in the table.

Warning

A correlated-k table with a variable gas has one more dimension than a standard Ktable (for the mixing ratio of the variable species; attribute self.xgrid).

Therefore, such tables are handled with a special class: Ktable5d.

Most methods available to Ktable can be used with Ktable5d (like plotting, writing, interpolating, binning down, remapping, etc.). Note that for some of these methods, a new argument (x_array) is sometimes needed.

Because they already represent a mix and never a pure species, they cannot be used exactly the same:

• They cannot be normalized by a table of volume mixing ratios.

• They cannot be combined with other Ktable5d.

• They cannot be used in a Kdatabase along with other regular Ktables

[56]:

start=time.time()

molecs=['H2O','CO2']
suffix='R300_0.3-50mu.ktable.SI'

database=xk.Kdatabase(molecs,suffix)
# we first remap on a smaller logP T grid for speed
database.remap_logPT(logp_array=np.arange(6.), t_array=[200.,300.,400.,500.,600.])
x_array=[1.e-7,0.1]
vgas={'H2O':'background'}
bg_gas={'CO2':4.e-4,'N2':'background'}
mix_var_gas=database.create_mix_ktable5d(vgas_comp=vgas,bg_comp=bg_gas,x_array=x_array)
end=time.time()
print("computation time ",end - start)

Do not have all the molecules in my database
Molecules to be treated:  {'CO2'}
only 1 molecule: CO2
I have all the requested molecules in my database
{'H2O'}
only 1 molecule: H2O
shape of the output Ktable5d (p,t,x,wn,g): [   6    5    2 1538   20]
computation time  6.644855976104736


Then we can bin down before writing to file

[57]:

R=10
IR_grid=xk.wavenumber_grid_R(220.,10000.,R)
VI_grid=xk.wavenumber_grid_R(5200.,33000.,R)

IR_mix_var_gas_ktab=mix_var_gas.copy()
IR_mix_var_gas_ktab.bin_down(IR_grid)

VI_mix_var_gas_ktab=mix_var_gas.copy()
VI_mix_var_gas_ktab.bin_down(VI_grid)

IR_mix_var_gas_ktab.write_LMDZ('data/corrk/lmdz/new_earth_h2o-var',band='IR')
VI_mix_var_gas_ktab.write_LMDZ('data/corrk/lmdz/new_earth_h2o-var',band='VI')

xk.finalize_LMDZ_dir('data/corrk/lmdz/new_earth_h2o-var',IR_mix_var_gas_ktab.Nw,VI_mix_var_gas_ktab.Nw)

Directory was already there: data/corrk/lmdz/new_earth_h2o-var
Directory was already there: data/corrk/lmdz/new_earth_h2o-var/IR39
Directory was already there: data/corrk/lmdz/new_earth_h2o-var
Directory was already there: data/corrk/lmdz/new_earth_h2o-var/VI19
Everything went ok. Your ktable is in: data/corrk/lmdz/new_earth_h2o-var/39x19
You'll probably need to add Q.dat before using it though!

[58]:

IR_mix_var_gas_ktab.write_hdf5('data/corrk/lmdz/new_earth_h2o.h5', kdata_unit='cm^2')


Let’s plot the result compared to the initial cross sections at two water mixing ratios

[59]:

ptest=100.
ttest=500.
fig,axs=plt.subplots(1,2,sharex=True,sharey=True)
database['H2O'].plot_spectrum(axs[0],p=ptest,t=ttest,g=1.,x=x_array[0],label='H2O')
database['CO2'].plot_spectrum(axs[0],p=ptest,t=ttest,g=1.,x=bg_gas['CO2'],label='CO2')
IR_mix_var_gas_ktab.plot_spectrum(axs[0],p=ptest,t=ttest,x=x_array[0],g=1.,label='mix/IR')
VI_mix_var_gas_ktab.plot_spectrum(axs[0],p=ptest,t=ttest,x=x_array[0],g=1.,label='mix/VI')
database['H2O'].plot_spectrum(axs[1],p=ptest,t=ttest,g=1.,x=x_array[-1],label='H2O')
database['CO2'].plot_spectrum(axs[1],p=ptest,t=ttest,g=1.,x=bg_gas['CO2'],label='CO2')
IR_mix_var_gas_ktab.plot_spectrum(axs[1],p=ptest,t=ttest,x=x_array[-1],g=1.,label='mix/IR')
VI_mix_var_gas_ktab.plot_spectrum(axs[1],p=ptest,t=ttest,x=x_array[-1],g=1.,label='mix/VI')
axs[0].set_title('$x_{H2O}=10^{-7}$')
axs[1].set_title('$x_{H2O}=10^{-1}$')
for ax in axs :
ax.set_ylim(bottom=1.e-40)
ax.set_yscale('log')
ax.set_xscale('log')
ax.legend()
fig.tight_layout()


## Adding the opacity of a trace species to an existing corr-k table¶

If you want to add the opacity of a given molecule to a corr-k table (whether it has variable gas or not), it is rather simple. Just load your table and add the species you want:

[60]:

IR_mix_var_gas_ktab=xk.Ktable5d(path='data/corrk/lmdz/new_earth_h2o-var',res='39x19',band='IR',
kdata_unit='m^2',p_unit='Pa')
co2=xk.Ktable('CO2_*R300_0.3-50mu.ktable.SI', remove_zeros=True)
co2.remap_logPT(t_array=IR_mix_var_gas_ktab.tgrid, logp_array=IR_mix_var_gas_ktab.logpgrid)
co2.bin_down(wnedges=IR_mix_var_gas_ktab.wnedges)
earth_plus_co2=IR_mix_var_gas_ktab+1.e-1*co2
print(earth_plus_co2)


file         : data/corrk/lmdz/new_earth_h2o-var
molecule     : H2O
p grid       : [1.e+00 1.e+01 1.e+02 1.e+03 1.e+04 1.e+05]
p unit       : Pa
t grid   (K) : [200. 300. 400. 500. 600.]
wn grid      : [ 231.56880099  255.92310439  282.83877223  312.58518557  345.46005652
381.79240782  421.94586586  466.32229996  515.36584436  569.56734336
629.46926376  695.67112414  768.83549494  849.69462979  939.05779409
1037.81936442 1146.96777977 1267.59543418 1400.90960974 1548.24455953
1711.07486126 1891.03017532 2089.91155497 2309.7094719  2552.62373755
2821.08551953 3117.78167359 3445.68163456 3808.06713546 4208.56505219
4651.18370251 5140.35296264 5680.96860295 6278.44128649 6938.75072067
7668.50550426 8475.00926841 9366.33377387 9917.13029426]
wn unit      : cm^-1
kdata unit   : m^2/molecule
weights      : [0.008807   0.02030071 0.03133602 0.04163837 0.05096506 0.05909727
0.06584432 0.07104805 0.07458649 0.07637669 0.07637669 0.07458649
0.07104805 0.06584432 0.05909727 0.05096506 0.04163837 0.03133602
0.02030071 0.008807  ]
x      (vmr) : [1.e-07 1.e-01]
data oredered following p, t, x, wn, g
shape        : [ 6  5  2 39 20]


[61]:

ptest=100.
ttest=500.
x_array=IR_mix_var_gas_ktab.xgrid
fig,ax=plt.subplots(1,1,sharex=True,sharey=True)
IR_mix_var_gas_ktab.plot_spectrum(ax,p=ptest,t=ttest,x=x_array[1],g=1.,label='mix')
co2.plot_spectrum(ax,p=ptest,t=ttest,x=0.1,g=1.,label='CO2')
earth_plus_co2.plot_spectrum(ax,p=ptest,t=ttest,x=x_array[1],g=1.,label='mix+0.01*CO2')
ax.set_ylim(bottom=1.e-40)
ax.set_yscale('log')
ax.set_xscale('log')
ax.legend()
fig.tight_layout()


Important

When starting from a mix with a variable gas:

• The opacity from the background and variable gases cannot be isolated.

• The values of the array for the vmr of the variable gas (self.xgrid) are not modified (diluted) during the mixing operation.

Therefore, this treatment is valid only if the mixing ratio of the other species is << 1 (you are adding a trace gas).

For this reason, it is always more accurate to start back from the data for the pure species when they are available.

# Creating corr-k from High_resolution spectra¶

Important

Many cells in this part of the tutorial are not executable because the data needed to run them are not in the sample data directory provided.

## A note on ascii vs hdf5 data¶

The native ascii format for many high-resolution spectra (e.g. Kspectrum) is highly inefficient and takes a lot of space. Some Kspectrum files take more than a minute to read where the hdf5 equivalent would take a second. Keeping it will be the bottleneck of all the following computation.

So we recommend you seriously consider converting your files to hdf5 before anything else.

If you have a bunch of file called k001, k002, k003, etc, k0N (let’s take N=63 for exemple) in path_in you might want to use something like

path_in='/Users/jleconte/root/data2/jaziri/exo_k/spectrum_data/o3_step_band/'
for i in np.arange(1,64):
file='k'+str(i).zfill(3)
xk.convert_kspectrum_to_hdf5(path_in+file,path_out+file+'.h5', wn_column=None,
kdata_column=None, data_type='abs_coeff',
file_kdata_unit='unspecified', kdata_unit='unspecified', skiprow=0)


The various available arguments are:

• file_kdata_unit: Specifies the unit for the opacity data in the file. The type of quantity may differ whether we are handling cross sections (surface) or absorption coefficients (inverse length)

• data_type: ‘xsec’ or ‘abs_coeff’. Whether the data read are cross-sections or absorption coefficients. Make sure you specify this if you changed the default kdata_column.

• wn_column/data_column: Number of column to be read for wavenumber and kdata in python convention (0 is first, 1 is second, etc.)

• kdata_unit: Unit to convert to.

• mult_factor: A multiplicative factor that can be applied to kdata (for example to correct for any dilution effect, or specific conversion).

• skiprows: Number of header lines to skip. For the latest Kspectrum format, the header is skipped automatically.

Important

here again, with HDF5 files, we are tracking the units and data type (xsec or abs_coeff). So conversion to the right units can be made here, or at the k-coeff creation stage.

For example, if you want to convert your absorption coefficients to a weird unit, you can do

path_in='/Users/jleconte/root/data2/jaziri/exo_k/spectrum_data/o3_step_band/'
for i in np.arange(1,64):
file='k'+str(i).zfill(3)
xk.convert_kspectrum_to_hdf5(path_in+file,path_out+file+'.h5', data_type='abs_coeff',
file_kdata_unit='m^-1', kdata_unit='cm^-1')


If you want to use the cross-section of a given species in a Ksectrum file (4nd column for example=>index=3) instead of the absorption coefficient and convert them to mks, you can do

path_in='/Users/jleconte/root/data2/jaziri/exo_k/spectrum_data/o3_step_band/'
for i in np.arange(1,64):
file='k'+str(i).zfill(3)
xk.convert_kspectrum_to_hdf5(path_in+file,path_out+file+'.h5', wn_column=None,
kdata_column=3, data_type='xsec', skiprow=0,
file_kdata_unit='cm^2/molecule', kdata_unit='m^2/molecule')


Note the last kdata_unit is superfluous if you have set exo_k.Settings().set_mks(True)

## Organizing your spectra into a grid¶

Now that you spent a few cpu x hours to compute some high resolution spectra, you would like to create a corr-k matrix with it.

For this you’ll need to provide the logpgrid and tgrid (and potentially xgrid) arrays onto which your spectra have been computed.

Warning

By default, log pressures are specified in Pa in logpgrid!!! If you want to use another unit, do not forget to specify it with the grid_p_unit keyword.

You’ll also need to provide the corresponding spectra filenames in an array of shape (logpgrid.size, tgrid.size (, xgrid.size)) with the following order:

[['spectrum for logP1,T1', 'spectrum for logP1,T2'],
['spectrum for logP2,T1', 'spectrum for logP2,T2'],
['spectrum for logP3,T1', 'spectrum for logP3,T2']]


or

[[['spectrum for logP1,T1,x1', 'spectrum for logP1,T1,x2',
'spectrum for logP1,T1,x3'],
['spectrum for logP1,T2,x1', 'spectrum for logP1,T2,x2',
'spectrum for logP1,T2,x3']],

[['spectrum for logP2,T1,x1', 'spectrum for logP2,T1,x2',
'spectrum for logP2,T1,x3'],
['spectrum for logP2,T2,x1', 'spectrum for logP2,T2,x2',
'spectrum for logP2,T2,x3']],

[['spectrum for logP3,T1,x1', 'spectrum for logP3,T1,x2',
'spectrum for logP3,T1,x3'],
['spectrum for logP3,T2,x1', 'spectrum for logP3,T2,x2',
'spectrum for logP3,T2,x3']]]


### Kspectrum files generated for LMDZ¶

If you are using Kspectrum files generated for LMDZ, your files should be called k001, k002, etc.

In the case of 4 pressure points, and 2 Temperature points, you can create the list of files in the right order as follows:

[62]:

xk.create_fname_grid_Kspectrum_LMDZ(4,2,suffix='')

[62]:

array([['k001', 'k002'],
['k003', 'k004'],
['k005', 'k006'],
['k007', 'k008']], dtype='<U4')


If you have a variable gas with 3 different volume mixing ratios, then use

[63]:

xk.create_fname_grid_Kspectrum_LMDZ(4, 2, 3, suffix='')

[63]:

array([[['k001', 'k009', 'k017'],
['k003', 'k011', 'k019'],
['k005', 'k013', 'k021'],
['k007', 'k015', 'k023']],

[['k002', 'k010', 'k018'],
['k004', 'k012', 'k020'],
['k006', 'k014', 'k022'],
['k008', 'k016', 'k024']]], dtype='<U4')


And if you have converted them to hdf5, as recommended above, you can then use

[64]:

xk.create_fname_grid_Kspectrum_LMDZ(4, 2, 3, suffix='.h5')

[64]:

array([[['k001.h5', 'k009.h5', 'k017.h5'],
['k003.h5', 'k011.h5', 'k019.h5'],
['k005.h5', 'k013.h5', 'k021.h5'],
['k007.h5', 'k015.h5', 'k023.h5']],

[['k002.h5', 'k010.h5', 'k018.h5'],
['k004.h5', 'k012.h5', 'k020.h5'],
['k006.h5', 'k014.h5', 'k022.h5'],
['k008.h5', 'k016.h5', 'k024.h5']]], dtype='<U7')


### Files named with information on P, T (and x) in the title¶

You can create your grid of filenames by hand. But if the names are regular enough, the following function may help you.

create_fname_grid() for details on arguments and options.

Warning

in create_fname_grid(), the walues in logpgrid, tgrid, and xgrid will be converted to string to be put in the filename. Using integers therefore provides more control to the way your numbers will look, but floats can be used anyway.

See the 3 examples below:

[65]:

logpgrid=[1,10]
tgrid=[100,150,200]
file_grid=xk.create_fname_grid('spectrum_CO2_{p}Pa_{t}K.h5',
logpgrid=logpgrid, tgrid=tgrid, p_kw='p', t_kw='t')
file_grid

[65]:

array([['spectrum_CO2_1Pa_100K.h5', 'spectrum_CO2_1Pa_150K.h5',
'spectrum_CO2_1Pa_200K.h5'],
['spectrum_CO2_10Pa_100K.h5', 'spectrum_CO2_10Pa_150K.h5',
'spectrum_CO2_10Pa_200K.h5']], dtype='<U25')

[66]:

logpgrid=[1,2]
tgrid=[100,150,200]
xgrid=[1.e-2,1.e-1]
file_grid=xk.create_fname_grid('O3-H2O_1e{logp}Pa_{t}K_vmr{x}.h5',
logpgrid=logpgrid, tgrid=tgrid, xgrid=xgrid, p_kw='logp', t_kw='t', x_kw='x')
file_grid

[66]:

array([[['O3-H2O_1e1Pa_100K_vmr0.01.h5', 'O3-H2O_1e1Pa_100K_vmr0.1.h5'],
['O3-H2O_1e1Pa_150K_vmr0.01.h5', 'O3-H2O_1e1Pa_150K_vmr0.1.h5'],
['O3-H2O_1e1Pa_200K_vmr0.01.h5', 'O3-H2O_1e1Pa_200K_vmr0.1.h5']],

[['O3-H2O_1e2Pa_100K_vmr0.01.h5', 'O3-H2O_1e2Pa_100K_vmr0.1.h5'],
['O3-H2O_1e2Pa_150K_vmr0.01.h5', 'O3-H2O_1e2Pa_150K_vmr0.1.h5'],
['O3-H2O_1e2Pa_200K_vmr0.01.h5', 'O3-H2O_1e2Pa_200K_vmr0.1.h5']]],
dtype='<U28')

[67]:

logpgrid=[1,10,100,1000]
tgrid=np.array([100.,150.,200.])
file_grid=xk.create_fname_grid('Kspec_{logp}Pa_{t}K.h5',
logpgrid=logpgrid, tgrid=tgrid, p_kw='logp', t_kw='t')
file_grid

[67]:

array([['Kspec_1Pa_100.0K.h5', 'Kspec_1Pa_150.0K.h5',
'Kspec_1Pa_200.0K.h5'],
['Kspec_10Pa_100.0K.h5', 'Kspec_10Pa_150.0K.h5',
'Kspec_10Pa_200.0K.h5'],
['Kspec_100Pa_100.0K.h5', 'Kspec_100Pa_150.0K.h5',
'Kspec_100Pa_200.0K.h5'],
['Kspec_1000Pa_100.0K.h5', 'Kspec_1000Pa_150.0K.h5',
'Kspec_1000Pa_200.0K.h5']], dtype='<U22')


## Actually making the Ktable¶

Once we have our grids of pressure (in Pa), temperature (, vmr), and the grid of corresponding filenames, producing the Ktable is done like this

temp=np.arange(100.,401.,50)
logpres=np.arange(-1,8)
file_grid=xk.create_fname_grid_Kspectrum_LMDZ(logpres.size, temp.size, suffix='.h5')
path='data/hires/yassin_O3_hdf5/'

#that's where you choose your wavenumber bin edges grid
wnedges=xk.wavenumber_grid_R(10.,10.1,980)
tmp=xk.hires_to_ktable(path=path, filename_grid=file_grid,
logpgrid=logpres, tgrid=temp, wnedges=wnedges, mol='O3')
tmp.remove_zeros() #if you do not want to have zeros in the table for interpolation purposes
print(tmp)


With these options, you create a table with 20 gauss point following the Gausse-Legendre quadrature.

exo_k.util.user_func.hires_to_ktable() is a function provided for the user that uses the exo_k.ktable.Ktable.hires_to_ktable() and exo_k.ktable.Ktable5d.hires_to_ktable() methods. A list of all the arguments and options can be found there.

Important

Whether data are in cross section or absorption coefficient form in the input data, they will be converted to cross sections to be put in the Ktable. As this conversion is done using the pressure grid you provide, it is important to specify its unit with the grid_p_unit keyword if you are NOT using Pa.

Hereafter are some variations around the same theme to show some commented examples of use. Pick the one that fits better to your use.

temp=np.arange(100.,401.,50)
logpres=np.arange(-1,8)
path='data/hires/yassin_O3_hdf5/'

#that's where you choose your wavenumber bin edges grid
wnedges=xk.wavenumber_grid_R(10.,10.1,980)
tmp=xk.hires_to_ktable(path=path,
# because files organized as k001, k002, etc.
#  is the default, the above result can be
#  obtained without specifying a filename_grid.
#  You'll need to provide suffix here, however
suffix='.h5',
logpgrid=logpres, tgrid=temp, wnedges=wnedges, mol='O3')
tmp.remove_zeros() #if you do not want to have zeros in the table for interpolation purposes
print(tmp)


Here starting from ascii files and doing directly the conversion from abs_coeff to cross sections.

temp=np.arange(100.,401.,50)
logpres=np.arange(-1,8)
path='data/hires/yassin_O3_hdf5/'
wnedges=xk.wavenumber_grid_R(10.,10.1,980)
tmp=xk.hires_to_ktable(path=path, data_type='abs_coeff',
logpgrid=logpres, tgrid=temp, wnedges=wnedges, mol='O3',
file_kdata_unit='m^-1', kdata_unit='m^2/molec')
tmp.remove_zeros() #if you do not want to have zeros in the table for interpolation purposes
print(tmp)


If you want to provide the g weights yourself, run this instead:

temp=np.arange(100.,401.,50)
logpres=np.arange(-1,8)
path='data/hires/yassin_O3_hdf5/'
wnedges=xk.wavenumber_grid_R(10.,10.1,980)

weights=[my g weights]

tmp=xk.hires_to_ktable(path=path, suffix='.h5',
logpgrid=logpres, tgrid=temp, wnedges=wnedges, weights=weights, mol='O3')
tmp.remove_zeros() #if you do not want to have zeros in the table for interpolation purposes


You can now save your brand new Ktable on disk with any of the writing routines (write_hdf5, write_LMDZ, etc.).

## Corr-k with a variable gas: Ktable5d¶

If you have a grid of spectra with a gas that is varying along one dimension, you just have to provide the grid of volume mixing ratios (xgrid) along with the other arguments.

hires_to_ktable will automatically detect the new dimension and create a Ktable5d instead of a regular Ktable (see 4.3).

# Dealing with cross sections: Xtable()objects¶

Cross sections are like corr-k tables, except for the fact that they do not have a “g-dimension”.

Thanks to the way numpy handles arrays and to the python inheritance system, handling a Xtable() is very similar to handling a Ktable() and you can do everything that you did before.

It is so similar that you can even load a Kdatabase() with Xtable() objects.

Just remember that the kdata array now has one less dimension.

[68]:

xsecs=xk.Xtable('H2O.R10000.TauREx.pickle', kdata_unit='m^2', p_unit='Pa')
xsecs

[68]:


molecule     : H2O
p grid       : [1.00000000e+00 2.15443469e+00 4.64158883e+00 1.00000000e+01
2.15443469e+01 4.64158883e+01 1.00000000e+02 2.15443469e+02
4.64158883e+02 1.00000000e+03 2.15443469e+03 4.64158883e+03
1.00000000e+04 2.15443469e+04 4.64158883e+04 1.00000000e+05
2.15443469e+05 4.64158883e+05 1.00000000e+06 2.15443469e+06
4.64158883e+06 1.00000000e+07]
p unit       : Pa
t grid   (K) : [ 100.  200.  300.  400.  500.  600.  700.  800.  900. 1000. 1100. 1200.
1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2200. 2400. 2600. 2800.
3000. 3200. 3400.]
wn grid      : [  666.61240713   666.67906837   666.74573628 ... 33326.66766653
33330.0003333  33333.33333333]
wn unit      : cm^-1
kdata unit   : m^2/molecule
data oredered following p, t, wn
shape        : [   22    27 39124]



Saving to a file is similar to Ktables.

[69]:

xsecs.write_hdf5('data/xsec/H2O.R10000.SI')

[70]:

ptest=1.e2
ttest=300.
fig,ax=plt.subplots(1,1,sharex=True,sharey=False)
xsecs.plot_spectrum(ax,p=ptest,t=ttest,yscale='log')
fig.tight_layout()


## Xtable: sample or bin_down¶

A unique feature of Xtables over Ktables is that there are two ways of changing the spectral grid:

* bin_down(wnedges= ): does an area-conserving binning
inside the bin edges provided as wnedges.
The final spectral grid has a dimension Nw=wnedges.size-1. Usually not recommanded (See Leconte, A&A, 2020).
* sample(wngrid= ): just interpolates the cross sections
on the given grid. The final spectral grid has a dimension Nw=wngrid.size.


## Loading Xtable with kspectrum-like spectra¶

You can also load your spectra to be used as cross sections.

For this you’ll need to provide the logpgrid and tgrid arrays onto which your spectra have been computed, as well as the corresponding spectra filenames in an array of shape (logpgrid.size,tgrid.size)

path='data/hires/yassin_O3_hdf5/'
temp=np.arange(100.,401.,50)
logpgrid=np.arange(-1,8)
hires_o3_xsec=xk.hires_to_xtable(path=path,filename_grid=None, suffix='.h5',
logpgrid=logpgrid,tgrid=tgrid,
mol='O3',kdata_unit='m^2')
hires_o3_xsec.remove_zeros()


You can then save it to a hdf5 file with write_hdf5().

## Creating a Ktable from an Xtable¶

As an Xtable represents monochromatic cross-sections, it is only natural that we should be able to produce Ktable with them. And we indeed can with the Ktable.xtable_to_ktable() method, or directly from the Ktable() constructor method with the xtable= keyword.

[71]:

R=20.
xsecs_h2o=xk.Xtable('H2O.R10000.TauREx.pickle', kdata_unit='m^2', p_unit='Pa')

ktab_h2o=xk.Ktable(xtable=xsecs_h2o, wnedges=xk.wavenumber_grid_R(1000.,20000.,R))

p_plot=1.e0
t_plot=300.
fig,ax=plt.subplots(1,1,sharex=False,sharey=False)
xsecs_h2o.plot_spectrum(ax, p=p_plot, t=t_plot, yscale='log', xscale='log', label='Original data')
ktab_h2o.plot_spectrum(ax, p=p_plot, t=t_plot, g=1., label='corr-k, g=1')
ktab_h2o.plot_spectrum(ax, p=p_plot, t=t_plot, g=0., label='corr-k, g=0')
ax.legend()
fig.tight_layout()


Warning

Although in theory there is no formal criterion separating a cross section table and a set of high-resolution spectra, in practice, cross section table often have a moderate resolution. Its resolution should be high enough to sample individual lines for a cross-section table to be suitable for the computation of accurate k-coefficients.

# Dealing with continua: Cia_table() and CIAdatabase() objects¶

Cia tables are like cross section tables, except for the fact that they do not have a pressure dimension.

They can still be loaded to be converted from one format to another (in particular HITRAN .cia files or water CKD files to HDF5), from one unit to another (our preferred unit is m^5 / molec^2 ).

Finally, they will also be usefull to account for continua in opacity calculations and forward models below

But first, as for Ktables, we need to setup the path to tell exo_k where to look for cia tables.

[104]:

xk.Settings().add_search_path('data/cia', path_type = 'cia')


Then you can load a CIAdatabase by providing the list of molecules you want to look for. It will try to load CIA data for as many pairs as it can find with these molecules.

[73]:

mols=['He','H2','N2','H2O']
cia_data

[73]:

The available molecule pairs are:
All tables do NOT have common spectral grid
You will need to run sample before using the database


You can filter the files you want to look for by adding a string filter as the first argument(s)

[74]:

mols=['He','H2','N2','H2O']
cia_data

[74]:

The available molecule pairs are:
All tables do NOT have common spectral grid
You will need to run sample before using the database


If you only want a subset of these pairs, you can directly specify the pairs you want to look for with the molecule_pairs keyword.

[75]:

mol_pairs=[['He','H2'],['N2','H2']]
cia_data

[75]:

The available molecule pairs are:
All tables do NOT have common spectral grid
You will need to run sample before using the database

[76]:

xk.Settings().set_mks(True)
mols=['He','H2','N2']
ttest=300.
fig,axs=plt.subplots(2,1,sharex=False,sharey=False)
cia_data['H2']['H2'].plot_spectrum(axs[0],t=ttest,x_axis='wns',yscale='log',label='H2-H2')
cia_data['H2']['He'].plot_spectrum(axs[0],t=ttest,x_axis='wns',label='H2-He')
cia_data['N2']['H2'].plot_spectrum(axs[1],t=ttest,x_axis='wns',yscale='log',label='N2-H2')
cia_data['N2']['N2'].plot_spectrum(axs[1],t=ttest,x_axis='wns',yscale='log',label='N2-N2')
for ax in axs:
ax.legend()
fig.tight_layout()


Because CIA are often on very inhomogeneous grids (as above), you will need to sample them on the grid of your Ktable or Xtable to use them.

[77]:

ref_wns_grid=h2o_ktab_SI.wns
cia_data.sample(ref_wns_grid)
fig,axs=plt.subplots(2,1,sharex=False,sharey=False)
cia_data['H2']['H2'].plot_spectrum(axs[0],t=ttest,x_axis='wns',yscale='log',label='H2-H2')
cia_data['H2']['He'].plot_spectrum(axs[0],t=ttest,x_axis='wns',label='H2-He')
cia_data['N2']['H2'].plot_spectrum(axs[1],t=ttest,x_axis='wns',yscale='log',label='N2-H2')
cia_data['N2']['N2'].plot_spectrum(axs[1],t=ttest,x_axis='wns',yscale='log',label='N2-N2')
for ax in axs:
ax.legend()
fig.tight_layout()


The method:

CIAdatabase.cia_cross_section(pressure,temperature,composition)


can be used to compute the effective CIA cross-section for the whole gas (in m^2 per average molecule)

[78]:

compo={'H2':0.2,'He':0.1,'N2':0.4}
raw=cia_data.cia_cross_section([1.,3.],[300.,500.],compo)
fig,axs=plt.subplots(2,1,sharex=False,sharey=False,figsize=(6,5))
cia_data['H2']['H2'].plot_spectrum(axs[0],t=ttest,x_axis='wns',yscale='log')
cia_data['H2']['He'].plot_spectrum(axs[0],t=ttest,x_axis='wns')
cia_data['N2']['H2'].plot_spectrum(axs[0],t=ttest,x_axis='wns',yscale='log')
axs[1].plot(cia_data.wns,raw[0])
axs[1].set_yscale('log')
axs[1].set_ylabel('Cross section (m^2/molec)')
axs[1].set_xlabel('Wavenumber')
fig.tight_layout()


# How to use exo_k as an opacity interpolator inside your code: The Gas_mix object¶

Exo_k can also be used to directly access and interpolate opacity within your code. For this, one possibility is to load and access individual K/X/CIA tables and add all the opacities by yourself.

Another is to use the Gas_mix object: it contains the information on the pressure/temperature conditions of the gas, its composition, and can be connected to radiative databases to compute the opacity of the gas from all contributions.

## Composition and background gas¶

An empty Gas_mix object can be instanciated as follows and filled later on.

[79]:

gas=xk.Gas_mix()
print(gas)

Volume mixing ratios in Gas_mix:
inactive_gas->1.0



The message shows that, at this time, the Gas_mix is empty. Effectively, it is filled with ‘inactive_gas’.

Composition is specified using a dictionary where molecules are the keys and the volume mixing ratios as values. The rest of the Gas_mix is filled with inactive gas so that Sum(vmr)=1.

[80]:

gas.set_composition({'CO2':0.01})
print(gas)
gas.set_composition({'N2':0.8})  # reset the whole composition
print(gas)

Volume mixing ratios in Gas_mix:
CO2->0.01
inactive_gas->0.99

Volume mixing ratios in Gas_mix:
N2->0.8
inactive_gas->0.19999999999999996



One of the gases can have a vmr equal to ‘background’. In this case, the vmr of this gas will always be updated to be 1.-Sum(other vmrs).

[81]:

gas=xk.Gas_mix({'CO2':0.01,'N2':'background'}) # can also directly instantiate with a composition
print(gas)
gas['O2']=0.2  # only adds a gas or changes its vmr.
print(gas)
gas['H2O']=[0.01, 0.001]  #  Accepts arrays
print(gas)

Volume mixing ratios in Gas_mix:
CO2->0.01
N2->0.99
inactive_gas->0.0

Volume mixing ratios in Gas_mix:
CO2->0.01
N2->0.79
inactive_gas->0.0
O2->0.2

Volume mixing ratios in Gas_mix:
CO2->0.01
N2->[0.78  0.789]
inactive_gas->0.0
O2->0.2
H2O->[0.01  0.001]



## Getting the molar mass¶

The molar mass (in kg/mol) of your mix can be computed simply using:

[82]:

gas.set_composition({'N2':'background','CO2':340.e-6,'H2O':0.001})
print(gas.molar_mass())
gas.set_composition({'CH3COOH':'background'})
print(gas.molar_mass())
gas.set_composition({'Mg3SiO4':'background'})
print(gas.molar_mass())
gas=xk.Gas_mix({'N2':'background','CO2':340.e-6,'H2O':[0.001,0.1]})
print(gas.molar_mass())

0.028008920548799998
0.06005256
0.1649981
[0.02800892 0.0270191 ]


## Connecting a gas with a radiative database¶

To connect with radiative data, first you need to load the databases you want to use. You can combine a CIAdatabase with a Kdatabase filled with either Ktables or Xtables (but not both).

Then these database can be added to a Gas_mix, or you can directly instantiate with these databases. This only needs to be done once.

Warning

With Gas_mix, radiative databases need to be converted in SI units to be used. The simplest way to achieve that is to use exo_k.Settings().set_mks(True).

[83]:

database=xk.Kdatabase(['CO2','H2O'],'R300_0.3-50mu.ktable.SI')

gas.set_k_database(database)

#Loads CIA and sample them on the right spectral grid.
cia_data.sample(database.wns)

gas.set_cia_database(cia_data)
print(gas)

Volume mixing ratios in Gas_mix:
N2->[0.99866 0.89966]
CO2->0.00034
H2O->[0.001 0.1  ]
inactive_gas->0.0



Or everything can be done at initialization

[84]:

database=xk.Kdatabase(['CO2','H2O'],'R300_0.3-50mu.ktable.SI')

#Loads CIA and sample them on the right spectral grid.
cia_data.sample(database.wns)

gas=xk.Gas_mix(composition={'N2':'background','CO2':340.e-6,'H2O':[0.001,0.1]},
k_database=database, cia_database=cia_data)
print(gas)

Volume mixing ratios in Gas_mix:
N2->[0.99866 0.89966]
CO2->0.00034
H2O->[0.001 0.1  ]
inactive_gas->0.0



Now, the effective cross section on the gas can be computed as many times as necessary by calling the Gas_mix.cross_section() method providing log pressure (in Pa), temperature, and composition.

[85]:

logp_array=6.
t_array=300.
gas.set_composition({'CO2':0.01,'N2':'background'})
# composition can be set once and for all and only logP and T or changed at call time.
kdata=gas.cross_section(logp_array=logp_array, t_array=t_array)
fig,ax=plt.subplots(1,1,sharey=False,figsize=(8,4))
ax.plot(10000./gas.wns,kdata[0,:,-1])
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylabel('Cross section (m^2/molec)')
ax.set_xlabel('Wavelength (micron)')
fig.tight_layout()

[86]:

logp_array=2.
t_array=300.
composition={'N2':'background','CO2':340.e-6,'H2O':1.e-3}
kdata=gas.cross_section(composition=composition, # composition can also be specified at call time
logp_array=logp_array, t_array=t_array)

fig,ax=plt.subplots(1,1,sharey=False,figsize=(8,4))
ax.plot(10000./gas.wns,kdata[0,:,-1])
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylabel('Cross section (m^2/molec)')
ax.set_xlabel('Wavelength (micron)')
fig.tight_layout()


If you need to call this function for many logP-T points, it is more efficient to call the function once, providing logP-T (and even vmr) arrays

[87]:

logp_array=np.linspace(1.,6.,100)
t_array=300.*np.ones(logp_array.size)
composition={'N2':'background','CO2':340.e-6,
'H2O':0.001*np.ones(logp_array.size)}
kdata=gas.cross_section(composition=composition,
logp_array=logp_array, t_array=t_array)
print(kdata.shape)

(100, 1538, 20)


## Choosing gas recombination method (when using Ktables)¶

Two recombination methods are available when several gases are present:

• Simple addition: Fast. It is ok at moderate resolution (and if you are not looking for ppm level accuracy), but can be inacurrate if you are using large bins that encompass molecular bands of several molecules. This is the default mode for Gas_mix.cross_section().

• Random overlap: Very accurate in many situations, but a bit longer, especially if you have many molecules. This is the default mode for regular Ktable recombinations.

For Gas_mix.cross_section(), the behavior can be controlled with the random_overlap=True/False keyword.

[88]:

logp_array=2.
t_array=300.
gas.set_composition({'N2':'background','CO2':340.e-6,'H2O':0.001})
gas.set_logPT(logp_array=logp_array, t_array=t_array)

kdata=gas.cross_section()
kdata_randoverlap=gas.cross_section(random_overlap=True)

%timeit -r 10 -n 10 gas.cross_section()
%timeit -r 10 -n 10 gas.cross_section(random_overlap=True)

fig,ax=plt.subplots(1,1,sharey=False,figsize=(8,4))
ax.plot(10000./gas.wns,kdata_randoverlap[0,:,-1], label='Rand. Overlap')
ax.set_yscale('log')
ax.set_xscale('log')
ax.legend()
ax.set_xlabel('Wavelength (micron)')
ax.set_ylabel('Cross section (m^2/molec)')
fig.tight_layout()

3.52 ms ± 100 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)
49 ms ± 6.47 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)


You can try any molecule with a regular Element Number of atoms notation.

# The Atmclass: An integrated atmospheric radiative-transfer model¶

Whether you want to test various assumptions (will my resolution be enough, are my initial data adequate, what weights should I use, etc.) or directly simulate the transmission/emission spectrum of a planet, we implemented a simple 1d atmosphere model so that you can do just that.

This revolves around the Atm() class, which defines an atmosphere with Pressure, Temperature, and composition profiles. In essence, we are using the structure of the Gas_mix object and add the additional attributes that make an atmosphere more than just some gas with a given pressure, temperature, and composition: a gravity, a planetary radius, etc.

As Xtable() objects can be used in a Kdatabase(), you can do the equivalent of Line by Line radiative transfer.

If you want to start from very high-res spectra, first load them into a Xtable(), and then into a Database.

Check out exo_k.atm.Atm and its base class exo_k.atm.Atm_profile for details on arguments and options.

## Emission¶

[89]:

database=xk.Kdatabase(['CO2'],'R300_0.3-50mu.ktable.SI')

T0=200.
atm_ck=xk.Atm(psurf=640.e0, ptop=1.e-1, Tsurf=T0, Tstrat=100., grav=3.4,
composition={'CO2':'background'}, Nlev=30, k_database=database)

FluxTop=atm_ck.emission_spectrum(integral=True, wl_range=[0.1,150.0])
SurfFlux=atm_ck.surf_bb(integral=True)
TopBB=atm_ck.top_bb(integral=True)
print("Forward model duration:")
%timeit -n 2 -r 10 atm_ck.emission_spectrum(integral=True, wl_range=[0.1,150.0])

Forward model duration:
132 ms ± 17.5 ms per loop (mean ± std. dev. of 10 runs, 2 loops each)


The integral=True keyword specifies that the source function is integrated inside each bin. If False, the planck function is just evaluated at the center of the bin. It is faster, and ok for thin bins, but inacurate for large ones.

The wl_range range keyword allows you to limit the wavelength range onto wich you want the calculation to be done to save time.

Important

The result of emission_spectrum() is a Spectrum object onto wich some methods can be used like: computing the total flux, plotting, binning, and some algeabric operations.

[90]:

print(FluxTop.total,'W/m^2 (BB at tsurf over the computed nu range is ',SurfFlux.total,'W/m^2)')
print('(sigma*T**4 is ',xk.SIG_SB*T0**4,'W/m^2)')

77.34575509180519 W/m^2 (BB at tsurf over the computed nu range is  82.97762550596835 W/m^2)
(sigma*T**4 is  90.725872 W/m^2)

[91]:

fig,axs=plt.subplots(1,2,sharey=False,figsize=(8,4))
SurfFlux.plot_spectrum(axs[0],label='Blackbody ($T_{surf}$)')
TopBB.plot_spectrum(axs[0],label='Blackbody ($T_{strat}$)')
FluxTop.plot_spectrum(axs[0],label='OLR')
SurfFlux.plot_spectrum(axs[1],per_wavenumber=False,label='Blackbody ($T_{surf}$)')
TopBB.plot_spectrum(axs[1],per_wavenumber=False,label='Blackbody ($T_{strat}$)')
FluxTop.plot_spectrum(axs[1],per_wavenumber=False,label='OLR')

axs[0].set_ylabel('Flux ($W/m^2/cm^{-1}$)')
axs[1].set_ylabel('Flux ($W/m^2/\mu m$)')
for ax in axs:
ax.legend()
fig.tight_layout()


You can also bin it down, and when spectra have the same wavenumber grid, perform algeabric operations)

[92]:

R=100
wavenumber_grid=xk.wavenumber_grid_R(200.,10000.,R)
binned_spec=FluxTop.bin_down_cp(wavenumber_grid)

fig,axs=plt.subplots(1,2,sharey=True)
FluxTop.plot_spectrum(axs[0],label='OLR')
binned_spec.plot_spectrum(axs[0],label='binned OLR')
(SurfFlux-FluxTop).plot_spectrum(axs[1],label='Surf Blackbody - OLR')

for ax in axs:
ax.set_ylabel('Flux ($W/m^2/cm^{-1}$)')
ax.legend()
fig.tight_layout()

[93]:

origin = 'upper'
delta = 0.025

x =database.wls
y =atm_ck.play[:]
X, Y = np.meshgrid(x, y)
Z = atm_ck.exp_minus_tau()

nr, nc = Z.shape

fig, ax = plt.subplots(constrained_layout=True)
ax.invert_yaxis()
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(left=0.3,right=20.)
CS = ax.contourf(X, Y, Z, 10, cmap=plt.cm.bone, origin=origin)
ax.set_title('Atm. Transmittance')
ax.set_xlabel('Wavelength (mu)')
ax.set_ylabel('Log(Pressure/Pa)')
# Make a colorbar for the ContourSet returned by the contourf call.
cbar = fig.colorbar(CS)
cbar.ax.set_ylabel('Exp(-tau)')
# Add the contour line levels to the colorbar

[93]:

Text(0, 0.5, 'Exp(-tau)')


## Emission with continuum and your own thermal profile¶

The Atm object can also be initialized with a custom logplev and tlev profile. The vmr of the molecules can also be specified as a profile.

Warning

The (log) pressure (in Pa) and temperature must be specified at the interfaces between the layers (the levels, with Nlev=Nlay+1). The surface is the bottom level. But opacities are computed at the layers midpoints. The pressure and temperature of a layer are computed using respectively the geometric and arithmetic means of the values at levels. The vmr must be specified at the layers (so logplev.size=vmr.size+1).

[94]:

database=xk.Kdatabase(['CO2','H2O'],'R300_0.3-50mu.ktable.SI')

T0=320.;Tstrat=180.
#create the level presure grid
logplev=np.linspace(-1,5,100)
#create the layer presure grid which is halfway in log
logplay=(logplev[1:]+logplev[:-1])/2.
tlev=T0*(10**logplev/1.e5)**0.28
tlev=np.where(tlev<Tstrat,Tstrat,tlev)
tlay=(tlev[1:]+tlev[:-1])/2.
# create your composition. Here, only water is variable and follows saturation.
# the composition must be given in the layers
xH2Olay=np.exp(-2257.e3*0.018/8.31* (1./tlay - 1./373.))*1.e5/10**logplay

#Loads CIA and sample them on the right spectral grid.
cia_data.sample(database.wns)

#Loads the atmosphere and computes some intermediate things
# without CIA
atm_ck=xk.Atm(logplev=logplev,tlev=tlev,grav=9.81,
composition={'CO2':4.e-4,'H2O':xH2Olay,'N2':'background'},
k_database=database)
# with CIA
atm_ck_cont=xk.Atm(logplev=logplev,tlev=tlev,grav=9.81,
composition={'CO2':4.e-4,'H2O':xH2Olay,'N2':'background'},
k_database=database, cia_database=cia_data)

# Computes emission flux for the two cases

# Computes surface BlackBody
SurfFlux=atm_ck.surf_bb(integral=True)

# Computes top of atmosphere BlackBody
TopBB=atm_ck.top_bb(integral=True)
print(FluxTop.total,'W/m^2 (BB at tsurf over the computed nu range is ',SurfFlux.total,'W/m^2)')
print('(sigma*T**4 is ',xk.SIG_SB*T0**4,'W/m^2)')

333.6292541916714 W/m^2 (BB at tsurf over the computed nu range is  566.7392171105968 W/m^2)
(sigma*T**4 is  594.5810747392 W/m^2)

[95]:

fig,axs=plt.subplots(1,2,sharey=False, figsize=(8,4))
SurfFlux.plot_spectrum(axs[0],label=r'Blackbody ($T_{surf}$)')
TopBB.plot_spectrum(axs[0],label=r'Blackbody ($T_{strat}$)')
FluxTop.plot_spectrum(axs[0],label=r'OLR w/o CIA')
FluxTop_cont.plot_spectrum(axs[0],label=r'OLR w/ CIA')
SurfFlux.plot_spectrum(axs[1],per_wavenumber=False,label=r'Blackbody ($T_{surf}$)')
TopBB.plot_spectrum(axs[1],per_wavenumber=False,label=r'Blackbody ($T_{strat}$)')
FluxTop.plot_spectrum(axs[1],per_wavenumber=False,label=r'OLR w/o CIA')
FluxTop_cont.plot_spectrum(axs[1],per_wavenumber=False,label=r'OLR w/ CIA')

axs[0].set_ylabel(r'Flux ($W/m^2/cm^{-1}$)')
axs[1].set_ylabel(r'Flux ($W/m^2/\mu m$)')
for ax in axs:
ax.legend()
#axs[0].invert_yaxis()
fig.tight_layout()


## Transmission¶

[96]:

T0=2000.;xH2O=1.e-5;Rs=1.*u.Rsun;grav=10.
database=xk.Kdatabase(['H2O'],'50mu','SI',remove_zeros=True)
cia_data.sample(database.wns)

atm_ck=xk.Atm(psurf=10.e5,ptop=1.e-4,Tsurf=T0,Tstrat=T0,grav=grav,
composition={'H2':'background','H2O':xH2O},Nlev=100,Rp=1.*u.Rjup,
k_database=database, cia_database=cia_data)
spec=atm_ck.transmission_spectrum(Rstar=Rs, rayleigh=True)
print("Forward model duration:")
%timeit -n 2 -r 5 atm_ck.transmission_spectrum(Rstar=Rs, rayleigh=True)

Forward model duration:
355 ms ± 53.2 ms per loop (mean ± std. dev. of 5 runs, 2 loops each)

[97]:

fig,ax=plt.subplots(1,1,sharex=True,sharey=True)
spec.plot_spectrum(ax,xscale='log',label='Native corr-k resolution')
spec.bin_down_cp(xk.wavenumber_grid_R(400.,30000.,200)).plot_spectrum(ax,xscale='log',label=r'$R=200$')
spec.bin_down_cp(xk.wavenumber_grid_R(400.,30000.,30)).plot_spectrum(ax,xscale='log',label=r'$R=30$')
ax.set_ylabel('Depth')
ax.legend()
fig.tight_layout()


## Heating rates¶

[98]:

T0=300.
Tstrat=300.
database=xk.Kdatabase(['CO2','H2O'],'_R300_0.3-50mu.ktable.SI')
cia_data.sample(database.wns,remove_zeros=True)

atm_ck=xk.Atm(psurf=1.e10,ptop=1.e-1,Tsurf=T0,Tstrat=Tstrat,grav=10.,
composition={'H2':'background','CO2':0.},
rcp=0.2183,Rp=1.*u.Rearth,
Nlev=30,
k_database=database,cia_database=cia_data)
HR=atm_ck.heating_rate(szangle=60.)
atm_ck.set_gas({'H2':'background','CO2':.05,'H2O':0.5})
HR2=atm_ck.heating_rate(szangle=60.)
HR3=atm_ck.heating_rate(szangle=89.99)

[99]:

fig,axs=plt.subplots(1,2,sharey=True)
axs[0].set_yscale('log')
axs[0].set_xscale('log')
axs[0].plot(HR,atm_ck.play)
axs[0].plot(HR2,atm_ck.play)
axs[0].plot(HR3,atm_ck.play)
axs[1].plot(atm_ck.tlev,atm_ck.plev)
axs[0].set_xlabel('Heating Rate (K/s)')
axs[1].set_xlabel('Temperature (K)')
axs[0].set_ylabel('Pressure (Pa)')
axs[0].invert_yaxis()
axs[0].set_xlim(left=1.e-20,right=1.)
fig.tight_layout()