Plotting orbit in 3D

%Question:  for the given initail state vectors: Plot the magnitude of the position vector versus time [x axis: Time (hours) – y axis: radius (km)] 
% Plot the magnitude of the velocity vector versus time [x axis: Time (hours) – y axis: velocity (km/s)] 
% Plot the orbit in 3D including a suitably dimensioned sphere for Earth.

% Issues: 1) 3D orbit is missing sphere. how do i get the spherre for earth to display?
% 2) Its an instinct that both X Y plot are not displayed correctly, can someone help verify?

function [f] = twobody (t,X)
% Designed to call two body eqautions of motion
% Author: Imran Shaikh, ISHA032
%x(1)= x position;
%x(2)= y position;
%x(3) = z position;
%x (4) = x velocity;
%x (5) = y velocity;
%x (6) = z velocity;
r = 6378.14 *1e3;
a = 15000 *1e3;             % in meters
e = 0.15;
mu = 398600;                % km^3/s^2
f = zeros (size(X));
Xdot(1:3) = X (4:6);
r=norm(X(1:3));
Xdot(4:6)= (-mu/r^3)*X(1:3);
end



%Script below to call the function


% Define initial positions
r = [6510.75986901532, 2676.16546382759, 333.334029373109]*1e3;  % in meters. Transpose form
v = [-2.23202460862428, 9.49860960555864, 1.18311436869621]*1e3; % m/s. Transpose form

options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12);              % stated options to improve efficiency

% Define time vector for one day with a step of 10 seconds
tspan = 0:10:24*3600;

% Call ode
x0 = [r, v];                                                     % Column vector form
[t, X] = ode45(@(t, X)twobody(t, X), tspan, x0, options);

x = X(:,1);
y = X(:,2);
z = X(:,3);
vx = X(:,4);
vy = X(:,5);
vz = X(:,6);

% Calculate radius
r_orbit = sqrt(x.^2 + y.^2 + z.^2);

% Plot results
figure
subplot(1,2,1)
plot(tspan/3600, r_orbit/1e3)
xlabel('Time (hours)')
ylabel('Radius (km)')
grid on
title('Radius vs. Time')

% Calculate velocity
v_orbit = sqrt(vx.^2 + vy.^2 + vz.^2);

subplot(1,2,2)
plot(tspan/3600, v_orbit/1e3)
xlabel('Time (hours)')
ylabel('Velocity (km/s)')
grid on
title('Velocity vs. Time')

% Plot 3D orbit
figure
plot3(x/1e3, y/1e3, z/1e3), grid on
xlabel('X (km)')
ylabel('Y (km)')
zlabel('Z (km)')

title('Orbit')

%Results obtained

Answer 1

[Walter Roberson]

You were initializing f = zeros(size(X)) and then assigning derivatives information into Xdot instead of into f.

The below is obviously still not right -- but why does your two-body function have a and e assigned but not use them?

I do not see anything in your twobody() function that would lead towards a curve.

% Define initial positions
r = [6510.75986901532, 2676.16546382759, 333.334029373109]*1e3;  % in meters. Transpose form
v = [-2.23202460862428, 9.49860960555864, 1.18311436869621]*1e3; % m/s. Transpose form

options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12);              % stated options to improve efficiency

% Define time vector for one day with a step of 10 seconds
tspan = 0:10:24*3600;

% Call ode
x0 = [r, v];                                                     % Column vector form
[t, X] = ode45(@(t, X)twobody(t, X), tspan, x0, options);

x = X(:,1);
y = X(:,2);
z = X(:,3);
vx = X(:,4);
vy = X(:,5);
vz = X(:,6);

% Calculate radius
r_orbit = sqrt(x.^2 + y.^2 + z.^2);

% Plot results
figure
subplot(1,2,1)
plot(tspan/3600, r_orbit/1e3)
xlabel('Time (hours)')
ylabel('Radius (km)')
grid on
title('Radius vs. Time')

% Calculate velocity
v_orbit = sqrt(vx.^2 + vy.^2 + vz.^2);

subplot(1,2,2)
plot(tspan/3600, v_orbit/1e3)
xlabel('Time (hours)')
ylabel('Velocity (km/s)')
grid on
title('Velocity vs. Time')

% Plot 3D orbit
figure
plot3(x/1e3, y/1e3, z/1e3), grid on
xlabel('X (km)')
ylabel('Y (km)')
zlabel('Z (km)')

title('Orbit')

min(x), max(x)

ans = -1.8634e+08

ans = 6.5108e+06

min(y), max(y)

ans = 2.6762e+06

ans = 8.2336e+08

min(z), max(z)

ans = 3.3333e+05

ans = 1.0255e+08

function Xdot = twobody (t,X)
% Designed to call two body eqautions of motion
% Author: Imran Shaikh, ISHA032
%x(1)= x position;
%x(2)= y position;
%x(3) = z position;
%x (4) = x velocity;
%x (5) = y velocity;
%x (6) = z velocity;
r = 6378.14 *1e3;
a = 15000 *1e3;             % in meters
e = 0.15;
mu = 398600;                % km^3/s^2
Xdot = zeros (size(X));
Xdot(1:3) = X (4:6);
r=norm(X(1:3));
Xdot(4:6)= (-mu/r^3)*X(1:3);
end

[imran shaikh]

Thanks Walter, I will review the codes. appreciate it