|
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; |