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