Some worked out use-cases for exo_k

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

Whereas the tutorial goes through all the concepts present in the library in a very progressive manner, here we propose some more complex examples that combine many of these concepts. The idea is to show some ‘real life’ examples of the use of the library that can be adapted by the users to their own needs.

Some of these examples use publicly available data: do not forget to acknowledge them if you use them in your work (see the where to find data section).

General initialization

[23]:
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
import time,sys,os
[24]:
# Uncomment the line below if you want to enable interactive plots
#%matplotlib notebook
plt.rcParams["figure.figsize"] = (7,4)
from matplotlib import cycler
font = {'color': 'dimgray', 'weight': 'bold', 'size': 10}
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)

Setting up global options for the library

[25]:
import exo_k as xk

datapath='data/'
#change to your own path

xk.Settings().set_search_path(datapath+'corrk', path_type='ktable')
xk.Settings().set_search_path(datapath+'xsec', path_type='xtable')
xk.Settings().set_search_path(datapath+'cia', path_type='cia')

xk.Settings().set_mks(True)
# to automatically convert to mks units

Loading ExoMol files and changing their resolution before saving them in different formats

For this, you will need to download the relevant x-sec of ktable files from the ExoMol website (http://exomol.com/data/data-types/opacity/) into the datapath/xsec or datapath/corrk.

[7]:
#let's decide on the spectral grid we want:
wn0=200.
wn1=10000.
Resolution=20.
new_wn_grid=xk.wavenumber_grid_R(wn0, wn1, Resolution)
# Here it is a grid a constant resolution, but any numpy array will do.

# a more focused logP-T grid
logp_array=np.arange(8)
t_array=np.arange(100,600,100)

dir_out='data/corrk/'

molecules = ['1H2-16O','12C-16O2','12C-1H4']
for mol in molecules:
    tmp_ktab=xk.Ktable(mol,'R1000_0.3-50mu.ktable.petitRADTRANS.h5', remove_zeros=True)

    # if for some reasons you need to limit the P-T grid (optional)
    tmp_ktab.remap_logPT(logp_array=logp_array, t_array=t_array)

    # the spectral binning phase
    tmp_ktab.bin_down(new_wn_grid)
    print(tmp_ktab)
    # choose any of the lines below for different formats
    tmp_ktab.write_hdf5(dir_out+mol+'R20.ktable') # hdf5 file with current units
    #tmp_ktab.write_hdf5(dir_out+mol+'R20.ktable', exomol_units=True) # hdf5 file with Exomol units
    #tmp_ktab.write_nemesis(dir_out+mol+'R20.ktable') # binary nemesis format
    #tmp_ktab.write_arcis(dir_out+mol+'R20.ktable') # fits ARCIS format

        file         : /Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/corrk/1H2-16O__POKAZATEL__R1000_0.3-50mu.ktable.petitRADTRANS.h5
        molecule     : H2O
        p grid       : [1.e+00 1.e+01 1.e+02 1.e+03 1.e+04 1.e+05 1.e+06 1.e+07]
        p unit       : Pa
        t grid   (K) : [100 200 300 400 500]
        wn grid      : [ 205.12710964  215.64420145  226.70051608  238.32370009  250.54281748
  263.38842243  276.89263562  291.08922462  306.01368831  321.70334562
  338.19742886  355.53718183  373.76596294  392.92935365  413.07527241
  434.25409451  456.51877804  479.92499631  504.53127705  530.39914878
  557.59329465  586.1817142   616.23589336  647.83098324  681.04598802
  715.96396251  752.67221983  791.26254975  831.8314482   874.48035855
  919.31592529  966.4502607  1016.00122516 1068.09272189 1122.85500677
 1180.42501404 1240.9466987  1304.57139649 1371.45820229 1441.77436795
 1515.69572052 1593.40710189 1675.10283097 1760.98718966 1851.27493358
 1946.19182912 2045.97521795 2150.87461054 2261.15230999 2377.08406799
 2498.95977434 2627.08418177 2761.77766804 2903.37703702 3052.23636091
 3208.72786553 3373.24286117 3546.192721   3728.00990977 3919.14906514
 4120.08813457 4331.3295704  4553.40158624 4786.85947781 5032.28701143
 5290.29788379 5561.53725644 5846.68336912 6146.44923561 6461.58442674
 6792.87694463 7141.15519313 7507.29004927 7892.19704091 8296.83863601
 8722.22664934 9169.42477249 9639.55123371 9940.24491055]
        wn unit      : cm^-1
        kdata unit   : m^2/molecule
        weights      : [0.04555284 0.10007147 0.14116799 0.1632077  0.1632077  0.14116799
 0.10007147 0.04555284 0.00506143 0.01111905 0.01568533 0.01813419
 0.01813419 0.01568533 0.01111905 0.00506143]
        data oredered following p, t, wn, g
        shape        : [ 8  5 79 16]


        file         : /Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/corrk/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5
        molecule     : CO2
        p grid       : [1.e+00 1.e+01 1.e+02 1.e+03 1.e+04 1.e+05 1.e+06 1.e+07]
        p unit       : Pa
        t grid   (K) : [100 200 300 400 500]
        wn grid      : [ 205.12710964  215.64420145  226.70051608  238.32370009  250.54281748
  263.38842243  276.89263562  291.08922462  306.01368831  321.70334562
  338.19742886  355.53718183  373.76596294  392.92935365  413.07527241
  434.25409451  456.51877804  479.92499631  504.53127705  530.39914878
  557.59329465  586.1817142   616.23589336  647.83098324  681.04598802
  715.96396251  752.67221983  791.26254975  831.8314482   874.48035855
  919.31592529  966.4502607  1016.00122516 1068.09272189 1122.85500677
 1180.42501404 1240.9466987  1304.57139649 1371.45820229 1441.77436795
 1515.69572052 1593.40710189 1675.10283097 1760.98718966 1851.27493358
 1946.19182912 2045.97521795 2150.87461054 2261.15230999 2377.08406799
 2498.95977434 2627.08418177 2761.77766804 2903.37703702 3052.23636091
 3208.72786553 3373.24286117 3546.192721   3728.00990977 3919.14906514
 4120.08813457 4331.3295704  4553.40158624 4786.85947781 5032.28701143
 5290.29788379 5561.53725644 5846.68336912 6146.44923561 6461.58442674
 6792.87694463 7141.15519313 7507.29004927 7892.19704091 8296.83863601
 8722.22664934 9169.42477249 9639.55123371 9940.24491055]
        wn unit      : cm^-1
        kdata unit   : m^2/molecule
        weights      : [0.04555284 0.10007147 0.14116799 0.1632077  0.1632077  0.14116799
 0.10007147 0.04555284 0.00506143 0.01111905 0.01568533 0.01813419
 0.01813419 0.01568533 0.01111905 0.00506143]
        data oredered following p, t, wn, g
        shape        : [ 8  5 79 16]


        file         : /Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/corrk/12C-1H4__YT34to10.R1000_0.3-50mu.ktable.petitRADTRANS.h5
        molecule     : CH4
        p grid       : [1.e+00 1.e+01 1.e+02 1.e+03 1.e+04 1.e+05 1.e+06 1.e+07]
        p unit       : Pa
        t grid   (K) : [100 200 300 400 500]
        wn grid      : [ 205.12710964  215.64420145  226.70051608  238.32370009  250.54281748
  263.38842243  276.89263562  291.08922462  306.01368831  321.70334562
  338.19742886  355.53718183  373.76596294  392.92935365  413.07527241
  434.25409451  456.51877804  479.92499631  504.53127705  530.39914878
  557.59329465  586.1817142   616.23589336  647.83098324  681.04598802
  715.96396251  752.67221983  791.26254975  831.8314482   874.48035855
  919.31592529  966.4502607  1016.00122516 1068.09272189 1122.85500677
 1180.42501404 1240.9466987  1304.57139649 1371.45820229 1441.77436795
 1515.69572052 1593.40710189 1675.10283097 1760.98718966 1851.27493358
 1946.19182912 2045.97521795 2150.87461054 2261.15230999 2377.08406799
 2498.95977434 2627.08418177 2761.77766804 2903.37703702 3052.23636091
 3208.72786553 3373.24286117 3546.192721   3728.00990977 3919.14906514
 4120.08813457 4331.3295704  4553.40158624 4786.85947781 5032.28701143
 5290.29788379 5561.53725644 5846.68336912 6146.44923561 6461.58442674
 6792.87694463 7141.15519313 7507.29004927 7892.19704091 8296.83863601
 8722.22664934 9169.42477249 9639.55123371 9940.24491055]
        wn unit      : cm^-1
        kdata unit   : m^2/molecule
        weights      : [0.04555284 0.10007147 0.14116799 0.1632077  0.1632077  0.14116799
 0.10007147 0.04555284 0.00506143 0.01111905 0.01568533 0.01813419
 0.01813419 0.01568533 0.01111905 0.00506143]
        data oredered following p, t, wn, g
        shape        : [ 8  5 79 16]

Creating k-coefficients for a new species not in ExoMol from high-resolution spectra from the petitRADTRANS database

Start by downloading some high resolution cross sections from petitRADTRANS: https://www.dropbox.com/sh/w7sa20v8qp19b4d/AABKF0GsjghsYLJMUJXDgrHma?dl=0

Here I am using data for Na. For the example to work with as few data as possible, I only consider a grid with 2 P and T points. So you will only need the following files in the Na directory:

  • ‘sigma_94_1215.K_0.001000bar.dat’

  • ‘sigma_94_1641.K_0.001000bar.dat’

  • ‘sigma_94_1215.K_1.000000bar.dat’

  • ‘sigma_94_1641.K_1.000000bar.dat’

  • wlen.dat

Feel free to extend to the full grid once you have downloaded the whole dataset.

[8]:
path_to_hires_spectra='/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/hires/PetitRADTRANS/Na'

press_grid_str=['0.001000','1.000000']  # notice the use of strings
logp_grid=[np.log10(float(p)) for p in press_grid_str]
t_grid=[1215, 1641]

file_grid=xk.create_fname_grid('sigma_94_{temp}.K_{press}bar.dat', logpgrid=press_grid_str, tgrid=t_grid,
        p_kw='press', t_kw='temp')
print(file_grid)

Hires_spectra=xk.hires_to_xtable(path=path_to_hires_spectra, filename_grid=file_grid, logpgrid=logp_grid, tgrid=t_grid,
                mol='Na', grid_p_unit='bar', binary=True, mass_amu=23.)

Hires_spectra
[['sigma_94_1215.K_0.001000bar.dat' 'sigma_94_1641.K_0.001000bar.dat']
 ['sigma_94_1215.K_1.000000bar.dat' 'sigma_94_1641.K_1.000000bar.dat']]
[8]:

        file         : /Users/jleconte/atmosphere/RadiativeTransfer/exo_k/data/hires/PetitRADTRANS/Na
        molecule     : Na
        p grid       : [   100. 100000.]
        p unit       : Pa
        t grid   (K) : [1215 1641]
        wn grid      : [  357.14265128   357.14300842   357.14336556 ... 33333.23696353
 33333.27029679 33333.30363008]
        wn unit      : cm^-1
        kdata unit   : m^2/molecule
        data oredered following p, t, wn
        shape        : [      2       2 4536178]

[9]:
fig,ax=plt.subplots(1,1,sharey=False, figsize=(8,4))
Hires_spectra.plot_spectrum(ax, p=1.e3, t=1300., xscale='log', yscale='log', label='p=1 mbar')
Hires_spectra.plot_spectrum(ax, p=1.e5, t=1300., label='p=1 bar')
ax.legend(loc='upper right')
[9]:
<matplotlib.legend.Legend at 0x18252d0860>
_images/examples-exo_k_14_1.png
[10]:
# Let's load a conventional ExoMOL k-table to use the same wavenumber and
#    g-grid
tmp_ktab=xk.Ktable('1H2-16O','R1000_0.3-50mu.ktable.petitRADTRANS.h5', remove_zeros=True)

wnedges=tmp_ktab.wnedges
weights=tmp_ktab.weights
ggrid=tmp_ktab.ggrid
## Or we can create a custom g-grid with 8 gauss legendre points between 0 and 0.95
##   and 8 points between 0.95 and 1 (as usual with petitRADTRANS data)
#weights, ggrid, gedges = xk.split_gauss_legendre(order=16, g_split=0.95)

ktab=xk.Ktable(xtable=Hires_spectra, wnedges=wnedges, weights=weights, ggrid=ggrid)

## choose any of the lines below for different formats
#full_path_to_write='data/corrk/Na_R1000.ktable'
#tmp_ktab.write_hdf5(full_path_to_write) # hdf5 file with current units
#tmp_ktab.write_hdf5(full_path_to_write, exomol_units=True) # hdf5 file with Exomol units
#tmp_ktab.write_nemesis(full_path_to_write) # binary nemesis format
#tmp_ktab.write_arcis(full_path_to_write) # fits ARCIS format
[11]:
fig,ax=plt.subplots(1,1,sharey=False, figsize=(8,4))
Hires_spectra.plot_spectrum(ax, p=1.e3, t=1300., xscale='log', yscale='log', label='p=1 mbar')
Hires_spectra.plot_spectrum(ax, p=1.e5, t=1300., label='p=1 bar')
ktab.plot_spectrum(ax, p=1.e5, t=1300., g=1, label='ktab, g=1, R=1000, p=1 bar')
ax.legend(loc='upper right')
/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/exo_k/data_table.py:235: RuntimeWarning: divide by zero encountered in log
  kc_p1t1=np.log(self.kdata[lpind[ii],tind[ii]][wngrid_filter].ravel())
/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/exo_k/data_table.py:236: RuntimeWarning: divide by zero encountered in log
  kc_p0t1=np.log(self.kdata[lpind[ii]-1,tind[ii]][wngrid_filter].ravel())
/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/exo_k/data_table.py:237: RuntimeWarning: divide by zero encountered in log
  kc_p1t0=np.log(self.kdata[lpind[ii],tind[ii]-1][wngrid_filter].ravel())
/Users/jleconte/atmosphere/RadiativeTransfer/exo_k/exo_k/data_table.py:238: RuntimeWarning: divide by zero encountered in log
  kc_p0t0=np.log(self.kdata[lpind[ii]-1,tind[ii]-1][wngrid_filter].ravel())
[11]:
<matplotlib.legend.Legend at 0x1831336550>
_images/examples-exo_k_16_2.png

Modelling transit spectra: sampled cross sections vs. k-coefficients (Leconte, A&A, 2020)

Here we will reproduce some of the results shown in Leconte (A&A, 2020) that demonstrate that the correlated-k method can be much more accurate than the sampled cross-section technique in many cases of interest. Although we use a different set of initial data (because these are publicly available and easily accessible), the results are mostly unaffected.

Start by downloading some high resolution cross sections for water from petitRADTRANS: https://www.dropbox.com/sh/w7sa20v8qp19b4d/AABKF0GsjghsYLJMUJXDgrHma?dl=0

[12]:
path_to_hires_spectra='data/hires/PetitRADTRANS/H2O/'

press_grid_str=['0.000001','0.000010','0.000100','0.001000','0.010000','0.100000','1.000000','10.000000','100.000000']
logp_grid=[np.log10(float(p)) for p in press_grid_str]
t_grid=[900, 1215]

file_grid=xk.create_fname_grid('sigma_01_{temp}.K_{press}bar.dat', logpgrid=press_grid_str, tgrid=t_grid,
        p_kw='press', t_kw='temp')
print(file_grid)

h2o_hires=xk.hires_to_xtable(path=path_to_hires_spectra, filename_grid=file_grid, logpgrid=logp_grid, tgrid=t_grid,
                mol='H2O', grid_p_unit='bar', binary=True, mass_amu=18.)

h2o_hires
[['sigma_01_900.K_0.000001bar.dat' 'sigma_01_1215.K_0.000001bar.dat']
 ['sigma_01_900.K_0.000010bar.dat' 'sigma_01_1215.K_0.000010bar.dat']
 ['sigma_01_900.K_0.000100bar.dat' 'sigma_01_1215.K_0.000100bar.dat']
 ['sigma_01_900.K_0.001000bar.dat' 'sigma_01_1215.K_0.001000bar.dat']
 ['sigma_01_900.K_0.010000bar.dat' 'sigma_01_1215.K_0.010000bar.dat']
 ['sigma_01_900.K_0.100000bar.dat' 'sigma_01_1215.K_0.100000bar.dat']
 ['sigma_01_900.K_1.000000bar.dat' 'sigma_01_1215.K_1.000000bar.dat']
 ['sigma_01_900.K_10.000000bar.dat' 'sigma_01_1215.K_10.000000bar.dat']
 ['sigma_01_900.K_100.000000bar.dat' 'sigma_01_1215.K_100.000000bar.dat']]
[12]:

        file         : data/hires/PetitRADTRANS/H2O/
        molecule     : H2O
        p grid       : [1.e-01 1.e+00 1.e+01 1.e+02 1.e+03 1.e+04 1.e+05 1.e+06 1.e+07]
        p unit       : Pa
        t grid   (K) : [ 900 1215]
        wn grid      : [  357.14265128   357.14300842   357.14336556 ... 33333.23696353
 33333.27029679 33333.30363008]
        wn unit      : cm^-1
        kdata unit   : m^2/molecule
        data oredered following p, t, wn
        shape        : [      9       2 4536178]

Notice the pressure grid has been automatically converted to Pa

[13]:
fig,ax=plt.subplots(1,1,sharey=False, figsize=(8,4))
h2o_hires.plot_spectrum(ax, p=1., t=1000., yscale='log', label='p=1 Pa')
h2o_hires.plot_spectrum(ax, p=1.e5, t=1000., yscale='log', label='p=1 bar')
ax.legend(loc='lower right')
[13]:
<matplotlib.legend.Legend at 0x1831502dd8>
_images/examples-exo_k_22_1.png

We clip the data to focus on a smaller wavenumber range (roughly, the WFC3 region).

[14]:
wn0=5000.;wn1=10000.;
h2o_hires.clip_spectral_range(wn_range=[wn0-2.,wn1+2.])
[15]:
Ri=15000.
wn_interm_grid=xk.wavenumber_grid_R(wn0, wn1, Ri)
sampled=h2o_hires.sample_cp(wn_interm_grid)

fig,axs=plt.subplots(1,2,sharey=False, figsize=(9,4))
h2o_hires.plot_spectrum(axs[0], p=1., t=1000., yscale='log', label='HR')
sampled.plot_spectrum(axs[0], p=1., t=1000., yscale='log', label='R='+str(Ri))
h2o_hires.plot_spectrum(axs[1], p=1., t=1000., yscale='log', label='HR')
sampled.plot_spectrum(axs[1], p=1., t=1000., yscale='log', label='R='+str(Ri))
axs[0].legend(loc='upper left')
axs[1].set_xlim(1.4,1.402)
fig.tight_layout()
_images/examples-exo_k_25_0.png

I can directly write this to an hdf5 file that is compatible with TauREX and petitRADTRANS (see above)

Now let’s compute a reference, high-resolution transmission spectrum:

[16]:
start=time.time()

database=xk.Kdatabase(None)
database.add_ktables(h2o_hires)

T0=1000.; xH2O=1.e-3; Rp=1*u.Rjup; Rs=1.*u.Rsun; grav=10.; nlev=100;

atm=xk.Atm(psurf=10.e5, ptop=1.e-4, Tsurf=T0, grav=grav,
                    composition={'H2':'background','H2O':xH2O}, Nlev=nlev, Rp=Rp,
                    k_database=database)
spec_ref=atm.transmission_spectrum(Rstar=Rs)

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

computation time= 15.350582122802734 s
[17]:
fig,ax=plt.subplots(1,1,sharex=True,figsize=(7,3))
spec_ref.plot_spectrum(ax)
ax.set_ylabel('Transit depth', fontdict=font)
fig.tight_layout()
_images/examples-exo_k_29_0.png

In the following, we sample our cross sections on a grid of intermediate resolution (Rinter) to compute the spectrum and then bin down to the desired final resolution

[18]:
timeit=False
#grid of final (instrumental data) resolution
Rfin=[10,20,50,100.,200.,500.,1000.]
#grid of intermediate resolution (at which computation will be performed)
Rinter=[1000.,3000.,10000.,30000.]

start=time.time()
res=[]
for Ri in Rinter:
    wn_interm_grid=xk.wavenumber_grid_R(wn0, wn1, Ri)
    database_LR=database.copy()
    database_LR.sample(wn_interm_grid)
    atm.set_k_database(k_database=database_LR)
    spec_LR=atm.transmission_spectrum(Rstar=Rs)
    if timeit:
        %timeit atm.transmission_spectrum(Rstar=Rs)
    for Rf in Rfin:
        print('intermediate resolution=',Ri,', final resolution=',Rf)
        wn_final_grid=xk.wavenumber_grid_R(wn0, wn1, Rf)
        spectmp=spec_LR.bin_down_cp(wn_final_grid)
        spec_comp=spec_ref.bin_down_cp(wn_final_grid)
        res.append([Ri,Rf,(spectmp-spec_comp).std(),spectmp,spectmp-spec_comp])

print('computation time=',time.time()-start,'s')
res2=np.array(res).reshape((len(Rinter),len(Rfin),5))

intermediate resolution= 1000.0 , final resolution= 10
intermediate resolution= 1000.0 , final resolution= 20
intermediate resolution= 1000.0 , final resolution= 50
intermediate resolution= 1000.0 , final resolution= 100.0
intermediate resolution= 1000.0 , final resolution= 200.0
intermediate resolution= 1000.0 , final resolution= 500.0
intermediate resolution= 1000.0 , final resolution= 1000.0
intermediate resolution= 3000.0 , final resolution= 10
intermediate resolution= 3000.0 , final resolution= 20
intermediate resolution= 3000.0 , final resolution= 50
intermediate resolution= 3000.0 , final resolution= 100.0
intermediate resolution= 3000.0 , final resolution= 200.0
intermediate resolution= 3000.0 , final resolution= 500.0
intermediate resolution= 3000.0 , final resolution= 1000.0
intermediate resolution= 10000.0 , final resolution= 10
intermediate resolution= 10000.0 , final resolution= 20
intermediate resolution= 10000.0 , final resolution= 50
intermediate resolution= 10000.0 , final resolution= 100.0
intermediate resolution= 10000.0 , final resolution= 200.0
intermediate resolution= 10000.0 , final resolution= 500.0
intermediate resolution= 10000.0 , final resolution= 1000.0
intermediate resolution= 30000.0 , final resolution= 10
intermediate resolution= 30000.0 , final resolution= 20
intermediate resolution= 30000.0 , final resolution= 50
intermediate resolution= 30000.0 , final resolution= 100.0
intermediate resolution= 30000.0 , final resolution= 200.0
intermediate resolution= 30000.0 , final resolution= 500.0
intermediate resolution= 30000.0 , final resolution= 1000.0
computation time= 1.5046379566192627 s

Here are the spectra for each intermediate and final resolution. The right column shows the difference with the reference, high-resolution spectrum.

[19]:
to_plot=res2[:,[0,2,3,6]]
nRfin=to_plot.shape[1]
fig,axs=plt.subplots(nRfin,2,sharex=True,figsize=(7,8))
color_idx = np.linspace(0.1, 0.9, to_plot.shape[0])
offset=11000.
for iRfin in range(nRfin):
    for i, data in enumerate(to_plot[:,iRfin]):
        (data[3]*1.e6-offset).plot_spectrum(axs[iRfin,0],label='$R_{sp}$='+str(int(data[0]/1.e3))+'k',
                    color=plt.cm.inferno(color_idx[i]))
        (data[4].abs()*1.e6).plot_spectrum(axs[iRfin,1],yscale='log',
                    label='$R_{sp}$='+str(int(data[0]/1.e3))+'k',
                    color=plt.cm.inferno(color_idx[i]))
    wn_final_grid=xk.wavenumber_grid_R(wn0,wn1,to_plot[0,iRfin,1])
    spec_comp=spec_ref.bin_down_cp(wn_final_grid)
    (spec_comp*1.e6-offset).plot_spectrum(axs[iRfin,0],marker='.', label='Reference')
    if iRfin==0:
        (spec_ref*1.e6-offset).plot_spectrum(axs[0,0],alpha=0.3)
        axs[iRfin,0].legend(fontsize='x-small', loc='upper right')
        axs[iRfin,1].legend(fontsize='x-small')
    axs[iRfin,0].text(0.01, .9, 'Final resolution: '+str(int(to_plot[0,iRfin,1])), transform=axs[iRfin,0].transAxes)
    axs[iRfin,0].set_ylabel('Depth (ppm)', fontdict=font)
    axs[iRfin,1].set_ylabel('Difference (ppm)', fontdict=font)
    axs[iRfin,0].set_xlabel(None)
    axs[iRfin,1].set_xlabel(None)
axs[-1,0].set_xlabel('Wavelength ($\mu m$)', fontdict=font)
axs[-1,1].set_xlabel('Wavelength ($\mu m$)', fontdict=font)
fig.tight_layout()
_images/examples-exo_k_33_0.png

Let’s summarize these results by looking only at the RMS error as a function of the final resolution for various intermediate resolutions.

[20]:
nRint=res2.shape[0]
color_idx = np.linspace(0.2, 0.9, nRint)
fig,ax=plt.subplots(1,1, sharex=True, figsize= (4.5,3))
for iRint in range(nRint):
    ax.plot(res2[iRint,:,1], res2[iRint,:,2]*1.e6, label='$R_{sp}$='+str(int(res2[iRint,0,0]/1.e3))+'k',
            color=plt.cm.inferno(color_idx[iRint]), linestyle='--')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel('RMS Error (ppm)', fontdict=font)
ax.set_xlabel('Final Resolution ($R_{fin}$)', fontdict=font)
ax.legend(fontsize='x-small', loc="upper left", ncol=2, frameon=False)
fig.tight_layout()
_images/examples-exo_k_35_0.png

Finaly, let’s compare with the result given by the correlated-k approach, computed directly at the final resolution.

[21]:
rescorrk=[]
order=8
for Rf in Rfin:
    print('Resolution=', Rf)
    wn_final_grid=xk.wavenumber_grid_R(wn0, wn1, Rf)
    database_corrk=xk.Kdatabase(None)
    database_corrk.add_ktables(xk.Ktable(xtable=h2o_hires, wnedges=wn_final_grid, order=order))
    atm.set_k_database(k_database=database_corrk)
    spec_corrk=atm.transmission_spectrum(Rstar=Rs)
    if timeit:
        %timeit atm.transmission_spectrum(Rstar=Rs)
    spectmp=spec_corrk.copy()
    spec_comp=spec_ref.bin_down_cp(wn_final_grid)
    rescorrk.append([1.,Rf,(spectmp-spec_comp).std(),spectmp,spectmp-spec_comp])

rescorrk=np.array(rescorrk).reshape((len(Rfin),5))
Resolution= 10
Resolution= 20
Resolution= 50
Resolution= 100.0
Resolution= 200.0
Resolution= 500.0
Resolution= 1000.0
[22]:
nRint=res2.shape[0]
color_idx = np.linspace(0.2, 0.9, nRint)
fig,ax=plt.subplots(1,1, sharex=True, figsize= (4.5,3))
for iRint in range(nRint):
    ax.plot(res2[iRint,:,1], res2[iRint,:,2]*1.e6, label='$R_{sp}$='+str(int(res2[iRint,0,0]/1.e3))+'k',
            color=plt.cm.inferno(color_idx[iRint]), linestyle='--')
ax.plot(rescorrk[:,1], rescorrk[:,2]*1.e6, marker='.', color='k', label='$k$-coeff')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel('RMS Error (ppm)', fontdict=font)
ax.set_xlabel('Final Resolution ($R_{fin}$)', fontdict=font)
ax.legend(fontsize='x-small', loc="upper left", ncol=2, frameon=False)
fig.tight_layout()
_images/examples-exo_k_38_0.png