1 #====================================================================== |
1 #====================================================================== |
2 # R D I _ U T I L S . P L |
2 # R D I _ U T I L S . P L |
3 # doc: Wed Feb 12 10:21:32 2003 |
3 # doc: Wed Feb 12 10:21:32 2003 |
4 # dlm: Fri Apr 12 09:22:10 2013 |
4 # dlm: Thu Jun 20 15:26:47 2013 |
5 # (c) 2003 A.M. Thurnherr |
5 # (c) 2003 A.M. Thurnherr |
6 # uE-Info: 44 68 NIL 0 0 72 2 2 4 NIL ofnI |
6 # uE-Info: 125 0 NIL 0 0 72 2 2 4 NIL ofnI |
7 #====================================================================== |
7 #====================================================================== |
8 |
8 |
9 # miscellaneous RDI-specific utilities |
9 # miscellaneous RDI-specific utilities |
10 |
10 |
11 # History: |
11 # History: |
40 # - immediately disabled this code becasue it does appear to make matters worse |
40 # - immediately disabled this code becasue it does appear to make matters worse |
41 # Sep 21, 2011: - added calculation of RMS heave acceleration |
41 # Sep 21, 2011: - added calculation of RMS heave acceleration |
42 # Mar 27, 2013: - BUG: 3-beam solutions were not used in ref_lr_w |
42 # Mar 27, 2013: - BUG: 3-beam solutions were not used in ref_lr_w |
43 # - disabled apparently unused code |
43 # - disabled apparently unused code |
44 # Apr 12, 2013: - added $min_pctg as optional parameter to mk_prof |
44 # Apr 12, 2013: - added $min_pctg as optional parameter to mk_prof |
|
45 # May 14, 2013: - added incident-velocity, w12 & w34 to mkProfile |
|
46 # Jun 5, 2013: - BUG: incident-flow warning was printed repeatedly |
|
47 # Jun 20, 2013: - BUG: warning had used &antsInfo() |
45 |
48 |
46 use strict; |
49 use strict; |
47 |
50 |
48 #====================================================================== |
51 #====================================================================== |
49 # fake_BT_RANGE(dta ptr) |
52 # fake_BT_RANGE(dta ptr) |
216 #====================================================================== |
219 #====================================================================== |
217 # ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz) = |
220 # ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz) = |
218 # mk_prof($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap); |
221 # mk_prof($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap); |
219 #====================================================================== |
222 #====================================================================== |
220 |
223 |
221 sub ref_lr_w($$$$$$$) # calc ref-level vert vels |
224 # calculate reference-layer vertical and incident velocities |
|
225 |
|
226 { my($warned); # static scope |
|
227 |
|
228 sub ref_lr_w($$$$$$$) |
222 { |
229 { |
223 my($dta,$ens,$rl_b0,$rl_b1,$min_corr,$max_e,$min_pctg) = @_; |
230 my($dta,$ens,$rl_b0,$rl_b1,$min_corr,$max_e,$min_pctg) = @_; |
224 my($i,@n,@bn,@v,@vel,@bv,@w); |
231 my($i,@n,@bn,@v,@vi,@vel,@veli,@bv,@w); |
|
232 my($w,$e,$nvi,$vi12,$vi43,@vbp,@velbp,@nbp,$w12,$w34,@w12,@w34); |
225 |
233 |
226 for ($i=$rl_b0; $i<=$rl_b1; $i++) { |
234 for ($i=$rl_b0; $i<=$rl_b1; $i++) { |
227 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0]) |
235 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0]) |
228 if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $min_corr); |
236 if ($dta->{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $min_corr); |
229 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1]) |
237 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][1]) |
239 if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < $min_pctg); |
247 if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < $min_pctg); |
240 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2]) |
248 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2]) |
241 if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < $min_pctg); |
249 if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < $min_pctg); |
242 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3]) |
250 undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3]) |
243 if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < $min_pctg); |
251 if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < $min_pctg); |
244 @v = velInstrumentToEarth($dta,$ens, |
252 @vi = velBeamToInstrument($dta,@{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]}); |
245 velBeamToInstrument($dta, |
253 @v = velInstrumentToEarth($dta,$ens,@vi); |
246 @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]})); |
254 @vbp = velBeamToBPEarth($dta,$ens,@{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]}); |
247 } else { |
255 } else { |
248 next if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 || |
256 next if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 || |
249 $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 || |
257 $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 || |
250 $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 || |
258 $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 || |
251 $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < $min_pctg); |
259 $dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < $min_pctg); |
252 @v = @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]}; |
260 @v = @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]}; |
253 # NB: no need to apply heading bias, as long as we only use w! |
261 unless ($warned) { |
|
262 print(STDERR "WARNING: incident-flow & beam-pair velocities not yet implemented for earth coordinates"); |
|
263 $warned = 1; |
|
264 } |
254 } |
265 } |
255 ### next if (!defined($v[3]) || abs($v[3]) > $max_e); # disallow 3-beam solutions |
266 ### next if (!defined($v[3]) || abs($v[3]) > $max_e); # disallow 3-beam solutions |
256 next if (defined($v[3]) && abs($v[3]) > $max_e); # allow 3-beam solutions |
267 next if (defined($v[3]) && abs($v[3]) > $max_e); # allow 3-beam solutions |
257 |
268 |
258 if (defined($v[2])) { # valid w |
269 if (defined($v[2])) { # valid velocity |
259 $vel[2] += $v[2]; $n[2]++; |
270 $vel[2] += $v[2]; $n[2]++; # vertical velocity |
260 $vel[3] += $v[3], $n[3]++ if defined($v[3]); |
271 $vel[3] += $v[3], $n[3]++ if defined($v[3]); # error velocity |
261 push(@w,$v[2]); # for stderr test |
272 push(@w,$v[2]); # save for stderr calculation |
|
273 } |
|
274 |
|
275 if (defined($vbp[1])) { # beam-pair vertical velocities |
|
276 $velbp[0] += $vbp[1]; $nbp[0]++; |
|
277 push(@w12,$vbp[1]); |
|
278 } |
|
279 if (defined($vbp[3])) { |
|
280 $velbp[1] += $vbp[3]; $nbp[1]++; |
|
281 push(@w34,$vbp[1]); |
262 } |
282 } |
263 |
283 |
|
284 if (defined($vi[0])) { # incident velocity |
|
285 $veli[0] += $vi[0]; |
|
286 $veli[1] += $vi[1]; |
|
287 $nvi++; |
|
288 } |
|
289 |
264 # The following code calculates beam-averaged ref-lr velocities. |
290 # The following code calculates beam-averaged ref-lr velocities. |
265 # I do not recall what this was implemented for. Disabled Mar 27, 2013. |
291 # I do not recall what this was implemented for. Disabled Mar 27, 2013. |
266 # |
292 # |
267 # if ($dta->{BEAM_COORDINATES}) { |
293 # if ($dta->{BEAM_COORDINATES}) { |
268 # $bv[0] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0], $bn[0]++ |
294 # $bv[0] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0], $bn[0]++ |
272 # $bv[2] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2], $bn[2]++ |
298 # $bv[2] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2], $bn[2]++ |
273 # if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2]); |
299 # if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][2]); |
274 # $bv[3] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3], $bn[3]++ |
300 # $bv[3] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3], $bn[3]++ |
275 # if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3]); |
301 # if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3]); |
276 # } |
302 # } |
277 } |
303 } # loop over ref-lr bins |
278 |
304 |
279 my($w) = $n[2] ? $vel[2]/$n[2] : undef; # w uncertainty |
305 $w = ($n[2] > 0) ? $vel[2]/$n[2] : undef; # calc means |
280 my($sumsq) = 0; |
306 $e = ($n[3] > 0) ? $vel[3]/$n[3] : undef; |
|
307 if ($nvi > 0) { |
|
308 $vi12 = $veli[0] / $nvi; |
|
309 $vi43 = $veli[1] / $nvi; |
|
310 } else { |
|
311 $vi12 = $vi43 = undef; |
|
312 } |
|
313 $w12 = ($nbp[0] > 0) ? $velbp[0]/$nbp[0] : undef; |
|
314 $w34 = ($nbp[1] > 0) ? $velbp[1]/$nbp[1] : undef; |
|
315 |
|
316 if (@w12) { # w uncertainty |
|
317 my($sumsq) = 0; |
|
318 for ($i=0; $i<=$#w12; $i++) { |
|
319 $sumsq += ($w12-$w12[$i])**2; |
|
320 } |
|
321 $dta->{ENSEMBLE}[$ens]->{W12} = $w12; |
|
322 $dta->{ENSEMBLE}[$ens]->{W12_ERR} = sqrt($sumsq)/($nbp[0]-1) |
|
323 if ($nbp[0]>=2); |
|
324 } |
|
325 |
|
326 if (@w34) { # w uncertainty |
|
327 my($sumsq) = 0; |
|
328 for ($i=0; $i<=$#w34; $i++) { |
|
329 $sumsq += ($w34-$w34[$i])**2; |
|
330 } |
|
331 $dta->{ENSEMBLE}[$ens]->{W34} = $w34; |
|
332 $dta->{ENSEMBLE}[$ens]->{W34_ERR} = sqrt($sumsq)/($nbp[1]-1) |
|
333 if ($nbp[1]>=2); |
|
334 } |
|
335 |
|
336 my($sumsq) = 0; # w uncertainty |
281 for ($i=0; $i<=$#w; $i++) { |
337 for ($i=0; $i<=$#w; $i++) { |
282 $sumsq += ($w-$w[$i])**2; |
338 $sumsq += ($w-$w[$i])**2; |
283 } |
339 } |
284 my($stderr) = $n[2]>=2 ? sqrt($sumsq)/($n[2]-1) : undef; |
340 my($stderr) = $n[2]>=2 ? sqrt($sumsq)/($n[2]-1) : undef; |
|
341 |
285 # The following stderr test introduces a huge gap near the bottom of |
342 # The following stderr test introduces a huge gap near the bottom of |
286 # the profiles. Without it, they seem more reasonable. |
343 # the profiles. Without it, they seem more reasonable. |
287 # next if ($stderr > 0.05); |
344 # next if ($stderr > 0.05); |
288 |
345 |
289 if (defined($w)) { # valid w |
346 if (defined($w)) { # valid velocity |
290 $dta->{ENSEMBLE}[$ens]->{W} = $w; |
347 $dta->{ENSEMBLE}[$ens]->{W} = $w; |
291 $dta->{ENSEMBLE}[$ens]->{W_ERR} = $stderr; |
348 $dta->{ENSEMBLE}[$ens]->{W_ERR} = $stderr; |
292 } |
349 } |
|
350 $dta->{ENSEMBLE}[$ens]->{ERR_VEL} = $e if (defined($e)); |
|
351 |
|
352 $dta->{ENSEMBLE}[$ens]->{W12} = $w12 if (defined($w12)); |
|
353 $dta->{ENSEMBLE}[$ens]->{W34} = $w34 if (defined($w34)); |
|
354 |
|
355 if (defined($vi12)) { |
|
356 $dta->{ENSEMBLE}[$ens]->{INCIDENT_VEL_T12} = $vi12; |
|
357 $dta->{ENSEMBLE}[$ens]->{INCIDENT_VEL_T43} = $vi43; |
|
358 } |
|
359 |
293 # The following code calculates beam-averaged ref-lr velocities. |
360 # The following code calculates beam-averaged ref-lr velocities. |
294 # I do not recall what this was implemented for. Disabled Mar 27, 2013. |
361 # I do not recall what this was implemented for. Disabled Mar 27, 2013. |
295 # |
362 # |
296 # if ($dta->{BEAM_COORDINATES}) { |
363 # if ($dta->{BEAM_COORDINATES}) { |
297 # $dta->{ENSEMBLE}[$ens]->{V1} = $bn[0]>=2 ? $bv[0]/$bn[0] : undef; |
364 # $dta->{ENSEMBLE}[$ens]->{V1} = $bn[0]>=2 ? $bv[0]/$bn[0] : undef; |
298 # $dta->{ENSEMBLE}[$ens]->{V2} = $bn[1]>=2 ? $bv[1]/$bn[1] : undef; |
365 # $dta->{ENSEMBLE}[$ens]->{V2} = $bn[1]>=2 ? $bv[1]/$bn[1] : undef; |
299 # $dta->{ENSEMBLE}[$ens]->{V3} = $bn[2]>=2 ? $bv[2]/$bn[2] : undef; |
366 # $dta->{ENSEMBLE}[$ens]->{V3} = $bn[2]>=2 ? $bv[2]/$bn[2] : undef; |
300 # $dta->{ENSEMBLE}[$ens]->{V4} = $bn[3]>=2 ? $bv[3]/$bn[3] : undef; |
367 # $dta->{ENSEMBLE}[$ens]->{V4} = $bn[3]>=2 ? $bv[3]/$bn[3] : undef; |
301 # } |
368 # } |
302 } |
369 |
|
370 } |
|
371 |
|
372 } # static scope |
303 |
373 |
304 |
374 |
305 sub mk_prof(...) # make profile |
375 sub mk_prof(...) # make profile |
306 { |
376 { |
307 my($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap,$min_pctg) = @_; |
377 my($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap,$min_pctg) = @_; |