--- 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: Sat Jun 7 10:43:54 2014
+# dlm: Fri Mar 6 15:51:49 2015
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 59 43 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 8 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: Tue Aug 5 14:39:26 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 569 0 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 410 82 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
# NOTES:
@@ -98,7 +98,12 @@
# 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
+# - added -v to allow calculation of package velocity
+# Aug 5, 2014: - BUG: (bad one, too): ref_lr_w called from mk_prof had edited some of the
+# horizontal velocity data, which were nevertheless used later on!!!
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -113,7 +118,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";
@@ -123,7 +128,7 @@
$antsSummary = "$version -- process LADCP data to get shear, time series";
$antsParseHeader = 0;
-&antsUsage('24a:b:c:df:g:i:kl:n:o:p:rs:t:u:w:z',2,
+&antsUsage('24a:b:c:df:g:i:kl:n:o:p:rs:t:u:v:w:z',2,
'[use -2)dary CTD sensor pair]',
'[require -4)-beam LADCP solutions]',
'[use -r)DI bottom-track data]',
@@ -133,6 +138,7 @@
'[-i)nitial LADCP time lag <guestimate>]',
'[-l)ag LADCP <by>] [auto-lag -w)indow <size[120s]>] [-n) <auto-lag windows[20]]',
'[correct echo amplitude -u)sing D[eines99]|V[isbeck04]|T[hurnherr11]|n[ocorr]',
+ '[ocean -v)elocity <file> for calculating package velocity]',
'[-d)iagnostic screen output] [-z)oom through problems]',
'output: [shear-p)rofile <file>] [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
' [-a)coustic backscatter <dts-file] [bottom-trac-k) profs]',
@@ -169,7 +175,7 @@
croak("$0: cannot decode -c $opt_c\n")
unless numberp($CC_offset) && numberp($CC_cos_fac) && numberp($CC_sin_fac);
}
-
+
$LADCP_file = &antsFileArg();
$CTD_file = &antsFileArg();
@@ -219,6 +225,34 @@
&antsAddParams('LADCPproc::vertical_resolution',$GRID_DZ);
#----------------------------------------------------------------------
+# Step 2a: Load Ocean Velocity Profile
+#----------------------------------------------------------------------
+
+if (defined($opt_v)) { # load velocity profile
+ print(STDERR "Reading ocean-velocity profile ($opt_v)...");
+ if (open(OVF,$opt_v)) {
+ @ovl = &antsFileLayout(OVF);
+ $ovdF = &localFnr('depth',@ovl);
+ $ovuF = &localFnr('u',@ovl);
+ $ovvF = &localFnr('v',@ovl);
+ my(@ov) = &antsFileIn(OVF);
+ my($first_depth) = $ov[$ovdF]; my($last_depth);
+ do {
+ push(@ovu,$ov[$ovuF]);
+ push(@ovv,$ov[$ovvF]);
+ $last_depth = $ov[$ovdF];
+ } while (@ov = &antsFileIn(OVF));
+ close(OVF);
+ croak("$opt_v: incompatible depth grid\n")
+ unless (($last_depth-$first_depth) == $#ovu*$GRID_DZ);
+ printf(STDERR "\n\t%d velocities",scalar(@ovv)) if ($opt_d);
+ } else {
+ printf(STDERR "\n\n\t\tWARNING: $opt_v: $!\n");
+ }
+ print(STDERR "\n");
+}
+
+#----------------------------------------------------------------------
# Step 3: Read CTD data
#----------------------------------------------------------------------
@@ -326,6 +360,7 @@
@{$LADCP{ENSEMBLE}[$ens]->{BEAM_VEL}[$bin]} = @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]};
@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
velInstrumentToEarth(\%LADCP,$ens,velBeamToInstrument(\%LADCP,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));
+# print(STDERR "<$ens,$bin>: @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}\n");
@{$LADCP{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$bin]} = # fake it to fool ref_lr_w
(0,0,0,defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]) ? 100 : 0);
}
@@ -371,8 +406,8 @@
print(STDERR "\t\tconstructing profile time series...")
if ($opt_d);
-($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
- mk_prof(\%LADCP,0,undef,1,6,70,0.1,$LADCP_max_gap);
+($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) = # NB: chose parameters to avoid editing by
+ mk_prof(\%LADCP,0,undef,1,6,0,0.1,$LADCP_max_gap,0); # ref_lr_w
croak("\n$LADCP_file: no good ensembles found\n")
unless defined($LADCP_start);
@@ -557,6 +592,8 @@
print(STDERR "Calculating shear profiles & producing time-depth-series (.tds) output...");
@antsNewLayout = ('ens','elapsed','CTD_depth','downcast','CTD_lat','CTD_lon','depth','u_z','v_z','w_z');
+push(@antsNewLayout,'pkg_u','pkg_v')
+ if defined($opt_v);
#--------------------------------------------------------------------------------
# callback routine to output .tds data, called once each for down-/upcasts after
@@ -578,23 +615,31 @@
'LADCPproc::min_depth',depthOfGI($mingi),'LADCPproc::max_depth',depthOfGI($#ens_vals));
for (my($gi)=0; $gi<@ush_vals; $gi++) {
for (my($i)=0; $i<@{$ush_vals[$gi]}; $i++){
- &antsOut($ens_vals[$gi][$i],
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
- depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ my(@out) = ($ens_vals[$gi][$i],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
+ depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ push(@out,$LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$U],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$V])
+ if defined($opt_v);
+ &antsOut(@out);
}
}
} else {
for (my($gi)=$#ush_vals; $gi>=0; $gi--) {
for (my($i)=0; $i<@{$ush_vals[$gi]}; $i++) {
- &antsOut($ens_vals[$gi][$i],
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
- $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
- depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ my(@out) = ($ens_vals[$gi][$i],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
+ depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ push(@out,$LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$U],
+ $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{PACKAGE_VELOCITY}[$V])
+ if defined($opt_v);
+ &antsOut(@out);
}
}
}
--- a/LADCPproc.BT
+++ b/LADCPproc.BT
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . B T
# doc: Wed Oct 20 21:05:37 2010
-# dlm: Thu Sep 19 13:17:41 2013
+# dlm: Tue Aug 5 14:38:08 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 119 96 NIL 0 0 72 10 2 4 NIL ofnI
+# uE-Info: 205 0 NIL 0 0 72 10 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -15,6 +15,10 @@
# - changed from echo amplitude to Sv
# May 16, 2012: - added support for -r)DI BT data
# Sep 19, 2013: - added support for $BT_range_method
+# Aug 5, 2014: - BUG: invalid velocities were used; see July 15 bug fix
+# in [LADCPproc.shearmethod]
+# - artifically removed BT profile data apparently below
+# seabed but that had passed previous tests
my($BEAM1) = 0;
my($BEAM2) = 1;
@@ -174,7 +178,7 @@
if ($opt_k) {
for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
- next if ($edit_flags[$ens][$bin]);
+ next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
printf(BTF "%d %d %d %f %f %f %f %f %f %f %f %f %f %f\n",
$LADCP{ENSEMBLE}[$ens]->{NUMBER},
depthOfBin($ens,$bin),$LADCP{ENSEMBLE}[$ens]->{DEPTH},
@@ -191,13 +195,14 @@
}
for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
- next if ($edit_flags[$ens][$bin]);
- my($gi) = depthOfBin($ens,$bin) / $GRID_DZ;
+ next if ($edit_flags[$ens][$bin] || $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]==0);
+ my($dob) = depthOfBin($ens,$bin);
+ next if ($dob > $water_depth);
+ my($gi) = int($dob / $GRID_DZ);
push(@{$BTu_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$CTD_u);
push(@{$BTv_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V]-$CTD_v);
push(@{$BTw_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]-$CTD_w);
}
-
}
sub getBTprof($$)
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.backscatter
+++ b/LADCPproc.backscatter
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . B A C K S C A T T E R
# doc: Wed Oct 20 13:02:27 2010
-# dlm: Sun May 25 21:38:44 2014
+# dlm: Fri Mar 6 15:52:50 2015
# (c) 2010 A.M. Thurnherr
-# uE-Info: 26 45 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 32 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -24,6 +24,11 @@
# Mar 21, 2014: - adapted to new [LADCPproc.utils]
# Mar 27, 2014: - adapted to depthOfBinAlongBeam()
# May 25, 2014: - made search_below integer
+# Jul 27, 2014: - moved depthOfGI() to [LADCPproc.utils]
+# Aug 5, 2014: - BUG: find_backscatter_seabed() discarded everything if
+# LADCP bin 1 had the backscatter max at the edge of
+# the search domain, which could easily happen when
+# the bottom stop was a long way from the seabed
my($BEAM1) = 0;
my($BEAM2) = 1;
@@ -244,8 +249,6 @@
} # for $end
} # sub
-sub depthOfGI($) { return $_[0]*$GRID_DZ + $GRID_DZ/2; } # depth corresponding to particular grid index
-
sub find_backscatter_seabed($)
{
my($water_depth) = @_;
@@ -253,7 +256,8 @@
my($search_below) = int(max(0,$water_depth-$BT_begin_search_above));
my($mdgi) = int($search_below/$GRID_DZ); # grid index to begin search
- print(STDERR "\n\t\tlooking for seabed below $search_below m (gi >= $mdgi)") if ($opt_d);
+ printf(STDERR "\n\t\tlooking for seabed below %d m (gi = [%d..%d])",$search_below,$mdgi,scalar(@nSv))
+ if ($opt_d);
print(STDERR "\n\t\tseabed-max grid indices:") if ($opt_d);
@@ -267,19 +271,13 @@
$min = $avg if ($avg < $min);
$max = $avg, $gimax = $gi if ($avg > $max);
}
- if ($gimax == $firstvalid) { # should be robust except maybe for huge bins
- printf(STDERR "\n\t\tdata from below-seabed bins (%d-%d) discarded",
- $bin+1,$LADCP{N_BINS}) if ($opt_d);
- last;
- }
- if ($max-$min>10 && $gimax!=$lastvalid) { # ignore boundary maxima & scatter
- printf(STDERR " %d",$gimax-$mdgi) if ($opt_d);
+ if ($max-$min>10 && $gimax!=$firstvalid && $gimax!=$lastvalid) { # ignore boundary maxima & scatter
+ printf(STDERR " %d",$gimax) if ($opt_d);
push(@wdepth_gi,$gimax);
}
}
return (depthOfGI(avg(@wdepth_gi)),stddev(@wdepth_gi)*$GRID_DZ);
-
}
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: Sun May 25 21:16:27 2014
+# dlm: Fri Mar 6 15:53:43 2015
# (c) 2010 A.M. Thurnherr
-# uE-Info: 31 72 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 33 0 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
# May 25, 2015: - added assertion to require numeric interpolated LADCP_w
# - BUG: interp_LADCP_w left gaps
# - added debug code to output bestLag input time series
@@ -198,7 +199,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,520 @@
+#======================================================================
+# 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;
--- a/LADCPproc.utils
+++ b/LADCPproc.utils
@@ -1,27 +1,23 @@
#======================================================================
# L A D C P P R O C . U T I L S
# doc: Fri Mar 21 15:16:59 2014
-# dlm: Thu Mar 27 19:36:29 2014
+# dlm: Sun Jul 27 17:00:57 2014
# (c) 2014 A.M. Thurnherr
-# uE-Info: 57 38 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 14 49 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Mar 21, 2014: - created
# - added rangeToBin()
# Mar 27, 2014: - added rangeToBinAlongBeam()
-
-sub rangeToBin($$)
-{
- 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}) / cos(rad($LADCP{BEAM_ANGLE}));
-}
+# Jul 27, 2014: - improved comments
+# - moved depthOfGI() here from [LADCPproc.backscatter]
#----------------------------------------------------------------------
-# - in the UH code
-# - the distance to the first bin is not soundspeed-corrected
-# - instrument tilt is not considered
+# calculate depth of particular bin in particular ensemble
+# in contrast to the original UH code:
+# - the distance to the first bin is soundspeed-corrected
+# - instrument tilt is considered
#----------------------------------------------------------------------
sub depthOfBin($$)
@@ -40,6 +36,24 @@
$LADCP{ENSEMBLE}[$ens]->{DEPTH} + &dzToBin($ens,$bin);
}
+#----------------------------------------------------------------------
+# calculate along-beam distance between transducer and center of
+# particular bin in particular ensemble
+# - used for acoustic backscatter correction
+#----------------------------------------------------------------------
+
+sub rangeToBin($$)
+{
+ 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}) / cos(rad($LADCP{BEAM_ANGLE}));
+}
+
+#----------------------------------------------------------------------------
+# calculate depth of particular bin of particular beam in particular ensemble
+# - used to map acoustic backscatter of different beams correctly when
+# instrument tilt is large
+#----------------------------------------------------------------------------
sub depthOfBinAlongBeam($$$)
{
@@ -67,6 +81,21 @@
$LADCP{ENSEMBLE}[$ens]->{DEPTH} + &dzToBinAlongBeam($ens,$bin,$beam);
}
+#---------------------------------------------------------------------------
+# return center depth corresponding to particular grid index (in shear grid)
+#---------------------------------------------------------------------------
+
+sub depthOfGI($) { return $_[0]*$GRID_DZ + $GRID_DZ/2; } # depth corresponding to particular grid index
+
+#----------------------------------------------------------------------
+# return ocean velocity (u,v,w) at a given depth
+#----------------------------------------------------------------------
+
+sub oceanVel($)
+{
+
+}
+
1;