% Patrick LeGresley
% Angel L. carrion
%
% Class II Weight Calculations
%
% March 2, 1999
%
% Last modified April 24, 1999



function [W_to,W_e,W_eng_single,x_cg_aft,x_cg_forward,z_cg_high,I_yy] =...
   wbfunction(tiltwing,b_wing,S_ref,taper_wing,sweep_wing_25c,tc_avg_wing,x_wing,...
   b_ht,Sht_ref,taper_ht,sweep_ht_25c,tc_avg_ht,x_ht,b_vt,Svt_ref,taper_vt,sweep_vt_25c,...
   x_vt,h,V,T_to_uninst,N_e,x_nosegear,x_maingear,W_to,W_payload,W_e,n_ult,V_C,V_D,V_H,D_eng,L_eng)

P_to_uninst = T_to_uninst;
N_p = 2;
N_bl = 3;
D_p = (b_wing/2) - 11;

% The moments of inertia are about the same point where the CG's were calculated about (the nose)
% tiltwing = 1 makes computations for tiltwing version

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define inputs

commercial = 0;                                 % commerical = 1 makes computations for commercial version

numbcomps = 21;                                 % number of components in weight breakdown

W_tfo = 0.005*W_to;
W_crew = 400;



% Wing inputs
AR_wing = b_wing^2/S_ref;                       % Wing aspect ratio (between 4 and 12 for GD equation)
tc_max_wing = 1.2* tc_avg_wing;                 % maximum wing thickness ratio

% wing sweep at leading edge (radians):
sweep_wing_le = 180/pi * atan(tan(sweep_wing_25c*pi/180) + 4*0.25*(1-taper_wing)/(1+taper_wing)/AR_wing);
% wing semi-chord sweep angle (deg):
sweep_wing_50c = 180/pi * atan(tan(sweep_wing_le*pi/180) - 4*0.5*(1-taper_wing)/(1+taper_wing)/AR_wing); % degrees

c_r_wing = 2*S_ref/(b_wing*(1+taper_wing));
c_t_wing = taper_wing * c_r_wing;
t_r_wing = tc_max_wing * c_r_wing;              % wing root thickness

% wing mean aerodynamic chord:
mac_wing = 2/3*(c_r_wing + c_t_wing - (c_r_wing*c_t_wing)/(c_r_wing + c_t_wing));
y_mac_wing = (mac_wing - c_r_wing)/((c_t_wing - c_r_wing) * 2/b_wing);


c_r_wing = 2*S_ref/(b_wing*(1+taper_wing));
c_t_wing = taper_wing*c_r_wing;
n_mac_wing = tan(sweep_wing_le/57.3)*y_mac_wing;
x_ac_wing = x_wing + n_mac_wing + 0.25*mac_wing;



W_fuel = W_to - W_e - W_crew - W_tfo - W_payload;						% fuel weight
W_MZF = W_to - W_fuel;									% maximum zero fuel weight




% Horizontal tail inputs
AR_ht = b_ht^2/Sht_ref;                         % horiz tail aspect ratio
c_r_ht = 2*Sht_ref/(b_ht*(1+taper_ht));		   % h tail root chord, ft
t_r_ht = 1.2*tc_avg_ht*c_r_ht;                  % h. tail maximum root thickness in feet
c_t_ht = c_r_ht * taper_ht;

K_h = 1.1;                                      % K_h = 1.0 for fixed incidence stabilizers
																% K_h = 1.1 for variable incidence stabilizers
                                                
sweep_ht_le = 180/pi * atan(tan(sweep_ht_25c*pi/180) + 4*0.25*(1-taper_ht)/(1+taper_ht)/AR_ht);
% horizontal tail semichord sweep angle:
sweep_ht_50c = 180/pi * atan(tan(sweep_ht_le*pi/180) - 4*0.5*(1-taper_ht)/(1+taper_ht)/AR_ht); % degrees


mac_ht = (2/3)*c_r_ht*((1+taper_ht+taper_ht^2)/(1+taper_ht));
y_mac_ht = b_ht*(1+2*taper_ht)/6/(1+taper_ht);
n_mac_ht = tan(sweep_ht_le/57.3)*y_mac_ht;
x_ac_ht = x_ht + n_mac_ht + 0.25*mac_ht;
x_h = x_ac_ht - x_ac_wing;


% Vertical tail inputs
N_vt = 2;                                       % Number of vertical tails
Svt_ref = 0.5*Svt_ref;
AR_vt = b_vt^2/Svt_ref;                         % vertical tail aspect ratio
c_r_vt = 2*Svt_ref/(b_vt*(1+taper_vt));
t_r_vt = 1.2*tc_avg_ht*c_r_vt;                  % v. tail maximum root thickness in feet
c_t_vt = c_r_vt * taper_vt;
z_h = 0;                                        % distance from v. tail root to where h. tail is mounted on v. tail
																% for fuselage mounted horizontal tails, set z_h = 0

S_r = 0.3*Svt_ref;                              % rudder area in sq ft

sweep_vt_le = 180/pi * atan(tan(sweep_vt_25c*pi/180) + 4*0.25*(1-taper_vt)/(1+taper_vt)/AR_vt);
% horizontal tail semichord sweep angle:
sweep_vt_50c = 180/pi * atan(tan(sweep_vt_le*pi/180) - 4*0.5*(1-taper_vt)/(1+taper_vt)/AR_vt); % degrees

