Merge pull request #546 from ut-issl/develop

Merge develop to main for Minor update v7.1.0

.github/workflows/google-test.yml

@@ -95,7 +95,7 @@ jobs:
           CC:  ${{ steps.compiler.outputs.CC  }}
           CXX: ${{ steps.compiler.outputs.CXX }}
         run: |
-          cmake . -DEXT_LIB_DIR=./ExtLibraries -DBUILD_64BIT=ON -DGOOGLE_TEST=ON
+          cmake . -DEXT_LIB_DIR=./ExtLibraries -DBUILD_64BIT=ON -DGOOGLE_TEST=ON -DCORE_DIR_FROM_EXE=../s2e-core
           cmake --build .
 
       - name: run test

.github/workflows/validate-renovate.yml

@@ -12,7 +12,7 @@ jobs:
 
     steps:
       - uses: actions/checkout@v4
-      - uses: actions/setup-node@v3
+      - uses: actions/setup-node@v4
 
       - name: install
         run: |

CMakeLists.txt

@@ -226,17 +226,14 @@ if(GOOGLE_TEST)
 
   # Unit test
   set(TEST_PROJECT_NAME ${PROJECT_NAME}_TEST)
-  set(TEST_FILES
-    src/library/math/test_quaternion.cpp
-    src/library/math/test_vector.cpp
-    src/library/math/test_matrix.cpp
-    src/library/math/test_matrix_vector.cpp
-    src/library/math/test_s2e_math.cpp
-    src/library/numerical_integration/test_runge_kutta.cpp
-    src/library/gravity/test_gravity_potential.cpp
-  )
+
+  # Add all test_*.cpp files as SOURCE_FILES
+  file(GLOB_RECURSE TEST_FILES ${CMAKE_CURRENT_LIST_DIR}/src/test_*.cpp)
+  # Uncomment the following line to exclude any files that match the REGEX from TEST_FILES
+  # list(FILTER TEST_FILES EXCLUDE REGEX ${CMAKE_CURRENT_LIST_DIR}/src/test_example.cpp)
+
   add_executable(${TEST_PROJECT_NAME} ${TEST_FILES})
-  target_link_libraries(${TEST_PROJECT_NAME} gtest gtest_main)
+  target_link_libraries(${TEST_PROJECT_NAME} gtest gtest_main gmock)
   target_link_libraries(${TEST_PROJECT_NAME} LIBRARY)
   include_directories(${TEST_PROJECT_NAME})
   add_test(NAME s2e-test COMMAND ${TEST_PROJECT_NAME})
@@ -246,6 +243,7 @@ if(GOOGLE_TEST)
   set_target_properties(${TEST_PROJECT_NAME} PROPERTIES LANGUAGE CXX)
   set_target_properties(${TEST_PROJECT_NAME} PROPERTIES CXX_STANDARD 17)
   set_target_properties(${TEST_PROJECT_NAME} PROPERTIES CXX_EXTENSIONS FALSE)
+  target_compile_definitions(${TEST_PROJECT_NAME} PRIVATE "CORE_DIR_FROM_EXE=\"${CORE_DIR_FROM_EXE}\"")
 
 endif()
 

data/sample/initialize_files/components/angular_velocity_observer.ini

@@ -0,0 +1,27 @@
+[COMPONENT_BASE]
+prescaler = 1
+
+[SENSOR_BASE_ANGULAR_VELOCITY_OBSERVER]
+// Constant bias noise at component frame [rad/s]
+constant_bias_c_rad_s(0) = 0.0
+constant_bias_c_rad_s(1) = 0.0
+constant_bias_c_rad_s(2) = 0.0
+
+// Standard deviation for random walk noise[rad/s]
+random_walk_standard_deviation_c_rad_s(0) = 0.0
+random_walk_standard_deviation_c_rad_s(1) = 0.0
+random_walk_standard_deviation_c_rad_s(2) = 0.0
+
+// Limit of random walk noise[rad/s]
+random_walk_limit_c_rad_s(0) = 0.0
+random_walk_limit_c_rad_s(1) = 0.0
+random_walk_limit_c_rad_s(2) = 0.0
+
+// Standard deviation of normal random noise[rad/s]
+normal_random_standard_deviation_c_rad_s(0) = 1e-3
+normal_random_standard_deviation_c_rad_s(1) = 1e-3
+normal_random_standard_deviation_c_rad_s(2) = 1e-3
+
+// Range [rad/s]
+range_to_constant_rad_s = 5.0  // smaller than range_to_zero
+range_to_zero_rad_s = 10.0

data/sample/initialize_files/components/attitude_observer.ini

