
TSPAN=0:.1:30;  % times to run takeoff for
%[row,col]=size(TSPAN)
ETAp=0.5       % propeller efficiency
ETAm=0.6       % Motor efficiency
rho=0.002377   % air density slugs per ft^3
Cd0=0.0300     % paracite drag coefficient
mu=0.05        % coefficient of rolling friction
WperS=0.2:.1:1;   % lbf/ft^2
WperBhp=0.1:.1:1; % lbf/Bhp
Vtakeoff=28       % Take-off velocity, ft/sec
RHOeNiCAD=  72900 % Energy density of NiCad (joule/lbf)
RHOeNiMH=   92500 % Energy density of NiMH (joule/lbf)
RHOeLiPoly=239000 % Energy density of Lithium Polymer (joule/lbf)


for i=1:length(WperS)
    for j=1:length(WperBhp)
        %str1=['wing loading= ',num2str(WperS(i)),', powerloading= ',num2str(WperBhp(j))];
        %disp(str1)
%[T,Y]=ode45(TakeOffEOM,TSPAN,Y0,[],wingloading,powerloading,ETAp,rho,Cd0,mu)
    end
end
wingloading=1
powerloading=10
Y0=[0;.1] % Initial conditions on x and xdot
[row,col]=size(Y0)
%[T,Y]=ode45(@TakeOffEOM,TSPAN,Y0,[],wingloading,powerloading,ETAp,rho,Cd0,mu);
[T,Y]=ode45(@TakeOffEOM,[0 .5],[0 .1],[],wingloading,powerloading,ETAp,rho,Cd0,mu);

% YI = INTERP1(X,Y,XI)
TakeOffTime=interp1(Y(:,2),T,Vtakeoff) % will be NaN if the velocity was ot attained.
if isnan(TakeOffTime); disp('ERROR: You have not integrated long enough to achieve Vtakeoff'); end
TakeOffDistance=interp1(T,Y(:,1),TakeOffTime)

subplot(311)
plot(T,Y(:,1))
ylabel('Position (ft)')
title('Takeoff ground roll')

subplot(312)
plot(T',Y(:,2))
xlabel('Time(sec)')
ylabel('Velocity (ft/sec)')

accelG=(550*ETAp./(Y(:,2)*powerloading)-.5*rho*Y(:,2).*Y(:,2)*Cd0/wingloading-mu);
subplot(313)
plot(T(40:length(T))',accelG(40:length(T)))
xlabel('Time(sec)')
ylabel('acceleration (g)')

% Battery weight fraction for takeoff
WbPerW=745.7*TakeOffTime/(ETAm*powerloading*RHOeNiCAD) % nondimensional

