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

#                    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

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

# 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

#	Jul  7, 2011: - use $BEAM? vars to clarify code

#				  - save SV values to use in BT code

#	Jul 14, 2011: - implemented Sv_V04

#	Jul 15, 2011: - modified Sv_V04 to take noise levels into account (Sv_T11)

#	Apr 10, 2012: - BUG: ensonified volume in Sv_T11 was wrong

#	May 16, 2012: - BUG: find_backscatter_seabed() used wrong start depth for

#						 search

#				  - BUG: same used bins entirely below seabed (only possible

#						 for shallow casts)

#	Mar 21, 2014: - adapted to new [LADCPproc.utils]

#	Mar 27, 2014: - adapted to depthOfBinAlongBeam()

#	May 25, 2014: - made search_below integer

#	Jul 27, 2014: - moved depthOfGI() to [LADCPproc.utils]

#	Aug  5, 2014: - BUG: find_backscatter_seabed() discarded everything if

#						 LADCP bin 1 had the backscatter max at the edge of

#						 the search domain, which could easily happen when

#						 the bottom stop was a long way from the seabed

my($BEAM1) = 0;

my($BEAM2) = 1;

my($BEAM3) = 2;

my($BEAM4) = 3;

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

# Volume Scattering Coefficient

#

# 1) There are currently 4 different methods for correcting echo

#	 amplitudes as a function of range. The method to calculate Sv

#    of Deines (1999) is the only published one, although Visbeck's

#    correction was presumably used by Visbeck and Thurnherr (2009).

# 	 The third method (Thurnherr 2011) improves on Visbeck's code.

#	 The fourth method consists in no correction.

# 2) The following observations stand out when plotting echo amplitudes

#	 against bin number in clear water (e.g. below 1400m in first

#	 downcast of IWISE yoyo cast 160; [Documentation/Sv_comparison.eps]):

#		- the corrections works approximately for the valid bins

#		  (2-13 in case of IWISE) but it amounts to a slight

#		  overcorrection, i.e. the resulting corrected backscatter

#		  profiles increase with increasing distance from the

#		  transducer

#		- in terms of remaining vertical structure, the different

#		  methods are nearly identical

#		- however, Visbeck (2004) shows fairly large instrument

#		  dependent biases because (instrument-specific) noise

#		  is not considered

#		- the vertical structure can be improved by decreasing the

#		  attenuation parameter to something like 0.05; I have

#		  verified this for the Deines and Thurnherr expressions

#		  with a single IWISE station; however, there remain

#		  small (but noticeable) instrument-specific biases

#	3) Therefore, the default code is Deines (1999); if this is

#	   not good enough, it is probably possible to derive an

#	   empirical correction from in-situ data.

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

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

# 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

# Results:

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

{

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

    my($Kc)     = 0.45;

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

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

}

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

# Volume Scattering Coefficient, following Visbeck (code 2004)

#

## function [ts,bcs]=targ($EA,$R,$alpha,$PL,$EAS,$ap)

## Target strength of EA for volume scatterer

## $EA = echoamp in  dB

## $R = distance in  m

## $alpha = attenuation dB/m

## $PL = pulse/bin legth in  m

## $EAS = source level

## $ap = aperature in degree

## M. Visbeck 2004

#

# Results:

#	- overall, correction with distance works similarly well to Deines (1999)

#	- constant bias, which could be taken care of by changing EAS

#	- however, much more serious UL/DL differences than Deines (1999)

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

sub Sv_V04($$$$$)

{

    my($temp,$PL,$Er,$R,$EA) = @_;	# only uses $R and $EA

#	my($alpha) = 0.039;		## attenuation dB/m for 150 kHz

	my($alpha) = 0.062;		## attenuation dB/m for 300 kHz

	my($EAS) = 100;			## source level in dB

	my($ap) = rad(2);		## beam aperature in DEGREE convert to radian

	my($r1) = tan($ap)*($R-$PL/2);	## radius of top and bottom of each bin

	my($r2) = tan($ap)*($R+$PL/2);

	my($V) = $PI*$PL/3 * ($r1**2+$r2**2+$r1*$r2);	## ensonified volume 

	my($TL) = 20*log10($R) + $alpha*$R;				## transmission loss

	my($TS) = 0.45*$EA - $EAS + 2*$TL - 10*log10($V);	## target strength

	my($BCS) = exp($TS/10);							## backscatter cross section

	return $TS;

}

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

