clear
clear functions
close all
ifig=0;
hold off
disp(' '); disp('>>> Start of script. <<<'); disp(' ')
echo on
% Script to design an end-to-end propulsion system for an
% electric-powered propeller-driven aircraft.
% Given:
%   drag polar, 
%   aircraft weight, air density, 
%   pitch to diameter ratio of the prop and prop data,
%   motor constants for a particular motor.
%
% Find:
%   speed for maximum endurance,
%   propeller diameter,
%   gear ratio,
%   voltage at which to operate the motor,
%   battery sizes to acheive the desited battery voltage,
%   endurance for single strand and dual strand batteries,
%   for an aircraft flying straight and level or turning with a 
%   specified turn radius.
%
% PHASE 1: AIRCRAFT SUBSYSTEM
%
% This data is roughly for the Boiler Express Aircraft
CD0=.0240; % drag coefficient when CL=0.
A=9; % aspect ratio span squared divided by reference area
e=.75; % Oswalds efficiency factor
V=10:1:40; % velocity in ft/sec
rho=.002377; % air density in slugs/ft^3
W= 10.25;    % aircraft weight in pounds (lbf)
S=   13.5;   % wing area
R=80;        % Radius of steady turn (ft)
echo off
[Pop,Vop,EtaAircraft,Pre,Ve,ifig]=DesignAircraft2(W,rho,S,CD0,A,e,V,ifig,R);

%
% PHASE 2: PROPELLER SUBSYSTEM
%
% Define hypothetical propeller for specified p/D=.631
string1=['Data for hypothetical propeler with tau=p/D=.631'];
disp(string1);
echo on
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];
echo off
[Jstar,CTstar,CPstar,Etastar,ifig]=PlotProp(J,CT,CP,string1,ifig);
% Propeller power equals the operating power of the aircraft
% The propeller is operating at a forward speed given by Vop
% Find the required propeller diameter
[Dactual,Jactual,nactual,RPMactual,EtaPropactual,Poutactual,Ppinactual,ifig]=...
DesignProp(rho,Pop,Vop,Jstar,CTstar,CPstar,Etastar,J,CT,CP,ifig); % <<< lots of work here
RPMprop=RPMactual;
PinPropWatts=Ppinactual*1.356; % convert to watts from ft-lbf/sec
PlotProp2(Jactual,EtaPropactual,J,CT,CP,ifig);
Pitch=tau*Dactual*12;
string2=['  You have selected a ',num2str(Dactual*12),'x',num2str(Pitch),' prop.'];
disp(string2); disp(' ');


%
% PHASE 3: GEAR BOX AND MOTOR SUBSYSTEMS
%
EtaGear=.96; % Efficiency of the gearbox
% At this point in the code the gear ratio is unknown.
% Required motor output taking into account the inefficiency of the grarbox
PoutMotorWatts=PinPropWatts/EtaGear; 


% 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

%
%
[Vinstar,Iinstar,Pinwattsstar,RPMstar,PoutHP,EtaMotorMax,ifig]=MotorMaxEff(PoutMotorWatts,Kv,Kt,R,Io,ifig); % <<<<<
disp(' '); disp('PRELIMINARY MOTOR AND GEARBOX ANALYSIS');
string20=['  The output power of this motor has to be ',num2str(PoutMotorWatts),' watts or ',num2str(PoutHP),' Hp.'];
string21=['  For maximum motor efficiency this motor must be provided with ',num2str(Vinstar),' input volts.'];
string22=['  Under these conditions, the input current will be ',num2str(Iinstar),' amperes'];
string23=['  and the input power will be ',num2str(Pinwattsstar),' watts.'];
string24=['  The motor shaft will be spinning at ',num2str(RPMstar),' RPM.'];
string25=['  The motor efficiency will be ',num2str(EtaMotorMax),'.'];
disp(string20); disp(string21); disp(string22); disp(string23); disp(string24); disp(string25); 
m=RPMstar/RPMprop;
string26=['  To match this motor to the propeller requires an optimal gear box ratio (m*) of ',num2str(m),'.'];
disp(string26); disp(' ');
mactual=input('What gear box ratio (m) do you want to use? Default=m*  >>>>');
disp(' ');
if isempty(mactual)
	mactual=m;                    % use default gear ratio for running the motor at maximum efficiency
	string20d=['You have selected the default gear ratio ',num2str(mactual),'.'];
	disp(string20d);
	RPMactual=mactual*RPMprop;
else
	RPMactual=mactual*RPMprop;    % Compute for specified gear ratio
	string20d=['You have selected the gear ratio ',num2str(mactual),'.'];
	disp(string20d);disp(' ');
	string20b=['At this point the motor RPM and output power of the motor are specified, so motor inputs can be found.'];
	string20c=['    Motor RPM= ',num2str(RPMactual),' RPM and motor output power= ',num2str(PoutMotorWatts),' watts'];
	disp(string20b); disp(string20c);
	% compute motor input properties
