function [a, b, c] = createCalculationTables(pRfb, pSoundsources, pTrackingParameter, pToaData, pSatData, pSoundvelocity, pLeapsecondsMatrix)
%CREATECALCULATIONMATRICES Summary of this function goes here
%   Detailed explanation goes here

%% Initialize required variables
%pSoundsources = artoa.controller.getSoundsourcesWithAppliedToa();
satPositions = [pSatData.lat_sat, pSatData.lon_sat];
satDates = artoa.convert.dmy2rd(pSatData.day_sat, pSatData.month_sat, pSatData.year_sat);
toaDates = [];
satToas = [];
satDistances = [];
floatDetails = pRfb.FLOAT;

%% Create table
results = struct();

%% Calculate SAT TOAs for every soundsource and get nearest neighbor measured toa
fnames = fieldnames(pSoundsources);
for i = 1:length(fnames)
    results.(fnames{i}) = table();
    % get start and end date of soundsource
    soundSourceBegin = artoa.convert.dmy2rd(pSoundsources.(fnames{i}).begemis(3), ...
        pSoundsources.(fnames{i}).begemis(2), ...
        pSoundsources.(fnames{i}).begemis(1));
    soundSourceEnd = artoa.convert.dmy2rd(pSoundsources.(fnames{i}).endemis(3), ...
        pSoundsources.(fnames{i}).endemis(2), ...
        pSoundsources.(fnames{i}).endemis(1));
    
    try
        soundvelocity = pSoundvelocity{fnames{i}, 1};
    catch
        soundvelocity = pSoundvelocity{1, 1};
    end
    
    [date, toa] = artoa.toa.predictFromGps( ...
        pRfb, ...
        pSoundsources.(fnames{i}), ...
        soundvelocity, ...
        pLeapsecondsMatrix ...
    );
    results.(fnames{i}).satDate = date;
    results.(fnames{i}).satToa = toa;
    results.(fnames{i}).daysSinceStart = ...
        date ...
        - artoa.convert.dmy2rd( ...
            pSoundsources.(fnames{i}).begemis(3), ...
            pSoundsources.(fnames{i}).begemis(2), ...
            pSoundsources.(fnames{i}).begemis(1) ...
        );
    tmpDistances = NaN(size(satPositions, 1), 1);
    % calculate distances
    for oDistance = 1:length(toa)
        if isnan(toa(oDistance)) | any(isnan(satPositions(oDistance, :)))
            %tmpDistances = [tmpDistances; NaN];
            continue;
        end
        tmpDistances(oDistance) = artoa.data.calculateGeodist( ...
                satPositions(oDistance, :), ...
                pSoundsources.(fnames{i}).position ...
        );
    end
    results.(fnames{i}).satDistances = tmpDistances;
    
    % set all sat dates to NaN that are not in soundsource mission range
    satDatesInRangeIndex = ...
        (results.(fnames{i}).satDate >= soundSourceBegin) ...
        & (results.(fnames{i}).satDate <= soundSourceEnd);
    results.(fnames{i}).satDate(~satDatesInRangeIndex) = NaN;
    results.(fnames{i}).satToa(~satDatesInRangeIndex) = NaN;
    results.(fnames{i}).daysSinceStart(~satDatesInRangeIndex) = NaN;
    results.(fnames{i}).satDistances(~satDatesInRangeIndex) = NaN;
    
    % find nearest neighbor toa
    tmpMeasuredToas = NaN(size(results.(fnames{i}).satDate));
    tmpIndexBelongingToSoso = strcmp(pToaData.soundSource, fnames{i});
    if ~any(tmpIndexBelongingToSoso) % no toa available so save NaN
        results.(fnames{i}).measuredToa = tmpMeasuredToas;
        continue;
    end
    tmpIndexBelongingToSoso = tmpIndexBelongingToSoso & (pToaData.status == 1);
    tmpToaDate = pToaData.toaDate(tmpIndexBelongingToSoso);
    tmpToa = pToaData.toa(tmpIndexBelongingToSoso);
    
    % interpolate toa data
    if ~strcmp(lower(pTrackingParameter.interpolationMethodString), 'none')
        [ ...
            tmpToaDate, ...
            tmpToa, ...
            ~ ...
        ] = artoa.data.interpolateRafosData( ...
            tmpToaDate, ...
            tmpToa, ...
            pTrackingParameter.interpolationInterval, ...
            pTrackingParameter.gapSize, ...
            lower(pTrackingParameter.interpolationMethodString) ...
        );
    end
    
    results.(fnames{i}).measuredToa = NaN(size(results.(fnames{i}).satDate));
    % find measured toa dates that belong to a sat date
    [values, satIndex, toaIndex] = intersect(results.(fnames{i}).satDate, ceil(tmpToaDate));
    results.(fnames{i}).measuredToa(satIndex) = tmpToa(toaIndex);

    % find nearest neighbor
