new file mode 100644
--- /dev/null
+++ b/INDEX
@@ -0,0 +1,28 @@
+#======================================================================
+# I N D E X
+# doc: Thu Feb 7 14:21:21 2008
+# dlm: Fri Feb 22 08:39:45 2008
+# (c) 2008 A.M. Thurnherr
+# uE-Info: 22 26 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+=General Utilities=
+
+[listHdr] list header info
+[listEns] list ensembles
+[listBT] list bottom-track data
+[listW] list vertical velocities; BROKEN
+[scanBins] list per-bin stats; BROKEN
+
+
+=Moored-ADCP Utilities=
+
+[RDI2grd] make GMT-compatible GRD from ADCP raw file
+[listBins] make per-bin velocity time series
+[meanProf] time-averaged mean profile
+
+
+=Lowered ADCP Utilities=
+
+[mkProfile] make depth-vs-time profile by integrating vertical velocities
+
new file mode 100755
--- /dev/null
+++ b/RDI2grd
@@ -0,0 +1,235 @@
+#!/usr/local/bin/perl
+#======================================================================
+# R D I 2 G R D
+# doc: Wed Aug 30 11:51:22 2006
+# dlm: Thu Jun 18 07:06:34 2009
+# (c) 2006 A.M. Thurnherr
+# uE-Info: 18 54 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# make GMT grd files from RDI file
+
+# HISTORY:
+# Aug 30, 2006: - created at end of GRAVILUCK cruise
+# Aug 31, 2006: - BUG: ensembles/bins were numbered from 0
+# - added -d)imensional coords
+# Sep 1, 2006: - fiddled with registration
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl]
+# Jun 18, 2009: - BUG: xysize had been called xyside
+
+# NOTES:
+# - regular grids only => no dimensional time axis for data collected
+# in multi-ensemble burst mode!
+
+# TODO:
+# - implement soundspeed corretion
+# - add temporal pre-averaging
+
+require "getopts.pl";
+$0 =~ m{(.*/)[^/]+};
+require "$1RDI_BB_Read.pl";
+require "$1RDI_Coords.pl";
+
+use NetCDF;
+
+sub dumpVar($$$$$)
+{
+ my($var,$units,$long_name,$fname,$dimnum) = @_;
+
+ my($ncId) = NetCDF::create("$opt_b$var.grd",NetCDF::CLOBBER);
+ NetCDF::setfill($ncId,NetCDF::NOFILL); # NetCDF library bug
+
+ my($sid) = NetCDF::dimdef($ncId,'side',2);
+ my($aid) = NetCDF::dimdef($ncId,'xysize',($le-$fe+1)*($lastBin-$firstBin+1));
+
+ my($xrid) = NetCDF::vardef($ncId,'x_range',NetCDF::DOUBLE,[$sid]);
+ my($yrid) = NetCDF::vardef($ncId,'y_range',NetCDF::DOUBLE,[$sid]);
+ my($zrid) = NetCDF::vardef($ncId,'z_range',NetCDF::DOUBLE,[$sid]);
+ my($spid) = NetCDF::vardef($ncId,'spacing',NetCDF::DOUBLE,[$sid]);
+ my($dmid) = NetCDF::vardef($ncId,'dimension',NetCDF::LONG,[$sid]);
+
+ my($zvid) = NetCDF::vardef($ncId,'z',NetCDF::FLOAT,[$aid]);
+
+ NetCDF::attput($ncId,NetCDF::GLOBAL,'title',NetCDF::CHAR,$ARGV[0]);
+ NetCDF::attput($ncId,NetCDF::GLOBAL,'source',NetCDF::CHAR,$usage);
+
+ NetCDF::attput($ncId,$xrid,'units',NetCDF::CHAR,
+ $opt_d ? 'day number' : 'ensemble number');
+ NetCDF::attput($ncId,$yrid,'units',NetCDF::CHAR,
+ $opt_d ? 'm' : 'bin number');
+ NetCDF::attput($ncId,$zrid,'units',NetCDF::CHAR,$units);
+
+ NetCDF::attput($ncId,$zvid,'long_name',NetCDF::CHAR,$long_name);
+ NetCDF::attput($ncId,$zvid,'scale_factor',NetCDF::DOUBLE,1);
+ NetCDF::attput($ncId,$zvid,'add_offset',NetCDF::DOUBLE,0);
+ NetCDF::attput($ncId,$zvid,'node_offset',NetCDF::LONG,0);
+
+ NetCDF::endef($ncId);
+
+ if ($opt_d) { # dimensional grid
+ my($ft) = $dta{ENSEMBLE}[$fe]->{DAYNO}; # x_range(side)
+ my($lt) = $dta{ENSEMBLE}[$le]->{DAYNO};
+ NetCDF::varput1($ncId,$xrid,0,$ft);
+ NetCDF::varput1($ncId,$xrid,1,$lt);
+
+ NetCDF::varput1($ncId,$yrid,0, # y_range(side)
+ $dta{DISTANCE_TO_BIN1_CENTER} + $firstBin*$dta{BIN_LENGTH});
+ NetCDF::varput1($ncId,$yrid,1,
+ $dta{DISTANCE_TO_BIN1_CENTER} + $lastBin*$dta{BIN_LENGTH});
+
+ NetCDF::varput1($ncId,$spid,0,($lt-$ft)/($le-$fe)); # spacing(side)
+ NetCDF::varput1($ncId,$spid,1,$dta{BIN_LENGTH});
+ } else { # non-dim grid
+ NetCDF::varput1($ncId,$xrid,0,$fe+1); # x_range(side)
+ NetCDF::varput1($ncId,$xrid,1,$le+1);
+ NetCDF::varput1($ncId,$yrid,0,$firstBin+1); # y_range(side)
+ NetCDF::varput1($ncId,$yrid,1,$lastBin+1);
+ NetCDF::varput1($ncId,$spid,0,1); # spacing(side)
+ NetCDF::varput1($ncId,$spid,1,1);
+ }
+
+ NetCDF::varput1($ncId,$dmid,0,$le-$fe+1); # dimension(side)
+ NetCDF::varput1($ncId,$dmid,1,$lastBin-$firstBin+1);
+
+ my($min) = 9e99; # z(xyside)
+ my($max) = -9e99;
+ my(@data);
+ for (my($b)=$lastBin; $b>=$firstBin; $b--) {
+ for (my($e)=$fe; $e<=$le; $e++) {
+ my($v) = $dta{ENSEMBLE}[$e]->{$fname}[$b][$dimnum];
+ $v = nan unless defined($v);
+ $min = $v if ($v < $min);
+ $max = $v if ($v > $max);
+ push(@data,$v);
+ }
+ }
+
+ my(@start) = (0);
+ my(@count) = (scalar(@data));
+
+ NetCDF::varput($ncId,$zvid,\@start,\@count,\@data);
+
+ NetCDF::varput1($ncId,$zrid,0,$min); # z_range(side)
+ NetCDF::varput1($ncId,$zrid,1,$max);
+
+ NetCDF::close($ncId);
+}
+
+#------
+# Usage
+#------
+
+$usage = "$0 @ARGV";
+die("Usage: $0 [-M)agnetic <declination>] [-r)ange <first_ens,last_ens>] " .
+ "[output -b)ase <name>] [-d)imensional coordinates] " .
+ "<RDI file>\n")
+ unless (&Getopts("b:dM:r:") && @ARGV == 1);
+
+print(STDERR "WARNING: magnetic declination not set!\n")
+ unless defined($opt_M);
+
+unless (defined($opt_b)) {
+ $opt_b = "$ARGV[0]_";
+ $opt_b =~ m{(.*)\.\d\d\d};
+}
+
+($first_ens,$last_ens) = split(',',$opt_r)
+ if defined($opt_r);
+
+#----------
+# Read Data
+#----------
+
+print(STDERR "Reading $ARGV[0]...");
+readData($ARGV[0],\%dta);
+printf(STDERR "%d complete ensembles\n",scalar(@{$dta{ENSEMBLE}}));
+
+#--------------------------------------------------
+# Find Good Ensemble Range & Make Earth Coordinates
+#--------------------------------------------------
+
+print(STDERR "Checking/transforming data...");
+$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
+
+if ($dta{BEAM_COORDINATES}) { # coords used
+ $beamCoords = 1;
+} elsif (!$dta{EARTH_COORDINATES}) {
+ die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
+}
+
+$lastGoodBin = 0;
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities
+ next if (defined($first_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
+ $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e
+ unless defined($P{first_ens});
+ last if (defined($last_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);
+ $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER};
+ $le = $e;
+
+ die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);
+ die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});
+ die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
+
+ for (my($b)=0; $b<$dta{N_BINS}; $b++) {
+ next unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]));
+ @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords
+ ? velInstrumentToEarth(\%dta,$e,
+ velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]})
+ )
+ : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
+ $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1;
+ $good_vels[$b]++;
+ $lastGoodBin = $b if ($b > $lastGoodBin);
+ $firstGoodEns = $e unless defined($firstGoodEns);
+ $lastGoodEns = $e;
+ }
+}
+
+unless (defined($opt_r)) {
+ $fe = $firstGoodEns;
+ $le = $lastGoodEns;
+}
+
+$firstBin = 0;
+$lastBin = $lastGoodBin;
+
+print(STDERR "\n");
+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);
+
+#-----------
+# Write Data
+#-----------
+
+print(STDERR "Writing NetCDF files");
+&dumpVar('u','[m/s]','eastward velocity' ,'VELOCITY',0); print(STDERR '.');
+&dumpVar('v','[m/s]','northward velocity','VELOCITY',1); print(STDERR '.');
+&dumpVar('w','[m/s]','vertical velocity' ,'VELOCITY',2); print(STDERR '.');
+&dumpVar('e','[m/s]','error velocity' ,'VELOCITY',3); print(STDERR '.');
+
+&dumpVar('ea1','[count]','beam-1 echo amplitude','ECHO_AMPLITUDE',0); print(STDERR '.');
+&dumpVar('ea2','[count]','beam-2 echo amplitude','ECHO_AMPLITUDE',1); print(STDERR '.');
+&dumpVar('ea3','[count]','beam-3 echo amplitude','ECHO_AMPLITUDE',2); print(STDERR '.');
+&dumpVar('ea4','[count]','beam-4 echo amplitude','ECHO_AMPLITUDE',3); print(STDERR '.');
+
+&dumpVar('corr1','[count]','beam-1 correlation','CORRELATION',0); print(STDERR '.');
+&dumpVar('corr2','[count]','beam-2 correlation','CORRELATION',1); print(STDERR '.');
+&dumpVar('corr3','[count]','beam-3 correlation','CORRELATION',2); print(STDERR '.');
+&dumpVar('corr4','[count]','beam-4 correlation','CORRELATION',3); print(STDERR '.');
+
+&dumpVar('pcg1','[count]','beam-1 %-good','PERCENT_GOOD',0); print(STDERR '.');
+&dumpVar('pcg2','[count]','beam-2 %-good','PERCENT_GOOD',1); print(STDERR '.');
+&dumpVar('pcg3','[count]','beam-3 %-good','PERCENT_GOOD',2); print(STDERR '.');
+&dumpVar('pcg4','[count]','beam-4 %-good','PERCENT_GOOD',3); print(STDERR '.');
+print(STDERR "\n");
+
+exit(0);
new file mode 100644
--- /dev/null
+++ b/RDI_BB_Read.pl
@@ -0,0 +1,809 @@
+#======================================================================
+# R D I _ B B _ R E A D . P L
+# doc: Sat Jan 18 14:54:43 2003
+# dlm: Wed Jun 4 09:43:15 2008
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 44 25 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# Read RDI BroadBand Binary Data Files (*.[0-9][0-9][0-9])
+
+# HISTORY:
+# Jan 18, 2003: - incepted aboard the Aurora Australis (KAOS)
+# Jan 19, 2003: - continued
+# Jan 20, 2003: - replaced INTENSITY by AMPLITUDE
+# Jan 21, 2003: - changed heading-correction field names
+# Jan 27, 2003: - cosmetics
+# Feb 14, 2003: - moved BT setup into header
+# Mar 15, 2003: - moved 10th xmit voltage into header as BATTERY field
+# - removed again, because values are not meaningful
+# Feb 24, 2004: - continued aboard Nathaniel B. Palmer (Anslope II)
+# - BUG: ensemble # was wrong on error messages
+# Feb 26, 2004: - removed ESW_ERROR and bitmasking of ERROR_STATUS_WORD
+# Feb 27, 2004: - removed some unused (commented-out) baggage
+# Mar 11, 2004: - BUG: renamed ACD -> ADC
+# Mar 18, 2004: - cosmetics
+# Mar 30, 2004: - rewrote to speed up reading; new version takes
+# ~40% less time
+# Sep 14, 2005: - made BT data optional (NUMBER_OF_DATA_TYPES)
+# - added DATA_FORMAT
+# Sep 15, 2005: - debugged
+# - implement checksum to robustly find EOF
+# - renamed to RDI_BB_Read.pl
+# - BUG: had used POSIX::mktime with wrong month def!
+# Oct 30, 2005: - added WH300 FW16.27 file format
+# - added DATA_FORMAT_VARIANT
+# - changed semantics so that first valid ensemble is
+# in E[0] (instead of E[$ensNo-1])
+# Nov 8, 2005: - replaced UNIXTIME by UNIX_TIME
+# - added SECNO
+# Aug 31: 2006: - added DAYNO
+# Aug 1, 2007: - BUG: typo in monthLength()
+# Sep 18, 2007: - modified readHeader() readDta() WBRhdr() WBRens() to
+# conserve memory (no large arrays as retvals)
+# Jun 4, 2008: - BUG: BB150 code was not considered on Sep 18, 2007
+
+# FIRMWARE VERSIONS:
+# It appears that different firmware versions generate different file
+# structures. Currently (Sep 2005) these routines have been tested
+# with the following firmware versions (as reported by [listhdr]):
+#
+# Firmw. DATA_FORMAT(_VARIANT) Owner Cruise FIXED_LEADER_LENGTH
+#------------------------------------------------------------
+# 05.52 BB150 (1) UH CLIVAR/P16S 42
+# 16.12 WH300 (1) FSU A0304 53
+# 16.21 WH300 (1) LDEO NBP0402 53
+# 16.27 WH300 (2) Nash ? 59
+
+# NOTES:
+# - RDI stores data in VAX/Intel byte order (little endian)
+# - the output data structure does not exactly mirror the file data
+# structure; the header is not stored at all and the fixed leader
+# data are not duplicated in every ensemble
+# - in the RDI files some fields that logically belong into the header
+# or the fixed leader (e.g. BT_MIN_CORRELATION) appear in the
+# ensemble data --- these are not read on input
+# - the field names are generally unabbreviated except for
+# BT (= Bottom Track), RL (= Reference Layer), MIN and MAX
+# - all arrays are 0-referenced, but the ensemble number is not!
+# - a list of filenames can be be passed to readData() so that
+# files split onto several memory cards (typically .000 .001 &c)
+# can be read --- not sure if this works, actually
+# - the RDI manuals are not entirely clear everywhere; I have made
+# guesses in some cases, but they should not affect the main
+# fields of interest
+# - some fields in the fixed leader are not really fixed in a LADCP
+# setting (e.g. xducer orientation); I'v made an educated guess
+# as to which fields to move to the ENS array
+# - all units except pressure are SI, i.e. in m and m/s
+# - I don't understand the ERROR_STATUS_WORD; here's what 3 different
+# instruments returned:
+# 0x88000100 FSU instrument during A0304 (Firmware 16.12)
+# 0x88008180 LDEO uplooker (slave) during NBP0402 (Firmware 16.21)
+# 0x00008180 LDEO downlooker (master) during NBP0402 (Firmware 16.21)
+# According to the manual (January 2001 version) this would, for example,
+# indicate power failures on both FSU and LDEO slave instruments...
+
+# &readData() returns perl obj (ref to anonymous hash) with the following
+# structure:
+#
+# DATA_FORMAT string BB150, WH300
+# DATA_FORMAT_VARIANT scalar ?
+# NUMBER_OF_DATA_TYPES scalar 6 (no BT) or 7
+# ENSEMBLE_BYTES scalar ?, number of bytes w/o checksum
+# HEADER_BYTES scalar ?
+# FIXED_LEADER_BYTES scalar ?
+# VARIABLE_LEADER_BYTES scalar ?
+# VELOCITY_DATA_BYTES scalar ?
+# CORRELATION_DATA_BYTES scalar ?
+# ECHO_INTENSITY_DATA_BYTES scalar ?
+# PERCENT_GOOD_DATA_BYTES scalar ?
+# BT_PRESENT bool NUMBER_OF_DATA_TYPES == 7
+# BT_DATA_BYTES scalar undefined, ? if BT_PRESENT
+# CPU_FW_VER scalar 0--255
+# CPU_FW_REV scalar 0--255
+# BEAM_FREQUENCY scalar 75, 150, 300, 600, 1200, 2400 [kHz]
+# CONVEX_BEAM_PATTERN bool undefined, 1
+# CONCAVE_BEAM_PATTERN bool undefined, 1
+# SENSOR_CONFIG scalar 1--3
+# XDUCER_HEAD_ATTACHED bool undefined, 1
+# BEAM_ANGLE scalar 15,20,30,undefined=other [deg]
+# N_BEAMS scalar 4--5
+# N_DEMODS scalar 2--3(???),undefined=n/a
+# N_BINS scalar 1--128
+# PINGS_PER_ENSEMBLE scalar 0--16384
+# BIN_LENGTH scalar 0.01--64 [m]
+# BLANKING_DISTANCE scalar 0-99.99 [m]
+# MIN_CORRELATION scalar 0--255
+# N_CODE_REPETITIONS scalar 0--255
+# MIN_PERCENT_GOOD scalar 1--100 [%]
+# MAX_ERROR_VELOCITY scalar 0--5 [m/s]
+# TIME_BETWEEN_PINGS scalar 0--? [s]
+# BEAM_COORDINATES bool undefined,1
+# INSTRUMENT_COORDINATES bool undefined,1
+# SHIP_COORDINATES bool undefined,1
+# EARTH_COORDINATES bool undefined,1
+# PITCH_AND_ROLL_USED bool undefined,1
+# USE_3_BEAM_ON_LOW_CORR bool undefined,1
+# BIN_MAPPING_ALLOWED bool undefined,1
+# HEADING_ALIGNMENT scalar -179.99..180 [deg]
+# HEADING_BIAS scalar -179.99..180 [deg]
+# CALCULATE_SPEED_OF_SOUND bool undefined,1
+# USE_PRESSURE_SENSOR bool undefined,1
+# USE_COMPASS bool undefined,1
+# USE_PITCH_SENSOR bool undefined,1
+# USE_ROLL_SENSOR bool undefined,1
+# USE_CONDUCTIVITY_SENSOR bool undefined,1
+# USE_TEMPERATURE_SENSOR bool undefined,1
+# SPEED_OF_SOUND_CALCULATED bool undefined,1
+# PRESSURE_SENSOR_AVAILABLE bool undefined,1
+# COMPASS_AVAILABLE bool undefined,1
+# PITCH_SENSOR_AVAILABLE bool undefined,1
+# ROLL_SENSOR_AVAILABLE bool undefined,1
+# CONDUCTIVITY_SENSOR_AVAILABLE bool undefined,1
+# TEMPERATURE_SENSOR_AVAILABLE bool undefined,1
+# DISTANCE_TO_BIN1_CENTER scalar 0--655.35 [m]
+# TRANSMITTED_PULSE_LENGTH scalar 0--655.35 [m]
+# RL_FIRST_BIN scalar 1--128
+# RL_LAST_BIN scalar 1--128
+# FALSE_TARGET_THRESHOLD scalar 0--254, undefined=disabled
+# LOW_LATENCY_SETTING scalar 0--5(???)
+# TRANSMIT_LAG_DISTANCE scalar 0--655.35 [m]
+# CPU_SERIAL_NUMBER scalar undefined, 0--65535 if WH300
+# NARROW_BANDWIDTH bool undefined,1 (only set if WH300)
+# WIDE_BANDWIDTH bool undefined,1 (only set if WH300)
+# TRANSMIT_POWER scalar undefined, 0--255(high) if WH300
+# TRANSMIT_POWER_HIGH bool undefined,1 (only set if WH300)
+# BT_PINGS_PER_ENSEMBLE scalar 0--999
+# BT_DELAY_BEFORE_REACQUIRE scalar 0--999
+# BT_MIN_CORRELATION scalar 0--255
+# BT_MIN_EVAL_AMPLITUDE scalar 0--255
+# BT_MIN_PERCENT_GOOD scalar 0--100 [%]
+# BT_MODE scalar 4,5,6(?)
+# BT_MAX_ERROR_VELOCITY scalar 0--5 [m/s], undef=not screened
+# BT_RL_MIN_SIZE scalar 0--99.9 [m]
+# BT_RL_NEAR scalar 0--999.9 [m]
+# BT_RL_FAR scalar 0--999.9 [m]
+# BT_MAX_TRACKING_DEPTH scalar 8--999.9 [m]
+# ENSEMBLE[ensemble_no-1] array ensemble info
+# XDUCER_FACING_UP bool undefined, 1
+# XDUCER_FACING_DOWN bool undefined, 1
+# N_BEAMS_USED scalar 3,4,5(?)
+# NUMBER scalar 1--16777215
+# BUILT_IN_TEST_ERROR scalar ?,undefined=none
+# SPEED_OF_SOUND scalar 1400--1600 [m/s]
+# XDUCER_DEPTH scalar 0.1--999.9 [m]
+# HEADING scalar 0--359.99 [deg]
+# PITCH scalar -20.00-20.00 [deg]
+# ROLL scalar -20.00-20.00 [deg]
+# SALINITY scalar 0-40 [psu]
+# TEMPERATURE scalar -5.00--40.00 [deg]
+# MIN_PRE_PING_WAIT_TIME scalar ? [s]
+# HEADING_STDDEV scalar 0-180 [deg]
+# PITCH_STDDEV scalar 0.0-20.0 [deg]
+# ROLL_STDDEV scalar 0.0-20.0 [deg]
+# ADC_XMIT_CURRENT scalar 0--255
+# ADC_XMIT_VOLTAGE scalar 0--255
+# ADC_AMBIENT_TEMPERATURE scalar 0--255
+# ADC_PRESSURE_PLUS scalar 0--255
+# ADC_PRESSURE_MINUS scalar 0--255
+# ADC_ATTITUDE_TEMPERATURE scalar 0--255
+# ADC_ATTITUDE scalar 0--255
+# ADC_CONTAMINATION scalar 0--255
+# ERROR_STATUS_WORD scalar undefined, ? (only set if WH300)
+# PRESSURE scalar undefined, ?-? [dbar] (only set if WH300)
+# PRESSURE_STDDEV scalar undefined, ?-? [dbar] (only set if WH300)
+# DATE string MM/DD/YYYY
+# YEAR scalar ?
+# MONTH scalar 1--12
+# DAY scalar 1--31
+# TIME string HH:MM:SS.hh
+# HOUR scalar 0--23
+# MINUTE scalar 0--59
+# SECONDS scalar 0--59.99
+# UNIX_TIME scalar 0--?
+# SECNO scalar 0--? (number of seconds since daystart)
+# DAYNO double fractional day number since start of current year (1.0 is midnight Jan 1st)
+# VELOCITY[bin][beam] scalars -32.767--32.768 [m/s], undef=bad
+# CORRELATION[bin][beam] scalars 1--255, undefined=bad
+# ECHO_AMPLITUDE[bin][beam] scalars 0--255
+# PERCENT_GOOD[bin][beam] scalars 0--255
+# BT_RANGE[beam] scalars tons [m]
+# BT_VELOCITY[beam] scalars see VELOCITY
+# BT_CORRELATION[beam] scalars see CORRELATION
+# BT_EVAL_AMPLITUDE[beam] scalars 0--255
+# BT_PERCENT_GOOD[beam] scalars see PERCENT_GOOD
+# BT_RL_VELOCITY[beam] scalars see VELOCITY
+# BT_RL_CORRELATION[beam] scalars see CORRELATION
+# BT_RL_ECHO_AMPLITUDE[beam] scalars see ECHO_AMPLITUDE
+# BT_RL_PERCENT_GOOD[beam] scalars see PERCENT_GOOD
+# BT_SIGNAL_STRENGTH[beam] scalars 0--255
+# HIGH_GAIN bool 1, undefined
+# LOW_GAIN bool 1, undefined
+
+use strict;
+use Time::Local; # timegm()
+
+#----------------------------------------------------------------------
+# Time Conversion Subroutines
+#----------------------------------------------------------------------
+
+sub monthLength($$) # of days in month
+{
+ my($Y,$M) = @_;
+
+ return 31 if ($M==1 || $M==3 || $M==5 || $M==7 ||
+ $M==8 || $M==10 || $M==12);
+ return 30 if ($M==4 || $M==6 || $M==9 || $M==11);
+ return 28 if ($Y%4 != 0);
+ return 29 if ($Y%100 != 0);
+ return 28 if ($Y%400 > 0);
+ return 29;
+}
+
+{ my($epoch,$lM,$lD,$lY,$ldn); # static scope
+
+ sub dayNo($$$$$$)
+ {
+ my($Y,$M,$D,$h,$m,$s) = @_;
+ my($dn);
+
+ if ($Y==$lY && $M==$lM && $D==$lD) { # same day as last samp
+ $dn = $ldn;
+ } else { # new day
+ $epoch = $Y unless defined($epoch); # 1st call
+ $lY = $Y; $lM = $M; $lD = $D; # store
+
+ for ($dn=0,my($cY)=$epoch; $cY<$Y; $cY++) { # multiple years
+ $dn += 337 + &monthLength($Y,$M);
+ }
+
+ $dn += $D; # day in month
+ while (--$M > 0) { # preceding months
+ $dn += &monthLength($Y,$M);
+ }
+
+ $ldn = $dn; # store
+ }
+ return $dn + $h/24 + $m/24/60 + $s/24/3600;
+ }
+
+} # static scope
+
+#----------------------------------------------------------------------
+# Read Data
+#----------------------------------------------------------------------
+
+my($WBRcfn); # current file name
+my(@WBRofs); # data type offsets
+
+my($FmtErr) = "%s: illegal %s Id 0x%04x at ensemble %d";
+
+sub WBRhdr($)
+{
+ my($dta) = @_;
+ my($buf,$hid,$did,$Ndt,$B,$W,$i,$dummy,$id);
+ my($B1,$B2,$B3,$B4,$B5,$B6,$B7,$W1,$W2,$W3,$W4,$W5);
+
+ #--------------------
+ # HEADER
+ #--------------------
+
+ read(WBRF,$buf,6) == 6 || die("$WBRcfn: $!\n");
+ ($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,"Data Source",$did,0));
+ die(sprintf("\n$WBRcfn: WARNING: unexpected number of data types (%d)\n",
+ $dta->{NUMBER_OF_DATA_TYPES}))
+ unless ($dta->{NUMBER_OF_DATA_TYPES} == 6 ||
+ $dta->{NUMBER_OF_DATA_TYPES} == 7);
+ $dta->{BT_PRESENT} = ($dta->{NUMBER_OF_DATA_TYPES} == 7);
+
+ read(WBRF,$buf,2*$dta->{NUMBER_OF_DATA_TYPES})
+ == 2*$dta->{NUMBER_OF_DATA_TYPES}
+ || die("$WBRcfn: $!\n");
+ @WBRofs = unpack("v$dta->{NUMBER_OF_DATA_TYPES}",$buf);
+
+ $dta->{HEADER_BYTES} = $WBRofs[0];
+ $dta->{FIXED_LEADER_BYTES} = $WBRofs[1] - $WBRofs[0];
+ $dta->{VARIABLE_LEADER_BYTES} = $WBRofs[2] - $WBRofs[1];
+ $dta->{VELOCITY_DATA_BYTES} = $WBRofs[3] - $WBRofs[2];
+ $dta->{CORRELATION_DATA_BYTES} = $WBRofs[4] - $WBRofs[3];
+ $dta->{ECHO_INTENSITY_DATA_BYTES} = $WBRofs[5] - $WBRofs[4];
+ if ($dta->{BT_PRESENT}) {
+ $dta->{PERCENT_GOOD_DATA_BYTES} = $WBRofs[6] - $WBRofs[5];
+ $dta->{BT_DATA_BYTES} = $dta->{ENSEMBLE_BYTES} - 4 - $WBRofs[6];
+ } else {
+ $dta->{PERCENT_GOOD_DATA_BYTES} = $dta->{ENSEMBLE_BYTES} - 4 - $WBRofs[5];
+ }
+
+# for ($i=0; $i<$dta->{NUMBER_OF_DATA_TYPES}; $i++) {
+# printf(STDERR "\nWBRofs[$i] = %d",$WBRofs[$i]);
+# }
+
+ $dta->{DATA_FORMAT_VARIANT} = 1;
+ if ($dta->{FIXED_LEADER_BYTES} == 53 || $dta->{FIXED_LEADER_BYTES} == 59) {
+ $dta->{DATA_FORMAT} = 'WH300';
+ $dta->{DATA_FORMAT_VARIANT} = 2 if ($dta->{FIXED_LEADER_BYTES} == 59);
+ } elsif ($dta->{FIXED_LEADER_BYTES} == 42) {
+ $dta->{DATA_FORMAT} = 'BB150';
+ } else {
+ printf(STDERR "\n$WBRcfn: WARNING: unknown data format (%d FIXED_LEADER_BYTES)\n",
+ $dta->{FIXED_LEADER_BYTES}
+ );
+ $dta->{DATA_FORMAT} = 'unknown';
+ }
+
+ #----------------------------------
+ # Check Data Format of 1st Ensemble
+ #----------------------------------
+
+ seek(WBRF,$WBRofs[1],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $id = unpack('v',$buf);
+ $id == 0x0080 || die(sprintf($FmtErr,$WBRcfn,"Variable Leader",$id,1));
+
+ seek(WBRF,$WBRofs[2],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $id = unpack('v',$buf);
+ $id == 0x0100 || die(sprintf($FmtErr,$WBRcfn,"Velocity Data",$id,1));
+
+ seek(WBRF,$WBRofs[3],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $id = unpack('v',$buf);
+ $id == 0x0200 || die(sprintf($FmtErr,$WBRcfn,"Correlation Data",$id,1));
+
+ seek(WBRF,$WBRofs[4],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $id = unpack('v',$buf);
+ $id == 0x0300 || die(sprintf($FmtErr,$WBRcfn,"Echo Intensity",$id,1));
+
+ seek(WBRF,$WBRofs[5],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $id = unpack('v',$buf);
+ $id == 0x0400 || die(sprintf($FmtErr,$WBRcfn,"Percent-Good Data",$id,1));
+
+ if ($dta->{BT_PRESENT}) {
+ seek(WBRF,$WBRofs[6],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $id = unpack('v',$buf);
+ $id == 0x0600 || die(sprintf($FmtErr,$WBRcfn,"Bottom Track",$id,1));
+ }
+
+ #--------------------
+ # FIXED LEADER
+ #--------------------
+
+ seek(WBRF,$WBRofs[0],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,42) == 42 || die("$WBRcfn: $!\n");
+ ($id,$dta->{CPU_FW_VER},$dta->{CPU_FW_REV},$B1,$B2,$dummy,$dummy,$dummy,
+ $dta->{N_BINS},$dta->{PINGS_PER_ENSEMBLE},$dta->{BIN_LENGTH},
+ $dta->{BLANKING_DISTANCE},$dummy,$dta->{MIN_CORRELATION},
+ $dta->{N_CODE_REPETITIONS},$dta->{MIN_PERCENT_GOOD},
+ $dta->{MAX_ERROR_VELOCITY},$dta->{TIME_BETWEEN_PINGS},$B3,$B4,$B5,
+ $dta->{HEADING_ALIGNMENT},$dta->{HEADING_BIAS},$B6,$B7,
+ $dta->{DISTANCE_TO_BIN1_CENTER},$dta->{TRANSMITTED_PULSE_LENGTH},
+ $dta->{REF_LAYER_FIRST_BIN},$dta->{REF_LAYER_LAST_BIN},
+ $dta->{FALSE_TARGET_THRESHOLD},$dta->{LOW_LATENCY_SETTING},
+ $dta->{TRANSMIT_LAG_DISTANCE}) =
+ unpack('vCCCCC3CvvvCCCCvCCCCvvCCvvCCCCv',$buf);
+
+ $id == 0x0000 || die(sprintf($FmtErr,$WBRcfn,"Fixed Leader",$id,0));
+
+ $dta->{BEAM_FREQUENCY} = 2**($B1 & 0x07) * 75;
+ $dta->{CONVEX_BEAM_PATTERN} = 1 if ($B1 & 0x08);
+ $dta->{CONCAVE_BEAM_PATTERN} = 1 if (!($B1 & 0x08));
+ $dta->{SENSOR_CONFIG} = ($B1 & 0x30) >> 4;
+ $dta->{XDUCER_HEAD_ATTACHED} = 1 if ($B1 & 0x40);
+
+ if (($B2 & 0x03) == 0x00) { $dta->{BEAM_ANGLE} = 15; }
+ elsif (($B2 & 0x03) == 0x01) { $dta->{BEAM_ANGLE} = 20; }
+ elsif (($B2 & 0x03) == 0x02) { $dta->{BEAM_ANGLE} = 30; }
+ if (($B2 & 0xF0) == 0x40) { $dta->{N_BEAMS} = 4; }
+ elsif (($B2 & 0xF0) == 0x50) { $dta->{N_BEAMS} = 5; $dta->{N_DEMODS} = 3; }
+ elsif (($B2 & 0xF0) == 0xF0) { $dta->{N_BEAMS} = 5; $dta->{N_DEMODS} = 2; }
+
+ $dta->{BIN_LENGTH} /= 100;
+ $dta->{BLANKING_DISTANCE} /= 100;
+
+ $dta->{MAX_ERROR_VELOCITY} /= 1000;
+ $dta->{TIME_BETWEEN_PINGS} *= 60;
+ $dta->{TIME_BETWEEN_PINGS} += $B3 + $B4/100;
+
+ $dta->{BEAM_COORDINATES} = 1 if (($B5 & 0x18) == 0x00);
+ $dta->{INSTRUMENT_COORDINATES} = 1 if (($B5 & 0x18) == 0x08);
+ $dta->{SHIP_COORDINATES} = 1 if (($B5 & 0x18) == 0x10);
+ $dta->{EARTH_COORDINATES} = 1 if (($B5 & 0x18) == 0x18);
+ $dta->{PITCH_AND_ROLL_USED} = 1 if ($B5 & 0x04);
+ $dta->{USE_3_BEAM_ON_LOW_CORR} = 1 if ($B5 & 0x02);
+ $dta->{BIN_MAPPING_ALLOWED} = 1 if ($B5 & 0x01);
+
+ $dta->{HEADING_ALIGNMENT} =
+ ($dta->{EARTH_COORDINATES} || $dta->{SHIP_COORDINATES}) ?
+ $dta->{HEADING_ALIGNMENT} / 100 : undef;
+ $dta->{HEADING_BIAS} =
+ ($dta->{EARTH_COORDINATES} || $dta->{SHIP_COORDINATES}) ?
+ $dta->{HEADING_BIAS} / 100 : undef;
+
+ $dta->{CALCULATE_SPEED_OF_SOUND} = 1 if ($B6 & 0x40);
+ $dta->{USE_PRESSURE_SENSOR} = 1 if ($B6 & 0x20);
+ $dta->{USE_COMPASS} = 1 if ($B6 & 0x10);
+ $dta->{USE_PITCH_SENSOR} = 1 if ($B6 & 0x08);
+ $dta->{USE_ROLL_SENSOR} = 1 if ($B6 & 0x04);
+ $dta->{USE_CONDUCTIVITY_SENSOR} = 1 if ($B6 & 0x02);
+ $dta->{USE_TEMPERATURE_SENSOR} = 1 if ($B6 & 0x01);
+
+ $dta->{SPEED_OF_SOUND_CALCULATED} = 1 if ($B7 & 0x40);
+ $dta->{PRESSURE_SENSOR_AVAILABLE} = 1 if ($B7 & 0x20);
+ $dta->{COMPASS_AVAILABLE} = 1 if ($B7 & 0x10);
+ $dta->{PITCH_SENSOR_AVAILABLE} = 1 if ($B7 & 0x08);
+ $dta->{ROLL_SENSOR_AVAILABLE} = 1 if ($B7 & 0x04);
+ $dta->{CONDUCTIVITY_SENSOR_AVAILABLE} = 1 if ($B7 & 0x02);
+ $dta->{TEMPERATURE_SENSOR_AVAILABLE} = 1 if ($B7 & 0x01);
+
+ $dta->{DISTANCE_TO_BIN1_CENTER} /= 100;
+ $dta->{TRANSMITTED_PULSE_LENGTH} /= 100;
+
+ $dta->{FALSE_TARGET_THRESHOLD} = undef
+ if ($dta->{FALSE_TARGET_THRESHOLD} == 255);
+ $dta->{TRANSMIT_LAG_DISTANCE} /= 100;
+
+ if ($dta->{DATA_FORMAT} eq 'WH300') {
+ read(WBRF,$buf,11) == 11 || die("$WBRcfn: $!\n");
+ ($W1,$W2,$W3,$W4,$W5,$dta->{TRANSMIT_POWER}) =
+ unpack('vvvvvC',$buf);
+
+ $dta->{CPU_SERIAL_NUMBER} = sprintf("%04X%04X%04X%04X",$W1,$W2,$W3,$W4);
+
+ $dta->{NARROW_BANDWIDTH} = ($W5 == 1);
+ $dta->{WIDE_BANDWIDTH} = ($W5 == 0);
+ $dta->{TRANSMIT_POWER_HIGH} = ($dta->{TRANSMIT_POWER} == 255);
+ }
+
+ #-----------------------
+ # 1st ENSEMBLE, BT Setup
+ #-----------------------
+
+ if ($dta->{BT_PRESENT}) {
+ seek(WBRF,$WBRofs[6],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,12) == 12 || die("$WBRcfn: $!\n");
+ ($id,$dta->{BT_PINGS_PER_ENSEMBLE},$dta->{BT_DELAY_BEFORE_REACQUIRE},
+ $dta->{BT_MIN_CORRELATION},$dta->{BT_MIN_EVAL_AMPLITUDE},
+ $dta->{BT_MIN_PERCENT_GOOD},$dta->{BT_MODE},
+ $dta->{BT_MAX_ERROR_VELOCITY}) = unpack('vvvCCCCv',$buf);
+
+ $id == 0x0600 ||
+ die(sprintf($FmtErr,$WBRcfn,"Bottom Track",$id,0,tell(WBRF)));
+
+ $dta->{BT_MAX_ERROR_VELOCITY} =
+ $dta->{BT_MAX_ERROR_VELOCITY} ? $dta->{BT_MAX_ERROR_VELOCITY} / 1000
+ : undef;
+
+ seek(WBRF,28,1) || die("$WBRcfn: $!");
+ read(WBRF,$buf,6) == 6 || die("$WBRcfn: $!\n");
+ ($dta->{BT_RL_MIN_SIZE},$dta->{BT_RL_NEAR},$dta->{BT_RL_FAR})
+ = unpack('vvv',$buf);
+
+ $dta->{BT_RL_MIN_SIZE} /= 10;
+ $dta->{BT_RL_NEAR} /= 10;
+ $dta->{BT_RL_FAR} /= 10;
+
+ seek(WBRF,20,1) || die("$WBRcfn: $!"); # skip data
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $dta->{BT_MAX_TRACKING_DEPTH} = unpack('v',$buf) / 10;
+ }
+
+ return $dta;
+}
+
+sub WBRens($$$$$)
+{
+ my($nbins,$ens_length,$BT_present,$data_format,$E) = @_;
+ my($start_ens,$B1,$B2,$B3,$B4,$I,$id,$bin,$beam,$buf,$dummy,@dta,$i,$cs);
+ my($ens,$ensNo,$dayStart);
+
+ for ($ens=$start_ens=0; 1; $ens++,$start_ens+=$ens_length+2) {
+# print(STDERR "start_ens = $start_ens\n");
+
+ #-------------------------------
+ # Make Sure Ensemble is Complete
+ #-------------------------------
+
+ # UH BB150 writes incomplete ensembles (i.e. short read
+ # indicates EOF). FSU WH300 has bogus data in incomplete
+ # final ensemble.
+
+ seek(WBRF,$start_ens,0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,$ens_length) == $ens_length || last;
+
+ read(WBRF,$cs,2) == 2 || last;
+ last unless (unpack('%16C*',$buf) == unpack('v',$cs));
+
+ #------------------------------
+ # Variable Leader
+ #------------------------------
+
+ seek(WBRF,$start_ens+$WBRofs[1],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,4) == 4 || die("$WBRcfn: $!\n");
+ ($id,$ensNo) = unpack("vv",$buf);
+
+ $id == 0x0080 ||
+ die(sprintf($FmtErr,$WBRcfn,"Variable Leader",$id,$ensNo+1));
+
+ if ($data_format eq 'BB150') { # non Y2K RTC
+ read(WBRF,$buf,7) == 7 || die("$WBRcfn: $!\n");
+ (${$E}[$ens]->{YEAR},${$E}[$ens]->{MONTH},
+ ${$E}[$ens]->{DAY},${$E}[$ens]->{HOUR},${$E}[$ens]->{MINUTE},
+ ${$E}[$ens]->{SECONDS},$B4) = unpack('CCCCCCC',$buf);
+ ${$E}[$ens]->{SECONDS} += $B4/100;
+ ${$E}[$ens]->{YEAR} += (${$E}[$ens]->{YEAR} > 80) ? 1900 : 2000;
+ } else {
+ seek(WBRF,7,1) || die("$WBRcfn: $!"); # use Y2K RTC
+ }
+
+ read(WBRF,$buf,1) == 1 || die("$WBRcfn: $!\n");
+ $ensNo += unpack('C',$buf) << 16;
+ ${$E}[$ens]->{NUMBER} = $ensNo;
+
+ read(WBRF,$buf,30) == 30 || die("$WBRcfn: $!\n");
+ (${$E}[$ens]->{BUILT_IN_TEST_ERROR},${$E}[$ens]->{SPEED_OF_SOUND},
+ ${$E}[$ens]->{XDUCER_DEPTH},${$E}[$ens]->{HEADING},
+ ${$E}[$ens]->{PITCH},${$E}[$ens]->{ROLL},
+ ${$E}[$ens]->{SALINITY},${$E}[$ens]->{TEMPERATURE},
+ ${$E}[$ens]->{MIN_PRE_PING_WAIT_TIME},$B1,$B2,
+ ${$E}[$ens]->{HEADING_STDDEV},${$E}[$ens]->{PITCH_STDDEV},
+ ${$E}[$ens]->{ROLL_STDDEV},${$E}[$ens]->{ADC_XMIT_CURRENT},
+ ${$E}[$ens]->{ADC_XMIT_VOLTAGE},${$E}[$ens]->{ADC_AMBIENT_TEMPERATURE},
+ ${$E}[$ens]->{ADC_PRESSURE_PLUS},${$E}[$ens]->{ADC_PRESSURE_MINUS},
+ ${$E}[$ens]->{ADC_ATTITUDE_TEMPERATURE},${$E}[$ens]->{ADC_ATTITUDE},
+ ${$E}[$ens]->{ADC_CONTAMINATION})
+ = unpack('vvvvvvvvCCCCCCCCCCCCCC',$buf);
+
+ ${$E}[$ens]->{BUILT_IN_TEST_ERROR} = undef
+ unless (${$E}[$ens]->{BUILT_IN_TEST_ERROR});
+ ${$E}[$ens]->{XDUCER_DEPTH} /= 10;
+ ${$E}[$ens]->{HEADING} /= 100;
+ ${$E}[$ens]->{PITCH} = unpack('s',pack('S',${$E}[$ens]->{PITCH})) / 100;
+ ${$E}[$ens]->{ROLL} = unpack('s',pack('S',${$E}[$ens]->{ROLL})) / 100;
+ ${$E}[$ens]->{TEMPERATURE} =
+ unpack('s',pack('S',${$E}[$ens]->{TEMPERATURE})) / 100;
+ ${$E}[$ens]->{MIN_PRE_PING_WAIT_TIME} *= 60;
+ ${$E}[$ens]->{MIN_PRE_PING_WAIT_TIME} += $B1 + $B2/100;
+ ${$E}[$ens]->{PITCH_STDDEV} /= 10;
+ ${$E}[$ens]->{ROLL_STDDEV} /= 10;
+
+ if ($data_format eq 'WH300') {
+ read(WBRF,$buf,23) == 23 || die("$WBRcfn: $!\n");
+ (${$E}[$ens]->{ERROR_STATUS_WORD},
+ $dummy,${$E}[$ens]->{PRESSURE},${$E}[$ens]->{PRESSURE_STDDEV},
+ $dummy,${$E}[$ens]->{YEAR},$B3,${$E}[$ens]->{MONTH},
+ ${$E}[$ens]->{DAY},${$E}[$ens]->{HOUR},${$E}[$ens]->{MINUTE},
+ ${$E}[$ens]->{SECONDS},$B4)
+ = unpack('VvVVCCCCCCCCC',$buf);
+
+ ${$E}[$ens]->{PRESSURE} /= 1000;
+ ${$E}[$ens]->{PRESSURE_STDDEV} /= 1000;
+ ${$E}[$ens]->{YEAR} *= 100; ${$E}[$ens]->{YEAR} += $B3;
+ ${$E}[$ens]->{SECONDS} += $B4/100;
+ }
+
+ ${$E}[$ens]->{DATE}
+ = sprintf("%02d/%02d/%d",${$E}[$ens]->{MONTH},
+ ${$E}[$ens]->{DAY},
+ ${$E}[$ens]->{YEAR});
+ ${$E}[$ens]->{TIME}
+ = sprintf("%02d:%02d:%05.02f",${$E}[$ens]->{HOUR},
+ ${$E}[$ens]->{MINUTE},
+ ${$E}[$ens]->{SECONDS});
+ ${$E}[$ens]->{DAYNO}
+ = &dayNo(${$E}[$ens]->{YEAR},${$E}[$ens]->{MONTH},${$E}[$ens]->{DAY},
+ ${$E}[$ens]->{HOUR},${$E}[$ens]->{MINUTE},${$E}[$ens]->{SECONDS});
+
+ ${$E}[$ens]->{UNIX_TIME}
+ = timegm(0,${$E}[$ens]->{MINUTE},
+ ${$E}[$ens]->{HOUR},
+ ${$E}[$ens]->{DAY},
+ ${$E}[$ens]->{MONTH}-1, # NB!!!
+ ${$E}[$ens]->{YEAR})
+ + ${$E}[$ens]->{SECONDS};
+
+ $dayStart = timegm(0,0,0,${$E}[$ens]->{DAY},
+ ${$E}[$ens]->{MONTH}-1, # NB!!!
+ ${$E}[$ens]->{YEAR})
+ unless defined($dayStart);
+ ${$E}[$ens]->{SECNO} = ${$E}[$ens]->{UNIX_TIME} - $dayStart;
+
+ seek(WBRF,$start_ens+$WBRofs[0]+4,0) # System Config / Fixed Leader
+ || die("$WBRcfn: $!");
+
+ read(WBRF,$buf,5) == 5 || die("$WBRcfn: $!\n");
+ ($B1,$dummy,$dummy,$dummy,${$E}[$ens]->{N_BEAMS_USED})
+ = unpack('CCCCC',$buf);
+ ${$E}[$ens]->{XDUCER_FACING_UP} = 1 if ($B1 & 0x80);
+ ${$E}[$ens]->{XDUCER_FACING_DOWN} = 1 unless ($B1 & 0x80);
+
+ #--------------------
+ # Velocity Data
+ #--------------------
+
+ my($ndata) = $nbins * 4;
+
+ seek(WBRF,$start_ens+$WBRofs[2],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2+$ndata*2) == 2+$ndata*2 || die("$WBRcfn: $!\n");
+ ($id,@dta) = unpack("vv$ndata",$buf);
+
+ $id == 0x0100 ||
+ die(sprintf($FmtErr,$WBRcfn,"Velocity Data",$id,$ens));
+
+ for ($i=0,$bin=0; $bin<$nbins; $bin++) {
+ for ($beam=0; $beam<4; $beam++,$i++) {
+ ${$E}[$ens]->{VELOCITY}[$bin][$beam] =
+ unpack('s',pack('S',$dta[$i])) / 1000
+ if ($dta[$i] != 0x8000);
+ }
+ }
+
+ #--------------------
+ # Correlation Data
+ #--------------------
+
+ seek(WBRF,$start_ens+$WBRofs[3],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2+$ndata) == 2+$ndata || die("$WBRcfn: $!\n");
+ ($id,@dta) = unpack("vC$ndata",$buf);
+
+ $id == 0x0200 ||
+ die(sprintf($FmtErr,$WBRcfn,"Correlation Data",$id,$ens));
+
+ for ($i=0,$bin=0; $bin<$nbins; $bin++) {
+ for ($beam=0; $beam<4; $beam++,$i++) {
+ ${$E}[$ens]->{CORRELATION}[$bin][$beam] = $dta[$i]
+ if ($dta[$i]);
+ }
+ }
+
+ #--------------------
+ # Echo Intensity Data
+ #--------------------
+
+ seek(WBRF,$start_ens+$WBRofs[4],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2+$ndata) == 2+$ndata || die("$WBRcfn: $!\n");
+ ($id,@dta) = unpack("vC$ndata",$buf);
+
+ $id == 0x0300 ||
+ die(sprintf($FmtErr,$WBRcfn,"Echo Intensity",$id,$ens));
+
+ for ($i=0,$bin=0; $bin<$nbins; $bin++) {
+ for ($beam=0; $beam<4; $beam++,$i++) {
+ ${$E}[$ens]->{ECHO_AMPLITUDE}[$bin][$beam] = $dta[$i];
+ }
+ }
+
+ #--------------------
+ # Percent Good Data
+ #--------------------
+
+ seek(WBRF,$start_ens+$WBRofs[5],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2+$ndata) == 2+$ndata || die("$WBRcfn: $!\n");
+ ($id,@dta) = unpack("vC$ndata",$buf);
+
+ $id == 0x0400 ||
+ die(sprintf($FmtErr,$WBRcfn,"Percent-Good Data",$id,$ens));
+
+ for ($i=0,$bin=0; $bin<$nbins; $bin++) {
+ for ($beam=0; $beam<4; $beam++,$i++) {
+ ${$E}[$ens]->{PERCENT_GOOD}[$bin][$beam] = $dta[$i];
+ }
+ }
+
+ #--------------------
+ # Bottom-Track Data
+ #--------------------
+
+ if ($BT_present) {
+ seek(WBRF,$start_ens+$WBRofs[6],0) || die("$WBRcfn: $!");
+ read(WBRF,$buf,2) == 2 || die("$WBRcfn: $!\n");
+ $id = unpack('v',$buf);
+
+ $id == 0x0600 ||
+ die(sprintf($FmtErr,$WBRcfn,"Bottom Track",$id,$ens));
+
+ seek(WBRF,14,1) || die("$WBRcfn: $!"); # BT config
+
+ read(WBRF,$buf,28) == 28 || die("$WBRcfn: $!\n");
+ @dta = unpack('v4v4C4C4C4',$buf);
+
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_RANGE}[$beam] = $dta[$beam] / 100
+ if ($dta[$beam]);
+ }
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_VELOCITY}[$beam] =
+ unpack('s',pack('S',$dta[4+$beam])) / 1000
+ if ($dta[4+$beam] != 0x8000);
+ }
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_CORRELATION}[$beam] = $dta[8+$beam]
+ if ($dta[8+$beam]);
+ }
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_EVAL_AMPLITUDE}[$beam] = $dta[12+$beam];
+ }
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_PERCENT_GOOD}[$beam] = $dta[16+$beam];
+ }
+
+ seek(WBRF,6,1) || die("$WBRcfn: $!"); # BT config
+
+ read(WBRF,$buf,20) == 20 || die("$WBRcfn: $!\n");
+ @dta = unpack('v4C4C4C4',$buf);
+
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_RL_VELOCITY}[$beam] =
+ unpack('s',pack('S',$dta[$beam])) / 1000
+ if ($dta[$beam] != 0x8000);
+ }
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_RL_CORRELATION}[$beam] = $dta[4+$beam]
+ if ($dta[4+$beam]);
+ }
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_RL_ECHO_AMPLITUDE}[$beam] = $dta[8+$beam];
+ }
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_RL_PERCENT_GOOD}[$beam] = $dta[12+$beam];
+ }
+
+ seek(WBRF,2,1) || die("$WBRcfn: $!"); # BT config
+
+ read(WBRF,$buf,9) == 9 || die("$WBRcfn: $!\n");
+ @dta = unpack('C4CC4',$buf);
+
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_SIGNAL_STRENGTH}[$beam] = $dta[$beam];
+ }
+ ${$E}[$ens]->{HIGH_GAIN} if ($dta[4]);
+ ${$E}[$ens]->{LOW_GAIN} unless ($dta[4]);
+ for ($beam=0; $beam<4; $beam++) {
+ ${$E}[$ens]->{BT_RANGE}[$beam] += $dta[5+$beam] * 655.36
+ if ($dta[5+$beam]);
+ }
+ } # BT present
+ } # ens loop
+}
+
+sub readHeader(@)
+{
+ my($WBRcfn,$dta) = @_;
+ open(WBRF,$WBRcfn) || die("$WBRcfn: $!\n");
+ WBRhdr($dta);
+}
+
+sub readData(@)
+{
+ my($WBRcfn,$dta) = @_;
+ open(WBRF,$WBRcfn) || die("$WBRcfn: $!\n");
+ WBRhdr($dta);
+ WBRens($dta->{N_BINS},$dta->{ENSEMBLE_BYTES},
+ $dta->{BT_PRESENT},$dta->{DATA_FORMAT},
+ \@{$dta->{ENSEMBLE}});
+}
+
+sub checkEnsemble($$)
+{
+ my($dta,$ens) = @_;
+ printf(STDERR "3 beams used in ensemble #$dta->{ENSEMBLE}[$ens]->{NUMBER}\n")
+ if ($dta->{ENSEMBLE}[$ens]->{N_BEAMS_USED} < 4);
+ die("BIT error in ensemble $dta->{ENSEMBLE}[$ens]->{NUMBER}\n")
+ if defined($dta->{ENSEMBLE}[$ens]->{BUILT_IN_TEST_ERROR});
+#
+# ERROR_STATUS_WORD CONTAINS APPARENTLY INSTRUMENT-SPECIFIC VALUES
+# => CHECK DISABLED
+#
+# die(sprintf("ESW = 0x%08lx in ensemble #%s\n",
+# $dta->{ENSEMBLE}[$ens]->{ERROR_STATUS_WORD},
+# $dta->{ENSEMBLE}[$ens]->{NUMBER}))
+# if ($dta->{ENSEMBLE}[$ens]->{ESW_ERROR});
+}
+
+1; # return true for all the world to see
new file mode 100644
--- /dev/null
+++ b/RDI_Coords.pl
@@ -0,0 +1,255 @@
+#======================================================================
+# R D I _ C O O R D S . P L
+# doc: Sun Jan 19 17:57:53 2003
+# dlm: Sun May 23 22:47:32 2010
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 28 74 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# RDI Workhorse Coordinate Transformations
+
+# HISTORY:
+# Jan 19, 2003: - written
+# Jan 21, 2003: - made it obey HEADING_BIAS (magnetic declination)
+# Jan 22, 3003: - corrected magnetic declination
+# Feb 16, 2003: - use pitch correction from RDI manual
+# Oct 11, 2003: - BUG: return value of atan() had been interpreted
+# as degrees instead of radians
+# Feb 27, 2004: - added velApplyHdgBias()
+# - changed non-zero HEADING_ALIGNMENT from error to warning
+# Sep 16, 2005: - added deg() for [mkprofile]
+# Aug 26, 2006: - BUG: incorrect transformation for uplookers
+# Nov 30, 2007: - optimized &velInstrumentToEarth(), velBeamToInstrument()
+# - added support for 3-beam solutions
+# Feb 12, 2008: - added threeBeamFlag
+# Mar 18, 2009: - added &gimbal_pitch(), &angle_from_vertical()
+# May 19, 2009: - added &velBeamToVertical()
+# May 23, 2009: - debugged & renamed to &velBeamToBPEarth
+# May 23, 2010: - changed prototypes of rad() & deg() to conform to ANTS
+
+use strict;
+use POSIX;
+
+my($PI) = 3.14159265358979;
+
+sub rad(@) { return $_[0]/180 * $PI; }
+sub deg(@) { return $_[0]/$PI * 180; }
+
+$RDI_Coords::minValidVels = 3; # 3-beam solutions ok
+
+$RDI_Coords::threeBeam_1 = 0; # stats
+$RDI_Coords::threeBeam_2 = 0;
+$RDI_Coords::threeBeam_3 = 0;
+$RDI_Coords::threeBeam_4 = 0;
+$RDI_Coords::fourBeam = 0;
+
+$RDI_Coords::threeBeamFlag = 0; # flag last transformation
+
+{ # STATIC SCOPE
+ my(@B2I);
+
+ sub velBeamToInstrument(@)
+ {
+ my($dta,$v1,$v2,$v3,$v4) = @_;
+ return undef unless (defined($v1) + defined($v2) +
+ defined($v3) + defined($v4)
+ >= $RDI_Coords::minValidVels);
+
+ unless (defined(@B2I)) {
+# print(STDERR "RDI_Coords::minValidVels = $RDI_Coords::minValidVels\n");
+ my($a) = 1 / (2 * sin(rad($dta->{BEAM_ANGLE})));
+ my($b) = 1 / (4 * cos(rad($dta->{BEAM_ANGLE})));
+ my($c) = $dta->{CONVEX_BEAM_PATTERN} ? 1 : -1;
+ my($d) = $a / sqrt(2);
+ @B2I = ([$c*$a, -$c*$a, 0, 0 ],
+ [0, 0, -$c*$a, $c*$a],
+ [$b, $b, $b, $b ],
+ [$d, $d, -$d, -$d ]);
+# print(STDERR "@{$B2I[0]}\n@{$B2I[1]}\n@{$B2I[2]}\n@{$B2I[3]}\n");
+ }
+
+ if (!defined($v1)) { # 3-beam solutions
+ $RDI_Coords::threeBeamFlag = 1;
+ $RDI_Coords::threeBeam_1++;
+ $v1 = -($v2*$B2I[3][1]+$v3*$B2I[3][2]+$v4*$B2I[3][3])/$B2I[3][0];
+ } elsif (!defined($v2)) {
+ $RDI_Coords::threeBeamFlag = 1;
+ $RDI_Coords::threeBeam_2++;
+ $v2 = -($v1*$B2I[3][0]+$v3*$B2I[3][2]+$v4*$B2I[3][3])/$B2I[3][1];
+ } elsif (!defined($v3)) {
+ $RDI_Coords::threeBeamFlag = 1;
+ $RDI_Coords::threeBeam_3++;
+ $v3 = -($v1*$B2I[3][0]+$v2*$B2I[3][1]+$v4*$B2I[3][3])/$B2I[3][2];
+ } elsif (!defined($v4)) {
+ $RDI_Coords::threeBeamFlag = 1;
+ $RDI_Coords::threeBeam_4++;
+ $v4 = -($v1*$B2I[3][0]+$v2*$B2I[3][1]+$v3*$B2I[3][2])/$B2I[3][3];
+ } else {
+ $RDI_Coords::threeBeamFlag = 0;
+ $RDI_Coords::fourBeam++;
+ }
+
+ return ($v1*$B2I[0][0]+$v2*$B2I[0][1],
+ $v3*$B2I[1][2]+$v4*$B2I[1][3],
+ $v1*$B2I[2][0]+$v2*$B2I[2][1]+$v3*$B2I[2][2]+$v4*$B2I[2][3],
+ $v1*$B2I[3][0]+$v2*$B2I[3][1]+$v3*$B2I[3][2]+$v4*$B2I[3][3]);
+ }
+} # STATIC SCOPE
+
+{ # STATIC SCOPE
+ my($hdg,$pitch,$roll,@I2E);
+
+ sub velInstrumentToEarth(@)
+ {
+ my($dta,$ens,$v1,$v2,$v3,$v4) = @_;
+ return undef unless (defined($v1) && defined($v2) &&
+ defined($v3) && defined($v4));
+
+ unless (@I2E &&
+ $hdg == $dta->{ENSEMBLE}[$ens]->{HEADING}
+ - $dta->{HEADING_BIAS} &&
+ $pitch == $dta->{ENSEMBLE}[$ens]->{PITCH} &&
+ $roll == $dta->{ENSEMBLE}[$ens]->{ROLL}) {
+ printf(STDERR "$0: warning HEADING_ALIGNMENT == %g ignored\n",
+ $dta->{HEADING_ALIGNMENT})
+ if ($dta->{HEADING_ALIGNMENT});
+ $hdg = $dta->{ENSEMBLE}[$ens]->{HEADING} - $dta->{HEADING_BIAS};
+ $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)));
+ my($sp,$cp) = (sin($rad_gimbal_pitch),cos($rad_gimbal_pitch));
+ my($sr,$cr) = (sin(rad($roll)), cos(rad($roll)));
+ @I2E = $dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}
+ ? (
+ [-$ch*$cr-$sh*$sp*$sr, $sh*$cp,-$ch*$sr+$sh*$sp*$cr],
+ [-$ch*$sp*$sr+$sh*$cr, $ch*$cp, $sh*$sr+$ch*$sp*$cr],
+ [+$cp*$sr, $sp, -$cp*$cr, ],
+ ) : (
+ [$ch*$cr+$sh*$sp*$sr, $sh*$cp, $ch*$sr-$sh*$sp*$cr],
+ [$ch*$sp*$sr-$sh*$cr, $ch*$cp,-$sh*$sr-$ch*$sp*$cr],
+ [-$cp*$sr, $sp, $cp*$cr, ],
+ );
+ }
+ return ($v1*$I2E[0][0]+$v2*$I2E[0][1]+$v3*$I2E[0][2],
+ $v1*$I2E[1][0]+$v2*$I2E[1][1]+$v3*$I2E[1][2],
+ $v1*$I2E[2][0]+$v2*$I2E[2][1]+$v3*$I2E[2][2],
+ $v4);
+
+ }
+} # STATIC SCOPE
+
+#======================================================================
+# velBeamToBPEarth3(@) 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.
+#======================================================================
+
+{ # STATIC SCOPE
+ my($TwoCosBAngle,$TwoSinBAngle);
+
+ sub velBeamToBPEarth(@)
+ {
+ my($dta,$ens,$b1,$b2,$b3,$b4) = @_;
+ my($v12,$w12,$v34,$w34);
+
+ return (undef,undef,undef,undef)
+ unless defined($b1) && defined($b2) && defined($b3) && defined($b4);
+
+ unless (defined($TwoCosBAngle)) {
+ $TwoCosBAngle = 2 * cos(rad($dta->{BEAM_ANGLE}));
+ $TwoSinBAngle = 2 * sin(rad($dta->{BEAM_ANGLE}));
+ }
+ my($roll) = rad($dta->{ENSEMBLE}[$ens]->{ROLL});
+ my($sr) = sin($roll); my($cr) = cos($roll);
+ my($pitch) = atan(tan(rad($dta->{ENSEMBLE}[$ens]->{PITCH})) * $cr); # gimbal pitch
+ 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($w12_ic) = ($b1+$b2)/$TwoCosBAngle;
+ $w12_ic *= -1 if ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP});
+ my($v34_ic) = ($b3-$b4)/$TwoSinBAngle;
+ my($w34_ic) = ($b3+$b4)/$TwoCosBAngle;
+ $w34_ic *= -1 if ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP});
+
+ if ($dta->{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;
+ }
+
+ return ($v12,$w12,$v34,$w34);
+ }
+}
+
+#======================================================================
+# velApplyHdgBias() applies the heading bias, which is used to correct
+# for magnetic declination for data recorded in Earth-coordinates ONLY.
+# Bias correction for beam-coordinate data is done in velInstrumentToEarth()
+#======================================================================
+
+{ # STATIC SCOPE
+ my($sh,$ch);
+
+ sub velApplyHdgBias(@)
+ {
+ my($dta,$ens,$v1,$v2,$v3,$v4) = @_;
+ return undef unless (defined($v1) && defined($v2));
+
+ unless (defined($sh)) {
+ printf(STDERR "$0: warning HEADING_ALIGNMENT == %g ignored\n",
+ $dta->{HEADING_ALIGNMENT})
+ if ($dta->{HEADING_ALIGNMENT});
+ $sh = sin(rad(-$dta->{HEADING_BIAS}));
+ $ch = cos(rad(-$dta->{HEADING_BIAS}));
+ }
+
+ return ( $v1*$ch + $v2*$sh,
+ -$v1*$sh + $v2*$ch,
+ $v3 ,
+ $v4 );
+ }
+} # STATIC SCOPE
+
+#----------------------------------------------------------------------
+# Pitch/Roll Functions
+#----------------------------------------------------------------------
+
+sub gimbal_pitch($$) # RDI coord trans manual
+{
+ my($tilt1,$tilt2) = @_;
+ return deg(atan(tan(rad($tilt1)) * cos(rad($tilt2))));
+}
+
+# - angle from vertical is home grown and should be treated with caution
+# - angle between two unit vectors given by acos(v1 dot v2)
+# - vertical unit vector v1 = (0 0 1) => dot product = z-component of v2
+# - when vertical unit vector is pitched in x direction, followed by
+# roll in y direction:
+# x = sin(pitch)
+# y = cos(pitch) * sin(roll)
+# z = cos(pitch) * cos(roll)
+# has been checked with sqrt(x^2+y^2+z^2) == 1
+# - for small angles, this is very similar to sqrt(pitch^2+roll^2)
+
+sub angle_from_vertical($$)
+{
+ my($tilt1,$tilt2) = @_;
+ my($rad_pitch) = atan(tan(rad($tilt1)) * cos(rad($tilt2)));
+ return deg(acos(cos($rad_pitch) * cos(rad($tilt2))));
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/RDI_Utils.pl
@@ -0,0 +1,360 @@
+#======================================================================
+# R D I _ U T I L S . P L
+# doc: Wed Feb 12 10:21:32 2003
+# dlm: Sun May 23 16:35:21 2010
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 156 42 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# miscellaneous RDI-specific utilities
+
+# History:
+# Feb 12, 2003: - created
+# Feb 14, 2003: - added check for short (bad) data files
+# Feb 26, 2004: - added Earth-coordinate support
+# - added ensure_BT_RANGE()
+# Mar 17, 2004: - set bad BT ranges to undef in ensure_BT_RANGE
+# - calc mean/stddev in ensure_BT_RANGE
+# Mar 20, 2004: - BUG: find_seabed() could bomb when not enough
+# bottom guesses passed the mode_width test
+# Mar 30, 2004: - added &soundSpeed()
+# Nov 8, 2005: - WATER_DEPTH => Z_BT
+# Dec 1, 2005: - renamed to RDI_Utils.pl
+# - folded in mk_prof from [mkProfile]
+# Nov 30, 2007: - adapted ref_lr_w() to 3-beam solutions
+# Feb 1, 2008: - added comment
+# Feb 22, 2008: - added ssCorr()
+# Apr 9, 2008: - BUG: duplicate line of code (without effect) in find_seabed()
+# - BUG: seabed < max depth was possible
+# Jan 2010: - fiddled with seabed detection params (no conclusion)
+# May 23, 2010: - renamed Z to DEPTH
+
+use strict;
+
+#======================================================================
+# fake_BT_RANGE(dta ptr)
+#======================================================================
+
+# During cruise NBP0204 it was found that one of Visbeck's RDIs consistently
+# returns zero as the bottom-track range, even thought the BT velocities
+# are good. This functions calculates the ranges if they are missing.
+
+sub cBTR($$$)
+{
+ my($d,$e,$b) = @_;
+ my($maxamp) = -9e99;
+ my($maxi);
+
+ for (my($i)=0; $i<$d->{N_BINS}; $i++) {
+ next unless ($d->{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b] > $maxamp);
+ $maxamp = $d->{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b];
+ $maxi = $i;
+ }
+ $d->{ENSEMBLE}[$e]->{BT_RANGE}[$b] =
+ $d->{DISTANCE_TO_BIN1_CENTER} + $maxi * $d->{BIN_LENGTH};
+}
+
+sub ensure_BT_RANGE($)
+{
+ my($d) = @_;
+ for (my($e)=0; $e<=$#{$d->{ENSEMBLE}}; $e++) {
+ my($sum) = my($n) = 0;
+ if (defined($d->{ENSEMBLE}[$e]->{BT_VELOCITY}[0])) {
+ for (my($b)=0; $b<4; $b++) {
+ cBTR($d,$e,$b)
+ unless defined($d->{ENSEMBLE}[$e]->{BT_RANGE}[$b]);
+ $sum += $d->{ENSEMBLE}[$e]->{BT_RANGE}[$b]; $n++;
+ }
+ } else {
+ for (my($b)=0; $b<4; $b++) {
+ $d->{ENSEMBLE}[$e]->{BT_RANGE}[$b] = undef;
+ }
+ }
+ $d->{ENSEMBLE}[$e]->{BT_MEAN_RANGE} = $sum/$n if ($n == 4);
+ }
+}
+
+#======================================================================
+# (seabed depth, stddev) = find_seabed(dta ptr, btm ensno, coord flag)
+#======================================================================
+
+# NOTE FOR YOYOS:
+# - this routine only detects the BT around the depeest depth!
+# - this is on purpose, because it is used for [listBT]
+
+# This is a PAIN:
+# if the instrument is too close to the bottom, the BT_RANGE
+# readings are all out; the only solution is to have a sufficiently
+# large search width (around the max(depth) ensemble) so that
+# a part of the up (oops, I forgot to finish this comment one year
+# ago during A0304 and now I don't understand it any more :-)
+
+my($search_width) = 200; # # of ensembles around bottom to search
+my($mode_width) = 10; # max range of bottom around mode
+my($min_dist) = 20; # min dist to seabed for good data
+my($z_offset) = 10000; # shift z to ensure +ve array indices
+
+sub find_seabed($$$)
+{
+ my($d,$be,$beamCoords) = @_;
+ my($i,$dd,$sd,$nd);
+ my(@guesses);
+
+ return undef unless ($be-$search_width >= 0 &&
+ $be+$search_width <= $#{$d->{ENSEMBLE}});
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless (defined($d->{ENSEMBLE}[$i]->{DEPTH}) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[0]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[1]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[2]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[3]));
+ my(@BT) = $beamCoords
+ ? velInstrumentToEarth($d,$i,
+ velBeamToInstrument($d,
+ @{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}}))
+ : velApplyHdgBias($d,$i,@{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}});
+ next unless (abs($BT[3]) < 0.05);
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} =
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[0]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[1]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[2]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[3]/4;
+ next unless ($d->{ENSEMBLE}[$i]->{DEPTH_BT} >= $min_dist);
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} *= -1
+ if ($d->{ENSEMBLE}[$i]->{XDUCER_FACING_UP});
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} += $d->{ENSEMBLE}[$i]->{DEPTH};
+ if ($d->{ENSEMBLE}[$i]->{DEPTH_BT} > $d->{ENSEMBLE}[$be]->{DEPTH}) {
+ $guesses[int($d->{ENSEMBLE}[$i]->{DEPTH_BT})+$z_offset]++;
+ $nd++;
+ } else {
+ undef($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ }
+ }
+ return undef unless ($nd>5);
+
+ my($mode,$nmax);
+ for ($i=0; $i<=$#guesses; $i++) { # find mode
+ $nmax=$guesses[$i],$mode=$i-$z_offset
+ if ($guesses[$i] > $nmax);
+ }
+
+ $nd = 0;
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ if (abs($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$mode) <= $mode_width) {
+ $dd += $d->{ENSEMBLE}[$i]->{DEPTH_BT};
+ $nd++;
+ } else {
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} = undef;
+ }
+ }
+ return undef unless ($nd >= 2);
+
+ $dd /= $nd;
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ $sd += ($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$dd)**2;
+ }
+
+ return ($dd, sqrt($sd/($nd-1)));
+}
+
+#======================================================================
+# c = soundSpeed($salin,$temp,$depth)
+#======================================================================
+
+# Taken from the RDI BroadBand primer manual. The reference given there
+# is Urick (1983).
+
+sub soundSpeed($$$)
+{
+ my($salin,$temp,$depth) = @_;
+ die("ERROR: soundSpeed($salin,$temp,$depth): non-numeric parameter\n")
+ unless numberp($salin) && numberp($temp) && numberp($depth);
+ return 1449.2 + 4.6*$temp -0.055*$temp**2 + 0.00029*$temp**3 +
+ (1.34 - 0.01*$temp) * ($salin - 35) + 0.016*$depth;
+}
+
+#======================================================================
+# fac = ssCorr($eRef,$salin,$temp,$depth)
+# $eRef : reference to current ensemble
+# $salin: * -> use instrument salinity
+# $temp : * -> use instrument temperature
+# $depth: * -> use instrument PRESSURE(!)
+#======================================================================
+
+{ my($warned);
+ sub ssCorr($$$$)
+ {
+ my($eRef,$S,$T,$D) = @_;
+ $S = $eRef->{SALINITY} if ($S eq '*');
+ $T = $eRef->{TEMPERATURE} if ($T eq '*');
+ if ($D eq '*') {
+ print(STDERR "WARNING: soundspeed correction using instrument pressure instead of depth!\n")
+ unless ($warned);
+ $warned = 1;
+ $D = $eRef->{PRESSURE};
+ }
+ return soundSpeed($S,$T,$D) / $eRef->{SPEED_OF_SOUND};
+ }
+}
+
+#======================================================================
+# ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz) =
+# mk_prof($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap);
+#======================================================================
+
+sub ref_lr_w($$$$$$) # calc ref-level vert vels
+{
+ my($dta,$ens,$rl_b0,$rl_b1,$min_corr,$max_e) = @_;
+ my($i,@n,@bn,@v,@vel,@bv,@w);
+
+ for ($i=$rl_b0; $i<=$rl_b1; $i++) {
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0])
+ if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $min_corr);
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1])
+ if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $min_corr);
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2])
+ if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $min_corr);
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3])
+ if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $min_corr);
+ if ($dta->{BEAM_COORDINATES}) {
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0])
+ if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100);
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1])
+ if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100);
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2])
+ if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100);
+ undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3])
+ if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+ @v = velInstrumentToEarth($dta,$ens,
+ velBeamToInstrument($dta,
+ @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
+ } else {
+ next if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||
+ $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||
+ $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||
+ $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+ @v = @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]};
+ # NB: no need to apply heading bias, as long as we only use w!
+ }
+ next if (!defined($v[3]) || abs($v[3]) > $max_e);
+
+ if (defined($v[2])) { # valid w
+ $vel[2] += $v[2]; $n[2]++;
+ $vel[3] += $v[3], $n[3]++ if defined($v[3]);
+ push(@w,$v[2]); # for stderr test
+ }
+
+ if ($dta->{BEAM_COORDINATES}) {
+ $bv[0] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0], $bn[0]++
+ if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0]);
+ $bv[1] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1], $bn[1]++
+ if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1]);
+ $bv[2] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2], $bn[2]++
+ if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2]);
+ $bv[3] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3], $bn[3]++
+ if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3]);
+ }
+ }
+
+ my($w) = $n[2] ? $vel[2]/$n[2] : undef; # w uncertainty
+ my($sumsq) = 0;
+ for ($i=0; $i<=$#w; $i++) {
+ $sumsq += ($w-$w[$i])**2;
+ }
+ my($stderr) = $n[2]>=2 ? sqrt($sumsq)/($n[2]-1) : undef;
+# The following stderr test introduces a huge gap near the bottom of
+# the profiles. Without it, they seem more reasonable.
+# next if ($stderr > 0.05);
+
+ if (defined($w)) { # valid w
+ $dta->{ENSEMBLE}[$ens]->{W} = $w;
+ $dta->{ENSEMBLE}[$ens]->{W_ERR} = $stderr;
+ }
+ if ($dta->{BEAM_COORDINATES}) {
+ $dta->{ENSEMBLE}[$ens]->{V1} = $bn[0]>=2 ? $bv[0]/$bn[0] : undef;
+ $dta->{ENSEMBLE}[$ens]->{V2} = $bn[1]>=2 ? $bv[1]/$bn[1] : undef;
+ $dta->{ENSEMBLE}[$ens]->{V3} = $bn[2]>=2 ? $bv[2]/$bn[2] : undef;
+ $dta->{ENSEMBLE}[$ens]->{V4} = $bn[3]>=2 ? $bv[3]/$bn[3] : undef;
+ }
+}
+
+
+sub mk_prof($$$$$$$$) # make profile
+{
+ my($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap) = @_;
+ my($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz);
+
+ for (my($z)=0,my($e)=0; $e<=$#{$dta->{ENSEMBLE}}; $e++) {
+ checkEnsemble($dta,$e) if ($check);
+
+ filterEnsemble($dta,$e)
+ if (defined($filter) &&
+ $dta->{ENSEMBLE}[$e]->{PERCENT_GOOD}[0][0] > 0);
+ ref_lr_w($dta,$e,$lr_b0,$lr_b1,$min_corr,$max_e); # ref. layer w
+
+ if (defined($firstgood)) {
+ $dta->{ENSEMBLE}[$e]->{ELAPSED_TIME} = # time since start
+ $dta->{ENSEMBLE}[$e]->{UNIX_TIME} -
+ $dta->{ENSEMBLE}[$firstgood]->{UNIX_TIME};
+ } else {
+ if (defined($dta->{ENSEMBLE}[$e]->{W})) { # start of prof.
+ $firstgood = $lastgood = $e;
+ $dta->{ENSEMBLE}[$e]->{ELAPSED_TIME} = 0;
+ $dta->{ENSEMBLE}[$e]->{DEPTH} = $dta->{ENSEMBLE}[$e]->{DEPTH_ERR} = 0;
+ }
+ next;
+ }
+
+ #--------------------------------------------------
+ # within profile: both $firstgood and $lastgood set
+ #--------------------------------------------------
+
+ if (!defined($dta->{ENSEMBLE}[$e]->{W})) { # gap
+ $w_gap_time += $dta->{ENSEMBLE}[$e]->{UNIX_TIME} -
+ $dta->{ENSEMBLE}[$e-1]->{UNIX_TIME};
+ next;
+ }
+
+ my($dt) = $dta->{ENSEMBLE}[$e]->{UNIX_TIME} - # time step since
+ $dta->{ENSEMBLE}[$lastgood]->{UNIX_TIME}; # ... last good ens
+
+ if ($dt > $max_gap) {
+ printf(STDERR "WARNING: %d-s gap too long, profile restarted at ensemble $e\n",$dt);
+ $firstgood = $lastgood = $e;
+ $dta->{ENSEMBLE}[$e]->{ELAPSED_TIME} = 0;
+ $z = $zErr = $maxz = 0;
+ $dta->{ENSEMBLE}[$e]->{DEPTH} = $dta->{ENSEMBLE}[$e]->{DEPTH_ERR} = 0;
+ $w_gap_time = 0;
+ next;
+ }
+
+ #-----------------------------------
+ # The current ensemble has a valid w
+ #-----------------------------------
+
+ $z += $dta->{ENSEMBLE}[$lastgood]->{W} * $dt; # integrate
+ $zErr += ($dta->{ENSEMBLE}[$lastgood]->{W_ERR} * $dt)**2;
+ $dta->{ENSEMBLE}[$e]->{DEPTH} = $z;
+ $dta->{ENSEMBLE}[$e]->{DEPTH_ERR} = sqrt($zErr);
+
+ $atbottom = $e, $maxz = $z if ($z > $maxz);
+
+ $lastgood = $e;
+ }
+
+ filterEnsembleStats() if defined($filter);
+
+ return ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz);
+}
+
+#----------------------------------------------------------------------
+# (true|false) = numberp(var)
+#----------------------------------------------------------------------
+
+sub numberp(@)
+{ return $_[0] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/; }
+
+
+1;
new file mode 100644
--- /dev/null
+++ b/WorkhorseUtils.pl
@@ -0,0 +1,162 @@
+#======================================================================
+# W O R K H O R S E U T I L S . P L
+# doc: Wed Feb 12 10:21:32 2003
+# dlm: Sun May 23 16:32:37 2010
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 142 42 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# miscellaneous Workhorse-specific utilities
+
+# History:
+# Feb 12, 2003: - created
+# Feb 14, 2003: - added check for short (bad) data files
+# Feb 26, 2004: - added Earth-coordinate support
+# - added ensure_BT_RANGE()
+# Mar 17, 2004: - set bad BT ranges to undef in ensure_BT_RANGE
+# - calc mean/stddev in ensure_BT_RANGE
+# Mar 20, 2004: - BUG: find_seabed() could bomb when not enough
+# bottom guesses passed the mode_width test
+# Mar 30, 2004: - added &soundSpeed()
+# Nov 8, 2005: - WATER_DEPTH => Z_BT
+# May 23, 2010: - Z* -> DEPTH*
+
+use strict;
+
+#======================================================================
+# fake_BT_RANGE(dta ptr)
+#======================================================================
+
+# During cruise NBP0204 it was found that one of Visbeck's RDIs consistently
+# returns zero as the bottom-track range, even thought the BT velocities
+# are good. This functions calculates the ranges if they are missing.
+
+sub cBTR($$$)
+{
+ my($d,$e,$b) = @_;
+ my($maxamp) = -9e99;
+ my($maxi);
+
+ for (my($i)=0; $i<$d->{N_BINS}; $i++) {
+ next unless ($d->{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b] > $maxamp);
+ $maxamp = $d->{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b];
+ $maxi = $i;
+ }
+ $d->{ENSEMBLE}[$e]->{BT_RANGE}[$b] =
+ $d->{DISTANCE_TO_BIN1_CENTER} + $maxi * $d->{BIN_LENGTH};
+}
+
+sub ensure_BT_RANGE($)
+{
+ my($d) = @_;
+ for (my($e)=0; $e<=$#{$d->{ENSEMBLE}}; $e++) {
+ my($sum) = my($n) = 0;
+ if (defined($d->{ENSEMBLE}[$e]->{BT_VELOCITY}[0])) {
+ for (my($b)=0; $b<4; $b++) {
+ cBTR($d,$e,$b)
+ unless defined($d->{ENSEMBLE}[$e]->{BT_RANGE}[$b]);
+ $sum += $d->{ENSEMBLE}[$e]->{BT_RANGE}[$b]; $n++;
+ }
+ } else {
+ for (my($b)=0; $b<4; $b++) {
+ $d->{ENSEMBLE}[$e]->{BT_RANGE}[$b] = undef;
+ }
+ }
+ $d->{ENSEMBLE}[$e]->{BT_MEAN_RANGE} = $sum/$n if ($n == 4);
+ }
+}
+
+#======================================================================
+# (seabed depth, stddev) = find_seabed(dta ptr, btm ensno, coord flag)
+#======================================================================
+
+# NOTE FOR YOYOS:
+# - this routine only detects the BT around the depeest depth!
+# - this is on purpose, because it is used for [listBT]
+
+# This is a PAIN:
+# if the instrument is too close to the bottom, the BT_RANGE
+# readings are all out; the only solution is to have a sufficiently
+# large search width (around the max(depth) ensemble) so that
+# a part of the up (oops, I forgot to finish this comment one year
+# ago during A0304 and now I don't understand it any more :-)
+
+my($search_width) = 200; # # of ensembles around bottom to search
+my($mode_width) = 10; # max range of bottom around mode
+my($z_offset) = 10000; # shift z to ensure +ve array indices
+
+sub find_seabed($$$)
+{
+ my($d,$be,$beamCoords) = @_;
+ my($i,$dd,$sd,$nd);
+ my(@guesses);
+
+ return undef unless ($be-$search_width >= 0 &&
+ $be+$search_width <= $#{$d->{ENSEMBLE}});
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless (defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[0]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[1]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[2]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[3]));
+ my(@BT) = $beamCoords
+ ? velInstrumentToEarth($d,$i,
+ velBeamToInstrument($d,
+ @{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}}))
+ : velApplyHdgBias($d,$i,@{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}});
+ next unless (abs($BT[3]) < 0.05);
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} =
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} =
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[0]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[1]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[2]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[3]/4;
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} *= -1
+ if ($d->{ENSEMBLE}[$i]->{XDUCER_FACING_UP});
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} += $d->{ENSEMBLE}[$i]->{DEPTH};
+ $guesses[int($d->{ENSEMBLE}[$i]->{DEPTH_BT})+$z_offset]++;
+ $nd++;
+ }
+ return undef unless ($nd>5);
+
+ my($mode,$nmax);
+ for ($i=0; $i<=$#guesses; $i++) { # find mode
+ $nmax=$guesses[$i],$mode=$i-$z_offset
+ if ($guesses[$i] > $nmax);
+ }
+
+ $nd = 0;
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ if (abs($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$mode) <= $mode_width) {
+ $dd += $d->{ENSEMBLE}[$i]->{DEPTH_BT};
+ $nd++;
+ } else {
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} = undef;
+ }
+ }
+ return undef unless ($nd >= 2);
+
+ $dd /= $nd;
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ $sd += ($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$dd)**2;
+ }
+
+ return ($dd, sqrt($sd/($nd-1)));
+}
+
+#======================================================================
+# c = soundSpeed()
+#======================================================================
+
+# Taken from the RDI BroadBand primer manual. The reference given there
+# is Urick (1983).
+
+sub soundSpeed($$$)
+{
+ my($salin,$temp,$depth) = @_;
+ return 1449.2 + 4.6*$temp -0.055*$temp**2 + 0.00029*$temp**3 +
+ (1.34 - 0.01*$temp) * ($salin - 35) + 0.016*$depth;
+}
+
+1;
new file mode 100755
--- /dev/null
+++ b/beamStats
@@ -0,0 +1,291 @@
+#!/usr/bin/perl
+#======================================================================
+# L I S T B I N S
+# doc: Fri Aug 25 15:57:05 2006
+# dlm: Thu Feb 21 15:05:40 2008
+# (c) 2006 A.M. Thurnherr
+# uE-Info: 32 57 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# Split data file into per-bin time series.
+
+# HISTORY:
+# Aug 25, 2006: - created from [listEns]
+# Aug 26, 2006: - added -M)agdecl
+# - changed -b to -f
+# Aug 27, 2006: - added %bin
+# Aug 28, 2006: - BUG: options were confused
+# Jan 4, 2007: - improved usage message for -a
+# - added %mag_decl
+# - BUG: roundoff error in %pct_good_vels
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl]
+# Feb 7, 2008: - added sound-speed correction
+# - enabled 3-beam solutions
+# Feb 8, 2008: - added -d)iscard <beam>
+# - added -b)eam coordinate output
+# Feb 12, 2008: - modified 3-beam output
+# - added -p)ct_good <min>
+# Feb 13, 2008: - various improvements
+# Feb 19, 2008: - BUG: division by zero
+# BUG: min() did not work with 1st elt undef
+# Feb 21, 2008: - BUG: had forgotten to undo debugging changes
+# - removed missing magdecl warning on -b
+
+# General Notes:
+# - everything (e.g. beams) is numbered from 1
+# - no support for BT data
+
+# Soundspeed Correction:
+# - applied as described in the RDI coord-trans manual
+# - sound-speed variation over range is ignored (valid for small gradients)
+# => - same simple correction for all velocity components
+# - simple correction for cell depths
+
+# Min %-good (min_pcg):
+# - nan for records w/o valid velocities
+# - min(%-good) of the beams used for the velocity solution
+# - min_pcg does not have to decrease monotonically with distance,
+# at least when 3-beam solutions are allowed and when -p is used to
+# edit the data
+# - non-monotonic min_pcg is particularly obvious with the DYNAMUCK BM_ADCP
+# 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";
+
+die("Usage: $0 [-r)ange <first_ens,last_ens>] " .
+ "[output -f)ile <pat[bin%d.raw]>] " .
+ "[-a)ll ens (not just those with good vels)] " .
+ "[-M)agnetic <declination>] " .
+ "[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
+ "[require -4)-beam solutions] [-d)iscard <beam#>] " .
+ "[-%)good <min>] " .
+ "[output -b)eam coordinates] " .
+ "<RDI file>\n")
+ unless (&Getopts("4abd:f:M:p:r:S:") && @ARGV == 1);
+
+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);
+
+$opt_f = 'bin%d.raw' unless defined($opt_f);
+$ifn = $ARGV[0];
+
+($first_ens,$last_ens) = split(',',$opt_r)
+ if defined($opt_r);
+
+($salin,$temp,$depth) = split(',',$opt_S)
+ if defined($opt_S);
+$variable_ssCorr = ($salin eq '*' || $temp eq '*' || $depth eq '*');
+
+#----------------------------------------------------------------------
+
+sub ssCorr($$$$) # sound velocity correction factor
+{
+ my($e,$S,$T,$D) = @_;
+ $S = $dta{ENSEMBLE}[$e]->{SALINITY} if ($S eq '*');
+ $T = $dta{ENSEMBLE}[$e]->{TEMPERATURE} if ($T eq '*');
+ if ($D eq '*') {
+ print(STDERR "WARNING: soundspeed correction using instrument pressure instead of depth!\n")
+ unless ($warned);
+ $warned = 1;
+ $D = $dta{ENSEMBLE}[$e]->{PRESSURE};
+ }
+ return soundSpeed($S,$T,$D) / $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
+}
+
+sub min(@) # return minimum
+{
+ my($min) = 99e99;
+ for (my($i)=0; $i<=$#_; $i++) {
+ $min = $_[$i] if defined($_[$i]) && ($_[$i] < $min);
+ }
+ return ($min == 99e99) ? nan : $min;
+}
+
+sub dumpBin($$$) # write time series of single bin
+{
+ my($b,$fe,$le) = @_;
+ my($file) = sprintf($opt_f,$b+1);
+
+ open(P,">$file") || die("$file: $!\n");
+ print(P "#ANTS#PARAMS# ");
+ foreach my $k (keys(%P)) {
+ print(P "$k\{$P{$k}\} ");
+ }
+ 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_good_vels{%.0f} ",100*$good_vels[$b]/($le-$fe+1));
+ printf(P "pct_3_beam{%.0f} ",$pct3b);
+ printf(P "bin{%d}",$b+1);
+ printf(P " soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!');
+ printf(P " dz{%g}",$dz[$b] *
+ (defined($opt_S) ? ssCorr($fe,$salin,$temp,$depth) : 1)
+ ) unless ($variable_ssCorr);
+ print( P "\n");
+
+ print(P "#ANTS#FIELDS# " .
+ "{ensemble} {date} {time} {elapsed} " .
+ "{heading} {pitch} {roll} " .
+ "{sig_heading} {sig_pitch} {sig_roll} " .
+ "{xmit_current} {xmit_voltage} " .
+ "{temp} " .
+ ($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") .
+ "{corr1} {corr2} {corr3} {corr4} " .
+ "{amp1} {amp2} {amp3} {amp4} " .
+ ($opt_b ? "{pcg1} {pcg2} {pcg3} {pcg4}" : "{3_beam} {min_pcg}")
+ );
+ print(P " {dz}") if ($variable_ssCorr);
+ print(P "\n");
+
+ my($t0) = $dta{ENSEMBLE}[$fe]->{UNIX_TIME};
+ for (my($e)=$fe; $e<=$le; $e++) {
+ next unless ($opt_a || $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]);
+
+ my($ssCorr) = defined($opt_S) ? ssCorr($e,$salin,$temp,$depth) : 1;
+
+ print(P "$dta{ENSEMBLE}[$e]->{NUMBER} ");
+ print(P "$dta{ENSEMBLE}[$e]->{DATE} ");
+ print(P "$dta{ENSEMBLE}[$e]->{TIME} ");
+ printf(P "%d ",$dta{ENSEMBLE}[$e]->{UNIX_TIME}-$t0);
+ print(P "$dta{ENSEMBLE}[$e]->{HEADING} ");
+ print(P "$dta{ENSEMBLE}[$e]->{PITCH} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ROLL} ");
+ print(P "$dta{ENSEMBLE}[$e]->{HEADING_STDDEV} ");
+ print(P "$dta{ENSEMBLE}[$e]->{PITCH_STDDEV} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ROLL_STDDEV} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_CURRENT} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE} ");
+ print(P "$dta{ENSEMBLE}[$e]->{TEMPERATURE} ");
+ if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$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);
+ printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] * $ssCorr);
+ } else {
+ print(P "nan nan nan nan ");
+ }
+ print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} ");
+ print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} ");
+ if ($opt_b) {
+ 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]}));
+ }
+ printf(P "%g ",$dz[$b]*$ssCorr) if ($variable_ssCorr);
+ print(P "\n");
+ }
+ close(P);
+}
+
+#----------------------------------------------------------------------
+# MAIN
+#----------------------------------------------------------------------
+
+$P{RDI_file} = $ifn;
+$P{mag_decl} = $opt_M if defined($opt_M);
+
+readData($ifn,\%dta);
+printf("%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}}));
+$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
+
+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")
+ if (!$dta{EARTH_COORDINATES});
+}
+
+$P{N_ensembles} = @{$dta{ENSEMBLE}};
+
+for (my($b)=0; $b<$dta{N_BINS}; $b++) { # calc dz
+ $dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
+}
+
+$lastGoodBin = 0;
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities
+ next if (defined($first_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
+ $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e
+ unless defined($P{first_ens});
+ last if (defined($last_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);
+ $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER};
+ $le = $e;
+
+ die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);
+ die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});
+ die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
+
+ for (my($b)=0; $b<$dta{N_BINS}; $b++) {
+ if (defined($opt_d)) {
+ undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
+ undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
+ }
+ 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);
+ $firstGoodEns = $e unless defined($firstGoodEns);
+ $lastGoodEns = $e;
+ }
+}
+
+unless (defined($opt_r)) {
+ $fe = $firstGoodEns;
+ $le = $lastGoodEns;
+}
+
+$firstBin = 0;
+$lastBin = $lastGoodBin;
+
+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);
+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: ");
+for ($b=$firstBin; $b<=$lastBin; $b++) { # generate output
+ dumpBin($b,$fe,$le);
+}
+print(STDERR "\n");
+
+exit(0);
new file mode 120000
--- /dev/null
+++ b/libRDI_Coords.pl
@@ -0,0 +1,1 @@
+RDI_Coords.pl
\ No newline at end of file
new file mode 100755
--- /dev/null
+++ b/listBT
@@ -0,0 +1,589 @@
+#!/usr/bin/perl
+#======================================================================
+# L I S T B T
+# doc: Sat Jan 18 18:41:49 2003
+# dlm: Sun May 23 16:33:33 2010
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 527 42 NIL 0 0 72 11 2 4 NIL ofnI
+#======================================================================
+
+# Extract Bottom-Track Data
+
+# NOTE: NO SOUND-SPEED CORRECTION APPLIED YET!!!
+
+# HISTORY:
+# Jan 18, 2003: - created
+# Jan 23, 2003: - added magnetic declination
+# Jan 25, 2003: - continued construction
+# Feb 11, 2003: - finally made it work
+# Feb 12, 2003: - added default profile output
+# Feb 13, 2003: - corrected raw output
+# Feb 14, 2003: - added errors if instrument-BT filters are more strict
+# than command-line values
+# Feb 18, 2003: - removed -d dependency on -W
+# Mar 3, 2003: - added -C)ompass correction
+# Mar 10, 2003: - added -f)orce to allow visbeck-style post processing
+# Mar 16, 2003: - added range comment
+# Feb 26, 2004: - added Earth-coordinate support
+# Feb 27, 2004: - made water-track calculation conditional (-E || -B)
+# Mar 9, 2004: - added magnetic_declination to %PARAMs
+# Apr 1, 2004: - added CTD u/v stats to %PARAMs
+# Apr 2, 2004: - added CTD_msf (mean square fluctuation) stat
+# Apr 3, 2004: - BUG: CTD vels were repeated for stats
+# - removed non-ANTS option
+# Nov 8, 2005: - UNIXTIME => UNIX_TIME
+# - adapted to new binary read library
+# - output editing statistics
+# Aug 15, 2006: - added -b
+# Aug 25, 2006: - fiddled
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
+# Nov 1, 2008: - BUG: sig(u) was reported instead of sig(v)
+# Jul 30, 2009: - NaN => nan
+
+# NOTES:
+# - the RDI BT data contains ranges that are greater than the
+# WT ping ranges. I don't know if those data are valid!
+# - there is a fair bit of heuristic used, especially in the
+# reference-layer calculation
+# - depth-correction (-m) is highly recommended because it allows
+# much better bad-BT detection and it is required for a valid
+# comparison with LADCP profiles
+# - the criterion for bottom-interference of the water-track data
+# is derived from Firing's [merge.c] (adding 1.5 bin lengths to
+# the calculated range), modified by taking the real beam angle
+# into account.
+# - from the RDI manuals it is not entirely clear whether the BT range
+# is given in vertical or in along-beam meters; comparison with the
+# WT range (calculated from the bin with the maximum echo amplitude)
+# shows that vertical meters are used
+
+# NOTES on quality checks:
+# -a minimum BT amplitude; setting this to 50 (RDI default is 30)
+# reduces the vertical range over which the bottom is detected but
+# not the quality of the bottom track data; therefore, this should
+# probably not be used.
+# -c minimum BT correlation; the RDI default for this parameter is 220,
+# which seems to work fine.
+# -e max error velocity (BT & WT); this is primarily used for detecting
+# good BT data, i.e. it should be set to a small value (Firing uses
+# 0.1m/s in merge); if too small a value is chosen too many good
+# data are discarded; note that the same error-velocity criterion
+# is used to separate good from bad data when mean profiles are
+# constructed.
+# -w max difference between reference-layer w and BT w; this is a
+# powerful criterion for determining good BT data; I like a value of
+# 0.03 m/s.
+# -d when the depth is corrected (-m) the...
+
+$0 =~ m{(.*)/[^/]+};
+require "$1/RDI_BB_Read.pl";
+require "$1/RDI_Coords.pl";
+require "$1/RDI_Utils.pl";
+use Getopt::Std;
+
+$USAGE = "$0 @ARGV";
+die("Usage: $0 " .
+ "[use -b)ins <1st,last>] " .
+ "[write -R)aw data] [write -B)T data] " .
+ "[write -E)nsembles <pref>] [-F)ilter ensembles <script>] " .
+ "[-C)ompass correction <amp/phase/bias>] " .
+ "[-w) <max-diff|0.03>] [-a)mp <min|30>] [-e)rr-vel <max|0.05>] " .
+ "[-c)orrelation <min|220>] " .
+ "[-W)ater <depth> [allowed -d)epth-diff <maxdiff|20>]] " .
+ "[-f)orce (no setup tests)] " .
+ "[-M)agnetic <declination>] " .
+ "<RDI file>\n")
+ unless (&getopts("BC:E:F:M:RW:a:b:c:d:e:fw:") && @ARGV == 1);
+
+print(STDERR "WARNING: magnetic declination not set!\n")
+ unless defined($opt_M);
+
+$opt_c = 220 unless defined($opt_c); # defaults
+$opt_a = 30 unless defined($opt_a);
+$opt_e = 0.05 unless defined($opt_e);
+$opt_w = 0.03 unless defined($opt_w);
+$opt_d = 20 unless defined($opt_d);
+
+if (defined($opt_C)) { # compass correction
+ ($CC_amp,$CC_phase,$CC_bias) = split('/',$opt_C);
+ die("$0: can't decode -C$opt_C\n")
+ unless defined($CC_bias);
+}
+
+unless ($opt_f) { # check BT setup
+ readHeader($ARGV[0],\%dta);
+ die("$0: minimum instrument BT correlation ($dta{BT_MIN_CORRELATION}) " .
+ "too large for selected criterion (-c $opt_c) --- use -f to override\n")
+ if ($dta{BT_MIN_CORRELATION} > $opt_c);
+ die("$0: minimum instrument BT amplitude ($dta{BT_MIN_EVAL_AMPLITUDE}) " .
+ "too large for selected criterion (-a $opt_a) --- use -f to override\n")
+ if ($dta{BT_MIN_EVAL_AMPLITUDE} > $opt_a);
+ die("$0: maximum instrument BT error velocity ($dta{BT_MAX_ERROR_VELOCITY}) " .
+ "too small for selected criterion (-e $opt_e) --- use -f to override\n")
+ if ($dta{BT_MAX_ERROR_VELOCITY} < $opt_e);
+}
+
+require $opt_F if defined($opt_F); # load filter
+
+print(STDERR "reading $ARGV[0]...");
+readData($ARGV[0],\%dta); # read data
+print(STDERR "done\n");
+
+$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
+ensure_BT_RANGE(\%dta); # make sure they're there
+
+$firstBin = $lastBin = '*'; # bins to use
+($firstBin,$lastBin) = split(',',$opt_b)
+ if defined($opt_b);
+$firstBin = 1 if ($firstBin eq '*');
+$lastBin = $dta{N_BINS} if ($lastBin eq '*');
+$firstBin--; $lastBin--;
+die("$ARGV[0]: not enough bins for ref layer\n")
+ unless ($lastBin-$firstBin >= 6);
+
+if ($dta{BEAM_COORDINATES}) { # coords used
+ $beamCoords = 1;
+} elsif (!$dta{EARTH_COORDINATES}) {
+ die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
+}
+
+#======================================================================
+# Calculate reference-layer w
+#======================================================================
+
+sub w($)
+{
+ my($ens) = @_;
+ my($i,$n,@v,$w);
+
+ for (my($b)=$firstBin; $b<=$lastBin; $b++) {
+ if ($beamCoords) {
+ next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][0] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][1] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][2] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][3] < 100);
+ @v = velInstrumentToEarth(\%dta,$ens,
+ velBeamToInstrument(\%dta,
+ @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$b]}));
+ } else {
+ next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][0] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][1] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][2] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][3] < 100);
+ @v = velApplyHdgBias(\%dta,$ens,
+ @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$b]});
+ }
+ next unless (defined($v[3]) && abs($v[3]) <= $opt_e);
+ $w += $v[2]; $n++;
+ }
+# printf(STDERR "$ens $n %.3f\n",$n>=1?$w/$n:-999);
+ return $n>=2 ? $w/$n : undef;
+}
+
+#======================================================================
+# Dump raw BT data from one ensemble
+#======================================================================
+
+sub dumpRaw($)
+{
+ my($e) = @_;
+
+ unless ($headerDone) {
+ print("#ANTS# [] $USAGE\n");
+ print("#ANTS#PARAMS# RDI_file{$ARGV[0]}\n");
+ print("#ANTS#FIELDS# {ens} {range1} {range2} {range3} {range4} " .
+ "{beamvel1} {beamvel2} {beamvel3} {beamvel4} {cor1} " .
+ "{cor2} {cor3} {cor4} {amp1} {amp2} {amp3} {amp4}\n");
+ $headerDone = 1;
+ }
+
+ printf("%d %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d\n",
+ $dta{ENSEMBLE}[$e]->{NUMBER},
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[0],
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[1],
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[2],
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[3],
+ $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0],
+ $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1],
+ $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2],
+ $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[3],
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0],
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1],
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2],
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3],
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[0],
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[1],
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[2],
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[3]
+ );
+}
+
+#======================================================================
+# Dump processed BT data from one ensemble
+#======================================================================
+
+sub dumpBT($)
+{
+ my($e) = @_;
+
+ unless ($headerDone) {
+ print("#ANTS# [] $USAGE\n");
+ printf("#ANTS#PARAMS# RDI_file{$ARGV[0]} bottom_time{%.1f}\n",
+ $dta{ENSEMBLE}[$maxz_e]->{ELAPSED});
+ print("#ANTS#FIELDS# {ens} {unix_time} {time} {depth} {BT_range} " .
+ "{WT_range} {u} {v} {w} {e} {w_ref} {corr} {amp}\n");
+ $headerDone = 1;
+ }
+
+ printf("%d %.2f %.2f %.1f %.1f %.1f %.4f %.4f %.4f %.4f %.4f %.1f %.1f\n",
+ $dta{ENSEMBLE}[$e]->{NUMBER},
+ $dta{ENSEMBLE}[$e]->{UNIX_TIME},
+ $dta{ENSEMBLE}[$e]->{ELAPSED},
+ $dta{ENSEMBLE}[$e]->{DEPTH},
+ $dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE},
+ $dta{ENSEMBLE}[$e]->{WT_MEAN_RANGE},
+ @{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}},
+ $dta{ENSEMBLE}[$e]->{W_REF},
+ $dta{ENSEMBLE}[$e]->{BT_MEAN_CORRELATION},
+ $dta{ENSEMBLE}[$e]->{BT_MEAN_EVAL_AMPLITUDE}
+ );
+}
+
+#======================================================================
+# Dump a single ensemble with valid BT data to separate file
+#======================================================================
+
+sub dumpEns(@) # write profile
+{
+ my($e) = @_;
+ my($b,$i);
+
+ open(P,">$opt_E.$e") || die("$opt_E.$e: $!\n");
+ print(P "#ANTS#PARAMS# " .
+ "depth{$dta{ENSEMBLE}[$e]->{DEPTH}} " .
+ "range{$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}} " .
+ "wt_range{$dta{ENSEMBLE}[$e]->{WT_MEAN_RANGE}} " .
+ "w_ref{$dta{ENSEMBLE}[$e]->{W_REF}} " .
+ "BT_u{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0]} " .
+ "BT_v{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1]} " .
+ "BT_w{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2]} " .
+ "BT_e{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[3]} " .
+ "BT_cor1{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0]} " .
+ "BT_cor2{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1]} " .
+ "BT_cor3{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2]} " .
+ "BT_cor4{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3]} " .
+ "BT_amp1{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[0]} " .
+ "BT_amp2{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[1]} " .
+ "BT_amp3{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[2]} " .
+ "BT_amp4{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[3]} " .
+ "BTFWT_u{$dta{ENSEMBLE}[$e]->{BTFWT_VELOCITY}[0]} " .
+ "BTFWT_v{$dta{ENSEMBLE}[$e]->{BTFWT_VELOCITY}[1]} " .
+ "BTFWT_w{$dta{ENSEMBLE}[$e]->{BTFWT_VELOCITY}[2]} " .
+ "\n"
+ );
+ print(P "#ANTS#FIELDS# " .
+ "{depth} {hab} {u} {v} {w} {e} {cor1} {cor2} {cor3} {cor4} " .
+ "{amp1} {amp2} {amp3} {amp4} {pcg1} {pcg2} {pcg3} {pcg4}\n"
+ );
+
+ my($slc) = (1-cos(rad($dta{BEAM_ANGLE})))*$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}
+ + 1.5*$dta{BIN_LENGTH}; # side-lobe contamination
+ for ($b=$firstBin; $b<=$lastBin; $b++) {
+ next unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]));
+ my($dz) = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
+ last if ($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}-$dz <= $slc);
+ my(@v) = $beamCoords
+ ? velInstrumentToEarth(\%dta,$e,
+ velBeamToInstrument(\%dta,
+ @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}))
+ : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
+ next unless defined($v[0]);
+ next if (abs($v[3]) > $opt_e ||
+ abs($v[2]-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2]) > 0.1);
+ $v[0] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0];
+ $v[1] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1];
+ my(@out) = (
+ $dta{ENSEMBLE}[$e]->{DEPTH}+$dz,
+ $dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}-$dz,
+ $v[0],$v[1],$v[2],$v[3],
+ @{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]},
+ @{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]},
+ @{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}
+ );
+ for ($i=0; $i<17; $i++) { $out[$i] = nan unless defined($out[$i]); }
+ print(P "@out\n");
+ }
+ close(P);
+}
+
+#======================================================================
+# Add Ensemble With Valid BT Data to Profile
+#======================================================================
+
+sub binEns($)
+{
+ my($e) = @_;
+ my($slc) = (1-cos(rad($dta{BEAM_ANGLE})))*$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}
+ + 1.5*$dta{BIN_LENGTH}; # side-lobe contamination
+ for (my($b)=$firstBin; $b<=$lastBin; $b++) {
+ next unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]));
+ my($dz) = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
+ last if ($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}-$dz <= $slc);
+ my(@v) = $beamCoords
+ ? velInstrumentToEarth(\%dta,$e,
+ velBeamToInstrument(\%dta,
+ @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}))
+ : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
+ next unless defined($v[0]);
+ next if (abs($v[3]) > $opt_e ||
+ abs($v[2]-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2]) > 0.1);
+
+ $v[0] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0];
+ $v[1] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1];
+
+ my($bin) = int(($dta{ENSEMBLE}[$e]->{DEPTH}+$dz) / $dta{BIN_LENGTH});
+ $minBin = $bin unless ($bin >= $minBin);
+ $maxBin = $bin unless ($bin <= $maxBin);
+ if (defined($BTn[$bin])) {
+ my($f1) = $BTn[$bin] / ($BTn[$bin]+1);
+ my($f2) = ($BTn[$bin]-1) / $BTn[$bin];
+ $BTsigu[$bin] =
+ $f2*$BTsigu[$bin] + $f1*($v[0]-$BTu[$bin])**2/$BTn[$bin];
+ $BTsigv[$bin] =
+ $f2*$BTsigv[$bin] + $f1*($v[1]-$BTv[$bin])**2/$BTn[$bin];
+ $BTn[$bin]++;
+ $BTu[$bin] = $f1*$BTu[$bin] + $v[0]/$BTn[$bin];
+ $BTv[$bin] = $f1*$BTv[$bin] + $v[1]/$BTn[$bin];
+ } else {
+ $BTu[$bin] = $v[0];
+ $BTv[$bin] = $v[1];
+ $BTn[$bin] = 1;
+ }
+ }
+}
+
+#======================================================================
+# Output Bottom-Referenced Profile
+#======================================================================
+
+sub dumpProf($$$)
+{
+ my($db,$wd,$md) = @_;
+ my(@sum,@mean);
+
+ for (my($i)=0; $i<=$#listCTDu; $i++) { # CTD vel mean
+ $sum[0] += $listCTDu[$i];
+ $sum[1] += $listCTDv[$i];
+ }
+ @mean = ($sum[0]/@listCTDu,$sum[1]/@listCTDv);
+ @sum = (0,0); # stddev
+ for (my($i)=0; $i<=$#listCTDu; $i++) {
+ $sum[0] += ($listCTDu[$i]-$mean[0])**2;
+ $sum[1] += ($listCTDv[$i]-$mean[1])**2;
+ }
+ @sigma = ($sum[0]/sqrt($#listCTDu),$sum[1]/sqrt($#listCTDv));
+ @sum = (0,0); # mean speed fluct
+ for (my($i)=1; $i<=$#listCTDu; $i++) { # also: list for median
+ push(@cfluc,sqrt(($listCTDu[$i]-$listCTDu[$i-1])**2 +
+ ($listCTDv[$i]-$listCTDv[$i-1])**2));
+ $sum[0] += $cfluc[$#cfluc];
+ }
+
+ printf("#ANTS#PARAMS# LADCP_depth_bias{%.1f} water_depth{%.1f} magnetic_declination{%.1f}\n",
+ $db,$wd,$md);
+ printf("#ANTS#PARAMS# CTD_u{%.3f} CTD_v{%.3f} CTD_sig_u{%.3f} CTD_sig_v{%.3f} CTD_mean_cfluc{%.4f} CTD_median_cfluc{%.4f}\n",
+ @mean,@sigma,$sum[0]/$#listCTDu,(sort{$a<=>$b}@cfluc)[@cfluc/2]);
+ printf("#ANTS#PARAMS# good_ens{$good}\n");
+ print("#ANTS#FIELDS# {depth} {u} {v} {sig_u} {sig_v} {n_data}\n");
+
+ for (my($bin)=$minBin; $bin<=$maxBin; $bin++) {
+ next unless defined($BTu[$bin]);
+ printf("%d %.3f %.3f %.3f %.3f %d\n",
+ ($bin+0.5)*$dta{BIN_LENGTH},
+ $BTu[$bin],$BTv[$bin],
+ sqrt($BTsigu[$bin]),sqrt($BTsigv[$bin]),
+ $BTn[$bin]);
+ }
+}
+
+#======================================================================
+# STEP 1: Calculate Depth (integrate w)
+#======================================================================
+
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {
+ checkEnsemble(\%dta,$e);
+ $dta{ENSEMBLE}[$e]->{W_REF} = w($e);
+ next unless (defined($start_e) ||
+ defined($dta{ENSEMBLE}[$e]->{W_REF}));
+ $start_e = $e unless defined($start_e);
+ $end_e = $e if defined($dta{ENSEMBLE}[$e]->{W_REF});
+ $lasttime = $curtime;
+ $curtime = $dta{ENSEMBLE}[$e]->{UNIX_TIME};
+ $dta{ENSEMBLE}[$e]->{ELAPSED} =
+ $curtime - $dta{ENSEMBLE}[$start_e]->{UNIX_TIME};
+ filterEnsemble(\%dta,$e)
+ if (defined($opt_F) &&
+ $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[0][3] > 0);
+ $z += $dta{ENSEMBLE}[$e]->{W_REF} *
+ (defined($lasttime) ? ($curtime - $lasttime) : 0);
+ $maxz_e=$e,$maxz = $z unless ($z < $maxz);
+ $dta{ENSEMBLE}[$e]->{DEPTH} = $z;
+}
+
+unless ($opt_R) {
+ ($w_depth,$swd) = find_seabed(\%dta,$maxz_e,$beamCoords);
+ die("$0: can't determine water depth (sigma = $swd)\n")
+ unless (defined($w_depth) && $swd < 10);
+
+ if (defined($opt_W)) { # adjust depth
+ $zbias = $w_depth - $opt_W; $w_depth = $opt_W;
+ for ($e=$start_e; $e<=$end_e; $e++) {
+ $dta{ENSEMBLE}[$e]->{DEPTH} -= $zbias;
+ }
+ }
+}
+
+# print(STDERR "maxz = $maxz, w_depth = $w_depth\n");
+
+#======================================================================
+# STEP 2: Process BT Data
+#======================================================================
+
+for ($e=$start_e; $e<=$end_e; $e++) {
+ next unless ($dta{ENSEMBLE}[$e]->{BT_RANGE}[0] && # BT data available
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[1] &&
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[2] &&
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[3]);
+# die("$0: don't know what to do with non-zero %-good and " .
+# "signal-strength values at ensemble " .
+# "#$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+# if ($dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[0] ||
+# $dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[1] ||
+# $dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[2] ||
+# $dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[3] ||
+# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[0] ||
+# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[1] ||
+# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[2] ||
+# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[3]);
+
+ if ($opt_R) { # dump raw data
+ dumpRaw($e);
+ next;
+ }
+
+ $dta{ENSEMBLE}[$e]->{HEADING} -= # compass correction
+ $CC_amp * sin(rad($dta{ENSEMBLE}[$e]->{HEADING} - $CC_phase))
+ + $CC_bias
+ if defined($opt_C);
+
+ @{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}} = $beamCoords # xform BT vel
+ ? velInstrumentToEarth(\%dta,$e,
+ velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}}))
+ : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}});
+
+ $dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE} = # mean vals
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[0]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[1]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[2]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_RANGE}[3]/4;
+ $dta{ENSEMBLE}[$e]->{BT_MEAN_CORRELATION} =
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3]/4;
+ $dta{ENSEMBLE}[$e]->{BT_MEAN_EVAL_AMPLITUDE} =
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[0]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[1]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[2]/4 +
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[3]/4;
+
+# next # could add this
+# if ($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE} < 50);
+
+ $bad_amp++,next
+ if ($dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[0] < $opt_a ||
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[1] < $opt_a ||
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[2] < $opt_a ||
+ $dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[3] < $opt_a);
+ $bad_corr++,next
+ if ($dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0] < $opt_c ||
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1] < $opt_c ||
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2] < $opt_c ||
+ $dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3] < $opt_c);
+
+ $bad_w_ref++,next # quality checks
+ if (abs($dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2] -
+ $dta{ENSEMBLE}[$e]->{W_REF}) > $opt_w);
+ $bad_e_vel++,next
+ if (abs($dta{ENSEMBLE}[$e]->{BT_VELOCITY}[3]) > $opt_e);
+ $bad_depth++,next
+ if (abs($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE} +
+ $dta{ENSEMBLE}[$e]->{DEPTH} - $w_depth) > $opt_d);
+
+ $good++;
+ push(@listCTDu,-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0]);
+ push(@listCTDv,-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1]);
+
+ if ($opt_E || $opt_B) {
+ my(@maxamp) = (0,0,0,0); # water-track range
+ my(@btm_e) = (0,0,0,0);
+ for ($b=$firstBin; $b<=$lastBin; $b++) {
+ for ($i=0; $i<4; $i++) {
+ if ($dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][$i] > $maxamp[$i]) {
+ $dta{ENSEMBLE}[$e]->{WT_RANGE}[$i] = $b;
+ $maxamp[$i] = $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][$i];
+ $btm_e[$i] = $e;
+ }
+ }
+ }
+ for ($i=0; $i<4; $i++) {
+ $dta{ENSEMBLE}[$e]->{WT_RANGE}[$i] *= $dta{BIN_LENGTH};
+ $dta{ENSEMBLE}[$e]->{WT_RANGE}[$i] += $dta{DISTANCE_TO_BIN1_CENTER};
+ }
+ $dta{ENSEMBLE}[$e]->{WT_MEAN_RANGE} =
+ $dta{ENSEMBLE}[$e]->{WT_RANGE}[0]/4 +
+ $dta{ENSEMBLE}[$e]->{WT_RANGE}[1]/4 +
+ $dta{ENSEMBLE}[$e]->{WT_RANGE}[2]/4 +
+ $dta{ENSEMBLE}[$e]->{WT_RANGE}[3]/4;
+
+ my($btm_e) = int($btm_e[0]/4+$btm_e[1]/4+$btm_e[2]/4+$btm_e[3]/4+0.5);
+ @{$dta{ENSEMBLE}[$btm_e]->{BTFWT_VELOCITY}} = $beamCoords # BT from WT
+ ? velInstrumentToEarth(\%dta,$e,
+ velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}}))
+ : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}});
+ }
+
+ dumpEns($e) if defined($opt_E); # output BT profiles
+ if ($opt_B) { dumpBT($e); }
+ else { binEns($e); }
+}
+
+filterEnsembleStats() if defined($opt_F);
+exit(0) if ($opt_R);
+
+printf(STDERR "%5d BT records removed due to bad w\n",$bad_w_ref)
+ if defined($bad_w_ref);
+printf(STDERR "%5d BT records removed due to bad err vel\n",$bad_e_vel)
+ if defined($bad_e_vel);
+printf(STDERR "%5d BT records removed due to bad echo amplitude\n",$bad_amp)
+ if defined($bad_amp);
+printf(STDERR "%5d BT records removed due to bad correlation\n",$bad_corr)
+ if defined($bad_corr);
+printf(STDERR "%5d BT records removed due to bad depth\n",$bad_depth)
+ if defined($bad_depth);
+
+die("$0: no good BT data\n") unless ($good);
+
+printf(STDERR "\n%5d BT records remaining\n",$good);
+
+dumpProf($zbias,$w_depth,-$dta{HEADING_BIAS})
+ unless ($opt_B);
+
+exit(0);
+
new file mode 100755
--- /dev/null
+++ b/listBins
@@ -0,0 +1,330 @@
+#!/usr/bin/perl
+#======================================================================
+# L I S T B I N S
+# doc: Fri Aug 25 15:57:05 2006
+# dlm: Sat May 23 16:34:30 2009
+# (c) 2006 A.M. Thurnherr
+# uE-Info: 269 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# Split data file into per-bin time series.
+
+# HISTORY:
+# Aug 25, 2006: - created from [listEns]
+# Aug 26, 2006: - added -M)agdecl
+# - changed -b to -f
+# Aug 27, 2006: - added %bin
+# Aug 28, 2006: - BUG: options were confused
+# Jan 4, 2007: - improved usage message for -a
+# - added %mag_decl
+# - BUG: roundoff error in %pct_good_vels
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl]
+# Feb 7, 2008: - added sound-speed correction
+# - enabled 3-beam solutions
+# Feb 8, 2008: - added -d)iscard <beam>
+# - added -b)eam coordinate output
+# Feb 12, 2008: - modified 3-beam output
+# - added -p)ct_good <min>
+# Feb 13, 2008: - various improvements
+# Feb 19, 2008: - BUG: division by zero
+# BUG: min() did not work with 1st elt undef
+# Feb 21, 2008: - BUG: had forgotten to undo debugging changes
+# - removed missing magdecl warning on -b
+# Feb 22, 2008: - moved ssCorr() to [RDI_Utils.pl]
+# - BUG: %d complete ensembles was written to STDOUT
+# - BUG: %N_ensembles was total number of ensembles, not
+# only the ones used
+# - BUG: 0 errvel was output for 3-beam solutions =>
+# wrong time-average statistics
+# May 19, 2009: - added -w to calculate vertical velocities from
+# two beam pairs separately
+# May 21, 2009: - added horizontal beampair velocities on -w
+# - -P)itchRoll <bias/bias>
+# May 22, 2009: - added -B) <bias/bias/bias/bias>
+# May 23, 2009: - adapted to changed beampair-velocity fun name
+
+# General Notes:
+# - everything (e.g. beams) is numbered from 1
+# - no support for BT data
+
+# Soundspeed Correction:
+# - applied as described in the RDI coord-trans manual
+# - sound-speed variation over range is ignored (valid for small gradients)
+# => - same simple correction for all velocity components
+# - simple correction for cell depths
+
+# Min %-good (min_pcg):
+# - nan for records w/o valid velocities
+# - min(%-good) of the beams used for the velocity solution
+# - min_pcg does not have to decrease monotonically with distance,
+# at least when 3-beam solutions are allowed and when -p is used to
+# edit the data
+# - non-monotonic min_pcg is particularly obvious with the DYNAMUCK BM_ADCP
+# 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";
+
+die("Usage: $0 [-r)ange <first_ens,last_ens>] " .
+ "[output -f)ile <pat[bin%d.raw]>] " .
+ "[-a)ll ens (not just those with good vels)] " .
+ "[-M)agnetic <declination>] " .
+ "[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
+ "[-P)itch/Roll <bias/bias>] [-B)eamvel <bias/bias/bias/bias>] " .
+ "[require -4)-beam solutions] [-d)iscard <beam#>] " .
+ "[-%)good <min>] " .
+ "[output -b)eam coordinates] [output two separate -w) estimates] " .
+ "<RDI file>\n")
+ unless (&Getopts("4abB:d:f:M:p:r:P:S:w") && @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);
+
+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);
+
+$opt_f = 'bin%d.raw' unless defined($opt_f);
+$ifn = $ARGV[0];
+
+($first_ens,$last_ens) = split(',',$opt_r)
+ if defined($opt_r);
+
+($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S)
+ if defined($opt_S);
+$variable_ssCorr = ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*');
+
+#----------------------------------------------------------------------
+
+sub min(@) # return minimum
+{
+ my($min) = 99e99;
+ for (my($i)=0; $i<=$#_; $i++) {
+ $min = $_[$i] if defined($_[$i]) && ($_[$i] < $min);
+ }
+ return ($min == 99e99) ? nan : $min;
+}
+
+sub dumpBin($$$) # write time series of single bin
+{
+ my($b,$fe,$le) = @_;
+ my($file) = sprintf($opt_f,$b+1);
+
+ open(P,">$file") || die("$file: $!\n");
+ print(P "#ANTS#PARAMS# ");
+ foreach my $k (keys(%P)) {
+ print(P "$k\{$P{$k}\} ");
+ }
+ if ($opt_b) {
+ printf(STDERR "%.0f%% ",100*$good_vels[$b]/($le-$fe+1));
+ } else {
+ 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(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!');
+ printf(P " dz{%g}",$dz[$b] *
+ (defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) : 1)
+ ) unless ($variable_ssCorr);
+ print( P "\n");
+
+ print(P "#ANTS#FIELDS# " .
+ "{ensemble} {date} {time} {elapsed} " .
+ "{heading} {pitch} {roll} " .
+ "{sig_heading} {sig_pitch} {sig_roll} " .
+ "{xmit_current} {xmit_voltage} " .
+ "{temp} " .
+ ($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") .
+ ($opt_w ? "{v12} {w12} {v34} {w34} " : "" ) .
+ "{corr1} {corr2} {corr3} {corr4} " .
+ "{amp1} {amp2} {amp3} {amp4} " .
+ ($opt_b ? "{pcg1} {pcg2} {pcg3} {pcg4}" : "{3_beam} {min_pcg}")
+ );
+ print(P " {dz}") if ($variable_ssCorr);
+ print(P "\n");
+
+ my($t0) = $dta{ENSEMBLE}[$fe]->{UNIX_TIME};
+ for (my($e)=$fe; $e<=$le; $e++) {
+ next unless ($opt_a || $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]);
+
+ my($ssCorr) = defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$e],$SS_salin,$SS_temp,$SS_depth) : 1;
+
+ print(P "$dta{ENSEMBLE}[$e]->{NUMBER} ");
+ print(P "$dta{ENSEMBLE}[$e]->{DATE} ");
+ print(P "$dta{ENSEMBLE}[$e]->{TIME} ");
+ printf(P "%d ",$dta{ENSEMBLE}[$e]->{UNIX_TIME}-$t0);
+ print(P "$dta{ENSEMBLE}[$e]->{HEADING} ");
+ print(P "$dta{ENSEMBLE}[$e]->{PITCH} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ROLL} ");
+ print(P "$dta{ENSEMBLE}[$e]->{HEADING_STDDEV} ");
+ print(P "$dta{ENSEMBLE}[$e]->{PITCH_STDDEV} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ROLL_STDDEV} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_CURRENT} ");
+ print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE} ");
+ print(P "$dta{ENSEMBLE}[$e]->{TEMPERATURE} ");
+ if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$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);
+ if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]) {
+ print(P "nan ");
+ } else {
+ printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] * $ssCorr);
+ }
+ if ($opt_w) {
+ printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]) ? "%g " : "nan ",
+ $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]);
+ printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]) ? "%g " : "nan ",
+ $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]);
+ printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]) ? "%g " : "nan ",
+ $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]);
+ printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]) ? "%g " : "nan ",
+ $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]);
+ }
+ } else {
+ print(P "nan nan nan nan ");
+ print(P "nan nan nan nan ") if ($opt_w);
+ }
+ print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} ");
+ print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} ");
+ if ($opt_b) {
+ 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]}));
+ }
+ printf(P "%g ",$dz[$b]*$ssCorr) if ($variable_ssCorr);
+ print(P "\n");
+ }
+ close(P);
+}
+
+#----------------------------------------------------------------------
+# MAIN
+#----------------------------------------------------------------------
+
+$P{RDI_file} = $ifn;
+$P{mag_decl} = $opt_M if defined($opt_M);
+
+readData($ifn,\%dta);
+printf(STDERR "%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}}));
+$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
+
+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")
+ if (!$dta{EARTH_COORDINATES});
+}
+
+for (my($b)=0; $b<$dta{N_BINS}; $b++) { # calc dz
+ $dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
+}
+
+$lastGoodBin = 0;
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities
+ next if (defined($first_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
+
+ $dta{ENSEMBLE}[$e]->{PITCH} -= $P{pitch_bias};
+ $dta{ENSEMBLE}[$e]->{ROLL} -= $P{roll_bias};
+
+ $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e
+ unless defined($P{first_ens});
+ last if (defined($last_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);
+ $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER};
+ $le = $e;
+
+ die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);
+ die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});
+ die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
+
+ for (my($b)=0; $b<$dta{N_BINS}; $b++) {
+ $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] -= $P{velbias_b1};
+ $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] -= $P{velbias_b2};
+ $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] -= $P{velbias_b3};
+ $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] -= $P{velbias_b4};
+ if ($opt_w) {
+ @{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =
+ velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
+ }
+
+ if (defined($opt_d)) {
+ undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
+ undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
+ }
+ 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);
+ $firstGoodEns = $e unless defined($firstGoodEns);
+ $lastGoodEns = $e;
+ }
+}
+
+unless (defined($opt_r)) {
+ $fe = $firstGoodEns;
+ $le = $lastGoodEns;
+}
+
+$P{N_ensembles} = $le - $fe + 1;
+
+$firstBin = 0;
+$lastBin = $lastGoodBin;
+
+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 {
+ 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: ");
+}
+for ($b=$firstBin; $b<=$lastBin; $b++) { # generate output
+ dumpBin($b,$fe,$le);
+}
+print(STDERR "\n");
+
+exit(0);
new file mode 100755
--- /dev/null
+++ b/listEns
@@ -0,0 +1,278 @@
+#!/usr/bin/perl
+#======================================================================
+# L I S T E N S
+# doc: Sat Jan 18 18:41:49 2003
+# dlm: Thu Jul 30 17:42:09 2009
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 247 47 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# Print useful info from the ensemble list or dump ensembles to
+# separate files.
+
+# HISTORY:
+# Jan 18, 2003: - created
+# Mar 18, 2004: - updated
+# Sep 15, 2005: - made ESW optional (BB150)
+# - change RDI binread library name
+# Aug 25, 2006: - added -r)ange
+# - added write -E)nsembles
+# Aug 26, 2006: - added -M)agdecl
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
+# Jan 26, 2008: - BUG: diagnostic output had been written to STDOUT
+# Feb 1, 2008: - BUG: still more diagnostic output written to STDOUT
+# - BUG: -E/-A combo had ignored -E
+# - changed %-good fieldnames for earth coordinates
+# - allowed for 3-beam solutions
+# Feb 7, 2008: - added -f)ields
+# Apr 4, 2008: - made -f output nan on undefined values
+# - BUG: -f fields did not allow array indices
+# - added in-w)ater data only
+# - restructured for simplicity
+# Mar 2, 2009: - added # of valid bin-1 vels to non-ANTS output
+# Jul 30, 2009: - NaN => nan
+
+# Notes:
+# - -E outputs data in earth coordinates
+# - -E output is always in ANTS format, ignoring -A
+# - no soundspeed correction
+
+require "getopts.pl";
+$0 =~ m{(.*/)[^/]+};
+require "$1RDI_BB_Read.pl";
+require "$1RDI_Coords.pl";
+
+die("Usage: $0 [-A)nts] [-Q)uiet (errcheck only)] " .
+ "[-f)ields <[name=]FIELD[,...]>] " .
+ "[write -E)nsemples <pref> [-M)agnetic <declination>] [min -p)ercent-good <#>]] " .
+ "[-r)ange <first_ens,last_ens>] [in-w)ater ensembles only]" .
+ "<RDI file...>\n")
+ unless (&Getopts("AE:f:M:p:Qr:w") && $#ARGV >= 0);
+
+print(STDERR "WARNING: magnetic declination not set!\n")
+ if defined($opt_E) && !defined($opt_M);
+
+die("$0: illegal option combination\n")
+ if ($opt_Q && $opt_A) || ((defined($opt_M) || defined($opt_p)) && !defined($opt_E));
+
+($first_ens,$last_ens) = split(',',$opt_r)
+ if defined($opt_r);
+
+undef($opt_A) if defined($opt_E);
+
+$opt_p = 0 unless defined($opt_p);
+
+if ($opt_f) { # additional fields
+ @addFields = split(',',$opt_f);
+ foreach my $f (@addFields) {
+ $f =~ s/\s//g; # remove spaces
+ @def = split('=',$f);
+ if (@def == 2) { # name=field
+ $addLayout .= $opt_A ? " {$def[0]}" : " $def[0]";
+ $f = $def[1];
+ } else { # field
+ $addLayout .= " {$f}";
+ }
+ }
+# print(STDERR "addLayout = $addLayout\n");
+# print(STDERR "\@addFields = @addFields\n");
+}
+
+#----------------------------------------------------------------------
+# MAIN
+#----------------------------------------------------------------------
+
+while (-f $ARGV[0]) {
+ if ($opt_A && !$opt_E) {
+ print("#ANTS#PARAMS# RDI_file{$ARGV[0]}\n");
+ } elsif (!$opt_Q) {
+ print(STDERR "$ARGV[0]: ");
+ }
+ readData(@ARGV,\%dta);
+ printf(STDERR "%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}}))
+ unless ($opt_Q);
+ $dta{HEADING_BIAS} = -$opt_M; # magnetic declination
+ shift;
+
+ if ($dta{BEAM_COORDINATES}) { # coords used
+ $beamCoords = 1;
+ } elsif (!$dta{EARTH_COORDINATES}) {
+ die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
+ }
+
+ if ($opt_A) { # select output fmt: ANTS
+ unless ($opt_Q) {
+ printf("#ANTS#PARAMS# N_ensembles{%d}\n",scalar(@{$dta{ENSEMBLE}}));
+ print('#ANTS#FIELDS# {ens} {time} {xducer_up} {temp} {hdg} {pitch} {roll} {XMIT_VOLTAGE} {XMIT_CURRENT}');
+ print(' {ESW}') if ($dta{DATA_FORMAT} eq 'WH300');
+ print("$addLayout\n");
+ }
+
+ $dumpEns = sub ($)
+ {
+ my($e) = @_;
+
+ printf('%d %lf %d %g %g %g %g %g %g',
+ $dta{ENSEMBLE}[$e]->{NUMBER},
+ $dta{ENSEMBLE}[$e]->{UNIX_TIME},
+ $dta{ENSEMBLE}[$e]->{XDUCER_FACING_UP} ? 1 : 0,
+ $dta{ENSEMBLE}[$e]->{TEMPERATURE},
+ $dta{ENSEMBLE}[$e]->{HEADING},
+ $dta{ENSEMBLE}[$e]->{PITCH},
+ $dta{ENSEMBLE}[$e]->{ROLL},
+ $dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE},
+ $dta{ENSEMBLE}[$e]->{ADC_XMIT_CURRENT},
+ );
+ printf(' %08X',$dta{ENSEMBLE}[$e]->{ERROR_STATUS_WORD})
+ if ($dta{DATA_FORMAT} eq 'WH300');
+ foreach my $f (@addFields) {
+ my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
+ $fn = $f unless defined($fn);
+ my($v) = eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi");
+ print(defined($v) ? " $v" : " nan");
+ }
+ print("\n");
+ }
+
+ } elsif ($opt_E) { # one file per ens
+
+ $dumpEns = sub ($)
+ {
+ my($e) = @_;
+ my($b,$i);
+ my($file) = "$opt_E.$dta{ENSEMBLE}[$e]->{NUMBER}";
+
+ open(P,">$file") || die("$file: $!\n");
+ print(P "#ANTS#PARAMS# " .
+ "BT_u{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0]} " .
+ "BT_v{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1]} " .
+ "BT_w{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2]} " .
+ "BT_e{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[3]} " .
+ "BT_cor1{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0]} " .
+ "BT_cor2{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1]} " .
+ "BT_cor3{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2]} " .
+ "BT_cor4{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3]} " .
+ "BT_amp1{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[0]} " .
+ "BT_amp2{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[1]} " .
+ "BT_amp3{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[2]} " .
+ "BT_amp4{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[3]} " .
+ "\n"
+ );
+ print(P "#ANTS#FIELDS# " .
+ "{dz} {u} {v} {w} {e} {cor1} {cor2} {cor3} {cor4} " .
+ "{amp1} {amp2} {amp3} {amp4} "
+ );
+ if ($beamCoords) {
+ print(P "{pcg1} {pcg2} {pcg3} {pcg4}");
+ } else {
+ print(P "{pc3beam} {pcBadErrVel} {pc1or2beam} {pc4beam}");
+ }
+ print(P "$addLayout\n");
+
+ for (my($b)=0; $b<$dta{N_BINS}; $b++) {
+ my(@v);
+ my($dz) = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
+
+ if ($beamCoords) {
+ undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])
+ if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0] < $opt_p);
+ undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1])
+ if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1] < $opt_p);
+ undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2])
+ if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2] < $opt_p);
+ undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])
+ if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3] < $opt_p);
+ @v = velInstrumentToEarth(\%dta,$e,
+ velBeamToInstrument(\%dta,
+ @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}));
+ } else {
+ @v = velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
+ }
+
+ $v[0] = nan unless defined($v[0]);
+ $v[1] = nan unless defined($v[1]);
+ $v[2] = nan unless defined($v[2]);
+ $v[3] = nan unless defined($v[3]);
+ my(@out) = (
+ $dz,$v[0],$v[1],$v[2],$v[3],
+ @{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]},
+ @{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]},
+ @{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}
+ );
+ foreach my $f (@addFields) {
+ my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
+ $fn = $f unless defined($fn);
+ push(@out,eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));
+ }
+ for ($i=0; $i<17+@addFields; $i++) {
+ $out[$i] = nan unless defined($out[$i]);
+ }
+ print(P "@out\n");
+ }
+ close(P);
+ }
+
+ } else { # neither ANTS nor ens files
+ unless ($opt_Q) {
+ if ($dta{DATA_FORMAT} eq 'WH300') {
+ printf(" # Date Time XD Temp Headng Pitch Roll vels(bin1) ESW$addLayout\n");
+ printf("-----------------------------------------------------------------------\n");
+ } else {
+ printf(" # Date Time XD Temp Headng Pitch Roll vels(bin1)$addLayout\n");
+ printf("-------------------------------------------------------------------\n");
+ }
+ }
+
+ $dumpEns = sub ($)
+ {
+ my($e) = @_;
+
+ printf('%5d %s %s %s %5.1f %6.1f %5.1f %5.1f %3d',
+ $dta{ENSEMBLE}[$e]->{NUMBER},
+ $dta{ENSEMBLE}[$e]->{DATE},
+ $dta{ENSEMBLE}[$e]->{TIME},
+ $dta{ENSEMBLE}[$e]->{XDUCER_FACING_UP} ? "UP" : "DN",
+ $dta{ENSEMBLE}[$e]->{TEMPERATURE},
+ $dta{ENSEMBLE}[$e]->{HEADING},
+ $dta{ENSEMBLE}[$e]->{PITCH},
+ $dta{ENSEMBLE}[$e]->{ROLL},
+ $dta{ENSEMBLE}[$e]->{BIN1VELS},
+ );
+ printf(' 0x%08X',$dta{ENSEMBLE}[$e]->{ERROR_STATUS_WORD})
+ if ($dta{DATA_FORMAT} eq 'WH300');
+ foreach my $f (@addFields) {
+ my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
+ $fn = $f unless defined($fn);
+ my($v) = eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi");
+ print(defined($v) ? " $v" : " nan");
+ }
+ print("\n");
+ }
+
+ }
+
+ for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {
+ next if (defined($first_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
+ last if (defined($last_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);
+ $dta{ENSEMBLE}[$e]->{BIN1VELS} =
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[1][0]) +
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[1][1]) +
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[1][2]) +
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[1][3]);
+ next if ($opt_w && $dta{ENSEMBLE}[$e]->{BIN1VELS}<3);
+
+ die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);
+ die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});
+ die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
+
+ &$dumpEns($e)
+ unless ($opt_Q);
+ }
+}
+
+exit(0);
new file mode 100755
--- /dev/null
+++ b/listHdr
@@ -0,0 +1,174 @@
+#!/usr/bin/perl
+#======================================================================
+# L I S T H D R
+# doc: Sat Jan 18 18:41:49 2003
+# dlm: Wed Jul 9 13:30:35 2008
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 66 0 NIL 0 0 72 8 2 4 NIL ofnI
+#======================================================================
+
+# Print useful info from the RDI BB header
+
+# HISTORY:
+# Jan 18, 2003: - incepted as a test for [WorkhorseBinRead.pl]
+# Jan 19, 2003: - continued
+# Feb 14, 2003: - added BT setup params
+# Mar 15, 2003: - added and removed BATTERY
+# Feb 24, 2004: - corrected TRANSMIT_LAG_DISTANCE units
+# Feb 26, 2004: - added ENSEMBLE_LENGTH
+# - added FIRMWARE
+# Mar 4, 2004: - added transducer orientation
+# Mar 30, 2004: - decified firmware output
+# Sep 14, 2005: - made BT data optional, dep. on NUMBER_OF_DATA_TYPES
+# - added DATA_FORMAT & related
+# Sep 15, 2005: - change BinRead library name
+# - compacted output format
+# Oct 30, 2005: - shuffled stuff, added DATA_FORMAT_VARIANT
+# Aug 21, 2006: - added CPU_SERIAL_NUMBER
+# - added usage error
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
+# Jul 9, 2008: - added output regarding available sensors
+
+$0 =~ m{(.*/)[^/]+};
+require "$1RDI_BB_Read.pl";
+
+die("Usage: $0 <RDI file[...]>\n")
+ unless (@ARGV);
+
+while (-f $ARGV[0]) {
+ print("$ARGV[0]:\n");
+ readHeader(@ARGV,\%hdr);
+ shift;
+
+ print(" Instrument Characteristics:\n");
+
+ printf("\tCPU_SERIAL_NUMBER\t\t= %s\n",$hdr{CPU_SERIAL_NUMBER});
+ printf("\tFIRMWARE\t\t\t= %d.%d\n",$hdr{CPU_FW_VER},$hdr{CPU_FW_REV});
+ printf("\tBEAM_FREQUENCY\t\t\t= %d kHz\n",$hdr{BEAM_FREQUENCY});
+ printf("\tBEAM_ANGLE\t\t\t= %d deg\n",$hdr{BEAM_ANGLE});
+ printf("\tN_BEAMS\t\t\t\t= %d\n",$hdr{N_BEAMS});
+ printf("\tN_DEMODS\t\t\t= %d\n",$hdr{N_DEMODS}) if defined($hdr{N_DEMODS});
+
+ printf("\tSensors\t\t\t\t: ");
+ printf("PRESSURE ") if ($hdr{PRESSURE_SENSOR_AVAILABLE});
+ printf("CONDUCTIVITY ") if ($hdr{CONDUCTIVITY_SENSOR_AVAILABLE});
+ printf("TEMPERATURE ") if ($hdr{TEMPERATURE_SENSOR_AVAILABLE});
+ printf("COMPASS ") if ($hdr{COMPASS_AVAILABLE});
+ printf("PITCH ") if ($hdr{PITCH_SENSOR_AVAILABLE});
+ printf("ROLL ") if ($hdr{ROLL_SENSOR_AVAILABLE});
+ print("\n");
+
+ printf("\tFlags\t\t\t\t: ");
+ printf("XDUCER_HEAD_ATTACHED ") if ($hdr{XDUCER_HEAD_ATTACHED});
+ printf("CONVEX_BEAM_PATTERN ") if ($hdr{CONVEX_BEAM_PATTERN});
+ printf("CONCAVE_BEAM_PATTERN ") if ($hdr{CONCAVE_BEAM_PATTERN});
+ print("\n");
+
+
+ print(" File Format:\n");
+
+ printf("\tDATA_FORMAT\t\t\t= %s (variant %d)\n",
+ $hdr{DATA_FORMAT},$hdr{DATA_FORMAT_VARIANT});
+ printf("\tNUMBER_OF_DATA_TYPES\t\t= %d\n",$hdr{NUMBER_OF_DATA_TYPES});
+ printf("\tENSEMBLE_BYTES\t\t\t= %3d bytes\n",$hdr{ENSEMBLE_BYTES});
+ printf("\tHEADER_BYTES\t\t\t= %3d bytes\n",$hdr{HEADER_BYTES});
+ printf("\tFIXED_LEADER_BYTES\t\t= %3d bytes\n",$hdr{FIXED_LEADER_BYTES});
+ printf("\tVARIABLE_LEADER_BYTES\t\t= %3d bytes\n",$hdr{VARIABLE_LEADER_BYTES});
+ printf("\tVELOCITY_DATA_BYTES\t\t= %3d bytes\n",$hdr{VELOCITY_DATA_BYTES});
+ printf("\tCORRELATION_DATA_BYTES\t\t= %3d bytes\n",$hdr{CORRELATION_DATA_BYTES});
+ printf("\tECHO_INTENSITY_DATA_BYTES\t= %3d bytes\n",$hdr{ECHO_INTENSITY_DATA_BYTES});
+ printf("\tPERCENT_GOOD_DATA_BYTES\t\t= %3d bytes\n",$hdr{PERCENT_GOOD_DATA_BYTES});
+ printf("\tBT_DATA_BYTES\t\t\t= %3d bytes\n",$hdr{BT_DATA_BYTES})
+ if ($dta->{BT_PRESENT});
+
+
+ print(" Coordinate System:\n");
+
+ printf("\tHEADING_ALIGNMENT_CORRECTION\t\t= %g deg\n",
+ $hdr{HEADING_ALIGNMENT_CORRECTION})
+ if defined($hdr{HEADING_ALIGNMENT_CORRECTION});
+ printf("\tHEADING_BIAS_CORRECTION\t\t= %g deg\n",
+ $hdr{HEADING_BIAS_CORRECTION})
+ if defined($hdr{HEADING_BIAS_CORRECTION});
+ print("\tFlags\t\t\t\t: ");
+ printf("BEAM_COORDINATES ") if ($hdr{BEAM_COORDINATES});
+ printf("INSTRUMENT_COORDINATES ") if ($hdr{INSTRUMENT_COORDINATES});
+ printf("SHIP_COORDINATES ") if ($hdr{SHIP_COORDINATES});
+ printf("EARTH_COORDINATES ") if ($hdr{EARTH_COORDINATES});
+ printf("PITCH_AND_ROLL_USED ") if ($hdr{PITCH_AND_ROLL_USED});
+ printf("BIN_MAPPING_ALLOWED ") if ($hdr{ALLOW_BIN_MAPPING});
+ print("\n");
+
+
+ if ($hdr{SPEED_OF_SOUND_CALCULATED}) {
+ print(" Speed-of-Sound Sensors Used:\n");
+ printf("\tPRESSURE_SENSOR_USED\n") if ($hdr{PRESSURE_SENSOR_USED});
+ printf("\tCOMPASS_USED\n") if ($hdr{COMPASS_USED});
+ printf("\tPITCH_SENSOR_USED\n") if ($hdr{PITCH_SENSOR_USED});
+ printf("\tROLL_SENSOR_USED\n") if ($hdr{ROLL_SENSOR_USED});
+ printf("\tCONDUCTIVITY_SENSOR_USED\n")
+ if ($hdr{CONDUCTIVITY_SENSOR_USED});
+ printf("\tTEMPERATURE_SENSOR_USED\n")
+ if ($hdr{TEMPERATURE_SENSOR_USED});
+ print("\n");
+ }
+
+
+ print(" Bin Setup:\n");
+ printf("\tN_BINS\t\t\t\t= %d\n", $hdr{N_BINS});
+ printf("\tBLANKING_DISTANCE\t\t= %g m\n", $hdr{BLANKING_DISTANCE});
+ printf("\tTRANSMIT_LAG_DISTANCE\t\t= %g m\n",
+ $hdr{TRANSMIT_LAG_DISTANCE});
+ printf("\tDISTANCE_TO_BIN1_CENTER\t\t= %g m\n",
+ $hdr{DISTANCE_TO_BIN1_CENTER});
+ printf("\tBIN_LENGTH\t\t\t= %g m\n", $hdr{BIN_LENGTH});
+ printf("\tTRANSMITTED_PULSE_LENGTH\t= %g m\n",
+ $hdr{TRANSMITTED_PULSE_LENGTH});
+ printf("\tRL_FIRST_BIN\t\t\t= %d\n", $hdr{RL_FIRST_BIN});
+ printf("\tRL_LAST_BIN\t\t\t= %d\n", $hdr{RL_LAST_BIN});
+
+
+ print(" Water-Track Setup:\n");
+ printf("\tPINGS_PER_ENSEMBLE\t\t= %d\n", $hdr{PINGS_PER_ENSEMBLE});
+ printf("\tTIME_BETWEEN_PINGS\t\t= %g s\n",$hdr{TIME_BETWEEN_PINGS});
+ printf("\tTRANSMIT_POWER\t\t\t= %d\n", $hdr{TRANSMIT_POWER});
+ printf("\tMIN_CORRELATION\t\t\t= %d\n", $hdr{MIN_CORRELATION});
+ printf("\tMIN_PERCENT_GOOD\t\t= %d %%\n", $hdr{MIN_PERCENT_GOOD});
+ printf("\tMAX_ERROR_VELOCITY\t\t= %g m/s\n",
+ $hdr{MAX_ERROR_VELOCITY});
+ printf("\tFALSE_TARGET_THRESHOLD\t\t= %d\n",
+ $hdr{FALSE_TARGET_THRESHOLD})
+ if defined($hdr{FALSE_TARGET_THRESHOLD});
+ printf("\tFlags\t\t\t\t: ");
+ printf("NARROW_BANDWIDTH ") if ($hdr{NARROW_BANDWIDTH});
+ printf("WIDE_BANDWIDTH ") if ($hdr{WIDE_BANDWIDTH});
+ printf("TRANSMIT_POWER_HIGH ") if ($hdr{TRANSMIT_POWER_HIGH});
+ printf("USE_3_BEAM_ON_LOW_CORR ") if ($hdr{USE_3_BEAM_ON_LOW_CORR});
+ print("\n");
+
+#----------------------------------------------------------------------
+
+ if ($hdr{NUMBER_OF_DATA_TYPES} == 7) {
+ print(" Bottom-Track Setup:\n");
+ printf("\tBT_MODE\t\t\t\t= %d\n", $hdr{BT_MODE});
+ printf("\tBT_PINGS_PER_ENSEMBLE\t\t= %d\n",
+ $hdr{BT_PINGS_PER_ENSEMBLE});
+ printf("\tBT_TIME_BEFORE_REACQUIRE\t= %g s\n",
+ $hdr{BT_TIME_BEFORE_REACQUIRE});
+ printf("\tBT_MIN_CORRELATION\t\t= %d\n",$hdr{BT_MIN_CORRELATION});
+ printf("\tBT_MIN_EVAL_AMPLITUDE\t\t= %d\n",
+ $hdr{BT_MIN_EVAL_AMPLITUDE});
+ printf("\tBT_MIN_PERCENT_GOOD\t\t= %d %%\n",
+ $hdr{BT_MIN_PERCENT_GOOD});
+ printf("\tBT_MAX_ERROR_VELOCITY\t\t= %g m/s\n",
+ $hdr{BT_MAX_ERROR_VELOCITY})
+ if defined($hdr{BT_MAX_ERROR_VELOCITY});
+ printf("\tBT_RL_MIN_SIZE\t\t\t= %g m\n", $hdr{BT_RL_MIN_SIZE});
+ printf("\tBT_RL_NEAR\t\t\t= %g m\n", $hdr{BT_RL_NEAR});
+ printf("\tBT_RL_FAR\t\t\t= %g m\n", $hdr{BT_RL_FAR});
+ printf("\tBT_MAX_TRACKING_DEPTH\t\t= %g m\n"
+ , $hdr{BT_MAX_TRACKING_DEPTH});
+ }
+}
+
+
new file mode 100755
--- /dev/null
+++ b/listW
@@ -0,0 +1,221 @@
+#!/usr/bin/perl
+#======================================================================
+# L I S T W
+# doc: Wed Mar 24 06:45:09 2004
+# dlm: Thu Jul 30 17:42:33 2009
+# (c) 2004 A.M. Thurnherr
+# uE-Info: 205 53 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# dump vertical velocities
+
+# NB: currently broken
+
+# HISTORY:
+# Mar 24, 2004: - created from [listens] and [mkprofile]
+# Mar 27, 2004: - added elapsed field
+# - floatized time
+# Apr 3, 2004: - cosmetics
+# Nov 8, 2005: - UNIXTIME => UNIX_TIME
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
+# Jul 30, 2009: - NaN => nan
+
+$0 =~ m{(.*)/[^/]+};
+require "$1/WorkhorseBinRead.pl";
+require "$1/WorkhorseCoords.pl";
+require "$1/WorkhorseUtils.pl";
+
+use Getopt::Std;
+
+$USAGE = "$0 @ARGV";
+die("Usage: $0 " .
+ "[-A)nts] " .
+ "[-F)ilter <script>] " .
+ "[bin -r)ange <bin|0,bin|*>] " .
+ "[-e)rr-vel <max|0.1>] [-c)orrelation <min|70>] " .
+ "[-S)alin <val|35>] [-t)emp <bias>] " .
+ "[output -f)ields <field[,...]> " .
+ "<RDI file>\n")
+ unless (&getopts("Ac:e:F:f:r:S:t:") && @ARGV == 1);
+
+$opt_e = 0.1 unless defined($opt_e); # defaults
+$opt_c = 70 unless defined($opt_c);
+$opt_S = 35 unless defined($opt_S);
+print(STDERR "WARNING: Using uncalibrated ADCP temperature!\n"),$opt_t = 0
+ unless defined($opt_t);
+
+require $opt_F if defined($opt_F); # load filter
+
+if ($opt_f) { # additional fields
+ @f = split(',',$opt_f);
+ foreach $f (@f) {
+ my($f) = $f; # copy it
+ $f =~ s{\[.*$}{}; # remove indices
+ $addFields .= " {$f}";
+ }
+}
+
+#----------------------------------------------------------------------
+
+print(STDERR "Reading $ARGV[0]..."); # read data
+readData($ARGV[0],\%dta);
+print(STDERR "done\n");
+
+if (defined($opt_r)) { # bin range
+ ($minb,$maxb) = split(',',$opt_r);
+ die("$0: can't decode -r $opt_r\n") unless defined($maxb);
+} else {
+ $minb = 0;
+ $maxb = $dta{N_BINS} - 1;
+}
+
+die("$ARGV[0]: not enough bins for choice of -r\n") # enough bins?
+ unless ($dta{N_BINS} >= $maxb);
+
+if ($dta{BEAM_COORDINATES}) { # coords used
+ $beamCoords = 1;
+} elsif (!$dta{EARTH_COORDINATES}) {
+ die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
+}
+
+#----------------------------------------------------------------------
+# Reference-Layer w (from [mkprofile])
+# - also sets W field when valid
+#----------------------------------------------------------------------
+
+sub ref_lr_w($)
+{
+ my($ens) = @_;
+ my($i,$n,$w);
+
+ for ($i=$minb; $i<=$maxb; $i++) {
+ next if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $opt_c ||
+ $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $opt_c ||
+ $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $opt_c ||
+ $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $opt_c);
+ if ($beamCoords) {
+ next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+ @v = velInstrumentToEarth(\%dta,$ens,
+ velBeamToInstrument(\%dta,
+ @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
+ } else {
+ next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+ @v = @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]};
+ }
+ next if (!defined($v[3]) || abs($v[3]) > $opt_e);
+
+ if (defined($v[2])) { # valid w
+ $dta{ENSEMBLE}[$ens]->{W}[$i] = $v[2];
+ $w += $v[2]; $n++;
+ }
+ }
+ return $n ? $w/$n : undef;
+}
+
+#----------------------------------------------------------------------
+
+print(STDERR "Generating profile by integrating w...");
+
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {
+ checkEnsemble(\%dta,$e); # sanity checks
+ filterEnsemble(\%dta,$e) # filter ensemble
+ if (defined($opt_F) &&
+ $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[0][0] > 0);
+
+ $dta{ENSEMBLE}[$e]->{REFW} = ref_lr_w($e);
+ next unless defined($dta{ENSEMBLE}[$e]->{REFW});
+
+ unless (defined($firstgood)) { # init profile
+ $firstgood = $lastgood = $e;
+ $depth = 0;
+ }
+
+ my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} - # time step since
+ $dta{ENSEMBLE}[$lastgood]->{UNIX_TIME}; # ... last good ens
+ if ($dt > 120) {
+ printf(STDERR "\nWARNING: %d-s gap too long, profile restarted at ensemble $e...",$dt);
+ $firstgood = $lastgood = $e;
+ $dt = $depth = 0;
+ }
+
+ $depth += $dta{ENSEMBLE}[$lastgood]->{REFW} * $dt # integrate depth
+ if ($dt > 0);
+ $dta{ENSEMBLE}[$e]->{DEPTH} = $depth;
+ $atbottom = $e, $maxdepth = $depth if ($depth > $maxdepth);
+
+ my($ss) = soundSpeed($opt_S, # sound-speed corr
+ $dta{ENSEMBLE}[$e]->{TEMPERATURE}-$opt_t,
+ $depth);
+ $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} =
+ $ss / $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
+ $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND} = $ss;
+ $dta{ENSEMBLE}[$e]->{REFW} *=
+ $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION};
+
+ $lastgood = $e;
+}
+
+printf(STDERR "done (max depth = %.1fm, depth at end of cast = %.1fm)\n",
+ $maxdepth,$depth);
+
+filterEnsembleStats() if defined($opt_F);
+
+#----------------------------------------------------------------------
+
+print(STDERR "Writing output...");
+
+if ($opt_A) {
+ print("#ANTS# [] $USAGE\n");
+ print("#ANTS#FIELDS# {ens} {unix_time} {time} {bin} {depth} {dz} {w} {ref_w} {dw} $addFields\n");
+ printf("#ANTS#PARAMS# maxdepth{$max_depth} bottom_time{%d}\n",
+ $dta{ENSEMBLE}[$atbottom]->{UNIX_TIME} -
+ $dta{ENSEMBLE}[$firstgood]->{UNIX_TIME});
+} else {
+ print("# ens-no time elapsed bin-no depth dz w ref-w dw $addFields\n");
+ print("#----------------------------------------------------------------------\n");
+}
+
+for ($e=$firstgood; $e<=$lastgood; $e++) {
+
+ for ($i=$minb; $i<=$maxb; $i++) { # dump valid
+ next unless defined($dta{ENSEMBLE}[$e]->{W}[$i]);
+
+ printf("%d %f %f %d %.1f %.1f %g %g %g ",
+ $e,$dta{ENSEMBLE}[$e]->{UNIX_TIME},
+ $dta{ENSEMBLE}[$e]->{UNIX_TIME} -
+ $dta{ENSEMBLE}[$firstgood]->{UNIX_TIME},$i,
+ $dta{ENSEMBLE}[$e]->{DEPTH} +
+ $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+ ($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),
+ $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+ ($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),
+ $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+ $dta{ENSEMBLE}[$e]->{W}[$i],
+ $dta{ENSEMBLE}[$e]->{REFW},
+ $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+ $dta{ENSEMBLE}[$e]->{W}[$i] - $dta{ENSEMBLE}[$e]->{REFW},
+ );
+
+ sub p($) { print(defined($_[0])?"$_[0] ":"nan "); }
+
+ if (defined(@f)) {
+ foreach $f (@f) {
+ my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
+ $fn = $f unless defined($fn);
+ p(eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));
+ }
+ }
+ print("\n");
+ }
+}
+
+print(STDERR "done\n");
+
+exit(0);
+
new file mode 100755
--- /dev/null
+++ b/meanProf
@@ -0,0 +1,319 @@
+#!/usr/bin/perl
+#======================================================================
+# M E A N P R O F
+# doc: Fri Feb 22 08:40:18 2008
+# dlm: Fri Feb 22 18:02:40 2008
+# (c) 2008 A.M. Thurnherr
+# uE-Info: 232 0 NIL 0 0 72 74 2 4 NIL ofnI
+#======================================================================
+
+# extract time-averaged mean profile from ADCP data
+
+# HISTORY:
+# Feb 22, 2008: - created from [listBins]
+
+# Soundspeed Correction:
+# - applied as described in the RDI coord-trans manual
+# - sound-speed variation over range is ignored (valid for small gradients)
+# => - same simple correction for all velocity components
+# - simple correction for cell depths
+
+require "getopts.pl";
+$0 =~ m{(.*/)[^/]+};
+require "$1RDI_BB_Read.pl";
+require "$1RDI_Coords.pl";
+require "$1RDI_Utils.pl";
+
+die("Usage: $0 [-r)ange <first_ens,last_ens>] " .
+ "[-Q)uiet (stats only)] " .
+ "[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
+ "[require -4)-beam solutions] [-d)iscard <beam#>] " .
+ "[-%)good <min>] " .
+ "[output -b)eam coordinates] " .
+ "[-M)agnetic <declination>] " .
+ "[-D)epth <depth>] " .
+ "<RDI file>\n")
+ unless (&Getopts("4bd:D:M:p:r:QS:") && @ARGV == 1);
+
+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);
+
+$ifn = $ARGV[0];
+
+($first_ens,$last_ens) = split(',',$opt_r)
+ if defined($opt_r);
+
+($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S)
+ if defined($opt_S);
+die("$0: Cannot do variable soundspeed correction (implementation restriction)\n")
+ if ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*');
+
+#----------------------------------------------------------------------
+# Read & Check Data, Transform Velocities
+#----------------------------------------------------------------------
+
+$P{RDI_file} = $ifn;
+$P{mag_decl} = $opt_M if defined($opt_M);
+
+print(STDERR "reading $ifn: ");
+readData($ifn,\%dta);
+printf(STDERR "%d complete ensembles.\n",scalar(@{$dta{ENSEMBLE}}));
+$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
+
+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")
+ if (!$dta{EARTH_COORDINATES});
+}
+
+for (my($b)=0; $b<$dta{N_BINS}; $b++) { # calc dz
+ $dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
+}
+
+$lastGoodBin = 0;
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities
+ next if (defined($first_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
+ $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e
+ unless defined($P{first_ens});
+ last if (defined($last_ens) &&
+ $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);
+ $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER};
+ $le = $e;
+
+ die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);
+ die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});
+ die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
+ if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
+
+ for (my($b)=0; $b<$dta{N_BINS}; $b++) {
+ if (defined($opt_d)) {
+ undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
+ undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
+ }
+ 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);
+
+ $sum_corr1[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0];
+ $sum_corr2[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][1];
+ $sum_corr3[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][2];
+ $sum_corr4[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][3];
+
+ $sum_amp1[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][0];
+ $sum_amp2[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][1];
+ $sum_amp3[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][2];
+ $sum_amp4[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][3];
+
+ 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);
+ $firstGoodEns = $e unless defined($firstGoodEns);
+ $lastGoodEns = $e;
+
+ $sum_u[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0];
+ $sum_v[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1];
+ $sum_w[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2];
+ $sum_e[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]
+ unless ($RDI_Coords::threeBeamFlag);
+
+ $sum_pcg1[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0];
+ $n_pcg1[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]);
+ $sum_pcg2[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1];
+ $n_pcg2[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]);
+ $sum_pcg3[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2];
+ $n_pcg3[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]);
+ $sum_pcg4[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3];
+ $n_pcg4[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]);
+ }
+}
+
+unless (defined($opt_r)) {
+ $fe = $firstGoodEns;
+ $le = $lastGoodEns;
+}
+$nEns = $le - $fe + 1;
+die("$0: insufficient data\n") if ($nEns < 2);
+$P{N_ensembles} = $nEns;
+
+$firstBin = 0;
+$lastBin = $lastGoodBin;
+
+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);
+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)
+ unless ($opt_b);
+
+#----------------------------------------------------------------------
+# Calculate Stddevs
+#----------------------------------------------------------------------
+
+for ($b=0; $b<=$lastGoodBin; $b++) {
+ $mean_corr1[$b] = $sum_corr1[$b] / $nEns; $mean_corr2[$b] = $sum_corr2[$b] / $nEns;
+ $mean_corr3[$b] = $sum_corr3[$b] / $nEns; $mean_corr4[$b] = $sum_corr4[$b] / $nEns;
+ $mean_amp1[$b] = $sum_amp1[$b] / $nEns; $mean_amp2[$b] = $sum_amp2[$b] / $nEns;
+ $mean_amp3[$b] = $sum_amp3[$b] / $nEns; $mean_amp4[$b] = $sum_amp4[$b] / $nEns;
+ $mean_pcg1[$b] = $sum_pcg1[$b] / $n_pcg1[$b]; $mean_pcg2[$b] = $sum_pcg2[$b] / $n_pcg2[$b];
+ $mean_pcg3[$b] = $sum_pcg3[$b] / $n_pcg3[$b]; $mean_pcg4[$b] = $sum_pcg4[$b] / $n_pcg4[$b];
+ $mean_u[$b] = $sum_u[$b] / $good_vels[$b]; $mean_v[$b] = $sum_v[$b] / $good_vels[$b];
+ $mean_w[$b] = $sum_w[$b] / $good_vels[$b];
+ $mean_e[$b] = $sum_e[$b] / ($good_vels[$b] - $three_beam[$b]);
+}
+
+for ($e=$fe; $e<=$le; $e++) {
+ for ($b=0; $b<=$lastGoodBin; $b++) {
+ $sumsq_corr1[$b] += ($mean_corr1[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0])**2;
+ $sumsq_corr2[$b] += ($mean_corr2[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][1])**2;
+ $sumsq_corr3[$b] += ($mean_corr3[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][2])**2;
+ $sumsq_corr4[$b] += ($mean_corr4[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][3])**2;
+
+ $sumsq_amp1[$b] += ($mean_amp1[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][0])**2;
+ $sumsq_amp2[$b] += ($mean_amp2[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][1])**2;
+ $sumsq_amp3[$b] += ($mean_amp3[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][2])**2;
+ $sumsq_amp4[$b] += ($mean_amp4[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][3])**2;
+
+ $sumsq_pcg1[$b] += ($mean_pcg1[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0])**2
+ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]);
+ $sumsq_pcg2[$b] += ($mean_pcg2[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1])**2
+ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]);
+ $sumsq_pcg3[$b] += ($mean_pcg3[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2])**2
+ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]);
+ $sumsq_pcg4[$b] += ($mean_pcg4[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3])**2
+ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]);
+
+ next unless ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]);
+
+ $sumsq_u[$b] += ($mean_u[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])**2;
+ $sumsq_v[$b] += ($mean_v[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1])**2;
+ $sumsq_w[$b] += ($mean_w[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2])**2;
+
+ next if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]);
+
+ $sumsq_e[$b] += ($mean_e[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])**2;
+ }
+}
+
+for ($b=0; $b<=$lastGoodBin; $b++) {
+ $var_corr1[$b] = $sumsq_corr1[$b] / ($nEns-1); $var_corr2[$b] = $sumsq_corr2[$b] / ($nEns-1);
+ $var_corr3[$b] = $sumsq_corr3[$b] / ($nEns-1); $var_corr4[$b] = $sumsq_corr4[$b] / ($nEns-1);
+ $var_amp1[$b] = $sumsq_amp1[$b] / ($nEns-1); $var_amp2[$b] = $sumsq_amp2[$b] / ($nEns-1);
+ $var_amp3[$b] = $sumsq_amp3[$b] / ($nEns-1); $var_amp4[$b] = $sumsq_amp4[$b] / ($nEns-1);
+ $var_pcg1[$b] = $sumsq_pcg1[$b] / ($n_pcg1[$b]-1) if ($n_pcg1[$b] > 1);
+ $var_pcg2[$b] = $sumsq_pcg2[$b] / ($n_pcg2[$b]-1) if ($n_pcg2[$b] > 1);
+ $var_pcg3[$b] = $sumsq_pcg3[$b] / ($n_pcg3[$b]-1) if ($n_pcg3[$b] > 1);
+ $var_pcg4[$b] = $sumsq_pcg4[$b] / ($n_pcg4[$b]-1) if ($n_pcg4[$b] > 1);
+ next unless ($good_vels[$b] > 1);
+ $var_u[$b] = $sumsq_u[$b] / ($good_vels[$b]-1);
+ $var_v[$b] = $sumsq_v[$b] / ($good_vels[$b]-1);
+ $var_w[$b] = $sumsq_w[$b] / ($good_vels[$b]-1);
+ next unless ($good_vels[$b] - $three_beam[$b] > 1);
+ $var_e[$b] = $sumsq_e[$b] / ($good_vels[$b] - $three_beam[$b] - 1);
+}
+
+#----------------------------------------------------------------------
+# Calculate Beam Statistics
+#----------------------------------------------------------------------
+
+#----------------------------------------------------------------------
+# Produce Output
+#----------------------------------------------------------------------
+
+unless ($opt_Q) {
+ my($ssCorr) = defined($opt_S)
+ ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth)
+ : 1;
+
+ print("#ANTS#PARAMS# ");
+ foreach my $k (keys(%P)) {
+ print("$k\{$P{$k}\} ");
+ }
+ printf("soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!');
+ print("\n");
+
+ print("#ANTS#FIELDS# " .
+ "{bin} {dz} " .
+ (defined($opt_D) ? "{depth} " : "") .
+ ($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") .
+ ($opt_b ? "{sig_v1} {sig_v2} {sig_v3} {sig_v4} " : "{sig_u} {sig_v} {sig_w} {sig_err_vel} ") .
+ "{corr1} {corr2} {corr3} {corr4} " .
+ "{sig_corr1} {sig_corr2} {sig_corr3} {sig_corr4} " .
+ "{amp1} {amp2} {amp3} {amp4} " .
+ "{sig_amp1} {sig_amp2} {sig_amp3} {sig_amp4} " .
+ "{pcg1} {pcg2} {pcg3} {pcg4} " .
+ "{sig_pcg1} {sig_pcg2} {sig_pcg3} {sig_pcg4}" .
+ "\n"
+ );
+
+ for ($b=$firstBin; $b<=$lastBin; $b++) {
+ printf("%d %.1f ",$b+1,$dz[$b]*$ssCorr);
+ printf("%.1f ",$opt_D - $dz[$b]*$ssCorr)
+ if defined($opt_D);
+
+ printf("%s ",defined($mean_u[$b]) ? $mean_u[$b] : nan);
+ printf("%s ",defined($mean_v[$b]) ? $mean_v[$b] : nan);
+ printf("%s ",defined($mean_w[$b]) ? $mean_w[$b] : nan);
+ printf("%s ",defined($mean_e[$b]) ? $mean_e[$b] : nan);
+
+ printf("%s ",defined($var_u[$b]) ? sqrt($var_u[$b]) : nan);
+ printf("%s ",defined($var_v[$b]) ? sqrt($var_v[$b]) : nan);
+ printf("%s ",defined($var_w[$b]) ? sqrt($var_w[$b]) : nan);
+ printf("%s ",defined($var_e[$b]) ? sqrt($var_e[$b]) : nan);
+
+ printf("%g %g %g %g ",$mean_corr1[$b],$mean_corr2[$b],
+ $mean_corr3[$b],$mean_corr4[$b]);
+ printf("%g %g %g %g ",sqrt($var_corr1[$b]),sqrt($var_corr2[$b]),
+ sqrt($var_corr3[$b]),sqrt($var_corr4[$b]));
+
+ printf("%g %g %g %g ",$mean_amp1[$b],$mean_amp2[$b],
+ $mean_amp3[$b],$mean_amp4[$b]);
+ printf("%g %g %g %g ",sqrt($var_amp1[$b]),sqrt($var_amp2[$b]),
+ sqrt($var_amp3[$b]),sqrt($var_amp4[$b]));
+
+ if ($good_vels[$b] > 0) {
+ printf("%g %g %g %g ",$mean_pcg1[$b],$mean_pcg2[$b],
+ $mean_pcg3[$b],$mean_pcg4[$b]);
+ printf("%g %g %g %g\n",sqrt($var_pcg1[$b]),sqrt($var_pcg2[$b]),
+ sqrt($var_pcg3[$b]),sqrt($var_pcg4[$b]));
+ } else {
+ print("nan nan nan nan ");
+ print("nan nan nan nan\n");
+ }
+ }
+}
+
+exit(0);
new file mode 100755
--- /dev/null
+++ b/mkProfile
@@ -0,0 +1,808 @@
+#!/usr/bin/perl
+#======================================================================
+# M K P R O F I L E
+# doc: Sun Jan 19 18:55:26 2003
+# dlm: Sun May 23 16:34:02 2010
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 788 36 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# Make an LADCP Profile by Integrating W (similar to Firing's scan*).
+
+# HISTORY:
+# Jan 19, 2003: - written in order to test the RDI libs
+# Jan 20, 2003: - added ensemble number
+# Jan 21, 2003: - added horizontal integration
+# Jan 22, 2003: - corrected magnetic declination
+# Jan 23, 2003: - added -F)ilter
+# Jan 24, 2003: - added more %PARAMs; started integration from 1st bin
+# - added -g, -f, battery status
+# Jan 25, 2003: - added more %PARAMs
+# Feb 1, 2003: - BUG: bottom-track quality checking was bad
+# Feb 8, 2003: - allowed for array-indices on -f
+# Feb 9, 2003: - added 50% goodvelbin
+# - removed unknown-field err on -f to allow -f W
+# Feb 10, 2003: - changed initialization depth to 0m
+# - changed %bottom_depth to %max_depth
+# Feb 11, 2003: - changed sign of magnetic declination
+# Feb 12, 2003: - corrected BT-range scaling
+# Feb 14, 2003: - added %pinging_hours, %min_range
+# - removed magnetic declination from default
+# Feb 26, 2004: - added earth coordinates
+# Mar 3, 2004: - removed requirement for -M on !-Q
+# - corrected range-stats on earth coordinates
+# Mar 4, 2004: - added number of ensebles to output
+# Mar 11, 2004: - BUG: rename ACD -> ADC
+# Mar 12, 2004: - added %bottom_xmit_{current|voltage}
+# Mar 16, 2004: - BUG: on -M u/v/x/y were wrong
+# Mar 17, 2004: - added error estimates on u/v/x/y
+# - removed battery stuff (has to be done btw casts)
+# Mar 18, 2004: - totally re-did u/v integration
+# Mar 19, 2004: - re-designed u/v uncertainty estimation
+# Mar 28, 2004: - added MEAN_CORRELATION, MEAN_ECHO_AMPLITUDE
+# Sep 15, 2005: - changed BinRead library name
+# - made max gap length variable
+# Sep 16, 2005: - re-did u,v,w uncertainties
+# Nov 8, 2005: - UNIXTIME => UNIX_TIME
+# - added unix_time, secno, z_BT to default output
+# Dec 1, 2005: - moved profile-building code to [RDI_utils.pl]
+# - changed -f syntax to allow name=FIELD
+# - added %bin1_dist, %bin_length
+# Dec 8, 2005: - remove spaces from -f argument to allow multiline
+# definitions in Makefiles
+# Nov 13, 2006: - BUG: end-of-cast depth had not been reported correctly
+# - cosmetics
+# Nov 30, 2007: - adapted to 3-beam solutions
+# Dec 11, 2007: - adapted to earlier modifications (Sep 2007) of
+# [RDI_BB_Read.pl]
+# Dec 14, 2007: - replaced z by depth
+# Dec 17, 2007: - BUG: downcast flag was set incorrectly
+# Jan 24, 2008: - rotation had been output as degrees/s; to make it more
+# consistent with pitch/roll, I changed it to simple degrees
+# - added net rotations [deployment]/down/up/[recovery]
+# Apr 9, 2008: - added profile -B)ottom depth
+# - BUG: depth of first bin was reported as beginning of cast
+# Oct 24, 2008: - added RANGE and RANGE_BINS fields
+# Mar 18, 2009: - BUG: pitch/roll calculation had typo
+# - calc pitch/roll separately for down-/upcasts
+# - removed approximations in pitch/roll calcs
+# Jul 30, 2009: - typo '<' removed from output
+# - NaN => nan
+
+# NOTES:
+# - the battery values are based on transmission voltages (different
+# from battery voltages) and reported without units (raw 8-bit a2d
+# values)
+# - -B with the CTD max depth can be used to linearly scale the depths;
+# even so, the profile can have negative depths, in particular when
+# the CTD is sent to a shallow depth first and then returned to the surface
+# before beginning the cast
+# - in one case that I looked at (Anslope ][, cast 82), there are large
+# depth errors, even when -B is used
+# - this utility works only approximately for uplookers (profile is
+# roughly ok, but apparently contaminated by surface reflection,
+# but stats are not ok; e.g. NBP0402 037U.prof)
+
+$0 =~ m{(.*)/[^/]+};
+require "$1/RDI_BB_Read.pl";
+require "$1/RDI_Coords.pl";
+require "$1/RDI_Utils.pl";
+require "getopts.pl";
+
+$USAGE = "$0 @ARGV";
+die("Usage: $0 " .
+ "[-A)nts] [-Q)uiet] [-F)ilter <script>] " .
+ "[-s)uppress checkensemble()] " .
+ "[require -4)-beam solutions] " .
+ "[-r)ef-layer <bin|1,bin|6>] [-n) vels <min|2>] " .
+ "[-e)rr-vel <max|0.1>] [-c)orrelation <min>] " .
+ "[-m)ax <gap>] " .
+ "[-d)rift <dx,dy>] [-g)ps <start lat,lon/end lat,lon>] " .
+ "[output -f)ields <field[,...]> " .
+ "[-M)agnetic <declination>] [profile -B)ottom <depth>] " .
+ "<RDI file>\n")
+ unless (&Getopts("4AB:F:M:Qd:r:n:e:c:g:f:m:s") && @ARGV == 1);
+
+die("$0: -Q and -A are mutually exclusive\n")
+ if ($opt_Q && $opt_A);
+
+$RDI_Coords::minValidVels = 4 if ($opt_4); # no 3-beam solutions
+
+require $opt_F if defined($opt_F); # load filter
+
+$opt_r = "1,6" unless defined($opt_r); # defaults
+$opt_n = 2 unless defined($opt_n);
+$opt_e = 0.1 unless defined($opt_e);
+$opt_c = 70 unless defined($opt_c);
+$opt_m = 120 unless defined($opt_m);
+
+($minb,$maxb) = split(',',$opt_r); # reference layer
+die("$0: can't decode -r $opt_r\n") unless defined($maxb);
+
+if ($opt_g) { # GPS info
+ ($s_lat,$s_lon,$e_lat,$e_lon) = gps_to_deg($opt_g);
+ $lat = $s_lat/2 + $e_lat/2;
+ $lon = $s_lon/2 + $e_lon/2;
+ $ddx = dist($lat,$s_lon,$lat,$e_lon);
+ $ddy = dist($s_lat,$lon,$e_lat,$lon);
+}
+
+if ($opt_d) { # ship drift
+ ($ddx,$ddy) = split(',',$opt_d);
+ die("$0: can't decode -d $opt_d\n") unless defined($ddy);
+}
+
+print(STDERR "Reading $ARGV[0]..."); # read data
+readData($ARGV[0],\%dta);
+print(STDERR "done\n");
+
+die("$ARGV[0]: not enough bins for choice of -r\n") # enough bins?
+ unless ($dta{N_BINS} >= $maxb);
+if ($dta{BEAM_COORDINATES}) { # coords used
+ $beamCoords = 1;
+} elsif (!$dta{EARTH_COORDINATES}) {
+ die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
+}
+if (defined($opt_M)) { # magnetic declination
+ $dta{HEADING_BIAS} = -1*$opt_M;
+} else {
+ $dta{HEADING_BIAS} = 0;
+}
+
+ensure_BT_RANGE(\%dta); # calc if missing
+
+if ($opt_f) { # additional fields
+ @f = split(',',$opt_f);
+ foreach $f (@f) {
+ $f =~ s/\s//g; # remove spaces
+ @def = split('=',$f);
+ if (@def == 2) { # name=field
+ $addFields .= " {$def[0]}";
+ $f = $def[1];
+ } else { # field
+ $addFields .= " {$f}";
+ }
+ }
+# print(STDERR "addFields = $addFields\n");
+# print(STDERR "\@f = @f\n");
+}
+
+#======================================================================
+# Misc funs used to decode options
+#======================================================================
+
+sub dist($$$$) # distance
+{
+ my($lat1,$lon1,$lat2,$lon2) = @_;
+ my($a) = 6378139; # Earth's radius
+
+ $lat1 = rad($lat1); $lon1 = rad($lon1);
+ $lat2 = rad($lat2); $lon2 = rad($lon2);
+ $ct1 = cos($lat1); $st1 = sin($lat1); $cp1 = cos($lon1);
+ $sp1 = sin($lon1); $ct2 = cos($lat2); $st2 = sin($lat2);
+ $cp2 = cos($lon2); $sp2 = sin($lon2);
+ $cos = $ct1*$cp1*$ct2*$cp2 + $ct1*$sp1*$ct2*$sp2 + $st1*$st2;
+ $cos = 1 if ($cos > 1);
+ $cos = -1 if ($cos < -1);
+ return $a * acos($cos);
+}
+
+sub deg_to_dec($) # parse degrees
+{
+ my($deg,$min) = split(':',$_[0]);
+ return $deg + $min/60;
+}
+
+sub gps_to_deg($) # decode lat/lon
+{
+ my($start,$end) = split('/',$_[0]);
+ my($sa,$so,$ea,$eo);
+
+ my($lat,$lon) = split(',',$start);
+ if ($lat =~ m{N$}) { $sa = deg_to_dec($`); }
+ elsif ($lat =~ m{S$}) { $sa = -deg_to_dec($`); }
+ else { $sa = $lat; }
+ if ($lon =~ m{E$}) { $so = deg_to_dec($`); }
+ elsif ($lon =~ m{W$}) { $so = -deg_to_dec($`); }
+ else { $so = $lon; }
+
+ my($lat,$lon) = split(',',$end);
+ if ($lat =~ m{N$}) { $ea = deg_to_dec($`); }
+ elsif ($lat =~ m{S$}) { $ea = -deg_to_dec($`); }
+ else { $ea = $lat; }
+ if ($lon =~ m{E$}) { $eo = deg_to_dec($`); }
+ elsif ($lon =~ m{W$}) { $eo = -deg_to_dec($`); }
+ else { $eo = $lon; }
+
+ return ($sa,$so,$ea,$eo);
+}
+
+#======================================================================
+# Step 1: Integrate w & determine water depth
+#======================================================================
+
+($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz) =
+ mk_prof(\%dta,!$opt_s,$opt_F,$minb,$maxb,$opt_c,$opt_e,$opt_m);
+
+die("$ARGV[0]: no good ensembles found\n")
+ unless defined($firstgood);
+
+if (defined($opt_B)) { # scale Z
+ my($zscale) = $opt_B / ($dta{ENSEMBLE}[$atbottom]->{DEPTH} -# downcast
+ $dta{ENSEMBLE}[$firstgood]->{DEPTH});
+# printf(STDERR "scaling downcast depths by %.2f\n",$zscale);
+ for (my($e)=$firstgood; $e<$atbottom; $e++) {
+ next unless defined($dta{ENSEMBLE}[$e]->{DEPTH});
+ $dta{ENSEMBLE}[$e]->{DEPTH} =
+ $dta{ENSEMBLE}[$firstgood]->{DEPTH} + $zscale *
+ ($dta{ENSEMBLE}[$e]->{DEPTH}-$dta{ENSEMBLE}[$firstgood]->{DEPTH});
+ }
+
+ $zscale = $opt_B / ($dta{ENSEMBLE}[$atbottom]->{DEPTH} - # upcast
+ $dta{ENSEMBLE}[$lastgood]->{DEPTH});
+# printf(STDERR "scaling upcast depths by %.2f\n",$zscale);
+ for (my($e)=$atbottom; $e<=$lastgood; $e++) {
+ next unless defined($dta{ENSEMBLE}[$e]->{DEPTH});
+ $dta{ENSEMBLE}[$e]->{DEPTH} =
+ $dta{ENSEMBLE}[$firstgood]->{DEPTH} + $zscale *
+ ($dta{ENSEMBLE}[$e]->{DEPTH}-$dta{ENSEMBLE}[$lastgood]->{DEPTH});
+ }
+}
+
+($water_depth,$sig_wd) = # sea bed
+ find_seabed(\%dta,$atbottom,$beamCoords);
+
+#======================================================================
+# Step 1a: determine alternate Z by using mean/sigma of w in gaps
+#======================================================================
+
+# This does not make much sense for w, because w is always very close
+# to zero. It might make sense for u and v, though, and it would
+# be more consistent with the way the displacement uncertainties are
+# calculated. However, the way the profiles are calculated at the
+# moment (using the last valid velocity across the gap) is probably
+# closer to the truth in most cases.
+
+#$dta{ENSEMBLE}[$firstgood]->{ALT_Z} = 0;
+#$dta{ENSEMBLE}[$firstgood]->{ALT_Z_ERR} = 0;
+#my($sumVar);
+#for ($e=$firstgood+1; $e<=$lastgood; $e++) {
+# my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} -
+# $dta{ENSEMBLE}[$e-1]->{UNIX_TIME};
+# $dta{ENSEMBLE}[$e]->{ALT_Z} =
+# $dta{ENSEMBLE}[$e-1]->{ALT_Z} +
+# $dt * (defined($dta{ENSEMBLE}[$e-1]->{W}) ?
+# $dta{ENSEMBLE}[$e-1]->{W} : $meanW);
+# $sumVar += defined($dta{ENSEMBLE}[$e-1]->{W}) ?
+# ($dta{ENSEMBLE}[$e-1]->{W_ERR} * $dt)**2 : ($dt**2)*$varW;
+# $dta{ENSEMBLE}[$e]->{ALT_Z_ERR} = sqrt($sumVar);
+#}
+
+#======================================================================
+# Step 2: Integrate u & v
+#======================================================================
+
+sub ref_lr_uv($$$) # calc ref-level u/v
+{
+ my($ens,$z,$water_depth) = @_;
+ my($i,$n,@v,@goodU,@goodV);
+
+ $water_depth = 99999 unless defined($water_depth);
+
+ for ($i=$minb; $i<=$maxb; $i++) {
+ next if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $opt_c ||
+ $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $opt_c ||
+ $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $opt_c ||
+ $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $opt_c);
+ if ($beamCoords) {
+ next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+ @v = velInstrumentToEarth(\%dta,$ens,
+ velBeamToInstrument(\%dta,
+ @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
+ } else {
+ next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||
+ $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+ @v = velApplyHdgBias(\%dta,$ens,
+ @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]});
+ }
+ next if (!defined($v[3]) || abs($v[3]) > $opt_e);
+
+# Martin's BT routines show strong shear just above sea bed
+# => skip lowest 20m.
+ if (defined($v[0])) { # valid u,v
+ if ($dta{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) {
+ if ($z - $dta{DISTANCE_TO_BIN1_CENTER}
+ - $i*$dta{BIN_LENGTH} > 0) {
+ push(@goodU,$v[0]); push(@goodV,$v[1]);
+ $dta{ENSEMBLE}[$ens]->{U} += $v[0];
+ $dta{ENSEMBLE}[$ens]->{V} += $v[1];
+ $n++;
+ }
+ } else {
+ if ($z + $dta{DISTANCE_TO_BIN1_CENTER}
+ + $i*$dta{BIN_LENGTH} < $water_depth-20) {
+ push(@goodU,$v[0]); push(@goodV,$v[1]);
+ $dta{ENSEMBLE}[$ens]->{U} += $v[0];
+ $dta{ENSEMBLE}[$ens]->{V} += $v[1];
+ $n++;
+ }
+ }
+ }
+ }
+
+ if ($n >= 2) {
+ my(@sumsq) = (0,0);
+ $dta{ENSEMBLE}[$ens]->{U} /= $n;
+ $dta{ENSEMBLE}[$ens]->{V} /= $n;
+ for ($i=0; $i<$n; $i++) {
+ $sumsq[0] += ($dta{ENSEMBLE}[$ens]->{U}-$goodU[$i])**2;
+ $sumsq[1] += ($dta{ENSEMBLE}[$ens]->{V}-$goodV[$i])**2;
+ }
+ $dta{ENSEMBLE}[$ens]->{U_ERR} = sqrt($sumsq[0])/($n-1);
+ $dta{ENSEMBLE}[$ens]->{V_ERR} = sqrt($sumsq[1])/($n-1);
+ } else {
+ $dta{ENSEMBLE}[$ens]->{U} = undef;
+ $dta{ENSEMBLE}[$ens]->{V} = undef;
+ }
+}
+
+#----------------------------------------------------------------------
+
+($x,$y) = (0,0); # init
+
+$dta{ENSEMBLE}[$firstgood]->{X} = $dta{ENSEMBLE}[$firstgood]->{X_ERR} = 0;
+$dta{ENSEMBLE}[$firstgood]->{Y} = $dta{ENSEMBLE}[$firstgood]->{Y_ERR} = 0;
+$prevgood = $firstgood;
+
+for ($e=$firstgood+1; defined($opt_M)&&$e<=$lastgood; $e++) {
+
+ #--------------------------------------------------
+ # within profile: both $firstgood and $prevgood set
+ #--------------------------------------------------
+
+ ref_lr_uv($e,$dta{ENSEMBLE}[$e]->{DEPTH},$water_depth) # instrument vel
+ if (defined($dta{ENSEMBLE}[$e]->{W}));
+
+ if (!defined($dta{ENSEMBLE}[$e]->{U})) { # gap
+ $uv_gap_time += $dta{ENSEMBLE}[$e]->{UNIX_TIME} -
+ $dta{ENSEMBLE}[$e-1]->{UNIX_TIME};
+ next;
+ }
+
+ my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} - # time step since
+ $dta{ENSEMBLE}[$prevgood]->{UNIX_TIME}; # ...last good ens
+
+ #-----------------------------------
+ # The current ensemble has valid u/v
+ #-----------------------------------
+
+ $x -= $dta{ENSEMBLE}[$prevgood]->{U} * $dt; # integrate
+ $xErr += ($dta{ENSEMBLE}[$prevgood]->{U_ERR} * $dt)**2;
+ $dta{ENSEMBLE}[$e]->{X} = $x;
+ $dta{ENSEMBLE}[$e]->{X_ERR} = sqrt($xErr);
+
+ $y -= $dta{ENSEMBLE}[$prevgood]->{V} * $dt;
+ $yErr += ($dta{ENSEMBLE}[$prevgood]->{V_ERR} * $dt)**2;
+ $dta{ENSEMBLE}[$e]->{Y} = $y;
+ $dta{ENSEMBLE}[$e]->{Y_ERR} = sqrt($yErr);
+
+ $prevgood = $e;
+}
+
+unless (defined($dta{ENSEMBLE}[$lastgood]->{X})) { # last is bad in u/v
+ my($dt) = $dta{ENSEMBLE}[$lastgood]->{UNIX_TIME} - # time step since
+ $dta{ENSEMBLE}[$prevgood]->{UNIX_TIME}; # ...last good ens
+
+ $x -= $dta{ENSEMBLE}[$prevgood]->{U} * $dt; # integrate
+ $xErr += ($dta{ENSEMBLE}[$prevgood]->{U_ERR} * $dt)**2;
+ $dta{ENSEMBLE}[$lastgood]->{X} = $x;
+ $dta{ENSEMBLE}[$lastgood]->{X_ERR} = sqrt($xErr);
+
+ $y -= $dta{ENSEMBLE}[$prevgood]->{V} * $dt;
+ $yErr += ($dta{ENSEMBLE}[$prevgood]->{V_ERR} * $dt)**2;
+ $dta{ENSEMBLE}[$lastgood]->{Y} = $y;
+ $dta{ENSEMBLE}[$lastgood]->{Y_ERR} = sqrt($yErr);
+}
+
+$firstgood++ if ($firstgood == 0); # centered diff
+$lastgood-- if ($lastgood == $#{$dta{ENSEMBLE}}); # in step 6
+
+#======================================================================
+# Step 3: Calculate Uncertainties
+#======================================================================
+
+# Time series of W_ERR indicate that errors are very large near the
+# surface and near the sea bed, perhaps because of reflections.
+# A reasonable estimate for typical uncertainty is therefore the mode
+# of the std errors.
+
+my(@histUErr,@histVErr,@histWErr);
+my($histRez) = 1e-4;
+
+for ($e=$firstgood; $e<=$lastgood; $e++) {
+ $histWErr[int($dta{ENSEMBLE}[$e]->{W_ERR}/$histRez+0.5)]++
+ if defined($dta{ENSEMBLE}[$e]->{W_ERR});
+ $histUErr[int($dta{ENSEMBLE}[$e]->{U_ERR}/$histRez+0.5)]++
+ if defined($dta{ENSEMBLE}[$e]->{U_ERR});
+ $histVErr[int($dta{ENSEMBLE}[$e]->{V_ERR}/$histRez+0.5)]++
+ if defined($dta{ENSEMBLE}[$e]->{V_ERR});
+}
+
+my($max) = 0; my($mode);
+for (my($i)=0; $i<=$#histWErr; $i++) {
+ next if ($histWErr[$i] < $max);
+ $max = $histWErr[$i]; $mode = $i;
+}
+$wErr = $mode * $histRez if defined($mode);
+
+$max = 0; $mode = undef;
+for (my($i)=0; $i<=$#histUErr; $i++) {
+ next if ($histUErr[$i] < $max);
+ $max = $histUErr[$i]; $mode = $i;
+}
+$uErr = $mode * $histRez if defined($mode);
+
+$max = 0; $mode = undef;
+for (my($i)=0; $i<=$#histVErr; $i++) {
+ next if ($histVErr[$i] < $max);
+ $max = $histVErr[$i]; $mode = $i;
+}
+$vErr = $mode * $histRez if defined($mode);
+
+#print(STDERR "u: mu = $meanU / sigma = $uErr\n");
+#print(STDERR "v: mu = $meanV / sigma = $vErr\n");
+#print(STDERR "w: mu = $meanW / sigma = $wErr\n");
+
+if (defined($opt_M)) { # displacement errors
+ $x_err = $uErr * $uv_gap_time + $dta{ENSEMBLE}[$lastgood]->{X_ERR};
+ $y_err = $vErr * $uv_gap_time + $dta{ENSEMBLE}[$lastgood]->{Y_ERR};
+}
+$z_err = $wErr * $w_gap_time + $dta{ENSEMBLE}[$lastgood]->{DEPTH_ERR};
+
+#printf(STDERR "x_err = $dta{ENSEMBLE}[$lastgood]->{X_ERR} + %g\n",
+# $uErr * $uv_gap_time);
+#printf(STDERR "y_err = $dta{ENSEMBLE}[$lastgood]->{Y_ERR} + %g\n",
+# $vErr * $uv_gap_time);
+#printf(STDERR "z_err = $dta{ENSEMBLE}[$lastgood]->{DEPTH_ERR} + %g\n",
+# $wErr * $w_gap_time);
+
+#======================================================================
+# Step 4: Calculate Beam Range Stats
+#======================================================================
+
+my($min_good_bins) = 999;
+my($worst_beam);
+
+sub count_good_vels($) # count good vels
+{
+ my($ens) = @_;
+ my($good) = -1; my($this_worst_beam);
+
+ if ($beamCoords) {
+ for (my($i)=0; $i<$dta{N_BINS}; $i++) {
+ for (my($b)=0; $b<4; $b++) {
+ $good=$i,$this_worst_beam=$b,$nVels[$i][$b]++
+ if defined($dta{ENSEMBLE}[$ens]->{VELOCITY}[$i][$b]);
+ }
+ }
+ } else {
+ for (my($i)=0; $i<$dta{N_BINS}; $i++) {
+ for (my($b)=0; $b<4; $b++) {
+ $good=$i,$this_worst_beam=$b,$nVels[$i][$b]++
+ if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][$b] >=
+ $dta{MIN_CORRELATION});
+ }
+ }
+ }
+ $min_good_ens=$ens, $min_good_bins=$good, $worst_beam=$this_worst_beam
+ if ((!defined($water_depth) ||
+ $dta{ENSEMBLE}[$ens]->{DEPTH} < $water_depth-200)
+ && $good >= 0 && $good < $min_good_bins);
+}
+
+#----------------------------------------------------------------------
+
+for ($e=$firstgood; $e<=$lastgood; $e++) { # range
+ my($i);
+ for ($i=0; $i<$dta{N_BINS}; $i++) {
+ last if (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][0]) +
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][1]) +
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][2]) +
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][3]) < 3);
+ }
+ $dta{ENSEMBLE}[$e]->{RANGE_BINS} = $i;
+ $dta{ENSEMBLE}[$e]->{RANGE} =
+ $dta{DISTANCE_TO_BIN1_CENTER} + $i * $dta{BIN_LENGTH};
+}
+
+for ($e=$firstgood; $e<=$lastgood; $e++) { # mean corr/amp
+ $sumcor = $sumamp = $ndata = 0;
+ for (my($i)=0; $i<$dta{N_BINS}; $i++) {
+ for (my($b)=0; $b<4; $b++) {
+ next unless ($dta{ENSEMBLE}[$e]->{CORRELATION}[$i][$b]);
+ $sumcor += $dta{ENSEMBLE}[$e]->{CORRELATION}[$i][$b];
+ $sumamp += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b];
+ $ndata++;
+ }
+ }
+ $dta{ENSEMBLE}[$e]->{MEAN_CORRELATION} = $sumcor/$ndata;
+ $dta{ENSEMBLE}[$e]->{MEAN_ECHO_AMPLITUDE} = $sumamp/$ndata;
+}
+
+for ($e=$firstgood+50; $e<=$lastgood-50; $e++) { # range stats
+ count_good_vels($e);
+}
+for ($i=0; $i<$dta{N_BINS}; $i++) {
+ for ($b=0; $b<4; $b++) {
+ $maxVels = $nVels[$i][$b] unless ($maxVels > $nVels[$i][$b]);
+ }
+}
+for ($i=0; $i<$dta{N_BINS}; $i++) {
+ for ($b=0; $b<4; $b++) {
+ $gb[$b] = $i if ($nVels[$i][$b] >= 0.8*$maxVels);
+ }
+}
+$gb = ($gb[0]+$gb[1]+$gb[2]+$gb[3]) / 4;
+
+#======================================================================
+# Step 5: Remove Ship Drift (probably not useful)
+#======================================================================
+
+if (defined($opt_M) && defined($ddx)) { # remove barotropic
+ $du = $ddx / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};# mean drift vel
+ $dv = $ddy / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
+ $iu = $dta{ENSEMBLE}[$lastgood]->{X} / # mean obs vel
+ $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
+ $iv = $dta{ENSEMBLE}[$lastgood]->{Y} /
+ $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
+
+ for ($e=$firstgood; $e<=$lastgood; $e++) {
+ next unless (defined($dta{ENSEMBLE}[$e]->{X}) &&
+ defined($dta{ENSEMBLE}[$e]->{Y}));
+ $dta{ENSEMBLE}[$e]->{U} -= $du;
+ $dta{ENSEMBLE}[$e]->{V} -= $dv;
+ $dta{ENSEMBLE}[$e]->{X} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($du-$iu);
+ $dta{ENSEMBLE}[$e]->{Y} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($dv-$iv);
+ }
+}
+
+#======================================================================
+# Step 6: Pitch, Roll, Rotation
+#======================================================================
+
+my($prrms,$dnprrms,$upprrms) = (0,0,0);
+my($rotrms,$prerot,$dnrot,$uprot,$postrot) = (0,0,0,0,0);
+
+sub rot($)
+{
+ my($e) = @_;
+ my($rot) = $dta{ENSEMBLE}[$e]->{HEADING} -
+ $dta{ENSEMBLE}[$e-1]->{HEADING};
+ $rot -= 360 if ($rot > 180);
+ $rot += 360 if ($rot < -180);
+ return $rot;
+}
+
+for ($e=1; $e<$firstgood; $e++) { # pre-deployment
+ $prerot += rot($e);
+}
+
+for (; $e<= $atbottom; $e++) { # downcast
+ $dta{ENSEMBLE}[$e]->{PITCHROLL} =
+ &angle_from_vertical($dta{ENSEMBLE}[$e]->{PITCH},
+ $dta{ENSEMBLE}[$e]->{ROLL});
+ $prrms += $dta{ENSEMBLE}[$e]->{PITCHROLL}**2;
+
+ $dta{ENSEMBLE}[$e]->{ROTATION} = rot($e);
+ $dnrot += $dta{ENSEMBLE}[$e]->{ROTATION};
+ $rotrms += $dta{ENSEMBLE}[$e]->{ROTATION}**2;
+}
+$dnprrms = $prrms;
+
+for (; $e<=$lastgood; $e++) { # upcast
+ $dta{ENSEMBLE}[$e]->{PITCHROLL} =
+ &angle_from_vertical($dta{ENSEMBLE}[$e]->{PITCH},
+ $dta{ENSEMBLE}[$e]->{ROLL});
+ $prrms += $dta{ENSEMBLE}[$e]->{PITCHROLL}**2;
+
+ $dta{ENSEMBLE}[$e]->{ROTATION} = rot($e);
+ $uprot += $dta{ENSEMBLE}[$e]->{ROTATION};
+ $rotrms += $dta{ENSEMBLE}[$e]->{ROTATION}**2;
+}
+$upprrms = $prrms - $dnprrms;
+
+for (; $e<=$#{$dta->{ENSEMBLE}}; $e++) { # post-recovery
+ $postrot += rot($e);
+}
+
+$prerot /= 360; # rotations, not degrees
+$dnrot /= 360;
+$uprot /= 360;
+$postrot /= 360;
+
+$prrms = sqrt($prrms/($lastgood-$firstgood));
+$dnprrms = sqrt($dnprrms/($atbottom-$firstgood));
+$upprrms = sqrt($upprrms/($lastgood-$atbottom));
+
+$rotrms = sqrt($rotrms/($lastgood-$firstgood));
+
+#======================================================================
+# PRODUCE OUTPUT
+#======================================================================
+
+printf(STDERR "# of ensembles : %d\n",scalar(@{$dta{ENSEMBLE}}));
+printf(STDERR "Start of cast : %s (#%5d) at %6.1fm\n",
+ $dta{ENSEMBLE}[$firstgood]->{TIME},
+ $dta{ENSEMBLE}[$firstgood]->{NUMBER},
+ $dta{ENSEMBLE}[$firstgood]->{DEPTH});
+printf(STDERR "Bottom of cast : %s (#%5d) at %6.1fm\n",
+ $dta{ENSEMBLE}[$atbottom]->{TIME},
+ $dta{ENSEMBLE}[$atbottom]->{NUMBER},
+ $dta{ENSEMBLE}[$atbottom]->{DEPTH});
+if (defined($water_depth)) {
+ printf(STDERR "Seabed : at %6.1fm (+-%dm)\n",$water_depth,$sig_wd);
+} else {
+ print(STDERR "Seabed : not found\n");
+}
+printf(STDERR "End of cast : %s (#%5d) at %6.1fm\n",
+ $dta{ENSEMBLE}[$lastgood]->{TIME},
+ $dta{ENSEMBLE}[$lastgood]->{NUMBER},
+ $dta{ENSEMBLE}[$lastgood]->{DEPTH});
+
+printf(STDERR "Rel. Displacement: x = %d(%d)m / y = %d(%d)m\n",
+ $dta{ENSEMBLE}[$lastgood]->{X}, $x_err,
+ $dta{ENSEMBLE}[$lastgood]->{Y}, $y_err,
+ ) if defined($opt_M);
+
+printf(STDERR "Cast Duration : %.1f hours (pinging for %.1f hours)\n",
+ $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME} / 3600,
+ ($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{UNIX_TIME} -
+ $dta{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);
+
+printf(STDERR "Minimum range : %dm at ensemble %d, beam %d\n",
+ $dta{DISTANCE_TO_BIN1_CENTER} +
+ $min_good_bins*$dta{BIN_LENGTH},
+ $dta{ENSEMBLE}[$min_good_ens]->{NUMBER},
+ $worst_beam);
+printf(STDERR "80%%-valid bins : %.1f\n",$gb+1);
+printf(STDERR "80%%-valid range : %dm\n",
+ $dta{DISTANCE_TO_BIN1_CENTER} + $gb*$dta{BIN_LENGTH});
+printf(STDERR "3-beam solutions : $RDI_Coords::threeBeam_1 " .
+ "$RDI_Coords::threeBeam_2 " .
+ "$RDI_Coords::threeBeam_3 " .
+ "$RDI_Coords::threeBeam_4\n")
+ unless ($opt_4);
+printf(STDERR "net rotations : [%d]/%d/%d/[%d]\n",$prerot,$dnrot,$uprot,$postrot);
+printf(STDERR "rms pitch/roll : %.1f/%.1f\n",$dnprrms,$upprrms);
+
+if ($opt_A) { # ANTS format
+ print("#ANTS# [] $USAGE\n");
+ $uFields = "{u} {u_err} {v} {v_err} {x} {x_err} {y} {y_err}"
+ if defined($opt_M);
+ print("#ANTS#FIELDS# {ens} {time} {elapsed} {secno} {downcast} " .
+ "{w} {w_err} {depth} {depth_err} {depth_BT} " .
+ "{pitchroll} {rotation} " .
+ "$uFields $addFields\n");
+
+ printf("#ANTS#PARAMS# date{$dta{ENSEMBLE}[$firstgood]->{DATE}} " .
+ "start_time{$dta{ENSEMBLE}[$firstgood]->{TIME}} " .
+ "bottom_time{$dta{ENSEMBLE}[$atbottom]->{TIME}} " .
+ "end_time{$dta{ENSEMBLE}[$lastgood]->{TIME}} " .
+ "bottom_xmit_voltage{$dta{ENSEMBLE}[$atbottom]->{ADC_XMIT_VOLTAGE}} " .
+ "bottom_xmit_current{$dta{ENSEMBLE}[$atbottom]->{ADC_XMIT_CURRENT}} " .
+ "pinging_duration{%.1f} " .
+ "cast_duration{%.1f} " .
+ "0.8_valid_bins{%.1f} " .
+ "0.8_valid_range{%.1f} " .
+ "max_depth{%.1f} " .
+ "depth_error{%.1f} " .
+ "min_range{%d} " .
+ "n_ensembles{%d} " .
+ "w_gap_time{%d} " .
+ "stderr_w{%.4f} " .
+ "rms_pitchroll{%.1f} " .
+ "downcast_rms_pitchroll{%.1f} " .
+ "upcast_rms_pitchroll{%.1f} " .
+ "rms_rotation{%.2f} " .
+ "deployment_rotations{%d} " .
+ "downcast_rotations{%d} " .
+ "upcast_rotations{%d} " .
+ "recovery_rotations{%d} " .
+ "bin1_dist{%.1f} " .
+ "bin_length{%.1f} " .
+ "\n",
+ ($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{UNIX_TIME} -
+ $dta{ENSEMBLE}[0]->{UNIX_TIME}),
+ $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
+ $gb+1,
+ $dta{DISTANCE_TO_BIN1_CENTER} + $gb*$dta{BIN_LENGTH},
+ $dta{ENSEMBLE}[$atbottom]->{DEPTH},
+ $dta{ENSEMBLE}[$lastgood]->{DEPTH} -
+ $dta{ENSEMBLE}[$firstgood]->{DEPTH},
+ $dta{DISTANCE_TO_BIN1_CENTER} +
+ $min_good_bins*$dta{BIN_LENGTH},
+ scalar(@{$dta{ENSEMBLE}}),
+ $w_gap_time,$wErr,$prrms,$dnprrms,$upprrms,$rotrms,
+ $prerot,$dnrot,$uprot,$postrot,
+ $dta{DISTANCE_TO_BIN1_CENTER},
+ $dta{BIN_LENGTH},
+ );
+ printf("#ANTS#PARAMS# magnetic_declination{$opt_M} " .
+ "uv_gap_time{%d} " .
+ "mean_u{%.4f} " .
+ "stderr_u{%.4f} " .
+ "dx{%d} " .
+ "dx_err{%d} " .
+ "mean_v{%.4f} " .
+ "stderr_v{%.4f} " .
+ "dy{%d} " .
+ "dy_err{%d}\n",
+ $uv_gap_time,
+ $dta{ENSEMBLE}[$lastgood]->{X} /
+ $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
+ $uErr, $dta{ENSEMBLE}[$lastgood]->{X}, $x_err,
+ $dta{ENSEMBLE}[$lastgood]->{Y} /
+ $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
+ $vErr, $dta{ENSEMBLE}[$lastgood]->{Y}, $y_err,
+ ) if defined ($opt_M);
+ print("#ANTS#PARAMS# start_lat{$s_lat} start_lon{$s_lon} " .
+ "end_lat{$e_lat} end_lon{$e_lon} " .
+ "lat{$lat} lon{$lon}\n")
+ if defined($lat);
+ if ($dta{TIME_BETWEEN_PINGS} == 0) {
+ print("#ANTS#PARAMS# pinging_rate{staggered}\n");
+ } else {
+ printf("#ANTS#PARAMS# pinging_rate{%.2f}\n",
+ 1/$dta{TIME_BETWEEN_PINGS});
+ }
+ printf("#ANTS#PARAMS# drift_x{%d} drift_y{%d} " .
+ "drift_u{%.3f} drift_v{%.3f} " .
+ "\n",$ddx,$ddy,$du,$dv) if defined($ddx);
+ if (defined($water_depth)) {
+ printf("#ANTS#PARAMS# water_depth{%d} sig-water_depth{%d}\n",
+ $water_depth,$sig_wd);
+ } else {
+ print("#ANTS#PARAMS# water_depth{nan} sig-water_depth{nan}\n");
+ }
+}
+
+sub p($) { print(defined($_[0])?"$_[0] ":"nan "); }
+sub pb($) { print($_[0]?"1 ":"0 "); }
+
+unless ($opt_Q) { # write profile
+ for ($e=$firstgood; $e<=$lastgood; $e++) {
+ p($dta{ENSEMBLE}[$e]->{NUMBER});
+ p($dta{ENSEMBLE}[$e]->{UNIX_TIME});
+ p($dta{ENSEMBLE}[$e]->{ELAPSED_TIME});
+ p($dta{ENSEMBLE}[$e]->{SECNO});
+ pb($dta{ENSEMBLE}[$e]->{UNIX_TIME} < $dta{ENSEMBLE}[$atbottom]->{UNIX_TIME});
+ p($dta{ENSEMBLE}[$e]->{W});
+ p($dta{ENSEMBLE}[$e]->{W_ERR});
+ p($dta{ENSEMBLE}[$e]->{DEPTH});
+ p($dta{ENSEMBLE}[$e]->{DEPTH_ERR});
+ p($dta{ENSEMBLE}[$e]->{DEPTH_BT});
+ p($dta{ENSEMBLE}[$e]->{PITCHROLL});
+ p($dta{ENSEMBLE}[$e]->{ROTATION});
+ if (defined($opt_M)) {
+ p($dta{ENSEMBLE}[$e]->{U}); p($dta{ENSEMBLE}[$e]->{U_ERR});
+ p($dta{ENSEMBLE}[$e]->{V}); p($dta{ENSEMBLE}[$e]->{V_ERR});
+ p($dta{ENSEMBLE}[$e]->{X}); p($dta{ENSEMBLE}[$e]->{X_ERR});
+ p($dta{ENSEMBLE}[$e]->{Y}); p($dta{ENSEMBLE}[$e]->{Y_ERR});
+ }
+ if (defined(@f)) {
+ foreach $f (@f) {
+ my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
+ $fn = $f unless defined($fn);
+ p(eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));
+ }
+ }
+ print("\n");
+ }
+}
+
+exit(0);
new file mode 100755
--- /dev/null
+++ b/scanBins
@@ -0,0 +1,74 @@
+#!/usr/bin/perl
+#======================================================================
+# S C A N B I N S
+# doc: Mon Jan 27 17:55:34 2003
+# dlm: Wed Sep 19 10:02:51 2007
+# (c) 2003 A.M. Thurnherr
+# uE-Info: 11 24 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# Collect Per-Bin Stats
+# NB: currently broken
+
+# HISTORY:
+# Jan 27, 2003: - created
+# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
+
+$0 =~ m{(.*)/[^/]+};
+require "$1/WorkhorseBinRead.pl";
+require "getopts.pl";
+
+die("Usage: $0 " .
+ "" .
+ "<RDI file>\n")
+ unless (&Getopts("") && @ARGV == 1);
+
+print(STDERR "Reading $ARGV[0]...");
+readData($ARGV[0],\%dta); # read data
+print(STDERR "done\n");
+
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {
+ checkEnsemble(\%dta,$e); # sanity checks
+ $nens++;
+
+ next unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[0][0]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[0][1]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[0][2]) &&
+ defined($dta{ENSEMBLE}[$e]->{VELOCITY}[0][3]));
+ $ngoodens++;
+
+ for ($b=0; $b<$dta{N_BINS}; $b++) { # collect stats
+ my($ngood) = defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])
+ + defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1])
+ + defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2])
+ + defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]);
+ if ($ngood == 4) { $ngood4[$b]++; }
+ elsif ($ngood == 3) { $ngood3[$b]++; }
+ else { $nbad[$b]++; }
+
+ for ($i=0; $i<4; $i++) {
+ if (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i])) {
+ $ngood[$b][$i]++;
+ $sumcor[$b][$i] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][$i];
+ $sumamp[$b][$i] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][$i];
+ }
+ }
+ }
+}
+
+printf("$ngoodens good ensembles out of $nens\n");
+for ($b=0; $b<$dta{N_BINS}; $b++) { # gen output
+ printf("%2d: vels: %3d%% 4-bin, %3d%% 3-bin, %3d%% bad; ",
+ $b+1,
+ 100*$ngood4[$b]/$ngoodens,100*$ngood3[$b]/$ngoodens,
+ 100*$nbad[$b]/$ngoodens);
+ printf("mean corr: %3d/%3d/%3d/%3d; mean amp: %3d/%3d/%3d/%3d",
+ $sumcor[$b][0]/$ngood[$b][0], $sumcor[$b][1]/$ngood[$b][1],
+ $sumcor[$b][2]/$ngood[$b][2], $sumcor[$b][3]/$ngood[$b][3],
+ $sumamp[$b][0]/$ngood[$b][0], $sumamp[$b][1]/$ngood[$b][1],
+ $sumamp[$b][2]/$ngood[$b][2], $sumamp[$b][3]/$ngood[$b][3]);
+ print("\n");
+}
+
+exit(0);
+
new file mode 100755
--- /dev/null
+++ b/test_velBeamToBeamPair
@@ -0,0 +1,97 @@
+#!/usr/bin/perl
+#======================================================================
+# T E S T _ V E L B E A M T O B E A M P A I R
+# doc: Thu May 21 13:40:27 2009
+# dlm: Thu May 21 15:03:47 2009
+# (c) 2009 A.M. Thurnherr
+# uE-Info: 89 24 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+require "RDI_Coords.pl";
+
+$dta{CONVEX_BEAM_PATTERN} = 1;
+$dta{BEAM_ANGLE} = 20;
+
+#---------------------------------------
+# Case 1: uplooker, bin 3 pointing south
+#---------------------------------------
+
+$dta{ENSEMBLE}[0]->{XDUCER_FACING_UP} = 1;
+$dta{ENSEMBLE}[0]->{HEADING} = 180;
+
+$dta{ENSEMBLE}[0]->{PITCH} = 0;
+$dta{ENSEMBLE}[0]->{ROLL} = 0;
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,2,1,0,0));
+@ant = velBeamToBeamPair(\%dta,0,2,1,0,0);
+die("1: $RDI[0] $ant[0]\n")
+ unless (abs($RDI[0]-$ant[0]) < 1e-6);
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,0,0,1,3));
+@ant = velBeamToBeamPair(\%dta,0,0,0,1,3);
+die("2: $RDI[1] $ant[2]\n")
+ unless (abs($RDI[1]-$ant[2]) < 1e-6);
+
+$dta{ENSEMBLE}[0]->{PITCH} = 20;
+$dta{ENSEMBLE}[0]->{ROLL} = 0;
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,2,1,2,1));
+@ant = velBeamToBeamPair(\%dta,0,2,1,2,1);
+die("3: $RDI[0] $ant[0]\n")
+ unless (abs($RDI[0]-$ant[0]) < 1e-6);
+die("4: $RDI[1] $ant[2]\n")
+ unless (abs($RDI[1]-$ant[2]) < 1e-6);
+
+$dta{ENSEMBLE}[0]->{PITCH} = 0;
+$dta{ENSEMBLE}[0]->{ROLL} = 15;
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,-2,-0.5,-2,-0.5));
+@ant = velBeamToBeamPair(\%dta,0,-2,-0.5,-2,-0.5);
+die("5: $RDI[0] $ant[0]\n")
+ unless (abs($RDI[0]-$ant[0]) < 1e-6);
+die("6: $RDI[1] $ant[2]\n")
+ unless (abs($RDI[1]-$ant[2]) < 1e-6);
+
+#-----------------------------------------
+# Case 1: downlooker, bin 3 pointing south
+#-----------------------------------------
+
+$dta{ENSEMBLE}[0]->{XDUCER_FACING_UP} = 0;
+$dta{ENSEMBLE}[0]->{HEADING} = 180;
+
+$dta{ENSEMBLE}[0]->{PITCH} = 0;
+$dta{ENSEMBLE}[0]->{ROLL} = 0;
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,2,1,0,0));
+@ant = velBeamToBeamPair(\%dta,0,2,1,0,0);
+die("7: $RDI[0] $ant[0]\n")
+ unless (abs($RDI[0]+$ant[0]) < 1e-6);
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,0,0,1,3));
+@ant = velBeamToBeamPair(\%dta,0,0,0,1,3);
+die("8: $RDI[1] $ant[2]\n")
+ unless (abs($RDI[1]-$ant[2]) < 1e-6);
+
+$dta{ENSEMBLE}[0]->{PITCH} = 20;
+$dta{ENSEMBLE}[0]->{ROLL} = 0;
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,2,1,2,1));
+@ant = velBeamToBeamPair(\%dta,0,2,1,2,1);
+die("9: $RDI[0] $ant[0]\n")
+ unless (abs($RDI[0]+$ant[0]) < 1e-6);
+die("10: $RDI[1] $ant[2]\n")
+ unless (abs($RDI[1]-$ant[2]) < 1e-6);
+
+$dta{ENSEMBLE}[0]->{PITCH} = 0;
+$dta{ENSEMBLE}[0]->{ROLL} = 15;
+
+@RDI = velInstrumentToEarth(\%dta,0,velBeamToInstrument(\%dta,-2,-0.5,-2,-0.5));
+@ant = velBeamToBeamPair(\%dta,0,-2,-0.5,-2,-0.5);
+die("11: $RDI[0] $ant[0]\n")
+ unless (abs($RDI[0]+$ant[0]) < 1e-6);
+die("12: $RDI[1] $ant[2]\n")
+ unless (abs($RDI[1]-$ant[2]) < 1e-6);
+
+#----------------------------------------------------------------------
+
+print(STDERR "ok\n");