added UrQMD v1.3cr files

Processes/UrQMD/Copyright

+
+Copyright
+
+
+UrQMD source and documentation are provided freely for the purpose of
+checking and reproducing published results of the authors.
+
+The Open Standard Codes and Routines (OSCAR)-Group has established -
+for good reasons - guidelines for reproducablity, usage and quality
+control of simlulations codes for pA and AA collsions.
+
+UrQMD is a complex model. In order to ensure that it is used correctly
+that all results are reproducible and that the proper credits are
+given we ask for your agreement to the following copyright and
+safeguard mechanisms in the OSCAR spirit.
+
+The UrQMD collaboration favors cooperation and joint projects with
+outside researchers. We encourage experimental collaborations to
+compare their results to UrQMD. We support you and/or cooperate on any
+sensible project related to UrQMD
+
+If you are interested in a project, please contact us.
+
+Projects without the participation of the UrQMD-Collaboration are
+accepted, if the project is not a current thesis topic of any
+UrQMD-Collaboration member.
+
+We expect that the code authors are informed about any changes and
+modifications made to the code. Any changes to the official version
+must be documented.
+
+The code or any fragments of it shall not be given away to third
+parties.  Similarily, events generated with UrQMD shall not be given
+to third parties without consent of the code authors.
+

Processes/UrQMD/addpart.f

+c $Id: addpart.f,v 1.4 2000/01/12 16:02:32 bass Exp $
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+C
+      subroutine addpart(index)
+c
+c     Revision : 1.0
+c
+cinput index : index for slot to create in particle arrays
+c
+c This subroutine creates an entry for a particle with index {\tt index} in all
+c particle arrays.
+c \danger{The {\tt nbar} and {\tt nmes} counters must be updated by the
+c calling routine !}
+c
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+      implicit none
+      include 'coms.f'
+      include 'newpart.f'
+      include 'freezeout.f'
+      integer ind,i,j,index
+      ind=index
+
+c     now shift vectors upwards
+      do 10 i=npart,ind,-1
+         r0(i+1)=r0(i)
+         rx(i+1)=rx(i)
+         ry(i+1)=ry(i)
+         rz(i+1)=rz(i)
+
+         p0(i+1)=p0(i)
+         px(i+1)=px(i)
+         py(i+1)=py(i)
+         pz(i+1)=pz(i)
+         fmass(i+1)=fmass(i)
+         ityp(i+1)=ityp(i)
+         iso3(i+1)=iso3(i)
+         ncoll(i+1)=ncoll(i)
+         lstcoll(i+1)=lstcoll(i)
+         charge(i+1)=charge(i)
+         spin(i+1)=spin(i)
+         dectime(i+1)=dectime(i)
+         tform(i+1)=tform(i)
+         xtotfac(i+1)=xtotfac(i)
+         origin(i+1)=origin(i)
+         strid(i+1)=strid(i)
+         uid(i+1)=uid(i)
+         frr0(i+1)=frr0(i)
+         frrx(i+1)=frrx(i)
+         frry(i+1)=frry(i)
+         frrz(i+1)=frrz(i)
+         frp0(i+1)=frp0(i)
+         frpx(i+1)=frpx(i)
+         frpy(i+1)=frpy(i)
+         frpz(i+1)=frpz(i)
+         ffermpx(i+1)=ffermpx(i)
+         ffermpy(i+1)=ffermpy(i)
+         ffermpz(i+1)=ffermpz(i)
+
+         r0_t(i+1)=r0_t(i)
+         rx_t(i+1)=rx_t(i)
+         ry_t(i+1)=ry_t(i)
+         rz_t(i+1)=rz_t(i)
+
+         do 11 j=1,2
+            p0td(j,i+1)=p0td(j,i)
+            pxtd(j,i+1)=pxtd(j,i)
+            pytd(j,i+1)=pytd(j,i)
+            pztd(j,i+1)=pztd(j,i)
+            fmasstd(j,i+1)=fmasstd(j,i)
+            ityptd(j,i+1)=ityptd(j,i)
+            iso3td(j,i+1)=iso3td(j,i)
+ 11      continue
+
+
+ 10      continue
+
+c     increment npart
+         npart=npart+1
+c
+         if(npart.ge.nmax) then
+            write(6,*)'*** error in addpart:too much particles>',nmax
+            write(6,*)' -> increase nmax in coms.f '
+            call file14out(999)
+            stop
+         endif
+
+c     initialize new slot
+         r0(ind)=0.0
+         rx(ind)=0.0
+         ry(ind)=0.0
+         rz(ind)=0.0
+
+         p0(ind)=0.0
+         px(ind)=0.0
+         py(ind)=0.0
+         pz(ind)=0.0
+         fmass(ind)=0.0
+         ityp(ind)=0
+         iso3(ind)=0
+         ncoll(ind)=0
+         lstcoll(ind)=0
+         charge(ind)=0
+         spin(ind)=-1
+         dectime(ind)=0.d0
+         tform(ind)=0.0d0
+         xtotfac(ind)=1.0d0
+         origin(ind)=0
+         strid(ind)=0
+         uid(ind)=0
+         frr0(ind)=0.d0
+         frrx(ind)=0.d0
+         frry(ind)=0.d0
+         frrz(ind)=0.d0
+         frp0(ind)=0.d0
+         frpx(ind)=0.d0
+         frpy(ind)=0.d0
+         frpz(ind)=0.d0
+         ffermpx(ind)=0.d0
+         ffermpy(ind)=0.d0
+         ffermpz(ind)=0.d0
+cpot
+         r0_t(ind)=0.d0
+         rx_t(ind)=0.d0
+         ry_t(ind)=0.d0
+         rz_t(ind)=0.d0
+ctd
+         do 12 j=1,2
+            p0td(j,ind)=0.d0
+            pxtd(j,ind)=0.d0
+            pytd(j,ind)=0.d0
+            pztd(j,ind)=0.d0
+            fmasstd(j,ind)=0.d0
+            ityptd(j,ind)=0
+            iso3td(j,ind)=0
+ 12      continue
+
+c     rescan collision table - all particle indices which have been
+c     shifted upwards must be modified in the collision tables, too.
+      call scantab(ind,1)
+
+c     the lstcoll array must also be shifted due to the new particle slot
+      do 15 i=1,npart
+         if(lstcoll(i).le.nmax) then
+            if(lstcoll(i).eq.ind) then
+               lstcoll(i)=0
+            elseif(lstcoll(i).gt.ind) then
+               lstcoll(i)=lstcoll(i) + 1
+            endif
+         endif
+ 15   continue
+
+
+      return
+      end

Processes/UrQMD/angdis.f

