Plotting orbit in 3D

100 ビュー (過去 30 日間)
imran shaikh
imran shaikh 2021 年 10 月 9 日
コメント済み: imran shaikh 2021 年 10 月 9 日
%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

回答 (1 件)

Walter Roberson
Walter Roberson 2021 年 10 月 9 日
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
  1 件のコメント
imran shaikh
imran shaikh 2021 年 10 月 9 日
Thanks Walter, I will review the codes. appreciate it

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by