path: root/Dragon/src/NXTTLC.f

*DECK NXTTLC
      SUBROUTINE NXTTLC(IPTRK ,IFTEMP,IPRINT,IGTRK ,NDIM  ,MAXMSH,
     >                  NFSUR ,NFREG ,NUCELL,NBUCEL,NBANGL,
     >                  LINMAX,MXGSUR,MXGREG,RENO  ,NBDR  ,ITYPBC,
     >                  IFMT  ,MAXPIN,AZMQUA,IPER  ,IUNFLD,SURVOL,
     >                  DGMESH,DANGLT,DDENWT,DNSANG,NBSANG,DEPART,
     >                  MAXSUB,MAXSGL,NTLINE,DVNOR ,DSNOR)
*
*-----------------------------------------------------------------------
*
*Purpose:
* To generate the cyclic tracking lines (specular tracking)
* for a geometry.
*
*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
* IPTRK   pointer to the TRACKING data structure in
*         update or creation mode.
* IFTEMP  pointer to a temporary TRACKING file in
*         update or creation mode.
* IPRINT  print level.
* IGTRK   flag to generate the tracking file. In the case where
*         IGTRK=1, the tracking is performed and
*         used to evaluate the track normalisation factor and the
*         tracking file is generated. When IGTRK=0, the tracking is
*         still performed and used to evaluate the
*         track normalisation factor but the tracking file is not
*         generated.
* NDIM    problem dimensions.
* MAXMSH  maximum number of elements in MESH array.
* NFSUR   number of surfaces.
* NFREG   number of regions.
* NUCELL  number of cell after unfolding in
*         $X$, $Y$ and $Z$ directions.
* NBUCEL  number of cells in unfolded geometry.
* NBANGL  number of angles.
* LINMAX  maximum number of segments in a track.
* MXGSUR  maximum number of surfaces for any geometry.
* MXGREG  maximum number of region for any geometry.
* RENO    track normalisation option. A value RENO=-1 implies
*         a direction dependent normalization of the tracks
*         for the volume while a value RENO=0, implies
*         a global normalisation.
* NBDR    number of directions for track normalization.
* ITYPBC  type of boundary conditions where:
*         =0 for geometry with Cartesianb oundaries;
*         =1 for geometry with annular boundary;
*         =2 for geometry with hexagonal boundary.
* IFMT    tracking file format:
*         =0 for short file;
*         =1 long file required by TLM:.
* MAXPIN  maximum number of pins in a cell.
* AZMQUA  tracking type.
* IPER    cell periodicity factor in each direction.
* IUNFLD  description of unfolded geometry.
* SURVOL  global surface volume vector.
* DGMESH  meshing vector for global geometry.
* DANGLT  angles.
* DDENWT  angular density for each angle.
* DNSANG  spatial density required.
* NBSANG  number of segments for each angles.
* DEPART  track starting point.
*
*Parameters: output
* MAXSUB  maximum number of subtracks in a line.
* MAXSGL  maximum number of segments in a line.
* NTLINE  total number of lines generated.
* DVNOR   ratio of analytic to tracked volume.
* DSNOR   ratio of analytic to tracked surface area.
*
*Reference:
*  G. Marleau,
*  New Geometries Processing in DRAGON: The NXT: Module,
*  Report IGE-260, Polytechnique Montreal,
*  Montreal, 2005.
*  \\\\
*  Based on the XELTI2 and XELTI3 routines.
*
*----------
*
      USE              GANLIB
      IMPLICIT         NONE
*----
*  Subroutine arguments
*----
      TYPE(C_PTR)      IPTRK
      INTEGER          IFTEMP
      INTEGER          IPRINT,IGTRK,NDIM,MAXMSH,NFSUR,NFREG,
     >                 NUCELL(3),NBUCEL,NBANGL,LINMAX,MXGSUR,MXGREG,
     >                 RENO,NBDR,ITYPBC,IFMT,MAXPIN,AZMQUA,IPER(3),
     >                 MAXSUB,MAXSGL,NTLINE
      INTEGER          IUNFLD(2,NBUCEL)
      DOUBLE PRECISION SURVOL(-NFSUR:NFREG),DGMESH(-1:MAXMSH,4),
     >                 DANGLT(NDIM,NBANGL),DDENWT(NBANGL),
     >                 DNSANG(NBANGL)
      INTEGER          NBSANG(5,NBANGL)
      DOUBLE PRECISION DEPART(NDIM,2,NBANGL)
      DOUBLE PRECISION DVNOR(NFREG,NBDR),
     >                 DSNOR(NFSUR)
