C
C $Id: bspline.F,v 1.2 1998/07/16 16:39:28 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
C     This file contains all the Cardinal B-spline functions needed
C     by the PME routines:
C
C     bspline      single point B-Spline algorithm
C     mbspline1    massive B-Spline algorithm, returns splines
C     mbspline     massive B-Spline algorithm, returns splines + derivatives
C     mkbspline    massive B-Spline algorithm for integers

CC----------------------------------------------------------------CC
CC----------------------------------------------------------------CC

#if 0
C     Note that the following routine is NOT used in the PME
C     algorithm. Instead of a multiple calls to bspline, massive
C     B-spline algorithms are used that make the computation
C     significantly faster (significant decrease of actual computation
C     time, and decrease of system-time by eliminating multiple calls)
C
C     Nevertheless, this algorithm is very instructive and
C     served as a model for the massive B-Spline algorithms.
C

      subroutine bsplinebis(u, n, smn, dsmn)

C     returns the Cardinal B-spline Mn(u) and
C     the derivative of the Cardinal B-spline at point u
C     (d Mn(u)/du)
C
C     Arjan van der Vaart, Oct. 1997
C
C     parameters:
C     *) u is the scaled fractional coordinate defined by
C     .  u(i) = K(i)*a(i)*r(i), i=1,3
C     .    with K integer representing the grid size
C     .         a the reciprocal box vector
C     .         r the (Cartesian) coordinate of the particle
C     *) n is the order of the B-spline, n > 2, n must be even
C     *) smn(1) will contain Mn(u), the dimension of smn
C     .  should be >= n-1 [used as workspace]
C     *) dsmn will contain dMn(u)

      implicit double precision(a-h,o-z)
      parameter (one=1.0)
      dimension smn(*)

CC----------------------------------------------------------------CC

      nm1 = n-1
      np1 = n+1
      xn = dble(n)
      xnm1 = xn-1.0

C     fill second level = M2(u), M2(u-1), M2(u-2), .. , M2(u-n+2)
      do 10 i=1,nm1
         x = abs(u-i)
         smn(i) = 1.0-min(one,x)
 10   continue

C     fill third till n-1 level = Mi(u), Mi(u-1) etc.
      do 30 i=3,nm1
         xi = dble(i)
         xim1 = xi-1.0
         je = np1-i
         do 20 j=1,je
            w = u-dble(j-1)
            smn(j) = (w*smn(j) + (xi-w)*smn(j+1))/xim1
 20      continue
 30   continue

C     calculate the derivative dMn(u)/du
      dsmn = smn(1)-smn(2)

C     calculate Mn(u)
      smn(1) = (u*smn(1)+(xn-u)*smn(2))/xnm1
      end

CC----------------------------------------------------------------CC
CC----------------------------------------------------------------CC

#endif

      subroutine mbspline1(upmex,upmey,upmez,xmn,ymn,zmn)

C     massive bspline algorithm
C     returns the Cardinal B-spline Mn(w) at point w
C     for w=u+int(n-u), u+int(n-u-1), .. , u, u-1,
C     .                  u-2, ... , u-n
C     in the x, y and z-direction.
C     [note that n=nspline is the order of the B-Spline]
C
C     Arjan van der Vaart, Oct. 1997
C
C     parameters:
C     *) iat is the atom number. This is needed for the scaled fractional
C     .  coordinate defined by
C     .  u(i) = K(i)*a(i)*r(i), i=1,3
C     .    with K integer representing the grid size
C     .         a the reciprocal box vector
C     .         r the (Cartesian) coordinate of the particle
C     *) xmn will contain the Cardinal B-spline of order n in for
C     .  the x coordinate;
C     .    xmn(1) = Mn(u+int(n-u))
C     .    xmn(int(u)+int(n-u)+1) = Mn(u-n)
C     .  xmn should be n+int(u)+int(n-u)+1-2 <= 2n-1 elements long,
C     .  this extra length of smn is used as workspace
C     *) ymn, zmn see xmn
C
C     note that nspline (= n) is the order of the B-spline,
C     nspline > 2, nspline must be even

      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

      dimension xmn(*), ymn(*), zmn(*)