+c $Id: angdis.f,v 1.15 2000/01/12 16:02:33 bass Exp $
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine angdisnew(sqrts,m1,m2,iline,costh,phi)
+c
+c     Revision : 1.1
+c
+c input:  sqrts, m1, m2, iline : characteristics of the ingoing channel
+coutput costh   : cos(theta) of theta-angle
+coutput phi     : phi-angle
+c
+c     {\tt angdisnew} delivers phi and cos(theta) scattering angles
+c     according to the angular distributions given by Mao et al.
+c     {\tt angdisnew} performs numerical inversion of the integral of the
+c     differential cross-section by means of a bisection method.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+
+      real*8 sqrts, m1, m2, costh, costhcoll, phi, ranf, pi, s
+      real*8 anginter
+      integer iline
+      logical symlog(38)
+
+      parameter (pi = 3.14159265358979312d0)
+
+c symmetrize or not angular distribution (depending on iline)
+c this data statements may be changed ...
+      data symlog /14*.true., 1*.false., 6*.true., 7*.false.,
+     &              7*.true., 1*.false., 2*.true./
+
+
+      s = sqrts*sqrts
+      phi = 2.0d0*pi*ranf(0)
+
+      goto (4, 4, 4, 4, 4, 4, 4, 4, 4, 9,
+     &      9, 4, 4, 4, 4, 4, 4, 4, 4, 9,
+c ISO-FB interpolation for MB iline 26/27/28
+     &      4, 4,10,10,10,11,11,11, 4, 4,
+     &      4, 4, 4, 4, 4, 9, 9, 4) abs(iline)
+
+c cross-sections NN (Mao et al.)
+   4  continue
+      costh = -costhcoll(s,m1,m2,symlog(abs(iline)))
+
+      return
+
+c isotropic decay
+   9  continue
+      costh = 1.0d0-2.0d0*ranf(0)
+      return
+
+c no deflection at all
+  10  continue
+      costh = 1.0d0
+      return
+
+c smooth interpolation between iso and f-b:
+  11  continue
+cbl only for intermediate masses, otherwise no deflection (zero degree scattering)
+      if(sqrts.gt.6d0) goto 10
+      costh = anginter()
+      return
+
+      end
+
+      function costhcoll(s,m1,m2,sym)
+      implicit none
+      real*8 costhcoll, s, m1, m2, x, dct, ct, dsigma, ranf
+      integer j
+      logical sym
+
+      x = ranf(0)
+      dct = 2.0d0
+      costhcoll = -1.0d0
+c
+c for jmax=12 the accuracy is better than 0.1 degree
+c
+      do j=1,12 ! accuracy 2**-jmax
+        dct = 0.5d0*dct
+        ct = costhcoll+dct
+        if (dsigma(s,m1,m2,sym,ct).le.x) costhcoll=ct
+      enddo
+c
+c randomize in final interval in order to avoid discrete angles
+c
+      costhcoll = costhcoll+ranf(0)*dct
+      return
+      end
+
+      function dsigma(s_in,m1_in,m2_in,sym,costh)
+c
+cc dsigma(s\_in,chosth) = int_-1^costh dsigma/dOmega(s_in,..) dOmega
+cc it is normalized such that dsigma(s\_in,-1) = 0 and
+cc                            dsigma(s\_in,1) = 1
+      implicit none
+
+      include 'coms.f'
+
+      real*8 dsigma
+      real*8 s_in, m1_in, m2_in, costh,
+     &       msi, cmsi, gsi, mom, cmom, gom, mpi, cmpi, gpi, m
+      real*8 m42, mpi2, cmpi2, d_pi1, d_pi2, cm6gp,
+     &       cpi_3, cpi_2, cpi_1, cpi_m, cpi_l, cpi_0,
+     &       msi2, cmsi2, cmsi4, cmsi6, d_si1, d_si2, d_si3, cm2gs,
+     &       csi_3, csi_2, csi_1, csi_m, csi_l, csi_0,
+     &       mom2, cmom2, cmom4, cmom6, d_om1, d_om2, d_om3, s_om1,
+     &       cm2go, com_3, com_2, com_1, com_m, com_l,
+     &       fac1, d_mx1, d_mx2, d_mx3,
+     &       cmx_o1, cmx_s1, cmx_om, cmx_sm, fac2, fac3,
+     &       cmx_olc, cmx_ols, cmx_slc, cmx_sls
+
+      real*8 sig, tp_pi, tp_si, tp_om, tm_pi, tm_si, tm_om,
+     &       bom_3, bom_2, bom_1, bom_m, bom_0, bom_l,
+     &       bmx_o1, bmx_s1, bmx_om, bmx_sm, bmx_ol, bmx_sl
+
+      real*8 s, tmax, tp, to, twos, brak1, norm,
+     &       t1_pi, t2_pi, t1_si, t2_si, t1_om, t2_om,
+     &       t3_pi, t4_pi, t3_si, t4_si, t3_om, t4_om
+
+      logical firstlog, sym
+      common /nn/ norm
+
+
+c define masses and coupling constants
+      data msi /0.550d0/ cmsi /1.200d0/ gsi  /9.40d0/
+      data mom /0.783d0/ cmom /0.808d0/ gom /10.95d0/
+      data mpi /0.138d0/ cmpi /0.510d0/ gpi  /7.27d0/
+      data m   /0.938d0/
+
+      save firstlog
+
+      save m42, mpi2, cmpi2, d_pi1, d_pi2, cm6gp,
+     &     cpi_3, cpi_2, cpi_1, cpi_m, cpi_l, cpi_0
+
+      save msi2, cmsi2, cmsi4, cmsi6, d_si1, d_si2, d_si3, cm2gs,
+     &     csi_3, csi_2, csi_1, csi_m, csi_l, csi_0
+
+      save mom2, cmom2, cmom4, cmom6, d_om1, d_om2, d_om3, s_om1,
+     &     cm2go, com_3, com_2, com_1, com_m, com_l
+
+      save fac1, d_mx1, d_mx2, d_mx3,
+     &     cmx_o1, cmx_s1, cmx_om, cmx_sm, fac2, fac3,
+     &     cmx_olc, cmx_ols, cmx_slc, cmx_sls
+
+      sig(tp_pi,tp_si,tp_om,tm_pi,tm_si,tm_om) =
+c pion
+     & +((cpi_3*tp_pi + cpi_2)*tp_pi + cpi_1)*tp_pi
+     & + cpi_m/tm_pi + cpi_0 + cpi_l*log(tp_pi*tm_pi)
+c sigma
+     & +((csi_3*tp_si + csi_2)*tp_si + csi_1)*tp_si
+     & + csi_m/tm_si + csi_0 + csi_l*log(tp_si*tm_si)
+c omega
+     & +((bom_3*tp_om + bom_2)*tp_om + bom_1)*tp_om
+     & + bom_m/tm_om + bom_0 + bom_l*log(tp_om*tm_om)
+c mix
+     & + bmx_o1*(tp_om - 1.0d0)
+     & + bmx_s1*(tp_si - 1.0d0)
+     & + bmx_om*log(tm_om)
+     & + bmx_sm*log(tm_si)
+     & + bmx_ol*log(tp_om)
+     & + bmx_sl*log(tp_si)
+
+c calculate constants only once!
+      if(firstlog) goto 1000
+      if (info) write(6,*)
+     $     '(info) dsigma: calculating constants for ang. dist.'
+
+c define constants for pion-Term (no s-dependence)
+      m42   = 4.0d0*m*m
+      mpi2  = mpi*mpi
+      cmpi2 = cmpi*cmpi
+      d_pi1 = cmpi2-mpi2
+      d_pi2 = d_pi1*d_pi1
+      cm6gp = 1.5d0*cmpi2**3*gpi**4*m42*m42/d_pi2
+
+      cpi_3 = -(cm6gp/3.0d0)
+      cpi_2 = -(cm6gp*mpi2/d_pi1)
+      cpi_1 = -(cm6gp*mpi2*(2.0d0*cmpi2 + mpi2)/d_pi2)
+      cpi_m = -(cm6gp*cmpi2*mpi2/d_pi2)
+      cpi_l = -(cm6gp*2.0d0*cmpi2*mpi2*(cmpi2 + mpi2)/d_pi2/d_pi1)
+      cpi_0 = -(cpi_3 + cpi_2 + cpi_1 + cpi_m)
+
+c define constants for sigma-Term (no s-dependence)
+      msi2  = msi*msi
+      cmsi2 = cmsi*cmsi
+      cmsi4 = cmsi2*cmsi2
+      cmsi6 = cmsi2*cmsi4
+      d_si1 = m42-cmsi2
+      d_si2 = m42-msi2
+      d_si3 = cmsi2-msi2
+      cm2gs = 0.5d0*cmsi2*gsi**4/d_si3**2
+
+      csi_3 = -(cm2gs*d_si1**2/3.0d0)
+      csi_2 = -(cm2gs*cmsi2*d_si1*d_si2/d_si3)
+      csi_1 = -(cm2gs*cmsi4*(2.0d0*d_si1 + d_si2)*d_si2/d_si3**2)
+      csi_m = -(cm2gs*cmsi6*d_si2**2/msi2/d_si3**2)
+      csi_l = -(cm2gs*cmsi6*d_si2*(d_si1 + d_si2)*2.0d0/d_si3**3)
+      csi_0 = -(csi_3 + csi_2 + csi_1 + csi_m)
+
+c define constants for omega-Term
+      mom2  = mom*mom
+      cmom2 = cmom*cmom
+      cmom4 = cmom2*cmom2
+      cmom6 = cmom2*cmom4
+      d_om1 = m42-cmom2
+      d_om2 = m42-mom2
+      d_om3 = cmom2-mom2
+      s_om1 = cmom2+mom2
+      cm2go = 0.5d0*cmom2*gom**4/d_om3**2
+
+      com_3 =  cm2go/3.0d0
+      com_2 = -(cm2go*cmom2/d_om3)
+      com_1 =  cm2go*cmom4/d_om3**2
+      com_m =  cm2go*cmom6/(d_om3**2*mom2)
+      com_l = -(cm2go*cmom6*4.0d0/d_om3**3)
+
+c define constants for mix-Term
+      fac1 = -((gsi*gom*cmsi2*cmom2)**2*m42)
+      d_mx1 = cmom2 - cmsi2
+      d_mx2 = cmom2 - msi2
+      d_mx3 = cmsi2 - mom2
+
+      cmx_o1 = fac1/(cmom2*d_mx1**2*d_mx2*d_om3)
+      cmx_s1 = fac1/(cmsi2*d_mx1**2*d_mx3*d_si3)
+      cmx_om = fac1/(d_om3**2*d_mx3**2*(mom2 - msi2))
+      cmx_sm = fac1/(d_si3**2*d_mx2**2*(msi2 - mom2))
+      fac2 = (-fac1)/(d_mx1**3*d_om3**2*d_mx2**2)
+      fac3 = (-fac1)/(d_mx1**3*d_mx3**2*d_si3**2)
+
+      cmx_olc =
+     &  fac2*(3.0d0*cmom2**3        - cmom2**2*cmsi2
+     &      - 2.0d0*cmom2**2*mom2   - 2.0d0*cmom2**2*msi2
+     &      + cmom2*mom2*msi2       + cmsi2*mom2*msi2
+     &      - 4.0d0*cmom2**2*m42    + 2.0d0*cmom2*cmsi2*m42
+     &      + 3.0d0*cmom2*mom2*m42  - cmsi2*mom2*m42
+     &      + 3.0d0*cmom2*msi2*m42  - cmsi2*msi2*m42
+     &      - 2.0d0*mom2*msi2*m42)
+
+      cmx_ols =
+     &  fac2*(8.0d0*cmom2**2        - 4.0d0*cmom2*cmsi2
+     &      - 6.0d0*cmom2*mom2      + 2.0d0*cmsi2*mom2
+     &      - 6.0d0*cmom2*msi2      + 2.0d0*cmsi2*msi2
+     &      + 4.0d0*mom2*msi2)
+
+      cmx_slc =
+     &  fac3*(cmom2*cmsi2**2        - 3.0d0*cmsi2**3
+     &      + 2.0d0*cmsi2**2*mom2   + 2.0d0*cmsi2**2*msi2
+     &      - cmom2*mom2*msi2       - cmsi2*mom2*msi2
+     &      - 2.0d0*cmom2*cmsi2*m42 + 4.0d0*cmsi2**2*m42
+     &      + cmom2*mom2*m42        - 3.0d0*cmsi2*mom2*m42
+     &      + cmom2*msi2*m42        - 3.0d0*cmsi2*msi2*m42
+     &      + 2.0d0*mom2*msi2*m42)
+
+      cmx_sls =
+     &  fac3*(4.0d0*cmom2*cmsi2     - 8.0d0*cmsi2**2
+     &      - 2.0d0*cmom2*mom2      + 6.0d0*cmsi2*mom2
+     &      - 2.0d0*cmom2*msi2      + 6.0d0*cmsi2*msi2
+     &      - 4.0d0*mom2*msi2)
+
+      firstlog = .true.
+      if (info) write(6,*) '(info) dsigma: calculation finished'
+
+c s-dependence beyond this point
+
+ 1000 continue
+      s = s_in - (m1_in+m2_in)**2 + m42
+      tmax = s-m42
+      tp = 0.5d0*(costh+1.0d0)*tmax
+      twos = 2.0d0*s
+
+c define s-dependent stuff for omega-Term
+      brak1 = (twos-m42)**2
+      bom_3 = com_3*(-(2.0d0*cmom2**2) - 2.0d0*cmom2*twos - brak1)
+      bom_2 = com_2*(2.0d0*cmom2*mom2 + s_om1*twos + brak1)
+      bom_1 = com_1*(-(4.0d0*cmom2*mom2) - 2.0d0*mom2**2 -
+     &               2.0d0*(cmom2+2*mom2)*twos - 3.0d0*brak1)
+      bom_m = com_m*(-(2.0d0*mom2**2)- 2.0d0*mom2*twos - brak1)
+      bom_l = com_l*(s_om1*mom2 + (cmom2 + 3.0d0*mom2)*s + brak1)
+      bom_0 = -(bom_3 + bom_2 + bom_1 + bom_m)
+
+c define s-dependent stuff for mix-Term
+      bmx_o1 = cmx_o1*(d_om1 - twos)
+      bmx_s1 = cmx_s1*(d_si1 - twos)
+      bmx_om = cmx_om*(d_om2 - twos)
+      bmx_sm = cmx_sm*(d_si2 - twos)
+      bmx_ol = cmx_olc + cmx_ols*s
+      bmx_sl = cmx_slc + cmx_sls*s
+
+      t1_pi = 1.0d0/(1.0d0+tmax/cmpi2)
+      t2_pi = 1.0d0+tmax/mpi2
+      t1_si = 1.0d0/(1.0d0+tmax/cmsi2)
+      t2_si = 1.0d0+tmax/msi2
+      t1_om = 1.0d0/(1.0d0+tmax/cmom2)
+      t2_om = 1.0d0+tmax/mom2
+
+      norm = sig(t1_pi,t1_si,t1_om,t2_pi,t2_si,t2_om)
+
+      t1_pi = 1.0d0/(1.0d0+tp/cmpi2)
+      t2_pi = 1.0d0+tp/mpi2
+      t1_si = 1.0d0/(1.0d0+tp/cmsi2)
+      t2_si = 1.0d0+tp/msi2
+      t1_om = 1.0d0/(1.0d0+tp/cmom2)
+      t2_om = 1.0d0+tp/mom2
+
+      if (sym) then
+        norm=2.0d0*norm
+        to = tmax-tp
+        t3_pi = 1.0d0/(1.0d0+to/cmpi2)
+        t4_pi = 1.0d0+to/mpi2
+        t3_si = 1.0d0/(1.0d0+to/cmsi2)
+        t4_si = 1.0d0+to/msi2
+        t3_om = 1.0d0/(1.0d0+to/cmom2)
+        t4_om = 1.0d0+to/mom2
+
+        dsigma = (sig(t1_pi,t1_si,t1_om,t2_pi,t2_si,t2_om)
+     &           -sig(t3_pi,t3_si,t3_om,t4_pi,t4_si,t4_om))/norm+0.5d0
+      else
+        dsigma = sig(t1_pi,t1_si,t1_om,t2_pi,t2_si,t2_om)/norm
+      end if
+      return
+      end
+
+      function anginter()
+      implicit none
+      real*8 ranf, anginter, a
+c*
+      a=8d0
+c*
+c the costheta distribution p(x) is proportional to exp(a*x)
+c x is chosen between -1 and +1
+c a=0  corresponds to isotropic distr.
+c a=infinity corresponds to exactly forward distr.
+c inverse transform method is used:
+
+      anginter=1d0/a*log( ranf(0)*(exp(a)-exp(-a))+exp(-a) )
+
+      if(anginter.lt.-1d0.or.anginter.gt.1d0)then
+       write(*,*)"#angdis# illegal costh value: ",anginter
+      endif
+
+      return
+      end
+

