In [1]:
import numpy as np
# from petitRADTRANS import Radtrans
import matplotlib.pyplot as plt
# from petitRADTRANS import nat_cst as nc
import JW_lib
import spectres
import astropy.io.fits
from astropy.io.votable import parse
import glob
# from poor_mans_nonequ_chem import poor_mans_nonequ_chem as pm
In [2]:
spec = JW_lib.readSpec(file_name="spec_combined_updated.txt")
plt.plot(spec.wavelength, spec.flux, marker="o", markersize=0.1, alpha=0.1)
plt.ylim(0.6, 1.2)
plt.xlim(2.29, 2.325)
Out[2]:
(2.29, 2.325)
In [3]:
from jax import config
config.update("jax_enable_x64", True)
from exojax.utils.grids import wavenumber_grid
from exojax.spec.multimol import MultiMol
nu_grid, wav, resolution = wavenumber_grid(4291.,
4366.,
100000,
unit="cm-1",
xsmode="premodit")
# 初始化 MultiMol 对象，确保使用正确的数据库名称
molmulti = [["H2O", "CO"]]
dbmulti = [["ExoMol", "ExoMol"]] # 假设所有分子都使用 ExoMol 数据库
mul = MultiMol(molmulti=molmulti, dbmulti=dbmulti)
# 获取多分子数据库
multimdb = mul.multimdb([nu_grid], crit=1.e-30, Ttyp=1500.)
# 生成多分子光谱吸收系数
auto_trange = [500., 2500.] # 自动温度范围，适用于大部分大气条件
dit_grid_resolution = 0.1 # DIT 网格分辨率，根据需要调整
multiopa = mul.multiopa_premodit(multimdb, [nu_grid], auto_trange, dit_grid_resolution=dit_grid_resolution)
xsmode = premodit xsmode assumes ESLOG in wavenumber space: mode=premodit ====================================================================== The wavenumber grid should be in ascending order. The users can specify the order of the wavelength grid by themselves. Your wavelength grid is in *** descending *** order ====================================================================== Multiple recommended databases found for H2O in ExoMol : ['POKAZATEL', 'POKAZATEL', 'POKAZATEL', 'POKAZATEL', 'POKAZATEL']. This is unexpected. Using the first Multiple recommended databases found for CO in ExoMol : ['Li2015', 'Li2015', 'Li2015']. This is unexpected. Using the first HITRAN exact name= H2(16O) Molecule: H2O Isotopologue: 1H2-16O Background atmosphere: H2 ExoMol database: None Local folder: .database\H2O\1H2-16O\POKAZATEL Transition files: => File 1H2-16O__POKAZATEL__04200-04300.trans => File 1H2-16O__POKAZATEL__04300-04400.trans # i_upper i_lower A nu_lines gup jlower jupper elower eupper Sij0 0 174444 167140 0.0026972 4200.000001 27 13 13 35369.749297 39569.749298 6.785075061528842e-98 1 69643 73288 0.4998 4200.000003000001 17 8 8 27489.345489 31689.345492 3.419571737978315e-79 2 570559 574315 0.00052419 4200.000006000002 61 30 30 32356.435604 36556.43561 6.841786739366168e-92 3 153954 158804 1.7935e-05 4199.648955999997 75 12 12 17973.809935 22173.458891 6.61845102635208e-63 4 282122 274199 1.6985e-05 4200.000009000003 105 17 17 35749.231065 39949.231074 2.626942731689608e-100 ... ... ... ... ... ... ... ... ... ... ... 84,947,935 533603 504505 0.0041694 4399.999992000005 57 27 28 31610.705314 36010.705306 1.7382864423065486e-89 84,947,936 703680 692571 0.062086 4399.999993000005 77 37 38 33591.900357 37991.90035 2.2981354451799078e-92 84,947,937 471458 498800 0.03739 4399.999993999998 153 26 25 31459.38811 35859.388104 8.730695316590409e-88 84,947,938 135402 118795 0.08965 4399.999995999999 69 10 11 31344.456996 35744.456992 1.650515787016003e-87 84,947,939 518271 543636 0.09221 4399.999997999999 165 28 27 34205.620029 38605.620027 3.703373088794551e-93
C:\Users\10923\anaconda3\envs\exojax_test1\lib\site-packages\exojax\utils\molname.py:178: FutureWarning: e2s will be replaced to exact_molname_exomol_to_simple_molname. warnings.warn(
Broadening code level: a1
C:\Users\10923\anaconda3\envs\exojax_test1\lib\site-packages\exojax\utils\molname.py:178: FutureWarning: e2s will be replaced to exact_molname_exomol_to_simple_molname. warnings.warn( c:\users\10923\radis\radis\api\exomolapi.py:607: AccuracyWarning: The default broadening parameter (alpha = 0.07 cm^-1 and n = 0.5) are used for J'' > 80 up to J'' = 93 warnings.warn( C:\Users\10923\anaconda3\envs\exojax_test1\lib\site-packages\exojax\spec\opacalc.py:171: UserWarning: dit_grid_resolution is not None. Ignoring broadening_parameter_resolution. warnings.warn(
HITRAN exact name= (12C)(16O) Molecule: CO Isotopologue: 12C-16O Background atmosphere: H2 ExoMol database: None Local folder: .database\CO\12C-16O\Li2015 Transition files: => File 12C-16O__Li2015.trans # i_upper i_lower A nu_lines gup jlower jupper elower Sij0 0 84 42 1.155e-06 2.405586 3 0 1 66960.7124 3.811968891483239e-164 1 83 41 1.161e-06 2.441775 3 0 1 65819.903 9.66302808612315e-162 2 82 40 1.162e-06 2.477774 3 0 1 64654.9206 2.743839242930895e-159 3 81 39 1.159e-06 2.513606 3 0 1 63465.8042 8.733228323835037e-157 4 80 38 1.152e-06 2.549292 3 0 1 62252.5793 3.1152203985525016e-154 ... ... ... ... ... ... ... ... ... ... 125,491 306 253 7.164e-10 22147.135424 15 6 7 80.7354 1.8282485560395954e-31 125,492 474 421 9.852e-10 22147.86595 23 10 11 211.4041 2.0425455628245774e-31 125,493 348 295 7.72e-10 22147.897299 17 7 8 107.6424 1.9589545214604644e-31 125,494 432 379 9.056e-10 22148.262711 21 9 10 172.978 2.0662209079393328e-31 125,495 390 337 8.348e-10 22148.273111 19 8 9 138.3903 2.03878272167021e-31 Broadening code level: a0 OpaPremodit: params automatically set. default elower grid trange (degt) file version: 2 Tref changed: 296.0K->2036.6563614826916K OpaPremodit: Tref_broadening is set to 1118.033988749894 K # of reference width grid : 10 # of temperature exponent grid : 11
uniqidx: 100%|████████████████████████████████████████████████████████████████████████| 20/20 [00:00<00:00, 245.13it/s]
Premodit: Twt= 603.4302162255675 K Tref= 2036.6563614826916 K Making LSD:|####################| 100% Making LSD:|####################| 100% Making LSD:|####################| 100%
C:\Users\10923\anaconda3\envs\exojax_test1\lib\site-packages\exojax\spec\opacalc.py:171: UserWarning: dit_grid_resolution is not None. Ignoring broadening_parameter_resolution. warnings.warn(
OpaPremodit: params automatically set. default elower grid trange (degt) file version: 2 Tref changed: 296.0K->2036.6563614826916K OpaPremodit: Tref_broadening is set to 1118.033988749894 K # of reference width grid : 5 # of temperature exponent grid : 3
uniqidx: 100%|███████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<?, ?it/s]
Premodit: Twt= 603.4302162255675 K Tref= 2036.6563614826916 K Making LSD:|####################| 100% Making LSD:|####################| 100% Making LSD:|####################| 100%
In [4]:
from exojax.spec.atmrt import ArtEmisPure
import numpy as np
art = ArtEmisPure(nu_grid=nu_grid, pressure_btm=1.e2, pressure_top=1.e-5, nlayer=75, rtsolver="ibased", nstream=8)
art.change_temperature_range(500.0, 2500.0)
Tarr = art.powerlaw_temperature(1600.0,0.07)
Parr = art.pressure
plt.plot(Tarr, Parr)
plt.yscale("log")
plt.ylim(np.max(Parr),np.min(Parr))
rtsolver: ibased Intensity-based n-stream solver, isothermal layer (e.g. NEMESIS, pRT like)
C:\Users\10923\anaconda3\envs\exojax_test1\lib\site-packages\exojax\spec\dtau_mmwl.py:14: FutureWarning: dtau_mmwl might be removed in future.
warnings.warn("dtau_mmwl might be removed in future.", FutureWarning)
Out[4]:
(100.0, 1e-05)
In [5]:
# 定义 Mixing Ratio
mmr_CO = art.constant_mmr_profile(0.01)
mmr_H2O = art.constant_mmr_profile(0.01)
mdb_CO = multimdb[0][0]
mdb_H2O = multimdb[0][1]
# multiopa 的第一个元素对应 H2O，第二个元素对应 CO
xsmatrix_H2O = multiopa[0][0].xsmatrix(Tarr, Parr)
xsmatrix_CO = multiopa[0][1].xsmatrix(Tarr, Parr)
from exojax.utils.astrofunc import gravity_jupiter
gravity = gravity_jupiter(1.0,70.0)
# 使用 H2O 和 CO 的 mixing ratio
dtau_H2O = art.opacity_profile_xs(xsmatrix_H2O, mmr_H2O, mdb_H2O.molmass, gravity)
dtau_CO = art.opacity_profile_xs(xsmatrix_CO, mmr_CO, mdb_CO.molmass, gravity)
from exojax.spec.contdb import CdbCIA
from exojax.spec.opacont import OpaCIA
cdb = CdbCIA(".database/H2-H2_2011.cia",nurange=nu_grid)
opacia = OpaCIA(cdb, nu_grid=nu_grid)
# for k in range(len(wavelength_grid_segments)):
# for i in range(len(mul.molmulti[k])):
# xsm = multiopa[k][i].xsmatrix(Tarr, Parr)
# xsmatrix = opa.xsmatrix(Tarr, Parr)
logacia_matrix = opacia.logacia_matrix(Tarr)
vmrH2 = 0.855 #VMR of H2
mmw = 2.33 # mean molecular weight of the atmosphere
dtaucia = art.opacity_profile_cia(logacia_matrix, Tarr, vmrH2, vmrH2, mmw, gravity)
dtau = dtau_CO + dtau_H2O + dtaucia
F = art.run(dtau, Tarr)
F_h2o = art.run(dtau_H2O, Tarr)
F_co = art.run(dtau_CO, Tarr)
H2-H2
In [12]:
f_norm = F/np.max(F)
f_norm_h2o = F_h2o/np.max(F_h2o)
f_norm_co = F_co/np.max(F_co)
wavegrid = 1e4/nu_grid
modeled_spec = JW_lib.Spectrum(wavegrid, f_norm)
window_size = 1500
window = np.ones(window_size) / window_size
modeled_spec.flux = np.convolve(modeled_spec.flux, window, mode='same')
window_size = 5
window = np.ones(window_size) / window_size
spec_copy = spec.copy()
spec_copy.flux = np.convolve(spec_copy.flux, window, mode='same')
fig=plt.figure(figsize=(15,4))
y_offset = 0.18
plt.plot(modeled_spec.wavelength,modeled_spec.flux + y_offset, color="red", alpha=0.5, label="modeled spectrum")
plt.plot(spec_copy.wavelength, spec_copy.flux, marker="o", markersize=0.8, alpha=0.2, color="black", label="observation")
plt.plot(wavegrid, f_norm_h2o, color="blue", alpha=0.5, label="H2O")
plt.plot(wavegrid, f_norm_co, color="orange", alpha=0.5, label="CO")
plt.ylim(0.0, 1.2)
plt.xlim(2.29, 2.325)
# plt.xlim(2.295, 2.305)
plt.xlabel("wavelength (um)")
plt.ylabel("flux normalized")
plt.legend()
plt.show()
