!
!                                **********************************************
!                                *             MODULE pbint                   *
!                                *  R. J. Purser, NOAA/NCEP/EMC 23rd May 2012 *
!                                *  jim.purser@noaa.gov                       *
!                                *                                            *
!                                **********************************************
!
! A suite of simple-to-use 1D interpolation routines for both regular and
! irregular monotonic grids.
!
! DIRECT DEPENDENCIES
! Libraries[their Modules]: pcomb[pcomb], pbend[pbend]
! Additional Modules      : pietc,pkind
!
! The methods provided here can all be considered to be generalizations, of
! various kinds, of the Lagrange polynomial interpolation method. Here, the
! interpolation is considered to be "applied" for a given target point when
! the stencil of interpolation weights is computed, the evaluation then
! being a simple "dot_product" operation over the designated sub-array of
! contiguous source values.
!
! For each kind of method, a uniform (unit-spaced) grid version and a 
! general non-uniform grid version is provided, but they then share the same
! module procedure interface names. The auxiliary arrays and constants
! for the uniform-grid versions are set up automatically, making them
! particularly easy to apply. For the generic-grid versions, these auxiliary
! arrays must be set up explicitly by the user before the application of
! the interpolation routines themselves.
!
! Broadly speaking, three kinds of interpolation are accommodated here:
! (1) Lagrange polynomial interpolation (without special measures to
!     preserve continuity of derivatives across source grid nodes).
!     Auxiliary arrays needed for the case of a general grid are precomputed
!     once by the appropriate "call getq"; the application to each particular
!     target point thereafter involves the appropriate "call getint" for
!     each target.
! (2) Lagrangian polynomial pairs combiined by complementary smooth blending
!     weights that go between 0 and 1 (and vice-versa) within the central
!     interval of the composite stencil. A "shape" parameter, <= 6, determines
!     the form of the sigmoidal blending function within the interval: for
!     positive integer shape, the blending function is a symmetric-parameter
!     incomplete Beta function, a polynomial of degree 2*shape-1; also
!     provided is a special shape=0 blending based on a Whittaker-function
!     smoothing formula which is formally continuous in ALL derivatives at
!     both ends of the interval. While this shape=0 and other relatively
!     large choices of shape ensure that several derivatives remain continuous,
!     the constraint that the blending is confined to within the single
!     central interval and only involves mixing two Lagrange formulae means
!     that, in many practical applications, the resulting interpolated
!     graphs appears to "lurch" in the middle of the interval from one
!     tentative solution to the next Lagrange-associated solution. For this
!     reason a third option is supplied.
!     Auxiliary arrays for the case of the general (i.e., nonuniform) grid 
!     are precomputed once by the appropriate "call getq". The application to 
!     each particular target requires the repeated appropriate "call getint".
! (3) Lagrange polynomials over possibly more than two consecutive ranges
!     are combined with smooth spline-associated weightings, which also can
!     span more than a single central interval of the composite stencil. Now,
!     in place of the "shape" parameter, a single "half-stencil" width
!     parameter, nh, is provided to determine both the order of formal
!     accuracy and the degree of smooth blending (which both increase with nh).
!     As the most trivial example, nh=1 denotes piecewise linear interpolation
!     (no blending). Nh=2, like shape=1 of the previous class, denotes 
!     linearly-weighted blending over each single central interval, but the
!     Lagrange polynomials blended are ALWAYS quadratics, i.e., the 
!     composite stencil ends up with nh*2=4 points (a cubic polynomial).
!     Generically, nh Lagrangian polynomials each of degree nh are blended
!     over the (nh-1) central intervals using blending weights obtained as
!     the partition of unity of piecewise (nh-1)-degree polynomial splines
!     that go as smoothly as possible beteen 0 and 1 over (nh-1) intervals.
!     This blending of Lagrange polynomials, being spread over possibly
!     several intervals, leads to very smooth and "natural" looking
!     interpolations, though possibly with more exuberant excursions between
!     source nodes that the simpler methods exhibit.
!     The auxiliary arrays for the general non-uniform grid require an initial
!     single "call getq" with the appropriate argument list. Then each
!     subsequent spline-weigted interpolation is accomplished by a
!     "call getsplint".
!
!     As a reminder, no "call getq" is needed when the source grid is uniform;
!     in this case we scale the coordinate to make the grid appear uniform
!     and the auxiliary arrays are automatically precomputed. However, as a
!     sop to the puritans of digital hygene, a "finalization" procedure
!     invoked by a "call finibint" after all the uniform-grid interpolations
!     have been complete deallocates the auxiliary arrays that were created
!     automatically. (Unfortunately, finibint cannot be invoked automatically
!     since there is no way to tell when the user has finished interpolating.)
!     
!=============================================================================
module pbint
!=============================================================================
use pkind, only: dp
use pietc, only: T,F,zero,one
implicit none
private
public:: finibint,getsplint,invalidgsplint,getint,getq
integer, allocatable,dimension(:,:,:):: tableaub
real(dp),allocatable,dimension(:)    :: xs,qq
integer                              :: nxs,nxsm,nxsp
data nxs/0/
interface finibint; module procedure finibint;                   end interface
interface inibint;  module procedure inibint;                    end interface
interface invalidgsplint; module procedure invalidgsplint;       end interface
interface getsplint
   module procedure getsplint,getsplintd,getgsplint,getgsplintd; end interface
