Plotting the tragectory of an asteroid in MATLAB

Click For Summary

Discussion Overview

The discussion revolves around plotting the trajectory of an asteroid in MATLAB using the ode23 solver. Participants are addressing issues related to the implementation of the code, specifically focusing on errors encountered during the numerical integration process. The scope includes technical debugging and code correction.

Discussion Character

  • Technical explanation
  • Debate/contested

Main Points Raised

  • One participant describes their attempt to plot an asteroid's trajectory using ode23, mentioning the gravitational influences of the Sun, Earth, Mars, and Jupiter.
  • Another participant requests a description of the problem to assist with debugging.
  • Concerns are raised about the lack of data sets necessary for debugging the code effectively.
  • A participant reports an error message indicating an undefined function 'R' when calling the Gravity function, suggesting a potential issue with variable naming or scope.
  • Further clarification is sought regarding the definition of 'R', with emphasis on the mismatch between input and output variable names in the Gravity function.

Areas of Agreement / Disagreement

Participants do not reach a consensus on the solution to the problem, as the discussion highlights multiple viewpoints regarding the error and the necessary corrections. The issue remains unresolved.

Contextual Notes

The discussion reveals limitations in the code related to variable definitions and function inputs, which may affect the execution of the ode solver. Specific assumptions about variable naming conventions and data availability are also noted.

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.