function DataTable = nappBackground(DataTable,varargin)
%nappBackground determines and subtracts the background from PE or XAS spectra
%
%   FileList = nappBackground(FileList,energyStart,energyEnd);
%
%   FileList = nappBackground(FileList,368,370,'tableSelection',tmpDataSelection,'energyScale','kinetic',...
%   'mode','mean','showFigures','off');
%
%
%   *input*
%   FileList: data table
%   energyStart: smallest value of energy range to use for background calculation
%       (can also be binding energy)
%   energyEnd: largest value of energy energy range to use for background calculation
%       (can also be binding energy)
%   tableSelection (optional): logical array to filter data table
%   tableVariable: choose to process either
%       'spectrum_counts','spectrum_countsPerSecond' (default), or
%       'spectrum_countsNorm'.
%   mode: mean(default) subtract the mean of the energy range from
%   spectrum; linfit subtract a linear regression through the energy range
%   from spectrum
%   showFigures, on(default): show fit in figure, off do not show fit in
%   figure
%   energyScale: kinetic(default): work with kinetic energy scale; binding
%   works with binding energy scale.
%
%
%   *output*
%       The table variable passed to the function as tableVariable
%       parameter will be overwritten with the background subtracted spectrum.
%
%       The function parameters (mode, energyStart, energyEnd) are written
%       to the table as log.



% Define conditions and defaults of input
defaultTableVariable        = 'spectrum_countsPerSecond';
expectedTableVariable       = {'spectrum_counts','spectrum_countsPerSecond','spectrum_countsNorm'};
defaultMode                 = 'mean';
expectedMode                = {'mean','linfit','exponential','mean, divide','linfit, divide','exponential, divide'};
defaultTableSelection       = true(size(DataTable,1),1);
defaultShowFigures          = 'on';
expectedShowFigures         = {'on','off'};
defaultEnergyScale          = 'kinetic';
expectedEnergyScale         = {'kinetic','binding','xRay'};

% Parse input
p = inputParser;
addRequired(p,'DataTable',@istable);
addRequired(p,'energyStart',@isnumeric)
addRequired(p,'energyEnd',@isnumeric)
addOptional(p,'tableSelection',defaultTableSelection,@islogical);
addParameter(p,'tableVariable',defaultTableVariable,...
    @(x) any(validatestring(x,expectedTableVariable)));
addParameter(p,'mode',defaultMode,...
    @(x) any(validatestring(x,expectedMode)));
addParameter(p,'showFigures',defaultShowFigures,...
    @(x) any(validatestring(x,expectedShowFigures)));
addParameter(p,'energyScale',defaultEnergyScale,...
    @(x) any(validatestring(x,expectedEnergyScale)));
parse(p,DataTable,varargin{:});

% check input (all elements of input variable need to have non-zero size)
if ~all(cell2mat(cellfun(@size,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable),'UniformOutput',false)),'all')
    disp('The table variable appears to be empty. I quit now')
    error('The table variable appears to be empty.')
end

