c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c
c     Program pwrvtx
c
c     This program will compute the Power diagram vertices and the
c     unbounded Power diagram edges for a set of points in 3-d
c     space from a Regular tetrahedralization for the set.
c
c     A Power diagram is essentially a Voronoi diagram with weights.
c
c     If no weights are present the program will compute the Voronoi
c     diagram vertices and the unbounded Voronoi diagram edges for
c     the set of points from a Delaunay tetrahedralization for the set.
c
c     It is assumed that a Regular/Delaunay tetrahedralization for
c     the set of points has already been computed using existing
c     program regtet/ragtat and that the output tetrahedron list
c     produced by regtet/ragtat contains only real tetrahedra (no
c     artificial ones), i. e. logical variable artfcl was set to
c     .false. during the execution of regtet/ragtat. The output
c     tetrahedron list from regtet/ragtat is then used as input below.
c     It is noted that because of degeneracies more than one
c     tetrahedron may correspond to the same Power/Voronoi vertex.
c     In addition, because of degeneracies, Power/Voronoi cells of
c     points may have facets that are not 2-dimensional, and if
c     weights are present there may be Power cells of points that are
c     not 3-dimensional.
c
c     Subroutine pwrvtx is the driver routine of the program and is
c     called from the main routine.
c
c     Description of output variables appears in driver routine pwrvtx.
c
c     Author: Javier Bernal
c
c     Disclaimer:
c
c     This software was developed at the National Institute of Standards
c     and Technology by employees of the Federal Government in the
c     course of their official duties. Pursuant to title 17 Section 105
c     of the United States Code this software is not subject to
c     copyright protection and is in the public domain. This software is
c     experimental. NIST assumes no responsibility whatsoever for its
c     use by other parties, and makes no guarantees, expressed or
c     implied, about its quality, reliability, or any other
c     characteristic. We would appreciate acknowledgement if the
c     software is used.
c
*MAIN
c
      program main
c
      integer nmax, nvmax
      parameter (nmax=150000, nvmax= 7*nmax)
c
      double precision x(nmax), y(nmax), z(nmax), w(nmax)
      integer ix(nmax), iy(nmax), iz(nmax), iw(nmax)
      integer ix2(nmax), iy2(nmax), iz2(nmax), iw2(nmax)
      double precision xp(nvmax), yp(nvmax), zp(nvmax)
      integer icon(8,nvmax), ifl(nvmax), is(nmax)
      integer nv, nw, nt, nq, naddl, icfig, iwfig
      double precision wlenx, wleny, wlenz, wlenw
      logical delaun, pntoff, flphis, artfcl, random, reccor, redchk
c
      double precision xcor, ycor, zcor, wght
      double precision xpmin, ypmin, zpmin, xpmax, ypmax, zpmax
      integer ideli, ipnti, iflpi, iarti, irani, ireci, iredi
      integer i, j, no
c
      integer ng, ifir, isec, ithi, ifou, larg
      double precision derro, dist1, dist2, dist3, dist4
      double precision dista, disti, diffd, dnom
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'This program is for computing the Power diagram ',
     *          'vertices and the'
      write(*,*)'unbounded Power diagram edges for a set of ',
     *          'points in 3-d space'
      write(*,*)'from a Regular tetrahedralization for the set.'
      write(*,*)' '
      write(*,*)'If no weights are present the Voronoi vertices ',
     *          'and the unbounded'
      write(*,*)'Voronoi edges are computed from a Delaunay ',
     *          'tetrahedralization for'
      write(*,*)'the set.'
      write(*,*)' '
      write(*,*)'It is assumed that a Regular/Delaunay ',
     *          'tetrahedralization for the'
      write(*,*)'set has already been computed using existing program ',
     *          'regtet/ragtat'
      write(*,*)'and that the output tetrahedron list produced by ',
     *          'regtet/ragtat'
      write(*,*)'contains only real tetrahedra (logical variable ',
     *          'artfcl was set to'
      write(*,*)'.false. during the execution of regtet/ragtat).'
      write(*,*)' '
      write(*,*)'The output tetrahedron list from regtet/ragtat ',
     *          'is then used as'
      write(*,*)'input by this program. Note that because of ',
     *          'degeneracies more than'
      write(*,*)'one tetrahedron may correspond to the same ',
     *          'Power/Voronoi vertex.'
      write(*,*)'In addition, because of degeneracies, Power/Voronoi ',
     *          'cells of'
      write(*,*)'points may have facets that are not 2-dimensional, ',
     *          'and if weights'
      write(*,*)'are present there may be Power cells of points that ',
     *          'are not'
      write(*,*)'3-dimensional.'
      write(*,*)' '
      write(*,*)'Arrays xp, yp, zp, icon, ifl will contain output ',
     *          'information.'
      write(*,*)'Description of these arrays appears in driver ',
     *          'routine pwrvtx.'
      write(*,*)' '
      write(*,*)'In the current main program arrays ifl, xp, yp, zp ',
     *          'are saved in'
      write(*,*)'output file power-vrts, and array icon in output ',
     *          'file tetra-unbd.'
c
c----------------------------------------------------------------------
c
c     open files
c
      open (11, file = 'pnts-wghts')
      open (12, file = 'tetrahedra')
      open (13, file = 'power-vrts')
      open (14, file = 'tetra-unbd')
c
c----------------------------------------------------------------------
c
c     read tetrahedralization information
c
      write(*,*)' '
      write(*,*)'Reading tetrahedralization information ...'
      read (12,80) ideli, ipnti, iflpi, iarti, irani, ireci, iredi
      delaun = .false.
      pntoff = .false.
      flphis = .false.
      artfcl = .false.
      random = .false.
      reccor = .false.
      redchk = .false.
      if(ideli.eq.1) delaun = .true.
      if(ipnti.eq.1) pntoff = .true.
      if(iflpi.eq.1) flphis = .true.
      if(iarti.eq.1) artfcl = .true.
      if(irani.eq.1) random = .true.
      if(ireci.eq.1) reccor = .true.
      if(iredi.eq.1) redchk = .true.
      if(.not.delaun) then
         read (12,*) nw, nt, icfig, iwfig
      else
         read (12,*) nw, nt, icfig
      endif
      if(nw.gt.nmax .or. nt.gt.nvmax) stop 10
      read (12,90) ((icon(i,j), i = 1, 8), j = 1, nt)
      read (12,95) (is(i), i = 1, nw)
      if(reccor .and. .not.delaun) then
         read (12,*) wlenx, wleny, wlenz, wlenw
         read (12,*) naddl
      elseif(reccor) then
         read (12,*) wlenx, wleny, wlenz
         read (12,*) naddl
      endif
c
   80 format (7(1x,i1))
   90 format (8i10)
   95 format (7i10)
c
c----------------------------------------------------------------------
c
      if(artfcl)then
         write(*,*)' '
         write(*,*)'Input tetrahedra contain artificial points.'
         write(*,*)'Such points are not allowed in this program.'
         write(*,*)'Program terminated.'
         stop 20
      endif
c
c----------------------------------------------------------------------
c
c     code for avoiding certain messages during compilation
c
      if(pntoff) ipnti = 1
      if(flphis) iflpi = 1
      if(random) irani = 1
      if(redchk) iredi = 1
c
c----------------------------------------------------------------------
c
c     read vertex data
c
      write(*,*)' '
      write(*,*)'Reading vertex data ...'
      nv = 0
      if(delaun) go to 130
c
c     read vertex data with weights
c
  100 continue
      read (11, *, end = 120) xcor, ycor, zcor, wght
      nv = nv + 1
      if (nv .gt. nmax) stop 30
      x(nv) = xcor
      y(nv) = ycor
      z(nv) = zcor
      w(nv) = wght
      go to 100
  120 continue
      go to 140
c
c     read vertex data without weights
c
  130 continue
      read (11, *, end = 140) xcor, ycor, zcor
      nv = nv + 1
      if (nv .gt. nmax) stop 40
      x(nv) = xcor
      y(nv) = ycor
      z(nv) = zcor
      w(nv) = 0.0d0
      go to 130
  140 continue
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'Computation of Power/Voronoi vertices to begin: ',
     *          'driver subroutine'
      write(*,*)'pwrvtx will be called from the main routine and ',
     *          'vertices will be'
      write(*,*)'computed.'
      write(*,*)' '
      write(*,*)'Please wait ...'
      write(*,*)' '
c
c----------------------------------------------------------------------
c 
c     call pwrvtx to compute Power/Voronoi vertices
c
      call pwrvtx(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2,
     *            xp, yp, zp, icon, ifl, is, nv, nw, nt, nq, nmax,
     *            nvmax, wlenx, wleny, wlenz, wlenw, naddl, icfig,
     *            iwfig, delaun, artfcl, reccor)
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'(Back to the main routine).'
      write(*,*)' '
      write(*,*)'Computation of Power/Voronoi vertices has been ',
     *          'completed.'
      write(*,*)' '
      write(*,*)'Number of points in set whose Power/Voronoi ',
     *          'diagram is'
      write(*,*)'under consideration = ',nw
      write(*,*)' '
      write(*,*)'Number of (bounded) tetrahedra in initial tetrahedron'
      write(*,*)'list, i. e. number of Power/Voronoi vertices = ',nt 
      write(*,*)' '
      write(*,*)'Number of bounded and unbounded tetrahedra in ',
     *          'final'
      write(*,*)'tetrahedron list, i. e. number of Power/Voronoi'
      write(*,*)'vertices and unbounded edges = ',nq
c
c----------------------------------------------------------------------
c
c     save Power/Voronoi vertices
c
      write(*,*)' '
      write(*,*)'Saving Power/Voronoi vertices and unbounded edges ',
     *          'information ...'
      write (13,*) nw, nt, nq
      write (13,400) (ifl(i),i=1,nq)
      do 300 i = 1, nq
         write (13,*) xp(i), yp(i), zp(i)
  300 continue
  400 format (40(1x,i1))
c
c----------------------------------------------------------------------
c
c     save tetrahedralization information including unbounded edges
c
      write(*,*)' '
      write(*,*)'Saving final tetrahedralization information ...'
      write (14,80) ideli, ipnti, iflpi, iarti, irani, ireci, iredi
      if(.not.delaun) then
         write (14,*) nw, nt, nq, icfig, iwfig
      else
         write (14,*) nw, nt, nq, icfig
      endif
      write (14,90) ((icon(i,j), i = 1, 8), j = 1, nq)
      write (14,95) (is(i), i = 1, nw)
      if(reccor .and. .not.delaun) then
         write (14,*) wlenx, wleny, wlenz, wlenw
         write (14,*) naddl
      elseif(reccor) then
         write (14,*) wlenx, wleny, wlenz
         write (14,*) naddl
      endif
