c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c
c     Program volare
c
c     This program will compute volumes of Power/Voronoi cells in the
c     Power/Voronoi diagram of a set of points in 3-dimensional space
c     and areas of facets of these cells.
c
c     It is assumed that a Regular/Delaunay tetrahedralization for
c     the set has already been computed using existing program regtet
c     and that the output tetrahedron list produced by regtet contains
c     only real tetrahedra (no artificial ones), i. e. logical variable
c     artfcl was set to .false. during the execution of regtet. It is
c     also assumed that the Power/Voronoi vertices and unbounded edges
c     of the Power/Voronoi diagram of the set have already been
c     computed using existing program pwrvtx with the output
c     tetrahedron list from regtet as input. The output Power/Voronoi
c     vertex and tetrahedron lists from pwrvtx are then used as input
c     below. It is noted that because of degeneracies there may be
c     facets that are not 2-dimensional any one of which will have zero
c     area if bounded. In addition, if weights are present, there
c     may be Power cells that are not 3-dimensional any one of which
c     will have zero volume if bounded.
c
c     Subroutine volare 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 volare.
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, nfmax
      parameter (nmax=150000, nvmax= 7*nmax, nfmax = 8*nmax)
c
      double precision xp(nvmax), yp(nvmax), zp(nvmax)
      integer ix(nvmax), iy(nvmax), iz(nvmax)
      integer ix2(nvmax), iy2(nvmax), iz2(nvmax)
      integer icon(8,nvmax), ifl(nvmax), id(nvmax)
      integer is(nmax), itl(nmax), ifn(nmax), ifc(nfmax)
      double precision vol(nmax), area(nfmax)
      integer nw, nt, nq, nb, icfig
c
      integer idumy, nw2, nq2, i, j
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'This program is for computing volumes of ',
     *          'Power/Voronoi cells'
      write(*,*)'in the Power/Voronoi diagram of a set of ',
     *          'points in 3-d space'
      write(*,*)'and for computing areas of facets ',
     *          'of these cells.'
      write(*,*)' '
      write(*,*)'It is assumed that a Regular/Delaunay ',
     *          'tetrahedralization of the'
      write(*,*)'set has already been computed using existing program ',
     *          'regtet and'
      write(*,*)'that the output tetrahedron list produced by regtet ',
     *          'contains only'
      write(*,*)'real tetrahedra (logical variable artfcl was set to ',
     *          '.false. during'
      write(*,*)'the execution of regtet).'
      write(*,*)' '
      write(*,*)'It is also assumed that the Power/Voronoi vertices ',
     *          'and the unbounded'
      write(*,*)'edges of the Power/Voronoi diagram of the set have ',
     *          'already been'
      write(*,*)'computed using existing program pwrvtx with the ',
     *          'output tetrahedron'
      write(*,*)'list from regtet as input.'
      write(*,*)' '
      write(*,*)'The output Power/Voronoi vertex and tetrahedron ',
     *          'lists from pwrvtx'
      write(*,*)'are then used as input by this program. Note that ',
     *          'because of'
      write(*,*)'degeneracies there may facets that are not ',
     *          '2-dimensional any one'
      write(*,*)'of which will have zero area if bounded. ',
     *          'In addition, if weights'
      write(*,*)'are present, there may be Power cells that are ',
     *          'not 3-dimensional'
      write(*,*)'any one of which will have zero volume if bounded.'
      write(*,*)' '
      write(*,*)'Arrays ifn, ifc, vol, area  will contain output ',
     *          'information.'
      write(*,*)'Description of these arrays appears in driver ',
     *          'routine volare.'
      write(*,*)' '
      write(*,*)'In the curret main program these arrays are saved ',
     *          'in output file'
      write(*,*)'areas-vols.'
c
c----------------------------------------------------------------------
c
c     open files
c
      open (12, file = 'tetra-unbd')
      open (13, file = 'power-vrts')
      open (14, file = 'areas-vols')
c
c----------------------------------------------------------------------
c
c     read tetrahedralization information
c
      write(*,*)' '
      write(*,*)'Reading tetrahedralization information ...'
      read (12,*) idumy
      read (12,*) nw, nq
      if(nw.gt.nmax .or. nq.gt.nvmax) stop 10
      read (12,90) ((icon(i,j), i = 1, 8), j = 1, nq)
      read (12,95) (is(i), i = 1, nw)
   90 format (8i10)
   95 format (7i10)
c
c----------------------------------------------------------------------
c
      if(idumy.lt.0) stop 20
c
c----------------------------------------------------------------------
c
c     set icfig for power vertices
c
      icfig = 9
