Exoplanet Detection: Radial Velocity Method

Astro 497: Week 4, Day 3

TableOfContents()

Logistics

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))$$

Markdown.parse("""
No closed form solution for \$E\$, so solve for \$E(t)\$ iteratively 
(see [code for `calc_ecc_anom(M, e)`]($(calc_ecc_anom_url))
""")

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

let
    ω_list = [0,π/4,π/2, 3π/4, π]
    ωstr_list = [L"0",L"\frac{\pi}{4}",L"\frac{\pi}{2}", L"\frac{3\pi}{4}", L"\pi"]
    e_list = [0.2, 0.4, 0.6, 0.8]
    plts = [ make_rv_vs_phase_panel(e,ω_list[ωi], panel_label = L"e=" * latexstring(e) * L"\;\; \omega=" * ωstr_list[ωi] ) for ωi in 1:length(ω_list),  e in e_list ]
    #annotate!(plts[3,4],0.5,-4,text(L"\mathrm{Time}",48))
    #annotate!(plts[1,3],-0.25,4,text(L"\mathrm{RV}",48,rotation=90))
    l = @layout [  b{0.04w} [ grid(4,5); c{0.03h} ]  ]
    plt_top = [ plot(link=:none,xlabel="", ylabel="",xticks=false,yticks=false,frame=:none) for e in e_list ]
    
    #annotate!(plt_top[1],0.5,0.5,text(L"\omega=0",36))
    plt_left = plot(link=:none,xlabel="", ylabel="",xticks=false,yticks=false,frame=:none)
    annotate!(plt_left,0.5,0.525,text(L"\mathrm{RV}",36,rotation=90,:center))
    plt_bottom = plot(link=:none,xlabel="", ylabel="",xticks=false,yticks=false,frame=:none)
    annotate!(plt_bottom,0.5,0.5,text(L"\mathrm{Time}",36,:center))
    plt = plot(plt_left,plts...,plt_bottom,  link=:none, layout = l, size=(2400,2400), thinkness_scaling=8)
end

function make_rv_vs_phase_panel(e, ω; P::Real=1, K::Real=1, M0::Real =0, panel_label="", t_max::Real = P, xticks=false, yticks=false )
    plt = plot(legend=:none, widen=false, xticks=xticks, yticks=yticks, margin=0mm, link=:none)
    t_plt = collect(range(0,stop=t_max,length=1000))
    rv_plt = calc_rv_keplerian.(t_plt, P, K, e, ω, M0)
    plot!(plt,t_plt,rv_plt, linecolor=:black) #, lineweight=4)
    xlims!(0,P)
    ylims!(-3.0,3.0)
    if length(panel_label)>0
        annotate!(plt, 0.505,2.1, text(panel_label,24,:center))
    end
    return plt
end

make_rv_vs_phase_panel (generic function with 1 method)

Interactive Keplerian RV Model

begin 
    plt_1pl = make_rv_vs_phase_panel(e_plt, ω_plt, P=P_plt, K=K_plt, M0=M0_plt, t_max = 100, xticks=true, yticks=true)
    xlabel!(plt_1pl, L"\mathrm{Time} \;\;\; (day)")
    ylabel!(plt_1pl, L"\Delta RV\;\;\;\; (m/s)")
    xlims!(plt_1pl,0,100)
    ylims!(plt_1pl,-6,6)
end

P: K: e: ω: M₀:

Pros & Cons

Strengths

Can detect planets with wide range of orbital inclinations
RV amplitude decreases slowly with orbital period
Provides information about planet mass
Effective tool for studying planets around sun-like stars
Can study stars at wide range of distances from Earth
Can characterize planetary systems with multiple planets
RV method has characterized hundreds of planetary systems

Weaknesses

