# HG changeset patch
# User A.M. Thurnherr <athurnherr@yahoo.com>
# Date 1470406228 14400
# Node ID 515b06dae59c32a9849b0190722d40ed393a9ba6
# Parent  7c394a2d1fc907143cdfb2e7be67f650dc4fd98e
version at end of ECOGIG EN586 cruise

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: Thu May 26 10:40:58 2016
+                    dlm: Fri Aug  5 10:10:06 2016
                     (c) 2012 A.M. Thurnherr
-                    uE-Info: 124 44 NIL 0 0 72 3 2 4 NIL ofnI
+                    uE-Info: 150 32 NIL 0 0 72 3 2 4 NIL ofnI
 ======================================================================
 
 --------------------------------------
@@ -112,9 +112,9 @@
 May 25, 2016:
 	- published for LADCP_w V1.3beta1
 
-------------------------
-V1.7 (bin interpolation)
-------------------------
+------------------------------------------------
+V1.7 (bin interpolation; better transformations)
+------------------------------------------------
 
 May 25, 2016:
 	- continue working on bin interpolation [RDI_Coords.pl]
@@ -122,3 +122,29 @@
 May 26, 2016:
 	- made it work
 	- updated version in [ADCP_tools_lib.pl]
+
+Jun  6, 2016:
+	- implemented coordinate transformations of Lohrman et al. (JAOT 1990)
+	- PREVIOUS 2-BEAM TRANSFORMATIONS WERE FAIRLY CRUDE APPROXIMATIONS
+	- [RDI_Coords.pl]: sign error in v34
+
+Jun  8, 2016:
+	- minor improvement in [RDI_Coords.pl]
+	- improvements to [editPD0]
+
+Jun  9, 2016:
+	- minor improvements to [listBins]	
+
+Jul  7, 2016:
+	- major BUG: velEarthToBPw() was wrong; new implementation
+	  debugged and verified by Paul Wanis from TRDI
+
+Jul 12, 2016:
+	- improvements to [editPD0]
+
+Jul 26, 2016:
+	- minor improvement to [splitPD0]
+
+Jul 30, 2016:
+	- minor bug in [RDI_PD0_IO.pl]
+	- improvements to [splitPD0]
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 26 17:40:20 2016
+#                    dlm: Sun Jul 31 13:54:26 2016
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 530 14 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 568 0 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # RDI Workhorse Coordinate Transformations
@@ -46,6 +46,13 @@
 #	May 19, 2016: - begin implemeting bin interpolation
 #	May 25, 2016: - continued
 #	May 26, 2016: - made it work
+#	May 30, 2016: - begin implementing 2nd order attitude transformations
+#	Jun  6, 2016: - toEarth transformation in beamToBPEarth was but crude approximation;
+#					updated with transformation taken from Lohrman et al. (JAOT 1990)
+#				  - BUG: v34 sign was inconsistent with RDI coord manual
+#	Jun  8, 2016: - added $ens as arg to velInstrumentToBeam() for consistency
+#	Jul  7, 2016: - added velEarthToBPw() with algorithm debugged and verified
+#					by Paul Wanis from TRDI
 
 use strict;
 use POSIX;
@@ -59,11 +66,12 @@
 # Tweakables
 #----------------------------------------------------------------------
 
-$RDI_Coords::minValidVels = 3;			# 3-beam solutions ok (velBeamToInstrument)
-$RDI_Coords::binMapping = 'linterp';	# 'linterp' or 'none' (earthVels, BPearthVels)
+$RDI_Coords::minValidVels = 3;				# 3-beam solutions ok (velBeamToInstrument)
+$RDI_Coords::binMapping = 'linterp';		# 'linterp' or 'none' (earthVels, BPearthVels)
+$RDI_Coords::beamTransformation = 'LHR90';	# set to 'RDI' to use 1st order transformations from RDI manual
 
 #----------------------------------------------------------------------
-# beam to earth transformation (from RDI coord trans manual)
+# beam to earth transformation 
 #----------------------------------------------------------------------
 
 $RDI_Coords::threeBeam_1 = 0;			# stats from velBeamToInstrument
@@ -123,6 +131,23 @@
 	}
 } # STATIC SCOPE
 
