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

#                    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

#                    (c) 2010 A.M. Thurnherr

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

# NOTES:

#	1) BT range is corrected for sound speed at the transducer. This is not

#	   accurate, but unlikely to be very wrong, at least for deep casts, 

#	   because the vertical sound speed variability near the seabed tends

#	   to be small. The seabed depth is only used for sidelobe editing,

#	   which is done with a generous safety margin (from the UH shear

#	   method implementation).

#	2) Acquisition sound speed of 1500m/s assumed.

#	3) To be reasonably accurate, DEPTH must be from the CTD at this stage.

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

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

sub find_seabed($$$)

{

	my($d,$be,$beamCoords) = @_;

	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++) {

		next unless (defined($d->{ENSEMBLE}[$i]->{CTD_DEPTH}) &&

					 defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[0]) &&

					 defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[1]) &&

					 defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[2]) &&

					 defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[3]));

		my(@BT) = $beamCoords

				? velInstrumentToEarth($d,$i,

					velBeamToInstrument($d,

						@{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}}))

				: velApplyHdgBias($d,$i,@{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}});

		next unless (abs($BT[3]) < 0.05);

		$d->{ENSEMBLE}[$i]->{DEPTH_BT} =

			 $d->{ENSEMBLE}[$i]->{BT_RANGE}[0]/4 +

			 $d->{ENSEMBLE}[$i]->{BT_RANGE}[1]/4 +

 			 $d->{ENSEMBLE}[$i]->{BT_RANGE}[2]/4 +

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

		$d->{ENSEMBLE}[$i]->{DEPTH_BT} *= -1

			if ($d->{ENSEMBLE}[$i]->{XDUCER_FACING_UP});

		$d->{ENSEMBLE}[$i]->{DEPTH_BT} += $d->{ENSEMBLE}[$i]->{CTD_DEPTH};

		if ($d->{ENSEMBLE}[$i]->{DEPTH_BT} > $d->{ENSEMBLE}[$be]->{CTD_DEPTH}) {

			$guesses[int($d->{ENSEMBLE}[$i]->{DEPTH_BT})+$z_offset]++;

			$nd++;

		} else {

			undef($d->{ENSEMBLE}[$i]->{DEPTH_BT});

		}

	}

	return undef unless ($nd>5);

	my($mode,$nmax);

	for ($i=0; $i<=$#guesses; $i++) {			# find mode

		$nmax=$guesses[$i],$mode=$i-$z_offset

			if ($guesses[$i] > $nmax);

	}

	$nd = 0;

	for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {

		next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});

		if (abs($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$mode) <= $mode_width) {

			$dd += $d->{ENSEMBLE}[$i]->{DEPTH_BT};

			$nd++;

		} else {

			undef($d->{ENSEMBLE}[$i]->{DEPTH_BT});

		}

	}

	return undef unless ($nd >= 2);

	$dd /= $nd;

	for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {

		next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});

		$sd += ($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$dd)**2;

	}

	return ($dd, sqrt($sd/($nd-1)));

}