--- a/LADCP_w	Tue Oct 11 17:05:36 2011 -0400
+++ b/LADCP_w	Sat Oct 15 23:23:57 2011 -0400
@@ -2,13 +2,12 @@
 #======================================================================
 #                    L A D C P _ W 
 #                    doc: Fri Dec 17 18:11:13 2010
-#                    dlm: Tue Oct 11 14:11:31 2011
+#                    dlm: Sat Oct 15 22:44:29 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 539 148 NIL 0 29 72 2 2 4 NIL ofnI
+#                    uE-Info: 810 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # TODO:
-#	- move (some? most? all?) data editing after profile begin/end
 #	  detection so that editing stats make sense
 #	- make timelagging work for short casts (make sure 10% is not more than half window size)
 #	- own seabed detection (P403)
@@ -79,14 +78,18 @@
 #	Oct 10, 2011: - BUG: LADCP_w in .tds output had not been soundspeed corrected
 #				  - removed -x
 #				  - added "false-positives" editing filter
-
-# PROCESSING PARAMETER FILE
-#	- # is comment character
-#	- invalid entries ignored
-# 	- valid entries begin with <ADCP-basename|-filename>: (NO INITIAL SPACE ALLOWED)
-#	- remainder of line is added to usage before ADCP file and LADCP_w is restarted
-#	- no argument with spaces allowed!
-#	- -0 suppresses acting on -l
+#	Oct 11, 2011: - adapted to the new parameter file, removing a few options
+#	Oct 12, 2011: - added surface-layer editing
+#				  - made %min/max_depth/elapsed more accurate
+#	Oct 13, 2011: - fiddled
+#	Oct 14, 2011: - renamed .prof fields from .N to .nsamp
+#				  - added %output_basename
+#				  - removed all ">" from open to allow use in pipelines
+#				  - made <run-label> argument optional
+#				  - replaced stdout output by $w_out
+#				  - renamed _out to out_; output_basename to out_basename
+#	Oct 15, 2011: - added editWOutliers()
+#				  - added step to remove single-ping noise
 
 # CTD REQUIREMENTS
 #	- elapsed		elapsed seconds; see note below
@@ -125,12 +128,6 @@
 #		1) 004DL.prof dc_w depth
 #		2) bindata -Sdowncast:1 -Fw.median,depth -n 20 depth 0 1000 10 004DL.w
 
-# TIME LAGGING
-#	- occasionally, the time lagging algorithm fails, in particular if the CTD is turned on
-#	  some time after the package enters the water
-#	- in this case, an initial guess can be provided with -i
-#	- e.g. plot 'CTD/24Hz/054.1Hz elapsed w','LADCP/raw/054UL.prof =$elapsed-850 w' => -i -850
-
 # VELOCITY AMBIGUITY ERRORS
 #	- quite extensive tests with DIMES US2 station 146, which has a lot of
 #	  ambiguity velocity errors, reveal that -m 1 catches those errors
@@ -138,14 +135,6 @@
 #	- even when the errors are not filtered with -m 1, they do not
 #	  affect the w profiles, as long as the median bin values are used
 
-# SCREEN LOGGING
-#	- there are 4 verbosity levels, selected by -v
-#		0:	only print errors
-#		1:	default, UNIX-like (warnings and info messages that are not produced for every cast)
-#		2:	progress messages and useful information
-#		>2:	debug messges
-#	- the most useful ones of these are 1 & 2
-
 ($WCALC) 	  = ($0                =~ m{^(.*)/[^/]*$});
 ($ANTS) 	  = (`which list`      =~ m{^(.*)/[^/]*$});
 ($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -162,16 +151,18 @@
 require "$PERL_TOOLS/RDI_BB_Read.pl";
 require "$PERL_TOOLS/RDI_Coords.pl";
 
+#------
+# Usage
+#------
+
 @ARGS = @ARGV;														# save opts & arguments
 
 $antsParseHeader = 0;
-&antsUsage('03:4a:b:c:d:e:f:g:h:i:k:l:m:n:o:p:r:s:t:uv:w:z:',1,
-	'[-l)oad processing parameters from <file>]',
+&antsUsage('3:4a:b:c:e:g:h:i:k:m:n:o:qs:t:uv:w:x:',1,
 	'[-v)erbosity <level[1]>]',
-    '[require -4)-beam solutions] [-d)iscard vels from <beam>]',
-	'[-r)ef-layer <bin[2],bin[6]>]',
+	'[-q)uick (no single-ping denoising)]',
+    '[require -4)-beam solutions]',
 	'[-c)orrelation <min[70]>] [-t)ilt <max[10deg]> [-e)rr-vel <max[0.1m/s]>]',
-	'[-z)ap (truncate) range by <n> bins]',
 	'[-h water <depth>]',
 	'[max LADCP time-series -g)ap <length[60s]>]',
 	'[-m)ax vertical <velocity[1m/s]>',
@@ -180,80 +171,36 @@
 	'[calculate -n) <lags,lags[10,100]>] [lag -w)indow <sz,sz[240s,20s]>]',
 	'[require top-3) lags to account for <frac[0.6]> of all]',
 	'[valid LADCP -b)ins <bin[2],bin[*]>',
-	'[-o)utput bin <resolution[50m]>] [-k) require <min[20]> samples]',
-	'[-f) write time-series <file>] [output binned -p)rofile to <file>]',
-	'<LADCP-file> [CTD-file]');
-
-&antsCardOpt(\$opt_v,1);											# suppress regular info
-
-$RDI_Coords::minValidVels = 4 if ($opt_4);							# suppress 3-beam solutions
-&antsCardOpt($opt_d);
-croak("$0: -d must be in the range 1-4, not $opt_d\n")
-	if defined($opt_d) && ($opt_d<1 || $opt_d>4);
-
-&antsFloatOpt(\$opt_c,70);											# min correlation
-&antsFloatOpt(\$opt_t,15);											# max tilt (pitch/roll)
-&antsFloatOpt(\$opt_e,0.1);											# max err vel
-&antsFloatOpt(\$opt_z,0);											# truncate final valid velocities
-&antsFloatOpt(\$opt_g,60);											# max LADCP gap length
-&antsFloatOpt(\$opt_m,1);											# max allowed vertical velocity
-
-&antsFloatOpt(\$opt_a,0);											# CTD depth adjustment
+	'[-o)utput bin <resolution[10m]>] [-k) require <min[20]> samples]',
+	'[e-x)ecute <perl-expr>]',
+	'<station-id> [run-label]');
 
 &antsFloatOpt($opt_i);												# externally supplied lag
 &antsUsageError() if ($opt_u && !defined($opt_i));
 
 &antsCardOpt(\$opt_s,0);											# skip # initial ensembles
 
-$opt_n = '10,100' unless defined($opt_n);							# number of time lags to carry out
-$opt_w = '240,20' unless defined($opt_w);							# time lag search window (full width)
-&antsFloatOpt(\$opt_3,0.6);
-
-&antsFloatOpt(\$opt_o,50);											# output bin size
-&antsCardOpt(\$opt_k,20);											# min samples required
-
-$opt_r = '2,6' unless defined($opt_r);								# reference layer bins for w for time matching
-$opt_b = '2,*' unless defined($opt_b);								# bins to use in w calculations
+$STN = &antsCardArg();												# station id
+$RUN = $ARGV[0],shift if (@ARGV);									# run label
 
-@n_lags = split(',',$opt_n);										# decode -n
-croak("$0: cannot decode -n $opt_n\n")
-	unless numberp($n_lags[0]) && numberp($n_lags[1]);
+#-----------------------------
+# Handle Processing Parameters
+#-----------------------------
 
-@w_size = split(',',$opt_w);										# decode -w
-croak("$0: cannot decode -w $opt_w\n")
-	unless numberp($w_size[0]) && numberp($w_size[1]);
-
-($refLr_firstBin,$refLr_lastBin) = split(',',$opt_r);				# decode -r
-croak("$0: cannot decode -r $opt_r\n")
-	unless numberp($refLr_firstBin) && numberp($refLr_lastBin);
+require "$WCALC/defaults.pl";										# load default/option parameters
+require "$processing_param_file";									# load processing parameters
+eval($opt_x) if defined($opt_x);									# eval cmd-line expression to override anything
 
-($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b);				# decode -b
-croak("$0: cannot decode -b $opt_b\n")
-	unless (numberp($LADCP_firstBin) &&
-			($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));
-
-$LADCP_file = &antsFileArg();
-($LADCP_basename) = ($LADCP_file =~ m{^.*/([^/]*)$});
+croak("$0: \$out_basename undefined\n")								# all plotting routines use this
+	unless defined($out_basename);
+&antsAddParams('out_basename',$out_basename);
 
-if (defined($opt_l) && !defined($opt_0)) {							# load per-cast processing parameters
-	my(@cast_params);
-	open(PF,$opt_l) || croak("$opt_l: $!\n");
-	while (<PF>) {
-		s/#.*//;
-		@cast_params=split(/\s+/),last if /^$LADCP_basename:/ || /^$LADCP_filename:/;
-	}
-	close(PF);
-
-	if (@cast_params) {												# found valid entry
-		if ($ARGS[$#ARGS] eq $LADCP_file) {							# CTD data on <stdin>
-			exec($0,@ARGS[0..$#ARGS-1],@cast_params[1..$#cast_params],'-0',$ARGS[$#ARGS]);
-			die("exec($0,@ARGS[0..$#ARGS-1],@cast_params[1..$#cast_params],'-0',$ARGS[$#ARGS]) failed!\n");
-    	} else {													# CTD file specified on cmdline
-			exec($0,@ARGS[0..$#ARGS-2],@cast_params[1..$#cast_params],'-0',$LADCP_file,$ARGS[$#ARGS]);
-			die("exec($0,@ARGS[0..$#ARGS-2],@cast_params[1..$#cast_params],'-0',$LADCP_file,$ARGS[$#ARGS]) failed!\n");
-    	}
-    }
-}
+$out_w 	    	= ">$out_w" 		 unless ($out_w =~ /^$|^\s*\|/);	# make sure O_CREAT is set unless it's a pipe
+$out_TL 	    = ">$out_TL" 		 unless ($out_TL =~ /^$|^\s*\|/);
+$out_TLhist    = ">$out_TLhist" 	 unless ($out_TLhist =~ /^$|^\s*\|/);
+$out_LADCPtis  = ">$out_LADCPtis"  unless ($out_LADCPtis =~ /^$|^\s*\|/);
+$out_profile    = ">$out_profile" 	 unless ($out_profile =~ /^$|^\s*\|/);
+$out_timeseries = ">$out_timeseries" unless ($out_timeseries =~ /^$|^\s*\|/);
 
 #----------------------------------------------------------------------
 # Screen Logging
@@ -264,32 +211,32 @@
 #----------------------------------------------------------------------
 
 sub progress(@)
-{  printf(STDERR @_) if ($opt_v > 1); }
+{  printf(STDERR @_) if ($verbosity_level > 1); }
 
 sub info(@)
 {
-	print(STDERR ($opt_v > 1) ? "\t" : "$LADCP_file: ");
-	printf(STDERR @_)  if ($opt_v > 0);
+	print(STDERR ($verbosity_level > 1) ? "\t" : "$LADCP_file: ");
+	printf(STDERR @_)  if ($verbosity_level > 0);
 }
 	
 sub warning(@)
 {
 	my($lvl,@msg) = @_;
-	return if ($opt_v == 0);
-	print(STDERR "\n") if ($opt_v > 1);
+	return if ($verbosity_level == 0);
+	print(STDERR "\n") if ($verbosity_level > 1);
 	print(STDERR "WARNING (level $lvl): ");
 	printf(STDERR @msg);
-	print(STDERR "\n") if ($opt_v > 1);
+	print(STDERR "\n") if ($verbosity_level > 1);
 }
 
 sub debugmsg(@)
-{ printf(STDERR @_) if ($opt_v > 2); }
+{ printf(STDERR @_) if ($verbosity_level > 2); }
 
 #----------------
 # Read LADCP data
 #----------------
 
-progress("Reading LADCP data ($LADCP_file)...\n");
+progress("Reading LADCP data from <$LADCP_file>...\n");
 readData($LADCP_file,\%LADCP);
 progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
 
@@ -331,23 +278,23 @@
 	$nvv = $cte = 0;
 	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
 		$nvv += countValidBeamVels($ens);
-		$cte += editCorr($ens,$opt_c);
-		$pte += editTilt($ens,$opt_t);
+		$cte += editCorr($ens,$min_correlation);
+		$pte += editTilt($ens,$max_tilt);
 	}
 	croak("$LADCP_file: no valid data\n") unless ($nvv > 0);
-	progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));
-	progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));
+	progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$min_correlation,$cte,round(100*$cte/$nvv));
+	progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$max_tilt,$pte,round(100*$pte/$nvv));
 } else {
 	progress("Editing velocity data...\n");
 	$nvv = $cte = 0;
 	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
 		$nvv += countValidBeamVels($ens);
-		$cte += editCorr_Earthcoords($ens,$opt_c);
-		$pte += editTilt($ens,$opt_t);
+		$cte += editCorr_Earthcoords($ens,$min_correlation);
+		$pte += editTilt($ens,$max_tilt);
 	}
 	croak("$LADCP_file: no valid data\n") unless ($nvv > 0);