% K_v = 1+0.15*(Sht_ref*z_h/(Svt_ref*b_vt)) for fin mounted horizontal tails
K_v = 1;                                        % K_v = 1.0 for fuselage mounted horizontal tails


mac_vt = (2/3)*c_r_vt*((1+taper_vt+taper_vt^2)/(1+taper_vt));
z_mac_vt = b_vt*(1+2*taper_vt)/6/(1+taper_vt);
n_mac_vt = tan(sweep_vt_le/57.3)*z_mac_vt;
x_ac_vt = x_vt + n_mac_vt + 0.25*mac_vt;

x_v = x_ac_vt - x_ac_wing;

% Fuselage inputs
w_max = 13;
l_nose = 13;												% length of nose
l_cargo = 27;												% cargo bay length
l_fus = 60;                                     % length of fuselage
w_fus = 13;                                     % maximum fuselage width
h_fus = 9.5;                                    % maximum fuselage height
K_inl = 1.0;                                    % factor for inlet location (see note under equation)
qbar_D = 136;                                   % design dive dynamic pressure in psf
K_fus = 1.4;                                    % factor for pressurization, gear location, cargo floor
																% see explanation with equation

S_fgs = 1700;                                   % fuselage gross shell area (approximately equal to fuselage wetted area)
qbar_L = qbar_D;                                % design dive dynamic pressure in psf
z_cg_fix = 0;                             % to move coord-sys to x-dir line through nose


% Nacelle inputs
N_inl = N_e;                                    % Number of inlets
A_inl = 0.9*pi*D_eng^2/4;                       % capture area per inlet
l_n = 1;                                        % nacelle length from inlet lip to compressor face
P2 = 30;                                        % maximum static pressure at engine compressor face in psi
																% typically 15 to 50 psi


% Landing gear inputs
z_ground2fus = 2;                               % distance from ground to bottom of fuselage, ft
l_s_m = 2.3;                                    % shock strut length for main gear in feet
W_L = 0.85*W_to;                                % landing weight
n_ult_L = 5.7;                                  % ultimate load factor for landing
K_g_r = 1.08;                                   % K_g_r = 1.0 for low wing airplanes
																% K_g_r = 1.08 for high wing airplanes

A_g_main = 40;                                  % constants A_g through D_g are defined in Table 5.1 of A/P Design Part V by Roskam
B_g_main = 0.16;
C_g_main = 0.019;
D_g_main = 1.5*10^-5;
A_g_nose = 20;
B_g_nose = 0.10;
C_g_nose = 0.0;
D_g_nose = 2.0*10^-6;


% Powerplant inputs
K_fsp = 6.55;                                   % JP-4 fuel: 6.55 lbs/gal
N_t = 2;                                        % Number of separate fuel tanks
fuelflow_TO = 0.5*T_to_uninst/3600;             % takeoff fuel flow in lbs/sec


%N_p = 4;                                       % number of propellers
%N_bl = 6;                                      % number of blades per propeller
%D_p = 15;                                      % propeller diameter
%P_to = 10500;                                  % takeoff power (hp)


% Fixed equipment inputs
K_fc = 0.64;                                    % K_fc = 0.44 for airplanes with unpowered flight controls
																% K_fc = 0.64 for airplanes with powered flight controls

K_hydr = 0.010;                                 % 0.0060 to 0.0120 for transports

V_pax = 2941;                                   % passenger cabin volume in cubic feet
N_pil = 2;                                      % number of pilots
R = 1500;                                       % range in nautical miles
N_cr = 3;                                       % number of crew members
N_pax = 24;                                     % number of passengers
N_pax_mil = 0;                                  % number of passengers in military configuration
l_pax = 26;                                     % length of passenger cabin
K_api = 887;                                    % for api weight calculations, see explanation with equation
V_pr = 5089;                                    % pressurized volume in cubic ft

K_apu = 0.008;                                  % typically 0.004 to 0.013

N_fdc = 2;                                      % number of flight deck crew
N_cc = 1;                                       % number of cabin crew
K_lav = 0.31;
K_buf = 1.02;
P_c = 12;                                       % design ultimate cabin pressure in psi

K_bc = 0.25;                                    % K_bc = 0.0646 without preload provisions
																% K_bc = 0.316 with preload provisions

K_pt = 0.0045;                                  % weight fraction of paint  (0.003 to 0.006)

W_e = 25500;                                    % airplane empty weight, first Guess

% end of inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[rho, Temp, a, mu] = atm(0);
M_h = 1.6878*V_H/a;


% Standard atmosphere calculations

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

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



% Class II Weights

weights = zeros(numbcomps,4);
iter = 0;
needtoiterate = 1;
W_to_ClassII = 0;

