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