%Ronald
%all alfas in DEGREES!
%by ascorrea, with help from Jacob Huffman
%"filler" variable just makes it easier to switch between Ronald and Grimace

function Ronald(paramfile, pranformfile, filler)

testrange={'20'};

for testnumber=1:length(testrange)
prefix='FDR0';
name=strcat(prefix, testrange{testnumber});
paramfile=strcat(name, '.param');
pranformfile=strcat(name, '.plan');


%DIRECTORY LOCATIONS--------------------------------------------------------------]
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';
outputsLOCATION='/Users/anthonyscorrea/Documents/MATLAB/AE441 - Home Edition/Ronaldplus Outputs';
%-----------------------------------------------------------------------]

n=50;
cd(paramLOCATION);
[variablenames, variablevalues]=textread(paramfile, '%s %f');
inputvalues=num2str(variablevalues);
input=strcat(variablenames, inputvalues);
[constantnames, constantvalues]=textread('constants.param', '%s %f');
cd(pranformLOCATION)
[spans, chords, offsets, airfoil]=textread(pranformfile, '%f %f %f %s');
cd(airfoilLOCATION)
[sref swet]=WhettedArea(spans, chords, offsets, airfoil, n);



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

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

Dragredux_BLI=(1-Dragredux_BLI/100);

wingspan=2*spans(end);
AR=wingspan^2/sref;
% wakeplot=findwake(airfoil, chords, spans, offsets, n, alt_cruise, M_cruise, laminarflow_cruise, fuse_width); 

%TAKEOFF------------------------------------------------------------------]
% cd0_to=findcd0(c_f, M_to, sref, swet);
[cd0_to D_wing_to D_fuse_to]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_to, alt_to, M_to, sref, c_f, D_engine, cd_landinggear, sfront_landinggear, fuse_width);
% finddragpolar(cd0_to, cLalpha_to, cL0, efficiency_to, AR, cd_correction, 'Takeoff')
% cLreq_to=findcLrequired(alt_to, sref, M_to, weight_to);
cLreq_to=.48;
alphareq_to=findalpharequired(cLalpha, cL0, cLreq_to);
cD_to=findcD(AR, efficiency_to, cd0_to, cLreq_to, cd_correction);
eloverdee_to=cLreq_to/cD_to;
cM_to=cMalpha*alphareq_to+cM0;
q_to=finddynpress(alt_to, M_to);
cD_to=cD_to*Dragredux_BLI;
cd0_to=cd0_to*Dragredux_BLI;
D_to=cD_to*q_to*sref;

%CRUISE------------------------------------------------------------------]
[cd0_cruise D_wing_cruise D_fuse_cruise]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_cruise, alt_cruise, M_cruise, sref, c_f, D_engine, 0, 0, fuse_width);
% finddragpolar(cd0_cruise, cLalpha, cL0, efficiency_cruise, AR, cd_correction, 'Cruise')
cLreq_cruise=findcLrequired(alt_cruise, sref, M_cruise, weight_cruise);
alphareq_cruise=findalpharequired(cLalpha, cL0, cLreq_cruise);
cD_cruise=findcD(AR, efficiency_cruise, cd0_cruise, cLreq_cruise, cd_correction);
cM_cruise=cMalpha*alphareq_cruise+cM0;
q_cruise=finddynpress(alt_cruise, M_cruise);
cD_cruise=cD_cruise*Dragredux_BLI;
cd0_cruise=cd0_cruise*Dragredux_BLI;
D_cruise=cD_cruise*q_cruise*sref;
eloverdee_cruise=cLreq_cruise/cD_cruise;

