c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c
c     This program is for locating points in a tetrahedralization by
c     moving in a straight line to each point from a point whose
c     location is known.
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, nlmax
      parameter (nmax=150000, nvmax= 7*nmax, nlmax=1000)
c
      double precision x(nmax), y(nmax), z(nmax)
      double precision xa(nlmax), ya(nlmax), za(nlmax)
      integer ix(nmax), iy(nmax), iz(nmax)
      integer ix2(nmax), iy2(nmax), iz2(nmax)
      integer icon(8,nvmax), ifl(nvmax), is(nmax)
      integer xc(8), yc(8), zc(8)
      integer nv, nw, nt, na, naddl, icfig, iwfig
      integer isclp(2), mhalf, mfull, nara, ibeg, isbeg, iftal
      integer nw0, ityp0, iscu0, sit01, sit02
      integer nw1, itype, iscur, site1, site2
      double precision wlenx, wleny, wlenz, wlenw, epz
      double precision deps, dscle, dfull, dfill
      logical delaun, pntoff, flphis, artfcl
      logical random, reccor, redchk
      logical regtet
c
      double precision xcor, ycor, zcor
      integer ideli, ipnti, iflpi, iarti, irani, ireci, iredi
      integer ico1, ico2, ico3, ico4, ico5, ico6, ico7, ico8
      integer i, j, nl
      character*1 answ
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'This program is for locating points in  ',
     *          'a tetrahedralization'
      write(*,*)'by moving in a straight line to each point ',
     *          'from a point'
      write(*,*)'whose location is known.'
c
   10 format(a1)
c
      write(*,*)' '
      write(*,*)'Do you want program to be run as if another program'
      write(*,*)'called regtet.f has already been run?(y/n)'
      write(*,*)'If that is the case then part of the input data for'
      write(*,*)'this program will be input/output data that is'
      write(*,*)'used/produced by regtet.f.'
      read(5,10) answ
      if(answ.eq.'y'.or. answ.eq.'Y') then
         regtet = .true.
         write(*,*)'Program will be run as if input tetrahedralization'
         write(*,*)'was computed with program regtet.f.'
      else
         regtet = .false.
         write(*,*)'Program will be run using input tetrahedralization'
         write(*,*)'provided by user not necessarily related to a'
         write(*,*)'previous execution of program regtet.f.'
         write(*,*)' '
         write(*,*)'Enter icfig, i. e. number of significant figures ',
     *             'of decimal part'
         write(*,*)'of point coordinates, -1 < icfig < 10, that will ',
     *             'be assumed to be'
         write(*,*)'compatible with input tetrahedralization during ',
     *             'current execution'
         write(*,*)'of program (total number of significant figures ',
     *             'figures should'
         write(*,*)'be at most 14 with at most 9 to either the left ',
     *             'or the right of'
         write(*,*)'the decimal point):'
         read(5,*) icfig
      endif
c
c----------------------------------------------------------------------
c
c     open files
c
c     Input file pnts-wghts contains the vertices of the
c     tetrahedralization, one vertex per line, in terms of their
c     coordinates, 3 numbers per line: the x-, y-, z-coordinates
c     of the vertex. Accordingly the number of lines in this file
c     equals the number of vertices, and the index of a vertex is
c     defined as the line number in this file that contains its
c     coordinates.
c
c     Input file tetahedra contains information about the
c     tetrahedalization in which points are to be located. If the
c     program is being run as if program regtet.f had been previously
c     executed then tetrahedra is just the output file by that name
c     obtained from the execution of regtet.f. Otherwise tetrahedra
c     contains information about the tetrahedralization, one
c     tetrahedron per line, in such a way that it can be read into
c     array icon with icon as described in program regtet.f. When the
c     execution of this program is not based on a previous execution
c     of regtet.f then the number of lines in this file equals the
c     number of tetrahedra in the tetrahedralization. In either case
c     given positive integer l less than or equal to the number of
c     tetrahedra in the tetrahedralization, once the data in input
c     file tetrahedra is read, icon(i,l), i = 1,...,8, will contain
c     information about neighbors and vertices of the lth tetrahedron
c     in the tetrahedralization. Again, a description of array icon
c     can be found in documention in program regtet.f.
c
c     Input file locatepnts contains the points to be located in the
c     tetrahedralization, one point per line, in terms of their
c     coordinates, 3 numbers per line: the x-, y-, z-coordinates
c     of the point.
c
c     Ouput file locatinfor contains information about the location
c     of the points in input file locatepnts, four integers per line:
c     itype, iscur, site1, site2. The number of lines in this file
c     equals the number of lines in file locatepnts. Accordingly
c     the information in a line of this file, say line l, corresponds
c     to the point in line l of file locatepnts. See below for
c     a description of the integers itype, iscur, site1, site2.
c     
      open (11, file = 'pnts-wghts')
      open (12, file = 'tetrahedra')
      open (13, file = 'locatepnts')
      open (14, file = 'locatinfor')
c
c----------------------------------------------------------------------
c
c     set tolerance
c
      epz = 0.01d0
c
c----------------------------------------------------------------------
c
c     read tetrahedralization information
c
      write(*,*)' '
      write(*,*)'Reading tetrahedralization information ...'
      if(.not.regtet) go to 100
      read (12,40) 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,50) ((icon(i,j), i = 1, 8), j = 1, nt)
      read (12,55) (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
   40 format (7(1x,i1))
   50 format (8i10)
   55 format (7i10)
      go to 110
c
  100 continue
      nt = 0
  105 continue
      read (12, *, end = 120) ico1, ico2, ico3, ico4,
     *                        ico5, ico6, ico7, ico8
      nt = nt + 1
      if (nt .gt. nvmax) stop 20
      icon(1,nt) = ico1
      icon(2,nt) = ico2
      icon(3,nt) = ico3
      icon(4,nt) = ico4
      icon(5,nt) = ico5
      icon(6,nt) = ico6
      icon(7,nt) = ico7
      icon(8,nt) = ico8
      go to 105
c
c----------------------------------------------------------------------
c
c     code for avoiding certain warning messages during compilation
c
  110 continue
      if(pntoff) ipnti = 1
      if(flphis) iflpi = 1
      if(random) irani = 1
      if(redchk) iredi = 1
      if(.not.delaun) then
         if(iwfig.lt.0 .or. iwfig.gt.9) stop 30
      endif
      if(reccor .and. .not.delaun) then
         wlenw = wlenw + 0.0d0
      endif
c
c----------------------------------------------------------------------
c
c     read vertex data
c
  120 continue
      write(*,*)' '
      write(*,*)'Reading vertex data ...'
      nv = 0
  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
      go to 130
  140 continue
c
c     read points to be located
c
      na = 0
  150 continue
      read (13, *, end = 160) xcor, ycor, zcor
      na = na + 1
      if (na .gt. nlmax) stop 50
      xa(na) = xcor
      ya(na) = ycor
      za(na) = zcor
      go to 150
  160 continue
      if(na.lt.1) stop 60
      if(.not.regtet) go to 200
c
c----------------------------------------------------------------------
c
c     call complt to complete and process vertex data structure
c
      call complt(x, y, z, xc, yc, zc, nv, nw, nt, nmax, nvmax,
     *            wlenx, wleny, wlenz, naddl, artfcl, reccor)
c
c----------------------------------------------------------------------
c
c     transform double precision coordinates into 2-integer
c     decomposition, find initial vertex, etc.
c
  200 continue
      call transf(x, y, z, ix, iy, iz, ix2, iy2, iz2, icon, is, ifl,
     *            nv, nw, nt, nmax, nvmax, icfig, nara, nw0, nw1,
     *            ibeg, isbeg, epz, mhalf, mfull, isclp, iftal,
     *            deps, dscle, dfull, dfill, artfcl)
      x(nw1) = x(ibeg)
      y(nw1) = y(ibeg)
      z(nw1) = z(ibeg)
      ix(nw1) = ix(ibeg)
      iy(nw1) = iy(ibeg)
      iz(nw1) = iz(ibeg)
      ix2(nw1) = ix2(ibeg)
      iy2(nw1) = iy2(ibeg)
      iz2(nw1) = iz2(ibeg)
      itype = 1
      iscur = isbeg
      site1 = ibeg
      site2 = 0
      nl = 0
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'Point location process to begin '
      write(*,*)' '
      write(*,*)'Please wait ...'
c
c----------------------------------------------------------------------
c 
c     call pntloc to locate current point
c
c
c     On input the following conventions are used when calling pntloc:
c
c     ityp0 = 1 then nw0 is vertex sit01 of iscu0
c     ityp0 = 2 then nw0 is in the interior of tetrahedron iscu0
c     ityp0 = 3 then nw0 is in edge of iscu0 with vertices sit01
c               and sit02
c     ityp0 = 4 then nw0 is in facet of iscu0 opposite vertex
c               sit01 of iscu0
c
c     On output the following conventions are used when calling pntloc:
c
c     itype = 1 then nw1 is vertex site1 of iscur
c     itype = 2 then nw1 is in the interior of tetrahedron iscur
c     itype = 3 then nw1 is in edge of iscur with vertices site1
c               and site2
c     itype = 4 then nw1 is in facet of iscur opposite vertex
c               site1 of iscur
c     itype =-1 then nw1 is outside union of tetrahedra
c
  300 continue
      nl = nl+1
      if(itype.eq.-1) go to 400
      x(nw0) = x(nw1)
      y(nw0) = y(nw1)
      z(nw0) = z(nw1)
      ix(nw0) = ix(nw1)
      iy(nw0) = iy(nw1)
      iz(nw0) = iz(nw1)
      ix2(nw0) = ix2(nw1)
      iy2(nw0) = iy2(nw1)
      iz2(nw0) = iz2(nw1)
      ityp0 = itype
      iscu0 = iscur
      sit01 = site1
      sit02 = site2
  400 continue
      xcor = xa(nl)
      ycor = ya(nl)
      zcor = za(nl)
      call pntloc(x, y, z, ix, iy, iz, ix2, iy2, iz2, icon, ifl, is,
     *            nt, nw0, ityp0, iscu0, sit01, sit02, nw1, itype,
     *            iscur, site1, site2, xc, yc, zc, epz, mhalf, mfull,
     *            isclp, iftal, dscle, dfull, dfill, xcor, ycor, zcor,
     *            nara)
c
      write(14,*) itype, iscur, site1, site2
c
      if(nl.lt.na) go to 300
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      write(*,*)'Point location process has been completed'
c
c----------------------------------------------------------------------
c
      write(*,*)' '
      stop
      end
*COMPLT
c**********************************************************************
c
c     Subroutine for completing and processing vertex data structure.
c
      subroutine complt(x, y, z, xc, yc, zc,  nv, nw, nt, nmax, nvmax,
     *                  wlenx, wleny, wlenz, naddl, artfcl, reccor)
c
      integer nmax, nvmax
      double precision x(*), y(*), z(*)
      integer xc(*), yc(*), zc(*)
      integer nv, nw, nt, irec, naddl
      double precision wlenx, wleny, wlenz
      logical artfcl, reccor
c
      double precision xd(8), yd(8), zd(8)
      double precision xmax, xmin, ymax, ymin, zmax, zmin
      double precision xlen, ylen, zlen, wlen, xctr, yctr, zctr
      double precision xmx, ymx, zmx, xmn, ymn, zmn
      double precision xlan, ylan, zlan, wlan, dfull
      double precision xint, yint, zint, xcor, ycor, zcor
      integer no, irec1, nu, i, ng, naddm, j, mfull
c
c     define mfull = 2^30
c
      mfull=1073741824
c
c     test parameters
c
      if(nv.lt.1 .or. nv.gt.nmax) stop 100
      if(nt.lt.1 .or. nt.gt.nvmax) stop 110
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 120
         if(naddl.lt.2) stop 130
      else
         wlenx = 0.0d0
         wleny = 0.0d0
         wlenz = 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)
      do 20 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)
   20 continue
