LADCPproc.BT
changeset 3 711dd29cb6dd
parent 1 54222c82435f
child 12 65582c172355
--- a/LADCPproc.BT
+++ b/LADCPproc.BT
@@ -1,52 +1,65 @@
 #======================================================================
 #                    L A D C P P R O C . B T 
 #                    doc: Wed Oct 20 21:05:37 2010
-#                    dlm: Mon Jan 10 12:07:17 2011
+#                    dlm: Fri Jul  8 03:24:56 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 11 28 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 148 69 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
 #	Oct 20, 2010: - created
 #	Jan 10, 2010: - -o => -k
+#	Jul  7, 2010: - added $DEBUG
+#				  - added BTrangeFlag
+#				  - added $BT processing parameters
+#				  - changed from echo amplitude to Sv
 
 my($BEAM1) = 0;
 my($BEAM2) = 1;
 my($BEAM3) = 2;
 my($BEAM4) = 3;
 
-my($nBTfound,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0);
+my($nBTfound,$nBTrangeFlag,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0,0);
+
+my($DEBUG) = 0;
 
 sub binBTprof($)
 {
 	my($ens) = @_;
 
-	my(@ea_max) = (0,0,0,0); my(@ea_max_bin) = (nan,nan,nan,nan);
-	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
-		$ea_max[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM1],
-		$ea_max_bin[$BEAM1] = $bin
-			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM1] > $ea_max[$BEAM1]);
-		$ea_max[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM2],
-		$ea_max_bin[$BEAM2] = $bin
-			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM2] > $ea_max[$BEAM2]);
-		$ea_max[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM3],
-		$ea_max_bin[$BEAM3] = $bin
-			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM3] > $ea_max[$BEAM3]);
-		$ea_max[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM4],
-		$ea_max_bin[$BEAM4] = $bin
-			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM4] > $ea_max[$BEAM4]);
+	my(@Sv_max) = (-9e99,-9e99,-9e99,-9e99); my(@Sv_max_bin) = (nan,nan,nan,nan);
+	for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
+		if (defined($BT_min_depth)) {								# manually supplied bottom depth range
+			my($dob) = &depthOfBin($ens,$bin);
+			next unless ($dob >= $BT_min_depth && $dob <= $BT_max_depth);
+		}
+		$Sv_max[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1],
+		$Sv_max_bin[$BEAM1] = $bin
+			if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1] > $Sv_max[$BEAM1]);
+		$Sv_max[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2],
+		$Sv_max_bin[$BEAM2] = $bin
+			if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2] > $Sv_max[$BEAM2]);
+		$Sv_max[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3],
+		$Sv_max_bin[$BEAM3] = $bin
+			if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3] > $Sv_max[$BEAM3]);
+		$Sv_max[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4],
+		$Sv_max_bin[$BEAM4] = $bin
+			if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4] > $Sv_max[$BEAM4]);
 	}
 
-#	print(STDERR "@ea_max | @ea_max_bin\n");
-
-	return unless (max(@ea_max_bin)-min(@ea_max_bin) <= 3);			# inconsistent range (&, impliclity, large tilt)
-		# SHOULD CONSIDER RELAXING LATER
+	print(STDERR "@Sv_max | @Sv_max_bin\n") if ($DEBUG);
 	$nBTfound++;
 
-	my($range_bin) = round(avg(@ea_max_bin));
-#	printf(STDERR "water_depth = $water_depth; BT peak depth = %d\n",depthOfBin($ens,$range_bin));
+	$nBTrangeFlag++,return											# inconsistent range (&, impliclity, large tilt)
+		unless (max(@Sv_max_bin)-min(@Sv_max_bin) <= $BT_max_bin_spread);
+
+	my($range_bin) = round(avg(@Sv_max_bin));
+	printf(STDERR "water_depth = $water_depth; BT peak depth = %d in bin $range_bin\n",depthOfBin($ens,$range_bin))
+		if ($DEBUG);
+
 	$nBTdepthFlag++,return											# BT range inconsistent with water depth
-		unless (abs($water_depth-depthOfBin($ens,$range_bin)) < 20+$sig_water_depth);
+		unless defined($BT_min_depth) ||
+			   (abs($water_depth-depthOfBin($ens,$range_bin)) < $sig_water_depth + $BT_max_depth_error);
 
 	# try bin of max plus one above and below
 	# this does not really work because, often, only one of the bins has valid velocities
@@ -54,7 +67,8 @@
 	my($w2) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
 	my($w3) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
 
-#	printf(STDERR "w123 = %.1f,%.1f,%.1f\n",$w1,$w2,$w3);
+	printf(STDERR "w123 = %.1f,%.1f,%.1f\n",$w1,$w2,$w3)
+		if ($DEBUG);
 
 	$w1 = 9e99 unless numberp($w1);									# no valid velocities
 	$w2 = 9e99 unless numberp($w1);
@@ -78,12 +92,13 @@
 			$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
 	}
 
-	$nBTwFlag++,return if (abs($CTD_w-$LADCP{ENSEMBLE}[$ens]->{W}) > 0.03);		# velocity error is too great
+	$nBTwFlag++,return if (abs($CTD_w-$LADCP{ENSEMBLE}[$ens]->{W}) > $BT_max_w_difference);
 
-#	printf(STDERR "good BT [%5.2f %5.2f %5.2f] found at ens $ens\n",$CTD_u,$CTD_v,$CTD_w);
+	printf(STDERR "good BT [%5.2f %5.2f %5.2f] found at ens $ens\n",$CTD_u,$CTD_v,$CTD_w)
+		if ($DEBUG);
 
 	if ($opt_k) {
-		for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+		for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
 			next if ($edit_flags[$ens][$bin]);
 			printf(BTF "%d %d %d %f %f %f %f %f %f %f %f %f %f %f\n",
 				$LADCP{ENSEMBLE}[$ens]->{NUMBER},
@@ -100,7 +115,7 @@
 	    print(BTF "nan nan nan nan\n");
 	}
 
-	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+	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;
 		push(@{$BTu_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$CTD_u);
@@ -120,15 +135,19 @@
 	}
 
 	for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
-		next unless ($water_depth-$LADCP{ENSEMBLE}[$ens]->{DEPTH} < 300);
+		next unless ($water_depth-$LADCP{ENSEMBLE}[$ens]->{DEPTH} < $BT_begin_search_above);
 		binBTprof($ens);
 	}
 
 	if ($opt_d) {
 		print(STDERR "\n\t$nBTfound BT ensembles found\n");
+	    print(STDERR "\t\t$nBTrangeFlag flagged bad because of inconsistent range to seabed\n");
 	    print(STDERR "\t\t$nBTdepthFlag flagged bad because of wrong bottom depth\n");
-	    print(STDERR "\t\t$nBTvalidVelFlag flagged bad because of no valid velocities\n");
-	    print(STDERR "\t\t$nBTwFlag flagged bad because of incorrect vertical velocity");
+	    print(STDERR "\t\t$nBTvalidVelFlag flagged bad because of lack of valid velocities\n");
+	    print(STDERR "\t\t$nBTwFlag flagged bad because of incorrect vertical velocities");
+	    printf(STDERR "\n\t=> %d velocities from %d BT ensembles used",
+	    				scalar(@BTu_vals),
+						$nBTfound-$nBTrangeFlag-$nBTdepthFlag-$nBTvalidVelFlag-$nBTwFlag);
 	}
 
 	@BTu  = @BTv  = @BTw  = ();