c=====================================================================
c
c	  PROGRAM DGDR
c                            Version Sept/1997 
c  based on the article by Bertulani, Canto, Hussein and Piza,
c  Phys. Rev. C53 (1996) 334
c                                   
c=====================================================================
c
c calculates 1 and 2 phonon excitation. Included are:
c     (a) Giant Dipole Resonance, IVGDR
c     (b) Isoscalar Giant Quadrupole Resonance, ISGQR
c     (c) Isovector Giant Quadrupole Resonance, IVGDR
c     (d) Double-GDR resonance with:
c                (d.1) L=2
c                (d.2) L=0
c
c In refs. [1,2] calculations are perturbative. In ref. [3] a
c coupled-chanels calculation was firstly done. Here the ideas
c of ref. [3] are extended to include:
c     (a) a correct description of the time dependence of the
c         electric dipole (E1) and quadrupole (E2) retarded
c         potentials, respectively
c     (b) all azimuthal substates properly
c     (c) the contributions of ISGQR and IVGQR states
c
c The calculation is non-perturbative, i.e., a coupled-channels
c calculation, where the strongest coupling, namely, the 
c                G.S.  <------>  GDR
c one, is calculated exactly. The time dependent occupation amplitude
c of the G.S. is used further to determine the transition amplitude 
c to the other states.
c
c Refs.
c         [1] PRL 71 (1993) 967, NPA 568 (1994) 931
c         [2] PRL 72 (1994) 1168
c         [3] PRL 72 (1994) 2147
c
c
c -------------------------------------------------------------------
c Here, we consider states labeled as:
c
c   i=0 (G.S.), 
c   i=1,2,3 (GDR states with l=1 and m=-1,0,1), 
c   i=4,5,6,7,8 (double-GDR states with l=2 and m=-2,-1,0,1,2), i=9 (l=0)
c   i=10-14 (ISGQR states with l=2 and m=-2,...,2)
c   i=15-19 (IVGQR states with l=2 and m=-2,...,2)
c          
c--------------------------------------------------------------------
c
      include 'dgdr_3d.dim'
c
c--------------------------------------------------------------------
c
      implicit real*8(a-h,o-z)
      complex*16 x(nmx),dlvc(3),y,a0(n1,nn),vc(0:n,0:n),za(nn),cint,
     &           zb(n,nn),z1(3,nn),zz,c2,c3
      dimension prob(n1,n1),xa(n1),xb(n1),p1e(n1),p2e(n1),xs(n1),
     &          xv(n1),pse(n1),pve(n1),probs(n1,n1),probv(n1,n1),
     &          prob2(n1,n1),x2r(n1) 
      common / par  / b,gamma,beta,s01,s12,s02,eps1,qsi1,pi
      common / sys  / zp,ap,zt,at,epn,homeg1,homeg2,gm1,e0
      common / gqr  / eisgqr,eivgqr,sqs,sqv
      common / w1 /  w011(3,3),w110(3,3),w022(4:8,4:8)
      common / w2 /  w112m1(3,4:8),w112ze(3,4:8),w112p1(3,4:8)
      common / spin / ll(0:n),mm(0:n)
      common / dlvc1 / y,c2,c3
      common / vcoup1 / c0,bb,fac02,facqs,facqv
      common / flag / iflag
      external dxdt
c
      open(11,file='dgdr_3d.in',status='old')
      open(12,file='dgdr_3d.out',status='unknown')
      open(13,file='dgdr_3d.dat',status='unknown')
      open(15,file='a0t.dat',status='unknown')
      open(17,file='lixo.dat',status='unknown')
c
      read(11,*) zp,ap,zt,at,epn
      read(11,*) homeg1,gm1,homeg2,emin,emax,de,bmin,bmax,db
      read(11,*) tcut,dt,frac01,frac12,frac02
      read(11,*) eisgqr,gms,eivgqr,gmv,frac2s,frac2v
