diff --git a/Function/geometry_calculations.m b/Function/geometry_calculations.m new file mode 100644 index 0000000000000000000000000000000000000000..15786ab4e3b88f6506e4cc2883152cdefd513e88 --- /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