function dat = Hbb_ReadFile(dirSpec,varargin)
%Hbb_ReadFile Read in Hyper-bb instrument data
%   DAT = readHyperBb('FILENAME') will read Hyper-bb instrument data
%   (ASCII-formatted) from the file FILENAME into a table DAT.
%
%   If 'FILENAME' is a path or includes wildcards (acceptable to dir
%   command) then multiple files are read and DAT is cell array of tables.
%   See also the 'Collect' parameter that can be used to combine all of the
%   files into a single table.
%
%   DAT = readHyperBb('FILENAME', 'PARAM',VALUE) allows specification of
%   parameter/value pairs separated by commas. Parameter/value combinations
%   are provided below.
%
%   Parameter/value combinations:
%
%   Parameter               Value (type) Description
%   ---------               ------------------------
%   IncludeRawData          (boolean) If true all raw data contained in the
%                           input filename is read into the table,
%                           including all LED on/of ADC values and stdevs.
%                           Default is false.
%
%   RemoveMultiplePmtGains  (boolean) If true then scans containing
%                           measurements with different PMT gain values
%                           (i.e, due to auto gain adjustment) are removed
%                           from the output table. Default is false.
%
%   Collect                 (boolean) If true, then all of the data files
%                           read in are combined into a single table.
%

% WHS 20191204 initial version
% WHS 20191216 changed multiple pmt gain to warning
% WHS 20200909 added ScatX and GainX to raw dat table
p = inputParser;
defaultCollect = false;
defaultRemoveMultiplePmtGains = false;
defaultIncludeRawData = false;

addRequired(p, 'dirSpec');
addParameter(p, 'Collect',defaultCollect);
addParameter(p, 'RemoveMultiplePmtGains',defaultRemoveMultiplePmtGains);
addParameter(p, 'IncludeRawData',defaultIncludeRawData);

parse(p,dirSpec,varargin{:});
   
   
dd = dir(dirSpec); 
datfiles = fullfile({dd.folder}',{dd.name}');

%                                                                              PMT   PMT   PMT   REF    REF    PMT   PMT    REF   REF   LOW   LOW   HIGH   HIGH   LOW    LOW   HIGH  HIGH
%                                           STEP           LED                 NET   ON    ON    ON     ON     OFF   OFF    OFF   OFF   ON    ON    ON     ON     OFF    OFF   OFF   OFF   LED   Water
%                              Date  TIME   POS        WL  PWR       GAIN      SIG   AVG   STD   AVG    STD    AVG   STD    AVG   STD   AVG   STD   AVG    STD    AVG    STD   AVG   STD   Temp  Temp   Depth

VarNames30={'ScanIdx','DataIdx','Date','Time','StepPos','wl','LedPwr','PmtGain','NetSig1','SigOn1','SigOn1Std','RefOn','RefOnStd','SigOff1','SigOff1Std','RefOff','RefOffStd', ...
  'SigOn2','SigOn2Std','SigOn3','SigOn3Std','SigOff2','SigOff2Std','SigOff3','SigOff3Std','LedTemp','WaterTemp','Depth','Debug1','zDistance'};


