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/laplac.f | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 src/labat/laplac.f (limited to 'src/labat/laplac.f') diff --git a/src/labat/laplac.f b/src/labat/laplac.f new file mode 100644 index 0000000..9a7dfd9 --- /dev/null +++ b/src/labat/laplac.f @@ -0,0 +1,41 @@ + subroutine laplac(r,gradup,graddn,i,imin,imax,done, + j lapl,laplup,lapldn) +c calculates laplace(density), i.e. (del^2)(density) +c spherical coordinates => (del^2)phi = 1/r^2*d/dr[r^2*d(phi)/dr] +c input r : radial coordinate (array) +c input gradup : grad(density), spin up (array) +c input graddn : grad(density), spin down (array) +c input i : counter, i.e. current position = r(i) +c input imin : min value of i +c input imax : max value of i +c in/output done : true if r^2*grad(density) is splined +c output lapl : laplace(density) +c output laplup,lapldn : laplace(density), spin up and down + implicit logical (a-z) + logical done + double precision r,gradup,graddn,lapl,laplup,lapldn,lapl2,temp + double precision rj2,r2,deriv + integer i,imin,imax,j + common/ggalap/temp,lapl2 + dimension r(561),gradup(561),graddn(561) + dimension lapl2(2,561),temp(2,561) + if (.not.done) then + do 10 j = imin,imax + rj2 = r(j)**2 + temp(1,j) = rj2*gradup(j) + temp(2,j) = rj2*graddn(j) + 10 continue + call spline(r,imin+1,imax,temp,lapl2) + done = .true. + endif + if (i.eq.imin) then + call laplac2(r,gradup,graddn,i+1,imin,imax,done,lapl,laplup, + j lapldn) + else + r2 = r(i)**2 + laplup = (deriv(r,temp,lapl2,1,i,imin,imax))/r2 + lapldn = (deriv(r,temp,lapl2,-1,i,imin,imax))/r2 + lapl = laplup + lapldn + endif + return + end -- cgit v1.2.3