-	progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));
-	progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));
+	progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$min_correlation,$cte,round(100*$cte/$nvv));
+	progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$max_tilt,$pte,round(100*$pte/$nvv));
 }
 
 #-------------------------------------------------------------------
@@ -362,13 +309,13 @@
 
 if ($LADCP{BEAM_COORDINATES}) {
 	progress("Calculating earth-coordinate velocities...\n");
-	progress("\tdiscarding velocities from beam $opt_d\n")
-		if defined($opt_d);
+	progress("\tdiscarding velocities from beam $discard_velocities_from_beam\n")
+		if ($discard_velocities_from_beam);
 	$nvw = 0;
 	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
 		for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
-			undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$opt_d-1])
-				if defined($opt_d);
+			undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$discard_velocities_from_beam-1])
+				if ($discard_velocities_from_beam);
 			($LADCP{ENSEMBLE}[$ens]->{U}[$bin],
 			 $LADCP{ENSEMBLE}[$ens]->{V}[$bin],
 			 $LADCP{ENSEMBLE}[$ens]->{W}[$bin],
@@ -413,23 +360,29 @@
 #----------------------------------------------
 # S1 STEP: Edit earth-coordinate -velocity data
 #	1) error-velocity threshold
+#	2) vertical-velocity outliers
+#	3) false-positives
+#	4) truncate range
 #----------------------------------------------
 
 progress("Editing earth-coordinate velocity data...\n");
 
-$ete = $zte = $fpte = 0;
+$evte = $trte = $fpte = $wote = 0;
 for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
-	$ete  += editErrVel($ens,$opt_e);
+	$evte += editErrVel($ens,$max_allowed_errvel);
+	$wote += editWOutliers($ens);
 	$fpte += editFalsePositives($ens);
-	$zte  += editTruncRange($ens,$opt_z) if ($opt_z > 0);
+	$trte += editTruncRange($ens,$truncate_farthest_valid_bins) if ($truncate_farthest_valid_bins > 0);
 }
 progress("\terror-velocity threshold (-e %.1f m/s): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
-	$opt_e,$ete,round(100*$ete/$nvw));
+	$max_allowed_errvel,$evte,round(100*$evte/$nvw));
+progress("\tvertical-velocity outliers: %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
+	$wote,round(100*$wote/$nvw));
 progress("\t\"false-positives\" filter: %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
 	$fpte,round(100*$fpte/$nvw));
 progress("\trange truncation (-z %d bins): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
-	$opt_z,$zte,round(100*$zte/$nvw))
-		if ($opt_z > 0);
+	$truncate_farthest_valid_bins,$trte,round(100*$trte/$nvw))
+		if ($truncate_farthest_valid_bins > 0);
 
 #----------------------------
 # Calculate LADCP time series
@@ -438,7 +391,7 @@
 progress("Calculating LADCP time-series...\n");
 
 ($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) =
-	calcLADCPts(\%LADCP,$opt_s,$refLr_firstBin,$refLr_lastBin,$opt_g);
+	calcLADCPts(\%LADCP,$opt_s,$refLr_firstBin,$refLr_lastBin,$max_LADCP_reflr_vel_gap);
 croak("$LADCP_file: no good ensembles\n")
 	unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0);
 
@@ -466,12 +419,29 @@
 						$LADCP{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);
 progress("\tMean ping interval: %.1f seconds\n",$LADCP{MEAN_DT});
 
+if (defined($out_LADCPtis)) {
+	progress("Writing LADCP-timeseries to <$out_LADCPtis>...\n");
+
+	@antsNewLayout = ('elapsed','reflr_w','reflr_w.stddev','reflr_w.nsamp','depth');
+	close(STDOUT);
+	open(STDOUT,"$out_LADCPtis") || croak("$out_LADCPtis: $!\n");
+	for (my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+		&antsOut($LADCP{ENSEMBLE}[$ens]->{ELAPSED},
+				 $LADCP{ENSEMBLE}[$ens]->{REFLR_W},
+				 $LADCP{ENSEMBLE}[$ens]->{REFLR_W_STDDEV},
+				 $LADCP{ENSEMBLE}[$ens]->{REFLR_W_NSAMP},
+				 $LADCP{ENSEMBLE}[$ens]->{DEPTH});
+	}
+	&antsOut('EOF'); close(STDOUT);
+}
+
 #--------------
 # Read CTD data
 #--------------
 
-progress("Reading CTD data...\n");
-croak("$0: no CTD data\n") unless (&antsIn());
+progress("Reading CTD data from <$CTD_file>...\n");
+open(STDIN,$CTD_file) || croak("$CTD_file: $!\n");
+croak("$CTD_file: no data\n") unless (&antsIn());
 undef($antsOldHeaders);
 ($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w');
 $CTD_temp = &fnrNoErr('temp');
@@ -482,12 +452,12 @@
 	croak("$0: cannot deal with non-numeric CTD elapsed time\n")
 		unless &antsNumbers($CTD_elapsed);
 	push(@{$CTD{ELAPSED}},$ants_[0][$CTD_elapsed]);
-	push(@{$CTD{DEPTH}},  $ants_[0][$CTD_depth]+$opt_a);
+	push(@{$CTD{DEPTH}},  $ants_[0][$CTD_depth]+$CTD_neg_press_offset);
 	push(@{$CTD{SVEL}},   $ants_[0][$CTD_svel]);
 	push(@{$CTD{W}},      $ants_[0][$CTD_w]);
 	push(@{$CTD{TEMP}},	  $ants_[0][$CTD_temp]) if defined($CTD_temp);
-	if ($ants_[0][$CTD_depth]+$opt_a > $CTD_maxdepth) {
-		$CTD_maxdepth = $ants_[0][$CTD_depth]+$opt_a;
+	if ($ants_[0][$CTD_depth]+$CTD_neg_press_offset > $CTD_maxdepth) {
+		$CTD_maxdepth = $ants_[0][$CTD_depth]+$CTD_neg_press_offset;
 		$CTD_atbottom = $#{$CTD{DEPTH}};
 	}
 } while (&antsIn());
@@ -536,8 +506,6 @@
 	$CTD{TIME_LAG} =
 		$CTD{ELAPSED}[$CTD_10pct_down] - $LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED};
 
-	print(STDERR "$CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0],  $CTD{ELAPSED}[$CTD_10pct_down],  $LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED}\n");
-	
 	progress("\telapsed(dz(CTD)=%.1fm) ~ elapsed(dz(LADCP)=%.1fm) + %.1fs\n",
 	    $CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH},$CTD{TIME_LAG});
 }
@@ -545,10 +513,10 @@
 if ($opt_u) {
 	progress("\tskipping time lagging (-u)\n");
 } else {
-	$CTD{TIME_LAG} = calc_lag($n_lags[0],$w_size[0],int(1/$CTD{DT}+0.5));
+	$CTD{TIME_LAG} = calc_lag($number_of_timelag_windows[0],$length_of_timelag_windows[0],int(1/$CTD{DT}+0.5));
 	progress("\telapsed(CTD) ~ elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG});
 	
-	$CTD{TIME_LAG} = calc_lag($n_lags[1],$w_size[1],1);
+	$CTD{TIME_LAG} = calc_lag($number_of_timelag_windows[1],$length_of_timelag_windows[1],1);
 	progress("\telapsed(CTD) = elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG});
 }
 
@@ -602,7 +570,7 @@
 		$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan];
 		$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan;
 		my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan];
