% *********************************************
% RunFirst - Runs each of the stability derivitives functions
% and sets up the array called constant that is used by other
% MATLAB programs to perform various dynamic analysis.
%
% *********************************************     
%
% A&AE 421 Fall 2001 - Purdue University
% 
% Note: This code is provided for a first order approximation of the dynamic 
%       stability and control derivatives of an airplane.
%
% Equations/Figures can be found in :
% 
% (Ref.1) Roskam, Jan. "Airplane Flight Dynamics and Automatic Flight
%         Controls"
%         Published by DARcorporation
%         120 E. Ninth St., Suite 2
%         Lawarence, KS 66044
%         Third Printing, 2001.
%
% (Ref.2) Roskam, Jan. "Methods for Estimating Satbility and
%         Control Derivatives of Conventional Subsonic Airplanes"
%         Published by the Author 
%         519 Boulder 
%         Laurance, Kansas 66044
%         Third Printing, 1997.
%  
% (Ref.3) Roskam, Jan. "Airplane Design: Part IV: Preliminary Calculation
%         of Aerodynamic, Thrust and Power Characteristics"
%         Published by Roskam Aviation and Engineering Corporation
%         Rt4, Box 274
%         Ottawa, Kansas 66067
%         Second Printing, 1990.
%
clear all
close all

format short g
% The next script gather data needed for the stability and 
% control derivative computations.
%BasicConstants_BoilerXP % Use this script to model the Boiler Xpress aircraft
BasicConstants_Cessna182 % Use this script to model the Cessna 182 aircraft
%BasicConstants_YourCraft % Write this script in the same format as above
%   to model your aircraft.

DerivedConstants  % This script computes some intermediate constants used below


% The 57 "constant" computed below are used by four dynamic and control
% software programs. Specifically, Simulink scripts FlatEarth.mdl and E_Earth.mdl
% use the first 547 "constant". These MATLAB programs perform 6-degree-of-freedom 
% aircraft simulation over a flat earth or an elliptical earth.
% Two other programs do simplified dynamic modeling (compute transfer functions)
% for the longitudinal or lateral-directional degrees of freedom.  They use all
% 67 of the "constant" defined below.
% These programs are called LatSC.m or LongSC.m. Examples of these programs are called 
% CessnaLongSC.m and CessnaLatSC.m. This software can be found at the 
% following web site.
%
% http://roger.ecn.purdue.edu/~andrisan/Courses/AAE421_S2001/Docs_Out/DC_Software/index.html
%
% Mass related inputs
constant(1)=W;      % W, Weight, pounds (lbf)
constant(2)=32.17405;  % g, Acceleration of gravity, ft/(sec*sec)
constant(3)=constant(1)/constant(2); % mass, slugs
constant(4)=Ixx;    % Ixx, slug*ft*ft
constant(5)=Iyy;   % Iyy, slug*ft*ft
constant(6)=Izz;   % Izz, slug*ft*ft
constant(7)=Ixz;     % Ixz, slug*ft*ft
constant(8)=eta_p;    % propeller efficiency, eta, nondimensional
constant(9)=0; % unassigned

% Derived constants from the inertia data
constant(10)=constant(4)*constant(6)-constant(7)*constant(7);  %gamma
constant(11)=((constant(5)-constant(6))*constant(6)-constant(7)*constant(7))/constant(10);% c1
constant(12)=(constant(4)-constant(5)+constant(6))*constant(7)/constant(10);% c2
constant(13)= constant(6)/constant(10);% c3 
constant(14)= constant(7)/constant(10);% c4 
constant(15)=(constant(6)-constant(4))/constant(5);% c5 
constant(16)= constant(7)/constant(5);% c6
constant(17)= 1/constant(5);% c7
constant(18)=  (constant(4)*(constant(4)-constant(5))+constant(7)*constant(7))/constant(10);% c8 
constant(19)=  constant(4)/constant(10);% c9 

% aircraft geometry
constant(20)=S_w;    % S, wing area, ft^2
constant(21)=c_w;    % cbar, mean geometric chord, ft
constant(22)=b_w;     % b, wing span, ft
constant(23)=0;      % phiT, thrust inclination angle, RADIANS
constant(24)=0;      % dT, thrust offset distance, ft

% Nondimensional Aerodynamic stability and control derivatives

