Skip to content

Commit

Permalink
Merge branch 'main' of https://github.com/kima-org/kima
Browse files Browse the repository at this point in the history
  • Loading branch information
@j-faria
j-faria committed Dec 10, 2024
2 parents 8663a61 + 66b51f1 commit 2f36f89
Show file tree
Hide file tree
Showing 6 changed files with 105 additions and 70 deletions.
6 changes: 3 additions & 3 deletions src/kima/ETVmodel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ void ETVmodel::calculate_mu()
ecc = components[j][3];
omega = components[j][4];

auto tau = brandt::keplerian_et(epochs, P, K, ecc, omega, phi, ephem1);
auto tau = brandt::keplerian_etv(epochs, P, K, ecc, omega, phi, ephem1);
for(size_t i=0; i<N; i++)
mu[i] += tau[i]/(24*3600);
}
Expand All @@ -173,7 +173,7 @@ void ETVmodel::remove_known_object()
auto epochs = data.get_epochs();
for (int j = 0; j < n_known_object; j++)
{
auto tau = brandt::keplerian_et(epochs, KO_P[j], KO_K[j], KO_e[j], KO_w[j], KO_phi[j], ephem1);
auto tau = brandt::keplerian_etv(epochs, KO_P[j], KO_K[j], KO_e[j], KO_w[j], KO_phi[j], ephem1);
for (size_t i = 0; i < data.N(); i++) {
mu[i] -= tau[i]/(24*3600);
}
Expand All @@ -187,7 +187,7 @@ void ETVmodel::add_known_object()
std::vector<double> ts;
for (int j = 0; j < n_known_object; j++)
{
auto tau = brandt::keplerian(epochs, KO_P[j], KO_K[j], KO_e[j], KO_w[j], KO_phi[j], ephem1);
auto tau = brandt::keplerian_etv(epochs, KO_P[j], KO_K[j], KO_e[j], KO_w[j], KO_phi[j], ephem1);
for (size_t i = 0; i < data.N(); i++) {
mu[i] += tau[i]/(24*3600);
}
Expand Down
21 changes: 9 additions & 12 deletions src/kima/kepler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -740,7 +740,7 @@ namespace brandt
return rv;
}

std::vector<double> keplerian_et(const std::vector<double> &epochs, const double &P,
std::vector<double> keplerian_etv(const std::vector<double> &epochs, const double &P,
const double &K, const double &ecc,
const double &w, const double &M0,
const double &ephem1)
Expand All @@ -754,9 +754,6 @@ namespace brandt
double sinw, cosw;
sincos(w, &sinw, &cosw);

// ecentricity factor for g, once per orbit
double g_e = sqrt((1 + ecc) / (1 - ecc));

// brandt solver calculations, once per orbit
double bounds[13];
double EA_tab[6 * 13];
Expand All @@ -765,14 +762,11 @@ namespace brandt
// std::cout << std::endl;
for (size_t i = 0; i < epochs.size(); i++)
{
double sinE, cosE;
double sinEf, cosEf;
double M = n * (epochs[i]*ephem1) + M0;
solver_fixed_ecc(bounds, EA_tab, M, ecc, &sinE, &cosE);
double g = g_e * ((1 - cosE) / sinE);
double g2 = g * g;
// std::cout << M << '\t' << ecc << '\t' << sinE << '\t' << cosE << std::endl;
// std::cout << '\t' << g << '\t' << g2 << std::endl;
ets[i] = K / (pow(1- pow(ecc * cosw,2.0),0.5)) * ((1-ecc*ecc)/(1+g2+ecc-ecc*g2) * (2*g*cosw + (1-g2)*sinw) + ecc*sinw);
solver_fixed_ecc(bounds, EA_tab, M, ecc, &sinEf, &cosEf);
to_f(ecc, 1 - ecc, &sinEf, &cosEf);
ets[i] = K / (pow(1- pow(ecc * cosw,2.0),0.5)) * ((1-ecc*ecc)/(1+ecc*cosEf) * (sinEf*cosw + cosEf*sinw) + ecc*sinw);
}

return ets;
Expand Down Expand Up @@ -1188,9 +1182,12 @@ NB_MODULE(kepler, m) {
m.def("keplerian2", &brandt::keplerian2,
"t"_a, "P"_a, "K"_a, "ecc"_a, "w"_a, "M0"_a, "M0_epoch"_a);

m.def("keplerian_etv", &brandt::keplerian_et,
m.def("keplerian_etv", &brandt::keplerian_etv,
"epochs"_a, "P"_a, "K"_a, "ecc"_a, "w"_a, "M0"_a, "ephem1"_a);

m.def("keplerian_gaia", &brandt::keplerian_gaia,
"t"_a, "psi"_a, "A"_a, "B"_a,"F"_a, "G"_a, "ecc"_a, "P"_a, "M0"_a, "M0_epoch"_a);

// [](nb::ndarray<double, nb::shape<nb::any>, nb::device::cpu> t,
// const double &P, const double &K, const double &ecc,
// const double &w, const double &M0, const double &M0_epoch)
Expand Down
2 changes: 1 addition & 1 deletion src/kima/kepler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ namespace brandt
const double &B, const double &F, const double &G,
const double &ecc, const double P, const double &M0,
const double &M0_epoch);
std::vector<double> keplerian_et(const std::vector<double> &t, const double &P,
std::vector<double> keplerian_etv(const std::vector<double> &t, const double &P,
const double &K, const double &ecc,
const double &w, const double &M0,
const double &M0_epoch);
Expand Down
89 changes: 41 additions & 48 deletions src/kima/postkepler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,39 +93,39 @@ namespace postKep
return K2;
}

inline double light_travel_time(double K1, double f, double w, double ecc)
inline double light_travel_time(double K1, double sinf, double cosf, double w, double ecc)
{

double delta_LT = pow(K1,2.0)*pow(sin(f + w), 2.0)*(1+ecc*cos(f))/c_light;
double delta_LT = pow(K1,2.0)*pow(sinf*cos(w) + cosf*sin(w), 2.0)*(1+ecc*cosf)/c_light;

return delta_LT;
}

inline double transverse_doppler(double K1, double f, double ecc, double cosi)
inline double transverse_doppler(double K1, double sinf, double cosf, double ecc, double cosi)
{
double sin2i = 1.0 - cosi*cosi;
double delta_TD = pow(K1,2.0)*(1 + ecc*cos(f) - (1-pow(ecc,2.0))/2)/(c_light*sin2i);
double delta_TD = pow(K1,2.0)*(1 + ecc*cosf - (1-pow(ecc,2.0))/2)/(c_light*sin2i);

return delta_TD;
}

inline double gravitational_redshift(double K1, double K2, double f, double ecc, double cosi)
inline double gravitational_redshift(double K1, double K2, double sinf, double cosf, double ecc, double cosi)
{
double sin2i = 1.0 - cosi*cosi;
double delta_GR = K1*(K1+K2)*(1+ecc*cos(f))/(c_light*sin2i);
double delta_GR = K1*(K1+K2)*(1+ecc*cosf)/(c_light*sin2i);

return delta_GR;
}

inline double v_tide(double R1, double M1, double M2, double P, double f, double w, double cosi)
inline double v_tide(double R1, double M1, double M2, double P, double sinf, double cosf, double w, double cosi)
{
double phi_0 = M_PI/2 - w;
double sin2i = 1.0 - cosi*cosi;

return 1184*M2/(M1*(M1+M2))*pow(R1,4.0)*pow(P,-3.0)*sin(2*(f-phi_0))*sin2i;
return 1184*M2/(M1*(M1+M2))*pow(R1,4.0)*pow(P,-3.0)*sin2i * 2*(sinf*cos(phi_0) - sin(phi_0)*cosf)*(cosf*cos(phi_0) + sinf*sin(phi_0));
}

double post_Newtonian(double K1, double f, double ecc, double w, double P, double cosi, double M1, double M2, double R1, bool GR, bool Tid)
double post_Newtonian(double K1, double sinf, double cosf, double ecc, double w, double P, double cosi, double M1, double M2, double R1, bool GR, bool Tid)
{
double K2;
double v = 0.0;
Expand All @@ -148,29 +148,29 @@ namespace postKep
{
R1 = pow(M1,0.8);
}
double delta_v_tide = v_tide(R1,M1,M2,P,f,w,cosi);
double delta_v_tide = v_tide(R1,M1,M2,P,sinf,cosf,w,cosi);
v = v + delta_v_tide;
}
if (GR)
{
double delta_LT = light_travel_time(K1, f, w, ecc);
double delta_TD = transverse_doppler(K1, f, ecc,cosi);
double delta_GR = gravitational_redshift(K1, K2, f, ecc,cosi);
double delta_LT = light_travel_time(K1, sinf,cosf, w, ecc);
double delta_TD = transverse_doppler(K1, sinf,cosf, ecc,cosi);
double delta_GR = gravitational_redshift(K1, K2, sinf,cosf, ecc,cosi);
v = v + delta_LT + delta_TD + delta_GR;
}


return v;
}

std::tuple <double,double> post_Newtonian_sb2(double K1,double K2, double f, double ecc, double w, double P, double cosi, double q, double R1, double R2, bool GR, bool Tid)
std::tuple <double,double> post_Newtonian_sb2(double K1,double K2, double sinf, double cosf, double ecc, double w, double P, double cosi, double q, double R1, double R2, bool GR, bool Tid)
{

double v1 = 0.0;
double v2 = 0.0;
double sin2i = 1.0 - cosi*cosi;
//calculate masses
double M1 = P * days * pow(K2,3.0) * pow((1 - pow(ecc,2.0)),1.5) * pow((1 + q),2.0) * pow(sin2i,1.5) / (TWO_PI * G);
double M1 = P * days * pow(K2,3.0) * pow((1 - pow(ecc,2.0)),1.5) * pow((1 + q),2.0) / pow(sin2i,1.5) / (TWO_PI * G)/Msun;
double M2 = M1 * q;
if (Tid)
{
Expand All @@ -182,21 +182,21 @@ namespace postKep
{
R2 = pow(M2,0.8);
}
double delta_v_tide1 = v_tide(R1,M1,M2,P,f,w,cosi);
double delta_v_tide1 = v_tide(R1,M1,M2,P,sinf,cosf,w,cosi);
v1 = v1 + delta_v_tide1;
double delta_v_tide2 = v_tide(R2,M2,M1,P,f,w-M_PI,cosi);
double delta_v_tide2 = v_tide(R2,M2,M1,P,sinf,cosf,w-M_PI,cosi);
v2 = v2 + delta_v_tide2;
}
if (GR)
{
double delta_LT = light_travel_time(K1, f, w, ecc);
double delta_TD = transverse_doppler(K1, f, ecc,cosi);
double delta_GR = gravitational_redshift(K1, K2, f, ecc,cosi);
double delta_LT = light_travel_time(K1, sinf,cosf, w, ecc);
double delta_TD = transverse_doppler(K1, sinf,cosf, ecc,cosi);
double delta_GR = gravitational_redshift(K1, K2, sinf,cosf, ecc,cosi);
v1 = v1 + delta_LT + delta_TD + delta_GR;

double delta_LT2 = light_travel_time(K2, f, w - M_PI, ecc);
double delta_TD2 = transverse_doppler(K2, f, ecc,cosi);
double delta_GR2 = gravitational_redshift(K2, K1, f, ecc,cosi);
double delta_LT2 = light_travel_time(K2, sinf,cosf, w - M_PI, ecc);
double delta_TD2 = transverse_doppler(K2, sinf,cosf, ecc,cosi);
double delta_GR2 = gravitational_redshift(K2, K1, sinf,cosf, ecc,cosi);
v2 = v2 + delta_LT2 + delta_TD2 + delta_GR2;
}

Expand All @@ -217,10 +217,6 @@ namespace postKep

// mean motion, once per orbit
double n = 2. * M_PI / P;


// ecentricity factor for g, once per orbit
double g_e = sqrt((1 + ecc) / (1 - ecc));

// brandt solver calculations, once per orbit
double bounds[13];
Expand All @@ -236,18 +232,14 @@ namespace postKep
double sinw, cosw;
sincos(w_t, &sinw, &cosw);

double sinE, cosE;
double sinEf, cosEf;
double M = n * (t[i] - M0_epoch) + M0;
brandt::solver_fixed_ecc(bounds, EA_tab, M, ecc, &sinE, &cosE);
double g = g_e * ((1 - cosE) / sinE);
double g2 = g * g;
brandt::solver_fixed_ecc(bounds, EA_tab, M, ecc, &sinEf, &cosEf);
brandt::to_f(ecc, 1-ecc, &sinEf, &cosEf);

double vrad = K * (cosw * ((1 - g2) / (1 + g2) + ecc) - sinw * ((2 * g) / (1 + g2)));


double f = acos((cosE - ecc) / (1 - ecc * cosE));
double vrad = K * (cosw * (cosEf + ecc) - sinw * sinEf);

double v_correction = postKep::post_Newtonian(K, f, ecc, w_t, P, cosi, M1, M2, R1, GR, Tid);
double v_correction = postKep::post_Newtonian(K, sinEf,cosEf, ecc, w_t, P, cosi, M1, M2, R1, GR, Tid);

rv[i] = vrad + v_correction;
}
Expand All @@ -268,10 +260,6 @@ namespace postKep

// mean motion, once per orbit
double n = 2. * M_PI / P;


// ecentricity factor for g, once per orbit
double g_e = sqrt((1 + ecc) / (1 - ecc));

// brandt solver calculations, once per orbit
double bounds[13];
Expand All @@ -287,18 +275,16 @@ namespace postKep
double sinw, cosw;
sincos(w_t, &sinw, &cosw);

double sinE, cosE;
double sinEf, cosEf;
double M = n * (t[i] - M0_epoch) + M0;
brandt::solver_fixed_ecc(bounds, EA_tab, M, ecc, &sinE, &cosE);
double g = g_e * ((1 - cosE) / sinE);
double g2 = g * g;
brandt::solver_fixed_ecc(bounds, EA_tab, M, ecc, &sinEf, &cosEf);
brandt::to_f(ecc, 1-ecc, &sinEf, &cosEf);

double vrad1 = K * (cosw * ((1 - g2) / (1 + g2) + ecc) - sinw * ((2 * g) / (1 + g2)));
double vrad2 = K/q * (sinw * ((2 * g) / (1 + g2)) - cosw * ((1 - g2) / (1 + g2) + ecc));
double vrad1 = K * (cosw * (cosEf + ecc) - sinw * sinEf);
double vrad2 = K/q * (-cosw * (cosEf + ecc) + sinw * sinEf);

double f = acos((cosE - ecc) / (1 - ecc * cosE));

auto [v_correction1,v_correction2] = postKep::post_Newtonian_sb2(K, K/q, f, ecc, w_t, P, cosi, q, R1, R2, GR, Tid);
auto [v_correction1,v_correction2] = postKep::post_Newtonian_sb2(K, K/q, sinEf,cosEf, ecc, w_t, P, cosi, q, R1, R2, GR, Tid);

rv1[i] = vrad1 + v_correction1;
rv2[i] = vrad2 + v_correction2;
Expand Down Expand Up @@ -328,3 +314,10 @@ namespace MassConv
}
}

