Classical Mechanics MatLab problem (oribtal motion)

[Rhabdovirus]

Hey guys, I've been fighting this problem all weekend with little avail. I only know a little Java, but we have this final project in MatLab, a program I don't really know.

Homework Statement

I need to compute the orbit of a comet using the Runge-Kutte, Euler and midpoint methods. The program should terminate after 1 orbit. I then need to compute eccentricity, a, period and r_{min}.

Homework Equations

I'm using the differential equation \mu \ddot{r} = F(r) + l²/\mu²*r²
F(r) is the central force, -GMm/r² in this case

The Attempt at a Solution

Here's my code so far, I'm having trouble figuring out how to plot it

%Differential equation:  r'' = F(r)/mu + ((l^2)/((mu^2)*(r^2)))
clear all; help orbit;
%Variables
%Fr = G*Msun*Mcomet/r^2;
%l = Mcomet*v*r;
%Constants
Msun = 2e30;
Mcomet = 2e6;
G = 6.67e-11;
mu = (Msun*Mcomet)/(Msun+Mcomet);
time = 0;

%initial conditions
r0 = 250000; %initial radius from Sun to comet in m
v0 = 8000; %initial tangential velocity of comet in m/s
r = r0; v = v0;
tau = 1000; %timestep in seconds
maxstep = 1000;
for iStep = 1:maxstep
   

  
%Euler Method v_n+1 = (v_n + tau*acc_n)
%				  r_n+1 = (r_n + tau*v_n)

acc = -G*Msun*r./norm(r)^3;
r = r + tau*v;
v = v + tau*acc;
time = time + tau;%Midpoint v_n+1 = v_n + tau*a_n
%			 r_n+1 = r_n + tau*v_n + .5*a_n*tau^2

acc = -G*Msun*r./norm(r)^3;
v = v + tau*acc;
r = r + tau*v + .5*acc*(tau^2);end

Homework Statement

Homework Equations

The Attempt at a Solution

 

[Rhabdovirus]

ugh, 75 views and no one knows!

here's an update: I figured out the algorithms, but I'm getting a straight line not the oscillatory one I should be getting.

%Complete the orbit of a comet around the Sun.  (1) Have the program stop when it has completed one orbit.
%(2) Have the program compute the period (period), eccentricity (sigma), semimajor axis (a) and perihelion distance (rmin). 


%Differential equation:  r'' = F(r)/mu + ((l^2)/((mu^2)*(r^2)))
clear all; help orbit;
%Variables
%Fr = G*Msun*Mcomet/r^2;
%l = Mcomet*v*r;
%Constants
Msun = 2e30;
Mcomet = 2e6;
G = 6.67e-11;
mu = (Msun*Mcomet)/(Msun+Mcomet);

%initial conditions

r(1) = 2500000; %initial radius from Sun to comet in m
v(1) = 8000; %initial tangential velocity of comet in m/s

l = Mcomet*r(1)*v(1);

tau = 1; %timestep in seconds
maxstep = 1000;
acc(1) = (((-1)*G*Msun*Mcomet)/((mu)*(r(1)^2))) + ((l^2)/((mu^2)*(r(1)^3)));
t(1)= tau
for i = 2:maxstep
 
  
%Euler Method v_n+1 = (v_n + tau*acc_n)
%				  r_n+1 = (r_n + tau*v_n)

r(i) = r(i-1) + tau*v(i-1);
acc(i) = (((-1)*G*Msun*Mcomet)/((mu)*(r(i)^2))) + ((l^2)/((mu^2)*(r(i)^3)));
v(i) = v(i-1) + tau*acc(i-1);
t(i) = (i-1)*tau;

plot(t,r);
end

 

[kuruman]

To get an orbit you need to solve a different diff-eq. Starting from the one you have, you need to change variables and write a second order diff-eq. where the independent variable is the azimuthal angle φ and the dependent variable is 1/r. This is standard stuff, treated in any intermediate level Classical Mechanics textbook.

If you can't find it, I can post the equation for you.