Merge pull request #217 from ut-issl/hotfix/planet-select

Fix bug in celestial information

data/SampleSat/ini/SampleSimBase.ini

@@ -93,8 +93,14 @@ num_of_selected_body = 3
 selected_body(0) = EARTH
 selected_body(1) = SUN
 selected_body(2) = MOON
-selected_body(3) = MARS
-
+selected_body(3) = MERCURY
+selected_body(4) = VENUS
+selected_body(5) = MARS
+selected_body(6) = JUPITER
+selected_body(7) = SATURN
+selected_body(8) = URANUS
+selected_body(9) = NEPTUNE
+selected_body(10) = PLUTO
 
 [FURNSH_PATH]
 // CSPICE Kernel files definition

src/Environment/Global/CelestialInformation.cpp

@@ -1,8 +1,10 @@
 #include "CelestialInformation.h"
 
-#include <Interface/LogOutput/LogUtility.h>
 #include <SpiceUsr.h>
 
+#include <Interface/LogOutput/LogUtility.h>
+#include <string.h>
+
 #include <iostream>
 #include <sstream>
 
@@ -16,6 +18,7 @@ CelestialInformation::CelestialInformation(string inertial_frame, string aber_co
       aber_cor_(aber_cor),
       center_obj_(center_obj),
       rotation_mode_(rotation_mode) {
+  // Initialize list
   int num_of_state = num_of_selected_body_ * 3;
   celes_objects_pos_from_center_i_ = new double[num_of_state];
   celes_objects_vel_from_center_i_ = new double[num_of_state];
@@ -29,10 +32,11 @@ CelestialInformation::CelestialInformation(string inertial_frame, string aber_co
     SpiceInt dim;
     SpiceDouble gravity_constant;
     bodvcd_c(planet_id, "GM", 1, &dim, &gravity_constant);
-    // CONVERT FROM [km^3/s^2] to [m^3/s^2]
+    // Convert unit [km^3/s^2] to [m^3/s^2]
     celes_objects_gravity_constant_[i] = gravity_constant * 1E+9;
   }
 
+  // Acquisition of radius
   for (int i = 0; i < num_of_selected_body_; i++) {
     SpiceInt planet_id = selected_body_[i];
     SpiceInt dim;
@@ -50,6 +54,7 @@ CelestialInformation::CelestialInformation(string inertial_frame, string aber_co
     celes_objects_mean_radius_m_[i] = pow(rx * ry * rz, 1.0 / 3.0);
   }
 
+  // Initialize rotation
   EarthRotation_ = new CelestialRotation(rotation_mode_, center_obj_);
 }
 
@@ -89,20 +94,24 @@ CelestialInformation::~CelestialInformation() {
 }
 
 void CelestialInformation::UpdateAllObjectsInfo(const double current_jd) {
-  SpiceDouble rv_buf[6], lt, et;
-  SpiceBoolean found;
-  const int maxlen = 100;
-  char namebuf[maxlen];
+  // Convert time
+  SpiceDouble et;
   string jd = "jd " + to_string(current_jd);
   str2et_c(jd.c_str(), &et);
 
   for (int i = 0; i < num_of_selected_body_; i++) {
     SpiceInt planet_id = selected_body_[i];
+
     // Acquisition of body name from id
+    SpiceBoolean found;
+    const int maxlen = 100;
+    char namebuf[maxlen];
     bodc2n_c(planet_id, maxlen, namebuf, (SpiceBoolean*)&found);
+
     // Acquisition of position and velocity
-    spkezr_c(namebuf, et, inertial_frame_.c_str(), aber_cor_.c_str(), center_obj_.c_str(), (SpiceDouble*)rv_buf, (SpiceDouble*)&lt);
-    // CONVERT [km], [km/s] to [m], [m/s]
+    SpiceDouble rv_buf[6];
+    GetPlanetOrbit(namebuf, et, (SpiceDouble*)rv_buf);
+    // Convert unit [km], [km/s] to [m], [m/s]
     for (int j = 0; j < 3; j++) {
       celes_objects_pos_from_center_i_[i * 3 + j] = rv_buf[j] * 1000.0;
       celes_objects_vel_from_center_i_[i * 3 + j] = rv_buf[j + 3] * 1000.0;
@@ -113,6 +122,7 @@ void CelestialInformation::UpdateAllObjectsInfo(const double current_jd) {
   EarthRotation_->Update(current_jd);
 }
 
+// Getters
 Vector<3> CelestialInformation::GetPosFromCenter_i(const int id) const {
   Vector<3> pos(0.0);
   if (id > num_of_selected_body_) return pos;
@@ -127,13 +137,6 @@ Vector<3> CelestialInformation::GetVelFromCenter_i(const int id) const {
   return vel;
 }
 
-Vector<3> CelestialInformation::GetRadii(const int id) const {
-  Vector<3> radii(0.0);
-  if (id > num_of_selected_body_) return radii;
-  for (int i = 0; i < 3; i++) radii[i] = celes_objects_planetographic_radii_m_[id * 3 + i];
-  return radii;
-}
-
 Vector<3> CelestialInformation::GetPosFromCenter_i(const char* body_name) const {
   int id = CalcBodyIdFromName(body_name);
   return GetPosFromCenter_i(id);
@@ -149,6 +152,15 @@ double CelestialInformation::GetGravityConstant(const char* body_name) const {
   return celes_objects_gravity_constant_[index];
 }
 
+double CelestialInformation::GetCenterBodyGravityConstant_m3_s2(void) const { return GetGravityConstant(center_obj_.c_str()); }
+
+Vector<3> CelestialInformation::GetRadii(const int id) const {
+  Vector<3> radii(0.0);
+  if (id > num_of_selected_body_) return radii;
+  for (int i = 0; i < 3; i++) radii[i] = celes_objects_planetographic_radii_m_[id * 3 + i];
+  return radii;
+}
+
 Vector<3> CelestialInformation::GetRadiiFromName(const char* body_name) const {
   int id = CalcBodyIdFromName(body_name);
   return GetRadii(id);
@@ -229,3 +241,19 @@ void CelestialInformation::DebugOutput(void) {
          << ": " << celes_objects_gravity_constant_[i] << "\n";
   }
 }
+
+void CelestialInformation::GetPlanetOrbit(const char* planet_name, double et, double orbit[6]) {
+  // Add `BARYCENTER` if needed
+  const int maxlen = 100;
+  char planet_name_[maxlen];
+  strcpy(planet_name_, planet_name);
+  if (strcmp(planet_name, "MARS") == 0 || strcmp(planet_name, "JUPITER") == 0 || strcmp(planet_name, "SATURN") == 0 ||
+      strcmp(planet_name, "URANUS") == 0 || strcmp(planet_name, "NEPTUNE") == 0 || strcmp(planet_name, "PLUTO") == 0) {
+    strcat(planet_name_, "_BARYCENTER");
+  }
+
+  // Get orbit
+  SpiceDouble lt;
+  spkezr_c((ConstSpiceChar*)planet_name_, (SpiceDouble)et, (ConstSpiceChar*)inertial_frame_.c_str(), (ConstSpiceChar*)aber_cor_.c_str(), (ConstSpiceChar*)center_obj_.c_str(), (SpiceDouble*)orbit, (SpiceDouble*)&lt);
+  return;
+}

src/Environment/Global/CelestialInformation.h

@@ -16,60 +16,67 @@ using libra::Vector;
 
 class CelestialInformation : public ILoggable {
  public:
-  // CONSTRUCTOR OF CELESTIAL INFORMATION
   CelestialInformation(std::string inertial_frame, std::string aber_cor, std::string center_obj, RotationMode rotation_mode, int num_of_selected_body,
                        int* selected_body);
   CelestialInformation(const CelestialInformation& obj);
   virtual ~CelestialInformation();
 
-  // UPDATE THE ALL SELECTED CELESTIAL OBJECTS INFORMATION
+  // ILoggable
+  virtual std::string GetLogHeader() const;
+  virtual std::string GetLogValue() const;
+
+  // Update the all selected celestial objects information
   void UpdateAllObjectsInfo(const double current_jd);
 
-  // GET FUNCTIONS
+  // Getters
+  // Orbit information
   Vector<3> GetPosFromCenter_i(const int id) const;
   Vector<3> GetVelFromCenter_i(const int id) const;
-  Vector<3> GetRadii(const int id) const;
   Vector<3> GetPosFromCenter_i(const char* body_name) const;
   Vector<3> GetVelFromCenter_i(const char* body_name) const;
+  // Gravity constants
   double GetGravityConstant(const char* body_name) const;
+  double GetCenterBodyGravityConstant_m3_s2(void) const;
+  // Shape information
+  Vector<3> GetRadii(const int id) const;
   Vector<3> GetRadiiFromName(const char* body_name) const;
   double GetMeanRadiusFromName(const char* body_name) const;
+  // Parameters
   inline int GetNumBody(void) const { return num_of_selected_body_; }
   inline int* GetSelectedBody(void) const { return selected_body_; }
-  int CalcBodyIdFromName(const char* body_name) const;
   inline std::string GetCenterBodyName(void) const { return center_obj_; }
-
+  // Members
   inline CelestialRotation GetEarthRotation(void) const { return *EarthRotation_; };
 
-  // FOR LOG OUTPUT
-  virtual std::string GetLogHeader() const;
-  virtual std::string GetLogValue() const;
-
-  // FOR DEBUG OUTPUT
+  // Calculation
+  int CalcBodyIdFromName(const char* body_name) const;
   void DebugOutput(void);
 
  private:
-  int num_of_selected_body_;
-  int* selected_body_;          // IDs of selected bodies.
-  std::string inertial_frame_;  // Definition of inertial frame. Default = "J2000"
-  std::string aber_cor_;        // stellar aberration correction. Default =
-                                // "NONE"（Ref：http://fermi.gsfc.nasa.gov/ssc/library/fug/051108/Aberration_Julie.ppt）
-  std::string center_obj_;      // center object. Default = "EARTH"
-  RotationMode rotation_mode_;  // designation of dynamics model. Default = "Full"
+  // Setting parameters
+  int num_of_selected_body_;    //!< number of selected body
+  int* selected_body_;          //!< SPICE IDs of selected bodies
+  std::string inertial_frame_;  //!< Definition of inertial frame
+  std::string aber_cor_;        //!< Stellar aberration correction
+                                //（Ref：http://fermi.gsfc.nasa.gov/ssc/library/fug/051108/Aberration_Julie.ppt）
+  std::string center_obj_;      //!< Center object of inertial frame
 
-  // Global Information. POS:[m], VEL:[m/s], GRAVITY CONSTANT (G*M):[m^3/s^2]
-  double* celes_objects_pos_from_center_i_;
-  double* celes_objects_vel_from_center_i_;
-  double* celes_objects_gravity_constant_;
-  // 3 axis planetographic radii.
-  // X-axis pass through the 0 degree latitude 0 degree longitude direction
-  // Z-axis pass through the 90 degree latitude direction
-  // Y-axis equal to the cross product of the unit Z-axis and X-axis vectors
-  double* celes_objects_planetographic_radii_m_;
-  // Mean radius: r = (rx * ry * rz)^(1/3)
-  double* celes_objects_mean_radius_m_;
+  // Calculated values
+  double* celes_objects_pos_from_center_i_;       //!< Position vector list at inertial frame [m]
+  double* celes_objects_vel_from_center_i_;       //!< Velocity vector list at inertial frame [m/s]
+  double* celes_objects_gravity_constant_;        //!< Gravity constant list [m^3/s^2]
+  double* celes_objects_mean_radius_m_;           //!< mean radius list [m] r = (rx * ry * rz)^(1/3)
+  double* celes_objects_planetographic_radii_m_;  //!< 3 axis planetographic radii.
+                                                  // X-axis pass through the 0 degree latitude 0 degree longitude direction
+                                                  // Z-axis pass through the 90 degree latitude direction
+                                                  // Y-axis equal to the cross product of the unit Z-axis and X-axis vectors
 
   // Rotational Motion of each planets
   CelestialRotation* EarthRotation_;
+  RotationMode rotation_mode_;  //!< Designation of rotation model
+
+  // Override function of SPICE
+  void GetPlanetOrbit(const char* planet_name, double et, double orbit[6]);
 };
+
 #endif  //__celestial_information_H__

src/Environment/Global/InitGlobalEnvironment.cpp

@@ -61,23 +61,10 @@ CelestialInformation* InitCelesInfo(std::string file_name) {
   const char* section = "PLANET_SELECTION";
   const char* furnsh_section = "FURNSH_PATH";
 
-  //各種定義読み込み
+  // Read SPICE setting
   std::string inertial_frame = ini_file.ReadString(section, "inertial_frame");
   std::string aber_cor = ini_file.ReadString(section, "aberration_correction");
   std::string center_obj = ini_file.ReadString(section, "center_object");
-  RotationMode rotation_mode;
-  std::string rotation_mode_temp = ini_file.ReadString(section, "rotation_mode");
-  if (rotation_mode_temp == "Idle") {
-    rotation_mode = Idle;
-  } else if (rotation_mode_temp == "Simple") {
-    rotation_mode = Simple;
-  } else if (rotation_mode_temp == "Full") {
-    rotation_mode = Full;
-  } else  // if rotation_mode is neither Idle, Simple, nor Full, set
-          // rotation_mode to Idle
-  {
-    rotation_mode = Idle;
-  }
 
   // SPICE Furnsh
   std::vector<std::string> keywords = {"TLS", "TPC1", "TPC2", "TPC3", "BSP"};
@@ -86,23 +73,39 @@ CelestialInformation* InitCelesInfo(std::string file_name) {
     furnsh_c(fname.c_str());
   }
 
-  //天体情報をinitialize
+  // Initialize celestial body list
   const int num_of_selected_body = ini_file.ReadInt(section, "num_of_selected_body");
-  int* selected_body = new int[num_of_selected_body];  // これのdeleteはCelestialInformationに任せる
-  SpiceInt planet_id;
-  SpiceBoolean found;
+  int* selected_body = new int[num_of_selected_body];
   for (int i = 0; i < num_of_selected_body; i++) {
+    // Convert body name to SPICE ID
     std::string selected_body_i = "selected_body(" + std::to_string(i) + ")";
     char selected_body_temp[30];
     ini_file.ReadChar(section, selected_body_i.c_str(), 30, selected_body_temp);
+    SpiceInt planet_id;
+    SpiceBoolean found;
     bodn2c_c(selected_body_temp, (SpiceInt*)&planet_id, (SpiceBoolean*)&found);
-    std::string body_name = selected_body_temp;
 
     // If the object specified in the ini file is not found, exit the program.
     assert(found == SPICETRUE);
 
     selected_body[i] = planet_id;
   }
+
+  // Read Rotation setting
+  RotationMode rotation_mode;
+  std::string rotation_mode_temp = ini_file.ReadString(section, "rotation_mode");
+  if (rotation_mode_temp == "Idle") {
+    rotation_mode = Idle;
+  } else if (rotation_mode_temp == "Simple") {
+    rotation_mode = Simple;
+  } else if (rotation_mode_temp == "Full") {
+    rotation_mode = Full;
+  } else  // if rotation_mode is neither Idle, Simple, nor Full, set
+          // rotation_mode to Idle
+  {
+    rotation_mode = Idle;
+  }
+
   CelestialInformation* celestial_info;
   celestial_info = new CelestialInformation(inertial_frame, aber_cor, center_obj, rotation_mode, num_of_selected_body, selected_body);