interface getint
   module procedure getbint,getbintd,getlint,getlintd,getlap,getlapd,&
                    getgbint,getgbintd
                                                                 end interface
interface getq
   module procedure getq,getgq1,getgq2,getgsplqy;                end interface
interface btableau;  module procedure btableau;                  end interface
interface bsigmoid;  module procedure bsigmoid;                  end interface
interface binprod;   module procedure binprod;                   end interface
interface evaltay;   module procedure evaltay,evaltayd;          end interface
contains

!=============================================================================
subroutine finibint!                                                [finibint]
!=============================================================================
if(allocated(xs)       )deallocate(xs)
if(allocated(qq)       )deallocate(qq)
if(allocated(tableaub) )deallocate(tableaub)
nxs=0
end subroutine finibint

!=============================================================================
subroutine inibint(nxsa)!                                            [inibint]
!=============================================================================
! nxs is the number of unit intervals in the grid xs(0:nxs) used for all
! the default uniform-grid Lagrange interpolation formulae.
! qq is the corresponding array of Lagrange denominators.
! tableaub is the integer array representation of the polynomial
! coefficients of the unit-grid B-splines.
! Note: these arrays are not used when the grid is non-uniform.
!=============================================================================
integer,intent(in):: nxsa
!-----------------------------------------------------------------------------
integer           :: i,j
real(dp)          :: p
!=============================================================================
call finibint
nxs=nxsa; nxsp=nxs+1; nxsm=nxs-1
allocate( xs(0:nxs), qq(0:nxs) )
allocate( tableaub(0:nxs-1,0:nxs-1,0:nxs-1) )

! Set up a unit grid of nxs intervals:
do i=0,nxs; xs(i)=i; enddo

! Unit-grid B-spline tableau for sigmoid spanning nxs-1 intervals:
call btableau(nxsm,tableaub)

! unit-grid q factors for Laplace interpolation over nxs intervals:
call getq(nxsp,xs,qq)
end subroutine inibint

!============================================================================
subroutine getsplint(nh,x,w)!                                     [getsplint]
!============================================================================
! Apply spline-weighted interpolation of half-stencil width nh to a target
! located at a proportion, x, of the distance within the central interval.
! Return the stencil of nh*2 weights, w.
!============================================================================
use pcomb, only: factorial
integer,                 intent(in ):: nh
real(dp),                intent(in ):: x
real(dp),dimension(nh*2),intent(out):: w
!----------------------------------------------------------------------------
real(dp),dimension(0:nh,0:nh-1)  :: wL
real(dp),dimension(0:nh-1)       :: ws
real(dp),dimension(0:nh-1,0:nh-1):: co
real(dp)                         :: facn,xt,faci
integer                          :: i,j,k,n,nhm,nhp
!============================================================================
if(nh/=nxs)call inibint(nh)
n=nh*2; nhm=nh-1; nhp=nh+1
facn=factorial(nhm)
do k=0,nhm
   faci=one
   do i=0,nhm
      co(i,k)=tableaub(i,0,k)*faci/facn
      faci=faci*(i+1)
   enddo
enddo

w=zero
do k=0,nhm
   xt=x+k
   call getint(nhp,qq,xs,xt,wL(:,k))
   call evaltay(nhm,co(:,k),x,ws(k)) 
   w(n-k-nh:n-k)=w(n-k-nh:n-k)+ws(k)*wL(:,k)
enddo
end subroutine getsplint

!=============================================================================
subroutine getsplintd(nh,x,w,wd)!                                  [getsplint]
!=============================================================================
! Like getsplint, but also return the derivative stencil wd, assuming the
! source grid has unit spacing.
!=============================================================================
use pcomb, only: factorial
integer,                 intent(in ):: nh
real(dp),                intent(in ):: x
real(dp),dimension(nh*2),intent(out):: w,wd
!----------------------------------------------------------------------------
real(dp),dimension(0:nh,0:nh-1)  :: wL,wLd
real(dp),dimension(0:nh-1)       :: ws,wsd
real(dp),dimension(0:nh-1,0:nh-1):: co
real(dp)                         :: facn,xt,faci
integer                          :: i,j,k,n,nhm,nhp
!============================================================================
if(nh/=nxs)call inibint(nh)
n=nh*2; nhm=nh-1; nhp=nh+1
facn=factorial(nhm)
do k=0,nhm
   faci=one
   do i=0,nhm
      co(i,k)=tableaub(i,0,k)*faci/facn
      faci=faci*(i+1)
   enddo
enddo

