Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab Plotting the tragectory of an asteroid in MATLAB

  1. Mar 2, 2017 #1
    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;

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

    DrClaude

    User Avatar

    Staff: Mentor

    You'll have to describe the problem.
     
  4. Mar 4, 2017 #3

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    And we can't help debug without the data sets.....
     
  5. Mar 6, 2017 #4
    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.
     

    Attached Files:

  6. Mar 6, 2017 #5

    DrClaude

    User Avatar

    Staff: Mentor

    Details, please! What is the error message?
     
  7. Mar 6, 2017 #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. Mar 6, 2017 #7

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    Where is R defined?
     
  9. Mar 6, 2017 #8

    DrClaude

    User Avatar

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Plotting the tragectory of an asteroid in MATLAB
  1. MATLAB plots (Replies: 1)

  2. Matlab Plot (Replies: 2)

  3. Plotting in MATLAB (Replies: 8)

  4. Plotting with MATLAB (Replies: 1)

  5. MatLab Plot (Replies: 6)

Loading...