%ORION-Fall 2000
%AAE451-sizing code
%Last Updated 20 SEP 2000

clear all
close all

%setting initial test values
gamma_test=0;
power_test=0;
weight_test=0;
energy_test=0;
thrust_test = 0;
test=0;
Power_thrust=0;

%Here initial estimates are made.

%control system weight estimates
n_servo=4; %number of servos required
weight_servo=2.0; %weight of a servoounces
weight_receiver=1.0; %weight of receiver in ounces
weight_esc=3.0; %weight of electronic spped control in ounces
weight_gyro=3.5; %oz
weight_tattletale8=15; %ounces

%Motor Parameters %MaxCim MaxN32-15D
Power_output=850; %rated output Watts
weight_motor=7.5; %estimate taken from MaxCim data
weight_gearbox=1.5; %oz
weight_propeller=2; %oz
cells_min=16; %the minimum number of cells the motor can handle
cells_max=16; %the maximum number of cells the motor can handle

%Battery Pack
cells=14; %number of cells-will change with power needed
cell_weight=2.8; %ounces-const once battery is chosen
cell_voltage=1.2; %Volts-const once battery is chosen
cell_capacity=2.1; %Amp-hours-const once battery is chosen

%Propulsion system efficiencies
eta_propeller=.7; %propeller efficiency
eta_gearbox=.95; %gear box efficiency
eta_motor=.9; %motor overall efficiency
eta_esc=.95; %electronic speed control efficency
eta_battery=.95; %battery efficency

%Cruise Speed
Vcruise=25; %ft/s. Desired Vcruise 
Vmax=30; %ft/s. Maximum flight speed in Mollenkop

%Stall Speed
Vstall_max=20; %(ft/s)
Vstall_desired=20; %(ft/s)

%Turn radius
radius_max=50; %ft
radius_desired=30; %ft

%Climb angle minimum
gamma_min=12; %degrees
%Wing Aspect Ratio
AR=2.2; %this value is estimated from historical data and subject to change

%Drag estimates
CLmax=1.1; %initial estimation of CLmax - 3D
Cf=.005; %this is an estimate taken from Brandt-Intro to Aero Design    
mu=1; %shorter wing span/longer wing span
r=1; %lift on shorter wing/lift on longer wing
sigma=0.5; %Prandtl's interference factor (gap/average span) ranges from .12 to 1

%Wetted Area estimates
CHT=0.45; %horizontal tail volume coefficient-based on homebuilt data
CVT=0.04; %vertical tail volume coefficient-based on homebuilt data
fuselage_height=.5; %ft
fuselage_width=.417; %ft
d_tire=.208; %ft-diamter of the tire
t_tire=.04167; %ft-thickness of the tire

%Take off distance
Dto_max=120; %ft
Dto=30; %ft
mhu=.9;%breaking coefficient:1.5*(dry concrete) was estimated for astroturf
%breaking coefficient:0.4-0.6:dry concrete,0.2-0.3:wet concrete,0.05-0.1:icy runway

%Climb inputs
cruise_alt=20; %ft

%Landing distance
distance_allow=120; %ft allowed to role before coming to a stop

%Atmospheric Data Estimates. Purdue is 607ft above sealevel
Patm=2041; %Press at 1000ft in lb/ft^2; input('Input a pressure for Mollenkopf in Dec.');
Tatm=515.1; %Temp at 1000ft in R.; input('Input a temperature for Mollenkopf in Dec.');
Ru=1545; %ft-lbf/lb-mol-R
mol=28.97; %lb/lb-mol air
g=32.174; %(ft/s^2)
Tsl=518.69; %R
Psl=2116.2; %lb/ft^2
rho=(.002308);

%Calculations begin

lbffts_Watt=1.3558; %conversion from lbf-ft/s to Watts

%calculate motor
eta_total=eta_motor*eta_gearbox*eta_propeller*eta_esc*eta_battery; %total porpulsion efficiency
Power_avail=Power_output*eta_total; %Watts
%Thrust_avail=(.1529*Power_avail+1.2473)/16; %lbf. (was divided by 16 b/c data was given in ounces)
   