+#----------------------------------------------------------------------
+# velInstrumentToEarth(\%ADCP,ens,v1,v2,v3,v4) => (u,v,w,e)
+#	- $RDI_Coords::beamTransformation = 'LHR90'
+#		- from Lohrmann, Hackett & Roet (J. Tech., 1990)
+#		- eq A1 maps to RDI matrix M (sec 5.6) with
+#			alpha = roll
+#			beta = gimball_pitch
+#			psi = calculation_pitch
+#			psi = asin{sin(beta) cos(alpha) / sqrt[1- sin(alpha)^2 sin(beta)^2]}
+#		- (I only checked for 0 heading, but this is sufficient)
+#	- $RDI_Coords::beamTransformation = 'RDI'
+#		- default prior to LADCP_w V1.3
+#		- from RDI manual
+#		- 99% accurate for p/r<8deg
+#			=> 1cm/s error for 1m/s winch speed!
+#----------------------------------------------------------------------
+
 { # STATIC SCOPE
 	my($hdg,$pitch,$roll,@I2E);
 
@@ -145,9 +170,12 @@
 			$pitch = $ADCP->{ENSEMBLE}[$ens]->{PITCH};
 			$roll  = $ADCP->{ENSEMBLE}[$ens]->{ROLL};
 			my($rad_gimbal_pitch) = atan(tan(rad($pitch)) * cos(rad($roll)));
+			my($rad_calc_pitch) = ($RDI_Coords::beamTransformation eq 'RDI') ? $rad_gimbal_pitch : 
+								  asin(sin($rad_gimbal_pitch)*cos(rad($roll)) /
+									   sqrt(1-sin(rad($roll))**2*sin($rad_gimbal_pitch)**2));
 			my($sh,$ch) = (sin(rad($hdg)),cos(rad($hdg)))
 				if defined($hdg);				
-			my($sp,$cp) = (sin($rad_gimbal_pitch),cos($rad_gimbal_pitch));
+			my($sp,$cp) = (sin($rad_calc_pitch),cos($rad_calc_pitch));
 			my($sr,$cr) = (sin(rad($roll)),	cos(rad($roll)));
 			@I2E = $ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}
 				 ? (
@@ -182,6 +210,9 @@
 #----------------------------------------------------------------------
 # velEarthToInstrument() transforms earth to instrument coordinates
 #	- based on manually inverted rotation matrix M (Sec 5.6 in coord-trans manual)
+#		- Paul Wanis from TRDI pointed out that M is orthonormal, which
+#		  implies that M^-1 = M' (where M' is the transpose), confirming
+#		  the (unnecessary) derivation
 #	- code was verified for both down- and uplookers
 #	- missing heading data (IMP) causes undef beam velocities
 #----------------------------------------------------------------------
@@ -218,7 +249,7 @@
 				  $ev)
 			   : (undef,undef,undef,undef);
 
-	}
+	} # velEarthToIntrument()
 } # STATIC SCOPE
 
 #----------------------------------------------------------------------
