]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | *PATCH,GEXAM1. |
2 | *CMZ : 3.21/02 29/03/94 15.41.35 by S.Giani | |
3 | *DECK, GTCKOV | |
4 | *CMZ : 12/09/95 11.07.14 by S.Ravndal | |
5 | *CMZ : 3.21/04 13/12/94 15.17.13 by S.Giani | |
6 | *-- Author : | |
7 | SUBROUTINE GTCKOV | |
8 | C. | |
9 | C. ****************************************************************** | |
10 | C. * * | |
11 | C. * This routine is called to follow the Cherenkov photons * | |
12 | C. * created during the tracking of charged particles and * | |
13 | C. * simulate the relevant processes along the way, until either * | |
14 | C. * the photon is absorbed or exits the detector. Processes * | |
15 | C. * currently simulated are absorption in-flight, and reflection * | |
16 | C. * /transmission/absorption at a medium boundary. There are two * | |
17 | C. * boundary types: dielectric-metal and dielectric-dielectric. * | |
18 | C. * For each of these there is a continuum of reflectivity * | |
19 | C. * and of surface quality from mirror finish to matte. The * | |
20 | C. * surface model is contained in routine GHSURF. * | |
21 | C. * * | |
22 | C. * ==>Called by : GTRACK * | |
23 | C. * Authors F.Carminati, R.Jones ************ * | |
24 | C. * * | |
25 | C. ****************************************************************** | |
26 | C. | |
27 | #include "geant321/gcbank.inc" | |
28 | #include "geant321/gccuts.inc" | |
29 | #include "geant321/gcjloc.inc" | |
30 | #include "geant321/gconsp.inc" | |
31 | #include "geant321/gcphys.inc" | |
32 | #include "geant321/gcstak.inc" | |
33 | #include "geant321/gctmed.inc" | |
34 | #include "geant321/gctrak.inc" | |
35 | #include "geant321/gcvolu.inc" | |
36 | #include "geant321/gcnum.inc" | |
37 | #include "geant321/gcvol1.inc" | |
38 | #include "geant321/gcunit.inc" | |
39 | ||
40 | #if !defined(CERNLIB_OLD) | |
41 | #include "geant321/gcvdma.inc" | |
42 | #endif | |
43 | ||
44 | * | |
45 | * ** The following common is in GTMEDI. LSAMVL is set to true if | |
46 | * ** we did not change volume yet | |
47 | * | |
48 | COMMON/GCCHAN/LSAMVL | |
49 | LOGICAL LSAMVL | |
50 | * | |
51 | DIMENSION R(3),D(3),U(3),QQ(3),vin(3) | |
52 | DIMENSION VBOU(3) | |
53 | LOGICAL LOLDTR | |
54 | #if !defined(CERNLIB_SINGLE) | |
55 | PARAMETER (EPSMAC=1.E-6) | |
56 | #endif | |
57 | #if defined(CERNLIB_SINGLE) | |
58 | PARAMETER (EPSMAC=1.E-11) | |
59 | #endif | |
60 | PARAMETER (MXPUSH=10) | |
61 | SAVE RIN1,EFFIC | |
62 | C. | |
63 | C. ------------------------------------------------------------------ | |
64 | * | |
65 | * *** Update local pointers if medium has changed | |
66 | * | |
67 | LOLDTR=.FALSE. | |
68 | IF(IUPD.EQ.0)THEN | |
69 | IUPD = 1 | |
70 | JTCKOV = LQ(JTM-3) | |
71 | IF(JTCKOV.EQ.0) THEN | |
72 | * | |
73 | * *** This Cerenkov photon has crossed into a black medium. | |
74 | * *** Just absorb it. | |
75 | IPROC = 101 | |
76 | STEP = 0. | |
77 | SLABS = 0. | |
78 | ISTOP = 2 | |
79 | DESTEP = 0. | |
80 | GOTO 110 | |
81 | ENDIF | |
82 | NPCKOV = Q(JTCKOV+1) | |
83 | JABSCO = LQ(JTCKOV-1) | |
84 | JEFFIC = LQ(JTCKOV-2) | |
85 | JINDEX = LQ(JTCKOV-3) | |
86 | JPOLAR = LQ(JSTAK-1) | |
87 | ENDIF | |
88 | IF(SLENG.LE.0.) THEN | |
89 | * | |
90 | * *** Calculate GEKRAT for the particle | |
91 | IF(VECT(7).GE.Q(JTCKOV+NPCKOV+1)) THEN | |
92 | GEKRAT=1. | |
93 | IEKBIN=NPCKOV-1 | |
94 | ELSEIF(VECT(7).LT.Q(JTCKOV+2)) THEN | |
95 | * | |
96 | * *** Particle below energy threshold ? Short circuit | |
97 | * *** This should never happen because the photons are generated | |
98 | * *** only above threshold | |
99 | * | |
100 | * GEKIN = 0. | |
101 | * GETOT = 0. | |
102 | * VECT(7)= 0. | |
103 | * ISTOP = 2 | |
104 | * NMEC = 1 | |
105 | * LMEC(1)= 30 | |
106 | * GO TO 110 | |
107 | ||
108 | GEKRAT=0. | |
109 | IEKBIN=1 | |
110 | ELSE | |
111 | JMIN = 1 | |
112 | JMAX = NPCKOV | |
113 | 10 JMED = (JMIN+JMAX)/2 | |
114 | IF(Q(JTCKOV+JMED+1).LT.VECT(7)) THEN | |
115 | JMIN = JMED | |
116 | ELSE | |
117 | JMAX = JMED | |
118 | ENDIF | |
119 | IF(JMAX-JMIN.GT.1) GO TO 10 | |
120 | IEKBIN = JMIN | |
121 | GEKRAT = (VECT(7) - Q(JTCKOV+IEKBIN+1))/ | |
122 | + (Q(JTCKOV+IEKBIN+2)-Q(JTCKOV+IEKBIN+1)) | |
123 | ENDIF | |
124 | GEKRT1=1.-GEKRAT | |
125 | RIN1=Q(JINDEX+IEKBIN)*GEKRT1+Q(JINDEX+IEKBIN+1)*GEKRAT | |
126 | EFFIC=Q(JEFFIC+IEKBIN)*GEKRT1+Q(JEFFIC+IEKBIN+1)*GEKRAT | |
127 | STEPLA=Q(JABSCO+IEKBIN)*GEKRT1+Q(JABSCO+IEKBIN+1)*GEKRAT | |
128 | ENDIF | |
129 | * | |
130 | * *** Compute current step size | |
131 | * | |
132 | IPROC = 103 | |
133 | STEP = STEMAX | |
134 | * | |
135 | * ** Step limitation due to in flight absorbtion ? | |
136 | * | |
137 | IF (ILABS.GT.0) THEN | |
138 | SLABS = STEPLA*ZINTLA | |
139 | IF (SLABS.LT.STEP) THEN | |
140 | STEP = SLABS | |
141 | IPROC = 101 | |
142 | ENDIF | |
143 | ENDIF | |
144 | * | |
145 | IF (STEP.LT.0.) STEP = 0. | |
146 | * | |
147 | * ** Step limitation due to geometry ? | |
148 | * | |
149 | STEPT=0. | |
150 | IF (STEP.GE.SAFETY) THEN | |
151 | CALL GTNEXT | |
152 | IF (IGNEXT.NE.0) THEN | |
153 | * | |
154 | * ** We are going to cross a boundary, so we need to simulate | |
155 | * ** boundary effects and to know what is on the other side. | |
156 | * ** For the moment save the current vector in the geometry tree. | |
157 | * | |
158 | #if !defined(CERNLIB_OLD) | |
159 | if(mycoun.gt.1.and.nfmany.gt.0)then | |
160 | nlevel=manyle(nfmany) | |
161 | do 99 i=1,nlevel | |
162 | names(i)=manyna(nfmany,i) | |
163 | number(i)=manynu(nfmany,i) | |
164 | 99 continue | |
165 | call glvolu(nlevel,names,number,ier) | |
166 | if(ier.ne.0)print *,'Fatal error in GLVOLU' | |
167 | ingoto=0 | |
168 | endif | |
169 | #endif | |
170 | NLEVL1 = NLEVEL | |
171 | DO 20 I=1,NLEVEL | |
172 | NAMES1(I) = NAMES(I) | |
173 | NUMBR1(I) = NUMBER(I) | |
174 | LVOLU1(I) = LVOLUM(I) | |
175 | 20 CONTINUE | |
176 | * | |
177 | * *** This is different from the other tracking routines. | |
178 | * *** We get to the boundary and then we just jump over it | |
179 | * *** So, linear transport till we are very near the boundary | |
180 | * | |
181 | #if !defined(CERNLIB_IBM) | |
182 | STEP = MAX(SNEXT-PREC,0.) | |
183 | #endif | |
184 | #if defined(CERNLIB_IBM) | |
185 | STEP = MAX(SNEXT-2.*PREC,0.) | |
186 | #endif | |
187 | C | |
188 | C In case of Cherenkovs beeing near a corner, particle holding is | |
189 | C prevented by a bigger step size. The value of VBOU is only used to | |
190 | C calculate the correct surface normal (by GLISUR and GGPERP). The | |
191 | C tracking of the Cherenkov is still using STEP | |
192 | C | |
193 | IF(SNEXT.GT.0.) THEN | |
194 | DO 25 I=1,3 | |
195 | VBOU(I)=VECT(I)+SNEXT*VECT(I+3) | |
196 | 25 CONTINUE | |
197 | END IF | |
198 | C | |
199 | IF(STEP.GT.0.) THEN | |
200 | DO 30 I=1,3 | |
201 | VECT(I)=VECT(I)+STEP*VECT(I+3) | |
202 | 30 CONTINUE | |
203 | ENDIF | |
204 | STEPT=STEP | |
205 | #if !defined(CERNLIB_IBM) | |
206 | STEP = SNEXT - STEP + PREC | |
207 | #endif | |
208 | #if defined(CERNLIB_IBM) | |
209 | STEP = SNEXT - STEP + 2.*PREC | |
210 | #endif | |
211 | IPROC = 0 | |
212 | INWVOL= 2 | |
213 | NMEC = 1 | |
214 | LMEC(1)=1 | |
215 | ENDIF | |
216 | * | |
217 | * Update SAFETY in stack companions, if any | |
218 | * This may well not work. | |
219 | IF (IQ(JSTAK+3).NE.0) THEN | |
220 | DO 40 IST = IQ(JSTAK+3),IQ(JSTAK+1) | |
221 | JST = JSTAK +3 +(IST-1)*NWSTAK | |
222 | Q(JST+11) = SAFETY | |
223 | 40 CONTINUE | |
224 | IQ(JSTAK+3) = 0 | |
225 | ENDIF | |
226 | * | |
227 | ELSE | |
228 | IQ(JSTAK+3) = 0 | |
229 | ENDIF | |
230 | * | |
231 | * *** Linear transport | |
232 | * | |
233 | VIN(1) = VECT(1) | |
234 | VIN(2) = VECT(2) | |
235 | VIN(3) = VECT(3) | |
236 | IF (INWVOL.EQ.2) THEN | |
237 | NBPUSH = 0 | |
238 | 50 DO 60 I = 1,3 | |
239 | VECTMP = VECT(I) +STEP*VECT(I+3) | |
240 | IF(VECTMP.EQ.VECT(I)) THEN | |
241 | * | |
242 | * *** Correct for machine precision | |
243 | * | |
244 | IF(VECT(I+3).NE.0.) THEN | |
245 | VECTMP = VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))* | |
246 | + EPSMAC | |
247 | IF(NMEC.GT.0) THEN | |
248 | IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1 | |
249 | ENDIF | |
250 | NMEC=NMEC+1 | |
251 | LMEC(NMEC)=104 | |
252 | #if defined(CERNLIB_DEBUG) | |
253 | WRITE(CHMAIL, 10000) | |
254 | CALL GMAIL(0,0) | |
255 | WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT | |
256 | CALL GMAIL(0,0) | |
257 | 10000 FORMAT(' Boundary correction in GTCKOV: ', | |
258 | + ' GEKIN NUMED STEP SNEXT') | |
259 | 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X) | |
260 | #endif | |
261 | ENDIF | |
262 | ENDIF | |
263 | VOUT(I) = VECTMP | |
264 | 60 CONTINUE | |
265 | NUMED1=NUMED | |
266 | CALL GTMEDI(VOUT,NUMED2) | |
267 | LOLDTR=.TRUE. | |
268 | IF(NUMED2.LE.0) THEN | |
269 | VECT(1) = VOUT(1) | |
270 | VECT(2) = VOUT(2) | |
271 | VECT(3) = VOUT(3) | |
272 | GO TO 110 | |
273 | ENDIF | |
274 | JVO=LQ(JVOLUM-LVOLUM(NLEVEL)) | |
275 | IF(LSAMVL.AND.Q(JVO+2).NE.12.) THEN | |
276 | * | |
277 | * *** In spite of our efforts we have not crossed the boundary | |
278 | * *** we increase the step size and try again | |
279 | * | |
280 | NBPUSH = NBPUSH + 1 | |
281 | IF (NBPUSH.LE.MXPUSH) THEN | |
282 | STEP = STEP + NBPUSH*PREC | |
283 | GOTO 50 | |
284 | ELSE | |
285 | INWVOL = 0 | |
286 | ENDIF | |
287 | ||
288 | ENDIF | |
289 | IF(NUMED1.EQ.NUMED2) THEN | |
290 | * | |
291 | * *** If we are in the same medium, nothing needs to be done! | |
292 | * | |
293 | VECT(1)=VOUT(1) | |
294 | VECT(2)=VOUT(2) | |
295 | VECT(3)=VOUT(3) | |
296 | IPROC=0 | |
297 | ELSE | |
298 | JTM2 = LQ(JTMED-NUMED2) | |
299 | IF(IQ(JTM2-2).GE.3) THEN | |
300 | JTCKV2 = LQ(JTM2-3) | |
301 | ELSE | |
302 | JTCKV2 = 0 | |
303 | ENDIF | |
304 | IF(JTCKV2.GT.0) THEN | |
305 | NPCKV2 = Q(JTCKV2+1) | |
306 | JABSC2 = LQ(JTCKV2-1) | |
307 | JEFFI2 = LQ(JTCKV2-2) | |
308 | JINDX2 = LQ(JTCKV2-3) | |
309 | IF(VECT(7).GE.Q(JTCKV2+NPCKV2+1)) THEN | |
310 | GEKRT2=1. | |
311 | IEKBI2=NPCKV2-1 | |
312 | ELSEIF(VECT(7).LT.Q(JTCKV2+2)) THEN | |
313 | GEKRT2=0. | |
314 | IEKBI2=1 | |
315 | ELSE | |
316 | JMIN = 1 | |
317 | JMAX = NPCKV2 | |
318 | 64 JMED = (JMIN+JMAX)/2 | |
319 | IF(Q(JTCKV2+JMED+1).LT.VECT(7)) THEN | |
320 | JMIN = JMED | |
321 | ELSE | |
322 | JMAX = JMED | |
323 | ENDIF | |
324 | IF(JMAX-JMIN.GT.1) GO TO 64 | |
325 | IEKBI2 = JMIN | |
326 | GEKRT2 = (VECT(7) - Q(JTCKV2+IEKBI2+1))/ | |
327 | + (Q(JTCKV2+IEKBI2+2)-Q(JTCKV2+IEKBI2+1)) | |
328 | ENDIF | |
329 | GEKRT1=1.-GEKRT2 | |
330 | ABSCO2=Q(JABSC2+IEKBI2)*GEKRT1+Q(JABSC2+IEKBI2+1)*GEKRT2 | |
331 | EFFIC2=Q(JEFFI2+IEKBI2)*GEKRT1+Q(JEFFI2+IEKBI2+1)*GEKRT2 | |
332 | IF(JINDX2.GT.0) THEN | |
333 | RIN2=Q(JINDX2+IEKBI2)*GEKRT1+Q(JINDX2+IEKBI2+1)*GEKRT2 | |
334 | ELSE | |
335 | RIN2=0. | |
336 | ENDIF | |
337 | IPROC = 102 | |
338 | ELSE | |
339 | ISTOP=2 | |
340 | DESTEP=0 | |
341 | NMEC = NMEC+1 | |
342 | LMEC(NMEC)=30 | |
343 | VECT(1)=VOUT(1) | |
344 | VECT(2)=VOUT(2) | |
345 | VECT(3)=VOUT(3) | |
346 | GOTO 110 | |
347 | ENDIF | |
348 | ENDIF | |
349 | ELSE | |
350 | DO 70 I = 1,3 | |
351 | VECT(I) = VECT(I) +STEP*VECT(I+3) | |
352 | 70 CONTINUE | |
353 | ENDIF | |
354 | * | |
355 | STEP = STEPT + STEP | |
356 | SLENG = SLENG +STEP | |
357 | * | |
358 | * *** Update time of flight | |
359 | * | |
360 | TOFG = TOFG +STEP*RIN1/CLIGHT | |
361 | * | |
362 | * *** Update interaction probabilities | |
363 | * | |
364 | IF (ILABS.GT.0) ZINTLA = ZINTLA -STEP/STEPLA | |
365 | * | |
366 | IF (IPROC.EQ.0) GO TO 110 | |
367 | NMEC = NMEC+1 | |
368 | LMEC(NMEC) = IPROC | |
369 | * | |
370 | * ** Absorbtion in flight ? | |
371 | * | |
372 | IF (IPROC.EQ.101) THEN | |
373 | ISTOP=2 | |
374 | CALL GRNDM(RNDM,1) | |
375 | IF(RNDM.LT.EFFIC) THEN | |
376 | * | |
377 | * *** Destep =/= 0 means that the photon has been detected | |
378 | * | |
379 | DESTEP=VECT(7) | |
380 | ELSE | |
381 | DESTEP=0. | |
382 | ENDIF | |
383 | * | |
384 | * ** Boundary action? | |
385 | * | |
386 | ELSE IF (IPROC.EQ.102) THEN | |
387 | IF(JINDX2.EQ.0) THEN | |
388 | * | |
389 | * *** Case dielectric -> metal | |
390 | * | |
391 | CALL GRNDM(RNDM,1) | |
392 | IF(RNDM.LT.ABSCO2) THEN | |
393 | * | |
394 | * *** Photon is absorbed in the next medium | |
395 | * | |
396 | NMEC=NMEC+1 | |
397 | LMEC(NMEC)=101 | |
398 | ISTOP = 2 | |
399 | CALL GRNDM(RNDM,1) | |
400 | IF(RNDM.LT.EFFIC2) THEN | |
401 | * | |
402 | * *** Destep =/= 0 means that the photon has been detected | |
403 | * | |
404 | DESTEP=VECT(7) | |
405 | ELSE | |
406 | DESTEP = 0. | |
407 | END IF | |
408 | VECT(1) = VOUT(1) | |
409 | VECT(2) = VOUT(2) | |
410 | VECT(3) = VOUT(3) | |
411 | GOTO 110 | |
412 | ELSE | |
413 | * | |
414 | * *** Photon is reflected (no polarization effects) | |
415 | * | |
416 | CALL GLISUR(VECT,VBOU,NUMED1,NUMED2,U,PDOTU,IERR) | |
417 | IF (IERR.NE.0) THEN | |
418 | WRITE(CHMAIL,10200) IERR | |
419 | CALL GMAIL(0,0) | |
420 | GO TO 110 | |
421 | END IF | |
422 | * | |
423 | * ** Restore old volume tree, the photon does not cross the boundary | |
424 | * | |
425 | * CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR) | |
426 | CALL GTMEDI(VIN,N) | |
427 | LOLDTR=.FALSE. | |
428 | * | |
429 | * | |
430 | * *** Reflect, but make sure we are in the old volume before | |
431 | * | |
432 | NBPUSH = 0 | |
433 | 80 NBPUSH = NBPUSH+1 | |
434 | IF(NBPUSH.GT.MXPUSH) THEN | |
435 | WRITE(CHMAIL,10300) NTMULT,NSTEP | |
436 | CALL GMAIL(0,0) | |
437 | ISTOP=1 | |
438 | GOTO 110 | |
439 | ELSE | |
440 | CALL GINVOL(VECT,ISAME) | |
441 | IF(ISAME.EQ.0) THEN | |
442 | PRECN = NBPUSH*PREC | |
443 | VECT(1) = VECT(1) - PRECN*VECT(4) | |
444 | VECT(2) = VECT(2) - PRECN*VECT(5) | |
445 | VECT(3) = VECT(3) - PRECN*VECT(6) | |
446 | GO TO 80 | |
447 | ENDIF | |
448 | ENDIF | |
449 | * | |
450 | NMEC=NMEC+1 | |
451 | LMEC(NMEC)=106 | |
452 | EDOTU = POLAR(1)*U(1) + POLAR(2)*U(2) + POLAR(3)*U(3) | |
453 | VECT(4) = +VECT(4) - 2*PDOTU*U(1) | |
454 | VECT(5) = +VECT(5) - 2*PDOTU*U(2) | |
455 | VECT(6) = +VECT(6) - 2*PDOTU*U(3) | |
456 | POLAR(1) = -POLAR(1) + 2*EDOTU*U(1) | |
457 | POLAR(2) = -POLAR(2) + 2*EDOTU*U(2) | |
458 | POLAR(3) = -POLAR(3) + 2*EDOTU*U(3) | |
459 | INWVOL = 0 | |
460 | ENDIF | |
461 | ELSE | |
462 | * | |
463 | * case dielectric-dielectric: | |
464 | * | |
465 | CALL GLISUR(VECT,VBOU,NUMED1,NUMED2,U,PDOTU,IERR) | |
466 | IF (IERR.NE.0) THEN | |
467 | WRITE(CHMAIL,10200) IERR | |
468 | CALL GMAIL(0,0) | |
469 | GO TO 110 | |
470 | END IF | |
471 | EDOTU = POLAR(1)*U(1) + POLAR(2)*U(2) + POLAR(3)*U(3) | |
472 | COST1 = -PDOTU | |
473 | IF (ABS(COST1).LT.1.) THEN | |
474 | SINT1 = SQRT((1-COST1)*(1+COST1)) | |
475 | SINT2 = SINT1*RIN1/RIN2 | |
476 | ELSE | |
477 | SINT1 = 0.0 | |
478 | SINT2 = 0.0 | |
479 | END IF | |
480 | IF (SINT2.GE.1) THEN | |
481 | * | |
482 | * *** Simulate total internal reflection | |
483 | * | |
484 | NMEC=NMEC+1 | |
485 | LMEC(NMEC)=106 | |
486 | * | |
487 | * ** Restore old volume tree, the photon does not cross the boundary | |
488 | * | |
489 | * CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR) | |
490 | CALL GTMEDI(VIN,N) | |
491 | LOLDTR=.FALSE. | |
492 | * | |
493 | * *** Reflect, but make sure we are in the old volume before | |
494 | * | |
495 | NBPUSH = 0 | |
496 | 90 NBPUSH = NBPUSH+1 | |
497 | IF(NBPUSH.GT.MXPUSH) THEN | |
498 | WRITE(CHMAIL,10300) NTMULT,NSTEP | |
499 | CALL GMAIL(0,0) | |
500 | ISTOP=1 | |
501 | GOTO 110 | |
502 | ELSE | |
503 | CALL GINVOL(VECT,ISAME) | |
504 | IF(ISAME.EQ.0) THEN | |
505 | PRECN = NBPUSH*PREC | |
506 | VECT(1) = VECT(1) - PRECN*VECT(4) | |
507 | VECT(2) = VECT(2) - PRECN*VECT(5) | |
508 | VECT(3) = VECT(3) - PRECN*VECT(6) | |
509 | GO TO 90 | |
510 | ENDIF | |
511 | ENDIF | |
512 | * | |
513 | VECT(4) = +VECT(4) - 2*PDOTU*U(1) | |
514 | VECT(5) = +VECT(5) - 2*PDOTU*U(2) | |
515 | VECT(6) = +VECT(6) - 2*PDOTU*U(3) | |
516 | POLAR(1) = -POLAR(1) + 2*EDOTU*U(1) | |
517 | POLAR(2) = -POLAR(2) + 2*EDOTU*U(2) | |
518 | POLAR(3) = -POLAR(3) + 2*EDOTU*U(3) | |
519 | INWVOL = 0 | |
520 | ELSE | |
521 | * | |
522 | * *** Calculate amplitude for transmission (Q = P x U) | |
523 | * | |
524 | COST2 = SIGN(1.0,COST1)*SQRT((1-SINT2)*(1+SINT2)) | |
525 | IF (SINT1.GT.1.E-4) THEN | |
526 | QQ(1) = VECT(5)*U(3) - VECT(6)*U(2) | |
527 | QQ(2) = VECT(6)*U(1) - VECT(4)*U(3) | |
528 | QQ(3) = VECT(4)*U(2) - VECT(5)*U(1) | |
529 | EPERP1 = (POLAR(1)*QQ(1) + POLAR(2)*QQ(2) + POLAR(3)* | |
530 | + QQ(3)) | |
531 | ENORM = SQRT(EPERP1**2+EDOTU**2) | |
532 | EPERP1 = EPERP1/ENORM | |
533 | EPARL1 = EDOTU/ENORM | |
534 | ELSE | |
535 | QQ(1) = POLAR(1) | |
536 | QQ(2) = POLAR(2) | |
537 | QQ(3) = POLAR(3) | |
538 | * | |
539 | * Here we follow Jackson's conventions and we set the parallel | |
540 | * component = 1 in case of a ray perpendicular to the surface | |
541 | EPERP1 = 0. | |
542 | EPARL1 = 1. | |
543 | END IF | |
544 | IF(COST1.NE.0.) THEN | |
545 | S1 = RIN1*COST1 | |
546 | EPERP2 = 2*S1*EPERP1/(RIN1*COST1+RIN2*COST2) | |
547 | EPARL2 = 2*S1*EPARL1/(RIN2*COST1+RIN1*COST2) | |
548 | E2 = EPERP2**2 + EPARL2**2 | |
549 | S2 = RIN2*COST2*E2 | |
550 | TCOEF = S2/S1 | |
551 | ELSE | |
552 | TCOEF = 0. | |
553 | ENDIF | |
554 | CALL GRNDM(RNDM,1) | |
555 | IF (RNDM.GT.TCOEF) THEN | |
556 | * | |
557 | * *** Simulate reflection | |
558 | * | |
559 | NMEC=NMEC+1 | |
560 | LMEC(NMEC)=106 | |
561 | * | |
562 | * *** Restore old volume tree, the photon does not cross the boundary | |
563 | * | |
564 | * CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR) | |
565 | CALL GTMEDI(VIN,N) | |
566 | LOLDTR=.FALSE. | |
567 | * | |
568 | * *** Reflect, but make sure we are in the old volume before | |
569 | * | |
570 | NBPUSH = 0 | |
571 | 100 NBPUSH = NBPUSH+1 | |
572 | IF(NBPUSH.GT.MXPUSH) THEN | |
573 | WRITE(CHMAIL,10300) NTMULT,NSTEP | |
574 | CALL GMAIL(0,0) | |
575 | ISTOP=1 | |
576 | GOTO 110 | |
577 | ELSE | |
578 | CALL GINVOL(VECT,ISAME) | |
579 | IF(ISAME.EQ.0) THEN | |
580 | PRECN = NBPUSH*PREC | |
581 | VECT(1) = VECT(1) - PRECN*VECT(4) | |
582 | VECT(2) = VECT(2) - PRECN*VECT(5) | |
583 | VECT(3) = VECT(3) - PRECN*VECT(6) | |
584 | GO TO 100 | |
585 | ENDIF | |
586 | ENDIF | |
587 | * | |
588 | VECT(4) = VECT(4) - 2*PDOTU*U(1) | |
589 | VECT(5) = VECT(5) - 2*PDOTU*U(2) | |
590 | VECT(6) = VECT(6) - 2*PDOTU*U(3) | |
591 | IF(SINT1.GT.1E-4) THEN | |
592 | EPARL2 = RIN2*EPARL2/RIN1-EPARL1 | |
593 | EPERP2 = EPERP2-EPERP1 | |
594 | E2 = EPERP2**2 + EPARL2**2 | |
595 | R(1) = U(1) + PDOTU*VECT(4) | |
596 | R(2) = U(2) + PDOTU*VECT(5) | |
597 | R(3) = U(3) + PDOTU*VECT(6) | |
598 | EABS = SQRT(E2)*SINT1 | |
599 | CPARL = EPARL2/EABS | |
600 | CPERP = EPERP2/EABS | |
601 | POLAR(1) = CPARL*R(1) - CPERP*QQ(1) | |
602 | POLAR(2) = CPARL*R(2) - CPERP*QQ(2) | |
603 | POLAR(3) = CPARL*R(3) - CPERP*QQ(3) | |
604 | ELSEIF(RIN2.GT.RIN1) THEN | |
605 | * | |
606 | * *** Case of ray perpendicular to the surface. No change or | |
607 | * *** an inversion of phase. | |
608 | POLAR(1) = -POLAR(1) | |
609 | POLAR(2) = -POLAR(2) | |
610 | POLAR(3) = -POLAR(3) | |
611 | ENDIF | |
612 | INWVOL = 0 | |
613 | ELSE | |
614 | * | |
615 | * *** Simulate transmission/refraction | |
616 | * | |
617 | NMEC=NMEC+1 | |
618 | LMEC(NMEC)=107 | |
619 | VECT(1) = VOUT(1) | |
620 | VECT(2) = VOUT(2) | |
621 | VECT(3) = VOUT(3) | |
622 | IF(SINT1.GT.1E-4) THEN | |
623 | ALPHA = COST1-COST2*(RIN2/RIN1) | |
624 | D(1) = VECT(4) + ALPHA*U(1) | |
625 | D(2) = VECT(5) + ALPHA*U(2) | |
626 | D(3) = VECT(6) + ALPHA*U(3) | |
627 | DABS = SQRT(D(1)**2+D(2)**2+D(3)**2) | |
628 | VECT(4) = D(1)/DABS | |
629 | VECT(5) = D(2)/DABS | |
630 | VECT(6) = D(3)/DABS | |
631 | PDOTU = -COST2 | |
632 | R(1) = U(1) - PDOTU*VECT(4) | |
633 | R(2) = U(2) - PDOTU*VECT(5) | |
634 | R(3) = U(3) - PDOTU*VECT(6) | |
635 | EABS = SQRT(E2) | |
636 | CPARL = EPARL2/(EABS*SINT2) | |
637 | CPERP = EPERP2/(EABS*SINT1) | |
638 | POLAR(1) = CPARL*R(1) + CPERP*QQ(1) | |
639 | POLAR(2) = CPARL*R(2) + CPERP*QQ(2) | |
640 | POLAR(3) = CPARL*R(3) + CPERP*QQ(3) | |
641 | ENDIF | |
642 | GEKRAT = GEKRT2 | |
643 | IEKBIN = IEKBI2 | |
644 | STEPLA = ABSCO2 | |
645 | EFFIC = EFFIC2 | |
646 | RIN1 = RIN2 | |
647 | END IF | |
648 | END IF | |
649 | END IF | |
650 | ENDIF | |
651 | * END GTCKOV | |
652 | 110 IF(LOLDTR) CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR) | |
653 | 10200 FORMAT(' **** GTCKOV: error from GLISUR = ',I10) | |
654 | 10300 FORMAT(' **** GTCKOV: unable to reflect at NTMULT = ', | |
655 | + I8,' step No. ',I8,' photon abandoned!') | |
656 | END |