#======================================================================
#                    . N M I N T E R P . L I N E A R 
#                    doc: Wed Nov 22 21:01:09 2000
#                    dlm: Fri Aug 30 13:26:58 2019
#                    (c) 2000 A.M. Thurnherr
#                    uE-Info: 101 9 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================

# linear interpolation for [fillgaps]

# HISTORY:
#	Apr  6, 2006: - created from [.interp.linear]
#	Jul 28, 2006: - added xf to ISInit() args
#	Jul 30, 2006: - added max-gap support
#	Aug  8, 2008: - documentation
#	Sep 23, 2011: - added support for %RECNO
#				  - removed xv from interp() args
#	Oct 19, 2011: - code did not work for %RECNO and trailing missing vals
#	Jan 23, 2014: - implemented huge speedup for sparse files

# NOTES:
#	- in contrast to the [.interp.*] routines, the [.nminterp.*] routines:
#		1) do not assume that x values increase monotonically
#		2) can handle nans
#		3) cannot be used with multiple buffers
#		4) use %RECNO when xfnr==-1

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
#
#	$ISOpts				string of allowed options
#	$ISOptsUsage		usage information string for options
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

$ISOpts = "";
$ISOptsUsage = "";

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &ISUsage() 	mangle parameters (options, really)
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub ISUsage() {}

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &ISInit(f,xf)				init interpolation of field f
#       f					field number
#		xf					x field number or -1 for %RECNO
#       <ret val>           none
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub ISInit($$) {}

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &interpolate(xf,mg,f,lvr,cr)		interpolate field f
#		xf							x field or -1 for %RECNO
#		mg							max allowed x gap
#		f							field number to interpolate
#		lvr							last record with valid y val
#		cr							current record
#		<ret val>					interpolated value
#
# NB:
#	- handle f == xf
#	- return undef if interpolation cannot be carried out
#	- x VALUES MAY NOT BE MONOTONIC
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

{ my(@nvr);		# static: next valid record (per field)

	sub interpolate($$$$$$)
	{
		my($xf,$mg,$f,$lvr,$cr) = @_;
	
		return $ants_[$cr][$xf] if ($xf == $f);
		return $ants_[$lvr][$f] 				# edge value is ok
			if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$lvr][$xf])) ||
			   ($xf<0 && $cr==$lvr);
	
		unless ($nvr[$f] > $cr) {
			$nvr[$f] = $cr + 1;
			while ($nvr[$f] <= $#ants_ && !numberp($ants_[$nvr[$f]][$f])) { $nvr[$f]++; }
		}
	
		return undef if ($nvr[$f] > $#ants_);
		return undef
			if ($xf>=0 && abs($ants_[$nvr[$f]][$xf]-$ants_[$lvr][$xf])>$mg);
		return $ants_[$nvr[$f]][$f] 				# edge value is ok
			if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$nvr[$f]][$xf])) ||
			   ($xf<0 && $cr==$nvr[$f]);
	
		if ($xf < 0) {
			die("$nvr[$f] - $lvr") unless ($nvr[$f] - $lvr > 0);
        }
		my($sc) = ($xf >= 0)
				? ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / ($ants_[$nvr[$f]][$xf] - $ants_[$lvr][$xf])
				: ($cr - $lvr) / ($nvr[$f] - $lvr);
		return $ants_[$lvr][$f] + $sc * ($ants_[$nvr[$f]][$f] - $ants_[$lvr][$f]);
	}
} # static scope

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

1;