Extended the data package in preparation for using with trajectory output window.

lib/+artoa/+data/addFloatDrift.m

+function [adjustedToaData] = addFloatDrift(pToaData, pFloatOffsetBegin, pFloatOffsetEnd)
+%ADDFLOATDRIFT Summary of this function goes here
+%   Detailed explanation goes here
+
+%% Prepare return variable
+adjustedToaData = pToaData;
+
+%% Parameter check
+if pFloatOffsetBegin == 0 && pFloatOffsetEnd == 0
+    return;
+end
+if pFloatOffsetBegin == pFloatOffsetEnd
+    adjustedToaData.toa = pToaData.toa + pFloatOffsetEnd;
+end
+
+%% Get required data
+maxToaDate = max(pToaData.toaDate);
+minToaDate = min(pToaData.toaDate);
+
+%% Calculate offset for every toa
+offset = ((pFloatOffsetEnd - pFloatOffsetBegin) / (maxToaDate - minToaDate)) ...
+    * (pToaData.toaDate - minToaDate) + pFloatOffsetBegin;
+
+%% Add calculated offset to toa data
+adjustedToaData.toa = pToaData.toa + offset;
+
+
+end
+

lib/+artoa/+data/calculateCombinationSegmentLeastSquares.m

+function [segmentPositions] = calculateCombinationSegmentLeastSquares(pCorrectedData, pCombinationDetails, pFloatReferenceTime, pSoundVelocity)
+%CALCULATECOMBINATIONSEGMENTLEASTSQUARES Summary of this function goes here
+%   Detailed explanation goes here
+
+%% Get sound sources
+soundsourceNames = strsplit(pCombinationDetails.soundsources{1}, ' ');
+
+%% 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];
+for i = 2:length(soundsourceNames)
+    [intersectedToaDates] = intersect(intersectedToaDates, pCorrectedData.(soundsourceNames{i}).toaDate);
+    soundsourcePositions = [soundsourcePositions; pCorrectedData.(soundsourceNames{i}).position];
+    
+end
+
+%% Remove all dates that are out of bounds
+intersectedToaDates = intersectedToaDates( ...
+    intersectedToaDates >= segmentStart ...
+    & intersectedToaDates <= segmentEnd ...
+);
+
+distanceToSoundsources = {};
+%segmentPositions = [pStartposition];
+segmentPositions = cellfun(@str2double, strsplit(pCombinationDetails.referencePosition{1}));
+for oDates = 1:length(intersectedToaDates)
+    currentDateValue = intersectedToaDates(oDates);
+    distanceToSoundsources{oDates} = [];
+    for i = 1:length(soundsourceNames)
+        currentName = soundsourceNames{i};
+        selection = pCorrectedData.(currentName).toaDate == currentDateValue;
+        filteredValue = pCorrectedData.(currentName).toa(selection);
+        %filteredData.(currentName).toa(oDates) = filteredValue(1);
+        %filteredValue = pCorrectedData.(currentName).toaDate(selection);
+        %filteredData.(currentName).toaDate(oDates) = filteredValue(1);
+        
+        relativeToa = ( ...
+            filteredValue(1) ...
+            + (pFloatReferenceTime(1) * 60 + pFloatReferenceTime(2)) * 60 ...
+            - (pCorrectedData.(currentName).soundsourceReferenceTime(1) * 3600 + pCorrectedData.(currentName).soundsourceReferenceTime(2) * 60 - pCorrectedData.(currentName).leapseconds) ...
+        );
+
+        windowStartTime = ...
+            pFloatReferenceTime(1) * 3600 ...
+            + pFloatReferenceTime(2) * 60 ...
+            - pCorrectedData.(currentName).soundsourceReferenceTime(1) * 3600 ...
+            - pCorrectedData.(currentName).soundsourceReferenceTime(2) * 60;
+        leapseconds = pCorrectedData.(currentName).leapseconds;
+        floatWindowStart = pCorrectedData.(currentName).floatWindowStart * 60;
+    
+        relativeToa = filteredValue(1) ...
+            + windowStartTime ...
+            + leapseconds ...
+            - floatWindowStart;
+        
+        % calculate distance to the source
+        distanceToSoundsources{oDates} = [ ...
+            distanceToSoundsources{oDates}; ...
+            relativeToa * (pSoundVelocity/1000)
+        ];
+%         distanceToSoundsources{oDates} = [ ...
+%             distanceToSoundsources{oDates}; ...
+%             (...
+%                 filteredValue(1) ... %filteredData.(currentName).toa(oDates) ...
+%                 + (pFloatReferenceTime(1) * 60 + pFloatReferenceTime(2)) * 60 ...
+%                 - (pCorrectedData.(currentName).soundsourceReferenceTime(1) * 3600 + pCorrectedData.(currentName).soundsourceReferenceTime(2) * 60 - pCorrectedData.(currentName).leapseconds) ...
+%             ) * (pSoundVelocity/1000)
+%         ];
+    
+    end
+end
+
+for oDates = 1:length(intersectedToaDates)
+    segmentPositions = [segmentPositions; ...
+        rad2deg(artoa.vendor.ls_converge(segmentPositions(end, :), length(soundsourceNames), soundsourcePositions, distanceToSoundsources{oDates}, pSoundVelocity/1000)) ...
+    ];
+end
+
+%% Remove reference position from trajectory
+segmentPositions = segmentPositions(2:end, :);
+
+end

