#======================================================================
# L A D C P P R O C . S H E A R M E T H O D
# doc: Fri Sep 17 20:27:53 2010
# dlm: Sun Jul 27 19:44:41 2014
# (c) 2010 A.M. Thurnherr & E. Firing
# uE-Info: 338 65 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# PERLified functions from E. Firing's [merge.c], modified by A.M. Thurnherr
# NOTES:
# - velocity integration removed
# - percent-good flag removed (no point for single-ping ensembles)
# - need the following per-ensemble fields from previous steps:
# DEPTH
# CTD_SVEL
# - negative depths are not allowed (should not happen given DEPTH)
# WEIRDNESSES IN ERIC'S CODE:
# - w reference layer in set_misc_flags is calculated without taking
# the ref layer bins into account
# - u,v ref lr vels in set_wake_flag use w ref layer bins
# - distance to bin 1 center is not sound-speed corrected [WEIRDNESS CORRECTED]
# - $tilt calculation is wrong. I do not understand this comment in 2014.
# HISTORY:
# Sep 17, 2010: - created
# Oct 13, 2010: - first working version
# Oct 14, 2010: - renamed from LADCPshear.UHcode
# Oct 15, 2010: - added support for -pPI edit suppresion
# Oct 19, 2010: - reversed semantics of -p
# Oct 25, 2010: - added W_CTD_BIT, renamed W_BIT to W_OUTLIER_BIT
# Oct 26, 2010: - added TILT_BIT
# Dec 10, 2010: - modified assertion to allow processing of UL data
# Jul 10, 2011: - added outTDseries() call
# Jul 12, 2011: - replaced -p by $PPI_editing_enabled flag
# Feb 19, 2012: - added elapsed time to binned shear output
# Apr 11, 2012: - added MISSING_CTD_DATA_BIT
# Sep 25, 2013: - added code to calc gridded lat/lon info
# Nov 12, 2013: - BUG: correlation editing removed most (all?) 3-beam
# solutions
# - BUG: set_shear_flag() calculated shdev (slightly?)
# wrongly
# Mar 4, 2014: - added support for missing PITCH/ROLL (TILT) & HEADING
# Mar 21, 2014: - moved depthOfBin() to [LADCPproc.utils]
# Jul 15, 2014: - BUG: occasional missing w values were used in w_ref calc;
# (potentially useless) tests for missing w were added wherever
# edit_flags is tested or w is used; this change greatly
# increase the number of shear samples available in case of
# DoMORE1 tow-yo#1 data
# - disabled two assertions about depths of shear bins, one of
# which was violated by DoMORE1 uplooker data (apparently
# valid data above the sea surface in the uplooker)
# Jul 27, 2014: - renamed from [LADCPproc.UHmethod] because the code has
# diverged quite a bit from the original mplementation. However,
# the shear editing core remains as in the UH code
# - added calculation of PACKAGE_VELOCITY (-v)
#======================================================================
# VELOCITY EDITING
#======================================================================
my($BADVEL_BIT) = 0x01;
my($ERRVEL_BIT) = 0x02;
my($CORREL_BIT) = 0x04;
my($W_OUTLIER_BIT) = 0x08;
my($SHEAR_BIT) = 0x10;
my($SIDELOBE_BIT) = 0x20;
my($WAKE_BIT) = 0x40;
my($PPI_BIT) = 0x80;
my($W_CTD_BIT) = 0x100;
my($TILT_BIT) = 0x200;
my($DELTA_TILT_BIT) = 0x400;
my($MISSING_CTD_DATA_BIT) = 0x800;
my(%flag_count);
sub set_wake_flag($$)
{
my($ens,$De) = @_;
my($n) = 0;
my($uref) = my($vref) = my($wref) = 0;
for (my($bin)=$wbin_start-1; $bin<$wbin_end; $bin++) { # calc crude ref lr vel from all data
next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
$uref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U];
$vref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V];
$wref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W];
$n++;
}
return if ($n==0);
$uref /= $n;
$vref /= $n;
$wref /= $n;
## if upward (=negative) velocity greater than minimum, calculate wake
## heading and inclination
if (defined($LADCP{ENSEMBLE}[$ens]->{HEADING}) && $wref<-$min_wake_w) {
my($wake_hd) = 180 / 3.14159265358979 * atan2($uref,$vref);
my($speed) = sqrt($uref*$uref + $vref*$vref);
my($wake_ang)= abs(180 / 3.14159265358979 * atan($speed/$wref));
$wake_hd += 360
if ($wake_hd < 0);
$LADCP{ENSEMBLE}[$ens]->{HEADING} += 360
if ($LADCP{ENSEMBLE}[$ens]->{HEADING} < 0);
my($wake_mod) = $wake_hd % 90; # % returns integer part of remainder, but that's sufficient
my($adcp_mod) = $LADCP{ENSEMBLE}[$ens]->{HEADING} % 90;
if (((abs($wake_mod - $adcp_mod) < $wake_hd_dif) ||
(abs($wake_mod - $adcp_mod) > (90 - $wake_hd_dif))) &&
($wake_ang > $wake_ang_min)) {
for (my($bin)=0; $bin<$n_wake_bins; $bin++) {
$flag_count{$WAKE_BIT}[$bin]++;
$edit_flags[$ens][$bin] |= $WAKE_BIT;
}
}
} ## if ($wref < -min_wake_w)
## This does not make a lot of sense, because it trims points
## on only one side of the wake error, and that side depends
## on whether the integration is forward or backward.
if (($edit_flags[$ens+$De][0]&$WAKE_BIT) &&
($edit_flags[$ens][0]&$WAKE_BIT == 0)) {
for (my($bin)=0; $bin<$n_wake_bins; $bin++) {
$flag_count{$WAKE_BIT}[$bin]++;
$edit_flags[$ens][$bin] |= $WAKE_BIT;
}
}
}
sub set_misc_flags($$)
{
my($ens,$De) = @_;
my($ww) = my($n) = 0;
my($SLIfac) = 1 - cos(rad($LADCP{BEAM_ANGLE}));
for (my($bin)=0; $bin<$w_ref_bin; $bin++) { # ref-lr w
next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
$ww += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W];
$n++;
}
$ww /= $n if ($n > 0);
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
## We use the standard criterion for bottom interference; e.g. for
## 30-degree beams, the last 15% of the profile is contaminated
## by the sidelobe bounce. 1.5 bin length is added to allow for
## the length of the bin and pulse, that is, contamination of part of a
## bin. Profiler tilt does not require a more stringent criterion.
if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN}) {
if (numberp($water_depth) &&
$water_depth - &depthOfBin($ens,$bin) <=
$SLIfac * ($water_depth - $LADCP{ENSEMBLE}[$ens]->{DEPTH})
+ 1.5 * $LADCP{BIN_LENGTH}) {
$edit_flags[$ens][$bin] |= $SIDELOBE_BIT;
$flag_count{$SIDELOBE_BIT}++;
}
} else { ## upward-looking
if (&depthOfBin($ens,$bin) <=
$SLIfac * $LADCP{ENSEMBLE}[$ens]->{DEPTH}
+ 1.5 * $LADCP{BIN_LENGTH}) {
$edit_flags[$ens][$bin] |= $SIDELOBE_BIT;
$flag_count{$SIDELOBE_BIT}++;
}
}
if ($ww != 0 &&
abs($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W] - $ww) > $w_dif) {
$flag_count{$W_OUTLIER_BIT}++;
$edit_flags[$ens][$bin] |= $W_OUTLIER_BIT;
}
if (abs($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$E]) > $e_max) {
$flag_count{$ERRVEL_BIT}++;
$edit_flags[$ens][$bin] |= $ERRVEL_BIT;
}
my($nBadCorr) = 0;
for (my($beam)=0; $beam<=3; $beam++) {
$nBadCorr++
if ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] < $min_cor);
}
if (abs($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$E]) > 0) { # 4-beam solution
if ($nBadCorr > 0) {
$flag_count{$CORREL_BIT}++;
$edit_flags[$ens][$bin] |= $CORREL_BIT;
}
} else { # 3-beam solution
if ($nBadCorr > 1) {
$flag_count{$CORREL_BIT}++;
$edit_flags[$ens][$bin] |= $CORREL_BIT;
}
}
if ($bin < $shbin_start-1 || $bin >= $shbin_end) { # manually remove vels outside shear bin range
$edit_flags[$ens][$bin] |= $BADVEL_BIT;
$flag_count{$BADVEL_BIT}++;
}
} # for ($bin=0...
}
## The following is for editing out the second bottom bounce.
# - in the UH code, tilt = max(pitch,roll)
# - using the real tilt (here) implies that PPI editing is too conservative
# in case of large tilts
# - since, however, the sound speed at the transducer is used instead
# of the mean soundspeed below the ADCP, the difference is unlikely
# to matter
sub set_PPI_flags($$)
{
my($ens,$De) = @_;
my($clip_z0,$clip_z1);
my($dt_ping) = $LADCP{ENSEMBLE}[$ens]->{UNIX_TIME} - $LADCP{ENSEMBLE}[$ens-1]->{UNIX_TIME};
if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN}) {
if (numberp($water_depth)) {
$clip_z1 = $water_depth
- $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
* cos(rad($LADCP{BEAM_ANGLE} + $LADCP{ENSEMBLE}[$ens]->{TILT}))
+ $clip_margin;
$clip_z0 = $water_depth
- $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
* cos(rad($LADCP{BEAM_ANGLE} - $LADCP{ENSEMBLE}[$ens]->{TILT}))
- $clip_margin;
}
} else { # upward-looking
$clip_z1 = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
* cos(rad($LADCP{BEAM_ANGLE} - $LADCP{ENSEMBLE}[$ens]->{TILT}))
+ $clip_margin;
$clip_z0 = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
* cos(rad($LADCP{BEAM_ANGLE} + $LADCP{ENSEMBLE}[$ens]->{TILT}))
- $clip_margin;
}
for (my($bin)=$first_clip_bin-1; $bin<$LADCP{N_BINS}; $bin++) {
next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
my($dob) = depthOfBin($ens,$bin);
if (!defined($LADCP{ENSEMBLE}[$ens]->{TILT}) || ($dob>=$clip_z0 && $dob<=$clip_z1)) {
$edit_flags[$ens][$bin] |= $PPI_BIT;
$flag_count{$PPI_BIT}++;
}
}
}
sub edit_velocity($$)
{
my($start,$end) = @_; # ensemble indices
my($De) = $start<$end ? 1 : -1; # downcast: De = 1; upcast: De = -1
$flag_count{$WAKE_BIT} = $flag_count{$W_OUTLIER_BIT} = $flag_count{$ERRVEL_BIT} =
$flag_count{$CORREL_BIT} = $flag_count{$SHEAR_BIT} = $flag_count{$BADVEL_BIT} =
$flag_count{$SIDELOBE_BIT} = $flag_count{$PPI_BIT} = $flag_count{$W_CTD_BIT} =
$flag_count{$TILT_BIT} = $flag_count{$DELTA_TILT_BIT} = $flag_count{$MISSING_CTD_DATA_BIT} = 0;
for (my($ens)=$start; $ens!=$end+$De; $ens+=$De) { # loop over all ens from start to end
next unless ($LADCP{ENSEMBLE}[$ens]->{W});
if (abs($LADCP{ENSEMBLE}[$ens]->{CTD_W}-$LADCP{ENSEMBLE}[$ens]->{W}) > $w_max_err) { # get rid of aliased vels (ambiguity)
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
$edit_flags[$ens][$bin] |= $W_CTD_BIT;
$flag_count{$W_CTD_BIT}++;
}
next;
}
if (!defined($LADCP{ENSEMBLE}[$ens]->{TILT}) ||
$LADCP{ENSEMBLE}[$ens]->{TILT}>$max_tilt) { # get rid ensembles with large tilt
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
$edit_flags[$ens][$bin] |= $TILT_BIT;
$flag_count{$TILT_BIT}++;
}
next;
}
unless (numberp($LADCP{ENSEMBLE}[$ens]->{DEPTH}) && # get rid of ensembles with insufficient CTD info
numberp($LADCP{ENSEMBLE}[$ens]->{CTD_SVEL})) {
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
$edit_flags[$ens][$bin] |= $MISSING_CTD_DATA_BIT;
$flag_count{$MISSING_CTD_DATA_BIT}++;
}
next;
} # get rid ensembles after large rotation
if (defined($LADCP{ENSEMBLE}[$ens]->{TILT}) &&
defined($LADCP{ENSEMBLE}[$ens-$De]->{TILT}) &&
abs($LADCP{ENSEMBLE}[$ens]->{TILT}-$LADCP{ENSEMBLE}[$ens-$De]->{TILT}) > $max_delta_tilt) {
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
$edit_flags[$ens][$bin] |= $DELTA_TILT_BIT;
$flag_count{$DELTA_TILT_BIT}++;
}
next;
}
for (my($bin)=$shbin_start-1; $bin<$shbin_end; $bin++) { # flag bad velocities
$edit_flags[$ens][$bin] |= $BADVEL_BIT,$flag_count{$BADVEL_BIT}++
unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]);
}
set_wake_flag($ens,$De);
set_misc_flags($ens,$De);
set_PPI_flags($ens,$De)
if $PPI_editing_enabled && ($clip_margin > 0); # PPI editing is off by default
}
}
#==============================================================================
# CALCULATE VELOCITY SHEAR
# - final output in @ush_mu,@vsh_mu,@wsh_mu,@ush_sig,@vsh_sig,@wsh_sig
# NEW (ant): elapsed time output in @esh_mu
# lat/lon output in @lash_mu, @losh_mu
# - @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh are defined "local" in &calc_shear
#==============================================================================
#----------------------------------------------------------------------
# uv_to_shear(ens,calc_pkg_vel)
# - sets @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh for a given ensemble
# - also sets $LADCP{ENSEMBLE}[$ens]->{PACKAGE_VELOCITY} on -v
#----------------------------------------------------------------------
sub uv_to_shear($$)
{
my($ens,$calc_pkg_vel) = @_;
my($nvel) = 0;
my(@pu,@pv);
@sh_i0 = @sh_i1 = ();
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { # select valid velocities
next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
$nvel++;
push(@sh_i1,$bin) if (@sh_i0);
push(@sh_i0,$bin) if ($bin < $LADCP{N_BINS}-1);
if (defined($calc_pkg_vel)) {
my($ovgi) = int(&depthOfBin($ens,$bin) / $GRID_DZ);
push(@pu,$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$ovu[$ovgi]);
push(@pv,$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V]-$ovv[$ovgi]);
}
}
if (defined($calc_pkg_vel)) { # calculate package velocity
$LADCP{ENSEMBLE}[$ens]->{PACKAGE_VELOCITY}[$U] = avg(@pu);
$LADCP{ENSEMBLE}[$ens]->{PACKAGE_VELOCITY}[$V] = avg(@pv);
}
@dsh = ();
for (my($i)=0; $i<@sh_i1; $i++) { # calc and bin shears
my($d0) = &depthOfBin($ens,$sh_i0[$i]);
my($d1) = &depthOfBin($ens,$sh_i1[$i]);
next unless ($d0>=0 && $d1>=0);
# die("$0: assertion failed (ens=$ens i=$i depth=$LADCP{ENSEMBLE}[$ens]->{DEPTH} sh_i0[$i]=$sh_i0[$i] sh_i1[$i]=$sh_i1[$i] d0=$d0 d1=$d1)")
# unless ($d0>=0 && $d1>=0);
# die("$0: assertion failed (ens=$ens i=$i sh_i0[$i]=$sh_i0[$i] sh_i1[$i]=$sh_i1[$i] d0=$d0 d1=$d1)")
# unless (($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN} && $d1-$d0>0) ||
# ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} && $d1-$d0<0));
$dsh[$i] = ($d1 + $d0) / 2;
$ush[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$U] -
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$U]) / ($d1-$d0);
$vsh[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$V] -
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$V]) / ($d1-$d0);
$wsh[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$W] -
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$W]) / ($d1-$d0);
$ens[$i] = $ens;
}
return $nvel;
}
# here ush_mu, sh_n, etc are still set from the pre-gridding pass
sub set_shear_flag($)
{
my($ens) = @_;
my(@ibad,@ush_dev,@vsh_dev);
for (my($i)=0; $i<@dsh; $i++) {
die("$0: assertion failed") unless numberp($dsh[$i]);
my($bsgi) = int($dsh[$i] / $SHEAR_PREGRID_DZ);
$ush_dev[$i] = $ush[$i] - $ush_mu[$bsgi];
$vsh_dev[$i] = $vsh[$i] - $vsh_mu[$bsgi];
push(@ibad,$i) if ($ush_sig[$i] > 0 &&
(abs($ush_dev[$i]/$ush_sig[$i]) > $max_shdev ||
abs($vsh_dev[$i]/$vsh_sig[$i]) > $max_shdev));
} ## end of loop through shears
## Look for internal glitches: a positive shear followed
## immediately by a compensating negative shear, for
## example. When one is found, flag the common velocity
## sample, and untag the two shears by setting each ibad
## to -1.
for (my($bi)=0; $bi<@ibad-1; $bi++) {
next unless ($ibad[$bi]+1 == $ibad[$bi+1]);
my($i) = $ibad[$bi];
my($bsgi) = int($dsh[$i] / $SHEAR_PREGRID_DZ);
if ($ush_sig[$bsgi] > 0 && $vsh_sig[$bsgi] > 0 &&
sqrt(($ush_dev[$i]+$ush_dev[$i+1])**2/$ush_sig[$bsgi]**2) < $max_shdev_sum &&
sqrt(($vsh_dev[$i]+$vsh_dev[$i+1])**2/$vsh_sig[$bsgi]**2) < $max_shdev_sum) {
$flag_count{$SHEAR_BIT}++;
$edit_flags[$ens][$sh_i1[$i]] |= $SHEAR_BIT;
$ibad[$bi] = $ibad[$bi+1] = -1;
}
} ## end of first loop through bad shears
## Now flag all remaining velocities involved in the shears
## listed by ibad.
for (my($bi)=0; $bi<@ibad; $bi++) {
next if ($ibad[$bi] < 0);
$flag_count{$SHEAR_BIT} += 2;
$edit_flags[$ens][$sh_i0[$ibad[$bi]]] |= $SHEAR_BIT;
$edit_flags[$ens][$sh_i1[$ibad[$bi]]] |= $SHEAR_BIT;
}
}
sub calc_shear($$$$)
{
my($start,$end,$grid_dz,$edit_shear) = @_;
my($De) = $start<$end ? 1 : -1; # downcast: De = 1; upcast: De = -1
local(@ush_vals,@vsh_vals,@wsh_vals,@ens_vals);
my($nvel,$nsh) = (0,0);
for (my($ens)=$start; $ens!=$end+$De; $ens+=$De) { # loop over all ens from start to end
local(@sh_i0,@sh_i1);
local(@dsh,@ush,@vsh,@wsh,@ens);
uv_to_shear($ens,0);
if ($edit_shear) {
set_shear_flag($ens);
$nvel += uv_to_shear($ens,defined($opt_v));
$nsh += @dsh;
}
for (my($i)=0; $i<@dsh; $i++) { # save shears for binning calculations
my($gi) = int($dsh[$i] / $grid_dz);
push(@{$ush_vals[$gi]},$ush[$i]);
push(@{$vsh_vals[$gi]},$vsh[$i]);
push(@{$wsh_vals[$gi]},$wsh[$i]);
push(@{$ens_vals[$gi]},$ens[$i]);
}
} # $ens loop
outTDseries($De==1) if ($edit_shear); # output depth-time time series
@ush_mu = @vsh_mu = @wsh_mu = @esh_mu = @lash_mu = @losh_mu = ();
@ush_sig = @vsh_sig = @wsh_sig = ();
for (my($gi)=0; $gi<@ush_vals; $gi++) { # calc grid means & stddev
my($sum_ush,$sum_vsh,$sum_wsh,$sum_esh,$sum_lash,$sum_losh);
$sh_n[$gi] = @{$ush_vals[$gi]};
for (my($vi)=0; $vi<$sh_n[$gi]; $vi++) {
$sum_ush += $ush_vals[$gi][$vi];
$sum_vsh += $vsh_vals[$gi][$vi];
$sum_wsh += $wsh_vals[$gi][$vi];
$sum_esh += $LADCP{ENSEMBLE}[$ens_vals[$gi][$vi]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l;
$sum_lash += $LADCP{ENSEMBLE}[$ens_vals[$gi][$vi]]->{CTD_LAT};
$sum_losh += $LADCP{ENSEMBLE}[$ens_vals[$gi][$vi]]->{CTD_LON};
}
$ush_mu[$gi] = $sh_n[$gi] ? $sum_ush/$sh_n[$gi] : nan;
$vsh_mu[$gi] = $sh_n[$gi] ? $sum_vsh/$sh_n[$gi] : nan;
$wsh_mu[$gi] = $sh_n[$gi] ? $sum_wsh/$sh_n[$gi] : nan;
$esh_mu[$gi] = $sh_n[$gi] ? $sum_esh/$sh_n[$gi] : nan;
$lash_mu[$gi] = $sh_n[$gi] ? $sum_lash/$sh_n[$gi] : nan;
$losh_mu[$gi] = $sh_n[$gi] ? $sum_losh/$sh_n[$gi] : nan;
}
for (my($gi)=0; $gi<@ush_vals; $gi++) { # calc & grid stddevs
my($sumsq_ush,$sumsq_vsh,$sumsq_wsh);
for (my($vi)=0; $vi<$sh_n[$gi]; $vi++) {
$sumsq_ush += ($ush_vals[$gi][$vi] - $ush_mu[$gi])**2;
$sumsq_vsh += ($vsh_vals[$gi][$vi] - $vsh_mu[$gi])**2;
$sumsq_wsh += ($wsh_vals[$gi][$vi] - $wsh_mu[$gi])**2;
}
$ush_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_ush/($sh_n[$gi]-1)) : nan;
$vsh_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_vsh/($sh_n[$gi]-1)) : nan;
$wsh_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_wsh/($sh_n[$gi]-1)) : nan;
}
if ($edit_shear && $opt_d) {
print(STDERR "\n\t\t$nvel valid velocities");
print(STDERR "\n\t\t$nsh valid shears");
print(STDERR "\n\t\tflag counts:");
print(STDERR "\n\t\t\tBADVEL_BIT = $flag_count{$BADVEL_BIT}")
if ($flag_count{$BADVEL_BIT});
print(STDERR "\n\t\t\tERRVEL_BIT = $flag_count{$ERRVEL_BIT}")
if ($flag_count{$ERRVEL_BIT});
print(STDERR "\n\t\t\tCORREL_BIT = $flag_count{$CORREL_BIT}")
if ($flag_count{$W_OUTLIER_BIT});
print(STDERR "\n\t\t\tW_OUTLIER_BIT = $flag_count{$W_OUTLIER_BIT}")
if ($flag_count{$W_OUTLIER_BIT});
print(STDERR "\n\t\t\tSHEAR_BIT = $flag_count{$SHEAR_BIT}")
if ($flag_count{$SIDELOBE_BIT});
print(STDERR "\n\t\t\tSIDELOBE_BIT = $flag_count{$SIDELOBE_BIT}")
if ($flag_count{$SIDELOBE_BIT});
print(STDERR "\n\t\t\tWAKE_BIT = $flag_count{$WAKE_BIT}")
if ($flag_count{$WAKE_BIT});
print(STDERR "\n\t\t\tPPI_BIT = $flag_count{$PPI_BIT}")
if ($flag_count{$PPI_BIT});
printf(STDERR "\n\t\t\tW_CTD_BIT = $flag_count{$W_CTD_BIT} (%d ensembles)",
$flag_count{$W_CTD_BIT}/$LADCP{N_BINS})
if ($flag_count{$W_CTD_BIT});
printf(STDERR "\n\t\t\tTILT_BIT = $flag_count{$TILT_BIT} (%d ensembles)",
$flag_count{$TILT_BIT}/$LADCP{N_BINS})
if ($flag_count{$TILT_BIT});
printf(STDERR "\n\t\t\tDELTA_TILT_BIT = $flag_count{$DELTA_TILT_BIT} (%d ensembles)",
$flag_count{$DELTA_TILT_BIT}/$LADCP{N_BINS})
if ($flag_count{$DELTA_TILT_BIT});
printf(STDERR "\n\t\t\tMISSING_CTD_DATA_BIT = $flag_count{$MISSING_CTD_DATA_BIT} (%d ensembles)",
$flag_count{$MISSING_CTD_DATA_BIT}/$LADCP{N_BINS})
if ($flag_count{$MISSING_CTD_DATA_BIT});
}
}
1;