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

#                    . I N T E R P . A D C P 

#                    doc: Fri Apr 16 16:07:48 2010

#                    dlm: Tue Aug  9 23:07:35 2011

#                    (c) 2010 A.M. Thurnherr

#                    uE-Info: 81 0 NIL 0 0 72 2 2 4 NIL ofnI

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

# interpolation scheme mimicking RDI ADCP response as described in

# Broadband primer, p.17

# HISTORY:

# 	Apr 16, 2010: - created

#	Aug  9, 2011: - added -u

# NOTES:

#	- interface is described in [.interp.linear]

$IS_opts = 'b:u';

$IS_optsUsage = '[pass -u)nfiltered] -b)in/pulse <length>';

sub IS_usage()

{

	&antsUsageError()

		unless defined($opt_b);

}

sub IS_init($$$$) {}

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#

# &IS_interpolate(br,idr,xf,xv,xi,f)	interpolate field f

#		br							data buffer reference

#		idr							init-data reference

#		xf							x field

#		xv							x value

#		xi							index of last record with x-value <= x

#		f							field number to interpolate

#		<ret val>					interpolated value

#

# NB:

#	- handle f == xf

#	- return NaN if any of the y values required is NaN

#

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub IS_interpolate($$$$$$)

{

	my($bR,$idR,$xf,$xv,$xi,$f) = @_;

	return $xv if ($xf == $f);							# return target x

	my($tow) = $xi;										# top of triangular sampling window

	while ($tow>=0 && (!numberp($bR->[$tow][$xf]) || $bR->[$tow][$xf] > $xv-$opt_b)) {

		$tow--;

	}

	if ($tow < 0) {										# incomplete window

		return nan unless ($opt_u);

		$tow = 0;

	} else {

		$tow++;

	}

	my($bow) = $xi+1;									# bottom of triangular sampling window

	while ($bow<=$#{$bR} && (!numberp($bR->[$bow][$xf]) || $bR->[$bow][$xf] < $xv+$opt_b)) {

		$bow++;

	}

	if ($bow > $#{$bR}) {								# incomplete window

		return nan unless ($opt_u);

		$bow = $#{$bR};

	} else {

		$bow--;

	}

	my($sweight) = 0;									# calculate weighted average

	my($sum) = 0;

	for (my($i)=$tow; $i<=$bow; $i++) {

		next unless (numberp($bR->[$i][$xf]) && numberp($bR->[$i][$f]));

		my($weight) = 1 - abs($bR->[$i][$xf]-$xv)/$opt_b;

		$sum += $weight * $bR->[$i][$f];

		$sweight += $weight;

	}

	return ($sweight>0) ? $sum/$sweight : nan;

}

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

1;