c
c----------------------------------------------------------------------
c
c     calculating min and max of Power/Voronoi vertices coordinates
c
      do 520 i = 1, nt
         if(ifl(i).eq.0) go to 520
         xpmax = xp(i)
         xpmin = xp(i)
         ypmax = yp(i)
         ypmin = yp(i)
         zpmax = zp(i)
         zpmin = zp(i)
         go to 530
  520 continue
      write(*,*)' '
      write(*,*)'All Power/Voronoi vertices marked as too large.'
      go to 900
  530 continue
      do 550 no = 1, nt
         if(ifl(i).eq.0) go to 550
         if (xp(no) .gt. xpmax) xpmax = xp(no)
         if (xp(no) .lt. xpmin) xpmin = xp(no)
         if (yp(no) .gt. ypmax) ypmax = yp(no)
         if (yp(no) .lt. ypmin) ypmin = yp(no)
         if (zp(no) .gt. zpmax) zpmax = zp(no)
         if (zp(no) .lt. zpmin) zpmin = zp(no)
  550 continue
c
      write(*,*)' '
      write(*,*)'min and max of Power/Voronoi vertices coordinates:'
      write(*,*)' '
      write(*,*)'xmin = ',xpmin
      write(*,*)'xmax = ',xpmax
      write(*,*)'ymin = ',ypmin
      write(*,*)'ymax = ',ypmax
      write(*,*)'zmin = ',zpmin
      write(*,*)'zmax = ',zpmax
c
c----------------------------------------------------------------------
c
c     test Power/Voronoi vertices
c
      ng = 0
      derro = 0.0d0
      larg = 0
      do 800 i = 1, nt
         if(ifl(i).eq.0) then
            larg = larg+1
            go to 800
         endif
         ifir = icon(5,i)
         isec = icon(6,i)
         ithi = icon(7,i)
         ifou = icon(8,i)
         if(ifir.le.ng .or. isec.le.ng .or. ithi.le.ng .or.
     *      ifou.le.ng) stop 50
c
         xcor = x(ifir)-xp(i)
         ycor = y(ifir)-yp(i)
         zcor = z(ifir)-zp(i)
         wght= w(ifir)
         dnom=dmax1(dabs(xcor),dabs(ycor),dabs(zcor),dsqrt(dabs(wght)))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            zcor = zcor/dnom
            wght = (wght/dnom)/dnom
            dist1 = dnom*dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         else
            dist1 = dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         endif
c
         xcor = x(isec)-xp(i)
         ycor = y(isec)-yp(i)
         zcor = z(isec)-zp(i)
         wght= w(isec)
         dnom=dmax1(dabs(xcor),dabs(ycor),dabs(zcor),dsqrt(dabs(wght)))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            zcor = zcor/dnom
            wght = (wght/dnom)/dnom
            dist2 = dnom*dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         else
            dist2 = dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         endif
c
         xcor = x(ithi)-xp(i)
         ycor = y(ithi)-yp(i)
         zcor = z(ithi)-zp(i)
         wght= w(ithi)
         dnom=dmax1(dabs(xcor),dabs(ycor),dabs(zcor),dsqrt(dabs(wght)))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            zcor = zcor/dnom
            wght = (wght/dnom)/dnom
            dist3 = dnom*dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         else
            dist3 = dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         endif
c
         xcor = x(ifou)-xp(i)
         ycor = y(ifou)-yp(i)
         zcor = z(ifou)-zp(i)
         wght= w(ifou)
         dnom=dmax1(dabs(xcor),dabs(ycor),dabs(zcor),dsqrt(dabs(wght)))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            zcor = zcor/dnom
            wght = (wght/dnom)/dnom
            dist4 = dnom*dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         else
            dist4 = dsqrt(dabs(xcor**2 + ycor**2 + zcor**2 - wght))
         endif
c
         dista = dmax1(dist1,dist2,dist3,dist4)
         disti = dmin1(dist1,dist2,dist3,dist4)
         diffd = dista-disti
         if(diffd.gt.derro) derro = diffd
  800 continue
      write(*,*)' '
      write(*,*)'Number of Power/Voronoi vertices marked as too ',
     *          'large = ',larg
      write(*,*)' '
      write(*,*)'Power/Voronoi distance error = ',derro
      write(*,*)'(This number should be close to zero).'
  900 continue
c
c     test unbounded edges
c
      derro = 0.0d0
      ifir = 0
      do 1000 i = nt+1, nq
         call sitord(icon, ifir, i)
         isec = icon(6,i)
         ithi = icon(7,i)
         ifou = icon(8,i)
         dnom = dabs((x(isec)-x(ifou))*xp(i) + (y(isec)-y(ifou))*yp(i) +
     *          (z(isec)-z(ifou))*zp(i))
         if(dnom.gt.derro) derro = dnom
         dnom = dabs((x(ithi)-x(ifou))*xp(i) + (y(ithi)-y(ifou))*yp(i) +
     *          (z(ithi)-z(ifou))*zp(i))
         if(dnom.gt.derro) derro = dnom
 1000 continue
      write(*,*)' '
      write(*,*)'Power/Voronoi unbounded edge error = ',derro
      write(*,*)'(This number should be close to zero).'
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      stop
      end
*PWRVTX
c**********************************************************************
c
c     Driver subroutine of Fortran 77 program PWRVTX for computing
c     the Power/Voronoi diagram vertices and unbounded edges for a
c     set of points in 3-dimensional space from a Regular/Delaunay
c     tetrahedralization for the set.
c
c     It is assumed that a Regular/Delaunay tetrahedralization for
c     the set of points has already been computed using existing
c     program regtet/ragtat and that the output tetrahedron list
c     produced by regtet/ragtat contains only real tetrahedra (no
c     artificial ones), i. e. logical variable artfcl was set to
c     .false. during the execution of regtet/ragtat. The output
c     tetrahedron list from regtet/ragtat is then used by this
c     program. It is noted that because of degeneracies more than
c     one tetrahedron may correspond to the same Power/Voronoi vertex.
c     In addition, because of degeneracies, Power/Voronoi cells of
c     points may have facets that are not 2-dimensional, and if
c     weights are present there may be Power cells of points that are
c     not 3-dimensional.
c
c     Arrays xp, yp, zp, icon, ifl will contain output information.
c     For i=1,...,nt, xp(i), yp(i), zp(i) will be the coordinates of
c     the Power/Voronoi vertices.
c     For i=nt+1,...,nq, xp(i), yp(i), zp(i) will be the coordinates
c     of the unbounded edges.
c     For j=1,...,8, i=1,...,nt, icon(j,i) will be the input
c     tetrahedron list modified to include unbounded edges neighbor
c     information.
c     For j=1,...,8, i=nt+1,...,nq, icon(j,i) will be the unbounded
c     edges list which is an extension of the modified input
c     tetrahedron list and is obtained by treating unbounded edges
c     as tetrahedra with one vertex, vertex 0, at infinity.
c     For i=1,...,nt, ifl(i) will equal 1 if the size of Power/Voronoi
c     vertex i is not considered to be too large.
c     If the size of vertex i is considered to be too large ifl(i)
c     will equal 0.
c     For i=nt+1,...,nq, ifl(i) will equal 2.
c
      subroutine pwrvtx(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2,
     *                  xp, yp, zp, icon, ifl, is, nv, nw, nt, nq, nmax,
     *                  nvmax, wlenx, wleny, wlenz, wlenw, naddl, icfig,
     *                  iwfig, delaun, artfcl, reccor)
c
      double precision x(*), y(*), z(*), w(*)
      integer ix(*), iy(*), iz(*), iw(*)
      integer ix2(*), iy2(*), iz2(*), iw2(*)
      double precision xp(*), yp(*), zp(*)
      integer icon(8,*), ifl(*), is(*)
      integer nv, nw, nt, nq, nmax, nvmax, naddl, icfig, iwfig
      double precision dscle, wlenx, wleny, wlenz, wlenw
      logical delaun, artfcl, reccor
c
      double precision xmin, ymin, zmin, wmin, xmax, ymax, zmax
      integer isclp(2), isclw(2), isclr(2), mhalf, mfull
c
c     test for artificial points (not allowed)
c
      if(artfcl) then
         write(*,*)' '
         write(*,*)'Input tetrahedra contain artificial points.'
         write(*,*)'Such points are not allowed in this program.'
         write(*,*)'Program terminated.'
         stop 100
      endif
c
c     set up data structures properly
c
      call setsup(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2,
     *            xmin, ymin, zmin, wmin, xmax, ymax, zmax,
     *            nv, nw, nt, nmax, nvmax, dscle,
     *            wlenx, wleny, wlenz, wlenw, naddl,
     *            delaun, reccor, icfig, iwfig,
     *            mhalf, mfull, isclp, isclw, isclr)
c
c     compute Power/Voronoi vertices
c
      call vrtcmp(ix, iy, iz, iw, ix2, iy2, iz2, iw2, xp, yp, zp,
     *            icon, ifl, is, nv, nt, nq, nvmax, dscle, delaun,
     *            mhalf, mfull, isclp, isclw, isclr)
c
      return
      end
*SETSUP
c
c     This subroutine completes information about tetrahedralization.
c
      subroutine setsup(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2,
     *                  xmin, ymin, zmin, wmin, xmax, ymax, zmax,
     *                  nv, nw, nt,  nmax, nvmax, dscle,
     *                  wlenx, wleny, wlenz, wlenw, naddl,
     *                  delaun, reccor, icsfig, iwsfig,
     *                  mhalf, mfull, isclp, isclw, isclr)
c
      double precision x(*), y(*), z(*), w(*)
      integer ix(*), iy(*), iz(*), iw(*)
      integer ix2(*), iy2(*), iz2(*), iw2(*)
      double precision xmin, ymin, zmin, wmin, xmax, ymax, zmax
      integer nv, nw, nt, nmax, nvmax, naddl
      double precision dscle, wlenx, wleny, wlenz, wlenw
      integer icsfig, iwsfig, mhalf, mfull
      logical delaun, reccor
      integer isclp(*), isclw(*), isclr(*)