Processes/UrQMD/boxinc.f

+c $Id: boxinc.f,v 1.3 1999/01/18 09:56:53 ernst Exp $
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     Unit     : all modules using the box
+c     Version  : 1.0
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c        max count of different species
+        integer bptmax
+         parameter(bptmax=20)
+
+c      counter
+        integer cbox
+c      Flags
+        integer boxflag
+        integer edensflag, mbflag
+        integer para,solid,mtest
+c      number of different species
+        integer mbox
+c      particle counters
+        integer bptpart(bptmax),bptityp(bptmax),bptiso3(bptmax)
+        real*8 bptpmax(bptmax)
+        real*8 edens
+c      edge length of a cube in fm
+        real*8  lbox
+c      half edge lengt of a cube
+        real*8  lboxhalbe
+c      double edge length of a cube
+        real*8  lboxd
+c      momenta
+        real*8 mbp0, mbpx, mbpy, mbpz
+
+        common /boxic/ cbox,boxflag, mbox, bptityp, bptiso3, bptpart
+        common /boxic/ edensflag, para, solid, mbflag, mtest
+        common /boxrc/ lbox, lboxhalbe, lboxd, bptpmax, edens
+        common /boxrc/ mbp0, mbpx, mbpy, mbpz
+

Processes/UrQMD/boxprg.f

+c $Id: boxprg.f,v 1.8 1999/01/18 09:56:53 ernst Exp $
+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
+
+

Processes/UrQMD/cascinit.f