*----
*  Local parameters
*----
      INTEGER          IOUT
      CHARACTER        NAMSBR*6
      PARAMETER       (IOUT=6,NAMSBR='NXTTLC')
      DOUBLE PRECISION DCUTOF,DZERO,DONE,DTWO
      PARAMETER       (DCUTOF=1.0D-8,DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
*----
*  Functions
*----
      DOUBLE PRECISION XDRCST,PI
      INTEGER          NXTLCA,NXTLHA,IRLA
*----
*  Local variables
*----
      INTEGER          LMAXT,ITST,IANGL,NGLINE,IDIR,ITXY(3),
     >                 NSCAN,NUMANG,NPOINT,NBPTS,ISCAN,
     >                 ISTART,ISUM,NBCOR(2),IFREF
      INTEGER          NBSINT,ISDIR,IBLIN,IELIN,ICI,IX,IY,IZ,
     >                 IOX,IOY,IOZ,IOC,ICEL,ITRN,ISURF
      DOUBLE PRECISION FACVOL,FACSUR,DWGT,WGTFAC,DAWGT,PMAX,PMIN,
     >                 DSTART,RADT2,ANGT2
      DOUBLE PRECISION ANGLED(3),ANGLES(3),ANGLEN(3),ANGLER(3),
     >                 TRKORI(3),TRKOR2(3),TRKORD(3),TRKORN(3),
     >                 TRKLIM(2),CELLPO(3,2)
      DOUBLE PRECISION DERR,DCERR,DSVERR,DMVERR,DAVERR,DSSERR,
     >                 DMSERR,DASERR,VCONTA,VCONT
      INTEGER          ISEG,JSEG,IREG,ILREG,ISUR,NBVERR,NBSERR,
     >                 NBV0,NBV1,NBS0,IPRINL,IPRINC
      INTEGER          ISD,NSDEB,ISF,NSFIN,NBREG,ISBL,NTSEG
      INTEGER          ISINT,NBSEG,IPT,JSDIR,ISDIRX,IPTN
      DOUBLE PRECISION DHALF
      DOUBLE PRECISION AAA
      INTEGER          IND,II,JJ
      LOGICAL          LNEW
*----
*  Allocatable arrays (local)
*----
      INTEGER, ALLOCATABLE, DIMENSION(:) :: NUMERO
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ICINT
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: LENGTH,DLENGT,
     >                                               DCINT
      INTEGER, ALLOCATABLE, DIMENSION(:) :: KANGL
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: TORIG
*----
*  Data
*----
      CHARACTER        CDIR(4)*1
      SAVE             CDIR
      DATA             CDIR /'X','Y','Z','R'/
*----
* Allocate: temporary storage (local)
* NUMERO  region/surface identification number
*         for segment.
* LENGTH  segment length.
* DLENGT  spatial location of each line segment.
* ICINT   identification of spatial position for each
*         line segment in cell description of geometry.
* DCINT   position of each intersection point for each
*         line segment in cell description of geometry.
* TORIG   track origin
* KANGL   angle index (quadrant)
*----
      ALLOCATE(NUMERO(LINMAX),LENGTH(LINMAX),DLENGT(LINMAX),
     > ICINT(0:5,2*(NBUCEL+4)),DCINT(2*(NBUCEL+4)))
*----
*  Processing starts:
*  print routine openning output header if required
*  and initialize various parameters.
*----
      IF(IPRINT .GE. 2) THEN
        WRITE(IOUT,6000) NAMSBR
      ENDIF
      PI=XDRCST('Pi',' ')
      DHALF=DONE/DTWO
      MAXSUB=0
      MAXSGL=0
      NTLINE=0
      LMAXT=2*(NBUCEL+4)
      ITST=RENO
      AAA=AZMQUA
      ITST=1
      DVNOR(:NFREG,:NBDR)=DZERO
      DSNOR(:NFSUR)=DZERO
      FACVOL=DTWO
      IF(NDIM .EQ. 2) THEN
        FACSUR=DTWO*PI
      ELSE
        FACSUR=DTWO*DTWO
      ENDIF
      IF(IPRINT .GT. 1000) THEN
        DO IDIR=1,NDIM
          WRITE(IOUT,6070) 'DGMESH '//CDIR(IDIR),NUCELL(IDIR)
          WRITE(IOUT,6071) (DGMESH(IX,IDIR),IX=1,NUCELL(IDIR))
        ENDDO
      ENDIF
*----
*  Loop over directions
*----
      IPRINL=IPRINT
      IPRINC=IPRINT
      NGLINE=0
      DO IANGL=1,NBANGL
        DO IDIR=1,NDIM
          ITXY(IDIR)=NBSANG(IDIR,IANGL)
        ENDDO
        ITXY(3)=ITXY(1)*ITXY(2)
        NSCAN=NBSANG(3,IANGL)
        NUMANG=NBSANG(4,IANGL)
        NPOINT=NBSANG(5,IANGL)
        ALLOCATE(TORIG(NDIM,NUMANG),KANGL(NUMANG))
*----
*  Loop over quadrant if both quadrant not scanned simultaneously
*  If NSCAN = 2  0 TO PI/2 AND PI/2 TO PI are scanned simultaneously
*  If NSCAN = 1  0 TO PI/2 AND PI/2 TO PI are scanned independently
*----
        NBPTS=NSCAN*NPOINT
        DO ISCAN=2,NSCAN,-1
*----
*  Initialize weight
*  Store track direction in ANGLES
*  and track starting point in TRKOR2
*----
          DWGT=DNSANG(IANGL)
          DAWGT=2.0D0*DWGT*DDENWT(IANGL)
          DO IDIR=1,NDIM
            ANGLED(IDIR)=DANGLT(IDIR,IANGL)
            TRKOR2(IDIR)=DEPART(IDIR,1,IANGL)
          ENDDO
          IF(ISCAN .EQ. 1) THEN
            IF(ITXY(1) .EQ. 0) THEN
              ANGLED(2)=-ANGLED(2)
            ELSE IF(ITXY(2) .EQ. 0) THEN
              ANGLED(1)=-ANGLED(1)
            ENDIF
          ENDIF
          DO IDIR=1,NDIM
            ANGLES(IDIR)=ANGLED(IDIR)
            ANGLEN(IDIR)=ANGLED(IDIR)
            ANGLER(IDIR)=ANGLED(IDIR)
            TRKORD(IDIR)=TRKOR2(IDIR)
            TRKORI(IDIR)=TRKORD(IDIR)
            TRKORN(IDIR)=TRKORD(IDIR)
          ENDDO
*----
*  Loop over points
*----
          IF(IPRINT .GE. 10) THEN
            WRITE(IOUT,6050) IANGL,ISCAN,ITXY,DWGT,
     >      (DANGLT(IDIR,IANGL),IDIR=1,NDIM),
     >      (DEPART(IDIR,1,IANGL),IDIR=1,NDIM)
          ENDIF
          ISTART=0
          ISUM=ISTART
          IFREF=1
          DO IPT=1,NBPTS
*----
*  Find if this track crosses geometry
*----
            IF(ISTART .EQ. 0) THEN
              IF(IFREF .EQ. 1) THEN
                DO IDIR=1,NDIM
                  ANGLES(IDIR)=ANGLER(IDIR)
                  ANGLED(IDIR)=ANGLER(IDIR)
                ENDDO
              ENDIF
            ENDIF
            ISINT=1
            DCINT(ISINT)=DZERO
            IRLA=-1
            IF(ITYPBC .EQ. 0) THEN
*----
*  Cartesian assembly
*----
              IRLA=NXTLCA(IPRINC,ITST  ,NDIM  ,MAXMSH,LMAXT,
     >                     NUCELL,TRKORI,ANGLES,DGMESH,
     >                     NBCOR ,NBSINT,ICINT ,DCINT)
            ELSE IF(ITYPBC .EQ. 1) THEN
*----
*  Annular assembly
*----
              CALL XABORT(NAMSBR//': Circular BC not implemented')
            ELSE IF(ITYPBC .EQ. 2) THEN
*----
*  Hexagonal assembly
*----
              IRLA=NXTLHA(IPRINC,ITST  ,NDIM  ,MAXMSH,LMAXT,
     >                     NUCELL,TRKORI,ANGLES,DGMESH,
     >                     NBCOR ,NBSINT,ICINT ,DCINT)
            ENDIF
            IF(IRLA .EQ. -1) CALL XABORT(NAMSBR//
     >': This type of cell cannot be tracked by NXT:')
*----
*  When no external face crossed go to next line
*----
            IF(IRLA .EQ. 0) GO TO 105
            ISDIRX=ABS(ICINT(0,1))
            ISDIR=ABS(ICINT(0,NBSINT))
            JSDIR=MOD(ISDIR,2)+1
            IF(ISTART .EQ. 0) THEN
*----
*  Define position of cyclic starting point
*----
              NGLINE=NGLINE+1
              DO IDIR=1,NDIM
                TRKORD(IDIR)=TRKORD(IDIR)+DCINT(ISINT)*ANGLES(IDIR)
                TRKORI(IDIR)=TRKORD(IDIR)
              ENDDO
              DO IDIR=1,NBCOR(1)
                DCINT(IDIR)=DZERO
              ENDDO
              IF(ISCAN .EQ. 1 .AND. ITXY(3) .NE. 0) THEN
                ANGLED(ISDIRX)=-ANGLED(ISDIRX)
                DO IDIR=1,NDIM
                  ANGLES(IDIR)=ANGLED(IDIR)
                ENDDO
                IF(ITYPBC .EQ. 0) THEN
*----
*  Cartesian assembly
*----
                  IRLA=NXTLCA(IPRINL,ITST  ,NDIM  ,MAXMSH,LMAXT,
     >                         NUCELL,TRKORI,ANGLES,DGMESH,
     >                         NBCOR ,NBSINT,ICINT ,DCINT)
                ELSE IF(ITYPBC .EQ. 1) THEN
*----
*  Annular assembly
*----
                  CALL XABORT(NAMSBR//': Circular BC not implemented')
                ELSE IF(ITYPBC .EQ. 2) THEN
*----
*  Hexagonal assembly
*----
                  IRLA=NXTLHA(IPRINL,ITST  ,NDIM  ,MAXMSH,LMAXT,
     >                         NUCELL,TRKORI,ANGLES,DGMESH,
     >                         NBCOR ,NBSINT,ICINT ,DCINT)
                ELSE
                  CALL XABORT(NAMSBR//
     >': This type of cell cannot be tracked by NXT:')
                ENDIF
                IF(IRLA .EQ. 0) GO TO 105
                DO IDIR=1,NDIM
                  TRKORD(IDIR)=TRKORD(IDIR)+DCINT(ISINT)*ANGLES(IDIR)
                  TRKORI(IDIR)=TRKORD(IDIR)
                ENDDO
                DO IDIR=1,NBCOR(1)
                  DCINT(IDIR)=DZERO
                ENDDO
                ISDIR=ABS(ICINT(0,NBSINT))
                JSDIR=MOD(ISDIR,2)+1
              ENDIF
              DO IDIR=1,NDIM
                ANGLES(IDIR)=ANGLED(IDIR)
                ANGLEN(IDIR)=ANGLED(IDIR)
              ENDDO
            ENDIF
*----
*  For track crossing geometry, end find surface direction
*  and end of track location
*----
            DO IDIR=1,NDIM
              AAA=DZERO
              DO ISINT=NBCOR(1),NBSINT-NBCOR(2)+1
                AAA=AAA+DCINT(ISINT)*ANGLES(IDIR)
              ENDDO
              TRKORN(IDIR)=TRKORI(IDIR)+AAA
              TORIG(IDIR,ISTART+1)=TRKORI(IDIR)
            ENDDO
            KANGL(ISTART+1)=0
            DO II=1,4*NBANGL
              IF((ANGLEN(1).EQ.DANGLT(1,II)).AND.
     >           (ANGLEN(2).EQ.DANGLT(2,II))) THEN
                  KANGL(ISTART+1)=II
                  GO TO 10
              ENDIF
            ENDDO
            CALL XABORT(NAMSBR//
     >': Unable to find an angular index for a subtrack')
 10         CONTINUE
            IPRINL=IPRINT
*            IF(NTLINE+1 .EQ. 1411
*     >    .OR. NTLINE+1 .EQ. 1414
*     >    .OR. NTLINE+1 .EQ. 1415
*     >                            ) IPRINL=IPRINT+4000
            IF(ISTART.EQ.0) THEN
              IF(IPRINL .GE. 500) THEN
                WRITE(IOUT,6051) NTLINE+1
              ENDIF
            ENDIF
            IF(IPRINL .GE. 500) THEN
              WRITE(IOUT,6052) TRKORI,ANGLEN
            ENDIF
*----
*  Find location with respect to end surfaces positions
*----
            PMAX=ABS(TRKORN(ISDIR)-DGMESH(NUCELL(ISDIR),ISDIR))
            PMIN=ABS(TRKORN(ISDIR)-DGMESH(0,ISDIR))
            IF(IPER(ISDIR) .EQ. 1) THEN
*----
*  Translate next initial starting point if possible
*----
              IF(PMAX .LT. DCUTOF) THEN
                TRKORN(ISDIR)=DGMESH(0,ISDIR)
              ELSE IF(PMIN .LT. DCUTOF) THEN
                TRKORN(ISDIR)=DGMESH(NUCELL(ISDIR),ISDIR)
              ELSE
                WRITE(IOUT,9100) NAMSBR,ISDIR,NUCELL(ISDIR),
     >          TRKORI(ISDIR),DGMESH(NUCELL(ISDIR),ISDIR),PMAX
                CALL XABORT(NAMSBR//
     >          ': Translation before cell limit attempted')
              ENDIF
            ELSE
*----
*  Change track direction upon reflexion if possible
*----
              IF(PMAX .GT. DCUTOF .AND. PMIN .GT. DCUTOF) THEN
                CALL XABORT(NAMSBR//
     >          ': Reflexion before cell limit attempted')
              ENDIF
              ANGLEN(ISDIR)=-ANGLEN(ISDIR)
            ENDIF
*----
*  Track seems ok, locate track segments for geometry
*----
            IF(ISTART .EQ. 0) THEN
*----
*  For first series of line segments initialize track information
*----
              IBLIN=1
              IELIN=0
            ELSE
*----
*  Next series of line segments, add at the end of
*  current line
*----
              IBLIN=IELIN+1
              DLENGT(IBLIN)=DZERO
            ENDIF
            ISTART=ISTART+1
            ISUM=ISUM+1
            NUMERO(IBLIN)=0
            DLENGT(IBLIN)=DCINT(1)
            DSTART=DCINT(1)
            IBLIN=IBLIN+1
*----
*  Find line segments for this track series
*----
            DO ICI=NBCOR(1)+1,NBSINT-NBCOR(2)+1
              IX=ICINT(1,ICI)
              IY=ICINT(2,ICI)
              IOX=IX
              CELLPO(1,2)=DGMESH(IX,1)
              CELLPO(1,1)=DGMESH(IX-1,1)
              CELLPO(2,2)=DGMESH(IY,2)
              CELLPO(2,1)=DGMESH(IY-1,2)
              IOY=(IY-1)*NUCELL(1)
              IOZ=0
              IF(NDIM .EQ. 3) THEN
                IZ=ICINT(3,ICI)
                IOZ=(IZ-1)*NUCELL(1)*NUCELL(2)
                CELLPO(3,2)=DGMESH(IZ,3)
                CELLPO(3,1)=DGMESH(IZ-1,3)
              ENDIF
              IOC=IOX+IOY+IOZ
              ICEL=IUNFLD(1,IOC)
              ITRN=IUNFLD(2,IOC)
              TRKLIM(1)=DSTART
              DSTART=DSTART+DCINT(ICI)
              TRKLIM(2)=DSTART
              IF(ICI .EQ. NBCOR(1)+1) THEN
*----
*  initial surfaces (at TRKLIM(1)) considered
*----
                ISURF=-1
              ELSE IF(ICI .EQ. NBSINT-NBCOR(2)) THEN
*----
*  final surfaces (at TRKLIM(2)) considered
*----
                ISURF=1
              ELSE
*----
*  no surface considered
*----
                ISURF=0
              ENDIF
*----
*  Track turned Cell
*----
              CALL NXTTCR(IPTRK ,IPRINL,ICEL  ,ITRN  ,ISURF ,
     >                    NDIM  ,MAXMSH,LINMAX,MXGSUR,MXGREG,
     >                    MAXPIN,CELLPO,TRKLIM,TRKORI,ANGLES,
     >                    IBLIN ,IELIN ,NUMERO,DLENGT)
              IPRINL=IPRINT
              AAA=MAX(ABS(TRKLIM(1)),ABS(TRKLIM(2)))
              DERR=ABS(DLENGT(IELIN)-DSTART)
              IF(AAA .GT. DONE) DERR=DERR/AAA
              IF(DERR .GT. DCUTOF) THEN
                IF(AAA .GT. DONE) THEN
                  WRITE(IOUT,9101) NAMSBR,NGLINE,
     >            DLENGT(IELIN),DSTART,DERR
                ELSE
                  WRITE(IOUT,9102) NAMSBR,NGLINE,
     >            DLENGT(IELIN),DSTART,DERR
                ENDIF
                CALL XABORT(NAMSBR//
     >': End of track does not coincide with end of cell')
              ENDIF
              IBLIN=IELIN+1
*----
*  End of tracking for this line segment
*----
            ENDDO
            IPRINL=IPRINT
*----
*  This series of line segments completed
*  Test if the cycle is back to starting point
*----
            RADT2=DZERO
            ANGT2=DZERO
            DO IDIR=1,NDIM
              RADT2=RADT2+(TRKORD(IDIR)-TRKORN(IDIR))**2
              ANGT2=ANGT2+(ANGLED(IDIR)-ANGLEN(IDIR))**2
            ENDDO
            RADT2=SQRT(RADT2)
            ANGT2=SQRT(ANGT2)
            IF(RADT2 .GT. DCUTOF .OR. ANGT2 .GT. DCUTOF) THEN
*----
*  Cycle is incomplete
*  Check if period is higher than predicted
*  reset starting point and direction
*----
              IF(ISTART .GE. NUMANG) CALL XABORT(NAMSBR//
     >        ': Cyclic period is too long')
                IF(IPRINL .GE. 500) THEN
                  IF(IPT .EQ. NBPTS) WRITE(IOUT,6064)
                ENDIF
              DO IDIR=1,NDIM
                TRKORI(IDIR)=TRKORN(IDIR)
                ANGLES(IDIR)=ANGLEN(IDIR)
              ENDDO
            ELSE
*----
*  Cycle is complete
*  Check if period is ok
*----
              NTLINE=NTLINE+1
*----
*  Process tracking line and save
*----
              NBSEG=IELIN
*----
*  Compress tracking vector for region with DLENGT=0.0,
*  for multiple surface of intersection and
*  for successive segments in the same region.
*----
              ISD=-1
              JSEG=0
              ILREG=-1
              DO ISEG=1,NBSEG
                IREG=NUMERO(ISEG)
                IF(DLENGT(ISEG) .GT. DZERO) THEN
                  IF(IREG .GT. 0) THEN
                    IF(IREG .EQ. ILREG) THEN
                      DLENGT(JSEG)=DLENGT(JSEG)+DLENGT(ISEG)
                    ELSE
                      JSEG=JSEG+1
                      DLENGT(JSEG)=DLENGT(ISEG)
                      NUMERO(JSEG)=IREG
                      ILREG=IREG
                    ENDIF
                    ISD=IREG
                  ELSE
                    ISF=NUMERO(ISEG+1)
                    IF(ISD .LT. 0) THEN
                      IF(ISF .GT. 0) THEN
                        JSEG=JSEG+1
                        DLENGT(JSEG)=DLENGT(ISEG)*DHALF
                        NUMERO(JSEG)=IREG
                        ILREG=-1
                      ENDIF
                    ELSE
                      JSEG=JSEG+1
                      DLENGT(JSEG)=DLENGT(ISEG)*DHALF
                      NUMERO(JSEG)=IREG
                      ILREG=-1
                    ENDIF
                    ISD=IREG
                  ENDIF
                ENDIF
              ENDDO
              NBSEG=JSEG
*----
*  Add contribution of track to volume integration for this angle
*  in this quadrant and compress tracking line by removing
*  segments with NUMERO=0.
*----
              JSEG=0
              NSDEB=0
              NBREG=0
              IND=0
              LNEW=.TRUE.
              DO ISEG=1,NBSEG
                IREG=NUMERO(ISEG)
                IF(DLENGT(ISEG) .GT. DZERO) THEN
                  IF(IREG .GT. 0) THEN
                    IF(LNEW) THEN
                      IND=IND+1
                      IF(IND.GT.ISTART) CALL XABORT(NAMSBR//
     >                ': ISTART overflow')
                      LNEW=.FALSE.
                    ENDIF
                    NBREG=NBREG+1
                    VCONT=DLENGT(ISEG)*DWGT*FACVOL
                    DVNOR(IREG,1)=DVNOR(IREG,1)+VCONT
                    IF(NBDR .GT. 1) THEN
                      II=KANGL(IND)
                      IF(II .GT. 2*NBANGL) THEN
                        JJ=II-2*NBANGL
                      ELSE
                        JJ=II+2*NBANGL
                      ENDIF
                      IF(DANGLT(1,II).EQ.DZERO) THEN
                        VCONTA=DHALF*DLENGT(ISEG)*DAWGT
                      ELSE
                        VCONTA=DLENGT(ISEG)*DAWGT
                      ENDIF
                      DVNOR(IREG,II+1)=DVNOR(IREG,II+1)+VCONTA
                      DVNOR(IREG,JJ+1)=DVNOR(IREG,JJ+1)+VCONTA
                    ENDIF
                    JSEG=JSEG+1
                    LENGTH(JSEG)=DLENGT(ISEG)
                    NUMERO(JSEG)=IREG
                  ELSE IF(IREG .LT. 0) THEN
                    LNEW=.TRUE.
                    IF(NBREG .EQ. 0) THEN
                      NSDEB=NSDEB+1
                    ENDIF
                    DSNOR(-IREG)=DSNOR(-IREG)+DLENGT(ISEG)*DWGT*FACSUR
                    JSEG=JSEG+1
                    LENGTH(JSEG)=DLENGT(ISEG)
                    NUMERO(JSEG)=IREG
                  ENDIF
                ENDIF
              ENDDO
              NBSEG=JSEG
              NSFIN=0
              DO ISEG=NBSEG,1,-1
                IREG=NUMERO(ISEG)
                IF(DLENGT(ISEG) .GT. DZERO) THEN
                  IF(IREG .GT. 0) THEN
                    GO TO 115
                  ELSE IF(IREG .LT. 0) THEN
                    NSFIN=NSFIN+1
                  ENDIF
                ENDIF
              ENDDO
 115          CONTINUE
*----
*  Matlab commands to create a simili-TLM plot
*----
              IF(IPRINL .GE. 2000) THEN
                WRITE(IOUT,6064)
                WRITE(IOUT,6063) NTLINE,'LineReg'
                WRITE(IOUT,6061) (NUMERO(ISEG),ISEG=1,NBSEG)
                WRITE(IOUT,6064)
                WRITE(IOUT,6063) NTLINE,'LinePos'
                WRITE(IOUT,6062) (LENGTH(ISEG),ISEG=1,NBSEG)
                WRITE(IOUT,6064)
                WRITE(IOUT,6065)
              ENDIF
*----
*  Store line on temporary tracking file if required
*----
              WGTFAC=DONE/DBLE(NSDEB*NSFIN)
              MAXSUB=MAX(MAXSUB,ISTART)
              MAXSGL=MAX(MAXSGL,NBSEG)
              NTSEG=NBSEG-NSDEB-NSFIN+2
              ISBL=0
              DO ISD=1,NSDEB
                DO ISF=0,NSFIN-1
                  ISBL=ISBL+1
                  IF(NSDEB*NSFIN .GT. 1 .AND. IPRINL .GE. 1000)
     >            WRITE(IOUT,6026) NTLINE-1+ISBL,NTLINE,ISBL
                  IF(IGTRK .EQ. 1) THEN
                    IF(IFMT == 1) THEN
                      WRITE(IFTEMP) ISTART,NTSEG,DWGT*WGTFAC,
     >                (KANGL(II),II=1,ISTART),
     >                NUMERO(ISD),
     >                (NUMERO(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                NUMERO(NBSEG-ISF),
     >                LENGTH(ISD),
     >                (LENGTH(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                LENGTH(NBSEG-ISF),
     >                NTLINE,1,1,1,
     >                ((TORIG(IDIR,II),IDIR=1,NDIM),II=1,ISTART)
                    ELSE
                      WRITE(IFTEMP) ISTART,NTSEG,DWGT*WGTFAC,
     >                (KANGL(II),II=1,ISTART),
     >                NUMERO(ISD),
     >                (NUMERO(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                NUMERO(NBSEG-ISF),
     >                LENGTH(ISD),
     >                (LENGTH(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                LENGTH(NBSEG-ISF)
                    ENDIF
                  ENDIF
                ENDDO
              ENDDO
              IF(ISBL .EQ. 0) THEN
                WRITE(IOUT,6027) NTLINE,IANGL,NBSEG,
     >                           NSDEB,NSFIN,NBREG,
     >                           NUMERO(1),DLENGT(1),
     >                           NUMERO(NBSEG),DLENGT(NBSEG)
                ISBL=1
                IF(IFMT == 1) THEN
                  WRITE(IFTEMP) ISTART,NTSEG,DWGT*WGTFAC,
     >            (KANGL(II),II=1,ISTART),
     >            NUMERO(1),
     >            (NUMERO(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >            NUMERO(NBSEG),
     >            LENGTH(1),
     >            (LENGTH(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >            LENGTH(NBSEG),
     >            NTLINE,1,1,1,
     >            ((TORIG(IDIR,II),IDIR=1,NDIM),II=1,ISTART)
                ELSE
                  WRITE(IFTEMP) ISTART,NTSEG,DWGT*WGTFAC,
     >            (KANGL(II),II=1,ISTART),
     >            NUMERO(1),
     >            (NUMERO(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >            NUMERO(NBSEG),
     >            LENGTH(1),
     >            (LENGTH(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >            LENGTH(NBSEG)
                ENDIF
              ENDIF
              NTLINE=NTLINE-1+ISBL
              IF(ISTART .EQ. NUMANG) THEN
*----
*  Displace line to next starting point
*----
                DO IDIR=1,NDIM
                  ANGLES(IDIR)=ANGLED(IDIR)
                  ANGLEN(IDIR)=ANGLED(IDIR)
                  TRKOR2(IDIR)=TRKOR2(IDIR)+DEPART(IDIR,2,IANGL)
                  TRKORD(IDIR)=TRKOR2(IDIR)
                  TRKORI(IDIR)=TRKORD(IDIR)
                  TRKORN(IDIR)=TRKORD(IDIR)
                ENDDO
              ELSE
                IF(NSCAN .EQ. 2) THEN
                  IF(IPER(ISDIR) .EQ. 1) THEN
                    IF(IFREF .EQ. -1) THEN
                      IFREF=1
                      DO IDIR=1,NDIM
                        ANGLES(IDIR)=ANGLED(IDIR)
                        ANGLEN(IDIR)=ANGLED(IDIR)
                        TRKOR2(IDIR)=TRKOR2(IDIR)+DEPART(IDIR,2,IANGL)
                        TRKORD(IDIR)=TRKOR2(IDIR)
                        TRKORI(IDIR)=TRKORD(IDIR)
                        TRKORN(IDIR)=TRKORD(IDIR)
                      ENDDO
                    ELSE
                      ANGLED(ISDIR)=-ANGLED(ISDIR)
                      IFREF=-1
                    ENDIF
                  ELSE IF(IPER(JSDIR) .EQ. 1) THEN
                    IF(IFREF .EQ. 1) THEN
                      DO IPTN=1,NBPTS
                        DO IDIR=1,NDIM
                          TRKORI(IDIR)=TRKORN(IDIR)
                          ANGLES(IDIR)=ANGLEN(IDIR)
                        ENDDO
*----
*  Find next surface intersection track crosses geometry
*----
                        IF(ITYPBC .EQ. 0) THEN
*----
*  Cartesian assembly
*----
                          IRLA=NXTLCA(IPRINL,ITST  ,NDIM  ,MAXMSH,LMAXT,
     >                               NUCELL,TRKORI,ANGLES,DGMESH,
     >                               NBCOR ,NBSINT,ICINT ,DCINT)
                        ELSE IF(ITYPBC .EQ. 1) THEN
                          CALL XABORT(NAMSBR//
     >                    ': Circular BC not implemented')
                        ELSE IF(ITYPBC .EQ. 2) THEN
*----
*  Hexagonal assembly
*----
                          IRLA=NXTLHA(IPRINL,ITST  ,NDIM  ,MAXMSH,LMAXT,
     >                               NUCELL,TRKORI,ANGLES,DGMESH,
     >                               NBCOR ,NBSINT,ICINT ,DCINT)
                        ELSE
                          CALL XABORT(NAMSBR//
     >': This type of cell cannot be tracked by NXT:')
                        ENDIF
*----
*  Define position of cyclic starting point
*----
                        DO IDIR=1,NDIM
                          AAA=DZERO
                          DO ISINT=NBCOR(1),NBSINT-NBCOR(2)+1
                            AAA=AAA+DCINT(ISINT)*ANGLES(IDIR)
                          ENDDO
                          TRKORN(IDIR)=TRKORI(IDIR)+AAA
                        ENDDO
                        ISDIRX=ABS(ICINT(0,NBSINT))
                        IF(ISDIRX .EQ. JSDIR) THEN
                          DO IDIR=1,NDIM
                            ANGLED(IDIR)=ANGLES(IDIR)
                          ENDDO
                          ANGLED(JSDIR)=-ANGLED(JSDIR)
                          GO TO 125
                        ENDIF
                        IF(IPER(ISDIRX) .EQ. 1) THEN
                          ANGLEN(ISDIRX)=-ANGLEN(ISDIRX)
                        ENDIF
                      ENDDO
 125                  CONTINUE
                      IFREF=-1
                    ELSE
                      DO IDIR=1,NDIM
                        ANGLES(IDIR)=ANGLED(IDIR)
                        ANGLEN(IDIR)=ANGLED(IDIR)
                        TRKOR2(IDIR)=TRKOR2(IDIR)+DEPART(IDIR,2,IANGL)
                        TRKORD(IDIR)=TRKOR2(IDIR)
                        TRKORI(IDIR)=TRKORD(IDIR)
                        TRKORN(IDIR)=TRKORD(IDIR)
                      ENDDO
                      IFREF=1
                    ENDIF
                  ENDIF
                ENDIF
                DO IDIR=1,NDIM
                  ANGLES(IDIR)=ANGLED(IDIR)
                  ANGLEN(IDIR)=ANGLED(IDIR)
                  TRKORD(IDIR)=TRKORN(IDIR)
                  TRKORI(IDIR)=TRKORD(IDIR)
                  TRKORN(IDIR)=TRKORD(IDIR)
                ENDDO
              ENDIF
              ISTART=0
            ENDIF
*----
*  End of tracking for this direction
*----
 105        CONTINUE
          ENDDO
*----
*  End of tracking for this set of quadrant
*----
        ENDDO
        DEALLOCATE(KANGL,TORIG)
*----
*  End of tracking for this set of angles
*----
      ENDDO
*----
*  Compute DVNOR and DSNOR by comparing ratio of analytical
*  to numerically integrated volume or surfaces.
*----
      NBVERR=0
      DSVERR=DZERO
      DMVERR=DZERO
      DAVERR=DZERO
      NBV0=0
      NBV1=0
      DO IREG=1,NFREG
        IF(IPRINT .GE. 20 ) THEN
          WRITE(IOUT,6030) IREG,SURVOL(IREG),DVNOR(IREG,1)
        ENDIF
        DO IDIR=1,NBDR
          IF(DVNOR(IREG,IDIR) .EQ. DZERO) THEN
            IF(IPRINT .GE. 10) THEN
              IANGL=NBDR-1
              WRITE(IOUT,9000) NAMSBR,IREG,IDIR
            ENDIF
            DVNOR(IREG,IDIR)=DONE
            IF(IDIR .EQ. 1) THEN
              NBV0=NBV0+1
            ELSE
              NBV1=NBV1+1
            ENDIF
          ELSE
            DVNOR(IREG,IDIR)=SURVOL(IREG)/DVNOR(IREG,IDIR)
            IF(IDIR .EQ. 1) THEN
              NBVERR=NBVERR+1
            ENDIF
          ENDIF
        ENDDO
        DCERR=100.0D0*(DONE-DVNOR(IREG,1))
        DMVERR=MAX(DMVERR,ABS(DCERR))
        DSVERR=DSVERR+DCERR*DCERR
        DAVERR=DAVERR+DCERR
      ENDDO
      IF(NBV0 .GT. 0) THEN
        WRITE(IOUT,9002) NAMSBR
      ENDIF
      IF(NBV1 .GT. 0) THEN
        WRITE(IOUT,9005) NAMSBR
      ENDIF
      DSVERR=SQRT(DSVERR/DBLE(NBVERR))
      DAVERR=DAVERR/DBLE(NBVERR)
      NBSERR=0
      DSSERR=DZERO
      DMSERR=DZERO
      DASERR=DZERO
      NBS0=0
      DO ISUR=1,NFSUR
        IF(IPRINT .GE. 20 ) THEN
          WRITE(IOUT,6031) ISUR,
     >    SURVOL(-ISUR),DSNOR(ISUR)
        ENDIF
        IF(DSNOR(ISUR) .EQ. DZERO) THEN
          IF(IPRINT .GE. 10) THEN
            WRITE(IOUT,9001) NAMSBR,-ISUR
          ENDIF
          NBS0=NBS0+1
          DSNOR(ISUR)=DONE
        ELSE
          DSNOR(ISUR)=SURVOL(-ISUR)
     >                /DSNOR(ISUR)
          NBSERR=NBSERR+1
        ENDIF
        DCERR=100.0D0*(DONE-DSNOR(ISUR))
        DMSERR=MAX(DMSERR,ABS(DCERR))
        DSSERR=DSSERR+DCERR*DCERR
        DASERR=DASERR+DCERR
      ENDDO
      IF(NBS0 .GT. 0) THEN
        WRITE(IOUT,9003) NAMSBR,NBS0
      ENDIF
      DSSERR=SQRT(DSSERR/DBLE(NBSERR))
      DASERR=DASERR/DBLE(NBSERR)
*----
*  Processing finished:
*  print track normalization vector.
*  and routine closing output header if required
*  and return
*----
      IF(ABS(IPRINT) .GE. 1) THEN
        WRITE(IOUT,6005) DSVERR,DMVERR,DAVERR
        IF(IPRINT .GE. 10) THEN
          DO IREG=1,NFREG
            WRITE(IOUT,6010) IREG,SURVOL(IREG)
            WRITE(IOUT,6012) DVNOR(IREG,1),
     >      100.0D0*(DONE-DVNOR(IREG,1))
          ENDDO
        ENDIF
        WRITE(IOUT,6006) DSSERR,DMSERR,DASERR
        IF(IPRINT .GE. 10) THEN
          DO ISUR=1,NFSUR
            WRITE(IOUT,6011) -ISUR,SURVOL(-ISUR)
            WRITE(IOUT,6012) DSNOR(ISUR),
     >      100.0D0*(DONE-DSNOR(ISUR))
          ENDDO
        ENDIF
        WRITE(IOUT,6001) NAMSBR
      ENDIF
*----
*  Save track normalisation vector
*----
      IF(MAXSGL .EQ. 0) THEN
        WRITE(IOUT,9004) NAMSBR
        MAXSGL=LINMAX
      ENDIF
*----
*  Deallocate temporary records
*----
      DEALLOCATE(DCINT,ICINT,DLENGT,LENGTH,NUMERO)
      RETURN
*----
*  Output formats
*----
 6000 FORMAT('(* Output from --',A6,'-- follows ')
 6001 FORMAT('   Output from --',A6,'-- completed *)')
 6005 FORMAT(' Global RMS, maximum and average errors (%) ',
     >       'on region volumes :',3(2X,F10.5))
 6006 FORMAT(' Global RMS, maximum and average errors (%) ',
     >       'on surface areas  :',3(2X,F10.5))
 6010 FORMAT(' Normalisation factors and relative errors (%) ',
     >       'for region ',I8,' with volume ',F19.10)
 6011 FORMAT(' Normalisation factors and relative error (%) ',
     >       'for surface ',I8,' with area   ',F19.10)
 6012 FORMAT((2X,F15.10,2X,F10.5))
 6026 FORMAT('Line',I10.10,'={',I10,',',I10,'};')
 6027 FORMAT('Problem with Line',I10.10,'={',5(I10,','),
     >       I10,',',F18.10,',',I10,',',F18.10,'};')
 6030 FORMAT(' Normalization volumes  =',I10,1P,2D20.10)
 6031 FORMAT(' Normalization surfaces =',I10,1P,2D20.10)
 6050 FORMAT(' Direction = ',I5,2X,'Scan =',4I5,2x,'Weight =',F20.15/
     >' Directions and starting point =',6F20.15)
 6051 FORMAT('% Torig',I10.10,/'Torig=[')
 6052 FORMAT(9F15.8)
 6061 FORMAT(I10)
 6062 FORMAT(F18.10)
 6063 FORMAT('% ',I10.10,/A7,'=[')
 6064 FORMAT(18X,'];')
 6065 FORMAT(
     >'xcol=jet(    2);'/
     >'lls=length(Torig);'/'lli=length(LinePos);'/
     >'is=-1;'/'nums=0;'/
     >'for i=1:lli'/'  if LineReg(i) < 0'/'    if is == -1'/
     >'      nums=nums+1;'/'      is=0;'/
     >'      xi=Torig(nums,1);'/'      yi=Torig(nums,2);'/
     >'      tx=Torig(nums,4);'/'      ty=Torig(nums,5);'/
     >'    elseif is == 0'/'      is=-1;'/'    end'/'  else'/
     >'    xf=xi+tx*LinePos(i);'/'    yf=yi+ty*LinePos(i);'/
     >'    xxx=line([xi,xf],[yi,yf]);'/
     >'    xi=xf;'/'    yi=yf;'/'    reg=LineReg(i);'/
     >'    set(xxx,''Color'',[xcol(reg,:)]);'/
     >'  end'/'end'/'clear Torig LinePosLineReg'/'pause ;')
 6070 FORMAT(1X,A8,5X,I10)
 6071 FORMAT(5F18.10)
 9000 FORMAT(1X,'***** Warning in ',A6,'*****'/
     >       7X,'For region ',I8,
     >       1X,'no crossing by angle ',I8)
 9001 FORMAT(1X,'***** Warning in ',A6,'*****'/
     >       7X,'For surface ',I8,
     >       1X,'no crossing by any angle ')
 9002 FORMAT(1X,'***** Warning in ',A6,'*****'/
     >       7X,' regions not tracked for any direction  ')
 9003 FORMAT(1X,' ***** Warning in ',A6,'*****'/
     >       7X,I8,' surfaces not tracked for direction ',I8)
 9004 FORMAT(1X,' ***** Warning in ',A6,'*****'/
     >       7X,' no line segments detected in tracking ')
 9005 FORMAT(1X,'***** Warning in ',A6,'*****'/
     >       7X,' regions not tracked for some directions')
 9100 FORMAT(1X,' ***** Error in ',A6,'***** for line ',2I10/
     >       3F20.15)
 9101 FORMAT(1X,' ***** Error in ',A6,'***** for line ',I8/
     >       7X,'Positions (current and reference ) =',1P,2D21.14/
     >       7X,'Relative error = ',D21.14)
 9102 FORMAT(1X,' ***** Error in ',A6,'***** for line ',I8/
     >       7X,'Positions (current and reference ) =',1P,2D20.12/
     >       7X,'Absolute error = ',D21.14)
      END