@@ -0,0 +1,7 @@
+[ATTITUDE_OBSERVER]
+// Standard deviation of force direction error [deg]
+error_angle_standard_deviation_deg = 1 
+
+[COMPONENT_BASE]
+// Prescaler with respect to the component update period
+prescaler = 1

data/sample/initialize_files/sample_gnss.ini

@@ -2,50 +2,66 @@
 directory_path = EXT_LIB_DIR_FROM_EXE/sp3/
 calculation = DISABLE
 
+// Choose from IGS, CODE_Final, JAXA_Final, QZSS_Final
 true_position_file_sort = IGS
-// choose from IGS, CODE_Final, JAXA_Final, QZSS_Final
+
 // As small as possible within the range.
 true_position_first = igs21610.sp3
 true_position_last = igs21613.sp3
+
+// Choose from 0: Lagrange, 1: Trigonometric : temporarily, only Trigonometric is valid, Lagrange is not yet implemented.
 true_position_interpolation_method = 1
-// choose from 0: Lagrange, 1: Trigonometric : temporarily, only Trigonometric is valid, Lagrange is not yet implemented.
+
+// If you chose Trigonometric as interpolation method, you must choose odds number here.
 true_position_interpolation_number = 9
-// if you chose Trigonometric as interpolation method, you must choose odds number here.
 
-// choose frome .sp3, .clk_30s or .clk
+// Choose frome .sp3, .clk_30s or .clk
 true_clock_file_extension = .clk_30s
+
+// Choose from IGS, CODE_Final, JAXA_Final, QZSS_Final
+// If you choose clk_30s, should choose the one equivalent to final, if you choose clk should choose the one equivalent to rapid.
 true_clock_file_sort = IGS
-// choose from IGS, CODE_Final, JAXA_Final, QZSS_Final
-//if you choose clk_30s, should choose the one equivalent to final, if you choose clk should choose the one equivalent to rapid.
+
 // As small as possible within the range.
-// The method is fixed with Lagrange interpolation, 3 (quadratic) recommended
 true_clock_first = igs21610.clk_30s
 true_clock_last = igs21613.clk_30s
+// The method is fixed with Lagrange interpolation, 3 (quadratic) recommended
 true_clock_interpolation_number = 3
 
-estimate_position_file_sort = madocaRT
+//
+// Estimated value
+//
 // choose from IGS, IGR, IGU,
 //            CODE_Final, CODE_Rapid
 //            JAXA_Final, JAXA_Rapid, JAXA_Ultra_rapid, madocaRT,
 //            QZSS_Final, QZSS_Rapid, QZSS_Ultra_rapid
+estimate_position_file_sort = madocaRT
+
 // As small as possible within the range.
 estimate_position_first = madoca21610.sp3
 estimate_position_last = madoca21613.sp3
-// when you use Ultra Rapid Product, you can choose "observe" or "predict"
-estimate_position_interpolation_method = 1
+
 // choose from 0: Lagrange, 1: Trigonometric : temporarily, only Trigonometric is valid, Lagrange is not yet implemented.
-estimate_position_interpolation_number = 9
+estimate_position_interpolation_method = 1
+
 // if you chose Trigonometric as interpolation method, you must choose odds number here.
-estimate_clock_file_sort = madocaRT
+estimate_position_interpolation_number = 9
+
 // choose from IGS, IGR, IGU,
 //            CODE_Final, CODE_Rapid
 //            JAXA_Final, JAXA_Rapid, JAXA_Ultra_rapid, madocaRT,
 //            QZSS_Final, QZSS_Rapid, QZSS_Ultra_rapid
+estimate_clock_file_sort = madocaRT
+
 // choose frome .sp3, .clk_30s or .clk
 estimate_clock_file_extension = .sp3
+
 // As small as possible within the range.
 estimate_clock_first = madoca21610.sp3
 estimate_clock_last = madoca21613.sp3
+
 // The method is fixed with Lagrange interpolation, 3 (quadratic) recommended
 estimate_clock_interpolation_number = 3
+
+// when you use Ultra Rapid Product, you can choose "observe" or "predict"
 estimate_ur_observe_or_predict = observe

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]

data/sample/initialize_files/sample_satellite.ini

@@ -144,5 +144,8 @@ rw_file = INI_FILE_DIR_FROM_EXE/components/reaction_wheel.ini
 thruster_file = INI_FILE_DIR_FROM_EXE/components/thruster.ini
 force_generator_file = INI_FILE_DIR_FROM_EXE/components/force_generator.ini
 torque_generator_file = INI_FILE_DIR_FROM_EXE/components/torque_generator.ini
+angular_velocity_observer_file = INI_FILE_DIR_FROM_EXE/components/angular_velocity_observer.ini
+attitude_observer_file = INI_FILE_DIR_FROM_EXE/components/attitude_observer.ini
 antenna_file = INI_FILE_DIR_FROM_EXE/components/spacecraft_antenna.ini
 component_interference_file = INI_FILE_DIR_FROM_EXE/components/component_interference.ini
