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

function MASTER441(paramfile, pranformfile, pranformConstraint)
diary
% datafile=input('Input data file: \nNote: You must delete the header
% before using \n:','s');
%DIRECTORY
%LOCATIONS--------------------------------------------------------------]
% rootLOCATION='\\ad.uiuc.edu\ae\correa2\Desktop\AE441\Matlab code\';
% airfoilLOCATION='\\ad.uiuc.edu\ae\correa2\Desktop\AE441\Airfoils';
% pranformLOCATION='\\ad.uiuc.edu\ae\correa2\Desktop\AE441\Planforms';
% pranformConstraintLOCATION='\\ad.uiuc.edu\engr\ews\homes\desktop\AE441\Pl
% anform Constraints';
rootLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Matlab code';
airfoilLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Airfoils';
pranformLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Planforms';
pranformConstraintLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Planform Constraints';
paramLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Parameter Files';
%-----------------------------------------------------------------------]

n=50;
cd(paramLOCATION);
[variablenames, variablevalues]=textread(paramfile, '%s %f');
cd(pranformLOCATION)
[spans, chords, offsets, airfoil]=textread(pranformfile, '%f %f %f %s');
cd(pranformConstraintLOCATION)
[spansCONST, chordsCONST, offsetsCONST, heightCONST]=textread(pranformConstraint, '%f %f %f %f');
cd(airfoilLOCATION)
[sref swet]=WhettedArea(spans, chords, offsets, airfoil, n);
cd(rootLOCATION)

h1=13.33;
h2=10.5;
h3=9.5;

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

wingspan=2*spans(end);
AR=wingspan^2/sref;
%CRUISE------------------------------------------------------------------]
cd0_cruise=findcd0(c_f, M_cruise, sref, swet);
% finddragpolar(cd0_cruise, cLalpha_cruise, cL0, efficiency, AR, cd_correction, 'cruise')
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;


hold on
cd(airfoilLOCATION)

[planformLE planformTE]=plotpranform(spans, chords, offsets, n);
[constraintLE constraintTE]=plotconstraint(spansCONST, chordsCONST, offsetsCONST, n);
plotheight=plotThick(spans, chords, offsets, airfoil, h1, n);
plotheight2=plotThick(spans, chords, offsets, airfoil, h2, n);
plotheight3=plotThick(spans, chords, offsets, airfoil, h3, n);

planformLE=[flipud(planformLE); -planformLE(:,1) planformLE(:,2)];
planformTE=[flipud(planformTE); -planformTE(:,1) planformTE(:,2)];
constraintLE=[flipud(constraintLE); -constraintLE(:,1) constraintLE(:,2)];
constraintTE=[flipud(constraintTE); -constraintTE(:,1) constraintTE(:,2)];
plotheight=[flipud(plotheight); -plotheight(:,1) plotheight(:,2)];
plotheight2=[flipud(plotheight2); -plotheight2(:,1) plotheight2(:,2)];
plotheight3=[flipud(plotheight3); -plotheight3(:,1) plotheight3(:,2)];


plot(planformLE (:,1), planformLE(:,2))
plot(planformTE (:,1), planformTE(:,2))
plot(constraintLE (:,1), constraintLE(:,2))
plot(constraintTE (:,1), constraintTE(:,2))
plot(plotheight(:,1),plotheight(:,2),'r')
plot(plotheight2(:,1),plotheight2(:,2),'r:')
plot(plotheight3(:,1),plotheight3(:,2),'r-.')
axis equal
hold off
fprintf('\n----------Cruise----------\n')
fprintf('alfa: %f, Cd0: %f, CL: %f, CD: %f, CM: %f, L/D: %f \n', alphareq_cruise, cd0_cruise, cLreq_cruise, cD_cruise, cM_cruise, eloverdee_cruise)
cd(rootLOCATION)
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, e, 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;
        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
function [pranformArea WetArea]=WhettedArea(spans, chords, ~, airfoil, n)
% fileName = 'AirData.mat';
wingSpan = 1;
% n = 50;

engineArea=600;
Area = 0; 
pranformArea = 0;

Chords = findChords(n, chords, spans);
for i = 1:(size(chords)-1)
    for j = 1:n
        Chords(i,j);
        dist = wingSpan*(spans(i+1)-spans(i))./n;
        Area = Area+findPerim( airfoil{i}, Chords(i,j), wingSpan).*dist;
    end
end

for i = 1:(size(chords)-1)
    for j = 1:n
        dist = wingSpan*(spans(i+1)-spans(i))./n;
        pranformArea = pranformArea+Chords(i,j).*dist;
%         fprintf('Chord: %f  pranformArea: %f \n', Chords(i,j), pranformArea)
    end
end

Area = Area *2+engineArea;

WetArea = Area;
pranformArea = pranformArea*2;

fprintf('Planform Area: %f, Wetted Area: %f \n', pranformArea, WetArea)
end
function Chords = findChords( n, Chord, X )
for( i = 1:(length(Chord)-1) )
    Xtemp = linspace(X(i),X(i+1),n);
    Chords(i,:) = (Chord(i+1)-Chord(i))./(X(i+1)-X(i)).*(Xtemp-X(i))+Chord(i);
    
end
end
function peri = findPerim(airFile, chord, span)
peri = 0;
Airfoil = load(airFile);
Airfoil = Airfoil*chord*span;
for i = 1:(size(Airfoil,1)-1)
    peri = peri + sqrt((Airfoil(i+1,1)-Airfoil(i,1))^2+(Airfoil(i+1,2)-Airfoil(i,2))^2);
end
end
function [plotLE plotTE]=plotpranform(spans, chords, offsets, n)
Chords = findChords(n, chords, spans);

%plot planform
for i = 1:(size(chords)-1)
    dist=0;
    for j = 1:n
        Chords(i,j);
        dist = (spans(i+1)-spans(i))./n;
        
        if i==1
            xplotLE(i,j+1)=spans(1)+dist*(j);
        else
            xplotLE(i,1)=xplotLE(i-1,end);
            xplotLE(i,j+1)=xplotLE(i-1,end)+dist*(j);
        end
        
        fitX(1)=spans(i);
        fitX(2)=spans(i+1);
        fitY(1)=offsets(i);
        fitY(2)=offsets(i+1);
        a = polyfit(fitX, fitY,1);
        yplotLE(i,j) = a(1)*xplotLE(i,j)+a(2);
        
        
        
    end
end


plotLE=zeros(1,2);
plotTE=zeros(1,2);

for i=1:(size(chords)-1)
    for j=1:n
        plotLE=[plotLE; xplotLE(i,j) -yplotLE(i,j)];
        plotTE=[plotTE; xplotLE(i,j) -(yplotLE(i,j)+Chords(i,j))];
    end
end
plotTE(1,:)=[];
plotLE(1,:)=[];
end
function [plotLE plotTE]=plotconstraint(spans, chords, offsets, n)

%plot constraint
for i = 1:(size(chords)-1)
    for j = 1:n
        dist = (spans(i+1)-spans(i))./n;
        if i==1
            xplotLE(i,j+1)=spans(1)+dist*(j);
        else
        xplotLE(i,1)=xplotLE(i-1,end);    
        xplotLE(i,j+1)=xplotLE(i-1,end)+dist*(j);
        end
        yplotLE(i,j) = offsets(i+1);
    end
end

plotLE=zeros(1,2);
plotTE=zeros(1,2);

for i=1:(size(chords)-1)
    for j=1:n
        plotLE=[plotLE; xplotLE(i,j) -yplotLE(i,j)];
        plotTE=[plotTE; xplotLE(i,j) (-yplotLE(i,j)-chords(i+1))];
    end
end

plotTE(1,:)=[];
plotLE(1,:)=[];
end
function plotThick=plotThick(spans, chords, offsets, airfoil, h, n)
n=n/25;
Chords = findChords(n, chords, spans);
plotThickLE=zeros(1,2);
plotThickTE=zeros(1,2);

