%======================================================================
% L O A D N A V . M
% doc: Thu Jun 17 18:01:50 2004
% dlm: Thu Feb 18 13:38:50 2016
% (c) 2004 ladcp@
% uE-Info: 182 7 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% MODIFICATIONS BY ANT:
% Jun 26, 2004: - totally re-wrote the file-reading, which reduces
% run time by over factor 90(!!!) --- 5.7 instead of
% 527 seconds in the test case --- for deep casts
% Jul 1, 2004: - added support for time bases
% - made input definition flexible
% - removed ipos argument
% Dec 18, 2004: - BUG: header_lines > 0 did not work
% Jun 9, 2005: - improved time-base comments
% Aug 9, 2006: - added support for elapsed-time
% Jan 5, 2007: - removed LADDER-1 specific code for IX_4 distribution
% - added Dan Torres' code to interpolate irregular
% GPS time series
% - added nav file layout into p structure
% Jan 17, 2007: - added support for new [geomag.m]
% Jan 26, 2007: - BUG: file layout default was in p (rather than f) structure
% Jun 30, 2008: - adapted to calculate magdev using external program geomag60
% Jul 2, 2008: - BUG: missing "" in geomag command to allow spcs in
% cof-file path
% Jul 17, 2008: - moved some code to [loadctd.m] to fix bug associated with
% adjusting of start/end positions in case of significant
% ADCP vs GPS/CTD clock offset (unclear whether this happened
% only when elapsed time was used)
% Jul 27, 2008: - nanmean() -> meannan()
% Oct 15, 2008: - replaced mean by median to get lat/lon (bad outliers in L1 data set)
% Dec 1, 2009: - BUG: geomag date check was wrong (Dec 1 2009 resulted in a date >= 2010)
% Jan 22, 2010: - adapted to Eric Firing's much simplified magdec utility
% Jan 3, 2011: - changed IGRF11 validity to end of 2015 (from 2010)
% Feb 18, 2016: - BUG: geomag year range check bombed in 2016
% - added p.interp_missing_GPS using code provided by Jay Hooper
function [d,p]=loadnav(f,d,p)
% function [d,p]=loadnav(f,d,p)
% This routine works for generic ASCII files, containing GPS time series,
% with fields time, lat and lon.
%====================
% TWEAKABLES
%====================
% PRE-AVERAGING OF GPS TIME (IN DAYS)
p = setdefv(p,'navtime_av',2/60/24);
% INTERPOLATE IRREGULAR NAV TIME SERIES (Dan Torres)
p=setdefv(p,'interp_nav_times',0);
% INTERPOLATE MISSING GPS VALUES (Jay Hooper)
p=setdefv(p,'interp_missing_GPS',1);
% FILE LAYOUT
f = setdefv(f,'nav_header_lines',0);
f = setdefv(f,'nav_fields_per_line',3);
f = setdefv(f,'nav_time_field',1);
f = setdefv(f,'nav_lat_field',2);
f = setdefv(f,'nav_lon_field',3);
% TIME BASE
% 0 for elapsed time in seconds
% 1 for year-day (1.0 = Jan 1, 00:00)
% 2 for Visbeck's Gregorian (see gregoria.m)
p = setdefv(p,'nav_time_base',0);
%======================================================================
% MODIFICATIONS BY ANT:
% Jun 26, 2004: - totally re-wrote the file-reading, which reduces
% run time by over factor 90(!!!) --- 5.7 instead of
% 527 seconds in the test case --- for deep casts
% Jul 1, 2004: - added support for time bases
% - made input definition flexible
% - removed ipos argument
% Dec 18, 2004: - BUG: header_lines > 0 did not work
% Jun 9, 2005: - improved time-base comments
% Aug 9, 2006: - added support for elapsed-time
% Jan 5, 2007: - removed LADDER-1 specific code for IX_4 distribution
% - added Dan Torres' code to interpolate irregular
% GPS time series
% - added nav file layout into p structure
% Jan 17, 2007: - added support for new [geomag.m]
% Jan 26, 2007: - BUG: file layout default was in p (rather than f) structure
% Jan 7, 2009: - tightened use of exist()
disp(['LOADNAV: load NAV time series ',f.nav])
if ~exist(f.nav,'file')
disp([' can not find ',f.nav])
return
end
% construct input format
cur_field = 1; input_format = '';
for i=1:f.nav_fields_per_line
switch i,
case f.nav_time_field,
i_time = cur_field; cur_field = cur_field + 1;
input_format = [input_format ' %g'];
case f.nav_lat_field
i_lat = cur_field; cur_field = cur_field + 1;
input_format = [input_format ' %g'];
case f.nav_lon_field
i_lon = cur_field; cur_field = cur_field + 1;
input_format = [input_format ' %g'];
otherwise
input_format = [input_format ' %*g'];
end
end
if cur_field ~= 4
error('File format definition error');
end
% open input & skip header
header_lines = f.nav_header_lines;
fp=fopen(f.nav);
while header_lines > 0
fgets(fp);
header_lines = header_lines - 1;
end
% read time series
[A,nread] = fscanf(fp,input_format,[3,inf]);
% close file
fclose(fp);
% NAV time
d.navtime_jul=A(i_time,:)';
switch f.nav_time_base
case 0 % elapsed time in seconds
d.navtime_jul = d.navtime_jul/24/3600 + julian(p.time_start);
case 1 % year-day
d.navtime_jul = d.navtime_jul + julian([p.time_start(1) 1 0 0 0 0]);
end
disp([' number of NAV scans: ',int2str(length(d.navtime_jul)),...
' delta t : ',num2str(median(diff(d.navtime_jul))*24*3600),' seconds'])
%----------------------------------------
% interpolate to regular time series
% code provided by Dan Torres
%----------------------------------------
if p.interp_nav_times
min_t = min(d.navtime_jul);
max_t = max(d.navtime_jul);
delta_t = median(diff(d.navtime_jul));
data = interp1q(d.navtime_jul,A([i_lat i_lon],:)',[min_t:delta_t:max_t]');
d.navtime_jul = [min_t:delta_t:max_t]';
disp(sprintf(' interpolated to %d NAV scans; delta_t = %.2f seconds',...
length(d.navtime_jul),median(diff(d.navtime_jul))*24*3600));
else
data=A([i_lat i_lon],:)';
end
p.navdata = 1;
d.slat = data(:,1);
d.slon = data(:,2);
%----------------------------------------
% interpolate missign GPS values
% code provided by Jay Hooper
%----------------------------------------
if p.interp_missing_GPS
bad_lon = find(d.slon == -9.990e-29);
if ~isempty(bad_lon),
good_lon = find(d.slon ~= -9.990e-29);
d.slon = interp1(d.navtime_jul(good_lon),d.slon(good_lon),d.navtime_jul);
end
bad_lat = find(d.slat == -9.990e-29);
if ~isempty(bad_lat),
good_lat = find(d.slat ~= -9.990e-29);
d.slat = interp1(d.navtime_jul(good_lat),d.slat(good_lat),d.navtime_jul);
end
end
% =================================================================
% - at this point nav data is in d.navtime_jul, d.slon, d.slat
% and p.navdata is set to 1
% - time shifting & extraction of begin/end position is handled
% in [loadctd.m]
% =================================================================
if ~isfinite(p.drot) % set magdecl
[s,o] = system('magdec');
if s == 1
p.drot = geomag(f,meannan(d.navtime_jul),medianan(d.slat),medianan(d.slon));
else
warn = sprintf('"magdec" not found; using old magdev code with IGRF00',f.IGRF);
disp(['WARNING: ' warn]);
p.warn(size(p.warn,1)+1,1:length(warn))=warn;
p.drot = magdev(medianan(d.slat),medianan(d.slon));
end
[d.ru,d.rv]=uvrot(d.ru,d.rv,p.drot);
[d.bvel(:,1),d.bvel(:,2)]=uvrot(d.bvel(:,1),d.bvel(:,2),p.drot);
disp(sprintf(' corrected for magnetic declination of %.1f deg',p.drot));
end
% ================================================================
function depth=p2z(p,lat)
% !!!!!! USES Z=0 AT P=0 (I.E. NOT 1ATM AT SEA SURFACE)
% pressure to depth conversion using
% saunders&fofonoff's method (deep sea res.,
% 1976,23,109-111)
% formula refitted for alpha(p,t,s) = eos80
% units:
% depth z meter
% pressure p dbars (original in bars, but below
% division by 10 is included)
% latitude lat deg
% checkvalue:
% depth = 9712.654 m
% for
% p = 1000. bars
% lat = 30. deg
% real lat,p
if nargin < 2, lat=54; end
p=p/10.;
x=sin(lat/57.29578);
x=x*x;
gr=9.780318*(1.0+(5.2788e-3+2.36e-5*x)*x)+1.092e-5*p;
depth=(((-1.82e-11*p+2.279e-7).*p-2.2512e-3).*p+97.2659).*p;
depth=depth./gr;
%
function a=datestrj(b)
%
% print julian date string
%
a=datestr(b-julian([0 1 0 0 0 0]));
%
%==============================================================
function [hours]=hms2h(h,m,s);
%HMS2H converts hours, minutes, and seconds to hours
%
% Usage: [hours]=hms2h(h,m,s); or [hours]=hms2h(hhmmss);
%
if nargin== 1,
hms=h;
h=floor(hms/10000);
ms=hms-h*10000;
m=floor(ms/100);
s=ms-m*100;
hours=h+m/60+s/3600;
else
hours=h+(m+s/60)/60;
end
%===============================================================
function dev=geomag(f,date,lat,lon);
% function dev=geomag(f,date,lat,lon);
%
% call SOEST magdec to compute magnetic deviation
% INSTALLATION INSTRUCTIONS:
% - src available as "geomag" at http://currents.soest.hawaii.edu/hg/hgwebdir.cgi
% - on UNIX systems (I tested Linux, MacOSX & FreeBSD),
% 1) compile by typing "make" or "gmake" in source directory
% 2) install by typing "make install" or "gmake install" as root
% 3) test by executing matlab command "system('magdec')"
% if this test produces an error and a return value of 127,
% the path is not set correctly
dstr = gregoria(date); % convert date (approx)
year = dstr(1); month = dstr(2); day = dstr(3);
if (year < 1980 || year > 2031)
error(sprintf('year = %d out of range',year));
end
% execute external program
CMD = sprintf('magdec %g %g %d %d %d',lon,lat,year,month,day);
disp(sprintf('executing %s',CMD));
[status,work] = system(CMD);
if status ~= 0
error(['cannot execute <' CMD '>']);
end
vals = sscanf(work,'%g'); % parse output
if length(vals) ~= 4
error(['unexpected output from <' CMD '>']);
end
dev = vals(1); % return result