I need help with MATLAB programmin. I am having a hard time with plotting in MATLAB. The following outputs a plot of the earth with a orbit inside it. The problem is that the calculations in the orbit do not include the radius of the earth. How do I increase the radius of the orbit by the radius of the earth to make the orbit outside of the earth. I tried adding the earth's radius to the x,y,z coordinates, but that only changes the place of the orbit in the plot not the radius.

I need help with MATLAB programmin. I am having a hard time with plotting in MATLAB. The following outputs a plot of the earth with a orbit inside it. The problem is that the calculations in the orbit do not include the radius of the earth. How do I increase the radius of the orbit by the radius of the earth to make the orbit outside of the earth. I tried adding the earth's radius to the x,y,z coordinates, but that only changes the place of the orbit in the plot not the radius.

% Given parameters
omega_earth = rad2deg(7.2921151467e-5);
period_of_repetition = 1 / 2;
ecc = 0.74;
inc = 63.4349;
raan = -86.915798;
argp = 270;
f = linspace(0, 360, 100);
mu = 398600.4418;
t = 0:99;

% Calculate semi-major axis
a = ((omega_earth * period_of_repetition) / (360))^(-2/3);

% Calculate periapsis and apoapsis distances
periapsis_distance = a * (1 - ecc);
apoapsis_distance = a * (1 + ecc);

% Initialize arrays to store Cartesian coordinates and velocities
x = zeros(1, 100);
y = zeros(1, 100);
z = zeros(1, 100);
vx = zeros(1, 100);
vy = zeros(1, 100);
vz = zeros(1, 100);

% Loop through each true anomaly
for i = 1:length(f)
    % Call kep2cart for each true anomaly
    [x(i), y(i), z(i), vx(i), vy(i), vz(i)] = kep2cart(a, ecc, inc, raan, argp, f(i));
end

% Plot 3D orbit
figure;
plot3(x, y, z, 'LineWidth', 2, 'Color', 'r');
hold on;

% Plot the Earth as a simple sphere
earth_radius = 6.371e3; % Earth's radius in meters
[x_earth, y_earth, z_earth] = sphere;
earth = surf(x_earth * earth_radius, y_earth * earth_radius, z_earth * earth_radius);
colormap([0 0 1]); % Set the color to blue for Earth
alpha(earth, 0.7); % Set transparency for Earth

% Set equal axis scale
axis equal;

% Set labels and title
xlabel('X (meters)');
ylabel('Y (meters)');
zlabel('Z (meters)');
title('3D Representation of the Orbit around Earth');

% Add grid and legend
grid on;
legend('Orbit', 'Earth');

% Display the plot
hold off;

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