#======================================================================
#                    A C O U S T I C _ B A C K S C A T T E R . P L 
#                    doc: Wed Oct 20 13:02:27 2010
#                    dlm: Tue May 24 16:34:24 2016
#                    (c) 2010 A.M. Thurnherr
#                    uE-Info: 212 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
#	Dec 30, 2010: - adapted for use with [LADCP_w]
#	Oct 19, 2011: - added support for $SS_{min,max}_allowed_range
#				  - BUG: acoustic-backscatter assumed 0 deg C
#				  - SV now saved in ensemble
#	Oct 21, 2011: - BUG: made code work for uplooker again
#	Mar  4, 2014: - added support for missing PITCH/ROLL (TILT)
#	Apr 17, 2014: - BUG: missing ;
#	Nov  7, 2014: - BUG: calc_binDepths() was called without valid CTD depth
#				  - added $SS_use_BT, $SS_min_signal, $SS_min_samp
#	Apr 20, 2015: - added comments
#				  - removed SS_{min,max}_allowed_range from calc_backscatter_profs()
#				  - added correct_backscatter() & linterp() from laptop
#	Apr 21, 2015: - added debug statements
#	May 14, 2015: - BUG: code did not work for partial-depth casts
#	Jun 18, 2015: - removed assertion marked by ##???, which bombed on P16N1#41 DL
#	Jan 26, 2016: - added %PARAMs
#	Mar 26, 2016: - BUG: nSv was declared local to this scope even though it is used outside
#	May 18, 2016: - improved logging
#   May 24, 2016: - calc_binDepths() -> binDepths()

#----------------------------------------------------------------------
# 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
# EMPIRICAL FINDINGS (after applying the empirical correction)
#  - Sv(WH300) ~ Sv(WH150)/1.4-34
#      - based on DoMORE-1: 004, 005
#      - identical instrument setup (0/8/8m)
#      - 003, with higher scattering is somewhat different
#----------------------------------------------------------------------

# 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 Sv($$$$$)
{
    my($temp,$PL,$Er,$R,$EA) = @_;
    my($C)      = -143;                 # RDI WHM300 (from Deines)
    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);
}

#----------------------------------------------------------------------
# Calculate per-bin backscatter profiles
#	input: 	first and last valid LADCP ensemble
#	output:	sSv[$depth][$bin]	sum of volume scattering coefficients
#			nSv[$depth][$bin]	number of samples in bin
#----------------------------------------------------------------------

# my(@sSv,@nSv);		BAD: nSv is referenced in [LADCP_w_ocean]

sub calc_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]);
    }
	debugmsg("per-beam noise levels = @Er\n");

	my($cosBeamAngle) = cos(rad($LADCP{BEAM_ANGLE}));
	for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
		next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
		my(@bd) = binDepths($ens);
		for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
			my($depth) = int($bd[$bin]);
			next if ($depth<0 || !defined($LADCP{ENSEMBLE}[$ens]->{TILT}));
			my($range_to_bin) = abs($bd[$bin] - $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})
									/ cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))
									/ $cosBeamAngle;
			my($temp) = defined($CTD_temp)
					  ? $CTD{TEMP}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]
					  : $LADCP{ENSEMBLE}[$ens]->{TEMPERATURE};
			$LADCP{ENSEMBLE}[$ens]->{SV}[$bin] =
				median(
					Sv($temp,$LADCP{TRANSMITTED_PULSE_LENGTH},$Er[0],$range_to_bin,
					   $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0]
				    ),
				    Sv($temp,$LADCP{TRANSMITTED_PULSE_LENGTH},$Er[1],$range_to_bin,
					   $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1]
					),
					Sv($temp,$LADCP{TRANSMITTED_PULSE_LENGTH},$Er[2],$range_to_bin,
					   $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2]
					),
					Sv($temp,$LADCP{TRANSMITTED_PULSE_LENGTH},$Er[3],$range_to_bin,
	     			   $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3]
					)
				);
    		$sSv[$depth][$bin] += $LADCP{ENSEMBLE}[$ens]->{SV}[$bin];
			$nSv[$depth][$bin]++;
		}
	}
}

#----------------------------------------------------------------------
# empirically adjust Sv in far bins for consistency with nearest
# valid bin
#	- based on bin-to-bin differences of 100m vertically averaged Sv
#	  profiles
#	- algorithm can leave artifacts in the near bins near bottom
#	  turn-around (DM1#001)
#----------------------------------------------------------------------

sub linterp($$$$$)
{
	my($x,$mix,$max,$miy,$may) = @_;
	return $miy + ($x-$mix)/($max-$mix)*($may-$miy);
}