-		if (abs($reflr_ocean_w) <= $opt_m) {
+		if (abs($reflr_ocean_w) <= $max_allowed_w) {
 			$sumWsq += &SQR($reflr_ocean_w);
 			$nWsq++;
 			if ($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} > 100 &&
@@ -686,13 +654,22 @@
 	progress("\t$nvrm velocities from $nerm ensembles removed\n");
 }
 
+#----------------------------------------------------------------------
+# Data Editing after LADCP and CTD data have been merged
+#	1) surface layer editing
+#----------------------------------------------------------------------
+
+progress("Removing data from instrument at surface...\n");
+$nerm = editSurfLayer($firstGoodEns,$lastGoodEns,$surface_layer_depth);
+progress("\t$nerm ensembles removed\n");
+
 #---------------------------------------------------------------------------
 # Depth-bin LADCP velocity data
 #
 # NOTES:
 #	1) ensemble and bin numbers are saved for maximum flexibility
 #	2) only ensemble/bins with valid vertical velocities are saved
-#	3) applying the full soundspeed correction to w is most likely pointless in
+#	3) applying the full soundspeed correction to w is pointless in
 #	   practice, but hey!, CPU cycles are cheap; [in a cast in the Gulf of Mexico
 #	   which has fairly pronounce soundspeed gradients, the max value of Kn
 #	   is 1.00004160558372, which gives rise to a correction of less than 0.2mm/s
@@ -708,7 +685,10 @@
 
 progress("\tdowncast...\n");
 for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) {						# downcast
-	next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+	unless (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})) {
+		$firstGoodEns++ if ($ens == $firstGoodEns);								# start has been edited away
+		next;
+	}
 	my(@bindepth) = calc_binDepths($ens);
 	for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
 		next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
@@ -721,7 +701,7 @@
 					 $bindepth[$bin]);
 		$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin] =
 			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] + $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}];
-		my($bi) = $bindepth[$bin]/$opt_o;
+		my($bi) = $bindepth[$bin]/$output_bin_size;
 		push(@{$DNCAST{ENSEMBLE}[$bi]},$ens);
 		push(@{$DNCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
 		push(@{$DNCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
@@ -742,8 +722,10 @@
 
 progress("\tupcast...\n");														# upcast
 
+my($realLastGoodEns);
 for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {		
 	next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+	$realLastGoodEns = $ens;
 	my(@bindepth) = calc_binDepths($ens);
 	for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
 		next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
@@ -756,7 +738,7 @@
 					 $bindepth[$bin]);
 		$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin] =
 			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] + $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}];
-		my($bi) = $bindepth[$bin]/$opt_o;
+		my($bi) = $bindepth[$bin]/$output_bin_size;
 		push(@{$UPCAST{ENSEMBLE}[$bi]},$ens);
 		push(@{$UPCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
 		push(@{$UPCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
@@ -774,9 +756,60 @@
 	$UPCAST{N_SAMP}[$bi]		= @{$UPCAST{W}[$bi]};
 }
 
-&antsAddParams('min_depth',$min_depth,'max_depth',$max_depth,
+&antsAddParams('min_depth',$min_depth,'max_depth',$max_depth,					# plot parameters
+			   'min_ens',$firstGoodEns,'max_ens',$realLastGoodEns,
 			   'min_elapsed',$LADCP{ENSEMBLE}[$firstGoodEns]->{CTD_ELAPSED},
-			   'max_elapsed',$LADCP{ENSEMBLE}[$lastGoodEns]->{CTD_ELAPSED});
+			   'max_elapsed',$LADCP{ENSEMBLE}[$realLastGoodEns]->{CTD_ELAPSED},
+			   'LADCP_bin_length',$LADCP{BIN_LENGTH});
+
+#------------------------------------------------------------------------------------------------------
+# remove single-ping noise
+#	- based on observation that without this step raw velocities tend to be strongly vertically striped
+#		=> residual CTD/LADCP mismatch errors & ADCP single-ping errors
+#	- remove ensemble-median residual from all vertical velocities; new residual calculated below
+#	- improves solution but can introduce errors in ensembles with strong outliers => median instead of
+#	  mean, and outlier filter above
+#------------------------------------------------------------------------------------------------------
+
+unless ($opt_q) {
+	progress("Removing single-ping residual noise...\n");
+	
+	for ($ens=$firstGoodEns; $ens<=$realLastGoodEns; $ens++) {
+		next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+		my(@residuals) = ();
+		my(@bindepth) = calc_binDepths($ens);
+		for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+			next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+			my($bi) = $bindepth[$bin]/$output_bin_size;
+			push(@residuals,
+				$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] -
+					(($ens < $LADCP_atbottom) ? $DNCAST{MEDIAN_W}[$bi]
+											  : $UPCAST{MEDIAN_W}[$bi]));
+		}
+		$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W} = median(@residuals);
+		for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+			next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] -=
+				$LADCP{ENSEMBLE}[$ens]->{MEDIAN_RESIDUAL_W};
+		}
+    }
+
+    progress("\tre-calculating profile mad...\n");
+    
+	for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) {						# bin data
+		for (my($i)=0; $i<@{$DNCAST{W}[$bi]}; $i++) {
+			$DNCAST{W}[$bi][$i] -= $LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
+		}
+		$DNCAST{MAD_W}[$bi] = mad2($DNCAST{MEDIAN_W}[$bi],@{$DNCAST{W}[$bi]});
+	}
+	for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
+		for (my($i)=0; $i<@{$UPCAST{W}[$bi]}; $i++) {
+			$UPCAST{W}[$bi][$i] -= $LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{MEDIAN_RESIDUAL_W};
+		}
+		$UPCAST{MAD_W}[$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
+	}
+
+} # unless ($opt_q);
 
 #--------------------------------------------------
 # Calculate BT-referenced vertical-velocity profile
@@ -784,12 +817,12 @@
 
 if (defined($water_depth)) {
 	progress("Calculating BT-referenced vertical velocities\n");
-	calc_BTprof($firstGoodEns,$lastGoodEns,$water_depth,$sig_water_depth);
+	calc_BTprof($firstGoodEns,$realLastGoodEns,$water_depth,$sig_water_depth);
 
 	my($sumSq) = my($n) = 0;
 	for (my($bi)=0; $bi<=$#{$BT{MEDIAN_W}}; $bi++) {
 		next unless defined($BT{MEDIAN_W}[$bi]);
-		next unless ($BT{N_SAMP}[$bi]>=$opt_k && $DNCAST{N_SAMP}[$bi]>=$opt_k && $UPCAST{N_SAMP}[$bi]>=$opt_k);
+		next unless ($BT{N_SAMP}[$bi]>=$min_w_nsamp && $DNCAST{N_SAMP}[$bi]>=$min_w_nsamp && $UPCAST{N_SAMP}[$bi]>=$min_w_nsamp);
 		$sumSq += ($BT{MEDIAN_W}[$bi] - $DNCAST{MEDIAN_W}[$bi]/2 - $UPCAST{MEDIAN_W}[$bi]/2)**2;
 		$n++;
 	}
@@ -806,122 +839,122 @@
 
 # NB: residual field is calculated with respect to down-/upcast medians in -o-size bins
 
-progress("Writing data...\n");
+if (defined($out_w)) {
+	progress("Writing vertical-velocity data to $out_w...\n");
 
-@antsNewLayout = ('ensemble','bin','elapsed','depth','CTD_depth','downcast',
-				  'w','residual','CTD_w','LADCP_w','errvel',
-				  'pitch','roll','tilt','heading','3_beam','svel');
+	@antsNewLayout = ('ensemble','bin','elapsed','depth','CTD_depth','downcast',
+					  'w','residual','CTD_w','LADCP_w','errvel',
+					  'pitch','roll','tilt','heading','3_beam','svel');
 
-for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) {						# downcast
-	next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
-	my(@bindepth) = calc_binDepths($ens);
-	for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
-		next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
-		my($bi) = $bindepth[$bin]/$opt_o;
-		&antsOut(
-			$ens,$bin,$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
-			$bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},1,
-			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
-			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $DNCAST{MEDIAN_W}[$bi],
-			$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
-			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
-			$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
-			$LADCP{ENSEMBLE}[$ens]->{PITCH},
-			$LADCP{ENSEMBLE}[$ens]->{ROLL},
-			$LADCP{ENSEMBLE}[$ens]->{TILT},
-			$LADCP{ENSEMBLE}[$ens]->{HEADING},
-			(defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + 			# only works for beam coords
-			 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
-			 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
-			 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
-			$CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
-		);
+	open(STDOUT,"$out_w") || croak("$out_w: $!\n");
+	for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) {						# downcast
+		next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+		my(@bindepth) = calc_binDepths($ens);
+		for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+			next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+			my($bi) = $bindepth[$bin]/$output_bin_size;
+			&antsOut(
+				$ens,$bin,$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+				$bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},1,
+				$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
+				$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $DNCAST{MEDIAN_W}[$bi],
+				$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+				$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
+				$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
+				$LADCP{ENSEMBLE}[$ens]->{PITCH},
+				$LADCP{ENSEMBLE}[$ens]->{ROLL},
+				$LADCP{ENSEMBLE}[$ens]->{TILT},
+				$LADCP{ENSEMBLE}[$ens]->{HEADING},
+				(defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + 			# only works for beam coords
+				 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
+				 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
+				 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
+				$CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+			);
+		}
 	}
-}
+	
+	for ($ens=$LADCP_atbottom; $ens<=$realLastGoodEns; $ens++) {						# upcast 
+		next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+		my(@bindepth) = calc_binDepths($ens);
+		for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+			next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+			my($bi) = $bindepth[$bin]/$output_bin_size;
+			&antsOut(
+				$ens,$bin,$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+				$bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},0,
+				$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
+				$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $UPCAST{MEDIAN_W}[$bi],
+				$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+				$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
+				$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
+				$LADCP{ENSEMBLE}[$ens]->{PITCH},
+				$LADCP{ENSEMBLE}[$ens]->{ROLL},
+				$LADCP{ENSEMBLE}[$ens]->{TILT},
+				$LADCP{ENSEMBLE}[$ens]->{HEADING},
+				(defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + 			# only works for beam coords
+				 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
+				 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
+				 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
+				$CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+			);
+		}
+	}
 
-for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {						# upcast 
-	next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
-	my(@bindepth) = calc_binDepths($ens);
-	for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
-		next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
-		my($bi) = $bindepth[$bin]/$opt_o;
-		&antsOut(
-			$ens,$bin,$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
-			$bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},0,
-			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
-			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $UPCAST{MEDIAN_W}[$bi],
-			$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
-			$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_W}[$bin],
-			$LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
-			$LADCP{ENSEMBLE}[$ens]->{PITCH},
-			$LADCP{ENSEMBLE}[$ens]->{ROLL},
-			$LADCP{ENSEMBLE}[$ens]->{TILT},
-			$LADCP{ENSEMBLE}[$ens]->{HEADING},
-			(defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + 			# only works for beam coords
-			 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
-			 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
-			 defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
-			$CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
-		);
-	}
+	&antsOut('EOF'); close(STDOUT);
 }
 	
 #---------------
 # Output profile
 #---------------
 
