Modeling a High Resolution Reflection Spectrum
Hajime Kawahara March 20th (2024)
from jax import config
config.update("jax_enable_x64", True)
Forget about the following code. I need this to run the code somewhere…
#username="exoplanet01"
username="kawahara"
if username == "exoplanet01":
import ssl
ssl._create_default_https_context = ssl._create_unverified_context
We need to get JoviSpec, the repository that contains the real
reflection spectra of Jupiter.
git clone https://github.com/HajimeKawahara/jovispec.git
python setup.py install
from jovispec import abcio
import pkg_resources
jupiter_data = pkg_resources.resource_filename("jovispec", "jupiter_data")
# red
rlambc, rspecc, rheadc = abcio.read_qfits("06033", jupiter_data, ext="q")
rlambw, rspecw, rheadw = abcio.read_qfits("06047", jupiter_data, ext="q")
rlambe, rspece, rheade = abcio.read_qfits("06049", jupiter_data, ext="q")
# rapid rotator
rlamb_ref, rspec_ref, rhead_ref = abcio.read_qfits(
"06031", jupiter_data, ext="q"
) # HD13041
# blue
# rlambc, rspecc, rheadc=abcio.read_qfits("06034",jupiter_data,ext="q")
# rlambw, rspecw, rheaadw=abcio.read_qfits("06048",jupiter_data,ext="q")
# rlambe, rspece, rheade=abcio.read_qfits("06050",jupiter_data,ext="q")
WARNING: VerifyWarning: Invalid 'BLANK' keyword in header. The 'BLANK' keyword is only applicable to integer data, and will be ignored in this HDU. [astropy.io.fits.hdu.image]
Sets the wavelength range
wavelength_start = 8500.0 #AA
wavelength_end = 8800.0 #AA
You find almost no telluric lines in this wavelength range
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10,2.5))
ax = fig.add_subplot(111)
plt.plot(rlamb_ref,rspec_ref,alpha=0.5,label="reference star")
plt.ylim(0.0,0.012)
plt.xlim(wavelength_start,wavelength_end)
plt.legend()
plt.show()
We observed the east edge, center, and west edge of the Jovian disk. The rapid spin rotation yields different Doppler shift for these observataions
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10,7.5))
ax = fig.add_subplot(311)
plt.plot(rlambc,rspecc,alpha=0.5,label="CENTER")
#plt.plot(rlambe,rspece,alpha=0.5,label="EAST")
#plt.plot(rlambw,rspecw*0.7,alpha=0.5,label="WEST")
plt.ylim(0.0,0.3)
plt.xlim(wavelength_start,wavelength_end)
plt.legend()
plt.title("Jupiter CH4 band")
#plt.xlim(7150.0,7200.0)
ax = fig.add_subplot(312)
plt.plot(rlambc,rspecc,alpha=0.5,label="CENTER")
plt.plot(rlambe,rspece,alpha=0.5,label="EAST")
plt.plot(rlambw,rspecw*0.7,alpha=0.5,label="WEST")
plt.legend()
plt.ylim(0.0,0.2)
plt.xlim(wavelength_start,wavelength_end)
plt.xlim(8625,8650)
ax = fig.add_subplot(313)
plt.plot(rlambc,rspecc,alpha=0.5,label="CENTER")
plt.plot(rlambe,rspece,alpha=0.5,label="EAST")
plt.plot(rlambw,rspecw*0.7,alpha=0.5,label="WEST")
plt.legend()
plt.ylim(0.0,0.2)
plt.xlim(wavelength_start,wavelength_end)
plt.xlabel("Wavelength (Angstrom)")
plt.savefig("ch4jupiter.png")
Some messy procedure for astronomers…
import numpy as np
rlamb = rlambc
rspec = rspecc
# mask some bad regions... as usual in astronomy
mask_wav = [
[8565.4, 8586.5],
[8689.5, 8713.0]
]
rlamb = np.array([float(d) for d in rlamb])
mask_index=np.digitize(mask_wav,rlamb)
for ind in mask_index:
rspec[ind[0]:ind[1]+1] = None
# None for outside wvelength start - end region
start_index=np.digitize(wavelength_start,rlamb)
end_index=np.digitize(wavelength_end,rlamb)
rspec[:start_index] = None
rspec[end_index:] = None
#Air-Vaccum correction
from specutils.utils.wcs_utils import refraction_index
import astropy.units as u
mask = rspec == rspec
rlamb = rlamb[mask]
rspec = rspec[mask]
nair = refraction_index(rlamb*u.AA,method="Ciddor1996")
rlamb = rlamb*nair
# ascending wavenumber form
from exojax.spec.unitconvert import wav2nu
rlamb = rlamb[::-1]
nus_obs = wav2nu(rlamb, unit="AA")
rspec = rspec[::-1]
The following spectrum is the target we will analyze (the center part of Jupiter)!
fig = plt.figure(figsize=(20,5))
plt.plot(nus_obs,rspec, label="masked spectrum")
plt.xlabel("wavenumber (cm-1)")
plt.legend()
plt.title("masked spectrum")
plt.show()
Telluric
Telluric lines stuff… ExoJAX can do that anyway.
from exojax.spec.api import MdbHitran
from exojax.spec.opacalc import OpaDirect
from exojax.spec.opacalc import OpaPremodit
from exojax.utils.grids import wavenumber_grid
from exojax.spec.unitconvert import wav2nu
N = 40000
margin = 10 # cm-1
nus_start = wav2nu(wavelength_end, unit="AA") - margin
nus_end = wav2nu(wavelength_start, unit="AA") + margin
mdb_water = MdbHitran("H2O", nurange=[nus_start, nus_end], isotope=1)
nus, wav, res = wavenumber_grid(nus_start, nus_end, N, xsmode="lpf", unit="cm-1")
opa_telluric = OpaPremodit(mdb_water, nu_grid=nus, allow_32bit=True, auto_trange=[80.0, 400.0])
xsmode = lpf xsmode assumes ESLOG in wavenumber space: mode=lpf ====================================================================== 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 ====================================================================== OpaPremodit: params automatically set. default elower grid trange (degt) file version: 2 Robust range: 79.45501192821337 - 740.1245313998245 K Change the reference temperature from 296.0K to 91.89455622053987 K.
/home/kawahara/exojax/src/exojax/spec/set_ditgrid.py:52: UserWarning: There exists negative or zero value.
warnings.warn("There exists negative or zero value.")
OpaPremodit: Tref_broadening is set to 178.88543819998304 K
OpaPremodit: gamma_air and n_air are used. gamma_ref = gamma_air/Patm
# of reference width grid : 22
# of temperature exponent grid : 7
uniqidx: 100%|██████████| 40/40 [00:00<00:00, 8397.43it/s]
Premodit: Twt= 328.42341041740974 K Tref= 91.89455622053987 K
Making LSD:|##########----------| 50%
Making LSD:|####################| 100%
Yes ah no, almost no telluric line confirmed…
import jax.numpy as jnp
T = 200.0
P = 0.5
xsv = opa_telluric.xsvector(T, P)
nl = 1.0e22
a = 0.25
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111)
plt.title("masked spectrum w/ telluric initial guess")
plt.plot(nus_obs, rspec)
#plt.ylim(0.0, 0.7)
plt.plot(nus, a * jnp.exp(-nl * xsv), alpha=0.5)
# plt.xlim(1.e8/wavelength_start, 1.e8/wavelength_end)
plt.xlabel("wavenumber (cm-1)")
plt.show()
Solar spectrum
This is the reflected “solar” spectrum by Jupiter! So, we need the solar spectrum template.
I found very good one: High-resolution solar spectrum taken from Meftar et al. (2023)
10.21413/SOLAR-HRS-DATASET.V1.1_LATMOS
import pandas as pd
#filename = "/home/kawahara/solar-hrs/Spectre_HR_LATMOS_Meftah_V1.txt"
filename = "/home/"+username+"/solar-hrs/Spectre_HR_LATMOS_Meftah_V1.txt"
dat = pd.read_csv(filename, names=("wav","flux"), comment=";", delimiter="\t")
dat["wav"] = dat["wav"]*10
wav_solar = dat["wav"][::-1]
solspec = dat["flux"][::-1]
nus_solar = wav2nu(wav_solar,unit="AA")
This is really the solar spectrum reflected by something…
vrv = (106.5-82.5)
vperc = vrv/300000
fig = plt.figure(figsize=(20,5))
ax = fig.add_subplot(211)
plt.plot(nus_obs,rspec,alpha=0.5,label="CENTER")
plt.plot(nus_solar*(1.0+vperc),solspec*0.17,lw=1,label="Solar")
plt.xlim(np.min(nus_obs),np.max(nus_obs))
plt.ylim(0.0,0.2)
plt.legend()
plt.title("Solar templeate vs obs")
ax = fig.add_subplot(212)
plt.plot(nus_obs,rspec,alpha=0.5,label="CENTER")
plt.plot(nus_solar*(1.0+vperc),solspec*0.17,lw=1,label="Solar")
plt.xlim(11400,11450)
plt.ylim(0.0,0.2)
plt.legend()
plt.xlabel("wavenumber (cm-1)")
Text(0.5, 0, 'wavenumber (cm-1)')
Atmospheric Model Setting
See Jupiter_cloud_model_using_amp.ipynb
from exojax.spec.atmrt import ArtReflectPure
from exojax.utils.astrofunc import gravity_jupiter
art = ArtReflectPure(nu_grid=nus, pressure_btm=1.0e2, pressure_top=1.0e-3, nlayer=100)
art.change_temperature_range(80.0, 400.0)
Tarr = art.powerlaw_temperature(150.0, 0.2)
Parr = art.pressure
mu = 2.3 # mean molecular weight
gravity = gravity_jupiter(1.0, 1.0)
fig = plt.figure()
ax = fig.add_subplot()
ax.plot(Tarr, art.pressure)
ax.invert_yaxis()
plt.yscale("log")
plt.xscale("log")
plt.show()
/home/kawahara/exojax/src/exojax/spec/dtau_mmwl.py:14: FutureWarning: dtau_mmwl might be removed in future.
warnings.warn("dtau_mmwl might be removed in future.", FutureWarning)
from exojax.spec.pardb import PdbCloud
from exojax.atm.atmphys import AmpAmcloud
pdb_nh3 = PdbCloud("NH3")
amp_nh3 = AmpAmcloud(pdb_nh3, bkgatm="H2")
amp_nh3.check_temperature_range(Tarr)
#condensate substance density
rhoc = pdb_nh3.condensate_substance_density #g/cc
from exojax.utils.zsol import nsol
from exojax.atm.mixratio import vmr2mmr
from exojax.spec.molinfo import molmass_isotope
n = nsol() #solar abundance
abundance_nh3 = n["N"]
molmass_nh3 = molmass_isotope("NH3", db_HIT=False)
fsed = 10.
sigmag = 2.0
Kzz = 1.e4
MMRbase_nh3 = vmr2mmr(abundance_nh3, molmass_nh3, mu)
rg_layer, MMRc = amp_nh3.calc_ammodel(Parr, Tarr, mu, molmass_nh3, gravity, fsed, sigmag, Kzz, MMRbase_nh3)
.database/particulates/virga/virga.zip exists. Remove it if you wanna re-download and unzip.
Refractive index file found: .database/particulates/virga/NH3.refrind
Miegrid file does not exist at .database/particulates/virga/miegrid_lognorm_NH3.mg.npz
Generate miegrid file using pdb.generate_miegrid if you use Mie scattering
/home/kawahara/exojax/src/exojax/atm/atmphys.py:50: UserWarning: min temperature 80.0 K is smaller than min(vfactor t range) 179.10000000000002 K
warnings.warn(
import matplotlib.pyplot as plt
# to convert MMR to g/L ...
from exojax.atm.idealgas import number_density
from exojax.utils.constants import m_u
fac = molmass_nh3*m_u*number_density(Parr,Tarr)*1.e3 #g/L
fig = plt.figure()
ax = fig.add_subplot(131)
plt.plot(rg_layer, Parr)
plt.xlabel("rg (cm)")
plt.ylabel("pressure (bar)")
plt.yscale("log")
ax.invert_yaxis()
ax = fig.add_subplot(132)
plt.plot(MMRc, Parr)
plt.xlabel("condensate MMR")
plt.yscale("log")
#plt.xscale("log")
ax.invert_yaxis()
ax = fig.add_subplot(133)
plt.plot(fac*MMRc, Parr)
plt.xlabel("cloud density g/L")
plt.yscale("log")
#plt.xscale("log")
ax.invert_yaxis()
rg = 1.e-5
from exojax.spec.opacont import OpaMie
opa = OpaMie(pdb_nh3, nus)
#sigma_extinction, sigma_scattering, asymmetric_factor = opa.mieparams_vector(rg,sigmag) # if using MieGrid
sigma_extinction, sigma_scattering, asymmetric_factor = opa.mieparams_vector_direct_from_pymiescatt(rg, sigmag)
100%|██████████| 5/5 [00:01<00:00, 4.10it/s]
# plt.plot(pdb_nh3.refraction_index_wavenumber, miepar[50,:,0])
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(311)
plt.plot(nus, asymmetric_factor, color="black")
plt.xscale("log")
plt.ylabel("$g$")
ax = fig.add_subplot(312)
plt.plot(nus, sigma_scattering/sigma_extinction, label="single scattering albedo", color="black")
plt.xscale("log")
plt.ylabel("$\\omega$")
ax = fig.add_subplot(313)
plt.plot(nus, sigma_extinction, label="ext", color="black")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("wavenumber (cm-1)")
plt.ylabel("$\\beta_0$")
plt.savefig("miefig_high.png")
Next, we consider the gas component, i.e. methane
from exojax.spec.api import MdbHitemp
Ignore the following unless you are interested in the elower max optimization
optimize = False #if you want to run the elower maximum optimization, turn True
if optimize:
mdb = MdbHitemp("CH4", nurange=[nus_start,nus_end], isotope=1)
from exojax.spec.optgrid import optelower
Tmax = 400.0
Pmin = 1.e-5
Eopt = optelower(mdb, nus, Tmax, Pmin)
Oh, we need HITEMP because our observation is in visible band, i.e., we need higher energy states.
Eopt = 3300.0 # this is from the above optimization result
mdb_reduced = MdbHitemp("CH4", nurange=[nus_start,nus_end], isotope=1, elower_max=Eopt)
The optimization stuff, ignore
if optimize:
plt.plot(mdb.elower,mdb.line_strength_ref,".",alpha=0.3)
plt.plot(mdb_reduced.elower,mdb_reduced.line_strength_ref,".",alpha=0.5,color="C2")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("E lower (cm-1)")
plt.ylabel("line strength at T="+str(mdb.Tref)+"K")
plt.axvline(Eopt, color="C2")
opa = OpaPremodit(mdb_reduced, nu_grid=nus, allow_32bit=True, auto_trange=[80.0, 300.0]) #this reduced the device memory use in 5/6.
OpaPremodit: params automatically set.
default elower grid trange (degt) file version: 2
Robust range: 79.45501192821337 - 740.1245313998245 K
Change the reference temperature from 296.0K to 91.89455622053987 K.
OpaPremodit: Tref_broadening is set to 154.91933384829665 K
OpaPremodit: gamma_air and n_air are used. gamma_ref = gamma_air/Patm
# of reference width grid : 5
# of temperature exponent grid : 2
uniqidx: 100%|██████████| 3/3 [00:00<00:00, 60.75it/s]
Premodit: Twt= 328.42341041740974 K Tref= 91.89455622053987 K
Making LSD:|####################| 100%
cloud opacity
dtau_cld = art.opacity_profile_cloud_lognormal(sigma_extinction, rhoc, MMRc, rg, sigmag, gravity)
gas opacity
mmr_ch4 = art.constant_mmr_profile(0.02) # (*_*)/~
molmass_ch4 = molmass_isotope("CH4", db_HIT=False)
xsmatrix = opa.xsmatrix(Tarr,Parr)
dtau_ch4 = art.opacity_profile_xs(xsmatrix, mmr_ch4,molmass_ch4,gravity)
dtau_cld_scat = art.opacity_profile_cloud_lognormal(sigma_scattering, rhoc, MMRc, rg, sigmag, gravity)
single_scattering_albedo = (dtau_cld_scat)/(dtau_cld + dtau_ch4)
sigma_scattering[None,:]/sigma_extinction[None,:] + np.zeros((len(art.pressure), len(nus)))
asymmetric_parameter = asymmetric_factor + np.zeros((len(art.pressure), len(nus)))
reflectivity_surface = np.zeros(len(nus))
dtau = dtau_cld + dtau_ch4
from exojax.plot.atmplot import plottau
plottau(np.log10(nus), dtau_ch4, Tarr, Parr)
plt.show()
plottau(np.log10(nus), dtau_cld, Tarr, Parr)
plt.show()
plottau(np.log10(nus), dtau, Tarr, Parr)
OK, uses the solar spectrum as an incomning spectrum!
# vv = 26.0
vv = 0.0
vpercp = (vrv + vv) / 300000
vpercp2 = vv / 300000
# incoming_flux = np.ones(len(nus))
incoming_flux = np.interp(nus, nus_solar * (1.0 + vpercp), solspec)
Fr = art.run(
dtau,
single_scattering_albedo,
asymmetric_parameter,
reflectivity_surface,
incoming_flux,
)
from exojax.spec.specop import SopInstProfile
from exojax.utils.instfunc import R2STD
sop = SopInstProfile(nus)
res_inst = 80000.0
std = R2STD(res_inst)
Fr_inst = sop.ipgauss(Fr, std)
Fr_samp = sop.sampling(Fr_inst, vv, nus)
Recall that this is the forward modelling. So, the following results
look not bad! If you have time change the MMR of CH4 (if you find the
line of (*_*)/~)?
fig = plt.figure(figsize=(30,5))
plt.plot(nus_obs,rspec,alpha=0.5,label="observation (Jupiter center)")
plt.plot(nus,Fr_samp,alpha=0.5,label="forward modelling")
plt.legend()
plt.xlabel("wavenumber (cm-1)")
plt.savefig("Jupiter_forward.png")
plt.show()
In microscopic view, it’s not so accurate (maybe Methane line list inaccuracy?)
fig = plt.figure(figsize=(30,15))
ax = fig.add_subplot(311)
plt.plot(nus_obs,rspec,alpha=0.5,label="observation (Jupiter center)")
plt.plot(nus,Fr_samp,alpha=0.5,label="forward modelling")
plt.plot(nus*(1-vpercp2),incoming_flux*0.2, label="incoming")
plt.legend()
plt.xlim(11550,11600)
ax = fig.add_subplot(312)
plt.plot(nus_obs,rspec,alpha=0.5,label="observation (Jupiter center)")
plt.plot(nus,Fr_samp,alpha=0.5,label="forward modelling")
plt.plot(nus*(1-vpercp2),incoming_flux*0.2, label="incoming")
plt.legend()
plt.xlim(11650,11700)
ax = fig.add_subplot(313)
plt.plot(nus_obs,rspec,alpha=0.5,label="observation (Jupiter center)")
plt.plot(nus,Fr_samp,alpha=0.5,label="forward modelling")
plt.plot(nus*(1-vpercp2),incoming_flux*0.2, label="incoming")
plt.legend()
plt.xlim(11675,11725)
plt.xlabel("wavenumber (cm-1)")
plt.savefig("Jupiter_forward_1.png")
plt.show()