c
      double precision dscli, dfull, dfill, decml
      integer no, irec, irec1, nu, ng, i, naddm, j, icsfi2, irsfig
      integer isgcl, isclu, ibsfig, itsfig
      double precision xint, yint, zint, xcor, ycor, zcor
      double precision wmen, wbig, epz, wmax
c
c     initialize Fortran 77 word lengths
c
      mhalf=32768
      mfull=1073741824
c
c     testing number of points or vertices
c
      if (nv.lt.1 .or. nv.gt.nmax) stop 110
      if (nt.lt.1 .or. nt.gt.nvmax) stop 120
c
c     test variables associated with a possible rectangular polyhedron
c
      if(reccor)then
         if(wlenx.le.0.0d0 .or. wleny.le.0.0d0 .or. wlenz.le.0.0d0)
     *      stop 130
         if(naddl.lt.2) stop 140
      else
         wlenx = 0.0d0
         wleny = 0.0d0
         wlenz = 0.0d0
         wlenw = 0.0d0
         naddl = 0
      endif
c
c     calculating min and max
c
      xmax = x(1)
      xmin = x(1)
      ymax = y(1)
      ymin = y(1)
      zmax = z(1)
      zmin = z(1)
      wmax = w(1)
      wmin = w(1)
      do 15 no = 1, nv
         if (x(no) .gt. xmax) xmax = x(no)
         if (x(no) .lt. xmin) xmin = x(no)
         if (y(no) .gt. ymax) ymax = y(no)
         if (y(no) .lt. ymin) ymin = y(no)
         if (z(no) .gt. zmax) zmax = z(no)
         if (z(no) .lt. zmin) zmin = z(no)
         if (w(no) .gt .wmax) wmax = w(no)
         if (w(no) .lt. wmin) wmin = w(no)
   15 continue
c
c     shift data
c
      irec = 0
      if(reccor) irec = irec + 6*(naddl**2) - 12*naddl + 8
      if(irec.eq.0) go to 25
      irec1 = irec + 1
      nv = nv + irec
      if(nv .gt. nmax) stop 150
      do 20 no = nv, irec1, -1
         nu = no - irec
         x(no) = x(nu)
         y(no) = y(nu)
         z(no) = z(nu)
         w(no) = w(nu)
   20 continue
   25 continue
      if(nv.ne.nw) stop 160
c
c     compute corners of rectangular polyhedron
c
      if(.not.reccor) go to 165
      ng = 0
      x(ng+1) = xmin - wlenx
      y(ng+1) = ymin - wleny
      z(ng+1) = zmax + wlenz
c
      x(ng+2) = xmin - wlenx
      y(ng+2) = ymax + wleny
      z(ng+2) = zmax + wlenz
c
      x(ng+3) = xmax + wlenx
      y(ng+3) = ymax + wleny
      z(ng+3) = zmax + wlenz
c
      x(ng+4) = xmax + wlenx
      y(ng+4) = ymin - wleny
      z(ng+4) = zmax + wlenz
c
      x(ng+5) = xmin - wlenx
      y(ng+5) = ymin - wleny
      z(ng+5) = zmin - wlenz
c
      x(ng+6) = xmin - wlenx
      y(ng+6) = ymax + wleny
      z(ng+6) = zmin - wlenz
c
      x(ng+7) = xmax + wlenx
      y(ng+7) = ymax + wleny
      z(ng+7) = zmin - wlenz
c
      x(ng+8) = xmax + wlenx
      y(ng+8) = ymin - wleny
      z(ng+8) = zmin - wlenz
c
      do 60 i=ng+1,ng+8
         if((x(i).ge.xmin.and.x(i).le.xmax).or.(y(i).ge.ymin.and.
     *       y(i).le.ymax).or.(z(i).ge.zmin.and.z(i).le.zmax)) stop 260
   60 continue
c
      xmin = xmin - wlenx
      ymin = ymin - wleny
      zmin = zmin - wlenz
      xmax = xmax + wlenx
      ymax = ymax + wleny
      zmax = zmax + wlenz
c
      if(naddl.eq.2) go to 155
c
c     compute other points in grid on surface of polyhedron
c
      naddm = naddl-2
      xint = (xmax-xmin)/dble(naddl-1)
      yint = (ymax-ymin)/dble(naddl-1)
      zint = (zmax-zmin)/dble(naddl-1)
      ng = ng+8
c
      do 119 i = 1, naddm
         xcor = xmin + dble(i)*xint
         ng = ng + 4
         x(ng-3) = xcor
         y(ng-3) = ymin
         z(ng-3) = zmin
         x(ng-2) = xcor
         y(ng-2) = ymin
         z(ng-2) = zmax
         x(ng-1) = xcor
         y(ng-1) = ymax
         z(ng-1) = zmin
         x(ng) = xcor
         y(ng) = ymax
         z(ng) = zmax
  119 continue
c
      do 120 i = 1, naddm
         ycor = ymin + dble(i)*yint
         ng = ng + 4
         y(ng-3) = ycor
         z(ng-3) = zmin
         x(ng-3) = xmin
         y(ng-2) = ycor
         z(ng-2) = zmin
         x(ng-2) = xmax
         y(ng-1) = ycor
         z(ng-1) = zmax
         x(ng-1) = xmin
         y(ng) = ycor
         z(ng) = zmax
         x(ng) = xmax
  120 continue
c
      do 121 i = 1, naddm
         zcor = zmin + dble(i)*zint
         ng = ng + 4
         z(ng-3) = zcor
         x(ng-3) = xmin
         y(ng-3) = ymin
         z(ng-2) = zcor
         x(ng-2) = xmin
         y(ng-2) = ymax
         z(ng-1) = zcor
         x(ng-1) = xmax
         y(ng-1) = ymin
         z(ng) = zcor
         x(ng) = xmax
         y(ng) = ymax
  121 continue
c
      do 130 i = 1, naddm
         xcor = xmin + dble(i)*xint
         do 125 j = 1, naddm
            ycor = ymin + dble(j)*yint
            ng = ng + 2
            x(ng-1) = xcor
            y(ng-1) = ycor
            z(ng-1) = zmin
            x(ng) = xcor
            y(ng) = ycor
            z(ng) = zmax
  125    continue
  130 continue
c
      do 140 i = 1, naddm
         ycor = ymin + dble(i)*yint
         do 135 j = 1, naddm
            zcor = zmin + dble(j)*zint
            ng = ng + 2
            y(ng-1) = ycor
            z(ng-1) = zcor
            x(ng-1) = xmin
            y(ng) = ycor
            z(ng) = zcor
            x(ng) = xmax
  135    continue
  140 continue
c
      do 150 i = 1, naddm
         zcor = zmin + dble(i)*zint
         do 145 j = 1, naddm
            xcor = xmin + dble(j)*xint
            ng = ng + 2
            z(ng-1) = zcor
            x(ng-1) = xcor
            y(ng-1) = ymin
            z(ng) = zcor
            x(ng) = xcor
            y(ng) = ymax
  145    continue
  150 continue
c
      if(ng.ne.irec) stop 280
c
  155 continue
      wmen = wmin
      if(.not.delaun .and. reccor) then
         wmen = wmen - wlenw
         do 160 i=irec1,irec
            w(i) = wmen
  160    continue
      endif
      if(wmen.lt.wmin) wmin = wmen
      if(wmen.gt.wmax) wmax = wmen
c
c     test # of significant figures of nondecimal part of coordinates
c
  165 continue
      write(*,*)' '
      write(*,*)'min and max of points coordinates:'
      write(*,*)' '
      write(*,*)'xmin = ',xmin
      write(*,*)'xmax = ',xmax
      write(*,*)'ymin = ',ymin
      write(*,*)'ymax = ',ymax
      write(*,*)'zmin = ',zmin
      write(*,*)'zmax = ',zmax
      epz = 0.01d0
      wbig = 0.0d0
      if(wbig .lt. dabs(xmax)) wbig = dabs(xmax)
      if(wbig .lt. dabs(xmin)) wbig = dabs(xmin)
      if(wbig .lt. dabs(ymax)) wbig = dabs(ymax)
      if(wbig .lt. dabs(ymin)) wbig = dabs(ymin)
      if(wbig .lt. dabs(zmax)) wbig = dabs(zmax)
      if(wbig .lt. dabs(zmin)) wbig = dabs(zmin)
      wbig = wbig + epz
c     WRITE(*,*)'COORDINATES WBIG=',WBIG
      ibsfig = 0
  180 continue
      ibsfig = ibsfig+1
      wbig = wbig/10.0d0
      if(wbig .ge. 1.0d0) go to 180
      if(ibsfig.gt.9) then
         write(*,*)'Number of significant figures of largest ',
     *             'nondecimal part of'
         write(*,*)'a point coordinate appears to be greater than 9.'
         write(*,*)'Program is terminated.'
         stop 282
      endif
      itsfig = ibsfig + icsfig
c     WRITE(*,*)'ITSFIG=',ITSFIG,' IBSFIG=',IBSFIG,' ICSFIG=',ICSFIG
      if(itsfig.gt.14) then
         write(*,*)' '
         write(*,*)'For this execution of the program the largest ',
     *             'total number of'
         write(*,*)'significant figures ',
     *             'that a coordinate requires appears to be ',itsfig
         write(*,*)'Program is terminated since the maximum ',
     *             'allowed is 14.'
         stop 284
      endif
c
c     if a Regular tetrahedralization test # of significant figures
c     of nondecimal part of weights
c
      if(delaun) go to 200
      wbig = 0.0d0
      if(wbig .lt. dabs(wmax)) wbig = dabs(wmax)
      if(wbig .lt. dabs(wmin)) wbig = dabs(wmin)
      wbig = wbig + epz
c     WRITE(*,*)'COORDINATES WBIG=',WBIG
      ibsfig = 0
  190 continue
      ibsfig = ibsfig+1
      wbig = wbig/10.0d0
      if(wbig .ge. 1.0d0) go to 190
      if(ibsfig.gt.9) then
         write(*,*)'Number of significant figures of largest ',
     *             'nondecimal part of'
         write(*,*)'a weight appears to be greater than 9.'
         write(*,*)'Program is terminated.'
         stop 286
      endif
      itsfig = ibsfig + iwsfig
