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

