#======================================================================
#                    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
#                    dlm: Wed Jun 15 15:31:43 2011
#                    (c) 2010 A.M. Thurnherr
#                    uE-Info: 67 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# HISTORY:
#	Oct 20, 2010: - created
#	Dec 10, 2010: - BUG: backscatter above sea surface made code bomb
#						 when run with uplooker data
#	Jun 15, 2011: - added calculation of backscatter profiles from
#				    subset of bins

#----------------------------------------------------------------------
# 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);
			next if ($gi < 0);
			my($range_to_bin) = &dzToBin($ens,$bin) / cos(rad($LADCP{BEAM_ANGLE}));
			my($Sv) = 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;

			$sSv[$gi][$bin] += $Sv;
			$nSv[$gi][$bin]++;

			if ($bin>=$Svbin_start && $bin<=$Svbin_end) {
				$sSv_prof[$gi] += $Sv;
				$nSv_prof[$gi]++;
			}
		}
	}
}

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;