diff --git a/src/dynamics/dynamics.cpp b/src/dynamics/dynamics.cpp
index 10c403e98..04d0886a7 100644
--- a/src/dynamics/dynamics.cpp
+++ b/src/dynamics/dynamics.cpp
@@ -29,7 +29,8 @@ void Dynamics::Initialize(const SimulationConfiguration* simulation_configuratio
                      local_celestial_information.GetGlobalInformation().GetCenterBodyGravityConstant_m3_s2(), "ORBIT", relative_information);
   attitude_ = InitAttitude(simulation_configuration->spacecraft_file_list_[spacecraft_id], orbit_, &local_celestial_information,
                            simulation_time->GetAttitudeRkStepTime_s(), structure->GetKinematicsParameters().GetInertiaTensor_b_kgm2(), spacecraft_id);
-  temperature_ = InitTemperature(simulation_configuration->spacecraft_file_list_[spacecraft_id], simulation_time->GetThermalRkStepTime_s());
+  temperature_ = InitTemperature(simulation_configuration->spacecraft_file_list_[spacecraft_id], simulation_time->GetThermalRkStepTime_s(),
+                                 &(local_environment_->GetSolarRadiationPressure()));
 
   // To get initial value
   orbit_->UpdateByAttitude(attitude_->GetQuaternion_i2b());
diff --git a/src/dynamics/thermal/heater.cpp b/src/dynamics/thermal/heater.cpp
index 2da6e158f..37b956121 100644
--- a/src/dynamics/thermal/heater.cpp
+++ b/src/dynamics/thermal/heater.cpp
@@ -1,10 +1,12 @@
 #include "heater.hpp"
 
+#include <cassert>
 #include <cmath>
 
 using namespace std;
 
 Heater::Heater(const unsigned int heater_id, const double power_rating_W) : heater_id_(heater_id), power_rating_W_(power_rating_W) {
+  AssertHeaterParams();
   heater_status_ = HeaterStatus::kOff;
   power_output_W_ = 0;
 }
@@ -26,3 +28,8 @@ void Heater::PrintParam(void) {
   cout << "  status            : " << static_cast<int>(heater_status_) << endl;
   cout << "  power output      : " << power_output_W_ << endl;
 }
+
+void Heater::AssertHeaterParams(void) {
+  assert(heater_id_ >= 1);       // Heater ID must be larger than 1
+  assert(power_rating_W_ >= 0);  // Power rating must be larger than 0
+}
diff --git a/src/dynamics/thermal/heater.hpp b/src/dynamics/thermal/heater.hpp
index a923fa52d..6472a6615 100644
--- a/src/dynamics/thermal/heater.hpp
+++ b/src/dynamics/thermal/heater.hpp
@@ -83,6 +83,12 @@ class Heater {
    * @brief Print power_rating_W, heater_status_, power_output_W_
    */
   void PrintParam(void);
+
+  /**
+   * @fn AssertHeaterParams
+   * @brief Check Heater Parameters
+   */
+  void AssertHeaterParams(void);
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_HEATER_HPP_
diff --git a/src/dynamics/thermal/heater_controller.cpp b/src/dynamics/thermal/heater_controller.cpp
index 8bb9e6db0..10280528d 100644
--- a/src/dynamics/thermal/heater_controller.cpp
+++ b/src/dynamics/thermal/heater_controller.cpp
@@ -1,11 +1,14 @@
 #include "heater_controller.hpp"
 
+#include <cassert>
 #include <cmath>
 
 using namespace std;
 
 HeaterController::HeaterController(const double lower_threshold_degC, const double upper_threshold_degC)
-    : lower_threshold_degC_(lower_threshold_degC), upper_threshold_degC_(upper_threshold_degC) {}
+    : lower_threshold_degC_(lower_threshold_degC), upper_threshold_degC_(upper_threshold_degC) {
+  AssertHeaterControllerParams();
+}
 
 HeaterController::~HeaterController() {}
 
@@ -21,3 +24,5 @@ void HeaterController::ControlHeater(Heater* p_heater, double temperature_degC)
     }
   }
 }
