Inclination fix (#84)

* fix to nu+omega range * added more inclinations to test range, tests now fail * version * removed abs() from Lvec in calculation of i_f, in both Turing.jl and Symbolics. Now the final inclination matches the initial one when > 90deg * added safety catch in i_f calculation, to make sure Lvec[3] is between -1 and 1 before calculating acos()
orlox · Oct 17, 2024 · 3da22aa · 3da22aa
1 parent cecb538
commit 3da22aa
Show file tree

Hide file tree

Showing 5 changed files with 5 additions and 5 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "SideKicks"
 uuid = "ae7baef6-f8d7-4188-9988-8bf5c4952615"
 authors = ["Pablo Marchant ", "Reinhold Willcox "]
-version = "0.3.7"
+version = "0.3.8"
 
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"

diff --git a/src/Orbits.jl b/src/Orbits.jl
@@ -337,7 +337,7 @@ function post_supernova_general_orbit_parameters(;m1_i, m2_i, a_i, e_i=0, m1_f=-
     L_n = R_n_par*L_par + R_n_per*L_per + R_n_z*L_z
     L_o = R_o_par*L_par + R_o_per*L_per + R_o_z*L_z
 
-    i_f = acos(min(1, abs(L_o)))
+    i_f = acos(max(-1,min(1, L_o)))
 
     # compute longitude of the ascending node
     # this is obtained by taking the vector Ωvec_f = Lvec_f × \hat{O} 

diff --git a/src/TuringModels.jl b/src/TuringModels.jl
@@ -363,7 +363,7 @@ function create_general_mcmc_model(;
         end
 
         dm2 = m2_i - m2_f
-        sum_ωi_νi = (ω_i + ν_i) % 360 # should be between 0 and 360 deg
+        sum_ωi_νi = (ω_i + ν_i) % 2π # should be between 0 and 2π
 
         return     ( m1_i,  m2_i,  P_i,  e_i,  a_i,  Ω_i,  ω_i,  i_i,  sum_ωi_νi,
                      vkick,  θ,  ϕ,  dm2,  frac,  ν_i, 

diff --git a/test/SymbolicFuncs.jl b/test/SymbolicFuncs.jl
@@ -112,7 +112,7 @@ function symbolic_kick_functions_vcm_and_orbital_elements()
         uncorrected_ω_final, 2*π+uncorrected_ω_final)
 
     #orbital inclination is an easy one
-    ι_final = acos(min(1,abs(Lvec_rot[3]/sqrt(L_x^2+L_y^2+L_z^2))))
+    ι_final = acos(max(-1,min(1,Lvec_rot[3]/sqrt(L_x^2+L_y^2+L_z^2))))
 
     #motion of the center of mass
     v_cm_rot = Rtotal*[v_xcm,v_ycm,v_zcm]

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -115,7 +115,7 @@ for vkick in [0, 30*km_per_s]
     for vimp in [0, 20*km_per_s]
         for Ω_i in [0.1π, 0.9π, 1.5π, 1.9π]
             for ω_i in [0.1π, 0.9π, 1.5π, 1.9π]
-                for i_i in [0.1π/2, 0.5π/2, 0.9π/2]
+                for i_i in [0.1π/2, 0.5π/2, 0.9π/2, 1.1π/2, 1.5π/2, 1.9π/2]
                     for ν_i in [0.1π, 0.9π, 1.5π, 1.9π]
                         for m1_f in [0.9*m1_i, 0.95*m1_i, m1_i]
                             for m2_f in [0.9*m2_i, 0.95*m2_i, m2_i]