-if (defined($opt_p)) {
-	progress("Writing vertical-velocity profiles to <$opt_p>...\n");
+if (defined($out_profile)) {
+	progress("Writing vertical-velocity profiles to <$out_profile>...\n");
 
-	@antsNewLayout = ('depth','dc_depth','dc_elapsed','dc_w','dc_w.mad','dc_w.N',
-							  'uc_depth','uc_elapsed','uc_w','uc_w.mad','uc_w.N',
-							  'elapsed','w','w.mad','w.N',
-	                          'BT_w','BT_w.mad','BT_w.N');
+	@antsNewLayout = ('depth','dc_depth','dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
+							  'uc_depth','uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp',
+							  'elapsed','w','w.mad','w.nsamp',
+	                          'BT_w','BT_w.mad','BT_w.nsamp');
 
-	&antsOut('EOF');
-	
-	close(STDOUT);
-	open(STDOUT,">$opt_p") || croak("$opt_p: $!\n");
+	open(STDOUT,"$out_profile") || croak("$out_profile: $!\n");
 
 	for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) {
-		&antsOut(($bi+0.5)*$opt_o,					# nominal depth
+		&antsOut(($bi+0.5)*$output_bin_size,					# nominal depth
 				 $DNCAST{MEAN_DEPTH}[$bi],$DNCAST{MEAN_ELAPSED}[$bi],
-				 $DNCAST{N_SAMP}[$bi]>=$opt_k?$DNCAST{MEDIAN_W}[$bi]:nan,
+				 $DNCAST{N_SAMP}[$bi]>=$min_w_nsamp?$DNCAST{MEDIAN_W}[$bi]:nan,
 				 $DNCAST{MAD_W}[$bi],$DNCAST{N_SAMP}[$bi],
 				 $UPCAST{MEAN_DEPTH}[$bi],$UPCAST{MEAN_ELAPSED}[$bi],
-				 $UPCAST{N_SAMP}[$bi]>=$opt_k?$UPCAST{MEDIAN_W}[$bi]:nan,
+				 $UPCAST{N_SAMP}[$bi]>=$min_w_nsamp?$UPCAST{MEDIAN_W}[$bi]:nan,
 				 $UPCAST{MAD_W}[$bi],$UPCAST{N_SAMP}[$bi],
 				 $DNCAST{MEAN_ELAPSED}[$bi]/2+$UPCAST{MEAN_ELAPSED}[$bi]/2,
-				 $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$opt_k ?
+				 $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$min_w_nsamp ?
 					($DNCAST{MEDIAN_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MEDIAN_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
 					nan,
-				 $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$opt_k ?
+				 $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$min_w_nsamp ?
 					 ($DNCAST{MAD_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MAD_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
 					 nan,
 				 $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi],
-				 $BT{N_SAMP}[$bi]>=$opt_k?$BT{MEDIAN_W}[$bi]:nan,
+				 $BT{N_SAMP}[$bi]>=$min_w_nsamp?$BT{MEDIAN_W}[$bi]:nan,
 				 $BT{MAD_W}[$bi],$BT{N_SAMP}[$bi]
 		);
 	}
+	&antsOut('EOF'); close(STDOUT);
 }
 
 #--------------------------------------
 # write time-series output if requested
 #--------------------------------------
 
-if (defined($opt_f)) {
-	progress("Writing time-series data to <$opt_f>...\n");
+if (defined($out_timeseries)) {
+	progress("Writing time-series data to <$out_timeseries>...\n");
 
 	@antsNewLayout = ('ens','elapsed',
 					  'depth','xducer_sound_speed','pitch','gimbal_pitch','roll','tilt','heading',
 					  'CTD_w','LADCP_reflr_w','LADCP_reflr_w_err','LADCP_reflr_w_nosscorr',
 					  'reflr_ocean_w');
-	&antsOut('EOF');
-
-	close(STDOUT);
-	open(STDOUT,">$opt_f") || croak("$opt_f: $!\n");
+	open(STDOUT,"$out_timeseries") || croak("$out_timeseries: $!\n");
 	 
-	for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+	for ($ens=$firstGoodEns; $ens<=$realLastGoodEns; $ens++) {
 		my($reflr_oc_w) = defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W})
 			   		    ? $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]
 			   		    : undef;
@@ -941,8 +974,7 @@
 				 $reflr_w);
 	}
 
-	close(STDOUT);
-	undef($antsHeadersPrinted);
+	&antsOut('EOF'); close(STDOUT);
 }
 
 &antsExit();

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LWplot_TL	Sat Oct 15 23:23:57 2011 -0400
@@ -0,0 +1,44 @@
+#!/bin/sh
+#======================================================================
+#                    L W P L O T _ T L 
+#                    doc: Thu Oct 13 10:51:49 2011
+#                    dlm: Sat Oct 15 22:24:18 2011
+#                    (c) 2011 A.M. Thurnherr
+#                    uE-Info: 15 31 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Oct 13, 2011: - created
+#	Oct 14, 2011: - fiddled
+
+[ -f .gmtdefaults4 ] ||
+	gmtset	PAPER_MEDIA letter+ \
+			LABEL_FONT_SIZE 14 ANNOT_FONT_SIZE_PRIMARY 14 \
+			WANT_EURO_FONT true \
+	        PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+while [ -z "$fields" ]
+do
+	read line
+	[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
+	[ -z "$min_elapsed" ] && min_elapsed=`expr -- "$line" : '#ANTS#PARAMS#.*elapsed.min{\([^}]*\)}`
+	[ -z "$max_elapsed" ] && max_elapsed=`expr -- "$line" : '#ANTS#PARAMS#.*elapsed.max{\([^}]*\)}`
+	[ -z "$best_scan_offset" ] && best_scan_offset=`expr -- "$line" : '#ANTS#PARAMS#.*best_scan_offset{\([^}]*\)}`
+	[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{// -e s/}//`
+done
+
+R=-R`echo $min_elapsed/60|bc`/`echo $max_elapsed/60+1|bc`/-24/24
+U=-R0/1/0/1
+J=-JX10
+
+while read elapsed scan_offset mad
+do
+	echo `echo "scale=2;$elapsed/60" | bc` $scan_offset
+done | psxy -P -K $J $R \
+			-Bf1a15:"Elapsed Time [min]":/f1a5:"Best Offset [scans]":WeSn \
+		    -Sc0.1 -G0
+	
+echo 0.02 0.98 12 0 0 TL $out_basename | pstext -O -K $U $J
+
+echo -e "`echo $min_elapsed/60|bc` $best_scan_offset\n`echo $max_elapsed/60+1|bc` $best_scan_offset" \
+	| psxy -O $R $J

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LWplot_TLhist	Sat Oct 15 23:23:57 2011 -0400
@@ -0,0 +1,43 @@
+#!/bin/sh
+#======================================================================
+#                    L W P L O T _ T L H I S T 
+#                    doc: Thu Oct 13 10:51:49 2011
+#                    dlm: Sat Oct 15 22:24:31 2011
+#                    (c) 2011 A.M. Thurnherr
+#                    uE-Info: 14 31 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Oct 13, 2011: - created
+
+[ -f .gmtdefaults4 ] ||
+	gmtset	PAPER_MEDIA letter+ \
+			LABEL_FONT_SIZE 14 ANNOT_FONT_SIZE_PRIMARY 14 \
+			WANT_EURO_FONT true \
+	        PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+while [ -z "$fields" ]
+do
+	read line
+	[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
+	[ -z "$n_windows" ] && n_windows=`expr -- "$line" : '#ANTS#PARAMS#.*n_windows{\([^}]*\)}`
+	[ -z "$best_scan_offset" ] && best_scan_offset=`expr -- "$line" : '#ANTS#PARAMS#.*best_scan_offset{\([^}]*\)}`
+	[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{// -e s/}//`
+done
+
+R=-R-24/24/0/60
+U=-R0/1/0/1
+J=-JX10
+
+while read scan_offset nsamp mad
+do
+	echo $scan_offset `echo 100*$nsamp/$n_windows | bc`
+done | psxy -P -K $J $R \
+			-Bf1a10:"Best Offset [scans]":/f1a5:"Fraction [%]":WeSn \
+		    -Sb1u -G127
+	
+echo 0.02 0.98 12 0 0 TL $out_basename | pstext -O -K $J $U
+
+echo -e "$best_scan_offset 0\n$best_scan_offset 60" \
+	| psxy -O $J $R
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LWplot_prof	Sat Oct 15 23:23:57 2011 -0400
@@ -0,0 +1,118 @@
+#!/bin/sh
+#======================================================================
+#					 L W P L O T _ P R O F 
+#                    doc: Fri Oct 14 09:42:36 2011
+#                    dlm: Sat Oct 15 23:20:48 2011
+#                    (c) 2011 A.M. Thurnherr
+#                    uE-Info: 96 50 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Oct 14, 2011: - created
+#	Oct 15, 2011: - added nsamp
+
+#--------------------------------------------------
+# Usage
+#--------------------------------------------------
+
+USAGE="Usage: $0 [-s)ave data to <file>] [in-file]"
+
+args=`getopt s: $*`
+[ $? != 0 ] && {
+	echo $USAGE >&2
+	exit 1
+}
+
+set -- $args
+for opt
+do
+	case "$opt" in
+		-s) opt_s="$2";
+			rm -f "$opt_s";
+			shift; shift;;
+		--) shift; break;;
+	esac
+done
+
+[ $# -gt 1 ] && {
+	echo $USAGE >&2
+	exit 1
+}
+
+if [ $# -eq 1 ]
+then
+	if [ -n "$opt_s" ]
+	then
+		exec <$1 $0 -s $opt_s
+	else
+		exec <$1 $0
+	fi
+fi
+
+#--------------------------------------------------
+# Read Header Data
+#--------------------------------------------------
+
+TMPFILE=/tmp/$$.LWplot_prof
+
+[ -f .gmtdefaults4 ] ||
+	gmtset	PAPER_MEDIA letter+ \
+			LABEL_FONT_SIZE 14 ANNOT_FONT_SIZE_PRIMARY 14 \
+			WANT_EURO_FONT true \
+	        PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+[ -n "$opt_s" ] && rm -f "$opt_s"
+
+while [ -z "$fields" ]
+do
+	read line
+	[ -n "$opt_s" ] && echo "$line" >> "$opt_s"
+	[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
+	[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
+	[ -z "$max_depth" ] && max_depth=`expr -- "$line" : '#ANTS#PARAMS#.*max_depth{\([^}]*\)}`
+	[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+done
+
+if [ -n "$opt_s" ]
+then
+	tee -a "$opt_s" > $TMPFILE
+else
+	cat > $TMPFILE
+fi
+
+#--------------------------------------------------
+# Plot Data
+#--------------------------------------------------
+
+R=-R-0.07/0.35/$min_depth/$max_depth
+R2=-R-200/200/$min_depth/$max_depth
+U=-R0/1/0/1
+X=-JX10/-10
+
+echo -e "0 $min_depth\n0 $max_depth" | psxy -P -K $R $X
+echo -e "0.07 $min_depth\n0.07 $max_depth\n0.18 $max_depth\n0.18 $min_depth" | psxy -O -K $R $X -L -G200
+echo -e "50 $min_depth\n50 $max_depth" | psxy -O -K $R2 $X
+echo -e "100 $min_depth\n100 $max_depth" | psxy -O -K $R2 $X
+echo -e "150 $min_depth\n150 $max_depth" | psxy -O -K $R2 $X
+
+awk '{print $4, $1}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -W4/red
+awk '{print $9, $1}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -W4/blue
+awk '{print $16,$1}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -W4/green
+	
+awk '{print  $5,$1, $4}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gred
+awk '{print $10,$1, $9}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gblue
+awk '{print $17,$1,$16}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Ggreen
+
+awk '{print  $6,$1, $4}' $TMPFILE | grep -vi nan | psxy -O -K -N $R2 $X -W1/red
+awk '{print $11,$1, $9}' $TMPFILE | grep -vi nan | psxy -O -K -N $R2 $X -W1/blue
+awk '{print $18,$1,$16}' $TMPFILE | grep -vi nan | psxy -O -K -N $R2 $X -W1/green
+
+echo 0.02 0.02 12 0 0 TL $out_basename | pstext -O -K $U $X
+echo 0.6 0.98 12 0 0 BR m.a.d. | pstext -O -K $U $X
+
+psbasemap -O -K $R $X -Bf0.01a0.1:"Vertical Velocity [m/s]":/f10a100:"Depth [m]":WeS
+psbasemap -O -K $R2 $X -Bf10a1000-950:"                                     # of Samples":N
+psbasemap -O -K $R2 $X -Ba1000-900N
+psbasemap -O $R2 $X -Ba1000-850N
+		 
+rm $TMPFILE

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/LWplot_w	Sat Oct 15 23:23:57 2011 -0400
@@ -0,0 +1,97 @@
+#!/bin/sh
+#======================================================================
+#					 L W P L O T _ W 
+#                    doc: Sat Oct 15 13:42:50 2011
+#                    dlm: Sat Oct 15 22:24:48 2011
+#                    (c) 2011 A.M. Thurnherr
+#                    uE-Info: 56 31 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Oct 15, 2011: - created
+
+#--------------------------------------------------
+# Usage
+#--------------------------------------------------
+
+USAGE="Usage: $0 [-s)ave data to <file>] [in-file]"
+
+args=`getopt s: $*`
+[ $? != 0 ] && {
+	echo $USAGE >&2
+	exit 1
+}
+
+set -- $args
+for opt
+do
+	case "$opt" in
+		-s) opt_s="$2";
+			rm -f "$opt_s";
+			shift; shift;;
+		--) shift; break;;
+	esac
+done
+
+[ $# -gt 1 ] && {
+	echo $USAGE >&2
+	exit 1
+}
+
+if [ $# -eq 1 ]
+then
+	if [ -n "$opt_s" ]
+	then
+		exec <$1 $0 -s $opt_s
+	else
+		exec <$1 $0
+	fi
+fi
+
+#--------------------------------------------------
+# Read Header Data
+#--------------------------------------------------
+
+[ -f .gmtdefaults4 ] ||
+	gmtset	PAPER_MEDIA letter+ \
+			LABEL_FONT_SIZE 14 ANNOT_FONT_SIZE_PRIMARY 14 \
+			WANT_EURO_FONT true \
+	        PLOT_DEGREE_FORMAT ddd:mm:ssF
+
+[ -n "$opt_s" ] && rm -f "$opt_s"
+
+while [ -z "$fields" ]
+do
+	read line
+	[ -n "$opt_s" ] && echo "$line" >> "$opt_s"
+	[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
+	[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
+	[ -z "$max_depth" ] && max_depth=`expr -- "$line" : '#ANTS#PARAMS#.*max_depth{\([^}]*\)}`
+	[ -z "$min_ens" ] && min_ens=`expr -- "$line" : '#ANTS#PARAMS#.*min_ens{\([^}]*\)}`
+	[ -z "$max_ens" ] && max_ens=`expr -- "$line" : '#ANTS#PARAMS#.*max_ens{\([^}]*\)}`
+	[ -z "$LADCP_bin_length" ] && LADCP_bin_length=`expr -- "$line" : '#ANTS#PARAMS#.*LADCP_bin_length{\([^}]*\)}`
+	[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+done
+
+#--------------------------------------------------
+# Plot Data
+#--------------------------------------------------
+
+R=-R`echo "scale=1;$min_ens-0.5"|bc`/`echo "scale=1;$max_ens+0.5"|bc`/`\
+	 echo "scale=1;$min_depth-$LADCP_bin_length/2"|bc`/`echo "scale=1;$max_depth+$LADCP_bin_length/2"|bc`
+U=-R0/1/0/1
+X=-JX10/-10
+
+ens_width=`echo "scale=5;10/($max_ens-$min_ens+1)"|bc`
+bin_length=`echo "scale=5;10*$LADCP_bin_length/($max_depth-$min_depth+$LADCP_bin_length)"|bc`
+
+while read ens bin elapsed depth CTD_depth downcast w extra
+do
+	[ -n "$opt_s" ] && echo "$ens $bin $elapsed $depth $CTD_depth $downcast $w $extra" >> "$opt_s"
+	echo $ens $depth $w $ens_width $bin_length
+done | psxy -P -K $R $X -C`which LWplot_w | sed 's@LWplot_w$@w.cpt@'` -Sr
+	
+echo 0.02 0.02 12 0 0 TL $out_basename | pstext -O -K $U $X
+
+psbasemap -O $R $X -Bf50a500:"Ensemble":/f10a100:"Depth [m]":WeSn
+

--- a/bottom_tracking.pl	Tue Oct 11 17:05:36 2011 -0400
+++ b/bottom_tracking.pl	Sat Oct 15 23:23:57 2011 -0400
@@ -1,14 +1,15 @@
 #======================================================================
 #                    B O T T O M _ T R A C K I N G . P L 
 #                    doc: Wed Oct 20 21:05:37 2010
-#                    dlm: Fri Dec 31 14:21:14 2010
+#                    dlm: Tue Oct 11 17:50:08 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 32 30 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 12 51 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
 #	Oct 20, 2010: - created
 #	Dec 30, 2010: - adapted for use with LADCP_w
+#	Oct 11, 2011: - moved defaults to [defaults.pl]
 
 # This code is essentially identical to the one used in LADCPproc. Differences:
 #	1) velocity editing is simpler: no wake editing, no PPI editing, no shear
@@ -16,30 +17,6 @@
 #	2) median/mad calculated instead of mean/stddev
 #	3) u,v not calculated
 
-# Don't look for BT-referenced velocities if package is more than $BT_max_range
-# above seabed. This parameter is frequency dependent and the current value is
-# appropriate (if rather high) for 300kHz Workhorse intruments.
-
-my($BT_max_range) = 300;
-
-
-# The code only tries to bin BT-referenced velocities if a consistent bottom
-# is available in all 4 beams. Ensembles where the range of bin numbers where
-# the maximum echo is found is greater than $max_BIT_bin_range_diff are rejected.
-# In addition to flukes this also rejects ensembles collected with large
-# instrument tilts. The value of 3 is a first guess that has not been explored.
-
-my($BT_max_bin_range_diff) = 3;
-
-
-# If the difference between measured vertical velocity of the seabed (i.e.
-# the package vertical velocity referenced by the seabed) and the vertical
-# velocity of the CTD (from dp/dt) si greater than $BT_max_w_error the current
-# ensemble is ignored and $nBTwFlag is increased. The value of
-# 3cm/s is taken from listBT developed on A0304 cruise.
-
-my($BT_max_w_error) = 0.03;
-
 #----------------------------------------------------------------------
 # bin valid BT-referenced velocities
 #	input:	ensemble number, water depth (with uncertainty)
@@ -125,7 +102,7 @@
 
 	for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
 		next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
-		my($gi) = int($bd[$bin]) / $opt_o;
+		my($gi) = int($bd[$bin]) / $output_bin_size;
 		push(@{$BTw[$gi]},$LADCP{ENSEMBLE}[$ens]->{W}[$bin]-$seafloor_w);
 	}
 }

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/defaults.pl	Sat Oct 15 23:23:57 2011 -0400
@@ -0,0 +1,241 @@
+#======================================================================
+#                    D E F A U L T S . P L 
+#                    doc: Tue Oct 11 17:11:21 2011
+#                    dlm: Sat Oct 15 21:01:56 2011
+#                    (c) 2011 A.M. Thurnherr
+#                    uE-Info: 144 26 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+#	Oct 11, 2011: - created
+#	Oct 12, 2011: - added $processing_param_file
+#	Oct 13, 2011: - added $out_basename, $TL_out, $TL_hist_out
+#	Oct 14, 2011: - added $w_out, $profile_out
+#				  - renamed _out to out_; out_basename to out_basename
+
+# Variable Names:
+#	- variables that are only used in a particular library are
+#	  prefixed with a 2-caps code
+
+#======================================================================
+# Data Input 
+#======================================================================
+
+# File to load cruise/cast-specific processing parameters from
+
+$processing_param_file = './ProcessingParams.pl';
+
+# CTD depth adjustment
+#	- set with -a
+#	- value is added to CTD pressure
+#	- use when CTD has -ve pressures
+
+$CTD_neg_press_offset = &antsFloatOpt($opt_a,0);
+
+
+# suppress 3-beam LADCP solutions
+
+$RDI_Coords::minValidVels = 4 if ($opt_4);							
+
+
+# bins to use in w calculations
+#	- set with -b
+#	- defaults to 2-last
+
+$opt_b = '2,*' unless defined($opt_b);
+($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b);
+croak("$0: cannot decode -b $opt_b\n")
+    unless (numberp($LADCP_firstBin) &&
+            ($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));
+
+#======================================================================
+# Logging and Output
+#======================================================================
+
+#	- there are 4 verbosity levels, selected by -v
+#		0:	only print errors
+#		1:	default, UNIX-like (warnings and info messages that are not produced for every cast)
+#		2:	progress messages and useful information
+#		>2:	debug messges
+#	- the most useful ones of these are 1 & 2
+
+$verbosity_level = &antsCardOpt($opt_v,1);
+
+# output base name
+
+$out_basename = sprintf('%03d%s',$STN,$RUN);
+
+
+# main w output
+
+$out_w = "| LWplot_w -s $out_basename.w > ${out_basename}_w.eps";
+
+
+# w profile output
+
+$out_profile = "| LWplot_prof -s $out_basename.prof > ${out_basename}_prof.eps";
+
+
+# output bin size in meters
+
+$output_bin_size = &antsFloatOpt($opt_o,10);
+
+
+# min w samples required for each vertical-velocity bin
+
+$min_w_nsamp = &antsCardOpt($opt_k,20);
+
+
+# diagnostic plots
+
+$out_TL 	= "| LWplot_TL     > ${out_basename}_TL.eps";
+$out_TLhist = "| LWplot_TLhist > ${out_basename}_TLhist.eps";
+
+
+#======================================================================
+# Data Editing
+#======================================================================
+
+# min correlation
+
+$min_correlation = &antsFloatOpt($opt_c,70);
+
+
+# max tilt (pitch/roll)
+
+$max_tilt = &antsFloatOpt($opt_t,15);
+
+
+# max err vel
+
+$max_allowed_errvel = &antsFloatOpt($opt_e,0.1);
+
+
+# truncate farthest valid velocities
+
+$truncate_farthest_valid_bins = 0;
+
+
+# discard velocities from chosen beam (1-4)
+
+$discard_velocities_from_beam = 0;
+
+
+# max LADCP gap length in seconds
+
+$max_LADCP_reflr_vel_gap = &antsFloatOpt($opt_g,60);
+
+
+# max allowed vertical velocity in m/s
+
+$max_allowed_w = &antsFloatOpt($opt_m,1);
+
+
+# all apparently valid velocities after a gap of at least this length
+# are deleted
+
+$DE_falsepositives_max_gap = 3;
+
+
+# in each ensemble, vertical velocities differing more than $DE_... times
+# mean absolute deviation from median, are considered outliers and
+# removed
+
+$DE_outliers_mad_limit = 2;
+
+
+# ensembles when instrument is shallower than 
+# $surface_layer_depth in meters are removed.
+# possible contamination: ship's hull, thrusters, bubble clouds
+# Inspired by 2011_IWISE station 8
+
+$surface_layer_depth = 25;
+
+#======================================================================
+# Time Lagging
+#======================================================================
+
+# reference layer bins for w for time matching
+
+($refLr_firstBin,$refLr_lastBin) = (2,6);
+
+
+# number of time lags during each of 2 lagging steps
+
+$opt_n = '10,100' unless defined($opt_n);
+@number_of_timelag_windows = split(',',$opt_n);
+croak("$0: cannot decode -n $opt_n\n")
+	unless numberp($number_of_timelag_windows[0]) && numberp($number_of_timelag_windows[1]);
+
+
+# time lag search window widths for each of 2 lagging steps
+#	- full width in seconds
+
+$opt_w = '240,20' unless defined($opt_w);
+@length_of_timelag_windows = split(',',$opt_w);
+croak("$0: cannot decode -w $opt_w\n")
+	unless numberp($length_of_timelag_windows[0]) && numberp($length_of_timelag_windows[1]);
+
+
+# if top 3 lags have spread greater than $TL_max_allowed_three_lag_spread
+# (in CTD scans) they must account for at least $TL_required_timelag_top_three_fraction
+# or there is an error
+#	- $TL_max_allowed_three_lag_spread default was initially set to 2 but found to be 
+#	  violated quite often during 2011_IWISE
+# 	- large spread may indicate dropped CTD scans
+# 	- the optimum value of $TL_max_allowed_three_lag_spread may be cast-duration dependent
+
+$TL_max_allowed_three_lag_spread = 3;
+$TL_required_top_three_fraction = &antsFloatOpt($opt_3,0.6);
+
+
+#======================================================================
+# Seabed Search
+#======================================================================
+
+# # of ensembles around bottom to search
+
+$SS_search_window_halfwidth = 200;	 
+
+
+# max allowed distance of seabed from mode of distribution
+
+$SS_max_allowed_depth_range = 10;
+
+
+# min allowed LADCP distance from seabed for good data
+
+$SS_min_allowed_hab = 20;
+
+
+#======================================================================
+# Bottom Tracking
+#======================================================================
+
+# Don't look for BT-referenced velocities if package is more than $BT_max_range
+# above seabed. This parameter is frequency dependent and the current value is
+# appropriate (if rather high) for 300kHz Workhorse intruments.
+
+$BT_max_range = 300;
+
+
+# The code only tries to bin BT-referenced velocities if a consistent bottom
+# is available in all 4 beams. Ensembles where the range of bin numbers where
+# the maximum echo is found is greater than $max_BIT_bin_range_diff are rejected.
+# In addition to flukes this also rejects ensembles collected with large
+# instrument tilts. The value of 3 is a first guess that has not been explored.
+
+$BT_max_bin_range_diff = 3;
+
+
+# If the difference between measured vertical velocity of the seabed (i.e.
+# the package vertical velocity referenced by the seabed) and the vertical
+# velocity of the CTD (from dp/dt) si greater than $BT_max_w_error the current
+# ensemble is ignored and $nBTwFlag is increased. The value of
+# 3cm/s is taken from listBT developed on A0304 cruise.
+
+$BT_max_w_error = 0.03;
+
+#======================================================================
+
+1;	# return true

--- a/edit_data.pl	Tue Oct 11 17:05:36 2011 -0400
+++ b/edit_data.pl	Sat Oct 15 23:23:57 2011 -0400
@@ -1,9 +1,9 @@
 #======================================================================
 #                    E D I T _ D A T A . P L 
 #                    doc: Sat May 22 21:35:55 2010
-#                    dlm: Tue Oct 11 13:48:21 2011
+#                    dlm: Sat Oct 15 20:59:38 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 173 0 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 163 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -18,6 +18,10 @@
 #	Oct 10, 2011: - added editFalsePositives()
 #				  - BUG: when earth velocities were edited, all were
 #						 counted, not just those between first and lastBin
+#	Oct 11, 2011: - moved defaults to [defaults.pl]
+#	Oct 12, 2011: - added &editSurfLayer()
+#				  - BUG: editSideLobes() was slightly loose
+#	Oct 15, 2011: - added editWOutliers()
 
 # NOTES:
 #	- all bins must be edited (not just the ones between $LADCP_firstBin
@@ -131,6 +135,7 @@
 
 	my($nrm) = 0;
 	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+		next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
 		next if (abs($LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) <= $lim);
 		undef($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
 		$nrm++ if ($bin>=$LADCP_firstBin-1 && $bin<=$LADCP_lastBin-1);
@@ -139,6 +144,31 @@
 }
 
 #======================================================================
+# $edited = editWOutliers($ens)
+#
+# NOTES:
+#	- call after Earth vels have been calculated
+#	- count only edited vels in selected bin range
+#	- $DE_outliers_mad_limit determines what is an outlier
+#======================================================================
+
+sub editWOutliers($)
+{
+	my($ens) = @_;
+	my($medw) = median(@{$LADCP{ENSEMBLE}[$ens]->{W}});
+	my($madw) = mad2($medw,@{$LADCP{ENSEMBLE}[$ens]->{W}});
+	
+	my($nrm) = 0;
+	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+		next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+		next if (abs($LADCP{ENSEMBLE}[$ens]->{W}[$bin]-$medw) <= $DE_outliers_mad_limit*$madw);
+		undef($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+		$nrm++ if ($bin>=$LADCP_firstBin-1 && $bin<=$LADCP_lastBin-1);
+	}
+	return $nrm;
+}
+
+#======================================================================
 # $edited = editTruncRange($ens,$nbins)
 #
 # NOTES:
@@ -169,8 +199,6 @@
 #		   length >= $FP_BAD_GAP; initial gap is not counted as such
 #======================================================================
 
-my($FP_BAD_GAP) = 3;
-
 sub editFalsePositives($)
 {
 	my($ens) = @_;
@@ -187,7 +215,7 @@
 	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
 		if ($s == 9) {													# skip initial gap
 			$s = 0 if defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
-		} elsif ($s == $FP_BAD_GAP) {									# gap too long => delete
+		} elsif ($s == $DE_falsepositives_max_gap) {					# gap too long => delete
 			next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
 			undef($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
 			$nrm++ if ($bin>=$LADCP_firstBin-1 && $bin<=$LADCP_lastBin-1);
@@ -217,7 +245,7 @@
 				   ? $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}
 				   : $wd - $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH};
 		my($sscorr) = $CTD{SVEL}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}] / 1500;
-		my($goodBins) =   ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER})
+		my($goodBins) =   ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER}) * cos(rad($LADCP{BEAM_ANGLE}))
 						/ ($sscorr*$LADCP{BIN_LENGTH})
 						- 1.5;
 
@@ -234,4 +262,26 @@
 	return ($nvrm,$nerm);
 }
 
+
+#======================================================================
+# $nerm = editSurfLayer($fromEns,$toEns,$surface_layer_depth)
+#
+# NOTES:
+#	1) When this code is executed the fully corrected instrument and
+#	   bin depths are known
+#	2) This code was inspired by 2011_IWISE station 8
+#	3) No point in counting the deleted velocities
+#======================================================================
+
+sub editSurfLayer($$$)
+{
+	my($fe,$te,$sld) = @_;		# first & last ens to process
+	my($nerm) = 0;				# of ensembles affected
+	for (my($e)=$fe; $e<=$te; $e++) {
+		undef($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}),$nerm++
+			if ($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH} <= $sld);
+	}
+	return $nerm;
+}
+
 1;

--- a/find_seabed.pl	Tue Oct 11 17:05:36 2011 -0400
+++ b/find_seabed.pl	Sat Oct 15 23:23:57 2011 -0400
@@ -1,14 +1,16 @@
 #======================================================================
 #                    F I N D _ S E A B E D . P L 
 #                    doc: Sun May 23 20:26:11 2010
-#                    dlm: Thu Dec 30 21:45:02 2010
+#                    dlm: Tue Oct 11 18:09:06 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 12 0 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 13 48 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
 #	May 23, 2010: - adapted from [perl-tools/RDI_Utils.pl]
 #	Dec 25, 2010: - adapted to changes in [LADCP_w]
+#	Oct 11, 2011: - moved defaults to [defaults.pl]
+#				  - increased z_offset from 10km to 15km
 
 # NOTES:
 #	1) BT range is corrected for sound speed at the transducer. This is not
@@ -24,10 +26,7 @@
 # (seabed median depth, mad) = find_seabed(dta ptr, btm ensno, coord flag)
 #======================================================================
 
-my($search_width) = 200;	# # of ensembles around bottom to search
-my($mode_width) = 10;		# max range of bottom around mode
-my($min_dist) = 20;			# min dist to seabed for good data
-my($z_offset) = 10000;		# shift z to ensure +ve array indices
+my($z_offset) = 15000;		# shift z to ensure +ve array indices
 
 sub find_seabed($$$)
 {
@@ -35,9 +34,9 @@
 	my($i,$dd,$sd,$nd);
 	my(@guesses);
 
-	return undef unless ($be-$search_width >= 0 &&
-						 $be+$search_width <= $#{$d->{ENSEMBLE}});
-	for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+	return undef unless ($be-$SS_search_window_halfwidth >= 0 &&
+						 $be+$SS_search_window_halfwidth <= $#{$d->{ENSEMBLE}});
+	for ($i=$be-$SS_search_window_halfwidth; $i<=$be+$SS_search_window_halfwidth; $i++) {
 		next unless (defined($d->{ENSEMBLE}[$i]->{CTD_DEPTH}) &&
 					 defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[0]) &&
 					 defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[1]) &&
@@ -56,7 +55,7 @@
 			 $d->{ENSEMBLE}[$i]->{BT_RANGE}[3]/4;
 		$d->{ENSEMBLE}[$i]->{DEPTH_BT} *= cos(rad($d->{BEAM_ANGLE}));
 		$d->{ENSEMBLE}[$i]->{DEPTH_BT} *= $CTD{SVEL}[$d->{ENSEMBLE}[$i]->{CTD_SCAN}]/1500;
-		next unless ($d->{ENSEMBLE}[$i]->{DEPTH_BT} >= $min_dist);
+		next unless ($d->{ENSEMBLE}[$i]->{DEPTH_BT} >= $SS_min_allowed_hab);
 		$d->{ENSEMBLE}[$i]->{DEPTH_BT} *= -1
 			if ($d->{ENSEMBLE}[$i]->{XDUCER_FACING_UP});
 		$d->{ENSEMBLE}[$i]->{DEPTH_BT} += $d->{ENSEMBLE}[$i]->{CTD_DEPTH};
@@ -76,9 +75,9 @@
 	}
 
 	$nd = 0;
-	for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+	for ($i=$be-$SS_search_window_halfwidth; $i<=$be+$SS_search_window_halfwidth; $i++) {
 		next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
-		if (abs($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$mode) <= $mode_width) {
+		if (abs($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$mode) <= $SS_max_allowed_depth_range) {
 			$dd += $d->{ENSEMBLE}[$i]->{DEPTH_BT};
 			$nd++;
 		} else {
@@ -88,7 +87,7 @@
 	return undef unless ($nd >= 2);
 
 	$dd /= $nd;
-	for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+	for ($i=$be-$SS_search_window_halfwidth; $i<=$be+$SS_search_window_halfwidth; $i++) {
 		next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
 		$sd += ($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$dd)**2;
 	}

--- a/time_lag.pl	Tue Oct 11 17:05:36 2011 -0400
+++ b/time_lag.pl	Sat Oct 15 23:23:57 2011 -0400
@@ -1,9 +1,9 @@
 #======================================================================
 #                    T I M E _ L A G . P L 
 #                    doc: Fri Dec 17 21:59:07 2010
-#                    dlm: Fri Sep 23 21:08:01 2011
+#                    dlm: Fri Oct 14 22:11:30 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 101 42 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 240 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -22,6 +22,20 @@
 #				    more than 2/3 of lags (instead of 50%)
 #	Sep 23, 2011: - added mad info to best lag guesses
 #				  - removed window-doubling heuristics
+#	Oct 12, 2011: - moved defaults to [defaults.pl]
+#				  - BUG: code did not work correctly when there were less than 3
+#						 valid offsets
+#				  - BUG: code did not work correctly unless all windows
+#						 returned valid offsets
+#	Oct 13, 2011: - added $TL_out and $TL_hist_out
+#				  - restricted 2nd-pass to time lags +-1.5s
+#				  - tightened allowed spread for weighted-average calc
+#				  - BUG: $le in &bestLag() could be > $lastGoodEns
+#				  - disabled weighted-average offset calculation
+#	Oct 14, 2011: - improved handling of %PARAMs
+#				  - BUG: last ens of window estimation was off, probably accounting
+#						 for Oct 13 BUG (fix disabled)
+#				  - renamed _out to out_
 
 # DIFFICULT STATIONS:
 #	NBP0901#005
@@ -29,12 +43,7 @@
 # TODO:
 #	- better seabed code (from LADCPproc)
 #	- intermediate-step timelagging guess
-#	- flip aliased ensembles
-
-my($MAX_ALLOWED_THREE_LAG_SPREAD) = 3;			# this was initially set to 2 but found to be
-												# violated quite often during 2011_IWISE. A
-												# large spread may indicate dropped CTD scans.
-												# The optimum value may be cast-duration dependent.
+#	? flip aliased ensembles
 
 sub mad_w($$$)									# mean absolute deviation
 {
@@ -50,7 +59,7 @@
 		) unless (abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2);
 		next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W});
 		my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w);
-		next unless (abs($dw) <= $opt_m);
+		next unless (abs($dw) <= $max_allowed_w);
 
 		$LADCP_mean_w += $LADCP{ENSEMBLE}[$e]->{REFLR_W};
 		$CTD_mean_w   += $CTD{W}[$s+$so];
@@ -64,7 +73,7 @@
 		my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
 		my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w);
 		next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W});
