path: root/Dragon/src/NXTTLS.f

*DECK NXTTLS
      SUBROUTINE NXTTLS(IPTRK ,IFTEMP,IPRINT,IGTRK ,NDIM  ,MAXMSH,
     >                  NFSUR ,NFREG ,NUCELL,NBUCEL,NBANGL,NQUAD ,
     >                  NPLANE,NPOINT,LINMAX,MXGSUR,MXGREG,RENO  ,
     >                  MAXPIN,NBTDIR,NBDR  ,ITYPBC,IFMT  ,
     >                  RCUTOF,SPACLN,WEIGHT,RADIUS,CENTER,
     >                  IUNFLD,SURVOL,DGMESH,DANGLT,DDENWT,
     >                  MAXSGL,NTLINE,DVNOR ,DSNOR)
*
*-----------------------------------------------------------------------
*
*Purpose:
* To generate the standard tracking lines (isotropic 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 data structure 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.
* NQUAD   number of quadrant (in 3-D) and quarter (in 2-D).
* NPLANE  number of normal planes considered.
* NPOINT  number of integration points along each axis
*         in a plane mormal to track direction.
* 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.
* MAXPIN  maximum number of pins in a cell.
* NBTDIR  number of tracks directions considered.
* 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:.
* RCUTOF  corner cutoff.
* SPACLN  linear track spacing in the plane.
* WEIGHT  weight associated with each line in the plane.
* RADIUS  radius of circle (2-D) or sphere (3-D) surrounding
*         the geometry.
* CENTER  center of circle (2-D) or sphere (3-D) surrounding
*         the geometry.
* IUNFLD  description of unfolded geometry.
* SURVOL  global surface volume vector.
* DGMESH  meshing vector for global geometry.
* DANGLT  angles.
* DDENWT  angular density for each angle.
*
*Parameters: output
* 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),NQUAD,NBANGL,NBUCEL,
     >                 NPLANE,NPOINT,LINMAX,MXGSUR,MXGREG,RENO,
     >                 MAXPIN,NTLINE,NBTDIR,NBDR,ITYPBC,IFMT,MAXSGL
      DOUBLE PRECISION RCUTOF,SPACLN,WEIGHT
      DOUBLE PRECISION RADIUS,CENTER(NDIM)
      INTEGER          IUNFLD(2,NBUCEL)
      DOUBLE PRECISION SURVOL(-NFSUR:NFREG),DGMESH(-1:MAXMSH,5),
     >                 DANGLT(NDIM,NQUAD,NBANGL),DDENWT(NQUAD,NBANGL),
     >                 DVNOR(NFREG,NBDR),
     >                 DSNOR(NFSUR)
*----
*  Local parameters
*----
      INTEGER          IOUT
      CHARACTER        NAMSBR*6
      PARAMETER       (IOUT=6,NAMSBR='NXTTLS')
      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          NPO2,NCUTOF,
     >                 ITDIR,IANGL,IQUAD,IPLANE,IDIR,
     >                 NPTA2,NPTA3,IPTA2,IPTA3,ITST,NBCOR(2),NBSINT,
     >                 ITRN,ICEL,ICI,JLINE,IBLIN,IELIN,NBSEG,ISURF,
     >                 IX,IY,IZ,IOX,IOY,IOZ,IOC,LMAXT
      DOUBLE PRECISION RAD2G,RAD2T,ANGLES(3,3),DNPDIR(3,2,3),
     >                 TRKORI(3),TRKOR2(3),TCUTOF(3,2,2),TRKLIM(2),
     >                 CELLPO(3,2),DSTART,DCERR,FACVOL,FACSUR,
     >                 DSVERR,DMVERR,DAVERR,DSSERR,DMSERR,DASERR,DERR
      DOUBLE PRECISION DWGT,DAWGT,VCONT,VCONTA,WGTFAC,TORIG(3)
      INTEGER          ISEG,JSEG,IREG,ILREG,ISUR,NBVERR,NBSERR,
     >                 NBV0,NBV1,NBS0,IPRINL
      INTEGER          ISD,NSDEB,ISF,NSFIN,NBREG,ISBL,NTSEG,IST
      INTEGER          IOFF,NGLINE
      DOUBLE PRECISION TOTVE,TOTVA,TOTVD,TOTVDR,DELV,
     >                 TOTSE,TOTSA,TOTSD,TOTSDR,DELS
      INTEGER          NICSS,ICISS,ICSR,ICSRR
