Geometry and orbital coefficient calculations for an orbit.

Function/geometry_calculations.m

+function [inclination, RAAN, argument_of_periapsis, eccentricity] = geometry_calculations(r, v, mu)
+    % Calculate the angular momentum vector
+    h = cross(r, v);
+    n = cross([0, 0, 1], h);
+    n_norm = norm(n);
+    h_norm = norm(h);
+    
+    % Calculate the inclination
+    inclination = acos(h(3) / h_norm);
+    inclination = rad2deg(inclination);  % Convert from radians to degrees
+    
+    % Calculate the Right Ascension of the Ascending Node (RAAN)
+    if n_norm ~= 0
+        RAAN = atan2(n(2), n(1));
+        RAAN = rad2deg(RAAN);  % Convert from radians to degrees
+        if RAAN < 0
+            RAAN = RAAN + 360;  % Ensure the angle is positive
+        end
+    else
+        RAAN = 0;  % RAAN is undefined for zero inclination
+    end
+    
+    % Calculate the eccentricity vector
+    r_norm = norm(r);
+    v_norm = norm(v);
+    e = (1/mu) * ((v_norm^2 - mu/r_norm) * r - dot(r, v) * v);
+    eccentricity = norm(e);  % Eccentricity is the magnitude of the eccentricity vector
+    
+    % Calculate the Argument of Periapsis
+    if n_norm ~= 0 && eccentricity > 0
+        argument_of_periapsis = acos(dot(n, e) / (n_norm * norm(e)));
+        if e(3) < 0
+            argument_of_periapsis = 2 * pi - argument_of_periapsis;  % Correct the angle if e_z is negative
+        end
+        argument_of_periapsis = rad2deg(argument_of_periapsis);  % Convert from radians to degrees
+    else
+        argument_of_periapsis = 0;  % Argument of periapsis is undefined if no inclination or eccentricity
+    end
+end