*module cstes MODULE cstes c Physical constants: hc=hbar*c; e2=e^2/4 pi epsilon_0;hm=hbar^2/2m_N; ci=i * real*8, parameter :: pi=3.141592653589793238d0 real*8, parameter :: pi=acos(-1d0) real*8 :: hc,e2,hm complex*16, parameter :: ci=(0,1) END MODULE *module reseaufo MODULE reseaufo c Wave function meshes integer :: Nr,Nt,Np real*8, parameter :: eps=1e-8 real*8 :: rm,hr END MODULE *module reseaueik MODULE reseaueik c Meshes for eikonal phase calculation integer :: Nb,NZ,Ns real*8 :: bmax,hbpol,v0,hZ,hs real*8, allocatable, dimension (:) :: xb END MODULE *module collision MODULE collision c Data about the collision: spin of the fragment, mass and charge of the c core(c), fragment(f), and target (T), reduced c-f mass,... integer :: iS,iSc,kZf,kZc,kZT real*8 :: e2ct,e2ft,Ac,Af,Apr,cc,cf,qc,qf,ATa,deuxm real*8 :: EAP,EKPT,etaPT END MODULE *module potentials MODULE potentials c Parameters of the two-body potentials: Vcf,VcT,VfT integer :: ntypo,Njpo,npotN,lamax,lamin,jelj integer, allocatable, dimension (:) :: Jpo,lpo real*8, allocatable, dimension (:) :: Vp,rp,ap,Vls,rls,als,rC real*8, allocatable, dimension (:,:) :: Vcf real*8 :: VcT,WcT,WDcT,rRcT,rIcT,rDcT,aRcT,aIcT,aDcT,rCcT real*8 :: VfT,WfT,WDfT,rRfT,rIfT,rDfT,aRfT,aIfT,aDfT,rCfT END MODULE *module boundstates MODULE boundstates c Physical bound states of the projectile at the mesh points integer :: Nel integer, allocatable, dimension (:) :: NRel,Jel,lel real*8, allocatable, dimension (:) :: Eel real*8, allocatable, dimension (:,:) :: foel END MODULE *module eli MODULE eli c Bound states of the c-f potential to be removed by SuSy integer :: Neli integer, allocatable, dimension (:) :: NReli,Jeli,leli real*8, allocatable, dimension (:) :: Eeli END MODULE *module phaseik MODULE phaseik c Eikonal phase, and parameters for their calculation: c lec determines the type of approximation c relcorr = 1 imply the use of gamma in the first order E1 phase c gaminv is the inverse of gamma => gaminv=sqrt(1-(v/c)^2) integer :: lec,Nl,relcorr real*8 :: gaminv complex*16, allocatable, dimension (:,:,:,:):: feik END MODULE *module Sbu MODULE Sbu c Breakup amplitude and angular amplitude c Dimension of these amplitudes integer :: NlJM,Nsed,Nq,NE real*8, allocatable, dimension (:) :: theta,Enrj complex*16, allocatable, dimension(:,:,:,:) :: SklJM,fbu END MODULE *pgm principal Chaconne program Chaconne c Pgm f90 c Programme that computes reactions involving one-nucleon halo nuclei c (2-body projectile) within the Coulomb-Corrected Eikonal approximation. c It can compute the elastic and inelastic scattering cross sections c (angular distributions) and the breakup cross sections c (angular distribution over a range of continuum energies, c energy distribution, parallel-momentum distribution) c Coulomb breakup reactions including both c The optical potentials between projectile and target c include both nuclear and Coulomb P-T interactions. c The Coulomb eikonal phase is renormalised as proposed by c Abu-Ibrahim and Suzuki [Prog. Theor. Phys. 112, 1013 (2004)], c following an idea of Margueron, Bonnaccorso, and Brink c [NPA 720, 337 (2003)]. c The programme is based on a subroutine written by D. Baye 'multeik'. USE cstes USE reseaufo USE reseaueik USE collision USE potentials USE boundstates USE eli USE phaseik USE Sbu implicit real*8 (a-h,o-z) parameter (lnom=80,nbes=2,epsrel=5d-2) character (len=2) typOL,lnb,c2l,Rel character (lnom) fich,fichres,fichs0,fichsk,fichcT,fichfT real*8 :: fp(30),yp(30),Pl(30),Ylm(466) complex*16 rint,csljp,csljm,csljpE,csljmE,cgamma,cg,cres,ctmp, a feik1,feik2 integer, allocatable, dimension(:):: Nob real*8, allocatable, dimension(:):: bmin,hb real*8, allocatable, dimension(:):: fBessJ,fBessY real*8, allocatable, dimension(:):: sruth,fibb real*8, allocatable, dimension(:):: S,sdE,sedbdl,selast real*8, allocatable, dimension(:,:):: sinel real*8, allocatable, dimension(:):: Si02 real*8, allocatable, dimension(:):: sedp real*8, allocatable, dimension (:,:,:):: foli real*8, allocatable, dimension(:,:):: ADlj real*8, allocatable, dimension(:,:,:,:):: ADfull complex*16, allocatable, dimension(:):: fcoul,ftild complex*16, allocatable, dimension(:,:):: phasOL complex*16, allocatable, dimension(:):: S20,Pstr complex*16, allocatable, dimension(:,:,:):: S0M,S0iM c Input reading c Physical constants: hc=hbar*c; e2=e^2/4 pi epsilon_0;hm=hbar^2/2m_N read(*,*)hc,e2,hm amn=hc**2/(2*hm) * write(*,*)' hbar^2/2m_N',hm,' m_N = ',amn c Mass, charge, and spin of core(c) and the fragment(f) read(*,*)Ac,Af,kZc,kZf,iSc,iS Apr=Af+Ac deuxm=(Af*Ac)/Apr/hm !=2\mu/\hbar^2 * write(*,*)' hbar^2/2mu',1/deuxm c Spin of f and c (up to now Sc is set equal to 0) nAc=Ac+0.5d0 nAf=Af+0.5d0 lzc=2-mod(nAc,2) if(iS==0)then c fragment spin is neglected lzf=2 else lzf=2-mod(nAf,2) endif iSc=iSc*lzc iS=iS*lzf c c-f potential read(*,*)ntypo,Njpo c ntypo= c 1,11: different potentials for Njpo-1 partial waves +1 for all the others c 2,12: different potentials for l-even and l-odd partial waves c 1,2: Woods Saxon potentials c 11,12: Gaussian potentials if(mod(ntypo,10)==2)then Njpo=2 elseif(mod(ntypo,10)/=1)then write(*,*)'Unprogrammed potential' goto 9999 endif allocate(Jpo(Njpo),lpo(Njpo)) allocate(Vp(Njpo),rp(Njpo),ap(Njpo),rC(Njpo)) allocate(Vls(Njpo),rls(Njpo),als(Njpo)) If(mod(ntypo,10)==1)then do ipo=1,Njpo-1 read(*,*)Jpo(ipo),lpo(ipo),Vp(ipo),rp(ipo),ap(ipo), a Vls(ipo),rls(ipo),als(ipo),rC(ipo) enddo Jpo=Jpo*lzf read(*,*)Vp(Njpo),rp(Njpo),ap(Njpo), a Vls(Njpo),rls(Njpo),als(Njpo),rC(Njpo) ElseIf(mod(ntypo,10)==2)then read(*,*)Vp(1),rp(1),ap(1),Vls(1),rls(1),als(1),rC(1)!1=odd l read(*,*)Vp(2),rp(2),ap(2),Vls(2),rls(2),als(2),rC(2)!2=even l EndIf c Bound states (etats lie: el) read(*,*)Nel allocate(NRel(Nel),Jel(Nel),lel(Nel),Eel(Nel)) read(*,*)(NRel(iel),Jel(iel),lel(iel),Eel(iel),iel=1,Nel) iJ0=Jel(1) l0=lel(1) Nef=iJ0+1 Nsed=(Nef+1)/2 Nelm=Nef do iel=2,Nel Nelm=Nelm+Jel(iel)+1 enddo ********************************************************************** c Forbidden states to be removed by SuSy * read(*,*)Neli Neli=0 * allocate(NReli(Neli),Jeli(Neli),leli(Neli),Eeli(Neli)) * read(*,*)(NReli(ieli),Jeli(ieli),leli(ieli),Eeli(ieli), * a ieli=1,Neli) * Jeli=Jeli*lzf ********************************************************************** c Collision data (target mass and charge, P energy/nucleon, trajectory type) read(*,*)ATa,kZT,EAP !,Rel c Rel= c RV: relativistic velocity c NR: non-relativistic velocity Rel='RV' select case (Rel) case('RV','Rv','rV','rv') !relativistic velocity Rel='RV' v0=1+2*EAP*hm/hc**2 v0=hc*sqrt(1-1/v0**2) !relativistic relative velocity*hbar case('NR','Nr','nR','nr') Rel='NR' v0=2*sqrt(EAP*hm) !non-relativistic velocity*hbar end select EKPT=abs(Apr*ATa/(Apr+ATa)/2/hm*v0) !incident wave vector e2ct=e2*kZc*kZT e2ft=e2*kZf*kZT cc=-Af/Apr cf=Ac/Apr etaPT=(e2ct+e2ft)/abs(v0) qc=2*e2ct*cc/v0 qf=2*e2ft*cf/v0 c Relativistic corrections (only if realcorr = 1) * read(*,*) relcorr relcorr=0 if(relcorr==1)then gaminv = 1.0/(1+2*EAP*hm/hc**2) !sqrt(1.0-(v0/hc)**2) ! else gaminv = 1.0 endif c Projectile-target potentials: npotN=V_PT type; bmaxN=max b with nuclear V_PT read(*,*)npotN,bmaxN,lamax c npotN= c 0: purely Coulomb (point-point) c 10: Woods-Saxon optical potential with point-spere Coulomb c 20: Gaussian potential (no Coulomb) c 100: chi phaseshift function given in a file If(lamax>=0)then !sum multipoles from 0 up to lamax lamax=lamax lamin=0 Else !only lamax term lamax=-lamax lamin=lamax Endif If(npotN<10)then !no potential to be read read(*,*) read(*,*) VcT=0 WcT=0 WDcT=0 VfT=0 WfT=0 WDfT=0 ElseIf(npotN==10)then !optical potentials c Potentiel c-T read(*,*)VcT,WcT,WDcT,rRcT,rIcT,rDcT,aRcT,aIcT,aDcT,rCcT c Potentiel f-T read(*,*)VfT,WfT,WDfT,rRfT,rIfT,rDfT,aRfT,aIfT,aDfT,rCfT ElseIf(npotN<30)then !Gaussian potentials WcT=0 WDcT=0 VfT=0 WfT=0 WDfT=0 c Potentiel c-T read(*,*)VcT,aRcT c Potentiel f-T read(*,*)VfT,aRfT ElseIf(npotN==100)then read(*,*)fichcT fichcT=trim(fichcT) read(*,*)fichfT fichfT=trim(fichfT) open(unit=101,status='unknown',file=fichcT) open(unit=102,status='unknown',file=fichfT) EndIf c Discretisation meshes c Radial and angular meshes: read(*,*)Nr,rm,Nt,Np if(mod(Nt,2)/=0)then write(*,*)'Nt must be even: Nt=',Nt STOP endif if(mod(Np,2)/=0)then write(*,*)'Np must be even: Np=',Np STOP endif hr=rm/Nr c Z and d meshes read(*,*)NZ,hZ,Ns,hs c b mesh (impact parameter) c Discretisation meshes b dependent: Nzb=Nb of b zones, bmax= max b read(*,*)Nzb,bmax,hbpol allocate(bmin(Nzb),hb(Nzb),Nob(Nzb)) Nb=0 !total Nb of b izb=1 DO read(*,*)bmin(izb),hb(izb) IF(bmin(izb)>=bmax)then if(izb==1)then Nob(izb)=1 hb=1 else Nob(izb)=0 hb(izb)=0 endif do izb2=izb+1,Nzb !skip next (useless) lines read(*,*)bmin(izb2),hb(izb2) Nob(izb2)=0 hb(izb2)=0 enddo Nzb=max(1,izb-1) EXIT ELSE If(izb>1)then Nob(izb-1)=(bmin(izb)-bmin(izb-1))/hb(izb-1) Nb=Nb+Nob(izb-1) if(Nob(izb-1)<0)then write(*,'(''Invalid bmin in b zone'',i3)')izb GOTO 9999 elseif(abs(bmin(izb-1)+Nob(izb-1)*hb(izb-1) a -bmin(izb))>1d-5)then write(*,'(''Invalid hb in b zone'',i3)')izb-1 GOTO 9999 endif EndIf ENDIF if(izb==Nzb) EXIT izb=izb+1 ENDDO If(bmin(Nzb)str c 100<=lec<110: only first-order Coulomb FF approximation if(lec>=100.and.lec<110)npotN=0 ilec=mod(lec,10) c Name of output files c wave function in 'fichfo' & additional calculations in 'fichres' read(*,'(a80)')fich if(ilec<7)then fichsk=trim(fich)//'.sk' else fichs0=trim(fich)//'.s0' endif select case (ilec) case(0) fichres=trim(fich)//'.sdb' NE=1 hE=1 case(1) fichres=trim(fich)//'.sdE' case(2) fichres=trim(fich)//'.sdp' case(3) fichres=trim(fich)//'.sdt' case(7) fichres=trim(fich)//'.str' NE=0 hE=1 case(8) fichres=trim(fich)//'.but' NE=0 hE=1 case(9) fichres=trim(fich)//'.sel' NE=0 hE=1 end select c Printing information about the calculation inputs c (value of the constants, potentials used, numerical parameters,...) write(*,7001)hc,e2,hm 7001 format(' hc=',3pd22.14,' e2=',1pd22.14,' hm=',2pd22.14) write(*,7002)Ac,Af,deuxm 7002 format('Ac=',f5.2,' Af=',f5.2,' deuxm=',1pd22.14) write(*,7003)kZc,kZf 7003 format('Zc=',i2,' Zf=',i2) write(*,7004)iSc,iS 7004 format('Sc=',i2,' Sf=',i2) if(ntypo<10)then write(*,*)' Woods-Saxon Potential' select case (ntypo) case(1) do ipo=1,Njpo-1 write(*,7006)Jpo(ipo),lpo(ipo), b Vp(ipo),rp(ipo),ap(ipo), a Vls(ipo),rls(ipo),als(ipo),rC(ipo) enddo 7006 format(' J=',i3,'l=',i3,' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3, a ' Vls=',f7.3,' rls=',f7.3,' als=',f7.3,' rC=',f7.3) write(*,*)'Others:' write(*,7005)Vp(Njpo),rp(Njpo),ap(Njpo), a Vls(Njpo),rls(Njpo),als(Njpo),rC(Njpo) 7005 format(' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,' Vls=',f7.3, a ' rls=',f7.3,' als=',f7.3,' rC=',f7.3) case(2) write(*,7035)Vp(1),rp(1),ap(1),Vls(1),rls(1),als(1),rC(1) 7035 format(' Odd-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3, a ' Vls=',f7.3,' rls=',f7.3,' als=',f7.3,' rC=',f7.3) write(*,7025)Vp(2),rp(2),ap(2),Vls(2),rls(2),als(2),rC(2) 7025 format('Even-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3, a ' Vls=',f7.3,' rls=',f7.3,' als=',f7.3,' rC=',f7.3) end select elseif(ntypo<20)then write(*,*)' Gaussian Potential' select case (ntypo) case(11) do ipo=1,Njpo-1 write(*,7206)Jpo(ipo),lpo(ipo), b Vp(ipo),rp(ipo),ap(ipo), a Vls(ipo),rls(ipo),als(ipo),rC(ipo) enddo 7206 format(' J=',i3,'l=',i3,' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3, a ' Vpp=',f7.3,' rpp=',f7.3,' app=',f7.3,' rC=',f7.3) write(*,*)'Others:' write(*,7205)Vp(Njpo),rp(Njpo),ap(Njpo), a Vls(Njpo),rls(Njpo),als(Njpo),rC(Njpo) 7205 format(' Vp=',f8.3,' rp=',f7.3,' ap=',f7.3,' Vpp=',f7.3, a ' rpp=',f7.3,' app=',f7.3,' rC=',f7.3) case(12) write(*,7235)Vp(1),rp(1),ap(1),Vls(1),rls(1),als(1),rC(1) 7235 format(' Odd-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3, a ' Vpp=',f7.3,' rpp=',f7.3,' app=',f7.3,' rC=',f7.3) write(*,7225)Vp(2),rp(2),ap(2),Vls(2),rls(2),als(2),rC(2) 7225 format('Even-l: Vp=',f8.3,' rp=',f7.3,' ap=',f7.3, a ' Vpp=',f7.3,' rpp=',f7.3,' app=',f7.3,' rC=',f7.3) end select endif write(*,7007)Nel 7007 format('Nel=',i2) write(*,7008)(iel,NRel(iel),Jel(iel),lel(iel),Eel(iel),iel=1,Nel) 7008 format('state',i2,': NR=',i2,' J= ',i2,' l= ',i2,' E0= ',f8.4) write(*,7011)ATa,kZT,EAP,v0/hc,VCK 7011 format('ATa=',f6.2,' ZT=',i2,' E/A=',f7.3,' v0/c=',f6.3, a ' VCK=',1pd13.6) write(*,7042)npotN * if(npotN==10.or.npotN==11.or.npotN==12)then if(npotN==10)then write(*,7040)VcT,WcT,WDcT,rRcT,rIcT,rDcT,aRcT,aIcT,aDcT,rCcT write(*,7041)VfT,WfT,WDfT,rRfT,rIfT,rDfT,aRfT,aIfT,aDfT,rCfT * if(npotN==11) write(*,*)'No f-T Coulomb potential' * if(npotN==12) write(*,*)'No c-T Coulomb potential' 7042 format('npotN= ',i2) 7040 format('VcT=',f9.4,' WcT=',f9.4,' WDcT=',f9.4,' rRcT=',f7.4, a ' rIcT=',f7.4,' rDcT=',f7.4,' aRcT=',f7.4,' aIcT=',f7.4, b ' aDcT=',f7.4,' rCcT=',f7.4) 7041 format('VfT=',f9.4,' WfT=',f9.4,' WDfT=',f9.4,' rRfT=',f7.4, a ' rIfT=',f7.4,' rDfT=',f7.4,' aRfT=',f7.4,' aIfT=',f7.4, b ' aDfT=',f7.4,' rCfT=',f7.4) elseif(npotN==20)then write(2,7044)VcT,aRcT write(2,7045)VfT,aRfT 7044 format('VcT=',f9.4,' aRcT=',f7.4) 7045 format('VfT=',f9.4,' aRfT=',f7.4) elseif(npotN==100)then write(*,*)' Numerical eikonal phases read externally' endif write(*,7012)Nr,rm 7012 format('Nr=',i7,' rm=',f9.2) write(*,7212)lamax,Nt,Np 7212 format('lamax=',i3,' Nt=',i3,' Np=',i3) write(*,7013)NZ,hZ,Ns,hs 7013 format('NZ= ',i5,' hZ=',f6.3,' Ns= ',i5,' hs=',f6.3) write(*,7014)bmin(1),bmax,Nb 7014 format('bmin=',f6.2,' bmax=',f6.2,' Nbtot=',i8) do izb=1,Nzb write(*,7015)izb,bmin(izb),Nob(izb),hb(izb) enddo 7015 format('b zone ',i3,' bmin=',f6.2,' Nb=',i4,' hb=',f7.4) write(*,7016)lec 7016 format('lec=',i3) Nl=lel(1)+lamax Nlj=(iS+1)*(Nl+1) NlJM=(iS+1)*(Nl+1)**2 ALLOCATE(foel(Nr,Nel),Vcf(Nr,Nlj),S(Nr),fibb(Nb)) fibb(1)=0 ib=0 Do izb=1,Nzb ib=ib+1 If(mod(Nob(izb),2)==1)then !Simpson fibb(ib)=fibb(ib)+xb(ib)*hb(izb)/300 fibb(ib+1)=4*xb(ib+1)*hb(izb)/300 do ib2=3,Nob(izb)-1,2 ib=ib+2 fibb(ib)=2*xb(ib)*hb(izb)/300 fibb(ib+1)=4*xb(ib+1)*hb(izb)/300 enddo ib=ib+1 fibb(ib)=xb(ib)*hb(izb)/300 Else !Trapezes fibb(ib)=fibb(ib)+xb(ib)*hb(izb)/200 do ib2=2,Nob(izb)-1 ib=ib+1 fibb(ib)=xb(ib)*hb(izb)/100 enddo ib=ib+1 fibb(ib)=xb(ib)*hb(izb)/200 EndIf EndDo SELECT CASE(ilec) CASE(0:3) !breakup cross sections ALLOCATE(Enrj(NE),phasOL(NE,Nlj),foli(Nr,Nlj,NE)) ALLOCATE(SklJM(NlJM,NE,Nsed,Nb)) if(ilec==3)then !angular distribution ALLOCATE(fbu(Nq,Nsed,NlJM,NE)) ALLOCATE(theta(Nq)) Do iq=1,Nq theta(iq)=pi/180*(th0+(iq-1)*hth) EndDo allocate(ADlj(Nq,0:Nlj)) endif Do iE=1,NE !defining energies at which SklJM are computed If(ilec==2)then !parallel-momentum distribution Enrj(iE)=(E0+(iE-1)*hE)**2/deuxm Else !energy or angular distributions Enrj(iE)=E0+(iE-1)*hE EndIf EndDo CASE(7) !total one-neutron stripping cross section ALLOCATE(Pstr(Nb)) CASE(8) !total breakup cross section ALLOCATE(S20(Nb),Si02(Nb)) CASE(9) !elastic scattering cross section ALLOCATE(theta(Nq),sruth(Nq),fcoul(Nq),ftild(Nq)) ALLOCATE(S0M(Nb,Nef,Nsed),selast(Nq)) ALLOCATE(S0iM(Nb,Nelm,Nsed),sinel(Nel,Nq)) Do iq=1,Nq theta(iq)=pi/180*(th0+(iq-1)*hth) EndDo END SELECT c Computes the c-f potential and the bound-state, and continuum wave functions jlj=0 Do l=0,Nl Do iJ=abs(2*l-iS),2*l+iS,2 jlj=jlj+1 do ir=1,Nr !computes Vcf=the core-fragment potential r=ir*hr Vcf(ir,jlj)=vJ(r,l,iJ)*deuxm enddo ****************************************************** c elimination of states by SuSy * U=Vcf(:,jlj) Do ieli=1,Neli if(l==leli(ieli).AND.iJ==Jeli(ieli))then ins=l+1 call super(Nr,hr,ins,Eeli(ieli),NReli(ieli), a Vcf(:,jlj),S,eps) endif EndDo !ieli ****************************************************** do iel=1,Nel !computes bound states If(l==lel(iel).and.iJ==Jel(iel))then S2=0 Ek=0 S=Vcf(:,jlj) call num1l(Nr,hr,Ek,S2,S,foel(:,iel),NRel(iel),eps) if(abs(1-Ek/deuxm/Eel(iel))>epsrel)then write(*,8001)iel,Eel(iel),Ek/deuxm goto 9998 else Eel(iel)=Ek/deuxm write(*,8002)iel,Eel(iel) endif EndIf enddo do iE=1,NE !computes continuum states Ek=sqrt(Enrj(iE)*deuxm) eta=kZc*kZf*e2*deuxm/2/Ek c sigma_l If(eta/=0)then if(l==0)then cg=1+ci*eta cg=cgamma(cg) phasOL(iE,:)=atan2(aimag(cg),real(cg)) elseif(iJ==abs(2*l-iS))then phasOL(iE,jlj:NlJ)=phasOL(iE,jlj:NlJ)+atan(eta/l) endif Else phasOL(iE,jlj)=0 EndIf c Calcul de l'onde distordue sur le reseau uniforme. 'dephase' normalise c la fct d'onde a cos(dl)*Fl(eta,k*r)+sin(dl)*Gl(eta,k*r) pour r->infini call dephase0(Nr,hr,Vcf(:,jlj),foli(:,jlj,iE),eps,dep, a l,eta,Ek,rm) c Ajout du dephasage nucleaire a la phase et passage en exp. phasOL(iE,jlj)=phasOL(iE,jlj)+dep-mod(l,4)*pi/2 phasOL(iE,jlj)=exp(ci*phasOL(iE,jlj)) enddo !iE EndDo !iJ EndDo !l 8001 Format('Error in binding energy of state',i2,2d16.8) 8002 Format('Binding energy of state',i2,'=',1pd16.8,'MeV') c Computes the eikonal phase with multipole expansion of the PT interaction call multeik if(npotN==100)then close(unit=101) close(unit=102) endif open(unit=2,status='unknown',file=fichres) SELECT CASE (ilec) CASE(0:3) c Computes the Skljm SklJM=0 cE=sqrt((1d0*iJ0+1d0)*(2d0*l0+1d0)/4/pi) DO 1 ib=1,Nb DO 2 ief=1,Nsed iM0=-iJ0+2*ief-2 DO 3 iE=1,NE Egv=(Enrj(iE)-Eel(1))*gaminv/v0 !relativistic correction x=Egv*xb(ib) if(x==0)then fp = 0 yp = 0 else call beskn(-x,nbes,fp,yp) endif If(lec<10)then !Usual Eikonal => no correction cor2=0 cor3=0 Else !Coulomb Corrected Eikonal cor2=(qc+qf)*Egv*yp(2) cor3=(qc+qf)*Egv*yp(1) EndIf jlJ=0 jlJM=0 Do l=0,Nl Do iJ=abs(2*l-iS),2*l+iS,2 jlJ=jlJ+1 Do iM=-iJ,iJ,2 jlJM=jlJM+1 iye=0 iyo=0 mu=(iM-iM0)/2 !fixed by Q3m if(mu<0.and.mod(mu,2)==1)then mfac=-1 else mfac=1 endif do la=0,lamax if(abs(mu)>la)goto 421 call troisj(2*l,2*la,2*l0,0,0,0,Q3) if(Q3==0)goto 421 call sixj(iJ,2*l,iS,2*l0,iJ0,2*la,Q6) if(Q6==0)goto 421 call troisj(iJ0,2*la,iJ,iM0,2*mu,-iM,Q3m) if(Q3m==0)goto 421 rint=0 If(mod(la+mu,2)==0)then !only feik 1 and 2 jye=iye+abs(mu)/2+1 !verif OK do jr=1,Nr !radial integral rint=rint+foel(jr,1)*foli(jr,jlj,iE)* a (feik(jye,jr,ib,1) b +ci*cor2*feik(jye,jr,ib,2)) enddo Else !only feik 3 jyo=iyo+abs(mu)/2+1 do jr=1,Nr !radial integral rint=rint-foel(jr,1)*foli(jr,jlj,iE)* a cor3*feik(jyo,jr,ib,3) enddo EndIf SklJM(jlJM,iE,ief,ib)=SklJM(jlJM,iE,ief,ib) a +rint*hr*Q3*Q6*Q3m*sqrt(2d0*la+1d0) 421 continue iye=iye+la/2+1 iyo=iyo+(la+1)/2 enddo !la if(mod(iS-iM,4)==0)then SklJM(jlJM,iE,ief,ib)=SklJM(jlJM,iE,ief,ib) a *mfac*cE*phasOL(iE,jlJ) b *sqrt((1d0*iJ+1d0)*(2d0*l+1d0)) else SklJM(jlJM,iE,ief,ib)=-SklJM(jlJM,iE,ief,ib) a *mfac*cE*phasOL(iE,jlJ) b *sqrt((1d0*iJ+1d0)*(2d0*l+1d0)) endif EndDo !iM EndDo !iJ EndDo !il 3 continue !iE 2 continue !m0 1 continue !b *ERSk c Printing of SklJM * open(unit=21,status='unknown',file=fichsk) * do ib=1,Nb * b=b0+(ib-1)*hb * write(21,*)b,(abs(SklJM(icomp,1,1,ib)), * a phase(SklJM(icomp,1,1,ib)),icomp=1,NlJM) * enddo * close(21) Select Case(ilec) Case(0,1) ALLOCATE(sedbdl(0:Nlj),sdE(0:Nlj)) c Energy distribution or Pbu(b) at 1 energy DO iE=1,NE Ek=sqrt(deuxm*Enrj(iE)) sdE=0 Do ib=1,Nb sedbdl=0 do ief=1,Nsed if(ief==Nsed.and.lzf==2)then !m0=0 c factor for the summation over m0 m0f=1 else m0f=2 endif jlJM=0 jlJ=0 Do l=0,Nl Do iJ=abs(2*l-iS),2*l+iS,2 jlj=jlj+1 slj=0 do iM=-iJ,iJ,2 jlJM=jlJM+1 sedbdl(jlj)=sedbdl(jlj)+ a m0f*abs(SklJM(jlJM,iE,ief,ib))**2 sedbdl(0)=sedbdl(0)+ a m0f*abs(SklJM(jlJM,iE,ief,ib))**2 enddo !iM EndDo !iJ EndDo !l enddo !m0 sedbdl=sedbdl*deuxm/pi/Ek/Nef *ERb If(ilec==0)then !thus NE=1 write(2,3003)xb(ib),(sedbdl(il),il=0,Nlj) Else c Integration over b sdE=sdE+fibb(ib)*sedbdl EndIf EndDo !ib *ERE If(ilec==1)then sdE=2*pi*sdE write(2,3003)Enrj(iE),(sdE(il),il=0,Nlj) EndIf ENDDO !iE DEALLOCATE(sedbdl,sdE) Case(2) c Parallel momentum distribution in the projectile CM rest frame ALLOCATE(sedp(-NE:NE)) sedp=0 do ib=1,Nb c Integration over b do ief=1,Nsed if(ief==Nsed.and.lzf==2)then c factor for the summation over m0 m0f=1 else m0f=2 endif DO ikp=0,NE if (ikp==0)then Ekp=0 else Ekp=E0+(ikp-1)*hE endif DO ik=max(ikp,1),NE Ek=E0+(ik-1)*hE costk=Ekp/Ek call HARSPH(costk,-Nl-1,Pl,Ylm) smmsp=0 smmsm=0 Do iM=-2*Nl-iS,2*Nl+iS,2 Do ims=-iS,iS,2 iml=iM-ims csljp=0 csljm=0 do il=abs(iml),2*Nl,2 l=il/2 do iJ=abs(il-iS),il+iS,2 If(iJ>=abs(iM))then jlJM=jfromlSJM(l,iS,iJ,iM) call clebs(il,iS,iJ,iml,ims,iM,CGO) jYlm=l*(l+1) jYlm=jYlm/2+abs(iml)/2+1 c NB: HARSPH gives the Ylm for m>=0 only, and phi=0 c values for m<0 are the same with a phase (-)^m c this phase is irrelevant here because of the |...|^2 in the sum over m csljp=csljp+CGO*Ylm(jYlm) a *SklJM(jlJM,ik,ief,ib) if(mod(il-iml,4)==0)then csljm=csljm+CGO*Ylm(jYlm) a *SklJM(jlJM,ik,ief,ib) else csljm=csljm-CGO*Ylm(jYlm) a *SklJM(jlJM,ik,ief,ib) endif EndIf enddo !iJ enddo !l smmsp=smmsp+abs(csljp)**2 smmsm=smmsm+abs(csljm)**2 EndDo !ims EndDo !iM c Integral over k by Simpson's rule c (breakup probabilties above Nk are assumed to be negligible) c Sum over m0 and integral over b (fibb) are included if(ikp==0)then if(ik==1)then !first step by trapezes sedp(0)=sedp(0) a +fibb(ib)*m0f*smmsp/Ek*(E0/2+hE/3) elseif(mod(ik,2)==0)then !rest by Simpson sedp(0)=sedp(0) a +fibb(ib)*m0f*4*smmsp*hE/Ek/3 else sedp(0)=sedp(0) a +fibb(ib)*m0f*2*smmsp*hE/Ek/3 endif elseif(ik==ikp)then sedp(ikp)=sedp(ikp) a +fibb(ib)*m0f*smmsp*hE/Ek/3 sedp(-ikp)=sedp(-ikp) a +fibb(ib)*m0f*smmsm*hE/Ek/3 elseif(mod(ik-ikp,2)==1)then sedp(ikp)=sedp(ikp) a +fibb(ib)*m0f*4*smmsp*hE/Ek/3 sedp(-ikp)=sedp(-ikp) a +fibb(ib)*m0f*4*smmsm*hE/Ek/3 else sedp(ikp)=sedp(ikp) a +fibb(ib)*m0f*2*smmsp*hE/Ek/3 sedp(-ikp)=sedp(-ikp) a +fibb(ib)*m0f*2*smmsm*hE/Ek/3 endif ENDDO !ik ENDDO !ikp enddo !m0 enddo !ib sedp=8*pi*sedp/Nef *ERp Do ikp=-NE,NE if(ikp<0)then write(2,'(1f8.4,1pd18.8)')-E0+(ikp+1)*hE,sedp(ikp) elseif(ikp==0)then write(2,'(1f8.4,1pd18.8)')ikp*hE,sedp(ikp) else write(2,'(1f8.4,1pd18.8)')E0+(ikp-1)*hE,sedp(ikp) endif EndDo DEALLOCATE (sedp) Case(3) c Angular distribution including integration over energy call calculFBU ADlj=0 do iq=1,Nq Do ief=1,Nsed if(ief==Nsed.and.lzf==2)then !m0=0 c factor for the summation over m0 m0f=1 else m0f=2 endif icomp=0 jlJ=0 Do l=0,Nl Do iJ=abs(2*l-iS),2*l+iS,2 jlj=jlj+1 do iM=-iJ,iJ,2 icomp=icomp+1 * if(NE==1)hE=1.0d0!si rectangles if(NE==1)hE=2 !si trapezes do iE=1,NE cbu=deuxm/pi/sqrt(deuxm*Enrj(iE))*EKPT**2/100 ! in barn ! Integration par trapezes if(iE==1.or.iE==NE)then ADlj(iq,jlj)=ADlj(iq,jlj)+m0f*hE/2*cbu a *(abs(fbu(iq,ief,icomp,iE)))**2 ADlj(iq,0)=ADlj(iq,0)+m0f*hE/2*cbu a *(abs(fbu(iq,ief,icomp,iE)))**2 else ADlj(iq,jlj)=ADlj(iq,jlj)+m0f*hE*cbu a *(abs(fbu(iq,ief,icomp,iE)))**2 ADlj(iq,0)=ADlj(iq,0)+m0f*hE*cbu a *(abs(fbu(iq,ief,icomp,iE)))**2 endif enddo !iE enddo !iM EndDo !iJ EndDo !l enddo !ief ADlj(iq,:)=ADlj(iq,:)/Nef *ERt * write(2,'(1pd10.4,60pd18.8)')180/pi*theta(iq),AD(iq), write(2,3003)180/pi*theta(iq),(ADlj(iq,jlj),jlj=0,Nlj) enddo DEALLOCATE (fbu,theta,ADlj) DEALLOCATE (Enrj,phasOL,foli,SklJM) End Select !ilec for breakup distributions CASE(7) c Computes the total one-neutron stripping cross section c ie total inelastic one-neutron removal cross section Pstr=0 il0=2*lel(1) iJ0=Jel(1) cstr=(1d0*il0+1d0)/2/sqrt(pi) if(mod(iJ0,2)/=0)cstr=-cstr iye=1 Do la=0,min(lamax,il0,iJ0),2 !only even la because of Q3 ila=2*la call troisj(il0,ila,il0,0,0,0,Q3) if(Q3==0)goto 827 call sixj(iJ0,il0,iS,il0,iJ0,ila,Q6) if(Q6==0)goto 827 sum0=0 do im0=-iJ0,iJ0,2 call troisj(iJ0,ila,iJ0,iM0,0,-iM0,Q3j) if(mod(iS+iM0,2)==0)then sum0=sum0+Q3j else sum0=sum0-Q3j endif enddo DO ib=1,Nb rint=0 !only feik 1 because no correction E1 do jr=1,Nr !radial integral if(la==0)then feik1=feik(iye,jr,ib,1)+sqrt(4*pi) else feik1=feik(iye,jr,ib,1) endif rint=rint+(foel(jr,1))**2*feik1 enddo !jr Pstr(ib)=Pstr(ib)+rint*hr*Q3*Q6*Q3j*cstr ENDDO !b 827 continue iye=iye+la+2 EndDo !la open(unit=20,status='unknown',file=fichs0) do ib=1,Nb write(20,'(1f10.4,2(1pd18.8))')xb(ib), a real(Pstr(ib)),aimag(Pstr(ib)) enddo close(unit=20) c Computes total one-neutron stripping cross section (integral over b) Str=0 do ib=2,Nb-1 Str=Str+fibb(ib)*real(Pstr(ib)) enddo Str=2*pi*Str write(*,'(''Total stripping cross section='',1pd18.8)')Str deallocate(Pstr) *ERbut write(2,'(1pd18.8)')Str CASE(8) c Computes the total breakup cross section through closing relation. c 1. Computes the norm of the wave function S20=0 il0=2*l0 CS2=sqrt((1d0*il0+1d0)/(1d0*iJ0+1d0))/2/pi if(mod(il0+iS+iJ0,4)/=0)Cs2=-Cs2 iLAmax=2*min(il0,iJ0) DO 800 iLA=0,iLAmax,4 !LA=2*Lambda, and Lambda is even by Q3l if(lzf==2)then !case iM0=0 call clebs(iJ0,iLA,iJ0,0,0,0,Q3jm) Q3jm=Q3jm/2 else !case iM0=1 call clebs(iJ0,iLA,iJ0,1,0,1,Q3jm) endif do iM0=iJ0,2,-2 call clebs(iJ0,iLA,iJ0,iM0,0,iM0,Q3) Q3jm=Q3jm+Q3 enddo call clebs(il0,iLA,il0,0,0,0,Q3l) if(Q3l==0)goto 800 call sixj(iJ0,il0,iS,il0,iJ0,iLA,Q6) if(Q6==0)goto 800 iye1=0 Do la1=0,lamax !lambda ila1=2*la1 iye2=0 Do la2=0,lamax !lambda prime ila2=2*la2 call clebs(ila2,ila1,iLA,0,0,0,Q3lala) if(Q3lala==0)goto 802 Q3lala=Q3lala*sqrt((1d0*ila1+1d0)*(1d0*ila1+1d0)) do 803 mu=-la1,la1 if(abs(mu)>la2)goto 803 call clebs(ila2,ila1,iLA,-2*mu,2*mu,0,Q3lamu) if(Q3lamu==0)goto 803 DO ib=1,Nb rint=0 If(mod(la1+mu,2)==0.and.mod(la2+mu,2)==0)then !only feik 1 because the E1 correction is neglected to compute norm of wf jye1=iye1+abs(mu)/2+1 jye2=iye2+abs(mu)/2+1 do jr=1,Nr if(la1==0.and.mu==0)then feik1=feik(jye1,jr,ib,1)+sqrt(4*pi) else feik1=feik(jye1,jr,ib,1) endif if(la2==0.and.mu==0)then feik2=conjg(feik(jye2,jr,ib,1)+sqrt(4*pi)) else feik2=conjg(feik(jye2,jr,ib,1)) endif rint=rint+foel(jr,1)**2*feik1*feik2 enddo if(mod(mu,2)==0)then S20(ib)=S20(ib) a +rint*hr*Q3lamu*Q3lala*Q6*Q3l*Q3jm else S20(ib)=S20(ib) a -rint*hr*Q3lamu*Q3lala*Q6*Q3l*Q3jm endif EndIf ENDDO !ib 803 continue !mu 802 Continue iye2=iye2+la2/2+1 EndDo !la2 iye1=iye1+la1/2+1 EndDo !la1 800 CONTINUE !LA S20=S20*CS2 c 2. Computes the S0 for the P-T elastic scattering (simple modif of CASE(9)) Si02=0 DO iel=1,Nel il=2*lel(iel) iJ=Jel(iel) cSi0=(1d0*il0+1d0)*(1d0*il+1d0)*(1d0*iJ+1d0)/4/pi iye=0 Do la=0,lamax ila=2*la if(ilail+il0)goto 820 call troisj(il,ila,il0,0,0,0,Q3) if(Q3==0)goto 820 call sixj(iJ,il,iS,il0,iJ0,ila,Q6) if(Q6==0)goto 820 do mu=-la,la DO ib=1,Nb rint=0 If(mod(la+mu,2)==0)then !only feik 1 because no correction E1 jye=iye+abs(mu)/2+1 do jr=1,Nr !radial integral if(la==0.and.mu==0)then feik1=feik(jye,jr,ib,1)+sqrt(4*pi) else feik1=feik(jye,jr,ib,1) endif rint=rint+foel(jr,iel)*foel(jr,1)*feik1 enddo EndIf Si02(ib)=Si02(ib)+(abs(rint*hr*Q3*Q6))**2*cSi0 ENDDO !ib enddo 820 continue iye=iye+la/2+1 EndDo !la ENDDO !iel open(unit=20,status='unknown',file=fichs0) do ib=1,Nb write(20,'(1f10.4,3(1pd18.8))')xb(ib),real(S20(ib)), a aimag(S20(ib)),Si02(ib) enddo close(unit=20) c Computes total breakup cross section (integral over b) Sbut=0 do ib=1,Nb Sbut=Sbut+fibb(ib)*(real(S20(ib))-Si02(ib)) enddo Sbut=2*pi*Sbut write(*,*)"computing bu" write(*,'(''Total breakup cross section='',1pd18.8)')Sbut deallocate(S20,Si02) *ERbut write(2,'(1pd18.8)')Sbut CASE (9) c Computes the S0 for the P-T elastic scattering S0M=0 cE=(1d0*iJ0+1d0)*(2d0*l0+1d0)/2/sqrt(pi) Do ief=1,Nsed iM0=-iJ0+2*ief-2 do ief2=1,Nef iM=-iJ0+2*ief2-2 mu=(iM-iM0)/2 !fixed by Q3m iye=0 if(mu<0.and.mod(mu,2)==1)then mfac=-1 else mfac=1 endif do la=0,lamax if(abs(mu)>la)goto 429 call troisj(2*l0,2*la,2*l0,0,0,0,Q3) if(Q3==0)goto 429 call sixj(iJ0,2*l0,iS,2*l0,iJ0,2*la,Q6) if(Q6==0)goto 429 call troisj(iJ0,2*la,iJ0,iM0,2*mu,-iM,Q3m) if(Q3m==0)goto 429 DO ib=1,Nb rint=0 If(mod(la+mu,2)==0)then !only feik 1 because no correction E1 jye=iye+abs(mu)/2+1 do jr=1,Nr !radial integral c No correction E1 at first order since elastic (l=l0, and J=J0) rint=rint+foel(jr,1)**2*feik(jye,jr,ib,1) enddo EndIf S0M(ib,ief2,ief)=S0M(ib,ief2,ief) a +rint*hr*Q3*Q6*Q3m*sqrt(2d0*la+1d0) ENDDO !b 429 continue iye=iye+la/2+1 enddo !la if(mod(iS-iM,4)==0)then S0M(:,ief2,ief)=S0M(:,ief2,ief)*mfac*cE else S0M(:,ief2,ief)=-S0M(:,ief2,ief)*mfac*cE endif enddo !ief2 (ie m) EndDo !ief (ie m0) *ERS0 c Printing of S0 open(unit=20,status='unknown',file=fichs0) do ib=1,Nb write(20,'(1f10.4,8(1pd18.8))')xb(ib),(abs(S0M(ib,icomp,1)), a phase(S0M(ib,icomp,1)),icomp=1,Nef) enddo !ib close(20) c Computes sigma_0 for PT motion cg=1+ci*etaPT cg=cgamma(cg) s0PT=atan2(aimag(cg),real(cg)) !OK c Computes Rutherford cross section 'sruth', c and Coulomb elastic scattering amplitude 'fcoul' Do iq=1,Nq sruth(iq)=(etaPT/EKPT)**2/sin(theta(iq)/2)**4/400 fcoul(iq)=-ci*etaPT/2.0d0/EKPT**2/sin(theta(iq)/2)**2 a *exp(-2*ci*etaPT*log(2*EKPT*sin(theta(iq)/2)) b +2*ci*S0PT) EndDo c Computes elastic scattering cross section open(unit=20,status='unknown',file=fichs0) nla=20 !nb of points for the Gauss integration in 'intbesm' NJmax=iJ0 Allocate(fBessJ(NJmax+1),fBessY(NJmax+1)) Do ief=1,Nsed iM0=-iJ0+(ief-1)*2 !2m0 Do ief2=1,Nef iM02=-iJ0+(ief2-1)*2 !2m0' ibess=abs(iM0-iM02)/2+1 ftild=0 DO ib=1,Nb-1 Npol=nint((xb(ib+1)-xb(ib))/hbpol) if(abs(Npol*1.0d0-(xb(ib+1)-xb(ib))/hbpol)>1d-3)then write(*,*)'Error in interpolation at ib=',ib write(*,*)'hb=',xb(ib+1)-xb(ib) write(*,*)'hbpol=',hbpol STOP endif do ipol=1,Npol bpol=xb(ib)+(ipol-1)*hbpol if(bpol<1d-3)then chic2=0 else * chic2=2*etaPT*log(bpol/2) c The Coulomb scattering phase is chosen equal to the c Coulomb eikonal phase that reproduces the exact Coulomb amplitude c cf. Bertulani & Danielewicz pg 423 Eq.(8.77). chic2=2*etaPT*log(EKPT*bpol) endif ctmp=S0M(ib,ief2,ief)+(S0M(ib+1,ief2,ief) a -S0M(ib,ief2,ief))*(ipol-1)/Npol do iq=1,Nq q=2*EKPT*sin(theta(iq)/2) * z=q*bpol !version w/ S17* x=q*bpol If(ibess==1)then !iM02==iM0 * BessJ0=S17AEF(q*bpol,IFAIL) * call BESSJ(NJmax,q*bpol,fBessJ,fBessY) call bessj0y0(x,fBessJ(1),fBessY(1)) ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp a *bpol*exp(ci*chic2)*hbpol ElseIf(ibess==2)then !iM02=iM0+-1 fBessJ(2)=bessj1(x) ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp a *bpol*exp(ci*chic2)*hbpol Else * z=q*bpol * call S17DEF(0.0d0,z,NJmax,'U',cBessJ,NZ,IFAIL) call BESSJ(NJmax,x,fBessJ,fBessY) * ftild(iq)=ftild(iq)+cBessJ(ibess+1)*ctmp * a *bpol*exp(ci*chic2)*hbpol ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp a *bpol*exp(ci*chic2)*hbpol Endif enddo !iq enddo !bpol ENDDO !b Do iq=1,Nq q=2*EKPT*sin(theta(iq)/2) call intbesm(ibess-1,bmax,q,etaPT,nla,cres) ftild(iq)=ftild(iq)+S0M(Nb,ief2,ief)*bmax**2*cres If(iM0==iM02)then selast(iq)=selast(iq)+2*abs(fcoul(iq)+ftild(iq))**2 Else selast(iq)=selast(iq)+2*abs(ftild(iq))**2 Endif EndDo EndDo !ief2 (ie m) EndDo !ief(ie m0) Deallocate(fBessJ,fBessY) close(unit=20) do iq=1,Nq selast(iq)=EKPT**2*selast(iq)/Nef/100 enddo ********************************************* c Computes the S0iM S0iM=0 cE=sqrt((1d0*iJ0+1d0)*(2d0*l0+1d0)/4/pi) DO ib=1,Nb DO ief=1,Nsed iM0=-iJ0+2*ief-2 ielm=Nef DO iel=2,Nel Egv=(Eel(iel)-Eel(1))*gaminv/v0 !relativistic correction x=Egv*xb(ib) if(x==0)then fp = 0 yp = 0 else call beskn(-x,nbes,fp,yp) endif If(lec<10)then !Usual Eikonal => no correction cor2=0 cor3=0 Else !Coulomb Corrected Eikonal cor2=(qc+qf)*Egv*yp(2) cor3=(qc+qf)*Egv*yp(1) EndIf l=lel(iel) iJ=Jel(iel) Do iM=-iJ,iJ,2 ielm=ielm+1 iye=0 iyo=0 mu=(iM-iM0)/2 !fixed by Q3m if(mu<0.and.mod(mu,2)==1)then mfac=-1 else mfac=1 endif do la=0,lamax if(abs(mu)>la)goto 422 call troisj(2*l,2*la,2*l0,0,0,0,Q3) if(Q3==0)goto 422 call sixj(iJ,2*l,iS,2*l0,iJ0,2*la,Q6) if(Q6==0)goto 422 call troisj(iJ0,2*la,iJ,iM0,2*mu,-iM,Q3m) if(Q3m==0)goto 422 rint=0 If(mod(la+mu,2)==0)then !only feik 1 and 2 jye=iye+abs(mu)/2+1 !verif OK do jr=1,Nr !radial integral rint=rint+foel(jr,1)*foel(jr,iel)* a (feik(jye,jr,ib,1) b +ci*cor2*feik(jye,jr,ib,2)) enddo Else !only feik 3 jyo=iyo+abs(mu)/2+1 !verif OK do jr=1,Nr !radial integral rint=rint-foel(jr,1)*foel(jr,iel)* a cor3*feik(jyo,jr,ib,3) enddo EndIf S0iM(ib,ielm,ief)=S0iM(ib,ielm,ief) a +rint*hr*Q3*Q6*Q3m*sqrt(2d0*la+1d0) 422 continue iye=iye+la/2+1 iyo=iyo+(la+1)/2 enddo !la if(mod(iS-iM,4)==0)then S0iM(ib,ielm,ief)=S0iM(ib,ielm,ief) a *mfac*cE b *sqrt((1d0*iJ+1d0)*(2d0*l+1d0)) else S0iM(ib,ielm,ief)=-S0iM(ib,ielm,ief) a *mfac*cE b *sqrt((1d0*iJ+1d0)*(2d0*l+1d0)) endif EndDo !iM ENDDO !iel ENDDO !iM0 ENDDO !ib c Calcul de Sinel write(*,*)'Calcul de la SINEL' sinel=0 hbpol=0.01d0 NJmax=Jel(1)+maxval(Jel)+1 Allocate(fBessJ(NJmax+1),fBessY(NJmax+1)) Do ief=1,Nsed iM0=-Jel(1)+(ief-1)*2 !2m0 ielm=Nef Do iel=2,Nel Do iM02=-Jel(iel),Jel(iel),2 !2m ielm=ielm+1 ibess=abs(iM0-iM02)/2+1 ftild=0 do ib=1,Nb-1 Npol=nint((xb(ib+1)-xb(ib))/hbpol) if(abs(Npol*1.0d0-(xb(ib+1)-xb(ib))/hbpol)>1d-3)then write(*,*)'Error in interpolation at ib=',ib write(*,*)'hb=',xb(ib+1)-xb(ib) write(*,*)'hbpol=',hbpol STOP endif do ipol=0,Npol-1 bpol=xb(ib)+ipol*hbpol ctmp=(S0iM(ib+1,ielm,ief)-S0iM(ib,ielm,ief)) a *ipol/Npol + S0iM(ib,ielm,ief) !phase menant a l'amplitude de diffusion Coulombienne exacte chic2=2*etaPT*log(EKPT*bpol) if(abs(bpol)<1e-3)chic2=0 do iq=1,Nq q=2*EKPT*sin(theta(iq)/2) * z=q*bpol * call S17DEF(0.0d0,z,NJmax,'U',BessJ,NZ,IFAIL) * ftild(iq)=ftild(iq)+BessJ(ibess)*ctmp * a *bpol*exp(ci*chic2)*hbpol * z=q*bpol !version w/ S17DEF * if(z==0)z=1d-15 x=q*bpol if(x==0)x=1d-15 If(ibess==1)then !iM02==iM0 call bessj0y0(x,fBessJ(1),fBessY(1)) ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp a *bpol*exp(ci*chic2)*hbpol ElseIf(ibess==2)then !iM02=iM0+-1 fBessJ(2)=bessj1(x) ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp a *bpol*exp(ci*chic2)*hbpol Else call BESSJ(NJmax,x,fBessJ,fBessY) ftild(iq)=ftild(iq)+fBessJ(ibess)*ctmp a *bpol*exp(ci*chic2)*hbpol Endif enddo !iq enddo !bpol enddo !ib do iq=1,Nq sinel(iel,iq)=sinel(iel,iq)+2*abs(ftild(iq))**2 enddo !iq EndDo !iM02 Enddo !iel Enddo !ief Deallocate(fBessJ,fBessY) Do iq=1,Nq sinel(:,iq)=EKPT**2*sinel(:,iq)/Nef/100 Enddo ********************************************* *ERel c Printing of elastic scattering cross section Do iq=1,Nq if(theta(iq).gt.1e-5)then write(2,'(1f8.4,10(1pd16.8))')180/pi*theta(iq), a sruth(iq),selast(iq),(sinel(iel,iq),iel=2,Nel) endif enddo DEALLOCATE(theta,sruth,fcoul,ftild) DEALLOCATE(S0M,S0iM,selast,sinel) END SELECT !ilect for reaction channel (el.scatt., total or distr. bu) 1001 format(/,' Breakup: ',a2,i2,' -> ',a2,i2,' + ',a2,i2,/) 1003 format(1h ,10f7.3) 1005 format(f11.6,4x,6(1pe14.5,2x)) 1006 format(1h ,f8.2,1p,5e13.4) 3000 format(2f10.6,2i3) 3001 format(3i3) 3003 format(f8.3,1p,60d16.8) 9998 continue DEALLOCATE(NRel,Jel,lel,Eel) DEALLOCATE(Jpo,lpo,Vp,rp,ap,rC,Vls,rls,als) DEALLOCATE(Vcf,foel,S,xb,fibb) close(1) close(2) 9999 continue *END END *multeik subroutine multeik c Computes the eikonal phase (including Coulomb correction) c Note the change in the definition of Vcs and Vfs from 'reneik12.f': c here they include i/hbar/v. In this way we can use directly the c chi phase-shift functions of Suzuki in the variable chi USE cstes USE collision USE potentials USE reseaufo USE reseaueik USE phaseik implicit real*8(a-h,o-z) real*8, parameter :: epslu=1d-4 real*8 :: xt(Nt),wt(Nt) complex*16 :: voptN,vopt,cpotct,cpotft,echin,echic,echi(3) complex*16 :: chincT,chinfT,chinPT complex*16 :: Vcs(0:Ns),Vfs(0:Ns),Vsn,Vsp real*8, allocatable, dimension (:):: yl,p real*8, allocatable, dimension (:,:):: ylm complex*16, allocatable, dimension (:,:):: f *** c Evaluate the dimension of the arrays, and allocate them if(mod(lamax,2)==0)then ny=(lamax/2+1)**2 else ny=(lamax/2+1)*(lamax/2+2) endif * print*,' ny =',ny nydim=(Nl+1)*(Nl+2)/2 allocate (ylm(nydim,Nt/2),feik(ny,Nr,Nb,3),f(ny,3)) allocate (p(Nl+1),yl(nydim)) *** c Computes the points and weights of the Gauss quadrature on theta ifail=0 zero=0 ! old: call d01bcf (0,-one,one,zero,zero,Nt,wt,xt,ifail) x1= -1.d0 x2= 1.d0 call gauleg(x1,x2,xt,wt,Nt) if(ifail.ne.0)print*,Nt,ifail c Computes the Ylm needed because of symmetries c (ie only phi=0, mu>=0, and la+mu even) Do it=1,Nt/2 !Nt is even to simplify the calculation x=xt(it) call harsph (x,-(lamax+1),p,yl) iy=0 do la=0,lamax do mu=0,la iy=iy+1 ylm(iy,it)=yl(iy) enddo !mu enddo !la enddo ! it IF(npotN>=10.and.npotN<100)then c Precalculation of the integral over Z of the P-T nuclear interaction c as a function of the perpendicular distance s c NB integrand symmetric in Z => *2 except for Z=0 c and invariant on the origin => no zz-z c Tested and OK for 11Be @70AMeV on lead for b=10-100fm, E=0-3MeV c with Ns=600 hs=0.05fm. Another interpolation could be better c eg 3pt instead of 2 to have at least the curve in the interpolation. Do js=0,Ns if(js==0)then s=1d-15 !at s=0 voptN diverges else s=js*hs endif if(npotN==10)then !optical potentials * Vcs(js)=vopt(VcT,rRcT,aRcT,WcT,rIcT,aIcT, * a WDcT,rDcT,aDcT,kZc,kZT,rCcT,s) * Vfs(js)=vopt(VfT,rRfT,aRfT,WfT,rIfT,aIfT, * a WDfT,rDfT,aDfT,kZf,kZT,rCfT,s) Vcs(js)=voptN(VcT,rRcT,aRcT,WcT,rIcT,aIcT, a WDcT,rDcT,aDcT,kZc,kZT,rCcT,s) Vfs(js)=voptN(VfT,rRfT,aRfT,WfT,rIfT,aIfT, a WDfT,rDfT,aDfT,kZf,kZT,rCfT,s) elseif(npotN==20)then !gaussian potentials Vcs(js)=vgauss(VcT,aRcT,s) Vfs(js)=vgauss(VfT,aRfT,s) endif do iz=1,Nz zz=iz*hz rxT=sqrt(s**2+zz**2) if(npotN==10)then !optical potentials * Vcs(js)=Vcs(js)+2*vopt(VcT,rRcT,aRcT,WcT,rIcT,aIcT, * a WDcT,rDcT,aDcT,kZc,kZT,rCcT,rxT) * Vfs(js)=Vfs(js)+2*vopt(VfT,rRfT,aRfT,WfT,rIfT,aIfT, * a WDfT,rDfT,aDfT,kZf,kZT,rCfT,rxT) Vcs(js)=Vcs(js)+2*voptN(VcT,rRcT,aRcT,WcT,rIcT,aIcT, a WDcT,rDcT,aDcT,kZc,kZT,rCcT,rxT) Vfs(js)=Vfs(js)+2*voptN(VfT,rRfT,aRfT,WfT,rIfT,aIfT, a WDfT,rDfT,aDfT,kZf,kZT,rCfT,rxT) elseif(npotN==20)then !gaussian potentials Vcs(js)=Vcs(js)+2*vgauss(VcT,aRcT,rxT) Vfs(js)=Vfs(js)+2*vgauss(VfT,aRfT,rxT) endif enddo !iz EndDo !js ! Old version where Vcs and Vfs contain only the integral of the potential * Vcs=Vcs*hz * Vfs=Vfs*hz ! New version where Vcs and Vfs contain also the -i/hbar/v factor Vcs=-ci*Vcs*hz/v0 Vfs=-ci*Vfs*hz/v0 * open(unit=27,file='Vopt.chi') * do js=0,Ns * write(27,*)js*hs,real(Vcs(js)),aimag(Vcs(is)), * a real(Vfs(js)),aimag(Vfs(is)) * enddo * close(unit=27) ELSEIF(npotN==100)then c The Chi phase-shift functions are precalculated from a GammaNN profile c function and read in separate files 101='fichcT' and 102='fichfT'. Do js=0,Ns s=js*hs * read(101,*)slu,Vsn,Vsp * write(*,*)s,Vsn,Vsp read(101,*)slu,VsnR,VsnI,VspR,VspI * write(*,*)s,VsnR,VsnI,VspR,VspI if(abs(slu-s)>epslu)then write(*,*)'Error in reading s in fichcT: ',s,slu STOP endif * Vcs(js)=-(Vsn+Vsp) !old version when there was a mistake in Suzuki's phaseshift functions * Vcs(js)=Vsn+Vsp Vcs(js)=VsnR+VspR+ci*(VsnI+VspI) * read(102,*)slu,Vsn,Vsp * write(*,*)s,Vsn,Vsp read(102,*)slu,VsnR,VsnI,VspR,VspI if(abs(slu-s)>epslu)then write(*,*)'Error in reading s in fichfT: ',s,slu STOP endif * Vfs(js)=-(Vsn+Vsp) !old version when there was a mistake in Suzuki's phaseshift functions * Vfs(js)=Vsn+Vsp Vfs(js)=VsnR+VspR+ci*(VsnI+VspI) EndDo ENDIF feik=0 c calculation of the eikonal phase corrected for Coulomb at first-order do 3 ib=1,nb !for each impact parameter if(xb(ib)==0)then b=1d-15 !at b=0 the Coulomb phase diverges else * b=b0+(ib-1)*hb b=xb(ib) endif do 4 ir=1,nr !for each radius r=ir*hr c f=0 do 5 it=1,Nt/2 !we only need 1/2 of the theta NB => Nt even? z=r*xt(it) sp=r*sqrt(1-xt(it)**2) do 6 ip=1,Np/2 !we only need 1/2 of the phi NB => Np even? phi=(2*ip-1)*pi/Np x=sp*cos(phi) y=sp*sin(phi) c Calculation of the nuclear eikonal phases IF(npotN==0)then !ie purely point-point Coulomb potential echin=1 !ie not nuclear potential * if(lec>=20)then * scT=sqrt((b+cc*x)**2+(cc*y)**2) * sfT=sqrt((b+cf*x)**2+(cf*y)**2) * endif ELSE scT=sqrt((b+cc*x)**2+(cc*y)**2) sfT=sqrt((b+cf*x)**2+(cf*y)**2) ics=scT/hs ifs=sfT/hs c 2-point interpolation * if(ics>=Ns)then * chi=0 * else !interpolation * chi=(ics+1-scT/hs)*Vcs(ics)+(scT/hs-ics)*Vcs(ics+1) * endif * if(ifs=Ns)then chincT=0 !beyond maximum impact parameter Else if(ics==0.or.scT/hs-ics>0.5d0) ics=ics+1 ps=scT/hs-ics if(ics==Ns)then chincT=ps*(ps-1)/2*Vcs(ics-1)+(1-ps**2)*Vcs(ics) else chincT=ps*(ps-1)/2*Vcs(ics-1)+(1-ps**2)*Vcs(ics) a +ps*(ps+1)/2*Vcs(ics+1) endif Endif c Interpolation of fT phase If(ifs>=Ns)then chinfT=0 !beyond maximum impact parameter Else if(ifs==0.or.sfT/hs-ifs>0.5d0) ifs=ifs+1 ps=sfT/hs-ifs if(ifs==Ns)then chinfT=ps*(ps-1)/2*Vfs(ifs-1)+(1-ps**2)*Vfs(ifs) else chinfT=ps*(ps-1)/2*Vfs(ifs-1)+(1-ps**2)*Vfs(ifs) a +ps*(ps+1)/2*Vfs(ifs+1) endif Endif ! Old version where Vcs and Vfs contain only the integral of the potential * chi=-chi/v0 * echin=exp(ci*chi) !nulear part of the phase ! New version where Vcs and Vfs contain also the -i/hbar/v factor ! this part is now below, where the Coulomb phase is evaluated * echin=exp(chi) * if(ib==200)write(*,*)b,chi,echin ! to include computation of stripping cross section if(lec==7)then !computes stripping cross section echin=abs(exp(chincT))**2*(1-abs(exp(chinfT))**2) else echin=exp(chincT+chinfT) endif ENDIF !npotN c Calculation of Coulomb eikonal phase c Exact Coulomb phase (pt-pt) chicct=1+2*cc*x/b+cc**2*(x**2+y**2)/b**2 chicct=e2ct/v0*log(chicct) chicft=1+2*cf*x/b+cf**2*(x**2+y**2)/b**2 chicft=e2ft/v0*log(chicft) c Coulomb with approximation of the logarithm * chicct=qc*x/b * chicft=qf*x/b c Coulomb (analytical integration up to Zm) (a=-monopole term) * Zm=10*Nz*hz * chicct=Zm/sqrt((b+cc*x)**2+cc**2*y**2) * chicct=-2*e2ct/v0*log(chicct+sqrt(chicct**2+1)) * a +2*e2ct/v0*log(Zm/b+sqrt((Zm/b)**2+1)) * chicft=Zm/sqrt((b+cf*x)**2+cf**2*y**2) * chicft=-2*e2ft/v0*log(chicft+sqrt(chicft**2+1)) * a +2*e2ft/v0*log(Zm/b+sqrt((Zm/b)**2+1)) c Coulomb V(RxT)-V(R) integrated up to Zm (in a different way) * Zm=10*Nz*hz * chicct=-2*e2ct/v0*log(b/sqrt((b+cc*x)**2+cc**2*y**2)) * a -2*e2ct/v0*log((Zm+sqrt(Zm**2+(b+cc*x)**2+cc**2*y**2)) * b /(Zm+sqrt(Zm**2+b**2))) * chicft=-2*e2ft/v0*log(b/sqrt((b+cf*x)**2+cf**2*y**2)) * a -2*e2ft/v0*log((Zm+sqrt(Zm**2+(b+cf*x)**2+cf**2*y**2)) * b /(Zm+sqrt(Zm**2+b**2))) chicPT=chicft+chicct *** c echic=exp(ci*(chicct+chicft)) ** echic=exp(ci*chicct)-ci*chicct * echic=ci*(chicct-chicft) !Coulomb phase * echic=exp(echic)-echic !Subtraction of the 1st-order c Coulomb phase and subtraction of the 1st-order of the eikonal If(lec<10)then !Usual Eikonal => no correction if(lec==7)then !computes stripping cross section echic=1 else echic=exp(ci*chicPT) endif c test without nuclear and/or Coulomb interaction * echin=1 * echic=1 echi(1)=echin*echic !total phase echi(2)=0 echi(3)=0 ElseIf(lec<20)then !Modified Eikonal echic=exp(ci*chicPT)-ci*chicPT * echic=exp(ci*chicPT) !sans correction !supprimer cor2 et cor3 * echic=exp(-ci*chicPT)+ci*chicPT !signes de Suzuki (ntypo=0) * echic=exp(ci*chicPT) !test pr verifier signes c test without nuclear and/or Coulomb interaction * echin=1 * echic=1 echi(1)=echin*echic !total phase echi(2)=echin*x !part to be corrected with K1 echi(3)=echin*z !part to be corrected with K0 * ElseIf(lec<30)then !Modified Eikonal with correction scT>str * if(scT<=str)then ! no correction * echic=exp(ci*chicPT) *c test without nuclear and/or Coulomb interaction ** echin=1 ** echic=1 * echi(1)=echin*echic !total phase * echi(2)=0 * echi(3)=0 * else * echic=exp(ci*chicPT)-ci*chicPT ** echic=exp(ci*chicPT) !sans correction !supprimer cor2 et cor3 ** echic=exp(-ci*chicPT)+ci*chicPT !signes de Suzuki (ntypo=0) ** echic=exp(ci*chicPT) !test pr verifier signes *c test without nuclear and/or Coulomb interaction ** echin=1 ** echic=1 * echi(1)=echin*echic !total phase * echi(2)=echin*x !part to be corrected with K1 * echi(3)=echin*z !part to be corrected with K0 * endif ElseIf(lec<110)then !First-order approximation c ie calculation without nuclear and/or Coulomb interaction echi(1)=1 !total phase echi(2)=x !part to be corrected with K1 echi(3)=z !part to be corrected with K0 EndIf * write(*,*)echi c test without any interaction * echi(1)=1 * echi(2)=0 * echi(3)=0 c Multipolar expansion of the eikonal phase iy=0 iye=0 iyo=0 do la=0,lamax !projection on spherical harmonics do mu=0,la iy=iy+1 !=(la,mu) if(mod(la+mu,2)==0)then iye=iye+1 * if(ib==1.and.ip==1.and.it==1)write(*,*)la,mu,iye feik(iye,ir,ib,1:2)=feik(iye,ir,ib,1:2) a +wt(it)*ylm(iy,it) b *cos(mu*phi)*echi(1:2)*8*pi/Np else iyo=iyo+1 * if(ib==1.and.ip==1.and.it==1)write(*,*)la,mu,iyo feik(iyo,ir,ib,3)=feik(iyo,ir,ib,3) a +wt(it)*ylm(iy,it) b *cos(mu*phi)*echi(3)*8*pi/Np endif enddo enddo 6 continue !ip 5 continue !it c feik(:,ir,ib,:)=f(:,:)*8*pi/Np !integral on the whole sphere feik(1,ir,ib,1)=feik(1,ir,ib,1)-sqrt(4*pi) c subtract the (unsignificant) part which'd correspond to echi=1 4 continue !ir 3 continue !ib deallocate(ylm,f,p,yl) return 1000 format(1h ,10f12.6) end *calculFBU subroutine calculFBU USE cstes USE reseaueik USE collision USE boundstates USE phaseik USE Sbu implicit real*8 (a-h,o-z) * real*8 :: fBessJ(30),fBessY(30) real*8, allocatable, dimension(:):: fBessJ,fBessY * complex*16 :: ctmp,z,cBessJ(30) !old version using S17DEF complex*16 :: ctmp fbu=0 * hbpol=0.05d0 * Nbesmax=(2*Nl+iS+Jel(1))/2 !old version with S17DEF Nbesmax=(2*Nl+iS+Jel(1))/2+1 !version using BESSJ from PHYSTM Allocate(fBessJ(Nbesmax+1),fBessY(Nbesmax+1)) do iq=1,Nq q=2*EKPT*sin(theta(iq)/2) do ib=1,Nb-1 * bleft=b0+(ib-1)*hb bleft=xb(ib) * Npol=hb/hbpol Npol=nint((xb(ib+1)-xb(ib))/hbpol) if(abs(Npol*1.0d0-(xb(ib+1)-xb(ib))/hbpol)>1d-3)then write(*,*)'Error in interpolation at ib=',ib write(*,*)'hb=',xb(ib+1)-xb(ib) write(*,*)'hbpol=',hbpol STOP endif do ipol=1,Npol bpol=bleft+(ipol-1)*hbpol if(abs(bpol)<1d-3)then chic2=0 else * tmpchic=2.0d0*bpol/abs(v0)/(Tout-Tin) * chic2=2*etaa*log((1+sqrt(1+tmpchi**2))/tmpchi) c The Coulomb scattering phase is chosen equal to the the c Coulomb eikonal phase that reproduces the exact Coulomb amplitude c cf. Bertulani & Danielewicz pg 423 Eq.(8.77). chic2=2*etaPT*log(EKPT*bpol) endif * z=q*bpol * call S17DEF(0,z,Nbesmax,'U',cBessJ,NZ,IFAIL) x=q*bpol call BESSJ(Nbesmax,x,fBessJ,fBessY) DO ief=1,Nsed iM0=-Jel(1)+(ief-1)*2 jlJM=0 Do l=0,Nl Do iJ=abs(2*l-iS),2*l+iS,2 Do iM=-iJ,iJ,2 jlJM=jlJM+1 ibess=abs(iM-iM0)/2+1 * Bessel=real(cBessJ(ibess)) !version w/S17DEF do iE=1,NE c Two point interpolation of Skljm ctmp=(SklJM(jlJM,iE,ief,ib+1) a -SklJM(jlJM,iE,ief,ib))*(ipol-1)/Npol b +SklJM(jlJM,iE,ief,ib) c Integral over b * fbu(iq,ief,jlJM,iE)= * a fbu(iq,ief,jlJM,iE) * b +Bessel*ctmp*bpol*hbpol !version w/S17DEF * c *exp(ci*chic2) fbu(iq,ief,jlJM,iE)= a fbu(iq,ief,jlJM,iE) b +fBessJ(ibess)*ctmp*bpol*hbpol c *exp(ci*chic2) enddo !iE EndDo !iM EndDo !iJ EndDo !l ENDDO !ief enddo !ipol enddo !ib enddo !iq Deallocate(fBessJ,fBessY) end *vJ function vJ(r,l,iJ) c Effective potential between the core and the fragment in the l J channel c It contains nuclear (Woods-Saxon or Gaussian), Coulomb (point-sphere) and centrifugal USE cstes USE collision USE potentials implicit real*8 (a-h,o-z) jls=iJ*(iJ+2)-4*l*(l+1)-iS*(iS+2) IF(mod(ntypo,10)==1)then !Different potentials for various partial waves ipo=1 do if(l==lpo(ipo).and.iJ==Jpo(ipo).or.ipo==Njpo)then c special channels or all others exit else ipo=ipo+1 endif enddo ELSEIF(mod(ntypo,10)==2)then !Different potentials for odd and even waves ipo=2-mod(l,2) ELSE write(*,*)'Uprogrammed potential' STOP ENDIF IF(ntypo<10)then !Woods-Saxon form eve=exp(-(r-rp(ipo))/ap(ipo)) vJ=Vp(ipo)*eve/(1+eve)+l*(l+1)/r**2/deuxm if(Vls(ipo)/=0)then eve=exp(-(r-rls(ipo))/als(ipo)) vJ=vJ-jls*Vls(ipo)*eve/(1+eve)**2/als(ipo)/r/8 endif ELSEIF(ntypo<20)then !Gaussian form eve=(r-rp(ipo))/ap(ipo) eve=exp(-eve**2/2) vJ=Vp(ipo)*eve+l*(l+1)/r**2/deuxm if(Vls(ipo)/=0)then !the ls term is used for the derivative eve=(r-rls(ipo))/als(ipo) eve=(r-rls(ipo))**2*exp(-eve**2/2) !(r-r0)^2*Gaussian form vJ=vJ+Vls(ipo)*eve endif ELSE write(*,*)'Unprogrammed potential' STOP ENDIF If (kZf/=0) then c Point-sphere Coulomb potential if(r>=rC(ipo))then c1=kZf*kZc*e2/r vJ=vJ+c1 else c1=kZf*kZc*e2/2/rC(ipo) c2=3-(r/rC(ipo))**2 vJ=vJ+c1*c2 endif Endif return end *vSW function vSW(V,r0,a,r) c Computes the value of a Woods-Saxon potential in r implicit none real*8 vSW,V,r0,a,r,eve eve=exp((r0-r)/a) vSW=V*eve/(1+eve) return end *dvSW function dvSW(V,r0,a,r) c Computes - the first derivative of a Woods-Saxon potential in r implicit none real*8 dvSW,V,r0,a,r,eve eve=exp((r0-r)/a) dvSW=V*eve/(1+eve)**2 return end *vopt function vopt(V,r0R,aR,W,r0I,aI,WD,r0D,aD,kZ1,kZ2,rC,r) c Computes the optical (complex) potential at r. c Woods-Saxon form factor is assumed, and the Coulomb potential is point-sphere USE cstes USE collision implicit real*8 (a-h,o-z) complex*16 vopt vopt=vSW(V,r0R,aR,r)+ci*vSW(W,r0I,aI,r) If(WD/=0) vopt=vopt+ci*dvSW(WD,r0D,aD,r) If(kZ1/=0.and.kZ2/=0)then c Coulomb point-sphere potential if(r>=rC)then c1=kZ1*kZ2*e2/r vopt=vopt+c1 else c1=kZ1*kZ2*e2/2/rC c2=3-(r/rC)**2 vopt=vopt+c1*c2 endif Endif return end *voptN function voptN(V,r0R,aR,W,r0I,aI,WD,r0D,aD,kZ1,kZ2,rC,r) c Computes the nuclear part of the optical (complex) potential at r. c Woods-Saxon form factor is assumed, and the Coulomb potential is point-sphere c The pure Coulomb is subtracted USE cstes USE collision implicit real*8 (a-h,o-z) complex*16 voptN voptN=vSW(V,r0R,aR,r)+ci*vSW(W,r0I,aI,r) If(WD/=0) voptN=voptN+ci*dvSW(WD,r0D,aD,r) If(kZ1/=0.and.kZ2/=0)then c Coulomb point-sphere potential minus pure Coulomb if(r=0.and.abs(1-e/e0)