New tracking method Singel Source plus Topography tracking.

artoa.ini

 % ARTOA ini file for dimes first three floats
 
 [files]
-%-soundsourcefile    C:\MyData\RAFOS\HAFOS.soso
--soundsourcefile    C:\MyData\RAFOS\TESTwithOFFSETS.soso
--floatfile          ..\HAFOS.flo  % not used in artoa4 yet
--bathyfile          ..\HAFOS.prj  % not used in artoa4 yet
--etopo              C:\MyData\TOPO\etopo1_ice_c.flt
+-soundsourcefile    C:\MyData\RAFOS\HAFOS.soso
+%-soundsourcefile C:\MyData\RAFOS\TESTwithOFFSETS.soso
+-floatfile                  ..\HAFOS.flo  % not used in artoa4 yet
+-bathyfile                ..\HAFOS.prj  % not used in artoa4 yet
+-etopo                    C:\MyData\TOPO\etopo1_ice_c.flt
+-floatProfiles           \\smb.isibhv.dmawi.de\projects\argoAWI\data\apex_analysis\mat\nemo_data.mat
+-trackingContours  C:\MyData\ARTOA\contours_Grid5000_Step50.mat
+
 
 [directory]
--interim        C:\MyData\RAFOS\
--rfb            C:\MyData\RAFOS\
--soundsource    C:\MyData\RAFOS\
+-interim            C:\MyData\RAFOS\
+-rfb                    C:\MyData\RAFOS\
+-soundsource   C:\MyData\RAFOS\
+-ancillary           C:\MyData\ARTOA\
 % the following parameters are not yet used in artoa4
 -argos          ..
 -rft            ..
@@ -23,6 +27,7 @@
 -surface        
 -ric             C:\MyData\RAFOS\
 
+
 [filemask]
 -interim        *.itm
 -rfb            *.rfb

plugins/tracking/SiSoTopo.m

-function [positions, clockError] = SiSoTopo(pData, pSoundsourcePositions, pCombinationDetails, pStartEndPosition)
+function [positions, clockError] = SiSoTopo_v2(pData, pSoundsourcePositions, pCombinationDetails, pStartEndPosition)
 %MYTRACKINGMETHOD A sample function to write your own tracking algorithms
 %
-%   Parameters:
+%   Parameters:plo
 %       pData (struct):
 %           Data structured by soundsource. The toa is already corrected
 %           (e.g. doppler correction etc.)
@@ -21,6 +21,8 @@ function [positions, clockError] = SiSoTopo(pData, pSoundsourcePositions, pCombi
 %       clockError (double):
 %           The estimated clock error.
 
+global artoaConfig;
+global artoaDataInput;
 %% Get reference position and all soundsource names that are involved
 referencePosition = cellfun(@str2double, strsplit(pCombinationDetails.referencePosition{1}));
 soundsourceNames = strsplit(pCombinationDetails.soundsources{1});
@@ -31,40 +33,215 @@ pCombinationDetails
 pSoundsourcePositions
 
 %% Return NaN because this is a sample
-positions = NaN(size(pData.(soundsourceNames{1}).date, 1), 2);
-clockError = [];
+positions  = NaN(size(pData.(soundsourceNames{1}).date, 1), 2);
+clockError = NaN(size(pData.(soundsourceNames{1}).date, 1),1);
+
+%% Additional data needed
+%%%  addpath('C:\Users\oboebel\Documents\MATLAB\FloatTools\lib\utils\seawater_library')  % for pressure to depth conversion
+%%%  wurde nach %%%%%   C:\Users\oboebel\Documents\MATLAB\artoa4\lib\vendor\seawater_library  verschoben 
+
+% the topography file with contour lines
+% % %
+% load('C:\Users\oboebel\Documents\MATLAB\artoa4development\tools\contours_Grid5000_Step50.mat')
+% replaced by next line
+load(artoaConfig.files.trackingContours)
+
+% get the corresponding profile depth.  Profile depths are spaced by 10
+% days.  To match daily rafos fixes, we use linear interpolation in-between
+% to obtain daily pressure readings.  NOte: the is intended for un-grounded
+% floats.  Grounded floats will actually have daily bottom pressure
+% readings. This might be queried by comparing the depth readings in the
+% RAFOS data block with the BOTTOM block pressure readings.
+% global artoaDataInput
+% profileInfo = artoaDataInput.rfb.PROFILE_DEPTH
+% profilePressure = profileInfo(:,8)
+% profileID       = profileInfo(:,1)
+% now the pressure date and the distance date need to be matched.
+%--> interpolate profile pressures based on distance dates
+
+% this needs to be recoded in line with artoa4 structure, rfb structure,
+% and apex to profile data,  currently hard coded to 8882
+%load('C:\Users\oboebel\Documents\MATLAB\artoa4development\floatdata.mat')
+% % % load('\\smb.isibhv.dmawi.de\projects\argoAWI\data\apex_analysis\mat\nemo_data.mat')
+floatProfilesAll = load(artoaConfig.files.floatProfiles);
+% added link from local nemo_data.mat to \\smb.isibhv.dmawi.de\projects\argoAWI\data\apex_analysis\mat\nemo_data.mat'
+
+[floatProfileDataSingle, index2Float] = nemo.float.find(floatProfilesAll.nemo_data,artoaDataInput.rfb.FLOAT.floatname);
+profileNumbers = 1:size(floatProfileDataSingle,2)
+launchposition  = artoaDataInput.rfb.FLOAT.launchpos;
+% % % %%% 8882
+% % % index2Float = 4;
+% % % profileNumbers = 1:77;
+% % % launchposition = [-74.5704, -29.1891];  % plot the launch position  % hardcoded, to be adjusted
+% % % %%% 8885
+% % % index2Float = 6;
+% % % profileNumbers = 1:36;
+% % % launchposition = [-74.5704, -29.1891];   %  dummy, achtung, falsche werte
+% % % %%% 9223
+% % % index2Float = 22;
+% % % profileNumbers = 1:36;
+% % % launchposition = [-70, -8.];  % plot the launch position  % hardcoded, to be adjusted
+cutOff = 1995;   % dbar,  if maxPressure less, float is considered grounded (black points), otherwise plot red points
+%%% Future:  add pop-up dialog to modify parameters if wanted. 
+
+
+% all data is loaded, here we go:
 
 % control plots %%
-figure
-lati = [-78 -60]
-loni = [-70  0]
-h = axesm('MapProjection','stereo','MapLatLimit',lati,'MapLonLimit',loni,'Frame','on','Grid','on') 
-MapHandle = worldmap(lati,loni);
-
-% plot the map
-load('C:\Users\oboebel\Documents\MATLAB\CPTCMAP\WeddellContours.mat')
-helper = fliplr(cArray{3}');
-latMap = helper(:,1);
-lonMap = helper(:,2);
-plotm(latMap,lonMap)   % index 3 = 2000m , then upwards in 50m increments
-hold on
-
-% program uses only first selected sound source
-centerSoSo = pSoundsourcePositions.W18a
-distance = [pData.W18a.distance{:}]';
-distance = distance/60/1.852  %--> distence in arch length
-
-for i = 1:length(distance)
-    [lat,lon] = scircle1(centerSoSo(1),centerSoSo(2),distance(i))
-    %plotm(lat,lon,'r');
-    [xi,yi] = polyxpoly(latMap,lonMap,lat,lon);   % hier muss die zum Druck passende latMap lonMap verwendung finden.
-    plotm(xi,yi,'xk')
-    drawnow
+figure(2)
+lati = [-75 -72];  % Markus Floats   --> move to popup
+loni = [-42 -20];
+
+lati = [-75 -60];   % Neumayer Floats   
+loni = [-35 -5];
+h = axesm('MapProjection','stereo','MapLatLimit',lati,'MapLonLimit',loni,'Frame','on','Grid','on');
+MapHandle = worldmap(lati,loni); hold on;
+
+% % prepare the map
+% helper = fliplr(cArray{3}');
+% latMap = helper(:,1);
+% lonMap = helper(:,2);
+% hold on
+
+sourceName = pCombinationDetails.soundsources;
+sourceName = sourceName{1}
+[token,remain] = strtok(sourceName,' ');
+sourceName = token;
+title(sourceName)
+
+plotm(launchposition(1),launchposition(2),'ko')  % plot the launch position  % hardcoded, to be adjusted
+
+eval(['centerSoSo = pSoundsourcePositions.',sourceName])
+plotm(centerSoSo(1),centerSoSo(2),'ob','MarkerSize',10)
+eval(['distanceVar = [pData.',sourceName,'.distance{:}]']); % these are the distances from the TOA, calculated on the fly by artoa4
+distanceVar = distanceVar';
+distanceVar = distanceVar/60/1.852;  %--> distence in arch length
+eval(['pongDate = [pData.',sourceName,'.date{:}]']);
+pongDate = fliplr(rd2dmy(pongDate'));
+pongDate = datetime((pongDate));  % this is the array of pong dates
+startDate = pongDate(1);
+
+for i = 1:length(profileNumbers)
+% % %     profileData      = nemo_data{index2Float}{i}.PROFILE_DATA;
+    profileData      =  floatProfileDataSingle{i}.PROFILE_DATA;
+    pressureMax(i,1) = profileData(1,7);
+    depthMax(i,1)    = sw_dpth(pressureMax(i,1),centerSoSo(1));
+    %   temperatureMax(i,1) = profileData(1,5);
+    profileTech          = floatProfileDataSingle{i}.PROFILE_TECHNICAL_DATA;
+    pressureMaxTime(i,1) = datetime(profileTech.xmit_ascent_start_time);  % this is when the pofile started
+    profileNumber(i,1)   =  profileTech.xmit_profile_number;
+end
+
+if ~pCombinationDetails.backwardTracking;
+    currentRefPosition = pCombinationDetails.referencePosition;
+elseif pCombinationDetails.backwardTracking
+    currentRefPosition = pCombinationDetails.referencePositionEnd;
+end
+currentRefPosition = cell2mat(currentRefPosition);
+[token,remain] = strtok(currentRefPosition,' ');
+currentRefPosition = [str2num(token), str2num(remain)]
+hp_select = plotm(currentRefPosition(1),currentRefPosition(2),'k*');
+textm(currentRefPosition(1),currentRefPosition(2),'  refPos');
+
+x =  datenum(pressureMaxTime) ;    % the profile timestamps
+v =  depthMax;          % the profile depths at the profile timestamps
+xq = datenum(pongDate);   % the RAFOS timestamps
+vq = interp1(x,v,xq);     % the interpolated profile depths at the RAFOS timestamps
+pongDepth = vq;
+
+%figure(3)   % pongDepth only covers the period of this sound source
+x =  datenum(pressureMaxTime) ;    % the profile timestamps
+v =  pressureMax;          % the profile depths at the profile timestamps
+xq = datenum(pongDate);   % the RAFOS timestamps
+vq = interp1(x,v,xq);     % the interpolated profile depths at the RAFOS timestamps
+% plot(x,v,'o'); hold on
+% plot(xq,vq,'-')
+% ylabel('pressure [dbar]')
+% plot( get(gca,'Xlim')', [cutOff cutOff]')
+
+for i = 1:length(distanceVar)  % loop over all distances / received pongs
+    [lat,lon] = scircle1(centerSoSo(1),centerSoSo(2),distanceVar(i));
+    figure(2)
+%     plotm(lat,lon,'r');
+    % find the contour fitting to this depth - latMap, lon Map
+    currentDepth = pongDepth(i);
+    currentDate  = pongDate(i);
+
+    if ~isnan(currentDepth)
+%       currentDepth = sw_dpth(currentDepth,centerSoSo(1));
+
+        helper = abs(abs(contourLevels)-currentDepth);
+        index2contour = find( min(abs(helper)) == helper) ;  % this is the contour closes to the current position
+        helper = fliplr(cArray{index2contour}');
+        latMap = helper(:,1);
+        lonMap = helper(:,2);
+        
+        index = (latMap > -68) | (latMap < -74.6) | (lonMap < -60);
+        latMap(index) = NaN;
+        lonMap(index) = NaN;
+                
+        figure(2)  
+        plotm(latMap(2:end),lonMap(2:end),'k:','MarkerSize',1)  % plot the contour line
+
+        [xi,yi] = polyxpoly(latMap,lonMap,lat,lon);
+%         P1 = polyshape(latMap,lonMap)
+%         P2 = polyshape(lat,lon)
+%         C = intersect(P1,P2)
+%         xi = C(1)
+%         yi = C(2)
+        
+        if ~isempty(xi) % continue/plot only if intersection exists
+
+            hp_temp = plotm(xi,yi,'.','MarkerSize',10,'Color',[.6 .6 .6]);  % all potential positions in grey
+            [ARCLEN, AZ] = distance(currentRefPosition(1),currentRefPosition(2),xi,yi);  % chose nearest Neighbor
+            indexNearby = find(ARCLEN == min(ARCLEN)) ;
+            xi = xi(indexNearby);
+            yi = yi(indexNearby);
+            %            currentRefPosition = [xi,yi];
+            hp_select = plotm(xi,yi,'.','MarkerSize',10,'Color',[.6 .6 .6]);  % overplot the closest position
+
+            [i length(distanceVar) currentDepth xi yi];
+            if currentDepth > sw_dpth(cutOff,centerSoSo(1)) % float touches bottom; cutoff[dbar]=1995, cutoff[m]=1965 
+                set(hp_select,'color','r'); % float did touch bottom.
+            else
+                if ~pCombinationDetails.backwardTracking;  % ergo forward tracking
+        %           if currentDate < (startDate+360)  % parking depth nominally 2000m
+                    if currentDate < (startDate)      % parking depth nominally 2000m
+                        set(hp_select,'color','k');
+                    else % parking depth nominally 800m
+                        set(hp_select,'color','g');
+                    end
+                elseif pCombinationDetails.backwardTracking;  % ergo backward tracking
+                    if currentDate >= (startDate)  % 
+                        set(hp_select,'color','g');
+                    else % parking depth nominally 800m   % der Teil hier ist noch unlogisch und muss überarbeitet werden
+                        set(hp_select,'color','g');
+                    end  
+                end
+            end
+            positions(i,:) = [xi, yi];   % see line 34
+            clockError(i,:)  = 0;
+            
+            delete(hp_temp) 
+
+            % check if profile label shall be added
+            helper = find(abs(datenum(currentDate)-x) <= 2 ); %0.75  current Pong date is within a day of profile date
+            if ~isempty(helper)
+                textm(xi,yi,num2str(profileNumber(helper)),'HorizontalAlignment','left','VerticalAlignment','bottom')
+            end
+
+        end
+        drawnow
+        pause
+    end
 end
 
 
 
 
+if 1 == 2
+    exportgraphics(gcf,'C:\MyData\APEX\8882.pdf','Append',true)
+end
 
 
 
@@ -72,3 +249,10 @@ end
 end
 
 
+
+
+
+
+
+
+