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
|
*DECK ALCACT
SUBROUTINE ALCACT(LCACT,NG,XG,WG)
*
*-----------------------------------------------------------------------
*
*Purpose:
* set the weights and base points for the different polar quadrature:
* - "Cactus" (Halsall)
* - "optimized" (Leonard and extension by Le Tellier)
*
*Copyright:
* Copyright (C) 1999 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): R. Roy and R. Le Tellier
*
*Reference:
* A. Leonard and C.T. Mc Daniel, "Optimal polar angles and weights for
* the characteristic method", Trans. Am. Nucl. Soc., 73, 172 (1995).
*
*Parameters: input
* LCACT type of quadrature (=1,2: values used by Halsall in
* Cactus (1980); =3: values proposed by Mc Daniel ng=2 only,
* extended to 3 and 4 with the same approach by R. Le Tellier
* in 06/2005)
* NG number of weights/base points.
*
*Parameters: output
* XG base points.
* WG weights.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER LCACT,NG
REAL XG(NG),WG(NG)
*----
* LOCAL VARIABLES
*----
INTEGER IG
DOUBLE PRECISION PI, PHIM, PHIP, DELPHI, WTEST, XTEST, YTEST,
> ZTEST, ZERO, ONE, HALF, QUART
PARAMETER( PI=3.14159265358979D0, ZERO=0.D0, ONE=1.D0,
> HALF=0.5D0, QUART=0.25D0 )
* For "Optimized" Quadrature
INTEGER I, I2, I3, I4, I5, IDEP, IFIN
PARAMETER ( I2= 1, I3= 3, I4= 6, I5= 10)
REAL XOP0(I5-1), WOP0(I5-1) ! Quadrature which minimizes the error on Ki2 without constraints.
REAL XOP1(I5-1), WOP1(I5-1) ! Quadrature which minimizes the error on Ki2 with P1 constraints.
REAL XGOP(I5-1), WGOP(I5-1) ! Gauss optimized quadrature.
SAVE XOP0, WOP0, XOP1, WOP1, XGOP, WGOP
* NG=2
DATA (XOP0(I),I=I2,I3-1) / 0.273658, 0.865714/
DATA (WOP0(I),I=I2,I3-1) / 0.139473, 0.860527/
DATA (XOP1(I),I=I2,I3-1) / 0.340183, 0.894215/
DATA (WOP1(I),I=I2,I3-1) / 0.194406, 0.805594/
DATA (XGOP(I),I=I2,I3-1) / 0.399374, 0.914448/
DATA (WGOP(I),I=I2,I3-1) / 0.250547, 0.749453/
* NG=3
DATA (XOP0(I),I=I3,I4-1) / 0.891439, 0.395534, 0.099812/
DATA (WOP0(I),I=I3,I4-1) / 0.793820, 0.188560, 0.017620/
DATA (XOP1(I),I=I3,I4-1) / 0.131209, 0.478170, 0.920079/
DATA (WOP1(I),I=I3,I4-1) / 0.029991, 0.250860, 0.719149/
DATA (XGOP(I),I=I3,I4-1) / 0.231156, 0.639973, 0.954497/
DATA (WGOP(I),I=I3,I4-1) / 0.085302, 0.341456, 0.573242/
* NG=4
DATA (XOP0(I),I=I4,I5-1) / 0.464167, 0.908274, 0.166004, 0.042181/
DATA (WOP0(I),I=I4,I5-1) / 0.218331, 0.746430, 0.032141, 0.003098/
DATA (XOP1(I),I=I4,I5-1) / 0.054819, 0.212313, 0.546065, 0.932318/
DATA (WOP1(I),I=I4,I5-1) / 0.005225, 0.051270, 0.272789, 0.670716/
DATA (XGOP(I),I=I4,I5-1) / 0.152641, 0.450820, 0.769181, 0.972320/
DATA (WGOP(I),I=I4,I5-1) / 0.037508, 0.167623, 0.338496, 0.456373/
! DATA (XGOP(I),I=I5,I6-1) / 0.159153, 0.450941, 0.724750, 0.868405,
! 1 0.980414/
! DATA (WGOP(I),I=I5,I6-1) / 0.039977, 0.156595, 0.220669, 0.204037,
! 1 0.378722/
*
IF( LCACT.EQ.1 )THEN
*---
* CACTUS 1:
*---
* Equal weight quadrature
DELPHI= ONE/DBLE(NG)
PHIM = ZERO
XTEST = ZERO
DO 10 IG= 1, NG
PHIP = DBLE(IG) * DELPHI
WTEST = PHIP - PHIM
WG(IG)= REAL( WTEST )
IF(IG.EQ.NG) THEN
YTEST = PI/2.0
ELSE
YTEST = SQRT(ONE - PHIP * PHIP) * PHIP + ASIN(PHIP)
ENDIF
ZTEST = HALF * (YTEST - XTEST) / WTEST
XG(IG)= REAL( SQRT(ONE - ZTEST*ZTEST) )
PHIM = PHIP
XTEST = YTEST
10 CONTINUE
ELSEIF( LCACT.EQ.2 )THEN
*---
* CACTUS 2:
*---
* Uniformly distributed angle quadrature
DELPHI= PI/DBLE(2*NG)
PHIM = ZERO
DO 20 IG= 1, NG
PHIP = DBLE(IG) * DELPHI
WTEST = SIN(PHIP) - SIN(PHIM)
WG(IG)= REAL( WTEST )
XTEST = HALF * (PHIP - PHIM)
YTEST = QUART * (SIN(PHIP+PHIP) - SIN(PHIM+PHIM))
ZTEST = (XTEST + YTEST) / WTEST
XG(IG)= REAL( SQRT(ONE - ZTEST*ZTEST) )
PHIM= PHIP
20 CONTINUE
ELSEIF(( LCACT.GE.3 ).AND.( LCACT.LE.5 ))THEN
*---
* OPTIMIZED:
*---
IDEP=0
IFIN=0
IF( NG.EQ.2 ) THEN
IDEP=I2
IFIN=I3-1
ELSE IF( NG.EQ.3 ) THEN
IDEP=I3
IFIN=I4-1
ELSE IF( NG.EQ.4 ) THEN
IDEP=I4
IFIN=I5-1
ELSE
CALL XABORT('ALCACT: LCACA=3 => NG= 2, 3 OR 4')
ENDIF
IF (LCACT.EQ.3) THEN
* Quadrature which minimizes the error on Ki2 without constraints.
DO 30 IG= IDEP, IFIN
XG(IG-IDEP+1)= REAL(SQRT(ONE - XOP0(IG)*XOP0(IG)))
WG(IG-IDEP+1)= WOP0(IG)
30 CONTINUE
ELSEIF (LCACT.EQ.4) THEN
* Quadrature which minimizes the error on Ki2 with P1 constraints.
DO 40 IG= IDEP, IFIN
XG(IG-IDEP+1)= REAL(SQRT(ONE - XOP1(IG)*XOP1(IG)))
WG(IG-IDEP+1)= WOP1(IG)
40 CONTINUE
ELSEIF (LCACT.EQ.5) THEN
* Gauss optimized quadrature.
DO 50 IG= IDEP, IFIN
XG(IG-IDEP+1)= REAL(SQRT(ONE - XGOP(IG)*XGOP(IG)))
WG(IG-IDEP+1)= WGOP(IG)
50 CONTINUE
ENDIF
ELSE
CALL XABORT('ALCACT: *LCACT* IN [1,5]')
ENDIF
RETURN
END
|