c
c     shift data
c
      irec = 0
      if(artfcl) irec = 8
      if(reccor) irec = irec + 6*(naddl**2) - 12*naddl + 8
      if(irec.eq.0) go to 40
      irec1 = irec + 1
      nv = nv + irec
      if(nv .gt. nmax) stop 150
      do 30 no = nv, irec1, -1
         nu = no - irec1 + 1
         x(no) = x(nu)
         y(no) = y(nu)
         z(no) = z(nu)
   30 continue
   40 continue
      if(nv.ne.nw) stop 160
c
c     construct cube
c
      if(.not.artfcl) go to 85
      xlen=xmax-xmin
      ylen=ymax-ymin
      zlen=zmax-zmin
      wlen=(dmax1(xlen,ylen,zlen)/2.0d0) + dmax1(wlenx,wleny,wlenz)
     *     + 4.0d0
c
      xctr=(xmax+xmin)/2.0d0
      yctr=(ymax+ymin)/2.0d0
      zctr=(zmax+zmin)/2.0d0
c     WRITE(*,*)'XYZCTR=',XCTR, YCTR, ZCTR,' WLEN=',WLEN
c
c     identify cube corner direction vectors
c
      xd(1) = - 1.0d0
      yd(1) = - 1.0d0
      zd(1) = + 1.0d0
c
      xd(2) = - 1.0d0
      yd(2) = + 1.0d0
      zd(2) = + 1.0d0
c
      xd(3) = + 1.0d0
      yd(3) = + 1.0d0
      zd(3) = + 1.0d0
c
      xd(4) = + 1.0d0
      yd(4) = - 1.0d0
      zd(4) = + 1.0d0
c
      xd(5) = - 1.0d0
      yd(5) = - 1.0d0
      zd(5) = - 1.0d0
c
      xd(6) = - 1.0d0
      yd(6) = + 1.0d0
      zd(6) = - 1.0d0
c
      xd(7) = + 1.0d0
      yd(7) = + 1.0d0
      zd(7) = - 1.0d0
c
      xd(8) = + 1.0d0
      yd(8) = - 1.0d0
      zd(8) = - 1.0d0
c
c     compute corners of cube
c
      do 50 i=1,8
         x(i)=xctr+wlen*xd(i)
         y(i)=yctr+wlen*yd(i)
         z(i)=zctr+wlen*zd(i)
         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 240
   50 continue
c
c     make coordinates of corners of cube into whole numbers
c
      dfull = dble(mfull)
      if(dabs(x(3)+1.0d0).ge.dfull) stop 242
      if(dabs(y(3)+1.0d0).ge.dfull) stop 243
      if(dabs(z(3)+1.0d0).ge.dfull) stop 244
      if(dabs(x(5)-1.0d0).ge.dfull) stop 245
      if(dabs(y(5)-1.0d0).ge.dfull) stop 246
      if(dabs(z(5)-1.0d0).ge.dfull) stop 247
c
      xmx = dble(idnint(x(3)+1.0d0))
      ymx = dble(idnint(y(3)+1.0d0))
      zmx = dble(idnint(z(3)+1.0d0))
      xmn = dble(idnint(x(5)-1.0d0))
      ymn = dble(idnint(y(5)-1.0d0))
      zmn = dble(idnint(z(5)-1.0d0))
c
      xlan = xmx - xmn
      ylan = ymx - ymn
      zlan = zmx - zmn
      wlan = dmax1(xlan,ylan,zlan) 
c
      x(1) = xmn
      y(1) = ymn
      z(1) = zmn + wlan
c
      x(2) = xmn
      y(2) = ymn + wlan
      z(2) = zmn + wlan
c
      x(3) = xmn + wlan
      y(3) = ymn + wlan
      z(3) = zmn + wlan
c
      x(4) = xmn + wlan
      y(4) = ymn
      z(4) = zmn + wlan
c
      x(5) = xmn
      y(5) = ymn
      z(5) = zmn
c
      x(6) = xmn
      y(6) = ymn + wlan
      z(6) = zmn
c
      x(7) = xmn + wlan
      y(7) = ymn + wlan
      z(7) = zmn
c
      x(8) = xmn + wlan
      y(8) = ymn
      z(8) = zmn
c
      do 55 i=1,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 250
   55 continue
c
c     get cube corner directions in their integer form
c
      do 58 i = 1,8
         xc(i) = idnint(xd(i))
         yc(i) = idnint(yd(i))
         zc(i) = idnint(zd(i))
   58 continue
c
c     compute corners of rectangular polyhedron
c
   85 continue
      if(.not.reccor) go to 165
      ng = 0
      if(artfcl) ng = 8
      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 100 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
  100 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 165
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
  165 continue
      return
      end
*TRANSF
c
c     Subroutine for transforming double precision coordinates into
c     their two-integer decompositions, finding initial vertex, etc.
c
      subroutine transf(x, y, z, ix, iy, iz, ix2, iy2, iz2, icon, is,
     *                 ifl, nv, nw, nt, nmax, nvmax, icfig, nara, nw0,
     *                 nw1, ibeg, isbeg, epz, mhalf, mfull, isclp,
     *                 iftal, deps, dscle, dfull, dfill, artfcl)
c
      integer nmax, nvmax
      double precision x(*), y(*), z(*)
      integer ix(*), iy(*), iz(*)
      integer ix2(*), iy2(*), iz2(*)
      integer icon(8,*), is(*), ifl(*)
      integer nv, nw, nt, icfig, nara, ibeg, isbeg
      logical artfcl
      double precision xmin, xmax, ymin, ymax, zmin, zmax, wbig, epz
      double precision deps, dscle, dfull, dfill
      integer isclp(*), mhalf, mfull, ibfig, itfig, nw0, nw1
      integer idmin, i, j, no, iftal
c
c     initialize for endpoints
c
      nw = nv
      nw0 = nw+1
      nw1 = nw+2
      if(nw1 .gt. nmax) stop 300
c
c     initialize for number of artificial points
c
      nara = 0
      if(artfcl) nara = 8
c
c     initialize Fortran 77 word lengths
c
      mhalf=32768
      mfull=1073741824
      deps = dble(0.9)
c
c     testing parameters
c
      if(nv .lt.1 .or. nv .gt. nmax) stop 310
      if(nt .lt.1 .or. nt .gt.nvmax) stop 320
c
c     initialize arrays
c
      do 100 i=1,nmax
         is(i) = 0
  100 continue
c
c     test array icon and set array is
c
      do 230 i=1,nt
         do 210 j=1,4
            if(icon(j,i).lt.0 .or. icon(j,i).gt.nt) stop 330
  210    continue
         do 220 j=5,8
            if(icon(j,i).le.0 .or. icon(j,i).gt.nv) stop 340
            is(icon(j,i))=i
  220    continue
  230 continue
c
c     find a vertex in tetrahedralization
c
      do 240 i = nara+1, nw
         if(is(i).gt.0) go to 250
  240 continue
      stop 350
  250 continue
      ibeg = i
      isbeg = is(ibeg)
      if(icon(5,isbeg).ne.ibeg .and. icon(6,isbeg).ne.ibeg .and.
     *   icon(7,isbeg).ne.ibeg .and. icon(8,isbeg).ne.ibeg) stop 360
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)
      do 260 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)
  260 continue
