%======================================================================

%                    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