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

#                    L A D C P P R O C . B A C K S C A T T E R 

#                    doc: Wed Oct 20 13:02:27 2010

#                    (c) 2010 A.M. Thurnherr

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

#----------------------------------------------------------------------

# Volume Scattering Coefficient, following Deines (IEEE 1999)

# NOTES:

#	- instrument specific! (300kHz Workhorse)

#   - no sound-speed correction applied

#   - R in bin center, instead of last quarter

#   - transmit power assumes 33V batteries

#----------------------------------------------------------------------

# NB:

#	- correction seems to work for a subset of bins (~bins 3-9 for 

#	  2010 P403 station 46) 

#	- this may imply that noise level depends on bin

# 	- far bins are important for seabed detection, i.e. cannot simply

#	  be discarded at this stage

sub log10 {

    my $n = shift;

    return log($n)/log(10);

}   

sub Sv($$$$$)

{

    my($temp,$PL,$Er,$R,$EA) = @_;

    my($C)      = -143;                 # RDI Workhorse monitor

    my($Ldbm)   = 10 * log10($PL);

    my($PdbW)   = 14.0;

    my($alpha)  = 0.069;

    my($Kc)     = 0.45;

    return $C + 10*log10(($temp+273)*$R**2) - $Ldbm - $PdbW

              + 2*$alpha*$R + $Kc*($EA-$Er);

}

sub mk_backscatter_profs($$)

{

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

	my(@Er) = (1e99,1e99,1e99,1e99);						# echo intensity reference level

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

		$Er[0] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0]

			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0] < $Er[0]);

		$Er[1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1]

			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1] < $Er[1]);

		$Er[2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2]

			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2] < $Er[2]);

		$Er[3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3]

			if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3] < $Er[3]);

    }

	print(STDERR "\n\t\@per-beam noise levels = @Er") if ($opt_d);

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

		for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

			my($gi) = int(&depthOfBin($ens,$bin) / $GRID_DZ);

			my($range_to_bin) = &dzToBin($ens,$bin) / cos(rad($LADCP{BEAM_ANGLE}));

			$sSv[$gi][$bin] += Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},

								  $LADCP{TRANSMITTED_PULSE_LENGTH},

								  $Er[0],$range_to_bin,

								  $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +

							   Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},

								  $LADCP{TRANSMITTED_PULSE_LENGTH},

								  $Er[1],$range_to_bin,

								  $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +

							   Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},

								  $LADCP{TRANSMITTED_PULSE_LENGTH},

								  $Er[2],$range_to_bin,

								  $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +

							   Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},

								  $LADCP{TRANSMITTED_PULSE_LENGTH},

								  $Er[3],$range_to_bin,

								  $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;

			$nSv[$gi][$bin]++;

		}

	}

}

sub depthOfGI($) { return $_[0]*$GRID_DZ + $GRID_DZ/2; }		# depth corresponding to particular grid index

sub find_backscatter_seabed($)

{

	my($search_below) = @_;

	my($mdgi) = int($search_below/$GRID_DZ);					# grid index to begin search

	my(@wdepth_gi);												# water_depth indices

	print(STDERR "\n\tseabed-max grid indices:") if ($opt_d);

	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { 			# find backscatter min/max below $search_below in each bin

		my($min,$max,$gimax,$lastvalid) = (1e99,-1e99,-1,-1);

		for (my($gi)=$mdgi; $gi<@nSv; $gi++) {

			next unless ($nSv[$gi][$bin] > 0);

			my($avg) = $sSv[$gi][$bin] / $nSv[$gi][$bin];

			$lastvalid = $gi;

			$min = $avg if ($avg < $min);

			$max = $avg, $gimax = $gi if ($avg > $max);

		}

		if ($max-$min>10 && $gimax!=$lastvalid) { 				# ignore boundary maxima & scatter

			printf(STDERR " %d",$gimax-$mdgi) if ($opt_d);

			push(@wdepth_gi,$gimax);

		}

	}

	return (depthOfGI(avg(@wdepth_gi)),stddev(@wdepth_gi)*$GRID_DZ);

}

1;