clear
clear functions
close all
ifig=0;
hold off
disp(' '); disp('>>> Start of script. <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'); disp(' ')
echo on
% Script to analyze an end-to-end propulsion system for an
% electric-powered propeller-driven aircraft.
% Dominick Andrisani
% September 20, 2005
%
% Given:
%   drag polar, 
%   aircraft weight, air density, 
%   pitch to diameter ratio of the prop and prop data,
%   motor constants for a particular motor.
%   propeller diameter,
%   gear ratio,
%   battery properties (number of cells in series, cell properties)

%
% Find:
%   Power available from the motor and power required as a function of aircraft velocity,
%   Speed of steady flight for an aircraft flying straight and level or turning with a 
%       specified turn radius,
%   Energy consumed for turning flight for a specified number of minutes and for climb.
%
% INPUTS: make changes here for different configurations and aircraft ****
nCells1=10; % Number of battery cells in series.
VoltsPerCell=1.2;         % volts per cell
mAmpsHoursPerCell=3000;   % milliamps hours per cell
gramspercell=57.24;   %grams per cell
PercentThrottle=100; % Percent throttle
EtaSpeed=1.0; % Efficiency of the gearbox
ii=8; % selected motor from motors listed in motordata
EtaGear=.96; % Efficiency of the gearbox
m=2;   % Gear ratio
Din=36;  % Propeller diameter inches
rho=.002377; % air density in slugs/ft^3
% Define hypothetical propeller for specified p/D=.631
string1q=['Data for hypothetical propeler with tau=p/D=.631'];
tau=.631;
J = [     0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000 0.7000    0.8000];
CT =[0.0933    0.0867    0.0773    0.0667    0.0553    0.0427    0.0293 0.0160    0.0013];
CP =[0.0400    0.0393    0.0387    0.0367    0.0333    0.0293    0.0227 0.0147    0.0027];
CD0=.0240; % drag coefficient when CL=0.
A=9; % aspect ratio span squared divided by reference area
e=.75; % Oswalds efficiency factor
W= 10.25;    % aircraft weight in pounds (lbf)
S=   13.5;   % wing area
Radius=80;        % Radius of steady turn (ft)
enduranceMin=5; % Endurance in minutes
Vclimb=32; % Vclimb: intenced climb speed, ft/sec
altitude=10; % altitude for cruise/turn 9ft)
%
% PHASE 1: BATTERY SUBSYSTEM
%
echo off
disp(' '); disp('BATTERY SUBSYSTEM'); disp('  The battery pack will be made up of individual cells with the following properties:')
slugspercell=(gramspercell/1000)/14.5939;
lbfpercell=32.17*slugspercell;
string54=['  ',num2str(VoltsPerCell),' volts per cell, and ',num2str(mAmpsHoursPerCell),' milliamp hours per cell, and ',num2str(gramspercell),' grams per cell.'];
disp(string54);
disp(' '); disp('  A single string battery pack designed for the above conditions will have the following properties.')
nVolts=nCells1*VoltsPerCell;
weight1=nCells1*lbfpercell;
BatteryEnergy1Joule=(mAmpsHoursPerCell*3600/1000)*nVolts; % Joules=watt*sec=ampere*volt*sec
string50=['  ',num2str(nCells1),' total cells, arranged in a 1x',num2str(nCells1),' array,'];
string51=['  producing ',num2str(nVolts),' volts,'];
string51a=['  containing ',num2str(BatteryEnergy1Joule),' joules of energy,'];
string54=['  and weighing ',num2str(weight1),' lbf.'];
disp(string50); disp(string51); disp(string51a); disp(string54);


disp(' '); disp('  A dual string battery pack designed for the above conditions will have the following properties.')
nCells2=2*nCells1;                    % Cells in both strands
weight2=nCells2*lbfpercell;
BatteryEnergy2Joule=2*BatteryEnergy1Joule; % Joules=watt*sec=ampere*volt*sec
string52=['  ',num2str(nCells2),' total cells, arranged in a 2x',num2str(nCells1),' array,'];
string53=['  producing ',num2str(nVolts),' volts,'];
string53a=['  containing ',num2str(BatteryEnergy2Joule),' joules of energy,'];
string55=['  and weighing ',num2str(weight2),' lbf.'];
disp(string52); disp(string53);  disp(string53a); disp(string55);