w =zero
wd=zero
do k=0,nhm
   xt=x+k
   call getint(nhp,qq,xs,xt,wL(:,k),wLd(:,k))
   call evaltayd(nhm,co(:,k),x,ws(k),wsd(k)) 
   w(n-k-nh:n-k)=w(n-k-nh:n-k)+ws(k)*wL(:,k)
   wd(n-k-nh:n-k)=wd(n-k-nh:n-k)+wsd(k)*wL(:,k)+ws(k)*wLd(:,k)
enddo
end subroutine getsplintd

!=============================================================================
subroutine getgsplint(L,m,nh,qss,yss,xss, x,Lw,mw,wt)!             [getsplint]
!=============================================================================
! Get general (nonuniform) grid spline-weighted interpolation coefficients, wt.
! L,m  : grid point index limits
! nh   : half the stencil width and degree of component Lagrange polynomials
! qss  : Table of Lagrange denominators (precomputed via getq)
! yss  : Table of Taylor coefficients of spline segments (also from getq)
! xss  : Coordinates of all the grid points.
! x    : Coordinate of interpolation target point.
! Lw,mw: range of the weight stencil needed for this target
! wt   : Array whose first mw+1-Lw elements contain the stencil weights.
!=============================================================================
integer,                              intent(in ):: L,m,nh
real(dp),dimension(0:nh,L:m),         intent(in ):: qss
real(dp),dimension(0:nh-1,0:nh-1,L:m),intent(in ):: yss
real(dp),dimension(L:m),              intent(in ):: xss
real(dp),                             intent(in ):: x
integer,                              intent(out):: Lw,mw
real(dp),dimension(0:nh*2-1),         intent(out):: wt
!-----------------------------------------------------------------------------
integer                 :: i,i0,i1,ia,ib,ic,mh,nhm,nhp,n
real(dp)                :: rx,wc,wspl,wct
real(dp),dimension(0:nh):: wtt
!=============================================================================
nhm=nh-1; nhp=nh+1; n=nh*2
ia=L; ib=m
do
   if(ib-ia==1)exit
   ic=(ia+ib)/2
   if(x<xss(ic))then; ib=ic; else; ia=ic; endif
enddo
Lw=max(L,ib-nh); Mw=min(m,ia+nh)
wt=zero

! Express location as distance to "right" of grid point ia:
rx=x-xss(ia)

! Loop over spline weights:
wc=one
i0=max(L,ib-nh); i1=min(m-nh,ia)
do i=i0,i1
   if(i==i1)then
      wspl=wc
   else
      call evaltay(nhm,yss(:,ia-i,i+1),rx,wct)
      wspl=wc-wct
      wc=wct
   endif
   call getint(nhp,qss(:,i),xss(i:i+nh),x,wtt)
   wt(i-i0:i-i0+nh)=wt(i-i0:i-i0+nh)+wspl*wtt
enddo
end subroutine getgsplint
   
!=============================================================================
subroutine getgsplintd(L,m,nh,qss,yss,xss, x,Lw,mw,wt,wtd)!        [getsplint]
!=============================================================================
! Same as getgsplint, but also compute the derivative weight stencil, wtd.
!=============================================================================
integer,                              intent(in ):: L,m,nh
real(dp),dimension(0:nh,L:m),         intent(in ):: qss
real(dp),dimension(0:nh-1,0:nh-1,L:m),intent(in ):: yss
real(dp),dimension(L:m),              intent(in ):: xss
real(dp),                             intent(in ):: x
integer,                              intent(out):: Lw,mw
real(dp),dimension(0:nh*2-1),         intent(out):: wt,wtd
!-----------------------------------------------------------------------------
integer                 :: i,i0,i1,ia,ib,ic,mh,nhm,nhp,n
real(dp)                :: rx,wc,wcd,wspl,wspld,wct,wctd
real(dp),dimension(0:nh):: wtt,wttd
!=============================================================================
nhm=nh-1; nhp=nh+1; n=nh*2
ia=L; ib=m
do
   if(ib-ia==1)exit
   ic=(ia+ib)/2
   if(x<xss(ic))then; ib=ic; else; ia=ic; endif
enddo
Lw=max(L,ib-nh); Mw=min(m,ia+nh)
wt =zero
wtd=zero

! Express location as distance to "right" of grid point ia:
rx=x-xss(ia)

! Loop over spline weights:
wc=one
wcd=zero
i0=max(L,ib-nh); i1=min(m-nh,ia)
do i=i0,i1
   if(i==i1)then
      wspl =wc
      wspld=wcd
   else
      call evaltayd(nhm,yss(:,ia-i,i+1),rx,wct,wctd)
      wspl =wc -wct
      wspld=wcd-wctd
      wc =wct
      wcd=wctd
   endif
   call getint(nhp,qss(:,i),xss(i:i+nh),x,wtt,wttd)
   wt (i-i0:i-i0+nh)=wt (i-i0:i-i0+nh)+wspl *wtt
   wtd(i-i0:i-i0+nh)=wtd(i-i0:i-i0+nh)+wspld*wtt+wspl*wttd
enddo
end subroutine getgsplintd
   
