diff --git a/data/sample/initialize_files/sample_local_environment.ini b/data/sample/initialize_files/sample_local_environment.ini
index 6c6e55a95..fce7453c5 100644
--- a/data/sample/initialize_files/sample_local_environment.ini
+++ b/data/sample/initialize_files/sample_local_environment.ini
@@ -10,6 +10,11 @@ magnetic_field_white_noise_standard_deviation_nT = 50.0
 [SOLAR_RADIATION_PRESSURE_ENVIRONMENT]
 calculation = ENABLE
 logging = ENABLE
+// The number of shadow generating bodies other than the central body
+number_of_third_shadow_source_ = 1
+// List of shadow generating bodies other than the central body
+// All these bodies must be included in the "selected_body_name" of "[CelestialInformation]"
+third_shadow_source_name(0) = MOON
 
 
 [ATMOSPHERE]
diff --git a/src/environment/local/solar_radiation_pressure_environment.cpp b/src/environment/local/solar_radiation_pressure_environment.cpp
index 846e42342..686b2a650 100644
--- a/src/environment/local/solar_radiation_pressure_environment.cpp
+++ b/src/environment/local/solar_radiation_pressure_environment.cpp
@@ -16,7 +16,7 @@
 SolarRadiationPressureEnvironment::SolarRadiationPressureEnvironment(LocalCelestialInformation* local_celestial_information)
     : local_celestial_information_(local_celestial_information) {
   solar_radiation_pressure_N_m2_ = solar_constant_W_m2_ / environment::speed_of_light_m_s;
-  shadow_source_name_ = local_celestial_information_->GetGlobalInformation().GetCenterBodyName();
+  shadow_source_name_list_.push_back(local_celestial_information_->GetGlobalInformation().GetCenterBodyName());
   sun_radius_m_ = local_celestial_information_->GetGlobalInformation().GetMeanRadiusFromName_m("SUN");
 }
 
@@ -24,7 +24,10 @@ void SolarRadiationPressureEnvironment::UpdateAllStates() {
   if (!IsCalcEnabled) return;
 
   UpdatePressure();
-  CalcShadowCoefficient(shadow_source_name_);
+  shadow_coefficient_ = 1.0;  // Initialize for multiple shadow source
+  for (auto shadow_source_name : shadow_source_name_list_) {
+    CalcShadowCoefficient(shadow_source_name);
+  }
 }
 
 void SolarRadiationPressureEnvironment::UpdatePressure() {
@@ -54,7 +57,7 @@ std::string SolarRadiationPressureEnvironment::GetLogValue() const {
 
 void SolarRadiationPressureEnvironment::CalcShadowCoefficient(std::string shadow_source_name) {
   if (shadow_source_name == "SUN") {
-    shadow_coefficient_ = 1.0;
+    shadow_coefficient_ *= 1.0;
     return;
   }
 
@@ -81,17 +84,17 @@ void SolarRadiationPressureEnvironment::CalcShadowCoefficient(std::string shadow
 
   if (c < fabs(a - b) && a <= b)  // The occultation is total (spacecraft is in umbra)
   {
-    shadow_coefficient_ = 0.0;
+    shadow_coefficient_ *= 0.0;
   } else if (c < fabs(a - b) && a > b)  // The occultation is partial but maximum
   {
     shadow_coefficient_ = 1.0 - (b * b) / (a * a);
   } else if (fabs(a - b) <= c && c <= (a + b))  // spacecraft is in penumbra
   {
     double A = a * a * acos(x / a) + b * b * acos((c - x) / b) - c * y;  // The area of the occulted segment of the apparent solar disk
-    shadow_coefficient_ = 1.0 - A / (libra::pi * a * a);
+    shadow_coefficient_ *= 1.0 - A / (libra::pi * a * a);
   } else {  // no occultation takes place
     assert(c > (a + b));
-    shadow_coefficient_ = 1.0;
+    shadow_coefficient_ *= 1.0;
   }
 }
 
@@ -104,5 +107,11 @@ SolarRadiationPressureEnvironment InitSolarRadiationPressureEnvironment(std::str
   srp_env.IsCalcEnabled = conf.ReadEnable(section, INI_CALC_LABEL);
   srp_env.is_log_enabled_ = conf.ReadEnable(section, INI_LOG_LABEL);
 
+  size_t number_of_third_shadow_source = conf.ReadInt(section, "number_of_third_shadow_source");
+  std::vector<std::string> list = conf.ReadVectorString(section, "third_shadow_source_name", number_of_third_shadow_source);
+  for (size_t i = 0; i < number_of_third_shadow_source; i++) {
+    srp_env.AddShadowSource(list[0]);
+  }
+
   return srp_env;
 }
diff --git a/src/environment/local/solar_radiation_pressure_environment.hpp b/src/environment/local/solar_radiation_pressure_environment.hpp
index e50c06189..8fbc73250 100644
--- a/src/environment/local/solar_radiation_pressure_environment.hpp
+++ b/src/environment/local/solar_radiation_pressure_environment.hpp
@@ -36,6 +36,16 @@ class SolarRadiationPressureEnvironment : public ILoggable {
    */
   void UpdateAllStates();
 
+  /**
+   * @fn AddShadowSource
+   * @brief Update pressure and shadow coefficients
+   * @param [in] shadow_source_name: Shadow source name
+   */
+  void AddShadowSource(const std::string shadow_source_name) {
+    // TODO: Add assertion
+    shadow_source_name_list_.push_back(shadow_source_name);
+  }
+
   // Getter
   /**
    * @fn GetPressure_N_m2
@@ -81,11 +91,11 @@ class SolarRadiationPressureEnvironment : public ILoggable {
   virtual std::string GetLogValue() const;
 
  private:
-  double solar_radiation_pressure_N_m2_;  //!< Solar radiation pressure [N/m^2]
-  double solar_constant_W_m2_ = 1366.0;   //!< Solar constant [W/m^2] TODO: We need to change the value depends on sun activity.
-  double shadow_coefficient_ = 1.0;       //!< Shadow function
-  double sun_radius_m_;                   //!< Sun radius [m]
-  std::string shadow_source_name_;        //!< Shadow source name
+  double solar_radiation_pressure_N_m2_;              //!< Solar radiation pressure [N/m^2]
+  double solar_constant_W_m2_ = 1366.0;               //!< Solar constant [W/m^2] TODO: We need to change the value depends on sun activity.
+  double shadow_coefficient_ = 1.0;                   //!< Shadow function
+  double sun_radius_m_;                               //!< Sun radius [m]
+  std::vector<std::string> shadow_source_name_list_;  //!< Shadow source name list
 
   LocalCelestialInformation* local_celestial_information_;  //!< Local celestial information