Commit df511007 authored by Rustam Gainutdinov's avatar Rustam Gainutdinov
Browse files

first omega-ms plot obtained

parent 1a4b1652
......@@ -20,11 +20,11 @@ module Cornerplot
end
function cornerplot(fchain, lnames, lnames_ed, hbinsize, cmap, errsize)
function cornerplot(fchain, lnames, lnames_ed, hbinsize, cmap, errsize; fsize=25)
println("Printing...")
n = size(fchain)[1]
figure(figsize = (25,25))
figure(figsize = (fsize,fsize))
fchmeds, f_lerr, f_rerr = chain_estimates(fchain)
fchmedsstr = string.(round.(fchmeds; sigdigits=3))
f_lerrstr = string.(round.(f_lerr; sigdigits=3))
......
......@@ -27,8 +27,11 @@ module DataSets
S55_ΔRV = sstars_df[!, 19]
sstars_art = sstars_df[!, 20]
sstars_data = [sstars_t, [S2_RA, S2_Dec, S2_RV, S38_RA, S38_Dec, S38_RV, S55_RA, S55_Dec, S55_RV], [S2_ΔRA, S2_ΔDec, S2_ΔRV, S38_ΔRA, S38_ΔDec, S38_ΔRV, S55_ΔRA, S55_ΔDec, S55_ΔRV]]
sstars_vals = [[S2_RA[i], S2_Dec[i], S2_RV[i], S38_RA[i], S38_Dec[i], S38_RV[i], S55_RA[i], S55_Dec[i], S55_RV[i]] for i in 1:length(S2_RA)]
sstars_errs = [[S2_ΔRA[i], S2_ΔDec[i], S2_ΔRV[i], S38_ΔRA[i], S38_ΔDec[i], S38_ΔRV[i], S55_ΔRA[i], S55_ΔDec[i], S55_ΔRV[i]] for i in 1:length(S2_RA)]
θ₀ = [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]
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]
end
......@@ -8,7 +8,7 @@ module Models
using Roots
using .EoMs
export model_massless
export model_massless, c_to_kms
rad_to_arcsec = 3600.0 * 1000.0 * 180.0 / π
c_to_kms = 299792.458
......@@ -62,27 +62,45 @@ module Models
end
function Rømer_delay(sol, Ω, i, t_obs)
function Rømer_delay_strict(sol, Ω, i, t_obs)
z = t -> turnorbit(sol(t)[1], sol(t)[2], 0.0, Ω, i)[3]
f = t_em_arg -> (t_em_arg - z(t_em_arg) - t_obs)
f = t_em_arg -> (t_em_arg + z(t_em_arg) - t_obs)
t_em = find_zero(f, t_obs)
return t_em
convkey = try
find_zero(f, t_obs)
catch e
e
end
if typeof(convkey) == Roots.ConvergenceFailed
return t_obs
else
return find_zero(f, t_obs)
end
end
function Rømer_delay(sol, Ω, i, 𝐭_obs)
z = t -> turnorbit(sol(t)[1], sol(t)[2], 0.0, Ω, i)[3]
return z.(𝐭_obs) .+ 𝐭_obs
end
function model_massless(θ, 𝐭_input)
𝐭_input = 𝐭_input ./ t̂_to_yr
S2_x0, S2_y0, S2_vx0, S2_vy0, S2_i, S2_Ω, S38_x0, S38_y0, S38_vx0, S38_vy0, S38_i, S38_Ω, S55_x0, S55_y0, S55_vx0, S55_vy0, S55_i, S55_Ω, M, R0, x_SgrA, y_SgrA, vx_SgrA, vy_SgrA, β, γ = θ
S2_x0, S2_y0, S2_vx0, S2_vy0, S2_i, S2_Ω, S38_x0, S38_y0, S38_vx0, S38_vy0, S38_i, S38_Ω, S55_x0, S55_y0, S55_vx0, S55_vy0, S55_i, S55_Ω, M, R0, vx_SgrA, vy_SgrA, β, γ = θ
S2_sol = solve_EoM_massless(S2_x0, S2_y0, S2_vx0 / c_to_kms, S2_vy0 / c_to_kms, M, β, γ)
S38_sol = solve_EoM_massless(S38_x0, S38_y0, S38_vx0 / c_to_kms, S38_vy0 / c_to_kms, M, β, γ)
S55_sol = solve_EoM_massless(S55_x0, S55_y0, S55_vx0 / c_to_kms, S55_vy0 / c_to_kms, 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)
𝐭_em_S2 = Rømer_delay(S2_sol, S2_Ω, S2_i, 𝐭_input)
𝐭_em_S38 = Rømer_delay(S38_sol, S38_Ω, S38_i, 𝐭_input)
𝐭_em_S55 = Rømer_delay(S55_sol, S55_Ω, S55_i, 𝐭_input)
S2_𝐱, S2_𝐲, S2_𝐯𝐱, S2_𝐯𝐲 = [hcat(S2_sol.(𝐭_em_S2)...)[i,:] for i in 1:4]
S38_𝐱, S38_𝐲, S38_𝐯𝐱, S38_𝐯𝐲 = [hcat(S38_sol.(𝐭_em_S38)...)[i,:] for i in 1:4]
......@@ -97,19 +115,19 @@ module Models
Dist = R0 / r̂_to_kpc
S2_RA = -((hcat(S2_𝐫...)[1,:] .+ vx_SgrA .* (𝐭_em_S2 .- t_beg)) ./ Dist) .* rad_to_arcsec
S2_Dec = ((hcat(S2_𝐫...)[2,:] .+ vy_SgrA .* (𝐭_em_S2 .- t_beg)) ./ Dist) .* rad_to_arcsec
S2_RA = ((hcat(S2_𝐫...)[1,:] .+ vx_SgrA .* (𝐭_em_S2 .- t_beg)) ./ Dist) .* rad_to_arcsec
S2_Dec = -((hcat(S2_𝐫...)[2,:] .+ vy_SgrA .* (𝐭_em_S2 .- t_beg)) ./ Dist) .* rad_to_arcsec
S2_RV = (-hcat(S2_𝐯...)[3,:] .- M ./ .(S2_𝐱.^2 .+ S2_𝐲.^2) .+ (S2_𝐯𝐱.^2 .+ S2_𝐯𝐲.^2) ./ 2.0) .* c_to_kms
S38_RA = -((hcat(S38_𝐫...)[1,:] .+ vx_SgrA .* (𝐭_em_S38 .- t_beg)) ./ Dist) .* rad_to_arcsec
S38_Dec = ((hcat(S38_𝐫...)[2,:] .+ vy_SgrA .* (𝐭_em_S38 .- t_beg)) ./ Dist) .* rad_to_arcsec
S38_RA = ((hcat(S38_𝐫...)[1,:] .+ vx_SgrA .* (𝐭_em_S38 .- t_beg)) ./ Dist) .* rad_to_arcsec
S38_Dec = -((hcat(S38_𝐫...)[2,:] .+ vy_SgrA .* (𝐭_em_S38 .- t_beg)) ./ Dist) .* rad_to_arcsec
S38_RV = (-hcat(S38_𝐯...)[3,:] .- M ./ .(S38_𝐱.^2 .+ S38_𝐲.^2) .+ (S38_𝐯𝐱.^2 .+ S38_𝐯𝐲.^2) ./ 2.0) .* c_to_kms
S55_RA = -((hcat(S55_𝐫...)[1,:] .+ vx_SgrA .* (𝐭_em_S55 .- t_beg)) ./ Dist) .* rad_to_arcsec
S55_Dec = ((hcat(S55_𝐫...)[2,:] .+ vy_SgrA .* (𝐭_em_S55 .- t_beg)) ./ Dist) .* rad_to_arcsec
S55_RA = ((hcat(S55_𝐫...)[1,:] .+ vx_SgrA .* (𝐭_em_S55 .- t_beg)) ./ Dist) .* rad_to_arcsec
S55_Dec = -((hcat(S55_𝐫...)[2,:] .+ vy_SgrA .* (𝐭_em_S55 .- t_beg)) ./ Dist) .* rad_to_arcsec
S55_RV = (-hcat(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]
return [[S2_RA[i], S2_Dec[i], S2_RV[i], S38_RA[i], S38_Dec[i], S38_RV[i], S55_RA[i], S55_Dec[i], S55_RV[i]] for i in 1:length(S2_RA)]
end
......
This diff is collapsed.
......@@ -15,16 +15,18 @@ module ParEst
function lnlike(θ, input_data)
t, x, Δx = input_data
χ² = ((x .- model_massless(θ, t)) ./ Δx) .^ 2
t, x_raw, Δx_raw = input_data
x = hcat(x_raw...)
Δx = hcat(Δx_raw...)
χ² = ((x .- hcat(model_massless(θ, t)...)) ./ Δx) .^ 2
return -0.5 * sum(χ²)
end
function lnprior(θ)
S2_x0, S2_y0, S2_vx0, S2_vy0, S2_i, S2_Ω, S38_x0, S38_y0, S38_vx0, S38_vy0, S38_i, S38_Ω, S55_x0, S55_y0, S55_vx0, S55_vy0, S55_i, S55_Ω, M, R0, x_SgrA, y_SgrA, vx_SgrA, vy_SgrA, β, γ = θ
bcond = 0.0 S2_i 180.0 && 0.0 S38_i 180.0 && 0.0 S55_i 180.0 && 0.0 S2_Ω 360.0 && 0.0 S38_Ω 360.0 && 0.0 S55_Ω 360.0 && √sum(S2_vx0^2 + S2_vy0^2) < 1.0 && √sum(S38_vx0^2 + S38_vy0^2) < 1.0 && √sum(S55_vx0^2 + S55_vy0^2) < 1.0 && M > 0 && 0.0 < β < 2.0 && 0.0 < γ < 2.0
S2_x0, S2_y0, S2_vx0, S2_vy0, S2_i, S2_Ω, S38_x0, S38_y0, S38_vx0, S38_vy0, S38_i, S38_Ω, S55_x0, S55_y0, S55_vx0, S55_vy0, S55_i, S55_Ω, M, R0, vx_SgrA, vy_SgrA, β, γ = θ
bcond = 0.0 S2_i 180.0 && 0.0 S38_i 180.0 && 0.0 S55_i 180.0 && 0.0 S2_Ω 360.0 && 0.0 S38_Ω 360.0 && 0.0 S55_Ω 360.0 && √sum(S2_vx0^2 + S2_vy0^2) < c_to_kms / 10.0 && √sum(S38_vx0^2 + S38_vy0^2) < c_to_kms / 10.0 && √sum(S55_vx0^2 + S55_vy0^2) < c_to_kms / 10.0 && M > 0 && 0.0 < β < 2.0 && 0.0 < γ < 2.0
if bcond
return 0.0
......
__precompile__()
module ScalarTensor
export ω_BD
AU_to_r_planck = 9.25701e+45
eV_to_m_planck = 1.78266e-36 / 2.17645e-8
function ω_BD_pl(mₛ_in, r_in, γ) #mₛ in eV, r in AU
mₛr = mₛ_in * eV_to_m_planck * r_in * AU_to_r_planck
return exp(-mₛr) / 2.0 * (1 + γ) / (1 - γ)
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