"""
Cross-correlation and radial velocity utility functions for GUIBRUSHR.
This module contains functions for computing planetary and stellar radial
velocities, including support for eccentric orbits and the Rossiter-McLaughlin
(Doppler shadow) effect.
"""
import numpy as np
from GUIBRUSHR.General_Constants.FunctionsAndConstants.Constant_Variables import ConstantVariables
from GUIBRUSHR.General_Constants.FunctionsAndConstants.radvel_g import rv_drive, timetrans_to_timeperi
[docs]
def rv_planet_and_star(
eccentricity, target, ecc_retrieval, opi_retrieval, kp_retrieval, night_obj, mode,
index_night_in_current_instrument, dVsys_arr, counter_derivative_params,
rv_dynamic=0
):
"""
Compute planetary and stellar radial velocities for a given night.
Args:
eccentricity (bool): Whether to use eccentric orbit model.
target: Target object with orbital parameters.
ecc_retrieval (float): Retrieved eccentricity value.
opi_retrieval (float): Retrieved argument of periastron.
kp_retrieval (float): Retrieved planetary radial velocity semi-amplitude (km/s).
night_obj: Night object containing phase and barycentric velocity data.
mode (str): Operating mode ('CC', 'Retrieval', or synthetic).
index_night_in_current_instrument (int): Index of the night within the instrument.
dVsys_arr: Array of systemic velocity offsets per night.
counter_derivative_params (int): Counter for derivative parameter indexing.
rv_dynamic (float, optional): Extra radial velocity (km/s) folded into the
planetary vtot. Used ONLY by the retrieval, where it carries the value of
the retrieval rv parameter; it replaces the former convolution-kernel rv
shift so that the radial velocity is computed in a single place. Defaults
to 0 (synthetic and shuffle pipelines). Forwarded to
calculate_rv_planet_and_star.
Returns:
tuple: (vtot, vtot_star, kpph) planetary velocity, stellar velocity, and Kp*phase arrays.
"""
phases = night_obj.phases_considered
hjdk = night_obj.hjd_ck_considered
barycentric = night_obj.barycentric_velocity_considered
if mode == "CC":
# Use the eccentricity-corrected semi-amplitude when the eccentric orbit
# model is selected; fall back to the circular kp otherwise (or if the
# target predates kp_ecc).
kp = getattr(target, "kp_ecc", None) if (eccentricity and getattr(target, "kp_ecc", None) is not None) else target.kp
ecc = target.eccentricity
opi = target.argument_of_periastron
coeff_vtot = -1
elif mode == "Retrieval":
kp = kp_retrieval
ecc = ecc_retrieval
opi = opi_retrieval
coeff_vtot = 1
else: # synthetic generation
kp = getattr(target, "kp_ecc", None) if (eccentricity and getattr(target, "kp_ecc", None) is not None) else target.kp
ecc = target.eccentricity
opi = target.argument_of_periastron
coeff_vtot = 1
phases = night_obj.phases
hjdk = night_obj.hjd_ck
barycentric = night_obj.barycentric_velocity
# SYSTEMIC VELOCITY CORRECTIONS
# Apply night-specific systemic velocity offset (first night is reference)
dvsys = 0 if (index_night_in_current_instrument == 0 or dVsys_arr is None or mode != "Retrieval") else dVsys_arr[counter_derivative_params]
if -999 in night_obj.hjd_ck:
hjd0 = 0
period = 1
times_of_observation = phases
else:
hjd0 = target.hjd0
period = target.period
times_of_observation = hjdk
systemic_velocity = target.systemic_velocity
ks = target.ks
if ks == 0:
ks = target.kp * target.mass * np.sin(np.deg2rad(target.inclination)) / (
target.stellar_mass * ConstantVariables.SOLAR_TO_JUPITER_MASSES)
if mode != "Retrieval":
# TODO REMOVE
print(f"Insert ks in yaml! ks_km_s: {ks}")
vtot, vtot_star, kpph = calculate_rv_planet_and_star(
hjd0, times_of_observation, eccentricity, period, ecc, opi, kp, ks, phases,
coeff_vtot, systemic_velocity, barycentric, dvsys, rv_dynamic=rv_dynamic
)
return vtot, vtot_star, kpph
[docs]
def calculate_rv_planet_and_star(
hjd0, times_of_observation, eccentricity, period, ecc, opi, kp, ks, phases,
coeff_vtot, systemic_velocity, barycentric, dvsys, vsys_offset_ccf=0, rv_dynamic=0
):
"""
Calculate planetary and stellar radial velocities from orbital parameters.
This is the single place where the planetary total velocity vtot is assembled
from all its additive components (orbital kpph, systemic velocity, per-night
dvsys, barycentric). Two optional additive terms are kept conceptually distinct:
- vsys_offset_ccf: extra systemic-velocity offset used ONLY by the
cross-correlation (CCF) pipeline to scan the Kp-Vsys plane. Default 0.
- rv_dynamic: extra radial velocity used ONLY by the retrieval, where it carries
the retrieval rv parameter. It replaces the former convolution-kernel rv shift,
keeping the radial velocity in one calculation. Default 0.
Sign convention: both extra terms sit inside the subtracted group, so with
coeff_vtot = 1 (Retrieval/synthetic) a positive rv_dynamic lowers vtot
(vtot = vtot_without_rv - rv_dynamic). This reproduces exactly the previous
behaviour where the kernel was translated by +rv. In CC mode coeff_vtot = -1,
so the CCF offset adds with the opposite sense, unchanged from before.
Args:
vsys_offset_ccf (float, optional): CCF-only systemic offset (km/s). Default 0.
rv_dynamic (float, optional): Retrieval-only extra radial velocity (km/s),
equal to the retrieval rv parameter. Default 0.
Returns:
tuple: (vtot, vtot_star, kpph) total velocity shift for planet, for star, and Kp*phase.
"""
if eccentricity:
# opi stores the PLANET's argument of periastron. timetrans_to_timeperi needs the
# STAR's omega (opi - pi); rv_drive needs the PLANET's omega (it yields the
# planet's RV directly). The pi offset between the two is intentional.
omega_planet = opi
omega_star = opi - np.pi
T_periastro = timetrans_to_timeperi(
hjd0, period, ecc, omega_star,
)
orbital_element = np.array([
period, T_periastro, ecc, omega_planet, kp
])
kpph = rv_drive(
times_of_observation, orbital_element,
)
else:
# Simple circular orbit approximation
kpph = kp * np.sin(2 * np.pi * phases)
# Calculate total velocity shift for planet spectrum.
# vsys_offset_ccf (CCF only) and rv_dynamic (retrieval only) are additive terms
# in the subtracted group; one of them is always 0 depending on the pipeline.
vtot = coeff_vtot * (
-kpph - (systemic_velocity + vsys_offset_ccf + rv_dynamic + dvsys) + barycentric
)
# Calculate stellar radial velocity for stellar spectrum correction
if ks is not None:
vtot_star = (
ks * np.sin(2 * np.pi * phases) - systemic_velocity + barycentric
)
else:
vtot_star = None
return vtot, vtot_star, kpph
[docs]
def rv_DopplerShadow(target, phases):
"""
Compute the Rossiter-McLaughlin (Doppler shadow) radial velocity anomaly.
Based on Cegla, H. M., et al. 2016, A&A, 588, A127.
Args:
target: Target object with projected_obliquity, a_Rs_ratio, inclination, v_sini,
and systemic_velocity attributes.
phases (np.ndarray): Orbital phase array.
Returns:
np.ndarray: Radial velocity anomaly plus systemic velocity (km/s).
"""
# formulae taken from
# Cegla, H. M., et al. 2016, A&A, 588, A127
# ---- center of the planet at any given orbital phase ----#
xp = target.a_Rs_ratio * np.sin(2 * np.pi * phases) # xp=a/R* sin(2*pi*phi)
yp = -target.a_Rs_ratio * np.cos(2 * np.pi * phases) * np.cos(np.deg2rad(target.inclination)) # yp=-a/R* con(2*pi*phi)*cos(ip)
l = np.deg2rad(target.projected_obliquity)
# ---- orthogonal distance from the spin-axis ----#
x = xp * np.cos(l) - yp * np.sin(l) # x_|_=xp * cos(lamda) - yp * sin(lamda)
# IF differential rotation and alpha is known
# diff. rot: Omega = Omega_eq*(1-alpha*sin2theta)
# y = xp*np.sin(l.radian) - yp*np.cos(l.radian)
# z = np.sqrt(1 - (x**2) - (y**2))
# beta = np.pi/2 + istar.radian.value
# z1 = z*np.cos(beta) - y*sin(beta)
# y1 = z*np.sin(beta) - y*cos(beta)
# v = x*vsini*(1-alpha*(y1**2))
v = x * target.v_sini
return v + target.systemic_velocity