Commit 04e1013a authored by Rustam Gainutdinov's avatar Rustam Gainutdinov
Browse files

Added two body problem modelling code

parent 04594aa4
__precompile__()
module EoMs
# This module contains the Equations of motion
# of PPN two-body problem
# Note: I assume G = c = 1
export EoM_massless
function EoM_massless(u, p, t)
x, y, vx, vy = u
M, β, γ = p
r = (x^2 + y^2)
φ = -M / r
= vx^2 + vy^2
rv = x * vx + y * vy
A = - (1 + 2 * (β + γ) * φ + γ * ) * M / r^3
B = M * (2 * γ + 2) * rv / r^3
acc_x = A * x + B * vx
acc_y = A * y + B * vy
return [vx, vy, acc_x, acc_y]
end
function EoM_general(u, p, t)
x1, y1, vx1, vy1, x2, y2, vx2, vy2 = u
M1, M2, β, γ = p
end
#S2 params: 98947.158163, 161544.112937, -0.0007324487518145469, 0.0020163968040454333, 4.31, 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, 0.0, 0.0, -3.49118611e-06, 1.47495749e-05, 1.0, 1.0]
end
__precompile__()
module Models
# This module stands for modelling the observational data
include("EoMs.jl")
using DifferentialEquations
using Roots
using .EoMs
export model_massless
rad_to_arcsec = 3600.0 * 1000.0 * 180.0 / π
c_to_kms = 299792.458
r̂_to_km = 1476625.03852
pc_to_km = 3600.0 * 180.0 / π * 149697870.7
r̂_to_kpc = r̂_to_km / pc_to_km / 1000.0
t̂_to_sec = 4.92549094921
yr_to_sec = 365.24219878 * 86400.0
t̂_to_yr = t̂_to_sec / yr_to_sec
t_beg = 1992.224 / t̂_to_yr
t_end = 2019.357 / t̂_to_yr
function turnorbit(x, y, ω_gr, Ω_gr, i_gr)
N = size(x)
r0 = [x, y, 0.0]
ω = ω_gr / 180.0 * π
Ω = Ω_gr / 180.0 * π + π / 2.0
i = i_gr / 180.0 * π
= cos(ω)
= sin(ω)
= cos(Ω)
= sin(Ω)
ci = cos(i)
si = sin(i)
Ω_turn = [ - 0.0; 0.0; 0.0 0.0 1.0]
i_turn = [1.0 0.0 0.0; 0.0 ci -si; 0.0 si ci]
ω_turn = [ - 0.0; 0.0; 0.0 0.0 1.0]
r = Ω_turn * i_turn * ω_turn * r0
return r
end
function solve_EoM_massless(x0, y0, vx0, vy0, M, β, γ; tol=(1e-7, 1e-10))
p = [M, β, γ]
u0 = [x0, y0, vx0, vy0]
tspan = (t_beg, t_end)
prob = ODEProblem(EoM_massless, u0, tspan, p, reltol=tol[1], abstol=tol[2])
sol = solve(prob)
return sol
end
function Rømer_delay(sol, Ω, i, t_obs)
z = t -> turnorbit(sol(t)[1], sol(t)[2], 0.0, Ω, i)
f = t_em_arg -> t_em_arg - z(t_em_arg) - t_obs
t_em = find_zero(f, t_obs)
return t_em
end
function model_massless(θ, 𝐭_input)
S2_orb, S38_orb, S55_orb, M, R0, x_SgrA, y_SgrA, vx_SgrA, vy_SgrA, β, γ = θ
S2_x0, S2_y0, S2_vx0, S2_vy0, S2_i, S2_Ω = S2_orb
S38_x0, S38_y0, S38_vx0, S38_vy0, S38_i, S38_Ω = S38_orb
S55_x0, S55_y0, S55_vx0, S55_vy0, S55_i, S55_Ω = S55_orb
S2_sol = solve_EoM_massless(S2_x0, S2_y0, S2_vx0 / 299792.458, S2_vy0 / 299792.458, M, β, γ)
S38_sol = solve_EoM_massless(S38_x0, S38_y0, S38_vx0 / 299792.458, S38_vy0 / 299792.458, M, β, γ)
S55_sol = solve_EoM_massless(S55_x0, S55_y0, S55_vx0 / 299792.458, S55_vy0 / 299792.458, M, β, γ)
𝐭_em_S2 = Rømer_delay.(Ref(S2_sol), Ref(S2_Ω), Ref(S2_i), 𝐭_input)
𝐭_em_S38 = Rømer_delay.(Ref(S38_sol), Ref(S38_Ω), Ref(S38_i), 𝐭_input)
𝐭_em_S55 = Rømer_delay.(Ref(S55_sol), Ref(S55_Ω), Ref(S55_i), 𝐭_input)
S2_𝐱, S2_𝐲, S2_𝐯𝐱, S2_𝐯𝐲 = S2_sol.(𝐭_em_S2)
S38_𝐱, S38_𝐲, S38_𝐯𝐱, S38_𝐯𝐲 = S38_sol.(𝐭_em_S38)
S55_𝐱, S55_𝐲, S55_𝐯𝐱, S55_𝐯𝐲 = S55_sol.(𝐭_em_S55)
S2_𝐫 = turnorbit.(S2_𝐱, S2_𝐲, Ref(0.0), Ref(S2_Ω), Ref(S2_i))
S2_𝐯 = turnorbit.(S2_𝐯𝐱, S2_𝐯𝐲, Ref(0.0), Ref(S2_Ω), Ref(S2_i))
S38_𝐫 = turnorbit.(S38_𝐱, S38_𝐲, Ref(0.0), Ref(S38_Ω), Ref(S38_i))
S38_𝐯 = turnorbit.(S38_𝐯𝐱, S38_𝐯𝐲, Ref(0.0), Ref(S38_Ω), Ref(S38_i))
S55_𝐫 = turnorbit.(S55_𝐱, S55_𝐲, Ref(0.0), Ref(S55_Ω), Ref(S55_i))
S55_𝐯 = turnorbit.(S55_𝐯𝐱, S55_𝐯𝐲, Ref(0.0), Ref(S55_Ω), Ref(S55_i))
Dist = R0 / r̂_to_kpc
S2_RA = -((S2_𝐫[:,1] + vx_SgrA*(𝐭_em_S2 .- t_beg)) ./ Dist) .* rad_to_arcsec
S2_Dec = ((S2_𝐫[:,2] + vy_SgrA*(𝐭_em_S2 .- t_beg)) ./ Dist) .* rad_to_arcsec
S2_RV = (-S2_𝐯[:,3] .- M ./ .(S2_𝐱.^2 + S2_𝐲.^2) .+ (S2_𝐯𝐱.^2 .+ S2_𝐯𝐲.^2) ./ 2.0) .* c_to_kms
S38_RA = -((S38_𝐫[:,1] + vx_SgrA*(𝐭_em_S38 .- t_beg)) ./ Dist) .* rad_to_arcsec
S38_Dec = ((S38_𝐫[:,2] + vy_SgrA*(𝐭_em_S38 .- t_beg)) ./ Dist) .* rad_to_arcsec
S38_RV = (-S38_𝐯[:,3] .- M / .(S38_𝐱.^2 + S38_𝐲.^2) .+ (S38_𝐯𝐱.^2 .+ S38_𝐯𝐲.^2) ./ 2.0) .* c_to_kms
S55_RA = -((S55_𝐫[:, 1] + vx_SgrA*(𝐭_em_S55 .- t_beg)) ./ Dist) .* rad_to_arcsec
S55_Dec = ((S55_𝐫[:, 2] + vy_SgrA*(𝐭_em_S55 .- t_beg)) ./ Dist) .* rad_to_arcsec
S55_RV = (-S55_𝐯[:, 3] .- M ./ .(S55_𝐱.^2 .+ S55_𝐲.^2) .+ (S55_𝐯𝐱.^2 .+ S55_𝐯𝐲.^2) ./ 2.0) .* c_to_kms
return (S2_RA, S2_Dec, S2_RV, S38_RA, S38_Dec, S38_RV, S55_RA, S55_Dec, S55_RV)
end
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment