function DataTable = nappScale(DataTable,varargin)
%NAPPSCALE Normalize XPS and NEXAFS spectra. 
%   The peak of the same region or of a different region with identical
%   SES file number (same sequence iteration) can be used for XPS data. 
%   NEXAFS spectra can be normed to their area or to the intensity
%   at an energy range of choice.
%
%   For XPS ('.ibw') data, NAPPSCALE divides spectra of one region
%   (regionOld, default is C1s) by the fitted peak area taken from a
%   spectrum of the same sequence but different region (regionNew, default
%   is O1s). The same sequence iteration is identified by the  SES file
%   number and acquisition time. For NEXAFS ('.h5') data, NAPPSCALE
%   normalizes spectra to the mean intensity or to the integrated area over
%   a chosen photon-energy interval.
%
%   DataTable = NAPPSCALE(DataTable) normalizes spectra in the input table
%   using default options and returns the updated table with
%   normalized spectra in the variable spectrum_countsNorm.
%
%   DataTable = NAPPSCALE(DataTable, tableSelection) restricts
%   normalization to rows where tableSelection is true. The length
%   of tableSelection must match the height of DataTable.
%
%   DataTable = NAPPSCALE(DataTable, tableSelection, energyStart, energyEnd)
%   additionally specifies an energy interval used to normalize NEXAFS
%   spectra ('.h5' data)
%
%   DataTable = NAPPSCALE(..., Name, Value) customizes the source variable,
%   normalization mode, line region mapping, and peak selection via
%   Name-Value arguments (see below):
%

%
%   Required input
%   --------------
%   DataTable
%       Table containing spectroscopy data and metadata.[web:2]
%       Depending on file type and options, the function uses:
%         - name                  : file name, used to distinguish 'ibw'
%                                   and 'h5' and to construct regionNew
%                                   file names (XPS)
%         - spectrum_counts       : 1D spectrum (cell of numeric)
%         - spectrum_countsPerSecond
%                                 : 1D spectrum normalized to rate
%         - spectrum_countsNorm   : (created/overwritten) normalized
%                                   spectrum
%         - peakArea              : numeric vector of fitted peak areas
%                                   (XPS, one row per region)
%         - lastModifiedDatenum   : datenum time stamp used to find the
%                                   matching regionNew spectrum
%         - xRayEkin_eV           : photon energy axis (NEXAFS)
%
%   Optional positional input
%   -------------------------
%   tableSelection
%       Logical vector selecting rows in DataTable to normalize.
%       Default: all rows true
%
%   energyStart, energyEnd
%       Photon-energy limits (in eV) that define the normalization
%       interval for NEXAFS spectra ('.h5').
%       Type   : numeric scalars
%       Use    : required for 'intensity' or 'area' mode on '.h5' data.
%
%   Name-Value arguments
%   --------------------
%   'tableVariable'
%       Name of the spectral variable to be normalized.[web:10]
%       Type   : char or string
%       Values : 'spectrum_countsPerSecond' (default)
%                'spectrum_counts'
%                'spectrum_countsNorm'
%
%   'regionOld'
%       XPS region name of the spectra to be normalized ('.ibw').
%       Used for file-name replacement when searching the reference
%       region
%       Type   : char or string
%       Values : 'C1s' (default) | 'O1s' | 'N1s' | 'Cl2p'
%
%   'regionNew'
%       XPS region name that provides the reference peak area used for
%       normalization.
%       Type   : char or string
%       Values : 'O1s' (default) | 'C1s' | 'N1s'
%
%   'peakNumber'
%       Index of the peak in the peakArea vector used as normalization
%       reference for XPS data.
%       Type   : positive integer scalar
%       Default: 1
%
%   'mode'
%       Normalization mode.
%       Type   : char or string
%       Values : 'area'      (default)
%                'intensity'
%
%       Behaviour:
%         - For '.ibw' (XPS) 
%             'area'      : divide selected spectra by the chosen peakArea
%                           of the corresponding regionNew spectrum.
%             'intensity' : currently not implemented (warning printed).
%         - For '.h5' (NEXAFS) 
%             'intensity' : divide spectra by the mean intensity over
%                           [energyStart, energyEnd].
%             'area'      : divide spectra by the integrated area under
%                           the curve over [energyStart, energyEnd]
%                           (trapz).
%
%   Output
%   ------
%   DataTable
%       Updated input table with normalized spectra and log variables.
%       The function ensures the following variables:
%
%       spectrum_countsNorm
%           Cell array containing the normalized spectrum for each
%           processed row. Existing values are overwritten.
%
%       logScaleRange
%           Numeric matrix with two columns. For '.h5' rows, stores the
%           [energyStart energyEnd] interval used for normalization.
%
%       logScaleMode
%           Cell array of character vectors with the mode ('area' or
%           'intensity') used for each normalized row.
%
%       logScaleDate
%           Cell array of character vectors containing the time stamp of
%           the normalization (datestr(now)).
%
%   Notes
%   -----
%   - All selected entries in the chosen tableVariable must contain
%     non-empty spectra; otherwise, the function throws an error.
%   - For XPS data, regionOld and the tableSelection filter must refer to
%     the same region; regionNew is found by name replacement and by
%     matching date and time to locate the correct reference spectrum.
%   - For NEXAFS data, the function supports both increasing and decreasing
%     photon-energy axes when locating the normalization interval.



