% *********************************************
% Make_Constant    Version 9.1 3/23/05
% This version uses the new Cl_beta code from Barnes McCormick
%
% OBJECTIVES: 1. Determine a set of derived constants that are 
%                used by the stability and control derivative functions.
%             2. 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.
%
% INPUTS: Many basic constants defined in the Matlab workspace.
% OUTPUTS: The array called constant of length 67.
% IMPORTANT NOTE: This script is INDEPENDENT of the particular aircraft being modeled.
%
%
% This version corrects Cl_beta, Cy_beta, and Cn_beta.
% This version corrects Cm_alpha, rudder control power, units, and comments.
% AAE565 Spring 2003
%
% *********************************************     
%
% A&AE 565 Spring 2003 - 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 VI: Preliminary Calculation
%         of Aerodynamic, Thrust and Power Characteristics"
%         Published by Roskam Aviation and Engineering Corporation
%         Rt4, Box 274
%         Ottawa, Kansas 66067
%         Second Printing, 1990.
%
% Ref: Barnes W. McCormick, Aerodynamics, Aeronautics, and Flight mechanics, Second Edition
%      John Wiley and Sons, Inc. 1995, pages 529-534 (used fror Cl_beta_Mc)

%disp(' '); disp(' Starting Make_Constant: The array called constant is now being defined.'); disp(' ')
% The derived constants below are defined from the basicConstants defined
% in the BasicConstants script.
% AR_v = b_v^2/S_v;               % Aspect Ratio of Vertical 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
% Xh=Xach-Xref;         			% Distance from airplane arbitrary reference point to the horizontal tail ac [ft]
% V_h = (Xh*S_h)/(c_h*S_w);       % Horizontal Tail Volume Coefficient.
% Xw=Xacw-Xref;					% Distance from arbitrary ref point to aero center of wing alone. [ft]
% l_h=Xach-.25*c_w;				% Distance from AC of Horizontal tail and wing c/4 (ft), positive aft.
% delta_f = delf;         		% Streamwise flap deflection [rad] - same as delf
% Gamma = dihedral;				% This is the geometric dihedral angle, positive for 
% %								dihedral, negative for anhedral [rad]
% lambda_w = lambda;     			% Taper ratio of wing again
% S_b = S_b_s;           			% Body side area [ft^2] again
% x_AC_vh=x_over_c_v*c_v;     	% X distance from LE of vertical tail to AC of horizontal tail [ft]
% U1 = U*1.687808;   	    		% Trim flight speed [ft/s]
% disp(['Trim airspeed= ',num2str(U1),' ft/sec'])
% nu=nufun(altitude);             % Kinematic viscosity, ft*ft/sec
% rho=rhofun(altitude);			% Air density at altitude [slug/ft^3]
% Rl_f= U1*l_f/nu;    			% Reynolds number of fuselage [non-dimensional]
% l_b = l_f;	           			% Total length of fuselage [ft]
% End of the derived constants section.

% echo on
% % 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 57 "constant". These MATLAB programs perform 6-degree-of-freedom 
% % aircraft simulation over a flat earth or an elliptical earth.
% % Two other programs (LongSC.m and LatSC.m) do simplified dynamic modeling (compute transfer functions)
% % for the longitudinal or lateral-directional degrees of freedom.  They addtionally use all
% % constant(58)-constant(67) of the "constant" defined below. 
% %
% echo off
% Mass related inputs
constant(1)=0.4126;      	% 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)=0.003106237;    	% Ixx, slug*ft*ft
constant(5)=0.001787763;   		% Iyy, slug*ft*ft
constant(6)=0.004624871;   		% Izz, slug*ft*ft
constant(7)=0.000081565;    	% Ixz, slug*ft*ft
constant(8)=0.5744;  	% 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)=1.37562;    % S, wing area, ft^2
constant(21)=0.184;    % cbar, mean geometric chord, ft
constant(22)=27.24/12;    % 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)=0.0385;   	% CDm, CD for minimum drag
constant(26)=0.0332; 		% k
constant(27)=0.0063;      	% CLdm, CL at the minimum drag point

