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});
 		}