c     WRITE(*,*)'ITSFIG=',ITSFIG,' IBSFIG=',IBSFIG,' IWSFIG=',IWSFIG
      if(itsfig.gt.14) then
         write(*,*)' '
         write(*,*)'For this execution of the program the largest ',
     *             'total number of'
         write(*,*)'significant figures ',
     *             'that a weight requires appears to be ',itsfig
         write(*,*)'Program is terminated since the maximum ',
     *             'allowed is 14.'
         stop 288
      endif
  200 continue
c
c     test number of significant figures of decimal part of coordinates
c
      if(icsfig.lt.0 .or. icsfig.gt.9) stop 290
      isclu = 1
      dscle = 1.0d0
      if(icsfig.eq.0) go to 220
      do 210 i = 1, icsfig
         isclu = 10*isclu
         dscle = 10.0d0*dscle
  210 continue
  220 continue
      if(iabs(isclu).ge.mfull) stop 295
      call decomp(isclp, isgcl, isclu, mhalf)
      if(isgcl.ne.1) stop 297
c
c     test lengths of x, y, z-coordinates, shift and make them integers
c
      dfull = dble(mfull)
      dfill=dfull/dscle
      do 235 i = 1, nv
         ix2(i) = 0
         iy2(i) = 0
         iz2(i) = 0
         if(dabs(x(i)).lt.dfill) then
            ix(i) = idnint(dscle*x(i))
            if(iabs(ix(i)).lt.mfull) then
               x(i) = dble(ix(i))/dscle
               go to 225
            endif
         endif
         if(dabs(x(i)).ge.dfull) stop 305
         ix(i) = idint(x(i))
         if(iabs(ix(i)).ge.mfull) stop 310
         decml = (x(i) - dint(x(i)))*dscle
         if(dabs(decml).ge.dfull) stop 312
         ix2(i) = idnint(decml)
         if(iabs(ix2(i)).ge.mfull) stop 315
         if(iabs(ix2(i)).eq.0) then
            x(i) = dble(ix(i))
            ix2(i) = mfull
         else
            x(i) = dble(ix(i)) + (dble(ix2(i))/dscle)
         endif
  225    continue
         if(dabs(y(i)).lt.dfill) then
            iy(i) = idnint(dscle*y(i))
            if(iabs(iy(i)).lt.mfull) then
               y(i) = dble(iy(i))/dscle
               go to 230
            endif
         endif
         if(dabs(y(i)).ge.dfull) stop 320
         iy(i) = idint(y(i))
         if(iabs(iy(i)).ge.mfull) stop 325
         decml = (y(i) - dint(y(i)))*dscle
         if(dabs(decml).ge.dfull) stop 327
         iy2(i) = idnint(decml)
         if(iabs(iy2(i)).ge.mfull) stop 330
         if(iabs(iy2(i)).eq.0) then
            y(i) = dble(iy(i))
            iy2(i) = mfull
         else
            y(i) = dble(iy(i)) + (dble(iy2(i))/dscle)
         endif
  230    continue
         if(dabs(z(i)).lt.dfill) then
            iz(i) = idnint(dscle*z(i))
            if(iabs(iz(i)).lt.mfull) then
               z(i) = dble(iz(i))/dscle
               go to 235
            endif
         endif
         if(dabs(z(i)).ge.dfull) stop 335
         iz(i) = idint(z(i))
         if(iabs(iz(i)).ge.mfull) stop 340
         decml = (z(i) - dint(z(i)))*dscle
         if(dabs(decml).ge.dfull) stop 342
         iz2(i) = idnint(decml)
         if(iabs(iz2(i)).ge.mfull) stop 345
         if(iabs(iz2(i)).eq.0) then
            z(i) = dble(iz(i))
            iz2(i) = mfull
         else
            z(i) = dble(iz(i)) + (dble(iz2(i))/dscle)
         endif
  235 continue
c
c     if a Regular tetrahedralization test number of significant figures
c     of decimal part of weights, test lengths of weights, shift and
c     make them integers
c
      if(delaun) go to 300
      icsfi2 = 2*icsfig
      irsfig = icsfi2 - iwsfig
      if(iwsfig.lt.0 .or. iwsfig.gt.9 .or. irsfig.lt.0 .or. irsfig.gt.9)
     *   stop 350
      isclu = 1
      dscli = 1.0d0
      if(iwsfig.eq.0) go to 250
      do 240 i = 1, iwsfig
         isclu = 10*isclu
         dscli = 10.0d0*dscli
  240 continue
  250 continue
      if(iabs(isclu).ge.mfull) stop 360
      call decomp(isclw, isgcl, isclu, mhalf)
      if(isgcl.ne.1) stop 363
      dfill = dfull/dscli
      do 260 i = 1, nv
         iw2(i) = 0
         if(dabs(w(i)).lt.dfill) then
            iw(i) = idnint(dscli*w(i))
            if(iabs(iw(i)).lt.mfull) then
               w(i) = dble(iw(i))/dscli
               go to 260
            endif
         endif
         if(dabs(w(i)).ge.dfull) stop 365
         iw(i) = idint(w(i))
         if(iabs(iw(i)).ge.mfull) stop 370
         decml = (w(i) - dint(w(i)))*dscli
         if(dabs(decml).ge.dfull) stop 372
         iw2(i) = idnint(decml)
         if(iabs(iw2(i)).ge.mfull) stop 375
         if(iabs(iw2(i)).eq.0) then
            w(i) = dble(iw(i))
            iw2(i) = mfull
         else
            w(i) = dble(iw(i)) + (dble(iw2(i))/dscli)
         endif
  260 continue
      isclu = 1
      if(irsfig.eq.0) go to 290
      do 270 i = 1, irsfig
         isclu = 10*isclu
  270 continue
  290 continue
      if(iabs(isclu).ge.mfull) stop 385
      call decomp(isclr, isgcl, isclu, mhalf)
      if(isgcl.ne.1) stop 390
  300 continue
c
      return
      end
*VRTCMP
c
c     This subroutine computes Power/Voronoi vertices.
c 
c
      subroutine vrtcmp(ix, iy, iz, iw, ix2, iy2, iz2, iw2, xp, yp, zp,
     *                  icon, ifl, is, nv, nt, nq, nvmax, dscle, delaun,
     *                  mhalf, mfull, isclp, isclw, isclr)
c
      integer ix(*), iy(*), iz(*), iw(*)
      integer ix2(*), iy2(*), iz2(*), iw2(*)
      double precision xp(*), yp(*), zp(*)
      integer icon(8,*), ifl(*), is(*)
      integer nv, nt, nq, nvmax
      integer mhalf, mfull, nkmax
      logical delaun
      integer isclp(*), isclw(*), isclr(*), ikon(8,1), a, b, c
      parameter(nkmax=30)
      integer io(nkmax), inx(nkmax), iny(nkmax), inz(nkmax)
      double precision r215, dnom, xnum, ynum, znum
      double precision deps, dscle, dscl2, dfull
c
      integer ng, i, j, ifir, isec, ithi, ifou, isgo, iko
      integer isgnx, isgny, isgnz, iknx, ikny, iknz
      integer ist1, ist2, ist3, icur, icar, now
c
c     compute vertices
c
      ng = 0
      dfull = dble(mfull)
      r215 = dble(mhalf)
      deps = dble(0.9)
      if(dscle.lt.deps) stop 400
      dscl2 = 2.0d0*dscle
      write(*,*)' '
      do 1000 i = 1, nt
         if(i.le.(i/10000)*10000)
     *   write(*,*)'Number of computed vertices = ',i
         xp(i) = 0.0d0
         yp(i) = 0.0d0
         zp(i) = 0.0d0
         ifir = icon(5,i)
         isec = icon(6,i)
         ithi = icon(7,i)
         ifou = icon(8,i)
         if(ifir.le.ng .or. isec.le.ng .or. ithi.le.ng .or.
     *      ifou.le.ng) stop 405
         call idenom(ix, iy, iz, ix2, iy2, iz2, ifir, isec, ithi,
     *               ifou, mhalf, mfull, isclp, io, isgo, iko)
         if(isgo.le.0) stop 410
         call xyznum(ix, iy, iz, iw, ix2, iy2, iz2, iw2, ifir, isec,
     *               ithi, ifou, mhalf, mfull, isclp, isclw, isclr,
     *               delaun, inx, iny, inz, isgnx, isgny, isgnz,
     *               iknx, ikny, iknz)
         call doubnm(io, isgo, iko, r215, dnom)
         if(dnom.lt.deps) stop 420
         call doubnm(inx, isgnx, iknx, r215, xnum)
         call doubnm(iny, isgny, ikny, r215, ynum)
         call doubnm(inz, isgnz, iknz, r215, znum)
         xp(i) = (xnum/dnom)/dscl2
         yp(i) = (ynum/dnom)/dscl2
         zp(i) = (znum/dnom)/dscl2
         ifl(i) = 1
 1000 continue
c
c     compute unbounded edges
c
      nq = nt
      do 1300 i = 1, nt
         do 1250 j = 1, 4
            if(icon(j,i).ne.0) go to 1250
            nq = nq + 1
            if(nq.gt.nvmax) stop 430
            icon(j,i) = nq
            if(j.eq.1)then
               ist1 = icon(6,i)
               ist2 = icon(7,i)
               ist3 = icon(8,i)
            elseif(j.eq.2)then
               ist1 = icon(7,i)
               ist2 = icon(5,i)
               ist3 = icon(8,i)
            elseif(j.eq.3)then
               ist1 = icon(8,i)
               ist2 = icon(5,i)
               ist3 = icon(6,i)
            else
               ist1 = icon(5,i)
               ist2 = icon(7,i)
               ist3 = icon(6,i)
            endif
            icon(1,nq) = i
            icon(5,nq) = 0
            icon(6,nq) = ist3
            icon(7,nq) = ist2
            icon(8,nq) = ist1
            call crossp(ix, iy, iz, ix2, iy2, iz2, ist1, ist2, ist3,
     *                  mhalf, mfull, isclp, inx, isgnx, iknx,
     *                  iny, isgny, ikny, inz, isgnz, iknz)
            call doubnm(inx, isgnx, iknx, r215, xnum)
            call doubnm(iny, isgny, ikny, r215, ynum)
            call doubnm(inz, isgnz, iknz, r215, znum)
            dnom = dmax1(dabs(xnum),dabs(ynum),dabs(znum))
            if(dnom.lt.deps) stop 435
            xnum = xnum/dnom
            ynum = ynum/dnom
            znum = znum/dnom
            dnom = dsqrt(xnum**2+ynum**2+znum**2)
            if(dnom.lt.deps) stop 440
            xp(nq) = xnum/dnom
            yp(nq) = ynum/dnom
            zp(nq) = znum/dnom
            ifl(nq) = 2
 1250    continue
 1300 continue
