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
|
*DECK QIJI3D
SUBROUTINE QIJI3D(NREG,NSOUT,NPIJ,NGRP,MXSEG,NCOR,SWVOID,LINE,
> WEIGHT,NUMERO,LENGHT,SIGTAL,NPSYS,DPR)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Integration for general 3D isotropic tracking.
*
*Copyright:
* Copyright (C) 1991 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
*
*Parameters: input
* NREG total number of regions.
* NSOUT number of outer surface.
* NPIJ number of probabilities in one group.
* NGRP number of groups.
* MXSEG number of segemnts on line.
* NCOR maximum number of corners.
* SWVOID flag to indicate if there are voids.
* LINE line number.
* WEIGHT line weight.
* NUMERO region crossed by track.
* LENGHT length of track.
* SIGTAL albedo-cross section vector.
* NPSYS non-converged energy group indices.
*
*Parameters: output
* DPR CP matrix.
*
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
INTEGER NREG,NSOUT,NPIJ,NGRP,MXSEG,NCOR, NUMERO(*),
> ISD(6),ISF(6),LIN2C,IL,JL,IG,LINE,NOIL,IND1,
> IND2,I,J,IN0,IN1,IN2,NPSYS(NGRP),NUNK
REAL WEIGHT, SIGTAL(-NSOUT:NREG,NGRP)
DOUBLE PRECISION LENGHT(*), DPR(NPIJ,NGRP), XSIL, XSIL2
LOGICAL SWVOID
REAL ZERO, ONE, HALF
PARAMETER ( ZERO=0.0E0, ONE=1.0E0, HALF=0.5E0)
REAL SIXT,CUTEXP
PARAMETER ( CUTEXP=0.02)
*
* Allocated arrays
INTEGER, ALLOCATABLE, DIMENSION(:) :: NRSEG
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DSCBEG, DSCEND,
> SEGLEN, PRODUC
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: STAYIN, GOSOUT
IN0(I) = ((I+NSOUT+1)*(I+NSOUT+2))/2
IN1(I,J) = ((I+NSOUT+1)*(I+NSOUT))/2 + (J+NSOUT+1)
IN2(I,J) = (MAX(I+NSOUT+1,J+NSOUT+1)*
> (MAX(I+NSOUT+1,J+NSOUT+1)-1))/2
> + MIN(I+NSOUT+1,J+NSOUT+1)
*
* Scratch storage allocation
ALLOCATE(NRSEG(MXSEG))
ALLOCATE(DSCBEG(NGRP), DSCEND(NGRP), SEGLEN(MXSEG), PRODUC(NGRP))
ALLOCATE(STAYIN(NGRP,MXSEG),GOSOUT(NGRP,MXSEG))
*
SIXT=HALF/3.0
NUNK=NSOUT+NREG+1
*
* 0) REFORMAT TRACKING LINE
LIN2C= LINE-2*NCOR
DO 10 JL= 1, NCOR
ISD(JL)= NUMERO(JL)
ISF(JL)= NUMERO(NCOR+LIN2C+JL)
10 CONTINUE
DO 20 IL= 1, LIN2C
NRSEG(IL)= NUMERO(NCOR+IL)
SEGLEN(IL)= LENGHT(NCOR+IL)
20 CONTINUE
DO 90 IG= 1, NGRP
PRODUC(IG)= WEIGHT
90 CONTINUE
*
IF( SWVOID )THEN
*
* 1) VOIDS ARE POSSIBLE
* PII CALCULATION AND ESCAPE
DO 240 IL = 1,LINE-2*NCOR
NOIL = NRSEG(IL)
IND1 = IN0(NOIL)
DO 100 IG= 1, NGRP
IF(NPSYS(IG).EQ.0) GO TO 100
XSIL=SIGTAL(NOIL,IG)*SEGLEN(IL)
IF( XSIL.EQ.ZERO )THEN
GOSOUT(IG,IL)= ONE
STAYIN(IG,IL)= SEGLEN(IL)
DPR(IND1,IG)= DPR(IND1,IG)
> + HALF*WEIGHT*SEGLEN(IL)*SEGLEN(IL)
ELSE IF(XSIL .LT. CUTEXP) THEN
XSIL2=XSIL*XSIL
STAYIN(IG,IL)=XSIL-XSIL2*(HALF-SIXT*XSIL)
GOSOUT(IG,IL)=ONE-STAYIN(IG,IL)
PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL)
DPR(IND1,IG)= DPR(IND1,IG)
> + WEIGHT*XSIL2*(HALF-SIXT*XSIL)
ELSE
GOSOUT(IG,IL)= EXP( -XSIL )
STAYIN(IG,IL)= (ONE - GOSOUT(IG,IL))
PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL)
DPR(IND1,IG)= DPR(IND1,IG)
> + WEIGHT*(XSIL-STAYIN(IG,IL))
ENDIF
100 CONTINUE
240 CONTINUE
ELSE
DO 241 IL = 1,LINE-2*NCOR
NOIL = NRSEG(IL)
IND1 = IN0(NOIL)
DO 101 IG= 1, NGRP
IF(NPSYS(IG).EQ.0) GO TO 101
XSIL=SIGTAL(NOIL,IG)*SEGLEN(IL)
IF(XSIL .LT. CUTEXP) THEN
XSIL2=XSIL*XSIL
STAYIN(IG,IL)=XSIL-XSIL2*(HALF-SIXT*XSIL)
GOSOUT(IG,IL)=ONE-STAYIN(IG,IL)
PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL)
DPR(IND1,IG)= DPR(IND1,IG)
> + WEIGHT*XSIL2*(HALF-SIXT*XSIL)
ELSE
GOSOUT(IG,IL)= EXP( -XSIL )
STAYIN(IG,IL)= (ONE - GOSOUT(IG,IL))
PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL)
DPR(IND1,IG)= DPR(IND1,IG)
> + WEIGHT*(XSIL-STAYIN(IG,IL))
ENDIF
101 CONTINUE
241 CONTINUE
ENDIF
* PIJ CALCULATION
DO 120 IG= 1, NGRP
DSCBEG(IG)= WEIGHT
120 CONTINUE
DO 260 IL = 1, LINE-2*NCOR
NOIL = NRSEG(IL)
DO 130 IG= 1, NGRP
IF(NPSYS(IG).NE.0) DSCEND(IG)= WEIGHT*STAYIN(IG,IL)
130 CONTINUE
DO 250 JL = IL+1, LINE-2*NCOR
IND2= IN2(NRSEG(JL),NOIL)
DO 140 IG= 1, NGRP
IF(NPSYS(IG).EQ.0) GO TO 140
DPR(IND2,IG)= DPR(IND2,IG) + STAYIN(IG,JL)*DSCEND(IG)
DSCEND(IG)= DSCEND(IG)*GOSOUT(IG,JL)
140 CONTINUE
250 CONTINUE
* PIS CALCULATION
DO 261 JL = 1, NCOR
IND1= IN1(NOIL,ISD(JL))
IND2= IN1(NOIL,ISF(JL))
DO 150 IG= 1, NGRP
IF(NPSYS(IG).EQ.0) GO TO 150
DPR(IND1,IG)= DPR(IND1,IG) + DSCBEG(IG)*STAYIN(IG,IL)
DPR(IND2,IG)= DPR(IND2,IG) + DSCEND(IG)
150 CONTINUE
261 CONTINUE
DO 160 IG= 1, NGRP
IF(NPSYS(IG).NE.0) DSCBEG(IG)= DSCBEG(IG)*GOSOUT(IG,IL)
160 CONTINUE
260 CONTINUE
* PSS CALCULATION
DO 265 IL = 1, NCOR
DO 264 JL = 1, NCOR
IND2= IN2(ISD(IL),ISF(JL))
DO 170 IG = 1, NGRP
IF(NPSYS(IG).NE.0) DPR(IND2,IG)= DPR(IND2,IG) + PRODUC(IG)
170 CONTINUE
264 CONTINUE
265 CONTINUE
*
* Scratch storage deallocation
DEALLOCATE(GOSOUT,STAYIN)
DEALLOCATE(PRODUC,SEGLEN,DSCEND,DSCBEG)
DEALLOCATE(NRSEG)
*
RETURN
END
|
