load engine.data
for i = 1:24,
  a(i,1) = 1;
  a(i,2) = engine(i,1);
  a(i,3) = a(i,2) * a(i,2);
  a(i,4) = engine(i,2)/10000.;
  a(i,5) = a(i,4) * a(i,4);
  a(i,6) = a(i,2) * a(i,4);
  b(i)=engine(i,3) - engine(i,4);
  c(i) = engine(i,5);
  TSFCtrue(i) = c(i)/b(i);
  mach(i) = engine(i,1);
  alt(i) = engine(i,2);
 end
end

aa= a'* a;
bb= a'* b';
cc=a'*c';
ainv = inv(aa);
coeff = ainv * bb
coeff2 = ainv*cc

Tcalc = a*coeff
TSFC = a*coeff2;

for i= 1:24
  TSFC2(i) = TSFC(i)/Tcalc(i);
end

TSFC2

for i = 1:17,
  for j = 1:17,
  machmsh(i) = 0.05*i;
  altmsh(j) = (j-1)/4.;
  dum(1) = 1;
  dum(2) = machmsh(i);
  dum(3) = machmsh(i)*machmsh(i);
  dum(4) = altmsh(j);
  dum(5) = altmsh(j) * altmsh(j);
  dum(6) = altmsh(j) * machmsh(i);
  Tcmsh(i,j) = dum*coeff;
  TSFCmsh(i,j) = dum*coeff2;
  TSFC2msh(i,j) = TSFCmsh(i,j)/Tcmsh(i,j);
 end
end

clf reset
hold on
view(0,90)
AXIS([0 4 0 1 0 0.5])
ylabel('Mach number')
xlabel('Altitude (x10^4)')
zlabel('TSFC (1/hr)')
mesh(altmsh,machmsh,TSFC2msh)