#                    dlm: Thu Oct 25 21:24:59 2012

#                    uE-Info: 22 41 NIL 0 0 70 2 2 4 NIL ofnI

#	Aug 22, 2012: - added documentation

#				  - added T_w()

#	Sep 24, 2012: - made "k" argument default in T_w()

#	Oct 25, 2012: - renamed T_SM to T_ASM

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

# Spectral corrections for LADCP data

#

# NOTES:

#	- see Polzin et al. (JAOT 2002), Thurnherr (JAOT 2012)

#	- to correct, multiply power densities (or power, I think) with corrections

#	- apply to down-/up-cast data only

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

# 1. Corrections for individual data acquisition and processing steps

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

# T_ravg(k,blen)

#	- correct for range averaging due to finite pulse and finite receive window

# 	- this version assumes bin-length == pulse-length

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

    my($k,$blen) =

    return 1 / sinc($k*$blen/2/$PI)**4;

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

# T_fdiff(k,dz)

#	- correct for first differencing on a grid with dz spacing

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

    my($k,$dz) =

    return 1 / sinc($k*$dz/2/$PI)**2;

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

# T_interp(k,blen,dz)

#	- correct for CODAS gridding-with-interpolation algorithm

#	- ONLY USED IN UH SOFTWARE

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

    my($k,$blen,$dz) =

    return 1 / sinc($k*$blen/2/$PI)**4 / sinc($k*$dz/2/$PI)**2;

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

# T_binavg(k,dz)

#	- correct for simple bin averaging

#	- Polzin et al. suggest to use blen instead of dz; this must be a typo

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

    my($k,$dz) =

    return 1 / sinc($k*$dz/2/$PI)**2;

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

# T_tilt(k,d')

#	- d' is a length scale that depends on tilt stats and range max

#	- on d' == 0, T_tilt() == 1, i.e. the correction is disabled

#	- d' = dprime(range_max)

#			- is from a quadratic fit to three data points given by Polzin et al.

#			- see Thurnherr (J. Tech. 2012) for notes

#			- on range_max == 0, d' == 0, i.e. the correction is disabled

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

    my($k,$dprime) =

    return $dprime ? 1 / sinc($k*$dprime/2/$PI)**2 : 1;

}

sub dprime($)

{

	return $_[0] ? -1.2 + 0.0857 * $_[0] - 0.000136 * $_[0]**2 : 0;

} 

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

# 2. Processing-Specific Corrections

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

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

# T_UH(k,blen,dz,range_max)

#	- UH implementation of the shear method (WOCE standard)

#	- range_max == 0 disables tilt correction

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

sub T_UH($$$$)

{

	my($k,$blen,$dz,$range_max) =

        &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);

    return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_interp($k,$blen,$dz) * T_tilt($k,dprime($range_max));

# T_ASM(k,blen,dz,range_max)

#	- re-implemented shear method with simple depth binning

#	- range_max == 0 disables tilt correction

sub T_ASM($$$$)

	my($k,$blen,$dz,$range_max) =

        &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);

    return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));

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

# T_VI(k,blen,preavg_dz,grid_dz,range_max)

# T_VI_alt(k,blen,preavg_dz,grid_dz,dprime)

#	- velocity inversion method of Visbeck (J. Tech., 2002)

#	- only valid if pre-averaging into superensembles is used

#	- range_max == 0 disables tilt correction

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

	my($k,$blen,$sel,$dz,$range_max) =

        &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <range max[m]>',@_);

	return T_VI_alt($k,$blen,$sel,$dz,dprime($range_max));        

	my($k,$blen,$sel,$dz,$dprime) =

    return ($sel>0) ? T_ravg($k,$blen) * T_binavg($k,$sel) * T_binavg($k,$dz) * T_tilt($k,$dprime)

    				: T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,$dprime);

}

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

# T_w(k,blen,dz,range_max)

#	- vertical-velocity method of Thurnherr (IEEE 2011)

#	- range_max == 0 disables tilt correction

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

{ my(@fc);

	sub T_w(@)

	{

		my($k,$blen,$dz,$range_max) =

			&antsFunUsage(4,'ffff',

				'[vertical wavenumber[rad/s]] <ADCP bin size[m]> <output grid resolution[m]> <range max[m]>',

				\@fc,'k',undef,undef,undef,@_);

		return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));

	}