c
c     test # of significant figures of nondecimal part of coordinates
c
      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(*,*)'WBIG=',WBIG
      ibfig = 0
  280 continue
      ibfig = ibfig+1
      wbig = wbig/10.0d0
      if(wbig .ge. 1.0d0) go to 280
      if(ibfig.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 405
      endif
      itfig = ibfig + icfig
c     WRITE(*,*)'ITFIG=',ITFIG,' IBFIG=',IBFIG,' ICFIG=',ICFIG
      if(itfig.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 ',itfig
         write(*,*)'Program is terminated since the maximum ',
     *             'allowed is 14.'
         stop 410
      endif
c
c     transform input double precision coordinates into their integer
c     decomposition according to the specified number of significant
c     figures on decimal part of coordinates
c
      call intran(x, y, z, ix, iy, iz, ix2, iy2, iz2, nv, mhalf,
     *            mfull, icfig, isclp, deps, dscle, dfull, dfill)
c
c     redefine number of significant figures of decimal part of
c     coordinates to be used during the current execution of program
c
c
c     test tetrahedralization
c
      call consis(icon, is, ifl, nw, nt)
c
      call orient(nt, icon, x, y, z, ix, iy, iz, ix2, iy2, iz2, idmin,
     *            mhalf, mfull, isclp, epz, artfcl)
      if(idmin.ne.0) then
         write(*,*)' '
         write(*,*)'In input tetrahedralization, ',
     *             'orientation violations detected.'
         write(*,*)'Number of violations = ',idmin
         write(*,*)'Program is terminated.'
         stop 415
      endif
c
c     initialize irray ifl
c
      iftal = 0
      do 150 i = 1, nvmax
         ifl(i) = 0
  150 continue
c
      return
      end
*INTRAN
c
c     subroutine intran to -
c
c     transform double precision coordinates into their integer
c     decomposition
c
      subroutine intran(x, y, z, ix, iy, iz, ix2, iy2, iz2, nv,
     *                  mhalf, mfull, icfig, isclp, deps, dscle,
     *                  dfull, dfill)
c
      double precision x(*), y(*), z(*)
      integer ix(*), iy(*), iz(*)
      integer ix2(*), iy2(*), iz2(*)
      integer nv, mhalf, mfull, icfig, isclp(*)
c
      double precision deps, dscle, dfull, dfill, decml
      integer isgcl, isclu, i
c
c     test # of significant figures of decimal part of coordinates
c
      if(icfig.lt.0 .or. icfig.gt.9) stop 510
      isclu = 1
      dscle = 1.0d0
      if(icfig.eq.0) go to 220
      do 210 i = 1, icfig
         isclu = 10*isclu
         dscle = 10.0d0*dscle
  210 continue
  220 continue
      if(iabs(isclu).ge.mfull) stop 520
      call decomp(isclp, isgcl, isclu, mhalf)
      if(isgcl.ne.1) stop 530
c
c     test lengths of x, y, z-coordinates, shift and make them integers
c
      dfull = dble(mfull)
      if(dscle.lt.deps) stop 540
      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 550
         ix(i) = idint(x(i))
         if(iabs(ix(i)).ge.mfull) stop 560
         decml = (x(i) - dint(x(i)))*dscle
         if(dabs(decml).ge.dfull) stop 570
         ix2(i) = idnint(decml)
         if(iabs(ix2(i)).ge.mfull) stop 580
         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 590
         iy(i) = idint(y(i))
         if(iabs(iy(i)).ge.mfull) stop 600
         decml = (y(i) - dint(y(i)))*dscle
         if(dabs(decml).ge.dfull) stop 610
         iy2(i) = idnint(decml)
         if(iabs(iy2(i)).ge.mfull) stop 620
         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 630
         iz(i) = idint(z(i))
         if(iabs(iz(i)).ge.mfull) stop 640
         decml = (z(i) - dint(z(i)))*dscle
         if(dabs(decml).ge.dfull) stop 650
         iz2(i) = idnint(decml)
         if(iabs(iz2(i)).ge.mfull) stop 660
         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
      return
      end
*PNTLOC
c 
c     Driver subroutine for locating a point in a tetrahedralization by
c     moving in a straight line from a point whose location is known.
c
c     Subroutine calls subroutine shishk for the purpose of finding
c     the location of target point.
c     The target point is represented by iend and it is found by
c     moving in a straight line from a point represented by ibeg;
c     itype, iscur, site1, site2 are input/output arguments for
c     subroutine shiskb.
c
c     On input/output the following conventions are used when
c     calling shishk:
c
c     itype = 1 then ibeg/iend is vertex site1 of iscur
c     itype = 2 then ibeg/iend is in the interior of tetrahedron iscur
c     itype = 3 then ibeg/iend is in edge of iscur with vertices site1
c               and site2
c     itype = 4 then ibeg/iend is in facet of iscur opposite vertex
c               site1 of iscur
c
c     On output only if itype equals -1 then iend is outside union of
c     tetrahedra.
c
      subroutine pntloc(xi, yi, zi, x, y, z, x2, y2, z2, icon, ifl, is,
     *                  nt, nw0, ityp0, iscu0, sit01, sit02, nw1,
     *                  itype, iscur, site1, site2, xc, yc, zc, epz,
     *                  mhalf, mfull, isclp, iftal, dscle, dfull,
     *                  dfill, xcor, ycor, zcor, nara)
c
c     declaration of variables
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*)
      integer icon(8,*), ifl(*), is(*)
      integer xc(*), yc(*), zc(*)
      integer isclp(*), mhalf, mfull, nara, nt
      integer site1, site2, sit01, sit02
      double precision dscle, dfull, dfill
      double precision xcor, ycor, zcor
      double precision epz, decml
      integer iftal, nw0, nw1
      integer ibeg, iend, iscur, ityp0, iscu0
      integer itype
c
c     initialize for destination point
c
      ibeg = nw0
      iend = nw1
      itype = ityp0
      iscur = iscu0
      site1 = sit01
      site2 = sit02
      xi(nw1) = xcor
      yi(nw1) = ycor
      zi(nw1) = zcor
      x2(nw1) = 0
      y2(nw1) = 0
      z2(nw1) = 0
      if(dabs(xi(nw1)).lt.dfill) then
         x(nw1) = idnint(dscle*xi(nw1))
         if(iabs(x(nw1)).lt.mfull) then
            xi(nw1) = dble(x(nw1))/dscle
            go to 322
         endif
      endif
      if(dabs(xi(nw1)).ge.dfull) stop 700
      x(nw1) = idint(xi(nw1))
      if(iabs(x(nw1)).ge.mfull) stop 710
      decml = (xi(nw1) - dint(xi(nw1)))*dscle
      if(dabs(decml).ge.dfull) stop 720
      x2(nw1) = idnint(decml)
      if(iabs(x2(nw1)).ge.mfull) stop 730
      if(iabs(x2(nw1)).eq.0) then
         xi(nw1) = dble(x(nw1))
         x2(nw1) = mfull
      else
         xi(nw1) = dble(x(nw1)) + (dble(x2(nw1))/dscle)
      endif
  322 continue
      if(dabs(yi(nw1)).lt.dfill) then
         y(nw1) = idnint(dscle*yi(nw1))
         if(iabs(y(nw1)).lt.mfull) then
            yi(nw1) = dble(y(nw1))/dscle
            go to 324
         endif
      endif
      if(dabs(yi(nw1)).ge.dfull) stop 740
      y(nw1) = idint(yi(nw1))
      if(iabs(y(nw1)).ge.mfull) stop 750
      decml = (yi(nw1) - dint(yi(nw1)))*dscle
      if(dabs(decml).ge.dfull) stop 760
      y2(nw1) = idnint(decml)
      if(iabs(y2(nw1)).ge.mfull) stop 770
      if(iabs(y2(nw1)).eq.0) then
         yi(nw1) = dble(y(nw1))
         y2(nw1) = mfull
      else
         yi(nw1) = dble(y(nw1)) + (dble(y2(nw1))/dscle)
      endif
  324 continue
      if(dabs(zi(nw1)).lt.dfill) then
         z(nw1) = idnint(dscle*zi(nw1))
         if(iabs(z(nw1)).lt.mfull) then
            zi(nw1) = dble(z(nw1))/dscle
            go to 326
         endif
      endif
      if(dabs(zi(nw1)).ge.dfull) stop 780
      z(nw1) = idint(zi(nw1))
      if(iabs(z(nw1)).ge.mfull) stop 790
      decml = (zi(nw1) - dint(zi(nw1)))*dscle
      if(dabs(decml).ge.dfull) stop 800
      z2(nw1) = idnint(decml)
      if(iabs(z2(nw1)).ge.mfull) stop 810
      if(iabs(z2(nw1)).eq.0) then
         zi(nw1) = dble(z(nw1))
         z2(nw1) = mfull
      else
         zi(nw1) = dble(z(nw1)) + (dble(z2(nw1))/dscle)
      endif
  326 continue
c
c     locate point iend (=nw1)
c
      call shishk(xi, yi, zi, x, y, z, x2, y2, z2, is, icon, ifl,
     *            ibeg, iend, itype, iscur, site1, site2, iftal, nt,
     *            nara, xc, yc, zc, mhalf, mfull, epz, isclp)
c
      return
      end
*SHISHK
c
c     shishkebab routines
c
c     subroutine shishk to -
c
c     move from a point whose location is known in the
c     tetrahedralization to a tetrahedron that contains
c     a second point, and identify the type of location of
c     the second point with respect to the tetrahedron
c
c     On input/output the following conventions are used when
c     calling shishk:
c
c     itype = 1 then ileft/k is vertex site1 of iscur
c     itype = 2 then ileft/k is in the interior of tetrahedron iscur
c     itype = 3 then ileft/k is in edge of iscur with vertices site1
c               and site2
c     itype = 4 then ileft/k is in facet of iscur opposite vertex
c               site1 of iscur
c
c     On output only if itype equals -1 then k is outside union of
c     tetrahedra.
c
      subroutine shishk(xi, yi, zi, x, y, z, x2, y2, z2, is, icon, id,
     *                  ileft, k, itype, iscur, site1, site2, iftal,
     *                  ivnxt, nara, xc, yc, zc, mhalf, mfull,
     *                  epz, isclp)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*)
      integer x2(*), y2(*), z2(*)
      integer is(*), icon(8,*), id(*)
      integer xc(*), yc(*), zc(*)
      integer ileft, k, itype, iscur, iftal, ivnxt, nara
      integer isclp(*), mhalf, mfull
      integer isite(4), a , b, c, d, side1, side2, site0, site1, site2
      integer i, imist, isyde, islst, isadj, ilift, isini
      double precision epz
c
      INTEGER ITUPE, ISID1, ISID2
c
c     reintialize array id if necessary
c
      if(iftal.gt.10000000) then
         iftal = 0
         do 50 i = 1, ivnxt
            id(i) = 0
   50    continue
      endif
c
c     test two points
c
c     WRITE(*,*)'ILEFT=',ILEFT,' ITYPE=',ITYPE,' IDENI=',IDENI
      if(ileft.le.nara .or. k.le.nara) stop 815
c
c
c     start shishk from ileft according to its position
c
      if(itype.eq.1) then
         if(site1.le.nara) stop 820
         a = site1
         go to 300
      elseif(itype.eq.2) then
c        WRITE(*,*)'ISCUR=',ISCUR
         itype = 0
         do 210 i = 5,8
            a = icon(i,iscur)
            if(a.gt.nara) go to 215
  210    continue
         stop 825
  215    continue
         call sitord(icon, a, iscur)
         b = icon(6,iscur)
         c = icon(7,iscur)
         d = icon(8,iscur)
         call trtest(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, itype,
     *               k, a, b, c, d, nara, side1, side2, mhalf, mfull,
     *               isclp, epz)
         if(itype.gt.0) go to 9000
         do 220 i = 1,4
            isite(i) = icon(i+4,iscur)
  220    continue
         do 230 i = 1,4
            a = isite(i)
            call sitord(icon, a, iscur)
            b = icon(6,iscur)
            c = icon(7,iscur)
            d = icon(8,iscur)