*----
*  Allocatable arrays (local)
*----
      INTEGER, ALLOCATABLE, DIMENSION(:) :: NUMERO
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ICINT
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: LENGTH,DLENGT,
     >                                               DCINT,DWGTRK
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DDIRET
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: DORITR
*----
*  Data
*----
      CHARACTER        CDIR(4)*1
      SAVE             CDIR
      DATA             CDIR /'X','Y','Z','R'/
*----
*  Processing starts:
*  print routine openning output header if required
*  and initialize various parameters.
*----
      IF(IPRINT .GE. 1) THEN
        WRITE(IOUT,6000) NAMSBR
      ENDIF
      LMAXT=4*(NBUCEL+4)
      WGTFAC=DONE
      PI=XDRCST('Pi',' ')
      MAXSGL=0
      ITST=RENO
      ITST=1
      NPTA2=NPOINT
      NPTA3=1
      RAD2G=RADIUS*RADIUS
      IF(NDIM .EQ. 3) NPTA3=NPOINT
      DVNOR(:NFREG,:NBDR)=DZERO
      DSNOR(:NFSUR)=DZERO
      NPO2=NPOINT/2
      IF(RCUTOF .EQ. DZERO) THEN
        NCUTOF= 1
      ELSE
        NCUTOF= 4
      ENDIF
      IOC=0
*----
* 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.
*   DDIRET  direction of tracking lines.
*   DWGTRK  weight of tracking lines.
*   DORITR  origin of tracking lines.
*----
      ALLOCATE(NUMERO(LINMAX),ICINT(0:5,LMAXT))
      ALLOCATE(LENGTH(LINMAX),DLENGT(LINMAX),DCINT(LMAXT),
     >         DDIRET(NDIM,NQUAD*NBANGL),DWGTRK(NQUAD*NBANGL),
     >         DORITR(NDIM*(NDIM+1),NPLANE,NQUAD*NBANGL))
*----
*  Compute number of track directions
*----
      NBTDIR=0
      DO IANGL=1,NBANGL
        DO IQUAD=1,NQUAD
          IF(DDENWT(IQUAD,IANGL) .GT. DZERO) THEN
            NBTDIR=NBTDIR+1
          ENDIF
        ENDDO
      ENDDO
      FACVOL=DTWO
      IF(NDIM .EQ. 2) THEN
        FACSUR=PI
      ELSE
        FACSUR=DTWO*DTWO
      ENDIF
      IF(IPRINT .GT. 1000) THEN
        IF(ITYPBC .EQ. 2) THEN
            WRITE(IOUT,6070) 'DGMESH X',NUCELL(1)
            WRITE(IOUT,6071) (DGMESH(IX,1),IX=-1,NUCELL(1))
            WRITE(IOUT,6070) 'DGMESH Y',NUCELL(1)
            WRITE(IOUT,6071) (DGMESH(IX,2),IX=-1,NUCELL(1))
            IF(NDIM .EQ.3) THEN
              WRITE(IOUT,6070) 'DGMESH Z',NUCELL(3)
              WRITE(IOUT,6071) (DGMESH(IX,3),IX=-1,NUCELL(3))
            ENDIF
        ELSE
          DO IDIR=1,NDIM
            WRITE(IOUT,6070) 'DGMESH '//CDIR(IDIR),NUCELL(IDIR)
            WRITE(IOUT,6071) (DGMESH(IX,IDIR),IX=-1,NUCELL(IDIR))
          ENDDO
        ENDIF
      ENDIF
      DDIRET(:NDIM,:NQUAD*NBANGL)=DZERO
      DWGTRK(:NQUAD*NBANGL)=DZERO
      DORITR(:NDIM*(NDIM+1),:NPLANE,:NQUAD*NBANGL)=DZERO