CC----------------------------------------------------------------CC

C     x coordinate
      iplusx = nspline-upmex
      iminx = upmex
      ikx = iminx+iplusx
      mx = ikx+1
      wx = upmex+iplusx

C     y coordinate
      iplusy = nspline-upmey
      iminy = upmey
      iky = iminy+iplusy
      my = iky+1
      wy = upmey+iplusy

C     z coordinate
      iplusz = nspline-upmez
      iminz = upmez
      ikz = iminz+iplusz
      mz = ikz+1
      wz = upmez+iplusz

      m = max(mx,my,mz)

C     fill M2
      do 10 i=1,nspline+m-2
         xmn(i) = 0.0
         ymn(i) = 0.0
         zmn(i) = 0.0
 10   continue
      xmn(ikx+1) = upmex-iminx
      xmn(ikx) = 1.0-xmn(ikx+1)
      ymn(iky+1) = upmey-iminy
      ymn(iky) = 1.0-ymn(iky+1)
      zmn(ikz+1) = upmez-iminz
      zmn(ikz) = 1.0-zmn(ikz+1)

C     fill M3 through M(nspline)
      do 50 i=3,nspline
         di = dble(i)
         dim1 = di-1.0
         do 20 j=1,nspline-i+mx
            x = wx-(j-1)
            xmn(j) = (x*xmn(j) + (di-x)*xmn(j+1))/dim1
 20      continue
         do 30 j=1,nspline-i+my
            y = wy-(j-1)
            ymn(j) = (y*ymn(j) + (di-y)*ymn(j+1))/dim1
 30      continue
         do 40 j=1,nspline-i+mz
            z = wz-(j-1)
            zmn(j) = (z*zmn(j) + (di-z)*zmn(j+1))/dim1
 40      continue
 50   continue

      end


CC----------------------------------------------------------------CC
CC----------------------------------------------------------------CC


      subroutine mbspline(upmex,upmey,upmez,xmn,ymn,zmn,dxmn,dymn,dzmn)

C     massive bspline algorithm
C     returns the derivative of the Cardinal B-spline at point w
C     (d Mn(w)/dw) for w=u+int(n-u), u+int(n-u-1), .. , u, u-1,
C     .                  u-2, ... , u-n
C     in the x, y and z-direction.
C     [note that n=nspline is the order of the B-Spline]
C
C     Arjan van der Vaart, Oct. 1997
C
C     parameters:
C     *) iat is the atom number. This is needed for the scaled fractional
C     .  coordinate defined by
C     .  u(i) = K(i)*a(i)*r(i), i=1,3
C     .    with K integer representing the grid size
C     .         a the reciprocal box vector
C     .         r the (Cartesian) coordinate of the particle
C     *) xmn will contain the Cardinal B-spline of order n in for
C     .  the x coordinate;
C     .    xmn(1) = Mn(u+int(n-u))
C     .    xmn(int(u)+int(n-u)+1) = Mn(u-n)
C     .  xmn should be n+int(u)+int(n-u)+1-2 <= 2n-1 elements long,
C     .  this extra length of smn is used as workspace
C     *) ymn, zmn see xmn
C     *) dxmn contains the derivatives dMn(w)/dw in the x direction;
C     .    dxmn(1) = dMn(u+int(n-u))/d(u+int(n-u))
C     .    dxmn(int(u)+int(n-u)+1) = dMn(u-n)/d(u-n)
C     .  dxmn should be int(u)+int(n-u)+1 <= n+1 elements long.
C     .  analogous for dymn and dzmn
C
C     note that nspline (= n) is the order of the B-spline,
C     nspline > 2, nspline must be even

      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     parameters:
      dimension xmn(*), ymn(*), zmn(*),
     &     dxmn(*), dymn(*), dzmn(*)

CC----------------------------------------------------------------CC

C     x coordinate
      iplusx = nspline-upmex
      iminx = upmex
      ikx = iminx+iplusx
      mx = ikx+1
      wx = upmex+iplusx

