Adding data saving + plotting from data files

.gitignore

@@ -1,4 +1,5 @@
 build/
+data/
 .DS_Store
 .vscode/
 *.png

.vscode/settings.json

@@ -82,6 +82,7 @@
         "semaphore": "cpp",
         "typeindex": "cpp",
         "dense": "cpp",
-        "mprealsupport": "cpp"
+        "mprealsupport": "cpp",
+        "hash_map": "cpp"
     }
 }

example_project/CMakeLists.txt

@@ -19,4 +19,4 @@ target_link_libraries(run PRIVATE nlohmann_json::nlohmann_json Eigen3::Eigen ${S
 set_target_properties(run PROPERTIES
     BUILD_RPATH "${CMAKE_INSTALL_PREFIX}/lib"
     INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib"    
-    )
+    )

example_project/sim_parameters_example.json

@@ -1,6 +1,6 @@
 {
   "Simulation duration": 25000,
-  "epsilon": 0.000000000001,
+  "epsilon": 0.00000000000001,
   "Initial timestep guess": 1,
   "Perturbation": true,
   "Drag": false

example_project/simulation_setup.cpp

@@ -1,14 +1,15 @@
-#include <iostream>
-
 #include "PhasedArrayGroundStation.h"
 #include "Satellite.h"
 #include "utils.h"
 
+#include <iostream>
+
 int main() {
   // This file demonstrates a few different ways you can run and
   // visualize data from satellite simulations.
   // First let's initialize the struct containing parameters for simulation and
   // plotting
+
   SimParameters sim_parameters("../sim_parameters_example.json");
   // Initialize satellite object from an input JSON file
   Satellite test_sat_1("../example_satellite_input_files/input.json");
@@ -22,8 +23,8 @@ int main() {
                                      t_thrust_start, t_thrust_end);
 
   std::array<double, 3> LVLH_thrust_vec_2 = {0.51, 20, -5};
-  double t_thrust_2_start = 5500;
-  double t_thrust_2_end = 6500;
+  double t_thrust_2_start = 3800;
+  double t_thrust_2_end = 4500;
   test_sat_1.add_LVLH_thrust_profile(LVLH_thrust_vec_2, t_thrust_2_start,
                                      t_thrust_2_end);
 
@@ -33,7 +34,6 @@ int main() {
 
   std::vector<Satellite> satellite_vector_1 = {test_sat_1, test_sat_2,
                                                test_sat_3};
-
   sim_parameters.initial_timestep_guess = 2;
   sim_parameters.total_sim_time = 25000;
   sim_parameters.epsilon = pow(10, -12);
@@ -69,6 +69,21 @@ int main() {
   sim_and_plot_attitude_evolution_gnuplot(satellite_vector_3, sim_parameters,
                                           "Pitch", file_name);
 
+  // Equivalent workflow saving data
+  std::string pitch_datafile_prefix = "Sim1-";
+  sim_and_write_to_file(satellite_vector_3, sim_parameters, pitch_datafile_prefix);
+  std::string plotted_parameter = "Pitch";
+  std::string output_filename_2D = "Pitch_from_file";
+  std::vector<std::string> datafile_name_vector = {};
+  for (Satellite sat_object : satellite_vector_3) {
+        std::string sat_name = sat_object.get_name();
+        std::string datafile_name = pitch_datafile_prefix + sat_name;
+        datafile_name_vector.push_back(datafile_name);
+    }
+  plot_2D_from_datafile(datafile_name_vector,
+                            plotted_parameter,
+                            output_filename_2D);
+
   // Now let's demonstrate effect of atmospheric drag approximation
   Satellite test_sat_8("../example_satellite_input_files/input_8.json");
   Satellite test_sat_9("../example_satellite_input_files/input_9.json");
@@ -78,11 +93,11 @@ int main() {
   double F_10 = 100;  // Solar radio ten centimeter flux
   double A_p = 120;   // Geomagnetic A_p index
 
-  // Collect drag parameters into a tuple with F_10 first and A_p second
-  std::pair<double, double> drag_elements = {F_10, A_p};
   sim_parameters.total_sim_time = 10000;
   sim_parameters.epsilon = pow(10, -14);
   sim_parameters.drag_bool = true;
+  sim_parameters.F_10 = F_10;
+  sim_parameters.A_p = A_p;
   file_name = "Eccentricity Plot";
   sim_and_plot_orbital_elem_gnuplot(satellite_vector_4, sim_parameters,
                                     "Eccentricity", file_name);
@@ -157,13 +172,6 @@ int main() {
   sim_and_draw_orbit_gnuplot(orbit_transfer_demo_vec, sim_parameters);
 
   // Now let's demonstrate changing the argument of periapsis
-  // Calibration strategy: there are some inherent oscillations in, e.g., arg of
-  // periapsis, even in absence of external forces or perturbations This
-  // oscillation magnitude can be different in initial and final orbits, so
-  // current strategy is to find the mean val of these oscillations at initial
-  // and final orbits, use those to get an offset, which can be applied to the
-  // final argument of periapsis to try to make the oscillations at the final
-  // orbit be close to the target value
 
   Satellite arg_periapsis_change_sat(
       "../example_satellite_input_files/arg_of_periapsis_test_input.json");
@@ -176,8 +184,7 @@ int main() {
   double final_arg_of_periapsis_deg = 25;
 
   thrust_magnitude = 0.1;  // N
-  // std::cout << "orbital period: " <<
-  // arg_periapsis_change_sat.calculate_orbital_period() << "\n";
+
   arg_periapsis_change_sat.add_maneuver(
       "Argument of Periapsis Change", t_thrust_start,
       final_arg_of_periapsis_deg, thrust_magnitude);
@@ -188,5 +195,28 @@ int main() {
                                     sim_parameters, "Argument of Periapsis",
                                     file_name);
 
+  // Now let's demonstrate the workflow if you choose to save simulation data to a data file then plot from the data file(s)
+  sim_parameters.epsilon = pow(10,-12);
+  sim_parameters.total_sim_time = 25000;
+  sim_parameters.x_increment = pow(10,7);
+  sim_parameters.y_increment = pow(10,7);
+  sim_parameters.z_increment = 5*pow(10,6);
+  sim_parameters.drag_bool = false;
+  std::string datafile_prefix = "Sim2-";
+  sim_and_write_to_file(satellite_vector_1, sim_parameters, datafile_prefix);
+  datafile_name_vector = {};
+  for (Satellite sat_object : satellite_vector_1) {
+        std::string sat_name = sat_object.get_name();
+        std::string datafile_name = datafile_prefix + sat_name;
+        datafile_name_vector.push_back(datafile_name);
+    }
+  std::string output_file_name = "sample_output_3D_plot";
+  std::string terminal_name = "qt";
+  plot_3D_from_datafile(datafile_name_vector,
+                                output_file_name,
+                                terminal_name,
+                                sim_parameters.x_increment,
+                                sim_parameters.y_increment,
+                                sim_parameters.z_increment);
   return 0;
 }

include/Satellite.h

@@ -328,16 +328,7 @@ class Satellite {
     return sqrt(pow(ECI_position_.at(0), 2) + pow(ECI_position_.at(1), 2) +
                 pow(ECI_position_.at(2), 2));
   }
-  double get_total_energy() {
-    double orbital_radius = get_radius();
-    double gravitational_potential_energy =
-        -G * mass_Earth * m_ / orbital_radius;
-
-    double orbital_speed = get_speed();
-    double kinetic_energy = (1.0 / 2.0) * m_ * (orbital_speed * orbital_speed);
-
-    return (gravitational_potential_energy + kinetic_energy);
-  }
+  double get_total_energy();
   double get_mass() { return m_; }
   double get_instantaneous_time() { return t_; }
   std::string get_name() { return name_; }

include/utils.h

@@ -522,4 +522,20 @@ int add_lowthrust_orbit_transfer(Satellite& input_satellite_object,
                                  const double final_orbit_semimajor_axis_km,
                                  const double thrust_magnitude,
                                  const double transfer_initiation_time = 0);
+
+int sim_and_write_to_file(std::vector<Satellite> input_satellite_vector,
+                                const SimParameters& input_sim_parameters,
+                                const std::string output_file_name_prefix,
+                                const double size_limit = INFINITY);
+
+void plot_3D_from_datafile(std::vector<std::string> input_datafile_name_vector,
+                                const std::string output_file_name = "3D plot",
+                                const std::string terminal_name = "qt",
+                                const double x_increment = 0,
+                                const double y_increment = 0,
+                                const double z_increment = 0);
+
+void plot_2D_from_datafile(std::vector<std::string> input_datafile_name_vector,
+                            const std::string plotted_parameter,
+                            const std::string output_file_name = "2D plot");
 #endif

src/Satellite.cpp

@@ -764,4 +764,24 @@ void Satellite::add_maneuver(const std::string maneuver_type,
   std::pair<std::string, std::array<double, 3>> maneuver_info = {maneuver_type,
                                                                  maneuver_vals};
   maneuvers_awaiting_initiation_.push_back(maneuver_info);
+}
+
+double Satellite::get_total_energy() {
+  // First, gravitational potential energy
+  double orbital_radius = get_radius();
+  double gravitational_potential_energy =
+      -G * mass_Earth * m_ / orbital_radius;
+
+  // Next, translational kinetic energy
+  double orbital_speed = get_speed();
+  double translational_kinetic_energy = (1.0 / 2.0) * m_ * (orbital_speed * orbital_speed);
+
+  // Next, rotational kinetic energy
+  double rotational_kinetic_energy = 0;
+  std::array<double,3> body_angular_velocity_vec_wrt_LVLH_in_body_frame = body_angular_velocity_vec_wrt_LVLH_in_body_frame_;
+  rotational_kinetic_energy += (0.5 * J_11_ * body_angular_velocity_vec_wrt_LVLH_in_body_frame.at(0)*body_angular_velocity_vec_wrt_LVLH_in_body_frame.at(0));
+  rotational_kinetic_energy += (0.5 * J_22_ * body_angular_velocity_vec_wrt_LVLH_in_body_frame.at(1)*body_angular_velocity_vec_wrt_LVLH_in_body_frame.at(1));
+  rotational_kinetic_energy += (0.5 * J_33_ * body_angular_velocity_vec_wrt_LVLH_in_body_frame.at(2)*body_angular_velocity_vec_wrt_LVLH_in_body_frame.at(2));
+
+  return (gravitational_potential_energy + translational_kinetic_energy + rotational_kinetic_energy);
 }