A look at the shadow of HD 332231 b
Here we will be using tracit to look at the shadow of HD 332231 b
First we import tracit. The two run commands are just to setup text rendering in the plots to LaTeX.
[1]:
import tracit
We want to create two dictionaries:
parin which we specify the values for the parameters, priors, boundaries, etc.datin which we give the filenames for our data, and we can also specify the de-trending/noise model
[2]:
par = tracit.par_struct(n_planets=1,n_phot=0,n_spec=1)
dat = tracit.dat_struct(n_phot=0,n_rvs=0,n_ls=1,n_sl=0)
For the parameter structure we specified that we have 1 planet in the system, 1 photometric system, and 1 spectroscopic system (RVs, shadow, or slope).
For the data we also define the number of photometric systems, while we further specify the type of spectroscopic measurements\(-\)n_rvs for RVs, n_ls for the shadow (lineshape), or n_sl for the subplanetary velocity (slope).
Here we want to look at the shadow, so we set n_ls=1.
Below we set the path to our shadow data, and give the indices (counting from 0) for the out-of-transit timestamps.
We also specify a range in the CCF, where there’s no bump. In this case, for instance, the bump of the CCF does not extend out to beyond 12 km/s. It’s of course important to ensure that the value you give is also inside your velocity grid. This is used to de-trend and estimate the scatter of the CCF.
[3]:
dat['LS filename_1'] = 'shadow.hdf5' # filename
dat['idxs_1'] = [39,40,41] # The out-of-transit indices
dat['No_bump_1'] = 12 # No bump in the CCF, i.e., from (vel < - 12 km/s) or (vel > 12 km/s)
dat['Velocity_range_1'] = 20
We then initialize the data, and the business module.
[4]:
tracit.ini_data(dat)
tracit.run_bus(par,dat)
Here we import some results from a previous run to make a nice plot.
[5]:
saved_results = 1 # set to 1 if you want to read in the results instead
if saved_results:
import pandas as pd
rdf = pd.read_csv('results_from_old_fit.csv')
tracit.update_pars(rdf,par,best_fit=False)
If you don’t have results from a previous run, you’ll have to set them individually as shown below.
[6]:
par['T0_b']['Value'] = 2458729.6940
par['lam_b']['Value'] = 0.0
Let’s take a look at the shadow of HD 332231 b. We do this by calling plot_shadow.
[7]:
tracit.run_exp(1)
tracit.plot_shadow(par,dat)
/home/emil/Desktop/PhD/tracit/src/tracit/shady.py:482: RuntimeWarning: divide by zero encountered in true_divide
tan = np.exp(-1*np.power(vel_1d/(zeta*y),2))/y
/home/emil/Desktop/PhD/tracit/src/tracit/shady.py:482: RuntimeWarning: invalid value encountered in true_divide
tan = np.exp(-1*np.power(vel_1d/(zeta*y),2))/y
Using indices [39, 40, 41] as out-of-transit spectra
We can also look at the out-of-transit CCFs. Here we’ll also look at all the in-transit CCFs after subtracting our average out-of-transit CCF, i.e., the 1D shadow.
[8]:
tracit.plot_oot_ccf(par,dat)
/home/emil/Desktop/PhD/tracit/src/tracit/shady.py:482: RuntimeWarning: divide by zero encountered in true_divide
tan = np.exp(-1*np.power(vel_1d/(zeta*y),2))/y
/home/emil/Desktop/PhD/tracit/src/tracit/shady.py:482: RuntimeWarning: invalid value encountered in true_divide
tan = np.exp(-1*np.power(vel_1d/(zeta*y),2))/y
Using indices [39, 40, 41] as out-of-transit spectra
## Spectroscopic system 1/Spectrograph\ 1 ##:
Reduced chi-squared for the oot CCF is:
12.089
Factor to apply to get a reduced chi-squared around 1.0 is:
3.477
Number of data points: 31
Number of fitting parameters: 0
#########################
/home/emil/Desktop/PhD/tracit/src/tracit/shady.py:482: RuntimeWarning: divide by zero encountered in true_divide
tan = np.exp(-1*np.power(vel_1d/(zeta*y),2))/y
/home/emil/Desktop/PhD/tracit/src/tracit/shady.py:482: RuntimeWarning: invalid value encountered in true_divide
tan = np.exp(-1*np.power(vel_1d/(zeta*y),2))/y
If you want to fit the shadow you can do as shown below.
[9]:
dat['Fit LS_1'] = 1 # whether to fit the shadow or not
par['lam_b']['Prior'] = 'uni'
par['lam_b']['Prior_vals'] = [90,2,-180,180]
par['FPs'] = ['lam_b']
mc = 0
if mc:
ndraws = 1500
nwalkers = 10
nproc = 1
ndf = tracit.mcmc(par,dat,ndraws,nwalkers,nproc=1,corner=True,chains=True,path_to_tracit=path)
You can also take a look at the slope. For that you also need to set the LC to true (for now at least).
[10]:
slope = 0
if slope:
par = tracit.par_struct(n_planets=1,n_phot=1,n_spec=1)
dat = tracit.dat_struct(n_phot=1,n_rvs=0,n_ls=0,n_sl=1)
dat['LC filename_1'] = 'lc_HD332231.txt'
import pandas as pd
rdf = pd.read_csv('results_from_old_fit.csv')
tracit.update_pars(rdf,par,best_fit=False)
par['T0_b']['Value'] = 2458729.6940
par['lam_b']['Value'] = 0.0
dat['SLs'] = 1
dat['SL filename_1'] = 'shadow.hdf5' # filename
dat['idxs_1'] = [39,40,41] # The out-of-transit indices
dat['No_bump_1'] = 12 # No bump in the CCF, i.e., from (vel < - 12 km/s) or (vel > 12 km/s)
tracit.ini_data(dat)
tracit.run_bus(par,dat)
tracit.plot_slope(par,dat)