Add files via upload

src/orbit/orbit_fitting.cpp

@@ -0,0 +1,104 @@
+#include <iostream>
+#include <cmath>
+#include <functional>
+#include <tuple>
+
+using namespace std;
+
+double b_const = 1.0;
+double x = 0;
+double y = 1;
+double dt = 0.01;
+double v0 = abs(y-x)/dt;
+
+double a(double t,double y) {
+   return b_const*cos(t);
+}
+
+double a_prime(double t,double y) {
+   return -1*b_const*sin(t);
+}
+
+double a_double_prime(double t,double y) {
+   return -1*b_const*cos(t);
+}
+
+double b(double t, double y) {
+   return b_const*sin(t);
+}
+
+double b_prime(double t,double y) {
+   return b_const*cos(t);
+}
+
+double b_double_prime(double t,double y) {
+   return -1*b_const*sin(t);
+}
+
+double f(double t,double y,double v) {
+   return a_prime(t,y);
+}
+
+double g(double t,double y,double v) {
+   return a_double_prime(t,y);
+}
+
+double f1(double t,double y,double v) {
+   return b_prime(t,y);
+}
+
+double g1(double t,double y,double v) {
+   return b_double_prime(t,y);
+}
+
+std::tuple<std::vector<double>, std::vector<double>, std::vector<double>> rungeKutta( double t, double y, double dt, double n, double v, std::function<double(double, double, double)> f, std::function<double(double, double, double)> g ){
+
+   vector<double> y_vals = {y};
+   vector<double> v_vals = {v};
+   vector<double> t_vals = {t};;
+   for(int i = 1; i <= n; i++) {
+
+      double k1 = f(t, y_vals[i-1], v_vals[i-1]);
+      double l1 = g(t, y_vals[i-1], v_vals[i-1]);
+      double k2 = f(t + dt/2, y_vals[i-1] + dt * (k1 / 2), v_vals[i-1] + dt * (l1 / 2));
+      double l2 = g(t + dt/2, y_vals[i-1] + dt * (k1 / 2), v_vals[i-1] + dt * (l1 / 2));
+      double k3 = f(t + dt/2, y_vals[i-1] + dt * (k2 / 2), v_vals[i-1] + dt * (l2 / 2));
+      double l3 = g(t + dt/2, y_vals[i-1] + dt * (k2 / 2), v_vals[i-1] + dt * (l2 / 2));
+      double k4 = f(t + dt, y_vals[i-1] + dt * k3, v_vals[i-1] + dt * l3);
+      double l4 = g(t + dt, y_vals[i-1] + dt * k3, v_vals[i-1] + dt * l3);
+      double y_new = y_vals[i-1] + (1.0 / 6) * dt * (k1 + 2*k2 + 2*k3 + k4);
+      y_vals.push_back(y_new);
+      double v_new = v_vals[i-1] + (1.0 / 6) * dt * (l1 + 2*l2 + 2*l3 + l4);
+      v_vals.push_back(v_new);
+      t += dt;
+      t_vals.push_back(t);
+
+   }
+   return std::make_tuple(t_vals, y_vals, v_vals);
+}
+
+
+
+int main() {
+
+   auto [t_vals, y_vals, v_vals] = rungeKutta(0,1,dt,1000,v0,f,g);
+   auto [t1_vals, x_vals, v1_vals] = rungeKutta(0,1,dt,1000,v0,f1,g1);
+   auto [t2_vals, z_vals, v2_vals] = rungeKutta(0,1,dt,1000,v0,f,g);
+
+   std::cout << "\n\nX values:\n";
+   for(double val : x_vals) {
+      std::cout << val << " ";
+   }
+   std::cout << "\n\nV1 values:\n";
+   for(double val : v1_vals) {
+      std::cout << val << " ";
+   }
+   std::cout << "\n\nZ values:\n";
+   for(double val : z_vals) {
+      std::cout << val << " ";
+   }
+   std::cout << std::endl;
+
+
+   return 0;
+}

src/orbit/orbit_parameters.cpp

@@ -0,0 +1,94 @@
+#include <iostream>
+#include <cmath>
+#include <array>
+
+const double MU = 3.986e14;
+
+std::array<double, 3> cross_product(const std::array<double, 3>& a, const std::array<double, 3>& b) {
+    return {a[1] * b[2] - a[2] * b[1],
+            a[2] * b[0] - a[0] * b[2],
+            a[0] * b[1] - a[1] * b[0]};
+}
+
+double norm(const std::array<double, 3>& vec) {
+    return std::sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
+}
+
+std::array<double, 3> angular_momentum(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    return cross_product(position, velocity);
+}
+
+std::array<double, 3> eccentricity_vector(const std::array<double, 3>& position,
+                                          const std::array<double, 3>& velocity) {
+    auto h = angular_momentum(position, velocity);
+    double r_norm = norm(position);
+    double v_norm = norm(velocity);
+
+    std::array<double, 3> term1;
+    for (int i = 0; i < 3; ++i)
+        term1[i] = (v_norm * v_norm - MU / r_norm) * position[i] - (position[i] * velocity[i]) * velocity[i];
+
+    std::array<double, 3> cross_vh = cross_product(velocity, h);
+    std::array<double, 3> e_vec;
+    for (int i = 0; i < 3; ++i)
+        e_vec[i] = (1.0 / MU) * (cross_vh[i] - (MU / r_norm) * position[i]);
+
+    return e_vec;
+}
+
+double eccentricity(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    auto e_vec = eccentricity_vector(position, velocity);
+    return norm(e_vec);
+}
+
+
+double semi_major_axis(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    double r_norm = norm(position);
+    double v_norm = norm(velocity);
+    double energy = (v_norm * v_norm) / 2 - MU / r_norm;
+    return -MU / (2 * energy);
+}
+
+double apoapsis(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    double a = semi_major_axis(position, velocity);
+    double e = eccentricity(position, velocity);
+    return a * (1 + e);
+}
+
+double periapsis(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    double a = semi_major_axis(position, velocity);
+    double e = eccentricity(position, velocity);
+    return a * (1 - e);
+}
+
+double inclination(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    auto h = angular_momentum(position, velocity);
+    return std::acos(h[2] / norm(h));
+}
+
+double raan(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    auto h = angular_momentum(position, velocity);
+    std::array<double, 3> K = {0, 0, 1};
+    auto N_vec = cross_product(K, h);
+    double N = norm(N_vec);
+    return 2 * M_PI - std::acos(N_vec[0] / N);
+}
+
+double argument_of_periapsis(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    auto h = angular_momentum(position, velocity);
+    std::array<double, 3> K = {0, 0, 1};
+    auto N_vec = cross_product(K, h);
+    auto e_vec = eccentricity_vector(position, velocity);
+    double e = norm(e_vec);
+    double N = norm(N_vec);
+    double e = eccentricity(position, velocity);
+    return 2 * M_PI - std::acos((N_vec[0] * e_vec[0] + N_vec[1] * e_vec[1] + N_vec[2] * e_vec[2]) / (N * e));
+}
+
+double true_anomaly(const std::array<double, 3>& position, const std::array<double, 3>& velocity) {
+    double r = norm(position);
+    auto e_vec = eccentricity_vector(position, velocity);
+    double e = norm(e_vec);
+    double e = eccentricity(position, velocity);
+    return std::acos((position[0] * e_vec[0] + position[1] * e_vec[1] + position[2] * e_vec[2]) / (r * e));
+}