# HG changeset patch
# User A.M. Thurnherr <athurnherr@yahoo.com>
# Date 1614801051 18000
# Node ID 21cf468fa8e0c0981c4806507a28a63d4f5fc63d
# Parent  51c5988a7f1f428ff50ba9616e4915487f94fe76
prior to A20

diff --git a/RDI_Coords.pl b/RDI_Coords.pl
--- a/RDI_Coords.pl
+++ b/RDI_Coords.pl
@@ -1,9 +1,9 @@
 #======================================================================
 #                    R D I _ C O O R D S . P L 
 #                    doc: Sun Jan 19 17:57:53 2003
-#                    dlm: Wed Mar 28 12:30:12 2018
+#                    dlm: Mon Jun 29 10:59:01 2020
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 109 22 NIL 0 0 72 10 2 4 NIL ofnI
+#                    uE-Info: 61 83 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # RDI Workhorse Coordinate Transformations
@@ -57,6 +57,8 @@
 #	Nov 26, 2017: - BUG: velBeamtoBPEarth() did not respect missing values
 #	Nov 27, 2017: - BUG: numbersp() from [antslib.pl] was used
 #	Mar 28, 2018: - added &loadInstrumentTransformation()
+#	Jun  5, 2020: - added sscorr_w & sscorr_w_mooring
+#   Jun 29, 2020: - added comments for sscorr_w, which conflicts with LADCP_w_ocean
 
 use strict;
 use POSIX;
@@ -673,5 +675,61 @@
 }
 
 #----------------------------------------------------------------------
+# Sound Speed Correction for Vertical Velocity
+#	- Usage: sscorr_w(<observed w>,<beam_angle>,<ADCP sVel setup>,salin,temp,press,<vertical temp. gradient>
+#----------------------------------------------------------------------
+
+sub sscorr_w($$$$$$$$)
+{
+	my($w,$beamAngle,$ssADCP,$salin,$temp,$press,$dz,$dtdz) = @_;
+	
+	return 'nan' unless numberp($w);
+	my($tanSqBeamAngle) = tan(rad($beamAngle))**2;
+	    
+	my($ssXD) 	= sVel($salin,$temp,$press);
+	my($ssBin) 	= sVel($salin,$temp+$dz*$dtdz,$press-$dz);
+	my($Kn) = sqrt(1 + (1 - $ssBin/$ssXD)**2 * $tanSqBeamAngle);
+	return $w * $ssBin/$ssADCP / $Kn;
+}
+
+#----------------------------------------------------------------------
+# Sound Speed Correction for Vertical Velocity
+#	- Usage: sscorr_w_mooring(\$ADCP,<ens-idx>,<bin-idx>,<press>,<salin>,<vertical temp. gradient>)
+#	- Notes:
+#		- RDI Coord. Trans. manual sec. 4.1
+#			- manual error: the ^2 applies to the []
+#		- difference between pressure and depth over instrument range ignored
+#	- Assumptions:
+#		- libEOS83.pl loaded
+#		- $ADCP{ENSEMBLE}[$ens-idx]->VELOCITY}[2] contains measured w
+#		- sound speed variation dominated by temperature
+#		- vertical temperature gradient is constant in time (violated
+#		  at least on superinertial time scales)
+#----------------------------------------------------------------------
+
+sub sscorr_w_mooring($$$$$)
+{												
+	my($ADCP,$ei,$bi,$press,$salin,$dtdz) = @_;
+
+	my($w) = $ADCP->{ENSEMBLE}[$ei]->{VELOCITY}[2];
+	return 'nan' unless numberp($w);
+
+	my($tanSqBeamAngle) = tan(rad($ADCP->{BEAM_ANGLE}))**2;
+	    
+	$Global::P{ITS} = 90;
+	my($ssXD) 	= sVel($salin,$ADCP->{ENSEMBLE}[$ei]->{TEMPERATURE},$press);
+	my($ssADCP)	= $ADCP->{ENSEMBLE}[$ei]->{SPEED_OF_SOUND};
+
+	my($dz)		= $ADCP->{ENSEMBLE}[$ei]->{BLANKING_DISTANCE} + $bi * $ADCP->{ENSEMBLE}[$ei]->{BIN_LENGTH};
+	$dz *= -1 if ($ADCP->{ENSEMBLE}[$ei]->{XDUCER_FACING_DOWN});			# z increases upward
+
+	my($ssBin) 	= sVel($salin,
+					   $ADCP->{ENSEMBLE}[$ei]->{TEMPERATURE} + $dz*$dtdz,
+					   $press-$dz);											# ignore press/depth difference across range
+
+	my($Kn) = sqrt(1 + (1 - $ssBin/$ssXD)**2 * $tanSqBeamAngle);		 	# RDI manual
+	return $w * $ssBin/$ssADCP / $Kn;
+}
+
 
 1;
