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