%Master 441 Code
%all alfas in DEGREES!
%by ascorrea

function master441()
diary
datafile=input('Input data file: \nNote: You must delete the header before using \n:','s');
[variablenames, variablevalues]=textread(datafile, '%s %f');


for k = 1:length(variablenames)
eval([variablenames{k} '= variablevalues(k);']);
end
i=1;
for AR=1:.1:10

%CRUISE------------------------------------------------------------------]
cd0_cruise=findcd0(c_f, M_cruise, sref, swet);
cLreq_cruise=findcLrequired(alt_cruise, sref, M_cruise, weight_cruise);
alphareq_cruise=findalpharequired(cLalpha_cruise, cL0, cLreq_cruise);
cD_cruise=findcD(AR, efficiency, cd0_cruise, cLreq_cruise, cd_correction);
eloverdee_cruise=cLreq_cruise/cD_cruise;
cM_cruise=cMalpha_cruise*alphareq_cruise+cM0;


plotAR(i)=AR;
plotLD(i)=eloverdee_cruise;
i=i+1;
end


plot(plotAR, plotLD)

fprintf('Max L/D: %f at AR: %f\n',max(plotLD), plotAR(max(plotLD)==plotLD));

end

function q=findcd0(c_f, M, sref, swet)
if M>.7
c_f=c_f*(1-.09*M^2);
end
q=c_f*(swet/sref);
end
function q=findcLalpha(AR, M, clalpha, sweep, sex, sref, doverb)
%ascorrea

betasquared=1-M^2;
beta=sqrt(betasquared);
eta=clalpha/(2*pi/beta);
numer=2*pi*AR;
denom1=(AR^2*beta^2)/eta^2;
denom2=(tand(sweep)^2/betasquared);
denom=2+sqrt(4+denom1*(1+denom2));
F=1.07*(1+doverb)^2;

cLalpha=(numer/denom)*(sex/sref)*F;

q=cLalpha;

end
function finddragpolar(cd0, cl_alfa, cL0, sweep_LE, ar, cd_correction, filename)
alphadeg=[-12:2];
% alpha=(pi/180)*(alphadeg);
alpha=alphadeg;
cL_alfa=cl_alfa*(ar/(ar+(2*(ar+4)/(ar+2))));

        for i=1:length(alpha)
        C_L(i)=cL_alfa*(alpha(i))+cL0;
        e=4.61*(1-.045*ar^.68)*cosd(sweep_LE)^.15-3.1;
        cd_i(i)=C_L(i)^2/(pi*ar*e);
        C_d(i)=cd0+cd_i(i)+cd_correction;
        end
        
        graphname='Drag Polar';

hold on
axisX0=linspace(0,max(cd));
axisY0=zeros(1,length(axisX0));
plot(axisX0, axisY0,'k--');
plot(C_d, C_L);
plot(C_d, C_L, '.');
title(graphname);
xlabel('C_d');
ylabel('C_L');
axis([0, .04, -1, 1.5]);
% for i=1:length(alphadeg)
%     alphalabel=num2str(alphadeg(i));
%     label=['  \alpha =' alphalabel];
% text(C_d(i), C_L(i), label)
% end
print ('-dtiff', filename)
close(1)
end
function q=findcLrequired(alt, s, v_mach, lift)
if alt < 36152
temp=59-.00356*alt;
press=2116*((temp+459.7)/518.6)^5.256;

elseif alt < 82345
    temp=-70;
    press=473.1*exp(1.73-.000048*alt);
    
else
    print('Program Not Valid for given altitude')
end

density=press/(1718*(temp+459.7));
v_fps=v_mach*sqrt(1.4*1718*(temp+459.7));
q=.5*density*v_fps^2;

q=lift/(q*s);

end
function q=findalpharequired(cLalpha, cL0, cLrequired)
q=(cLrequired-cL0)/cLalpha;
end
function q=findcD(ar, e, cd0, C_L, cd_correction)

%         e=4.61*(1-.045*ar^.68)*cosd(sweep_LE)^.15-3.1;
        cd_i=C_L^2/(pi*ar*e);
        cd_i=C_L^2/(pi*ar*e);
        q=cd0+cd_i+cd_correction;
end