diff --git a/listBins b/listBins
--- a/listBins
+++ b/listBins
@@ -2,9 +2,9 @@
 #======================================================================
 #                    L I S T B I N S 
 #                    doc: Fri Aug 25 15:57:05 2006
-#                    dlm: Thu Feb 13 10:36:24 2020
+#                    dlm: Fri Jun  5 13:45:23 2020
 #                    (c) 2006 A.M. Thurnherr
-#                    uE-Info: 133 33 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 367 103 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # Split data file into per-bin time series.
@@ -74,6 +74,7 @@
 #	Jun 13, 2018: - adpated to RTI files (disabled BIT error check)
 #				  - BUG: dn did not have sufficient digits
 #	Feb 13, 2020: - added -z
+#	May 11, 2020: - removed -z, added -t -m
 
 # General Notes:
 #	- everything (e.g. beams) is numbered from 1
@@ -118,19 +119,18 @@
 require "$ANTS/ants.pl";
 require "$ANTS/libconv.pl";
 
-die("Usage: $0 [-z) progress dots] [-r)ange <first_ens,last_ens>] [-l)ast <bin>] [-R)enumber ensembles from 1] " .
+die("Usage: $0 [-r)ange <first_ens,last_ens>] [-l)ast <bin>] [-R)enumber ensembles from 1] " .
 			  "[-o)utput <redirection[>bin%d.raw]>] " .
 			  "[output -a)ll ens (not just those with good vels)] " .
 			  "[-M)agnetic <declination>] " .
 			  "[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
 			  "[Instrument -T)ransformation Matrix <file>] " .
+			  '[disable bin -m)apping] [use TRDI beam-to-earth -t)ransformation] ' .
 			  "[-P)itch/Roll <bias/bias>] [-B)eamvel <bias/bias/bias/bias>] " .
 		 	  "[require -4)-beam solutions] [-d)iscard <beam#>] " .
 		 	  "[-p)ct-good <min>] " .
 			  "<RDI file>\n")
-	unless (&getopts("4aB:d:l:M:o:p:r:P:RS:T:z") && @ARGV == 1);
-
-$global::readDataProgress = 10000 if defined($opt_z);
+	unless (&getopts("4aB:d:l:mM:o:p:r:P:RS:tT:") && @ARGV == 1);
 
 ($P{pitch_bias},$P{roll_bias}) = split('[,/]',$opt_P);
 ($P{velbias_b1},$P{velbias_b2},$P{velbias_b3},$P{velbias_b4}) = split('[,/]',$opt_B);
@@ -140,7 +140,9 @@
 
 $opt_p = 0 unless defined($opt_p);
 
-$RDI_Coords::minValidVels = 4 if ($opt_4);			# no 3-beam solutions
+$RDI_Coords::minValidVels 		= 4 if ($opt_4);			# no 3-beam solutions
+$RDI_Coords::binMapping 		= 'none' if ($opt_m);  		# 'linterp' is default
+$RDI_Coords::beamTransformation = 'RDI'  if ($opt_t);		# 'LHR90' is default
 
 print(STDERR "WARNING: magnetic declination not set!\n")
 	unless defined($opt_M);
@@ -362,7 +364,7 @@
 			@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]} =					# save beam velocities
 				@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]};
 
-			@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =				# calculate w12, w34
+			@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =				# calculate v12, w12, v34, w34
 				velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
 
 			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = 						# calculate earth velocities
