%%%%%%%%%%%%%%%%%%%%
%%%%Appendix F-3
%%%%Team Orion
%%%%Structure Analysis
%%%%Last Updated 10/25/00
%%%%%%%%%%%%%%%%%%%%

clear all
close all

%Input parameters

b=6.63;  		 	    								   % span(ft)
b_half = b/2;        								   % half span  
b_htail = 3.1634;		 								   % (ft)              
b_half_htail = b_htail/2;								% (ft)
c=b/4.4; 		 	 	  								   % chord(ft)
tc=.14;  		 	 	  									% thickness to chord ratio
c_root_htail = 1.3181;								   % (ft)
c_tip_htail = .7908;   									% (ft)
c_ave_htail = (c_root_htail + c_tip_htail)/2;   % (ft)
w=10.2 ;			 	     			  						% weight of plane (lb) 
lambda=1;	 	 	     			  						% taper ratio
rho_balsa = 11 *.031080950172 ;					   % density of balsa(slug/ft^3)
rho_spruce = 34 *.031080950172 ;						% density of spruce(slug/ft^3)
E_balsa=625; 	 										   % Youngs Modulus for balsa(ksi)
E_spruce=1500;	 										   % Youngs Modulus for balsa(ksi)
v=1 ;			  	 	 										% Poissons ratio
th_ribs=1/16/12;   										% rib thickness(ft)
cg_plane = 1.071;										   % (ft)
x_sid5 = 6.38;									         % (ft) 
n_ribs_wing = 29;
n_ribs_tail = 14;

%wing dimensions

x_sparcap=-1/8/12; 									   % (ft)
z_sparcap=-3/8/12;										% (ft)
z_shearweb=-c*tc;										   % (ft)
x_shearweb=-1/16/12;  									% (ft)
A_sparcap=x_sparcap * z_sparcap;						% Area of spar caps(ft^2)
V_cap = A_sparcap * b;									% (ft^3) 
A_shearweb=x_shearweb * z_shearweb;					% Area of shear web(ft^2)
V_shearweb = A_shearweb * b;							% (ft^3)
z_betweenwings = -.15*b;     							% (ft)
A_spar = A_shearweb + A_sparcap;					   % (ft^2)
V_spar = A_spar * A_sparcap * b;					   % (ft^3)
V_rib = (7.4126e-004); 									% (ft^3)
A_wing = .20097;											% (ft^2)
V_wing = A_wing * b;									   % (ft^3)
V_wings = V_wing * 2;

%horizontal tail dimensions
x_sparcap_htail = -1/8/12;  							            % (ft)
z_sparcap_htail = -1/8/12;  							            % (ft)
x_shearweb_htail = -1/16/12;    					    			   % (ft)   
z_shearweb_htail = -0.1186; 							 			   % (ft)
A_sparcap_htail = x_sparcap_htail * z_sparcap_htail;		   % (ft^2)
V_cap_htail = A_sparcap_htail * b_htail; 						   % (ft^3)
A_shearweb_htail = x_shearweb_htail * z_shearweb_htail; 	   % (ft^2) 
V_web_htail = A_shearweb_htail * b_htail;               	   % (ft^3)
V_htailspar = A_sparcap_htail * A_shearweb_htail * b_htail; % (ft^3) 
A_htail = c_ave_htail * .04 * c_ave_htail;					   % (ft)
V_htail = (.5 * (c_root_htail - c_tip_htail) * b_htail * .04 * c_ave_htail) + (c_tip_htail * b_htail * .04 * c_ave_htail);% (ft^3)
V_rib_htail = (c_ave_htail * (c_ave_htail * .04) * th_ribs);% (ft^3)

%vertical tail dimensions
z_vtail = -1.24;      												 % (ft)
x_vtail_root = -1.2522;												 % (ft)
x_vtail_tip  = -0.7513;												 % (ft) 
x_vtail_ave = (x_vtail_root + x_vtail_tip)/2;				 % (ft)
b_vtail = x_vtail_root - x_vtail_tip;							 % (ft) 
A_vtail_tri = .5 * b_vtail * z_vtail;     					 % (ft^2)
A_vtail_rec = z_vtail * x_vtail_tip; 			         	 % (ft^2)
A_vtail = A_vtail_tri + A_vtail_rec;			  				 % (ft^2)
V_vtail = A_vtail_tri * A_vtail_rec * (.04 * x_vtail_ave);% (ft^2)

%fuselage dimensions
z_fuselage =-.5;     	       % (ft)
y_fuselage = .417;   	       % (ft)
y_fuselage_half = y_fuselage/2;% (ft)
x_fuselage_front = -.91;       % (ft)
z_fuselage1a = - .0883;        % (ft)
z_fuselage1b = - .25;	       % (ft)
z_fuselage1c = - .292;   		 % (ft)
A_fuselage1a = .5 * z_fuselage1a * x_fuselage_front;				  % (ft^2)
A_fuselage1b = x_fuselage_front* z_fuselage1b;	          		  % (ft^2)
A_fuselage1c = .5 * x_fuselage_front * z_fuselage1c;				  % (ft^2)
V_fuselage1a = .5 * z_fuselage1a * x_fuselage_front * y_fuselage;% (ft^3)
V_fuselage1b = x_fuselage_front * z_fuselage1b * y_fuselage;	  % (ft^3)
V_fuselage1c = .5 * x_fuselage_front * z_fuselage1c * y_fuselage;% (ft^3)
x_fuselage2 = -1.5;  	 													  %( ft) length of front part of fuselage 
x_fuselage3 = -4.35; 	 													  %( ft) length of back part of fuselage
x_fuselagetotal = x_fuselage_front + x_fuselage2 + x_fuselage3;  %( ft)
A_fuselage2 = z_fuselage * x_fuselage2;    						     % (ft^2)
V_fuselage2 = z_fuselage * x_fuselage2 * y_fuselage;     		  % (ft^3)
A_fuselage3 = .5*(z_fuselage * x_fuselage3);					 		  % (ft^3)
V_fuselage3 = .5*(z_fuselage * x_fuselage3 * y_fuselage);		  % (ft^3)
A_fuselage = A_fuselage1a + A_fuselage1b + A_fuselage1c + A_fuselage2 + A_fuselage3;% (ft^2)
V_fuselage = V_fuselage1a + V_fuselage1b + V_fuselage1c + V_fuselage2 + V_fuselage3;% (ft^3) 

A_sid5 = A_fuselage + A_vtail_tri + A_vtail_rec + A_sparcap_htail + A_shearweb_htail + A_spar; % (ft^2) from side view
V_sid5 = V_fuselage + V_vtail + V_htail + V_wings% (ft^3)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%control system weight estimates
n_servo=4; 						   %number of servos required
weight_servo=2.0/16 * n_servo;%weight of a servo (lbs)
weight_receiver=1.0/16; 	   %weight of receiver in (lbs)
weight_speedcontrol=3.0/16; 	%weight of electronic spped control in (lbs)
weight_gyro=3.5/16; 			   %(lbs)
weight_datarecord=8/16; 	   %(lbs)
weight_interface = 7/16;      % (lbs)


