# HG changeset patch
# User A.M. Thurnherr <athurnherr@yahoo.com>
# Date 1586789566 14400
# Node ID 56bdfe65a697932358ca7acfdedc66c48f483488
# Parent  15c603bc4f70920b5f10e9f0fd8da61b9df1f799
.

diff --git a/.interp.ADCP b/.interp.ADCP
new file mode 100644
--- /dev/null
+++ b/.interp.ADCP
@@ -0,0 +1,90 @@
+#======================================================================
+#                    . I N T E R P . A D C P 
+#                    doc: Fri Apr 16 16:07:48 2010
+#                    dlm: Tue Aug  9 23:07:35 2011
+#                    (c) 2010 A.M. Thurnherr
+#                    uE-Info: 81 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# interpolation scheme mimicking RDI ADCP response as described in
+# Broadband primer, p.17
+
+# HISTORY:
+# 	Apr 16, 2010: - created
+#	Aug  9, 2011: - added -u
+
+# NOTES:
+#	- interface is described in [.interp.linear]
+
+$IS_opts = 'b:u';
+$IS_optsUsage = '[pass -u)nfiltered] -b)in/pulse <length>';
+
+sub IS_usage()
+{
+	&antsUsageError()
+		unless defined($opt_b);
+}
+
+sub IS_init($$$$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_interpolate(br,idr,xf,xv,xi,f)	interpolate field f
+#		br							data buffer reference
+#		idr							init-data reference
+#		xf							x field
+#		xv							x value
+#		xi							index of last record with x-value <= x
+#		f							field number to interpolate
+#		<ret val>					interpolated value
+#
+# NB:
+#	- handle f == xf
+#	- return NaN if any of the y values required is NaN
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_interpolate($$$$$$)
+{
+	my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+	
+	return $xv if ($xf == $f);							# return target x
+
+	my($tow) = $xi;										# top of triangular sampling window
+	while ($tow>=0 && (!numberp($bR->[$tow][$xf]) || $bR->[$tow][$xf] > $xv-$opt_b)) {
+		$tow--;
+	}
+	if ($tow < 0) {										# incomplete window
+		return nan unless ($opt_u);
+		$tow = 0;
+	} else {
+		$tow++;
+	}
+
+	my($bow) = $xi+1;									# bottom of triangular sampling window
+	while ($bow<=$#{$bR} && (!numberp($bR->[$bow][$xf]) || $bR->[$bow][$xf] < $xv+$opt_b)) {
+		$bow++;
+	}
+
+	if ($bow > $#{$bR}) {								# incomplete window
+		return nan unless ($opt_u);
+		$bow = $#{$bR};
+	} else {
+		$bow--;
+	}
+
+	my($sweight) = 0;									# calculate weighted average
+	my($sum) = 0;
+	for (my($i)=$tow; $i<=$bow; $i++) {
+		next unless (numberp($bR->[$i][$xf]) && numberp($bR->[$i][$f]));
+		my($weight) = 1 - abs($bR->[$i][$xf]-$xv)/$opt_b;
+		$sum += $weight * $bR->[$i][$f];
+		$sweight += $weight;
+	}
+
+	return ($sweight>0) ? $sum/$sweight : nan;
+}
+
+#----------------------------------------------------------------------
+
+1;
diff --git a/.interp.fillnan b/.interp.fillnan
new file mode 100644
--- /dev/null
+++ b/.interp.fillnan
@@ -0,0 +1,37 @@
+#======================================================================
+#                    . I N T E R P . F I L L N A N 
+#                    doc: Tue Sep 18 16:36:28 2012
+#                    dlm: Tue Sep 18 16:43:36 2012
+#                    (c) 2012 A.M. Thurnherr
+#                    uE-Info: 34 5 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+# fill missing with nans
+
+# HISTORY:
+#	Sep 18, 2012: - adapted from [.interp.nnbr]
+
+# see [.interp.linear] for documentation of interface
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+sub IS_usage() {}
+sub IS_init($$$$) {}
+
+sub IS_interpolate($$$$$$)
+{
+	my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+
+	if ($f == $xf) {
+		return $xv;
+    } elsif ($bR->[$xi][$xf] == $xv) {
+		return $bR->[$xi][$f];
+    } elsif ($bR->[$xi+1][$xf] == $xv) {
+		return $bR->[$xi+1][$f];
+	} else {
+		return nan;
+	}
+}
+
+1;
diff --git a/.interp.linear b/.interp.linear
new file mode 100644
--- /dev/null
+++ b/.interp.linear
@@ -0,0 +1,87 @@
+#======================================================================
+#                    . I N T E R P . L I N E A R 
+#                    doc: Wed Nov 22 21:01:09 2000
+#                    dlm: Tue Aug  5 14:13:14 2008
+#                    (c) 2000 A.M. Thurnherr
+#                    uE-Info: 23 47 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# linear interpolation
+
+# HISTORY:
+# 	Nov 22, 2000: - created
+#	Jul 28, 2006: - Version 3.3 [HISTORY]
+#				  - added xf to ISInit() args
+#	Aug 22, 2006: - modified to allow use with [match]
+#	Aug  5, 2008: - added idr param to IS_init()
+
+# NOTES:
+#	- the [.interp.*] routines assume that x increases strictly monotonically;
+#	  for utilities dealing with non-monotonic data, use [.nminterp.*] instead
+#	- the [.interp.*] routines are written to work on multiple buffers,
+#	  rather than just @ants_; this implies that data created by IS_init()
+#	  must be stored separately for each buffer
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$IS_opts			string of allowed options
+#	$IS_optsUsage		usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_usage() 	mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_usage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_init(br,f,xf)			init interpolation of field f
+#		br					data buffer reference
+#		idr					init-data reference
+#       f					field number
+#		xf					x field number
+#       <ret val>           none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_init($$$$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_interpolate(br,idr,xf,xv,xi,f)	interpolate field f
+#		br							data buffer reference
+#		idr							init-data reference
+#		xf							x field
+#		xv							x value
+#		xi							index of last record with x-value <= x
+#		f							field number to interpolate
+#		<ret val>					interpolated value
+#
+# NB:
+#	- handle f == xf
+#	- return NaN if any of the y values required is NaN
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_interpolate($$$$$$)
+{
+	my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+	return $xv if ($xf == $f);
+	return nan unless (numberp($bR->[$xi][$f]) && numberp($bR->[$xi+1][$f]));
+
+	my($sc) = ($xv - $bR->[$xi][$xf]) / ($bR->[$xi+1][$xf] - $bR->[$xi][$xf]);
+	return $bR->[$xi][$f] + $sc * ($bR->[$xi+1][$f] - $bR->[$xi][$f]);
+}
+
+#----------------------------------------------------------------------
+
+1;
diff --git a/.interp.linear.orig b/.interp.linear.orig
new file mode 100644
--- /dev/null
+++ b/.interp.linear.orig
@@ -0,0 +1,90 @@
+#======================================================================
+#                    . I N T E R P . L I N E A R 
+#                    doc: Wed Nov 22 21:01:09 2000
+#                    dlm: Fri Sep 23 16:32:20 2011
+#                    (c) 2000 A.M. Thurnherr
+#                    uE-Info: 62 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# linear interpolation
+
+# HISTORY:
+# 	Nov 22, 2000: - created
+#	Jul 28, 2006: - Version 3.3 [HISTORY]
+#				  - added xf to ISInit() args
+#	Aug 22, 2006: - modified to allow use with [match]
+#	Aug  5, 2008: - added idr param to IS_init()
+#	Sep 23, 2011: - added support for xfnr==-1 (%RECNO)
+
+# NOTES:
+#	- the [.interp.*] routines assume that x increases strictly monotonically;
+#	  for utilities dealing with non-monotonic data, use [.nminterp.*] instead
+#	- the [.interp.*] routines are written to work on multiple buffers,
+#	  rather than just @ants_; this implies that data created by IS_init()
+#	  must be stored separately for each buffer
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$IS_opts			string of allowed options
+#	$IS_optsUsage		usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_usage() 	mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_usage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_init(br,f,xf)			init interpolation of field f
+#		br					data buffer reference
+#		idr					init-data reference
+#       f					field number
+#		xf					x field number or -1 for %RECNO
+#       <ret val>           none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_init($$$$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_interpolate(bR,idR,xf,xv,xi,f)	interpolate field f
+#		bR							data buffer reference
+#		idR							init-data reference
+#		xf							x field or -1 for %RECNO
+#		xv							x value
+#		xi							index of last record with x-value <= x
+#		f							field number to interpolate
+#		<ret val>					interpolated value
+#
+# NB:
+#	- handle f == xf
+#	- return NaN if any of the y values required is NaN
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_interpolate($$$$$$)
+{
+	my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+	return $xv if ($xf == $f);
+	return nan unless (numberp($bR->[$xi][$f]) && numberp($bR->[$xi+1][$f]));
+
+	my($sc) = ($xf < 0)
+			? ($xv - $bR->[$xi][$xf]) / ($bR->[$xi+1][$xf] - $bR->[$xi][$xf]);
+			: ($xv - $bR->[$xi][$xf]) / ($bR->[$xi+1][$xf] - $bR->[$xi][$xf]);
+	return $bR->[$xi][$f] + $sc * ($bR->[$xi+1][$f] - $bR->[$xi][$f]);
+}
+
+#----------------------------------------------------------------------
+
+1;
diff --git a/.interp.nnbr b/.interp.nnbr
new file mode 100644
--- /dev/null
+++ b/.interp.nnbr
@@ -0,0 +1,37 @@
+#======================================================================
+#                    . I N T E R P . N N B R 
+#                    doc: Wed Nov 22 21:01:09 2000
+#                    dlm: Wed Aug  3 11:48:59 2011
+#                    (c) 2000 A.M. Thurnherr
+#                    uE-Info: 32 0 NIL 0 0 72 8 2 4 NIL ofnI
+#======================================================================
+
+# nearest neighbor resampling
+
+# HISTORY:
+# 	Nov 22, 2000: - adapted from [.interp.linear]
+#	Apr  3, 2004: - added nominal -x
+#	Jan 12, 2006: - renamed from [.interp.subsample]
+#	Jul 28, 2006: - added xf to ISInit() args
+#	Aug 22, 2006: - adapted to work with [match]
+#   Aug  5, 2008: - added idr param to IS_init()
+#	Aug  3, 2011: - removed -x, which is now handled by [resample]
+
+# see [.interp.linear] for documentation of interface
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+sub IS_usage() {}
+sub IS_init($$$$) {}
+
+sub IS_interpolate($$$$$$)
+{
+	my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+
+	return $bR->[$xi+1][$f]
+		if ($bR->[$xi+1][$xf] - $xv < $xv - $bR->[$xi][$xf]);
+	return $bR->[$xi][$f];
+}
+
+1;
diff --git a/.interp.poly b/.interp.poly
new file mode 100644
--- /dev/null
+++ b/.interp.poly
@@ -0,0 +1,94 @@
+#======================================================================
+#                    . I N T E R P . P O L Y 
+#                    doc: Thu Nov 23 21:30:25 2000
+#                    dlm: Tue Aug  5 14:20:19 2008
+#                    (c) 2000 A.M. Thurnherr
+#                    uE-Info: 58 0 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+# polynomial interpolation
+
+# HISTORY:
+# 	Nov 23, 2000: - created
+#	Nov 24, 2000: - added -c)enter x-val in window
+#	Dec 13, 2002: - replaced -n by -o
+#	Oct 22, 2003: - changed -e to -r (conflicts with -e option of [resample])
+#	Jul 28, 2006: - removed debugging code
+#				  - added xf to ISInit() args
+#				  - BUG: sufficient-input data test was buggy
+#   Aug 22, 2006: - adapted to work with [match]
+#	Sep 11, 2006: - BUG: sufficient-input ermesg was wrong
+#   Aug  5, 2008: - added idr param to IS_init()
+
+# see [.interp.linear] for documentation of interface
+
+# NOTES:
+#	- $opt_o <poly order> is required
+#	- $opt_r is used to return the error estimates instead of the
+#	  interpolated values
+#	- -o 0 is almost like using -s subsample, the only difference
+#	  being that subsample returns the real x values while -o 0
+#	  returns the subsampled x values
+#	- on -c the window for each interpolation estimate is centered
+#	  around the interpolated x value; this implies window-sliding
+#	  between tabulated x values resulting in functional disconti-
+#	  nuities
+
+require	"$ANTS/polint.pl";
+
+$IS_opts = "co:r";
+$IS_optsUsage = "[e-r)ror estimates] [-c)enter x-val] -o <poly order>";
+
+sub IS_usage()
+{
+	die("$0 (.interp.poly): ERROR! -o required\n")
+		unless (defined($opt_o) && $opt_o >= 0);
+	&antsInfo("WARNING: -c generates discontinuous function")
+		if ($opt_c);
+	&antsInfo("WARNING: high-order polynomial interpolation deprecated")
+		if ($opt_o > 5);
+	&antsInfo("order-1 polynomial; -s linear is more efficient")
+		if ($opt_o == 1);
+	&antsInfo("giving error estimates") if ($opt_r);
+}
+
+sub IS_init($$$$)
+{
+	my($bR,$idR,$f,$xf) = @_;
+	die(sprintf("$0 (.interp.poly): ERROR! " .
+				"not enough data points (%d) for order-$opt_o poly\n",
+				 $#{$bR}+1))
+		if ($opt_o > $#{$bR});
+}
+
+sub findWindow($$$$$)							# find window of size n=opt_o+1
+{
+	my($bR,$f,$v,$i,$n) = @_;
+	my($cn);									# current win size
+
+	return ($bR->[$i+1][$f] - $v < $v - $bR->[$i][$f]) ? $i+1 : $i
+		if ($n == 1);							# nearest neighbor
+		
+	if ($opt_c) {								# center window around x value
+		for ($cn=2; $cn<$n; $cn++) {			# grow window
+			$i-- if (($i>0 && $bR->[$i+$cn-1][$f] - $v >= $v - $bR->[$i][$f]) ||
+					 ($i+$cn-1 == $#{$bR}));
+		}
+	} else {
+		$i -= int(($n-1)/2);						# ideally
+		$i  = 0 if ($i < 0);						# off leftmost column
+		$i  = $#{$bR}+1 - $n if ($i+$n > $#{$bR}+1);# off rightmost column
+	}
+
+	return $i;
+}
+
+sub IS_interpolate($$$$$$)
+{
+	my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+	return $xv if ($xf == $f);
+
+	$xi = &findWindow($bR,$xf,$xv,$xi,$opt_o+1);
+	my($y,$dy) = &polint($bR,$xf,$xv,$xi,$opt_o+1,$f);
+	return $opt_r ? $dy : $y;
+}
diff --git a/.interp.spline b/.interp.spline
new file mode 100644
--- /dev/null
+++ b/.interp.spline
@@ -0,0 +1,125 @@
+#======================================================================
+#                    . I N T E R P . S P L I N E 
+#                    doc: Wed Nov 22 21:01:09 2000
+#                    dlm: Tue Aug  5 14:20:39 2008
+#                    (c) 2000 A.M. Thurnherr
+#                    uE-Info: 31 0 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+# spline interpolation
+
+# HISTORY:
+# 	May 27, 2003: - adapted from [.interp.linear]
+#	Jul 19, 2004: - BUG: i instead of $i --- 'orrible!
+#	May 15, 2006: - BUG: assertion did not work on -e
+#	Jul 28, 2006: - Version 3.3 [HISTORY]
+#				  - made 100% compatible with [.interp.spline]
+#				  - added xf to ISInit() args
+#	Jul 30, 2006: - continued debugging
+#   Aug 22, 2006: - adapted to work with [match]
+#   Aug  5, 2008: - added idr param to IS_init()
+
+# see [.interp.linear] for documentation of interface
+
+# EXAMPLE:
+#	Matlab: x = 0:10;  y = sin(x); xx = 0:.25:10;
+#		yy = spline(x,y,xx); plot(x,y,'o',xx,yy)
+# 	gnuplot/ANTS:
+#		set style function points;
+#		set style data lines;
+#		plot [0:10] sin(x), \
+#		'<Cat -f x:0,1,10 | list x y=sin($x) | resample -s spline x \#0-10:0.25'
+                        
+$IS_opts = "B:";
+$IS_optsUsage = "[-B)oundary cond <1st/last>]";
+
+#----------------------------------------------------------------------
+
+my($yp1,$ypn);
+
+sub IS_usage()
+{
+	if (defined($opt_B)) {								# NR p.115
+		($yp1,$ypn) = split('/',$opt_B);
+		croak("$0: can't decode -B $opt_B\n")
+			unless (numberp($yp1) && numberp($ypn));
+	}
+}
+
+#----------------------------------------------------------------------
+
+sub IS_init($$$$) {
+	my($bR,$idR,$f,$xf) = @_;
+	my($i,$k,$p,$qn,$sig,$un,@u);
+	my($n) = scalar(@{$bR});
+
+	if (defined($opt_B)) {								# handle boundary cond
+		$idR->[1][$f] = -0.5;
+		$u[1] = (3/($bR->[1][$xf]-$bR->[0][$xf])) *
+				    (($bR->[1][$f]-$bR->[0][$f]) /
+					 ($bR->[1][$xf]-$bR->[0][$xf]) - $yp1);
+	} else {
+		$idR->[1][$f] = $u[1] = 0;
+	}
+
+	for ($i=2; $i<=$n-1; $i++) {
+		$sig = ($bR->[$i-1][$xf]-$bR->[$i-2][$xf]) /
+				($bR->[$i][$xf]-$bR->[$i-2][$xf]);
+		$p = $sig*$idR->[$i-1][$f] + 2;
+		$idR->[$i][$f] = ($sig-1)/$p;
+		$u[$i] = ($bR->[$i][$f]-$bR->[$i-1][$f]) /
+					($bR->[$i][$xf]-$bR->[$i-1][$xf])
+						- ($bR->[$i-1][$f]-$bR->[$i-2][$f]) /
+							($bR->[$i-1][$xf]-$bR->[$i-2][$xf]);
+		$u[$i] = (6*$u[$i]/($bR->[$i][$xf]-$bR->[$i-2][$xf]) -
+					$sig * $u[$i-1]) / $p;
+	}
+	
+	if (defined($opt_B)) {
+		$qn = 0.5;
+		$un = (3/($bR->[$n-1][$xf]-$bR->[$n-2][$xf])) *
+				($ypn-($bR->[$n-1][$f]-$bR->[$n-2][$f]) /
+					($bR->[$n-1][$xf]-$bR->[$n-2][$xf]));
+    } else {
+		$qn = $un = 0;
+	}
+	
+	$idR->[$n][$f] = ($un-$qn*$u[$n-1]) / ($qn*$idR->[$n-1][$f]+1);
+	for ($k=$n-1;$k>=1;$k--) {
+		$idR->[$k][$f] = $idR->[$k][$f] * $idR->[$k+1][$f] + $u[$k];
+	}
+}
+
+#----------------------------------------------------------------------
+
+sub IS_interpolate($$$$$$)
+{
+	my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+
+	return $xv if ($xf == $f);
+
+	return $bR->[$xi][$f]						# edge values are ok
+		if equal($xv,$bR->[$xi][$xf]);
+	return $bR->[$xi+1][$f]
+		if equal($xv,$bR->[$xi+1][$xf]);
+
+	return nan unless (numberp($bR->[$xi][$f]) &&
+					   numberp($bR->[$xi+1][$f]));
+
+	croak("$0: assertion $bR->[$xi+1][$xf] >= $xv >=  $bR->[$xi][$xf] failed")
+		unless (defined($opt_e) || $bR->[$xi+1][$xf] >= $xv && $xv >= $bR->[$xi][$xf]);
+		
+	my($h) = $bR->[$xi+1][$xf] - $bR->[$xi][$xf];
+	croak("$0: assertion #2 failed") unless ($h > 0);
+	my($a) = ($bR->[$xi+1][$xf] - $xv) / $h;
+	my($b) = ($xv - $bR->[$xi][$xf]) / $h;
+
+	return $a*$bR->[$xi][$f] + $b*$bR->[$xi+1][$f] +
+			(($a*$a*$a-$a) * $idR->[$xi+1][$f] +
+			 ($b*$b*$b-$b) * $idR->[$xi+2][$f]) *
+			 	($h*$h)/6;
+}
+
+#----------------------------------------------------------------------
+
+1;
diff --git a/.isopycnal_TS.gamma_n b/.isopycnal_TS.gamma_n
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.gamma_n
@@ -0,0 +1,36 @@
+#======================================================================
+#                    . I S O P Y C N A L _ T S . G A M M A _ N 
+#                    doc: Tue Dec 13 21:50:18 2005
+#                    dlm: Mon Dec 19 13:04:15 2005
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 27 30 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Dec 14, 2005: - created
+#	Dec 19, 2005: - finalized
+
+# NOTES:
+#	- requires %lat/%lon PARAMs
+
+require "$ANTS/libgamma.pl";				# load equation of state
+
+unless (defined($P{ITS})) {
+	&antsInfo("using default %ITS=90");
+	&antsAddParams(ITS,90);
+}
+$gamma::temp_scale = $P{ITS};
+
+croak("$0: need %lat/%lon\n")
+	unless defined($P{lat}) && defined($P{lon});
+croak("$0: need %press\n")
+	unless defined($P{press});
+
+sub density($$)
+{
+	my($S,$T) = @_;
+	my($gamma) = gamma::gamma_n($S,$T,$P{press},$P{lat},$P{lon});
+	return $gamma > 0 ? $gamma : undef;
+}
+
+1;
diff --git a/.isopycnal_TS.sigma0 b/.isopycnal_TS.sigma0
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.sigma0
@@ -0,0 +1,29 @@
+#======================================================================
+#                    . I S O P Y C N A L _ T S . S I G M A 0 
+#                    doc: Mon Dec 19 12:28:01 2005
+#                    dlm: Mon Dec 19 13:27:24 2005
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 19 56 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl";				# load equation of state
+
+unless (defined($P{ITS})) {
+	&antsInfo("using default %ITS=90");
+	&antsAddParams(ITS,90);
+}
+
+&antsInfo("Warning: ignoring non-zero %press = $P{press}")
+	if ($P{press} != 0);
+$P{press} = 0;
+
+sub density($$)
+{
+	my($S,$T) = @_;
+	return sigma($S,$T,$P{press},0);
+}
+
+1;
diff --git a/.isopycnal_TS.sigma2 b/.isopycnal_TS.sigma2
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.sigma2
@@ -0,0 +1,29 @@
+#======================================================================
+#                    . I S O P Y C N A L _ T S . S I G M A 2 
+#                    doc: Mon Dec 19 12:28:01 2005
+#                    dlm: Mon Dec 19 13:29:23 2005
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 26 32 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl";				# load equation of state
+
+unless (defined($P{ITS})) {
+	&antsInfo("using default %ITS=90");
+	&antsAddParams(ITS,90);
+}
+
+&antsInfo("Warning: ignoring non-zero %press = $P{press}")
+	if ($P{press} != 0);
+$P{press} = 0;
+
+sub density($$)
+{
+	my($S,$T) = @_;
+	return sigma($S,$T,$P{press},2000);
+}
+
+1;
diff --git a/.isopycnal_TS.sigma4 b/.isopycnal_TS.sigma4
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.sigma4
@@ -0,0 +1,29 @@
+#======================================================================
+#                    . I S O P Y C N A L _ T S . S I G M A 4 
+#                    doc: Mon Dec 19 12:28:01 2005
+#                    dlm: Wed Sep 16 10:19:10 2009
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 26 34 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Sep 16, 2009: created
+
+require "$ANTS/libEOS83.pl";				# load equation of state
+
+unless (defined($P{ITS})) {
+	&antsInfo("using default %ITS=90");
+	&antsAddParams(ITS,90);
+}
+
+&antsInfo("Warning: ignoring non-zero %press = $P{press}")
+	if ($P{press} != 0);
+$P{press} = 0;
+
+sub density($$)
+{
+	my($S,$T) = @_;
+	return sigma($S,$T,$P{press},4000);
+}
+
+1;
diff --git a/.isopycnal_TS.theta0 b/.isopycnal_TS.theta0
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.theta0
@@ -0,0 +1,22 @@
+#======================================================================
+#                    . I S O P Y C N A L _ T S . T H E T A 0 
+#                    doc: Mon Dec 19 09:43:02 2005
+#                    dlm: Mon Dec 19 13:38:50 2005
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 21 42 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# calculate in-situ temperature at surface given theta2
+
+# HISTORY:
+#	Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl";				# load equation of state
+
+$temperature_fname = 'theta0';
+
+sub temperature($$)
+{
+	my($salin,$potemp) = @_;
+	return &temp($salin,$potemp,$P{press},0);
+}
diff --git a/.isopycnal_TS.theta2 b/.isopycnal_TS.theta2
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.theta2
@@ -0,0 +1,22 @@
+#======================================================================
+#                    . I S O P Y C N A L _ T S . T H E T A 2 
+#                    doc: Mon Dec 19 09:43:02 2005
+#                    dlm: Mon Dec 19 13:43:18 2005
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 21 32 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# calculate in-situ temperature at surface given theta2
+
+# HISTORY:
+#	Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl";				# load equation of state
+
+$temperature_fname = 'theta2';
+
+sub temperature($$)
+{
+	my($salin,$potemp) = @_;
+	return &temp($salin,$potemp,$P{press},2000);
+}
diff --git a/.lmfit.exp b/.lmfit.exp
new file mode 100644
--- /dev/null
+++ b/.lmfit.exp
@@ -0,0 +1,117 @@
+#======================================================================
+#                    . L M F I T . E X P 
+#                    doc: Wed Feb 24 09:40:06 1999
+#                    dlm: Fri Jul 28 13:40:56 2006
+#                    (c) 1999 A.M. Thurnherr
+#                    uE-Info: 30 41 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+#	- a number of global variables to be set during loading
+#	- a number of subs to perform admin tasks (usage, init, ...)
+#	- a sub to evaluate the model function which is to be fitted using
+#	  a number of pararams which are all stored in @A (beginning at
+#	  A[1]!!!). You also need to return the partial derivatives of
+#	  the model function wrt all params.
+#	- the interface is documented between +++++++ lines
+
+# fit exponential A[3]+A[2]*exp(A[1]*x) to data
+#
+# NOTES:
+#	- initial parameter estimates are crucial
+#	- there is currently no heuristics
+
+# HISTORY:
+#	Mar 11, 1999: - created from [./.mfit.poly] & [./.mfit.gauss]
+#	Jul 31, 1999: - typecheck usage
+#   Mar 17, 2001: - param->arg
+#	Jan 16, 2006: - added notes
+#	Jul 28, 2006: - Version 3.3 [HISTORY]
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$modelOpts			string of allowed options
+#	$modelOptsUsage		usage information string for options
+#	$modelMinArgs		min # of arguments of model
+#	$modelArgsUsage		usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+#	$modelNFit			number of params to fit in model
+#	@nameA				symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "";
+$modelOptsUsage = "";
+$modelMinArgs = 0;
+$modelArgsUsage = "[exp [mul [add guess]]]";
+$modelNFit = 3;
+$nameA[1] = "exp";
+$nameA[2] = "mul";		
+$nameA[3] = "add";	
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
+#					filenames after model arguments; this usually sets
+#					@A (the model parameters) but these can later be
+#					calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+	$A[1] = nan; $A[2] = nan; $A[3] = nan;				# usage
+	$A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	$A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	$A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit()		initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+	$A[1] = 1 unless (numberp($A[1]));
+	$A[2] = 1 unless (numberp($A[2]));
+	$A[3] = 0 unless (numberp($A[3]));
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda)	evaluate polynom and derivatives
+#		x					x value (NOT xfnr)
+#		A					reference to @A
+#		dyda				reference to array for partial derivatives
+#							(wrt individaul params in @A)
+#		<ret val>			y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)					# y = A[3]+A[2]*exp(A[1]*x)
+{
+	my($x,$AR,$dydaR) = @_;
+	my($e) = exp($AR->[1]*$x);
+
+	$dydaR->[1] = $AR->[2]*$x*$e;
+	$dydaR->[2] = $e;
+	$dydaR->[3] = 1;					# partial derivatives
+
+	return $AR->[3] + $AR->[2]*$e;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup()	cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+}
diff --git a/.lmfit.gauss b/.lmfit.gauss
new file mode 100644
--- /dev/null
+++ b/.lmfit.gauss
@@ -0,0 +1,176 @@
+#======================================================================
+#                    . L M F I T . G A U S S 
+#                    doc: Wed Feb 24 09:40:06 1999
+#                    dlm: Fri Jul 28 13:32:35 2006
+#                    (c) 1999 A.M. Thurnherr
+#                    uE-Info: 35 51 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+#	- a number of global variables to be set during loading
+#	- a number of subs to perform admin tasks (usage, init, ...)
+#	- a sub to evaluate the model function which is to be fitted using
+#	  a number of pararams which are all stored in @A (beginning at
+#	  A[1]!!!). You also need to return the partial derivatives of
+#	  the model function wrt all params.
+#	- the interface is documented between +++++++ lines
+
+# Gauss data model (i.e. fit Gaussian curve)
+# NB: - fitting is rather sensitive to the input parameters, thus
+# 	    a heuristic has been added to guess them (by setting them
+# 	    to NaN)
+# 	  - another fickle parameter is the y-offset (zero line); thus
+# 	    a heuristics has been added for this one as well
+#	  - the parameters are peak, mean, standard deviation
+
+# HISTORY:
+#	Feb 24, 1999: - created together with [./cfit]
+#	Feb 25, 1999: - cosmetic changes
+#	Jul 31, 1999: - parameter typecheck
+#	Oct 04, 1999: - changed param names
+#	Oct 05, 1999: - improved heuristics
+#				  - changed e-scale to sigma
+#   Mar 17, 2001: - param->arg
+#	Jul 28, 2006: - Version 3.3 [HISTORY]; &isnan()
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$modelOpts			string of allowed options
+#	$modelOptsUsage		usage information string for options
+#	$modelMinArgs		min # of arguments of model
+#	$modelArgsUsage		usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+#	$modelNFit			number of params to fit in model
+#	@nameA				symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "y";
+$modelOptsUsage = "[-y)shift]";
+$modelMinArgs = 0;
+$modelArgsUsage = "[peak guess [mean guess [sigma guess]]]";
+$modelNFit = 3;			
+$nameA[1] = "peak";
+$nameA[2] = "mean";		
+$nameA[3] = "sigma";	
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
+#					filenames after model arguments; this usually sets
+#					@A (the model parameters) but these can later be
+#					calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+	$A[1] = nan; $A[2] = nan; $A[3] = nan;				# usage
+	$A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	$A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	$A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit()		initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+	my($i,$j,$ymin,$ymax,$xatymax);
+
+#	--------------------------------------------------
+#	heuristics for initial model param values
+#	--------------------------------------------------
+
+	$ymin = 1e33, $ymax = -1e33, $xatymax = 0;
+	for ($i=0; $i<=$#ants_; $i++) {
+		next if ($antsFlagged[$i]);
+		$ymin = $ants_[$i][$yfnr]
+			if ($ants_[$i][$yfnr] < $ymin);
+		$ymax = $ants_[$i][$yfnr], $xatymax = $ants_[$i][$xfnr]
+			if ($ants_[$i][$yfnr] > $ymax);
+	}
+	$A[1] = $ymax - $ymin if isnan($A[1]);			# peak guess
+	$A[2] = $xatymax if isnan($A[2]);				# mean guess
+	if (isnan($A[3])) {								# sigma guess
+		for ($i=1;
+			 $i<=$#ants_ && !$antsFlagged[$i]
+						 && $ants_[$i][$yfnr]-$ymin<0.36*$A[1];
+			 $i++) {}				   
+		for ($j=$#ants_;
+			 $j>=1 && !$antsFlagged[$i]
+				   && $ants_[$j][$yfnr]-$ymin < 0.36*$A[1];
+			 $j--) {}
+		$A[3] = abs($ants_[$i][$xfnr]-$ants_[$j][$xfnr]) / 2.0;
+		$A[3] *= 0.71;								# scale by 1/sqrt(2)
+		if ($A[3] == 0) {
+			&antsInfo("$model: sigma heuristic failed (set to 1)!");
+			$A[3] = 1;
+	    }
+	}
+
+#	--------------------------------------------------
+#	y shift (-y option)
+#	--------------------------------------------------
+
+	if ($opt_y) {
+		for ($i=1; $i<=$#ants_; $i++) {
+			$ants_[$i][$yfnr] -= $ymin;
+		}
+	}
+
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda)	evaluate sum of Gaussians (p.528) at x
+#		x					x value (NOT xfnr)
+#		A					reference to @A
+#		dyda				reference to array for partial derivatives
+#							(wrt individaul params in @A)
+#		<ret val>			y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+	my($x,$AR,$dydaR) = @_;
+	my($i,$fac,$ex,$arg,$sqrt2sig);
+	my($y) = 0;
+
+	for ($i=1; $i < $#{$AR}; $i+=3) {
+		$sqrt2sig = (1.4142135623731*$AR->[$i+2]);
+		$arg = ($x - $AR->[$i+1]) / $sqrt2sig;
+		$ex  = exp(-$arg*$arg);
+		$fac = $AR->[$i] * $ex * 2*$arg;
+		$y += $AR->[$i] * $ex;
+		
+		$dydaR->[$i]   = $ex;
+		$dydaR->[$i+1] = $fac / $sqrt2sig;
+		$dydaR->[$i+2] = $fac * $arg / $sqrt2sig;
+	}
+	return $y;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup()	cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+	if ($opt_y) {
+		$A[1] += $ymin;
+		for ($i=1; $i<=$#ants_; $i++) {
+			$ants_[$i][$yfnr] += $ymin;
+		}
+	}
+}
diff --git a/.lmfit.normal b/.lmfit.normal
new file mode 100644
--- /dev/null
+++ b/.lmfit.normal
@@ -0,0 +1,175 @@
+#======================================================================
+#                    . L M F I T . N O R M A L 
+#                    doc: Wed Feb 24 09:40:06 1999
+#                    dlm: Fri Jul 28 13:35:24 2006
+#                    (c) 1999 A.M. Thurnherr
+#                    uE-Info: 34 51 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+#	- a number of global variables to be set during loading
+#	- a number of subs to perform admin tasks (usage, init, ...)
+#	- a sub to evaluate the model function which is to be fitted using
+#	  a number of pararams which are all stored in @A (beginning at
+#	  A[1]!!!). You also need to return the partial derivatives of
+#	  the model function wrt all params.
+#	- the interface is documented between +++++++ lines
+
+# check if a given distribution is normal
+# NB:
+#	- fitting is based on gauss curve fitting [.lmfit.gauss]
+#	- heuristics are taken from there and scaled for the normal
+#	  parameter choices
+#	- simplified, e.g. y-shift is removed (does not make sense for
+#	  distribution)
+#	- added chi^2 significance testing to &modelCleanup() on -x
+
+# HISTORY:
+#	Oct 04, 1999: - created from [.lmfit.gauss]
+#	Oct 05, 1999: - added chi^2 significance test
+#				  - removed -y
+#				  - improved heuristics
+#   Mar 17, 2001: - param->arg
+#	Jul 28, 2006: - Version 3.3 [HISTORY]; &isnan()
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$modelOpts			string of allowed options
+#	$modelOptsUsage		usage information string for options
+#	$modelMinArgs		min # of arguments of model
+#	$modelArgsUsage		usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+#	$modelNFit			number of params to fit in model
+#	@nameA				symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "x";
+$modelOptsUsage = "[-x chi^2 test]";
+$modelMinArgs = 0;
+$modelArgsUsage = "[area guess [mean guess [sigma guess]]]";
+$modelNFit = 3;			
+$nameA[1] = "area";
+$nameA[2] = "mean";		
+$nameA[3] = "sigma";	
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
+#					filenames after model arguments; this usually sets
+#					@A (the model parameters) but these can later be
+#					calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+	$A[1] = nan; $A[2] = nan; $A[3] = nan;				# usage
+	$A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	$A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	$A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit()		initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+	my($i,$j,$ymin,$ymax,$xatymax);
+
+#	--------------------------------------------------
+#	heuristics for initial model param values
+#	--------------------------------------------------
+
+	$ymin = 1e33, $ymax = -1e33, $xatymax = 0;
+	for ($i=0; $i<=$#ants_; $i++) {
+		next if ($antsFlagged[$i]);
+		$ymin = $ants_[$i][$yfnr]
+			if ($ants_[$i][$yfnr] < $ymin);
+		$ymax = $ants_[$i][$yfnr], $xatymax = $ants_[$i][$xfnr]
+			if ($ants_[$i][$yfnr] > $ymax);
+	}
+	$A[1] = $ymax - $ymin if isnan($A[1]);			# peak guess
+	$A[2] = $xatymax if isnan($A[2]);				# mean guess
+	if (isnan($A[3])) {								# e-scale guess
+		for ($i=1;
+			 $i<=$#ants_ && !$antsFlagged[$i]
+						 && $ants_[$i][$yfnr]-$ymin<0.36*$A[1];
+			 $i++) {}
+		for ($j=$#ants_;
+			 $j>=1 && !$antsFlagged[$i]
+				   && $ants_[$j][$yfnr]-$ymin < 0.36*$A[1];
+			 $j--) {}
+		$A[3] = abs($ants_[$i][$xfnr]-$ants_[$j][$xfnr]) / 2.0;
+		if ($A[3] == 0.0) {
+			&antsInfo("$model: sigma heuristic failed (set to 1)!");
+			$A[3] = 1.0;
+	    }
+	}
+
+	$A[1] *= 1.77 * $A[3]							# gauss -> normal
+		unless (isnan($A[1]) || isnan($A[3]));
+	$A[3] *= 0.71 unless isnan($A[3]);
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda)	evaluate Normal distribution curve at x
+#		x					x value (NOT xfnr)
+#		A					reference to @A
+#		dyda				reference to array for partial derivatives
+#							(wrt individual params in @A)
+#		<ret val>			y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+	my($x,$AR,$dydaR) = @_;
+
+	my($peak) = $AR->[1] / (2.506628274631 * $AR->[3]);
+	my($dx  ) = $x - $AR->[2];
+	my($sig2) = $AR->[3] * $AR->[3];
+	my($expo) = exp(-$dx*$dx/(2*$sig2));
+	my($norm) = $peak * $expo;
+
+	if (defined($dydaR)) {
+		$dydaR->[1] = $norm / $AR->[1];
+		$dydaR->[2] = $norm * $dx / $sig2;
+	    $dydaR->[3] = $norm/$AR->[3] * ($dx*$dx/$sig2 - 1);
+	}
+
+	return $norm;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup()	cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+	return unless ($opt_x);
+
+	require "$ANTS/libfuns.pl";
+	my($chisq) = 0;
+	my($nval,$prob,$sign);
+
+	for ($i=0; $i<=$#ants_; $i++) {
+		next if ($antsFlagged[$i]);
+#		next if ($ants_[$i][$yfnr] <= 1);	# IGNORE TAIL HEURISTICS
+		$nval = &modelEvaluate($ants_[$i][$xfnr],\@A);
+		$chisq += ($ants_[$i][$yfnr] - $nval)**2 / $nval;
+	}
+	$prob = &gammq(($ndata-3)/2,$chisq/2);
+	$sign = int($prob*100);
+	&antsInfo("$model: normal-distr. hypothesis disproved at $sign%% sign. level");
+}
diff --git a/.lmfit.poly b/.lmfit.poly
new file mode 100644
--- /dev/null
+++ b/.lmfit.poly
@@ -0,0 +1,129 @@
+#======================================================================
+#                    . L M F I T . P O L Y 
+#                    doc: Wed Feb 24 09:40:06 1999
+#                    dlm: Fri Jul 28 13:35:50 2006
+#                    (c) 1999 A.M. Thurnherr
+#                    uE-Info: 28 41 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+#	- a number of global variables to be set during loading
+#	- a number of subs to perform admin tasks (usage, init, ...)
+#	- a sub to evaluate the model function which is to be fitted using
+#	  a number of pararams which are all stored in @A (beginning at
+#	  A[1]!!!). You also need to return the partial derivatives of
+#	  the model function wrt all params.
+#	- the interface is documented between +++++++ lines
+
+# fit polynomial (sum of A_i x^i) to data
+# NB:
+
+# HISTORY:
+#	Feb 25, 1999: - created
+#	Mar 14, 1999: - cosmetic changes
+#	Jul 31, 1999: - argument typechecking
+#   Mar 17, 2001: - param->arg
+#	Jan 12, 2006: - specify order with -o as in [.lsfit.poly]
+#	Jul 28, 2006: - Version 3.3 [HISTORY]
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$modelOpts			string of allowed options
+#	$modelOptsUsage		usage information string for options
+#	$modelMinArgs		min # of arguments of model
+#	$modelArgsUsage		usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+#	$modelNFit			number of params to fit in model
+#	@nameA				symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "o:";
+$modelOptsUsage = "-o)rder <n>";
+$modelMinArgs = 0;
+$modelArgsUsage = "[c0 [c1 [...]]]";
+
+&antsInfo("non-linear method deprecated; use `lsfit' instead");
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
+#					filenames after model arguments; this usually sets
+#					@A (the model parameters) but these can later be
+#					calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+	my($c);
+	
+	die("$0 (.lmfit.poly): ERROR! -o required\n")	# order of polynomial
+		unless (defined($opt_o) && $opt_o >= 0);
+	$modelNFit = &antsCardOpt($opt_o)+1;			
+	
+	for ($c=0; $c<$modelNFit; $c++) {				# init coefficients
+		if ($#ARGV >= 0 && ! -r $ARGV[0]) {
+			$A[$c+1] = &antsFloatArg();
+		} else {
+			$A[$c+1] = nan;
+		}
+		$nameA[$c+1] = "c$c";						# and names
+	}
+	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit()		initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+	my($c);
+
+	for ($c=0; $c<$modelNFit; $c++) {				# init coefficients
+		$A[$c+1] = 10**-$c unless (numberp($A[$c+1]));
+	}
+
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda)	evaluate polynomial and derivatives
+#		x					x value (NOT xfnr)
+#		A					reference to @A
+#		dyda				reference to array for partial derivatives
+#							(wrt individaul params in @A)
+#		<ret val>			y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+	my($x,$AR,$dydaR) = @_;
+	my($i);
+	my($pow) = 1;
+	my($y) = 0;
+
+	for ($i=1; $i<=$modelNFit; $i++) {
+		$y += $AR->[$i]*$pow;
+		$dydaR->[$i] = $pow;
+		$pow *= $x;
+	}
+	return $y;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup()	cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+}
diff --git a/.lmfit.sqrt b/.lmfit.sqrt
new file mode 100644
--- /dev/null
+++ b/.lmfit.sqrt
@@ -0,0 +1,105 @@
+#======================================================================
+#                    . L M F I T . S Q R T 
+#                    doc: Fri Oct 10 15:50:42 2014
+#                    dlm: Sat Oct 11 09:30:48 2014
+#                    (c) 2014 A.M. Thurnherr
+#                    uE-Info: 27 32 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+#	- a number of global variables to be set during loading
+#	- a number of subs to perform admin tasks (usage, init, ...)
+#	- a sub to evaluate the model function which is to be fitted using
+#	  a number of pararams which are all stored in @A (beginning at
+#	  A[1]!!!). You also need to return the partial derivatives of
+#	  the model function wrt all params.
+#	- the interface is documented between +++++++ lines
+
+# fit square root A[1]*sqrt(x) to data
+#
+# NOTES:
+#	- initial parameter estimates may be important
+#	- there is currently no heuristics
+
+# HISTORY:
+#	Oct 10, 2014: - created from [.lmfit.exp]
+#	Oct 11, 2014: - made it work
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$modelOpts			string of allowed options
+#	$modelOptsUsage		usage information string for options
+#	$modelMinArgs		min # of arguments of model
+#	$modelArgsUsage		usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+#	$modelNFit			number of params to fit in model
+#	@nameA				symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "";
+$modelOptsUsage = "";
+$modelMinArgs = 0;
+$modelArgsUsage = "[scale guess]";
+$modelNFit = 1;
+$nameA[1] = "scale";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
+#					filenames after model arguments; this usually sets
+#					@A (the model parameters) but these can later be
+#					calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+	$A[1] = nan; 			
+	$A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit()		initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+	$A[1] = 1 unless (numberp($A[1]));		# scale
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda)	evaluate polynom and derivatives
+#		x					x value (NOT xfnr)
+#		A					reference to @A
+#		dyda				reference to array for partial derivatives
+#							(wrt individaul params in @A)
+#		<ret val>			y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)					# y  = A[1]*sqrt(*x)
+{										
+	my($x,$AR,$dydaR) = @_;
+	my($v) = sqrt($x);
+
+	$dydaR->[1] = $v;					# dy/dA[1] = sqrt(x)
+	return $AR->[1]*$v;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup()	cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+}
diff --git a/.lsfit.bilin b/.lsfit.bilin
new file mode 100644
--- /dev/null
+++ b/.lsfit.bilin
@@ -0,0 +1,141 @@
+#======================================================================
+#                    . L S F I T . B I L I N 
+#                    doc: Wed Feb 24 09:40:06 1999
+#                    dlm: Fri Jul 28 13:36:36 2006
+#                    (c) 1999 A.M. Thurnherr
+#                    uE-Info: 32 41 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# linear model function to your data:
+#	- a number of global variables to be set during loading
+#	- a number of subs to perform admin tasks (usage, init, ...)
+#	- a sub to evaluate the basis funs at a given x value; each
+#	  y value must be stored in @A (beginning with
+#	  A[1]!!!).
+#	- the interface is documented between +++++++ lines
+
+# fit bi-linear function to data, i.e. y = A + B*x1 + C*x2
+
+# HISTORY:
+#	Aug 01, 1999: - adapted from [.lsfit.poly]
+#	Aug 02, 1999: - added &antsDescription()
+#	Sep 26, 1999: - cosmetics
+#				  - added vars & covars
+#	Sep 27, 1999: - changed from covar to sigmas
+#	Oct 01, 1999: - cosmetics
+#	Oct 06, 1999: - added -l
+#   Mar 17, 2001: - param->arg
+#	May  2, 2001: - updated doc
+#	Nov 17, 2005: - commented out antsDescription()
+#				  - updated stats on -p
+#	Jul 28, 2006: - Version 3.3 [HISTORY]
+
+# NOTES:
+#	- could be easily extended to multidimensional linear fit but
+#	  what's the use?
+#	- -p zeroes param C for bilinear spice method (T = A + B*sig + C*neph)
+#	  so that residual field becomes spice anomaly (linearly correlated
+#	  with neph)
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$modelOpts			string of allowed options
+#	$modelOptsUsage		usage information string for options
+#	$modelMinArgs		min # of arguments of model
+#	$modelArgsUsage		usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+#	$modelNFit			number of params to fit in model
+#	@nameA				symbolic names of model parameters
+#
+# You should call &antsDescription() for the -ct options here
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "pl:";
+$modelOptsUsage = "[-p)artial f/r] [-l)imit r_BC <min,max>]";
+$modelMinArgs = 1;
+$modelArgsUsage = "<2nd x-field>";
+$modelNFit = 3;
+$nameA[1] = "A"; $nameA[2] = "B"; $nameA[3] = "C";
+$A[1] = nan; $A[2] = nan; $A[3] = nan;
+#&antsDescription("c","bilin_$nameA[1]",
+#				 "c","bilin_$nameA[2]",
+#				 "c","bilin_$nameA[3]");
+#&antsDescription("t","bilin_sigma_$nameA[1]",
+#				 "t","bilin_sigma_$nameA[2]",
+#				 "t","bilin_sigma_$nameA[3]",
+#				 "t","bilin_ccc_$nameA[1]_$nameA[2]",
+#				 "t","bilin_ccc_$nameA[1]_$nameA[3]",
+#				 "t","bilin_ccc_$nameA[2]_$nameA[3]");
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
+#					filenames after model arguments; this usually sets
+#					@A (the model parameters) but these can later be
+#					calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+	if (defined($opt_l)) {
+		($minBC,$maxBC) = split(',',$opt_l);
+		&antsUsageError("\n>>> error with -l")
+			unless (defined($maxBC) && $maxBC > $minBC);
+	}
+	$x2fnr = &antsFieldArg();
+	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit()		initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(idx,xfnr,vals)	evaluate basis funs at x (NB: x1, x2)
+#		idx				       	current index in @ants_
+#		xfnr					field number of x field
+#		vals					reference to return values (1-relative!)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+	my($idx,$xfnr,$valsR) = @_;
+	my($i);
+
+	$valsR->[1] = 1;
+	$valsR->[2] = $ants_[$idx][$xfnr];
+	$valsR->[3] = $ants_[$idx][$x2fnr];
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup()	cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup() 
+{
+	if (defined($opt_l)) {
+		my($ccc) = $covar[2][3]/(sqrt($covar[2][2])*sqrt($covar[3][3]));
+		if ($ccc < $minBC || $ccc > $maxBC) {
+			&antsInfo("CCC B/C = %.3g out of range, fit discarded",$ccc);
+			$suppressFit = 1;
+		}
+	}
+	if ($opt_p) {
+		&antsInfo("parameter $nameA[3] = %.3g discarded",$A[3]);
+		$A[3] = 'discarded';
+		$covar[3][3] = $RMS = $sig = nan;
+	}
+}
diff --git a/.match_minimize.MAD b/.match_minimize.MAD
new file mode 100644
--- /dev/null
+++ b/.match_minimize.MAD
@@ -0,0 +1,54 @@
+#======================================================================
+#                    . M A T C H _ M I N I M I Z E . M A D 
+#                    doc: Tue Aug 22 18:31:22 2006
+#                    dlm: Wed Aug 23 23:15:13 2006
+#                    (c) 2006 A.M. Thurnherr
+#                    uE-Info: 39 26 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# sum of absolute deviations
+
+#======================================================================
+
+sub MM_usage(@)
+{
+	@MM_fname = @_;
+
+	croak("match-minimize usage: -m SAD(<field>[,...])\n")
+		unless (@MM_fname > 0);
+	for (my($fi)=0; $fi<@MM_fname; $fi++) {
+		$MM_fnr[$fi]  =  &fnr($MM_fname[$fi]);
+		$MM_wfnr[$fi] = &wfnr($MM_fname[$fi]);
+		&IS_init($MM_wfnr[$fi],$MM_cwfnr);			# prepare interpolation
+	}
+	for (my($r)=0; $r<@ants_; $r++) {
+		for (my($fi)=0; $fi<@MM_fname; $fi++) {
+			croak("$0: $MM_fname[$fi] must be numeric\n")
+				unless numberp($ants_[$r][$MM_fnr[$fi]]) &&
+					   numberp($ants_[$r][$MM_wfnr[$fi]]);
+		}
+	}
+}
+
+#======================================================================
+
+sub MM_eval(@)
+{
+	&MW_warp(@_);
+	
+	my($sad) = my($nad) = 0;
+	for (my($r)=0; $r<@ants_; $r++) {
+		my($xi) = &bSearch($MW_cwfnr,$ants_[$r][$MW_cfnr]);
+		next unless defined($xi);
+		$nad++;
+		for (my($fi)=0; $fi<@MM_fnr; $fi++) {
+			my($wv) = &IS_interpolate(\@wf_,$MW_cwfnr,$ants_[$r][$MW_cfnr],$xi,$MM_wfnr[$fi]);
+			$sad += abs($ants_[$r][$MM_fnr[$fi]] - $wv);
+		}
+	}
+	return $sad/$nad;
+}
+
+#======================================================================
+
+1;
diff --git a/.match_warp.stretch b/.match_warp.stretch
new file mode 100644
--- /dev/null
+++ b/.match_warp.stretch
@@ -0,0 +1,59 @@
+#======================================================================
+#                    . M A T C H _ W A R P . S T R E T C H 
+#                    doc: Tue Aug 22 18:31:22 2006
+#                    dlm: Wed Aug 23 23:51:25 2006
+#                    (c) 2006 A.M. Thurnherr
+#                    uE-Info: 52 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# stretch (and shift!) monotonically increasing coordinate
+
+#======================================================================
+
+sub MW_usage(@)
+{
+	croak("match-warp usage: -w stretch(<field>)\n")
+		unless (@_ == 1);
+	($MW_cfname) = @_;
+	
+	$MW_cfnr  =  fnr($MW_cfname);
+	$MW_cwfnr = wfnr($MW_cfname);
+
+	croak("$0: field $MW_cfname must be numeric\n")
+		unless numberp($ants_[0][$MW_cfnr]);
+	for (my($r)=1; $r<@ants_; $r++) {
+		croak("$0: field $MW_cfname of rec $r must be numeric ($ants_[$r][$MW_cfnr])\n")
+			unless numberp($ants_[$r][$MW_cfnr]);
+		croak("$0: field $MW_cfname must be monotonically increasing (" .
+				  "$ants_[$r][$MW_cfnr]=>$ants_[$r-1][$MW_cfnr]" .
+			  	  ")\n")
+			unless ($ants_[$r][$MW_cfnr] >= $ants_[$r-1][$MW_cfnr]);
+	}
+
+	croak("$0: warp-file field $MW_cfname must be numeric\n")
+		unless numberp($wf_[0][$MW_cwfnr]);
+	for (my($r)=1; $r<@wf_; $r++) {
+		croak("$0: warp-file field $MW_cfname of rec $r must be numeric ($wf_[$r][$MW_cwfnr])\n")
+			unless numberp($wf_[$r][$MW_cwfnr]);
+		croak("$0: warp-file field $MW_cfname must be monotonically increasing (" .
+				  "$wf_[$r][$MW_cwfnr]=>$wf_[$r-1][$MW_cwfnr]" .
+			  	  ")\n")
+			unless ($wf_[$r][$MW_cwfnr]  >= $wf_[$r-1][$MW_cwfnr]);
+		$MW_unwarped[$r] = $wf_[$r][$MW_cwfnr];
+	}
+}
+
+#======================================================================
+
+sub MW_warp(@)
+{
+	my($offset,$scale) = @_;
+
+	for (my($r)=0; $r<@wf_; $r++) {
+		$wf_[$r][$MW_cwfnr] = $MW_unwarped[$r]*$scale + $offset;
+	}
+}
+
+#======================================================================
+
+1;
diff --git a/.nminterp.linear b/.nminterp.linear
--- a/.nminterp.linear
+++ b/.nminterp.linear
@@ -1,9 +1,9 @@
 #======================================================================
 #                    . N M I N T E R P . L I N E A R 
 #                    doc: Wed Nov 22 21:01:09 2000