c
c----------------------------------------------------------------------
c
c     read power vertices
c
      write(*,*)' '
      write(*,*)'Reading power vertices ...'
      read (13,*) nw2, nt, nq2
      if(nw2.ne.nw .or. nt.ge.nq2 .or. nq2.ne.nq) stop 30
      read (13,*) (ifl(i),i=1,nq)
      do 300 i = 1, nt
         read (13,*) xp(i), yp(i), zp(i)
  300 continue
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'Computation of areas and volumes to begin: driver ',
     *          'subroutine volare'
      write(*,*)'will be called from the main routine and areas and ',
     *          'volumes will be'
      write(*,*)'computed.'
      write(*,*)' '
      write(*,*)'Please wait ...'
      write(*,*)' '
c
c----------------------------------------------------------------------
c 
c     call volare to compute areas and volumes
c
      call volare(xp, yp, zp, ix, iy, iz, ix2, iy2, iz2, icon, ifl, id,
     *            is, itl, ifn, ifc, vol, area, nmax, nvmax, nfmax,
     *            nw, nt, nq, nb, icfig)
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'(Back to the main routine).'
      write(*,*)'Computation of areas and volumes has been completed.'
c
c----------------------------------------------------------------------
c
c     save areas and volumes information
c
      write(*,*)' '
      write(*,*)'Saving area and volume information ...'
      write (14,*) nw, nb
      write (14,600) (ifn(i),i=1,nw)
      write (14,700) (ifc(i),i=1,nb)
      write (14,*) (vol(i),i=1,nw)
      write (14,*) (area(i),i=1,nb)
  600 format (7i10)
  700 format (8i10)
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      stop
      end
*VOLARE
c**********************************************************************
c
c     Driver subroutine of Fortran 77 program VOLARE for computing
c     volumes of Power/Voronoi cells in the Power/Voronoi diagram
c     of a set of points in 3-dimensional space and areas of facets
c     of these cells.
c
c     It is assumed that a Regular/Delaunay tetrahedralization for
c     the set has already been computed using existing program regtet
c     and that the output tetrahedron list produced by regtet contains
c     only real tetrahedra (no artificial ones), i. e. logical variable
c     artfcl was set to .false. during the execution of regtet. It is
c     also assumed that the Power/Voronoi vertices and unbounded edges
c     of the Power/Voronoi diagram of the set have already been
c     computed using existing program pwrvtx with the output
c     tetrahedron list from regtet as input. The output Power/Voronoi
c     vertex and tetrahedron lists from pwrvtx are then used as input
c     by this program. It is noted that because of degeneracies there
c     may be facets that are not 2-dimensional any one of which will
c     have zero area if bounded. In addition, if weights are present,
c     there may be Power cells that are not 3-dimensional any one of
c     which will have zero volume if bounded.
c
c     Arrays ifn, ifc, vol, area will contain output information.
c     For i=1,...,nw, vol(i) will be the volume of the Power/Voronoi
c     cell of point i if point i is active, non-redundant, and its
c     cell is bounded and none of its vertices was marked as being
c     too large by program pwrvtx. Otherwise vol(i) will equal -40.0
c     if point i is not active or was marked as redundant by program
c     regtet; -30.0 if the Power/Voronoi cell of point i is unbounded
c     and has at least one vertex marked as being too large; -20.0 if
c     the cell is unbounded but has no vertices marked as being too
c     large; and -10.0 if the cell is bounded but has at least one
c     vertex marked as being too large.
c     For l=1,...,nb, area(l) will be the area of the lth facet of
c     the Power/Voronoi diagram of the set if the facet is bounded
c     and none of its vertices was marked as being too large by
c     program pwrvtx. Otherwise area(l) will equal -30.0 if the facet
c     is unbounded and has at least one vertex marked as being too
c     large; -20.0 if the facet is unbounded but has no vertices
c     marked as being too large; and -10.0 if the facet is bounded
c     but has at least one vertex marked as being too large. Which
c     facet is the lth facet is determined using arrays ifn and ifc
c     as described below.
c
c     For l=1,...,nb, ifc(l) will equal an integer, 0 < ifc(l) < nw+1.
c     For i=1,...,nw, ifn(i) will equal zero if point i is not active
c     or was marked as redundant by program regtet. Otherwise if l2
c     equals ifn(i) and l1 equals the largest among 0 and ifn(j),
c     j < i, then -1 < l1 < l2 < nb+1, i < ifc(l) for l = l1+1,...,l2, 
c     and point ifc(l) is a Power/Voronoi neighbor of point i for each
c     l, l = l1+1,...,l2. For l = l1+1,...,l2, the facet shared by
c     point i and point ifc(l) is the lth facet of the Power/Voronoi
c     diagram of the set. It is noted that because of degeneracies
c     this facet may be not be 2-dimensional. Also it is noted that
c     nb equals the largest among ifn(i), i = 1,...,nw.
c
c     The idea is to have a way of identifying each facet in terms of
c     the two points that define it, to do this only one time, and then
c     keep a list of the facets identified this way in a way that is
c     systematic. Essentially for each integer i, i = 1, ..., nw, if
c     point i is active and not redundant, two integers l1 and l2 are
c     identified so that the points ifc(l), l = l1+1, ..., l2, are the
c     Voronoi/Power neighbors of point i for which i < ifc(l). The last
c     inequality guarantees that each facet is identified only once.
c     Since point i is active and not redundant, ifn(i) is not zero and
c     has been set in the program so that the desired value of l2 is
c     actually ifn(i). If point (i-1) is active and not redundant, then
c     ifn(i-1) is not zero and the desired value of l1 is then ifn(i-1).
c     If point (i-1) is either not active or redundant then ifn(i-1) is
c     zero and ifn(i-2) then has to be checked in the same way, and so
c     on. This is why l1 is the largest of 0 and ifn(j), j<i. The last
c     inequality means that if i is greater than 1 then l1 will be 0 or
c     the first of ifn(i-1), ifn(i-2), ..., ifn(1) that is not zero. If
c     i equals 1 then l1 is definitely 0.
c    
      subroutine volare(xp, yp, zp, ix, iy, iz, ix2, iy2, iz2, icon,
     *                  ifl, id, is, itl, ifn, ifc, vol, area,
     *                  nmax, nvmax, nfmax, nw, nt, nq, nb, icfig)
c
      double precision xp(*), yp(*), zp(*)
      integer ix(*), iy(*), iz(*)
      integer ix2(*), iy2(*), iz2(*)
      integer icon(8,*), ifl(*), id(*), is(*), itl(*), ifn(*), ifc(*)
      double precision vol(*), area(*)
      integer nmax, nvmax, nfmax, nw, nt, nq, nb, icfig, mhalf, mfull
c
c     initialize Fortran 77 word lengths
c
      mhalf=32768
      mfull=1073741824
c
c     test parameters
c
      if(nw.gt.nmax .or. nt.ge.nq .or. nq.gt.nvmax) stop 110
c
c     write(*,*)' '
c     write(*,*)'Entering vacomp ...'
      call vacomp(xp, yp, zp, ix, iy, iz, ix2, iy2, iz2, icon, ifl, id,
     *            is, itl, ifn, ifc, vol, area, nmax, nvmax, nfmax,
     *            nw, nt, nq, nb, icfig, mhalf, mfull)
c     write(*,*)' '
c     write(*,*)'Leaving vacomp ...'
c
      return
      end
*VACOMP
c
c     This subroutine computes volumes and areas
c
      subroutine vacomp(x, y, z, ix, iy, iz, ix2, iy2, iz2, icon, ifl,
     *                  id, is, itl, ifn, ifc, vol, area, nmax, nvmax,
     *                  nfmax, nw, nt, nq, nb, icsfig, mhalf, mfull)
c
      double precision x(*), y(*), z(*)
      integer ix(*), iy(*), iz(*)
      integer ix2(*), iy2(*), iz2(*)
      integer icon(8,*), ifl(*), id(*), is(*), itl(*), ifn(*), ifc(*)
      double precision vol(*), area(*)
      integer nmax, nvmax, nfmax, nw, nt, nq, nb, icsfig, mhalf, mfull
      double precision r215, deps, dscle, dfull, dfill, decml
      integer i, ng, isclu, isclp(2), isgcl
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
      deps = 0.9d0
      if(dscle.lt.deps) stop 293
      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, nt
         if(ifl(i).eq.0) go to 235
         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     identify Power/Voronoi neighbor pairs
c     compute areas and volumes
c
      ng = 0
      r215 = dble(mhalf)
c
      call nghpair(ix, iy, iz, ix2, iy2, iz2, icon, ifl, is, ifn, ifc,
     *             id, itl, vol, area, nmax, nvmax, nfmax, nw, nt, nq,
     *             nb, ng, mhalf, mfull, isclp, r215, dscle, deps)
c
      return
      end
*NGHPAIR
c
c     subroutine nghpair to -
c
c     obtain neighbor pairs in power diagram
c     compute areas and volumes
c
      subroutine nghpair(x, y, z, x2, y2, z2, icon, ifl, is, ifn, ifc,
     *                   id, itl, vol, area, nmax, nvmax, nfmax,
     *                   nw, nt, nq, nb, ng, mhalf, mfull, isclp,
     *                   r215, dscle, deps)