%Weight calculations
weight_battery=cells*cell_weight; %ounces-initial estimate
weight_payload=n_servo*weight_servo+weight_receiver+weight_esc+weight_tattletale8+weight_gyro; %ounces-estimated will change with more concrete component weights
weight_empty=1.5*(weight_motor+weight_gearbox+weight_propeller+weight_battery); %inital calculation of empty weight based on AstroFlight
W=1.2263*weight_empty-0.1574; %initial take-off weight based on WTA database

while test == 0 
   
  
   while weight_test == 0
      
	   %Wing Area
	   %Calculating Wing Area for Stall
   	S_stall=(W/16)/(CLmax*rho*(Vstall_desired^2)*.5); %ft^2
	   %Calculating Wing Area as a function of turn
   	S_turn=(W/16)/(CLmax*rho*(Vmax^2)*.5); %ft^2
   
   		if S_stall < S_turn
		      S=S_turn;
		   else
      		S=S_stall;
		   end
     
      %Wing_span
   	b=sqrt(S*AR); %ft
   
		%   if b > 8
      %      disp('span too large, decrease AR')
      %      pause
		%   end
   
   	%chord
	   c=(S/2)/b; %ft
   	t=c*.14; %ft-establishes a t/c of .14.-subject to change
   
   	%Wetted Area Calculations
	   fuselage_length=.88*b; %ft-historical S vs fuselage length data
   	LVT=0.7*fuselage_length; %moment arm of vertical tail
	   LHT=0.7*fuselage_length; %moment arm of horizontal tail
   	SVT=CVT*b*S/LVT; %ft^2-area of vertical tail
	   SHT=CHT*c*S/LHT; %ft^2-area of the horizontal tail
      
      %Horizontal Tail Parameters
      lamdah=0.6; %found from Raymer p85
      ARh=3; %found from Raymer p85
      b_htail=sqrt(ARh*SHT);
      Croot_htail=2*SHT/(b_htail*(1+lamdah));
      Ctip_htail=lamdah*Croot_htail;
      Cbar_htail=(2/3)*Croot_htail*(1+lamdah+(lamdah^2))/(1+lamdah);
      sweep_leh=atan((Croot_htail-Ctip_htail)/(b_htail/2));%leading edge sweep (rad)
      sweep_qch=atan((tan(sweep_leh))- ((1-lamdah)/(ARh*(1+lamdah)))); %sweep for the quarter chord (rad)
       
      %Vertical Tail Paramters
      ARv=1.3; %from Raymer 185
      lamdav=0.6;
      b_vtail=sqrt(ARv*SVT);
      Croot_vtail=2*SVT/(b_vtail*(1+lamdav));
      Ctip_vtail=lamdav*Croot_vtail;
      Cbar_vtail=(2/3)*Croot_vtail*(1+lamdav+(lamdav^2))/(1+lamdav);
      sweep_lev=atan((Croot_vtail-Ctip_vtail)/(b_vtail));%leading edge sweep
      sweep_qcv=atan((tan(sweep_lev))- ((1-lamdav)/(ARv*(1+lamdav)))); %sweep for the quarter chord
      
      Swet_fuse=2*fuselage_height*fuselage_width+2*(fuselage_height+fuselage_width)*fuselage_length; %ft^2
	   Swet_wing=(S-c*fuselage_width)*(1.977+.52*(t/c)); %ft^2
   	Swet_tail=2*SHT+2*SVT; %note thickness of the tail is negligable
	   Swet_landgear=2*(pi*((d_tire/2)^2)*t_tire); %ft^2
   	Swet_struts=11*(2*(.00833*.04167)+(.04167*c)+(.00833*c)); %ft^2
   
   	%Calculate weight based on area
      weight_wing=2.75*S-1.9564; %ounces Orion Database
      weight_tail=2.75*(SHT+SVT)-1.9564; %ounces Orion Database
      thickness_fuselage=(1/8)/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)*16; %(lb)->oz
      weight_struts=(Swet_struts/2*0.00833*spruce_density)*16; %(lb)->oz
      weight_mainwheel=1.28; %oz 2.5" tire
      weight_tailwheel=.5; %oz
      weight_landinggear=5; %oz
      Wnew=weight_wing+weight_tail+weight_fuselage+2*weight_mainwheel+weight_tailwheel+weight_landinggear+weight_payload+weight_battery; %oz
   
   	if abs(W-Wnew) < 0.001 %ounce difference
     	   weight_test=1;
	   else
