+ enum prottypes, + rnb interpolation

utils/PropInfoBase.cpp

@@ -9,7 +9,7 @@ PropInfoBase::PropInfoBase(
     const bool sun, 
     const bool moon, 
     const IntegratorCfg& icfg, 
-    const int type)
+    const proptypes type)
     :
     m_atmModel(std::move(atmModel)),
     m_earthModel(std::move(earthModel)),
@@ -29,11 +29,11 @@ bool PropInfoBase::Load(myXML::Contain* parent)
     std::string type;
     LoadXML(con, "Propagator", type);
 
-    if ((m_type = TypeConvert(type)) == -1) {
-        throw (Exp() << "Unsupported model type:" << m_type);
+    if ((m_type = TypeConvert(type)) == proptypes::UNKNOWN) {
+        throw Exp() << "Unsupported model type:" << TypeConvert(m_type);
     }
 
-    if (m_type == SGP_PROP || m_type == PROGNOZ_PROP || m_type == KEPLER_PROP) {
+    if (m_type == proptypes::SGP_PROP || m_type == proptypes::PROGNOZ_PROP || m_type == proptypes::KEPLER_PROP) {
         return true;
     }
 
@@ -54,11 +54,11 @@ void PropInfoBase::load(const jsonval& parent)
     std::string type;
     loadjson(parent, "prop", type);
 
-    if ((m_type = TypeConvert(type)) == -1) {
-        throw (Exp() << "Unsupported model type:" << m_type);
+    if ((m_type = TypeConvert(type)) == proptypes::UNKNOWN) {
+        throw Exp() << "Unsupported model type:" << TypeConvert(m_type);
     }
 
-    if (m_type == SGP_PROP || m_type == PROGNOZ_PROP || m_type == KEPLER_PROP) {
+    if (m_type == proptypes::SGP_PROP || m_type == proptypes::PROGNOZ_PROP || m_type == proptypes::KEPLER_PROP) {
         return;
     }
 
@@ -77,7 +77,7 @@ std::string PropInfoBase::GetAsXML() const
 
     PutValue(outs, "Propagator", TypeConvert(m_type));
 
-    if (m_type == SGP_PROP || m_type == PROGNOZ_PROP || m_type == KEPLER_PROP) {
+    if (m_type == proptypes::SGP_PROP || m_type == proptypes::PROGNOZ_PROP || m_type == proptypes::KEPLER_PROP) {
         return outs.str();
     }
     else {
@@ -106,18 +106,18 @@ std::string PropInfoBase::GetAsXML() const
     return outs.str();
 }
 
-std::string PropInfoBase::TypeConvert(int value)
+std::string PropInfoBase::TypeConvert(proptypes value)
 {
-    if (value == ST_PROP) {
+    if (value == proptypes::ST_PROP) {
         return "ST";
     }
-    else if (value == KEPLER_PROP) {
+    else if (value == proptypes::KEPLER_PROP) {
         return "KEPLER";
     }
-    else if (value == PROGNOZ_PROP) {
+    else if (value == proptypes::PROGNOZ_PROP) {
         return "PROGNOZ";
     }
-    else if (value == SGP_PROP) {
+    else if (value == proptypes::SGP_PROP) {
         return "SGP";
     }
     else {
@@ -125,22 +125,22 @@ std::string PropInfoBase::TypeConvert(int value)
     }
 }
 
-int PropInfoBase::TypeConvert(const std::string& value)
+proptypes PropInfoBase::TypeConvert(const std::string& value)
 {
     if (value == "ST") {
-        return ST_PROP;
+        return proptypes::ST_PROP;
     }
     else if (value == "KEPLER") {
-        return KEPLER_PROP;
+        return proptypes::KEPLER_PROP;
     }
     else if (value == "PROGNOZ") {
-        return PROGNOZ_PROP;
+        return proptypes::PROGNOZ_PROP;
     }
     else if (value == "SGP") {
-        return SGP_PROP;
+        return proptypes::SGP_PROP;
     }
     else {
-        return -1;
+        return proptypes::UNKNOWN;
     }
 }
 

utils/PropInfoBase.h

@@ -8,12 +8,13 @@
 #include "json_functions.h"
 #include "integratorcfg.h"
 
-enum proptypes
+enum class proptypes
 {
     ST_PROP = 0,
     KEPLER_PROP,
     PROGNOZ_PROP,
-    SGP_PROP
+    SGP_PROP,
+    UNKNOWN
 };
 
 class PropInfoBase
@@ -29,13 +30,13 @@ public:
         const bool sun,
         const bool moon,
         const IntegratorCfg& icfg,
-        const int type);
+        const proptypes type);
 
     bool Load(myXML::Contain* parent);
     void load(const jsonval& parent);
     bool IsSun() const { return m_sun; }
     bool IsMoon() const { return m_moon; }
-    int GetType() const { return m_type; }
+    proptypes GetType() const { return m_type; }
     //bool is_prec_itrf() const { return m_prec_itrf; }
     const AtmModel& GetAtmModel() const { return m_atmModel; }
     const EarthModel& GetEarthModel() const { return m_earthModel; }
@@ -45,15 +46,15 @@ public:
     const IntegratorCfg& get_integrator_cfg() const { return m_integrator_cfg; }
 private:
 
-    static std::string TypeConvert(int value) ;
-    static int TypeConvert(const std::string& value);
+    static std::string TypeConvert(proptypes value) ;
+    static proptypes TypeConvert(const std::string& value);
 
     AtmModel m_atmModel;
     EarthModel m_earthModel;
     std::string m_lightModel;
     bool m_sun;
     bool m_moon;
-    int m_type;
+    proptypes m_type;
     int m_id;
     IntegratorCfg m_integrator_cfg;
 };

utils/Propagator.cpp

@@ -40,7 +40,7 @@ bool Propagator::Init(const PropInfoBase& propInfo, const SInitOrbit& initOrbit,
 {
     m_variates = variates;
 
-    if (propInfo.GetType() == ST_PROP) {
+    if (propInfo.GetType() == proptypes::ST_PROP) {
         if (!m_variates) {
             InitStPredictor(m_stPredictor, initOrbit, propInfo);
             m_stPredictor.GetPredictor()->SetAccuracy(
@@ -60,7 +60,7 @@ bool Propagator::Init(const PropInfoBase& propInfo, const SInitOrbit& initOrbit,
             m_predictor = &m_stvPredictor;
         }
     }
-    else if (propInfo.GetType() == KEPLER_PROP) {
+    else if (propInfo.GetType() == proptypes::KEPLER_PROP) {
         //	if (!variates){
         InitKepPredictor(m_kpPredictor, initOrbit);
         m_predictor = &m_kpPredictor;
@@ -68,7 +68,7 @@ bool Propagator::Init(const PropInfoBase& propInfo, const SInitOrbit& initOrbit,
         //		TAssert(0);
         //	}
     }
-    else if (propInfo.GetType() == PROGNOZ_PROP) {
+    else if (propInfo.GetType() == proptypes::PROGNOZ_PROP) {
         //	if (!variates){
         InitPrognozPredictor(m_prPredictor, initOrbit);
         m_predictor = &m_prPredictor;
@@ -76,7 +76,7 @@ bool Propagator::Init(const PropInfoBase& propInfo, const SInitOrbit& initOrbit,
         //		TAssert(0);
         //	}
     }
-    else if (propInfo.GetType() == SGP_PROP) {
+    else if (propInfo.GetType() == proptypes::SGP_PROP) {
         //	if (!variates){
         InitSGPPredictor(m_sgpPredictor, initOrbit);
         m_predictor = &m_sgpPredictor;

utils/interpolate.cpp

@@ -55,7 +55,7 @@ void XInterpolator::propagate_dir(
     std::vector<IRes>& xyz
 )
 {
-    const bool variates = orbit.GetCovMtx().HasMtx() && orbit.GetMotionModel().GetType() == ST_PROP && calc_rnb;
+    const bool variates = orbit.GetCovMtx().HasMtx() && orbit.GetMotionModel().GetType() == proptypes::ST_PROP && calc_rnb;
     Propagator prop(orbit, variates);
     const double step = get_step(orbit);
 
@@ -80,10 +80,10 @@ void XInterpolator::propagate_dir(
             prop.Propagate(date, norb);
             std::map<std::string, double> covrnb = norb.GetCovMtx().GetCovRnbMap(norb.GetPhasePointJ2000().CoordsVel);
             xyz.emplace_back(
-                norb.GetPhasePointJ2000().CoordsVel, 
-                Vect3(covrnb["RR"], covrnb["NN"], covrnb["BB"]), 
-                date, 
-                age, 
+                norb.GetPhasePointJ2000().CoordsVel,
+                Vect3(sqrt(covrnb["RR"]), sqrt(covrnb["NN"]), sqrt(covrnb["BB"])),
+                date,
+                age,
                 orbit.GetId());
         }
 
@@ -177,6 +177,23 @@ void fill(const Vect6& vec, const double jdi, const double jdf, ei& ephi)
     }
 }
 
+void fill_rnb(const Vect3& vec, const double jdi, const double jdf, ei& ephi)
+{
+    ephi.nv++;
+
+    ephi.outv[0][ephi.nv] = vec.x * 1E6;
+    ephi.outv[1][ephi.nv] = vec.y * 1E6;
+    ephi.outv[2][ephi.nv] = vec.z * 1E6;
+
+    ephi.jdi[ephi.nv] = jdi;
+    ephi.jdf[ephi.nv] = jdf;
+
+    if (ephi.jdf[ephi.nv] > 1.0e0) {
+        ephi.jdi[ephi.nv] += 1.e0;
+        ephi.jdf[ephi.nv] -= 1.e0;
+    }
+}
+
 void finalize(ei& ephi)
 {
     (ephi.nv)++;
@@ -211,7 +228,7 @@ void XInterpolator::create(const std::vector<IRes>& xyz)
         init(m_rnb_ephi, xyz.size());
 
         for (const auto& p : xyz) {
-            fill(Vect6(p.get_rnb(), Vect3()), p.get_date().GetDays(), p.get_date().GetFraction(), m_rnb_ephi);
+            fill_rnb(p.get_rnb(), p.get_date().GetDays(), p.get_date().GetFraction(), m_rnb_ephi);
         }
 
         finalize(m_rnb_ephi);
@@ -262,3 +279,42 @@ Vect6 XInterpolator::get_pos(const TimeJD& date, const bool tks) const
         outvec[5] / 1E3
     };
 }
+
+Vect3 XInterpolator::get_rnb(const TimeJD& date, const bool tks) const
+{
+    if (!m_rnb_ephi.jdi) {
+        throw Exp() << "interpolation not initialised\n";
+    }
+
+    if (!(date >= m_beg && date <= m_end)) {
+        throw Exp()
+            << "interpolation date out of range\n"
+            << "date is " << DateToStr(date) << "\n"
+            << "inerval is " << DateToStr(m_beg) << " - " << DateToStr(m_end) << "\n";
+    }
+
+    double outvec[6];
+    const int fvel = 1;
+    int ircode;
+    const double jdi = date.GetDays();
+    const double sec = date.GetFraction() * 86400;
+    cinterp_pvac(m_rnb_ephi, jdi, sec, outvec, -1, &ircode);
+
+    if (ircode != 0) {
+        throw Exp() << "interpolation fail\n";
+    }
+
+    if (!tks) {
+        return {
+            outvec[0],
+            outvec[1],
+            outvec[2],
+        };
+    }
+
+    return {
+        outvec[0] / 1E6,
+        outvec[1] / 1E6,
+        outvec[2] / 1E6,
+    };
+}

utils/interpolate.h

@@ -13,6 +13,7 @@ public:
     XInterpolator(const SInitOrbit& orbit, const TimeJD& beg, const TimeJD& end, bool calc_rnb);
     XInterpolator(std::vector<IRes> xyz);
     Vect6 get_pos(const TimeJD& date, bool tks = false) const;
+    Vect3 get_rnb(const TimeJD& date, bool tks = false) const;
     TimeJD get_beg() const { return m_beg; }
     TimeJD get_end() const { return m_end; }
     static void calc_positions(