%LOITER------------------------------------------------------------------]
[cd0_loiter D_wing_loiter D_fuse_loiter]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_loiter, alt_loiter, M_loiter, sref, c_f, D_engine, 0, 0, fuse_width);
% finddragpolar(cd0_loiter, cLalpha_loiter, cL0, efficiency_loiter, AR, cd_correction, 'Loiter')
cLreq_loiter=findcLrequired(alt_loiter, sref, M_loiter, weight_loiter);
alphareq_loiter=findalpharequired(cLalpha, cL0, cLreq_loiter);
cD_loiter=findcD(AR, efficiency_loiter, cd0_loiter, cLreq_loiter, cd_correction);
cM_loiter=cMalpha*alphareq_loiter+cM0;
q_loiter=finddynpress(alt_loiter, M_loiter);
cD_loiter=cD_loiter*Dragredux_BLI;
cd0_loiter=cd0_loiter*Dragredux_BLI;
D_loiter=cD_loiter*q_loiter*sref;
eloverdee_loiter=cLreq_loiter/cD_loiter;

%LANDING------------------------------------------------------------------]
[cd0_landing D_wing_landing D_fuse_landing]=findcd0(spans, chords, offsets, airfoil, n, laminarflow_landing, alt_landing, M_landing, sref, c_f, D_engine, cd_landinggear, sfront_landinggear, fuse_width);
% finddragpolar(cd0_landing, cLalpha_landing, cL0, efficiency_landing, AR, cd_correction, 'Landing')
% cLreq_landing=findcLrequired(alt_landing, sref, M_landing, weight_landing);
cLreq_landing=.401;
alphareq_landing=findalpharequired(cLalpha, cL0, cLreq_landing);
cD_landing=findcD(AR, efficiency_landing, cd0_landing, cLreq_landing, cd_correction);
cM_landing=cMalpha*alphareq_landing+cM0;
q_landing=finddynpress(alt_landing, M_landing);
cD_landing=cD_landing*Dragredux_BLI;
cd0_landing=cd0_landing*Dragredux_BLI;
D_landing=cD_landing*q_landing*sref;
eloverdee_landing=cLreq_landing/cD_landing;

cd(outputsLOCATION)
filenameplan=regexprep(pranformfile, '.plan', '_planform.tiff');
filenameoutput=regexprep(pranformfile, '.plan', '.txt');
plottitle=regexprep(pranformfile, '.plan', '');

fid = fopen(filenameoutput, 'w');
plottitle

fprintf(fid, plottitle);



fprintf(fid,'\nOutputs--------------------------------------------------\n');

fprintf(fid, 'Wingspan: %f, AR: %f, Planform Area: %f, Wetted Area: %f \n', wingspan, AR, sref, swet);

