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
|
*DECK SNQU02
SUBROUTINE SNQU02(NLF,JOP,U,W,TPQ,UPQ,VPQ,WPQ)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Set the level-symmetric (type 2) quadratures.
*
*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): A. Hebert
*
*Parameters: input
* NLF order of the SN approximation (even number).
*
*Parameters: output
* JOP number of base points per axial level in one octant.
* U base points in $\\xi$ of the axial quadrature. Used with
* zero-weight points.
* W weights for the axial quadrature in $\\xi$.
* TPQ base points in $\\xi$ of the 2D SN quadrature.
* UPQ base points in $\\mu$ of the 2D SN quadrature.
* VPQ base points in $\\eta$ of the 2D SN quadrature.
* WPQ weights of the 2D SN quadrature.
*
*-----------------------------------------------------------------------
*
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NLF,JOP(NLF/2)
REAL U(NLF/2),W(NLF/2),TPQ(NLF*(NLF/2+1)/4),UPQ(NLF*(NLF/2+1)/4),
1 VPQ(NLF*(NLF/2+1)/4),WPQ(NLF*(NLF/2+1)/4)
*----
* LOCAL VARIABLES
*----
PARAMETER(PI=3.141592654,MAXNLF=24,MAXEQ=64,MAXNBA=78,MAXW=16)
INTEGER INWEI(MAXNBA)
DOUBLE PRECISION WSUM2,WEI(MAXW),ZMAT(MAXEQ,MAXW+1),UD(MAXW)
*----
* SET THE UNIQUE QUADRATURE VALUES.
*----
IF(NLF.GT.MAXNLF) CALL XABORT('SNQU02: MAXNLF OVERFLOW.')
M2=NLF/2
NPQ=M2*(M2+1)/2
ZMU1=1.0D0/(3.0D0*DBLE(NLF-1))
NW=1+(NLF*(NLF+8)-1)/48
IF(NW.GT.MAXW) CALL XABORT('SNQU02: MAXW OVERFLOW.')
IF(NLF.EQ.2) THEN
ZMU1=0.33333333
ELSE IF(NLF.EQ.4) THEN
ZMU1=0.12251480
ELSE IF(NLF.EQ.6) THEN
ZMU1=0.07109447
ELSE IF(NLF.EQ.8) THEN
ZMU1=0.04761903
ELSE IF(NLF.EQ.10) THEN
ZMU1=0.03584310
ELSE IF(NLF.EQ.12) THEN
ZMU1=0.02796615
ELSE IF(NLF.EQ.14) THEN
ZMU1=0.02310250
ELSE IF(NLF.EQ.16) THEN
ZMU1=0.01931398
ELSE IF(NLF.EQ.18) THEN
ZMU1=0.01692067
ELSE IF(NLF.EQ.20) THEN
ZMU1=0.01455253
ELSE
CALL XABORT('SNQU02: 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(NLF-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('SNQU02: 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('SNQU02: INVALID VALUE OF NW.')
IF(IPR.NE.NPQ) CALL XABORT('SNQU02: BAD VALUE ON NPQ.')
*----
* SET THE RECTANGULAR SYSTEM AND SOLVE IT USING THE QR METHOD.
*----
NEQ=0
DO IPL=0,NLF,2
DO IPK=IPL,NLF-IPL,2
IF(MOD(IPL+IPK,2).EQ.1) CYCLE
NEQ=NEQ+1
IF(NEQ.GT.MAXEQ) CALL XABORT('SNQU02: 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.
*----
IPQ=0
WSUM=0.0
DO IP=1,M2
WSUM2=0.0D0
DO IQ=1,JOP(IP)
IPQ=IPQ+1
WPQ(IPQ)=REAL(WEI(INWEI(IPQ))*PI/2.0)
WSUM2=WSUM2+WEI(INWEI(IPQ))
ENDDO
W(IP)=REAL(WSUM2)
WSUM=WSUM+REAL(WSUM2*PI/2.0)
ENDDO
RETURN
END
|
