--- 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/antsusage.pl
+++ b/antsusage.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S U S A G E . P L
# doc: Fri Jun 19 13:43:05 1998
-# dlm: Tue Apr 2 22:26:55 2013
+# dlm: Sat Jan 31 15:49:23 2015
# (c) 1998 A.M. Thurnherr
-# uE-Info: 157 44 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 652 28 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -155,6 +155,8 @@
# Apr 2, 2013: - BUG: pref{}suff special args did sometimes produce unexpanded as well
# as expanded output (unexpanded should be produced only if the
# expansion is empty)
+# Jan 30, 2015: - added &antsFunOpt()
+# Jan 31, 2015: - made it work
# NOTES:
# - ksh expands {}-arguments with commas in them!!! Use + instead
@@ -632,4 +634,23 @@
}
}
+#----------------------------------------------------------------------
+# antsFunOpt(\$opt_x) parses the contents of $opt_x as follows:
+#
+# name(value) => $opt_x = 'name'; $name = "value";
+# name => no change
+#----------------------------------------------------------------------
+
+sub antsFunOpt(@)
+{
+ my($opt) = @_;
+ return unless defined($opt);
+ croak("antsusage.pl: antsFunOpt(@_): argument is not a ref\n")
+ unless ref($opt);
+ my($name,$param) = (${$opt} =~ m{^([^\)]+)\(([^\)]+)\)$});
+ return unless defined($param);
+ eval(sprintf('$%s = "%s";',$name,$param));
+ ${$opt} = $name;
+}
+
1; # return true
--- 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: Mon Oct 6 09:48:22 2014
+# dlm: Tue Nov 18 12:42:30 2014
# (c) 2011 A.M. Thurnherr
-# uE-Info: 21 10 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 22 53 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -19,6 +19,7 @@
# 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";
@@ -91,24 +92,24 @@
sub Sw(@)
{
- my($omega,$m,$lat,$N) = &antsFunUsage(-3,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_);
- my($GM_b) = 1300; #m # pycnocline lengthscale
+ 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*>',@_);
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($GM_jstar,$N,$omega);
- return $GM_E0 * $GM_b * 2 * $f**2/$omega**2/B($omega) * $GM_jstar / ($m+$mstar)**2;
+ 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($GM_jstar,$N);
- return $pi * $GM_E0 * $GM_b * $N * $f * $GM_jstar / ($m+$mstar)**2;
+ 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;