lib/+artoa/+data/calculateOffsetsAndDates.m

+function [offsets, dates] = calculateOffsetsAndDates(pFloatRafosDates, pSoundsourceRafosDates, pSoundsourceOffset, pDrift, pSchedule)
+%CALCULATEOFFSETSANDDATES Summary of this function goes here
+%   Detailed explanation goes here
+
+
+%% Prepare variables
+
+if isnan(pDrift)
+    pDrift = 0;
+end
+
+schedule = pSchedule / 24;
+
+soundsourceBegin = pSoundsourceRafosDates(1);
+soundsourceEnd = pSoundsourceRafosDates(2);
+soundsourceOffsetTotal = [];
+soundsourceDaysTotal = [];
+[m, ~] = size(pSoundsourceRafosDates);
+if m == 1
+    soundsourceEnd = pFloatRafosDates(2);
+end
+
+floatBegin = pFloatRafosDates(1);
+floatEnd = pFloatRafosDates(2);
+
+offsetDate = soundsourceBegin;
+offset = 0;
+if ~isempty(pSoundsourceOffset)
+    offsetDate = artoa.convert.dmy2rd( ...
+        pSoundsourceOffset(:, 3), ...
+        pSoundsourceOffset(:, 2), ...
+        pSoundsourceOffset(:, 1) ...
+    );
+    offset = pSoundsourceOffset(:, 4);
+end
+
+%% Calculations
+
+for i = 1:length(offset) - 1
+    tmpDates = [offsetDate(i):schedule:offsetDate(i + 1)]';
+    tmpOffset = ( ...
+        (offset(i + 1) - offset(i)) ...
+         / (offsetDate(i + 1) - offsetDate(i)) ...
+    ) * [tmpDates - offsetDate(i) - 1] + offset(i);
+    soundsourceOffsetTotal = [soundsourceOffsetTotal; tmpOffset];
+    soundsourceDaysTotal = [soundsourceDaysTotal; tmpDates];
+    clear tmpDates tmpOffset;
+end
+
+floatDays = [offsetDate(end):schedule:soundsourceEnd]';
+tmpOffset = (floatDays - offsetDate(length(offsetDate))) * pDrift + offset(length(offset));
+soundsourceOffsetTotal = [soundsourceOffsetTotal;tmpOffset];
+soundsourceDaysTotal = [soundsourceDaysTotal; floatDays];
+clear tmpOffset floatDays;
+
+%% Determine start and end date
+startDate = soundsourceBegin;
+if floatBegin >= soundsourceBegin
+    startDate = floatBegin;
+end
+endDate = soundsourceEnd;
+if floatEnd <= soundsourceEnd
+    endDate = floatEnd;
+end
+
+totalDays = [startDate:schedule:endDate]';
+index = find((soundsourceDaysTotal <= endDate) & (soundsourceDaysTotal >= startDate));
+
+offsets = soundsourceOffsetTotal(index);
+dates = soundsourceDaysTotal(index);
+
+end
+

lib/+artoa/+data/doppcorr.m