-		next unless (abs($dw) <= $opt_m);
+		next unless (abs($dw) <= $max_allowed_w);
 		$sad += abs($dw);
 		$n++;
 	}
@@ -75,6 +84,10 @@
 sub bestLag($$$$)								# find best lag in window
 {
 	my($fe,$le,$ww,$soi) = @_;					# first/last LADCP ens, window width, scan-offset increment
+#	$le = $lastGoodEns
+#		unless ($le <= $lastGoodEns);
+	die("assertion failed\n\tfe = $fe, le = $le, firstGoodEns = $firstGoodEns, lastGoodEns = $lastGoodEns")
+		unless ($fe>=$firstGoodEns && $le<=$lastGoodEns);
 
 	my($bestso) = 0;							# error at first-guess offset
 	my($bestmad) = mad_w($fe,$le,0);
@@ -98,6 +111,7 @@
 	my($n_windows,$w_size,$scan_increment) = @_;
 
 RETRY:
+	my($failed) = undef;
 	progress("Calculating $n_windows time lags from ${w_size}s-long windows at %dHz resolution...\n",
 		int(1/$scan_increment/$CTD{DT}+0.5));
 
@@ -113,31 +127,48 @@
 
 	my($skip_ens) = int(($approx_joint_profile_end_ens - $approx_joint_profile_start_ens) / 10 + 0.5);
 
-	my(%nBest,%madBest);
+	my(@elapsed,@so,@mad,%nBest,%madBest);
+	my($n_valid_windows) = 0;
 	for (my($wi)=0; $wi<$n_windows; $wi++) {
-		my($fe) = $approx_joint_profile_start_ens + $skip_ens + $wi*int(($approx_joint_profile_end_ens-$approx_joint_profile_start_ens-2*$skip_ens)/$n_windows+0.5);
-		my($so,$mad) = bestLag($fe,$fe+int($w_size/$LADCP{MEAN_DT}+0.5),$w_size,$scan_increment);
-		debugmsg("%.1f cm/s mad(w) at %3d scans offset\n",100*$mad,$so);
+		my($fe) = $approx_joint_profile_start_ens + $skip_ens +
+					int(($approx_joint_profile_end_ens-$approx_joint_profile_start_ens-2*$skip_ens)*$wi/$n_windows+0.5);
+		my($so,$mad) = bestLag($fe,$fe+int($w_size/$LADCP{MEAN_DT}+0.5),$scan_increment==1?3:$w_size,$scan_increment);
+		$elapsed[$wi] = $LADCP{ENSEMBLE}[$fe+int($w_size/2/$LADCP{MEAN_DT}+0.5)]->{ELAPSED};
+		die("assertion failed\nfe=$fe, lastGoodEns=$lastGoodEns, w_size=$w_size") unless ($elapsed[$wi]);
+		next unless ($mad < 9e99);
+		$so[$wi] = $so; $mad[$wi] = $mad;
+		$n_valid_windows++;
 		$nBest{$so}++; $madBest{$so} += $mad;
 	}
-	
-	my(@best_lag);
 	foreach my $i (keys(%nBest)) {
 		$madBest{$i} /= $nBest{$i};
-		$best_lag[0] = $i if ($nBest{$i} > $nBest{$best_lag[0]});
+	}
+
+	my(@best_lag);
+	foreach my $lag (keys(%nBest)) {
+		$best_lag[0] = $lag if ($nBest{$lag} > $nBest{$best_lag[0]});
 	}
-	foreach my $i (keys(%nBest)) {
-		next if ($i == $best_lag[0]);
-		$best_lag[1] = $i if ($nBest{$i} > $nBest{$best_lag[1]});
+	foreach my $lag (keys(%nBest)) {
+		next if ($lag == $best_lag[0]);
+		$best_lag[1] = $lag if ($nBest{$lag} > $nBest{$best_lag[1]});
+	}
+	foreach my $lag (keys(%nBest)) {
+		next if ($lag == $best_lag[0] || $lag == $best_lag[1]);
+		$best_lag[2] = $lag if ($nBest{$lag} > $nBest{$best_lag[2]});
 	}
-	foreach my $i (keys(%nBest)) {
-		next if ($i == $best_lag[0] || $i == $best_lag[1]);
-		$best_lag[2] = $i if ($nBest{$i} > $nBest{$best_lag[2]});
+	if ($nBest{$best_lag[2]}) {
+		progress("\t3 most popular offsets: %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad)\n",
+			$best_lag[0],int(($nBest{$best_lag[0]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[0]},
+			$best_lag[1],int(($nBest{$best_lag[1]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[1]},
+	        $best_lag[2],int(($nBest{$best_lag[2]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[2]});
+	} elsif ($nBest{$best_lag[1]}) {
+		progress("\toffsets: %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad)\n",
+			$best_lag[0],int(($nBest{$best_lag[0]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[0]},
+			$best_lag[1],int(($nBest{$best_lag[1]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[1]});
+	} else {
+		progress("\toffset: %d (%d%% %.1fcm/s mad)\n",
+			$best_lag[0],int(($nBest{$best_lag[0]}/$n_valid_windows)*100+0.5),100*$madBest{$best_lag[0]});
 	}
-	progress("\t3 most popular offsets: %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad), %d (%d%% %.1fcm/s mad)\n",
-		$best_lag[0],int(($nBest{$best_lag[0]}/$n_windows)*100+0.5),100*$madBest{$best_lag[0]},
-		$best_lag[1],int(($nBest{$best_lag[1]}/$n_windows)*100+0.5),100*$madBest{$best_lag[1]},
-		$best_lag[2],int(($nBest{$best_lag[2]}/$n_windows)*100+0.5),100*$madBest{$best_lag[2]});
 
 # BETTER HEURISTIC NEEDED!
 ###	if ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} <= 6) {
@@ -147,28 +178,27 @@
 ###		goto RETRY;
 ###	}
 	
-	unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $opt_3*$n_windows) {
-		if (max(@best_lag)-min(@best_lag) > $MAX_ALLOWED_THREE_LAG_SPREAD) {
-			croak(sprintf("$0: cannot determine a valid lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
-				int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5)))
+	unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $TL_required_top_three_fraction*$n_valid_windows) {
+		if (max(@best_lag)-min(@best_lag) > $TL_max_allowed_three_lag_spread) {
+			$failed = sprintf("$0: cannot determine a valid lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
+				int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_valid_windows+0.5));
 		} else {
 			warning(1,"top 3 tags account for only %d%% of total\n",
-				int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5));
+				int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_valid_windows+0.5));
 		}
 	}
 
-	my($bmo);
-	if (max(@best_lag)-min(@best_lag) > 5 || $nBest{$best_lag[0]}/$n_windows >= 2/3) {
-		$bmo = $best_lag[0];
-		progress("\tunambiguously best offset = %d scans\n",$bmo);
-	} else {
-		$bmo = ($nBest{$best_lag[0]}*$best_lag[0] +
-				$nBest{$best_lag[1]}*$best_lag[1] +
-				$nBest{$best_lag[2]}*$best_lag[2]) / ($nBest{$best_lag[0]} +
-													  $nBest{$best_lag[1]} +
-													  $nBest{$best_lag[2]});
-		progress("\tweighted-mean offset = %.1f scans\n",$bmo);
-	}
+	my($bmo) = $best_lag[0];
+#	if (max(@best_lag)-min(@best_lag) > 3 || $nBest{$best_lag[0]}/$n_valid_windows >= 2/3) {
+#		progress("\tunambiguously best offset = %d scans\n",$bmo);
+#	} else {
+#		$bmo = ($nBest{$best_lag[0]}*$best_lag[0] +
+#				$nBest{$best_lag[1]}*$best_lag[1] +
+#				$nBest{$best_lag[2]}*$best_lag[2]) / ($nBest{$best_lag[0]} +
+#													  $nBest{$best_lag[1]} +
+#													  $nBest{$best_lag[2]});
+#		progress("\tweighted-mean offset = %.1f scans\n",$bmo);
+#	}
 
 	if ($bmo > 0.9*$w_size/2/$CTD{DT}) {
 		warning(0,"lag too close to the edge of the window --- trying again after adjusting the guestimated offset\n");
@@ -183,6 +213,44 @@
 		goto RETRY;
 	}
 
+	if (defined($out_TL) && $scan_increment==1) {
+		progress("\tsaving/plotting time-lagging time series...\n");
+	
+		my($saveParams) = $antsCurParams;
+		@antsNewLayout = ('elapsed','scan_offset','mad');
+		open(STDOUT,"$out_TL") || croak("$out_TL: $!\n");
+
+		&antsAddParams('best_scan_offset',$bmo);
+		&antsAddParams('elapsed.min',$elapsed[0]);
+		&antsAddParams('elapsed.max',$elapsed[$#elapsed]);
+
+		for (my($wi)=0; $wi<@elapsed; $wi++) {
+			&antsOut($elapsed[$wi],$so[$wi],$mad[$wi]);
+		}
+
+		&antsOut('EOF'); close(STDOUT);
+		$antsCurParams = $saveParams;
+	}
+	
+	if (defined($out_TLhist) && $scan_increment==1) {
+		progress("\tsaving/plotting time-lagging histogram...\n");
+	
+		my($saveParams) = $antsCurParams;
+		@antsNewLayout = ('scan_offset','nsamp','mad.avg');
+		open(STDOUT,"$out_TLhist") || croak("$out_TLhist: $!\n");
+	
+		&antsAddParams('n_windows',$n_windows);
+		&antsAddParams('best_scan_offset',$bmo);
+	
+		for (my($so)=-24; $so<=24; $so++) {
+			&antsOut($so,$nBest{$so}?$nBest{$so}:0,$madBest{$so});
+		}
+	
+		&antsOut('EOF'); close(STDOUT);
+		$antsCurParams = $saveParams;
+	}
+
+	croak($failed) if defined($failed);
 	return $CTD{TIME_LAG}+$bmo*$CTD{DT};
 }
 

--- a/time_series.pl	Tue Oct 11 17:05:36 2011 -0400
+++ b/time_series.pl	Sat Oct 15 23:23:57 2011 -0400
@@ -1,9 +1,9 @@
 #======================================================================
 #                    T I M E _ S E R I E S . P L 
 #                    doc: Sun May 23 16:40:53 2010
-#                    dlm: Tue Oct 11 14:08:55 2011
+#                    dlm: Wed Oct 12 12:36:46 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 17 69 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 102 23 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -15,6 +15,8 @@
 #	Jul  2, 2011: - tightened gap-detection code
 #	Jul  4, 2011: - added support for $skip_ens
 #	Oct 11, 2011: - BUG: {DEPTH} had not been set at start of profile
+#	Oct 12, 2011: - re-worked ref_lr_w()
+#				  - stopped depth integration across gaps >= 5s
 
 # NOTES:
 #	- resulting DEPTH field based on integrated w without any sound speed correction
@@ -24,27 +26,16 @@
 sub ref_lr_w($$$$)										# calc ref-layer vert vels
 {
 	my($dta,$ens,$rl_b0,$rl_b1) = @_;
-	my($i,@n,@bn,@v,@vel,@bv,@w);
-
-	for ($i=$rl_b0-1; $i<=$rl_b1-1; $i++) {
-		if (defined($dta->{ENSEMBLE}[$ens]->{W}[$i])) {							# valid w
-			$vel[2] += $dta->{ENSEMBLE}[$ens]->{W}[$i]; $n[2]++;
-			$vel[3] += $dta->{ENSEMBLE}[$ens]->{ERRVEL}[$i], $n[3]++ if defined($dta->{ENSEMBLE}[$ens]->{ERRVEL}[$i]);
-			push(@w,$dta->{ENSEMBLE}[$ens]->{W}[$i]); 							# for stderr test
-		}
-	}
+	my(@w);
 
-	my($w) = $n[2] ? $vel[2]/$n[2] : undef;				# w uncertainty
-	my($sumsq) = 0;
-	for ($i=0; $i<=$#w; $i++) {
-		$sumsq += ($w-$w[$i])**2;
+	for (my($bin)=$rl_b0-1; $bin<=$rl_b1-1; $bin++) {
+		push(@w,$dta->{ENSEMBLE}[$ens]->{W}[$bin])
+			if defined($dta->{ENSEMBLE}[$ens]->{W}[$bin]);
 	}
-	my($stderr) = $n[2]>=2 ? sqrt($sumsq)/($n[2]-1) : undef;
-
-	if (defined($w)) {									# valid w
-		$dta->{ENSEMBLE}[$ens]->{REFLR_W} = $w;
-		$dta->{ENSEMBLE}[$ens]->{REFLR_W_ERR} = $stderr;
-	}
+	return unless (@w);
+	$dta->{ENSEMBLE}[$ens]->{REFLR_W} = avg(@w);
+	$dta->{ENSEMBLE}[$ens]->{REFLR_W_STDDEV} = stddev2($dta->{ENSEMBLE}[$ens]->{REFLR_W},@w);
+	$dta->{ENSEMBLE}[$ens]->{REFLR_W_NSAMP} = @w;
 }
 
 #======================================================================
@@ -107,13 +98,14 @@
 			}
 		}
 	
-		$depth += $dta->{ENSEMBLE}[$lastgood]->{REFLR_W} * $dt;			# integrate
+		$depth += $dta->{ENSEMBLE}[$lastgood]->{REFLR_W} * $dt			# integrate
+			if ($dt < 5);
 		$dta->{ENSEMBLE}[$e]->{DEPTH} = $depth;
 	
 		$atbottom = $e, $max_depth = $depth if ($depth > $max_depth); 
 		$lastgood = $e;
 	}
-	
+
 	return ($firstgood,$lastgood,$atbottom,$w_gap_time);
 }
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/w.cpt	Sat Oct 15 23:23:57 2011 -0400
@@ -0,0 +1,16 @@
+#ANTS# [10/15/11 ...vertical_velocity] mkCPT '-oc' 'polar' '--' '-0.07' '-0.05' '-0.04' '-0.03' '-0.02' '-0.01' '0.01' '0.02' '0.03' '0.04' '0.05' '0.07'
+#ANTS#FIELDS# {from_z} {from_R} {from_G} {from_B} {to_z} {to_R} {to_G} {to_B} 
+-0.07	64	64	255	-0.05	64	64	255
+-0.05	112	112	255	-0.04	112	112	255
+-0.04	143	143	255	-0.03	143	143	255
+-0.03	175	175	255	-0.02	175	175	255
+-0.02	207	207	255	-0.01	207	207	255
+-0.01	255	255	255	0.01	255	255	255
+0.01	255	207	207	0.02	255	207	207
+0.02	255	175	175	0.03	255	175	175
+0.03	255	143	143	0.04	255	143	143
+0.04	255	112	112	0.05	255	112	112
+0.05	255	64	64	0.07	255	64	64
+B	0	0	255
+F	255	0	0
+N	127	127	127