while abs(W_to - W_to_ClassII)/W_to > 0.001,
   iter = iter + 1;
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Wing Weights
   
   W_w_USAF = 96.948*((W_to*n_ult/10^5)^0.65*(AR_wing/cos(sweep_wing_25c/57.3))^0.57*(S_ref/100)^0.61*...
      ((1+taper_wing)/2/tc_max_wing)^0.36*(1+V_H/500)^0.5)^0.993;
   
   %W_w_GD = (0.00428*S_ref^0.48*AR_wing*M_h^0.43*(W_to*n_ult)^0.84*taper_wing^0.14)/...
   %   ((100*tc_max_wing)^0.76*cos(sweep_wing_50c/57.3)^1.54);
   
   W_w_Tor = 0.0017*W_MZF*(b_wing/cos(sweep_wing_50c/57.3))^0.75*(1+(6.3*cos(sweep_wing_50c/57.3)/b_wing))^0.5*...
      n_ult^0.55*(b_wing*S_ref/t_r_wing/W_MZF/cos(sweep_wing_50c/57.3))^0.3;
   
   % special notes for Torenbeek wing weight equation:
   % 2% should be added for spoilers and speed brakes
   % 2 wing mounted engines reduce wing weight by 5%
   % 4 wing mounted engines reduce wing weight by 10%
   % if LG not mounted to wing reduce the wing weight by 5%
   % for braced wings reduce the wing weight by 30% (does not include weight of stut)
   % weight of wing strut is approximately 10% of the wing weight
   % for fowler flaps add 2% to the wing weight
   
   %W_w = (W_w_USAF + W_w_GD + W_w_Tor)/3;
   W_w = (W_w_USAF + W_w_Tor)/2;
   
   % Wing folding stuff:
   W_w_fold = 0.12 * W_w;
   
   % Wing pivoting mechanism  (for fixed wing only)
   if tiltwing == 1,
      W_w_piv = 0.15 * W_w;
   else
      W_w_piv = 0;
   end
   
   W_w_tot = W_w + W_w_fold + W_w_piv;
   
   I_xx_on_wing = 0;		% Ixx for the component about its own cg
   I_yy_on_wing = 0;		% Iyy for the component about its own cg
   I_zz_on_wing = 0;		% Izz for the component about its own cg
   x_cg_wing = x_wing + 0.3 * mac_wing + y_mac_wing*tan(sweep_wing_le*pi/180);  % 30% chord at MAC
   z_cg_wing = h_fus + z_ground2fus - z_cg_fix + (tc_avg_wing * c_r_wing)/2;
   y_cg_wing = 0;
   
   weights(1,1) = W_w_tot;
   weights(1,2) = x_cg_wing * weights(1,1);
   weights(1,4) = z_cg_wing * weights(1,1);
   
   % end of wing weights
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Horizontal Tail Weights
   
   W_ht_USAF = 127*((W_to*n_ult/10^5)^0.87*(Sht_ref/100)^1.2*0.289*(x_h/10)^0.483*(b_ht/t_r_ht)^0.5)^0.458;
   
   W_ht_GD = 0.0034*((W_to*n_ult)^0.813*Sht_ref^0.584*(b_ht/t_r_ht)^0.033*(mac_wing/x_h)^0.28)^0.915;
   
   W_ht_Tor = K_h*Sht_ref*(3.81*(Sht_ref)^0.2*V_D/(1000*(cos(sweep_ht_50c/57.3)^0.5))-0.287);
   % K_h = 1.0 for fixed incidence stabilizers
   % K_h = 1.1 for variable incidence stabilizers
   
   W_ht = (W_ht_USAF + W_ht_GD + W_ht_Tor)/3;
   
   I_yy_on_ht = 0;		% Iyy for the component about its own cg
   I_xx_on_ht = 0;		% Iyy for the component about its own cg
   I_zz_on_ht = 0;		% Iyy for the component about its own cg

	x_cg_ht = x_ht + 0.3 * c_r_ht + 0.35*b_ht/2*tan(sweep_ht_le*pi/180);  % 30% chord at 35% semi-span
	y_cg_ht = 0;
   z_cg_ht = h_fus + z_ground2fus - z_cg_fix;
   
   weights(2,1) = W_ht;
   weights(2,2) = x_cg_ht * weights(2,1);
   weights(2,4) = z_cg_ht * weights(2,1);
   
   % end of horizontal tail weights
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Vertical Tail(s) Weights
   
   W_vt_USAF = 98.5*((W_to*n_ult/10^5)^0.87*(Svt_ref/100)^1.2*0.289*(b_vt/t_r_vt)^0.5)^0.458;
   
   W_vt_GD = 0.19*((1+z_h/b_vt)^0.5*(W_to*n_ult)^0.363*Svt_ref^1.089*M_h^0.601*x_v^-0.726*...
      (1+S_r/Svt_ref)^0.217*AR_vt^0.337*(1+taper_vt)^0.363*(cos(sweep_vt_25c/57.3))^-0.484)^1.014;
   
   W_vt_Tor = K_v*Svt_ref*(3.81*(Svt_ref^0.2*V_D/1000/cos(sweep_vt_50c/57.3)^0.5)-0.287);
   % K_v = 1.0 for fuselage mounted horizontal tails
   % K_v = 1+0.15*(Sht_ref*z_h/(Svt_ref*b_vt))   for fin mounted horizontal tails
   
   W_vt = N_vt*(W_vt_USAF + W_vt_GD + W_vt_Tor)/3;
   
   % Vtail folding (only for fixed wing):
   if tiltwing ~= 1,
      W_vt_fold = 0.13 * W_vt/2;  % only one v. tail will be folded
   else
      W_vt_fold = 0;
   end
   
   W_vt_tot = W_vt + W_vt_fold;
   
   I_yy_on_vt = 0;		% Iyy for the component about its own cg
   I_xx_on_vt = 0;		% Iyy for the component about its own cg
   I_zz_on_vt = 0;		% Iyy for the component about its own cg
   x_cg_vt = x_vt + 0.3 * c_r_vt + 0.35*b_vt/2*tan(sweep_vt_le*pi/180);  % 30% chord at 35% semi-span
   y_cg_vt = 0;
   z_cg_vt = h_fus + z_ground2fus - z_cg_fix;
   
   weights(3,1) = W_vt_tot;
   weights(3,2) = x_cg_vt * weights(3,1);
   weights(3,4) = z_cg_vt * weights(3,1);
   % end of vertical tail(s) weights
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Fuselage Weights
   
   W_fus_USAF = 200*((W_to*n_ult/10^5)^0.286*(l_fus/10)^0.857*((w_fus+h_fus)/10)*(V_C/100)^0.338)^1.1;
   
   W_fus_GD = 2*10.43*(K_inl)^1.42*(qbar_D/100)^0.283*(W_to/1000)^0.95*(l_fus/h_fus)^0.71;
   % K_inl = 1.25 for airplanes with inlets in or on the fuselage for a buried engine installation
   % K_inl = 1.0 for inlets located elsewhere
   
   W_fus_Tor = 0.0218*K_fus*(V_D*x_h/(w_fus+h_fus))^0.5*S_fgs^1.2;
   % K_fus = 1.08 for a pressurized fuselage
   % K_fus = 1.07 for a main gear attached to the fuselage
   % K_fus = 1.10 for a cargo airplane with a cargo floor
   % effects are multiplicative for airplanes equipped with all of the above
   
   W_fus_USN = 11.03*K_inl^1.23*(qbar_L/100)^0.245*(W_to/1000)^0.98*(l_fus/h_fus)^0.61;
   
   W_fus = (W_fus_USAF + W_fus_GD + W_fus_Tor + W_fus_USN)/4;
   
   I_yy_on_fus = 0;		% Iyy for the component about its own cg
   I_xx_on_fus = 0;		% Iyy for the component about its own cg
   I_zz_on_fus = 0;		% Iyy for the component about its own cg
   
   x_cg_fus = 35;
   y_cg_fus = 0;
   z_cg_fus = h_fus/2 + z_ground2fus - z_cg_fix;
   
   weights(4,1) = W_fus;
   weights(4,2) = x_cg_fus * W_fus;
   weights(4,4) = z_cg_fus * W_fus;
   % end of fuselage weights
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Nacelle Weights
   
   % USAF method - the nacelle weight is included in the powerplant weight
   
   W_nac_GD = 7.435*N_inl*(A_inl^0.5*l_n*P2)^0.731;
   % for turbofan engines
   
   W_nac_Tor = 0.065*(0.9*T_to_uninst);
   
   W_nac = (W_nac_GD + W_nac_Tor)/2;
   
   W_nac_single = W_nac/N_inl;
   
   [x_cg_eng_in,x_cg_eng_out,x_cg_prop,y_cg_eng_in,y_cg_eng_out,y_cg_prop,z_cg_eng_topwing]...
      =eng_location(tiltwing,w_max,c_r_wing,tc_avg_wing,sweep_wing_25c,taper_wing,AR_wing,b_wing,D_eng,L_eng);
   
   I_yy_on_nac = 0;		% Iyy for the component about its own cg
   I_zz_on_nac = 0;		% Iyy for the component about its own cg
   I_xx_on_nac = 0;		% Iyy for the component about its own cg
   
   x_cg_eng_in = x_cg_eng_in + x_wing;
   x_cg_eng_out = x_cg_eng_out + x_wing;
   x_cg_nac =  (x_cg_eng_in + x_cg_eng_out)/2;
   z_cg_nac = h_fus + z_ground2fus - z_cg_fix + (tc_avg_wing * c_r_wing) + z_cg_eng_topwing; % z_cg_eng ft above wing
   
   weights(5,1) = W_nac;
   weights(5,2) = x_cg_nac*W_nac;
   weights(5,4) = z_cg_nac*W_nac;
   % end of nacelle weights
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Landing Gear Weights
   
   W_gear_USAF = 0.054*(l_s_m)^0.501*(W_L*n_ult_L)^0.684;
   % includes nose gear weight
   % n_ult_L may be taken as 5.7
   
   W_gear_GD = 62.61*(W_to/1000)^0.84;
   
   W_gear_Tor_nose = K_g_r*(A_g_nose + B_g_nose*W_to^0.75 + C_g_nose*W_to + D_g_nose*W_to^1.5);
   W_gear_Tor_main = K_g_r*(A_g_main + B_g_main*W_to^0.75 + C_g_main*W_to + D_g_main*W_to^1.5);
   W_gear_Tor = W_gear_Tor_nose + W_gear_Tor_main;
   % K_g_r = 1.0 for low wing airplanes
   % K_g_r = 1.08 for high wing airplanes
   % the constants A_g through D_g are defined in Table 5.1 of A/P Design Part V by Roskam
   
   W_gear_USN = 129.1*(W_to/1000)^0.66;
   
   
   W_gear = (W_gear_USAF + W_gear_GD + W_gear_Tor + W_gear_USN)/4;
   
   W_gear_nose = (W_gear_Tor_nose/W_gear_Tor)*W_gear;
   W_gear_main = (W_gear_Tor_main/W_gear_Tor)*W_gear;
   
   weights(6,1) = W_gear_nose;
   weights(6,2) = W_gear_nose*x_nosegear;
   weights(6,4) = W_gear_nose*1.5;
   weights(7,1) = W_gear_main;
   weights(7,2) = W_gear_main*x_maingear;
   weights(7,4) = W_gear_main*1.5;
   % end of landing gear weights
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Powerplant Weights
   
   if tiltwing == 1
      W_eng_single = 10^((log10(0.95*P_to_uninst/N_e)+0.1486)/1.1931);
      W_eng = N_e*W_eng_single;
      if (0.95*P_to_uninst/N_p) < 1500
         K_prop_1 = 31.92;
      else
         K_prop_1 = 24.0;
      end
      W_prop = K_prop_1*N_p*N_bl^0.391*((D_p*(0.95*P_to_uninst/N_e))/1000)^0.782;
      W_eng = W_eng + W_prop;
   else   
      W_eng_single = 10^((log10(0.9*(T_to_uninst)/N_e)-0.6072)/1.0229);
      W_eng = N_e*W_eng_single;
   end
   
   W_fs = 80*(N_e + N_t - 1) + 15*N_t^0.5*(W_fuel/K_fsp)^0.333;
   
   W_ec = 88.46*((l_fus+b_wing)*N_e/100)^0.294;
   % for wing mounted engines
   
   W_ess = 38.93*(W_eng/1000)^0.918;
   % for jet engines using electric starting systems
   
   W_osc = 0;   % ?????????????????
   
   W_p_GD = W_ec + W_ess + W_osc;
   
   
   W_p_Tor = 36*N_e*fuelflow_TO + W_osc;
   
   W_p = (W_p_GD + W_p_Tor)/2;
   
   I_xx_on_eng = 0;		% Iyy for the component about its own cg
   I_yy_on_eng = 0;		% Iyy for the component about its own cg
   I_zz_on_eng = 0;		% Iyy for the component about its own cg
   
   [x_cg_eng_in,x_cg_eng_out,x_cg_prop,y_cg_eng_in,y_cg_eng_out,y_cg_prop,z_cg_eng_topwing]...
      =eng_location(tiltwing,w_max,c_r_wing,tc_avg_wing,sweep_wing_25c,taper_wing,AR_wing,b_wing,D_eng,L_eng);
   
   % z_cg_eng_topwing = z-CG position measured from top of wing
   
   x_cg_eng_in = -x_cg_eng_in + x_wing;
   x_cg_eng_out = -x_cg_eng_out + x_wing;
   
   x_cg_eng = (x_cg_eng_in + x_cg_eng_out)/2;
   z_cg_eng = h_fus + z_ground2fus - z_cg_fix + (tc_avg_wing * c_r_wing) + z_cg_eng_topwing; % z_cg_eng_topwing ft above wing
   
   I_yy_on_fs = 0;		% Iyy for the component about its own cg
   x_cg_fs = x_cg_wing + 0.2 * x_cg_eng;
   z_cg_fs = z_cg_wing;
   
   I_xx_on_p = 0;
   I_yy_on_p = 0;			% Iyy for the component about its own cg
   I_zz_on_p = 0;
   x_cg_p = x_cg_eng;
   y_cg_p = 0;
   z_cg_p = z_cg_eng;
   
   weights(8,1) = W_eng;
   weights(8,2) = x_cg_eng * W_eng;
   weights(8,4) = z_cg_eng * W_eng;
   
   weights(9,1) = W_fs;
   weights(9,2) = x_cg_fs * W_fs;
   weights(9,4) = z_cg_fs * W_fs;
   
   weights(10,1) = W_p;
   weights(10,2) = x_cg_p * W_p;
   weights(10,4) = z_cg_p * W_p;
   
   % end of powerplant weights
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Fixed Equipment Weights
   
   
   % Flight controls
   
   W_fc_GD_ct = 56.01*(W_to*qbar_D/100000)^0.576;
   W_fc_Tor_ct = K_fc*W_to^(2/3);
   W_fc_GD_mil = 15.96*(W_to*qbar_L/100000)^0.815;
   % includes the weight of all associated hydraulic and/or pneumatic systems
   
   I_xx_on_fc = 0;
   I_yy_on_fc = 0;		% Iyy for the component about its own cg
   I_zz_on_fc = 0;
   x_cg_fc = 0.4 * (x_wing + y_mac_wing*tan(sweep_wing_le*pi/180));  % 40% chord at MAC
   y_cg_fc = 0;
   z_cg_fc = z_cg_wing;
   
   W_fc = (W_fc_GD_ct + W_fc_Tor_ct + W_fc_GD_mil)/3;
   weights(11,1) = W_fc;
   weights(11,2) = x_cg_fc * W_fc;
   weights(11,4) = z_cg_fc * W_fc;
   
   % Hydraulic system
   W_hydr = K_hydr*W_to;
   
   I_xx_on_hydr = 0;		% Iyy for the component about its own cg
   I_yy_on_hydr = 0;		% Iyy for the component about its own cg
   I_zz_on_hydr = 0;		% Iyy for the component about its own cg
   x_cg_hydr = 0.75 * x_cg_wing + 0.25 * x_cg_ht;			% 75% at wing cg and 25% at tail cg
   y_cg_hydr = 0;
   z_cg_hydr = 0.75 * z_cg_wing + 0.25 * z_cg_ht;
   
   weights(12,1) = W_hydr;
   weights(12,2) = W_hydr * x_cg_hydr;
   weights(12,4) = W_hydr * z_cg_hydr;
   
   
   % Instrumentation, Avionics, and Electronics (iae)
   
   W_iae_GD = N_pil*(15+0.032*(W_to/1000))+N_e*(5+0.006*(W_to/1000))+0.15*(W_to/1000)+0.012*W_to;
   W_iae_Tor = 0.575*W_e^0.556*R^0.25;
   
   I_zz_on_iae = 0;		% Iyy for the component about its own cg
   I_yy_on_iae = 0;		% Iyy for the component about its own cg
   I_xx_on_iae = 0;		% Iyy for the component about its own cg
   x_cg_iae = 0.4 * l_nose;
   y_cg_iae = 0;
   z_cg_iae = 8;
   
   W_iae = (W_iae_GD + W_iae_Tor)/2;
   weights(13,1) = W_iae;
   weights(13,2) = W_iae * x_cg_iae;
   weights(13,4) = W_iae * z_cg_iae;
   
   % Air-conditioning, Pressurization, and Anti/De-icing Systems (api)
   
   W_api_GD = 469*(V_pax*(N_cr+N_pax)/10000)^0.419;
   W_api_Tor = 6.75*(l_pax)^1.28;
   W_api_GD_mil = K_api*(V_pr/100)^0.242;
   % K_api = 887 for subsonic airplanes with wing and tail anti-icing
   
   W_api = (W_api_GD + W_api_Tor + W_api_GD_mil)/3;
   
   I_zz_on_api = 0;		% Iyy for the component about its own cg
   I_yy_on_api = 0;		% Iyy for the component about its own cg
   I_xx_on_api = 0;		% Iyy for the component about its own cg
   x_cg_api = l_nose;
   y_cg_api = 0;
   z_cg_api = 8;
   
   weights(14,1) = W_api;
   weights(14,2) = x_cg_api * W_api;
   weights(14,4) = z_cg_api * W_api;
   
   
   % Electrical System (need weights for fuel system and electronics to do this)
   
   W_els_GD = 1163*((W_fs + W_iae)/1000)^0.506;
   W_els_Tor = 10.8*V_pax^0.7*(1-0.018*V_pax^0.35);
   
   W_els = (W_els_GD + W_els_Tor)/2;
   
   I_yy_on_els = 0;		% Iyy for the component about its own cg
   I_zz_on_els = 0;		% Iyy for the component about its own cg
   I_xx_on_els = 0;		% Iyy for the component about its own cg
   x_cg_els = 0.75 * x_cg_fus + 0.25 * x_cg_p;		% 75% at fuselage CG and 25% at prop CG
   y_cg_els = 0;
   z_cg_els = 0.75 * z_cg_fus + 0.25 * z_cg_p;
   
   weights(15,1) = W_els;
   weights(15,2) = x_cg_els * W_els;
   weights(15,4) = z_cg_els * W_els;
   
   % Oxygen System
   
   W_ox_GD = 7*(N_cr + N_pax)^0.702;
   W_ox_Tor = 30 + 1.2*N_pax;
   
   W_ox = (W_ox_GD + W_ox_Tor)/2;
   
   x_cg_ox = 30;
   z_cg_ox = 8;
   weights(16,1) = W_ox;
   weights(16,2) = W_ox*x_cg_ox;
   weights(16,4) = W_ox*z_cg_ox;
   
   % APU Weight
   W_apu = K_apu*W_to;
   
   I_xx_on_apu = 0;		% Iyy for the component about its own cg
   I_yy_on_apu = 0;		% Iyy for the component about its own cg
   I_zz_on_apu = 0;		% Iyy for the component about its own cg
   x_cg_apu = 35;
   y_cg_apu = 0;
   z_cg_apu = 8;
   
   weights(17,1) = W_apu;
   weights(17,2) = x_cg_apu * W_apu;
   weights(17,4) = z_cg_apu * W_apu;
   
   % Furnishings weight
   
   x_cg_fur = x_cg_fus;
   
   % commercial transport
   W_fur_GD = 55*N_fdc + 32*N_pax + 15*N_cc + K_lav*N_pax^1.33 + K_buf*N_pax^1.12 +...
      109*(N_pax*(1+P_c)/100)^0.505 + 0.771*(W_to/1000);
   W_fur_Tor = 0.211*(W_to - W_fuel)^0.91;
   
   W_fur_comm = (W_fur_GD + W_fur_Tor)/2;
   
   % military transport
   W_fur_GD_mil = 83*N_cr^0.726 + 32*N_pax_mil + 1.11*N_pax_mil^1.33 + 0.771*(W_to/1000);
   W_fur_mil = W_fur_GD_mil;
   
   if commercial == 1
      weights(18,1) = W_fur_comm;
   else
      weights(18,1) = W_fur_mil;
   end
   weights(18,2) = x_cg_fur * weights(18,1);
   weights(18,4) = 6*weights(18,1);
   
   % baggage and cargo handling equipment
   W_bc = K_bc*(N_pax)^1.456;
   
   x_cg_bc = 45;
   weights(19,1) = W_bc;
   weights(19,2) = W_bc*x_cg_bc;
   weights(19,4) = W_bc*6;
   
   % auxiliary items
   W_aux = 0.01*W_e;
   
   x_cg_aux = 30;
   weights(20,1) = W_aux;
   weights(20,2) = W_aux*x_cg_aux;
   weights(20,4) = W_aux*8;
   
   % paint
   W_pt = K_pt*W_to;
   
   x_cg_pt = 30;
   weights(21,1) = W_pt;
   weights(21,2) = W_pt*x_cg_pt;
   weights(21,4) = W_pt*9;
   
   % end of weight calculations
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   
   sumweights = [sum(weights,1)];
   
   W_e_ClassII = sumweights(1);
   x_cg_empty = sumweights(2)/W_e_ClassII;
   y_cg_empty = sumweights(3)/W_e_ClassII;
   z_cg_empty = sumweights(4)/W_e_ClassII;
   
   cg_empty = [x_cg_empty y_cg_empty z_cg_empty]';
   
   %payload = zeros(4,3);
   %payload(1,1) = W_tfo;
   %payload(2,1) = W_crew;
   %payload(3,1) = W_cargo;
   
      
   x_cg_crew = 0.45 * l_nose;
   y_cg_crew = 0;
   z_cg_crew = 0;
   
   x_cg_fuel = x_cg_wing;
   y_cg_fuel = 0;
   z_cg_fuel = z_cg_wing;
   
   x_cg_cargo = l_nose + l_cargo/2;
   y_cg_cargo = 0;
   z_cg_cargo = 0;
   
   x_cg_TO = (sumweights(2) + x_cg_crew*W_crew + x_cg_fuel*W_fuel + x_cg_cargo*W_payload)/(W_e_ClassII+W_crew+W_fuel+W_payload);
   y_cg_TO = (sumweights(3) + y_cg_crew*W_crew + y_cg_fuel*W_fuel + y_cg_cargo*W_payload)/(W_e_ClassII+W_crew+W_fuel+W_payload);
   z_cg_TO = (sumweights(4) + z_cg_crew*W_crew + z_cg_fuel*W_fuel + z_cg_cargo*W_payload)/(W_e_ClassII+W_crew+W_fuel+W_payload);
   cg_TO = [x_cg_TO y_cg_TO z_cg_TO];
   
   x_cg_nocargo = (sumweights(2) + x_cg_crew*W_crew + x_cg_fuel*W_fuel)/(W_e_ClassII+W_crew+W_fuel);
   y_cg_nocargo = (sumweights(3) + y_cg_crew*W_crew + y_cg_fuel*W_fuel)/(W_e_ClassII+W_crew+W_fuel);
   z_cg_nocargo = (sumweights(4) + z_cg_crew*W_crew + z_cg_fuel*W_fuel)/(W_e_ClassII+W_crew+W_fuel);
   cg_nocargo = [x_cg_nocargo y_cg_nocargo z_cg_nocargo];
   
   x_cg_aft = max([x_cg_empty x_cg_TO x_cg_nocargo]);
   x_cg_forward = min([x_cg_empty x_cg_TO x_cg_nocargo]);
   z_cg_high = max([z_cg_empty z_cg_TO z_cg_nocargo]);
   
   W_to_ClassII = W_e_ClassII + W_tfo + W_fuel + W_payload;
   
   W_to_old = W_to;
   W_to = W_to_ClassII;
   W_e = W_e_ClassII;
   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ixx calculation:
