- #1

- 8

- 0

Ive got a problem in MATLAB and have searched the internet a lot 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 each other in terms of their gravity. I've put what I've 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, I am 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