% Drag Polar CD=k(CLstatic-CLdm)^2 + CDm
constant(25)=Cd_0;   % CDm, CD for minimum drag
constant(26)=k; % k
constant(27)=0;      % CLdm, CL at the minimum drag point

% Lift Force
constant(28)=CL_0(S_w,S_h,M,tc_w,alpha_0,epsilon_t,i_w,i_h,epsilon_0_h,AR_w,Lambda_c4,Lambda_c2,lambda_w,kappa,beta,b_w,d,AR_h,eta_h);  % CL0
constant(29)=CL_alpha(AR_w,AR_h,Lambda_c2,lambda_w,l_h,h_h,b_w,d,eta_h,S_h,S_w,kappa_h,Lambda_c2_h,beta,kappa);   % CLalpha
constant(30)=CL_de(S_w,S_h,AR_h,ce_ch,eta_oe,eta_ie,beta,kappa_h,lambda_h,Lambda_c2_h,tc_h,delta_e,Cl_alpha_h);   % CLdeltaE
constant(31)=CL_alpha_dot(l_h, h_h, b_w, lambda, AR_w, AR_h, Lambda_c4, Lambda_c4_h, beta, kappa, kappa_h, V_h, eta_h);    % CLalphadot
constant(32)=CL_q(Xw,b_w,c_w,c_h,AR_w,Lambda_c4,Lambda_c2,Lambda_c2_h,Xh,S_h,S_w,eta_h,AR_h,beta, V_h,b_h, kappa, kappa_h);    % CLQ
% Side Force
constant(33)=0;      % CY0
constant(34)=Cy_beta(two_r_one,eta_v,beta,AR_v,b_v,Z_h,x_over_c_v,lambda_v,S_v,S_w,Lambda_c4,Lambda_c4_v,Z_w,d,dihedral,wingloc,Z_w1,S_h,S_o); % CYbeta
constant(35)=Cy_da(S_w);      % CYdeltaA
constant(36)=Cy_dr(S_w,b_w,S_h,S_v,b_v,c_v,x_AC_vh,two_r_one,AR_v,l_v,Z_v,eta_or,eta_ir,cr_cv,beta,kappa_v,Lambda_c2_v,lambda_v,delta_r,alpha);  % CydeltaR
constant(37)=Cy_p(b_w,l_v,Z_v,alpha,two_r_one,eta_v,beta,AR_v,b_v,Z_h,x_over_c_v,lambda_v,S_v,S_w,Lambda_c4_v,Z_w,d,dihedral,wingloc,Z_w1,S_h,S_o); % Cy_p
constant(38)=Cy_r(alpha,two_r_one,eta_v,beta,AR_v,b_v,b_w,Z_h,x_over_c_v,lambda_v,S_v,S_w,Lambda_c4,Z_w,d,dihedral,wingloc,Z_w1,S_h,S_o,l_v,Z_v);  % CYR
% Rolling Moment
constant(39)=0;       % Cl0
constant(40)=Cl_beta(CL_wb,theta,theta_h,Lambda_c2, Lambda_c4,Lambda_c4_h,S_w,b_w,AR_w,AR_h,Z_v,alpha,l_v, M,two_r_one,eta_v,beta,AR_v,b_v,Z_h,x_over_c_v,lambda_v,S_v,Z_w,d,wingloc,Z_w1,S_h,S_o, dihedral,dihedral_h,Lambda_c2_h,l_b,b_h,CL_hb); % Clbeta
constant(41)=Cl_da(S_w,AR_w,ca_cw,eta_ia,eta_oa,beta,kappa,Lambda_c4,lambda_w,tc_w,Cl_alpha_w);   % CldeltaA
constant(42)=Cl_dr(S_w,b_w,S_h,S_v,b_v,c_v,x_AC_vh,two_r_one,AR_v,l_v,Z_v,eta_or,eta_ir,cr_cv,beta,kappa_v,Lambda_c2_v,lambda_v,delta_r,alpha);  % CldeltaR
constant(43)=Cl_p(b_h,b_w,Z_v,two_r_one,eta_v,AR_v,b_v,Z_h,x_over_c_v,lambda_v,S_v,S_w,Lambda_c4,Lambda_c4_h,Z_w,d,lambda_w,wingloc,Z_w1,S_h,S_o,beta,Cl_alpha,Cl_alpha_h,AR_w,lambda_h, AR_h);  % ClP
constant(44)=Cl_r(l_v,alpha,Z_v,S_h,S_v,x_over_c_v,b_v,Z_h,two_r_one,lambda_v,S_w,Z_w,d,AR_w,beta,Lambda_c4,c_w,Xw,AR_v,eta_v,b_w,cf,delta_f,theta,lambda_w,Gamma);  % ClR
% Pitching Moment
constant(45)=Cm_0(S_h_slip,S_w,S_h,M,tc_w,alpha_0,epsilon_t,i_w,i_h,epsilon_0_h,AR_w,Lambda_c4,lambda_w,beta,Cm_0_r,Cm_o_t,Lambda_c2_h,kappa_h,AR_h,Xw);   % CM0
constant(46)=Cm_alpha(AR_w,AR_h,Lambda_c2,lambda_w,l_h,h_h,b_w,c_w,d,eta_h,S_h,S_w,kappa_h,Lambda_c2_h,beta,kappa,Xw); % Cmalpha
constant(47)=Cm_de(S_w,S_h,AR_h,ce_ch,eta_oe,eta_ie,beta,kappa_h,lambda_h,Lambda_c2_h,tc_h,delta_e,Cl_alpha_h,V_h); % CMdeltaE
constant(48)=Cm_a_dot(l_h,h_h,b_w,lambda_w,AR_w,Lambda_c4,M,Cl_alpha,q_bar_h,q_bar,S_h,Xh,c_w,S_w);  % CMalphadot
constant(49)=Cm_q(Xw,b_w,c_w,c_h,AR_w,Lambda_c4,Lambda_c2,Xh,S_h,S_w,eta_h,AR_h,beta,V_h,b_h,Cl_alpha,B);  % CMQ
% Yawing Moment
constant(50)=0;         % CN0
constant(51)=Cn_beta(S_w,b_w,alpha,l_v,Z_v,l_f,S_b_s,Rl_f,x_m,h1_fuse,h2_fuse,hmax_fuse,wmax_fuse,two_r_one,eta_v,M,AR_v,b_v,Z_h,x_over_c_v,lambda_v,S_v,Z_w,d,Z_w1,S_h,Lambda_c4,Lambda_c2_v);  % CNbeta
constant(52)=Cn_da(S_w,AR_w,ca_cw,eta_ia,eta_oa,beta,kappa,Lambda_c4,lambda_w,tc_w,Cl_alpha_w,CL); % CNdeltaA
constant(53)=Cn_dr(S_w,b_w,S_h,S_v,b_v,c_v,x_AC_vh,two_r_one,AR_v,l_v,Z_v,eta_or,eta_ir,cr_cv,beta,kappa_v,Lambda_c2_v,lambda_v,delta_r,alpha); % CNdeltaR
constant(54)=Cn_p(c_w,B,theta,adelf,delf,b_w,l_v,b_h,Z_v,two_r_one,eta_v,AR_v,b_v,Z_h,x_over_c_v,lambda_v,S_v,S_w,Lambda_c4,Lambda_c2,Lambda_c4_h,Z_w,d,lambda_w,wingloc,Z_w1,S_h,S_o,beta,Cl_alpha,Cl_alpha_h,AR_w,lambda_h, AR_h, b_f,Xw,CL,Lambda_c4_v,alpha); % CNP
constant(55)=Cn_r(CL,C_bar_D_o,l_v,alpha,Z_v,S_h,S_v,x_over_c_v,b_v,Z_h,two_r_one,lambda_v,S_w,Z_w,d,AR_w,beta,Lambda_c4,c_w,Xw,AR_v,eta_v,b_w); % CNR
% Reference positions
constant(56)=.264;  % XbarRef, nondimensional
constant(57)=.264;  % XbarCG,  nondimensional

