I need help with programming in MATLAB. I am trying to convert kepler elemnts to cartesian coordinates. The code works perfectly if my true anomaly only has one value. The problem occurs when true anomaly f is equal to 100 different points within 0 to 360 degrees. I get size errors. How do I get 100 different cartesian coordinates using the 100 f values.
I need help with programming in MATLAB. I am trying to convert kepler elemnts to cartesian coordinates. The code works perfectly if my true anomaly only has one value. The problem occurs when true anomaly f is equal to 100 different points within 0 to 360 degrees. I get size errors. How do I get 100 different cartesian coordinates using the 100 f values.
% Given parameters
omega_earth = rad2deg(7.2921151467e-5); % Earth's rotational rate in deg/s
period_of_repetition = 1 / 2; % Two orbits per day
% Orbital parameters
ecc = 0.74;
inc = 63.4349; % Inclination (deg)
raan = -86.915798; % RAAN (deg)
argp = 270; % Arg_of_Perigee (deg)
f = linspace(0, 360, 100); % True anomaly (deg)
mu = 398600.4418; % gravitational parameter of Earth (km^3/s^2)
t = 0:99;
% Calculate semi-major axis
a = ((omega_earth * period_of_repetition) / (360))^(-2/3) % km
% Calculate periapsis and apoapsis distances
periapsis_distance = a * (1 - ecc); % km
apoapsis_distance = a * (1 + ecc); % km
% Calculate orbital radius evolution with time
%r = (a * (1 - ecc^2) ./ (1 + ecc * cosd(f))) + 6731;
[x, y, z, vx, vy, vz] = kep2cart(a, ecc, inc, raan, argp, f)
function [x, y, z, vx, vy, vz] = kep2cart(a, ecc, inc, raan, argp, f)
%=========================================================================%
% FUNCTION: kep2cart
% DESCRIPTION: Converts kepler orbital elements to cartesian coordiantes
% INPUTS:
% > a = semi-major axis [km]
% > ecc = eccentricity
% > inc = inclination angle [degs]
% > raan = right ascension of the ascending node [degs]
% > argp = arguement of periapsis [degs]
% > f = true anomaly [degs]
%
% OUPUTS:
% > x = 1x1 position value in the x-direction [km]
% > y = 1x1 position value in the y-direction [km]
% > z = 1x1 position value in the z-direction [km]
% > vx = 1x1 velocity value in the x-direction [km/s]
% > vy = 1x1 velocity value in the y-direction [km/s]
% > vz = 1x1 velocity value in the z-direction [km/s]
%=========================================================================%
% Gravitational parameter
mu = 398600.4418; % (km^3/s^2)
% Calculate angular momentum
h = sqrt(a * (1-ecc^2)*mu); % (km^2/s)
% Calculate position and velocity in the periforcal frame
r_w = ((h^2 / mu) / (1 + ecc * cosd(f))) .* [cosd(f), sind(f), 0];
v_w = (mu / h) .* [-sind(f), ecc + cosd(f), 0];
% Calculate the Rotational Matrix
R1 = [cosd(-argp) -sind(-argp) 0;
sind(-argp) cosd(-argp) 0;
0 0 1];
R2 = [1 0 0;
0 cosd(-inc) -sind(-inc);
0 sind(-inc) cosd(-inc)];
R3 = [cosd(-raan) -sind(-raan) 0;
sind(-raan) cosd(-raan) 0;
0 0 1];
% Calculate position
r_rot = r_w * R1 * R2 * R3;
% Calculate velocity vector
v_rot = v_w * R1 * R2 * R3;
% Define the cartesian coordinates
x = r_rot(1);
y = r_rot(2);
z = r_rot(3);
vx = v_rot(1);
vy = v_rot(2);
vz = v_rot(3);
end
Step by step
Solved in 4 steps with 4 images