C     y coordinate
      iplusy = nspline-upmey
      iminy = upmey
      iky = iminy+iplusy
      my = iky+1
      wy = upmey+iplusy

C     z coordinate
      iplusz = nspline-upmez
      iminz = upmez
      ikz = iminz+iplusz
      mz = ikz+1
      wz = upmez+iplusz

      m = max(mx,my,mz)

C     fill M2
      do 10 i=1,nspline+m-2
         xmn(i) = 0.0
         ymn(i) = 0.0
         zmn(i) = 0.0
 10   continue
      xmn(ikx+1) = upmex-iminx
      xmn(ikx) = 1.0-xmn(ikx+1)
      ymn(iky+1) = upmey-iminy
      ymn(iky) = 1.0-ymn(iky+1)
      zmn(ikz+1) = upmez-iminz
      zmn(ikz) = 1.0-zmn(ikz+1)

C     fill M3 through M(nspline-1)
      do 50 i=3,nspline-1
         di = dble(i)
         dim1 = di-1.0
         do 20 j=1,nspline-i+mx
            x = wx-(j-1)
            xmn(j) = (x*xmn(j) + (di-x)*xmn(j+1))/dim1
 20      continue
         do 30 j=1,nspline-i+my
            y = wy-(j-1)
            ymn(j) = (y*ymn(j) + (di-y)*ymn(j+1))/dim1
 30      continue
         do 40 j=1,nspline-i+mz
            z = wz-(j-1)
            zmn(j) = (z*zmn(j) + (di-z)*zmn(j+1))/dim1
 40      continue
 50   continue

C     fill dMn and Mn
      dnm1 = dnspline-1.0
      xnmw = dnspline-wx
      do 60 i=1,mx
         im1 = i-1
         dxmn(i) = xmn(i)-xmn(i+1)
         xmn(i) = ((wx-im1)*xmn(i)+(xnmw+im1)*xmn(i+1))/dnm1
 60   continue
      ynmw = dnspline-wy
      do 70 i=1,my
         im1 = i-1
         dymn(i) = ymn(i)-ymn(i+1)
         ymn(i) = ((wy-im1)*ymn(i)+(ynmw+im1)*ymn(i+1))/dnm1
 70   continue
      znmw = dnspline-wz
      do 80 i=1,mz
         im1 = i-1
         dzmn(i) = zmn(i)-zmn(i+1)
         zmn(i) = ((wz-im1)*zmn(i)+(znmw+im1)*zmn(i+1))/dnm1
 80   continue

      end

CC----------------------------------------------------------------CC
CC----------------------------------------------------------------CC

      subroutine mkbspline(smn)

C     massive bspline algorithm
C     returns the Cardinal B-spline Mn(w) at point k+1 for
C     k=0, 1, .. , n-2
C
C     Arjan van der Vaart, Oct. 1997
C
C     parameters:
C     *) smn will contain the Cardinal B-spline;
C     .    smn(1) = Mn(1)
C     .    smn(k) = Mn(k+1)
C
C     nspline is the order of the B-Spline
C
      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

      dimension smn(*), wmn(maxspline2m3)

CC----------------------------------------------------------------CC

      m = nspline-1
      n2m1 = m+nspline
      n2m3 = n2m1-2

C     fill Mn2
      do 10 i=1,n2m3
         wmn(i) = 0.0
 10   continue
      wmn(m) = 1.0

C     fill Mn3 through Mn(nspline-1)
      do 30 i=3,m
         xi = dble(i)
         xim1 = xi-1.0
         do 20 j=1,n2m1-i
            x = dnspline-j
            wmn(j) = (x*wmn(j) + (xi-x)*wmn(j+1))/xim1
 20      continue
 30   continue

C     fill Mn(nspline)
      xim1 = dnspline-1.0
      do 40 j=1,m
         nj = nspline-j
         x = dble(nj)
         smn(nj) = (x*wmn(j) + (dnspline-x)*wmn(j+1))/xim1
 40   continue

      end