% Trim conditions. These may or not be used by subsequent programs. Small 
% variations in these trim flight conditions are OK.
constant(58)=U;     % trim speed, Vt, ft/sec
constant(59)=5000;  % Trim altitude, ft
constant(60)=0;     % Trim alpha, >>>DEGREES<<<This is not used by CessnaLongSC
% The constants below are used only by CessnaLongSC and CessnaLatSC
constant(61)=0;     % CLu=0
constant(62)=0;     % CDu=0
constant(63)=-.096; % CTxu
constant(64)=0;     % Cmu
constant(65)=0;     % CmTu
constant(66)=0;     % CmTalpha
constant(67)=0;     % CDdeltae

% #########################################################################

% *********************************************
% BasicConstants_Cessna182
%
% *********************************************
% BasicConstants - Identifies, describes, and assigns all of the 
%                  the most basic variables for analyzing the control
%                  and stability of a generic aircraft.
% *********************************************     
%
% A&AE 421 Fall 2001 - Purdue University
% 
% Note: This code is provided for a first order approximation of the dynamic 
%       analysis of an airplane and is not intended for final designs.
%
% Equations/Figures can be found in :
% 
% (Ref.1) Roskam, Jan. "Airplane Flight Dynamics and Automatic Flight
%         Controls"
%         Published by DARcorporation
%         120 E. Ninth St., Suite 2
%         Lawarence, KS 66044
%         Third Printing, 2001.
%
% (Ref.2) Roskam, Jan. "Methods for Estimating Satbility and
%         Control Derivatives of Conventional Subsonic Airplanes"
%         Published by the Author 
%         519 Boulder 
%         Laurance, Kansas 66044
%         Third Printing, 1997.
%  
% (Ref.3) Roskam, Jan. "Airplane Design: Part IV: Preliminary Calculation
%         of Aerodynamic, Thrust and Power Characteristics"
%         Published by Roskam Aviation and Engineering Corporation
%         Rt4, Box 274
%         Ottawa, Kansas 66067
%         Second Printing, 1990.
%

