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
|
*DECK SYBT1D
SUBROUTINE SYBT1D(NPIJ,RAD,LGSPH,NGAUSS,Z)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Tracking for pij calculation using the method of Kavenoky in annular
* or spherical geometry. The tracking is used by SYBALC.
*
*Copyright:
* Copyright (C) 2002 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
* NPIJ number of regions.
* RAD radius of regions array.
* LGSPH geometry flag (=.TRUE. for spherical geometry).
* NGAUSS number of Gauss points.
*
*Parameters: output
* Z tracking information.
*
*-----------------------------------------------------------------------
*
* REFERENCE:
* A. KAVENOKY, 'CALCUL ET UTILISATION DES PROBABILITES DE PREMIERE
* COLLISION POUR LES MILIEUX HETEROGENES A UNE DIMENSION: LES PROGRAMMES
* ALCOLL ET CORTINA', NOTE CEA-N-1077, COMMISSARIAT A L'ENERGIE
* ATOMIQUE, MARS 1969.
*----
* SUBROUTINE ARGUMENTS
*----
LOGICAL LGSPH
INTEGER NPIJ,NGAUSS
REAL RAD(0:NPIJ),Z(*)
*----
* LOCAL VARIABLES
*----
PARAMETER(MAXGAU=64,PI=3.1415926535)
REAL ALPR(MAXGAU),PWR(MAXGAU)
*
CALL ALGJP(NGAUSS,ALPR,PWR)
SUM=0.0
IZ=1
RIK1=RAD(0)
RI=1.0/RAD(NPIJ)
DO 60 IK=1,NPIJ
RIK2=RAD(IK)
RD=RIK2-RIK1
DO 50 IL=1,NGAUSS
R=RIK2-RD*ALPR(IL)**2
R2=R*R
IZ=IZ+2
IZ1=IZ
*
*----
* STORE INTEGRATION WEIGHT (TIMES 2) FOR THIS LINE
*----
AUX=RD*PWR(IL)*4.0
IF(LGSPH) AUX=AUX*R*PI
Z(IZ)=AUX
CT1=0.
CT2=0.
DO 40 I=IK,NPIJ
CT2=SQRT(RAD(I)**2-R2)
IZ=IZ+1
*----
* STORE LENGTH OF PATH
*----
Z(IZ)=CT2-CT1
CT1=CT2
40 CONTINUE
*----
* STORE COS(PHI)*INTEGRATION WEIGHT ( TIMES 2 )
*----
Z(IZ1-1)=Z(IZ1)*CT2*RI
SUM=SUM+AUX/CT2
50 CONTINUE
RIK1=RIK2
60 CONTINUE
IZ=IZ+1
IF(LGSPH) THEN
Z(1)=SUM/(PI*2.0*RAD(NPIJ))
ELSE
Z(1)=SUM/PI
ENDIF
RETURN
END
|