@@ -234,7 +265,7 @@
 
 	sub velInstrumentToBeam(@)
 	{
-		my($ADCP,$x,$y,$z,$ev) = @_;
+		my($ADCP,$ens,$x,$y,$z,$ev) = @_;
 		return undef unless (defined($x) + defined($y) +
 					   		 defined($z) + defined($ev) == 4);
 
@@ -260,15 +291,42 @@
 sub velEarthToBeam(@)
 {
 	my($ADCP,$ens,$u,$v,$w,$ev) = @_;
-	return velInstrumentToBeam($ADCP,
+	return velInstrumentToBeam($ADCP,$ens,
 				velEarthToInstrument($ADCP,$ens,$u,$v,$w,$ev));
 }
 
+#----------------------------------------------------------------------
+# velEarthToBPw() returns w12 and w34 for beam-coordinate data
+#	- I am grateful for Paul Wanis from TRDI who corrected a
+#	  bug in my transformation (fixed in V1.3). [The bug did not
+#	  affect the final w profiles significantly, because w12 and w34
+#	  are used only as diagnostics.]
+#	- algorithm:
+#		1) rotate into instrument coordinates
+#		2) w12 = w + e*tan(beam_angle)/sqrt(2)
+#		   w34 = w - e*tan(beam_angle)/sqrt(2)
+#		3) rotate into horizontal coords (earth coords w/o
+#		   considering heading, i.e. same as earth coords
+#		   in case of w
+#----------------------------------------------------------------------
+
+sub velEarthToBPw(@)
+{
+	my($ADCP,$ens,$u,$v,$w,$ev) = @_;
+	my(@iv) = velEarthToInstrument(@_);
+	my(@iv12) = my(@iv34) = @iv;
+	$iv12[2] += $iv[3] * tan(rad($ADCP->{BEAM_ANGLE}))/sqrt(2);
+	$iv34[2] -= $iv[3] * tan(rad($ADCP->{BEAM_ANGLE}))/sqrt(2);
+	my(@ev12) = velInstrumentToEarth($ADCP,$ens,@iv12);
+	my(@ev34) = velInstrumentToEarth($ADCP,$ens,@iv34);
+	return ($ev12[2],$ev34[2]);
+}
+
 #======================================================================
 # velBeamToBPEarth(@) calculates the vertical- and horizontal vels
 # from the two beam pairs separately. Note that (w1+w2)/2 is 
-# identical to the w estimated according to RDI without 3-beam 
-# solutions.
+# identical to the w estimated according to RDI (ignoring 3-beam 
+# solutions).
 #======================================================================
 
 { # STATIC SCOPE
@@ -287,34 +345,35 @@
 			$TwoCosBAngle = 2 * cos(rad($ADCP->{BEAM_ANGLE}));
 			$TwoSinBAngle = 2 * sin(rad($ADCP->{BEAM_ANGLE}));
 		}
-		my($roll)  = rad($ADCP->{ENSEMBLE}[$ens]->{ROLL});							
-		my($sr) = sin($roll); my($cr) = cos($roll);
-		my($pitch) = atan(tan(rad($ADCP->{ENSEMBLE}[$ens]->{PITCH})) * $cr);	# gimbal pitch
-		my($sp) = sin($pitch); my($cp) = cos($pitch);
+		my($rad_roll)  = rad($ADCP->{ENSEMBLE}[$ens]->{ROLL});							
+		my($sr) = sin($rad_roll); my($cr) = cos($rad_roll);
+		my($rad_gimbal_pitch) = atan(tan(rad($ADCP->{ENSEMBLE}[$ens]->{PITCH})) * $cr);	# gimbal pitch
+		my($rad_calc_pitch) = ($RDI_Coords::beamTransformation eq 'RDI') ? $rad_gimbal_pitch :
+							  asin(sin($rad_gimbal_pitch)*cos($rad_roll) /
+								   sqrt(1-sin($rad_roll)**2*sin($rad_gimbal_pitch)**2));
+		my($sp) = sin($rad_calc_pitch); my($cp) = cos($rad_calc_pitch);
 
 		# Sign convention:
 		#	- refer to Coord manual Fig. 3
 		#	- v12 is horizontal velocity from beam1 to beam2, i.e. westward for upward-looking ADCP
 		#	  with beam 3 pointing north (heading = 0)
-		#	- w is +ve upward, regardless of instrument orientation
 
-		my($v12_ic) = ($b1-$b2)/$TwoSinBAngle;							# instrument coords with w vertical up
+		my($v12_ic) = ($b1-$b2)/$TwoSinBAngle;									# instrument coords
+		my($v34_ic) = ($b4-$b3)/$TwoSinBAngle;									# consistent with RDI Coords
 		my($w12_ic) = ($b1+$b2)/$TwoCosBAngle;
-		$w12_ic *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP});
-		my($v34_ic) = ($b3-$b4)/$TwoSinBAngle;
 		my($w34_ic) = ($b3+$b4)/$TwoCosBAngle;
-		$w34_ic *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP});
+		my($w_ic) = ($w12_ic+$w34_ic) / 2;
 	    
-		if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) {				# beampair Earth coords
-			$w12 = $w12_ic*$cr + $v12_ic*$sr - $v34_ic*$sp;
-			$v12 = $v12_ic*$cr - $w12_ic*$sr + $w34_ic*$sp;
-			$w34 = $w34_ic*$cp - $v34_ic*$sp + $v12_ic*$sr;
-    	    $v34 = $v34_ic*$cp + $w34_ic*$sp - $w12_ic*$sr;
-		} else {
-			$w12 = $w12_ic*$cr - $v12_ic*$sr - $v34_ic*$sp;
-			$v12 = $v12_ic*$cr + $w12_ic*$sr + $w34_ic*$sp;
-			$w34 = $w34_ic*$cp - $v34_ic*$sp - $v12_ic*$sr;
-        	$v34 = $v34_ic*$cp + $w34_ic*$sp + $w12_ic*$sr;
+		if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN}) {					# beampair Earth coords
+			$v12 = $v12_ic*$cr 		+ $v34_ic*0			+ $w_ic*$sr;			# Lohrman et al. (1990) A1
+			$v34 = $v12_ic*$sp*$sr	+ $v34_ic*$cp		- $w_ic*$sp*$cr;		#	- defined for z upward => DL
+			$w12 =-$v12_ic*$cp*$sr	+ $v34_ic*$sp		+ $w12_ic*$cp*$cr;
+			$w34 =-$v12_ic*$cp*$sr	+ $v34_ic*$sp		+ $w34_ic*$cp*$cr;
+		} else {																# RDI Coord trans manual
+			$v12 =-$v12_ic*$cr 		+ $v34_ic*0			- $w_ic*$sr;			# 	- as above with 1st & 3rd cols negated
+			$v34 =-$v12_ic*$sp*$sr	+ $v34_ic*$cp		+ $w_ic*$sp*$cr;
+			$w12 = $v12_ic*$cp*$sr	+ $v34_ic*$sp		- $w12_ic*$cp*$cr;
+			$w34 = $v12_ic*$cp*$sr	+ $v34_ic*$sp		- $w34_ic*$cp*$cr;
 		}
 
 		$v12=$w12=undef unless (defined($b1) && defined($b2));