c
      integer x(*), y(*), z(*)
      integer x2(*), y2(*), z2(*)
      integer icon(8,*), ifl(*), is(*), ifn(*), ifc(*), id(*), itl(*)
      double precision vol(*), area(*)
      integer nkmax
      parameter (nkmax = 30)
      integer io(nkmax), iox(nkmax), ioy(nkmax), ioz(nkmax)
      integer nmax, nvmax, nfmax, nw, nt, nq, nb, ng, mhalf, mfull
      double precision r215, dscle, deps, dnum, dnux, dvol, dare
      double precision xnum, ynum, znum
      integer ifval, i, ileft, ine1, iscur, isini, ilift
      integer ivoli, ivuli, iboli, ibuli
      integer isadj, isbeg, ine, nb1, isclp(*)
      integer isgox, isgoy, isgoz, ikox, ikoy, ikoz
      integer isgo, iko
c
c     initialize
c
      do 30 i=1,nmax
         ifn(i)=0
   30 continue
c
      ifval=0
      do 50 i=1,nvmax
         id(i)=0
   50 continue
c
      do 70 i=1,nmax
         itl(i)=0
   70 continue
c
      do 80 i=1,nmax
         vol(i) = -40.0d0
   80 continue
c
      ileft = ng
      nb = 0
c
c     process next point
c
  100 continue
      ine1 = nb + 1
      ivoli = 1
      ivuli = 1
      ileft = ileft + 1
      if(ileft.gt.nw) go to 2000
      if(ileft.le.(ileft/1000)*1000)
     *   write(*,*)'Number of processed cells =',ileft
      ifval = ifval + 1
      iscur = is(ileft)
      if(iscur.le.0) go to 100
      if(iscur.gt.nt) stop 410
      isini = iscur
c
c     reorder isini so that ileft equals icon(5,isini)
c
      call site1ordr(icon, ileft, isini)
c
      do 300 i=6,8
         ilift=icon(i,isini)
         if(ilift.le.ng) stop 420
         itl(ilift) = ifval
         nb = nb + 1
         if(nb.gt.nfmax) stop 425
         ifc(nb) = ilift
         is(ilift) = isini
  300 continue
c
c     mark current tetrahedron
c
  400 continue
      if(ifl(iscur).eq.0) ivuli = -1
      id(iscur) = ifval
c
c     obtain next tetrahedron with ileft as a vertex
c
      isadj = icon(2,iscur)
      if(isadj.le.0.or.isadj.gt.nq) stop 430
      if(ifl(isadj).eq.2) then
         if(isadj.le.nt) stop 432
         ivoli =-1
         go to 600
      endif
      if(isadj.gt.nt) stop 434
      if(id(isadj).eq.ifval) go to 600
      ilift = icon(8,iscur)
      go to 900
  600 continue
      isadj = icon(3,iscur)
      if(isadj.le.0.or.isadj.gt.nq) stop 440
      if(ifl(isadj).eq.2) then
         if(isadj.le.nt) stop 442
         ivoli =-1
         go to 700
      endif
      if(isadj.gt.nt) stop 444
      if(id(isadj).eq.ifval) go to 700
      ilift = icon(6,iscur)
      go to 900
  700 continue
      isadj = icon(4,iscur)
      if(isadj.le.0.or.isadj.gt.nq) stop 450
      if(iscur .eq. isini) go to 800
      if(isadj.gt.nt) stop 452
      if(icon(3,isadj) .eq. iscur) then
          iscur = isadj
          go to 700
      elseif(icon(2,isadj) .eq. iscur) then
          iscur = isadj
          go to 600
      elseif(icon(4,isadj) .eq. iscur) then
          if(isadj .ne. isini) stop 460
          go to 1000
      else
          stop 465
      endif
  800 continue
      if(ifl(isadj).eq.2) then
         if(isadj.le.nt) stop 467
         ivoli =-1
         go to 1000
      endif
      if(isadj.gt.nt) stop 468
      if(id(isadj).eq.ifval) go to 1000
      ilift = icon(7,iscur)
c
c     reorder isadj so that ileft equals icon(5,isadj) and ilift
c     equals icon(6,isadj)
c
  900 continue
      call reordr(icon, ileft, ilift, isadj)
      iscur = isadj
      ilift=icon(8,iscur)
      if(ilift.le.ng) stop 470
      if(itl(ilift).eq.ifval) go to 400
      itl(ilift) = ifval
      nb = nb + 1
      if(nb.gt.nfmax) stop 475
      ifc(nb) = ilift
      is(ilift) = iscur
      go to 400
