Extended float package.

lib/+artoa/+float/calculateTrajectory.m

+function [trajectory] = calculateTrajectory(pFloatDetails, pToaData, pSoundsources, pTrackingParameter, pLeapsecondsMatrix)
+%CALCULATETRAJECTORY Summary of this function goes here
+%   Detailed explanation goes here
+
+
+%% Get required data
+soundsourceCombinations = pTrackingParameter.soundsourceCombinations;
+uniqueToaDates = unique(pToaData.toaDate);
+
+%% Get all sound sources that are involved
+involvedSoundsourceNames = ...
+    unique(flatten( ...
+        cellfun( ...
+            @strsplit, soundsourceCombinations.soundsources, 'UniformOutput', false ...
+        ) ...
+    ));
+involvedSoundsources = artoa.data.extractSoundsourcesFromStruct( ...
+    involvedSoundsourceNames, ...
+    pSoundsources ...
+);
+
+%% Add float drift to TOA
+pToaData = artoa.data.addFloatDrift( ...
+        pToaData, ...
+        pTrackingParameter.floatOffsetBegin, ...
+        pTrackingParameter.floatOffsetEnd ...
+);
+
+%% Prepare data for every sound source
+preparedData = struct();
+for i = 1:length(involvedSoundsources)
+    currentSoundsource = involvedSoundsources{i};
+    preparedData.(currentSoundsource.sourcename) = struct();
+    % get toadata that is applied to soundsource and not deleted
+    currentToaSelection = ...
+        strcmp(pToaData.soundSource, {currentSoundsource.sourcename}) ...
+        & pToaData.status ~= 2;
+    
+    if ~any(currentToaSelection)
+        trajectory = false;
+        return
+    end
+    
+    % get toa and date of current sound source
+    currentToa = pToaData.toa(currentToaSelection);
+    currentToaDate = pToaData.toaDate(currentToaSelection);
+    
+    
+    % get rafos start and enddate of soundsource
+    rafosBeginEmission = artoa.convert.dmy2rd( ...
+        currentSoundsource.begemis(3), ...
+        currentSoundsource.begemis(2), ...
+        currentSoundsource.begemis(1) ...
+    );
+    rafosEndEmission = artoa.convert.dmy2rd( ...
+        currentSoundsource.endemis(3), ...
+        currentSoundsource.endemis(2), ...
+        currentSoundsource.endemis(1) ...
+    );
+
+    % get rafos date of float based on cycle (expecting only 1 cycle)
+    rafosFloatDateCycleBegin = artoa.convert.dmy2rd( ...
+        pFloatDetails.cycle(6), ...
+        pFloatDetails.cycle(5), ...
+        pFloatDetails.cycle(4) ...
+    );
+    rafosFloatDateCycleEnd = artoa.convert.dmy2rd( ...
+        pFloatDetails.cycle(13), ...
+        pFloatDetails.cycle(12), ...
+        pFloatDetails.cycle(11) ...
+    );
+
+    % calculate soundsource offsets and dates
+    [offsets, dates] = artoa.data.calculateOffsetsAndDates( ...
+        [rafosFloatDateCycleBegin, rafosFloatDateCycleEnd], ...
+        [rafosBeginEmission, rafosEndEmission], ...
+        currentSoundsource.offset, ...
+        currentSoundsource.drift, ...
+        currentSoundsource.schedule ...
+    );
+
+    % check overlap between float and sound source deployment and add sound
+    % source offset
+    for oDates = 1:length(dates)
+        kx = find(abs(currentToaDate - dates(oDates)) < 1/2);
+        currentToa(kx) = currentToa(kx) - offsets(oDates);
+    end
+    
+    % subtract offset (corrective value) of sound source
+    tmpFnames = fieldnames(pSoundsources);
+    soundsourceOffsetCorrectiveValue = ...
+        pTrackingParameter.soundsourceOffsets.(currentSoundsource.sourcename);
+    currentToa = currentToa - soundsourceOffsetCorrectiveValue(1);
+        
+    % interpolate toa data
+    [ ...
+        preparedData.(currentSoundsource.sourcename).toaDate, ...
+        preparedData.(currentSoundsource.sourcename).toa, ...
+        preparedData.(currentSoundsource.sourcename).toaIndex ...
+    ] = artoa.data.interpolateRafosData( ...
+        currentToaDate, ...
+        currentToa, ...
+        pTrackingParameter.interpolationInterval, ...
+        pTrackingParameter.gapSize, ...
+        lower(pTrackingParameter.interpolationMethodString) ...
+    );
+    preparedData.(currentSoundsource.sourcename).toaIndex = preparedData.(currentSoundsource.sourcename).toaIndex';
+    
+    % combine interpolated and original data and sort them
+    preparedData.(currentSoundsource.sourcename).toaDate = [preparedData.(currentSoundsource.sourcename).toaDate; currentToaDate];
+    preparedData.(currentSoundsource.sourcename).toa = [preparedData.(currentSoundsource.sourcename).toa; currentToa];
+    
+    % Find the nearest neighbor rafos date
+    for oRafosDate = 1:length(preparedData.(currentSoundsource.sourcename).toaDate)
+        neighbor = artoa.data.findNearestNeighbor(uniqueToaDates, preparedData.(currentSoundsource.sourcename).toaDate(oRafosDate));
+        preparedData.(currentSoundsource.sourcename).toaDate(oRafosDate) = neighbor;
+        clear neighbor;
+    end
+    
+    % sort toa by date
+    [preparedData.(currentSoundsource.sourcename).toaDate, sortIndices] = sort(preparedData.(currentSoundsource.sourcename).toaDate);
+    preparedData.(currentSoundsource.sourcename).toa = preparedData.(currentSoundsource.sourcename).toa(sortIndices);
+    
+    % shift toa data
+    [preparedData.(currentSoundsource.sourcename).toaDate, preparedData.(currentSoundsource.sourcename).toa] = ...
+        artoa.data.simulshift( ...
+            preparedData.(currentSoundsource.sourcename).toaDate, ...
+            preparedData.(currentSoundsource.sourcename).toa, ...
+            pFloatDetails.schedule, ...
+            pFloatDetails.windowstart, ...
+            currentSoundsource.reftime, ...
+            pFloatDetails.phasereftime ...
+    );
+
+    % save position
+    preparedData.(currentSoundsource.sourcename).position = currentSoundsource.position;
+    % save float window start
+    preparedData.(currentSoundsource.sourcename).floatWindowStart = pFloatDetails.windowstart;
+    
+    
+    % save soundsource reference time
+    preparedData.(currentSoundsource.sourcename).soundsourceReferenceTime = currentSoundsource.reftime;
+    
+    
+    % get leapseconds for soundsource
+    preparedData.(currentSoundsource.sourcename).leapseconds = artoa.data.calculateLeapseconds( ...
+        artoa.convert.dmy2rd(pFloatDetails.launchtime(3), pFloatDetails.launchtime(2), pFloatDetails.launchtime(1)), ...
+        artoa.convert.dmy2rd(currentSoundsource.begemis(3), currentSoundsource.begemis(2), currentSoundsource.begemis(1)), ...
+        pLeapsecondsMatrix ...
+    );
+    if isempty(preparedData.(currentSoundsource.sourcename).leapseconds)
+        preparedData.(currentSoundsource.sourcename).leapseconds = 0;
+    end
+    
+    if ~pTrackingParameter.dopplerCorrection
+        continue
+    end
+    
+    % apply doppler correction if selected
+    [preparedData.(currentSoundsource.sourcename).toaDate, preparedData.(currentSoundsource.sourcename).toa] = ...
+        artoa.data.doppcorr( ...
+            preparedData.(currentSoundsource.sourcename).toaDate, preparedData.(currentSoundsource.sourcename).toa, pFloatDetails.schedule ...
+    );
+    
+    %hms = artoa.convert.rd2hms(preparedData.(currentSoundsource.sourcename).toaDate);
+    %preparedData.(currentSoundsource.sourcename).toa = preparedData.(currentSoundsource.sourcename).toa - 60*hms(2) - hms(3);
+    %preparedData.(currentSoundsource.sourcename).toa = preparedData.(currentSoundsource.sourcename).toa + pFloatDetails.windowstart * 60;
+    %preparedData.(currentSoundsource.sourcename).toa = preparedData.(currentSoundsource.sourcename).toa - 3600;
+end
+
+%% Get static data
+floatReferenceTime = pFloatDetails.phasereftime;
+soundVelocity = soundsourceCombinations(:, 5);
+%startPosition = [pSatData.lat_sat(1) pSatData.lon_sat(1)];
+trackingMethod = lower(pTrackingParameter.trackingMethodString);
+
+trajectory = [];
+
+switch (trackingMethod)
+    case 'least square'
+        trajectory = artoa.float.calculateTrajectoryLeastSquares( ...
+            preparedData, ...
+            soundsourceCombinations, ...startPosition, ...
+            floatReferenceTime, ...
+            soundVelocity ...
+        );
+end
+
+return
+
+%% Calculate positions for every combination segment
+trajectory = startPosition;
+
+for oCombination = 1:size(soundsourceCombinations, 1)
+    trajectory = [ ...
+        trajectory; ...
+        artoa.data.calculateCombinationSegmentLeastSquares( ...
+            preparedData, ...
+            soundsourceCombinations(oCombination, :), ...
+            trajectory(end, :), ...
+            floatReferenceTime, ...
+            soundVelocity(oCombination) ...
+        ) ...
+    ];
+end
+
+% x_new=ls_converge(refpos,length(sosopos),sosopos,dist,svsoso((jx)));
+% if (sum(~isnan(dist)) ~= length(sosopos)) x_new=nan*size(x_new);end
+% pnt=x_new(1,1:2);
+% %disp([num2str(ix),' ',num2str(mx),num2str(pnt(1)*RAD),' ',num2str(pnt(2)*RAD)]);
+% clkerror = NaN;
+
+
+return
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+%% Get required data
+soundsourceCombinations = pTrackingParameter.soundsourceCombinations;
+trackingMethod = pTrackingParameter.trackingMethodString;
+startPosition = [pSatData.lat_sat(1) pSatData.lon_sat(1)];
+
+
+%% Run through all combinations
+
+% get the current combination
+currentCombination = soundsourceCombinations(1, :);
+% isolate the sound source data of this combination
+currentSources = artoa.data.extractSoundsourcesFromStruct( ...
+    strsplit(currentCombination{3}), pSoundsources ...
+);
+% get the rafos date boundaries for this combination
+beginRafosDate = str2double(currentCombination{1});
+endRafosDate = str2double(currentCombination{2});
+% create a selection based on the information above
+
+% equals to true if the dates are in range
+selectionOfDates = (pToaData.toaDate >= beginRafosDate) & (pToaData.toaDate <= endRafosDate);
+
+soundsourceData = struct();
+soundsourcePositions = [];
+for oSoundsource = 1:length(currentSources)
+    % collect the toa points that have been applied to the sound source and
+    % is in date range
+    selectionSoundsource = ...
+        strcmp(pToaData.soundSource, {currentSources{oSoundsource}.sourcename}) ...
+        & selectionOfDates;
+    [soundsourceData(oSoundsource).rafosDate, sortIndices] = sort(pToaData.toaDate(selectionSoundsource));
+    soundsourceData(oSoundsource).toa = pToaData.toa(selectionSoundsource);
+    % sort toa
+    soundsourceData(oSoundsource).toa = soundsourceData(oSoundsource).toa(sortIndices);
+    % get position
+    soundsourcePositions = [soundsourcePositions; currentSources{oSoundsource}.position];
+end
+
+% search all correspondences
+sameDateIndices = [];
+for oDate = 1:length(soundsourceData(1).rafosDate)
+    tmpDate = soundsourceData(1).rafosDate(oDate);
+    tmpDateExists = true;
+    tmpIndices = oDate;
+    for oSoundsource = 2:length(currentSources)
+        tmpDateExists = tmpDateExists & any(soundsourceData(oSoundsource).rafosDate == tmpDate);
+        tmpIndices = [tmpIndices, find(soundsourceData(oSoundsource).rafosDate == tmpDate)];
+    end
+    if tmpDateExists
+        sameDateIndices(end + 1, :) = tmpIndices;
+    end
+end
+
+
+for oSoundsource = 1:length(currentSources)
+    soundsourceData(oSoundsource).distanceToSoundsource = ...
+        (soundsourceData(oSoundsource).toa ...
+        + (pFloatDetails.phasereftime(1)*60 + pFloatDetails.phasereftime(2)) * 60 ...
+        - (currentSources{oSoundsource}.reftime(:, 1)*3600 + currentSources{oSoundsource}.reftime(:, 2) * 60) ... % - leapseconds missing
+        ) * 1464;
+end
+
+calculatedPositions = startPosition;
+
+for oIndex = 1:size(sameDateIndices, 1)
+    sv = 1.464;
+    calculatedPositions = [ ...
+        calculatedPositions; ...
+        artoa.vendor.ls_converge( ...
+            calculatedPositions(end, :), ...
+            length(currentSources), ...
+            soundsourcePositions, ...
+            [soundsourceData(1).distanceToSoundsource(sameDateIndices(oIndex));soundsourceData(2).distanceToSoundsource(sameDateIndices(oIndex))], ...
+            sv ...
+        )* 180 / pi ...
+    ];
+end
+
+% run through all dates
+indicesDates = find(selectionOfDates);
+for oDates = 1:length(indicesDates)
+    % get the toa
+    currentToa = pToaData.toa(indicesDates(oDates));
+    
+end
+
+
+
+
+switch (trackingMethod)
+    case 'Least Square'
+%         gpsRafosDates(o) = artoa.convert.dmy2rd( ...
+%             gpsData(o, pRfb.SAT_FORMAT.day_sat), ...
+%             gpsData(o, pRfb.SAT_FORMAT.month_sat), ...
+%             gpsData(o, pRfb.SAT_FORMAT.year_sat) ...
+%         );
+        %distanceSoundSourceToFloat = artoa.data.calculateGeodist(currentGpsPosition, soundSourcePosition);
+end
+
+
+trajectory = [];
+
+end
\ No newline at end of file