NB_MODULE(postkepler, m) {
m.def("keplerian_prec", &postKep::keplerian_prec,
"t"_a, "P"_a, "K"_a, "ecc"_a, "w"_a, "wdot"_a, "M0"_a, "M0_epoch"_a,"cosi"_a, "M1"_a, "M2"_a, "R1"_a, "GR"_a, "Tid"_a);
m.def("keplerian_prec_sb2", &postKep::keplerian_prec_sb2,
"t"_a, "P"_a, "K"_a, "q"_a, "ecc"_a, "w"_a, "wdot"_a, "M0"_a, "M0_epoch"_a,"cosi"_a, "R1"_a, "R2"_a, "GR"_a, "Tid"_a);

}
12 changes: 6 additions & 6 deletions src/kima/postkepler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,12 @@ namespace postKep
inline double f_dash_M(double K,double M0, double M1, double P, double ecc);
inline double get_K2_v1(double K1, double M, double P, double ecc);
inline double get_K2_v2(double K1, double M, double P, double ecc);
inline double light_travel_time(double K1, double f, double w, double ecc);
inline double transverse_doppler(double K1, double f, double ecc, double cosi);
inline double gravitational_redshift(double K1, double K2, double f, double ecc, double cosi);
inline double v_tide(double R1, double M1, double M2, double P, double f, double w);
double post_Newtonian(double K1, double f, double ecc, double w, double P, double M1, double M2, double R1, bool GR, bool Tid);
double post_Newtonian_sb2(double K1, double K2, double f, double ecc, double w, double P, double q, double R1, double R2, bool GR, bool Tid);
inline double light_travel_time(double K1, double sinf, double cosf, double w, double ecc);
inline double transverse_doppler(double K1, double sinf, double cosf, double ecc, double cosi);
inline double gravitational_redshift(double K1, double K2, double sinf, double cosf, double ecc, double cosi);
inline double v_tide(double R1, double M1, double M2, double P, double sinf, double cosf, double w);
double post_Newtonian(double K1, double sinf, double cosf, double ecc, double w, double P, double M1, double M2, double R1, bool GR, bool Tid);
double post_Newtonian_sb2(double K1, double K2, double sinf, double cosf, double ecc, double w, double P, double q, double R1, double R2, bool GR, bool Tid);

