function [a,CL,CD,CM] = xfoil(varargin)
% XFOIL   Run XFOIL sequences and return the results
% Alex Stoll / 28 February 2012
%
%    [A,CL,CD,CM] = XFOIL(X,Y,RE,M,PARAMETER,VALUES) runs XFOIL with the
%    provided airfoil coordinates (X,Y). These coordinates start at the trailing
%    edge of the upper surface, trace the upper surface to the leading edge, and
%    then continue along the lower surface to end at the trailing edge. RE = 0
%    runs XFOIL inviscidly. PARAMETER is either 'a' for alpha (in degrees) or
%    'c' for Cl, and VALUES is the value(s) of this parameter to run.
%    (Regardless of the value of the vector VALUES, the vector actually used is
%    linspace(min(VALUES),max(VALUES),length(VALUES)).) A, CL, CD, and CM are
%    vectors of the XFOIL results; if XFOIL does not converge at one of the
%    parameter values, the corresponding results will either be missing from
%    these vectors or will be Inf.
%    GDES\CADD is run on the airfoil to smooth it out, and PANE is run to
%    repanel it. If this behavior is unwanted, it must be removed/commented out
%    in this file. Temporary files (script, airfoil coordinates, output, and
%    polar) are written in C:\ - this behavor can be changed by editing this
%    file.
%
%    For example, [A,CL,CD,CM] = XFOIL(X,Y,1e6,0.3,'a',linspace(-5,5,11)) runs
%    XFOIL with the provided coordinates X and Y at a Reynolds number of 1e6 and
%    a Mach number of 0.3 at a sequence of 11 alphas from -5 to 5.
%
%    XFOIL(X,Y,RE,M,PARAMETER,VALUES,ID) suffixes the value of ID to the ends of
%    the file names of all the temporary files written in this function. This is
%    useful when running it in parallel, so that the same temporary files are
%    not attempted to be written to simultaneously. ID can be an integer, float,
%    or string.
%
%    XFOIL(AIRFOIL,RE,M,PARAMETER,VALUES) uses the airfoil file in the string
%    AIRFOIL instead of explicitly providing coordinates.
%
%    XFOIL(AIRFOIL,RE,M,PARAMETER,VALUES,ID) uses the airfoil file in the string
%    AIRFOIL and suffixes the value of ID to the filenames of the temporary
%    files.
%
%    For example, [A,CL,CD,CM] = XFOIL('n0012.dat',0,0.5,'c',0) runs the airfoil
%    file 'n0012.dat' inviscidly at Mach 0.5 at Cl = 0 to find the zero-lift
%    angle of attack (in degrees).
%
%    This script works on Windows. Does it work on other operating systems? Heck
%    if I know. I imagine you will at least have to change the "C:\" part.

% Set this to your XFOIL directory. This is the directory containing "bin"
% which contais the "xfoil" executable.
path = '.\XFOIL';

id = [];
if ischar(varargin{1})
    writeairfoil = false;
    afname = varargin{1};
    Re = varargin{2};
    M = varargin{3};
    varname = varargin{4};
    vars = varargin{5};
    if nargin>5, id = varargin{6}; end
    if ~ischar(id), id = num2str(id); end
    n = 1;
else
    writeairfoil = true;
    x = varargin{1};
    y = varargin{2};
    if size(y,1) == 1, y = y'; end
    n = size(y,2);
    Re = varargin{3};
    M = varargin{4};
    varname = varargin{5};
    vars = varargin{6};
    if nargin>6, id = varargin{7}; end
    if ~ischar(id), id = num2str(id); end
    afname = ['C:\\airfoil' id '.tmp'];
    
    % Write airfoil coordinates file in MSES format
    fp = fopen(afname,'w');
    fprintf(fp,['xfoil.m ' datestr(now) '\r\n']);
    fprintf(fp,' -2.0  3.0  -2.5  3.5\r\n');
    for i=1:n
        for j=1:size(y,1)
            fprintf(fp,['     ' num2str(x(j)) '    ' num2str(y(j,i)) '\r\n']);
        end
        if i<n, fprintf(fp,'999.0 999.0\r\n'); end
    end
    fclose(fp);
end
lv = length(vars);
zspan = vars(lv)/vars(1) < 0;   % If VARS spans 0, run XFOIL differently

