clear
clf
global N1 N2 N3

f = 0;
z0 = 500; % initial position (m)
w0 = 5; % initial wind velocity (m/s)
TT  = 3600; % run time (s)
dt = 1; % step size (s)
%adiobatic lapse rates
G_A = 9.76e-3; % adiabatic lapse rate (K/m)
G_AJ = 1.963e-3; %Adiobatic lapse rate of Jupiter
G_AM = 4.5e-3; %adiobatic Lapse rate of Mars
%measured lapse rate of tropospheres
G = 6.4e-3; % measured lapse rate (K/m) Earth
G_J = 1.2e-3; %Measured lapse rate of Jupiter
G_M = 0.000998; %measured lapse rate of Mars
%gravities of planets
g = 9.806; % gravitational acceleration
g_J = 25.95; %Jupiter Gravity
g_M = 3.77; %marshan gravity
T0 = 287; % surface temperature
T0M = 242; %surface temp of Mars
T0J = 170; %surface temprature of Jupiter

n=floor(TT/dt); % number of steps
T = T0 - G*z0; % temperature at initial position (K)
TM = T0M - G_M*z0;
TJ = T0J - G_J*z0;
N1 = sqrt((g*(G_A - G))/T); % Brunt-Vaisala frequency if Earth
N2 = sqrt((g_M*(G_AM - G_M))/TM); % Brunt-Vaisala frequency of Mars
N3 = sqrt((g_J*(G_AJ - G_J))/TJ); % Brunt-Vaisala frequency of Jupiter
Y = [z0 ; w0; z0 ; w0 ; z0 ; w0];
t = 0;
for i = 1:n
    f1 = spring_eq(t     , Y                );
    f2 = spring_eq(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = spring_eq(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = spring_eq(t+dt  , Y+f3 .*   dt     );
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    plot( t , Y(1), 'b', t , Y(3) , 'r', t , Y(5) , 'g');
    axis([0 TT -1e3 1e3]);
    hold on
    pause(0.01);
    if mod(t + 35 , 36) == 0;
        f = f + 1;
        saveas(gcf, ['planet_',sprintf('%03d',f)],'png');
    end
    t = t + dt;
end