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

Matlab Planetary motion

  1. May 24, 2009 #1
    Hello,

    Ive got a problem in matlab and have searched the internet alot for help with this but nothing really addresses my issue specifically.\. I wish to use the ode45 solver to calculate and plot the orbits of all the planets in the solar system. I can manage just fine to get them to go around the sun but am struggling with even understanding how to get them to interact with eachother in terms of their gravity. Ive put what ive done for earth and venus below, please forgive the sloppy work I am not very good at it yet.



    function Gravity1 ()

    clear;
    format long eng;

    StartTime = 0;
    TimeStep = 24 * 3600 * 5;
    EndTime = 24 * 3600 * 365;


    %Venus initial conditions

    VX0 = 1.08e11; %Position X
    VY0 = 0; %Position Y
    VVX0 = 0; %Velocity X
    VVY0 = 35000; %Velocity Y

    %Earth intial conditions

    EX0 = 1.496e11; %Position X
    EY0 = 0; %Position Y
    EVX0 = 0; %Velocity X
    EVY0 = 29500; %Velocity Y


    TVec = StartTime:TimeStep:EndTime;

    %Column Vectors for initial conditions
    EarthVec = [VX0; VY0; VVX0; VVY0];
    VenusVec = [EX0; EY0; EVX0; EVY0];

    %Built in solver
    [TimeVec, EarthMat] = ode45(@GravityDE, TVec, EarthVec);
    [TimeVec, VenusMat] = ode45(@GravityDE, TVec, VenusVec);


    EarthXVec = EarthMat (:,1);
    EarthYVec = EarthMat (:,2);
    VenusXVec = VenusMat (:,1);
    VenusYVec = VenusMat (:,2);

    hold on;
    plot (EarthXVec,EarthYVec, '.');
    plot (VenusXVec,VenusYVec,'.');



    function dYVec = GravityDE (TVec, YVec)


    MStar = 1.99e30;
    GravConst = 6.67259e-11;

    X = YVec(1);
    Y = YVec(2);
    VX = YVec(3);
    VY = YVec(4);

    RangeSq = X.^2 + Y.^2;
    Theta = atan2(Y,X);
    AccMag = (GravConst .* MStar ./ RangeSq);
    AccX = -AccMag .* cos (Theta);
    AccY = -AccMag .* sin (Theta);

    dYVec = [VX; VY; AccX; AccY];



    Sorry about the minimal comments, Im stressed out and have butchered things here and there, I hope you can follow this. Any help on improvements would also be appreciated. I really want to understand this and figure out how to deal with more than 2 bodies.

    Thanks
     
  2. jcsd
  3. May 24, 2009 #2
    That sounds great, "certainly I'm going to learn MATLAB someday"! But I don't think its possible to make planets revolve while considering gravitational interaction b/w them. This is so because Ive learnt that nobody has ever been able to solve the 3-BODY problem.
     
  4. May 24, 2009 #3

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    You need to
    • Make a *big* vector that incorporates all of the planetary states in one structure.
    • Consider the Sun to be a planet. It orbits, too.
    • Use a better form for the gravitational acceleration (i.e., use vectors only. That atan2 stuff forces an oversimplistic form.)
    • Compute the gravitational interactions among all bodies. (i.e., you need a loop).

    If you want to be at all realistic you need to go to 3D. The planets do not orbit in the same plane.
     
  5. May 24, 2009 #4

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    That is just wrong. Of course the N-body problem is solvable. It just can't be expressed in terms of elementary functions. One has to resort to numerical integration techniques to find a solution.
     
  6. May 24, 2009 #5
    Thanks for the reply. I think I can manage most of that however what I don't know how to do is to calculate the gravitational attraction between a planet at a position and a planet at another position. Would I use a similar function to the one I had before? the mind boggles.

    My vector knowledge isn't hugely advanced unfortunately.
     
  7. May 24, 2009 #6
    In general, the gravitational force on a body of mass m1 at position r1 due to a second body, of mass m2 at position r2 is given by

    F = Gm1m2 (r2 - r1)/|r2 - r1|3

    The vector part is just the inverse square dependence on the separation distance, in the direction from body 1 to body 2.

    If this level of vector math is unfamiliar to you, then I think you might find a detailed simulation of this problem to be kind of difficult. Vectors are used precisely because they make problems like this easier.
     
  8. May 24, 2009 #7
    That does look familiar, thanks. I could do all this stuff on paper just not terribly well in matlab. I'll get there eventually.
     
  9. May 25, 2009 #8
    Ok now the next problem is that when i put the state vectors into a column as below:


    %Venus initial conditions

    VX0 = 1.08e11; %Position X
    VY0 = 0; %Position Y
    VVX0 = 0; %Velocity X
    VVY0 = 35000; %Velocity Y

    %Earth intial conditions

    EX0 = 1.496e11; %Position X
    EY0 = 0; %Position Y
    EVX0 = 0; %Velocity X
    EVY0 = 29500; %Velocity Y


    InitialStateVec = [VX0; VY0; VVX0; VVY0; EX0; EY0; EVX0; EVY0]


    I get an error in matlab when it tries to do the calculations with the GravityDE function:


    ??? Error using ==> odearguments at 116
    Solving GRAVITYDE requires an initial condition vector of length 8.


    It works fine when I do just the 4 states for the 1 planet, but it doesn't like it for 2.


    Once again, thanks for taking the time to read this.
     
  10. May 25, 2009 #9

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Did you rewrite your GRAVITYDE function to work with vectors of length 8?
     
  11. May 25, 2009 #10
    Ok that was a problem, I thought I did rewrite it and saved it but matlab kept using an older version for some reason. Figured that bit out now.

    I've tried to incorporate loops to calculate how each body affects the acceleration of the planet, in the case of 3 bodies there's a for loop with 2 nested for loops each of which calculated the acceleration due to the star and the other for the other planet. I am not sure if I did this correctly nor do I know how to set the for loop condition, which I assume is the time period over which this is all taking place...not sure!

    I'll copy what I've changed so far below, this assuming the input for the function is the column containing the 4 states of the 2 planets. Any suggestions are welcome.



    function dYVec = GravityDE (TVec, YVec)




    %Calculate the individual accelerations due to bodies then add them
    %together


    %Calculations for Venus
    for

    X = YVec(1);
    Y = YVec(2);
    VX = YVec(3);
    VY = YVec(4);
    X2 = YVec(5);
    Y2 = YVec(6);
    VX2 = YVec(7);
    VY2 = YVec(8);

    %Acceleration due to star

    for

    MBody = 2.00e30;
    [AccX,AccY] = Acceleration(X,X2,Y,Y2);
    AccX1 = AccX;
    AccY1=AccY;

    end

    %Acceleration due to Earth

    for

    Mbody = 5.98e24
    [AccX,AccY] = Acceleration (X,X2,Y,Y2);
    AccX2=AccX;
    AccY2=AccY;

    end

    %Calculate sum of Acceleration Vectors
    VenusAccX = AccX1 + AccX2;
    VenusAccY = AccY1 + AccY2;
    %Venus Velocity Vectors
    VenusVX = VX;
    VenusVY = VY;


    end


    %Calculations for Earth
    for

    X = YVec(5);
    Y = YVec(6);
    VX = YVec(7);
    VY = YVec(8);
    X2 = YVec(1);
    Y2 = YVec(2);
    VX2 = YVec(3);
    VY2 = YVec(4);

    %Acceleration due to star

    for

    MBody = 2.00e30;
    [AccX,AccY] = Acceleration(X,X2,Y,Y2);
    AccX1 = AccX;
    AccY1=AccY;

    end

    %Acceleration due to Venus

    for

    Mbody = 5.98e24
    [AccX,AccY] = Acceleration (X,X2,Y,Y2);
    AccX2=AccX;
    AccY2=AccY;

    end

    %Calculate sum of Acceleration Vectors

    EarthAccX = AccX1 + AccX2;
    EarthAccY = AccY1 + AccY2;
    %Earth Velocity vectors

    EarthVX = VX;
    EarthVY = VY;


    end



    dYVec = [VenusVX; VenusVY; VenusAccX; VenusAccY; EarthVX; EarthVY; EarthAccX; EarthAccY];




    function [AccX,AccY] = Acceleration(X,X2,Y,Y2)


    GravConst = 6.67259e-11;
    RangeSq = (X2-X).^2 + (Y2 -Y).^2;
    Theta = atan2 ((X2-X),(Y2-Y));
    AccMag = GravConst .* MBody ./ RangeSq;
    PlntAccX = -AccMag .* cos (Theta);
    PlntAccY = -AccMag .* sin (Theta);


    I also thought in future I'd include the state of the star in the input matrix for this function so that it is easier to calculate. Though, I think I'm going to just make it so that the star is stationary in space as my head is already hurting from this alone. Thanks very much so far.
    I'm still using the atan function. Once I figure this whole thing out I'll work on getting that sorted. Thanks very much so far.
     
  12. May 25, 2009 #11

    Thank You D H.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Matlab Planetary motion
  1. Matlab and tables (Replies: 0)

  2. Simple Matlab (Replies: 2)

  3. Arduino and Matlab (Replies: 1)

Loading...