aircraft='Cessna182, cruise configuration';

adelf = 0;           % Two dimensional lift effectiveness parameter Ref.(2), Equ(8.7)alpha = 0;           % Angle of attack [deg]
alpha_0 = -.2;	    	% Airfoil zero-lift AOA [deg]
AR_h = 3.93;         % Aspect ratio of the horizontal tail
AR_w = 7.45;         % Aspect ratio of the wing
b_f = 20;            % Span of the flap [ft]
b_h = 11.33;     		% Span of the horizontal tail [ft]
b_h_oe = 5.8333;     % Elevator outboard position [ft]
b_h_ie = 0.5;     	% Elevator inboard position [ft]
b_w = 36;            % Span of the wing [ft]
b_v = 5;             % Vertical tail span measured from fuselage centerline [ft]
b_v_or = 4.95;       % Outboard position of rudder [ft]
b_v_ir = 0.0;        % Inboard position of rudder [ft]
c_a = 1.0;           % Chord of aileron [ft]
C_bar_D_o = .0270;   % Parasite drag
Cd_0 = 0.027;        % Drag coefficient at zero lift (parasite drag)
c_e = 2.70;	         % Elevator flap chord [ft]
cf = 1.2;            % Length of the wing flap [ft]
c_h = 6.3;           % Mean aerodynamic chord of the horizontal tail [ft]
CL = 0.307;          % Lift coefficient (3-D)
CL_hb=.307;          % Lift coefficient of the horzontal tail/body
CL_wb=.307; 	      % Lift coefficient of the wing/body
Cl_alpha_h  = 2*pi;	% 2-D Lift curve slope of wing
Cl_alpha_v	= 2*pi;  % 2-D Lift curve slope of vertical tail
Cl_alpha = 6;        % Two-dimensional lift curve slope
Cl_alpha_w=Cl_alpha; % Two-dimensional lift curve slope
Cm_0_r = -.2;        % Zero lift pitching moment coefficient of the wing root
Cm_o_t = -.2;        % Zero lift pitching moment coefficient of the wing tip
c_r = 1.52;          % Chord of the rudder [ft]
c_w = 4.9;           % Mean aerodynamic chord of the wing [ft]
c_v = 3.8;           % Mean aerodynamic chord of the vertical tail [ft]
D_p = 4;			   	% Diamter of propeller [ft]
d = 4.5;					% Average diameter of the fuselage [ft]
delf = 0;            % Streamwise flap deflection [deg] 
delta_e = 0;	   	% Elevator deflection [deg]
delta_f = 0;         % Streamwise flap deflection [deg]
delta_r = 0;         % Rudder deflection [deg]
dihedral = 0;        % Geometric dihedral angle of the wing [deg]
dihedral_h = 0;      % Geometric dihedral angle of the horizontal tail [deg]
e = 0.9;             % Oswald efficiency factor
epsilon_t = -.421;	% Horizontal tail twist angle [deg]
epsilon_0_h=epsilon_t;% Horizontal tail twist angle [deg]
eta_h = 0.85;	      % Ratio of dynamic pressure at the tail to that of the free stream eta_ia = .06;        % Percent span position of inboard edge of aileron
eta_oa = 1.0;        % Percent span position of outboard edge of aileron
eta_p = 0.70;			% Propeller Efficiency
eta_v = 1;				% Ratio of the dynmaic pressure at the vertical tail to that of the freestream
Gamma = 2*pi/180;    % This is the geometric dihedral angle, positive for dihedral, negative for anhedral [rad]
h1_fuse =4;          % Heighth of the fuselage at 1/4 of the its length 
h2_fuse =3;          % Heighth of the fuselage at 3/4 of the its length 
h_h = 1.7;           % Height from chord plane of wing to chord plane of horizontal tail [ft] - Fig 3.7, Ref. 2
hmax_fuse = 5;       % Maximum heighth of the fuselage [ft]  
Ixx = 948;
Iyy = 1346;
Izz = 1967;
Ixz = 0;
i_h = 2;             % Horiz tail incidence angle [deg]  
i_w	= -.2;			% Wing incidence angle [deg]  
k = 0.0554;				% k
Lambda=0;            % Wing sweep angle [deg]
Lambda_c4=0;         % Sweep angle at the quarter-chord of the wing [deg]      
Lambda_c2=0;         % Sweep angle at the half-chord of the wing [deg]
Lambda_c2_v = 0.35;  % Sweep angle at the half-chord of the vertical tail [deg]
Lambda_c4_v = 0.35;  % Sweep angle at the quarter-chord of the vertical tail [deg]
Lambda_c2_h = 0;     % Sweep angle at the half-chord of the horizontal tail [deg]
Lambda_c4_h = 0;     % Sweep angle at the quarter-chord of the vertical tail [deg]
lambda = 0.75;       % Taper ratio of the wing
lambda_h = 0.75;  	% Horizontal tail taper ratio
lambda_v	= 0.5;      % Vertical tail taper ratio
lambda_w = lambda;   % Taper ratio of the wing
l_b = 25;	         % This is the total length of the fuselage [ft]
l_f =12;             % The horizontal length of the fuselage [ft]  (*******)
l_h = 12;            % Distance from c/4 of wing to c/4 of horizontal tail [ft] 
l_v = 13;    			% Horizontal distance from the aircraft CG to the vertical tail aero center [ft]
M = 0.2;             % Mach number
q_bar = 1; 				% Dynamic pressure ratio (free stream)
q_bar_h = 1;         % Dynamic pressure ratio at the tail
Rl_f=10^6;           % Fuselage Reynold's number
rho = 0.00204;			% Air density at 5000ft [slugs/ft^3]
S_b_s =60;           % Body side area [ft^2]
S_b = S_b_s;         % Body side area [ft^2] 
S_h	= 48;				% Area of the horizontal tail [ft^2]
S_h_slip	= 8;	   	% Area of horiz. tail that is covered in prop-wash - See fig.(8.64) Ref.(3) [ft^2]
S_o = 6;	            % Take X1/l_b and plug into: .378+.527*(X1/l_b)=(Xo/l_b), S_o is the fuselage x-sectional area at Xo. (ft^2) Ref.(2), Fig. 7.2  
S_w = 174;           % Aera of the wing [ft^2]
S_v = 25;   			% Surface area of vertical tail [ft^2]
T = 0;					% Temperature [F]
tc_w = .10;         	% Thickness to chord ratio of the wing 
tc_h = 0.10;     	   % Thickness to chord ratio of horizontal tail
theta = 0;           % This is the wing twist in degrees, negative for washout [deg]
theta_h = 0;         % Horizontal tail twist between the root and tip stations, negative for washout [deg]
two_r_one = 2;       % Fuselage depth in region of vertical tail [ft] Ref.(2), Figure 7.5
U = 129.73;				% Free Stream Velocity [knots]
U1  = 220.1;			% Cruise flight speed [ft/s]
W = 2650;            % Weight of Airplane [lbf]
wingloc = 1;         % If the aircraft is a highwing: (wingloc=1), low-wing: (wingloc=0) 
wmax_fuse = 5;       % Maximum fuselage width [ft]
X1 = 12.5;	         % Distance from the front of the fuselage where the x-sectional area decrease (dS_x/dx) is greatest (most negative) [ft] - Ref.(2), Fig. 7.2  
x_AC_vh =1;          % X distance from LE of vertical tail to AC of horizontal tail [ft]
Xh = 18;             % Distance from airplane cg to the horizontal tail ac [ft] 
x_m = 5;             % Center of gravity location from the leading edge [ft]
x_over_c_v = .5;     % X distance from leading edge of vertical fin mean chord to horizontal aero center [ft] 
Xref = 6;            % Arbitrary reference point located on the airplane's axis of symmetry. 
							% Measured as positive aft, starting from the leading edge of the mean aero. chord. [ft]