+function [toa_date,toa]=doppcorr(toa_date,toa,sampl_int);
+
+% doppcorr.m
+% 10-NOV-2000, HFurey.  This version based on the centered differencing 
+% method recommended by Phil Richardson and Amy Bower.  
+% Constants in both the old code, and Olaf's code is based
+% on a 250 Hz transmission frequency and 80 second signal length, which should
+% not be assumed constant.  There are currently no input variable that
+% can be used for signal strength and transmission frequency, perhaps
+% add to rfb-file header?
+% Also, if sound velocity is to be used, please make sure chosen sound velocity
+% is used, not the default (usually del grosso).
+
+% ALGORITHM:
+% deltaT = (delta * velocity_towards_soso * signal_start_frequency)/
+%          (signal sweep * sound speed);
+% dTOA = ( diff(TOA) * sv )./(diff(TOA date) * secperday);
+% doppler corrected TOAs = TOAs - (deltaT * dTOAs);
+% Note: sound velocity does not enter into code below because it is cancelled 
+% out in algorithm above. 
+%------------------------------------------------------------------------------
+
+% Doppler Constant:
+% signal_length = 80 seconds
+% signal_start_frequency = 259.375 Hz
+% signal_sweep = 1.523 Hz;
+% sec_per_day = 86400;
+%
+% doppcon = (signal_length*signal_start_frequency)/(signal_sweep*sec_per_day);
+doppcon = 0.151990;
+
+% innocuous.
+[toa_date,ix]=sort(toa_date);
+toa=toa(ix);
+
+% I really need to know sampling interval to correctly apply the time.
+% I am not going to use the diff routine because will not do centered difference.
+% this makes for a long loop - but ebbene.
+dctoa = toa;
+diffdate = diff(toa_date); % n-1 long
+% doppcorrTOA = toa - doppcon*(d(toa)/d(t)), provided d(t) is sampling interval
+% I want to avoid calculating the correction over a period of time > samplerate.
+sampl_int = sampl_int./24; % convert to days
+
+for m = 1:length(toa_date)
+   % if first toa
+	if m==1
+		if diffdate(m)==sampl_int % else a single point, and leave TOA unchanged
+		dctoa(m) = toa(m) - doppcon*((toa(m+1)-toa(m))/sampl_int);
+		end
+	
+	% if last toa
+	elseif m==length(toa_date)
+		if diffdate(m-1)==sampl_int % else a single point, and leave TOA unchanged
+		dctoa(m) = toa(m) - doppcon*((toa(m)-toa(m-1))/sampl_int);
+		end
+
+	% if middle toas
+   else
+		if diffdate(m)==sampl_int & diffdate(m-1)==sampl_int % centered difference
+		dctoa(m) = toa(m) - doppcon*((toa(m+1)-toa(m-1))/(2*sampl_int));
+		elseif diffdate(m)==sampl_int & diffdate(m-1)~=sampl_int % forward difference
+		dctoa(m) = toa(m) - doppcon*((toa(m+1)-toa(m))/sampl_int);
+		elseif diffdate(m)~=sampl_int & diffdate(m-1)==sampl_int % backwards difference
+		dctoa(m) = toa(m) - doppcon*((toa(m)-toa(m-1))/sampl_int);
+		end
+		% if point does not qualify for date restrictions then a stand-alone point,
+		% and should remain unchanged.
+	end
+end
+
+toa=dctoa;
\ No newline at end of file

lib/+artoa/+data/findNearestNeighbor.m

+function [nearestneighbor] = findNearestNeighbor(vector,value)
+
+% small function to find the indices of the neighboring values
+% that would surround 'value'
+% 31july98, HDH.
+% INPUT: vector - the vector of values in which one hopes to find neighbors,
+%        value - the value for which one wishes to find neighbors,
+% OUTPUT: smallneighbor - the neighbor value just less than the 'value',
+%         bigneighbor - the neighbor value just greater than the 'value'.
+% 09 January 2010, HHF. Modified to return nearest neighbor.
+
+
+diffs = vector - value;
+diffs = diffs(:);
+% need to find the two smallest diffs, pos and neg.
+[a,b] = find(diffs<0);
+[c,d] = find(diffs>=0);
+maxneg = max(diffs(a));
+minpos = min(diffs(c));
+if isempty(maxneg) % no value in 'vector' less than 'value'
+    nearestneighbor= vector(1);
+elseif isempty(minpos)  % no value in 'vector' greater than 'value'
+    nearestneighbor= vector(end);
+else % 'value' is in middle of 'vector'
+    smallneighbor = vector(find(diffs==maxneg));
+    bigneighbor = vector(find(diffs==minpos));
+
+    if abs(value-smallneighbor) < abs(value-bigneighbor)
+        nearestneighbor = smallneighbor;
+    else
+        nearestneighbor = bigneighbor;
+    end
+end
\ No newline at end of file

