loadnav.m

%======================================================================
%                    L O A D N A V . M 
%                    doc: Thu Jun 17 18:01:50 2004
%                    dlm: Tue Jan 28 13:24:01 2020
%                    (c) 2004 ladcp@
%                    uE-Info: 43 64 NIL 0 0 72 2 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
%   Sep 14, 2018: - BUG: 2008 code move broke working with single GPS file for entire cruise
%   Sep  4, 2019: - adapted to GK new magdev.m
%		  - added p.magdec_source
%   Jan 28, 2020: - 2008 code move is required because at the loadnav stage p.time_start and
%		    end reflect the LADCP turn-on and -off times

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);

% MAGNETIC DECLINATION
%    There are three different sources for magnetic declination
%	- specify p.drot manually
%	- p.magdec_source = 1	% use external magdec program, if available, otherwise ...
%	- p.magdec_source = 2	% use [magdev.m] from Gerd Krahmann
p=setdefv(p,'magdec_source',2);

% 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

%----------------------------------------------------------------------
% The following code was moved into [loadctd.m] in July 2008 to fix a
% bug that most likely only occurs with elapsed times. In 2018, it was
% discovered that the bug fix introduced another bug that only affects
% data sets with a single GPS time series from which the per-station
% information must be extracted. Therefore, the code was moved back
% here but only called when not working with elapsed time.
% In Jan 2020 I realized that the code cannot be here at all, because 
% p.time_start and end are still the LADCP turn-on and -off times.
% Therefore, I disabled the code here and re-enabled it for all
% time bases in [loadctd.m]. I do not remember what bug this now causes.
%----------------------------------------------------------------------

if 0
  if min(d.navtime_jul)>max(d.time_jul) | max(d.navtime_jul)<min(d.time_jul)
    disp('NAV timeseries does not overlap WRONG STATION????')
    disp(['  NAV-data : ',datestrj(d.navtime_jul(1)),'	to ',datestrj(d.navtime_jul(end))])
    disp([' ADCP-data : ',datestrj(d.time_jul(1)),'  to ',datestrj(d.time_jul(end))])
    disp(' will ignore nav data')
    p.navdata = 0;
    d.navtime_jul = d.time_jul; % make sure same length
    d.slat = d.navtime_jul*NaN;
    d.slon = d.slat;
  else
    if d.navtime_jul(1) > d.time_jul(1)
      disp('NAV timeseries starts after ADCP timeseries: used first NAV value to patch ')
      disp(['  NAV-data : ',datestrj(d.navtime_jul(1)),'  to ',datestrj(d.navtime_jul(end))])
      disp([' ADCP-data : ',datestrj(d.time_jul(1)),'  to ',datestrj(d.time_jul(end))])
      d.navtime_jul(1) = d.time_jul(1);
    end
    
    if d.navtime_jul(end) < d.time_jul(end)
      disp('NAV timeseries ends before ADCP timeseries: used last NAV value to patch ')
      disp(['  NAV-data : ',datestrj(d.navtime_jul(1)),'  to ',datestrj(d.navtime_jul(end))])
      disp([' ADCP-data : ',datestrj(d.time_jul(1)),'  to ',datestrj(d.time_jul(end))])
      d.navtime_jul(end) = d.time_jul(end);
    end
  
    % find valid
    ii=find(diff(d.navtime_jul)>0);
    ii=[ii;length(d.navtime_jul)];
  
    % average over p.navtime_av days
    dt2m=[-p.navtime_av:(1/3600/24):p.navtime_av]';
    slon=medianan(interp1q(d.navtime_jul(ii),d.slon(ii),julian(p.time_start)+dt2m));
    slat=medianan(interp1q(d.navtime_jul(ii),d.slat(ii),julian(p.time_start)+dt2m));
keyboard
    p.nav_start=[fix(slat), (slat-fix(slat))*60, fix(slon), (slon-fix(slon))*60];
    elon=medianan(interp1q(d.navtime_jul(ii),d.slon(ii),julian(p.time_end)+dt2m));
    elat=medianan(interp1q(d.navtime_jul(ii),d.slat(ii),julian(p.time_end)+dt2m));
    p.nav_end=[fix(elat), (elat-fix(elat))*60, fix(elon), (elon-fix(elon))*60];
  
    % interpolate on RDI data
    % this also shortens vectors to length of d.time_jul, which may be the
    % only thing that is really needed
    d.slon=interp1q(d.navtime_jul(ii),d.slon(ii),d.time_jul')';
    d.slat=interp1q(d.navtime_jul(ii),d.slat(ii),d.time_jul')';

    p.poss = [NaN NaN NaN NaN];
    p.pose = [NaN NaN NaN NaN];
  end

  p=setdefv(p,'poss',[NaN NaN NaN NaN]);
  [slat,slon] = pos2str(p.poss(1)+p.poss(2)/60,p.poss(3)+p.poss(4)/60);
  disp([' update start pos  from:',slat,'  ',slon])
  p.poss=p.nav_start;
  [slat,slon] = pos2str(p.poss(1)+p.poss(2)/60,p.poss(3)+p.poss(4)/60);
  disp(['			to:',slat,'  ',slon])
  
  p=setdefv(p,'pose',[NaN NaN NaN NaN]);
  [slat,slon] = pos2str(p.pose(1)+p.pose(2)/60,p.pose(3)+p.pose(4)/60);
  disp([' update end pos    from:',slat,'  ',slon])
  p.pose=p.nav_end;
  [slat,slon] = pos2str(p.pose(1)+p.pose(2)/60,p.pose(3)+p.pose(4)/60);
  disp(['			to:',slat,'  ',slon])

end % if p.navdata

% =================================================================
% - 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
 if p.magdec_source == 1					% external magdec program
   [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 magdev Matlab code');
     disp(['WARNING: ' warn]);
     p.warn(size(p.warn,1)+1,1:length(warn))=warn;
   end
 end

 if ~isfinite(p.drot)
   p.drot = magdev(medianan(d.slat),medianan(d.slon),0,p.time_start(1));
 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