echo on
%
% Phase 2: Speeed controller
%
echo off
Vspeed=nVolts*PercentThrottle/100; % scale battery voltage by percent throttle

echo on
%
% Phase 3: Motor
%
echo off

% Get array of motor data
motordata
[nmotor,col]=size(mdata);
% Select motor from the available choices
% string20a=['Which of the  ',num2str(nmotor),' available motors do you want to use? Default=1  >>>>'];
% disp(' ');
% ii=input(string20a);
% if isempty(ii)
% 	ii=1; % Select motor #1
% end


titstr=mname(ii,:);
string29=['You have selected the motor ',num2str(ii),', the ',titstr,' motor.']; disp(' ');
disp(string29);
Kv=mdata(ii,1);      % RPM/volt
Kt=mdata(ii,2);      % inch-ounce per ampere
R= mdata(ii,3);      % Ohms
Io=mdata(ii,4);      % amperes
Vin=mdata(ii,6) ;    % nominal motor voltage, volts
HPnom=mdata(ii,7);   % nominal horsepower at the conditions for max efficiency

Poutwatts=HPnom*745.7; % watts
ifig=1; ifigMotor=ifig; % Plot of motor data
[Vinstar,Iinstar,Pinwattsstar,RPMstar,PoutHP,EtaMotorMax,ifig]=MotorMaxEff(Poutwatts,Kv,Kt,R,Io,ifig);

Iin=1.5*Io:Iinstar/50:3*Iinstar;
ni=length(Iin);
% the code below takes vector Iin and makes all the outputs vectors of the
% same size
[PoutMotorWatts,Poutfpps,PoutHp,Pin,Tout,Toutfp,RPMmotor,eta]=motor(Vspeed,Iin,Kv,Kt,R,Io); % Vspeed is the battery/speed controller voltage
%[PoutMotorWatts,Poutfpps,PoutHp,Pin,Tout,Toutfp,RPM,eta]=motor(Vinstar,Iin,Kv,Kt,R,Io); % Vinstar is best input voltage
subplot(311); plot(Iin,PoutMotorWatts)
title('Motor Properties')
xlabel('Current (amp)')
ylabel('Output Power (watt)')
subplot(312); plot(Iin,Pin)
xlabel('Current (amp)')
ylabel('Input Power (watt)')
subplot(313); plot(Iin,eta)
xlabel('Current (amp)')
ylabel('Efficiency')


echo on
%
% Phase 4: Gearbox
%
echo off


PinPropWatts=PoutMotorWatts*EtaGear; 
RPMprop=RPMmotor/m;
ifig=ifig+1; figure(ifig) % Plot of gearbox data
subplot(211); plot(Iin,PinPropWatts)
title('Gearbox Properties')
xlabel('Motor Current (amp)')
ylabel('Gearbox Output Power(watts)')
subplot(212); plot(Iin,RPMprop)
xlabel('Motor Current (amp)')
ylabel('Gearbox Output RPM')

echo on
%
% Phase 5: Propeller
%
echo off
disp(string1q);
% Plot of propeller data
[Jstar,CTstar,CPstar,Etastar,ifig]=PlotProp(J,CT,CP,string1q,ifig);
% Propeller input power and RPM is specified for each Iin.
% We can compute the required velocity of the aircraft that will consume
% that much input power. Then we compute the advance ratio and the 
% resulting output power (power available) and thrust.
Dft=Din/12; % feet
% Ppin=rho*n^3Dft^5*Cp
PinPropfpps=PinPropWatts/1.355818; % ft-lbf/sec
n=RPMprop/60; % revolutions per second
CPmat=PinPropfpps./(rho*n.^3*Dft^5);
CPmax=max(CP);
ic=find(CPmat<=CPmax);
if isempty(ic)
    disp(' <<< This prop is not big enough to absord this much power. Try a larger diameter prop.')
    break
