Transit & Absorption Time-Series¶
In this notebook, we simulate a Sun-like star hosting a Jupiter-sized planet, representing the stellar local spectra with a PHOENIX model. We compute a time series of spectra during the planetary transit within the wavelength range 5882–5902 Å, and from these obtain the corresponding absorption time series in the stellar and planetary rest frames. We then illustrate the evolution of the absorption spectra as the planet moves across the stellar disk, culminating in the presentation of the resulting mean absorption spectrum.
In [1]:
Copied!
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import TwoSlopeNorm
from SOAP.utils import transgrad, compute_planet_doppler_shift, transit_durations
from SOAP.visualizer import plot_absorption_map
import SOAP
from SOAP import PHOENIX
from SOAP.fast_starspot import doppler_shift
from astropy.constants import M_earth, M_sun
from IPython.display import HTML
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import TwoSlopeNorm
from SOAP.utils import transgrad, compute_planet_doppler_shift, transit_durations
from SOAP.visualizer import plot_absorption_map
import SOAP
from SOAP import PHOENIX
from SOAP.fast_starspot import doppler_shift
from astropy.constants import M_earth, M_sun
from IPython.display import HTML
In [2]:
Copied!
# Star-like star, rigid rotation
radius=1 # Stellar radius [solar radii]
mass=1 # Stellar mass [solar masses]
prot=24.47 # Stellar rotation period [days]
incl=45 # Stellar inclination [degrees] -> 90 is equator-on, 0 is pole-on
logg= 4.4 # Stellar surface gravity [cgs]
feh= 0.00 # Stellar metallicity [dex]
Teff=5777 # Stellar effective temperature [K]
ldcn=[[0.5,0.16]] # Linear and quadratic limb-darkening coefficients
start_psi=0.0 # Starting phase of the star
# Mock jupiter-like planet close to the star
Pp=3 # Orbital period [days]
e=0.0 # Orbital eccentricity
w=90.0 # Argument of periastron [degrees]
ip=90 # Orbital inclination [degrees]
lbda=0 # Sky-projected spin-orbit misalignement [degrees]
a=9 # Semi-major axis [stellar radii]
Rp=0.15 # Planetary radius [stellar radii]
Mp=317.8 # Planetary mass [Earth masses]
# Observational parameters
λ =[5882, 5902] # Wavelength range [Angstrom]
Res=140000 # Instrumental resolution
# Star-like star, rigid rotation
radius=1 # Stellar radius [solar radii]
mass=1 # Stellar mass [solar masses]
prot=24.47 # Stellar rotation period [days]
incl=45 # Stellar inclination [degrees] -> 90 is equator-on, 0 is pole-on
logg= 4.4 # Stellar surface gravity [cgs]
feh= 0.00 # Stellar metallicity [dex]
Teff=5777 # Stellar effective temperature [K]
ldcn=[[0.5,0.16]] # Linear and quadratic limb-darkening coefficients
start_psi=0.0 # Starting phase of the star
# Mock jupiter-like planet close to the star
Pp=3 # Orbital period [days]
e=0.0 # Orbital eccentricity
w=90.0 # Argument of periastron [degrees]
ip=90 # Orbital inclination [degrees]
lbda=0 # Sky-projected spin-orbit misalignement [degrees]
a=9 # Semi-major axis [stellar radii]
Rp=0.15 # Planetary radius [stellar radii]
Mp=317.8 # Planetary mass [Earth masses]
# Observational parameters
λ =[5882, 5902] # Wavelength range [Angstrom]
Res=140000 # Instrumental resolution
Where we use will import Phoenix spectra that match the stellar properties. These will be used to compute the RVs in the simulation
In [3]:
Copied!
spQ = PHOENIX(wave_range=λ, normalize=True, teff=Teff, logg=logg, Z=feh)
tplt={"wave":spQ.wave, "flux":spQ.flux} # Template spectrum matching the input local spectra to cumpute the RVs
spQ = PHOENIX(wave_range=λ, normalize=True, teff=Teff, logg=logg, Z=feh)
tplt={"wave":spQ.wave, "flux":spQ.flux} # Template spectrum matching the input local spectra to cumpute the RVs
Start the simulation¶
In [4]:
Copied!
sim = SOAP.Simulation(pixel=spQ, pixel_spot=None, wlll=(λ[1]+λ[0])/2, inst_reso=Res, grid=600, active_regions=[], ring=None)
sim = SOAP.Simulation(pixel=spQ, pixel_spot=None, wlll=(λ[1]+λ[0])/2, inst_reso=Res, grid=600, active_regions=[], ring=None)
In [5]:
Copied!
# Set the properties of the star and planet
sim.star.set(prot=prot,incl=incl,u1=ldcn[0][0],u2= ldcn[0][1], start_psi=start_psi,radius=radius,mass=mass,teff=Teff)
sim.planet.set(P=Pp,t0=start_psi,e=e,w=w,ip=ip,lbda=lbda,a=a,Rp=Rp,Mp=Mp)
# Set the properties of the star and planet
sim.star.set(prot=prot,incl=incl,u1=ldcn[0][0],u2= ldcn[0][1], start_psi=start_psi,radius=radius,mass=mass,teff=Teff)
sim.planet.set(P=Pp,t0=start_psi,e=e,w=w,ip=ip,lbda=lbda,a=a,Rp=Rp,Mp=Mp)
In [6]:
Copied!
# Provide the range of stellar phases to simulate
psi = np.linspace(-0.1,0.1,50) / sim.star.prot
# Provide the range of stellar phases to simulate
psi = np.linspace(-0.1,0.1,50) / sim.star.prot
In [7]:
Copied!
# Calculate the transit signal
out=sim.calculate_signal(psi, renormalize_rv=True, template=tplt)
# Calculate the transit signal
out=sim.calculate_signal(psi, renormalize_rv=True, template=tplt)
In [8]:
Copied!
# Convert from transmission to absorption [%]
absorption_spec=(1-sim.pixel_trans)*100
# Convert from transmission to absorption [%]
absorption_spec=(1-sim.pixel_trans)*100
In [9]:
Copied!
# Plot the absorption map in the stellar rest frame
fig, ax = plot_absorption_map(sim,psi, absorption_spec,λ)
plt.show()
# Plot the absorption map in the stellar rest frame
fig, ax = plot_absorption_map(sim,psi, absorption_spec,λ)
plt.show()
In [10]:
Copied!
# Doppler shift of the map to the planetary rest frame
absorpt_prest=compute_planet_doppler_shift(sim,psi,absorption_spec)
# Doppler shift of the map to the planetary rest frame
absorpt_prest=compute_planet_doppler_shift(sim,psi,absorption_spec)
In [11]:
Copied!
# Plot the absorption map in the planet rest frame
fig, ax = plot_absorption_map(sim,psi,absorpt_prest,λ)
plt.show()
# Plot the absorption map in the planet rest frame
fig, ax = plot_absorption_map(sim,psi,absorpt_prest,λ)
plt.show()
In [12]:
Copied!
# Create an animation of the transit with the absorption map
ani = sim.visualize_animation(
output=out,
plot_type="tr",
lim=[5889.950-1.5, 5889.950+1.5],
ref_wave=5889.950,
plot_lims=[-1, 1]
)
HTML(ani.to_jshtml())
# Create an animation of the transit with the absorption map
ani = sim.visualize_animation(
output=out,
plot_type="tr",
lim=[5889.950-1.5, 5889.950+1.5],
ref_wave=5889.950,
plot_lims=[-1, 1]
)
HTML(ani.to_jshtml())
Out[12]: