#! /usr/bin/env julia

htab = [3280.84 * x for x in (0.0, 11.0, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852)]
ttab = [1.8 * x for x in (288.15, 216.65, 216.65, 228.65, 270.65, 270.65,
                          214.65, 189.946)]
ptab = [2116.21662 * x for x in  (1.0, 2.233611e-1, 5.403295e-2, 8.5666784e-3,
                        1.0945601e-3, 6.6063531e-4, 3.9046834e-5, 3.68501e-6)]

gtab = fill(0.0, NTAB) 

for i = 1:NTAB-1
    gtab[i] = ((ttab[i+1] - ttab[i]) / (htab[i+1] - htab[i]))
end

gtab[NTAB] = 0.0

function atmosphere(alt::Float64)
    h::Float64 = alt * REARTH / (alt + REARTH)      # convert geometric to
                                                    # geopotential altitude
    i = 1
    j = NTAB
    while true
        k = (i + j) / 2

        if h < htab[k]
            j = k
        else
            i = k
        end

        if j <= i+1
            break
        end
    end

    tgrad = gtab[i]                                 # i will be in 1:NTAB-1
    tbase = ttab[i]
    deltah = h - htab[i]
    tlocal = tbase + tgrad * deltah
    theta = tlocal / ttab[i]                        # temperature ratio

    if tgrad == 0.0
        delta = ptab[i] * exp(-GMR * deltah / tbase)
    else
        delta = ptab[i] * (tbase / tlocal)^^(GMR/tgrad)
    end

    sigma = delta / theta                           # density ratio

    return delta, theta, sigma
end