!=============================================================================
subroutine invalidgsplint(L,m,nh,xss,flag)!                   [invalidgsplint]
!=============================================================================
! Flag this grid xss(L:m) if there is any feature about it that precludes
! its use in the spline-weighted interpolation scheme of index nh.
!=============================================================================
integer,                intent(in ):: L,m,nh
real(dp),dimension(L:m),intent(in ):: xss
Logical,                intent(out):: flag
!-----------------------------------------------------------------------------
integer:: i
!=============================================================================
flag=nh<1;     if(flag)return ! nh must be positive
flag=m-L<1;    if(flag)return ! Interval count must be positive
flag=m-L<nh+1; if(flag)return ! Must be elbow-room for spline weighting
do i=L,m-1
   flag=(xss(i)>=xss(i+1)); if(flag)return ! monotonicity and increasing xss
enddo
end subroutine invalidgsplint

!=============================================================================
subroutine getbint(shape,npth,rx, wt)!                                [getint]
!=============================================================================
! Compute the blended-interpolation weights for a target located at a fraction
! rx along the central interval of npth*2 uniformly spaced source points
! indexed [1-npth,npth]. The weights are returned in array wt.
! The integer "shape" parameter decides the mode of blending; 0 implies
! Whittaker function blending (arbitrarily differentiable); 1--6 imply
! beta function polynomial blending modes at successively higher orders.
!=============================================================================
use pbend, only: wbend
integer,                          intent(IN ):: shape,npth
real(dp),                         intent(IN ):: rx
real(dp),dimension(1-npth:npth),  intent(OUT):: wt
!-----------------------------------------------------------------------------
real(dp)                         :: wb
real(dp),dimension(1-npth:npth-1):: wta
real(dp),dimension(2-npth:npth  ):: wtb
integer                          :: npthm,npthmm
!=============================================================================
if(npth*2-2/=nxs)call inibint(npth*2-2)
select case(shape)
case(0); call wbend(rx,wb)
case(1); wb=rx
case(2); wb=rx**2*(3-2*rx)
case(3); wb=rx**3*(10+rx*(-15+rx*6))
case(4); wb=rx**4*(35+rx*(-84+rx*(70-rx*20)))
case(5); wb=rx**5*(126+rx*(-420+rx*(540+rx*(-315+rx*70))))
case(6); wb=rx**6*(462+rx*(-1980+rx*(3465+rx*(-3080+rx*(1386-rx*252)))))
case default; stop 'In getbint; shape code not recognized'
end select
npthm=npth-1; npthmm=npth-2
call getlap(nxsp,qq,xs,rx+npthm ,wta)
call getlap(nxsp,qq,xs,rx+npthmm,wtb)
wt=zero
wt( -npthm:npthm)=           (one-wb)*wta
wt(-npthmm:npth )=wt(-npthmm:npth)+wb*wtb
end subroutine getbint

!=============================================================================
subroutine getbintd(shape,npth,rx, wt,wtd)!                           [getint]
!=============================================================================
! Like getbint, but also provides the weights wtd for the derivative of the
! target (assuming the source grid to have unit spacing)
!=============================================================================
use pbend, only: wbend
integer,                          intent(IN ):: shape,npth
real(dp),                         intent(IN ):: rx
real(dp),dimension(1-npth:npth),  intent(OUT):: wt,wtd
!-----------------------------------------------------------------------------
real(dp)                         :: wb,wbd,b
real(dp),dimension(1-npth:npth-1):: wta,wtad
real(dp),dimension(2-npth:npth)  :: wtb,wtbd
integer                          :: npthm,npthmm
!=============================================================================
if(npth*2-2/=nxs)call inibint(npth*2-2)
b=rx*(one-rx)
select case(shape)
case(0); call wbend(rx,wb,wbd)
case(1); wb=rx;                                             wbd=one
case(2); wb=rx**2*(3-2*rx);                                 wbd=6*b
case(3); wb=rx**3*(10+rx*(-15+rx*6));                       wbd=30*b**2
case(4); wb=rx**4*(35+rx*(-84+rx*(70-rx*20)));              wbd=140*b**3
case(5); wb=rx**5*(126+rx*(-420+rx*(540+rx*(-315+rx*70)))); wbd=630*b**4
case(6); wb=rx**6*(462+rx*(-1980+rx*(3465+rx*(-3080+rx*(1386-rx*252)))))
                                                            wbd=2772*b**5
