adding mass calculations to MotorModel

gui/MainWindow.cpp

@@ -146,8 +146,8 @@ void MainWindow::on_loadRSE_button_clicked()
    const std::vector<model::MotorModel>& 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()));
    }
 }
 

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);

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

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<int> 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<int> 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<std::pair<double, double>> massCurve;
    
 };
 

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<std::pair<double, double>>& 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)
    {

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
     */

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 += \

sim/Propagator.cpp

@@ -47,9 +47,9 @@ Propagator::Propagator(Rocket* r)
        /* dpitch/dt    */ [](double, const std::vector<double>& s) -> double { return s[9]; },
        /* dyaw/dt      */ [](double, const std::vector<double>& s) -> double { return s[10]; },
        /* droll/dt     */ [](double, const std::vector<double>& s) -> double { return s[11]; },
-       /* dpitchRate/dt */ [this](double, const std::vector<double>& s) -> double { (getTorqueP() - s[7] * s[8] * (getIroll() - getIyaw())) / getIpitch(); },
-       /* dyawRate/dt   */ [this](double, const std::vector<double>& s) -> double { (getTorqueQ() - s[6] * s[9] * (getIpitch() - getIroll())) / getIyaw(); },
-      /* drollRate/dt   */ [this](double, const std::vector<double>& s) -> double { (getTorqueR() - s[6] * s[7] * (getIyaw() - getIpitch())) / getIroll(); }));
+       /* dpitchRate/dt */ [this](double, const std::vector<double>& s) -> double { return (getTorqueP() - s[7] * s[8] * (getIroll() - getIyaw())) / getIpitch(); },
+       /* dyawRate/dt   */ [this](double, const std::vector<double>& s) -> double { return (getTorqueQ() - s[6] * s[9] * (getIpitch() - getIroll())) / getIyaw(); },
+      /* drollRate/dt   */ [this](double, const std::vector<double>& s) -> double { return (getTorqueR() - s[6] * s[7] * (getIyaw() - getIpitch())) / getIroll(); }));
 
 
    integrator->setTimeStep(timeStep);

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<double>("<xmlattr>.avgThrust", 0.0);
    mm.burnTime  = v.get<double>("<xmlattr>.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));
 }
 
 

utils/ThrustCurveAPI.cpp

@@ -148,7 +148,8 @@ std::vector<model::MotorModel> 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<model::MotorModel> 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)

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

utils/math/UtilityMathFunctions.h

@@ -0,0 +1,40 @@
+#ifndef UTILITYMATHFUNCTIONS_H
+#define UTILITYMATHFUNCTIONS_H
+
+/// \cond
+// C headers
+// C++ headers
+#include <cmath>
+#include <limits>
+// 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<double>::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<typename T>
+bool floatingPointEqual(T a, T b, int ulp = 4)
+{
+   return std::fabs(a - b) <= std::numeric_limits<T>::epsilon() * std::fabs(a + b) * ulp
+       // unless the result is subnormal
+       || std::fabs(a - b) < std::numeric_limits<T>::min();
+}
+
+} // namespace math
+} // namespace utils
+#endif // UTILITYMATHFUNCTIONS_H