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
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
|
*DECK XELTI2
SUBROUTINE XELTI2( IPRT,IFTEMP,NANGLE,DENUSR,ISYMM,ANGLES,DENSTY,
> NTOTCL,MAXR,REMESH,LINMAX,RCUTOF,
> NSUR,NVOL,INDEL,MINDIM,
> MAXDIM,ICOORD,INCR,ICUR,TRKBEG,CONV,TRKDIR,
> LENGHT,NUMERO,DDENWT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Construct the sequential tape that will contain tracks for
* isotropic BC for 2-D calculation.
*
*Copyright:
* Copyright (C) 1989 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
* IPRT intermediate printing level for output.
* IFTEMP tracking file number.
* NANGLE number of angles used in the tracking process.
* DENUSR density of tracks in the plane perpendicular
* to the tracking angles.
* ISYMM flag for symetry:
* = 2 reflection plane normal to X axis;
* = 4 reflection plane normal to Y axis;
* = 8 reflection plane normal to X and Y axis.
* ANGLES 3D angle values.
* DENSTY density of tracks angle by angle.
* NTOTCL number of cylindres of a type + 2.
* MAXR max number of real mesh values in REMESH.
* REMESH real mesh values (rect/cyl).
* LINMAX max. number of track segments in a single track.
* RCUTOF cutof for corner tracking.
* NSUR number of surfaces.
* NVOL number of zones.
* INDEL numbering of surfaces & zones.
* MINDIM min index values for all axes (rect/cyl).
* MAXDIM max index values for all axes (rect/cyl).
* ICOORD principal axes direction (X/Y/Z) for meshes.
* ICUR current zonal location for a track segment.
* INCR increment direction for next track segment.
* TRKBEG position where a track begins.
* CONV segments of tracks.
* TRKDIR direction of a track in all axes.
* LENGHT relative lenght of each segment in a track.
* NUMERO material identification of each track segment.
* DDENWT density of tracks angle by angle.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*
* DECLARE DUMMY ARGUMENTS
INTEGER IPRT,IFTEMP,NANGLE,NTOTCL,MAXR,LINMAX
REAL TRKBEG(NTOTCL), TRKDIR(NTOTCL), CONV(NTOTCL),
> REMESH(MAXR),DENUSR
DOUBLE PRECISION DENSTY(NANGLE),ANGLES(2,NANGLE),LENGHT(LINMAX),
> DDENWT(NANGLE)
INTEGER MINDIM(NTOTCL), MAXDIM(NTOTCL), ICUR(NTOTCL),
> ICOORD(NTOTCL), INCR(NTOTCL), NUMERO(LINMAX),
> INDEL(4,*)
*
* DECLARE LOCAL VARIABLES
REAL TRKEND(3), TRKORI(3), BARY(2), PROJC2(3), RCUTOF,
> ANGEQN(2,2), ANGLE2(2), TRKCUT(3,2),
> TCUTOF(2,2), TORIC(3), DENLIN, R2CIRC, RCIRC, DP,
> ODDNXT, SURTOT, VOLTOT, DDENST, ANN, TOTLEN, TOTXXX
DOUBLE PRECISION WEIGHT,WTA,WTD
INTEGER NSBEG(4), NSEND(4), NSCUT(2), NDIM, I, NPOINT,
> NPO2, NCUTOF, NOTRAK, IANGL, IDIM, NANGLS, NESTIM,
> NSOLMX, IREF1, IZZ, IANG, NTTRK, NDEBS, I2, J,
> NSUR, NVOL, NCROS, K, K3, LINE, ISB, ISE, JC
INTEGER ISYMM
REAL WGA,WLA,WGD,WLD
EQUIVALENCE ( ANGEQN(1,2), ANGLE2 )
INTEGER IOUT
PARAMETER ( IOUT=6 )
CHARACTER TEDATA*13
*
* ONE WEIGHT FOR ALL TRACKS
* DENLIN= # OF TRACKS / CM
DENLIN= DENUSR
NDIM= 2
PROJC2(1)= 0.0
PROJC2(2)= 0.0
PROJC2(3)= 1.0
TRKBEG(3)= 0.0
TRKDIR(3)= 0.0
TRKEND(3)= 0.0
*
* COMPUTE THE CIRCUMSCRIBED RADIUS
* THE COORDINATE FOR THE TRUE CENTER OF THE CELL
R2CIRC= 0.0
DO 10 I = 1, 2
BARY(I)= 0.5 * (REMESH(MAXDIM(I)) + REMESH(MINDIM(I)))
R2CIRC= R2CIRC
> + (REMESH(MAXDIM(I)) - REMESH(MINDIM(I)))**2
10 CONTINUE
R2CIRC= 0.25 * R2CIRC
RCIRC = SQRT(R2CIRC)
*
* NPOINT= # OF TRACKS ALONG THE PERPENDICULAR AXIS
NPOINT= INT( 2. * RCIRC * DENLIN )
****** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE *****
****** CHANGE THIS "NPOINT" PARAMETER HAS TREMENDOUS EFFECTS ON TRACKING
****** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE *****
*
* POSSIBLE OTHER CHOICES (EXPLORED WITHOUT SUCCESS) ARE ==>
*1-) NPOINT= INT( 2. * RCIRC * DENLIN ) + 1
*2-) NPOINT= NINT( 2. * RCIRC * DENLIN )
*3-) NPOINT= NINT( 2. * RCIRC * DENLIN ) + 1
*
* KEEP "NPOINT" ODD & CORRECT DENSITY
NPO2 = NPOINT / 2
NPOINT= 2 * NPO2 + 1
DP = 2. * RCIRC / NPOINT
DENLIN= 1. / DP
DENUSR= DENLIN
ODDNXT= (2*NPO2+3) / (2.*RCIRC)
IF( RCUTOF.EQ.0.0 )THEN
NCUTOF= 1
ELSE
NCUTOF= 2
ENDIF
*
NOTRAK= 0
SURTOT= 0.0
VOLTOT= 0.0
DDENST= 1.0 / DENUSR
WEIGHT= 0.5*DDENST/DBLE(NANGLE)
CALL XELEQN( 2, 0, ANGEQN )
DO 15 IANGL= 1, NANGLE
CALL XELEQN( 2, NANGLE, ANGEQN )
DENSTY(IANGL)= REAL(2*NANGLE)
ANGLES(1,IANGL)= ANGEQN(1,1)
ANGLES(2,IANGL)= ANGEQN(2,1)
DDENWT(IANGL)=WEIGHT
15 CONTINUE
*----
* Optimize tracking by taking into account
* symmetry of geometry
* 1) Do not track symmetric lines
* 2) Nodify weight as required for tracks droped
*----
WGA=1.0
WLA=1.0
WGD=1.0
WLD=1.0
NANGLS= NANGLE
IF(ISYMM .EQ. 2 .OR.
> ISYMM .EQ. 4 .OR.
> ISYMM .EQ. 8 ) THEN
NANGLS= (NANGLE+1)/2
WGA=2.0
DO 200 IANGL=1,NANGLS-1
DENSTY(IANGL)=0.5*DENSTY(IANGL)
200 CONTINUE
IF(2*NANGLS .EQ. NANGLE) THEN
WLA=WGA
DENSTY(NANGLS)=0.5*DENSTY(NANGLS)
ENDIF
ENDIF
IF(ISYMM .EQ. 8 ) THEN
WGD=2.0
NPOINT=(NPOINT+1)/2
ENDIF
*
* COPY ANGLES AND DENSITIES ON TEMPORARY TRACKING FILE
WRITE(IFTEMP) ((ANGLES(IDIM,IANGL),IDIM=1,NDIM),IANGL=1,NANGLS)
WRITE(IFTEMP) (DENSTY(IANGL) ,IANGL=1,NANGLS)
*
* TO REINITIATE THE EQN ANGLES ROUTINE
CALL XELEQN( 2, 0, ANGEQN )
*
NSOLMX= 0
NDEBS= 0
IF( IPRT.GT.1 )THEN
*
* PREPARE & PRINT THE ESTIMATED NUMBER OF TRACKS
NESTIM= NPOINT * NANGLS
WRITE(IOUT,'(1H )')
WRITE(IOUT,'( 8H ECHO = ,I8,15H TRACKS/ANGLE )')
> NPOINT
WRITE(IOUT,'(30H ECHO = NEXT ODD DENSITY > ,F15.7,3H/CM )')
> ODDNXT
WRITE(IOUT,'( 8H ECHO = ,28H ESTIMATED NUMBER OF TRACKS= ,I8)')
> NESTIM
*
* PREPARE PRINTING WITH VARIABLE FORMAT
WRITE(IOUT,'(1H )')
WRITE(IOUT,6003) NANGLE
IF(NANGLS .NE. NANGLE) THEN
WRITE(IOUT,6004) NANGLS
ENDIF
NSOLMX= MIN(9, NANGLS/10)
IREF1=0
WRITE(IOUT,'( 1H0,10(I1,9X))') (IREF1, IZZ=0,NSOLMX)
WRITE(IOUT,'( 1H ,10(I1,9X))') (MOD(IZZ,10), IZZ=0,NSOLMX)
WRITE(IOUT,'( 2H 0)')
TEDATA= '(1H+,TXXX,I1)'
ENDIF
IANG=0
NTTRK=0
DO 290 IANGL= 1, NANGLS
IANG=IANG+1
WTA=WGA
IF(IANGL .EQ. NANGLS) THEN
WTA=WLA
ENDIF
IF( IPRT.GT.1 )THEN
IF( MOD(IANG,100) .EQ. 0 )THEN
IREF1=IREF1+1
NDEBS= NSOLMX+1
NSOLMX=MIN(NDEBS+9, NANGLS/10)
WRITE(IOUT,'( 1H0,10(I1,9X))')(IREF1,IZZ=NDEBS,NSOLMX)
WRITE(IOUT,'( 1H ,10(I1,9X))')
> (MOD(IZZ,10),IZZ=NDEBS,NSOLMX)
WRITE(IOUT,'( 2H 0)')
ELSE
IF( IPRT.GT.10000.AND.MOD(IANG,100).NE.0 )THEN
WRITE(IOUT,'( 1H ,10(I1,9X))')(IREF1,IZZ=NDEBS,NSOLMX)
WRITE(IOUT,'( 1H ,10(I1,9X))')
> (MOD(IZZ,10),IZZ=NDEBS,NSOLMX)
WRITE(IOUT,'( 2H )')
ENDIF
WRITE(TEDATA(7:9),'(I3.3)') MOD(IANG,100) + 2
WRITE(IOUT,TEDATA) MOD(IANG,10)
ENDIF
ENDIF
CALL XELEQN( 2, NANGLE, ANGEQN )
*
DO 40 I = 1, 2
TRKDIR(I)= ANGEQN(I,1)
*
* DETERMINE THE DIRECTION OF INCREASE FOR TRKDIR
INCR(I)= +1
IF( TRKDIR(I) .LT. 0.0 ) INCR(I)= -1
*
* MODIFY PERPENDICULAR ANGLES TO TAKE INTO ACCOUNT DP
ANGLE2(I)= DP * ANGLE2(I)
IF( NCUTOF.NE.1 )THEN
TCUTOF(I,1)= RCUTOF*ANGLE2(I)
TCUTOF(I,2)= -TCUTOF(I,1)
ENDIF
*
* DETERMINE THE ORIGINE OF ALL TRACKS
TRKORI(I)= BARY(I) - (NPO2+1) * ANGLE2(I)
40 CONTINUE
DO 180 I2 = 1, NPOINT
WTD=WGD
IF(I2 .EQ. NPOINT) THEN
WTD=WLD
ENDIF
ANN = 0.0
DO 50 J = 1, 2
TRKORI(J)= TRKORI(J) + ANGLE2(J)
ANN= ANN + (TRKORI(J)-BARY(J))**2
50 CONTINUE
*
* ELIMINATE TRACKS OUTSIDE CIRCUMSCRIBED CIRCLE
IF( ANN.GT.R2CIRC ) GO TO 180
*----
* Start test print
*
* WRITE(IOUT,7002) I2,I3,(TRKORI(JJ),JJ=1,3)
*7002 FORMAT(' ORIGINE MESH:',I10,5X,I10,5X,3(F11.5))
* Finish test print
*----
*
* WHICH EXTERNAL SURFACES DO THIS TRACK CROSS ?
NTTRK=NTTRK+1
CALL XELLSR( NDIM, NTOTCL, NSUR, MAXR, REMESH,
> INDEL, MINDIM, MAXDIM, ICOORD, ICUR, INCR,
> TRKORI, TRKDIR, TRKCUT, NSCUT, NCROS,
> TOTLEN)
*
* WHEN NOT SURFACES ARE CROSSED, ELIMINATE THE TRACK
IF( NCROS.LT.2 ) GO TO 180
DO 70 K= 1, NDIM
TRKBEG(K)= TRKCUT(K,1)
TRKEND(K)= TRKCUT(K,2)
70 CONTINUE
DO 75 K= 1, 4
NSBEG(K)= NSCUT(1)
NSEND(K)= NSCUT(2)
75 CONTINUE
IF( NCUTOF.NE.1 )THEN
DO 77 K= 1, 2
DO 76 K3= 1, 2
TORIC(K3)= TRKORI(K3)+TCUTOF(K3,K)
76 CONTINUE
CALL XELLSR( NDIM, NTOTCL, NSUR, MAXR, REMESH,
> INDEL, MINDIM, MAXDIM, ICOORD, ICUR, INCR,
> TORIC, TRKDIR, TRKCUT, NSCUT, NCROS,
> TOTXXX)
IF(NSCUT(1).NE.0) NSBEG(K)= NSCUT(1)
IF(NSCUT(2).NE.0) NSEND(K)= NSCUT(2)
77 CONTINUE
ENDIF
CALL XELLIN( NDIM, NTOTCL, MAXR, REMESH,
> NSUR, NVOL, INDEL, MINDIM, MAXDIM,
> ICOORD, ICUR, INCR, TRKBEG, TRKEND, TRKDIR,
> PROJC2, TOTLEN,
> CONV, LINMAX, LENGHT, NUMERO, LINE)
NOTRAK= NOTRAK+1
*
WRITE(IFTEMP) 1,LINE+2*NCUTOF,WEIGHT*WTA*WTD,IANG,
> (NSBEG(ISB),ISB=1,NCUTOF),
> (NUMERO(I),I=1,LINE),
> (NSEND(ISE),ISE=1,NCUTOF),
> ( DBLE(1.0/NCUTOF),ISB=1,NCUTOF),
> (LENGHT(I),I=1,LINE),
> ( DBLE(1.0/NCUTOF),ISE=1,NCUTOF)
IF( IPRT.GT.10000)THEN
WRITE(IOUT,6001) NOTRAK,
> NCUTOF,(TRKBEG(JC),JC=1,3),
> NCUTOF,(TRKEND(JC),JC=1,3),
> (TRKDIR(JC),JC=1,3)
WRITE(IOUT,6002) (1.0/NCUTOF,NSBEG(ISB),ISB=1,NCUTOF),
> (LENGHT(I), NUMERO(I),I=1,LINE),
> (1.0/NCUTOF,NSEND(ISE),ISE=1,NCUTOF)
ENDIF
180 CONTINUE
290 CONTINUE
IF( IPRT.GT.1 )THEN
WRITE(IOUT,'(1H )')
WRITE(IOUT,'(27H0ECHO = TRACKING PROPERTIES )')
WRITE(IOUT,'( 8H0ECHO = ,I3,20H ANGLES AND DENSITY:,
> F9.3,4H/CM )')
> NANGLS, DENUSR
WRITE(IOUT,'( 8H0ECHO = ,I10,3H / ,I10,
> 23H TRACKS STORED ON TAPE ,I2/)')
> NOTRAK,NTTRK,IFTEMP
ENDIF
RETURN
*
6001 FORMAT(' #',I8,1P,' B',I1,'(',2(E10.2,','),E10.2,')',
> ' E',I1,'(',2(E10.2,','),E10.2,')',
> ' D(',2(E10.2,','),E10.2,')' )
6002 FORMAT(1P,5(1X,E15.7,1X,I6))
6003 FORMAT(' '/
> ' ECHO = ',I3,' ANGLES TO BE TRACKED IN RANGE 0 TO PI')
6004 FORMAT(' ',I3,' ANGLES IN RANGE 0 TO PI/2 AFTER SYMMETRY')
END
|