%         W
%         pause
  	  	   weight_test=0;
         W=Wnew;
	   end
 
   end %End structure loop
   
   disp('Passed Structure')
   pause  
   


   %Total wetted area
   Swet=Swet_fuse+Swet_wing+Swet_tail+Swet_landgear+Swet_struts; %ft^2
   
   %Vehicle Drag Estimate
   e=((mu^2)*(1+r^2))/((mu^2)+(2*sigma*mu*r)+(r^2)); %Oswald efficiency factor for a bi-plane
   k=1/(pi*AR*e);%taken from Roskam Airplane Aerodynamics and performance
   	%Calculation of Cdo
      Vto=1.2*sqrt((2*((W/16)/S))/(rho*CLmax));
      
      V=[Vto Vmax Vcruise]; %This will cause Cdo to be calculated for the three velocities 
      							 %seen during the mission.  	
      Re_perlength=.002308*V./3.717E-7; %Reynolds number divided by length

		%Wing
   
			%Skin Friction   
			Re_wing=Re_perlength.*c;

			Cf_wing=(0.455./((log10(Re_wing)).^2.58));
   
   		%Form Factor
      	Zw=2;
      
      	Kw=1+Zw.*(t./c)+100.*((t./c).^4);
      
      	%Interference Factor
      	Qw=1.5;	%Estimated 
      
      
      	%Drag Coefficient
      	Cd_w=Kw.*Qw.*Cf_wing.*Swet_wing./S;
   
   
		%Fuselage
   
			%Skin Friction   
			Re_fuselage=Re_perlength.*fuselage_length;

			Cf_fuselage=(0.455./((log10(Re_fuselage)).^2.58));

			%Form Factor
			lamdaf=(fuselage_length./fuselage_height); %fineness ratio
      
      	Kf=1+(60./(lamdaf.^3))+(lamdaf./400);
      
      	%Interference Factor
      	Qf=1; %estimated
      
      	%Drag Coefficient
      	Cd_f=Kf.*Qf.*Cf_fuselage.*Swet_fuse./S;
      
   	%Vertical Tail
   
   		%Skin Friction
			Re_vtail=Re_perlength.*Cbar_vtail;

			Cf_vtail=(0.455./((log10(Re_vtail)).^2.58));
      
   	   %Form Factor
      	Kvtail=(2.*cos(sweep_qcv));
              
      	%Interference Factor
      	Qvtail=1;	%Estimated 
      
      
      	%Drag Coefficient
      	Cd_vtail=Kvtail.*Qvtail.*Cf_vtail.*2.*SVT./S;
      
		
		%Horizontal Tail

			%Skin Friction

			Re_htail=Re_perlength.*Cbar_htail;

			Cf_htail=(0.455./((log10(Re_htail)).^2.58));
   
   	   %Form Factor
      	Khtail=(2.*cos(sweep_qch));
              
      	%Interference Factor
      	Qhtail=1;	%Estimated 
      
      
      	%Drag Coefficient
	      Cd_htail=Khtail.*Qhtail.*Cf_htail.*2.*SHT./S;

	%Miscellaneous Drag Components Calculated in accordance with Raymer p.346-348

		%Struts
   
   	Frontal_Area_St=11.*(c.*8.33E-3);
   
   	Factor_St=0.4;	%D/q / Frontal Area for rectangular strut estimated from Raymer p. 349
   
   	Cd_St=(Frontal_Area_St.*Factor_St)./S;
   
   	Cd_St=1.2.*Cd_St; %Raymer p. 347, stimate that accounts for mutual interference
   
   	%Landing Gear
   
   	Frontal_Area_Lg=2.*(d_tire.*t_tire);
   
   	Factor_Lg=0.25; %D/q / Frontal Area for regular tire from Raymer p. 349

		Cd_Lg=(Frontal_Area_Lg.*Factor_Lg)./S;

		Cd_Lg=1.2.*Cd_Lg; %%Raymer p. 347, stimate that accounts for mutual interference
   
   	%Total Drag excluding Induced Drag
   
   	CDo=Cd_w+Cd_f+Cd_vtail+Cd_htail+Cd_Lg+Cd_St;
   
   
   	%Extra Drag that accounts for Leakage and Protuberance Drag and Upswept fuselage
   
   	Cd_extra=0.15.*CDo;
   
   	%Revaluate Total Drag excluding induced drag to account for extra misc. drag
   
   	CDo=CDo+Cd_extra;
      CDo_Vto=CDo(1);
      CDo_Vmax=CDo(2);
      CDo_Vcruise=CDo(3);
   
   
   
   CD=CDo_Vmax+(CLmax^2)/(pi*AR*e); %Roskam
   Drag_max=(CD*S)*rho*(Vmax^2)*.5;
   Lift_max=CLmax*rho*(Vmax^2)*(S/2);
   
   %Take-Off Calculations
   
   %Take off Speed
   Vto=1.2*sqrt((2*((W/16)/S))/(rho*CLmax));
   Drag_to=(CDo_Vto*S)*.5*rho*(.7*Vto)^2;
   L_to=0;
   %L_to=CLmax*rho*((.7*Vto)^2)*(S/2);
   %Take off Thrust 
   
   Thrust_to=1.44*((W/16)^2)/(g*CLmax*Dto*rho*S)+0.05*((W/16)-L_to)+Drag_to; %lbf
   %Need to find Power from Thrust
   Power_to=(Thrust_to*Vto)*lbffts_Watt; %Watts
   %Need to find Energy from Power
   time_to=(1.2*sqrt((2*(W/16))/(rho*S*CLmax)))/((g*(Thrust_to-Drag_to-0.05*((W/16)-L_to)))/(W/16)); %sec
   Energy_to=(Power_to*time_to)/3600; %Watt-hr   
   
   %Climb Calculations
   
   hdot_min=Vto*sin(gamma_min*pi/180); %ft/s-climb rate minimum
   CL_climb=sqrt(3*CDo_Vto/k); %Roskam 3.25
   CD_climb=4*CDo_Vto; %Roskam 3.26                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
   Drag_climb=(CD_climb*S)/(2*rho*(Vto)^2); %currently it is assumed that climb is at Vto
   %gammanew=asin((Thrust_avail-Drag_climb)/(W/16))*(180/pi);
   Roc=Vto*sin(gamma_min*pi/180); %ft/s-rate of climb  
   L_climb=(W/16)*cos(gamma_min*pi/180); %lbf
   %Climb Thrust
   Thrust_climb=(W/16)*sin(gamma_min*pi/180)+Drag_climb; %lbf
   %Need to find Power from Thrust
   
   %Power_climb=((eta_propeller./(((Vto.*sin(gamma_min.*pi/180))./550)+(((W/16)/S)^(1/2))./(19.*(CL_climb.^3/2)./(CD_climb))))^(-1))*(W/16);

   Power_climb=(Thrust_climb*Vto)*lbffts_Watt; %Watts
   Power_climb_poo=2*Power_climb/eta_total;
   %Need to find Energy from Power
   time_climb=cruise_alt/Roc; %sec
   Energy_climb=(Power_climb_poo*time_climb)/3600; %Watt-hr
   
   %Cruise/Loiter Calculations   
   
   %Cruise Thrust
   CL_cruise=(W/16)/(.5*rho*(Vcruise^2)*S);
   CD_cruise=CDo_Vcruise+k*CL_cruise^2;
   Thrust_cruise=((CD_cruise*S)*(.5*rho*(Vcruise)^2)); %lbf
   %Need to find Power from Thrust
   Power_cruise=(Thrust_cruise*Vcruise)*lbffts_Watt; %Watts
   %Need to find Energy from Power
   time_cruise=648; %sec
   Energy_cruise=(Power_cruise*time_cruise)/3600; %Watt-hr
   
   %Turning Calculations
   
   %Turn Rate
   psidot_desired=Vcruise/radius_desired; %radians/s
   psidot_max=Vcruise/radius_max; %radians/s
   n=sqrt((((psidot_desired*Vcruise)/g)^2)+1); %load factor-Raymer
   CL_turn=((n)*(W/16))/(.5*rho*(Vcruise^2)*S); %Brand 5.11.1
   CD_turn=CDo_Vcruise+k*(CL_turn^2);
   L_turn=n*(W/16); %lbf-assuming that is cruise L=W
   Drag_turn=(0.5*rho*S*CD_turn*Vcruise^2); %lbf
   Thrust_turn=Drag_turn; %lbf-Thrust=Drag in a sustained turn Brandt
   %Need to find Power from Thrust
   Power_turn=(Drag_turn*Vcruise)*lbffts_Watt; %Watts
   %Need to find Energy from Thrust
   time_turn=254.5; %sec. will need to change is Vcruise ne 25ft/s or radius ne 35.7ft.
   %this was calculated assuming about 15 minutes of flight
   Energy_turn=(Power_turn*time_turn)/3600; %Watt-hr
   
   %Calculating total energy, thrust required
   Energy_total=(Energy_to+Energy_climb+Energy_cruise+Energy_turn);
   
   Thrust_total=max([Thrust_to Thrust_climb Thrust_cruise Thrust_turn]);
   
   %check power
   Power_req=max([Power_to Power_climb Power_cruise Power_turn Power_thrust]);
   %Excess Power
   Power_excess=(Power_avail-Power_req)/(W/16); 
   
   %######################
   %	Times above calculaetd based on 37.5ft radius turn
   %	Fligh path is in Keith's note book
   %######################
   %Compute battery energy
   
   cell_energy=cell_voltage*cell_capacity; %Watt hours per cell
   battery_voltage=cells*cell_voltage; %for cells in series
   battery_energy=battery_voltage*cell_capacity*eta_total; %Energy capacity
   cells_new=ceil(Energy_total/cell_energy);
   
   cells_test=mod(cells_new,cells);
   
   if cells_test == 0
      energy_test=1;
   else if  cells == cells_min
         energy_test=1;
      else 
         if cells_new < cells_min
            disp('Cells min reached')
            energy_test=0;
            cells=cells_min;    
         else if cells_new > cells_max
               disp('Cells max reached')