c           WRITE(*,*)'FCTEST 1'
            call fctest(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc,
     *                  itype, k, imist, ileft, b, c, d, nara, side1,
     *                  side2, mhalf, mfull, isclp, epz)
            if(itype .gt. 0) stop 830
            if(itype .eq. 0) go to 230
            if(itype .eq. -2) go to 1100
            if(itype .eq. -3) then
               a = imist
               go to 300
            elseif(itype .eq. -4) then
               site0 = a
               site1 = imist
               go to 2000
            else
               stop 840
            endif
  230    continue
c        WRITE(13,*)'ISCUR=',ISCUR,' IFL=',IFL(ISCUR),' ABCD=',A,B,C,D
c        WRITE(13,*)'ILEFT=',ILEFT,' K=',K,' IDENI=',IDENI
         stop 850
      elseif(itype.eq.3) then
         call reordr(icon, site1, site2, iscur)
         site0 = icon(7,iscur)
         go to 2000
      elseif(itype.eq.4) then
         itype = 0
         call sitord(icon, site1, iscur)
         b = icon(6,iscur)
         c = icon(7,iscur)
         d = icon(8,iscur)
         call fcsign(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, isyde,
     *               k, b, c, d, nara, mhalf, mfull, isclp, epz)
         if(isyde.lt.0) then
            islst = iscur
            isadj = iabs(icon(1,iscur))
            if(isadj.le.0) then
               itype = -1
               go to 9000
            endif
            if(isadj.gt.ivnxt) stop 860
            iscur = isadj
            if(iabs(icon(1,iscur)) .eq. islst) then
               ilift = icon(5,iscur)
            elseif(iabs(icon(2,iscur)) .eq. islst) then
               ilift = icon(6,iscur)
            elseif(iabs(icon(3,iscur)) .eq. islst) then
               ilift = icon(7,iscur)
            elseif(iabs(icon(4,iscur)) .eq. islst) then
               ilift = icon(8,iscur)
            else
               stop 870
            endif
            isyde = -isyde
         else
            ilift = site1
            c = icon(8,iscur)
            d = icon(7,iscur)
         endif
c        WRITE(*,*)'FTFIND'
         call ftfind(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, itype,
     *               ileft, k, ilift, imist, b, c, d, nara, side1,
     *               side2, mhalf, mfull, isclp, epz, isyde)
         if(itype .gt. 0) then
            call reordr(icon, ilift, b, iscur)
            go to 9000
         elseif(itype .eq. -2) then
            a = imist
            call sitord(icon, a, iscur)
            b = icon(6,iscur)
            c = icon(7,iscur)
            d = icon(8,iscur)
            go to 1100
         elseif(itype .eq. -3) then
            a = imist
            go to 300
         elseif(itype .eq. -4) then
            if(isyde.gt.0) then
               if(imist.eq.b)then
                  site0 = c
               elseif(imist.eq.c)then
                  site0 = d
               else
                  site0 = b
               endif
            else
               site0 = ilift
            endif
            site1 = imist
            go to 2000
         else
            stop 880
         endif
      endif
      stop 890
c
  300 continue
c
c     find tetrahedron with point a as a vertex for which the ray with
c     origin point a and through point k intersects the facet of the
c     tetrahedron opposite to point a
c
      itype = 0
      iftal = iftal + 1
      iscur = is(a)
      if(iscur.le.0.or.iscur.gt.ivnxt) stop 910
      isini = iscur
c
c     reorder isini so that vertex a equals icon(5,isini)
c
      call sitord(icon, a, isini)
c
c     test current facet
c
  400 continue
      b = icon(6,iscur)
      c = icon(7,iscur)
      d = icon(8,iscur)
      id(iscur) = iftal
c     WRITE(*,*)'FCTEST 2'
      call fctest(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, itype, k,
     *            imist, a, b, c, d, nara, side1, side2, mhalf, mfull,
     *            isclp, epz)
      if(itype .eq. 0) go to 500
      if(itype .gt. 0) go to 9000
      if(itype .eq. -2) go to 1100
      if(itype .eq. -3) then
         a = imist
         go to 300
      elseif(itype .eq. -4) then
         site0 = a
         site1 = imist
         go to 2000
      else
         stop 920
      endif
c
c     obtain next tetrahedron with point a as a vertex
c
  500 continue
      isadj = iabs(icon(2,iscur))
      if(isadj.le.0) go to 600
      if(isadj.gt.ivnxt) stop 930
      if(id(isadj) .eq. iftal) go to 600
      ilift = icon(8,iscur)
      go to 900
  600 continue
      isadj = iabs(icon(3,iscur))
      if(isadj.le.0) go to 700
      if(isadj.gt.ivnxt) stop 940
      if(id(isadj) .eq. iftal) go to 700
      ilift = icon(6,iscur)
      go to 900
  700 continue
      isadj = iabs(icon(4,iscur))
      if((iscur.ne.isini.and.isadj.le.0).or.isadj.gt.ivnxt) stop 950
      if(iscur .eq. isini) go to 800
      if(iabs(icon(3,isadj)) .eq. iscur) then
          iscur = isadj
          go to 700
      elseif(iabs(icon(2,isadj)) .eq. iscur) then
          iscur = isadj
          go to 600
      elseif(iabs(icon(4,isadj)) .eq. iscur) then
          if(isadj .ne. isini) stop 960
          go to 1000
      else
          stop 970
      endif
  800 continue
      if(isadj.le.0) go to 1000
      if(id(isadj) .eq. iftal) go to 1000
      ilift = icon(7,iscur)
c
c     reorder isadj so that a equals icon(5,isadj) and ilift
c     equals icon(6,isadj)
c
  900 continue
      call reordr(icon, a, ilift, isadj)
      iscur = isadj
      go to 400
c
c     can not find intersected tetrahedron
c
 1000 continue
      itype = -1
      go to 9000
c
c     obtain next tetrahedron along line segment as it crosses a facet
c
 1100 continue
      isadj = iabs(icon(1,iscur))
      if(isadj.le.0) then
         itype = -1
         go to 9000
      endif
      if(isadj.gt.ivnxt) stop 1005
      islst = iscur
      iscur = isadj
      if(iabs(icon(1,iscur)) .eq. islst) then
          ilift = icon(5,iscur)
      elseif(iabs(icon(2,iscur)) .eq. islst) then
          ilift = icon(6,iscur)
      elseif(iabs(icon(3,iscur)) .eq. islst) then
          ilift = icon(7,iscur)
      elseif(iabs(icon(4,iscur)) .eq. islst) then
          ilift = icon(8,iscur)
      else
          stop 1010
      endif
c
c     obtain opposite facet of tetrahedron intersected by line
c     segment
c
c     WRITE(*,*)'FCFIND'
      call fcfind(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, itype,
     *            ileft, k, ilift, imist, b, c, d, nara, side1, side2,
     *            mhalf, mfull, isclp, epz)
      if(itype .gt. 0) then
         call reordr(icon, ilift, b, iscur)
         go to 9000
      elseif(itype .eq. -2) then
         call sitord(icon, imist, iscur)
         b = icon(6,iscur)
         c = icon(7,iscur)
         d = icon(8,iscur)
         go to 1100
      elseif(itype .eq. -3) then
         a = ilift
         go to 300
      elseif(itype .eq. -4) then
         if(imist.eq.b)then
            site0 = c
         elseif(imist.eq.c)then
            site0 = d
         else
            site0 = b
         endif
         site1 = imist
         go to 2000
      else
         stop 1074
      endif
c
c     obtain next tetrahedron along line segment as it crosses an edge
c
 2000 continue
c     WRITE(*,*)'FCEDGE'
      call fcedge(x, y, z, x2, y2, z2, xc, yc, zc, itype, ileft, k,
     *            icon, iscur, imist, ivnxt, site0, site1, side1, side2,
     *            nara, mhalf, mfull, isclp)
      if(itype .gt. 0 .or. itype.eq.-1) go to 9000
      if(itype .eq. -2) then
         call sitord(icon, imist, iscur)
         b = icon(6,iscur)
         c = icon(7,iscur)
         d = icon(8,iscur)
         go to 1100
      elseif(itype .eq. -3) then
         a = imist
         go to 300
      elseif(itype.eq.-4) then
         if(imist.gt.0) then
            if(imist.eq.site1)then
               site0 = icon(7,iscur)
            elseif(imist.eq.site0)then
               site0 = site1
            endif
         else
            imist = -imist
            if(imist.eq.site1)then
               site0 = icon(8,iscur)
            else
               site0 = site1
            endif
         endif
         site1 = imist
         go to 2000
      else
         stop 1090
      endif
c
 9000 continue
c
c     test position of k
c
      if(itype.eq.-1) go to 9060
      itupe = 0
      do 9020 i = 5,8
         a = icon(i,iscur)
         if(a.gt.nara) go to 9040
 9020 continue
      stop 1095
 9040 continue
      if(a.eq.icon(5,iscur)) then
         b = icon(6,iscur)
         c = icon(7,iscur)
         d = icon(8,iscur)
      elseif(a.eq.icon(6,iscur)) then
         b = icon(5,iscur)
         c = icon(8,iscur)
         d = icon(7,iscur)
      elseif(a.eq.icon(7,iscur)) then
         b = icon(5,iscur)
         c = icon(6,iscur)
         d = icon(8,iscur)
      else
         b = icon(5,iscur)
         c = icon(7,iscur)
         d = icon(6,iscur)
      endif
      call trtest(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, itupe, k,
     *            a, b, c, d, nara, isid1, isid2, mhalf, mfull, isclp,
     *            epz)
      if(itupe.ne.itype) then
         WRITE(*,*)'ILEFT=',ILEFT,' K=',K,' ISCUR=',ISCUR
         WRITE(*,*)'ITUPE=',ITUPE,' ITYPE=',ITYPE
         WRITE(*,*)'ISID1=',ISID1,' ISID2=',ISID2
         WRITE(*,*)'ABCD=',A,B,C,D
         stop 1097
      endif
 9060 continue
c
      if(itype .eq. 1) then
         site1 = icon(side1+4,iscur)
         site2 = 0
      elseif(itype .eq. 2) then
         site1 = 0
         site2 = 0
      elseif(itype .eq. 3) then
         site1 = icon(side1+4,iscur)
         site2 = icon(side2+4,iscur)
         call reordr(icon, site1, site2, iscur)
         site1 = icon(7,iscur)
         site2 = icon(8,iscur)
      elseif(itype .eq. 4) then
         site1 = icon(side1+4,iscur)
         site2 = 0
      elseif(itype .eq. -1) then
         iscur = 0
         site1 = 0
         site2 = 0
      else
         stop 1120
      endif
