% Hbb_Process compute corrected scattering and Beta, given a data file,
%   plaque calibration file, and temperature calibration file.
% 
% Usage:
%   dat = Hbb_Process(Data_FilePath, Plaque_Cal_FilePath, Temperature_Cal_FilePath)
%
% Inputs:
%   Inputs are the file path to the Hyper-bb data file you wish to process,
%   the calibration file (.mat) generated using the Hbb_ReadCal, and the
%   temperature calibration (.mat) provided by Sequoia Scientific.
%
% Optional Parameters:
%
%   Parameter               Value (type) Description
%   ---------               ------------------------
%   Reshape                 (boolean) If true the data will be reshaped so
%                           each row contains an entire wavelength scan.
%                           Some of the metadata will be averaged for each
%                           scan, such as temperature and depth. Much of
%                           the raw signal data are excluded from the 
%                           reshaped table (such as sig on/off signals).
%
% Output:
%   Hbb_ReadCal outputs a structure containing both the raw and processed data. 
%   If the 'Reshape' option is chosen, only processed data is returned.
%
%   ScatCor1,2,3 are the low, medium, and high gain channel net scattering 
%   signals after applying the gain factors Gfe12 and Gfe23, as well as 
%   correcting for the PMT gain setting, Gpmt. ScatTempCor is after
%   applying the temperature correction to ScatCor.
%
%   ScatTempCorX is the highest gain channel scattering that is not saturated.
%   
%   Beta_u is the VSF calculated from ScatCorX using the mu calibration
%   factors determined from reflective plaque measurements.
%
% Sequoia Scientific, Inc. - 11/30/2022

function dat = Hbb_Process(datafile, calfile_plaque, calfile_temp, varargin)

p = inputParser;
addRequired(p, 'datafile');
addRequired(p, 'calfile_plaque');
addRequired(p, 'calfile_temp');
addParameter(p, 'Reshape',false);
parse(p,datafile,calfile_plaque,calfile_temp,varargin{:});

% read data file
dat = Hbb_ReadFile(datafile, 'IncludeRawData',true);

% read cal files
load(calfile_plaque);
if ~all(isfield(cal, {'gain12', 'gain23', 'pmtGamma', 'pmtRefGain', 'H', 'rho', 'muWavelengths', 'muFactors'}))
    error('Input cal struct does not contain required fields.')
end

%load temperature cal file
load(calfile_temp);

if istable(dat)
    if sum(ismember(dat.Properties.VariableNames,{'Scat1','Scat2','Scat3','PmtGain'})) == 4
        
        dat = processDataTable(dat, cal, cal_temp);
        
    else
        warning(['dat{' num2str(kd) '} does not contain proper data.'])
    end
else
    warning(['dat{' num2str(kd) '} is not a data table.'])
end

if p.Results.Reshape
    [dat, ~] = reshapeHyperBb(dat);
end

% save calibration information in table properties
dat = addprop(dat,{'PlaqueCal','TemperatureCal'}, ...
              {'table','table'});
          
dat.Properties.CustomProperties.PlaqueCal = cal;
dat.Properties.CustomProperties.TemperatureCal = cal_temp;

end 

function dat = processDataTable(dat, cal_plaque, cal_temp)

    % Interpolate wl and pmt to find dark offset
    darkOffset_scat1 = interp2(cal_plaque.darkCalPmtGain,cal_plaque.darkCalWavelength,cal_plaque.darkCalScat1,dat.PmtGain,dat.wl, 'linear');
    darkOffset_scat2 = interp2(cal_plaque.darkCalPmtGain,cal_plaque.darkCalWavelength,cal_plaque.darkCalScat2,dat.PmtGain,dat.wl, 'linear');
    darkOffset_scat3 = interp2(cal_plaque.darkCalPmtGain,cal_plaque.darkCalWavelength,cal_plaque.darkCalScat3,dat.PmtGain,dat.wl, 'linear');

    % Subtract Dark Offset
    scat1_darkRemoved = dat.Scat1 - darkOffset_scat1;
    scat2_darkRemoved = dat.Scat2 - darkOffset_scat2;
    scat3_darkRemoved = dat.Scat3 - darkOffset_scat3;

    % Apply PMT and front end gain factors
    Gpmt = (dat.PmtGain./cal_plaque.pmtRefGain).^cal_plaque.pmtGamma;
    scat1_gainCorrected = scat1_darkRemoved.*cal_plaque.gain12.*cal_plaque.gain23.*Gpmt;
    scat2_gainCorrected = scat2_darkRemoved.*cal_plaque.gain23.*Gpmt;
    scat3_gainCorrected = scat3_darkRemoved.*Gpmt;
   
    % Add data to output table
    dat = addvars(dat,scat1_gainCorrected,'After','Scat3', 'NewVariableNames', 'ScatCor1');
    dat = addvars(dat,scat2_gainCorrected,'After','ScatCor1', 'NewVariableNames', 'ScatCor2');
    dat = addvars(dat,scat3_gainCorrected,'After','ScatCor2', 'NewVariableNames', 'ScatCor3');

    % Apply temperature correction
     tempCoeff = GetTemperatureCoefficients(cal_temp,dat.wl,dat.LedTemp);
     scat1_tempCorrected = dat.ScatCor1.*tempCoeff;
     scat2_tempCorrected = dat.ScatCor2.*tempCoeff;
     scat3_tempCorrected = dat.ScatCor3.*tempCoeff;
    
    % Add data to output table
    dat = addvars(dat,scat1_tempCorrected, 'After','ScatCor3', 'NewVariableNames', 'ScatTempCor1');
    dat = addvars(dat,scat2_tempCorrected, 'After','ScatTempCor1', 'NewVariableNames', 'ScatTempCor2');
    dat = addvars(dat,scat3_tempCorrected, 'After','ScatTempCor2', 'NewVariableNames', 'ScatTempCor3');
    dat = addvars(dat,tempCoeff, 'After','ScatTempCor3', 'NewVariableNames', 'TempCorrCoeff'); % save the calculated temp coeff.

    
    % Select highest non-saturated gain channel
    ScatTempCorX = dat.ScatTempCor3;
    ScatTempCorX(isnan(ScatTempCorX)) = dat.ScatTempCor2(isnan(ScatTempCorX)); % If high gain is saturated, replace with low gain
    ScatTempCorX(isnan(ScatTempCorX)) = dat.ScatTempCor1(isnan(ScatTempCorX)); % it low gain is saturated, replace with raw pmt
    dat = addvars(dat,ScatTempCorX, 'After','ScatTempCor3', 'NewVariableNames', 'ScatTempCorX');

    % Temperature correct mu calibration
     tempCoeff_mu = GetTemperatureCoefficients(cal_temp,cal_plaque.muWavelengths,cal_plaque.muLedTemp);
     muFactors_tempCorrected = cal_plaque.muFactors.*tempCoeff_mu;

    % Calculate Beta
    wlDat = sort(unique(dat.wl))';
    muFactors = interp1(cal_plaque.muWavelengths,muFactors_tempCorrected, wlDat, 'pchip');
    dat = addvars(dat,NaN(height(dat),1), 'After','ScatTempCorX', 'NewVariableNames', 'Beta_u'); % preallocate array
    for kw = 1:length(wlDat)
        dat.Beta_u(dat.wl==wlDat(kw)) = dat.ScatTempCorX(dat.wl==wlDat(kw)).*muFactors(kw);
    end
           