case default; stop 'In getbint; shape code not recognized'
end select
npthm=npth-1; npthmm=npth-2
call getlapd(nxsp,qq,xs,rx+npthm ,wta,wtad)
call getlapd(nxsp,qq,xs,rx+npthmm,wtb,wtbd)
wt =zero
wtd=zero
wt( -npthm:npthm) =              (one-wb)*wta
wt(-npthmm:npth ) =wt(-npthmm:npth)   +wb*wtb
wtd( -npthm:npthm)=                  -wbd*wta + (one-wb)*wtad
wtd(-npthmm:npth )=wtd(-npthmm:npth )+wbd*wtb +       wb*wtbd
end subroutine getbintd
!=============================================================================
subroutine getgbint(shape,L,M,npth,qss,xss,x, Lw,Mw,wt)!              [getint]
!=============================================================================
! General blended interpolation for a nonuniform grid and a single target.
!
! Assume the denominators for all the interior length-(2*npth-1) 
! Lagrange interpolation formulae are given in array qss, together with the 
! grid points, xss(L:M), themselves. 
! By blending consecutive pairs of the Lagrange formulae whose stencil
! centers staddle the nominal interval assigned to the target, a wider
! weight stencils, wt, of width 2*npth can be constructed, that ensures
! continuity of 1st derivative of the target value wrt the target location, x.
! The lowest and highest grid index for this wider weight stencil are Lw,Mw.
! The integer "shape" parameter decides the mode of blending; 0 implies
! Whittaker function blending (arbitrarily differentiable); 1--6 imply
! beta function polynomial blending modes at successly higher orders.
!=============================================================================
use pbend, only: wbend
integer,                                       intent(IN ):: shape,L,M,npth
real(dp),dimension(npth*2-1,L+npth-1:M+1-npth),intent(IN ):: qss
real(dp),dimension(L:M),                       intent(IN ):: xss
real(dp),                                      intent(IN ):: x
integer,                                       intent(OUT):: Lw,Mw
real(dp),dimension(npth*2),                    intent(OUT):: wt
!-----------------------------------------------------------------------------
real(dp)                         :: wa,wb,rx
real(dp),dimension(1-npth:npth-1):: wta
real(dp),dimension(2-npth:npth  ):: wtb
integer                          :: ia,ib,ic,npt,nptm
!=============================================================================
! Check for elbow room:
if(M+1-L<2*npth)stop 'In getgbint; Grid too small to accommodate stencil'

! Determine the appropriate interval to assign to target x:
ia=L+npth-1; ib=m+1-npth
do
   if(ib-ia==1)exit
   ic=(ia+ib)/2
   if(x<xss(ic))then; ib=ic; else; ia=ic; endif
enddo
Lw=ib-npth
Mw=ia+npth

! Express the location as a fractional distance along the nominal interval:
rx=(x-xss(ia))/(xss(ib)-xss(ia))

! Obtain the blending weight pair, wa and wb, according to the value of rx
! and the choice of shape parameter:

if(rx<=zero)then;    wa=one;  wb=zero ! <- extreme left
elseif(rx>=one)then; wa=zero; wb=one  ! <- extreme right
else                                  ! Interior:
   select case(shape)
   case(0); call wbend(rx,wb)
   case(1); wb=rx
   case(2); wb=rx**2*(3-2*rx)
   case(3); wb=rx**3*(10+rx*(-15+rx*6))
   case(4); wb=rx**4*(35+rx*(-84+rx*(70-rx*20)))
   case(5); wb=rx**5*(126+rx*(-420+rx*(540+rx*(-315+rx*70))))
   case(6); wb=rx**6*(462+rx*(-1980+rx*(3465+rx*(-3080+rx*(1386-rx*252)))))
   case default; stop 'In getgbint; shape code not recognized'
   end select
   wa=one-wb
endif

! Obtain the pair of Lagrange weight stencils, wta and wtb, centered on
! grid points, ia and ib:
npt=npth*2; nptm=npt-1
call getint(nptm,qss(:,ia),xss(Lw:Mw-1),x,wta)
call getint(nptm,qss(:,ib),xss(Lw+1:Mw),x,wtb)

! Combine the Lagrange weight stencils with blending weights, wa and wb,
! to obtain the final weight stencil, wt:
wt=zero; wt(1:nptm)=wa*wta; wt(2:npt )=wt(2:npt)+wb*wtb
end subroutine getgbint

!=============================================================================
subroutine getgbintd(shape,L,M,npth,qss,xss,x, Lw,Mw,wt,wtd)!         [getint]
!=============================================================================
! Like getgbint, but also output the derivative, wtd, at the target
!=============================================================================
use pbend, only: wbend
integer,                                       intent(IN ):: shape,L,M,npth
real(dp),dimension(npth*2-1,L+npth-1:M+1-npth),intent(IN ):: qss
real(dp),dimension(L:M),                       intent(IN ):: xss
real(dp),                                      intent(IN ):: x
integer,                                       intent(OUT):: Lw,Mw
real(dp),dimension(npth*2),                    intent(OUT):: wt,wtd
!-----------------------------------------------------------------------------
real(dp)                         :: wa,wb,rx,wad,wbd,b
real(dp),dimension(1-npth:npth-1):: wta,wtad
real(dp),dimension(2-npth:npth)  :: wtb,wtbd
integer                          :: ia,ib,ic,npt,nptm
!=============================================================================
! Check for elbow room:
if(M+1-L<2*npth)stop 'In getgbint; Grid too small to accommodate stencil'

! Determine the appropriate interval to assign to target x:
ia=L+npth-1; ib=m+1-npth
do
   if(ib-ia==1)exit
   ic=(ia+ib)/2
   if(x<xss(ic))then; ib=ic; else; ia=ic; endif
