RDI_Coords.pl
 changeset 36 515b06dae59c parent 35 7c394a2d1fc9 child 39 3bddaa514ef5
```--- 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

+#----------------------------------------------------------------------
+#	- \$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 @@
if defined(\$hdg);
? (
@@ -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(@)
{
return undef unless (defined(\$x) + defined(\$y) +
defined(\$z) + defined(\$ev) == 4);

@@ -260,15 +291,42 @@
sub velEarthToBeam(@)
{
}

+#----------------------------------------------------------------------
+# 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(@iv) = velEarthToInstrument(@_);
+	my(@iv12) = my(@iv34) = @iv;
+	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 @@
}
-		my(\$sr) = sin(\$roll); my(\$cr) = cos(\$roll);
-		my(\$sp) = sin(\$pitch); my(\$cp) = cos(\$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 @@