% FILE: Constraint3.m
% Script to generate constraint diagram: 
%
%
disp(' '); disp('*** Start here ***'); disp(' ')
% DataSection
WperSmin=.2   % Limit on the axes of the constraint diagram (lbf/ft^2)
WperSmax=.8   % Limit on the axes of the constraint diagram (lbf/ft^2)
WperPmin= 0   % Limit on the axes of the constraint diagram (lbf/hp)
WperPmax=200   % Limit on the axes of the constraint diagram (lbf/hp)
Vs=20 % Stall speed (ft/sec)
CLmax=[1.3,1.4,1.5] % Possible values of maximum lift coefficient (nondimensional), use 3 of them
rho=0.002377 % air density (slugs/ft^3)
Vcr=28       % cruise speed (ft/sec)
EtaP=.5 % propeller efficiency (nondimensional)
CD0=[.022,.026,.030] % Possible values of CD0, use 3 of them
LoverDmax=[10,12] % estimated maximum lift to drag ratio
gamma=20/57.3 % Take-off flight path angle


% Stall speed constraint
WperS1=.5*rho*Vs^2*CLmax; % wing loading constraints
numCLmax=length(CLmax);
ifig=0;
ifig=ifig+1; figure(ifig)
clf
axis([WperSmin,WperSmax,WperPmin,WperPmax])
WperSdat=[WperS1(1),WperS1(1)];
WperPdat=[WperPmin,WperPmax];
plot(WperSdat,WperPdat)
hold on
title('Constraint Diagram')
hash_right(WperSdat,WperPdat)
%hash_left(WperSdat,WperPdat,30)
WperSdat=[WperS1(2),WperS1(2)];
plot(WperSdat,WperPdat)
hash_right(WperSdat,WperPdat)
%hash_left(WperSdat,WperPdat)
WperSdat=[WperS1(3),WperS1(3)];
plot(WperSdat,WperPdat)
hash_right(WperSdat,WperPdat)
%hash_left(WperSdat,WperPdat,-30)
%string1=['Stall: CLmax=[', num2str(CLmax),'], Vs= ', num2str(Vs), ' f/s'];
%text2(.2,.2,string1)


% Cruise speed constraint
slopes=((.75)*(550)/(.5*rho*1.1))*EtaP./(Vcr^3*CD0);
inc=(WperSmax-WperSmin)/10;
WperSdat=WperSmin:inc:WperSmax;
plot(WperSdat,slopes(1)*WperSdat)
hash_left(WperSdat,slopes(1)*WperSdat,0)
plot(WperSdat,slopes(2)*WperSdat)
hash_left(WperSdat,slopes(2)*WperSdat,0)
plot(WperSdat,slopes(3)*WperSdat)
hash_left(WperSdat,slopes(3)*WperSdat,0)
%string2=['Cruise Constraint: CD0=[', num2str(CD0),'], Vcruise= ', num2str(Vcr), ' ft/sec'];
%text2(.05,.9,string2)
%text2(.3,.03,'DESIGN SPACE')

% Climb constraint
WperPclimb=550*EtaP./(Vcr./(.866*LoverDmax)+Vcr*sin(gamma))
plot([WperSmin WperSmax],[WperPclimb(1) WperPclimb(1)])
plot([WperSmin WperSmax],[WperPclimb(2) WperPclimb(2)])
hash_left([WperSmin WperSmax],[WperPclimb(1) WperPclimb(1)])
hash_left([WperSmin WperSmax],[WperPclimb(2) WperPclimb(2)])
%string3=['Climb constraint, gamma= ',num2str(gamma*57.3),' deg, Vclimb= ',num2str(Vcr),' ft/sec. Lower L/D gives lower line.']
%text2(.05,.08,string3)
%string4=['L/D max= ',num2str(LoverDmax)]
%text2(.05,.15,string4)

xlabel('Wing loading lbf/ft^2')
ylabel('Power loding (lbf/hp)')
axis([WperSmin,WperSmax,WperPmin,WperPmax])
%disp('Click twice on the desired design point')
% [X,Y] = GINPUT(N)
%[WperSin,WperHPin]=ginput(1)
%plot(WperSin,WperHPin,'rx')
%weight=2.5
%S=weight/WperSin
%Bhp=weight/WperHPin

% Add text to figur4es after drawing last plot because scaling changes.
string1=['Stall Constraint: CLmax=[', num2str(CLmax),'], Vs= ', num2str(Vs), ' f/s'];
text2(.4,.4,string1)

string2=['Cruise Constraint: CD0=[', num2str(CD0),'], Vcruise= ', num2str(Vcr), ' ft/sec'];
text2(.05,.9,string2)
text2(.3,.03,'DESIGN SPACE')

string3=['Climb constraint, gamma= ',num2str(gamma*57.3),' deg, Vclimb= ',num2str(Vcr),' ft/sec. Lower L/D gives lower line.']
text2(.05,.08,string3)
string4=['L/D max= ',num2str(LoverDmax)]
text2(.05,.15,string4)

hold off
