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

#                    . I N T E R P . S P L I N E 

#                    doc: Wed Nov 22 21:01:09 2000

#                    dlm: Tue Aug  5 14:20:39 2008

#                    (c) 2000 A.M. Thurnherr

#                    uE-Info: 31 0 NIL 0 0 72 10 2 4 NIL ofnI

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

# spline interpolation

# HISTORY:

# 	May 27, 2003: - adapted from [.interp.linear]

#	Jul 19, 2004: - BUG: i instead of $i --- 'orrible!

#	May 15, 2006: - BUG: assertion did not work on -e

#	Jul 28, 2006: - Version 3.3 [HISTORY]

#				  - made 100% compatible with [.interp.spline]

#				  - added xf to ISInit() args

#	Jul 30, 2006: - continued debugging

#   Aug 22, 2006: - adapted to work with [match]

#   Aug  5, 2008: - added idr param to IS_init()

# see [.interp.linear] for documentation of interface

# EXAMPLE:

#	Matlab: x = 0:10;  y = sin(x); xx = 0:.25:10;

#		yy = spline(x,y,xx); plot(x,y,'o',xx,yy)

# 	gnuplot/ANTS:

#		set style function points;

#		set style data lines;

#		plot [0:10] sin(x), \

#		'<Cat -f x:0,1,10 | list x y=sin($x) | resample -s spline x \#0-10:0.25'

$IS_opts = "B:";

$IS_optsUsage = "[-B)oundary cond <1st/last>]";

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

my($yp1,$ypn);

sub IS_usage()

{

	if (defined($opt_B)) {								# NR p.115

		($yp1,$ypn) = split('/',$opt_B);

		croak("$0: can't decode -B $opt_B\n")

			unless (numberp($yp1) && numberp($ypn));

	}

}

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

sub IS_init($$$$) {

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

	my($i,$k,$p,$qn,$sig,$un,@u);

	my($n) = scalar(@{$bR});

	if (defined($opt_B)) {								# handle boundary cond

		$idR->[1][$f] = -0.5;

		$u[1] = (3/($bR->[1][$xf]-$bR->[0][$xf])) *

				    (($bR->[1][$f]-$bR->[0][$f]) /

					 ($bR->[1][$xf]-$bR->[0][$xf]) - $yp1);

	} else {

		$idR->[1][$f] = $u[1] = 0;

	}

	for ($i=2; $i<=$n-1; $i++) {

		$sig = ($bR->[$i-1][$xf]-$bR->[$i-2][$xf]) /

				($bR->[$i][$xf]-$bR->[$i-2][$xf]);

		$p = $sig*$idR->[$i-1][$f] + 2;

		$idR->[$i][$f] = ($sig-1)/$p;

		$u[$i] = ($bR->[$i][$f]-$bR->[$i-1][$f]) /

					($bR->[$i][$xf]-$bR->[$i-1][$xf])

						- ($bR->[$i-1][$f]-$bR->[$i-2][$f]) /

							($bR->[$i-1][$xf]-$bR->[$i-2][$xf]);

		$u[$i] = (6*$u[$i]/($bR->[$i][$xf]-$bR->[$i-2][$xf]) -

					$sig * $u[$i-1]) / $p;

	}

	if (defined($opt_B)) {

		$qn = 0.5;

		$un = (3/($bR->[$n-1][$xf]-$bR->[$n-2][$xf])) *

				($ypn-($bR->[$n-1][$f]-$bR->[$n-2][$f]) /

					($bR->[$n-1][$xf]-$bR->[$n-2][$xf]));

    } else {

		$qn = $un = 0;

	}

	$idR->[$n][$f] = ($un-$qn*$u[$n-1]) / ($qn*$idR->[$n-1][$f]+1);

	for ($k=$n-1;$k>=1;$k--) {

		$idR->[$k][$f] = $idR->[$k][$f] * $idR->[$k+1][$f] + $u[$k];

	}

}

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

sub IS_interpolate($$$$$$)

{

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

	return $xv if ($xf == $f);

	return $bR->[$xi][$f]						# edge values are ok

		if equal($xv,$bR->[$xi][$xf]);

	return $bR->[$xi+1][$f]

		if equal($xv,$bR->[$xi+1][$xf]);

	return nan unless (numberp($bR->[$xi][$f]) &&

					   numberp($bR->[$xi+1][$f]));

	croak("$0: assertion $bR->[$xi+1][$xf] >= $xv >=  $bR->[$xi][$xf] failed")

		unless (defined($opt_e) || $bR->[$xi+1][$xf] >= $xv && $xv >= $bR->[$xi][$xf]);

	my($h) = $bR->[$xi+1][$xf] - $bR->[$xi][$xf];

	croak("$0: assertion #2 failed") unless ($h > 0);

	my($a) = ($bR->[$xi+1][$xf] - $xv) / $h;

	my($b) = ($xv - $bR->[$xi][$xf]) / $h;

	return $a*$bR->[$xi][$f] + $b*$bR->[$xi+1][$f] +

			(($a*$a*$a-$a) * $idR->[$xi+1][$f] +

			 ($b*$b*$b-$b) * $idR->[$xi+2][$f]) *

			 	($h*$h)/6;

}

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

1;