Skip to content
Snippets Groups Projects
Commit 6c63433f authored by leprob001's avatar leprob001
Browse files

Extended data package in preparation for GPS fixes TOA plot.

parent 3c7d9352
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
lib/+artoa/+data/calculateEllipk2.m 0 → 100644
+ 59
− 0
View file @ 6c63433f
function [abstand,ang] = calculateEllipk2(st,sc)
%function [abstand,ang] = ellipk2(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
% Version 2.0 14 Nov 1996 R.Goldsmith
% 2.0 fixed angle computation as is was not done correctly
% if the angle approached 270 degrees. 270.1 flipped
% over to 90.1. Might still be wrong for other cases
% but works for my case.
%
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;
if z > 0
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);
% compute bearing...
ang1 = cos(sc(:,1)) .* tan(st(:,1)) - sin(sc(:,1)) .* cos(sc(:,2)-st(:,2));
if ang1 ~= 0,
ang = atan(sin(sc(:,2)-st(:,2))./ang1);
for ii=1:length(st(:,1)),
if ang1 < 0,
if sc(ii,1) > st(ii,1),
ang = ang + pi;
end
if sc(ii,1) < st(ii,1),
ang = ang - pi;
end
end
ang = 2 * pi - ang;
ang = ang / pi * 180;
end
end
else
z=-1.0;
% fprintf('error because distance <= 0')
end
abstand=z;
%the end
lib/+artoa/+data/calculateGeodist.m 0 → 100644
+ 49
− 0
View file @ 6c63433f
% GEODIST calculate the distance and bearing between two positions
%
% SYNOPSIS:
% function [dist,ang] = geodist(pos1,pos2,unit)
%
% DESCRIPTION:
% Calculates the distance and bearing between pos1 and pos2. The optional
% parameter 'unit' specifies if the positions are in degrees ('degree') or
% in radiant ('rad'). If unit is not given, 'degree' would be the default.
%
% Goedist returns distance from pos1 to pos2 (dist) and the bearing from
% pos1 to pos2 (ang) in degrees rel. to north.
%
%
% rlat=[pos1(1),pos2(1)]
% rlon=[pos1(2),pos2(2)]
%
% calls function ellipk2 to calculate distance and bearing
function [dist,ang] = calculateGeodist(pos1,pos2,unit)
if nargin < 3
unit='degree';
end
rlat=[pos1(1),pos2(1)];
rlon=[pos1(2),pos2(2)];
if strcmp(unit,'degree')
% rad = degree/180*pi !
rlat = rlat / 180 * pi;
rlon = rlon / 180 * pi;
end
s=size(rlat);
if s(1) > 1
% make sure that rlat and rlon have more rows than columns
if s(1) < s(2)
rlat=rlat';
rlon=rlon';
end
end
st=[rlat(:,1) rlon(:,1)];
sc=[rlat(:,2) rlon(:,2)];
[dist,ang] = artoa.data.calculateEllipk2(st,sc);
% dist is the distance from sc to st and ang is the bearing from sc to
% st in degrees rel. to north.
lib/+artoa/+data/calculateLeapseconds.m 0 → 100644
+ 39
− 0
View file @ 6c63433f
function lsec=calculateLeapseconds(date1,date2)
% LEAPSEC return the leapseconds between two dates
% (dates must be in rafos day (use dmy2rd)
% lsec=lsec(date1)-lsec(date2)!
lsec=[];
dmys=[ 1 1 1972 10;
1 7 1972 11;
1 1 1973 12;
1 1 1974 13;
1 1 1975 14;
1 1 1976 15;
1 1 1977 16;
1 1 1978 17;
1 1 1979 18;
1 1 1980 19;
1 7 1981 20;
1 7 1982 21;
1 7 1983 22;
1 7 1985 23;
1 1 1988 24;
1 1 1990 25;
1 1 1991 26;
1 7 1992 27;
1 7 1993 28;
1 7 1994 29;
1 1 1996 30;
1 7 1997 31;
1 1 1999 32;
1 1 2006 33];
rd=dmy2rd(dmys(:,1:3));
leap=dmys(:,4);
ax=find(rd <= date1);
bx=find(rd <= date2);
lsec=leap(max(ax))-leap(max(bx));
lib/+artoa/+data/calculateSoundVelocity.m 0 → 100644
+ 35
− 0
View file @ 6c63433f
function [calculatedSoundVelocity] = calculateSoundVelocity(pTemperature, pPressure, pMethod, pSoundSource)
%CALCULATESOUNDVELOCITY Calculates the sound velocity based on the given values.
% Parameters:
% pTemperature The temperature array.
% pPressure The pressure array.
% pMethod The chosen method. Available methods are:
% del grosso
% linear
% soundsource
% pSoundSource The sound source data. Only required if pMethod is
% set to soundsource.
%
% Returns:
% The calculated sound velocity in m/s.
switch lower(pMethod) % del grosso is default
case 'linear'
indices = (pTemperature > 0 & pTemperature < 20);
meanTemperature = mean(pTemperature(indices));
if meanTemperature > 0
calculatedSoundVelocity = (meanTemperature - 7) * 0.0011 + 1490;
else
calculatedSoundVelocity = 1490;
end
clear indices;
case 'soundsource'
calculatedSoundVelocity = pSoundSource.sound_speed;
otherwise
calculatedSoundVelocity = artoa.vendor.oceanstoolbox.sndspd( ...
35, mean(pTemperature), mean(pPressure), 'del grosso' ...
);
end
end
lib/+artoa/+data/predictToaFromGps.m 0 → 100644
+ 98
− 0
View file @ 6c63433f
function [gpsRafosDates, predictedToa] = predictToaFromGps(pRfb, pSoundSource, pSoundSpeedParameter)
%ESTIMATETOAFROMGPS Summary of this function goes here
% Detailed explanation goes here
%% Collect Float data
gpsData = pRfb.SAT_DATA;
floatSchedule = pRfb.FLOAT.schedule;
floatPhasereftime = pRfb.FLOAT.phasereftime(1) + pRfb.FLOAT.phasereftime(2)/60;
floatWindowStart = pRfb.FLOAT.windowstart;
floatWindowLength = pRfb.FLOAT.windowlength;
floatLaunchTime = artoa.convert.dmy2rd( ...
pRfb.FLOAT.launchtime(3), ...
pRfb.FLOAT.launchtime(2), ...
pRfb.FLOAT.launchtime(1) ...
) + pRfb.FLOAT.launchtime(4)/24 + pRfb.FLOAT.launchtime(5)/24/60;
% floatSurfaceEnd = artoa.convert.dmy2rd( ...
% pRfb.FLOAT.cycle(13), pRfb.FLOAT.cycle(12), pRfb.FLOAT.cycle(11) ...
% ) + pRfb.FLOAT.cycle(14)/24 + pRfb.FLOAT.cycle(15)/24/60;
%% Collect sound source data
soundSourceTime = pSoundSource.reftime(1) + pSoundSource.reftime(2)/60;
%soundSourceSchedule = pSoundSource.schedule;
%soundSourceLaunchDate = [pSoundSource.launchtime(3) pSoundSource.launchtime(2) pSoundSource.launchtime(1)];
soundSourceEmissionBegin = artoa.convert.dmy2rd( ...
[pSoundSource.begemis(3) pSoundSource.begemis(2) pSoundSource.begemis(1)] ...
) + pSoundSource.begemis(4)/24 + pSoundSource.begemis(5)/24/60;
%soundSourceOffset = 0;
% if (~isempty(pSoundSource.offset))
% soundSourceOffset = pSoundSource.offset(4);
% end
soundSourceDrift = 0;
if ~isempty(pSoundSource.drift) && ~isnan(pSoundSource.drift)
soundSourceDrift = pSoundSource.drift;
end
soundSourcePosition = pSoundSource.position;
%% Collect generic variables
leapSeconds = artoa.data.calculateLeapseconds( ...
floatLaunchTime, ...
soundSourceEmissionBegin ...
);
%% Calculations
loopLength = floor(24/floatSchedule);
floatTime = NaN(loopLength, 1);
for o = 1:loopLength
floatTime(o) = mod(floatPhasereftime + (o-1) * floatSchedule, 24);
end
floatTime = sort(floatTime);
travelCorrection = soundSourceTime - floatTime;
maxTime = (floatWindowStart + floatWindowLength) / 60;
minTime = floatWindowStart / 60;
indicesInRange = (travelCorrection < maxTime) & (travelCorrection >= minTime);
windowStartTime = min(travelCorrection(indicesInRange));
clear maxTime minTime indicesInRange;
% sound xmit time (seconds) = (hours*60 plus minutes)*60 -leapseconds - sosooffset
%soundXmitTime = windowStartTime * 3600 - leapSeconds + soundSourceOffset + pAdditionalOffset;
soundSourceDriftSinceLaunch = soundSourceDrift ...
* (floatLaunchTime - soundSourceEmissionBegin);
%soundSourceDriftSinceSurfaceEnd = soundSourceDrift ...
% * (floatSurfaceEnd - artoa.convert.dmy2rd(soundSourceEmissionBegin));
%predictedToaAtLaunch = soundXmitTime + soundSourceDriftSinceLaunch;
%predictedToaAtSurface = soundXmitTime + soundSourceDriftSinceSurfaceEnd;
%% Predict toa for every gps position
gpsMeasurementCount = size(gpsData, 1);
predictedToa = zeros(gpsMeasurementCount, 1);
gpsRafosDates = NaN(gpsMeasurementCount, 1);
for o = 1:gpsMeasurementCount
currentGpsPosition = [str2double(gpsData{o, pRfb.SAT_FORMAT.lat_sat}), str2double(gpsData{o, pRfb.SAT_FORMAT.lon_sat})];
gpsRafosDates(o) = artoa.convert.dmy2rd( ...
str2double(gpsData{o, pRfb.SAT_FORMAT.day_sat}), ...
str2double(gpsData{o, pRfb.SAT_FORMAT.month_sat}), ...
str2double(gpsData{o, pRfb.SAT_FORMAT.year_sat}) ...
);
currentSoundsourceDrift = soundSourceDrift * (gpsRafosDates(o) - soundSourceEmissionBegin);
% calculate distance between sound source and gps
distanceSoundSourceToFloat = artoa.data.calculateGeodist(currentGpsPosition, soundSourcePosition);
if (isempty(windowStartTime))
predictedToa(o) = NaN;
else
predictedToa(o) = windowStartTime * 3600 - leapSeconds + currentSoundsourceDrift ...
+ distanceSoundSourceToFloat / (artoa.data.calculateSoundVelocity( ...
pSoundSpeedParameter.temperature, pSoundSpeedParameter.pressure, ...
pSoundSpeedParameter.method, pSoundSpeedParameter.soundSource ...
)/1000); % TODO: CHECK IF /1000 IS CORRECT
end
end
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment