Improve this site

Intro to Analyzing RV Timeseries

Astro 589: Week 4

Circular RV Model

Velocity of Planet:

$$v_{pl} = 2\pi a_{pl} / P$$

$$a = a_{pl} + a_\star$$

Velocity of Star:

$$M_\star v_\star = m_{pl} v_{pl}$$

$$v_\star = \frac{m_{pl}}{M_\star+m_{pl}} \frac{2\pi a}{P}$$

Doppler effective is only sensitive to motion projected onto observer's line of sight

$$RV_{\star} = \frac{m_{pl}}{M_\star+m_{pl}} \frac{2\pi a}{P} \sin i$$

Keplerian Orbit

Credit: Lasunncty via Wikipedia (CC BY-SA 3.0 license)

Keplerian RV Model

$$\Delta RV(t) = \frac{K}{\sqrt{1-e^2}} \left[\cos(\omega+ν(t)) + e \cos(\omega) \right]$$

$$\begin{eqnarray} K & = & \frac{2\pi a \sin i}{P} \frac{m_{pl}}{M_\star+m_{pl}} \\ % & = & \frac{2\pi a \sin i}{P} \\ \end{eqnarray}$$

$$\tan\left(\frac{\nu(t)}{2}\right) = \sqrt{\frac{1+e}{1-e}} \tan\left(\frac{E(t)}{2}\right)$$

  • Mean anomaly (\(M(t)\)) increases linear with time

  • Eccentric anomaly (\(E(t)\)) specifies position in orbit using angle from center of elipse

  • True anomaly (\(\nu(t)\) or \(f(t)\) or \(T(t)\)) specifies position in orbit using angle from focus of ellipse

Credit: CheCheDaWaff via Wikipedia (CC BY-SA 4.0 license)

Keplerian Orbit in Time

$$M(t) = \frac{2π t}{P} + M_0$$

Kepler Equation

$$M(t) = E(t) - e \sin(E(t))$$