c
      return
      end
*TRTEST
c
c     This subroutine will test whether a point k is in tetrahedron
c     with vertices a, b, c, d. Vertex a is assumed not to be
c     artificial.
c
      subroutine trtest(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc,
     *                  itype, k, a, b, c, d, nara, side1, side2, mhalf,
     *                  mfull, isclp, epz)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer itype, k, nara, mhalf, mfull, isclp(*)
      integer iside(4), a, b, c, d, side1, side2
      double precision epz
      integer iasign, ipout
c
c     determine whether ray with origin point a and through point k 
c     intersects facet of tetrahedron opposite to point a
c
      if(b.le.nara .or. c.le.nara .or. d.le.nara) go to 100
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, a, c, d,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(2) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, a, d, b,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(3) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, a, b, c,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(4) = ipout
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 1000
c
c     determine whether point k is in tetrahedron
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, b, d, c,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(1) = ipout
      if(iside(1).lt.0) go to 1000
c
   50 continue
      call pntype(iside, itype, side1, side2)
      go to 1000
c
c     there is at least one artificial vertex
c
  100 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, a, c, d, k,
     *            mhalf, mfull, isclp, iasign)
      iside(2) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, a, d, b, k,
     *            mhalf, mfull, isclp, iasign)
      iside(3) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, a, b, c, k,
     *            mhalf, mfull, isclp, iasign)
      iside(4) = iasign
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 1000
c
c     determine whether point k is in tetrahedron
c
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, b, d, c, k,
     *            mhalf, mfull, isclp, iasign)
      iside(1) = iasign
      if(iside(1).lt.0) go to 1000
      go to 50
c
 1000 continue
      return
      end
*FCSIGN
c
c     This subroutine will determine where with respect to facet
c     with vertices b, c, d, point k is.
c
      subroutine fcsign(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc,
     *                  isyde, k, b, c, d, nara, mhalf, mfull, isclp,
     *                  epz)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer b, c, d, isyde, k, nara, mhalf, mfull, isclp(*)
      double precision epz
      integer iasign
c
      if(b.le.nara .or. c.le.nara .or. d.le.nara) go to 100
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, b, d, c,
     *            mhalf, mfull, isclp, epz, isyde)
      go to 1000
c
c     there is at least one artificial vertex
c
  100 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, b, d, c, k,
     *            mhalf, mfull, isclp, iasign)
      isyde = iasign
c
 1000 continue
      return
      end
*FTFIND
c
c     This subroutine tests whether a point on a ray with origin in
c     the interior of a facet of a tetrahedron is in the tetrahedron
c     and if not finds other facet of the tetrahedron intersected by
c     the ray
c
      subroutine ftfind(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc,
     *                  itype, ileft, k, ilift, imist, b, c, d, nara,
     *                  side1, side2, mhalf, mfull, isclp, epz, isyde)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer iside(4), b, c, d, side1, side2, isclp(*)
      integer itype, ileft, k, ilift, imist, nara, mhalf, mfull
      double precision epz
      integer isyde, iasign, idot1, idot2, idot3, ipout
c
c     determine whether point k is in tetrahedron
c
      if(b.le.nara .or. c.le.nara .or. d.le.nara .or. ilift.le.nara)
     *   go to 100
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ilift, d, c,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(2) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ilift, c, b,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(3) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ilift, b, d,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(4) = ipout
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 200
c
   50 continue
      iside(1) = isyde
      call pntype(iside, itype, side1, side2)
      go to 1000
c
c     there is at least one artificial vertex
c
  100 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, d, c, k,
     *            mhalf, mfull, isclp, iasign)
      iside(2) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, b, d, k,
     *            mhalf, mfull, isclp, iasign)
      iside(4) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, c, b, k,
     *            mhalf, mfull, isclp, iasign)
      iside(3) = iasign
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 500
      go to 50
c
c     k is not in tetrahedron
c
  200 continue
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, b,
     *            ilift, mhalf, mfull, isclp, epz, idot1)
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, c,
     *            ilift, mhalf, mfull, isclp, epz, idot2)
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, d,
     *            ilift, mhalf, mfull, isclp, epz, idot3)
      go to 600
c
  500 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, b, ilift, k,
     *            mhalf, mfull, isclp, idot1)
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, c, ilift, k,
     *            mhalf, mfull, isclp, idot2)
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, d, ilift, k,
     *            mhalf, mfull, isclp, idot3)
c
  600 continue
      if(isyde.eq.0) go to 700
      itype = -2
      if(idot1 .lt. 0 .and. idot2 .gt. 0) then
          imist = d
      elseif(idot2 .lt. 0 .and. idot3 .gt. 0) then
          imist = b
      elseif(idot3 .lt. 0 .and. idot1 .gt. 0) then
          imist = c
      elseif(idot1 .eq. 0 .and. idot2 .eq. 0 .and. idot3 .eq.0) then
          itype = -3
          imist = ilift
      elseif(idot1 .eq. 0) then
          itype = -4
          imist = b
      elseif(idot2 .eq. 0) then
          itype = -4
          imist = c
      elseif(idot3 .eq. 0) then
          itype = -4
          imist = d
      else
          stop 1180
      endif
      go to 1000
c
  700 continue
      itype = -4
      if(idot1 .lt. 0 .and. idot2 .gt. 0) then
          imist = c
      elseif(idot2 .lt. 0 .and. idot3 .gt. 0) then
          imist = d
      elseif(idot3 .lt. 0 .and. idot1 .gt. 0) then
          imist = b
      elseif(idot1 .eq. 0 .and. idot2 .eq. 0 .and. idot3 .eq.0) then
          stop 1190
      elseif(idot1 .eq. 0) then
          itype = -3
          imist = b
      elseif(idot2 .eq. 0) then
          itype = -3
          imist = c
      elseif(idot3 .eq. 0) then
          itype = -3
          imist = d
      else
          stop 1200
      endif
c
 1000 continue
      return
      end
*FCTEST
c
c     This subroutine will test whether a ray with origin a vertex of
c     a tetrahedron intersects the facet opposite the vertex of the
c     tetrahedron and whether a point in the interior of the ray is
c     contained in the tetrahedron
c
      subroutine fctest(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc,
     *                  itype, k, imist, a, b, c, d, nara, side1, side2,
     *                  mhalf, mfull, isclp, epz)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      double precision epz
      integer isclp(*), itype, k, imist, mhalf, mfull, iasign, ipout
      integer iside(4), a, b, c, d, side1, side2, nara
c
c     determine whether ray with origin point a and through point k 
c     intersects facet of current tetrahedron opposite to point a
c
      if(b.le.nara .or. c.le.nara .or. d.le.nara) go to 100
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, a, c, d,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(2) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, a, d, b,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(3) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, a, b, c,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(4) = ipout
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 1000
c
c     determine whether point k is in tetrahedron
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, b, d, c,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(1) = ipout
      if(iside(1).lt.0) go to 500
c
   50 continue
      call pntype(iside, itype, side1, side2)
      go to 1000
c
c     there is at least one artificial vertex
c
  100 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, a, c, d, k,
     *            mhalf, mfull, isclp, iasign)
      iside(2) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, a, d, b, k,
     *            mhalf, mfull, isclp, iasign)
      iside(3) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, a, b, c, k,
     *            mhalf, mfull, isclp, iasign)
      iside(4) = iasign
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 1000
c
c     determine whether point k is in tetrahedron
c
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, b, d, c, k,
     *            mhalf, mfull, isclp, iasign)
      iside(1) = iasign
      if(iside(1).lt.0) go to 500
      go to 50
c
c     ray intersects facet but point k is not in tetrahedron
c
  500 continue
c
c     ray intersects interior of facet
c
      if(iside(2).gt.0 .and. iside(3).gt.0 .and. iside(4).gt.0) then
         itype = -2
         go to 1000
      endif
c
      if(iside(2).eq.0 .and. iside(3).eq.0 .and. iside(4).eq.0)
     *   stop 1220
c
c     ray intersects a vertex of facet
c
      if (iside(2).eq.0 .and. iside(3).eq.0) then
         itype = -3
         imist = d
         go to 1000
      elseif (iside(2).eq.0 .and. iside(4).eq.0) then
         itype = -3
         imist = c
         go to 1000
      elseif (iside(3).eq.0 .and. iside(4).eq.0) then
         itype = -3
         imist = b
         go to 1000
      endif
c
c     ray intersects the interior of an edge of facet
c
      itype = -4
      if (iside(2) .eq. 0) then
         imist = c
      elseif (iside(3) .eq. 0) then
         imist = d
      elseif (iside(4) .eq. 0) then 
         imist = b
      else
         stop 1240
      endif
c
 1000 continue
      return
      end
*FCFIND
c
c     This subroutine tests whether a point on a ray that intersects
c     the interior of a facet of a tetrahedron is in the tetrahedron
c     and if not finds other facet of the tetrahedron intersected by
c     the ray
c
      subroutine fcfind(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc,
     *                  itype, ileft, k, ilift, imist, b, c, d, nara,
     *                  side1, side2, mhalf, mfull, isclp, epz)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      double precision epz
      integer isclp(*), mhalf, mfull, nara
      integer itype, ileft, k, ilift, imist, iasign, ipout
      integer idut1, idut2, idut3, idot1, idot2, idot3
      integer iside(4), b, c, d, side1, side2
c
c     determine whether point k is in tetrahedron
c
      if(b.le.nara .or. c.le.nara .or. d.le.nara .or. ilift.le.nara)
     *   go to 100
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ilift, d, c,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(2) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ilift, c, b,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(3) = ipout
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ilift, b, d,
     *            mhalf, mfull, isclp, epz, ipout)
      iside(4) = ipout
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 200
c
   50 continue
      iside(1) = 1
      call pntype(iside, itype, side1, side2)
      go to 1000
c
c     there is at least one artificial vertex
c
  100 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, d, c, k,
     *            mhalf, mfull, isclp, iasign)
      iside(2) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, b, d, k,
     *            mhalf, mfull, isclp, iasign)
      iside(4) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, c, b, k,
     *            mhalf, mfull, isclp, iasign)
      iside(3) = iasign
c
      if(iside(2).lt.0 .or. iside(3).lt.0 .or. iside(4).lt.0) go to 300
      go to 50
c
c     k is not in tetrahedron
c
c     determine position of ilift with repect to current situation
c
  200 continue
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, ileft, ilift, c,
     *            b, mhalf, mfull, isclp, epz, idut1)
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, ileft, ilift, d,
     *            c, mhalf, mfull, isclp, epz, idut2)
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, ileft, ilift, b,
     *            d, mhalf, mfull, isclp, epz, idut3)
