1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
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
|