c
c     compute areas and volume for current point
c
 1000 continue
      if(ine1.gt.nb) stop 480
      isbeg = is(ifc(ine1))
      dvol = 0.0d0
      do 1200 ine = ine1, nb
         dare = 0.0d0
         iboli = 1
         ibuli = 1
         ilift = ifc(ine)
         isini = is(ilift)
         if(isini.le.0 .or. isini.gt.nt) stop 490
         if(ifl(isini).eq.0) ibuli =-1
         call reordr(icon, ileft, ilift, isini)
         iscur = icon(3,isini)
         if(iscur.le.0 .or. iscur.gt.nq) stop 510
         if(ifl(iscur).eq.2) then
            if(iscur.le.nt) stop 512
            iboli =-1
         endif
         if(ifl(iscur).ne.2 .and. iscur.gt.nt) stop 514
         if(ifl(iscur).eq.0) ibuli =-1
         call reordr(icon, ileft, ilift, iscur)
         isadj = icon(3,iscur)
         if(isadj .eq. isini) stop 520
 1150    continue
         if(isadj.le.0 .or. isadj.gt.nq) stop 530
         if(ifl(isadj).eq.2) then
            if(isadj.le.nt) stop 532
            iboli =-1
         endif
         if(ifl(isadj).ne.2 .and. isadj.gt.nt) stop 534
         if(ifl(isadj).eq.0) ibuli =-1
         call reordr(icon, ileft, ilift, isadj)
c
         if(iboli.ne.-1 .and. ibuli.ne.-1) then
c
c        here compute area of isini-iscur-isadj triangle and
c        add to area(ine)
c
            call crossp(x, y, z, x2, y2, z2, isini, iscur, isadj,
     *                  mhalf, mfull, isclp, iox, isgox, ikox,
     *                  ioy, isgoy, ikoy, ioz, isgoz, ikoz)
            if(isgox.eq.0 .and. isgoy.eq.0 .and. isgoz.eq.0) go to 1160
            call doubnm(iox, isgox, ikox, r215, xnum)
            call doubnm(ioy, isgoy, ikoy, r215, ynum)
            call doubnm(ioz, isgoz, ikoz, r215, znum)
            dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum))
            if(dnux.lt.deps) dnux = 1.0d0
            xnum = xnum/dnux
            ynum = ynum/dnux
            znum = znum/dnux
            dnum = dsqrt(xnum**2+ynum**2+znum**2)
            dare = dare + ((dnux/dscle)/dscle)*dnum
         endif
c
         if(ivoli.ne.-1 .and. ivuli.ne.-1) then
c
c           here compute volume of isbeg-isini-iscur-isadj tetrahedron
c           and add to vol(ileft)
c
            if(iboli.eq.-1 .or. ibuli.eq.-1) stop 540
            call tetvol(x, y, z, x2, y2, z2, isbeg, isini, iscur, isadj,
     *                  mhalf, mfull, isclp, io, isgo, iko)
            if(isgo.le.0) go to 1160
            call doubnm(io, isgo, iko, r215, dnum)
            dvol = dvol + (((dnum/dscle)/dscle)/dscle)
         endif
c
 1160    continue
         iscur = isadj
         isadj = icon(3,iscur)
         if(isadj .ne. isini) go to 1150
         if(iboli.eq.-1 .and. ibuli.eq.-1) then
            area(ine) = -30.0
         elseif(iboli.eq.-1) then
            area(ine) = -20.0
         elseif(ibuli.eq.-1) then
            area(ine) = -10.0
         else
            area(ine) = dare/2.0d0
         endif
 1200 continue
c
      if(ivoli.eq.-1 .and. ivuli.eq.-1) then
         vol(ileft) = -30.0
      elseif(ivoli.eq.-1) then
         vol(ileft) = -20.0
      elseif(ivuli.eq.-1) then
         vol(ileft) = -10.0
      else
         vol(ileft) = dvol/6.0d0
      endif
c
      nb1 = nb
      nb = ine1 - 1
      do 1600 ine = ine1, nb1
         if(ifc(ine).lt.ileft) go to 1600
         nb = nb+1
         ifc(nb) = ifc(ine)
         area(nb) = area(ine)
 1600 continue
c
      if(nb.lt.ine1) go to 100
      ifn(ileft) = nb
      go to 100
 2000 continue
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
*SITE1ORDR
c
c     subroutine site1ordr 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 site1ordr(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
*TETVOL
c
c     Routine for determining volume of tetrahedron with vertices
c     ifir, isec, ithi, ifou times 6.
c
      subroutine tetvol(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
*DOUBNM
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