disp(' '); disp(' ');
% disp('Below this point in Make_Constant there may be multiple occurrances of warnings like the following')
% disp('    Warning: 0.25 exceeded interpolation limits')
% disp('    Warning: 0.25 was converted to 0.5 (plus some traceback information)')
% disp('Each warning means that a variable was found to be out of limits')
% disp('for using an empirical data found in a figure. No extrapolation was performed.')
% disp('Instead, the variable was limited to the value shown and the program proceded.')
% disp('The careful user will check the that figure in the appropriate reference for each warning to see ')
% disp('if the action of the computer was reasonable.')
% disp(' ')

% Lift Force
constant(28)=0.5082;   %CL_0
constant(29)=4.977;   % CLalpha
constant(30)=0.0924;   % CLdeltaE
constant(31)=1.102935989;    % CLalphadot
constant(32)=1.951635944;    % CLQ
% Side Force
constant(33)=0;      % CY0
constant(34)=-0.341381675; % CYbeta
constant(35)=0;      % CYdeltaA
constant(36)=0;  % CydeltaR
constant(37)=0; % Cy_p
constant(38)=0.173001337;  % CYR
% Rolling Moment
constant(39)=0;       % Cl0

constant(40)=-0.175871037;% Clbeta
%disp(' '); disp('Cl_beta using the original code from Roskam, per radian')
%constant(40)
% disp(' '); disp('Cl_beta per degree broken down into individual components using the new code from McCormick')
% constant(40)=Cl_beta_Mc(low_wing,Cl_alpha_v,S_v,Z_v,S_w,b_w,AR_v,AR_w,Lambda_c2,lambda_w,dihedral,W,U1,rho); % Clbeta
% disp(' '); disp('Cl_beta using the new code from McCormick, per radian')
% constant(40)
constant(41)=0.234568921;   % CldeltaA
constant(42)=0;  % CldeltaR
constant(43)=-.340164879;  % ClP
constant(44)=0;  % ClR
% Pitching Moment
constant(45)=-0.0094;   % CM0
constant(46)=-0.89586; % Cmalpha
constant(47)=-0.0447; % CMdeltaE
constant(48)=-3.5722;  % CMalphadot
constant(49)=-6.609;  % CMQ
% Yawing Moment
constant(50)=0;         % CN0
constant(51)=0.0846105;  % CNbeta
constant(52)=-0.074853516; % CNdeltaA
constant(53)=0; % CNdeltaR
constant(54)=0; % CNP
constant(55)=-0.084117928; % CNR
% Reference positions
constant(56)=0.5793;  % XbarRef, nondimensional, measured aft from leading edge of wing mean aerodynamic chord.
constant(57)=0.388276;  	% XbarCG,  nondimensional, measured aft from leading edge of wing mean aerodynamic chord.
% constant(57 can be changed at a later time if you use program FlatEarth or E_Earth.

% Trim conditions. These may or not be used by subsequent programs LongSC.m and LatSC.m. Small 
% variations in these trim flight conditions are OK.
constant(58)=17;    	% Trim speed, Vt, ft/sec aka U1
constant(59)=10;  % Trim altitude, ft
constant(60)=6;    		% Trim alpha, >>>DEGREES<<<This is not used by LongSC
% 						The constants below are used only by LongSC.m and LatSC.m MATLAB scripts for liear small perturbation
% 						longitudinal and lateral-directional analysis as defined in Reference 1.
constant(61)=0;     	% CLu
constant(62)=0;     	% CDu
constant(63)=0; 		% CTxu
constant(64)=0;     	% Cmu
constant(65)=0;     	% CmTu
constant(66)=0;     	% CmTalpha
constant(67)=0;     	% CDdeltae

disp(' '); disp(' ');
% disp('Above this point in Make_Constant there may be multiple occurrances of warnings like the following')
% disp('    Warning: 0.25 exceeded interpolation limits')
% disp('    Warning: 0.25 was converted to 0.5 (plus some traceback information)')
% disp('Each warning means that a variable was found to be out of limits')
% disp('for using an empirical data found in a figure. No extrapolation was performed.')
% disp('Instead, the variable was limited to the value shown and the program proceded.')
% disp('The careful user will check the that figure in the appropriate reference for each warning to see ')
% disp('if the action of the computer was reasonable.')
% disp(' '); disp('Done with Make_Constant: The array called constant is now defined.'); disp(' ')