@@ -327,6 +386,8 @@
 #===================================================================
 # velBeamToBPInstrument(@) calculates the instrument-coordinate vels
 # from the two beam pairs separately.
+#	- in spite of the function name, the output is in ship
+#	  coordinates (instr coords with w up)
 #===================================================================
 
 { # STATIC SCOPE
@@ -357,7 +418,7 @@
 			$w12 *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP});
 		}
 		if (defined($b3) && defined($b4)) {
-			$v34 = ($b3-$b4)/$TwoSinBAngle;
+			$v34 = ($b4-$b3)/$TwoSinBAngle;
 			$w34 = ($b3+$b4)/$TwoCosBAngle;
 			$w34 *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP});
 		}
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 _ P D 0 _ I O . P L 
+#                    R D I _ B B _ R E A D . P L 
 #                    doc: Sat Jan 18 14:54:43 2003
-#                    dlm: Mon Feb 29 12:30:04 2016
+#                    dlm: Sat Jul 30 18:34:46 2016
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 1232 63 NIL 0 0 72 10 2 4 NIL ofnI
+#                    uE-Info: 402 62 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # Read RDI BroadBand Binary Data Files (*.[0-9][0-9][0-9])
@@ -75,7 +75,9 @@
 #				  - BUG: most WBPens() error messages used wrong file name
 #	Feb 23, 2016: - changed WBRhdr() to use 2nd ensemble (with correct data-source id)
 #	Feb 26, 2016: - added basic BT data to WBPens(); not BT_RL_* and BT_SIGNAL_STRENGTH
-#	Feb 29, 2016: - LEAP DAY: actually got BT data to work
+#	Feb 29, 2016: - LEAP DAY: actually got BT data patching to work
+#	Jul 30, 2016: - BUG: incomplete last ensemble with garbage content was returned on reading
+#						 WH300 data
 
 # FIRMWARE VERSIONS:
 #	It appears that different firmware versions generate different file
@@ -390,14 +392,14 @@
 	sysread(WBRF,$buf,6) == 6 || die("$WBRcfn: $!");
 	($hid,$did,$dta->{ENSEMBLE_BYTES},$dummy,$dta->{NUMBER_OF_DATA_TYPES})
 		= unpack('CCvCC',$buf);
-	$hid == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Header",$hid,0));
-	$did == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Header",$did,0));
+	$hid == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Header (hid)",$hid,0));
+	$did == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Header (did)",$did,0));
 
 	$start_ens = sysseek(WBRF,$dta->{ENSEMBLE_BYTES}-6+2,1) || die("$WBRcfn: $!");
 	sysread(WBRF,$buf,6) == 6 || die("$WBRcfn: $!");
 	($hid,$did,$dta->{ENSEMBLE_BYTES},$dummy,$dta->{NUMBER_OF_DATA_TYPES})
 		= unpack('CCvCC',$buf);
-	$hid == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Header",$hid,0));
+	$hid == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Header (hid2)",$hid,0));
 	$dta->{DATA_SOURCE_ID} = $did;
 	if ($did == 0x7f) {
 		$dta->{PRODUCER} = 'TRDI ADCP';
@@ -696,10 +698,13 @@
 		# final ensemble.
 
 		sysseek(WBRF,$start_ens,0) || die("$WBRcfn: $!");
-		sysread(WBRF,$buf,$ens_length) == $ens_length || last;
+		unless ((sysread(WBRF,$buf,$ens_length) == $ens_length) &&
+				(sysread(WBRF,$cs,2) == 2)) {
+			pop(@{$E});
+			last;
+		}
 
-		sysread(WBRF,$cs,2) == 2 || last;
-		last unless (unpack('%16C*',$buf) == unpack('v',$cs));
+		pop(@{$E}),last unless (unpack('%16C*',$buf) == unpack('v',$cs));
 
 		#------------------------------
 		# Variable Leader
@@ -782,7 +787,9 @@
 			${$E}[$ens]->{SECONDS} += $B4/100;
 		}
 