*----
*  Loop over angles in a quarter (2-D) or a quadrant (3-D)
*----
      NGLINE=0
      NTLINE=0
      ITDIR=0
      DO IANGL=1,NBANGL
*----
*  Loop over 2 quarters (2-D) or 4 quadrants (3-D)
*----
        DO IQUAD=1,NQUAD
*----
*  Do not track angle with 0 density
*  because of the problem symmetry
*----
          IF(DDENWT(IQUAD,IANGL) .EQ. DZERO) GO TO 105
          DWGT=WEIGHT/DDENWT(IQUAD,IANGL)
          DAWGT=WEIGHT
*----
*  Track this angle
*----
          ITDIR=ITDIR+1
          DWGTRK(ITDIR)=DWGT
*----
*  Find planes mormal to selected direction
*----
          TRKORI(:3)=DZERO
          TRKOR2(:3)=DZERO
          ANGLES(:3,:3)=DZERO
          DNPDIR(:3,:2,:3)=DZERO
          DO IDIR=1,NDIM
            ANGLES(IDIR,1)=DANGLT(IDIR,IQUAD,IANGL)
            DDIRET(IDIR,ITDIR)=ANGLES(IDIR,1)
          ENDDO
          CALL NXTQPS(NDIM,ANGLES,DNPDIR)
*----
*  Loop over planes normal to direction
*  1 in 2-D and 3 in 3-D
*----
          DO IPLANE=1,NPLANE
            DO IDIR=1,NDIM
              DO JLINE=2,NDIM
                ANGLES(IDIR,JLINE)=SPACLN*DNPDIR(IDIR,JLINE-1,IPLANE)
              ENDDO
              IF(NCUTOF .NE. 1)THEN
                TCUTOF(IDIR,1,1)=RCUTOF*(ANGLES(IDIR,2)+ANGLES(IDIR,3))
                TCUTOF(IDIR,1,2)=RCUTOF*(ANGLES(IDIR,2)-ANGLES(IDIR,3))
                TCUTOF(IDIR,2,1)=-TCUTOF(IDIR,1,2)
                TCUTOF(IDIR,2,2)=-TCUTOF(IDIR,1,1)
              ENDIF
              TRKOR2(IDIR)=CENTER(IDIR)
     >                    -DBLE(NPO2+1)*(ANGLES(IDIR,2)+ANGLES(IDIR,3))
            ENDDO
*----
*  Fill array for localisation of integration lines
*----
            DO IDIR=1,NDIM
              DORITR(IDIR,IPLANE,ITDIR)=TRKOR2(IDIR)
            ENDDO
            IOFF=NDIM
            DO JLINE=1,NDIM
              DO IDIR=1,NDIM
                DORITR(IDIR+IOFF,IPLANE,ITDIR)=ANGLES(IDIR,JLINE)
              ENDDO
              IOFF=IOFF+NDIM
            ENDDO
*----
*  Loop over lines on first normal axis (direction 2)
*----
            DO IPTA2=1,NPTA2
*----
*  Displace starting point by an additional value of SPALLN
*  along first normal axis (direction 2)
*----
              DO IDIR=1,NDIM
                TRKOR2(IDIR)=TRKOR2(IDIR)+ANGLES(IDIR,2)
                TRKORI(IDIR)=TRKOR2(IDIR)
              ENDDO
*----
*  For 3-D models:
*  Loop over lines on second normal axis (direction 3)
*  For 2-D models:
*  NPTA3=1
*----
              DO IPTA3=1,NPTA3
                RAD2T=0.0
