Refactoring part of argument of periapsis change maneuver workflow, rerunning clang-format, updating test coverage reports

.gitignore

@@ -1,4 +1,8 @@
 build/
 .DS_Store
 .vscode/
-*.png
+*.png
+src/files_to_format.txt
+include/files_to_format.txt
+tests/files_to_format.txt
+files_to_format.txt

README.md

@@ -32,7 +32,7 @@
 
    - Continuous (low) thrust maneuvers to change argument of perigee
 
-      - Note: Thrust profiles of argument of periapsis change maneuvers are calculated assuming a "continuous" maneuver takes place (i.e., a burn arc angle $\alpha = \pi/2$ radians in https://apps.dtic.mil/sti/tr/pdf/ADA384536.pdf). Therefore it is not recommended to use this feature outside of parameter ranges where this assumption does not hold.
+      - Note: Thrust profiles of argument of periapsis change maneuvers are calculated assuming a "continuous" maneuver takes place (i.e., a burn arc angle $\alpha = \pi/2$ radians in https://apps.dtic.mil/sti/tr/pdf/ADA384536.pdf). Therefore it is not recommended to use this feature outside of parameter ranges where this assumption does not hold. Maneuvers that take place over a fraction of an orbit are not recommended.
 
 
 ## TO DO
@@ -72,6 +72,8 @@ Simulation of zero-inclination orbits isn't currently supported, as the magnitud
 
 Calculation of instantaneous orbital angular acceleration does not currently include any contribution from the time derivative of the magnitude of the orbital angular momentum vector.
 
+Loss of propellant mass during burns is not yet accounted for.
+
 Visualization/plotting is done via Gnuplot. The copyright and permission notice of Gnuplot is shown below:
 
 Copyright 1986 - 1993, 1998, 2004   Thomas Williams, Colin Kelley

example_satellite_input_files/arg_of_periapsis_calibration_input.json

@@ -1,7 +1,7 @@
 {
   "Inclination": 5,
   "RAAN": 0,
-  "Argument of Periapsis": 20,
+  "Argument of Periapsis": 30,
   "Eccentricity": 0.1,
   "Semimajor Axis": 15000,
   "True Anomaly": 0,

include/PhasedArrayGroundStation.h

@@ -1,41 +1,49 @@
 #ifndef GROUNDSTATION_HEADER
 #define GROUNDSTATION_HEADER
 
-#include <iostream>
 #include <Eigen/Dense>
+#include <iostream>
 #include <stdexcept>
+
 #include "Satellite.h"
 
 using Eigen::Vector3d;
 
 class PhasedArrayGroundStation {
-    // Phased array ground station class
-    public:
-      double latitude = {0};
-      double longitude = {0};
-      double height = {0};
-      Vector3d ECEF_position_;
-      double max_beam_angle_from_normal_ = {0}; //deg
-      int num_beams_ = {0};
-      // To enforce the maximum number of satellites a ground station can communicate with at any one time,
-      // going to create a map where each key is the satellite index of a given satellite (i.e., index in a satellite vector being plotted)
-      // Each key has a val which is a vector of pairs, where the first element is the start of a time range where the
-      // ground station is communicating with that satellite, and the second element is the end time for that range
-      // So communications with satellites are characterized by vectors of time ranges
+  // Phased array ground station class
+ public:
+  double latitude = {0};
+  double longitude = {0};
+  double height = {0};
+  Vector3d ECEF_position_;
+  double max_beam_angle_from_normal_ = {0};  // deg
+  int num_beams_ = {0};
+  // To enforce the maximum number of satellites a ground station can
+  // communicate with at any one time, going to create a map where each key is
+  // the satellite index of a given satellite (i.e., index in a satellite vector
+  // being plotted) Each key has a val which is a vector of pairs, where the
+  // first element is the start of a time range where the ground station is
+  // communicating with that satellite, and the second element is the end time
+  // for that range So communications with satellites are characterized by
+  // vectors of time ranges
+
+  std::map<size_t, std::vector<std::pair<double, double>>> linked_sats_map_;
 
-      std::map<size_t,std::vector<std::pair<double,double>>> linked_sats_map_;
-
-      PhasedArrayGroundStation(double latitude, double longitude, double height, double max_beam_angle_from_normal, int num_beams = 1) {
-        generate_ECEF_position();
-        num_beams_ = num_beams;
-        max_beam_angle_from_normal_ = max_beam_angle_from_normal;
-      }
-      void generate_ECEF_position();
-      Vector3d get_ECI_position(double input_time);
-      double distance_to_satellite(Satellite input_satellite);
-      double angle_to_satellite_from_normal(Satellite input_satellite);
-      int num_sats_connected_at_this_time(double input_time);
-      void update_linked_sats_map(const size_t satellite_index, const double connection_time, const double previous_time);
-  };
+  PhasedArrayGroundStation(double latitude, double longitude, double height,
+                           double max_beam_angle_from_normal,
+                           int num_beams = 1) {
+    generate_ECEF_position();
+    num_beams_ = num_beams;
+    max_beam_angle_from_normal_ = max_beam_angle_from_normal;
+  }
+  void generate_ECEF_position();
+  Vector3d get_ECI_position(double input_time);
+  double distance_to_satellite(Satellite input_satellite);
+  double angle_to_satellite_from_normal(Satellite input_satellite);
+  int num_sats_connected_at_this_time(double input_time);
+  void update_linked_sats_map(const size_t satellite_index,
+                              const double connection_time,
+                              const double previous_time);
+};
 
 #endif

include/Satellite.h

@@ -19,13 +19,16 @@ using json = nlohmann::json;
 
 class ThrustProfileLVLH {
   // Note: for now, thrust forces are assumed to act through center of mass of
-  // satellite.
+  // the satellite.
  public:
   double t_start_ = {0};
   double t_end_ = {0};
   bool arg_of_periapsis_change_thrust_profile = false;
+  // =======================================================================
+  // Only used for argument of periapsis change thrust profiles
   double sign_of_delta_omega = {1};
-  double thrust_magnitude_ = {0}; // Only used for argument of periapsis change thrust profiles
+  double thrust_magnitude_ = {0};  
+  // =======================================================================
   std::array<double, 3> LVLH_force_vec_ = {0, 0, 0};
   ThrustProfileLVLH(const double t_start, const double t_end,
                     const std::array<double, 3> LVLH_force_vec) {
@@ -45,24 +48,33 @@ class ThrustProfileLVLH {
     }
   }
   ThrustProfileLVLH(const double t_start, const double final_arg_of_periapsis,
-                  const  double thrust_magnitude, const double initial_arg_of_periapsis,
-                  const double satellite_eccentricity, const double satellite_a,
-                  const double satellite_mass) {
+                    const double thrust_magnitude,
+                    const double initial_arg_of_periapsis,
+                    const double satellite_eccentricity,
+                    const double satellite_a, const double satellite_mass) {
     // Thrust profile for continuous-thrust argument of periapsis change
     // Main ref: https://apps.dtic.mil/sti/tr/pdf/ADA384536.pdf
-    // 
-    // Used Eq. 67 in https://link.springer.com/article/10.1007/s10569-021-10033-9#Sec15 to determine burn angles
+    // Used Eq. 67 in
+    // https://link.springer.com/article/10.1007/s10569-021-10033-9#Sec15 to
+    // determine burn angles Arguments of periapsis are input in radians
     arg_of_periapsis_change_thrust_profile = true;
     t_start_ = t_start;
-    const double alpha = M_PI/2.0; // continous thrust
-    (final_arg_of_periapsis > initial_arg_of_periapsis) ? sign_of_delta_omega = 1.0 : sign_of_delta_omega = -1.0;
+    const double alpha = M_PI / 2.0;  // continous thrust
+    (final_arg_of_periapsis > initial_arg_of_periapsis)
+        ? sign_of_delta_omega = 1.0
+        : sign_of_delta_omega = -1.0;
     thrust_magnitude_ = thrust_magnitude;
-    const double mu_Earth = G*mass_Earth;
-    const double delta_omega_mag = abs(final_arg_of_periapsis - initial_arg_of_periapsis); // = delta_omega / sgn(delta_omega)
-    const double delta_V = (2.0/3.0)*sqrt(mu_Earth/satellite_a)*satellite_eccentricity*delta_omega_mag/sqrt(1-satellite_eccentricity*satellite_eccentricity);
-    const double acceleration = thrust_magnitude/satellite_mass;
-    t_end_ = t_start + delta_V/acceleration;
-    std::cout << "Transfer time: " << t_end_ << "\n";
+    const double mu_Earth = G * mass_Earth;
+    const double delta_omega_mag =
+        abs(final_arg_of_periapsis -
+            initial_arg_of_periapsis);  // = delta_omega / sgn(delta_omega)
+    const double delta_V =
+        (2.0 / 3.0) * sqrt(mu_Earth / satellite_a) * satellite_eccentricity *
+        delta_omega_mag /
+        sqrt(1 - satellite_eccentricity * satellite_eccentricity);
+    const double acceleration = thrust_magnitude / satellite_mass;
+    t_end_ = t_start + delta_V / acceleration;
+    std::cout << "Maneuver duration: " << t_end_ << "\n";
   }
   bool operator==(const ThrustProfileLVLH& input_profile) {
     return ((t_start_ == input_profile.t_start_) &&
@@ -106,23 +118,24 @@ class BodyframeTorqueProfile {
 class Satellite {
  private:
   double inclination_ = {0};
-  double raan_ = {
-      0};  // Assuming RAAN can be used interchangeably with longitude of
-           // ascending node for the Earth-orbiting satellites simulated here
+  double raan_ = {0};  
+  // Assuming RAAN can be used interchangeably with longitude of
+  // ascending node for the Earth-orbiting satellites simulated here
   double arg_of_periapsis_ = {0};
   double eccentricity_ = {0};
   double a_ = {0};
   double true_anomaly_ = {0};
   double orbital_period_ = {0};
-  double m_ = {1};  // default value to prevent infinities in acceleration
-                    // calculations from a=F/m
+  double m_ = {1};  
+  // default value to prevent infinities in acceleration
+  // calculations from a=F/m
   // double I_={1}; //moment of inertia, taken to be same for all 3 principal
   // axes, set to default value for same reasons as mass
   double t_ = {0};
 
   double orbital_rate_ = {0};
-  double orbital_angular_acceleration_ = {0};  // Time derivative of orbital
-                                               // rate
+  // Time derivative of orbital rate
+  double orbital_angular_acceleration_ = {0};  
 
   // Now body-frame attributes
   // Assuming diagonal J matrix
@@ -141,8 +154,8 @@ class Satellite {
 
   // body-frame angular velocities relative to the LVLH frame, represented in
   // the body frame
-  std::array<double, 3> body_angular_velocity_vec_wrt_LVLH_in_body_frame_ = {
-      0, 0, 0};
+  std::array<double, 3> body_angular_velocity_vec_wrt_LVLH_in_body_frame_ 
+    = {0, 0, 0};
 
   // quaternion representing attitude of satellite body frame with respect to
   // the LVLH frame
@@ -163,7 +176,8 @@ class Satellite {
   std::vector<std::array<double, 3>> list_of_ECI_forces_at_this_time_ = {};
   std::vector<std::array<double, 3>> list_of_body_frame_torques_at_this_time_ =
       {};
-
+  std::vector<std::pair<std::string, std::array<double, 3>>>
+      maneuvers_awaiting_initiation_ = {};
   double drag_surface_area = {0};  // Surface area of satellite used for
   // atmospheric drag calculations
 
@@ -354,8 +368,9 @@ class Satellite {
   void add_LVLH_thrust_profile(
       const std::array<double, 3> input_LVLH_thrust_vector,
       const double input_thrust_start_time, const double input_thrust_end_time);
-  void add_LVLH_thrust_profile(const double input_thrust_start_time, const double final_arg_of_periapsis,
-                                const double input_thrust_magnitude);
+  void add_LVLH_thrust_profile(const double input_thrust_start_time,
+                               const double final_arg_of_periapsis,
+                               const double input_thrust_magnitude);
   void add_bodyframe_torque_profile(
       const std::array<double, 3> input_bodyframe_direction_unit_vec,
       const double input_bodyframe_torque_magnitude,
@@ -381,6 +396,11 @@ class Satellite {
   double get_attitude_val(const std::string input_attitude_val_name);
   double calculate_orbital_period();
   int add_arg_of_periapsis_change_thrust();
+  void check_for_maneuvers_to_initialize();
+  void add_maneuver(const std::string maneuver_type,
+                    const double maneuver_start_time,
+                    const double final_parameter_val,
+                    const double thrust_magnitude);
 };
 
 #endif

include/utils.h

@@ -6,8 +6,8 @@
 #include <iostream>
 #include <nlohmann/json.hpp>
 
-#include "Satellite.h"
 #include "PhasedArrayGroundStation.h"
+#include "Satellite.h"
 
 using Eigen::Matrix3d;
 using Eigen::MatrixXd;
@@ -16,7 +16,6 @@ using Eigen::Vector4d;
 
 using json = nlohmann::json;
 
-
 struct SimParameters {
   double initial_timestep_guess = 1;
   double total_sim_time = 10;
@@ -88,9 +87,11 @@ std::array<double, 6> RK4_deriv_function_orbit_position_and_velocity(
 // template <int T>
 // std::array<double, T> RK4_step(
 //     const std::array<double, T> y_n, const double input_step_size,
-//     std::function<std::array<double, T>(const std::array<double, T> input_y_vec,
+//     std::function<std::array<double, T>(const std::array<double, T>
+//     input_y_vec,
 //                                         const double input_spacecraft_mass,
-//                                         const std::vector<std::array<double, 3>>
+//                                         const std::vector<std::array<double,
+//                                         3>>
 //                                             input_vec_of_force_vectors_in_ECI)>
 //         input_derivative_function,
 //     const double input_spacecraft_mass,
@@ -101,7 +102,8 @@ std::array<double, 6> RK4_deriv_function_orbit_position_and_velocity(
 //     const std::vector<std::array<double, 3>>
 //         input_vec_of_force_vectors_in_ECI_at_t_and_step = {}) {
 //   // ref:
-//   // https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method
+//   //
+//   https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method
 //   std::array<double, T> y_nplus1 = y_n;
 
 //   // first, k=1;
@@ -148,11 +150,9 @@ std::array<double, 6> RK4_deriv_function_orbit_position_and_velocity(
 //   return y_nplus1;
 // }
 
-void sim_and_draw_orbit_gnuplot(
-    std::vector<Satellite> input_satellite_vector, 
-    const SimParameters& input_sim_parameters,
-    const std::string output_file_name = "output"
-);
+void sim_and_draw_orbit_gnuplot(std::vector<Satellite> input_satellite_vector,
+                                const SimParameters& input_sim_parameters,
+                                const std::string output_file_name = "output");
 
 template <int T>
 std::pair<std::array<double, T>, std::pair<double, double>> RK45_step(
@@ -166,6 +166,8 @@ std::pair<std::array<double, T>, std::pair<double, double>> RK45_step(
     double input_t_n, double input_epsilon, double input_inclination,
     double input_arg_of_periapsis, double input_true_anomaly,
     bool perturbation) {
+  // std::cout << "In RK45_step, received perturbation bool: " << perturbation
+  // << "\n";
   // Version for satellite orbital motion time evolution
   // Implementing RK4(5) method for its adaptive step size
   // Refs:https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
@@ -205,8 +207,8 @@ std::pair<std::array<double, T>, std::pair<double, double>> RK45_step(
     std::array<double, T> k_vec_at_this_s;
     for (size_t s_ind = 0; s_ind < k_ind; s_ind++) {
       for (size_t y_val_ind = 0; y_val_ind < y_n.size(); y_val_ind++) {
-        y_n_evaluated_value.at(y_val_ind) += RK_matrix(k_ind, s_ind) *
-                                             k_vec_vec.at(s_ind).at(y_val_ind);
+        y_n_evaluated_value.at(y_val_ind) +=
+            RK_matrix(k_ind, s_ind) * k_vec_vec.at(s_ind).at(y_val_ind);
       }
     }
     std::array<double, 6> derivative_function_output =
@@ -290,14 +292,14 @@ std::array<double, 3> convert_cylindrical_to_cartesian(
     const double input_r_comp, const double input_theta_comp,
     const double input_z_comp, const double input_theta);
 void sim_and_plot_orbital_elem_gnuplot(
-    std::vector<Satellite> input_satellite_vector, 
+    std::vector<Satellite> input_satellite_vector,
     const SimParameters& input_sim_parameters,
     const std::string input_orbital_element_name,
     const std::string file_name = "output");
 void sim_and_plot_attitude_evolution_gnuplot(
     std::vector<Satellite> input_satellite_vector,
     const SimParameters& input_sim_parameters,
-    const std::string input_plotted_val_name, 
+    const std::string input_plotted_val_name,
     const std::string file_name = "output");
 
 Matrix3d rollyawpitch_bodyframe_to_LVLH(
@@ -499,22 +501,25 @@ std::array<double, 3> convert_array_from_LVLH_to_bodyframe(
     const std::array<double, 3> input_LVLH_frame_array, const double input_roll,
     const double input_yaw, const double input_pitch);
 
+Vector3d convert_lat_long_to_ECEF(const double latitude, const double longitude,
+                                  const double height);
+Vector3d convert_ECEF_to_ECI(const Vector3d input_ECEF_position,
+                             const double input_time);
 
-Vector3d  convert_lat_long_to_ECEF(const double latitude, const double longitude, const double height);
-Vector3d  convert_ECEF_to_ECI(const Vector3d input_ECEF_position, const double input_time);
-
-void sim_and_plot_gs_connectivity_distance_gnuplot(PhasedArrayGroundStation input_ground_station,
-  std::vector<Satellite> input_satellite_vector, 
-  const SimParameters& input_sim_parameters,
-  const std::string file_name = "output");
-
-void sim_and_plot_gs_connectivity_gnuplot(PhasedArrayGroundStation input_ground_station,
-  std::vector<Satellite> input_satellite_vector, 
-  const SimParameters& input_sim_parameters,
-  const std::string file_name = "output");
+void sim_and_plot_gs_connectivity_distance_gnuplot(
+    PhasedArrayGroundStation input_ground_station,
+    std::vector<Satellite> input_satellite_vector,
+    const SimParameters& input_sim_parameters,
+    const std::string file_name = "output");
 
-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);
+void sim_and_plot_gs_connectivity_gnuplot(
+    PhasedArrayGroundStation input_ground_station,
+    std::vector<Satellite> input_satellite_vector,
+    const SimParameters& input_sim_parameters,
+    const std::string file_name = "output");
 
-double calibrate_mean_val(Satellite satellite_object, const SimParameters& input_sim_parameters, const std::string input_parameter_name);
+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);
 #endif