%     for oSatDates = 1:length(results.(fnames{i}).satDate)
%         tmpSatDate = results.(fnames{i}).satDate(oSatDates);
%         neighborDate = artoa.data.findNearestNeighbor( ...
%             tmpToaDate, ...
%             tmpSatDate ...
%         );
%         if abs(neighborDate - tmpSatDate) <= pTrackingParameter.gapSize
%             tmpMeasuredToas(oSatDates) = tmpToa(tmpToaDate == neighborDate);
%         end
%     end
%     results.(fnames{i}).measuredToa = tmpMeasuredToas;
    clear tmpDistances;
end

%% Construct matrices A and B

rowCount = length(satDates) * length(fnames);

aSoundsources = table();
%aCore = zeros(rowCount, 2 * length(fnames));
%b = zeros(rowCount, 1);
b = table(zeros(rowCount, 1), 'VariableNames', {'toa'});
c = table(zeros(rowCount, 1), 'VariableNames', {'ReferenceTime'});
Distances = zeros(rowCount, 1);
daysSinceFloatStart = NaN(rowCount, 1);

% disable warning because table produces a lot when adding the soundsources
warning('off');

for i = 1:length(fnames)
    aSoundsources{:, (fnames{i})} = 0;
    rowIndices = ((i - 1) * length(satDates) + 1):i * length(satDates);
    %startColIndex = (2 * (i - 1)) + 1;
    Distances(rowIndices, 1) = results.(fnames{i}).satDistances;
    %aCore(rowIndices, startColIndex) = 1;
    %aCore(rowIndices, startColIndex + 1) = results.(fnames{i}).daysSinceStart;
    
    aSoundsources{rowIndices, (fnames{i})} = results.(fnames{i}).daysSinceStart;
    b{rowIndices, 'toa'} = results.(fnames{i}).measuredToa;
    c{rowIndices, 'ReferenceTime'} = ...
        (pSoundsources.(fnames{i}).reftime(1) * 3600 + pSoundsources.(fnames{i}).reftime(2) * 60) ...
        - (floatDetails.phasereftime(1) * 3600 + floatDetails.phasereftime(2) * 60 + floatDetails.windowstart * 60);
    
    daysSinceFloatStart(rowIndices, 1) = ...
        results.(fnames{i}).satDate ...
        - artoa.convert.dmy2rd( ...
            floatDetails.launchtime(3), ...
            floatDetails.launchtime(2), ...
            floatDetails.launchtime(1) ...
        );
    % construct b
%     foundToas = findCorrespondingToas(fnames{i}, results.(fnames{i}).satDate);
end

a = addvars(aSoundsources, Distances, 'Before', 1); % add distances to soundsource
a.Float = daysSinceFloatStart; % append float days since start
% a = [ ...
%     Distances, aCore{:, :}, ones(size(daysSinceFloatStart)), daysSinceFloatStart ...
% ];

% remove all NaN from matrix
indicesToUse = all(~isnan(a{:, :}), 2) & all(aSoundsources{:, :} >= 0, 2) & ~isnan(b{:, :});

a = a(indicesToUse, :);
b = b(indicesToUse, :);
c = c(indicesToUse, :);

% reenable warning
warning('on');

end