+telescope_file = INI_FILE_DIR_FROM_EXE/components/telescope.ini

data/sample/initialize_files/sample_simulation_base.ini

@@ -1,6 +1,6 @@
 [TIME]
 // Simulation start time [UTC]
-simulation_start_time_utc = 2020/01/01 12:00:00.0
+simulation_start_time_utc = 2020/04/01 12:00:00.0
 
 // Simulation duration [sec]
 simulation_duration_s = 200

scripts/Plot/plot_angular_velocity_observer.py

@@ -0,0 +1,94 @@
+#
+# Plot Angular Velocity Observer
+#
+# arg[1] : read_file_tag : time tag for default CSV output log file. ex. 220627_142946
+#
+
+#
+# Import
+#
+# plots
+import matplotlib.pyplot as plt
+# local function
+from common import find_latest_log_tag
+from common import add_log_file_arguments
+from common import read_3d_vector_from_csv
+from common import read_scalar_from_csv
+# arguments
+import argparse
+
+# Arguments
+aparser = argparse.ArgumentParser()
+aparser = add_log_file_arguments(aparser)
+aparser.add_argument('--no-gui', action='store_true')
+args = aparser.parse_args()
+
+#
+# Read Arguments
+#
+# log file path
+path_to_logs = args.logs_dir
+
+read_file_tag = args.file_tag
+if read_file_tag == None:
+  print("file tag does not found. use latest.")
+  read_file_tag = find_latest_log_tag(path_to_logs)
+
+print("log: " + read_file_tag)
+
+#
+# CSV file name
+#
+read_file_name  = path_to_logs + '/' + 'logs_' + read_file_tag + '/' + read_file_tag + '_default.csv'
+
+#
+# Data read and edit
+#
+# Read S2E CSV
+time = read_scalar_from_csv(read_file_name, 'elapsed_time[s]')
+
+measured_angular_velocity_b_rad_s = read_3d_vector_from_csv(read_file_name, 'angular_velocity_observer_measured_value_b', 'rad/s')
+true_angular_velocity_b_rad_s = read_3d_vector_from_csv(read_file_name, 'spacecraft_angular_velocity_b', 'rad/s')
+
+# Statistics
+# We assume that the component frame and the body frame is same
+error_rad_s = measured_angular_velocity_b_rad_s - true_angular_velocity_b_rad_s
+average = [0.0, 0.0, 0.0]
+standard_deviation = [0.0, 0.0, 0.0]
+for i in range(3):
+  average[i] = error_rad_s[i].mean()
+  standard_deviation[i] = error_rad_s[i].std()
+
+#
+# Plot
+#
+unit = ' rad/s'
+
+fig, axis = plt.subplots(3, 1, squeeze = False, tight_layout = True, sharex = True)
+axis[0, 0].plot(time[0], measured_angular_velocity_b_rad_s[0], marker=".", c="red",    label="MEASURED-X")
+axis[0, 0].plot(time[0], true_angular_velocity_b_rad_s[0], marker=".", c="orange",  label="TRUE-X")
+axis[0, 0].legend(loc = 'upper right')
+axis[0, 0].text(0.01, 0.99, "Error average:" + format(average[0], '+.2e') + unit, verticalalignment = 'top', transform = axis[0, 0].transAxes)
+axis[0, 0].text(0.01, 0.89, "Standard deviation:" + format(standard_deviation[0], '+.2e') + unit, verticalalignment = 'top', transform = axis[0, 0].transAxes)
+
+axis[1, 0].plot(time[0], measured_angular_velocity_b_rad_s[1], marker=".", c="green",  label="MEASURED-Y")
+axis[1, 0].plot(time[0], true_angular_velocity_b_rad_s[1], marker=".", c="yellow",  label="TRUE-Y")
+axis[1, 0].legend(loc = 'upper right')
+axis[1, 0].text(0.01, 0.99, "Error average:" + format(average[1], '+.2e') + unit, verticalalignment = 'top', transform = axis[1, 0].transAxes)
+axis[1, 0].text(0.01, 0.89, "Standard deviation:" + format(standard_deviation[1], '+.2e') + unit, verticalalignment = 'top', transform = axis[1, 0].transAxes)
+
+axis[2, 0].plot(time[0], measured_angular_velocity_b_rad_s[2], marker=".", c="blue",   label="MEASURED-Z")
+axis[2, 0].plot(time[0], true_angular_velocity_b_rad_s[2], marker=".", c="purple",  label="TRUE-Z")
+axis[2, 0].legend(loc = 'upper right')
+axis[2, 0].text(0.01, 0.99, "Error average:" + format(average[2], '+.2e') + unit, verticalalignment = 'top', transform = axis[2, 0].transAxes)
+axis[2, 0].text(0.01, 0.89, "Standard deviation:" + format(standard_deviation[2], '+.2e') + unit, verticalalignment = 'top', transform = axis[2, 0].transAxes)
+
+fig.suptitle("Angular Velocity Observer")
+fig.supylabel("Angular Velocity" + unit)
+fig.supxlabel("Time [s]")
+
+# Data save
+if args.no_gui:
+  plt.savefig(read_file_tag + "_angular_velocity_observer.png") # save last figure only
+else:
+  plt.show()

