Skip to content

Commit

Permalink
redo postkepler dependency
Browse files Browse the repository at this point in the history
  • Loading branch information
@ThomasBaycroft
ThomasBaycroft committed Jan 15, 2024
1 parent c226446 commit 8ca9e5f
Show file tree
Hide file tree
Showing 6 changed files with 61 additions and 59 deletions.
9 changes: 5 additions & 4 deletions src/kima/BINARIESmodel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ using namespace std;
// using namespace Eigen;
using namespace DNest4;
using namespace nijenhuis;
using namespace brandt;
using namespace postKep;

#define TIMING false
Expand Down Expand Up @@ -383,7 +384,7 @@ void BINARIESmodel::remove_known_object()
for(int j=0; j<n_known_object; j++)
{
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
auto v = brandt::keplerian_prec(data.t, P_anom, KO_K[j], KO_e[j], KO_w[j], KO_wdot[j], KO_phi[j], data.M0_epoch);
auto v = postKep::keplerian_prec(data.t, P_anom, KO_K[j], KO_e[j], KO_w[j], KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
mu[i] -= v[i];
Expand Down Expand Up @@ -414,7 +415,7 @@ void BINARIESmodel::remove_known_object_secondary()
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
K2 = KO_K[j]/KO_q[j];
w = KO_w[j] - M_PI;
auto v = brandt::keplerian_prec(data.t, P_anom, K2, KO_e[j], w, KO_wdot[j], KO_phi[j], data.M0_epoch);
auto v = postKep::keplerian_prec(data.t, P_anom, K2, KO_e[j], w, KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
mu_2[i] -= v[i];
Expand Down Expand Up @@ -442,7 +443,7 @@ void BINARIESmodel::add_known_object()
for(int j=0; j<n_known_object; j++)
{
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
auto v = brandt::keplerian_prec(data.t, P_anom, KO_K[j], KO_e[j], KO_w[j], KO_wdot[j], KO_phi[j], data.M0_epoch);
auto v = postKep::keplerian_prec(data.t, P_anom, KO_K[j], KO_e[j], KO_w[j], KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
mu[i] += v[i];
Expand Down Expand Up @@ -471,7 +472,7 @@ void BINARIESmodel::add_known_object_secondary()
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
K2 = KO_K[j]/KO_q[j];
w = KO_w[j] - M_PI;
auto v = brandt::keplerian_prec(data.t, P_anom, K2, KO_e[j], w, KO_wdot[j], KO_phi[j], data.M0_epoch);
auto v = postKep::keplerian_prec(data.t, P_anom, K2, KO_e[j], w, KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
mu_2[i] += v[i];
Expand Down
2 changes: 1 addition & 1 deletion src/kima/examples/Kepler16/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ def Kepler16(run=False, **kwargs):
return model

if __name__ == '__main__':
model = Kepler16(run=True, steps=100000)
model = Kepler16(run=True, steps=1000)
res = kima.load_results()


49 changes: 0 additions & 49 deletions src/kima/kepler.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
#include "kepler.h"
using namespace postKep;

const double TWO_PI = M_PI * 2;
const double PI_D_2 = M_PI / 2;
Expand Down Expand Up @@ -753,54 +752,6 @@ namespace brandt
}


//
std::vector<double> keplerian_prec(const std::vector<double> &t, const double &P,
const double &K, const double &ecc,
const double &w, const double &wdot, const double &M0,
const double &M0_epoch)
{
// allocate RVs
std::vector<double> rv(t.size());


// 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];
double EA_tab[6 * 13];
get_bounds(bounds, EA_tab, ecc);

// std::cout << std::endl;
for (size_t i = 0; i < t.size(); i++)
{
double Tp = M0_epoch-(P*M0)/(TWO_PI);
double w_t = postKep::change_omega(w, wdot, t[i], Tp);
// sin and cos of argument of periastron
double sinw, cosw;
sincos(w_t, &sinw, &cosw);

double sinE, cosE;
double M = n * (t[i] - M0_epoch) + 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;
rv[i] = K * (cosw * ((1 - g2) / (1 + g2) + ecc) - sinw * ((2 * g) / (1 + g2)));

// double f ;

// double v_correction = postKep::post_Newtonian(K, f, ecc, w_t, P, double M1, double M2, double R1, bool GR, bool Tid);
}

return rv;
}

std::vector<double> keplerian_gaia(const std::vector<double> &t, const std::vector<double> &psi, const double &A,
const double &B, const double &F, const double &G,
const double &ecc, const double P, const double &M0,
Expand Down
5 changes: 0 additions & 5 deletions src/kima/kepler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
#include <algorithm>
#include <iostream>
#include <tuple>
#include "postkepler.h"

// sincos on apple
#ifdef __APPLE__
Expand Down Expand Up @@ -71,10 +70,6 @@ namespace brandt
const double &K, const double &ecc,
const double &w, const double &M0,
const double &M0_epoch);
std::vector<double> keplerian_prec(const std::vector<double> &t, const double &P,
const double &K, const double &ecc,
const double &w, const double &wdot, const double &M0,
const double &M0_epoch);
std::vector<double> keplerian_gaia(const std::vector<double> &t, const std::vector<double> &psi, const double &A,
const double &B, const double &F, const double &G,
const double &ecc, const double P, const double &M0,
Expand Down
49 changes: 49 additions & 0 deletions src/kima/postkepler.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#include "postkepler.h"

using namespace brandt;

const double TWO_PI = M_PI * 2;

namespace postKep
Expand Down Expand Up @@ -160,4 +162,51 @@ namespace postKep

return v;
}

//
std::vector<double> keplerian_prec(const std::vector<double> &t, const double &P,
const double &K, const double &ecc,
const double &w, const double &wdot, const double &M0,
const double &M0_epoch)
{
// allocate RVs
std::vector<double> rv(t.size());


// 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];
double EA_tab[6 * 13];
brandt::get_bounds(bounds, EA_tab, ecc);

// std::cout << std::endl;
for (size_t i = 0; i < t.size(); i++)
{
double Tp = M0_epoch-(P*M0)/(TWO_PI);
double w_t = change_omega(w, wdot, t[i], Tp);
// sin and cos of argument of periastron
double sinw, cosw;
sincos(w_t, &sinw, &cosw);

double sinE, cosE;
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;

rv[i] = K * (cosw * ((1 - g2) / (1 + g2) + ecc) - sinw * ((2 * g) / (1 + g2)));

// double f ;

// double v_correction = postKep::post_Newtonian(K, f, ecc, w_t, P, double M1, double M2, double R1, bool GR, bool Tid);
}

return rv;
}
}
6 changes: 6 additions & 0 deletions src/kima/postkepler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

#include <cmath>
#include <iostream>
#include "kepler.h"

namespace postKep
{
Expand All @@ -17,4 +18,9 @@ namespace postKep
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);

std::vector<double> keplerian_prec(const std::vector<double> &t, const double &P,
const double &K, const double &ecc,
const double &w, const double &wdot, const double &M0,
const double &M0_epoch);
}

0 comments on commit 8ca9e5f

Please sign in to comment.