coload.f

c $Id: coload.f,v 1.11 2001/04/06 21:48:16 weber Exp $
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function sigtot(ind1,ind2,sqrts)
c
c     Revision : 1.0
C
cinput ind1 : index of particle 1
cinput ind2 : index of particle 2
cinput sqrts: $\sqrt{s}$ of collision between part. 1 and 2
c
c     {\tt sigtot} returns the total cross section (in mbarn) for the collision
c     between the particles with the indices {\tt ind1} and {\tt ind2}.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       implicit none
       include 'coms.f'
       include 'comres.f'
       include 'newpart.f'
c
       integer isigline
       integer ind1,ind2,collclass
       integer ityp1,ityp2,iso31,iso32
       real*8 sqrts,mminit
c for detailed-balance
       integer nCh,ii
       real*8 e1,e2,sigma


c     reset sigtot, store often needed values in scalars
       sigtot=0.d0
       ityp1=ityp(ind1)
       ityp2=ityp(ind2)
       iso31=iso3(ind1)
       iso32=iso3(ind2)

c     now get collision class (line-number for sigmaLN array in blockres.f)
c     isigline classifies the collision (pp,pn,Delta N, Meson-Baryon etc)
       isigline=collclass(ityp1,iso31,ityp2,iso32)

       if(isigline.eq.0) return ! zero cross section for collclass=0