# Volume Scattering Coefficient, Thurnherr [2011]

#

# in an attempt to improved on Visbeck's code, I made several

# modifications:

#	- add 5cm to beam radius (transducer radius)

#	- take beam noise into account (as in Deines [1999])

#	- use different equation for ensonified volume

# results:

#	- inclusion of beam noise greatly reduces UL/DL differences in IWISE data

#	- constant offset from Deines values, which could be taken care of by changing EAS

#	- remaining changes don't have great effects

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

sub Sv_T11($$$$$)

{

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

#	my($alpha) = 0.067;				# Deines [1999]

	my($alpha) = 0.062;				# Visbeck [2004]

#	my($alpha) = 0.048;				# attenuation dB/m for 300 kHz

	my($EAS) = 100;					# source level in dB

	my($ap) = rad(2);				# beam aperature in DEGREE convert to radian

	my($r0) = 0.05;					# beam radius at source (transducer radius)

	my($d1) = $R - $PL/2;			# distance to top and bottom of bin

	my($d2) = $R + $PL/2;

	my($r1) = $r0 + $d1*tan($ap);	# radius of top and bottom of each bin

	my($r2) = $r0 + $d2*tan($ap);

	my($V) = $PI/3 * ($d2*$r2**2 - $d1*$r1**2);	# ensonified volume

	my($TL) = 20*log10($R) + $alpha*$R;		# transmission loss

	my($TS) = 0.45*($EA-$Er) - $EAS + 2*$TL - 10*log10($V);	# target strength

	return $TS;

}

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

# no correction for attenuation losses

#

# however, noise level is still taken into consideration

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

sub Sv_nocorr($$$$$)

{

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

	return 0.45*($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[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM1]

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

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

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

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

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

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

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

    }

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

	my($Svfunc);											# Sv method

	if ($opt_u =~ /^[dD]/) {

		$Svfunc = \&Sv_D99;

	} elsif ($opt_u =~ /^[vV]/) {

		$Svfunc = \&Sv_V04;

	} elsif ($opt_u =~ /^[tT]/) {

		$Svfunc = \&Sv_T11;

	} else {

		$Svfunc = \&Sv_nocorr;

	}

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

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

			my($range_to_bin) = &rangeToBin($ens,$bin);

			for (my($beam)=0; $beam<4; $beam++) {

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

				next if ($gi < 0);

				$LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$beam] = &$Svfunc($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},

							 							    		 $LADCP{TRANSMITTED_PULSE_LENGTH},

														    		 $Er[$beam],$range_to_bin,

							            			        		 $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$beam]);

				$sSv[$gi][$bin] += $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$beam];

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

				if ($bin>=$Svbin_start && $bin<=$Svbin_end) {

					$sSv_prof[$gi] += $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$beam];

					$nSv_prof[$gi]++;

				} # if $bin

			} # for $beam

		} # for $bin

	} # for $end

} # sub

sub find_backscatter_seabed($)

{

	my($water_depth) = @_;

	my(@wdepth_gi);												# water_depth indices

	my($search_below) = int(max(0,$water_depth-$BT_begin_search_above));

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

	printf(STDERR "\n\t\tlooking for seabed below %d m (gi = [%d..%d])",$search_below,$mdgi,scalar(@nSv))

		if ($opt_d);

	print(STDERR "\n\t\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,$firstvalid,$lastvalid) = (1e99,-1e99,-1,-1,-1);

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

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

			$firstvalid = $gi if ($firstvalid < 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!=$firstvalid && $gimax!=$lastvalid) { 				# ignore boundary maxima & scatter

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

			push(@wdepth_gi,$gimax);

		}

	}

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

}

1;