C+
C NAME:
C	FLINT
C PURPOSE:
C	Linear interpolation in N dimensions
C CALLING SEQUENCE:
	function FLINT(NN,ND,F,P,EPS)
C INPUTS:
C	NN	integer		N=abs(NN): # dimensions. If NN<0 then values in F equal
C				to BadR4() are considered 'invalid' (see PROCEDURE).
C	ND(N)	integer		size of each dimension. Up to 10 degenerate dimensions
C				of 1 element are allowed (e.g. F(10,1,10,1,10))
C	F(*)	real		array with function values on which to interpolate
C	P(N)	real		point in N-dim cube, using the array indices as coordinate axes
C	EPS	real		controls interpolation very close to grid lines (see PROCEDURE)
C OUTPUTS:
C	Result	real		the interpolated value
C CALLS:
C	BadR4, ArrI4Copy, ArrR4Copy, ArrR8Copy
C SEE ALSO:
C	FLINT8
C RESTRICTIONS:
C >	The coordinate values of P are best kept inside the range for each
C	dimension: 1<=P(i)<=ND(i), i.e. the point should lie 'inside'
C	the 'volume' covered by the array F.
C	Strictly speaking extrapolation is possible for one unit in each index:
C	as long as 0<=P(i)<=ND(i)+1 the result will be useful.
C	(in particular, for P(I)<0 the calculation of the principal corner
C	will not work properly). In this case the value BadR4() is returned.
C >	Currently 10 degerate dimensions are allowed. This is set as parameter N1MAX.
C	If there are more than N1MAX degenerate dimensions then the value BadR4()
C	is returned.
C PROCEDURE:
C >	The interpolating function is a linear function in all dimensions which
C	matches the function values F in all corners of the N-dimensional cube.
C	1-dim:	f = a+bx
C	2-dim:	f = a+bx+cy+dxy
C	3-dim:	f = a+bx+cy+dz+exy+fyz+gxz+hxyz
C	etc.
C	The constants are uniquely determined by the condition that the
C	function values in the corners must be matched:
C	1-dim:	f = (1-wx)F(0)+wx*F(1)
C	2-dim:	f = (1-wx)(1-wy)F(0,0)+(1-wx)wy*F(0,1)+
C		       wx*(1-wy)F(1,0)+   wx*wy*F(1,1)
C	3-dim:	f = (1-wx)(1-wy)(1-wz)F(0,0,0)+
C		    (1-wx)(1-wy)   wz*F(0,0,1)+
C		    (1-wx)   wy*(1-wz)F(0,1,0)+
C		    (1-wx)   wy*   wz*F(0,1,1)+
C		       wx*(1-wy)(1-wz)F(1,0,0)+
C		       wx*(1-wy)   wz*F(1,0,1)+
C		       wx*   wy*(1-wz)F(1,1,0)+
C		       wx*   wy*   wz*F(1,1,1)
C	etc. for higher dimensions.
C >	The array F is declared in the calling program as
C		real	F(ND(1),ND(2),..,ND(N))
C >	The array P defines a point in the N-dim cube where the interpolated
C	function value is needed.
C >	The principal corner NP is the array element nearest to P.
C	with all indices smaller than the coordinate values of P.
C >	The weights W are defined as P-NP. By definition W >= 0.
C >	Array NP selects an element in the N-dim array F:
C		F(NP(1)  ,NP(2)  ,..,NP(N))
C	This is one corner (the principal corner) of an N-dimensional cube.
C	The other corners are:
C		F(NP(1)+n1,NP(2)+n2,..,NP(N)+nn)
C	where n1,n2,..,nn take on values 0 and 1 (the principal corner
C	has n1=n2=..=nn=0). The N-dimensional cube has a total of 2^N corners.
C
C	Use the array indices as coordinate system:
C	Point (1,1,..,1) has function value F(1,1,..,1).
C	Point (ND(1),ND(2),..,ND(N)) has function value F(ND(1),ND(2),..,ND(N)).
C	etc. for all other elements in array F.
C
C	The weights W define a point P inside the N-dimensional cube.
C		P(NP(1)+W(1),NP(2)+W(P),..,NP(N)+W(N))
C	i.e. the weights W are the distances of P to the principal corner
C	along each of the coordinate axes.
C >	If NN<0 then, if one or more of the corners of the N-dim cube contain
C	the value BadR4() (defined in math.h), then the value BadR4() is
C	returned. Typically BadR4() is used to flag entries in F as 'invalid'.
C >	For degenerate dimensions the corresponding coordinate in P is ignored.
C >	If the condition 0 <= P(i) <= ND(i)+1 is violated on any of the
C	non-degenerate dimensions (ND(i) > 1) than the result BadR4() is
C	returned.
C
C >	The constant EPS influences how interpolation near grid lines is done.
C	If the point P is closer than EPS to a side of an N-dim cube
C	then P is assumed to lie exactly on the side of the cube.
C	This argument was added primarily to handle differences between
C	compilers (specifically Absoft and G77) when P is close to the side
C	of a cube. Typically EPS will be no larger than the machine precision,
C	e.g. as returned by href=TinyR4=. In all cases I have seen so far
C	setting EPS to a non-zero small value was unnecessary, and could have
C	been avoided with a little bit of careful coding. Preferably EPS is
C	set to zero.
C
C MODIFICATION HISTORY:
C	APR-1996, Paul Hick (UCSD/CASS)
C	MAY-1997, Paul Hick (UCSD/CASS)
C	    added BadR4() check
C	DEC-1999, Paul Hick (UCSD/CASS)
C	    added capacity to deal with degenerate dimensions; added FLINT8
C	DEC-2000, Paul Hick (UCSD/CASS)
C	    removed stop statement used to stop execution if there are more
C	    then N1MAX degenerate dimensions. Now the value BadR4() is returned.
C	AUG-2001, Paul Hick (UCSD/CASS)
C	    added a check for zero weight of any of the corners involved in the
C	    interpolation. If the weight is zero the corner now is ignored if
C	    the function value is bad (previously a bad corner would result in
C	    a bad return value). As a result if the point P is on a face of the
C	    N-dim cube only the corners of the face matter (i.e. an (N-1)-dim
C	    interpolation is done); if P is on an edge of the N-dim cube then an
C	    (N-2)-dim interpolation is done; etc.
C	AUG-2001, Paul Hick (UCSD/CASS)
C	    A bad value would be returned if any of the coordinates P(I) was
C	    more than one unit outside the range [1,ND(I)], even for degenerate
C	    dimensions. Now the coordinate P(I) for a degerate dimension is
C	    completely ignored (effectively it is set to P(I) = 1; the value is
C	    not explicitly used in the calculation).
C	SEP-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    Added argument EPS
C-
	    integer	NN
	    integer	ND(*)
	    real	F(*)
	    real	P(*)
	    real	EPS

	logical		bCheckBad

	parameter	(N1MAX=10)		! Max # of degenerate dimensions
	integer		ND1(N1MAX)
	real		P1 (N1MAX)

	real		Bad
	real		FI
	real		FP
	real		FF

	NP(I) = max(1,min(ND(I)-1,int(P(I))))	! Index I of principal corner

	N = abs(NN)
	bCheckBad = NN .lt. 0

	Bad = BadR4()

	!-------
	! Don't extrapolate more than one unit. Make sure that
	! 0 <= P(I) <= ND(I)+1 or that the dimension I is degenerate
	! (if ND(I) = 1 then the coordinate P(I) is not used).

	I = 1
	do while (I .le. N .and. ((0 .le. P(I) .and. P(I) .le. ND(I)+1) .or. ND(I) .eq. 1))
	    I = I+1
	end do

	if (I .le. N) then	! Invalid coordinate P(I)
	    FF = Bad

	else
	    N1 = 0
	    do I=1,N		! Check for degenerate dimensions (of 1 element)
		if (ND(I) .eq. 1) then
		    N1 = N1+1	! Count # degenerate dimensions
		    if (N1 .le. N1MAX) then
			ND1(N1) = I	! Position of degenerate dimension
			P1 (N1) = P(I)	! Degenerate coordinate
		    end if
		end if
	    end do

	    if (N1 .gt. N1MAX) then! Too many degenerate dimensions
		FF = Bad

	    else if (N .eq. N1) then! All dimensions degenerate:
		FF = F(1)	! .. only 1 element in F

	    else

		!-----
		! Remove the degenerate dimensions from the ND and P arrays, by
		! shifting the non-degenerate dimensions to a lower index.

		do I=N1,1,-1
		    IP = ND1(I)
		    if (IP .ne. N) then
			call ArrI4Copy(N-IP,ND(IP+1),ND(IP))
			call ArrR4Copy(N-IP,P (IP+1),P (IP))
		    end if
		    ND(N) = 0	! Not really necessary
		    P (N) = 0.0
		end do

		N = N-N1	! # non-degenerate dimensions

		!-----
		! ND(1:N) now contains the non-degenerate dimensions and NP(1:N) the
		! the matching coordinates. ND(N+1:*) and P(N+1:*) will not be used.

		!-----
		! In the calling program F is declared as an n-dimensional array.
		! Internally F is addressed as a one dimensional array.
		!
		! If F is an n-dimensional array F(ND1,ND2,..,NDn) then element F(I1,I2,..,In)
		! is located at one-dimensional index IP:
		!
		! (1-dim array)	 IP = 1+ I1-1
		! (2-dim array)	 IP = 1+ I1-1 +ND1*( I2-1 )
		! (3-dim array)	 IP = 1+ I1-1 +ND1*( I2-1 +ND2*(I3-1) )
		! (n-dim array)	 IP = 1+ I1-1 +ND1*( I2-1 +ND2*( + ... ND[n-1]*(In-1))...)

		IP = NP(N)-1	! Calculate offset for the principal corner:
		do I=N-1,1,-1
		    IW = NP(I)-1
		    IP = IW+ND(I)*IP
		end do

		IP0 = IP+1	! One-dimensional index of principal corner

		!-----
		! The n-dim cube has 2^n corners, identified by ICN=0,2^n-1.
		! The array index of the principal corner in the i-th dimension: NP(I)
		! The index for the other corners in the i-th dimension  : NP(I)+0 or NP(I)+1
		! The offset 0 or 1 is calculated according to the following scheme
		! (example is for a 3-dim array).
		!
		!           0               1          x = 1st dim
		!         /   \           /   \
		!       0       1       0       1      y = 2nd dim
		!     /   \   /   \   /   \   /   \
		!     0   1   0   1   0   1   0   1    z = 3rd dim
		!
		!     0   1   2   3   4   5   6   7      = ICN
		!
		! E.g. corner ICN=5 has array indices NP(1)+1,NP(2)+0,NP(3)+1
		! The principal corner ICN=0 has indices NP(1),NP(2),NP(3)

		FF  = 0.0
		ICN = -1
		NCN = 2**N-1

		do while (ICN .lt. NCN)	! Loop over all ICN corners of the N-dim cube around P
		    ICN = ICN+1

		    IC = 2*ICN
		    FI = 1.0
		    IP = 0

		    do I=N,1,-1		! Loop over dimensions N-1 to 1
			IC = IC/2
			IW = mod(IC,2)	! Offset 0/1 for I-th dim
					! Distance to principal corner =
			DI = P(I)-NP(I)	! .. weight 1-DI for offset 0; DI for offset 1

			if (abs(DI  ) .lt. EPS) DI = 0.0
			if (abs(DI-1) .lt. EPS) DI = 1.0

			FI = FI*((1-IW)*(1-DI)+IW*DI)	! Total weight for corner ICN

			IP = IW+ND(I)*IP! Total 1-dim of corner ICN relative
					! .. to the principal corner ICN=0
		    end do

		    !------
		    ! If weight for corner ICN is zero, then the function value
		    ! does not matter and we can skip the corner. Note that we
		    ! check for FI = 0.0 before we check for bad values.

		    if (FI .ne. 0.0) then
			FP = F(IP0+IP)	! If bad value, return with Bad
			if (bCheckBad .and. FP .eq. Bad) then
			    ICN = NCN	! Break do while loop
			    FF  = Bad	! Return bad value
			else
			    FF = FF+FI*FP	! Accumulate interpolation
			end if
		    end if
		end do

		!-------
		! If there were degenerate dimensions we still need to restore
		! the arrays ND and P to their input values, using the arrays
		! ND1 and P1.

		N = N+N1		! Original # dimensions

		do I=1,N1		! Loop over degenerate dimensions
		    IP = ND1(I)
		    if (IP .ne. N) then	! Make room for insertion by shifting right
			call ArrI4Copy(N-IP,ND(IP),ND(IP+1))
			call ArrR4Copy(N-IP,P (IP),P (IP+1))
		    end if
		    ND(IP) = 1		! Restore input values
		    P (IP) = P1(I)
		end do

	    end if

	end if

	FLINT = FF		! Interpolated value (maybe bad)

	return
	end