+
+void HeaterController::AssertHeaterControllerParams() { assert(upper_threshold_degC_ > lower_threshold_degC_); }
diff --git a/src/dynamics/thermal/heater_controller.hpp b/src/dynamics/thermal/heater_controller.hpp
index b678d6da1..db271ec84 100644
--- a/src/dynamics/thermal/heater_controller.hpp
+++ b/src/dynamics/thermal/heater_controller.hpp
@@ -10,6 +10,7 @@
 #include <library/logger/logger.hpp>
 #include <string>
 #include <vector>
+
 #include "heater.hpp"
 
 class HeaterController {
@@ -72,6 +73,12 @@ class HeaterController {
    * @param[in] temperature_degC
    */
   void ControlHeater(Heater* p_heater, double temperature_degC);
+
+  /**
+   * @fn AssertHeaterControllerParams
+   * @brief Check Parameters of HeaterController
+   */
+  void AssertHeaterControllerParams(void);
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_HEATER_CONTROLLER_HPP_
diff --git a/src/dynamics/thermal/heatload.cpp b/src/dynamics/thermal/heatload.cpp
index 6aa33b67f..b0ac5e106 100644
--- a/src/dynamics/thermal/heatload.cpp
+++ b/src/dynamics/thermal/heatload.cpp
@@ -1,5 +1,6 @@
 #include "heatload.hpp"
 
+#include <cassert>
 #include <cmath>
 #include <environment/global/physical_constants.hpp>
 
@@ -15,8 +16,8 @@ Heatload::Heatload(int node_id, std::vector<double> time_table_s, std::vector<do
   heater_heatload_W_ = 0.0;
   total_heatload_W_ = 0.0;
 
-  // [FIXME] Include assertion for size of time_table_s_ and internal_heatload_table_W_ to be larger than 2
-  // [FIXME] Include assertion for size of time_table_s_ and internal_heatload_table_W_ to be the same
+  AssertHeatloadParams();
+
   time_table_period_s_ = time_table_s_[time_table_s_.size() - 1];
   residual_elapsed_time_s_ = fmod(elapsed_time_s_, time_table_period_s_);
 }
@@ -48,3 +49,11 @@ void Heatload::SetElapsedTime_s(double elapsed_time_s) {
     }
   }
 }
