1 files changed, 65 insertions, 0 deletions

diff --git a/src/labat/deriv.f b/src/labat/deriv.f
new file mode 100644
index 0000000..b440070
--- /dev/null
+++ b/src/labat/deriv.f
@@ -0,0 +1,65 @@
+      function deriv(r,x,x2,spin,i,imin,imax)
+c  deriv calculates the derivative of x with respect to r, at
+c  position r(i).
+c  input r        : radial coordinate (array)
+c  input x        : function to be derivated  (array)
+c  input x2       : second derivative of x (array)
+c  input spin     : -1 if spin-down, 1 if spin-up, 0 otherwise
+c  input i        : counter, i.e. current position = r(i)
+c  input imin     : min value of i
+c  input imax     : max value of i
+      implicit logical (a-z)
+      double precision deriv,r,x,x2
+      double precision h,d,dforw,dback,dforup,dfordn,dbacup,dbacdn
+      integer spin,i,imin,imax
+      dimension r(561),x(2,561),x2(2,561)
+      if (i.eq.imin) then
+         h = (r(i+1)-r(i))/2.d0
+         if (spin.eq.-1) then
+            d = x(2,i)
+            call splint(r,x,x2,-1,imin,imax,r(i)+h,dforw)
+         elseif (spin.eq.1) then
+            d = x(1,i)
+            call splint(r,x,x2,1,imin,imax,r(i)+h,dforw)
+         else
+            d = (x(1,i) + x(2,i))
+            call splint(r,x,x2,1,imin,imax,r(i)+h,dforup)
+            call splint(r,x,x2,-1,imin,imax,r(i)+h,dfordn)
+            dforw = dforup + dfordn
+         endif   
+         deriv = (dforw-d)/h
+      elseif (i.eq.imax) then
+         h = (r(i)-r(i-1))/2.d0
+         if (spin.eq.-1) then
+            call splint(r,x,x2,-1,imin,imax,r(i)-h,dback)
+            d = x(2,i)
+         elseif (spin.eq.1) then
+            call splint(r,x,x2,1,imin,imax,r(i)-h,dback)
+            d = x(1,i)
+         else
+            call splint(r,x,x2,1,imin,imax,r(i)-h,dbacup)
+            call splint(r,x,x2,-1,imin,imax,r(i)-h,dbacdn)
+            dback = dbacup + dbacdn
+            d = x(1,i) + x(2,i)
+         endif   
+         deriv = (d-dback)/h
+      else
+         h = (r(i)-r(i-1))/2.d0
+         if (spin.eq.-1) then
+            call splint(r,x,x2,-1,imin,imax,r(i)+h,dforw)
+            call splint(r,x,x2,-1,imin,imax,r(i)-h,dback)
+         elseif (spin.eq.1) then
+            call splint(r,x,x2,1,imin,imax,r(i)+h,dforw)
+            call splint(r,x,x2,1,imin,imax,r(i)-h,dback)
+         else
+            call splint(r,x,x2,-1,imin,imax,r(i)+h,dfordn)
+            call splint(r,x,x2,-1,imin,imax,r(i)-h,dbacdn)
+            call splint(r,x,x2,1,imin,imax,r(i)+h,dforup)
+            call splint(r,x,x2,1,imin,imax,r(i)-h,dbacup)
+            dback = dbacup + dbacdn
+            dforw = dforup + dfordn
+         endif   
+         deriv = (dforw-dback)/(2.d0*h)
+      endif
+      return
+      end