--- a/HISTORY
+++ b/HISTORY
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Tue May 15 18:04:39 2012
- dlm: Tue Aug 6 20:59:38 2013
+ dlm: Sun Jul 27 16:21:46 2014
(c) 2012 A.M. Thurnherr
- uE-Info: 61 37 NIL 0 0 72 3 2 8 NIL ofnI
+ uE-Info: 61 7 NIL 0 0 72 3 2 8 NIL ofnI
======================================================================
May 15, 2012:
@@ -58,5 +58,5 @@
Aug 6, 2013:
- V1.2 [.hg/hgrc] [LADCPproc.version]
- - but fixed in [LADCPproc.defaults]
+ - bug fixed in [LADCPproc.defaults]
--- a/LADCPintsh
+++ b/LADCPintsh
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P I N T S H
# doc: Thu Oct 14 21:22:50 2010
-# dlm: Thu Mar 20 12:14:53 2014
+# dlm: Tue Jul 15 09:14:27 2014
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 266 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 526 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
$antsSummary = 'integrate LADCP shear';
--- a/LADCPproc
+++ b/LADCPproc
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P P R O C
# doc: Thu Sep 16 20:36:10 2010
-# dlm: Tue May 20 10:46:41 2014
+# dlm: Sun Jul 27 16:23:50 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 569 0 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 103 69 NIL 0 0 72 11 2 4 NIL ofnI
#======================================================================
# NOTES:
@@ -98,7 +98,9 @@
# Mar 24, 2014: - added $pitch_offset, $roll_offset
# Mar 27, 2014: - BUG: Sv output did not reflect true bin depths
# Apr 25, 2014: - added LADCP_errvel to -t output
-# May 20, 2014: - merged laptop with whoosher versions (folded in Feb 22 change)
+# May 20, 2014: - merged laptop with whoosher versions (folded in Feb 22 change)
+# Jul 27, 2014: - renamed LADCPproc.UHcode to LADCPproc.shearmethod, because the code has
+# diverged more and more from the UH implementation
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -113,7 +115,7 @@
require "$LADCPPROC/LADCPproc.bestLag";
require "$LADCPPROC/LADCPproc.BT";
require "$LADCPPROC/LADCPproc.backscatter";
-require "$LADCPPROC/LADCPproc.UHcode";
+require "$LADCPPROC/LADCPproc.shearmethod";
require "$PERL_TOOLS/RDI_BB_Read.pl";
require "$PERL_TOOLS/RDI_Coords.pl";
require "$PERL_TOOLS/RDI_Utils.pl";
deleted file mode 100644
--- a/LADCPproc.UHcode
+++ /dev/null
@@ -1,495 +0,0 @@
-#======================================================================
-# L A D C P P R O C . U H C O D E
-# doc: Fri Sep 17 20:27:53 2010
-# dlm: Tue May 20 10:47:15 2014
-# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 404 0 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]
-
-#======================================================================
-# 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]);
- $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
- if (!$edit_flags[$ens][$bin]) {
- $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]);
-
- ## 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]);
-
- 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
- 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)
-# - sets @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh for a given ensemble
-#----------------------------------------------------------------------
-
-sub uv_to_shear($)
-{
- my($ens) = @_;
- my($nvel) = 0;
-
- @sh_i0 = @sh_i1 = ();
- for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
- next if ($edit_flags[$ens][$bin]);
- $nvel++;
- push(@sh_i1,$bin) if (@sh_i0);
- push(@sh_i0,$bin) if ($bin < $LADCP{N_BINS}-1);
- }
-
- @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]);
- 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);
- if ($edit_shear) {
- set_shear_flag($ens);
- $nvel += uv_to_shear($ens);
- $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;
--- a/LADCPproc.bestLag
+++ b/LADCPproc.bestLag
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . B E S T L A G
# doc: Tue Sep 28 21:58:48 2010
-# dlm: Wed Mar 19 21:33:30 2014
+# dlm: Sat Jul 19 18:32:26 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 28 45 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 29 44 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
@@ -26,6 +26,7 @@
# Oct 19, 2012: - BUG: opt_i had wrong sign!
# Jun 25, 2013: - adapted to :: %PARAM convention
# Mar 19, 2014: - moved %PARAM to LADCPproc
+# Jul 19, 2014: - made lagging obey -z)oom
sub interp_LADCP_w($$)
{
@@ -181,7 +182,7 @@
$best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
}
croak("\n$0: cannot determine a valid lag\n")
- unless ($nBest{$best_lag} > $opt_n/3);
+ unless ($opt_z || $nBest{$best_lag}>$opt_n/3);
print(STDERR "\n\n\t\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")
if ($nBest{$best_lag} < $opt_n/2);
new file mode 100644
--- /dev/null
+++ b/LADCPproc.shearmethod
@@ -0,0 +1,507 @@
+#======================================================================
+# 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 16:25:56 2014
+# (c) 2010 A.M. Thurnherr & E. Firing
+# uE-Info: 314 0 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
+
+#======================================================================
+# 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)
+# - sets @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh for a given ensemble
+#----------------------------------------------------------------------
+
+sub uv_to_shear($)
+{
+ my($ens) = @_;
+ my($nvel) = 0;
+
+ @sh_i0 = @sh_i1 = ();
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ 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);
+ }
+
+ @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);
+ if ($edit_shear) {
+ set_shear_flag($ens);
+ $nvel += uv_to_shear($ens);
+ $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;