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

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

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

#                    dlm: Wed Jan 18 18:46:33 2012

#                    (c) 2011 A.M. Thurnherr

#                    uE-Info: 103 19 NIL 0 0 72 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

require "$ANTS/libvec.pl";

require "$ANTS/libfuns.pl";

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

# Polzin et al., JAOT 2002 LADCP shear corrections

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

# NOTES:

#	- apply to downcast data only

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

# individual corrections

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

# NB: Dzb = (Dzt == Dzr) assumed

sub T_ravg($$)

{

    my($kz,$Dzb) =

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

    return 1 / sinc($kz*$Dzb/2/$PI)**4;

}

sub T_fdiff($$)

{

    my($kz,$Dzd) =

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

    return 1 / sinc($kz*$Dzd/2/$PI)**2;

}

sub T_interp($$$)

{

    my($kz,$Dzb,$Dzg) =

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

    return 1 / sinc($kz*$Dzb/2/$PI)**4 / sinc($kz*$Dzg/2/$PI)**2;

}

# NB: Polzin et al claim that Dz should be ADCP bin size, which does not seem to make sense

sub T_binavg($$)

{

    my($kz,$Dzg) =

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

    return 1 / sinc($kz*$Dzg/2/$PI)**2;

}

sub T_tilt($$)

{

    my($kz,$dprime) =

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

    return 1 / sinc($kz*$dprime/2/$PI)**2;

}

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

# convenient combinations

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

sub LADCP_tilt_dprime($)

{

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

} 

sub T_UH($$$$)

{

	my($kz,$blen,$grez,$maxrange) =

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

    return T_ravg($kz,$blen) * T_fdiff($kz,$blen) * T_interp($kz,$blen,$grez) * T_tilt($kz,LADCP_tilt_dprime($maxrange));

}

sub T_SM($$$$)

{

	my($kz,$blen,$grez,$maxrange) =

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

    return T_ravg($kz,$blen) * T_fdiff($kz,$blen) * T_binavg($kz,$grez) * T_tilt($kz,LADCP_tilt_dprime($maxrange));

}

sub T_VI($$$$$)

{

	my($kz,$blen,$sel,$grez,$maxrange) =

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

	return T_VI_alt($kz,$blen,$sel,$grez,LADCP_tilt_dprime($maxrange));        

}

sub T_VI_alt($$$$$)

{

	my($kz,$blen,$sel,$grez,$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("$0: tilt-dprime outside range [0..$blen]\n")

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

    return ($sel>0) ? T_ravg($kz,$blen) * T_binavg($kz,$sel) * T_binavg($kz,$grez) * T_tilt($kz,$dprime)

    				: T_ravg($kz,$blen) * T_binavg($kz,$grez) * T_tilt($kz,$dprime);

}

1;