%======================================================================
%                    L O A D N A V . M 
%                    doc: Thu Jun 17 18:01:50 2004
%                    dlm: Mon Jan  3 02:01:17 2011
%                    (c) 2004 ladcp@
%                    uE-Info: 252 25 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)

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

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

% =================================================================
% - 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 < 1900 || year > 2015)
	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