for i = 1:(size(chords)-1)
    for j = 1:n
        Chords(i,j);
        dist = (spans(i+1)-spans(i))./n;
        
        if i==1
            xplotLE(i,j+1)=spans(1)+dist*(j);
        else
            xplotLE(i,1)=xplotLE(i-1,end);
            xplotLE(i,j+1)=xplotLE(i-1,end)+dist*(j);
        end
        
        fitX(1)=spans(i);
        fitX(2)=spans(i+1);
        fitY(1)=offsets(i);
        fitY(2)=offsets(i+1);
        a = polyfit(fitX, fitY,1);
        yplotLE(i,j) = a(1)*xplotLE(i,j)+a(2);
        
        chordloc=linspace(0, Chords(i,j));
        
        for k=1:length(chordloc)-1
            thickness(k)=findthickness(airfoil{i}, Chords(i,j), chordloc(k));
            if thickness(k)>=h
                plotThickLE=[plotThickLE; xplotLE(i,j) -(chordloc(k)+yplotLE(i,j))];
                break
            end
        end
        
        for k=1:length(chordloc)-1
            thicknessLE(k)=findthickness(airfoil{i}, chords(i), chordloc(length(chordloc)-k));
            if thickness(k)>=h
                plotThickTE=[plotThickTE; xplotLE(i,j) -(chordloc(length(chordloc)-k)+yplotLE(i,j))];
                break
            end
        end
        
    end
end


plotLE=zeros(1,2);
plotTE=zeros(1,2);

for i=1:(size(chords)-1)
    for j=1:n
        plotLE=[plotLE; xplotLE(i,j) -yplotLE(i,j)];
        plotTE=[plotTE; xplotLE(i,j) -(yplotLE(i,j)+Chords(i,j))];
    end
end

plotTE(1,:)=[];
plotLE(1,:)=[];
plotThickLE(1,:)=[];
plotThickTE(1,:)=[];
plotThick=[plotThickLE; flipud(plotThickTE)];
end
function thickness=findthickness(airfoil, Chord, xLoc)
Airfoil = load(airfoil);
Airfoil = Airfoil*Chord;
% blah....goodnight
    
    AirfoilX=Airfoil(:,1);
    [min_diff, k]=min(abs(AirfoilX));
    
    for j=1:k
    AirfoilYupper(j)=Airfoil(j,2);
    AirfoilXupper(j)=Airfoil(j,1);
    end
    
    for j=1:length(Airfoil)-k
    AirfoilYlower(j)=Airfoil(length(Airfoil)-(j-1), 2);
    AirfoilXlower(j)=Airfoil(length(Airfoil)-(j-1), 1);
    end
    
    AirfoilYupper=fliplr(AirfoilYupper);
    AirfoilXupper=fliplr(AirfoilXupper);
    AirfoilYlower=fliplr(AirfoilYlower);
    AirfoilXlower=fliplr(AirfoilXlower);
    AirfoilYlower=[0 AirfoilYlower];
    AirfoilXlower=[0 AirfoilXlower];
    
    [min_diff, k]=min(abs(AirfoilX));
    
    [min_diff, kU]=min(abs(AirfoilXupper-xLoc));
    [min_diff, kL]=min(abs(AirfoilXlower-xLoc));
    
    if xLoc<AirfoilXupper(kU)
        xFitU(1)=AirfoilXupper(kU-1);
        xFitU(2)=AirfoilXupper(kU);
        yFitU(1)=AirfoilYupper(kU-1);
        yFitU(2)=AirfoilYupper(kU);
    else
        xFitU(1)=AirfoilXupper(kU);
        xFitU(2)=AirfoilXupper(kU+1);
        yFitU(1)=AirfoilYupper(kU);
        yFitU(2)=AirfoilYupper(kU+1);
    end
      
    if xLoc<AirfoilXlower(kL)
        xFitL(1)=AirfoilXlower(kL-1);
        xFitL(2)=AirfoilXlower(kL);
        yFitL(1)=AirfoilYlower(kL-1);
        yFitL(2)=AirfoilYlower(kL);
    else
        xFitL(1)=AirfoilXlower(kL);
        xFitL(2)=AirfoilXlower(kL+1);
        yFitL(1)=AirfoilYlower(kL);
        yFitL(2)=AirfoilYlower(kL+1);
    end
            aU = polyfit(xFitU, yFitU,1);
            aL = polyfit(xFitL, yFitL,1);
%             yplotLE(i,j) = a(1)*xplotLE(i,j)+a(2);
%     thickness =  AirfoilYupper(xLoc_indu)-AirfoilYlower(xLoc_indl);
    thickness =  (aU(1)*xLoc+aU(2))-(aL(1)*xLoc+aL(2));
end