orbital_model7

Orbital mechanics--with RK4 method.
Written by R.Sonnenfeld, New Mexico Tech Physics
 Version 1.0  11/7/2008.

USAGE orbital_model7(v2,150E9,1.99E30,1000);
The Runge-Kutta notation here is the same as used by Weisstein in
Weisstein, Eric W. "Runge-Kutta Method." 
From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Runge-KuttaMethod.html
================
SET INITIAL CONDITIONS
x(1)=z0(1);y(1)=z0(2);
vx(1)=z0(3);vy(1)=z0(4);

function [time, x,y ]=orbital_model7(v0,R0,M,npts)
%Orbital mechanics--with RK4 method.
%Written by R.Sonnenfeld, New Mexico Tech Physics
% Version 1.0  11/7/2008.
%
%USAGE orbital_model7(v2,150E9,1.99E30,1000);
%The Runge-Kutta notation here is the same as used by Weisstein in
%Weisstein, Eric W. "Runge-Kutta Method."
%From MathWorld--A Wolfram Web Resource.
%http://mathworld.wolfram.com/Runge-KuttaMethod.html
%================
%SET INITIAL CONDITIONS
%x(1)=z0(1);y(1)=z0(2);
%vx(1)=z0(3);vy(1)=z0(4);
z0=zeros(1,4);
z0(1)=R0;
z0(4)=v0;
%================
%SET CONSTANTS
global GM
G=6.67E-11;
GM=G*M;   %Calculate constants
%================
%SET PARAMETERS
%npts=10000;

%================
%SET TIME RANGE
t0=0;
tf=18*365*86400;  %18 years
%================
%ACTUALLY SOLVE ODE
t=linspace(0,tf,npts);
    [t,z]=RK4(t,z0);

% Note: The initial conditions go in z0.
% z0 is an array, not just a single #
% For this problem (a second order ODE) z0 contains initial x, y positions and vx, vy
% velocities.
%================
%UNPACK THE SOLUTIONS INTO CONVENIENTLY NAMED COLUMN VECTORS
x=z(:,1); y=z(:,2);
vx=z(:,3); vy=z(:,4);
time=t;
end

function [t,z]=RK4(t,z0)
global GM
% Note: The initial conditions go in z0.
% z0 is an array, not just a single #
% For this problem (a second order ODE)
%z0 contains initial x, y positions and vx, vy velocities.
% The array t is not assumed to be evenly spaced (but it can be);
% The first returned data point corresponds to t(1) and the last
% to t(length(t)).
len=length(t); 
z=zeros(len,4); %Set up a matrix for x, y, vx, vy
x=zeros(len,1);y=zeros(len,1);vx=zeros(len,1);vy=zeros(len,1);
%Plug initial conditions into convenient variables
x(1)=z0(1);y(1)=z0(2);
vx(1)=z0(3);vy(1)=z0(4);
for n=1:len-1
    dt=t(n+1)-t(n);
    %k1
    magr1=sqrt((x(n)^2+y(n)^2))^3;
    k1vx=dt*(-GM*x(n)/magr1);
    k1vy=dt*(-GM*y(n)/magr1);
    k1x=dt*vx(n);
    k1y=dt*vy(n);
    %k2
    magr2=sqrt(((x(n)+k1x/2)^2+(y(n)+k1y/2)^2))^3;
    k2vx=dt*(-GM*(x(n)+k1x/2)/magr2);
    k2vy=dt*(-GM*(y(n)+k1y/2)/magr2);
    k2x=dt*(vx(n)+k1vx/2);
    k2y=dt*(vy(n)+k1vy/2);
    %k3
    magr3=sqrt(((x(n)+k2x/2)^2+(y(n)+k2y/2)^2))^3;
    k3vx=dt*(-GM*(x(n)+k2x/2)/magr3);
    k3vy=dt*(-GM*(y(n)+k2y/2)/magr3);
    k3x=dt*(vx(n)+k2vx/2);
    k3y=dt*(vy(n)+k2vy/2);
    %k4
    magr4=sqrt(((x(n)+k3x)^2+(y(n)+k3y)^2))^3;
    k4vx=dt*(-GM*(x(n)+k3x)/magr4);
    k4vy=dt*(-GM*(y(n)+k3y)/magr4);
    k4x=dt*(vx(n)+k3vx);
    k4y=dt*(vy(n)+k3vy);
    %%
    vx(n+1)=vx(n)+k1vx/6+k2vx/3+k3vx/3+k4vx/6;
    x(n+1) = x(n)+k1x/6 +k2x/3 +k3x/3 +k4x/6;
    vy(n+1)=vy(n)+k1vy/6+k2vy/3+k3vy/3+k4vy/6;      
    y(n+1) = y(n)+k1y/6 +k2y/3 +k3y/3 +k4y/6;    
end
%keyboard
z(:,1)=x;  %1st column
z(:,2)=y;  %2nd column
z(:,3)=vx; 
z(:,4)=vy;
end