c
      if(idut1.le.0 .or. idut2.le.0 .or. idut3.le.0) go to 700
      go to 400
c
c     there is at least one artificial vertex
c
  300 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, c, b, ileft,
     *            mhalf, mfull, isclp, idut1)
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, d, c, ileft,
     *            mhalf, mfull, isclp, idut2)
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ilift, b, d, ileft,
     *            mhalf, mfull, isclp, idut3)
c
      if(idut1.le.0 .or. idut2.le.0 .or. idut3.le.0) go to 700
      go to 500
c
c     ilift, ileft, b, d, c, form a strictly convex set
c
  400 continue
c
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, b,
     *            ilift, mhalf, mfull, isclp, epz, idot1)
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, c,
     *            ilift, mhalf, mfull, isclp, epz, idot2)
      call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, d,
     *            ilift, mhalf, mfull, isclp, epz, idot3)
      go to 600
c
  500 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, b, ilift, k,
     *            mhalf, mfull, isclp, idot1)
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, c, ilift, k,
     *            mhalf, mfull, isclp, idot2)
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, d, ilift, k,
     *            mhalf, mfull, isclp, idot3)
c
  600 continue
      itype = -2
      if(idot1 .lt. 0 .and. idot2 .gt. 0) then
          imist = d
      elseif(idot2 .lt. 0 .and. idot3 .gt. 0) then
          imist = b
      elseif(idot3 .lt. 0 .and. idot1 .gt. 0) then
          imist = c
      elseif(idot1 .eq. 0 .and. idot2 .eq. 0 .and. idot3 .eq.0) then
          itype = -3
      elseif(idot1 .eq. 0) then
          itype = -4
          imist = b
      elseif(idot2 .eq. 0) then
          itype = -4
          imist = c
      elseif(idot3 .eq. 0) then
          itype = -4
          imist = d
      else
          stop 1260
      endif
      go to 1000
c
  700 continue
      if(idut1.le.0 .and. idut2.le.0 .and. idut3.le.0) stop 1280
      itype = -2
      if(idut1.le.0 .and. idut2.le.0)then
         imist = c
      elseif(idut2.le.0 .and. idut3.le.0)then
         imist = d
      elseif(idut3.le.0 .and. idut1.le.0)then
         imist = b
      elseif(idut1.le.0)then
         if(d.gt.nara .and. ilift.gt.nara) then
            call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, d,
     *                  ilift, mhalf, mfull, isclp, epz, idot3)
         else
            call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, d,
     *                  ilift, k, mhalf, mfull, isclp, idot3)
         endif
         if(idot3.gt.0)then
            imist = b
         elseif(idot3.lt.0)then
            imist = c
         else
            itype = -4
            imist = d
         endif
      elseif(idut2.le.0)then
         if(b.gt.nara .and. ilift.gt.nara) then
            call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, b,
     *                  ilift, mhalf, mfull, isclp, epz, idot1)
         else
            call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, b,
     *                  ilift, k, mhalf, mfull, isclp, idot1)
         endif
         if(idot1.gt.0)then
            imist = c
         elseif(idot1.lt.0)then
            imist = d
         else
            itype = -4
            imist = b
         endif
      else
         if(c.gt.nara .and. ilift.gt.nara) then
            call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, ileft, c,
     *                  ilift, mhalf, mfull, isclp, epz, idot2)
         else
            call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, c,
     *                  ilift, k, mhalf, mfull, isclp, idot2)
         endif
         if(idot2.gt.0)then
            imist = d
         elseif(idot2.lt.0)then
            imist = b
         else
            itype = -4
            imist = c
         endif
      endif
c
 1000 continue
      return
      end
*FCEDGE
c
c     This subroutine will test whether a ray through an edge of a
c     tetrahedron intersects either of the facets of the tetrahedron
c     opposite the edge and whether a point in the interior of the
c     ray is contained in the tetrahedron
c
      subroutine fcedge(x, y, z, x2, y2, z2, xc, yc, zc, itype, ileft,
     *                  k, icon, iscur, imist, ivnxt, site0, site1,
     *                  side1, side2, nara, mhalf, mfull, isclp)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer icon(8,*)
      integer iside(4), site0, site1, site2, site3, side1, side2
      integer itype, ileft, k, iscur, imist, ivnxt, nara, mhalf, mfull
      integer isclp(*), isnow, idat, idut, iasign, idot0, ipout
c
c     find intersecting facet
c
      itype = 0
      call reordr(icon, site0, site1, iscur)
      site2 = icon(7,iscur)
      site0 = icon(8,iscur)
      isnow = iabs(icon(1,iscur))
c
  300 continue
      if(isnow.le.0) go to 350
      if(isnow.gt.ivnxt) stop 1380
      call reordr(icon, site0, site1, isnow)
      site3 = icon(8,isnow)
      if(site1.gt.nara .and. site2.gt.nara .and. site0.gt.nara) then
         call ipsign(x, y, z, x2, y2, z2, site1, site0, site2, k,
     *                 mhalf, mfull, isclp, idat)
      else
         call artsig(x, y, z, x2, y2, z2, xc, yc, zc, site1, site0,
     *               site2, k, mhalf, mfull, isclp, idat)
      endif
      if(site1.gt.nara .and. site2.gt.nara .and. site3.gt.nara) then
         call ipsign(x, y, z, x2, y2, z2, site1, site3, site2, k,
     *               mhalf, mfull, isclp, idut)
      else
         call artsig(x, y, z, x2, y2, z2, xc, yc, zc, site1, site3,
     *               site2, k, mhalf, mfull, isclp, idut)
      endif
      if(idat.le.0 .and. idut.ge.0) go to 400
      if(isnow.eq.iscur) stop 1420
      site0 = site3
      isnow = iabs(icon(1,isnow))
      go to 300
c
  350 continue
      site0 = icon(5,iscur)
      site3 = icon(8,iscur)
      isnow = iscur
  380 continue
      if(site1.gt.nara .and. site2.gt.nara .and. site0.gt.nara) then
         call ipsign(x, y, z, x2, y2, z2, site1, site0, site2, k,
     *               mhalf, mfull, isclp, idat)
      else
         call artsig(x, y, z, x2, y2, z2, xc, yc, zc, site1, site0,
     *               site2, k, mhalf, mfull, isclp, idat)
      endif
      if(site1.gt.nara .and. site2.gt.nara .and. site3.gt.nara) then
         call ipsign(x, y, z, x2, y2, z2, site1, site3, site2, k,
     *               mhalf, mfull, isclp, idut)
      else
         call artsig(x, y, z, x2, y2, z2, xc, yc, zc, site1, site3,
     *               site2, k, mhalf, mfull, isclp, idut)
      endif
      if(idat.le.0 .and. idut.ge.0) go to 400
      isnow = iabs(icon(4,isnow))
      if(isnow.le.0) then
         itype = -1
         go to 1000
      endif
      if(isnow.gt.ivnxt .or. isnow.eq.iscur) stop 1450
      site3 = site0
      call reordr(icon, site3, site1, isnow)
      site0 = icon(7,isnow)
      call reordr(icon, site0, site1, isnow)
      go to 380
c
  400 continue
      iscur = isnow
c
c     determine whether point k is in tetrahedron
c
      if(site0.le.nara .or. site1.le.nara .or. site2.le.nara .or.
     *   site3.le.nara) go to 500
c
      call ipsign(x, y, z, x2, y2, z2, site1, site0, site3, k,
     *            mhalf, mfull, isclp, ipout)
      iside(3) = ipout
      call ipsign(x, y, z, x2, y2, z2, site2, site3, site0, k,
     *            mhalf, mfull, isclp, ipout)
      iside(2) = ipout
c
      if(iside(2).lt.0 .or. iside(3).lt.0) go to 600
c
  450 continue
      iside(1) = idut
      iside(4) = -idat
      call pntype(iside, itype, side1, side2)
      go to 1000
c
c     there is at least one artificial vertex
c
  500 continue
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, site1, site0, site3,
     *            k, mhalf, mfull, isclp, iasign)
      iside(3) = iasign
      call artsig(x, y, z, x2, y2, z2, xc, yc, zc, site2, site3, site0,
     *            k, mhalf, mfull, isclp, iasign)
      iside(2) = iasign
c
      if(iside(2).lt.0 .or. iside(3).lt.0) go to 600
      go to 450
c
c     k is not in tetrahedron but ray intersects one of the facets
c     of the tetrahedron opposite the edge
c
  600 continue
      if(site0.gt.nara .and. site3.gt.nara) then
         call ipsign(x, y, z, x2, y2, z2, ileft, site0, site3, k,
     *               mhalf, mfull, isclp, idot0)
      else
         call artsig(x, y, z, x2, y2, z2, xc, yc, zc, ileft, site0,
     *               site3, k, mhalf, mfull, isclp, idot0)
      endif
      if(idat.eq.0) go to 700
      if(idot0.gt.0)then
         if(idut.gt.0) then
            itype = -2
            imist = site1
         else
            itype = -4
            imist = site2
         endif
      elseif(idot0.lt.0)then
         if(idut.gt.0) then
            itype = -2
            imist = site2
         else
            itype = -4
            imist = site1
         endif
      else
         if(idut.gt.0) then
            itype = -4
            imist = site0
         else
            itype = -3
            imist = site3
         endif
      endif
      go to 1000
  700 continue
      if(idot0.gt.0)then
         if(idut.gt.0) then
            itype = -4
            imist = -site2
         else
            itype = -3
            imist = site2
         endif
      elseif(idot0.lt.0)then
         if(idut.gt.0) then
            itype = -4
            imist = -site1
         else
            itype = -3
            imist = site1
         endif
      else
         if(idut.gt.0)then
            itype = -3
            imist = site0
         else
            stop 1480
         endif
      endif
c
 1000 continue
      return
      end
*PNTYPE
c
c     This subroutine determines point type with respect to a
c     tetrahedron that contains a point
c
      subroutine pntype(iside, itype, side1, side2)
c
      integer iside(*), side1, side2, itype
c
c     point is in the interior of tetrahedron
c
      if(iside(1).gt.0 .and. iside(2).gt.0 .and. iside(3).gt.0 .and.
     *   iside(4).gt.0) then
         itype = 2
         go to 1000
      endif