enddo
Lw=ib-npth
Mw=ia+npth

! Express the location as a fractional distance along the nominal interval:
rx=(x-xss(ia))/(xss(ib)-xss(ia))

! Obtain the blending weight pair, wa and wb, according to the value of rx
! and the choice of shape parameter:
if(rx<=0)then;       wa=one; wb=zero; wad=zero; wbd=zero ! <- extreme left
elseif(rx>=one)then; wa=zero; wb=one; wad=zero; wbd=zero ! <- extreme right
else                                                     ! Interior:
   b=rx*(one-rx)
   select case(shape)
   case(0); call wbend(rx,wb,wbd)
   case(1); wb=rx;                                             wbd=one
   case(2); wb=rx**2*(3-2*rx);                                 wbd=6*b
   case(3); wb=rx**3*(10+rx*(-15+rx*6));                       wbd=30*b**2
   case(4); wb=rx**4*(35+rx*(-84+rx*(70-rx*20)));              wbd=140*b**3
   case(5); wb=rx**5*(126+rx*(-420+rx*(540+rx*(-315+rx*70)))); wbd=630*b**4
   case(6); wb=rx**6*(462+rx*(-1980+rx*(3465+rx*(-3080+rx*(1386-rx*252)))))
      wbd=2772*b**5
   case default; stop 'In getgbintd; shape code not recognized'
   end select
   wbd=wbd/(xss(ib)-xss(ia)); wa=one-wb; wad=-wbd
endif

! Obtain the pair of Lagrange weight stencils, wta and wtb, centered on
! grid points, ia and ib:
npt=npth*2; nptm=npt-1
call getint(nptm,qss(:,ia),xss(Lw:Mw-1),x,wta,wtad)
call getint(nptm,qss(:,ib),xss(Lw+1:Mw),x,wtb,wtbd)

! Combine the Lagrange weight stencils with blending weights, wa and wb,
! to obtain the final weight stencil, wt:
wt=zero; wtd=zero; wt(1:nptm)=wa*wta; wt(2:npt )=wt(2:npt)+wb*wtb
wtd(1:nptm)=            wad*wta+wa*wtad
wtd(2:npt )=wtd(2:npt )+wbd*wtb+wb*wtbd
end subroutine getgbintd

!=============================================================================
subroutine getlint(npth,rx, wt)!                                      [getint]
!=============================================================================
! Unblended Lagrange polynomial interpolation weights to a target located a 
! fraction rx along the central interval of a set of 2*npth uniformly-
! spaced source grid points here indexed [1-npth,npth]. The weights are
! returned in array wt.
!=============================================================================
integer,                         intent(IN ):: npth
real(dp),                        intent(IN ):: rx
real(dp),dimension(1-npth:npth), intent(OUT):: wt
integer:: npthm
!=============================================================================
if(npth*2-1/=nxs)call inibint(npth*2-1)
npthm=npth-1
call getlap(nxsp,qq,xs, rx+npthm,wt)
end subroutine getlint

!=============================================================================
subroutine getlintd(npth,rx, wt,wtd)!                                 [getint]
!=============================================================================
! Like getlint, but also provide the weights wtd for the derivative of the
! target, assuming the source grid to have unit spacing.
!=============================================================================
integer,                         intent(IN ):: npth
real(dp),                        intent(IN ):: rx
real(dp),dimension(1-npth:npth), intent(OUT):: wt,wtd
integer:: npthm
!=============================================================================
if(npth*2-1/=nxs)call inibint(npth*2-1)
npthm=npth-1
call getlapd(nxsp,qq,xs, rx+npthm,wt,wtd)
end subroutine getlintd

!==============================================================================
subroutine getlap(n,qs,xss,x,w)!                                       [getint]
!==============================================================================
! Assuming that the Q-factors, qs,  have already been calculated, get the
! N general Lagrange interpolation weights, w, for the local span of the grid
! given by the nodes, XS, for an interpolation to the target, X.
!==============================================================================
integer,              intent(IN ):: n
real(dp),dimension(n),intent(IN ):: qs,xss
real(dp),             intent(IN ):: x
real(dp),dimension(n),intent(OUT):: w
!------------------------------------------------------------------------------
real(dp),dimension(n):: r
integer              :: i
real(dp)             :: L
!==============================================================================
L=one
do i=1,n
   r(i)=x-xss(i)
   L=L*r(i)
   if(r(i)/=zero)r(i)=one/r(i)
enddo
w=one
do i=1,n
   if(r(i)==zero)cycle
   w(i)=qs(i)*L*r(i)
enddo
end subroutine getlap

