Feature/add earth ir and fix albedo to thermal calc

example/settings/sample_satellite/local_environment.ini

@@ -22,6 +22,11 @@ calculation = DISABLE
 // Earth albedo factor: Percentage of sunlight reflected off the Earth surface. Value between 0.0 and 1.0
 earth_albedo_factor = 0.3
 
+[EARTH_INFRARED]
+calculation = ENABLE
+// Earth infrared temperature hot/cold side [K]. Default value is Global effective temperature of 250K
+earth_infrared_temperature_hot_side = 250.0
+earth_infrared_temperature_cold_side = 250.0
 
 [ATMOSPHERE]
 calculation = ENABLE

example/settings/sample_satellite/thermal_csv_files/heatload.csv

@@ -1,4 +1,4 @@
 NodeID/Times[s],0,500,501,1000,1001,1500
 0,10,10,20,20,10,10
 1,0,0,0,0,0,0
-2,0,0,0,0,0,0
+2,0,0,0,0,0,0

example/settings/sample_satellite/thermal_csv_files/node.csv

@@ -1,4 +1,4 @@
-Node_id,Node_label,"node_type (0: diffusive, 1: boundary, 2: arithmetic)",heater_node_id,capacity[J/K],alpha,area[m^2],normal_v_b_x,normal_v_b_y,normal_v_b_z,initial_temperature[K]
-0,BUS,0,1,880,0.2,0.06,1,0,0,300
-1,SAP,0,0,100,0,0.02,0,0,1,250
-2,SPACE,1,0,0,0,0,0,0,0,2.73
+Node_id,Node_label,node_type (0: diffusive 1: boundary 2: arithmetic),heater_node_id,capacity[J/K],alpha,epsilon,area[m^2],normal_v_b_x,normal_v_b_y,normal_v_b_z,initial_temperature[K]
+0,BUS,0,1,880,0.2,0.7,0.06,1,0,0,300
+1,SAP,0,0,100,0,0,0.02,0,0,1,250
+2,SPACE,1,0,0,0,0,0,0,0,0,2.73

src/dynamics/dynamics.cpp

@@ -31,7 +31,8 @@ void Dynamics::Initialize(const simulation::SimulationConfiguration* simulation_
                                      simulation_time->GetAttitudeRkStepTime_s(), structure->GetKinematicsParameters().GetInertiaTensor_b_kgm2(),
                                      spacecraft_id);
   temperature_ = thermal::InitTemperature(simulation_configuration->spacecraft_file_list_[spacecraft_id], simulation_time->GetThermalRkStepTime_s(),
-                                          &(local_environment_->GetSolarRadiationPressure()), &(local_environment_->GetEarthAlbedo()));
+                                          &(local_environment_->GetSolarRadiationPressure()), &(local_environment_->GetEarthAlbedo()),
+                                          &(local_environment_->GetEarthInfrared()));
 
   // To get initial value
   orbit_->UpdateByAttitude(attitude_->GetQuaternion_i2b());

src/dynamics/thermal/heatload.hpp

