Code:

function calcCdi (Cl, Ar, k)  {
        return k * (Cl^2) / (Pi * Ar)
}

# ------------------------------------------------
function getAirspeed (cl, area, weight)  {
    return sqrt(weight / (0.5 * Rho * area * cl))
}

function getCl    (airspeed, area, weight) {
    return weight / (0.5 * Rho * area * airspeed^2)
}

# ------------------------------------------------
function getLift  (cl, airspeed, area)  {
    return 0.5 * Rho * area * cl * airspeed^2
}

# ------------------------------------------------
function gen (Weight, Area, Ar, Cd, col)  {
        v0     = int(getAirspeed(CLmax, Area, Weight))
        v1     = int(getAirspeed(CLmin, Area, Weight))
        N      = v1 - v0

        printf "## airspeed  %8.4f - %8.4f\n", v0, v1

    for (n = 0; n <= N; n++)  {
            airspeed  [n]  = v0 + n
            cl        [n]  = getCl(airspeed [n], Area, Weight)
            cdi       [n]  = calcCdi(cl [n], Ar, 1)

            induced   [n]  = getLift(cdi [n], airspeed [n], Area)
            parasitic [n]  = getLift(Cd,  airspeed [n], Area)

            lift      [n]  = getLift(cl [n], airspeed [n], Area)
            drag      [n]  = induced [n] + parasitic [n]

            Ld        [n]  = lift [n] / drag [n]
            sink      [n]  = 4 * airspeed [n] / Ld [n]
             
            airspeed  [n] /= 5280/3600
        }
}

# ----------------------------------------------