1 files changed, 41 insertions, 0 deletions

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