reset
#set parametric
set style function lines
set size 1.0, 1.0
set origin 0.0, 0.0
set multiplot
set grid
unset key
set angles radians
set samples 250
# constants
k = 1.4
M = 0.02896
R = 8.314
L = 0.0065
g = 9.81
# standard
T0 = 288.15
rho0 = 1.225
p0 = 101325.0
#
Tr(h) = (T0-L*h)/T0
Ec = g*M/L/R

# functions
MachToMS(Ma,h) = 20.046 * sqrt(T0 - L*h) * Ma
FeetToMeter(x) = x * 0.3048
KnotToMS(x) = x * 0.51444 
MSToKnot(x) = x / 0.51444

rho(h) = rho0 * (Tr(h) ** (Ec-1.0))
p(h) = p0 * (Tr(h) ** Ec)

qi(vTAS,h) = 0.5 * rho(h) * vTAS*vTAS
qc(vTAS,h) = p(h) * (1.0 + ((k-1.0)*rho(h)*vTAS*vTAS)/(k*2.0*p(h)))**(k/(k-1.0)) - p(h)

TAS(vMAS,h) = sqrt((((0.5*rho0*vMAS*vMAS+p(h))/p(h))**((k-1.0)/k) - 1.0) / (((k-1.0)*rho(h))/(k*2.0*p(h))))

EAS(vTAS,h) = sqrt(2.0*qi(vTAS,h)/rho0)
MAS(vTAS,h) = sqrt(2.0*qc(vTAS,h)/rho0) 
CAS(vMAS,h) = TAS(MAS(vMAS,h),0)

vc(vTAS,h) = EAS(vTAS,h)-CAS(vTAS,h)

MASkt(vTASkt,hft) = MSToKnot(MAS(KnotToMS(vTASkt),FeetToMeter(hft)))
CASkt(vTASkt,hft) = MSToKnot(CAS(KnotToMS(vTASkt),FeetToMeter(hft)))
EASkt(vTASkt,hft) = MSToKnot(EAS(KnotToMS(vTASkt),FeetToMeter(hft)))

vckt(vMa,hft) = MSToKnot(vc(MachToMS(vMa,FeetToMeter(hft)),FeetToMeter(hft)))

#  Plot Magnitude Response
set title "MAS, CAS, EAS as functions of TAS"
set xlabel "TAS [kt]"
set ylabel "MAS, CAS, EAS [kt]" 1,0

set xrange [0.0:650.0]
set yrange [0.0:650.0]

plot \
  MASkt(x,00000.0), \
  MASkt(x,20000.0), \
  MASkt(x,40000.0)

plot \
  CASkt(x,00000.0), \
  CASkt(x,20000.0), \
  CASkt(x,40000.0)

plot \
  EASkt(x,00000.0), \
  EASkt(x,20000.0), \
  EASkt(x,40000.0)