-		pop(@{$E}),last if (${$E}[$ens]->{MONTH}>12);						# 10/15/2014; IWISE#145 UL ???
+# 		THE FOLLOWING LINE OF CODE WAS REMOVED 7/30/2016 WHEN I ADDED A POP
+#		TO THE last STATEMENT ABOVE (INCOMPLETE ENSEMBLE)
+#		pop(@{$E}),last if (${$E}[$ens]->{MONTH}>12);						# 10/15/2014; IWISE#145 UL ???
 
 		if ($fixed_leader_bytes == 58) {									# Explorer DVL
 			sysread(WBRF,$buf,14) == 14 || die("$WBRcfn: $!");
@@ -1187,7 +1194,7 @@
 		my($nxt);
 		for ($nxt=6; $nxt<$ndt; $nxt++) {										# scan until BT found
 			sysseek(WBPF,$start_ens+$WBPofs[$nxt],0) || die("$WBPcfn: $!");
-			sysread(WBPF,$buf,2) == 2 || die("$WBPcfn: $!");
+			sysread(WBPF,$buf,2) == 2 || die("$WBPcfn: $! [ens=$ens]");
 			$id = unpack('v',$buf);
 			last if ($id == 0x0600);
 		}
diff --git a/RDI_Utils.pl b/RDI_Utils.pl
--- a/RDI_Utils.pl
+++ b/RDI_Utils.pl
@@ -1,9 +1,9 @@
 #======================================================================
 #                    R D I _ U T I L S . P L 
 #                    doc: Wed Feb 12 10:21:32 2003
-#                    dlm: Thu May 19 10:23:48 2016
+#                    dlm: Sat Jul 30 09:46:59 2016
 #                    (c) 2003 A.M. Thurnherr
-#                    uE-Info: 56 62 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 438 0 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # miscellaneous RDI-specific utilities
diff --git a/editPD0 b/editPD0
--- a/editPD0
+++ b/editPD0
@@ -2,9 +2,9 @@
 #======================================================================
 #                    E D I T P D 0 
 #                    doc: Mon Nov 25 20:24:31 2013
-#                    dlm: Tue Apr 12 21:45:31 2016
+#                    dlm: Tue Jul 12 18:56:55 2016
 #                    (c) 2013 A.M. Thurnherr
-#                    uE-Info: 235 46 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 118 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # edit RDI PD0 file, e.g. to replace pitch/roll/heading with external values
@@ -21,20 +21,24 @@
 #		h(<heading>)			set heading alue value of current ensemble
 #
 #		swap_beams(<b1>,<b2>)	swap data from beams b1 and b2
-#								input in beam coords required
-#								beam rotation is equivalent to 3 consecutive beam swaps
-#								basic BT data are swapped as well (not RL and not SIGNAL_STRENGTH)
+#									- input in beam coords required
+#									- beam rotation is equivalent to 3 consecutive beam swaps
+#									- basic BT data are swapped as well (not RL and not SIGNAL_STRENGTH)
 #
 #		earth2beam()			transform beam to earth coordinates
-#								does not handle bin-remapping
-#								input in beam coords required
+#									- does not handle bin-remapping
+#									- input in earth coords required
+#
+#		beam2earth()			transform earth to beam coordinates
+#									- does not handle bin-remapping
+#									- input in beam coords required
 #
 #		instrument2beam()		transform instrument to earth coordinates
-#								does not handle bin-remapping
-#								input in instrument coords required
+#									- does not handle bin-remapping
+#									- input in instrument coords required
 #
-#		ensure_UL()				overwrite transducer-orientation flag in data file
-#		ensure_DL()
+#		ensure_UL()				correct data for wrong transducer orientation
+#		ensure_DL()					- sets correct flag & negates roll value
 #
 #	- -x notes:
 #		- multiple perl expressions can be combined with ,
@@ -60,6 +64,12 @@
 #	Feb 26, 2016: - added basic BT data to swap_beams()
 #				  - added earth2beam()
 #	Apr 12, 2016: - added instrument2beam()
+#	Jun  3, 2016: - added beam2earth()
+#				  - BUG: instrument2earth() set wrong flag
+#	Jun  8, 2016: - adapted to new interface of velInstrumentToBeam()
+#				  - added %-good to beam2earth and earth2beam
+#				  - made single-ping ensemble requirement for most routines
+#	Jul 12, 2016: - updated ensure_{DL,UL} routines
 
 use Getopt::Std;
 
@@ -99,16 +109,31 @@
 #
 # override transducer orientation
 #
+#	These routines are intended to correct ADCP data for
+#	erroneous orientation switch readings, primarily because
+#	of a stuck switch. While not fully debugged, negating
+#	the roll value greatly improves the vertical velocity
+#	solutions of 2007 CLIVAR I08S profile #1. (#2-#7 could
+#	also be used for testing)
+#
+
 sub ensure_DL()
 {
-	$dta{ENSEMBLE}[$e]->{XDUCER_FACING_DOWN} = 1;
-	$dta{ENSEMBLE}[$e]->{XDUCER_FACING_UP} = undef;
+	if ($dta{ENSEMBLE}[$e]->{XDUCER_FACING_UP}) {
+		$dta{ENSEMBLE}[$e]->{ROLL} *= -1;
+		$dta{ENSEMBLE}[$e]->{XDUCER_FACING_DOWN} = 1;
+		$dta{ENSEMBLE}[$e]->{XDUCER_FACING_UP} = undef;
+	}
 	return 1;
 }
+
 sub ensure_UL()
 {
-	$dta{ENSEMBLE}[$e]->{XDUCER_FACING_DOWN} = undef;
-	$dta{ENSEMBLE}[$e]->{XDUCER_FACING_UP} = 1;
+	if ($dta{ENSEMBLE}[$e]->{XDUCER_FACING_DOWN}) {
+		$dta{ENSEMBLE}[$e]->{ROLL} *= -1;
+		$dta{ENSEMBLE}[$e]->{XDUCER_FACING_UP} = 1;
+		$dta{ENSEMBLE}[$e]->{XDUCER_FACING_DOWN} = undef;
+	}
 	return 1;
 }
 
@@ -124,6 +149,8 @@
 
 	die("$ARGV[0]: beam-coordinate data required\n")
 		unless ($dta{BEAM_COORDINATES});
+	die("$ARGV[0]: single-ping ensembles required\n")
+		unless ($dta{PINGS_PER_ENSEMBLE} == 1);
 
 	if ($dta{BT_PRESENT}) {
 		$tmp = $dta{ENSEMBLE}[$e]->{BT_RANGE}[$b1-1];
@@ -178,13 +205,21 @@
 		unless ($checked) {
 			die("$ARGV[0]: earth-coordinate data required\n")
 				unless ($dta{EARTH_COORDINATES});
+			die("$ARGV[0]: single-ping ensembles required\n")
+				unless ($dta{PINGS_PER_ENSEMBLE} == 1);
 			$dta{BEAM_COORDINATES} = 1; undef($dta{EARTH_COORDINATES});
 			$checked = 1;
 		}
 	    
 		for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
-			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]} =
-				velEarthToBeam(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]});
+			if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][3] == 100) {			# 4-beam solution
+               	@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]} =
+					velEarthToBeam(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]});
+				@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin]} = (100,100,100,100);
+			} else {															# 3-beam solution or no solution
+				undef(@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]});
+				@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin]} = (0,0,0,0);
+			}
 		}
 	
 		return 1;
