c $Id$ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine bptinit(ibox) c c Unit : Initis all the particles setted by the bpt command c Version : 1.0 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none include 'coms.f' include 'comres.f' include 'boxinc.f' include 'options.f' c Var c counter, spin integer i,ibox,fchg,getspin c randomnumbergenerator, particlemass, deacy times real*8 ranf,massit,dectim c momentum, angle distribution real*8 P,cost,sint,phi ecm=0.d0 ebeam=0.d0 c main program c loop over all particles do 42 i=npart+1,npart+bptpart(ibox) c configuration space r0(i)=0.d0 rx(i)=lboxhalbe*(1-2*ranf(0)) ry(i)=lboxhalbe*(1-2*ranf(0)) rz(i)=lboxhalbe*(1-2*ranf(0)) c isospin and ityp iso3(i)=bptiso3(ibox) ityp(i)=bptityp(ibox) c set baryon and meson numbers if(abs(ityp(i)).le.maxbar) then nbar=nbar+1 else nmes=nmes+1 endif c charge charge(i)=fchg(iso3(i),ityp(i)) c massarray fmass(i)=massit(ityp(i)) c Spin spin(i)=getspin(ityp(i),-1) c decaytime dectime(i)=dectim(i,1) 42 continue c End of loop if (edensflag.le.0) then c homogenious momentum distribution, randomly distributed c max momentum is a parameter do 45 i=npart+1,npart+bptpart(ibox) P=bptpmax(ibox)*ranf(0)**(1./3.) cost = 1.-2.*ranf(0) sint = sqrt(1.-cost**2) phi = 2.*Pi*ranf(0) px(i) = P*sint*cos(phi) py(i) = P*sint*sin(phi) pz(i) = P*cost call setonshell(i) 45 continue elseif (edensflag.ge.1) then c energiedensity c loop over all particles do 60 i=npart+1,npart+bptpart(ibox) P=bptpmax(ibox)/bptpart(ibox)*ranf(0)**(1./3.) cost = 1.-2.*ranf(0) sint = sqrt(1.-cost**2) phi = 2.*Pi*ranf(0) c different initialisations if (para.eq.0) then c Boxmode if (i.eq.1) write(*,*) 'Boxmode' px(i) = P*sint*cos(phi) py(i) = P*sint*sin(phi) pz(i) = P*cost elseif(para.eq.1) then c stream over stream if (i.eq.1) write(*,*) 'streammode' px(i) = 0.d0 py(i) = 0.d0 pz(i) = bptpmax(ibox)/bptpart(ibox)*(-1.d0)**i elseif(para.eq.2) then c slab on slab if (i.eq.1) write(*,*) 'slabmode' px(i)=0.d0 py(i)=0.d0 if (rz(i).gt.0) then pz(i)=(-1.0d0)*bptpmax(ibox)/bptpart(ibox) else pz(i)=bptpmax(ibox)/bptpart(ibox) endif endif 60 continue endif c sum of particles npart=npart+bptpart(ibox) Write(*,*) 'Particles = ',npart cc return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function swapi(x,dx) c c Version: 1.0 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 swapi, x, dx swapi = x 1 if (swapi.lt.-dx) then swapi = swapi + 2.0d0*dx goto 1 end if 2 if (swapi.gt.dx) then swapi = swapi - 2.0d0*dx goto 2 end if end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Function Energie(alpha,max) c c Unit : calculate the energy c Version : 1.0 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none include 'coms.f' integer i real*8 alpha,max real*8 Energie, E E=0 Do 42 i=1,npart E=E+sqrt((alpha**2)*(px(i)*px(i)+py(i)*py(i)+pz(i)*pz(i))+ $ fmass(i)**2) 42 continue Energie=E-max Return End ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Function Regula(me) c c Unit : Searches for the zero of the function c Version : 1.0 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none include 'coms.f' real*8 Regula,xn,xu,x0,me,Energie integer i real*8 E1,E2,E3 c init xu=0.d0 x0=3.d0 c main Write(*,*) 'Regula is running!' i=0 10 Continue i=i+1 E1=Energie(x0,me) E2=Energie(xu,me) xn=x0-(E1*(x0-xu))/(E1-E2) E3=Energie(xn,me) IF ((E2*E3).LE.0) then x0=xn else xu=xn EndIF IF ((ABS(x0-xu).GE.1.D-12).and.(i.le.1000).and.( & ((E3.ge.1.D-12).or.(-E3.ge.1.D-12)))) goto 10 Regula=xn End ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 function wallcoll(ipart,wall) c c Unit : Collisions with an imaginary wall c Version : 1.0 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none include 'coms.f' include 'boxinc.f' include 'options.f' c var real*8 betax,betay,betaz real*8 ty,tz real*8 tn integer wall,ipart integer wally,wallz c Mainprogram c init the variables wall=0 tn=0 c velocity betax=px(ipart)/p0(ipart) betay=py(ipart)/p0(ipart) betaz=pz(ipart)/p0(ipart) c check which wall is reached next and sort by impact time c wall presents the wall and tn the time if (betax.lt.0) then wall=-4 tn=(-lboxhalbe-rx(ipart))/(-max(-betax,1.d-13)) else wall=-1 tn=((lboxhalbe-rx(ipart))/max(betax,1.d-13)) endif if (betay.lt.0) then wally=-5 ty=(-lboxhalbe-ry(ipart))/(-max(-betay,1.d-13)) else wally=-2 ty=((lboxhalbe-ry(ipart))/max(betay,1.d-13)) endif if(ty.lt.tn) then tn=ty wall=wally endif if (betaz.lt.0) then wallz=-6 tz=(-lboxhalbe-rz(ipart))/(-max(-betaz,1.d-13)) else wallz=-3 tz=((lboxhalbe-rz(ipart))/max(betaz,1.d-13)) endif if(tz.lt.tn) then tn=tz wall=wallz endif c sets the time for the earliest wall collision wallcoll=tn return End