# HG changeset patch
# User A.M. Thurnherr <athurnherr@yahoo.com>
# Date 1425578194 18000
# Node ID 4097040ad31c1868f946ae349eed0880d192e506
# Parent  bdd9251865629bd058c367d1a76cb3327f8205d6
?

diff --git a/ANTSlib b/ANTSlib
--- a/ANTSlib
+++ b/ANTSlib
@@ -2,34 +2,42 @@
 #======================================================================
 #                    A N T S L I B 
 #                    doc: Wed May 16 06:19:16 2012
-#                    dlm: Wed May 16 06:28:16 2012
+#                    dlm: Thu Oct 30 09:27:54 2014
 #                    (c) 2012 A.M. Thurnherr
-#                    uE-Info: 11 36 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 14 39 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
 #	May 16, 2012: - created for V5.0
+#	Oct 29, 2014: - made it work again for V6
+#	Oct 30, 2014: - changed output
+#				  - added version check
 
-($ANTSlib) = ($0 =~ m{^(.*)/[^/]*$});
+($ANTSLIB) = ($0 =~ m{^(.*)/[^/]*$});
+$antsMinLibVersion = 6.0;
 
-require "$ANTSlib/ants.pl";
-require "$ANTSlib/libCPT.pl";
-require "$ANTSlib/libEOS83.pl";
-require "$ANTSlib/libGM.pl";
-require "$ANTSlib/libLADCP.pl";
-require "$ANTSlib/libNODC.pl";
-require "$ANTSlib/libPOSIX.pl";
-require "$ANTSlib/libRWalk.pl";
-require "$ANTSlib/libWOCE.pl";
-require "$ANTSlib/libWOCE_oldstyle.pl";
-require "$ANTSlib/libconv.pl";
-require "$ANTSlib/libfuns.pl";
-require "$ANTSlib/libgamma.pl";
-require "$ANTSlib/libstats.pl";
-require "$ANTSlib/libtides.pl";
-require "$ANTSlib/libubtest.pl";
-require "$ANTSlib/libvec.pl";
+require "$ANTSLIB/ants.pl";
+require "$ANTSLIB/libCPT.pl";
+require "$ANTSLIB/libEOS83.pl";
+require "$ANTSLIB/libGM.pl";
+require "$ANTSLIB/libLADCP.pl";
+require "$ANTSLIB/libNODC.pl";
+require "$ANTSLIB/libPOSIX.pl";
+require "$ANTSLIB/libRWalk.pl";
+require "$ANTSLIB/libWOCE.pl";
+require "$ANTSLIB/libWOCE_oldstyle.pl";
+require "$ANTSLIB/libconv.pl";
+require "$ANTSLIB/libfuns.pl";
+require "$ANTSLIB/libgamma.pl";
+require "$ANTSLIB/libstats.pl";
+require "$ANTSLIB/libtides.pl";
+require "$ANTSLIB/libubtest.pl";
+require "$ANTSLIB/libvec.pl";
 
-chomp($about = `sed -n '/^description =/s/description = //p' $ANTSlib/.hg/hgrc`);
-print(STDERR "$about installed in $ANTSlib\n");
+chomp($about = `sed -n '/^description =/s/description = //p' $ANTSLIB/.hg/hgrc`);
+($aboutVer) = ($about =~ /V(.*)$/);
+die(sprintf("$0: inconsistent version numbers (.hg/hgrc: $aboutVer; \$antsLibVersion: %.1f)\n",$antsLibVersion))
+	unless ($aboutVer == $antsLibVersion);
+
+print("$ANTSLIB: $about\n");
 exit(0);
diff --git a/INDEX b/INDEX
--- a/INDEX
+++ b/INDEX
@@ -1,9 +1,9 @@
 ======================================================================
                     I N D E X 
                     doc: Wed Jun 18 09:46:58 1997
-                    dlm: Tue Nov 19 10:23:18 2013
+                    dlm: Mon Nov  3 12:42:00 2014
                     (c) 1997 Andreas Thurnherr
-                    uE-Info: 49 58 NIL 0 0 72 2 2 4 NIL ofnI
+                    uE-Info: 186 41 NIL 0 0 72 2 2 4 NIL ofnI
 ======================================================================
 
 NOTES:
@@ -173,6 +173,7 @@
 
 =Libraries (for -L)=
 
+[libconv.pl]		conversions
 [libCPT.pl]			GMT cpt file parsing
 [libEOS83.pl]		equation of state
 [libfuns.pl]		chapter 6 from NR: special functions; own stuff
@@ -182,11 +183,11 @@
 [libLADCP.pl]		LADCP-related funs
 [libNODC.pl]!!!		NODC SD2 conversions
 [libPOSIX.pl]		POSIX functions
+[libSBE.pl]			Seabird CTD utilities
 [libstats.pl]		significance tests
 [libubtest.pl]		testing stuff
 [libvec.pl]			vector stuff, including distance on globe
 [libWOCE.pl]		WOCE conversions
-[libconv.pl]		conversions
 
 =Numerical Recipes Routines=
 
diff --git a/README b/README
--- a/README
+++ b/README
@@ -1,13 +1,13 @@
 ======================================================================
                     R E A D M E 
                     doc: Tue May 15 18:10:40 2012
-                    dlm: Wed May 16 06:30:40 2012
+                    dlm: Mon Nov  3 12:40:57 2014
                     (c) 2012 A.M. Thurnherr
-                    uE-Info: 19 0 NIL 0 0 72 3 2 4 NIL ofnI
+                    uE-Info: 10 58 NIL 0 0 72 3 2 4 NIL ofnI
 ======================================================================
 
 This directory contains a set of perl libraries written and copyrighted
-by A.M. Thurnherr. The software is written in perl (V5.12.4 or later
+by A.M. Thurnherr. The software is written in perl (V5.8.8 or later
 should work) and assumed to run in a UN*X environment.
 
 To find out which version you have, run "make" in current directory.
diff --git a/ants.pl b/ants.pl
--- a/ants.pl
+++ b/ants.pl
@@ -2,9 +2,9 @@
 #======================================================================
 #                    A N T S . P L 
 #                    doc: Fri Jun 19 14:01:06 1998
-#                    dlm: Mon Sep 24 12:41:50 2012
+#                    dlm: Thu Oct 30 09:33:41 2014
 #                    (c) 1998 A.M. Thurnherr
-#                    uE-Info: 25 76 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 22 55 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -13,10 +13,16 @@
 #  Jul  5, 2006: - removed `basename`
 #  Jul 19, 2006: - added error if exec($ANTS_PERL) fails
 #  Sep 24, 2012: - added support for $ANTSLIB
+#  Oct 29, 2014: - added $antsLibVersion with compile-time version check
 
 exec($ENV{ANTS_PERL},$0,@ARGV),die("$ENV{ANTS_PERL}: $!")
     if (defined($ENV{ANTS_PERL}) && $^X ne $ENV{ANTS_PERL});
 
+$antsLibVersion = 6.0;
+die(sprintf("$0: obsolete library V%.1f; V%.1f required\n",
+	$antsLibVersion,$antsMinLibVersion))
+		if (!defined($antsMinLibVersion) || $antsMinLibVersion>$antsLibVersion);
+
 if (defined($ANTSLIB)) {							# new style (V5)
 	require "$ANTSLIB/antsusage.pl";
 	require "$ANTSLIB/antsio.pl";
diff --git a/libGM.pl b/libGM.pl
--- a/libGM.pl
+++ b/libGM.pl
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L I B G M . P L 
 #                    doc: Sun Feb 20 14:43:47 2011
-#                    dlm: Sat Apr 20 20:05:51 2013
+#                    dlm: Tue Nov 18 12:42:30 2014
 #                    (c) 2011 A.M. Thurnherr
-#                    uE-Info: 62 0 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 22 53 NIL 0 0 70 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -18,6 +18,8 @@
 #	Aug 23, 2012: - cosmetics?
 #	Sep  7, 2012: - made N0, E0, b, jstar global
 #	Dec 28, 2012: - added allowance for small roundoff error to Sw()
+#	Oct  6, 2014: - made omega optional in Sw()
+#	Nov 18, 2014: - made b & jstar mandatory for Sw()
 
 require "$ANTS/libEOS83.pl";
 
@@ -53,6 +55,21 @@
 # velocity.  Since w WKB-scales inversely with N, the largest signals should
 # be in the abyss where you therefore likely have the best chance of
 # measuring it.
+#
+# E. Kunze (email, Sep 19, 2013):
+#
+# S[w](omega, m) = PI*E0*b*[f*sqrt(omega^2-f^2)/omega]*[jstar/(m+mstar)^2] with
+#
+# S[w](m) = PI*E0*b*N*f*[jstar/(m+mstar)^2]
+#
+# where the nondimensional spectral energy level E0 = 6.3e-5, stratification
+# lengthscale b = 1500 m, jstar = 3, mstar = jstar*PI*N/b/N0, and N0 =
+# 5.3e-3 rad/s.
+#
+# NOTES:
+#	- b=1500m is a likely typo, as Gregg & Kunze (1991) have b=1300m
+#	- k_z == m
+#    
 #======================================================================
 
 sub m($$)	# vertical wavenumber as a function of mode number & stratification params
@@ -61,7 +78,7 @@
 
 	return defined($omega) && ($omega <= $GM_N0)
 		   ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j
-		   : $pi * $j * $N / ($GM_b * $GM_N0);			# valid, except in vicinity of buoyancy turning frequency (p. 285)
+		   : $pi * $j * $N / ($GM_b * $GM_N0);		# valid, except in vicinity of buoyancy turning frequency (Munk 1981, p.285)
 }
 
 sub B($)											# structure function (omega dependence)
@@ -73,20 +90,27 @@
 }
 
 
-sub Sw($$$$)
+sub Sw(@)
 {
-	my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_);
+	my($omega,$m,$lat,$b,$jstar,$N) =
+		&antsFunUsage(-5,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]> <b[m]> <j*>',@_);
 
-	local($f) = abs(&f($lat));
-	$omega += $PRACTICALLY_ZERO if ($omega < $f);
-	$omega -= $PRACTICALLY_ZERO if ($omega > $N);
-	return nan if ($omega < $f || $omega > $N);
-
-	my($GM_b) = 1300; #m								# pycnocline lengthscale
-
-	my($mstar) = &m($GM_jstar,$N,$omega);
-
-	return $GM_E0 * $GM_b * 2 * $f**2/$omega**2/B($omega) * $GM_jstar / ($m+$mstar)**2;
+	if (defined($N)) {									# Sw(omega,m)
+		local($f) = abs(&f($lat));
+		$omega += $PRACTICALLY_ZERO if ($omega < $f);
+		$omega -= $PRACTICALLY_ZERO if ($omega > $N);
+		return nan if ($omega < $f || $omega > $N);
+		my($mstar) = &m($jstar,$N,$omega);
+		return $GM_E0 * $b * 2 * $f**2/$omega**2/B($omega) * $jstar / ($m+$mstar)**2;
+	} else {											# Sw(m)
+		$N = $lat;										# shift arguments to account for missing omega
+		$lat = $m;
+		local($f) = abs(&f($lat));
+		$m = $omega;
+		undef($omega);
+		my($mstar) = &m($jstar,$N);
+		return $pi * $GM_E0 * $b * $N * $f * $jstar / ($m+$mstar)**2;
+	}
 }
 
 #----------------------------------------------------------------------
diff --git a/libSBE.pl b/libSBE.pl
new file mode 100755
--- /dev/null
+++ b/libSBE.pl
@@ -0,0 +1,306 @@
+#======================================================================
+#                    L I B S B E . P L 
+#                    doc: Mon Nov  3 12:42:14 2014
+#                    dlm: Mon Nov  3 19:35:53 2014
+#                    (c) 2014 A.M. Thurnherr
+#                    uE-Info: 83 13 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Nov  3, 2014: - exported from [importCNV]
+
+#----------------------------------------------------------------------
+# fname_SBE2std($)
+#	- standardize field names (also adds correct unit %PARAMs)
+#----------------------------------------------------------------------
+
+sub fname_SBE2std($)
+{
+	$_ = $_[0];
+
+	return 'lat' 		if /^lat/;
+	return 'lon' 		if /^lon/;
+	return 'press'		if /^prDM/;
+	return 'depth'		if /^depSM/;
+	return 'O2' 		if /^sbeox0/;
+	return 'alt_O2' 	if /^sbeox1/;
+	return 'salin' 		if /^sal00/;
+	return 'alt_salin' 	if /^sal11/;
+	return 'elapsed'	if /^timeS/;
+	
+	if (/^t090/) {											# temperatures with different scales
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 90);
+		&antsAddParams('ITS',90); $P{ITS} = 90;
+		return 'temp';
+	} elsif (/^t068/) {
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 68);
+		&antsAddParams('ITS',68); $P{ITS} = 68;
+		return 'temp';
+	}
+		
+	if (/^t190/) {
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 90);
+		&antsAddParams('ITS',90); $P{ITS} = 90;
+		return 'alt_temp';
+	} elsif (/^t168/) {
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 68);
+		&antsAddParams('ITS',68); $P{ITS} = 68;
+		return 'alt_temp';
+	}
+
+	if (m{^c0S/m}) {										# conductivity with different units
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} ne 'S/m');
+		&antsAddParams('cond.unit','S/m'); $P{cond.unit} = 'S/m';
+		return 'cond';
+	} elsif (m{^c0mS/cm}) {
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} != 'mS/cm');
+		&antsAddParams('cond.unit','mS/cm'); $P{cond.unit} = 'mS/cm';
+		return 'cond';
+	}
+		
+	if (m{^c1S/m}) {
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} != 'S/m');
+		&antsAddParams('cond.unit','S/m'); $P{cond.unit} = 'S/m';
+		return 'alt_cond';
+	} elsif (m{^c1mS/cm}) {
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} != 'mS/cm');
+		&antsAddParams('cond.unit','mS/cm'); $P{cond.unit} = 'mS/cm';
+		return 'alt_cond';
+	}
+
+	return $_;
+}
+
+# same as above but leaving names in place (only setting %PARAMs)
+sub fname_SBE($)
+{
+	$_ = $_[0];
+
+	if (/^t090/) {											# temperatures with different scales
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 90);
+		&antsAddParams('ITS',90); $P{ITS} = 90;
+	} elsif (/^t068/) {
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 68);
+		&antsAddParams('ITS',68); $P{ITS} = 68;
+	}
+		
+	if (/^t190/) {
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 90);
+		&antsAddParams('ITS',90); $P{ITS} = 90;
+	} elsif (/^t168/) {
+		croak("$0: inconsistent temperature scales\n")
+			if defined($P{ITS}) && ($P{ITS} != 68);
+		&antsAddParams('ITS',68); $P{ITS} = 68;
+	}
+
+	if (m{^c0S/m}) {										# conductivity with different units
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} ne 'S/m');
+		&antsAddParams('cond.unit','S/m'); $P{cond.unit} = 'S/m';
+	} elsif (m{^c0mS/cm}) {
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} != 'mS/cm');
+		&antsAddParams('cond.unit','mS/cm'); $P{cond.unit} = 'mS/cm';
+	}
+		
+	if (m{^c1S/m}) {
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} != 'S/m');
+		&antsAddParams('cond.unit','S/m'); $P{cond.unit} = 'S/m';
+	} elsif (m{^c1mS/cm}) {
+		croak("$0: inconsistent conductivity units\n")
+			if defined($P{cond.unit}) && ($P{cond.unit} != 'mS/cm');
+		&antsAddParams('cond.unit','mS/cm'); $P{cond.unit} = 'mS/cm';
+	}
+
+	return $_;
+}
+
+#----------------------------------------------------------------------
+# SBE_checkTime($$)
+# 	- make sure all times are (roughly) the same
+#----------------------------------------------------------------------
+
+{ # static scope
+	my($target_month,$target_day,$target_year,$target_time);
+
+sub SBE_checkTime($$)
+{
+	return unless $_[1];
+	my($mo,$dy,$yr,$tm) = split('\s+',$_[0]);
+
+	unless (defined($target_month)) {
+		$target_month = $mo;
+		$target_day   = $dy;
+		$target_year  = $yr;
+		$target_time  = $tm;
+		return;
+	}
+
+	croak("$0: inconsistent dates in header ($target_month $target_day $target_year vs $mo $dy $yr)\n")
+		unless ($target_month eq $mo && $target_day == $dy && $target_year == $yr);
+	croak("$0: inconsistent times in header ($target_time vs $tm)\n")
+		unless (abs(frac_day(split(':',$target_time))-frac_day(split(':',$tm))) < 1/60/24);
+}
+
+} # static scope
+
+#----------------------------------------------------------------------
+# sub SBE_parseHeader(FP,std-field-names,time-check)
+#	- parse header information
+#----------------------------------------------------------------------
+
+sub SBE_parseHeader($$$)
+{
+	my($FP,$sfn,$tc) = @_;
+	my($hdr,$nfields,$nrecs,$deg,$min,$NS,$EW,$lat,$lon,$sampint,$badval,$ftype);
+
+	while (1) { 										# parse header
+		chomp($hdr = <$FP>);
+		$hdr =~ s/\r*$//;
+		die("$0: unexpected EOF (format error)\n") unless defined($hdr);
+		last if ($hdr eq '*END*');
+	    
+		$nfields = $',next if ($hdr =~ /nquan = /); 	# Layout
+		$nrecs = $',next if ($hdr =~ /nvalues = /);
+		if ($hdr =~ /name (\d+) = ([^:]+):/) {
+			$antsNewLayout[$1] = $sfn ? fname_SBE2std($2) : fname_SBE($2);
+			next;
+		}
+		    
+		SBE_checkTime($1,$tc),next 							# sanity time check
+			if ($hdr =~ /NMEA UTC \(Time\) = (.*)/);
+		SBE_checkTime($1,$tc),next
+			if ($hdr =~ /System UpLoad Time = (.*)/);
+	
+		&antsAddParams('CNV_File',$1),next				# selected metadata
+			if ($hdr =~ /FileName = (.*)$/);
+		SBE_checkTime($1,$tc),&antsAddParams('start_time',$1),next
+			if ($hdr =~ /start_time = (.*)/);
+	
+		&antsAddParams('station',$1),next
+			if ($hdr =~ /Station\s*:\s*(.*)/);
+		&antsAddParams('ship',$1),next
+			if ($hdr =~ /Ship\s*:\s*(.*)/);
+		&antsAddParams('cruise',$1),next
+			if ($hdr =~ /Cruise\s*:\s*(.*)/);
+		&antsAddParams('time',$1),next
+			if ($hdr =~ /Time\s*:\s*(.*)/);
+		&antsAddParams('date',$1),next
+			if ($hdr =~ /Date\s*:\s*(.*)/);
+	
+		if (($hdr =~ /Latitude\s*[:=]\s*/) && !defined($lat)) {
+			($deg,$min,$NS) = split(/\s+/,$');
+			croak("$0: cannot decode latitude ($')\n")
+				unless ($NS eq 'N' || $NS eq 'S');
+			$lat = $deg + $min/60;
+			$lat *= -1 if ($NS eq 'S');
+			&antsAddParams('lat',$lat);
+			next;
+		}
+		if (($hdr =~ /Longitude\s*[:=]\s*/) && !defined($lon)) {
+			($deg,$min,$EW) = split(/\s+/,$');
+			croak("$0: cannot decode longitude ($')\n")
+				unless ($EW eq 'E' || $EW eq 'W');
+			$lon = $deg + $min/60;
+			$lon *= -1 if ($EW eq 'W');
+			&antsAddParams('lon',$lon);
+			next;
+		}
+	    
+		if ($hdr =~ /interval = seconds: /) {
+			$sampint = 1*$';
+			&antsAddParams('sampling_frequency',1/$sampint);
+			next;
+		}
+		if ($hdr =~ /interval = decibars: /) {
+			$sampint = 1*$';
+			&antsAddParams('sampling_press_interval',$sampint);
+			next;
+		}
+	    
+		$badval = $',next
+			if ($hdr =~ /bad_flag = /); 
+		$ftype = $',next
+			if ($hdr =~ /file_type = /);    
+	}
+
+	croak("$0: cannot determine file layout\n")
+		unless (@antsNewLayout && defined($nfields) && defined($nrecs));
+	croak("$0: cannot determine missing value\n")
+	    unless defined($badval);
+
+	@antsLayout = @antsNewLayout;
+	return ($nfields,$nrecs,$sampint,$badval,$ftype,$lat,$lon);
+}
+
+#----------------------------------------------------------------------
+# SBEin($$)
+#	- read SBE CTD data
+#----------------------------------------------------------------------
+
+{ my(@dta); my($nextR)=0;										# static scope
+
+sub SBEin($$$$$)
+{
+	my($FP,$ftype,$nf,$nr,$bad) = @_;
+	my(@add);
+
+	splice(@ants_,0,$antsBufSkip);								# shift buffers
+
+	if ($ftype eq 'ascii') {
+		until ($#ants_>=0 && &antsBufFull()) {
+			return undef unless (@add = &antsFileIn($FP));
+			for (my($f)=0; $f<$nf; $f++) {
+				$add[$f] = nan if ($add[$f] == $bad);
+			}
+			push(@ants_,[@add]);
+		}
+	} elsif ($ftype eq 'binary') {
+		unless (@dta) {											# read binary data once
+			my($fbits) = 8 * length(pack('f',0));
+			croak(sprintf("$0: incompatible native CPU float representation (%d instead of 32bits)\n",fbits))
+				unless ($fbits == 32);  
+			my($dta);
+			croak("$0: can't read binary data\n")
+				unless (read($FP,$dta,4*$nf*$nr) == 4*$nf*$nr);
+			print(STDERR "WARNING: extraneous data at EOF\n")
+				unless eof($FP);
+			$dta = pack('V*',unpack('N*',$dta)) 				# big-endian CPU
+				if (unpack('h*', pack('s', 1)) =~ /01/);		# c.f. perlport(1)
+			@dta = unpack("f*",$dta);
+			for ($r=0; $r<$nr; $r++) {
+				for ($f=0; $f<$nf; $f++) {
+					@add[$f] = $dta[$r*$nf+$f] == $bad ? nan : $dta[$r*$nf+$f];
+				}
+	        }
+	    }
+		until ($#ants_>=0 && &antsBufFull()) {					# copy next out
+			return undef unless ($nextR < $nr);
+			@add = $dta[$nextR*$nf..($nextR+1)*$nr+$nf];
+			push(@ants_,[@add]);
+			$nextR++;
+        }
+    } else {
+		croak("$0: unknown file type $ftype\n");
+    }
+
+	return $#ants_+1;											# ok
+}
+
+} # static scope
+
+#----------------------------------------------------------------------
+
+1;																# return true
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 
 #                    doc: Wed Mar 24 13:59:27 1999
-#                    dlm: Sun Nov 24 23:23:11 2013
+#                    dlm: Sat Jan 31 15:51:14 2015
 #                    (c) 1999 A.M. Thurnherr
-#                    uE-Info: 192 1 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 150 0 NIL 0 0 70 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -34,6 +34,8 @@
 #	Oct 15, 2012: - added max_i(), min_i()
 #	Nov 24, 2013: - renamed N to ndata
 #				  - added fdiff()
+#	Jan 30, 2015: - added log_avg()
+#				  - added noise_avg()
 
 require "$ANTS/libfuns.pl";
 
@@ -117,6 +119,50 @@
 	return ($N>0)?$sum/$N:nan;
 }
 
+#----------------------------------------------------------------------
+# log_avg does the average in log space, which may be useful
+# for log-normal distributions. However, when I tired this
+# in response to a reviewer suggestion, the results were
+# noticably but not significantly worse, even though at least
+# some of the distributions looked quite log-normal.
+# It appears, however, that 32 samples are enough to sample
+# the distribution of dissipation adequately.
+#----------------------------------------------------------------------
+
+sub log_avg(@)	# exp(avg(log(x)))
+{
+	my($N) = my($sum) = 0;
+	for (my($i)=0; $i<=$#_; $i++) { $N++,$sum+=log($_[$i]) if (numberp($_[$i])); }
+	return ($N>0)?exp($sum/$N):nan;
+}
+
+#----------------------------------------------------------------------
+# noise_avg(noise_lvl,frac)  calculates an arithmetic average after
+# setting all samples below noise level to a certain fraction of the
+# noise level. 66% was used for the VKE paper.
+# Requires $noise_avg to be set, e.g. with &antsFunOpt()
+#----------------------------------------------------------------------
+
+sub noise_avg(@)
+{
+	my($nlv,$frac) = split(',',$noise_avg);
+	croak("libstats.pl: noise_avg: cannot decode \$noise_avg = <$noise_avg> <$nlv,$frac> (noise level, fraction)\n")
+		unless ($nlv > 0) && ($frac < 1);
+	my($nsamp) = my($sum) = 0;
+	for (my($i)=0; $i<=$#_; $i++) {
+		if (numberp($_[$i])) {
+			$nsamp++;
+			if ($_[$i]  > $nlv) {
+				$sum += $_[$i];
+			} else {
+				$sum += $frac*$nlv;
+			}
+		}
+	}
+	return ($nsamp > 0) ? $sum/$nsamp : nan;
+}
+
+
 sub stddev2(@)		# avg, val, val, val, ...
 {
 	my($N) = my($sum) = 0;