-#                    dlm: Thu Jan 23 13:39:13 2014
+#                    dlm: Fri Aug 30 13:26:58 2019
 #                    (c) 2000 A.M. Thurnherr
-#                    uE-Info: 19 61 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 101 9 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # linear interpolation for [fillgaps]
@@ -96,6 +96,9 @@
 			if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$nvr[$f]][$xf])) ||
 			   ($xf<0 && $cr==$nvr[$f]);
 	
+		if ($xf < 0) {
+			die("$nvr[$f] - $lvr") unless ($nvr[$f] - $lvr > 0);
+        }
 		my($sc) = ($xf >= 0)
 				? ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / ($ants_[$nvr[$f]][$xf] - $ants_[$lvr][$xf])
 				: ($cr - $lvr) / ($nvr[$f] - $lvr);
diff --git a/.nminterp.nan b/.nminterp.nan
new file mode 100644
--- /dev/null
+++ b/.nminterp.nan
@@ -0,0 +1,71 @@
+#======================================================================
+#                    . N M I N T E R P . N A N 
+#                    doc: Sun Jul  6 20:37:31 2014
+#                    dlm: Sun Jul  6 20:39:26 2014
+#                    (c) 2014 A.M. Thurnherr
+#                    uE-Info: 67 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# filling gaps with nan [fillgaps]
+
+# HISTORY:
+#	Jul  6, 2014: - created from [.nminterp.linear]
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$ISOpts				string of allowed options
+#	$ISOptsUsage		usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$ISOpts = "";
+$ISOptsUsage = "";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISUsage() 	mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub ISUsage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISInit(f,xf)				init interpolation of field f
+#       f					field number
+#		xf					x field number or -1 for %RECNO
+#       <ret val>           none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub ISInit($$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &interpolate(xf,mg,f,lvr,cr)		interpolate field f
+#		xf							x field or -1 for %RECNO
+#		mg							max allowed x gap
+#		f							field number to interpolate
+#		lvr							last record with valid y val
+#		cr							current record
+#		<ret val>					interpolated value
+#
+# NB:
+#	- handle f == xf
+#	- return undef if interpolation cannot be carried out
+#	- x VALUES MAY NOT BE MONOTONIC
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub interpolate($$$$$$)
+{
+	my($xf,$mg,$f,$lvr,$cr) = @_;
+	
+	return ($xf == $f) ? $ants_[$cr][$xf] : nan;
+} # static scope
+
+#----------------------------------------------------------------------
+
+1;
diff --git a/.nminterp.spline b/.nminterp.spline
new file mode 100644
--- /dev/null
+++ b/.nminterp.spline
@@ -0,0 +1,153 @@
+#======================================================================
+#                    . N M I N T E R P . S P L I N E 
+#                    doc: Wed Nov 22 21:01:09 2000
+#                    dlm: Fri Sep 23 16:52:31 2011
+#                    (c) 2000 A.M. Thurnherr
+#                    uE-Info: 19 50 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# spline interpolation for [fillgaps]
+
+# HISTORY:
+#	Apr  3, 2006: - created from [.interp.spline]
+#	Jul 28, 2006: - made to work
+#	Jul 30, 2006: - added max-gap support
+#				  - BUG: extended interpolation into last interval
+#				  - BUG: assertions made it incompatible with decreasing data
+#	Aug  8, 2008: - renamed
+#	Sep 23, 2011: - removed $xv from interpolate() args
+#				  - added error when xf<0 (%RECNO)
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+#	$ISOpts				string of allowed options
+#	$ISOptsUsage		usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$ISOpts = '';
+$ISOptsUsage = '';
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISUsage() 	mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub ISUsage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISInit(f,xf)				init interpolation of field f
+#       f					field number
+#		xf					x field number
+#       <ret val>           none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+my(@y2);
+
+sub ISInit($$) {
+	my($f,$xf) = @_;
+	my($k,$i,$p,$qn,$sig,$un,@u);
+
+	croak(".nminterp.spline: cannot use %RECNO as x-field (implementation restriction)\n")
+		unless ($xf >= 0);
+
+	my($pi,$ppi,$li);				# find idx of first 3 & last valid recs
+	for (my($r)=0; $r<=$#ants_; $r++) {
+		if (numberp($ants_[$r][$f])) {
+			unless (defined($ppi)) {
+				$ppi = $pi; $pi = $i; $i = $r;
+			}
+			$li = $r;
+		}
+	}
+	unless (defined($ppi)) {
+		&antsInfo("WARNING: not enough <$antsLayout[$f]> data for spline interpolation");
+		return;
+	}
+		
+
+	my($fpi) = $pi;
+	$y2[$pi][$f] = $u[$pi] = 0;
+
+	for (; $i<=$li; $i++) {
+		next unless numberp($ants_[$i][$f]);
+		
+		$sig = ($ants_[$pi][$xf]-$ants_[$ppi][$xf]) /
+				($ants_[$i][$xf]-$ants_[$ppi][$xf]);
+		$p = $sig*$y2[$pi][$f] + 2;
+		$y2[$i][$f] = ($sig-1)/$p;
+		$u[$i] = ($ants_[$i][$f]-$ants_[$pi][$f]) /
+					($ants_[$i][$xf]-$ants_[$pi][$xf])
+						- ($ants_[$pi][$f]-$ants_[$ppi][$f]) /
+							($ants_[$pi][$xf]-$ants_[$ppi][$xf]);
+		$u[$i] = (6*$u[$i]/($ants_[$i][$xf]-$ants_[$ppi][$xf]) -
+					$sig * $u[$pi]) / $p;
+
+		$ppi = $pi; $pi = $i;
+	}
+
+	$qn = $un = 0;
+
+	$y2[$li][$f] = ($un-$qn*$u[$pi]) / ($qn*$y2[$pi][$f]+1);
+	my($pk) = $li;
+	for ($k=$pi; $k>=$fpi; $k--) {
+		next unless numberp($ants_[$k][$f]);
+		$y2[$k][$f] = $y2[$k][$f] * $y2[$pk][$f] + $u[$k];
+		$pk = $k;
+	}
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &interpolate(xf,mg,f,lvr,cr)	interpolate field f
+#		xf							x field
+#		mg							max allowed x gap
+#		f							field number to interpolate
+#		lvr							last record with valid y val
+#		cr							current record
+#		<ret val>					interpolated value
+#
+# NB:
+#	- handle f == xf
+#	- return undef if interpolation cannot be carried out
+#	- x VALUES MAY NOT BE MONOTONIC
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub interpolate($$$$$$)
+{
+	my($xf,$mg,$f,$lvr,$cr) = @_;
+	
+	return $ants_[$cr][$xf] if ($xf == $f);
+
+	return $ants_[$lvr][$f]							# edge values are ok
+		if equal($ants_[$cr][$xf],$ants_[$lvr][$xf]);
+
+	my($nvr1) = $cr + 1;							# find next
+	while ($nvr1 <= $#y2 && !numberp($y2[$nvr1][$f])) { $nvr1++; }
+	return undef									# none or gap too large
+		if ($nvr1>$#ants_ || abs($ants_[$nvr1][$xf]-$ants_[$lvr][$xf])>$mg);
+	return $ants_[$nvr1][$f]						# edge values are ok
+		if equal($ants_[$cr][$xf],$ants_[$nvr1][$xf]);
+
+	my($nvr2) = $nvr1 + 1;
+	while ($nvr2 <= $#y2 && !numberp($y2[$nvr2][$f])) { $nvr2++; }
+
+	my($h) = $ants_[$nvr1][$xf] - $ants_[$lvr][$xf];
+	my($a) = ($ants_[$nvr1][$xf] - $ants_[$cr][$xf]) / $h;
+	my($b) = ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / $h;
+
+	return $a*$ants_[$lvr][$f] + $b*$ants_[$nvr1][$f] +
+			(($a*$a*$a-$a) * $y2[$nvr1][$f] +
+			 ($b*$b*$b-$b) * $y2[$nvr2][$f]) *
+			 	($h*$h)/6;
+}
+
+#----------------------------------------------------------------------
+
+1;
diff --git a/.xycontour.sigma0 b/.xycontour.sigma0
new file mode 100644
--- /dev/null
+++ b/.xycontour.sigma0
@@ -0,0 +1,47 @@
+#======================================================================
+#                    . I S O P L E T H . S I G M A 0 
+#                    doc: Tue Dec 13 21:50:18 2005
+#                    dlm: Wed Dec 14 09:40:50 2005
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 44 43 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Dec 14, 2005: - created
+
+# NOTES:
+#	- assumes T/S fields to be called 'temp' and 'salin'
+
+require "$ANTS/libEOS83.pl";				# load equation of state
+
+unless (defined($P{ITS})) {
+	&antsInfo("using default %ITS=90");
+	&antsAddParams(ITS,90);
+}
+
+sub yfield_isopleth($)						# return y-field based on x-field
+{
+	my($xf) = @_;
+	return $xf eq 'salin' ? 'temp' : 'salin';
+}
+
+sub bracket_isopleth($$$$)					# heuristically bracket y field
+{
+	my($xf,$xv,$iv,$prev_y) = @_;
+
+	if ($xf eq 'salin') {
+		return (-10,30);
+	} else {
+		return (34,37);
+	}
+}
+
+sub eval_isopleth($$$)						# evaluate function for isopleth
+{
+	my($xf,$x,$y) = @_;
+
+	return $xf eq 'salin' ? sigma($x,$y,0,0)
+						  : sigma($y,$x,0,0)
+}
+
+1;
diff --git a/.xycontour.sigma2 b/.xycontour.sigma2
new file mode 100644
--- /dev/null
+++ b/.xycontour.sigma2
@@ -0,0 +1,47 @@
+#======================================================================
+#                    . I S O P L E T H . S I G M A 2 
+#                    doc: Tue Dec 13 21:50:18 2005
+#                    dlm: Wed Dec 14 09:11:01 2005
+#                    (c) 2005 A.M. Thurnherr
+#                    uE-Info: 19 26 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Dec 13, 2005: - created
+
+# NOTES:
+#	- assumes T/S fields to be called 'temp' and 'salin'
+
+require "$ANTS/libEOS83.pl";				# load equation of state
+
+unless (defined($P{ITS})) {
+	&antsInfo("using default %ITS=90");
+	&antsAddParams(ITS,90);
+}
+
+sub yfield_isopleth($)						# return y-field based on x-field
+{
+	my($xf) = @_;
+	return $xf eq 'salin' ? 'temp' : 'salin';
+}
+
+sub bracket_isopleth($$$$)					# heuristically bracket y field
+{
+	my($xf,$xv,$iv,$prev_y) = @_;
+
+	if ($xf eq 'salin') {
+		return (-10,30);
+	} else {
+		return (34,37);
+	}
+}
+
+sub eval_isopleth($$$)						# evaluate function for isopleth
+{
+	my($xf,$x,$y) = @_;
+
+	return $xf eq 'salin' ? sigma($x,$y,0,2000)
+						  : sigma($y,$x,0,2000)
+}
+
+1;
diff --git a/OLD b/OLD
new file mode 100644
--- /dev/null
+++ b/OLD
@@ -0,0 +1,721 @@
+#======================================================================
+#                    O L D 
+#                    doc: Tue Nov 26 21:59:40 2013
+#                    dlm: Mon Jul  1 15:15:19 2019
+#                    (c) 2017 A.M. Thurnherr
+#                    uE-Info: 698 42 NIL 0 0 70 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Nov 26, 2013: - created
+#	Dec  1, 2013: - removed HDG_ROT
+#				  - added support for IMP data gaps
+#	Mar  4, 2014: - added support for DATA_SOURCE_ID
+#	May 22, 2014: - use +/-2deg to assess quality of heading offset
+#	May 23, 2014: - added $dhist_binsize, $dhist_min_pirom
+#	Jul 27, 2016: - updated variable names for consistency
+#	Jul 28, 2016: - major re-write if merging routines
+#	Jul 29, 2016: - cosmetics
+#				  - increased heading-offset resolution from 2 to 1 degrees
+#				  - BUG: inconsistent heading definition used (from old IMP with
+#						 confused coordinates)
+#				  - modified initial timelag guess (there was a bug and it is 
+#				    likely more robust based on end time rather than start time)
+#	Aug  5, 2016: - BUG: weird statement accessing LADCP_begin-1
+#				  - BUG: DSID of first ensemble was not left original
+#	Aug 22, 2016: - major changes to timelagging
+#	Aug 23, 2016: - changed semantics for removing ensembles with bad attitudes:
+#					instead of setting attitude to undef (or large pitch/roll),
+#					clearEns() is used
+#	Aug 24, 2016: - overhauled time-lagging
+#	Aug 25, 2016: - significant code cleanup
+#	Aug 26, 2016: - added _hdg_err.ps output plot
+#	Oct 13, 2016: - made hdg nan for invalid records (BUG with current versions of IMP+LADCP, IMPatchPD0)
+#	Nov 22, 2016: - added heading-offset plot
+#				  - added sensor info to plots
+#	Nov 29, 2016: - added stats to compass error plot
+#	Dec 29, 2016: - improved histogram plot
+#	Nov 16, 2017: - adapted rot_vecs() to KVM coordinates
+#				  - made sensor information optional in
+#	Nov 20, 2017: - major code cleanup
+#	Nov 22, 2017: - replaced "IMP" output in routines used by KVH by "IMU"
+#	Dec  8, 2017: - replaced remaing "IMP" output (e.g. in plot) by "IMU"
+#	May 22, 2018: - added data trimming to rot_vec
+#	May 23, 2018: - added horizontal field strength to mag calib plot
+#				  - added support for S/N 8 board inside UL WH300 (neg_piro == 2)
+#	May 24, 2018: - continued working (coord_trans == 2)
+#	May 25, 2018: - continued working
+#	May 30, 2018: - BUG: trimming did not re-calculate elapsed time
+#	Jun  5, 2018: - BUG: on -c main magnetometer and accelerometer vectors were
+#						 modified twice
+#	Jun  7, 2018: - relaxed definition of -e
+#	Jun  9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does not matter in practice)
+
+#----------------------------------------------------------------------
+# gRef() library
+#----------------------------------------------------------------------
+
+sub pl_mag_calib_begin($$$)															# initialize mag_calib plot
+{
+	my($pfn,$plotsize,$axlim) = @_;
+	my($xmin,$xmax) = (-$axlim,$axlim);
+	my($ymin,$ymax) = (-$axlim,$axlim);
+	GMT_begin($pfn,"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+	GMT_psxy('-Sc0.05');
+}
+
+sub pl_mag_calib_plot($$$)															# plot data point
+{
+	my($valid,$magX,$magY) = @_;
+	if ($valid)	{ print(GMT "> -Wred -Gred\n$magX $magY\n"); }
+	else 		{ print(GMT "> -Wgreen -Ggreen\n$magX $magY\n"); }
+}
+
+sub pl_mag_calib_end($$)															# finish mag_calib plot
+{
+	my($axlim,$HF_mag,$sensor_info) = @_;
+	
+	GMT_psxy('-Sc0.1 -Gblue');														# calibration circle
+	for (my($a)=0; $a<2*$pi; $a+=0.075) {
+		printf(GMT "%f %f\n",$HF_mag*sin($a),$HF_mag*cos($a));
+	}
+	
+	if ($axlim <= 0.1) {															# axes labels
+		GMT_psbasemap('-Bg1a.04f.001:"X Magnetic Field [Gauss]":/g1a0.02f0.001:"Y Magnetic Field [Gauss]":WeSn');
+	} else {
+		GMT_psbasemap('-Bg1a.1f.01:"X Magnetic Field [Gauss]":/g1a0.1f0.01:"Y Magnetic Field [Gauss]":WeSn');
+	}
+    GMT_unitcoords();																# horizontal field strength
+    GMT_pstext('-F+f12,Helvetica,blue+jTR -N');
+   	printf(GMT "0.98 0.98 HF = %.2f Gauss\n",$HF_mag);
+   	
+   	printf(GMT "0.98 0.94 $sensor_info\n")											# sensor info
+		if ($sensor_info ne '');
+
+    GMT_pstext('-F+f14,Helvetica,blue+jTL -N');										# profile id
+    printf(GMT "0.01 1.06 $P{profile_id}\n");
+	GMT_end();
+}
+
+sub rot_vecs(@) 																	# rotate & output IMU vector data 
+{
+	my($coord_trans,$min_elapsed,$max_elapsed,$plot_milapsed,$plot_malapsed) = @_;		# negate KVH pitch/roll data if first arg set to 1
+	$min_elapsed = 0 unless defined($min_elapsed);
+	$max_elapsed = 9e99 unless defined($max_elapsed);
+	$plot_milapsed = $min_elapsed unless defined($plot_milapsed);
+	$plot_malapsed = $max_elapsed unless defined($plot_malapsed);
+
+	while (&antsIn()) {
+		next if ($ants_[0][$elapsedF] < $min_elapsed);								# trim data
+		last if ($ants_[0][$elapsedF] > $max_elapsed);
+		
+		my($cpiro) = -1;															# current pitch/roll accelerometer
+		my(@R); 																	# rotation matrix
+		for (my($i)=0; $i<@vecs; $i++) {											# rotate vector data
+			if ($piro[$i][0] != $cpiro) {											# next sensor chip
+				$cpiro = $piro[$i][0];
+
+				my($accX) = $ants_[0][$vecs[$cpiro][0]];
+				my($accY) = $ants_[0][$vecs[$cpiro][1]];
+				my($accZ) = $ants_[0][$vecs[$cpiro][2]];
+
+				if ($coord_trans == 2) {											# S/N 8 inside UL WH300
+					$accY *= -1; $accZ *= -1;
+				}
+				
+				my($roll)  = atan2($accY,$accZ); 									# eqn 25 from Freescale AN3461
+				my($pitch) = atan2($accX,sqrt($accY**2+$accZ**2));     				# eqn 26 from <Freescale AN3461
+				if ($coord_trans == 1) {											# KVH
+					$pitch *= -1;
+					$roll  *= -1;
+                }
+				$ants_[0][$piro[$i][1]] = deg($pitch);								# add pitch/roll to data
+				$ants_[0][$piro[$i][2]] = deg($roll);
+
+				my($sp) = sin($pitch); my($cp) = cos($pitch);						# define rotation matrix
+				my($sr) = sin($roll);  my($cr) = cos($roll);
+				@R = ([ $cp,	 0,   -$sp	  ],
+					  [-$sp*$sr, $cr, -$cp*$sr],
+					  [ $sp*$cr, $sr,  $cp*$cr]);
+			}
+			my($xval) = $ants_[0][$vecs[$i][0]];
+			my($yval) = $ants_[0][$vecs[$i][1]];
+			my($zval) = $ants_[0][$vecs[$i][2]];
+			
+			if ($coord_trans == 2) {												# S/N 8 inside UL WH300
+				$yval *= -1; $zval *= -1;
+			}
+			
+			$ants_[0][$vecs[$i][0]] = ($xval-$bias[$i][0]) * $R[0][0] +
+									  ($yval-$bias[$i][1]) * $R[0][1] +
+									  ($zval-$bias[$i][2]) * $R[0][2];
+			$ants_[0][$vecs[$i][1]] = ($xval-$bias[$i][0]) * $R[1][0] +
+			  						  ($yval-$bias[$i][1]) * $R[1][1] +
+									  ($zval-$bias[$i][2]) * $R[1][2];
+			$ants_[0][$vecs[$i][2]] = ($xval-$bias[$i][0]) * $R[2][0] +
+									  ($yval-$bias[$i][1]) * $R[2][1] +
+									  ($zval-$bias[$i][2]) * $R[2][2];
+		}
+	
+		my($magX) = $ants_[0][$magXF];
+		my($magY) = $ants_[0][$magYF];
+		my($magZ) = $ants_[0][$magZF];
+		my($accX) = $ants_[0][$accXF];
+		my($accY) = $ants_[0][$accYF];
+		my($accZ) = $ants_[0][$accZF];
+
+		my($HF)   = sqrt($magX**2+$magY**2);
+		my($valid)= ($HF >= $minfac*$HF_mag) && ($HF <= $maxfac*$HF_mag);
+		my($hdg)  = $valid ? mag_heading($magX,$magY) : nan;
+
+		&antsOut($ants_[0][$elapsedF]-$min_elapsed,$ants_[0][$tempF],
+				 RDI_pitch($ants_[0][$pitchF],$ants_[0][$rollF]),$ants_[0][$rollF],
+				 $hdg,$accX,$accY,$accZ,$magX,$magY,$magZ,
+				 vel_u($HF,$hdg),vel_v($HF,$hdg),$valid);
+
+		pl_mag_calib_plot($valid,$magX,$magY)
+			if defined($P{profile_id}) &&
+				($ants_[0][$elapsedF] >= $plot_milapsed) &&
+				($ants_[0][$elapsedF] <= $plot_malapsed);
+	}
+}
+
+#----------------------------------------------------------------------
+# LADCP merging library
+#----------------------------------------------------------------------
+
+#-----------------------------------------------------
+# Instrument Offset Estimation
+#
+#	1: resolution of histogram in deg
+#		1 deg okay for good sensors
+#		2 deg sometimes required (2016 P18 003 2nd sensor)
+#	2: min tilt anom to consider for offset estimation
+#		0.3 deg works even for calm casts
+#		increased values improve histogram
+#		2.0 deg is too high for quiet casts (2016 P18 003)
+#	3: minimum fraction of hist mode required
+#		10% default
+#		decreasing histogram resolution is better than
+#			decreasing this value, I think
+#-----------------------------------------------------
+
+#----------------------------------------------------------------------
+# trim_out_of_water()
+#	- attempt to remove out-of-water records from time-series of
+#	  horizontal acceleration
+#----------------------------------------------------------------------
+
+sub trim_out_of_water($)
+{
+	my($verbose) = @_;
+
+	#--------------------------------------------------------------------------
+	# first-difference horizontal acceleration at full resolution to pass-filter
+	# dangling motion
+	#--------------------------------------------------------------------------
+
+	$IMP{Ah}[0] 	= sqrt($ants_[0][$accXF]**2+$ants_[0][$accYF]**2);
+	$IMP{dAhdt}[0]	= nan;
+	for (my($r)=1; $r<@ants_; $r++) {
+		$IMP{Ah}[$r] = sqrt($ants_[$r][$accXF]**2+$ants_[$r][$accYF]**2);
+		$IMP{dAhdt}[$r] = ($IMP{Ah}[$r]-$IMP{Ah}[$r-1]) / ($ants_[$r][$elapsedF]-$ants_[$r-1][$elapsedF]);
+	}
+
+	#--------------------------------------------------------------------------------------
+	# create 10-s binned time series to calculate rms values of this quantity (dAhdt), and
+	# scale this with cos(sqrt($$pitch**2+$$roll**2)) to dampen underwater peaks (when the
+	# instrument has a large tilt because it is being dragged)
+	#	NB: dAhdt, pitch and roll are only set up to last full bin (not so, sum and n)
+	#--------------------------------------------------------------------------------------
+
+	my(@dAhdt,@pitch,@roll,@dAhdt_rms);
+	my(@sum) = my(@sume) = my(@sump) = my(@sumr) = my(@n) = (0);
+	
+	my($bin_start) = $ants_[0][$elapsedF];
+	for (my($r)=1; $r<@ants_; $r++) {   
+	
+		if ($ants_[$r][$elapsedF] - $bin_start <= 10) { 						# within 10-s bin
+			$sum[$#sum] += $IMP{dAhdt}[$r]; 									# sums
+			$sume[$#sume] += $ants_[$r][$elapsedF];
+			$sump[$#sump] += $ants_[$r][$pitchF];
+			$sumr[$#sumr] += $ants_[$r][$rollF];
+			$n[$#n]++;
+			next;
+		}
+	
+		$dAhdt[$#sum] = $sum[$#sum] / $n[$#n];								# bin done => means
+		$elapsed[$#sum] = $sume[$#sume] / $n[$#n];
+		$pitch[$#sum] = $sump[$#sump] / $n[$#n];
+		$roll[$#sum] = $sumr[$#sumr] / $n[$#n];
+	
+		my($sumsq) = 0; 													# sum of squares for rms(accel)
+		for (my($rr)=$r-$n[$#n]; $rr<$r; $rr++) {
+			$sumsq += ($IMP{dAhdt}[$rr] - $dAhdt[$#sum])**2;
+		}
+		$dAhdt_rms[$#sum] = sqrt($sumsq / $n[$#n]);
+	    
+		push(@sum,$IMP{dAhdt}[$r]); 										# begin next bin
+		push(@sume,$ants_[$r][$elapsedF]);
+		push(@sump,$ants_[$r][$pitchF]);
+		push(@sumr,$ants_[$r][$rollF]);
+		push(@n,1);
+		$bin_start = $ants_[$r][$elapsedF];
+	}
+
+	#--------------------------------------------
+	# trim beginning/end when IMP is out of water
+	#--------------------------------------------
+
+	my($i,$si);
+	for ($i=int(@dAhdt_rms/2); $i>0; $i--) {
+		last if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0);
+	}
+	if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0) {
+		for ($si=0; $ants_[$si][$elapsedF]<=$elapsed[$i]; $si++) {}
+		splice(@ants_,0,$si);
+		printf(STDERR "\n\t\t%5d  leading out-of-water IMU records removed",$si)
+			if ($si>0 && $verbose);
+	} else {
+		print(STDERR "\n\t\tWARNING: no leading out-of-water IMU records detected/removed") if $verbose;
+	}
+	
+	for ($i=int(@dAhdt_rms/2); $i<@dAhdt_rms; $i++) {
+		last if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0);
+	}
+	if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0) {
+		for ($si=$#ants_; $ants_[$si][$elapsedF]>=$elapsed[$i]; $si--) {}
+		my($rem) = @ants_ - $si;
+		splice(@ants_,$si);
+		printf(STDERR "\n\t\t%5d trailing out-of-water IMU records removed",$rem)
+			if ($rem>0 && $verbose);
+	} else {
+		print(STDERR "\n\t\tWARNING: no trailing out-of-water IMU records detected/removed") if $verbose;
+	}
+	
+	printf(STDERR "\n\t\tcast duration		  : %.1f min",
+		($ants_[$#ants_][$elapsedF] - $ants_[0][$elapsedF]) / 60)
+	        if $verbose;
+}
+
+#--------------------------------------------------------------------------------
+# prep_piro_IMP()
+# 	- calculate pitch/roll offsets & tilt azimuth for IMP
+#	- during an attempt to improve time lagging for the 2015 IOPAS data set,
+#	  it was noticed that one particular instrument, WHM300#12973 maxes out
+#	  pitch at 27.36 degrees, whereas the roll may not be maxed out at 28.76 deg,
+#	  the max observed during the cruise.
+#	- therefore, IMP{TILT_AZIM} and IMP{TILT_ANOM} are calculated here, first,
+#	  with a pitch/roll cutoff value of 29 degrees
+#	- after the time lagging, when the LADCP start and end times are known,
+#	  the TILT values are re-calculated without the pitch/roll limit, and
+#	  using only the correct time range
+#---------------------------------------------------------------------------------
+
+sub prep_piro_IMP($)
+{
+	my($verbose) = @_;
+	my($RDI_max_tilt) = 29; 
+	my($IMP_pitch_mean,$IMP_roll_mean,$nPR) = (0,0,0);
+	
+	for (my($r)=0; $r<@ants_; $r++) {
+		next unless numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF]);
+		$nPR++;
+		$IMP_pitch_mean += min($ants_[$r][$pitchF],$RDI_max_tilt);
+		$IMP_roll_mean	+= min($ants_[$r][$rollF],$RDI_max_tilt);
+	}
+	$IMP_pitch_mean /= $nPR;
+	$IMP_roll_mean /= $nPR;
+	printf(STDERR "\n\t\tIMU mean pitch/roll	  : %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
+			if $verbose;
+	
+	for (my($r)=0; $r<@ants_; $r++) {
+		next unless numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF]);
+		$IMP{TILT_AZIMUTH}[$r] = tilt_azimuth(min($ants_[$r][$pitchF],$RDI_max_tilt)-$IMP_pitch_mean,
+											  min($ants_[$r][$rollF],$RDI_max_tilt) -$IMP_roll_mean);
+		$IMP{TILT_ANOM}[$r] = angle_from_vertical(min($ants_[$r][$pitchF],$RDI_max_tilt)-$IMP_pitch_mean,
+												  min($ants_[$r][$rollF],$RDI_max_tilt) -$IMP_roll_mean);
+	}
+	return ($IMP_pitch_mean,$IMP_roll_mean);
+}
+
+#----------------------------------------------------------------------
+# prep_piro_LADCP()
+#	- calculate pitch/roll offsets & tilt azimuth of LADCP
+#----------------------------------------------------------------------
+
+sub prep_piro_LADCP($)
+{
+	my($verbose) = @_;
+
+	my($LADCP_pitch_mean,$LADCP_roll_mean) = (0,0);
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+			gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+		$LADCP_pitch_mean += $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH};
+		$LADCP_roll_mean  += $LADCP{ENSEMBLE}[$ens]->{ROLL};
+	}
+	$LADCP_pitch_mean /= ($LADCP_end-$LADCP_begin+1);
+	$LADCP_roll_mean  /= ($LADCP_end-$LADCP_begin+1);
+	printf(STDERR "\n\t\tLADCP mean pitch/roll	  : %.1f/%.1f deg",$LADCP_pitch_mean,$LADCP_roll_mean)
+			if $verbose;
+	
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		$LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH} =
+			tilt_azimuth($LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$LADCP_pitch_mean,
+						 $LADCP{ENSEMBLE}[$ens]->{ROLL}-$LADCP_roll_mean);
+		$LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} =
+			angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$LADCP_pitch_mean,
+								$LADCP{ENSEMBLE}[$ens]->{ROLL}-$LADCP_roll_mean);
+	}
+	
+	print(STDERR "\n") if $verbose;
+	return ($LADCP_pitch_mean,$LADCP_roll_mean);
+}
+
+#------------------------------------------------------------------
+# sub calc_hdg_offset()
+#	- estimate heading offset from tilt time series
+#	- returns heading offset and updated IMP mean tilts
+#	- also creates diagnostic plot with pl_hdg_offset()
+#------------------------------------------------------------------
+
+sub pl_hdg_offset($@)
+{
+	my($dhist_binsize,$modefrac,@dhist) = @_;
+	
+	my($plotsize) = '13c';
+	my($xmin,$xmax) = (-180.5,180.5);
+	my($ymin) = 0;
+	my($ymax) = 1.05 * $dhist[$HDG_offset];
+	    
+	GMT_begin("$P{profile_id}${opt_a}_hdg_offset.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+	GMT_psxy("-Sb${dhist_binsize}u -GCornFlowerBlue");
+	for (my($i)=0; $i<@dhist; $i++) {
+		next unless $dhist[$i];
+		printf(GMT "%f $dhist[$i]\n",$i*$dhist_binsize>180 ? $i*$dhist_binsize-360 : $i*$dhist_binsize);
+	}
+	GMT_psbasemap('-Bg45a90f15:"IMU Heading Offset [\260]":/ga100f10:"Frequency":WeSn');
+	GMT_unitcoords();
+	GMT_pstext('-F+f14,Helvetica,CornFlowerBlue+jTR -N');
+	printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+	GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+	printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
+    GMT_end();
+}
+
+sub calc_hdg_offset($)
+{
+	my($verbose) = @_;
+
+	print(STDERR "\n\tRe-calculating IMU pitch/roll anomalies") if $verbose;
+	($IMP_pitch_mean,$IMP_roll_mean,$nPR) = (0,0,0);
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+		if ($r < 0 && $ens == $LADCP_begin) {
+			$r = int(($LADCP{ENSEMBLE}[++$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT})
+				while ($r < 0);
+			printf(STDERR "\n\tIMU data begin with instrument already in water => skipping %ds of LADCP data",
+				$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME})
+					if ($verbose);
+			$LADCP_begin = $ens;
+		}
+		if ($r > $#ants_) {
+			printf(STDERR "\n\tIMU data end while instrument is still in water => truncating %ds of LADCP data",
+				$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME})
+					if ($verbose);
+			$LADCP_end = $ens - 1;
+			last;
+		}
+		next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+		$nPR++;
+		$IMP_pitch_mean += $ants_[$r][$pitchF];
+		$IMP_roll_mean	+= $ants_[$r][$rollF];
+	}
+	$IMP_pitch_mean /= $nPR;
+	$IMP_roll_mean /= $nPR;
+	
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+		next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+		$LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} =
+			$IMP{TILT_AZIMUTH}[$r] = tilt_azimuth($ants_[$r][$pitchF]-$IMP_pitch_mean,
+												  $ants_[$r][$rollF] -$IMP_roll_mean);
+		$LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} =
+			$IMP{TILT_ANOM}[$r] = angle_from_vertical($ants_[$r][$pitchF]-$IMP_pitch_mean,
+													  $ants_[$r][$rollF] -$IMP_roll_mean);
+	}
+	
+	printf(STDERR "\n\t\tIMU mean pitch/roll: %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
+			if $verbose;
+	
+	if (defined($opt_s)) {
+		$HDG_offset = $opt_s;
+		printf(STDERR "\n\tHEADING_OFFSET = %.1f (set by -s)\n",$HDG_offset) if $verbose;
+	} else {
+		my($dhist_binsize,$dhist_min_pirom,$dhist_min_mfrac) = split(/,/,$opt_e);
+		croak("$0: cannot decode -e $opt_e\n")
+			unless ($dhist_binsize > 0 && $dhist_min_pirom > 0 && $dhist_min_mfrac >= 0);
+	
+		my(@dhist); my($nhist) = my($modeFreq) = 0;
+		for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+			my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+			next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+			next unless (abs($ants_[$r][$pitchF]-$IMP_pitch_mean) >= $dhist_min_pirom &&
+						 abs($ants_[$r][$rollF] -$IMP_roll_mean) >= $dhist_min_pirom);
+			$dhist[int(angle_pos($LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH}-$IMP{TILT_AZIMUTH}[$r])/$dhist_binsize+0.5)]++;
+			$nhist++;
+		}
+		croak("$0: empty histogram\n")
+			unless ($nhist);
+	
+		$HDG_offset = 0;
+		for (my($i)=1; $i<@dhist-1; $i++) { 							# make sure mode is not on edge
+			$HDG_offset = $i if ($dhist[$i] >= $dhist[$HDG_offset]);
+		}
+		my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
+	    
+		$HDG_offset *= $dhist_binsize;
+		pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
+	    
+		if ($opt_f) {
+			printf(STDERR "\n\nIGNORED WARNING (-f): Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100)
+				if ($modefrac < $dhist_min_mfrac);
+		} else {
+			croak(sprintf("\n$0: Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100))
+				if ($modefrac < $dhist_min_mfrac);
+		}
+	
+		printf(STDERR "\n\t") if $verbose;
+		printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac) if $verbose;
+	}
+	return ($HDG_offset,$IMP_pitch_mean,$IMP_roll_mean);
+}
+
+#-----------------------------------------------------------
+# rot_IMP()
+# 	- rotate IMP Data Into LADCP Instrument Coords
+#	- also replaced pitch/roll by corresponding anomalies!!!
+#-----------------------------------------------------------
+
+sub rot_IMP($)
+{
+	my($verbose) = @_;
+	my($crho) = cos(rad($HDG_offset));
+	my($srho) = sin(rad($HDG_offset));
+	
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+	
+		if (numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF])) { 				# pitch/roll
+			my($rot_p) = (($ants_[$r][$pitchF]-$IMP_pitch_mean)*$crho +
+						  ($ants_[$r][$rollF]-$IMP_roll_mean)*$srho);
+			my($rot_r) = (-($ants_[$r][$pitchF]-$IMP_pitch_mean)*$srho +
+						   ($ants_[$r][$rollF]-$IMP_roll_mean)*$crho);
+			$ants_[$r][$pitchF] = $rot_p;
+			$ants_[$r][$rollF]	= $rot_r;
+		}
+	    
+		$ants_[$r][$hdgF] = angle_pos($ants_[$r][$hdgF] - $HDG_offset)
+			if numberp($ants_[$r][$hdgF]);
+	}
+	
+	my($rot_p) =  $IMP_pitch_mean * $crho + $IMP_roll_mean * $srho; 			# mean pitch roll
+	my($rot_r) = -$IMP_pitch_mean * $srho + $IMP_roll_mean * $crho;
+	$IMP_pitch_mean = $rot_p;
+	$IMP_roll_mean	= $rot_r;
+	
+	print(STDERR "\n") if $verbose;
+	return ($IMP_pitch_mean,$IMP_roll_mean);
+}
+
+#----------------------------------------------------------------------
+# create_merge_plots()
+#   - tilt time series (*_time_lag.ps)
+#----------------------------------------------------------------------
+
+sub create_merge_plots($$$)
+{
+    my($basename,$plotsize,$verbose) = @_;
+
+	#---------------------------------
+	# Tilt Time Series (*_time_lag.ps)
+	#---------------------------------
+
+	my($mint,$maxt) = (99,-99);
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		if (($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 180) ||
+			($LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 180)) {
+				$mint = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} < $mint);
+				$maxt = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} > $maxt);
+				$mint = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < $mint);
+				$maxt = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} > $maxt);
+		}
+	}
+
+	my($xmin,$xmax) = (-90,90);
+	my($ymin) = round($mint-0.5);
+	my($ymax) = round($maxt+0.5);
+	
+	GMT_begin("${basename}_time_lag.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+
+	GMT_psxy('-W2,coral');
+	for (my($ens) = $LADCP_begin + 5; 
+		 $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 90;
+		 $ens++) {
+			printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME},
+								 $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM});
+	}
+	GMT_psxy('-W1');
+	for (my($ens) = $LADCP_begin + 5; 
+		 $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 90;
+		 $ens++) {
+			printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME},
+								 $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM});
+	}
+	GMT_psxy('-W2,SeaGreen');
+	for (my($ens) = $LADCP_end - 5; 
+		 $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 90;
+		 $ens--) {
+			printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
+								 $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM});
+	}
+	GMT_psxy('-W1');
+	for (my($ens) = $LADCP_end - 5; 
+		 $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 90;
+		 $ens--) {
+			printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
+								 $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM});
+	}
+
+    GMT_psbasemap('-Bg30a30f10:"Elapsed Time [sec]":/g5a5f1:"Tilt Magnitude [\260]":WeSn');
+	GMT_unitcoords();
+    GMT_pstext('-F+f14,Helvetica,Coral+jTL');
+	    printf(GMT "0.52 0.98 downcast\n");
+    GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR');
+	    printf(GMT "0.48 0.98 upcast\n");	    
+	GMT_psxy('-W4,LightSkyBlue');
+		printf(GMT "0.5 0\n0.5 1\n");
+    GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+	    printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
+	GMT_end();
+
+	#------------------------------
+	# Heading Errors (*_hdg_err.ps)
+	#------------------------------
+
+	my(@err_binned,@err_nsamp);
+	my($sumErr) = 0; my($nErr) = $LADCP_end - $LADCP_begin + 1;
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		next unless ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < 10);
+		my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+		my($bi) = int($ants_[$r][$hdgF]/5);
+		my($err) = angle_diff($ants_[$r][$hdgF],$LADCP{ENSEMBLE}[$ens]->{HEADING});
+		next unless numberp($err);
+		$err_binned[$bi] += $err; $sumErr += $err;
+		$err_nsamp[$bi]++;
+	}
+	for (my($bi)=0; $bi<@err_nsamp; $bi++) {
+		$err_binned[$bi] = ($err_nsamp[$bi] >= 5)
+						 ? $err_binned[$bi]/$err_nsamp[$bi]
+						 : undef;
+	}
+	my(@err_dssq);
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		next unless ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < 10);
+		my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+		my($bi) = int($ants_[$r][$hdgF]/5);
+		my($err) = angle_diff($ants_[$r][$hdgF],$LADCP{ENSEMBLE}[$ens]->{HEADING});
+		next unless numberp($err);
+		$err_dssq[$bi] += ($err-$err_binned[$bi])**2;
+	}
+
+	my($xmin,$xmax) = (0,360);
+	my($ymin,$ymax) = (-45,45);
+	GMT_begin("${basename}_hdg_err.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+	GMT_psxy('-Ey/3,CornFlowerBlue');
+	my($sumSq,$sumBe) = my($nSq,$nBe) = (0,0);
+	for (my($bi)=0; $bi<@err_binned; $bi++) {
+		next unless ($err_nsamp[$bi] >= 2);
+		next unless numberp($err_binned[$bi]);
+		$sumSq += $err_binned[$bi]**2; $nSq++;
+		$sumBe += $err_binned[$bi]; $nBe++;
+#		printf(GMT "%f %f\n",2.5+5*$bi,$err_binned[$bi]);
+		printf(GMT "%f %f %f\n",2.5+5*$bi,$err_binned[$bi],sqrt($err_dssq[$bi]/($err_nsamp[$bi]-1)));
+	}
+    GMT_psbasemap('-Bg90a45f5:"ADCP Heading [\260]":/g15a15f5:"ADCP Compass Error [\260]":WeSn');
+	GMT_unitcoords();
+    GMT_pstext('-F+f12,Helvetica,CornFlowerBlue+jTR -Gwhite -C25%');
+    printf(GMT "0.98 0.98 rms error = %7.1f \260\n",sqrt($sumSq/$nSq));
+    printf(GMT "0.98 0.94 time-averaged error = %7.1f \260\n",$sumErr/$nErr);
+    printf(GMT "0.98 0.90 heading-averaged error = %7.1f \260\n",$sumBe/$nBe);
+    GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+    printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
+    GMT_end();
+	                        
+	print(STDERR "\n") if $verbose;
+}
+
+#----------------------------------------------------------------------
+# output_merged()
+#	- output merged data
+#----------------------------------------------------------------------
+
+sub output_merged($)
+{
+	my($verbose) = @_;
+	
+	my($tazimF)		= &antsNewField('tilt_azimuth');
+	my($tanomF)		= &antsNewField('tilt_magnitude');
+	my($L_tazimF) 	= &antsNewField('LADCP_tilt_azimuth');
+	my($L_tanomF) 	= &antsNewField('LADCP_tilt_magnitude');
+	my($L_elapsedF) = &antsNewField('LADCP_elapsed');
+	my($L_ensF) 	= &antsNewField('LADCP_ens');
+	my($L_depthF) 	= &antsNewField('LADCP_depth_estimate');
+	my($L_pitchF)	= &antsNewField('LADCP_pitch');
+	my($L_rollF)	= &antsNewField('LADCP_roll');
+	my($L_hdgF)		= &antsNewField('LADCP_hdg');
+	my($dcF)		= &antsNewField('downcast');
+
+	for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+		my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+#		print(STDERR "$r\[$ens\,$LADCP{ENSEMBLE}[$ens]->{NUMBER}] = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});\n");
+		if ($r<0 || $r>$#ants_) {												# ensemble beyond limits of IMU data
+			my(@out);
+			if ($ens >= $LADCP_begin && $ens <= $LADCP_end) {
+				$out[$elapsedF] 	= $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG};
+				$out[$L_elapsedF] 	= $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME};
+			}
+			$out[$L_ensF]			= $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+			$out[$dcF]				= ($ens <= $LADCP_bottom);
+			&antsOut(@out);
+		} elsif ($ens < $LADCP_begin || $ens > $LADCP_end) {					# pre deplyment or post recovery
+			$ants_[$r][$L_elapsedF] = undef;
+			$ants_[$r][$L_ensF] 	= $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+			$ants_[$r][$L_pitchF]	= $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
+			$ants_[$r][$L_rollF]	= $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
+			$ants_[$r][$L_hdgF] 	= $LADCP{ENSEMBLE}[$ens]->{HEADING};							    
+			$ants_[$r][$dcF]        = nan;
+			&antsOut(@{$ants_[$r]});
+		} else {
+			$ants_[$r][$tazimF] 	= angle($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} + $HDG_offset);
+			$ants_[$r][$tanomF] 	= $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM};
+			$ants_[$r][$L_tazimF]	= $LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH};
+			$ants_[$r][$L_tanomF]	= $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM};
+			$ants_[$r][$L_elapsedF] = $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME};
+			$ants_[$r][$L_ensF] 	= $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+			$ants_[$r][$L_depthF]	= $LADCP{ENSEMBLE}[$ens]->{DEPTH};
+			$ants_[$r][$L_pitchF]	= $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
+			$ants_[$r][$L_rollF]	= $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
+			$ants_[$r][$L_hdgF] 	= $LADCP{ENSEMBLE}[$ens]->{HEADING};							    
+			$ants_[$r][$dcF]        = ($ens <= $LADCP_bottom);
+			&antsOut(@{$ants_[$r]});
+		}
+	}
+	
+	print(STDERR "\n") if $verbose;
+}
+
+#----------------------------------------------------------------------
+
+1; # return true for all the world to see
diff --git a/antsio.pl b/antsio.pl
--- a/antsio.pl
+++ b/antsio.pl
@@ -2,9 +2,9 @@
 #======================================================================
 #                    A N T S I O . P L 
 #                    doc: Fri Jun 19 19:22:51 1998