for kf = 1:length(datfiles)

    disp(['Reading ' num2str(kf) '/' num2str(length(datfiles)) '   ' datfiles{kf}])
    dat{kf} = readtable(datfiles{kf},'ReadVariableNames',0, ...
                    'HeaderLines',0, 'Delimiter',{'\t',' '}, 'MultipleDelimsAsOne',true);

    if size(dat{kf},2)==30 % normal data file (last column is zDist if cal)
        dat{kf}.Properties.VariableNames = VarNames30;
    else
        error('Improper number of columns in data file')
    end
    
    for ks = unique(dat{kf}.ScanIdx')
        scanSel = dat{kf}.ScanIdx==ks;
        if length(unique(dat{kf}.PmtGain(scanSel)))~=1
            if p.Results.RemoveMultiplePmtGains
                warning(['  Scan ' num2str(ks) ' has multiple PmtGain, REMOVING'])
            	dat{kf}(scanSel,:) = [];
            else
                warning(['  Scan ' num2str(ks) ' has multiple PmtGain'])
            end
        end
    end
    
    
    
    satLevel = 4000;
    dat{kf}.SigOn3(dat{kf}.SigOn3>satLevel) = NaN;
    dat{kf}.SigOn2(dat{kf}.SigOn2>satLevel) = NaN;
    dat{kf}.SigOn1(dat{kf}.SigOn1>satLevel) = NaN;
    dat{kf}.SigOff3(dat{kf}.SigOff3>satLevel) = NaN;
    dat{kf}.SigOff2(dat{kf}.SigOff2>satLevel) = NaN;
    dat{kf}.SigOff1(dat{kf}.SigOff1>satLevel) = NaN;
    % dat{kf}.SigOn3(dat{kf}.SigOn3>satLevel | dat{kf}.SigOn3<0) = NaN;
    % dat{kf}.SigOn2(dat{kf}.SigOn2>satLevel | dat{kf}.SigOn2<0) = NaN;
    % dat{kf}.SigOn1(dat{kf}.SigOn1>satLevel | dat{kf}.SigOn1<0) = NaN;
    
    % figure, hold on, plot(dat{kf}.SigOn3,'.-')    
    % plot(dat{kf}.SigOn1,'.-')
    % plot(dat{kf}.SigOn1,'.-')

    % Calculate the net high and low gain, net ref
    dat{kf} = addvars(dat{kf},dat{kf}.SigOn3-dat{kf}.SigOff3, 'After','NetSig1', 'NewVariableNames', 'NetSig3');
    dat{kf} = addvars(dat{kf},dat{kf}.SigOn2-dat{kf}.SigOff2, 'After','NetSig1', 'NewVariableNames', 'NetSig2');
    dat{kf} = addvars(dat{kf},dat{kf}.RefOn-dat{kf}.RefOff, 'Before','NetSig1', 'NewVariableNames', 'NetRef');
    
    
    dat{kf} = addvars(dat{kf},dat{kf}.NetSig1./dat{kf}.NetRef, 'Before','NetSig1', 'NewVariableNames', 'Scat1');
    dat{kf} = addvars(dat{kf},dat{kf}.NetSig2./dat{kf}.NetRef, 'Before','NetSig1', 'NewVariableNames', 'Scat2');
    dat{kf} = addvars(dat{kf},dat{kf}.NetSig3./dat{kf}.NetRef, 'Before','NetSig1', 'NewVariableNames', 'Scat3');
    
    % ScatX is the highest raw net signal from the front end that is not
    % saturated. GainX indicates the front end channel (1,2,3) assigned to
    % ScatX.
    % Start by assigning ScatX and GainX to highest front end gain level (3)
    ScatX = dat{kf}.Scat3;
    GainX = repmat(3,size(ScatX));
    if any(isnan(ScatX)) % If high gain is saturated, replace with low gain (2)
        GainX(isnan(ScatX)) = 2;
        ScatX(isnan(ScatX)) = dat{kf}.Scat2(isnan(ScatX)); 
    end
    if any(isnan(ScatX)) % it low gain is still saturated, replace with raw pmt no gain (1)
        GainX(isnan(ScatX)) = 1;
        ScatX(isnan(ScatX)) = dat{kf}.Scat1(isnan(ScatX)); 
    end
    
    dat{kf} = addvars(dat{kf},ScatX,GainX, 'After','Scat3', 'NewVariableNames', {'ScatX','GainX'});
    
    
    dat{kf} = addvars(dat{kf}, datetime(dat{kf}.Date, 'InputFormat','yyyy/MM/dd') + dat{kf}.Time, ...
                'NewVariableNames','Timestamp', 'Before','Date');
    dat{kf} = removevars(dat{kf}, {'Date','Time'});
    
    % dat{kf} = mergevars(dat{kf}, {'StepPos','PmtOn','PmtOnStd','RefOn','RefOnStd','PmtOff','PmtOffStd','RefOff','RefOffStd', ...
    %       'LowOn','LowOnStd','HighOn','HighOnStd','LowOff','LowOffStd','HighOff','HighOffStd','Debug1','Debug2'}, 'NewVariableName','RawData');
    
%     if dat{kf}.zDistance(1) == -1
%         dat{kf} = removevars(dat{kf}, 'zDistance');
%     end
    
    
    
    % dat{kf} = mergevars(dat{kf}, {'StepPos','NetRef','RefOn','RefOnStd','RefOff','RefOffStd', 'NetSig1','SigOn1','SigOn1Std','SigOff1','SigOff1Std', ...
    %         'NetSig2','SigOn2','SigOn2Std','SigOff2','SigOff2Std', 'NetSig3','SigOn3','SigOn3Std','SigOff3','SigOff3Std','Debug1'}, 'NewVariableName','RawData');
    % dat{kf} = movevars(dat{kf}, 'RawData', 'After','Depth');

    
    if ~p.Results.IncludeRawData 
        % dat{kf} = removevars(dat{kf}, 'RawData');
        dat{kf} = removevars(dat{kf}, {'StepPos','NetRef','RefOn','RefOnStd','RefOff','RefOffStd', 'NetSig1','SigOn1','SigOn1Std','SigOff1','SigOff1Std', ...
                'NetSig2','SigOn2','SigOn2Std','SigOff2','SigOff2Std', 'NetSig3','SigOn3','SigOn3Std','SigOff3','SigOff3Std','Debug1'});
    end
      
end

if length(dat)==1
    dat = dat{1};
elseif length(dat)>1 && p.Results.Collect==true

    disp(' ')
    disp('Collecting data into single table...')
    
    datArr = dat;
    dat = datArr{1};

    dat.Properties.VariableNames{'ScanIdx'} = 'FileScanIdx';
    
    dat = addvars(dat, repmat(1,height(dat),1), ...
                'NewVariableNames','FileIdx', 'Before',1); %#ok<REPMAT>
            
    dat = addvars(dat, dat.FileScanIdx, ...
                'NewVariableNames','ScanIdx', 'Before',1);

    % lastScanIdx = max(dat.FileScanIDx)
    
    for kf = 2:length(datfiles)

        tmpDat = datArr{kf};

        tmpDat.Properties.VariableNames{'ScanIdx'} = 'FileScanIdx';

        tmpDat = addvars(tmpDat, repmat(kf,height(tmpDat),1), ...
                    'NewVariableNames','FileIdx', 'Before',1);
        tmpDat = addvars(tmpDat, max(dat.ScanIdx) + tmpDat.FileScanIdx, ...
                    'NewVariableNames','ScanIdx', 'Before',1);
        
        dat = [dat; tmpDat];        
                
        % ScanIdx = (lastScanIdx+1):(lastScanIdx+1+height(datArr{kf}))
        
    end    
    
    
end