@@ -202,13 +237,15 @@
 		unless ($checked) {
 			die("$ARGV[0]: instrument-coordinate data required\n")
 				unless ($dta{INSTRUMENT_COORDINATES});
+			die("$ARGV[0]: single-ping ensembles required\n")
+				unless ($dta{PINGS_PER_ENSEMBLE} == 1);
 			$dta{BEAM_COORDINATES} = 1; undef($dta{INSTRUMENT_COORDINATES});
 			$checked = 1;
 		}
 	    
 		for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
 			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]} =
-				velInstrumentToBeam(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]});
+				velInstrumentToBeam(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]});
 		}
 	
 		return 1;
@@ -226,7 +263,9 @@
 		unless ($checked) {
 			die("$ARGV[0]: instrument-coordinate data required\n")
 				unless ($dta{INSTRUMENT_COORDINATES});
-			$dta{BEAM_COORDINATES} = 1; undef($dta{INSTRUMENT_COORDINATES});
+			die("$ARGV[0]: single-ping ensembles required\n")
+				unless ($dta{PINGS_PER_ENSEMBLE} == 1);
+			$dta{EARTH_COORDINATES} = 1; undef($dta{INSTRUMENT_COORDINATES});
 			$checked = 1;
 		}
 	    
@@ -240,6 +279,38 @@
 
 }
 
