# HG changeset patch # User A.M. Thurnherr # Date 1275638608 14400 # Node ID 229a0d72d2abfa76ec9a15857d3edbc1a316a5fe first hg version diff --git a/INDEX b/INDEX 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 + diff --git a/RDI2grd b/RDI2grd 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 ] [-r)ange ] " . + "[output -b)ase ] [-d)imensional coordinates] " . + "\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); diff --git a/RDI_BB_Read.pl b/RDI_BB_Read.pl 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 diff --git a/RDI_Coords.pl b/RDI_Coords.pl 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; diff --git a/RDI_Utils.pl b/RDI_Utils.pl 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; diff --git a/WorkhorseUtils.pl b/WorkhorseUtils.pl 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; diff --git a/beamStats b/beamStats 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 +# - added -b)eam coordinate output +# Feb 12, 2008: - modified 3-beam output +# - added -p)ct_good +# 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 ] " . + "[output -f)ile ] " . + "[-a)ll ens (not just those with good vels)] " . + "[-M)agnetic ] " . + "[-S)oundspeed correction " . + "[require -4)-beam solutions] [-d)iscard ] " . + "[-%)good ] " . + "[output -b)eam coordinates] " . + "\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); diff --git a/libRDI_Coords.pl b/libRDI_Coords.pl new file mode 120000 --- /dev/null +++ b/libRDI_Coords.pl @@ -0,0 +1,1 @@ +RDI_Coords.pl \ No newline at end of file diff --git a/listBT b/listBT 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 ] [-F)ilter ensembles