scripts/Plot/plot_attitude_observer.py

@@ -0,0 +1,95 @@
+#
+# Plot Attitude Observer
+#
+# arg[1] : read_file_tag : time tag for default CSV output log file. ex. 220627_142946
+#
+
+#
+# Import
+#
+# plots
+import matplotlib.pyplot as plt
+# local function
+from common import find_latest_log_tag
+from common import add_log_file_arguments
+from common import read_quaternion_from_csv
+from common import read_scalar_from_csv
+from common import calc_error_angle_from_quaternions
+# arguments
+import argparse
+
+# Arguments
+aparser = argparse.ArgumentParser()
+aparser = add_log_file_arguments(aparser)
+aparser.add_argument('--no-gui', action='store_true')
+args = aparser.parse_args()
+
+#
+# Read Arguments
+#
+# log file path
+path_to_logs = args.logs_dir
+
+read_file_tag = args.file_tag
+if read_file_tag == None:
+  print("file tag does not found. use latest.")
+  read_file_tag = find_latest_log_tag(path_to_logs)
+
+print("log: " + read_file_tag)
+
+#
+# CSV file name
+#
+read_file_name  = path_to_logs + '/' + 'logs_' + read_file_tag + '/' + read_file_tag + '_default.csv'
+
+#
+# Data read and edit
+#
+# Read S2E CSV
+time = read_scalar_from_csv(read_file_name, 'elapsed_time[s]')
+
+measured_quaternion_i2b = read_quaternion_from_csv(read_file_name, 'attitude_observer_quaternion_i2b')
+true_quaternion_i2b = read_quaternion_from_csv(read_file_name, 'spacecraft_quaternion_i2b')
+
+# Statistics
+error_angle_rad = calc_error_angle_from_quaternions(measured_quaternion_i2b, true_quaternion_i2b)
+error_average = error_angle_rad.mean()
+standard_deviation = error_angle_rad.std()
+
+#
+# Plot
+#
+
+fig, axis = plt.subplots(4, 1, squeeze = False, tight_layout = True, sharex = True)
+axis[0, 0].plot(time[0], measured_quaternion_i2b[0], marker=".", c="red",    label="MEASURED-X")
+axis[0, 0].plot(time[0], true_quaternion_i2b[0], marker=".", c="orange",  label="TRUE-X")
+axis[0, 0].legend(loc = 'upper right')
+
+axis[1, 0].plot(time[0], measured_quaternion_i2b[1], marker=".", c="green",    label="MEASURED-Y")
+axis[1, 0].plot(time[0], true_quaternion_i2b[1], marker=".", c="yellow",  label="TRUE-Y")
+axis[1, 0].legend(loc = 'upper right')
+
+axis[2, 0].plot(time[0], measured_quaternion_i2b[2], marker=".", c="blue",    label="MEASURED-Z")
+axis[2, 0].plot(time[0], true_quaternion_i2b[2], marker=".", c="purple",  label="TRUE-Z")
+axis[2, 0].legend(loc = 'upper right')
+
+axis[3, 0].plot(time[0], measured_quaternion_i2b[3], marker=".", c="black",    label="MEASURED-W")
+axis[3, 0].plot(time[0], true_quaternion_i2b[3], marker=".", c="gray",  label="TRUE-W")
+axis[3, 0].legend(loc = 'upper right')
+
+fig.suptitle("Attitude Observer Quaternion")
+fig.supylabel("Quaternion")
+fig.supxlabel("Time [s]")
+
+unit = 'rad'
+plt.figure(0)
+plt.plot(time[0], error_angle_rad, marker=".", c="red")
+plt.title("Error angle \n" + "Error average:" + format(error_average, '+.2e') + unit + "\n Standard deviation:" + format(standard_deviation, '+.2e') + unit)
+plt.xlabel("Time [s]")
+plt.ylabel("Angle [rad]")
+
+# Data save
+if args.no_gui:
+  plt.savefig(read_file_tag + "_attitude_observer.png") # save last figure only
+else:
+  plt.show()