*----
*  For 3-D models:
*  Displace starting point by an additional value of SPALLN
*  along second normal axis (direction 3)
*  For 2-D models:
*  No displacement since ANGLES(IDIR,3)=0
*----
                NGLINE=NGLINE+1
                DO IDIR=1,NDIM
                  TRKORI(IDIR)=TRKORI(IDIR)+ANGLES(IDIR,3)
                  RAD2T=RAD2T+(TRKORI(IDIR)-CENTER(IDIR))**2
                ENDDO
*----
*  Eliminate tracks outside circle or sphere surrounding geometry
*----
                IF(RAD2T .GT. RAD2G) GO TO 115
*----
*  Find cells crossed by track
*----
                IPRINL=IPRINT
*                IF(NTLINE+1 .EQ. 676
*     >        .OR. NTLINE+1 .EQ. 1814
*     >        .OR. NTLINE+1 .EQ. 1831
*     >                           ) IPRINL=IPRINT+4000
                IRLA=-1
                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)
                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 115
*----
*  Test for multiple line segments in hexagonal assemblies
*
*----
                NICSS=1                  
                IF(ITYPBC .EQ. 2) THEN
                  DO ICI=2,NBSINT
                    IF(ICINT(0,ICI) .EQ. 0) THEN
                      NICSS=NICSS+1                  
                    ENDIF
                  ENDDO
                  ICSRR=NBCOR(1)+1
                ELSE
                  ICSRR=NBCOR(1)+1
                ENDIF
*----
*  For each region crossed loop track geometry
*  present in this region
*----
*                write(6,*) 'NICSS,ICSRR=',NICSS,ICSRR
                DO ICISS=1,NICSS
                ICSR=ICSRR
                NTLINE=NTLINE+1
                IF(IPRINL .GE. 500) THEN
                  WRITE(IOUT,6039)
                  WRITE(IOUT,6040) NTLINE,ITDIR,IPLANE,
     >            IPTA2,IPTA3,WEIGHT,SPACLN
                  WRITE(IOUT,6041) NTLINE,TRKORI
                  WRITE(IOUT,6042) NTLINE,ANGLES
                  WRITE(IOUT,6043) NTLINE,DCINT(ICSR-1)
                ELSE IF(IPRINL .GE. 50) THEN
                  WRITE(IOUT,6040) NTLINE,ITDIR,IPLANE,
     >            IPTA2,IPTA3,WEIGHT,SPACLN
                ENDIF
                IBLIN=1
                IELIN=0
                NUMERO(IBLIN)=0
                DLENGT(IBLIN)=DCINT(ICSR-1)
                DSTART=DCINT(ICSR-1)
                IBLIN=IBLIN+1
                DO ICI=ICSR,NBSINT-NBCOR(2)+1
                  ICSRR=ICSRR+1
                  IF(ITYPBC .EQ. 0) THEN
                    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
                  ELSE IF(ITYPBC .EQ. 2) THEN
                    IOC=ICINT(0,ICI+1)
*----
*  For multiple track segment IOC=0 indicates that the current
*  segment is completed and that a new track segment should be
*  started at ICI+2
*----
                    IF(IOC .EQ. 0) THEN
                      ICSRR=ICSRR+2
                      GO TO 125
                    ENDIF
                    IOC=ICINT(0,ICI)
                    IX=ICINT(1,ICI)
                    IY=ICINT(1,ICI)
                    CELLPO(1,2)=DGMESH(IX,1)
                    CELLPO(1,1)=DGMESH(IX,1)
                    CELLPO(2,2)=DGMESH(IY,2)
                    CELLPO(2,1)=DGMESH(IY,2)
                    IF(NDIM .EQ. 3) THEN
                      IZ=ICINT(3,ICI)
                      IOZ=(IZ-1)*NUCELL(1)
                      CELLPO(3,2)=DGMESH(IZ,3)
                      CELLPO(3,1)=DGMESH(IZ-1,3)
                      IOC=IOZ+IX                  
                    ENDIF
                  ENDIF
                  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)
                  DERR=MAX(ABS(TRKLIM(1)),ABS(TRKLIM(2)))
                  DERR=(DLENGT(IELIN)-DSTART)/DERR
                  IF(DERR .GT. DCUTOF) THEN
                    WRITE(IOUT,9100) NAMSBR,NTLINE,
     >              DLENGT(IELIN),DSTART,DERR
                    CALL XABORT(NAMSBR//
     >': End of track does not coincide with end of cell')
                  ENDIF
                  IBLIN=IELIN+1
                ENDDO
 125            CONTINUE               
                NBSEG=IELIN