Less sensitive to planets with nearly face-on orbits
Less sensitive to planets with large orbital periods
Characterizing long-period planets requires long-term observing programs
Mass information subject to uncertainties in $\sin i$ and $M_\star$.
Little Doppler information in spectra of stars with large $T_{\mathrm{eff}}$ or $v \sin i$
Applying RV to faint stars requires large telescopes
Characterizing planetary systems with low-mass planets and/or multiple planets can require hundreds of observations
Interpretting observations can be tricky:
- Unknown number of planets
- Potential for degeneracies between model parameters
Requires advanced high-resolution spectrographs & calibration system
Relatively few observatories have precision RV capabilities

Example RV Data

Ingest data

begin
    fn = joinpath("../_assets/week4/legacy_data.csv")
    if !isfile(fn) || !(filesize(fn)>0)
        path = joinpath(pwd(),"data")
        mkpath(path)
        fn = joinpath(path,"legacy_data.csv")
        fn = Downloads.download("https://github.com/leerosenthalj/CLSI/raw/master/legacy_tables/legacy_data.csv", fn)
    end
    if filesize(fn) > 0
        df_all = CSV.read(fn, DataFrame)
        select!(df_all,["name","jd","mnvel","errvel", "cts","sval","tel"])
        # Rename columns to match labels from table in original paper
        rename!(df_all, "name"=>"Name","jd"=>"d","mnvel"=>"RVel","errvel"=>"e_RVel","tel"=>"Inst", "sval"=>"SVal")	
        star_names = unique(df_all.Name)
        md"Read RVs for $(length(star_names)) stars from California Legacy Survey from https://github.com/leerosenthalj/CLSI & from [Rosenthal et al. (2021)](https://doi.org/10.3847/1538-4365/abe23c) into `df_all`."
    else
        df_all = DataFrame()
        star_names = String[""]

        danger(md"Error reading data file with RVs.  Expect empty plots below.")
    end
    
end

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

begin
    star_name = star_names[starid]
    df_star = df_all |> @filter( _.Name == star_name ) |> DataFrame
end;

begin 
    df_star_by_inst = DataFrame()
    try
    df_star_by_inst = df_star |> @groupby( _.Inst ) |> @map( {bjd = _.d, rv = _.RVel, σrv = _.e_RVel, inst= key(_), nobs_inst=length(_) }) |> DataFrame;
    catch
    end
end;

t_offset = 2450000;  # So labels on x-axis are more digestable

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

Star ID to plot: HD

starid = searchsortedfirst(star_names,star_name_to_plt);

let
    #upscale
    plt = plot() #legend=:none, widen=true)
    num_inst = size(df_star_by_inst,1)
    rvoffset = zeros(4) # [result2.minimizer[11], result2.minimizer[12], 0, 0]
    #rvoffset[1:2] .= result1.minimizer[6:7]
    #slope = result1.minimizer[8]
    #rvoffset[1:2] .= result2.minimizer[11:12]
    for inst in 1:num_inst
        lab = df_star_by_inst[inst,:inst]
        if lab == "lick" continue end
        if lab == "k" lab = "Keck (pre)" end
        if lab == "j" lab = "Keck (post)" end
        if lab == "apf" lab = "APF" end
        scatter!(df_star_by_inst[inst,:bjd].-t_offset,
                df_star_by_inst[inst,:rv].-rvoffset[inst],
                yerr=collect(df_star_by_inst[inst,:σrv]),
                label=lab)
                #markersize=4*upscale, legendfontsize=upscale*12
        #plot!(t_plt,model_1pl)
        #plot!(t_plt.-2450000,model_2pl)
        #scatter!(df_star_by_inst[2,:bjd].-2450000,df_star_by_inst[2,:rv], yerr=collect(df_star_by_inst[2,:σrv]),label=:none)
        #scatter!(df_star_by_inst[3,:bjd].-2450000,df_star_by_inst[3,:rv], yerr=collect(df_star_by_inst[3,:σrv]),label=:none)
        #scatter!(df_star_by_inst[4,:bjd].-2450000,df_star_by_inst[4,:rv], yerr=collect(df_star_by_inst[4,:σrv]),label=:none)
    end
    xlabel!("Time (d)")
    ylabel!("RV (m/s)")
    title!("HD " * star_name )
    #savefig(plt,joinpath(homedir(),"Downloads","RvEx.pdf"))
    plt
