Changeset 3633 for trunk

Timestamp:

Feb 19, 2025, 2:27:45 PM (4 months ago)

Author:

afalco

Message:

exo_k: updated scripts for better info and removed useless binning.
AF

Location:

trunk/UTIL/corrk_exo_k

Files:

• bin_corrk.py (modified) (2 diffs)
• dace_to_corrk.py (modified) (4 diffs)
• lmdz_corrk_using_exo_k.py (moved) (moved from trunk/UTIL/corrk_exo_k/lmdz_corrk_from_exomol.py) (1 diff)
• xsec_to_corrk.py (modified) (2 diffs)

trunk/UTIL/corrk_exo_k/bin_corrk.py

@@ -7 +7 @@
 R=100
 wnedges=xk.wavenumber_grid_R(10000/100, 10000/.3, R)
-print("wnedges range: ", wnedges.min(), wnedges.max())
+print("Spectral resolution asked: wn=", wnedges.min(), wnedges.max(), ", R =", R)
 
 ### folder where corrk will be saved
@@ -66 +66 @@
 
         tmp_ktab=xk.Ktable(mol=mol, remove_zeros=True)
+        print(f"Initial spectral resolution = {tmp_ktab.Nw}")
 
-        if logpgrid is None: # use first molecule to set p and t grid
+        if logpgrid is None or tgrid is None: # use first molecule to set p and t grid
             logpgrid = tmp_ktab.logpgrid
             tgrid = tmp_ktab.tgrid
 
-        print(f"Remapping: {mol}")
-        tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid)
         print(f"Binning: {mol}")
         tmp_ktab.bin_down(wnedges)
+
+        if (tmp_ktab.logpgrid != logpgrid).any() or (tmp_ktab.tgrid != tgrid).any():
+            print("Initial P grid:", tmp_ktab.pgrid)
+            print("Initial T grid:", tmp_ktab.tgrid)
+            print(f"Remapping {mol}")
+            print("New P grid:", np.power(logpgrid,10))
+            print("New T grid:", tgrid)
+            tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid)
+
+        print(tmp_ktab)
         print(f"Writing: {mol}")
         tmp_ktab.write_hdf5(dir_out+mol+f"_R{R}.corrk.h5")

trunk/UTIL/corrk_exo_k/dace_to_corrk.py

@@ -7 +7 @@
 
 ### choose spectral resolution for corrk
-R=500
+R=20000
 wnedges=xk.wavenumber_grid_R(10000/100, 10000/.3, R)
-print(wnedges.min(), wnedges.max())
+print("Spectral resolution asked: wn=", wnedges.min(), wnedges.max(), ", R =", R)
 
 ### set data paths
@@ -15 +15 @@
 corrk_dir = "datadir/corrk_data/"
 dir_out = corrk_dir+f"/R{R}_from_dace/"
-ref=xk.Ktable(filename=corrk_dir + '/R500_from_R15000/CH4_R500.corrk.h5', mol="CH4")
+
+### you can set your desired p and t grid here (otherwise, first molecule is used)
+logpgrid = None
+tgrid = None
+# ref=xk.Ktable(filename=corrk_dir + '/R500_from_R15000/CH4_R500.corrk.h5', mol="CH4")
+# logpgrid = ref.logpgrid
+# tgrid = ref.tgrid
+
+
 
 ### select (output) molecules and their corresponding (input) isotopologues
 molecules = {
-    "CO":  "12C-16O",
-    "C2H2":"12C2-1H2",
-    "C2H4":"12C2-1H4",
+    # "CO":  "12C-16O",
+    # "CO2":  "12C-16O2",
+    # "H2O":  "1H2-16O",
+    # # "C2H2":"12C2-1H2",
+    # # "C2H4":"12C2-1H4",
     "CH4": "12C-1H4",
+    # "TiO": "48Ti-16O",
+    # "VO": "51V-16O",
 }
 
-xk.Settings().set_mks(True)
+# xk.Settings().set_mks(True)
 xk.Settings().set_log_interp(False)
 xk.Settings().set_case_sensitive(True)
 
-def dace_to_corrk(molecule=None, isotopologue=None, order=17):
+def dace_to_corrk(molecule=None, isotopologue=None, order=17, logpgrid=logpgrid, tgrid=tgrid):
     if molecule is None:
         molecule=isotopologue_to_species(isotopologue)
@@ -43 +55 @@
         return
 
-    tmp=xk.hires_to_ktable(path=mol_corrk, helios=True, wnedges=wnedges, mol=molecule, order=order) # test conversion from a whole folder to tmp
+    print(f"Reading: {mol}")
+    hr = xk.Hires_spectrum(glob.glob(mol_corrk+"/*.bin")[0], helios=True)
+    print(f"Initial spectral resolution = {hr.Nw}")
 