*----
*  Compress tracking vector for region with DLENGT=0.0 and
*  for successive segments in the same region.
*----
                JSEG=0
                ILREG=-1
                NBREG=0
                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
                    ELSE
                      JSEG=JSEG+1
                      DLENGT(JSEG)=DLENGT(ISEG)
                      NUMERO(JSEG)=IREG
                      ILREG=-1
                    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
                NSFIN=0
                NBREG=0
                DO ISEG=1,NBSEG
                  IREG=NUMERO(ISEG)
                  IF(DLENGT(ISEG) .GT. DZERO) THEN
                    IF(IREG .GT. 0) THEN
                      NBREG=NBREG+1
                      VCONT=DLENGT(ISEG)*DWGT*FACVOL
                      DVNOR(IREG,1)=DVNOR(IREG,1)+VCONT
                      IF(NBDR .GT. 1) THEN
                        VCONTA=DLENGT(ISEG)*DAWGT
                        DVNOR(IREG,ITDIR+1)=DVNOR(IREG,ITDIR+1)+VCONTA
                      ENDIF
                      JSEG=JSEG+1
                      LENGTH(JSEG)=DLENGT(ISEG)
                      NUMERO(JSEG)=IREG
                    ELSE IF(IREG .LT. 0) THEN
                      IF(NBREG .EQ. 0) THEN
                        NSDEB=NSDEB+1
                      ELSE
                        NSFIN=NSFIN+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
                IF(NSDEB .GT. 1 .OR. NSFIN .GT. 1) THEN
                  IF(IPRINL .GE. 500) THEN
                    WRITE(IOUT,6025) NTLINE,ITDIR,NBSEG,
     >                               NSDEB,NSFIN,NBREG,
     >                               DLENGT(1),DLENGT(NBSEG)
                    WRITE(IOUT,6023) 'LineReg',NTLINE
                    WRITE(IOUT,6021) (NUMERO(ISEG),ISEG=1,NBSEG)
                    WRITE(IOUT,6024)
                    WRITE(IOUT,6023) 'LinePos',NTLINE
                    WRITE(IOUT,6022) (LENGTH(ISEG),ISEG=1,NBSEG)
                    WRITE(IOUT,6024)
                  ENDIF
                ENDIF
                MAXSGL=MAX(MAXSGL,NBSEG)
                IF(NSDEB*NSFIN .EQ. 0) THEN
*----
*  Missing outer or inner surface
*  Skip track and return warning
*----
                  WRITE(IOUT,9006) NAMSBR
                  WRITE(IOUT,9026) NTLINE,ITDIR,NBSEG,NSDEB,NSFIN,
     >                             DLENGT(1),DLENGT(NBSEG)
                ELSE
*----
*  Store line on temporary tracking file if required
*----
                  WGTFAC=DONE/DBLE(NSDEB*NSFIN)
                  IF(IPRINL .GE. 500) THEN
                    WRITE(IOUT,6020) NTLINE,ITDIR,NBSEG,DWGT*WGTFAC,
     >                               DLENGT(1),DLENGT(NBSEG)