c
c--------------------------------------------------------------------
c zp,ap,zt,at = atomic and mass numbers of projectile and target. 
c               (We assume that the excitation is in the target)
c         epn = lab energy per nucleon
c      homeg1 = GDR excitation energy (hbar omega)
c         gm1 = width of the resonance
c      homeg2 = 2 phonon GDR energy (hbar omega)
c         gm2 = 2 phonon GDR width
c      eisgqr = energy of the isoscalar giant quadrupole resonance
c         gms = width of the ISGQR  
c      eivgqr = energy of the isovector giant quadrupole resonance
c         gmv = width of the IVGQR  
c      erange = energy range around the GDR energy taken into account
c          de = energy mesh spacing
c        bmax = maximal impact parameter in the cross sections 
c          db = impact parameter mesh spacing
c        tcut = dimensionless time cut-off in for the time integration
c          dt = dimensionless mesh spacing for the time integration
c      frac01 = fraction of the dipole sum rule exhausted 0-->1
c      frac12 = fraction of the sum rule exhausted for 1-->2
c      frac02 = fraction of the quadrupole sum rule exhausted 
c               by the double-GDR 0-->2xGDR (L=2)
c      frac2s = fraction of the quadrupole sum rule exhausted by
c               the isoscalar usual quadrupole resonance
c      frac2v = fraction of the quadrupole sum rule exhausted by
c               the isovector usual quadrupole resonance
c --------------------------------------------------------------------
c
c L and M values for the states:
c
      ll(0)=0
      mm(0)=0
      do i=1,3
	ll(i)=1
	mm(i)=i-2
      end do
      do i=4,8
	ll(i)=2
	mm(i)=i-6
      end do
      ll(9)=0
      mm(9)=0
      do i=10,14
	ll(i)=2
	mm(i)=i-12
      enddo
      do i=15,19
	ll(i)=2 
	mm(i)=i-17
      enddo
c
c Wigner 3j symbols
c
      ze=0.d0
      on=1.d0
      tw=2.d0
      do i=1,3
	 fi=dfloat(mm(i)) 
	 do j=1,3
	    fj=dfloat(mm(j))
	    w011(i,j)=wigner(ze,on,on,ze,fi,fj)
	    w110(i,j)=wigner(on,on,ze,fi,fj,ze)
	 enddo
      enddo
      do i=1,3
	 fi=dfloat(mm(i))
	 do j=4,8
	    fj=dfloat(mm(j))
	    w112m1(i,j)=wigner(on,on,tw,-on,fi,fj)
	    w112ze(i,j)=wigner(on,on,tw,ze,fi,fj)
	    w112p1(i,j)=wigner(on,on,tw,on,fi,fj)
	 enddo
      enddo
      do i=4,8
	 fi=dfloat(mm(i))
	 do j=4,8
	    fj=dfloat(mm(j))
	    w022(i,j)=wigner(ze,tw,tw,ze,fi,fj)
	 enddo
      enddo
c
c constants      
c
      pi=4.*atan(1.d0)
      y=dcmplx(0.d0,1.d0)
      hc=197.33
c
c width of the double GDR 
c
      gm2=sqrt(2.)*gm1
c      
c coupling strengths as given by sum rules (in units of e^2):      
c sum rules are described in ref. [2]
c
      sd=9./4./pi*20.7*zt*(at-zt)/at
      s01=sqrt(frac01*sd/homeg1)
      s12=sqrt(frac12*2.*sd/(homeg2-homeg1))
      rt=1.2*at**(1./3.)
      sq=2.*15./4./pi*20.7*rt**2*zt**2/at
      sqs=sqrt(frac2s*sq/eisgqr)
      sqv=sqrt(frac2v*sq*(at-zt)/zt/eivgqr)
      s02=sqrt(frac02*sq/homeg2)
      c0=zp/137.
c      
      gamma=1.+epn/931.5
      beta=sqrt(gamma*gamma-1.)/gamma