fprintf(fid, '\n----------Takeoff--------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_to, cd0_to, cLreq_to, cD_to, D_to, cM_to, eloverdee_to);

fprintf(fid, '\n----------Cruise---------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f,  CM: %f, L/D: %f \n D_Wing:%f, D_fuselage:%f', alphareq_cruise, cd0_cruise, cLreq_cruise, cD_cruise, D_cruise, cM_cruise, eloverdee_cruise, D_wing_cruise, D_fuse_cruise);

fprintf(fid, '\n----------Loiter---------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_loiter, cd0_loiter, cLreq_loiter, cD_loiter, D_loiter, cM_loiter, eloverdee_loiter);

fprintf(fid, '\n----------Landing--------------------------------\n alfa: %f, Cd0: %f, CL: %f, CD: %f, D: %f, CM: %f, L/D: %f \n', alphareq_landing, cd0_landing, cLreq_landing, cD_landing, D_landing, cM_landing, eloverdee_landing);
fprintf(fid, '\nInputs--------------------------------------------------\n');
for k=1:length(input)-1
   fprintf(fid, '%s      %s \n',input{k}, input{k+1});
end


fclose(fid);
cd(rootLOCATION)
end
end

% function cd0=findcd0(c_f, M, sref, swet)
% if M>.7
% c_f=c_f*(1-.09*M^2);
% end
% cd0=c_f*(swet/sref);
% end
function [cd0 D_wing D_fuse]=findcd0(spans, chords, offsets, airfoil, n, laminarflow, alt, v_mach, sref, c_f, D_engine, cd_lg, Sfront_lg, width)
wingSpan = 1;
% n = 50;
width=width/2;

[q v_fps dynvisc]=findconstants(alt, v_mach);
drag=0;
drag_fuse=0;
drag_wing=0;

Chords = findChords(n, chords, spans);
for i = 1:(size(chords)-1)
    for j = 1:n
        dist = wingSpan*(spans(i+1)-spans(i))./n;
                chordloc=linspace(0, Chords(i,j));
                delta_x = Chords(i,j)/100;
                
                if spans(i)>width
                    for k = 1:100
                        x(k)= k * Chords(i,j)/100;
                        Re_x(k)= v_fps * x(k)/dynvisc;
                        if (k <= 100*laminarflow )
                            cf_x(k) = 0.664 /(sqrt(Re_x(k)));
                        else % Turbulent
                            cf_x(k) = 0.455 /(log(0.06 * Re_x(k))^2);
                        end
                        xcf_x(k)=x(k)*cf_x(k);
                    end
                    t = trapz(cf_x);
                    Cf = 2 * delta_x * t /Chords(i,j);
                    drag_section = q*Cf*findPerim(airfoil{i}, Chords(i,j), wingSpan).*dist;
                    drag_wing=drag_wing+drag_section;
                else
                    for k = 1:100
                        x(k)= k * Chords(i,j)/100;
                        Re_x(k)= v_fps * x(k)/dynvisc;
                        if (k <= 100*laminarflow -10)
                            cf_x(k) = 0.664 /(sqrt(Re_x(k)));
                        else % Turbulent
                            cf_x(k) = 0.455 /(log(0.06 * Re_x(k))^2);
                        end
                        xcf_x(k)=x(k)*cf_x(k);
                    end
                    t = trapz(cf_x);
                    Cf = 2 * delta_x * t /Chords(i,j);
                    drag_section = q*Cf*findPerim(airfoil{i}, Chords(i,j), wingSpan).*dist;
                    drag_fuse=drag_fuse+drag_section;
                end
                t = trapz(cf_x);
                Cf = 2 * delta_x * t /Chords(i,j);
                drag_section = q*Cf*findPerim(airfoil{i}, Chords(i,j), wingSpan).*dist;
                drag=drag+drag_section;
         
    end
end

if v_mach>.7
c_f=c_f*(1-.09*v_mach)^2;
end
% cd0engine=c_f*(s_engine/sref);
D_lg=Sfront_lg*q*cd_lg;
D_wing=drag_wing*2;
D_fuse=drag_fuse*2;
drag_tot=D_wing+D_engine+D_lg+D_fuse;
cd0=drag_tot/(q*sref);
cD_lg=D_lg/(q*sref);
cD_engine=D_engine/(q*sref);
cD_wing=D_wing/(q*sref);
cD_fuse=D_fuse/(q*sref);
end
function finddragpolar(cd0, cl_alfa, cL0, e, ar, cd_correction, filename)
alphadeg=[-12:10];
% 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
grid on
axisX0=linspace(0,max(cd));
axisY0=zeros(1,length(axisX0));
plot(axisX0, axisY0,'k--');
plot(C_d, C_L, 'k', 'LineWidth', 2);
plot(C_d, C_L, 'k.');
title(graphname);
xlabel('C_d', 'FontSize', 16);
ylabel('C_L', 'FontSize', 16);
axis([0, .02, -.3, .4]);
set(gca, 'FontSize',16)

% for i=1:length(alphadeg)
%     alphalabel=num2str(alphadeg(i));
%     label=['  \alpha =' alphalabel];
% text(C_d(i), C_L(i), label)
% end

print ('-dtiff', filename)
grid off
close(1)
end
function q=finddynpress(alt, v_mach)
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

%slpcuft: clugs ber cubic foot
density_slpcuft=press/(1718*(temp+459.7));
v_fps=v_mach*sqrt(1.4*1718*(temp+459.7));
q=.5*density_slpcuft*v_fps^2;
end
function [q v_fps dynvisc]=findconstants(alt, v_mach)
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

%slpcuft: clugs ber cubic foot
density_slpcuft=press/(1718*(temp+459.7));
v_fps=v_mach*sqrt(1.4*1718*(temp+459.7));
q=finddynpress(alt, v_mach);

S = 110.4;
beta = 1.458e-6;
temp_k=[(temp-32) / (1.8)] + 273 ;
viscosity_pascalseconds = (beta*temp_k^(1.5))/(temp_k+S);
%kgpmcu:kilograms per meter cubed
%msqpse:meters squared per second
%sqftpsec:square feet per second
density_kgpmcu=density_slpcuft*515.4;
dynamic_viscosity_msqpsec=viscosity_pascalseconds/density_kgpmcu;
dynamic_viscosity_sqftpsec=dynamic_viscosity_msqpsec*10.76;
dynvisc=dynamic_viscosity_sqftpsec;
density=density_slpcuft;
end
function cLreq=findcLrequired(alt, s, v_mach, lift)
q=finddynpress(alt, v_mach);

cLreq=lift/(q*s);

end
function alpha=findalpharequired(cLalpha, cL0, cLrequired)
alpha=(cLrequired-cL0)/cLalpha;
end
function cD=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);
        cD=cd0+cd_i+cd_correction;
end
function [pranformArea WetArea]=WhettedArea(spans, chords, offsets, airfoil, n)
% fileName = 'AirData.mat';
wingSpan = 1;
% n = 50;

engineArea=0;
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;
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 [wakeplot]=findwake(airfoil, chords, spans, offsets, n, alt, v_mach, laminarflow, width) 
[q v_fps dynvisc]=findconstants(alt, v_mach);
Chords = findChords(n, chords, spans);
wakeplot=zeros(100,1);
wakeplotX=zeros(100,1);
wakeplotY=zeros(100,1);
width=width/2;
% for i = 1:(size(chords)-1)
for i=1:1:(size(chords)-1)
    for j = 1:n
%         wakecoord(j,:)=linspace(0, Chords(i,j), n);
dist = (spans(i+1)-spans(i))./n;
y(j)=spans(i)+j*dist;

        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);

        if spans(i)>width
            lam=laminarflow;
        else
            lam=laminarflow-.15;
        end
        x=linspace(0, Chords(i,j));
        for k = 1:length(x)
            if (k <= 100*lam)
                delta_x=( 5.2 * x(k) )*(v_fps*x(k)/dynvisc)^-0.5;
                wake(k,j)= delta_x;
            else %if turbulent
                x_tr=x(100*lam);
                Re_xtr = v_fps * x_tr / dynvisc;
                delta_xtr=(5.2*x_tr)*(v_fps*x(k)/dynvisc)^-(0.5);
                delta_l = ((delta_xtr*((v_fps)^0.2))/(0.37*((dynvisc)^0.2)))^(1/0.8);
                lt = x(k) - x_tr + delta_l;
                Re_lt = v_fps * lt / dynvisc;
                delta = 0.37 * lt /(Re_lt)^0.2;
                wake(k,j) = delta;
            end

wakeX(k,j)=x(k)+yplotLE(i,j);
wakeY(k,j)=y(j);
        end
        
        for k = 1:100

            wakeX(k,j)=x(k)+yplotLE(i,j);
            wakeY(k,j)=y(j);
        end
        
        

wakeplot=[wakeplot, wake];
wakeplotX=[wakeplotX, wakeX];
wakeplotY=[wakeplotY, wakeY];
    end

end
wakeplot(:,1)=[ ];
wakeplotX(:,1)=[ ];
wakeplotY(:,1)=[ ];
wakeplot(1,:)=[0];
flipud(wakeplot);
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