No closed form solution for \(E\), so solve for \(E(t)\) iteratively (see code for calc_ecc_anom(M, e)

make_rv_vs_phase_panel (generic function with 1 method)

Interactive Keplerian RV Model

P: K: e: ω: M₀:

Example RV Data

# Examples to show: 114783, 119850, 13931, 187123

Star ID to plot: HD

Fitting RV Model to data

Fit one planet

Fit 1-planet model: 

θinit1 = [496.9, 10, 0.1, 0.1, 0.1, 1, 2, 1e-4, 2.0] # for 114783
9-element Vector{Float64}:
 496.9
  10.0
   0.1
   0.1
   0.1
   1.0
   2.0
   0.0001
   2.0

Fit 2nd planet to residuals

Fit 2nd planet to residuals: 

θinit_resid = [4319.0, 5.0, 0.1, 0.1, π/4, 1.0, -1.0 , 0.0001, 1.0];  # for 114783

Fit 2-planet model

Fit 2nd planet model:

Setup

begin
    using CSV, DataFrames, Query
    using Optim
    using Plots, LaTeXStrings, Plots.Measures
    using PlutoUI, PlutoTeachingTools
    using Downloads
    using ParameterHandling

    default(size=(1200,800))
end
     

Keplerian Radial Velocity Model Code

begin 
    """ Calculate RV from t, P, K, e, ω and M0	"""
    function calc_rv_keplerian end 
    #calc_rv(t, p::Vector) = calc_rv(t, p[1],p[2],p[3],p[4],p[5])
    calc_rv_keplerian(t, p::Vector) = calc_rv_keplerian(t, p...)
    function calc_rv_keplerian(t, P,K,e,ω,M0) 
        mean_anom = t*2π/P-M0
        ecc_anom = calc_ecc_anom(mean_anom,e)
        true_anom = calc_true_anom(ecc_anom,e)
        rv = K/sqrt((1-e)*(1+e))*(cos(ω+true_anom)+e*cos(ω))
    end
end
calc_rv_keplerian (generic function with 2 methods)
function calc_true_anom(ecc_anom::Real, e::Real)
    true_anom = 2*atan(sqrt((1+e)/(1-e))*tan(ecc_anom/2))
end
calc_true_anom (generic function with 1 method)
begin
    calc_ecc_anom_cell_id = PlutoRunner.currently_running_cell_id[] |> string
    calc_ecc_anom_url = "#$(calc_ecc_anom_cell_id)"
    """
       calc_ecc_anom( mean_anomaly, eccentricity )
       calc_ecc_anom( param::Vector )
    
    Estimates eccentric anomaly for given 'mean_anomaly' and 'eccentricity'.
    If passed a parameter vector, param[1] = mean_anomaly and param[2] = eccentricity. 
    
    Optional parameter `tol` specifies tolerance (default 1e-8)
    """
    function calc_ecc_anom end
    
    function calc_ecc_anom(mean_anom::Real, ecc::Real; tol::Real = 1.0e-8)
      	if !(0 <= ecc <= 1.0)
            println("mean_anom = ",mean_anom,"  ecc = ",ecc)
        end
        @assert 0 <= ecc <= 1.0
        @assert 1e-16 <= tol < 1
      	M = rem2pi(mean_anom,RoundNearest)
        E = ecc_anom_init_guess_danby(M,ecc)
        local E_old
        max_its_laguerre = 200
        for i in 1:max_its_laguerre
           E_old = E
           E = update_ecc_anom_laguerre(E_old, M, ecc)
           if abs(E-E_old) < tol break end
        end
        return E
    end
    
    function calc_ecc_anom(param::Vector; tol::Real = 1.0e-8)
        @assert length(param) == 2
        calc_ecc_anom(param[1], param[2], tol=tol)
    end;
end
calc_ecc_anom (generic function with 2 methods)
"""
   ecc_anom_init_guess_danby(mean_anomaly, eccentricity)

Returns initial guess for the eccentric anomaly for use by itterative solvers of Kepler's equation for bound orbits.  

Based on "The Solution of Kepler's Equations - Part Three"
Danby, J. M. A. (1987) Journal: Celestial Mechanics, Volume 40, Issue 3-4, pp. 303-312 (1987CeMec..40..303D)
"""
function ecc_anom_init_guess_danby(M::Real, ecc::Real)
    @assert -2π<= M <= 2π
    @assert 0 <= ecc <= 1.0
    if  M < zero(M)
        M += 2π
    end
    E = (M<π) ? M + 0.85*ecc : M - 0.85*ecc
end;
"""
   update_ecc_anom_laguerre(eccentric_anomaly_guess, mean_anomaly, eccentricity)

Update the current guess for solution to Kepler's equation
  
Based on "An Improved Algorithm due to Laguerre for the Solution of Kepler's Equation"
   Conway, B. A.  (1986) Celestial Mechanics, Volume 39, Issue 2, pp.199-211 (1986CeMec..39..199C)
"""
function update_ecc_anom_laguerre(E::Real, M::Real, ecc::Real)
  #es = ecc*sin(E)
  #ec = ecc*cos(E)
  (es, ec) = ecc .* sincos(E)  # Does combining them provide any speed benefit?
  F = (E-es)-M
  Fp = one(M)-ec
  Fpp = es
  n = 5
  root = sqrt(abs((n-1)*((n-1)*Fp*Fp-n*F*Fpp)))
  denom = Fp>zero(E) ? Fp+root : Fp-root
  return E-n*F/denom
end;

Fitting Keplerian RV model

begin 
    """ Calculate RV from t, P, K, e, ω, M0	and C"""
    function calc_rv_keplerian_plus_const end 
    calc_rv_keplerian_plus_const(t, p::Vector) = calc_rv_keplerian_plus_const(t, p...)
    
    function calc_rv_keplerian_plus_const(t, P,K,e,ω,M0,C) 
        calc_rv_keplerian(t, P,K,e,ω,M0) + C
    end
end
calc_rv_keplerian_plus_const (generic function with 2 methods)
begin
    """ Calculate RV from t, P1, K1, e1, ω1, M01, P2, K2, e2, ω2, M02 and C"""
    calc_rv_2keplerians_plus_const(t, p::Vector) = calc_rv_2keplerians_plus_const(t, p...)
    function calc_rv_2keplerians_plus_const(t, P1,K1,e1,ω1,M01,P2,K2,e2,ω2,M02,C) 
        calc_rv_keplerian(t, P1,K1,e1,ω1,M01) + calc_rv_keplerian(t, P2,K2,e2,ω2,M02) + C
    end
end
calc_rv_2keplerians_plus_const (generic function with 2 methods)
""" Calculate RV from t, P, K, e, ω, M0	and C with optional slope and t_mean"""
function model_1pl(t, P, K, e, ω, M, C; slope=0.0, t_mean = 0.0)
    calc_rv_keplerian(t-t_mean,P,K,e,ω,M) + C + slope * (t-t_mean)
end
""" Calculate RV from t, P1, K1, e1, ω1, M01, P2, K2, e2, ω2, M02 and C with optional slope and t_mean"""
function model_2pl(t, P1, K1, e1, ω1, M1, P2, K2, e2, ω2, M2, C; slope=0.0, t_mean = 0.0)
    rv = calc_rv_keplerian(t-t_mean,P1,K1,e1,ω1,M1) + 
         calc_rv_keplerian(t-t_mean,P2,K2,e2,ω2,M2) + 
         C + slope * (t-t_mean)
end
""" Convert vector of (P,K,h,k,ω+M0) to vector of (P, K, e, ω, M0) """
function PKhkωpM_to_PKeωM(x::Vector) 
    (P, K, h, k, ωpM) = x
    ω = atan(h,k)
    return [P, K, sqrt(h^2+k^2), ω, ωpM-ω]
end

Loss functions

The functions below:

  • assume observations from exactly two instruments.

  • make use of the global variables data1, data2 and t_mean.

loss_1pl (generic function with 1 method)
loss_1pl_resid (generic function with 1 method)
loss_2pl (generic function with 1 method)

Ingest data

Read RVs for 719 stars from California Legacy Survey from https://github.com/leerosenthalj/CLSI & from Rosenthal et al. (2021) into df_all.

Code for parsing machine readable tables from AAS Journals

(Not actually used here, since found CSV version of files at https://github.com/leerosenthalj/CLSI. But potentially useful for student projects.)

read_apj_mrt (generic function with 1 method)
extract_entry (generic function with 1 method)

Built with Julia 1.11.1 and

CSV 0.10.15
DataFrames 1.3.6
Downloads 1.6.0
LaTeXStrings 1.3.1
Optim 1.7.8
ParameterHandling 0.4.10
Plots 1.26.1
PlutoTeachingTools 0.1.7
PlutoUI 0.7.60
Query 1.0.0

To run this notebook locally, download this file and open it with Pluto.jl.