+c $Id: cascinit.f,v 1.11 2002/05/14 12:33:50 balazs Exp $
+      subroutine cascinit(ZZ,AA,nucleus)
+
+      implicit none
+
+      include 'coms.f'
+      include 'options.f'
+      include 'inputs.f'
+
+      integer Z,A,i,getspin,fchg,ZZ,AA,antia, j,k, nrjc,izcnt,nucleus
+      real*8 R,P,R2,ranf,sint,phi,nucrad,cost,drr
+      real*8 xcm,ycm,zcm,pxcm,pycm,pzcm,densnorm,add,avd,ws
+      real*8 r_long, r_short, rdef, ratio, alpha
+c
+c mass correction according to binding energy
+      real*8 meff
+ce bcorr to prevent recycling odd nuclei
+         logical bcorr
+        common /ini/ bcorr
+        bcorr=.false.
+
+C averaged density of one Gaussian in units of central value
+      avd = 0.5**1.5
+
+      pxcm=0.d0
+      pycm=0.d0
+      pzcm=0.d0
+
+      A=abs(AA)
+      Z=abs(ZZ)
+      antia=A/AA
+
+      PT_AA(nucleus)=A
+
+      if(A.gt.AAmax) then
+         write(6,*)'***(E): Mass ',A,' exceeds initialization limit!'
+         write(6,*)'        -> increase parameter AAmax=',AAmax
+         write(6,*)'           in include-file inputs.f '
+         write(6,*)'        -> uqmd aborting ... '
+         stop
+      endif
+
+      if (A.eq.1) then  ! proton or neutron
+        i = 1
+        PT_iso3(i,nucleus) = antia*(-1+2*Z)
+        PT_ityp(i,nucleus) = antia*1
+        uid_cnt=uid_cnt+1
+        PT_uid(i,nucleus)=uid_cnt
+        PT_spin(i,nucleus) = 1
+        PT_dectime(i,nucleus)=1.d31
+        PT_charge(i,nucleus)=fchg(PT_iso3(i,nucleus),
+     &                            PT_ityp(i,nucleus))
+        PT_fmass(i,nucleus) = EMNUC
+        PT_p0(i,nucleus)=sqrt(PT_fmass(i,nucleus)**2)
+        return
+      end if
+
+C
+C Initialisation of A nucleons in configuration space
+C
+      xcm=0.d0
+      ycm=0.d0
+      zcm=0.d0
+      R2=0.d0
+      if (CTOption(24).eq.1) then
+        R2 = nucrad(A) + 10.0
+      elseif(CTOption(24).eq.0)then
+        R2 = nucrad(A)
+      endif
+      nrjc = 0
+      if(CTOption(24).eq.2)then ! use fast method
+        call nucfast(A,nucleus)
+        do j=1,A
+          PT_dectime(j,nucleus)=1.d31
+          xcm=xcm+PT_rx(j,nucleus)
+          ycm=ycm+PT_ry(j,nucleus)
+          zcm=zcm+PT_rz(j,nucleus)
+        enddo
+      else
+      do 1 j=1,A
+        PT_dectime(j,nucleus)=1.d31
+111     nrjc = nrjc+1
+        R = R2*ranf(0)**(1./3.)
+        cost = 1.-2.*ranf(0)
+        sint = sqrt(1.-cost**2)
+        phi = 2.*Pi*ranf(0)
+        PT_r0(j,nucleus) = 0.d0
+        PT_rx(j,nucleus) = R*sint*cos(phi)
+        PT_ry(j,nucleus) = R*sint*sin(phi)
+        PT_rz(j,nucleus) = R*cost
+        if (CTParam(21).gt.0.0) then
+           ratio = sqrt((1 + 4.0*CTParam(21)/3.0) /
+     $                  (1 - 2.0*CTParam(21)/3.0) )
+           alpha = atan( sqrt(PT_rx(j,nucleus)*PT_rx(j,nucleus)
+     $     + PT_ry(j,nucleus)*PT_ry(j,nucleus))/PT_rz(j,nucleus))
+           r_short = nucrad(A)*(ratio**(-(1.0/3.0)))
+           r_long  = nucrad(A)*(ratio**(2.0/3.0))
+           rdef = r_short/sqrt(1-((r_long**2-r_short**2)/r_long**2)
+     $            *cos(alpha)*cos(alpha))
+        else
+           rdef = nucrad(A)
+        endif
+
+
+        if (CTOption(24).eq.1) then
+           WS = 1 / ( 1 + dexp( ( R - rdef ) / 0.545 ) )
+cdh        if (ranf(0) > WS ) then
+           if (ranf(0) .gt. WS ) then
+c               write (6,*) "rejected: ",R
+              goto 111
+           endif
+        else
+           do 11 k=1,j-1
+              drr=(PT_rx(j,nucleus)-PT_rx(k,nucleus))**2
+     &             +(PT_ry(j,nucleus)-PT_ry(k,nucleus))**2
+     &             +(PT_rz(j,nucleus)-PT_rz(k,nucleus))**2
+              if (drr.lt.2.6.and.nrjc.lt.CTParam(46)) goto 111
+ 11        continue
+        endif
+        xcm=xcm+PT_rx(j,nucleus)
+        ycm=ycm+PT_ry(j,nucleus)
+        zcm=zcm+PT_rz(j,nucleus)
+1     continue
+      endif
+
+      if (nrjc.ge.CTParam(46)) then
+c       write(6,*)'*** warning: initialisation corrupt '
+         bcorr=.true.
+      end if
+
+      xcm = xcm/dble(A)
+      ycm = ycm/dble(A)
+      zcm = zcm/dble(A)
+      do 13 j=1,A
+        if (CTOption(24).eq.0) then
+          PT_rx(j,nucleus) = PT_rx(j,nucleus)-xcm
+          PT_ry(j,nucleus) = PT_ry(j,nucleus)-ycm
+          PT_rz(j,nucleus) = PT_rz(j,nucleus)-zcm
+        endif
+        PT_rho(j,nucleus) = avd
+13    continue
+
+C local proton density in nucleus A,Z
+      do 14 j=1,Z
+        do 15 k=j+1,Z
+          drr=(PT_rx(j,nucleus)-PT_rx(k,nucleus))**2
+     &       +(PT_ry(j,nucleus)-PT_ry(k,nucleus))**2
+     &       +(PT_rz(j,nucleus)-PT_rz(k,nucleus))**2
+          add=exp(-(2.0*gw*drr))
+          PT_rho(j,nucleus) = PT_rho(j,nucleus)+add
+          PT_rho(k,nucleus) = PT_rho(k,nucleus)+add
+15      continue
+14    continue
+
+C local neutron density in nucleus A,Z
+      do 16 j=Z+1,A
+        do 17 k=j+1,A
+          drr=(PT_rx(j,nucleus)-PT_rx(k,nucleus))**2
+     &       +(PT_ry(j,nucleus)-PT_ry(k,nucleus))**2
+     &       +(PT_rz(j,nucleus)-PT_rz(k,nucleus))**2
+          add=exp(-(2.0*gw*drr))
+          PT_rho(j,nucleus) = PT_rho(j,nucleus)+add
+          PT_rho(k,nucleus) = PT_rho(k,nucleus)+add
+17      continue
+16    continue
+
+      densnorm = (2.0*gw/pi)**1.5
+      do 18 j=1,A
+        PT_rho(j,nucleus) = PT_rho(j,nucleus)*densnorm
+        PT_pmax(j,nucleus) = hqc*(3.0*pi*pi*PT_rho(j,nucleus))**(1./3.)
+18    continue
+
+      izcnt=0
+      do 12 j=1,A
+         P = PT_pmax(j,nucleus)*ranf(0)**(1./3.)
+         cost = 1.-2.*ranf(0)
+         sint = sqrt(1.-cost**2)
+         phi = 2.*Pi*ranf(0)
+         PT_px(j,nucleus) = P*sint*cos(phi)
+         PT_py(j,nucleus) = P*sint*sin(phi)
+         PT_pz(j,nucleus) = P*cost
+         pxcm=pxcm+PT_px(j,nucleus)
+         pycm=pycm+PT_py(j,nucleus)
+         pzcm=pzcm+PT_pz(j,nucleus)
+         if (j.le.Z) then
+            PT_iso3(j,nucleus)= 1*antia
+            PT_charge(j,nucleus)=1*antia
+         else
+            PT_iso3(j,nucleus)= -(1*antia)
+            PT_charge(j,nucleus)=0
+         endif
+
+         PT_spin(j,nucleus) = getspin(1,-1)
+         PT_ityp(j,nucleus) = 1*antia
+        uid_cnt=uid_cnt+1
+        PT_uid(j,nucleus)=uid_cnt
+
+12    continue
+
+c perform CM-correction
+      pxcm=pxcm/A
+      pycm=pycm/A
+      pzcm=pzcm/A
+      do 2 i=1,A
+         PT_px(i,nucleus)=PT_px(i,nucleus)-pxcm
+         PT_py(i,nucleus)=PT_py(i,nucleus)-pycm
+         PT_pz(i,nucleus)=PT_pz(i,nucleus)-pzcm
+c effective masses for initial energy corr. (CTOption(11).eq.0)
+         r=sqrt(PT_rx(i,nucleus)**2+PT_ry(i,nucleus)**2
+     &         +PT_rz(i,nucleus)**2)
+         p=sqrt(PT_px(i,nucleus)**2+PT_py(i,nucleus)**2
+     &         +PT_pz(i,nucleus)**2)
+ctp060202         PT_fmass(i,nucleus) = meff(z,a,r,p)
+         PT_fmass(i,nucleus) = meff(z,a,p)
+         PT_p0(i,nucleus)=sqrt(PT_px(i,nucleus)**2+PT_py(i,nucleus)**2
+     &                    +PT_pz(i,nucleus)**2+PT_fmass(i,nucleus)**2)
+ 2    continue
+c end of CM-correction
+
+      return
+      end
+
+      subroutine nucfast(IA,JJ)
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+
+      include 'inputs.f'
+      real*8 nucrad
+
+      am=0.545
+      rad=nucrad(ia)/am
+      cr1=1.+3./rad+6./rad**2+6./rad**3
+      cr2=3./rad
+      cr3=3./rad+6./rad**2
+      rimmax=0.0
+
+      DO I=1,IA
+ 1      zuk=ranf(0)*cr1-1.
+        if(zuk.le.0.)then
+          tt=rad*(ranf(0)**.3333-1.)
+        elseif(zuk.le.cr2 )then
+          tt=-log(ranf(0))
+        elseif(zuk.lt.cr3 )then
+          tt=-log(ranf(0))-log(ranf(0))
+        else
+          tt=-log(ranf(0))-log(ranf(0))-log(ranf(0))
+        endif
+        if(ranf(0).gt.1./(1.+exp(-abs(tt))))goto 1
+        rim=tt+rad
+        z=rim*(2D0*ranf(0)-1D0)
+        rim=dsqrt(rim*rim-z*z)
+cc        if(rimmax.lt.rim*AM)rimmax=rim*AM
+        PT_r0(I,JJ)=0d0
+        PT_rz(I,JJ)=Z * AM
+
+ 2      s1=2d0*ranf(0)-1d0
+        s2=2d0*ranf(0)-1d0
+        s3=s1*s1+s2*s2
+        if(s3.gt.1d0)goto 2
+        s3=dsqrt(s3)
+        c=s1/s3
+        s=s2/s3
+        PT_rx(I,JJ)=rim * C * AM
+        PT_ry(I,JJ)=rim * S * AM
+      enddo
+cc      print *,"rimmax",rimmax
+cc      if(debug.ge.3)then
+cc        write (*,*) "nucleons"
+cc        do i=1,ia
+cc          write (*,'(i3,3g12.6)')i,PT_rx(i,JJ),PT_ry(i,JJ),PT_rz(i,JJ)
+cc        enddo
+cc        write (*,*)
+cc      endif
+      return
+      END
+
+
+      function nucrad(AA)
+      implicit none
+      real*8 nucrad, r_0
+      integer A,AA
+      include 'coms.f'
+      include 'options.f'
+
+      A=abs(AA)
+c root mean square radius of nucleus of mass A
+c r_0 corresponding to rho0
+      if (CTOption(24).ge.1) then
+c root mean square radius of nucleus of mass A (Mayer-Kuckuck)
+         nucrad = 1.128 * a**(1./3.) - 0.89 * a**(-(1./3.))
+      else
+         r_0 = (0.75/pi/rho0)**(1./3.)
+c subtract gaussian tails, for distributing centroids correctly
+         nucrad = r_0*(0.5*(a + (a**(1./3.)-1.)**3.))**(1./3.)
+      endif
+
+      return
+      end
+
+      subroutine boostnuc(i1,i2,pin,b,dst)
+      implicit none
+      include 'coms.f'
+      include 'options.f'
+      integer i1,i2,i
+      real*8 b,dst,ei,ti
+      real*8 pin,beta,gamma
+
+      do 1 i=i1,i2
+
+      beta = pin/sqrt(pin**2+fmass(i)**2)
+      gamma = 1.d0/sqrt(1.d0-beta**2)
+
+c  Gallilei-Trafo in x-direction (impact parameter)
+c  projectile hits at POSITIVE x
+         rx(i) = rx(i) + b
+c  distance between nuclei: projectile at NEGATIVE z for dst < 0
+         if(CTOption(23).eq.0)then
+           ti = r0(i)
+           rz(i) = rz(i)/gamma+dst/gamma
+         else
+          rz(i) = (rz(i) + dst)
+         end if
+
+
+         Ei = p0(i)
+         p0(i) = gamma*(p0(i) - beta*pz(i))
+         pz(i) = gamma*(pz(i) - beta*Ei)
+
+ 1    continue
+      return
+      end
+
+ctp060202      real*8 function meff(z,a,r,p)
+      real*8 function meff(z,a,p)
+c mean binding energy of a nucleon in a nucleus according to weizaecker
+      implicit none
+      include 'options.f'
+      real*8 av,as,ac,aa,ap,mdef,p,e,EMNUC
+      integer z,a
+      parameter (av=0.01587,as=0.01834,ac=0.00071)
+      parameter (aa=0.09286,ap=11.46,EMNUC=0.938)
+      if(CTOption(11).ne.0.or.a.eq.1)then
+        meff=EMNUC
+        return
+      end if
+c...mass defect
+      mdef=-(av*A)+as*A**0.66667+ac*z*z/a**0.33333+aa*(z-a/2.)**2/a
+c...energy per nucleon = binding energy + nucleon mass
+      e=min(0d0,mdef/a)+EMNUC
+      meff=sqrt(e**2-p**2)
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine getnucleus(nucleus,offset)
+c
+c     Revision: 1.0
+c
+cinput nucleus : 1=projectile, 2=target
+cinput offset  : offset for location of nucleus in particle vectors
+c
+c output : via common blocks
+c
+c This subroutine read in a nucleus which has been initialized
+c by {\tt cascinit} and stored in the {\tt PT\_ *(i,nucleus)} arrays.
+c The respective nucleus is then rotated randomly in configuration
+c and momentum space to yield a new initial state.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'coms.f'
+      include 'options.f'
+      include 'inputs.f'
+
+c local variables
+      real*8 eul1,eul2,eul3,ceul1,seul1,ceul2,seul2,ceul3,seul3
+      real*8 vecx0,vecy0,vecz0,vecx1,vecy1,vecz1,vecx2
+      real*8 vecy2,vecz2
+      integer k,nucleus,offset
+
+c functions
+      real*8 ranf
+
+c ***
+c *** rotation around Euler-angles
+c ***
+      if (CTParam(21).gt.0.0) then
+         eul1 = 0.0
+         eul2 = 0.0
+         eul3 = 0.0
+      else
+         eul1 = ranf(0)*2.0*pi
+         eul2 = ranf(0)*2.0*pi
+         eul3 = ranf(0)*2.0*pi
+      endif
+
+      ceul1 = cos(eul1)
+      seul1 = sin(eul1)
+      ceul2 = cos(eul2)
+      seul2 = sin(eul2)
+      ceul3 = cos(eul3)
+      seul3 = sin(eul3)
+
+      do 178 k=1,PT_AA(nucleus)
+
+c rotate in configuration space
+
+         vecx0 = PT_rx(k,nucleus)
+         vecy0 = PT_ry(k,nucleus)
+         vecz0 = PT_rz(k,nucleus)
+
+
+         vecx1 =   ceul1*vecx0 + seul1*vecy0
+         vecy1 = -(seul1*vecx0)+ ceul1*vecy0
+         vecz1 =                                     vecz0
+
+         vecx2 =   ceul2*vecx1               + seul2*vecz1
+         vecy2 =                       vecy1
+         vecz2 = -(seul2*vecx1)              + ceul2*vecz1
+
+         vecx0 =   ceul3*vecx2 + seul3*vecy2
+         vecy0 = -(seul3*vecx2)+ ceul3*vecy2
+         vecz0 =                                     vecz2
+
+         rx(k+offset) = vecx0
+         ry(k+offset) = vecy0
+         rz(k+offset) = vecz0
+
+c rotate in momentum space
+
+         vecx0 = PT_px(k,nucleus)
+         vecy0 = PT_py(k,nucleus)
+         vecz0 = PT_pz(k,nucleus)
+
+
+         vecx1 =   ceul1*vecx0 + seul1*vecy0
+         vecy1 = -(seul1*vecx0)+ ceul1*vecy0
+         vecz1 =                                     vecz0
+
+         vecx2 =   ceul2*vecx1               + seul2*vecz1
+         vecy2 =                       vecy1
+         vecz2 = -(seul2*vecx1)              + ceul2*vecz1
+
+         vecx0 =   ceul3*vecx2 + seul3*vecy2
+         vecy0 = -(seul3*vecx2) + ceul3*vecy2
+         vecz0 =                                     vecz2
+
+         px(k+offset) = vecx0
+         py(k+offset) = vecy0
+         pz(k+offset) = vecz0
+
+c initialize the other quantum numbers
+
+         iso3(k+offset)=PT_iso3(k,nucleus)
+         ityp(k+offset)=PT_ityp(k,nucleus)
+         uid(k+offset)=PT_uid(k,nucleus)
+         spin(k+offset)=PT_spin(k,nucleus)
+         dectime(k+offset)=PT_dectime(k,nucleus)
+         charge(k+offset)=PT_charge(k,nucleus)
+         fmass(k+offset)=PT_fmass(k,nucleus)
+         r0(k+offset)=PT_r0(k,nucleus)
+         p0(k+offset)=PT_p0(k,nucleus)
+
+ 178  continue
+
+      return
+      end
+
+C####C##1#########2#########3#########4#########5#########6#########7##
+      real*8 function rnfWSX(AA,zmin,zmax)
+c  yields a $x^n$ distributet value for $x$ between mmin and mmax
+cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
+
+      implicit none
+
+      real*8 zmin,zmax,ranf,a,rr,yf,z,nucrad
+      integer AA
+      parameter (a=0.54d0)
+
+      rr=nucrad(aa)
+
+      if(zmax.lt.rr)then
+        write(6,*)'rnfwsx: maximum radius seems too low'
+        stop
+      end if
+
+ 108  continue
+      z=zmin+ranf(0)*(zmax-zmin)
+      yf=z*z/((zmax-zmin)**3)*0.5d0/(1d0+exp(z-rr)/a)
+      rnfWSX=z
+      if(yf.gt.1d0)stop'rnfWSX: wrong normalisaton:'
+      if(yf.gt.ranf(0))    return
+      goto 108
+
+      end
+
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+      subroutine setonshell(i)
+c
+c     Revision : 1.0
+c     This subroutine set particle i on-shell
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      include 'coms.f'
+      integer i
+
+      p0(i) = sqrt(px(i)**2+py(i)**2+pz(i)**2+fmass(i)**2)
+
+      return
+      end
+