!==============================================================================
subroutine getlapd(n,qs,xss,x,w,wd)!                                   [getint]
!==============================================================================
! Like getlap, but also provide the weights, wd, for the derivative of the
! target point.
!==============================================================================
integer,              intent(IN ):: n
real(dp),dimension(n),intent(IN ):: qs,xss
real(dp),             intent(IN ):: x
real(dp),dimension(n),intent(OUT):: w,wd
!------------------------------------------------------------------------------
real(dp),dimension(n):: r
integer              :: i,j
real(dp)             :: L,Ld
!==============================================================================
L =one
Ld=zero
do i=1,n
   r(i)=x-xss(i)
   Ld=Ld*r(i)+L
   L =L *r(i)
   if(r(i)/=zero)then
      r(i)=one/r(i)
   endif
enddo
w =one
wd=zero
j=0
do i=1,n
   if(r(i)==zero)then
      j=i
      cycle
   endif
   w(i)=qs(i)*L*r(i)
   wd(i)=qs(i)*r(i)*(Ld-L*r(i))
enddo
if(j/=0)wd(j)=-sum(wd)
end subroutine getlapd

!==============================================================================
subroutine getq(n,xss,qs)!                                               [getq]
!==============================================================================
! Get the Lagrange interpolation Q-factors, qs, for a grid span that need not
! be regular. The span of the grid corresponding to the stencil is given
! in the array xss.
!==============================================================================
integer,              intent(IN ):: n
real(dp),dimension(n),intent(IN ):: xss
real(dp),dimension(n),intent(OUT):: qs
!-----------------------------------------------------------------------------
integer:: i,j
!==============================================================================
do j=1,n
   qs(j)=1
   do i=1,n
      if(i==j)cycle
      qs(j)=qs(j)*(xss(j)-xss(i))
   enddo
   qs(j)=1/qs(j)
enddo
end subroutine getq
!==============================================================================
subroutine getgq2(L,M,npth,xss,qss)!                                     [getq]
!==============================================================================
integer,                                       intent(IN ):: L,M,npth
real(dp),dimension(L:M),                       intent(IN ):: xss
real(dp),dimension(npth*2-1,L+npth-1:M+1-npth),intent(OUT):: qss
!=============================================================================
call getq(1+M-L,npth,xss,qss)
end subroutine getgq2
!==============================================================================
subroutine getgq1(m,npth,xss,qss)!                                       [getq]
!==============================================================================
integer,                                   intent(IN ):: m,npth
real(dp),dimension(m),                     intent(IN ):: xss
real(dp),dimension(npth*2-1,npth:m+1-npth),intent(OUT):: qss
!-----------------------------------------------------------------------------
integer:: j,npthm,nptm
!=============================================================================
npthm=npth-1; nptm=npth*2-1
do j=1+npthm,m-npthm; call getq(nptm,xss(j-npthm:j+npthm),qss(:,j)); enddo
end subroutine getgq1

!=============================================================================
subroutine getgsplqy(L,m,nh,xss, qss,yss)!                              [getq]
!=============================================================================
! For spline-weighted interpolation, get both the q-factors (1/denominators) 
! for component Lagrange polynomials, and the Taylor coefficients yss, for the
! splines that, collectively, partition unity in a smooth way.
! The following conventions are adopted:
! Grid points are indexed between L and m;
! Intervals are indexed between 1 and m;
! q-sets for each polynomial are each nominally associated with the grid point
! at which the stencil begins;
! yss-sets are each associated with the grid point at which this sigmoid
! begins. 
! Lagrange polynomials are all of degree nh-1.
! At the outer two intervals, just one polynomial components is combined,
! and, if the grid is sufficiently large, the innermost intervals involve
! combinations of nh Lagrange polynomials.
!=============================================================================
integer,                              intent(in ):: L,m,nh
real(dp),dimension(L:m),              intent(in ):: xss
real(dp),dimension(0:nh,L:m),         intent(out):: qss
real(dp),dimension(0:nh-1,0:nh-1,L:m),intent(out):: yss
!-----------------------------------------------------------------------------
integer                        :: i,j,k,nhm,nhp
real(dp),dimension(0:nh-1,nh-1):: yy
!=============================================================================
nhm=nh-1
nhp=nh+1
qss=zero
do k=L,m-nh
   call getq(nhp,xss(k:k+nh),qss(:,k))
enddo

yss=zero
do i=L+1,m-nh
   call bsigmoid(nhm,xss(i:i+nhm),yy)
   do k=1,nhm
      yss(:,k,i)=yy(:,k)
   enddo
enddo
end subroutine getgsplqy

!=============================================================================
subroutine btableau(n,tableau)!                                     [btableau]
!=============================================================================
! Construct the integer tableau of values proportional to the polynomial
! coefficients of the unit-grid B-splines and their derivatives. The constant 
! of proportionality is n!. 
! The derivative coefficients are progressively offset. The origin of each
! polynomial is the "left" of its respective interval, assuming the axis
! points to the right.
! Collectively, the spline segments form a bell shape spanning intervals
! 0 to n, and the integral of their graph is equal to one. Also, for any
! given 0<r<1, the sum of the values of the spline segements at position
! r in their respective intervals is unity, regardless of r, so their
! values constitute a convenient smooth partition of unity as a function of r.
!=============================================================================
use pcomb, only: factorial,binomial
integer,                       intent(in ):: n
integer,dimension(0:n,0:n,0:n),intent(out):: tableau
!-----------------------------------------------------------------------------
integer,dimension(0:n):: vec,facs,bincos
integer               :: alt,i,j,jm,k,np,facn
!=============================================================================
tableau=0
np=n+1
alt=1
do i=0,n
   facs(i)=factorial(i)
   bincos(i)=alt*binomial(np,i) ! Alternating-sign binomial coefficients
   alt=-alt
