% FILE: Weight_3.m
% Preliminary weight estimator for electric powereed aircraft
% Revised 9/5/06
disp(' '); disp('>>>>>>>>>Start here <<<<<<<<<'); disp(' ')
LoverDmax=12           % for fixed gear GA aircraft (Skyhawk) (See Raymer p. 22)
LoverD=.866*LoverDmax  % for loiter (See Raymer p. 22)
Vloiter=40             % ft/sec, Estimated loiter speed
ETAmotor=0.75
ETAprop= 0.60
%RHOb=72900    % battery energy density for NiCad joule per pound
%RHOb=9.25E+04 % battery energy density for NiMH joule per pound
RHOb=2.39E+05 % battery energy density for Lithium polymer joule per pound
disp('Battery energy density for NiCad batteries, joules per pound')
EnduranceMIN=7
Wpayload=1 % payload weight pounds

EnduranceSEC=EnduranceMIN*60
TimeLoiterStraight=EnduranceSEC/2   % Loiter time in straight flight (sec)
TimeLoiterTurn=EnduranceSEC/2       % Loiter time in turning flight (sec)
g=32.17 % acceleration of gravity ft/sec^2


% For loiter in straight flight
WlsperW=Vloiter*1.356*TimeLoiterStraight/(ETAmotor*ETAprop*RHOb*LoverD)

% For loiter in turning flight
R=50 % Turn radius at loiter from mission spec.
phi=atan(Vloiter*Vloiter/(R*g)) % bank angle in the turn (rad)
WltperW=Vloiter*1.356*TimeLoiterTurn/(ETAmotor*ETAprop*RHOb*LoverD*cos(phi))

% For climbing flight

gamma=20/57.3 % climb angle (rad)
TimeClimb=12/(Vloiter*sin(gamma)) % time to climb to 12 feet
WclimbperW=Vloiter*1.356*TimeClimb*(cos(gamma)/LoverD+sin(gamma))/(ETAmotor*ETAprop*RHOb)


% For Takeoff
disp('From integration of eoms at takeoff, assume that the battery')
disp(' weight fraction is .002.')
WtoperW=.002

% For warm-up assume takeoff times aree about 3 sec and 
% warm-up times are about 30 seconds.
disp('Assume that the warmup weight fraction is 10 times the ')
disp(' takeoff weight fraction.')
WwarmperW=10*WtoperW


% Assemble the complete battery weight fraction.
WbperW=WlsperW+WltperW+WclimbperW+WtoperW+WwarmperW


Weight=0:1:10; %weight in pounds
echo on
WminusWe=.2103*Weight+.1243; % formula for historical data (pounds)
echo off
disp('Your weight estimate will only be as good at that historical data represented in the equation above')
Wbattery=WbperW*Weight;
WbplusWpay=Wbattery+Wpayload;
plot(Weight,WminusWe,Weight,WbplusWpay)
xlabel('Weight~lbf')
ylabel('W-We and Wb+Wp~lbf')

% Determination of aircraft weight
delta=WminusWe-WbplusWpay;
% YI = INTERP1(X,Y,XI)
Waircraft=interp1(delta,Weight,0)
y=.2103*Waircraft+.1243;
string1=['Estimated aircraft weight is ',num2str(Waircraft),' pounds.']
text2(.25,.2,['           ',string1])
title('Weight estimation using historical weight data')
legend('Historical data','Estimated weight')
hold on; plot(Waircraft,y,'o'); hold off
Wb=WbperW*Waircraft
string2=['Estimated battery weight is ',num2str(Wb),' pounds.']
text2(.25,.15,['           ',string2])
string2=['Payload weight is ',num2str(Wpayload),' pounds.']
text2(.25,.1,['           ',string2])