%--------------------------------------------------------------------------
%
%    VMCAv1.m  Program to Calculate the Minimum Control Speed in Air
%
%               Mike Cavanaugh  3/08/04
%
%--------------------------------------------------------------------------
clear

%------Initial Conditions for a Right Engine Failure
phi = -5;                   % bank angle in degrees, FARs allow 5 deg bank toward good engine
wmax = 640000;              % maximum aircraft weight in pounds for VMCA calculation
wmin = 440000;              % minimum aircraft weight in pounds for VMCA calculation
dw = 2000;                  % weight increment in pounds for VMCA calculation
drmax =  15;                % maximum rudder deflection in degrees, must be a positive number   
damax =  -25;               % maximum aileron deflection in degrees, must be a negative number
Tmaxsl = 50000;             % maximum asymmetric thrust in pounds

%------Aircraft Geometry-------------------------
lt = 68.5 ;                 % asymmetric thrust moment arm in feet
s = 5500;                   % reference wing area in square feet
b = 195.7;                  % wing span in feet
di = 8.4;                   % engine inlet diameter in feet
dsl = .002377;              % air density at sea level in slugs per cubic foot

%------Stability Derivatives---------------------
cnbeta =  0.002618;         % yawing moment coefficient per degree of sideslip
clbeta = -0.003857;         % rolling moment coefficient per degree of sideslip
cybeta = -0.016756;         % side force coefficient per degree of sideslip
cnda =    0.000112;         % yawing moment coefficient per degree of aileron deflection
clda =    0.000805;         % rolling moment coefficient per degree of aileron deflection
cyda =    0.0;              % side force coefficient per degree of aileron deflection
cndr =   -0.001902;         % yawing moment coefficient per degree of rudder deflection
cldr =    0.000122;         % rolling moment coefficient per degree of rudder deflection
cydr =    0.003054;         % side force coefficient per degree of rudder deflection
clmax = 1.6;                % aircraft maximum lift coefficient for this test condition

%-----Calculate drag of windmilling engine (Ref. Torenbeek G-8)
cdeng = (.1934*di^2)/s;

%-----Loop in Weight for Rudder Limited VMCA-----
i = 1;
for W = wmin:dw:wmax
  
 %-----Find Rudder Limited Solution--------------
  dr = drmax;
  A11 = cybeta;   
  A12 = cyda;   
  A13 = 0;   
  A14 = ((2*W)/(dsl*s))*sin(phi/57.3);
  A21 = cnbeta;
  A22 = cnda;
  A23 = cdeng*lt/b;
  A24 = (2*Tmaxsl*lt)/(s*b*dsl);
  A31 = clbeta;
  A32 = clda;
  A33 = 0;
  A34 = 0;
  A41 = 0;
  A42 = 0;
  A43 = 1;
  A44 = 0;

  B11 = -cydr*dr;
  B21 = -cndr*dr;
  B31 = -cldr*dr;
  B41 = 1;

  A = [A11 A12 A13 A14; A21 A22 A23 A24; A31 A32 A33 A34; A41 A42 A43 A44];
  B = [B11; B21; B31; B41];

  c = inv(A)*B;

  beta = c(1,1);
  da = c(2,1);
  vfts = c(4,1)^-.5;
  q = .5 * dsl * vfts^2;
  vkts = vfts/1.15*(3600/5280);
  vskts = 0.5928854*((2*W)/(dsl*s*clmax))^.5;
  
  rudlim(i,1) = W;
  rudlim(i,2) = beta;
  rudlim(i,3) = dr;
  rudlim(i,4) = da;
  rudlim(i,5) = vkts;
  rudlim(i,6) = vskts;
 
  if abs(da) >= abs(damax) 
      break
  end   
  
  %-----plot results-----------------------------
  if i>1
   hold on
   plot(W,vkts,'-ob',W,vskts,'or')
   grid on
   box on
   xlabel('WEIGHT')
   ylabel('VMCA (ktas)')
   legend ('vmca','vstall', 4)
  end
  i=i+1;   
end

%----Loop in Weight for Aileron Limited Solution-
i = 1;
for W = wmax:-dw:wmin
    
  %-----Find Aileron Limited Solution------------
  da = damax;
  AA11 = cybeta;   
  AA12 = cydr;   
  AA13 = 0;   
  AA14 = ((2*W)/(dsl*s))*sin(phi/57.3);
  AA21 = cnbeta;
  AA22 = cndr;
  AA23 = cdeng*lt/b;
  AA24 = (2*Tmaxsl*lt)/(s*b*dsl);
  AA31 = clbeta;
  AA32 = cldr;
  AA33 = 0;
  AA34 = 0;
  AA41 = 0;
  AA42 = 0;
  AA43 = 1;
  AA44 = 0;

  BB11 = -cyda*da;
  BB21 = -cnda*da;
  BB31 = -clda*da;
  BB41 = 1;

  AA = [AA11 AA12 AA13 AA14; AA21 AA22 AA23 AA24; AA31 AA32 AA33 AA34; AA41 AA42 AA43 AA44];
  BB = [BB11; BB21; BB31; BB41];

  cc = inv(AA)*BB;

  beta = cc(1,1);
  dr = cc(2,1);
  vfts = cc(4,1)^-.5;
  q = .5*dsl*vfts^2;
  vkts = vfts/1.15*(3600/5280);
  vskts = 0.5928854*((2*W)/(dsl*s*clmax))^.5;
  
  aillim(i,1) = W;
  aillim(i,2) = beta;
  aillim(i,3) = dr;
  aillim(i,4) = da;
  aillim(i,5) = vkts;
  aillim(i,6) = vskts;
  
  if abs(dr) >= abs(drmax)
     break
  end  
  
  %-----plot results-----------------------------
  if i>1
   hold on
   plot(W,vkts,'-ob',W,vskts,'-or')    
  end
  i=i+1;
end