c     get pointer for total cross section (column #2 in SigmaLN array)
       iline=SigmaLn(2,1,isigline) ! flag of total cross section

c     if zero, cross section is zero, return
       if(iline.eq.0) return

       if(iline.gt.0) then
c
c     if gt zero we have a tabulated or parameterized total cross section,
c     all table-lookups and parametrizations are accessed via crossx
c
             call crossx(iline,sqrts,ityp1,iso31,
     &                         max(fmass(ind1),mminit(ityp1)),
     &                         ityp2,iso32,
     &                         max(fmass(ind2),mminit(ityp2)),
     &                         sigtot)

       else
c
c     total cross section via sum of partial cross sections
c
c     get number of exit-channels:
          if (isigline.gt.maxreac) then
             write (6,*) '4isigline: ',isigline
          endif
          nCh=SigmaLn(1,1,isigline)
c
c transformation quantities  into NN system for proper kinematics
c (necessary for detailed balance cross sections)
c     first compute transformation betas
          e1=sqrt(fmass(ind1)**2+px(ind1)**2
     &         +py(ind1)**2+pz(ind1)**2)
          e2=sqrt(fmass(ind2)**2+px(ind2)**2
     &         +py(ind2)**2+pz(ind2)**2)
          betax=(px(ind1)+px(ind2))/(e1+e2)
          betay=(py(ind1)+py(ind2))/(e1+e2)
          betaz=(pz(ind1)+pz(ind2))/(e1+e2)
c     now transform momenta
          pxnn=px(ind1)
          pynn=py(ind1)
          pznn=pz(ind1)
          p0nn=e1
c     call to Lorentz transformation
          call rotbos(0d0,0d0,-betax,-betay,-betaz,
     &         pxnn,pynn,pznn,p0nn)
          pnn=sqrt(pxnn*pxnn+pynn*pynn+pznn*pznn)
c     end of transform part
c
c     loop over exit channels for sum of partial cross sections
c     partial cross sections start in column #3 of SigmaLN in blockres.f
          do 10 ii=3,nCh+2
c     get pointer for partial cross section
             iline=SigmaLn(ii,1,isigline)
c     normal partial cross sections
             if(iline.gt.0) then
                   call crossx(iline,sqrts,ityp1,iso31,fmass(ind1),
     &                         ityp2,iso32,fmass(ind2),sigma)
             else
c     detailed balance partial cross section (must be computed now)
c     (crossz delivers inelastic channel for SINGLE resonance)
                call crossz(iline,sqrts,ityp1,iso31,fmass(ind1),
     &                      ityp2,iso32,fmass(ind2),sigma)
             endif
c     perform summation of partial cross sections
             sigtot=sigtot+sigma
c     end of loop for partial cross sections
 10       continue
c     end of total / sum of partial cross sections if
       endif

c     return to caller
       return

       end


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine getnext(k)
c
c     Revision : 1.0
C
coutput k : index of next collison
c
c {\tt getnext} returns the index of the next collision or decay
c to be performed.
c If no further collisons occur in the timestep, {\tt getnext} returns
c {\tt k}=0 .
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      implicit none
      integer k
      include 'coms.f'
      include 'options.f'
      include 'colltab.f'

c     shrink collision tables in case 80% load to counter possible overflow
       if(dble(nct)/ncollmax.ge.0.8) call collshrink

c
c     set k to current collision/decay (which has already been performed and
c     and is outdated by the time of this call)
      k = actcol

 1    continue
c     increment counter, next entry in table is now the current one
      k = k+1
c     if the end of the collision table has been reached, return with k=0
      if (k.gt.nct) then
         k = 0
         actcol = k
         return
      endif

c     if the current entry in the collision table is marked "F" - false -
c     (due to previous interaction of one of the collision partners)
c     then find new collision partners for the particle(s) via calls
c     to collupd
      if (.not.ctvalid(k)) then
         call collupd(cti1(k),1)
c     second call only if not a decay entry
         if(cti2(k).gt.0) call collupd(cti2(k),1)
c     if the current collision is now marked "T" - true - return
         if(ctvalid(k)) then
            actcol = k
            return
         endif
c     the current collision is still marked false, go to top of loop
c     (and increment counter)
         goto 1
      endif
c
c     the current entry is marked "T" - true - this is the next
c     collision to be perfomed
      actcol = k
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine collshrink
c
c     Revision : 1.0
c
c     This subroutine deletes all entries in the collision tables between
c     1 and {\tt actcol}-1. It's purpose is to counter a possible overflow
c     of the collision tables.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer i
      include 'colltab.f'

      do 104 i=actcol,nct
         cttime(1+i-actcol) = cttime(i)
         ctsqrts(1+i-actcol) = ctsqrts(i)
         ctsigtot(1+i-actcol) = ctsigtot(i)
         cti1(1+i-actcol) = cti1(i)
         cti2(1+i-actcol) = cti2(i)
         ctvalid(1+i-actcol) = ctvalid(i)
         ctcolfluc(1+i-actcol) = ctcolfluc(i)
 104  continue
c     recalculate number of collisions in tables
      nct=1+nct-actcol
c     reset pointer to current collision
      actcol=1

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine colload
c
c     Revision : 1.0
c
c
c This routine fills the collision tables with all collisions and decays
c to be performed in the current timestep. Within the timestep,
c particle propagation is assumed on straight lines. This routine
c actually only performs the outer of the double particle loop and
c calls {\tt collupd} for the inner loop.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      implicit none
      integer i

      include 'coms.f'
      include 'colltab.f'

c     reset number of collisions in table and current collision pointer
      nct = 0
      actcol = 0

c reset all collision arrays
      do 10 i=1,ncollmax
         cttime(i)=0.d0
         ctsqrts(i)=0.d0
         ctsigtot(i)=0.d0
         cti1(i)=0
         cti2(i)=0
         ctvalid(i)=.false.
         ctcolfluc(i)=1.d0
 10   continue

c     outer loop over all particles
      do 20 i=1,npart
c     call collupd for inner loop, -1: inner loop only from i+1 to npart
         call collupd(i,-1)
 20   continue
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine collupd(i,all)
c
c     Revision : 1.0
C
cinput i   : particle to be checked for collison or decay
cinput all : flag for update mode
c
c {\tt collupd} checks whether particle {\tt i} will collide or decay
c in the time interval between the current time {\tt acttime}
c and the end of the time step.
c {\tt collupd} uses the variable {\tt tmin} to find the {\bf earliest}
c interaction/decay of particle {\tt i} and store it in the
c collision arrays via a call to {\tt ctupdate}.
c For {\tt all}$>0$ all other particles from 1 to {\tt npart} are checked
c (necessary for update after a collision/decay), whereas for
c {\tt all}$<0$ only the particles with the indices from {\tt i+1} to
c {\tt npart} are checked (for calls via {\tt colload}).
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      implicit none

      include 'coms.f'
      include 'options.f'
      include 'colltab.f'
      include 'boxinc.f'

      real*8 dst2
      integer i,all
      integer j,imin,jmin,j0,A
      integer wall
      integer stidx
      logical isstable
      real*8  tn, wallcoll
cc
      real*8 tcoll,sqrs,sigt,sqrts,sigtot
      real*8 smin,sigmin,sigfac

      real*8 colfaci,colfacj,colfac,colfacmin

c     number of "initial" particles in event
      A = At+Ap
c     initially tmin is set to the time-interval until the end of the timestep
c     tmin is then minimized to the time of the first interaction of
c     particle i
c
c     acttime: current time
c     time: time at beginning of timestep
c     dtimestep: length of timestep
      tmin=dtimestep-acttime+time
c
c     other information to be stored together with tmin
      smin = 0 ! sqrt(s) of interaction
      sigmin = 0 ! total cross section of interaction
      imin = 0 ! index of first particle
      jmin = 0 ! index of second particle
c
c  check, if particle i is a resonance, that might decay within the remaining
c  part of the timestep.
c  If so, than treat decay as collision (with particle 2= 0)
      if (dectime(i)-acttime.lt.tmin) then
         isstable = .false.
         do 132 stidx=1,nstable
            if (ityp(i).eq.stabvec(stidx)) then
c               write (6,*) 'no decay of particle ',ityp(i)
               isstable = .true.
            endif
 132     enddo
         if (.not.isstable) then
            tmin = dectime(i) - acttime
            smin=fmass(i)
            imin = i
            jmin = 0
         endif
      endif


ccccccccccccccccccccccccccccccccc
c new walls selected if mbflag equals 2
      if ((mbox.gt.0).and.(mbflag.eq.2)) then

c acttime is the ACTUAL time
c time is the time at the BEGINNING of the timestep
c tmin is being minimized RELATIVE to the beginning of the timestep

c tn is the relativ time to the next wall colision
        tn=wallcoll(i,wall)

c comparison wheather a particle decays before a wall colision
            if (tn.lt.tmin) then
c set the time
                    tmin=tn
c set the particle number
                    imin=i
c set the wall
                    jmin=wall
            endif
        endif
cc


c default setting: loop does only go from i+1 to npart
         j0 = i+1
c in case of "update mode" let loop run starting from 1
         if (all.gt.0) j0 = 1
c
c  Now check, which is the earliest collision of particle i
c
         do 101 j=j0,npart

c check for some exclusion cases
            if (
c 1. avoid "self-interaction"
     &           i.ne.j
c 2. particles which have interacted with each other in the past
c    are only allowed to interact with each other if at least one
c    of them has had an interaction in between.
c    (due to string-decays the structure is a little complicated, since
c    one particle can have multiple partners of it's last interaction)
     &          .AND.
     &          ((lstcoll(i).ne.j.and.lstcoll(i).ne.(nmax+strid(j)+1))
     &           .or.
     &           (lstcoll(j).ne.i.and.lstcoll(j).ne.(nmax+strid(i)+1))
     &           .or.
     &           (lstcoll(i).eq.0.and.lstcoll(j).eq.0))
c 3. Particles within the projectile or target respectively are per
c    default only allowed to interact with each other in case they
c    have already had an interaction with a particle of the target
c    or projectile respectively. This can be turned off by setting
c    CTOption(6) to a nonzero value.
c    This rule of course must not not apply to produced particles.
c    In the case of a meson nucleus collision, projectile and
c    target may be swapped in the particle vectors (therefore the
c    use of Apt instead of Ap, because Apt has been then swapped
c    accordingly)
     &          .AND.
     &          (CToption(6).ne.0.or.i.gt.A.or.j.gt.A.or.
     &           ncoll(i)+ncoll(j).gt.0.or.
     &          (i.le.Apt.and.j.gt.Apt).or.(j.le.Apt.and.i.gt.Apt))
     &         ) then

c
c  determine time of minimal approach of particles i and j
c  relative to current time
                  call nxtcoll(i,j,dst2,tcoll)
c
c  are the particles close enough - check cut off
c (default: 250 mbarn, defined in coms.f)
               if (dst2.lt.hit_sphere) then
c does the collision occur in the current time step?
                  if (tcoll.gt.0.d0.and.tcoll.lt.tmin) then
c     get sqrt(s) and total cross section
                     sqrs = sqrts(i,j)
c
c reduced cross section for leading hadrons of string fragmentation
c within their formation time
c the scaling factor is sigfac which is determined by the
c individual particles scaling factors xtotfac
                     if(tform(i).le.acttime+tcoll
     &                       .and.tform(j).le.acttime+tcoll) then
                       sigfac=1.d0
                     else if(tform(i).le.acttime+tcoll
     &                       .and.tform(j).gt.acttime+tcoll) then
                       sigfac=xtotfac(j)
                     else if(tform(j).le.acttime+tcoll
     &                       .and.tform(i).gt.acttime+tcoll) then
                       sigfac=xtotfac(i)
                     else
                       sigfac=xtotfac(i)*xtotfac(j)
                     endif
c     get total cross section via call to sigtot and rescale
                     sigt = sigfac*sigtot(i,j,sqrs)
c     rescale sigtot due to color fluctuations
ctp060202                     call colorfluc(ityp(i),ityp(j),sqrs,
                     call colorfluc(ityp(i),ityp(j),
     &                              colfaci,colfacj)
                     colfac=colfaci*colfacj
                     sigt = sigt*colfac
c     are we within the geometrical cross section
                     if (sigt.gt.max(1.d-8,(10.d0*pi*dst2))) then
c     this collision is to beat, now
                        tmin = tcoll
                        smin = sqrs
                        sigmin = sigt
                        imin = i
                        jmin = j
                        colfacmin=colfac
                     endif
                  endif
               endif
            endif
 101     continue
         if (imin.gt.0) then
c  if we found something, then update table via call to ctupdate
c  (keep in mind: tmin is relative to actual time!)
            call ctupdate(imin,jmin,acttime+tmin,smin,sigmin,colfacmin)
c     in case of "full load mode" after every collision
c     only the first entry in the collision table is relevant
            if(CTOption(17).ne.0) nct=1
         endif

       return
       end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine ctupdate(i,j,t,s,sig,cfac)
c
cinput i : index of 1st colliding particle
cinput j : index of 2nd colliding particle (0 for decay)
cinput t : (absolut) time  of collision/decay
cinput s : $sqrt{s}$ (GeV) of collison
cinput sig : total cross section (mbarn) of collision
cinput cfac :  scaling factor for color fluctuation
c
c This subroutine updates the collision arrays.
c It determines the (chonologically) correct position for the new
c entry in the collision arrays, creates the respective slot and
c inserts the new entry (via a call to {\tt ctset}.
c Then the arrays are scanned to tag the chonologically first collision
c of particle {\tt i} or {\tt j} {\em true} and all subsequent ones
c as {\em false}.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       implicit none

       integer i,j
       real*8 t,s,sig,cfac
       include 'colltab.f'
       integer k,tfound,ncoll1

       ncoll1 = nct + 1
       tfound = ncoll1
c loop over all collisions in table
       do 101 k=actcol+1,nct
c     get the correct position for the new entry (chronologically sorted)
          if (tfound.eq.ncoll1.and.t.le.cttime(k)) tfound = k
c     the above construct works as follows:
c     as long as t > cttime(k), the first term is true and teh
c     the second is false. At the correct position for the new
c     entry BOTH terms are true, for all later times the first
c     term is false in order not to overwrite the value of tfound
 101   continue

c make sure that collision is not already listed in the table
       if (.not.(t.eq.cttime(tfound).and.(i.eq.cti1(tfound).and.
     &      j.eq.cti2(tfound).or.i.eq.cti2(tfound).and.
     &      j.eq.cti1(tfound)))) then
c then create slot for  new entry
             do 102 k=nct,tfound,-1
                cttime(k+1) = cttime(k)
                ctsqrts(k+1) = ctsqrts(k)
                ctsigtot(k+1) = ctsigtot(k)
                cti1(k+1) = cti1(k)
                cti2(k+1) = cti2(k)
                ctvalid(k+1) = ctvalid(k)
                ctcolfluc(k+1) = ctcolfluc(k)
 102         continue
c     increment number of collisions/decays in table
          nct = ncoll1
       endif
c insert new entry into the created slot via call to ctset
       call ctset(tfound,i,j,t,s,sig,cfac)
c
c     only the chonologically FIRST collision of particle i is set to true
c
       do 103 k=actcol+1,nct,1
c     the newly found collision must be omitted in the following sequence
          if (k.eq.tfound) goto 103
c  is there already a collision with i or j ?
          if (cti1(k).eq.i.or.cti2(k).eq.i) then
c     if the other collision is at an earlier time (the table is time-ordered,
c     therefore k < tfound corresponds to an earlier time)
c     then set the new one to false or vice versa
             if (k.lt.tfound.and.ctvalid(k)) then
                ctvalid(tfound) = .false.
             else
                ctvalid(k) = .false.
             endif
          endif
c     do likewise for the second particle
          if (j.gt.0.and.(cti1(k).eq.j.or.cti2(k).eq.j)) then
             if (k.lt.tfound.and.ctvalid(k)) then
                ctvalid(tfound) = .false.
             else
                ctvalid(k) = .false.
             endif
          endif
 103   continue

       return
       end



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine ctset(tfound,i,j,t,s,sig,cfac)
c
c     Revision : 1.0
C
cinput tfound : index in coll. array for coll. to be inserted
cinput i : index of first colliding particle
cinput j : index of second colliing particle (0 for decay, negative: wall)
cinput t : (absolute) time of collision
cinput s : $sqrt{s}$ (GeV) of collison
cinput sig : total cross section (mbarn) of collsion
cinput cfac :  scaling factor for color fluctuation
c
c {\tt ctset} enters the collision of particles {\tt i} and {\tt j}
c into the collision arrays at index {\tt tfound}.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      implicit none

      integer tfound,i,j

      real*8 t,s,sig,cfac
      include 'colltab.f'
c
      cttime(tfound) = t
      ctsqrts(tfound) = s
      ctsigtot(tfound) = sig
      cti1(tfound) = i
      cti2(tfound) = j
      ctvalid(tfound) = .true.
      ctcolfluc(tfound) = cfac
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       real*8 function sqrts(i,j)
c
c     Revision : 1.0
c
c     input: i,j : numbers of colliding particles
c     output: $\sqrt{s}$ of collision as return value
c
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       implicit none

       include 'coms.f'
       integer i,j
       real*8 p10,p20
       p10 = sqrt((px(i)+ffermpx(i))**2
     +           +(py(i)+ffermpy(i))**2
     +           +(pz(i)+ffermpz(i))**2+fmass(i)**2)
       p20 = sqrt((px(j)+ffermpx(j))**2
     +           +(py(j)+ffermpy(j))**2
     +           +(pz(j)+ffermpz(j))**2+fmass(j)**2)
       sqrts = sqrt((p10+p20)**2
     +               -(px(i)+ffermpx(i)+px(j)+ffermpx(j))**2
     +               -(py(i)+ffermpy(i)+py(j)+ffermpy(j))**2
     +               -(pz(i)+ffermpz(i)+pz(j)+ffermpz(j))**2)
       return
       end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine nxtcoll(j,k,dst,dtauc)
c
c     Revision : 1.0
c
c input:  j,k   : indices of colliding particles
coutput dst   : impact parameter squared
coutput dtauc : collisiontime in the computational system
c
c     {\tt nxtcoll}  is the heart  of the collision term. It determines
c     the time  in the  computional system, when the  collision between
c     j and k took or will take place ({\tt dtauc}). The squared impact
c     parameter of the collision is returned in {\tt dst}. {\tt dst} is
c     independent of the computational system  in which the coordinates
c     of j and k are given.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none
      integer j,k, i
      include 'coms.f'
      real*8 u1(0:3), u2(0:3), p1(0:3), p2(0:3)
      real*8 dp(0:3), du(0:3), bt_com(3), bt(3)
      real*8 du2, dp2, dudp, bt_du, bt_dp, btdr, bt2
      real*8 dst, gam_com, gam_com2, bt_com2, dtauc

c
c do some assignments first
c
      u1(0) = r0(j)
      u1(1) = rx(j)
      u1(2) = ry(j)
      u1(3) = rz(j)
      u2(0) = r0(k)
      u2(1) = rx(k)
      u2(2) = ry(k)
      u2(3) = rz(k)

      p1(0) = sqrt(px(j)**2+py(j)**2+pz(j)**2+fmass(j)**2)
      p1(1) = px(j)
      p1(2) = py(j)
      p1(3) = pz(j)
      p2(0) = sqrt(px(k)**2+py(k)**2+pz(k)**2+fmass(k)**2)
      p2(1) = px(k)
      p2(2) = py(k)
      p2(3) = pz(k)
c
c -velocity and gamma-factor of the two particle center of momentum
c frame (com) measured in the computational system
c
      bt_com2 = 0.0d0
      do 1 i=1,3
        bt_com(i) = -((p1(i)+p2(i))/(p1(0)+p2(0)))
        bt_com2 = bt_com2 + bt_com(i)**2
    1 continue
      gam_com = 1.0d0/sqrt(1.0d0-bt_com2)
      gam_com2 = gam_com**2/(1.0d0+gam_com)
c
c calculate some numbers which are needed for the Lorentz-transformation
c
      bt_du = 0.0d0
      bt_dp = 0.0d0
      do 2 i=1,3
        bt_du = bt_du + bt_com(i)*(u1(i) - u2(i))
        bt_dp = bt_dp + bt_com(i)*(p1(i) - p2(i))
    2 continue
c
c calculate bt_com square and the dotproduct bt_com*(r(j)-r(k)),
c where the r's are given in the computational frame
c
c perform Lorentz-transformation of relative distance and relative
c momentum vectors of particles j and k into com-frame
c
c use the resulting 3-vectors du and dp to obtain du and dp squared
c and the dotproduct du*dp
c
      du2  = 0.0d0
      dp2  = 0.0d0
      dudp = 0.0d0
      bt2  = 0.0d0
      btdr = 0.0d0
      do 3 i=1,3
        bt(i) = p1(i)/p1(0) - p2(i)/p2(0) ! Rechensystem rel. vel.
        bt2 = bt2 + bt(i)**2
        btdr = btdr + bt(i)*(u1(i)-u2(i)) ! Rechensystem
        du(i) = u1(i)-u2(i) + bt_com(i)*
     *      (gam_com2*bt_du+gam_com*(u1(0)-u2(0)))
        dp(i) = p1(i)-p2(i) + bt_com(i)*
     *      (gam_com2*bt_dp+gam_com*(p1(0)-p2(0)))
        du2  = du2  + du(i)*du(i)
        dp2  = dp2  + dp(i)*dp(i)
        dudp = dudp + du(i)*dp(i)
    3 continue
c
c  obtain collision time and impact parameter squared
c
      dtauc = -(btdr/bt2)
      dst   = du2 - dudp*dudp/dp2

      end


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine scantab(ind,offs)
c
c     Revision : 1.0
c
cinput ind : index of particle
cinput offs: offset
c
c     {\tt scantab} adjusts the collision table to changed particle
c     indices due to calls to {\tt addpart} or {\tt delpart}.
c
c     In case of an annihilation a list is created of those particles
c     which have lost their collision partner and have to be rechecked
c     for possible collisions/decays.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none
      include 'colltab.f'
      include 'coms.f'
      integer ind,offs,i,k
      logical rescan

c     reset counters for the rechecking of particles due to lost collision
c     partners
      nsav=0
      rescan=.false.

c     call from addpart
      if (offs.gt.0) then
c     shift upwards, if necessary
         do 1 i=1,nct
            if (cti1(i).ge.ind) cti1(i) = cti1(i) + offs
            if (cti2(i).ge.ind) cti2(i) = cti2(i) + offs
 1       continue
c
c     call from delpart
      else

c     omit scan, if last collision in table
         if(actcol.eq.nct) return
c     start loop with next collision in table
         i=actcol+1
 2       continue
         if((cti1(i).eq.ind).or.(cti2(i).eq.ind)) then
c     a dubious entry has been found, now
c     look for particles with 'lost' collision partners
            if(cti1(i).eq.ind.and.cti2(i).gt.0) then
c     save particles with 'lost' partners in the ctsav array
               nsav=nsav+1
               ctsav(nsav)=cti2(i)
               if(ctsav(nsav).gt.ind) ctsav(nsav)=ctsav(nsav) + offs
            elseif(cti2(i).eq.ind.and.cti1(i).gt.0) then
               nsav=nsav+1
               ctsav(nsav)=cti1(i)
               if(ctsav(nsav).gt.ind) ctsav(nsav)=ctsav(nsav) + offs
            endif
c     delete obsolete collision
            do 4 k=i+1,nct
               cttime(k-1) = cttime(k)
               ctsqrts(k-1) = ctsqrts(k)
               ctsigtot(k-1) = ctsigtot(k)
               cti1(k-1) = cti1(k)
               cti2(k-1) = cti2(k)
               ctvalid(k-1) = ctvalid(k)
               ctcolfluc(k-1) = ctcolfluc(k)
 4          continue
c     decrement collision counter
            nct = nct-1
            rescan=.true.
         else
c     else entry is OK
            rescan=.false.
         endif
c     shift rest of table, in case entry has not been deleted
         if(.not.rescan) then
            if (cti1(i).gt.ind) cti1(i) = cti1(i) + offs
            if (cti2(i).gt.ind) cti2(i) = cti2(i) + offs
            i=i+1
         endif
c     condinue procedure until all collisions have been scanned
         if(i.le.nct) goto 2
      endif

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine printtab
c
c     Revision: 1.0
c
c     {\tt printtab} prints the contents of the collision arrays on
c     unit 6 and marks the current collision with an *. This subroutine
c     is used for debugging purposes only.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        implicit none

        include 'colltab.f'
        integer i
        character*1 c
c
        write(6,*) 'colltab:'
        do 101 i=1,nct
           c = ' '
           if (i.eq.actcol) c = '*'
           write(6,'(i4,1x,L1,A1,2(i4,1x),4(f6.3,1x))')
     &            i,ctvalid(i),c,cti1(i),cti2(i),cttime(i),ctsqrts(i)
 101    continue
        return
        end


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ctp060202      subroutine colorfluc(it1,it2,ws,fac1,fac2)
      subroutine colorfluc(it1,it2,fac1,fac2)
c
c     Revision : 1.0
c
cinput it1 : ityp particle 1
cinput it2 : ityp particle 2
cinput ws  : $\sqrt{s}$ of collision
coutput fac1 : x-section scaling factor of particle 1
coutput fac2 : x-section scaling factor of particle 2
c
c     Modifies the hadron cross section due to color fluctuations,
c     ref. L. Frankfurt et al.: Ann. Rev. Nucl. Part. Sc. 44 (1994)501
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none
      include 'options.f'

ctp060202      real*8 ws,fac1,fac2
      real*8 fac1,fac2
      real*8 x,y,Pmes,Pbar,ranf
      integer it1,it2


      fac1=1d0
      fac2=1d0

c switched on ?
      if (ctoption(42).eq.0) return


c particle 1 is baryon
      if (abs(it1).lt.100) then
 1     x=ranf(0)*3d0
       y=ranf(0)

c probability disrtibution for x-section factor
       Pbar=(1.19+32.11*x-15.65*x**2-1.24*x**3+0.94*x**4)/18.d0
       Pbar=max(0d0,Pbar)
       if(y.gt.Pbar) goto 1
       fac1=x
c particle 1 is meson
      else
 2     x=ranf(0)*5d0
       y=ranf(0)

c probability disrtibution for x-section factor
       Pmes=(21.76+4.41*x-3-79*x**2+0.40*x**3)/25.d0
       Pmes=max(0d0,Pmes)
       if(y.gt.Pmes) goto 2
       fac1=x
      endif

c particle 2 is baryon
      if (abs(it2).lt.100) then
 3     x=ranf(0)*3d0
       y=ranf(0)

c probability disrtibution for x-section factor
       Pbar=(1.19+32.11*x-15.65*x**2-1.24*x**3+0.94*x**4)/18.d0
       Pbar=max(0d0,Pbar)
       if(y.gt.Pbar) goto 3
       fac2=x
c particle 2 is meson
      else
 4     x=ranf(0)*5d0
       y=ranf(0)

c probability disrtibution for x-section factor
       Pmes=(21.76+4.41*x-3-79*x**2+0.40*x**3)/25.d0
       Pmes=max(0d0,Pmes)
       if(y.gt.Pmes) goto 4
       fac2=x
      endif

      return
      end