c
c     connect unbounded edges
c
      if(nq.le.nt) stop 445
      do 1500 icur = nt+1, nq
         icar = icon(1,icur)
         a = icon(8,icur)
         b = icon(7,icur)
         c = icon(6,icur)
         do 1350 j = 1, 8
            ikon(j,1) = icon(j,icar)
 1350    continue
         ist1 = a
         ist2 = b
         ist3 = c
         call reordr(ikon,ist1,ist2,1)
         do 1375 j = 2, 4
            now = ikon(4,1)
 1360       continue
            if(now.eq.icur .or. now.eq.icar .or. now.eq.0) stop 450
            if(ifl(now).eq.2) then
               icon(j,icur) = now
            elseif(ifl(now).eq.1) then
               call reordr(icon,ist1,ist2,now)
               now = icon(4,now)
               go to 1360
            else
               stop 480
            endif
            if(j.eq.2) then
               ist1 = c
               ist2 = a
               ist3 = b
               call reordr(ikon,ist1,ist2,1)
            endif
            if(j.eq.3) then
               ist1 = b
               ist2 = c
               ist3 = a
               call reordr(ikon,ist1,ist2,1)
            endif
 1375    continue
 1500 continue
c
c     test consistency of tetrahedralization
c
      call consis(icon, is, ifl, nv, nq)
c
c     redefine array ifl
c
      do 1600 i = 1, nt
         ifl(i) = 1
         dnom = dmax1(dabs(xp(i)),dabs(yp(i)),dabs(zp(i)))
         if(dnom.gt.dfull) ifl(i) = 0
 1600 continue
      do 1700 i = nt+1, nq
         ifl(i) = 2
 1700 continue
c
      return
      end
*CONSIS
c
c     subroutine consis to -
c
c     test consistency of diagram
c
c     May 1, 1989
c
      subroutine consis(icon, is, ifl, n, ivnxt)
c
      integer icon(8,*), is(*), ifl(*), ikon(8,1)
      integer site0, site1, site2, site3, n, ivnxt
      integer i, iscur, isone, islst, isini, indx
c
c     test initial tetrahedron for each site
c
      do 50 i = 1, n
          iscur = is(i)
          if (iscur .le. 0) goto 50
          if(icon(5,iscur) .ne. i .and. icon(6,iscur) .ne. i .and.
     *       icon(7,iscur) .ne. i .and. icon(8,iscur) .ne. i) stop 2810
   50 continue
c
c     initialize
c
      isone = 1
      do 60 i = 1, n
         if(is(i) .gt. 0) go to 80
   60 continue
      stop 2820
   80 continue
      islst = is(i)
      isini = islst
c
      do 100 i = 1, ivnxt
          ifl(i) = 0
  100 continue
c
      ifl(isini) = 1
      indx = 1
      iscur = icon(1,isini)
      if(iscur.eq.0) go to 500
      site0 = icon(5,isini)
      site1 = icon(6,isini)
      site2 = icon(7,isini)
      site3 = icon(8,isini)
c
c     reorder iscur relative to site1 and site2, and test
c
  200 continue
      if(site0.eq.site1 .or. site0.eq.site2 .or. site0.eq.site3 .or.
     *   site1.eq.site2 .or. site1.eq.site3 .or. site2.eq.site3)
     *   stop 2830
      call reordr(icon, site1, site2, iscur)
      if(icon(7,iscur) .ne. site3) stop 2840
      if(icon(4,iscur) .ne. islst) stop 2850
      if(icon(8,iscur) .eq. site0) stop 2855
      ifl(iscur) = 1
c
c     obtain next tetrahedron
c
      islst = iscur
      indx = 1
      iscur = icon(1,islst)
      if(iscur.eq.0) go to 500
      site0 = icon(5,islst)
      site1 = icon(6,islst)
      site2 = icon(7,islst)
      site3 = icon(8,islst)
      if(ifl(iscur) .ne. 1) go to 200
c
c     reorder iscur relative to site1 and site2, and test
c
  300 continue
      if(site0.eq.site1 .or. site0.eq.site2 .or. site0.eq.site3 .or.
     *   site1.eq.site2 .or. site1.eq.site3 .or. site2.eq.site3)
     *   stop 2860
      do 400 i = 1, 8
          ikon(i,1) = icon(i,iscur)
  400 continue
      call reordr(ikon, site1, site2, isone)
      if(ikon(7,1) .ne. site3) stop 2865
      if(ikon(4,1) .ne. islst) stop 2870
      if(ikon(8,1) .eq. site0) stop 2875
c
c     obtain next tetrahedron
c
  500 continue
      if(indx.eq.1) then
          indx = 2
          iscur = icon(2,islst)
          if(iscur.eq.0) go to 500
          site0 = icon(6,islst)
          site1 = icon(5,islst)
          site2 = icon(8,islst)
          site3 = icon(7,islst)
          if(ifl(iscur) .ne. 1) go to 200
          go to 300
      elseif(indx.eq.2) then
          indx = 3
          iscur = icon(3,islst)
          if(iscur.eq.0) go to 500
          site0 = icon(7,islst)
          site1 = icon(5,islst)
          site2 = icon(6,islst)
          site3 = icon(8,islst)
          if(ifl(iscur) .ne. 1) go to 200
          go to 300
      elseif(indx.eq.3) then
          if(islst .ne. isini) then
              iscur = islst
              islst = icon(4,iscur)
              if(islst. eq. 0) stop 2880
              if(icon(1,islst) .eq. iscur) then
                  indx = 1
              elseif(icon(2,islst) .eq. iscur) then
                  indx = 2
              elseif(icon(3,islst) .eq. iscur) then
                  indx = 3
              elseif(icon(4,islst) .eq. iscur) then
                  indx = 4
              else
                  stop 2885
              endif
              go to 500
          else
              indx = 4
              iscur = icon(4,islst)
              if(iscur.eq.0) go to 500
              site0 = icon(8,islst)
              site1 = icon(5,islst)
              site2 = icon(7,islst)
              site3 = icon(6,islst)
              if(ifl(iscur) .ne. 1) go to 200
              go to 300
          endif
      endif
      if(islst .ne. isini) stop 2890
c
c     write (*,*) ' '
c     write (*,*) '**************************************'
c     write (*,*) 'consistency check satisfied'
c     write (*,*) '**************************************'
c     write (*,*) ' '
c
      return
      end
*REORDR
c
c     subroutine reordr to -
c
c     reorder icon(i,iscur), i = 1, ..., 8, so that site1 equals
c     icon(5,iscur) and site2 equals icon(6,iscur)
c
c     July 22, 1988
c
      subroutine reordr(icon, site1, site2, iscur)
c
      integer icon(8,*), site1, site2, iscur, itemp
c
      if(icon(5,iscur) .eq. site1) go to 200
      if(icon(6,iscur) .eq. site1) then
          itemp = icon(1,iscur)
          icon(1,iscur) = icon(2,iscur)
          icon(2,iscur) = icon(4,iscur)
          icon(4,iscur) = itemp
          itemp = icon(5,iscur)
          icon(5,iscur) = icon(6,iscur)
          icon(6,iscur) = icon(8,iscur)
          icon(8,iscur) = itemp
      elseif(icon(7,iscur) .eq. site1) then
          itemp = icon(1,iscur)
          icon(1,iscur) = icon(3,iscur)
          icon(3,iscur) = icon(2,iscur)
          icon(2,iscur) = itemp
          itemp = icon(5,iscur)
          icon(5,iscur) = icon(7,iscur)
          icon(7,iscur) = icon(6,iscur)
          icon(6,iscur) = itemp
      elseif(icon(8,iscur) .eq. site1) then
          itemp = icon(1,iscur)
          icon(1,iscur) = icon(4,iscur)
          icon(4,iscur) = icon(3,iscur)
          icon(3,iscur) = itemp
          itemp = icon(5,iscur)
          icon(5,iscur) = icon(8,iscur)
          icon(8,iscur) = icon(7,iscur)
          icon(7,iscur) = itemp
      else
          stop 2910
      endif
  200 continue
c
      if(icon(6,iscur) .eq. site2) go to 300
      if(icon(7,iscur) .eq. site2) then
          itemp = icon(2,iscur)
          icon(2,iscur) = icon(3,iscur)
          icon(3,iscur) = icon(4,iscur)
          icon(4,iscur) = itemp
          itemp = icon(6,iscur)
          icon(6,iscur) = icon(7,iscur)
          icon(7,iscur) = icon(8,iscur)
          icon(8,iscur) = itemp
      elseif(icon(8,iscur) .eq. site2) then
          itemp = icon(2,iscur)
          icon(2,iscur) = icon(4,iscur)
          icon(4,iscur) = icon(3,iscur)
          icon(3,iscur) = itemp
          itemp = icon(6,iscur)
          icon(6,iscur) = icon(8,iscur)
          icon(8,iscur) = icon(7,iscur)
          icon(7,iscur) = itemp
      else
          stop 2920
      endif
  300 continue
c
      return
      end
*SITORD
c
c     subroutine sitord to -
c
c     reorder icon(i,iscur), i = 1, ..., 8, so that site1 equals
c     icon(5,iscur) 
c
c     July 22, 1988
c
      subroutine sitord(icon, site1, iscur)
c
      integer icon(8,*), site1, iscur, itemp
c
      if(icon(5,iscur) .eq. site1) go to 200
      if(icon(6,iscur) .eq. site1) then
          itemp = icon(1,iscur)
          icon(1,iscur) = icon(2,iscur)
          icon(2,iscur) = icon(4,iscur)
          icon(4,iscur) = itemp
          itemp = icon(5,iscur)
          icon(5,iscur) = icon(6,iscur)
          icon(6,iscur) = icon(8,iscur)
          icon(8,iscur) = itemp
      elseif(icon(7,iscur) .eq. site1) then
          itemp = icon(1,iscur)
          icon(1,iscur) = icon(3,iscur)
          icon(3,iscur) = icon(2,iscur)
          icon(2,iscur) = itemp
          itemp = icon(5,iscur)
          icon(5,iscur) = icon(7,iscur)
          icon(7,iscur) = icon(6,iscur)
          icon(6,iscur) = itemp
      elseif(icon(8,iscur) .eq. site1) then
          itemp = icon(1,iscur)
          icon(1,iscur) = icon(4,iscur)
          icon(4,iscur) = icon(3,iscur)
          icon(3,iscur) = itemp
          itemp = icon(5,iscur)
          icon(5,iscur) = icon(8,iscur)
          icon(8,iscur) = icon(7,iscur)
          icon(7,iscur) = itemp
      else
          stop 3010
      endif
  200 continue
      return
      end
