%Climb preformance calculation for dhc6 twin otter.

clear all;
clf;

% Aircraft data and ISA
Weight=5760.0*9.8;
Power_sl=2.0*690.0*745.;
T0=288.;
g=9.8;
R=287.;
S=39.0;
rho0=1.225;
L=0.0065;
P0=101300.;
gamma=1.4;
eff=0.8;
vals=[0.01:1.:50.];

%Altitude and Velocity range
h=[0.:500.:40000.];
V=[0.:2.5:200.];

for i=1:81,
 alt=h(i)*0.3048;
 temp=T0-alt*L;
 press=P0*(temp/T0)^(g/L/R);
 rho=press/(R*temp);
 Power=Power_sl*rho/rho0;
 a=sqrt(gamma*R*temp); 
 for j=1:81,
  Vel=V(j);
  Thrust=Power*eff/Vel;
  Cl=Weight/(0.5*rho*Vel*Vel*S);
  M=Vel/a;
  if (M<0.99),
% compressibility correction
   beta=sqrt(1.-M*M);
   Cd=(0.02+0.0425*Cl*Cl)/beta;
   D=Cd*(0.5*rho*Vel*Vel*S);
% Specific Excess Power
   Ps(i,j)=Vel/Weight*(Thrust-D);
  else
   Ps(i,j)=0;
  end;
% Specific Energy
  he(i,j)=alt+Vel*Vel/2./g;
 end;
end;

% Contour Specific Excess Power and Specific Energy
contour(V,h,Ps,vals,'-w');
hold on;
contour(V,h,he,h,'-g');
xlabel('Velocity (m/s)');
ylabel('Altitude (ft)');
title('Specific Excess Power Contours');
disp('---Paused---');
pause;

clf;
hold on;
k=1;
n=0;
hestart=250;
for i=1:81,
 clear x;
 clear y;
 ymin0=2.0;
 xmin0=0.;
 for j=1:81,
  x(j)=he(i,j);
  y(j)=Ps(i,j);
  if y(j)<0.5,
   y(j)=0.5;
  end;
  y(j)=1./y(j);
  if i==1,
   if x(j)>hestart,
    if y(j)<ymin0,
     n=n+1;
     ymin(n)=y(j)
     xmin(n)=x(j)
    end;
   end;
  end;
  if y(j)<ymin0,
   ymin0=y(j);
   xmin0=x(j);
  end;
 end;
 if i>1,
  n=n+1;
  ymin(n)=ymin0;
  xmin(n)=xmin0;
 end;
 k=k+1;
 if k==2,
  k=0;
 else
  plot(x,y,'-r');
 end;
end;
plot(xmin,ymin,'-w');
xlabel('he');
ylabel('1/Ps');
title('Path Time Prediction');

sum=0.0;
for i=1:(n-1),
 sum=sum+(xmin(i+1)-xmin(i))*(ymin(i+1)+ymin(i))*0.5
end;


disp('---Paused---');
pause;