std::vector<double> keplerian_prec(const std::vector<double> &t, const double &P,
const double &K, const double &ecc,
Expand Down
45 changes: 45 additions & 0 deletions tests/test_keplerian.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,51 @@ def test_keplerian():
kepler.keplerian2(t, P, K, ecc, w, 0.0, 0.0)
)

def test_keplerian_etv():
from kima.kepler import keplerian_etv

epochs = np.linspace(0, 1, 10)
P, K, ecc, w, eph1 = 1.0, 1.0, 0.0, np.pi/2, 0.01

assert np.isfinite(keplerian_etv([0.0], P, K, ecc, w, 0.0, eph1))
assert np.isfinite(keplerian_etv(epochs, P, K, ecc, w, 0.0, eph1)).all()
np.testing.assert_allclose(
keplerian_etv([0.0, P/2/eph1, P/eph1], P, K, ecc, w, 0.0, eph1),
[K, -K, K]
)

ecc = 0.1
cnu1 = -0.19738041725338823
snu1 = 0.9803269714155979
np.testing.assert_allclose(
keplerian_etv([0.0, P/4/eph1, P/2/eph1], P, K, ecc, w, 0.0, eph1),
[K, K*((1-ecc**2)/(1+ecc*cnu1)*cnu1+ecc), -K]
)

w=0
np.testing.assert_allclose(
keplerian_etv([0.0, P/4/eph1, P/2/eph1], P, K, ecc, w, 0.0, eph1),
[0, K*((1-ecc**2)**0.5)/(1+ecc*cnu1)*snu1, 0]
)

def test_keplerian_gaia():
from kima.kepler import keplerian_gaia

t = np.linspace(0, 1, 10)
psi = np.linspace(0,2*np.pi,10)
P, ecc, A, B, F, G = 1.0, 0.0, 1.0, -1.0, 1.0, -1.0

assert np.isfinite(keplerian_gaia([0.0],[0.0],A,B,F,G,ecc,P,0.0,0.0))
assert np.isfinite(keplerian_gaia(t,psi,A,B,F,G,ecc,P,0.0,0.0)).all()


ecc=0.1

np.testing.assert_allclose(
keplerian_gaia([0.0,P/4 - ecc*P/2/np.pi,P/2],[0.0,0.0,np.pi/2],A,B,F,G,ecc,P,0.0,0.0),
[1-ecc, (1-ecc**2)**(0.5) - ecc, 1+ecc]
)


@pytest.mark.skip(reason="only for testing")
def test_speed():
Expand Down

0 comments on commit 2f36f89

Please sign in to comment.