JayFlynn
- 7
- 0
I am trying to plot the trajectory of an asteroid in MATLAB using ode23. The only bodies in the system are The Sun, Earth, Mars and Jupiter and their orbital data has been loaded from data files. I have picked arbitrary initial conditions for the asteroid and believe my forces are correct. My only issue is getting the ode23 to work. Any idea why the ode solver doesn't work?
<Moderator's note: please use code tags>
Matlab:
function SEMJ_5
global T X_J Y_J Z_J X_M Y_M Z_M X_E Y_E Z_E mu_E mu_M mu_J mu_S
G = 6.674e-11 * 1e+9 * 149.6e6; % Gravitational Constant (km^3.s^-2.kg^-1)
M_E = 5.972e+24; % Mass of the Earth (kg)
M_M = 0.64171e+24; % Mass of Mars (kg)
M_J = 1898.19e+24; % Mass of Jupiter (kg)
M_S = 1988500e+24; % Mass of Jupiter (kg)
mu_E = G*M_E;
mu_M = G*M_M;
mu_J = G*M_J;
mu_S = G*M_S;
% Load astronomical data
A = dlmread('1996-2017Earth.txt');
B = dlmread('1996-2017Mars.txt');
C = dlmread('1996-2017Jupiter.txt');
T = 0:(size(A,1)-1);
% Convert data into cartesian coordiantes for Earth
R_E = A(:,3);
theta_E = A(:,4);
Phi_E = A(:,5);
X_E = R_E.*cosd(theta_E).*cosd(Phi_E)*149.6e6;
Y_E = R_E.*cosd(theta_E).*sind(Phi_E)*149.6e6;
Z_E = R_E.*sind(theta_E)*149.6e6;
% Convert data into cartesian coordiantes for Mars
R_M = B(:,3);
theta_M = B(:,4);
Phi_M = B(:,5);
X_M = R_M.*cosd(theta_M).*cosd(Phi_M)*149.6e6;
Y_M = R_M.*cosd(theta_M).*sind(Phi_M)*149.6e6;
Z_M = R_M.*sind(theta_M)*149.6e6;
% Convert data into cartesian coordiantes for Jupiter
R_J= C(:,3);
theta_J = C(:,4);
Phi_J = C(:,5);
X_J = R_J.*cosd(theta_J).*cosd(Phi_J)*149.6e6;
Y_J = R_J.*cosd(theta_J).*sind(Phi_J)*149.6e6;
Z_J = R_J.*sind(theta_J)*149.6e6;
% Define initial conditions of the Asteroid
rX0 = 149.6e6; % Initial x coordinate
rY0 = 149.6e6; % Initial y coordinate
rZ0 = 149.6e6; % Initial z coordinate
V_AX0 = 10000 % Initial Velocity in x
V_AY0 = 10000 % Initial Velocity in y
V_AZ0 = 10000 % Initial Velocity in z
r_vecIC = [rX0 rY0 rZ0]; % Initial possition vector
Vel_vecIC = [V_AX0 V_AY0 V_AZ0]; % Initial velocity vector
init_val = [r_vecIC Vel_vecIC]; %Initial values for ode23
%% Numerical Integration
% ode23
[T,RV] = ode23(@(t,RV) Gravity(t, RV),T', init_val);
% :,_ describes what row/collum that value relates to
x_1 = RV(:,1);
y_1 = RV(:,2);
z_1 = RV(:,3);
u_1 = RV(:,4);
v_1 = RV(:,5);
w_1 = RV(:,6);%%
% Plot orbits
% figure(1)
% p0 = plot3(x_1 y_1 z_1);
figure(1)
p1 = plot3(X_E,Y_E,Z_E);
figure(1)
p2 = plot3(X_M,Y_M,Z_M);
figure(1)
p3 = plot3(X_J,Y_J,Z_J);
legend([p1 p2 p3],'Earth','Mars','Jupiter')
%Plot the Sun to check scale
figure(1)
r_S = 69570000;
colormap(autumn(5));
[x,y,z] = sphere(50);
r = r_S;
surf( r*x, r*y, r*z,'EdgeColor','y') % sphere with radius r_S centred at (0,0,0)
title('Orbits of Earth, Mars and Jupiter')
set(gca,'Color','k');
axis equal
hold on
function dRV = Gravity(t, RV)
global T X_J Y_J Z_J X_M Y_M Z_M X_E Y_E Z_E mu_E mu_M mu_J mu_S
X = R(1);
Y = R(2);
Z = R(3);
VX = R(4);
VY = R(5);
VZ = R(6);
% Interpolate Astronomical Data
XE = spline(T, X_E, t);
YE = spline(T, Y_E, t);
ZE = spline(T, Z_E, t);
R_E = [XE YE ZE];
XM = spline(T, X_M, t);
YM = spline(T, Y_M, t);
ZM = spline(T, Z_M, t);
R_M = [XM YM ZM];
XJ = spline(T, X_J, t);
YJ = spline(T, Y_J, t);
ZJ = spline(T, Z_J, t);
R_J = [XJ YJ ZJ];
R_S = [0 0 0];
F_E = mu_E*(R_E - R_A)/norm(R_E - R_A)^3;
F_M = mu_M*(R_M - R_A)/norm(R_M - R_A)^3;
F_J = mu_J*(R_J - R_A)/norm(R_J - R_A)^3;
F_S = mu_S*(R_S - R_A)/norm(R_S - R_A)^3;
dV = F_E + F_M + F_J + F_S;
% dR = [VX VY VZ];
% dRV = [dV dR].';
dR(1) = VX
dR(2) = VY
dR(3) = VZ
dR(4) = dV * X
dR(5) = dV * Y
dR(6) = dV * Z
Last edited by a moderator: