--- 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);
--- 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=
--- 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.
--- 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";
--- 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;
+ }
}
#----------------------------------------------------------------------
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
--- 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;