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