c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c c Program pwrtet c c This program is for doing 3D nearest point searches. c c Given a set of points S in 3-dimensional space, this program c will identify Power cells in the Power diagram of S that c contain points in a set R in 3-dimensional space, c i.e. for each point in R finds a point in S that is as Power c close to the point in R as any other point in S. c c If no weights are present then Voronoi cells in the c Voronoi diagram of S that contain points in R are identified, c i.e. for each point in R finds a point in S that is as close c to the point in R as any other point in S. c c It is assumed that a Regular/Delaunay tetrahedralization for S c has already been computed using existing program regtet and that c the output tetrahedron list produced by regtet contains real c as well as artificial tetrahedra, i. e. logical variable artfcl c was set to .true. during the execution of regtet. The output c tetrahedron list from regtet is then used as input below. It is c noted that the number of significant figures of decimal part of c point coordinates for points in R is the same as for points in S c which is inherited from the run of regtet. c c Array icl will contain ouput information. If point r is the mth c point in R and n is the absolute value of icl(m) then the Power/ c Voronoi cell of the nth point in S contains point r. c If icl(m) is positive then point r is in the covex hull of S. c In the exterior otherwise. c c The program is based on an algorithm for constructing Regular c tetrahedralizations with incremental topological flipping. c However no actual flippings take place. Given a point in R, c if weights are present a weight is assigned to the point so that c the Power cell of the point in the Power diagram of S together c with the point contains the point. The program then goes through c the motions of inserting the point into the Regular/Delaunay c tetrahedralization for S without actually doing it. This way a c subset of points in S is identified that contains what would be c the Voronoi/Power neighbors of the point in the Voronoi/Power c diagram of S together with the point. The desired point in S c is then a point in this subset closest to the the point in R c in the Power sense. If weights are present care is taken in c choosing the weight assigned to the point in R so that it is as c small as possible. c c Subroutine pwrtet is the driver routine of the program and is c called from the main routine. c c Author: Javier Bernal c National Institute of Standards and Technology (NIST) 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 c Setting of parameters: c c 1. Flipping history used for locating points: c c integer nmax, nvmax, nlmax c parameter (nmax=150000, nvmax=55*nmax, nlmax=150000) c c 2. Flipping history not used for locating points: c integer nmax, nvmax, nlmax parameter (nmax=150000, nvmax= 7*nmax, nlmax=150000) c double precision x(nmax), y(nmax), z(nmax), w(nmax) double precision xa(nlmax), ya(nlmax), za(nlmax) integer ix(nmax), iy(nmax), iz(nmax), iw(nmax) integer ix2(nmax), iy2(nmax), iz2(nmax), iw2(nmax) integer icon(8,nvmax), is(nmax), icl(nlmax), id(nvmax) integer nv, nw, nt, na, nb, naddl, icfig, iwfig double precision wlenx, wleny, wlenz, wlenw, epz logical delaun, pntoff, flphis, artfcl, random, reccor, redchk c double precision xcor, ycor, zcor, wght integer irec, i, j integer ideli, ipnti, iflpi, iarti, irani, ireci, iredi c c---------------------------------------------------------------------- c write(*,*)' ' write(*,*)'This program is for doing 3D nearest point searches.' write(*,*)' ' write(*,*)'It identifies Power cells in Power diagram of a' write(*,*)'set S in 3-d space that contain points in a set R.' write(*,*)' ' write(*,*)'If no weights are present then Voronoi replaces ', * 'Power above.' write(*,*)' ' write(*,*)'It is assumed that a Regular/Delaunay ', * 'tetrahedralization for S' write(*,*)'has already been computed using existing program ', * 'regtet and that' write(*,*)'the output tetrahedron list produced by regtet ', * 'contains real as' write(*,*)'well as artificial tetrahedra (logical variable ', * 'artfcl was set to' write(*,*)'.true. during the execution of regtet).' write(*,*)' ' write(*,*)'The output tetrahedron list from regtet ', * 'is then used as input' write(*,*)'by this program. Note that the number of significant ', * 'figures' write(*,*)'of decimal part of point coordinates for points in R ', * 'is the same' write(*,*)'as for points in S which is inherited from the run ', * 'of regtet.' write(*,*)' ' write(*,*)'Array icl will contain ouput information.' write(*,*)'If point r is the mth point in R and n is the ', * 'absolute value of' write(*,*)'icl(m) then the Power/Voronoi cell of the ', * 'nth point in S ' write(*,*)'contains point r.' write(*,*)'If icl(m) is positive then point r is in the convex ', * 'hull of S.' write(*,*)'In the exterior otherwise.' write(*,*)' ' write(*,*)'In the current main program array icl is saved in ', * 'output file' write(*,*)'targetpoly.' c c---------------------------------------------------------------------- c c open files c open (11, file = 'pnts-wghts') open (12, file = 'tetrahedra') open (13, file = 'targetpnts') open (14, file = 'targetpoly') c c---------------------------------------------------------------------- c c read Regular/Delaunay 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 c test for presence of artificial points and rectangular grid c if(.not.artfcl) then write(*,*)' ' write(*,*)'Input tetrahedra do not contain artificial points.' write(*,*)'Program terminated.' stop 20 endif c c---------------------------------------------------------------------- c c code for avoiding certain warning messages during compilation c if(pntoff) ipnti = 1 if(random) irani = 1 if(redchk) iredi = 1 c c---------------------------------------------------------------------- c c set tolerance c epz = 0.001d0 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 irec = 8 if(reccor) irec = irec + 6*(naddl**2) - 12*naddl + 8 if(nv+irec .ne. nw) stop 50 c c read target points, i.e. points for which Power/Delaunay c cells are to be identified c na = 0 150 continue read (13, *, end = 160) xcor, ycor, zcor na = na + 1 if (na .gt. nlmax) stop 60 icl(na) = 0 xa(na) = xcor ya(na) = ycor za(na) = zcor go to 150 160 continue c if(na.lt.1) stop 70 c c---------------------------------------------------------------------- c write(*,*)' ' write(*,*)'Identification of cells to begin: driver ', * 'subroutine pwrtet will' write(*,*)'be called from the ', * 'main routine and target points will be processed.' write(*,*)' ' write(*,*)'Please wait ...' write(*,*)' ' c c---------------------------------------------------------------------- c c call pwrtet to identify cells c call pwrtet(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2, * xa, ya, za, icon, is, icl, id, nv, nw, nt, na, nb, * nmax, nvmax, nlmax, wlenx, wleny, wlenz, wlenw, naddl, * icfig, iwfig, epz, delaun, flphis, artfcl, reccor) c c---------------------------------------------------------------------- c write(*,*)' ' write(*,*)'(Back to the main routine).' c c---------------------------------------------------------------------- c c save cells target information c write(*,*)' ' write(*,*)'Saving cells target information in file ', * 'targetpoly ...' write (14,*) nb write (14,500) (icl(i), i=1,nb) 500 format (8i10) c c---------------------------------------------------------------------- c write(*,*)' ' stop end *PWRTET c********************************************************************** c c Driver subroutine of Fortran 77 program PWRTET for doing 3D c nearest point searches. It identifies Power cells in the c Power diagram of a set S of 3-d points that contain points in c a set R, i.e. for each point in R finds a point in S that is c as Power close to the point in R as any other point in S. c c If no weights are present then Voronoi cells in the c Voronoi diagram of S that contain points in R are identified, c i.e. for each point in R finds a point in S that is as close c to the point in R as any other point in S. c c It is assumed that a Regular/Delaunay tetrahedralization for S c has already been computed using existing program regtet and that c the output tetrahedron list produced by regtet contains real c as well as artificial tetrahedra, i. e. logical variable artfcl c was set to .true. during the execution of regtet. The output c tetrahedron list from regtet is then used as input for pwrtet. c It is noted that the number of significant figures of decimal c part of point coordinates for points in R is the same as for c points in S which is inherited from the run of regtet. c c Array icl will contain ouput information. If point r is the mth c point in R and n is the absolute value of icl(m) then the Power/ c Voronoi cell of the nth point in S contains point r. c If icl(m) is positive then point r is in the covex hull of S. c In the exterior otherwise. c c The program is based on an algorithm for constructing Regular c tetrahedralizations with incremental topological flipping. c However no actual flippings take place. Given a point in R, c if weights are present a weight is assigned to the point so that c the Power cell of the point in the Power diagram of S together c with the point contains the point. The program then goes through c the motions of inserting the point into the Regular/Delaunay c tetrahedralization for S without actually doing it. This way a c subset of points in S is identified that contains what would be c the Voronoi/Power neighbors of the point in the Voronoi/Power c diagram of S together with the point. The desired point in S c is then a point in this subset closest to the the point in R c in the Power sense. If weights are present care is taken in c choosing the weight assigned to the point in R so that it is as c small as possible. c c c Author: Javier Bernal c c********************************************************************** c subroutine pwrtet(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, * iw2, xa, ya, za, icon, is, icl, id, nv, nw, * nt, na, nb, nmax, nvmax, nlmax, wlenx, wleny, * wlenz, wlenw, naddl, icfig, iwfig, epz, * delaun, flphis, artfcl, reccor) c integer nmax, nvmax, nlmax double precision x(*), y(*), z(*), w(*) double precision xa(*), ya(*), za(*) integer ix(*), iy(*), iz(*), iw(*) integer ix2(*), iy2(*), iz2(*), iw2(*) integer icon(8,*), is(*), icl(*), id(*) integer nv, nw, nt, na, nb, naddl, icfig, iwfig double precision wlenx, wleny, wlenz, wlenw, epz double precision dscle, dfull, dfell logical delaun, flphis, artfcl, reccor c integer isclp(2), isclw(2), isclr(2) double precision xmin, xmax, ymin, ymax, zmin, zmax, wmin, wmax integer mhalf, mfull, i, iftal, no integer irec, irec1, nv1, nu, mxlook, k c c double precision dmin, dcur c integer lmin, lcur, iviol, ipout c c test for presence of artificial points and rectangular grid c if(.not.artfcl) then write(*,*)' ' write(*,*)'Input tetrahedra do not contain artificial points.' write(*,*)'Program terminated.' stop 80 endif c c initialize Fortran 77 word lengths c mhalf=32768 mfull=1073741824 c c testing parameters and number of input points or vertices c if(nv.lt.1 .or. nw.gt.nmax .or. nt.lt.12 .or. nt.gt.nvmax) stop 90 if(na.lt.1 .or. na.gt.nlmax) stop 100 nb = na c c initialize arrays w and id c if(delaun) then do 80 i = 1, nmax w(i) = 0.0d0 80 continue endif c iftal = 0 do 100 i = 1, nvmax id(i) = 0 100 continue 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 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 150 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) 150 continue c WRITE(*,*)' ' c WRITE(*,*)'WMAX=',WMAX c irec = 8 if(reccor) irec = irec + 6*(naddl**2) - 12*naddl + 8 if(nv+irec .ne. nw) stop 140 irec1 = irec + 1 nv1 = nv nv = nv + irec if(nv .gt. nmax) stop 150 do 190 no = nv, irec1, -1 nu = no - irec x(no) = x(nu) y(no) = y(nu) z(no) = z(nu) w(no) = w(nu) 190 continue c c write(*,*)' ' c write(*,*)'Entering poltri ...' call poltri(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2, * xa, ya, za, icon, is, icl, id, xmin, ymin, zmin, * wmin, xmax, ymax, zmax, wmax, iftal, nv, na, nmax, * nvmax, wlenx, wleny, wlenz, wlenw, naddl, irec, nt, * mxlook, delaun, flphis, icfig, iwfig, mhalf, mfull, * isclp, isclw, isclr, epz, dscle, dfull, dfell) c WRITE(*,*)' ' c WRITE(*,*)'MAXLOOK=',MAXLOOK c write (*,*)' ' c write (*,*)'Leaving poltri ...' c nv = nv1 c c---------------------------------------------------------------------- c c test solution the hard way c c write(*,*)' ' c write(*,*)'Testing solution the hard way ...' c write(*,*)' ' c iviol = 0 c k = nw+1 c if(k.gt.nmax) stop 160 c do 900 i = 1, nb c if(i.le.(i/1000)*1000) c * write(*,*)'Number of target points processed = ',i c x(k) = xa(i) c y(k) = ya(i) c z(k) = za(i) c call intran(x, y, z, ix, iy, iz, ix2, iy2, iz2, k, mfull, c * dscle, dfull, dfell) c lmin = iabs(icl(i))+8 c dmin = (x(lmin)-x(k))**2 + (y(lmin)-y(k))**2 + c * (z(lmin)-z(k))**2 - w(lmin) c do 800 lcur = 9, nw c if(is(lcur).eq.0) go to 800 c dcur = (x(lcur)-x(k))**2 + (y(lcur)-y(k))**2 + c * (z(lcur)-z(k))**2 - w(lcur) c if(dsign(dsqrt(dabs(dcur)),dcur) .ge. c * dsign(dsqrt(dabs(dmin)),dmin)+epz) go to 800 c if(dsign(dsqrt(dabs(dcur)),dcur) .le. c * dsign(dsqrt(dabs(dmin)),dmin)-epz) then c iviol = iviol+1 c go to 900 c endif c call ddsign(ix, iy, iz, iw, ix2, iy2, iz2, iw2, k, lmin, c * lcur, mhalf, mfull, isclp, isclw, isclr, c * delaun, ipout) c if(ipout.eq.-1) then c iviol = iviol+1 c go to 900 c endif c IF(IPOUT.EQ.-1) THEN c WRITE(*,*)'I=',I,' LMIN=',LMIN,' LCUR=',LCUR c WRITE(*,*)'DMIN=',DMIN,' DCUR=',DCUR c WRITE(*,*)'XYZW(LMIN)=',X(LMIN),Y(LMIN),Z(LMIN),W(LMIN) c WRITE(*,*)'XYZW(LCUR)=',X(LCUR),Y(LCUR),Z(LCUR),W(LCUR) c ENDIF c 800 continue c 900 continue c write(*,*)' ' c write(*,*)'# of proximity violations = ',iviol c c update target points c k = nw+1 if(k.gt.nmax) stop 170 do 1000 i = 1, nb x(k) = xa(i) y(k) = ya(i) z(k) = za(i) call intran(x, y, z, ix, iy, iz, ix2, iy2, iz2, k, mfull, * dscle, dfull, dfell) xa(i) = x(k) ya(i) = y(k) za(i) = z(k) 1000 continue c return end *POLTRI c c This subroutine will process target points, etc. c subroutine poltri(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2, * xa, ya, za, icon, is, icl, id, xmin, ymin, zmin, * wmin, xmax, ymax, zmax, wmax, iftal, nv, na, * nmax, nvmax, wlenx, wleny, wlenz, wlenw, naddl, * irec, tetra, mxlook, delaun, flphis, icsfig, * iwsfig, mhalf, mfull, isclp, isclw, isclr, epz, * dscle, dfull, dfell) c double precision x(*), y(*), z(*), w(*) double precision xa(*), ya(*), za(*) integer ix(*), iy(*), iz(*), iw(*) integer ix2(*), iy2(*), iz2(*), iw2(*) integer icon(8,*), is(*), icl(*), id(*) double precision xc(8), yc(8), zc(8) integer ixc(8), iyc(8), izc(8) integer nmax, nvmax integer nv, na, icsfig, iwsfig, mxlook, naddl, irec integer mhalf, mfull, ibsfig, itsfig double precision wlenx, wleny, wlenz, wlenw, epz, wbig, wmen, wman logical delaun, flphis integer isclp(*), isclw(*), isclr(*), tetra double precision xmin, xmax, ymin, ymax, zmin, zmax, wmin, wmax double precision xmin2, xmax2, ymin2, ymax2, zmin2, zmax2 double precision xmin3, xmax3, ymin3, ymax3, zmin3, zmax3 double precision xint, yint, zint, xcor, ycor, zcor double precision xlen, ylen, zlen, wlen, xctr, yctr, zctr double precision xlan, ylan, zlan, wlan double precision xmx, ymx, zmx, xmn, ymn, zmn double precision dscle, dscli, dfull, dfell, dfill double precision decml, r215, deps integer iftal, i, j, i9, icsfi2, irsfig, isclu, isgcl integer issin, k, ka, no, naddm, ng, lmin c integer lman INTEGER NTET REAL WTET c c initialize c mxlook = 0 c c update min and max c xmin2 = xmin xmax2 = xmax ymin2 = ymin ymax2 = ymax zmin2 = zmin zmax2 = zmax do 30 no = 1, na if (xa(no) .gt. xmax) xmax = xa(no) if (xa(no) .lt. xmin) xmin = xa(no) if (ya(no) .gt. ymax) ymax = ya(no) if (ya(no) .lt. ymin) ymin = ya(no) if (za(no) .gt. zmax) zmax = za(no) if (za(no) .lt. zmin) zmin = za(no) 30 continue c c test parameters c if (nv .gt. nmax) stop 210 if (nv .lt. 9) stop 220 if (nvmax .lt. 12) stop 230 c c construct cube c xlen=xmax-xmin ylen=ymax-ymin zlen=zmax-zmin wlan=(dmax1(xlen,ylen,zlen)/2.0d0) + dmax1(wlenx,wleny,wlenz) wlen=wlan + 4.0d0 if(wlen.le.wlan) stop 235 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 xc(1) = - 1.0d0 yc(1) = - 1.0d0 zc(1) = + 1.0d0 c xc(2) = - 1.0d0 yc(2) = + 1.0d0 zc(2) = + 1.0d0 c xc(3) = + 1.0d0 yc(3) = + 1.0d0 zc(3) = + 1.0d0 c xc(4) = + 1.0d0 yc(4) = - 1.0d0 zc(4) = + 1.0d0 c xc(5) = - 1.0d0 yc(5) = - 1.0d0 zc(5) = - 1.0d0 c xc(6) = - 1.0d0 yc(6) = + 1.0d0 zc(6) = - 1.0d0 c xc(7) = + 1.0d0 yc(7) = + 1.0d0 zc(7) = - 1.0d0 c xc(8) = + 1.0d0 yc(8) = - 1.0d0 zc(8) = - 1.0d0 c c compute corners of cube c do 50 i=1,8 x(i)=xctr+wlen*xc(i) y(i)=yctr+wlen*yc(i) z(i)=zctr+wlen*zc(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 if(irec.eq.8) go to 155 c xmin3 = xmin xmax3 = xmax ymin3 = ymin ymax3 = ymax zmin3 = zmin zmax3 = zmax xmin = xmin2 xmax = xmax2 ymin = ymin2 ymax = ymax2 zmin = zmin2 zmax = zmax2 c c compute corners of rectangular polyhedron c x( 9) = xmin - wlenx y( 9) = ymin - wleny z( 9) = zmax + wlenz c x(10) = xmin - wlenx y(10) = ymax + wleny z(10) = zmax + wlenz c x(11) = xmax + wlenx y(11) = ymax + wleny z(11) = zmax + wlenz c x(12) = xmax + wlenx y(12) = ymin - wleny z(12) = zmax + wlenz c x(13) = xmin - wlenx y(13) = ymin - wleny z(13) = zmin - wlenz c x(14) = xmin - wlenx y(14) = ymax + wleny z(14) = zmin - wlenz c x(15) = xmax + wlenx y(15) = ymax + wleny z(15) = zmin - wlenz c x(16) = xmax + wlenx y(16) = ymin - wleny z(16) = zmin - wlenz c do 60 i=9,16 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 255 60 continue if(x(1).ge.x( 9) .or. y(1).ge.y( 9) .or. z(1).le.z( 9) .or. * x(2).ge.x(10) .or. y(2).le.y(10) .or. z(2).le.z(10) .or. * x(3).le.x(11) .or. y(3).le.y(11) .or. z(3).le.z(11) .or. * x(4).le.x(12) .or. y(4).ge.y(12) .or. z(4).le.z(12) .or. * x(5).ge.x(13) .or. y(5).ge.y(13) .or. z(5).ge.z(13) .or. * x(6).ge.x(14) .or. y(6).le.y(14) .or. z(6).ge.z(14) .or. * x(7).le.x(15) .or. y(7).le.y(15) .or. z(7).ge.z(15) .or. * x(8).le.x(16) .or. y(8).ge.y(16) .or. z(8).ge.z(16)) stop 260 c xmin = xmin - wlenx ymin = ymin - wleny zmin = zmin - wlenz xmax = xmax + wlenx ymax = ymax + wleny zmax = zmax + wlenz if (xmax .gt. xmax3) xmax3 = xmax if (xmin .lt. xmin3) xmin3 = xmin if (ymax .gt. ymax3) ymax3 = ymax if (ymin .lt. ymin3) ymin3 = ymin if (zmax .gt. zmax3) zmax3 = zmax if (zmin .lt. zmin3) zmin3 = zmin c if(naddl.eq.2) go to 153 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 = 16 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 265 c 153 continue xmin = xmin3 xmax = xmax3 ymin = ymin3 ymax = ymax3 zmin = zmin3 zmax = zmax3 c 155 continue c c find i9 c do 160 i = 9, nv if(is(i).gt.0) go to 165 160 continue stop 270 165 continue i9 = i c wmen = wmin if(.not.delaun .and. irec.gt.8) then wmen = wmen - wlenw do 167 i=9,irec w(i) = wmen 167 continue endif if(wmen.lt.wmin) wmin = wmen if(wmen.gt.wmax) wmax = wmen wman = wmin - 10.0d0 do 170 i=1,8 w(i) = wman 170 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(*,*)'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 273 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 275 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 278 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 280 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 300 c c test lengths of x, y, z-coordinates, shift and make them integers c dfull = dble(mfull) dfell=dfull/dscle r215 = dble(mhalf) deps = dble(0.9) if(dscle.lt.deps) stop 305 do 235 i = 1, nv call intran(x, y, z, ix, iy, iz, ix2, iy2, iz2, i, mfull, * dscle, dfull, dfell) 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(nv+1 .gt. nmax) stop 310 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 c get cube corner directions in their integer form c do 320 i = 1,8 ixc(i) = idnint(xc(i)) iyc(i) = idnint(yc(i)) izc(i) = idnint(zc(i)) 320 continue c c find power cells that contain target points c write(*,*)'Identifying cells that contain target points ...' write(*,*)' ' k = nv+1 issin = i9 lmin = issin NTET = 0 do 380 ka = 1, na if(ka.le.(ka/1000)*1000) * write(*,*)'Number of target points processed = ',ka x(k) = xa(ka) y(k) = ya(ka) z(k) = za(ka) call intran(x, y, z, ix, iy, iz, ix2, iy2, iz2, k, mfull, * dscle, dfull, dfell) call pntins(x, y, z, w, ix, iy, iz, iw, ix2, iy2, iz2, iw2, * icon, is, icl, id, k, ka, tetra, mxlook, * ixc, iyc, izc, issin, lmin, iftal, delaun, * flphis, mhalf, mfull, isclp, isclw, isclr, * dscle, dscli, dfill, dfull, r215, deps, epz, NTET) 380 continue WTET = REAL(NTET)/REAL(NA) WRITE(*,*)' ' WRITE(*,*)'AVERAGE # OF CHECKED TETRAHEDRA PER TARGET PNT =',WTET c return end *PNTINS c c This subroutine will find location of new point c c This routine also calls routine 'sphere' for the purpose of optimizing c for the locally Regular/Delaunay property c subroutine pntins(xi, yi, zi, wi, x, y, z, w, x2, y2, z2, w2, * icon, is, icl, id, k, ka, tetra, mxlook, * xc, yc, zc, issin, lmin, iftal, delaun, * flphis, mhalf, mfull, isclp, isclw, isclr, * dscle, dscli, dfill, dfull, r215, deps, epz, * NTET) c double precision xi(*), yi(*), zi(*), wi(*) integer x(*), y(*), z(*), w(*) integer x2(*), y2(*), z2(*), w2(*) integer icon(8,*), is(*), icl(*), id(*) integer xc(*), yc(*), zc(*) integer k, ka, mxlook, issin, lmin, iftal integer isclp(*), isclw(*), isclr(*), curr, side1, side2, tetra integer site1, site2, ihull double precision dscle, dscli, dfill, dfull, r215, deps, epz logical delaun, flphis integer mhalf, mfull, itype, look, i INTEGER NTET c side1 = 0 side2 = 0 if(flphis) then itype = -1 look = 1 call gettet(itype, k, xi, yi, zi, x, y, z, x2, y2, z2, icon, * curr, side1, side2, xc, yc, zc, mhalf, mfull, * isclp, epz) c 50 continue if (icon(5,curr) .lt. 0) then call lkdown(icon, curr, xi, yi, zi, x, y, z, x2, y2, z2, * itype, k, side1, side2, xc, yc, zc, mhalf, * mfull, isclp, epz) look = look + 1 if(itype.eq.-1) stop 410 goto 50 endif if (look .gt. mxlook) mxlook = look else call shishk(xi, yi, zi, x, y, z, x2, y2, z2, is, icon, id, * issin, k, side1, side2, curr, iftal, itype, tetra, * xc, yc, zc, mhalf, mfull, isclp, epz) if(itype.eq.-1) stop 420 endif c ihull = 0 if(itype .eq. 1) then site1 = icon(side1+4,curr) if(site1.le.8) stop 430 elseif(itype .eq. 2) then if(icon(5,curr).le.8 .or. icon(6,curr).le.8 .or. * icon(7,curr).le.8 .or. icon(8,curr).le.8) ihull = 1 elseif(itype .eq. 3) then site1 = icon(side1+4,curr) site2 = icon(side2+4,curr) call reordr(icon, site1, site2, curr) site1 = icon(7,curr) site2 = icon(8,curr) if(site1.le.8 .or. site2.le.8) ihull = 1 elseif(itype .eq. 4) then site1 = icon(side1+4,curr) call sitord(icon, site1, curr) if(icon(6,curr).le.8 .or. icon(7,curr).le.8 .or. * icon(8,curr).le.8) ihull = 1 else stop 440 endif c if(delaun .and. itype.eq.1) then lmin = site1 icl(ka) = lmin-8 issin = site1 go to 1000 endif c do 100 i = 5, 8 issin = icon(i,curr) if(issin.gt.8) go to 120 100 continue stop 450 120 continue c c optimize for Regular/Delaunay property c c write(*,*)'Entering sphere ...' call sphere(icon, icl, id, xi, yi, zi, wi, x, y, z, w, x2, y2, * z2, w2, k, ka, curr, tetra, iftal, xc, yc, zc, * delaun, lmin, mhalf, mfull, isclp, isclw, isclr, * dscle, dscli, dfill, dfull, r215, deps, epz, NTET) if(ihull.eq.1) icl(ka) = -icl(ka) c 1000 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, i, * mfull, dscle, dfull, dfill) c double precision x(*), y(*), z(*) integer ix(*), iy(*), iz(*) integer ix2(*), iy2(*), iz2(*) integer i, mfull double precision dscle, dfull, dfill, decml c c test lengths of x, y, z-coordinates, shift and make them integers c 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 *GETTET c c This subroutine will test each of the 1st 12 tetrahedra to find c where new point is located. It'll do so by calling 'gette2'. c subroutine gettet(itype, k, xi, yi, zi, x, y, z, x2, y2, z2, icon, * curr, side1, side2, xc, yc, zc, mhalf, mfull, * isclp, epz) c double precision xi(*), yi(*), zi(*) integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*) integer icon(8,*), curr, a, b, c, d, side1, side2, flag integer isclp(*), mhalf, mfull double precision epz integer itype, k c do 100 curr = 1, 12 a = iabs(icon(5, curr)) b = icon(6, curr) c = icon(7, curr) d = icon(8, curr) c flag = icon(5,curr) icon(5,curr)=a c call vrtord(icon, curr, a, b, c, d) if(flag.lt.0) icon(5,curr)=-a if(a.le.8.or.b.gt.8.or.c.gt.8.or.d.gt.8) stop 680 c call gette2(a, b, c, d, xi, yi, zi, x, y, z, x2, y2, z2, * itype, k, side1, side2, flag, xc, yc, zc, * mhalf, mfull, isclp, epz) if (itype .ne. -1) goto 200 c 100 continue stop 690 200 continue return end *LKDOWN c c This subroutine will traverse thru children of curr until point c is found to be in one of them c subroutine lkdown(icon, curr, xi, yi, zi, x, y, z, x2, y2, z2, * itype, k, side1, side2, xc, yc, zc, mhalf, * mfull, isclp, epz) c double precision xi(*), yi(*), zi(*) integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer icon(8,*), xc(*), yc(*), zc(*) double precision epz integer isclp(*), itype, k, mhalf, mfull, i, newcur integer curr, side1, side2, a, b, c, d, flag c c test children of current tetrahedron c itype = -1 do 100 i = 1, 4 newcur = icon(i,curr) if (newcur .le. 0) goto 100 c a = iabs(icon(5,newcur)) b = icon(6,newcur) c = icon(7,newcur) d = icon(8,newcur) c flag = icon(5,newcur) icon(5,newcur)=a c call vrtord(icon, newcur, a, b, c, d) if(flag.lt.0) icon(5,newcur)=-a if(a.le.8) stop 710 c call gette2(a, b, c, d, xi, yi, zi, x, y, z, x2, y2, z2, * itype, k, side1, side2, flag, xc, yc, zc, * mhalf, mfull, isclp, epz) if (itype .eq. -1) goto 100 curr = newcur goto 1000 100 continue c 1000 continue return end *GETTE2 c c This subroutine will check for each tetra, if the point is equal to an c existing vertex, inside (interior, edge, side), or outside curr tetra. c subroutine gette2(a, b, c, d, xi, yi, zi, x, y, z, x2, y2, z2, * itype, k, side1, side2, flag, xc, yc, zc, * mhalf, mfull, isclp, epz) c double precision xi(*), yi(*), zi(*) integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer xc(*), yc(*), zc(*) double precision epz integer isclp(*), itype, k, mhalf, mfull, ifn, ipout integer iside(4), a, b, c, d, side1, side2, flag c c determine position of point k relative to facets of tetrahedron c if(b.le.8 .or. c.le.8 .or. d.le.8) go to 100 c call irsign(xi, yi, zi, x, y, z, x2, y2, z2, k, b, d, c, * mhalf, mfull, isclp, epz, ipout) iside(1) = ipout 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 c point k is not in tetrahedron c 50 continue if(iside(1).lt.0 .or. iside(2).lt.0 .or. iside(3).lt.0 .or. * iside(4).lt.0) go to 1000 c if(flag.lt.0) then itype = 0 go to 1000 endif c c point k 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 805 c c point k 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 k 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 k 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 807 endif go to 1000 c c there is at least one artificial vertex c 100 continue if(b.le.8) then iside(1) = 1 go to 120 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, ipout) iside(1) = ipout if(iside(1).ne.0) go to 120 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * d, c, k, b, mhalf, mfull, isclp, ipout) iside(1) = ipout go to 120 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, ipout) iside(1) = ipout else ifn = 1 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * b, d, k, c, mhalf, mfull, isclp, ifn, ipout) iside(1) = ipout endif if(iside(1).ne.0) go to 120 call ipsign(x, y, z, x2, y2, z2, b, d, c, k, mhalf, * mfull, isclp, ipout) iside(1) = ipout 120 continue c if(c.le.8.and.d.le.8) then call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * c, d, k, a, mhalf, mfull, isclp, ipout) iside(2) = ipout if(iside(2).ne.0) go to 140 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * c, d, k, a, mhalf, mfull, isclp, ipout) iside(2) = ipout go to 140 elseif(c.le.8) then ifn = 0 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * a, d, k, c, mhalf, mfull, isclp, ifn, ipout) iside(2) = ipout else ifn = 1 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * a, c, k, d, mhalf, mfull, isclp, ifn, ipout) iside(2) = ipout endif if(iside(2).ne.0) go to 140 call ipsign(x, y, z, x2, y2, z2, a, c, d, k, mhalf, * mfull, isclp, ipout) iside(2) = ipout 140 continue c if(d.gt.8) then call ipsign(x, y, z, x2, y2, z2, a, d, b, k, mhalf, * mfull, isclp, ipout) iside(3) = ipout go to 160 elseif(b.le.8) then call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * d, b, k, a, mhalf, mfull, isclp, ipout) iside(3) = ipout if(iside(3).ne.0) go to 160 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * d, b, k, a, mhalf, mfull, isclp, ipout) iside(3) = ipout go to 160 else ifn = 0 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * a, b, k, d, mhalf, mfull, isclp, ifn, ipout) iside(3) = ipout endif if(iside(3).ne.0) go to 160 call ipsign(x, y, z, x2, y2, z2, a, d, b, k, mhalf, * mfull, isclp, ipout) iside(3) = ipout 160 continue c if(c.gt.8) then call ipsign(x, y, z, x2, y2, z2, a, b, c, k, mhalf, * mfull, isclp, ipout) iside(4) = ipout go to 180 elseif(b.le.8) then call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * b, c, k, a, mhalf, mfull, isclp, ipout) iside(4) = ipout if(iside(4).ne.0) go to 180 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * b, c, k, a, mhalf, mfull, isclp, ipout) iside(4) = ipout go to 180 else ifn = 1 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * a, b, k, c, mhalf, mfull, isclp, ifn, ipout) iside(4) = ipout endif if(iside(4).ne.0) go to 180 call ipsign(x, y, z, x2, y2, z2, a, b, c, k, mhalf, * mfull, isclp, ipout) iside(4) = ipout 180 continue go to 50 c 1000 continue return end *SHISHK c c shishkebab routines c c subroutine shishk to - c c move from a vertex in the tetrahedralization to a tetrahedron c that contains a point, and identify the type of location of c the point with respect to the tetrahedron c subroutine shishk(xi, yi, zi, x, y, z, x2, y2, z2, is, icon, id, * ileft, k, side1, side2, iscur, iftal, itype, * ivnxt, xc, yc, zc, mhalf, mfull, isclp, epz) c double precision xi(*), yi(*), zi(*) integer x(*), y(*), z(*) integer x2(*), y2(*), z2(*) integer is(*), icon(8,*), id(*) double precision epz integer xc(*), yc(*), zc(*) integer isclp(*), ileft, k, iftal, itype, ivnxt, mhalf, mfull integer a, b, c, d, side1, side2, site0, site1 integer iscur, isini, imist, isadj, ilift, islst, i c c reinitialize array id if necessary c if(iftal.gt.10000000) then iftal = 0 do 50 i = 1, ivnxt id(i) = 0 50 continue endif c if(ileft .le. 8) stop 911 a = ileft 100 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 912 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 call fctest(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, itype, k, * imist, a, b, c, d, side1, side2, mhalf, mfull, isclp, * epz) if(itype .gt. 0) go to 9000 if(itype .eq. 0) go to 500 if(itype .eq. -2) go to 1100 if(itype .eq. -3) then a = imist go to 100 elseif(itype .eq. -4) then site0 = a site1 = imist go to 2000 else stop 913 endif c c obtain next tetrahedron with point a as a vertex c 500 continue isadj = iabs(icon(2,iscur)) if(isadj.le.0.or.isadj.gt.ivnxt) stop 914 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.or.isadj.gt.ivnxt) stop 915 if(id(isadj) .eq. iftal) go to 700 ilift = icon(6,iscur) go to 900 700 continue isadj = iabs(icon(4,iscur)) if(isadj.le.0.or.isadj.gt.ivnxt) stop 916 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 917 go to 1000 else stop 918 endif 800 continue 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 stop 919 c c obtain next tetrahedron along line segment as it crosses a facet c 1100 continue islst = iscur isadj = iabs(icon(1,iscur)) if(isadj.le.0.or.isadj.gt.ivnxt) stop 921 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 922 endif c c obtain opposite facet of tetrahedron intersected by line c segment c call fcfind(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, itype, * ileft, k, ilift, imist, b, c, d, 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 100 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 923 endif c c obtain next tetrahedron along line segment as it crosses an edge c 2000 continue call fcedge(x, y, z, x2, y2, z2, xc, yc, zc, itype, ileft, k, * icon, iscur, imist, ivnxt, site0, site1, side1, * side2, mhalf, mfull, isclp) if(itype .gt. 0) 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 100 elseif(itype.eq.-4) then if(imist.eq.site1)then site0 = icon(7,iscur) elseif(imist.eq.site0)then site0 = site1 endif site1 = imist go to 2000 else stop 924 endif c 9000 continue 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 *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, 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 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.8 .or. c.le.8 .or. d.le.8) 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 931 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 932 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, * 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 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.8 .or. c.le.8 .or. d.le.8 .or. ilift.le.8) 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 200 go to 50 c c k is not in tetrahedron c c determine position of ilift with repect to current situation c 200 continue if(b.le.8 .or. c.le.8 .or. d.le.8 .or. ilift.le.8) go to 300 c 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 c c ilift, ileft, b, d, c, form a strictly convex set c 400 continue if(b.le.8 .or. c.le.8 .or. d.le.8 .or. ilift.le.8) go to 500 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 951 endif go to 1000 c 700 continue if(idut1.le.0 .and. idut2.le.0 .and. idut3.le.0) stop 952 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.8 .and. ilift.gt.8) 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.8 .and. ilift.gt.8) 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.8 .and. ilift.gt.8) 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, mhalf, mfull, isclp) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*) integer icon(8,*) integer isclp(*), mhalf, mfull integer itype, ileft, k, iscur, imist, ivnxt integer iside(4), site0, site1, site2, site3, side1, side2 integer isnow, idut, iasign, idot0, ipout c c find intersecting facet c 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.or.isnow.gt.ivnxt) stop 961 if(isnow.eq.iscur) stop 963 call reordr(icon, site0, site1, isnow) site3 = icon(8,isnow) if(site1.gt.8 .and. site2.gt.8 .and. site3.gt.8) 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(idut.ge.0) go to 400 site0 = site3 isnow = iabs(icon(1,isnow)) go to 300 c 400 continue iscur = isnow c c determine whether point k is in tetrahedron c if(site0.le.8 .or. site1.le.8 .or. site2.le.8 .or. site3.le.8) * 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) = 1 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.8 .and. site3.gt.8) 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(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 c 1000 continue return end *FCSIGN c c This subroutine will determine where point k is relative to c facet with vertices b, c, d. IF b, c, d are not artificial it c will compute either the perpendicular distance or the shortest c distance from the facet to the point. c subroutine fcsign(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, * isyde, k, b, c, d, mhalf, mfull, isclp, epz, * dist, r215, deps, dscle, idype, delaun) c double precision xi(*), yi(*), zi(*) integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*) integer b, c, d, isyde, k, mhalf, mfull, isclp(*), ipossi, idype double precision epz, dist, r215, deps, dscle logical delaun c if(b.le.8 .or. c.le.8 .or. d.le.8) go to 100 c if(delaun .and. idype.eq.-1) then call tstnce(xi, yi, zi, b, d, c, epz, k, dist, ipossi, idype) else call dstanc(xi, yi, zi, b, d, c, epz, k, dist, ipossi) endif if(ipossi.eq.0) then isyde = 1 if(dist.lt.0.0d0) isyde = -1 else if(delaun .and. idype.eq.-1) then call tstexa(x, y, z, x2, y2, z2, b, d, c, k, mhalf, mfull, * isclp, r215, deps, dscle, dist, isyde, idype) else call dstexa(x, y, z, x2, y2, z2, b, d, c, k, mhalf, mfull, * isclp, r215, deps, dscle, dist, isyde) endif endif 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, isyde) 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(*), itype, side1, side2 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 971 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 972 endif c 1000 continue 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 *SPHERE c c This subroutine will optimize locally at point k for the c Regular/Delaunay property c subroutine sphere(icon, icl, idl, xi, yi, zi, wi, x, y, z, w, * x2, y2, z2, w2, k, ka, isini, tetra, iftal, * xc, yc, zc, delaun, lmin, mhalf, mfull, isclp, * isclw, isclr, dscle, dscli, dfill, dfull, * r215, deps, epz, NTET) c double precision xi(*), yi(*), zi(*), wi(*) integer x(*), y(*), z(*), w(*), x2(*), y2(*), z2(*), w2(*) integer icon(8,*), icl(*), idl(*) double precision tdist, xctr, yctr, zctr, dmin, dcur, dmun double precision dscle, dscli, dfill, dfull double precision r215, deps, epz, dist, wmax, decml integer iftal, mhalf, mfull, idype, iftai integer xc(*), yc(*), zc(*) integer ikon(8,1), opvert, tetra, a, b, c, d, oddsid, ipout integer isclp(*), isclw(*), isclr(*), iside(4), sidist(4) logical delaun integer k, ka, i, ipossi, itide, j, ib, ic, id, ifn, istt, isyde integer isodd, isite, isini, islst, iscur, indx, ix integer lmun, lmin, lcur INTEGER NTET c c reinitialize array idl if necessary c if(iftal.gt.10000000) then iftal = 0 do 50 i = 1, tetra idl(i) = 0 50 continue endif c iftai = iftal dmin = (xi(lmin) - xi(k))**2 + (yi(lmin) - yi(k))**2 + * (zi(lmin) - zi(k))**2 - wi(lmin) wmax = -dmin dmun = dmin lmun = lmin c NTET = NTET+1 do 150 i = 5, 8 lcur = icon(i,isini) if(lcur.le.8) go to 150 dcur = (xi(lcur) - xi(k))**2 + (yi(lcur) - yi(k))**2 + * (zi(lcur) - zi(k))**2 - wi(lcur) if(dsign(dsqrt(dabs(dcur)),dcur) .ge. dsign(dsqrt(dabs(dmin)), * dmin)+epz) go to 150 if(dsign(dsqrt(dabs(dcur)),dcur) .le. dsign(dsqrt(dabs(dmin)), * dmin)-epz) then lmin = lcur dmin = dcur wmax = -dcur go to 150 endif call ddsign(x, y, z, w, x2, y2, z2, w2, k, lmin, lcur, * mhalf, mfull, isclp, isclw, isclr, delaun, ipout) if(ipout.eq.-1) then lmin = lcur dmin = dcur wmax = -dcur endif 150 continue c if(delaun) go to 180 160 continue wi(k) = wmax + epz w2(k) = 0 if(dabs(wi(k)).lt.dfill) then w(k) = idnint(dscli*wi(k)) if(iabs(w(k)).lt.mfull) then wi(k) = dble(w(k))/dscli go to 180 endif endif if(dabs(wi(k)).ge.dfull) stop 1165 w(k) = idint(wi(k)) if(iabs(w(k)).ge.mfull) stop 1170 decml = (wi(k) - dint(wi(k)))*dscli if(dabs(decml).ge.dfull) stop 1172 w2(k) = idnint(decml) if(iabs(w2(k)).ge.mfull) stop 1175 if(iabs(w2(k)).eq.0) then wi(k) = dble(w(k)) w2(k) = mfull else wi(k) = dble(w(k)) + (dble(w2(k))/dscli) endif 180 continue c a = k iftal = iftal+1 idl(isini) = iftal indx = 1 islst = isini iscur = icon(1,isini) if(iscur.eq.0) go to 1100 b = icon(6,isini) c = icon(7,isini) d = icon(8,isini) c c reorder iscur relative to b and c, and test c 200 continue ix = max0(b,c,d) if(b .eq. ix) go to 220 if(c .eq. ix) then ix = b b = c c = d d = ix else ix = b b = d d = c c = ix endif 220 continue if(c.gt.b .or. d.gt.b) stop 1210 call reordr(icon, b, c, iscur) if(icon(7,iscur) .ne. d) stop 1220 if(icon(4,iscur) .ne. islst) stop 1230 opvert = icon(8,iscur) if(delaun) go to 230 if(idl(iscur).le.iftai) go to 225 idl(iscur) = iftal islst = iscur indx = 1 iscur = icon(1,islst) if(iscur.eq.0) go to 1100 b = icon(6,islst) c = icon(7,islst) d = icon(8,islst) if(idl(iscur) .ne. iftal) go to 200 go to 1100 225 continue if(opvert.le.8) go to 230 dcur = (xi(opvert) - xi(k))**2 + (yi(opvert) - yi(k))**2 + * (zi(opvert) - zi(k))**2 - wi(opvert) if(wmax .ge. -dcur) go to 230 wmax = -dcur go to 160 c 230 continue NTET = NTET+1 idype = -1 call fcsign(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, * isyde, k, b, c, d, mhalf, mfull, isclp, * epz, dist, r215, deps, dscle, idype, delaun) if(isyde.lt.0) go to 1100 if(isyde.gt.0) go to 235 idype = 0 call fcsign(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, * isyde, k, b, c, opvert, mhalf, mfull, isclp, * epz, dist, r215, deps, dscle, idype, delaun) if(isyde.lt.0) go to 1100 call fcsign(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, * isyde, k, c, d, opvert, mhalf, mfull, isclp, * epz, dist, r215, deps, dscle, idype, delaun) if(isyde.lt.0) go to 1100 call fcsign(xi, yi, zi, x, y, z, x2, y2, z2, xc, yc, zc, * isyde, k, d, b, opvert, mhalf, mfull, isclp, * epz, dist, r215, deps, dscle, idype, delaun) if(isyde.lt.0) go to 1100 go to 990 235 continue if(.not.delaun .or. b.le.8 .or. c.le.8 .or. d.le.8) go to 240 c dist = (dist-epz)**2 c if(dist .gt. dmin) go to 1000 if(dist-epz .gt. dsqrt(dabs(dmin))) go to 1000 if(idype .eq. lmin) go to 1000 240 continue ikon(5,1) = a ikon(6,1) = b ikon(7,1) = c ikon(8,1) = d c if(opvert.le.8.or.c.le.8.or.d.le.8) go to 300 c c test for Regular/Delaunay configuration c call ctrad(xi, yi, zi, wi, xctr, yctr, zctr, b, c, d, opvert, * epz, delaun, ipossi) if(ipossi.eq.1) go to 250 call bisphr(xi, yi, zi, wi, a, b, epz, xctr, yctr, zctr, tdist, * delaun, ipossi) if(ipossi.eq.1) go to 250 if(tdist.le.0.0d0) go to 1000 go to 990 250 continue call iqsign(x, y, z, w, x2, y2, z2, w2, a, b, c, d, opvert, * mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 go to 990 c 300 continue if(opvert.le.8.and.c.gt.8.and.d.gt.8) go to 1000 c c determine signs of distances from opvert to facets of c tetrahedron with vertices a, b, c, d c if(opvert.gt.8)then if(c.le.8.and.d.le.8)then call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * d, c, a, opvert, mhalf, mfull, isclp, ipout) iside(2) = ipout if(iside(2).ne.0) go to 310 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * d, c, a, opvert, mhalf, mfull, isclp, ipout) iside(2) = ipout go to 310 endif call vrtarr(c,d,opvert,ib,ic,id) if(ic.gt.8)then ifn = 0 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * ib, ic, a, id, mhalf, mfull, isclp, ifn, ipout) iside(2) = ipout else ifn = 1 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * ib, id, a, ic, mhalf, mfull, isclp, ifn, ipout) iside(2) = ipout endif if(iside(2).ne.0) go to 310 call ipsign(x, y, z, x2, y2, z2, ib, id, ic, a, mhalf, * mfull, isclp, ipout) iside(2) = ipout 310 continue c call vrtarr(b,opvert,d,ib,ic,id) if(d.gt.8)then call ipsign(x, y, z, x2, y2, z2, ib, id, ic, a, mhalf, * mfull, isclp, ipout) iside(3) = ipout go to 320 endif if(ic.gt.8)then ifn = 0 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * ib, ic, a, id, mhalf, mfull, isclp, ifn, ipout) iside(3) = ipout else ifn = 1 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * ib, id, a, ic, mhalf, mfull, isclp, ifn, ipout) iside(3) = ipout endif if(iside(3).ne.0) go to 320 call ipsign(x, y, z, x2, y2, z2, ib, id, ic, a, mhalf, * mfull, isclp, ipout) iside(3) = ipout 320 continue c call vrtarr(b,c,opvert,ib,ic,id) if(c.gt.8)then call ipsign(x, y, z, x2, y2, z2, ib, id, ic, a, mhalf, * mfull, isclp, ipout) iside(4) = ipout go to 330 endif if(ic.gt.8)then ifn = 0 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * ib, ic, a, id, mhalf, mfull, isclp, ifn, ipout) iside(4) = ipout else ifn = 1 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * ib, id, a, ic, mhalf, mfull, isclp, ifn, ipout) iside(4) = ipout endif if(iside(4).ne.0) go to 330 call ipsign(x, y, z, x2, y2, z2, ib, id, ic, a, mhalf, * mfull, isclp, ipout) iside(4) = ipout 330 continue c else if(c.le.8.and.d.le.8) then iside(2) = 1 elseif(c.gt.8) then call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * opvert, d, a, c, mhalf, mfull, isclp, ipout) iside(2) = ipout if(iside(2).ne.0) go to 340 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * opvert, d, a, c, mhalf, mfull, isclp, ipout) iside(2) = ipout else call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * c, opvert, a, d, mhalf, mfull, isclp, ipout) iside(2) = ipout if(iside(2).ne.0) go to 340 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * c, opvert, a, d, mhalf, mfull, isclp, ipout) iside(2) = ipout endif 340 continue c if(d.gt.8)then ifn = 1 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * b, d, a, opvert, mhalf, mfull, isclp, ifn, ipout) iside(3) = ipout else call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * d, opvert, a, b, mhalf, mfull, isclp, ipout) iside(3) = ipout if(iside(3).ne.0) go to 350 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * d, opvert, a, b, mhalf, mfull, isclp, ipout) iside(3) = ipout go to 350 endif if(iside(3).ne.0) go to 350 call ipsign(x, y, z, x2, y2, z2, b, d, opvert, a, mhalf, * mfull, isclp, ipout) iside(3) = ipout 350 continue c if(c.gt.8)then ifn = 0 call ipsig3(x, y, z, x2, y2, z2, xc, yc, zc, * b, c, a, opvert, mhalf, mfull, isclp, ifn, ipout) iside(4) = ipout else call ipsig4(x, y, z, x2, y2, z2, xc, yc, zc, * opvert, c, a, b, mhalf, mfull, isclp, ipout) iside(4) = ipout if(iside(4).ne.0) go to 360 call ipsig6(x, y, z, x2, y2, z2, xc, yc, zc, * opvert, c, a, b, mhalf, mfull, isclp, ipout) iside(4) = ipout go to 360 endif if(iside(4).ne.0) go to 360 call ipsign(x, y, z, x2, y2, z2, b, opvert, c, a, mhalf, * mfull, isclp, ipout) iside(4) = ipout 360 continue endif c c set sidist array c do 400 j = 2, 4 if(iside(j) .gt. 0) then sidist(j) = 0 elseif(iside(j) .lt. 0) then sidist(j) = -1 else sidist(j) = 1 endif 400 continue c c check for Regular/Delaunay property c if ((sidist(2). eq. 0) .and. (sidist(3). eq .0) .and. * (sidist(4). eq. 0)) then if(opvert.gt.8) go to 990 if(c.le.8.and.d.le.8)then call ipsig2(x, y, z, x2, y2, z2, xc, yc, zc, opvert, d, * opvert, c, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) elseif(c.le.8)then call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, c, opvert, * b, d, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) else call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, d, opvert, * c, b, a, c, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) endif if(itide .ge. 0) go to 1000 go to 990 elseif ((sidist(2).eq.-1) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.0)) then if(c.le.8.and.d.le.8) stop 1660 if(opvert.gt.8) go to 990 if(c.le.8)then call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, c, opvert, * b, d, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) else call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, d, opvert, * c, b, a, c, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) endif if(itide .ge. 0) go to 1000 go to 990 elseif ((sidist(2).eq.0) .and. (sidist(3).eq.-1) .and. * (sidist(4).eq.0)) then if(d.gt.8) go to 1000 if(opvert.gt.8.and.c.gt.8) go to 990 if(opvert.gt.8)then call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, d, c, b, * opvert, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) elseif(c.le.8)then call ipsig2(x, y, z, x2, y2, z2, xc, yc, zc, opvert, d, * opvert, c, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) else call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, d, opvert, * c, b, a, c, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) endif if(itide .ge. 0) go to 1000 go to 990 elseif ((sidist(2).eq.0) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.-1)) then if(c.gt.8) go to 1000 if(opvert.gt.8.and.d.gt.8) go to 990 if(opvert.gt.8)then call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, d, c, b, * opvert, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) elseif(d.le.8)then call ipsig2(x, y, z, x2, y2, z2, xc, yc, zc, opvert, d, * opvert, c, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) else call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, c, opvert, * b, d, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) endif if(itide .ge. 0) go to 1000 go to 990 elseif ((sidist(2).eq.1) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.0)) then if(c.le.8.and.d.le.8) stop 1670 if(opvert.gt.8) go to 990 if(c.le.8)then call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, c, opvert, * b, d, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) else call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, d, opvert, * c, b, a, c, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) endif if(itide .ge. 0) go to 1000 go to 990 elseif ((sidist(2).eq.0) .and. (sidist(3).eq.1) .and. * (sidist(4).eq.0)) then if(opvert.gt.8.and.d.gt.8)then call iqsig2(x, y, z, w, x2, y2, z2, w2, b, d, opvert, a, * mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide.ge.0) go to 1000 go to 990 endif if(opvert.gt.8) go to 990 if(d.gt.8) go to 1000 if(c.gt.8)then call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, d, opvert, * c, b, a, c, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) else call ipsig2(x, y, z, x2, y2, z2, xc, yc, zc, opvert, d, * opvert, c, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) endif if(itide .ge. 0) go to 1000 go to 990 elseif ((sidist(2).eq.0) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.1)) then if(opvert.gt.8.and.c.gt.8)then call iqsig2(x, y, z, w, x2, y2, z2, w2, b, opvert, c, a, * mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide.ge.0) go to 1000 go to 990 endif if(opvert.gt.8) go to 990 if(c.gt.8) go to 1000 if(d.gt.8)then call ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, c, opvert, * b, d, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) else call ipsig2(x, y, z, x2, y2, z2, xc, yc, zc, opvert, d, * opvert, c, a, b, mhalf, mfull, isclp, istt) if(istt.gt.0) go to 990 if(istt.lt.0) go to 1000 call iqsign(x, y, z, w, x2, y2, z2, w2, opvert, d, c, b, * a, mhalf, mfull, isclp, isclw, isclr, delaun, itide) endif if(itide .ge. 0) go to 1000 go to 990 elseif(delaun) then go to 1000 endif c if ((sidist(2).eq.0) .and. (sidist(3).eq.-1) .and. * (sidist(4).eq.-1)) then oddsid = 2 go to 900 elseif ((sidist(2).eq.-1) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.-1)) then oddsid = 3 go to 900 elseif ((sidist(2).eq.-1) .and. (sidist(3).eq.-1) .and. * (sidist(4).eq.0)) then oddsid = 4 go to 900 elseif ((sidist(2).eq.0) .and. (sidist(3).eq.1) .and. * (sidist(4).eq.1)) then oddsid = 2 if(opvert.le.8)go to 900 isodd = ikon(6,1) if(isodd.le.8) stop 1680 call iqsig1(x, y, z, w, x2, y2, z2, w2, isodd, opvert, k, * mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 elseif ((sidist(2).eq.1) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.1)) then oddsid = 3 if(opvert.le.8)go to 900 isodd = ikon(7,1) if(isodd.le.8) stop 1690 call iqsig1(x, y, z, w, x2, y2, z2, w2, isodd, opvert, k, * mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 elseif ((sidist(2).eq.1) .and. (sidist(3).eq.1) .and. * (sidist(4).eq.0)) then oddsid = 4 if(opvert.le.8)go to 900 isodd = ikon(8,1) if(isodd.le.8) stop 1710 call iqsig1(x, y, z, w, x2, y2, z2, w2, isodd, opvert, k, * mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 elseif ((sidist(2).eq.0) .and. (sidist(3).eq.-1) .and. * (sidist(4).eq.1)) then oddsid = 2 isite = ikon(7,1) if(opvert.le.8 .or. isite.le.8)go to 900 isodd = ikon(6,1) if(isodd.le.8) stop 1720 call iqsig2(x, y, z, w, x2, y2, z2, w2, isodd, isite, opvert, * k, mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 elseif ((sidist(2).eq.1) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.-1)) then oddsid = 3 isite = ikon(8,1) if(opvert.le.8 .or. isite.le.8)go to 900 stop 1730 elseif ((sidist(2).eq.-1) .and. (sidist(3).eq.1) .and. * (sidist(4).eq.0)) then oddsid = 4 isite = ikon(6,1) if(opvert.le.8 .or. isite.le.8)go to 900 isodd = ikon(8,1) if(isodd.le.8) stop 1740 call iqsig2(x, y, z, w, x2, y2, z2, w2, isodd, isite, opvert, * k, mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 elseif ((sidist(2).eq.0) .and. (sidist(3).eq.1) .and. * (sidist(4).eq.-1)) then oddsid = -2 isite = ikon(8,1) if(opvert.le.8 .or. isite.le.8)go to 800 isodd = ikon(6,1) if(isodd.le.8) stop 1750 call iqsig2(x, y, z, w, x2, y2, z2, w2, isodd, isite, opvert, * k, mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 elseif ((sidist(2).eq.-1) .and. (sidist(3).eq.0) .and. * (sidist(4).eq.1)) then oddsid = -3 isite = ikon(6,1) if(opvert.le.8 .or. isite.le.8)go to 800 isodd = ikon(7,1) if(isodd.le.8) stop 1760 call iqsig2(x, y, z, w, x2, y2, z2, w2, isodd, isite, opvert, * k, mhalf, mfull, isclp, isclw, isclr, delaun, itide) if(itide .ge. 0) go to 1000 elseif ((sidist(2).eq.1) .and. (sidist(3).eq.-1) .and. *(sidist(4).eq.0)) then oddsid = -4 isite = ikon(7,1) if(opvert.le.8 .or. isite.le.8) go to 800 stop 1770 endif go to 990 c 800 continue oddsid=-oddsid 900 continue isite = ikon(oddsid+4,1) if(isite.le.8) stop 1780 go to 1000 c c obtain next tetrahedron c 990 continue idl(iscur) = iftal if(opvert.gt.8) then dcur = (xi(opvert) - xi(k))**2 + (yi(opvert) - yi(k))**2 + * (zi(opvert) - zi(k))**2 - wi(opvert) if(dsign(dsqrt(dabs(dcur)),dcur) .ge. dsign(dsqrt(dabs(dmin)), * dmin)+epz) go to 995 if(dsign(dsqrt(dabs(dcur)),dcur) .le. dsign(dsqrt(dabs(dmin)), * dmin)-epz) then lmin = opvert dmin = dcur go to 995 endif call ddsign(x, y, z, w, x2, y2, z2, w2, k, lmin, opvert, * mhalf, mfull, isclp, isclw, isclr, delaun, ipout) if(ipout.eq.-1) then lmin = opvert dmin = dcur endif endif 995 continue islst = iscur indx = 1 iscur = icon(1,islst) if(iscur.eq.0) go to 1100 b = icon(6,islst) c = icon(7,islst) d = icon(8,islst) if(idl(iscur) .ne. iftal) go to 200 go to 1100 c c obtain next tetrahedron c 1000 continue if(opvert.gt.8) then dcur = (xi(opvert) - xi(k))**2 + (yi(opvert) - yi(k))**2 + * (zi(opvert) - zi(k))**2 - wi(opvert) if(dsign(dsqrt(dabs(dcur)),dcur) .ge. dsign(dsqrt(dabs(dmun)), * dmun)+epz) go to 1100 if(dsign(dsqrt(dabs(dcur)),dcur) .le. dsign(dsqrt(dabs(dmun)), * dmun)-epz) then lmun = opvert dmun = dcur go to 1100 endif call ddsign(x, y, z, w, x2, y2, z2, w2, k, lmun, opvert, * mhalf, mfull, isclp, isclw, isclr, delaun, ipout) if(ipout.eq.-1) then lmun = opvert dmun = dcur endif endif 1100 continue if(indx.eq.1) then indx = 2 iscur = icon(2,islst) if(iscur.eq.0) go to 1100 b = icon(5,islst) c = icon(8,islst) d = icon(7,islst) if(idl(iscur) .ne. iftal) go to 200 go to 1100 elseif(indx.eq.2) then indx = 3 iscur = icon(3,islst) if(iscur.eq.0) go to 1100 b = icon(5,islst) c = icon(6,islst) d = icon(8,islst) if(idl(iscur) .ne. iftal) go to 200 go to 1100 elseif(indx.eq.3) then if(islst .ne. isini) then iscur = islst islst = icon(4,iscur) if(islst .le. 0) stop 1810 if(idl(islst) .ne. iftal) stop 1815 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 1820 endif go to 1100 else indx = 4 iscur = icon(4,islst) if(iscur.eq.0) go to 1100 b = icon(5,islst) c = icon(7,islst) d = icon(6,islst) if(idl(iscur) .ne. iftal) go to 200 go to 1100 endif endif if(islst .ne. isini) stop 1830 c if(dsign(dsqrt(dabs(dmun)),dmun) .ge. dsign(dsqrt(dabs(dmin)), * dmin)+epz) go to 1300 if(dsign(dsqrt(dabs(dmun)),dmun) .le. dsign(dsqrt(dabs(dmin)), * dmin)-epz) go to 1200 call ddsign(x, y, z, w, x2, y2, z2, w2, k, lmin, lmun, * mhalf, mfull, isclp, isclw, isclr, delaun, ipout) if(ipout.ne.-1) go to 1300 1200 continue write(*,*)'Polyhedron found does not appear to contain point.' write(*,*)'Program terminated.' WRITE(*,*)'DMUN=',DMUN,' DMIN=',DMIN stop 1850 1300 continue icl(ka) = lmin-8 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 *VRTORD c c This routine will order vertices a, b, c, d of a tetrahedron c so that a>b, b>c, b>d. Data structure is updated c subroutine vrtord(icon, curr, a, b, c, d) c integer icon(8,*), a, b, c, d, curr, it c if(a.lt.b)then it=a a=b b=it endif if(a.lt.c)then it=a a=c c=it endif if(a.lt.d)then it=a a=d d=it endif if(b.lt.c) b=c if(b.lt.d) b=d call reordr(icon,a,b,curr) c = icon(7,curr) d = icon(8,curr) if(b.gt.a.or.c.gt.b.or.d.gt.b) stop 3110 c 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 *TSTNCE c c This subroutine will compute the shortest distance from a point c to a facet of a tetrahedron. c subroutine tstnce(x, y, z, p, q, r, epz, k, dist, ipossi, idype) c integer p, q, r, k, ipossi, idype double precision x(*), y(*), z(*) double precision epz, dist, dust, dest, dust1, dust2, dust3 double precision xvec1, yvec1, zvec1, xvec2, yvec2, zvec2 double precision xvec3, yvec3, zvec3, dst1, dst2, dst3 double precision dotx, doty, dotz, dutx, duty, dutz double precision dmax, dmux, dlan, dlen, dlun, dotw double precision xvecp, yvecp, zvecp, dstp, dstq, dstr double precision xvecq, yvecq, zvecq, xvecr, yvecr, zvecr 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 dlan = dsqrt (dotx**2 + doty**2 + dotz**2) if(dlan.lt.epz .or. dlan/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(dlan,dlun) dotw=xvecp*dotx+yvecp*doty+zvecp*dotz dist=dotw/dlun if(dist.gt.-epz .and. dist.lt.epz)then ipossi = 1 go to 1000 endif dist=dotw/dlan if(dist.gt.-epz .and. dist.lt.epz)then ipossi = 1 go to 1000 endif if(dist.le.-epz) go to 1000 c dutx = yvec1 * dotz - doty * zvec1 duty = - xvec1 * dotz + dotx * zvec1 dutz = xvec1 * doty - dotx * yvec1 dmux = dst1*dmax dmux = dmax1(dlan,dmux) dlen = dsqrt (dutx**2 + duty**2 + dutz**2) if(dlen.lt.epz .or. dlen/dmux.lt.epz)then ipossi = 1 go to 1000 endif c dlun=dstp*dmux dlun=dmax1(dlen,dlun) dotw=xvecp*dutx+yvecp*duty+zvecp*dutz dust=dotw/dlun if(dust.gt.-epz .and. dust.lt.epz)then ipossi = 1 go to 1000 endif dust=dotw/dlen if(dust.gt.-epz .and. dust.lt.epz)then ipossi = 1 go to 1000 endif dust1 = dust if(dust1.le.-epz) go to 100 c dlun = dmax1(dstp,dst1) dest = (xvecp*xvec1 + yvecp*yvec1 + zvecp*zvec1)/dlun if(dest.gt.-epz .and. dest.lt.epz)then ipossi = 1 go to 1000 endif if(dest.le.-epz) go to 100 c xvecq = x(k) - x(q) yvecq = y(k) - y(q) zvecq = z(k) - z(q) dstq=dsqrt(xvecq**2+yvecq**2+zvecq**2) if(dstq.lt.epz) then ipossi = 1 go to 1000 endif c dlun = dmax1(dstq,dst1) dest = (xvecq*xvec1 + yvecq*yvec1 + zvecq*zvec1)/dlun if(dest.gt.-epz .and. dest.lt.epz)then ipossi = 1 go to 1000 endif if(dest.le.-epz) go to 900 dust2 = -1.0d0 go to 200 c 100 continue dutx = doty * zvec2 - yvec2 * dotz duty = - dotx * zvec2 + xvec2 * dotz dutz = dotx * yvec2 - xvec2 * doty dmux = dst2*dmax dmux = dmax1(dlan,dmux) dlen = dsqrt (dutx**2 + duty**2 + dutz**2) if(dlen.lt.epz .or. dlen/dmux.lt.epz)then ipossi = 1 go to 1000 endif c dlun=dstp*dmux dlun=dmax1(dlen,dlun) dotw=xvecp*dutx+yvecp*duty+zvecp*dutz dust=dotw/dlun if(dust.gt.-epz .and. dust.lt.epz)then ipossi = 1 go to 1000 endif dust=dotw/dlen if(dust.gt.-epz .and. dust.lt.epz)then ipossi = 1 go to 1000 endif dust2 = dust if(dust1.le.-epz .and. dust2.le.-epz) go to 200 if(dust2.ge.epz) go to 125 dist = dstp idype = p go to 1000 c 125 continue dlun = dmax1(dstp,dst2) dest = (xvecp*xvec2 + yvecp*yvec2 + zvecp*zvec2)/dlun if(dest.gt.-epz .and. dest.lt.epz)then ipossi = 1 go to 1000 endif if(dest.ge.epz) go to 150 dist = dstp idype = p go to 1000 c 150 continue xvecr = x(k) - x(r) yvecr = y(k) - y(r) zvecr = z(k) - z(r) dstr=dsqrt(xvecr**2+yvecr**2+zvecr**2) if(dstr.lt.epz) then ipossi = 1 go to 1000 endif c dlun = dmax1(dstr,dst2) dest = (xvecr*xvec2 + yvecr*yvec2 + zvecr*zvec2)/dlun if(dest.gt.-epz .and. dest.lt.epz)then ipossi = 1 go to 1000 endif if(dest.le.-epz) go to 900 c 200 continue dutx = doty * zvec3 - yvec3 * dotz duty = - dotx * zvec3 + xvec3 * dotz dutz = dotx * yvec3 - xvec3 * doty dmux = dst3*dmax dmux = dmax1(dlan,dmux) dlen = dsqrt (dutx**2 + duty**2 + dutz**2) if(dlen.lt.epz .or. dlen/dmux.lt.epz)then ipossi = 1 go to 1000 endif c xvecr = x(k) - x(r) yvecr = y(k) - y(r) zvecr = z(k) - z(r) dstr=dsqrt(xvecr**2+yvecr**2+zvecr**2) if(dstr.lt.epz) then ipossi = 1 go to 1000 endif c dlun=dstr*dmux dlun=dmax1(dlen,dlun) dotw=xvecr*dutx+yvecr*duty+zvecr*dutz dust=dotw/dlun if(dust.gt.-epz .and. dust.lt.epz)then ipossi = 1 go to 1000 endif dust=dotw/dlen if(dust.gt.-epz .and. dust.lt.epz)then ipossi = 1 go to 1000 endif dust3 = dust if(dust1.le.-epz .and. dust2.le.-epz .and. * dust3.le.-epz) go to 1000 if(dust2.ge.epz .or. dust3.ge.epz) go to 220 dist = dstq idype = q go to 1000 c 220 continue if(dust3.ge.epz) go to 225 dist = dstr idype = r go to 1000 c 225 continue dlun = dmax1(dstr,dst3) dest = (xvecr*xvec3 + yvecr*yvec3 + zvecr*zvec3)/dlun if(dest.gt.-epz .and. dest.lt.epz)then ipossi = 1 go to 1000 endif if(dest.ge.epz) go to 250 dist = dstr idype = r go to 1000 c 250 continue xvecq = x(k) - x(q) yvecq = y(k) - y(q) zvecq = z(k) - z(q) dstq=dsqrt(xvecq**2+yvecq**2+zvecq**2) if(dstq.lt.epz) then ipossi = 1 go to 1000 endif c dlun = dmax1(dstq,dst3) dest = (xvecq*xvec3 + yvecq*yvec3 + zvecq*zvec3)/dlun if(dest.gt.-epz .and. dest.lt.epz)then ipossi = 1 go to 1000 endif if(dest.le.-epz) go to 900 dist = dstq idype = q go to 1000 c 900 continue dist = dsqrt(dist**2 + dust**2) c 1000 continue return end *DSTNCE c c This subroutine will compute the perpendicular distance from a c point to a facet of a tetrahedron with emphasis on the sign. 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 *DSTANC c c This subroutine will compute the perpendicular distance from a point c to a facet of a tetrahedron with emphasis on both the sign and c magnitude. c subroutine dstanc(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, dotw 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) dotw=xvecp*dotx+yvecp*doty+zvecp*dotz dist=dotw/dlun if(dist.gt.-epz .and. dist.lt.epz)then ipossi = 1 go to 1000 endif dist=dotw/dlen if(dist.gt.-epz .and. dist.lt.epz)then ipossi = 1 endif c 1000 continue return end *CTRAD c c This subroutine will compute the orthogonal center of a tetrahedron c subroutine ctrad(x, y, z, w, xctr, yctr, zctr, a, b, c, d, * epz, delaun, ipossi) c double precision x(*), y(*), z(*), w(*) double precision epz, xctr, yctr, zctr double precision xm, ym, zm, xn, yn, zn, xu, yu, zu, xv, yv, zv double precision xw, yw, zw, xq, yq, zq, xe, ye, ze, xl, yl, zl double precision xt, yt, zt double precision norm, lambda, normu, normv, denom, dmax integer a, b, c, d, ipossi logical delaun c c initialize c ipossi = 0 c c find midpoints of edges ac and ab c xm = (x(a) + x(c)) / 2.0d0 ym = (y(a) + y(c)) / 2.0d0 zm = (z(a) + z(c)) / 2.0d0 c xn = (x(a) + x(b)) / 2.0d0 yn = (y(a) + y(b)) / 2.0d0 zn = (z(a) + z(b)) / 2.0d0 c c compute edge vectors u and v for edges ac and ab c xu = x(c) - x(a) yu = y(c) - y(a) zu = z(c) - z(a) c xv = x(b) - x(a) yv = y(b) - y(a) zv = z(b) - z(a) c c compute lengths of u and v c normu = dsqrt(xu**2 + yu**2 + zu**2) normv = dsqrt(xv**2 + yv**2 + zv**2) if(normu.lt.epz .or. normv.lt.epz) then ipossi = 1 go to 1000 endif dmax = dmax1(normu,normv) c c find perpendicular to facet abc of tetrahedron c xw = yu * zv - zu * yv yw = -xu * zv + zu * xv zw = xu * yv - yu * xv c c test whether edges ac, ab are colinear c norm = dsqrt(xw**2 + yw**2 + zw**2)/dmax if(norm .lt. epz)then ipossi = 1 go to 1000 endif xw = xw/normu yw = yw/normu zw = zw/normu c c normalize u and v c xu = xu / normu yu = yu / normu zu = zu / normu xv = xv / normv yv = yv / normv zv = zv / normv c c compute orthogonal center of edge ac c if(.not.delaun)then lambda = ((w(a)-w(c))/normu)/2.0d0 xm = xm + lambda*xu ym = ym + lambda*yu zm = zm + lambda*zu endif c c compute orthogonal center of edge ab c if(.not.delaun)then lambda = ((w(a)-w(b))/normv)/2.0d0 xn = xn + lambda*xv yn = yn + lambda*yv zn = zn + lambda*zv endif c c find perpendicular to edge v in plane that contains facet abc c xq = yw * zv - zw * yv yq = -xw * zv + zw * xv zq = xw * yv - yw * xv norm = dsqrt(xq**2 + yq**2 + zq**2) if(norm.lt.epz) then ipossi = 1 go to 1000 endif c c compute orthogonal center of facet abc c denom = xu*xq + yu*yq + zu*zq if(denom .gt. -epz .and. denom .lt. epz) then ipossi = 1 go to 1000 endif lambda = (xu*(xm-xn) + yu*(ym-yn) + zu*(zm-zn)) / denom c xe = xn + lambda*xq ye = yn + lambda*yq ze = zn + lambda*zq c c compute edge vector t for edge ad c xl = (x(a) + x(d)) / 2.0d0 yl = (y(a) + y(d)) / 2.0d0 zl = (z(a) + z(d)) / 2.0d0 c xt = x(d) - x(a) yt = y(d) - y(a) zt = z(d) - z(a) norm = dsqrt(xt**2 + yt**2 + zt**2) if(norm .lt. epz) then ipossi = 1 go to 1000 endif xt = xt / norm yt = yt / norm zt = zt / norm c c compute orthogonal center of edge ad c if(.not.delaun)then lambda = ((w(a)-w(d))/norm)/2.0d0 xl = xl + lambda*xt yl = yl + lambda*yt zl = zl + lambda*zt endif c c compute orthogonal center of tetrahedron c denom = xt*xw + yt*yw + zt*zw if(denom .gt. -epz .and. denom .lt. epz) then ipossi = 1 go to 1000 endif lambda = (xt*(xl-xe) + yt*(yl-ye) + zt*(zl-ze)) / denom c xctr = xe + lambda*xw yctr = ye + lambda*yw zctr = ze + lambda*zw c 1000 continue return end *BISPHR c c This subroutine will compute the distance from a point c (xctr,yctr,zctr) to the chordale plane between two points c opvert and ivrt c subroutine bisphr(x, y, z, w, opvert, ivrt, epz, * xctr, yctr, zctr, tdist, delaun, ipossi) c double precision x(*), y(*), z(*), w(*), norm integer opvert, ivrt, ipossi double precision epz, tdist, wambda, xctr, yctr, zctr, dif double precision xm, ym, zm, xu, yu, zu, xd, yd, zd, dmax double precision xu2, yu2, zu2 logical delaun c c find midpoint of edge from opvert to ivrt c xm = (x(opvert) + x(ivrt)) / 2.0d0 ym = (y(opvert) + y(ivrt)) / 2.0d0 zm = (z(opvert) + z(ivrt)) / 2.0d0 c c find vector from ivrt to opvert c xu = x(opvert) - x(ivrt) yu = y(opvert) - y(ivrt) zu = z(opvert) - z(ivrt) c norm = dsqrt(xu**2 + yu**2 + zu**2) if(norm .lt. epz) then ipossi = 1 go to 1000 endif xu2 = xu/norm yu2 = yu/norm zu2 = zu/norm c c compute orthogonal center of edge ivrt-opvert c if(.not.delaun)then wambda = ((w(ivrt)-w(opvert))/norm)/2.0d0 xm = xm + wambda*xu2 ym = ym + wambda*yu2 zm = zm + wambda*zu2 endif c c compute distance c xd = xctr - xm yd = yctr - ym zd = zctr - zm dif = dsqrt(xd**2 + yd**2 + zd**2) dmax = dmax1(norm,dif) tdist = (xd*xu + yd*yu + zd*zu) / dmax if(tdist.gt. -epz .and. tdist.lt. epz) then ipossi = 1 endif c 1000 continue return end *DSTEXA c c Routine to compute almost exact perpendicular distance from c a plane spanned by points ifir, isec, ithi to a point ifou c subroutine dstexa(x, y, z, x2, y2, z2, ifir, isec, ithi, ifou, * mhalf, mfull, isclp, r215, deps, dscle, dist, * ipout) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer nkmax parameter (nkmax = 30) integer io(nkmax), iox(nkmax), ioy(nkmax), ioz(nkmax) integer ifir, isec, ithi, ifou, mhalf, mfull, isclp(*) integer isgo, iko, isgox, ikox, isgoy, ikoy, isgoz, ikoz, ipout double precision dist, dnux double precision r215, dnom, xnum, ynum, znum, dnum, deps, dscle c call crsinn(x, y, z, x2, y2, z2, ifir, isec, ithi, ifou, mhalf, * mfull, isclp, io, isgo, iko, iox, isgox, ikox, * ioy, isgoy, ikoy, ioz, isgoz, ikoz) ipout = isgo call doubnm(io, isgo, iko, r215, dnum) call doubnm(iox, isgox, ikox, r215, xnum) call doubnm(ioy, isgoy, ikoy, r215, ynum) call doubnm(ioz, isgoz, ikoz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3310 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) if(dnom.lt.deps) stop 3320 dist = ((dnum/dnux)/dnom)/dscle c return end *TSTEXA c c Routine to compute either almost exact perpendicular distance c or almost exact shortest distance from a plane spanned c by points ifir, isec, ithi to a point ifou c subroutine tstexa(x, y, z, x2, y2, z2, ifir, isec, ithi, ifou, * mhalf, mfull, isclp, r215, deps, dscle, dist, * ipout, idype) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer nkmax parameter (nkmax = 30) integer io(nkmax), iox(nkmax), ioy(nkmax), ioz(nkmax) integer ia(nkmax), iax(nkmax), iay(nkmax), iaz(nkmax) integer ifir, isec, ithi, ifou, mhalf, mfull, isclp(*) integer isgo, iko, isgox, ikox, isgoy, ikoy, isgoz, ikoz integer isga, ika, isgax, ikax, isgay, ikay, isgaz, ikaz integer ipout, idype, idst1, idst2, idst3, isig1, isig2 double precision dist, dust, dnux double precision r215, dnom, xnum, ynum, znum, dnum, deps, dscle c call crsinn(x, y, z, x2, y2, z2, ifir, isec, ithi, ifou, mhalf, * mfull, isclp, io, isgo, iko, iox, isgox, ikox, * ioy, isgoy, ikoy, ioz, isgoz, ikoz) ipout = isgo if(ipout.le.0) go to 1000 call doubnm(io, isgo, iko, r215, dnum) call doubnm(iox, isgox, ikox, r215, xnum) call doubnm(ioy, isgoy, ikoy, r215, ynum) call doubnm(ioz, isgoz, ikoz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3410 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) if(dnom.lt.deps) stop 3420 dist = ((dnum/dnux)/dnom)/dscle c call crsunn(x, y, z, x2, y2, z2, ifir, isec, ifou, mhalf, mfull, * isclp, iox, isgox, ikox, ioy, isgoy, ikoy, * ioz, isgoz, ikoz, ia, isga, ika, iax, isgax, ikax, * iay, isgay, ikay, iaz, isgaz, ikaz) idst1 = isga if(idst1.le.0) go to 100 c call innunn(x, y, z, x2, y2, z2, ifir, isec, ifou, mhalf, mfull, * isclp, isig1, isig2) if(isig1.lt.0) go to 100 c if(isig2.gt.0) go to 50 call doubnm(ia, isga, ika, r215, dnum) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3430 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) if(dnom.lt.deps) stop 3440 dust = ((dnum/dnux)/dnom)/dscle go to 900 c 50 continue idst2 = -1 go to 200 c 100 continue call crsunn(x, y, z, x2, y2, z2, ithi, ifir, ifou, mhalf, mfull, * isclp, iox, isgox, ikox, ioy, isgoy, ikoy, * ioz, isgoz, ikoz, ia, isga, ika, iax, isgax, ikax, * iay, isgay, ikay, iaz, isgaz, ikaz) idst2 = isga if(idst1.le.0 .and. idst2.le.0) go to 200 if(idst2.gt.0) go to 125 call vectrd(x, y, z, x2, y2, z2, ifou, ifir, mhalf, mfull, isclp, * iax, isgax, ikax, iay, isgay, ikay, iaz, isgaz, ikaz) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3450 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) dist = (dnux/dscle)*dnom idype = ifir go to 1000 c 125 continue call innunn(x, y, z, x2, y2, z2, ifir, ithi, ifou, mhalf, mfull, * isclp, isig1, isig2) if(isig1.ge.0) go to 150 call vectrd(x, y, z, x2, y2, z2, ifou, ifir, mhalf, mfull, isclp, * iax, isgax, ikax, iay, isgay, ikay, iaz, isgaz, ikaz) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3460 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) dist = (dnux/dscle)*dnom idype = ifir go to 1000 c 150 continue if(isig2.gt.0) go to 200 call doubnm(ia, isga, ika, r215, dnum) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3470 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) if(dnom.lt.deps) stop 3480 dust = ((dnum/dnux)/dnom)/dscle go to 900 c 200 continue call crsunn(x, y, z, x2, y2, z2, isec, ithi, ifou, mhalf, mfull, * isclp, iox, isgox, ikox, ioy, isgoy, ikoy, * ioz, isgoz, ikoz, ia, isga, ika, iax, isgax, ikax, * iay, isgay, ikay, iaz, isgaz, ikaz) idst3 = isga if(idst1.le.0 .and. idst2.le.0 .and. idst3.le.0) go to 1000 if(idst2.gt.0 .or. idst3.gt.0) go to 220 call vectrd(x, y, z, x2, y2, z2, ifou, isec, mhalf, mfull, isclp, * iax, isgax, ikax, iay, isgay, ikay, iaz, isgaz, ikaz) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3490 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) dist = (dnux/dscle)*dnom idype = isec go to 1000 c 220 continue if(idst3.gt.0) go to 225 call vectrd(x, y, z, x2, y2, z2, ifou, ithi, mhalf, mfull, isclp, * iax, isgax, ikax, iay, isgay, ikay, iaz, isgaz, ikaz) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3500 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) dist = (dnux/dscle)*dnom idype = ithi go to 1000 c 225 continue call innunn(x, y, z, x2, y2, z2, ithi, isec, ifou, mhalf, mfull, * isclp, isig1, isig2) if(isig1.ge.0) go to 250 call vectrd(x, y, z, x2, y2, z2, ifou, ithi, mhalf, mfull, isclp, * iax, isgax, ikax, iay, isgay, ikay, iaz, isgaz, ikaz) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3510 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) dist = (dnux/dscle)*dnom idype = ithi go to 1000 c 250 continue if(isig2.gt.0) go to 275 call doubnm(ia, isga, ika, r215, dnum) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3520 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) if(dnom.lt.deps) stop 3530 dust = ((dnum/dnux)/dnom)/dscle go to 900 c 275 continue call vectrd(x, y, z, x2, y2, z2, ifou, isec, mhalf, mfull, isclp, * iax, isgax, ikax, iay, isgay, ikay, iaz, isgaz, ikaz) call doubnm(iax, isgax, ikax, r215, xnum) call doubnm(iay, isgay, ikay, r215, ynum) call doubnm(iaz, isgaz, ikaz, r215, znum) dnux = dmax1(dabs(xnum),dabs(ynum),dabs(znum)) if(dnux.lt.deps) stop 3540 xnum = xnum/dnux ynum = ynum/dnux znum = znum/dnux dnom = dsqrt(xnum**2+ynum**2+znum**2) dist = (dnux/dscle)*dnom idype = isec go to 1000 c 900 continue dist = dsqrt(dist**2 + dust**2) c 1000 continue c return end *DDSIGN c c subroutine for determining which is larger between c the power distance of point ifir from point isec and the c power distance of point ifir from point ithi using exact c arithmetic: if computation results in a positive c number then distance of ifir from ithi is larger, etc. c subroutine ddsign(x, y, z, w, x2, y2, z2, w2, ifir, isec, ithi, * mhalf, mfull, isclp, isclw, isclr, * delaun, ipout) c integer x(*), y(*), z(*), w(*), x2(*), y2(*), z2(*), w2(*) integer ifir, isec, ithi integer isclp(*), isclw(*), isclr(*), mhalf, mfull, nkmax, ipout logical delaun parameter (nkmax = 30) integer io(nkmax), iu(nkmax) integer ix2(nkmax), iy2(nkmax), iz2(nkmax) integer ix3(nkmax), iy3(nkmax), iz3(nkmax) integer ixf(nkmax), iyf(nkmax), izf(nkmax) integer iw2(nkmax), iw3(nkmax) integer ix22(nkmax), iy22(nkmax), iz22(nkmax) integer ix32(nkmax), iy32(nkmax), iz32(nkmax) integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew integer ixthw, iythw, izthw integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2 integer ixth2, iyth2, izth2 integer iwsew, iwthw, iwse2, iwth2 integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2 integer isgx3, isgy3, isgz3, ikx3, iky3, ikz3 integer isgx22, isgy22, isgz22, ikx22, iky22, ikz22 integer isgx32, isgy32, isgz32, ikx32, iky32, ikz32 integer isgw2, isgw3, ikw2, ikw3, isgcl, ikcl integer isgo, isgu, iko, iku c if(delaun) then isgw2 = 0 isgw3 = 0 else iwsew = w(isec) iwthw = w(ithi) c iwse2 = w2(isec) iwth2 = w2(ithi) c isgcl = 1 ikcl = 2 call decmp2(io, isgo, iko, iwsew, iwse2, mhalf, mfull, isclw) call mulmul(io, isclr, iw2, isgo, isgcl, isgw2, iko, ikcl, * ikw2, nkmax, mhalf) call decmp2(io, isgo, iko, iwthw, iwth2, mhalf, mfull, isclw) call mulmul(io, isclr, iw3, isgo, isgcl, isgw3, iko, ikcl, * ikw3, nkmax, mhalf) endif 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(ix2, ix2, ix22, isgx2, isgx2, isgx22, ikx2, ikx2, * ikx22, nkmax, mhalf) call mulmul(iy2, iy2, iy22, isgy2, isgy2, isgy22, iky2, iky2, * iky22, nkmax, mhalf) call mulmul(iz2, iz2, iz22, isgz2, isgz2, isgz22, ikz2, ikz2, * ikz22, nkmax, mhalf) if(isgx22.lt.0 .or. isgy22.lt.0 .or. isgz22.lt.0) stop 3610 c call mulmul(ix3, ix3, ix32, isgx3, isgx3, isgx32, ikx3, ikx3, * ikx32, nkmax, mhalf) call mulmul(iy3, iy3, iy32, isgy3, isgy3, isgy32, iky3, iky3, * iky32, nkmax, mhalf) call mulmul(iz3, iz3, iz32, isgz3, isgz3, isgz32, ikz3, ikz3, * ikz32, nkmax, mhalf) if(isgx32.lt.0 .or. isgy32.lt.0 .or. isgz32.lt.0) stop 3620 c isgy32 = -isgy32 call muldif(ix32, iy32, io, isgx32, isgy32, isgo, * ikx32, iky32, iko, nkmax, mhalf) isgz32 = -isgz32 call muldif(io, iz32, iu, isgo, isgz32, isgu, * iko, ikz32, iku, nkmax, mhalf) call muldif(iu, iw3, io, isgu, isgw3, isgo, * iku, ikw3, iko, nkmax, mhalf) call muldif(io, ix22, iu, isgo, isgx22, isgu, * iko, ikx22, iku, nkmax, mhalf) call muldif(iu, iy22, io, isgu, isgy22, isgo, * iku, iky22, iko, nkmax, mhalf) call muldif(io, iz22, iu, isgo, isgz22, isgu, * iko, ikz22, iku, nkmax, mhalf) isgw2 = -isgw2 call muldif(iu, iw2, io, isgu, isgw2, isgo, * iku, ikw2, iko, nkmax, mhalf) c ipout = isgo c return end *CRSINN c c Routine for determining cross product of two vectors c and , and inner product of this cross product with a c third vector c c subroutine crsinn(x, y, z, x2, y2, z2, ifir, isec, ithi, ifou, * mhalf, mfull, isclp, io, isgo, iko, iox, isgox, * ikox, ioy, isgoy, ikoy, ioz, isgoz, ikoz) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer io(*), iox(*),ioy(*), ioz(*) integer ifir, isec, ithi, ifou integer isclp(*), mhalf, mfull, nkmax parameter (nkmax = 30) integer 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, iko, isgox, ikox, isgoy, ikoy, isgoz, ikoz integer isgu, isgv, isgw, 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, iox, isgv, isgu, isgox, ikv, iku, ikox, * nkmax, mhalf) call mulmul(iox, ix4, io, isgox, isgx4, isgo, ikox, 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, ioy, isgv, isgu, isgoy, ikv, iku, ikoy, * nkmax, mhalf) call mulmul(ioy, iy4, iu, isgoy, isgy4, isgu, ikoy, 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, ioz, isgv, isgu, isgoz, ikv, iku, ikoz, * nkmax, mhalf) call mulmul(ioz, iz4, iu, isgoz, isgz4, isgu, ikoz, ikz4, iku, * nkmax, mhalf) isgu =-isgu call muldif(iw, iu, io, isgw, isgu, isgo, ikw, iku, iko, * nkmax, mhalf) c return end *CRSUNN c c Routine for determining cross product of two vectors c and (iox,ioy,ioz) and inner product of this cross product with a c third vector c c subroutine crsunn(x, y, z, x2, y2, z2, ifir, isec, ifou, * mhalf, mfull, isclp, iox, isgox, ikox, * ioy, isgoy, ikoy, ioz, isgoz, ikoz, * ia, isga, ika, iax, isgax, ikax, * iay, isgay, ikay, iaz, isgaz, ikaz) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer iox(*),ioy(*), ioz(*) integer ia(*), iax(*),iay(*), iaz(*) integer ifir, isec, ifou integer isclp(*), mhalf, mfull, nkmax parameter (nkmax = 30) integer io(nkmax), iu(nkmax), iv(nkmax), iw(nkmax) integer ix2(nkmax), iy2(nkmax), iz2(nkmax) integer ix4(nkmax), iy4(nkmax), iz4(nkmax) integer ixf(nkmax), iyf(nkmax), izf(nkmax) integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew integer ixfow, iyfow, izfow integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2 integer ixfo2, iyfo2, izfo2 integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2 integer isgx4, isgy4, isgz4, ikx4, iky4, ikz4 integer isgox, ikox, isgoy, ikoy, isgoz, ikoz integer isga, ika, isgax, ikax, isgay, ikay, isgaz, ikaz 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) 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) 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, 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, ioz, iv, isgy2, isgoz, isgv, iky2, ikoz, ikv, * nkmax, mhalf) call mulmul(iz2, ioy, iu, isgz2, isgoy, isgu, ikz2, ikoy, iku, * nkmax, mhalf) call muldif(iv, iu, iax, isgv, isgu, isgax, ikv, iku, ikax, * nkmax, mhalf) call mulmul(iax, ix4, ia, isgax, isgx4, isga, ikax, ikx4, ika, * nkmax, mhalf) c call mulmul(iz2, iox, iv, isgz2, isgox, isgv, ikz2, ikox, ikv, * nkmax, mhalf) call mulmul(ix2, ioz, iu, isgx2, isgoz, isgu, ikx2, ikoz, iku, * nkmax, mhalf) call muldif(iv, iu, iay, isgv, isgu, isgay, ikv, iku, ikay, * nkmax, mhalf) call mulmul(iay, iy4, iu, isgay, isgy4, isgu, ikay, iky4, iku, * nkmax, mhalf) isgu =-isgu call muldif(ia, iu, iw, isga, isgu, isgw, ika, iku, ikw, * nkmax, mhalf) c call mulmul(ix2, ioy, iv, isgx2, isgoy, isgv, ikx2, ikoy, ikv, * nkmax, mhalf) call mulmul(iy2, iox, iu, isgy2, isgox, isgu, iky2, ikox, iku, * nkmax, mhalf) call muldif(iv, iu, iaz, isgv, isgu, isgaz, ikv, iku, ikaz, * nkmax, mhalf) call mulmul(iaz, iz4, iu, isgaz, isgz4, isgu, ikaz, ikz4, iku, * nkmax, mhalf) isgu =-isgu call muldif(iw, iu, ia, isgw, isgu, isga, ikw, iku, ika, * nkmax, mhalf) c return end *INNUNN c c Routine for determining signs of inner products of vectors c and , and and c subroutine innunn(x, y, z, x2, y2, z2, ifir, isec, ifou, * mhalf, mfull, isclp, isga, isgo) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer ifir, isec, ifou integer isclp(*), mhalf, mfull, nkmax parameter (nkmax = 30) integer ia(nkmax), io(nkmax), iax(nkmax), iay(nkmax), iaz(nkmax) integer iu(nkmax), iw(nkmax) integer ix2(nkmax), iy2(nkmax), iz2(nkmax) integer ix4(nkmax), iy4(nkmax), iz4(nkmax) integer ixf(nkmax), iyf(nkmax), izf(nkmax) integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew integer ixfow, iyfow, izfow integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2 integer ixfo2, iyfo2, izfo2 integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2 integer isgx4, isgy4, isgz4, ikx4, iky4, ikz4 integer isga, ika, isgo, iko integer isgax, ikax, isgay, ikay, isgaz, ikaz integer isgu, isgw, iku, ikw c ixfiw = x(ifir) iyfiw = y(ifir) izfiw = z(ifir) ixsew = x(isec) iysew = y(isec) izsew = z(isec) 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) 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, iax, isgo, isgxf, isgax, iko, ikxf, ikax, * nkmax, mhalf) call decmp2(iu, isgu, iku, ixfow, ixfo2, mhalf, mfull, isclp) call muldif(iu, ixf, ix4, isgu, isgxf, isgx4, iku, ikxf, ikx4, * nkmax, mhalf) call muldif(iu, io, ix2, isgu, isgo, isgx2, iku, iko, ikx2, * nkmax, mhalf) call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp) call muldif(io, iyf, iay, isgo, isgyf, isgay, iko, ikyf, ikay, * nkmax, mhalf) call decmp2(iu, isgu, iku, iyfow, iyfo2, mhalf, mfull, isclp) call muldif(iu, iyf, iy4, isgu, isgyf, isgy4, iku, ikyf, iky4, * nkmax, mhalf) call muldif(iu, io, iy2, isgu, isgo, isgy2, iku, iko, iky2, * nkmax, mhalf) call decmp2(io, isgo, iko, izsew, izse2, mhalf, mfull, isclp) call muldif(io, izf, iaz, isgo, isgzf, isgaz, iko, ikzf, ikaz, * nkmax, mhalf) call decmp2(iu, isgu, iku, izfow, izfo2, mhalf, mfull, isclp) call muldif(iu, izf, iz4, isgu, isgzf, isgz4, iku, ikzf, ikz4, * nkmax, mhalf) call muldif(iu, io, iz2, isgu, isgo, isgz2, iku, iko, ikz2, * nkmax, mhalf) c call mulmul(iax, ix4, ia, isgax, isgx4, isga, ikax, ikx4, ika, * nkmax, mhalf) c call mulmul(iay, iy4, iu, isgay, isgy4, isgu, ikay, iky4, iku, * nkmax, mhalf) isgu =-isgu call muldif(ia, iu, iw, isga, isgu, isgw, ika, iku, ikw, * nkmax, mhalf) c call mulmul(iaz, iz4, iu, isgaz, isgz4, isgu, ikaz, ikz4, iku, * nkmax, mhalf) isgu =-isgu call muldif(iw, iu, ia, isgw, isgu, isga, ikw, iku, ika, * nkmax, mhalf) c call mulmul(iax, ix2, io, isgax, isgx2, isgo, ikax, ikx2, iko, * nkmax, mhalf) c call mulmul(iay, iy2, iu, isgay, isgy2, isgu, ikay, iky2, iku, * nkmax, mhalf) isgu =-isgu call muldif(io, iu, iw, isgo, isgu, isgw, iko, iku, ikw, * nkmax, mhalf) c call mulmul(iaz, iz2, iu, isgaz, isgz2, isgu, ikaz, ikz2, iku, * nkmax, mhalf) isgu =-isgu call muldif(iw, iu, io, isgw, isgu, isgo, ikw, iku, iko, * nkmax, mhalf) c return end *VECTRD c c Routine for determining vector c subroutine vectrd(x, y, z, x2, y2, z2, ifir, isec, mhalf, * mfull, isclp, iax, isgax, ikax, * iay, isgay, ikay, iaz, isgaz, ikaz) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*) integer iax(*),iay(*), iaz(*) integer ifir, isec integer isclp(*), mhalf, mfull, nkmax parameter (nkmax = 30) integer io(nkmax) integer ixf(nkmax), iyf(nkmax), izf(nkmax) integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2 integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf integer isgax, ikax, isgay, ikay, isgaz, ikaz integer isgo, iko c ixfiw = x(ifir) iyfiw = y(ifir) izfiw = z(ifir) ixsew = x(isec) iysew = y(isec) izsew = z(isec) c ixfi2 = x2(ifir) iyfi2 = y2(ifir) izfi2 = z2(ifir) ixse2 = x2(isec) iyse2 = y2(isec) izse2 = z2(isec) 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, iax, isgo, isgxf, isgax, iko, ikxf, ikax, * nkmax, mhalf) call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp) call muldif(io, iyf, iay, isgo, isgyf, isgay, iko, ikyf, ikay, * nkmax, mhalf) call decmp2(io, isgo, iko, izsew, izse2, mhalf, mfull, isclp) call muldif(io, izf, iaz, isgo, isgzf, isgaz, iko, ikzf, ikaz, * nkmax, mhalf) c return end *DOUBNM c subroutine doubnm(io, isgo, iko, r215, dnum) c integer io(*) double precision dnum, r215, rpwr integer isgo, iko, i c if(isgo.eq.0) then dnum = dble(0) go to 1000 else if(iko .lt. 2) stop 3710 if(iko .gt. 68) stop 3720 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 *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 *IPSIG1 c subroutine ipsig1(x, y, z, x2, y2, z2, xc, yc, zc, ifir, * isec, ithi, ifou, ifif, isix, mhalf, * mfull, isclp, ipout) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*) integer ifir, isec, ithi, ifou, ifif, isix 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, ixsiw, iysiw, izsiw integer ixfuw, iyfuw, izfuw, ixsuw, iysuw, izsuw integer ixth2, iyth2, izth2, ixfo2, iyfo2, izfo2 integer ixfi2, iyfi2, izfi2, ixsi2, iysi2, izsi2 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 ixthw = x(ithi) iythw = y(ithi) izthw = z(ithi) ixfow = x(ifou) iyfow = y(ifou) izfow = z(ifou) ixfiw = x(ifif) iyfiw = y(ifif) izfiw = z(ifif) ixsiw = x(isix) iysiw = y(isix) izsiw = z(isix) c ixth2 = x2(ithi) iyth2 = y2(ithi) izth2 = z2(ithi) ixfo2 = x2(ifou) iyfo2 = y2(ifou) izfo2 = z2(ifou) ixfi2 = x2(ifif) iyfi2 = y2(ifif) izfi2 = z2(ifif) ixsi2 = x2(isix) iysi2 = y2(isix) izsi2 = z2(isix) c ixfuw = xc(ifir) iyfuw = yc(ifir) izfuw = zc(ifir) ixsuw = xc(isec) iysuw = yc(isec) izsuw = zc(isec) c ikxf = 2 ikyf = 2 ikzf = 2 call decomp(ixf, isgxf, ixfuw, mhalf) call decomp(iyf, isgyf, iyfuw, mhalf) call decomp(izf, isgzf, izfuw, mhalf) c call decmp2(ix3, isgx3, ikx3, ixthw, ixth2, mhalf, mfull, isclp) call decmp2(iy3, isgy3, iky3, iythw, iyth2, mhalf, mfull, isclp) call decmp2(iz3, isgz3, ikz3, izthw, izth2, mhalf, mfull, isclp) call decmp2(ix5, isgx5, ikx5, ixfiw, ixfi2, mhalf, mfull, isclp) call decmp2(iy5, isgy5, iky5, iyfiw, iyfi2, mhalf, mfull, isclp) call decmp2(iz5, isgz5, ikz5, izfiw, izfi2, mhalf, mfull, isclp) c iko = 2 call decomp(io, isgo, ixsuw, mhalf) call muldif(io, ixf, ix2, isgo, isgxf, isgx2, iko, ikxf, ikx2, * nkmax, mhalf) call decomp(io, isgo, iysuw, mhalf) call muldif(io, iyf, iy2, isgo, isgyf, isgy2, iko, ikyf, iky2, * nkmax, mhalf) call decomp(io, isgo, izsuw, mhalf) call muldif(io, izf, iz2, isgo, isgzf, isgz2, iko, ikzf, ikz2, * nkmax, mhalf) c call decmp2(io, isgo, iko, ixfow, ixfo2, mhalf, mfull, isclp) call muldif(io, ix3, ix4, isgo, isgx3, isgx4, iko, ikx3, ikx4, * nkmax, mhalf) call decmp2(io, isgo, iko, iyfow, iyfo2, mhalf, mfull, isclp) call muldif(io, iy3, iy4, isgo, isgy3, isgy4, iko, iky3, iky4, * nkmax, mhalf) call decmp2(io, isgo, iko, izfow, izfo2, mhalf, mfull, isclp) call muldif(io, iz3, iz4, isgo, isgz3, isgz4, iko, ikz3, ikz4, * nkmax, mhalf) call decmp2(io, isgo, iko, ixsiw, ixsi2, mhalf, mfull, isclp) call muldif(io, ix5, ix6, isgo, isgx5, isgx6, iko, ikx5, ikx6, * nkmax, mhalf) call decmp2(io, isgo, iko, iysiw, iysi2, mhalf, mfull, isclp) call muldif(io, iy5, iy6, isgo, isgy5, isgy6, iko, iky5, iky6, * nkmax, mhalf) call decmp2(io, isgo, iko, izsiw, izsi2, mhalf, mfull, isclp) call muldif(io, iz5, iz6, isgo, isgz5, isgz6, iko, ikz5, ikz6, * nkmax, mhalf) c call mulmul(iy2, iz4, io, isgy2, isgz4, isgo, iky2, ikz4, iko, * nkmax, mhalf) call mulmul(io, ix6, iv, isgo, isgx6, isgv, iko, ikx6, ikv, * nkmax, mhalf) c call mulmul(iz2, ix4, io, isgz2, isgx4, isgo, ikz2, ikx4, iko, * nkmax, mhalf) call mulmul(io, iy6, iu, isgo, isgy6, isgu, iko, iky6, iku, * nkmax, mhalf) isgu =-isgu call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko, * nkmax, mhalf) c call mulmul(ix2, iy4, iv, isgx2, isgy4, isgv, ikx2, iky4, ikv, * nkmax, mhalf) call mulmul(iv, iz6, iu, isgv, isgz6, isgu, ikv, ikz6, iku, * nkmax, mhalf) isgu =-isgu call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv, * nkmax, mhalf) c call mulmul(iz2, iy4, io, isgz2, isgy4, isgo, ikz2, iky4, iko, * nkmax, mhalf) call mulmul(io, ix6, iu, isgo, isgx6, isgu, iko, ikx6, iku, * nkmax, mhalf) call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko, * nkmax, mhalf) c call mulmul(ix2, iz4, iv, isgx2, isgz4, isgv, ikx2, ikz4, ikv, * nkmax, mhalf) call mulmul(iv, iy6, iu, isgv, isgy6, isgu, ikv, iky6, iku, * nkmax, mhalf) call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv, * nkmax, mhalf) c call mulmul(iy2, ix4, io, isgy2, isgx4, isgo, iky2, ikx4, iko, * nkmax, mhalf) call mulmul(io, iz6, iu, isgo, isgz6, isgu, iko, ikz6, iku, * nkmax, mhalf) call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko, * nkmax, mhalf) c ipout = isgo c return end *IPSIG2 c subroutine ipsig2(x, y, z, x2, y2, z2, xc, yc, zc, ifir, * isec, ithi, ifou, ifif, isix, mhalf, * mfull, isclp, ipout) c integer x(*), y(*), z(*), x2(*), y2(*), z2(*), xc(*), yc(*), zc(*) integer ifir, isec, ithi, ifou, ifif, isix 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, ixsiw, iysiw, izsiw integer ixfuw, iyfuw, izfuw, ixsuw, iysuw, izsuw integer ixfi2, iyfi2, izfi2, ixsi2, iysi2, izsi2 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(ifif) iyfiw = y(ifif) izfiw = z(ifif) ixsiw = x(isix) iysiw = y(isix) izsiw = z(isix) c ixfi2 = x2(ifif) iyfi2 = y2(ifif) izfi2 = z2(ifif) ixsi2 = x2(isix) iysi2 = y2(isix) izsi2 = z2(isix) c ixfuw = xc(ifir) iyfuw = yc(ifir) izfuw = zc(ifir) ixsuw = xc(isec) iysuw = yc(isec) izsuw = zc(isec) ixthw = xc(ithi) iythw = yc(ithi) izthw = zc(ithi) ixfow = xc(ifou) iyfow = yc(ifou) izfow = zc(ifou) c ikxf = 2 ikyf = 2 ikzf = 2 ikx3 = 2 iky3 = 2 ikz3 = 2 call decomp(ixf, isgxf, ixfuw, mhalf) call decomp(iyf, isgyf, iyfuw, mhalf) call decomp(izf, isgzf, izfuw, mhalf) call decomp(ix3, isgx3, ixthw, mhalf) call decomp(iy3, isgy3, iythw, mhalf) call decomp(iz3, isgz3, izthw, mhalf) c call decmp2(ix5, isgx5, ikx5, ixfiw, ixfi2, mhalf, mfull, isclp) call decmp2(iy5, isgy5, iky5, iyfiw, iyfi2, mhalf, mfull, isclp) call decmp2(iz5, isgz5, ikz5, izfiw, izfi2, mhalf, mfull, isclp) c iko = 2 call decomp(io, isgo, ixsuw, mhalf) call muldif(io, ixf, ix2, isgo, isgxf, isgx2, iko, ikxf, ikx2, * nkmax, mhalf) call decomp(io, isgo, iysuw, mhalf) call muldif(io, iyf, iy2, isgo, isgyf, isgy2, iko, ikyf, iky2, * nkmax, mhalf) call decomp(io, isgo, izsuw, mhalf) call muldif(io, izf, iz2, isgo, isgzf, isgz2, iko, ikzf, ikz2, * nkmax, mhalf) call decomp(io, isgo, ixfow, mhalf) call muldif(io, ix3, ix4, isgo, isgx3, isgx4, iko, ikx3, ikx4, * nkmax, mhalf) call decomp(io, isgo, iyfow, mhalf) call muldif(io, iy3, iy4, isgo, isgy3, isgy4, iko, iky3, iky4, * nkmax, mhalf) call decomp(io, isgo, izfow, mhalf) call muldif(io, iz3, iz4, isgo, isgz3, isgz4, iko, ikz3, ikz4, * nkmax, mhalf) c call decmp2(io, isgo, iko, ixsiw, ixsi2, mhalf, mfull, isclp) call muldif(io, ix5, ix6, isgo, isgx5, isgx6, iko, ikx5, ikx6, * nkmax, mhalf) call decmp2(io, isgo, iko, iysiw, iysi2, mhalf, mfull, isclp) call muldif(io, iy5, iy6, isgo, isgy5, isgy6, iko, iky5, iky6, * nkmax, mhalf) call decmp2(io, isgo, iko, izsiw, izsi2, mhalf, mfull, isclp) call muldif(io, iz5, iz6, isgo, isgz5, isgz6, iko, ikz5, ikz6, * nkmax, mhalf) c call mulmul(iy2, iz4, io, isgy2, isgz4, isgo, iky2, ikz4, iko, * nkmax, mhalf) call mulmul(io, ix6, iv, isgo, isgx6, isgv, iko, ikx6, ikv, * nkmax, mhalf) c call mulmul(iz2, ix4, io, isgz2, isgx4, isgo, ikz2, ikx4, iko, * nkmax, mhalf) call mulmul(io, iy6, iu, isgo, isgy6, isgu, iko, iky6, iku, * nkmax, mhalf) isgu =-isgu call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko, * nkmax, mhalf) c call mulmul(ix2, iy4, iv, isgx2, isgy4, isgv, ikx2, iky4, ikv, * nkmax, mhalf) call mulmul(iv, iz6, iu, isgv, isgz6, isgu, ikv, ikz6, iku, * nkmax, mhalf) isgu =-isgu call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv, * nkmax, mhalf) c call mulmul(iz2, iy4, io, isgz2, isgy4, isgo, ikz2, iky4, iko, * nkmax, mhalf) call mulmul(io, ix6, iu, isgo, isgx6, isgu, iko, ikx6, iku, * nkmax, mhalf) call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko, * nkmax, mhalf) c call mulmul(ix2, iz4, iv, isgx2, isgz4, isgv, ikx2, ikz4, ikv, * nkmax, mhalf) call mulmul(iv, iy6, iu, isgv, isgy6, isgu, ikv, iky6, iku, * nkmax, mhalf) call muldif(io, iu, iv, isgo, isgu, isgv, iko, iku, ikv, * nkmax, mhalf) c call mulmul(iy2, ix4, io, isgy2, isgx4, isgo, iky2, ikx4, iko, * nkmax, mhalf) call mulmul(io, iz6, iu, isgo, isgz6, isgu, iko, ikz6, iku, * nkmax, mhalf) call muldif(iv, iu, io, isgv, isgu, isgo, ikv, 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 *IQSIGN c c subroutine for determining position of point ifif with respect c to sphere determined by (weighted) points ifir, isec, ithi, ifou c if positive then ifif is outside the sphere c if negative then ifif is inside the sphere c if zero then ifif is in the surface of the sphere c subroutine iqsign(x, y, z, w, x2, y2, z2, w2, ifir, isec, * ithi, ifou, ifif, mhalf, mfull, isclp, * isclw, isclr, delaun, ipout) c integer x(*), y(*), z(*), w(*), x2(*), y2(*), z2(*), w2(*) integer ifir, isec, ithi, ifou, ifif, mhalf, mfull, nkmax, ipout integer isclp(*), isclw(*), isclr(*) parameter (nkmax = 30) integer io(nkmax), iu(nkmax), iv(nkmax), ip(nkmax) integer iq2(nkmax), iq3(nkmax), iq4(nkmax), iq5(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 ixf(nkmax), iyf(nkmax), izf(nkmax) integer ixf2(nkmax), iyf2(nkmax), izf2(nkmax) integer iwf(nkmax), iw2(nkmax), iw3(nkmax), iw4(nkmax), iw5(nkmax) logical delaun integer iwfuw, iwsew, iwthw, iwfow, iwfiw integer ixthw, iythw, izthw, ixfow, iyfow, izfow integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew integer ixfuw, iyfuw, izfuw integer iwfu2, iwse2, iwth2, iwfo2, iwfi2 integer ixth2, iyth2, izth2, ixfo2, iyfo2, izfo2 integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2 integer ixfu2, iyfu2, izfu2 integer isgw2, isgw3, isgw4, isgw5, ikw2, ikw3, ikw4, ikw5 integer isgq2, isgq3, isgq4, isgq5, ikq2, ikq3, ikq4, ikq5 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 isgx5, isgy5, isgz5, ikx5, iky5, ikz5 integer isgo, isgu, isgv, isgp, iko, iku, ikv, ikp integer isgwf, isgcl, ikwf, ikcl c if(delaun) then isgw2 = 0 isgw3 = 0 isgw4 = 0 isgw5 = 0 else iwfuw = w(ifir) iwsew = w(isec) iwthw = w(ithi) iwfow = w(ifou) iwfiw = w(ifif) c iwfu2 = w2(ifir) iwse2 = w2(isec) iwth2 = w2(ithi) iwfo2 = w2(ifou) iwfi2 = w2(ifif) 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) call decmp2(io, isgo, iko, iwfiw, iwfi2, mhalf, mfull, isclw) call muldif(io, iwf, iu, isgo, isgwf, isgu, iko, ikwf, iku, * nkmax, mhalf) call mulmul(iu, isclr, iw5, isgu, isgcl, isgw5, iku, ikcl, * ikw5, 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) ixfiw = x(ifif) iyfiw = y(ifif) izfiw = z(ifif) 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) ixfi2 = x2(ifif) iyfi2 = y2(ifif) izfi2 = z2(ifif) 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 frterm(ixfiw, iyfiw, izfiw, ixf, iyf, izf, isgxf, isgyf, * isgzf, ikxf, ikyf, ikzf, ixf2, iyf2, izf2, * isgxf2, isgyf2, isgzf2, ikxf2, * ikyf2, ikzf2, iw5, ix5, iy5, iz5, iq5, isgw5, isgx5, * isgy5, isgz5, isgq5, ikw5, ikx5, iky5, ikz5, ikq5, * mhalf, mfull, ixfi2, iyfi2, izfi2, isclp) c call mulmul(iq5, ix2, iv, isgq5, isgx2, isgv, ikq5, ikx2, ikv, * nkmax, mhalf) call mulmul(iq5, ix3, iu, isgq5, isgx3, isgu, ikq5, ikx3, iku, * nkmax, mhalf) call mulmul(iq5, ix4, ip, isgq5, isgx4, isgp, ikq5, ikx4, ikp, * nkmax, mhalf) call detrm3(iv, iy2, iz2, isgv, isgy2, isgz2, * iu, iy3, iz3, isgu, isgy3, isgz3, * ip, iy4, iz4, isgp, isgy4, isgz4, * ikv, iku, ikp, iky2, iky3, iky4, * ikz2, ikz3, ikz4, io, isgo, iko, mhalf) 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, iu, isgu, iku, mhalf) call mulmul(iu, ix5, ip, isgu, isgx5, isgp, iku, ikx5, ikp, * nkmax, mhalf) call muldif(io, ip, iv, isgo, isgp, isgv, iko, ikp, ikv, * nkmax, 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, iu, isgu, iku, mhalf) call mulmul(iu, iy5, io, isgu, isgy5, isgo, iku, iky5, iko, * nkmax, mhalf) call muldif(iv, io, ip, isgv, isgo, isgp, ikv, iko, ikp, * nkmax, 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, iu, isgu, iku, mhalf) call mulmul(iu, iz5, iv, isgu, isgz5, isgv, iku, ikz5, ikv, * nkmax, mhalf) call muldif(ip, iv, io, isgp, isgv, isgo, ikp, ikv, iko, * nkmax, mhalf) c ipout = isgo c return end *IQSIG1 c c subroutine for determining position of point ifif with respect to c line segment determined by (weighted) points ifir, isec, assuming c ifir, isec and ifif are colinear c if positive then ifif is in the exterior of the line segment c if negative then ifif is in the interior of the line segment c if zero then ifif is one of the endpoints c subroutine iqsig1(x, y, z, w, x2, y2, z2, w2, ifir, isec, ifif, * mhalf, mfull, isclp, isclw, isclr, delaun, ipout) c integer x(*), y(*), z(*), w(*), x2(*), y2(*), z2(*), w2(*) integer ifir, isec, ifif, mhalf, mfull, nkmax, ipout integer isclp(*), isclw(*), isclr(*) parameter (nkmax = 30) integer io(nkmax), iu(nkmax), iv(nkmax) integer iq2(nkmax), iq5(nkmax) integer ix2(nkmax), iy2(nkmax), iz2(nkmax) integer ix5(nkmax), iy5(nkmax), iz5(nkmax) integer ixf(nkmax), iyf(nkmax), izf(nkmax) integer ixf2(nkmax), iyf2(nkmax), izf2(nkmax) integer iwf(nkmax), iw2(nkmax), iw5(nkmax) logical delaun integer iwfuw, iwsew, iwfiw integer ixfiw, iyfiw, izfiw, ixsew, iysew, izsew integer ixfuw, iyfuw, izfuw integer iwfu2, iwse2, iwfi2 integer ixfi2, iyfi2, izfi2, ixse2, iyse2, izse2 integer ixfu2, iyfu2, izfu2 integer isgw2, isgw5, ikw2, ikw5 integer isgq2, isgq5, ikq2, ikq5 integer isgxf, isgyf, isgzf, ikxf, ikyf, ikzf integer isgxf2, isgyf2, isgzf2, ikxf2, ikyf2, ikzf2 integer isgx2, isgy2, isgz2, ikx2, iky2, ikz2 integer isgx5, isgy5, isgz5, ikx5, iky5, ikz5 integer isgo, isgu, isgv, iko, iku, ikv integer isgwf, isgcl, ikwf, ikcl c if(delaun) then isgw2 = 0 isgw5 = 0 else iwfuw = w(ifir) iwsew = w(isec) iwfiw = w(ifif) c iwfu2 = w2(ifir) iwse2 = w2(isec) iwfi2 = w2(ifif) 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, iwfiw, iwfi2, mhalf, mfull, isclw) call muldif(io, iwf, iu, isgo, isgwf, isgu, iko, ikwf, iku, * nkmax, mhalf) call mulmul(iu, isclr, iw5, isgu, isgcl, isgw5, iku, ikcl, * ikw5, nkmax, mhalf) endif c ixfuw = x(ifir) iyfuw = y(ifir) izfuw = z(ifir) ixsew = x(isec) iysew = y(isec) izsew = z(isec) ixfiw = x(ifif) iyfiw = y(ifif) izfiw = z(ifif) c ixfu2 = x2(ifir) iyfu2 = y2(ifir) izfu2 = z2(ifir) ixse2 = x2(isec) iyse2 = y2(isec) izse2 = z2(isec) ixfi2 = x2(ifif) iyfi2 = y2(ifif) izfi2 = z2(ifif) 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 5205 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(ixfiw, iyfiw, izfiw, ixf, iyf, izf, isgxf, isgyf, * isgzf, ikxf, ikyf, ikzf, ixf2, iyf2, izf2, * isgxf2, isgyf2, isgzf2, ikxf2, * ikyf2, ikzf2, iw5, ix5, iy5, iz5, iq5, isgw5, isgx5, * isgy5, isgz5, isgq5, ikw5, ikx5, iky5, ikz5, ikq5, * mhalf, mfull, ixfi2, iyfi2, izfi2, isclp) c call detrm1(iq5, ix2, iy2, iz2, ix2, iy2, iz2, isgq5, isgx2, * isgy2, isgz2, isgx2, isgy2, isgz2, ikq5, ikx2, iky2, * ikz2, ikx2, iky2, ikz2, iu, isgu, iku, mhalf) c call detrm1(iq2, ix2, iy2, iz2, ix5, iy5, iz5, isgq2, isgx2, * isgy2, isgz2, isgx5, isgy5, isgz5, ikq2, ikx2, iky2, * ikz2, ikx5, iky5, ikz5, iv, isgv, ikv, mhalf) c call muldif(iu, iv, io, isgu, isgv, isgo, iku, ikv, iko, * nkmax, mhalf) c ipout = isgo c return end *IQSIG2 c c subroutine for determining position of point ifif with respect c to circle determined by (weighted) points ifir, isec, ithi, c if positive then ifif is outside the circle c if negative then ifif is inside the circle c if zero then ifif is in the circle c subroutine iqsig2(x, y, z, w, x2, y2, z2, w2, ifir, isec, ithi, * ifif, mhalf, mfull, isclp, isclw, isclr, delaun, ipout) c integer x(*), y(*), z(*), w(*), x2(*), y2(*), z2(*), w2(*) integer ifir, isec, ithi, ifif, mhalf, mfull, nkmax, ipout integer isclp(*), isclw(*), isclr(*) parameter (nkmax = 30) integer io(nkmax), iu(nkmax), iv(nkmax), ip(nkmax) integer iq2(nkmax), iq3(nkmax), iq5(nkmax) integer ix2(nkmax), iy2(nkmax), iz2(nkmax) integer ix3(nkmax), iy3(nkmax), iz3(nkmax) integer ix5(nkmax), iy5(nkmax), iz5(nkmax) integer ixf(nkmax), iyf(nkmax), izf(nkmax) integer ixf2(nkmax), iyf2(nkmax), izf2(nkmax) integer iwf(nkmax), iw2(nkmax), iw3(nkmax), iw5(nkmax) logical delaun integer iwfuw, iwsew, iwthw, iwfiw integer ixthw, iythw, izthw, ixfiw, iyfiw, izfiw integer ixfuw, iyfuw, izfuw, ixsew, iysew, izsew integer iwfu2, iwse2, iwth2, iwfi2 integer ixth2, iyth2, izth2, ixfi2, iyfi2, izfi2 integer ixfu2, iyfu2, izfu2, ixse2, iyse2, izse2 integer isgw2, isgw3, isgw5, ikw2, ikw3, ikw5 integer isgq2, isgq3, isgq5, ikq2, ikq3, ikq5 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 isgx5, isgy5, isgz5, ikx5, iky5, ikz5 integer isgo, isgu, isgv, isgp, iko, iku, ikv, ikp integer isgwf, isgcl, ikwf, ikcl c if(delaun) then isgw2 = 0 isgw3 = 0 isgw5 = 0 else iwfuw = w(ifir) iwsew = w(isec) iwthw = w(ithi) iwfiw = w(ifif) c iwfu2 = w2(ifir) iwse2 = w2(isec) iwth2 = w2(ithi) iwfi2 = w2(ifif) 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, iwfiw, iwfi2, mhalf, mfull, isclw) call muldif(io, iwf, iu, isgo, isgwf, isgu, iko, ikwf, iku, * nkmax, mhalf) call mulmul(iu, isclr, iw5, isgu, isgcl, isgw5, iku, ikcl, * ikw5, 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) ixfiw = x(ifif) iyfiw = y(ifif) izfiw = z(ifif) 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) ixfi2 = x2(ifif) iyfi2 = y2(ifif) izfi2 = z2(ifif) 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 5305 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(ixfiw, iyfiw, izfiw, ixf, iyf, izf, isgxf, isgyf, * isgzf, ikxf, ikyf, ikzf, ixf2, iyf2, izf2, * isgxf2, isgyf2, isgzf2, ikxf2, * ikyf2, ikzf2, iw5, ix5, iy5, iz5, iq5, isgw5, isgx5, * isgy5, isgz5, isgq5, ikw5, ikx5, iky5, ikz5, ikq5, * mhalf, mfull, ixfi2, iyfi2, izfi2, isclp) c call detrm0(iq5, iy2, iz2, iy3, iz3, isgq5, isgy2, isgz2, * isgy3, isgz3, ikq5, iky2, ikz2, iky3, ikz3, * iv, isgv, ikv, mhalf) c call detrm0(iq5, iz2, ix2, iz3, ix3, isgq5, isgz2, isgx2, * isgz3, isgx3, ikq5, ikz2, ikx2, ikz3, ikx3, * iu, isgu, iku, mhalf) c call detrm0(iq5, ix2, iy2, ix3, iy3, isgq5, isgx2, isgy2, * isgx3, isgy3, ikq5, ikx2, iky2, ikx3, iky3, * ip, isgp, ikp, mhalf) c call detrm3(iv, ix2, ix3, isgv, isgx2, isgx3, * iu, iy2, iy3, isgu, isgy2, isgy3, * ip, iz2, iz3, isgp, isgz2, isgz3, * ikv, iku, ikp, ikx2, iky2, ikz2, * ikx3, iky3, ikz3, io, isgo, iko, mhalf) c call detrm2(iq2, ix2, iy2, iz2, isgq2, isgx2, isgy2, isgz2, * iq3, ix3, iy3, iz3, isgq3, isgx3, isgy3, isgz3, * ikq2, ikx2, iky2, ikz2, ikq3, ikx3, iky3, ikz3, * iu, isgu, iku, mhalf) call mulmul(iu, ix5, ip, isgu, isgx5, isgp, iku, ikx5, ikp, * nkmax, mhalf) call muldif(io, ip, iv, isgo, isgp, isgv, iko, ikp, ikv, * nkmax, mhalf) c call detrm2(iq2, iy2, iz2, ix2, isgq2, isgy2, isgz2, isgx2, * iq3, iy3, iz3, ix3, isgq3, isgy3, isgz3, isgx3, * ikq2, iky2, ikz2, ikx2, ikq3, iky3, ikz3, ikx3, * iu, isgu, iku, mhalf) call mulmul(iu, iy5, io, isgu, isgy5, isgo, iku, iky5, iko, * nkmax, mhalf) call muldif(iv, io, ip, isgv, isgo, isgp, ikv, iko, ikp, * nkmax, mhalf) c call detrm2(iq2, iz2, ix2, iy2, isgq2, isgz2, isgx2, isgy2, * iq3, iz3, ix3, iy3, isgq3, isgz3, isgx3, isgy3, * ikq2, ikz2, ikx2, iky2, ikq3, ikz3, ikx3, iky3, * iu, isgu, iku, mhalf) call mulmul(iu, iz5, iv, isgu, isgz5, isgv, iku, ikz5, ikv, * nkmax, mhalf) call muldif(ip, iv, io, isgp, isgv, isgo, ikp, ikv, iko, * nkmax, mhalf) c ipout = isgo 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 *DETRM0 c subroutine detrm0(iq, ix2, iy2, ix3, iy3, isgq, isgx2, isgy2, * isgx3, isgy3, ikq, ikx2, iky2, ikx3, iky3, * io, isgo, iko, mhalf) c integer nkmax parameter (nkmax = 30) integer io(*), iu(nkmax), iv(nkmax), iw(nkmax) integer iq(*), ix2(*), iy2(*), ix3(*), iy3(*) integer isgq, isgx2, isgy2, isgx3, isgy3 integer ikq, ikx2, iky2, ikx3, iky3, isgo, iko, mhalf integer isgu, isgv, isgw, iku, ikv, ikw c call mulmul(iq, ix2, iv, isgq, isgx2, isgv, ikq, ikx2, ikv, * nkmax, mhalf) call mulmul(iv, iy3, iu, isgv, isgy3, isgu, ikv, iky3, iku, * nkmax, mhalf) call mulmul(iq, iy2, iv, isgq, isgy2, isgv, ikq, iky2, ikv, * nkmax, mhalf) call mulmul(iv, ix3, iw, isgv, isgx3, isgw, ikv, ikx3, ikw, * nkmax, mhalf) call muldif(iu, iw, io, isgu, isgw, isgo, iku, ikw, iko, * nkmax, mhalf) c return end *DETRM1 c subroutine detrm1(iq, ix2, iy2, iz2, ix3, iy3, iz3, * isgq, isgx2, isgy2, isgz2, isgx3, isgy3, isgz3, * ikq, ikx2, iky2, ikz2, ikx3, iky3, ikz3, * io, isgo, iko, mhalf) c integer nkmax parameter (nkmax = 30) integer io(*), iv(nkmax), iw(nkmax) integer iq(*), ix2(*), iy2(*), iz2(*), ix3(*), iy3(*), iz3(*) integer isgq, isgx2, isgy2, isgz2, isgx3, isgy3, isgz3 integer ikq, ikx2, iky2, ikz2, ikx3, iky3, ikz3 integer isgo, isgv, isgw, iko, ikv, ikw, mhalf c call mulmul(iq, ix2, iv, isgq, isgx2, isgv, ikq, ikx2, ikv, * nkmax, mhalf) call mulmul(iv, ix3, io, isgv, isgx3, isgo, ikv, ikx3, iko, * nkmax, mhalf) c call mulmul(iq, iy2, iv, isgq, isgy2, isgv, ikq, iky2, ikv, * nkmax, mhalf) call mulmul(iv, iy3, iw, isgv, isgy3, isgw, ikv, iky3, ikw, * nkmax, mhalf) isgw =-isgw call muldif(io, iw, iv, isgo, isgw, isgv, iko, ikw, ikv, * nkmax, mhalf) c call mulmul(iq, iz2, io, isgq, isgz2, isgo, ikq, ikz2, iko, * nkmax, mhalf) call mulmul(io, iz3, iw, isgo, isgz3, isgw, iko, ikz3, ikw, * nkmax, mhalf) isgw =-isgw call muldif(iv, iw, io, isgv, isgw, isgo, ikv, ikw, iko, * nkmax, mhalf) c return end *DETRM2 c subroutine detrm2(iq2, ix2, iy2, iz2, isgq2, isgx2, isgy2, isgz2, * iq3, ix3, iy3, iz3, isgq3, isgx3, isgy3, isgz3, * ikq2, ikx2, iky2, ikz2, ikq3, ikx3, iky3, ikz3, * io, isgo, iko, mhalf) c integer nkmax parameter (nkmax = 30) integer io(*), ip(nkmax), ir(nkmax), iv(nkmax), iw(nkmax) integer iq2(*), iq3(*) integer ix2(*), iy2(*), iz2(*), ix3(*), iy3(*), iz3(*) integer isgq2, isgx2, isgy2, isgz2, isgq3, isgx3, isgy3, isgz3 integer ikq2, ikx2, iky2, ikz2, ikq3, ikx3, iky3, ikz3 integer isgo, isgp, isgr, isgv, isgw, iko, ikp, ikr, ikv, ikw integer mhalf c call detrm0(iq2, ix2, iy2, ix3, iy3, isgq2, isgx2, isgy2, * isgx3, isgy3, ikq2, ikx2, iky2, ikx3, iky3, * iv, isgv, ikv, mhalf) call mulmul(iv, iy3, io, isgv, isgy3, isgo, ikv, iky3, iko, * nkmax, mhalf) c call detrm0(iq2, iz2, ix2, iz3, ix3, isgq2, isgz2, isgx2, * isgz3, isgx3, ikq2, ikz2, ikx2, ikz3, ikx3, * iv, isgv, ikv, mhalf) call mulmul(iv, iz3, ip, isgv, isgz3, isgp, ikv, ikz3, ikp, * nkmax, mhalf) c call muldif(io, ip, ir, isgo, isgp, isgr, iko, ikp, ikr, * nkmax, mhalf) c call detrm0(iq3, ix2, iy2, ix3, iy3, isgq3, isgx2, isgy2, * isgx3, isgy3, ikq3, ikx2, iky2, ikx3, iky3, * iv, isgv, ikv, mhalf) call mulmul(iv, iy2, io, isgv, isgy2, isgo, ikv, iky2, iko, * nkmax, mhalf) c call detrm0(iq3, iz2, ix2, iz3, ix3, isgq3, isgz2, isgx2, * isgz3, isgx3, ikq3, ikz2, ikx2, ikz3, ikx3, * iv, isgv, ikv, mhalf) call mulmul(iv, iz2, ip, isgv, isgz2, isgp, ikv, ikz2, ikp, * nkmax, mhalf) c call muldif(io, ip, iw, isgo, isgp, isgw, iko, ikp, ikw, * nkmax, mhalf) call muldif(ir, iw, io, isgr, isgw, isgo, ikr, ikw, iko, * 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