SUBROUTINE InterpolateIRC( x, f, g, iPower, xTab, fTab, gTab, iTab, NkTab )
! Internal reflection coefficient interpolator.
! Returns f, g, iPower for given x using tabulated values.
! Uses polynomial interpolation to approximate the function between the tabulated values
USE PolyMod
IMPLICIT NONE
INTEGER, PARAMETER :: N = 3 ! order of the polynomial for interpolation
INTEGER, INTENT( IN ) :: NkTab, iTab( NkTab )
REAL ( KIND=8 ), INTENT( IN ) :: xTab( NkTab )
COMPLEX ( KIND=8 ), INTENT( IN ) :: fTab( NkTab ), gTab( NkTab ), x
INTEGER, INTENT( OUT ) :: iPower
COMPLEX ( KIND=8 ), INTENT( OUT ) :: f, g
INTEGER :: i, j, iDel, iLeft, iMid, iRight, NAct
REAL ( KIND=8 ) :: xReal
COMPLEX ( KIND=8 ) :: xT( N ), fT( N ), gT( N )
xReal = DBLE( x ) ! taking the real part
iLeft = 1
iRight = NkTab
! Three cases: x left, in, or right of tabulated interval
IF ( xReal < xTab( iLeft ) ) THEN
f = fTab( iLeft )
g = gTab( iLeft )
iPower = iTab( iLeft )
ELSE IF( xReal > xTab( iRight ) ) THEN
f = fTab( iRight )
g = gTab( iRight )
iPower = iTab( iRight )
ELSE
! Search for bracketing abscissas:
! Log base 2 (NPts) stabs required for a bracket
DO WHILE ( iLeft /= iRight-1 )
iMid = ( iLeft + iRight )/2
IF ( xTab( iMid ) > xReal ) THEN
iRight = iMid
ELSE
iLeft = iMid
ENDIF
END DO
! Extract the subset for interpolation and scale
iLeft = MAX( iLeft - ( N - 2 ) / 2, 1 )
iRight = MiN( iRight + ( N - 1 ) / 2, NkTab )
NAct = iRight - iLeft + 1
DO i = 1, NAct
j = i + iLeft - 1
iDel = iTab( j ) - iTab( iLeft )
xT( i ) = xTab( j )
fT( i ) = fTab( j ) * 10.0D0 ** iDel
gT( i ) = gTab( j ) * 10.0D0 ** iDel
END DO
! Construct the interpolate
f = Poly( x, xT, fT, NAct )
g = Poly( x, xT, gT, NAct )
iPower = iTab( iLeft )
ENDIF
END SUBROUTINE InterpolateIRC