Xw = 0.2;            % Distance from the airplane cg to wing ac (positive aftward) [ft]
Z_h = -2.5;      		% Negative of the hoizontal distance from the fuselage centerline to the horizontal tail aero center (negative number) - Ref.(2), Fig. 7.6
Z_v = .5; 				% Vertical distance from the aircraft CG to the vertical tail aero center - Ref.(2), Fig. 7.18
Z_w = 2;					% This is the vertical distance from the wing root c/4 to the fuselage centerline, positive downward - Ref.(2), Equ(7.5)
Z_w1 = -.5;	         % Distance from body centerline to c/4 of wing root chord, positive for c/4 point below body centerline (ft) - Ref.(2), Fig. 7.1  



% ******************************************************************
% DerivedConstants - Identifies, describes, and assigns all of the 
%                  the basic derived constants used for analyzing the 
%                  control and stability of a generic aircraft.
% ******************************************************************     
%
% A&AE 421 Fall 2001
% 
% Equations/Figures can be found in :
%
% (Ref.1) Roskam, Jan. "Airplane Flight Dynamics and Automatic Flight
%         Controls"
%         Published by DARcorporation
%         120 E. Ninth St., Suite 2
%         Lawarence, KS 66044
%         Third Printing, 2001.
%
% (Ref.2) Roskam, Jan. "Methods for Estimating Satbility and
%         Control Derivatives of Conventional Subsonic Airplanes"
%         Published by the Author 
%         519 Boulder 
%         Laurance, Kansas 66044
%         Third Printing, 1997.
%  
% (Ref.3) Roskam, Jan. "Airplane Design: Part IV: Preliminary Calculation
%         of Aerodynamic, Thrust and Power Characteristics"
%         Published by Roskam Aviation and Engineering Corporation
%         Rt4, Box 274
%         Ottawa, Kansas 66067
%         Second Printing, 1990.
%

