ravest.model ============ .. py:module:: ravest.model .. autoapi-nested-parse:: Radial velocity models for planets, stars, and trends. Provides classes for modelling radial velocity signals from planetary orbits and stellar trends (constant, linear, quadratic). Classes ------- .. autoapisummary:: ravest.model.Planet ravest.model.Instrument ravest.model.Trend ravest.model.Star Functions --------- .. autoapisummary:: ravest.model._solve_kepler ravest.model._true_anomaly ravest.model._radial_velocity_from_f ravest.model._njit_kepler_rv ravest.model._compute_rv ravest.model.calculate_mpsini ravest.model.fold_time_series Module Contents --------------- .. py:function:: _solve_kepler(Mi: float, e: float) -> tuple[float, float] Solve Kepler's equation for a single mean anomaly value. Kepler's equation :math:`E - e\sin E = M` must be solved iteratively for the eccentric anomaly :math:`E`. Uses Halley's method (cubic convergence) with the same tolerance and iteration limits as scipy's newton solver. :param Mi: Mean anomaly at a single time (radians). :type Mi: float :param e: Eccentricity of the orbit, 0 < e < 1 (dimensionless). :type e: float :returns: * **cos_E** (*float*) -- Cosine of the eccentric anomaly. * **sin_E** (*float*) -- Sine of the eccentric anomaly. .. py:function:: _true_anomaly(cos_E: float, sin_E: float, e: float, sqrt_1me2: float) -> tuple[float, float] Compute cos(f) and sin(f) of the true anomaly from the eccentric anomaly. Calculates the true anomaly directly from cos(E) and sin(E) rather than via the arctan formula, avoiding a trig call per element. :param cos_E: Cosine of the eccentric anomaly. :type cos_E: float :param sin_E: Sine of the eccentric anomaly. :type sin_E: float :param e: Eccentricity of the orbit (dimensionless). :type e: float :param sqrt_1me2: Precomputed :math:`\sqrt{1 - e^2}`. :type sqrt_1me2: float :returns: * **cos_f** (*float*) -- Cosine of the true anomaly. * **sin_f** (*float*) -- Sine of the true anomaly. .. rubric:: Notes Rather than computing :math:`f` as an angle via the arctan formula :math:`f = 2\arctan\left(\sqrt{\frac{1+e}{1-e}} \tan\frac{E}{2}\right)`, we obtain :math:`\cos f` and :math:`\sin f` directly from :math:`\cos E` and :math:`\sin E`: .. math:: \cos f = \frac{\cos E - e}{1 - e\cos E} \sin f = \frac{\sqrt{1 - e^2} \sin E}{1 - e\cos E} The :math:`\cos f` identity is a standard result from the geometry of the auxiliary circle (Perryman, Exoplanet Handbook, eq. 2.6). :math:`\sin f` follows from :math:`\sin^2 f = 1 - \cos^2 f`, which simplifies to :math:`(1 - e^2)\sin^2 E / (1 - e\cos E)^2`. Since :math:`\cos E` and :math:`\sin E` are already available from the Kepler solver, this avoids any additional trig calls. .. py:function:: _radial_velocity_from_f(cos_f: float, sin_f: float, K: float, cos_w: float, sin_w: float, e_cos_w: float) -> float Compute radial velocity from the true anomaly of a single observation. Uses the trig identity :math:`\cos(f + \omega) = \cos f \cos \omega - \sin f \sin \omega` to avoid computing :math:`f` explicitly. :param cos_f: Cosine of the true anomaly. :type cos_f: float :param sin_f: Sine of the true anomaly. :type sin_f: float :param K: Semi-amplitude of the stellar radial velocity (m/s). :type K: float :param cos_w: Precomputed cos(w). :type cos_w: float :param sin_w: Precomputed sin(w). :type sin_w: float :param e_cos_w: Precomputed e*cos(w). :type e_cos_w: float :returns: Radial velocity (m/s). :rtype: float .. rubric:: Notes The RV equation is: .. math:: v_r = K \left[ \cos(f + \omega) + e\cos\omega \right] Rather than computing the angle :math:`f` and evaluating :math:`\cos(f + \omega)`, we expand via the addition identity :math:`\cos(f + \omega) = \cos f \cos \omega - \sin f \sin \omega`, using :math:`\cos f` and :math:`\sin f` from ``_true_anomaly`` directly. :math:`\cos\omega`, :math:`\sin\omega`, and :math:`e\cos\omega` are constant across all observations and precomputed by the caller. .. py:function:: _njit_kepler_rv(M: numpy.ndarray, e: float, K: float, w: float) -> numpy.ndarray Solve Kepler's equation and compute radial velocity for an array of times. Loops over mean anomaly values element-by-element, calling the scalar Kepler solver, true anomaly, and RV functions for each. Numba inlines these calls at compile time, so performance is identical to a single monolithic loop while keeping each step readable. :param M: Mean anomaly at each time (radians). :type M: np.ndarray :param e: Eccentricity of the orbit, 0 < e < 1 (dimensionless). :type e: float :param K: Semi-amplitude of the stellar radial velocity (m/s). :type K: float :param w: Argument of periapsis of the star (radians). :type w: float :returns: Radial velocity of the star due to the planet at each time (m/s). :rtype: np.ndarray .. py:function:: _compute_rv(M: numpy.ndarray, e: float, K: float, w: float) -> numpy.ndarray Compute radial velocity from mean anomaly, dispatching by eccentricity. For circular orbits (e=0) the RV simplifies to K*cos(M + w) and can be computed directly with numpy. For eccentric orbits, delegates to the numba-compiled solver. :param M: Mean anomaly at each time (radians). :type M: np.ndarray :param e: Eccentricity of the orbit, 0 <= e < 1 (dimensionless). :type e: float :param K: Semi-amplitude of the stellar radial velocity (m/s). :type K: float :param w: Argument of periapsis of the star (radians). :type w: float :returns: Radial velocity of the star due to the planet at each time (m/s). :rtype: np.ndarray .. py:class:: Planet(letter: str, parameterisation: ravest.param.Parameterisation, params: dict[str, float]) Planet defined by its orbital parameters. :param letter: The label of the planet, e.g. "b", "c". Must be a single letter. :type letter: `str` :param parameterisation: The set of planetary parameters used to define the planet. :type parameterisation: `parameterisation` :param params: The orbital parameters, matching the parameterisation. :type params: `dict` .. py:attribute:: letter .. py:attribute:: parameterisation .. py:attribute:: params .. py:attribute:: _rvparams .. py:method:: __repr__() -> str .. py:method:: __str__() -> str .. py:method:: _calculate_mean_motion(period: float) -> float Calculate mean motion (mean angular rate of orbit in radians/day). This is the mean angular rate calculated by dividing the orbital period by 2*pi radians. :param period: The orbital period in days. :type period: `float` :returns: The mean motion, the average angular rate of orbit (radians/day). :rtype: `float` .. py:method:: _calculate_mean_anomaly(t: numpy.ndarray, n: float, time_peri: float) -> numpy.ndarray Calculate mean anomaly (radians). For an eccentric orbit with period P, mean anomaly is a fictitious angle that increases linearly with time ``t`` for an angular rate ``n``, the mean motion. It is the angle from the point of periapsis passage that would have been swept in a circular orbit of period ``P`` with fixed angular rate ``n``. :param n: The mean motion (radians/day). :type n: `float` :param t: The time to evaluate the mean anomaly at (day). :type t: `float` :param time_peri: The time of periapsis of the star's orbit (day). :type time_peri: `float` :returns: The mean anomaly (radians). :rtype: `float` .. py:method:: radial_velocity(t: numpy.ndarray) -> numpy.ndarray Calculate radial velocity of the star due to the planet, at time t. Uses a numba-compiled Kepler solver that solves Kepler's equation and computes the radial velocity in a single pass. :param t: The time to calculate the radial velocity at (day). :type t: np.ndarray :returns: Radial velocity of the reflex motion of star due to the planet (m/s). :rtype: np.ndarray .. py:method:: mpsini(mass_star: float, unit: str = 'kg') -> float Calculate the minimum mass of the planet. :param mass_star: The mass of the star in solar masses. :type mass_star: `float` :param unit: The unit to return the planetary minimum mass in. Options are "kg", "M_earth", "M_jupiter". :type unit: `str` :returns: The minimum mass of the planet (solar masses). :rtype: `float` .. py:class:: Instrument(name: str, g: float, jit: float) Instrument-specific parameters for RV observations. Represents the measurement characteristics of a specific instrument/telescope, including its velocity offset (gamma) and additional noise (jitter). :param name: The name/identifier of the instrument (e.g., "HARPS"). Must match the labels used in the data's instrument column. :type name: str :param g: Gamma offset - the constant RV offset for this instrument [m/s]. :type g: float :param jit: Jitter - additional noise added in quadrature to uncertainties [m/s]. Must be >= 0. :type jit: float .. rubric:: Examples >>> harps = Instrument("HARPS", g=5.0, jit=2.0) >>> print(harps) Instrument HARPS: γ=5.0 m/s, jitter=2.0 m/s >>> hires = Instrument("HIRES", g=-3.6, jit=1.5) >>> hires.g -3.6 .. py:attribute:: name .. py:attribute:: g .. py:attribute:: jit .. py:method:: __repr__() -> str .. py:method:: __str__() -> str .. py:class:: Trend(t0: float, params: dict[str, float]) System-wide trend in the radial velocity of the star. Represents long-term velocity trends that apply across all instruments, such as acceleration from a distant companion. The constant offset (gamma) is now handled per-instrument via the Instrument class. :param t0: The reference zero-point time for the linear and quadratic trend. Recommended to be the mean of the input times. :type t0: float :param params: The parameters of the trend: linear and quadratic components. These must be named "gd" and "gdd" respectively (gamma-dot, gamma-dot-dot). These are in units of m/s/day and m/s/day^2 respectively. :type params: dict :returns: The radial velocity of the star due to the trend (m/s). :rtype: float .. rubric:: Notes The radial velocity contribution from the trend is calculated as the sum of the linear and quadratic components: `gd*(t-t0)` and `gdd*((t-t0)**2)`. The constant offset (gamma) is now instrument-specific and handled via the Instrument class. Use `star.gamma_offsets(instrument)` to get per-instrument offsets. In general the trend is used to account for long-term effects such as acceleration from a distant companion. If you see a strong linear or quadratic trend in the data, it is worth investigating. .. rubric:: Examples >>> trend = Trend(t0=2458000.0, params={"gd": 0.001, "gdd": 0.0}) >>> rv_trend = trend.radial_velocity(times) .. py:attribute:: gammadot .. py:attribute:: gammadotdot .. py:method:: __str__() -> str .. py:method:: __repr__() -> str .. py:method:: _linear(t: numpy.ndarray, t0: float) -> numpy.ndarray .. py:method:: _quadratic(t: numpy.ndarray, t0: float) -> numpy.ndarray .. py:method:: radial_velocity(t: numpy.ndarray) -> numpy.ndarray Calculate radial velocity contribution from trend components. :param t: Time values :type t: array_like :returns: RV trend values (linear + quadratic terms) :rtype: array_like .. py:class:: Star(name: str, mass: float) Star with orbiting planet(s). :param name: The name of the star system. :type name: `str` :param mass: The mass of the star (solar masses) :type mass: `float` .. attribute:: planets The dict storing the `Planet` objects, the key is the `planet.letter` attribute. :type: `dict` .. attribute:: num_planets The number of `Planet` objects stored in the `Star` object. :type: int .. py:attribute:: name .. py:attribute:: mass .. py:attribute:: planets .. py:attribute:: instruments .. py:attribute:: num_planets :value: 0 .. py:method:: __repr__() -> str .. py:method:: __str__() -> str .. py:method:: add_planet(planet: Planet) -> None Store `Planet` object in `planets` dict with the key `Planet.letter`. Planets cannot share letters; if two planets have the same letter then the second one will overwrite the first. :param planet: A `ravest.model.Planet` object :type planet: `Planet` .. py:method:: add_trend(trend: Trend) -> None Store `Trend` object in `trend` attribute. :param trend: A `ravest.model.Trend` object :type trend: `Trend` .. py:method:: add_instrument(instrument: Instrument) -> None Store `Instrument` object in `instruments` dict with key `instrument.name`. Instruments cannot share names; if two instruments have the same name then the second one will overwrite the first. :param instrument: A `ravest.model.Instrument` object :type instrument: `Instrument` .. py:method:: gamma_offsets(instrument: numpy.ndarray) -> numpy.ndarray Return gamma offset for each observation based on instrument column. :param instrument: Array of instrument names for each observation. :type instrument: np.ndarray :returns: Array of gamma offsets (m/s) for each observation. :rtype: np.ndarray :raises ValueError: If an instrument in the data is not found in the Star's instruments. .. py:method:: jitter_values(instrument: numpy.ndarray) -> numpy.ndarray Return jitter value for each observation based on instrument column. :param instrument: Array of instrument names for each observation. :type instrument: np.ndarray :returns: Array of jitter values (m/s) for each observation. :rtype: np.ndarray .. py:method:: radial_velocity(t: numpy.ndarray) -> numpy.ndarray Calculate the radial velocity of the star at time ``t`` due to the planets and trend. Calculate the radial velocity of the star at time ``t`` due to the planets and the trend (constant velocity, linear, and quadratic.) This is a linear sum of the RVs of each of the `Planet` objects stored in the `Star`, and the RV of the `Trend` object. :param t: The time at which to calculate the radial velocity (day) :type t: `float` :returns: The radial velocity at time `t` (m/s) :rtype: `float` .. py:method:: mpsini(planet_letter: str, unit: str = 'kg') -> float Calculate the minimum mass of a planet in the star system. :param planet_letter: The letter of the planet to calculate the minimum mass of. :type planet_letter: `str` :param unit: The unit to return the planetary minimum mass in. Options are "kg", "M_earth", "M_jupiter". :type unit: `str` :returns: The minimum mass of the planet. :rtype: `float` .. py:method:: phase_plot(t: numpy.ndarray, ydata: numpy.ndarray, yerr: numpy.ndarray, instrument: numpy.ndarray) -> None Generate a phase plot for each planet. Given RV ``ydata`` at time ``t`` with errorbars ``yerr``, generates a phase plot for each planet. Data is coloured by instrument and gamma offsets are subtracted before plotting. :param t: The time of the observations ``ydata`` (day) :type t: np.ndarray :param ydata: The observed radial velocity at time ``t`` (m/s) :type ydata: np.ndarray :param yerr: The measurement error of the datapoint ``ydata`` (m/s) :type yerr: np.ndarray :param instrument: The instrument/telescope name for each observation. :type instrument: np.ndarray .. py:function:: calculate_mpsini(mass_star: float, period: float, semi_amplitude: float, eccentricity: float, unit: str = 'kg') -> float Calculate the minimum mass of the planet. :param mass_star: The mass of the star in solar masses. :type mass_star: `float` :param period: The orbital period of the planet in days. :type period: `float` :param semi_amplitude: The semi-amplitude of the radial velocity of the star in m/s. :type semi_amplitude: `float` :param eccentricity: The eccentricity of the orbit, 0 <= e < 1 (dimensionless). :type eccentricity: `float` :param unit: The unit to return the planetary minimum mass in. Options are "kg", "M_earth", "M_jupiter". :type unit: `str` :returns: The minimum mass mpsini of the planet. :rtype: `float` :raises ValueError: If the unit is not one of "kg", "M_earth", "M_jupiter". .. py:function:: fold_time_series(times: numpy.ndarray, period: float, t_ref: float) -> tuple[numpy.ndarray, numpy.ndarray] Fold time series to orbital phase and return sorted arrays. Converts times to orbital phases in the range [-0.5, 0.5] and returns both the sorted phases and the indices needed to sort other arrays consistently. :param times: Time values to fold :type times: np.ndarray :param period: Orbital period :type period: float :param t_ref: Reference time (usually Tc - time of transit/conjunction) :type t_ref: float :returns: * **phases_sorted** (*np.ndarray*) -- Phase-folded times in range [-0.5, 0.5], sorted ascending * **sort_indices** (*np.ndarray*) -- Indices that sort the original times by phase .. rubric:: Examples >>> times = np.array([0, 1, 2, 3, 4]) >>> period = 2.0 >>> t_ref = 0.5 # Reference time >>> phases, indices = fold_time_series(times, period, t_ref) >>> # phases will be sorted from -0.5 to 0.5 >>> # indices can be used to sort other arrays consistently