end

Fitting RV Model to data

if size(df_star_by_inst,1)>0  # Warning: Assume exactly 2 instruments providing RV data
    data1 = (t=collect(df_star_by_inst[1,:bjd]).-t_offset, rv=collect(df_star_by_inst[1,:rv]), σrv=collect(df_star_by_inst[1,:σrv]))
    if size(df_star_by_inst,1) > 1
        data2 = (t=collect(df_star_by_inst[2,:bjd]).-t_offset, rv=collect(df_star_by_inst[2,:rv]), σrv=collect(df_star_by_inst[2,:σrv]))
    else
        data2 = (t=Float64[], rv=Float64[], σrv=Float64[])
    end
    t_mean = (sum(data1.t)+sum(data2.t))/(length(data1.t).+length(data2.t))
    t_plt = range(minimum(vcat(data1.t,data2.t)), stop=maximum(vcat(data1.t,data2.t)), step=1.0)
end;

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

if try_fit_1pl && @isdefined data1 
    result1 = Optim.optimize(loss_1pl, θinit1, BFGS(), autodiff=:forward);
end

if @isdefined result1 
    result1.minimizer[1:5]
end

if try_fit_1pl && @isdefined data1
    pred_1pl = map(t->model_1pl(t,PKhkωpM_to_PKeωM(result1.minimizer[1:5])...,0.0,slope=0.0, t_mean=t_mean),t_plt);
end

if @isdefined result1
let
    #upscale
    plt = plot() #legend=:none, widen=true)
    num_inst = size(df_star_by_inst,1)
    rvoffset = zeros(4) # [result2.minimizer[11], result2.minimizer[12], 0, 0]
    rvoffset[1:2] .= result1.minimizer[6:7]
    slope = result1.minimizer[8]
    #rvoffset[1:2] .= result2.minimizer[11:12]
    for inst in 1:num_inst
        lab = df_star_by_inst[inst,:inst]
        if lab == "lick" continue end
        if lab == "k" lab = "Keck (pre)" end
        if lab == "j" lab = "Keck (post)" end
        if lab == "apf" lab = "APF" end
        scatter!(df_star_by_inst[inst,:bjd].-t_offset,
                df_star_by_inst[inst,:rv].-rvoffset[inst],
                yerr=collect(df_star_by_inst[inst,:σrv]),
                label=lab)
        #plot!(t_plt,model_1pl)
        #plot!(t_plt.-2450000,model_2pl)
        #scatter!(df_star_by_inst[2,:bjd].-2450000,df_star_by_inst[2,:rv], yerr=collect(df_star_by_inst[2,:σrv]),label=:none)
        #scatter!(df_star_by_inst[3,:bjd].-2450000,df_star_by_inst[3,:rv], yerr=collect(df_star_by_inst[3,:σrv]),label=:none)
        #scatter!(df_star_by_inst[4,:bjd].-2450000,df_star_by_inst[4,:rv], yerr=collect(df_star_by_inst[4,:σrv]),label=:none)
    end
    plot!(t_plt,pred_1pl, label=:none)
    xlabel!("Time (d)")
    ylabel!("RV (m/s)")
    title!("HD " * star_name )
    #savefig(plt,joinpath(homedir(),"Downloads","RvEx.pdf"))
    plt
end
end

if @isdefined result1
    resid1 = vcat(
     data1.rv .- model_1pl.(data1.t,PKhkωpM_to_PKeωM(result1.minimizer[1:5])...,result1.minimizer[6],slope=result1.minimizer[8], t_mean=t_mean),
    data2.rv .- model_1pl.(data2.t,PKhkωpM_to_PKeωM(result1.minimizer[1:5])...,result1.minimizer[7],slope=result1.minimizer[8], t_mean=t_mean) )
end;

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

if try_fit_2nd_pl_to_resid && (@isdefined data1) && (@isdefined resid1)
result_resid = Optim.optimize(loss_1pl_resid, θinit_resid, BFGS(),autodiff=:forward ) ;
end