Processes/UrQMD/colltab.f

+c $Id: colltab.f,v 1.5 1999/11/24 19:47:48 ssoff Exp $
+c
+cdes  This file contains the  uqmd collision tables
+c
+      integer ncollmax
+      parameter (ncollmax = 100) ! maximum number of entries in collision table
+      integer nct,actcol,nsav,apt
+      real*8 cttime(0:ncollmax),ctsqrts(ncollmax),ctsigtot(ncollmax)
+      real*8 ctcolfluc(ncollmax)
+      logical ctvalid(ncollmax)
+      real*8 tmin
+      integer cti1(ncollmax),cti2(ncollmax),ctsav(ncollmax)
+c      integer updi1(ncollmax),updi2(ncollmax)
+c
+c     cttime  : collision time
+c     ctsqrts : $sqrt{s}$ of collision
+c     ctsigtot: total cross section in mbarn
+c     tmin    : paramteter for {\tt collupd}
+c     cti1    : index of particle 1
+c     cti2    : index of particle 2
+c     nct     : number of collisions in the table
+c     actcol  : current collision
+c     ctvalid : tag whether collision is {\em true} or {\em false}
+c     ctsav   : list of particles which lost their collision partner
+c     nsav    : number of entries in {\tt ctsav}
+c     apt     : mass of first particle/composite in the part. arrays
+
+      common /colltab/cttime,ctsqrts,ctsigtot,tmin,cti1,cti2,nct,actcol,
+     &     ctvalid,ctsav,nsav,apt,ctcolfluc

