rough estimate of 2-33 polar

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
}
}
# ---------------------------------------------------------
BEGIN {
Pi = atan2(0, -1)
Rho = 0.002378 # slugs / cu.ft.
# general
CLmax = 1.5
CLmin = 0.175
# 2-33 parameters
Empty = 600 # lbs
Cargo = 300
Weight = Empty + Cargo # lbs
Span = 51 # ft
Area = 220 # sq.ft.
Ar = 11.85
# estimates
Cd = 0.007
gen(Weight, Area, Ar, Cd)
}

a 2nd polar plot is added. It shows polars with weight increases equivalent to

lift increases necessitated by bank angle increases of 15, 30, 45 and 60.

It also shows AOA/Cl corresponding to min-sink speed in each plot.