c
c details of the time, energy and impact-parameter meshes
c
      ncut=tcut/dt
      erange=emax-emin
      ne=erange/de
      if (ne.lt.1)ne=1
      nb=(bmax-bmin)/db
      if(nb.le.0)nb=1
      nt=2*ncut+1
c
c conversion to odd number of integration points (required by Simpson 
c routines)
c
      if((-1)**ne.eq.1)ne=ne+1
      if((-1)**nb.eq.1)nb=nb+1
c
c********        OUTPUT  # 1      
c      
      write(12,100)zp,ap,zt,at,epn,frac01,frac12,frac02,
     &             homeg1,gm1,homeg2,gm2,frac2s,frac2v,eisgqr,gms
     &             ,eivgqr,gmv,emax,ne,nb,nt
c
c
c _______________________________    start mesh in impact parameter 
c
      do jb=1,nb
	b=bmin+db*(jb-1)
	e0=gamma*beta*hc/b
	qsi1=gm1/e0
c       
c eps1 is the real part of the pole of the Lorenzian, which contributes in
c energy integrations
c
	eps1=sqrt(homeg1**2-gm1**2/4.)/e0
c
	c2= y*(eps1-y*qsi1/2.)
	c3=dcmplx(1.d0,-qsi1/eps1/2.)
	bb=beta*b
	fac=sqrt(pi/30.)*c0/bb/b
	fac02=fac*s02
	facqs=fac*sqs
	facqv=fac*sqv
c
c
	ndif=6
	tau=-tcut
c       
c initialization of the variables x(i)
c
	do ijk=1,6
	  x(ijk)=cmplx(0.,0.)
	end do
	a0(jb,1)=cmplx(1.,0.)
c
c
	iflag=1
	call vcoup(tau,vc,dlvc)
	do jmu=1,3
	   zb(jmu,1)=conjg(vc(0,jmu))
     &     *exp(y*dcmplx(eps1,-qsi1/2.)*tau)
	   z1(jmu,1)=dcmplx(0.d0,0.d0)
	end do
c

c ______________________________   start coupled-channels calculation
c
	do jt= 1,nt-1
	  call rk4(x,ndif,tau,dt,dxdt)
	  a0(jb,jt+1)=dcmplx(1.d0,0.d0)+x(1)+x(2)+x(3)
	  call vcoup(tau,vc,dlvc)
c
	  do jmu=1,3
	    zb(jmu,jt+1)=conjg(vc(0,jmu))
     &      *exp(y*dcmplx(eps1,-qsi1/2.)*tau)*a0(jb,jt+1)
c
c zb(jmu,jt) is the integrand corresponding to the second integral
c to calculate the DGDR population, appearing in eq. (64) of Phys.
c Rev. C57 (1996) 334. Notice that there is a misprint in this eq.
c There should be a complex conjugate in the potential V_mu^{(01)}.
c
	  end do
	end do
c
c__________________________________ end coupled-channels calculation
c
c
	p0b=abs(a0(jb,nt))**2