+
+void Heatload::AssertHeatloadParams() {
+  // size of time_table_s_ and internal_heatload_table_W_ must be larger than 1
+  assert(time_table_s_.size() >= 1);
+  assert(internal_heatload_table_W_.size() >= 1);
+  // size of time_table_s_ and internal_heatload_table_W_ must be the same
+  assert(time_table_s_.size() == internal_heatload_table_W_.size());
+}
diff --git a/src/dynamics/thermal/heatload.hpp b/src/dynamics/thermal/heatload.hpp
index ad8272d70..805b1456b 100644
--- a/src/dynamics/thermal/heatload.hpp
+++ b/src/dynamics/thermal/heatload.hpp
@@ -99,6 +99,12 @@ class Heatload {
    * @param[in] heater_heatload_W
    */
   inline void SetHeaterHeatload_W(double heater_heatload_W) { heater_heatload_W_ = heater_heatload_W; }
+
+  /**
+   * @fn AssertHeatloadParams
+   * @brief Check if Heatload Parameters are Correct
+   */
+  void AssertHeatloadParams();
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_HEATLOAD_HPP_
diff --git a/src/dynamics/thermal/initialize_heater.cpp b/src/dynamics/thermal/initialize_heater.cpp
index c00c7fc19..fb488ad54 100644
--- a/src/dynamics/thermal/initialize_heater.cpp
+++ b/src/dynamics/thermal/initialize_heater.cpp
@@ -5,12 +5,26 @@
 
 #include "initialize_heater.hpp"
 
-// [FIXME] Write discription of csv file for heater initialization
+#include <cassert>
+
+/* Import heater properties by reading CSV File (heaters.csv)
+
+[File Formats of heater.csv]
+column 1: Heater_id(int, Use values larger than or equal to 1)
+column 2: Power Rating (double, [W])
+column 3: Lower threshold of control (double, [degC])
+column 4: Upper threshold of control (double, [degC])
+
+First row is for Header, data begins from the second row
+*/
 
 Heater InitHeater(const std::vector<std::string>& heater_str) {
   using std::stod;
   using std::stoi;
 
+  size_t heater_str_size_defined = 4;                    // Correct size of heater_str
+  assert(heater_str.size() == heater_str_size_defined);  // Check if size of heater_str is correct
+
   int heater_id = 0;
   double power_rating_W = 0;  // [W]
 
@@ -28,6 +42,9 @@ Heater InitHeater(const std::vector<std::string>& heater_str) {
 HeaterController InitHeaterController(const std::vector<std::string>& heater_str) {
   using std::stod;
 
+  size_t heater_str_size_defined = 4;                    // Correct size of heater_str
+  assert(heater_str.size() == heater_str_size_defined);  // Check if size of heater_str is correct
+
   // Index to read from heater_str for each parameter
   int index_lower_threshold = 2;
   int index_upper_threshold = 3;
diff --git a/src/dynamics/thermal/initialize_heatload.cpp b/src/dynamics/thermal/initialize_heatload.cpp
index c76b79be5..a0e16176e 100644
--- a/src/dynamics/thermal/initialize_heatload.cpp
+++ b/src/dynamics/thermal/initialize_heatload.cpp
@@ -5,17 +5,40 @@
 
 #include "initialize_heatload.hpp"
 
-// [FIXME] Write discription of csv file for heatload initialization
+#include <cassert>
+
+/* Import heatload properties by reading CSV File (heatload.csv)
+
+[File Formats of heatload.csv]
+First Row: Array of Simulation Times [s]
+After Second Row:
+  First Column: Node ID
+  After Second Column: Array of Heatload Values [W]
+
+Heatload object calculates heatload value of certain time by interpolating values of the heatload array
+
+Data Example:
+NodeID/Time, 0, 99, 100, 300, 500
+0,          10, 10,  20,  10, 20
+
+In this example, heatload at time 200s is calculated as 15W from interpolation
+
+The number of nodes and heatloads must be the same
+
+*/
 
 Heatload InitHeatload(const std::vector<std::string>& time_str, const std::vector<std::string>& internal_heatload_str) {
   using std::stod;
   using std::stoi;
 
-  // [FIXME] Include assertion for size of time_str and internal_heatload_str to be larger than 2
+  assert(time_str.size() >= 2);                             // Number of columns must be larger than 2
+  assert(internal_heatload_str.size() >= 2);                // Number of columns must be larger than 2
+  assert(time_str.size() == internal_heatload_str.size());  // Number of columns must be same
+
   std::vector<double> time_table_s(time_str.size() - 1);                            // exclude index
   std::vector<double> internal_heatload_table_W(internal_heatload_str.size() - 1);  // exclude index
 
-  int node_id = stoi(internal_heatload_str[0]);                                     // First data of internal_heatload_str is node id
+  int node_id = stoi(internal_heatload_str[0]);  // First data of internal_heatload_str is node id
 
   // read table
   for (unsigned int i = 0; i < time_str.size() - 1; ++i) {
diff --git a/src/dynamics/thermal/initialize_node.cpp b/src/dynamics/thermal/initialize_node.cpp
index 920a10225..25b51775c 100644
--- a/src/dynamics/thermal/initialize_node.cpp
+++ b/src/dynamics/thermal/initialize_node.cpp
@@ -5,15 +5,16 @@
 
 #include "initialize_node.hpp"
 
+#include <cassert>
 #include <iostream>
 #include <library/initialize/initialize_file_access.hpp>
 #include <string>
 #include <typeinfo>
 #include <vector>
 
-/* Import node properties by reading CSV File (Node.csv)
+/* Import node properties by reading CSV File (node.csv)
 
-[File Formats of Node.csv]
+[File Formats of node.csv]
 column 1: Node_id(int)
 column 2: Node_label(string)
 column 3: Node_type(int, 0: diffusive, 1: boundary), Arithmetic node to be implemented as future work
@@ -38,6 +39,9 @@ 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
+  assert(node_str.size() == node_str_size_defined);  // Check if size of node_str is correct
+
   int node_id = 0;                  // node number
   std::string node_label = "temp";  // node name
   int node_type_int = 0;            // node type
@@ -69,6 +73,15 @@ Node InitNode(const std::vector<std::string>& node_str) {
   for (int i = 0; i < 3; i++) {
     normal_v_b[i] = stod(node_str[index_normal_v_b_head + i]);
   }
+
+  // Normalize Norm Vector (Except for Boundary and Arithmetic Nodes)
+  if (node_type_int == 0) {
+    double norm = normal_v_b.CalcNorm();
+    for (int i = 0; i < 3; i++) {
+      normal_v_b[i] = normal_v_b[i] / norm;
+    }
+  }
+
   temperature_K = stod(node_str[index_temperature]);
 
   NodeType node_type = NodeType::kDiffusive;
diff --git a/src/dynamics/thermal/initialize_temperature.cpp b/src/dynamics/thermal/initialize_temperature.cpp
index 9aede7736..ca68bb419 100644
--- a/src/dynamics/thermal/initialize_temperature.cpp
+++ b/src/dynamics/thermal/initialize_temperature.cpp
@@ -5,6 +5,7 @@
 
 #include "initialize_temperature.hpp"
 
+#include <cassert>
 #include <environment/global/simulation_time.hpp>
 #include <library/initialize/initialize_file_access.hpp>
 #include <string>
@@ -57,7 +58,7 @@ First row is time data
 using std::string;
 using std::vector;
 
-Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s) {
+Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s, const SolarRadiationPressureEnvironment* srp_environment) {
   auto mainIni = IniAccess(file_name);
 
   vector<Node> node_list;
@@ -70,6 +71,7 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
   vector<vector<string>> heater_str_list;    // string vector of heater property data
   vector<vector<string>> heatload_str_list;  // string vector of heatload property data
   size_t node_num = 1;
+  size_t heatload_num = 1;
   size_t heater_num = 1;
 
   bool is_calc_enabled = mainIni.ReadEnable("THERMAL", "calculation");
@@ -100,6 +102,7 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
   /*since we don't know the number of node_list yet, set node_num=100 temporary.
     Recall that Nodes_num are given to this function only to reseve memory*/
 
+  heatload_num = heatload_str_list.size() - 1;
   auto times_itr = heatload_str_list.begin();  // First Row is Time Data
   for (auto itr = heatload_str_list.begin() + 1; itr != heatload_str_list.end(); ++itr) {
     heatload_list.push_back(InitHeatload(*times_itr, *itr));
@@ -118,6 +121,8 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
     node_list.push_back(InitNode(*itr));
   }
 
+  assert(node_num == heatload_num);  // Number of nodes and heatload lists must be the same
+
   // Read Heater Properties from CSV File
   string filepath_heater = file_path + "heaters.csv";
   IniAccess conf_heater(filepath_heater);
@@ -141,8 +146,13 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s
   conf_cij.ReadCsvDoubleWithHeader(conductance_matrix, node_num, 1, 1);
   conf_rij.ReadCsvDoubleWithHeader(radiation_matrix, node_num, 1, 1);
 
+  assert(conductance_matrix.size() == node_num);                      // Dimension must be same as node_num
+  assert(radiation_matrix.size() == node_num);                        // Dimension must be same as node_num
+  assert(conductance_matrix.size() == conductance_matrix[0].size());  // Must be square matrix
+  assert(radiation_matrix.size() == radiation_matrix[0].size());      // Must be square matrix
+
   Temperature* temperature;
   temperature = new Temperature(conductance_matrix, radiation_matrix, node_list, heatload_list, heater_list, heater_controller_list, node_num,
-                                rk_prop_step_s, is_calc_enabled, solar_calc_setting, debug);
+                                rk_prop_step_s, srp_environment, is_calc_enabled, solar_calc_setting, debug);
   return temperature;
 }
diff --git a/src/dynamics/thermal/initialize_temperature.hpp b/src/dynamics/thermal/initialize_temperature.hpp
index 1b248a1e7..ee841f0bd 100644
--- a/src/dynamics/thermal/initialize_temperature.hpp
+++ b/src/dynamics/thermal/initialize_temperature.hpp
@@ -15,6 +15,6 @@
  * @param[in] rk_prop_step_s: time step interval for temperature propagation integration
  * @return Temperature*
  */
-Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s);
+Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s, const SolarRadiationPressureEnvironment* srp_environment);
 
 #endif  // S2E_DYNAMICS_THERMAL_INITIALIZE_TEMPERATURE_HPP_
diff --git a/src/dynamics/thermal/node.cpp b/src/dynamics/thermal/node.cpp
index 6eef4c865..586bd6c51 100644
--- a/src/dynamics/thermal/node.cpp
+++ b/src/dynamics/thermal/node.cpp
@@ -5,6 +5,7 @@
 
 #include "node.hpp"
 
+#include <cassert>
 #include <cmath>
 #include <environment/global/physical_constants.hpp>
 
@@ -20,25 +21,20 @@ Node::Node(const int node_id, const string node_name, const NodeType node_type,
       capacity_J_K_(capacity_J_K),
       alpha_(alpha),
       area_m2_(area_m2),
-      normal_vector_b_(normal_vector_b),
-      node_type_(node_type) {
+      node_type_(node_type),
+      normal_vector_b_(normal_vector_b) {
+  AssertNodeParams();
   solar_radiation_W_ = 0;
 }
 
 Node::~Node() {}
 
-double Node::CalcSolarRadiation_W(libra::Vector<3> sun_direction_b) {
-  // FIXME: constants
-  double R = 6.96E+8;                              // Distance from sun
-  double T = 5778;                                 // sun surface temperature [K]
-  double sigma = 5.67E-8;                          // stephan-boltzman constant
-  double a = sun_direction_b.CalcNorm();           // distance from sun
-  double S = pow((R / a), 2) * sigma * pow(T, 4);  // Solar Constant at s/c position S[W/m^2]
-  double cos_theta = InnerProduct(sun_direction_b, normal_vector_b_) / a / normal_vector_b_.CalcNorm();
+double Node::CalcSolarRadiation_W(libra::Vector<3> sun_direction_b, double solar_flux_W_m2) {
+  double cos_theta = InnerProduct(sun_direction_b, normal_vector_b_);
 
   // calculate sun_power
   if (cos_theta > 0)
-    solar_radiation_W_ = S * area_m2_ * alpha_ * cos_theta;
+    solar_radiation_W_ = solar_flux_W_m2 * area_m2_ * alpha_ * cos_theta;
   else
     solar_radiation_W_ = 0;
   return solar_radiation_W_;
@@ -66,3 +62,12 @@ void Node::PrintParam(void) {
   }
   cout << endl;
 }
+
+void Node::AssertNodeParams(void) {
+  assert(node_id_ >= 0);                                      // Node ID should be larger than 0
+  assert(heater_id_ >= 0);                                    // Heater ID should be larger than 0
+  assert(temperature_K_ >= environment::absolute_zero_degC);  // Temperature must be larger than zero kelvin
+  assert(capacity_J_K_ >= 0);                                 // Capacity must be larger than 0, use 0 when node is boundary or arithmetic
+  assert(0 <= alpha_ && alpha_ <= 1);                         // alpha must be between 0 and 1
+  assert(area_m2_ >= 0);                                      // Area must be larger than 0
+}
diff --git a/src/dynamics/thermal/node.hpp b/src/dynamics/thermal/node.hpp
index 0b96e506e..6999b09ce 100644
--- a/src/dynamics/thermal/node.hpp
+++ b/src/dynamics/thermal/node.hpp
@@ -70,7 +70,7 @@ class Node {
    * @param sun_direction_b: Sun direction in body frame
    * @return double: Solar Radiation [W]
    */
-  double CalcSolarRadiation_W(libra::Vector<3> sun_direction_b);  // 太陽入射熱を計算
+  double CalcSolarRadiation_W(libra::Vector<3> sun_direction_b, double solar_flux_W_m2);
 
   // Getter
   /**
@@ -137,6 +137,12 @@ class Node {
    * @brief Print parameters of node in debug output
    */
   void PrintParam(void);
+
+  /**
+   * @fn AssertNodeParams
+   * @brief Check if Node Parameters are Correct
+   */
+  void AssertNodeParams(void);
 };
 
 #endif  // S2E_DYNAMICS_THERMAL_NODE_HPP_
diff --git a/src/dynamics/thermal/temperature.cpp b/src/dynamics/thermal/temperature.cpp
index 74cccba2c..1c31f7398 100644
--- a/src/dynamics/thermal/temperature.cpp
+++ b/src/dynamics/thermal/temperature.cpp
@@ -15,7 +15,8 @@ using namespace std;
 
 Temperature::Temperature(const vector<vector<double>> conductance_matrix_W_K, const vector<vector<double>> radiation_matrix_m2, vector<Node> nodes,
                          vector<Heatload> heatloads, vector<Heater> heaters, vector<HeaterController> heater_controllers, const int node_num,
-                         const double propagation_step_s, const bool is_calc_enabled, const SolarCalcSetting solar_calc_setting, const bool debug)
+                         const double propagation_step_s, const SolarRadiationPressureEnvironment* srp_environment, const bool is_calc_enabled,
+                         const SolarCalcSetting solar_calc_setting, const bool debug)
     : conductance_matrix_W_K_(conductance_matrix_W_K),
       radiation_matrix_m2_(radiation_matrix_m2),
       nodes_(nodes),
@@ -23,7 +24,8 @@ Temperature::Temperature(const vector<vector<double>> conductance_matrix_W_K, co
       heaters_(heaters),
       heater_controllers_(heater_controllers),
       node_num_(node_num),
-      propagation_step_s_(propagation_step_s),  // ルンゲクッタ積分時間刻み幅
+      propagation_step_s_(propagation_step_s),
+      srp_environment_(srp_environment),
       is_calc_enabled_(is_calc_enabled),
       solar_calc_setting_(solar_calc_setting),
       debug_(debug) {
@@ -44,8 +46,13 @@ Temperature::Temperature() {
 
 Temperature::~Temperature() {}
 
-void Temperature::Propagate(libra::Vector<3> sun_direction_b, const double time_end_s) {
+void Temperature::Propagate(libra::Vector<3> sun_position_b_m, const double time_end_s) {
   if (!is_calc_enabled_) return;
+  double sun_distance_m = sun_position_b_m.CalcNorm();
+  libra::Vector<3> sun_direction_b;
+  for (int i = 0; i < 3; i++) {
+    sun_direction_b[i] = sun_position_b_m[i] / sun_distance_m;
+  }
   while (time_end_s - propagation_time_s_ - propagation_step_s_ > 1.0e-6) {
     CalcRungeOneStep(propagation_time_s_, propagation_step_s_, sun_direction_b, node_num_);
     propagation_time_s_ += propagation_step_s_;
@@ -65,9 +72,8 @@ void Temperature::Propagate(libra::Vector<3> sun_direction_b, const double time_
       cout << setprecision(4) << itr->GetSolarRadiation_W() << "  ";
     }
     cout << "SunDir:  ";
-    double norm_sun_direction_b = sun_direction_b.CalcNorm();
     for (int i = 0; i < 3; i++) {
-      cout << setprecision(3) << sun_direction_b[i] / norm_sun_direction_b << "  ";
+      cout << setprecision(3) << sun_direction_b[i] << "  ";
     }
     cout << "Heatload:  ";
     for (auto itr = heatloads_.begin(); itr != heatloads_.end(); ++itr) {
@@ -113,7 +119,7 @@ void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra:
   }
 }
 
-vector<double> Temperature::CalcTemperatureDifferentials(vector<double> temperatures_K, double t, libra::Vector<3> sun_direction, int node_num) {
+vector<double> Temperature::CalcTemperatureDifferentials(vector<double> temperatures_K, double t, libra::Vector<3> sun_direction_b, int node_num) {
   // TODO: consider the following unused arguments are really needed
   UNUSED(temperatures_K);
 
@@ -121,8 +127,9 @@ vector<double> Temperature::CalcTemperatureDifferentials(vector<double> temperat
   for (int i = 0; i < node_num; i++) {
     heatloads_[i].SetElapsedTime_s(t);
     if (nodes_[i].GetNodeType() == NodeType::kDiffusive) {
+      double solar_flux_W_m2 = srp_environment_->GetPowerDensity_W_m2();
       if (solar_calc_setting_ == SolarCalcSetting::kEnable) {
-        double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_direction);
+        double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_direction_b, solar_flux_W_m2);
         heatloads_[i].SetSolarHeatload_W(solar_radiation_W);
       }
       double heater_power_W = GetHeaterPower_W(i);
@@ -157,7 +164,6 @@ double Temperature::GetHeaterPower_W(int node_id) {
 }
 
 void Temperature::UpdateHeaterStatus(void) {
-  // [FIXME] Heater status doesn't get updated...
   for (auto itr = nodes_.begin(); itr != nodes_.end(); ++itr) {
     int heater_id = itr->GetHeaterId();
     if (heater_id > 0) {
diff --git a/src/dynamics/thermal/temperature.hpp b/src/dynamics/thermal/temperature.hpp
index 72962a41b..ca0af4265 100644
--- a/src/dynamics/thermal/temperature.hpp
+++ b/src/dynamics/thermal/temperature.hpp
@@ -6,6 +6,7 @@
 #ifndef S2E_DYNAMICS_THERMAL_TEMPERATURE_HPP_
 #define S2E_DYNAMICS_THERMAL_TEMPERATURE_HPP_
 
+#include <environment/local/solar_radiation_pressure_environment.hpp>
 #include <library/logger/loggable.hpp>
 #include <string>
 #include <vector>
@@ -38,10 +39,11 @@ class Temperature : public ILoggable {
   std::vector<HeaterController> heater_controllers_;         // vector of heater controllers
   int node_num_;                                             // number of nodes
   double propagation_step_s_;                                // propagation step [s]
-  double propagation_time_s_;            // Incremented time inside class Temperature [s], finish propagation when reaching end_time
-  bool is_calc_enabled_;                 // Whether temperature calculation is enabled
-  SolarCalcSetting solar_calc_setting_;  // setting for solar calculation
-  bool debug_;                           // Activate debug output or not
+  double propagation_time_s_;  // Incremented time inside class Temperature [s], finish propagation when reaching end_time
+  const SolarRadiationPressureEnvironment* srp_environment_;  // SolarRadiationPressureEnvironment for calculating solar flux
+  bool is_calc_enabled_;                                      // Whether temperature calculation is enabled
+  SolarCalcSetting solar_calc_setting_;                       // setting for solar calculation
+  bool debug_;                                                // Activate debug output or not
 
   /**
    * @fn CalcRungeOneStep
@@ -49,7 +51,7 @@ class Temperature : public ILoggable {
    *
    * @param[in] time_now_s: Current elapsed time [s]
    * @param[in] time_step_s: Time step of RK4 [s]
-   * @param[in] sun_direction_b: Sun direction in body frame
+   * @param[in] sun_direction_b: Sun position in body frame [m]
    * @param[in] node_num: Number of nodes
    */
   void CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_direction_b, int node_num);
@@ -59,7 +61,7 @@ class Temperature : public ILoggable {
    *
    * @param temperatures_K: [UNUSED] Temperatures of each node [K]
    * @param time_now_s: Current elapsed time [s]
-   * @param sun_direction_b: Sun direction in body frame
+   * @param[in] sun_direction_b: Sun direction in body frame
    * @param node_num: Number of nodes
    * @return std::vector<double>: Differential of thermal equilibrium equation at time now
    */
@@ -85,8 +87,8 @@ class Temperature : public ILoggable {
    */
   Temperature(const std::vector<std::vector<double>> conductance_matrix_W_K, const std::vector<std::vector<double>> radiation_matrix_m2,
               std::vector<Node> nodes, std::vector<Heatload> heatloads, std::vector<Heater> heaters, std::vector<HeaterController> heater_controllers,
-              const int node_num, const double propagation_step_s, const bool is_calc_enabled, const SolarCalcSetting solar_calc_setting,
-              const bool debug);
+              const int node_num, const double propagation_step_s, const SolarRadiationPressureEnvironment* srp_environment,
+              const bool is_calc_enabled, const SolarCalcSetting solar_calc_setting, const bool debug);
   /**
    * @fn Temperature
    * @brief Construct a new Temperature object, used when thermal calculation is disabled.
@@ -103,10 +105,10 @@ class Temperature : public ILoggable {
    * @fn Propagate
    * @brief Propagate thermal calculation until time_end_s
    *
-   * @param sun_direction_b: Sun direction in body frame
+   * @param[in] sun_position_b_m: Sun position in body frame [m]
    * @param time_end_s: Time to finish propagation [s]
    */
-  void Propagate(libra::Vector<3> sun_direction_b, const double time_end_s);
+  void Propagate(libra::Vector<3> sun_position_b_m, const double time_end_s);
 
   // Getter
   /**
diff --git a/src/environment/global/physical_constants.hpp b/src/environment/global/physical_constants.hpp
index 40a79daea..c2acd7de0 100644
--- a/src/environment/global/physical_constants.hpp
+++ b/src/environment/global/physical_constants.hpp
@@ -23,6 +23,7 @@ DEFINE_PHYSICAL_CONSTANT(speed_of_light_m_s, 299792458.0L)                   //!
 DEFINE_PHYSICAL_CONSTANT(boltzmann_constant_J_K, 1.380649e-23L)              //!< Boltzmann constant [J/K]
 DEFINE_PHYSICAL_CONSTANT(stefan_boltzmann_constant_W_m2K4, 5.670374419e-8L)  //!< Stefan-Boltzmann constant [W/m2K4]
 DEFINE_PHYSICAL_CONSTANT(absolute_zero_degC, -273.15L)                       //!< Absolute zero [degC]
+DEFINE_PHYSICAL_CONSTANT(solar_constant_W_m2, 1.366e3L)                      //!< Solar constant [W/m2]
 }  // namespace physics
 
 inline namespace astronomy {
diff --git a/src/environment/local/solar_radiation_pressure_environment.cpp b/src/environment/local/solar_radiation_pressure_environment.cpp
index 846e42342..449ecd56a 100644
--- a/src/environment/local/solar_radiation_pressure_environment.cpp
+++ b/src/environment/local/solar_radiation_pressure_environment.cpp
@@ -20,6 +20,8 @@ SolarRadiationPressureEnvironment::SolarRadiationPressureEnvironment(LocalCelest
   sun_radius_m_ = local_celestial_information_->GetGlobalInformation().GetMeanRadiusFromName_m("SUN");
 }
 
+SolarRadiationPressureEnvironment::SolarRadiationPressureEnvironment() {}
+
 void SolarRadiationPressureEnvironment::UpdateAllStates() {
   if (!IsCalcEnabled) return;
 
diff --git a/src/environment/local/solar_radiation_pressure_environment.hpp b/src/environment/local/solar_radiation_pressure_environment.hpp
index c886a5364..0e9a26fab 100644
--- a/src/environment/local/solar_radiation_pressure_environment.hpp
+++ b/src/environment/local/solar_radiation_pressure_environment.hpp
@@ -23,6 +23,13 @@ class SolarRadiationPressureEnvironment : public ILoggable {
    * @param [in] local_celestial_information: Local celestial information
    */
   SolarRadiationPressureEnvironment(LocalCelestialInformation* local_celestial_information);
+
+  /**
+   * @fn SolarRadiationPressureEnvironment
+   * @brief Constructor (Use for constructing Temperature with no params)
+   */
+  SolarRadiationPressureEnvironment();
+
   /**
    * @fn ~SolarRadiationPressureEnvironment
    * @brief Destructor