new file mode 100644
--- /dev/null
+++ b/LADCPproc.UHcode
@@ -0,0 +1,435 @@
+#======================================================================
+# 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 Oct 26 14:40:39 2010
+# (c) 2010 A.M. Thurnherr & E. Firing
+# uE-Info: 431 79 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# PERLified functions from Eric's [merge.c]; with mods
+
+# NOTES:
+# - velocity integration removed
+# - percent-good flag removed (no point for single-ping ensembles)
+# - need the following per-ensemble fields from previous steps:
+# CTD_DEPTH
+# CTD_SVEL
+# - negative depths are not allowed (should not happen given CTD_DEPTH)
+
+# WEIRDNESSES:
+# - 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 [CORRECTED]
+# - $tilt calculation is wrong
+
+# 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
+
+#======================================================================
+# 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(%flag_count);
+
+# apparently, in Eric's code, DISTANCE_TO_BIN1_CENTER is not sound-speed corrected
+sub dzToBin($$)
+{
+ my($ens,$bin) = @_;
+ my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
+ return $sscorr * ($LADCP{DISTANCE_TO_BIN1_CENTER} + $bin*$LADCP{BIN_LENGTH});
+}
+
+sub depthOfBin($$)
+{
+ my($ens,$bin) = @_;
+ return $LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} ?
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH} - &dzToBin($ens,$bin) :
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH} + &dzToBin($ens,$bin);
+}
+
+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 ($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;
+ }
+
+ for (my($beam)=0; $beam<=3; $beam++) {
+ if ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] < $min_cor) {
+ $flag_count{$CORREL_BIT}++;
+ $edit_flags[$ens][$bin] |= $CORREL_BIT;
+ }
+ }
+
+ if ($bin < $shbin_start-1 || $bin >= $shbin_end) {
+ $edit_flags[$ens][$bin] |= $BADVEL_BIT;
+ $flag_count{$BADVEL_BIT}++;
+ }
+ } # for ($bin=0...
+}
+
+## The following is for editing out the second bottom bounce.
+sub set_PPI_flags($$)
+{
+ my($ens,$De) = @_;
+
+ 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 ($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} = 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 ($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;
+ } # get rid ensembles after large rotation
+ if (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 $opt_p && ($clip_margin > 0); # PPI editing is off by default
+ }
+}
+
+#======================================================================
+# CALCULATE VELOCITY SHEAR
+# - output in @ush_mu,@vsh_mu,@wsh_mu,@ush_sig,@vsh_sig,@wsh_sig
+# - @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh are defined "local" in calc_shear
+#======================================================================
+
+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 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 ($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);
+ }
+
+ 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[$bsgi]) > $max_shdev ||
+ abs($vsh_dev[$i]/$vsh_sig[$bsgi]) > $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
+
+ my(@ush_vals,@vsh_vals,@wsh_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);
+
+ 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]);
+ }
+ } # $ens loop
+
+ @ush_mu = @vsh_mu = @wsh_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);
+
+ $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];
+ }
+ $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;
+ }
+
+ 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}");
+ print(STDERR "\n\t\t\tERRVEL_BIT = $flag_count{$ERRVEL_BIT}");
+ print(STDERR "\n\t\t\tCORREL_BIT = $flag_count{$CORREL_BIT}");
+ print(STDERR "\n\t\t\tW_OUTLIER_BIT = $flag_count{$W_OUTLIER_BIT}");
+ print(STDERR "\n\t\t\tSHEAR_BIT = $flag_count{$SHEAR_BIT}");
+ print(STDERR "\n\t\t\tSIDELOBE_BIT = $flag_count{$SIDELOBE_BIT}");
+ print(STDERR "\n\t\t\tWAKE_BIT = $flag_count{$WAKE_BIT}");
+ print(STDERR "\n\t\t\tPPI_BIT = $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});
+ printf(STDERR "\n\t\t\tTILT_BIT = $flag_count{$TILT_BIT} (%d ensembles)",
+ $flag_count{$TILT_BIT}/$LADCP{N_BINS});
+ printf(STDERR "\n\t\t\tDELTA_TILT_BIT = $flag_count{$DELTA_TILT_BIT} (%d ensembles)",
+ $flag_count{$DELTA_TILT_BIT}/$LADCP{N_BINS});
+ }
+}
+
+1;