else
    Iinpossible=Iin(ic);
    CPpossible=CPmat(ic);
    npossible=n(ic);
    ifig=ifig+1; figure(ifig); plot(Iinpossible,CPpossible)
    title('Possible CP for possible Iin')
    xlabel('Motor Current (amps)')
    ylabel('Propeller CP')
    text2(.05,.9,'The propeller is not capable of absorbing any more power then shown here.')
    text2(.05,.85,'If you need to absord more power (go to higher currents) make the propeller bigger')
end

% YI = INTERP1(X,Y,XI)
Jpossible=interp1(CP,J,CPpossible);
ip=find(isfinite(Jpossible)); % some J will be NaN. Get rid of those here.
Jpossible2=Jpossible(ip);
CPpossible2=CPpossible(ip);
npossible2=npossible(ip);
Vpossible2=Dft*npossible2.*Jpossible2;
CTpossible2=interp1(J,CT,Jpossible2);
Iinpossible2=Iinpossible(ip);
Thrustpossible2=rho*npossible2.^2*Dft.^4.*CTpossible2;
PowerAvailable=Thrustpossible2.*Vpossible2;


ifig=ifig+1; figure(ifig) % Plot of Power available and power required.
plot(Vpossible2,PowerAvailable,'k-')
hold on; plot(Vpossible2,PowerAvailable,'kx'); 
text2(.05,.9,'Black x is Power Available from the propeller')
text2(.05,.85,'Red + is Power Required by the aircraft')
hold off
title('Power Available and Power Required')
xlabel('Velocity (ft/sec)')
ylabel('Power (ft-lbf/sec)')

echo on
%
% Phase 6: Aircraft
%
echo off

tanphi=(Vpossible2.*Vpossible2)/(32.17*Radius); %tangent of bank angle Roskam eqn 4.89
phirad=atan(tanphi); % bank angle, radians
n=1./cos(phirad); % load factor, g, Roskam eqn 4.92
k=1/(pi*A*e); % k from the drag polar CD=CD0+k*CL^2
CL=2*n*W./(rho*S*Vpossible2.*Vpossible2); % Lift cofficient as a function of velocity
CD=CD0+k*CL.*CL; % Drag coefficient as a function of velocity
PowerRequired=.5*rho*S*Vpossible2.^3.*CD; % Power required ft/lbf/sec
hold on; plot(Vpossible2,PowerRequired,'r-'); hold off %% More plot of Power available and power required.
hold on; plot(Vpossible2,PowerRequired,'r+'); hold off

PowerExcess=PowerAvailable-PowerRequired;
% YI = INTERP1(X,Y,XI)

% Take powers that are on the front-side of the power required curve.
PaMax=max(PowerAvailable);
% YI = INTERP1(X,Y,XI)
Vmax=interp1(PowerAvailable,Vpossible2,PaMax);
kk=find(Vpossible2>=Vmax);
Vpossible3=Vpossible2(kk);
PowerAvailable3=PowerAvailable(kk);
PowerRequired3=PowerRequired(kk);
Vsteady=interp1(PowerAvailable3-PowerRequired3,Vpossible3,0);

%Vsteady=interp1(PowerExcess,Vpossible2,0);
Psteady=interp1(Vpossible2,PowerAvailable,Vsteady);
Isteady=interp1(PowerAvailable,Iinpossible2,Psteady);
Jsteady=interp1(Vpossible2,Jpossible2,Vsteady);
PmotorSteady=Isteady*Vspeed; % watts
PbatterySteady=PmotorSteady/EtaSpeed; % watts
Esteady=PbatterySteady*enduranceMin*60; %Energy to fly at this condition for the specified amount ot time.
% Isteady=interp1(PowerRequired,Iinpossible2,Psteady) % This would work as
% well as the line above.
tanphi=(Vsteady.*Vsteady)/(32.17*Radius); %tangent of bank angle Roskam eqn 4.89
phirad=atan(tanphi); % bank angle, radians
n=1./cos(phirad); % load factor, g, Roskam eqn 4.92
CLsteady=2*n*W./(rho*S*Vsteady.*Vsteady); % Lift cofficient as a function of velocity
string1a=['Lift coefficient at cruise/turn is ',num2str(CLsteady,4),'.'];


