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