C==================================================================== C PROGRAM RELEX C C Bertulani multiple Coulomb+nuclear relativistic excitation program C C Solves Coupled-Channels equations for the relativistic C Coulomb+nuclear excitation of E0 (only nuclear), E1, E2 and M1 modes. C C.A. Bertulani - version 13/Jun/98 C C-------------------------------------------------------------------- C HOW TO MAKE AN INPUT C-------------------------------------------------------------------- C INPUT and OUTPUT: C Inputs are in files RELEX.DIM, RELEX.IN and OPTW.IN (input of C optical potential: optional). C Units are in MeV and fm. Output cross sections in milibarns. C Output files: RELEX.OUT --> Probabilities and cross sections C RELEX2.OUT --> Statistical tensors C-------------------------------------------------------------------- C C RELEX.DIM (dimension of vectors and matrices) C --------- C NMAX = maximum number of levels. C NBMAX = maximum # of impact parameters. C NSTMAX = maximum # of total magnetic substates. C NGRID= maximum # of points in optical potential grid (EVEN). C C RELEX.IN (unit 10) C --------- C FIRST ROW: ap,zp,at,zt,eca,iw,iout C input of system charges and masses (AP,ZP,AT,ZT), bombarding energy C per nucleon (ECA) in MeV. Option (IW): IW=0 (1) for projectile (target) C excitation. IOUT=(0)1 for (no) output of statistical tensors. C C SECOND ROW: nb,accur,bmin,itot C NB=number of points in impact parameter mesh. C ACCUR=accuracy required for time integration at each impact parameter. C BMIN= minimum impact parameter for integration of cross section C ( if BMIN=0 integration starts B = 1.2(AP^1/3+AT^1/3)/2 ). C ITOT=1(0) prints (does not print) out impact parameter probabilities. C C THIRD ROW: iopw,iopnuc C IOPW is an option to enter (or not) the optical potential: IOPW=1(0). C If optical potential is provided, it should be stored in OPTW.IN in C rows of R x Real[U(R)] x Imag[U(R)]. C The first line in this file is the number of rows (maximum=NGRID). C The program makes an interpolation for intermediate values. C IOPNUC=1(0) is an option to compute (or not) nuclear excitation. C C FOURTH ROW: nst C NST is the number of nuclear levels. C C FIFTH, AND FOLLOWING ROWS: i, ex(i), spin(i) C Input of state label (I), energy (EX), and spin (SPIN). C I ranges from 1 to NST. C C FOLLOWING ROWS: STATE1,STATE2,MATE1,MATE2,MATM1 C Reduced matrix elements for E1, E2 and M1 excitations: C , j > i , C for the electromagnetic transitions: MATE1(i -> j) (electric dipole), C MATE2(i -> j) (ele. quadrupole), and MATM1(i -> j) (magnetic dip.). C To end reading add a row of zeros. C C Notation of Edmonds, "Angular Momentum in Quantum Mechanics", C Princeton University Press (1957), is used. C C FOLLOWING ROWS: j,delte(0,j),delte(1,j),delte(2,j) C If IOPNUC=1 enter values of deformation parameters for monopole, dipole, C and quadrupole nuclear excitations (DELTE0,DELTE1,DELTE2) for each excited C state J. C _______ C EXCITATION OF GIANT RESONANCES: C A good guess for the electric reduced matrix elements is obtained from C C ||^2 = (2I_i+1) . C 9/(4 pi) . hbar^2/(2 m_N) . e^2 . NZ/A / E_x C C where E_x = E_j - E_i . For E2 transitions C C ||^2 = (2I_i+1) . C 15/(4 pi) . hbar^2/m_N . e^2 . R^2 . Z^2/A / E_x (for isoscalar) C 15/(4 pi) . hbar^2/m_N . e^2 . R^2 . NZ/A / E_x (for isovector) C C The excitation of GR's exhaust a fraction of these sum rules, which is C close to unity. C C For the nuclear deformation parameters, the sum rules are C Delte0^2 = 2 pi . hbar^2/m_N / ( A E_x) C Delte1^2 = chi . pi / 2 . hbar^2/m_N . A / (NZ E_x) C Delte2^2 = 20 pi / 3 . hbar^2/m_N / (A E_x) C where chi= 2 (R_n - R_p)/(R_n + R_p) is the neutron skin fraction. C C-------------------------------------------------------------------- C START PROGRAM C-------------------------------------------------------------------- C implicit real*8(a-h,o-z) include 'relex.dim' integer*4 intr(0:ngrid),mj(nstmax),ii(0:nstmax),istop(0:nmax) real*8 mm(nstmax),jj(nstmax),intb(nbmax),b(nbmax),mnkap, & mat(3,0:nmax,0:nmax),spin(nmax),fd(2),delte(0:2,nmax), & exx(nstmax),pcc(nmax,nbmax),ex(nmax),dens1(0:ngrid), & dens2(0:ngrid),dru(0:ngrid),r(0:ngrid),z(0:ngrid),tb(nbmax), & ru(0:ngrid),psitot(3,nstmax,nstmax),psi(3,nbmax,nstmax,nstmax), & xi(nstmax,nstmax),spinave(nmax,nbmax),pint(nbmax),pint_2(nbmax), & bcorr(nbmax) complex*16 ei,ca(nstmax),dca(nstmax),wn(0:ngrid),wn2(0:ngrid), & phasen(nbmax),psinuc(0:2,0:ngrid,nstmax,nstmax),wn1(0:ngrid), & fact,u(0:ngrid),ubm(0:2,0:ngrid),stat(nmax,0:10,-10:10) real*8 iv(0:130),fak(0:130),fad(0:130) common /gf1 / iv common /gf2 / fak common /gf3 / fad common /a1/ beta,gamma,ei,pi,hc common/a3/xi common/a4/mm,n common/a5/psitot common/a6/psinuc common/a7/bb common/a8/iopnuc common/a9/delr common/int1/ intr C external dcadt,rkqs C C initialize factorial routine call gfv(30) C open(10,file='relex.in', status='old',readonly) open(11,file='optw.in', status='old',readonly) open(12,file='relex.out', status='unknown') open(15,file='relex2.out', status='unknown') C read(10,*)ap,zp,at,zt,eca,iw,iout read(10,*)nb,accur,bmin,itot read(10,*)iopw,iopnuc if(nb.ge.nbmax) then write(12,*)' WARNING - number of impact parameters is too big' write(12,*)' Decrease value of NB, or increase NBMAX' go to 2000 endif C write(12,*)' ' write(12,'('' **** Output of RELEX ****'')') write(12,*)' ' write(12,109) 109 format(' Input parameters:',/) write(12,1)ap,zp,at,zt,eca,iw,ngrid,nb,accur,iopw,iopnuc,iout C 1 format(2x,' Nuclear mass and charge numbers: ',4(F4.0,2x),/, 1 2x, ' Laboratory energy/per nucleon [MeV]: ',F6.0,/,2x, 2 ' Option for the excited nucleus [0(1) for proj (targ)] : ', 3 I2,/,2x,' Integration parameters (NGRID,NB,ACCUR):',I4,2x, 4 I4,2x,F9.4,/,2x,' Enter optical potential 0(no) 1(yes):',I2,/, 5 2x,' Compute nuclear excitation 0(no) 1(yes):',I2,/, 5 2x,' Output of statistical 0(no) 1(yes):',I2,/, 6 I5,/) C read(10,*)nst if(nst.gt.10)goto 1500 write(12,2)nst 2 format(2x,' Number of nuclear states = ',I2,//, 1 ' State E [MeV] Spin ') C do i=1,nst read(10,*)j,ex(j),spin(j) write(12,3)j,ex(j),spin(j) 3 format(2x,I3,3x,F10.3,4x,F10.3) enddo C write(12,4) 4 format(/,' Red. matrix elem. (E1, M1: fm.e), (E2: fm^2.e)' 1 ,//,2x,16x,'E1',10x,'E2',10x,'M1') C C read matrices (E/M;L, I_i-->I_j) values do k=1,nstmax read(10,*)i,j,MAT(1,i,j),MAT(2,i,j),MAT(3,i,j) if(i.eq.0)goto 51 write(12,5)i,j,mat(1,i,j),mat(2,i,j),mat(3,i,j) 5 format(2x,i2,' --> ',i2,2X,3(f10.3,2x)) enddo C C construct magnetic quantum numbers and catalogue states 51 j=1 ii(0)=0 istop(0)=0 do i=1,nst mj(i)=nint(2.*spin(i)+1.) do m=0,mj(i)-1 ii(j)=i jj(j)=spin(i) mm(j)=-spin(i)+dfloat(m) exx(j)=ex(i) j=j+1 enddo istop(i)=istop(i-1)+mj(i) enddo n=j-1 if(n.gt.nstmax) then write(12,*)' WARNING - number of states is too big' write(12,*)' Increase parameter NSTMAX in file RELEX.DIM' go to 2000 endif C C Index T is reserved for the excited nucleus (depends on option IW) if(iw.eq.0) then dum1=ap dum2=zp ap=at zp=zt at=dum1 zt=dum2 endif C C constants (signn=nucleon-nucleon cross section) ei=cmplx(0.,1.) pi=acos(-1.) e2=1.44 amu=931.5 hc=197.33 signn=4. C C nuclear deformation parameters if(iopnuc.eq.1) then write(12,6) 6 format(/,' Deform. Par. DELTE0 (dimensionless),', 1 ' DELTE1^2 (fm), and ','DELTE2^2 (fm) ' 2 ,//,2x,' Transit.',7x,'DELTE0',6x,'DELTE1',6x, 3 'DELTE2') do j=2,nst read(10,*)idum,delte(0,j),delte(1,j),delte(2,j) write(12,7)idum,delte(0,j),delte(1,j),delte(2,j) 7 format(2x,i2,9X,3(f10.3,2x)) enddo endif C C Lorentz variables gamma=1.+eca/amu beta=sqrt((gamma-1.)*(gamma+1.))/gamma C C grazing impact parameter rp13=1.2*ap**(1./3.) rt13=1.2*at**(1./3.) b13=rp13+rt13 C C mesh in impact parameter and integration factors bmin1=b13/2. bmax1=1.5*b13 nb1=int(nb/2) if(iv(nb1).ge.0.)nb1=nb1+1 delb1=(bmax1-bmin1)/(nb1-1) do ib=1,nb1 b(ib)=bmin1+(ib-1)*delb1 enddo nb2=nb-nb1 bmax=200. delb2=(bmax-bmax1)/nb2 do ib=nb1+1,nb b(ib)=bmax1+(ib-nb1)*delb2 enddo C C optical potential (rmax=maximum range of the potential) rmax=30. delr=rmax/ngrid do i=0,ngrid r(i)=i*delr z(i)=i*delr enddo C C factors for Simpson's integration in nuclear interactions intr(0)=1 intr(ngrid)=1 ig=4 do i=1,ngrid-1 intr(i)=ig ig=6-ig enddo C C read optical potential from OPTW.IN if(iopw.eq.0)goto 81 read(11,*)nr ru(0)=0. do i=1,nr read(11,*)ru(i),vdum1,vdum2 u(i)=vdum1+ei*vdum2 dru(i)=ru(i)-ru(i-1) enddo u(0)=u(1) C C interpolate optical potential ir=0 ncount=0 8 ir=ir+1 if(ir.gt.nr-1)goto 81 do i=ncount,ngrid wn(i)=0. aux=ru(ir)-r(i) if(aux.ge.0.) then al=1.-aux/dru(ir) wn(i)=al*u(ir)+(1.-al)*u(ir-1) ncount=ncount+1 else goto 8 endif enddo C 81 if(iopnuc.eq.0)goto 9 C C Nuclear interaction in the Bohr-Mottelson model call derivative(wn,wn1,wn2,delr,ngrid) do i=0,ngrid ubm(0,i)=3.*wn(i)+r(i)*wn1(i) ubm(1,i)=1.5*(wn1(i)+r(i)*wn2(i)/3.) ubm(2,i)=wn1(i) enddo C 9 if(iopw.eq.1)goto 10 C C if optical potential is not known: compute liquid-drop densities do i=0,ngrid dens1(i)=rhonp(r(i),ap,zp)+rhopp(r(i),ap,zp) dens2(i)=rhonp(r(i),at,zt)+rhopp(r(i),at,zt) enddo fact=-ei*signn*hc*beta/2. call twofold(fact,r,rmax,delr,dens1,dens2,wn) C C nuclear eikonal phase and absorption factor 10 call phnuc(wn,b,nb,z,phasen) do ib=1,nb tb(ib)=exp(-2.*dimag(phasen(ib))) enddo C C computation of psi matrix and nuclear excitation matrices C notation: 1=E1, 2=E2, 3=M1 fac=zp*e2/beta/hc fd(1)=fac*sqrt(2.*pi/3.) fd(2)=-fac*sqrt(2.*pi/15.)/2. do i=1,n-1 do j=i+1,n mu=nint(mm(j)-mm(i)) do l=1,3 psi(l,ib,i,j)=0. ll=l if(l.eq.3)ll=1 bmat=mat(l,ii(i),ii(j)) iv1=iv(nint(jj(j)+dfloat(ll)-mm(j)+1)) three= threej(jj(j),dfloat(ll),jj(i),-mm(j),dfloat(mu),mm(i)) aux=bmat*iv1*three do ib=1,nb psi(l,ib,i,j)=fd(ll)*aux/b(ib)**ll enddo enddo c nuclear part do ln=0,2 bmatn=delte(ln,ii(j)) iv2=iv(nint(jj(j)-mm(j))) zero=0.d0 three= threej(jj(j),dfloat(ln),jj(i),-mm(j),dfloat(mu),mm(i)) three0=threej(jj(j),dfloat(ln),jj(i),zero,zero,zero) auxn=bmatn*sqrt((2.*jj(i)+1.)/ & (2.*jj(j)+1.)/4./pi)*iv2*three*three0*gamma do ir=0,ngrid psinuc(ln,ir,i,j)=psinuc(ln,ir,i,j)+auxn*ubm(ln,ir) enddo enddo enddo enddo write(15,110) C 110 format(1H1,9x, & 49HThe angular distribution tensors STAT(J,KA,KAPPA)) write(15,111) 111 format(7x,1HJ,5x,2HKA,4x,5HKAPPA,9x,9HReal STAT,11x,9HImag STAT) c C ========================== C start integration in time for each impact parameter C C impact parameter loop icc=0 do ib=1,nb e0=hc*beta*gamma/b(ib) bb=b(ib) do i=1,n do j=1,n xi(i,j)=(exx(i)-exx(j))/e0 if(j.lt.nst)spinave(j,ib)=0. do l=1,3 psitot(l,i,j)=psi(l,ib,i,j) enddo enddo enddo C initialize statistical tensors do ka=0,4,2 do kappa=-ka,ka do j=2,nst stat(j,ka,kappa)=0. enddo enddo enddo C C initialize amplitudes and average over initial polarization mj1=1 1001 do k=1,n ca(k)=0. dca(k)=0. enddo ca(mj1)=cmplx(1.,0.) C C ************** integrate coupled-channels equations C C initial guess for the time mesh tmax=15. nt=50 dtau=2.*tmax/nt t1=-tmax t2=tmax nvar=n h1=dtau hmin=0. C accuracy eps=accur C Numerical Recipes routine for adaptive Runge-Kutta integration C call odeint(ca,nvar,t1,t2,eps,h1,hmin,nok,nbad,dcadt,rkqs) C C ************ end time integration C C excitation probabilities averaged over initial spin do j=1,nst sum=0. do i=istop(j-1)+1,istop(j) sum=sum+abs(ca(i))**2 enddo pcc(j,ib)=sum*tb(ib) enddo C C check if norm is conserved xn=0. do j=1,nst xn=xn+pcc(j,ib) enddo xn=xn+1.-tb(ib) eps10=eps*10. if(abs(xn-1.).gt.eps10) then write(15,*)' Error in probability exceeds 10 x ACCUR' go to 2000 endif C C sum excitation probabilities over initial spin do j=1,nst spinave(j,ib)=spinave(j,ib)+pcc(j,ib)/mj(1) enddo if(iout.eq.0)goto 112 C C statistical tensors do ka=0,4,2 do kappa=-ka,ka do i=1,nst jini=istop(ii(i-1))+1 jfin=istop(ii(i)) do j=jini,jfin mnkap=mm(j)+kappa iv3=iv(nint(abs(jj(j)+mnkap))) three= threej(jj(j),jj(j),dfloat(ka),-mnkap, & mm(j),dfloat(kappa)) jkap=j+kappa if(jkap.ge.jini.and.jkap.le.jfin) then stat(i,ka,kappa)=stat(i,ka,kappa)+iv3*three* & sqrt(2.*jj(j)+1.)*conjg(ca(j))*ca(jkap)/mj(1) endif enddo enddo enddo enddo C 112 mj1=mj1+1 if(mj1.le.mj(1))goto 1001 C end sum over intial orientation C C output of statistical tensors if(iout.eq.1) then write(15,*)'b=',b(ib) do j=2,nst do ka=0,4,2 do kappa=0,ka if(stat(j,ka,kappa).ne.0.) then write(15,113)j,ka,kappa,stat(j,ka,kappa) endif 113 format(1x,3I7,2E22.4) enddo enddo enddo endif C enddo C C end impact parameter loop C ========================= C if(itot.eq.1) then write(12,116) 116 format(//, & ' Impact parameter, state, and occupation probabilities',/) do ib=1,nb do j=2,nst if(j.eq.2) then write(12,118)b(ib),j,spinave(j,ib) else write(12,119)j,spinave(j,ib) endif 118 format(2x,g12.4,i2,2x,g12.4) 119 format(14x,i2,2x,g12.4) enddo enddo endif c c **** Compute total cross sections write(12,160),bmin 160 format(/,' Cross sections in mb, with BMIN=',f9.4) c factor for recoil correction ared=ap*at/(ap+at) aclose=zp*zt*e2/(ared*amu*beta*beta) bc=pi*aclose/2./gamma C do k=2,nst sum1=0. C correct for Coulomb recoil do ib=1,nb bcorr(ib)=b(ib)-bc pint(ib)=20.*pi*b(ib)*spinave(k,ib) enddo C interpolate and integrate call spline(bcorr,pint,nb,1.d30,1.d30,pint_2) call qsimp(bcorr,pint,pint_2,nb,max(bmin1,bmin),bmax,cross) C output WRITE(12,'('' State '',i3,'' = '','// & 'g12.4)')k,cross enddo C **** end cross section C goto 2000 1500 write(12,1501) 1501 format(26h0ERROR: NST EXCEEDS NSTMAX) 2000 stop end C-------------------------------------------------------------------- C END MAIN PROGRAM C-------------------------------------------------------------------- C C==================================================================== C subroutine dcadt(tau,ca,dca) C C time derivative of the occupation amplitudes C include 'relex.dim' implicit real*8(a-h,o-z) complex*16 ei,ca(nstmax),dca(nstmax),expt,vmat(nstmax,nstmax) real*8 xi(nstmax,nstmax),mm(nstmax) common /a1/ beta,gamma,ei,pi,hc common/a3/xi common/a4/mm,n call vint(tau,vmat) do i=1,n dca(i)=0. do j=1,n expt=exp(ei*tau*xi(j,i)) dca(i)=dca(i)-ei*vmat(i,j)*expt*ca(j) enddo enddo return end C C===================================================================== C subroutine vint(t,vmat) C C time-dependent fields for E1, E2 and M1 transitions C implicit real*8(a-h,o-z) include 'relex.dim' complex*16 ei,q,ve1,ve2,vm1,vmat(nstmax,nstmax),Ylm, & psinuc(0:2,0:ngrid,nstmax,nstmax) real*8 mm(nstmax) real*8 psitot(3,nstmax,nstmax),xi(nstmax,nstmax) common /a1/ beta,gamma,ei,pi,hc common/a3/xi common/a4/mm,n common/a5/psitot common/a6/psinuc common/a7/bb common/a8/iopnuc common/a9/delr t2=t*t phi=1./sqrt(1.+t2) phi3=phi**3 phi5=phi**5 beta2=beta*beta gamma2=gamma*gamma do i=1,n-1 do j=i+1,n mu=nint(mm(j)-mm(i)) vmat(i,j)=0. x=xi(j,i) C C electric dipole t.d. field if(mu.eq.0) then q=sqrt(2.)*gamma*(t*phi3-ei*x*beta2*phi) else q=-mu*phi3 endif ve1=psitot(1,i,j)*q C C electric quadrupole t.d. field if(iabs(mu).eq.2) then q=3.*phi5 else if(iabs(mu).eq.1) then q=mu*gamma*(6.*t*phi5-ei*x*beta2*phi3) else q=sqrt(6.)*gamma2*((2.*t2-1.)*phi5- & ei*beta2*t*phi3) endif endif ve2=psitot(2,i,j)*q C C magnetic dipole t.d. field q=iabs(mu)*ei*beta*phi3 vm1=psitot(3,i,j)*q C C add E1, E2 and M1 vmat(i,j)=ve1+ve2+vm1 if(iopnuc.ne.1)goto 20 c goto 20 C C add nuclear interaction rr=bb/phi ir=int(rr/delr) h=1.-(rr-delr*ir)/delr theta=t*phi zero=0.d0 vnuc=0. if(ir+1.lt.ngrid.and.ir.gt.0.) then do ln=0,2 if(abs(mu).le.ln)vnuc=vnuc+(h*psinuc(ln,ir,i,j)+ & (1.-h)*psinuc(ln,ir+1,i,j))*Ylm(ln,mu,theta,zero) enddo vmat(i,j)=vmat(i,j)+vnuc endif C C interaction is hermitian 20 vmat(j,i)=conjg(vmat(i,j)) enddo vmat(i,i)=0. enddo return end C C===================================================================== C real*8 FUNCTION rhopp(rg,ap,zp) C C Droplet densities for protons - See Myers and Swiatecki, Ann. of C Phys. A55 (1969) 395 and C Ann. of Phys. A84 (1974) 186 implicit real*8 (a-h,o-z) pi=3.141597265 tpp = 2.4 bpp = 0.413 * tpp deltap = ((ap-2.*zp)/ap + 8.076E-3*zp*ap**(-.66666)) # / ( 1. + 4.871 * ap**(-.33333) ) epsp = -0.1724*ap**(-.33333) + 3.051E-3 * zp**2 * ap**(-1.33333) # + .4166666 * deltap**2 rp = 1.18 * ap**0.33333 * ( 1. + epsp ) dp = 0.66666 * ((ap-2.*zp)/ap - deltap ) * rp rpp = rp - ((ap-zp)/ap) * dp cpp = rpp * ( 1. - (bpp/rpp)**2 ) rho0pp = 3. * zp / ( 12.5664 * cpp**3 * # ( 1. + pi**2 * tpp**2 / (19.36 * cpp**2))) rhopp = rho0pp / ( 1. + exp((rg-cpp)/(tpp/4.4)) ) return end C C===================================================================== C real*8 FUNCTION rhonp(rg,ap,zp) C C Droplet densities for neutrons - From Myers and Swiatecki, Ann. of C Phys. A55 (1969) 395 and C Ann. of Phys. A84 (1974) 186 implicit real*8 (a-h,o-z) pi=3.141597265 tnp = 2.4 bnp = 0.413 * tnp deltap = ((ap-2.*zp)/ap + 8.076E-3*zp*ap**(-.66666)) # / ( 1. + 4.871 * ap**(-.33333) ) epsp = -0.1724*ap**(-.33333) + 3.051E-3 * zp**2 * ap**(-1.33333) # + .4166666 * deltap**2 rp = 1.18 * ap**0.33333 * ( 1. + epsp ) dp = 0.66666 * ((ap-2.*zp)/ap - deltap ) * rp rnp = rp + (zp/ap) * dp cnp = rnp * ( 1. - (bnp/rnp)**2 ) rho0np = 3. * ( ap - zp )/ ( 12.5664 * cnp**3 * # ( 1. + pi**2 * tnp**2 / (19.36 * cnp**2))) rhonp = rho0np / ( 1. + exp((rg-cnp)/(tnp/4.4)) ) return end C C===================================================================== C SUBROUTINE phnuc(wn,b,nb,z,phasen) C C nuclear eikonal phase C implicit real*8 (a-h,o-z) include 'relex.dim' real*8 b(nbmax),z(0:ngrid) integer*4 intr(0:ngrid) complex*16 ei,phasen(nbmax),vnucl,wn(0:ngrid) common /a1/ beta,gamma,ei,pi,hc common/a9/ delr common /int1/ intr do k=1,nb phasen(k)=0. C integrate in z coordinate - Simpson's rule do i=0,ngrid vnuc=0. rv=sqrt(b(k)**2+z(i)**2) l=int(rv/delr) al=1.-(rv-delr*l)/delr if(l+1.lt.ngrid.and.l.gt.0)vnucl=al*wn(l)+(1.-al)*wn(l+1) phasen(k)=phasen(k)+intr(i)*vnucl enddo phasen(k)=-2./hc/beta*delr/3.*phasen(k) enddo return end C C===================================================================== C SUBROUTINE twofold(fact,r,rmax,delr,dens1,dens2,v) C C Folding of densities C implicit real*8 (a-h,o-z) include 'relex.dim' real*8 dens1(0:ngrid),dens2(0:ngrid),r(0:ngrid),x(0:ngrid), & rint(0:ngrid) complex*16 v(0:ngrid),fact integer*4 intr(0:ngrid) common /int1/ intr pi=3.141597265 delx=2./ngrid do i=0,ngrid x(i)=-1.+i*delx rint(i)=r(i) enddo do i=0,ngrid sum2=0. do j=0,ngrid sum1=0. do k=0,ngrid aux=sqrt(abs(r(i)**2+rint(j)**2-2.*r(i)*rint(j)*x(k))) ifrac=int(aux/delr) al=1.-(aux-delr*ifrac)/delr if(ifrac+1.lt.ngrid.and.ifrac.gt.0) & sum1=sum1+intr(k)*(al*dens1(ifrac)+ & (1.-al)*dens1(ifrac+1)) enddo ifrac=int(rint(j)/delr+1.) al=1.-(rint(j)-delr*ifrac)/delr if(ifrac+1.lt.ngrid.and.ifrac.gt.0) & sum2=sum2+intr(j)*rint(j)**2*(al*dens2(ifrac) & +(1.-al)*dens2(ifrac+1))*delx/3.*sum1 enddo v(i)=2.*pi*fact*sum2/3.*delr enddo return end C C===================================================================== C real*8 function threej(a1,a2,a3,b1,b2,b3) C C Returns WIGNER 3J-COEFFICIENTS (for Clebsh coeficients change C name to clebsh and delete lines where b3=-b3) C a1,a2,a3 are j1,j2,j3 and b1,b2,b3 are m1,m2,m3. C array f(n+1)=n!/a**n, where a is an arbitrary number C that cancels in the calculation of clebsh; its purpose C is to prevent under/overflows in array f C implicit real*8 (a-h,o-z) dimension f(85) save f data ifirst /1/ C if (ifirst .eq. 1) then f(1)=1.0 do 10 i=1,84 10 f(i+1)=f(i)*dfloat(i)/15. ifirst=0 endif C C delete this line to run this program for clebsh coefficients b3=-b3 C clebsh=0. threej=0. if(abs(b1)-a1 .gt. 0.01) goto 100 if(abs(b2)-a2 .gt. 0.01) goto 100 if(abs(b3)-a3 .gt. 0.01) goto 100 if(abs(b1+b2-b3).gt.0.01) goto 100 if(a3-a1-a2 .gt. 0.01) goto 100 if(abs(a1-a2)-a3 .gt. 0.01) goto 100 l1=a1+a2-a3+1.01 l2=a1+b1+1.01 l3=a2+b2+1.01 l4=a2-b2+1.01 l5=a3+a2-a1+1.01 l6=a1-b1+1.01 l7=a3+b3+1.01 l8=a3-b3+1.01 l9=a3+a1-a2+1.01 l10=a3+a1+a2+1.01 m1=l1 m2=l6 m3=l3 m4=(l7+l8-l3-l4+l2-l6)/2+1 m5=(l7+l8-l2-l6-l3+l4)/2+1 m6=1 zw=sqrt(f(l1)*f(l2)*f(l6)) tw=sqrt(f(l9)/f(l10)*(2.0*a3+1.)/float(l10)) sw=tw*sqrt(f(l3)*f(l7)*f(l8)*f(l4)*f(l5)) min=max0(-m4,-m5,-m6)+2 max=min0(m1,m2,m3) fr=-(-1.)**min if(max.lt.min) return do m=min,max mw=m-1 clebsh=clebsh+sw*fr/f(m4+mw)/(f(m2-mw)/zw)/f(m3-mw)/ 1 f(m1-mw)/f(m5+mw)/f(m6+mw) fr=-fr enddo threej=(-1)**nint(a1-a2+b3)/sqrt(2.*a3+1.)*clebsh 100 continue C C delete this line to run this program for clebsh coefficients b3=-b3 C return end C C===================================================================== C SUBROUTINE Derivative(Wf,Wf1,Wf2,Dx,Npnts) C C Calculates first (WF1) and second derivative (WF2) of function WF C Initialize Ncount=0 in calling program C implicit real*8(a-h,o-z) real*8 A(6,6),B(6,6) complex*16 Wf(0:Npnts),Wf1(0:Npnts),Wf2(0:Npnts) common/coeff/ Ncount if(Ncount.gt.0) go to 1 Ncount=1 ds=120.0d0*dx A(1,1)=-274.0d0/ds A(1,2)=600.0d0/ds A(1,3)=-600.0d0/ds A(1,4)=400.0d0/ds A(1,5)=-150.0d0/ds A(1,6)=24.0d0/ds A(2,1)=-24.0d0/ds A(2,2)=-130.0d0/ds A(2,3)=240.0d0/ds A(2,4)=-120.0d0/ds A(2,5)=40.0d0/ds A(2,6)=-6.0d0/ds A(3,1)=6.0d0/ds A(3,2)=-60.0d0/ds A(3,3)=-40.0d0/ds A(3,4)=120.0d0/ds A(3,5)=-30.0d0/ds A(3,6)=4.0d0/ds A(4,1)=-4.0d0/ds A(4,2)=30.0d0/ds A(4,3)=-120.0d0/ds A(4,4)=40.0d0/ds A(4,5)=60.0d0/ds A(4,6)=-6.0d0/ds A(5,1)=6.0d0/ds A(5,2)=-40.0d0/ds A(5,3)=120.0d0/ds A(5,4)=-240.0d0/ds A(5,5)=130.0d0/ds A(5,6)=24.0d0/ds A(6,1)=-24.0d0/ds A(6,2)=150.0d0/ds A(6,3)=-400.0d0/ds A(6,4)=600.0d0/ds A(6,5)=-600.0d0/ds A(6,6)=274.0d0/ds C ds=60.0d0*Dx**2 C B(1,1)=225.0d0/ds B(1,2)=-770.0d0/ds B(1,3)=1070.0d0/ds B(1,4)=-780.0d0/ds B(1,5)=305.0d0/ds B(1,6)=-50.0d0/ds B(2,1)=50.0d0/ds B(2,2)=-75.0d0/ds B(2,3)=-20.0d0/ds B(2,4)=70.0d0/ds B(2,5)=-30.0d0/ds B(2,6)=5.0d0/ds B(3,1)=-5.0d0/ds B(3,2)=80.0d0/ds B(3,3)=-150.0d0/ds B(3,4)=80.0d0/ds B(3,5)=-5.0d0/ds B(3,6)=0.0d0 B(4,1)=0.0d0 B(4,2)=-5.0d0/ds B(4,3)=80.0d0/ds B(4,4)=-150.0d0/ds B(4,5)=80.0d0/ds B(4,6)=-5.0d0/ds B(5,1)=5.0d0/ds B(5,2)=-30.0d0/ds B(5,3)=70.0d0/ds B(5,4)=-20.0d0/ds B(5,5)=-75.0d0/ds B(5,6)=50.0d0/ds B(6,1)=-50.0d0/ds B(6,2)=305.0d0/ds B(6,3)=-780.0d0/ds B(6,4)=1070.0d0/ds B(6,5)=-770.0d0/ds B(6,6)=225.0d0/ds 1 continue C do 2 Ir=3,NPNTS-3 WF1(Ir)=A(3,1)*Wf(Ir-2)+A(3,2)*Wf(Ir-1)+A(3,3)*Wf(Ir) & +A(3,4)*Wf(Ir+1)+A(3,5)*Wf(Ir+2)+A(3,6)*Wf(Ir+3) WF2(Ir)=B(3,1)*Wf(Ir-2)+B(3,2)*Wf(Ir-1)+B(3,3)*Wf(Ir) & +B(3,4)*Wf(Ir+1)+B(3,5)*Wf(Ir+2)+B(3,6)*Wf(Ir+3) 2 continue Ir=1 WF1(IR)=A(1,1)*Wf(Ir)+A(1,2)*Wf(Ir+1)+A(1,3)*Wf(Ir+2) & +A(1,4)*Wf(Ir+3)+A(1,5)*Wf(Ir+4)+A(1,6)*Wf(Ir+5) C WF2(IR)=B(1,1)*Wf(Ir)+B(1,2)*Wf(Ir+1)+B(1,3)*Wf(Ir+2) & +B(1,4)*Wf(Ir+3)+B(1,5)*Wf(Ir+4)+B(1,6)*Wf(Ir+5) C Ir=2 WF1(IR)=A(2,1)*Wf(Ir-1)+A(2,2)*Wf(Ir)+A(2,3)*Wf(Ir+1) & +A(2,4)*Wf(Ir+2)+A(2,5)*Wf(Ir+3)+A(2,6)*Wf(Ir+4) C WF2(IR)=B(2,1)*Wf(Ir-1)+B(2,2)*Wf(Ir)+B(2,3)*Wf(Ir+1) & +B(2,4)*Wf(Ir+2)+B(2,5)*Wf(Ir+3)+B(2,6)*Wf(Ir+4) C Ir=NPNTS-2 WF1(IR)=A(4,1)*Wf(Ir-3)+A(4,2)*Wf(Ir-2)+A(4,3)*Wf(Ir-1) & +A(4,4)*Wf(Ir)+A(4,5)*Wf(Ir+1)+A(4,6)*Wf(Ir+2) WF2(IR)=B(4,1)*Wf(Ir-3)+B(4,2)*Wf(Ir-2)+B(4,3)*Wf(Ir-1) & +B(4,4)*Wf(Ir)+B(4,5)*Wf(Ir+1)+B(4,6)*Wf(Ir+2) C Ir=NPNTS-1 WF1(IR)=A(5,1)*Wf(Ir-4)+A(5,2)*Wf(Ir-3)+A(5,3)*Wf(Ir-2) & +A(5,4)*Wf(Ir-1)+A(5,5)*Wf(Ir)+A(5,6)*Wf(Ir+1) WF2(IR)=B(5,1)*Wf(Ir-4)+B(5,2)*Wf(Ir-3)+B(5,3)*Wf(Ir-2) & +B(5,4)*Wf(Ir-1)+B(5,5)*Wf(Ir)+B(5,6)*Wf(Ir+1) C Ir=NPNTS WF1(IR)=A(6,1)*Wf(Ir-5)+A(6,2)*Wf(Ir-4)+A(6,3)*Wf(Ir-3) & +A(6,4)*Wf(Ir-2)+A(6,5)*Wf(Ir-1)+A(6,6)*Wf(Ir) WF2(IR)=B(6,1)*Wf(Ir-5)+B(6,2)*Wf(Ir-4)+B(6,3)*Wf(Ir-3) & +B(6,4)*Wf(Ir-2)+B(6,5)*Wf(Ir-1)+B(6,6)*Wf(Ir) C return end C C===================================================================== C complex*16 function Ylm(l,mm,theta,phi) C C This routine computes spherical harmonics (theta, phi in radians) C implicit real*8 (a-h,o-z) complex*16 cephi real*8 fak(0:130) common /gf2 / fak m=abs(mm) pi=3.141597265 x=cos(theta) sinp = sin(phi) cosp = cos(phi) ci=cmplx(0.,1.) cephi = cosp + ci*sinp if(m.gt.l.or.abs(x).gt.1.) then write(12,*)'bad arguments in Ylm' stop endif pmm=1. if(m.gt.0) then somx2=sqrt((1.-x)*(1.+x)) fact=1. do 11 i=1,m pmm=-pmm*fact*somx2 fact=fact+2. 11 continue endif if(l.eq.m) then plgndr=pmm else pmmp1=x*(2*m+1)*pmm if(l.eq.m+1) then plgndr=pmmp1 else do 12 ll=m+2,l pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m) pmm=pmmp1 pmmp1=pll 12 continue plgndr=pll endif endif Ylm=sqrt((2.*l+1.)/4./pi*fak(l-m)/fak(l+m))*plgndr*cephi**m if(mm.lt.0)Ylm=(-1.)**m*conjg(Ylm) return end C C===================================================================== C subroutine gfv(ie) C C Calculates sign and factorials of integers and half int. C C iv(n) = (-1)**n C fak(n) = n! C fad(n) = n!! C implicit real*8 (a-h,o-z) real*8 iv(0:130),fak(0:130),fad(0:130) common /gf1 / iv common /gf2 / fak common /gf3 / fad C if (ie.gt.130) then write(12,*) 'STOP IN GFV: IE LARGER 130 ' stop endif iv(0) = +1 fak(0) = 1.d0 fad(0)=1. fad(1)=1. do 10 i=1,ie iv(i) = -iv(i-1) 10 fak(i) = i*fak(i-1) do 11 i=3,ie,2 fad(i) = i*fad(i-2) 11 continue do 12 i=2,ie,2 fad(i) = i*fad(i-2) 12 continue return end C C The following routines are taken from the Numerical Recipes(c) C (C) Copr. 1986-92 Numerical Recipes Software. C****************** C Routines for Runge-Kutta integration begin here C C======================================================================= C SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs, *rkqs) C C Complex Runge-Kutta integration routine with adaptive stepsize control C implicit real*8(a-h,o-z) COMPLEX*16 ystart(nvar) PARAMETER (MAXSTP=10000,KMAXX=200,TINY=1.e-30) include 'relex.dim' COMPLEX*16 dydx(NSTMAX),y(NSTMAX),yp(NSTMAX,KMAXX), & yscal(NSTMAX) EXTERNAL derivs,rkqs REAL*8 xp(KMAXX) COMMON /path/ kmax,kount,dxsav,xp,yp x=x1 h=sign(h1,x2-x1) nok=0 nbad=0 kount=0 do 11 i=1,nvar y(i)=ystart(i) 11 continue if (kmax.gt.0) xsav=x-2.*dxsav do 16 nstp=1,MAXSTP call derivs(x,y,dydx) do 12 i=1,nvar yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY 12 continue if(kmax.gt.0)then if(abs(x-xsav).gt.abs(dxsav)) then if(kount.lt.kmax-1)then kount=kount+1 xp(kount)=x do 13 i=1,nvar yp(i,kount)=y(i) 13 continue xsav=x endif endif endif if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs) if(hdid.eq.h)then nok=nok+1 else nbad=nbad+1 endif if((x-x2)*(x2-x1).ge.0.)then do 14 i=1,nvar ystart(i)=y(i) 14 continue if(kmax.ne.0)then kount=kount+1 xp(kount)=x do 15 i=1,nvar yp(i,kount)=y(i) 15 continue endif return endif if(abs(hnext).lt.hmin) then write(12,*)'stepsize smaller than minimum in odeint' stop endif h=hnext 16 continue write(12,*)'too many steps in odeint' stop return END C C===================================================================== C SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs) C implicit real*8(a-h,o-z) COMPLEX*16 dydx(n),y(n),yscal(n) EXTERNAL derivs include 'relex.dim' CU USES derivs,rkck COMPLEX*16 yerr(NSTMAX),ytemp(NSTMAX) PARAMETER (SAFETY=0.9,PGROW=-.2,PSHRNK=-.25,ERRCON=1.89e-4) h=htry 1 call rkck(y,dydx,n,x,h,ytemp,yerr,derivs) errmax=0. do 11 i=1,n errmax=max(errmax,abs(yerr(i)/yscal(i))) 11 continue errmax=errmax/eps if(errmax.gt.1.)then h=SAFETY*h*(errmax**PSHRNK) if(h.lt.0.1*h)then h=.1*h endif xnew=x+h if(xnew.eq.x) then write(12,*)'stepsize underflow in rkqs' stop endif goto 1 else if(errmax.gt.ERRCON)then hnext=SAFETY*h*(errmax**PGROW) else hnext=5.*h endif hdid=h x=x+h do 12 i=1,n y(i)=ytemp(i) 12 continue return endif END C C===================================================================== C SUBROUTINE rkck(y,dydr,n,x,h,yout,yerr,derivs) C implicit real*8(a-h,o-z) COMPLEX*16 dydr(n),y(n),yerr(n),yout(n) EXTERNAL derivs include 'relex.dim' CU USES derivs COMPLEX*16 ak2(NSTMAX),ak3(NSTMAX),ak4(NSTMAX),ak5(NSTMAX), & ak6(NSTMAX),ytemp(NSTMAX) PARAMETER (A2=.2,A3=.3,A4=.6,A5=1.,A6=.875,B21=.2,B31=3./40., *B32=9./40.,B41=.3,B42=-.9,B43=1.2,B51=-11./54.,B52=2.5, *B53=-70./27.,B54=35./27.,B61=1631./55296.,B62=175./512., *B63=575./13824.,B64=44275./110592.,B65=253./4096.,C1=37./378., *C3=250./621.,C4=125./594.,C6=512./1771.,DC1=C1-2825./27648., *DC3=C3-18575./48384.,DC4=C4-13525./55296.,DC5=-277./14336., *DC6=C6-.25) do 11 i=1,n ytemp(i)=y(i)+B21*h*dydr(i) 11 continue call derivs(x+A2*h,ytemp,ak2) do 12 i=1,n ytemp(i)=y(i)+h*(B31*dydr(i)+B32*ak2(i)) 12 continue call derivs(x+A3*h,ytemp,ak3) do 13 i=1,n ytemp(i)=y(i)+h*(B41*dydr(i)+B42*ak2(i)+B43*ak3(i)) 13 continue call derivs(x+A4*h,ytemp,ak4) do 14 i=1,n ytemp(i)=y(i)+h*(B51*dydr(i)+B52*ak2(i)+B53*ak3(i)+B54*ak4(i)) 14 continue call derivs(x+A5*h,ytemp,ak5) do 15 i=1,n ytemp(i)=y(i)+h*(B61*dydr(i)+B62*ak2(i)+B63*ak3(i)+B64*ak4(i)+ *B65*ak5(i)) 15 continue call derivs(x+A6*h,ytemp,ak6) do 16 i=1,n yout(i)=y(i)+h*(C1*dydr(i)+C3*ak3(i)+C4*ak4(i)+C6*ak6(i)) 16 continue do 17 i=1,n yerr(i)=h*(DC1*dydr(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i)+DC6* *ak6(i)) 17 continue return END C C===================================================================== C SUBROUTINE rk4(y,dydx,n,x,h,yout,derivs) C implicit real*8(a-h,o-z) COMPLEX*16 dydx(n),y(n),yout(n) EXTERNAL derivs include 'relex.dim' COMPLEX*16 dym(NSTMAX),dyt(NSTMAX),yt(NSTMAX) hh=h*0.5 h6=h/6. xh=x+hh 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 yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i)) 14 continue return END C C End routines for Runge-Kutta integration C****************** C Routines for spline interpolation begin here C C===================================================================== C SUBROUTINE spline(x,y,n,yp1,ypn,y2) C implicit real*8(a-h,o-z) REAL*8 x(n),y(n),y2(n) PARAMETER (NMAX=500) REAL*8 u(NMAX) if (yp1.gt..99e30) then y2(1)=0. u(1)=0. else y2(1)=-0.5 u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) endif do 11 i=2,n-1 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) p=sig*y2(i-1)+2. y2(i)=(sig-1.)/p u(i)=(6.*((y(i+1)-y(i))/(x(i+ *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig* *u(i-1))/p 11 continue if (ypn.gt..99e30) then qn=0. un=0. else qn=0.5 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) do 12 k=n-1,1,-1 y2(k)=y2(k)*y2(k+1)+u(k) 12 continue return END C===================================================================== C SUBROUTINE splint(xa,ya,y2a,n,x,y) C implicit real*8(a-h,o-z) REAL*8 xa(n),y2a(n),ya(n) klo=1 khi=n 1 if (khi-klo.gt.1) then k=(khi+klo)/2 if(xa(k).gt.x)then khi=k else klo=k endif goto 1 endif h=xa(khi)-xa(klo) if (h.eq.0.) then write(12,*)'bad xa input in splint' stop endif a=(xa(khi)-x)/h b=(x-xa(klo))/h y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** *2)/6. return END C C End routines for spline interpolation C****************** C Routines for Simpson integration begin here C======================================================================= c SUBROUTINE trapzd(X,Y,Y2,Nu,a,b,s,n) c c Integration of a real function by trapezoidal rule. C N = integration points, A (B) is lower (upper) limit, Y (Y2) is c a function (its second derivative) calculated at Nu array points X. C Integral result = S c implicit real*8(a-h,o-z) REAL*8 x(nu),y(nu),Y2(nu) IF (n.EQ.1) THEN CALL splint(X,Y,Y2,nu,a,funca) CALL splint(X,Y,Y2,nu,b,funcb) s = .5*(b-a) * (funca + funcb) ELSE it = 2**(n-2) tnm = it del = (b-a)/tnm xx = a+.5*del sum=.0 DO j= 1,it CALL splint(x,y,y2,nu,xx,funcx) sum = sum + funcx xx = xx + del ENDDO s = .5 * (s + (b-a)*sum/tnm) ENDIF RETURN END c C======================================================================= c SUBROUTINE qsimp(x,y,y2,n,a,b,s) c c Integration of a real function by Simpson's rule. Use routine TRAPZD. C N = integration points, A (B) is lower (upper) limit, Y (Y2) is c a function (its second derivative) calculated at N array points X. C Integral result = S c implicit real*8(a-h,o-z) REAL*8 x(n),y(n),y2(n) PARAMETER( eps = 1.e-6,jmax = 20) ost = -1.e30 os = -1.e30 DO j = 1,jmax CALL trapzd(x,y,y2,n,a,b,st,j) s = (4.*st-ost)/3. IF (abs(s-os) .lt. eps*abs(os)) RETURN os = s ost = st ENDDO WRITE(12,*)' too many steps in qsimp' STOP END C End Routines for Simpson integration C-------------------------------------------------------------------- C END PROGRAM C--------------------------------------------------------------------