%I_xx_off = zeros(numbcomps,1);
%I_xx_off(1) = I_xx_on_wing + (y_cg_wing^2 + z_cg_wing^2) * W_w_tot;
%I_xx_off(2) = I_xx_on_ht + (y_cg_ht^2 + z_cg_ht^2) * W_ht;
%I_xx_off(3) = I_xx_on_vt + ((b_ht/2)^2 + z_cg_vt^2) * W_vt_tot;
%I_xx_off(4) = I_xx_on_fus + (y_cg_fus^2 + z_cg_fus^2) * W_fus;
%I_xx_off(5) = I_xx_on_nac + 2*(y_cg_eng_in^2 + z_cg_eng^2 + y_cg_eng_out^2 + z_cg_eng^2) * W_nac_single;
%I_xx_off(8) = I_xx_on_eng + 2*(y_cg_eng_in^2 + z_cg_eng^2 + y_cg_eng_out^2 + z_cg_eng^2) * W_eng_single;
%I_xx_off(9) = I_xx_on_fus + (y_cg_fus^2 + z_cg_fs^2) * W_fs;
%I_xx_off(10) = I_xx_on_p + (y_cg_p^2 + z_cg_p^2) * W_p;
%I_xx_off(11) = I_xx_on_fc + (y_cg_fc^2 + z_cg_fc^2) * W_fc;
%I_xx_off(12) = I_xx_on_hydr + (y_cg_hydr^2 + z_cg_hydr^2) * W_hydr;
%I_xx_off(13) = I_xx_on_iae + (y_cg_iae^2 + z_cg_iae^2) * W_iae;
%I_xx_off(14) = I_xx_on_api + (y_cg_api^2 + z_cg_api^2) * W_api;
%I_xx_off(15) = I_xx_on_els + (y_cg_els^2 + z_cg_els^2) * W_els;
%I_xx_off(17) = I_xx_on_apu + (y_cg_apu^2 + z_cg_apu^2) * W_apu;

