pcubillos / pyratbay

Python Radiative Transfer in a Bayesian Framework

Pyrat Bay: Python Radiative Transfer in a Bayesian framework

A Forward-modeling and retrieval code to model exoplanet atmospheres and spectra.

Install as:

pip install pyratbay

or:

conda install -c conda-forge pyratbay

Docs at:

https://pyratbay.readthedocs.io/en/latest/

Cite as (ADS):

@ARTICLE{CubillosBlecic2021mnrasPyratBay,
       author = {{Cubillos}, Patricio E. and {Blecic}, Jasmina},
        title = "{The PYRAT BAY framework for exoplanet atmospheric modelling: a population study of Hubble/WFC3 transmission spectra}",
      journal = {\mnras},
     keywords = {radiative transfer, methods: statistical, planets and satellites: atmosphere, Astrophysics - Earth and Planetary Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics},
         year = 2021,
        month = aug,
       volume = {505},
       number = {2},
        pages = {2675-2702},
          doi = {10.1093/mnras/stab1405},
archivePrefix = {arXiv},
       eprint = {2105.05598},
 primaryClass = {astro-ph.EP},
       adsurl = {https://ui.adsabs.harvard.edu/abs/2021MNRAS.505.2675C},
      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

pyratbay's Issues

Emission integration bug.

Bug: When integrating B*exp(-tau), if tau doesn't properly sample the optical depth, the Simpson integration return bogus results. This happens when there is a layer with very large extinction coefficient, making the optical depth jump directly into tau > taumax.

I need to find a better way to integrate f(x)exp(-x)dx.
So far, adding extra points and interpolating is my best idea.
Can I at least optimize when to add these points, how many do I need, and where to put them?

Re-factor the cia reader as a general cross-section reader

Change function and variable names.
Accept a single species.
Distinguish CIA from line-absorption bands.
Etc.

Implement VALD LBL database reader

Python3

Add compatibility with Python3.

House keeping

Minor bug fixes:

• pass nproc from PB to MC3.
• warning when wavenumber sampling does not match EC table wn sampling.
• Skip path check when mode=opacity.
• skip tmin, tmax check if extfile exists.
• skip bulk check if retflag does not have mol flag.
• Use backslash for long line commands in documentation (for pdf).
• Remove checks for rpars, hpars, and cpars when parameters are in MCMC params input.
• Remove system parameters check for opacity run.

Re-factor cross-section interpolation

Include the wavenumber interpolation in the initialization step, leaving only the temperature interpolation for the RT calculation. This will make the code much faster.

Add warning for sum(q) != 1

Implemenet info() methods for the Pyrat-Bay objects

Implement a method that prints the values of an object's variables.

Add scripts to format the cross-section files

• Borysow CIA
• HITRAN CIA
• HITRAN cross section (CSX)
• ExoMol cross section (xsecs)

Add mplanet argument

Use this variable to calculate the surface gravity (which will depend on the planetary radius) when gplanet is not given.

Post Processing I: Best Fit

Figure out if I need to rerun Pyrat to get the best-fit result. Can I save the best run from the MCMC?

Implement line-transition reader for Kurucz/Schwenke TiO

Retrieval IV: MC3

Connect MC3 to Pyrat. This involve several, potentially complicated, tasks:

• Pre-processing.
• Shared memory for large files?
• MCMC hook as in the BART code? No.

Post Processing III: Plots

Incorporate all useful output plots from BART.
Think what else would be good to show.

Update binsearch in Lineread

See pylineread in the Transit package.

Replace Alkali Voigt calculation with astropy routine

Astropy's Voigt1D routine seems to be fast enough to run during the retrieval. This avoids using the coarse Voigt widths sampling, which was producing funky contribution functions.

Make a script to generate a table out of the Pyrat classes

Add error and warning messages for required and defaulted arguments

Including (update as necessary):

• atmfile not found (this must be an error).
• Defaulted atm units.
• wllow or wnup not found.
• wlup or wnlow not found.
• Defaulted path.
• rstar not found.
• Defaulted outspec.
• Default molfile.
• Default radunits, punits.
• Default gplanet.
• Isotope's species not in atm file.
• Not-found error for tmin, tmax, tstep if extgrid.
• Default vextent, Dmin, Dmax, nDop, Lmin, Lmax, nLor, DLratio.
• Default maxdepth.
• Default raygrid.
• Default wlunits, wnunits.
• wnosamp, wnstep.

VO partition

Update VO partition function with values from Barklem & Collet (2016).
(These cover temperatures below 1000 K).

Post Processing II: Contribution Functions

Incorporate Jasmina's code. I've never seen it in detail, might need several edits.

Cloud models: III

Implement constant cross-section opacity cloud model (as in Sing et al. 2016).
Define cloud top and bottom?
To which particles relate this cross section? total? large-particles?

create a pickle/npz file output?

optional?
always?
how to parse?

Fix timestamps when init() and run() are run separately.

Add Rayleigh-scattering models

• for H2 from Dalgarno & Williams (1962).
• fitting model from Lecavelier des Etangs et al. (2008).
• H and He from Dalgarno (see Kurucz 1970).

See also references in Benneke & Seager (2012).

Optimize PT posterior plot

Less memory, faster computation:

• Compute profiles from unique() view of the posterior.
• Compute percentiles per layer.

Implement CH4 database reader for Hargreaves et al.

This issue involves two tasks:

• process cross-section files.
• process line-by-line files.

Add temperature-profile model

Check that the C implementation works and add to the repo.

Retrieval II: Kurucz

Implement code to read/process the Kurucz stellar models.

Describe better makeradius

Think of a better name.
Fix the docstring to be more precise of what's going on.

Implement database opacities for Na and K

See references in Freedman et al. (2008), section 2.4.

Retrieval III: Filter Passband Files

Implement the code to process the filter passbands.

Make demo runs

For both the project's front page and wiki.

Create a separate file with the warnings

Separate from the log output so the user can see that something might have gone wrong.

Add a message at the end saying: "There were [N] warnings during the execution, see 'warnings.txt'."

Retrieval I: reloadatm

Implement the reloadatm() function. This entails:

• recalculate the temperature, abundances (balance() function). Feed it to pyrat.
• Recalculate the other atmospheric properties (radius, densities, mean mass).

enable opacity run from pressure params instead of atmfile

Currently, opacity runs require an atmospheric file. But it only takes the pressure info from it.
Enable pressure array input by defining: {pmin, pmax, nlayers} in the config file.

Parallelize the opacity grid routine

Should be simple, right?

Cloud models: II

Implement cloud deck model: a wavelength-independent sharp opaque layer.

retrieval mmm

In the retrieval atmospheric update step, I need to update the mean molecular mass after recomputing the abundances, such that the recomputed radius has the correct value.

Implement line-transition reader for Partridge and Schwenke

Cloud models: I

Implement Ackerman and Marley (2001) cloud model.

Post Processing IV: Logs

Decide between a unique log for the whole run (full MCMC pyrat-bay, pyrat, pre-processing, post-processing)? exclude/mute pyrat log at all?

Implement line-transition reader for B. Plez VO

Adjust pbay inputs' default warning/error messages

Do the calls only when the defaulted variables are requested (depending on the specific requested step: temperature, atm, pyrat, MCMC).

Implement dual screen/file writing

Enable the option to write to file at the same time that we are writing to screen.

Update the HITRAN data file.

TLI bug with multiple molecules

current version of the tli code has some bug. When you try to write several species in the same time it writes the file in some weird order and the opacity table code breaks.
So the code does not break while writing the tli file, rather when writing the opacity file.

Here is the bug that appeared to me several times with different species.
When writing one species per tli file the opacity code does not give any error.

Cheers
Jass

==========================================================
Reading line transition info.
  Read TLI file:
  './PyratBay_24July2018_HD189_HITRAN_HITEMP_H20_CO_CO2_CH4_NH3_N2_03_30um.tli'.
  TLI data storage: Little endian
  TLI version: 6.4.10.
  TLI wavenumber range (cm-1): [333.3, 30303.0]
  Number of data bases: 6
  Data base name: 'HITRAN CO'
  Molecule name: 'CO'
  Temperature range: 70.0 -- 3000.0 K.
   Isotope: 26,  mass: 27.9949 u,  isotopic ratio: 0.9865
   Isotope: 36,  mass: 28.9983 u,  isotopic ratio: 0.01108
   Isotope: 28,  mass: 29.9992 u,  isotopic ratio: 0.001978
   Isotope: :  4.446e-08
         P,  mass: 0.0000 u,  isotopic ratio: 8.9034e-10
---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-4-e0048d0c7465> in <module>()
----> 1 pyrat = pb.pbay.run('/home/jasmina/Work/3D-retrieval/Ian-HD189733b/transit-runs/re-re-run/HD189_HITEMP_HITRAN_opacity.cfg')

/home/jasmina/Work/3D-retrieval/Ian-HD189733b/pyratbay/pyratbay/pbay/driver.pyc in run(argv, main)
    118 
    119     # Initialize pyrat object:
--> 120     pyrat = py.init(args.cfile, log=log)
    121 
    122     # Compute spectrum and return pyrat object if requested:

/home/jasmina/Work/3D-retrieval/Ian-HD189733b/pyratbay/pyratbay/pyrat/driver.pyc in init(argv, main, log)
     74 
     75   # Read line database:
---> 76   rl.readlinedb(pyrat)
     77   timestamps.append(time.time())
     78 

/home/jasmina/Work/3D-retrieval/Ian-HD189733b/pyratbay/pyratbay/pyrat/readlinedb.pyc in readlinedb(pyrat)
     34            format(pyrat.linedb[n]), pyrat.log, 2)
     35     # Read headers info:
---> 36     dbindex.append(dbindex[-1] + readheader(pyrat, TLI[n]))
     37 
     38   # Set link to molecules' indices:

/home/jasmina/Work/3D-retrieval/Ian-HD189733b/pyratbay/pyratbay/pyrat/readlinedb.pyc in readheader(pyrat, linefile)
    124       dbindex[j] = i + pyrat.lt.ndb
    125       lenIsoName = pt.unpack(linefile, 1,          "h")
--> 126       name[j]    = pt.unpack(linefile, lenIsoName, "s")
    127       mass[j]    = pt.unpack(linefile, 1,          "d")
    128       ratio[j]   = pt.unpack(linefile, 1,          "d")

/home/jasmina/Work/3D-retrieval/Ian-HD189733b/pyratbay/pyratbay/tools/ptools.pyc in unpack(file, n, dtype)
    284   fmt  = "{:d}{:s}".format(n, dtype)
    285   # Calculate the number of bytes to read:
--> 286   size = struct.calcsize(fmt)
    287   # Read:
    288   output = struct.unpack(fmt, file.read(size))

error: bad char in struct format

Rename pconstants.units to pc.u?

Enable tophat passbands by wavelength and width

Instead of giving filter files as input. This will make spectra runs easier to set up.

Save tli or opacity files as npz?

Test if it is any faster or uses less resources (mem).

isotopes without PF create NaNs

When including a molecule with isotopes without PF, the code (apparently) still includes them, then the spectra will produce NaNs, because the line strengths divide by the PF, which is zero.

E.g., HITRAN N2, isotope '45'.

The fix: do not include in the TLI file these 'bad' line transitions.