# Matlab Plotting the tragectory of an asteroid in MATLAB

Tags:
1. Mar 2, 2017

### JayFlynn

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?

Code (Matlab M):

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;

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: Mar 3, 2017
2. Mar 3, 2017

### Staff: Mentor

You'll have to describe the problem.

3. Mar 4, 2017

### Dr Transport

And we can't help debug without the data sets.....

4. Mar 6, 2017

### JayFlynn

Okay thanks guys, im 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.

File size:
447.4 KB
Views:
31
File size:
447.4 KB
Views:
33
File size:
447.4 KB
Views:
26
5. Mar 6, 2017

### Staff: Mentor

Details, please! What is the error message?

6. Mar 6, 2017

### JayFlynn

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);

7. Mar 6, 2017

### Dr Transport

Where is R defined?

8. Mar 6, 2017

### Staff: Mentor

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.