%I_xx = sum(I_xx_off)/32.2;		% total aircraft Ixx

% Iyy calculation:

I_yy_off = zeros(numbcomps + 1,1);
I_yy_off(1) = I_yy_on_wing + (x_cg_wing^2 + z_cg_wing^2) * W_w_tot;
I_yy_off(2) = I_yy_on_ht + (x_cg_ht^2 + z_cg_ht^2) * W_ht;
I_yy_off(3) = I_yy_on_vt + (x_cg_vt^2 + z_cg_vt^2) * W_vt_tot;
I_yy_off(4) = I_yy_on_fus + (x_cg_fus^2 + z_cg_fus^2) * W_fus;
I_yy_off(5) = I_yy_on_nac + 2*(x_cg_eng_in^2 + z_cg_eng^2 + x_cg_eng_out^2 + z_cg_eng^2) * W_nac_single;
I_yy_off(8) = I_yy_on_eng + 2*(x_cg_eng_in^2 + z_cg_eng^2 + x_cg_eng_out^2 + z_cg_eng^2) * W_eng_single;
I_yy_off(9) = I_yy_on_fs + (x_cg_fs^2 + z_cg_fs^2) * W_fs;
I_yy_off(10) = I_yy_on_p + (x_cg_p^2 + z_cg_p^2) * W_p;
I_yy_off(11) = I_yy_on_fc + (x_cg_fc^2 + z_cg_fc^2) * W_fc;
I_yy_off(12) = I_yy_on_hydr + (x_cg_hydr^2 + z_cg_hydr^2) * W_hydr;
I_yy_off(13) = I_yy_on_iae + (x_cg_iae^2 + z_cg_iae^2) * W_iae;
I_yy_off(14) = I_yy_on_api + (x_cg_api^2 + z_cg_api^2) * W_api;
I_yy_off(15) = I_yy_on_els + (x_cg_els^2 + z_cg_els^2) * W_els;
I_yy_off(17) = I_yy_on_apu + (x_cg_apu^2 + z_cg_apu^2) * W_apu;

