Added functions to reproduce Berti's plot

code/SubTasks/TwoBody/DataSets.jl

@@ -32,6 +32,8 @@ module DataSets
 
     sstars_data = (sstars_t, sstars_vals, sstars_errs)
 
-    θ₀ = [98947.158163, 161544.112937, -2.19582537e+02,  6.04500473e+02, 1.33902287e+02,  2.26126780e+02,  215932.83605999997, 29862.5963615, 6.71724897e+01,  5.87446040e+02,  1.70792316e+02,  9.34728478e+01, 139381.778738, -37591.510541999996, -3.45931227e+02,  1.00503663e+03, 1.53023909e+02,  3.23354746e+02, 4.31, 8.3, -3.49118611e-06,  1.47495749e-05, 1.0, 1.0]
+    #θ₀ = [98947.158163, 161544.112937, -2.19582537e+02,  6.04500473e+02, 1.33902287e+02,  2.26126780e+02,  215932.83605999997, 29862.5963615, 6.71724897e+01,  5.87446040e+02,  1.70792316e+02,  9.34728478e+01, 139381.778738, -37591.510541999996, -3.45931227e+02,  1.00503663e+03, 1.53023909e+02,  3.23354746e+02, 4.31, 8.3, -3.49118611e-06,  1.47495749e-05, 1.0, 1.0]
+
+    θ₀ = [101242.16764299569, 163115.15856316627, -209.32745199708657, 609.2274529754067, 135.5210562812287, 227.04515423639108, 19583.21818679856, 27162.349994682932, 57.93021035914549, 597.6757305313353, 176.60480011425918, 92.95428843303475, 139784.05807133473, -52028.52456219171, -247.73569673346708, 988.8216275239129, 148.80442935896926, 320.44376068418524, 4.506242786705499, 8.505059736146732, -5.0792651134477405e-6, 1.0045033094256844e-5, 1.1940393286488078, 0.751667430842911]
 
 end

code/SubTasks/TwoBody/ScalarTensor.jl

 __precompile__()
 module ScalarTensor
 
-    export ω_BD
+    export Θ, ω_BD_pl, dot_P, residual_2σ, residual_2σ_1, residual_2σ_2
 
+    # unit transformation constants
+    # for ex, to transform x AUs into planck lengths, one should multiply
+    # x by AU_to_r_planck:
     AU_to_r_planck = 9.25701e+45
     eV_to_m_planck = 1.78266e-36 / 2.17645e-8
+    M_sun_to_m_planck = 9.136366e+37
+    days_to_t_planck = 1.6026232e+48
+    km_to_r_planck = 6.187240576e+37
+    gradbyyr_to_radbyt_planck = π / 180.0 / 31556925.216 * 5.39116e-44
+
+    function Θ(x)
+    #Heavyside
+
+        if x ≤ 0.0
+            return 0.0
+        else
+            return 1.0
+        end
+
+    end
 
     function ω_BD_pl(mₛ_in, r_in, γ)   #mₛ in eV, r in AU
+    # This function computes ω_BD+1.5 for Shapiro time delay γ_PPN experiments
 
         mₛr = mₛ_in * eV_to_m_planck * r_in * AU_to_r_planck
-    
+
         return exp(-mₛr) / 2.0 * (1 + γ) / (1 - γ)
 
     end
 