enddo
facn=factorial(n)
bincos=facn*bincos ! Weight bincos by n!
vec=0
do k=0,n !<- loop over the intervals spanned by this B-spline
   vec(n)=vec(n)+bincos(k)
   tableau(:,0,k)=vec/facs
   vec(0)=sum(tableau(:,0,k))
   do j=1,n
      jm=j-1
      do i=1,n-jm
         tableau(i+jm,j,k)=i*tableau(i+jm,jm,k)
      enddo
      vec(j)=sum(tableau(:,j,k))
   enddo
enddo
end subroutine btableau

!=============================================================================
subroutine bsigmoid(n,x,yy)!                                        [bsigmoid]
!=============================================================================
! Given a non-uniform grid x(0:n) of n intervals, find the Taylor coefficients
! at the start of each interval of the n-th degree polynomial within that
! interval so the, collectively, these polynomial pieces form a sigmoid
! function continuous at interval boundaries until the n-th derivative,
! joining y=0 at <x(0) with y=1 at >x(n). The polynomials are related to
! the B-splines of degree n.
!=============================================================================
use pcomb, only: flag,next_bin,factorial
integer,                  intent(in ):: n
real(dp),dimension(0:n),  intent(in ):: x
real(dp),dimension(0:n,n),intent(out):: yy
!-----------------------------------------------------------------------------
real(dp)                   :: p,rr
real(dp),dimension(0:n-1,n):: xx
integer                    :: i,j,k
integer,dimension(n-1)     :: b
integer,dimension(n)       :: alt
!=============================================================================
alt(1)=1; do i=2,n; alt(i)=-alt(i-1); enddo
do j=1,n; do i=0,j-1; xx(i,j)=x(j)-x(i); enddo; enddo ! (Upper triangular)
yy=zero
if(n==1)then
   yy(1,1)=one/xx(0,1)
else
   flag=T
   do
      call next_bin(b)
      if(flag)then; flag=F; exit; endif
      k=1+sum(b)
      call binprod(n,b,xx,p)
      yy(n,k)=yy(n,k)+alt(k)/p
   enddo
endif
yy(n,:)=yy(n,:)*factorial(n)
! Coefficients yy(:,1) are already complete; now find coefficients in the
! remaining intervals:
do k=2,n
   j=k-1
   rr=xx(j-1,j)! <- width of preceding interval
   do i=0,n-1
      call evaltay(n-i,yy(i:n,j),rr,yy(i,k))
   enddo
enddo
end subroutine bsigmoid

!=============================================================================
subroutine binprod(n,b,xx,p)!                                        [binprod]
!=============================================================================
! Construct the binary-coded product term, p, from the pool of potential
! factors supplied by xx.
!=============================================================================
integer,                    intent(in ):: n
integer ,dimension(n-1),    intent(in ):: b
real(dp),dimension(0:n-1,n),intent(in ):: xx
real(dp),                   intent(out):: p
!-----------------------------------------------------------------------------
integer:: i,j,k
!=============================================================================
j=sum(b)
k=j+1
p=xx(j,k)
do i=1,n-1
   if(b(i)==0)then; k=k+1
   else           ; j=j-1
   endif
   p=p*xx(j,k)
enddo
end subroutine binprod
   
!=============================================================================
subroutine evaltay(n,co,x,y)!                                        [evaltay]
!=============================================================================
! Evaluate a Taylor-Maclaurin expansion of degree n.
! co(0:n) are the Taylor coefficients.
!=============================================================================
integer,                intent(in ):: n
real(dp),dimension(0:n),intent(in ):: co
real(dp),               intent(in ):: x
real(dp),               intent(out):: y
!----------------------------------------------------------------------------
integer:: i
!============================================================================
y=zero
do i=n,1,-1
   y=(co(i)+y)*x/i
enddo
y=y+co(0)
end subroutine evaltay

!=============================================================================
subroutine evaltayd(n,co,x,y,yd)!                                    [evaltay]
!=============================================================================
! Evaluate a Taylor-Maclaurin expansion of degree n and its derivative.
! co(0:n) are the Taylor coefficients.
!=============================================================================
integer,                intent(in ):: n
real(dp),dimension(0:n),intent(in ):: co
real(dp),               intent(in ):: x
real(dp),               intent(out):: y,yd
!----------------------------------------------------------------------------
integer:: i
!============================================================================
yd=zero
y =zero
do i=n,1,-1
   yd=(yd*x+co(i)+y)/i
   y =(co(i)+y)*x/i
enddo
y=y+co(0)
end subroutine evaltayd

end module pbint