% Write XFOIL script
inpname = ['C:\xfoil' id '.inp'];
fp = fopen(inpname,'w');
fprintf(fp, 'y\r\n');       % Yes, load xfoil.def

% Turn off plotting
fprintf(fp, 'plop\r\n');
fprintf(fp, 'G\r\n\r\n');

for i=1:n % For each airfoil
    % Load airfoil
    fprintf(fp, ['load ' afname '\r\n']);
    if n>1, fprintf(fp, [int2str(i) '\r\n']); end
    
    % Smooth out the airfoil for higher chance of convergence
    fprintf(fp, 'gdes\r\n');
    fprintf(fp, 'cadd\r\n');
    fprintf(fp,'\r\n'); % accept default corner angle criterion
    fprintf(fp,'1\r\n'); % uniform spline parameter type
    fprintf(fp,'\r\n'); % accept default x limits
    fprintf(fp, '\r\npane\r\n');
    
    % Set operating point
    fprintf(fp, 'oper\r\n');
    if Re(i) > 0, fprintf(fp, ['visc ' int2str(round(Re(i))) '\r\n']); end
    fprintf(fp, ['Mach ' num2str(M(i)) '\r\n']);
    fprintf(fp, 'iter 100\r\n');
    fprintf(fp, 'pacc\r\n');
    polarfile = ['C:\\clcd' int2str(i) '-' id '.txt'];
    fprintf(fp, [polarfile '\r\n\r\n']);
    
    % Do the calculation
    if lv == n || lv == 1
        fprintf(fp, [varname ' ' num2str(vars(i)) '\r\n']);
    else
        dv = mean(diff(vars));
        [~, iv0] = min(abs(vars));
        v0 = vars(iv0);
        if zspan
            fprintf(fp, [varname(1) 'seq ' num2str(v0) ' ' num2str(vars(1)) ' ' num2str(dv) '\r\n']);
            fprintf(fp, 'vpar\r\n');
            fprintf(fp, 'init\r\n\r\n');
            fprintf(fp, [varname(1) 'seq ' num2str(v0) ' ' num2str(vars(lv)) ' ' num2str(dv) '\r\n']);
        else
            fprintf(fp, [varname(1) 'seq ' num2str(vars(1)) ' ' num2str(vars(lv)) ' ' num2str(dv) '\r\n']);
        end
    end
    
    % Reset pacc and visc
    fprintf(fp, 'pacc\r\n');
    if Re(i) > 0, fprintf(fp, 'visc\r\n'); end
    fprintf(fp, '\r\n');
    
    % Before running, delete the polar file if it already exists.
    if exist(polarfile,'file'), delete(polarfile); end
end
fprintf(fp, 'quit\r\n');
fclose(fp);

% Run XFOIL and find CL, CD, etc.
outname = ['C:\xfoilout' id '.txt'];  % Save the output for troubleshooting
dos([path '\bin\xfoil < ' inpname ' > ' outname]);
% Read in the output from the polar file(s)
if n>1
    CL = zeros(1,n);
    CD = CL;
    a = CL;
    CM = CL;
    for i=1:n
        coeffs = textread(['C:\clcd' int2str(i) '-' id '.txt'], '%s');
        if length(coeffs)<53,
            a(i) = Inf;
            CL(i) = Inf;
            CD(i) = Inf;
            CM(i) = Inf;
        else
            a(i) = str2double(coeffs(49));
            CL(i) = str2double(coeffs(50));
            CD(i) = str2double(coeffs(51));
            CM(i) = str2double(coeffs(53));
        end
        delete(['clcd' int2str(i) '.txt'])
    end
else
    coeffs = reshape(textread(['C:\clcd1-' id '.txt'], '%n','headerlines',12),7,[])';
    [a,o] = unique(coeffs(:,1));
    CL = coeffs(o,2);
    CD = coeffs(o,3);
    CM = coeffs(o,5);
end

% Delete temporary files
if writeairfoil, delete(afname); end    % Delete airfoil coordinates file
delete(inpname);                        % Delete script
delete(outname);                        % Delete output file
for i=1:n                               % Delete polar files
    delete(['C:\clcd' int2str(i) '-' id '.txt'])
end