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 vector
    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