Commit 047a69bd authored by Paul Mollière's avatar Paul Mollière
Browse files

Merge with retrieval package complete.

parent 1370e7ab
......@@ -11,7 +11,7 @@ necessary folder structure if you want to install high-resolution
opacities later (:math:`\lambda/\Delta\lambda=10^6`).
Thus, to get started download the `opacity and input data
<https://keeper.mpdl.mpg.de/f/4b9409d9d17d443cb6ee/?dl=1>`_
<https://keeper.mpdl.mpg.de/f/f5aba635d3a244adb3c0/?dl=1>`_
(12.1 GB), unzip them, and put the "input_data" folder somewhere on
your computer (it does not matter where).
......
......@@ -34,26 +34,26 @@ The linelists available on the `Exomol website <http://www.exomol.com>`_ were re
2. Navigate to your petitRADTRANS opacity folder, which is the ``path to input_data/opacities/lines/corr_k/`` folder.
Make a new folder called, for example, ``H2O_Chubb`` and place the ``1H2-16O__POKAZATEL__R1000_0.3-50mu.ktable.petitRADTRANS.h5`` file inside.
You're done! pRT will check automatically if an ``*.h5`` file is in the opacity folder, so as long as you keep the ``.h5`` extension, the file in the ``H2O_Chubb`` folder can be called whatever you like. The name of the folder, which we chose to be ``H2O_Chubb`` here, can also be anything. After dropping the ``.h5`` file in the ``H2O_Chubb`` folder, the opacity species ``H2O_Chubb`` is now ready for use! Note that also the abundances must be specified for ``H2O_Chubb`` when using petitRADTRANS with this opacity table.
You're done! pRT will check automatically if an ``*.h5`` file is in the opacity folder, so as long as you keep the ``.h5`` extension, the file in the ``H2O_Chubb`` folder can be called whatever you like, as long as it begins with the correct chemical formula for that species. The name of the folder, which we chose to be ``H2O_Chubb`` here, can also be anything. After dropping the ``.h5`` file in the ``H2O_Chubb`` folder, the opacity species ``H2O_Chubb`` is now ready for use! Note that also the abundances must be specified for ``H2O_Chubb`` when using petitRADTRANS with this opacity table.
.. _OWtopRT:
-------------------------------------------------------------------
Converting cross-section grids from `opacity.world`_
Converting cross-section grids from `DACE`_
-------------------------------------------------------------------
Pre-computed opacities are also available from `opacity.world`_,
Pre-computed opacities are also available from `DACE`_,
which have been generated using the method presented in `Grimm & Heng
(2015) <https://iopscience.iop.org/article/10.1088/0004-637X/808/2/182>`_ .
The website allows to download the cross-section tables as a function
of pressure and temperature.
Simply decide on any P-T and wavelength range that you are interested
Simply decide on any P-T range line list that you are interested
in. Note that their spectral coordinate is wavenumber, in units of
:math:`{\rm cm}^{-1}`.
Then transform the resulting binary files, that you downloaded from
`opacity.world`_, like described in the following. This script was
`DACE`_, like described in the following. This script was
developed together with Mantas Zilinskas, for `Ziliniskas et
al. (2020) <https://arxiv.org/abs/2003.05354>`_.
......@@ -61,25 +61,25 @@ First we load the relevant packages.
.. code-block:: python
import numpy as np
import struct
import os
from scipy.interpolate import interp1d
import math
import sys
import numpy as np
import struct
import glob
from scipy.interpolate import interp1d
import math
import sys
Then we define where the input file are and where ouput files are
supposed to be put. This example here is for the HDO opacity.
supposed to be put. This example here is for the H2O opacity.
.. code-block:: python
# Paths to opacity world files and output directory (ADJUST ACCORDINGLY)
path_to_files = 'opacity.world.s3/generateddata/1H-2H-16O__VTT_e2b/'
path_to_output = 'sigmas_adapted_pRT/'
filelist = os.listdir(path_to_files) # Find all opacity world files in the directory
# Paths to DACE files and output directory (ADJUST ACCORDINGLY)
path_to_files = '../1H2-16O__POKAZATEL_e2b/'
path_to_output = '../1H2-16O__POKAZATEL_e2b/'
filelist = glob.glob(path_to_files+'Out*') # Find all opacity world files in the directory
Opacity.world saves the opacities in units of :math:`{\rm cm}^{2}{\rm
DACE saves the opacities in units of :math:`{\rm cm}^{2}{\rm
g}^{-1}`, but the petitRADTRANS conversions scripts need :math:`{\rm
cm}^{2}`. So we will have to convert below. For this it is important
that the mass of the absorber species is defined, in units of
......@@ -88,18 +88,18 @@ amu. **Do not forget to adapt this for every new species!**
.. code-block:: python
# Properties of chosen species
species_mass = 19.
species_mass = 18.
This function below will read the binary files downloaded from
`opacity.world`_:
`DACE`_:
.. code-block:: python
def read_bin_single(filename):
def read_bin_single(filename):
""" Read a binary opacity world file.
"""
# Open file
file = open(filename,mode='rb')
# Read content
......@@ -111,15 +111,15 @@ This function below will read the binary files downloaded from
points = int(len(cont)/4)
# Create array of the appropriate length
x = np.ones(points)
# Read the binary data into the array
for i in range(int(points)):
test = struct.unpack('f',cont[i*4:(i+1)*4])
x[i] = test[0]
return x
Finally we define the function that reads the binary `opacity.world`_
Finally we define the function that reads the binary `DACE`_
files, and saves them in the format that can be used by the opacity input
generating scripts of petitRADTRANS. For this you also need the file
that defines the petitRADTRANS wavelength grid, which can be
......@@ -127,57 +127,59 @@ downloaded here: `wlen_petitRADTRANS.dat`_
.. code-block:: python
def convert():
def convert():
""" Converts opacity.world binary files for further pRT processing
"""
# Read the fiducial petitRADTRANS wavelength grid
wavelength_petit = np.genfromtxt('wlen_petitRADTRANS.dat')
for file in filelist:
# Reads oworld file
opa = read_bin_single(path_to_files + file)
# Temp and pressure for naming files
t = str(int(file.split('/')[-1].split('_')[3]))
p = str(file.split('/')[-1].split('_')[4].split('.bin')[0].replace('n','-').replace('p',' '))
p = p[:2] + '.' + p[2:]
p = str(np.round(1e1**float(p), 6))
print (t,p)
# Wavenumber points from range given in the file names
wl_start = int(file.split('/')[-1].split('_')[1])
wl_end = int(file.split('/')[-1].split('_')[2])
wlen = np.linspace(wl_start, wl_end, len(opa))
# Convert to cm or [micron]
wavelength = 1./wlen#/1e-4
# Invert them to go from a accending wavenumber ordering
# to an accending wavelength ordering.
wavelength = wavelength[::-1]
sigma = opa[::-1]
# OW opacities cm^2/g, convert to cm^2 by *species_mass*amu
sigma = sigma*species_mass*1.66053892e-24
# Interpolate
sig_interp = interp1d(wavelength, sigma,bounds_error=False,fill_value=0.0)
sig_interpolated_petit = sig_interp(wavelength_petit)
# Check if interp values are below 0 or NaN
for i in sig_interpolated_petit:
if i < 0.:
print (i)
elif math.isnan(i):
print (i)
#### SAVING REBINNED #### Around 300 MB per grid point
# New file name is 'sigma_+ temp + .K_ + Pressure + bar.dat'
np.savetxt(path_to_output + 'sigma_' + str(t) + '.K_' + str(p) + 'bar.dat',
np.column_stack((wavelength_petit, sig_interpolated_petit)))
""" Converts opacity.world binary files for further pRT processing
"""
# Read the fiducial petitRADTRANS wavelength grid
wavelength_petit = np.genfromtxt('wlen_petitRADTRANS.dat')
for file in filelist:
# Reads oworld file
opa = read_bin_single(path_to_files + file)
# Temp and pressure for naming files
t = file[16:21]
t = t.lstrip('0')
p = file[22:26]
print (t,p)
# Wavenumber points from range given in the file names
wl_start = int(file[4:9])
wl_end = int(file[10:15])
wlen = np.linspace(wl_start, wl_end, len(opa))
# Convert to cm or [micron]
wavelength = 1./wlen#/1e-4
# Invert them to go from a accending wavenumber ordering
# to an accending wavelength ordering.
wavelength = wavelength[::-1]
sigma = opa[::-1]
# OW opacities cm^2/g, convert to cm^2 by *species_mass*amu
sigma = sigma*species_mass*1.66053892e-24
# Interpolate
sig_interp = interp1d(wavelength, sigma,bounds_error=False,fill_value=0.0)
sig_interpolated_petit = sig_interp(wavelength_petit)
# Check if interp values are below 0 or NaN
for i in sig_interpolated_petit:
if i < 0.:
print (i)
elif math.isnan(i):
print (i)
#### SAVING REBINNED #### Around 300 MB per grid point
# New file name is 'sigma_+ temp + .K_ + Pressure + bar.dat'
np.savetxt(path_to_output + 'sigma_' + str(t) + '.K_' + str(p) + 'bar.dat', \
np.column_stack((wavelength_petit, \
sig_interpolated_petit)))
Then you just need to start the conversion:
......@@ -190,8 +192,8 @@ k-tables. This is done in an analogous way as explained in Section
:ref:`EXtopPRT` below. When doing this, note that you can omit the step rebinning the cross-section
files to the petitRADTRANS wavelength grid, because this was already
done in ``convert()`` above!
.. _opacity.world: http://opacity.world/
.. _DACE: https://dace.unige.ch/opacityDatabase/
The opacities can then be installed as described in Section
:ref:`install` below.
......@@ -207,7 +209,7 @@ In this example we will show you how to do this using ExoCross, the
open-source opacity calculator of the `Exomol`_ database.
ExoCross can be downloaded `here <https://github.com/Trovemaster/exocross>`_, is described in
`Yurchenko et al. (2018)`_ and documented `here
<https://exocross.readthedocs.io>`_.
<https://exocross.readthedocs.io>`_.
.. _Exomol: http://www.exomol.com
.. _Yurchenko et al. (2018): https://arxiv.org/abs/1801.09803
......@@ -221,7 +223,7 @@ gfortran. The relevant lines in "makefile" should look like this:
.. code-block:: bash
#FOR = ifort
#FFLAGS = -O3 -qopenmp -traceback -ip
#FFLAGS = -O3 -qopenmp -traceback -ip
FOR = gfortran
FFLAGS = -O2 -fopenmp -std=f2008
......@@ -319,7 +321,7 @@ This calculates the opacity of NaH with the following settings
Line (2018)`_. Also see the text around Equation 12 in `Sharp &
Burrows (2007)`_ for more information. Sometimes more detailed
broadening information is available on Exomol, `see here`_.
.. _Hartmann et al. (2002): http://adsabs.harvard.edu/abs/2002JQSRT..72..117H
.. _Grimm & Heng (2015): https://arxiv.org/abs/1503.03806
.. _Amundsen et al. (2014): https://arxiv.org/abs/1402.0814
......@@ -331,7 +333,7 @@ If more detailed broadening information is avaiable (not for NaH) you can replac
the lines below ``species`` with something like
.. code-block:: bash
species
0 gamma 0.06 n 0.5 t0 296 file path_toH2_broadening_information_file model J ratio 0.860000
1 gamma 0.06 n 0.5 t0 296 file path_toHe_broadening_information_file model J ratio 0.140000
......@@ -412,7 +414,7 @@ calculated above.
import numpy as np
from scipy.interpolate import interp1d
# Read the opacity file from ExoCross
dat = np.genfromtxt('NaH_1000K_1em5bar.out.xsec')
wavelength = 1./dat[:,0]
......@@ -461,7 +463,7 @@ folder where the opacity files are. In our simple example (just one
NaH file at 1000 K and :math:`10^{-5}` bar, its content looks like this:
.. code-block:: bash
NaH_1000K_1em5bar_petit_grid.dat
Let's start with the k-table calculation, for the low-resolution
......@@ -472,15 +474,15 @@ are interested in (NaH has 24 amu, so just put 24, like below):
.. code-block:: fortran
! (c) Paul Molliere 2014
program calc_k_g
implicit none
!-----------------------------------------------------------
! ||| ||| ||| !
! \|||/ \|||/ \|||/ !
! v v v !
! v v v !
!----------------------------------------------------------!
!----------------------------------------------------------!
! DO NOT FORGET TO CHANGE THE MASS OF THE MOLECULE
......@@ -492,11 +494,11 @@ are interested in (NaH has 24 amu, so just put 24, like below):
! /|||\ /|||\ /|||\ !
! ||| ||| ||| !
!----------------------------------------------------------!
Next, compile the Fortran source:
.. code-block:: bash
gfortran -o calc_k_g_r1000_ptrad calc_k_g_r1000_ptrad.f90
Lastly, create a folder called kappa_gs_r1000. Now, take care that the opacity files, the compiled Fortran routine,
......@@ -504,7 +506,7 @@ sigma_list.ls, retrieval_NP_16_ggrid.dat and the kappa_gs_r1000 folder
are all in the same folder. And that you are in this folder. Type
.. code-block:: bash
./calc_k_g_r1000_ptrad
and all k-tables will be generated and placed into the kappa_gs_r1000
......@@ -516,7 +518,7 @@ to have the correct molecule mass. **Do not change the wavelength boundary value
For NaH, with mass 24, it should look like this:
.. code-block:: bash
# Minimum wavelength in cm
0.3d-4
# Maximum wavelength in cm
......@@ -527,7 +529,7 @@ For NaH, with mass 24, it should look like this:
Next, compile the high-resolution opacity conversion routine:
.. code-block:: bash
gfortran -o make_short make_short.f90
Now, again take care that the opacity files, the compiled Fortran routine,
......@@ -535,7 +537,7 @@ sigma_list.ls, short_stream_lambs_mass.dat and the short_stream folder
are all in the same folder. And that you are in this folder. Type
.. code-block:: bash
./make_short
and all high resolution opacity tables will be generated and placed into the short_stream
......@@ -550,7 +552,7 @@ The new opacity files are now ready to be installed. Before that
create a file called "molparam_id.txt" with the following content
.. code-block:: bash
#### Species ID (A2) format
06
#### molparam value
......
Analysis of results
===================
OLD (from Mollière+19): Analysis of results
===========================================
Here we give some example code for how to analyse one's results after
having carried out a retrieval. Example code is given for posterior
......
Emission spectrum retrieval
===========================
OLD (from Mollière+19): Emission spectrum retrieval
===================================================
.. toctree::
:maxdepth: 2
......
Transmission spectrum retrieval
===============================
OLD (from Mollière+19): Transmission spectrum retrieval
=======================================================
.. toctree::
:maxdepth: 2
......
Retrieval Examples
==================
In this section we will give two fully imlemented examples of
retrievals, using synthetic JWST observations. One retrieval for
emission and one for transmission spectra. These are the retrievals
reported on in the `petitRADTRANS paper
<https://arxiv.org/abs/1904.11504>`_.
In addition we also give some implemented examples on how to analyze
the results, by producing contribution functions, corner plots and
pressure-temperature confidene envelopes.
In this section we will give fully imlemented examples of
retrievals. **Please use the** `petitRADTRANS Retrieval Tutorial <notebooks/pRT_Retrieval_Example.html>`_
**from now on, we still give our old setups from the** `petitRADTRANS paper
<https://arxiv.org/abs/1904.11504>`_ **though**.
.. toctree::
:maxdepth: 2
......
......@@ -12,7 +12,8 @@ of exoplanets for clear and cloudy planets. pRT also incorporates (**new!**)
an easy subpackage for running retrievals with nested sampling.
**To get started with some examples on how to run pRT immediately,
see** `"Getting started" <content/notebooks/getting_started.html>`_. **Otherwise read on for some more general info.**
see** `"Getting started" <content/notebooks/getting_started.html>`_. **Otherwise read on for some more general info,
also on how to best run retrievals.**
pRT has two different opacity treatment modes. The low resolution mode runs calculations
at :math:`\lambda/\Delta\lambda\leq 1000` using the so-called correlated-k treatment for opacities.
......@@ -29,16 +30,17 @@ spectral calculations. The user should verify whether this leads to solutions wh
described in the `API <autoapi/petitRADTRANS/radtrans/index.html#petitRADTRANS.radtrans.Radtrans>`_ here. An example is
given `here <content/notebooks/Rebinning_opacities.html>`_.
pRT can calculate transmission and emission spectra of exoplanets, for clear or cloudy atmospheres. The different cloud treatment
(gray clouds, power law clouds, "real" clouds using optical constants) are described in the tutorial `here <content/notebooks/clouds.html>`_. Scattering is
included in pRT, but must be specifically turned on for emission spectra (note that scattering increases the runtime),
pRT can calculate transmission and emission spectra of exoplanets, for clear or cloudy atmospheres. The different cloud treatments
(gray clouds, power law clouds, "real" clouds using optical constants) are described in the tutorial `here <content/notebooks/clouds.html>`_.
Scattering is included in pRT, but must be specifically turned on for emission spectra (note that scattering increases the runtime),
see `Scattering for Emission Spectra <content/notebooks/emis_scat.html>`_. pRT can also calculate the reflection of light
at the surface of rocky planets, for which the use can specify wavelength-dependent albedos and emissivities. This is likewise
at the surface of rocky planets, for which the user can specify wavelength-dependent albedos and emissivities. This is likewise
explained in `Scattering for Emission Spectra <content/notebooks/emis_scat.html>`_.
The newly added retrieval subpackage is documented `here <content/notebooks/pRT_Retrieval_Example.html>`_.
At the moment pRT retrievals are making use of the `PyMultiNest <https://johannesbuchner.github.io/PyMultiNest/>`_
package for parameter inference. Of course you are free to use pRT spectral synthesis routines with any other inference tool of your liking.
For this you will have to setup your own retrieval framework, however (you can modify / check our source code for inspiration).
petitRADTRANS is available under the MIT License, and documented in
`Mollière et al. (2019) <https://arxiv.org/abs/1904.11504>`_, for the general code, and `Mollière et al. (2020) <https://arxiv.org/abs/2006.09394>`_, Alei et al. (in prep.), for the scattering implementation. Please cite these papers if you make use of petitRADTRANS in your work.
......
......@@ -209,7 +209,7 @@ low-resolution version (<span class="math notranslate nohighlight">\(\lambda/\De
provides all relevant input files for pRT to run, and contains the
necessary folder structure if you want to install high-resolution
opacities later (<span class="math notranslate nohighlight">\(\lambda/\Delta\lambda=10^6\)</span>).</p>
<p>Thus, to get started download the <a class="reference external" href="https://keeper.mpdl.mpg.de/f/4b9409d9d17d443cb6ee/?dl=1">opacity and input data</a>
<p>Thus, to get started download the <a class="reference external" href="https://keeper.mpdl.mpg.de/f/f5aba635d3a244adb3c0/?dl=1">opacity and input data</a>
(12.1 GB), unzip them, and put the “input_data” folder somewhere on
your computer (it does not matter where).</p>
<p>Next, please add the following environment variable to your
......
......@@ -108,6 +108,7 @@
</ul>
</li>
<li class="toctree-l2"><a class="reference internal" href="pRT_Retrieval_Example.html">petitRADTRANS Retrieval Tutorial</a></li>
<li class="toctree-l2"><a class="reference internal" href="Rebinning_opacities.html">Rebinning opacities</a></li>
<li class="toctree-l2"><a class="reference internal" href="nat_cst_utility.html">Utility functions</a></li>
<li class="toctree-l2"><a class="reference internal" href="poor_man.html">Interpolating chemical equilibrium abundances</a></li>
</ul>
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment