From eafcf6fe89312e649e0b8cfc72d92f63fe281dd2 Mon Sep 17 00:00:00 2001
From: Lewin Probst <info@emirror.de>
Date: Sat, 2 Nov 2019 11:53:29 +0100
Subject: [PATCH] Improved trajectory output window.
+ Added checkbox to enable mercator projection
+ Added checkbox to show/hide position dates
+ Extended position dates by rafos day
+ Added colors for trajectories
/ Improved soundvelocity caluclation
+ Added clock error calculation
---
.../checkboxEnableMercatorProjection.m | 13 ++
.../checkboxShowPositionDates.m | 23 +++
.../+track/+trajectoryOutput/open.m | 5 +-
.../+track/+trajectoryOutput/plot.m | 45 ++++-
.../+track/+trajectoryOutput/plotTrajectory.m | 2 +-
.../tableGeneratedTracksSelect.m | 9 +-
lib/+artoa/+controller/+track/run.m | 52 ++++-
.../+controller/setPlotHandleVisibility.m | 49 +++++
lib/+artoa/+convert/+mercator/imerc.m | 19 ++
lib/+artoa/+convert/+mercator/merc.m | 15 ++
lib/+artoa/+convert/axes2mercator.m | 63 ++++++
lib/+artoa/+data/calculateEllipk4.m | 32 +++
lib/+artoa/+data/interpolateRafosData.m | 23 ++-
.../+float/calculateCombinationSegment.m | 184 ++++++++++++++++++
lib/+artoa/+float/calculateTrajectory.m | 98 +++++-----
lib/+artoa/+gui/+track/trajectoryOutput.m | 15 +-
lib/+artoa/+vendor/xnavai.m | 180 +++++++++++++++++
17 files changed, 751 insertions(+), 76 deletions(-)
create mode 100644 lib/+artoa/+controller/+track/+trajectoryOutput/checkboxEnableMercatorProjection.m
create mode 100644 lib/+artoa/+controller/+track/+trajectoryOutput/checkboxShowPositionDates.m
create mode 100644 lib/+artoa/+controller/setPlotHandleVisibility.m
create mode 100644 lib/+artoa/+convert/+mercator/imerc.m
create mode 100644 lib/+artoa/+convert/+mercator/merc.m
create mode 100644 lib/+artoa/+convert/axes2mercator.m
create mode 100644 lib/+artoa/+data/calculateEllipk4.m
create mode 100644 lib/+artoa/+float/calculateCombinationSegment.m
create mode 100644 lib/+artoa/+vendor/xnavai.m
diff --git a/lib/+artoa/+controller/+track/+trajectoryOutput/checkboxEnableMercatorProjection.m b/lib/+artoa/+controller/+track/+trajectoryOutput/checkboxEnableMercatorProjection.m
new file mode 100644
index 0000000..e842888
--- /dev/null
+++ b/lib/+artoa/+controller/+track/+trajectoryOutput/checkboxEnableMercatorProjection.m
@@ -0,0 +1,13 @@
+function [] = checkboxEnableMercatorProjection(~, event)
+%CHECKBOXUPDATETRACKPARAMETERWINDOW Summary of this function goes here
+% Detailed explanation goes here
+
+global artoaWorkspace;
+
+artoaWorkspace.trajectoryOutput.enableMercatorProjection = logical(event.Source.Value);
+
+%% Replot
+artoa.controller.track.trajectoryOutput.plot();
+
+end
+
diff --git a/lib/+artoa/+controller/+track/+trajectoryOutput/checkboxShowPositionDates.m b/lib/+artoa/+controller/+track/+trajectoryOutput/checkboxShowPositionDates.m
new file mode 100644
index 0000000..b71e3e8
--- /dev/null
+++ b/lib/+artoa/+controller/+track/+trajectoryOutput/checkboxShowPositionDates.m
@@ -0,0 +1,23 @@
+function [] = checkboxShowPositionDates(~, event)
+%CHECKBOXUPDATETRACKPARAMETERWINDOW Summary of this function goes here
+% Detailed explanation goes here
+
+global artoaWorkspace artoaGui;
+
+artoaWorkspace.trajectoryOutput.showPositionDates = logical(event.Source.Value);
+
+%% Get trajectory handles
+trajectoryHandles = artoaGui.trajectoryOutput.trajectoryHandles;
+trajectoryWorkspace = artoaWorkspace.trajectoryOutput.trajectories;
+
+%% Hide or show text
+for i = 1:length(trajectoryHandles)
+ artoa.controller.setPlotHandleVisibility( ...
+ trajectoryHandles{i}.textPositions, ...
+ (artoaWorkspace.trajectoryOutput.showPositionDates && ~trajectoryWorkspace{i}.hidden) ...
+ );
+end
+
+
+end
+
diff --git a/lib/+artoa/+controller/+track/+trajectoryOutput/open.m b/lib/+artoa/+controller/+track/+trajectoryOutput/open.m
index 53dcc34..ed34452 100644
--- a/lib/+artoa/+controller/+track/+trajectoryOutput/open.m
+++ b/lib/+artoa/+controller/+track/+trajectoryOutput/open.m
@@ -17,12 +17,13 @@ callbacks = struct();
callbacks.CloseRequestFcn = @artoa.controller.track.trajectoryOutput.close;
callbacks.buttonDeleteSingle = @artoa.controller.track.trajectoryOutput.buttonDeleteSelectedTrajectory;
callbacks.buttonDeleteAll = '';
-callbacks.checkboxMercator = '';
+callbacks.checkboxMercator = @artoa.controller.track.trajectoryOutput.checkboxEnableMercatorProjection;
callbacks.checkboxTrackVisible = '';
callbacks.buttonTrackInfos = '';
callbacks.buttonVelocities = '';
callbacks.tableGeneratedTracksSelect = @artoa.controller.track.trajectoryOutput.tableGeneratedTracksSelect;
callbacks.checkboxUpdateTrackParameterWindow = @artoa.controller.track.trajectoryOutput.checkboxUpdateTrackParameterWindow;
+callbacks.checkboxShowPositionDates = @artoa.controller.track.trajectoryOutput.checkboxShowPositionDates;
%% Open the gui
artoa.gui.track.trajectoryOutput( ...
@@ -41,6 +42,8 @@ artoa.controller.track.trajectoryOutput.updateTableGeneratedTracks();
%% Update gui
artoaGui.trajectoryOutput.checkboxUpdateTrackParameterWindow.Value = ...
artoaWorkspace.trajectoryOutput.updateTrackParameterWindow;
+artoaGui.trajectoryOutput.checkboxShowPositionDates.Value = ...
+ artoaWorkspace.trajectoryOutput.showPositionDates;
%% Plot everything available
artoa.controller.track.trajectoryOutput.plot();
diff --git a/lib/+artoa/+controller/+track/+trajectoryOutput/plot.m b/lib/+artoa/+controller/+track/+trajectoryOutput/plot.m
index e29cb8d..ed95101 100644
--- a/lib/+artoa/+controller/+track/+trajectoryOutput/plot.m
+++ b/lib/+artoa/+controller/+track/+trajectoryOutput/plot.m
@@ -37,10 +37,10 @@ end
for i = 1:length(trajectories)
% initialize handle container
artoaGui.trajectoryOutput.trajectoryHandles{i} = struct();
- if trajectories{i}.hidden
- continue;
- end
- trajectoryColor = artoaWorkspace.soundsourceColors(mod(trajectories{i}.id, size(artoaWorkspace.soundsourceColors, 1)), :);
+ trajectoryColor = artoaWorkspace.soundsourceColors( ...
+ mod(trajectories{i}.id, size(artoaWorkspace.soundsourceColors, 1)) + 1, ...
+ : ...
+ );
% plot the trajectory itself
[ ...
artoaGui.trajectoryOutput.trajectoryHandles{i}.scatterPositions, ...
@@ -61,7 +61,44 @@ for i = 1:length(trajectories)
trajectories{i}, ...
trajectoryColor ...
);
+
+
+ artoa.controller.setPlotHandleVisibility( ...
+ { ...
+ artoaGui.trajectoryOutput.trajectoryHandles{i}.scatterPositions, ...
+ artoaGui.trajectoryOutput.trajectoryHandles{i}.linePositions, ...
+ artoaGui.trajectoryOutput.trajectoryHandles{i}.textPositions, ...
+ artoaGui.trajectoryOutput.trajectoryHandles{i}.scatterReferencePoints, ...
+ artoaGui.trajectoryOutput.trajectoryHandles{i}.textReferencePoints ...
+ }, ...
+ ~trajectories{i}.hidden ...
+ );
+ if ~trajectories{i}.hidden
+ artoa.controller.setPlotHandleVisibility( ...
+ { ...
+ artoaGui.trajectoryOutput.trajectoryHandles{i}.textPositions ...
+ }, ...
+ artoaWorkspace.trajectoryOutput.showPositionDates ...
+ );
+ end
+end
+
+%% Change trajectory plot to mercator if required
+trajectoryAxes = artoaGui.trajectoryOutput.axesTrajectoryOutput;
+
+if ~artoaWorkspace.trajectoryOutput.enableMercatorProjection
+ trajectoryAxes.UserData = [];
+ trajectoryAxes.XLimMode = 'auto';
+ trajectoryAxes.YLimMode = 'auto';
+ return;
end
+trajectoryAxes = artoaGui.trajectoryOutput.axesTrajectoryOutput;
+artoa.convert.axes2mercator( ...
+ trajectoryAxes, ...
+ false ...
+);
+
+
end
diff --git a/lib/+artoa/+controller/+track/+trajectoryOutput/plotTrajectory.m b/lib/+artoa/+controller/+track/+trajectoryOutput/plotTrajectory.m
index dab4895..2fa166c 100644
--- a/lib/+artoa/+controller/+track/+trajectoryOutput/plotTrajectory.m
+++ b/lib/+artoa/+controller/+track/+trajectoryOutput/plotTrajectory.m
@@ -32,7 +32,7 @@ for i = 1:15:length(t)
textHandles{i} = text( ...
pTrajectory.longitude(i), ...
pTrajectory.latitude(i), ...
- {' '; [num2str(t(i, 3)) '-' pad(num2str(t(i, 2)), 2, 'left', '0') '-' pad(num2str(t(i, 1)), 2, 'left', '0')]}, ...
+ {' '; [num2str(t(i, 3)) '-' pad(num2str(t(i, 2)), 2, 'left', '0') '-' pad(num2str(t(i, 1)), 2, 'left', '0') ', ' num2str(pTrajectory.date(i))]}, ...
'FontSize', 6, ...
'Color', pColor ...
);
diff --git a/lib/+artoa/+controller/+track/+trajectoryOutput/tableGeneratedTracksSelect.m b/lib/+artoa/+controller/+track/+trajectoryOutput/tableGeneratedTracksSelect.m
index 959afb2..02ee753 100644
--- a/lib/+artoa/+controller/+track/+trajectoryOutput/tableGeneratedTracksSelect.m
+++ b/lib/+artoa/+controller/+track/+trajectoryOutput/tableGeneratedTracksSelect.m
@@ -2,7 +2,7 @@ function [] = tableGeneratedTracksSelect(~, event)
%TABLEGENERATEDTRACKSSELECT Summary of this function goes here
% Detailed explanation goes here
-global artoaWorkspace;
+global artoaWorkspace artoaGui;
%% Get data
try
@@ -24,8 +24,13 @@ if selectedColumn == 1
~artoaWorkspace.trajectoryOutput.trajectories{selectedRow}.hidden;
% Update table
artoa.controller.track.trajectoryOutput.updateTableGeneratedTracks();
+ % Hide trajectory
+ artoa.controller.setPlotHandleVisibility( ...
+ artoaGui.trajectoryOutput.trajectoryHandles{selectedRow}, ...
+ ~artoaWorkspace.trajectoryOutput.trajectories{selectedRow}.hidden ...
+ );
% Plot again
- artoa.controller.track.trajectoryOutput.plot();
+ %artoa.controller.track.trajectoryOutput.plot();
end
%% Replace track parameter
diff --git a/lib/+artoa/+controller/+track/run.m b/lib/+artoa/+controller/+track/run.m
index 228595e..495ad66 100644
--- a/lib/+artoa/+controller/+track/run.m
+++ b/lib/+artoa/+controller/+track/run.m
@@ -19,13 +19,53 @@ if ~artoa.data.hasMember(artoaWorkspace, {'trajectoryOutput', 'trajectories'})
artoaWorkspace.trajectoryOutput.trajectories = {};
end
+%% Prepare sound velocity
+soundVelocity = NaN( ...
+ size(artoaWorkspace.trackParameter.soundsourceCombinations, 1), ...
+ 3 ...
+);
+
+%% Get sound velocity
+for i = 1:size(soundVelocity, 1)
+ if strcmpi( ...
+ artoaWorkspace.trackParameter.soundspeedMethodString, ...
+ 'manual' ...
+ ) ...
+ && strcmpi( ...
+ artoaWorkspace.trackParameter.trackingMethodString, ...
+ 'hyperbolic' ...
+ )
+ soundVelocity(i, :) = [ ...
+ artoaWorkspace.trackParameter.soundsourceCombinations.soundspeed1(i), ...
+ artoaWorkspace.trackParameter.soundsourceCombinations.soundspeed2(i), ...
+ artoaWorkspace.trackParameter.soundsourceCombinations.soundspeed3(i) ...
+ ];
+ elseif strcmpi( ...
+ artoaWorkspace.trackParameter.soundspeedMethodString, ...
+ 'manual' ...
+ )
+ soundVelocity(i, :) = ...
+ artoaWorkspace.trackParameter.soundsourceCombinations.soundspeed1(i);
+ else
+ fnames = fieldnames(artoaWorkspace.filteredSoundsources);
+ soundVelocity(i, :) = artoa.data.calculateSoundVelocity( ...
+ artoaWorkspace.temperature(artoaWorkspace.statusTemperature == 1), ...
+ artoaWorkspace.pressure(artoaWorkspace.statusPressure == 1), ...
+ artoaWorkspace.trackParameter.soundspeedMethodString, ...
+ artoaWorkspace.filteredSoundsources.(fnames{1}) ... % WHICH SOUNDSOURCE TO USE IF "soundsource" IS SELECTED??
+ );
+ % TODO: IMPLEMENT OTHER SOUNDVELOCITY OPTIONS
+ end
+end
+
%% Calculate and append track
-[trajectory, trajectoryDates] = artoa.float.calculateTrajectory( ...
+[trajectory, trajectoryDates, trajectoryClockError] = artoa.float.calculateTrajectory( ...
artoaWorkspace.float, ...
artoaWorkspace.toaData, ...
artoaWorkspace.filteredSoundsources, ...
artoaWorkspace.trackParameter, ...
+ soundVelocity, ...
artoaDataInput.ini.leapseconds ...
);
@@ -46,7 +86,12 @@ trajectoryObject.trackParameter = artoaWorkspace.trackParameter;
trajectoryObject.latitude = trajectory(:, 1);
trajectoryObject.longitude = trajectory(:, 2);
trajectoryObject.date = trajectoryDates;
-trajectoryObject.id = artoaWorkspace.trajectoryOutput.trajectories{end}.id + 1;
+trajectoryObject.clockError = trajectoryClockError;
+if isempty(artoaWorkspace.trajectoryOutput.trajectories)
+ trajectoryObject.id = 1;
+else
+ trajectoryObject.id = artoaWorkspace.trajectoryOutput.trajectories{end}.id + 1;
+end
trajectoryObject.hidden = false;
% store to workspace
@@ -58,8 +103,5 @@ artoa.controller.track.trajectoryOutput.open();
%% Plot
artoa.controller.track.trajectoryOutput.plot();
-%% Plot trajectory
-%artoa.controller.track.trajectoryOutput.plotTrajectory(trajectory);
-
end
diff --git a/lib/+artoa/+controller/setPlotHandleVisibility.m b/lib/+artoa/+controller/setPlotHandleVisibility.m
new file mode 100644
index 0000000..f20f38f
--- /dev/null
+++ b/lib/+artoa/+controller/setPlotHandleVisibility.m
@@ -0,0 +1,49 @@
+function [] = setPlotHandleVisibility(pHandles, pVisible)
+%HIDEORSHOWHANDLES Summary of this function goes here
+% Detailed explanation goes here
+
+%% Parameter check
+if ~islogical(pVisible)
+ warning('Parameter pVisible needs to be logical. Skipping...');
+ return;
+end
+
+%% Check data structure
+if iscell(pHandles)
+ processCell(pHandles);
+ return;
+end
+
+if isstruct(pHandles)
+ processStruct(pHandles);
+ return;
+end
+
+%% Check if handle is valid
+if isempty(pHandles) || ~isvalid(pHandles)
+ return;
+end
+
+%% Set visibility
+if pVisible
+ pHandles.Visible = 'on';
+else
+ pHandles.Visible = 'off';
+end
+
+ function processCell(pCell)
+ for i = 1:length(pCell)
+ artoa.controller.setPlotHandleVisibility(pCell{i}, pVisible);
+ end
+ end
+
+ function processStruct(pHandles)
+ fnames = fieldnames(pHandles);
+ for oFnames = 1:length(fnames)
+ artoa.controller.setPlotHandleVisibility(pHandles.(fnames{oFnames}), pVisible);
+ end
+ end
+
+
+end
+
diff --git a/lib/+artoa/+convert/+mercator/imerc.m b/lib/+artoa/+convert/+mercator/imerc.m
new file mode 100644
index 0000000..44a5547
--- /dev/null
+++ b/lib/+artoa/+convert/+mercator/imerc.m
@@ -0,0 +1,19 @@
+function yo = imerc(y);
+
+%IMERC inverse mercator projection
+
+%-----------
+%yo = 2 * atan(exp(y)) - pi/2;
+%yo = yo/pi*180;
+%-----------
+
+c1 = 0.19231594611859;
+c2 = 0.3773144885e-3;
+
+yo = 2 * atan(exp(y)) - pi/2;
+yo = yo/pi*180;
+
+yo = yo + sin(yo/90*pi)*c1 + sin(yo/45*pi)*c2;
+
+
+
diff --git a/lib/+artoa/+convert/+mercator/merc.m b/lib/+artoa/+convert/+mercator/merc.m
new file mode 100644
index 0000000..bbb6848
--- /dev/null
+++ b/lib/+artoa/+convert/+mercator/merc.m
@@ -0,0 +1,15 @@
+function [yo] = merc(y)
+% MERC transfer coordinates according Mercator projection.
+% Yo = MERC(Y) returns the mercator-transformed of Y.
+
+a = 6378137.000;
+b = 6356752.314;
+e = sqrt(a^2-b^2)/a;
+
+rad = pi/180;
+
+s = e * sin(y*rad);
+
+yo = log(tan(pi./4 + (y*rad)/2) .* (((1-s)./(1+s)) .^ (e/2) ));
+
+yo = real(yo);
diff --git a/lib/+artoa/+convert/axes2mercator.m b/lib/+artoa/+convert/axes2mercator.m
new file mode 100644
index 0000000..27bcebc
--- /dev/null
+++ b/lib/+artoa/+convert/axes2mercator.m
@@ -0,0 +1,63 @@
+function [] = axes2mercator(pAxesHandle, pInverse)
+%AXES2MERCATOR Converts the given axes and all its children to mercator projection.
+% Detailed explanation goes here
+
+%% Get required variables
+children = pAxesHandle.Children;
+
+if nargin == 1
+ pInverse = false;
+end
+
+%% Check if already converted
+if strcmp(pAxesHandle.UserData, 'mercator')
+ warning([mfilename ': Given axes has already been converted to mercator! Skipping...']);
+ return;
+end
+
+%% Set units
+pAxesHandle.Units = 'norm';
+
+%% Convert every children
+for i = 1:length(children)
+ currentChild = children(i);
+ % get position, depending on type
+ switch (currentChild.Type)
+ case 'image'
+ disp('Image discovered, skipping...');
+ continue;
+ case 'text'
+ if pInverse
+ currentChild.Position(2) = artoa.convert.mercator.imerc( ...
+ currentChild.Position(2) ...
+ );
+ else
+ currentChild.Position(2) = artoa.convert.mercator.merc( ...
+ currentChild.Position(2) ...
+ );
+ end
+ otherwise
+ if pInverse
+ currentChild.YData = artoa.convert.mercator.imerc( ...
+ currentChild.YData ...
+ );
+ else
+ currentChild.YData = artoa.convert.mercator.merc( ...
+ currentChild.YData ...
+ );
+ end
+ end
+end
+
+%% Set corrected limits
+if pInverse
+ pAxesHandle.YLim = artoa.convert.mercator.imerc(pAxesHandle.YLim);
+else
+ pAxesHandle.YLim = artoa.convert.mercator.merc([pAxesHandle.YLim(1) - 1, pAxesHandle.YLim(2) + 1]);
+end
+
+%% Save state to axes
+pAxesHandle.UserData = 'mercator';
+
+end
+
diff --git a/lib/+artoa/+data/calculateEllipk4.m b/lib/+artoa/+data/calculateEllipk4.m
new file mode 100644
index 0000000..3be6e72
--- /dev/null
+++ b/lib/+artoa/+data/calculateEllipk4.m
@@ -0,0 +1,32 @@
+function z = calculateEllipk4(st,sc)
+%
+% NO ANGLE!
+%
+%function [abstand,ang] = ellipk4(st,sc)
+%
+% routine computes great arc distance z in kilometers on a spheroid
+% with flattening =1/295=.00339 using the andoyer-lambert approximation.
+% z is the distance from sc to st and ang is the bearing from sc to
+% st in true radians rel. to north.
+% the coordinates for sc and st(lat,long) must be in radians. this and
+% other details make this routine different from an earlier version
+% ellips used at yale.
+% 6377. is the equatorial radius of the earth minus the depth of the
+% sound channel(=1200meters).
+% 6378 is the equatorial radius of the earth, approximately.
+%
+% changed to matlab by c.schmid march 94
+%
+
+snst=sin(st(:,1));
+snsc=sin(sc(:,1));
+cfz = snst * snsc + cos(st(:,1)) * cos(sc(:,1)) * cos(st(:,2)-sc(:,2));
+z = 1 - cfz * cfz;
+sfz=sqrt(z);
+z=atan2(sfz,cfz);
+a=snst-snsc;
+b=snst+snsc;
+delz = (z+3*sfz) / (1.-cfz) *a *a;
+delz = .84752e-3 * (delz + (z-3.*sfz) / (1.+cfz) *b *b);
+z=6378.0*(z-delz);
+%the end
diff --git a/lib/+artoa/+data/interpolateRafosData.m b/lib/+artoa/+data/interpolateRafosData.m
index 04e4dd0..d95eee8 100644
--- a/lib/+artoa/+data/interpolateRafosData.m
+++ b/lib/+artoa/+data/interpolateRafosData.m
@@ -60,6 +60,10 @@ if (strcmpi(method, 'none')) % if the interpolation method was none just return
return
end
+if (strcmpi(method, 'cubic'))
+ method = 'pchip';
+end
+
[x, indx] = sort(x); % sort the time axis, which seems to consist of at least 2 sebments
y = y(indx);
@@ -97,6 +101,9 @@ yint = interp1(x, y, xint, lower(char(method))); % interpolate the data
xover = ismember(xint, int_index_array); % HHF, 17JUL03: Ahhhh! This is the problem!
if isempty(isgap)
+ xint = toColumn(xint);
+ yint = toColumn(yint);
+ xover = toColumn(xover);
return
end
@@ -107,8 +114,18 @@ for ix = 1:length(isgap)
end
% and now cut them out:
jx = find(xint ~= (-1));
-xint = xint(jx);
-yint = yint(jx);
-xover = xover(jx);
+xint = toColumn(xint(jx));
+yint = toColumn(yint(jx));
+xover = toColumn(xover(jx));
+
+
+ function [columnVector] = toColumn(vector)
+ if ~iscolumn(vector)
+ columnVector = vector';
+ else
+ columnVector = vector;
+ end
+ end
+
end
diff --git a/lib/+artoa/+float/calculateCombinationSegment.m b/lib/+artoa/+float/calculateCombinationSegment.m
new file mode 100644
index 0000000..b7206f8
--- /dev/null
+++ b/lib/+artoa/+float/calculateCombinationSegment.m
@@ -0,0 +1,184 @@
+function [segmentPositions, segmentDates, clockError] = calculateCombinationSegment(pCorrectedData, pCombinationDetails, pFloatReferenceTime, pSoundVelocity, pTrackingMethod)
+%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]
+% pTrackingMethod (string/char):
+% Contains the tracking algorithm that should be used.
+% Available methods:
+% least square
+% exclusive least square
+% circular
+% hyperbolic
+
+%% Initialize return variables
+segmentPositions = [];
+segmentDates = [];
+clockError = [];
+
+%% Get sound sources
+soundsourceNames = unique(strsplit(pCombinationDetails.soundsources{1}, ' '));
+
+%% Get tracking method
+trackingMethod = lower(pTrackingMethod);
+
+%% Get start and end point
+segmentStart = pCombinationDetails{1, 1};
+segmentEnd = pCombinationDetails{1, 2};
+
+%% Get reference position
+referencePosition = cellfun(@str2double, strsplit(pCombinationDetails.referencePosition{1}));
+
+%% 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 = referencePosition;
+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(i)/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
+
+if strcmp(trackingMethod, 'hyperbolic') && (length(soundsourceNames) < 3)
+ trackingMethod = 'circular';
+end
+
+switch trackingMethod
+ case {'least square', 'exclusive least square'}
+ for oDates = 1:length(intersectedToaDates)
+ if strcmp(pTrackingMethod, 'exclusive least square') ...
+ && (sum(~isnan(distanceToSoundsources{oDates})) ~= length(soundsourcePositions))
+ segmentPositions = [segmentPositions; [NaN, NaN]];
+ clockError = [clockError; NaN(1, length(soundsourceNames))];
+ continue;
+ end
+ segmentPositions = [segmentPositions; ...
+ rad2deg(artoa.vendor.ls_converge(segmentPositions(end, :), length(soundsourceNames), soundsourcePositions, distanceToSoundsources{oDates}, pSoundVelocity(1)/1000)) ...
+ ];
+ clockError = [clockError; NaN(1, length(soundsourceNames))];
+ end
+ case 'circular'
+ for oDates = 1:length(intersectedToaDates)
+ distances = distanceToSoundsources{oDates};
+ segmentPositions = [segmentPositions; ...
+ rad2deg(artoa.vendor.xnavai(deg2rad(soundsourcePositions), distances, deg2rad(referencePosition), 2)) ...
+ ];
+ clockError = [clockError; NaN(1, length(soundsourceNames))];
+ end
+ case 'hyperbolic'
+ for oDates = 1:length(intersectedToaDates)
+ distances = distanceToSoundsources{oDates};
+ if any(isnan(distances))
+ segmentPositions = [segmentPositions; [NaN NaN]];
+ clockError = [clockError; NaN(1, length(soundsourceNames))];
+ continue;
+ end
+ segmentPositions = [segmentPositions; ...
+ rad2deg(artoa.vendor.xnavai(deg2rad(soundsourcePositions), distances, deg2rad(referencePosition), 3)) ...
+ ];
+ currentClockError = NaN(1, 3);
+ for oClockError = 1:3
+ currentClockError(oClockError) = ...
+ ( ...
+ artoa.data.calculateEllipk4(segmentPositions(end, :), ...
+ deg2rad(soundsourcePositions(oClockError,:))) - distances(oClockError) ...
+ ) / (pSoundVelocity(oClockError)/1000);
+ end
+ clockError = [clockError; currentClockError];
+ end
+
+% if length(unique(sosonr)) < 3; % check if 3 sources have been specified in track panel
+% fprintf(1,'Hyperbolic tracking not possible: only 2 sound sources specified, using circular tracking\r');
+% pnt=xnavai(sosopos/RAD,dist,refpos/RAD,2); % use circular by default
+% else
+% if (isnan(dist(1)) | isnan(dist(2)) | isnan(dist(3))); % check if data for all 3 sources is availble
+% tpnt=[tpnt;[NaN,NaN]];
+% else
+% pnt=xnavai(sosopos/RAD,dist,refpos/RAD,3);
+% clkerror(1) = (artoa.data.calculateEllipk4(pnt,sosopos(1,:)/RAD)-dist(1))/svsoso(1);
+% clkerror(2) = (artoa.data.calculateEllipk4(pnt,sosopos(2,:)/RAD)-dist(2))/svsoso(2);
+% clkerror(3) = (artoa.data.calculateEllipk4(pnt,sosopos(3,:)/RAD)-dist(3))/svsoso(3);
+% end
+% end
+
+end
+
+%% Store toa dates
+segmentDates = intersectedToaDates;
+
+%% Remove reference position from trajectory
+segmentPositions = segmentPositions(2:end, :);
+
+end
diff --git a/lib/+artoa/+float/calculateTrajectory.m b/lib/+artoa/+float/calculateTrajectory.m
index c92b94e..ed78297 100644
--- a/lib/+artoa/+float/calculateTrajectory.m
+++ b/lib/+artoa/+float/calculateTrajectory.m
@@ -1,4 +1,4 @@
-function [trajectory, trajectoryDates] = calculateTrajectory(pFloatDetails, pToaData, pSoundsources, pTrackingParameter, pLeapsecondsMatrix)
+function [trajectory, trajectoryDates, trajectoryClockError] = calculateTrajectory(pFloatDetails, pToaData, pSoundsources, pTrackingParameter, pSoundVelocity, pLeapsecondsMatrix)
%CALCULATETRAJECTORY Summary of this function goes here
% Detailed explanation goes here
@@ -87,7 +87,6 @@ for i = 1:length(involvedSoundsources)
end
% subtract offset (corrective value) of sound source
- tmpFnames = fieldnames(pSoundsources);
soundsourceOffsetCorrectiveValue = ...
pTrackingParameter.soundsourceOffsets.(currentSoundsource.sourcename);
currentToa = currentToa - soundsourceOffsetCorrectiveValue(1);
@@ -104,26 +103,34 @@ for i = 1:length(involvedSoundsources)
pTrackingParameter.gapSize, ...
lower(pTrackingParameter.interpolationMethodString) ...
);
- preparedData.(currentSoundsource.sourcename).toaIndex = preparedData.(currentSoundsource.sourcename).toaIndex';
+ %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];
+ 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;
+ preparedData.(currentSoundsource.sourcename).toaDate(oRafosDate) = ...
+ artoa.data.findNearestNeighbor( ...
+ uniqueToaDates, ...
+ preparedData.(currentSoundsource.sourcename).toaDate(oRafosDate) ...
+ );
end
% sort toa by date
- [preparedData.(currentSoundsource.sourcename).toaDate, sortIndices] = sort(preparedData.(currentSoundsource.sourcename).toaDate);
- preparedData.(currentSoundsource.sourcename).toa = preparedData.(currentSoundsource.sourcename).toa(sortIndices);
+ [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 ...
+ ] = artoa.data.simulshift( ...
preparedData.(currentSoundsource.sourcename).toaDate, ...
preparedData.(currentSoundsource.sourcename).toa, ...
pFloatDetails.schedule, ...
@@ -136,12 +143,9 @@ for i = 1:length(involvedSoundsources)
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)), ...
@@ -159,59 +163,45 @@ for i = 1:length(involvedSoundsources)
% 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 ...
+ 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);
+%trackingMethod = lower(pTrackingParameter.trackingMethodString);
trajectory = [];
+trajectoryDates = [];
+trajectoryClockError = [];
-switch (trackingMethod)
- case 'least square'
- [trajectory, trajectoryDates] = artoa.float.calculateTrajectoryLeastSquares( ...
+%% Calculate for every combination
+for oCombination = 1:size(soundsourceCombinations, 1)
+ [segmentPositions, segmentDates, segmentClockError] = artoa.float.calculateCombinationSegment( ...
preparedData, ...
- soundsourceCombinations, ...startPosition, ...
+ soundsourceCombinations(oCombination, :), ...%trajectory(end, :), ...
floatReferenceTime, ...
- soundVelocity ...
- );
- case 'exclusive least square'
- trajectory = artoa.float.calculateTrajectoryExclusiveLeastSquares( ...
- );
-end
-
-return
-
-%% Calculate positions for every combination segment
-trajectory = startPosition;
-
-for oCombination = 1:size(soundsourceCombinations, 1)
+ pSoundVelocity(oCombination, :), ...
+ pTrackingParameter.trackingMethodString ...
+ );
+ if islogical(segmentPositions) && segmentPositions == false
+ trajectory = false;
+ trajectoryDates = false;
+ trajectoryClockError = false;
+ return
+ end
trajectory = [ ...
trajectory; ...
- artoa.data.calculateCombinationSegmentLeastSquares( ...
- preparedData, ...
- soundsourceCombinations(oCombination, :), ...
- trajectory(end, :), ...
- floatReferenceTime, ...
- soundVelocity(oCombination) ...
- ) ...
+ segmentPositions ...
];
+ trajectoryDates = [trajectoryDates; segmentDates];
+ if size(segmentClockError, 2) == 2 && size(trajectoryClockError, 2) == 3
+ segmentClockError(:, 3) = NaN(size(segmentClockError, 1), 1);
+ end
+ trajectoryClockError = [trajectoryClockError; segmentClockError];
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
\ No newline at end of file
diff --git a/lib/+artoa/+gui/+track/trajectoryOutput.m b/lib/+artoa/+gui/+track/trajectoryOutput.m
index b74ac84..2571e52 100644
--- a/lib/+artoa/+gui/+track/trajectoryOutput.m
+++ b/lib/+artoa/+gui/+track/trajectoryOutput.m
@@ -17,7 +17,7 @@ availableCallbacks = { ...
'buttonTrackInfos', ...
'buttonVelocities', ...
'tableGeneratedTracksSelect', ...
- 'buttonGenerateNewTrajectoryId' ...
+ 'checkboxShowPositionDates' ...
};
@@ -47,6 +47,8 @@ set( ...
%% Setup axes
artoaGui.figures.trajectoryOutput.CurrentAxes = axes();
artoaGui.figures.trajectoryOutput.CurrentAxes.Position = [.1 .1 .55 .85];
+artoaGui.figures.trajectoryOutput.CurrentAxes.Color = [.2 .2 .2];
+artoaGui.figures.trajectoryOutput.CurrentAxes.GridColor = [.5 .5 .5];
titleVal = [ ...
'Project: ' artoaWorkspace.float.projectname ...
' - Float ID: ' num2str(artoaWorkspace.float.floatname) ...
@@ -77,7 +79,7 @@ artoaGui.trajectoryOutput.frameTrackList = uipanel( ...
columns = { ...
'', ...
'No.', ...
- 'Soundspeed Method' ...
+ 'Tracking Method' ...
};
artoaGui.trajectoryOutput.tableTrackList = uitable( ...
artoaGui.trajectoryOutput.frameTrackList, ...
@@ -163,14 +165,15 @@ artoaGui.trajectoryOutput.checkboxMercator = uicontrol( ...
'CallBack', pCallbacks.checkboxMercator ...
);
-artoaGui.trajectoryOutput.buttonGenerateNewTrajectoryId = uicontrol( ...
+
+artoaGui.trajectoryOutput.checkboxShowPositionDates = uicontrol( ...
'Parent', artoaGui.trajectoryOutput.frameDisplayOptions, ...
- 'String', 'Generate new trajectory id', ...
- 'Style', 'PushButton', ...
+ 'String', 'Show date positions', ...
+ 'Style', 'CheckBox', ...
'FontSize', 8, ...
'Units', 'normalized', ...
'Position', [.1 .6 .8 .175], ...
- 'CallBack', pCallbacks.buttonGenerateNewTrajectoryId ...
+ 'CallBack', pCallbacks.checkboxShowPositionDates ...
);
artoaGui.trajectoryOutput.buttonTrackInfos = uicontrol( ...
diff --git a/lib/+artoa/+vendor/xnavai.m b/lib/+artoa/+vendor/xnavai.m
new file mode 100644
index 0000000..efdd062
--- /dev/null
+++ b/lib/+artoa/+vendor/xnavai.m
@@ -0,0 +1,180 @@
+% XNAVAI locate a point on a spheroid using circular or hyperbolic navigation
+%
+% SYNOPSIS:
+% function [pnt]=xnavai(sosopos,dist,refpoint,mindex)
+%
+% DESCRIPTION:
+% XNAVAI locates a point from a spheroid using circular or hyperbolic
+% navigation.
+%
+% ARGUMETNS:
+% INPUT: sosopos matrix of the size 2 by 2 which contains the
+% positions (latitude, longitude) of the sound
+% sources used for this point in radiants.
+% The matrix should look like:
+% [soso_lat1, soso_long1;
+% soso_lat2, soso_long2]
+% dist vector which contains the distance of each
+% sound source in km used for this point.
+% refpoint vector with length 2, which contains the
+% position of the refpoint in radiants.
+% The vector should look like:
+% [ref_lat, ref_lon]
+% mindex =2 --> circular navigation
+% =3 --> hyperbolic navigation
+%
+% OUTPUT: pnt vector with length 2, which contains the
+% calculated floatposition in radiants.#
+% The vector will look like:
+% [pnt_lat, pnt_lon]
+%
+% Maximum number of iteration steps set to 20
+
+function [pnt]=xnavai(sosopos,dist,refpoint,mindex)
+
+a(1:3)=999;
+dt1=999;
+dt2=999;
+dt3=999;
+dt4=999;
+dt5=999;
+dt6=999;
+
+TOL=0.1;
+INCR=0.000156788;
+ERAD=6377;
+
+imax=20;
+icnt=0;
+stop=0;
+pnt=[NaN, NaN];
+
+if mindex==2 % the circular case
+ if ~(isnan(dist(1)) || isnan(dist(2))) % if neither dist(1) nor dist(2) = NaN
+ d1=dist(1);
+ d2=dist(2);
+ sosoposition(1,:) = sosopos(1,:);
+ sosoposition(2,:) = sosopos(2,:);
+ elseif ~(isnan(dist(1)) || isnan(dist(3))) % if neither dist(1) nor dist(3) = NaN
+ d1=dist(1);
+ d2=dist(3);
+ sosoposition(1,:) = sosopos(1,:);
+ sosoposition(2,:) = sosopos(3,:);
+ elseif ~(isnan(dist(2)) || isnan(dist(3))) % if neither dist(1) nor dist(3) = NaN
+ d1=dist(2);
+ d2=dist(3);
+ sosoposition(1,:) = sosopos(2,:);
+ sosoposition(2,:) = sosopos(3,:);
+ end
+elseif mindex==3 % the hyperbolic case
+ d1=dist(2)-dist(1);
+ d2=dist(3)-dist(1);
+ sosoposition(1,:) = sosopos(1,:);
+ sosoposition(2,:) = sosopos(2,:);
+ sosoposition(3,:) = sosopos(3,:);
+end
+
+istop=0;
+
+while ((icnt < imax) && ~stop)
+
+ istop=istop+1;
+
+ if istop > 20
+ pnt=[NaN, NaN];
+ stop=1;
+ end
+
+ ps=refpoint;
+
+ a(1)=artoa.data.calculateEllipk4(ps,sosoposition(1,:));
+ a(2)=artoa.data.calculateEllipk4(ps,sosoposition(2,:));
+ if mindex==3
+ a(3)=artoa.data.calculateEllipk4(ps,sosoposition(3,:));
+ end
+
+
+ if ((a(1)<0) || (a(2)<0) || (a(3)<0))
+ stop=1;
+ else
+
+ if mindex==2
+ cr1=d1-a(1);
+ cr2=d2-a(2);
+ elseif mindex==3
+ cr1=d1-a(2)+a(1);
+ cr2=d2-a(3)+a(1);
+ end
+ % check if the error is small enough
+ err=abs(cr1)+abs(cr2);
+
+ if (err <= TOL)
+ pnt=refpoint;
+ stop=1;
+ else
+ % if not, go on incrementing position
+ ps(1)=refpoint(1);
+ ps(2)=refpoint(2)+INCR;
+ dt1=artoa.data.calculateEllipk4(ps,sosoposition(1,:));
+ dt2=artoa.data.calculateEllipk4(ps,sosoposition(2,:));
+ if mindex==3
+ dt3=artoa.data.calculateEllipk4(ps,sosoposition(3,:));
+ end
+
+
+ if ((dt1<0) || (dt2<0) || (dt3<0))
+ stop=1;
+ else
+
+ ps(1)=refpoint(1)+INCR;
+ ps(2)=refpoint(2);
+ dt4=artoa.data.calculateEllipk4(ps,sosoposition(1,:));
+ dt5=artoa.data.calculateEllipk4(ps,sosoposition(2,:));
+ if mindex==3
+ dt6=artoa.data.calculateEllipk4(ps,sosoposition(3,:));
+ end
+
+ if ((dt4<0) || (dt5<0) || (dt6<0))
+ stop=1;
+ else
+ if mindex==2
+ b1=dt1-a(1);
+ b2=dt2-a(2);
+ b3=dt4-a(1);
+ b4=dt5-a(2);
+ elseif mindex==3
+ b1=(dt2-dt1)-(a(2)-a(1));
+ b2=(dt3-dt1)-(a(3)-a(1));
+ b3=(dt5-dt4)-(a(2)-a(1));
+ b4=(dt6-dt4)-(a(3)-a(1));
+ end
+
+ dnom=(b1*b4-b2*b3)*ERAD;
+
+ c(2)=(cr1*b4-b3*cr2)/dnom;
+ c(1)=(cr2*b1-b2*cr1)/dnom;
+
+ im=max(abs(c))==abs(c);
+ cc=c(im);
+
+ if mindex==2
+ im=min(abs(a(1:2)))==abs(a(1:2));
+ aa=a(im);
+ elseif mindex==3
+ im=min(abs(a(1:3)))==abs(a(1:3));
+ aa=a(im);
+ end
+
+ cp=2*cc*ERAD/aa;
+
+ if cp < 1
+ cp=1;
+ end
+
+ refpoint=refpoint+c/cp;
+
+ end
+ end
+ end
+ end
+end
--
GitLab