function [T,p,rho,cs]=atmosphere(h)
    g=9.80665;
    R=287.053;    
    if (h>=0)&(h<=11000) then
        T0=288.15;
        a1=-0.0065;
        T=T0-0.0065*h;
        p0=101325;
        p=p0*(T/T0)^(-g/(a1*R));
    elseif (h>11000)&(h<=20000)
        T=216.65;
        p1=22632.1;
        p=p1*exp(-g/(R*T)*(h-11000))
    elseif (h>20000)&(h<=32000)
        T2=216.65;
        a2=0.001;
        T=T2+0.001*(h-20000);
        p2=5474.89;
        p=p2*(T/T2)^(-g/(a2*R));
    elseif (h>32000)
        T3=228.65;
        a3=0.0028;
        T=T3+0.0028*(h-32000);
        p3=868.019;
        p=p3*(T/T3)^(-g/(a3*R));
    end
    rho=p/(R*T);
    cs=sqrt(1.4*R*T)
endfunction

function chute(dt)
g=9.6;
hp=39000;
vp=0;
i=1;
while (hp>0)
    [T,p,rho,cs]=atmosphere(hp)
    vs=vp+dt*(-g+(0.5*1.3*0.4/110)*vp^2*rho);
    hs=hp+dt*vs;
    //disp(-vs/cs)
    hp=hs;
    vp=vs;
    h(i)=hp;
    c(i)=cs;
    v(i)=vp;
    i=i+1;
end
clf
scf(0);
plot(h,-3.6*v, "r")
plot(h,3.6*c, "g")
scf(1);
plot2d(h, -v);
//plot(t, f)
endfunction