%motor weight estimates
weight_motor=7.5/16;   %estimate taken from MaxCim data (lbs)
weight_gearbox=1.5/16; %(lbs)
weight_propeller=2/16; %(lbs)

%Battery Pack
cells=18; 		    						 %number of cells-will change with power needed
cell_weight=2.8/16;                  %lbs-const once battery is chosen
weight_battery = cell_weight * cells;% (lbs)
%Calculate weight based on area
      %weight_wing=(2.75*S-1.9564)/16;			        %(lbs) Orion Database
      %weight_tail=(2.75*(SHT+SVT)-1.9564)/16; 		  %(lbs) Orion Database
      %thickness_fuselage=(1/4)/12;      %(in)->ft
      %balsa_density=10;  					  %(lbs/ft^3)
      %spruce_density=28; 					  %(lbs/ft^3)
      %weight_fuselage=((Swet_fuse*thickness_fuselage*balsa_density)*1.15); %(lbs)
      %weight_struts=(Swet_struts/2*0.00833*spruce_density); 					%(lbs)
      %weight_mainwheel=1.28/16; %(lbs) 2.5" tire
      %weight_tailwheel=.5/16;   %(lbs)
      %weight_landinggear=5/16;  %(lbs)
      %Wnew=weight_wing+weight_tail+weight_fuselage+2*weight_mainwheel+weight_tailwheel+weight_landinggear+weight_payload; %(lbs)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Centroid of Spar
xc_cap = x_sparcap/2; % (ft)
zc_cap = z_sparcap/2; % (ft)
xc_web = x_shearweb/2;% (ft)
zc_web = z_shearweb/2;% (ft)

%zc_leadingedge = width_leadingedge/2;   % (ft)
%yc_leadingedge = height_leadingedge/2;  % (ft)
%zc_trailingedge = width_trailingedge/2; % (ft)
%yc_trailingedge = height_trailingedge/2;% (ft)

xc_spar = ((2*A_sparcap * xc_cap) + (A_shearweb * xc_web))/(2*A_sparcap + A_shearweb); 				     % (ft)
zc_spar = ((A_sparcap * zc_cap) + (A_sparcap * zc_cap) + (A_shearweb * zc_web))/(2*A_sparcap + A_shearweb);% (ft)

xc_rib = -7.4835/12;% (ft)
yc_rib = .0312/12;  % (ft)
zc_rib = -.5692/12; % (ft)

%Centroid of fuselage
xc_fuselage1a = x_fuselage_front * (2/3);% (ft)
zc_fuselage1a = z_fuselage1a/3; 	 	 	  % (ft)
xc_fuselage1b = x_fuselage_front/2;   	  % (ft)
zc_fuselage1b = z_fuselage1a/2;			  % (ft)
xc_fuselage1c = x_fuselage_front * (2/3);% (ft)
zc_fuselage1c = z_fuselage1c * (2/3);    % (ft)

xc_fuselage2 = x_fuselage2/2;            % (ft)
zc_fuselage2 = z_fuselage/2;       	     % (ft)

xc_fuselage3 = x_fuselage3/3;     	     % (ft)
zc_fuselage3 = z_fuselage * (2/3);	     % (ft)

xc_fuselage = (A_fuselage1a * xc_fuselage1a) + (A_fuselage1b * xc_fuselage1b)+ (A_fuselage1c * xc_fuselage1c)+(A_fuselage2 *...
   xc_fuselage2) + (A_fuselage3 *xc_fuselage3)/(A_fuselage2 + A_fuselage3 + A_fuselage1a + A_fuselage1b + A_fuselage1c); % (ft)
zc_fuselage = (A_fuselage1a * zc_fuselage1a) + (A_fuselage1b * zc_fuselage1b)+ (A_fuselage1c * zc_fuselage1c)+(A_fuselage2 * ...
   zc_fuselage3) + (A_fuselage3 *zc_fuselage3)/(A_fuselage2 + A_fuselage3 + A_fuselage1a + A_fuselage1b + A_fuselage1c); % (ft)

%Centroid of H-tail
xc_cap_htail = x_sparcap_htail/2; % (ft)
zc_cap_htail = z_sparcap_htail/2; % (ft)
xc_web_htail = x_shearweb_htail/2;% (ft)
zc_web_htail = z_shearweb_htail/2;% (ft)

xc_spar_htail = ((2*A_sparcap_htail * xc_cap_htail) + (A_shearweb_htail * xc_web_htail))/...
   (2*A_sparcap_htail + A_shearweb_htail);% (ft)
zc_spar_htail = ((A_sparcap_htail * zc_cap_htail) + (A_sparcap_htail * zc_cap_htail) + ...
   (A_shearweb_htail * zc_web_htail))/(2*A_sparcap_htail + A_shearweb_htail);% (ft)


%Centroid of V-tail
xc_vtail_tri = b_vtail * (2/3);% (ft)
zc_vtail_tri = z_vtail/3; % (ft) 

xc_vtail_rec = x_vtail_tip/2;% (ft)
zc_vtail_rec = z_vtail/2;    % (ft)

xc_vtail = (A_vtail_tri * xc_vtail_tri) + (A_vtail_rec * xc_vtail_rec)/(A_vtail_rec + A_vtail_tri);% (ft)
zc_vtail = (A_vtail_tri * zc_vtail_tri) + (A_vtail_rec * zc_vtail_rec)/(A_vtail_rec + A_vtail_tri);% (ft)

%Centroid of internal equipment 



A_gearbox_side = (1*3)/144;           % (ft^2)
A_gearbox_top = (1*1.5)/144;			  % (ft^2)
A_gearbox_front = (1.5*3)/144;        % (ft^2)
xc_gearbox = -.5/12;  		  		     % (ft)
zc_gearbox = -5.5/12;		           % (ft)

A_motorside_side = (1.5*2.25)/144;    % (ft^2)
A_motorside_top = (1.5 * 2.25)/144;   % (ft^2)
A_motorside_front = (1.5*1.5)/144;    % (ft^2)
xc_motor = -2.125/12; 		 		     % (ft)
zc_motor = -5.75/12;        		     % (ft)

A_speedcontrol_side = (1.5*1)/144;    % (ft^2)
A_speedcontrol_top = (1.25*1.5);      % (ft^2) 
A_speedcontrol_front = (1.5*1)/144;   % (ft^2)
xc_speedcontrol = -6.875/12;		     % (ft)
zc_speedcontrol = -5.5/12;	 		     % (ft)

A_receiver_side = (.75*1.75)/144;     % (ft^2)
A_receiver_top = (.75*1.25)/144;      % (ft^2)
A_receiver_front = (1.25*1.75)/144;   % (ft^2)
xc_receiver = -8.375/12;	 		     % (ft)
zc_receiver = -6.5/12;		 		     % (ft)

