]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/gtrak/gtckov.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gtrak / gtckov.F
CommitLineData
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
8C.
9C. ******************************************************************
10C. * *
11C. * This routine is called to follow the Cherenkov photons *
12C. * created during the tracking of charged particles and *
13C. * simulate the relevant processes along the way, until either *
14C. * the photon is absorbed or exits the detector. Processes *
15C. * currently simulated are absorption in-flight, and reflection *
16C. * /transmission/absorption at a medium boundary. There are two *
17C. * boundary types: dielectric-metal and dielectric-dielectric. *
18C. * For each of these there is a continuum of reflectivity *
19C. * and of surface quality from mirror finish to matte. The *
20C. * surface model is contained in routine GHSURF. *
21C. * *
22C. * ==>Called by : GTRACK *
23C. * Authors F.Carminati, R.Jones ************ *
24C. * *
25C. ******************************************************************
26C.
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
62C.
63C. ------------------------------------------------------------------
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
187C
188C In case of Cherenkovs beeing near a corner, particle holding is
189C prevented by a bigger step size. The value of VBOU is only used to
190C calculate the correct surface normal (by GLISUR and GGPERP). The
191C tracking of the Cherenkov is still using STEP
192C
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
198C
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)
25710000 FORMAT(' Boundary correction in GTCKOV: ',
258 + ' GEKIN NUMED STEP SNEXT')
25910100 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)
65310200 FORMAT(' **** GTCKOV: error from GLISUR = ',I10)
65410300 FORMAT(' **** GTCKOV: unable to reflect at NTMULT = ',
655 + I8,' step No. ',I8,' photon abandoned!')
656 END