if @isdefined result_resid
    model_resid = map(t->calc_rv_keplerian(t.-t_mean,PKhkωpM_to_PKeωM(result_resid.minimizer[1:5])...),t_plt);
end;

if @isdefined result_resid
let
    #upscale
    plt = plot(legend=:none, widen=true)
    num_inst = size(df_star_by_inst,1)
    rvoffset = result1.minimizer[6:7]
    slope = result1.minimizer[8]

    scatter!(data1.t,
                data1.rv.-
                model_1pl.(data1.t,PKhkωpM_to_PKeωM(result1.minimizer[1:5])...,rvoffset[1],slope=result1.minimizer[8], t_mean=t_mean),
                yerr=data1.σrv) #,
                #markersize=4*upscale, legendfontsize=upscale*12)
    scatter!(data2.t,
                data2.rv.-
                model_1pl.(data2.t,PKhkωpM_to_PKeωM(result1.minimizer[1:5])...,rvoffset[2],slope=result1.minimizer[8], t_mean=t_mean),
                yerr=data2.σrv),
                #markersize=4*upscale, legendfontsize=upscale*12)
    plot!(t_plt, model_resid)
    xlabel!("Time (d)")
    ylabel!("RV (m/s)")
    title!("HD " * star_name * " (residuals to 1 planet model)")
    #savefig(plt,joinpath(homedir(),"Downloads","RvEx.pdf"))
    plt
end
end

Fit 2-planet model

Fit 2nd planet model:

if try_fit_2pl && (@isdefined result1) && (@isdefined result_resid)
    θinit2 = [result1.minimizer[1:5]..., result_resid.minimizer[1:5]..., result1.minimizer[6]+result_resid.minimizer[6],result1.minimizer[7]+result_resid.minimizer[7], result1.minimizer[8]+result_resid.minimizer[8], result_resid.minimizer[9]];
    result2 = Optim.optimize(loss_2pl, θinit2, BFGS(), autodiff=:forward)
end

if (@isdefined result1) && (@isdefined result2)
    pred_2pl = map(t->model_2pl(t,PKhkωpM_to_PKeωM(result2.minimizer[1:5])...,PKhkωpM_to_PKeωM(result2.minimizer[6:10])...,0.0,slope=result2.minimizer[13], t_mean=t_mean),t_plt)
    resid2 = vcat(
     data1.rv .- model_2pl.(data1.t,PKhkωpM_to_PKeωM(result2.minimizer[1:5])...,PKhkωpM_to_PKeωM(result2.minimizer[6:10])...,result2.minimizer[11],slope=result2.minimizer[13], t_mean=t_mean),
    data2.rv .- model_2pl.(data2.t,PKhkωpM_to_PKeωM(result2.minimizer[1:5])...,PKhkωpM_to_PKeωM(result2.minimizer[6:10])...,result2.minimizer[12],slope=result2.minimizer[13], t_mean=t_mean) )
end;

if @isdefined result2
    #upscale
    plt_fit = plot(widen=true, xticks=false)
    num_inst = size(df_star_by_inst,1)
    rvoffset = result2.minimizer[11:12]
    slope = result2.minimizer[13]
    #rvoffset[1:2] .= result2.minimizer[11:12]
    for inst in 1:num_inst
        lab = df_star_by_inst[inst,:inst]
        if lab == "lick" continue end
        if lab == "k" lab = "Keck (pre)" end
        if lab == "j" lab = "Keck (post)" end
        if lab == "apf" lab = "APF" end
        scatter!(df_star_by_inst[inst,:bjd].-t_offset,
                df_star_by_inst[inst,:rv].-rvoffset[inst],
                yerr=collect(df_star_by_inst[inst,:σrv]),
                label=lab)
                #markersize=3*upscale, legendfontsize=upscale*12,
        #plot!(t_plt,model_1pl)
        #plot!(t_plt.-2450000,model_2pl)
        #scatter!(df_star_by_inst[2,:bjd].-2450000,df_star_by_inst[2,:rv], yerr=collect(df_star_by_inst[2,:σrv]),label=:none)
        #scatter!(df_star_by_inst[3,:bjd].-2450000,df_star_by_inst[3,:rv], yerr=collect(df_star_by_inst[3,:σrv]),label=:none)
        #scatter!(df_star_by_inst[4,:bjd].-2450000,df_star_by_inst[4,:rv], yerr=collect(df_star_by_inst[4,:σrv]),label=:none)
    end
    plot!(t_plt,pred_2pl,linecolor=:black, label=:none)
    xlabel!("Time (d)")
    ylabel!("RV (m/s)")
    #title!("HD " * star_name )
    #savefig(plt_fit,joinpath(homedir(),"Downloads","RvEx.pdf"))
    plt_fit
