# HG changeset patch
# User A.M. Thurnherr <athurnherr@yahoo.com>
# Date 1452117593 0
# Node ID 7c7da52363c234825b83a8ac68d29318b81c5484
# Parent  1c128aaaca6f3642cd7a7da483b10a03b8889bf6
transport version

diff --git a/ADCP_tools_lib.pl b/ADCP_tools_lib.pl
new file mode 100644
--- /dev/null
+++ b/ADCP_tools_lib.pl
@@ -0,0 +1,21 @@
+#======================================================================
+#                    A D C P _ T O O L S _ L I B . P L 
+#                    doc: Tue Jan  5 10:45:47 2016
+#                    dlm: Tue Jan  5 10:50:18 2016
+#                    (c) 2016 A.M. Thurnherr
+#                    uE-Info: 21 0 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Jan  5, 2015: - created
+
+$ADCP_tools_version = 1.4;		# Jan  5, 2016
+
+die(sprintf("$0: obsolete ADCP_tools V%.1f; V%.1f required\n",
+    $ADCP_tools_version,$ADCP_tools_minVersion))
+        if (!defined($ADCP_tools_minVersion) || $ADCP_tools_minVersion>$ADCP_tools_version);
+
+require "$ADCP_TOOLS/RDI_Coords.pl";
+require "$ADCP_TOOLS/RDI_PD0_IO.pl";
+require "$ADCP_TOOLS/RDI_Utils.pl";
+
diff --git a/HISTORY b/HISTORY
--- a/HISTORY
+++ b/HISTORY
@@ -1,9 +1,9 @@
 ======================================================================
                     H I S T O R Y 
                     doc: Tue May 15 18:04:39 2012
-                    dlm: Wed Nov  4 11:33:11 2015
+                    dlm: Tue Jan  5 13:22:10 2016
                     (c) 2012 A.M. Thurnherr
-                    uE-Info: 29 32 NIL 0 0 72 3 2 8 NIL ofnI
+                    uE-Info: 35 46 NIL 0 0 72 3 2 8 NIL ofnI
 ======================================================================
 
 May 15, 2012:
@@ -26,4 +26,10 @@
 
 Nov  4, 2015:
   - merged with Oct 2 version on Studio desktop, which ignores
-    initial garbage in PD0 files  
+    initial garbage in PD0 files
+
+Jan  5, 2015: V1.4
+	- added [ADCP_tools_lib.pl] with compile-time version control
+	- [RDI_Coords.pl] added &velEarthToBeam()
+	- updated [listBins] to use versioned libs and calc w12 & w34
+	  from earth-coordinate data correctly
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: Thu May 29 09:19:54 2014
+#                    dlm: Tue Jan  5 13:56:55 2016
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 282 0 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 185 37 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # RDI Workhorse Coordinate Transformations
@@ -39,6 +39,7 @@
 #						 heading
 #				  - removed some old debug statements
 #				  - removed unused code from &velBeamToBPInstrument
+#	Jan  5, 2016: - added &velEarthToInstrument(@), &velInstrumentToBeam(@)
 
 use strict;
 use POSIX;
@@ -155,6 +156,79 @@
 	}
 } # STATIC SCOPE
 