c
c
c ****************** Carrying out the first integration in eq. (64): 
c
c z1(\mu,t') = \int_{-infty}^t dt' V_\mu^{(01)}(t''} .........
c
c ******************
	do jmu=1,3
c
c the integration is carried out for odd number of points 
c (requirement of the Simpson routine). For even numbers 
c it is taken as the average of neighbouring values
c
	   do jt= 3,nt,2
	      tau=-tcut+dt*dble(jt-1)
	      call vcoup(tau,vc,dlvc)
	      do jj=1,jt
		 za(jj)=zb(jmu,jj)
	      end do
	      call csimp(jt,za,dt,cint)
	      z1(jmu,jt)=cint
	      z1(jmu,jt-1)=(z1(jmu,jt-2)+z1(jmu,jt))/2.
	   end do
	end do
c
c----------------------------------------------------------------
c
c
c**********       OUTPUT  # 2
c
c        write(*,*)'p0b=',p0b
c
c**************************  energy spectrum of one-phonon states
c
	do je=1,ne
	  prob(jb,je)=0.
	  probs(jb,je)=0.
	  probv(jb,je)=0.
	  e=emin+de*(je-1)
	  ea2=e*e-homeg1**2
	  eas=e*e-eisgqr**2
	  eav=e*e-eivgqr**2
	  efac = 2.*gm1/pi*e*e/(ea2**2+(gm1*e)**2)
	  efacs = 2.*gms/pi*e*e/(eas**2+(gms*e)**2)
	  efacv = 2.*gmv/pi*e*e/(eav**2+(gmv*e)**2)
	  ee=e/e0
c          
c giant dipole          
c
	  do i=1,3                           
	    do jt=1,nt
	      tau=-tcut+dt*(jt-1)
	      iflag=1
	      call vcoup(tau,vc,dlvc)
	      za(jt)=conjg(vc(0,i))*exp(y*ee*tau)*a0(jb,jt)
	    end do
	    call csimp(nt,za,dt,cint)
	    prob(jb,je)=prob(jb,je)+efac*abs(cint)**2
	  end do
c
c giant isoscalar quadrupole          
c          
	  do i=10,14                           
	    do jt=1,nt
	      tau=-tcut+dt*(jt-1)
	      iflag=4
	      call vcoup(tau,vc,dlvc)
	      za(jt)=conjg(vc(0,i))*exp(y*ee*tau)*a0(jb,jt)
	    end do
	    call csimp(nt,za,dt,cint)
	    probs(jb,je)=probs(jb,je)+efacs*abs(cint)**2
	  end do
c
c giant isovector quadrupole          
c          
	  do i=15,19                           
	    do jt=1,nt
	      tau=-tcut+dt*(jt-1)
	      iflag=5
	      call vcoup(tau,vc,dlvc)
	      za(jt)=conjg(vc(0,i))*exp(y*ee*tau)*a0(jb,jt)
	    end do
	    call csimp(nt,za,dt,cint)
	    probv(jb,je)=probv(jb,je)+efacv*abs(cint)**2
	  end do
	end do
c
c*************      end spectrum of one-phonon states        
c
c***************************  energy spectrum of two-phonon states
c
	do je=1,ne
	  prob2(jb,je)=0.
	  auxx=0.
	  e=emin+de*(je-1)
	  ea2=e*e-homeg2**2
	  efac = 2.*gm2/pi*e*e/(ea2**2+(gm2*e)**2)
	  ee=e/e0
	  do i=4,9                          
	    do jt=1,nt
	      tau=-tcut+dt*(jt-1)
	      iflag=2
	      call vcoup(tau,vc,dlvc)
	      zz=dcmplx(0.,0.)
	      do j=1,3
		zz=zz+conjg(vc(j,i))*dcmplx(1.d0,-qsi1/eps1/2.)*
     &                 exp(-y*tau*dcmplx(eps1,-qsi1/2.))*z1(j,jt)
	      end do
	      iflag=3
	      call vcoup(tau,vc,dlvc)
c
c WARNING. Direct excitation is set to ZERO ! If  two-step excitation
c          is wanted, comment the next 2 lines and remove the comment
c          in the line below the next warning.
c              
	      vc(0,i)=0. 
	      za(jt)=zz*exp(y*ee*tau)
     &              +conjg(vc(0,i))*a0(jb,jt)*exp(y*ee*tau)
c
c WARNING. Two-step excitation is set to zero !
c             za(jt)=(conjg(vc(0,i))*a0(jb,jt))*exp(y*ee*tau)
c            
	    end do
	    call csimp(nt,za,dt,cint)
	    auxx=auxx+abs(cint)**2
	  end do
	  prob2(jb,je)=efac*auxx
	end do      
c          
c************     end spectrum of double-GDR states
c
c*************     Checking unitarity
c
	if(ne.eq.1)then
	    p1b=prob(jb,ne)
	    p2b=prob2(jb,ne)
	else
	  do je=1,ne
	    xa(je)=prob(jb,je)
	    x2r(je)=prob2(jb,je)
	  end do
	  p1b=arsimp(ne,xa,de)
	  p2b=arsimp(ne,x2r,de)
	end if
c
c************      OUTPUT  # 3
c
	write(12,155)b,p0b,p1b,p0b+p1b,p2b
c
c*********       end check of unitarity
c
      end do
c
c___________________________________    end mesh in impact parameter
c
c*********   integration over impact parameter
c
      do je=1,ne
	e=emin+de*(je-1)
	if(nb.eq.1)then
	    p1e(je)=20.*pi*b*prob(1,je)
	    p2e(je)=20.*pi*b*prob2(1,je)
	    pse(je)=20.*pi*b*probs(1,je)
	    pve(je)=20.*pi*b*probv(1,je)
	else
	  do jb=1,nb
	    b=bmin+db*(jb-1)
	    xa(jb)=20.*pi*b*prob(jb,je)
	    xb(jb)=20.*pi*b*prob2(jb,je)
	    xs(jb)=20.*pi*b*probs(jb,je)
	    xv(jb)=20.*pi*b*probv(jb,je)
	  end do
	  p1e(je)=arsimp(nb,xa,db)
	  p2e(je)=arsimp(nb,xb,db)
	  pse(je)=arsimp(nb,xs,db)
	  pve(je)=arsimp(nb,xv,db)
	end if
	ptot=p1e(je)+p2e(je)+pse(je)+pve(je)
c
c**********   OUTPUT # 4 : final spectrum (d sigma / dE)
c
	write(13,160)e,p1e(je),p2e(je),pse(je),pve(je),ptot
      end do
c
c************    end integration over impact parameter
c
c************     Total cross sections
c
      sig1=arsimp(ne,p1e,de)        
      sig2=arsimp(ne,p2e,de)        
      sigs=arsimp(ne,pse,de)        
      sigv=arsimp(ne,pve,de)        
c      
c*********     OUTPUT # 5 : cross sections
c
      write(12,170)sig1,sig2,sigs,sigv
c
100   format(' *******************************************************'
     &,//, ' checking input data:',//,' zp,ap,zt,at,epn =',5f6.0,/
     @     ' frac01,frac12,frac02 = ',3(g10.2,1x),/,
     & ' hw1,gm1,hw2,gm2 = ',4(f8.1),/,
     & ' frac2s,frac2v = ',2(g10.2,2x),/,
     & ' eisgqr,gms,eivgqr,gmv = ',4(f8.1),/,' emax,ne,nb,nt = ',
     & f6.2,3(2x,i4),/,
     &' *******************************************************',//,
     &'   b(fm) ',9x,'P_0 ',10x,'P_1 ',10x,'norm',10x,'P_2'/)
155   format(g12.4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4)
160   format(6(1x,e12.4))     
170   format(' Total GDR   cross section (mb) = ',g12.4,/,
     &       ' Total 2xGDR cross section (mb) = ',g12.4,/,
     &       ' Total ISGQR cross section (mb) = ',g12.4,/,
     &       ' Total IVGQR cross section (mb) = ',g12.4)
      stop
      end         
c         
c         
c======================================================================
c
      subroutine dxdt(tau,x,dx)
c
c=====================================================================
c
      include 'dgdr_3d.dim'
c
c---------------------------------------------------------------------      
      implicit real*8(a-h,o-z)
      complex*16 y,x(nmx),dx(nmx),vc(0:n,0:n),dlvc(3),c1,c2,c3
      common / par  / b,gamma,beta,s01,s12,s02,eps1,qsi1,pi
      common / dlvc1 / y,c2,c3
c
      call vcoup(tau,vc,dlvc)
      c1=(cmplx(1.,0.)+x(1)+x(2)+x(3))*c3
      do i=1,3
	dx(i)=x(i+3)
	dx(i+3)=(dlvc(i)-c2)*x(i+3)-abs(vc(0,i))**2*c1
      end do
      return
      end
c
c======================================================================
c
	subroutine rk4(y,m,x,h,derivs)
c
c======================================================================
c * complex runge-kutta integration routine (from numerical recipes) *
c----------------------------------------------------------------------
c
	include 'dgdr_3d.dim'
c
c----------------------------------------------------------------------
	implicit real*8(a-h,o-z)
	complex*16 y,dydx,yt,dyt,dym
	dimension y(nmx),dydx(nmx),yt(nmx),dyt(nmx),dym(nmx)
	external derivs
c
	hh=h*.5
	h6=h/6.
	xh=x+hh
c
	call derivs(x,y,dydx)
	do 11 i=1,m
	yt(i)=y(i)+hh*dydx(i)
 11     continue
	call derivs(xh,yt,dyt)
	do 12 i=1,m
	yt(i)=y(i)+hh*dyt(i)
 12     continue
	call derivs(xh,yt,dym)
	do 13 i=1,m
	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,m
	y(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
 14     continue
	x=x+h
	return
	end
c
c
c----------------------------------------------------------------------
c
      subroutine vcoup(t,vc,dlvc)
c
c  for IFLAG =
c               1     returns V(0 --> GDR)
c               2     returns V(GDR --> 2xGDR, L=0, and L=2)
c               3     returns V(0 --> 2XGDR, L=0, and L=2)
c               4     returns V(0 --> ISGQR)
c               5     returns V(0 --> IVGDR)
c---------------------------------------------------------------------
c
	include 'dgdr_3d.dim'
c
c--------------------------------------------------------------------      
      implicit real*8(a-h,o-z)
      complex*16 vc(0:n,0:n),dlvc(3),ab,ac,xn,xd
c
c
      common / par  / b,gamma,beta,s01,s12,s02,eps1,qsi1,pi
      common / sys  / zp,ap,zt,at,epn,homeg1,homeg2,gm1,e0
      common / spin / ll(0:n),mm(0:n)
      common / w1 /  w011(3,3),w110(3,3),w022(4:8,4:8)
      common / w2 /  w112m1(3,4:8),w112ze(3,4:8),w112p1(3,4:8)
      common / gqr  / eisgqr,eivgqr,sqs,sqv
      common / vcoup1 / c0,bb,fac02,facqs,facqv
      common / flag / iflag
c
c******  FLAG
c
      sqt=sqrt(1.d0+t*t)
      sqt5=sqt**5
      go to (11,12,13,14,15)iflag    
c
c*** dipole interaction: 0-->1 (GDR)
c
11    do i=1,3
	m=4-i
	dx=s01*sqrt(2.*pi/3.)*(w011(1,m)-
     &         w011(3,m))
	dz=s01*sqrt(4.*pi/3.)*w011(2,m)
	ab=c0*dz/bb/gamma
	ac=c0*dx/bb
	xn=ab*t+ac
	xd=sqt**3
	vc(0,i)=(-1.)**(1-mm(i))*xn/xd
	vc(i,0)=conjg(vc(0,i))
c        
c*** logarithimic derivative of the interaction
c
	dlvc(i)=dcmplx(0.d0,0.d0)
	if (abs(xn).gt.1.d-8) then
	  dlvc(i)=ab/xn - 3.*t/sqt**2
	end if
      enddo
      goto 16
c        
c*** dipole interaction: 1-->2xGDR (L=0)
c
12      continue
	do i=1,3  
	  dx=s12*sqrt(2.*pi/3.)*(w110(i,1)-
     &         w110(i,3))
	  dz=s12*sqrt(4.*pi/3.)*w110(i,2)
	  ab=c0*dz/bb/gamma
	  ac=c0*dx/bb
	  xn=ab*t+ac
	  xd=sqt**3
	  vc(i,9)=xn/xd
	  vc(9,i)=conjg(vc(i,9))
      enddo
c        
c*** dipole interaction: 1-->2xGDR (L=2)
c
      do i=1,3  
	do k=4,8
	  m=12-k 
	  dx=s12*sqrt(2.*pi/3.)*(w112m1(i,m)-
     &         w112p1(i,m))
	  dz=s12*sqrt(4.*pi/3.)*w112ze(i,m)
	  ab=c0*dz/bb/gamma
	  ac=c0*dx/bb
	  xn=ab*t+ac
	  xd=sqt**3
	  vc(i,k)=(-1.)**(2-mm(k))*xn/xd
	  vc(k,i)=conjg(vc(i,k))
	enddo
      enddo
      goto 16
c        
c*** quadrupole interaction for the transition 0-->2xGDR (L=2)
c
13    continue
      do i=4,8
	m=12-i
	q1=w022(8,m)+w022(4,m)
	q2=w022(7,m)-w022(5,m)
	q3=sqrt(6.)*w022(6,m)
	vc(0,i)=-(-1.)**(2-mm(i))*fac02/sqt5*(q1*3.
     &     +3.*(2.-beta**2)*q2*gamma*t
     &     +q3*(2.*t*t-1.))
	vc(i,0)=conjg(vc(0,i))
      end do
      goto 16
c
c*** quadrupole interaction for the transition 0-->ISGQR
c
14    continue
      do i=10,14
	m=18-i 
	q1=w022(8,m)+w022(4,m)
	q2=w022(7,m)-w022(5,m)
	q3=sqrt(6.)*w022(6,m)
	vc(0,i)=-(-1.)**(2-mm(i))*facqs/sqt5*(q1*3.
     &     +3.*(2.-beta**2)*q2*gamma*t
     &     +q3*(2.*t*t-1.))
	vc(i,0)=conjg(vc(0,i))
      end do
      goto 16
c
c*** quadrupole interaction for the transition 0-->IVGQR
c
15    continue
      do i=15,19
	m=23-i
	q1=w022(8,m)+w022(4,m)
	q2=w022(7,m)-w022(5,m)
	q3=sqrt(6.)*w022(6,m)
	vc(0,i)=-(-1.)**(2-mm(i))*facqv/sqt5*(q1*3.
     &     +3.*(2.-beta**2)*q2*gamma*t
     &     +q3*(2.*t*t-1.))
	vc(i,0)=conjg(vc(0,i))
      end do
c
16    return
      end
c
c
C=======================================================================
c
       real*8 function arsimp(nm,a,h) 
c
c======================================================================
c
c Integration of a real function by Sympson's rule. N (must be odd)  is 
c the number mesh points, A is a vector containing the integrand values
c on the mesh and H is the step size.
c
c----------------------------------------------------------------------
c
       implicit real*8(a-h,o-z)                
       include 'dgdr_3d.dim'
       dimension a(n1)
       l=nm                                     
       s=a(1)-a(l)                           
       do 1 i=    2,l,2                            
 1     s=s+4.d0*a(i)+2.d0*a(i+1)           
       sum=h*s/3.d0                       
       arsimp = sum
       return                                  
       end
c
c
c---------------------------------------------------------------------
       subroutine csimp(nm,a,h,sum) 
       implicit real*8(a-h,o-z)                
       include 'dgdr_3d.dim'
       complex*16 a(nn),s,sum                        
       l=nm                                     
       s=a(1)-a(l)                           
       do 1 i=2,l,2                            
 1     s=s+4.d0*a(i)+2.d0*a(i+1)           
       sum=h*s/3.d0                       
       return                                  
       end
c
C
C***********************************************************************
C
C     Returns Wigner 3j-coefficients (for Clebsh coeficients change
C     name to clebsh and delete line 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
C                               Written by C.A. Bertulani, Dec. 1994
C
C***********************************************************************
C      
      real*8 function wigner(a1,a2,a3,b1,b2,b3)
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. 
      wigner=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
      wigner=(-1)**nint(a1-a2+b3)/sqrt(2.*a3+1.)*clebsh
100   b3=-b3
      return
      end