A_gyro_side = (1.25*1.5)/144; 	     % (ft^2) 
A_gyro_top = (1.5*1.25)/144;
A_gyro_front = (1.5*1.25)/144; 
xc_gyro = -6.875/12;			 		     % (ft)
zc_gyro = -6.25/12;			 		  	  % (ft)

A_datarecord_side = (2.25*1.75)/144;  % (ft^2)
A_datarecord_top = (2.25*3.25)/144;   % (ft^2)
A_datarecord_front = (1.75*3.25)/144; % (ft^2)  
xc_datarecord = - 6.875/12;           % (ft)
zc_datarecord = - 4.625/12;           % (ft)

A_battery_side = (2*5)/144;			  % (ft^2)
A_battery_top = (2*4)/144;            % (ft^2)
A_battery_front = (5*2)/144; 
xc_battery = -4.75/12;		 			  % (ft)
zc_battery = -3.75/12;		 			  % (ft)

A_servo_side = (1.25*1.5)/144;		  % (ft^2)
A_servo_top = (.75*3)/144; 			  % (ft^2)
A_servo_front = (1.25*3)/144; 
xc_servo = -24/12;					  	  % (ft)
zc_servo = -2.5/12;			 			  % (ft)

A_interface_side = (1.25*5.75)/144;	  % (ft^2)
A_interface_top = (5.75*3.25)/144;	  % (ft^2)
A_interface_front = (1.25*3.5)/144;   % (ft^2)
xc_interface = -9.675/12; 	 			  % (ft)
zc_interface = -2.125/12;				  % (ft)




%Centroid of sid5
xc_sid5 = ((A_vtail * xc_vtail) + ( A_htail * xc_spar_htail) + (A_fuselage * xc_fuselage) + (A_wing * xc_spar) + (A_interface_side * xc_interface)...
   + (A_gearbox_side * xc_gearbox) + (A_servo_side * xc_servo) + (A_battery_side * xc_battery) + (A_datarecord_side * xc_datarecord) + ...
   (A_gyro_side * xc_gyro) + (A_receiver_side * xc_receiver)+ (A_speedcontrol_side * xc_speedcontrol) + (A_motorside_side * xc_motor))/...
   ((A_vtail+ A_htail + A_fuselage + A_wing + A_interface_side + A_gearbox_side + A_servo_side + A_battery_side + A_datarecord_side + A_gyro_side +...
   A_receiver_side + A_speedcontrol_side + A_motorside_side))% (ft) 



zc_sid5 = ((A_vtail * zc_vtail) + ( A_htail * zc_spar_htail) + (A_fuselage * zc_fuselage) + (A_wing *zc_spar) + (A_interface_side * zc_interface)...
   + (A_gearbox_side * zc_gearbox) + (A_servo_side * zc_servo) + (A_battery_side * zc_battery) + (A_datarecord_side * zc_datarecord) + ...
   (A_gyro_side * zc_gyro) + (A_receiver_side * zc_receiver)+ (A_speedcontrol_side * zc_speedcontrol) + (A_motorside_side * zc_motor))/...
   ((A_vtail+ A_htail + A_fuselage + A_wing + A_interface_side + A_gearbox_side + A_servo_side + A_battery_side + A_datarecord_side + A_gyro_side +...
   A_receiver_side + A_speedcontrol_side + A_motorside_side))% (ft) 

cgx_sid5 = 1.071;  % (ft)
cgz_sid5 = zc_sid5;% (ft)


%Distances from center of gravity

x_spar_top = abs((.91+1.5*.25) - cgx_sid5);				 	 % (ft)
z_spar_top = abs((.15 * b) + cgz_sid5);				     	 % (ft)
x_spar_bot = abs((.91+1.5*.25) - cgx_sid5);				 	 % (ft)
z_spar_bot = abs(-cgz_sid5 + zc_spar);     					 % (ft)
x_rib_top = abs(1.5336 - cgx_sid5);		   				 	 % (ft) xc of rib + .91 - xc_sid5
z_rib_top = abs((.15 * b) + cgz_sid5);  				       % (ft)
x_rib_bot = abs(1.5336 - cgx_sid5);		   			    	 % (ft)
z_rib_bot = abs(cgz_sid5 - zc_rib);       				 	 % (ft) 
x_vtail = abs((x_sid5 - zc_vtail)- cgz_sid5);      	    % (ft)
z_vtail = abs((z_fuselage + zc_vtail) - cgz_sid5); 	    % (ft)
x_htail = abs((x_sid5 - xc_spar_htail) - cgx_sid5);	    % (ft)
z_htail = abs(z_fuselage - cgz_sid5);

y_rib1 = 3/12;
y_rib2 = 6/12;
y_rib3 = 9/12;
y_rib4 = 12/12;
y_rib5 = 15/12;
y_rib6 = 18/12;
y_rib7 = 21/12;
y_rib8 = 24/12;
y_rib9 = 27/12;
y_rib10 = 30/12;
y_rib11 = 33/12;
y_rib12 = 36/12;
y_rib13 = 39/12;
y_rib14 = -3/12;
y_rib15 = -6/12;
y_rib16 = -9/12;
y_rib17 = -12/12;
y_rib18 = -15/12;
y_rib19 = -18/12;
y_rib20 = -21/12;
y_rib21 = -24/12;
y_rib22 = -27/12;
y_rib23 = -30/12;
y_rib24 = -33/12;
y_rib25 = -36/12;
y_rib26 = -39/12;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Moment of inertia for wing
Ixc_cap = (x_sparcap * z_sparcap^3)/12;		   % (ft^4)
Izc_cap = (z_sparcap * x_sparcap^3)/12;		   % (ft^4)
Ixc_web = (x_shearweb * z_shearweb^3)/12;       % (ft^4)
Izc_web = (z_shearweb * x_shearweb^3)/12;       % (ft^4)

%Mass moment of inertia for wing
Izm_shearweb_top = ((x_spar_top)^2)* rho_balsa * V_shearweb; 						 % (slug ft^2)
Ixm_shearweb_top = ((z_spar_top)^2)* rho_balsa * V_shearweb; 						 % (slug ft^2)
Iym_shearweb_top = ((x_spar_top)^2 + (z_spar_top)^2 )* rho_balsa * V_shearweb; % (slug ft^2)
Izm_cap_top = ((x_spar_top)^2)* rho_spruce * V_cap;		    					    % (slug ft^2)
Ixm_cap_top = ((z_spar_top)^2)* rho_spruce * V_cap;		    					    % (slug ft^2)
Iym_cap_top = ((x_spar_top)^2 + (z_spar_top)^2 )* rho_spruce * V_cap;		    % (slug ft^2)

