c c________________________________________________________________________ PROGRAM ELAST C C.A. Bertulani - version 14/jul/99 c________________________________________________________________________ c***** c This program calculates the elastic scattering cross section c with the optical limit to the Glauber model including the c Coulomb correction to the trajectories. C See, Nucl. Phys. A588 (1995) 667 c***** c________________________________________________________________________ C Input should be placed in file ELAST.IN, parameters in ELAST.DIM. C For a description of how to make an input see the end of this program. C A reasonable set of parameters is C N1=2000, N2=2000, N3=6000, N4=7000, N5=400 C C Units are in MeV and fm. Cross section in mb/sr. c***** c Parameters included in file 'elast.dim'. c***** implicit real*8(a-h,o-z) include 'elast.dim' common simq(0:n1),simr(0:n2),simb(0:n3), * qmax,dq,thr,pcm,alfnn,eta,pi,q1,pr0,tr0,nn,nq dimension tfp(0:n1),tft(0:n1),dif(n5),xl(0:n3), * b(0:n3),sigma(n4),fqr(0:n3),fqi(0:n3),delc(n5) open(5,file='elast.in',status='old') open(2,file='elast.out',status='unknown') open(8,file='elast.dat',status='unknown') c iii = 0 c c***** Inputs c read(5,*)ap,zp,at,zt read(5,*)elab,signn,alfnn,af read(5,*)ndenp,ndent read(5,*)qmax,dq,bmin,bmax,db read(5,*)thmin,thmax,dth c c***** OUTPUT #1 c write(2,700)ap,zp,at,zt,elab,signn,alfnn,af,ndenp,ndent, & qmax,dq,bmin,bmax,db,thmin,thmax,dth 700 format(' ***************************************************' &,//, ' checking input data:',//,' ap,zp,at,zt =',4f6.0,/ @ ' elab,signn,alfnn,af =',4(f8.2,1x),/, & ' ndenp,ndent =',2(2x,i4),/, & ' qmax,dq,bmin,bmax,db =',5(f8.2,1x),/, & ' thmin,thmax,dth =',3(f8.2,1x),/, &' ***************************************************',//) c qmin=0.01 nth=(thmax-thmin)/dth+1.0001 c constants pi=acos(-1.) amu=938. hbarc=197.3 c momenta and energies with relativistic kinematics c for the nuclei redm=ap*at/(ap+at) apmass=ap*amu atmass=at*amu eta=zp*zt*1.44/hbarc*sqrt(apmass*0.5/elab) plab=sqrt(elab*(elab+2.*apmass)) beta=plab/(elab+apmass) gamma=1./sqrt((1.-beta)*(1.+beta)) ecm=sqrt(atmass**2+apmass**2+2.*(elab+apmass)*atmass) pcm=plab*atmass/ecm pcm=pcm/hbarc c for the nucleons plab_nn=plab/ap ecm_nn=sqrt(2.*amu**2+2.*(elab/ap+amu)*amu) pcm_nn=plab_nn*amu/ecm_nn pcm_nn=pcm_nn/hbarc C mesh in impact parameter and momentum transfer nq=qmax/dq+1.0001 nb=(bmax-bmin)/db+1.0001 nbmin=bmin/db+0.00001 nbmax=bmax/db+0.00001 C factors for nucleon-nucleon amplitudes fnnr=alfnn*signn*pcm_nn/(4.*pi) fnni=signn*pcm_nn/(4.*pi) C variables used in Simpson's integration C ** for space integration d=1./3.0 d2=2.0*d d4=4.0*d simr(0)=d*0.1 nsimp=1999 do1i=1,nsimp,2 simr(i)=d4*0.1 1 simr(i+1)=d2*0.1 simr(2000)=d*0.1 C ** for momentum integration simq(0)=d*dq do 11 iq=1,nq-1,2 simq(iq)=d4*dq 11 simq(iq+1)=d2*dq simq(nq)=d*dq C ** for impact parameter integration simb(nbmin)=d*db simb(nbmax)=d*db do 12 ib=nbmin+1,nbmax-1,2 simb(ib)=d4*db 12 simb(ib+1)=d2*db c c** projectile data and Fourier transform of densities if(ndenp.eq.2)go to 2 read(5,*)dum1,a0p,dum2 r0p=ap/(pi**1.5*a0p**3) call fdgau(r0p,a0p,tfp) go to 3 2 continue read(5,*)pc,pw,pz bp=20. call fnor(bp,pc,pw,pz,pnor) pr0=ap/pnor q1=qmin-dq do 4 iq=0,nq q1=q1+dq call fdws(bp,pc,pw,pz,pr0,ftri) tfp(iq)=ftri 4 continue 3 continue c c** projectile data and Fourier transform of densities if(ndent.eq.2)go to 5 read(5,*)dum1,a0t,dum2 r0t=at/(pi**1.5*a0t**3) call fdgau(r0t,a0t,tft) go to 6 5 continue read(5,*)tc,tw,tz bt=20. call fnor(bt,tc,tw,tz,tnor) tr0=at/tnor q1=qmin-dq do 7 iq=0,nq q1=q1+dq call fdws(bt,tc,tw,tz,tr0,ftri) tft(iq)=ftri 7 continue 6 continue c c** Coulomb phase shift c delc(1)=atan(eta) do 49 l=2,n5 delc(l)=delc(l-1)+atan(eta/float(l)) fase=eta*log(l+0.5) dif(l)=delc(l)-fase dife=dif(l)-dif(l-1) if(abs(dife).lt.1.0e-04)go to 490 49 continue 490 fase=dif(l) print*,fase c c**** elastic phase shift c bs=bmin-db do 51 jb=nbmin,nbmax c c correction for Coulomb deflection bs=bs+db xl(jb)=eta+dsqrt(eta**2.+(pcm*bs)**2.) b(jb)=xl(jb)/pcm c fesum=0. qs=qmin-dq do 50 jq=0,nq qs=qs+dq qb=qs*b(jb) fj0=bessj0(qb) ef=exp(-0.5*af*qs**2) felast=tfp(jq)*tft(jq)*qs*fj0*ef/pcm_nn*simq(jq) fesum=fesum+felast 50 continue ans=fesum ansr=ans*fnnr ansi=ans*fnni eansi=exp(-ansi) c c*** these are the glauber phases: c sr = eansi*cos(ansr) si = eansi*sin(ansr) sigc=eta*log(pcm*bs)+fase C and their exponentials fqr(jb)=-sin(2*sigc)+eansi*sin(ansr+2*sigc) fqi(jb)=cos(2*sigc)-eansi*cos(ansr+2*sigc) 51 continue c c*** cross section c th=thmin-dth do 32 jth=1,nth th=th+dth thr=th*pi/180. fir=0. fii=0. bs=bmin-db do 31 jb=nbmin,nbmax bs=bs+db arg=bs*2.*pcm*sin(thr/2.) fj0t=bessj0(arg) fj0b=fj0t*bs fintqr=fqr(jb)*fj0b*simb(jb) fintqi=fqi(jb)*fj0b*simb(jb) fir=fir+fintqr fii=fii+fintqi 31 continue anbr=fir anbi=fii C Coulomb amplitudes a0=-eta/pcm fcr=a0/2./(sin(thr/2.)**2)*cos(eta*log(sin(thr/2.)**2)) fci=-a0/2./(sin(thr/2.)**2)*sin(eta*log(sin(thr/2.)**2)) C Rutherford cross section sigmam=(eta/(2*pcm))**2/(sin(thr/2))**4 ampr=fcr+pcm*anbr ampi=fci+pcm*anbi sigma(jth)=(ampr**2+ampi**2) sipi=sigma(jth)/sigmam sigma(jth) = log(sigma(jth)) c c*** OUTPUT #2 c this is the Cross Section on a log scale: c write(2,*)th,sigma(jth) csig to coul ratio write(2,*)th,sipi write(8,*)th,sipi 32 continue stop end c c________________________________________________________________________ c________________________________________________________________________ c c subroutine fdgau(r0,a0,tf) c c**** This routine calculates the Fourier transform of the density c implicit real*8(a-h,o-z) include 'elast.dim' common simq(0:n1),simr(0:n2),simb(0:n3), * qmax,dq,thr,pcm,alfnn,eta,pi,q1,pr0,tr0,nn,nq dimension tf(n4) qmin=0.01 q=qmin-dq do 1 i=1,nq q=q+dq tf(i)=(pi*a0**2)**1.5*r0*exp((-1)*(q*a0)**2/4) 1 continue return end c c________________________________________________________________________ c________________________________________________________________________ c subroutine fnor(bx,c,w,z,xnor) C c*** subroutine for the calculation of the normalization of the **** density in the case of Fermi distrib. Normalizes to the mass. C implicit real*8(a-h,o-z) include 'elast.dim' common simq(0:n1),simr(0:n2),simb(0:n3), * qmax,dq,thr,pcm,alfnn,eta,pi,q1,pr0,tr0,nn,nq xnor=0. ib=bx/0.01 do 1 i=0,ib x=i*.1 temp=(x-c)/z if(temp.gt.20.)temp=20. den=(1+w*(x/c)**2)/(1.+exp(temp)) xfnor=den*pi*4.*x**2.*simr(i) xnor=xnor+xfnor 1 continue return end c c________________________________________________________________________ c________________________________________________________________________ c subroutine fdws(bx,c,w,z,r0,fd) C c*** subroutine for the calculation of the density in the momentum c*** space using a Fermi shape for the density C implicit real*8(a-h,o-z) include 'elast.dim' common simq(0:n1),simr(0:n2),simb(0:n3), * qmax,dq,thr,pcm,alfnn,eta,pi,q1,pr0,tr0,nn,nq fd=0. ib=bx/0.1 do 1 i=0,ib x=i*0.1 if(abs(q1).le.1.e-20)go to 25 fac=4.*pi*r0*x*sin(x*q1)/q1 go to 26 25 continue fac=4.*pi*r0*x**2 26 continue temp=(x-c)/z if(temp.gt.20.)temp=20. fac1=exp(temp) sum=fac/(1.+fac1)*simr(i) if(abs(w).lt.1.e-20)go to 3 sum=sum+fac*w*(x/c)**2/(1.+fac1)*simr(i) 3 continue fd=fd+sum 1 continue return end c c________________________________________________________________________ c________________________________________________________________________ c function bessj0(x) c*** returns the Bessel function of order 0 implicit real*8(a-h,o-z) data p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, * -.2073370639d-5,.2093887211d-6/,q1,q2,q3,q4,q5/-.1562499995d-1, * .1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/ data r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0, * 651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/, * s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0, * 9494680.718d0,59272.64853d0,267.8532712d0,1.d0/ if(abs(x).lt.8.)then y=x**2 bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))) * /(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))) else ax=abs(x) z=8./ax y=z**2 xx=ax-.785398164 bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+ * y*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) endif return end C C**************************************************************************** C* C************ Explaining the inputs of program ELAST *************** C ***** The actual input program is denoted ELAST.IN ******* C ***** Integration parameters are in file ELAST.DIM ******* C ***** Units are in MeV and fm ***** C------------------------------------------------------------------------------ C 1st row: AP,ZP,AT,ZT c - AT and AP = target and projectile masses c - ZT and ZP = target and projectile number of protons c------------- C 2nd row: ELAB,SIGNN,ALFNN,AF C **** SIGNN, ALFNN, and AF are parameters for the elementary C scattering amplitude c - ELAB = incident energy in the laboratory c - SIGNN = nucleon-nucleon cross section in fm^2 c - ALPFNN = ratio Real-to-Imaginary part of the nucleon-nucleon c scattering amplitude c f_nn = (k/4/pi) . signn . (i+alfnn) .exp(-q^2.af/2) c-------------- C 3rd row: NDENP,NDENT C - NDEN* = 1 (2) = gaussian (Woods-Saxon) density c-------------- C 4th row: QMAX,DQ,BMIN,BMAX,DB C 5th row: THMIN,THMAX,DTH c ***** These are integration data c - QMAX = approx. pi/R, where R is the nuclear size c - DQ = integration step in q c - BMIN = minimum impact parameter for the integral c which calculates the eikonal phase c - BMAX = corresponding maximum value of the b c - DB = integration step in b C - THMIN = minimum scattering angle (degrees) C - THMAX = maximum scattering angle C - DTH = step in the scattering angle c-------------- C 6th row: DUM1,A0P,DUM2 C 6th row: PC,PW,PZ c **** parameters for the density of the projectile c*** (Gaussian, nden*=1): (first choice for 6th row) c rho(r)=rho(0) . exp(- r^2 / a0*^2), with rho(0)=a*/(pi^1.5 . a0*^3) c In this case set 1st and last entry of 6th row (dum1 and dum2) to c anything. c*** (Woods-Saxon, nden*=2): (second choice for 6th row) c rho(r)=rho(0) . (1 + *w . r^2 / *c^2) / ( 1+exp(-(r-*c)/*z) ) c In the case of the WS the program calculates rho(0) via the suboutine c fnorn. c-------------- C 7th row: DUM1,A0T,DUM2 C 7th row: TC,TW,TZ c **** parameters for the density of the target C same rules as above applies c***************