AR_v = b_h^2/S_v;                 % Aspect Ratio of Horizontal Tail
AR_h = b_h^2/S_h;					    % Aspect Ratio of Horizontal Tail
AR_w = b_w^2/S_w;                 % Aspect Ratio of Horizontal Tail
B=sqrt(1-M^2*(cos(Lambda_c4))^2); % Compressibility correction factor
beta = sqrt(1-M^2);               % Compressibility correction factor
Beta = beta;                      % Compressibility correction factor
ca_cw = c_a/c_w;                  % Ratio of aileron chord to wing chord
Cd = Cd_0 + (CL^2/(pi*AR_w*e));   % Drag Coefficient
ce_ch = c_e/c_h;				     	 % Elevator flap chord 
cr_cv			= c_r/c_v;            % Rudder flap chord 
eta_ie = b_h_ie/(b_h/2);	  	    % Percent span position of inboard edge of elevator
eta_ir		= b_v_ir/b_v;         % Percent span position of inboard edge of rudder
eta_oe = b_h_oe/(b_h/2);		    % Percent span position of outboard edge of elevator
eta_or		= b_v_or/b_v;         % Percent span position of outboard edge of rudder
kappa=Cl_alpha/(2*pi);            % Ratio of 2D lift coefficient to 2pi for the wing
kappa_h = Cl_alpha_h/(2*pi);	    % Ratio of 2D lift coefficient to 2pi for the horiz. tail
kappa_v		= Cl_alpha_v/(2*pi);  % Ratio of 2D lift coefficient to 2pi for the vert. tail
V_h = (Xh*S_h)/(c_h*S_w);         % Horizontal Tail Volume Coefficient