+#
+# transform beam to earth coordinates
+#
+{ my($checked);
+
+	sub beam2earth()
+	{
+		unless ($checked) {
+			die("$ARGV[0]: beam-coordinate data required\n")
+				unless ($dta{BEAM_COORDINATES});
+			die("$ARGV[0]: single-ping ensembles required\n")
+				unless ($dta{PINGS_PER_ENSEMBLE} == 1);
+			$dta{EARTH_COORDINATES} = 1; undef($dta{BEAM_COORDINATES});
+			$checked = 1;
+		}
+	    
+		for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
+			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]} =
+				velBeamToEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]});
+			$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][0] = 100*$RDI_Coords::threeBeamFlag;	# 3-beam solution
+			$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][1] = 0;								# error velocity not checked
+			$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][2] =									# no solution -> more than 1 bad beam
+								@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$bin]} ? 0 : 100;
+			$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][2] =									# 4-beam solution
+								100 -  $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$bin][0];
+		}
+	
+		return 1;
+	}
+
+}
+
 #--------------------------------------------------
 # Main Routine
 #--------------------------------------------------
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 Mar 17 07:39:43 2016
+#                    dlm: Thu Jun  9 19:09:47 2016
 #                    (c) 2006 A.M. Thurnherr
-#                    uE-Info: 61 0 NIL 0 0 72 10 2 4 NIL ofnI
+#                    uE-Info: 332 0 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # Split data file into per-bin time series.
@@ -58,6 +58,8 @@
 #	Jan 31, 2016: - started debugging the obviously wrong Earth2Beam() transformation
 #	Feb 29, 2016: - continued debugging; removed debugging code
 #   Mar 17, 2016: - adapted to new Getopt library
+#	Jun  9, 2016: - minor improvements
+#				  - BUG: velBeamToEarth() has new interface
 
 # General Notes:
 #	- everything (e.g. beams) is numbered from 1
@@ -149,12 +151,13 @@
 	foreach my $k (keys(%P)) {
 		print(P "$k\{$P{$k}\} ");
 	}
-	if ($beamCoords) {
-		my($pct3b) = ($good_vels[$b] > 0) ? 100*$three_beam[$b]/$good_vels[$b] : nan;
+	my($pct3b);
+#	if ($beamCoords) {
+		$pct3b = ($good_vels[$b] > 0) ? 100*$three_beam[$b]/$good_vels[$b] : nan;
 		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));
-	}
+#	} 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));
@@ -293,7 +296,7 @@
 
 	for (my($b)=0; $b<$dta{N_BINS}; $b++) {
 		if ($beamCoords) {
-			for (my($i)=0; $i<4; $i++) {										# percent-good editing (-p)
+			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]);
@@ -321,12 +324,12 @@
 				velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
 
 			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = 						# calculate earth velocities
-				velBeamToEarth(\%dta,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});
+				velBeamToEarth(\%dta,$e,@{$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
+			unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {		# not a valid transformation
+				undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});
 				next;
 			}
 		} else { 															# Earth coordinates
@@ -339,8 +342,9 @@
 			@{$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
+			$three_beam[$b] += ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]/100 * $dta{PINGS_PER_ENSEMBLE});
+			unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {		# no valid velocity
+				undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});
 				next;
 			}
 
diff --git a/splitPD0 b/splitPD0
--- a/splitPD0
+++ b/splitPD0
@@ -1,10 +1,10 @@
 #!/usr/bin/perl
 #======================================================================
-#                    S P L I T R D I 
+#                    S P L I T P D 0 
 #                    doc: Sat Aug 21 22:20:27 2010
-#                    dlm: Sun Sep 14 17:47:48 2014
+#                    dlm: Sat Jul 30 18:50:43 2016
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 54 0 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 69 60 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # split RDI files based on list of ensemble numbers (e.g. from yoyo -t)
@@ -14,7 +14,11 @@
 #	Jun 24, 2011: - replaced -b, -n by -o
 #   Feb 13, 2014: - updated doc
 #	Mar 19, 2014: - added -s)kip small files