+#----------------------------------------------------------------------
+# velEarthToInstrument() transforms earth to instrument coordinates
+#	- based on manually inverted rotation matrix M (Sec 5.6 in coord-trans manual)
+#	- missing heading data (IMP) causes undef beam velocities
+#----------------------------------------------------------------------
+
+{ # STATIC SCOPE
+	my($hdg,$pitch,$roll,@E2I);
+
+	sub velEarthToInstrument(@)
+	{
+		my($dta,$ens,$u,$v,$w,$ev) = @_;
+
+		unless (@E2I) {
+			$hdg   = $dta->{ENSEMBLE}[$ens]->{HEADING} - $dta->{HEADING_BIAS} if defined($dta->{ENSEMBLE}[$ens]->{HEADING});
+			$pitch = $dta->{ENSEMBLE}[$ens]->{PITCH};
+			$roll  = $dta->{ENSEMBLE}[$ens]->{ROLL};
+			my($rad_gimbal_pitch) = atan(tan(rad($pitch)) * cos(rad($roll)));
+			my($sh,$ch) = (sin(rad($hdg)),cos(rad($hdg)))
+				if defined($hdg);				
+			my($sp,$cp) = (sin($rad_gimbal_pitch),cos($rad_gimbal_pitch));
+			my($sr,$cr) = (sin(rad($roll)),	cos(rad($roll)));
+			@E2I = $dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}
+				 ? ([$ch*-$cr+$sh*$sp*-$sr,	 $ch*$sp*-$sr-$sh*-$cr,	$cp*-$sr],
+				    [$sh*$cp,				 $ch*$cp,				$sp		],
+				    [$ch*-$sr-$sh*$sp*-$cr,	-$sh*-$sr-$ch*$sp*-$cr,	$cp*-$cr])
+				 : ([$ch*$cr+$sh*$sp*$sr,	 $ch*$sp*$sr-$sh*$cr,	$cp*$sr	],
+				    [$sh*$cp,				 $ch*$cp,				$sp		],
+				    [$ch*$sr-$sh*$sp*$cr,	-$sh*$sr-$ch*$sp*$cr,	$cp*$cr	]);
+		}
+
+		return defined($dta->{ENSEMBLE}[$ens]->{HEADING})
+			   ? ($u*$E2I[0][0]+$v*$E2I[0][1]+$w*$E2I[0][2],
+				  $u*$E2I[1][0]+$v*$E2I[1][1]+$w*$E2I[1][2],
+				  $u*$E2I[2][0]+$v*$E2I[2][1]+$w*$E2I[2][2],
+				  $ev)
+			   : (undef,undef,undef,undef);
+
+	}
+} # STATIC SCOPE
+
+#----------------------------------------------------------------------
+# velInstrumentToBeam() transforms instrument to beam coordinates
+#	- based on manually solved eq system in sec 5.3 of coord manual
+#	- does not implement bin-remapping
+#	- does not work for 3-beam solutions, as it is not known which
+#	  beam was bad
+#----------------------------------------------------------------------
+
+{ # STATIC SCOPE
+	my($a,$b,$c,$d);
+
+	sub velInstrumentToBeam(@)
+	{
+		my($dta,$x,$y,$z,$ev) = @_;
+		return undef unless (defined($x) + defined($y) +
+					   		 defined($z) + defined($ev) == 4);
+
+		unless (defined($a)) {
+			$a = 1 / (2 * sin(rad($dta->{BEAM_ANGLE})));
+			$b = 1 / (4 * cos(rad($dta->{BEAM_ANGLE})));
+			$c = $dta->{CONVEX_BEAM_PATTERN} ? 1 : -1;
+			$d = $a / sqrt(2);
+		}
+
+		return ( $x/(2*$a*$c) + $z/(4*$b) + $ev/(4*$d),
+				-$x/(2*$a*$c) + $z/(4*$b) + $ev/(4*$d),
+				-$y/(2*$a*$c) + $z/(4*$b) - $ev/(4*$d),
+				 $y/(2*$a*$c) + $z/(4*$b) - $ev/(4*$d));
+
+	}
+} # STATIC SCOPE
+
 #======================================================================
 # velBeamToBPEarth(@) calculates the vertical- and horizontal vels
 # from the two beam pairs separately. Note that (w1+w2)/2 is 
diff --git a/RDI_PD0_IO.pl b/RDI_PD0_IO.pl
--- a/RDI_PD0_IO.pl
+++ b/RDI_PD0_IO.pl
@@ -1,9 +1,9 @@
 #======================================================================
-#                    R D I _ B B _ R E A D . P L 
+#                    R D I _ P D 0 _ I O . P L 
 #                    doc: Sat Jan 18 14:54:43 2003
-#                    dlm: Fri Oct  2 19:17:15 2015
+#                    dlm: Fri Dec 18 18:01:45 2015
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 66 47 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 1106 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # Read RDI BroadBand Binary Data Files (*.[0-9][0-9][0-9])
@@ -64,6 +64,8 @@
 #					incomplete ensemble at the end, which seems to imply that there is
 #				    a garbage final ensemble that passes the checksum test???
 #	Oct  2, 2015: - added &skip_initial_trash()