diff --git a/listEns b/listEns
--- a/listEns
+++ b/listEns
@@ -2,9 +2,9 @@
 #======================================================================
 #                    L I S T E N S 
 #                    doc: Sat Jan 18 18:41:49 2003
-#                    dlm: Thu Feb 13 10:36:50 2020
+#                    dlm: Sun Feb 28 15:46:21 2021
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 96 54 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 351 1 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 $synopsis = 'list ensemble summaries (default), dump ensembles (-E), time-average ensembles (-T)';
@@ -61,6 +61,7 @@
 #	Apr 10, 2018: - added -l)ast bin
 #	May 31, 2018: - BUG: -A was disabled by default
 #	Feb 13, 2020: - added support for $readDataProgress
+#	Feb 19, 2021: - BUG: -T did not handle new years correctly
 
 # Notes:
 #	- -E/-B outputs data in earth coordinates, unless -b is set also
@@ -311,19 +312,18 @@
 } elsif (defined($opt_T)) { 									# time-average ensembles
 
 	my(@tmp) = split(',',$opt_T);								# decode -T 
-	my($Tstart,$deltaT,$Tend,$month,$day,$year);
+	my($Tstart,$deltaT,$Tend,$month,$day);						# NB: $yearbase needs to be global!
 	if (@tmp == 3) {
 		($Tstart,$deltaT,$Tend) = @tmp;
 	} elsif (@tmp == 2) {
 		($Tstart,$deltaT) = @tmp;
-		($month,$day,$year) = split('/',$dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{DATE});
-		$Tend = str2dec_time($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{DATE},$dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{TIME},$year);
+		($month,$day,$yearbase) = split('/',$dta{ENSEMBLE}[0]->{DATE});
+		$Tend = str2dec_time($dta{ENSEMBLE}[0]->{DATE},$dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{TIME},$yearbase);
 	} else {
-		($month,$day,$year) = split('/',$dta{ENSEMBLE}[0]->{DATE});
-		$Tstart = str2dec_time($dta{ENSEMBLE}[0]->{DATE},$dta{ENSEMBLE}[0]->{TIME},$year);
+		($month,$day,$yearbase) = split('/',$dta{ENSEMBLE}[0]->{DATE});
+		$Tstart = str2dec_time($dta{ENSEMBLE}[0]->{DATE},$dta{ENSEMBLE}[0]->{TIME},$yearbase);
 		($deltaT) = @tmp;
-		($month,$day,$year) = split('/',$dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{DATE});
-		$Tend = str2dec_time($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{DATE},$dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{TIME},$year);
+		$Tend = str2dec_time($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{DATE},$dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{TIME},$yearbase);
 	}
 	$Tstart = &{&antsCompileConstExpr($')} if ($Tstart =~ m{^=});
 	$deltaT = &{&antsCompileConstExpr($')} if ($deltaT =~ m{^=});
@@ -346,11 +346,12 @@
 		my($e) = @_;
 		my($b,$i);
 
-		my($month,$day,$year) = ($e >= 0) ? split('/',$dta{ENSEMBLE}[$e]->{DATE}) : (undef,undef,undef);
-		my($dn) = ($e >= 0) ? str2dec_time($dta{ENSEMBLE}[$e]->{DATE},$dta{ENSEMBLE}[$e]->{TIME},$year) : undef;
+		my($dn) = ($e >= 0) ? str2dec_time($dta{ENSEMBLE}[$e]->{DATE},$dta{ENSEMBLE}[$e]->{TIME},$yearbase) : undef;
 
+#		print(STDERR "ens#$e at $dn (cbin = $cbin)\n");
 		if ($e<0 || $dn>=$Tstart+$cbin*$deltaT) {				# dump full bin
 			my($file) = sprintf("$basename.T$fnfmt",$cbin); 	# file name: <basename>.T0001
+#			print(STDERR "dumping average to $file...\n");
 
 			do {												# skip empy bins
 				$cbin++;
@@ -530,6 +531,8 @@
 	}
 } # define output format
 
+die("$yearbase") unless ($yearbase > 1900);
+
 #----------------------------------------------------------------------
 # Loop Over Ensembles
 #----------------------------------------------------------------------