% add empty table variables (if they don't exist)
existColumn = strcmp('logBackgroundRange',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.logBackgroundRange = NaN(size(DataTable,1),2);
end
existColumn = strcmp('logBackgroundMode',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.logBackgroundMode = cellstr(repmat(' ',size(p.Results.DataTable,1),1));
end
existColumn = strcmp('logBackgroundDate',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.logBackgroundDate = cellstr(repmat(' ',size(p.Results.DataTable,1),1));
end
% add empty variable to table
existColumn = strcmp('spectrum_countsNorm',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.spectrum_countsNorm = cell(size(DataTable,1),1);
end


% repeat this for all columns in input.
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)

    if startsWith(p.Results.mode,'linfit')
        ft = fittype( 'poly1' );
    elseif contains(p.Results.mode,'exponential')
        ft = fittype( 'exp2' );
    end

    if contains(p.Results.energyScale,'kinetic')

        if ~startsWith(p.Results.mode,'mean')
            [xData, yData] = prepareCurveData(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii},p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});

            excludedPoints = excludedata( xData, yData, 'domain',[p.Results.energyStart(1) p.Results.energyEnd(1)] );
            if size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==1
                excludedPointsEins = excludedPoints;
                clearvars excludedPoints
                excludedPointsZwei = excludedata( xData, yData, 'domain',[xData(1) p.Results.energyStart(2)] );
                excludedPoints = excludedPointsEins & excludedPointsZwei;
            elseif size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==2
                excludedPointsEins = excludedPoints;
                clearvars excludedPoints
                excludedPointsZwei = excludedata( xData, yData, 'domain',[xData(1) p.Results.energyStart(2)] );
                excludedPointsDrei = excludedata( xData, yData, 'domain',[p.Results.energyEnd(2) xData(end)] );
                excludedPoints = excludedPointsEins & excludedPointsZwei & excludedPointsDrei;
            end

            opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', ~excludedPoints );
            opts.Display = 'Off';

            % Fit model to data.
            [fitresult, ~] = fit( xData, yData, ft, opts );
            if contains(p.Results.showFigures,'on')
                % Plot fit with data.
                figure( 'Name', 'untitled fit 1' );
                plot( fitresult, xData, yData, ~excludedPoints );
                % legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
                % Label axes
                xlabel( 'kinetic energy [eV]', 'Interpreter', 'none' );
                ylabel( 'PE intensity', 'Interpreter', 'none' );
                grid on
                title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
            end
            if endsWith(p.Results.mode,'divide')
                tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}./fitresult(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});

            else
                tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});
            end
        elseif startsWith(p.Results.mode,'mean')
            b_1         = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= p.Results.energyStart,1);
            b_2         = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= p.Results.energyEnd,1);

            if b_2 == 1
                b_2=10;
            end

            b_mean   = mean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(b_1:b_2),'omitnan');
            tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false); %#ok<*AGROW>


        end

    elseif contains(p.Results.energyScale,'binding')
        if ~startsWith(p.Results.mode,'mean')
            [xData, yData] = prepareCurveData(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii},p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});

            % Set up fittype and options.

            excludedPoints = excludedata( xData, yData, 'domain',[p.Results.energyStart(1) p.Results.energyEnd(1)] );
            if size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==1
                excludedPointsEins = excludedPoints;
                clearvars excludedPoints
                excludedPointsZwei = excludedata( xData, yData, 'domain',[xData(1) p.Results.energyStart(2)] );
                excludedPoints = excludedPointsEins & excludedPointsZwei;
            elseif size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==2
                excludedPointsEins = excludedPoints;
                clearvars excludedPoints
                excludedPointsZwei = excludedata( xData, yData, 'domain',[p.Results.energyStart(2) p.Results.energyEnd(2)] );
                excludedPoints = excludedPointsEins & excludedPointsZwei;
            end

            opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', ~excludedPoints );
            opts.Display = 'Off';



            % Fit model to data.
            [fitresult, ~] = fit( xData, yData, ft, opts );
            if contains(p.Results.showFigures,'on')
                % Plot fit with data.
                figure( 'Name', 'untitled fit 1' );
                plot( fitresult, xData, yData, excludedPoints );
                % legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
                % Label axes
                xlabel( 'binding energy [eV]', 'Interpreter', 'none' );
                ylabel( 'PE intensity', 'Interpreter', 'none' );
                grid on
                title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
            end
            if endsWith(p.Results.mode,'divide')
                tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}./fitresult(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii});

            else
                tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii});
            end
        elseif startsWith(p.Results.mode,'mean')
            b_1         = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= p.Results.energyStart,1);
            b_2         = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= p.Results.energyEnd,1);

            if b_2 == 1
                b_2=10;
            end

            b_mean   = mean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(b_2:b_1),'omitnan');
            tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false);


        end
    elseif contains(p.Results.energyScale,'xRay')
        if contains(p.Results.mode,'linfit') % LINEAR FIT

            % Fit: 'untitled fit 1'.
            [xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
                p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
            % Set up fittype and options.
            ft = fittype( 'poly1' );
            excludedPoints = excludedata( xData, yData, 'domain',[p.Results.energyStart p.Results.energyEnd] );
            fo = fitoptions('Exclude', excludedPoints);
            % Fit model to data.
            [fitresult, ~] = fit( xData, yData, ft, fo );
            if contains(p.Results.showFigures,'on')
                % Plot fit with data.
                figure( 'Name', 'fit pre-edge' );
                plot( fitresult, xData, yData, excludedPoints );
                xlabel('Photon Energy eV')
                ylabel('partial A-M.E.Y.')
                grid on
                title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
            end
            tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii});

        elseif contains(p.Results.mode,'mean') % mean

            if p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(1) < p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(end)
                % prepare Data

                [xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
                    p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
                % MEAN
                b_1 = find(xData >= p.Results.energyStart,1);
                b_2 = find(xData >= p.Results.energyEnd,1);

                b_mean = mean(yData(b_1:b_2));
            elseif p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(1) > p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(end)
                % prepare Data
                [xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
                    p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
                % MEAN
                b_1 = find(xData <= p.Results.energyStart,1);
                b_2 = find(xData <= p.Results.energyEnd,1);

                b_mean = mean(yData(b_2:b_1));
            else
                disp('Warning: Issue with NEXAFS data or energies')
            end


            tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false);
        end


    end
end

% fill selected rows
if endsWith(p.Results.mode,'divide')
    DataTable(p.Results.tableSelection,:).spectrum_countsNorm                        = tmp;
else
    DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)                  = tmp;
end
% write log to table
DataTable.logBackgroundMode(p.Results.tableSelection,:)                     = cellstr(repmat(p.Results.mode,size(p.Results.DataTable(p.Results.tableSelection,:),1),1));
DataTable.logBackgroundRange(p.Results.tableSelection,:)                    = repmat([p.Results.energyStart(1) p.Results.energyEnd(1)],size(p.Results.DataTable(p.Results.tableSelection,:),1),1);
DataTable.logBackgroundDate(p.Results.tableSelection,:)                     = cellstr(repmat(datestr(now),size(p.Results.DataTable(p.Results.tableSelection,:),1),1));

end