summaryrefslogtreecommitdiff
path: root/Dragon/src/NXTLCA.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/NXTLCA.f')
-rw-r--r--Dragon/src/NXTLCA.f809
1 files changed, 809 insertions, 0 deletions
diff --git a/Dragon/src/NXTLCA.f b/Dragon/src/NXTLCA.f
new file mode 100644
index 0000000..9469b2c
--- /dev/null
+++ b/Dragon/src/NXTLCA.f
@@ -0,0 +1,809 @@
+*DECK NXTLCA
+ FUNCTION NXTLCA(IPRINT,ITST ,NDIM ,MXMESH,LINMAX,
+ > MESH ,ORITRK,DIRTRK,DCMESH,
+ > NBCOR ,NBSINT,ISINT ,TRKLSI)
+*
+*----------
+*
+*Purpose:
+* To track a Cartesian 2-D or 3-D geometry
+* using the NXT tracking procedure.
+*
+*Copyright:
+* Copyright (C) 2005 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+*Author(s): G. Marleau.
+*
+*Parameters: input
+* IPRINT print level.
+* ITST type of tracking, where:
+* =-1 only the exact geometry
+* is considered taking into account the
+* submesh in each direction;
+* = 0 only the global geometry
+* is considered without taking into account the
+* submesh in each direction;
+* = 1 both the global
+* geometry (as a first step) and the exact geometry
+* are considered taking into account the
+* submesh in each direction.
+* NDIM dimension of problem.
+* MXMESH maximum number of spatial subdivision in
+* $X$, $Y$ or $Z$.
+* LINMAX maximum number of segments in a track.
+* MESH effective number of spatial subdivision in
+* each direction ($X$, $Y$ and $Z$).
+* ORITRK a point on the track (origin).
+* DIRTRK the track direction (director cosines).
+* DCMESH spatial description of the parallepiped.
+*
+*Parameters: output
+* NXTLCA number of surfaces intersected.
+* NBCOR number of corner found for each external faces.
+* NBSINT number of surface crossed by track.
+* ISINT direction of plane intersected and
+* the surfaces crossed by the track.
+* TRKLSI the surface intersection distance.
+*
+*Reference:
+* G. Marleau,
+* New Geometries Processing in DRAGON: The NXT: Module,
+* Report IGE-260, Polytechnique Montreal,
+* Montreal, 2005.
+*
+*----------
+*
+ IMPLICIT NONE
+*----
+* Subroutine arguments
+*----
+ INTEGER IPRINT,ITST,NDIM,MXMESH,LINMAX
+ INTEGER MESH(NDIM)
+ DOUBLE PRECISION ORITRK(NDIM),
+ > DIRTRK(NDIM),
+ > DCMESH(-1:MXMESH,5)
+ INTEGER NBCOR(2),NBSINT
+ INTEGER ISINT(0:5,LINMAX)
+ DOUBLE PRECISION TRKLSI(LINMAX)
+ INTEGER NXTLCA
+*----
+* Local parameters
+*----
+ INTEGER IOUT
+ CHARACTER NAMSBR*6
+ PARAMETER (IOUT=6,NAMSBR='NXTLCA')
+ DOUBLE PRECISION DCUTOF,DZERO,DONE,DTWO
+ PARAMETER (DCUTOF=1.0D-8,DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
+*----
+* Local variables
+*----
+ INTEGER IDIR,INEXT,IFRST,ILAST,IFACE,IOFF,NBTINT
+ INTEGER JDIR,JNEXT,JFRST,JLAST,JFACE,JOFF,JSUR
+ INTEGER KDIR,KNEXT,KFRST,KLAST,KFACE,KOFF,KSUR
+ INTEGER IC1,IDIRS,IDIRF,IDIRB,IRFIN(5)
+ INTEGER ISDIR,JSDIR,KSDIR
+ DOUBLE PRECISION DP1,DP2,TRKK,TRKJ
+ DOUBLE PRECISION TRKDIS,TRKOLD,DELTKD
+*----
+* Data
+*----
+ CHARACTER CDIR(4)*1
+ SAVE CDIR
+ DATA CDIR /'X','Y','Z','R'/
+*----
+* Verify ITST option and reset to default value if invalid
+*----
+ IF(ITST .LT. -1 .OR. ITST .GT. 1) THEN
+*----
+* Reset ITST=1 (complete analysis) if the value of ITST is invalid
+*----
+ ITST=1
+ ENDIF
+ KFRST=1
+ KLAST=1
+ KNEXT=1
+ TRKK=DZERO
+ TRKOLD=DZERO
+ KOFF=1
+*----
+* Initialise output vectors
+*----
+ NBCOR(1)=0
+ NBCOR(2)=0
+ NBSINT=0
+ ISINT(0:5,:LINMAX)=0
+ TRKLSI(:LINMAX)=DZERO
+*----
+* Print header if required
+*----
+ IF(IPRINT .GT. 1000) THEN
+ WRITE(IOUT,6000) NAMSBR
+ WRITE(IOUT,6011) 'meshx={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,1),IC1=0,MESH(1))
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'meshy={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,2),IC1=0,MESH(2))
+ WRITE(IOUT,6013)
+ IF(NDIM .EQ. 3) THEN
+ WRITE(IOUT,6011) 'meshz={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,3),IC1=0,MESH(3))
+ WRITE(IOUT,6013)
+ ENDIF
+ WRITE(IOUT,6011) 'trackorigin={ '
+ WRITE(IOUT,6012) ORITRK
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'trackdirection={ '
+ WRITE(IOUT,6012) DIRTRK
+ WRITE(IOUT,6013)
+ ENDIF
+*----
+* For option ITST=0, 1, consider global geometry
+* and test if line crosses this geometry.
+*----
+ IF(ITST .EQ. 0 .OR. ITST .EQ. 1) THEN
+*----
+* Scan over Cartesian directions
+*----
+ DO 100 IDIR=1,NDIM
+*----
+* IDIR is main planes direction ($X$, $Y$ or $Z$)
+* JDIR is first direction of plane perpendicular
+* to main direction ($Y, $Z$ or $X$).
+* KDIR is second direction of plane perpendicular
+* to main direction ($Z$, $X$ or $Y$).
+*----
+ JDIR=MOD(IDIR,NDIM)+1
+ KDIR=MOD(IDIR+1,NDIM)+1
+*----
+* Select planes order in direction IDIR
+*----
+ IF(DIRTRK(IDIR) .LT. 0) THEN
+ INEXT=-1
+ IFRST=MESH(IDIR)+1
+ ILAST=1
+ ELSE
+ INEXT=1
+ IFRST=1
+ ILAST=MESH(IDIR)+1
+ ENDIF
+*----
+* Select planes order in direction JDIR
+*----
+ IF(DIRTRK(JDIR) .LT. 0) THEN
+ JNEXT=-1
+ JFRST=MESH(JDIR)+1
+ JLAST=1
+ ELSE
+ JNEXT=1
+ JFRST=1
+ JLAST=MESH(JDIR)+1
+ ENDIF
+*----
+* Select planes order in direction KDIR (if required)
+*----
+ IF(NDIM .EQ. 3) THEN
+ IF(DIRTRK(KDIR) .LT. 0) THEN
+ KNEXT=-1
+ KFRST=MESH(KDIR)+1
+ KLAST=1
+ ELSE
+ KNEXT=1
+ KFRST=1
+ KLAST=MESH(KDIR)+1
+ ENDIF
+ ENDIF
+*----
+* Scan over planes in direction IDIR
+*----
+ DO 101 IFACE=IFRST,ILAST,ILAST-IFRST
+*----
+* Compute track length required to reach a face
+* and intersection point with infinite plane
+*----
+ TRKDIS=(DCMESH(IFACE-1,IDIR)-ORITRK(IDIR))/DIRTRK(IDIR)
+ TRKJ=DBLE(JNEXT)*(TRKDIS*DIRTRK(JDIR)+ORITRK(JDIR))
+ IF(NDIM .EQ. 3) THEN
+ TRKK=DBLE(KNEXT)*(TRKDIS*DIRTRK(KDIR)+ORITRK(KDIR))
+ ENDIF
+*----
+* Test if point is in finite Cartesian plane (j,k)
+*----
+ DP1=DBLE(JNEXT)*DCMESH(JFRST-1,JDIR)
+ DP2=DBLE(JNEXT)*DCMESH(JLAST-1,JDIR)
+ IF(TRKJ .GE. DP1-DCUTOF
+ > .AND. TRKJ .LE. DP2+DCUTOF ) THEN
+ IF(NDIM .EQ. 3) THEN
+*----
+* For 3-D, consider intersection with plane
+*----
+ DP1=DBLE(KNEXT)*DCMESH(KFRST-1,KDIR)
+ DP2=DBLE(KNEXT)*DCMESH(KLAST-1,KDIR)
+ IF(TRKK .GE. DP1-DCUTOF
+ > .AND. TRKK .LE. DP2+DCUTOF ) THEN
+*----
+* If no point, save distance set NBSINT to 1 and continue
+* if one point already known, check if this point is
+* at a different distance exit, otherwise, this point
+* is already identified (intersection of 2 or more planes)
+* just continue
+*----
+ IF(NBSINT .EQ. 0) THEN
+ NBSINT=1
+ TRKOLD=TRKDIS
+ ELSE IF(TRKOLD .NE. TRKDIS) THEN
+ NBSINT=2
+ GO TO 105
+ ENDIF
+ ENDIF
+ ELSE
+*----
+* For 2-D, intersection with a line is sufficient
+*----
+ IF(NBSINT .EQ. 0) THEN
+ NBSINT=1
+ TRKOLD=TRKDIS
+ ELSE IF(TRKOLD .NE. TRKDIS) THEN
+ NBSINT=2
+ GO TO 105
+ ENDIF
+ ENDIF
+ ENDIF
+ 101 CONTINUE
+ 100 CONTINUE
+ 105 CONTINUE
+ ELSE
+*----
+* If ITST=-1, assume that 2 surface are intersected
+* and continue
+*----
+ NBSINT=2
+ ENDIF
+ NXTLCA=NBSINT
+*----
+* If line does not intersect any surfaces, exit
+*----
+ IF(NXTLCA .EQ. 0) THEN
+ IF(IPRINT .GT. 1000) THEN
+ WRITE(IOUT,6002) NAMSBR
+ ENDIF
+ RETURN
+ ENDIF
+*----
+* If line does not intersect 2 surfaces, abort
+*----
+ IF(NBSINT .NE. 2) THEN
+ IF(NBSINT .EQ. 1) CALL XABORT(NAMSBR//
+ > ': Invalid line since it intersects only 1 surface')
+ IF(IPRINT .GT. 1000)
+ > WRITE(IOUT,9000) NBSINT
+ ENDIF
+*----
+* If ITST=0, return number of surfaces intersected by line
+*----
+ IF(ITST .EQ. 0) RETURN
+*----
+* If ITST=1 and NXTLCA=2 or ITST=-1
+* track geometry with submesh
+*----
+ NBSINT=0
+*----
+* Scan over Cartesian directions
+*----
+ DO 110 IDIR=1,NDIM
+*----
+* IDIR is main planes direction ($X$, $Y$ or $Z$)
+* JDIR is first direction of plane perpendicular
+* to main direction ($Y, $Z$ or $X$).
+* KDIR is second direction of plane perpendicular
+* to main direction ($Z$, $X$ or $Y$).
+*----
+ JDIR=MOD(IDIR,NDIM)+1
+ KDIR=MOD(IDIR+1,NDIM)+1
+*----
+* Select planes order in direction IDIR
+*----
+ IF(DIRTRK(IDIR) .LT. 0) THEN
+ INEXT=-1
+ IFRST=MESH(IDIR)+1
+ ILAST=1
+ IOFF=-1
+ ELSE
+ INEXT=1
+ IFRST=1
+ ILAST=MESH(IDIR)+1
+ IOFF=0
+ ENDIF
+*----
+* Select planes order in direction JDIR
+*----
+ IF(DIRTRK(JDIR) .LT. 0) THEN
+ JNEXT=-1
+ JFRST=MESH(JDIR)+1
+ JLAST=2
+ JOFF=-1
+ ELSE
+ JNEXT=1
+ JFRST=1
+ JLAST=MESH(JDIR)
+ JOFF=0
+ ENDIF
+*----
+* Select planes order in direction KDIR (if required)
+*----
+ IF(NDIM .EQ. 3) THEN
+ IF(DIRTRK(KDIR) .LT. 0) THEN
+ KNEXT=-1
+ KFRST=MESH(KDIR)+1
+ KLAST=2
+ KOFF=-1
+ ELSE
+ KNEXT=1
+ KFRST=1
+ KLAST=MESH(KDIR)
+ KOFF=0
+ ENDIF
+ ENDIF
+*----
+* Scan over planes in direction IDIR
+*----
+ DO 111 IFACE=IFRST,ILAST,INEXT
+ ISDIR=IDIR
+ IF(IFACE .EQ. IFRST .OR. IFACE .EQ. ILAST)
+ > ISDIR=-ISDIR
+*----
+* Compute track length required to reach a face
+* and intersection point with infinite plane
+*----
+ TRKDIS=(DCMESH(IFACE-1,IDIR)-ORITRK(IDIR))/DIRTRK(IDIR)
+ TRKJ=DBLE(JNEXT)*(TRKDIS*DIRTRK(JDIR)+ORITRK(JDIR))
+ IF(NDIM .EQ. 3) THEN
+ TRKK=DBLE(KNEXT)*(TRKDIS*DIRTRK(KDIR)+ORITRK(KDIR))
+ ENDIF
+*----
+* Test if point is in finite Cartesian plane (j,k) and add to number
+* of surfaces crossed if this is the case
+* (j=1 or NPLJ+1 and k=1 or NPLK+1)
+* For values of j=1,NPLJ and k=1,NPLK a volume was crossed
+*----
+ DO 112 JFACE=JFRST,JLAST,JNEXT
+ JSDIR=JDIR
+ IF(JFACE .EQ. JFRST .OR. JFACE .EQ. JLAST)
+ > JSDIR=-JSDIR
+ DP1=DBLE(JNEXT)*DCMESH(JFACE-1,JDIR)
+ DP2=DBLE(JNEXT)*DCMESH(JFACE+JNEXT-1,JDIR)
+ IF(TRKJ .GE. DP1-DCUTOF
+ > .AND. TRKJ .LE. DP2+DCUTOF ) THEN
+ IF(NDIM .EQ. 3) THEN
+*----
+* For 3-D, consider intersection with plane
+*----
+ DO 113 KFACE=KFRST,KLAST,KNEXT
+ KSDIR=KDIR
+ IF(KFACE .EQ. KFRST .OR. KFACE .EQ. KLAST)
+ > KSDIR=-KSDIR
+ DP1=DBLE(KNEXT)*DCMESH(KFACE-1,KDIR)
+ DP2=DBLE(KNEXT)*DCMESH(KFACE+KNEXT-1,KDIR)
+ IF(TRKK .GE. DP1-DCUTOF
+ > .AND. TRKK .LE. DP2+DCUTOF ) THEN
+ NBSINT=NBSINT+1
+ IF(IPRINT .GT. 1000) THEN
+ WRITE(IOUT,6011) 'Intersection at '
+ WRITE(IOUT,6014) IDIR,IFACE,
+ > JDIR,JFACE+JOFF,KDIR,KFACE+KOFF
+ WRITE(IOUT,6015) TRKDIS,DIRTRK
+ ENDIF
+ IF(NBSINT .EQ. 1) THEN
+ ISINT(0,NBSINT)=ISDIR
+ ISINT(IDIR,NBSINT)=-IFACE
+ ISINT(JDIR,NBSINT)=JFACE+JOFF
+ ISINT(KDIR,NBSINT)=KFACE+KOFF
+ TRKLSI(NBSINT)=TRKDIS
+ ELSE
+*----
+* Scan previous distances and reorder if necessary
+*----
+ NBTINT=NBSINT
+ DO 120 JSUR=NBTINT-1,1,-1
+ DELTKD=TRKDIS-TRKLSI(JSUR)
+ IF(DELTKD .GE. -DCUTOF) THEN
+*----
+* For more than one intersection at a given distance:
+* a) if this is an external surface, store the distance
+* for corner duplication
+* b) if this is an internal surface select adequate next region
+* to cross and reset NBSINT=NBSINT-1
+*----
+ IF(ABS(DELTKD) .LE. DCUTOF) THEN
+ IF(ISDIR .GT. 0 ) THEN
+ IF(ISINT(0,JSUR) .GT. 0) THEN
+ ISINT(0,JSUR)=ISINT(0,JSUR)+10*ISDIR
+ ISINT(IDIR,JSUR)=-IFACE
+ ENDIF
+ NBSINT=NBSINT-1
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'SurInt A =',JSUR,IFACE,JFACE,KFACE,
+* >ISINT(0,JSUR),ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),
+* >ISINT(KDIR,JSUR),
+* >TRKDIS,TRKLSI(JSUR)
+ GO TO 125
+ ELSE
+ IF(ISINT(0,JSUR) .GT. 0) THEN
+ ISINT(0,JSUR)=ISDIR
+ ISINT(IDIR,JSUR)=-IFACE
+ ISINT(JDIR,JSUR)=JFACE+JOFF
+ ISINT(KDIR,JSUR)=KFACE+KOFF
+ NBSINT=NBSINT-1
+ GO TO 125
+ ENDIF
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'SurInt B =',JSUR,IFACE,JFACE,KFACE,
+* >ISINT(0,JSUR),ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),
+* >ISINT(KDIR,JSUR),
+* >TRKDIS,TRKLSI(JSUR)
+ ENDIF
+ ENDIF
+ DO 121 KSUR=NBTINT,JSUR+2,-1
+ ISINT(0,KSUR)=ISINT(0,KSUR-1)
+ ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+ ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+ ISINT(KDIR,KSUR)=ISINT(KDIR,KSUR-1)
+ TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 121 CONTINUE
+ ISINT(0,JSUR+1)=ISDIR
+ ISINT(IDIR,JSUR+1)=-IFACE
+ ISINT(JDIR,JSUR+1)=JFACE+JOFF
+ ISINT(KDIR,JSUR+1)=KFACE+KOFF
+ TRKLSI(JSUR+1)=TRKDIS
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'SurInt C =',JSUR,IFACE,JFACE,KFACE,
+* >ISINT(0,JSUR),ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),
+* >ISINT(KDIR,JSUR),
+* >TRKDIS,TRKLSI(JSUR)
+ GO TO 125
+ ENDIF
+ 120 CONTINUE
+ DO 122 KSUR=NBTINT,2,-1
+ ISINT(0,KSUR)=ISINT(0,KSUR-1)
+ ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+ ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+ ISINT(KDIR,KSUR)=ISINT(KDIR,KSUR-1)
+ TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 122 CONTINUE
+ ISINT(0,1)=ISDIR
+ ISINT(IDIR,1)=-IFACE
+ ISINT(JDIR,1)=JFACE+JOFF
+ ISINT(KDIR,1)=KFACE+KOFF
+ TRKLSI(1)=TRKDIS
+ 125 CONTINUE
+ ENDIF
+ ENDIF
+ 113 CONTINUE
+ ELSE
+*----
+* For 2-D, intersection with a line is sufficient
+*----
+ NBSINT=NBSINT+1
+ IF(IPRINT .GT. 1000) THEN
+ WRITE(IOUT,6011) 'Intersection at '
+ WRITE(IOUT,6014) IDIR,IFACE,JDIR,JFACE
+ WRITE(IOUT,6015) TRKDIS,DIRTRK
+ ENDIF
+ IF(NBSINT .EQ. 1) THEN
+ ISINT(0,NBSINT)=ISDIR
+ ISINT(IDIR,NBSINT)=-IFACE
+ ISINT(JDIR,NBSINT)=JFACE+JOFF
+ TRKLSI(NBSINT)=TRKDIS
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'Int sur A =',NBSINT,
+* >ISINT(IDIR,NBSINT),ISINT(JDIR,NBSINT),TRKLSI(NBSINT)
+ ELSE
+*----
+* Scan previous distances and reorder if necessary
+*----
+ NBTINT=NBSINT
+ DO 130 JSUR=NBTINT-1,1,-1
+ DELTKD=TRKDIS-TRKLSI(JSUR)
+ IF(DELTKD .GE. -DCUTOF) THEN
+*----
+* For more than one intersection at a given distance:
+* a) if this is an external surface, store the distance
+* for corner duplication
+* b) if this is an internal surface select adequate next region
+* to cross and reset NBSINT=NBSINT-1
+*----
+ IF(ABS(DELTKD) .LE. DCUTOF) THEN
+ IF(ISDIR .GT. 0 ) THEN
+ IF(ISINT(0,JSUR) .GT. 0) THEN
+ ISINT(0,JSUR)=ISINT(0,JSUR)+10*ISDIR
+ ISINT(IDIR,JSUR)=-IFACE
+ ENDIF
+ NBSINT=NBSINT-1
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'Int sur A =',JSUR,
+* >ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),TRKLSI(JSUR)
+ GO TO 135
+ ELSE
+ IF(ISINT(0,JSUR) .GT. 0) THEN
+ ISINT(0,JSUR)=ISDIR
+ ISINT(IDIR,JSUR)=-IFACE
+ ISINT(JDIR,JSUR)=JFACE+JOFF
+ NBSINT=NBSINT-1
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'Int sur B =',JSUR,
+* >ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),TRKLSI(JSUR)
+ GO TO 135
+ ENDIF
+ ENDIF
+ ENDIF
+ DO 131 KSUR=NBTINT,JSUR+2,-1
+ ISINT(0,KSUR)=ISINT(0,KSUR-1)
+ ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+ ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+ TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 131 CONTINUE
+ ISINT(0,JSUR+1)=ISDIR
+ ISINT(IDIR,JSUR+1)=-IFACE
+ ISINT(JDIR,JSUR+1)=JFACE+JOFF
+ TRKLSI(JSUR+1)=TRKDIS
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'Int sur C =',JSUR+1,
+* >ISINT(IDIR,JSUR+1),ISINT(JDIR,JSUR+1),TRKLSI(JSUR+1)
+ GO TO 135
+ ENDIF
+ 130 CONTINUE
+ DO 132 KSUR=NBTINT,2,-1
+ ISINT(0,KSUR)=ISINT(0,KSUR-1)
+ ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+ ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+ TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 132 CONTINUE
+ ISINT(0,1)=ISDIR
+ ISINT(IDIR,1)=-IFACE
+ ISINT(JDIR,1)=JFACE+JOFF
+ TRKLSI(1)=TRKDIS
+* IF(IPRINT .GT. 1000)
+* >write(6,*) 'Int sur D =',1,
+* >ISINT(IDIR,1),ISINT(JDIR,1),TRKLSI(1)
+ 135 CONTINUE
+ ENDIF
+ ENDIF
+ ENDIF
+ 112 CONTINUE
+ 111 CONTINUE
+ 110 CONTINUE
+ IF(IPRINT .GT. 1000) THEN
+ WRITE(IOUT,6011) 'Intersections'
+ DO 210 JSUR=1,NBSINT
+ WRITE(IOUT,6010) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5)
+ 210 CONTINUE
+ ENDIF
+ IF(NBSINT .GE. 2) THEN
+*----
+* Find number of corners for initial surface intersections
+*----
+ TRKDIS=TRKLSI(1)
+ NBCOR(1)=1
+ DO 140 JSUR=2,NBSINT
+ DELTKD=TRKLSI(JSUR)-TRKDIS
+ IF(DELTKD .GT. DCUTOF) GO TO 145
+ NBCOR(1)=NBCOR(1)+1
+ 140 CONTINUE
+ 145 CONTINUE
+*----
+* Find number of corners for final surface intersections
+*----
+ TRKDIS=TRKLSI(NBSINT)
+ NBCOR(2)=1
+ DO 150 JSUR=NBSINT-1,1,-1
+ DELTKD=TRKDIS-TRKLSI(JSUR)
+ IF(DELTKD .GT. DCUTOF) GO TO 155
+ NBCOR(2)=NBCOR(2)+1
+ 150 CONTINUE
+ 155 CONTINUE
+ ENDIF
+*----
+* Print intersection points
+*----
+ IF(IPRINT .GT. 1000) THEN
+ WRITE(IOUT,6011) 'Initial planes '
+ DO JSUR=1,NBCOR(1)
+ WRITE(IOUT,6010) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5)
+ ENDDO
+ IF(NBCOR(1)+1 .LE. NBSINT-NBCOR(2)) THEN
+ WRITE(IOUT,6011) 'Intermediate planes '
+ DO JSUR=NBCOR(1)+1,NBSINT-NBCOR(2)
+ WRITE(IOUT,6010) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5)
+ ENDDO
+ ENDIF
+ WRITE(IOUT,6011) 'Final planes '
+ DO JSUR=NBSINT-NBCOR(2)+1,NBSINT
+ WRITE(IOUT,6010) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5)
+ ENDDO
+ ENDIF
+*----
+* Identify final faces
+*----
+ DO JSUR=NBSINT,NBSINT-NBCOR(2)+1,-1
+ DO IDIRS=0,5
+ ISINT(IDIRS,JSUR+1)=ISINT(IDIRS,JSUR)
+ ENDDO
+ TRKLSI(JSUR+1)=TRKLSI(JSUR)
+ IDIRS=ABS(ISINT(0,JSUR+1))
+ IF(IDIRS .GT. 10) CALL XABORT(NAMSBR//
+ >': Final outer surface identified as an internal corner')
+ IDIRF=ABS(ISINT(IDIRS,JSUR+1))
+ IF(IDIRF .EQ. 1) THEN
+ ISINT(IDIRS,JSUR+1)=-1
+ ELSE IF(IDIRF .EQ. MESH(IDIRS)+1) THEN
+ ISINT(IDIRS,JSUR+1)=-2
+ ELSE
+*----
+* Print information for debug
+*----
+ WRITE(IOUT,6011) 'meshx={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,1),IC1=0,MESH(1))
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'meshy={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,2),IC1=0,MESH(2))
+ WRITE(IOUT,6013)
+ IF(NDIM .EQ. 3) THEN
+ WRITE(IOUT,6011) 'meshz={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,3),IC1=0,MESH(3))
+ WRITE(IOUT,6013)
+ ENDIF
+ WRITE(IOUT,6011) 'trackorigin={ '
+ WRITE(IOUT,6012) ORITRK
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'trackdirection={ '
+ WRITE(IOUT,6012) DIRTRK
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'Initial planes '
+ DO KSUR=1,NBCOR(1)
+ WRITE(IOUT,6010) TRKLSI(KSUR),
+ > (ISINT(IDIR,KSUR),IDIR=0,5)
+ ENDDO
+ IF(NBCOR(1)+1 .LE. NBSINT-NBCOR(2)) THEN
+ WRITE(IOUT,6011) 'Intermediate planes '
+ DO KSUR=NBCOR(1)+1,NBSINT-NBCOR(2)
+ WRITE(IOUT,6010) TRKLSI(KSUR),
+ > (ISINT(IDIR,KSUR),IDIR=0,5)
+ ENDDO
+ ENDIF
+ WRITE(IOUT,6011) 'Final planes '
+ DO KSUR=NBSINT-NBCOR(2)+1,NBSINT
+ WRITE(IOUT,6010) TRKLSI(KSUR),
+ > (ISINT(IDIR,KSUR),IDIR=0,5)
+ ENDDO
+ WRITE(IOUT,9001) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5),IDIRF
+ CALL XABORT(NAMSBR//': Invalid final '//CDIR(IDIRS)//
+ >' directed surface')
+ ENDIF
+ ENDDO
+*----
+* Identify regions
+*----
+ DO JSUR=NBSINT-NBCOR(2)+1,NBCOR(1)+1,-1
+ DO IDIRS=1,5
+ IRFIN(IDIRS)=ABS(ISINT(IDIRS,JSUR-1))
+ IF(ISINT(IDIRS,JSUR) .LT. 0 ) THEN
+ IF(ISINT(IDIRS,JSUR-1) .LT. 0 ) THEN
+ IRFIN(IDIRS)=-MAX(ISINT(IDIRS,JSUR-1),ISINT(IDIRS,JSUR))
+ ENDIF
+ ELSE
+ IF(ISINT(IDIRS,JSUR-1) .LT. 0 ) THEN
+ IRFIN(IDIRS)=ISINT(IDIRS,JSUR)
+ ENDIF
+ ENDIF
+ ENDDO
+ TRKLSI(JSUR)=TRKLSI(JSUR)-TRKLSI(JSUR-1)
+ DO IDIRS=1,5
+ IF(ISINT(IDIRS,JSUR) .LT. 0 ) THEN
+ ISINT(IDIRS,JSUR)=IRFIN(IDIRS)
+ ELSE
+ ISINT(IDIRS,JSUR)=ISINT(IDIRS,JSUR)
+ ENDIF
+ ENDDO
+ ENDDO
+*----
+* Identify initial face
+*----
+ DO JSUR=1,NBCOR(1)
+ IDIRS=ABS(ISINT(0,JSUR))
+ IF(IDIRS .GT. 10) CALL XABORT(NAMSBR//
+ >': Initial outer surface identified as an internal corner')
+ IDIRB=ABS(ISINT(IDIRS,JSUR))
+ IF(IDIRB .EQ. 1) THEN
+ ISINT(IDIRS,JSUR)=-1
+ ELSE IF(IDIRB .EQ. MESH(IDIRS)+1) THEN
+ ISINT(IDIRS,JSUR)=-2
+ ELSE
+*----
+* Print information for debug
+*----
+ WRITE(IOUT,6011) 'meshx={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,1),IC1=0,MESH(1))
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'meshy={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,2),IC1=0,MESH(2))
+ WRITE(IOUT,6013)
+ IF(NDIM .EQ. 3) THEN
+ WRITE(IOUT,6011) 'meshz={ '
+ WRITE(IOUT,6012) (DCMESH(IC1,3),IC1=0,MESH(3))
+ WRITE(IOUT,6013)
+ ENDIF
+ WRITE(IOUT,6011) 'trackorigin={ '
+ WRITE(IOUT,6012) ORITRK
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'trackdirection={ '
+ WRITE(IOUT,6012) DIRTRK
+ WRITE(IOUT,6013)
+ WRITE(IOUT,6011) 'Initial planes '
+ DO KSUR=1,NBCOR(1)
+ WRITE(IOUT,6010) TRKLSI(KSUR),
+ > (ISINT(IDIR,KSUR),IDIR=0,5)
+ ENDDO
+ IF(NBCOR(1)+1 .LE. NBSINT-NBCOR(2)) THEN
+ WRITE(IOUT,6011) 'Intermediate planes '
+ DO KSUR=NBCOR(1)+1,NBSINT-NBCOR(2)
+ WRITE(IOUT,6010) TRKLSI(KSUR),
+ > (ISINT(IDIR,KSUR),IDIR=0,5)
+ ENDDO
+ ENDIF
+ WRITE(IOUT,6011) 'Final planes '
+ DO KSUR=NBSINT-NBCOR(2)+1,NBSINT
+ WRITE(IOUT,6010) TRKLSI(KSUR),
+ > (ISINT(IDIR,KSUR),IDIR=0,5)
+ ENDDO
+ WRITE(IOUT,9001) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5),IDIRB
+ CALL XABORT(NAMSBR//': Invalid initial '//CDIR(IDIRS)//
+ >' directed surface')
+ ENDIF
+ ENDDO
+*----
+* Print final track information
+*----
+ IF(IPRINT .GT. 1000) THEN
+ WRITE(IOUT,6011) 'Initial face '
+ DO JSUR=1,NBCOR(1)
+ WRITE(IOUT,6010) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5)
+ ENDDO
+ WRITE(IOUT,6011) 'Regions '
+ DO JSUR=NBCOR(1)+1,NBSINT-NBCOR(2)+1
+ WRITE(IOUT,6010) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5)
+ ENDDO
+ WRITE(IOUT,6011) 'Final face '
+ DO JSUR=NBSINT-NBCOR(2)+2,NBSINT+1
+ WRITE(IOUT,6010) TRKLSI(JSUR),
+ > (ISINT(IDIR,JSUR),IDIR=0,5)
+ ENDDO
+ WRITE(IOUT,6001) NAMSBR
+ ENDIF
+ RETURN
+*----
+* Output formats
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows ')
+ 6001 FORMAT(' Output from --',A6,'-- completed *)')
+ 6002 FORMAT(' No region crossed'/
+ > ' Output from --',A6,'-- completed *)')
+ 6010 FORMAT(1X,F25.16,6I10)
+ 6011 FORMAT(A20)
+ 6012 FORMAT(6(1X,F25.16,:,','))
+ 6013 FORMAT('};')
+ 6014 FORMAT(6I10)
+ 6015 FORMAT(6(F25.16,2X))
+ 9000 FORMAT(' Warning : ',I10,' surfaces crossed')
+ 9001 FORMAT(1X,F25.16,7I10)
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:24:09 +0000