]>
Commit | Line | Data |
---|---|---|
9ef1c2d9 | 1 | c----------------------------------------------------------------------- |
2 | c activate making inicon table by the follwing commands in optns file: | |
3 | c----------------------------------------------------------------------- | |
4 | c set nevent 1 | |
5 | c set bminim 3 (for example) | |
6 | c set bmaxim 5 (for example) | |
7 | c set ninicon 1000 (for example) | |
8 | c set icocore 1 !or 2 | |
9 | c make icotable | |
10 | c switch splitting on | |
11 | c switch fusion on | |
12 | c switch clusterdecay off | |
13 | c----------------------------------------------------------------------- | |
14 | c inicon table means storage of IcoE, IcoV, IcoF | |
15 | c IcoE ... energy density in the comoving frame | |
16 | c IcoV ... flow velocity in hyperbolic coordinates | |
17 | c v=(vx,vy,veta=vz/tau) | |
18 | c with vz being the velocity in the frame moving with rapidity eta | |
19 | c IcoV ... net flavor density in the comoving frame | |
20 | c the indices of these fields ix,iy,iz also refer to hyperbolic coordinates | |
21 | c (x,y,eta) at given tau | |
22 | c the corresponding grid is defined in epos.incico | |
23 | c----------------------------------------------------------------------- | |
24 | ||
25 | c----------------------------------------------------------------------- | |
26 | subroutine IniCon(nin) | |
27 | c----------------------------------------------------------------------- | |
28 | include 'epos.inc' | |
29 | include 'epos.incico' | |
30 | double precision seedf | |
31 | character*80 fn | |
32 | logical table | |
33 | ||
34 | ii=index(fnhi(1:nfnhi),".histo")-1 | |
35 | fn=fnhi(1:ii)//".ico" | |
36 | inquire(file=fn(1:ii+4),exist=table) | |
37 | ||
38 | if(icotabr.eq.1)then !read from file | |
39 | ||
40 | if(.not.table)then | |
41 | write(ifmt,'(//10x,3a//)')'STOP: file ' | |
42 | * ,fn(1:ii+4),' not found !!!' | |
43 | stop | |
44 | endif | |
45 | write(ifmt,'(3a)')'read from ',fn(1:ii+4),' ...' | |
46 | open(97,file=fn,status='old') | |
47 | read(97,*) laprojx,maprojx,latargx,matargx | |
48 | read(97,*) engyx | |
49 | read(97,*) bminimx,bmaximx | |
50 | read(97,*) tauicox | |
51 | read(97,*) iabs_ninicon | |
52 | read(97,*) nxicox,nyicox,nzicox | |
53 | read(97,*)xminicox,xmaxicox,yminicox,ymaxicox,zminicox,zmaxicox | |
54 | if(laprojx.ne.laproj.or.maprojx.ne.maproj | |
55 | * .or.latargx.ne.latarg.or.matargx.ne.matarg)stop'(1)' | |
56 | if(engyx.ne.engy) stop'variables dont match (2) ' | |
57 | if(bminimx.ne.bminim.or.bmaximx.ne.bmaxim) | |
58 | * stop'variables dont match (3) ' | |
59 | if(tauicox.ne.tauico)stop'variables dont match (4) ' | |
60 | if(nxicox.ne.nxico.or.nyicox.ne.nyico | |
61 | * .or.nzicox.ne.nzico)stop'(5)' | |
62 | if( xminicox.ne.xminico.or.xmaxicox.ne.xmaxico | |
63 | * .or.yminicox.ne.yminico.or.ymaxicox.ne.ymaxico | |
64 | * .or.zminicox.ne.zminico.or.zmaxicox.ne.zmaxico)stop'(6)' | |
65 | read(97,*) IcoE,IcoV,IcoF | |
66 | close(97) | |
67 | ||
68 | elseif(icotabr.eq.0)then ! calculate | |
69 | ||
70 | !if(icotabm.eq.1.and.table)then | |
71 | !write(ifmt,'(//10x,3a//)')'STOP: file ' | |
72 | !* ,fn(1:ii+4),' already exists !!!' | |
73 | !stop | |
74 | !endif | |
75 | ||
76 | if(nin.eq.1)then | |
77 | write(ifmt,'(2a,i7,a)')' calculate inicon table ', | |
78 | & 'based on',iabs(ninicon),' ico-events ...' | |
79 | call IcoStr(1) | |
80 | endif | |
81 | ||
82 | call ranfgt(seedf) | |
83 | if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0) | |
84 | & write(ifmt,'(a,i7,a,f6.2,a,d27.16)')' ico-event ',nin | |
85 | & ,' bimevt:',bimevt,' seedf:',seedf | |
86 | !if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0) | |
87 | !& print*,'+++++ time: ',timefin-timeini | |
88 | seedc=seedf | |
89 | call IcoStr(2) | |
90 | ||
91 | if(nin.eq.iabs(ninicon))then | |
92 | write(ifmt,*)'normalize and diagonalize engy-mom tensor' | |
93 | call IcoStr(3) | |
94 | ||
95 | if(icotabm.eq.1)then ! make table | |
96 | write(ifmt,'(2a)')' write inicon table into ',fn(1:ii+4) | |
97 | open(97,file=fn,status='unknown') | |
98 | write(97,*) laproj,maproj,latarg,matarg | |
99 | write(97,*) engy | |
100 | write(97,*) bminim,bmaxim | |
101 | write(97,*) tauico | |
102 | write(97,*) iabs(ninicon) | |
103 | write(97,*) nxico,nyico,nzico | |
104 | write(97,*) xminico,xmaxico,yminico,ymaxico,zminico,zmaxico | |
105 | write(97,*) IcoE,IcoV,IcoF | |
106 | close(97) | |
107 | endif | |
108 | endif | |
109 | ||
110 | else | |
111 | ||
112 | stop'IniCon: wrong choice' | |
113 | ||
114 | endif | |
115 | ||
116 | end | |
117 | ||
118 | c----------------------------------------------------------------------- | |
119 | subroutine IcoStr(iflag) | |
120 | c----------------------------------------------------------------------- | |
121 | c energy momentum tensor and flavor current | |
122 | c and corresponding contractions from string method | |
123 | c | |
124 | c iflag = 1 initialization | |
125 | c = 2 sum up density | |
126 | c = 3 normalization | |
127 | c | |
128 | c output: common blocks /Ico3/ /Ico4/ /Ico5/ | |
129 | c----------------------------------------------------------------------- | |
130 | include 'epos.inc' | |
131 | include 'epos.incico' | |
132 | double precision IcoT(4,4,nxico,nyico,nzico) | |
133 | * ,IcoC(3,4,nxico,nyico,nzico) | |
134 | common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl) | |
135 | *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3) | |
136 | common /Ico10/ntot,ish0,irs | |
137 | double precision p(4),v(4),u(4),eps,T(4,4),gamma,rap,rapx | |
138 | integer ic(2),jc(nflav,2) | |
139 | logical ok | |
140 | goto (1,10,100),iflag | |
141 | return | |
142 | ||
143 | c...................................................................... | |
144 | c initialization | |
145 | c...................................................................... | |
146 | ||
147 | 1 do ix=1,nxico | |
148 | do iy=1,nyico | |
149 | do iz=1,nzico | |
150 | IcoE(ix,iy,iz)=0. | |
151 | do i=1,3 | |
152 | IcoV(i,ix,iy,iz)=0. | |
153 | enddo | |
154 | do i=1,3 | |
155 | IcoF(i,ix,iy,iz)=0. | |
156 | enddo | |
157 | do i=1,4 | |
158 | do j=1,4 | |
159 | IcoT(i,j,ix,iy,iz)=0. | |
160 | enddo | |
161 | enddo | |
162 | do i=1,3 | |
163 | do j=1,4 | |
164 | IcoC(i,j,ix,iy,iz)=0. | |
165 | enddo | |
166 | enddo | |
167 | enddo | |
168 | enddo | |
169 | enddo | |
170 | ntot=0 | |
171 | return | |
172 | ||
173 | c...................................................................... | |
174 | c fill arrays | |
175 | c...................................................................... | |
176 | ||
177 | 10 ntot=ntot+1 | |
178 | dxico=(xmaxico-xminico)/nxico | |
179 | dyico=(ymaxico-yminico)/nyico | |
180 | dzico=(zmaxico-zminico )/nzico | |
181 | ||
182 | do j=1,nptl | |
183 | ok=.true. | |
184 | if(icocore.eq.2 | |
185 | . .and.(ityptl(j).lt.20.or.ityptl(j).ge.40))ok=.false. | |
186 | if(ok | |
187 | . .and.dezptl(j).lt.1e3.and.j.le.maxfra.and.istptl(j).eq.2)then | |
188 | aa=cos(phievt) | |
189 | bb=sin(phievt) | |
190 | x=xptl(j)*aa+yptl(j)*bb | |
191 | cc=-sin(phievt) | |
192 | dd=cos(phievt) | |
193 | y=xptl(j)*cc+yptl(j)*dd | |
194 | z=dezptl(j) | |
195 | ix=1+(x-xminico)/dxico | |
196 | iy=1+(y-yminico)/dyico | |
197 | iz=1+(z-zminico)/dzico | |
198 | if(ix.ge.1.and.ix.le.nxico.and. | |
199 | . iy.ge.1.and.iy.le.nyico.and. | |
200 | . iz.ge.1.and.iz.le.nzico)then | |
201 | !~~~~~~~~~~~~ | |
202 | ! T^i^k = \int d3p p^i p^k /E f(\vec x,\vec p) | |
203 | ! N^k = \int d3p p^k /E f_net(\vec x,\vec p) | |
204 | !~~~~~~~~~~~~ | |
205 | rapx=zminico+(iz-0.5)*dzico | |
206 | amt2=pptl(5,j)**2+pptl(1,j)**2+pptl(2,j)**2 | |
207 | pp=pptl(4,j)+abs(pptl(3,j)) | |
208 | if(amt2.gt.0..and.pp.gt.0.)then | |
209 | amt=sqrt(amt2) | |
210 | rap=sign(1.,pptl(3,j))*alog(pp/amt) | |
211 | p(1)=pptl(1,j) | |
212 | p(2)=pptl(2,j) | |
213 | p(3)=amt*sinh(rap-rapx) | |
214 | p(4)=amt*cosh(rap-rapx) | |
215 | do i=1,4 | |
216 | do k=1,4 | |
217 | IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)+p(i)*p(k)/p(4) | |
218 | enddo | |
219 | enddo | |
220 | id=idptl(j) | |
221 | ida=iabs(id/10) | |
222 | ids=id/iabs(id) | |
223 | if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10 | |
224 | if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids | |
225 | if(ida.eq.213)id=1230*ids | |
226 | ic(1)=idtrai(1,id,1) | |
227 | ic(2)=idtrai(2,id,1) | |
228 | call iddeco(ic,jc) | |
229 | do i=1,3 | |
230 | fi=jc(i,1)-jc(i,2) | |
231 | do k=1,4 | |
232 | IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)+p(k)/p(4)*fi | |
233 | enddo | |
234 | enddo | |
235 | endif | |
236 | !~~~~~~~~~~~~ | |
237 | endif | |
238 | endif | |
239 | enddo | |
240 | ||
241 | return | |
242 | ||
243 | c............................................................................ | |
244 | c normalization and diagonalization | |
245 | c............................................................................ | |
246 | ||
247 | 100 continue | |
248 | ||
249 | !~~~normalize | |
250 | ||
251 | vol= (xmaxico-xminico)/nxico | |
252 | . *(ymaxico-yminico)/nyico | |
253 | . *(zmaxico-zminico)/nzico *tauico | |
254 | fac=1./ntot/vol | |
255 | ||
256 | do ix=1,nxico | |
257 | do iy=1,nyico | |
258 | do iz=1,nzico | |
259 | do i=1,4 | |
260 | do k=1,4 | |
261 | IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)*fac | |
262 | enddo | |
263 | enddo | |
264 | do i=1,3 | |
265 | do k=1,4 | |
266 | IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)*fac | |
267 | enddo | |
268 | enddo | |
269 | enddo | |
270 | enddo | |
271 | enddo | |
272 | ||
273 | !~~~diagonalize T | |
274 | ||
275 | do ix=1,nxico | |
276 | do iy=1,nyico | |
277 | do iz=1,nzico | |
278 | do i=1,4 | |
279 | do k=1,4 | |
280 | T(i,k)=IcoT(i,k,ix,iy,iz) | |
281 | enddo | |
282 | enddo | |
283 | call DiagTmunu(T,u,v,gamma,eps,nonz,ix,iy,iz) | |
284 | if(nonz.eq.1)then | |
285 | IcoE(ix,iy,iz)=eps | |
286 | IcoV(1,ix,iy,iz)=v(1) | |
287 | IcoV(2,ix,iy,iz)=v(2) | |
288 | !rapx=zminico+(iz-0.5)*dzico | |
289 | !rap=rapx+log(gamma*(1+v(3))) | |
290 | !IcoV(3,ix,iy,iz)=tanh(rap) | |
291 | IcoV(3,ix,iy,iz)=v(3) / tauico | |
292 | do i=1,3 | |
293 | IcoF(i,ix,iy,iz)=IcoC(i,4,ix,iy,iz)*u(4) | |
294 | do k=1,3 | |
295 | IcoF(i,ix,iy,iz)=IcoF(i,ix,iy,iz) | |
296 | . -IcoC(i,k,ix,iy,iz)*u(k) | |
297 | enddo | |
298 | enddo | |
299 | endif | |
300 | enddo | |
301 | enddo | |
302 | enddo | |
303 | ||
304 | return | |
305 | end | |
306 | ||
307 | c------------------------------------------------------------------------------------ | |
308 | subroutine DiagTmunu(Tmunu,u,v,gamma,eps,nonz,ix,iy,iz) | |
309 | c------------------------------------------------------------------------------------ | |
310 | ! solve lambda * T * lambda = diag(eps,P,P,P), for symmetric T | |
311 | ! | |
312 | ! lambda = g g*vx g*vy g*vz =symmetric! | |
313 | ! g*vx a*vx*vx+1 a*vy*vx a*vz*vx | |
314 | ! g*vy a*vx*vy a*vy*vy+1 a*vz*vy | |
315 | ! g*vz a*vx*vz a*vy*vz a*vz*vz+1 | |
316 | ! with g=gamma, a=g*g/(g+1) | |
317 | ! | |
318 | ! so T*lambda(v)=lambda(-v)*diag(eps,P,P,P) | |
319 | ! first column: four equations | |
320 | ! eps=T00+T0k*vk | |
321 | ! -eps*vi=Ti0+Tik*vk | |
322 | ! solved iteratively to get eps, vi | |
323 | ! returns u,v,eps if nonz=1 (otherwise zero tensor) | |
324 | ! (by T.Kodama) | |
325 | ! (K. Werner: modifs) | |
326 | ||
327 | include 'epos.incico' | |
328 | double precision Tmunu(4,4), u(4),eps,gamma,g,a, Lor(4,4),w(4),sum | |
329 | double precision v(4),vx(3),tt(4,4),err,sg(4) | |
330 | ||
331 | sg(4)=1d0 | |
332 | v(4)=1d0 | |
333 | do i=1,3 | |
334 | sg(i)=-1d0 | |
335 | enddo | |
336 | sum=0d0 | |
337 | do i=1,4 | |
338 | do k=1,4 | |
339 | sum=sum+dabs(Tmunu(i,k)) | |
340 | end do | |
341 | end do | |
342 | nonz=0 | |
343 | if(sum.eq.0.0d0)return | |
344 | nonz=1 | |
345 | ||
346 | do k=1,3 | |
347 | v(k)=0. | |
348 | end do | |
349 | eps=0 | |
350 | ||
351 | do lrep=1,100 | |
352 | epsx=eps | |
353 | do i=1,3 | |
354 | vx(i)=v(i) | |
355 | end do | |
356 | eps=Tmunu(4,4) | |
357 | do k=1,3 | |
358 | eps=eps-Tmunu(4,k)*vx(k) | |
359 | end do | |
360 | if(eps.le.0d0)then | |
361 | print*,'DiagTmunu: sum(abs(Tmunu))=',sum,' eps=',eps | |
362 | stop'negative energy density. ' | |
363 | endif | |
364 | do i=1,3 | |
365 | Tv=0d0 | |
366 | do k=1,3 | |
367 | Tv=Tv+Tmunu(i,k)*vx(k) | |
368 | end do | |
369 | v(i)=(Tmunu(i,4)-Tv)/eps | |
370 | end do | |
371 | if(lrep.gt.60)then | |
372 | do i=1,3 | |
373 | v(i)=0.5d0*(vx(i)+v(i)) | |
374 | enddo | |
375 | endif | |
376 | !print*,'Tmunu: ',lrep,abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3) | |
377 | err=1d-6 | |
378 | if(lrep.gt.50)err=1d-5 | |
379 | if(lrep.gt.89)err=1d-4 | |
380 | if(abs(eps-epsx).lt.err.and.abs(v(1)-vx(1)).lt.err | |
381 | . .and.abs(v(2)-vx(2)).lt.err.and.abs(v(3)-vx(3)).lt.err)goto1 | |
382 | do i=1,4 | |
383 | w(i)=0 | |
384 | do k=1,4 | |
385 | w(i)=w(i)+Tmunu(i,k)*v(k)*sg(k) | |
386 | enddo | |
387 | w(i)=w(i)-eps*v(i) | |
388 | enddo | |
389 | if(lrep.gt.95 | |
390 | ..and.w(1)*w(1)+w(2)*w(2)+w(3)*w(3)+w(4)*w(4).lt.err)goto1 | |
391 | end do | |
392 | ||
393 | 1 v2=0.d0 | |
394 | do i=1,3 | |
395 | v2=v2+v(i)**2 | |
396 | end do | |
397 | if(v2.ge.1.)stop'DiagTmunu: v2 ge 1. ' | |
398 | gamma=1./sqrt(abs(1.-v2)) | |
399 | u(4)=gamma | |
400 | u(1)=v(1)*gamma | |
401 | u(2)=v(2)*gamma | |
402 | u(3)=v(3)*gamma | |
403 | ||
404 | !~~~check | |
405 | g=gamma | |
406 | a=g*g/(g+1d0) | |
407 | Lor(4,4)=g | |
408 | do k=1,3 | |
409 | Lor(4,k)=-g*v(k) | |
410 | Lor(k,4)=Lor(4,k) | |
411 | enddo | |
412 | do i=1,3 | |
413 | do k=1,3 | |
414 | Lor(i,k)=a*v(i)*v(k) | |
415 | enddo | |
416 | enddo | |
417 | do k=1,3 | |
418 | Lor(k,k)=Lor(k,k)+1 | |
419 | enddo | |
420 | do i=1,4 | |
421 | do k=1,4 | |
422 | tt(i,k)=0d0 | |
423 | do m=1,4 | |
424 | do n=1,4 | |
425 | tt(i,k)=tt(i,k)+Lor(i,m)*Tmunu(m,n)*Lor(n,k) | |
426 | enddo | |
427 | enddo | |
428 | enddo | |
429 | enddo | |
430 | err=err*1d2 | |
431 | if(tt(1,4).gt.err.or.tt(2,4).gt.err.or.tt(3,4).gt.err)then | |
432 | print*,'************** DiagTmunu: nonzero T14 or T24 or T34' | |
433 | . ,' after ',lrep,' iterations' | |
434 | print*,'d_eps or d_v(i): ',abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3) | |
435 | print*,'Tmunu ( ix=',ix-nxico/2-1,' iy=',iy-nyico/2-1,' ) :' | |
436 | . ,' iz=',iz-nzico/2-1 | |
437 | do i=1,4 | |
438 | write(*,'(4f10.5,2x,4f10.5)') | |
439 | . (sngl(tmunu(i,k)),k=1,4),(sngl(tt(i,k)),k=1,4) | |
440 | enddo | |
441 | endif | |
442 | ||
443 | end | |
444 |