*IDENOM
c
c     Routine for determining denominator used in formulas for
c     computing x-,y-z-coordinates of orthogonal center of
c     tetrahedron with vertices ifir, isec, ithi, ifou.
c
      subroutine idenom(x, y, z, x2, y2, z2, ifir, isec, ithi,
     *                  ifou, mhalf, mfull, isclp, io, isgo, iko)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*)
      integer ifir, isec, ithi, ifou
      integer isclp(*), mhalf, mfull, nkmax
      parameter (nkmax = 30)
      integer io(*), iu(nkmax), iv(nkmax), iw(nkmax)
      integer ix2(nkmax), iy2(nkmax), iz2(nkmax)
      integer ix3(nkmax), iy3(nkmax), iz3(nkmax)
      integer ix4(nkmax), iy4(nkmax), iz4(nkmax)
      integer ixf(nkmax), iyf(nkmax), izf(nkmax)
      integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew
      integer ixthw, iythw, izthw, ixfow, iyfow, izfow
      integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2
      integer ixth2, iyth2, izth2, ixfo2, iyfo2, izfo2
      integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf
      integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2
      integer isgx3, isgy3, isgz3, ikx3, iky3, ikz3
      integer isgx4, isgy4, isgz4, ikx4, iky4, ikz4
      integer isgo, isgu, isgv, isgw, iko, iku, ikv, ikw
c
      ixfiw = x(ifir)
      iyfiw = y(ifir)
      izfiw = z(ifir)
      ixsew = x(isec)
      iysew = y(isec)
      izsew = z(isec)
      ixthw = x(ithi)
      iythw = y(ithi)
      izthw = z(ithi)
      ixfow = x(ifou)
      iyfow = y(ifou)
      izfow = z(ifou)
c
      ixfi2 = x2(ifir)
      iyfi2 = y2(ifir)
      izfi2 = z2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
      izse2 = z2(isec)
      ixth2 = x2(ithi)
      iyth2 = y2(ithi)
      izth2 = z2(ithi)
      ixfo2 = x2(ifou)
      iyfo2 = y2(ifou)
      izfo2 = z2(ifou)
c
      call decmp2(ixf, isgxf, ikxf, ixfiw, ixfi2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfiw, iyfi2, mhalf, mfull, isclp)
      call decmp2(izf, isgzf, ikzf, izfiw, izfi2, mhalf, mfull, isclp)
c
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix2, isgo, isgxf, isgx2, iko, ikxf, ikx2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy2, isgo, isgyf, isgy2, iko, ikyf, iky2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izsew, izse2, mhalf, mfull, isclp)
      call muldif(io, izf, iz2, isgo, isgzf, isgz2, iko, ikzf, ikz2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, ixthw, ixth2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix3, isgo, isgxf, isgx3, iko, ikxf, ikx3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iythw, iyth2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy3, isgo, isgyf, isgy3, iko, ikyf, iky3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izthw, izth2, mhalf, mfull, isclp)
      call muldif(io, izf, iz3, isgo, isgzf, isgz3, iko, ikzf, ikz3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, ixfow, ixfo2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix4, isgo, isgxf, isgx4, iko, ikxf, ikx4,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iyfow, iyfo2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy4, isgo, isgyf, isgy4, iko, ikyf, iky4,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izfow, izfo2, mhalf, mfull, isclp)
      call muldif(io, izf, iz4, isgo, isgzf, isgz4, iko, ikzf, ikz4,
     *              nkmax, mhalf)
c
      call mulmul(iy2, iz3, iv, isgy2, isgz3, isgv, iky2, ikz3, ikv,
     *              nkmax, mhalf)
      call mulmul(iz2, iy3, iu, isgz2, isgy3, isgu, ikz2, iky3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, iw, isgv, isgu, isgw, ikv, iku, ikw,
     *              nkmax, mhalf)
      call mulmul(iw, ix4, io, isgw, isgx4, isgo, ikw, ikx4, iko,
     *              nkmax, mhalf)
