Support for importing RadVel posteriors

OctofitterRadialVelocity/src/OctofitterRadialVelocity.jl

@@ -153,6 +153,7 @@ jd2mjd(jd) = jd + 2400000.5
 
 include("harps.jl")
 include("hires.jl")
+include("radvel.jl")
 
 
 # Plot recipe for astrometry data

OctofitterRadialVelocity/src/radvel.jl

@@ -0,0 +1,161 @@
+#=
+This file contains utilities for working
+with posteriors generated by RadVel
+=#
+using MCMCChains
+
+"""
+Given a Table / DataFrame containing a RadVel 
+posterior, convert it into Octofitter format.
+This allows it to be analyzed using all the 
+available tools of these packages.
+
+Only works with one planet.
+
+Supply priors for the stellar mass and parallax
+using Distributions.jl
+"""
+function radvel_posterior(
+    table,
+    M_prior,
+    plx_prior,
+
+)
+
+    # Calculations adapted from Sarah Blunt's
+    # https://github.com/sblunt/orbitize_radvel_utils/blob/main/compute_sep.py
+
+    tau_ref_epoch = 58849
+
+    semiamp = table.k1
+    period_days = table.per1
+    period_yrs = period_days/365.25
+    if hasproperty(table, :e1)
+        eccentricity = table.e1
+    else
+        eccentricity = fill(0.0, size(semiamp))
+    end
+
+    M_star = rand(M_prior, size(semiamp))
+
+    msini = Msini.(semiamp, period_days, M_star, eccentricity)
+
+    cosi = (2 .* rand(length(semiamp))) .- 1
+    inc = acos.(cosi)
+    m_pl = msini ./ sin.(inc)
+    mtot = M_star .+ m_pl*Octofitter.mjup2msol
+    sma = @. (period_yrs^2 * mtot)^(1/3)
+    if hasproperty(table, :e1)
+        omega_st_rad = table.w1
+    else
+        omega_st_rad = fill(0.0, size(semiamp))
+    end 
+
+    omega_pl_rad = omega_st_rad .+ π
+    parallax = rand(plx_prior, length(semiamp))
+
+    Ω = 2π .* rand(length(semiamp))
+    if hasproperty(table, :tp1)
+        tp_mjd = table.tp1 .- 2400000.5
+    else # have to convert from time of transit to time of to time of periastron 
+        tp_mjd = timetrans_to_timeperi.(table.tc1, period_days, eccentricity, omega_st_rad) .- 2400000.5
+    end
+    tau_ref_epoch = 58849
+    tau = @. (tp_mjd - tau_ref_epoch)/(period_days)
+    tau = mod.(tau, 1.0)
+
+    # [:loglike, :logpost, :tree_depth, :numerical_error, :M, :rv0_1, :rv0_2, :rv0_3, :rv0_4, :jitter_1, :jitter_2, :jitter_3, :jitter_4, :rv0_5, :jitter_5, :rv0_6, :jitter_6, :plx, :b_a, :b_τy, :b_τx, :b_mass, :b_i, :b_Ωy, :b_Ωx, :b_e, :b_ω, :b_τ, :b_Ω]
+
+    return MCMCChains.Chains(
+        # (data),
+        # (names)
+        hcat(
+            M_star,
+            parallax,
+            m_pl./Octofitter.mjup2msol,
+            sma,
+            inc,
+            Ω,
+            eccentricity,
+            omega_pl_rad,
+            tau,
+            msini
+        ),
+        [:M, :plx, :b_mass, :b_a, :b_i, :b_Ω, :b_e, :b_ω, :b_τ, :b_msini],
+    )
+    # return map(eachindex(semiamp)) do i
+    #     orbit(
+    #         a=sma[i],
+    #         M=M_star[i],
+    #         e=eccentricity[i],
+    #         τ=tau[i],
+    #         i=inc[i],
+    #         Ω=Ω[i],
+    #         ω=omega_pl_rad[i],
+    #         plx=parallax[i],
+    #     )
+    # end
+
+end
+
+
+
+# Adapted from 
+# https://radvel.readthedocs.io/en/latest/_modules/radvel/utils.html#Msini
+
+"""Calculate Msini
+
+Calculate Msini for a given K, P, stellar mass, and e
+
+Semi-amplitude K is in units of **earth masses**.
+Mstar is in solar mass units.
+
+Period must be in days.
+
+Returns Msini in *jupiter masses*
+"""
+function Msini(K, P, Mstar, e)
+
+    K_Msun =  K * 0.0031Octofitter.mjup2msol
+
+    # Normalization.
+    # RV m/s of a 1.0 Jupiter mass planet tugging on a 1.0
+    # solar mass star on a 1.0 year orbital period
+    K_0 = 28.4329
+
+
+    P_year = P / 365.25
+
+    # assumes that Mp << Mstar
+    Msini = K_Msun / K_0 * sqrt(1.0 - e^2) * Mstar^(2/3) * P_year^(1/3)
+
+    return Msini/Octofitter.mjup2msol
+end
+
+
+
+
+
+"""
+Convert Time of Transit to Time of Periastron Passage
+
+Args:
+    tc (float): time of transit
+    per (float): period [days]
+    ecc (float): eccentricity
+    omega (float): longitude of periastron (radians)
+
+Returns:
+    float: time of periastron passage
+
+"""
+function timetrans_to_timeperi(tc, per, ecc, omega)
+
+
+    f = π/2 - omega
+    ee = 2 * atan(tan(f/2) * sqrt((1-ecc)/(1+ecc)))  # eccentric anomaly
+    tp = tc - per/(2π) * (ee - ecc*sin(ee))      # time of periastron
+
+    return tp
+
+end