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'};
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


% repeat this for all columns in input.
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
    if contains(p.Results.energyScale,'kinetic')
        % LINEAR FIT
        if contains(p.Results.mode,'linfit')
            disp('linfit not yet supported')
            % MEAN
        else
            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);


        end

    elseif contains(p.Results.energyScale,'binding')
        if contains(p.Results.mode,'linfit')
            disp('linfit not yet supported')
            % MEAN
        else
            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' );
                h = 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
DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)                  = tmp;
% 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 p.Results.energyEnd],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