diff options
Diffstat (limited to 'Dragon/src/NXTQSS.f')
| -rw-r--r-- | Dragon/src/NXTQSS.f | 173 |
1 files changed, 173 insertions, 0 deletions
diff --git a/Dragon/src/NXTQSS.f b/Dragon/src/NXTQSS.f new file mode 100644 index 0000000..6336ec9 --- /dev/null +++ b/Dragon/src/NXTQSS.f @@ -0,0 +1,173 @@ +*DECK NXTQSS + SUBROUTINE NXTQSS(IPRINT,NDIM ,ITYPBC,MAXMSH,NUCELL,DENUSR, + > DGMESH,NPLANE,NPOINT,DENLIN,SPACLN, + > WEIGHT,RADIUS,CENTER) +* +*----------------------------------------------------------------------- +* +*Purpose: +* To define standard spatial quadrature. +* +*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, R.Roy +* +*Parameters: input +* IPRINT print level. +* NDIM number of dimensions for geometry. +* ITYPBC type of boundary conditions ( +* +/- 2 for hexagones; +* +/- 1 for annular; +* 0 for Cartesian). +* MAXMSH maximum number of elements in mesh vector for +* each directions. +* NUCELL number of cell after unfolding in +* $X$, $Y$ and $Z$ directions. +* DENUSR user defined track density. +* DGMESH meshing vector for global geometry. +* +*Parameters: output +* NPLANE number of normal planes considered. +* NPOINT number of integration points along each axis +* in a plane mormal to track direction. +* DENLIN effective track density in plane. +* 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. +* +*Reference: +* G. Marleau, +* New Geometries Processing in DRAGON: The NXT: Module, +* Report IGE-260, Polytechnique Montreal, +* Montreal, 2005. +* \\\\ +* Extracted from the subroutine XELTI2 and XELTI3. +* +*---------- +* + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + INTEGER IPRINT,NDIM,ITYPBC,MAXMSH,NUCELL(3) + DOUBLE PRECISION DENUSR,DGMESH(-1:MAXMSH,3) + INTEGER NPLANE,NPOINT + DOUBLE PRECISION DENLIN,SPACLN,WEIGHT + DOUBLE PRECISION RADIUS,CENTER(NDIM) +*---- +* Local parameters +*---- + INTEGER IOUT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='NXTQSS') + DOUBLE PRECISION DZERO,DONE,DTWO + PARAMETER (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0) +*---- +* Local variables +*---- + INTEGER IDIR,LSTCEL,NPO2,IX + DOUBLE PRECISION DM,XD +*---- +* Processing starts: +* print routine openning output header if required +* and initialize various parameters. +*---- + IF(IPRINT .GE. 10) THEN + WRITE(IOUT,6000) NAMSBR + ENDIF +*---- +* Find the radius of the sphere (3-D) or circle surrounding +* the cell. +* Also find the true center of the cell +*---- + IF(NDIM .EQ. 3) THEN + NPLANE=3 + DENLIN=SQRT(DENUSR) + ELSE IF(NDIM .EQ. 2) THEN + NPLANE=1 + DENLIN=DENUSR + ENDIF + IF(ABS(ITYPBC) .EQ. 2) THEN +*---- +* HEXAGONAL +*---- + RADIUS=DZERO + DO IDIR=1,2 + LSTCEL=NUCELL(IDIR) + CENTER(IDIR)=DGMESH(1,IDIR) + DM=DZERO + DO IX=1,LSTCEL + XD=DGMESH(IX,IDIR)-DGMESH(1,IDIR)+DGMESH(0,IDIR) + DM=MAX(DM,XD*XD) + ENDDO + RADIUS=RADIUS+DM + ENDDO + DM=DZERO + DO IDIR=3,NDIM + LSTCEL=NUCELL(IDIR) + CENTER(IDIR)=(DGMESH(LSTCEL,IDIR)+DGMESH(0,IDIR))/DTWO + DM=((DGMESH(LSTCEL,IDIR)-DGMESH(0,IDIR))/DTWO)**2 + ENDDO + RADIUS=SQRT(RADIUS+DM) + ELSE IF(ABS(ITYPBC) .EQ. 1) THEN +*---- +* RADIAL +*---- + RADIUS=DZERO + DO IDIR=1,NDIM + LSTCEL=NUCELL(IDIR) + CENTER(IDIR)=(DGMESH(LSTCEL,IDIR)+DGMESH(0,IDIR))/DTWO + RADIUS=RADIUS+(DGMESH(LSTCEL,IDIR)-DGMESH(0,IDIR))**2 + ENDDO + RADIUS=SQRT(RADIUS)/DTWO + ELSE +*---- +* CARTESIAN +*---- + RADIUS=DZERO + DO IDIR=1,NDIM + LSTCEL=NUCELL(IDIR) + CENTER(IDIR)=(DGMESH(LSTCEL,IDIR)+DGMESH(0,IDIR))/DTWO + RADIUS=RADIUS+(DGMESH(LSTCEL,IDIR)-DGMESH(0,IDIR))**2 + ENDDO + RADIUS=SQRT(RADIUS)/DTWO + ENDIF + NPOINT=INT(DTWO*RADIUS*DENLIN)+1 + NPO2=NPOINT/2 + NPOINT=2*NPO2+1 + SPACLN=DTWO*RADIUS/DBLE(NPOINT) + DENLIN=DONE/SPACLN + IF(NDIM .EQ. 3) DENLIN=DENLIN*DENLIN + WEIGHT=DONE/(DBLE(NPLANE)*DENLIN) + IF(IPRINT .GE. 10) THEN + WRITE(IOUT,6002) DENUSR,DENLIN,WEIGHT,NPOINT,SPACLN + WRITE(IOUT,6010) RADIUS,(CENTER(IDIR),IDIR=1,NDIM) + WRITE(IOUT,6001) NAMSBR + ENDIF +*---- +* Processing finished: return +*---- + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + 6002 FORMAT(' Spatial tracking information :'/ + > ' Minimum tracking density requested =',1P,E15.7/ + > ' Effective tracking density selected=',E15.7/ + > ' Effective tracking weight =',E15.7/ + > ' Number of points per direction =',I15/ + > ' Linear track spacing =',E15.7) + 6010 FORMAT(' RADIUS = ',1P,E20.12/ + > ' CENTER = ',3E20.12) + END |