+    function dot_P_GR_nopershift(pars, Δpars)
+    # This function computes Ṗ_GR with eq. 96
+
+        m₁_in, m₂_in, q, P_in = pars #with units as stated in paper
+        Δm₁_in, Δm₂_in, Δq, ΔP_in = Δpars #uncertanties of pars
+
+        P, ΔP = (P_in, ΔP_in) .* days_to_t_planck
+
+        m₁, Δm₁, m₂, Δm₂ = (m₁_in, Δm₁_in, m₂_in, Δm₂_in) .* M_sun_to_m_planck
+        m = m₁ + m₂
+        Δm = √(Δm₁^2 + Δm₂^2)
+
+        dot_P_GR = -192.0 * π / 5.0 * q / (1.0 + q)^2 * (2 * π * m / P)^(5.0/3.0) # eq. 96
+        ∂dot_P_GR_m = -64.0 * π * q / (1.0 + q)^2 * (2 * π / P)^(5.0/3.0) * m^(2.0/3.0) * Δm # That is ∂Ṗ_GR / ∂m * Δm
+        ∂dot_P_GR_q = -192.0 * π / 5.0 * (2 * π * m / P)^(5.0/3.0) * (1.0 - q) / (1.0 + q)^4 * Δq # That is ∂Ṗ_GR / ∂q * Δq
+        ∂dot_P_GR_P = 64 * π * q / (1.0 + q)^2 / P * (2 * π * m / P) ^ (5.0/3.0) * ΔP # That is ∂Ṗ_GR / ∂P * ΔP
+        Δdot_P_GR = √(∂dot_P_GR_m^2 + ∂dot_P_GR_q^2 + ∂dot_P_GR_P^2)
+
+        return dot_P_GR, Δdot_P_GR
+
+    end
+
+    function dot_P_GR_pershift(pars, Δpars)
+    # This function computes Ṗ_GR with eq. 98
+
+        q, P_in, dot_ω_in = pars
+        Δq, ΔP_in, Δdot_ω_in = Δpars
+
+        P, ΔP = (P_in, ΔP_in) .* days_to_t_planck
+        dot_ω, Δdot_ω = (dot_ω_in, Δdot_ω_in) .* gradbyyr_to_radbyt_planck
+
+        term1 = -4.0 / (1.0 + q)^2 * 8.0 / 15.0 / √(3.0) * (P / 2.0 / π)^1.5 * dot_ω^1.5
+        dot_P_GR = term1 * q * P * dot_ω
+        ∂dot_P_GR_q = term1 * (1.0 - q) / (1.0 + q) * P * dot_ω * Δq
+        ∂dot_P_GR_dot_ω = 2.5 * term1 * q * P * Δdot_ω
+        ∂dot_P_GR_P = 2.5 * term1 * q * dot_ω * ΔP
+        Δdot_P_GR = √(∂dot_P_GR_q^2 + ∂dot_P_GR_dot_ω^2 + ∂dot_P_GR_P^2)
+
+        return (dot_P_GR, Δdot_P_GR)
+
+    end
+
+    function dot_P(ω_BD, mₛ_in, pars, Δpars; s₁ = 0.2, s₂ = 1e-4)
+    # This function computes Ṗ with eq. 95
+
+        m₁_in, m₂_in, P_in, dot_P_GR = pars #with units as stated in paper
+        Δm₁_in, Δm₂_in, ΔP_in, Δdot_P_GR = Δpars #uncertanties of pars
+
+        m₁, Δm₁, m₂, Δm₂ = (m₁_in, Δm₁_in, m₂_in, Δm₂_in) .* M_sun_to_m_planck
+        m = m₁ + m₂
+        Δm = √(Δm₁^2 + Δm₂^2)
+
+        P, ΔP = (P_in, ΔP_in) .* days_to_t_planck
+        ω = 2 * π / P
+        Δω = 2 * π / P^2 * ΔP
+
+        mₛ = mₛ_in * eV_to_m_planck
+
+        r = (m * P^2 / 4.0 / π^2)^(1.0/3.0)
+        Δr = √(((P / (2 * π * m))^(2.0/3.0) / 3.0 * Δm)^2 + (2.0 / 3.0 * (m / (4 * π^2 * P))^(1.0/3.0) * ΔP)^2)
+
+        ξ = 1.0 / (2.0 + ω_BD)
+        Γ = 1.0 - 2.0 * (s₁ * m₂ + m₁ * s₂) / m
+        ∂Γ_m₁ = 2.0 * m₂ / m^2 * (s₂ - s₁) * Δm₁
+        ∂Γ_m₂ = 2.0 * m₁ / m^2 * (s₁ - s₂) * Δm₂
+        ΔΓ = √(∂Γ_m₁^2 + ∂Γ_m₂^2)
+
+        𝓖 = 1.0 - ξ * (s₁ + s₂ - 2 * s₁ * s₂)
+        term1 = (4 * ω^2 - mₛ^2) / (4 * ω^2)  # I defined this 2 terms
+        term2 = ξ * Γ * term1 * Θ(2 * ω - mₛ) # for ease of formulas
+
+        κ1 = 𝓖^2 * (12.0 - 6.0 * ξ + Γ * term1 * term2) # eq. 2
+        ∂κ1_ω = 2 * 𝓖^2 * Γ * term2 * mₛ^2 / ω^3 * Δω
+        ∂κ1_Γ = 2 * 𝓖^2 * term1 * term2 * ΔΓ
+        Δκ1 = √(∂κ1_ω^2 + ∂κ1_Γ^2)
+
+        κD = 2 * 𝓖 * ξ * (ω^2 - mₛ^2) / ω^2 * Θ(ω - mₛ) # eq. 2 too
+        ΔκD = abs(4 * 𝓖 * ξ * mₛ^2 / ω^3 * Θ(ω - mₛ) * Δω)
+
+        𝓢 = abs(s₁ - s₂)
+        term3 = (2 * π * m / P)^(-2.0/3.0)
+        Δterm3 = 2.0 / 3.0 * abs(term3) * √((Δm / m)^2 + (ΔP / P)^2)
+
+        ∂dot_P_dot_P_GR = (𝓖 ^ (-4.0/3.0) * κ1 / 12.0 + 5.0 / 48.0 * (2 * π * m / P)^(-2.0/3.0) * 𝓢^2 * κD)
+        dot_P = dot_P_GR * ∂dot_P_dot_P_GR # eq. 95
+        ∂dot_P_dot_P_GR = ∂dot_P_dot_P_GR * Δdot_P_GR
+        ∂dot_P_term3 = 5.0 / 48.0 * dot_P_GR * 𝓢^2 * κD * Δterm3
+        ∂dot_P_κ1 = dot_P_GR * 𝓖^(-4.0/3.0) / 12.0 * Δκ1
+        ∂dot_P_κD = dot_P_GR * 5.0 / 48.0 * term3 * 𝓢^2 * ΔκD
+        Δdot_P = √(∂dot_P_dot_P_GR^2 + ∂dot_P_term3^2 + ∂dot_P_κ1^2 + ∂dot_P_κD^2)
+        return (dot_P, Δdot_P)
+
+    end
+
+    function residual_2σ(ω_BD, mₛ_in, pars, Δpars, dot_P_obs, Δdot_P_obs; s₁ = 0.2, s₂ = 1e-4)
+    # this function computes |Ṗ_obs - Ṗ| - 2σ
+
+        dot_P_model, Δdot_P_model = dot_P(ω_BD, mₛ_in, pars, Δpars; s₁ = s₁, s₂ = s₂)
+        σ = √(Δdot_P_model^2 + Δdot_P_obs^2) # σ is combined uncertainty of Ṗ_obs and Ṗ_model
+        return abs(dot_P_obs - dot_P_model) - 2*σ
+
+    end
+
+    function get_residual_func(m₁_in, m₂_in, q, P_in, dot_ω_in, Δm₁_in, Δm₂_in, Δq, ΔP_in, Δdot_ω_in, dot_P_obs, Δdot_P_obs)
+
+        if dot_ω_in == 0.0
+            dot_P_GR, Δdot_P_GR = dot_P_GR_nopershift((m₁_in, m₂_in, q, P_in), (Δm₁_in, Δm₂_in, Δq, ΔP_in))
+        else
+            dot_P_GR, Δdot_P_GR = dot_P_GR_pershift((q, P_in, dot_ω_in), (Δq, ΔP_in, Δdot_ω_in))
+        end
+
+        return (ω_BD, mₛ_in) -> residual_2σ(ω_BD, mₛ_in, (m₁_in, m₂_in, P_in, dot_P_GR), (Δm₁_in, Δm₂_in, ΔP_in, Δdot_P_GR), dot_P_obs, Δdot_P_obs)
+
+    end
+
+    residual_2σ_1 = get_residual_func(1.27, 1.02, 1.245, 0.1976509593, 5.3096, 0.01, 0.01, 0.014, 1e-10, 4e-4, -4.01e-13, 0.25e-13)
+    residual_2σ_2 = get_residual_func(1.64, 0.16, 10.5, 0.60467271355, 0.0, 0.22, 0.02, 0.5, 3e-11, 0.0, -1.5e-14, 1.5e-14)
+
 end