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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
|
*DECK NXTLSN
SUBROUTINE NXTLSN(NDIM ,ORDRE ,NQUAD ,NBANGL,DQUAD ,DANGLT,
> DDENWT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* To define level-symmetric (type 2) quadrature angles.
*
*
*Copyright:
* Copyright (C) 2006 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):
* A. Hebert and R. Le Tellier
*
*Parameters: input
* NDIM number of dimensions for geometry.
* ORDRE 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 director cosines of angles.
* DDENWT angular density for each angle.
*
*----------
*
IMPLICIT NONE
*----
* Subroutine arguments
*----
INTEGER NDIM,ORDRE,NQUAD,NBANGL
DOUBLE PRECISION DQUAD(NQUAD),DANGLT(NDIM,NQUAD,NBANGL),
> DDENWT(NQUAD,NBANGL)
*----
* Local parameters
*----
INTEGER IOUT
CHARACTER NAMSBR*6
PARAMETER (IOUT=6,NAMSBR='NXTLSN')
DOUBLE PRECISION DZERO,DONE,DTWO
PARAMETER (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
INTEGER MAXORD,MAXNBA,MAXEQ,MAXW
PARAMETER (MAXORD=20,MAXNBA=55,MAXEQ=64,MAXW=12)
*----
* Functions
*----
DOUBLE PRECISION XDRCST,PI
*----
* Local variables
*----
INTEGER I,IP,IQ,IS,IW,IANG,JOP(MAXORD/2),M2,NPQ,NW,IPR,
> IPL,IPK,INMAX,IPQ,NW0,II,KK,LL,NEQ
REAL U(MAXORD/2),TPQ(MAXNBA),UPQ(MAXNBA),
> VPQ(MAXNBA),WPQ(MAXNBA)
DOUBLE PRECISION FAC,ZMU,ZETA,ZMU1,ZMU2,DDA,X,Y,Z,REF
INTEGER INWEI(MAXNBA)
DOUBLE PRECISION WEI(MAXW),ZMAT(MAXEQ,MAXW+1),UD(MAXW)
*----
* Set the unique quadrature values
*----
IF(ORDRE.GT.MAXORD) CALL XABORT(NAMSBR//': MAXORD OVERFLOW.')
M2=ORDRE/2
NPQ=M2*(M2+1)/2
ZMU1=1.0D0/(3.0D0*DBLE(ORDRE-1))
NW=1+(ORDRE*(ORDRE+8)-1)/48
IF(NW.GT.MAXW) CALL XABORT('NXTLSN: MAXW OVERFLOW.')
IF(ORDRE.EQ.2) THEN
ZMU1=0.33333333
ELSE IF(ORDRE.EQ.4) THEN
ZMU1=0.12251480
ELSE IF(ORDRE.EQ.6) THEN
ZMU1=0.07109447
ELSE IF(ORDRE.EQ.8) THEN
ZMU1=0.04761903
ELSE IF(ORDRE.EQ.10) THEN
ZMU1=0.03584310
ELSE IF(ORDRE.EQ.12) THEN
ZMU1=0.02796615
ELSE IF(ORDRE.EQ.14) THEN
ZMU1=0.02310250
ELSE IF(ORDRE.EQ.16) THEN
ZMU1=0.01931398
ELSE IF(ORDRE.EQ.18) THEN
ZMU1=0.01692067
ELSE IF(ORDRE.EQ.20) THEN
ZMU1=0.01455253
ELSE
CALL XABORT(NAMSBR//': ORDER NOT AVAILABLE.')
ENDIF
U(1)=REAL(SQRT(ZMU1))
DO I=2,M2
ZMU2=ZMU1+2.0D0*DBLE(I-1)*(1.0D0-3.0D0*ZMU1)/DBLE(ORDRE-2)
U(I)=REAL(SQRT(ZMU2))
ENDDO
*----
* Compute the position of weights
*----
IPR=0
INMAX=0
DO IP=1,M2
JOP(IP)=M2-IP+1
DO IQ=1,JOP(IP)
IPR=IPR+1
IF(IPR.GT.MAXNBA) CALL XABORT('NXTLSN: MAXNBA OVERFLOW.')
TPQ(IPR)=U(IP)
UPQ(IPR)=U(M2+2-IP-IQ)
VPQ(IPR)=U(IQ)
IS=MIN(IP,IQ,M2+2-IP-IQ)
NW0=0
DO II=1,IS-1
NW0=NW0+(M2-3*(II-1)+1)/2
ENDDO
KK=IP-IS+1
LL=IQ-IS+1
IF(KK.EQ.1)THEN
INWEI(IPR)=NW0+MIN(LL,M2-3*(IS-1)+1-LL)
ELSEIF(LL.EQ.1)THEN
INWEI(IPR)=NW0+MIN(KK,M2-3*(IS-1)+1-KK)
ELSE
INWEI(IPR)=NW0+MIN(KK,LL)
ENDIF
INMAX=MAX(INMAX,INWEI(IPR))
ENDDO
ENDDO
IF(INMAX.NE.NW) CALL XABORT(NAMSBR//': INVALID VALUE OD NW.')
IF(IPR.NE.NPQ) CALL XABORT(NAMSBR//': BAD VALUE ON NPQ.')
*----
* Set the rectangular system and solve it using the QR method
*----
NEQ=0
DO IPL=0,ORDRE,2
DO IPK=IPL,ORDRE-IPL,2
IF(MOD(IPL+IPK,2).EQ.1) CYCLE
NEQ=NEQ+1
IF(NEQ.GT.MAXEQ) CALL XABORT(NAMSBR//': MAXEQ OVERFLOW.')
DO IW=1,NW
ZMAT(NEQ,IW)=0.0D0
ENDDO
DO IPQ=1,NPQ
ZMU=TPQ(IPQ)
ZETA=UPQ(IPQ)
IW=INWEI(IPQ)
ZMAT(NEQ,IW)=ZMAT(NEQ,IW)+(ZMU**IPK)*(ZETA**IPL)
ENDDO
REF=1.0D0/DBLE(IPK+IPL+1)
DO I=1,IPL-1,2
REF=REF*DBLE(I)/DBLE(IPK+I)
ENDDO
ZMAT(NEQ,NW+1)=REF
ENDDO
ENDDO
CALL ALST2F(MAXEQ,NEQ,NW,ZMAT,UD)
CALL ALST2S(MAXEQ,NEQ,NW,ZMAT,UD,ZMAT(1,NW+1),WEI)
*----
* Set the level-symmetric quadratures
*----
PI=XDRCST('Pi',' ')
IPQ=0
DO IP=1,M2
DO IQ=1,JOP(IP)
IPQ=IPQ+1
WPQ(IPQ)=REAL(WEI(INWEI(IPQ))*PI)/2.0
ENDDO
ENDDO
*----
* Fill-in DANGLT and DDENWT array
*----
FAC=4.0*PI
DO IANG=1,NBANGL
X = DBLE(TPQ(IANG))
Y = DBLE(UPQ(IANG))
Z = DBLE(VPQ(IANG))
DDA=DBLE(FAC/WPQ(IANG))
DANGLT(1,1,IANG)=X
DANGLT(2,1,IANG)=Y
DANGLT(3,1,IANG)=Z
DDENWT(1,IANG)=DQUAD(1)*DDA
DANGLT(1,2,IANG)=-X
DANGLT(2,2,IANG)=Y
DANGLT(3,2,IANG)=Z
DDENWT(2,IANG)=DQUAD(2)*DDA
DANGLT(1,3,IANG)=X
DANGLT(2,3,IANG)=-Y
DANGLT(3,3,IANG)=Z
DDENWT(3,IANG)=DQUAD(3)*DDA
DANGLT(1,4,IANG)=-X
DANGLT(2,4,IANG)=-Y
DANGLT(3,4,IANG)=Z
DDENWT(4,IANG)=DQUAD(4)*DDA
ENDDO
*
RETURN
END
|
