% This script computes basic lift and pitching momewnt data using the
% Roskam methods

disp(' '); disp(' >>>>>>>>>Start here <<<<<<<'); disp(' ')
disp(' INPUTS')
clear
format compact
% Basic constants from which the stability and control derivatives are computed

% Properties of the aircraft
% M             - Mach number
M=.3
% h_h           - the vertical distance from the wing root trailing edge to the chord line of the horizontal tail (ft)
h_h=0
% d             - Fuselage diameter (ft)
d=5
% epsilon_0_h 	- Downwash angle of tail at zero AOA [rad]
epsilon_0_h=0
% eta_h         - Ratio of dynamic pressure at the tail to that of the free stream
%				  usually assumed to be in the range of .90<eta_h<1.0
eta_h=1
% Xref			Distance from the leading edge of the wing mean aerodynamic chord
%                               to the arbitrary reference point. The equivalent force system
% 								for the aerodynamic force system is given about this point.
%						        Measured as positive aft, starting from the leading edge of the mean aero. chord. [ft]
Xref=.264*4.9
% Xacw			Distance from the leading edge of the wing mean aerodynamic chord
%                               to the aerodynamic center of the wing. 
%						        Measured as positive aft, starting from the
%						        leading edge of the mean aero. chord. [ft]
Xacw=.25*4.9
% Xacwb			Distance from the leading edge of the wing mean aerodynamic chord
%                               to the aerodynamic center of the wing and body. 
%						        Measured as positive aft, starting from the
%						        leading edge of the mean aero. chord. [ft]
Xacwb=.25*4.9

% Xach			Distance from the leading edge of the wing mean aerodynamic chord
%                               to the aerodynamic center of the horizontal tail. 
%						        Measured as positive aft, starting from the
%						        leading edge of the mean aero. chord. [ft]
Xach=3*4.9

% Properties of the Wing
% S_w           - Area of the wing (ft^2)
S_w=174
% b_w           - span of the wing (ft)
b_w=36
% c_w           - chord of the wing (ft)
c_w=4.9
% Lambda_c4	    - Wing quarter chord sweep angle [rad]
Lambda_c4=0
% Lambda_c2     - Wing Half chord sweep angle (rad)
Lambda_c2=0
% lambda_w      - Wing Taper ratio
lambda_w=.9
% tc_w			- Wing thickness to chord ratio
tc_w=.15
% alpha_0		- Wing airfoil zero-lift AOA [rad] 2D or 3D?????
alpha_0=0
% i_w			- Wing incidence angle [rad]
i_w=0
% Cl_alpha      - Wing lift curve slope
Cl_alpha=.9*2*pi

% Properties of the horizontal tail
% S_h           - Area of the horizontal tail (ft^2)
S_h=41.4
% S_h_slip		- Area of horizontal tail that is covered in prop-wash - see fig.(8.64) Ref.(3) [ft^2]
S_h_slip=41.4
% c_h           - Horizontal tail chord (ft)
c_h=2.5
% b_h           - Span of the horizontal tail (ft)
b_h=11+8/12
% Lambda_c2     - Horizontal tail half chord sweep angle (rad)
Lambda_c2_h=0
% lambda_h		- Horizontal tail taper ratio
lambda_h=.9
% epsilon_t		- Horizontal tail twist angle [rad]
epsilon_t=0
% i_h			- Horizontal tail incidence angle [rad]
i_h=0
% tc_h			- Thickness to chord ratio of horizontal tail
tc_h=.15
%		delta_e 		- Elevator deflection [rad]
delta_e=0
% Cl_alpha_h    - 2D lift curve slope of the airfoil of the horizontal tail
Cl_alpha_h=.9*2*pi



% Properties of the elevator 
% c_e           - Elevator chord (ft)
c_e=1.25

b_h_ie=0

b_h_oe=1

% delta_e 		- Elevator deflection [rad]
delta_e=0




% Cm_0_wf         The zero-lift pitching moment coefficient of the wing-fuselage combination (Pg. 320 of Ref(3))
% Cm_0_c          The zero-lift pitching moment coefficient of the canard
% Cm_0_h          The zero-lift pitching moment coefficient of the horizontal tail
% Cm_0_w          The zero-lift pitching moment coefficient of the wing
% Cm_0_f          The zero-lift pitching moment coefficient of the fuselage
Cm_0_r=0
Cm_o_t=0



%		ce_ch			- Ratio of elevator flap chord to horizontal tail chord
%		eta_oe		- Percent span position of outboard edge of elevator
%		eta_ie		- Percent span position of inboard edge of elevator
%		beta			- Compressibility correction factor (sqrt(1-Mach^2))
%		Lambda_c2_h	- Horizontal tail sweep angle at 1/2chord [rad]
%		lambda_h		- Horizontal tail taper ratio
%		tc_h			- Thickness to chord ratio of horizontal tail

%		Cl_alpha_h	- Horizontal tail 2D lift curve slope
%		V_h			- Horizontal tail volume coefficient

% Derived constants
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
ce_ch = c_e/c_h;				% Elevator flap chord to horizontal tail chord ratio 
eta_ie = b_h_ie/(b_h/2);	  	% Percent span position of inboard edge of elevator
eta_oe = b_h_oe/(b_h/2);		% Percent span position of outboard edge of elevator
kappa=Cl_alpha/(2*pi);          % Ratio of 2D lift curve slope to 2pi for the wing
kappa_h = Cl_alpha_h/(2*pi);	% Ratio of 2D lift curve slope to 2pi for the horiz. 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.
% End of the derived constants section.

%tempprint(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)
%t1=[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]


disp(' '); disp(' OUTPUTS'); disp(' ')
% Compute stability and control derivatives
% Lift 

%[CL_0f]=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)
[CL_0c]=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)
[CL_alphac]=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)
[CL_dec]=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)
                   
% Pitching moment where Xref is the moment reference point
[Cm_0c]=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,Xh,c_w,eta_h)
[Cm_alphac]=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,Xref,Xacwb)
[Cm_dec]=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)
%t2=[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]
%tdif=t2-t1

%
% Applicable References (NB: the first number in parentheses is the number by 
% which the given reference will be referred to in the code)
% (3)[Ref. ] Roskam's Airplane Design Series: Book VI


