From 5df79c53745fde5d6c3340a2979b1429cd5892c1 Mon Sep 17 00:00:00 2001 From: Henrik Rydberg Date: Sat, 8 Oct 2011 20:30:28 +0200 Subject: Initial import of htcd system 1.0 Signed-off-by: Henrik Rydberg --- src/labat/splint.f | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 src/labat/splint.f (limited to 'src/labat/splint.f') 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 -- cgit v1.2.3