%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Patrick LeGresley
% Placard Diagram & V-n Diagram
% Last modified: April 26, 1999
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Known problems with this code:
%
% none :)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs:
%  b_wing = wing span (ft)
%  S_ref = wing reference area (ft^2)
%  taper_wing = wing taper ratio
%  sweep_wing_25c = wing quarter chord sweep angle (deg)
%  h_cruise = cruise altitude (ft)
%  V_cruise = cruise speed (KTAS)
%  W_to = takeoff weight (lbs)
%  CL_max = clean wing maximum lift coefficient at zero Mach number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs:
%  n_ult = ultimate load factor for weight calculations (g's)
%  V_C = structural design cruise speed in knots equivalent airspeed (keas)
%  V_D = design dive speed in knots equivalent airspeed (keas)
%  V_H = maximum speed at sealevel (keas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [n_ult, V_C, V_D, V_H] = vn(b_wing,S_ref,taper_wing,sweep_wing_25c,h_cruise,...
   V_cruise,W_to,CL_max)



% Misc. stuff
h_Vc = 26000;								% Vc altitude (ft)
h_max = h_cruise + 5000;			   % Maximum altitude for calculations (ft)
h_vn = 20000;								% Altitude for computing V-n diagram (ft)

% Atmospheric constants
T_o = 518.7;
rho_o = 0.002377;
R = 1716;
gamma = 1.4;
g=32.2;

% Curve fit to data from course notes to correct CLmax for Mach number effects
Mach_corr = polyfit([0 .1 .2 .3 .4 .5 .6 .7 .8 .9],...
   [1.3 1.3 1.3 1.25 1.15 1 0.85 0.8 0.75 0.85],3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computations of needed variables from the inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find speed of sound (in ft/sec) at cruise altitude
if h_cruise < 36089
      theta = 1- 6.875E-6*h_cruise;
      T = T_o*theta;
else
      T = -69.7 + 459.67;
end
a = sqrt(gamma*R*T);

% compute cruise Mach number
M_cruise = 1.6878*V_cruise/a;

% compute wing aspect ratio
AR_wing = b_wing^2/S_ref;

% find wing mean geometric chord
c_r_wing = 2*S_ref/(b_wing*(1+taper_wing));
c_t_wing = taper_wing*c_r_wing;
mgc_wing = 0.5*(c_r_wing + c_t_wing);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Placard Diagram Calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculate the structural design cruise Mach number
M_c = 1.06*M_cruise;

% Calculate the design dive Mach number
M_D = 1.07*M_c;

i=0;
for h=h_max:-100:0
   i=i+1;
   hplot(i) = h;
   
   % Calculate temperature and density at given altitude
   if h < 36089
      theta = 1- 6.875E-6*h;
      T = T_o*theta;
      rho = rho_o*theta^4.2561;
   else
      T = -69.7 + 459.67;
      rho = rho_o*0.2971*exp(-((h-36089)/20806.7));
   end
   
   % Calculate speed of sound
   a = sqrt(gamma*R*T);					
   
   if h >= h_Vc
      
      % Calculations for altitudes at or above Vc altitude
      
      % Calculate design cruise speed and dive speed in KTAS
      V_c(i) = 0.592*M_c*a;
      V_D(i) = 0.592*M_D*a;
      
      % Calculate corresponding dynamic pressures
      q_c = 0.5*rho*V_c(i)^2;
      q_D = 0.5*rho*V_D(i)^2;
      
   else
      
      % Calculations for altitudes below Vc altitude
      
      % Calculate design cruise speed and dive speed in KTAS
      V_c(i) = sqrt(2*q_c/rho);
      V_D(i) = 1.15*V_c(i);
      
      % Extend the constant Mach number line to meet the constant EAS line
      if 0.592*M_D*a < V_D(i)
         V_D(i) = 0.592*M_D*a;
      end
      
   end
   
end



% Plot the placard diagram
%figure(1)
%plot(V_c,hplot,'b-',V_D,hplot,'r--')
%title('Placard Diagram')
%xlabel('True Airspeed (knots)')
%ylabel('Altitude (feet)')
%legend('Cruise','Dive',2)
%grid on

V_H = V_c(length(hplot));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% V-n Diagram Calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculate atmospheric properties
if h_vn < 36089
   theta = 1- 6.875E-6*h_vn;
   T = T_o*theta;
   rho = rho_o*theta^4.2561;
   sigma=rho/rho_o;
else
   T = -69.7 + 459.67;
   rho = rho_o*0.2971*exp(-((h_vn-36089)/20806.7));
end
sigma=rho/rho_o;
a=sqrt(gamma*R*T);

% Find Vc and Vd at the current altitude
V_c_keas = V_c(h_vn/100)*sqrt(sigma);
V_D_keas = V_D(h_vn/100)*sqrt(sigma);

% Compute the semi-chord sweep angle from the quarter-chord sweep angle
sweep_wing_0c = tan(sweep_wing_25c/57.3) + 4*(0.25)*...
   (1-taper_wing)/AR_wing/(1+taper_wing);
sweep_wing_50c = tan(sweep_wing_0c) - 4*(0.5)*...
   (1-taper_wing)/AR_wing/(1+taper_wing);

i=0;
for V=1:V_D_keas;
   i=i+1;
   
   V_keas(i) = V;
   
   % Calculate true airspeed in ft/sec corresponding to the current EAS
   v=1.6878*V/sqrt(sigma);
   
   % Calculate Mach number
   M=v/a;
   
   % Correct the low speed CLmax for Mach number effects
   Mprime = M*cos(sweep_wing_25c/57.3)/cos(24.5/57.3);
   CL_max_Mach = CL_max*polyval(Mach_corr,Mprime)/polyval(Mach_corr,0);
   
   % Calculate the stall line
   n_max(i) = v^2*rho*S_ref*CL_max_Mach/(2*W_to);
   
   % Calculate CLalpha, including Mach number effects
   kappa = 0.97;
   beta = sqrt(1-M^2);
   C_L_alpha = 2*pi*AR_wing/(2+sqrt(AR_wing^2*beta^2/kappa^2*...
      (1+tan(sweep_wing_50c/57.3)^2/beta^2)+4));
   
   % Find values needed for gust load factors
   mu_g_gust = 2*(W_to/S_ref)/(rho*mgc_wing*g*C_L_alpha);
   K_g = 0.88*mu_g_gust/(5.3+mu_g_gust);
   
   % Get the gust velocities corresponding to the current altitude
   if h < 20000
      U_de_V_B = 66;
      U_de_V_C = 50;
      U_de_V_D = 25;
   elseif h> 50000
      U_de_V_B = 38;
      U_de_V_C = 25;
      U_de_V_D = 12.5;
   else
      U_de_V_B = 84.67 - 0.000933*h;
      U_de_V_C = 66.67 - 0.000833*h;
      U_de_V_D = 33.34 - 0.000417*h;
   end
   
   % Compute the gust lines
   n_lim_gust_V_B(V) = 1 + V*(K_g*U_de_V_B*C_L_alpha)/(498*W_to/S_ref);
   n_lim_gust_V_C(V) = 1 + V*(K_g*U_de_V_C*C_L_alpha)/(498*W_to/S_ref);
   n_lim_gust_V_D(V) = 1 + V*(K_g*U_de_V_D*C_L_alpha)/(498*W_to/S_ref);
   
end


% Calculate positive limit maneuvering load factor, n_lim_pos
n_lim_pos = 2.1 + (24000/(W_to + 10000));

% n_lim_pos must be greater than 2.5 at all times
if n_lim_pos < 2.5
   n_lim_pos = 2.5;
end

% n_lim_pos need not be greater than 3.8 at W_to
if n_lim_pos > 3.8
   n_lim_pos = 3.8;
end

% create vector of n_lim_pos for plot
n_lim_pos = n_lim_pos*ones(1,(length(V_keas)));

% find the minimum load factor at each velocity for plotting
for j=1:length(V_keas)
   n_pos(j) = min([n_max(j) n_lim_pos(j)]);
end

% find V_B
k_VB = 10^9;
for k = 2:length(n_max);
   if n_lim_gust_V_B(k) < n_max(k)
      if k < k_VB
         k_VB = k;
      end
   end   
end

if k_VB > length(n_max)
   k_VB = length(k);
end

V_B_keas = V_keas(k);

if V_B_keas > sqrt(sigma)*sqrt(n_lim_gust_V_B(round(V_c_keas)))*sqrt(2*W_to/(rho*S_ref*CL_max))/1.6878
   V_B_keas = sqrt(sigma)*sqrt(n_lim_gust_V_B(round(V_c_keas)))*sqrt(2*W_to/(rho*S_ref*CL_max))/1.6878;
end

% find the stall line for the gust diagram
for k=1:V_B_keas
   v_gust(k) = k;
   n_lim_gust(k) = n_max(k);
end

% Add the data points at V_C and V_D for plotting
v_gust(k+1) = V_c_keas;
n_lim_gust(k+1) = n_lim_gust_V_C(floor(V_c_keas));
v_gust(k+2) = floor(V_D_keas);
n_lim_gust(k+2) = n_lim_gust_V_D(floor(V_D_keas));
v_gust(k+3) = floor(V_D_keas);
n_lim_gust(k+3) = 0;
V_keas(length(V_keas)+1) = floor(V_D_keas);
n_pos(length(V_keas))= 0;


% Plot the V-n diagram
%figure(2)
%plot(V_keas,n_pos,'b-',v_gust,n_lim_gust,'r--')
%grid on
%title('V-n Diagram')
%xlabel('Speed (keas)')
%ylabel('Load Factor')
%legend('Maneuver','Gust',2)


% Define variables to be passed out of the function
n_ult = 1.5*max([max(n_pos) max(n_lim_gust)]);
V_C = V_c_keas;
V_D = V_D_keas;