From cd95c112d3c66035cb4f79dda64bda43a200127d Mon Sep 17 00:00:00 2001 From: Travis Hunter Date: Sun, 23 Apr 2023 16:43:55 -0600 Subject: [PATCH] adding mass calculations to MotorModel --- gui/MainWindow.cpp | 4 +- gui/ThrustCurveMotorSelector.cpp | 4 +- model/MotorModel.cpp | 84 +++++++++++++++++++++++++++++ model/MotorModel.h | 78 ++++++++++++++++++--------- model/ThrustCurve.cpp | 9 +++- model/ThrustCurve.h | 2 + qtrocket.pro | 2 + sim/Propagator.cpp | 6 +-- utils/RSEDatabaseLoader.cpp | 6 ++- utils/ThrustCurveAPI.cpp | 6 ++- utils/math/UtilityMathFunctions.cpp | 15 ++++++ utils/math/UtilityMathFunctions.h | 40 ++++++++++++++ 12 files changed, 220 insertions(+), 36 deletions(-) create mode 100644 utils/math/UtilityMathFunctions.cpp create mode 100644 utils/math/UtilityMathFunctions.h diff --git a/gui/MainWindow.cpp b/gui/MainWindow.cpp index 68245d1..4b7ec5a 100644 --- a/gui/MainWindow.cpp +++ b/gui/MainWindow.cpp @@ -146,8 +146,8 @@ void MainWindow::on_loadRSE_button_clicked() const std::vector& motors = loader.getMotors(); for(const auto& motor : motors) { - std::cout << "Adding: " << motor.commonName << std::endl; - engineSelector->addItem(QString(motor.commonName.c_str())); + std::cout << "Adding: " << motor.data.commonName << std::endl; + engineSelector->addItem(QString(motor.data.commonName.c_str())); } } diff --git a/gui/ThrustCurveMotorSelector.cpp b/gui/ThrustCurveMotorSelector.cpp index aee1681..d19228b 100644 --- a/gui/ThrustCurveMotorSelector.cpp +++ b/gui/ThrustCurveMotorSelector.cpp @@ -69,7 +69,7 @@ void ThrustCurveMotorSelector::on_searchButton_clicked() for(const auto& i : motors) { - ui->motorSelection->addItem(QString::fromStdString(i.commonName)); + ui->motorSelection->addItem(QString::fromStdString(i.data.commonName)); } } @@ -87,7 +87,7 @@ void ThrustCurveMotorSelector::on_setMotor_clicked() std::end(motorModels), [&commonName](const auto& item) { - return item.commonName == commonName; + return item.data.commonName == commonName; }); QtRocket::getInstance()->getRocket()->setMotorModel(mm); diff --git a/model/MotorModel.cpp b/model/MotorModel.cpp index 8e79192..985d71e 100644 --- a/model/MotorModel.cpp +++ b/model/MotorModel.cpp @@ -1,4 +1,5 @@ #include "model/MotorModel.h" +#include "utils/math/Constants.h" namespace model { @@ -13,4 +14,87 @@ MotorModel::~MotorModel() } +double MotorModel::getMass(double simTime) const +{ + // the current mass is the emptyMass + the current prop mass + + // If ignition hasn't occurred, return the totalMass + if(!ignitionOccurred) + { + return data.totalWeight; + } + else if(simTime - ignitionTime <= data.burnTime) + { + double thrustTime = simTime - ignitionTime; + // Find the right interval in the massCurve + auto i = massCurve.cbegin(); + while(i->first <= thrustTime) + { + // If thrustTime is equal to a data point that we have, then just return + // the mass at that time. Otherwise it fell between two points and we + // will interpolate + if(i->first == thrustTime) + { + // return empty mass + remaining propellant mass + return emptyMass + i->second; + } + else + { + i++; + } + } + // linearly interpolate the propellant mass. Then return the empty mass + remaining prop mass + double tStart = std::prev(i)->first; + double tEnd = i->first; + double propMassStart = std::prev(i)->second; + double propMassEnd = i->second; + double slope = (propMassEnd - propMassStart) / (tEnd - tStart); + double currentMass = emptyMass + propMassStart + (thrustTime - tStart) * slope; + return currentMass; + + } + else + { + return emptyMass; + } +} + +void MotorModel::setMetaData(const MetaData& md) +{ + data = md; + emptyMass = data.totalWeight - data.propWeight; + + // Calculate the Isp for the motor, as we'll need this for the computing the mass flow rate. + // This will be the total impulse in Newton-seconds over + // the propellant weight. The prop mass is in grams, hence the division by 1000.0 to get kg + isp = data.totalImpulse / (utils::math::Constants::g0 * data.propWeight / 1000.0); + + // Precompute the mass curve. Having this precomputed will ensure multiple calls to getMass() + // or getThrust() during the same time step don't accidentally decrement the mass multiple times. + // Having a lookup table will ensure consistent mass values, as well as speed up the simulation, + // just at the cost of some extra space + + // Most motor data in the RASP format has a limitation of 32 data points. We're not going to + // match that, so we can pick whatever we want and just interpolate values. We can have 128 + // for example + + massCurve.reserve(128); + double timeStep = data.burnTime / 127.0; + double t = 0.0; + double propMass{data.propWeight}; + for(std::size_t i = 0; i < 127; ++i) + { + massCurve.push_back(std::make_pair(t + i*timeStep, propMass)); + propMass -= thrust.getThrust(t + i*timeStep) * timeStep * data.propWeight / data.totalImpulse; + } + massCurve.push_back(std::make_pair(data.burnTime, 0.0)); + + + +} + +void MotorModel::moveMetaData(MetaData&& md) +{ + data = std::move(md); +} } // namespace model diff --git a/model/MotorModel.h b/model/MotorModel.h index b5df4cb..17524c7 100644 --- a/model/MotorModel.h +++ b/model/MotorModel.h @@ -255,11 +255,12 @@ public: if(name == "SU" || name == "Single Use") return MOTORTYPE::SU; - else if(name == "Hybrid") - return MOTORTYPE::HYBRID; else if(name == "reload" || name == "Reload") return MOTORTYPE::RELOAD; + else // It's a hybrid + return MOTORTYPE::HYBRID; + } }; @@ -349,35 +350,64 @@ public: } }; + /// TODO: make these MotorModel members private. Public just for testing //private: - MotorAvailability availability{AVAILABILITY::REGULAR}; /// Motor Availability - double avgThrust{0.0}; /// Average thrust in Newtons - double burnTime{0.0}; /// Burn time in seconds - CertOrg certOrg{CERTORG::UNC}; /// The certification organization, defaults to Uncertified - std::string commonName{""}; /// Common name, e.g. A8 or J615 - // int dataFiles - std::vector delays; /// 1000 delay means no ejection charge - std::string designation{""}; /// Other name, usually includes prop type, e.g. H110W - double diameter{0}; /// motor diameter in mm - std::string impulseClass; /// Motor letter, e.g. 'A', 'B', '1/2A', 'M', etc - std::string infoUrl{""}; /// TODO: ??? - double length{0.0}; /// motor length in mm - MotorManufacturer manufacturer{MOTORMANUFACTURER::UNKNOWN}; /// Motor Manufacturer + struct MetaData + { + MetaData() = default; + ~MetaData() = default; + MetaData(const MetaData&) = default; + MetaData(MetaData&&) = default; + MetaData& operator=(const MetaData&) = default; + MetaData& operator=(MetaData&&) = default; + + MotorAvailability availability{AVAILABILITY::REGULAR}; /// Motor Availability + double avgThrust{0.0}; /// Average thrust in Newtons + double burnTime{0.0}; /// Burn time in seconds + CertOrg certOrg{CERTORG::UNC}; /// The certification organization, defaults to Uncertified + std::string commonName{""}; /// Common name, e.g. A8 or J615 + // int dataFiles + std::vector delays; /// 1000 delay means no ejection charge + std::string designation{""}; /// Other name, usually includes prop type, e.g. H110W + double diameter{0}; /// motor diameter in mm + std::string impulseClass; /// Motor letter, e.g. 'A', 'B', '1/2A', 'M', etc + std::string infoUrl{""}; /// TODO: ??? + double length{0.0}; /// motor length in mm + MotorManufacturer manufacturer{MOTORMANUFACTURER::UNKNOWN}; /// Motor Manufacturer + + double maxThrust{0.0}; /// Max thrust in Newtons + std::string motorIdTC{""}; /// 24 character hex string used by thrustcurve.org to ID a motor + std::string propType{""}; /// Propellant type, e.g. black powder + double propWeight{0.0}; /// Propellant weight in grams + bool sparky{false}; /// true if the motor is "sparky", false otherwise + double totalImpulse{0.0}; /// Total impulse in Newton-seconds + double totalWeight{0.0}; /// Total weight in grams + MotorType type{MOTORTYPE::SU}; /// Motor type, e.g. single-use, reload, or hybrid + std::string lastUpdated{""}; /// Date last updated on ThrustCurve.org + }; + + double getMass(double simTime) const; + double getThrust(double simTime); + + void setMetaData(const MetaData& md); + void moveMetaData(MetaData&& md); + + void startMotor(double startTime) { ignitionOccurred = true; ignitionTime = startTime; } - double maxThrust{0.0}; /// Max thrust in Newtons - std::string motorIdTC{""}; /// 24 character hex string used by thrustcurve.org to ID a motor - std::string propType{""}; /// Propellant type, e.g. black powder - double propWeight{0.0}; /// Propellant weight in grams - bool sparky{false}; /// true if the motor is "sparky", false otherwise - double totalImpulse{0.0}; /// Total impulse in Newton-seconds - double totalWeight{0.0}; /// Total weight in grams - MotorType type{MOTORTYPE::SU}; /// Motor type, e.g. single-use, reload, or hybrid - std::string lastUpdated{""}; /// Date last updated on ThrustCurve.org // Thrust parameters + MetaData data; +private: + bool ignitionOccurred{false}; + double emptyMass; + double isp; + double maxTime; + double ignitionTime; ThrustCurve thrust; /// The measured motor thrust curve + + std::vector> massCurve; }; diff --git a/model/ThrustCurve.cpp b/model/ThrustCurve.cpp index a5c4539..a19cb5f 100644 --- a/model/ThrustCurve.cpp +++ b/model/ThrustCurve.cpp @@ -7,7 +7,8 @@ // 3rd party headers /// \endcond -#include "ThrustCurve.h" +#include "model/ThrustCurve.h" +#include "utils/Logger.h" ThrustCurve::ThrustCurve(std::vector>& tc) @@ -54,6 +55,12 @@ void ThrustCurve::setIgnitionTime(double t) double ThrustCurve::getThrust(double t) { // calculate t relative to the start time of the motor + static bool burnout{false}; + if(!burnout && t >= (maxTime + ignitionTime)) + { + burnout = true; + utils::Logger::getInstance()->info("Motor burnout at time: " + std::to_string(t)); + } t -= ignitionTime; if(t < 0.0 || t > maxTime) { diff --git a/model/ThrustCurve.h b/model/ThrustCurve.h index 9bf6129..c50c064 100644 --- a/model/ThrustCurve.h +++ b/model/ThrustCurve.h @@ -41,6 +41,8 @@ public: void setIgnitionTime(double t); + double getMaxTime() const { return maxTime; } + /** * TODO: Get rid of this. This is for temporary testing */ diff --git a/qtrocket.pro b/qtrocket.pro index 02a5458..c43ceff 100644 --- a/qtrocket.pro +++ b/qtrocket.pro @@ -37,6 +37,7 @@ SOURCES += \ utils/ThreadPool.cpp \ utils/ThrustCurveAPI.cpp \ utils/math/Quaternion.cpp \ + utils/math/UtilityMathFunctions.cpp \ utils/math/Vector3.cpp HEADERS += \ @@ -75,6 +76,7 @@ HEADERS += \ utils/Triplet.h \ utils/math/Constants.h \ utils/math/Quaternion.h \ + utils/math/UtilityMathFunctions.h \ utils/math/Vector3.h FORMS += \ diff --git a/sim/Propagator.cpp b/sim/Propagator.cpp index c204f6a..d812ee8 100644 --- a/sim/Propagator.cpp +++ b/sim/Propagator.cpp @@ -47,9 +47,9 @@ Propagator::Propagator(Rocket* r) /* dpitch/dt */ [](double, const std::vector& s) -> double { return s[9]; }, /* dyaw/dt */ [](double, const std::vector& s) -> double { return s[10]; }, /* droll/dt */ [](double, const std::vector& s) -> double { return s[11]; }, - /* dpitchRate/dt */ [this](double, const std::vector& s) -> double { (getTorqueP() - s[7] * s[8] * (getIroll() - getIyaw())) / getIpitch(); }, - /* dyawRate/dt */ [this](double, const std::vector& s) -> double { (getTorqueQ() - s[6] * s[9] * (getIpitch() - getIroll())) / getIyaw(); }, - /* drollRate/dt */ [this](double, const std::vector& s) -> double { (getTorqueR() - s[6] * s[7] * (getIyaw() - getIpitch())) / getIroll(); })); + /* dpitchRate/dt */ [this](double, const std::vector& s) -> double { return (getTorqueP() - s[7] * s[8] * (getIroll() - getIyaw())) / getIpitch(); }, + /* dyawRate/dt */ [this](double, const std::vector& s) -> double { return (getTorqueQ() - s[6] * s[9] * (getIpitch() - getIroll())) / getIyaw(); }, + /* drollRate/dt */ [this](double, const std::vector& s) -> double { return (getTorqueR() - s[6] * s[7] * (getIyaw() - getIpitch())) / getIroll(); })); integrator->setTimeStep(timeStep); diff --git a/utils/RSEDatabaseLoader.cpp b/utils/RSEDatabaseLoader.cpp index 47851fc..cc4f766 100644 --- a/utils/RSEDatabaseLoader.cpp +++ b/utils/RSEDatabaseLoader.cpp @@ -35,7 +35,7 @@ RSEDatabaseLoader::~RSEDatabaseLoader() void RSEDatabaseLoader::buildAndAppendMotorModel(boost::property_tree::ptree& v) { - model::MotorModel mm; + model::MotorModel::MetaData mm; mm.availability = model::MotorModel::MotorAvailability(model::MotorModel::AVAILABILITY::REGULAR); mm.avgThrust = v.get(".avgThrust", 0.0); mm.burnTime = v.get(".burn-time", 0.0); @@ -70,7 +70,9 @@ void RSEDatabaseLoader::buildAndAppendMotorModel(boost::property_tree::ptree& v) thrustData.push_back(std::make_pair(tdata, fdata)); } - motors.emplace_back(std::move(mm)); + model::MotorModel motorModel; + motorModel.moveMetaData(std::move(mm)); + motors.emplace_back(std::move(motorModel)); } diff --git a/utils/ThrustCurveAPI.cpp b/utils/ThrustCurveAPI.cpp index 8778d92..5039c6e 100644 --- a/utils/ThrustCurveAPI.cpp +++ b/utils/ThrustCurveAPI.cpp @@ -148,7 +148,8 @@ std::vector ThrustCurveAPI::searchMotors(const SearchCriteria iter != jsonResult["results"].end(); ++iter) { - model::MotorModel mm; + model::MotorModel motorModel; + model::MotorModel::MetaData mm; mm.commonName = (*iter)["commonName"].asString(); std::string availability = (*iter)["availability"].asString(); @@ -185,7 +186,8 @@ std::vector ThrustCurveAPI::searchMotors(const SearchCriteria else mm.type = model::MotorModel::MotorType(model::MotorModel::MOTORTYPE::HYBRID); - retVal.push_back(mm); + motorModel.moveMetaData(std::move(mm)); + retVal.push_back(motorModel); } } catch(const std::exception& e) diff --git a/utils/math/UtilityMathFunctions.cpp b/utils/math/UtilityMathFunctions.cpp new file mode 100644 index 0000000..21d5dfa --- /dev/null +++ b/utils/math/UtilityMathFunctions.cpp @@ -0,0 +1,15 @@ + +/// \cond +// C headers +// C++ headers +// 3rd party headers +/// \endcond + + +namespace utils +{ +namespace math +{ + +} // namespace math +} // namespace utils diff --git a/utils/math/UtilityMathFunctions.h b/utils/math/UtilityMathFunctions.h new file mode 100644 index 0000000..bf47777 --- /dev/null +++ b/utils/math/UtilityMathFunctions.h @@ -0,0 +1,40 @@ +#ifndef UTILITYMATHFUNCTIONS_H +#define UTILITYMATHFUNCTIONS_H + +/// \cond +// C headers +// C++ headers +#include +#include +// 3rd party headers +/// \endcond + +namespace utils +{ +namespace math +{ + +/** + * @brief doublesEqual compares two doubles for equality, using a value for number of smallest + * places to ignore + * + * This is derived from the example on cppreference.com: https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon + * The default ulp of 8 with a numeric_limits::epsilon() of ~2e-16, so ulp of 8 yields 12 + * significant figures. A ulp of 10 yields 6 significant figures. + * @param a the first double to compare + * @param b the second double to compare + * @param ulp number of smallest decimal places to ignore + * @return + */ + +template +bool floatingPointEqual(T a, T b, int ulp = 4) +{ + return std::fabs(a - b) <= std::numeric_limits::epsilon() * std::fabs(a + b) * ulp + // unless the result is subnormal + || std::fabs(a - b) < std::numeric_limits::min(); +} + +} // namespace math +} // namespace utils +#endif // UTILITYMATHFUNCTIONS_H