1 files changed, 34 insertions, 0 deletions

diff --git a/src/labat/spline.f b/src/labat/spline.f
new file mode 100644
index 0000000..aa523be
--- /dev/null
+++ b/src/labat/spline.f
@@ -0,0 +1,34 @@
+      subroutine spline(r,start,n,db,db2)
+c  spline calculates the cubic spline interpolation of the density
+c  together with subroutine splint.
+c  the main ideas are "stolen" from "numerical recipes", cambridge
+c  university press (1992).
+c  input r     : position coordinate (array)
+c  input start : index, r(start) = starting point of the spline
+c  input n     : index, r(n) = r_max
+c  input db    : spin up and down density. (array)
+c  output db2  : second derivative of db (array)
+      implicit logical (a-z)
+      double precision r,db,db2,u,sig,p,qn,un
+      integer spin,start,n,i
+      dimension r(561),db(2,561),db2(2,561),u(561)
+      do 100 spin = 1,2
+         db2(spin,start) = 0.d0
+         u(start) = 0.d0
+         do 10 i = start+1,n-1
+            sig = (r(i)-r(i-1))/(r(i+1)-r(i-1))
+            p = sig*db2(spin,i-1)+2.d0
+            db2(spin,i) = (sig-1.d0)/p
+            u(i) = (6.d0*((db(spin,i+1)-db(spin,i))/(r(i+1)-r(i))-
+     j           (db(spin,i)-db(spin,i-1))/(r(i)-r(i-1)))/
+     j           (r(i+1)-r(i-1))-sig*u(i-1))/p
+ 10      continue
+         qn = 0.d0
+         un = 0.d0
+         db2(spin,n) = (un-qn*u(n-1))/(qn*db2(spin,n-1)+1.d0)
+         do 20 i = n-1,start,-1
+            db2(spin,i) = db2(spin,i)*db2(spin,i+1)+u(i)
+ 20      continue
+ 100  continue
+      return
+      end