c
c     unacceptable situation
c
      if(iside(1).eq.0 .and. iside(2).eq.0 .and. iside(3).eq.0 .and.
     *   iside(4).eq.0) stop 2010
c
c     point is a vertex of tetrahedron
c
      if(iside(1).eq.0 .and. iside(2).eq.0 .and. iside(3).eq.0) then
         itype = 1
         side1 = 4
         go to 1000
      elseif(iside(1).eq.0 .and. iside(2).eq.0 .and. iside(4).eq.0) then
         itype = 1
         side1 = 3
         go to 1000
      elseif(iside(1).eq.0 .and. iside(3).eq.0 .and. iside(4).eq.0) then
         itype = 1
         side1 = 2
         go to 1000
      elseif(iside(2).eq.0 .and. iside(3).eq.0 .and. iside(4).eq.0) then
         itype = 1
         side1 = 1
         go to 1000
      endif
c
c     point is in the interior of an edge of tetrahedron
c
      if (iside(1).eq.0 .and. iside(2).eq.0) then
         itype = 3 
         side1 = 1
         side2 = 2
         go to 1000
      elseif (iside(1).eq.0 .and. iside(3).eq.0) then
         itype = 3 
         side1 = 1
         side2 = 3
         go to 1000
      elseif (iside(1).eq.0 .and. iside(4).eq.0) then
         itype = 3 
         side1 = 1
         side2 = 4
         go to 1000
      elseif (iside(2).eq.0 .and. iside(3).eq.0) then
         itype = 3 
         side1 = 2
         side2 = 3
         go to 1000
      elseif (iside(2).eq.0 .and. iside(4).eq.0) then
         itype = 3 
         side1 = 2
         side2 = 4
         go to 1000
      elseif (iside(3).eq.0 .and. iside(4).eq.0) then
         itype = 3 
         side1 = 3
         side2 = 4
         go to 1000
      endif
c
c     point is in the interior of a facet of tetrahedron
c
      itype = 4
      if (iside(1) .eq. 0) then
         side1 = 1
      elseif (iside(2) .eq. 0) then
         side1 = 2
      elseif (iside(3) .eq. 0) then
         side1 = 3
      elseif (iside(4) .eq. 0) then 
         side1 = 4
      else
         stop 2020
      endif
c
 1000 continue
      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, i
      integer 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 2815
   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
      return
      end
*ORIENT
c
c     This subroutine will test the orientation of the tetrahedra
c
      subroutine orient(tetra, icon, xi, yi, zi, x, y, z, x2, y2, z2,
     *                  idmin, mhalf, mfull, isclp, epz, artfcl)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*)
      integer tetra, icon(8,*), a, b, c, d, idmin, nara, i
      integer mhalf, mfull, isclp(*), iside
      logical artfcl
      double precision epz
c
c     test all tetrahedra
c
      nara = 0
      if(artfcl) nara = 8
      idmin = 0
      do 200 i=1,tetra
         if(icon(5,i).lt.0) go to 200
         if(icon(5,i) .le. nara .or. icon(6,i) .le. nara .or.
     *     icon(7,i) .le. nara .or. icon(8,i) .le. nara)
     *   go to 200
         a=icon(5,i)
         b=icon(6,i)
         c=icon(7,i)
         d=icon(8,i)
         call irsign(xi, yi, zi, x, y, z, x2, y2, z2, d, a, b, c,
     *               mhalf, mfull, isclp, epz, iside)
         if(iside .le. 0) idmin = idmin+1
  200 continue
c
      return
      end
*ARTSIG
c
c     This subroutine determines the position of a point with respect
c     to a plane when artificial points may be involved
c
      subroutine artsig(x, y, z, x2, y2, z2, xc, yc, zc, ib, id, ic, k,
     *                  mhalf, mfull, isclp, iasign)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer isclp(*), mhalf, mfull
      integer b, d, c, ib, id, ic, k, iasign, ifn
c
      call vrtarr(ib, id, ic, b, d, c)
c
      if(d.gt.8 .and. c.gt.8) then
         call ipsign(x, y, z, x2, y2, z2, b, d, c, k, mhalf,
     *               mfull, isclp, iasign)
         go to 100
      elseif(b.le.8) then
         iasign = 1
         go to 100
      elseif(d.le.8.and.c.le.8) then
         call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, d, c, k, b,
     *               mhalf, mfull, isclp, iasign)
         if(iasign.ne.0) go to 100
         call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, d, c, k, b,
     *               mhalf, mfull, isclp, iasign)
         go to 100
      elseif(d.le.8) then
         ifn = 0
         call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, b, c, k, d,
     *               mhalf, mfull, isclp, ifn, iasign)
      else
         ifn = 1
         call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, b, d, k, c,
     *               mhalf, mfull, isclp, ifn, iasign)
      endif
      if(iasign.ne.0) go to 100
      call ipsign(x, y, z, x2, y2, z2, b, d, c, k, mhalf,
     *            mfull, isclp, iasign)
c
  100 continue
      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
*VRTARR
c
c This routine will arrange vertices b, c, d of a tetrahedron
c so that b>c, b>d. Data structure is not updated
c
      subroutine vrtarr(i2,i3,i4,b,c,d)
c
      integer b, c, d, i2, i3, i4, ix
c
      b=i2
      c=i3
      d=i4
      ix = max0(b,c,d)
      if(b .eq. ix) go to 100
      if(c .eq. ix) then
         ix = b
         b = c
         c = d
         d = ix
      else
         ix = b
         b = d
         d = c
         c = ix
      endif
  100 continue
      if(c.gt.b. or. d.gt.b) stop 3210
c
      return
      end
*IRSIGN
c
c     subroutine to determine position of point site0 with respect
c     to the plane spanned by points site1, site2, site3
c
      subroutine irsign(xi, yi, zi, x, y, z, x2, y2, z2, site0, site1,
     *                  site2, site3, mhalf, mfull, isclp, epz, ipout)
c
      double precision xi(*), yi(*), zi(*)
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*)
      double precision epz, dist
      integer isclp(*), mhalf, mfull, ipossi
      integer site0, site1, site2, site3, ipout
c
      call dstnce(xi, yi, zi, site1, site2, site3, epz, site0, dist,
     *            ipossi)
      if(ipossi.eq.0) then
         ipout = 1
         if(dist.lt.0.0d0) ipout = -1
      else
         call ipsign(x, y, z, x2, y2, z2, site1, site2, site3, site0,
     *               mhalf, mfull, isclp, ipout)
      endif
c
      return
      end
*DSTNCE
c
c This subroutine will compute the distance from a point to a facet of
c a tetrahedron.
c
      subroutine dstnce(x, y, z, p, q, r, epz, k, dist, ipossi)
c
      integer p, q, r, k, ipossi
      double precision x(*), y(*), z(*)
      double precision epz, dist
      double precision xvec1, yvec1, zvec1, xvec2, yvec2, zvec2
      double precision xvec3, yvec3, zvec3, dst1, dst2, dst3
      double precision dotx, doty, dotz, dmax, dlun
      double precision xvecp, yvecp, zvecp, dstp, dlen
c
      ipossi = 0
      xvec1 = x(q) - x(p)
      yvec1 = y(q) - y(p)
      zvec1 = z(q) - z(p)
      xvec2 = x(r) - x(p)
      yvec2 = y(r) - y(p)
      zvec2 = z(r) - z(p)
      xvec3 = x(q) - x(r)
      yvec3 = y(q) - y(r)
      zvec3 = z(q) - z(r)
      dst1=dsqrt(xvec1**2+yvec1**2+zvec1**2)
      dst2=dsqrt(xvec2**2+yvec2**2+zvec2**2)
      dst3=dsqrt(xvec3**2+yvec3**2+zvec3**2)
      if(dst1.lt.epz .or. dst2.lt.epz .or. dst3.lt.epz) then
         ipossi = 1
         go to 1000
      endif
      dmax = dmax1(dst1,dst2,dst3)
c
      dotx = yvec1 * zvec2 - yvec2 * zvec1
      doty = - xvec1 * zvec2 + xvec2 * zvec1
      dotz = xvec1 * yvec2 - xvec2 * yvec1
      dlen = dsqrt (dotx**2 + doty**2 + dotz**2)
      if(dlen.lt.epz .or. dlen/dmax.lt.epz)then
         ipossi = 1
         go to 1000
      endif
c
      xvecp = x(k) - x(p)
      yvecp = y(k) - y(p)
      zvecp = z(k) - z(p)
      dstp=dsqrt(xvecp**2+yvecp**2+zvecp**2)
      if(dstp.lt.epz) then
         ipossi = 1
         go to 1000
      endif
c
      dlun=dstp*dmax
      dlun=dmax1(dlen,dlun)
      dist=(xvecp*dotx+yvecp*doty+zvecp*dotz)/dlun
      if(dist.gt.-epz .and. dist.lt.epz)then
         ipossi = 1
      endif
c
 1000 continue
      return
      end
*IPSIGN
c
c     subroutine for determining position of point ifou with respect
c     to plane that contains points ifir, isec, ithi
c     if positive then ifou is on positive side of plane
c     if negative then ifou is on negative side of plane
c     if zero then ifou is in plane
c
      subroutine ipsign(x, y, z, x2, y2, z2, ifir, isec, ithi,
     *                  ifou, mhalf, mfull, isclp, ipout)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*)
      integer ifir, isec, ithi, ifou
      integer isclp(*), mhalf, mfull, nkmax, ipout
      parameter (nkmax = 30)
      integer io(nkmax), 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
      ipout = isgo
c
      return
      end
