0
|
1 |
#======================================================================
|
|
2 |
# R D I _ U T I L S . P L
|
|
3 |
# doc: Wed Feb 12 10:21:32 2003
|
|
4 |
# dlm: Sun May 23 16:35:21 2010
|
|
5 |
# (c) 2003 A.M. Thurnherr
|
|
6 |
# uE-Info: 156 42 NIL 0 0 72 2 2 4 NIL ofnI
|
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# miscellaneous RDI-specific utilities
|
|
10 |
|
|
11 |
# History:
|
|
12 |
# Feb 12, 2003: - created
|
|
13 |
# Feb 14, 2003: - added check for short (bad) data files
|
|
14 |
# Feb 26, 2004: - added Earth-coordinate support
|
|
15 |
# - added ensure_BT_RANGE()
|
|
16 |
# Mar 17, 2004: - set bad BT ranges to undef in ensure_BT_RANGE
|
|
17 |
# - calc mean/stddev in ensure_BT_RANGE
|
|
18 |
# Mar 20, 2004: - BUG: find_seabed() could bomb when not enough
|
|
19 |
# bottom guesses passed the mode_width test
|
|
20 |
# Mar 30, 2004: - added &soundSpeed()
|
|
21 |
# Nov 8, 2005: - WATER_DEPTH => Z_BT
|
|
22 |
# Dec 1, 2005: - renamed to RDI_Utils.pl
|
|
23 |
# - folded in mk_prof from [mkProfile]
|
|
24 |
# Nov 30, 2007: - adapted ref_lr_w() to 3-beam solutions
|
|
25 |
# Feb 1, 2008: - added comment
|
|
26 |
# Feb 22, 2008: - added ssCorr()
|
|
27 |
# Apr 9, 2008: - BUG: duplicate line of code (without effect) in find_seabed()
|
|
28 |
# - BUG: seabed < max depth was possible
|
|
29 |
# Jan 2010: - fiddled with seabed detection params (no conclusion)
|
|
30 |
# May 23, 2010: - renamed Z to DEPTH
|
|
31 |
|
|
32 |
use strict;
|
|
33 |
|
|
34 |
#======================================================================
|
|
35 |
# fake_BT_RANGE(dta ptr)
|
|
36 |
#======================================================================
|
|
37 |
|
|
38 |
# During cruise NBP0204 it was found that one of Visbeck's RDIs consistently
|
|
39 |
# returns zero as the bottom-track range, even thought the BT velocities
|
|
40 |
# are good. This functions calculates the ranges if they are missing.
|
|
41 |
|
|
42 |
sub cBTR($$$)
|
|
43 |
{
|
|
44 |
my($d,$e,$b) = @_;
|
|
45 |
my($maxamp) = -9e99;
|
|
46 |
my($maxi);
|
|
47 |
|
|
48 |
for (my($i)=0; $i<$d->{N_BINS}; $i++) {
|
|
49 |
next unless ($d->{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b] > $maxamp);
|
|
50 |
$maxamp = $d->{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b];
|
|
51 |
$maxi = $i;
|
|
52 |
}
|
|
53 |
$d->{ENSEMBLE}[$e]->{BT_RANGE}[$b] =
|
|
54 |
$d->{DISTANCE_TO_BIN1_CENTER} + $maxi * $d->{BIN_LENGTH};
|
|
55 |
}
|
|
56 |
|
|
57 |
sub ensure_BT_RANGE($)
|
|
58 |
{
|
|
59 |
my($d) = @_;
|
|
60 |
for (my($e)=0; $e<=$#{$d->{ENSEMBLE}}; $e++) {
|
|
61 |
my($sum) = my($n) = 0;
|
|
62 |
if (defined($d->{ENSEMBLE}[$e]->{BT_VELOCITY}[0])) {
|
|
63 |
for (my($b)=0; $b<4; $b++) {
|
|
64 |
cBTR($d,$e,$b)
|
|
65 |
unless defined($d->{ENSEMBLE}[$e]->{BT_RANGE}[$b]);
|
|
66 |
$sum += $d->{ENSEMBLE}[$e]->{BT_RANGE}[$b]; $n++;
|
|
67 |
}
|
|
68 |
} else {
|
|
69 |
for (my($b)=0; $b<4; $b++) {
|
|
70 |
$d->{ENSEMBLE}[$e]->{BT_RANGE}[$b] = undef;
|
|
71 |
}
|
|
72 |
}
|
|
73 |
$d->{ENSEMBLE}[$e]->{BT_MEAN_RANGE} = $sum/$n if ($n == 4);
|
|
74 |
}
|
|
75 |
}
|
|
76 |
|
|
77 |
#======================================================================
|
|
78 |
# (seabed depth, stddev) = find_seabed(dta ptr, btm ensno, coord flag)
|
|
79 |
#======================================================================
|
|
80 |
|
|
81 |
# NOTE FOR YOYOS:
|
|
82 |
# - this routine only detects the BT around the depeest depth!
|
|
83 |
# - this is on purpose, because it is used for [listBT]
|
|
84 |
|
|
85 |
# This is a PAIN:
|
|
86 |
# if the instrument is too close to the bottom, the BT_RANGE
|
|
87 |
# readings are all out; the only solution is to have a sufficiently
|
|
88 |
# large search width (around the max(depth) ensemble) so that
|
|
89 |
# a part of the up (oops, I forgot to finish this comment one year
|
|
90 |
# ago during A0304 and now I don't understand it any more :-)
|
|
91 |
|
|
92 |
my($search_width) = 200; # # of ensembles around bottom to search
|
|
93 |
my($mode_width) = 10; # max range of bottom around mode
|
|
94 |
my($min_dist) = 20; # min dist to seabed for good data
|
|
95 |
my($z_offset) = 10000; # shift z to ensure +ve array indices
|
|
96 |
|
|
97 |
sub find_seabed($$$)
|
|
98 |
{
|
|
99 |
my($d,$be,$beamCoords) = @_;
|
|
100 |
my($i,$dd,$sd,$nd);
|
|
101 |
my(@guesses);
|
|
102 |
|
|
103 |
return undef unless ($be-$search_width >= 0 &&
|
|
104 |
$be+$search_width <= $#{$d->{ENSEMBLE}});
|
|
105 |
for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
|
|
106 |
next unless (defined($d->{ENSEMBLE}[$i]->{DEPTH}) &&
|
|
107 |
defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[0]) &&
|
|
108 |
defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[1]) &&
|
|
109 |
defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[2]) &&
|
|
110 |
defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[3]));
|
|
111 |
my(@BT) = $beamCoords
|
|
112 |
? velInstrumentToEarth($d,$i,
|
|
113 |
velBeamToInstrument($d,
|
|
114 |
@{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}}))
|
|
115 |
: velApplyHdgBias($d,$i,@{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}});
|
|
116 |
next unless (abs($BT[3]) < 0.05);
|
|
117 |
$d->{ENSEMBLE}[$i]->{DEPTH_BT} =
|
|
118 |
$d->{ENSEMBLE}[$i]->{BT_RANGE}[0]/4 +
|
|
119 |
$d->{ENSEMBLE}[$i]->{BT_RANGE}[1]/4 +
|
|
120 |
$d->{ENSEMBLE}[$i]->{BT_RANGE}[2]/4 +
|
|
121 |
$d->{ENSEMBLE}[$i]->{BT_RANGE}[3]/4;
|
|
122 |
next unless ($d->{ENSEMBLE}[$i]->{DEPTH_BT} >= $min_dist);
|
|
123 |
$d->{ENSEMBLE}[$i]->{DEPTH_BT} *= -1
|
|
124 |
if ($d->{ENSEMBLE}[$i]->{XDUCER_FACING_UP});
|
|
125 |
$d->{ENSEMBLE}[$i]->{DEPTH_BT} += $d->{ENSEMBLE}[$i]->{DEPTH};
|
|
126 |
if ($d->{ENSEMBLE}[$i]->{DEPTH_BT} > $d->{ENSEMBLE}[$be]->{DEPTH}) {
|
|
127 |
$guesses[int($d->{ENSEMBLE}[$i]->{DEPTH_BT})+$z_offset]++;
|
|
128 |
$nd++;
|
|
129 |
} else {
|
|
130 |
undef($d->{ENSEMBLE}[$i]->{DEPTH_BT});
|
|
131 |
}
|
|
132 |
}
|
|
133 |
return undef unless ($nd>5);
|
|
134 |
|
|
135 |
my($mode,$nmax);
|
|
136 |
for ($i=0; $i<=$#guesses; $i++) { # find mode
|
|
137 |
$nmax=$guesses[$i],$mode=$i-$z_offset
|
|
138 |
if ($guesses[$i] > $nmax);
|
|
139 |
}
|
|
140 |
|
|
141 |
$nd = 0;
|
|
142 |
for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
|
|
143 |
next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
|
|
144 |
if (abs($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$mode) <= $mode_width) {
|
|
145 |
$dd += $d->{ENSEMBLE}[$i]->{DEPTH_BT};
|
|
146 |
$nd++;
|
|
147 |
} else {
|
|
148 |
$d->{ENSEMBLE}[$i]->{DEPTH_BT} = undef;
|
|
149 |
}
|
|
150 |
}
|
|
151 |
return undef unless ($nd >= 2);
|
|
152 |
|
|
153 |
$dd /= $nd;
|
|
154 |
for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
|
|
155 |
next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
|
|
156 |
$sd += ($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$dd)**2;
|
|
157 |
}
|
|
158 |
|
|
159 |
return ($dd, sqrt($sd/($nd-1)));
|
|
160 |
}
|
|
161 |
|
|
162 |
#======================================================================
|
|
163 |
# c = soundSpeed($salin,$temp,$depth)
|
|
164 |
#======================================================================
|
|
165 |
|
|
166 |
# Taken from the RDI BroadBand primer manual. The reference given there
|
|
167 |
# is Urick (1983).
|
|
168 |
|
|
169 |
sub soundSpeed($$$)
|
|
170 |
{
|
|
171 |
my($salin,$temp,$depth) = @_;
|
|
172 |
die("ERROR: soundSpeed($salin,$temp,$depth): non-numeric parameter\n")
|
|
173 |
unless numberp($salin) && numberp($temp) && numberp($depth);
|
|
174 |
return 1449.2 + 4.6*$temp -0.055*$temp**2 + 0.00029*$temp**3 +
|
|
175 |
(1.34 - 0.01*$temp) * ($salin - 35) + 0.016*$depth;
|
|
176 |
}
|
|
177 |
|
|
178 |
#======================================================================
|
|
179 |
# fac = ssCorr($eRef,$salin,$temp,$depth)
|
|
180 |
# $eRef : reference to current ensemble
|
|
181 |
# $salin: * -> use instrument salinity
|
|
182 |
# $temp : * -> use instrument temperature
|
|
183 |
# $depth: * -> use instrument PRESSURE(!)
|
|
184 |
#======================================================================
|
|
185 |
|
|
186 |
{ my($warned);
|
|
187 |
sub ssCorr($$$$)
|
|
188 |
{
|
|
189 |
my($eRef,$S,$T,$D) = @_;
|
|
190 |
$S = $eRef->{SALINITY} if ($S eq '*');
|
|
191 |
$T = $eRef->{TEMPERATURE} if ($T eq '*');
|
|
192 |
if ($D eq '*') {
|
|
193 |
print(STDERR "WARNING: soundspeed correction using instrument pressure instead of depth!\n")
|
|
194 |
unless ($warned);
|
|
195 |
$warned = 1;
|
|
196 |
$D = $eRef->{PRESSURE};
|
|
197 |
}
|
|
198 |
return soundSpeed($S,$T,$D) / $eRef->{SPEED_OF_SOUND};
|
|
199 |
}
|
|
200 |
}
|
|
201 |
|
|
202 |
#======================================================================
|
|
203 |
# ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz) =
|
|
204 |
# mk_prof($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap);
|
|
205 |
#======================================================================
|
|
206 |
|
|
207 |
sub ref_lr_w($$$$$$) # calc ref-level vert vels
|
|
208 |
{
|
|
209 |
my($dta,$ens,$rl_b0,$rl_b1,$min_corr,$max_e) = @_;
|
|
210 |
my($i,@n,@bn,@v,@vel,@bv,@w);
|
|
211 |
|
|
212 |
for ($i=$rl_b0; $i<=$rl_b1; $i++) {
|
|
213 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0])
|
|
214 |
if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $min_corr);
|
|
215 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1])
|
|
216 |
if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $min_corr);
|
|
217 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2])
|
|
218 |
if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $min_corr);
|
|
219 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3])
|
|
220 |
if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $min_corr);
|
|
221 |
if ($dta->{BEAM_COORDINATES}) {
|
|
222 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0])
|
|
223 |
if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100);
|
|
224 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1])
|
|
225 |
if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100);
|
|
226 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2])
|
|
227 |
if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100);
|
|
228 |
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3])
|
|
229 |
if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
|
|
230 |
@v = velInstrumentToEarth($dta,$ens,
|
|
231 |
velBeamToInstrument($dta,
|
|
232 |
@{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
|
|
233 |
} else {
|
|
234 |
next if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||
|
|
235 |
$dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||
|
|
236 |
$dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||
|
|
237 |
$dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
|
|
238 |
@v = @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]};
|
|
239 |
# NB: no need to apply heading bias, as long as we only use w!
|
|
240 |
}
|
|
241 |
next if (!defined($v[3]) || abs($v[3]) > $max_e);
|
|
242 |
|
|
243 |
if (defined($v[2])) { # valid w
|
|
244 |
$vel[2] += $v[2]; $n[2]++;
|
|
245 |
$vel[3] += $v[3], $n[3]++ if defined($v[3]);
|
|
246 |
push(@w,$v[2]); # for stderr test
|
|
247 |
}
|
|
248 |
|
|
249 |
if ($dta->{BEAM_COORDINATES}) {
|
|
250 |
$bv[0] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0], $bn[0]++
|
|
251 |
if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0]);
|
|
252 |
$bv[1] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1], $bn[1]++
|
|
253 |
if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1]);
|
|
254 |
$bv[2] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2], $bn[2]++
|
|
255 |
if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2]);
|
|
256 |
$bv[3] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3], $bn[3]++
|
|
257 |
if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3]);
|
|
258 |
}
|
|
259 |
}
|
|
260 |
|
|
261 |
my($w) = $n[2] ? $vel[2]/$n[2] : undef; # w uncertainty
|
|
262 |
my($sumsq) = 0;
|
|
263 |
for ($i=0; $i<=$#w; $i++) {
|
|
264 |
$sumsq += ($w-$w[$i])**2;
|
|
265 |
}
|
|
266 |
my($stderr) = $n[2]>=2 ? sqrt($sumsq)/($n[2]-1) : undef;
|
|
267 |
# The following stderr test introduces a huge gap near the bottom of
|
|
268 |
# the profiles. Without it, they seem more reasonable.
|
|
269 |
# next if ($stderr > 0.05);
|
|
270 |
|
|
271 |
if (defined($w)) { # valid w
|
|
272 |
$dta->{ENSEMBLE}[$ens]->{W} = $w;
|
|
273 |
$dta->{ENSEMBLE}[$ens]->{W_ERR} = $stderr;
|
|
274 |
}
|
|
275 |
if ($dta->{BEAM_COORDINATES}) {
|
|
276 |
$dta->{ENSEMBLE}[$ens]->{V1} = $bn[0]>=2 ? $bv[0]/$bn[0] : undef;
|
|
277 |
$dta->{ENSEMBLE}[$ens]->{V2} = $bn[1]>=2 ? $bv[1]/$bn[1] : undef;
|
|
278 |
$dta->{ENSEMBLE}[$ens]->{V3} = $bn[2]>=2 ? $bv[2]/$bn[2] : undef;
|
|
279 |
$dta->{ENSEMBLE}[$ens]->{V4} = $bn[3]>=2 ? $bv[3]/$bn[3] : undef;
|
|
280 |
}
|
|
281 |
}
|
|
282 |
|
|
283 |
|
|
284 |
sub mk_prof($$$$$$$$) # make profile
|
|
285 |
{
|
|
286 |
my($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap) = @_;
|
|
287 |
my($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz);
|
|
288 |
|
|
289 |
for (my($z)=0,my($e)=0; $e<=$#{$dta->{ENSEMBLE}}; $e++) {
|
|
290 |
checkEnsemble($dta,$e) if ($check);
|
|
291 |
|
|
292 |
filterEnsemble($dta,$e)
|
|
293 |
if (defined($filter) &&
|
|
294 |
$dta->{ENSEMBLE}[$e]->{PERCENT_GOOD}[0][0] > 0);
|
|
295 |
ref_lr_w($dta,$e,$lr_b0,$lr_b1,$min_corr,$max_e); # ref. layer w
|
|
296 |
|
|
297 |
if (defined($firstgood)) {
|
|
298 |
$dta->{ENSEMBLE}[$e]->{ELAPSED_TIME} = # time since start
|
|
299 |
$dta->{ENSEMBLE}[$e]->{UNIX_TIME} -
|
|
300 |
$dta->{ENSEMBLE}[$firstgood]->{UNIX_TIME};
|
|
301 |
} else {
|
|
302 |
if (defined($dta->{ENSEMBLE}[$e]->{W})) { # start of prof.
|
|
303 |
$firstgood = $lastgood = $e;
|
|
304 |
$dta->{ENSEMBLE}[$e]->{ELAPSED_TIME} = 0;
|
|
305 |
$dta->{ENSEMBLE}[$e]->{DEPTH} = $dta->{ENSEMBLE}[$e]->{DEPTH_ERR} = 0;
|
|
306 |
}
|
|
307 |
next;
|
|
308 |
}
|
|
309 |
|
|
310 |
#--------------------------------------------------
|
|
311 |
# within profile: both $firstgood and $lastgood set
|
|
312 |
#--------------------------------------------------
|
|
313 |
|
|
314 |
if (!defined($dta->{ENSEMBLE}[$e]->{W})) { # gap
|
|
315 |
$w_gap_time += $dta->{ENSEMBLE}[$e]->{UNIX_TIME} -
|
|
316 |
$dta->{ENSEMBLE}[$e-1]->{UNIX_TIME};
|
|
317 |
next;
|
|
318 |
}
|
|
319 |
|
|
320 |
my($dt) = $dta->{ENSEMBLE}[$e]->{UNIX_TIME} - # time step since
|
|
321 |
$dta->{ENSEMBLE}[$lastgood]->{UNIX_TIME}; # ... last good ens
|
|
322 |
|
|
323 |
if ($dt > $max_gap) {
|
|
324 |
printf(STDERR "WARNING: %d-s gap too long, profile restarted at ensemble $e\n",$dt);
|
|
325 |
$firstgood = $lastgood = $e;
|
|
326 |
$dta->{ENSEMBLE}[$e]->{ELAPSED_TIME} = 0;
|
|
327 |
$z = $zErr = $maxz = 0;
|
|
328 |
$dta->{ENSEMBLE}[$e]->{DEPTH} = $dta->{ENSEMBLE}[$e]->{DEPTH_ERR} = 0;
|
|
329 |
$w_gap_time = 0;
|
|
330 |
next;
|
|
331 |
}
|
|
332 |
|
|
333 |
#-----------------------------------
|
|
334 |
# The current ensemble has a valid w
|
|
335 |
#-----------------------------------
|
|
336 |
|
|
337 |
$z += $dta->{ENSEMBLE}[$lastgood]->{W} * $dt; # integrate
|
|
338 |
$zErr += ($dta->{ENSEMBLE}[$lastgood]->{W_ERR} * $dt)**2;
|
|
339 |
$dta->{ENSEMBLE}[$e]->{DEPTH} = $z;
|
|
340 |
$dta->{ENSEMBLE}[$e]->{DEPTH_ERR} = sqrt($zErr);
|
|
341 |
|
|
342 |
$atbottom = $e, $maxz = $z if ($z > $maxz);
|
|
343 |
|
|
344 |
$lastgood = $e;
|
|
345 |
}
|
|
346 |
|
|
347 |
filterEnsembleStats() if defined($filter);
|
|
348 |
|
|
349 |
return ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz);
|
|
350 |
}
|
|
351 |
|
|
352 |
#----------------------------------------------------------------------
|
|
353 |
# (true|false) = numberp(var)
|
|
354 |
#----------------------------------------------------------------------
|
|
355 |
|
|
356 |
sub numberp(@)
|
|
357 |
{ return $_[0] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/; }
|
|
358 |
|
|
359 |
|
|
360 |
1;
|