import numpy as np
import scipy.optimize as optimize
from ForMoSA.utils_spec import continuum_estimate
[docs]
def hc_model(global_params, obs_dict, flx_mod_spectro, indobs):
'''
For high-contrast companions, where the star speckles signal contaminate the data
Args:
global_params (object): Class containing each parameter used in ForMoSA
obs_dict (dict): Dictionay containing all the observationnal entries (photometry, spectroscopy and/or optional)
flx_mod_spectro (array): Model flux of the companion
indobs (int): Index of the current observation loop
Returns:
- contributions (array): Results of the high-constrast model
- flx_mod_spectro (array): Model of the high-constrast model
- flx_mod_spectro (array): Model of the high-constrast model
- speckles (array): Speckles contribution
Authors: Allan Denis
'''
# Initialize contributions
contributions = 1
# Estimate the continuum of the model
flx_mod_spectro *= obs_dict['transm']
flx_mod_spectro_cont = continuum_estimate(obs_dict['wav_spectro'],
flx_mod_spectro,
obs_dict['res_spectro'],
global_params.wav_cont[indobs % len(global_params.wav_cont)],
float(global_params.res_cont[indobs % len(global_params.res_cont)]))
if global_params.hc_bounds_lsq[indobs % len(global_params.hc_bounds_lsq)] != 'NA':
bounds =(float(global_params.hc_bounds_lsq[indobs*2 % len(global_params.hc_bounds_lsq)]),
float(global_params.hc_bounds_lsq[indobs*2 + 1 % len(global_params.hc_bounds_lsq)]))
else:
bounds = (0, np.inf)
# Select model
if global_params.logL_type[indobs % len(global_params.logL_type)].startswith('chi2'):
contributions, flx_mod_spectro, speckle, systematics = hc_model_estimate_speckles(obs_dict['flx_spectro'],
obs_dict['flx_spectro_cont'],
obs_dict['star_flx'],
obs_dict['star_flx_cont'],
flx_mod_spectro,
flx_mod_spectro_cont,
obs_dict['err_spectro'],
bounds,
obs_dict['system'])
else:
flx_mod_spectro, speckle, systematics = hc_model_remove_speckles(obs_dict['flx_spectro'],
obs_dict['flx_spectro_cont'],
obs_dict['star_flx'],
obs_dict['star_flx_cont'],
flx_mod_spectro,
flx_mod_spectro_cont,
obs_dict['err_spectro'])
return contributions, flx_mod_spectro, speckle, systematics
# ----------------------------------------------------------------------------------------------------------------------
[docs]
def hc_model_remove_speckles(flx_obs_spectro, flx_cont_obs_spectro, star_flx_obs_spectro, star_flx_cont_obs_spectro, flx_mod_spectro, flx_cont_mod_spectro, err_obs_spectro):
'''
high-constrast
Args:
flx_obs_spectro (array): Flux of the data
flx_cont_obs_spectro (array): Continuum of the data
star_flx_obs_spectro (array): Flux of the star data
star_flx_cont_obs_spectro (array): Continuum of the star data
flx_mod_spectro (array): Model of the companion
flx_cont_mod_spectro (array): Continuum of the model of the companion
err_obs_spectro (array): Noise of the data
Returns:
- flx_mod_spectro (array): High-resolution content of planet model
- speckles (array): Speckles contribution
- systematics (float): Not taken into account here
Authors: Allan Denis
'''
speckles = star_flx_obs_spectro[:, len(star_flx_obs_spectro[0]) // 2] / star_flx_cont_obs_spectro * flx_cont_obs_spectro # Speckles modulation (Bidot et al. 2023, Landman et al. 2023)
alpha = np.sum(((flx_obs_spectro - speckles) * flx_mod_spectro) / err_obs_spectro**2) / np.sum((flx_mod_spectro**2) / err_obs_spectro**2)
flx_mod_spectro = alpha * (flx_mod_spectro - flx_cont_mod_spectro)
systematics = 0
return flx_mod_spectro, speckles, systematics
[docs]
def hc_model_estimate_speckles(flx_obs_spectro, flx_cont_obs_spectro, star_flx_obs_spectro, star_flx_cont_obs_spectro, flx_mod_spectro, flx_cont_mod_spectro, err, bounds, system_obs_spectro):
'''
high-constrast model of planet and star contributions
Args:
flx_obs_spectro (array): Flux of the data
flx_cont_obs_spectro (array): Continuum of the data
star_flx_obs_spectro (array): Flux of the star data
star_flx_cont_obs_spectro (array): Continuum of the star data
flx_mod_spectro (array): Model of the companion
flx_cont_mod_spectro (array): Continuum of the model of the companion
weights (array): Weights to apply to the data
bounds (tuple): Bounds to be applied to the estimated parameters
system_obs_spectro (array): Systematics
Returns:
- results.x (array): Results of the high-constrast model
- flx_mod_spectro (array): Model of the high-constrast model
- speckles (array): Speckles contribution
Authors: Allan Denis
'''
speckles = np.ones((star_flx_obs_spectro.shape))
for star_i in range(len(star_flx_obs_spectro[0])):
speckles[:,star_i] = star_flx_obs_spectro[:,star_i] / star_flx_cont_obs_spectro
weights = 1 / (err**2)
ind_star = 1 + len(speckles[0])
if len(system_obs_spectro) > 0:
ind_system = ind_star + len(system_obs_spectro[0])
else:
ind_system = ind_star
# # # # # # Solve linear Least Squares A.x = b
# Build matrix A
A = np.zeros([np.size(flx_obs_spectro), ind_system])
A[:, 0] = weights * (flx_mod_spectro - flx_cont_mod_spectro * speckles[:, len(speckles[0]) // 2])
for star_i in range(len(speckles[0])):
A[:, star_i + 1] = weights * (speckles[:,star_i] * flx_cont_obs_spectro)
for system_i in range(ind_system - ind_star):
A[:, system_i + ind_star] = weights * system_obs_spectro[:, system_i]
# Build vector b
b = weights * flx_obs_spectro
# Solve linear Least Squares
results = optimize.lsq_linear(A, b, bounds=bounds)
# Model
flx_mod_spectro = np.dot(A[:, 0], results.x[0]) / weights
# Speckles
speckles = np.dot(A[:, 1:ind_star], results.x[1:ind_star]) / weights
# Systematics
systematics = np.dot(A[:,ind_star:], results.x[ind_star:]) / weights
return results.x, flx_mod_spectro, speckles, systematics