c
      call mulmul(iz2, ix3, iv, isgz2, isgx3, isgv, ikz2, ikx3, ikv,
     *              nkmax, mhalf)
      call mulmul(ix2, iz3, iu, isgx2, isgz3, isgu, ikx2, ikz3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, iw, isgv, isgu, isgw, ikv, iku, ikw,
     *              nkmax, mhalf)
      call mulmul(iw, iy4, iu, isgw, isgy4, isgu, ikw, iky4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(io, iu, iw, isgo, isgu, isgw, iko, iku, ikw,
     *              nkmax, mhalf)
c
      call mulmul(ix2, iy3, iv, isgx2, isgy3, isgv, ikx2, iky3, ikv,
     *              nkmax, mhalf)
      call mulmul(iy2, ix3, iu, isgy2, isgx3, isgu, iky2, ikx3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
      call mulmul(io, iz4, iu, isgo, isgz4, isgu, iko, ikz4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iw, iu, io, isgw, isgu, isgo, ikw, iku, iko,
     *              nkmax, mhalf)
c
      return
      end
*DECMP2
c
c     subroutine decmp2 
c
c     to decompose a regular or non-regular length integer
c
      subroutine decmp2(ia, isga, ika, iwi, iwi2, mhalf, mfull, isclp)
c
      integer ia(*), isga, ika, iwi, iwi2, mhalf, mfull, isclp(*)
      integer nkmax
      parameter (nkmax = 30)
      integer iu(nkmax), io(nkmax), isgu, isgo, iku, iko, isgcl, ikcl
c
      call decomp(ia, isga, iwi, mhalf)
      ika = 2
      if(iwi2.ne.0) then
         isgcl = 1
         ikcl = 2
         call mulmul(ia, isclp, iu, isga, isgcl, isgu, ika, ikcl,
     *                 iku, nkmax, mhalf)
         if(iwi2.eq.mfull) iwi2 = 0
         call decomp(io, isgo, iwi2, mhalf)
         isgo = -isgo
         iko = 2
         call muldif(iu, io, ia, isgu, isgo, isga, iku, iko, ika,
     *                 nkmax, mhalf)
      endif
c
      return
      end
*DECOMP
c
c     subroutine decomp
c
c     to decompose a regular length integer
c
c     iwi = isga*(ia(1) + ia(2) * mhalf)
c
c     iwi  is a regular length integer
c     isga is a sign integer (-1, 0, 1)
c     ia(1) and ia(2) are integers less than mhalf
c
      subroutine decomp(ia, isga, iwi, mhalf)
c
      integer ia(*), isga, iwi, mhalf, ivi
c
      if(iwi.gt.0) then
         isga = 1
         ivi = iwi
      elseif(iwi.lt.0) then
         isga =-1
         ivi = -iwi
      else
         isga = 0
         ia(1) = 0
         ia(2) = 0
         return
      endif
      ia(2) = ivi/mhalf
      ia(1) = ivi - ia(2)*mhalf
c
      return
      end
*MULMUL
c
c     subroutine mulmul
c
c     to perform a multiple precision integer multiplication
c     (for multiplying 2 or more integers)
c
c     io = ia * ib
c
c     ia represents a decomposed integer
c     ib represents a decomposed integer
c     io is the product of ia and ib in its decomposed form
c
      subroutine mulmul(ia, ib, io, isga, isgb, isgo, ika, ikb, iko,
     *                    nkmax, mhalf)
c
      integer ia(*), ib(*), io(*)
      integer isga, isgb, isgo, ika, ikb, iko, nkmax, mhalf
      integer i, ipt, ipr, iko1, k, j
c
      if(isga.eq.0.or.isgb.eq.0)then
         isgo=0
         iko = 2
         io(1) = 0
         io(2) = 0
         return
      endif
c
      iko = ika + ikb
      if(iko.gt.nkmax) stop 4710
c
      if(isga.gt.0)then
         if(isgb.gt.0)then
            isgo = 1
         else
            isgo =-1
         endif
      else
         if(isgb.gt.0)then
            isgo =-1
         else
            isgo = 1
         endif
      endif
c
      iko1 = iko - 1
      ipr = 0
c
      do 200 i = 1, iko1
         ipt = ipr
         k = i
         do 180 j = 1, ikb
            if(k .lt. 1) go to 190
            if(k .gt. ika) go to 150
            ipt = ipt + ia(k)*ib(j)
  150       continue
            k = k - 1
  180    continue
  190    continue
         ipr = ipt/mhalf
         io(i) = ipt - ipr*mhalf
  200 continue
c
      io(iko) = ipr
      if(ipr.ge.mhalf) stop 4720
c
      iko1 = iko
      do 300 i = iko1, ika+1, -1
         if(io(i) .ne. 0) go to 400
         iko = iko - 1
  300 continue
  400 continue
c
      return
      end
*MULDIF
c
c     subroutine muldif
c
c     to perform a multiple precision integer subtraction
c     (for subtracting a decomposed product from another)
c
c     io = ia - ib
c
c     ia represents a decomposed regular length integer or the 
c        decomposed product of two or more regular length integers
c     ib is similarly described
c     io is a decomposed integer which represents ia - ib
c
      subroutine muldif(ia, ib, io, isga, isgb, isgo, ika, ikb, iko,
     *                    nkmax, mhalf)
c
      integer ia(*), ib(*), io(*)
      integer isga, isgb, isgo, ika, ikb, iko, nkmax, mhalf
      integer i, iko1, irel
c
      if(isgb.eq.0)then
         if(isga.eq.0)then
            isgo=0
            iko = 2
            io(1) = 0
            io(2) = 0
            return
         endif
         isgo = isga
         iko = ika
         do 100 i=1,iko
            io(i) = ia(i)
  100    continue
      elseif(isga.eq.0)then
         isgo =-isgb
         iko = ikb
         do 200 i=1,iko
            io(i) = ib(i)
  200    continue
      else
         iko = ika
         if(ikb.lt.ika) then
            do 300 i=ikb+1,ika
               ib(i) = 0
  300       continue
         elseif(ika.lt.ikb) then
            iko = ikb
            do 400 i=ika+1,ikb
               ia(i) = 0
  400       continue
         endif
         if(isga*isgb.gt.0)then
            irel = 0
            do 500 i = iko, 1, -1
               if(ia(i).gt.ib(i))then
                  irel = 1
                  go to 600
               elseif(ia(i).lt.ib(i))then
                  irel = -1
                  go to 600
               endif
  500       continue
  600       continue
            if(irel.eq.0)then
               isgo = 0
               do 700 i=1,iko
                  io(i) = 0
  700          continue
            else
               isgo=isga*irel
               io(1) = (ia(1)-ib(1))*irel
               do 800 i=2,iko
                  if(io(i-1).lt.0) then
                     io(i) =-1
                     io(i-1) = io(i-1) + mhalf
                  else
                     io(i) = 0
                  endif
                  io(i) = io(i) + (ia(i)-ib(i))*irel
  800          continue
               if(io(iko).lt.0) stop 4810
            endif
         else
            isgo=isga
            io(1) = ia(1)+ib(1)
            do 900 i=2,iko
               if(io(i-1).ge.mhalf) then
                  io(i) = 1
                  io(i-1) = io(i-1) - mhalf
               else
                  io(i) = 0
               endif
               io(i) = io(i) + ia(i)+ib(i)
  900       continue
            if(io(iko).ge.mhalf) then
               iko = iko+1
               if(iko.gt.nkmax) stop 4820
               io(iko) = 1
               io(iko-1) = io(iko-1) - mhalf
            endif
         endif
      endif
c
      if(iko .eq. 2) go to 1400
      iko1 = iko
      do 1300 i = iko1, 3, -1
         if(io(i) .ne. 0) go to 1400
         iko = iko - 1
 1300 continue
 1400 continue
c
      return
      end
*XYZNUM
c
c     Routine for determining numerator used in formulas for
c     computing x-,y-z-coordinates of orthogonal center of
c     tetrahedron with vertices ifir, isec, ithi, ifou.
c
      subroutine xyznum(x, y, z, w, x2, y2, z2, w2, ifir, isec, ithi,
     *                  ifou, mhalf, mfull, isclp, isclw, isclr,
     *                  delaun, inx, iny, inz, isgnx, isgny, isgnz,
     *                  iknx, ikny, iknz)
c
      integer x(*), y(*), z(*), w(*), x2(*), y2(*), z2(*), w2(*)
      integer inx(*), iny(*), inz(*)
      integer isgnx, isgny, isgnz, iknx, ikny, iknz
      integer ifir, isec, ithi, ifou, mhalf, mfull, nkmax
      integer isclp(*), isclw(*), isclr(*)
      parameter (nkmax = 30)
      integer io(nkmax), iu(nkmax)
      integer iq2(nkmax), iq3(nkmax), iq4(nkmax)
      integer ix2(nkmax), iy2(nkmax), iz2(nkmax)
      integer ix3(nkmax), iy3(nkmax), iz3(nkmax)
      integer ix4(nkmax), iy4(nkmax), iz4(nkmax)
      integer ixf(nkmax), iyf(nkmax), izf(nkmax)
      integer ixf2(nkmax), iyf2(nkmax), izf2(nkmax)
      integer iwf(nkmax), iw2(nkmax), iw3(nkmax), iw4(nkmax)
      logical delaun
      integer iwfuw, iwsew, iwthw, iwfow
      integer ixfuw, iyfuw, izfuw, ixsew, iysew, izsew
      integer ixthw, iythw, izthw, ixfow, iyfow, izfow
      integer iwfu2, iwse2, iwth2, iwfo2
      integer ixfu2, iyfu2, izfu2, ixse2, iyse2, izse2
      integer ixth2, iyth2, izth2, ixfo2, iyfo2, izfo2
      integer isgw2, isgw3, isgw4, ikw2, ikw3, ikw4
      integer isgq2, isgq3, isgq4, ikq2, ikq3, ikq4
      integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf
      integer isgxf2, isgyf2, isgzf2, ikxf2, ikyf2, ikzf2
      integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2
      integer isgx3, isgy3, isgz3, ikx3, iky3, ikz3
      integer isgx4, isgy4, isgz4, ikx4, iky4, ikz4
      integer isgo, isgu, iko, iku
      integer isgwf, isgcl, ikwf, ikcl
c
      if(delaun) then
         isgw2 = 0
         isgw3 = 0
         isgw4 = 0
      else
         iwfuw = w(ifir)
         iwsew = w(isec)
         iwthw = w(ithi)
         iwfow = w(ifou)
c
         iwfu2 = w2(ifir)
         iwse2 = w2(isec)
         iwth2 = w2(ithi)
         iwfo2 = w2(ifou)
c
         call decmp2(iwf,isgwf,ikwf, iwfuw,iwfu2, mhalf, mfull, isclw)
         isgcl = 1
         ikcl = 2
         call decmp2(io, isgo, iko, iwsew, iwse2, mhalf, mfull, isclw)
         call muldif(io, iwf, iu, isgo, isgwf, isgu, iko, ikwf, iku,
     *                 nkmax, mhalf)
         call mulmul(iu, isclr, iw2, isgu, isgcl, isgw2, iku, ikcl,
     *                 ikw2, nkmax, mhalf)
         call decmp2(io, isgo, iko, iwthw, iwth2, mhalf, mfull, isclw)
         call muldif(io, iwf, iu, isgo, isgwf, isgu, iko, ikwf, iku,
     *                 nkmax, mhalf)
         call mulmul(iu, isclr, iw3, isgu, isgcl, isgw3, iku, ikcl,
     *                 ikw3, nkmax, mhalf)
         call decmp2(io, isgo, iko, iwfow, iwfo2, mhalf, mfull, isclw)
         call muldif(io, iwf, iu, isgo, isgwf, isgu, iko, ikwf, iku,
     *                 nkmax, mhalf)
         call mulmul(iu, isclr, iw4, isgu, isgcl, isgw4, iku, ikcl,
     *                 ikw4, nkmax, mhalf)
      endif
c
      ixfuw = x(ifir)
      iyfuw = y(ifir)
      izfuw = z(ifir)
      ixsew = x(isec)
      iysew = y(isec)
      izsew = z(isec)
      ixthw = x(ithi)
      iythw = y(ithi)
      izthw = z(ithi)
      ixfow = x(ifou)
      iyfow = y(ifou)
      izfow = z(ifou)
c
      ixfu2 = x2(ifir)
      iyfu2 = y2(ifir)
      izfu2 = z2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
      izse2 = z2(isec)
      ixth2 = x2(ithi)
      iyth2 = y2(ithi)
      izth2 = z2(ithi)
      ixfo2 = x2(ifou)
      iyfo2 = y2(ifou)
      izfo2 = z2(ifou)
c
      call decmp2(ixf, isgxf, ikxf, ixfuw, ixfu2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfuw, iyfu2, mhalf, mfull, isclp)
      call decmp2(izf, isgzf, ikzf, izfuw, izfu2, mhalf, mfull, isclp)
      call mulmul(ixf, ixf, ixf2, isgxf, isgxf, isgxf2, ikxf, ikxf,
     *              ikxf2, nkmax, mhalf)
      call mulmul(iyf, iyf, iyf2, isgyf, isgyf, isgyf2, ikyf, ikyf,
     *              ikyf2, nkmax, mhalf)
      call mulmul(izf, izf, izf2, isgzf, isgzf, isgzf2, ikzf, ikzf,
     *              ikzf2, nkmax, mhalf)
      if(isgxf2.lt.0 .or. isgyf2.lt.0 .or. isgzf2.lt.0) stop 5105
c
      call frterm(ixsew, iysew, izsew, ixf, iyf, izf, isgxf, isgyf,
     *            isgzf, ikxf, ikyf, ikzf, ixf2, iyf2, izf2,
     *            isgxf2, isgyf2, isgzf2, ikxf2,
     *            ikyf2, ikzf2, iw2, ix2, iy2, iz2, iq2, isgw2, isgx2,
     *            isgy2, isgz2, isgq2, ikw2, ikx2, iky2, ikz2, ikq2,
     *            mhalf, mfull, ixse2, iyse2, izse2, isclp)
c
      call frterm(ixthw, iythw, izthw, ixf, iyf, izf, isgxf, isgyf,
     *            isgzf, ikxf, ikyf, ikzf, ixf2, iyf2, izf2,
     *            isgxf2, isgyf2, isgzf2, ikxf2,
     *            ikyf2, ikzf2, iw3, ix3, iy3, iz3, iq3, isgw3, isgx3,
     *            isgy3, isgz3, isgq3, ikw3, ikx3, iky3, ikz3, ikq3,
     *            mhalf, mfull, ixth2, iyth2, izth2, isclp)
c
      call frterm(ixfow, iyfow, izfow, ixf, iyf, izf, isgxf, isgyf,
     *            isgzf, ikxf, ikyf, ikzf, ixf2, iyf2, izf2,
     *            isgxf2, isgyf2, isgzf2, ikxf2,
     *            ikyf2, ikzf2, iw4, ix4, iy4, iz4, iq4, isgw4, isgx4,
     *            isgy4, isgz4, isgq4, ikw4, ikx4, iky4, ikz4, ikq4,
     *            mhalf, mfull, ixfo2, iyfo2, izfo2, isclp)
c
      call detrm3(iq2, iy2, iz2, isgq2, isgy2, isgz2,
     *            iq3, iy3, iz3, isgq3, isgy3, isgz3,
     *            iq4, iy4, iz4, isgq4, isgy4, isgz4,
     *            ikq2, ikq3, ikq4, iky2, iky3, iky4,
     *            ikz2, ikz3, ikz4, inx, isgnx, iknx, mhalf)
c
      call detrm3(iq2, iz2, ix2, isgq2, isgz2, isgx2,
     *            iq3, iz3, ix3, isgq3, isgz3, isgx3,
     *            iq4, iz4, ix4, isgq4, isgz4, isgx4,
     *            ikq2, ikq3, ikq4, ikz2, ikz3, ikz4,
     *            ikx2, ikx3, ikx4, iny, isgny, ikny, mhalf)
c
      call detrm3(iq2, ix2, iy2, isgq2, isgx2, isgy2,
     *            iq3, ix3, iy3, isgq3, isgx3, isgy3,
     *            iq4, ix4, iy4, isgq4, isgx4, isgy4,
     *            ikq2, ikq3, ikq4, ikx2, ikx3, ikx4,
     *            iky2, iky3, iky4, inz, isgnz, iknz, mhalf)
c
      return
      end
*FRTERM
c
      subroutine frterm(ixsew, iysew, izsew, ixf, iyf, izf, isgxf,
     *                  isgyf, isgzf, ikxf, ikyf, ikzf, ixf2, iyf2,
     *                  izf2, isgxf2, isgyf2, isgzf2,
     *                  ikxf2, ikyf2, ikzf2, iw2, ix2, iy2, iz2,
     *                  iq2, isgw2, isgx2, isgy2, isgz2, isgq2, ikw2,
     *                  ikx2, iky2, ikz2, ikq2, mhalf, mfull,
     *                  ixse2, iyse2, izse2, isclp)
c
      integer ixsew, iysew, izsew, isgxf, isgyf, isgzf
      integer isgxf2, isgyf2, isgzf2
      integer ikxf, ikyf, ikzf, ikxf2, ikyf2, ikzf2
      integer isgw2, isgx2, isgy2, isgz2, isgq2
      integer ikw2, ikx2, iky2, ikz2, ikq2, mhalf, mfull
      integer ixse2, iyse2, izse2, isclp(*)
      integer isgo, isgu, isgv, isgp, iko, iku, ikv, ikp
      integer ixf(*), iyf(*), izf(*), ixf2(*), iyf2(*), izf2(*)
      integer iw2(*), ix2(*), iy2(*), iz2(*), iq2(*)
      integer nkmax
      parameter (nkmax = 30)
      integer io(nkmax), iu(nkmax), iv(nkmax), ip(nkmax)
c
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix2, isgo, isgxf, isgx2, iko, ikxf, ikx2,
     *              nkmax, mhalf)
      call mulmul(io, io, iu, isgo, isgo, isgu, iko, iko, iku,
     *              nkmax, mhalf)
      call muldif(iu, ixf2, iv, isgu, isgxf2, isgv, iku, ikxf2, ikv,
     *              nkmax, mhalf)
      call muldif(iv, iw2, ip, isgv, isgw2, isgp, ikv, ikw2, ikp,
     *              nkmax, mhalf)
c
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy2, isgo, isgyf, isgy2, iko, ikyf, iky2,
     *              nkmax, mhalf)
      call mulmul(io, io, iu, isgo, isgo, isgu, iko, iko, iku,
     *              nkmax, mhalf)
      call muldif(iu, iyf2, iv, isgu, isgyf2, isgv, iku, ikyf2, ikv,
     *              nkmax, mhalf)
      isgv=-isgv
      call muldif(ip, iv, iu, isgp, isgv, isgu, ikp, ikv, iku,
     *              nkmax, mhalf)