end;

if @isdefined resid2
    plt_resid = plot(legend=:none)
    
    scatter!(data1.t, resid2[1:length(data1.t)], yerr=data1.σrv)
    scatter!(data2.t, resid2[length(data1.t)+1:end], yerr=data2.σrv)
    plot!(t_plt,t_plt.*0, linestyle=:dot, linecolor=:black)
    xlabel!("Time (d)")
    ylabel!("ΔRV (m/s)")
end;

if @isdefined plt_fit 
    l = @layout [a{0.7h} ; b ]
    plt_combo = 
    plot(plt_fit, plt_resid, layout = l)
    plt_combo
end

Setup

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

notebook_id = PlutoRunner.notebook_id[] |> string

"6b35a506-5e66-11ed-3176-fb3a73788cb8"

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.

function loss_1pl(θ) 
    (P1, K1, h1, k1, Mpω1, C1, C2, slope, σj ) = θ
    ( P1, K1, e1, ω1, M1 ) = PKhkωpM_to_PKeωM([P1, K1, h1, k1, Mpω1])
    if e1>1 return 1e6*e1 end
    rv_model1 = model_1pl.(data1.t,P1,K1,e1,ω1,M1,C1, slope=slope, t_mean=t_mean)
    loss = 0.5*sum(((rv_model1.-data1.rv)./(data1.σrv.+σj^2)).^2)
    rv_model2 = model_1pl.(data2.t,P1,K1,e1,ω1,M1,C2, slope=slope, t_mean=t_mean)
    loss += 0.5*sum((rv_model2.-data2.rv).^2 ./(data2.σrv.^2 .+σj^2))
    loss += 0.5*sum(log.(2π*(data1.σrv.^2 .+σj^2)))
    loss += 0.5*sum(log.(2π*(data2.σrv.^2 .+σj^2)))
    return loss
end

loss_1pl (generic function with 1 method)

function loss_1pl_resid(θ) 
    (P1, K1, h1, k1, Mpω1, C1, C2, slope, σj  ) = θ
    ( P1, K1, e1, ω1, M1 ) = PKhkωpM_to_PKeωM([P1, K1, h1, k1, Mpω1])
    if e1>1 return 1e6*e1 end
    rv_model1 = model_1pl.(data1.t,P1,K1,e1,ω1,M1,C1, slope=slope, t_mean=t_mean)
    loss = 0.5*sum(((rv_model1.-resid1[1:length(data1.t)])./(data1.σrv.+σj^2)).^2)
    rv_model2 = model_1pl.(data2.t,P1,K1,e1,ω1,M1,C2, slope=slope, t_mean=t_mean)
    loss += 0.5*sum(((rv_model2.-resid1[length(data1.t)+1:end])./(data2.σrv.+σj^2)).^2)
    loss += 0.5*sum(log.(2π*(data1.σrv.^2 .+σj^2)))
    loss += 0.5*sum(log.(2π*(data2.σrv.^2 .+σj^2)))
    return loss
end

loss_1pl_resid (generic function with 1 method)

