MATLAB Plotting the tragectory of an asteroid in MATLAB

AI Thread Summary
The discussion focuses on troubleshooting an issue with plotting an asteroid's trajectory in MATLAB using the ode23 solver. The user has correctly set up the gravitational forces from the Sun, Earth, Mars, and Jupiter, but encounters an error indicating that the variable 'R' is undefined in the Gravity function. The problem arises because the function is using 'R' instead of the provided input 'RV', leading to confusion in the code. To resolve the issue, the variable 'R' should be replaced with 'RV' to ensure the correct input is utilized. This adjustment should allow the ode23 solver to function properly and plot the asteroid's trajectory.
JayFlynn
Messages
7
Reaction score
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?

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
<Moderator's note: please use code tags>
 
Last edited by a moderator:
Physics news on Phys.org
JayFlynn said:
Any idea why the ode solver doesn't work?
You'll have to describe the problem.
 
And we can't help debug without the data sets...
 
Okay thanks guys, I am new to all this. Heres the data files. When i run the code it says there is an error with my input arguments in the line where i call the ode solver.
 

Attachments

JayFlynn said:
When i run the code it says there is an error with my input arguments in the line where i call the ode solver.
Details, please! What is the error message?
 
Undefined function 'R' for input arguments
of type 'double'.

Error in SEMJ_5>Gravity (line 99)
X = R(1);

Error in SEMJ_5>@(t,RV)Gravity(t,RV) (line
60)
[T,RV] = ode23(@(t,RV) Gravity(t, RV),T',
init_val);

Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:}); % ODE15I
sets args{1} to yp0.

Error in ode23 (line 114)
odearguments(FcnHandlesUsed,
solver_name, ode, tspan, y0, options,
varargin);

Error in SEMJ_5 (line 60)
[T,RV] = ode23(@(t,RV) Gravity(t, RV),T',
init_val);
 
JayFlynn said:
Undefined function 'R' for input arguments
of type 'double'.

Error in SEMJ_5>Gravity (line 99)
X = R(1);
Wasn't that error message clear enough?!?

You are using R and dR in your function Gravity, while the input is RV and the output dRV.
 

Hot Threads

Recent Insights

Back
Top