whoosher version
authorA.M. Thurnherr <athurnherr@yahoo.com>
Thu, 05 Mar 2015 12:51:48 +0000
changeset 10 3dfa16523886
parent 9 1a7983cbb82a
child 15 ebd8a4ddd7f2
whoosher version
INDEX
README
antsusage.pl
libGM.pl
libSBE.pl
libstats.pl
--- 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;