Changeset 3671 for trunk/UTIL

Timestamp:

Mar 5, 2025, 5:54:17 PM (4 months ago)

Author:

afalco

Message:

corrk: updated scripts to match API of dace and exo_k.
Added missing atoms in exomol species.
AF

Location:

trunk/UTIL/corrk_exo_k

Files:

• Exomol_species.dat (modified) (1 diff)
• dace_to_corrk.py (modified) (3 diffs)
• download_dace.py (modified) (3 diffs)
• exomol2.py (added)
• lmdz_corrk_using_exo_k.py (modified) (3 diffs)

trunk/UTIL/corrk_exo_k/Exomol_species.dat ¶

@@ +1 @@
+K/               39K                      39K__Kurucz-Allard-K             K//39K/Kurucz-Allard-K                          0        1        0
+Na/              23Na                     23Na__Kurucz-Allard-Na           Na//23Na/Kurucz-Allard-Na                       0        1        0
 MgH/             24Mg-1H                  24Mg-1H__Yadin                   MgH//24Mg-1H/Yadin                              0        1        0
 MgH/             24Mg-1H                  24Mg-1H__UGAMOP                  MgH//24Mg-1H/UGAMOP                             0        1        0

trunk/UTIL/corrk_exo_k/dace_to_corrk.py ¶

@@ -3 +3 @@
 import matplotlib.pyplot as plt
 import astropy.units as u
+import astropy.constants as cst
 import glob
 from chemistry import isotopologue_to_species, species_name_to_common_isotopologue_name
 
 ### choose spectral resolution for corrk
-R=20000
+R=10000
 wnedges=xk.wavenumber_grid_R(10000/100, 10000/.3, R)
 print("Spectral resolution asked: wn=", wnedges.min(), wnedges.max(), ", R =", R)
 
 ### set data paths
-dace_path = "datadir/dace/"
+# dace_path = "datadir/cold/"
+dace_path = "datadir/hot/"
 corrk_dir = "datadir/corrk_data/"
 dir_out = corrk_dir+f"/R{R}_from_dace/"
@@ -26 +28 @@
 
 ### select (output) molecules and their corresponding (input) isotopologues
+### set value to 'None' to get common isotopologue
 molecules = {
-    # "CO":  "12C-16O",
-    # "CO2":  "12C-16O2",
+    "CO":  "12C-16O",
+    "CO2":  None,
     # "H2O":  "1H2-16O",
     # # "C2H2":"12C2-1H2",
     # # "C2H4":"12C2-1H4",
-    "CH4": "12C-1H4",
+    # "CH4": "12C-1H4",
     # "TiO": "48Ti-16O",
     # "VO": "51V-16O",
@@ -61 +64 @@
     tmp_ktab=xk.hires_to_ktable(path=mol_corrk, helios=True, wnedges=wnedges, mol=molecule, order=order) # test conversion from a whole folder to tmp
 
-    tmp_ktab.kdata = tmp_ktab.kdata * 10 * xk.Molar_mass().fetch(molecule) / xk.N_A
+    tmp_ktab.kdata = tmp_ktab.kdata * 10 * xk.Molar_mass().fetch(molecule) / cst.N_A
     # cm^2/g * 10000/1000 * kg/mol / (molecule/mol) -> m^2/molecule
     tmp_ktab.kdata_unit = "m^2/molecule"

trunk/UTIL/corrk_exo_k/download_dace.py ¶

@@ -1 +1 @@
 from dace import *
 
+# output_dir = "datadir/hot"
+# temperature_range=[50,3000]
+# pressure_range=[-8,5]
+
+output_dir = "datadir/cold"
 temperature_range=[50,300]
-pressure_range=[1e-8,1e2]
+pressure_range=[-8,2]
 
 ### select linelist for each molecule
@@ -9 +14 @@
 
 molecules = {
-# "C2H2": default_linelist,
-# "C2H4": default_linelist,
+"C2H2": default_linelist,
+"C2H4": default_linelist,
 "CH4":  default_linelist,
+"CO":   default_linelist,
 # "CO2":  default_linelist,
-# "CO":   default_linelist,
 # "H2O":  default_linelist,
 # "H2S":  default_linelist,
@@ -33 +38 @@
 }
 for mol in molecules:
-    download_molecule(molecule=mol,temperature_range=temperature_range,pressure_range=pressure_range, line_list=molecules[mol], output_dir="datadir/dace")
+    download_molecule(molecule=mol,temperature_range=temperature_range,pressure_range=pressure_range, line_list=molecules[mol], output_dir=output_dir)
 
 

trunk/UTIL/corrk_exo_k/lmdz_corrk_using_exo_k.py ¶

@@ -1 +1 @@
 import exo_k as xk
 import numpy as np
+import matplotlib.pyplot as plt
 
 ### folder with corrk data for each molecule (see 'xsec_to_corrk.py' to convert cross-sections to correlated-k tables)
@@ -42 +43 @@
 
 VI_grid = np.array([ 1670.,  2080.,  2220.,  2450.,  2750.,  2920.,  3200.,  3750., 4000.,  4180.,  4260.,  4460.,  4700.,  5000.,  5400.,  5730., 5930.,  6300.,  6700.,  7700.,  8200.,  8750.,  9250., 35000.])
+res=f"{len(IR_grid)-1}x{len(VI_grid)-1}"
+print(res)
 
 wnedges = np.unique(np.concatenate([IR_grid,VI_grid]))
@@ -62 +65 @@
 xk.Settings().set_case_sensitive(True)
 
-kdata = xk.Kdatabase(gases)
-for gas in kdata.ktables:
-    # convert Xtables to Ktables
-    if type(kdata.ktables[gas]) is xk.Xtable:
-        kdata.ktables[gas] = xk.Ktable(wnedges=wnedges, xtable=kdata.ktables[gas])
+### select pressure, temperature, and mixing ratio to plot
+p=1
+t=300
+x=1e-4
+
+fig, ax = plt.subplots()
+kdata = xk.Kdatabase(molecules=None)
+ktables = []
+for gas in gases:
+    try:
+        ktable = xk.Ktable(mol=gas)
+    except:
+        ktable = xk.Ktable(wnedges=wnedges, xtable=xk.Xtable(mol=gas))
+    ktable.plot_spectrum(ax=ax,p=p,t=t,x=x,g=1,label=gas)
+    ktable.bin_down(wnedges)
+    ktables.append(ktable)
+kdata.add_ktables(*ktables)
 
 mix_var_gas=kdata.create_mix_ktable5d(vgas_comp=vgas,bg_comp=bg_gas,x_array=x_array)