%               pause
            else
               disp('Cells fault')
               energy_test=0;
%               weight_test=0;
               cells=cells_new;
%               %Weight calculations
%					weight_battery=cells*cell_weight; %ounces-initial estimate
%					weight_payload=n_servo*weight_servo+weight_receiver+weight_esc+weight_tattletale8+weight_gyro; %ounces-estimated will change with more concrete component weights
%					weight_empty=1.5*(weight_motor+weight_gearbox+weight_propeller+weight_battery); %inital calculation of empty weight based on AstroFlight
%					W=1.2263*weight_empty-0.1574; %initial take-off weight based on WTA database
            end
         end
      end
   end
   
   
	if Power_avail < Power_req
      disp('The motor is too small') 
	   pause   
   end 
   
   %check on Climb angle
%   if gammanew < gamma_min
%      disp('Climb angle fault')
%      pause
%   end %if loop ends

   %Check if the plane passes the tests
   if energy_test == 1
      test=1;
   else
      test=0;
      weight_test=0;
   end 
         
         disp('completed tests')
         pause
         
end %end while
      
disp('Passed All Tests!')
w_lb=W/16;
w_lb
weight_battery
cells
S
b
c
Power_req
Power_avail
Energy_total
battery_energy

%AAE451-ORION Constraint Diagram
%Plotting Constraint Equations to find allowable design space 
%Plotting Wing Loading (W/S) vs Power Loading (P/S)