function loss_2pl(θ) 
    (P1, K1, h1, k1, Mpω1, P2, K2, h2, k2, Mpω2, C1, C2, slope, σj ) = θ
    (P1, K1, e1, ω1, M1 ) = PKhkωpM_to_PKeωM([P1, K1, h1, k1, Mpω1])
    (P2, K2, e2, ω2, M2 ) = PKhkωpM_to_PKeωM([P2, K2, h2, k2, Mpω2])
    if e1>1 return 1e6*e1 end
    if e2>1 return 1e6*e2 end
    rv_model1 = model_2pl.(data1.t,P1,K1,e1,ω1,M1,P2,K2,e2,ω2,M2,C1, slope=slope, t_mean=t_mean)
    loss = 0.5*sum((rv_model1.-data1.rv).^2 ./(data1.σrv.^2 .+σj^2))
    loss += 0.5*sum(log.(2π*(data1.σrv.^2 .+σj^2)))
    rv_model2 = model_2pl.(data2.t,P1,K1,e1,ω1,M1,P2,K2,e2,ω2,M2,C2, slope=slope, t_mean=t_mean)
    loss += 0.5*sum((rv_model2.-data2.rv).^2 ./(data2.σrv.^2 .+σj^2))
    loss += 0.5*sum(log.(2π*(data2.σrv.^2 .+σj^2)))
    return loss
end

loss_2pl (generic function with 1 method)

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.)

#=begin
    fn = joinpath("../_assets/week4/","apjsabe23ct6_mrt.txt");
    #=
    if !isfile(fn) || !(filesize(fn)>0)
        fn = Downloads.download("https://psuastro497.github.io/Fall2022/assets/week4/apjsabe23ct6_mrt.txt", fn)
    end
    =#
    if filesize(fn)>0
        df_all = read_apj_mrt(fn)
        star_names = unique(df_all.Name)
        md"Read machine readable version of Table 6 from [Rosenthal et al. (2021)](https://doi.org/10.3847/1538-4365/abe23c) into `df_all`."
    else 
        df_all = DataFrame()
        star_names = String[""]
        danger(md"Error reading data file with RVs.  Expect empty plots below.")
    end
end
=#

function read_apj_mrt(fn::AbstractString)
    lines = readlines(fn)
    line_start_fmt_specs = findfirst(l->occursin("Bytes",l),lines)+2
    line_stop_fmt_specs = findfirst(l->occursin(r"^\-+$",l),lines[line_start_fmt_specs:end]) + (line_start_fmt_specs-1) -1
    line_data_start = findlast(l->occursin(r"^\-+$",l),lines) +1
    fmt_info = map(l->match(r"^\s*(\d+)-\s+(\d+)\s+(\w)(\d\.?\d*)\s*(\S+)\s+(\S+)\s+(.*)$",l).captures,lines[line_start_fmt_specs:line_stop_fmt_specs]) 
    colnames = map(f->f[6],fmt_info)
    data = map(fmt->Base.Fix2(extract_entry,fmt).(lines[line_data_start:length(lines)]),fmt_info)
    df = DataFrame(data,colnames )
end

read_apj_mrt (generic function with 1 method)

function extract_entry(line::AbstractString, fmt_info)
    substr = line[parse(Int,fmt_info[1]):parse(Int,fmt_info[2])]
    if occursin("--",substr)
        return missing
    end
    if  fmt_info[3] == "A" 
        return strip(substr)
    end
    if !occursin(r"\d",substr) 
        return missing
    end
    if fmt_info[3] == "I" 
        return parse(Int,substr)
    elseif  fmt_info[3] == "F" 
        return parse(Float64,substr)
    else
        @warn "Need to add instructions for parsing type " fmt_info[1][3]
        return nothing
    end
end

extract_entry (generic function with 1 method)

Built with Julia 1.8.2 and

CSV 0.10.3
DataFrames 1.3.2
Downloads 1.6.0
LaTeXStrings 1.3.0
Optim 1.7.2
ParameterHandling 0.4.5
Plots 1.26.0
PlutoTeachingTools 0.1.7
PlutoUI 0.7.37
Query 1.0.0

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