% AXIS([XMIN XMAX YMIN YMAX])
limits=axis;
hold on; plot([Vsteady Vsteady],[limits(3) Psteady],':'); 
plot([limits(1) Vsteady],[Psteady Psteady],':');
hold off
string1=['This aircraft will fly in steady level flight at V= ',num2str(Vsteady,4),' ft/sec.'];
string2=['At this speed Power Required=Power Available= ',num2str(Psteady,4),' ft-lbf/sec.'];
string2a=['At this condition the motor is consuming ',num2str(Isteady,4),' amps of power,'];
string2b=['and will consume ',num2str(Esteady,6),' joules of energy from battery to fly for ',num2str(enduranceMin),' minutes.'];

text2(.05,.30,string1)
text2(.05,.24,string2)
text2(.05,.18,string2a)
text2(.05,.12,string2b)
text2(.05,.06,string1a)

disp(' '); disp(string1); disp(string2); disp(string2a); disp(string2b); disp(string1a)

% Lets compute some propeller conditions at Vsteady.
% Jpossible2=Jpossible(ip);
% CPpossible2=CPpossible(ip);
% npossible2=npossible(ip);
% Vpossible2=Dft*npossible2.*Jpossible2;
% CTpossible2=interp1(J,CT,Jpossible2);
% Iinpossible2=Iinpossible(ip);
% Thrustpossible2=rho*npossible2.^2*Dft.^4.*CTpossible2;
% PowerAvailable=Thrustpossible2.*Vpossible2;
% YI = INTERP1(X,Y,XI)
CTsteady=interp1(Vpossible2,CTpossible2,Vsteady);
CPsteady=interp1(Vpossible2,CPpossible2,Vsteady);
Etapossible2=CTpossible2.*Jpossible2./CPpossible2;
etaPropSteady=interp1(Vpossible2,Etapossible2,Vsteady);
% % Plot of propeller data with cruise/turn and climb conditions highlighted
% [Jstar,CTstar,CPstar,Etastar,ifig]=PlotProp(J,CT,CP,string1q,ifig);
% subplot(311); hold on; 
% limits=axis;
% plot([Jsteady Jsteady],[limits(3) CTsteady],'k:')
% plot([limits(1) Jsteady],[CTsteady CTsteady],'k:')
% string1b=['CT for cruise/turn = ',num2str(CTsteady,3),' for J= ',num2str(Jsteady,3)];
% text(.1,CTsteady+.006,string1b)
% plot(Jpossible2,CTpossible2,'rx')
% hold off
% subplot(312); hold on; 
% limits=axis;
% plot([Jsteady Jsteady],[limits(3) CPsteady],'k:')
% plot([limits(1) Jsteady],[CPsteady CPsteady],'k:')
% string1c=['CP for cruise/turn = ',num2str(CPsteady,3),' for J= ',num2str(Jsteady,3)];
% text(.1,CPsteady+.005,string1c)
% plot(Jpossible2,CPpossible2,'rx')
% hold off
% subplot(313); hold on; 
% limits=axis;
% plot([Jsteady Jsteady],[limits(3) etaPropSteady],'k:')
% plot([limits(1) Jsteady],[etaPropSteady etaPropSteady],'k:')
% string1d=['Eta for cruise/turn = ',num2str(etaPropSteady,3),' for J= ',num2str(Jsteady,3)];
% text(.1,Etastar-.2,string1d)
% plot(Jpossible2,Etapossible2,'rx')
% 
% hold off



% ifig=ifig+1; % Plot of turn conditions
% figure(ifig); 
% subplot(311); plot(Vpossible2,phirad*180/pi)
% ylabel('Bank angle, phi (deg)')
% string40=['Steady turn with turn radius of ',num2str(Radius),' ft'];
% title(string40)
% subplot(312); plot(Vpossible2,n)
% ylabel('Load factor n (g)')
% subplot(313); plot(Vpossible2,CL)
% ylabel('Lift coefficient')
% xlabel('Velocity (ft/sec)')

ifig=ifig+1; figure(ifig) % Plot of power available and power required.
plot(Iinpossible2,PowerAvailable,'k-')
hold on; plot(Iinpossible2,PowerAvailable,'kx'); 
text2(.05,.9,'Black x is Power Available from the propeller')
text2(.05,.85,'Red + is Power Required by the aircraft')
hold off
title('Power Available and Power Required')
xlabel('Motor current (amps)')
ylabel('Power (ft-lbf/sec)')
hold on; plot(Iinpossible2,PowerRequired,'r-'); hold off
hold on; plot(Iinpossible2,PowerRequired,'r+'); hold off