%Equations: Endurance,Climb,Turn,Takeoff,Stall 

%Given info

q=.5*rho*S*(Vcruise^2);

WingLoading=0.2:.01:1.2;
PowerLoading=0.01:.1:90;
      
%V-stall Constraint
WingLoading_stall=((Vstall_desired).^2)*rho*CLmax./2;
figure(1)
subplot(221)     
plot(WingLoading_stall,PowerLoading)
title('Vstall Constraint')
xlabel('Wing Loading [lb/ft^2]')
ylabel('Power Loading [lb/hp]')

%Endurance Constraint
PowerLoading_end=((((S.*rho./2)*(WingLoading)).^(-.5))*((4*CDo_Vcruise).^(-1))*((3*CDo_Vcruise./k).^(3/4)))*eta_total;
subplot(222)
plot(WingLoading,PowerLoading_end)
title('Endurance Constraint')
xlabel('Wing Loading [lb/ft^2]')
ylabel('Power Loading [lb/hp]')
   
%Cruise: 
LDcr=1./((q*CDo_Vcruise./WingLoading)+(WingLoading./(q*pi*AR*e)));
TWcr = 1./LDcr;
PLcr = 1./(Vcruise.*TWcr)*550*eta_total;
subplot(223)
plot(WingLoading,PLcr)
title('Cruise Speed (V=25 ft/s) Constraint')
xlabel('Wing Loading [lb/ft^2]')
ylabel('Power Loading [lb/hp]')

