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
|
*DECK DETCTL
SUBROUTINE DETCTL(NX,NY,NZ,NEL,VECT,RESP,NDET,XCT,YCT,ZCT,COR,
1 KEYF,IPRT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Call the subroutines that perform the parabolic interpolation.
*
*Copyright:
* Copyright (C) 2010 Ecole Polytechnique de Montreal.
*
*Author(s):
* E. Varin, M. Guyot
*
*Parameters:
* NX number of x mesh-splitted elements
* NY number of y mesh-splitted elements
* NZ number of z mesh-splitted elements
* NEL number of finite elements
* VECT
* RESP flux reads by the detector
* NDET number of detectors
* XCT center coordinates of each mesh-splitted elements for x
* YCT center coordinates of each mesh-splitted elements for y
* ZCT center coordinates of each mesh-splitted elements for z
* COR coordinates of the center of the detector
* KEYF keyflux recover from L_TRACK object
* IPRT printing index
*
*----------------------------------------------------------------------- *
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NX,NY,NZ,NEL,NDET,KEYF(NEL),IPRT
REAL VECT(*),COR(*),XCT,YCT,ZCT,RESP(NDET)
*----
* LOCAL VARIABLES
*----
REAL D1,D2,D3,X1,X2,X3,XX1,XX2,XX3,Y1,Y2,
1 Y3,YY1,YY2,YY3,Z1,Z2,Z3,ZZ1,ZZ2,ZZ3,PD1,PD2,PD3,PPD1,
2 PPD2,PPD3,CE,BE,AH
INTEGER I,III,NIJK,I1,I2,I3,IP1,IP2,IP3,J1,J2,J3,JP1,JP2,JP3,
1 K1,K2,K3,KP1,KP2,KP3,K0,JJJ
IF(IPRT.GT.4) WRITE(6,1000)
IF (NDET.LE.0) RETURN
NIJK = NX*NY
DO 10 III=1,NDET
I = (III-1)*3
D1 = COR(I+1)
D2 = COR(I+2)
D3 = COR(I+3)
*----
* DETERMINE CENTER OF INTERPOLATE RANGE
*----
CALL DETRTR(D1,XCT,NX,XX1,XX2,XX3,IP1,IP2,IP3)
X1 = XX1
X2 = XX2
X3 = XX3
I1 = IP1
I2 = IP2
I3 = IP3
CALL DETRTR(D2,YCT,NY,YY1,YY2,YY3,JP1,JP2,JP3)
Y1 = YY1
Y2 = YY2
Y3 = YY3
J1 = JP1
J2 = JP2
J3 = JP3
CALL DETRTR(D3,ZCT,NZ,ZZ1,ZZ2,ZZ3,KP1,KP2,KP3)
Z1 = ZZ1
Z2 = ZZ2
Z3 = ZZ3
K1 = KP1
K2 = KP2
K3 = KP3
IF (IPRT.GT.4) THEN
IF (MOD(III,25).EQ.0) WRITE(6,1000)
ENDIF
IF(IPRT.GT.4) THEN
WRITE(6,2000) III,D1,X1,X2,X3,D2,Y1,Y2,Y3,D3,Z1,Z2,Z3,
> I1,I2,I3, J1,J2,J3, K1,K2,K3
ENDIF
*----
* INTERPOLATION IN X AT PLANE Z=K1
*----
K0 = (K1-1)*NIJK
*----
* INTERPOLATION IN X AT PLANE Y=J1,Z=K1
*----
JJJ = NX*(J1-1)
PD1 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN X AT PLANE Y=J2,Z=K1
*----
JJJ = NX*(J2-1)
PD2 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN X AT PLANE Y=J3,Z=K1
*----
JJJ = NX*(J3-1)
PD3 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN Y AT PLANE Z=K1
*----
CALL DETPAR(Y1,Y2,Y3,PD1,PD2,PD3,AH,BE,CE)
PPD1 = AH*D2*D2 + BE*D2 + CE
*----
* INTERPOLATION IN X AT PLANE Z=K2
*----
K0 = (K2-1)*NIJK
*----
* INTERPOLATION IN X AT PLANE Y=J1,Z=K2
*----
JJJ = NX*(J1-1)
PD1 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN X AT PLANE Y=J2,Z=K2
*----
JJJ = NX*(J2-1)
PD2 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN X AT PLANE Y=J3,Z=K2
*----
JJJ = NX*(J3-1)
PD3 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN Y AT PLANE Z=K2
*----
CALL DETPAR2(Y1,Y2,Y3,PD1,PD2,PD3,AH,BE,CE)
PPD2 = AH*D2*D2 + BE*D2 + CE
*----
* INTERPOLATION IN X AT PLANE Z=K3
*----
K0 = (K3-1)*NIJK
*----
* INTERPOLATION IN X AT PLANE Y=J1,Z=K3
*----
JJJ = NX*(J1-1)
PD1 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN X AT PLANE Y=J2,Z=K3
*----
JJJ = NX*(J2-1)
PD2 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN X AT PLANE Y=J3,Z=K3
*----
JJJ = NX*(J3-1)
PD3 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
*----
* INTERPOLATION IN Y AT PLANE Z=K3
*----
CALL DETPAR2(Y1,Y2,Y3,PD1,PD2,PD3,AH,BE,CE)
PPD3 = AH*D2*D2 + BE*D2 + CE
*----
* INTERPOLATION IN Z
*----
CALL DETPAR2(Z1,Z2,Z3,PPD1,PPD2,PPD3,AH,BE,CE)
RESP(III) = AH*D3*D3 + BE*D3 + CE
10 CONTINUE
RETURN
1000 FORMAT(//,57X,'BRACKETING PROCESS',
> /,57X,'******************',
> //,5X,'DET',4X,'X ',8X,'X1',8X,'X2',8X,'X3',4X,
> 4X,'Y ',8X,'Y1',8X,'Y2',8X,'Y3',4X,
> 4X,'Z ',8X,'Z1',8X,'Z2',8X,'Z3',4X,/)
2000 FORMAT(5X,I3.3,12F10.4,/,5X,3X,3(10X,3(2X,I6.6,2X)))
END
|