+#	Dec 18, 2015: - added most data types to WBPofs()
+#				  - BUG: WBPens() requires round() for scaled values
 
 # FIRMWARE VERSIONS:
 #	It appears that different firmware versions generate different file
@@ -942,10 +944,15 @@
 }
 
 #----------------------------------------------------------------------
-# WBPens(nbins,fixed_leader_bytes,^data)
-# 	- patch PD0 file with new data
-#	- file must already exist and have correct structure
-#	- currently only does part of variable leader data, esp. attitude
+# writeData(output_file_name,^data) WBPens(nbins,fixed_leader_bytes,^data)
+#	- writeData() copies file previously read with readData() to output_file_name
+# 	- WBPens() patches new PD0 file with ^data
+#		- ^data is modified!!!!
+#		- output file must already exist and have correct structure
+#		- only subset of data structure is patched:
+#			- Header: Data Source Id
+#			- Var Ldr: Soundspeed, Depth, Heading, Pitch, Roll, Temp, Salin
+#			- Data: Velocity, Correlation, Echo Amp, %-Good, 
 #----------------------------------------------------------------------
 
 sub writeData(@)
@@ -962,6 +969,13 @@
 	                   \@{$dta->{ENSEMBLE}});
 }
 
+sub round(@)
+{
+	return $_[0] >= 0 ? int($_[0] + 0.5)
+					  : int($_[0] - 0.5);
+}
+
+
 sub WBPens($$$)
 {
 	my($nbins,$fixed_leader_bytes,$E) = @_;
@@ -998,22 +1012,24 @@
 	
 		sysseek(WBPF,$start_ens+$WBPofs[1]+12,0) || die("$WBPcfn: $!");
 		
-		${$E}[$ens]->{XDUCER_DEPTH} *= 10;
+		${$E}[$ens]->{XDUCER_DEPTH} = round(${$E}[$ens]->{XDUCER_DEPTH}*10);
 
-		# IMP EXTENSIONS
-		#---------------
+		#-----------------------------
+		# IMP allows for missing value
+		#-----------------------------
+
 		${$E}[$ens]->{HEADING} = defined(${$E}[$ens]->{HEADING})
-							   ? ${$E}[$ens]->{HEADING} * 100
+							   ? round(${$E}[$ens]->{HEADING}*100)
 							   : 0xF000;
 		${$E}[$ens]->{PITCH} = defined(${$E}[$ens]->{PITCH})
-							 ? unpack('S',pack('s',${$E}[$ens]->{PITCH}*100))
+							 ? unpack('S',pack('s',round(${$E}[$ens]->{PITCH}*100)))
 							 : 0x8000;
 		${$E}[$ens]->{ROLL} = defined(${$E}[$ens]->{ROLL})
-						    ? unpack('S',pack('s',${$E}[$ens]->{ROLL}*100))
+						    ? unpack('S',pack('s',round(${$E}[$ens]->{ROLL}*100)))
 						    : 0x8000;
 
 		${$E}[$ens]->{TEMPERATURE} =
-			unpack('S',pack('s',${$E}[$ens]->{TEMPERATURE}*100));
+			unpack('S',pack('s',round(${$E}[$ens]->{TEMPERATURE}*100)));
 
 		sysseek(WBPF,2,1);			# skip built-in test which reads as 0 but is usually undef		
 									# this was found not to matter, but there is no reason to edit
@@ -1031,6 +1047,83 @@
 		my($nw) = syswrite(WBPF,$buf,14);
 		$nw == 14 || die("$WBPcfn: $nw bytes written ($!)");
 
+
+		#--------------------
+		# Velocity Data
+		#--------------------
+
+		sysseek(WBPF,$start_ens+$WBPofs[2]+2,0) || die("$WBRcfn: $!");	# skip velocity data id (assume it is correct)
+		for ($bin=0; $bin<$nbins; $bin++) {
+			for ($beam=0; $beam<4; $beam++) {
+				${$E}[$ens]->{VELOCITY}[$bin][$beam] = defined(${$E}[$ens]->{VELOCITY}[$bin][$beam])
+							   						 ? round(${$E}[$ens]->{VELOCITY}[$bin][$beam]*1000)
+							   						 : 0x8000;
+				$buf = pack('v',unpack('S',pack('s',${$E}[$ens]->{VELOCITY}[$bin][$beam])));
+				my($nw) = syswrite(WBPF,$buf,2);
+				$nw == 2 || die("$WBPcfn: $nw bytes written ($!)");
+			}
+		}
+
+		#--------------------
+		# Correlation Data
+		#--------------------
+
+		sysseek(WBPF,$start_ens+$WBPofs[3]+2,0) || die("$WBRcfn: $!");
+		for ($bin=0; $bin<$nbins; $bin++) {
+			for ($beam=0; $beam<4; $beam++) {
+				$buf = pack('C',${$E}[$ens]->{CORRELATION}[$bin][$beam]);
+				my($nw) = syswrite(WBPF,$buf,1);
+				$nw == 1 || die("$WBPcfn: $nw bytes written ($!)");
+			}
+		}
+
+		#--------------------
+		# Echo Intensity Data
+		#--------------------
+
+		sysseek(WBPF,$start_ens+$WBPofs[4]+2,0) || die("$WBRcfn: $!");
+
+		for ($bin=0; $bin<$nbins; $bin++) {
+			for ($beam=0; $beam<4; $beam++) {
+				$buf = pack('C',${$E}[$ens]->{ECHO_AMPLITUDE}[$bin][$beam]);
+				my($nw) = syswrite(WBPF,$buf,1);
+				$nw == 1 || die("$WBPcfn: $nw bytes written ($!)");
+			}
+		}
+
+		#--------------------
+		# Percent Good Data
+		#--------------------
+
+		sysseek(WBPF,$start_ens+$WBPofs[5]+2,0) || die("$WBRcfn: $!");
+
+		for ($i=0,$bin=0; $bin<$nbins; $bin++) {
+			for ($beam=0; $beam<4; $beam++,$i++) {
+				$buf = pack('C',${$E}[$ens]->{PERCENT_GOOD}[$bin][$beam]);
+				my($nw) = syswrite(WBPF,$buf,1);
+				$nw == 1 || die("$WBPcfn: $nw bytes written ($!)");
+			}
+		}
+
+		#-----------------------------------------
+		# Bottom-Track Data
+		#	- scan through remaining data types
+		#-----------------------------------------
+
+		my($nxt);
+		for ($nxt=6; $nxt<$ndt; $nxt++) {										# scan until BT found
+			sysseek(WBPF,$start_ens+$WBPofs[$nxt],0) || die("$WBRcfn: $!");
+			sysread(WBPF,$buf,2) == 2 || die("$WBRcfn: $!");
+			$id = unpack('v',$buf);
+			last if ($id == 0x0600);
+		}
+
+		unless ($nxt == $ndt) {													# BT found
+			sysseek(WBPF,14,1) || die("$WBRcfn: $!");							# skip BT config
+			# NOT YET IMPLEMENTED
+		}
+
+
 		#----------------
 		# Update Checksum
 		#----------------
diff --git a/editPD0 b/editPD0
--- a/editPD0
+++ b/editPD0
@@ -2,15 +2,17 @@
 #======================================================================
 #                    E D I T P D 0 
 #                    doc: Mon Nov 25 20:24:31 2013
-#                    dlm: Tue Nov 26 10:31:05 2013
+#                    dlm: Fri Dec 18 20:56:45 2015
 #                    (c) 2013 A.M. Thurnherr
-#                    uE-Info: 70 1 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 104 1 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # edit RDI PD0 file, e.g. to replace pitch/roll/heading with external values
 
 # HISTORY:
 #	Nov 25, 2013: - created
+#	Dec 18, 2015: - added switch_beams()
+#				  - added -x
 
 $0 =~ m{(.*)/[^/]+}; 
 require "$1/RDI_PD0_IO.pl";
@@ -18,12 +20,12 @@
 
 $USAGE = "$0 @ARGV";
 die("Usage: $0 " .
-	'-e) <edit-file> ' .
+	'-e) <edit-file> | -x) <expr> ' .
 	"<input file> <output file>\n")
-		unless (&Getopts('e:') && @ARGV == 2);
+		unless (&Getopts('e:x:') && @ARGV == 2);
 
-die("$0: -e <edit-file> is required\n")
-	unless (-r $opt_e);
+die("$0: -e <edit-file> or -x <expr> required\n")
+	unless (defined($opt_x) || -r $opt_e);
 
 print(STDERR "Reading $ARGV[0]...");				# read data
 readData($ARGV[0],\%dta);
@@ -31,46 +33,84 @@
 
 #----------------------------------------------------------------------
 
-print(STDERR "Editing according to <$opt_e>...");
+print(STDERR "Editing Data...");				
 
 #--------------------------------------------------
 # INTERFACE
-#	- example <edit_file> entry:
+#	- example <edit_file> entry for external attitude data
 #		162     p(3), r(4), h(3.14)
+#	- example <edit_file> entry for beam switching
+#		*		switch_beams(3,4)
 #--------------------------------------------------
 
 sub p($) { $dta{ENSEMBLE}[$e]->{PITCH} = $_[0]; }
 sub r($) { $dta{ENSEMBLE}[$e]->{ROLL} = $_[0]; }
 sub h($) { $dta{ENSEMBLE}[$e]->{HEADING} = $_[0]; }
 
+sub switch_beams($$)
+{
+	my($b1,$b2) = @_;
+
+#	print(STDERR "\n entering switch_beams($b1,$b2) for ens = $e...");
+
+	for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
+		my($tmp) = $dta{ENSEMBLE}[$e]->{VELOCITY}[$bin][$b1-1];
+		$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin][$b1-1] = $dta{ENSEMBLE}[$e]->{VELOCITY}[$bin][$b2-1];
+		$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin][$b2-1] = $tmp;
+
+		$tmp = $dta{ENSEMBLE}[$e]->{CORRELATION}[$bin][$b1-1];
+		$dta{ENSEMBLE}[$e]->{CORRELATION}[$bin][$b1-1] = $dta{ENSEMBLE}[$e]->{CORRELATION}[$bin][$b2-1];
+		$dta{ENSEMBLE}[$e]->{CORRELATION}[$bin][$b2-1] = $tmp;
+
+		$tmp = $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$bin][$b1-1];
+		$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$bin][$b1-1] = $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$bin][$b2-1];
+		$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$bin][$b2-1] = $tmp;
+
+		$tmp = $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][$b1-1];
+		$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][$b1-1] = $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][$b2-1];
+		$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][$b2-1] = $tmp;
+	}
+
+	return 1;
+}
+
 #--------------------------------------------------
 # Main Routine
 #--------------------------------------------------
 