end

function tempCoeff = GetTemperatureCoefficients(cal_temp,wavelength,temperature) 

% Generate temperature correction grid
LEDTempRange = [min(temperature):0.1:max(temperature)+0.1]; % need to make sure the max value is included
TempCorrGrid = NaN(length(cal_temp.wl),length(LEDTempRange));
for n = 1:length(cal_temp.wl)
   TempCorrGrid(n,:) = polyval(cal_temp.coeff(n,:),LEDTempRange);
end

% 2D interpolate (wavelength x temperature) to find temperature
% correction factor
tempCoeff = interp2(LEDTempRange, cal_temp.wl, TempCorrGrid, temperature, wavelength, 'linear');

end


function [rdat, wl] = reshapeHyperBb(dat)
%RESHAPEHYPERBB Reformat Hyper-bb data into spectral-column data
%   DAT2 = reshapeHyperBb(DAT1) will reshape Hyper-bb instrument data
%   contained in table DAT1 (i.e., output of readHyperBb) and convert to an
%   output format where each table row is a multi-spectral measurement.

wl = sort(unique(dat.wl))';
scanIdx = sort(unique(dat.ScanIdx));

ns = length(scanIdx);
nwl = length(wl);

%            ScanIdx    Timestamp                      LedPwr     PmtGain
rdat = table(NaN(ns,1), repmat(dat.Timestamp(1),ns,1), NaN(ns,1), NaN(ns,1), ...
            NaN(ns,nwl), NaN(ns,nwl), ...
            NaN(ns,1), NaN(ns,1), NaN(ns,1), ...% LedTemp    WaterTemp    Depth
        'VariableNames',{'ScanIdx','Timestamp','LedPwr','PmtGain','wl','Beta_u','LedTemp','WaterTemp','Depth'});
    
for ks = 1:length(scanIdx)
    
    thisScan = dat.ScanIdx==scanIdx(ks);
    
    rdat.ScanIdx(ks) = scanIdx(ks);
    rdat.Timestamp(ks) = min(dat.Timestamp(thisScan)); % start of scan
    if length(unique(dat.LedPwr(thisScan)))~=1
        rdat.LedPwr(ks) = NaN; % multiple values indicated gain/power switch during scan
    else
        rdat.LedPwr(ks) = unique(dat.LedPwr(thisScan));
    end
    if length(unique(dat.PmtGain(thisScan)))~=1
        rdat.PmtGain(ks) = NaN; % multiple values indicated gain/power switch during scan
    else
        rdat.PmtGain(ks) = unique(dat.PmtGain(thisScan));
    end
    
    rdat.LedTemp(ks) = nanmean(dat.LedTemp(thisScan)); 
    rdat.WaterTemp(ks) = nanmean(dat.WaterTemp(thisScan)); 
    rdat.Depth(ks) = nanmean(dat.Depth(thisScan)); 
    
    % [dat.wl(thisScan) dat.(p.Results.Variable)(thisScan,:)]
    
    for kw = 1:nwl
    	thisWavelength = dat.wl==wl(kw);
        if any(thisScan&thisWavelength)
             rdat.wl(ks,kw) = wl(kw);
            rdat.Beta_u(ks,kw) = dat.Beta_u(thisScan&thisWavelength);
        end
    end
    
end
end
