1c
2c $Id: hnd_dervxyz.F 19696 2010-10-29 16:53:42Z d3y133 $
3c
4c Taken and modified from HONDO
5c
6      subroutine dim_dervxyz(mder)
7      implicit none
8c
9c     Routine can do ddvxyz and dvxyz by setting mder to 2 or 1
10c
11c     ----- gauss-hermite quadrature using minimum point formula -----
12c
13#include "hnd_whermt.fh"
14c
15      double precision xint, yint, zint, t, x0, y0, z0, xi, yi, zi
16      double precision xj, yj, zj, cx, cy, cz, zero, ptx, pty, ptz
17      double precision dum, px, py, pz, ax, ay, az, bx, by, bz
18      integer ni, nj, npts, imin, imax, i, mder, ii, jj
19      common/hnd_xyzder/xint,yint,zint,t,x0,y0,z0,xi,yi,zi,xj,yj,zj
20     1                  ,ni,nj,cx,cy,cz
21      data zero /0.0d+00/
22      xint = zero
23      yint = zero
24      zint = zero
25      npts = (ni+nj+mder-2)/2+1
26      imin = hermin(npts)
27      imax = hermax(npts)
28      do 11 i = imin,imax
29         dum = h(i)*t
30         ptx = dum+x0
31         pty = dum+y0
32         ptz = dum+z0
33         px = ptx-cx
34         py = pty-cy
35         pz = ptz-cz
36         px = px**mder
37         py = py**mder
38         pz = pz**mder
39         ax = ptx-xi
40         ay = pty-yi
41         az = ptz-zi
42         bx = ptx-xj
43         by = pty-yj
44         bz = ptz-zj
45         do ii = 1, ni-1
46            px=px*ax
47            py=py*ay
48            pz=pz*az
49         enddo
50         do jj = 1, nj-1
51            px=px*bx
52            py=py*by
53            pz=pz*bz
54         enddo
55         dum = w(i)
56         xint = xint+dum*px
57         yint = yint+dum*py
58         zint = zint+dum*pz
59   11 continue
60      return
61      end
62