Processes/UrQMD/comnorm.f

+c $Id: comnorm.f,v 1.5 1999/01/18 09:56:55 ernst Exp $
+      integer n
+      parameter (n = 400)
+      real*8 x_norm(0:3,1:n),y_norm(0:3,1:n)
+      real*8 y2a(0:3,1:n),y2b(0:3,1:n),dx
+      common /normsplin/ x_norm,y_norm,y2a,y2b,dx

Processes/UrQMD/comres.f

+c $Id: comres.f,v 1.15 2003/06/29 14:26:36 weber Exp $
+c
+cdes This file contains definitions for the collision term
+c
+      integer maxbar,maxbra,minbar
+      integer offmeson,maxmeson,pimeson,maxbrm,minnuc,mindel
+      integer maxbrs1,maxbrs2
+      integer numnuc,numdel,nucleon,maxnuc,maxdel
+      integer minmes,maxmes
+
+
+      parameter (minnuc=1) ! lowest baryon particle ID
+      parameter (minmes=100) ! lowest meson particle ID
+      parameter (maxmes=132) ! hightest meson particle ID
+
+c number of resonances of a kind
+      parameter (numnuc=16) ! number of nucleon resonances
+      parameter (numdel=10) ! number of delta resonances
+c indices of minimal and maximal itype of a kind (redundant but nice)
+      parameter (maxnuc=minnuc+numnuc-1) ! highest nucleon ID
+      parameter (mindel=minnuc+maxnuc)   ! lowest delta ID
+      parameter (maxdel=mindel+numdel-1) ! highest delta ID
+
+c minres & maxres define the range of nonstable & nonstrange baryons
+      integer minres,maxres
+      parameter (minres=minnuc+1) ! lowest baryon resonance ID
+      parameter (maxres=maxdel)   ! highest (nonstrange) baryon
+                                  ! resonance ID
+
+c strangenes.ne.0 baryon resonances
+      integer minlam,minsig,mincas,minome
+      integer numlam,numsig,numcas,numome
+      integer maxlam,maxsig,maxcas,maxome
+      parameter (numlam=13) ! number of lambda states
+      parameter (numsig=9)  ! number of sigma states
+      parameter (numcas=6)  ! number of cascade states
+      parameter (numome=1)  ! number of omega states
+      parameter (minlam=mindel+numdel)   ! ID of lowest lambda state
+      parameter (maxlam=minlam+numlam-1) ! ID of highest lambda state
+      parameter (minsig=minlam+numlam)   ! ID of lowest sigma state
+      parameter (maxsig=minsig+numsig-1) ! ID of highest sigma state
+      parameter (mincas=minsig+numsig)   ! ID of lowest cascade state
+      parameter (maxcas=mincas+numcas-1) ! ID of highest cascade state
+      parameter (minome=mincas+numcas)   ! ID of lowest omega state
+      parameter (maxome=minome+numome-1) ! ID of highest omega state
+
+c minbar & maxbar define the range of all baryons
+      parameter (minbar=minnuc) ! ID of lowest baryon state
+      parameter (maxbar=maxome) ! ID of highest baryon state
+
+      parameter (offmeson=minmes) ! offset between zero and lowest
+                                  ! meson state
+      parameter (maxmeson=maxmes) ! ID of highest meson state
+c... these variables are in principal obsolete and should be exchanged
+c were referenced
+
+c... avoid hard coded itypes
+      integer itrho,itome,iteta,itkaon,itphi,itetapr
+      parameter (itkaon=106)   ! ID of kaon
+      parameter (itrho=104)    ! ID of rho meson
+      parameter (itome=103)    ! ID of omega meson
+      parameter (iteta=102)    ! ID of eta
+      parameter (itphi=109)    ! ID of phi
+      parameter (itetapr=107)  ! ID of eta'
+      parameter (pimeson=101)  ! ID of $\pi$
+      parameter (nucleon=minnuc) ! ID of nucleon
+
+      integer itmin,itmax
+      parameter (itmin=minnuc)  ! lowest defined ID
+      parameter (itmax=maxmes)  ! highest defined ID
+c
+      parameter (maxbra=11)  ! decay channels for $s=0$ baryon resonances
+      parameter (maxbrm=25) ! decay channels for meson resonances
+      parameter (maxbrs1=10)! decay channels for $s=1$ baryon resonances
+      parameter (maxbrs2=3) ! decay channels for $s=2$ baryon resonances
+
+c
+       integer mlt2it(maxmes-minmes) ! meson IDs sorted by multipletts
+
+
+      real*8 massoff,mresmin,mresmax
+      parameter (massoff=1d-4)      ! offset for mass generation
+      parameter (mresmin=1.0765d0)  ! minimum baryon resonance mass
+      parameter (mresmax=5d0)       ! maximum baryon resonance mass
+
+      character*45 versiontag
+      common /versioning/ versiontag
+
+      real*8 massres(minbar:maxbar),widres(minbar:maxbar)
+      real*8 branmes(0:maxbrm,minmes+1:maxmes)
+      real*8 branres(0:maxbra,minnuc+1:maxdel)
+      real*8 branbs1(0:maxbrs1,minlam+1:maxsig)
+      real*8 branbs2(0:maxbrs2,mincas+1:maxcas)
+      integer Jres(minbar:maxbar)
+      integer Jmes(minmes:maxmes)
+      integer pares(minbar:maxbar),pames(minmes:maxmes)
+      integer Isores(minbar:maxbar), Isomes(minmes:maxmes)
+      integer brtype(4,0:maxbra),bmtype(4,0:maxbrm)
+      integer bs1type(4,0:maxbrs1),bs2type(4,0:maxbrs2)
+      real*8 massmes(minmes:maxmes)
+      real*8 mmesmn(minmes:maxmes)
+      real*8 widmes(minmes:maxmes)
+      integer strres(minbar:maxbar),strmes(minmes:maxmes)
+
+      integer lbr(0:maxbra,minnuc+1:maxdel)
+      integer lbs1(0:maxbrs1,minlam+1:maxsig)
+      integer lbs2(0:maxbrs2,mincas+1:maxcas)
+      integer lbm(0:maxbrm,minmes+1:maxmes)
+
+      common /resonances/ massres,widres,massmes,widmes,mmesmn,
+     ,                    branres,branmes,branbs1,branbs2,
+     ,                    bs1type,bs2type,lbs1,lbs2,lbm,
+     ,                    jres,jmes,lbr,brtype,pares,pames,
+     ,                    bmtype,
+     ,                    Isores,Isomes,strres,strmes,mlt2it
+
+c     massres   : baryon mass table
+c     widres    : baryon decay width table
+c     massmes   : meson mass table
+c     widmes    : meson decay width table
+c     mmesmn    : table of minimum masses for meson resonances
+c     branres   : branching ratios for $s=0$ baryon resonances
+c     branmes   : branching ratios for meson resonances
+c     branbs1   : branching ratios for $s=1$ baryon resonances
+c     branbs2   : branching ratios for $s=2$ baryon resonances
+c     brtype    : definitions of the decay branches for $s=0$ baryon resonances
+c     bmtype    : definitions of the decay branches for meson resonances
+c     bs1type   : definitions of the decay branches for $s=1$ baryon resonances
+c     bs2type   : definitions of the decay branches for $s=2$ baryon resonances
+c     lbr       : decay angular momenta for $s=0$ baryon decays
+c     lbm       : decay angular momenta for meson decays
+c     lbs1      : decay angular momenta for $s=1$ baryon decays
+c     lbs2      : decay angular momenta for $s=2$ baryon decays
+c     jres      : spin table for baryons
+c     jmes      : spin table for mesons
+c     pares     : parity table for baryons
+c     pames     : parity table for mesons
+c     isores    : isospin table for baryons
+c     isomes    : isospin table for mesons
+c     strres    : strangeness table for baryons
+c     strmes    : strangeness table for mesons
+c
+c
+ccccccccccccccccccccc sigtab-declarations cccccccccccccccccccccccccccccccccc
+
+      integer itblsz,nsigs,maxreac,maxpsig,sigver
+c     VERSION NUMBER of SIGTAB
+      parameter (sigver = 1000) ! version number for collision
+                                ! term and tables
+ccccccccccccccccccccccccccccccccccccccc
+c
+
+      parameter (maxreac = 13) ! maximum number of collision classes
+      parameter (maxpsig = 12) ! maximum number of cross
+                               ! sections per class
+      parameter (nsigs = 10)   ! number of tabulated cross sections
+
+      PARAMETER (ITBLSZ= 100)  ! table size of cross section tables
+
+
+c
+      integer sigmaLN(maxpsig,2,maxreac)
+      integer sigmainf(nsigs,20)
+      real*8 sigmascal(nsigs,5), sigmas(nsigs,itblsz)
+c
+      common/sigtabi/sigmaln,sigmainf
+      common/sigtabr/sigmas,sigmascal
+
+c     sigmaln   : pointer array connecting collision classes with cross sections
+c     sigmainf  : information related to tabulated cross sections
+c     sigmascal : information related to tabulated cross sections
+c     sigmas    : tabulated cross sections