lib/+artoa/+float/calculateTrajectoryLeastSquares.m

+function [trajectory] = calculateTrajectoryLeastSquares(pPreparedData, pSoundsourceCombinations, pFloatReferenceTime, pSoundVelocity)
+%CALCULATETRAJECTORYLEASTSQUARES Summary of this function goes here
+%   Detailed explanation goes here
+
+trajectory = [];
+%trajectory = pStartPosition;
+
+for oCombination = 1:size(pSoundsourceCombinations, 1)
+    segmentPositions = artoa.data.calculateCombinationSegmentLeastSquares( ...
+            pPreparedData, ...
+            pSoundsourceCombinations(oCombination, :), ...%trajectory(end, :), ...
+            pFloatReferenceTime, ...
+            pSoundVelocity{1, oCombination} ...
+    );
+    trajectory = [ ...
+        trajectory; ...
+        segmentPositions ...
+    ];
+end
+
+end
+

lib/+artoa/+float/scatterSatDataPositions.m

+function [scatterHandle, textHandles] = scatterSatDataPositions(pAxesHandle, pSatDataPositions)
+%SCATTERSATDATAPOSITIONS Summary of this function goes here
+%   Detailed explanation goes here
+
+
+%% Remove NaN
+x = pSatDataPositions.lon_sat(~isnan(pSatDataPositions.lon_sat));
+y = pSatDataPositions.lat_sat(~isnan(pSatDataPositions.lat_sat));
+year = pSatDataPositions.year_sat(~isnan(pSatDataPositions.year_sat));
+month = pSatDataPositions.month_sat(~isnan(pSatDataPositions.month_sat));
+day = pSatDataPositions.day_sat(~isnan(pSatDataPositions.day_sat));
+
+%% Calculate step size for text
+steps = 1:3:(floor(length(x)/3)*3);
+
+%% Prepare text cell
+textCell = {};
+for i = steps
+    textCell{end + 1} = {'', ['  ' num2str(year(i)) '-' num2str(month(i)) '-' num2str(day(i))]};
+end
+
+%% Plot
+hold(pAxesHandle, 'on');
+scatterHandle = scatter(pAxesHandle, x, y);
+textHandles = text(pAxesHandle, x(steps), y(steps), textCell, 'FontSize', 7);
+hold(pAxesHandle, 'off');
+
+end
+