]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TAmpt/AMPT/zpc.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
randomize reaction plane
[u/mrichter/AliRoot.git] / TAmpt / AMPT / zpc.f
CommitLineData
0119ef9a 1c.................... zpc.f
2c PROGRAM ZPC
3 SUBROUTINE ZPCMN
4c Version: 1.0.1
5c Author: Bin Zhang
6c (suggestions, problems -> bzhang@nt1.phys.columbia.edu)
7cms
8cms dlw & gsfs Comments our writing of output files
9cms
10 implicit double precision (a-h, o-z)
11clin-4/20/01 PARAMETER (NMAXGL = 16000)
12 parameter (MAXPTN=400001)
13 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
14cc SAVE /para3/
15 SAVE
16c
17c loop over events
18 do 1000 i = 1, nevnt
19 ievt = i
20c generation of the initial condition for one event
21 call inievt
22c loop over many runs of the same event
23 do 2000 j = 1, nsbrun
24 isbrun = j
25c initialization for one run of an event
26 call inirun
27clin-4/2008 not used:
28c CALL HJAN1A
29 3000 continue
30c do one collision
31 call zpcrun(*4000)
32 call zpca1
33 goto 3000
34 4000 continue
35 call zpca2
36 2000 continue
37 1000 continue
38 call zpcou
39clin-5/2009 ctest off
40c 5/17/01 calculate v2 for parton already frozen out:
41c call flowp(3)
42c.....to get average values for different strings
43 CALL zpstrg
44 RETURN
45 end
46
47******************************************************************************
48******************************************************************************
49
50 block data zpcbdt
51c set initial values in block data
52
53 implicit double precision (a-h, o-z)
54 parameter (MAXPTN=400001)
55 PARAMETER (MAXSTR=150001)
56 common /para1/ mul
57cc SAVE /para1/
58 common /para2/ xmp, xmu, alpha, rscut2, cutof2
59cc SAVE /para2/
60 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
61cc SAVE /para3/
62 common /para4/ iftflg, ireflg, igeflg, ibstfg
63cc SAVE /para4/
64 common /para5/ iconfg, iordsc
65cc SAVE /para5/
66 common /para6/ centy
67cc SAVE /para6/
68clin-6/2009 nsmbbbar and nsmmeson respectively give the total number of
69c baryons/anti-baryons and mesons for each event:
70c common /para7/ ioscar
71 common /para7/ ioscar,nsmbbbar,nsmmeson
72cc SAVE /para7/
73 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
74 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
75 & XMASS0(MAXPTN), ITYP0(MAXPTN)
76cc SAVE /prec1/
77 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
78 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
79 & xmass(MAXPTN), ityp(MAXPTN)
80cc SAVE /prec2/
81 common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
82 & pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
83 & xmasss(MAXPTN), ityps(MAXPTN)
84cc SAVE /prec3/
85 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
86cc SAVE /prec4/
87 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
88cc SAVE /prec5/
89 common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
90cc SAVE /prec6/
91 common /aurec1/ jxa, jya, jza
92cc SAVE /aurec1/
93 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
94cc SAVE /aurec2/
95 common /ilist1/
96 & iscat, jscat, next(MAXPTN), last(MAXPTN),
97 & ictype, icsta(MAXPTN),
98 & nic(MAXPTN), icels(MAXPTN)
99cc SAVE /ilist1/
100 common /ilist2/ icell, icel(10,10,10)
101cc SAVE /ilist2/
102 common /ilist3/ size1, size2, size3, v1, v2, v3, size
103cc SAVE /ilist3/
104 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
105cc SAVE /ilist4/
106c 6/07/02 initialize in ftime to expedite compiling:
107c common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
108cc SAVE /ilist5/
109 common /ilist6/ t, iopern, icolln
110cc SAVE /ilist6/
111 COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
112cc SAVE /ilist7/
113 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
114cc SAVE /ilist8/
115 common /rndm1/ number
116cc SAVE /rndm1/
117 common /rndm2/ iff
118cc SAVE /rndm2/
119 common /rndm3/ iseedp
120cc SAVE /rndm3/
121 common /ana1/ ts(12)
122cc SAVE /ana1/
123 common /ana2/
124 & det(12), dn(12), detdy(12), detdn(12), dndy(12),
125 & det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
126 & det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
127cc SAVE /ana2/
128 common /ana3/ em(4, 4, 12)
129cc SAVE /ana3/
130 common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
131cc SAVE /ana4/
132 SAVE
133 data centy/0d0/
134c 6/07/02 initialize in ftime to expedite compiling:
135c data (ct(i), i = 1, MAXPTN)/MAXPTN*0d0/
136c data (ot(i), i = 1, MAXPTN)/MAXPTN*0d0/
137c data tlarge/1000000.d0/
138 data number/0/
139 data ts/0.11d0, 0.12d0, 0.15d0, 0.2d0, 0.3d0, 0.4d0, 0.6d0,
140 & 0.8d0, 1d0, 2d0, 4d0, 6d0/
141c
142 end
143
144******************************************************************************
145******************************************************************************
146
147 subroutine inizpc
148
149 implicit double precision (a-h, o-z)
150 SAVE
151
152 call readpa
153
154 call inipar
155
156 call inian1
157
158 return
159 end
160
161 subroutine readpa
162
163 implicit double precision (a-h, o-z)
164
165cc external ran1
166
167 character*50 str
168
169 common /para2/ xmp, xmu, alpha, rscut2, cutof2
170cc SAVE /para2/
171 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
172cc SAVE /para3/
173 common /para4/ iftflg, ireflg, igeflg, ibstfg
174cc SAVE /para4/
175 common /para5/ iconfg, iordsc
176cc SAVE /para5/
177 common /para7/ ioscar,nsmbbbar,nsmmeson
178cc SAVE /para7/
179 common /ilist3/ size1, size2, size3, v1, v2, v3, size
180cc SAVE /ilist3/
181 common /rndm1/ number
182cc SAVE /rndm1/
183 common /rndm2/ iff
184cc SAVE /rndm2/
185 common /rndm3/ iseedp
186cc SAVE /rndm3/
187 SAVE
188
189 str=str
190 iseed=iseedp
191c this is the initialization file containing the initial values of
192c the parameters
193cbz1/31/99
194c open (5, file = 'zpc.ini', status = 'unknown')
195cbz1/31/99end
196
197c this is the final data file containing general info about the cascade
198cbz1/31/99
199c open (6, file = 'zpc.res', status = 'unknown')
200cms open (25, file = 'ana/zpc.res', status = 'unknown')
201cbz1/31/99end
202
203c this is the input file containing initial particle records
204cbz1/25/99
205c open (7, file = 'zpc.inp', status = 'unknown')
206cbz1/25/99end
207
208c this gives the optional OSCAR standard output
209cbz1/31/99
210c open (8, file = 'zpc.oscar', status = 'unknown')
211 if(ioscar.eq.1) then
212cms open (26, file = 'ana/parton.oscar', status = 'unknown')
213cms open (19, file = 'ana/hadron.oscar', status = 'unknown')
214 endif
215cbz1/31/99end
216
217c 2/11/03 combine zpc initialization into ampt.ini:
218c open (29, file = 'zpc.ini', status = 'unknown')
219c read (29, *) str, xmp
220 xmp=0d0
221c read (29, *) str, xmu
222c read (29, *) str, alpha
223 cutof2 = 4.5d0 * (alpha / xmu) ** 2
224c read (29, *) str, rscut2
225 rscut2=0.01d0
226c read (29, *) str, nsevt
227 nsevt=1
228c read (29, *) str, nevnt
229 nevnt=1
230c read (29, *) str, nsbrun
231 nsbrun=1
232c read (29, *) str, iftflg
233 iftflg=0
234c read (29, *) str, ireflg
235 ireflg=1
236cbz1/31/99
237 IF (ireflg .EQ. 0) THEN
238cms OPEN (27, FILE = 'zpc.inp', STATUS = 'UNKNOWN')
239 END IF
240cbz1/31/99end
241c read (29, *) str, igeflg
242 igeflg=0
243c read (29, *) str, ibstfg
244 ibstfg=0
245c read (29, *) str, iconfg
246 iconfg=1
247c read (29, *) str, iordsc
248 iordsc=11
249c read (29, *) str, ioscar
250c read (29, *) str, v1, v2, v3
251 v1=0.2d0
252 v2=0.2d0
253 v3=0.2d0
254c read (29, *) str, size1, size2, size3
255 size1=1.5d0
256 size2=1.5d0
257 size3=0.7d0
258 if (size1 .eq. 0d0 .or. size2 .eq. 0d0 .or.
259 & size3 .eq. 0d0) then
260 if (size1 .ne. 0d0 .or. size2 .ne. 0d0 .or. size3 .ne. 0d0
261 & .or. v1 .ne. 0d0 .or. v2 .ne. 0d0 .or. v3 .ne. 0d0) then
262 print *, 'to get rid of space division:'
263 print *, 'set all sizes and vs to 0'
264 stop 'chker'
265 end if
266 end if
267 size = min(size1, size2, size3)
268c read (29, *) str, iff
269 iff=-1
270c read (29, *) str, iseed
271
272c 10/24/02 get rid of argument usage mismatch in ran1():
273 isedng=-iseed
274c a = ran1(-iseed)
275 a = ran1(isedng)
276c read (29, *) str, irused
277 irused=2
278 do 1001 i = 1, irused - 1
279c a = ran1(2)
280 iseed2=2
281 a = ran1(iseed2)
282 1001 continue
283c 10/24/02-end
284
285 if (iconfg .eq. 2 .or. iconfg .eq. 3) then
286 v1 = 0d0
287 v2 = 0d0
288 end if
289
290 if (iconfg .eq. 4 .or. iconfg .eq. 5) then
291 v1 = 0d0
292 v2 = 0d0
293 v3 = 0d0
294 end if
295
296c close(5)
297
298 return
299 end
300
301 subroutine inipar
302
303 implicit double precision (a-h,o-z)
304
305 common /para4/ iftflg, ireflg, igeflg, ibstfg
306cc SAVE /para4/
307 common /para6/ centy
308cc SAVE /para6/
309 SAVE
310
311 if (ibstfg .ne. 0) then
312 centy = -6d0
313 end if
314
315 return
316 end
317
318 subroutine inian1
319
320 implicit double precision (a-h,o-z)
321
322 common /para4/ iftflg, ireflg, igeflg, ibstfg
323cc SAVE /para4/
324 common /ana1/ ts(12)
325cc SAVE /ana1/
326 SAVE
327 if (ibstfg .ne. 0) then
328 a = cosh(6d0)
329 do 1001 i = 1, 12
330 ts(i) = ts(i) * a
331 1001 continue
332 end if
333
334 return
335 end
336
337******************************************************************************
338
339 subroutine inievt
340
341 implicit double precision (a-h, o-z)
342
343 COMMON /para1/ mul
344cc SAVE /para1/
345 common /para4/ iftflg, ireflg, igeflg, ibstfg
346cc SAVE /para4/
347 SAVE
348
349cbz1/25/99
350c mul = 0
351cbz1/25/99
352 if (ireflg .eq. 0) call readi
353 if (igeflg .ne. 0) call genei
354 if (ibstfg .ne. 0) call boosti
355
356 return
357 end
358
359 subroutine readi
360
361 implicit double precision (a-h, o-z)
362 parameter (MAXPTN=400001)
363 double precision field(9)
364 common /para1/ mul
365cc SAVE /para1/
366 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
367cc SAVE /para3/
368 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
369 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
370 & XMASS0(MAXPTN), ITYP0(MAXPTN)
371cc SAVE /prec1/
372 SAVE
373 do 1001 i = 1, MAXPTN
374 if (ievt .ne. 1 .and. i .eq. 1) then
375 ityp0(i) = ntyp
376 gx0(1) = field(1)
377 gy0(1) = field(2)
378 gz0(1) = field(3)
379 ft0(1) = field(4)
380 px0(1) = field(5)
381 py0(1) = field(6)
382 pz0(1) = field(7)
383 e0(1) = field(8)
384 xmass0(i) = field(9)
385 mul = 1
386 else
387 900 read (27, *, end = 1000) neve, ntyp, field
388 if (neve .lt. nsevt) goto 900
389 if (neve .gt.
390 & nsevt + ievt - 1) goto 1000
391 ityp0(i) = ntyp
392 gx0(i) = field(1)
393 gy0(i) = field(2)
394 gz0(i) = field(3)
395 ft0(i) = field(4)
396 px0(i) = field(5)
397 py0(i) = field(6)
398 pz0(i) = field(7)
399 e0(i) = field(8)
400 xmass0(i) = field(9)
401 mul = mul + 1
402 end if
403 1001 continue
404
405 1000 continue
406
407 return
408 end
409
410 subroutine genei
411
412 implicit double precision (a-h, o-z)
413 parameter (MAXPTN=400001)
414 common /para1/ mul
415cc SAVE /para1/
416 common /para2/ xmp, xmu, alpha, rscut2, cutof2
417cc SAVE /para2/
418 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
419cc SAVE /para3/
420 common /para5/ iconfg, iordsc
421cc SAVE /para5/
422 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
423 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
424 & XMASS0(MAXPTN), ITYP0(MAXPTN)
425cc SAVE /prec1/
426 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
427cc SAVE /prec5/
3006c44b 428 common /lora/ enenew, pxnew, pynew, pznew
429cc SAVE /lora/
0119ef9a 430 common /rndm3/ iseedp
431cc SAVE /rndm3/
432 SAVE
433
434cc external ran1
435
436 iseed=iseedp
437 incmul = 4000
438 temp = 0.5d0
439 etamin = -5d0
440 etamax = 5d0
441 r0 = 5d0
442 tau0 = 0.1d0
443 deta = etamax - etamin
444
445 do 1001 i = mul + 1, mul + incmul
446 ityp0(i) = 21
447 xmass0(i) = xmp
e2670fbe 448 call energya(e, temp)
0119ef9a 449 call momntm(px, py, pz, e)
450c 7/20/01:
451c e = sqrt(e ** 2 + xmp ** 2)
452 e = dsqrt(e ** 2 + xmp ** 2)
453 if (iconfg .le. 3) then
454 eta(i) = etamin + deta * ran1(iseed)
455 bex = 0d0
456 bey = 0d0
457 bez = -tanh(eta(i))
e2670fbe 458 call lorenza(e, px, py, pz, bex, bey, bez)
0119ef9a 459 px0(i) = pxnew
460 py0(i) = pynew
461 pz0(i) = pznew
462 e0(i) = enenew
463 else
464 px0(i) = px
465 py0(i) = py
466 pz0(i) = pz
467 e0(i) = e
468 end if
469 1001 continue
470
471 do 1002 i = mul + 1, mul + incmul
472 if (iconfg .le. 3) then
473 gz0(i) = tau0 * sinh(eta(i))
474 ft0(i) = tau0 * cosh(eta(i))
475 if (iconfg .eq. 1) then
476 call posit1(x, y, r0)
477 gx0(i) = x + px0(i) * ft0(i)/e0(i)
478 gy0(i) = y + py0(i) * ft0(i)/e0(i)
479 else if (iconfg .eq. 2 .or. iconfg .eq. 3) then
480 call posit2(x, y)
481 gx0(i) = x
482 gy0(i) = y
483 end if
484 else
485 ft0(i) = 0d0
486 call posit3(x, y, z)
487 gx0(i) = x
488 gy0(i) = y
489 gz0(i) = z
490 end if
491 1002 continue
492
493 mul = mul + incmul
494
495c check if it's necessary to adjust array size 'adarr'
496 if (mul .ge. MAXPTN .or. mul .eq. 0) then
497 print *, 'event',ievt,'has',mul,'number of gluon',
498 & 'adjusting counting is necessary'
499 stop 'adarr'
500 end if
501
502 return
503 end
504
505 subroutine posit1(x, y, r0)
506
507 implicit double precision (a-h, o-z)
508
509cc external ran1
510 common /rndm3/ iseedp
511cc SAVE /rndm3/
512 SAVE
513
514 iseed=iseedp
515 10 x = 2d0 * ran1(iseed) - 1d0
516 y = 2d0 * ran1(iseed) - 1d0
517 if (x ** 2 + y ** 2 .gt. 1d0) goto 10
518 x = x * r0
519 y = y * r0
520
521 return
522 end
523
524 subroutine posit2(x, y)
525
526 implicit double precision (a-h, o-z)
527
528c external ran1
529
530 common /ilist3/ size1, size2, size3, v1, v2, v3, size
531cc SAVE /ilist3/
532 common /rndm3/ iseedp
533cc SAVE /rndm3/
534 SAVE
535 iseed=iseedp
536 x = 2d0 * ran1(iseed) - 1d0
537 y = 2d0 * ran1(iseed) - 1d0
538 x = x * 5d0 * size1
539 y = y * 5d0 * size2
540
541 return
542 end
543
544 subroutine posit3(x, y, z)
545
546 implicit double precision (a-h, o-z)
547
548cc external ran1
549
550 common /ilist3/ size1, size2, size3, v1, v2, v3, size
551cc SAVE /ilist3/
552 common /rndm3/ iseedp
553cc SAVE /rndm3/
554 SAVE
555
556 iseed=iseedp
557 x = 2d0 * ran1(iseed) - 1d0
558 y = 2d0 * ran1(iseed) - 1d0
559 z = 2d0 * ran1(iseed) - 1d0
560 x = x * 5d0 * size1
561 y = y * 5d0 * size2
562 z = z * 5d0 * size3
563
564 return
565 end
566
e2670fbe 567 subroutine energya(e, temp)
0119ef9a 568
569c to generate the magnitude of the momentum e,
570c knowing the temperature of the local thermal distribution temp
571
572 implicit double precision (a-h, o-z)
573
574cc external ran1
575
576 common /para2/ xmp, xmu, alpha, rscut2, cutof2
577cc SAVE /para2/
578 common /rndm3/ iseedp
579cc SAVE /rndm3/
580 SAVE
581
582 iseed=iseedp
583 1000 continue
584
585 e = ran1(iseed)
586 e = e * ran1(iseed)
587 e = e * ran1(iseed)
588
589 if (e .le. 0d0) goto 1000
590 e = - temp * log(e)
591 if (ran1(iseed) .gt.
592 & exp((e - dsqrt(e ** 2 + xmp ** 2))/temp)) then
593 goto 1000
594 end if
595
596 return
597 end
598
599 subroutine momntm(px, py, pz, e)
600
601c to generate the 3 components of the momentum px, py, pz,
602c from the magnitude of the momentum e
603
604 implicit double precision (a-h,o-z)
605
606cc external ran1
607
608 parameter (pi = 3.14159265358979d0)
609 common /rndm3/ iseedp
610cc SAVE /rndm3/
611 SAVE
612
613 iseed=iseedp
614 cost = 2d0 * ran1(iseed) - 1d0
615c 7/20/01:
616c sint = sqrt(1d0 - cost ** 2)
617 sint = dsqrt(1d0 - cost ** 2)
618 phi = 2d0 * pi * ran1(iseed)
619
620 px = e * sint * cos(phi)
621 py = e * sint * sin(phi)
622 pz = e * cost
623
624 return
625 end
626
627 subroutine boosti
628
629 implicit double precision (a-h, o-z)
630 parameter (MAXPTN=400001)
631 common /para1/ mul
632cc SAVE /para1/
633 common /para6/ centy
634cc SAVE /para6/
635 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
636 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
637 & XMASS0(MAXPTN), ITYP0(MAXPTN)
638cc SAVE /prec1/
3006c44b 639 common /lora/ enenew, pxnew, pynew, pznew
640cc SAVE /lora/
0119ef9a 641 SAVE
642
e2670fbe 643 external lorenza
0119ef9a 644
645 bex = 0d0
646 bey = 0d0
647 bez = - tanh(centy)
648
649c save data for many runs of the same initial condition
650 do 1001 i = 1, mul
651 px1 = gx0(i)
652 py1 = gy0(i)
653 pz1 = gz0(i)
654 e1 = ft0(i)
e2670fbe 655 call lorenza(e1, px1, py1, pz1, bex, bey, bez)
0119ef9a 656 gx0(i) = pxnew
657 gy0(i) = pynew
658 gz0(i) = pznew
659 ft0(i) = enenew
660 px1 = px0(i)
661 py1 = py0(i)
662 pz1 = pz0(i)
663 e1 = e0(i)
e2670fbe 664 call lorenza(e1, px1, py1, pz1, bex, bey, bez)
0119ef9a 665 px0(i) = pxnew
666 py0(i) = pynew
667 pz0(i) = pznew
668 e0(i) = enenew
669 1001 continue
670
671 return
672 end
673
674******************************************************************************
675
676 subroutine inirun
677 SAVE
678
679c sort prec2 according to increasing formation time
680 call ftime
681 call inirec
682 call iilist
683 call inian2
684
685 return
686 end
687
688 subroutine ftime
689c this subroutine generates formation time for the particles
690c indexing ft(i)
691c input e(i)
692c output ft(i), indx(i)
693
694 implicit double precision (a-h, o-z)
695
696 external ftime1
697 parameter (MAXPTN=400001)
698 common /para1/ mul
699cc SAVE /para1/
700 common /para2/ xmp, xmu, alpha, rscut2, cutof2
701cc SAVE /para2/
702 common /para4/ iftflg, ireflg, igeflg, ibstfg
703cc SAVE /para4/
704 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
705 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
706 & XMASS0(MAXPTN), ITYP0(MAXPTN)
707cc SAVE /prec1/
708 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
709cc SAVE /prec4/
710 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
711cc SAVE /ilist4/
712 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
713cc SAVE /ilist5/
714 common /par1/ formt
715cc SAVE /par1/
716 common/anim/nevent,isoft,isflag,izpc
717cc SAVE /anim/
718 common /rndm3/ iseedp
719cc SAVE /rndm3/
720 SAVE
721
722 iseed=iseedp
723clin-6/07/02 initialize here to expedite compiling, instead in zpcbdt:
724 do 1001 i = 1, MAXPTN
725 ct(i)=0d0
726 ot(i)=0d0
727 1001 continue
728 tlarge=1000000.d0
729clin-6/07/02-end
730
731 if (iftflg .eq. 0) then
732c 5/01/01 different prescription for parton initial formation time:
733 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
734 do 1002 i = 1, mul
735 if (ft0(i) .gt. tlarge) ft0(i) = tlarge
736 1002 continue
737 goto 150
738 else
739c 5/01/01-end
740
741 do 1003 i = 1, MAXPTN
742 ft0(i) = tlarge
743 1003 continue
744 do 1004 i = 1, mul
745 xmt2 = px0(i) ** 2 + py0(i) ** 2 + xmp ** 2
746 formt = xmt2 / e0(i)
747 ft0(i) = ftime1(iseed)
748 if (ft0(i) .gt. tlarge) ft0(i) = tlarge
749 1004 continue
750c 5/01/01:
751 endif
752
753 end if
754
755c 5/01/01:
756 150 continue
757
758c call index1(MAXPTN, mul, ft0, indx)
759 if (mul .gt. 1) then
760 call index1(MAXPTN, mul, ft0, indx)
761 else
762clin-7/09/03: need to set value for mul=1:
763 indx(1)=1
764 end if
765c
766 return
767 end
768
769 subroutine inirec
770
771 implicit double precision (a-h, o-z)
772cc external ran1
773 parameter (MAXPTN=400001)
774 common /para1/ mul
775cc SAVE /para1/
776 common /para4/ iftflg, ireflg, igeflg, ibstfg
777cc SAVE /para4/
778 common /para5/ iconfg, iordsc
779cc SAVE /para5/
780 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
781 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
782 & XMASS0(MAXPTN), ITYP0(MAXPTN)
783cc SAVE /prec1/
784 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
785 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
786 & xmass(MAXPTN), ityp(MAXPTN)
787cc SAVE /prec2/
788 common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
789 & pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
790 & xmasss(MAXPTN), ityps(MAXPTN)
791cc SAVE /prec3/
792 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
793cc SAVE /prec4/
794 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
795cc SAVE /prec5/
796 common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
797cc SAVE /prec6/
798 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
799cc SAVE /ilist4/
800cbz1/25/99
801 COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
802cc SAVE /ilist7/
803 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
804cc SAVE /ilist8/
805cbz1/25/99end
806 COMMON /smearz/smearp,smearh
807cc SAVE /smearz/
808 dimension vxp(MAXPTN), vyp(MAXPTN), vzp(MAXPTN)
809 common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN)
810cc SAVE /precpa/
811 common/anim/nevent,isoft,isflag,izpc
812cc SAVE /anim/
813clin-6/06/02 local parton freezeout:
814 common /frzprc/
815 & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
816 & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
817 & xmfrz(MAXPTN),
818 & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
819cc SAVE /frzprc/
820 common /rndm3/ iseedp
821cc SAVE /rndm3/
822 common /para7/ ioscar,nsmbbbar,nsmmeson
823 COMMON /AREVT/ IAEVT, IARUN, MISS
824 SAVE
825 iseed=iseedp
826clin-6/06/02 local freezeout initialization:
827 if(isoft.eq.5) then
828 itlast=0
829 call inifrz
830 endif
831
832 do 1001 i = 1, mul
833clin-7/09/01 define indx(i) to save time:
834c ityp(i) = ityp0(indx(i))
835c gx(i) = gx0(indx(i))
836c gy(i) = gy0(indx(i))
837c gz(i) = gz0(indx(i))
838c ft(i) = ft0(indx(i))
839c px(i) = px0(indx(i))
840c py(i) = py0(indx(i))
841c pz(i) = pz0(indx(i))
842c e(i) = e0(indx(i))
843c xmass(i) = xmass0(indx(i))
844ccbz1/25/99
845c LSTRG1(I) = LSTRG0(INDX(I))
846c LPART1(I) = LPART0(INDX(I))
847ccbz1/25/99end
848 indxi=indx(i)
849 ityp(i) = ityp0(indxi)
850 gx(i) = gx0(indxi)
851 gy(i) = gy0(indxi)
852 gz(i) = gz0(indxi)
853 ft(i) = ft0(indxi)
854 px(i) = px0(indxi)
855 py(i) = py0(indxi)
856 pz(i) = pz0(indxi)
857 e(i) = e0(indxi)
858 xmass(i) = xmass0(indxi)
859 LSTRG1(I) = LSTRG0(INDXI)
860 LPART1(I) = LPART0(INDXI)
861 vxp(I) = vxp0(INDXI)
862 vyp(I) = vyp0(INDXI)
863 vzp(I) = vzp0(INDXI)
864clin-7/09/01-end
865c
866clin-6/06/02 local freezeout initialization:
867 if(isoft.eq.5) then
868 idfrz(i)=ityp(i)
869 gxfrz(i)=gx(i)
870 gyfrz(i)=gy(i)
871 gzfrz(i)=gz(i)
872 ftfrz(i)=ft(i)
873 pxfrz(i)=px(i)
874 pyfrz(i)=py(i)
875 pzfrz(i)=pz(i)
876 efrz(i)=e(i)
877 xmfrz(i)=xmass(i)
878 ifrz(i)=0
879 endif
880clin-6/06/02-end
881 1001 continue
882
883c save particle info for fixed time analysis
884 do 1002 i = 1, mul
885 ityps(i) = ityp(i)
886 gxs(i) = gx(i)
887 gys(i) = gy(i)
888 gzs(i) = gz(i)
889 fts(i) = ft(i)
890 pxs(i) = px(i)
891 pys(i) = py(i)
892 pzs(i) = pz(i)
893 es(i) = e(i)
894 xmasss(i) = xmass(i)
895 1002 continue
896
897clin-6/2009
898cms if(isoft.eq.1.and.(ioscar.eq.2.or.ioscar.eq.3))
899cms 1 write(92,*) iaevt,miss,mul
900
901 do 1003 i = 1, mul
902 energy = e(i)
903 vx(i) = px(i) / energy
904 vy(i) = py(i) / energy
905 vz(i) = pz(i) / energy
906 if (iftflg .eq. 0) then
907 formt = ft(i)
908c 7/09/01 propagate partons with parent velocity till formation
909c so that partons in same hadron have 0 distance:
910c gx(i) = gx(i) + vx(i) * formt
911c gy(i) = gy(i) + vy(i) * formt
912c gz(i) = gz(i) + vz(i) * formt
913 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
914 gx(i) = gx(i) + vxp(i) * formt
915 gy(i) = gy(i) + vyp(i) * formt
916 gz(i) = gz(i) + vzp(i) * formt
917 else
918 gx(i) = gx(i) + vx(i) * formt
919 gy(i) = gy(i) + vy(i) * formt
920 gz(i) = gz(i) + vz(i) * formt
921 endif
922c 7/09/01-end
923c
924c 3/27/00-ctest off no smear z on partons to avoid eta overflow:
925c gz(i) = gz(i)+smearp*(2d0 * ran1(iseed) - 1d0)
926c to give eta=y +- smearp*random:
927c smeary=smearp*(2d0 * ran1(iseed) - 1d0)
928c smearf=dexp(2*smeary)*(1+vz(i))/(1-vz(i)+1.d-8)
929c gz(i) = gz(i)+formt*(smearf-1)/(smearf+1)
930c 3/27/00-end
931 end if
932
933clin-6/2009 write out initial parton information after string melting
934c and after propagating to its format time:
935 if(ioscar.eq.2.or.ioscar.eq.3) then
936 if(dmax1(abs(gx(i)),abs(gy(i)),
937 1 abs(gz(i)),abs(ft(i))).lt.9999) then
938cms write(92,200) ityp(i),px(i),py(i),pz(i),xmass(i),
939cms 1 gx(i),gy(i),gz(i),ft(i)
940 else
941cms write(92,201) ityp(i),px(i),py(i),pz(i),xmass(i),
942cms 1 gx(i),gy(i),gz(i),ft(i)
943 endif
944 endif
945cyy 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
946cyy 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
947c
948 1003 continue
949
950 if (iconfg .le. 3) then
951 do 1004 i = 1, mul
952 if (ft(i) .le. abs(gz(i))) then
953 eta(i) = 1000000.d0
954 else
955 eta(i) = 0.5d0 * log((ft(i) + gz(i)) / (ft(i) - gz(i)))
956 end if
957 if (e(i) .le. abs(pz(i))) then
958 rap(i) = 1000000.d0
959 else
960 rap(i) = 0.5d0 * log((e(i) + pz(i)) / (e(i) - pz(i)))
961 end if
962 tau(i) = ft(i) / cosh(eta(i))
963 1004 continue
964
965 do 1005 i = 1, mul
966 etas(i) = eta(i)
967 raps(i) = rap(i)
968 taus(i) = tau(i)
969 1005 continue
970 end if
971
972 return
973 end
974
975 subroutine iilist
976
977 implicit double precision (a-h, o-z)
978 parameter (MAXPTN=400001)
979 common /para1/ mul
980cc SAVE /para1/
981 common /ilist1/
982 & iscat, jscat, next(MAXPTN), last(MAXPTN),
983 & ictype, icsta(MAXPTN),
984 & nic(MAXPTN), icels(MAXPTN)
985cc SAVE /ilist1/
986 common /ilist2/ icell, icel(10,10,10)
987cc SAVE /ilist2/
988 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
989cc SAVE /ilist4/
990 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
991cc SAVE /ilist5/
992 common /ilist6/ t, iopern, icolln
993cc SAVE /ilist6/
994 SAVE
995
996 iscat = MAXPTN
997 jscat = MAXPTN
998
999 do 1001 i = 1, mul
1000 next(i) = 0
1001 last(i) = 0
1002 icsta(i) = 0
1003 nic(i) = 0
1004 icels(i) = 0
1005 1001 continue
1006
1007 icell = 0
1008 do 1004 i1 = 1, 10
1009 do 1003 i2 = 1, 10
1010 do 1002 i3 = 1, 10
1011 icel(i1, i2, i3) = 0
1012 1002 continue
1013 1003 continue
1014 1004 continue
1015
1016 ichkpt = 0
1017 ifmpt = 1
1018
1019 do 1005 i = 1, mul
1020 ct(i) = tlarge
1021 ot(i) = tlarge
1022 1005 continue
1023
1024 iopern = 0
1025 icolln = 0
1026 t = 0.d0
1027
1028 return
1029 end
1030
1031 subroutine inian2
1032
1033 implicit double precision (a-h, o-z)
1034
1035 common /para5/ iconfg, iordsc
1036cc SAVE /para5/
1037 common /ana2/
1038 & det(12), dn(12), detdy(12), detdn(12), dndy(12),
1039 & det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
1040 & det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
1041cc SAVE /ana2/
1042 SAVE
1043
1044 if (iconfg .le. 3) then
1045 do 1001 i = 1, 12
1046 det(i) = 0d0
1047 dn(i) = 0d0
1048 det1(i) = 0d0
1049 dn1(i) = 0d0
1050 det2(i) = 0d0
1051 dn2(i) = 0d0
1052 1001 continue
1053 end if
1054
1055 return
1056 end
1057
1058******************************************************************************
1059******************************************************************************
1060
1061 subroutine zpcrun(*)
1062
1063 implicit double precision (a-h, o-z)
1064 parameter (MAXPTN=400001)
1065 parameter (tend1 = 250d0)
1066 parameter (tend2 = 6.1d0)
1067 common /para1/ mul
1068cc SAVE /para1/
1069 common /para5/ iconfg, iordsc
1070 common /para7/ ioscar,nsmbbbar,nsmmeson
1071cc SAVE /para5/
1072 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1073 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1074 & xmass(MAXPTN), ityp(MAXPTN)
1075cc SAVE /prec2/
1076 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1077cc SAVE /prec4/
1078 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1079cc SAVE /prec5/
1080 common /ilist1/
1081 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1082 & ictype, icsta(MAXPTN),
1083 & nic(MAXPTN), icels(MAXPTN)
1084cc SAVE /ilist1/
1085 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1086cc SAVE /ilist4/
1087 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1088cc SAVE /ilist5/
1089 common /ilist6/ t, iopern, icolln
1090cc SAVE /ilist6/
1091 common /ana1/ ts(12)
1092cc SAVE /ana1/
1093 common/anim/nevent,isoft,isflag,izpc
1094cc SAVE /anim/
1095 COMMON /AREVT/ IAEVT, IARUN, MISS
1096 SAVE
1097
1098c save last collision info
1099 if (mod(ictype, 2) .eq. 0) then
1100 call savrec(iscat)
1101 call savrec(jscat)
1102 end if
1103
1104c1 get operation type
1105 call getict(t1)
1106c2 check freezeout condition
1107 if (iconfg .eq. 1 .and. t1 .gt. tlarge / 2d0) return 1
1108 if (iconfg .eq. 2 .or. iconfg .eq. 3) then
1109 if (t1 .gt. tend1) return 1
1110c if (ichkpt .eq. mul) then
1111c ii = 0
1112c do i = 1, mul
1113c gztemp = gz(i) + vz(i) * (t1 - ft(i))
1114c if (sqrt(t1 ** 2 - gztemp ** 2) .lt. tend) then
1115c ii = 1
1116c goto 1000
1117c end if
1118c end do
1119c 1000 continue
1120c if (ii .eq. 0) return 1
1121c end if
1122 end if
1123 if (iconfg .eq. 4 .or. iconfg .eq. 5) then
1124 if (t1 .gt. tend2) return 1
1125 end if
1126
1127clin-6/06/02 local freezeout for string melting,
1128c decide what partons have frozen out at time t1:
1129 if(isoft.eq.5) then
1130 call local(t1)
1131 endif
1132
1133c3 update iopern, t
1134
1135 iopern = iopern + 1
1136 t = t1
1137 if (mod(ictype, 2) .eq. 0) then
1138 icolln = icolln + 1
1139
1140c 4/18/01-ctest off
1141c write (2006, 1233) 'iscat=', iscat, 'jscat=', jscat,
1142c write (2006, *) 'iscat=', iscat, ' jscat=', jscat,
1143c 1 ityp(iscat), ityp(jscat)
1144c write (2006, 1233) 'iscat=', max(indx(iscat), indx(jscat)),
1145c & 'jscat=', min(indx(iscat), indx(jscat))
1146
1147c write (2006, 1234) ' icolln=', icolln, 't=', t
1148
1149c 1233 format (a10, i10, a10, i10)
1150c 1234 format (a15, i10, a5, f23.17, a5, f23.17)
1151 end if
1152
1153c4.1 deal with formation
1154 if (iconfg .eq. 1
1155 & .or. iconfg .eq. 2
1156 & .or. iconfg .eq. 4) then
1157 if (ictype .eq. 1 .or. ictype .eq. 2 .or.
1158 & ictype .eq. 5 .or. ictype .eq. 6) then
1159 call celasn
1160 end if
1161 end if
1162
1163c4.2 deal with collisions
1164
1165 if (ictype .ne. 1) then
1166
1167 iscat0 = iscat
1168 jscat0 = jscat
1169
1170c iscat is the larger one so that if it's a wall collision,
1171c it's still ok
1172 iscat = max0(iscat0, jscat0)
1173 jscat = min0(iscat0, jscat0)
1174
1175ctest off check icsta(i): 0 with f77 compiler
1176c write(9,*) 'BB:ictype,t1,iscat,jscat,icsta(i)=',
1177c 1 ictype,t1,iscat,jscat,icsta(iscat)
1178
1179c check collision time table error 'tterr'
1180clin-4/2008 to avoid out-of-bound error in next():
1181c if (jscat .ne. 0 .and. next(jscat) .ne. iscat)
1182c & then
1183c print *, 'iscat=', iscat, 'jscat=', jscat,
1184c & 'next(', jscat, ')=', next(jscat)
1185c
1186c if (ct(iscat) .lt. tlarge / 2d0) stop 'tterr'
1187c if (ct(jscat) .lt. tlarge / 2d0) stop 'tterr'
1188c end if
1189 if (jscat .ne. 0) then
1190 if(next(jscat) .ne. iscat) then
1191 print *, 'iscat=', iscat, 'jscat=', jscat,
1192 & 'next(', jscat, ')=', next(jscat)
1193 if (ct(iscat) .lt. tlarge / 2d0) stop 'tterr'
1194 if (ct(jscat) .lt. tlarge / 2d0) stop 'tterr'
1195 endif
1196 end if
1197clin-4/2008-end
1198
1199c4.2.1 collisions with wall
1200
1201c 8/19/02 avoid actual argument in common blocks of cellre:
1202 niscat=iscat
1203 njscat=jscat
1204c if (icsta(iscat) .ne. 0) call cellre(iscat, t)
1205c if (jscat .ne. 0) then
1206c if (icsta(jscat) .ne. 0) call cellre(jscat, t)
1207c end if
1208 if (icsta(iscat) .ne. 0) call cellre(niscat, t)
1209 if (jscat .ne. 0) then
1210 if (icsta(jscat) .ne. 0) call cellre(njscat, t)
1211 end if
1212
1213c4.2.2 collision between particles
1214
1215clin-6/2009 write out info for each collision:
1216c if (mod(ictype, 2) .eq. 0) call scat(t, iscat, jscat)
1217 if (mod(ictype, 2) .eq. 0) then
1218 if(ioscar.eq.3) then
1219cms write(95,*) 'event,miss,iscat,jscat=',iaevt,miss,iscat,jscat
1220 if(dmax1(abs(gx(iscat)),abs(gy(iscat)),
1221 1 abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)),
1222 2 abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat)))
1223 3 .lt.9999) then
1224cms write(95,200) ityp(iscat),px(iscat),py(iscat),
1225cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1226cms 2 gz(iscat),ft(iscat)
1227cms write(95,200) ityp(jscat),px(jscat),py(jscat),
1228cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1229cms 2 gz(jscat),ft(jscat)
1230 else
1231cms write(95,201) ityp(iscat),px(iscat),py(iscat),
1232cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1233cms 2 gz(iscat),ft(iscat)
1234cms write(95,201) ityp(jscat),px(jscat),py(jscat),
1235cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1236cms 2 gz(jscat),ft(jscat)
1237 endif
1238 endif
1239c
1240 call scat(t, iscat, jscat)
1241c
1242 if(ioscar.eq.3) then
1243 if(dmax1(abs(gx(iscat)),abs(gy(iscat)),
1244 1 abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)),
1245 2 abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat)))
1246 3 .lt.9999) then
1247cms write(95,200) ityp(iscat),px(iscat),py(iscat),
1248cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1249cms 2 gz(iscat),ft(iscat)
1250cms write(95,200) ityp(jscat),px(jscat),py(jscat),
1251cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1252cms 2 gz(jscat),ft(jscat)
1253 else
1254cms write(95,201) ityp(iscat),px(iscat),py(iscat),
1255cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1256cms 2 gz(iscat),ft(iscat)
1257cms write(95,201) ityp(jscat),px(jscat),py(jscat),
1258cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1259cms 2 gz(jscat),ft(jscat)
1260 endif
1261 endif
1262 endif
1263cyy format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
1264cyy format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
1265
1266 end if
1267
1268c5 update the interaction list
1269 call ulist(t)
1270
1271c6 update ifmpt. ichkpt
1272c old ichkpt and ifmpt are more conveniently used in ulist
1273 if (ifmpt .le. mul) then
1274 if (ictype .ne. 0 .and. ictype .ne. 3
1275 & .and. ictype .ne. 4) then
1276 ichkpt = ichkpt + 1
1277 ifmpt = ifmpt + 1
1278 end if
1279 end if
1280
1281 return
1282 end
1283
1284 subroutine savrec(i)
1285
1286 implicit double precision (a-h, o-z)
1287 parameter (MAXPTN=400001)
1288 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1289 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1290 & xmass(MAXPTN), ityp(MAXPTN)
1291cc SAVE /prec2/
1292 common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
1293 & pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
1294 & xmasss(MAXPTN), ityps(MAXPTN)
1295cc SAVE /prec3/
1296 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1297cc SAVE /prec5/
1298 common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
1299cc SAVE /prec6/
1300 SAVE
1301
1302 ityps(i) = ityp(i)
1303 gxs(i) = gx(i)
1304 gys(i) = gy(i)
1305 gzs(i) = gz(i)
1306 fts(i) = ft(i)
1307 pxs(i) = px(i)
1308 pys(i) = py(i)
1309 pzs(i) = pz(i)
1310 es(i) = e(i)
1311 xmasss(i) = xmass(i)
1312 etas(i) = eta(i)
1313 raps(i) = rap(i)
1314 taus(i) = tau(i)
1315
1316 return
1317 end
1318
1319 subroutine getict(t1)
1320 implicit double precision (a-h, o-z)
1321 parameter (MAXPTN=400001)
1322 common /para1/ mul
1323cc SAVE /para1/
1324 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1325 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1326 & xmass(MAXPTN), ityp(MAXPTN)
1327cc SAVE /prec2/
1328 common /ilist1/
1329 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1330 & ictype, icsta(MAXPTN),
1331 & nic(MAXPTN), icels(MAXPTN)
1332cc SAVE /ilist1/
1333 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1334cc SAVE /ilist4/
1335 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1336cc SAVE /ilist5/
1337 SAVE
1338
1339c neglect possibility of 2 collisions at the same time
1340c0 set initial conditions
1341
1342 t1 = tlarge
1343 iscat = 0
1344 jscat = 0
1345
1346c1 get next collision between particles
1347 do 1001 i = 1, ichkpt
1348 if (ot(i) .lt. t1) then
1349 t1 = ot(i)
1350 iscat = i
1351 end if
1352 1001 continue
1353 if (iscat .ne. 0) jscat = next(iscat)
1354
1355c2 get ictype
1356c 10/30/02 ictype=0:collision; 1:parton formation
1357 if (iscat .ne. 0 .and. jscat .ne. 0) then
1358 if (icsta(iscat) .eq. 0 .and. icsta(jscat) .eq. 0) then
1359 ictype = 0
1360 else
1361 ictype = 4
1362 end if
1363 else if (iscat .ne. 0 .or. jscat .ne. 0) then
1364 ictype = 3
1365 end if
1366c
1367 if (ifmpt .le. mul) then
1368 if (ft(ifmpt) .lt. t1) then
1369 ictype = 1
1370 t1 = ft(ifmpt)
1371 else if (ft(ifmpt) .eq. t1) then
1372 if (ictype .eq. 0) ictype = 2
1373 if (ictype .eq. 3) ictype = 5
1374 if (ictype .eq. 4) ictype = 6
1375 end if
1376 end if
1377
1378 return
1379 end
1380
1381 subroutine celasn
1382c this subroutine is used to assign a cell for a newly formed particle
1383c output: nic(MAXPTN) icels(MAXPTN) in the common /ilist1/
1384c icell, and icel(10,10,10) in the common /ilist2/
1385
1386 implicit double precision (a-h, o-z)
1387 parameter (MAXPTN=400001)
1388 common /para1/ mul
1389cc SAVE /para1/
1390 common /para5/ iconfg, iordsc
1391cc SAVE /para5/
1392 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1393 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1394 & xmass(MAXPTN), ityp(MAXPTN)
1395cc SAVE /prec2/
1396 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1397cc SAVE /prec4/
1398 common /ilist1/
1399 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1400 & ictype, icsta(MAXPTN),
1401 & nic(MAXPTN), icels(MAXPTN)
1402cc SAVE /ilist1/
1403 common /ilist2/ icell, icel(10,10,10)
1404cc SAVE /ilist2/
1405 common /ilist3/ size1, size2, size3, v1, v2, v3, size
1406cc SAVE /ilist3/
1407 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1408cc SAVE /ilist4/
1409 SAVE
1410
1411 external integ
1412
1413 i = ifmpt
1414 tt = ft(i)
1415 td = tt - size
1416 if (iconfg .eq. 1 .and. (size1 .eq. 0d0 .or.
1417 & size2 .eq. 0d0 .or. size3 .eq. 0d0)) then
1418 i1 = 11
1419 i2 = 11
1420 i3 = 11
1421 else if (iconfg .eq. 4 .or. td .le. 0d0) then
1422 i1 = integ(gx(i) / size1) + 6
1423 i2 = integ(gy(i) / size2) + 6
1424 i3 = integ(gz(i) / size3) + 6
1425 if (integ(gx(i) / size1) .eq. gx(i) / size1 .and.
1426 & vx(i) .lt. 0d0)
1427 & i1 = i1 - 1
1428 if (integ(gy(i) / size2) .eq. gy(i) / size2 .and.
1429 & vy(i) .lt. 0d0)
1430 & i2 = i2 - 1
1431 if (integ(gz(i) / size3) .eq. gz(i) / size3 .and.
1432 & vz(i) .lt. 0d0)
1433 & i3 = i3 - 1
1434 else
1435 i1 = integ(gx(i) / (size1 + v1 * td)) + 6
1436 i2 = integ(gy(i) / (size2 + v2 * td)) + 6
1437 i3 = integ(gz(i) / (size3 + v3 * td)) + 6
1438 if (integ(gx(i) / (size1 + v1 * td)) .eq. gx(i) /
1439 & (size1 + v1 * td) .and. vx(i) .lt. (i1 - 6) * v1)
1440 & i1 = i1 - 1
1441 if (integ(gy(i) / (size2 + v2 * td)) .eq. gy(i)/
1442 & (size2 + v2 * td) .and. vy(i) .lt. (i2 - 6) * v2)
1443 & i2 = i2 - 1
1444 if (integ(gz(i) / (size3 + v3 * td)) .eq. gz(i)/
1445 & (size3 + v3 * td) .and. vz(i) .lt. (i3 - 6) * v3)
1446 & i3 = i3 - 1
1447 end if
1448
1449 if (i1 .le. 0 .or. i1 .ge. 11 .or. i2 .le. 0 .or.
1450 & i2 .ge. 11 .or. i3 .le. 0 .or. i3 .ge. 11) then
1451 i1 = 11
1452 i2 = 11
1453 i3 = 11
1454 end if
1455
1456 if (i1 .eq. 11) then
1457 j = icell
1458 call newcre(i, j)
1459 icell = j
1460 icels(i) = 111111
1461 else
1462 j = icel(i1, i2, i3)
1463 call newcre(i, j)
1464 icel(i1, i2, i3) = j
1465 icels(i) = i1 * 10000 + i2 * 100 + i3
1466 end if
1467
1468 return
1469 end
1470
1471 integer function integ(x)
1472c this function is used to get the largest integer that is smaller than
1473c x
1474
1475 implicit double precision (a-h, o-z)
1476 SAVE
1477
1478 if (x .lt. 0d0) then
1479 integ = int(x - 1d0)
1480 else
1481 integ = int( x )
1482 end if
1483
1484 return
1485 end
1486
1487 subroutine cellre(i, t)
1488c this subroutine is used for changing the cell of a particle that
1489c collide with the wall
1490
1491 implicit double precision (a-h, o-z)
1492 parameter (MAXPTN=400001)
1493 common /para5/ iconfg, iordsc
1494cc SAVE /para5/
1495 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1496 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1497 & xmass(MAXPTN), ityp(MAXPTN)
1498cc SAVE /prec2/
1499 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1500cc SAVE /prec4/
1501 common /aurec1/ jxa, jya, jza
1502cc SAVE /aurec1/
1503 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
1504cc SAVE /aurec2/
1505 common /ilist1/
1506 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1507 & ictype, icsta(MAXPTN),
1508 & nic(MAXPTN), icels(MAXPTN)
1509cc SAVE /ilist1/
1510 common /ilist2/ icell, icel(10,10,10)
1511cc SAVE /ilist2/
1512 common /ilist3/ size1, size2, size3, v1, v2, v3, size
1513cc SAVE /ilist3/
1514 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1515cc SAVE /ilist4/
1516 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1517cc SAVE /ilist5/
1518 SAVE
1519
1520 logical good
1521
1522 external integ
1523
1524c this happens before update the /prec2/ common; in contrast with
1525c scat which happens after updating the glue common
1526
1527 t0 = t
1528
1529 1000 continue
1530
1531 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
1532 k = mod(icsta(i), 10)
1533
1534 if (k .eq. 1) then
1535 gx(i) = gx(i) - 10d0 * size1
1536 dgxa(i) = dgxa(i) + 10d0 * size1
1537 do 1001 ii = 1, ichkpt
1538 if (next(ii) .eq. i) then
1539 dgxa(ii) = dgxa(ii) - 10d0 * size1
1540 end if
1541 1001 continue
1542 end if
1543 if (k .eq. 2) then
1544 gx(i) = gx(i) + 10d0 * size1
1545 dgxa(i) = dgxa(i) - 10d0 * size1
1546 do 1002 ii = 1, ichkpt
1547 if (next(ii) .eq. i) then
1548 dgxa(ii) = dgxa(ii) + 10d0 * size1
1549 end if
1550 1002 continue
1551 end if
1552 if (k .eq. 3) then
1553 gy(i) = gy(i) - 10d0 * size2
1554 dgya(i) = dgya(i) + 10d0 * size2
1555 do 1003 ii = 1, ichkpt
1556 if (next(ii) .eq. i) then
1557 dgya(ii) = dgya(ii) - 10d0 * size2
1558 end if
1559 1003 continue
1560 end if
1561 if (k .eq. 4) then
1562 gy(i) = gy(i) + 10d0 * size2
1563 dgya(i) = dgya(i) - 10d0 * size2
1564 do 1004 ii = 1, ichkpt
1565 if (next(ii) .eq. i) then
1566 dgya(ii) = dgya(ii) + 10d0 * size2
1567 end if
1568 1004 continue
1569 end if
1570 if (iconfg .eq. 5) then
1571 if (k .eq. 5) then
1572 gz(i) = gz(i) - 10d0 * size3
1573 dgza(i) = dgza(i) + 10d0 * size3
1574 do 1005 ii = 1, ichkpt
1575 if (next(ii) .eq. i) then
1576 dgza(ii) = dgza(ii) - 10d0 * size3
1577 end if
1578 1005 continue
1579 end if
1580 if (k .eq. 6) then
1581 gz(i) = gz(i) + 10d0 * size3
1582 dgza(i) = dgza(i) - 10d0 * size3
1583 do 1006 ii = 1, ichkpt
1584 if (next(ii) .eq. i) then
1585 dgza(ii) = dgza(ii) + 10d0 * size3
1586 end if
1587 1006 continue
1588 end if
1589 end if
1590 else
1591 icels0 = icels(i)
1592
1593 i1 = icels0 / 10000
1594 i2 = (icels0 - i1 * 10000) / 100
1595 i3 = icels0 - i1 * 10000 - i2 * 100
1596
1597cc for particle inside the cube
1598 if (i1 .ge. 1 .and. i1 .le. 10
1599 & .and. i2 .ge. 1 .and. i2 .le. 10
1600 & .and. i3 .ge. 1 .and. i3 .le. 10) then
1601
1602c this assignment takes care of nic(i)=0 automatically
1603 if (icel(i1, i2, i3) .eq. i) icel(i1, i2, i3) = nic(i)
1604
1605c1 rearrange the old cell
1606
1607 call oldcre(i)
1608
1609c2 rearrange the new cell
1610
1611 k = mod(icsta(i), 10)
1612
1613c2.1 particle goes out of the cube
1614 if (iconfg .eq. 1) then
1615 good = (i1 .eq. 1 .and. k .eq. 2)
1616 & .or. (i1 .eq. 10 .and. k .eq. 1)
1617 & .or. (i2 .eq. 1 .and. k .eq. 4)
1618 & .or. (i2 .eq. 10 .and. k .eq. 3)
1619 & .or. (i3 .eq. 1 .and. k .eq. 6)
1620 & .or. (i3 .eq. 10 .and. k .eq. 5)
1621 end if
1622 if (iconfg .eq. 2) then
1623 good = (i3 .eq. 1 .and. k .eq. 6)
1624 & .or. (i3 .eq. 10 .and. k .eq. 5)
1625 end if
1626 if (good) then
1627
1628c j = icell
1629 call newcre(i, icell)
1630c icell = j
1631
1632 icels(i) = 111111
1633
1634c2.2 particle moves inside the cube
1635 else
1636
1637 if (k .eq. 1) i1 = i1 + 1
1638 if (k .eq. 2) i1 = i1 - 1
1639 if (k .eq. 3) i2 = i2 + 1
1640 if (k .eq. 4) i2 = i2 - 1
1641 if (k .eq. 5) i3 = i3 + 1
1642 if (k .eq. 6) i3 = i3 - 1
1643
1644 if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1645 if (i1 .eq. 0) then
1646 i1 = 10
1647 gx(i) = gx(i) + 10d0 * size1
1648 end if
1649 if (i1 .eq. 11) then
1650 i1 = 1
1651 gx(i) = gx(i) - 10d0 * size1
1652 end if
1653 if (i2 .eq. 0) then
1654 i2 = 10
1655 gy(i) = gy(i) + 10d0 * size2
1656 end if
1657 if (i2 .eq. 11) then
1658 i2 = 1
1659 gy(i) = gy(i) - 10d0 * size2
1660 end if
1661 if (iconfg .eq. 4) then
1662 if (i3 .eq. 0) then
1663 i3 = 10
1664 gz(i) = gz(i) + 10d0 * size3
1665 end if
1666 if (i3 .eq. 11) then
1667 i3 = 1
1668 gz(i) = gz(i) - 10d0 * size3
1669 end if
1670 end if
1671 end if
1672
1673 j = icel(i1, i2, i3)
1674
1675 call newcre(i, j)
1676c in case icel changes
1677
1678 icel(i1 ,i2, i3) = j
1679
1680 icels(i) = i1 * 10000 + i2 * 100 + i3
1681
1682 end if
1683
1684cc for particles outside the cube
1685 else
1686
1687 if (icell .eq. i) icell = nic(i)
1688
1689 call oldcre(i)
1690
1691 k = mod(icsta(i), 10)
1692
1693 ddt = t - ft(i)
1694 dtt = t - size
1695 if (dtt .le. 0d0) then
1696 i1 = integ((gx(i) + vx(i) * ddt) / size1) + 6
1697 i2 = integ((gy(i) + vy(i) * ddt) / size2) + 6
1698 i3 = integ((gz(i) + vz(i) * ddt) / size3) + 6
1699 else
1700 i1 = integ((gx(i) + vx(i) * ddt) /
1701 & (size1 + v1 * dtt)) + 6
1702 i2 = integ((gy(i) + vy(i) * ddt) /
1703 & (size2 + v2 * dtt)) + 6
1704 i3 = integ((gz(i) + vz(i) * ddt) /
1705 & (size3 + v3 * dtt)) + 6
1706 end if
1707
1708
1709 if (k .eq. 1) i1 = 1
1710 if (k .eq. 2) i1 = 10
1711 if (k .eq. 3) i2 = 1
1712 if (k .eq. 4) i2 = 10
1713 if (k .eq. 5) i3 = 1
1714 if (k .eq. 6) i3 = 10
1715
1716 j = icel(i1, i2, i3)
1717 call newcre(i, j)
1718 icel(i1, i2, i3) = j
1719
1720 icels(i) = i1 * 10000 + i2 * 100 + i3
1721
1722 end if
1723 end if
1724
1725 if (next(i) .ne. 0) then
1726 otmp = ot(next(i))
1727 ctmp = ct(next(i))
1728 end if
1729
1730 if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
1731 call dchout(i, k, t)
1732 else
1733 if (iconfg .eq. 1) then
1734 call dchin1(i, k, i1, i2, i3, t)
1735 else if (iconfg .eq. 2) then
1736 call dchin2(i, k, i1, i2, i3, t)
1737 else if (iconfg .eq. 4) then
1738 call dchin3(i, k, i1, i2, i3, t)
1739 end if
1740 end if
1741
1742 if (icsta(i) / 10 .eq. 11) then
1743 ot(next(i)) = otmp
1744 ct(next(i)) = ctmp
1745 next(next(i)) = i
1746 call wallc(i, i1, i2, i3, t0, tmin1)
1747 if (tmin1 .lt. ct(i)) then
1748 icsta(i) = icsta(i) + 10
1749 t0 = tmin1
1750 goto 1000
1751 end if
1752 end if
1753
1754 return
1755 end
1756
1757 subroutine oldcre(i)
1758c this subroutine is used to rearrange the old cell nic when a particle
1759c goes out of the cell
1760
1761 implicit double precision (a-h, o-z)
1762 parameter (MAXPTN=400001)
1763 common /ilist1/
1764 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1765 & ictype, icsta(MAXPTN),
1766 & nic(MAXPTN), icels(MAXPTN)
1767cc SAVE /ilist1/
1768 SAVE
1769
1770 if (nic(i) .eq. 0) return
1771
1772 j = nic(i)
1773
1774 if (nic(j) .eq. i) then
1775 nic(j) = 0
1776 return
1777 end if
1778
1779 do 10 while (nic(j) .ne. i)
1780 j = nic(j)
1781 10 continue
1782
1783 nic(j) = nic(i)
1784
1785 return
1786 end
1787
1788
1789 subroutine newcre(i, k)
1790c this subroutine is used to mk rearrange of the new cell a particle
1791c enters,
1792c input i
1793c output nic(i)
1794
1795 implicit double precision (a-h, o-z)
1796 parameter (MAXPTN=400001)
1797 common /ilist1/
1798 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1799 & ictype, icsta(MAXPTN),
1800 & nic(MAXPTN), icels(MAXPTN)
1801cc SAVE /ilist1/
1802 SAVE
1803
1804 if (k .eq. 0) then
1805 k = i
1806 nic(i) = 0
1807 else if (nic(k) .eq. 0) then
1808 nic(k) = i
1809 nic(i) = k
1810 else
1811 j = k
1812 do 10 while (nic(j) .ne. k)
1813 j = nic(j)
1814 10 continue
1815
1816 nic(j) = i
1817 nic(i) = k
1818
1819 end if
1820
1821 return
1822 end
1823
1824 subroutine scat(t, iscat, jscat)
1825
1826c this subroutine is used to calculate the 2 particle scattering
1827
1828 implicit double precision (a-h, o-z)
1829 SAVE
1830
1831 call newpos(t, iscat)
1832 call newpos(t, jscat)
1833 call newmom(t)
1834
1835 return
1836 end
1837
1838 subroutine newpos(t, i)
1839
1840c this subroutine is used to calculate the 2 particle scattering
1841c get new position
1842
1843 implicit double precision (a-h, o-z)
1844 parameter (MAXPTN=400001)
1845 common /para5/ iconfg, iordsc
1846cc SAVE /para5/
1847 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1848 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1849 & xmass(MAXPTN), ityp(MAXPTN)
1850cc SAVE /prec2/
1851 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1852cc SAVE /prec4/
1853 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1854cc SAVE /prec5/
1855 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1856cc SAVE /ilist5/
1857 SAVE
1858
1859 t=t
1860 dt1 = ct(i) - ft(i)
1861
1862 gx(i) = gx(i) + vx(i) * dt1
1863 gy(i) = gy(i) + vy(i) * dt1
1864 gz(i) = gz(i) + vz(i) * dt1
1865 ft(i) = ct(i)
1866
1867 if (iconfg .le. 3) then
1868 if (ft(i) .le. abs(gz(i))) then
1869 eta(i) = 1000000.d0
1870 else
1871 eta(i) = 0.5d0 * log((ft(i) + gz(i)) / (ft(i) - gz(i)))
1872 end if
1873 tau(i) = ft(i) / cosh(eta(i))
1874 end if
1875
1876 return
1877 end
1878
1879 subroutine newmom(t)
1880
1881c this subroutine is used to calculate the 2 particle scattering
1882
1883 implicit double precision (a-h, o-z)
1884
1885 parameter (hbarc = 0.197327054d0)
1886 parameter (MAXPTN=400001)
1887 parameter (pi = 3.14159265358979d0)
1888 COMMON /para1/ mul
1889cc SAVE /para1/
1890 common /para2/ xmp, xmu, alpha, rscut2, cutof2
1891cc SAVE /para2/
1892 common /para5/ iconfg, iordsc
1893cc SAVE /para5/
1894ctrans
1895 common /para6/ centy
1896cc SAVE /para6/
1897ctransend
1898 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1899 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1900 & xmass(MAXPTN), ityp(MAXPTN)
1901cc SAVE /prec2/
1902 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1903cc SAVE /prec4/
1904 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1905cc SAVE /prec5/
1906 common /aurec1/ jxa, jya, jza
1907cc SAVE /aurec1/
1908 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
1909cc SAVE /aurec2/
1910 common /ilist1/
1911 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1912 & ictype, icsta(MAXPTN),
1913 & nic(MAXPTN), icels(MAXPTN)
1914cc SAVE /ilist1/
1915 common /ilist3/ size1, size2, size3, v1, v2, v3, size
1916cc SAVE /ilist3/
3006c44b 1917 common /lora/ enenew, pxnew, pynew, pznew
1918cc SAVE /lora/
0119ef9a 1919 common /cprod/ xn1, xn2, xn3
1920cc SAVE /cprod/
1921 common /rndm2/ iff
1922cc SAVE /rndm2/
1923 common/anim/nevent,isoft,isflag,izpc
1924cc SAVE /anim/
1925 common /frzprc/
1926 & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1927 & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1928 & xmfrz(MAXPTN),
1929 & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1930cc SAVE /frzprc/
1931 SAVE
1932
1933 t=t
1934clin-6/06/02 no momentum change for partons already frozen out,
1935c however, spatial upgrade is needed to ensure overall system freezeout:
1936 if(isoft.eq.5) then
1937 if(ifrz(iscat).eq.1.or.ifrz(jscat).eq.1) then
1938 last(iscat) = jscat
1939 last(jscat) = iscat
1940 return
1941 endif
1942 endif
1943clin-6/06/02-end
1944
1945c iff is used to randomize the interaction to have both attractive and
1946c repulsive
1947
1948 iff = - iff
1949
1950 if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1951 icels1 = icels(iscat)
1952 i1 = icels1 / 10000
1953 j1 = (icels1 - i1 * 10000) / 100
1954 icels2 = icels(jscat)
1955 i2 = icels2 / 10000
1956 j2 = (icels2 - i2 * 10000) / 100
1957 if (iconfg .eq. 4) then
1958 k1 = icels1 - i1 * 10000 - j1 * 100
1959 k2 = icels2 - i2 * 10000 - j2 * 100
1960 end if
1961 end if
1962
1963 px1 = px(iscat)
1964 py1 = py(iscat)
1965 pz1 = pz(iscat)
1966 e1 = e(iscat)
1967 x1 = gx(iscat)
1968 y1 = gy(iscat)
1969 z1 = gz(iscat)
1970 t1 = ft(iscat)
1971 px2 = px(jscat)
1972 py2 = py(jscat)
1973 pz2 = pz(jscat)
1974 e2 = e(jscat)
1975
1976 if (iconfg .eq. 1) then
1977 x2 = gx(jscat)
1978 y2 = gy(jscat)
1979 z2 = gz(jscat)
1980 else if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1981 if (i1 - i2 .gt. 5) then
1982 x2 = gx(jscat) + 10d0 * size1
1983 else if (i1 - i2 .lt. -5) then
1984 x2 = gx(jscat) - 10d0 * size1
1985 else
1986 x2 = gx(jscat)
1987 end if
1988 if (j1 - j2 .gt. 5) then
1989 y2 = gy(jscat) + 10d0 * size2
1990 else if (j1 - j2 .lt. -5) then
1991 y2 = gy(jscat) - 10d0 * size2
1992 else
1993 y2 = gy(jscat)
1994 end if
1995 if (iconfg .eq. 4) then
1996 if (k1 - k2 .gt. 5) then
1997 z2 = gz(jscat) + 10d0 * size3
1998 else if (k1 - k2 .lt. -5) then
1999 z2 = gz(jscat) - 10d0 * size3
2000 else
2001 z2 = gz(jscat)
2002 end if
2003 else
2004 z2 = gz(jscat)
2005 end if
2006 else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2007 x2 = gx(jscat) + dgxa(jscat)
2008 y2 = gy(jscat) + dgya(jscat)
2009 if (iconfg .eq. 5) then
2010 z2 = gz(jscat) + dgza(jscat)
2011 else
2012 z2 = gz(jscat)
2013 end if
2014 end if
2015 t2 = ft(jscat)
2016ctrans
2017 rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2 -
2018 & (py1 + py2) ** 2 - (pz1 + pz2) ** 2
2019ctransend
2020 bex = (px1 + px2) / (e1 + e2)
2021 bey = (py1 + py2) / (e1 + e2)
2022 bez = (pz1 + pz2) / (e1 + e2)
2023
e2670fbe 2024 call lorenza(e1, px1, py1, pz1, bex, bey, bez)
0119ef9a 2025cc SAVE pxnew, ..., values for later use.
2026 px1 = pxnew
2027 py1 = pynew
2028 pz1 = pznew
2029 e1 = enenew
2030
2031 pp2 = pxnew ** 2 + pynew ** 2 + pznew ** 2
2032 call getht(iscat, jscat, pp2, that)
2033 theta = dacos(that / (2d0 * pp2) + 1d0)
2034 theta = dble(iff) * theta
2035
2036c we boost to the cm frame, get rotation axis, and rotate 1 particle
2037c momentum
2038
e2670fbe 2039 call lorenza(t1, x1, y1, z1, bex, bey, bez)
0119ef9a 2040
2041 x1 = pxnew
2042 y1 = pynew
2043 z1 = pznew
2044
e2670fbe 2045 call lorenza(t2, x2, y2, z2, bex, bey, bez)
0119ef9a 2046
2047 x2 = pxnew
2048 y2 = pynew
2049 z2 = pznew
2050
2051c notice now pxnew, ..., are new positions
2052 call cropro(x1-x2, y1-y2, z1-z2, px1, py1, pz1)
2053
2054 call xnormv(xn1, xn2, xn3)
2055
2056cbz1/29/99
2057c call rotate(xn1, xn2, xn3, theta, px1, py1, pz1)
2058 call zprota(xn1, xn2, xn3, theta, px1, py1, pz1)
2059cbz1/29/99end
2060
2061c we invert the momentum to get the other particle's momentum
2062 px2 = -px1
2063 py2 = -py1
2064 pz2 = -pz1
2065clin-4/13/01: modify in case m1, m2 are different:
2066c e2 = e1
2067 e2 = dsqrt(px2**2+py2**2+pz2**2+xmass(jscat)**2)
2068
2069c boost the 2 particle 4 momentum back to lab frame
e2670fbe 2070 call lorenza(e1, px1, py1, pz1, -bex, -bey, -bez)
0119ef9a 2071 px(iscat) = pxnew
2072 py(iscat) = pynew
2073 pz(iscat) = pznew
2074 e(iscat) = enenew
e2670fbe 2075 call lorenza(e2, px2, py2, pz2, -bex, -bey, -bez)
0119ef9a 2076 px(jscat) = pxnew
2077 py(jscat) = pynew
2078 pz(jscat) = pznew
2079 e(jscat) = enenew
2080
2081 vx(iscat) = px(iscat) / e(iscat)
2082 vy(iscat) = py(iscat) / e(iscat)
2083 vz(iscat) = pz(iscat) / e(iscat)
2084 vx(jscat) = px(jscat) / e(jscat)
2085 vy(jscat) = py(jscat) / e(jscat)
2086 vz(jscat) = pz(jscat) / e(jscat)
2087
2088 last(iscat) = jscat
2089 last(jscat) = iscat
2090
2091 if (iconfg .le. 3) then
2092 if (e(iscat) .le. abs(pz(iscat))) then
2093 rap(iscat) = 1000000.d0
2094 else
2095 rap(iscat) = 0.5d0 * log((e(iscat) + pz(iscat)) /
2096 & (e(iscat) - pz(iscat)))
2097 end if
2098
2099 if (e(jscat) .le. abs(pz(jscat))) then
2100 rap(jscat) = 1000000.d0
2101 else
2102 rap(jscat) = 0.5d0 * log((e(jscat) + pz(jscat)) /
2103 & (e(jscat) - pz(jscat)))
2104 end if
2105
2106ctrans
2107 rap1 = rap(iscat)
2108 rap2 = rap(jscat)
2109
2110 if ((rap1 .lt. centy + 0.5d0 .and.
2111 & rap1 .gt. centy - 0.5d0)) then
2112c write (9, *) sqrt(ft(iscat) ** 2 - gz(iscat) ** 2), rts2
2113 end if
2114 if ((rap2 .lt. centy + 0.5d0 .and.
2115 & rap2 .gt. centy - 0.5d0)) then
2116c write (9, *) sqrt(ft(jscat) ** 2 - gz(jscat) ** 2), rts2
2117 end if
2118ctransend
2119 end if
2120
2121 return
2122 end
2123
2124 subroutine getht(iscat, jscat, pp2, that)
2125
2126c this subroutine is used to get \hat{t} for a particular processes
2127
2128 implicit double precision (a-h, o-z)
2129
2130 parameter (hbarc = 0.197327054d0)
2131 parameter (MAXPTN=400001)
2132 common /para2/ xmp, xmu, alpha, rscut2, cutof2
2133cc SAVE /para2/
2134 common/anim/nevent,isoft,isflag,izpc
2135cc SAVE /anim/
2136cc external ran1
2137 common /rndm3/ iseedp
2138cc SAVE /rndm3/
2139 SAVE
2140
2141 iscat=iscat
2142 jscat=jscat
2143 iseed=iseedp
2144 xmu2 = (hbarc * xmu) ** 2
2145 xmp2 = xmp ** 2
2146 xm2 = xmu2 + xmp2
2147 rx=ran1(iseed)
2148 that = xm2*(1d0+1d0/((1d0-xm2/(4d0*pp2+xm2))*rx-1d0))
2149ctest off isotropic scattering:
2150c & + 1d0/((1d0 - xm2 / (4d0 * pp2 + xm2)) * ran1(2) - 1d0))
2151c if(izpc.eq.100) that=-4d0*pp2*ran1(2)
2152 if(izpc.eq.100) that=-4d0*pp2*rx
2153
2154 return
2155 end
2156
2157 subroutine ulist(t)
2158c this subroutine is used to update a new collision time list
2159c notice this t has been updated
2160
2161 implicit double precision (a-h, o-z)
2162 parameter (MAXPTN=400001)
2163 common /ilist1/
2164 & iscat, jscat, next(MAXPTN), last(MAXPTN),
2165 & ictype, icsta(MAXPTN),
2166 & nic(MAXPTN), icels(MAXPTN)
2167cc SAVE /ilist1/
2168 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
2169cc SAVE /ilist4/
2170 SAVE
2171
2172 if (ictype .eq. 1 .or. ictype .eq. 2 .or. ictype .eq. 5
2173 & .or. ictype .eq. 6) then
2174 l = ifmpt
2175 call ulist1(l, t)
2176 end if
2177 if (ictype .ne. 1) then
2178 l = iscat
2179 call ulist1(l, t)
2180 if (jscat .ne. 0) then
2181 l = jscat
2182 call ulist1(l, t)
2183 end if
2184 end if
2185
2186 return
2187 end
2188
2189 subroutine ulist1(l, t)
2190c this subroutine is used to update the interaction list when particle
2191c l is disturbed.
2192
2193 implicit double precision (a-h, o-z)
2194 parameter (MAXPTN=400001)
2195 common /para5/ iconfg, iordsc
2196cc SAVE /para5/
2197 common /ilist1/
2198 & iscat, jscat, next(MAXPTN), last(MAXPTN),
2199 & ictype, icsta(MAXPTN),
2200 & nic(MAXPTN), icels(MAXPTN)
2201cc SAVE /ilist1/
2202 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2203cc SAVE /ilist5/
2204 SAVE
2205
2206 icels0 = icels(l)
2207 i1 = icels0 / 10000
2208 i2 = (icels0 - i1 * 10000) / 100
2209 i3 = icels0 - i1 * 10000 - i2 * 100
2210c save collision info for use when the collision is a collision with wall
2211c otherwise wallc will change icsta
2212 k = mod(icsta(l), 10)
2213
2214 call wallc(l, i1, i2, i3, t, tmin1)
2215 tmin = tmin1
2216 nc = 0
2217
2218 if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
2219 call chkout(l, t, tmin, nc)
2220 else
2221 if (iconfg .eq. 1) then
2222 call chkin1(l, i1, i2, i3, t, tmin, nc)
2223 else if (iconfg .eq. 2) then
2224 call chkin2(l, i1, i2, i3, t, tmin, nc)
2225 else if (iconfg .eq. 4) then
2226 call chkin3(l, i1, i2, i3, t, tmin, nc)
2227 else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2228 call chkcel(l, i1, i2, i3, t, tmin, nc)
2229 end if
2230 end if
2231
2232 call fixtim(l, t, tmin1, tmin, nc)
2233
2234 return
2235 end
2236
2237 subroutine wallc(i, i1, i2, i3, t, tmin)
2238c this subroutine calculates the next time for collision with wall
2239c for particle i
2240c input particle label i,t
2241c output tmin collision time with wall, icsta(i) wall collision
2242c information
2243
2244 implicit double precision (a-h, o-z)
2245 parameter (MAXPTN=400001)
2246 common /para5/ iconfg, iordsc
2247cc SAVE /para5/
2248 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2249cc SAVE /ilist5/
2250 SAVE
2251
2252 tmin = tlarge
2253
2254 if (iconfg .le. 2 .or. iconfg .eq. 4) then
2255c if particle is inside the cube
2256 if ((i1 .ge. 1 .and. i1 .le. 10)
2257 & .or. (i2 .ge. 1 .and. i2 .le. 10)
2258 & .or. (i3 .ge. 1 .and. i3 .le. 10)) then
2259 call wallc1(i, i1, i2, i3, t, tmin)
2260c if particle is outside the cube
2261 else
2262 call wallcb(i, t, tmin)
2263 end if
2264 else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2265 call wallc2(i, i1, i2, i3, t, tmin)
2266 end if
2267
2268 return
2269 end
2270
2271 subroutine wallc1(i, i1, i2, i3, t, tmin)
2272c this subroutine is used to get wall collision time
2273c when particle is inside the cube, it sets the icsta at the same time
2274c input i,i1,i2,i3,t
2275c output tmin, icsta(i)
2276c note the icsta is not finally set. we need further judgement in
2277c fixtim
2278
2279 implicit double precision (a-h, o-z)
2280 parameter (MAXPTN=400001)
2281 common /para5/ iconfg, iordsc
2282cc SAVE /para5/
2283 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2284 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2285 & xmass(MAXPTN), ityp(MAXPTN)
2286cc SAVE /prec2/
2287 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2288cc SAVE /prec4/
2289 common /ilist1/
2290 & iscat, jscat, next(MAXPTN), last(MAXPTN),
2291 & ictype, icsta(MAXPTN),
2292 & nic(MAXPTN), icels(MAXPTN)
2293cc SAVE /ilist1/
2294 common /ilist3/ size1, size2, size3, v1, v2, v3, size
2295cc SAVE /ilist3/
2296 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2297cc SAVE /ilist5/
2298 SAVE
2299
2300 x1p = gx(i)
2301 x2p = gy(i)
2302 x3p = gz(i)
2303 tf = ft(i)
2304 v1p = vx(i)
2305 v2p = vy(i)
2306 v3p = vz(i)
2307
2308 if (t .lt. size .and. tf .lt. size) then
2309
2310 if (v1p .gt. 0d0) then
2311 t1 = ((dble(i1) - 5d0) * size1 - x1p) / v1p + tf
2312 else if (v1p .lt. 0d0) then
2313 t1 = ((dble(i1) - 6d0) * size1 - x1p) / v1p + tf
2314 else
2315 t1 = tlarge
2316 end if
2317
2318 if (v2p .gt. 0d0) then
2319 t2 = ((dble(i2) - 5d0) * size2 - x2p) / v2p + tf
2320 else if (v2p .lt. 0d0) then
2321 t2 = ((dble(i2) - 6d0) * size2 - x2p) / v2p + tf
2322 else
2323 t2 = tlarge
2324 end if
2325
2326 if (v3p .gt. 0d0) then
2327 t3 = ((dble(i3) - 5d0) * size3 - x3p) / v3p + tf
2328 else if (v3p .lt. 0d0) then
2329 t3 = ((dble(i3) - 6d0) * size3 - x3p) / v3p + tf
2330 else
2331 t3 = tlarge
2332 end if
2333
2334c if a particle is on the wall, we don't collide it on the same wall
2335
2336c if (t1 .eq. 0d0) t1 = tlarge
2337c if (t2 .eq. 0d0) t2 = tlarge
2338c if (t3 .eq. 0d0) t3 = tlarge
2339
2340 tmin = min(t1, t2, t3)
2341
2342c set icsta,
2343c after checking this is not an earlier collision comparing with
2344c a collision with another particle, we need to set icsta=0
2345c after checking whether there is also a particle collision
2346c at the same time, we need to reset the second bit of icsta
2347
2348 if (tmin .eq. t1) then
2349 if (v1p .gt. 0d0) then
2350 icsta(i) = 101
2351 else
2352 icsta(i) = 102
2353 end if
2354 end if
2355
2356 if (tmin .eq. t2) then
2357 if (v2p .gt. 0d0) then
2358 icsta(i) = 103
2359 else
2360 icsta(i) = 104
2361 end if
2362 end if
2363
2364 if (tmin .eq. t3) then
2365 if (v3p .gt. 0d0) then
2366 icsta(i) = 105
2367 else
2368 icsta(i) = 106
2369 end if
2370 end if
2371
2372 if (tmin .le. size) return
2373
2374 end if
2375
2376 if (v1p .gt. (i1 - 5) * v1) then
2377 t1 = ((i1 - 5) * (size1 - v1 * size) +
2378 & v1p * tf - x1p) / (v1p - (i1 - 5) * v1)
2379 else if (v1p .lt. (i1 - 6) * v1) then
2380 t1 = ((i1 - 6) * (size1 - v1 * size) +
2381 & v1p * tf - x1p) / (v1p - (i1 - 6) * v1)
2382 else
2383 t1 = tlarge
2384 end if
2385
2386 if (v2p .gt. (i2 - 5) * v2) then
2387 t2 = ((i2 - 5) * (size2 - v2 * size) +
2388 & v2p * tf - x2p) / (v2p - (i2 - 5) * v2)
2389 else if (v2p .lt. (i2 - 6) * v2) then
2390 t2 = ((i2 - 6) * (size2 - v2 * size) +
2391 & v2p * tf - x2p) / (v2p - (i2 - 6) * v2)
2392 else
2393 t2 = tlarge
2394 end if
2395
2396 if (v3p .gt. (i3 - 5) * v3) then
2397 t3 = ((i3 - 5) * (size3 - v3 * size) +
2398 & v3p * tf - x3p) / (v3p - (i3 - 5) * v3)
2399 else if (v3p .lt. (i3 - 6) * v3) then
2400 t3 = ((i3 - 6) * (size3 - v3 * size) +
2401 & v3p * tf - x3p) / (v3p - (i3 - 6) * v3)
2402 else
2403 t3 = tlarge
2404 end if
2405
2406c if a particle is on the wall, we don't collide it on the same wall
2407
2408c if (t1 .eq. 0d0) t1 = tlarge
2409c if (t2 .eq. 0d0) t2 = tlarge
2410c if (t3 .eq. 0d0) t3 = tlarge
2411
2412 tmin = min(t1, t2, t3)
2413
2414c set icsta,
2415c after checking this is not an earlier collision comparing with
2416c a collision with another particle, we need to set icsta=0
2417c after checking whether there is also a particle collision
2418c at the same time, we need to reset the second bit of icsta
2419
2420 if (tmin .eq. t1) then
2421 if (v1p .gt. (i1 - 5) * v1) then
2422 icsta(i) = 101
2423 else
2424 icsta(i) = 102
2425 end if
2426 end if
2427
2428 if (tmin .eq. t2) then
2429 if (v2p .gt. (i2 - 5) * v2) then
2430 icsta(i) = 103
2431 else
2432 icsta(i) = 104
2433 end if
2434 end if
2435
2436 if (tmin .eq. t3) then
2437 if (v3p .gt. (i3 - 5) * v3) then
2438 icsta(i) = 105
2439 else
2440 icsta(i) = 106
2441 end if
2442 end if
2443
2444 return
2445 end
2446
2447 subroutine wallc2(i, i1, i2, i3, t, tmin)
2448c this subroutine is used to get wall collision time
2449c when particle is inside the cube, it sets the icsta at the same time
2450c input i,i1,i2,i3,t
2451c output tmin, icsta(i)
2452c note the icsta is not finally set. we need further judgement in
2453c fixtim
2454
2455 implicit double precision (a-h, o-z)
2456 parameter (MAXPTN=400001)
2457
2458 common /para5/ iconfg, iordsc
2459cc SAVE /para5/
2460 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2461 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2462 & xmass(MAXPTN), ityp(MAXPTN)
2463cc SAVE /prec2/
2464 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2465cc SAVE /prec4/
2466 common /ilist1/
2467 & iscat, jscat, next(MAXPTN), last(MAXPTN),
2468 & ictype, icsta(MAXPTN),
2469 & nic(MAXPTN), icels(MAXPTN)
2470cc SAVE /ilist1/
2471 common /ilist3/ size1, size2, size3, v1, v2, v3, size
2472cc SAVE /ilist3/
2473 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2474cc SAVE /ilist5/
2475 SAVE
2476
2477 i1=i1
2478 i2=i2
2479 i3=i3
2480 t=t
2481 x1p = gx(i)
2482 x2p = gy(i)
2483 x3p = gz(i)
2484 tf = ft(i)
2485 v1p = vx(i)
2486 v2p = vy(i)
2487 v3p = vz(i)
2488
2489 if (v1p .gt. 0d0) then
2490 t1 = (5d0 * size1 - x1p) / v1p + tf
2491 else if (v1p .lt. 0d0) then
2492 t1 = (-5d0 * size1 - x1p) / v1p + tf
2493 else
2494 t1 = tlarge
2495 end if
2496
2497 if (v2p .gt. 0d0) then
2498 t2 = (5d0 * size2 - x2p) / v2p + tf
2499 else if (v2p .lt. 0d0) then
2500 t2 = (- 5d0 * size2 - x2p) / v2p +tf
2501 else
2502 t2 = tlarge
2503 end if
2504
2505 if (iconfg .eq. 5) then
2506 if (v3p .gt. 0d0) then
2507 t3 = (5d0 * size3 - x3p) / v3p + tf
2508 else if (v3p .lt. 0d0) then
2509 t3 = (- 5d0 * size3 - x3p) / v3p +tf
2510 else
2511 t3 = tlarge
2512 end if
2513 else
2514 t3 = tlarge
2515 end if
2516
2517 tmin = min(t1, t2, t3)
2518
2519c set icsta,
2520c after checking this is not an earlier collision comparing with
2521c a collision with another particle, we need to set icsta=0
2522c after checking whether there is also a particle collision
2523c at the same time, we need to reset the second bit of icsta
2524
2525 if (tmin .eq. t1) then
2526 if (v1p .gt. 0d0) then
2527 icsta(i) = 101
2528 else
2529 icsta(i) = 102
2530 end if
2531 end if
2532
2533 if (tmin .eq. t2) then
2534 if (v2p .gt. 0d0) then
2535 icsta(i) = 103
2536 else
2537 icsta(i) = 104
2538 end if
2539 end if
2540
2541 if (tmin .eq. t3) then
2542 if (v3p .gt. 0d0) then
2543 icsta(i) = 105
2544 else
2545 icsta(i) = 106
2546 end if
2547 end if
2548
2549 return
2550 end
2551
2552 subroutine wallcb(i, t, tmin)
2553c this subroutine is used to calculate the wall collision time
2554c when the particle is outside the cube
2555c input i,t
2556c output tmin,icsta(i)
2557c note the icsta is not finally set. we need further judgement in
2558c fixtim
2559
2560 implicit double precision (a-h, o-z)
2561 parameter (MAXPTN=400001)
2562
2563 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2564 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2565 & xmass(MAXPTN), ityp(MAXPTN)
2566cc SAVE /prec2/
2567 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2568cc SAVE /prec4/
2569 common /ilist1/
2570 & iscat, jscat, next(MAXPTN), last(MAXPTN),
2571 & ictype, icsta(MAXPTN),
2572 & nic(MAXPTN), icels(MAXPTN)
2573cc SAVE /ilist1/
2574 common /ilist3/ size1, size2, size3, v1, v2, v3, size
2575cc SAVE /ilist3/
2576 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2577cc SAVE /ilist5/
2578 SAVE
2579
2580c check if there is a collision by looking at the closest approach point
2581c and see if it's inside the cube
2582
2583 if (size1 .eq. 0d0 .or. size2 .eq. 0d0 .or.
2584 & size3 .eq. 0d0) return
2585
2586 x1p = gx(i)
2587 x2p = gy(i)
2588 x3p = gz(i)
2589 v1p = vx(i)
2590 v2p = vy(i)
2591 v3p = vz(i)
2592 tf = ft(i)
2593
2594 if (t .lt. size .and. tf .lt. size) then
2595 if (x1p .lt. - 5d0 * size1 .and. v1p .gt. 0d0) then
2596 t1 = (- 5d0 * size1 - x1p) / v1p + tf
2597 else if(x1p .gt. 5d0 * size1 .and. v1p .lt. 0d0) then
2598 t1 = - (x1p - 5d0 * size1) / v1p + tf
2599 else
2600 t1 = tlarge
2601 end if
2602
2603 if (t1 .ne. tlarge) then
2604 x2pp = x2p + v2p * (t1 - tf)
2605 x3pp = x3p + v3p * (t1 - tf)
2606 if (x2pp .le. - 5d0 * size2 .or. x2pp .ge. 5d0 * size2
2607 & .or. x3pp .le. - 5d0 * size3
2608 & .or. x3pp .ge. 5d0 * size3)
2609 & t1 = tlarge
2610 end if
2611
2612 if (x2p .lt. - 5d0 * size2 .and. v2p .gt. 0d0) then
2613 t2 = (- 5d0 * size2 - x2p) / v2p + tf
2614 else if(x2p .gt. 5d0 * size2 .and. v2p .lt. 0d0) then
2615 t2 = - (x2p - 5d0 * size2) / v2p + tf
2616 else
2617 t2 = tlarge
2618 end if
2619
2620 if (t2 .ne. tlarge) then
2621 x1pp = x1p + v1p * (t2 - tf)
2622 x3pp = x3p + v3p * (t2 - tf)
2623 if (x1pp .le. - 5d0 * size1 .or. x1pp .ge. 5d0 * size1
2624 & .or. x3pp .le. - 5d0 * size3 .or. x3pp .ge. 5d0 * size3)
2625 & t2 = tlarge
2626 end if
2627
2628 if (x3p .lt. - 5d0 * size3 .and. v3p .gt. 0d0) then
2629 t3 = (- 5d0 * size3 - x3p) / v3p + tf
2630 else if(x3p .gt. 5d0 * size3 .and. v3p .lt. 0d0) then
2631 t3 = - (x3p - 5d0 * size3) / v3p + tf
2632 else
2633 t3 = tlarge
2634 end if
2635
2636 if (t3 .ne. tlarge) then
2637 x1pp = x1p + v1p * (t3 - tf)
2638 x2pp = x2p + v2p * (t3 - tf)
2639 if (x1pp .le. - 5d0 * size1 .or. x1pp .ge. 5d0 * size1
2640 & .or. x2pp .le. - 5d0 * size2 .or. x2pp .ge. 5d0 * size2)
2641 & t3 = tlarge
2642 end if
2643
2644 tmin = min(t1, t2, t3)
2645
2646c set icsta,
2647c after checking this is not an earlier collision comparing with
2648c a collision with another particle, we need to set icsta=0
2649c after checking whether there is also a particle collision
2650c at the same time, we need to reset the second bit of icsta
2651
2652 if (tmin .eq. t1) then
2653 if (v1p .gt. 0d0) then
2654 icsta(i) = 101
2655 else
2656 icsta(i) = 102
2657 end if
2658 end if
2659
2660 if (tmin .eq. t2) then
2661 if (v2p .gt. 0d0) then
2662 icsta(i) = 103
2663 else
2664 icsta(i) = 104
2665 end if
2666 end if
2667
2668 if (tmin .eq. t3) then
2669 if (v3p .gt. 0d0) then
2670 icsta(i) = 105
2671 else
2672 icsta(i) = 106
2673 end if
2674 end if
2675
2676 if (tmin .le. size) return
2677
2678 end if
2679
2680c notice now x1q, x2q, x3q are coordinates at time t
2681 x1q = x1p + v1p * (t - tf)
2682 x2q = x2p + v2p * (t - tf)
2683 x3q = x3p + v3p * (t - tf)
2684
2685 if (x1q .lt. - 5d0 * (size1 + v1 * (t - size)) .and.
2686 & v1p .gt. - 5d0 * v1) then
2687 t1 = (- 5d0 * (size1 - v1 * size) + v1p * tf - x1p) /
2688 & (v1p - (- 5d0) * v1)
2689 icsta1 = 101
2690 else if (x1q .gt. 5d0 * (size1 + v1 * (t-size)) .and.
2691 & v1p .lt. 5d0 * v1) then
2692 t1 = (5d0 * (size1 - v1 * size) + v1p * tf - x1p) /
2693 & (v1p - 5d0 * v1)
2694 icsta1 = 102
2695 else
2696 t1 = tlarge
2697 end if
2698
2699 if (t1 .ne. tlarge) then
2700 x2pp = x2p + v2p * (t1 - tf)
2701 x3pp = x3p + v3p * (t1 - tf)
2702 if (x2pp .le. - 5d0 * (size2 + v2 * (t1 - size))
2703 & .or. x2pp .ge. 5d0 * (size2 + v2 * (t1 - size))
2704 & .or. x3pp .le. - 5d0 * (size3 + v3 * (t1 - size))
2705 & .or. x3pp .ge. 5d0 * (size3 + v3 * (t1 - size)))
2706 & t1 = tlarge
2707 end if
2708
2709 if (x2q .lt. - 5d0 * (size2 + v2 * (t - size)) .and.
2710 & v2p .gt. - 5d0 * v2) then
2711 t2 = (- 5d0 * (size2 - v2 * size) + v2p * tf - x2p) /
2712 & (v2p - (- 5d0) * v2)
2713 icsta2 = 103
2714 else if (x2q .gt. 5d0 * (size2 + v2 * (t - size)) .and.
2715 & v2p .lt. 5d0 * v2) then
2716 t2 = (5d0 * (size2 - v2 * size) + v2p * tf - x2p) /
2717 & (v2p - 5d0 * v2)
2718 icsta2 = 104
2719 else
2720 t2 = tlarge
2721 end if
2722
2723 if (t2 .ne. tlarge) then
2724 x1pp = x1p + v1p * (t2 - tf)
2725 x3pp = x3p + v3p * (t2 - tf)
2726 if (x1pp .le. - 5d0 * (size1 + v1 * (t2 - size))
2727 & .or. x1pp .ge. 5d0 * (size1 + v1 * (t2 - size))
2728 & .or. x3pp .le. - 5d0 * (size3 + v3 * (t2 - size))
2729 & .or. x3pp .ge. 5d0 * (size3 + v3 * (t2 - size)))
2730 & t2 = tlarge
2731 end if
2732
2733 if (x3q .lt. - 5d0 * (size3 + v3 * (t - size)) .and.
2734 & v3p .gt. - 5d0 * v3) then
2735 t3 = (- 5d0 * (size3 - v3 * size) + v3p * tf - x3p) /
2736 & (v3p - (- 5d0) * v3)
2737 icsta3 = 105
2738 else if (x3q .gt. 5d0 * (size3 + v3 * (t - size)) .and.
2739 & v3p .lt. 5d0 * v3) then
2740 t3 = (5d0 * (size3 - v3 * size) + v3p * tf - x3p) /
2741 & (v3p - 5d0 * v3)
2742 icsta3 = 106
2743 else
2744 t3 = tlarge
2745 end if
2746
2747 if (t3 .ne. tlarge) then
2748 x2pp = x2p + v2p * (t3 - tf)
2749 x1pp = x1p + v1p * (t3 - tf)
2750 if (x2pp .le. - 5d0 * (size2 + v2 * (t3 - size))
2751 & .or. x2pp .ge. 5d0 * (size2 + v2 * (t3 - size))
2752 & .or. x1pp .le. - 5d0 * (size1 + v1 * (t3 - size))
2753 & .or. x1pp .ge. 5d0 * (size1 + v1 * (t3 - size)))
2754 & t3 = tlarge
2755 end if
2756
2757 tmin = min(t1, t2, t3)
2758
2759c set icsta,
2760c after checking this is not an earlier collision comparing with
2761c a collision with another particle, we need to set icsta=0
2762c after checking whether there is also a particle collision
2763c at the same time, we need to reset the second bit of icsta
2764
2765 if (tmin .eq. t1) then
2766 icsta(i) = icsta1
2767 else if (tmin .eq. t2) then
2768 icsta(i) = icsta2
2769 else if (tmin .eq. t3) then
2770 icsta(i) = icsta3
2771 end if
2772
2773 return
2774 end
2775
2776 subroutine chkout(l, t, tmin, nc)
2777c this subroutine is used to check the collisions with particles in
2778c surface cells to see if we can get a smaller collision time than tmin
2779c with particle nc, when the colliding particle is outside the cube
2780c input l,t,tmin,nc
2781c output tmin, nc
2782
2783 implicit double precision (a-h, o-z)
2784 parameter (MAXPTN=400001)
2785
2786 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2787 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2788 & xmass(MAXPTN), ityp(MAXPTN)
2789cc SAVE /prec2/
2790 SAVE
2791
2792 m1 = 11
2793 m2 = 11
2794 m3 = 11
2795 call chkcel(l, m1, m2, m3, t, tmin, nc)
2796
2797 do 1003 i = 1, 10
2798 do 1002 j = 1, 10
2799 do 1001 k = 1, 10
2800 if (i .eq. 1 .or. i .eq. 10 .or. j .eq. 1
2801 & .or. j .eq. 10 .or. k .eq. 1 .or. k .eq. 10)
2802 & call chkcel(l, i, j, k, t, tmin, nc)
2803 1001 continue
2804 1002 continue
2805 1003 continue
2806
2807 return
2808 end
2809
2810 subroutine chkin1(l, i1, i2, i3, t, tmin, nc)
2811c this subroutine is used to check collisions for particle inside
2812c the cube
2813c and update the afftected particles through chkcel
2814
2815 implicit double precision (a-h, o-z)
2816 SAVE
2817
2818c itest is a flag to make sure the 111111 cell is checked only once
2819 itest = 0
2820
2821 do 1003 i = i1 - 1, i1 + 1
2822 do 1002 j = i2 - 1, i2 + 1
2823 do 1001 k = i3 - 1, i3 + 1
2824 if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
2825 & j .le. 10 .and. k .ge. 1 .and. k .le. 10) then
2826 call chkcel(l, i, j, k, t, tmin, nc)
2827 else if (itest .eq. 0) then
2828 m1 = 11
2829 m2 = 11
2830 m3 = 11
2831 call chkcel(l, m1, m2, m3, t, tmin, nc)
2832 itest = 1
2833 end if
2834 1001 continue
2835 1002 continue
2836 1003 continue
2837
2838 return
2839 end
2840
2841 subroutine chkin2(l, i1, i2, i3, t, tmin, nc)
2842c this subroutine is used to check collisions for particle inside
2843c the cube
2844c and update the afftected particles through chkcel
2845
2846 implicit double precision (a-h, o-z)
2847 SAVE
2848
2849c itest is a flag to make sure the 111111 cell is checked only once
2850 itest = 0
2851
2852 do 1003 i = i1 - 1, i1 + 1
2853 do 1002 j = i2 - 1, i2 + 1
2854 do 1001 k = i3 - 1, i3 + 1
2855 ia = i
2856 ib = j
2857 ic = k
2858 if (k .ge. 1 .and. k .le. 10) then
2859 if (i .eq. 0) ia = 10
2860 if (i .eq. 11) ia = 1
2861 if (j .eq. 0) ib = 10
2862 if (j .eq. 11) ib = 1
2863 call chkcel(l, ia, ib, ic, t, tmin, nc)
2864 end if
2865 1001 continue
2866 1002 continue
2867 1003 continue
2868
2869 return
2870 end
2871
2872 subroutine chkin3(l, i1, i2, i3, t, tmin, nc)
2873c this subroutine is used to check collisions for particle inside
2874c the cube
2875c and update the afftected particles through chkcel
2876
2877 implicit double precision (a-h, o-z)
2878 SAVE
2879
2880c itest is a flag to make sure the 111111 cell is checked only once
2881 itest = 0
2882
2883 do 1003 i = i1 - 1, i1 + 1
2884 do 1002 j = i2 - 1, i2 + 1
2885 do 1001 k = i3 - 1, i3 + 1
2886 if (i .eq. 0) then
2887 ia = 10
2888 else if (i .eq. 11) then
2889 ia = 1
2890 else
2891 ia = i
2892 end if
2893 if (j .eq. 0) then
2894 ib = 10
2895 else if (j .eq. 11) then
2896 ib = 1
2897 else
2898 ib = j
2899 end if
2900 if (k .eq. 0) then
2901 ic = 10
2902 else if (k .eq. 11) then
2903 ic = 1
2904 else
2905 ic = k
2906 end if
2907 call chkcel(l, ia, ib, ic, t, tmin, nc)
2908 1001 continue
2909 1002 continue
2910 1003 continue
2911
2912 return
2913 end
2914
2915 subroutine chkcel(il, i1, i2, i3, t, tmin, nc)
2916c this program is used to check through all the particles
2917c in the cell (i1,i2,i3) and see if we can get a particle collision
2918c with time less than the original input tmin ( the collision time of
2919c il with the wall
2920c and update the affected particles
2921
2922 implicit double precision (a-h, o-z)
2923 parameter (MAXPTN=400001)
2924
2925 common /para5/ iconfg, iordsc
2926cc SAVE /para5/
2927 common /ilist1/
2928 & iscat, jscat, next(MAXPTN), last(MAXPTN),
2929 & ictype, icsta(MAXPTN),
2930 & nic(MAXPTN), icels(MAXPTN)
2931cc SAVE /ilist1/
2932 common /ilist2/ icell, icel(10, 10, 10)
2933cc SAVE /ilist2/
2934 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
2935cc SAVE /ilist4/
2936 SAVE
2937
2938 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2939 jj = ichkpt
2940 do 1001 j = 1, jj
2941 call ck(j, ick)
2942c 10/24/02 get rid of argument usage mismatch in ud2():
2943 jud2=j
2944c if (ick .eq. 1) call ud2(j, il, t, tmin, nc)
2945 if (ick .eq. 1) call ud2(jud2, il, t, tmin, nc)
2946 1001 continue
2947 return
2948 end if
2949
2950 if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
2951 l = icell
2952 else
2953 l = icel(i1, i2, i3)
2954 end if
2955
2956c if there is no particle
2957 if (l .eq. 0) then
2958 return
2959 end if
2960 j = nic(l)
2961c if there is only one particle
2962 if (j .eq. 0) then
2963 call ck(l, ick)
2964 if (ick .eq. 1) call ud2(l, il, t, tmin, nc)
2965
2966c if there are many particles
2967 else
2968
2969c we don't worry about the other colliding particle because it's
2970c set in last(), and will be checked in ud2
2971
2972 call ck(l, ick)
2973 if (ick .eq. 1) call ud2(l, il, t, tmin, nc)
2974
2975 do 10 while(j .ne. l)
2976 call ck(j, ick)
2977 if (ick .eq. 1) call ud2(j, il, t, tmin, nc)
2978 j = nic(j)
2979 10 continue
2980 end if
2981
2982 return
2983 end
2984
2985 subroutine ck(l, ick)
2986c this subroutine is used for chcell to check whether l should be
2987c checked or not for updating tmin, nc
2988c input l
2989c output ick
2990c if ick=1, l should be checked, otherwise it should not be.
2991
2992 implicit double precision (a-h, o-z)
2993 parameter (MAXPTN=400001)
2994
2995 common /ilist1/
2996 & iscat, jscat, next(MAXPTN), last(MAXPTN),
2997 & ictype, icsta(MAXPTN),
2998 & nic(MAXPTN), icels(MAXPTN)
2999cc SAVE /ilist1/
3000 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
3001cc SAVE /ilist4/
3002 SAVE
3003
3004 ick = 1
3005 if (ictype .eq. 1) then
3006 if (l .eq. ifmpt) ick = 0
3007 else if (ictype .eq. 0 .or. ictype .eq. 3 .or.
3008 & ictype .eq. 4) then
3009 if (l .eq. iscat .or. l .eq. jscat) ick = 0
3010 else
3011 if (l .eq. iscat .or. l .eq. jscat .or.
3012 & l .eq. ifmpt) ick = 0
3013 end if
3014c notice il is either iscat or jscat, or ifmpt, we deal with them
3015c seperately according to ictype
3016
3017 return
3018 end
3019
3020 subroutine dchout(l, ii, t)
3021c this subroutine is used to check collisions of l with particles when
3022c l is outside the cube and the collision just happened is a collision
3023c including a collision with wall (hence we need to use dcheck to throw
3024c away old collisions that are not in the new neighboring cells.
3025
3026c input l,t
3027
3028 implicit double precision (a-h, o-z)
3029 parameter (MAXPTN=400001)
3030
3031 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3032 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3033 & xmass(MAXPTN), ityp(MAXPTN)
3034cc SAVE /prec2/
3035 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3036cc SAVE /prec4/
3037 common /ilist3/ size1, size2, size3, v1, v2, v3, size
3038cc SAVE /ilist3/
3039 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3040cc SAVE /ilist5/
3041 SAVE
3042
3043 external integ
3044
3045 tt = ft(l)
3046 td = t - size
3047 x1 = gx(l) + vx(l) * (t - tt)
3048 x2 = gy(l) + vy(l) * (t - tt)
3049 x3 = gz(l) + vz(l) * (t - tt)
3050 if (td .le. 0d0) then
3051 i1 = integ(x1 / size1) + 6
3052 i2 = integ(x2 / size2 ) + 6
3053 i3 = integ(x3 / size3 ) + 6
3054 if (integ(x1 / size1) .eq. x1 / size1 .and. vx(l) .lt. 0d0)
3055 & i1 = i1 - 1
3056 if (integ(x2 / size2) .eq. x2 / size2 .and. vy(l) .lt. 0d0)
3057 & i2 = i2 - 1
3058 if (integ(x3 / size3) .eq. x3 / size3 .and. vz(l) .lt. 0d0)
3059 & i3 = i3 - 1
3060 else
3061 i1 = integ(x1 / (size1 + v1 * td)) + 6
3062 i2 = integ(x2 / (size2 + v2 * td)) + 6
3063 i3 = integ(x3 / (size3 + v3 * td)) + 6
3064c 10/24/02 (i) below should be (l):
3065 if (integ(x1 / (size1 + v1 * td)) .eq.
3066 & x1 / (size1 +v1 * td) .and.
3067 & vx(l) .lt. (i1 - 6) * v1) i1 = i1 - 1
3068c & vx(i) .lt. (i1 - 6) * v1) i1 = i1 - 1
3069 if (integ(x2 / (size2 + v2 * td)) .eq.
3070 & x2 / (size2 + v2 * td) .and.
3071 & vy(l) .lt. (i2 - 6) * v2) i2 = i2 - 1
3072c & vy(i) .lt. (i2 - 6) * v2) i2 = i2 - 1
3073 if (integ(x3 / (size3 + v3 * td)) .eq.
3074 & x3 / (size3 + v3 * td) .and.
3075 & vz(l) .lt. (i3 - 6) * v3) i3 = i3 - 1
3076c & vz(i) .lt. (i3 - 6) * v3) i3 = i3 - 1
3077 end if
3078
3079 if (ii .eq. 1) then
3080 i = 9
3081 do 1002 j = i2 - 1, i2 + 1
3082 do 1001 k = i3 - 1, i3 + 1
3083 if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3084 & k .le. 10) then
3085 call dchcel(l, i, j, k, t)
3086 end if
3087 1001 continue
3088 1002 continue
3089 end if
3090
3091 if (ii .eq. 2) then
3092 i = 2
3093 do 1004 j = i2 - 1, i2 + 1
3094 do 1003 k = i3 - 1, i3 + 1
3095 if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3096 & k .le. 10) then
3097 call dchcel(l, i, j, k, t)
3098 end if
3099 1003 continue
3100 1004 continue
3101 end if
3102
3103 if (ii .eq. 3) then
3104 j = 9
3105 do 1006 i = i1 - 1, i1 + 1
3106 do 1005 k = i3 - 1, i3 + 1
3107 if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3108 & k .le. 10) then
3109 call dchcel(l, i, j, k, t)
3110 end if
3111 1005 continue
3112 1006 continue
3113 end if
3114
3115 if (ii .eq. 4) then
3116 j = 2
3117 do 1008 i = i1 - 1, i1 + 1
3118 do 1007 k = i3 - 1, i3 + 1
3119 if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3120 & k .le. 10) then
3121 call dchcel(l, i, j, k, t)
3122 end if
3123 1007 continue
3124 1008 continue
3125 end if
3126
3127 if (ii .eq. 5) then
3128 k = 9
3129 do 1010 i = i1 - 1, i1 + 1
3130 do 1009 j = i2 - 1, i2 + 1
3131 if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3132 & j .le. 10) then
3133 call dchcel(l, i, j, k, t)
3134 end if
3135 1009 continue
3136 1010 continue
3137 end if
3138
3139 if (ii .eq. 6) then
3140 k = 2
3141 do 1012 i = i1 - 1, i1 + 1
3142 do 1011 j = i2 - 1, i2 + 1
3143 if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3144 & j .le. 10) then
3145 call dchcel(l, i, j, k, t)
3146 end if
3147 1011 continue
3148 1012 continue
3149 end if
3150
3151 return
3152 end
3153
3154 subroutine dchin1(l, ii, i1, i2, i3, t)
3155c this subroutine is used to check collisions for particle inside
3156c the cube when the collision just happened is a collision including
3157c collision with wall
3158c and update the afftected particles through chkcel
3159
3160c input l,ii(specifying the direction of the wall collision),
3161c i1,i2,i3, (specifying the position of the cell
3162c we are going to check)
3163c t
3164
3165 implicit double precision (a-h, o-z)
3166 SAVE
3167
3168c itest is a flag to make sure the 111111 cell is checked only once
3169 itest = 0
3170
3171 if (ii .eq. 1) then
3172 if (i1 .eq. 1) goto 100
3173 if (i1 .eq. 2) then
3174 if (i2 .ge. 2 .and. i2 .le. 9 .and. i3 .ge. 2 .and.
3175 & i3 .le. 9) then
3176 i = 11
3177 j = 11
3178 k = 11
3179 call dchcel(l, i, j, k, t)
3180 end if
3181 goto 100
3182 end if
3183 i = i1 - 2
3184 do 1002 j = i2 - 1, i2 + 1
3185 do 1001 k = i3 - 1, i3 + 1
3186 if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3187 & k .le. 10)
3188 & call dchcel(l, i, j, k, t)
3189 1001 continue
3190 1002 continue
3191 end if
3192
3193 if (ii .eq. 2) then
3194 if (i1 .eq. 10) goto 100
3195 if (i1 .eq. 9) then
3196 if (i2 .ge. 2 .and. i2 .le. 9 .and. i3 .ge. 2 .and.
3197 & i3 .le. 9) then
3198 i = 11
3199 j = 11
3200 k = 11
3201 call dchcel(l, i, j, k, t)
3202 end if
3203 goto 100
3204 end if
3205 i = i1 + 2
3206 do 1004 j = i2 - 1, i2 + 1
3207 do 1003 k = i3 - 1, i3 + 1
3208 if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3209 & k .le. 10)
3210 & call dchcel(l, i, j, k, t)
3211 1003 continue
3212 1004 continue
3213 end if
3214
3215 if (ii .eq. 3) then
3216 if (i2 .eq. 1) goto 100
3217 if (i2 .eq. 2) then
3218 if (i1 .ge. 2 .and. i1 .le. 9 .and. i3 .ge. 2 .and.
3219 & i3 .le. 9) then
3220 i = 11
3221 j = 11
3222 k = 11
3223 call dchcel(l, i, j, k, t)
3224 end if
3225 goto 100
3226 end if
3227 j = i2 - 2
3228 do 1006 i = i1 - 1, i1 + 1
3229 do 1005 k = i3 - 1, i3 + 1
3230 if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3231 & k .le. 10)
3232 & call dchcel(l, i, j, k, t)
3233 1005 continue
3234 1006 continue
3235 end if
3236
3237 if (ii .eq. 4) then
3238 if (i2 .eq. 10) goto 100
3239 if (i2 .eq. 9) then
3240 if (i1 .ge. 2 .and. i1 .le. 9 .and. i3 .ge. 2 .and.
3241 & i3 .le. 9) then
3242 i = 11
3243 j = 11
3244 k = 11
3245 call dchcel(l, i, j, k, t)
3246 end if
3247 goto 100
3248 end if
3249 j = i2 + 2
3250 do 1008 i = i1 - 1, i1 + 1
3251 do 1007 k = i3 - 1, i3 + 1
3252 if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3253 & k .le. 10)
3254 & call dchcel(l, i, j, k, t)
3255 1007 continue
3256 1008 continue
3257 end if
3258
3259 if (ii .eq. 5) then
3260 if (i3 .eq. 1) goto 100
3261 if (i3 .eq. 2) then
3262 if (i1 .ge. 2 .and. i1 .le. 9 .and. i2 .ge. 2 .and.
3263 & i2 .le. 9) then
3264 i = 11
3265 j = 11
3266 k = 11
3267 call dchcel(l, i, j, k, t)
3268 end if
3269 goto 100
3270 end if
3271 k = i3 - 2
3272 do 1010 i = i1 - 1, i1 + 1
3273 do 1009 j = i2 - 1, i2 + 1
3274 if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3275 & j .le. 10)
3276 & call dchcel(l, i, j, k, t)
3277 1009 continue
3278 1010 continue
3279 end if
3280
3281 if (ii .eq. 6) then
3282 if (i3 .eq. 10) goto 100
3283 if (i3 .eq. 9) then
3284 if (i1 .ge. 2 .and. i1 .le. 9 .and. i2 .ge. 2 .and.
3285 & i2 .le. 9) then
3286 i = 11
3287 j = 11
3288 k = 11
3289 call dchcel(l, i, j, k, t)
3290 end if
3291 goto 100
3292 end if
3293 k = i3 + 2
3294 do 1012 i = i1 - 1, i1 + 1
3295 do 1011 j = i2 - 1, i2 + 1
3296 if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3297 & j .le. 10)
3298 & call dchcel(l, i, j, k, t)
3299 1011 continue
3300 1012 continue
3301 end if
3302
3303 100 continue
3304
3305 return
3306 end
3307
3308 subroutine dchin2(l, ii, i1, i2, i3, t)
3309c this subroutine is used to check collisions for particle inside
3310c the cube when the collision just happened is a collision including
3311c collision with wall
3312c and update the afftected particles through chkcel
3313
3314c input l,ii(specifying the direction of the wall collision),
3315c i1,i2,i3, (specifying the position of the cell
3316c we are going to check)
3317c t
3318
3319 implicit double precision (a-h, o-z)
3320 SAVE
3321
3322 if (ii .eq. 1) then
3323 i = i1 - 2
3324 if (i .le. 0) i = i + 10
3325 ia = i
3326 do 1002 j = i2 - 1, i2 + 1
3327 do 1001 k = i3 - 1, i3 + 1
3328 ib = j
3329 ic = k
3330 if (j .eq. 0) ib = 10
3331 if (j .eq. 11) ib = 1
3332 if (k .ge. 1 .and. k .le. 10) then
3333 call dchcel(l, ia, ib, ic, t)
3334 end if
3335 1001 continue
3336 1002 continue
3337 end if
3338
3339 if (ii .eq. 2) then
3340 i = i1 + 2
3341 if (i .ge. 11) i = i - 10
3342 ia = i
3343 do 1004 j = i2 - 1, i2 + 1
3344 do 1003 k = i3 - 1, i3 + 1
3345 ib = j
3346 ic = k
3347 if (j .eq. 0) ib = 10
3348 if (j .eq. 11) ib = 1
3349 if (k .ge. 1 .and. k .le. 10) then
3350 call dchcel(l, ia, ib, ic, t)
3351 end if
3352 1003 continue
3353 1004 continue
3354 end if
3355
3356 if (ii .eq. 3) then
3357 j = i2 - 2
3358 if (j .le. 0) j = j + 10
3359 ib = j
3360 do 1006 i = i1 - 1, i1 + 1
3361 do 1005 k = i3 - 1, i3 + 1
3362 ia = i
3363 ic = k
3364 if (i .eq. 0) ia = 10
3365 if (i .eq. 11) ia = 1
3366 if (k .ge. 1 .and. k .le. 10) then
3367 call dchcel(l, ia, ib, ic, t)
3368 end if
3369 1005 continue
3370 1006 continue
3371 end if
3372
3373 if (ii .eq. 4) then
3374 j = i2 + 2
3375 if (j .ge. 11) j = j - 10
3376 ib = j
3377 do 1008 i = i1 - 1, i1 + 1
3378 do 1007 k = i3 - 1, i3 + 1
3379 ia = i
3380 ic = k
3381 if (i .eq. 0) ia = 10
3382 if (i .eq. 11) ia = 1
3383 if (k .ge. 1 .and. k .le. 10) then
3384 call dchcel(l, ia, ib, ic, t)
3385 end if
3386 1007 continue
3387 1008 continue
3388 end if
3389
3390 if (ii .eq. 5) then
3391 if (i3 .eq. 2) goto 100
3392 k = i3 - 2
3393 ic = k
3394 do 1010 i = i1 - 1, i1 + 1
3395 do 1009 j = i2 - 1, i2 + 1
3396 ia = i
3397 ib = j
3398 if (i .eq. 0) ia = 10
3399 if (i .eq. 11) ia = 1
3400 if (j .eq. 0) ib = 10
3401 if (j .eq. 11) ib = 1
3402 call dchcel(l, ia, ib, ic, t)
3403 1009 continue
3404 1010 continue
3405 end if
3406
3407 if (ii .eq. 6) then
3408 if (i3 .eq. 9) goto 100
3409 k = i3 + 2
3410 ic = k
3411 do 1012 i = i1 - 1, i1 + 1
3412 do 1011 j = i2 - 1, i2 + 1
3413 ia = i
3414 ib = j
3415 if (i .eq. 0) ia = 10
3416 if (i .eq. 11) ia = 1
3417 if (j .eq. 0) ib = 10
3418 if (j .eq. 11) ib = 1
3419 call dchcel(l, ia, ib, ic, t)
3420 1011 continue
3421 1012 continue
3422 end if
3423
3424 100 continue
3425
3426 return
3427 end
3428
3429 subroutine dchin3(l, ii, i1, i2, i3, t)
3430c this subroutine is used to check collisions for particle inside
3431c the cube when the collision just happened is a collision including
3432c collision with wall
3433c and update the afftected particles through chkcel
3434
3435c input l,ii(specifying the direction of the wall collision),
3436c i1,i2,i3, (specifying the position of the cell
3437c we are going to check)
3438c t
3439
3440 implicit double precision (a-h, o-z)
3441 SAVE
3442
3443 if (ii .eq. 1) then
3444 i = i1 - 2
3445 if (i .le. 0) i = i + 10
3446 ia = i
3447 do 1002 j = i2 - 1, i2 + 1
3448 do 1001 k = i3 - 1, i3 + 1
3449 ib = j
3450 ic = k
3451 if (j .eq. 0) ib = 10
3452 if (j .eq. 11) ib = 1
3453 if (k .eq. 0) ic = 10
3454 if (k .eq. 11) ic = 1
3455 call dchcel(l, ia, ib, ic, t)
3456 1001 continue
3457 1002 continue
3458 end if
3459
3460 if (ii .eq. 2) then
3461 i = i1 + 2
3462 if (i .ge. 11) i = i - 10
3463 ia = i
3464 do 1004 j = i2 - 1, i2 + 1
3465 do 1003 k = i3 - 1, i3 + 1
3466 ib = j
3467 ic = k
3468 if (j .eq. 0) ib = 10
3469 if (j .eq. 11) ib = 1
3470 if (k .eq. 0) ic = 10
3471 if (k .eq. 11) ic = 1
3472 call dchcel(l, ia, ib, ic, t)
3473 1003 continue
3474 1004 continue
3475 end if
3476
3477 if (ii .eq. 3) then
3478 j = i2 - 2
3479 if (j .le. 0) j = j + 10
3480 ib = j
3481 do 1006 i = i1 - 1, i1 + 1
3482 do 1005 k = i3 - 1, i3 + 1
3483 ia = i
3484 ic = k
3485 if (i .eq. 0) ia = 10
3486 if (i .eq. 11) ia = 1
3487 if (k .eq. 0) ic = 10
3488 if (k .eq. 11) ic = 1
3489 call dchcel(l, ia, ib, ic, t)
3490 1005 continue
3491 1006 continue
3492 end if
3493
3494 if (ii .eq. 4) then
3495 j = i2 + 2
3496 if (j .ge. 11) j = j - 10
3497 ib = j
3498 do 1008 i = i1 - 1, i1 + 1
3499 do 1007 k = i3 - 1, i3 + 1
3500 ia = i
3501 ic = k
3502 if (i .eq. 0) ia = 10
3503 if (i .eq. 11) ia = 1
3504 if (k .eq. 0) ic = 10
3505 if (k .eq. 11) ic = 1
3506 call dchcel(l, ia, ib, ic, t)
3507 1007 continue
3508 1008 continue
3509 end if
3510
3511 if (ii .eq. 5) then
3512 k = i3 - 2
3513 if (k .le. 0) k = k + 10
3514 ic = k
3515 do 1010 i = i1 - 1, i1 + 1
3516 do 1009 j = i2 - 1, i2 + 1
3517 ia = i
3518 ib = j
3519 if (i .eq. 0) ia = 10
3520 if (i .eq. 11) ia = 1
3521 if (j .eq. 0) ib = 10
3522 if (j .eq. 11) ib = 1
3523 call dchcel(l, ia, ib, ic, t)
3524 1009 continue
3525 1010 continue
3526 end if
3527
3528 if (ii .eq. 6) then
3529 k = i3 + 2
3530 if (k .ge. 11) k = k - 10
3531 ic = k
3532 do 1012 i = i1 - 1, i1 + 1
3533 do 1011 j = i2 - 1, i2 + 1
3534 ia = i
3535 ib = j
3536 if (i .eq. 0) ia = 10
3537 if (i .eq. 11) ia = 1
3538 if (j .eq. 0) ib = 10
3539 if (j .eq. 11) ib = 1
3540 call dchcel(l, ia, ib, ic, t)
3541 1011 continue
3542 1012 continue
3543 end if
3544c
3545 return
3546 end
3547
3548 subroutine dchcel(l, i, j, k, t)
3549c this subroutine is used to recalculate next collision time for
3550c particles in the cell i,j,k if the next collision partener is l
3551
3552 implicit double precision (a-h, o-z)
3553 parameter (MAXPTN=400001)
3554
3555 common /ilist1/
3556 & iscat, jscat, next(MAXPTN), last(MAXPTN),
3557 & ictype, icsta(MAXPTN),
3558 & nic(MAXPTN), icels(MAXPTN)
3559cc SAVE /ilist1/
3560 common /ilist2/ icell, icel(10, 10, 10)
3561cc SAVE /ilist2/
3562 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3563cc SAVE /ilist5/
3564 SAVE
3565
3566 if (i .eq. 11 .or. j .eq. 11 .or. k .eq. 11) then
3567 if ( .not. (i .eq. 11 .and. j .eq. 11 .and.
3568 & k .eq. 11)) stop 'cerr'
3569 m = icell
3570 else
3571 m = icel(i, j, k)
3572 end if
3573
3574 if (m .eq. 0) return
3575 if (next(m) .eq. l) then
3576 tm = tlarge
3577 last0 = 0
3578 call reor(t, tm, m, last0)
3579 end if
3580 n = nic(m)
3581 if (n .eq. 0) return
3582 do 10 while(n .ne. m)
3583 if (next(n) .eq. l) then
3584 tm = tlarge
3585 last0 = 0
3586 call reor(t, tm, n, last0)
3587 end if
3588 n = nic(n)
3589 10 continue
3590
3591 return
3592 end
3593
3594 subroutine fixtim(l, t, tmin1, tmin, nc)
3595c this subroutine is used to compare the collision time with wall tmin1
3596c and new collision time with particles for particle l
3597c when used in ulist, input nc may be 0, which indicates no particle
3598c collisions happen before wall collision, of course, then tmin=tmin1
3599
3600 implicit double precision (a-h, o-z)
3601 parameter (MAXPTN=400001)
3602
3603 common /ilist1/
3604 & iscat, jscat, next(MAXPTN), last(MAXPTN),
3605 & ictype, icsta(MAXPTN),
3606 & nic(MAXPTN), icels(MAXPTN)
3607cc SAVE /ilist1/
3608 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3609cc SAVE /ilist5/
3610 SAVE
3611
3612 t=t
3613 k = nc
3614 if (tmin .lt. tmin1) then
3615 ot(l) = tmin
3616 if (ct(l) .lt. tmin1) then
3617 icsta(l) = 0
3618 else
3619 icsta(l) = icsta(l) + 10
3620 end if
3621 next(l) = k
3622 else if (tmin .eq. tmin1) then
3623 ot(l) = tmin
3624 if (nc .eq. 0) then
3625 next(l) = 0
3626 else
3627 icsta(l) = icsta(l) + 10
3628 next(l) = k
3629 end if
3630 else
3631 ot(l) = tmin1
3632 next(l) = 0
3633 end if
3634
3635 return
3636 end
3637
3638 subroutine ud2(i, j, t, tmin, nc)
3639c this subroutine is used to update next(i), ct(i), ot(i),
3640c and get tmin, nc for j
3641
3642 implicit double precision (a-h, o-z)
3643 parameter (MAXPTN=400001)
3644
3645 common /para5/ iconfg, iordsc
3646cc SAVE /para5/
3647 common /aurec1/ jxa, jya, jza
3648cc SAVE /aurec1/
3649 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
3650cc SAVE /aurec2/
3651 common /ilist1/
3652 & iscat, jscat, next(MAXPTN), last(MAXPTN),
3653 & ictype, icsta(MAXPTN),
3654 & nic(MAXPTN), icels(MAXPTN)
3655cc SAVE /ilist1/
3656 common /ilist3/ size1, size2, size3, v1, v2, v3, size
3657cc SAVE /ilist3/
3658 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3659cc SAVE /ilist5/
3660 SAVE
3661
3662 logical allok
3663
3664 call isco(i, j, allok, tm, t1, t2)
3665
3666 if (allok) then
3667c tm eq tmin, change nc to make sure fixtime get the collision with both
3668c wall and particle
3669
3670 if (tm .lt. tmin) then
3671 tmin = tm
3672 ct(j) = t2
3673 nc = i
3674 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
3675 dgxa(j) = jxa * 10d0 * size1
3676 dgya(j) = jya * 10d0 * size2
3677 if (iconfg .eq. 5) then
3678 dgza(j) = jza * 10d0 * size3
3679 end if
3680 end if
3681 end if
3682
3683 if (tm .le. ot(i)) then
3684 ct(i) = t1
3685 icels0 = icels(i)
3686 i1 = icels0 / 10000
3687 i2 = (icels0 - i1 * 10000) / 100
3688 i3 = icels0 - i1 * 10000 - i2 * 100
3689 call wallc(i, i1, i2, i3, t, tmin1)
3690 call fixtim(i, t, tmin1, tm, j)
3691 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
3692 dgxa(i) = - jxa * 10d0 * size1
3693 dgya(i) = - jya * 10d0 * size2
3694 if (iconfg .eq. 5) then
3695 dgza(i) = - jza * 10d0 * size3
3696 end if
3697 end if
3698 end if
3699
3700 if (tm .gt. ot(i) .and. next(i) .eq. j) then
3701 ct(i) = t1
3702 call reor(t, tm, i, j)
3703 end if
3704
3705 else if (next(i) .eq. j) then
3706
3707 tm = tlarge
3708
3709 call reor(t, tm, i, j)
3710
3711 end if
3712
3713 return
3714 end
3715
3716 subroutine isco(i, j, allok, tm, t1, t2)
3717
3718 implicit double precision (a-h, o-z)
3719
3720 common /para5/ iconfg, iordsc
3721cc SAVE /para5/
3722 SAVE
3723
3724 logical allok
3725
3726 iorder = iordsc / 10
3727 if (iconfg .eq. 1) then
3728 if (iorder .eq. 1) then
3729 call isco1(i, j, allok, tm, t1, t2)
3730 else if (iorder .eq. 2) then
3731 call isco2(i, j, allok, tm, t1, t2)
3732 else if (iorder .eq. 3) then
3733 call isco3(i, j, allok, tm, t1, t2)
3734 end if
3735 else if (iconfg .eq. 2 .or. iconfg .eq. 4) then
3736 if (iorder .eq. 1) then
3737 call isco4(i, j, allok, tm, t1, t2)
3738 else if (iorder .eq. 2) then
3739 call isco5(i, j, allok, tm, t1, t2)
3740 else if (iorder .eq. 3) then
3741 call isco6(i, j, allok, tm, t1, t2)
3742 end if
3743 else if (iconfg .eq. 3) then
3744 if (iorder .eq. 1) then
3745 call isco7(i, j, allok, tm, t1, t2)
3746 else if (iorder .eq. 2) then
3747 call isco8(i, j, allok, tm, t1, t2)
3748 else if (iorder .eq. 3) then
3749 call isco9(i, j, allok, tm, t1, t2)
3750 end if
3751 else if (iconfg .eq. 5) then
3752 if (iorder .eq. 1) then
3753 call isco10(i, j, allok, tm, t1, t2)
3754 else if (iorder .eq. 2) then
3755 call isco11(i, j, allok, tm, t1, t2)
3756 else if (iorder .eq. 3) then
3757 call isco12(i, j, allok, tm, t1, t2)
3758 end if
3759 end if
3760
3761 return
3762 end
3763
3764 subroutine isco1(i, j, allok, tm, t1, t2)
3765c this subroutine is used to decide whether there is a collision between
3766c particle i and j, if there is one allok=1, and tm gives the
3767c collision time, t1 the collision time for i,
3768c t2 the collision time for j
3769
3770 implicit double precision (a-h, o-z)
3771 parameter (MAXPTN=400001)
3772
3773 common /para2/ xmp, xmu, alpha, rscut2, cutof2
3774cc SAVE /para2/
3775 common /para5/ iconfg, iordsc
3776cc SAVE /para5/
3777 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3778 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3779 & xmass(MAXPTN), ityp(MAXPTN)
3780cc SAVE /prec2/
3781 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3782cc SAVE /prec4/
3783 common /ilist1/
3784 & iscat, jscat, next(MAXPTN), last(MAXPTN),
3785 & ictype, icsta(MAXPTN),
3786 & nic(MAXPTN), icels(MAXPTN)
3787cc SAVE /ilist1/
3788 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3789cc SAVE /ilist5/
3790 SAVE
3791
3792 logical allok
3793
3794c preventing consecutive collisions
3795 allok = last(i) .ne. j .or. last(j) .ne. i
3796
3797c set up numbers for later calculations
3798 i1 = i
3799 i2 = j
3800
3801 p4 = ft(i2) - ft(i1)
3802 p1 = gx(i2) - gx(i1)
3803 p2 = gy(i2) - gy(i1)
3804 p3 = gz(i2) - gz(i1)
3805
3806 q4 = e(i1)
3807 q1 = px(i1)
3808 q2 = py(i1)
3809 q3 = pz(i1)
3810
3811 r4 = e(i2)
3812 r1 = px(i2)
3813 r2 = py(i2)
3814 r3 = pz(i2)
3815
3816 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
3817 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
3818 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
3819 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
3820 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
3821 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
3822
3823c make sure particle 2 formed early
3824 h = a + b
3825 if (h .gt. 0d0) then
3826 g = a
3827 a = -b
3828 b = -g
3829
3830 g = c
3831 c = d
3832 d = g
3833
3834 i1 = j
3835 i2 = i
3836 end if
3837
3838c check the approaching criteria
3839 if (allok) then
3840
3841 vp = a * d - b * ee
3842
3843 allok = allok .and. vp .lt. 0d0
3844
3845 end if
3846
3847c check the closest approach distance criteria
3848 if (allok) then
3849
3850 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
3851 & (ee ** 2 - c * d)
3852
3853 allok = allok .and. dm2 .lt. cutof2
3854
3855 end if
3856
3857c check the time criteria
3858 if (allok) then
3859
3860 tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
3861 tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
3862 tm = 0.5d0 * (tc1 + tc2)
3863
3864 allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
3865
3866 end if
3867
3868c check rts cut
3869 if (allok) then
3870
3871 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
3872 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
3873
3874 allok = allok .and. rts2 .gt. rscut2
3875 end if
3876
3877 if (.not. allok) then
3878 tm = tlarge
3879 t1 = tlarge
3880 t2 = tlarge
3881 else if (h .gt. 0d0) then
3882 t1 = tm
3883 t2 = tm
3884 else
3885 t1 = tm
3886 t2 = tm
3887 end if
3888
3889 return
3890 end
3891
3892 subroutine isco2(i, j, allok, tm, t1, t2)
3893c this subroutine is used to decide whether there is a collision between
3894c particle i and j, if there is one allok=1, and tm gives the
3895c collision time, t1 the collision time for i,
3896c t2 the collision time for j
3897
3898 implicit double precision (a-h, o-z)
3899 parameter (MAXPTN=400001)
3900
3901 common /para2/ xmp, xmu, alpha, rscut2, cutof2
3902cc SAVE /para2/
3903 common /para5/ iconfg, iordsc
3904cc SAVE /para5/
3905 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3906 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3907 & xmass(MAXPTN), ityp(MAXPTN)
3908cc SAVE /prec2/
3909 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3910cc SAVE /prec4/
3911 common /ilist1/
3912 & iscat, jscat, next(MAXPTN), last(MAXPTN),
3913 & ictype, icsta(MAXPTN),
3914 & nic(MAXPTN), icels(MAXPTN)
3915cc SAVE /ilist1/
3916 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3917cc SAVE /ilist5/
3918 SAVE
3919
3920 logical allok
3921
3922c preventing consecutive collisions
3923 allok = last(i) .ne. j .or. last(j) .ne. i
3924
3925c set up numbers for later calculations
3926 i1 = i
3927 i2 = j
3928
3929 p4 = ft(i2) - ft(i1)
3930 p1 = gx(i2) - gx(i1)
3931 p2 = gy(i2) - gy(i1)
3932 p3 = gz(i2) - gz(i1)
3933
3934 q4 = e(i1)
3935 q1 = px(i1)
3936 q2 = py(i1)
3937 q3 = pz(i1)
3938
3939 r4 = e(i2)
3940 r1 = px(i2)
3941 r2 = py(i2)
3942 r3 = pz(i2)
3943
3944 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
3945 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
3946 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
3947 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
3948 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
3949 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
3950
3951c make sure particle 2 formed early
3952 h = a + b
3953 if (h .gt. 0d0) then
3954 g = a
3955 a = -b
3956 b = -g
3957
3958 g = c
3959 c = d
3960 d = g
3961
3962 i1 = j
3963 i2 = i
3964 end if
3965
3966c check the approaching criteria
3967 if (allok) then
3968
3969 vp = a * d - b * ee
3970
3971 allok = allok .and. vp .lt. 0d0
3972
3973 end if
3974
3975c check the closest approach distance criteria
3976 if (allok) then
3977
3978 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
3979 & (ee ** 2 - c * d)
3980
3981 allok = allok .and. dm2 .lt. cutof2
3982
3983 end if
3984
3985c check the time criteria
3986 if (allok) then
3987
3988 tc1 = ft(i1) - e(i1) * (a * d - b * ee)/(ee ** 2 - c * d)
3989 tc2 = ft(i2) + e(i2) * (b * c - a * ee)/(ee ** 2 - c * d)
3990 if (iordsc .eq. 20) then
3991 tm = min(tc1, tc2)
3992 else if (iordsc .eq. 21) then
3993 tm = 0.5d0 * (tc1 + tc2)
3994 else
3995 tm = max(tc1, tc2)
3996 end if
3997
3998 allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
3999
4000 end if
4001
4002c check rts cut
4003 if (allok) then
4004
4005 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4006 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
4007
4008 allok = allok .and. rts2 .gt. rscut2
4009 end if
4010
4011 if (.not. allok) then
4012 tm = tlarge
4013 t1 = tlarge
4014 t2 = tlarge
4015 else if (h .gt. 0d0) then
4016 t1 = tc2
4017 t2 = tc1
4018 else
4019 t1 = tc1
4020 t2 = tc2
4021 end if
4022
4023 return
4024 end
4025
4026 subroutine isco3(i, j, allok, tm, t1, t2)
4027c this subroutine is used to decide whether there is a collision between
4028c particle i and j, if there is one allok=1, and tm gives the
4029c collision time, t1 the collision time for i,
4030c t2 the collision time for j
4031
4032 implicit double precision (a-h, o-z)
4033 parameter (MAXPTN=400001)
4034
4035 common /para2/ xmp, xmu, alpha, rscut2, cutof2
4036cc SAVE /para2/
4037 common /para5/ iconfg, iordsc
4038cc SAVE /para5/
4039 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4040 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4041 & xmass(MAXPTN), ityp(MAXPTN)
4042cc SAVE /prec2/
4043 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4044cc SAVE /prec4/
4045 common /ilist1/
4046 & iscat, jscat, next(MAXPTN), last(MAXPTN),
4047 & ictype, icsta(MAXPTN),
4048 & nic(MAXPTN), icels(MAXPTN)
4049cc SAVE /ilist1/
4050 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4051cc SAVE /ilist5/
4052 SAVE
4053
4054 logical allok
4055
4056c preventing consecutive collisions
4057 allok = last(i) .ne. j .or. last(j) .ne. i
4058
4059 if (ft(i) .ge. ft(j)) then
4060 i1 = j
4061 i2 = i
4062 else
4063 i1 = i
4064 i2 = j
4065 end if
4066
4067 if (allok) then
4068
4069 t1 = ft(i1)
4070 vx1 = vx(i1)
4071 vy1 = vy(i1)
4072 vz1 = vz(i1)
4073
4074 t2 = ft(i2)
4075
4076 dvx = vx(i2) - vx1
4077 dvy = vy(i2) - vy1
4078 dvz = vz(i2) - vz1
4079
4080 dt = t2 - t1
4081
4082 dx = gx(i2) - gx(i1) - vx1 * dt
4083 dy = gy(i2) - gy(i1) - vy1 * dt
4084 dz = gz(i2) - gz(i1) - vz1 * dt
4085
4086 vp = dvx * dx + dvy * dy + dvz * dz
4087
4088 allok = allok .and. vp .lt. 0d0
4089
4090 end if
4091
4092 if (allok) then
4093
4094 v2= dvx * dvx + dvy * dvy + dvz * dvz
4095
4096 if (v2 .eq. 0d0) then
4097 tm = tlarge
4098 else
4099 tm = t2 - vp / v2
4100 end if
4101
4102c note now tm is the absolute time
4103
4104 allok = allok .and. tm .gt. t1 .and. tm .gt. t2
4105
4106 end if
4107
4108 if (allok) then
4109
4110 dgx = dx - dvx * t2
4111 dgy = dy - dvy * t2
4112 dgz = dz - dvz * t2
4113
4114 dm2 = - v2 * tm ** 2 + dgx * dgx + dgy * dgy + dgz * dgz
4115
4116 allok = allok .and. dm2 .lt. cutof2
4117
4118 end if
4119
4120 if (allok) then
4121
4122 e1 = e(i1)
4123 px1 = px(i1)
4124 py1 = py(i1)
4125 pz1 = pz(i1)
4126 e2 = e(i2)
4127 px2 = px(i2)
4128 py2 = py(i2)
4129 pz2 = pz(i2)
4130
4131 rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
4132 & - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
4133
4134 allok = allok .and. rts2 .gt. rscut2
4135 end if
4136
4137 if (.not. allok) then
4138 tm = tlarge
4139 t1 = tlarge
4140 t2 = tlarge
4141 else
4142 t1 = tm
4143 t2 = tm
4144 end if
4145
4146 return
4147 end
4148
4149 subroutine isco4(i, j, allok, tm, t1, t2)
4150c this subroutine is used to decide whether there is a collision between
4151c particle i and j, if there is one allok=1, and tm gives the
4152c collision time, t1 the collision time for i,
4153c t2 the collision time for j
4154
4155 implicit double precision (a-h, o-z)
4156 parameter (MAXPTN=400001)
4157
4158 common /para2/ xmp, xmu, alpha, rscut2, cutof2
4159cc SAVE /para2/
4160 common /para5/ iconfg, iordsc
4161cc SAVE /para5/
4162 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4163 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4164 & xmass(MAXPTN), ityp(MAXPTN)
4165cc SAVE /prec2/
4166 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4167cc SAVE /prec4/
4168 common /ilist1/
4169 & iscat, jscat, next(MAXPTN), last(MAXPTN),
4170 & ictype, icsta(MAXPTN),
4171 & nic(MAXPTN), icels(MAXPTN)
4172cc SAVE /ilist1/
4173 common /ilist3/ size1, size2, size3, v1, v2, v3, size
4174cc SAVE /ilist3/
4175 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4176cc SAVE /ilist5/
4177 SAVE
4178
4179 logical allok
4180
4181c preventing consecutive collisions
4182 allok = last(i) .ne. j .or. last(j) .ne. i
4183
4184c set up numbers for later calculations
4185
4186 icels1 = icels(i)
4187 ii1 = icels1 / 10000
4188 jj1 = (icels1 - ii1 * 10000) / 100
4189 kk1 = icels1 - ii1 * 10000 - jj1 * 100
4190 icels2 = icels(j)
4191 ii2 = icels2 / 10000
4192 jj2 = (icels2 - ii2 * 10000) / 100
4193 kk2 = icels2 - ii2 * 10000 - jj2 * 100
4194
4195 i1 = i
4196 i2 = j
4197
4198 p4 = ft(i2) - ft(i1)
4199 p1 = gx(i2) - gx(i1)
4200 p2 = gy(i2) - gy(i1)
4201 p3 = gz(i2) - gz(i1)
4202
4203 if (ii1 - ii2 .gt. 5) then
4204 p1 = p1 + 10d0 * size1
4205 else if (ii1 - ii2 .lt. -5) then
4206 p1 = p1 - 10d0 * size1
4207 end if
4208 if (jj1 - jj2 .gt. 5) then
4209 p2 = p2 + 10d0 * size2
4210 else if (jj1 - jj2 .lt. -5) then
4211 p2 = p2 - 10d0 * size2
4212 end if
4213 if (kk1 - kk2 .gt. 5) then
4214 p3 = p3 + 10d0 * size3
4215 else if (kk1 - kk2 .lt. -5) then
4216 p3 = p3 - 10d0 * size3
4217 end if
4218
4219 q4 = e(i1)
4220 q1 = px(i1)
4221 q2 = py(i1)
4222 q3 = pz(i1)
4223
4224 r4 = e(i2)
4225 r1 = px(i2)
4226 r2 = py(i2)
4227 r3 = pz(i2)
4228
4229 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4230 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4231 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4232 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4233 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4234 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4235
4236c make sure particle 2 formed early
4237 h = a + b
4238 if (h .gt. 0d0) then
4239 g = a
4240 a = -b
4241 b = -g
4242
4243 g = c
4244 c = d
4245 d = g
4246
4247 i1 = j
4248 i2 = i
4249 end if
4250
4251c check the approaching criteria
4252 if (allok) then
4253
4254 vp = a * d - b * ee
4255
4256 allok = allok .and. vp .lt. 0d0
4257
4258 end if
4259
4260c check the closest approach distance criteria
4261 if (allok) then
4262
4263 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4264 & (ee ** 2 - c * d)
4265
4266 allok = allok .and. dm2 .lt. cutof2
4267
4268 end if
4269
4270c check the time criteria
4271 if (allok) then
4272
4273 tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4274 tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4275 tm = 0.5d0 * (tc1 + tc2)
4276
4277 allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
4278
4279 end if
4280
4281c check rts cut
4282 if (allok) then
4283
4284 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4285 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
4286
4287 allok = allok .and. rts2 .gt. rscut2
4288 end if
4289
4290 if (.not. allok) then
4291 tm = tlarge
4292 t1 = tlarge
4293 t2 = tlarge
4294 else if (h .gt. 0d0) then
4295 t1 = tm
4296 t2 = tm
4297 else
4298 t1 = tm
4299 t2 = tm
4300 end if
4301
4302 return
4303 end
4304
4305 subroutine isco5(i, j, allok, tm, t1, t2)
4306c this subroutine is used to decide whether there is a collision between
4307c particle i and j, if there is one allok=1, and tm gives the
4308c collision time, t1 the collision time for i,
4309c t2 the collision time for j
4310
4311 implicit double precision (a-h, o-z)
4312 parameter (MAXPTN=400001)
4313
4314 common /para2/ xmp, xmu, alpha, rscut2, cutof2
4315cc SAVE /para2/
4316 common /para5/ iconfg, iordsc
4317cc SAVE /para5/
4318 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4319 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4320 & xmass(MAXPTN), ityp(MAXPTN)
4321cc SAVE /prec2/
4322 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4323cc SAVE /prec4/
4324 common /ilist1/
4325 & iscat, jscat, next(MAXPTN), last(MAXPTN),
4326 & ictype, icsta(MAXPTN),
4327 & nic(MAXPTN), icels(MAXPTN)
4328cc SAVE /ilist1/
4329 common /ilist3/ size1, size2, size3, v1, v2, v3, size
4330cc SAVE /ilist3/
4331 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4332cc SAVE /ilist5/
4333 SAVE
4334
4335 logical allok
4336
4337c preventing consecutive collisions
4338 allok = last(i) .ne. j .or. last(j) .ne. i
4339
4340c set up numbers for later calculations
4341
4342 icels1 = icels(i)
4343 ii1 = icels1 / 10000
4344 jj1 = (icels1 - ii1 * 10000) / 100
4345 kk1 = icels1 - ii1 * 10000 - jj1 * 100
4346 icels2 = icels(j)
4347 ii2 = icels2 / 10000
4348 jj2 = (icels2 - ii2 * 10000) / 100
4349 kk2 = icels2 - ii2 * 10000 - jj2 * 100
4350
4351 i1 = i
4352 i2 = j
4353
4354 p4 = ft(i2) - ft(i1)
4355 p1 = gx(i2) - gx(i1)
4356 p2 = gy(i2) - gy(i1)
4357 p3 = gz(i2) - gz(i1)
4358
4359 if (ii1 - ii2 .gt. 5) then
4360 p1 = p1 + 10d0 * size1
4361 else if (ii1 - ii2 .lt. -5) then
4362 p1 = p1 - 10d0 * size1
4363 end if
4364 if (jj1 - jj2 .gt. 5) then
4365 p2 = p2 + 10d0 * size2
4366 else if (jj1 - jj2 .lt. -5) then
4367 p2 = p2 - 10d0 * size2
4368 end if
4369 if (kk1 - kk2 .gt. 5) then
4370 p3 = p3 + 10d0 * size3
4371 else if (kk1 - kk2 .lt. -5) then
4372 p3 = p3 - 10d0 * size3
4373 end if
4374
4375 q4 = e(i1)
4376 q1 = px(i1)
4377 q2 = py(i1)
4378 q3 = pz(i1)
4379
4380 r4 = e(i2)
4381 r1 = px(i2)
4382 r2 = py(i2)
4383 r3 = pz(i2)
4384
4385 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4386 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4387 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4388 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4389 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4390 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4391
4392c make sure particle 2 formed early
4393 h = a + b
4394 if (h .gt. 0d0) then
4395 g = a
4396 a = -b
4397 b = -g
4398
4399 g = c
4400 c = d
4401 d = g
4402
4403 i1 = j
4404 i2 = i
4405 end if
4406
4407c check the approaching criteria
4408 if (allok) then
4409
4410 vp = a * d - b * ee
4411
4412 allok = allok .and. vp .lt. 0d0
4413
4414 end if
4415
4416c check the closest approach distance criteria
4417 if (allok) then
4418
4419 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4420 & (ee ** 2 - c * d)
4421
4422 allok = allok .and. dm2 .lt. cutof2
4423
4424 end if
4425
4426c check the time criteria
4427 if (allok) then
4428
4429 tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4430 tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4431 if (iordsc .eq. 20) then
4432 tm = min(tc1, tc2)
4433 else if (iordsc .eq. 21) then
4434 tm = 0.5d0 * (tc1 + tc2)
4435 else
4436 tm = max(tc1, tc2)
4437 end if
4438
4439 allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
4440
4441 end if
4442
4443c check rts cut
4444 if (allok) then
4445
4446 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4447 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
4448
4449 allok = allok .and. rts2 .gt. rscut2
4450 end if
4451
4452 if (.not. allok) then
4453 tm = tlarge
4454 t1 = tlarge
4455 t2 = tlarge
4456 else if (h .gt. 0d0) then
4457 t1 = tc2
4458 t2 = tc1
4459 else
4460 t1 = tc1
4461 t2 = tc2
4462 end if
4463
4464 return
4465 end
4466
4467 subroutine isco6(i, j, allok, tm, t1, t2)
4468c this subroutine is used to decide whether there is a collision between
4469c particle i and j, if there is one allok=1, and tm gives the
4470c collision time, t1 the collision time for i,
4471c t2 the collision time for j
4472
4473 implicit double precision (a-h, o-z)
4474 parameter (MAXPTN=400001)
4475
4476 common /para2/ xmp, xmu, alpha, rscut2, cutof2
4477cc SAVE /para2/
4478 common /para5/ iconfg, iordsc
4479cc SAVE /para5/
4480 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4481 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4482 & xmass(MAXPTN), ityp(MAXPTN)
4483cc SAVE /prec2/
4484 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4485cc SAVE /prec4/
4486 common /ilist1/
4487 & iscat, jscat, next(MAXPTN), last(MAXPTN),
4488 & ictype, icsta(MAXPTN),
4489 & nic(MAXPTN), icels(MAXPTN)
4490cc SAVE /ilist1/
4491 common /ilist3/ size1, size2, size3, v1, v2, v3, size
4492cc SAVE /ilist3/
4493 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4494cc SAVE /ilist5/
4495 SAVE
4496
4497 logical allok
4498
4499c preventing consecutive collisions
4500 allok = last(i) .ne. j .or. last(j) .ne. i
4501
4502 if (ft(i) .ge. ft(j)) then
4503 i1 = j
4504 i2 = i
4505 else
4506 i1 = i
4507 i2 = j
4508 end if
4509
4510 icels1 = icels(i1)
4511 ii1 = icels1 / 10000
4512 jj1 = (icels1 - ii1 * 10000) / 100
4513 kk1 = icels1 - ii1 * 10000 - jj1 * 100
4514 icels2 = icels(i2)
4515 ii2 = icels2 / 10000
4516 jj2 = (icels2 - ii2 * 10000) / 100
4517 kk2 = icels2 - ii2 * 10000 - jj2 * 100
4518
4519 if (allok) then
4520
4521 t1 = ft(i1)
4522 vx1 = vx(i1)
4523 vy1 = vy(i1)
4524 vz1 = vz(i1)
4525
4526 t2 = ft(i2)
4527
4528 dvx = vx(i2) - vx1
4529 dvy = vy(i2) - vy1
4530 dvz = vz(i2) - vz1
4531
4532 dt = t2 - t1
4533
4534 dx = gx(i2) - gx(i1) - vx1 * dt
4535 dy = gy(i2) - gy(i1) - vy1 * dt
4536 dz = gz(i2) - gz(i1) - vz1 * dt
4537
4538 if (ii1 - ii2 .gt. 5) then
4539 dx = dx + 10d0 * size1
4540 else if (ii1 - ii2 .lt. - 5) then
4541 dx = dx - 10d0 * size1
4542 end if
4543
4544 if (jj1 - jj2 .gt. 5) then
4545 dy = dy + 10d0 * size2
4546 else if (jj1 - jj2 .lt. - 5) then
4547 dy = dy - 10d0 * size2
4548 end if
4549
4550 if (kk1 - kk2 .gt. 5) then
4551 dz = dz + 10d0 * size3
4552 else if (kk1 - kk2 .lt. -5) then
4553 dz = dz - 10d0 * size3
4554 end if
4555
4556 vp = dvx * dx + dvy * dy + dvz * dz
4557
4558 allok = allok .and. vp .lt. 0d0
4559
4560 end if
4561
4562 if (allok) then
4563
4564 v2p = dvx * dvx + dvy * dvy + dvz * dvz
4565
4566 if (v2p .eq. 0d0) then
4567 tm = tlarge
4568 else
4569 tm = t2 - vp / v2p
4570 end if
4571
4572c note now tm is the absolute time
4573
4574 allok = allok .and. tm .gt. t1 .and. tm .gt. t2
4575
4576 end if
4577
4578 if (allok) then
4579
4580 dgx = dx - dvx * t2
4581 dgy = dy - dvy * t2
4582 dgz = dz - dvz * t2
4583
4584 dm2 = - v2p * tm ** 2 + dgx * dgx + dgy * dgy + dgz * dgz
4585
4586 allok = allok .and. dm2 .lt. cutof2
4587
4588 end if
4589
4590 if (allok) then
4591
4592 e1 = e(i1)
4593 px1 = px(i1)
4594 py1 = py(i1)
4595 pz1 = pz(i1)
4596 e2 = e(i2)
4597 px2 = px(i2)
4598 py2 = py(i2)
4599 pz2 = pz(i2)
4600
4601 rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
4602 & - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
4603
4604 allok = allok .and. rts2 .gt. rscut2
4605 end if
4606
4607 if (.not. allok) then
4608 tm = tlarge
4609 t1 = tlarge
4610 t2 = tlarge
4611 else
4612 t1 = tm
4613 t2 = tm
4614 end if
4615
4616 return
4617 end
4618
4619 subroutine isco7(i, j, allok, tm, t1, t2)
4620c this subroutine is used to decide whether there is a collision between
4621c particle i and j, if there is one allok=1, and tm gives the
4622c collision time, t1 the collision time for i,
4623c t2 the collision time for j
4624
4625 implicit double precision (a-h, o-z)
4626 parameter (MAXPTN=400001)
4627
4628 common /para2/ xmp, xmu, alpha, rscut2, cutof2
4629cc SAVE /para2/
4630 common /para5/ iconfg, iordsc
4631cc SAVE /para5/
4632 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4633 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4634 & xmass(MAXPTN), ityp(MAXPTN)
4635cc SAVE /prec2/
4636 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4637cc SAVE /prec4/
4638 common /aurec1/ jxa, jya, jza
4639cc SAVE /aurec1/
4640 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4641cc SAVE /aurec2/
4642 common /ilist1/
4643 & iscat, jscat, next(MAXPTN), last(MAXPTN),
4644 & ictype, icsta(MAXPTN),
4645 & nic(MAXPTN), icels(MAXPTN)
4646cc SAVE /ilist1/
4647 common /ilist3/ size1, size2, size3, v1, v2, v3, size
4648cc SAVE /ilist3/
4649 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4650cc SAVE /ilist5/
4651 SAVE
4652
4653 logical allok, allokp
4654
4655c preventing consecutive collisions
4656 allok = last(i) .ne. j .or. last(j) .ne. i
4657
4658c set up numbers for later calculations
4659
4660 tm = tlarge
4661
4662 if (allok) then
4663 do 1000 ii = - 1, 1
4664 do 2000 jj = - 1, 1
4665
4666 allokp = .true.
4667
4668 i1 = i
4669 i2 = j
4670
4671 p4 = ft(j) - ft(i)
4672 p1 = gx(j) - gx(i)
4673 p2 = gy(j) - gy(i)
4674 p3 = gz(j) - gz(i)
4675
4676 p1 = p1 + ii * 10d0 * size1
4677 p2 = p2 + jj * 10d0 * size2
4678
4679 q4 = e(i)
4680 q1 = px(i)
4681 q2 = py(i)
4682 q3 = pz(i)
4683
4684 r4 = e(j)
4685 r1 = px(j)
4686 r2 = py(j)
4687 r3 = pz(j)
4688
4689 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4690 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4691 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4692 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4693 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4694 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4695
4696c make sure particle 2 formed early
4697 h = a + b
4698 if (h .gt. 0d0) then
4699 g = a
4700 a = -b
4701 b = -g
4702 g = c
4703 c = d
4704 d = g
4705 i1 = j
4706 i2 = i
4707 end if
4708
4709c check the approaching criteria
4710 if (allokp) then
4711 vp = a * d - b * ee
4712 allokp = allokp .and. vp .lt. 0d0
4713 end if
4714
4715c check the closest approach distance criteria
4716 if (allokp) then
4717 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4718 & (ee ** 2 - c * d)
4719 allokp = allokp .and. dm2 .lt. cutof2
4720 end if
4721
4722c check the time criteria
4723 if (allokp) then
4724 tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4725 tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4726 tmp = 0.5d0 * (tc1 + tc2)
4727 allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
4728 end if
4729
4730 if (allokp .and. tmp .lt. tm) then
4731 tm = tmp
4732 jxa = ii
4733 jya = jj
4734cd dgxa(j) = ii * 10d0 * size1
4735cd dgya(j) = jj * 10d0 * size2
4736cd dgxa(i) = - dgxa(j)
4737cd dgya(i) = - dgya(j)
4738 end if
4739
4740 2000 continue
4741 1000 continue
4742
4743 if (tm .eq. tlarge) then
4744 allok = .false.
4745 end if
4746
4747 end if
4748
4749c check rts cut
4750 if (allok) then
4751
4752 q4 = e(i1)
4753 q1 = px(i1)
4754 q2 = py(i1)
4755 q3 = pz(i1)
4756
4757 r4 = e(i2)
4758 r1 = px(i2)
4759 r2 = py(i2)
4760 r3 = pz(i2)
4761
4762 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4763 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
4764
4765 allok = allok .and. rts2 .gt. rscut2
4766 end if
4767
4768 if (.not. allok) then
4769 tm = tlarge
4770 t1 = tlarge
4771 t2 = tlarge
4772 else if (h .gt. 0d0) then
4773 t1 = tm
4774 t2 = tm
4775 else
4776 t1 = tm
4777 t2 = tm
4778 end if
4779
4780 return
4781 end
4782
4783 subroutine isco8(i, j, allok, tm, t1, t2)
4784c this subroutine is used to decide whether there is a collision between
4785c particle i and j, if there is one allok=1, and tm gives the
4786c collision time, t1 the collision time for i,
4787c t2 the collision time for j
4788
4789 implicit double precision (a-h, o-z)
4790 parameter (MAXPTN=400001)
4791
4792 common /para2/ xmp, xmu, alpha, rscut2, cutof2
4793cc SAVE /para2/
4794 common /para5/ iconfg, iordsc
4795cc SAVE /para5/
4796 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4797 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4798 & xmass(MAXPTN), ityp(MAXPTN)
4799cc SAVE /prec2/
4800 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4801cc SAVE /prec4/
4802 common /aurec1/ jxa, jya, jza
4803cc SAVE /aurec1/
4804 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4805cc SAVE /aurec2/
4806 common /ilist1/
4807 & iscat, jscat, next(MAXPTN), last(MAXPTN),
4808 & ictype, icsta(MAXPTN),
4809 & nic(MAXPTN), icels(MAXPTN)
4810cc SAVE /ilist1/
4811 common /ilist3/ size1, size2, size3, v1, v2, v3, size
4812cc SAVE /ilist3/
4813 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4814cc SAVE /ilist5/
4815 SAVE
4816
4817 logical allok, allokp
4818
4819c preventing consecutive collisions
4820 allok = last(i) .ne. j .or. last(j) .ne. i
4821
4822c set up numbers for later calculations
4823
4824 tm = tlarge
4825
4826 if (allok) then
4827 do 1000 ii = - 1, 1
4828 do 2000 jj = - 1, 1
4829
4830 allokp = .true.
4831
4832 i1 = i
4833 i2 = j
4834
4835 p4 = ft(j) - ft(i)
4836 p1 = gx(j) - gx(i)
4837 p2 = gy(j) - gy(i)
4838 p3 = gz(j) - gz(i)
4839
4840 p1 = p1 + ii * 10d0 * size1
4841 p2 = p2 + jj * 10d0 * size2
4842
4843 q4 = e(i)
4844 q1 = px(i)
4845 q2 = py(i)
4846 q3 = pz(i)
4847
4848 r4 = e(j)
4849 r1 = px(j)
4850 r2 = py(j)
4851 r3 = pz(j)
4852
4853 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4854 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4855 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4856 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4857 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4858 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4859
4860c make sure particle 2 formed early
4861 h = a + b
4862 if (h .gt. 0d0) then
4863 g = a
4864 a = -b
4865 b = -g
4866 g = c
4867 c = d
4868 d = g
4869 i1 = j
4870 i2 = i
4871 end if
4872
4873c check the approaching criteria
4874 if (allokp) then
4875 vp = a * d - b * ee
4876 allokp = allokp .and. vp .lt. 0d0
4877 end if
4878
4879c check the closest approach distance criteria
4880 if (allokp) then
4881 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4882 & (ee ** 2 - c * d)
4883 allokp = allokp .and. dm2 .lt. cutof2
4884 end if
4885
4886c check the time criteria
4887 if (allokp) then
4888 tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4889 tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4890 if (iordsc .eq. 20) then
4891 tmp = min(tc1, tc2)
4892 else if (iordsc .eq. 21) then
4893 tmp = 0.5d0 * (tc1 + tc2)
4894 else
4895 tmp = max(tc1, tc2)
4896 end if
4897 allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
4898 end if
4899
4900 if (allokp .and. tmp .lt. tm) then
4901 tm = tmp
4902 jxa = ii
4903 jya = jj
4904 ha = h
4905 tc1a = tc1
4906 tc2a = tc2
4907cd dgxa(j) = ii * 10d0 * size1
4908cd dgya(j) = jj * 10d0 * size2
4909cd dgxa(i) = - dgxa(j)
4910cd dgya(i) = - dgya(j)
4911 end if
4912
4913 2000 continue
4914 1000 continue
4915
4916 if (tm .eq. tlarge) then
4917 allok = .false.
4918 end if
4919
4920 end if
4921
4922c check rts cut
4923 if (allok) then
4924
4925 q4 = e(i1)
4926 q1 = px(i1)
4927 q2 = py(i1)
4928 q3 = pz(i1)
4929
4930 r4 = e(i2)
4931 r1 = px(i2)
4932 r2 = py(i2)
4933 r3 = pz(i2)
4934
4935 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4936 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
4937
4938 allok = allok .and. rts2 .gt. rscut2
4939 end if
4940
4941 if (.not. allok) then
4942 tm = tlarge
4943 t1 = tlarge
4944 t2 = tlarge
4945 else if (ha .gt. 0d0) then
4946 t1 = tc2a
4947 t2 = tc1a
4948 else
4949 t1 = tc1a
4950 t2 = tc2a
4951 end if
4952
4953 return
4954 end
4955
4956 subroutine isco9(i, j, allok, tm, t1, t2)
4957c this subroutine is used to decide whether there is a collision between
4958c particle i and j, if there is one allok=1, and tm gives the
4959c collision time, t1 the collision time for i,
4960c t2 the collision time for j
4961
4962 implicit double precision (a-h, o-z)
4963 parameter (MAXPTN=400001)
4964
4965 common /para2/ xmp, xmu, alpha, rscut2, cutof2
4966cc SAVE /para2/
4967 common /para5/ iconfg, iordsc
4968cc SAVE /para5/
4969 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4970 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4971 & xmass(MAXPTN), ityp(MAXPTN)
4972cc SAVE /prec2/
4973 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4974cc SAVE /prec4/
4975 common /aurec1/ jxa, jya, jza
4976cc SAVE /aurec1/
4977 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4978cc SAVE /aurec2/
4979 common /ilist1/
4980 & iscat, jscat, next(MAXPTN), last(MAXPTN),
4981 & ictype, icsta(MAXPTN),
4982 & nic(MAXPTN), icels(MAXPTN)
4983cc SAVE /ilist1/
4984 common /ilist3/ size1, size2, size3, v1, v2, v3, size
4985cc SAVE /ilist3/
4986 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4987cc SAVE /ilist5/
4988 SAVE
4989
4990 logical allok, allokp
4991
4992c preventing consecutive collisions
4993 allok = last(i) .ne. j .or. last(j) .ne. i
4994
4995 if (ft(i) .ge. ft(j)) then
4996 i1 = j
4997 i2 = i
4998 isign = -1
4999 else
5000 i1 = i
5001 i2 = j
5002 isign = 1
5003 end if
5004
5005 if (allok) then
5006 tm = tlarge
5007
5008 t1 = ft(i1)
5009 vx1 = vx(i1)
5010 vy1 = vy(i1)
5011 vz1 = vz(i1)
5012
5013 t2 = ft(i2)
5014
5015 dvx = vx(i2) - vx1
5016 dvy = vy(i2) - vy1
5017 dvz = vz(i2) - vz1
5018
5019 dt = t2 - t1
5020
5021 do 1000 ii = - 1, 1
5022 do 2000 jj = - 1, 1
5023
5024 allokp = .true.
5025
5026 dx = gx(i2) - gx(i1) - vx1 * dt
5027 dy = gy(i2) - gy(i1) - vy1 * dt
5028 dz = gz(i2) - gz(i1) - vz1 * dt
5029
5030 dx = dx + ii * 10d0 * size1
5031 dy = dy + jj * 10d0 * size2
5032
5033 vp = dvx * dx + dvy * dy + dvz * dz
5034
5035 allokp = allokp .and. vp .lt. 0d0
5036
5037 if (allokp) then
5038
5039 v2 = dvx * dvx + dvy * dvy + dvz * dvz
5040
5041 if (v2 .eq. 0d0) then
5042 tmp = tlarge
5043 else
5044 tmp = t2 - vp / v2
5045 end if
5046
5047c note now tm is the absolute time
5048
5049 allokp = allokp .and. tmp .gt. t1 .and.
5050 & tmp .gt. t2
5051
5052 end if
5053
5054 if (allokp) then
5055
5056 dgx = dx - dvx * t2
5057 dgy = dy - dvy * t2
5058 dgz = dz - dvz * t2
5059
5060 dm2 = - v2 * tmp ** 2 + dgx * dgx +
5061 & dgy * dgy + dgz * dgz
5062
5063 allokp = allokp .and. dm2 .lt. cutof2
5064
5065 end if
5066
5067 if (allokp .and. tmp .lt. tm) then
5068 tm = tmp
5069 jxa = isign * ii
5070 jya = isign * jj
5071 end if
5072
5073 2000 continue
5074 1000 continue
5075
5076 if (tm .eq. tlarge) then
5077 allok = .false.
5078 end if
5079 end if
5080
5081 if (allok) then
5082
5083 e1 = e(i1)
5084 px1 = px(i1)
5085 py1 = py(i1)
5086 pz1 = pz(i1)
5087 e2 = e(i2)
5088 px2 = px(i2)
5089 py2 = py(i2)
5090 pz2 = pz(i2)
5091
5092 rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
5093 & - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
5094
5095 allok = allok .and. rts2 .gt. rscut2
5096 end if
5097
5098 if (.not. allok) then
5099 tm = tlarge
5100 t1 = tlarge
5101 t2 = tlarge
5102 else
5103 t1 = tm
5104 t2 = tm
5105 end if
5106
5107 return
5108 end
5109
5110 subroutine isco10(i, j, allok, tm, t1, t2)
5111c this subroutine is used to decide whether there is a collision between
5112c particle i and j, if there is one allok=1, and tm gives the
5113c collision time, t1 the collision time for i,
5114c t2 the collision time for j
5115
5116 implicit double precision (a-h, o-z)
5117 parameter (MAXPTN=400001)
5118
5119 common /para2/ xmp, xmu, alpha, rscut2, cutof2
5120cc SAVE /para2/
5121 common /para5/ iconfg, iordsc
5122cc SAVE /para5/
5123 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5124 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5125 & xmass(MAXPTN), ityp(MAXPTN)
5126cc SAVE /prec2/
5127 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5128cc SAVE /prec4/
5129 common /aurec1/ jxa, jya, jza
5130cc SAVE /aurec1/
5131 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5132cc SAVE /aurec2/
5133 common /ilist1/
5134 & iscat, jscat, next(MAXPTN), last(MAXPTN),
5135 & ictype, icsta(MAXPTN),
5136 & nic(MAXPTN), icels(MAXPTN)
5137cc SAVE /ilist1/
5138 common /ilist3/ size1, size2, size3, v1, v2, v3, size
5139cc SAVE /ilist3/
5140 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5141cc SAVE /ilist5/
5142 SAVE
5143
5144 logical allok, allokp
5145
5146c preventing consecutive collisions
5147 allok = last(i) .ne. j .or. last(j) .ne. i
5148
5149c set up numbers for later calculations
5150
5151 tm = tlarge
5152
5153 if (allok) then
5154 do 1000 ii = - 1, 1
5155 do 2000 jj = - 1, 1
5156 do 3000 kk = -1, 1
5157 allokp = .true.
5158
5159 i1 = i
5160 i2 = j
5161
5162 p4 = ft(j) - ft(i)
5163 p1 = gx(j) - gx(i)
5164 p2 = gy(j) - gy(i)
5165 p3 = gz(j) - gz(i)
5166
5167 p1 = p1 + ii * 10d0 * size1
5168 p2 = p2 + jj * 10d0 * size2
5169 p3 = p3 + kk * 10d0 * size3
5170
5171 q4 = e(i)
5172 q1 = px(i)
5173 q2 = py(i)
5174 q3 = pz(i)
5175
5176 r4 = e(j)
5177 r1 = px(j)
5178 r2 = py(j)
5179 r3 = pz(j)
5180
5181 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
5182 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
5183 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
5184 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
5185 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
5186 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
5187
5188c make sure particle 2 formed early
5189 h = a + b
5190 if (h .gt. 0d0) then
5191 g = a
5192 a = -b
5193 b = -g
5194 g = c
5195 c = d
5196 d = g
5197 i1 = j
5198 i2 = i
5199 end if
5200
5201c check the approaching criteria
5202 if (allokp) then
5203 vp = a * d - b * ee
5204 allokp = allokp .and. vp .lt. 0d0
5205 end if
5206
5207c check the closest approach distance criteria
5208 if (allokp) then
5209 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
5210 & (ee ** 2 - c * d)
5211 allokp = allokp .and. dm2 .lt. cutof2
5212 end if
5213
5214c check the time criteria
5215 if (allokp) then
5216 tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
5217 tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
5218 tmp = 0.5d0 * (tc1 + tc2)
5219 allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
5220 end if
5221
5222 if (allokp .and. tmp .lt. tm) then
5223 tm = tmp
5224 jxa = ii
5225 jya = jj
5226 jza = kk
5227cd dgxa(j) = ii * 10d0 * size1
5228cd dgya(j) = jj * 10d0 * size2
5229cd dgxa(i) = - dgxa(j)
5230cd dgya(i) = - dgya(j)
5231 end if
5232
5233 3000 continue
5234 2000 continue
5235 1000 continue
5236
5237 if (tm .eq. tlarge) then
5238 allok = .false.
5239 end if
5240
5241 end if
5242
5243c check rts cut
5244 if (allok) then
5245
5246 q4 = e(i1)
5247 q1 = px(i1)
5248 q2 = py(i1)
5249 q3 = pz(i1)
5250
5251 r4 = e(i2)
5252 r1 = px(i2)
5253 r2 = py(i2)
5254 r3 = pz(i2)
5255
5256 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
5257 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
5258
5259 allok = allok .and. rts2 .gt. rscut2
5260 end if
5261
5262 if (.not. allok) then
5263 tm = tlarge
5264 t1 = tlarge
5265 t2 = tlarge
5266 else if (h .gt. 0d0) then
5267 t1 = tm
5268 t2 = tm
5269 else
5270 t1 = tm
5271 t2 = tm
5272 end if
5273
5274 return
5275 end
5276
5277 subroutine isco11(i, j, allok, tm, t1, t2)
5278c this subroutine is used to decide whether there is a collision between
5279c particle i and j, if there is one allok=1, and tm gives the
5280c collision time, t1 the collision time for i,
5281c t2 the collision time for j
5282
5283 implicit double precision (a-h, o-z)
5284 parameter (MAXPTN=400001)
5285
5286 common /para2/ xmp, xmu, alpha, rscut2, cutof2
5287cc SAVE /para2/
5288 common /para5/ iconfg, iordsc
5289cc SAVE /para5/
5290 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5291 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5292 & xmass(MAXPTN), ityp(MAXPTN)
5293cc SAVE /prec2/
5294 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5295cc SAVE /prec4/
5296 common /aurec1/ jxa, jya, jza
5297cc SAVE /aurec1/
5298 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5299cc SAVE /aurec2/
5300 common /ilist1/
5301 & iscat, jscat, next(MAXPTN), last(MAXPTN),
5302 & ictype, icsta(MAXPTN),
5303 & nic(MAXPTN), icels(MAXPTN)
5304cc SAVE /ilist1/
5305 common /ilist3/ size1, size2, size3, v1, v2, v3, size
5306cc SAVE /ilist3/
5307 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5308cc SAVE /ilist5/
5309 SAVE
5310
5311 logical allok, allokp
5312
5313c preventing consecutive collisions
5314 allok = last(i) .ne. j .or. last(j) .ne. i
5315
5316c set up numbers for later calculations
5317
5318 tm = tlarge
5319
5320 if (allok) then
5321 do 1000 ii = - 1, 1
5322 do 2000 jj = - 1, 1
5323 do 3000 kk = - 1, 1
5324
5325 allokp = .true.
5326
5327 i1 = i
5328 i2 = j
5329
5330 p4 = ft(j) - ft(i)
5331 p1 = gx(j) - gx(i)
5332 p2 = gy(j) - gy(i)
5333 p3 = gz(j) - gz(i)
5334
5335 p1 = p1 + ii * 10d0 * size1
5336 p2 = p2 + jj * 10d0 * size2
5337 p3 = p3 + kk * 10d0 * size3
5338
5339 q4 = e(i)
5340 q1 = px(i)
5341 q2 = py(i)
5342 q3 = pz(i)
5343
5344 r4 = e(j)
5345 r1 = px(j)
5346 r2 = py(j)
5347 r3 = pz(j)
5348
5349 a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
5350 b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
5351 c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
5352 d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
5353 ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
5354 f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
5355
5356c make sure particle 2 formed early
5357 h = a + b
5358 if (h .gt. 0d0) then
5359 g = a
5360 a = -b
5361 b = -g
5362 g = c
5363 c = d
5364 d = g
5365 i1 = j
5366 i2 = i
5367 end if
5368
5369c check the approaching criteria
5370 if (allokp) then
5371 vp = a * d - b * ee
5372 allokp = allokp .and. vp .lt. 0d0
5373 end if
5374
5375c check the closest approach distance criteria
5376 if (allokp) then
5377 dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
5378 & (ee ** 2 - c * d)
5379 allokp = allokp .and. dm2 .lt. cutof2
5380 end if
5381
5382c check the time criteria
5383 if (allokp) then
5384 tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
5385 tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
5386 if (iordsc .eq. 20) then
5387 tmp = min(tc1, tc2)
5388 else if (iordsc .eq. 21) then
5389 tmp = 0.5d0 * (tc1 + tc2)
5390 else
5391 tmp = max(tc1, tc2)
5392 end if
5393 allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
5394 end if
5395
5396 if (allokp .and. tmp .lt. tm) then
5397 tm = tmp
5398 jxa = ii
5399 jya = jj
5400 jza = kk
5401 ha = h
5402 tc1a = tc1
5403 tc2a = tc2
5404cd dgxa(j) = ii * 10d0 * size1
5405cd dgya(j) = jj * 10d0 * size2
5406cd dgxa(i) = - dgxa(j)
5407cd dgya(i) = - dgya(j)
5408 end if
5409
5410 3000 continue
5411 2000 continue
5412 1000 continue
5413
5414 if (tm .eq. tlarge) then
5415 allok = .false.
5416 end if
5417
5418 end if
5419
5420c check rts cut
5421 if (allok) then
5422
5423 q4 = e(i1)
5424 q1 = px(i1)
5425 q2 = py(i1)
5426 q3 = pz(i1)
5427
5428 r4 = e(i2)
5429 r1 = px(i2)
5430 r2 = py(i2)
5431 r3 = pz(i2)
5432
5433 rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
5434 & - (q2 + r2) ** 2 - (q3 + r3) ** 2
5435
5436 allok = allok .and. rts2 .gt. rscut2
5437 end if
5438
5439 if (.not. allok) then
5440 tm = tlarge
5441 t1 = tlarge
5442 t2 = tlarge
5443 else if (ha .gt. 0d0) then
5444 t1 = tc2a
5445 t2 = tc1a
5446 else
5447 t1 = tc1a
5448 t2 = tc2a
5449 end if
5450
5451 return
5452 end
5453
5454 subroutine isco12(i, j, allok, tm, t1, t2)
5455c this subroutine is used to decide whether there is a collision between
5456c particle i and j, if there is one allok=1, and tm gives the
5457c collision time, t1 the collision time for i,
5458c t2 the collision time for j
5459
5460 implicit double precision (a-h, o-z)
5461 parameter (MAXPTN=400001)
5462
5463 common /para2/ xmp, xmu, alpha, rscut2, cutof2
5464cc SAVE /para2/
5465 common /para5/ iconfg, iordsc
5466cc SAVE /para5/
5467 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5468 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5469 & xmass(MAXPTN), ityp(MAXPTN)
5470cc SAVE /prec2/
5471 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5472cc SAVE /prec4/
5473 common /aurec1/ jxa, jya, jza
5474cc SAVE /aurec1/
5475 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5476cc SAVE /aurec2/
5477 common /ilist1/
5478 & iscat, jscat, next(MAXPTN), last(MAXPTN),
5479 & ictype, icsta(MAXPTN),
5480 & nic(MAXPTN), icels(MAXPTN)
5481cc SAVE /ilist1/
5482 common /ilist3/ size1, size2, size3, v1, v2, v3, size
5483cc SAVE /ilist3/
5484 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5485cc SAVE /ilist5/
5486 SAVE
5487
5488 logical allok, allokp
5489
5490c preventing consecutive collisions
5491 allok = last(i) .ne. j .or. last(j) .ne. i
5492
5493 if (ft(i) .ge. ft(j)) then
5494 i1 = j
5495 i2 = i
5496 isign = -1
5497 else
5498 i1 = i
5499 i2 = j
5500 isign = 1
5501 end if
5502
5503 if (allok) then
5504 tm = tlarge
5505
5506 t1 = ft(i1)
5507 vx1 = vx(i1)
5508 vy1 = vy(i1)
5509 vz1 = vz(i1)
5510
5511 t2 = ft(i2)
5512
5513 dvx = vx(i2) - vx1
5514 dvy = vy(i2) - vy1
5515 dvz = vz(i2) - vz1
5516
5517 dt = t2 - t1
5518
5519 do 1000 ii = - 1, 1
5520 do 2000 jj = - 1, 1
5521 do 3000 kk = -1, 1
5522
5523 allokp = .true.
5524
5525 dx = gx(i2) - gx(i1) - vx1 * dt
5526 dy = gy(i2) - gy(i1) - vy1 * dt
5527 dz = gz(i2) - gz(i1) - vz1 * dt
5528
5529 dx = dx + ii * 10d0 * size1
5530 dy = dy + jj * 10d0 * size2
5531 dz = dz + kk * 10d0 * size3
5532
5533 vp = dvx * dx + dvy * dy + dvz * dz
5534
5535 allokp = allokp .and. vp .lt. 0d0
5536
5537 if (allokp) then
5538
5539 v2 = dvx * dvx + dvy * dvy + dvz * dvz
5540
5541 if (v2 .eq. 0d0) then
5542 tmp = tlarge
5543 else
5544 tmp = t2 - vp / v2
5545 end if
5546
5547c note now tm is the absolute time
5548
5549 allokp = allokp .and. tmp .gt. t1 .and.
5550 & tmp .gt. t2
5551
5552 end if
5553
5554 if (allokp) then
5555
5556 dgx = dx - dvx * t2
5557 dgy = dy - dvy * t2
5558 dgz = dz - dvz * t2
5559
5560 dm2 = - v2 * tmp ** 2 + dgx * dgx +
5561 & dgy * dgy + dgz * dgz
5562
5563 allokp = allokp .and. dm2 .lt. cutof2
5564
5565 end if
5566
5567 if (allokp .and. tmp .lt. tm) then
5568 tm = tmp
5569 jxa = isign * ii
5570 jya = isign * jj
5571 jza = isign * kk
5572 end if
5573
5574 3000 continue
5575 2000 continue
5576 1000 continue
5577
5578 if (tm .eq. tlarge) then
5579 allok = .false.
5580 end if
5581 end if
5582
5583 if (allok) then
5584
5585 e1 = e(i1)
5586 px1 = px(i1)
5587 py1 = py(i1)
5588 pz1 = pz(i1)
5589 e2 = e(i2)
5590 px2 = px(i2)
5591 py2 = py(i2)
5592 pz2 = pz(i2)
5593
5594 rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
5595 & - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
5596
5597 allok = allok .and. rts2 .gt. rscut2
5598 end if
5599
5600 if (.not. allok) then
5601 tm = tlarge
5602 t1 = tlarge
5603 t2 = tlarge
5604 else
5605 t1 = tm
5606 t2 = tm
5607 end if
5608
5609 return
5610 end
5611
5612 subroutine reor(t, tmin, j, last0)
5613c this subroutine is used to fix ct(i) when tm is greater than ct(i)
5614c next(i) is last1 or last2
5615
5616 implicit double precision (a-h, o-z)
5617 parameter (MAXPTN=400001)
5618
5619 common /para5/ iconfg, iordsc
5620cc SAVE /para5/
5621 common /ilist1/
5622 & iscat, jscat, next(MAXPTN), last(MAXPTN),
5623 & ictype, icsta(MAXPTN),
5624 & nic(MAXPTN), icels(MAXPTN)
5625cc SAVE /ilist1/
5626cd common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5627cc SAVE /ilist5/
5628 SAVE
5629
5630 icels0 = icels(j)
5631
5632 i1 = icels0 / 10000
5633 i2 = (icels0 - i1 * 10000) / 100
5634 i3 = icels0 - i1 * 10000 - i2 * 100
5635
5636 call wallc(j, i1, i2, i3, t, tmin1)
5637
5638 if (tmin .le. tmin1) then
5639 nc = last0
5640 else
5641 tmin = tmin1
5642 nc = 0
5643 end if
5644
5645 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5646 call chcell(j, i1, i2, i3, last0, t, tmin, nc)
5647 else
5648 if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
5649 call chout(j, last0, t, tmin, nc)
5650 else
5651 if (iconfg .eq. 1) then
5652 call chin1(j, i1, i2, i3, last0, t, tmin, nc)
5653 else if (iconfg .eq. 2) then
5654 call chin2(j, i1, i2, i3, last0, t, tmin, nc)
5655 else if (iconfg .eq. 4) then
5656 call chin3(j, i1, i2, i3, last0, t, tmin, nc)
5657 end if
5658 end if
5659 end if
5660
5661 call fixtim(j, t, tmin1, tmin, nc)
5662
5663 return
5664 end
5665
5666 subroutine chout(l, last0, t, tmin, nc)
5667c this subroutine is used to check the surface when the colliding
5668c particle is outside the cube
5669
5670 implicit double precision (a-h, o-z)
5671 parameter (MAXPTN=400001)
5672
5673 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5674 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5675 & xmass(MAXPTN), ityp(MAXPTN)
5676cc SAVE /prec2/
5677 SAVE
5678
5679 m1 = 11
5680 m2 = 11
5681 m3 = 11
5682 call chcell(l, m1, m2, m3, last0, t, tmin, nc)
5683
5684 do 1003 i = 1, 10
5685 do 1002 j = 1, 10
5686 do 1001 k = 1, 10
5687 if (i .eq. 1 .or. i .eq. 10 .or. j .eq. 1 .or.
5688 & j .eq. 10 .or. k .eq. 1 .or. k. eq. 10)
5689 & call chcell(l, i, j, k, last0, t, tmin, nc)
5690 1001 continue
5691 1002 continue
5692 1003 continue
5693
5694 return
5695 end
5696
5697 subroutine chin1(l, i1, i2, i3, last0, t, tmin, nc)
5698c this subroutine is used to check collisions for particle inside
5699c the cube
5700c and update the afftected particles through chcell
5701
5702 implicit double precision (a-h, o-z)
5703 SAVE
5704
5705c itest is a flag to make sure the 111111 cell is checked only once
5706 itest = 0
5707
5708 do 1003 i = i1 - 1, i1 + 1
5709 do 1002 j = i2 - 1, i2 + 1
5710 do 1001 k = i3 - 1, i3 + 1
5711 if (i .ge. 1 .and. i .le. 10
5712 & .and. j .ge. 1 .and. j .le. 10
5713 & .and. k .ge. 1 .and. k .le. 10) then
5714 call chcell(l, i, j, k, last0, t, tmin, nc)
5715 else if (itest .eq. 0) then
5716 m1 = 11
5717 m2 = 11
5718 m3 = 11
5719 call chcell(l, m1, m2, m3, last0, t, tmin, nc)
5720 itest = 1
5721 end if
5722 1001 continue
5723 1002 continue
5724 1003 continue
5725
5726 return
5727 end
5728
5729 subroutine chin2(l, i1, i2, i3, last0, t, tmin, nc)
5730c this subroutine is used to check collisions for particle inside
5731c the cube
5732c and update the afftected particles through chcell
5733
5734 implicit double precision (a-h, o-z)
5735 SAVE
5736
5737c itest is a flag to make sure the 111111 cell is checked only once
5738 itest = 0
5739
5740 do 1003 i = i1 - 1, i1 + 1
5741 do 1002 j = i2 - 1, i2 + 1
5742 do 1001 k = i3 - 1, i3 + 1
5743 ia = i
5744 ib = j
5745 ic = k
5746 if (k .ge. 1 .and. k .le. 10) then
5747 if (i .eq. 0) ia = 10
5748 if (i .eq. 11) ia = 1
5749 if (j .eq. 0) ib = 10
5750 if (j .eq. 11) ib = 1
5751 call chcell(l, ia, ib, ic, last0, t, tmin, nc)
5752 end if
5753 1001 continue
5754 1002 continue
5755 1003 continue
5756
5757 return
5758 end
5759
5760 subroutine chin3(l, i1, i2, i3, last0, t, tmin, nc)
5761c this subroutine is used to check collisions for particle inside
5762c the cube
5763c and update the afftected particles through chcell
5764
5765 implicit double precision (a-h, o-z)
5766 SAVE
5767
5768c itest is a flag to make sure the 111111 cell is checked only once
5769 itest = 0
5770
5771 do 1003 i = i1 - 1, i1 + 1
5772 do 1002 j = i2 - 1, i2 + 1
5773 do 1001 k = i3 - 1, i3 + 1
5774 if (i .eq. 0) then
5775 ia = 10
5776 else if (i .eq. 11) then
5777 ia = 1
5778 else
5779 ia = i
5780 end if
5781 if (j .eq. 0) then
5782 ib = 10
5783 else if (j .eq. 11) then
5784 ib = 1
5785 else
5786 ib = j
5787 end if
5788 if (k .eq. 0) then
5789 ic = 10
5790 else if (k .eq. 11) then
5791 ic = 1
5792 else
5793 ic = k
5794 end if
5795 call chcell(l, ia, ib, ic, last0, t, tmin, nc)
5796 1001 continue
5797 1002 continue
5798 1003 continue
5799
5800 return
5801 end
5802
5803 subroutine chcell(il, i1, i2, i3, last0, t, tmin, nc)
5804c this program is used to check through all the particles, except last0
5805c in the cell (i1,i2,i3) and see if we can get a particle collision
5806c with time less than the original input tmin ( the collision time of
5807c il with the wall
5808c last0 cas be set to 0 if we don't want to exclude last0
5809
5810 implicit double precision (a-h, o-z)
5811 parameter (MAXPTN=400001)
5812
5813 common /para5/ iconfg, iordsc
5814cc SAVE /para5/
5815 common /ilist1/
5816 & iscat, jscat, next(MAXPTN), last(MAXPTN),
5817 & ictype, icsta(MAXPTN),
5818 & nic(MAXPTN), icels(MAXPTN)
5819cc SAVE /ilist1/
5820 common /ilist2/ icell, icel(10, 10, 10)
5821cc SAVE /ilist2/
5822 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
5823cc SAVE /ilist4/
5824 SAVE
5825
5826 t=t
5827 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5828 jj = ichkpt
5829 do 1001 j = 1, jj
5830c 10/24/02 get rid of argument usage mismatch in mintm():
5831 jmintm=j
5832 if (j .ne. il .and. j .ne. last0)
5833 & call mintm(il, jmintm, tmin, nc)
5834c & call mintm(il, j, tmin, nc)
5835
5836 1001 continue
5837 return
5838 end if
5839
5840c set l
5841 if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
5842 l = icell
5843 else
5844 l = icel(i1 ,i2, i3)
5845 end if
5846
5847 if (l .eq. 0) return
5848
5849 j = nic(l)
5850
5851c if there is only one particle
5852 if (j .eq. 0) then
5853
5854c if it's not il or last0,when last is not wall
5855 if (l .eq. il .or. l .eq. last0) return
5856 call mintm(il, l, tmin, nc)
5857
5858c if there are many particles
5859 else
5860 if (l .ne. il .and. l .ne. last0)
5861 & call mintm(il, l, tmin, nc)
5862 do 10 while(j .ne. l)
5863 if (j .ne. il .and. j .ne. last0)
5864 & call mintm(il, j, tmin, nc)
5865 j = nic(j)
5866 10 continue
5867 end if
5868
5869 return
5870 end
5871
5872 subroutine mintm(i, j, tmin, nc)
5873c this subroutine is used to check whether particle j has smaller
5874c collision time with particle i than other particles
5875c or in other words, update next(i)
5876
5877c input i,j,tmin,nc
5878c output tmin,nc
5879
5880 implicit double precision (a-h, o-z)
5881 parameter (MAXPTN=400001)
5882
5883 common /para5/ iconfg, iordsc
5884cc SAVE /para5/
5885 common /aurec1/ jxa, jya, jza
5886cc SAVE /aurec1/
5887 common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5888cc SAVE /aurec2/
5889 common /ilist1/
5890 & iscat, jscat, next(MAXPTN), last(MAXPTN),
5891 & ictype, icsta(MAXPTN),
5892 & nic(MAXPTN), icels(MAXPTN)
5893cc SAVE /ilist1/
5894 common /ilist3/ size1, size2, size3, v1, v2, v3, size
5895cc SAVE /ilist3/
5896 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5897cc SAVE /ilist5/
5898 SAVE
5899
5900 logical allok
5901
5902 call isco(i, j, allok, tm, t1, t2)
5903
5904 if (allok .and. tm .lt. tmin) then
5905 tmin = tm
5906 ct(i) = t1
5907 nc = j
5908 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5909 dgxa(i) = - jxa * 10d0 * size1
5910 dgya(i) = - jya * 10d0 * size2
5911 if (iconfg .eq. 5) then
5912 dgza(i) = - jza * 10d0 * size3
5913 end if
5914 end if
5915 end if
5916
5917 return
5918 end
5919
5920******************************************************************************
5921******************************************************************************
5922
5923 subroutine zpca1
5924
5925 implicit double precision (a-h, o-z)
5926 parameter (MAXPTN=400001)
5927
5928 common /ilist1/
5929 & iscat, jscat, next(MAXPTN), last(MAXPTN),
5930 & ictype, icsta(MAXPTN),
5931 & nic(MAXPTN), icels(MAXPTN)
5932cc SAVE /ilist1/
5933 SAVE
5934
5935 if (mod(ictype,2) .eq. 0) then
5936 call zpca1a(iscat)
5937 call zpca1a(jscat)
5938clin-5/2009 ctest off v2 for parton:
5939c call flowp(1)
5940 end if
5941
5942 return
5943 end
5944
5945 subroutine zpca1a(i)
5946
5947 implicit double precision (a-h, o-z)
5948 parameter (MAXPTN=400001)
5949
5950 common /para2/ xmp, xmu, alpha, rscut2, cutof2
5951cc SAVE /para2/
5952 common /para5/ iconfg, iordsc
5953cc SAVE /para5/
5954 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5955 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5956 & xmass(MAXPTN), ityp(MAXPTN)
5957cc SAVE /prec2/
5958 common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
5959 & pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
5960 & xmasss(MAXPTN), ityps(MAXPTN)
5961cc SAVE /prec3/
5962 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
5963cc SAVE /prec5/
5964 common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
5965cc SAVE /prec6/
5966 common /ana1/ ts(12)
5967cc SAVE /ana1/
5968 SAVE
5969
5970 if (iconfg .eq. 1) then
5971 t1 = fts(i)
5972 t2 = ft(i)
5973 ipic = 11
5974 else if (iconfg .eq. 2 .or.
5975 & iconfg .eq. 3) then
5976cd t1 = fts(i)
5977cd t2 = ft(i)
5978 t1 = taus(i)
5979 t2 = tau(i)
5980 ipic = 12
5981 else if (iconfg .eq. 4 .or.
5982 & iconfg .eq. 5) then
5983 t1 = fts(i)
5984 t2 = ft(i)
5985 ipic = 12
5986 end if
5987
5988 if (iconfg .le. 3) then
5989 do 1002 ian = 1, ipic
5990 if (t1 .le. ts(ian) .and.
5991 & t2 .gt. ts(ian)) then
5992 rapi = raps(i)
5993c 7/20/01:
5994c et = sqrt(pxs(i) ** 2 + pys(i) ** 2 + xmp ** 2)
5995 et = dsqrt(pxs(i) ** 2 + pys(i) ** 2 + xmp ** 2)
5996 call zpca1b(rapi, et, ian)
5997 end if
5998 1002 continue
5999 else
6000 do 1003 ian = 1, ipic
6001 if (t1 .le. ts(ian) .and.
6002 & t2 .gt. ts(ian)) then
6003 p0 = es(i)
6004 p1 = pxs(i)
6005 p2 = pys(i)
6006 p3 = pzs(i)
6007 call zpca1c(p0, p1, p2, p3, ian)
6008 end if
6009 1003 continue
6010 end if
6011
6012 return
6013 end
6014
6015 subroutine zpca1b(rapi, et, ian)
6016
6017 implicit double precision (a-h, o-z)
6018 parameter (MAXPTN=400001)
6019
6020 common /para6/ centy
6021cc SAVE /para6/
6022 common /ilist6/ t, iopern, icolln
6023cc SAVE /ilist6/
6024 common /ana2/
6025 & det(12), dn(12), detdy(12), detdn(12), dndy(12),
6026 & det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6027 & det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6028cc SAVE /ana2/
6029 SAVE
6030
6031 if (rapi .gt. centy - 0.5d0 .and.
6032 & rapi .lt. centy + 0.5d0) then
6033 det2(ian) = det2(ian) + et
6034 dn2(ian) = dn2(ian) + 1d0
6035cdtrans
6036 if (ian .eq. 10) then
6037cd write (10, *) t, det2(ian)
6038 end if
6039 if (ian .eq. 11) then
6040cd write (11, *) t, det2(ian)
6041 end if
6042 if (ian .eq. 12) then
6043cd write (12, *) t, det2(ian)
6044 end if
6045cdtransend
6046 if (rapi .gt. centy - 0.25d0 .and.
6047 & rapi .lt. centy + 0.25d0) then
6048 det1(ian) = det1(ian) + et
6049 dn1(ian) = dn1(ian) + 1d0
6050 if (rapi .gt. centy - 0.1d0 .and.
6051 & rapi .lt. centy + 0.1d0) then
6052 det(ian) = det(ian) + et
6053 dn(ian) = dn(ian) + 1d0
6054 end if
6055 end if
6056 end if
6057
6058 return
6059 end
6060
6061 subroutine zpca1c(p0, p1, p2, p3, ian)
6062
6063 implicit double precision (a-h, o-z)
6064
6065 common /ana3/ em(4, 4, 12)
6066cc SAVE /ana3/
6067
6068 dimension en(4)
6069 SAVE
6070
6071 en(1) = p0
6072 en(2) = p1
6073 en(3) = p2
6074 en(4) = p3
6075
6076 do 1002 i = 1, 4
6077 do 1001 j = 1, 4
6078 em(i, j, ian) = em(i, j, ian) + en(i) * en(j) / p0
6079 1001 continue
6080 1002 continue
6081
6082 return
6083 end
6084
6085******************************************************************************
6086
6087 subroutine zpca2
6088
6089 implicit double precision (a-h, o-z)
6090 parameter (MAXPTN=400001)
6091
6092 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6093cc SAVE /para3/
6094 common /para5/ iconfg, iordsc
6095cc SAVE /para5/
6096 common /para7/ ioscar,nsmbbbar,nsmmeson
6097cc SAVE /para7/
6098 common /ilist6/ t, iopern, icolln
6099cc SAVE /ilist6/
6100 common /rndm1/ number
6101cc SAVE /rndm1/
6102 common /rndm2/ iff
6103cc SAVE /rndm2/
6104 common /rndm3/ iseedp
6105cc SAVE /rndm3/
6106 COMMON /AREVT/ IAEVT, IARUN, MISS
6107cc SAVE /AREVT/
6108 SAVE
6109
6110 if (iconfg .le. 3) then
6111 call zpca2a
6112 else
6113 call zpca2b
6114 end if
6115
6116 if (ioscar .eq. 1) then
6117 call zpca2c
6118 end if
6119
6120cbzdbg2/17/99
6121c write (25, *) 'Event', nsevt - 1 + ievt,
6122c & ', run', isbrun,
6123c WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN,
6124c & ',\n\t number of operations = ', iopern,
6125c & ',\n\t number of collisions between particles = ',
6126c & icolln,
6127c & ',\n\t freezeout time=', t,
6128c & ',\n\t ending at the ', number, 'th random number',
6129c & ',\n\t ending collision iff=', iff
6130cms WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN
6131cms WRITE (25, *) ' number of operations = ', iopern
6132cms WRITE (25, *) ' number of collisions between particles = ',
6133cms & icolln
6134cms WRITE (25, *) ' freezeout time=', t
6135cms WRITE (25, *) ' ending at the ', number, 'th random number'
6136cms WRITE (25, *) ' ending collision iff=', iff
6137
6138 return
6139 end
6140
6141 subroutine zpca2a
6142
6143 implicit double precision (a-h, o-z)
6144 parameter (MAXPTN=400001)
6145
6146 common /para1/ mul
6147cc SAVE /para1/
6148 common /para2/ xmp, xmu, alpha, rscut2, cutof2
6149cc SAVE /para2/
6150 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6151cc SAVE /para3/
6152 common /para5/ iconfg, iordsc
6153cc SAVE /para5/
6154 common /para6/ centy
6155cc SAVE /para6/
6156 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6157 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6158 & xmass(MAXPTN), ityp(MAXPTN)
6159cc SAVE /prec2/
6160 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
6161cc SAVE /prec5/
6162 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
6163cc SAVE /ilist4/
6164 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
6165cc SAVE /ilist5/
6166 common /ilist6/ t, iopern, icolln
6167cc SAVE /ilist6/
6168 common /rndm1/ number
6169cc SAVE /rndm1/
6170 common /rndm2/ iff
6171cc SAVE /rndm2/
6172 common /rndm3/ iseedp
6173cc SAVE /rndm3/
6174 common /ana1/ ts(12)
6175cc SAVE /ana1/
6176 common /ana2/
6177 & det(12), dn(12), detdy(12), detdn(12), dndy(12),
6178 & det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6179 & det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6180cc SAVE /ana2/
6181 common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
6182cc SAVE /ana4/
6183 SAVE
6184
6185 do 1004 i = 1, ichkpt
6186 rapi = rap(i)
6187c 7/20/01:
6188c et = sqrt(px(i) ** 2 + py(i) ** 2 + xmp ** 2)
6189 et = dsqrt(px(i) ** 2 + py(i) ** 2 + xmp ** 2)
6190
6191 do 1001 j = 1, 24
6192 if (rapi .gt. j + centy - 13d0
6193 & .and. rapi .lt. j + centy - 12d0) then
6194 fdetdy(j) = fdetdy(j) + et
6195 fdndy(j) = fdndy(j) + 1d0
6196 end if
6197 1001 continue
6198
6199 do 1002 j = 1, 12
6200 if (et .gt. 0.5d0 * (j - 1) .and.
6201 & et .lt. 0.5d0 * j ) then
6202 fdndpt(j) = fdndpt(j) + 1d0
6203 end if
6204 1002 continue
6205
6206 if (iconfg .eq. 1) then
6207 t1 = ft(i)
6208 t2 = tlarge
6209 ipic = 11
6210 else
6211 t1 = tau(i)
6212 t2 = tlarge
6213 ipic = 12
6214 end if
6215
6216 do 1003 ian = 1, ipic
6217 if (t1 .le. ts(ian) .and.
6218 & t2 .gt. ts(ian)) then
6219 call zpca1b(rapi, et, ian)
6220 end if
6221 1003 continue
6222
6223 if (iconfg .eq. 1) then
6224 call zpca1b(rapi, et, 12)
6225 end if
6226 1004 continue
6227
6228 do 1005 ian = 1, 12
6229c if (dn(ian) .eq. 0d0 .or. dn1(ian) .eq. 0d0 .or.
6230c & dn2(ian) .eq. 0d0) then
6231c print *, 'event=', ievt
6232c print *, 'dn(', ian, ')=', dn(ian), 'dn1(', ian,
6233c & ')=', dn1(ian), 'dn2(', ian, ')=', dn2(ian)
6234c end if
6235 detdy(ian) = detdy(ian) + det(ian)
6236 if (dn(ian) .ne. 0) then
6237 detdn(ian) = detdn(ian) + det(ian) / dn(ian)
6238 end if
6239 dndy(ian) = dndy(ian) + dn(ian)
6240 detdy1(ian) = detdy1(ian) + det1(ian)
6241 if (dn1(ian) .ne. 0) then
6242 detdn1(ian) = detdn1(ian) + det1(ian) / dn1(ian)
6243 end if
6244 dndy1(ian) = dndy1(ian) + dn1(ian)
6245 detdy2(ian) = detdy2(ian) + det2(ian)
6246 if (dn2(ian) .ne. 0) then
6247 detdn2(ian) = detdn2(ian) + det2(ian) / dn2(ian)
6248 end if
6249 dndy2(ian) = dndy2(ian) + dn2(ian)
6250 1005 continue
6251
6252 return
6253 end
6254
6255 subroutine zpca2b
6256
6257 implicit double precision (a-h, o-z)
6258 parameter (MAXPTN=400001)
6259
6260 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6261 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6262 & xmass(MAXPTN), ityp(MAXPTN)
6263cc SAVE /prec2/
6264 common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
6265cc SAVE /ilist4/
6266 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
6267cc SAVE /ilist5/
6268 common /ana1/ ts(12)
6269cc SAVE /ana1/
6270 SAVE
6271
6272 do 1002 i = 1, ichkpt
6273 t1 = ft(i)
6274 t2 = tlarge
6275 ipic = 12
6276
6277 do 1001 ian = 1, ipic
6278 if (t1 .le. ts(ian) .and.
6279 & t2 .gt. ts(ian)) then
6280 p0 = e(i)
6281 p1 = px(i)
6282 p2 = py(i)
6283 p3 = pz(i)
6284 call zpca1c(p0, p1, p2, p3, ian)
6285 end if
6286 1001 continue
6287 1002 continue
6288
6289 return
6290 end
6291
6292 subroutine zpca2c
6293
6294 implicit double precision (a-h, o-z)
6295
6296 character*8 code, versn
6297 character*4 reffra
6298 integer aproj, zproj, atarg, ztarg, event
6299 parameter (MAXPTN=400001)
6300
6301 common /para1/ mul
6302cc SAVE /para1/
6303 common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6304 & px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6305 & xmass(MAXPTN), ityp(MAXPTN)
6306cc SAVE /prec2/
6307 SAVE
6308 data nff/0/
6309
6310c file header
6311 if (nff .eq. 0) then
6312cms write (26, 101) 'OSCAR1997A'
6313cms write (26, 101) 'final_id_p_x'
6314 code = 'ZPC'
6315 versn = '1.0.1'
6316 aproj = -1
6317 zproj = -1
6318 atarg = -1
6319 ztarg = -1
6320 reffra = 'cm'
6321 ebeam = 0d0
6322 ntestp = 1
6323cms write (26, 102) code, versn, aproj, zproj, atarg, ztarg,
6324cms & reffra, ebeam, ntestp
6325 nff = 1
6326 event = 1
6327 bimp = 0d0
6328 phi = 0d0
6329 end if
6330
6331c comment
6332
6333c event header
6334cms write (26, 103) event, mul, bimp, phi
6335
6336c particles
6337 do 99 i = 1, mul
6338cms write (26, 104) i, ityp(i),
6339cms & px(i), py(i), pz(i), e(i), xmass(i),
6340cms & gx(i), gy(i), gz(i), ft(i)
6341 99 continue
6342
6343 event = event + 1
6344
6345cyy 101 format (a12)
6346cyy 102 format (2(a8, 2x), '(',i3, ',',i6, ')+(',i3, ',', i6, ')',
6347cyy & 2x, a4, 2x, e10.4, 2x, i8)
6348cyy 103 format (i10, 2x, i10, 2x, f8.3, 2x, f8.3)
6349cyy 104 format (i10, 2x, i10, 2x, 9(e12.6, 2x))
6350
6351 return
6352 end
6353
6354******************************************************************************
6355
6356 subroutine zpcou
6357
6358 implicit double precision (a-h, o-z)
6359
6360 common /para5/ iconfg, iordsc
6361cc SAVE /para5/
6362 SAVE
6363
6364 if (iconfg .le. 3) then
6365 call zpcou1
6366 else
6367 call zpcou2
6368 end if
6369
6370 return
6371 end
6372
6373 subroutine zpcou1
6374
6375 implicit double precision (a-h, o-z)
6376
6377 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6378cc SAVE /para3/
6379 common /ana1/ ts(12)
6380cc SAVE /ana1/
6381 common /ana2/
6382 & det(12), dn(12), detdy(12), detdn(12), dndy(12),
6383 & det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6384 & det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6385cc SAVE /ana2/
6386 common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
6387cc SAVE /ana4/
6388 SAVE
6389c
6390 dpt = 0.5d0
6391 dy2 = 1d0
6392 dy1 = 0.5d0
6393 dy = 0.2d0
6394 ntotal = nevnt * nsbrun
6395c
6396 return
6397 end
6398
6399 subroutine zpcou2
6400
6401 implicit double precision (a-h, o-z)
6402
6403 common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6404cc SAVE /para3/
6405 common /ilist3/ size1, size2, size3, v1, v2, v3, size
6406cc SAVE /ilist3/
6407 common /ana1/ ts(12)
6408cc SAVE /ana1/
6409 common /ana3/ em(4, 4, 12)
6410cc SAVE /ana3/
6411 SAVE
6412c
6413cms open (28, file = 'ana4/em.dat', status = 'unknown')
6414 vol = 1000.d0 * size1 * size2 * size3
6415 ntotal = nevnt * nsbrun
6416
6417 do 1002 ian = 1, 12
6418cms write (28, *) '*** for time ', ts(ian), 'fm(s)'
6419 do 1001 i = 1, 4
6420cms write (28, *) em(i, 1, ian) / vol / ntotal,
6421cms & em(i, 2, ian) / vol / ntotal,
6422cms & em(i, 3, ian) / vol / ntotal,
6423cms & em(i, 4, ian) / vol / ntotal
6424 1001 continue
6425 1002 continue
6426
6427 return
6428 end
6429
6430******************************************************************************
6431
e2670fbe 6432 subroutine lorenza(energy, px, py, pz, bex, bey, bez)
0119ef9a 6433
6434c add in a cut for beta2 to prevent gam to be nan (infinity)
6435
6436 implicit double precision (a-h, o-z)
6437
3006c44b 6438 common /lora/ enenew, pxnew, pynew, pznew
6439cc SAVE /lora/
0119ef9a 6440 SAVE
6441
6442 beta2 = bex ** 2 + bey ** 2 + bez ** 2
6443 if (beta2 .eq. 0d0) then
6444 enenew = energy
6445 pxnew = px
6446 pynew = py
6447 pznew = pz
6448 else
6449 if (beta2 .gt. 0.999999999999999d0) then
6450 beta2 = 0.999999999999999d0
6451 print *,'beta2=0.999999999999999'
6452 end if
6453clin-7/20/01:
6454c gam = 1.d0 / sqrt(1.d0 - beta2)
6455 gam = 1.d0 / dsqrt(1.d0 - beta2)
6456 enenew = gam * (energy - bex * px - bey * py - bez * pz)
6457 pxnew = - gam * bex * energy + (1.d0
6458 & + (gam - 1.d0) * bex ** 2 / beta2) * px
6459 & + (gam - 1.d0) * bex * bey/beta2 * py
6460 & + (gam - 1.d0) * bex * bez/beta2 * pz
6461 pynew = - gam * bey * energy
6462 & + (gam - 1.d0) * bex * bey / beta2 * px
6463 & + (1.d0 + (gam - 1.d0) * bey ** 2 / beta2) * py
6464 & + (gam - 1.d0) * bey * bez / beta2 * pz
6465 pznew = - gam * bez * energy
6466 & + (gam - 1.d0) * bex * bez / beta2 * px
6467 & + (gam - 1.d0) * bey * bez / beta2 * py
6468 & + (1.d0 + (gam - 1.d0) * bez ** 2 / beta2) * pz
6469 endif
6470
6471 return
6472 end
6473
6474 subroutine index1(n, m, arrin, indx)
6475c indexes the first m elements of ARRIN of length n, i.e., outputs INDX
6476c such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m
6477
6478 implicit double precision (a-h, o-z)
6479
6480 dimension arrin(n), indx(n)
6481 SAVE
6482 do 1001 j = 1, m
6483 indx(j) = j
6484 1001 continue
6485 l = m / 2 + 1
6486 ir = m
6487 10 continue
6488 if (l .gt. 1) then
6489 l = l - 1
6490 indxt = indx(l)
6491 q = arrin(indxt)
6492 else
6493 indxt = indx(ir)
6494 q = arrin(indxt)
6495 indx(ir) = indx(1)
6496 ir = ir - 1
6497 if (ir .eq. 1) then
6498 indx(1) = indxt
6499 return
6500 end if
6501 end if
6502 i = l
6503 j = l + l
6504 20 if (j .le. ir) then
6505 if (j .lt. ir) then
6506 if (arrin(indx(j)) .lt. arrin(indx(j + 1))) j = j + 1
6507 end if
6508 if (q .lt. arrin(indx(j))) then
6509 indx(i) = indx(j)
6510 i = j
6511 j = j + j
6512 else
6513 j = ir + 1
6514 end if
6515 goto 20
6516 end if
6517 indx(i) = indxt
6518 goto 10
6519
6520 end
6521
6522
6523 double precision function ftime1(iseed)
6524
6525c this program is used to generate formation time
6526c the calling program needs a common /par1/
6527c and declare external ftime1
6528
6529clin-8/19/02
6530 implicit double precision (a-h, o-z)
6531
6532cc external ran1
6533
6534 parameter (hbarc = 0.197327054d0)
6535
6536 common /par1/ formt
6537cc SAVE /par1/
6538 SAVE
6539
6540 aa = hbarc / formt
6541
6542clin7/20/01:
6543c ftime1 = aa * sqrt(1d0 / ran1(iseed) - 1d0)
6544 ftime1 = aa * dsqrt(1d0 / ran1(iseed) - 1d0)
6545 return
6546 end
6547
6548
6549 subroutine cropro(vx1, vy1, vz1, vx2, vy2, vz2)
6550
6551c this subroutine is used to calculate the cross product of
6552c (vx1,vy1,vz1) and (vx2,vy2,vz2) and get the result (vx3,vy3,vz3)
6553c and put the vector into common /cprod/
6554
6555 implicit double precision (a-h, o-z)
6556
6557 common/cprod/ vx3, vy3, vz3
6558cc SAVE /cprod/
6559 SAVE
6560
6561 vx3 = vy1 * vz2 - vz1 * vy2
6562 vy3 = vz1 * vx2 - vx1 * vz2
6563 vz3 = vx1 * vy2 - vy1 * vx2
6564
6565 return
6566 end
6567
6568 subroutine xnormv(vx, vy, vz)
6569
6570c this subroutine is used to get a normalized vector
6571
6572 implicit double precision (a-h, o-z)
6573 SAVE
6574
6575clin-7/20/01:
6576c vv = sqrt(vx ** 2 + vy ** 2 + vz ** 2)
6577 vv = dsqrt(vx ** 2 + vy ** 2 + vz ** 2)
6578 vx = vx / vv
6579 vy = vy / vv
6580 vz = vz / vv
6581
6582 return
6583 end
6584
6585cbz1/29/99
6586c subroutine rotate(xn1, xn2, xn3, theta, v1, v2, v3)
6587 subroutine zprota(xn1, xn2, xn3, theta, v1, v2, v3)
6588cbz1/29/99end
6589
6590c this subroutine is used to rotate the vector (v1,v2,v3) by an angle theta
6591c around the unit vector (xn1, xn2, xn3)
6592
6593 implicit double precision (a-h, o-z)
6594 SAVE
6595
6596 vx = v1
6597 vy = v2
6598 vz = v3
6599 c = cos(theta)
6600 omc = 1d0 - c
6601 s = sin(theta)
6602 a11 = xn1 ** 2 * omc + c
6603 a12 = xn1 * xn2 * omc - s * xn3
6604 a13 = xn1 * xn3 * omc + s * xn2
6605 a21 = xn1 * xn2 * omc + s * xn3
6606 a22 = xn2 **2 * omc + c
6607 a23 = xn2 * xn3 * omc - s * xn1
6608 a31 = xn1 * xn3 * omc - s * xn2
6609 a32 = xn3 * xn2 * omc + s * xn1
6610 a33 = xn3 ** 2 * omc + c
6611 v1 = vx * a11 + vy * a12 + vz * a13
6612 v2 = vx * a21 + vy * a22 + vz * a23
6613 v3 = vx * a31 + vy * a32 + vz * a33
6614
6615 return
6616 end
6617
6618c double precision function ran1(idum)
6619
6620c* return a uniform random deviate between 0.0 and 1.0. set idum to
6621c* any negative value to initialize or reinitialize the sequence.
6622
6623c implicit double precision (a-h, o-z)
6624
6625c dimension r(97)
6626
6627c common /rndm1/ number
6628cc SAVE /rndm1/
6629c parameter (m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 1d0 / m1)
6630c parameter (m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 1d0 / m2)
6631c parameter (m3 = 243000, ia3 = 4561, ic3 = 51349)
6632clin-6/23/00 save ix1-3:
6633clin-10/30/02 r unsaved, causing wrong values for ran1 when compiled with f77:
6634cc SAVE ix1,ix2,ix3,r
6635c SAVE
6636c data iff/0/
6637
6638c if (idum .lt. 0 .or. iff .eq. 0) then
6639c iff = 1
6640c ix1 = mod(ic1 - idum, m1)
6641c ix1 = mod(ia1 * ix1 + ic1, m1)
6642c ix2 = mod(ix1, m2)
6643c ix1 = mod(ia1 * ix1 + ic1, m1)
6644c ix3 = mod(ix1, m3)
6645c do 11 j = 1, 97
6646c ix1 = mod(ia1 * ix1 + ic1, m1)
6647c ix2 = mod(ia2 * ix2 + ic2, m2)
6648c r(j) = (dble(ix1) + dble(ix2) * rm2) * rm1
6649c 11 continue
6650c idum = 1
6651c end if
6652c ix1 = mod(ia1 * ix1 + ic1, m1)
6653c ix2 = mod(ia2 * ix2 + ic2, m2)
6654c ix3 = mod(ia3 * ix3 + ic3, m3)
6655clin-7/01/02 j = 1 + (97 * i x 3) / m3
6656c j=1+(97*ix3)/m3
6657clin-4/2008:
6658c if (j .gt. 97 .or. j .lt. 1) pause
6659c if (j .gt. 97 .or. j .lt. 1) print *, 'In zpc ran1, j<1 or j>97',j
6660c ran1 = r(j)
6661c r(j) = (dble(ix1) + dble(ix2) * rm2) * rm1
6662
6663clin-6/23/00 check random number generator:
6664c number = number + 1
6665c if(number.le.100000) write(99,*) 'number, ran1=', number,ran1
6666
6667c return
6668c end