Izm_shearweb_bot = ((x_spar_bot)^2)* rho_balsa * V_shearweb; 					    % (slug ft^2)
Ixm_shearweb_bot = ((z_spar_bot)^2)* rho_balsa * V_shearweb; 					    % (slug ft^2)
Iym_shearweb_bot = ((x_spar_bot)^2 + (z_spar_bot)^2)* rho_balsa * V_shearweb;  % (slug ft^2)
Izm_cap_bot = ((x_spar_bot)^2)* rho_spruce * V_cap;		    					    % (slug ft^2)
Ixm_cap_bot = ((z_spar_bot)^2)* rho_spruce * V_cap;		    					    % (slug ft^2)
Iym_cap_bot = ((x_spar_bot)^2 + (z_spar_bot)^2)* rho_spruce * V_cap;			    % (slug ft^2)

Izm_spar_top = Izm_shearweb_top + Izm_cap_top;					 % (slug ft^2)
Ixm_spar_top = Ixm_shearweb_top + Ixm_cap_top;					 % (slug ft^2)
Iym_spar_top = Iym_shearweb_top + Iym_cap_top;					 % (slug ft^2)
Izm_spar_bot = Izm_shearweb_bot + Izm_cap_bot;					 % (slug ft^2)
Ixm_spar_bot = Ixm_shearweb_bot + Ixm_cap_bot;					 % (slug ft^2)
Iym_spar_bot = Iym_shearweb_bot + Iym_cap_bot;					 % (slug ft^2)

