From e88991933c226eb7a27a822b6f73dc4aff24c1e0 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Mon, 26 Aug 2024 17:05:38 -0500 Subject: [PATCH 01/19] earth GDL; continuous events --- grss/version.txt | 2 +- include/gr15.h | 25 +++++-- include/grss.h | 3 +- include/simulation.h | 53 ++++++++++----- src/force.cpp | 36 ++++++++++ src/gr15.cpp | 137 ++++++++++++++++++++++++++----------- src/grss.cpp | 61 +++++++++-------- src/observe.cpp | 159 ++++++++++++++++++++++++++----------------- src/simulation.cpp | 135 +++++++++++++++++++++++++++--------- 9 files changed, 421 insertions(+), 190 deletions(-) diff --git a/grss/version.txt b/grss/version.txt index ef8d7569..d87edbfc 100644 --- a/grss/version.txt +++ b/grss/version.txt @@ -1 +1 @@ -4.2.0 \ No newline at end of file +4.2.1 \ No newline at end of file diff --git a/include/gr15.h b/include/gr15.h index 68846328..19aedd54 100644 --- a/include/gr15.h +++ b/include/gr15.h @@ -79,11 +79,28 @@ void refine_b(std::vector &b, std::vector &e, const real &dtRatio, const size_t &dim); /** - * @brief Check if any events need to be applied after the current timestep. + * @brief Check if any impulsive events need to be applied after the current timestep. */ -void check_and_apply_events(PropSimulation *propSim, const real &t, - real &tNextEvent, size_t &nextEventIdx, - std::vector &xInteg); +void check_and_apply_impulsive_events(PropSimulation *propSim, const real &t, + std::vector &xInteg); + +/** + * @brief Check if any continuous events need to be applied after the current timestep. + */ +void check_continuous_events(PropSimulation *propSim, const real &t); + +/** + * @brief Top level function to check for events after the current timestep. + */ +void check_events(PropSimulation *propSim, const real &t, std::vector &xInteg); + +/** + * @brief Check whether timestep is too large that it would skip an event. + * + * @param[in] propSim PropSimulation object for the integration. + * @param[in] dt Current timestep. + */ +void event_timestep_check(PropSimulation *propSim, real &dt); /** * @brief 15th-order Gauss-Radau integrator for the PropSimulation. diff --git a/include/grss.h b/include/grss.h index 5809aa96..960082dd 100644 --- a/include/grss.h +++ b/include/grss.h @@ -20,7 +20,8 @@ void PropSimulation::integrate() { // flip vectors if integration is backwards in time if (this->integParams.t0 > this->integParams.tf) { - std::reverse(this->events.begin(), this->events.end()); + std::reverse(this->eventMngr.impulsiveEvents.begin(), this->eventMngr.impulsiveEvents.end()); + std::reverse(this->eventMngr.continuousEvents.begin(), this->eventMngr.continuousEvents.end()); std::reverse(this->xObserver.begin(), this->xObserver.end()); std::reverse(this->observerInfo.begin(), this->observerInfo.end()); std::reverse(this->tEval.begin(), this->tEval.end()); diff --git a/include/simulation.h b/include/simulation.h index 472f6eb1..e3e40d23 100644 --- a/include/simulation.h +++ b/include/simulation.h @@ -222,6 +222,13 @@ class IntegBody : public Body { * @param t Time of the event. * @param bodyName Name of the IntegBody for the event * @param bodyIndex Index of the IntegBody in the PropSimulation. + * @param xIntegIndex Starting index of the body's state in the flattened state vector. + * @param deltaV Delta-V for the event. + * @param multiplier Multiplier for the Delta-V. + * @param dt Time duration for the continuous event. + * @param threshold Threshold for the continuous impulse to reach in time tEvent+dt. + * @param isContinuous Flag to indicate if the event is continuous. + * @param isHappening Flag to indicate if the continuous event is happening. */ class Event { private: @@ -229,27 +236,36 @@ class Event { real t; std::string bodyName; size_t bodyIndex; -}; + size_t xIntegIndex; -/** - * @brief Class for impulse events in a PropSimulation. - * - * @param deltaV Delta-V for the impulse. - * @param multiplier Multiplier for the Delta-V. - */ -class ImpulseEvent : public Event { - private: - public: + // for impulsive events std::vector deltaV = {0.0L, 0.0L, 0.0L}; real multiplier = 1.0L; + + // for continuous ejecta events + real dt = 1.0L; + real threshold = 0.999L; + bool isContinuous = false; + bool isHappening = false; + real c = 3.8L; /** * @brief Apply the impulse event to the body. - * - * @param[in] t Time of the event. - * @param[inout] xInteg State of the body. - * @param[in] propDir Direction of propagation. */ - void apply(const real &t, std::vector &xInteg, const real &propDir); + void apply_impulsive(const real &t, std::vector &xInteg, const real &propDir); +}; + +class EventManager { + private: + public: + std::vector impulsiveEvents = {}; + std::vector continuousEvents = {}; + size_t nextImpEventIdx = 0; + size_t nextConEventIdx = 0; + real tNextImpEvent = std::numeric_limits::quiet_NaN(); + real tNextConEvent = std::numeric_limits::quiet_NaN(); + size_t nImpEvents = 0; + size_t nConEvents = 0; + bool allConEventDone = true; }; /** @@ -461,7 +477,7 @@ class PropSimulation { IntegrationParameters integParams; std::vector spiceBodies; std::vector integBodies; - std::vector events; + EventManager eventMngr; std::vector caParams; std::vector impactParams; real t; @@ -506,6 +522,11 @@ class PropSimulation { */ void add_event(IntegBody body, real tEvent, std::vector deltaV, real multiplier = 1.0L); + /** + * @brief Add an ejecta event to the simulation. + */ + void add_event(IntegBody body, real tEvent, std::vector deltaV, + real multiplier = 1.0L, real dt = 1.0L, real threshold = 0.999L); /** * @brief Set the values of the PropSimulation Constants object. */ diff --git a/src/force.cpp b/src/force.cpp index 3004a4fd..8179b740 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -37,6 +37,12 @@ static void force_nongrav(const PropSimulation *propSim, std::vector &accI */ static void force_thruster(const PropSimulation *propSim, std::vector &accInteg); +/** + * @brief Compute the acceleration of the system due to a continuous event. + */ +static void force_continuous_event(const real &t, const PropSimulation *propSim, + std::vector &accInteg); + /** * @param[in] t Time [TDB MJD] * @param[in] xInteg State vector @@ -115,6 +121,7 @@ void get_state_der(PropSimulation *propSim, const real &t, force_J2(propSim, accInteg, allSTMs); force_nongrav(propSim, accInteg, allSTMs); force_thruster(propSim, accInteg); + force_continuous_event(t, propSim, accInteg); #ifdef PRINT_FORCES forceFile.open("cpp.11", std::ios::app); forceFile << std::setw(10) << "total_acc" << std::setw(25) << accInteg[0] @@ -636,3 +643,32 @@ static void force_thruster(const PropSimulation *propSim, forceFile.close(); #endif } + +/** + * @param[in] propSim PropSimulation object. + * @param[inout] accInteg State derivative vector. + */ +static void force_continuous_event(const real &t, const PropSimulation *propSim, + std::vector &accInteg) { + if (!propSim->eventMngr.allConEventDone){ + const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; + const real propDir = forwardProp ? 1.0 : -1.0; + for (size_t i = 0; i < propSim->eventMngr.continuousEvents.size(); i++){ + if (propSim->eventMngr.continuousEvents[i].isHappening){ + const Event event = propSim->eventMngr.continuousEvents[i]; + const size_t starti = event.xIntegIndex / 2; + const real tPastEvent = t - event.t; + // f = np.sin(np.pi*x/(2*time_for_val)) + // df = np.pi/(2*time_for_val)*np.cos(np.pi*x/(2*time_for_val)) + const real preFac = PI/(2*event.dt); + real accFac = preFac*cos(preFac*tPastEvent); + if (accFac > preFac || accFac < 0.0){ + accFac = 0.0; + } + accInteg[starti + 0] += accFac * event.deltaV[0] * propDir; + accInteg[starti + 1] += accFac * event.deltaV[1] * propDir; + accInteg[starti + 2] += accFac * event.deltaV[2] * propDir; + } + } + } +} diff --git a/src/gr15.cpp b/src/gr15.cpp index ad6213de..0805d524 100644 --- a/src/gr15.cpp +++ b/src/gr15.cpp @@ -225,30 +225,104 @@ void refine_b(std::vector &b, std::vector &e, const real &dtRatio, c /** * @param[in] propSim PropSimulation object for the integration. * @param[in] t Current time. - * @param[out] tNextEvent Time of the next event. - * @param[out] nextEventIdx Index of the next event. * @param[in] xInteg Vector of integrated states after the current timestep. */ -void check_and_apply_events(PropSimulation *propSim, const real &t, - real &tNextEvent, size_t &nextEventIdx, - std::vector &xInteg) { - while (nextEventIdx < propSim->events.size() && t == tNextEvent) { +void check_and_apply_impulsive_events(PropSimulation *propSim, const real &t, + std::vector &xInteg) { + size_t *nextImpEventIdx = &propSim->eventMngr.nextImpEventIdx; + real *tNextImpEvent = &propSim->eventMngr.tNextImpEvent; + while (*nextImpEventIdx < propSim->eventMngr.nImpEvents && t == *tNextImpEvent) { // apply events for the state just reached by the integrator - real propDir; - if (propSim->integParams.t0 < propSim->integParams.tf) { - propDir = 1.0L; + bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; + real propDir = forwardProp ? 1.0L : -1.0L; + propSim->eventMngr.impulsiveEvents[*nextImpEventIdx].apply_impulsive(t, xInteg, propDir); + // update next event index and time + (*nextImpEventIdx)++; + if (*nextImpEventIdx < propSim->eventMngr.nImpEvents) { + *tNextImpEvent = propSim->eventMngr.impulsiveEvents[*nextImpEventIdx].t; } else { - propDir = -1.0L; + *tNextImpEvent = propSim->integParams.tf; } - propSim->events[nextEventIdx].apply(t, xInteg, propDir); - // update next event index and time - nextEventIdx++; - if (nextEventIdx < propSim->events.size()) { - tNextEvent = propSim->events[nextEventIdx].t; + } +} + +/** + * @param[in] propSim PropSimulation object for the integration. + * @param[in] t Current time. + */ +void check_continuous_events(PropSimulation *propSim, const real &t) { + if (!propSim->eventMngr.allConEventDone){ + bool allDone = true; + bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; + for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) { + // const bool eventStarted = t >= propSim->eventMngr.continuousEvents[i].t; + // const bool eventEnded = t > propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt; + bool eventStarted, eventEnded; + if (forwardProp) { + eventStarted = t >= propSim->eventMngr.continuousEvents[i].t; + eventEnded = t >= (propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt); + } else { + eventStarted = t <= (propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt); + eventEnded = t < propSim->eventMngr.continuousEvents[i].t; + } + if (eventStarted && !eventEnded) { + propSim->eventMngr.continuousEvents[i].isHappening = true; + } else { + propSim->eventMngr.continuousEvents[i].isHappening = false; + } + if (!eventEnded) { + allDone = false; + } + } + if (allDone) { + propSim->eventMngr.allConEventDone = true; + propSim->eventMngr.tNextConEvent = propSim->integParams.tf; + } + } +} + +/** + * @param[in] propSim PropSimulation object for the integration. + * @param[in] t Current time. + * @param[in] xInteg Vector of integrated states after the current timestep. + */ +void check_events(PropSimulation *propSim, const real &t, std::vector &xInteg){ + check_and_apply_impulsive_events(propSim, t, xInteg); + check_continuous_events(propSim, t); +} + +/** + * @param[in] propSim PropSimulation object for the integration. + */ +void event_timestep_check(PropSimulation *propSim, real &dt) { + real tNextEvent = propSim->integParams.tf; + const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; + if (propSim->eventMngr.nImpEvents > 0) { + const real tNextImpEvent = propSim->eventMngr.tNextImpEvent; + if (forwardProp) { + tNextEvent = fmin(tNextEvent, tNextImpEvent); } else { - tNextEvent = propSim->integParams.tf; + tNextEvent = fmax(tNextEvent, tNextImpEvent); + } + } + if (!propSim->eventMngr.allConEventDone) { + for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) { + if (propSim->eventMngr.continuousEvents[i].isHappening) { + // if any continuous event is happening, set dt to 1/10th of the event duration + dt = propSim->eventMngr.continuousEvents[i].dt/10.0L; + dt = forwardProp ? dt : -dt; + } + if (forwardProp && propSim->t < propSim->eventMngr.continuousEvents[i].t) { + tNextEvent = fmin(tNextEvent, propSim->eventMngr.continuousEvents[i].t); + } else if (!forwardProp && propSim->t > propSim->eventMngr.continuousEvents[i].t+propSim->eventMngr.continuousEvents[i].dt) { + tNextEvent = fmax(tNextEvent, propSim->eventMngr.continuousEvents[i].t+propSim->eventMngr.continuousEvents[i].dt); + } } } + if ((forwardProp && propSim->t + dt > tNextEvent) || + (!forwardProp && propSim->t + dt < tNextEvent)) { + dt = tNextEvent - propSim->t; + } } /** @@ -342,21 +416,8 @@ void gr15(PropSimulation *propSim) { std::vector g(7 * dim, 0.0); std::vector e(7 * dim, 0.0); real dtReq; - real tNextEvent = propSim->integParams.tf; - size_t nextEventIdx = 0; - if (t == propSim->integParams.t0) { - nextEventIdx = 0; - } - if (propSim->events.size() != 0) { - tNextEvent = propSim->events[0].t; - } - check_and_apply_events(propSim, t, tNextEvent, nextEventIdx, xInteg); - if ((propSim->integParams.tf > propSim->integParams.t0 && - t + dt > tNextEvent) || - (propSim->integParams.tf < propSim->integParams.t0 && - t + dt < tNextEvent)) { - dt = tNextEvent - t; - } + check_events(propSim, t, xInteg); + event_timestep_check(propSim, dt); propSim->interpParams.tStack.push_back(t); propSim->interpParams.xIntegStack.push_back(xInteg); std::vector accInteg0(dim, 0.0); @@ -430,11 +491,12 @@ void gr15(PropSimulation *propSim) { approx_xInteg(xInteg0, accInteg0, dt, 1.0, b, dim, propSim->integBodies, xInteg, xIntegCompCoeffs); t += dt; + check_events(propSim, t, xInteg); get_state_der(propSim, t, xInteg, accInteg0); - check_and_apply_events(propSim, t, tNextEvent, nextEventIdx, - xInteg); propSim->interpParams.tStack.push_back(t); propSim->interpParams.xIntegStack.push_back(xInteg); + propSim->t = t; + propSim->xInteg = xInteg; propSim->integParams.timestepCounter++; refine_b(b, e, dtReq/dt, dim); check_ca_or_impact(propSim, t-dt, xInteg0, t, xInteg); @@ -445,14 +507,7 @@ void gr15(PropSimulation *propSim) { keepStepping = 0; } dt = dtReq; - if ((propSim->integParams.tf > propSim->integParams.t0 && - t + dt > tNextEvent) || - (propSim->integParams.tf < propSim->integParams.t0 && - t + dt < tNextEvent)) { - dt = tNextEvent - t; - } - propSim->t = t; - propSim->xInteg = xInteg; + event_timestep_check(propSim, dt); oneStepDone = 1; } } diff --git a/src/grss.cpp b/src/grss.cpp index b1f54fa9..53660a32 100644 --- a/src/grss.cpp +++ b/src/grss.cpp @@ -636,31 +636,6 @@ PYBIND11_MODULE(libgrss, m) { Reconstructed state transition matrix. )mydelimiter"); - py::class_(m, "Event", R"mydelimiter( - The Event class contains the properties of an integration event. - )mydelimiter") - .def(py::init<>()) - .def_readwrite("t", &ImpulseEvent::t, R"mydelimiter( - MJD TDB Time of the event. - )mydelimiter") - .def_readwrite("bodyName", &ImpulseEvent::bodyName, R"mydelimiter( - Name of the integration body to apply the event to. - )mydelimiter") - .def_readwrite("bodyIndex", &ImpulseEvent::bodyIndex, R"mydelimiter( - Index of the integration body to apply the event to. - )mydelimiter"); - - py::class_(m, "ImpulseEvent", R"mydelimiter( - The ImpulseEvent class contains the properties of an impulsive delta-V event. - )mydelimiter") - .def(py::init<>()) - .def_readwrite("deltaV", &ImpulseEvent::deltaV, R"mydelimiter( - Delta-V of the event. - )mydelimiter") - .def_readwrite("multiplier", &ImpulseEvent::multiplier, R"mydelimiter( - Multiplier on the delta-V the event. - )mydelimiter"); - m.def("propSim_parallel_omp", &propSim_parallel_omp, py::arg("refSim"), py::arg("isCometary"), py::arg("allBodies"), py::arg("maxThreads") = 128, R"mydelimiter( @@ -734,9 +709,6 @@ PYBIND11_MODULE(libgrss, m) { R"mydelimiter( Integration bodies of the simulation. List of PropSimulation.IntegBody objects. )mydelimiter") - .def_readwrite("events", &PropSimulation::events, R"mydelimiter( - Events of the simulation. List of PropSimulation.Event objects. - )mydelimiter") .def_readwrite("caParams", &PropSimulation::caParams, R"mydelimiter( Close approach parameters of the simulation. List of PropSimulation.CloseApproachParameters objects. )mydelimiter") @@ -888,8 +860,12 @@ PYBIND11_MODULE(libgrss, m) { name : str Name of the body to remove. )mydelimiter") - .def("add_event", &PropSimulation::add_event, py::arg("body"), - py::arg("tEvent"), py::arg("deltaV"), py::arg("multiplier") = 1.0L, + .def("add_event", + static_cast, real)>( + &PropSimulation::add_event), + py::arg("body"), py::arg("tEvent"), py::arg("deltaV"), + py::arg("multiplier") = 1.0L, R"mydelimiter( Adds an impulsive delta-V event to the simulation. @@ -904,6 +880,31 @@ PYBIND11_MODULE(libgrss, m) { multiplier : real Multiplier to apply to the delta-V. )mydelimiter") + .def("add_event", + static_cast, real, real, real)>( + &PropSimulation::add_event), + py::arg("bodyName"), py::arg("tEvent"), py::arg("deltaV"), + py::arg("multiplier") = 1.0L, py::arg("dt") = 1.0L, + py::arg("threshold") = 0.999L, + R"mydelimiter( + Adds an ejecta event to the simulation. + + Parameters + ---------- + body : str + Name of the body to apply the delta-V to. + tEvent : real + MJD Epoch of the event. Must be in TDB. + deltaV : list of real + Delta-V to apply to the body. + multiplier : real + Multiplier to apply to the delta-V. + dt : real + Time duration for the continuous ejecta event. + threshold : real + Threshold for the ejecta impulse to reach in time tEvent+dt. + )mydelimiter") .def("set_sim_constants", &PropSimulation::set_sim_constants, py::arg("du2m") = 149597870700.0L, py::arg("tu2s") = 86400.0L, py::arg("G") = 6.6743e-11L / diff --git a/src/observe.cpp b/src/observe.cpp index dbc94a84..730d4990 100644 --- a/src/observe.cpp +++ b/src/observe.cpp @@ -9,85 +9,116 @@ void get_glb_correction(PropSimulation *propSim, const size_t &interpIdx, const real &tInterpGeom, std::vector &xInterpApparentBary) { - double sunState[9]; - get_spk_state(10, tInterpGeom, propSim->spkEphem, sunState); - std::vector earthState = propSim->xObserver[interpIdx]; - - std::vector sunEarthPos = {earthState[0] - sunState[0], - earthState[1] - sunState[1], - earthState[2] - sunState[2]}; - real sunEarthDist; - vnorm(sunEarthPos, sunEarthDist); - std::vector sunTargetPos = {xInterpApparentBary[0] - sunState[0], - xInterpApparentBary[1] - sunState[1], - xInterpApparentBary[2] - sunState[2]}; - real sunTargetDist; - vnorm(sunTargetPos, sunTargetDist); - std::vector earthTargetPos = {xInterpApparentBary[0] - earthState[0], - xInterpApparentBary[1] - earthState[1], - xInterpApparentBary[2] - earthState[2]}; - real earthTargetDist; - vnorm(earthTargetPos, earthTargetDist); - - real G = propSim->consts.G; + const real G = propSim->consts.G; real sunGM = 0; + real earthGM = 0; for (size_t i = 0; i < propSim->integParams.nSpice; i++) { if (propSim->spiceBodies[i].spiceId == 10) { sunGM = G * propSim->spiceBodies[i].mass; } + if (propSim->spiceBodies[i].spiceId == 399) { + earthGM = G * propSim->spiceBodies[i].mass; + } } - if (sunGM == 0) { + if (sunGM == 0 || earthGM == 0) { throw std::runtime_error( - "Sun GM not found in get_delta_delay_relativistic"); + "Sun GM or Earth GM not found in get_delta_delay_relativistic"); } - real c = propSim->consts.clight; + const real c = propSim->consts.clight; - // from section 7.2.4 in the Explanatory Supplement to the Astronomical - // Almanac, 3rd edition - std::vector e(3, 0.0); - vunit(sunEarthPos, e); - std::vector q(3, 0.0); - vunit(sunTargetPos, q); + double sunState[9]; + double earthState[9]; + get_spk_state(10, tInterpGeom, propSim->spkEphem, sunState); + get_spk_state(399, tInterpGeom, propSim->spkEphem, earthState); + std::vector observerState = propSim->xObserver[interpIdx]; + std::vector observerTargetPos = {xInterpApparentBary[0] - observerState[0], + xInterpApparentBary[1] - observerState[1], + xInterpApparentBary[2] - observerState[2]}; + real observerTargetDist; + vnorm(observerTargetPos, observerTargetDist); std::vector p(3, 0.0); - vunit(earthTargetPos, p); + vunit(observerTargetPos, p); + std::vector deltaP1Targ(3, 0.0); std::vector deltaP1Star(3, 0.0); - std::vector p1(3, 0.0); - - real pDotQ, eDotP, qDotE; - vdot(p, q, pDotQ); - vdot(e, p, eDotP); - vdot(q, e, qDotE); - - real g1 = 2 * sunGM / c / c / sunEarthDist; - real g2 = 1.0L + qDotE; - - deltaP1Targ[0] = g1 * (pDotQ * e[0] - eDotP * q[0]) / g2; - deltaP1Targ[1] = g1 * (pDotQ * e[1] - eDotP * q[1]) / g2; - deltaP1Targ[2] = g1 * (pDotQ * e[2] - eDotP * q[2]) / g2; - - deltaP1Star[0] = g1 * (e[0] - eDotP * p[0]) / (1 + eDotP); - deltaP1Star[1] = g1 * (e[1] - eDotP * p[1]) / (1 + eDotP); - deltaP1Star[2] = g1 * (e[2] - eDotP * p[2]) / (1 + eDotP); - - // do absolute correction first - p1[0] = p[0] + deltaP1Targ[0]; - p1[1] = p[1] + deltaP1Targ[1]; - p1[2] = p[2] + deltaP1Targ[2]; - // if not gaia obs, do full relative correction - if (propSim->obsType[interpIdx] != 3) { - p1[0] -= deltaP1Star[0]; - p1[1] -= deltaP1Star[1]; - p1[2] -= deltaP1Star[2]; + std::vector p1(p); + + size_t numBendingBodies = 1; // sun (1) or sun+earth (2) + for (size_t i = 0; i < numBendingBodies; i++){ + double centralBodyState[6]; + real centralBodyGM; + switch (i) { + case 0: + for (size_t j = 0; j < 6; j++) { + centralBodyState[j] = sunState[j]; + } + centralBodyGM = sunGM; + break; + case 1: + for (size_t j = 0; j < 6; j++) { + centralBodyState[j] = earthState[j]; + } + centralBodyGM = earthGM; + break; + default: + throw std::runtime_error( + "get_glb_correction: central body index must be 0 or 1"); + break; + } + + std::vector centralBodyObserverPos = {observerState[0] - centralBodyState[0], + observerState[1] - centralBodyState[1], + observerState[2] - centralBodyState[2]}; + real centralBodyObserverDist; + vnorm(centralBodyObserverPos, centralBodyObserverDist); + std::vector centralBodyTargetPos = {xInterpApparentBary[0] - centralBodyState[0], + xInterpApparentBary[1] - centralBodyState[1], + xInterpApparentBary[2] - centralBodyState[2]}; + real centralBodyTargetDist; + vnorm(centralBodyTargetPos, centralBodyTargetDist); + + // from section 7.2.4 in the Explanatory Supplement to the Astronomical + // Almanac, 3rd edition + std::vector e(3, 0.0); + vunit(centralBodyObserverPos, e); + std::vector q(3, 0.0); + vunit(centralBodyTargetPos, q); + + real pDotQ, eDotP, qDotE; + vdot(p, q, pDotQ); + vdot(e, p, eDotP); + vdot(q, e, qDotE); + + real g1 = 2 * centralBodyGM / c / c / centralBodyObserverDist; + real g2 = 1.0L + qDotE; + + deltaP1Targ[0] += g1 * (pDotQ * e[0] - eDotP * q[0]) / g2; + deltaP1Targ[1] += g1 * (pDotQ * e[1] - eDotP * q[1]) / g2; + deltaP1Targ[2] += g1 * (pDotQ * e[2] - eDotP * q[2]) / g2; + + deltaP1Star[0] += g1 * (e[0] - eDotP * p[0]) / (1 + eDotP); + deltaP1Star[1] += g1 * (e[1] - eDotP * p[1]) / (1 + eDotP); + deltaP1Star[2] += g1 * (e[2] - eDotP * p[2]) / (1 + eDotP); + + // do absolute correction first + p1[0] += deltaP1Targ[0]; + p1[1] += deltaP1Targ[1]; + p1[2] += deltaP1Targ[2]; + // if not gaia obs, do full relative correction + if (propSim->obsType[interpIdx] != 3) { + p1[0] -= deltaP1Star[0]; + p1[1] -= deltaP1Star[1]; + p1[2] -= deltaP1Star[2]; + } } - earthTargetPos[0] = earthTargetDist * p1[0]; - earthTargetPos[1] = earthTargetDist * p1[1]; - earthTargetPos[2] = earthTargetDist * p1[2]; + observerTargetPos[0] = observerTargetDist * p1[0]; + observerTargetPos[1] = observerTargetDist * p1[1]; + observerTargetPos[2] = observerTargetDist * p1[2]; - xInterpApparentBary[0] = earthState[0] + earthTargetPos[0]; - xInterpApparentBary[1] = earthState[1] + earthTargetPos[1]; - xInterpApparentBary[2] = earthState[2] + earthTargetPos[2]; + xInterpApparentBary[0] = observerState[0] + observerTargetPos[0]; + xInterpApparentBary[1] = observerState[1] + observerTargetPos[1]; + xInterpApparentBary[2] = observerState[2] + observerTargetPos[2]; } /** diff --git a/src/simulation.cpp b/src/simulation.cpp index b015a212..f432261e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -339,14 +339,14 @@ void IntegBody::prepare_stm(){ * @param[inout] xInteg State of the body. * @param[in] propDir Direction of propagation. */ -void ImpulseEvent::apply(const real& t, std::vector& xInteg, +void Event::apply_impulsive(const real& t, std::vector& xInteg, const real& propDir) { if (t != this->t) { throw std::runtime_error( - "ImpulseEvent::apply: Integration time does " + "Event::apply_impulsive: Integration time does " "not match event time. Cannot apply impulse."); } - size_t velStartIdx = 6 * this->bodyIndex + 3; + size_t velStartIdx = this->xIntegIndex + 3; for (size_t i = 0; i < 3; i++) { xInteg[velStartIdx + i] += propDir * this->multiplier * this->deltaV[i]; } @@ -985,29 +985,23 @@ void PropSimulation::remove_body(std::string name) { std::cout << "Error: Body " << name << " not found." << std::endl; } -/** - * @param[in] body IntegBody object to apply the ImpulseEvent to. - * @param[in] tEvent Time at which to apply the ImpulseEvent. - * @param[in] deltaV Delta-V for the impulse. - * @param[in] multiplier Multiplier for the Delta-V. - */ -void PropSimulation::add_event(IntegBody body, real tEvent, - std::vector deltaV, real multiplier) { +size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, + const real &tEvent) { // check if tEvent is valid - const bool forwardProp = this->integParams.tf > this->integParams.t0; - const bool backwardProp = this->integParams.tf < this->integParams.t0; + const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; + const bool backwardProp = propSim->integParams.tf < propSim->integParams.t0; if ((forwardProp && - (tEvent < this->integParams.t0 || tEvent >= this->integParams.tf)) || + (tEvent < propSim->integParams.t0 || tEvent >= propSim->integParams.tf)) || (backwardProp && - (tEvent > this->integParams.t0 || tEvent <= this->integParams.tf))) { + (tEvent > propSim->integParams.t0 || tEvent <= propSim->integParams.tf))) { throw std::invalid_argument("Event time " + std::to_string(tEvent) + " is not within simulation time bounds."); } // check if body exists bool bodyExists = false; size_t bodyIndex; - for (size_t i = 0; i < this->integParams.nInteg; i++) { - if (this->integBodies[i].name == body.name) { + for (size_t i = 0; i < propSim->integParams.nInteg; i++) { + if (propSim->integBodies[i].name == body.name) { bodyExists = true; bodyIndex = i; break; @@ -1016,30 +1010,93 @@ void PropSimulation::add_event(IntegBody body, real tEvent, if (!bodyExists) { throw std::invalid_argument("Integration body with name " + body.name + " does not exist in simulation " + - this->name); + propSim->name); } - ImpulseEvent event; - event.t = tEvent; - event.deltaV = deltaV; - event.multiplier = multiplier; - event.bodyName = body.name; - event.bodyIndex = bodyIndex; - // add event to this->events sorted by time - if (this->events.size() == 0) { - this->events.push_back(event); + return bodyIndex; +} + +void event_postprocess(PropSimulation *propSim, Event &event) { + // assign event xIntegIndex + event.xIntegIndex = 0; + for (size_t i = 0; i < event.bodyIndex; i++) { + event.xIntegIndex += 2*propSim->integBodies[i].n2Derivs; + } + // add event to either continuous or impulse events based on isContinuous + std::vector *eventsList; + if (event.isContinuous) { + eventsList = &propSim->eventMngr.continuousEvents; + } else { + eventsList = &propSim->eventMngr.impulsiveEvents; + } + if (eventsList->size() == 0) { + eventsList->push_back(event); } else { - for (size_t i = 0; i < this->events.size(); i++) { - if (event.t < this->events[i].t) { - this->events.insert(this->events.begin() + i, event); + for (size_t i = 0; i < eventsList->size(); i++) { + if (event.t < eventsList->at(i).t) { + eventsList->insert(eventsList->begin() + i, event); break; - } else if (i == this->events.size() - 1) { - this->events.push_back(event); + } else if (i == eventsList->size() - 1) { + eventsList->push_back(event); break; } } } } +/** + * @param[in] body IntegBody object to apply the ImpulseEvent to. + * @param[in] tEvent Time at which to apply the ImpulseEvent. + * @param[in] deltaV Delta-V for the impulse. + * @param[in] multiplier Multiplier for the Delta-V. + */ +void PropSimulation::add_event(IntegBody body, real tEvent, + std::vector deltaV, real multiplier) { + size_t bodyIndex = event_preprocess(this, body, tEvent); + Event event; + event.t = tEvent; + event.bodyName = body.name; + event.bodyIndex = bodyIndex; + event.deltaV = deltaV; + event.multiplier = multiplier; + + event.dt = 0.0; + event.threshold = 1.0; + event.isContinuous = false; + event.isHappening = false; + // event.c is infinity for impulse events + event.c = std::numeric_limits::infinity(); + event_postprocess(this, event); + this->eventMngr.nImpEvents++; +} + +/** + * @param[in] body IntegBody object to apply the EjectaEvent to. + * @param[in] tEvent Time at which to apply the EjectaEvent. + * @param[in] deltaV Delta-V for the ejecta. + * @param[in] multiplier Multiplier for the Delta-V. + * @param[in] dt Time duration for the ejecta. + * @param[in] threshold Threshold for the ejecta impulse to reach in time tEvent+dt. + */ +void PropSimulation::add_event(IntegBody body, real tEvent, + std::vector deltaV, real multiplier, + real dt, real threshold) { + size_t bodyIndex = event_preprocess(this, body, tEvent); + Event event; + event.t = tEvent; + event.bodyName = body.name; + event.bodyIndex = bodyIndex; + event.deltaV = deltaV; + event.multiplier = multiplier; + + event.dt = dt; + event.threshold = threshold; + event.isContinuous = true; + event.isHappening = false; + event.c = tanh(threshold)/dt; + event_postprocess(this, event); + this->eventMngr.nConEvents++; +} + /** * @param[in] du2m Distance unit conversion factor (default: AU to meters). * @param[in] tu2s Time unit conversion factor (default: days to seconds). @@ -1155,7 +1212,19 @@ void PropSimulation::preprocess() { } bool backwardProp = this->integParams.t0 > this->integParams.tf; if (backwardProp) { - std::reverse(this->events.begin(), this->events.end()); + std::reverse(this->eventMngr.impulsiveEvents.begin(), this->eventMngr.impulsiveEvents.end()); + std::reverse(this->eventMngr.continuousEvents.begin(), this->eventMngr.continuousEvents.end()); + } + this->eventMngr.nextImpEventIdx = 0; + this->eventMngr.nextConEventIdx = 0; + if (this->eventMngr.nImpEvents > 0) { + this->eventMngr.tNextImpEvent = this->eventMngr.impulsiveEvents[0].t; + } + if (this->eventMngr.nConEvents > 0) { + this->eventMngr.tNextConEvent = this->eventMngr.continuousEvents[0].t; + this->eventMngr.allConEventDone = false; + } else { + this->eventMngr.allConEventDone = true; } this->isPreprocessed = true; } From bbbe2dc1ad7e99b3d090e6a4abdbb52bddc67f3b Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Mon, 26 Aug 2024 17:31:27 -0500 Subject: [PATCH 02/19] remove unused threshold for tanh continuous events --- include/simulation.h | 7 ++----- src/grss.cpp | 9 +++------ src/simulation.cpp | 8 +------- 3 files changed, 6 insertions(+), 18 deletions(-) diff --git a/include/simulation.h b/include/simulation.h index e3e40d23..44e1cd35 100644 --- a/include/simulation.h +++ b/include/simulation.h @@ -226,7 +226,6 @@ class IntegBody : public Body { * @param deltaV Delta-V for the event. * @param multiplier Multiplier for the Delta-V. * @param dt Time duration for the continuous event. - * @param threshold Threshold for the continuous impulse to reach in time tEvent+dt. * @param isContinuous Flag to indicate if the event is continuous. * @param isHappening Flag to indicate if the continuous event is happening. */ @@ -244,10 +243,8 @@ class Event { // for continuous ejecta events real dt = 1.0L; - real threshold = 0.999L; bool isContinuous = false; bool isHappening = false; - real c = 3.8L; /** * @brief Apply the impulse event to the body. */ @@ -521,12 +518,12 @@ class PropSimulation { * @brief Add an impulse event to the simulation. */ void add_event(IntegBody body, real tEvent, std::vector deltaV, - real multiplier = 1.0L); + real multiplier); /** * @brief Add an ejecta event to the simulation. */ void add_event(IntegBody body, real tEvent, std::vector deltaV, - real multiplier = 1.0L, real dt = 1.0L, real threshold = 0.999L); + real multiplier, real dt); /** * @brief Set the values of the PropSimulation Constants object. */ diff --git a/src/grss.cpp b/src/grss.cpp index 53660a32..f2bedbe1 100644 --- a/src/grss.cpp +++ b/src/grss.cpp @@ -865,7 +865,7 @@ PYBIND11_MODULE(libgrss, m) { std::vector, real)>( &PropSimulation::add_event), py::arg("body"), py::arg("tEvent"), py::arg("deltaV"), - py::arg("multiplier") = 1.0L, + py::arg("multiplier"), R"mydelimiter( Adds an impulsive delta-V event to the simulation. @@ -882,11 +882,10 @@ PYBIND11_MODULE(libgrss, m) { )mydelimiter") .def("add_event", static_cast, real, real, real)>( + IntegBody, real, std::vector, real, real)>( &PropSimulation::add_event), py::arg("bodyName"), py::arg("tEvent"), py::arg("deltaV"), - py::arg("multiplier") = 1.0L, py::arg("dt") = 1.0L, - py::arg("threshold") = 0.999L, + py::arg("multiplier"), py::arg("dt"), R"mydelimiter( Adds an ejecta event to the simulation. @@ -902,8 +901,6 @@ PYBIND11_MODULE(libgrss, m) { Multiplier to apply to the delta-V. dt : real Time duration for the continuous ejecta event. - threshold : real - Threshold for the ejecta impulse to reach in time tEvent+dt. )mydelimiter") .def("set_sim_constants", &PropSimulation::set_sim_constants, py::arg("du2m") = 149597870700.0L, py::arg("tu2s") = 86400.0L, diff --git a/src/simulation.cpp b/src/simulation.cpp index f432261e..7aa038d6 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1060,11 +1060,8 @@ void PropSimulation::add_event(IntegBody body, real tEvent, event.multiplier = multiplier; event.dt = 0.0; - event.threshold = 1.0; event.isContinuous = false; event.isHappening = false; - // event.c is infinity for impulse events - event.c = std::numeric_limits::infinity(); event_postprocess(this, event); this->eventMngr.nImpEvents++; } @@ -1075,11 +1072,10 @@ void PropSimulation::add_event(IntegBody body, real tEvent, * @param[in] deltaV Delta-V for the ejecta. * @param[in] multiplier Multiplier for the Delta-V. * @param[in] dt Time duration for the ejecta. - * @param[in] threshold Threshold for the ejecta impulse to reach in time tEvent+dt. */ void PropSimulation::add_event(IntegBody body, real tEvent, std::vector deltaV, real multiplier, - real dt, real threshold) { + real dt) { size_t bodyIndex = event_preprocess(this, body, tEvent); Event event; event.t = tEvent; @@ -1089,10 +1085,8 @@ void PropSimulation::add_event(IntegBody body, real tEvent, event.multiplier = multiplier; event.dt = dt; - event.threshold = threshold; event.isContinuous = true; event.isHappening = false; - event.c = tanh(threshold)/dt; event_postprocess(this, event); this->eventMngr.nConEvents++; } From 95f6a63be0a093556c2377434ed6f80645ff64ca Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Tue, 27 Aug 2024 00:23:00 -0500 Subject: [PATCH 03/19] add continuous event estimation --- grss/fit/fit_simulation.py | 24 ++++++++++++++++++--- src/force.cpp | 26 ++++++++++++++-------- src/gr15.cpp | 3 ++- src/interpolate.cpp | 44 ++++++++++++++++---------------------- 4 files changed, 58 insertions(+), 39 deletions(-) diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index f22b5243..b9e1b328 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -987,11 +987,14 @@ def _x_dict_to_events(self, x_dict): for i in range(len(self.fixed_propsim_params['events'])): fixed_event = tuple(self.fixed_propsim_params['events'][i]) event = [None]*5 - event[0] = fixed_event[0] + event[0] = fixed_event[0] # time of delta-v event[1] = x_dict[f"dvx{i}"] if f"dvx{i}" in x_dict.keys() else fixed_event[1] event[2] = x_dict[f"dvy{i}"] if f"dvy{i}" in x_dict.keys() else fixed_event[2] event[3] = x_dict[f"dvz{i}"] if f"dvz{i}" in x_dict.keys() else fixed_event[3] event[4] = x_dict[f"mult{i}"] if f"mult{i}" in x_dict.keys() else fixed_event[4] + fixed_event_continuous = len(fixed_event) == 6 and fixed_event[5] != 0.0 + if fixed_event_continuous or f"dt{i}" in x_dict.keys(): + event.append(x_dict[f"dt{i}"] if f"dt{i}" in x_dict.keys() else fixed_event[5]) events.append(tuple(event)) return events @@ -1023,10 +1026,14 @@ def _check_and_add_events(self, prop_sim_past, prop_sim_future, integ_body, even dvy = event[2] dvz = event[3] multiplier = event[4] + event_args = (integ_body, t_event, [dvx, dvy, dvz], multiplier) + if len(event) == 6: + dt = event[5] + event_args += (dt,) if t_event < self.t_sol: - prop_sim_past.add_event(integ_body, t_event, [dvx, dvy, dvz], multiplier) + prop_sim_past.add_event(*event_args) else: - prop_sim_future.add_event(integ_body, t_event, [dvx, dvy, dvz], multiplier) + prop_sim_future.add_event(*event_args) return prop_sim_past, prop_sim_future def _get_perturbed_state(self, key): @@ -1055,6 +1062,7 @@ def _get_perturbed_state(self, key): fd_delta : float Finite difference perturbation. """ + fd_pert = np.inf if self.fit_cartesian: if key in {'x', 'y', 'z'}: fd_pert = 1e-8 @@ -1074,6 +1082,11 @@ def _get_perturbed_state(self, key): fd_pert = 1e0 if key[:3] in {'dvx', 'dvy', 'dvz'}: fd_pert = 1e-11 + if key[:2] == 'dt': + fd_pert = 1.0 + if fd_pert == np.inf: + raise ValueError("Finite difference perturbation not defined" + " for the given state parameter.") x_plus = self.x_nom.copy() x_minus = self.x_nom.copy() # fd_pert = finite difference perturbation to nominal state for calculating derivatives @@ -1720,6 +1733,11 @@ def filter_lsq(self, verbose=True): print("WARNING: Eccentricity is negative per least squares state correction. " "Setting to 0.0. This solution may not be trustworthy.") self.x_nom = dict(zip(self.x_nom.keys(), next_state)) + # make sure any keys starting with dt are positive and stay in the range [0, 180] + for key in self.x_nom: + if key[:2] == 'dt': + self.x_nom[key] = abs(self.x_nom[key]) + self.x_nom[key] = min(self.x_nom[key], 180.0) # add iteration self._add_iteration(i+1, rms_u, rms_w, chi_sq) if verbose: diff --git a/src/force.cpp b/src/force.cpp index 8179b740..87ab2556 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -655,19 +655,27 @@ static void force_continuous_event(const real &t, const PropSimulation *propSim, const real propDir = forwardProp ? 1.0 : -1.0; for (size_t i = 0; i < propSim->eventMngr.continuousEvents.size(); i++){ if (propSim->eventMngr.continuousEvents[i].isHappening){ - const Event event = propSim->eventMngr.continuousEvents[i]; - const size_t starti = event.xIntegIndex / 2; - const real tPastEvent = t - event.t; + const size_t starti = + propSim->eventMngr.continuousEvents[i].xIntegIndex / 2; + const real tPastEvent = + t - propSim->eventMngr.continuousEvents[i].t; // f = np.sin(np.pi*x/(2*time_for_val)) // df = np.pi/(2*time_for_val)*np.cos(np.pi*x/(2*time_for_val)) - const real preFac = PI/(2*event.dt); - real accFac = preFac*cos(preFac*tPastEvent); - if (accFac > preFac || accFac < 0.0){ + const real preFac = + PI / (2 * propSim->eventMngr.continuousEvents[i].dt); + real accFac = preFac * cos(preFac * tPastEvent); + if (accFac > preFac || accFac < 0.0) { accFac = 0.0; } - accInteg[starti + 0] += accFac * event.deltaV[0] * propDir; - accInteg[starti + 1] += accFac * event.deltaV[1] * propDir; - accInteg[starti + 2] += accFac * event.deltaV[2] * propDir; + accInteg[starti + 0] += accFac * + propSim->eventMngr.continuousEvents[i].deltaV[0] * + propSim->eventMngr.continuousEvents[i].multiplier * propDir; + accInteg[starti + 1] += accFac * + propSim->eventMngr.continuousEvents[i].deltaV[1] * + propSim->eventMngr.continuousEvents[i].multiplier * propDir; + accInteg[starti + 2] += accFac * + propSim->eventMngr.continuousEvents[i].deltaV[2] * + propSim->eventMngr.continuousEvents[i].multiplier * propDir; } } } diff --git a/src/gr15.cpp b/src/gr15.cpp index 0805d524..61a6dd89 100644 --- a/src/gr15.cpp +++ b/src/gr15.cpp @@ -309,7 +309,8 @@ void event_timestep_check(PropSimulation *propSim, real &dt) { for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) { if (propSim->eventMngr.continuousEvents[i].isHappening) { // if any continuous event is happening, set dt to 1/10th of the event duration - dt = propSim->eventMngr.continuousEvents[i].dt/10.0L; + // dt is the min of the current dt and the one dictated by the event + dt = fmin(1.0, propSim->eventMngr.continuousEvents[i].dt/10.0L); dt = forwardProp ? dt : -dt; } if (forwardProp && propSim->t < propSim->eventMngr.continuousEvents[i].t) { diff --git a/src/interpolate.cpp b/src/interpolate.cpp index bae8ee7d..88fd9dc5 100644 --- a/src/interpolate.cpp +++ b/src/interpolate.cpp @@ -216,38 +216,30 @@ void get_interpIdxInWindow(const PropSimulation *propSim, bool &interpIdxInWindow) { interpIdxInWindow = false; const size_t interpIdx = propSim->interpIdx; + real tCheck = propSim->tEval[interpIdx]; + if (propSim->tEvalUTC) { + tCheck += delta_et_utc(tCheck)/86400.0; + } bool fwInWindow, fwInWindowMarginStart, fwInWindowMarginEnd; bool bwInWindow, bwInWindowMarginStart, bwInWindowMarginEnd; - fwInWindow = - (forwardProp && propSim->tEval[interpIdx] >= tWindowStart && - propSim->tEval[interpIdx] < tNext); + fwInWindow = (forwardProp && tCheck >= tWindowStart && tCheck < tNext); fwInWindowMarginStart = - (forwardProp && - propSim->tEval[interpIdx] <= propSim->integParams.t0 && - propSim->tEval[interpIdx] + propSim->tEvalMargin >= - propSim->integParams.t0 && - propSim->tEval[interpIdx] + propSim->tEvalMargin >= tWindowStart); + (forwardProp && tCheck <= propSim->integParams.t0 && + tCheck + propSim->tEvalMargin >= propSim->integParams.t0 && + tCheck + propSim->tEvalMargin >= tWindowStart); fwInWindowMarginEnd = - (forwardProp && - propSim->tEval[interpIdx] >= propSim->integParams.tf && - propSim->tEval[interpIdx] - propSim->tEvalMargin <= - propSim->integParams.tf && - propSim->tEval[interpIdx] - propSim->tEvalMargin <= tNext); - bwInWindow = - (backwardProp && propSim->tEval[interpIdx] <= tWindowStart && - propSim->tEval[interpIdx] > tNext); + (forwardProp && tCheck >= propSim->integParams.tf && + tCheck - propSim->tEvalMargin <= propSim->integParams.tf && + tCheck - propSim->tEvalMargin <= tNext); + bwInWindow = (backwardProp && tCheck <= tWindowStart && tCheck > tNext); bwInWindowMarginStart = - (backwardProp && - propSim->tEval[interpIdx] >= propSim->integParams.t0 && - propSim->tEval[interpIdx] - propSim->tEvalMargin <= - propSim->integParams.t0 && - propSim->tEval[interpIdx] - propSim->tEvalMargin <= tWindowStart); + (backwardProp && tCheck >= propSim->integParams.t0 && + tCheck - propSim->tEvalMargin <= propSim->integParams.t0 && + tCheck - propSim->tEvalMargin <= tWindowStart); bwInWindowMarginEnd = - (backwardProp && - propSim->tEval[interpIdx] <= propSim->integParams.tf && - propSim->tEval[interpIdx] + propSim->tEvalMargin >= - propSim->integParams.tf && - propSim->tEval[interpIdx] + propSim->tEvalMargin >= tNext); + (backwardProp && tCheck <= propSim->integParams.tf && + tCheck + propSim->tEvalMargin >= propSim->integParams.tf && + tCheck + propSim->tEvalMargin >= tNext); interpIdxInWindow = fwInWindow || fwInWindowMarginStart || fwInWindowMarginEnd || bwInWindow || bwInWindowMarginStart || bwInWindowMarginEnd; From e5edc8bf9b2d21dab65c90456ae609498e02eb8f Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Wed, 28 Aug 2024 15:02:27 -0500 Subject: [PATCH 04/19] added chi value logging and postfit summary --- grss/fit/fit_ades.py | 1 + grss/fit/fit_simulation.py | 79 +++++++++++++++++++++++--------------- grss/version.txt | 2 +- 3 files changed, 51 insertions(+), 31 deletions(-) diff --git a/grss/fit/fit_ades.py b/grss/fit/fit_ades.py index cc7e04f9..fcdbd04b 100644 --- a/grss/fit/fit_ades.py +++ b/grss/fit/fit_ades.py @@ -32,6 +32,7 @@ } ades_grss_columns = { 'obsTimeMJD': 'float', 'obsTimeMJDTDB': 'float', 'cosDec': 'float', + 'resChi': 'float', } ades_column_types = ades_keep_columns | ades_add_columns | ades_grss_columns diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index b9e1b328..e36da36d 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -46,7 +46,7 @@ def __init__(self, iter_number, x_nom, covariance, obs, rms_u, rms_w, chi_sq): self.covariance = covariance variance = np.sqrt(np.diag(self.covariance)) self.variance = dict(zip([f'var_{k}' for k in self.x_nom.keys()], variance)) - self.obs = obs + self.obs = obs.copy() force_del_idx = self.obs[self.obs['selAst'] == 'd'].index auto_del_idx = self.obs[self.obs['selAst'] == 'D'].index self.rejected_idx = np.concatenate((force_del_idx, auto_del_idx)) @@ -549,7 +549,7 @@ def __init__(self, x_init, obs_df, cov_init=None, n_iter_max=10, self._compute_obs_weights() self.n_iter = 0 self.n_iter_max = n_iter_max - self.iters = [] + self.iters = [[]] self.de_kernel = de_kernel self.de_kernel_path = default_kernel_path self.analytic_partials = False @@ -1546,6 +1546,7 @@ def _get_rms_and_reject_outliers(self, partials, residuals, start_rejecting): chi_sq = 0 j = 0 sel_ast = self.obs['selAst'].values + res_chisq_vals = np.nan*np.ones(len(self.observer_info_lengths)) for i, obs_info_len in enumerate(self.observer_info_lengths): if obs_info_len in {4, 7}: size = 2 @@ -1554,8 +1555,8 @@ def _get_rms_and_reject_outliers(self, partials, residuals, start_rejecting): else: raise ValueError("Observer info length not recognized.") resid = residuals[i] - # calculate chi-squared for each residual - if start_rejecting: + # calculate chi-squared for each residual if after the first iteration + if self.n_iter > 1: obs_cov = self.obs_cov[i] obs_partials = partials[j:j+size, :] if sel_ast[i] in {'D', 'd'}: @@ -1567,22 +1568,24 @@ def _get_rms_and_reject_outliers(self, partials, residuals, start_rejecting): else: resid_cov_det = resid_cov[0,0]*resid_cov[1,1] - resid_cov[0,1]*resid_cov[1,0] resid_cov_inv = np.array([[resid_cov[1,1], -resid_cov[0,1]], - [-resid_cov[1,0], resid_cov[0,0]]])/resid_cov_det + [-resid_cov[1,0], resid_cov[0,0]]])/resid_cov_det residual_chi_squared = resid @ resid_cov_inv @ resid.T - # outlier rejection, only reject RA/Dec measurements - if size == 2: - if residual_chi_squared > chi_reject**2 and sel_ast[i] not in {'a', 'd'}: - if sel_ast[i] == 'A': - self.num_rejected += 1 - sel_ast[i] = 'D' - elif residual_chi_squared < chi_recover**2 and sel_ast[i] == 'D': - sel_ast[i] = 'A' - self.num_rejected -= 1 + res_chisq_vals[i] = residual_chi_squared + # outlier rejection, only reject RA/Dec measurements + if start_rejecting and size == 2: + if residual_chi_squared > chi_reject**2 and sel_ast[i] not in {'a', 'd'}: + if sel_ast[i] == 'A': + self.num_rejected += 1 + sel_ast[i] = 'D' + elif residual_chi_squared < chi_recover**2 and sel_ast[i] == 'D': + sel_ast[i] = 'A' + self.num_rejected -= 1 if sel_ast[i] not in {'D', 'd'}: rms_u += resid @ resid.T chi_sq += resid @ self.obs_weight[i] @ resid.T j += size self.obs['selAst'] = sel_ast + self.obs['resChi'] = res_chisq_vals**0.5 rejected_fraction = self.num_rejected/len(self.obs) if rejected_fraction > 0.25: print("WARNING: More than 25% of observations rejected. Consider changing the", @@ -1721,9 +1724,6 @@ def filter_lsq(self, verbose=True): # get new covariance self.covariance = cov # get new state - if i == 0: - # add prefit iteration - self._add_iteration(0, rms_u, rms_w, chi_sq) if self.constraint_dir is not None: constr_hat = self.constraint_dir/np.linalg.norm(self.constraint_dir) delta_x -= np.dot(delta_x, constr_hat)*constr_hat @@ -1762,11 +1762,17 @@ def filter_lsq(self, verbose=True): print(msg) if not self.reject_outliers: break + # add postfit iteration if converged + if self.converged: + residuals, partials = self._get_residuals_and_partials() + rms_u, rms_w, chi_sq = self._get_rms_and_reject_outliers(partials, residuals, + start_rejecting) + self._add_iteration(self.n_iter+1, rms_u, rms_w, chi_sq) if self.n_iter == self.n_iter_max and not self.converged: print("WARNING: Maximum number of iterations reached without converging.") return None - def print_summary(self, iter_idx=-1, out_stream=sys.stdout): + def print_summary(self, iter_idx=None, out_stream=sys.stdout): """ Prints a summary of the orbit fit calculations at a given iteration. @@ -1780,10 +1786,14 @@ def print_summary(self, iter_idx=-1, out_stream=sys.stdout): None : NoneType None """ + iter_idx = len(self.iters)-1 if iter_idx is None else iter_idx data = self.iters[iter_idx] arc = self.obs.obsTimeMJD.max() - self.obs.obsTimeMJD.min() - str_to_print = "Summary of the orbit fit calculations at iteration " - str_to_print += f"{data.iter_number} (of {self.n_iter}):\n" + str_to_print = "Summary of the orbit fit calculations" + if data.iter_number <= self.n_iter: + str_to_print += f" after iteration {data.iter_number} (of {self.n_iter}):\n" + else: + str_to_print += " after postfit pass:\n" str_to_print += "==============================================================\n" str_to_print += f"RMS unweighted: {data.unweighted_rms}\n" str_to_print += f"RMS weighted: {data.weighted_rms}\n" @@ -1799,7 +1809,7 @@ def print_summary(self, iter_idx=-1, out_stream=sys.stdout): str_to_print += "\t\t\tUncertainty\t\t\tChange\t\t\t\tChange (sigma)\n" init_variance = np.sqrt(np.diag(self.covariance_init)) final_variance = np.sqrt(np.diag(self.covariance)) - init_sol = self.iters[0].x_nom + init_sol = self.x_init final_sol = data.x_nom with np.errstate(divide='ignore'): for i, key in enumerate(init_sol.keys()): @@ -1835,7 +1845,9 @@ def plot_summary(self, auto_close=False): """ plt.rcParams['font.size'] = 12 plt.rcParams['axes.labelsize'] = 12 - ticks = np.arange(1, self.n_iter+1, 1) + ticks = np.arange(1, self.n_iter+2, 1) + labels = list(ticks) + labels[-1] = "Postfit" start_idx = 1 iters_for_plot = self.iters[start_idx:] final_iter = iters_for_plot[-1] @@ -1845,7 +1857,8 @@ def plot_summary(self, auto_close=False): [iteration.unweighted_rms for iteration in iters_for_plot], label=f"Final Unweighted RMS={final_iter.unweighted_rms:.3e}") plt.xticks(ticks) - plt.xlabel("Iteration #") + plt.gca().set_xticklabels(labels) + plt.xlabel("Iteration") plt.ylabel("Unweighted RMS") plt.legend() plt.subplot(2, 2, 2) @@ -1853,7 +1866,8 @@ def plot_summary(self, auto_close=False): [iteration.weighted_rms for iteration in iters_for_plot], label=f"Final Weighted RMS={final_iter.weighted_rms:.3e}") plt.xticks(ticks) - plt.xlabel("Iteration #") + plt.gca().set_xticklabels(labels) + plt.xlabel("Iteration") plt.ylabel("Weighted RMS") plt.legend() plt.subplot(2, 2, 3) @@ -1861,7 +1875,8 @@ def plot_summary(self, auto_close=False): [iteration.chi_squared for iteration in iters_for_plot], label=fr"Final $\chi^2$={final_iter.chi_squared:.3e}") plt.xticks(ticks) - plt.xlabel("Iteration #") + plt.gca().set_xticklabels(labels) + plt.xlabel("Iteration") plt.ylabel(r"$\chi^2$") plt.legend() plt.subplot(2, 2, 4) @@ -1869,7 +1884,8 @@ def plot_summary(self, auto_close=False): [iteration.reduced_chi_squared for iteration in iters_for_plot], label=fr"Final Reduced $\chi^2$={final_iter.reduced_chi_squared:.3e}") plt.xticks(ticks) - plt.xlabel("Iteration #") + plt.gca().set_xticklabels(labels) + plt.xlabel("Iteration") plt.ylabel(r"Reduced $\chi^2$") plt.legend() plt.tight_layout() @@ -2033,8 +2049,11 @@ def save(self, filename): f.write(f'{iter_summary_buff}{iter_summary_str}{iter_summary_buff}\n') f.write(subsection_full + '\n') for itrn in self.iters[1:]: - f.write("Summary of the orbit fit calculations at iteration") - f.write(f" {itrn.iter_number} (of {self.n_iter}):\n") + f.write("Summary of the orbit fit calculations") + if itrn.iter_number <= self.n_iter: + f.write(f" after iteration {itrn.iter_number} (of {self.n_iter}):\n") + else: + f.write(" after postfit pass:\n") f.write(f"RMS unweighted: {itrn.unweighted_rms}\n") f.write(f"RMS weighted: {itrn.weighted_rms}\n") f.write(f"chi-squared: {itrn.chi_squared}\n") @@ -2043,7 +2062,7 @@ def save(self, filename): f.write(f'{"Variable":<8}{"Initial Value":>20}') f.write(f'{"Uncertainty":>20}{"Fitted Value":>20}') f.write(f'{"Uncertainty":>20}{"Change":>20}{"Change":>8}'+' (\u03C3)\n') - init_sol = self.iters[0].x_nom + init_sol = self.x_init final_sol = itrn.x_nom with np.errstate(divide='ignore'): for i, key in enumerate(init_sol.keys()): @@ -2064,7 +2083,7 @@ def save(self, filename): f.write(f'{final_val-init_val:+18.11e}{half_tab}') f.write(f'{(final_val-init_val)/init_unc:+10.3f}\n') f.write("\n") - if itrn.iter_number != self.n_iter: + if itrn.iter_number <= self.n_iter: f.write(subsection_full + '\n') mean_0 = np.array(list(self.x_init.values())) diff --git a/grss/version.txt b/grss/version.txt index d87edbfc..078bf8b7 100644 --- a/grss/version.txt +++ b/grss/version.txt @@ -1 +1 @@ -4.2.1 \ No newline at end of file +4.2.2 \ No newline at end of file From dce336c96023975e17ba1e1cb9d1d92a2fda20b2 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Wed, 28 Aug 2024 16:36:43 -0500 Subject: [PATCH 05/19] new PCK files just dropped --- grss/fit/fit_radar.py | 9 --------- grss/kernels/get_kernels.py | 20 ++++++++++---------- grss/version.txt | 2 +- src/simulation.cpp | 13 +++++++------ tests/cpp/prop/pck_map.cpp | 10 +++++----- 5 files changed, 23 insertions(+), 31 deletions(-) diff --git a/grss/fit/fit_radar.py b/grss/fit/fit_radar.py index 935a4e4c..b2806820 100644 --- a/grss/fit/fit_radar.py +++ b/grss/fit/fit_radar.py @@ -71,7 +71,6 @@ def add_radar_obs(obs_df, t_min_tdb=None, t_max_tdb=None, verbose=False): print(f"Read in {num_obs} radar observations from JPL radar API.") data = raw_data['data'] time_range_count = 0 - pre1972_count = 0 for i in range(num_obs): obs = data[num_obs-i-1] date = Time(obs[1], format='iso', scale='utc') @@ -105,14 +104,6 @@ def add_radar_obs(obs_df, t_min_tdb=None, t_max_tdb=None, verbose=False): obs_df.loc[idx,'frq'] = freq obs_df.loc[idx,'sigDelay'] = obs_sigma if delay else np.nan obs_df.loc[idx,'sigDoppler'] = obs_sigma if doppler else np.nan - if date.utc.mjd < 41317: - pre1972_count += 1 - obs_df.loc[idx,'selAst'] = 'd' - if pre1972_count > 0: - print(f"\tWARNING: {pre1972_count} radar observations were taken before 1972.\n" - "\t\tThese residuals cannot be reliably computed because high-accuracy Earth " - "orientation kernels are not available for this time period.\n" - "\t\tForce deleting them...") if verbose: print(f"\tFiltered to {num_obs-time_range_count} observations that satisfy the " "time range constraints.") diff --git a/grss/kernels/get_kernels.py b/grss/kernels/get_kernels.py index 817af211..c9c03490 100644 --- a/grss/kernels/get_kernels.py +++ b/grss/kernels/get_kernels.py @@ -37,26 +37,26 @@ # get the earth orientation binary spice kernels and their comments if they are not already present # latest earth pck -latest_earth_pck_exists = os.path.exists(f'{script_dir}/earth_latest_high_prec.bpc') +latest_earth_pck_exists = os.path.exists(f'{script_dir}/earth_latest.bpc') latest_earth_pck_old = (latest_earth_pck_exists and - time.time() - os.path.getmtime(f'{script_dir}/earth_latest_high_prec.bpc') + time.time() - os.path.getmtime(f'{script_dir}/earth_latest.bpc') > 86400) if not latest_earth_pck_exists or latest_earth_pck_old: print('Downloading latest Earth binary PCK...') - os.system((f'curl --silent --show-error -o {script_dir}/earth_latest_high_prec.bpc ' + os.system((f'curl --silent --show-error -o {script_dir}/earth_latest.bpc ' f'{NAIF_SITE}/pck/earth_latest_high_prec.bpc')) # historical earth pck -historical_earth_pck_exists = os.path.exists(f'{script_dir}/earth_720101_230601.bpc') +historical_earth_pck_exists = os.path.exists(f'{script_dir}/earth_historic.bpc') if not historical_earth_pck_exists: - print('Downloading historical Earth binary PCK...') - os.system((f'curl --silent --show-error -o {script_dir}/earth_720101_230601.bpc ' - f'{NAIF_SITE}/pck/earth_720101_230601.bpc')) + print('Downloading historic Earth binary PCK...') + os.system((f'curl --silent --show-error -o {script_dir}/earth_historic.bpc ' + f'{NAIF_SITE}/pck/earth_620120_240827.bpc')) # predicted earth pck -predicted_earth_pck_exists = os.path.exists(f'{script_dir}/earth_200101_990825_predict.bpc') +predicted_earth_pck_exists = os.path.exists(f'{script_dir}/earth_predict.bpc') if not predicted_earth_pck_exists: print('Downloading predicted Earth binary PCK...') - os.system((f'curl --silent --show-error -o {script_dir}/earth_200101_990825_predict.bpc ' - f'{NAIF_SITE}/pck/earth_200101_990825_predict.bpc')) + os.system((f'curl --silent --show-error -o {script_dir}/earth_predict.bpc ' + f'{NAIF_SITE}/pck/earth_200101_990827_predict.bpc')) # generic frame kernels generic_frame_kernel_exists = os.path.exists(f'{script_dir}/pck00011.tpc') if not generic_frame_kernel_exists: diff --git a/grss/version.txt b/grss/version.txt index 078bf8b7..ec87108d 100644 --- a/grss/version.txt +++ b/grss/version.txt @@ -1 +1 @@ -4.2.2 \ No newline at end of file +4.2.3 \ No newline at end of file diff --git a/src/simulation.cpp b/src/simulation.cpp index 7aa038d6..39dcd48c 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -21,9 +21,10 @@ void get_baseBodyFrame(const int &spiceId, const real &tMjdTDB, break; case 399: baseBodyFrame = "ITRF93"; - // High precision frame is not defined before 1972 JAN 01 - // 00:00:42.183 TDB or after 2099 AUG 25 00:01:09.182 TDB - if (tMjdTDB < 41317.000488239L || tMjdTDB > 87940.000800717L) { + // High precision frame is not defined + // before 1962 JAN 20 00:00:41.184 TDB or + // after 2099 AUG 27 00:01:09.182 TDB + if (tMjdTDB < 37684.0004767L || tMjdTDB > 87942.0008007L) { baseBodyFrame = "IAU_EARTH"; } break; @@ -691,9 +692,9 @@ PropSimulation::PropSimulation(std::string name, real t0, } this->spkEphem.mbPath = kernel_mb; this->spkEphem.sbPath = kernel_sb; - this->pckEphem.histPckPath = DEkernelPath + "earth_720101_230601.bpc"; - this->pckEphem.latestPckPath = DEkernelPath + "earth_latest_high_prec.bpc"; - this->pckEphem.predictPckPath = DEkernelPath + "earth_200101_990825_predict.bpc"; + this->pckEphem.histPckPath = DEkernelPath + "earth_historic.bpc"; + this->pckEphem.latestPckPath = DEkernelPath + "earth_latest.bpc"; + this->pckEphem.predictPckPath = DEkernelPath + "earth_predict.bpc"; } /** diff --git a/tests/cpp/prop/pck_map.cpp b/tests/cpp/prop/pck_map.cpp index f1facd3a..301c6d53 100644 --- a/tests/cpp/prop/pck_map.cpp +++ b/tests/cpp/prop/pck_map.cpp @@ -12,8 +12,8 @@ int main(){ std::random_device rd; // Will be used to obtain a seed for the random number engine std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd() // limits of distribution are dictated by the time span of - // earth_720101_230601.bpc and earth_200101_990825_predict.bpc - std::uniform_real_distribution<> dis(41317.0005, 87940.0008); + // earth_620120_240827.bpc and earth_200101_990827_predict.bpc + std::uniform_real_distribution<> dis(37684.0004767L, 87942.0008007L); std::cout << "/////////////////////// PCK map accuracy test ///////////////////////" @@ -24,9 +24,9 @@ int main(){ std::string from = "ITRF93"; std::string to = "J2000"; std::string pck_text = "../../../grss/kernels/pck00011.tpc"; - std::string pck_hist = "../../../grss/kernels/earth_720101_230601.bpc"; - std::string pck_latest = "../../../grss/kernels/earth_latest_high_prec.bpc"; - std::string pck_predict = "../../../grss/kernels/earth_200101_990825_predict.bpc"; + std::string pck_hist = "../../../grss/kernels/earth_historic.bpc"; + std::string pck_latest = "../../../grss/kernels/earth_latest.bpc"; + std::string pck_predict = "../../../grss/kernels/earth_predict.bpc"; //////////// GRSS //////////// PckEphemeris pckEphem; From b6743470a148424b2c717e37ebe88361f508716b Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Sun, 1 Sep 2024 19:00:16 -0500 Subject: [PATCH 06/19] updated light deflection --- src/observe.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/observe.cpp b/src/observe.cpp index 730d4990..41ab3185 100644 --- a/src/observe.cpp +++ b/src/observe.cpp @@ -39,10 +39,7 @@ void get_glb_correction(PropSimulation *propSim, const size_t &interpIdx, std::vector p(3, 0.0); vunit(observerTargetPos, p); - std::vector deltaP1Targ(3, 0.0); - std::vector deltaP1Star(3, 0.0); std::vector p1(p); - size_t numBendingBodies = 1; // sun (1) or sun+earth (2) for (size_t i = 0; i < numBendingBodies; i++){ double centralBodyState[6]; @@ -89,13 +86,15 @@ void get_glb_correction(PropSimulation *propSim, const size_t &interpIdx, vdot(e, p, eDotP); vdot(q, e, qDotE); - real g1 = 2 * centralBodyGM / c / c / centralBodyObserverDist; - real g2 = 1.0L + qDotE; + const real g1 = 2 * centralBodyGM / c / c / centralBodyObserverDist; + const real g2 = 1.0L + qDotE; + std::vector deltaP1Targ(3, 0.0); deltaP1Targ[0] += g1 * (pDotQ * e[0] - eDotP * q[0]) / g2; deltaP1Targ[1] += g1 * (pDotQ * e[1] - eDotP * q[1]) / g2; deltaP1Targ[2] += g1 * (pDotQ * e[2] - eDotP * q[2]) / g2; + std::vector deltaP1Star(3, 0.0); deltaP1Star[0] += g1 * (e[0] - eDotP * p[0]) / (1 + eDotP); deltaP1Star[1] += g1 * (e[1] - eDotP * p[1]) / (1 + eDotP); deltaP1Star[2] += g1 * (e[2] - eDotP * p[2]) / (1 + eDotP); @@ -110,6 +109,11 @@ void get_glb_correction(PropSimulation *propSim, const size_t &interpIdx, p1[1] -= deltaP1Star[1]; p1[2] -= deltaP1Star[2]; } + // re-normalize + const real p1Norm = sqrt(p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2]); + for (size_t j = 0; j < 3; j++) { + p1[j] /= p1Norm; + } } observerTargetPos[0] = observerTargetDist * p1[0]; From 49db4b23c3f8e5a72f3654fe5dca0c2c9bfda4ee Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Wed, 11 Sep 2024 10:43:01 -0500 Subject: [PATCH 07/19] added a prioris --- grss/fit/fit_simulation.py | 127 +++++++++++++++++++++++++++++++------ 1 file changed, 106 insertions(+), 21 deletions(-) diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index e36da36d..347f2182 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -561,6 +561,12 @@ def __init__(self, x_init, obs_df, cov_init=None, n_iter_max=10, for key in nongrav_info: self.fixed_propsim_params[key] = nongrav_info[key] self.fixed_propsim_params['events'] = events if events is not None else [] + self.prior_est = None + self.prior_sigs = None + self._priors_given = False + self._xbar0 = None + self._info0 = None + self._prior_constant = None self.reject_outliers = True self.reject_criteria = [3.0, 2.8] # number of rejected obs is the count of ['D', 'd'] in selAst @@ -621,6 +627,38 @@ def _check_initial_solution(self, x_init, cov_init): raise ValueError(msg) return None + def _check_priors(self): + """ + Check the prior estimates and sigmas provided by the user + and make sure they are valid. + + Returns + ------- + None : NoneType + None + + Raises + ------ + ValueError + If the prior estimates and sigmas do not have the same length. + """ + # check that the keys in self.x_init are the same as in self.prior_est and self.prior_sigs + if self.prior_est is not None and self.prior_sigs is not None: + self._priors_given = True + if len(self.prior_est) != len(self.prior_sigs): + raise ValueError("Prior estimates and sigmas must have the same length.") + for key in self.prior_est.keys(): + if key not in self.x_init: + raise ValueError(f"Key {key} not found in initial state vector.") + self._info0 = np.zeros_like(self.covariance_init) + self._xbar0 = np.zeros(self.n_fit) + for key in self.prior_est.keys(): + idx = list(self.x_init.keys()).index(key) + self._info0[idx, idx] = 1/self.prior_sigs[key]**2 + self._xbar0[idx] = self.prior_est[key] - self.x_init[key] + self._prior_constant = list(self.x_nom.values())+self._xbar0 + return None + def _flatten_and_clean(self, arr): """ Flatten an array and remove any NaN values. @@ -1584,6 +1622,11 @@ def _get_rms_and_reject_outliers(self, partials, residuals, start_rejecting): rms_u += resid @ resid.T chi_sq += resid @ self.obs_weight[i] @ resid.T j += size + # # write res_chisq_vals to file if any values are negative + # if np.any(res_chisq_vals < 0): + # print("WARNING: Negative outlier rejection chi-squared values detected. " + # "Writing to file 'warn_neg_chisq_vals.txt'.") + # np.savetxt("warn_neg_chisq_vals.txt", res_chisq_vals, fmt="%.16f") self.obs['selAst'] = sel_ast self.obs['resChi'] = res_chisq_vals**0.5 rejected_fraction = self.num_rejected/len(self.obs) @@ -1667,8 +1710,12 @@ def _get_lsq_state_correction(self, partials, residuals): delta_x : array The state correction. """ - atwa = np.zeros((self.n_fit, self.n_fit)) - atwb = np.zeros(self.n_fit) + if self._priors_given: + atwa = self._info0.copy() + atwb = (self._info0 @ self._xbar0).copy() + else: + atwa = np.zeros((self.n_fit, self.n_fit)) + atwb = np.zeros(self.n_fit) j = 0 sel_ast = self.obs['selAst'].values for i, obs_info_len in enumerate(self.observer_info_lengths): @@ -1686,11 +1733,11 @@ def _get_lsq_state_correction(self, partials, residuals): j += size # use pseudo-inverse if the data arc is less than 7 days if self.obs.obsTimeMJD.max() - self.obs.obsTimeMJD.min() < 7.0: - cov = np.linalg.pinv(atwa, rcond=1e-20, hermitian=True) + self.covariance = np.linalg.pinv(atwa, rcond=1e-20, hermitian=True) else: - cov = np.array(libgrss.matrix_inverse(atwa)) - delta_x = cov @ atwb - return delta_x.ravel(), cov + self.covariance = np.array(libgrss.matrix_inverse(atwa)) + delta_x = self.covariance @ atwb + return delta_x.ravel() def filter_lsq(self, verbose=True): """ @@ -1706,38 +1753,52 @@ def filter_lsq(self, verbose=True): None : NoneType None """ + self._check_priors() start_rejecting = False if verbose: print("Iteration\t\tUnweighted RMS\t\tWeighted RMS", "\t\tChi-squared\t\tReduced Chi-squared") + delta_x = np.zeros(self.n_fit) for i in range(self.n_iter_max): + if self._priors_given: + self._xbar0 -= delta_x + a_priori_constant = list(self.x_nom.values())+self._xbar0 + a_priori_constant_violation = np.linalg.norm(a_priori_constant-self._prior_constant) + if a_priori_constant_violation >= 1e-11: + print(f"a priori constraint constant violation before iteration {i+1}: " + f"{a_priori_constant_violation}. Solution may not be trustworthy. " + "Violation for each estimated parameter:") + print(dict(zip(self.x_nom.keys(), a_priori_constant-self._prior_constant))) self.n_iter = i+1 # get residuals and partials residuals, partials = self._get_residuals_and_partials() # calculate rms and reject outliers here if desired rms_u, rms_w, chi_sq = self._get_rms_and_reject_outliers(partials, residuals, start_rejecting) - # get initial guess - curr_state = np.array(list(self.x_nom.values())) + # get current state + curr_state = list(self.x_nom.values()) # get state correction - delta_x, cov = self._get_lsq_state_correction(partials, residuals) - # get new covariance - self.covariance = cov - # get new state + delta_x = self._get_lsq_state_correction(partials, residuals) if self.constraint_dir is not None: constr_hat = self.constraint_dir/np.linalg.norm(self.constraint_dir) delta_x -= np.dot(delta_x, constr_hat)*constr_hat - next_state = curr_state + delta_x - if self.fit_cometary and next_state[0] < 0.0: - next_state[0] = 0.0 - print("WARNING: Eccentricity is negative per least squares state correction. " - "Setting to 0.0. This solution may not be trustworthy.") - self.x_nom = dict(zip(self.x_nom.keys(), next_state)) - # make sure any keys starting with dt are positive and stay in the range [0, 180] + # make sure eccentricity is non-negative and + # any keys starting with dt are 0.1 at minimum for key in self.x_nom: + if self.fit_cometary and key == 'e': + idx = list(self.x_nom.keys()).index(key) + if curr_state[idx] + delta_x[idx] < 0.0: + delta_x[idx] = -curr_state[idx] + print("WARNING: Eccentricity is negative per least squares " + "state correction. Setting to 0.0. This solution " + "may not be trustworthy.") if key[:2] == 'dt': - self.x_nom[key] = abs(self.x_nom[key]) - self.x_nom[key] = min(self.x_nom[key], 180.0) + idx = list(self.x_nom.keys()).index(key) + if curr_state[idx] + delta_x[idx] < 0.1: + delta_x[idx] = 0.1 - curr_state[idx] + # get new nominal state + next_state = curr_state + delta_x + self.x_nom = dict(zip(self.x_nom.keys(), next_state)) # add iteration self._add_iteration(i+1, rms_u, rms_w, chi_sq) if verbose: @@ -1968,6 +2029,18 @@ def save(self, filename): f.write(f"{key.upper()}: {val:.11f} \u00B1 {sig:.11f}") if key in units_dict: f.write(f" {units_dict[key]}") + if self._priors_given and key in self.prior_est: + p_val = self.prior_est[key] + p_sig = self.prior_sigs[key] + if key in ['om', 'w', 'i']: + p_val *= 180/np.pi + p_sig *= 180/np.pi + if key in ['a1', 'a2', 'a3']: + f.write(f" (a priori: {p_val:.11e} \u00B1 {p_sig:.11e})") + else: + f.write(f" (a priori: {p_val:.11f} \u00B1 {p_sig:.11f})") + if key in units_dict: + f.write(f" {units_dict[key]}") f.write("\n") f.write("\n") f.write("Initial Covariance Matrix:\n") @@ -1995,6 +2068,18 @@ def save(self, filename): f.write(f"{key.upper()}: {val:.11f} \u00B1 {sig:.11f}") if key in units_dict: f.write(f" {units_dict[key]}") + if self._priors_given and key in self.prior_est: + p_val = self.prior_est[key] + p_sig = self.prior_sigs[key] + if key in ['om', 'w', 'i']: + p_val *= 180/np.pi + p_sig *= 180/np.pi + if key in ['a1', 'a2', 'a3']: + f.write(f" (a priori: {p_val:.11e} \u00B1 {p_sig:.11e})") + else: + f.write(f" (a priori: {p_val:.11f} \u00B1 {p_sig:.11f})") + if key in units_dict: + f.write(f" {units_dict[key]}") f.write("\n") f.write("\n") f.write("Final Covariance Matrix:\n") From bffa7ac7ab128ceae87d92ac962fbc6c9e0ed40d Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Sat, 14 Sep 2024 19:39:46 -0500 Subject: [PATCH 08/19] OD cleanup --- grss/fit/fit_optical.py | 3 ++- grss/fit/fit_simulation.py | 14 +++++++++----- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/grss/fit/fit_optical.py b/grss/fit/fit_optical.py index a396b8bb..06789485 100644 --- a/grss/fit/fit_optical.py +++ b/grss/fit/fit_optical.py @@ -479,7 +479,8 @@ def apply_debiasing_scheme(obs_df, lowres, verbose): unbias_counter += len(group) continue if cat_code not in biased_catalogs: - print(f"\tUnknown star catalog: {cat}") + if verbose: + print(f"\tUnknown star catalog: {cat}") obs_df.loc[group.index, 'biasRA'] = 0.0 obs_df.loc[group.index, 'biasDec'] = 0.0 continue diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index 347f2182..710b4bf9 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -691,14 +691,18 @@ def _add_simulated_obs(self): weights = self.simulated_obs['weights'] if not len(times) == len(modes) == len(weights): raise ValueError("Simulated observation data must have the same length.") + perm_id = self.obs['permID'][0] + prov_id = self.obs['provID'][0] for i, time in enumerate(times): idx = len(self.obs) mode = modes[i] if mode not in {'SIM_CCD', 'SIM_OCC', 'SIM_RAD_DEL', 'SIM_RAD_DOP'}: raise ValueError(f"Unknown simulated observation mode {mode}.") weight = weights[i] - self.obs.loc[idx, 'permID'] = 'SIM_'+str(self.obs['permID'][0]) - self.obs.loc[idx, 'provID'] = 'SIM_'+str(self.obs['provID'][0]) + if isinstance(perm_id, str): + self.obs.loc[idx, 'permID'] = f'SIM_{perm_id}' + if isinstance(prov_id, str): + self.obs.loc[idx, 'provID'] = f'SIM_{prov_id}' self.obs.loc[idx, 'obsTime'] = f'{time.utc.isot}Z' self.obs.loc[idx, 'obsTimeMJD'] = time.utc.mjd self.obs.loc[idx, 'obsTimeMJDTDB'] = time.tdb.mjd @@ -1347,7 +1351,7 @@ def _inflate_uncertainties(self, prop_sim_past, prop_sim_future): | special_codes["spacecraft"] ) for i in range(len(self.observer_info_lengths)): - if modes[i] not in {'RAD', 'SIM_RAD_DEL', 'SIM_RAD_DOP', 'SIM_OCC'}: + if modes[i] not in {'RAD', 'SIM_RAD_DEL', 'SIM_RAD_DOP', 'SIM_CCD', 'SIM_OCC'}: sig_ra = sig_ra_vals[i] sig_dec = sig_dec_vals[i] sig_corr = sig_corr_vals[i] @@ -1611,11 +1615,11 @@ def _get_rms_and_reject_outliers(self, partials, residuals, start_rejecting): res_chisq_vals[i] = residual_chi_squared # outlier rejection, only reject RA/Dec measurements if start_rejecting and size == 2: - if residual_chi_squared > chi_reject**2 and sel_ast[i] not in {'a', 'd'}: + if abs(residual_chi_squared) > chi_reject**2 and sel_ast[i] not in {'a', 'd'}: if sel_ast[i] == 'A': self.num_rejected += 1 sel_ast[i] = 'D' - elif residual_chi_squared < chi_recover**2 and sel_ast[i] == 'D': + elif abs(residual_chi_squared) < chi_recover**2 and sel_ast[i] == 'D': sel_ast[i] = 'A' self.num_rejected -= 1 if sel_ast[i] not in {'D', 'd'}: From fb6ed6ecf19a9cc5eac66270617115b42a569e3e Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Sat, 14 Sep 2024 20:09:21 -0500 Subject: [PATCH 09/19] change continuous events to exponential; no longer letting them change timestep --- grss/version.txt | 2 +- include/simulation.h | 7 ++++--- src/force.cpp | 24 ++++++++---------------- src/gr15.cpp | 23 ++++++++--------------- src/grss.cpp | 14 +++++++------- src/simulation.cpp | 17 +++++++++-------- 6 files changed, 37 insertions(+), 50 deletions(-) diff --git a/grss/version.txt b/grss/version.txt index ec87108d..81911389 100644 --- a/grss/version.txt +++ b/grss/version.txt @@ -1 +1 @@ -4.2.3 \ No newline at end of file +4.3.0 \ No newline at end of file diff --git a/include/simulation.h b/include/simulation.h index 44e1cd35..b5404ae9 100644 --- a/include/simulation.h +++ b/include/simulation.h @@ -242,7 +242,8 @@ class Event { real multiplier = 1.0L; // for continuous ejecta events - real dt = 1.0L; + std::vector expAccel0 = {0.0L, 0.0L, 0.0L}; + real tau = 1.0L; bool isContinuous = false; bool isHappening = false; /** @@ -522,8 +523,8 @@ class PropSimulation { /** * @brief Add an ejecta event to the simulation. */ - void add_event(IntegBody body, real tEvent, std::vector deltaV, - real multiplier, real dt); + void add_event(IntegBody body, real tEvent, std::vector expAccel0, + real multiplier, real tau); /** * @brief Set the values of the PropSimulation Constants object. */ diff --git a/src/force.cpp b/src/force.cpp index 87ab2556..7ef4df60 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -657,24 +657,16 @@ static void force_continuous_event(const real &t, const PropSimulation *propSim, if (propSim->eventMngr.continuousEvents[i].isHappening){ const size_t starti = propSim->eventMngr.continuousEvents[i].xIntegIndex / 2; - const real tPastEvent = - t - propSim->eventMngr.continuousEvents[i].t; - // f = np.sin(np.pi*x/(2*time_for_val)) - // df = np.pi/(2*time_for_val)*np.cos(np.pi*x/(2*time_for_val)) - const real preFac = - PI / (2 * propSim->eventMngr.continuousEvents[i].dt); - real accFac = preFac * cos(preFac * tPastEvent); - if (accFac > preFac || accFac < 0.0) { - accFac = 0.0; - } - accInteg[starti + 0] += accFac * - propSim->eventMngr.continuousEvents[i].deltaV[0] * + const real tPastEvent = t - propSim->eventMngr.continuousEvents[i].t; + const real postFac = exp(-tPastEvent/propSim->eventMngr.continuousEvents[i].tau); + accInteg[starti + 0] += propSim->eventMngr.continuousEvents[i].expAccel0[0] * + postFac * propSim->eventMngr.continuousEvents[i].multiplier * propDir; - accInteg[starti + 1] += accFac * - propSim->eventMngr.continuousEvents[i].deltaV[1] * + accInteg[starti + 1] += propSim->eventMngr.continuousEvents[i].expAccel0[1] * + postFac * propSim->eventMngr.continuousEvents[i].multiplier * propDir; - accInteg[starti + 2] += accFac * - propSim->eventMngr.continuousEvents[i].deltaV[2] * + accInteg[starti + 2] += propSim->eventMngr.continuousEvents[i].expAccel0[2] * + postFac * propSim->eventMngr.continuousEvents[i].multiplier * propDir; } } diff --git a/src/gr15.cpp b/src/gr15.cpp index 61a6dd89..00f4265b 100644 --- a/src/gr15.cpp +++ b/src/gr15.cpp @@ -255,14 +255,12 @@ void check_continuous_events(PropSimulation *propSim, const real &t) { bool allDone = true; bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) { - // const bool eventStarted = t >= propSim->eventMngr.continuousEvents[i].t; - // const bool eventEnded = t > propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt; bool eventStarted, eventEnded; if (forwardProp) { eventStarted = t >= propSim->eventMngr.continuousEvents[i].t; - eventEnded = t >= (propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt); + eventEnded = t >= propSim->integParams.tf; } else { - eventStarted = t <= (propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt); + eventStarted = t <= propSim->integParams.t0; eventEnded = t < propSim->eventMngr.continuousEvents[i].t; } if (eventStarted && !eventEnded) { @@ -306,17 +304,12 @@ void event_timestep_check(PropSimulation *propSim, real &dt) { } } if (!propSim->eventMngr.allConEventDone) { - for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) { - if (propSim->eventMngr.continuousEvents[i].isHappening) { - // if any continuous event is happening, set dt to 1/10th of the event duration - // dt is the min of the current dt and the one dictated by the event - dt = fmin(1.0, propSim->eventMngr.continuousEvents[i].dt/10.0L); - dt = forwardProp ? dt : -dt; - } - if (forwardProp && propSim->t < propSim->eventMngr.continuousEvents[i].t) { - tNextEvent = fmin(tNextEvent, propSim->eventMngr.continuousEvents[i].t); - } else if (!forwardProp && propSim->t > propSim->eventMngr.continuousEvents[i].t+propSim->eventMngr.continuousEvents[i].dt) { - tNextEvent = fmax(tNextEvent, propSim->eventMngr.continuousEvents[i].t+propSim->eventMngr.continuousEvents[i].dt); + for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) { + const real tNextConEvent = propSim->eventMngr.continuousEvents[i].t; + if (forwardProp && propSim->t < tNextConEvent) { + tNextEvent = fmin(tNextEvent, tNextConEvent); + } else if (!forwardProp && propSim->t > tNextConEvent) { + tNextEvent = fmax(tNextEvent, tNextConEvent); } } } diff --git a/src/grss.cpp b/src/grss.cpp index f2bedbe1..a873b3df 100644 --- a/src/grss.cpp +++ b/src/grss.cpp @@ -884,8 +884,8 @@ PYBIND11_MODULE(libgrss, m) { static_cast, real, real)>( &PropSimulation::add_event), - py::arg("bodyName"), py::arg("tEvent"), py::arg("deltaV"), - py::arg("multiplier"), py::arg("dt"), + py::arg("bodyName"), py::arg("tEvent"), py::arg("expAccel0"), + py::arg("multiplier"), py::arg("tau"), R"mydelimiter( Adds an ejecta event to the simulation. @@ -895,12 +895,12 @@ PYBIND11_MODULE(libgrss, m) { Name of the body to apply the delta-V to. tEvent : real MJD Epoch of the event. Must be in TDB. - deltaV : list of real - Delta-V to apply to the body. + expAccel0 : list of real + Initial exponentiallly decaying acceleration. multiplier : real - Multiplier to apply to the delta-V. - dt : real - Time duration for the continuous ejecta event. + Multiplier for the exponential decay. + tau : real + Time constant for the exponentiallly decaying acceleration. )mydelimiter") .def("set_sim_constants", &PropSimulation::set_sim_constants, py::arg("du2m") = 149597870700.0L, py::arg("tu2s") = 86400.0L, diff --git a/src/simulation.cpp b/src/simulation.cpp index 39dcd48c..2c75429c 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1060,7 +1060,7 @@ void PropSimulation::add_event(IntegBody body, real tEvent, event.deltaV = deltaV; event.multiplier = multiplier; - event.dt = 0.0; + event.tau = 0.0; event.isContinuous = false; event.isHappening = false; event_postprocess(this, event); @@ -1070,22 +1070,23 @@ void PropSimulation::add_event(IntegBody body, real tEvent, /** * @param[in] body IntegBody object to apply the EjectaEvent to. * @param[in] tEvent Time at which to apply the EjectaEvent. - * @param[in] deltaV Delta-V for the ejecta. - * @param[in] multiplier Multiplier for the Delta-V. - * @param[in] dt Time duration for the ejecta. + * @param[in] expAccel0 Initial exponentiallly decaying acceleration. + * @param[in] multiplier Multiplier for the exponential decay. + * @param[in] tau Time constant for the exponentiallly decaying acceleration. */ void PropSimulation::add_event(IntegBody body, real tEvent, - std::vector deltaV, real multiplier, - real dt) { + std::vector expAccel0, real multiplier, + real tau) { size_t bodyIndex = event_preprocess(this, body, tEvent); Event event; event.t = tEvent; event.bodyName = body.name; event.bodyIndex = bodyIndex; - event.deltaV = deltaV; + event.deltaV = std::vector(3, 0.0); event.multiplier = multiplier; - event.dt = dt; + event.expAccel0 = expAccel0; + event.tau = tau; event.isContinuous = true; event.isHappening = false; event_postprocess(this, event); From 69756139fefb0f2f9c18ffd903e7b22da3fc7e08 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Tue, 17 Sep 2024 17:21:34 -0500 Subject: [PATCH 10/19] changed default initial dt; added margin for binary PCK reads --- include/simulation.h | 4 ++-- src/gr15.cpp | 6 ++---- src/grss.cpp | 2 +- src/pck.cpp | 5 +++-- 4 files changed, 8 insertions(+), 9 deletions(-) diff --git a/include/simulation.h b/include/simulation.h index b5404ae9..8ad65c8e 100644 --- a/include/simulation.h +++ b/include/simulation.h @@ -249,7 +249,7 @@ class Event { /** * @brief Apply the impulse event to the body. */ - void apply_impulsive(const real &t, std::vector &xInteg, const real &propDir); + void apply_impulsive(PropSimulation *propSim, const real &t, std::vector &xInteg); }; class EventManager { @@ -543,7 +543,7 @@ class PropSimulation { bool convergedLightTime = false, std::vector> observerInfo = std::vector>(), - bool adaptiveTimestep = true, real dt0 = 0.0L, + bool adaptiveTimestep = true, real dt0 = 1.0L, real dtMin = 1.0e-4L, real dtChangeFactor = 0.25L, real tolInteg = 1.0e-11L, real tolPC = 1.0e-16L); /** diff --git a/src/gr15.cpp b/src/gr15.cpp index 00f4265b..eb19fc83 100644 --- a/src/gr15.cpp +++ b/src/gr15.cpp @@ -5,7 +5,7 @@ * @return real Initial timestep. */ real get_initial_timestep(PropSimulation *propSim){ - real dt0 = propSim->integParams.dtMin; + real dt0 = 1.0L; if (propSim->integParams.dt0 != 0.0) { dt0 = fabs(propSim->integParams.dt0); } @@ -233,9 +233,7 @@ void check_and_apply_impulsive_events(PropSimulation *propSim, const real &t, real *tNextImpEvent = &propSim->eventMngr.tNextImpEvent; while (*nextImpEventIdx < propSim->eventMngr.nImpEvents && t == *tNextImpEvent) { // apply events for the state just reached by the integrator - bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; - real propDir = forwardProp ? 1.0L : -1.0L; - propSim->eventMngr.impulsiveEvents[*nextImpEventIdx].apply_impulsive(t, xInteg, propDir); + propSim->eventMngr.impulsiveEvents[*nextImpEventIdx].apply_impulsive(propSim, t, xInteg); // update next event index and time (*nextImpEventIdx)++; if (*nextImpEventIdx < propSim->eventMngr.nImpEvents) { diff --git a/src/grss.cpp b/src/grss.cpp index a873b3df..2a51cdd1 100644 --- a/src/grss.cpp +++ b/src/grss.cpp @@ -928,7 +928,7 @@ PYBIND11_MODULE(libgrss, m) { py::arg("tEvalUTC") = false, py::arg("evalApparentState") = false, py::arg("convergedLightTime") = false, py::arg("observerInfo") = std::vector>(), - py::arg("adaptiveTimestep") = true, py::arg("dt0") = 0.0L, + py::arg("adaptiveTimestep") = true, py::arg("dt0") = 1.0L, py::arg("dtMin") = 1.0e-4L, py::arg("dtChangeFactor") = 0.25L, py::arg("tolInteg") = 1.0e-11L, py::arg("tolPC") = 1.0e-16L, R"mydelimiter( diff --git a/src/pck.cpp b/src/pck.cpp index 9a458507..f15f2d9a 100644 --- a/src/pck.cpp +++ b/src/pck.cpp @@ -750,10 +750,11 @@ void get_pck_rotMat(const std::string &from, const std::string &to, if (t0_mjd < ephem.histPck->targets[0].beg || t0_mjd > ephem.predictPck->targets[0].end){ throw std::runtime_error("get_pck_rotMat: The epoch is outside the range of the binary PCK files."); } - if (t0_mjd <= ephem.histPck->targets[0].end){ + const real tSwitchMargin = 0.1; + if (t0_mjd <= ephem.histPck->targets[0].end-tSwitchMargin){ bpc = ephem.histPck; // std::cout << "Using histPck" << std::endl; - } else if (t0_mjd <= ephem.latestPck->targets[0].end){ + } else if (t0_mjd <= ephem.latestPck->targets[0].end-tSwitchMargin){ bpc = ephem.latestPck; // std::cout << "Using latestPck" << std::endl; } else { From 0f70402c10d53316e731227d29bfee984d081a28 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Tue, 17 Sep 2024 17:24:22 -0500 Subject: [PATCH 11/19] add delta-v vector estimation --- grss/fit/fit_simulation.py | 15 +++++++++------ src/simulation.cpp | 37 ++++++++++++++++++++++++++++++++----- 2 files changed, 41 insertions(+), 11 deletions(-) diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index 710b4bf9..6aa92dc8 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -843,7 +843,7 @@ def _get_prop_sim_past(self, name, t_eval_utc, eval_apparent_state, """ # pylint: disable=no-member t_eval_past = self.obs.obsTimeMJD.values[self.past_obs_idx] - tf_past = np.min(t_eval_past) + tf_past = np.min(self.obs.obsTimeMJDTDB.values[self.past_obs_idx]) prop_sim_past = libgrss.PropSimulation(name, self.t_sol, self.de_kernel, self.de_kernel_path) prop_sim_past.tEvalMargin = 1.0 @@ -887,7 +887,7 @@ def _get_prop_sim_future(self, name, t_eval_utc, eval_apparent_state, """ # pylint: disable=no-member t_eval_future = self.obs.obsTimeMJD.values[self.future_obs_idx] - tf_future = np.max(t_eval_future) + tf_future = np.max(self.obs.obsTimeMJDTDB.values[self.future_obs_idx]) prop_sim_future = libgrss.PropSimulation(name, self.t_sol, self.de_kernel, self.de_kernel_path) prop_sim_future.tEvalMargin = 1.0 @@ -1340,7 +1340,7 @@ def _inflate_uncertainties(self, prop_sim_past, prop_sim_future): radius_nonzero = self.fixed_propsim_params['radius'] != 0.0 fac = 0.0 if radius_nonzero: - rel_dists = np.linalg.norm(np.array(state_eval)[:, :3], axis=1) + rel_dists = np.linalg.norm([state[:3] for state in state_eval], axis=1) lmbda = 0.3 # from fuentes-munoz et al. 2024 au2m = 1.495978707e11 fac = (lmbda*self.fixed_propsim_params['radius']/au2m/rel_dists*180/np.pi*3600)**2 @@ -1412,7 +1412,8 @@ def _get_analytic_partials(self, prop_sim_past, prop_sim_future): size = 1 else: raise ValueError("Observer info length not recognized.") - part = np.zeros((size, self.n_fit)) + # part is partial of observable with respect to state at observation time + part = np.zeros((size, 6)) if self.past_obs_exist and i < len_past_idx: sim_idx = i optical_partials = past_optical_partials @@ -1423,14 +1424,16 @@ def _get_analytic_partials(self, prop_sim_past, prop_sim_future): optical_partials = future_optical_partials radar_partials = future_radar_partials state = future_states - stm = libgrss.reconstruct_stm(state[sim_idx][6:]) + # stm is partial of state at observation time with respect to nominal state + stm = libgrss.reconstruct_stm(state[sim_idx][6:])[:6] if is_optical: part[0, :6] = optical_partials[sim_idx, :6]*cos_dec[i] part[1, :6] = optical_partials[sim_idx, 6:12] else: part[0, :6] = radar_partials[sim_idx, :6] + # partial is partial of observable with respect to nominal state partial = part @ stm - partials[partials_idx:partials_idx+size, :] = partial + partials[partials_idx:partials_idx+size, :partial.shape[1]] = partial partials_idx += size return partials diff --git a/src/simulation.cpp b/src/simulation.cpp index 2c75429c..b4ce10c7 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -340,17 +340,26 @@ void IntegBody::prepare_stm(){ * @param[inout] xInteg State of the body. * @param[in] propDir Direction of propagation. */ -void Event::apply_impulsive(const real& t, std::vector& xInteg, - const real& propDir) { +void Event::apply_impulsive(PropSimulation *propSim, const real &t, std::vector& xInteg) { if (t != this->t) { throw std::runtime_error( "Event::apply_impulsive: Integration time does " "not match event time. Cannot apply impulse."); } size_t velStartIdx = this->xIntegIndex + 3; + bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; + real propDir = forwardProp ? 1.0L : -1.0L; for (size_t i = 0; i < 3; i++) { xInteg[velStartIdx + i] += propDir * this->multiplier * this->deltaV[i]; } + if (propSim->integBodies[this->bodyIndex].propStm) { + IntegBody *body = &propSim->integBodies[this->bodyIndex]; + const size_t numNongravs = body->ngParams.a1Est + body->ngParams.a2Est + body->ngParams.a3Est; + const size_t DeltaVStmIdx = this->xIntegIndex + 6 + 36 + 6*numNongravs; + xInteg[DeltaVStmIdx+3] = forwardProp ? 1.0L : -1.0L; + xInteg[DeltaVStmIdx+10] = forwardProp ? 1.0L : -1.0L; + xInteg[DeltaVStmIdx+17] = forwardProp ? 1.0L : -1.0L; + } } /** @@ -748,7 +757,7 @@ void PropSimulation::prepare_for_evaluation( // sort tEval into ascending order or descending order based on the // integration direction (not during orbit fits) if (observerInfo.size() == 0) { - std::vector tEvalSortedIdx(tEval.size()); + std::vector tEvalSortedIdx(tEval.size()); sort_vector(tEval, forwardProp, tEvalSortedIdx); } if (forwardProp) { @@ -986,7 +995,7 @@ void PropSimulation::remove_body(std::string name) { std::cout << "Error: Body " << name << " not found." << std::endl; } -size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, +static size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, const real &tEvent) { // check if tEvent is valid const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; @@ -1016,7 +1025,21 @@ size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, return bodyIndex; } -void event_postprocess(PropSimulation *propSim, Event &event) { +static void event_stm_handling(PropSimulation *propSim, IntegBody *body, const Event &event){ + body->n2Derivs += 3*3; + propSim->integBodies[event.bodyIndex].n2Derivs += 3*3; + propSim->integParams.n2Derivs += 3*3; + + const bool backwardProp = propSim->integParams.tf < propSim->integParams.t0; + std::vector extraVec = std::vector(6*3, 0.0); + // add extra vector at the end of stm vector + for (size_t i = 0; i < extraVec.size(); i++) { + body->stm.push_back(extraVec[i]); + propSim->integBodies[event.bodyIndex].stm.push_back(extraVec[i]); + } +} + +static void event_postprocess(PropSimulation *propSim, Event &event) { // assign event xIntegIndex event.xIntegIndex = 0; for (size_t i = 0; i < event.bodyIndex; i++) { @@ -1059,6 +1082,10 @@ void PropSimulation::add_event(IntegBody body, real tEvent, event.bodyIndex = bodyIndex; event.deltaV = deltaV; event.multiplier = multiplier; + // modify stm if estimating Delta-V + if (body.propStm) { + event_stm_handling(this, &body, event); + } event.tau = 0.0; event.isContinuous = false; From 02da849bba61dcfb0a2626b5c94ea31548500eed Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Wed, 18 Sep 2024 12:01:50 -0500 Subject: [PATCH 12/19] c++ docs update --- include/gr15.h | 5 ----- src/force.cpp | 1 + src/gr15.cpp | 1 + src/simulation.cpp | 2 +- 4 files changed, 3 insertions(+), 6 deletions(-) diff --git a/include/gr15.h b/include/gr15.h index 19aedd54..808b93db 100644 --- a/include/gr15.h +++ b/include/gr15.h @@ -96,16 +96,11 @@ void check_events(PropSimulation *propSim, const real &t, std::vector &xIn /** * @brief Check whether timestep is too large that it would skip an event. - * - * @param[in] propSim PropSimulation object for the integration. - * @param[in] dt Current timestep. */ void event_timestep_check(PropSimulation *propSim, real &dt); /** * @brief 15th-order Gauss-Radau integrator for the PropSimulation. - * - * @param[inout] propSim PropSimulation object for the integration. */ void gr15(PropSimulation *propSim); diff --git a/src/force.cpp b/src/force.cpp index 7ef4df60..f84400ef 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -645,6 +645,7 @@ static void force_thruster(const PropSimulation *propSim, } /** + * @param[in] t Time [TDB MJD] * @param[in] propSim PropSimulation object. * @param[inout] accInteg State derivative vector. */ diff --git a/src/gr15.cpp b/src/gr15.cpp index eb19fc83..4cc73b5c 100644 --- a/src/gr15.cpp +++ b/src/gr15.cpp @@ -289,6 +289,7 @@ void check_events(PropSimulation *propSim, const real &t, std::vector &xIn /** * @param[in] propSim PropSimulation object for the integration. + * @param[in] dt Current timestep. */ void event_timestep_check(PropSimulation *propSim, real &dt) { real tNextEvent = propSim->integParams.tf; diff --git a/src/simulation.cpp b/src/simulation.cpp index b4ce10c7..9862e15a 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -336,9 +336,9 @@ void IntegBody::prepare_stm(){ } /** + * @param[in] propSim PropSimulation object. * @param[in] t Time of the event. * @param[inout] xInteg State of the body. - * @param[in] propDir Direction of propagation. */ void Event::apply_impulsive(PropSimulation *propSim, const real &t, std::vector& xInteg) { if (t != this->t) { From 67a827b1a169e2005e52192c3f2595de5abc44ac Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Thu, 19 Sep 2024 13:21:48 -0500 Subject: [PATCH 13/19] joss paper update --- joss/joss_paper.bib | 28 ++++++++++++++++++++++++++++ joss/joss_paper.md | 8 ++++---- 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/joss/joss_paper.bib b/joss/joss_paper.bib index b1981578..c1f8f1b1 100644 --- a/joss/joss_paper.bib +++ b/joss/joss_paper.bib @@ -83,6 +83,34 @@ @article{Milani1999 author = {Andrea Milani and Giovanni B. Valsecchi}, } +@Article{Farnocchia2019, + author = {Farnocchia, Davide and Eggl, Siegfried and Chodas, Paul W. and Giorgini, Jon D. and Chesley, Steven R.}, + journal = {Celestial Mechanics and Dynamical Astronomy}, + title = {Planetary encounter analysis on the B-plane: a comprehensive formulation}, + year = {2019}, + issn = {1572-9478}, + month = aug, + number = {8}, + volume = {131}, + doi = {10.1007/s10569-019-9914-4}, + publisher = {Springer Science and Business Media LLC}, +} + +@article{Park_2021, +doi = {10.3847/1538-3881/abd414}, +url = {https://dx.doi.org/10.3847/1538-3881/abd414}, +year = {2021}, +month = {feb}, +publisher = {The American Astronomical Society}, +volume = {161}, +number = {3}, +pages = {105}, +author = {Ryan S. Park and William M. Folkner and James G. Williams and Dale H. Boggs}, +title = {The JPL Planetary and Lunar Ephemerides DE440 and DE441}, +journal = {The Astronomical Journal}, +abstract = {The planetary and lunar ephemerides called DE440 and DE441 have been generated by fitting numerically integrated orbits to ground-based and space-based observations. Compared to the previous general-purpose ephemerides DE430, seven years of new data have been added to compute DE440 and DE441, with improved dynamical models and data calibration. The orbit of Jupiter has improved substantially by fitting to the Juno radio range and Very Long Baseline Array (VLBA) data of the Juno spacecraft. The orbit of Saturn has been improved by radio range and VLBA data of the Cassini spacecraft, with improved estimation of the spacecraft orbit. The orbit of Pluto has been improved from use of stellar occultation data reduced against the Gaia star catalog. The ephemerides DE440 and DE441 are fit to the same data set, but DE441 assumes no damping between the lunar liquid core and the solid mantle, which avoids a divergence when integrated backward in time. Therefore, DE441 is less accurate than DE440 for the current century, but covers a much longer duration of years −13,200 to +17,191, compared to DE440 covering years 1550–2650.} +} + @article{Holman2023, doi = {10.3847/PSJ/acc9a9}, year = {2023}, diff --git a/joss/joss_paper.md b/joss/joss_paper.md index bd568555..8995b9c1 100644 --- a/joss/joss_paper.md +++ b/joss/joss_paper.md @@ -34,15 +34,15 @@ aas-journal: Planetary Science Journal # Statement of Need -Understanding the motion of small solar system bodies is of utmost importance when looking at the problem through the lens of planetary defense. The ability to compute the orbit of an asteroid or a comet from various observations and then predicting its trajectory in the future is critical to understanding its impact hazard. The NASA Center for Near-Earth Object Studies (CNEOS) at the Jet Propulsion Laboratory (JPL) has developed a suite of state-of-the-art tools over the course of decades for this specific purpose. However, these tools are not publicly available. With the expected increase in the number of Near-Earth Object (NEO) observations as well as discoveries when new observatories such as the Vera C. Rubin Observatory come online [@Schwamb2023], there is a need for a publicly available library that can reliably perform both orbit determination and propagation for small bodies in the solar system. +Modeling the motion of small solar system bodies is of utmost importance when looking at the problem in the context of planetary defense. Reliably computing the orbit of an asteroid or a comet from various observations, and then predicting its trajectory in the future is critical in assessing the associated Earth-impact hazard. The NASA Center for Near-Earth Object Studies at the Jet Propulsion Laboratory (JPL) has developed a suite of state-of-the-art tools over the course of decades for this specific purpose. However, these tools are not publicly available. With the expected increase in the number of Near-Earth Object (NEO) observations as well as discoveries when new observatories such as the Vera C. Rubin Observatory come online [@Schwamb2023], there is a need for a publicly available library that can reliably perform both orbit determination and propagation for small bodies in the solar system. Such a publicly available library will enable community efforts in planetary defense research. # Summary -In this paper, we present ``GRSS``, the Gauss-Radau Small-body Simulator, an open-source library for orbit determination and propagation of small bodies in the solar system. ``GRSS`` is an open-source, MIT licensed software library with a C++11 foundation and a Python binding. The propagator is based on the ``IAS15`` algorithm, a 15th order integrator based on Gauss-Radau quadrature [@Rein2014]. Only the particles of interest are integrated within ``GRSS`` to reduce computational cost. The states for the planets and Big16 main-belt asteroids are computed using memory-mapped SPICE digital ephemeris kernels as done by @Holman2023 in the ``ASSIST`` orbit propagator library. In addition to the propagator, the C++ portion of the library also has the ability to predict impacts and calculate close encounter circumstances using various formulations of the B-plane [@Kizner1961; @Opik1976; @Chodas1999; @Milani1999]. +In this paper, we present ``GRSS``, the Gauss-Radau Small-body Simulator, an open-source library for orbit determination and propagation of small bodies in the solar system. ``GRSS`` is an open-source, MIT licensed software library with a C++11 foundation and a Python binding. The propagator is based on the ``IAS15`` algorithm, a 15th order integrator based on Gauss-Radau quadrature [@Rein2014]. Only the particles of interest are integrated within ``GRSS`` to reduce computational cost. The states for the planets and 16 largest main-belt asteroids are computed using memory-mapped JPL digital ephemeris kernels [@Park2021] as done in the ``ASSIST`` orbit propagator library [@Holman2023]. In addition to the propagator, the C++ portion of the library also has the ability to predict impacts and calculate close encounter circumstances using various formulations of the B-plane [@Kizner1961; @Opik1976; @Chodas1999; @Milani1999; @Farnocchia2019]. The C++ functionality is exposed to Python through a binding generated using ``pybind11``[^1]. The Python layer of ``GRSS`` uses the propagator as the foundation to compute the orbits of small bodies from a given set of optical and/or radar astrometry from the Minor Planet Center[^2], the JPL Small Body Radar Astrometry database[^3], and the Gaia Focused Product Release (FPR) solar system observations database[^4]. Additionally, the orbit determination modules also have the ability to fit especially demanding measurements such as stellar occultations. These capabilities of the ``GRSS`` library have already been used to study the the heliocentric changes in the orbit of the (65803) Didymos binary asteroid system as a result of the DART impact [@Makadia2022; @Makadia2024]. -``GRSS`` will continue to be developed in the future, with anticipated contributions including the ability to perform mission studies for future asteroid deflections. ``GRSS`` is publicly available to the community through the Python Package Index (PyPI) and the source code is available on GitHub[^5]. Therefore, ``GRSS`` will allow the research community to have access to a reliable and efficient tool for studying the dynamics of small bodies in the solar system. +``GRSS`` will continue to be developed in the future, with anticipated contributions including the ability to perform mission studies for future asteroid deflections. ``GRSS`` is publicly available to the community through the Python Package Index (PyPI) and the source code is available on GitHub[^5]. Therefore, ``GRSS`` is a reliable and efficient tool that the community has access to for studying the dynamics of small bodies in the solar system. [^1]: [^2]: @@ -52,6 +52,6 @@ The C++ functionality is exposed to Python through a binding generated using ``p # Acknowledgements -R.M. acknowledges funding from a NASA Space Technology Graduate Research Opportunities (NSTGRO) award, NASA contract No. 80NSSC22K1173. The work of S.R.C. and D.F. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (No. 80NM0018D0004). This work has made use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium. This library has also extensively used data and services provided by the International Astronomical Union Minor Planet Center. Data from the MPC's database is made freely available to the public. Funding for this data and the MPC's operations comes from a NASA PDCO grant (80NSSC22M0024), administered via a University of Maryland - SAO subaward (106075-Z6415201). The MPC's computing equipment is funded in part by the above award, and in part by funding from the Tamkin Foundation. +R.M. acknowledges funding from a NASA Space Technology Graduate Research Opportunities (NSTGRO) award, NASA contract No. 80NSSC22K1173. The work of S.R.C. and D.F. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (No. 80NM0018D0004). This work has made use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium. This library has also extensively used data and services provided by the International Astronomical Union Minor Planet Center (MPC). Data from the MPC's database is made freely available to the public. Funding for this data and the MPC's operations comes from a NASA PDCO grant (80NSSC22M0024), administered via a University of Maryland - SAO subaward (106075-Z6415201). The MPC's computing equipment is funded in part by the above award, and in part by funding from the Tamkin Foundation. # References From 542eeea2761b9dc85644ec8750661e173d7725e9 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Thu, 19 Sep 2024 13:25:20 -0500 Subject: [PATCH 14/19] joss bib update --- joss/joss_paper.bib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/joss/joss_paper.bib b/joss/joss_paper.bib index c1f8f1b1..7418da9f 100644 --- a/joss/joss_paper.bib +++ b/joss/joss_paper.bib @@ -96,7 +96,7 @@ @Article{Farnocchia2019 publisher = {Springer Science and Business Media LLC}, } -@article{Park_2021, +@article{Park2021, doi = {10.3847/1538-3881/abd414}, url = {https://dx.doi.org/10.3847/1538-3881/abd414}, year = {2021}, From d80ff24d8d4317dbeb7905b9e2cc5e778375f572 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Thu, 19 Sep 2024 14:24:54 -0500 Subject: [PATCH 15/19] added continuous exp accel estimation; joss paper edits --- grss/fit/fit_simulation.py | 30 ++++++++++++++++++++---------- joss/joss_paper.md | 10 +++++----- 2 files changed, 25 insertions(+), 15 deletions(-) diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index 6aa92dc8..1490b090 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -1028,15 +1028,23 @@ def _x_dict_to_events(self, x_dict): events = [] for i in range(len(self.fixed_propsim_params['events'])): fixed_event = tuple(self.fixed_propsim_params['events'][i]) + fixed_event_impulsive = len(fixed_event) == 5 + fixed_event_continuous = len(fixed_event) == 6 + if fixed_event_continuous and fixed_event[5] == 0.0: + raise ValueError("Continuous event time constant cannot be 0.") + assert fixed_event_impulsive or fixed_event_continuous event = [None]*5 event[0] = fixed_event[0] # time of delta-v - event[1] = x_dict[f"dvx{i}"] if f"dvx{i}" in x_dict.keys() else fixed_event[1] - event[2] = x_dict[f"dvy{i}"] if f"dvy{i}" in x_dict.keys() else fixed_event[2] - event[3] = x_dict[f"dvz{i}"] if f"dvz{i}" in x_dict.keys() else fixed_event[3] - event[4] = x_dict[f"mult{i}"] if f"mult{i}" in x_dict.keys() else fixed_event[4] - fixed_event_continuous = len(fixed_event) == 6 and fixed_event[5] != 0.0 - if fixed_event_continuous or f"dt{i}" in x_dict.keys(): + if fixed_event_impulsive: + event[1] = x_dict[f"dvx{i}"] if f"dvx{i}" in x_dict.keys() else fixed_event[1] + event[2] = x_dict[f"dvy{i}"] if f"dvy{i}" in x_dict.keys() else fixed_event[2] + event[3] = x_dict[f"dvz{i}"] if f"dvz{i}" in x_dict.keys() else fixed_event[3] + elif fixed_event_continuous: + event[1] = x_dict[f"ax{i}"] if f"ax{i}" in x_dict.keys() else fixed_event[1] + event[2] = x_dict[f"ay{i}"] if f"ay{i}" in x_dict.keys() else fixed_event[2] + event[3] = x_dict[f"az{i}"] if f"az{i}" in x_dict.keys() else fixed_event[3] event.append(x_dict[f"dt{i}"] if f"dt{i}" in x_dict.keys() else fixed_event[5]) + event[4] = x_dict[f"mult{i}"] if f"mult{i}" in x_dict.keys() else fixed_event[4] events.append(tuple(event)) return events @@ -1064,11 +1072,11 @@ def _check_and_add_events(self, prop_sim_past, prop_sim_future, integ_body, even """ for event in events: t_event = event[0] - dvx = event[1] - dvy = event[2] - dvz = event[3] + dvx_or_ax = event[1] + dvy_or_ay = event[2] + dvz_or_az = event[3] multiplier = event[4] - event_args = (integ_body, t_event, [dvx, dvy, dvz], multiplier) + event_args = (integ_body, t_event, [dvx_or_ax, dvy_or_ay, dvz_or_az], multiplier) if len(event) == 6: dt = event[5] event_args += (dt,) @@ -1124,6 +1132,8 @@ def _get_perturbed_state(self, key): fd_pert = 1e0 if key[:3] in {'dvx', 'dvy', 'dvz'}: fd_pert = 1e-11 + if key[:2] in {'ax', 'ay', 'az'}: + fd_pert = 1e-11 if key[:2] == 'dt': fd_pert = 1.0 if fd_pert == np.inf: diff --git a/joss/joss_paper.md b/joss/joss_paper.md index 8995b9c1..d1795413 100644 --- a/joss/joss_paper.md +++ b/joss/joss_paper.md @@ -34,15 +34,15 @@ aas-journal: Planetary Science Journal # Statement of Need -Modeling the motion of small solar system bodies is of utmost importance when looking at the problem in the context of planetary defense. Reliably computing the orbit of an asteroid or a comet from various observations, and then predicting its trajectory in the future is critical in assessing the associated Earth-impact hazard. The NASA Center for Near-Earth Object Studies at the Jet Propulsion Laboratory (JPL) has developed a suite of state-of-the-art tools over the course of decades for this specific purpose. However, these tools are not publicly available. With the expected increase in the number of Near-Earth Object (NEO) observations as well as discoveries when new observatories such as the Vera C. Rubin Observatory come online [@Schwamb2023], there is a need for a publicly available library that can reliably perform both orbit determination and propagation for small bodies in the solar system. Such a publicly available library will enable community efforts in planetary defense research. +Modeling the motion of small solar system bodies is of utmost importance when looking at the problem in the context of planetary defense. Reliably computing the orbit of an asteroid or a comet from various observations, and then predicting its trajectory in the future is critical in assessing the associated Earth-impact hazard. The NASA Center for Near-Earth Object Studies at the Jet Propulsion Laboratory (JPL) has developed a suite of state-of-the-art tools over the course of decades for this specific purpose. However, these tools are not publicly available. With the expected increase in the number of Near-Earth Object observations as well as discoveries when new observatories such as the Vera C. Rubin Observatory come online [@Schwamb2023], there is a need for a publicly available library that can reliably perform both orbit determination and propagation for small bodies in the solar system. Such a publicly available library will enable community efforts in planetary defense research. # Summary -In this paper, we present ``GRSS``, the Gauss-Radau Small-body Simulator, an open-source library for orbit determination and propagation of small bodies in the solar system. ``GRSS`` is an open-source, MIT licensed software library with a C++11 foundation and a Python binding. The propagator is based on the ``IAS15`` algorithm, a 15th order integrator based on Gauss-Radau quadrature [@Rein2014]. Only the particles of interest are integrated within ``GRSS`` to reduce computational cost. The states for the planets and 16 largest main-belt asteroids are computed using memory-mapped JPL digital ephemeris kernels [@Park2021] as done in the ``ASSIST`` orbit propagator library [@Holman2023]. In addition to the propagator, the C++ portion of the library also has the ability to predict impacts and calculate close encounter circumstances using various formulations of the B-plane [@Kizner1961; @Opik1976; @Chodas1999; @Milani1999; @Farnocchia2019]. +In this paper, we present ``GRSS``, the Gauss-Radau Small-body Simulator, an open-source library for orbit determination and propagation of small bodies in the solar system. ``GRSS`` is an open-source software library with a C++11 foundation and a Python binding. The propagator is based on the ``IAS15`` algorithm, a 15th order integrator based on Gauss-Radau quadrature [@Rein2014]. Only the particles of interest are integrated within ``GRSS`` to reduce computational cost. The states for the planets and 16 largest main-belt asteroids are computed using memory-mapped JPL digital ephemeris kernels [@Park2021] as done in the ``ASSIST`` orbit propagator library [@Holman2023]. In addition to the propagator, the C++ portion of the library also has the ability to predict impacts and calculate close encounter circumstances using various formulations of the B-plane [@Kizner1961; @Opik1976; @Chodas1999; @Milani1999; @Farnocchia2019]. -The C++ functionality is exposed to Python through a binding generated using ``pybind11``[^1]. The Python layer of ``GRSS`` uses the propagator as the foundation to compute the orbits of small bodies from a given set of optical and/or radar astrometry from the Minor Planet Center[^2], the JPL Small Body Radar Astrometry database[^3], and the Gaia Focused Product Release (FPR) solar system observations database[^4]. Additionally, the orbit determination modules also have the ability to fit especially demanding measurements such as stellar occultations. These capabilities of the ``GRSS`` library have already been used to study the the heliocentric changes in the orbit of the (65803) Didymos binary asteroid system as a result of the DART impact [@Makadia2022; @Makadia2024]. +The C++ functionality is exposed to Python through a binding generated using ``pybind11``[^1]. The Python layer of ``GRSS`` uses the propagator as the foundation to compute the orbits of small bodies from a given set of optical and/or radar astrometry from the Minor Planet Center[^2], the JPL Small Body Radar Astrometry database[^3], and the Gaia Focused Product Release solar system observations database[^4]. Additionally, the orbit determination modules also have the ability to fit especially demanding measurements such as stellar occultations. These capabilities of the ``GRSS`` library have already been used to study the the heliocentric changes in the orbit of the (65803) Didymos binary asteroid system as a result of the DART impact [@Makadia2022; @Makadia2024]. -``GRSS`` will continue to be developed in the future, with anticipated contributions including the ability to perform mission studies for future asteroid deflections. ``GRSS`` is publicly available to the community through the Python Package Index (PyPI) and the source code is available on GitHub[^5]. Therefore, ``GRSS`` is a reliable and efficient tool that the community has access to for studying the dynamics of small bodies in the solar system. +``GRSS`` will continue to be developed in the future, with anticipated contributions including the ability to perform mission studies for future asteroid deflections. ``GRSS`` is publicly available to the community through the Python Package Index and the source code is available on GitHub[^5]. Therefore, ``GRSS`` is a reliable and efficient tool that the community has access to for studying the dynamics of small bodies in the solar system. [^1]: [^2]: @@ -52,6 +52,6 @@ The C++ functionality is exposed to Python through a binding generated using ``p # Acknowledgements -R.M. acknowledges funding from a NASA Space Technology Graduate Research Opportunities (NSTGRO) award, NASA contract No. 80NSSC22K1173. The work of S.R.C. and D.F. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (No. 80NM0018D0004). This work has made use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium. This library has also extensively used data and services provided by the International Astronomical Union Minor Planet Center (MPC). Data from the MPC's database is made freely available to the public. Funding for this data and the MPC's operations comes from a NASA PDCO grant (80NSSC22M0024), administered via a University of Maryland - SAO subaward (106075-Z6415201). The MPC's computing equipment is funded in part by the above award, and in part by funding from the Tamkin Foundation. +R.M. acknowledges funding from a NASA Space Technology Graduate Research Opportunities award, NASA contract No. 80NSSC22K1173. The work of S.R.C. and D.F. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (No. 80NM0018D0004). This work has made use of data from the European Space Agency mission Gaia, processed by the Gaia Data Processing and Analysis Consortium. This library has also extensively used data and services provided by the International Astronomical Union Minor Planet Center (MPC). Data from the MPC's database is made freely available to the public. Funding for this data and the MPC's operations comes from a NASA Planetary Defense Coordination Office grant (80NSSC22M0024), administered via a University of Maryland - Smithsonian Astrophysical Observatory subaward (106075-Z6415201). The MPC's computing equipment is funded in part by the above award, and in part by funding from the Tamkin Foundation. # References From ffc8868bd43b9f21c755fd367da985e907902cf6 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Fri, 20 Sep 2024 11:58:49 -0500 Subject: [PATCH 16/19] added analytic multiplier estimation for impulsive events --- include/simulation.h | 1 + src/simulation.cpp | 35 ++++++++++++++++++++++++++++------- 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/include/simulation.h b/include/simulation.h index 8ad65c8e..6f91e1a4 100644 --- a/include/simulation.h +++ b/include/simulation.h @@ -240,6 +240,7 @@ class Event { // for impulsive events std::vector deltaV = {0.0L, 0.0L, 0.0L}; real multiplier = 1.0L; + bool multiplierEst = false; // for continuous ejecta events std::vector expAccel0 = {0.0L, 0.0L, 0.0L}; diff --git a/src/simulation.cpp b/src/simulation.cpp index 9862e15a..6d56b838 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -356,9 +356,15 @@ void Event::apply_impulsive(PropSimulation *propSim, const real &t, std::vector< IntegBody *body = &propSim->integBodies[this->bodyIndex]; const size_t numNongravs = body->ngParams.a1Est + body->ngParams.a2Est + body->ngParams.a3Est; const size_t DeltaVStmIdx = this->xIntegIndex + 6 + 36 + 6*numNongravs; - xInteg[DeltaVStmIdx+3] = forwardProp ? 1.0L : -1.0L; - xInteg[DeltaVStmIdx+10] = forwardProp ? 1.0L : -1.0L; - xInteg[DeltaVStmIdx+17] = forwardProp ? 1.0L : -1.0L; + if (this->multiplierEst){ + xInteg[DeltaVStmIdx+3] = propDir * this->deltaV[0]; + xInteg[DeltaVStmIdx+4] = propDir * this->deltaV[1]; + xInteg[DeltaVStmIdx+5] = propDir * this->deltaV[2]; + } else { + xInteg[DeltaVStmIdx+3] = propDir; + xInteg[DeltaVStmIdx+10] = propDir; + xInteg[DeltaVStmIdx+17] = propDir; + } } } @@ -1026,12 +1032,16 @@ static size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, } static void event_stm_handling(PropSimulation *propSim, IntegBody *body, const Event &event){ - body->n2Derivs += 3*3; - propSim->integBodies[event.bodyIndex].n2Derivs += 3*3; - propSim->integParams.n2Derivs += 3*3; + int numEventParams = 3; + if (event.multiplierEst) { + numEventParams = 1; + } + body->n2Derivs += 3*numEventParams; + propSim->integBodies[event.bodyIndex].n2Derivs += 3*numEventParams; + propSim->integParams.n2Derivs += 3*numEventParams; const bool backwardProp = propSim->integParams.tf < propSim->integParams.t0; - std::vector extraVec = std::vector(6*3, 0.0); + std::vector extraVec = std::vector(6*numEventParams, 0.0); // add extra vector at the end of stm vector for (size_t i = 0; i < extraVec.size(); i++) { body->stm.push_back(extraVec[i]); @@ -1081,6 +1091,17 @@ void PropSimulation::add_event(IntegBody body, real tEvent, event.bodyName = body.name; event.bodyIndex = bodyIndex; event.deltaV = deltaV; + // if multiplier is numeric, estimate it otherwise estimate full deltaV vector + if (!std::isfinite(multiplier)) { + if (body.propStm){ + event.multiplierEst = false; + } + multiplier = 1.0; + } else { + if (body.propStm) { + event.multiplierEst = true; + } + } event.multiplier = multiplier; // modify stm if estimating Delta-V if (body.propStm) { From 83018e23daba1d7c2d2cce7d6cc5ec81041655ba Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Fri, 20 Sep 2024 12:06:03 -0500 Subject: [PATCH 17/19] added internet connectivity check for downloading files --- grss/utils.py | 13 +++++++++++++ grss/version.txt | 2 +- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/grss/utils.py b/grss/utils.py index 734a7e13..1dad90f1 100644 --- a/grss/utils.py +++ b/grss/utils.py @@ -2,6 +2,7 @@ import os import time import gzip +import requests __all__ = [ 'default_kernel_path', 'grss_path', @@ -11,6 +12,14 @@ grss_path = os.path.dirname(os.path.abspath(__file__)) default_kernel_path = f'{grss_path}/kernels/' +def _connected_to_internet(url='http://www.google.com/', timeout=5): + try: + _ = requests.head(url, timeout=timeout) + return True + except requests.ConnectionError: + print("No internet connection available.") + return False + def _download_ades_files(file): """ Download a text file of the observatory code info from the Minor Planet Center. @@ -66,6 +75,10 @@ def initialize(): None : NoneType None """ + internet = _connected_to_internet() + if not internet: + print("Please connect to the internet to download/update the necessary data files.") + return None # run the get_debiasing_data.py script os.system(f'python {grss_path}/debias/get_debiasing_data.py') # run the get_kernels.py script diff --git a/grss/version.txt b/grss/version.txt index 81911389..ecedc98d 100644 --- a/grss/version.txt +++ b/grss/version.txt @@ -1 +1 @@ -4.3.0 \ No newline at end of file +4.3.1 \ No newline at end of file From 0f1090bf5c4f39179aa84ae1dfb4307b3af3ca55 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Fri, 20 Sep 2024 18:51:23 -0500 Subject: [PATCH 18/19] con. event timestep control; log file minor updates --- grss/fit/fit_simulation.py | 39 +++++++++++++++++++++----------------- src/gr15.cpp | 12 ++++++++++++ 2 files changed, 34 insertions(+), 17 deletions(-) diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index 1490b090..b8f352fe 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -1135,7 +1135,8 @@ def _get_perturbed_state(self, key): if key[:2] in {'ax', 'ay', 'az'}: fd_pert = 1e-11 if key[:2] == 'dt': - fd_pert = 1.0 + # 10 percent of the time constant + fd_pert = 0.1*self.x_nom[key] if fd_pert == np.inf: raise ValueError("Finite difference perturbation not defined" " for the given state parameter.") @@ -2002,6 +2003,8 @@ def save(self, filename): 'w': 'deg', 'i': 'deg', 'x': 'AU', 'y': 'AU', 'z': 'AU', 'vx': 'AU/day', 'vy': 'AU/day', 'vz': 'AU/day', 'a1': 'AU/day^2', 'a2': 'AU/day^2', 'a3': 'AU/day^2', + 'dv': 'AU/day', 'ax': 'AU/day^2', 'ay': 'AU/day^2', 'az': 'AU/day^2', + 'dt': 'day' } with open(filename, 'x', encoding='utf-8') as f: f.write(section_full + '\n') @@ -2040,24 +2043,25 @@ def save(self, filename): if key in ['om', 'w', 'i']: val *= 180/np.pi sig *= 180/np.pi - if key in ['a1', 'a2', 'a3']: + if key[:2] in {'a1', 'a2', 'a3', 'dv', 'ax', 'ay', 'az'}: f.write(f"{key.upper()}: {val:.11e} \u00B1 {sig:.11e}") else: f.write(f"{key.upper()}: {val:.11f} \u00B1 {sig:.11f}") - if key in units_dict: - f.write(f" {units_dict[key]}") + if key[:2] in units_dict: + f.write(f" {units_dict[key[:2]]}") if self._priors_given and key in self.prior_est: p_val = self.prior_est[key] p_sig = self.prior_sigs[key] if key in ['om', 'w', 'i']: p_val *= 180/np.pi p_sig *= 180/np.pi - if key in ['a1', 'a2', 'a3']: - f.write(f" (a priori: {p_val:.11e} \u00B1 {p_sig:.11e})") + if key[:2] in {'a1', 'a2', 'a3', 'dv', 'ax', 'ay', 'az'}: + f.write(f" (a priori: {p_val:.11e} \u00B1 {p_sig:.11e}") else: - f.write(f" (a priori: {p_val:.11f} \u00B1 {p_sig:.11f})") - if key in units_dict: - f.write(f" {units_dict[key]}") + f.write(f" (a priori: {p_val:.11f} \u00B1 {p_sig:.11f}") + if key[:2] in units_dict: + f.write(f" {units_dict[key[:2]]}") + f.write(")") f.write("\n") f.write("\n") f.write("Initial Covariance Matrix:\n") @@ -2079,24 +2083,25 @@ def save(self, filename): if key in ['om', 'w', 'i']: val *= 180/np.pi sig *= 180/np.pi - if key in ['a1', 'a2', 'a3']: + if key[:2] in {'a1', 'a2', 'a3', 'dv', 'ax', 'ay', 'az'}: f.write(f"{key.upper()}: {val:.11e} \u00B1 {sig:.11e}") else: f.write(f"{key.upper()}: {val:.11f} \u00B1 {sig:.11f}") - if key in units_dict: - f.write(f" {units_dict[key]}") + if key[:2] in units_dict: + f.write(f" {units_dict[key[:2]]}") if self._priors_given and key in self.prior_est: p_val = self.prior_est[key] p_sig = self.prior_sigs[key] if key in ['om', 'w', 'i']: p_val *= 180/np.pi p_sig *= 180/np.pi - if key in ['a1', 'a2', 'a3']: - f.write(f" (a priori: {p_val:.11e} \u00B1 {p_sig:.11e})") + if key[:2] in {'a1', 'a2', 'a3', 'dv', 'ax', 'ay', 'az'}: + f.write(f" (a priori: {p_val:.11e} \u00B1 {p_sig:.11e}") else: - f.write(f" (a priori: {p_val:.11f} \u00B1 {p_sig:.11f})") - if key in units_dict: - f.write(f" {units_dict[key]}") + f.write(f" (a priori: {p_val:.11f} \u00B1 {p_sig:.11f}") + if key[:2] in units_dict: + f.write(f" {units_dict[key[:2]]}") + f.write(")") f.write("\n") f.write("\n") f.write("Final Covariance Matrix:\n") diff --git a/src/gr15.cpp b/src/gr15.cpp index 4cc73b5c..6fd4fd20 100644 --- a/src/gr15.cpp +++ b/src/gr15.cpp @@ -310,6 +310,18 @@ void event_timestep_check(PropSimulation *propSim, real &dt) { } else if (!forwardProp && propSim->t > tNextConEvent) { tNextEvent = fmax(tNextEvent, tNextConEvent); } + // break up the timestep if early in a continuous event + if (propSim->eventMngr.continuousEvents[i].isHappening) { + const real timeConstant = propSim->eventMngr.continuousEvents[i].tau; + const real timeConstantFac = 5.0L; + const real numSegments = 100.0L; + const real dtSegment = timeConstantFac*timeConstant/numSegments; + if (forwardProp && propSim->t < tNextConEvent+timeConstantFac*timeConstant) { + dt = fmin(dt, dtSegment); + } else if (!forwardProp && propSim->t < tNextConEvent+timeConstantFac*timeConstant) { + dt = fmax(dt, -dtSegment); + } + } } } if ((forwardProp && propSim->t + dt > tNextEvent) || From 89dac9d4c04f7af0ef1e272420b3e9d20e75f692 Mon Sep 17 00:00:00 2001 From: Rahil Makadia Date: Mon, 23 Sep 2024 17:08:13 -0500 Subject: [PATCH 19/19] updated method for adding events; continuous events analytic estimation --- grss/fit/fit_simulation.py | 53 ++++++++++---- grss/version.txt | 2 +- include/simulation.h | 30 ++++---- include/stm.h | 6 ++ src/force.cpp | 26 ++++--- src/gr15.cpp | 6 +- src/grss.cpp | 113 ++++++++++++++++------------- src/simulation.cpp | 143 +++++++++++++++++++++---------------- src/stm.cpp | 51 +++++++++++++ 9 files changed, 273 insertions(+), 157 deletions(-) diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py index b8f352fe..7e926f52 100644 --- a/grss/fit/fit_simulation.py +++ b/grss/fit/fit_simulation.py @@ -1030,7 +1030,7 @@ def _x_dict_to_events(self, x_dict): fixed_event = tuple(self.fixed_propsim_params['events'][i]) fixed_event_impulsive = len(fixed_event) == 5 fixed_event_continuous = len(fixed_event) == 6 - if fixed_event_continuous and fixed_event[5] == 0.0: + if fixed_event_continuous and fixed_event[4] == 0.0: raise ValueError("Continuous event time constant cannot be 0.") assert fixed_event_impulsive or fixed_event_continuous event = [None]*5 @@ -1039,12 +1039,13 @@ def _x_dict_to_events(self, x_dict): event[1] = x_dict[f"dvx{i}"] if f"dvx{i}" in x_dict.keys() else fixed_event[1] event[2] = x_dict[f"dvy{i}"] if f"dvy{i}" in x_dict.keys() else fixed_event[2] event[3] = x_dict[f"dvz{i}"] if f"dvz{i}" in x_dict.keys() else fixed_event[3] + event[4] = x_dict[f"mult{i}"] if f"mult{i}" in x_dict.keys() else fixed_event[4] elif fixed_event_continuous: event[1] = x_dict[f"ax{i}"] if f"ax{i}" in x_dict.keys() else fixed_event[1] event[2] = x_dict[f"ay{i}"] if f"ay{i}" in x_dict.keys() else fixed_event[2] event[3] = x_dict[f"az{i}"] if f"az{i}" in x_dict.keys() else fixed_event[3] - event.append(x_dict[f"dt{i}"] if f"dt{i}" in x_dict.keys() else fixed_event[5]) - event[4] = x_dict[f"mult{i}"] if f"mult{i}" in x_dict.keys() else fixed_event[4] + event[4] = x_dict[f"dt{i}"] if f"dt{i}" in x_dict.keys() else fixed_event[4] + event.append(fixed_event[5]) events.append(tuple(event)) return events @@ -1070,20 +1071,42 @@ def _check_and_add_events(self, prop_sim_past, prop_sim_future, integ_body, even prop_sim_future : libgrss.PropSimulation object PropSimulation object for the future. """ - for event in events: - t_event = event[0] - dvx_or_ax = event[1] - dvy_or_ay = event[2] - dvz_or_az = event[3] - multiplier = event[4] - event_args = (integ_body, t_event, [dvx_or_ax, dvy_or_ay, dvz_or_az], multiplier) - if len(event) == 6: - dt = event[5] - event_args += (dt,) + est_keys = self.x_nom.keys() + for i, event_list in enumerate(events): + t_event = event_list[0] + dvx_or_ax = event_list[1] + dvy_or_ay = event_list[2] + dvz_or_az = event_list[3] + multiplier_or_dt = event_list[4] + continuous_flag = event_list[5] if len(event_list) == 6 else False + event = libgrss.Event() + event.t = t_event + event.bodyName = integ_body.name + event.isContinuous = continuous_flag + if continuous_flag: + event.expAccel0 = [dvx_or_ax, dvy_or_ay, dvz_or_az] + event.tau = multiplier_or_dt + event.expAccel0Est = ( + f"ax{i}" in est_keys and + f"ay{i}" in est_keys and + f"az{i}" in est_keys + ) + event.tauEst = f"dt{i}" in est_keys + event.eventEst = event.expAccel0Est or event.tauEst + else: + event.deltaV = [dvx_or_ax, dvy_or_ay, dvz_or_az] + event.multiplier = multiplier_or_dt + event.deltaVEst = ( + f"dvx{i}" in est_keys and + f"dvy{i}" in est_keys and + f"dvz{i}" in est_keys + ) + event.multiplierEst = f"mult{i}" in est_keys + event.eventEst = event.deltaVEst or event.multiplierEst if t_event < self.t_sol: - prop_sim_past.add_event(*event_args) + prop_sim_past.add_event(event) else: - prop_sim_future.add_event(*event_args) + prop_sim_future.add_event(event) return prop_sim_past, prop_sim_future def _get_perturbed_state(self, key): diff --git a/grss/version.txt b/grss/version.txt index ecedc98d..7e961f9e 100644 --- a/grss/version.txt +++ b/grss/version.txt @@ -1 +1 @@ -4.3.1 \ No newline at end of file +4.3.2 \ No newline at end of file diff --git a/include/simulation.h b/include/simulation.h index 6f91e1a4..9fe18d45 100644 --- a/include/simulation.h +++ b/include/simulation.h @@ -234,19 +234,27 @@ class Event { public: real t; std::string bodyName; + bool isContinuous = false; + bool eventEst = false; size_t bodyIndex; size_t xIntegIndex; + bool hasStarted = false; // for impulsive events - std::vector deltaV = {0.0L, 0.0L, 0.0L}; - real multiplier = 1.0L; + std::vector deltaV = std::vector(3, std::numeric_limits::quiet_NaN()); + real multiplier = std::numeric_limits::quiet_NaN(); + bool deltaVEst = false; bool multiplierEst = false; // for continuous ejecta events - std::vector expAccel0 = {0.0L, 0.0L, 0.0L}; - real tau = 1.0L; - bool isContinuous = false; - bool isHappening = false; + std::vector expAccel0 = std::vector(3, std::numeric_limits::quiet_NaN()); + real tau = std::numeric_limits::quiet_NaN(); + bool expAccel0Est = false; + bool tauEst = false; + /** + * @brief Empty constructor for the Event class. + */ + Event() {}; /** * @brief Apply the impulse event to the body. */ @@ -517,15 +525,9 @@ class PropSimulation { */ void remove_body(std::string name); /** - * @brief Add an impulse event to the simulation. - */ - void add_event(IntegBody body, real tEvent, std::vector deltaV, - real multiplier); - /** - * @brief Add an ejecta event to the simulation. + * @brief Add an event to the simulation. */ - void add_event(IntegBody body, real tEvent, std::vector expAccel0, - real multiplier, real tau); + void add_event(Event event); /** * @brief Set the values of the PropSimulation Constants object. */ diff --git a/include/stm.h b/include/stm.h index acbfaf05..22791863 100644 --- a/include/stm.h +++ b/include/stm.h @@ -71,4 +71,10 @@ void stm_nongrav(STMParameters &stmParams, const real &g, const real &dy, const real &dz, const real &dvx, const real &dvy, const real &dvz, real *rVec, real *nVec); +/** + * @brief Compute the derivatives of a continuous event with respect to position, velocity, and parameters. + */ +void stm_continuous_event(STMParameters &stmParams, + const PropSimulation *propSim, const size_t &eventIdx, + const real &tPastEvent, const real &postFac); #endif diff --git a/src/force.cpp b/src/force.cpp index f84400ef..f34ea7f9 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -41,7 +41,7 @@ static void force_thruster(const PropSimulation *propSim, std::vector &acc * @brief Compute the acceleration of the system due to a continuous event. */ static void force_continuous_event(const real &t, const PropSimulation *propSim, - std::vector &accInteg); + std::vector &accInteg, STMParameters* allSTMs); /** * @param[in] t Time [TDB MJD] @@ -121,7 +121,7 @@ void get_state_der(PropSimulation *propSim, const real &t, force_J2(propSim, accInteg, allSTMs); force_nongrav(propSim, accInteg, allSTMs); force_thruster(propSim, accInteg); - force_continuous_event(t, propSim, accInteg); + force_continuous_event(t, propSim, accInteg, allSTMs); #ifdef PRINT_FORCES forceFile.open("cpp.11", std::ios::app); forceFile << std::setw(10) << "total_acc" << std::setw(25) << accInteg[0] @@ -648,27 +648,31 @@ static void force_thruster(const PropSimulation *propSim, * @param[in] t Time [TDB MJD] * @param[in] propSim PropSimulation object. * @param[inout] accInteg State derivative vector. + * @param[in] allSTMs STMParameters vector for IntegBodies in the simulation. */ static void force_continuous_event(const real &t, const PropSimulation *propSim, - std::vector &accInteg) { + std::vector &accInteg, STMParameters* allSTMs) { if (!propSim->eventMngr.allConEventDone){ const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; const real propDir = forwardProp ? 1.0 : -1.0; for (size_t i = 0; i < propSim->eventMngr.continuousEvents.size(); i++){ - if (propSim->eventMngr.continuousEvents[i].isHappening){ + if (propSim->eventMngr.continuousEvents[i].hasStarted){ const size_t starti = propSim->eventMngr.continuousEvents[i].xIntegIndex / 2; const real tPastEvent = t - propSim->eventMngr.continuousEvents[i].t; - const real postFac = exp(-tPastEvent/propSim->eventMngr.continuousEvents[i].tau); + const real postFac = exp(-tPastEvent/propSim->eventMngr.continuousEvents[i].tau) * propDir; accInteg[starti + 0] += propSim->eventMngr.continuousEvents[i].expAccel0[0] * - postFac * - propSim->eventMngr.continuousEvents[i].multiplier * propDir; + postFac; accInteg[starti + 1] += propSim->eventMngr.continuousEvents[i].expAccel0[1] * - postFac * - propSim->eventMngr.continuousEvents[i].multiplier * propDir; + postFac; accInteg[starti + 2] += propSim->eventMngr.continuousEvents[i].expAccel0[2] * - postFac * - propSim->eventMngr.continuousEvents[i].multiplier * propDir; + postFac; + const size_t bodyIdx = propSim->eventMngr.continuousEvents[i].bodyIndex; + if (propSim->integBodies[bodyIdx].propStm && + propSim->eventMngr.continuousEvents[i].eventEst) { + stm_continuous_event(allSTMs[bodyIdx], propSim, i, + tPastEvent, postFac); + } } } } diff --git a/src/gr15.cpp b/src/gr15.cpp index 6fd4fd20..a2c8d322 100644 --- a/src/gr15.cpp +++ b/src/gr15.cpp @@ -262,9 +262,9 @@ void check_continuous_events(PropSimulation *propSim, const real &t) { eventEnded = t < propSim->eventMngr.continuousEvents[i].t; } if (eventStarted && !eventEnded) { - propSim->eventMngr.continuousEvents[i].isHappening = true; + propSim->eventMngr.continuousEvents[i].hasStarted = true; } else { - propSim->eventMngr.continuousEvents[i].isHappening = false; + propSim->eventMngr.continuousEvents[i].hasStarted = false; } if (!eventEnded) { allDone = false; @@ -311,7 +311,7 @@ void event_timestep_check(PropSimulation *propSim, real &dt) { tNextEvent = fmax(tNextEvent, tNextConEvent); } // break up the timestep if early in a continuous event - if (propSim->eventMngr.continuousEvents[i].isHappening) { + if (propSim->eventMngr.continuousEvents[i].hasStarted) { const real timeConstant = propSim->eventMngr.continuousEvents[i].tau; const real timeConstantFac = 5.0L; const real numSegments = 100.0L; diff --git a/src/grss.cpp b/src/grss.cpp index 2a51cdd1..1c082cee 100644 --- a/src/grss.cpp +++ b/src/grss.cpp @@ -555,7 +555,7 @@ PYBIND11_MODULE(libgrss, m) { Radius of the body. cometaryState : list of real Initial Heliocentric Ecliptic Cometary state of the body. - ngParams : PropSimulation.NongravParameters + ngParams : libgrss.NongravParameters Non-gravitational parameters of the body. )mydelimiter") @@ -577,7 +577,7 @@ PYBIND11_MODULE(libgrss, m) { Initial barycentric Cartesian position of the body. vel : list of real Initial barycentric Cartesian velocity of the body. - ngParams : PropSimulation.NongravParameters + ngParams : libgrss.NongravParameters Non-gravitational parameters of the body. )mydelimiter") .def_readwrite("spiceId", &IntegBody::spiceId, R"mydelimiter( @@ -622,6 +622,52 @@ PYBIND11_MODULE(libgrss, m) { None. )mydelimiter"); + // expose event class and its empty constructor + py::class_(m, "Event", R"mydelimiter( + The Event class contains the properties of an impulsive or continuous event. + Impulsive events are instantaneous changes in the velocity of the body. + Continuous events are an exponentially decaying acceleration acting on the body. + )mydelimiter") + .def(py::init<>()) + .def_readwrite("t", &Event::t, R"mydelimiter( + Time of the event (MJD TDB). + )mydelimiter") + .def_readwrite("bodyName", &Event::bodyName, R"mydelimiter( + Name of the body the event is acting on. Throws an error if the body + is not found in the simulation. + )mydelimiter") + .def_readwrite("isContinuous", &Event::isContinuous, R"mydelimiter( + Whether the event is a continuous event. False by default. + )mydelimiter") + .def_readwrite("eventEst", &Event::eventEst, R"mydelimiter( + Flag for whether to estimate anything for the event. False by default. + )mydelimiter") + .def_readwrite("deltaV", &Event::deltaV, R"mydelimiter( + Delta-V vector of the event. + )mydelimiter") + .def_readwrite("multiplier", &Event::multiplier, R"mydelimiter( + Multiplier of the event. For impulsive events, this is a scalar + multiplier for the delta-V vector. Not used for continuous events. + )mydelimiter") + .def_readwrite("deltaVEst", &Event::deltaVEst, R"mydelimiter( + Flag for whether to estimate the delta-V of the event. False by default. + )mydelimiter") + .def_readwrite("multiplierEst", &Event::multiplierEst, R"mydelimiter( + Flag for whether to estimate the multiplier of the event. False by default. + )mydelimiter") + .def_readwrite("expAccel0", &Event::expAccel0, R"mydelimiter( + Initial acceleration vector for the continuous event. + )mydelimiter") + .def_readwrite("tau", &Event::tau, R"mydelimiter( + Time constant for the exponential decay of the continuous event. + )mydelimiter") + .def_readwrite("expAccel0Est", &Event::expAccel0Est, R"mydelimiter( + Flag for whether to estimate the initial acceleration of the continuous event. False by default. + )mydelimiter") + .def_readwrite("tauEst", &Event::tauEst, R"mydelimiter( + Flag for whether to estimate the time constant of the continuous event. False by default. + )mydelimiter"); + m.def("reconstruct_stm", &reconstruct_stm, py::arg("stm"), R"mydelimiter( Reconstruct the state transition matrix from the flattened vector. @@ -695,26 +741,26 @@ PYBIND11_MODULE(libgrss, m) { Path to the SPICE DE kernel. )mydelimiter") .def_readwrite("consts", &PropSimulation::consts, R"mydelimiter( - Constants of the simulation. PropSimulation.Constants object. + Constants of the simulation. libgrss.Constants object. )mydelimiter") .def_readwrite("integParams", &PropSimulation::integParams, R"mydelimiter( - Integration parameters of the simulation. PropSimulation.IntegParams object. + Integration parameters of the simulation. libgrss.IntegParams object. )mydelimiter") .def_readwrite("spiceBodies", &PropSimulation::spiceBodies, R"mydelimiter( - SPICE bodies of the simulation. List of PropSimulation.SpiceBodies objects. + SPICE bodies of the simulation. List of libgrss.SpiceBodies objects. )mydelimiter") .def_readwrite("integBodies", &PropSimulation::integBodies, R"mydelimiter( - Integration bodies of the simulation. List of PropSimulation.IntegBody objects. + Integration bodies of the simulation. List of libgrss.IntegBody objects. )mydelimiter") .def_readwrite("caParams", &PropSimulation::caParams, R"mydelimiter( - Close approach parameters of the simulation. List of PropSimulation.CloseApproachParameters objects. + Close approach parameters of the simulation. List of libgrss.CloseApproachParameters objects. )mydelimiter") .def_readwrite("impactParams", &PropSimulation::impactParams, R"mydelimiter( - Impact parameters of the simulation. List of PropSimulation.ImpactParameters objects. + Impact parameters of the simulation. List of libgrss.ImpactParameters objects. )mydelimiter") .def_readwrite("t", &PropSimulation::t, R"mydelimiter( Current time of the simulation. @@ -724,7 +770,7 @@ PYBIND11_MODULE(libgrss, m) { )mydelimiter") .def_readwrite("interpParams", &PropSimulation::interpParams, R"mydelimiter( - Interpolation parameters of the simulation. PropSimulation.InterpolationParameters object. + Interpolation parameters of the simulation. libgrss.InterpolationParameters object. )mydelimiter") .def_readwrite("tEvalUTC", &PropSimulation::tEvalUTC, R"mydelimiter( Whether the MJD evaluation time is in UTC for each value in PropSimulation.tEval, @@ -812,7 +858,7 @@ PYBIND11_MODULE(libgrss, m) { Parameters ---------- - body : PropSimulation.SpiceBody + body : libgrss.SpiceBody SPICE body to add to the simulation. )mydelimiter") .def("map_ephemeris", &PropSimulation::map_ephemeris, R"mydelimiter( @@ -848,7 +894,7 @@ PYBIND11_MODULE(libgrss, m) { Parameters ---------- - body : PropSimulation.IntegBody + body : libgrss.IntegBody Integration body to add to the simulation. )mydelimiter") .def("remove_body", &PropSimulation::remove_body, py::arg("name"), @@ -860,47 +906,14 @@ PYBIND11_MODULE(libgrss, m) { name : str Name of the body to remove. )mydelimiter") - .def("add_event", - static_cast, real)>( - &PropSimulation::add_event), - py::arg("body"), py::arg("tEvent"), py::arg("deltaV"), - py::arg("multiplier"), - R"mydelimiter( - Adds an impulsive delta-V event to the simulation. - - Parameters - ---------- - body : str - Name of the body to apply the delta-V to. - tEvent : real - MJD Epoch of the event. Must be in TDB. - deltaV : list of real - Delta-V to apply to the body. - multiplier : real - Multiplier to apply to the delta-V. - )mydelimiter") - .def("add_event", - static_cast, real, real)>( - &PropSimulation::add_event), - py::arg("bodyName"), py::arg("tEvent"), py::arg("expAccel0"), - py::arg("multiplier"), py::arg("tau"), - R"mydelimiter( - Adds an ejecta event to the simulation. + .def("add_event", &PropSimulation::add_event, py::arg("event"), + R"mydelimiter( + Adds an event to the simulation. Parameters ---------- - body : str - Name of the body to apply the delta-V to. - tEvent : real - MJD Epoch of the event. Must be in TDB. - expAccel0 : list of real - Initial exponentiallly decaying acceleration. - multiplier : real - Multiplier for the exponential decay. - tau : real - Time constant for the exponentiallly decaying acceleration. + event : libgrss.Event + Event to add to the simulation. )mydelimiter") .def("set_sim_constants", &PropSimulation::set_sim_constants, py::arg("du2m") = 149597870700.0L, py::arg("tu2s") = 86400.0L, @@ -943,7 +956,7 @@ PYBIND11_MODULE(libgrss, m) { Can be TDB or UTC based on tEvalUTC. tEvalUTC : bool Whether the MJD evaluation time is in UTC for each value in - PropSimulation.tEval, as opposed to TDB. + libgrss.tEval, as opposed to TDB. evalApparentState : bool Whether to evaluate the apparent state of the integration bodies. convergedLightTimes : bool diff --git a/src/simulation.cpp b/src/simulation.cpp index 6d56b838..27962a94 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -352,10 +352,33 @@ void Event::apply_impulsive(PropSimulation *propSim, const real &t, std::vector< for (size_t i = 0; i < 3; i++) { xInteg[velStartIdx + i] += propDir * this->multiplier * this->deltaV[i]; } - if (propSim->integBodies[this->bodyIndex].propStm) { + if (propSim->integBodies[this->bodyIndex].propStm && this->eventEst) { IntegBody *body = &propSim->integBodies[this->bodyIndex]; const size_t numNongravs = body->ngParams.a1Est + body->ngParams.a2Est + body->ngParams.a3Est; - const size_t DeltaVStmIdx = this->xIntegIndex + 6 + 36 + 6*numNongravs; + size_t DeltaVStmIdx = this->xIntegIndex + 6 + 36 + 6*numNongravs; + for (size_t i = 0; i < propSim->eventMngr.impulsiveEvents.size(); i++) { + const bool sameBody = propSim->eventMngr.impulsiveEvents[i].bodyIndex == this->bodyIndex; + const bool hasStarted = propSim->eventMngr.impulsiveEvents[i].hasStarted; + const bool deltaVEst = propSim->eventMngr.impulsiveEvents[i].deltaVEst; + const bool multiplierEst = propSim->eventMngr.impulsiveEvents[i].multiplierEst; + if (sameBody && hasStarted && multiplierEst) { + DeltaVStmIdx += 6*1; // only multiplier + } else if (sameBody && hasStarted && deltaVEst) { + DeltaVStmIdx += 6*3; // full deltaV vector + } + } + for (size_t i = 0; i < propSim->eventMngr.continuousEvents.size(); i++) { + const bool sameBody = propSim->eventMngr.continuousEvents[i].bodyIndex == this->bodyIndex; + const bool hasStarted = propSim->eventMngr.continuousEvents[i].hasStarted; + const bool expAccel0Est = propSim->eventMngr.continuousEvents[i].expAccel0Est; + const bool tauEst = propSim->eventMngr.continuousEvents[i].tauEst; + if (sameBody && hasStarted && expAccel0Est) { + DeltaVStmIdx += 6*3; // expAccel0 vector + } + if (sameBody && hasStarted && tauEst) { + DeltaVStmIdx += 6*1; // time constant + } + } if (this->multiplierEst){ xInteg[DeltaVStmIdx+3] = propDir * this->deltaV[0]; xInteg[DeltaVStmIdx+4] = propDir * this->deltaV[1]; @@ -366,6 +389,7 @@ void Event::apply_impulsive(PropSimulation *propSim, const real &t, std::vector< xInteg[DeltaVStmIdx+17] = propDir; } } + this->hasStarted = true; } /** @@ -1001,7 +1025,7 @@ void PropSimulation::remove_body(std::string name) { std::cout << "Error: Body " << name << " not found." << std::endl; } -static size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, +static size_t event_preprocess(PropSimulation *propSim, const std::string &eventBodyName, const real &tEvent) { // check if tEvent is valid const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0; @@ -1017,14 +1041,14 @@ static size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, bool bodyExists = false; size_t bodyIndex; for (size_t i = 0; i < propSim->integParams.nInteg; i++) { - if (propSim->integBodies[i].name == body.name) { + if (propSim->integBodies[i].name == eventBodyName) { bodyExists = true; bodyIndex = i; break; } } if (!bodyExists) { - throw std::invalid_argument("Integration body with name " + body.name + + throw std::invalid_argument("Integration body with name " + eventBodyName + " does not exist in simulation " + propSim->name); } @@ -1032,11 +1056,27 @@ static size_t event_preprocess(PropSimulation *propSim, const IntegBody &body, } static void event_stm_handling(PropSimulation *propSim, IntegBody *body, const Event &event){ - int numEventParams = 3; - if (event.multiplierEst) { + int numEventParams = -1; + if (!event.isContinuous && !event.deltaVEst && event.multiplierEst) { + // only multiplier is estimated for an impulsive delta-V + numEventParams = 1; + } else if (!event.isContinuous && event.deltaVEst && !event.multiplierEst) { + // only delta-V vector is estimated for an impulsive delta-V + numEventParams = 3; + } else if (event.isContinuous && !event.expAccel0Est && event.tauEst) { + // only time constant is estimated for a continuous event numEventParams = 1; + } else if (event.isContinuous && event.expAccel0Est && !event.tauEst) { + // only exponential initial acceleration is estimated for a continuous event + numEventParams = 3; + } else if (event.isContinuous && event.expAccel0Est && event.tauEst) { + // exponential initial acceleration and time constant are estimated for a continuous event + numEventParams = 4; + } else { + throw std::invalid_argument( + "event_stm_handling: Invalid event estimation configuration for analytic STM " + "propagation."); } - body->n2Derivs += 3*numEventParams; propSim->integBodies[event.bodyIndex].n2Derivs += 3*numEventParams; propSim->integParams.n2Derivs += 3*numEventParams; @@ -1044,12 +1084,15 @@ static void event_stm_handling(PropSimulation *propSim, IntegBody *body, const E std::vector extraVec = std::vector(6*numEventParams, 0.0); // add extra vector at the end of stm vector for (size_t i = 0; i < extraVec.size(); i++) { - body->stm.push_back(extraVec[i]); propSim->integBodies[event.bodyIndex].stm.push_back(extraVec[i]); } } static void event_postprocess(PropSimulation *propSim, Event &event) { + // modify stm if estimating Delta-V + if (propSim->integBodies[event.bodyIndex].propStm && event.eventEst) { + event_stm_handling(propSim, &propSim->integBodies[event.bodyIndex], event); + } // assign event xIntegIndex event.xIntegIndex = 0; for (size_t i = 0; i < event.bodyIndex; i++) { @@ -1078,67 +1121,41 @@ static void event_postprocess(PropSimulation *propSim, Event &event) { } /** - * @param[in] body IntegBody object to apply the ImpulseEvent to. - * @param[in] tEvent Time at which to apply the ImpulseEvent. - * @param[in] deltaV Delta-V for the impulse. - * @param[in] multiplier Multiplier for the Delta-V. + * @param[in] event Event object to add to the simulation. */ -void PropSimulation::add_event(IntegBody body, real tEvent, - std::vector deltaV, real multiplier) { - size_t bodyIndex = event_preprocess(this, body, tEvent); - Event event; - event.t = tEvent; - event.bodyName = body.name; +void PropSimulation::add_event(Event event) { + size_t bodyIndex = event_preprocess(this, event.bodyName, event.t); event.bodyIndex = bodyIndex; - event.deltaV = deltaV; - // if multiplier is numeric, estimate it otherwise estimate full deltaV vector - if (!std::isfinite(multiplier)) { - if (body.propStm){ - event.multiplierEst = false; + event.hasStarted = false; + if (event.isContinuous) { + // make sure exponential accelerations are not nan and time constant is postive + if (!std::isfinite(event.expAccel0[0]) || !std::isfinite(event.expAccel0[1]) || + !std::isfinite(event.expAccel0[2]) || !std::isfinite(event.tau) || event.tau <= 0.0) { + throw std::invalid_argument( + "add_event: Exponential acceleration must be finite and time " + "constant must be finite and positive."); + } + // make sure deltaV is nan and multiplier is nan + if (std::isfinite(event.deltaV[0]) || std::isfinite(event.deltaV[1]) || + std::isfinite(event.deltaV[2]) || std::isfinite(event.multiplier)) { + throw std::invalid_argument("add_event: Delta-V and multiplier must be nan."); } - multiplier = 1.0; + this->eventMngr.nConEvents++; } else { - if (body.propStm) { - event.multiplierEst = true; + // make sure deltaV is not nan and multiplier is not nan + if (!std::isfinite(event.deltaV[0]) || !std::isfinite(event.deltaV[1]) || + !std::isfinite(event.deltaV[2]) || !std::isfinite(event.multiplier)) { + throw std::invalid_argument("add_event: Delta-V and multiplier must be finite."); } + // make sure exponential acceleration is nan + if (std::isfinite(event.expAccel0[0]) || std::isfinite(event.expAccel0[1]) || + std::isfinite(event.expAccel0[2]) || std::isfinite(event.tau)) { + throw std::invalid_argument( + "add_event: Exponential acceleration must be zero."); + } + this->eventMngr.nImpEvents++; } - event.multiplier = multiplier; - // modify stm if estimating Delta-V - if (body.propStm) { - event_stm_handling(this, &body, event); - } - - event.tau = 0.0; - event.isContinuous = false; - event.isHappening = false; - event_postprocess(this, event); - this->eventMngr.nImpEvents++; -} - -/** - * @param[in] body IntegBody object to apply the EjectaEvent to. - * @param[in] tEvent Time at which to apply the EjectaEvent. - * @param[in] expAccel0 Initial exponentiallly decaying acceleration. - * @param[in] multiplier Multiplier for the exponential decay. - * @param[in] tau Time constant for the exponentiallly decaying acceleration. - */ -void PropSimulation::add_event(IntegBody body, real tEvent, - std::vector expAccel0, real multiplier, - real tau) { - size_t bodyIndex = event_preprocess(this, body, tEvent); - Event event; - event.t = tEvent; - event.bodyName = body.name; - event.bodyIndex = bodyIndex; - event.deltaV = std::vector(3, 0.0); - event.multiplier = multiplier; - - event.expAccel0 = expAccel0; - event.tau = tau; - event.isContinuous = true; - event.isHappening = false; event_postprocess(this, event); - this->eventMngr.nConEvents++; } /** diff --git a/src/stm.cpp b/src/stm.cpp index 664e7c7f..99bc60ee 100644 --- a/src/stm.cpp +++ b/src/stm.cpp @@ -525,3 +525,54 @@ void stm_nongrav(STMParameters &stmParams, const real &g, stmParams.dfdpar[8] += g*nHat[2]; } } + +/** + * @param[inout] stmParams Structure containing the STM submatrices and their derivatives. + * @param[in] propSim PropSimulation object. + * @param[in] eventIdx Index of the continuous event. + * @param[in] tPastEvent Time since the continuous event started [days]. + * @param[in] postFac Exponential decay factor for the continuous event. + */ +void stm_continuous_event(STMParameters &stmParams, + const PropSimulation *propSim, const size_t &eventIdx, + const real &tPastEvent, const real &postFac){ + const Event *thisEvent = &propSim->eventMngr.continuousEvents[eventIdx]; + const IntegBody *body = &propSim->integBodies[thisEvent->bodyIndex]; + const size_t numNongravs = body->ngParams.a1Est + body->ngParams.a2Est + body->ngParams.a3Est; + size_t dfdparIdx = 3*numNongravs; + for (size_t i = 0; i < propSim->eventMngr.impulsiveEvents.size(); i++) { + const bool sameBody = propSim->eventMngr.impulsiveEvents[i].bodyIndex == thisEvent->bodyIndex; + const bool hasStarted = propSim->eventMngr.impulsiveEvents[i].hasStarted; + const bool deltaVEst = propSim->eventMngr.impulsiveEvents[i].deltaVEst; + const bool multiplierEst = propSim->eventMngr.impulsiveEvents[i].multiplierEst; + if (sameBody && hasStarted && multiplierEst) { + dfdparIdx += 3*1; // only multiplier + } else if (sameBody && hasStarted && deltaVEst) { + dfdparIdx += 3*3; // full deltaV vector + } + } + for (size_t i = 0; i < eventIdx; i++) { + const bool sameBody = propSim->eventMngr.continuousEvents[i].bodyIndex == thisEvent->bodyIndex; + const bool hasStarted = propSim->eventMngr.continuousEvents[i].hasStarted; + const bool expAccel0Est = propSim->eventMngr.continuousEvents[i].expAccel0Est; + const bool tauEst = propSim->eventMngr.continuousEvents[i].tauEst; + if (sameBody && hasStarted && expAccel0Est) { + dfdparIdx += 3*3; // expAccel0 vector + } + if (sameBody && hasStarted && tauEst) { + dfdparIdx += 3*1; // time constant + } + } + if (thisEvent->isContinuous && thisEvent->expAccel0Est) { + stmParams.dfdpar[dfdparIdx + 0] = postFac; + stmParams.dfdpar[dfdparIdx + 4] = postFac; + stmParams.dfdpar[dfdparIdx + 8] = postFac; + dfdparIdx += 9; + } + if (thisEvent->isContinuous && thisEvent->tauEst) { + const real postFac2 = tPastEvent/thisEvent->tau/thisEvent->tau; + for (size_t i = 0; i < 3; i++) { + stmParams.dfdpar[dfdparIdx + i] = -thisEvent->expAccel0[i]*postFac*postFac2; + } + } +}