% Define conditions and defaults of input
defaultTableVariable    = 'spectrum_countsPerSecond';
expectedTableVariable   = {'spectrum_counts','spectrum_countsPerSecond','spectrum_countsNorm'};
defaultMode             = 'area';
expectedMode            = {'intensity','area'};
defaultRegionOld        = 'C1s';
expectedRegionOld       = {'O1s','C1s','N1s','Cl2p'};
defaultRegionNew        = 'O1s';
expectedRegionNew       = {'O1s','C1s','N1s'};
defaultPeakNumber       = 1;
defaultTableSelection   = true(size(DataTable,1));


% Parse input
p = inputParser;
addRequired(p,'DataTable',@istable);
addOptional(p,'tableSelection',defaultTableSelection,@islogical);
addOptional(p,'energyStart',@isnumeric)
addOptional(p,'energyEnd',@isnumeric)
addParameter(p,'tableVariable',defaultTableVariable,...
    @(x) any(validatestring(x,expectedTableVariable)));
addParameter(p,'regionOld',defaultRegionOld,...
    @(x) any(validatestring(x,expectedRegionOld)));
addParameter(p,'regionNew',defaultRegionNew,...
    @(x) any(validatestring(x,expectedRegionNew)));
addParameter(p,'peakNumber',defaultPeakNumber);
addParameter(p,'mode',defaultMode,...
    @(x) any(validatestring(x,expectedMode)));
parse(p,DataTable,varargin{:});

% add empty variable to table
existColumn = strcmp('spectrum_countsNorm',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.spectrum_countsNorm = cell(size(DataTable,1),1);
end

existColumn = strcmp('logScaleRange',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.logBackgroundRange = NaN(size(DataTable,1),2);
end
existColumn = strcmp('logScaleMode',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.logBackgroundMode = cellstr(repmat(' ',size(p.Results.DataTable,1),1));
end
existColumn = strcmp('logScaleDate',DataTable.Properties.VariableNames);
if ~any(existColumn)
    DataTable.logBackgroundDate = cellstr(repmat(' ',size(p.Results.DataTable,1),1));
end

% 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


% repeat this for all columns in input.
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
    if endsWith(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'ibw') || endsWith(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'xy')
        % get name and date of input file and replace old region with new
        % region. This will be the files where the peak area is taken from.
        tmpName  = replace(p.Results.DataTable(p.Results.tableSelection,:).name{ii}, ...
            p.Results.regionOld,p.Results.regionNew);
        tmpDate  = p.Results.DataTable(p.Results.tableSelection,:).lastModifiedDatenum(ii);
        % select spectra based on the new names and dates of original file
        tmpDataSelection =  contains(p.Results.DataTable.name,tmpName)...
            & ismember(year(p.Results.DataTable.lastModifiedDatenum),year(tmpDate))...
            & ismember(month(p.Results.DataTable.lastModifiedDatenum),month(tmpDate))...
            & ismember(day(p.Results.DataTable.lastModifiedDatenum),day(tmpDate))...
            & ismember(hour(p.Results.DataTable.lastModifiedDatenum),hour(tmpDate))...
            & ismember(minute(p.Results.DataTable.lastModifiedDatenum),minute(tmpDate));
        % check peak area 
        existColumn = strcmp('peakArea',DataTable.Properties.VariableNames);
        if any(existColumn)
                if isnan(DataTable( tmpDataSelection,:).peakArea(p.Results.peakNumber))
                    disp('The peak area appears to be empty. I quit now')
                    disp(DataTable( tmpDataSelection,:).name{ii})
                    error('The peak area appears to be empty.')
                end
        else
              disp('The data table appear to have not peak area column. I quit now')
                    error('No peak area column')  
        end
        % do the norming
        if contains(p.Results.mode,'area')
            DataTable(p.Results.tableSelection,:).spectrum_countsNorm{ii} = ...
                DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}...
                ./  DataTable( tmpDataSelection,:).peakArea(p.Results.peakNumber);
        elseif contains(p.Results.mode,'intensity')
            disp('Intensity mode not yet programmed for XPS data')
        end
    
    elseif endsWith(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'h5')
        % if ~exist(p.Results.energyStart) | ~exist(p.Results.energyEnd)
        %     disp('The energy region used for normalisation appears to be empty. I quit now')
        %    error('The energy region used to normalise appears to be empty.')
        % end
        if contains(p.Results.mode,'intensity')
            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_col_scaling')
            end

            tmp_3(:,1) = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}...
                ./b_mean;

            DataTable(p.Results.tableSelection,:).spectrum_countsNorm{ii}  = tmp_3;
            clearvars tmp_3
        elseif contains(p.Results.mode,'area')
            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});
                % Area
                b_1 = find(xData >= p.Results.energyStart,1);
                b_2 = find(xData >= p.Results.energyEnd,1);

                b_area = trapz(xData(b_1:b_2),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_area = trapz(xData(b_2:b_1),yData(b_2:b_1));
            else
                disp('Warning: Issue with NEXAFS_col_scaling')
            end

            tmp_3(:,1) = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}...
                ./b_area;

             DataTable(p.Results.tableSelection,:).spectrum_countsNorm{ii} = tmp_3;
            clearvars tmp_3

        end
    end
end
% write log to table
DataTable.logScaleMode(p.Results.tableSelection,:)                  = cellstr(repmat(p.Results.mode,size(p.Results.DataTable(p.Results.tableSelection,:),1),1));
if endsWith(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'h5')
    DataTable.logScaleRange(p.Results.tableSelection,:)               = repmat([p.Results.energyStart p.Results.energyEnd],size(p.Results.DataTable(p.Results.tableSelection,:),1),1);
end
DataTable.logScaleDate(p.Results.tableSelection,:)                   = cellstr(repmat(datestr(now),size(p.Results.DataTable(p.Results.tableSelection,:),1),1));

end