ETVmodel M0_epoch publish

kima-org · Jul 30, 2024 · 977cfdd · 977cfdd
1 parent 00e52e3
commit 977cfdd
Show file tree

Hide file tree

Showing 3 changed files with 43 additions and 20 deletions.
diff --git a/src/kima/Data.cpp b/src/kima/Data.cpp
@@ -1120,6 +1120,8 @@ time  vrad  error  quant  error
         .def_prop_ro("et", [](ETVData &d) { return d.get_et(); }, "The observed mid-eclipse times")
         .def_prop_ro("etsig", [](ETVData &d) { return d.get_etsig(); }, "The uncertainties in the eclipse times");
 
+        .def_rw("M0_epoch", &ETVData::M0_epoch, "reference epoch for the mean anomaly")
+
         //
         //.def("load", &GAIAData::load, "filename"_a, "units"_a, "skip"_a, "max_rows"_a, "delimiter"_a)
 }
diff --git a/src/kima/ETVmodel.cpp b/src/kima/ETVmodel.cpp
@@ -166,10 +166,9 @@ void ETVmodel::remove_known_object()
     auto epochs = data.get_epochs();
     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);
-        for (size_t i = 0; i < data.N(); i++)
-        {
-            mu[i] -= tau[i] / (24 * 3600);
+        auto tau = brandt::keplerian_et(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);
         }
     }
 }
@@ -182,9 +181,8 @@ void ETVmodel::add_known_object()
     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);
-        for (size_t i = 0; i < data.N(); i++)
-        {
-            mu[i] += tau[i] / (24 * 3600);
+        for (size_t i = 0; i < data.N(); i++) {
+            mu[i] += tau[i]/(24*3600);
         }
     }
 }

diff --git a/src/kima/kepler.cpp b/src/kima/kepler.cpp
@@ -1185,17 +1185,40 @@ NB_MODULE(kepler, m) {
     m.def("keplerian2", &brandt::keplerian,
           "t"_a, "P"_a, "K"_a, "ecc"_a, "w"_a, "M0"_a, "M0_epoch"_a,
           KEPLERIAN_DOC);
-
-    m.def("keplerian", [](const std::vector<double> &t, 
-                           const double &P, const double &K, const double &ecc, 
-                           const double &w, const double &M0, const double &M0_epoch) {
-        size_t size = t.size();
-        struct Temp { std::vector<double> v; };
-        Temp *temp = new Temp();
-        temp->v = brandt::keplerian(t, P, K, ecc, w, M0, M0_epoch);
-        nb::capsule owner(temp, [](void *p) noexcept { delete (Temp *) p; });
-        return nb::ndarray<nb::numpy, double>(temp->v.data(), {size}, owner);
-    }, "t"_a, "P"_a, "K"_a, "ecc"_a, "w"_a, "M0"_a, "M0_epoch"_a, KEPLERIAN_DOC
-    // nb::rv_policy::copy
-    );
+    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,
+          "epochs"_a, "P"_a, "K"_a, "ecc"_a, "w"_a, "M0"_a, "ephem1"_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) 
+        //      {
+        //         using array_type = nb::ndarray<nb::numpy, double, nb::shape<1, nb::any>>;
+        //         std::vector<double> rv(t.size());
+        //         // 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));
+
+        //         // brandt solver calculations, once per orbit
+        //         double bounds[13];
+        //         double EA_tab[6 * 13];
+        //         brandt::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;
+        //             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)));
+        //         }
+        //         return rv;
+        //         // return nb::ndarray<nb::numpy, double, nb::shape<nb::any>>(rv.data(), 1, shape);
+        //   }, "keplerian function");
 }