Izm_rib_top1  = ((x_rib_top^2)+ (y_rib1^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top2  = ((x_rib_top^2)+ (y_rib2^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top3  = ((x_rib_top^2)+ (y_rib3^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top4  = ((x_rib_top^2)+ (y_rib4^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top5  = ((x_rib_top^2)+ (y_rib5^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top6  = ((x_rib_top^2)+ (y_rib6^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top7  = ((x_rib_top^2)+ (y_rib7^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top8  = ((x_rib_top^2)+ (y_rib8^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top9  = ((x_rib_top^2)+ (y_rib9^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top10 = ((x_rib_top^2)+ (y_rib10^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top11 = ((x_rib_top^2)+ (y_rib11^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top12 = ((x_rib_top^2)+ (y_rib12^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top13 = ((x_rib_top^2)+ (y_rib13^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top14 = ((x_rib_top^2)+ (y_rib14^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top15 = ((x_rib_top^2)+ (y_rib15^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top16 = ((x_rib_top^2)+ (y_rib16^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top17 = ((x_rib_top^2)+ (y_rib17^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top18 = ((x_rib_top^2)+ (y_rib18^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top19 = ((x_rib_top^2)+ (y_rib19^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top20 = ((x_rib_top^2)+ (y_rib20^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top21 = ((x_rib_top^2)+ (y_rib21^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top22 = ((x_rib_top^2)+ (y_rib22^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top23 = ((x_rib_top^2)+ (y_rib23^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top24 = ((x_rib_top^2)+ (y_rib24^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top25 = ((x_rib_top^2)+ (y_rib25^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_top26 = ((x_rib_top^2)+ (y_rib26^2)) * V_rib * rho_balsa; 				 % (slug ft^2)


Izm_rib_bot1  = ((x_rib_bot^2)+ (y_rib1^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot2  = ((x_rib_bot^2)+ (y_rib2^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot3  = ((x_rib_bot^2)+ (y_rib3^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot4  = ((x_rib_bot^2)+ (y_rib4^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot5  = ((x_rib_bot^2)+ (y_rib5^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot6  = ((x_rib_bot^2)+ (y_rib6^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot7  = ((x_rib_bot^2)+ (y_rib7^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot8  = ((x_rib_bot^2)+ (y_rib8^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot9  = ((x_rib_bot^2)+ (y_rib9^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot10 = ((x_rib_bot^2)+ (y_rib10^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot11 = ((x_rib_bot^2)+ (y_rib11^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot12 = ((x_rib_bot^2)+ (y_rib12^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot13 = ((x_rib_bot^2)+ (y_rib13^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot14 = ((x_rib_bot^2)+ (y_rib14^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot15 = ((x_rib_bot^2)+ (y_rib15^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot16 = ((x_rib_bot^2)+ (y_rib16^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot17 = ((x_rib_bot^2)+ (y_rib17^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot18 = ((x_rib_bot^2)+ (y_rib18^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot19 = ((x_rib_bot^2)+ (y_rib19^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot20 = ((x_rib_bot^2)+ (y_rib20^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot21 = ((x_rib_bot^2)+ (y_rib21^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot22 = ((x_rib_bot^2)+ (y_rib22^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot23 = ((x_rib_bot^2)+ (y_rib23^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot24 = ((x_rib_bot^2)+ (y_rib24^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot25 = ((x_rib_bot^2)+ (y_rib25^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Izm_rib_bot26 = ((x_rib_bot^2)+ (y_rib26^2)) * V_rib * rho_balsa; 				 % (slug ft^2)

Ixm_rib_top1  = ((z_rib_top^2)+ (y_rib1^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top2  = ((z_rib_top^2)+ (y_rib2^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top3  = ((z_rib_top^2)+ (y_rib3^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top4  = ((z_rib_top^2)+ (y_rib4^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top5  = ((z_rib_top^2)+ (y_rib5^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top6  = ((z_rib_top^2)+ (y_rib6^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top7  = ((z_rib_top^2)+ (y_rib7^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top8  = ((z_rib_top^2)+ (y_rib8^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top9  = ((z_rib_top^2)+ (y_rib9^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top10 = ((z_rib_top^2)+ (y_rib10^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top11 = ((z_rib_top^2)+ (y_rib11^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top12 = ((z_rib_top^2)+ (y_rib12^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top13 = ((z_rib_top^2)+ (y_rib13^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top14 = ((z_rib_top^2)+ (y_rib14^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top15 = ((z_rib_top^2)+ (y_rib15^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top16 = ((z_rib_top^2)+ (y_rib16^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top17 = ((z_rib_top^2)+ (y_rib17^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top18 = ((z_rib_top^2)+ (y_rib18^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top19 = ((z_rib_top^2)+ (y_rib19^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top20 = ((z_rib_top^2)+ (y_rib20^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top21 = ((z_rib_top^2)+ (y_rib21^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top22 = ((z_rib_top^2)+ (y_rib22^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top23 = ((z_rib_top^2)+ (y_rib23^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top24 = ((z_rib_top^2)+ (y_rib24^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top25 = ((z_rib_top^2)+ (y_rib25^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_top26 = ((z_rib_top^2)+ (y_rib26^2)) * V_rib * rho_balsa; 				 % (slug ft^2)


Ixm_rib_bot1  = ((z_rib_bot^2)+ (y_rib1^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot2  = ((z_rib_bot^2)+ (y_rib2^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot3  = ((z_rib_bot^2)+ (y_rib3^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot4  = ((z_rib_bot^2)+ (y_rib4^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot5  = ((z_rib_bot^2)+ (y_rib5^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot6  = ((z_rib_bot^2)+ (y_rib6^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot7  = ((z_rib_bot^2)+ (y_rib7^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot8  = ((z_rib_bot^2)+ (y_rib8^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot9  = ((z_rib_bot^2)+ (y_rib9^2))  * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot10 = ((z_rib_bot^2)+ (y_rib10^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot11 = ((z_rib_bot^2)+ (y_rib11^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot12 = ((z_rib_bot^2)+ (y_rib12^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot13 = ((z_rib_bot^2)+ (y_rib13^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot14 = ((z_rib_bot^2)+ (y_rib14^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot15 = ((z_rib_bot^2)+ (y_rib15^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot16 = ((z_rib_bot^2)+ (y_rib16^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot17 = ((z_rib_bot^2)+ (y_rib17^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot18 = ((z_rib_bot^2)+ (y_rib18^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot19 = ((z_rib_bot^2)+ (y_rib19^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot20 = ((z_rib_bot^2)+ (y_rib20^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot21 = ((z_rib_bot^2)+ (y_rib21^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot22 = ((z_rib_bot^2)+ (y_rib22^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot23 = ((z_rib_bot^2)+ (y_rib23^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot24 = ((z_rib_bot^2)+ (y_rib24^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot25 = ((z_rib_bot^2)+ (y_rib25^2)) * V_rib * rho_balsa; 				 % (slug ft^2)
Ixm_rib_bot26 = ((z_rib_bot^2)+ (y_rib26^2)) * V_rib * rho_balsa; 				 % (slug ft^2)





Izm_rib_top = Izm_rib_top1 + Izm_rib_top2 + Izm_rib_top3 + Izm_rib_top4 + Izm_rib_top5 + Izm_rib_top6 + Izm_rib_top7 + ...
   Izm_rib_top20 + Izm_rib_top21 + Izm_rib_top22 + Izm_rib_top23 + Izm_rib_top24 + Izm_rib_top25 + Izm_rib_top26 +  ...
Izm_rib_top8 + Izm_rib_top9 + Izm_rib_top10 + Izm_rib_top11 + Izm_rib_top12 + Izm_rib_top13 + Izm_rib_top14 + Izm_rib_top15 + ...
Izm_rib_top16 + Izm_rib_top17 + Izm_rib_top18 + Izm_rib_top19;

Izm_rib_bot = Izm_rib_bot1 + Izm_rib_bot2 + Izm_rib_bot3 + Izm_rib_bot4 + Izm_rib_bot5 + Izm_rib_bot6 + Izm_rib_bot7 + ...
   Izm_rib_bot20 + Izm_rib_bot21 + Izm_rib_bot22 + Izm_rib_bot23 + Izm_rib_bot24 + Izm_rib_bot25 + Izm_rib_bot26 +  ...
Izm_rib_bot8 + Izm_rib_bot9 + Izm_rib_bot10 + Izm_rib_bot11 + Izm_rib_bot12 + Izm_rib_bot13 + Izm_rib_bot14 + Izm_rib_bot15 + ...
Izm_rib_bot16 + Izm_rib_bot17 + Izm_rib_bot18 + Izm_rib_bot19;

Ixm_rib_top = Ixm_rib_top1 + Ixm_rib_top2 + Ixm_rib_top3 + Ixm_rib_top4 + Ixm_rib_top5 + Ixm_rib_top6 + Ixm_rib_top7 + ...
   Ixm_rib_top20 + Ixm_rib_top21 + Ixm_rib_top22 + Ixm_rib_top23 + Ixm_rib_top24 + Ixm_rib_top25 + Ixm_rib_top26 +  ...
Ixm_rib_top8 + Ixm_rib_top9 + Ixm_rib_top10 + Ixm_rib_top11 + Ixm_rib_top12 + Ixm_rib_top13 + Ixm_rib_top14 + Ixm_rib_top15 + ...
Ixm_rib_top16 + Ixm_rib_top17 + Ixm_rib_top18 + Ixm_rib_top19;

Ixm_rib_bot = Ixm_rib_bot1 + Ixm_rib_bot2 + Ixm_rib_bot3 + Ixm_rib_bot4 + Ixm_rib_bot5 + Ixm_rib_bot6 + Ixm_rib_bot7 + ...
   Ixm_rib_bot20 + Ixm_rib_bot21 + Ixm_rib_bot22 + Ixm_rib_bot23 + Ixm_rib_bot24 + Ixm_rib_bot25 + Ixm_rib_bot26 +  ...
Ixm_rib_bot8 + Ixm_rib_bot9 + Ixm_rib_bot10 + Ixm_rib_bot11 + Ixm_rib_bot12 + Ixm_rib_bot13 + Ixm_rib_bot14 + Ixm_rib_bot15 + ...
Ixm_rib_bot16 + Ixm_rib_bot17 + Ixm_rib_bot18 + Ixm_rib_bot19;


Iym_rib_top = (x_rib_top^2 + z_rib_top^2 ) * V_rib * rho_balsa; 				 % (slug ft^2)
Iym_rib_bot = (x_rib_bot^2 + x_rib_bot^2) * V_rib * rho_balsa; 				 % (slug ft^2)


Izm_wing = Izm_spar_top + Izm_rib_top + Izm_spar_bot + Izm_rib_bot% (slug ft^2)
Iym_wing = Iym_spar_top + Iym_rib_top + Iym_spar_bot + Iym_rib_bot% (slug ft^2)
Ixm_wing = Ixm_spar_top + Ixm_rib_top + Ixm_spar_bot + Ixm_rib_bot% (slug ft^2)

%Moment of inertia for fuselage

Izc_fuselage1a = (x_fuselage_front * z_fuselage1a^3)/36;% (ft^4) 
Iyc_fuselage1a = (z_fuselage1a * x_fuselage_front^3)/48;% (ft^4) 
Izc_fuselage1b = (x_fuselage_front * z_fuselage1b^3)/12;% (ft^4)
Iyc_fuselage1b = (z_fuselage1b * x_fuselage_front^3)/12;% (ft^4)
Izc_fuselage1c = (x_fuselage_front * z_fuselage1c^3)/36;% (ft^4)
Iyc_fuselage1c = (z_fuselage1c * x_fuselage_front^3)/48;% (ft^4)

Izc_fuselage2 = (x_fuselage2 * z_fuselage^3)/12; % (ft^4)
Iyc_fuselage2 = (z_fuselage * x_fuselage2^3)/12; % (ft^4)

Izc_fuselage3 = (x_fuselage3 * z_fuselage^3)/36;  % (ft^4)
Iyc_fuselage3 = (z_fuselage * x_fuselage3^3)/48;  % (ft^4) 


%Mass moment of inertia for the fuselage

Izm_fuselage1a = ((xc_fuselage1a - xc_sid5)^2) * rho_balsa* V_fuselage1a;% (slug ft^2) 
Ixm_fuselage1a = ((zc_fuselage1a - zc_sid5)^2) * rho_balsa* V_fuselage1a;% (slug ft^2) 
Iym_fuselage1a = ((xc_fuselage1a - xc_sid5)^2 + (zc_fuselage1a - zc_sid5)^2) * rho_balsa* V_fuselage1a;% (slug ft^2)
Izm_fuselage1b = ((xc_fuselage1b - xc_sid5)^2) * rho_balsa* V_fuselage1b;% (slug ft^2)
Ixm_fuselage1b = ((zc_fuselage1b - zc_sid5)^2) * rho_balsa* V_fuselage1b;% (slug ft^2)
Iym_fuselage1b = ((xc_fuselage1b - xc_sid5)^2 + (zc_fuselage1b - zc_sid5)^2) * rho_balsa* V_fuselage1b;% (slug ft^2)
Izm_fuselage1c = ((xc_fuselage1c - xc_sid5)^2) * rho_balsa* V_fuselage1c;% (slug ft^2)
Ixm_fuselage1c = ((zc_fuselage1c - zc_sid5)^2) * rho_balsa* V_fuselage1c;% (slug ft^2)
Iym_fuselage1c = ((xc_fuselage1c - xc_sid5)^2 + (zc_fuselage1c - zc_sid5)^2) * rho_balsa* V_fuselage1c;% (slug ft^2)

Izm_fuselage2 = ((xc_fuselage2 - xc_sid5)^2) * rho_balsa * V_fuselage2;% (slug ft^2)
Ixm_fuselage2 = ((zc_fuselage2 - zc_sid5)^2) * rho_balsa * V_fuselage2;% (slug ft^2)
Iym_fuselage2 = ((xc_fuselage2 - xc_sid5)^2 + (zc_fuselage2 - zc_sid5)^2) * rho_balsa * V_fuselage2;% (slug ft^2)

Izm_fuselage3 =((xc_fuselage3 - xc_sid5)^2) * rho_balsa * V_fuselage3;% (slug ft^2)
Ixm_fuselage3 =((zc_fuselage3 - zc_sid5)^2) * rho_balsa * V_fuselage3;% (slug ft^2)
Iym_fuselage3 =((xc_fuselage3 - xc_sid5)^2 + (zc_fuselage3 - zc_sid5)^2) * rho_balsa * V_fuselage3;% (slug ft^2)

Izm_fuselage = Izm_fuselage1a + Izm_fuselage1b + Izm_fuselage1c + Izm_fuselage2 + Izm_fuselage3;% (slug ft^2)
Ixm_fuselage = Ixm_fuselage1a + Ixm_fuselage1b + Ixm_fuselage1c + Ixm_fuselage2 + Ixm_fuselage3;% (slug ft^2)
Iym_fuselage = Iym_fuselage1a + Iym_fuselage1b + Iym_fuselage1c + Iym_fuselage2 + Iym_fuselage3;% (slug ft^2)


%Moment of inertia for H-tail
Izc_cap_htail = (x_sparcap_htail * z_sparcap_htail^3)/12;% (ft^4)
Iyc_cap_htail = (z_sparcap_htail * x_sparcap_htail^3)/12;% (ft^4)
Izc_web_htail = (x_shearweb_htail * z_shearweb_htail^3)/12; % (ft^4)
Iyc_web_htail = (z_shearweb_htail * x_shearweb_htail^3)/12; % (ft^4)


%Mass moment of inertia for H-tail

Izm_web_htail = ((x_htail)^2)* rho_balsa * V_web_htail; % (slug ft^2)
Ixm_web_htail = ((z_htail)^2)* rho_balsa * V_web_htail; % (slug ft^2)
Iym_web_htail = ((x_htail)^2 + (z_htail)^2)* rho_balsa * V_web_htail; % (slug ft^2)
Izm_cap_htail = ((x_htail)^2)* rho_spruce * V_cap_htail;% (slug ft^2)
Ixm_cap_htail = ((z_htail)^2)* rho_spruce * V_cap_htail;% (slug ft^2)
Iym_cap_htail = ((x_htail)^2 + (z_htail)^2)* rho_spruce * V_cap_htail;% (slug ft^2)

Izm_spar_htail = Izm_web_htail + Izm_cap_htail;			  % (slug ft^2)
Iym_spar_htail = Iym_web_htail + Iym_cap_htail;			  % (slug ft^2)
Ixm_spar_htail = Ixm_web_htail + Ixm_cap_htail;			  % (slug ft^2)

y_htrib1 = 3/12;
y_htrib2 = 6/12;
y_htrib3 = 9/12;
y_htrib4 = 12/12;
y_htrib5 = 15/12;
y_htrib6 = 18/12;
y_htrib7 = 21/12;
y_htrib8 = -3/12;
y_htrib9 = -6/12;
y_htrib10 = -9/12;
y_htrib11= -12/12;
y_htrib12 = -15/12;
y_htrib13 = -18/12;
y_htrib14 = -21/12;

Izm_rib_htail = (((x_htail)^2) * (rho_balsa * (V_rib_htail * n_ribs_tail)));% (slug ft^2)

Ixm_rib_htail1 = (((z_htail)^2 + (y_htrib1)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail2 = (((z_htail)^2 + (y_htrib2)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail3 = (((z_htail)^2 + (y_htrib3)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail4 = (((z_htail)^2 + (y_htrib4)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail5 = (((z_htail)^2 + (y_htrib5)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail6 = (((z_htail)^2 + (y_htrib6)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail7 = (((z_htail)^2 + (y_htrib7)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail8 = (((z_htail)^2 + (y_htrib8)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail9 = (((z_htail)^2 + (y_htrib9)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail10 = (((z_htail)^2 + (y_htrib10)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail11= (((z_htail)^2 + (y_htrib11)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail12 = (((z_htail)^2 + (y_htrib12)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail13 = (((z_htail)^2 + (y_htrib13)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)
Ixm_rib_htail14 = (((z_htail)^2 + (y_htrib14)^2) * (rho_balsa * V_rib_htail));% (slug ft^2)

Ixm_rib_htail = Ixm_rib_htail1 + Ixm_rib_htail2 + Ixm_rib_htail3 + Ixm_rib_htail4 + Ixm_rib_htail5 + Ixm_rib_htail6 +...
   Ixm_rib_htail7 + Ixm_rib_htail8 + Ixm_rib_htail9 + Ixm_rib_htail10 + Ixm_rib_htail11 + Ixm_rib_htail12 +...
Ixm_rib_htail13 + Ixm_rib_htail14;
Iym_rib_htail = (((x_htail)^2 + (z_htail)^2) * rho_balsa * V_rib_htail);% (slug ft^2)

%Moment of inertia for V-tail
Izc_vtail_tri = (b_vtail * z_vtail^3)/36;    % (ft^4)
Iyc_vtail_tri = (z_vtail * b_vtail^3)/48;    % (ft^4)
Izc_vtail_rec = (x_vtail_tip * z_vtail^3)/12;% (ft^4)
Iyc_vtail_rec = (z_vtail * x_vtail_tip^3)/12;% (ft^4)

%Mass moment of inertia for V-tail
Izm_vtail = ((x_vtail)^2)* rho_balsa * V_vtail;              % (slug ft^2)
Ixm_vtail = ((z_vtail)^2)* rho_balsa * V_vtail;              % (slug ft^2)
Iym_vtail = ((x_vtail)^2 + (z_vtail)^2)* rho_balsa * V_vtail;% (slug ft^2)

%Mass moment of inertia for internal equipment

Izm_gearbox =  ((xc_gearbox - xc_sid5)^2)  * weight_gearbox;				                 				   % (slug ft^2)
Ixm_gearbox =  ((zc_gearbox - zc_sid5)^2)  * weight_gearbox;				                 				   % (slug ft^2)
Iym_gearbox =  ((xc_gearbox - xc_sid5)^2 + (zc_gearbox - zc_sid5)^2) * weight_gearbox;					   % (slug ft^2)
Izm_motor = ((xc_motor - xc_sid5)^2) * weight_motor;																   % (slug ft^2)
Ixm_motor = ((zc_motor - zc_sid5)^2) * weight_motor;																   % (slug ft^2)
Iym_motor = ((xc_motor - xc_sid5)^2+ (zc_motor - zc_sid5)^2) * weight_motor;							 	   % (slug ft^2)
Izm_speedcontrol = ((xc_speedcontrol - xc_sid5)^2) * weight_speedcontrol;									   % (slug ft^2)
Ixm_speedcontrol = ((zc_speedcontrol - zc_sid5)^2) * weight_speedcontrol;									   % (slug ft^2)
Iym_speedcontrol = ((xc_speedcontrol - xc_sid5)^2 + (zc_speedcontrol - zc_sid5)^2) * weight_speedcontrol;% (slug ft^2)
Izm_receiver = ((xc_receiver - xc_sid5)^2) * weight_receiver;				 									   % (slug ft^2)
Ixm_receiver = ((zc_receiver - zc_sid5)^2) * weight_receiver;				 									   % (slug ft^2)
Iym_receiver = ((xc_receiver - xc_sid5)^2 + (zc_receiver - zc_sid5)^2) * weight_receiver;				   % (slug ft^2)
Izm_gyro = ((xc_gyro - xc_sid5)^2) * weight_gyro;								 								 	   % (slug ft^2)
Ixm_gyro = ((zc_gyro - zc_sid5)^2) * weight_gyro;								 								 	   % (slug ft^2)
Iym_gyro = ((xc_gyro - xc_sid5)^2 + (zc_gyro - zc_sid5)^2) * weight_gyro;							   	   % (slug ft^2)
Izm_datarecord = ((xc_datarecord - xc_sid5)^2) * weight_datarecord;		 									   % (slug ft^2)
Ixm_datarecord = ((zc_datarecord - zc_sid5)^2) * weight_datarecord;		 									   % (slug ft^2)
Iym_datarecord = ((xc_datarecord - xc_sid5)^2 + (zc_datarecord - zc_sid5)^2) * weight_datarecord;	   % (slug ft^2)
Izm_battery = ((xc_battery - xc_sid5)^2) * weight_battery;				    									   % (slug ft^2)
Ixm_battery = ((zc_battery - zc_sid5)^2) * weight_battery;				    									   % (slug ft^2)
Iym_battery = ((xc_battery - xc_sid5)^2 + (zc_battery - zc_sid5)^2) * weight_battery;					   % (slug ft^2)
Izm_servo = ((xc_servo - xc_sid5)^2) * weight_servo;							 									   % (slug ft^2)
Ixm_servo = ((zc_servo - zc_sid5)^2) * weight_servo;							 									   % (slug ft^2)
Iym_servo = ((xc_servo - xc_sid5)^2 + (zc_servo - zc_sid5)^2) * weight_servo;								   % (slug ft^2)
Izm_interface = ((xc_interface - xc_sid5)^2) * weight_interface;			 									   % (slug ft^2)
Ixm_interface = ((zc_interface - zc_sid5)^2) * weight_interface;			 									   % (slug ft^2)
Iym_interface = ((xc_interface - xc_sid5)^2 + (zc_interface - zc_sid5)^2) * weight_interface;		   % (slug ft^2)


%Mass moment of inertia for Sid5

Izm_sid5 = (Izm_wing + Izm_fuselage + Izm_spar_htail + Izm_vtail + Izm_rib_htail + Izm_gearbox + Izm_motor + Izm_speedcontrol +...
   Izm_receiver + Izm_gyro + Izm_datarecord + Izm_battery + Izm_servo + Izm_interface)   % slug ft^2

Ixm_sid5 = (Ixm_wing + Ixm_fuselage + Ixm_spar_htail + Ixm_vtail + Ixm_rib_htail + Ixm_gearbox + Ixm_motor + Ixm_speedcontrol +...
   Ixm_receiver + Ixm_gyro + Ixm_datarecord + Ixm_battery + Ixm_servo + Ixm_interface)   % slug ft^2

Iym_sid5 = (Iym_wing + Iym_fuselage + Iym_spar_htail + Iym_vtail + Izm_rib_htail + Iym_gearbox + Iym_motor + Iym_speedcontrol + ...
   Iym_receiver + Iym_gyro + Iym_datarecord + Iym_battery + Iym_servo + Iym_interface)    % slug ft^2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Parallel-axis theorem for wing
Iz_capnew = Izc_cap + (A_sparcap)*((zc_cap - zc_spar)^2); % (ft^4)
Ix_capnew = Ixc_cap + (A_sparcap)*((xc_cap - xc_spar)^2); % (ft^4)
Iz_webnew = Izc_web + (A_shearweb)*((zc_web - zc_spar)^2);% (ft^4)
Ix_webnew = Ixc_web + (A_shearweb)*((xc_web - xc_spar)^2);% (ft^4)

Iz_total = 2 * Iz_capnew + Iz_webnew;% (ft^4)
Ix_total = 2 * Ix_capnew + Ix_webnew;% (ft^4)

%Parallel-axis theorem relating wing moment to plane's cg

Iz_wing1 = Iz_total + (A_spar) * ((zc_sid5 - zc_spar)^2);% (ft^4)
Ix_wing1 = Ix_total + (A_spar) * ((xc_spar - xc_sid5)^2);% (ft^4)
Iz_wing2 = Iz_total + (A_spar) * ((zc_sid5 - zc_spar)^2);% (ft^4)
Ix_wing2 = Ix_total + (A_spar) * ((xc_sid5 - xc_spar)^2);% (ft^4)

Iz_wing = Iz_wing1 + Iz_wing2;% (ft^4)
Ix_wing = Ix_wing1 + Ix_wing2;% (ft^4)

%Parallel-axis theorem for fuselage

Iz_fuselage1a_new = Izc_fuselage1a + (A_fuselage1a)*((zc_sid5 - zc_fuselage1a)^2); 			% (ft^4)
Ix_fuselage1a_new = Ixc_fuselage1a + (A_fuselage1a)*((xc_fuselage1a - xc_sid5 )^2);		   % (ft^4)
Iz_fuselage1b_new = Izc_fuselage1b + (A_fuselage1b)*((zc_sid5 - zc_fuselage1b)^2); 			% (ft^4)
Ix_fuselage1b_new = Ixc_fuselage1b + (A_fuselage1b)*((xc_fuselage1b - xc_sid5 )^2);			% (ft^4)
Iz_fuselage1c_new = Izc_fuselage1c + (A_fuselage1c)*((zc_sid5 - zc_fuselage1c)^2); 			% (ft^4)
Ix_fuselage1c_new = Ixc_fuselage1c + (A_fuselage1c)*((xc_sid5 - xc_fuselage1c)^2); 			% (ft^4)
Iz_fuselagemid_new = Izc_fuselagemid + (Area_fuselage2)*((zc_fuselagemid - zc_sid5)^2); 	% (ft^4)
Ix_fuselagemid_new = Ixc_fuselagemid + (Area_fuselage2)*((xc_fuselagemid - xc_sid5)^2); 	% (ft^4)
Iz_fuselageback_new = Izc_fuselageback + (Area_fuselage3)*((zc_fuselageback - zc_sid5)^2);% (ft^4)
Ix_fuselageback_new = Ixc_fuselageback + (Area_fuselage3)*((xc_fuselageback - xc_sid5)^2);% (ft^4)

Iz_fuselage = Iz_fuselage1a_new + Iz_fuselage1b_new + Iz_fuselage1c_new + Iz_fuselagemid_new + Iz_fuselageback_new;% (ft^4)
Ix_fuselage = Ix_fuselage1a_new + Ix_fuselage1b_new + Ix_fuselage1c_new + Ix_fuselagemid_new + Ix_fuselageback_new;% (ft^4)


%Parallel-axis theorem for H-tail
Iz_capnew_htail = Izc_cap_htail + (A_sparcap_htail)*((zc_cap_htail - zc_spar_htail)^2);  % (ft^4)
Iy_capnew_htail = Iyc_cap_htail + (A_sparcap_htail)*((yc_cap_htail - yc_spar_htail)^2);  % (ft^4)
Iz_webnew_htail = Izc_web_htail + (A_shearweb_htail)*((zc_web_htail - zc_spar_htail)^2); % (ft^4)
Iy_webnew_htail = Iyc_web_htail + (A_shearweb_htail)*((yc_web_htail - yc_spar_htail)^2); % (ft^4)

Iz_htail_total = Iz_capnew_htail + Iz_webnew_htail;  % (ft^4)
Iy_htail_total = Iy_capnew_htail + Iy_webnew_htail;  % (ft^4)

%Parallel-axis theorem relating H-tail moment to plane's cg
Iz_htail = Iz_htail_total + (A_htail)*((zc_sid5 - zc_spar_htail)^2);% (ft^4)
Iy_htail = Iy_htail_total + (A_htail)*((yc_sid5 - yc_spar_htail)^2);% (ft^4)

%Parallel-axis theorem for V-tail relating V-tail moment to plane's cg

Iz_vtailnew_tri = Izc_vtail_tri + (A_vtail_tri) * ((zc_vtail_tri - zc_sid5)^2);% (ft^4)
Iy_vtailnew_tri = Iyc_vtail_tri + (A_vtail_tri) * ((yc_vtail_tri - yc_sid5)^2);% (ft^4)
Iz_vtailnew_rec = Izc_vtail_rec + (A_vtail_rec) * ((zc_vtail_rec - zc_sid5)^2);% (ft^4)
Iy_vtailnew_rec = Iyc_vtail_rec + (A_vtail_rec) * ((yc_vtail_rec - yc_sid5)^2);% (ft^4)

Iz_vtail = Iz_vtailnew_tri + Iz_vtailnew_rec;% (ft^4)
Iy_vtail = Iy_vtailnew_tri + Iy_vtailnew_rec;% (ft^4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Moment of inertia for sid5

Iz_sid5 = (Iz_vtail + Iz_htail + Iz_fuselage + Iz_wing)  % (ft^4)
Iy_sid5 = (Iy_vtail + Iy_htail + Iy_fuselage + Iy_wing)  % (ft^4)


%Wing loading
dl=0.01;% (ft)
l=0:dl:.5*b;% (ft) - l is the postion along the span

%Trapezoidal method
q1o = c;											 % (lbf/ft)
q1 = q1o*(1-(2.*l./b).*(1-lambda))*1.89;% (lbf/ft)
%q1o=w/4*3.75/(b/2);
%q1=q1o/(b/2);

%Elliptical method
q2o = (4*c)/(pi);						  % (lbf/ft)
q2=q2o*sqrt(1-(2.*(l)./b).^2)*1.89;% (lbf/ft)

%Average of Trapezoidal and Elliptical wing loadings
qave=(q1+q2)./2;% (lbf/ft)

figure(1)
plot(l,q1,'m',l,q2,'k',l,qave,'b')
title('Wing loading using Trapezoidal and Elliptical methods')
xlabel('Distance from wing root (ft)')
ylabel('Wing loading (lbf/ft)')
legend('Trapezoidal','Elliptical','Average')
axis([0 b/2 0 max(q2)+.5])

%Setup indecies for integration
d = size(l);
n = d(2);

%Shear force

V=cumtrapz(l,qave);% (lbf)

figure(2)
plot(b/2-l,V)
title('Shear force vs distance from wing root')
xlabel('Distance from wing root (ft)')
ylabel('Shear force (lbf)')
axis([0 b/2 0 max(V)])

%found from integrating the wing loading


%Moment

M=cumtrapz(l,V);% (lbf * ft)

figure(3)
plot(b/2-l,M)
title('Moment vs distance from wing root')
xlabel('Distance from wing root (ft)')
ylabel('Moment (lbf * ft)')
axis([0 b/2 0 max(M)])


%Wing root bending analysis
dy = .01;% (ft)
y = -c*tc/2:dy:c*tc/2;% (ft)
sigma_spar = (M(d(2))*y./Iz_total)./144;% (lbf/in^2)

figure(4)
plot(sigma_spar,y*12,sigma_spar,(c*tc/2-z_sparcap)*12,'*r',sigma_spar,(-c*tc/2+z_sparcap)*12,'*r')
title('Height of spar vs normal stress')
xlabel('Normal stress (lbf/in^2)')
ylabel('Height (in)')
axis([min(sigma_spar) max(sigma_spar) min(y)*12 max(y)*12])