I_yy_off(22) = (x_cg_cargo^2 + z_cg_cargo^2) * W_payload;
I_yy = sum(I_yy_off)/32.2;		% total aircraft Iyy

% Izz calculation:
%I_zz_off = zeros(numbcomps + 1,1);
%I_zz_off(1) = I_zz_on_wing + (y_cg_wing^2 + x_cg_wing^2) * W_w_tot;
%I_zz_off(2) = I_zz_on_ht + (y_cg_ht^2 + x_cg_ht^2) * W_ht;
%I_zz_off(3) = I_zz_on_vt + ((b_ht/2)^2 + x_cg_vt^2) * W_vt_tot;
%I_zz_off(4) = I_zz_on_fus + (y_cg_fus^2 + x_cg_fus^2) * W_fus;
%I_zz_off(5) = I_zz_on_nac + 2*(x_cg_eng_in^2 + y_cg_eng_in^2 + x_cg_eng_out^2 + y_cg_eng_out^2) * W_nac_single;
%I_zz_off(8) = I_zz_on_eng + 2*(x_cg_eng_in^2 + y_cg_eng_in^2 + x_cg_eng_out^2 + y_cg_eng_out^2) * W_eng_single;
%I_zz_off(9) = I_zz_on_fus + (y_cg_fus^2 + x_cg_fus^2) * W_fus;
%I_zz_off(10) = I_zz_on_p + (y_cg_p^2 + x_cg_p^2) * W_p;
%I_zz_off(11) = I_zz_on_fc + (y_cg_fc^2 + x_cg_fc^2) * W_fc;
%I_zz_off(12) = I_zz_on_hydr + (y_cg_hydr^2 + x_cg_hydr^2) * W_hydr;
%I_zz_off(13) = I_zz_on_iae + (y_cg_iae^2 + x_cg_iae^2) * W_iae;
%I_zz_off(14) = I_zz_on_api + (y_cg_api^2 + x_cg_api^2) * W_api;
%I_zz_off(15) = I_zz_on_els + (y_cg_els^2 + x_cg_els^2) * W_els;
%I_zz_off(17) = I_zz_on_apu + (y_cg_apu^2 + x_cg_apu^2) * W_apu;