-open(EF,$opt_e) || die("$opt_e: $!\n");
-while (<EF>) {
-	s/\#.*//;
-	next if m/^\s+$/;
-	my($ens,$expr) = m/^\s*(\d+)\s+(.*)$/;
-
-	my($id) = ($expr =~ m/^([A-Z]+)\s/);							# e.g. PITCH, ROLL, HEADING
-	$expr = sprintf('$dta{ENSEMBLE}[$e]->{%s}',$id)
+if (defined($opt_x)) {
+	push(@EE,'*');
+	my($id) = ($opt_x =~ m/^([A-Z]+)\s/);										# e.g. PITCH, ROLL, HEADING
+	$opt_x = sprintf('$dta{ENSEMBLE}[$e]->{%s}',$id)
 		if defined($id);
-		
-	push(@EE,$ens);
-	push(@EX,$expr);
-}
-close(EF);
+	push(@EX,$opt_x);
+}		
 
-for (local($e)=my($eei)=0; $e<@{$dta{ENSEMBLE}}; $e++) {			# local() needed for p(), r(), h()
-	if ($EE[$eei] == $dta{ENSEMBLE}[$e]->{NUMBER}) {				# match => edit
+if (defined($opt_e)) {
+	open(EF,$opt_e) || die("$opt_e: $!\n");
+	while (<EF>) {
+		s/\#.*//;
+		next if m/^\s+$/;
+		my($ens,$expr) = m/^\s*(\*|\d+)\s+(.*)$/;
+	
+		my($id) = ($expr =~ m/^([A-Z]+)\s/);										# e.g. PITCH, ROLL, HEADING
+		$expr = sprintf('$dta{ENSEMBLE}[$e]->{%s}',$id)
+			if defined($id);
+		    
+		push(@EE,$ens);
+		push(@EX,$expr);
+	}
+	close(EF);
+}
+
+for (local($e)=my($eei)=0; $e<@{$dta{ENSEMBLE}}; $e++) {						# local() needed for p(), r(), h()
+#	print(STDERR "\n\@ens=$EE[$eei]: $EX[$eei]\n");
+	if ($EE[$eei] eq '*' || $EE[$eei] == $dta{ENSEMBLE}[$e]->{NUMBER}) {		# match => edit
 #		print(STDERR "\n\@ens=$EE[$eei]: $EX[$eei]");
-#		print(STDERR "\n\tp($dta{ENSEMBLE}[$e]->{PITCH}), r($dta{ENSEMBLE}[$e]->{ROLL}), h($dta{ENSEMBLE}[$e]->{HEADING})");
 		eval($EX[$eei]) || die("$@ while executing <$EX[$eei]>\n");
-#		print(STDERR "\n\tp($dta{ENSEMBLE}[$e]->{PITCH}), r($dta{ENSEMBLE}[$e]->{ROLL}), h($dta{ENSEMBLE}[$e]->{HEADING})...");
-	} elsif ($EE[$eei] > $dta{ENSEMBLE}[$e]->{NUMBER}) {			# next edit later in file => skip
+	} elsif ($EE[$eei] > $dta{ENSEMBLE}[$e]->{NUMBER}) {						# next edit later in file => skip
 		next;
-	} else {														# need next edit
+	} else {																	# need next edit
 		$eei++;
 		last if ($eei >= @EE);
 		redo;
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: Tue Jun 16 09:28:17 2015
+#                    dlm: Wed Jan  6 16:08:54 2016
 #                    (c) 2006 A.M. Thurnherr
-#                    uE-Info: 53 74 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 153 28 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # Split data file into per-bin time series.
@@ -51,6 +51,10 @@
 #	Nov 24, 2014: - enabled -w always
 #	Mar 22, 2015: - replaced -f by -o (allowing for pipes)
 #	Jun 16, 2015: - BUG: velocity bias code did not respect bad velocities
+#	Jan  5, 2016: - adapted to [ANTS_tools_lib.pl]
+#				  - adapted to calculation of w12, w34 from earth-coordinate data
+#				  - several other changes to the code that should not affect the results
+#	Jan  6, 2015: - -b removed (always output beamvels)
 
 # General Notes:
 #	- everything (e.g. beams) is numbered from 1
@@ -79,10 +83,10 @@
 #	  data, where one of the beams performed much worse than the others
 
 require "getopts.pl";
-$0 =~ m{(.*/)[^/]+};
-require "$1RDI_BB_Read.pl";
-require "$1RDI_Coords.pl";
-require "$1RDI_Utils.pl";
+
+$ADCP_tools_minVersion = 1.4;
+($ADCP_TOOLS) = ($0 =~ m{(.*/)[^/]+});
+require "$ADCP_TOOLS/ADCP_tools_lib.pl";
 
 die("Usage: $0 [-r)ange <first_ens,last_ens>] [-R)enumber ensembles] " .
 			  "[-o)utput <redirection[>bin%d.raw]>] " .
@@ -92,9 +96,8 @@
 			  "[-P)itch/Roll <bias/bias>] [-B)eamvel <bias/bias/bias/bias>] " .
 		 	  "[require -4)-beam solutions] [-d)iscard <beam#>] " .
 		 	  "[-p)ct-good <min>] " .
-		 	  "[output -b)eam coordinates] " .
 			  "<RDI file>\n")
-	unless (&Getopts("4abB:d:M:o:p:r:P:RS:") && @ARGV == 1);
+	unless (&Getopts("4aB:d:M:o:p:r:P:RS:") && @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);
@@ -102,15 +105,12 @@
 die("$0: -4 and -d are mutually exclusive\n")
 	if ($opt_4 && defined($opt_d));
 
-die("$0: -p and -b are mutually exclusive\n")
-	if ($opt_b && defined($opt_p));
-
 $opt_p = 0 unless defined($opt_p);
 
 $RDI_Coords::minValidVels = 4 if ($opt_4);			# no 3-beam solutions
 
 print(STDERR "WARNING: magnetic declination not set!\n")
-	unless defined($opt_M) || defined($opt_b);
+	unless defined($opt_M);
 
 $opt_o = '>bin%d.raw' unless defined($opt_o);
 $ifn = $ARGV[0];
@@ -146,13 +146,14 @@
 	foreach my $k (keys(%P)) {
 		print(P "$k\{$P{$k}\} ");
 	}
-	if ($opt_b) {
-		printf(STDERR "%.0f%% ",100*$good_vels[$b]/($le-$fe+1));
-	} else {
+	if ($beamCoords) {
 		my($pct3b) = ($good_vels[$b] > 0) ? 100*$three_beam[$b]/$good_vels[$b] : nan;
-		printf(STDERR "%.0f%%/%.0f%% ",100*$good_vels[$b]/($le-$fe+1),$pct3b);
-		printf(P "pct_3_beam{%.0f} ",$pct3b);
+		printf(STDERR "%02d:%.0f%%/%.0f%% ",$b+1,100*$good_vels[$b]/($le-$fe+1),$pct3b);
+	} else {
+		printf(STDERR "%02d:%.0f%% ",$b+1,100*$good_vels[$b]/($le-$fe+1));
 	}
+
+	printf(P "pct_3_beam{%.0f} ",$pct3b);
 	printf(P "pct_good_vels{%.0f} ",100*$good_vels[$b]/($le-$fe+1));
 	printf(P "bin{%d}",$b+1);
 	printf(P " soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!');
@@ -167,11 +168,11 @@
 			"{sig_heading} {sig_pitch} {sig_roll} " .
 			"{xmit_current} {xmit_voltage} " .
 			"{temp} " .
-			($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") .
+			"{bv1} {bv2} {bv3} {bv4} {u} {v} {w} {err_vel} " .
 			"{v12} {w12} {v34} {w34} " .
 			"{corr1} {corr2} {corr3} {corr4} " .
 			"{amp1} {amp2} {amp3} {amp4} " .
-			($opt_b ? "{pcg1} {pcg2} {pcg3} {pcg4}" : "{3_beam} {min_pcg}")
+			"{pcg1} {pcg2} {pcg3} {pcg4} {3_beam} {min_pcg}"
 	);
 	print(P " {dz}") if ($variable_ssCorr);
 	print(P "\n");
@@ -196,6 +197,7 @@
 		print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE} ");
 		print(P "$dta{ENSEMBLE}[$e]->{TEMPERATURE} ");
 		if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {
+			printf(P "%g %g %g %g ",@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
 			printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] * $ssCorr);
 			printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] * $ssCorr);
 			printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] * $ssCorr);
@@ -218,11 +220,14 @@
 		}
 		print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} ");
 		print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} ");
-		if ($opt_b) {
+		if ($beamCoords) {
 			print(P "@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]} ");
-		} else {
 			printf(P "%d ",$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]);
 			printf(P "%s ",min(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}));
+		} else {
+			print(P "nan nan nan nan ");
+			print(P "$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0] ");
+			print(P "nan ");
 		}
 		printf(P "%g ",$dz[$b]*$ssCorr) if ($variable_ssCorr);
 		print(P "\n");
@@ -244,10 +249,14 @@
 if ($dta{BEAM_COORDINATES}) {						# coords
 	$beamCoords = 1;
 } else {
-	die("$0: -b requires input in beam coordinates\n")
-		if ($opt_b);
-	die("$ifn: only beam and earth coordinates implemented so far\n")
+	die("$ifn: only beam and earth coordinates supported\n")
 		if (!$dta{EARTH_COORDINATES});
+	die("$ifn: -p requires beam-coordinate data\n")
+		if ($opt_p > 0);
+	die("$ifn: -d requires beam-coordinate data\n")
+		if defined($opt_d);
+	die("$ifn: -B requires beam-coordinate data\n")
+		if defined($opt_B);
 }
 
 for (my($b)=0; $b<$dta{N_BINS}; $b++) {				# calc dz
@@ -280,34 +289,63 @@
         if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
 
 	for (my($b)=0; $b<$dta{N_BINS}; $b++) {
-		$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] -= $P{velbias_b1} if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]);
-		$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] -= $P{velbias_b2} if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]);
-		$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] -= $P{velbias_b3} if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]);
-		$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] -= $P{velbias_b4} if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]);
-		if (defined($opt_d)) {
-			undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
-			undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
+		if ($beamCoords) {
+			for (my($i)=0; $i<4; $i++) {										# percent-good editing (-p)
+				if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i] < $opt_p) {
+					undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i]);
+					undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]);
+				}
+	        }
+
+			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] -= $P{velbias_b1}			# beam-velocity biases (-B)
+				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]);
+			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] -= $P{velbias_b2}
+				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]);
+			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] -= $P{velbias_b3}
+				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]);
+			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] -= $P{velbias_b4}
+				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]);
+
+			if (defined($opt_d)) {											# discard data from given beam (-d)
+				undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
+				undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
+			}
+
+			@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]} =					# save beam velocities
+				@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]};
+
+			@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =				# calculate w12, w34
+				velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
+
+			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = 						# calculate earth velocities
+				velInstrumentToEarth(\%dta,$e,
+					velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]})
+			  	);
+			$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag;
+			$three_beam[$b] += $RDI_Coords::threeBeamFlag;
+
+			unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {
+				undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});			# not sure when this can happen
+				next;
+			}
+		} else { 															# Earth coordinates
+			@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]} =					# calculate beam velocities
+				velInstrumentToBeam(\%dta,
+					&velEarthToInstrument(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}));
+				                                            
+			@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =				# calculate w12, w34
+				velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
+
+			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = 						# correct for heading bias
+				velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
+
+			unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {
+				undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});			# not sure when/if this can happen
+				next;
+			}
+
 		}
