From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/NXTTLS.f | 901 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 901 insertions(+) create mode 100644 Dragon/src/NXTTLS.f (limited to 'Dragon/src/NXTTLS.f') diff --git a/Dragon/src/NXTTLS.f b/Dragon/src/NXTTLS.f new file mode 100644 index 0000000..8613301 --- /dev/null +++ b/Dragon/src/NXTTLS.f @@ -0,0 +1,901 @@ +*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 -- cgit v1.2.3