-#	Sep 14, 2014: - added -f)irst 
+#	Sep 14, 2014: - added -f)irst
+#	Jul 26, 2016: - changed file numbering to 1-relative
+#	Jul 30, 2016: - modified -o default
+#				  - added code to set DSID of first ensemble of each
+#					output file to 0x7f7f
 
 # NOTES:
 #   - it is assumed that the input file begins with ensemble #1
@@ -23,7 +27,7 @@
 
 # FILE NAME CONVENTION:
 #   - in order to assign individual yoyo casts numerical station numbers,
-#     by default, the yoyo cast number is inserted after the station number
+#     by default, an underscore and the yoyo cast number is added to the basename
 
 # EXAMPLES:
 #   splitRDI 017DL000.000 `mkProfile 017DL000.000 | yoyo -QFens -ut`
@@ -33,41 +37,50 @@
 use Getopt::Std;
 
 die("Usage: $0 " .
-	"[-o)ut-file <fmt[e.g. 017%02dDL000.000]>] " .
+	"[-o)ut-file <fmt[e.g. 017DL_%02d.000]>] " .
 	"[-f)irst <cast #>] " .
 	"[require -m)in <ens> to produce output] " .
 	"<RDI file> <ens> <ens[...]>\n")
 		unless (&getopts('f:o:m:') && @ARGV>=3);
 
-$opt_o = substr($ARGV[0],0,3)					# default output filename format
-		 . '%02d'
-		 . substr($ARGV[0],-9)
-			unless defined($opt_o);
+unless (defined($opt_o)) {
+	my($bn,$extn) = ($ARGV[0] =~ m{([^/]+)\.([^\.]+)$});
+	$opt_o = "${bn}_%02d.$extn";
+}
 
-$opt_m = 0 unless defined($opt_m);				# default: produce tiny files as well
+$opt_m = 0 unless defined($opt_m);								# default: produce tiny files as well
 	
-readHeader($ARGV[0],\%hdr); shift;				# get length of ensembles
+readHeader($ARGV[0],\%hdr); shift;								# get length of ensembles
 $ens_len = $hdr{ENSEMBLE_BYTES} + 2;
 
-$first_ens = $ARGV[0]+1; shift;					# initialize loop
+$first_ens = $ARGV[0]+1; shift;									# initialize loop
 $last_ens  = $ARGV[0]; shift;
-$cnr = 0;
+$cnr = 1;
 
-do {											# split data
-	sysseek(WBRF,($first_ens-2)*$ens_len,0) ||	# read next block of data
+do {															# split data
+	sysseek(WBRF,($first_ens-2)*$ens_len,0) ||					# begin next block of ensembles
 		die("$WBRcfn: $!");
 	$last_ens++ unless defined($ARGV[0]);
-	$nBytes = ($last_ens-$first_ens+1) * $ens_len;
+
+	sysread(WBRF,$ids,2) || die("$WBRcfn: file truncated");		# read 1st ensemble & ensure DSID is 0x7f
+	$ids = pack('H4','7f7f');									# reset DSID
+	sysread(WBRF,$febuf,$ens_len-4) == $ens_len-4 ||
+		die("$WBRcfn: file truncated");
+	sysread(WBRF,$csum,2) || die("$WBRcfn: file truncated");
+	$csum = pack('v',unpack('%16C*',$ids.$febuf));				# re-calculate checksum
+	
+	$nBytes = ($last_ens-$first_ens) * $ens_len;				# read remaining ensembles in block
 	sysread(WBRF,$buf,$nBytes) == $nBytes ||
 		die("$WBRcfn: file truncated");
 
-	if ($last_ens-$first_ens+1 > $opt_m) {		# produce file only if sufficient # of ensembles
+	if ($last_ens-$first_ens+1 >= $opt_m) {						# write output only if sufficient # of ensembles
 		$fn = sprintf($opt_o,$opt_f+$cnr++);
 		open(F,">$fn") || die("$fn: $!\n");
-		syswrite(F,$buf,$nBytes) == $nBytes ||
+		syswrite(F,$ids.$febuf.$csum.$buf,$nBytes+$ens_len) == $nBytes+$ens_len ||
 			die("$fn: $!\n");
 		close(F);
-	    printf(STDERR "$fn: %d ensembles ($nBytes bytes)\n",$last_ens-$first_ens+1);
+	    printf(STDERR "$fn: %d ensembles (%d bytes)\n",
+						$last_ens-$first_ens+1,$nBytes+$ens_len);
 	}
 	
 	$first_ens = $last_ens+1;