commit f3f7677108051b1d1bf873700eb19b61f5a670fb Author: Milo Priegnitz Date: Fri Jun 21 18:20:06 2024 +0200 Add USSA 1976 model diff --git a/README.md b/README.md new file mode 100644 index 0000000..e98e052 --- /dev/null +++ b/README.md @@ -0,0 +1,4 @@ +# Models +Collecion of Altitude models. + +- U.S. Standard Atmosphere 1976 \ No newline at end of file diff --git a/include/sta/models/altitude.hpp b/include/sta/models/altitude.hpp new file mode 100644 index 0000000..1cd283a --- /dev/null +++ b/include/sta/models/altitude.hpp @@ -0,0 +1,58 @@ +#ifndef INC_MODELS_ALTITUDE_ +#define INC_MODELS_ALTITUDE_ + +#include + +namespace sta +{ +namespace models { + class BarometricFormula + { + /** + * @brief Get the altitude from the pressure using the barometric formula. The model is valid for altitudes up to leat 9 km. + * + * @param p Pressure in hPa + * @return Altitude in meters + */ + public: + /* + * @brief Constructor for the Barometric Formula + * @param pRef Reference pressure in hPa + * @param hRef Reference altitude in meters + */ + BarometricFormula(float pRef = 1013.25, float hRef = 0); + + /* + * @brief Get the altitude from the pressure using the barometric formula + * @param p Pressure in hPa + * @return Altitude in meters + */ + float getAltitude(float p); + + private: + float p0_; + }; + + namespace USSA1976 + { + /** + * @brief Get the altitude from the pressure using the U.S. Standard Atmosphere 1976 model (https://ntrs.nasa.gov/search.jsp?R=19770009539) + * and the barometric formula for multiple layers. The model is valid for altitudes up to at leat 71 km. See https://en.wikipedia.org/wiki/Barometric_formula for more information. + * + * @param p Pressure in hPa + * @return Altitude in meters + */ + float getAltitude(float p); + /** + * @brief Get the altitude as Gaussian Random Variable from the pressure as Gaussian Random Variable using the U.S. Standard Atmosphere 1976 model (https://ntrs.nasa.gov/search.jsp?R=19770009539) + * and the barometric formula for multiple layers. The model is valid for altitudes up to at leat 71 km. See https://en.wikipedia.org/wiki/Barometric_formula for more information. + * + * @param p Pressure in hPa as a Gaussian Random Variable + * @return Altitude in meters as a Gaussian Random Variable + */ + sta::math::probability::grv getAltitude(sta::math::probability::grv p); + } // namespace USSA1976 +} // namespace models +} // namespace sta + +#endif //INC_MODELS_ALTITUDE_ diff --git a/include/sta/models/constants.hpp b/include/sta/models/constants.hpp new file mode 100644 index 0000000..b7b15c0 --- /dev/null +++ b/include/sta/models/constants.hpp @@ -0,0 +1,12 @@ +#ifndef STA_MODELS_CONSTANTS_HPP +#define STA_MODELS_CONSTANTS_HPP + +//Altitude Constants + +#define STA_CONSTANT_R 8.3144598 //J/(mol*K) + +#define STA_CONSTANT_G_0 9.80665 //m/s^2 + +#define STA_CONSTANT_M 0.028964425278793993 //kg/mol + +#endif //STA_MODELS_CONSTANTS_HPP \ No newline at end of file diff --git a/library.json b/library.json new file mode 100644 index 0000000..cd3cac0 --- /dev/null +++ b/library.json @@ -0,0 +1,11 @@ +{ + "owner" : "sta", + "name": "sta-models", + "version": "0.1.0", + "dependencies": [ + { + "url": "git@gitlab.com:sta-git/alpaka/sta-peak.git", + "ref": "main" + } + ] +} \ No newline at end of file diff --git a/src/models/altitude.cpp b/src/models/altitude.cpp new file mode 100644 index 0000000..0bc1e47 --- /dev/null +++ b/src/models/altitude.cpp @@ -0,0 +1,153 @@ +#include +#include +#include +namespace sta +{ +namespace models { + +BarometricFormula::BarometricFormula(float pRef, float hRef) + { + p0_ = pRef / std::pow(1 - hRef / 44330, 5.255); + } + +float BarometricFormula::getAltitude(float p){ + return 44330 * (1 - std::pow(p / p0_, 1/5.255)); +} + +float USSA1976::getAltitude(float p){ + float T_b = 288.15; + float L_b = 0.0065; + float P_b = 1013.25; + float exp_b = 5.2561; + float H_b = 0; + + if(p > 226.321){ + T_b = 288.15; + L_b = 0.0065; + P_b = 1013.25; + exp_b = 5.2561; + H_b = 0; + } + else if(p > 54.7489){ + T_b = 216.65; + L_b = 0; + P_b = 226.321; + exp_b = 0; + H_b = 11000; + } + else if(p > 8.6802){ + T_b = 216.65; + L_b = -0.001; + P_b = 54.7489; + exp_b = -34.1626; + H_b = 20000; + } + else if(p > 1.1091){ + T_b = 228.65; + L_b = -0.0028; + P_b = 8.6802; + exp_b = -12.2009; + H_b = 32000; + } + else if(p > 0.6699){ + T_b = 270.65; + L_b = 0; + P_b = 1.1091; + exp_b = 0; + H_b = 47000; + } + else if(p > 0.0399){ + T_b = 270.65; + L_b = 0.0028; + P_b = 0.6699; + exp_b = 12.2009; + H_b = 51000; + } + else{ + T_b = 214.65; + L_b = 0.002; + P_b = 0.0399; + exp_b = 17.0816; + H_b = 71000; + } + + if(L_b == 0){ + return H_b - STA_CONSTANT_R*T_b*std::log(p/P_b) / (STA_CONSTANT_G_0* STA_CONSTANT_M); + } + return T_b/L_b * (-std::pow((p/P_b),(1/exp_b)) + 1) + H_b; +} + +sta::math::probability::grv USSA1976::getAltitude(sta::math::probability::grv p){ + float T_b = 288.15; + float L_b = 0.0065; + float P_b = 1013.25; + float exp_b = 5.2561; + float H_b = 0; + + if(p.mean > 226.321){ + T_b = 288.15; + L_b = 0.0065; + P_b = 1013.25; + exp_b = 5.2561; + H_b = 0; + } + else if(p.mean > 54.7489){ + T_b = 216.65; + L_b = 0; + P_b = 226.321; + exp_b = 0; + H_b = 11000; + } + else if(p.mean > 8.6802){ + T_b = 216.65; + L_b = -0.001; + P_b = 54.7489; + exp_b = -34.1626; + H_b = 20000; + } + else if(p.mean > 1.1091){ + T_b = 228.65; + L_b = -0.0028; + P_b = 8.6802; + exp_b = -12.2009; + H_b = 32000; + } + else if(p.mean > 0.6699){ + T_b = 270.65; + L_b = 0; + P_b = 1.1091; + exp_b = 0; + H_b = 47000; + } + else if(p.mean > 0.0399){ + T_b = 270.65; + L_b = 0.0028; + P_b = 0.6699; + exp_b = 12.2009; + H_b = 51000; + } + else{ + T_b = 214.65; + L_b = 0.002; + P_b = 0.0399; + exp_b = 17.0816; + H_b = 71000; + } + + float h; + float variance_h; + if(L_b == 0){ + h = H_b - STA_CONSTANT_R*T_b*std::log(p.mean/P_b) / (STA_CONSTANT_G_0* STA_CONSTANT_M); + variance_h = std::pow(( STA_CONSTANT_R * T_b / (STA_CONSTANT_G_0* STA_CONSTANT_M*p.mean)),2) * p.variance; + }else{ + float inner = -std::pow((p.mean/P_b),(1/exp_b)); + h = T_b/L_b * (inner + 1) + H_b; + variance_h = std::pow(inner*(T_b/L_b)*(1/exp_b) / p.mean ,2) * p.variance; + } + return {h, variance_h}; +} + + + +} // namespace models +} // namespace sta