lib/+artoa/+data/interpolateRafosData.m

+function [xint, yint, xover] = interpolateRafosData(x, y, interval, gapsize, method)
+%INTERPOLATERAFOSDATA  interpolate the rafos data over a specified gapsize
+%
+% SYNOPSIS:
+%   function [xint,yint,xover]=interpolateRafosData(x,y,interval,gapsize,method)
+%
+% DESCRIPTION:
+%   INTPRAF4 interpolates data using linear or spline interpolation.
+%   If the gap over which was interpolated is greater than the gapsize
+%   specified the interpolated points over this gap will be set to zero.
+%   Only data points greater than zero will be interpolated.
+%
+% ARGUMENTS:
+%     INPUT:    x,y             Data points to be interpolated over x
+%               gapsize         Maximum gap over which to interpolate. The associated unit is that of the x-vector
+%   I think this comment is incorrect: (Olaf)   gapsize         maximum number of data points over which interpolation shall be executed
+%               interval        
+%               method          String containing the method to use
+%                               for interpolation. If method = 'spline'
+%                               the points will be interpolated using
+%                               spline interpolation. If method = 'linear'
+%                               the points will be interpolated using
+%                               linear interpolation.
+%    OUTPUT:    yint            Interpolated y-data points
+%               xint            Interpolated y-data points.
+%               xover           x-values of interpolated data.
+%   The function return the interpolated y-values in yint. All data points
+%   which have been interpolated are given back by there indicies xint.
+%   The interpolation coverd xover
+%	This is the call
+%   [itoa_date,itoa]=intpraf4(TOA_DATE(dx),toa(dx),interval,gapsize,imethod);
+%  interval = time step of data
+%
+%  17JUL03, HHF:  This is not returning the correct values for 'xover', the flag which
+%  tells whether point has been interpolated or not ... I think that it is returning a
+%  '1' for interped only once every day (24hrs), no matter if all cycles < 24 hrs have been 
+%  interpolated.
+%  17 September 2012 HHF: this program is the problem - with infinite
+%  division, e.g. 8./24, it returns time steps that do not match the
+%  original msgdate.  I have performed '(round(().*100))./100' to clip to teh
+%  hunreths digit all incoming time stamps.  For example, the 8-hourly
+%  shedule with yields infinite divisiton is no clipped to 0.08, 0.42, 0.75
+%  in day fraction, which include the 2-hr init offset fo rthe listening
+%  window.  I will put fix at this function call, but not inside this
+%  function.
+
+% if input arguments are missing....
+if nargin < 5
+  method = 'spline';
+end
+if nargin < 4
+  gapsize = 1;
+end
+
+
+if (strcmpi(method, 'none')) % if the interpolation method was none just return the input
+    yint = y;
+    xint = x;
+    xover = [];
+    return
+end
+
+[x, indx] = sort(x);   % sort the time axis, which seems to consist of at least 2 sebments
+y = y(indx); 
+
+gaps = diff(x);
+% this finds the gaps larger than gapsize, over which no interpolation shall be perfomed
+isgap = find(gaps > gapsize + 1);
+gapstime = NaN(length(isgap), 2);
+for ix = 1:length(isgap)
+    gapstime(ix, :) = [x(isgap(ix)) x(isgap(ix) + 1)];
+end
+
+% this finds the gaps <=  gapsize, over which  interpolation shall be perfomed
+isint = find(gaps <= gapsize + 1 & gaps > interval / 24);
+
+% 17JULY03, HHF: Below finds times at which interp is preformed, BUT, decimates by 1 (i.e., 24 hours).
+% This means that in a 6-hourly float, only one interp per day happens.
+%%$for ix=1:length(isint)
+%%$   int_index_array = [int_index_array, [x(isint(ix))+1:1:x(isint(ix)+1)-1]]; % these are the times at which interp is performed
+%%$end
+
+% 17JUL03, HHF: This is correct, accounts for vacriable cycle length, and will flag interpolated
+% points correctly.
+int_index_array = [];
+for ix = 1:length(isint)
+    int_index_array = [ ...
+        int_index_array, ...
+        [x(isint(ix))+(interval/24):(interval/24):x(isint(ix)+1)-(interval/24)] ...
+        ]; % these are the times at which interp is performed
+end
+
+xint = [min(x):interval / 24:max(x)];           % interpolate time axis
+yint = interp1(x, y, xint, lower(char(method))); % interpolate the data
+
+%generate array with 1 where interpolated, 0 otherwise
+xover = ismember(xint, int_index_array); % HHF, 17JUL03: Ahhhh! This is the problem!
+
+if isempty(isgap)
+    return
+end
+
+% cut out segments which shall not be interpolated:
+for ix = 1:length(isgap)
+  gap = find((xint > gapstime(ix, 1)) & (xint < gapstime(ix, 2)));
+  xint(gap) = ones(size(gap)) * (-1);   %set xint to -1 for gaps
+end
+% and now cut them out:
+jx = find(xint ~= (-1));
+xint = xint(jx);
+yint = yint(jx);
+xover = xover(jx);
+
+end