sub correct_backscatter($$)
{
	my($LADCP_start,$LADCP_end) = @_;
	my($bin) = $LADCP_firstBin-1;
	my(@refSvProf,@refSvSamp,$depth,$i);

	&antsAddParams('Sv_ref_bin',$Sv_ref_bin);

RETRY:
	for ($depth=0; $depth<@nSv; $depth++) {						# create reference profile
		next unless ($nSv[$depth][$Sv_ref_bin-1] > 0);
		$refSvProf[int($depth/100)] += $sSv[$depth][$Sv_ref_bin-1] / $nSv[$depth][$Sv_ref_bin-1];
		$refSvSamp[int($depth/100)]++;
    }
    $Sv_ref_bin++,goto RETRY
    	unless (@refSvSamp);

	for ($i=0; $i<@refSvSamp; $i++) {
		next unless ($refSvSamp[$i] > 0);
		$refSvProf[$i] /= $refSvSamp[$i];
	}
	for ($i=0; $i<5; $i++) {									# extrapolate bottom value by 500m
		push(@refSvProf,$refSvProf[$#refSvProf]);
		push(@refSvSamp,$refSvSamp[$#refSvSamp]);
	}
	for ($i=$#refSvProf-1; $i>=0; $i--) {						# extrapolate upward to surface
		next if ($refSvSamp[$i] > 0);
		$refSvProf[$i] = $refSvProf[$i+1];
		$refSvSamp[$i] = $refSvSamp[$i+1];
	}
	info("\tusing bin %d as reference\n",$Sv_ref_bin);

	my(@dSvProf);												# create profiles for all bins
	for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {	
		my(@dSvSamp);									
		for ($depth=0; $depth<@nSv; $depth++) {					# create Sv-difference profile for current bin
			next unless ($nSv[$depth][$bin] > 0);
			$dSvProf[int($depth/100)][$bin] += $sSv[$depth][$bin] / $nSv[$depth][$bin];
			$dSvSamp[int($depth/100)]++;
		}
#		print(STDERR "dSvProf[bin$bin] = ");
		for ($i=0; $i<@dSvSamp; $i++) {
			next unless ($dSvSamp[$i] > 0);
			die("assertion failed (refSvSamp[$i] = $refSvSamp[$i])") unless ($refSvSamp[$i] > 0);
			$dSvProf[$i][$bin] = $dSvProf[$i][$bin]/$dSvSamp[$i] - $refSvProf[$i];
			$dSvProf[$i][$bin] = $dSvProf[$i-1][$bin]
				if (abs($dSvProf[$i][$bin]) > 10);
#			printf(STDERR "%.1f ",$dSvProf[$i][$bin]);
		}
		$i--;													# trim deepest, often anomalous, value
#		print(STDERR "[delete] ");
		while ($i <= int(scalar(@nSv)/100)+1) {
			$dSvProf[$i][$bin] = $dSvProf[$i-1][$bin];			# extrapolate 100m
#			printf(STDERR "[%.1f] ",$dSvProf[$i][$bin]);
			$i++;
		}
#		print(STDERR "\n");
	}

    for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {		# correct Sv data
		next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
		my(@bd) = binDepths($ens);
		for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
			next unless numberp($LADCP{ENSEMBLE}[$ens]->{SV}[$bin]);
			my($depth) = int($bd[$bin]);
			if (numberp($dSvProf[int($depth/100)][$bin]) && numberp($dSvProf[int($depth/100)+1][$bin])) {
#				print(STDERR "\n$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]? $dSvProf[int($depth/100)][$bin] $dSvProf[int($depth/100)+1][$bin] int($depth/100)+1[$bin]");
				$LADCP{ENSEMBLE}[$ens]->{SV}[$bin] -= # $dSvProf[int($depth/100)][$bin];
					linterp($depth,100*int($depth/100),100*int($depth/100)+100,
							$dSvProf[int($depth/100)][$bin],$dSvProf[int($depth/100)+1][$bin]);
			} else {
				$LADCP{ENSEMBLE}[$ens]->{SV}[$bin] = nan;
			}
		}
	}
}

#----------------------------------------------------------------------
# determine location of seabed from backscatter profiles
#	input:	depth below which seabed can possibly be (e.g. max CTD depth)
#	output:	median/mad of estimated water depth
#----------------------------------------------------------------------

sub find_backscatter_seabed($)
{
	my($search_below) = int($_[0]);										# grid index to begin search
	my(@wdepth,@Sv_rng);												# list of water_depth indices

	&antsAddParams('SS_min_signal',$SS_min_signal,'SS_min_samp',$SS_min_samp,
				   'SS_max_allowed_depth_range',$SS_max_allowed_depth_range);
				   
	for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { 	# find backscatter min/max below $search_below in each bin
		my($minSv,$maxSv,$depthmaxSv,$lastvalid) = (1e99,-1e99,-1,-1);
		for (my($depth)=$search_below; $depth<@nSv; $depth++) {
			next unless ($nSv[$depth][$bin] > 0);
			my($Sv) = $sSv[$depth][$bin] / $nSv[$depth][$bin];
			$lastvalid = $depth;
			$minSv = $Sv if ($Sv < $minSv);
			$maxSv = $Sv, $depthmaxSv = $depth if ($Sv > $maxSv);
		}
		if ($maxSv-$minSv >= $SS_min_signal) { 							# ignore scatter
			push(@Sv_rng,round($maxSv-$minSv));
			push(@wdepth,$depthmaxSv);
		}
	}

	if (@wdepth) {
		info("\t%d bins with seabed signatures found (\@Sv_rng: @Sv_rng)\n",scalar(@wdepth));
    } else {
		info("\tno bins with seabed signatures found\n");
    }
	return (undef,undef) if (scalar(@wdepth) < $SS_min_samp);			# require min number of samples
	
	my($wd) = median(@wdepth);
	return ($wd,mad2($wd,@wdepth));

}

1;