diff options
Diffstat (limited to 'Dragon/src/NXTQAS.f')
| -rw-r--r-- | Dragon/src/NXTQAS.f | 150 |
1 files changed, 150 insertions, 0 deletions
diff --git a/Dragon/src/NXTQAS.f b/Dragon/src/NXTQAS.f new file mode 100644 index 0000000..7bebadf --- /dev/null +++ b/Dragon/src/NXTQAS.f @@ -0,0 +1,150 @@ +*DECK NXTQAS + SUBROUTINE NXTQAS(IPRINT,NDIM ,AZMQUA,NANGL ,NQUAD ,NBANGL, + > DQUAD ,DANGLT,DDENWT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* To define quadrature angles for a given tracking option. +* +*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, M. Hampartzounian +* +*Parameters: input +* IPRINT print level. +* NDIM number of dimensions for geometry. +* AZMQUA quadrature type. +* NANGL quadrature order. +* NQUAD number of quadrant (in 3-D) and quarter (in 2-D). +* NBANGL number of angles. +* DQUAD relative density of each quadrant. +* +*Parameters: output +* DANGLT angles. +* DDENWT angular density for each angle. +* +*Reference: +* G. Marleau, +* New Geometries Processing in DRAGON: The NXT: Module, +* Report IGE-260, Polytechnique Montreal, +* Montreal, 2005. +* Extracted from the subroutine XELTS2 of EXCELL. +* +*---------- +* + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + INTEGER IPRINT,NDIM,AZMQUA,NANGL,NQUAD,NBANGL + DOUBLE PRECISION DQUAD(NQUAD) + DOUBLE PRECISION DANGLT(NDIM,NQUAD,NBANGL),DDENWT(NQUAD,NBANGL) +*---- +* Local parameters +*---- + INTEGER IOUT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='NXTQAS') + DOUBLE PRECISION DZERO,DONE,DTWO + PARAMETER (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0) +*---- +* Functions +*---- + DOUBLE PRECISION XDRCST,PI +*---- +* Local variables +*---- + INTEGER IANG,IQUAD,IDIR + DOUBLE PRECISION DTHETA,THETA,DDA,COST,SINT +*---- +* Processing starts: +* print routine openning output header if required +* and initialize various parameters. +*---- + IF(IPRINT .GE. 10) THEN + WRITE(IOUT,6000) NAMSBR + WRITE(IOUT,6010) NDIM,AZMQUA,NANGL,NQUAD + ENDIF + PI=XDRCST('Pi',' ') + IF(NDIM .EQ. 2) THEN + IF(AZMQUA .EQ. 1) THEN +*---- +* Trapezoidal quadrature +* NBANGL point quadrature for both: +* (1/2*Pi)*Integral(0,Pi/2) and (1/2*Pi)*Integral(Pi/2,Pi) +* Quadrature weight = Pi/(2*NBANGL) +* DENSITY=(2*Pi)/weight=4*NBANGL +*---- + THETA =PI/DBLE(4*NANGL) + DTHETA=DTWO*THETA + DDA=DBLE(4*NBANGL) + DO IANG=1,NBANGL + COST=COS(THETA) + SINT=SIN(THETA) + DANGLT(1,1,IANG)=COS(THETA) + DANGLT(2,1,IANG)=SIN(THETA) + DDENWT(1,IANG)=DQUAD(1)*DDA + DANGLT(1,2,IANG)=-SIN(THETA) + DANGLT(2,2,IANG)=COS(THETA) + DDENWT(2,IANG)=DQUAD(2)*DDA + THETA=THETA+DTHETA + ENDDO + ELSE + WRITE(IOUT,9000) NAMSBR,AZMQUA + CALL XABORT(NAMSBR//': INVALID QUADRATURE OPTION IN 2D') + ENDIF + ELSE IF(NDIM .EQ. 3) THEN + IF(AZMQUA .EQ. 1) THEN + CALL NXTQEW(NDIM ,NANGL ,NQUAD ,NBANGL,DQUAD , + > DANGLT,DDENWT) + ELSE IF(AZMQUA .EQ. 4) THEN + CALL NXTQLC(NDIM ,NANGL ,NQUAD ,NBANGL,DQUAD , + > DANGLT,DDENWT) + ELSE IF(AZMQUA .EQ. 5) THEN + CALL NXTQLT(NDIM ,NANGL ,NQUAD ,NBANGL,DQUAD , + > DANGLT,DDENWT) + ELSE IF(AZMQUA .EQ. 6) THEN + CALL NXTLSN(NDIM ,NANGL ,NQUAD ,NBANGL,DQUAD , + > DANGLT,DDENWT) + ELSE IF(AZMQUA .EQ. 7) THEN + CALL NXTQRN(NDIM ,NANGL ,NQUAD ,NBANGL,DQUAD , + > DANGLT,DDENWT) + ENDIF + ENDIF +*---- +* Processing finished: return +*---- + IF(IPRINT .GE. 10) THEN + WRITE(IOUT,6011) + DO IANG=1,NBANGL + DO IQUAD=1,NQUAD + IF(DDENWT(IQUAD,IANG) .GT. DZERO) THEN + WRITE(IOUT,6012) IANG,IQUAD, + > (DANGLT(IDIR,IQUAD,IANG),IDIR=1,NDIM), + > DDENWT(IQUAD,IANG) + ENDIF + ENDDO + ENDDO + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + 6010 FORMAT(1X,'NDIM =',I8,1X,'AZMQUA=',I8, + > 1X,'NANGL =',I8,1X,'NQUAD =',I8) + 6011 FORMAT(' Tracking directions and weights '/ + > 1X,' Angle',1X,' Quadrant', + > 1X,' Directions and weight') + 6012 FORMAT(2(1X,I10),4(2X,F24.14)) + 9000 FORMAT(A6,': AZMQUA=',I5,' is invalid in 2D') + END |