-		@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =
-			velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
-		for (my($i)=0; $i<4; $i++) {
-			if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i] < $opt_p) {
-				undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i]);
-				undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]);
-			}
-        }
-		@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords
-			? velInstrumentToEarth(\%dta,$e,
-				  velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]})
-			  )
-			: velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]})
-				unless ($opt_b);
-		unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {
-			undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});
-			next;
-		}
-		$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag;
-		$three_beam[$b] += $RDI_Coords::threeBeamFlag;
+
 		$dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1;
 		$good_vels[$b]++; 
 		$lastGoodBin = $b if ($b > $lastGoodBin);
@@ -329,15 +367,15 @@
 print( STDERR "Start      : $dta{ENSEMBLE}[$fe]->{DATE} $dta{ENSEMBLE}[$fe]->{TIME}\n");
 print( STDERR "End        : $dta{ENSEMBLE}[$le]->{DATE} $dta{ENSEMBLE}[$le]->{TIME}\n");
 printf(STDERR "Bins       : %d-%d\n",$firstBin+1,$lastBin+1);
-if ($opt_b) {
-	print(STDERR "Good       : ");
-} else {
+if ($beamCoords) {
 	printf(STDERR "3-Beam     : %d %d %d %d\n",$RDI_Coords::threeBeam_1,
 											   $RDI_Coords::threeBeam_2,
 											   $RDI_Coords::threeBeam_3,
 											   $RDI_Coords::threeBeam_4);
 	
 	print(STDERR "Good/3-beam: ");
+} else {
+	print(STDERR "Good       : ");
 }
 for ($b=$firstBin; $b<=$lastBin; $b++) {				# generate output
 	dumpBin($b,$fe,$le);