#======================================================================

#                    L A D C P P R O C . B T 

#                    doc: Wed Oct 20 21:05:37 2010

#                    dlm: Tue Aug  5 14:38:08 2014

#                    (c) 2010 A.M. Thurnherr

#                    uE-Info: 205 0 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

#	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;

my($BEAM3) = 2;

my($BEAM4) = 3;

my($nBTfound,$nBTrangeFlag,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0,0);

my($DEBUG) = 0;

sub CTDvelFromSv($)

{

	my($ens) = @_;

	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 "@Sv_max | @Sv_max_bin\n") if ($DEBUG);

	$nBTfound++;

	$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 defined($BT_min_depth) ||

			   (abs($water_depth-depthOfBin($ens,$range_bin)) < $sig_water_depth + $BT_max_depth_error);

	my($CTD_u,$CTD_v,$CTD_w);

	if ($BT_range_method == 0) {										# take BT vel from bin with Sv max

		$nBTvalidVelFlag++,return unless numberp($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W]);

		$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$U];

		$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$V];

		$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];

	} elsif ($BT_range_method == 1) {									# take "best-fit" with w_reflr

		# try bin of max plus one above and below

		# this does not really work because, often, only one of the bins has valid velocities

		my($w1) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];

		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)

			if ($DEBUG);

		$w1 = 9e99 unless numberp($w1); 								# no valid velocities

		$w2 = 9e99 unless numberp($w1);

		$w3 = 9e99 unless numberp($w1);

		if (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2) &&

			abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w3)) {

				$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$U];

				$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$V];

				$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];

		} elsif (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2)) {

				$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$U];

				$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$V];

				$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];

		} else {

				$nBTvalidVelFlag++,return if ($w2 == 9e99); 			# none of 3 bins has valid velocity

				$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$U];

				$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$V];

				$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];

	    }

	} elsif ($BT_range_method == 2) {									# Visbeck method (median from 3 bins)

		croak("$0: need \$BT_range_Visbeck_center\n")

			unless defined($BT_range_Visbeck_center);

		$CTD_u = median($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center+1][$U],

						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center][$U],

						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center-1][$U]);

		$nBTvalidVelFlag++,return unless numberp($CTD_u);

		$CTD_v = median($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center+1][$V],

						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center][$V],

						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center-1][$V]);

		$CTD_w = median($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center+1][$W],

						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center][$W],

						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center-1][$W]);

	} else {

		croak("$0: unknown \$BT_range_method == $BT_range_method\n");

	}

	return ($CTD_u,$CTD_v,$CTD_w);

}

sub depthAtRange($$)

{

	my($ens,$range) = @_;

	return $LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} ?

	       $LADCP{ENSEMBLE}[$ens]->{DEPTH} - $range :

	       $LADCP{ENSEMBLE}[$ens]->{DEPTH} + $range;

}

sub CTDvelFromBT($)

{

	my($ens) = @_;

	return (undef,undef,undef)

		unless defined($LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$W]);

	$nBTfound++;

	$nBTrangeFlag++,return											# inconsistent range (&, impliclity, large tilt)

		unless (max(@{$LADCP{ENSEMBLE}[$ens]->{BT_RANGE}})-min(@{$LADCP{ENSEMBLE}[$ens]->{BT_RANGE}})

					<= $BT_max_bin_spread*$LADCP{BIN_LENGTH});

	my($range) = round(avg(@{$LADCP{ENSEMBLE}[$ens]->{BT_RANGE}}));

	$nBTdepthFlag++,return											# BT range inconsistent with water depth

		unless defined($BT_min_depth) ||

			   (abs($water_depth-depthAtRange($ens,$range)) < $sig_water_depth + $BT_max_depth_error);

	return ($LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$U],

			$LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$V],

			$LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$W]);

}

sub binBTprof($)

{

	my($ens) = @_;

	my($CTD_u,$CTD_v,$CTD_w);

	($CTD_u,$CTD_v,$CTD_w) = $opt_r ? CTDvelFromBT($ens) : CTDvelFromSv($ens);

	return unless defined($CTD_w);

	$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)

		if ($DEBUG);

	if ($opt_k) {

		for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $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},

				$LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL},

				$CTD_u,$CTD_v,$CTD_w,

				$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U],

				$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V],

				$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W],

				$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$CTD_u,

				$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V]-$CTD_v,

				$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]-$CTD_w);

		}

	    print(BTF "nan nan nan nan\n");

	}

	for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {

		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($$)

{

	my($LADCP_start,$LADCP_end) = @_;

	if ($opt_k) {

		open(BTF,">BT.profs");

	    print(BTF "#ANTS#FIELDS# {ens} {depth} {CTD_depth} {pitch} {roll} {CTD_u} {CTD_v} {CTD_w} {u} {v} {w} {BT_u} {BT_v} {BT_w}\n");

	}

	for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {

		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 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  = ();

	@BTu_sig = @BTv_sig = @BTw_sig = ();

	@BT_nsamp = ();

	for (my($gi)=0; $gi<@BTu_vals; $gi++) {								# calc grid means & stddev

		my($sum_u,$sum_v,$sum_w);

		$BT_nsamp[$gi] = @{$BTu_vals[$gi]};

		for (my($vi)=0; $vi<$BT_nsamp[$gi]; $vi++) {

			$sum_u += $BTu_vals[$gi][$vi];

			$sum_v += $BTv_vals[$gi][$vi];

			$sum_w += $BTw_vals[$gi][$vi];

		}

		$BTu[$gi] = $BT_nsamp[$gi] ? $sum_u/$BT_nsamp[$gi] : nan;

		$BTv[$gi] = $BT_nsamp[$gi] ? $sum_v/$BT_nsamp[$gi] : nan;

		$BTw[$gi] = $BT_nsamp[$gi] ? $sum_w/$BT_nsamp[$gi] : nan;

	}

	for (my($gi)=0; $gi<@BTu_vals; $gi++) {								# calc & grid stddevs

		my($sumsq_u,$sumsq_v,$sumsq_w);

		for (my($vi)=0; $vi<$BT_nsamp[$gi]; $vi++) {

			$sumsq_u += ($BTu_vals[$gi][$vi] - $BTu[$gi])**2;

			$sumsq_v += ($BTv_vals[$gi][$vi] - $BTv[$gi])**2;

			$sumsq_w += ($BTw_vals[$gi][$vi] - $BTw[$gi])**2;

		}

		$BTu_sig[$gi] = $BT_nsamp[$gi]>1 ? sqrt($sumsq_u/($BT_nsamp[$gi]-1)) : nan;

		$BTv_sig[$gi] = $BT_nsamp[$gi]>1 ? sqrt($sumsq_v/($BT_nsamp[$gi]-1)) : nan;

		$BTw_sig[$gi] = $BT_nsamp[$gi]>1 ? sqrt($sumsq_w/($BT_nsamp[$gi]-1)) : nan;

	}

	close(BTF) if ($opt_k);

}

1;