c **** smatx **** c calculates s-matrix and phase shifts c C.A. Bertulani, August, 1990 c implicit real*8(a-h,o-z) common/snbar/ snbar,pshift common/der/am2,aux1 common/inte/dr,ri,nax common/et/eta,xk common/pot/v0r,r0r,ar,v0i,r0i,ai,am,a13,ecm,all,l common/cou/z1,z2,a1,a2,rc,rc2,aux6 common/aux/aux4,aux5 c include 'crossec.dim' complex*16 psi(nmesh),sn(0:ldim),fcou,fn,ci dimension cph(0:ldim),pl(0:ldim) c open(1,file='crossec.in',status='old') open(10,file='crossec.out',status='unknown') open(11,file='smat.dat',status='unknown') open(12,file='angdist.dat',status='unknown') c read(1,*)v0r,r0r,ar,v0i,r0i,ai,rc0 c only coulomb? if yes, inuc=0 inuc=1 checkn=abs(v0r)+abs(v0i) if(check.lt.1d-8)inuc=0 c read(1,*)z1,a1,z2,a2,ecm,lmax read(1,*)dr,nax read(1,*)tmin,tmax,dt c ident=0 check=abs(z1-z2)+abs(a1-a2) if(check.lt.1.) ident=1 nt=int((tmax-tmin)/dt)+1 ci=cmplx(0.d0,1.d0) am0=a1*a2/(a1+a2) a13=a1**.333333+a2**.333333 rc=rc0*a13 rc2=rc*2. am=931.5*am0 am2=2.*am h2=197.33**2 aux4=r0r*a13 aux5=r0i*a13 aux6=1.44*z1*z2 pi=acos(-1.d0) c aux1=am2*ecm/h2 xk=dsqrt(aux1) aux1=aux1/ecm eta=1.44*z1*z2*am/xk/h2 aaa=eta/xk ruth90=aaa**2 c write(10,100)z1,a1,z2,a2,ecm,ldim,eta,xk,10.*ruth90 c c-------------- angular momentum mesh -------------------------------- c do l=0,lmax all=dble(l*(l+1)) sn(l)=cmplx(1.d0,0.d0) if(inuc.eq.1)then call wave(psi) sn(l)=snbar*cmplx(cos(2.*pshift),sin(2.*pshift)) end if end do c c c-------------- Angular distribution ---------------------------------- c c call coul(lmax,eta,cph) c do l=0,lmax cps=cph(l)/pi snr=real(sn(l)) sni=dimag(sn(l)) psn=atan2(sni,snr)/pi/2. write(10,200)l,abs(sn(l)),psn,cps write(11,200)l,abs(sn(l)),psn,cps write(* ,200)l,abs(sn(l)),psn,cps end do c write(10,300) do it=1,nt tdeg=tmin+dble(it-1)*dt t=tdeg*pi/180. cost=cos(t) call legpol(lmax,cost,pl) fn=cmplx(0.d0,0.d0) if(inuc.eq.1)then do l=0,ldim delta=1.d0 if(ident.eq.1)delta=dble(1.d0+(-1.d0)**l) fn=fn+delta*dble(2*l+1)*exp(2.*ci*cph(l))* & (sn(l)-cmplx(1.d0,0.d0))*pl(l)/(2.*ci*xk) end do end if s2=sin(t/2.)**2 au=-eta*log(s2)+pi+2.*cph(0) fcou=eta*exp(ci*au)/(2.*xk*s2) if(ident.eq.1) then c2=cos(t/2.)**2 au=-eta*log(c2)+pi+2.*cph(0) fcou=fcou+eta*exp(ci*au)/(2.*xk*c2) end if sigma=abs(fn+fcou)**2/ruth90 c multiplication by 10 transforms the cross section in mb c sigruth=(aaa/(2.*s2))**2/ruth90 if(ident.eq.1) sigruth=sigruth+10.*(aaa/(2.*c2))**2/ruth90 ratio=sigma/sigruth write(10,400)tdeg,sigma,sigruth,ratio write(12,400)tdeg,sigma,sigruth,ratio end do c 100 format('SYSTEM DATA:',//,' Z1=',f3.0,2x,'Z2=',f3.0,2x,'A1=',f3.0, & 2x,'A2=',f3.0,//,' Ecm(MeV)=',g14.4,'lmax=',i3,//,' eta=', & g12.4,2x,'k(fm)=',g12.4,5x,'SIG_ruth(90)=',g12.4,'mb', &///,3x,'l',5x,'|S_l|',5x,'n-pshift',5x,'C-pshift',//) 200 format(2x,i3,3(2x,g12.4)) 300 format(///,' Angular Distribution',//, & ' Theta(deg)',3x,'sigma/a**2',5x,'sigruth/a**a',5x, & 'Ratio',//) 400 format(2x,f8.2,4x,g12.4,2x,g12.4,2x,g12.4) c stop end c c function vnr(r) implicit real*8(a-h,o-z) common/aux/aux4,aux5 common/pot/v0r,r0r,ar,v0i,r0i,ai,am,a13,e,all,l cr=dexp((r-aux4)/ar) vnr=v0r/(1.+cr) vnr=v0r*exp(-(r/aux4)**2) return end c function vni(r) implicit real*8(a-h,o-z) common/aux/aux4,aux5 common/pot/v0r,r0r,ar,v0i,r0i,ai,am,a13,e,all,l ci=dexp((r-aux5)/ai) vni=v0i/(1.+ci) return end c function vcou(r) implicit real*8(a-h,o-z) common/cou/z1,z2,a1,a2,rc,rc2,aux6 if(r.gt.rc)go to 1 vcou=aux6*(3.-(r/rc)**2)/rc2 return 1 vcou=aux6/r return end c function deriv(r) implicit real*8(a-h,o-z) complex *16 deriv common/der/am2,aux1 common/cou/z1,z2,a1,a2,rc,rc2,aux6 common/pot/v0r,r0r,ar,v0i,r0i,ai,am,a13,e,all,l derivr=all/(r*r)-aux1*(e-vcou(r)-vnr(r)) derivi=aux1*vni(r) deriv = dcmplx(derivr,derivi) return end c c subroutine wave(psi) implicit real*8(a-h,o-z) common/snbar/snbar,pshift common/cou/z1,z2,a1,a2,rc,rc2,aux6 common/inte/dr,ri,nax common/et/eta,xk common/pot/v0r,r0r,ar,v0i,r0i,ai,am,a13,e,all,l c include 'crossec.dim' dimension f(ldim),fp(ldim),g(ldim),gp(ldim),cnn(nmesh) complex*16 hp,hm,dhp,dhm,dlg,anorm,snuc complex*16 psi(nmesh) complex*16 yy(2) c dr1=dr/10. n=nax+1 psi(1)=cmplx(0.d0,0.d0) yy(1)=cmplx(0.d0,0.d0) yy(2)=cmplx(1.d0,0.d0) r=1.d-10 cnn(1)=0.d0 do 1000 i=2,n cnn(i)=cnn(i-1) call rk4(yy,2,r,dr) chn=abs(yy(1)) if(chn.gt.1.d30)then cnn(i)=cnn(i)+log(chn) yy(1)=yy(1)/chn yy(2)=yy(2)/chn end if psi(i)=yy(1) 1000 continue rmax=dr*dble(nax) lp=l+1 rom=rmax*xk accur=1.d0 step=100.d0 call rcwfn(rom,eta,l,l,f,fp,g,gp,accur,step) devi=r-rmax if (dabs(devi).gt.dr1) write(6,707)r,rmax,devi 707 format(/,' deu zebra!! rmatch (=',f8.4,' fm) diferente 1 de rmax (=',f8.4,' fm) !! delta r=',g10.2//) hp=dcmplx(g(lp),f(lp)) hm=dconjg(hp) dhp=dcmplx(gp(lp),fp(lp)) dhm=dconjg(dhp) dlg=yy(2)/yy(1) snuc=(xk*dhm-dlg*hm)/(xk*dhp-dlg*hp) anorm=.5*(hm-snuc*hp)*dcmplx(cexp(cmplx(0.,-1.57)))/yy(1) do 2000 i=2,n reno=exp(cnn(n)-cnn(i)) psi(i)=psi(i)*anorm/reno r=dr*dble(i-1) 2000 continue snbar=abs(snuc) pshift=datan2(dimag(snuc),dble(real(snuc)))/2. return end c subroutine derivs(x,y,dydx) implicit real*8(a-h,o-z) complex*16 y(2),dydx(2),deriv dydx(1)=y(2) dydx(2)=y(1)*deriv(x) return end c c subroutine rk4(y,n,x,h) parameter (nmax=10) implicit real*8(a-h,o-z) complex*16 y,dydx,yt,dyt,dym dimension y(n),dydx(nmax),yt(nmax),dyt(nmax),dym(nmax) hh=h*.5 h6=h/6. xh=x+hh c call derivs(x,y,dydx) do 11 i=1,n yt(i)=y(i)+hh*dydx(i) 11 continue call derivs(xh,yt,dyt) do 12 i=1,n yt(i)=y(i)+hh*dyt(i) 12 continue call derivs(xh,yt,dym) do 13 i=1,n yt(i)=y(i)+h*dym(i) dym(i)=dyt(i)+dym(i) 13 continue call derivs(x+h,yt,dyt) do 14 i=1,n y(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i)) 14 continue x=x+h return end c c ---------------------------------------------------------------------- subroutine rcwfn(rho,eta,minl,maxl,fc,fcp,gc,gcp,accur,step) c ---------------------------------------------------------------------- implicit real*8 (a-h,o-z) c include 'crossec.dim' c real*8 k,k1,k2,k3,k4,m1,m2,m3,m4 dimension fc(ldim),fcp(ldim),gc(ldim),gcp(ldim) c*** coulomb wavefunctions calculated at r=rho by the c*** continued fraction method of steed minl,maxl are actual l-values c*** see barnett feng steed and goldfarb computer physics commun 1974 pace=step acc=accur if(pace.lt.100.0)pace=100.0 if(acc.lt.1.0e-15.or.acc.gt.1.0e-6)acc=1.0e-6 r=rho ktr=1 lmax=maxl lmin1=minl+1 xll1=dfloat(minl*lmin1) eta2=eta*eta turn=eta+dsqrt(eta2+xll1) if(r.lt.turn.and.dabs(eta).ge.1.0e-6)ktr=-1 ktrp=ktr go to 2 1 r=turn tf=f tfp=fp lmax=minl ktrp=1 2 etar=eta*r rho2=r*r pl=dfloat(lmax+1) pmx=pl+0.5 c*** continued fraction for fp(maxl)/f(maxl) xl is f xlprime is fp fp=eta/pl+pl/r dk=etar*2.0 del=0.0 d=0.0 f=1.0 k=(pl*pl-pl+etar)*(2.0*pl-1.0) if(pl*pl+pl+etar.ne.0.0)go to 3 r=r+1.0e-6 go to 2 3 h=(pl*pl+eta2)*(1.0-pl*pl)*rho2 k=k+dk+pl*pl*6.0 d=1.0/(d*h+k) del=del*(d*k-1.0) if(pl.lt.pmx)del=-r*(pl*pl+eta2)*(pl+1.0)*d/pl pl=pl+1.0 fp=fp+del if(d.lt.0.0)f=-f if(pl.gt.20000.)go to 11 if(dabs(del/fp).ge.acc)go to 3 fp=f*fp if(lmax.eq.minl)go to 5 fc(lmax+1)=f fcp(lmax+1)=fp c*** downward recursion to minl for f and fp,arrays gc,gcp are storage l=lmax do 4 lp=lmin1,lmax pl=dfloat(l) gc(l+1)=eta/pl+pl/r gcp(l+1)=dsqrt(eta2+pl*pl)/pl fc(l)=(gc(l+1)*fc(l+1)+fcp(l+1))/gcp(l+1) fcp(l)=gc(l+1)*fc(l)-gcp(l+1)*fc(l+1) 4 l=l-1 f=fc(lmin1) fp=fcp(lmin1) 5 if(ktrp.eq.-1)go to 1 c*** repeat for r=turn if rho lt turn c*** now obtain p+i.q for minl from continued fraction (32) c*** real arithmetic to facilitate conversion to ibm using real*8 p=0.0 q=r-eta pl=0.0 ar=-(eta2+xll1) ai=eta br=2.0*q bi=2.0 wi=2.0*eta dr=br/(br*br+bi*bi) di=-bi/(br*br+bi*bi) dp=-(ar*di+ai*dr) dq=(ar*dr-ai*di) 6 p=p+dp q=q+dq pl=pl+2.0 ar=ar+pl ai=ai+wi bi=bi+2.0 d=ar*dr-ai*di+br di=ai*dr+ar*di+bi t=1.0/(d*d+di*di) dr=t*d di=-t*di h=br*dr-bi*di-1.0 k=bi*dr+br*di t=dp*h-dq*k dq=dp*k+dq*h dp=t if(pl.gt.46000.)go to 11 if(dabs(dp)+dabs(dq).ge.(dabs(p)+dabs(q))*acc)go to 6 p=p/r q=q/r c*** solve for fp,g,gp and normalise f at l=minl g=(fp-f*p)/q gp =p*g-q*f w=1.0/dsqrt(fp*g-f*gp) g=w*g gp=w*gp if(ktr.eq.1)go to 8 f=tf fp=tfp lmax=maxl c*** runge-kutta integration of g(minl) and gp(minl) inwards from turn c*** see fox and mayers 1968 pg 202 if(rho.lt.0.2*turn)pace=999.0 r3=1.0/3.0 h=(rho-turn)/(pace+1.0) h2=0.5*h i2=(pace+0.001) etah=eta*h h2ll=h2*xll1 s=(etah+h2ll/r)/r-h2 7 rh2=r+h2 t=(etah+h2ll/rh2)/rh2-h2 k1=h2*gp m1=s*g k2=h2*(gp+m1) m2=t*(g+k1) k3=h*(gp+m2) m3=t*(g+k2) m3=m3+m3 k4=h2*(gp+m3) rh=r+h s=(etah+h2ll/rh)/rh-h2 m4=s*(g+k3) g=g+(k1+k2+k2+k3+k4)*r3 gp=gp+(m1+m2+m2+m3+m4)*r3 r=rh i2=i2-1 if(i2.ge.0)go to 7 w=1.0/(fp*g-f*gp) c*** upward recursion from gc(minl) and gcp(minl),stored values are r,s c*** renormalise fc,fcp for each l-value 8 gc(lmin1)=g gcp(lmin1)=gp if(lmax.eq.minl)go to10 do 9 l=lmin1,lmax t=gc(l+1) gc(l+1)=(gc(l)*gc(l+1)-gcp(l))/gcp(l+1) gcp(l+1)=gc(l)*gcp(l+1)-gc(l+1)*t fc(l+1)=w*fc(l+1) 9 fcp(l+1)=w*fcp(l+1) fc(lmin1)=fc(lmin1)*w fcp(lmin1)=fcp(lmin1)*w return 10 fc(lmin1)=w*f fcp(lmin1)=w*fp return 11 w=0.0 g=0.0 gp=0.0 go to 8 end c c c c===================================================================== subroutine legpol(maxl,cosbet,pl) c===================================================================== c c Calculates Legenrde Polynomials PL(L) = P_L (COSBET) c c--------------------------------------------------------------------- c include 'crossec.dim' implicit real*8 (a-h,o-z) dimension pl(0:ldim) c pl(0) = 1.d0 pl(1) = cosbet do 1 l = 2,maxl 1 pl(l) = ((2*l-1)*cosbet*pl(l-1)-(l-1)*pl(l-2))/l c return c-end-legpol end c c====================================================================== subroutine coul(lmax,eta,cph) c---------------------------------------------------------------------- c calculates Coulomb phase-shifts c---------------------------------------------------------------------- implicit real*8(a-h,o-z) c include 'crossec.dim' dimension cph(0:ldim) c sto=16.d0+eta*eta cph(0)=-eta+(eta/2.d0)*dlog(sto)+3.5d0*datan(eta/4.d0)-( 1 datan(eta)+datan(eta/2.d0)+datan(eta/3.d0))-(eta/(12.d0* 2 sto))*(1.d0+(1.d0/30.d0)*(eta**2-48.d0)/sto**2+(1.d0/105.d0) 3 *(eta**4-160.d0*eta**2+1280.d0)/sto**4) c do 1 ii=1,lmax fi=dble(ii) cph(ii)=cph(ii-1)+datan(eta/fi) 1 continue return end c