#======================================================================
#                    P O L I N T . P L 
#                    doc: Thu Nov 23 20:38:46 2000
#                    dlm: Tue Aug  5 14:06:31 2008
#                    (c) 2000 A.M. Thurnherr
#                    uE-Info: 17 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# 2nd edition NR polint.c adapted to ANTS

# HISTORY:
#	Nov 23, 2000: - created for [.interp.poly]
#	Jan 12, 2006: - BUG: higher-order polynomials could not be used
#						 to interpolate linear function
#   Jul  1, 2006: - Version 3.3 [HISTORY]
#	Jul 28, 2006: - cosmetics
#	Aug  5, 2008: - BUG: [.interp.poly] takes data from ref, not @ants_

# NOTES:
#	- &vector()-allocated arrays are numbered from 1
#	- (nan,nan) is returned on non-numeric required @ants_ values
#	- in contrast to the NR routine, the error value returned is +ve

require "$ANTS/nrutil.pl";

sub polint($$$$$$)								# ($y,$dy) = &polint(...)
{
	my($dR,$xf,$xv,$ti,$n,$yf) = @_;
	my($y,$dy);

	my($i,$m); my($ns) = 1;
	my($den,$dif,$dift,$ho,$hp,$w);
	my(@c,@d);

	for ($i=0; $i<$n; $i++) {					# check for nans
		return (nan,nan)
			unless (numberp($dR->[$ti+$i][$xf]) &&
					numberp($dR->[$ti+$i][$yf]));
	}

	$dif = abs($xv - $dR->[$ti][$xf]);
	&vector(\@c,1,$n);
	&vector(\@d,1,$n);
	for ($i=1; $i<=$n; $i++) {
		$dift = abs($xv - $dR->[$ti+$i-1][$xf]);
		if ($dift < $dif) {
			$ns  = $i;
			$dif = $dift;
		}
		$c[$i] = $dR->[$ti+$i-1][$yf];
		$d[$i] = $dR->[$ti+$i-1][$yf];
	}
	$y = $dR->[$ti+$ns---1][$yf];				# WHAT A CONSTRUCT :-)
	for ($m=1; $m<$n; $m++) {
		for ($i=1; $i<=$n-$m; $i++) {
			$ho  = $dR->[$ti+$i-1][$xf] - $xv;
			$hp  = $dR->[$ti+$i+$m-1][$xf] - $xv;
			$w   = $c[$i+1] - $d[$i];
			$den = $ho - $hp;
### The following two lines of code are the original, which makes polint
### fail when interpolating a linear function with a higher-order polynomial,
### as is done in [ubtest/resample.TF]. 
###			croak("$0 (polint.pl): ERROR!") if ($den == 0);
###			$den   = $w / $den;
### The following line of code is the replacement that solves the bug.
			$den   = $w / $den unless ($den == 0);
			$d[$i] = $hp * $den;
			$c[$i] = $ho * $den;
		}
		$dy = (2*$ns < ($n-$m)) ? $c[$ns+1] : $d[$ns--];
		$y += $dy;
	}
	return ($y,abs($dy));
}

1;