echo on
%
% Phase 7: Climb Performance
%
echo off

% Climb is done with wings level so there is no turn to consider
% As a result there is a defferent power required.
k=1/(pi*A*e); % k from the drag polar CD=CD0+k*CL^2
CL=2*W./(rho*S*Vpossible2.*Vpossible2); % Lift cofficient as a function of velocity
CD=CD0+k*CL.*CL; % Drag coefficient as a function of velocity
PowerRequiredLevel=.5*rho*S*Vpossible2.^3.*CD; % Power required ft/lbf/sec
PowerExcessLevel=PowerAvailable-PowerRequiredLevel;


gammaRad=asin(PowerExcessLevel./(W*Vpossible2));
gammaDeg=gammaRad*180/pi;

fig=ifig+1; figure(ifig) % Plot of climb data with climb at specified speed highlighted.
plot(Vpossible2,gammaDeg,'k-')
hold on; plot(Vpossible2,gammaDeg,'kx'); 
%text2(.05,.9,'Black x is Power Available from the propeller')
%text2(.05,.85,'Red + is Power Required by the aircraft')
hold off
title('Climb Performance with Wings Level')
xlabel('Velocity (ft/sec)')
ylabel('Flight Path Angle (deg)')
% YI = INTERP1(X,Y,XI)
gammaCLIMBdeg=interp1(Vpossible2,gammaDeg,Vclimb);


limits=axis;
hold on; plot([Vclimb Vclimb],[limits(3) gammaCLIMBdeg],':')
plot([limits(1) Vclimb],[gammaCLIMBdeg gammaCLIMBdeg],':'); hold off
Pclimb=interp1(Vpossible2,PowerAvailable,Vclimb);
Iclimb=interp1(PowerAvailable,Iinpossible2,Pclimb);
PclimbMotor=Vspeed*Iclimb; % watts
PclimbBattery=PclimbMotor/EtaSpeed;
deltaT=altitude/(Vclimb*sin(gammaCLIMBdeg*pi/180)); % Time to climb specified altitude.
EnergyClimb=PclimbBattery*deltaT; % joules to climb to specified altitude
CLclimb=2*W./(rho*S*Vclimb.*Vclimb); % Lift cofficient as a function of velocity
string7a=['Lift coefficient at climb is ',num2str(CLclimb,4),'.'];
Jclimb=interp1(Vpossible2,Jpossible2,Vclimb);
% Lets compute some propeller conditions at Vclimb.
% Jpossible2=Jpossible(ip);
% CPpossible2=CPpossible(ip);
% npossible2=npossible(ip);
% Vpossible2=Dft*npossible2.*Jpossible2;
% CTpossible2=interp1(J,CT,Jpossible2);
% Iinpossible2=Iinpossible(ip);
% Thrustpossible2=rho*npossible2.^2*Dft.^4.*CTpossible2;
% PowerAvailable=Thrustpossible2.*Vpossible2;
% YI = INTERP1(X,Y,XI)
CTclimb=interp1(Vpossible2,CTpossible2,Vclimb);
CPclimb=interp1(Vpossible2,CPpossible2,Vclimb);
etaPropClimb=interp1(Vpossible2,Etapossible2,Vclimb);




string3=['This aircraft has a climb angle of ',num2str(gammaCLIMBdeg,4),' deg at a speed of ',num2str(Vclimb,4),' ft/sec.'];
string4=['At this condition the propeller is producing power equaling ',num2str(Pclimb,4),' ft-lbf/sec, '];
string5=['and the motor is consuming ',num2str(Iclimb,4),' amps of current at ',num2str(Vspeed,4),' volts.'];
string6=['The motor is consuming ',num2str(PclimbMotor,4),' watts of power and will take ',num2str(deltaT,4),' seconds to climb to ',num2str(altitude),' ft.'];
string7=['Climb will consume ',num2str(EnergyClimb,4),' joules of energy from the battery.'];