-#                    dlm: Mon Apr 23 14:20:58 2018
+#                    dlm: Wed Apr 10 16:57:59 2019
 #                    (c) 1998 A.M. Thurnherr
-#                    uE-Info: 216 82 NIL 0 0 70 2 2 4 NIL ofnI
+#                    uE-Info: 217 48 NIL 0 0 70 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -214,6 +214,7 @@
 #	Mar 10, 2017: - BUG: antsCheckDeps() used ctime instead of mtime!!!
 #	Apr  5, 2017: - BUG: stale file mtime dependency info was not printed correctly
 #	Apr 23, 2018: - BUG: @antsLayout was not kept up-to-date when layout-changes are allowed
+#	Apr 10, 2019: - disabled dependency warnings
 
 # GENERAL NOTES:
 #	- %P was named without an ants-prefix because associative arrays
@@ -292,7 +293,7 @@
 #	- by default, tests current <> file
 #----------------------------------------------------------------------
 
-{ my($warned);
+{ my($warned) = 1;									# disable dependency warning
 
   sub antsCheckDeps()
   {
diff --git a/antsutils.pl b/antsutils.pl
--- a/antsutils.pl
+++ b/antsutils.pl
@@ -2,9 +2,9 @@
 #======================================================================
 #                    A N T S U T I L S . P L 
 #                    doc: Fri Jun 19 23:25:50 1998
-#                    dlm: Wed May 25 12:16:39 2016
+#                    dlm: Fri Aug 30 13:39:50 2019
 #                    (c) 1998 A.M. Thurnherr
-#                    uE-Info: 104 31 NIL 0 0 70 10 2 4 NIL ofnI
+#                    uE-Info: 108 62 NIL 0 0 70 10 2 4 NIL ofnI
 #======================================================================
 
 # Miscellaneous auxillary functions
@@ -102,6 +102,10 @@
 #	May 18, 2015: - added antsFindParam()
 #	Jun 21, 2015: - added antsParam(), modified antsRequireParam()
 #	May 12, 2016: - added &div2() to prevent division by zero errors
+#	Apr  5, 2019: - disabled weird line of code in antsFunUsage() (see comment)
+#				  - improved error messages in antsFunUsage()
+#				  - BUG: antsFunUsage did not work with -ve argc (variable argument funs)
+#	Aug 30, 2019: - BUG: antsLoadModel() did not respect $ANTS
 
 # fnr notes:
 #	- matches field names starting with the string given, i.e. "sig" is
@@ -446,8 +450,7 @@
 		require "$pref.$name";
 		return $name;
 	} else {											# load from ANTSlib 
-		my($path) = ($0 =~ m{^(.*)/[^/]*$});
-		require "$path/$pref.$name";
+		require "$ANTS/$pref.$name";
 		return $name;
     }
 }
@@ -467,8 +470,7 @@
 			require "$pref.$name";
 			return ($name,split(',',$args));
 		} else {
-			my($path) = ($0 =~ m{^(.*)/[^/]*$});
-			require "$path/$pref.$name";
+			require "$ANTS/$pref.$name";
 			return ($name,split(',',$args));
 		}
 	}
@@ -538,7 +540,7 @@
 
 	if (ref($params[0]) && @antsLayout>0 && @params<2*$argc+1) {  # default params
 		my(@newparams);									# 2nd test is for abc
-		my($npi) = $argc+1;
+		my($npi) = abs($argc)+1;
 
 		$listAllRecs = 1;								# special flag for list(1)
 
@@ -551,7 +553,7 @@
 			return(@newparams);
 		}
 	    
-		for (my($i)=1; $i<=$argc; $i++) {				# fill cache & do tests
+		for (my($i)=1; $i<=abs($argc); $i++) {				# fill cache & do tests
 			if (defined($params[$i])) {
 				push(@{$params[0]},&fnr($params[$i]));
 				push(@newparams,&antsVal($params[0]->[$#{$params[0]}]));
@@ -564,24 +566,27 @@
 		croak("usage: $msg\n") unless ($npi > $#params);
 		
 		@params = @newparams;
-	} elsif (ref($params[0])) {
-		splice(@params,0,$argc+1);
+	} elsif (ref($params[0])) {								# remove array ref & list of field names
+		splice(@params,0,abs($argc)+1);
 	}
 
 	if ($argc >= 0) {									# argument count
-		croak("usage: $msg\n") unless (@params == $argc);
+		croak("usage: $msg [params = @params]\n") unless (@params == $argc);
 	} else {
-		croak("usage: $msg\n") unless (@params >= -$argc);
+		croak("usage: $msg [params = @params])\n") unless (@params >= -$argc);
 	}
     
 	for (my($i)=0; $i<length($types); $i++) {			# type checking
 		$_ = substr($types,$i,1);
 		SWITCH: {
-			last unless defined($params[$i]);
-			&antsNoCardErr("",$params[$i]),last SWITCH if (/c/);
-			&antsNoIntErr("",$params[$i]),last SWITCH if (/i/);
-			&antsNoFloatErr("",$params[$i]),last SWITCH if (/f/);
-			&antsNoFileErr("",$params[$i]),last SWITCH if (/F/);
+# 4/5/19: The following line of code prevents proper type checking when one of the
+#		  arguments is undefined. I do not know under what circumstances the code
+#	      is required. Therfore I disabled it temporarily.
+#			last unless defined($params[$i]);
+			&antsNoCardErr(sprintf("argument #%d in $msg (params = @params)",$i+1),$params[$i]),last SWITCH if (/c/);
+			&antsNoIntErr(sprintf("argument #%d in $msg",$i+1),$params[$i]),last SWITCH if (/i/);
+			&antsNoFloatErr(sprintf("argument #%d in $msg (params = @params)",$i+1),$params[$i]),last SWITCH if (/f/);
+			&antsNoFileErr(sprintf("argument #%d in $msg",$i+1),$params[$i]),last SWITCH if (/F/);
 			if (/\d/) {
 				croak("$0: $params[$i] is not a string of length $_\n")
 					unless ($_ == length($params[$i]));
diff --git a/libGM76.pl b/libGM76.pl
--- a/libGM76.pl
+++ b/libGM76.pl
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L I B G M 7 6 . P L 
 #                    doc: Sun Feb 20 14:43:47 2011
-#                    dlm: Sat Feb  2 02:33:11 2019
+#                    dlm: Sat Apr  6 19:33:05 2019
 #                    (c) 2011 A.M. Thurnherr
-#                    uE-Info: 144 43 NIL 0 0 70 2 2 4 NIL ofnI
+#                    uE-Info: 34 48 NIL 0 0 70 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -27,7 +27,11 @@
 #				  - added GM_strain
 #				  - BUG: Sw usage message had wrong parameter order
 #				  - renamed Sw => GM_VKE; Su_z => GM_shear
-
+#	Mar 31, 2019: - updated doc for shear and strain spectra
+#	Apr  1, 2019: - added GM_vdiv
+#	Apr  5, 2019: - BUG: GM_VKE was erroneous (argument shifting was wrong)
+#				  - adapted to improved antsFunUsage()
+#	Apr  6, 2019: - fiddled during debugging of fs_finestructure issue
 
 require "$ANTS/libEOS83.pl";
 
@@ -100,7 +104,7 @@
 sub GM_VKE(@)
 {
 	my($omega,$m,$lat,$b,$jstar,$N) =
-		&antsFunUsage(-5,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+		&antsFunUsage(-5,'fffff','GM_VKE([frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
 
 	if (defined($N)) {									# Sw(omega,m)
 		local($f) = abs(&f($lat));
@@ -110,8 +114,8 @@
 		my($mstar) = &m($jstar,$N,$omega);
 		return $GM_E0 * $b * 2 * $f**2/$omega**2/B($omega) * $jstar / ($m+$mstar)**2;
 	} else {											# Sw(m), i.e. integrated over all omega; as in Thurnherr et al., GRL 2015
-		$N = $lat;										# shift arguments to account for missing omega
-		$lat = $m;
+		$N = $jstar;  $jstar = $b;						# shift arguments to account for missing omega
+		$b = $lat; $lat = $m;
 		local($f) = abs(&f($lat));
 		$m = $omega;
 		undef($omega);
@@ -121,15 +125,38 @@
 }
 
 #----------------------------------------------------------------------
+# Vertical Divergence (normalized by N)
+#	- implemented to investigate shear-to-vertical divergence ratio
+#	- S[w_z/N] = S[w_z] / N^2 = S[w] * m^2 / N^2
+#	- GM shear-to-vd variance ratio = 3N/2f
+#	- GM strain-to-vd variance ratio = N/2f
+#----------------------------------------------------------------------
+
+sub GM_vdiv(@)
+{
+	my($m,$lat,$b,$jstar,$N) =
+		&antsFunUsage(5,'fffff','GM_vdiv(<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
+	return GM_VKE($m,$lat,$b,$jstar,$N) * $m**2 / $N**2;
+}
+
+#----------------------------------------------------------------------
 # Shear and Strain m-spectra (i.e. integrated over f)
+#	- shear is buoyancy-frequency normalized
 #	- spectral density
-#	- from Thurnherr et al. (GRL 2015)
+#	- from Thurnherr et al. (GRL 2015), which is from email info
+#	  Eric Kunze, 04/22/2015: First, the standard GHP parameterization that I
+#	  implement uses either the N-normalized shear or strain spectra. The GM
+#	  versions of these are
+#		S[V_z/N](m) = (3*PI/2)*E0*b*jstar*m^2/(m+mstar)^2 and
+#		S[Z_z](m) = (PI/2)*E0*b*jstar*m^2/(m+mstar)^2
+#	  though in practice mstar is not particularly important because the
+#	  variance of both these quantities is dominated by high m.
 #----------------------------------------------------------------------
 
 sub GM_shear(@)
 {
 	my($m,$lat,$b,$jstar,$N) =
-		&antsFunUsage(5,'fffff','<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+		&antsFunUsage(5,'fffff','GM_shear(<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
 	local($f) = abs(&f($lat));
 	my($mstar) = &m($jstar,$N);
 	return 3 * $pi/2 * $GM_E0 * $b * $jstar * $m**2 / ($m+$mstar)**2;
@@ -138,7 +165,7 @@
 sub GM_strain(@)
 {
 	my($m,$lat,$b,$jstar,$N) =
-		&antsFunUsage(5,'fffff','<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+		&antsFunUsage(5,'fffff','GM_strain(<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
 	local($f) = abs(&f($lat));
 	my($mstar) = &m($jstar,$N);
 	return $pi/2 * $GM_E0 * $b * $jstar * $m**2 / ($m+$mstar)**2;
diff --git a/libIMP.pl b/libIMP.pl
--- a/libIMP.pl
+++ b/libIMP.pl
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L I B I M P . P L 
 #                    doc: Tue Nov 26 21:59:40 2013
-#                    dlm: Sat Jun  9 12:07:15 2018
+#                    dlm: Mon Jul  1 15:47:57 2019
 #                    (c) 2017 A.M. Thurnherr
-#                    uE-Info: 52 109 NIL 0 0 70 2 2 4 NIL ofnI
+#                    uE-Info: 689 0 NIL 0 0 70 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -49,7 +49,16 @@
 #	Jun  5, 2018: - BUG: on -c main magnetometer and accelerometer vectors were
 #						 modified twice
 #	Jun  7, 2018: - relaxed definition of -e
-#	Jun  9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does not matter in practice)
+#	Jun  9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does
+#						 not matter in practice)
+#	Jun 30, 2019: - renamed -s in IMP+LADCP to -o (option not used in KVH+LADCP)
+#				  - changed semantics so that offset histogram is plotted even on -o
+#				  - made create_merge_plots return errors
+#	Jul  1, 2019: - BUG: pre-deployment output cleared IMU pitch/roll fields (1 record
+#						 affected)
+#				  - modified output of pre-/post-cast records to include all LADCP
+#				    but no IMU data
+#				  - BUG: GIMBAL_PITCH had not only been defined for in-water records
 
 #----------------------------------------------------------------------
 # gRef() library
@@ -349,6 +358,10 @@
 	my($verbose) = @_;
 
 	my($LADCP_pitch_mean,$LADCP_roll_mean) = (0,0);
+	for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+		$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+			gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+	}
 	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
 		$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
 			gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
@@ -398,7 +411,11 @@
 	GMT_psbasemap('-Bg45a90f15:"IMU Heading Offset [\260]":/ga100f10:"Frequency":WeSn');
 	GMT_unitcoords();
 	GMT_pstext('-F+f14,Helvetica,CornFlowerBlue+jTR -N');
-	printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+	if (defined($opt_o)) {
+		printf(GMT "0.99 1.06 -o %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+	} else {
+		printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+	}
 	GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
 	printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
     GMT_end();
@@ -449,35 +466,36 @@
 	printf(STDERR "\n\t\tIMU mean pitch/roll: %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
 			if $verbose;
 	
-	if (defined($opt_s)) {
-		$HDG_offset = $opt_s;
-		printf(STDERR "\n\tHEADING_OFFSET = %.1f (set by -s)\n",$HDG_offset) if $verbose;
-	} else {
-		my($dhist_binsize,$dhist_min_pirom,$dhist_min_mfrac) = split(/,/,$opt_e);
-		croak("$0: cannot decode -e $opt_e\n")
-			unless ($dhist_binsize > 0 && $dhist_min_pirom > 0 && $dhist_min_mfrac >= 0);
+	my($dhist_binsize,$dhist_min_pirom,$dhist_min_mfrac) = split(/,/,$opt_e);
+	croak("$0: cannot decode -e $opt_e\n")
+		unless ($dhist_binsize > 0 && $dhist_min_pirom > 0 && $dhist_min_mfrac >= 0);
 	
-		my(@dhist); my($nhist) = my($modeFreq) = 0;
-		for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
-			my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
-			next unless numberp($IMP{TILT_AZIMUTH}[$r]);
-			next unless (abs($ants_[$r][$pitchF]-$IMP_pitch_mean) >= $dhist_min_pirom &&
-						 abs($ants_[$r][$rollF] -$IMP_roll_mean) >= $dhist_min_pirom);
-			$dhist[int(angle_pos($LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH}-$IMP{TILT_AZIMUTH}[$r])/$dhist_binsize+0.5)]++;
-			$nhist++;
-		}
-		croak("$0: empty histogram\n")
-			unless ($nhist);
+	my(@dhist); my($nhist) = my($modeFreq) = 0;
+	for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+		my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+		next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+		next unless (abs($ants_[$r][$pitchF]-$IMP_pitch_mean) >= $dhist_min_pirom &&
+					 abs($ants_[$r][$rollF] -$IMP_roll_mean) >= $dhist_min_pirom);
+		$dhist[int(angle_pos($LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH}-$IMP{TILT_AZIMUTH}[$r])/$dhist_binsize+0.5)]++;
+		$nhist++;
+	}
+	croak("$0: empty histogram\n")
+		unless ($nhist);
 	
+	if (defined($opt_o)) {												# set heading offset, either with -o
+		$HDG_offset = $opt_o;
+	} else {															# or from the data
 		$HDG_offset = 0;
 		for (my($i)=1; $i<@dhist-1; $i++) { 							# make sure mode is not on edge
 			$HDG_offset = $i if ($dhist[$i] >= $dhist[$HDG_offset]);
 		}
-		my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
-	    
 		$HDG_offset *= $dhist_binsize;
-		pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
-	    
+	}
+
+	my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
+	pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
+
+	unless (defined($opt_o)) {	    									# make sure data are consistent (unless -o is used)
 		if ($opt_f) {
 			printf(STDERR "\n\nIGNORED WARNING (-f): Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100)
 				if ($modefrac < $dhist_min_mfrac);
@@ -485,10 +503,17 @@
 			croak(sprintf("\n$0: Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100))
 				if ($modefrac < $dhist_min_mfrac);
 		}
+	}
 	
-		printf(STDERR "\n\t") if $verbose;
-		printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac) if $verbose;
+	printf(STDERR "\n\t") if $verbose;
+	if ($opt_o) {
+		printf(STDERR "IMU heading offset = -o %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac)
+			if $verbose;
+	} else {
+		printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac)
+			if $verbose;
 	}
+
 	return ($HDG_offset,$IMP_pitch_mean,$IMP_roll_mean);
 }
 
@@ -532,6 +557,7 @@
 #----------------------------------------------------------------------
 # create_merge_plots()
 #   - tilt time series (*_time_lag.ps)
+#	- heading errors (*_hdg_err.ps)
 #----------------------------------------------------------------------
 
 sub create_merge_plots($$$)
@@ -654,11 +680,13 @@
     GMT_end();
 	                        
 	print(STDERR "\n") if $verbose;
+	return (sqrt($sumSq/$nSq),$sumErr/$nErr,$sumBe/$nBe);
 }
 
 #----------------------------------------------------------------------
 # output_merged()
 #	- output merged data
+#	- there must be exactly one record for each PD0 ensemble
 #----------------------------------------------------------------------
 
 sub output_merged($)
@@ -690,16 +718,13 @@
 			$out[$dcF]				= ($ens <= $LADCP_bottom);
 			&antsOut(@out);
 		} elsif ($ens < $LADCP_begin || $ens > $LADCP_end) {					# pre deplyment or post recovery
-			$ants_[$r][$L_elapsedF] = undef;
-			$ants_[$r][$L_ensF] 	= $LADCP{ENSEMBLE}[$ens]->{NUMBER};
-			$ants_[$r][$L_pitchF]	= undef;
-			$ants_[$r][$L_rollF]	= undef;
-			$ants_[$r][$L_hdgF]		= undef;
-			$ants_[$r][$pitchF]		= undef;			
-			$ants_[$r][$rollF]		= undef;			
-			$ants_[$r][$hdgF]		= undef;			
-			$ants_[$r][$dcF]        = ($ens <= $LADCP_bottom);
-			&antsOut(@{$ants_[$r]});
+			my(@out);															# correct IMP record NOT known
+			$out[$L_elapsedF] 		= undef;
+			$out[$L_ensF] 			= $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+			$out[$L_pitchF]			= $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
+			$out[$L_rollF]			= $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
+			$out[$L_hdgF]			= $LADCP{ENSEMBLE}[$ens]->{HEADING};							    
+			&antsOut(@out);
 		} else {
 			$ants_[$r][$tazimF] 	= angle($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} + $HDG_offset);
 			$ants_[$r][$tanomF] 	= $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM};
@@ -710,8 +735,8 @@
 			$ants_[$r][$L_depthF]	= $LADCP{ENSEMBLE}[$ens]->{DEPTH};
 			$ants_[$r][$L_pitchF]	= $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
 			$ants_[$r][$L_rollF]	= $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
-			$ants_[$r][$L_hdgF] 	= $LADCP{ENSEMBLE}[$ens]->{HEADING};
-			$ants_[$r][$dcF]        = ($ens <= $LADCP_bottom);
+			$ants_[$r][$L_hdgF] 	= $LADCP{ENSEMBLE}[$ens]->{HEADING};							    
+			$ants_[$r][$dcF]        = ($ens <= $LADCP_bottom) ? 1 : 0;
 			&antsOut(@{$ants_[$r]});
 		}
 	}
diff --git a/libLADCP.pl b/libLADCP.pl
--- a/libLADCP.pl
+++ b/libLADCP.pl
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L I B L A D C P . P L 
 #                    doc: Wed Jun  1 20:38:19 2011
-#                    dlm: Wed Apr 25 17:41:36 2018
+#                    dlm: Sat Apr  6 19:42:07 2019
 #                    (c) 2011 A.M. Thurnherr
-#                    uE-Info: 46 27 NIL 0 0 70 2 2 4 NIL ofnI
+#                    uE-Info: 211 34 NIL 0 0 70 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -24,6 +24,8 @@
 #				  - added parameter checks to processing-specific corrections
 #	May 18, 2015: - added pulse length to T_w() and T_w_z()
 #	Apr 25, 2018: - added eps_VKE() parameterization
+#	Apr  5, 2018: - adapted to improved antsFunUsage()
+#				  - BUG eps_VKE() had erroneous string
 
 require "$ANTS/libvec.pl";
 require "$ANTS/libfuns.pl";
@@ -43,7 +45,7 @@
 sub eps_VKE(@)
 {
     my($p0,$c) =
-        &antsFunUsage(-1,'.','<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]]',@_);
+        &antsFunUsage(-1,'.','eps_VKE(<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]])',@_);
 	$c = 0.021 unless defined($c);
     return numberp($p0) ? ($p0 / $c) ** 2 : nan;
 }
@@ -69,7 +71,7 @@
 sub T_ravg(@)
 {
     my($k,$blen,$plen) =
-        &antsFunUsage(-2,'ff','<vertical wavenumber[rad/s]> <bin-length[m]> [pulse-length[m]]',@_);
+        &antsFunUsage(-2,'ff','T_ravg(<vertical wavenumber[rad/s]> <bin-length[m]> [pulse-length[m]])',@_);
 	$plen = $blen unless defined($plen);        
     return 1 / sinc($k*$blen/2/$PI)**2 / sinc($k*$plen/2/$PI)**2;
 }
@@ -82,7 +84,7 @@
 sub T_fdiff($$)
 {
     my($k,$dz) =
-        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <differencing interval[m]>',@_);
+        &antsFunUsage(2,'ff','T_fdiff(<vertical wavenumber[rad/s]> <differencing interval[m]>)',@_);
     return 1 / sinc($k*$dz/2/$PI)**2;
 }
 
@@ -95,7 +97,7 @@
 sub T_interp($$$)
 {
     my($k,$blen,$dz) =
-        &antsFunUsage(3,'fff','<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>',@_);
+        &antsFunUsage(3,'fff','T_interp(<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>)',@_);
     return 1 / sinc($k*$blen/2/$PI)**4 / sinc($k*$dz/2/$PI)**2;
 }
 
@@ -108,7 +110,7 @@
 sub T_binavg($$)
 {
     my($k,$dz) =
-        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <grid resolution[m]>',@_);
+        &antsFunUsage(2,'ff','T_binavg(<vertical wavenumber[rad/s]> <grid resolution[m]>)',@_);
     return 1 / sinc($k*$dz/2/$PI)**2;
 }
 
@@ -125,7 +127,7 @@
 sub T_tilt($$)
 {
     my($k,$dprime) =
-        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <d-prime[m]>',@_);
+        &antsFunUsage(2,'ff','T_tilt(<vertical wavenumber[rad/s]> <d-prime[m]>)',@_);
     return $dprime ? 1 / sinc($k*$dprime/2/$PI)**2 : 1;
 }
 
@@ -147,7 +149,7 @@
 sub T_UH($$$$)
 {
 	my($k,$blen,$dz,$range_max) =
-        &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+        &antsFunUsage(-3,'fff','T_UH(<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>)',@_);
 	croak("T_UH($k,$blen,$dz,$range_max): bad parameters\n")
 		unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				
     return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_interp($k,$blen,$dz) * T_tilt($k,dprime($range_max));
@@ -162,7 +164,7 @@
 sub T_ASM($$$$)
 {
 	my($k,$blen,$dz,$range_max) =
-        &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+        &antsFunUsage(-3,'fff','T_ASM(<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>)',@_);
 	croak("T_ASM($k,$blen,$dz,$range_max): bad parameters\n")
 		unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				
     return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
@@ -180,14 +182,14 @@
 sub T_VI($$$$$)
 {
 	my($k,$blen,$sel,$dz,$range_max) =
-        &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <range max[m]>',@_);
+        &antsFunUsage(-4,'ff.f','T_VI(<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <range max[m]>)',@_);
 	return T_VI_alt($k,$blen,$sel,$dz,dprime($range_max));        
 }
 
 sub T_VI_alt($$$$$)
 {
 	my($k,$blen,$sel,$dz,$dprime) =
-        &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <tilt d-prime[m]>',@_);
+        &antsFunUsage(-4,'ff.ff','T_VI_alt(<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <tilt d-prime[m]>)',@_);
 	croak("T_VI_alt($k,$blen,$sel,$dz,$range_max): bad parameters\n")
 		unless ($k>=0 && $blen>0 && $sel ne '' && $dz>0 && $range_max>=0);				
 	croak("$0: tilt-dprime outside range [0..$blen]\n")
@@ -206,9 +208,9 @@
 	sub T_w(@)
 	{
 		my($k,$blen,$plen,$dz,$range_max) =
-			&antsFunUsage(5,'fffff',
-				'[vertical wavenumber[rad/s]] <ADCP bin length[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>',
-				\@fc,'k',undef,undef,undef,undef,@_);
+			&antsFunUsage(-4,'ffff',
+				'T_w([vertical wavenumber[rad/s]] <ADCP bin length[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>)',
+				\@fc,'k',undef,undef,undef,@_);
 		croak("T_w($k,$blen,$plen,$dz,$range_max): bad parameters\n")
 			unless ($k>=0 && $blen>0 && $plen>0 && $dz>0 && $range_max>=0);				
 		return T_ravg($k,$blen,$plen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
@@ -227,9 +229,9 @@
 	sub T_w_z(@)
 	{
 		my($k,$blen,$plen,$dz,$range_max) =
-			&antsFunUsage(5,'fffff',
-				'[vertical wavenumber[rad/s]] <ADCP bin size[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>',
-				\@fc,'k',undef,undef,undef,undef,@_);
+			&antsFunUsage(-4,'ffff',
+				'T_w_z([vertical wavenumber[rad/s]] <ADCP bin size[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>)',
+				\@fc,'k',undef,undef,undef,@_);
 		croak("T_w_z($k,$blen,$plen,$dz,$range_max): bad parameters\n")
 			unless ($k>=0 && $blen>0 && $plen>0 && $dz>0 && $range_max>=0);				
 		return T_ravg($k,$blen,$plen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)) * T_fdiff($k,$dz);
diff --git a/libLADCP.pl.orig b/libLADCP.pl.orig
new file mode 100644
--- /dev/null
+++ b/libLADCP.pl.orig
@@ -0,0 +1,240 @@
+#======================================================================
+#                    L I B L A D C P . P L 
+#                    doc: Wed Jun  1 20:38:19 2011
+#                    dlm: Mon Apr  1 16:13:46 2019
+#                    (c) 2011 A.M. Thurnherr
+#                    uE-Info: 27 65 NIL 0 0 70 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Jun  1, 2011: - created
+#	Jul 29, 2011: - improved
+#	Aug 10, 2011: - made correct
+#	Aug 11, 2011: - added "convenient combinations"
+#	Aug 18, 2011: - made buoyancy frequency non-constant in S()
+#	Jan  4, 2012: - improved T_VI to allow correcting even without superensembles
+#	Jan  5, 2012: - removed S(), which is just pwrdens/N^2 (rather than
+#				    pwrdens/N^2/(2pi) as I erroneously thought)
+#	Jan 18, 2012: - added T_VI_alt() to allow assessment of tilt correction extrema
+#	Aug 22, 2012: - added documentation
+#				  - added T_w()
+#	Sep 24, 2012: - made "k" argument default in T_w()
+#	Oct 25, 2012: - renamed T_SM() to T_ASM()
+#	Jun 26. 2013: - added T_w_z()
+#				  - added parameter checks to processing-specific corrections
+#	May 18, 2015: - added pulse length to T_w() and T_w_z()
+#	Apr 25, 2018: - added eps_VKE() parameterization
+#	Apr  1, 2019: - removed plen parameter from T_w() and T_w_z()
+
+require "$ANTS/libvec.pl";
+require "$ANTS/libfuns.pl";
+
+#------------------------------------------------------------------------------
+# VKE parameterization for epsilon
+#
+# NOTES:
+#	- see Thurnherr et al. (GRL 2015)
+#	- calculate eps from p0
+#	- optional second argument allows free choice of parameterization constant
+#	- default value is from paper, which is slightly lower than the current value
+#	  used in [LADCP_VKE], which applies the parameterization only to spectra
+#	  passing a few tests
+#------------------------------------------------------------------------------
+
+sub eps_VKE(@)
+{
+    my($p0,$c) =
+        &antsFunUsage(-1,'.','<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]]',@_);
+	$c = 0.021 unless defined($c);
+    return numberp($p0) ? ($p0 / $c) ** 2 : nan;
+}
+
+#------------------------------------------------------------------------------
+# Spectral corrections for LADCP data
+#
+# NOTES:
+#	- see Polzin et al. (JAOT 2002), Thurnherr (JAOT 2012)
+#	- to correct, multiply power densities (or power, I think) with corrections
+#	- apply to down-/up-cast data only
+#------------------------------------------------------------------------------
+#----------------------------------------------------------------------
+# 1. Corrections for individual data acquisition and processing steps
+#----------------------------------------------------------------------
+
+#------------------------------------------------------------------------------
+# T_ravg(k,blen[,plen])
+#	- correct for range averaging due to finite pulse and finite receive window
+# 	- when called with 2 arguments, bin-length == pulse-length assumed
+#------------------------------------------------------------------------------
+
+sub T_ravg(@)
+{
+    my($k,$blen,$plen) =
+        &antsFunUsage(-2,'ff','<vertical wavenumber[rad/s]> <bin-length[m]> [pulse-length[m]]',@_);
+	$plen = $blen unless defined($plen);        
+    return 1 / sinc($k*$blen/2/$PI)**2 / sinc($k*$plen/2/$PI)**2;
+}
+
+#-------------------------------------------------------------
+# T_fdiff(k,dz)
+#	- correct for first differencing on a grid with dz spacing
+#-------------------------------------------------------------
+
+sub T_fdiff($$)
+{
+    my($k,$dz) =
+        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <differencing interval[m]>',@_);
+    return 1 / sinc($k*$dz/2/$PI)**2;
+}
+
+#------------------------------------------------------------
+# T_interp(k,blen,dz)
+#	- correct for CODAS gridding-with-interpolation algorithm
+#	- ONLY USED IN UH SOFTWARE
+#------------------------------------------------------------
+
+sub T_interp($$$)
+{
+    my($k,$blen,$dz) =
+        &antsFunUsage(3,'fff','<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>',@_);
+    return 1 / sinc($k*$blen/2/$PI)**4 / sinc($k*$dz/2/$PI)**2;
+}
+
+#-------------------------------------------------------------------------
+# T_binavg(k,dz)
+#	- correct for simple bin averaging
+#	- Polzin et al. suggest to use blen instead of dz; this must be a typo
+#-------------------------------------------------------------------------
+
+sub T_binavg($$)
+{
+    my($k,$dz) =
+        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <grid resolution[m]>',@_);
+    return 1 / sinc($k*$dz/2/$PI)**2;
+}
+
+#--------------------------------------------------------------------------------
+# T_tilt(k,d')
+#	- d' is a length scale that depends on tilt stats and range max
+#	- on d' == 0, T_tilt() == 1, i.e. the correction is disabled
+#	- d' = dprime(range_max)
+#			- is from a quadratic fit to three data points given by Polzin et al.
+#			- see Thurnherr (J. Tech. 2012) for notes
+#			- on range_max == 0, d' == 0, i.e. the correction is disabled
+#--------------------------------------------------------------------------------
+
+sub T_tilt($$)
+{
+    my($k,$dprime) =
+        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <d-prime[m]>',@_);
+    return $dprime ? 1 / sinc($k*$dprime/2/$PI)**2 : 1;
+}
+
+sub dprime($)
+{
+	return $_[0] ? -1.2 + 0.0857 * $_[0] - 0.000136 * $_[0]**2 : 0;
+} 
+
+#======================================================================
+# 2. Processing-Specific Corrections
+#======================================================================
+
+#----------------------------------------------------------------------
+# T_UH(k,blen,dz,range_max)
+#	- UH implementation of the shear method (WOCE standard)
+#	- range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+sub T_UH($$$$)
+{
+	my($k,$blen,$dz,$range_max) =
+        &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+	croak("T_UH($k,$blen,$dz,$range_max): bad parameters\n")
+		unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				
+    return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_interp($k,$blen,$dz) * T_tilt($k,dprime($range_max));
+}
+
+#----------------------------------------------------------------------
+# T_ASM(k,blen,dz,range_max)
+#	- re-implemented shear method with simple depth binning
+#	- range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+sub T_ASM($$$$)
+{
+	my($k,$blen,$dz,$range_max) =
+        &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+	croak("T_ASM($k,$blen,$dz,$range_max): bad parameters\n")
+		unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				
+    return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
+}
+
+#------------------------------------------------------------
+# T_VI(k,blen,preavg_dz,grid_dz,range_max)
+# T_VI_alt(k,blen,preavg_dz,grid_dz,dprime)
+#	- velocity inversion method of Visbeck (J. Tech., 2002)
+#	- only valid if pre-averaging into superensembles is used
+#	- range_max == 0 disables tilt correction
+#	- sel == nan disables pre-averaging correction
+#------------------------------------------------------------
+
+sub T_VI($$$$$)
+{
+	my($k,$blen,$sel,$dz,$range_max) =
+        &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <range max[m]>',@_);
+	return T_VI_alt($k,$blen,$sel,$dz,dprime($range_max));        
+}
+
+sub T_VI_alt($$$$$)
+{
+	my($k,$blen,$sel,$dz,$dprime) =
+        &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <tilt d-prime[m]>',@_);
+	croak("T_VI_alt($k,$blen,$sel,$dz,$range_max): bad parameters\n")
+		unless ($k>=0 && $blen>0 && $sel ne '' && $dz>0 && $range_max>=0);				
+	croak("$0: tilt-dprime outside range [0..$blen]\n")
+		unless ($dprime>=0 && $dprime<=$blen);
+    return ($sel>0) ? T_ravg($k,$blen) * T_binavg($k,$sel) * T_binavg($k,$dz) * T_tilt($k,$dprime)
+    				: T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,$dprime);
+}
+
+#----------------------------------------------------------------------
+# T_w(k,blen,dz,range_max)
+#	- vertical-velocity method of Thurnherr (IEEE 2011)
+#	- range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+{ my(@fc);
+	sub T_w(@)
+	{
+		my($k,$blen,$dz,$range_max) =
+			&antsFunUsage(4,'ffff',
+				'[vertical wavenumber[rad/s]] <ADCP bin length[m]> <output grid resolution[m]> <range max[m]>',
+				\@fc,'k',undef,undef,undef,undef,@_);
+		croak("T_w($k,$blen,$dz,$range_max): bad parameters\n")
+			unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				
+		return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
+	}
+}
+
+#----------------------------------------------------------------------
+# T_w_z(k,blen,dz,range_max)
+#	- vertical-velocity method of Thurnherr (IEEE 2011)
+#	- first differencing of gridded shear to calculate dw/dz
+#	- NB: grid-scale differentiation assumed
+#	- range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+{ my(@fc);
+	sub T_w_z(@)
+	{
+		my($k,$blen,$dz,$range_max) =
+			&antsFunUsage(4,'ffff',
+				'[vertical wavenumber[rad/s]] <ADCP bin size[m]> <output grid resolution[m]> <range max[m]>',
+				\@fc,'k',undef,undef,undef,undef,@_);
+		croak("T_w_z($k,$blen,$dz,$range_max): bad parameters\n")
+			unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				
+		return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)) * T_fdiff($k,$dz);
+	}
+}
+
+1;
diff --git a/libstats.pl b/libstats.pl
--- a/libstats.pl
+++ b/libstats.pl
@@ -1,9 +1,9 @@
 #======================================================================
-#                    L I B S T A T S . P L 
+#                    . . / L I B / L I B S T A T S . P L 
 #                    doc: Wed Mar 24 13:59:27 1999
-#                    dlm: Sat Jan 31 15:51:14 2015
+#                    dlm: Thu Mar 12 15:11:24 2020
 #                    (c) 1999 A.M. Thurnherr
-#                    uE-Info: 150 0 NIL 0 0 70 2 2 4 NIL ofnI
+#                    uE-Info: 45 44 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -36,10 +36,37 @@
 #				  - added fdiff()
 #	Jan 30, 2015: - added log_avg()
 #				  - added noise_avg()
+#	Mar 26, 2019: - added regress()
 
 require "$ANTS/libfuns.pl";
 
 #----------------------------------------------------------------------
+# simple linear regression
+#	- returns regression coefficient (slope)
+#----------------------------------------------------------------------
+
+sub regress(@)
+{
+	my($n) = @_/2;
+	return nan unless ($n > 1);
+
+	my(@X) = @_[0..$n-1];
+	my(@Y) = @_[$n..2*$n-1];
+
+	my($sumx) = 0;
+	for (my($i)=0; $i<$n; $i++) { $sumx += $X[$i]; }
+    my($midx) = $sumx / @Y;
+    my($rcoeff) = my($sumxdsq) = 0;
+	for (my($i)=0; $i<$n; $i++) {
+		my($xd) = $X[$i] - $midx;
+		$sumxdsq += $xd**2;
+		$rcoeff += $xd * $Y[$i];
+	}
+    $rcoeff /= $sumxdsq;
+	return $rcoeff;
+}
+
+#----------------------------------------------------------------------
 # estimate stderr given stddev & degrees of freedom
 #	- return nan for dof <= 0
 #----------------------------------------------------------------------