Processes/UrQMD/coms.f

+c $Id: coms.f,v 1.21 2003/06/29 14:26:36 weber Exp $
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c standard common block for uQMD
+c
+cdes This file contains the standard commom blocks of UrQMD
+c
+      real*8 emnuc
+      parameter (emnuc = 0.938) ! nucleon mass
+
+      integer nmax, nspl
+      real*8 hit_sphere
+      parameter (nmax = 500) ! maximum number of particles
+      parameter (nspl = 500)  ! dimension of spline arrays
+      parameter (hit_sphere = 8.d0)  ! hard collision cutoff: 251 mbarn
+
+c     debug and validity range
+      logical check, info, warn
+      parameter (check=.true., info=.false., warn=.false.)
+
+      integer Ap, At, Zp, Zt, npart, nbar, nmes, ctag
+      integer nsteps,ranseed,event,eos,dectag,uid_cnt
+      integer NHardRes,NSoftRes,NDecRes,NElColl,NBlColl
+      real*8  time,  acttime, bdist, ebeam, bimp,bmin,ecm
+c 7 integer
+
+      common /sys/ npart, nbar, nmes, ctag,nsteps,uid_cnt,
+     +             ranseed,event,Ap,At,Zp,Zt,eos,dectag,
+     +             NHardRes,NSoftRes,NDecRes,NElColl,NBlColl
+      common /rsys/ time,acttime,bdist,bimp,bmin,ebeam,ecm
+
+c     Ap       : projectile mass
+c     Zp       : projectile charge
+c     At       : target mass
+c     Zt       : target charge
+c     npart    : total number of particles
+c     nbar     : number of baryons AND antibaryons
+c     nmes     : number of mesons
+c     ctag     : counter of All interactions (coll. and dec.)
+c     nsteps   : number of timesteps
+c     uid_cnt  : counter for assigning unique particle ID-tags
+c     ranseed  : random number generator seed of event
+c     event    : event counter
+c     eos      : flag for the EoS chosen
+c     dectag   : counter for decays
+c     NHardRes : counter for resonance excitation in BB collisions
+c     NSoftRes : counter for resonance excitation in MB collisions
+c     NElColl  : counter for elastic collisions
+c     NBlColl  : counter for Pauli-blocked collisions
+c     time     : system time at beginning of timestep
+c     acttime  : current system time
+c     bdist    : maximum impact parameter (of event sample)
+c     bimp     : actual impact parameter
+c     bmin     : minimum impact parameter (of event sample)
+c     ebeam    : incident beam energy (lab frame)
+c     ecm      : initial projectile-hadron target-hadron c.m. energy
+
+      logical firstseed
+      common /comseed/firstseed
+
+      logical lsct(nmax),
+     +        logSky, logYuk, logCb, logPau
+      common /logic/ lsct, logSky, logYuk, logCb, logPau
+c 2*nmax*nmax logical
+
+      integer spin(nmax),ncoll(nmax),charge(nmax),strid(nmax),
+     +        ityp(nmax),lstcoll(nmax),iso3(nmax),origin(nmax),uid(nmax)
+c 6*nmax integer
+
+
+
+      real*8 eps, er0, pi, rho0,hqc
+      parameter (eps  = 1.0E-12)  ! small number
+      parameter (er0  = 1.128379167) ! used for error function
+      parameter (pi   = 3.1415926535) ! useful constant
+      parameter (rho0 = 0.16) ! nuclear matter ground state density
+      parameter (hqc  = 0.197327) ! value of $\hbar c$
+
+c IMPORTANT: when you change the version number please change also
+c            the versiontag in blockres.f !
+      integer version, laires
+      parameter ( version = 10035) ! version number
+      parameter ( laires  = 10002) ! additional tag
+
+c MD temporary arrays
+       real*8 r0_t(nmax), rx_t(nmax), ry_t(nmax), rz_t(nmax)
+
+       common/mdprop/ r0_t, rx_t, ry_t, rz_t
+
+      real*8
+     +     gw, sgw, delr, fdel, dt,
+     +     da, db,
+     +     Cb0, Yuk0, Pau0, Sky20, Sky30, gamSky, gamYuk, drPau, dpPau,
+     +     dtimestep
+c 19 real*8
+
+      real*8 cutmax, cutPau, cutCb, cutYuk, cutSky, cutdww
+      common /cuts/ cutmax, cutPau, cutCb, cutYuk, cutSky, cutdww
+
+      real*8 spx(nspl), spPauy(nspl), outPau(nspl),
+     +                spCby(nspl),  outCb(nspl),
+     +                spYuky(nspl), outYuk(nspl),
+     +                spSkyy(nspl), outSky(nspl),
+     +                spdwwy(nspl), outdww(nspl)
+
+      common /spdata/ spx, spPauy, outPau, spCby,  outCb,
+     +                     spYuky, outYuk, spSkyy, outSky,
+     +                     spdwwy, outdww
+
+      real*8
+     +     r0(nmax), rx(nmax), ry(nmax), rz(nmax),
+     +     p0(nmax), px(nmax), py(nmax), pz(nmax),
+     +     airx(nmax), airy(nmax), airz(nmax),
+     +     aipx(nmax), aipy(nmax), aipz(nmax),
+     +     aorx(nmax,4), aory(nmax,4), aorz(nmax,4),
+     +     aopx(nmax,4), aopy(nmax,4), aopz(nmax,4),
+     +     fmass(nmax), rww(nmax),
+     +     dectime(nmax), tform(nmax), xtotfac(nmax)
+
+      common/isys/spin,ncoll,charge,ityp,lstcoll,iso3,origin,strid,
+     +            uid
+      common /coor/ r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime
+      common /frag/ tform, xtotfac
+
+c     spin    : particle spin
+c     ncoll   : particle number of collisions
+c     charge  : particle charge
+c     strid   : ID of last string, the particle was contained in
+c     ityp    : particle ID
+c     lstcoll : tag of last interaction of particle
+c     iso3    : $2 \cdot I_3$ of particle
+c     origin  : ID of last interaction of particle
+c     uid     : unique particle ID tag
+c     r0      : particle time
+c     rx      : $x$ coordinate
+c     ry      : $y$ coordinate
+c     rz      : $z$ coordinate
+c     p0      : particle energy
+c     px      : $p_x$ momentum
+c     py      : $p_y$ momentum
+c     pz      : $p_z$ momentum
+c     fmass   : mass of particle
+c     rww     : ??
+c     dectime : decay time of particle
+c     tform   : formation time of particle
+c     xtotfac : cross section scaling factor during formation time
+c
+      common /aios/ airx, airy, airz, aipx, aipy, aipz,
+     +              aorx, aory, aorz, aopx, aopy, aopz
+
+      common /pots/ Cb0, Yuk0, Pau0, Sky20, Sky30, gamSky,
+     +              gamYuk, drPau, dpPau, gw, sgw, delr, fdel,
+     +              dt,da, db,dtimestep
+
+
+c spectator arrays
+        integer smax
+        parameter(smax=500)  ! maximum number of spectators
+        real*8 r0s(smax), rxs(smax), rys(smax), rzs(smax),
+     +               p0s(smax), pxs(smax), pys(smax), pzs(smax),
+     +               sfmass(smax)
+
+
+        integer sspin(smax), scharge(smax), sityp(smax), siso3(smax),
+     +          suid(smax)
+
+        integer nspec
+
+        common /scoor/ r0s, rxs, rys, rzs, p0s, pxs ,pys, pzs, sfmass
+
+        common /sisys/ sspin, scharge, sityp, siso3, suid
+
+        common /ssys/ nspec
+
+c     sspin    : spectator particle spin
+c     scharge  : spectator particle charge
+c     sityp    : spectator particle ID
+c     siso3    : $2 \cdot I_3$ of spectator particle
+c     r0s      : spectator particle time
+c     rxs      : spectator $x$ coordinate
+c     rys      : spectator $y$ coordinate
+c     rzs      : spectator $z$ coordinate
+c     p0s      : spectator particle energy
+c     pxs      : spectator $p_x$ momentum
+c     pys      : spectator $p_y$ momentum
+c     pzs      : spectator $p_z$ momentum
+c     sfmass   : mass of spectator particle
+
+
+c
+        real*8 p0td(2,nmax),pxtd(2,nmax),pytd(2,nmax),pztd(2,nmax),
+     +         fmasstd(2,nmax)
+        integer ityptd(2,nmax),iso3td(2,nmax)
+        integer itypt(2),uidt(2),origint(2),iso3t(2)
+
+        common /rtdelay/p0td,pxtd,pytd,pztd,fmasstd
+        common /itdelay/ityptd,iso3td
+        common /svinfo/itypt,uidt,origint,iso3t
+
+c     p0td    : energy of parent particles of resonace (DP formalism)
+c     pxtd    : $p_x$ of parent particles of resonace (DP formalism)
+c     pytd    : $p_y$ of parent particles of resonace (DP formalism)
+c     pztd    : $p_z$ of parent particles of resonace (DP formalism)
+c     fmasstd : mass of parent particles of resonace (DP formalism)
+c     ityptd  : ID of parent particles of resonace (DP formalism)
+c     iso3td  : $2\cdot I_3$ of parent particles of resonace (DP formalism)
+
+        real*8 ffermpx(nmax), ffermpy(nmax), ffermpz(nmax)
+        real*8 peq1, peq2
+        common /ffermi/ ffermpx, ffermpy, ffermpz
+        common /peq/ peq1,peq2
+
+c     ffermpx  : fermi momentum in $x$ direction
+c     ffermpy  : fermi momentum in $y$ direction
+c     ffermpz  : fermi momentum in $z$ direction

