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
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
|
*DECK PIJI2D
SUBROUTINE PIJI2D(NREG,NSOUT,NSLINE,NCOR,SWVOID,SIGTAL,WEIGHT,
> SEGLEN,NRSEG,SEGPAT,DPR,
> MKI0,BIN0,PAS0,L0,
> MKI1,BIN1,PAS1,XLM1,L1,
> MKI2,BIN2,PAS2,XLM2)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Integration for general 2D 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.
* NSLINE number of segemnts on line.
* NCOR maximum number of corners.
* SWVOID flag to indicate if there are voids.
* SIGTAL albedo-cross section vector.
* WEIGHT line weight.
* SEGLEN length of track.
* NRSEG region crossed by track.
* MKI0 nb element quadratic BICKLEY table order N.
* BIN0 elements quadratic BICKLEY table order N.
* PAS0 step quadratic BICKLEY table order N.
* L0 log divergence quadratic BICKLEY table order N.
* MKI1 nb element quadratic BICKLEY table order N+1.
* BIN1 elements quadratic BICKLEY table order N+1.
* PAS1 step quadratic BICKLEY table order N+1.
* XLM1 upper limit quadratic BICKLEY table order N+1.
* L1 log divergence quadratic BICKLEY table order N+1.
* MKI2 nb element quadratic BICKLEY table order N+2.
* BIN2 elements quadratic BICKLEY table order N+2.
* PAS2 step quadratic BICKLEY table order N+2.
* XLM2 upper limit quadratic BICKLEY table order N+2.
*
*Parameters: output
* DPR CP matrix.
*
*Parameters: scratch
* SEGPAT optical path.
*
*Comments:
* PIJ => WITH BICKLEY FUNCTIONS OF ORDER 1,2,3
* PIJK => WITH BICKLEY FUNCTIONS OF ORDER 3,4,5
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* PARAMETERS
*----
REAL PI
PARAMETER (PI=3.1415926535897932)
*----
* INTERNAL FUNCTIONS
*----
DOUBLE PRECISION CBIN0,CBIN1,CBIN2
*----
* VARIABLES
*----
INTEGER NREG,NSOUT,NSLINE,NCOR,NRSEG(NSLINE)
LOGICAL SWVOID
REAL SIGTAL(-NSOUT:NREG),SEGPAT(NSLINE)
DOUBLE PRECISION WEIGHT,SEGLEN(NSLINE)
DOUBLE PRECISION DPR(-NSOUT:NREG,-NSOUT:NREG)
INTEGER MKI0,L0,MKI1,L1,MKI2
REAL BIN0(0:MKI0,3),PAS0,
> BIN1(0:MKI1,3),PAS1,XLM1,
> BIN2(0:MKI2,3),PAS2,XLM2
*----
* Local variables
*----
INTEGER IL,ISD,ISF,NSEG,ISEG,IREG,LPOS,ICSEG,JCSEG,
> JSEG,IR1
DOUBLE PRECISION FLOG0,FLOG1,ZREF1,ZREF2,ZREF3,XPOS,ZINTP,
> X3,ZI,X
*----
* QUADRATIC INTERPOLATION IN BICKLEY TABLES
*----
CBIN0(IL,X)= BIN0(IL,1)+X*(BIN0(IL,2)+X*BIN0(IL,3))
CBIN1(IL,X)= BIN1(IL,1)+X*(BIN1(IL,2)+X*BIN1(IL,3))
CBIN2(IL,X)= BIN2(IL,1)+X*(BIN2(IL,2)+X*BIN2(IL,3))
*----
* LOGARITHMIC DIVERGENCE FOR KI1 AND KI2
* FOR KI1 -> X*LOG(X) FOR K<L0 AND L0>0
* FOR KI2 -> -X**2*LOG(X)/2. FOR K<L1 AND L1>0
*----
FLOG0=FLOAT(L0/MAX(1,L0))
FLOG1=-0.5*FLOAT(L1/MAX(1,L1))
*----
* Process track required
*----
NSEG=NSLINE-2*NCOR
ZI=0.0D0
ZREF3=CBIN2(0,0.0D0)*WEIGHT
ZREF2=CBIN1(0,0.0D0)*WEIGHT
ZREF1=CBIN0(0,0.0D0)*WEIGHT
*----
* CONVERT SEGMENT LENGHT TO PATH LENGTH
*----
ICSEG=NCOR
DO 210 ISEG=1,NSEG
ICSEG=ICSEG+1
IREG=NRSEG(ICSEG)
SEGPAT(ISEG)=REAL(SEGLEN(ICSEG)*SIGTAL(IREG))
DPR(IREG,IREG)=DPR(IREG,IREG)+ZREF2*SEGPAT(ISEG)
210 CONTINUE
*----
* INTEGRATION
*----
XPOS=0.0D0
ZINTP=ZREF3
*----
* INTEGRATE OVER FIRST REGION AND COMPUTE PVS, PSS
*----
ICSEG=NCOR+1
DO 220 ISEG=1,NSEG
IREG=NRSEG(ICSEG)
IR1=NRSEG(NCOR+1)
XPOS=XPOS+SEGPAT(ISEG)
LPOS=MIN(NINT(XPOS*PAS2),MKI2)
ZI=WEIGHT*CBIN2(LPOS,XPOS)
X3=ZI-ZINTP
DPR(IR1,IREG)=DPR(IR1,IREG)+X3
DO 225 ISD=1,NCOR
DPR(IREG,NRSEG(ISD))=DPR(IREG,NRSEG(ISD))-X3
225 CONTINUE
IF(XPOS.GT.XLM2) GO TO 221
ZINTP=ZI
ICSEG=ICSEG+1
220 CONTINUE
221 CONTINUE
IR1=NRSEG(NCOR+1)
DO 226 ISF=NSLINE-NCOR+1,NSLINE
DPR(IR1,NRSEG(ISF))=DPR(IR1,NRSEG(ISF))-ZI
DO 227 ISD=1,NCOR
DPR(NRSEG(ISD),NRSEG(ISF))=
> DPR(NRSEG(ISD),NRSEG(ISF))+ZI
227 CONTINUE
226 CONTINUE
ICSEG=NCOR+2
DO 230 ISEG=2,NSEG
XPOS=0.0D0
ZINTP=ZREF3
JCSEG=ICSEG
DO 240 JSEG=ISEG,NSEG
XPOS=XPOS+SEGPAT(JSEG)
LPOS=MIN(NINT(XPOS*PAS2),MKI2)
ZI=WEIGHT*CBIN2(LPOS,XPOS)
X3=ZI-ZINTP
DPR(NRSEG(ICSEG-1),NRSEG(JCSEG))=
> DPR(NRSEG(ICSEG-1),NRSEG(JCSEG))-X3
DPR(NRSEG(ICSEG),NRSEG(JCSEG))=
> DPR(NRSEG(ICSEG),NRSEG(JCSEG))+X3
IF(XPOS.GT.XLM2) GO TO 241
ZINTP=ZI
JCSEG=JCSEG+1
240 CONTINUE
241 CONTINUE
DO 235 ISF=NSLINE-NCOR+1,NSLINE
DPR(NRSEG(ICSEG),NRSEG(ISF))=
> DPR(NRSEG(ICSEG),NRSEG(ISF))-ZI
DPR(NRSEG(ICSEG-1),NRSEG(ISF))=
> DPR(NRSEG(ICSEG-1),NRSEG(ISF))+ZI
235 CONTINUE
ICSEG=ICSEG+1
230 CONTINUE
ICSEG=NCOR+NSEG
DO 236 ISF=NSLINE-NCOR+1,NSLINE
DPR(NRSEG(ICSEG),NRSEG(ISF))=
> DPR(NRSEG(ICSEG),NRSEG(ISF))+ZREF3
236 CONTINUE
*----
* FOR VOID REGIONS RESET PROBABILITIES
*----
IF(SWVOID) THEN
ICSEG=NCOR+1
DO 300 ISEG=1,NSEG
IF(SIGTAL(NRSEG(ICSEG)).NE.0.0) GO TO 301
XPOS=0.0D0
DPR(NRSEG(ICSEG),NRSEG(ICSEG))=
> DPR(NRSEG(ICSEG),NRSEG(ICSEG))
> + 0.5*ZREF1*SEGLEN(ICSEG)*SEGLEN(ICSEG)
ZINTP=ZREF2*SEGLEN(ICSEG)
JCSEG=ICSEG+1
DO 310 JSEG=ISEG+1,NSEG
IF(SIGTAL(NRSEG(JCSEG)).EQ.0.0) THEN
LPOS=MIN(NINT(XPOS*PAS0),MKI0)
IF(LPOS.LT.L0.AND.XPOS.NE.0.0) THEN
ZI=WEIGHT*SEGLEN(ICSEG)*SEGLEN(JCSEG)*
> (CBIN0(LPOS,XPOS)+FLOG0*XPOS*LOG(XPOS))
ELSE
ZI=WEIGHT*SEGLEN(ICSEG)*SEGLEN(JCSEG)*CBIN0(LPOS,XPOS)
ENDIF
X3=0.5*ZI
ELSE
XPOS=XPOS+SEGPAT(JSEG)
LPOS=MIN(NINT(XPOS*PAS1),MKI1)
IF(LPOS.LT.L1.AND.XPOS.NE.0.0) THEN
ZI=WEIGHT*SEGLEN(ICSEG)*
> (CBIN1(LPOS,XPOS)+FLOG1*XPOS*XPOS*LOG(XPOS))
ELSE
ZI=WEIGHT*SEGLEN(ICSEG)*CBIN1(LPOS,XPOS)
ENDIF
X3=ZINTP-ZI
ZINTP=ZI
ENDIF
DPR(NRSEG(ICSEG),NRSEG(JCSEG))=
> DPR(NRSEG(ICSEG),NRSEG(JCSEG))+X3
IF(XPOS.GT.XLM1) GO TO 311
JCSEG=JCSEG+1
310 CONTINUE
311 CONTINUE
DO 320 ISF=NSLINE-NCOR+1,NSLINE
DPR(NRSEG(ICSEG),NRSEG(ISF))=
> DPR(NRSEG(ICSEG),NRSEG(ISF))+ZINTP
320 CONTINUE
XPOS=0.0D0
ZINTP=ZREF2*SEGLEN(ICSEG)
JCSEG=ICSEG-1
DO 330 JSEG=ISEG-1,1,-1
IF(SIGTAL(NRSEG(JCSEG)).EQ.0.0) THEN
LPOS=MIN(NINT(XPOS*PAS0),MKI0)
IF(LPOS.LT.L0.AND.XPOS.NE.0.0) THEN
ZI=WEIGHT*SEGLEN(ICSEG)*SEGLEN(JCSEG)*
> (CBIN0(LPOS,XPOS)+FLOG0*XPOS*LOG(XPOS))
ELSE
ZI=WEIGHT*SEGLEN(ICSEG)*SEGLEN(JCSEG)*CBIN0(LPOS,XPOS)
ENDIF
X3=0.5*ZI
ELSE
XPOS=XPOS+SEGPAT(JSEG)
LPOS=MIN(NINT(XPOS*PAS1),MKI1)
IF(LPOS.LT.L1.AND.XPOS.NE.0.0) THEN
ZI=WEIGHT*SEGLEN(ICSEG)*
> (CBIN1(LPOS,XPOS)+FLOG1*XPOS*XPOS*LOG(XPOS))
ELSE
ZI=WEIGHT*SEGLEN(ICSEG)*CBIN1(LPOS,XPOS)
ENDIF
X3=ZINTP-ZI
ZINTP=ZI
ENDIF
DPR(NRSEG(ICSEG),NRSEG(JCSEG))=
> DPR(NRSEG(ICSEG),NRSEG(JCSEG))+X3
IF(XPOS.GT.XLM1) GO TO 331
JCSEG=JCSEG-1
330 CONTINUE
331 CONTINUE
DO 340 ISD=1,NCOR
DPR(NRSEG(ICSEG),NRSEG(ISD))=
> DPR(NRSEG(ICSEG),NRSEG(ISD))+ZINTP
340 CONTINUE
301 CONTINUE
ICSEG=ICSEG+1
300 CONTINUE
ENDIF
*----
* RETURN
*----
RETURN
END
|