@@ -41,7 +41,9 @@ class Heatload {
    * @fn UpdateTotalHeatload
    * @brief Update total heatload value by summing up all factors
    */
-  void UpdateTotalHeatload(void) { total_heatload_W_ = solar_heatload_W_ + albedo_heatload_W_ + internal_heatload_W_ + heater_heatload_W_; }
+  void UpdateTotalHeatload(void) {
+    total_heatload_W_ = solar_heatload_W_ + albedo_heatload_W_ + earth_infrared_heatload_W_ + internal_heatload_W_ + heater_heatload_W_;
+  }
 
   // Getter
   /**
@@ -54,6 +56,11 @@ class Heatload {
    * @brief Return Albedo Heatload
    */
   inline double GetAlbedoHeatload_W(void) const { return albedo_heatload_W_; }
+  /**
+   * @fn GetEarthInfraredHeatload_W
+   * @brief Return Earth Infrared Heatload
+   */
+  inline double GetEarthInfraredHeatload_W(void) const { return earth_infrared_heatload_W_; }
   /**
    * @fn GetInternalHeatload_W
    * @brief Return Internal Heatload
@@ -91,6 +98,11 @@ class Heatload {
    * @param[in] albedo_heatload_W
    */
   inline void SetAlbedoHeatload_W(const double albedo_heatload_W) { albedo_heatload_W_ = albedo_heatload_W; }
+  /**
+   * @brief Set Earth Infrared Heatload [W]
+   * @param[in] earth_infrared_heatload_W
+   */
+  inline void SetEarthInfraredHeatload_W(const double earth_infrared_heatload_W) { earth_infrared_heatload_W_ = earth_infrared_heatload_W; }
   /**
    * @brief Set Heater Heatload [W]
    * @param[in] heater_heatload_W
@@ -103,12 +115,13 @@ class Heatload {
   std::vector<double> time_table_s_;               //!< Times that internal heatload values are defined [s]
   std::vector<double> internal_heatload_table_W_;  //!< Defined internal heatload values [W]
 
-  unsigned int elapsed_time_idx_;  //!< index of time_table_s_ that is closest to elapsed_time_s_
-  double solar_heatload_W_;        //!< Heatload from solar flux [W]
-  double albedo_heatload_W_;       //!< Heatload from albedo flux [W]
-  double internal_heatload_W_;     //!< Heatload from internal dissipation [W]
-  double heater_heatload_W_;       //!< Heatload from heater [W]
-  double total_heatload_W_;        //!< Total heatload [W]
+  unsigned int elapsed_time_idx_;     //!< index of time_table_s_ that is closest to elapsed_time_s_
+  double solar_heatload_W_;           //!< Heatload from solar flux [W]
+  double albedo_heatload_W_;          //!< Heatload from albedo flux [W]
+  double earth_infrared_heatload_W_;  //!< Heatload from earth infrared flux [W]
+  double internal_heatload_W_;        //!< Heatload from internal dissipation [W]
+  double heater_heatload_W_;          //!< Heatload from heater [W]
+  double total_heatload_W_;           //!< Total heatload [W]
 
   double time_table_period_s_;      //!< Value of last element of time_table_s_, which represents the period of the heatload table [s]
   double residual_elapsed_time_s_;  //!< Residual of dividing elapsed_time_s_ by time_table_period_s_ [s]

src/dynamics/thermal/node.cpp

@@ -8,26 +8,29 @@
 #include <cassert>
 #include <cmath>
 #include <environment/global/physical_constants.hpp>
+#include <math_physics/math/constants.hpp>
 
 using namespace std;
 using namespace s2e::math;
 
 namespace s2e::dynamics::thermal {
 
 Node::Node(const size_t node_id, const string node_name, const NodeType node_type, const size_t heater_id, const double temperature_ini_K,
-           const double capacity_J_K, const double alpha, const double area_m2, math::Vector<3> normal_vector_b)
+           const double capacity_J_K, const double alpha, const double epsilon, const double area_m2, math::Vector<3> normal_vector_b)
     : node_id_(node_id),
       node_name_(node_name),
       heater_id_(heater_id),
       temperature_K_(temperature_ini_K),
       capacity_J_K_(capacity_J_K),
       alpha_(alpha),
+      epsilon_(epsilon),
       area_m2_(area_m2),
       node_type_(node_type),
       normal_vector_b_(normal_vector_b) {
   AssertNodeParams();
   solar_radiation_W_ = 0.0;
   albedo_radiation_W_ = 0.0;
+  earth_infrared_W_ = 0.0;
 }
 
 Node::~Node() {}
@@ -43,20 +46,91 @@ double Node::CalcSolarRadiation_W(math::Vector<3> sun_direction_b, double solar_
   return solar_radiation_W_;
 }
 
-double Node::CalcAlbedoRadiation_W(math::Vector<3> earth_position_b_m, double earth_albedo_W_m2) {
-  math::Vector<3> earth_direction_b = earth_position_b_m.CalcNormalizedVector();
+double Node::CalcAlbedoRadiation_W(math::Vector<3> earth_position_b_m, math::Vector<3> sun_direction_b, double earth_albedo_W_m2) {
+  math::Vector<3> sat2earth_direction_b = earth_position_b_m.CalcNormalizedVector();
+  math::Vector<3> sat2sun_direction_b = sun_direction_b;
 
-  double cos_theta_albedo = InnerProduct(earth_direction_b, normal_vector_b_);
+  math::Vector<3> vec_a = -sat2earth_direction_b;
+  math::Vector<3> vec_b = (sat2sun_direction_b - sat2earth_direction_b).CalcNormalizedVector();
 
-  // albedo radiation calculation; earth_albedo_W_m2 reflects the shadow coefficient.
-  if (cos_theta_albedo > 0.0) {
-    albedo_radiation_W_ = earth_albedo_W_m2 * area_m2_ * alpha_ * cos_theta_albedo;
+  double cos_theta = InnerProduct(vec_a, vec_b);
+  double lamda = acos(InnerProduct(sat2earth_direction_b, normal_vector_b_));
+  double h = earth_position_b_m.CalcNorm() - environment::earth_equatorial_radius_m;
+  double H = earth_position_b_m.CalcNorm() / environment::earth_equatorial_radius_m;
+  double phi_m = asin(1.0 / H);
+  double b = sqrt(H * H - 1.0);
+  double view_factor;
+
+  // Calc view factor
+  // ref)POWER INPUT TO A SMALL FLAT PLATE FROM A DIFFUSELY RADIATING SPHERE WITH APPLICATION TO EARTH SATELLITES: THE SPINNING PLATE
+  if (h < 1732e3) {
+    if (lamda <= math::pi / 2.0 - phi_m) {
+      view_factor = cos(lamda) / H / H;
+    } else if (lamda <= math::pi / 2.0 + phi_m) {
+      view_factor = 1.0 / 2.0 - 1.0 / math::pi * asin(b / (H * sin(lamda))) +
+                    1.0 / (math::pi * pow(H, 2.0)) * (cos(lamda) * acos(-b / tan(lamda)) - b * sqrt(1.0 - pow(H * cos(lamda), 2.0)));
+    } else {
+      view_factor = 0.0;
+    }
   } else {
-    albedo_radiation_W_ = 0.0;
+    if (lamda < math::pi / 2.0) {
+      // Consider in terms of steradian
+      view_factor = 0.25 / H / H;
+    } else {
+      view_factor = 0.0;
+    }
   }
+
+  // Banister's approximation. ref) RADIATION GEOMETRY FACTOR BETWEEN THE EARTH AND A SATELLITE
+  if (cos_theta > 0.0) {
+    // TODO: correlate the value of the exponent with the view factor
+    view_factor *= pow(cos_theta, 3.0);
+  } else {
+    view_factor = 0.0;
+  }
+
+  // albedo radiation calculation; earth_albedo_W_m2 reflects the shadow coefficient.
+  albedo_radiation_W_ = earth_albedo_W_m2 * area_m2_ * alpha_ * view_factor;
+
   return albedo_radiation_W_;
 }
 
+double Node::CalcEarthInfraredRadiation_W(math::Vector<3> earth_position_b_m, double earth_infrared_W_m2) {
+  math::Vector<3> earth_direction_b = earth_position_b_m.CalcNormalizedVector();
+
+  double lamda = acos(InnerProduct(earth_direction_b, normal_vector_b_));
+  double h = earth_position_b_m.CalcNorm() - environment::earth_equatorial_radius_m;
+  double H = earth_position_b_m.CalcNorm() / environment::earth_equatorial_radius_m;
+  double phi_m = asin(1.0 / H);
+  double b = sqrt(H * H - 1.0);
+  double view_factor;
+
+  // Calc view factor
+  // ref)POWER INPUT TO A SMALL FLAT PLATE FROM A DIFFUSELY RADIATING SPHERE WITH APPLICATION TO EARTH SATELLITES: THE SPINNING PLATE
+  if (h < 1732e3) {
+    if (lamda <= math::pi / 2.0 - phi_m) {
+      view_factor = cos(lamda) / H / H;
+    } else if (lamda <= math::pi / 2.0 + phi_m) {
+      view_factor = 1.0 / 2.0 - 1.0 / math::pi * asin(b / (H * sin(lamda))) +
+                    1.0 / (math::pi * pow(H, 2.0)) * (cos(lamda) * acos(-b / tan(lamda)) - b * sqrt(1.0 - pow(H * cos(lamda), 2.0)));
+    } else {
+      view_factor = 0.0;
+    }
+  } else {
+    if (lamda < math::pi / 2.0) {
+      // Consider in terms of steradian
+      view_factor = 0.25 / H / H;
+    } else {
+      view_factor = 0.0;
+    }
+  }
+
+  // earth infrared radiation calculation
+  earth_infrared_W_ = earth_infrared_W_m2 * area_m2_ * epsilon_ * view_factor;
+
+  return earth_infrared_W_;
+}
+
 void Node::PrintParam(void) {
   string node_type_str = "";
   if (node_type_ == NodeType::kDiffusive) {
@@ -70,6 +144,7 @@ void Node::PrintParam(void) {
   cout << "  node_name    : " << node_name_ << endl;
   cout << "  temperature  : " << temperature_K_ << endl;
   cout << "  alpha        : " << alpha_ << endl;
+  cout << "  epsilon      : " << epsilon_ << endl;
   cout << "  capacity     : " << capacity_J_K_ << endl;
   cout << "  node type    : " << node_type_str << endl;
   cout << "  heater id    : " << heater_id_ << endl;
@@ -99,6 +174,14 @@ void Node::AssertNodeParams(void) {
     std::cerr << "The value is set as 0.0." << std::endl;
     alpha_ = 0.0;
   }
+
+  // epsilon must be between 0 and 1
+  if (epsilon_ < 0.0 || epsilon_ > 1.0) {
+    std::cerr << "[WARNING] node: epsilon is over the range [0, 1]." << std::endl;
+    std::cerr << "The value is set as 0.0." << std::endl;
+    epsilon_ = 0.0;
+  }
+
   // Area must be larger than 0
   if (area_m2_ < 0.0) {
     std::cerr << "[WARNING] node: area is less than zero [m2]." << std::endl;
@@ -116,16 +199,17 @@ column 3: Node_type(int, 0: diffusive, 1: boundary), Arithmetic node to be imple
 column 4: Heater No
 column 5: capacity
 column 6: alpha
-column 7: area
-column 8,9,10: normal vector of surface(body frame)
-column 11: initial temperature(K)
+column 7: epsilon
+column 8: area
+column 9,10,11: normal vector of surface(body frame)
+column 12: initial temperature(K)
 
 First row is for Header, data begins from the second row
 Ex.
-Node_id,Node_label,node_type,heater_id,capacity,alpha,area,normal_v_b_x,normal_v_b_y,normal_v_b_z,initial_temperature
-0,BUS,0,1,880,0.2,0.06,1,0,0,300
-1,SAP,0,0,100,0,0.02,0,0,1,250
-2,SPACE,1,0,0,0,0,0,0,0,2.73
+Node_id,Node_label,node_type,heater_id,capacity,alpha,epsilon,area,normal_v_b_x,normal_v_b_y,normal_v_b_z,initial_temperature
+0,BUS,0,1,880,0.2,0.1,0.06,1,0,0,300
+1,SAP,0,0,100,0,0.8,0.02,0,0,1,250
+2,SPACE,1,0,0,0,0,0,0,0,0,2.73
 
 Be sure to include at least one boundary node to avoid divergence
 */
@@ -134,7 +218,7 @@ Node InitNode(const std::vector<std::string>& node_str) {
   using std::stod;
   using std::stoi;
 
-  size_t node_str_size_defined = 11;                 // Correct size of node_str
+  size_t node_str_size_defined = 12;                 // Correct size of node_str
   assert(node_str.size() == node_str_size_defined);  // Check if size of node_str is correct
 
   size_t node_id = 0;               // node number
@@ -144,6 +228,7 @@ Node InitNode(const std::vector<std::string>& node_str) {
   double temperature_K = 0.0;       // [K]
   double capacity_J_K = 0.0;        // [J/K]
   double alpha = 0.0;               // []
+  double epsilon = 0.0;             // []
   double area_m2 = 0.0;             // [m^2]
 
   // Index to read from node_str for each parameter
@@ -153,16 +238,18 @@ Node InitNode(const std::vector<std::string>& node_str) {
   size_t index_heater_id = 3;
   size_t index_capacity = 4;
   size_t index_alpha = 5;
-  size_t index_area = 6;
-  size_t index_normal_v_b_head = 7;
-  size_t index_temperature = 10;
+  size_t index_epsilon = 6;
+  size_t index_area = 7;
+  size_t index_normal_v_b_head = 8;
+  size_t index_temperature = 11;
 
   node_id = stoi(node_str[index_node_id]);
   node_label = node_str[index_node_label];
   node_type_int = stoi(node_str[index_node_type]);
   heater_id = stoi(node_str[index_heater_id]);
   capacity_J_K = stod(node_str[index_capacity]);
   alpha = stod(node_str[index_alpha]);
+  epsilon = stod(node_str[index_epsilon]);
   area_m2 = stod(node_str[index_area]);
   math::Vector<3> normal_v_b;
   for (size_t i = 0; i < 3; i++) {
@@ -185,7 +272,7 @@ Node InitNode(const std::vector<std::string>& node_str) {
   } else if (node_type_int == 2) {
     node_type = NodeType::kArithmetic;
   }
-  Node node(node_id, node_label, node_type, heater_id, temperature_K, capacity_J_K, alpha, area_m2, normal_v_b);
+  Node node(node_id, node_label, node_type, heater_id, temperature_K, capacity_J_K, alpha, epsilon, area_m2, normal_v_b);
   return node;
 }
 

src/dynamics/thermal/node.hpp

@@ -43,11 +43,12 @@ class Node {
    * @param[in] temperature_ini_K: Initial temperature of node [K]
    * @param[in] capacity_J_K: Heat capacity of node [J/K]
    * @param[in] alpha: Solar absorptivity of face with possibility of solar incidence
+   * @param[in] epsilon: Emissivity of face with possibility of infrared radiation
    * @param[in] area_m2: Area of face with possibility of solar incidence [m^2]
    * @param[in] normal_vector_b: Normal vector of face with possibility of solar incidence (Body frame)
    */
   Node(const size_t node_id, const std::string node_name, const NodeType node_type, const size_t heater_id, const double temperature_ini_K,
-       const double capacity_J_K, const double alpha, const double area_m2, math::Vector<3> normal_vector_b);
+       const double capacity_J_K, const double alpha, const double epsilon, const double area_m2, math::Vector<3> normal_vector_b);
   /**
    * @fn ~Node
    * @brief Destroy the Node object
@@ -67,10 +68,20 @@ class Node {
    * @brief Calculate albedo radiation [W] from earth direction, albedo factor, area, and normal vector
    *
    * @param earth_position_b_m: Earth position in body frame
+   * @param sun_direction_b: Sun direction in body frame
    * @param earth_albedo_W_m2: Earth albedo [W/m^2]
    * @return double: Albedo Radiation [W]
    */
-  double CalcAlbedoRadiation_W(math::Vector<3> earth_position_b_m, double earth_albedo_W_m2);
+  double CalcAlbedoRadiation_W(math::Vector<3> earth_position_b_m, math::Vector<3> sun_direction_b, double earth_albedo_W_m2);
+  /**
+   * @fn CalcEarthInfraredRadiation_W
+   * @brief Calculate Earth Infrared Radiation [W] from earth direction, earth infrared radiation, area, and normal vector
+   *
+   * @param earth_position_b_m: Earth position in body frame
+   * @param earth_infrared_W_m2: Earth Infrared Radiation [W/m^2]
+   * @return double: Earth Infrared Radiation [W]
+   */
+  double CalcEarthInfraredRadiation_W(math::Vector<3> earth_position_b_m, double earth_infrared_W_m2);
 
   // Getter
   /**
@@ -121,6 +132,12 @@ class Node {
    * @return double: Albedo Radiation [W]
    */
   inline double GetAlbedoRadiation_W(void) const { return albedo_radiation_W_; }
+  /**
+   * @fn GetEarthInfraredRadiation_W
+   * @brief Return Earth Infrared Radiation [W]
+   * @return double: Earth Infrared Radiation [W]
+   */
+  inline double GetEarthInfraredRadiation_W(void) const { return earth_infrared_W_; }
   /**
    * @fn GetNodeType
    * @brief Return Node Type
@@ -151,9 +168,11 @@ class Node {
   double temperature_K_;
   double capacity_J_K_;
   double alpha_;
+  double epsilon_;
   double area_m2_;
   double solar_radiation_W_;
   double albedo_radiation_W_;
+  double earth_infrared_W_;
   NodeType node_type_;
   math::Vector<3> normal_vector_b_;