EPOS/epos167/epos-ico-158.f

   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