--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/loadnav.m Tue Oct 20 16:23:49 2009 -0400
@@ -0,0 +1,267 @@
+%======================================================================
+% L O A D N A V . M
+% doc: Thu Jun 17 18:01:50 2004
+% dlm: Fri Mar 5 15:49:00 2010
+% (c) 2004 ladcp@
+% uE-Info: 167 13 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
+
+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 >= 2011)
+ 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