summaryrefslogtreecommitdiff
path: root/Dragon/src/NXTTLC.f
diff options
context:
space:
mode:
authorstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
committerstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
commit7dfcc480ba1e19bd3232349fc733caef94034292 (patch)
tree03ee104eb8846d5cc1a981d267687a729185d3f3 /Dragon/src/NXTTLC.f
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Dragon/src/NXTTLC.f')
-rw-r--r--Dragon/src/NXTTLC.f1046
1 files changed, 1046 insertions, 0 deletions
diff --git a/Dragon/src/NXTTLC.f b/Dragon/src/NXTTLC.f
new file mode 100644
index 0000000..51dc74c
--- /dev/null
+++ b/Dragon/src/NXTTLC.f
@@ -0,0 +1,1046 @@
+*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
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:19:39 +0000