Plotting the tragectory of an asteroid in MATLAB

In summary, the conversation discusses the process of plotting 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. The individual has chosen initial conditions for the asteroid and believes their forces are correct, but they are having issues getting the ode23 solver to work. They ask for help with debugging and provide the necessary data files.
  • #1
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?

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
  • #2
JayFlynn said:
Any idea why the ode solver doesn't work?
You'll have to describe the problem.
 
  • #3
And we can't help debug without the data sets...
 
  • #4
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

  • 1996-2017Mars.txt
    447.4 KB · Views: 682
  • 1996-2017Jupiter.txt
    447.4 KB · Views: 630
  • 1996-2017Earth.txt
    447.4 KB · Views: 571
  • #5
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?
 
  • #6
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);
 
  • #8
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.
 

1. What is the purpose of plotting the trajectory of an asteroid in MATLAB?

The purpose of plotting the trajectory of an asteroid in MATLAB is to visually represent the path of the asteroid as it moves through space. This can help scientists make predictions about its future path and potential impact on Earth.

2. How is MATLAB used to plot the trajectory of an asteroid?

MATLAB is used to plot the trajectory of an asteroid by utilizing mathematical equations and data points to calculate and graph the asteroid's position over time. This allows for a more accurate and detailed representation of the asteroid's path.

3. What data is needed to plot the trajectory of an asteroid in MATLAB?

The data needed to plot the trajectory of an asteroid in MATLAB includes the asteroid's initial position, velocity, and mass, as well as any external forces acting on the asteroid. This data can be obtained from observations, simulations, or other scientific studies.

4. Can MATLAB predict the exact path of an asteroid?

No, MATLAB cannot predict the exact path of an asteroid. While it can accurately calculate and plot the trajectory based on given data, there are many factors that can affect the asteroid's path, such as gravitational forces from other objects and collisions with debris.

5. How can plotting the trajectory of an asteroid in MATLAB contribute to asteroid impact prevention?

By plotting the trajectory of an asteroid in MATLAB, scientists can better understand the potential path and impact of an asteroid and make more accurate predictions. This information can then be used to develop strategies and technologies for potentially deflecting or destroying an asteroid before it reaches Earth.

Hot Threads

Recent Insights

Back
Top