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;