calculateCombinationSegment.m

function [ ...
    segmentPositions, ...
    segmentDates, ...
    segmentClockError, ...
    segmentResiduals ...
    ] = ...
    calculateCombinationSegment( ...
        pCorrectedData, ...
        pCombinationDetails, ...
        pFloatReferenceTime, ...
        pSoundVelocity, ...
        pBackwardTracking, ...
        pStartEndPosition ...
    )
%CALCULATECOMBINATIONSEGMENT Combination segment calculation.
%   Using the input data, the function calculates the segment positions,
%   its corresponding dates as well as the clock error (in case of hyperbolic calculation)
%
%   Parameters:
%       pCorrectedData (struct):
%           Data structured by soundsource. The toa is already corrected
%           (e.g. doppler correction etc.)
%       pCombinationDetails (table):
%           The combination table that is produced by the track parameter
%           window.
%       pFloatReferenceTime (double):
%           Vector containing the float reference time, in the
%           format [hour, minute]
%       pSoundVelocity (double):
%           Sound velocity next to the sound sources that are used.
%           If hyperbolic is being used, it needs to be a vector containing
%           at least as many entries as the count of sound sources being
%           used.
%           For all other methods, only one value needs to be set.
%           Unit: [m/s]
%       pBackwardTracking (bool):
%           If true, the data will be prepared for backward tracking
%       pStartEndPosition (double):
%           Contains the start and end position of the segment as a matrix
%           consisting of row vectors. It is passed to the tracking methods.
%
%   Returns:
%       segmentPositions (double):
%           The calculated list of segment positions as a column vector.
%       segmentDates (double):
%           The corresponding rafos dates of the segment positions.
%       clockError (double):
%           The estimated clock error.
%       residuals (table):
%           The calculated residuals corresponding to segmentDates for
%           every soundsource.

%% Get sound sources
soundsourceNames = unique(strsplit(strtrim(pCombinationDetails.soundsources{1}), ' '));

%% Get tracking method
trackingMethod = pCombinationDetails.trackingMethod{:};

%% Get start and end point
segmentStart = pCombinationDetails{1, 1};
segmentEnd = pCombinationDetails{1, 2};

%% Filter toa data
intersectedToaDates = pCorrectedData.(soundsourceNames{1}).toaDate;
soundsourcePositions = [pCorrectedData.(soundsourceNames{1}).position];
pluginSoundsourcePositions = struct();
pluginSoundsourcePositions.(soundsourceNames{1}) = pCorrectedData.(soundsourceNames{1}).position;
for i = 2:length(soundsourceNames)
    [intersectedToaDates] = intersect(intersectedToaDates, pCorrectedData.(soundsourceNames{i}).toaDate);
    soundsourcePositions = [soundsourcePositions; pCorrectedData.(soundsourceNames{i}).position];
    pluginSoundsourcePositions.(soundsourceNames{i}) = pCorrectedData.(soundsourceNames{i}).position;
end

%% Remove all dates that are out of bounds
intersectedToaDates = intersectedToaDates( ...
    intersectedToaDates >= segmentStart ...
    & intersectedToaDates <= segmentEnd ...
);


%% Create a table for every soundsource (required by tracking plugin system)
pluginTables = struct();
for i = 1:length(soundsourceNames)
    pluginTables.(soundsourceNames{i}) = table();
end

%% Prepare data

lenIntersectedToaDates = length(intersectedToaDates);
distanceToSoundsources = cell(1, lenIntersectedToaDates);

warning('off'); % disable warnings because creating the tables are not consistent

for oDates = 1:lenIntersectedToaDates
    currentDateValue = intersectedToaDates(oDates);
    distanceToSoundsources{oDates} = [];
    for i = 1:length(soundsourceNames)
        currentName = soundsourceNames{i};
        currentData = pCorrectedData.(currentName);
        selection = currentData.toaDate == currentDateValue;
        filteredValue = currentData.toa(selection);

        windowStartTime = ...
            pFloatReferenceTime(1) * 3600 ...
            + pFloatReferenceTime(2) * 60 ...
            - currentData.soundsourceReferenceTime(1) * 3600 ...
            - currentData.soundsourceReferenceTime(2) * 60;
        leapseconds = currentData.leapseconds;
        floatWindowStart = currentData.floatWindowStart * 60;
    
        relativeToa = filteredValue(1) ...
            + windowStartTime ...
            + leapseconds ...
            + floatWindowStart;
        
        % calculate distance to the source
        distanceToSoundsources{oDates} = [ ...
            distanceToSoundsources{oDates}; ...
            relativeToa * (pSoundVelocity(i)/1000)
        ];
    
        % store in pluginTables
        pluginTables.(currentName).date{oDates} = currentDateValue;
        pluginTables.(currentName).toa{oDates} = relativeToa;
        pluginTables.(currentName).distance{oDates} = distanceToSoundsources{oDates}(i);

    end

end
warning('on');

%% Backward tracking
% if enabled, everything needs to be upside down for the calculation
if pBackwardTracking
    for i = 1:length(soundsourceNames)
        pluginTables.(soundsourceNames{i}) = flipud(pluginTables.(soundsourceNames{i}));
    end
end

%% Call the corresponding tracking method from the plugin folder
[segmentPositions, segmentClockError] = artoa.plugins.tracking.callTrackingMethod( ...
    trackingMethod, ...
    pluginTables, ...
    pluginSoundsourcePositions, ...
    pCombinationDetails, ...
    pStartEndPosition ...
    );

%% Restore backward tracking changes
% if enabled, the calculated positions need to be flipped upside down to be
% able to display them correctly
if pBackwardTracking
    segmentPositions = flipud(segmentPositions);
    segmentClockError = flipud(segmentClockError);
end

%% Store toa dates
segmentDates = intersectedToaDates;

%% Calculate residuals

% initalize variables
segmentResiduals = struct();
for i = 1:length(soundsourceNames)
    segmentResiduals.(soundsourceNames{i}) = [];
end

% calculate residuals for every timestep
for oDates = 1:length(intersectedToaDates)
    [tmpResiduals, ~] = artoa.data.calculateResiduals( ...
        segmentPositions(1, :), ...
        distanceToSoundsources{oDates}, ...
        soundsourcePositions, ...
        pSoundVelocity ...
    );
    % split up for storing it as table
    for i = 1:length(soundsourceNames)
        segmentResiduals.(soundsourceNames{i}) = [ ...
            segmentResiduals.(soundsourceNames{i}); ...
            tmpResiduals(i) ...
        ];
    end
end


end