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

#                    L I B L A D C P . P L 

#                    doc: Wed Jun  1 20:38:19 2011

#                    dlm: Mon Apr  1 16:13:46 2019

#                    (c) 2011 A.M. Thurnherr

#                    uE-Info: 27 65 NIL 0 0 70 2 2 4 NIL ofnI

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

# HISTORY:

#	Jun  1, 2011: - created

#	Jul 29, 2011: - improved

#	Aug 10, 2011: - made correct

#	Aug 11, 2011: - added "convenient combinations"

#	Aug 18, 2011: - made buoyancy frequency non-constant in S()

#	Jan  4, 2012: - improved T_VI to allow correcting even without superensembles

#	Jan  5, 2012: - removed S(), which is just pwrdens/N^2 (rather than

#				    pwrdens/N^2/(2pi) as I erroneously thought)

#	Jan 18, 2012: - added T_VI_alt() to allow assessment of tilt correction extrema

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

#	Jun 26. 2013: - added T_w_z()

#				  - added parameter checks to processing-specific corrections

#	May 18, 2015: - added pulse length to T_w() and T_w_z()

#	Apr 25, 2018: - added eps_VKE() parameterization

#	Apr  1, 2019: - removed plen parameter from T_w() and T_w_z()

require "$ANTS/libvec.pl";

require "$ANTS/libfuns.pl";

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

# VKE parameterization for epsilon

#

# NOTES:

#	- see Thurnherr et al. (GRL 2015)

#	- calculate eps from p0

#	- optional second argument allows free choice of parameterization constant

#	- default value is from paper, which is slightly lower than the current value

#	  used in [LADCP_VKE], which applies the parameterization only to spectra

#	  passing a few tests

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

sub eps_VKE(@)

{

    my($p0,$c) =

        &antsFunUsage(-1,'.','<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]]',@_);

	$c = 0.021 unless defined($c);

    return numberp($p0) ? ($p0 / $c) ** 2 : nan;

}

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

# 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[,plen])

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

# 	- when called with 2 arguments, bin-length == pulse-length assumed

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

sub T_ravg(@)

{

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

        &antsFunUsage(-2,'ff','<vertical wavenumber[rad/s]> <bin-length[m]> [pulse-length[m]]',@_);

	$plen = $blen unless defined($plen);        

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

}

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

# T_fdiff(k,dz)

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

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

sub T_fdiff($$)

{

    my($k,$dz) =

        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <differencing interval[m]>',@_);

    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

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

sub T_interp($$$)

{

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

        &antsFunUsage(3,'fff','<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>',@_);

    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

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

sub T_binavg($$)

{

    my($k,$dz) =

        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <grid resolution[m]>',@_);

    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

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

sub T_tilt($$)

{

    my($k,$dprime) =

        &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <d-prime[m]>',@_);

    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]>',@_);

	croak("T_UH($k,$blen,$dz,$range_max): bad parameters\n")

		unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				

    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]>',@_);

	croak("T_ASM($k,$blen,$dz,$range_max): bad parameters\n")

		unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				

    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

#	- sel == nan disables pre-averaging correction

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

sub T_VI($$$$$)

{

	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));        

}

sub T_VI_alt($$$$$)

{

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

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

	croak("T_VI_alt($k,$blen,$sel,$dz,$range_max): bad parameters\n")

		unless ($k>=0 && $blen>0 && $sel ne '' && $dz>0 && $range_max>=0);				

	croak("$0: tilt-dprime outside range [0..$blen]\n")

		unless ($dprime>=0 && $dprime<=$blen);

    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 length[m]> <output grid resolution[m]> <range max[m]>',

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

		croak("T_w($k,$blen,$dz,$range_max): bad parameters\n")

			unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				

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

	}

}

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

# T_w_z(k,blen,dz,range_max)

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

#	- first differencing of gridded shear to calculate dw/dz

#	- NB: grid-scale differentiation assumed

#	- range_max == 0 disables tilt correction

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

{ my(@fc);

	sub T_w_z(@)

	{

		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,undef,@_);

		croak("T_w_z($k,$blen,$dz,$range_max): bad parameters\n")

			unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);				

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

	}

}

1;