Processes/UrQMD/comstr.f

+c $Id: comstr.f,v 1.2 1999/01/18 09:56:58 ernst Exp $
+
+      parameter (njspin=8)
+
+      real*8  PJSPNS,PMIX1S(3,njspin),PMIX2S(3,njspin),PBARS,
+     *PARQLS,PARRS
+      real*8  PJSPNC,PMIX1C(3,njspin),PMIX2C(3,njspin),PBARC
+
+
+
+      COMMON/FRGSPA/ PJSPNS,PMIX1S,PMIX2S,PBARS,
+     *PARQLS,PARRS
+      COMMON/FRGCPA/ PJSPNC,PMIX1C,PMIX2C,PBARC
+
+c parm gives the probability for different meson multiplets according
+c to spin degeneracy and average mass ratios
+c spin-parity 0- : 1- : 0+ : 1+ : 2+ = parm(1):parm(2)...:parm(njspin)
+      real*8 parm(njspin)
+
+      common/coparm/parm
+
+      real*8 pi
+      COMMON/CONST/ PI

Processes/UrQMD/comwid.f

+c $Id: comwid.f,v 1.13 2003/06/29 14:26:36 weber Exp $
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     Revision : 1.0
+c
+c     common block for the tabulated branching ratios
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c number of points for the interpolation
+      integer widnsp
+c the branching ratios are interpolated with high precision in the range
+c from mintab to maxtab1 and from maxtab1 to maxtab2 with a smaller precision
+
+c version number of the table
+      integer tabver
+      real*8 mintab, maxtab1, maxtab2
+
+c set the default values
+      parameter (widnsp=120)
+      parameter (mintab=0.1d0)
+      parameter (maxtab1=5.0d0)
+      parameter (maxtab2=50.0d0)
+c increase this parameter, if you make changes, which require a new table
+      parameter (tabver=9)
+
+c tabulated x-values (i.e. sqrt(s) of the collision)
+      real*8 tabx (1:widnsp)
+c tabulated y-values (i.e. the branching ratios) and the second
+c derivatives of the function.
+c full baryon ratio
+      real*8 fbtaby (1:widnsp,minbar:maxbar,1:2)
+c partial baryon ratios
+cdh   real*8 pbtaby (1:widnsp,1:2,minbar:maxbar,0:maxbrs1)
+      real*8 pbtaby (1:widnsp,1:2,minbar:maxbar,0:maxbra)
+c full meson ratio
+      real*8 fmtaby (1:widnsp,minmes:maxmes,1:2)
+c partial meson ratios
+      real*8 pmtaby (1:widnsp,1:2,minmes:maxmes,0:maxbrm)
+
+c Breit-Wigner norms
+c norm of Breit-Wigner with mass dependent widths baryons/mesons
+        real*8 bwbarnorm(minbar:maxbar),bwmesnorm(minmes:maxmes)
+
+c tabulated fppfit()
+c tabulated x-values (i.e. sqrt(s) of the collision)
+      real*8 tabxnd (1:widnsp)
+c 2-resonance channels
+c                      x     deriv ND N*..D*
+      real*8 frrtaby(1:widnsp,1:2,1:2,2:maxdel)
+
+c this flag indicates the progress of tabulating the function
+      integer wtabflg
+c name of file containing the tables
+      character*77 tabname
+
+      common /decaywidth/ tabx,fbtaby,pbtaby,fmtaby,pmtaby,wtabflg
+      common /brwignorm/ bwbarnorm,bwmesnorm
+      common /xsections/ tabxnd,frrtaby
+      common /tabnames/ tabname

Processes/UrQMD/dectim.f

+c $Id: dectim.f,v 1.7 1999/01/18 09:56:59 ernst Exp $
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      real*8 function dectim(ind,iproc)
+c
+c     Revision : 1.0
+C
+cinput ind : ID of particle
+cinput iproc: process ID for resonance creation
+couput dectim: time of decay
+c
+c     This function computes a random choice for the time at which
+c     a resonance will decay and transformes it to the computational
+c     frame.
+c
+c
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+      implicit none
+      include 'coms.f'
+      include 'options.f'
+      include 'colltab.f'
+
+      integer ind,iproc
+      real*8 gg,wid,tau,ranf,fwidth,widit,fbwnorm,mr,massit
+      real*8 tmp,factor
+c
+c first determine width of resonace
+c
+      if(CTOption(1).eq.0.and.CTOption(34).ne.1) then
+c     mass dependent width
+         wid=fwidth(ityp(ind),iso3(ind),fmass(ind))*CTParam(1)
+      else
+c     fixed width
+         wid=widit(ityp(ind))*CTParam(1)
+      end if
+
+C  ...  REST FRAME DECAY TIME
+      if(WID .GT. 1.d-10) then
+         if(CTOption(34).lt.2) then
+c     "normal" life time tau=1/gamma
+            TAU=-(dLOG(1.d0-RANF(0))/WID)
+         else
+c     use Danielewicz delay
+            if(iproc.ne.36.and.iproc.ne.37) then
+c     delay for scattering wave
+               TAU=-(dLOG(1.d0-RANF(0))*
+     &              fbwnorm(fmass(ind),ityp(ind),iso3(ind))*pi/2.d0)
+            else
+c     delay for forward wave
+               if(CTOption(34).eq.2) then
+                  factor=1.d0/CTParam(58)
+               elseif(CTOption(34).eq.3) then
+                  factor=(ctsigtot(actcol)-CTParam(58))/CTParam(58)
+               elseif(CTOption(34).eq.4) then
+                  tmp=dsqrt(2.d0/(wid*3.14d0*
+     &                 fbwnorm(fmass(ind),ityp(ind),iso3(ind))))
+                  factor=1.d0/(CTParam(58)*tmp)
+               else
+                  factor=1.d0
+               endif
+               mr=massit(ityp(ind))
+               TAU=-(dLOG(1.d0-RANF(0))*factor*
+     &              2.d0*pi*fbwnorm(fmass(ind),ityp(ind),iso3(ind))*
+     &              (fmass(ind)-mr)**2/wid**2)
+            endif
+         endif
+      ELSE
+c     stable particle
+         DECTIM=1.d34
+         RETURN
+      END IF
+C  ...  APPLY TIME DILATION
+C  ...  GAMMA FOR THE RESONANCE RESTFRAME <-> COMP. FRAME TRAFO
+      gg=p0(ind)/fmass(ind)
+      DECTIM=TAU*GG*hqc
+
+      RETURN
+      END

Processes/UrQMD/error.f

+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine error (function_name, message, value, level)
+c
+c     Revision : 1.0
+c
+c     output of errors and warnings
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+
+      character function_name*(*)
+      character message*(*)
+      real*8 value
+      integer level
+
+      include 'inputs.f'
+
+      integer errdev
+
+      errdev=6
+
+      if ((level.lt.1).or.(level.gt.3)) then
+         write (errdev,*) '*** Error in Subroutine error:'
+         write (errdev,*) '*** Message: Wrong Errorlevel'
+         write (errdev,*) '*** Related Value: ', level
+      endif
+
+      if (level.eq.1) then
+         write (errdev,*) '*** Warning in Subroutine ',function_name,':'
+      elseif (level.eq.2) then
+         write (errdev,*) '*** Error in Subroutine ',function_name,':'
+      else
+         write (errdev,*) '*** Fatal Error in Subroutine ',
+     $        function_name,':'
+      endif
+      write (errdev,*) '*** Message: ',message
+      write (errdev,*) '*** Related Value: ',value
+
+      if (level.ge.3) then
+         write (errdev,*)
+         write (errdev,*) '*** Program stopped.'
+         stop
+      endif
+
+      return
+      end