subroutine splint(r,db,db2,spin,start,n,ri,dbi)
c  splint calculates the cubic spline interpolation of the density
c  together with subroutine spline.
c  the main ideas are "stolen" from "numerical recipes", cambridge
c  university press (1992).
c  input r     : position coordinate (array)
c  input db    : spin up and down density. (array)
c  input db2   : second derivative of db (array)
c  input spin  : 1 if spin up, -1 if spin down
c  input start : starting point of spline = r(start)
c  input n     : index, r(n) = r_max
c  input ri    : position of interpolation
c  output dbi  : interpolated value of db
      implicit logical (a-z)
      integer spin,spinn,n,klo,khi,k,start
      double precision r,db,db2,ri,dbi,h,a,b
      dimension r(561),db(2,561),db2(2,561)
      spinn = spin
      if (spinn.eq.-1) spinn = 2
      klo = start
      khi = n
 10   if (khi-klo.gt.1) then
         k = (khi+klo)/2
         if (r(k).gt.ri) then
            khi = k
         else
            klo = k
         endif
      goto 10
      endif
      h = r(khi)-r(klo)
      a = (r(khi)-ri)/h
      b = (ri-r(klo))/h
      dbi = a*db(spinn,klo)+b*db(spinn,khi)+((a**3-a)*
     j     db2(spinn,klo)+(b**3-b)*db2(spinn,khi))*(h**2)/6.d0
      return
      end