*IPSIG3
c
      subroutine ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, ifir, isec,
     *                  ithi, ifou, mhalf, mfull, isclp, ifn, ipout)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer ifir, isec, ithi, ifou
      integer isclp(*), mhalf, mfull, ifn, nkmax, ipout
      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 ix4(nkmax), iy4(nkmax), iz4(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
      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, 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
      ixfow = xc(ifou)
      iyfow = yc(ifou)
      izfow = zc(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
      ikx4 = 2
      iky4 = 2
      ikz4 = 2
      call decomp(ix4, isgx4, ixfow, mhalf)
      call decomp(iy4, isgy4, iyfow, mhalf)
      call decomp(iz4, isgz4, izfow, mhalf)
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(iy4, iz2, io, isgy4, isgz2, isgo, iky4, ikz2, iko,
     *              nkmax, mhalf)
      call mulmul(io, ix3, iv, isgo, isgx3, isgv, iko, ikx3, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz4, ix2, io, isgz4, isgx2, isgo, ikz4, ikx2, iko,
     *              nkmax, mhalf)
      call mulmul(io, iy3, iu, isgo, isgy3, isgu, iko, iky3, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix4, iy2, iv, isgx4, isgy2, isgv, ikx4, iky2, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iz3, iu, isgv, isgz3, isgu, ikv, ikz3, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz4, iy2, io, isgz4, isgy2, isgo, ikz4, iky2, iko,
     *              nkmax, mhalf)
      call mulmul(io, ix3, iu, isgo, isgx3, isgu, iko, ikx3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix4, iz2, iv, isgx4, isgz2, isgv, ikx4, ikz2, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iy3, iu, isgv, isgy3, isgu, ikv, iky3, iku,
     *              nkmax, mhalf)
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iy4, ix2, io, isgy4, isgx2, isgo, iky4, ikx2, iko,
     *              nkmax, mhalf)
      call mulmul(io, iz3, iu, isgo, isgz3, isgu, iko, ikz3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      ipout = isgo
      if(ifn.eq.1) ipout = -isgo
c
      return
      end
*IPSIG4
c
      subroutine ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, ifir, isec,
     *                  ithi, ifou, mhalf, mfull, isclp, ipout)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer ifir, isec, ithi, ifou
      integer isclp(*), mhalf, mfull, nkmax, ipout
      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 ix4(nkmax), iy4(nkmax), iz4(nkmax)
      integer ixf(nkmax), iyf(nkmax), izf(nkmax)
      integer ixthw, iythw, izthw, ixfow, iyfow, izfow
      integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew
      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, iko, iku, ikv
c
      ixthw = x(ithi)
      iythw = y(ithi)
      izthw = z(ithi)
      ixfow = x(ifou)
      iyfow = y(ifou)
      izfow = z(ifou)
c
      ixth2 = x2(ithi)
      iyth2 = y2(ithi)
      izth2 = z2(ithi)
      ixfo2 = x2(ifou)
      iyfo2 = y2(ifou)
      izfo2 = z2(ifou)
c
      ixfiw = xc(ifir)
      iyfiw = yc(ifir)
      izfiw = zc(ifir)
      ixsew = xc(isec)
      iysew = yc(isec)
      izsew = zc(isec)
c
      call decmp2(ixf, isgxf, ikxf, ixfow, ixfo2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfow, iyfo2, mhalf, mfull, isclp)
      call decmp2(izf, isgzf, ikzf, izfow, izfo2, mhalf, mfull, isclp)
c
      ikx2 = 2
      iky2 = 2
      ikz2 = 2
      ikx3 = 2
      iky3 = 2
      ikz3 = 2
      call decomp(ix2, isgx2, ixfiw, mhalf)
      call decomp(iy2, isgy2, iyfiw, mhalf)
      call decomp(iz2, isgz2, izfiw, mhalf)
      call decomp(ix3, isgx3, ixsew, mhalf)
      call decomp(iy3, isgy3, iysew, mhalf)
      call decomp(iz3, isgz3, izsew, mhalf)
c
      call decmp2(io, isgo, iko, ixthw, ixth2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix4, isgo, isgxf, isgx4, iko, ikxf, ikx4,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iythw, iyth2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy4, isgo, isgyf, isgy4, iko, ikyf, iky4,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izthw, izth2, mhalf, mfull, isclp)
      call muldif(io, izf, iz4, isgo, isgzf, isgz4, iko, ikzf, ikz4,
     *              nkmax, mhalf)
c
      call mulmul(iy2, iz3, io, isgy2, isgz3, isgo, iky2, ikz3, iko,
     *              nkmax, mhalf)
      call mulmul(io, ix4, iv, isgo, isgx4, isgv, iko, ikx4, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz2, ix3, io, isgz2, isgx3, isgo, ikz2, ikx3, iko,
     *              nkmax, mhalf)
      call mulmul(io, iy4, iu, isgo, isgy4, isgu, iko, iky4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix2, iy3, iv, isgx2, isgy3, isgv, ikx2, iky3, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iz4, iu, isgv, isgz4, isgu, ikv, ikz4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz2, iy3, io, isgz2, isgy3, isgo, ikz2, iky3, iko,
     *              nkmax, mhalf)
      call mulmul(io, ix4, iu, isgo, isgx4, isgu, iko, ikx4, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix2, iz3, iv, isgx2, isgz3, isgv, ikx2, ikz3, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iy4, iu, isgv, isgy4, isgu, ikv, iky4, iku,
     *              nkmax, mhalf)
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iy2, ix3, io, isgy2, isgx3, isgo, iky2, ikx3, iko,
     *              nkmax, mhalf)
      call mulmul(io, iz4, iu, isgo, isgz4, isgu, iko, ikz4, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      ipout = isgo
c
      return
      end
*IPSIG6
c
      subroutine ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, ifir, isec,
     *                  ithi, ifou, mhalf, mfull, isclp, ipout)
c
      integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*)
      integer ifir, isec, ithi, ifou
      integer isclp(*), mhalf, mfull, nkmax, ipout
      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 ix4(nkmax), iy4(nkmax), iz4(nkmax)
      integer ix5(nkmax), iy5(nkmax), iz5(nkmax)
      integer ix6(nkmax), iy6(nkmax), iz6(nkmax)
      integer ixf(nkmax), iyf(nkmax), izf(nkmax)
      integer ixthw, iythw, izthw, ixfow, iyfow, izfow
      integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew
      integer ixfuw, iyfuw, izfuw, ixsuw, iysuw, izsuw
      integer ixth2, iyth2, izth2, ixfo2, iyfo2, izfo2
      integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2
      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 isgx5, isgy5, isgz5, ikx5, iky5, ikz5
      integer isgx6, isgy6, isgz6, ikx6, iky6, ikz6
      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)
      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
      ixfuw = xc(ifir)
      iyfuw = yc(ifir)
      izfuw = zc(ifir)
      ixsuw = xc(isec)
      iysuw = yc(isec)
      izsuw = zc(isec)
c
      call decmp2(ixf, isgxf, ikxf, ixfow, ixfo2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfow, iyfo2, mhalf, mfull, isclp)
      call decmp2(izf, isgzf, ikzf, izfow, izfo2, mhalf, mfull, isclp)
c
      ikx5 = 2
      iky5 = 2
      ikz5 = 2
      ikx6 = 2
      iky6 = 2
      ikz6 = 2
      call decomp(ix5, isgx5, ixfuw, mhalf)
      call decomp(iy5, isgy5, iyfuw, mhalf)
      call decomp(iz5, isgz5, izfuw, mhalf)
      call decomp(ix6, isgx6, ixsuw, mhalf)
      call decomp(iy6, isgy6, iysuw, mhalf)
      call decomp(iz6, isgz6, izsuw, mhalf)
c
      call decmp2(io, isgo, iko, ixfiw, ixfi2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix2, isgo, isgxf, isgx2, iko, ikxf, ikx2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iyfiw, iyfi2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy2, isgo, isgyf, isgy2, iko, ikyf, iky2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izfiw, izfi2, mhalf, mfull, isclp)
      call muldif(io, izf, iz2, isgo, isgzf, isgz2, iko, ikzf, ikz2,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix3, isgo, isgxf, isgx3, iko, ikxf, ikx3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy3, isgo, isgyf, isgy3, iko, ikyf, iky3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izsew, izse2, mhalf, mfull, isclp)
      call muldif(io, izf, iz3, isgo, isgzf, isgz3, iko, ikzf, ikz3,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, ixthw, ixth2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix4, isgo, isgxf, isgx4, iko, ikxf, ikx4,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iythw, iyth2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy4, isgo, isgyf, isgy4, iko, ikyf, iky4,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, izthw, izth2, mhalf, mfull, isclp)
      call muldif(io, izf, iz4, isgo, isgzf, isgz4, iko, ikzf, ikz4,
     *              nkmax, mhalf)
c
      call mulmul(iy5, iz3, io, isgy5, isgz3, isgo, iky5, ikz3, iko,
     *              nkmax, mhalf)
      call mulmul(io, ix4, iv, isgo, isgx4, isgv, iko, ikx4, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz5, ix3, io, isgz5, isgx3, isgo, ikz5, ikx3, iko,
     *              nkmax, mhalf)
      call mulmul(io, iy4, iu, isgo, isgy4, isgu, iko, iky4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix5, iy3, iv, isgx5, isgy3, isgv, ikx5, iky3, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iz4, iu, isgv, isgz4, isgu, ikv, ikz4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz5, iy3, io, isgz5, isgy3, isgo, ikz5, iky3, iko,
     *              nkmax, mhalf)
      call mulmul(io, ix4, iu, isgo, isgx4, isgu, iko, ikx4, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix5, iz3, iv, isgx5, isgz3, isgv, ikx5, ikz3, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iy4, iu, isgv, isgy4, isgu, ikv, iky4, iku,
     *              nkmax, mhalf)
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iy5, ix3, io, isgy5, isgx3, isgo, iky5, ikx3, iko,
     *              nkmax, mhalf)
      call mulmul(io, iz4, iu, isgo, isgz4, isgu, iko, ikz4, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(iy6, iz2, iv, isgy6, isgz2, isgv, iky6, ikz2, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, ix4, iu, isgv, isgx4, isgu, ikv, ikx4, iku,
     *              nkmax, mhalf)
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz6, ix2, io, isgz6, isgx2, isgo, ikz6, ikx2, iko,
     *              nkmax, mhalf)
      call mulmul(io, iy4, iu, isgo, isgy4, isgu, iko, iky4, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix6, iy2, iv, isgx6, isgy2, isgv, ikx6, iky2, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iz4, iu, isgv, isgz4, isgu, ikv, ikz4, iku,
     *              nkmax, mhalf)
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iz6, iy2, io, isgz6, isgy2, isgo, ikz6, iky2, iko,
     *              nkmax, mhalf)
      call mulmul(io, ix4, iu, isgo, isgx4, isgu, iko, ikx4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix6, iz2, iv, isgx6, isgz2, isgv, ikx6, ikz2, ikv,
     *              nkmax, mhalf)
      call mulmul(iv, iy4, iu, isgv, isgy4, isgu, ikv, iky4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv,
     *              nkmax, mhalf)
c
      call mulmul(iy6, ix2, io, isgy6, isgx2, isgo, iky6, ikx2, iko,
     *              nkmax, mhalf)
      call mulmul(io, iz4, iu, isgo, isgz4, isgu, iko, ikz4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
c
      ipout = isgo
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