end
disp(' ');

% [Vin,Iin,Pinwatts,PoutHP,EtaMotor,ifig]=MotorInputs(PoutWatts,RPM,Kv,Kt,R,Io,ifig)
[Vinactual,Iinactual,Pinwattsactual,PoutHP,EtaMotoractual,ifig]=MotorInputs(PoutMotorWatts,RPMactual,Kv,Kt,R,Io,ifig);

string30a=['MOTOR AND GEAR BOX DESIGN SUMMARY'];
string30= ['  The output power of this motor is ',num2str(PoutMotorWatts),' watts or ',num2str(PoutHP),' Hp.'];
string31= ['  The motor input voltage is        ',num2str(Vinactual),' volts.'];
string32= ['  The motor input current is        ',num2str(Iinactual),' amperes.'];
string33= ['  and the motor input power is      ',num2str(Pinwattsactual),' watts.'];
string34= ['  The electric motor shaft is spinning at  ',num2str(RPMactual),' RPM.'];
string34a=['  The gear box output shaft is spinning at ',num2str(RPMprop),' RPM.'];
string35= ['  The gear box ratio (m) is         ',num2str(mactual),'.'];
string36= ['  The motor efficiency is           ',num2str(EtaMotoractual),'.'];
string37= ['  The gear box efficiency is        ',num2str(EtaGear),'.'];
disp(string30a); disp(string30); disp(string31); disp(string32); disp(string33); disp(string34); disp(string34a);
disp(string35); disp(string36); disp(string37)

% Comments on these motor operating conditions
string40e=['COMMENTS ON THE MOTOR OPERATING CONDITIONS'];
string40= ['  You are asking a motor with nominal horsepower of ',num2str(HPnom),' to produce a horsepower of ',num2str(PoutHP),'.'];
string40a=['  Since the produced power is greater than the nominal power you are overdriving the motor.'];
string40b=['  Consider a larger motor!'];
string40c=['  Since the produced power is less than the nominal power you are underdriving the motor.'];
string40d=['  Consider a smaller motor!'];
disp(' '); disp(string40e); disp(string40);
if PoutHP>HPnom; 
	disp(string40a); disp(string40b);
else 
	disp(string40c); disp(string40d);
end

% Plot the motor operating properties for the given voltage.
Iintest=Io:30;
[Pouttest,Poutfppstest,PoutHptest,Pintest,Touttest,Toutfptest,RPMtest,etatest]=motor(Vinactual,Iintest,Kv,Kt,R,Io);
[ifig]=MotorPlot2(Vinactual,Iinactual,Iintest,Pouttest,Touttest,RPMtest,etatest,titstr,ifig);


%
% PHASE 4: BATTERY SUBSYSTEM
%
disp(' '); disp('BATTERY SUBSYSTEM'); disp('  The battery pack will be made up of individual cells with the following properties:')
VoltsPerCell=1.2;         % volts per cell
mAmpsHoursPerCell=3000;   % milliamps hours per cell
gramspercell=57.24;   %grams per cell
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.')
nCells1=ceil(Vinactual/VoltsPerCell);
nVolts=nCells1*VoltsPerCell;
weight1=nCells1*lbfpercell;
BatteryEnergy1Joule=(mAmpsHoursPerCell*3600/1000)*nVolts; % Joules=watt*sec=ampere*volt*sec
ActualEnduranceMin1=(1/60)*(BatteryEnergy1Joule/Pinwattsactual); % predicted endurance for single strand battery, minutes
string50=['  ',num2str(nCells1),' total cells, arranged in a 1x',num2str(nCells1),' array.'];
string51=['  producing ',num2str(nVolts),' volts, and giving a predicted endurance of ',num2str(ActualEnduranceMin1),' minutes,'];
string54=['  and weighing ',num2str(weight1),' lbf.'];
disp(string50); disp(string51); 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
ActualEnduranceMin2=2*ActualEnduranceMin1; % predicted endurance for dual strand battery, minutes
string52=['  ',num2str(nCells2),' total cells, arranged in a 2x',num2str(nCells1),' array.'];
string53=['  producing ',num2str(nVolts),' volts, and giving a predicted endurance of ',num2str(ActualEnduranceMin2),' minutes'];
string55=['  and weighing ',num2str(weight2),' lbf.'];
disp(string52); disp(string53); disp(string55);
disp(' '); disp('  These endurance numbers do not include energy spent in other mission phases, (e.g., take off, climb, turning).')


%
% SUMMARY 
%
disp(' '); disp('SYSTEM SUMMARY');
EtaOverall=EtaMotoractual*EtaGear*EtaPropactual*EtaAircraft;
string60=['  The overall efficiency of your design in the product of the individual subsystem efficiencies'];
string61=['  The overall efficiency is ',num2str(EtaOverall),'.'];
disp(string60); disp(string61);





