function D=density(y)
H=y*3.28084;
N=length(H);
for i=1:N
    h=H(i);
    if h<36152
        D(i)=2116*(1.000192827-0.6864635557e-5*h)^5.256/...
            (8.911266*10^5-6.11608*h);
    elseif (h>36152 & h<82345) 
        D(i)=0.7066418961e-3*exp(1.73-0.48e-4*h);
    elseif h>82345
        D(i)=51.97/((.6529822044+0.4205343863e-5*h)^11.388...
            *(4.3748870*10^5+2.81752*h));
    else
        D(i)=0;
    end
end
D=D*515.3788;