text2(.02,.22,string3)
text2(.02,.18,string4)
text2(.02,.14,string5)
text2(.02,.10,string6)
text2(.02,.06,string7)
text2(.02,.02,string7a)
disp(' '); disp(string3); disp(string4); disp(string5);  disp(string6);  disp(string7); disp(string7a);








ifig=ifig+1; figure(ifig) % Plot of motor porperties with two flight conditions highlighted.
% use function AddDottedLine(x,y,string,stringXfraction,stringYfraction)
subplot(311); plot(Iin,PoutMotorWatts)
title('Motor Properties')
ylabel('Output Power (watt)')

subplot(312); plot(Iin,Pin)
title('Motor Properties')
ylabel('Input Power (watt)')
string8=['Aircraft climbing at ',num2str(Vclimb),' ft/sec will consume ',num2str(Iclimb,4),' amps of current.'];
AddDottedLine(Iclimb,PclimbMotor,string8,.05,.9)
string8a=['Aircraft turning/cruising at ',num2str(Vsteady,4),' ft/sec will consume ',num2str(Isteady,4),' amps of current.'];
AddDottedLine(Isteady,PmotorSteady,string8a,.05,.75)
hold on; text2(.05,.6,string8); hold off
disp(' '); disp(string8); disp(string8a); 
hold off

subplot(313); plot(Iin,eta)
xlabel('Current (amp)')
ylabel('Efficiency')
limits=axis; hold on
% YI = INTERP1(X,Y,XI)
etaClimb=interp1(Iin,eta,Iclimb);
string9=['Motor efficiency at climb is ',num2str(etaClimb,2),'.'];
AddDottedLine(Iclimb,etaClimb,[' '],.05,.75)
% YI = INTERP1(X,Y,XI)
etaSteady=interp1(Iin,eta,Isteady);
string9a=['Aircraft turning/cruising at ',num2str(Vsteady,4),' ft/sec will consume ',num2str(Isteady,4),' amps of current.'];
string9b=['Motor efficiency at cruise/turn is ',num2str(etaSteady,2),'.'];
AddDottedLine(Isteady,etaSteady,[' '],.05,.75)
limits=axis; hold on
text2(.05,.44,string8)
text2(.05,.32,string9)
text2(.05,.20,string9a)
text2(.05,.08,string9b)
disp(' '); disp(string8); disp(string9); disp(string9a); disp(string9b);
hold off



% Plot of propeller data with cruise/turn and climb conditions highlighted
% use function AddDottedLine(x,y,string,stringXfraction,stringYfraction)
[Jstar,CTstar,CPstar,Etastar,ifig]=PlotProp(J,CT,CP,string1q,ifig);
subplot(311);
string1b=['CT for cruise/turn = ',num2str(CTsteady,3),' for J= ',num2str(Jsteady,3)];
AddDottedLine(Jsteady,CTsteady,string1b,.1,.9)
hold on; plot(Jpossible2,CTpossible2,'rx')
string10b=['CT for climb = ',num2str(CTclimb,4),' for J = ',num2str(Jclimb,4)];
AddDottedLine(Jclimb,CTclimb,string10b,.1,.75)
hold off % Just to be sure.

subplot(312);
string1c=['CP for cruise/turn = ',num2str(CPsteady,3),' for J= ',num2str(Jsteady,3)];
AddDottedLine(Jsteady,CPsteady,string1c,.1,.2)
hold on; plot(Jpossible2,CPpossible2,'rx')
string10b=['CP for climb = ',num2str(CPclimb,4),' for J = ',num2str(Jclimb,4)];
AddDottedLine(Jclimb,CPclimb,string10b,.1,.1)
hold off % Just to make sure

subplot(313);
string1d=['Eta for cruise/turn = ',num2str(etaPropSteady,3),' for J= ',num2str(Jsteady,3)];
% text(.1,Etastar-.2,string1d)
AddDottedLine(Jsteady,etaPropSteady,string1d,.1,.2)
hold on; plot(Jpossible2,Etapossible2,'rx')
string10a=['Eta for climb = ',num2str(etaPropClimb,2),' for J = ',num2str(Jclimb,4)];
AddDottedLine(Jclimb,etaPropClimb,string10a,.1,.1)
hold off % Just to make sure hold is off.

