nespinoza / juliet

A versatile modelling tool for transiting and non-transiting (single and multiple) exoplanetary systems
MIT License
55 stars 30 forks source link

Integrating Scalpels #122

Open LucaNap opened 1 month ago

LucaNap commented 1 month ago

Hi everyone,

I am exploring the possibility of adding Scalpels to juliet (for who's not familiar with it, check Cameron+2020), and I am wondering how should it be done. From what I understand, juliet should first read the base Scalpels vector (U0) from a file, and then simply calculate the residuals during the fit as:

res = (RV_values - matrix_N@matrix_C) * (RV_values - matrix_N@matrix_C)

where

matrix_N = np.array([ np.ones(len(RV_times)), rv_drive(RV_times, p1), ..., U0]).T
matrix_C = np.array([offset, 1, ..., alpha]) 

rv_drive is the planetary signal model evaluated with Radvel for planet 1 (and so on, if there are others), and p1 are the parameters of the first planet. Of course, for each planet there should be a new rv_drive in matrix_N and "1" in matric_C. Finally, "offset" and "alpha" are the new priors related to the Scalpels vector (if U0 is mono-dimensional, otherwise there should be an "alpha" for each dimension of the vector). For example, if we have N radial-velocity points, one planet and a mono-dimensional U0, then:

matrix_N@matrix_C = 1*offset + K*1 + U0*alpha = offset + K + U0*alpha.

where matrix_N dimensions are Nx3 and matrix_C 3x1. So, I guess the main change should be the addition of a correction factor to K during the sampling. @nespinoza, what do you think?