%ROC
%PowerLoading_ROC=(((Roc+4*(2*WingLoading./rho).^(.5))*CDo*(3*CDo./k).^(-3/4)).^-1)*550;
PowerLoading_ROC=eta_total./(((Vto.*sin(gamma_min.*pi/180))./550)+((WingLoading.^(1/2))./(19.*(CL_climb.^3/2)./(CD_climb))));
PowerLoading_ROC2=1./(((Vto.*sin(gamma_min.*pi/180))./550)+((WingLoading.^(1/2))./(19.*(CL_climb.^3/2)./(CD_climb))));

%PwrRoc=(19*(CL_climb^(3/2))/CD_climb)./(WingLoading.^0.5);
subplot(224)
plot(WingLoading,PowerLoading_ROC)
hold on
plot(WingLoading,PowerLoading_ROC2)

%plot(WingLoading,PwrRoc,'--')
title('Rate of Climb Constraint')
xlabel('Wing Loading [lb/ft^2]')
ylabel('Power Loading [lb/hp]')

%Take Off Constaint;Ground roll
Vto=1.2*sqrt((2*(WingLoading))/(rho*CLmax));
PowerLoading_to=((1.008*Vto.*WingLoading./(rho*Dto*CLmax*g)+(0.1715*CDo_Vto*rho*(Vto.^3))./WingLoading+(0.7*.05*Vto)).^(-1))*550*eta_total;
figure (2)
subplot(211)
plot(WingLoading,PowerLoading_to,'--')
title('TO/Ground roll Constraint w/ 120ft To distance')
xlabel('Wing Loading [lb/ft^2]')
ylabel('Power Loading [lb/hp]')
%Turn Constraint
PowerLoading_turn=((Vcruise*((q*CDo_Vcruise./WingLoading)+((1.2.^2)*k*WingLoading./q))).^(-1))*550*eta_total;
subplot(212)
plot(WingLoading,PowerLoading_turn,'g')
title('Turn Contraint')
xlabel('Wing Loading [lb/ft^2]')
ylabel('Power Loading [lb/hp]')
hold off

figure(3)
hold on
axis([.2 1.2 0 90])
plot(WingLoading,PowerLoading_end,'m',WingLoading,PLcr,'g',WingLoading,PowerLoading_ROC,'b',WingLoading,PowerLoading_to,'c',WingLoading,PowerLoading_turn,'k',WingLoading_stall,PowerLoading,'r')
legend('Endurance','Cruise','Climb','Take-off','Turn','Vstall')
xlabel('Wing Loading [lb/ft^2]')
ylabel('Power Loading [lb/hp]')
title('Constraint Diagram')
gtext('Design Space')
eta=eta_motor*eta_gearbox*eta_propeller;
plot(w_lb/S,w_lb/(eta*Power_output/745.7),'p')
gtext('Design Point')
