1 files changed, 37 insertions, 0 deletions

diff --git a/src/labat/splint.f b/src/labat/splint.f
new file mode 100644
index 0000000..a25a942
--- /dev/null
+++ b/src/labat/splint.f
@@ -0,0 +1,37 @@
+      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