-    tmp.kdata = tmp.kdata * 10 * xk.Molar_mass().fetch(molecule) / xk.N_A
+    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
     # cm^2/g * 10000/1000 * kg/mol / (molecule/mol) -> m^2/molecule
-    tmp.kdata_unit = "m^2/molecule"
-    tmp.remap_logPT(ref.logpgrid, ref.tgrid)
-    tmp.remove_zeros()
+    tmp_ktab.kdata_unit = "m^2/molecule"
+    tmp_ktab.remove_zeros()
 
-    print(tmp)
-    tmp.write_hdf5(dir_out+f"{molecule}_corrk_dace_{R}.h5")
+
+    if logpgrid is None or tgrid is None: # use first molecule to set p and t grid
+            logpgrid = tmp_ktab.logpgrid
+            tgrid = tmp_ktab.tgrid
+
+    if (tmp_ktab.logpgrid != logpgrid).any() or (tmp_ktab.tgrid != tgrid).any():
+        print("Initial P grid:", tmp_ktab.pgrid)
+        print("Initial T grid:", tmp_ktab.tgrid)
+        print(f"Remapping {mol}")
+        print("New P grid:", np.power(logpgrid,10))
+        print("New T grid:", tgrid)
+        tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid)
+
+    print(tmp_ktab)
+    print(f"Writing: {mol}")
+    tmp_ktab.write_hdf5(dir_out+f"{molecule}_corrk_dace_{R}.h5")
 
     # %matplotlib notebook
@@ -58 +87 @@
     t_plot=300
     fig,ax=plt.subplots(figsize=(9,6.5))
-    tmp.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='g=1')
+    tmp_ktab.plot_spectrum(ax,p=p_plot,t=t_plot,g=1.,yscale='log',xscale='log',label='g=1')
     fig.tight_layout()
     plt.savefig(dir_out+f"corrk_dace_R{R}_{molecule}_p{p_plot}_t{t_plot}.pdf")

trunk/UTIL/corrk_exo_k/lmdz_corrk_using_exo_k.py

@@ -4 +4 @@
 ### folder with corrk data for each molecule (see 'xsec_to_corrk.py' to convert cross-sections to correlated-k tables)
 xk.Settings().set_search_path("datadir/exomol/taurex_R15000/")
+# xk.Settings().set_search_path("datadir/corrk_data/R20000_from_dace/")
 
 corrk_dir = "datadir/corrk_data/"

trunk/UTIL/corrk_exo_k/xsec_to_corrk.py

@@ -7 +7 @@
 R=500
 wnedges=xk.wavenumber_grid_R(10000/100, 10000/.3, R)
-print("wnedges range: ", wnedges.min(), wnedges.max())
+print("Spectral resolution asked: wn=", wnedges.min(), wnedges.max(), ", R =", R)
 
 ### folder where corrk will be saved
@@ -63 +63 @@
             print(f"Reading: {mol}")
 
-        tmp_ktab=xk.Ktable(wnedges=wnedges, xtable=xk.Xtable(pattern, mol=mol, remove_zeros=True))
+        xtable=xk.Xtable(pattern, mol=mol, remove_zeros=True)
+        print(f"Initial spectral resolution = {xtable.Nw}")
 
-        if logpgrid is None: # use first molecule to set p and t grid
+        tmp_ktab=xk.Ktable(wnedges=wnedges, xtable=xtable)
+
+        # print(f"Binning: {mol}")
+        # tmp_ktab.bin_down(wnedges)
+
+
+        if logpgrid is None or tgrid is None: # use first molecule to set p and t grid
             logpgrid = tmp_ktab.logpgrid
             tgrid = tmp_ktab.tgrid
 
-        print(f"Remapping: {mol}")
-        tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid)
+        if (tmp_ktab.logpgrid != logpgrid).any() or (tmp_ktab.tgrid != tgrid).any():
+            print("Initial P grid:", tmp_ktab.pgrid)
+            print("Initial T grid:", tmp_ktab.tgrid)
+            print(f"Remapping {mol}")
+            print("New P grid:", np.power(logpgrid,10))
+            print("New T grid:", tgrid)
+            tmp_ktab.remap_logPT(logp_array=logpgrid, t_array=tgrid)
+
+        print(tmp_ktab)
         print(f"Writing: {mol}")
         tmp_ktab.write_hdf5(dir_out+mol+f"_R{R}.corrk.h5")