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/spline.f | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 src/labat/spline.f (limited to 'src/labat/spline.f') 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 -- cgit v1.2.3