diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml index 8abcc53..1096b00 100644 --- a/.github/workflows/python-publish.yml +++ b/.github/workflows/python-publish.yml @@ -5,13 +5,14 @@ on: release: types: [published] -permissions: - contents: read - jobs: deploy: - runs-on: ubuntu-latest + environment: + name: deploy-to-pypi + url: https://pypi.org/p/kima + permissions: + id-token: write steps: - uses: actions/checkout@v3 @@ -31,7 +32,3 @@ jobs: - name: Publish package uses: pypa/gh-action-pypi-publish@release/v1 - with: - repository-url: https://test.pypi.org/legacy/ - user: __token__ - password: ${{ secrets.PYPIT_API_TOKEN }} diff --git a/src/kima/BINARIESmodel.cpp b/src/kima/BINARIESmodel.cpp index 73a1cba..38f23ad 100644 --- a/src/kima/BINARIESmodel.cpp +++ b/src/kima/BINARIESmodel.cpp @@ -244,21 +244,21 @@ void BINARIESmodel::calculate_mu() omega = components[j][4]; -// auto v = brandt::keplerian(data.t, P, K, ecc, omega, phi, data.M0_epoch); -// for(size_t i=0; i(0, 0.5); + a0prior = make_shared(0.01, 10); if (!omegaprior) omegaprior = make_shared(0, 2*M_PI); if (!cosiprior) cosiprior = make_shared(-1, 1); if (!Omegaprior) - Omegaprior = make_shared(0, 2*M_PI); + Omegaprior = make_shared(0, M_PI); } } diff --git a/src/kima/Data.cpp b/src/kima/Data.cpp index bbc7922..ab01ac2 100644 --- a/src/kima/Data.cpp +++ b/src/kima/Data.cpp @@ -876,8 +876,8 @@ GAIAdata::GAIAdata() {}; pf = data[4]; - // epoch for the mean anomaly, by default the time of the first observation - M0_epoch = t[0]; + // epoch for the mean anomaly, by default the gaia reference time + M0_epoch = 57388.5; // How many points did we read? if (VERBOSE) diff --git a/src/kima/GAIAmodel.cpp b/src/kima/GAIAmodel.cpp index 1e21e4f..d501ee9 100644 --- a/src/kima/GAIAmodel.cpp +++ b/src/kima/GAIAmodel.cpp @@ -4,6 +4,7 @@ using namespace std; // using namespace Eigen; using namespace DNest4; using namespace nijenhuis; +using namespace brandt; #define TIMING false @@ -42,13 +43,13 @@ void GAIAmodel::setPriors() // BUG: should be done by only one thread! Jprior = make_prior(0.1,100.); if (!da_prior) - da_prior = make_prior(0.0,pow(10,1)); + da_prior = make_prior(0.0,pow(10,2)); if (!dd_prior) - dd_prior = make_prior(0.0,pow(10,1)); + dd_prior = make_prior(0.0,pow(10,2)); if (!mua_prior) - mua_prior = make_prior(0.0,pow(10,1)); + mua_prior = make_prior(0.0,pow(10,2)); if (!mud_prior) - mud_prior = make_prior(0.0,pow(10,1)); + mud_prior = make_prior(0.0,pow(10,2)); if (!plx_prior) plx_prior = make_prior(0.01,1000.); @@ -181,26 +182,35 @@ void GAIAmodel::calculate_mu() omega = components[j][4]; cosi = components[j][5]; Omega = components[j][6]; - } - - for(size_t i=0; i rv(t.size()); - //Need to change P -> Panom, and have w move, maybe this slows the process down to much to be worth it? // mean motion, once per orbit double n = 2. * M_PI / P; - // sin and cos of argument of periastron, once per orbit - double sinw, cosw; - sincos(w, &sinw, &cosw); + // ecentricity factor for g, once per orbit double g_e = sqrt((1 + ecc) / (1 - ecc)); @@ -780,18 +778,66 @@ namespace brandt // 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; + 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 keplerian_gaia(const std::vector &t, const std::vector &psi, const double &A, + const double &B, const double &F, const double &G, + const double &ecc, const double P, const double &M0, + const double &M0_epoch) + { + // allocate wks + std::vector wk(t.size()); + + // mean motion, once per orbit + double n = 2. * M_PI / P; + + // brandt solver calculations, once per orbit + double bounds[13]; + double EA_tab[6 * 13]; + get_bounds(bounds, EA_tab, ecc); + + for (size_t i = 0; i < t.size(); i++) + { + + double sinE, cosE; + double M = n * (t[i] - M0_epoch) + M0; + + solver_fixed_ecc(bounds, EA_tab, M, ecc, &sinE, &cosE); + + double X = cosE - ecc; + double Y = sinE * sqrt(1-ecc*ecc); + + double sinpsi, cospsi; + sincos(psi[i], &sinpsi, &cospsi); + // std::cout << M << '\t' << ecc << '\t' << sinE << '\t' << cosE << std::endl; + // std::cout << '\t' << g << '\t' << g2 << std::endl; + wk[i] = (B*X + G*Y)*sinpsi + (A*X + F*Y)*cospsi; + } + + return wk; + + } std::vector keplerian2(const std::vector &t, const double &P, const double &K, const double &ecc, diff --git a/src/kima/kepler.h b/src/kima/kepler.h index 09c16e5..bb47ce9 100644 --- a/src/kima/kepler.h +++ b/src/kima/kepler.h @@ -82,6 +82,10 @@ namespace brandt const double &K, const double &ecc, const double &w, const double &wdot, const double &M0, const double &M0_epoch); + std::vector keplerian_gaia(const std::vector &t, const std::vector &psi, const double &A, + const double &B, const double &F, const double &G, + const double &ecc, const double P, const double &M0, + const double &M0_epoch); } diff --git a/src/kima/postkepler.cpp b/src/kima/postkepler.cpp index 1db2452..7cd04de 100644 --- a/src/kima/postkepler.cpp +++ b/src/kima/postkepler.cpp @@ -96,19 +96,19 @@ namespace postKep return delta_LT; } - inline double transverse_doppler(double K1, double f, double ecc) + inline double transverse_doppler(double K1, double f, double ecc, double cosi=0) { //here assume inclination of 90 (eclipsing) so sin(i) = 1 - double sini = 1.0; + double sini = 1.0 - cosi*cosi; double delta_TD = pow(K1,2.0)*(1 + ecc*cos(f) - (1-pow(ecc,2.0))/2)/(c_light*pow(sini,2.0)); return delta_TD; } - inline double gravitational_redshift(double K1, double K2, double f, double ecc) + inline double gravitational_redshift(double K1, double K2, double f, double ecc, double cosi=0) { //again assume inclination of 90 (eclipsing) so sin(i) = 1 - double sini = 1.0; + double sini = 1.0 - cosi*cosi; double delta_GR = K1*(K1+K2)*(1+ecc*cos(f))/(c_light*pow(sini,2.0)); return delta_GR; diff --git a/src/kima/postkepler.h b/src/kima/postkepler.h index 99f9aef..88dcf54 100644 --- a/src/kima/postkepler.h +++ b/src/kima/postkepler.h @@ -13,8 +13,8 @@ namespace postKep 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); - inline double gravitational_redshift(double K1, double K2, double f, 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); } \ No newline at end of file