%I_zz_off(22) = (x_cg_cargo^2 + y_cg_cargo^2) * W_payload;
%I_zz = sum(I_zz_off)/32.2;		% total aircraft Izz


return

% Outputs
fprintf('\n');
fprintf('Weights (lbs)\n');
fprintf('-----------------------------------------\n');
fprintf('Wing =                  %6.0f\n',weights(1,1));
fprintf('H. Tail =               %6.0f\n',weights(2,1));
fprintf('V. Tail(s) =            %6.0f\n',weights(3,1));
fprintf('Fuselage =              %6.0f\n',weights(4,1));
fprintf('Nacelles =              %6.0f\n',weights(5,1));
fprintf('Nose Gear =             %6.0f\n',weights(6,1));
fprintf('Main Gear =             %6.0f\n',weights(7,1));
fprintf('Engines =               %6.0f\n',weights(8,1));
fprintf('Fuel System =           %6.0f\n',weights(9,1));
fprintf('Propulsion =            %6.0f\n',weights(10,1));
fprintf('Flight Controls =       %6.0f\n',weights(11,1));
fprintf('Hydraulic =             %6.0f\n',weights(12,1));
fprintf('Instr, Avionics, Elec = %6.0f\n',weights(13,1));
fprintf('Air Cond, Press =       %6.0f\n',weights(14,1));
fprintf('Electrical =            %6.0f\n',weights(15,1));
fprintf('Oxygen =                %6.0f\n',weights(16,1));
fprintf('APU =                   %6.0f\n',weights(17,1));
fprintf('Furnishings =           %6.0f\n',weights(18,1));
fprintf('Bag Handling =          %6.0f\n',weights(10,1));
fprintf('Auxiliary =             %6.0f\n',weights(20,1));
fprintf('Paint =                 %6.0f\n',weights(21,1));
fprintf('\n');

fprintf('Empty =                 %6.0f\n',W_e_ClassII);
fprintf('\n');

return

fprintf('Trap. Fuel & Oil =      %6.0f\n',payload(1,1));
fprintf('Crew =                  %6.0f\n',payload(2,1));
fprintf('Fuel =                  %6.0f\n',payload(3,1));
fprintf('Cargo =                 %6.0f\n',payload(4,1));
fprintf('\n');
fprintf('Takeoff =               %6.0f\n',W_to);
fprintf('\n');

