From de208b1d3207554568815e3387f607758e8b509b Mon Sep 17 00:00:00 2001 From: Thomas West <tw00956@surrey.ac.uk> Date: Tue, 3 Sep 2024 15:23:56 +0000 Subject: [PATCH] Geometry and orbital coefficient calculations for an orbit. --- Function/geometry_calculations.m | 39 ++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 Function/geometry_calculations.m diff --git a/Function/geometry_calculations.m b/Function/geometry_calculations.m new file mode 100644 index 0000000..15786ab --- /dev/null +++ b/Function/geometry_calculations.m @@ -0,0 +1,39 @@ +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 -- GitLab