*                    WRITE(6,*) IPLANE,IPTA2,IPTA3,
*     >                    (TRKORI(IST),IST=1,NDIM),DCINT(ICSR-1)
                    WRITE(IOUT,6023) 'LineReg',NTLINE
                    WRITE(IOUT,6021) (NUMERO(ISEG),ISEG=1,NBSEG)
                    WRITE(IOUT,6024)
                    WRITE(IOUT,6023) 'LinePos',NTLINE
                    WRITE(IOUT,6022) (LENGTH(ISEG),ISEG=1,NBSEG)
                    WRITE(IOUT,6024)
                  ENDIF
                  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. 500)
     >                WRITE(IOUT,6026) NTLINE-1+ISBL,NTLINE,ISBL
                      IF(IGTRK .EQ. 1) THEN
                        IF(IFMT .EQ. 1) THEN
                          DO IST=1,NDIM
                            TORIG(IST)=TRKORI(IST)+DCINT(ICSR-1)*
     >                      DANGLT(IST,IQUAD,IANGL)
                          ENDDO
                          WRITE(IFTEMP) 1,NTSEG,DWGT*WGTFAC,ITDIR,
     >                    NUMERO(ISD),
     >                    (NUMERO(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                    NUMERO(NBSEG-ISF),
     >                    DONE,
     >                    (LENGTH(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                    DONE,
     >                    NTLINE-1+ISBL,IPLANE,IPTA2,IPTA3,
     >                    (TORIG(IST),IST=1,NDIM)
                        ELSE
                          WRITE(IFTEMP) 1,NTSEG,DWGT*WGTFAC,ITDIR,
     >                    NUMERO(ISD),
     >                    (NUMERO(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                    NUMERO(NBSEG-ISF),
     >                    DONE,
     >                    (LENGTH(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >                    DONE
                        ENDIF
                      ENDIF
                    ENDDO
                  ENDDO
                ENDIF
                IF(ISBL .EQ. 0) THEN
                  WRITE(IOUT,6027) NTLINE,ITDIR,NBSEG,
     >                             NSDEB,NSFIN,NBREG,
     >                             NUMERO(1),DLENGT(1),
     >                             NUMERO(NBSEG),DLENGT(NBSEG)
                  ISBL=1
                  IF(IGTRK .EQ. 1) THEN
                    WRITE(IFTEMP) 1,NTSEG,DWGT*WGTFAC,ITDIR,
     >              NUMERO(1),
     >              (NUMERO(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >              NUMERO(NBSEG),
     >              DONE,
     >              (LENGTH(ISEG),ISEG=NSDEB+1,NBSEG-NSFIN),
     >              DONE
                  ENDIF
                ENDIF
                NTLINE=NTLINE-1+ISBL
                ENDDO
*----
*  Exit because line is outside circle or sphere surrounding geometry
*----
 115            CONTINUE
*----
*  END loop over points on second normal axis
*----
              ENDDO
*----
*  END loop over points on first normal axis
*----
            ENDDO
*----
*  END loop over planes
*----
          ENDDO
*----
*  Exit because angle with 0 density not tracked
*----
 105      CONTINUE
*----
*  END loop over quarter or quadrant
*----
        ENDDO
*----
*  END loop over angles
*----
      ENDDO
*----
*  Save general tracking information
*----
      CALL LCMPUT(IPTRK,'TrackingOrig',
     >            NDIM*(NDIM+1)*NPLANE*NBTDIR,4,DORITR)
      CALL LCMPUT(IPTRK,'TrackingWgtD',NBTDIR,4,DWGTRK)
      CALL LCMPUT(IPTRK,'TrackingDirc',NDIM*NBTDIR,4,DDIRET)
*----
*  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
      TOTVE=DZERO
      TOTVA=DZERO
      TOTVD=DZERO
      TOTVDR=DZERO
      DO IREG=1,NFREG
        IF(IPRINT .GE. 20 ) THEN
          WRITE(IOUT,6030) IREG,
     >    SURVOL(IREG),DVNOR(IREG,1)
        ENDIF
        TOTVE=TOTVE+SURVOL(IREG)
        TOTVA=TOTVA+DVNOR(IREG,1)
        DELV=SURVOL(IREG)-DVNOR(IREG,1)
        TOTVDR=TOTVDR+DELV*DELV
        DO IDIR=1,NBDR
          IF(DVNOR(IREG,IDIR) .EQ. DZERO) THEN
            IF(IPRINT .GE. 10) THEN
              WRITE(IOUT,9000) NAMSBR,IREG,ITDIR
            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
      TOTVD=100.0D0*(TOTVE-TOTVA)/TOTVE
      TOTVDR=100.0D0*SQRT(TOTVDR/DBLE(NBVERR))/TOTVE
      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
      TOTSE=DZERO
      TOTSA=DZERO
      TOTSD=DZERO
      TOTSDR=DZERO
      DO ISUR=1,NFSUR
        IF(IPRINT .GE. 20 ) THEN
          WRITE(IOUT,6031) ISUR,
     >    SURVOL(-ISUR),DSNOR(ISUR)
        ENDIF
        TOTSE=TOTSE+SURVOL(-ISUR)
        TOTSA=TOTSA+DSNOR(ISUR)
        DELS=SURVOL(-ISUR)-DSNOR(ISUR)
        TOTSDR=TOTSDR+DELS*DELS
        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
      TOTSD=100.0D0*(TOTSE-TOTSA)/TOTSE
      TOTSDR=100.0D0*SQRT(TOTSDR/DBLE(NBSERR))/TOTSE
      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(IPRINT .GE. 1) THEN
        WRITE(IOUT,6007) TOTVE,TOTVA,TOTVD,TOTVDR
        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,6008) TOTSE,TOTSA,TOTSD,TOTSDR
        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
*----
*  Scratch storage deallocation
*----
      DEALLOCATE(DORITR,DWGTRK,DDIRET,DCINT,DLENGT,LENGTH)
      DEALLOCATE(ICINT,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))
 6007 FORMAT(' Total exact volume       = ',F19.10/
     >       ' Total approximate volume = ',F19.10/
     >       ' Error on total volume    = ',4X,F10.5/
     >       ' RMS Error on total volume= ',4X,F10.5)
 6008 FORMAT(' Total exact surface area       = ',F19.10/
     >       ' Total approximate surface area = ',F19.10/
     >       ' Error on total surface area    = ',4X,F10.5/
     >       ' RMS Error on total surface area= ',4X,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))
 6020 FORMAT('Line',I10.10,'={',2(I10,','),
     >       F18.10,',',F18.10,',',F18.10,'};')
 6021 FORMAT(6(I10,:,',',9X))
 6022 FORMAT(6(F18.10,:,','))
 6023 FORMAT(A7,I10.10,'={')
 6024 FORMAT(18X,'};')
 6025 FORMAT('Line',I10.10,'={',5(I10,','),
     >       F18.10,',',F18.10,'};')
 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)
 6039 FORMAT(1X)
 6040 FORMAT('Track',I10.10,'={',4(I10,','),
     >F15.8,',',F15.8,'};')
 6041 FORMAT('Torig',I10.10,'={',2(F15.8,','),F15.8,'};')
 6042 FORMAT('Tdire',I10.10,'={',8(F15.8,','),F15.8,'};')
 6043 FORMAT('Tstrt',I10.10,'={',F15.8,'};')
 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')
 9006 FORMAT(1X,'***** Warning in ',A6,'*****'/
     >       7X,' Final or initial surface could not be identified')
 9026 FORMAT('Line',I10.10,'={',4(I10,','),
     >       F18.10,',',F18.10,'};')
 9100 FORMAT(1X,' ***** Error in ',A6,'***** for line ',I8/
     >       7X,'Positions (current and reference ) =',1P,2D21.14/
     >       7X,'Relative error = ',D21.14)
      END