lib/+artoa/+data/simulshift.m

+function [toa_date,toa]=simulshift(toa_date,toa,sampl_int,windowstart,sosoreftime,FloatReferenceTime);
+
+% 10-NOV-2000, HFurey.  
+% Written for ARTOA2 input, based on code from ARTOA_ARTRK.
+% INPUT:    toa_date    - truncated julian date of each toa (a non-continuous time vector)
+%			toa         - toas of a source that fall on toa_date
+%			sampl_int   - the sample interval: need this so can tell if diff(toa_date) 
+%                         greater than one cycle
+%			windowstart - the start times of every window within a float cycle (e.g., 30 60 90)
+%			sosoreftime - the time the source ponged, in hour minute (e.g., 01 30)
+% OUTPUT:   toa         - toa string adjusted for simultaneity
+%			toa_date    - sorted toa date string (innocuous? relevant? necessary?  keeping cuz
+%				          was in the doppler shift code from Olaf/ARTOA2
+
+
+%------------------------------------------------------------------------------
+
+% this seems OK.  innocuous.
+[toa_date,ix]=sort(toa_date);
+toa=toa(ix);
+
+% I really need to know sampling interval to correctly apply the time.
+sstoa = toa;   % simultaneity shifted toa
+diffdate = diff(toa_date); % n-1 long, in days
+difftoa = diff(toa);
+sampl_int = sampl_int./24;
+
+FloatReferenceTime = FloatReferenceTime(1)*60 + FloatReferenceTime(2);  % OB 2017 03 05 convert [hh mm] to [minutes]
+shifttime = sum(windowstart)./length(windowstart) + FloatReferenceTime; % in minutes  % OB 2017 03 05 added last summand
+sosoreftime = sosoreftime(1)*60 + sosoreftime(2); % in minutes
+timeshift = (shifttime - sosoreftime); % positive if soso ponged before reftime, 
+									   % negative if soso ponged after reftime.
+adjfraction = timeshift./1440; % timeshift/num minutes perday
+
+% simulshiftTOA = toa + adjfraction*(d(toa)/d(t)), provided d(t) is sampling interval
+% I want to avoid calculating the correction over a period of time > samplerate.
+if timeshift ~= 0 % if shift necessary, then proceed ...
+	for m = 1:(length(toa_date)-1) % forward differencing method - no shift to end points.
+        %next line changed by G. Roudaut 04/08/2005 to get ride of roundoff errors
+		%if diffdate(m)==sampl_int % else a single or segment end point; leave TOA unchanged
+		if abs(diffdate(m)-sampl_int) < (sampl_int/10) % else a single or segment end point; leave TOA unchanged
+            sstoa(m) = toa(m) + adjfraction*(difftoa(m)/sampl_int);
+		end
+		% if point does not qualify for date restrictions then a stand-alone point,
+		% and should remain unchanged.
+    
+    end
+end
+
+toa=sstoa;