c
      call decmp2(io, isgo, iko, izsew, izse2, mhalf, mfull, isclp)
      call muldif(io, izf, iz2, isgo, isgzf, isgz2, iko, ikzf, ikz2,
     *              nkmax, mhalf)
      call mulmul(io, io, iv, isgo, isgo, isgv, iko, iko, ikv,
     *              nkmax, mhalf)
      call muldif(iv, izf2, ip, isgv, isgzf2, isgp, ikv, ikzf2, ikp,
     *              nkmax, mhalf)
      isgp=-isgp
      call muldif(iu, ip, iq2, isgu, isgp, isgq2, iku, ikp, ikq2,
     *              nkmax, mhalf)
c
      return
      end
*DETRM3
c
      subroutine detrm3(ix2, iy2, iz2, isgx2, isgy2, isgz2,
     *                  ix3, iy3, iz3, isgx3, isgy3, isgz3,
     *                  ix4, iy4, iz4, isgx4, isgy4, isgz4,
     *                  ikx2, ikx3, ikx4, iky2, iky3, iky4,
     *                  ikz2, ikz3, ikz4, io, isgo, iko, mhalf)
c
      integer nkmax
      parameter (nkmax = 30)
      integer io(*), iu(nkmax), iv(nkmax), iw(nkmax)
      integer ix2(*), iy2(*), iz2(*), ix3(*), iy3(*), iz3(*)
      integer ix4(*), iy4(*), iz4(*)
      integer isgx2, isgy2, isgz2, isgx3, isgy3, isgz3
      integer isgx4, isgy4, isgz4, ikx2, ikx3, ikx4, iky2, iky3, iky4
      integer ikz2, ikz3, ikz4, isgo, iko, mhalf
      integer isgu, isgv, isgw, iku, ikv, ikw
c
      call mulmul(ix3, iy4, iv, isgx3, isgy4, isgv, ikx3, iky4, ikv,
     *              nkmax, mhalf)
      call mulmul(ix4, iy3, iu, isgx4, isgy3, isgu, ikx4, iky3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, iw, isgv, isgu, isgw, ikv, iku, ikw,
     *              nkmax, mhalf)
      call mulmul(iw, iz2, io, isgw, isgz2, isgo, ikw, ikz2, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix2, iy4, iv, isgx2, isgy4, isgv, ikx2, iky4, ikv,
     *              nkmax, mhalf)
      call mulmul(ix4, iy2, iu, isgx4, isgy2, isgu, ikx4, iky2, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, iw, isgv, isgu, isgw, ikv, iku, ikw,
     *              nkmax, mhalf)
      call mulmul(iw, iz3, iu, isgw, isgz3, isgu, ikw, ikz3, iku,
     *              nkmax, mhalf)
      call muldif(io, iu, iw, isgo, isgu, isgw, iko, iku, ikw,
     *              nkmax, mhalf)
c
      call mulmul(ix3, iy2, iv, isgx3, isgy2, isgv, ikx3, iky2, ikv,
     *              nkmax, mhalf)
      call mulmul(ix2, iy3, iu, isgx2, isgy3, isgu, ikx2, iky3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
      call mulmul(io, iz4, iu, isgo, isgz4, isgu, iko, ikz4, iku,
     *              nkmax, mhalf)
      call muldif(iw, iu, io, isgw, isgu, isgo, ikw, iku, iko,
     *              nkmax, mhalf)
c
      return
      end
*DOUBNUM
c
      subroutine doubnm(io, isgo, iko, r215, dnum)
c
      integer io(*)
      double precision dnum, r215, rpwr
      integer isgo, iko, i
c
      if(isgo.eq.0) then
         dnum = dble(0)
         go to 1000
      else
         if(iko .lt. 2) stop 6010
         if(iko .gt. 68) stop 6020
         rpwr = dble(1)
         dnum = dble(io(1))
         do 100 i = 2, iko
            rpwr = rpwr*r215
            dnum = dnum + dble(io(i))*rpwr
  100    continue
      endif
      if(isgo.lt.0) dnum = -dnum
c
 1000 continue
      return
      end
*CROSSP
c
c     Routine for determining cross product of two vectors <ifir,isec>
c     and <ifir,ithi>.
c
      subroutine crossp(x, y, z, x2, y2, z2, ifir, isec, ithi,
     *                  mhalf, mfull, isclp, iox, isgox, ikox,
     *                  ioy, isgoy, ikoy, ioz, isgoz, ikoz)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*)
      integer iox(*),ioy(*), ioz(*)
      integer ifir, isec, ithi
      integer isclp(*), mhalf, mfull, nkmax
      parameter (nkmax = 30)
      integer io(nkmax), iu(nkmax), iv(nkmax)
      integer ix2(nkmax), iy2(nkmax), iz2(nkmax)
      integer ix3(nkmax), iy3(nkmax), iz3(nkmax)
      integer ixf(nkmax), iyf(nkmax), izf(nkmax)
      integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew
      integer ixthw, iythw, izthw
      integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2
      integer ixth2, iyth2, izth2
      integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf
      integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2
      integer isgx3, isgy3, isgz3, ikx3, iky3, ikz3
      integer isgox, ikox, isgoy, ikoy, isgoz, ikoz
      integer isgo, isgu, isgv, iko, iku, ikv
c
      ixfiw = x(ifir)
      iyfiw = y(ifir)
      izfiw = z(ifir)
      ixsew = x(isec)
      iysew = y(isec)
      izsew = z(isec)
      ixthw = x(ithi)
      iythw = y(ithi)
      izthw = z(ithi)
c
      ixfi2 = x2(ifir)
      iyfi2 = y2(ifir)
      izfi2 = z2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
      izse2 = z2(isec)
      ixth2 = x2(ithi)
      iyth2 = y2(ithi)
      izth2 = z2(ithi)
c
      call decmp2(ixf, isgxf, ikxf, ixfiw, ixfi2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfiw, iyfi2, mhalf, mfull, isclp)
      call decmp2(izf, isgzf, ikzf, izfiw, izfi2, mhalf, mfull, isclp)
c
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix2, isgo, isgxf, isgx2, iko, ikxf, ikx2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy2, isgo, isgyf, isgy2, iko, ikyf, iky2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izsew, izse2, mhalf, mfull, isclp)
      call muldif(io, izf, iz2, isgo, isgzf, isgz2, iko, ikzf, ikz2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, ixthw, ixth2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix3, isgo, isgxf, isgx3, iko, ikxf, ikx3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iythw, iyth2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy3, isgo, isgyf, isgy3, iko, ikyf, iky3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izthw, izth2, mhalf, mfull, isclp)
      call muldif(io, izf, iz3, isgo, isgzf, isgz3, iko, ikzf, ikz3,
     *              nkmax, mhalf)
c
      call mulmul(iy2, iz3, iv, isgy2, isgz3, isgv, iky2, ikz3, ikv,
     *              nkmax, mhalf)
      call mulmul(iz2, iy3, iu, isgz2, isgy3, isgu, ikz2, iky3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, iox, isgv, isgu, isgox, ikv, iku, ikox,
     *              nkmax, mhalf)
c
      call mulmul(iz2, ix3, iv, isgz2, isgx3, isgv, ikz2, ikx3, ikv,
     *              nkmax, mhalf)
      call mulmul(ix2, iz3, iu, isgx2, isgz3, isgu, ikx2, ikz3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, ioy, isgv, isgu, isgoy, ikv, iku, ikoy,
     *              nkmax, mhalf)
c
      call mulmul(ix2, iy3, iv, isgx2, isgy3, isgv, ikx2, iky3, ikv,
     *              nkmax, mhalf)
      call mulmul(iy2, ix3, iu, isgy2, isgx3, isgu, iky2, ikx3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, ioz, isgv, isgu, isgoz, ikv, iku, ikoz,
     *              nkmax, mhalf)
c
      return
      end
