cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This program generates the results discussed in: 
c  "Eigenvalue problem for a coupled channel Schrodinger equation
c  with application to the description of deformed nuclear systems"
c  V. Ledoux, M. Rizea, M. Van Daele, G. Vanden Berghe, I. Silisteanu
c
c  Program to calculate the eigenvalues and eigenfunctions of a
c  system of Schrodinger equations of the form:
c  Y'' = (U(r) + \lambda V(r)) Y
c  \lambda is the spectral parameter (depth of the nuclear potential)
c
c  A constant perturbation method is used for propagation 
c
c  The program uses some subroutines from the CPC Library 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c   declariations/initializations:
c
c     set dimensions 
      parameter(ndim=2,ndim2=ndim*ndim,ndd=2*ndim,ndd1=ndd-1)       
      parameter(kmax=20000)
      parameter(neq=2*ndim2,neq2=2*neq)
      parameter(maxlog=14,maxit=15,mxd=3*maxit+1,lwrk=mxd*maxlog*ndim)
      parameter(ndeg=3,npoints=kmax)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
      real*8 hbar
      common /valene/ ce(ndim)
	common /valp/ ck(ndim)
	common /eigend/ depth
	common /indexes/ i,j
      common /valq/ cq
      common /slopes/ qq(ndim,ndim)
      common /radii/ rinit,rmatch,rfin,rdx
      common /radial/ raz(0:kmax),kfin
      common /matchn/ kmatc
      common /tol/ rtol
      common /dprec/ rprec
      common /solsd/ sy(neq),syp(neq),dy(neq),dyp(neq),wght(ndd)
      common /solut/ sol(kmax,ndim2),dsol(kmax,ndim2)
	common /pnuclear/ pnuc(kmax,ndim2)
      common /totsol/ tsol(kmax,ndim),tdsol(kmax,ndim)
      common /solutp/ solp(0:kmax,ndim),dsolp(0:kmax,ndim)
	common /workn/ w(4*ndim),wm(4*ndim)
	common/poti1/v1c(ndim,ndim),v10(ndim,ndim),v11(ndim,ndim),
     &  v12(ndim,ndim),lo(ndim)   
	common/worki1/aw(lwrk),yw(lwrk),adw(lwrk),ydw(lwrk),cw(ndim),
     ' reps,iterat(ndim),ml
	common /fituv/ pu(ndim2,0:ndeg),pv(ndim2,0:ndeg)
	common /fituv0/ pu0(ndim2,kmax),pv0(ndim2,kmax)
	common /fituv1/ pu1(ndim2,kmax),pv1(ndim2,kmax)
	common /fituv2/ pu2(ndim2,kmax),pv2(ndim2,kmax)
	external fctu,fctv
	common /diagonal/ amatr(ndim,ndim),dval(ndim),zval(ndim,ndim),
     +  wval(ndim),tval(ndim,ndim),pval(ndim,ndim)
	common /fitmat/ fmatr(ndim2,0:ndeg),fuval(ndim,ndim),
     +  fvval(ndim,ndim),fval20c(ndim,ndim) 
	common /comat/ fval10(ndim,ndim),fval11(ndim,ndim),fval22(ndim,
     +  ndim),fval20(ndim,ndim),fval200(ndim,ndim),fval12(ndim,ndim)
	dimension rdelta(4)
	common /partit/ rpart(kmax),npart
      common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
	common /parcoul/ rz,rz2,r1,r2
	common /lorb/ lang(ndim)
      common /legend/ r,rh,k	
      common /jwrite/ ifin
	common /rowcol/ irow,jcol 
	common /match_sys/ amatch(neq,neq),bmatch(neq)
     +                    ,imatch(neq,neq)
	  
	dimension valq(20),valphi(20)
	dimension rp(neq),rxi(neq,neq)
      common/hion/rap,rzp,rrp,rat,rzt,rrt
	common /coulrad/ radc
      common/rott/rbeta2t,rbeta4t,ntarg
      common/dimen/nlevel
      common/valpi/rpi,hbar
      dimension renerg(ndim)
	character*4 name,namen
	character*2 numb(25)/
     +  '01','02','03','04','05','06','07','08','09','10',
     +  '11','12','13','14','15','16','17','18','19','20',
     +  '21','22','23','24','25'/
	data rhbarc/197.327053d0/,ruma/931.49432d0/,resq/1.4399748d0/ 
	data rxl/0.d0/,rxu/1.d0/
	rhbar=6.5821220d-22
	hbar=rhbarc
	reps=1.d-12
	rpi=4.d0*atan(1.d0)

c      input data file
       open(1,file='sys.inp')  

c      output files with potential values 
c       open(11,file='potu.11')
c       open(12,file='potu.12')
c       open(21,file='potu.21')
c       open(22,file='potu.22')
c       open(31,file='potv.11')
c       open(32,file='potv.12')
c       open(41,file='potv.21')
c       open(42,file='potv.22')

c     
c    read input data from sys.inp:
c         
c    'give masses target /projectile:'
      read(1,*) rat,rap 
      rmass=rat*rap/(rat+rap)*ruma  
      rm=rat*rap/(rat+rap)
      rhb22mu = (rhbarc * rhbarc)/ (2.d0 * rmass)
	r2muhb2=1.d0/rhb22mu
	rw0=0.d0
c    'give nuclear radius and diffuseness'
      read(1,*) rv0,rdif
      radp=rv0*rat**(1.d0/3.d0)
      r0=rv0*(rap**(1.d0/3.d0)+rat**(1.d0/3.d0))
      rrt=radp
c    'give deformations'
	read(1,*) rbeta2t,rbeta4t
c    'give l_i'
	do i=1,ndim
	read(1,*) lang(i)
	enddo
c    'give atomic masses: z=zt+zp, z2=zp' 
	read(1,*) rz,rz2
	rzt=rz-rz2
	rzp=rz2 
c    'give coulomb radii - r2 is used for the present application'
	read(1,*) r1,r2
      rrp=r2*rap**(1.d0/3.d0)
      radc=r2*(rap**(1.d0/3.d0)+rat**(1.d0/3.d0))
c    'give r_init and r_step'
	read(1,*) r_in,r_dx
      rdx=r_dx
c    'give (complex) energies:'
	do i=1,ndim
      read(1,*) rer,rei 
      renerg(i)=rer
      ce(i)=r2muhb2*dcmplx(rer,rei)
	ck(i)=cdsqrt(ce(i))
	enddo
	nlevel=ndim
	ntarg=ndim-1
c     input tolerance:              
      rtol=1.d-10 
	rprec=1.d-10
	itermax=50
c    'give matching point, endpoint'
      read(1,*) rmatc,rfin
	rfin=rfin+1.d-3*rfin
	kfin=rfin/rdx
	rfin=kfin*rdx
	
	np=kfin  
	
c    'give (complex) v0_init, initial value for depth:'

      read(1,*) rrv0,riv0  
	depth=dcmplx(rrv0,riv0)
       
c     output data file
      open(3,file='sys.rep')

c          'cmat0' generates the matrix $\hat O$ and obtain its eigenvalues 
c          and eigenfunctions (by 'mdiag'). To obtain the Wigner 3-j 
c          coefficients 3-j the subroutine 'cg' is used. cg calculates in fact 
c          other coefficients (Clebsch-Gordan) and then uses a relation 
c          between Clebsch-Gordan and Wigner coefficients
      call cmat0
	 
	rinit=r_in
      raz(0)=rinit
	rdx=r_dx 
	 
c     this sbr prepares the integration on I1: 
	call pi1(ndim,ndim,rinit,lo,v1c,v10,v11,v12,ml)

c
c     stepsize selection procedure (from rinit to radc):
c
      write(*,*)'constructing mesh...'
	r=rinit
	npn=1
	rpart(npn)=rinit
	radfin=radc
	rh=rdx
507	r=rpart(npn)
	rnext=r+rh
	if(rnext.ge.radfin) rh=radfin-r
	numiter=0
505	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	do k=0,ndeg
      call quadg16(rxl,rxu,fctu,cofk)
	if(k.eq.0) then
	pu(l+i,k)=(2*k+1)*cofk
	else
	pu(l+i,k)=(2*k+1)*rh**2*cofk
	endif
      call quadg16(rxl,rxu,fctv,cofk)
	if(k.eq.0) then
	pv(l+i,k)=(2*k+1)*cofk
	else
	pv(l+i,k)=(2*k+1)*rh**2*cofk
	endif
	enddo
	enddo
	enddo

c     compute {\bf F}_0
      l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
	amatr(i,j)=pu(l+i,0)+depth*pv(l+i,0)
	enddo
	enddo

c     dchouse computes the eigenvalues(dval) and 
c     eigenvectors of amatr
c     amatr= {\bf F}_0
	call dchouse(amatr,ndim,ndim,dval,zval,wval)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	if(i.eq.j) then
	fmatr(l+i,0)=dval(i)*rh**2
	else
	fmatr(l+i,0)=0.d0
	endif
	enddo
	enddo

c     compute fmatr = {\bf F}^D_i 
	call dcgmtr(zval,tval,ndim,ndim)
	do k=1,ndeg
	call dcgmpr(tval,pu(1,k),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fuval,ndim,ndim,ndim)
	call dcgmpr(tval,pv(1,k),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fvval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	fmatr(l+i,k)=fuval(i,j)+depth*fvval(i,j)
	enddo
	enddo
	enddo

c     compute estimate of local error:
c       f1*f2 in tval
	call dcgmpr(fmatr(1,1),fmatr(1,2),tval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=-0.5d0*fmatr(l+i,3)-(tval(i,j)-tval(j,i))/120.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax1,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,3),fmatr(1,0),fuval,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,1),fmatr(1,0),fvval,ndim,ndim,ndim)
	do j=1,ndim
	do i=1,ndim
	fval10(i,j)=fvval(i,j)-fvval(j,i)
	enddo
	enddo
	call dcgmpr(fval10,fmatr(1,0),fvval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=2.5d0*fmatr(l+i,3)+(fuval(i,j)-fuval(j,i))/8.d0-
     +  (tval(i,j)-tval(j,i))/120.d0-(fvval(i,j)-fvval(j,i))/96.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax2,ndim,ndim,ndim)
	
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=-0.75d0*(fuval(i,j)-fuval(j,i))+(tval(i,j)+
     +  tval(j,i))/4.d0-(tval(i,j)-tval(j,i))/40.d0-
     +  (fvval(i,j)-fvval(j,i))/32.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdelu=rpmax1/3.d0+rpmax2/15.d0+rpmax3/105.d0
	rdelta(1)=dabs(rdelu)
	call dcgmpr(fmatr(1,2),fmatr(1,2),fval22,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,2),fmatr(1,0),fval20,ndim,ndim,ndim)
	do j=1,ndim
	do i=1,ndim
	fval20c(i,j)=fval20(i,j)-fval20(j,i)
	enddo
	enddo
	call dcgmpr(fval20c,fmatr(1,0),fval200,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,1),fval10,fval11,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=(fuval(i,j)-fuval(j,i))/8.d0-fval22(i,j)/40.d0+
     +  (fval200(i,j)-fval200(j,i))/480.d0+(fval11(i,j)-fval11(j,i))/
     +  480.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax1,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=-5.d0*(fuval(i,j)-fuval(j,i))/8.d0-3.d0*fval22(i,j)/
     +  20.d0+(fval200(i,j)-fval200(j,i))/30.d0+
     +  (11.d0*fval11(i,j)-fval11(j,i))/480.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax2,ndim,ndim,ndim)
	l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
      pval(i,j)=57.d0*fval22(i,j)/40.d0+11.d0*fval11(j,i)/240.d0+
     +  fval11(i,j)/10.d0-9.d0*(fval200(i,j)-fval200(j,i))/160.d0
      enddo
      enddo
      call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdeldu=rpmax1/3.d0+rpmax2/15.d0+rpmax3/105.d0
	rdelta(2)=dabs(rdeldu*rh) 
	call dcgmpr(fmatr(1,1),fmatr(1,1),fval12,ndim,ndim,ndim)
	l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
	pval(i,j)=-(fuval(i,j)-fuval(j,i))/8.d0-fval12(i,j)/24.d0+
     +  fval20c(i,j)/8.d0
	enddo
	enddo
      call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdelv=rpmax3/105.d0
	rdelta(3)=dabs(rdelv/rh)
	l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
	pval(i,j)=fmatr(l+i,3)/2.d0-(tval(i,j)-tval(j,i))/120.d0
	enddo
	enddo
      call dcmlmx(pval,rpmax1,ndim,ndim,ndim)
	l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
	pval(i,j)=-2.5d0*fmatr(l+i,3)-(fuval(i,j)-fuval(j,i))/8.d0-
     +	(tval(i,j)-tval(j,i))/120.d0-(fvval(i,j)-fvval(j,i))/96.d0
	enddo
	enddo
      call dcmlmx(pval,rpmax2,ndim,ndim,ndim)
	l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
	pval(i,j)=(fuval(i,j)-fuval(j,i))/2.d0-(tval(i,j)+tval(j,i))/
     +  4.d0-(tval(i,j)-tval(j,i))/120.d0+(fvval(i,j)-fvval(j,i))/
     +  96.d0
	enddo
	enddo
      call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdeldv=rpmax1/3.d0+rpmax2/15.d0+rpmax3/105.d0
	rdelta(4)=dabs(rdeldv)
	repsil=0.d0
	do i=1,4
	if(rdelta(i).gt.repsil) repsil=rdelta(i)
	enddo
c     construct new stepsize: 
	rhnew=rh*(rtol/repsil)**(1.d0/7.d0)
	numiter=numiter+1
	rdifer=dabs(rhnew/rh-1.d0)
	if(rdifer.gt.0.1d0) then
	rh=rhnew
	if(numiter.gt.20) then
	rh=0.1*rh
	go to 509
	endif
	go to 505
	else
509	npn=npn+1
	if(npn.gt.npoints) pause
	rnext=r+rh
	if(rnext.ge.radfin) then
	rpart(npn)=radfin
	rh=radfin-r
	go to 506
	endif
	rpart(npn)=r+rh
	go to 507
	endif
506	continue
	npart=npn
**************************************************
c     stepsize selection to rfin
	npn=npart
	radfin=rfin
	rh=rdx
407	r=rpart(npn)
	rnext=r+rh
	if(rnext.ge.radfin) rh=radfin-r
405	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	do k=0,ndeg
      call quadg16(rxl,rxu,fctu,cofk)
	if(k.eq.0) then
	pu(l+i,k)=(2*k+1)*cofk
	else
	pu(l+i,k)=(2*k+1)*rh**2*cofk
	endif
      call quadg16(rxl,rxu,fctv,cofk)
	if(k.eq.0) then
	pv(l+i,k)=(2*k+1)*cofk
	else
	pv(l+i,k)=(2*k+1)*rh**2*cofk
	endif
	enddo
	enddo
	enddo

      l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
	amatr(i,j)=pu(l+i,0)+depth*pv(l+i,0)
	enddo
	enddo

c     compute {\bf F}^D_0 = fmatr(,0)
	call dchouse(amatr,ndim,ndim,dval,zval,wval)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	if(i.eq.j) then
	fmatr(l+i,0)=dval(i)*rh**2
	else
	fmatr(l+i,0)=0.d0
	endif
	enddo
	enddo

c     compute {\bf F}^D_i ,i>0
	call dcgmtr(zval,tval,ndim,ndim)
	do k=1,ndeg
	call dcgmpr(tval,pu(1,k),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fuval,ndim,ndim,ndim)
	call dcgmpr(tval,pv(1,k),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fvval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	fmatr(l+i,k)=fuval(i,j)+depth*fvval(i,j)
	enddo
	enddo
	enddo

c       f1*f2 in tval
	call dcgmpr(fmatr(1,1),fmatr(1,2),tval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=-0.5d0*fmatr(l+i,3)-(tval(i,j)-tval(j,i))/120.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax1,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,3),fmatr(1,0),fuval,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,1),fmatr(1,0),fvval,ndim,ndim,ndim)
	do j=1,ndim
	do i=1,ndim
	fval10(i,j)=fvval(i,j)-fvval(j,i)
	enddo
	enddo
	call dcgmpr(fval10,fmatr(1,0),fvval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=2.5d0*fmatr(l+i,3)+(fuval(i,j)-fuval(j,i))/8.d0-
     +  (tval(i,j)-tval(j,i))/120.d0-(fvval(i,j)-fvval(j,i))/96.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax2,ndim,ndim,ndim)
	
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=-0.75d0*(fuval(i,j)-fuval(j,i))+(tval(i,j)+
     +  tval(j,i))/4.d0-(tval(i,j)-tval(j,i))/40.d0-
     +  (fvval(i,j)-fvval(j,i))/32.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdelu=rpmax1/3.d0+rpmax2/15.d0+rpmax3/105.d0
	rdelta(1)=dabs(rdelu)
	call dcgmpr(fmatr(1,2),fmatr(1,2),fval22,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,2),fmatr(1,0),fval20,ndim,ndim,ndim)
	do j=1,ndim
	do i=1,ndim
	fval20c(i,j)=fval20(i,j)-fval20(j,i)
	enddo
	enddo
	call dcgmpr(fval20c,fmatr(1,0),fval200,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,1),fval10,fval11,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=(fuval(i,j)-fuval(j,i))/8.d0-fval22(i,j)/40.d0+
     +  (fval200(i,j)-fval200(j,i))/480.d0+(fval11(i,j)-fval11(j,i))/
     +  480.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax1,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=-5.d0*(fuval(i,j)-fuval(j,i))/8.d0-3.d0*fval22(i,j)/
     +  20.d0+(fval200(i,j)-fval200(j,i))/30.d0+
     +  (11.d0*fval11(i,j)-fval11(j,i))/480.d0
	enddo
	enddo
	call dcmlmx(pval,rpmax2,ndim,ndim,ndim)
	l=-ndim
        do j=1,ndim
        l=l+ndim
        do i=1,ndim
        pval(i,j)=57.d0*fval22(i,j)/40.d0+11.d0*fval11(j,i)/240.d0+
     +  fval11(i,j)/10.d0-9.d0*(fval200(i,j)-fval200(j,i))/160.d0
        enddo
        enddo
        call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdeldu=rpmax1/3.d0+rpmax2/15.d0+rpmax3/105.d0
	rdelta(2)=dabs(rdeldu*rh) 
	call dcgmpr(fmatr(1,1),fmatr(1,1),fval12,ndim,ndim,ndim)
	l=-ndim
        do j=1,ndim
        l=l+ndim
        do i=1,ndim
	pval(i,j)=-(fuval(i,j)-fuval(j,i))/8.d0-fval12(i,j)/24.d0+
     +  fval20c(i,j)/8.d0
	enddo
	enddo
        call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdelv=rpmax3/105.d0
	rdelta(3)=dabs(rdelv/rh)
	l=-ndim
        do j=1,ndim
        l=l+ndim
        do i=1,ndim
	pval(i,j)=fmatr(l+i,3)/2.d0-(tval(i,j)-tval(j,i))/120.d0
	enddo
	enddo
        call dcmlmx(pval,rpmax1,ndim,ndim,ndim)
	l=-ndim
        do j=1,ndim
        l=l+ndim
        do i=1,ndim
	pval(i,j)=-2.5d0*fmatr(l+i,3)-(fuval(i,j)-fuval(j,i))/8.d0-
     +	(tval(i,j)-tval(j,i))/120.d0-(fvval(i,j)-fvval(j,i))/96.d0
	enddo
	enddo
        call dcmlmx(pval,rpmax2,ndim,ndim,ndim)
	l=-ndim
        do j=1,ndim
        l=l+ndim
        do i=1,ndim
	pval(i,j)=(fuval(i,j)-fuval(j,i))/2.d0-(tval(i,j)+tval(j,i))/
     +  4.d0-(tval(i,j)-tval(j,i))/120.d0+(fvval(i,j)-fvval(j,i))/
     +  96.d0
	enddo
	enddo
        call dcmlmx(pval,rpmax3,ndim,ndim,ndim)
	rdeldv=rpmax1/3.d0+rpmax2/15.d0+rpmax3/105.d0
	rdelta(4)=dabs(rdeldv)
	repsil=0.d0
	do i=1,4
	if(rdelta(i).gt.repsil) repsil=rdelta(i)
	enddo
	rhnew=rh*(rtol/repsil)**(1.d0/7.d0)
	rdifer=dabs(rhnew/rh-1.d0)
c	write(*,*) 'eps =',repsil,' hnew =',rhnew,' dif =',rdifer
	if(rdifer.gt.0.1d0) then
	rh=rhnew
	go to 405
	else
	npn=npn+1
	if(npn.gt.npoints) pause
	rnext=r+rh
	if(rnext.ge.radfin) then
	rpart(npn)=radfin
	rh=radfin-r
	go to 406
	endif
	rpart(npn)=r+rh
	go to 407
	endif
406	continue
	npart=npn
c
c     mesh is now constructed
c     compute data associated to mesh
c
c     compute matching point:
	do i=1,npart
	if(rpart(i).gt.rmatc) then
	kmatc=i
	rmatch=rpart(i)
	go to 408
	endif
	enddo
408	continue
	write(*,*) '  number of intervals in mesh =',npart,
     c' matching point =',rmatch
	do kk=1,npart
	r=rpart(kk)
      l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
      pnuc(kk,l+i)=potnuc(i,j,r)
      enddo
      enddo
      enddo

	do m=1,npart-1
	r=rpart(m)
	rh=rpart(m+1)-r
	rhp2=rh*rh
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	do k=0,ndeg-1
      call quadg16(rxl,rxu,fctu,cofku)
      call quadg16(rxl,rxu,fctv,cofkv)
	if(k.eq.0) then
	pu0(l+i,m)=(2*k+1)*cofku
	pv0(l+i,m)=(2*k+1)*cofkv
	endif
	if(k.eq.1) then
	pu1(l+i,m)=(2*k+1)*rhp2*cofku
	pv1(l+i,m)=(2*k+1)*rhp2*cofkv
	endif
	if(k.eq.2) then
	pu2(l+i,m)=(2*k+1)*rhp2*cofku
	pv2(l+i,m)=(2*k+1)*rhp2*cofkv
	endif
	enddo
	enddo
	enddo
	enddo

c     check, approximations of the functions U and V are ok:
c	do m1=2,npart
c	r=rpart(m1)
c	rd=rpart(m1)-rpart(m1-1)
c	rh=rd
c	m2=m1-1
c	rr=rd/rh
c	l=-ndim
c	do j=1,ndim
c	l=l+ndim
c	do i=1,ndim
c     valpu=potu(i,j,r)
c      valpv=potnuc(i,j,r)
c      rpol1=2.d0*rr-1.d0
c	rpol2=1.d0+rr*(-6.d0+6.d0*rr)
c	valpua=pu0(l+i,m2)+pu1(l+i,m2)*rpol1/rh**2+pu2(l+i,m2)*
c     +  rpol2/rh**2
c	valpva=pv0(l+i,m2)+pv1(l+i,m2)*rpol1/rh**2+pv2(l+i,m2)*
c     +  rpol2/rh**2
c	if(valpu.ne.(0.d0,0.d0)) then
c	releru=cdabs((valpu-valpua)/valpu)
c	else
c	releru=cdabs(valpu-valpua)
c	endif
c	if(valpv.ne.(0.d0,0.d0)) then
c	relerv=cdabs((valpv-valpva)/valpv)
c	else
c	relerv=cdabs(valpv-valpva)
c	endif
c	if(i.eq.1.and.j.eq.1) then	
c      write(11,2212) r,dreal(valpu),dreal(valpua),releru
c      write(31,2212) r,dreal(valpv),dreal(valpva),relerv
c2212	format(1x,f8.4,1x,2d24.14,2x,d18.10)
c      endif
c	if(i.eq.2.and.j.eq.1) then	
c        write(21,2212) r,dreal(valpu),dreal(valpua),releru
c        write(41,2212) r,dreal(valpv),dreal(valpva),relerv
c        endif
c	if(i.eq.1.and.j.eq.2) then	
c        write(12,2212) r,dreal(valpu),dreal(valpua),releru
c        write(32,2212) r,dreal(valpv),dreal(valpva),relerv
c        endif
c	if(i.eq.2.and.j.eq.2) then	
c        write(22,2212) r,dreal(valpu),dreal(valpua),releru
c        write(42,2212) r,dreal(valpv),dreal(valpva),relerv
c        endif
c	enddo
c	enddo
c	enddo
c	close(11)
c	close(12)
c	close(21)
c	close(22)
c	close(31)
c	close(32)
c	close(41)
c	close(42)
c
c       open(11,file='petu.11')
c       open(12,file='petu.12')
c       open(21,file='petu.21')
c       open(22,file='petu.22')
c       open(31,file='petv.11')
c       open(32,file='petv.12')
c       open(41,file='petv.21')
c       open(42,file='petv.22')
c	rstep=0.01d0
c	jmax=(rpart(npart)-rpart(1))/rstep
c	do j1=0,jmax
c	r=rpart(1)+j1*rstep
c	l=-ndim
c	do j=1,ndim
c	l=l+ndim
c	do i=1,ndim
c        valpu=potu(i,j,r)
c        valpv=potnuc(i,j,r)
c	if(i.eq.1.and.j.eq.1) then	
c        write(11,2212) r,dreal(valpu)
c        write(31,2212) r,dreal(valpv)
c        endif
c	if(i.eq.2.and.j.eq.1) then	
c        write(21,2212) r,dreal(valpu)
c        write(41,2212) r,dreal(valpv)
c        endif
c	if(i.eq.1.and.j.eq.2) then	
c        write(12,2212) r,dreal(valpu)
c        write(32,2212) r,dreal(valpv)
c        endif
c	if(i.eq.2.and.j.eq.2) then	
c        write(22,2212) r,dreal(valpu)
c        write(42,2212) r,dreal(valpv)
c        endif
c	enddo
c	enddo
c	enddo
c	close(11)
c	close(12)
c	close(21)
c	close(22)
c	close(31)
c	close(32)
c	close(41)
c	close(42)

c 
c    Compute eigenvalue
c
	rrv0=dreal(depth)
	riv0=dimag(depth)
      iter=0
	itermax=200
c     'corr' propagates left- and right hand solution
c     and computes correction for eigenvalue
52    call corr(rrv0,riv0,del)
      write(*,11) iter,dreal(depth),dreal(del)
11	format(1x,'it=',i3,2x,'V_0=',g21.12,2x,'dv=',1pd12.5)
      write(3,11) iter,dreal(depth),dreal(del)
	iter=iter+1
      depth=depth+del
	if(cdabs(del).lt.rtol) go to 53
	if(iter.gt.itermax) go to 53
	rrv0=dreal(depth)
	riv0=dimag(depth)
	go to 52
53	write(3,15) depth
	write(*,15) depth
15	format(1x,'resulting depth =  ',2g22.12)
 	ranc=0.d0
	do kk=1,npart-1
	rquadl=0.d0
	rquadr=0.d0
      do i=1,ndim
	rquadl=rquadl+tsol(kk,i)*dconjg(tsol(kk,i))
	rquadr=rquadr+tsol(kk+1,i)*dconjg(tsol(kk+1,i))
      enddo
	rquadp=0.5d0*(rquadl+rquadr)*(rpart(kk+1)-rpart(kk))
	ranc=ranc+rquadp
	enddo

      rcn=1.d0/dsqrt(ranc)
      do kk=1,npart
      do i=1,ndim
      tsol(kk,i)=tsol(kk,i)*rcn
      tdsol(kk,i)=tdsol(kk,i)*rcn
      enddo
      enddo
	bsol=0.d0
	bdsol=1.d0
c  wp** - files containing the eigenfunctions
      name(1:2)='wp'
      do n=1,ndim
      name(3:4)=numb(n)
      open(18,file=name)
c1515 format(f9.5,1x,2g17.9,1x,2g17.9)
1515  format(1x,f9.5,3x,g17.9,2x,g17.9)
	r=0.d0
      write(18,1515) r,dreal(bsol),dreal(bdsol)
      do k=1,npart
      r=rpart(k)
      write(18,1515) r,dreal(tsol(k,n)),dreal(tdsol(k,n))
      enddo
      close(18)
      enddo
	stop
      end

c******************************************************************	     
        subroutine quadg16(xl,xu,fct,y) 
c   this subroutine performs the gaussian 
c   quadrature by a 16-point rule.
c     input parameters:
c     xl = the lower limit of the integration. 
c     xu = the upper limit of the integration. 
c     fct= the function (must be declared external in calling program)
c     output parameter: 
c     y = the integral
c
	implicit double precision (a-h,o-z)
	complex*16 s,fct,y
	parameter (n=8)
	dimension x(n),w(n)
	data x/
     +  0.095012509837637440185d0,0.281603550779258913230d0,
     +  0.458016777657227386342d0,0.617876244402643748447d0,
     +  0.755404408355003033895d0,0.865631202387831743880d0,
     +  0.944575023073232576078d0,0.989400934991649932596d0/
        data w/
     +  0.189450610455068496285d0,0.182603415044923588867d0,
     +  0.169156519395002538189d0,0.149595988816576732081d0,
     +  0.124628971255533872052d0,0.095158511682492784810d0,
     +  0.062253523938647892863d0,0.027152459411754094852d0/
c
c executable statements
c
 	a=0.5d0*(xu-xl)
        b=0.5d0*(xu+xl)    
	s=0.d0
	do i=1,n
	s=s+w(i)*(fct(-a*x(i)+b)+fct(a*x(i)+b))
	enddo
	y=a*s
     	return
      end

	function fctu(t)
	implicit double precision (a-h,o-z)
	double complex fctu,potu
	common /legend/ x,h,k
	common /indexes/ i,j
c     poleg= shifted Legendre polynomial
	if(k.eq.0) then
	poleg=1.d0
	endif
	if(k.eq.1) then
	poleg=2.d0*t-1.d0
	endif
	if(k.eq.2) then
	poleg=1.d0+t*(-6.d0+6.d0*t)
	endif
	if(k.eq.3) then
	poleg=-1.d0+t*(12.d0+t*(-30.d0+t*20.d0))
	endif
	if(k.eq.4) then
	poleg=1.d0+t*(-20.d0+t*(90.d0+t*(-140.d0+t*70.d0)))
	endif
	if(k.eq.5) then
	poleg=-1.d0+t*(30.d0+t*(-210.d0+t*(560.d0+t*(-630.d0+
     #  t*252.d0))))
	endif
	if(k.eq.6) then
	z=2.d0*t-1.d0
	poleg=(231.d0*z**6-315.d0*z**4+105.d0*z**2-5.d0)/16.d0
	endif
	fctu=poleg*potu(i,j,x+t*h)
	return
	end

	function fctv(t)
	implicit double precision (a-h,o-z)
	double complex fctv,potnuc
	common /legend/ x,h,k
	common /indexes/ i,j
	if(k.eq.0) then
	poleg=1.d0
	endif
	if(k.eq.1) then
	poleg=2.d0*t-1.d0
	endif
	if(k.eq.2) then
	poleg=1.d0+t*(-6.d0+6.d0*t)
	endif
	if(k.eq.3) then
	poleg=-1.d0+t*(12.d0+t*(-30.d0+t*20.d0))
	endif
	if(k.eq.4) then
	poleg=1.d0+t*(-20.d0+t*(90.d0+t*(-140.d0+t*70.d0)))
	endif
	if(k.eq.5) then
	poleg=-1.d0+t*(30.d0+t*(-210.d0+t*(560.d0+t*(-630.d0+
     #  t*252.d0))))
	endif
	if(k.eq.6) then
	z=2.d0*t-1.d0
	poleg=(231.d0*z**6-315.d0*z**4+105.d0*z**2-5.d0)/16.d0
	endif
	fctv=poleg*potnuc(i,j,x+t*h)
	return
	end

c..................................................

	subroutine dchouse(a,n,np,d,z,e)
c   computes the eigenvalues and eigenvectors of a complex symmetric 
c   n by n matrix a, stored in an np by np physical array 
c   by householder reduction and q-l algorithm with implicit shifts.
c   d is a vector of length np. on output, it returns the
c   eigenvalues.
c   e is a working vector of length n.
c   the matrix z (n by n matrix stored in np by np array) contains
c   on output the eigenvectors. the k-th column of z
c   returns the normalized eigenvector corresponding to d(k).
	implicit double complex(a-h,o-z)
	logical loop1,loop2
	dimension a(np,np),z(np,np),d(np),e(np)
	data zero/(0.d0,0.d0)/,one/(1.d0,0.d0)/
	do  i=1,n
	  do  j=1,n
	    z(i,j)=a(i,j)
	  enddo
	enddo
	if(n.gt.1) then
	  do i=n,2,-1
	    l=i-1
	    h=zero
	    scale=zero
	    if(l.gt.1) then
	      do k=1,l
	        scale=scale+cdabs(z(i,k))
	      enddo
	      if(scale.eq.zero) then
	        e(i)=z(i,l)
	      else
	        do k=1,l
	          z(i,k)=z(i,k)/scale
	          h=h+z(i,k)**2
	        enddo
	        f=z(i,l)
	        sqh=cdsqrt(h)
	        if(dreal(f).lt.0.d0) sqh=-sqh
	        g=-sqh	
	        e(i)=scale*g
	        h=h-f*g
	        z(i,l)=f-g
	        f=zero
	        do  j=1,l
	          z(j,i)=z(i,j)/h
	          g=zero
	          do  k=1,j
	            g=g+z(j,k)*z(i,k)
	          enddo
	          if(l.gt.j) then
	            do  k=j+1,l
	              g=g+z(k,j)*z(i,k)
	            enddo
	          endif
	          e(j)=g/h
	          f=f+e(j)*z(i,j)
	        enddo
	        hh=f/(h+h)
	        do j=1,l
	          f=z(i,j)
	          g=e(j)-hh*f
	          e(j)=g
	          do  k=1,j
	            z(j,k)=z(j,k)-f*e(k)-g*z(i,k)
                    enddo
	        enddo
	      endif
	    else
	      e(i)=z(i,l)
	    endif
	    d(i)=h
	  enddo
          endif
	d(1)=zero
	e(1)=zero
	do  i=1,n
	  l=i-1
	  if(d(i).ne.zero) then
	    do  j=1,l
	      g=zero
	      do  k=1,l
	        g=g+z(i,k)*z(k,j)
	      enddo
	      do  k=1,l
	        z(k,j)=z(k,j)-g*z(k,i)
	      enddo
	    enddo
 	  endif
	  d(i)=z(i,i)
	  z(i,i)=one
	  if(l.ge.1) then
	    do  j=1,l
	      z(i,j)=zero
	      z(j,i)=zero
	    enddo
	  endif
	enddo
	if(n.gt.1) then
	  do  i=2,n
	    e(i-1)=e(i)
	  enddo
	  e(n)=zero
	  do l=1,n
	    iter=0
	    loop2=.true.
	    do while(loop2)
	      m=l-1
	      id=0
	      loop1=.true.
	      do while(loop1)
	        m=m+1
	        dd=cdabs(d(m))+cdabs(d(m+1))
	        if(cdabs(e(m))+dd.eq.dd)then
	          loop1=.false.
	          id=1
 	        endif
	        if(m.eq.n-1)loop1=.false.
	      enddo
	      if(m.eq.n-1.and.id.eq.0)m=n
	      if(m.ne.l) then
	        if(iter.eq.30) then
	          write(*,*) 'too  many  iterations'
	          return
	        endif
	        iter=iter+1
	        g=(d(l+1)-d(l))/(2.d0*e(l))
	        r=cdsqrt(g**2+one)
	        if(dreal(g).lt.0.d0) r=-r
	        g=d(m)-d(l)+e(l)/(g+r)
	        s=one
	        c=one
	        p=zero
	        do i=m-1,l,-1
	          f=s*e(i)
	          b=c*e(i)
	          if(cdabs(f).ge.cdabs(g)) then
	            c=g/f
	            r=cdsqrt(c**2+one)
	            e(i+1)=f*r
	            s=one/r
	            c=c*s
	          else
	            s=f/g
	            r=cdsqrt(s**2+one)
	            e(i+1)=g*r
	            c=one/r
	            s=s*c
	          endif
	          g=d(i+1)-p
	          r=(d(i)-g)*s+2.d0*c*b
	          p=s*r
	          d(i+1)=g+p
	          g=c*r-b
	          do k=1,n
	            f=z(k,i+1)
	            z(k,i+1)=s*z(k,i)+c*f
	            z(k,i)=c*z(k,i)-s*f
	          enddo
                  enddo
	        d(l)=d(l)-p
	        e(l)=g
	        e(m)=zero
	      else
	        loop2=.false.
	      endif
	    enddo
	  enddo
	endif
	return
	end
c     ..................................................................                
      subroutine dcgmpr(a,b,w,n,m,l)
c  this multiplies two double complex matrices, w=a*b. 
c        description of parameters                                      
c           a - name of first input matrix                              
c           b - name of second input matrix                             
c           w - name of output matrix                                   
c           n - number of rows in a                                     
c           m - number of columns in a and rows in b                    
c           l - number of columns in b                                  
c                                                                       
c        remarks                                                        
c           all matrices must be stored as general matrices             
c           matrix w cannot be in the same location as matrix a or b         
c           number of columns of matrix a must be equal to number of row
c           of matrix b                                
      implicit double complex (a-h,o-q,s-z)                                 
      dimension a(*),b(*),w(*)                                                                                       
      czero=dcmplx(0.d0,0.d0)                                           
      ir=0                                                              
      ik=-m                                                             
      do k=1,l                                                       
        ik=ik+m                                                           
        do j=1,n                                                       
          ir=ir+1                                                           
          ji=j-n                                                            
          ib=ik                                                             
          w(ir)=czero                                                       
          do i=1,m                                                       
            ji=ji+n                                                           
            ib=ib+1                                                           
            w(ir)=w(ir)+a(ji)*b(ib)                                           
	enddo
       enddo
      enddo
      return                                                            
      end                                                               

c                                                                       
c     ..................................................................            
      subroutine dcgmtr(a,w,n,m)
c  this transposes a double complex general matrix         
c        description of parameters                                      
c           a - name of matrix to be transposed                         
c           w - name of resultant matrix                                
c           n - number of rows of a and columns of w                    
c           m - number of columns of a and rows of w                    
c                                                                       
c        remarks                                                        
c           matrix w cannot be in the same location as matrix a         
c           matrices a and w must be stored as general matrices         
c                                                           
      implicit double complex (a-h,o-q,s-z)	                                
      dimension a(*),w(*)                   
      ir=0                                                              
      do i=1,n                                                       
        ij=i-n                                                            
        do j=1,m                                                       
          ij=ij+n                                                           
          ir=ir+1                                                           
          w(ir)=a(ij)                
        enddo
      enddo	
      return                                                            
      end                                                               
c                                                                       
c        ...............................................................                
      subroutine dcmcpy(a,w,n,m,ms)      
c  this copies the input double complex matrix a into w
c        description of parameters                                      
c           a - name of input matrix                                    
c           w - name of output matrix                                   
c           n - number of rows in a or w                                
c           m - number of columns in a or w                             
c           ms  - one digit number for storage mode of matrix a (and w) 
c                  0 - general                                          
c                  1 - symmetric                                        
c                  2 - diagonal                         

      implicit double complex (a-h,o-q,s-z)	                                
      dimension a(*),w(*)                                               
c                                                                       
c        compute vector length, it                                      
c                                                                       
      call lloc(n,m,it,n,m,ms)                                          
c                                                                       
c        copy matrix                                                    
c                                                                       
      do  i=1,it                                                       
        w(i)=a(i)              
      enddo
      return                                                            
      end                                                               
c     ..................................................................
      subroutine lloc(i,j,ir,n,m,ms)                                    
c   this computes a vector subscript for an element in a matrix of    
c           specified storage mode                                      
c                                                                       
c        description of parameters                                      
c           i   - row number of element                                 
c           j   - column number  of element                             
c           ir  - resultant vector subscript                            
c           n   - number of rows in matrix                              
c           m   - number of columns in matrix                           
c           ms  - one digit number for storage mode of matrix           
c                  0 - general                                          
c                  1 - symmetric                                        
c                  2 - diagonal                                         
c                                                                       
c        method                                                         
c           ms=0   subscript is computed for a matrix with n*m elements 
c                  in storage (general matrix)                          
c           ms=1   subscript is computed for a matrix with n*(n+1)/2 in 
c                  storage (upper triangle of symmetric matrix). if     
c                  element is in lower triangular portion, subscript is 
c                  corresponding element in upper triangle.             
c           ms=2   subscript is computed for a matrix with n elements   
c                  in storage (diagonal elements of diagonal matrix).   
c                  if element is not on diagonal (and therefore not in  
c                  storage), ir is set to zero.                         
c                                                                      
c                                                                       
      ix=i                                                              
      jx=j 
      if(ms.eq.0)ir=n*(jx-1)+ix                                                   
      if(ms.eq.1)then
        if(ix.lt.jx)then
          ir=ix+(jx*jx-jx)/2                                               
        else
          ir=jx+(ix*ix-ix)/2                                               
        endif
      endif
      if(ms.eq.2)then
        if(ix.eq.jx)then
	ir=ix
        else
	ir=0
        endif
      endif
      return
      end	

c.................................................
                                                                        
      subroutine dcmlmx(a,ramax,n,m,np)                                 
c this computes ramax, the absolute value of the maximal matrix element
c of matrix a
      implicit double complex (a-h,o-q,s-z)                             
      implicit double precision (r)
      dimension a(np,*)                                                 
      ramax=0.d0                                                        
      do i=1,n                                                        
        do j=1,m                                                        
          ra=cdabs(a(i,j))                                                  
          if(ra.gt.ramax) ramax=ra                                          
        enddo
      enddo
      return                                                            
      end                                                               
                                                                        
c..................................................

                                                                        
      subroutine dcunitm(a,n,m,np)                                      
c  this generates the unit matrix a
      implicit double complex (a-h,o-q,s-z)	                                
      dimension a(np,*)                                                 
      call dczerom(a,n,m,np)                                            
      n0=n                                                              
      if(n.gt.m) n0=m                                                   
      do  i=1,n0                                                       
        a(i,i)=dcmplx(1.0d0,0.d0)                                        
      enddo
      return                                                            
      end                                                               
                                                                        
c..................................................                                                                        

      subroutine dczerom(a,n,m,np)                                      
c  this generates the zero matrix a
      implicit double complex (a-h,o-q,s-z)	                                
      dimension a(np,*)                  
      zero=dcmplx(0.d0,0.d0)                                            
      do 1 i=1,n                                                        
      do 1 j=1,m                                                        
      a(i,j)=zero                                                       
    1 continue                                                          
      return                                                            
      end                                          

c.................................................                     
                                                                                                                                                                                                                       
      subroutine dccopy(a,b,n,m,np)                                     
c  this copies matrix a into b
      implicit double complex (a-h,o-q,s-z)	                                
      dimension a(np,*),b(np,*)                                         
      do 1 i=1,n                                                        
      do 1 j=1,m                                                        
      b(i,j)=a(i,j)                                                     
    1 continue                                                          
      return                                                            
      end                                                               
                                                  
c..................................................
                                                                                                                                                                      
      subroutine dcmsum(a,b,c,x,n,m,np)                                 
c  this performs the matrix operation c=a+xb, where x is a complex number
      implicit double complex (a-h,o-q,s-z)	                                
      dimension a(np,*),b(np,*),c(np,*)                                 
      do 1 i=1,n                                                        
      do 1 j=1,m                                                        
      c(i,j)=a(i,j)+x*b(i,j)                                            
    1 continue                                                          
      return                                                            
      end                                                               
c....................................................
	subroutine renorm(n,y,yp,yk,ykp,cw1,cw2,cw3,cw4,back,indx)
c
c   the input n by n matrices y,yp,yk and yp, stored as vectors, are converted
c   by using the permuted lu representation of y, i.e. y=p*l*u where  p is the n 
c   by n permutation matrix. in output the 
c   mentioned matrices are replaced by y*u^-1=p*l, yp*u^-1, yk*u^-1 and ykp*u^-1.
c   the output n by n matrix back, also stored as a vector, helps reconverting   
c   these data in the original representation. output back = u times input 
c   back.	
c   version january 2001
c
	implicit  double complex (a-h,o-q,s-z)
	implicit double precision(r)
	dimension y(*),yp(*),yk(*),ykp(*),
     !    cw1(*),cw2(*),cw3(*),cw4(*),back(*),indx(*)
	data czero /(0.d0,0.d0)/, cone /(1.d0,0.d0)/
	call dccopy(y,cw1,n,n,n)
      call ludcmpn(cw1,n,n,indx,d,ierr)
	call dccopy(cw1,cw3,n,n,n)
	do i=1,n
	  do j=1,n
	    ivec=i+(j-1)*n
	    if(j.gt.i)cw1(ivec)=czero
	    if(j.eq.i)cw1(ivec)=cone
	    if(j.lt.i)cw3(ivec)=czero
 	  enddo
	enddo
c l is now in vector cw1. matrix pl is constructed in cw2 by permuting the lines 
c of l as it follows:
	do i=1,n
	  do m=1,n
	    cw2(n*(m-1)+indx(i))=cw1(n*(m-1)+i)
	  enddo
	enddo
c  at this moment the pl and u matrices are in vectors cw2 and cw3, respectively.
	call dccopy(cw3,cw4,n,n,n)
	call dcludcmp(cw4,n,n,indx,d,ierr)
      call dcunitm(cw1,n,n,n)                                      
	do j=1,n
	 call dclubksb(cw4,n,n,indx,cw1((j-1)*n+1))
	enddo
c  at this moment cw1 is u^-1.
	call dccopy(cw2,y,n,n,n)
	call dcgmpr(yp,cw1,cw4,n,n,n)
	call dccopy(cw4,yp,n,n,n)
	call dcgmpr(yk,cw1,cw4,n,n,n)
	call dccopy(cw4,yk,n,n,n)
	call dcgmpr(ykp,cw1,cw4,n,n,n)
	call dccopy(cw4,ykp,n,n,n)
	call dcgmpr(cw3,back,cw4,n,n,n)
	call dccopy(cw4,back,n,n,n)
	return
	end
c..................................................

      subroutine dcludcmp(a,n,np,indx,d,ierror)
c  this is the double complex version of the homonyme sbr. from numerical
c  recipes
      implicit double complex(a-h,o-q,s-z)
      implicit double precision (r)
      parameter (tiny=(1.0d-30,1.0d-30))
      dimension a(np,np),indx(n),rvv(n)
      d=dcmplx(1.d0,0.d0)
	ierror=0
      do 12 i=1,n
        raamax=0.d0
        do 11 j=1,n
          if (cdabs(a(i,j)).gt.raamax) raamax=cdabs(a(i,j))
11      continue
        if (raamax.eq.0.d0)then
	write (*,*) 'singular matrix in dcludcmp'
	ierror=1
	return
	endif
        rvv(i)=1.d0/raamax
12    continue
      do 19 j=1,n
        if (j.gt.1) then
          do 14 i=1,j-1
            sum=a(i,j)
            if (i.gt.1)then
              do 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
13            continue
              a(i,j)=sum
            endif
14        continue
        endif
        raamax=0.d0
        do 16 i=j,n
          sum=a(i,j)
          if (j.gt.1)then
            do 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
15          continue
            a(i,j)=sum
          endif
          rdum=rvv(i)*cdabs(sum)
          if (rdum.ge.raamax) then
            imax=i
            raamax=rdum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          rvv(imax)=rvv(j)
        endif
        indx(j)=imax
        if(j.ne.n)then
          if(cdabs(a(j,j)).eq.0.d0) a(j,j)=tiny
          dum=1.d0/a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      if(cdabs(a(n,n)).eq.0.d0) a(n,n)=tiny
      return
      end

c........................................................

      subroutine ludcmpn(a,n,np,nord,d,ierror)
c  this is a slightly modified double complex version of sbr. ludcmp from 
c  numerical recipes. the main difference is the argument nord, instead 
c  of indx.
c  
      implicit double complex(a-h,o-q,s-z)
      implicit double precision (r)
      parameter (tiny=(1.0d-30,1.0d-30))
      dimension a(np,np),rvv(n),nord(*)
      d=dcmplx(1.d0,0.d0)
	do j=1,n
	nord(j)=j
	enddo
	ierror=0
      do 12 i=1,n
        raamax=0.d0
        do 11 j=1,n
          if (abs(a(i,j)).gt.raamax) raamax=abs(a(i,j))
11      continue
        if (raamax.eq.0.d0)then
	write (*,*) 'singular matrix in ludcmpn'
	ierror=1
	return
	endif
        rvv(i)=1.d0/raamax
12    continue
      do 19 j=1,n
        if (j.gt.1) then
          do 14 i=1,j-1
            sum=a(i,j)
            if (i.gt.1)then
              do 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
13            continue
              a(i,j)=sum
            endif
14        continue
        endif
        raamax=0.d0
        do 16 i=j,n
          sum=a(i,j)
          if (j.gt.1)then
            do 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
15          continue
            a(i,j)=sum
          endif
          rdum=rvv(i)*abs(sum)
          if (rdum.ge.raamax) then
            imax=i
            raamax=rdum
          endif
16      continue
        if (j.ne.imax)then
	iord=nord(imax)
	jord=nord(j)
	nord(imax)=jord
	nord(j)=iord
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)	
            a(j,k)=dum
17        continue
          d=-d
          rvv(imax)=rvv(j)
        endif
c        indx(j)=imax
        if(j.ne.n)then
          if(abs(a(j,j)).eq.0.d0) a(j,j)=tiny
          dum=1.d0/a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      if(abs(a(n,n)).eq.0.d0) a(n,n)=tiny
      return
      end

c..................................................

      subroutine dclubksb(a,n,np,indx,b)
c  this is the double complex version of the homonyme sbr. from numerical
c  recipes
      implicit double complex(a-h,o-q,s-z)
      implicit double precision (r)
      dimension a(np,np),indx(n),b(n)
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.d0) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        if(i.lt.n)then
          do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
13        continue
        endif
        b(i)=sum/a(i,i)
14    continue
      return
      end

c...............................................
      subroutine prop(r,rh,int,if,sol,dsol)
c     propagates solution (sol) and its first 
c     derivative (dsol)
      include 'params'
      parameter(ndeg=3)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
	common /fituv0/ pu0(ndim2,kmax),pv0(ndim2,kmax)
	common /fituv1/ pu1(ndim2,kmax),pv1(ndim2,kmax)
	common /fituv2/ pu2(ndim2,kmax),pv2(ndim2,kmax)
	common /diagonal/ amatr(ndim,ndim),dval(ndim),zval(ndim,ndim),
     +  wval(ndim),tval(ndim,ndim),pval(ndim,ndim)
	common /fitmat/ fmatr(ndim2,0:ndeg),fuval(ndim,ndim),
     +  fvval(ndim,ndim),fval20c(ndim,ndim) 
	common /comat/ fval10(ndim,ndim),fval11(ndim,ndim),fval22(ndim,
     +  ndim),fval20(ndim,ndim),fval200(ndim,ndim),fval12(ndim,ndim)
	common /eigend/ depth
	dimension cu(ndim2,2),cud(ndim2,0:2),cv(ndim2,2),cvd(ndim2,2)
	dimension csi(ndim,ndim),eta(ndim2,0:2)
	dimension sol(ndim,ndim),dsol(ndim,ndim)
      l=-ndim
      do j=1,ndim
      l=l+ndim
      do i=1,ndim
	amatr(i,j)=pu0(l+i,int)+depth*pv0(l+i,int)
	enddo
	enddo
c     diagonalize F_0:
	call dchouse(amatr,ndim,ndim,dval,zval,wval)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	if(i.eq.j) then
	fmatr(l+i,0)=dval(i)*rh**2
	else
	fmatr(l+i,0)=0.d0
	endif
	enddo
	enddo

	call dcgmtr(zval,tval,ndim,ndim)
	do k=1,ndeg-1
	if(k.eq.1) then
	call dcgmpr(tval,pu1(1,int),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fuval,ndim,ndim,ndim)
	call dcgmpr(tval,pv1(1,int),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fvval,ndim,ndim,ndim)
	endif
	if(k.eq.2) then
	call dcgmpr(tval,pu2(1,int),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fuval,ndim,ndim,ndim)
	call dcgmpr(tval,pv2(1,int),pval,ndim,ndim,ndim)
	call dcgmpr(pval,zval,fvval,ndim,ndim,ndim)
	endif
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	fmatr(l+i,k)=fuval(i,j)+depth*fvval(i,j)
	enddo
	enddo
	enddo
	call dcgmpr(fmatr(1,1),fmatr(1,0),fval10,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
c     {\bf C}^{(u)}_1  
	cu(l+i,1)=-0.5d0*fmatr(l+i,1)+(fval10(i,j)-fval10(j,i))/24.d0
	enddo
	enddo
	call dcgmpr(fmatr(1,1),fmatr(1,1),fval12,ndim,ndim,ndim)
	call dcgmpr(fmatr(1,2),fmatr(1,0),fval20,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
c     {\bf C}^{(u)}_2  
	cu(l+i,2)=-fval12(i,j)/24.d0+(fval10(i,j)-fval10(j,i))/6.d0-
     +  (fval20(i,j)-fval20(j,i))/8.d0
	enddo
	enddo
	call dcgmpr(fmatr(1,1),fmatr(1,2),fuval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
c     {\bf C}^{(u')}_0  
	cud(l+i,0)=fmatr(l+i,2)/2.d0+(fval10(i,j)-fval10(j,i))/24.d0-
     +  (fuval(i,j)-fuval(j,i))/120.d0
	enddo
	enddo
	do j=1,ndim
	do i=1,ndim
	fval20c(i,j)=fval10(i,j)-fval10(j,i)
	enddo
	enddo
	call dcgmpr(fval20c,fmatr(1,0),fvval,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
c     {\bf C}^{(u')}_1  
	cud(l+i,1)=-1.5d0*fmatr(l+i,2)-fval12(i,j)/24.d0+
     +  fval20c(i,j)/8.d0-(fval20(i,j)-fval20(j,i))/8.d0-
     +  (fvval(i,j)-fvval(j,i))/96.d0
c     {\bf C}^{(u')}_2  
	cud(l+i,2)=-7.d0*fval12(i,j)/24.d0+(fuval(i,j)-fuval(j,i))/
     +  8.d0+3.d0*(fval20(i,j)-fval20(j,i))/8.d0-
     +  (fvval(i,j)-fvval(j,i))/32.d0
c     {\bf C}^{(v)}_2  
	cv(l+i,2)=-fmatr(l+i,2)/2.d0+fval20c(i,j)/24.d0
c     {\bf C}^{(v')}_1  
	cvd(l+i,1)=fmatr(l+i,1)/2.d0+fval20c(i,j)/24.d0
c     {\bf C}^{(v')}_2  
	cvd(l+i,2)=-fval12(i,j)/24.d0-fval20c(i,j)/12.d0+
     +  (fval20(i,j)-fval20(j,i))/8.d0
	enddo
	enddo

	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	if(i.ne.j) then
	csi(i,j)=0.d0
	do k=0,2
	eta(l+i,k)=0.d0
	enddo
	else
	z=dval(i)*rh**2
c     compute eta-functions:
	call cgebas(z,csival,eta0,eta1,eta2,eta3,eta4,eta5,eta6)
	csi(i,i)=csival
	eta(l+i,0)=eta0
	eta(l+i,1)=eta1
	eta(l+i,2)=eta2
	endif
	enddo
	enddo

	call dcgmpr(cu(1,1),eta(1,1),fval11,ndim,ndim,ndim)
	call dcgmpr(cu(1,2),eta(1,2),fval22,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	fuval(i,j)=csi(i,j)+fval11(i,j)+fval22(i,j)
	enddo
	enddo	
	call dcgmpr(cud(1,0),eta(1,0),fval10,ndim,ndim,ndim)
	call dcgmpr(cud(1,1),eta(1,1),fval11,ndim,ndim,ndim)
	call dcgmpr(cud(1,2),eta(1,2),fval12,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	if(i.ne.j) then
	zeta=0.d0
	else
	z=dval(i)*rh**2
	zeta=z*eta(l+i,0)
	endif
	tval(i,j)=zeta+fval10(i,j)+fval11(i,j)+fval12(i,j)
	tval(i,j)=tval(i,j)/rh
	enddo
	enddo	
	call dcgmpr(cv(1,2),eta(1,2),fval12,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	fvval(i,j)=(eta(l+i,0)+fval12(i,j))*rh
	enddo
	enddo
	call dcgmpr(cvd(1,1),eta(1,1),fval11,ndim,ndim,ndim)
	call dcgmpr(cvd(1,2),eta(1,2),fval12,ndim,ndim,ndim)
	l=-ndim
	do j=1,ndim
	l=l+ndim
	do i=1,ndim
	pval(i,j)=csi(i,j)+fval11(i,j)+fval12(i,j)
	enddo
	enddo

	call dcgmtr(zval,fval20c,ndim,ndim)
      call dcgmpr(zval,fuval,fval11,ndim,ndim,ndim)
      call dcgmpr(fval11,fval20c,fuval,ndim,ndim,ndim)
      call dcgmpr(zval,fvval,fval11,ndim,ndim,ndim)
      call dcgmpr(fval11,fval20c,fvval,ndim,ndim,ndim)
      call dcgmpr(zval,tval,fval11,ndim,ndim,ndim)
      call dcgmpr(fval11,fval20c,tval,ndim,ndim,ndim)
      call dcgmpr(zval,pval,fval11,ndim,ndim,ndim)
      call dcgmpr(fval11,fval20c,pval,ndim,ndim,ndim)

      if(if.eq.1) then
      call dcgmpr(fuval,sol,fval11,ndim,ndim,ndim)
      call dcgmpr(fvval,dsol,fval12,ndim,ndim,ndim)
	call dcgmpr(tval,sol,fval10,ndim,ndim,ndim)
	call dcgmpr(pval,dsol,fval20,ndim,ndim,ndim)
	do j=1,ndim
	do i=1,ndim
	sol(i,j)=fval11(i,j)+fval12(i,j)
	dsol(i,j)=fval10(i,j)+fval20(i,j)
	enddo
	enddo
	else
	call dcgmtr(pval,zval,ndim,ndim)
      call dcgmpr(zval,sol,fval11,ndim,ndim,ndim)
	call dcgmtr(fvval,zval,ndim,ndim)
      call dcgmpr(zval,dsol,fval12,ndim,ndim,ndim)
	call dcgmtr(tval,zval,ndim,ndim)
      call dcgmpr(zval,sol,fval10,ndim,ndim,ndim)
	call dcgmtr(fuval,zval,ndim,ndim)
      call dcgmpr(zval,dsol,fval20,ndim,ndim,ndim)
      do j=1,ndim
      do i=1,ndim
      sol(i,j)=fval11(i,j)-fval12(i,j)
      dsol(i,j)=-fval10(i,j)+fval20(i,j)
      enddo
      enddo
	endif
	return
	end

c..................................................
      subroutine cgebas(z,csi,eta0,eta1,eta2,eta3,eta4,eta5,eta6)
c  for given complex z, this computes the basic functions csi(z), eta0(z)
c  eta1(z),...,eta6(z).
c  this is version april 1990.
      implicit double complex(a-h,o-q,s-z)
      implicit double precision(r)
      rtz=1.5d0
      ru=cdabs(z)
      if(ru.ge.rtz) then
        squ=cdsqrt(z)*dcmplx(0.d0,1.d0)
        eta0=cdsin(squ)/squ
        csi=cdcos(squ)
        eta1=(csi-eta0)/z
        eta2=(eta0-3.d0*eta1)/z
        eta3=(eta1-5.d0*eta2)/z
        eta4=(eta2-7.d0*eta3)/z
        eta5=(eta3-9.d0*eta4)/z
        eta6=(eta4-11.d0*eta5)/z
        return
      endif 	
      eta6=(269583552.d4+z*(89861184.d3+z*(1321488.d3+z*(11592.d3+
     !z*(69.d3+z*(3.d2+z))))))/3643017329952.d5
      eta5=(14018344704.d2+z*(539167104.d2+z*(89861184.d1+z*(
     !880992.d1+z*(5796.d1+z*(276.d0+z))))))/14572069319808.d3
      eta4=z*eta6+11.d0*eta5
      eta3=z*eta5+9.d0*eta4
      eta2=z*eta4+7.d0*eta3
      eta1=z*eta3+5.d0*eta2
      eta0=z*eta2+3.d0*eta1
      csi=z*eta1+eta0
      return
      end


c******************************************************************
        subroutine cmat0
c
c preparation for the nuclear coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      parameter (nlevelmax=25)
      dimension a(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common /coulrad/ radc
      common/rott/beta2t,beta4t,ntarg
      common/dimen/nlevel
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/valpi/pi,hbar
      external cg
      frot2t=rt*beta2t
      frot4t=rt*beta4t
      do i=1,nlevelmax
      do j=1,nlevelmax
      a(i,j)=0.d0
      enddo
      enddo
      mdim=nlevel
      i=0
      do 32 it=0,ntarg
      i=i+1
      j=0
      do 42 jt=0,ntarg
      j=j+1  
      if(i.gt.j) then
      a(i,j)=a(j,i)
      goto 42
      endif
      c=frot2t*sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
     &/4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)
      c=c+frot4t*sqrt((2.d0*2*it+1.d0)*9.d0*(2.d0*2*jt+1.d0)
     &/4.d0/pi)*cg(2*it,0,4,0,2*jt,0)**2/(2.d0*2*jt+1.d0)
      a(i,j)=c 
42    continue
32    continue
      call mdiag(a,mdim)
      return
      end
c******************************************************************
        subroutine cmat(r,cpot)
c
c coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      complex*16 cpot,vncc,vn
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common /coulrad/ radc
      common/rott/beta2t,beta4t,ntarg
      common/dimen/nlevel
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/valpi/pi,hbar
        do 10 i=1,nlevelmax
        do 20 j=1,nlevelmax
        cpot(i,j)=0.d0
 20     continue
 10     continue
                                                                                
        ndim=nlevel
        i=0
        do 32 it=0,ntarg
        i=i+1
        j=0
        do 42 jt=0,ntarg
        j=j+1
        if(i.gt.j) then
        cpot(i,j)=cpot(j,i)
        goto 42
        endif
        c=0.d0
c nuclear coupling                                                                               
        do 50 k=1,ndim
        c=c+vncc(r,ev(k))*evec(i,k)*evec(j,k)
 50     continue
c coulomb coupling
        c=c+sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct2(r)
     &  +sqrt((2.d0*2*it+1.d0)*9.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,4,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct4(r)
        cpot(i,j)=c
 42     continue
 32     continue
        do 51 i=1,ndim
        cpot(i,i)=cpot(i,i)-vn(r)
 51     continue
	return
	end

c******************************************************************
        subroutine cmatu(r,cpot)
c
c coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      complex*16 cpot
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common /coulrad/ radc
      common/rott/beta2t,beta4t,ntarg
      common/dimen/nlevel
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/valpi/pi,hbar
        do 10 i=1,nlevelmax
        do 20 j=1,nlevelmax
        cpot(i,j)=0.d0
 20     continue
 10     continue
       ndim=nlevel
        i=0
        do 32 it=0,ntarg
        i=i+1
        j=0
        do 42 jt=0,ntarg
        j=j+1
        if(i.gt.j) then
        cpot(i,j)=cpot(j,i)
        goto 42
        endif
        c=0.d0
c coulomb coupling
        c=c+sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct2(r)
     &  +sqrt((2.d0*2*it+1.d0)*9.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,4,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct4(r)
        cpot(i,j)=c
 42     continue
 32     continue
        return
        end

c****************************************************************
        function vncc(r,rxt)
c
c woods-saxon potential for nuclear potential
c
c****************************************************************
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
        common /eigend/ depth 
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
                                                                                
        if(r.lt.50.d0) then
           rr=r-r0-rxt
        vncc=-depth/(1.d0+exp(rr/rdif))
        else
        vncc=0.d0
        endif
        return
        end
c****************************************************************
        function vn(r)
c
c woods-saxon potential for nuclear potential
c
c****************************************************************
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
        common /eigend/ depth 
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
                                                                                
        if(r.lt.50.d0) then
           rr=r-r0
        vn=-depth/(1.d0+exp(rr/rdif))
        else
        vn=0.d0
        endif
        return
        end

c******************************************************************
        subroutine cmatd(r,cpot)
c
c coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      complex*16 cpot,vnccd,vnd
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common /coulrad/ radc
      common/rott/beta2t,beta4t,ntarg
      common/dimen/nlevel
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/valpi/pi,hbar
        do 10 i=1,nlevelmax
        do 20 j=1,nlevelmax
        cpot(i,j)=0.d0
 20     continue
 10     continue
                                                                                
        ndim=nlevel
        i=0
        do 32 it=0,ntarg
        i=i+1
        j=0
        do 42 jt=0,ntarg
        j=j+1
        if(i.gt.j) then
        cpot(i,j)=cpot(j,i)
        goto 42
        endif
        c=0.d0
c nuclear coupling                                                                               
        do 50 k=1,ndim
        c=c+vnccd(r,ev(k))*evec(i,k)*evec(j,k)
 50     continue
        cpot(i,j)=c
 42     continue
 32     continue
        do 51 i=1,ndim
        cpot(i,i)=cpot(i,i)-vnd(r)
 51     continue
	return
	end
c****************************************************************
        function vnccd(r,rxt)
c
c woods-saxon potential for nuclear potential
c
c****************************************************************
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
        common /eigend/ depth 
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
                                                                                
        if(r.lt.50.d0) then
           rr=r-r0-rxt
        vnccd=-1.d0/(1.d0+exp(rr/rdif))
        else
        vnccd=0.d0
        endif
        return
        end
c****************************************************************
        function vnd(r)
c
c woods-saxon potential for nuclear potential
c
c****************************************************************
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
        common /eigend/ depth 
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
                                                                                
        if(r.lt.50.d0) then
           rr=r-r0
        vnd=-1.d0/(1.d0+exp(rr/rdif))
        else
        vnd=0.d0
        endif
        return
        end

c*****************************************************************
        function fct2(r)
c
c coulomb coupling form factor for target excitation
c (rotational e2 coupling)
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
        lambda=2
        rc=rt
        if(r.gt.rc) then
        fct2=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rt/r)**lambda
        else
        fct2=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
        endif
                                                                                
        fct2=fct2*(beta2t+2.d0*sqrt(5.d0/pi)*beta2t**2/7.d0)
                                                                                
        return
        end
c*****************************************************************
        function fct4(r)
c
c coulomb coupling form factor for target excitation
c (rotational e4 coupling)
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
                                                                              
        lambda=4
        rc=rt
        if(r.gt.rc) then
        fct4=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/r
     &                                         *(rt/r)**lambda
        else
        fct4=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
        endif
                                                                                
        fct4=fct4*(beta4t+9.d0*beta2t**2/7.d0/sqrt(pi))
        return
        end

c*****************************************************************
        function pcoul(r)
c
c spherical coulomb for target excitation
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
        rc=radc
        if(r.gt.rc) then
        pcoul=zp*zt/137.d0*hbar/r
        else
        pcoul=zp*zt/137.d0*hbar*(3.d0-(r/rc)*(r/rc))/(2.d0*rc)
        endif
        return
        end

c*****************************************************************
        function cpcoul(r)
c
c spherical coulomb for target excitation
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
        rc=radc
        if(r.gt.rc) then
        cpcoul=0.d0
        else
        cpcoul=zp*zt/137.d0*hbar*(3.d0-(r/rc)*(r/rc))/(2.d0*rc)
        endif
        return
        end

c*****************************************************************
        function spcoul(r)
c
c spherical coulomb for target excitation
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
        rc=radc
        if(r.gt.rc) then
        spcoul=zp*zt/137.d0*hbar
        else
        spcoul=0.d0
        endif
        return
        end

c******************************************************************
        subroutine cmats(r,cpot)
c
c coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      complex*16 cpot
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common /coulrad/ radc
      common/rott/beta2t,beta4t,ntarg
      common/dimen/nlevel
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/valpi/pi,hbar
        do 10 i=1,nlevelmax
        do 20 j=1,nlevelmax
        cpot(i,j)=0.d0
 20     continue
 10     continue
                                                                                
        ndim=nlevel
        i=0
        do 32 it=0,ntarg
        i=i+1
        j=0
        do 42 jt=0,ntarg
        j=j+1
        if(i.gt.j) then
        cpot(i,j)=cpot(j,i)
        goto 42
        endif
        c=0.d0
c coulomb coupling
        c=c+sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct2s(r)
     &  +sqrt((2.d0*2*it+1.d0)*9.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,4,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct4s(r)
        cpot(i,j)=c
 42     continue
 32     continue
	return
	end

c*****************************************************************
        function fct2s(r)
c
c coulomb coupling form factor for target excitation
c (rotational e2 coupling)
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
        lambda=2
        rc=rt
        if(r.gt.rc) then
        fct2s=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar
     &                                         *(rt/r)**lambda
        else
        fct2s=0.d0
        endif
                                                                                
        fct2s=fct2s*(beta2t+2.d0*sqrt(5.d0/pi)*beta2t**2/7.d0)
                                                                                
        return
        end
c*****************************************************************
        function fct4s(r)
c
c coulomb coupling form factor for target excitation
c (rotational e4 coupling)
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
                                                                              
        lambda=4
        rc=rt
        if(r.gt.rc) then
        fct4s=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar
     &                                         *(rt/r)**lambda
        else
        fct4s=0.d0
        endif
                                                                                
        fct4s=fct4s*(beta4t+9.d0*beta2t**2/7.d0/sqrt(pi))
        return
        end

c******************************************************************
        subroutine cmatc(r,cpot)
c
c coupling matrix
c
c*****************************************************************
      implicit real*8 (a-h,o-z)
      complex*16 cpot,vncc,vn
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
      common/hion/ap,zp,rp,at,zt,rt
      common /coulrad/ radc
      common/rott/beta2t,beta4t,ntarg
      common/dimen/nlevel
      common/eigen/ev(nlevelmax),evec(nlevelmax,nlevelmax)
      common/valpi/pi,hbar
        do 10 i=1,nlevelmax
        do 20 j=1,nlevelmax
        cpot(i,j)=0.d0
 20     continue
 10     continue
                                                                                
        ndim=nlevel
        i=0
        do 32 it=0,ntarg
        i=i+1
        j=0
        do 42 jt=0,ntarg
        j=j+1
        if(i.gt.j) then
        cpot(i,j)=cpot(j,i)
        goto 42
        endif
        c=0.d0
c nuclear coupling                                                                               
        do 50 k=1,ndim
        c=c+vncc(r,ev(k))*evec(i,k)*evec(j,k)
 50     continue
c coulomb coupling
        c=c+sqrt((2.d0*2*it+1.d0)*5.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,2,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct2c(r)
     &  +sqrt((2.d0*2*it+1.d0)*9.d0*(2.d0*2*jt+1.d0)
     &  /4.d0/pi)*cg(2*it,0,4,0,2*jt,0)**2/(2.d0*2*jt+1.d0)*fct4c(r)
        cpot(i,j)=c
 42     continue
 32     continue
        do 51 i=1,ndim
        cpot(i,i)=cpot(i,i)-vn(r)
 51     continue
	return
	end
c*****************************************************************
        function fct2c(r)
c
c coulomb coupling form factor for target excitation
c (rotational e2 coupling)
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
        lambda=2
        rc=rt
        if(r.gt.rc) then
        fct2c=0.d0
        else
        fct2c=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
        endif
                                                                                
        fct2c=fct2c*(beta2t+2.d0*sqrt(5.d0/pi)*beta2t**2/7.d0)
                                                                                
        return
        end
c*****************************************************************
        function fct4c(r)
c
c coulomb coupling form factor for target excitation
c (rotational e4 coupling)
c
c*****************************************************************
        implicit real*8(a-h,o-z)
        common/hion/ap,zp,rp,at,zt,rt
        common /coulrad/ radc
        common/rott/beta2t,beta4t,ntarg
        common/valpi/pi,hbar
                                                                              
        lambda=4
        rc=rt
        if(r.gt.rc) then
        fct4c=0.d0
        else
        fct4c=3.d0/(2.d0*lambda+1.d0)*zp*zt/137.d0*hbar/rc
     &                                         *(r/rt)**lambda
        endif
                                                                                
        fct4c=fct4c*(beta4t+9.d0*beta2t**2/7.d0/sqrt(pi))
        return
        end

c****************************************************************
      function cg(j1,m1,j2,m2,j3,m3)
c
c  clebsh-gordan coefficient <j1 m1 j2 m2 | j3 m3>
c****************************************************************
      implicit real*8(a-h,o-z)
	external fact

	if(m1+m2.ne.m3) then
	cg=0.d0
	return
	endif

	if(j3.lt.abs(j1-j2)) then
	cg=0.d0
	return
	endif

	if(j3.gt.j1+j2) then
	cg=0.d0
	return
	endif

      ka=j1+j2-j3
      kb=j3+j1-j2
      kc=j2+j3-j1
      kd=j1+j2+j3+1
      del=sqrt(fact(ka)*fact(kb)*fact(kc)/fact(kd))

      cg=0.d0

      do 10 n=0,max(j1+j2-j3,j1-m1,j2+m2)

      ka1=j1+j2-j3-n
      if(ka1.lt.0.d0) go to 10
      ka2=j3-j2+m1+n
      if(ka2.lt.0.d0) go to 10
      ka3=j3-j1-m2+n
      if(ka3.lt.0.d0) go to 10
      ka4=j1-m1-n
      if(ka4.lt.0.d0) go to 10
      ka5=j2+m2-n
      if(ka5.lt.0.d0) go to 10
      an=n*1.d0

      cg=cg+(-1.d0)**n
     &  /(fact(n)*fact(ka1)*fact(ka2)*fact(ka3)*fact(ka4)*fact(ka5))

 10   continue

 20   cg=cg*sqrt(fact(j1+m1)*fact(j1-m1))
      cg=cg*sqrt(fact(j2+m2)*fact(j2-m2))
      cg=cg*sqrt(fact(j3+m3)*fact(j3-m3))

      cg=cg*sqrt(2.d0*j3+1.d0)*del

      return
      end


                         function fact(n)
                         implicit real*8(a-h,o-z)
			 if(n.lt.0) then
			 fact=0.d0
			 return
			 endif
                         if (n.eq.0) then
                            fact=1.d0
                            go to 10
                            endif
                            fact=1.d0
                            do 20 i=1,n
                               fact=fact*i*1.d0
 20                            continue
 10                            return
                               end

c*************************************************************
	subroutine mdiag(a,n)
c
c the jacobi method to compute all eigenvalues and eigenvectors 
c of a real symmetric matrix a, which is of size n by n, stored 
c in a physical np by np array. on output, the elements of a 
c above the diagonal are destroyed. d returns the eigenvalues of 
c a in its first n elements. v is a matrix with the same logical 
c and physical dimensions as a whose columns contain, on output, 
c the normalized eigenvectors of a. nrot returns the number of 
c jacobi rotations which were required. 
c
c n.b.; the i-th component of the eigenvector for the k-th 
c       eigenvalue is given by v(i,k). 
c
c************************************************************
	implicit real*8(a-h,o-z)
	parameter(nmax=100)
        parameter (nlevelmax=25)
        dimension a(nlevelmax,nlevelmax)
	dimension b(nmax),z(nmax)
        common/eigen/d(nlevelmax),v(nlevelmax,nlevelmax)

	do 12 ip=1,n
	do 11 iq=1,n
	v(ip,iq)=0.d0
 11	continue
	v(ip,ip)=1.d0
 12	continue

	do 13 ip=1,n
	b(ip)=a(ip,ip)
	d(ip)=b(ip)
	z(ip)=0.d0
 13	continue

	nrot=0

	do 24 i=1,50
	sm=0.d0

	do 15 ip=1,n-1
	do 14 iq=ip+1,n
	sm=sm+abs(a(ip,iq))
 14	continue
 15     continue

	if(sm.eq.0.d0) return

	if(i.lt.4) then
	tresh=0.2d0*sm/n**2
	else
	tresh=0.d0
	endif

	do 22 ip=1,n-1
	do 21 iq=ip+1,n
	g=100.d0*abs(a(ip,iq))
	
	if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip)))
     &     .and.(abs(d(iq))+g.eq.abs(d(iq)))) then
	   a(ip,iq)=0.d0
	elseif(abs(a(ip,iq)).gt.tresh) then
	   h=d(iq)-d(ip)

	   if(abs(h)+g.eq.abs(h)) then
	     t=a(ip,iq)/h
	   else
	     theta=0.5d0*h/a(ip,iq)
	     t=1.d0/(abs(theta)+sqrt(1.d0+theta**2))
	     if(theta.lt.0.d0) t=-t
	   endif

	c=1.d0/sqrt(1.d0+t**2)
	s=t*c
	tau=s/(1.d0+c)
	h=t*a(ip,iq)
	z(ip)=z(ip)-h
	z(iq)=z(iq)+h
	d(ip)=d(ip)-h
	d(iq)=d(iq)+h
	a(ip,iq)=0.d0

	do 16 j=1,ip-1
	g=a(j,ip)
	h=a(j,iq)
	a(j,ip)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
 16	continue

	do 17 j=ip+1,iq-1
	g=a(ip,j)
	h=a(j,iq)
	a(ip,j)=g-s*(h+g*tau)
	a(j,iq)=h+s*(g-h*tau)
 17	continue

	do 18 j=iq+1,n
	g=a(ip,j)
	h=a(iq,j)
	a(ip,j)=g-s*(h+g*tau)
	a(iq,j)=h+s*(g-h*tau)
 18	continue

	do 19 j=1,n
	g=v(j,ip)
	h=v(j,iq)
	v(j,ip)=g-s*(h+g*tau)
	v(j,iq)=h+s*(g-h*tau)
 19	continue
	
	nrot=nrot+1
	endif

 21	continue
 22	continue

	do 23 ip=1,n
	b(ip)=b(ip)+z(ip)
	d(ip)=b(ip)
	z(ip)=0.d0

 23	continue
 24	continue

	pause '50 iterations should never happen'

	return
	end


       subroutine fcn(neq,r,y,ydp)
      include 'params'
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      dimension  y(neq), ydp(neq)
      common /valene/ ce(ndim)
        common /eigend/ depth
	common /fituv0/ pu0(ndim2,kmax),pv0(ndim2,kmax)
	common /fituv1/ pu1(ndim2,kmax),pv1(ndim2,kmax)
	common /fituv2/ pu2(ndim2,kmax),pv2(ndim2,kmax)
	common /partit/ rpart(kmax),npart
	if(r.eq.rpart(npart)) then
	rd=rpart(npart)-rpart(npart-1)
	rh=rd
	m2=npart-1
	go to 346
	endif
	do m1=1,npart-1
	if(r.ge.rpart(m1).and.r.lt.rpart(m1+1)) then
	m2=m1
	go to 345
	endif
	enddo
345	rd=r-rpart(m2)
	rh=rpart(m2+1)-rpart(m2)
346	rr=rd/rh
        rpol1=2.d0*rr-1.d0
	rpol2=1.d0+rr*(-6.d0+6.d0*rr)
	m=-ndim
	do j=1,ndim
	m=m+ndim
	do i=1,ndim
	k=m+i
	if(i.eq.j) then
	alam=ce(i)
        sum = -alam*y(k)
        suma = -alam*y(k)
	else
	sum=0.d0
	suma=0.d0
	endif
	l1=-ndim
      do l=1,ndim
	l1=l1+ndim
	valpua=pu0(l1+i,m2)+pu1(l1+i,m2)*rpol1/rh**2+
     +  pu2(l1+i,m2)*rpol2/rh**2
	valpva=pv0(l1+i,m2)+pv1(l1+i,m2)*rpol1/rh**2+
     +  pv2(l1+i,m2)*rpol2/rh**2
      suma=suma+(valpua+depth*valpva)*y(m+l)
      enddo
      ydp(k)=suma
      enddo
      enddo
      return
      end

******************************************************************
	subroutine pi1(n,np,rx1,lo,vc,v0,v1,v2,m)
c****************************************************************** 
c     this sbr prepares the integration on i1. it gives the
c     coefficients of the approximation of the potential.
c     input parameters:
c     n = number of equations
c     np= the physical dimension of the arrays rl, vc, v0, v1, v2.
c     rx1 = the radius at which the integration on i1 is performed.
c     output parameters:
c     rl = vector containing the l-values
c     vc = array containing the coefficient of 1/r
c     v0 = array containing the constant term of the potential
c     v1 = array containing the coefficient of r
c     v2 = array containing the coefficient of r**2
c     m  = dimension of the working arrays needed in persys.
c****************************************************************** 
	implicit double complex (a-h,o-q,s-z)
	implicit double precision (r)
	dimension cpfs(6),cpf(6)
	external spot,pot
	dimension vc(np,1),v0(np,1),v1(np,1),v2(np,1),lo(1)
	data maxit/10/
	rx2=rx1*rx1
	rx3=rx2*rx1
	rx4=rx2*rx2
	do 10 i=1,n
	do 10 j=1,n
	call sixval(spot,i,j,0.d0,rx1,cpfs)
	call fitpot(cpfs,rx1,vs0,vs1,vs2)
	call sixval(pot,i,j,0.d0,rx1,cpf)
	call fitpot(cpf,rx1,vn0,vn1,vn2)
	call copot(vs0,vs1,vs2,vn0,vn1,vn2,rx1,vts,vt0,vt1,vt2)
	vc(i,j)=vts/rx1
	v0(i,j)=vt0/rx2
	v1(i,j)=vt1/rx3
	v2(i,j)=vt2/rx4
10	continue
	do 20 i=1,n
	call ident(i,li,mi)
	lo(i)=li
20	continue
	call rangel(n,lo,maxit,m)
c	write(*,*) 'the dim m in rutine persys should be at least:',m
	return
	end
*********************************************************************
	subroutine copot(vs0,vs1,vs2,v0,v1,v2,r,vts,vt0,vt1,vt2)
	implicit double complex (a-h,o-q,s-z)
	implicit double precision (r)
	data two/(2.d0,0.d0)/,thr/(3.d0,0.d0)/,fou/(4.d0,0.d0)/,
     *  twe/(12.d0,0.d0)/
	vts=(vs0-two*(vs1-vs2))/r
	vt0=v0-two*(v1-v2)+fou*(vs1-thr*vs2)/r
	vt1=fou*(v1-thr*v2)+twe*vs2/r
	vt2=twe*v2
	return
	end

***********************************************************************
	subroutine fitpot(v,r,v0,v1,v2)
	implicit double complex (a-h,o-q,s-z)
	implicit double precision (r)
	dimension v(6)
	rr=r*r
      v0=rr*(19.d0*(v(1)+v(6))+75.d0*(v(2)+v(5))+5.d1*(v(3)+v(4)))/
     1288.d0
      v1=rr*(111.d0*(v(6)-v(1))+425.d0*(v(5)-v(2))
     1-15.d1*(v(4)-v(3)))/1344.d0
      v2=rr*25.d0*(5.d0*(v(1)+v(6))+6.d0*(v(2)+v(5))
     1-11.d0*(v(3)+v(4)))/1008.d0
      return
        end        
**************************************************************	
        subroutine sixval(func,i,j,r,rh,pf)
        implicit double precision (r)
	implicit double complex (a-h,o-q,s-z)
        dimension pf(6)
        external func
        do 1 k=1,6
        pf(k)=func(i,j,r+.2d0*rh*(k-1))
1	continue
	return
	end
**************************************************************	
	subroutine rangel(n,lo,maxit,m)
c	this subroutine produces a parameter needed in subroutine perth.
c	the value of this parameter depends on the n values of angular
c	momenta al and represents the maximum number of logarithmic terms
c	in the perturbative solution. the solution may be purely polynomial
c	or combination of the form: sum(c*x**p*log(x)**s). at each iteration
c	the maximum value of s either increases or remains unchanged.
	implicit double precision(a-h,o-z)
	dimension lo(*)
	m=1
	do i=1,n
	  iter=1
	  ism1=1
	    do j=1,n
	      do ip=1,4
	        lip=lo(i)+ip
	        ieta=lo(j)-lip
	        if(ieta.eq.0) then
	        ism1=2
	        endif
	      enddo
	    enddo
	  ism=ism1
	  ipmin=1
	  ipmax=4
	  ipm=ipmax-ipmin+1
	    do iter=2,maxit
	      ipm=ipm+3
	      ipmin=ipmin+1
	      ipmax=ipmax+4
	      ism1=ism
	        do j=1,n
	          do ip=ipmin,ipmax
	            lip=lo(i)+ip
	            ieta=lo(j)-lip
	            if(ieta.eq.0) then
	            ism1=ism+1
	            endif
	          enddo
	        enddo
	      ism=ism1
	    enddo
	  if(ism.gt.m) m=ism
	enddo
	return
	end

      subroutine soli1(rx,sy,syp,syk,sykp) 
      include 'params'
      parameter(neq=2*ndim2,neq2=2*neq)
      parameter(maxlog=14,maxit=15,mxd=3*maxit+1,lwrk=mxd*maxlog*ndim)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
	common/poti1/v1c(ndim,ndim),v10(ndim,ndim),v11(ndim,ndim),
     &  v12(ndim,ndim),lo(ndim)   
	common/worki1/aw(lwrk),yw(lwrk),adw(lwrk),ydw(lwrk),cw(ndim),
     &	reps,iter(ndim),ml
        common /valp/ ck(ndim)
	dimension sy(1),syp(1),syk(1),sykp(1)
        call persys(ndim,ndim,rx,lo,v1c,v10,v11,v12,reps,maxit,
     *ml,lwrk,aw,yw,adw,ydw,cw,sy(1),syp(1),syk(1),sykp(1),iter)
	do i=1,ndim2
	syk(i)=0.d0
	sykp(i)=0.d0
	enddo
c      write(*,*) (iter(id),id=1,ndim)
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine persys(n,np,rx,lo,vc,v0,v1,v2,reps,maxit,
     *maxlog,lwrk,wa,wy,wad,wyd,wu,sol,dsol,zol,dzol,iter)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	this subroutine computes the solution of a system of coupled
c	schrodinger equations in the vicinity of the origin. a pertur-
c	bative method is used. 
c	parameters:
c	-n     = number of equations.
c       -np    = physical dimension of the matrices.
c	-rx    = the point in which the solution is computed.
c       -ak    = value of k.
c	-lo    = array of dimension n containing the angular momenta.
c	-vc    = matrix of dimension n*n containing the coefficients 
c		     of r**(-1) (the coulombian term). 
c	-v0    = matrix of dimension n*n containing the coefficients 
c		     of r**0. 
c	-v1    = matrix of dimension n*n containing the coefficients 
c		     of r**1. 
c	-v2    = matrix of dimension n*n containing the coefficients 
c		     of r**2.
c	-reps  = the required accuracy.
c                its minimum value is 1.d-14.
c       -maxit = maximum number of iterations.
c                its value depends on the desired accuracy. 
c                a reasonable value is 15. 
c       -maxlog= number of logarithmic terms as given by the
c		     subroutine range.
c       -lwrk  = dimension of the working arrays wa and wy.
c	         lwrk = (3*maxit+1)*maxlog*n
c	-wa,wad= working arrays of dimension lwrk.
c	-wy,wyd= working arrays of dimension lwrk.
c       -wu    = working vector of dimension n.
c	-sol   = matrix of dimension n*n containing the solution.
c	-dsol  = matrix of dimension n*n containing the first derivative
c			 of the solution.
c	-zol   = matrix of dimension n*n containing the derivative of 
c                the solution with respect with k.
c	-dzol  = matrix of dimension n*n containing the first derivative
c		 of zol.
c	-iter  = integer vector of dimension n. on output iter(j) contains
c                the number of iterations required for the column j of
c                solution. if the desired accuracy is not reached after
c                maxit iterations, this parameter receives the value -maxit.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double complex(a-h,o-q,s-z)
      implicit double precision (r)
      include 'params'
      dimension lo(*),vc(np,*),v0(np,*),v1(np,*),v2(np,*),sol(np,*),dsol
     !(np,*),wu(*),wa(*),wy(*),zol(np,*),dzol(np,*),wad(*),wyd(*)
        common /valp/ ak(ndim)
      dimension iter(*)
c      dimension ddsol(10,10),ddzol(10,10),dif(10,10),diz(10,10)
      data czero/(0.d0,0.d0)/,cone/(1.d0,0.d0)/,rzero/0.d0/
c     rx2=rx*rx
      maxpol= 3*maxit+1
      lwrk1 = maxpol*maxlog*n
      if(lwrk1.gt.lwrk) then
      write(*,*) ' dimension of working vectors should be at least ',lwrk1
      stop
      endif
      clx=dlog(rx)
      do k=1,n
      ener=ak(k)*ak(k)
      da=-(ak(k)+ak(k))
        v0(k,k) = v0(k,k) - ener
      enddo
c	compute each column of solution
      do i=1,n
        kont=0
        ierc=0
        iacc=0
        it=0
        ism1=1
c	test for maximum number of logarithmic terms
	do j=1,n
	  do ip=1,4
	    lip=lo(i)+ip
	    ieta=lip*(lip+1)-lo(j)*(lo(j)+1)
	    if(ieta.eq.0) then
	      ism1=2
	    endif
	  enddo
	enddo
	ism=ism1
c     starting values of sol, dsol, zol and dzol
	jm=-maxlog
	do j=1,n
	  zol(j,i)=czero
	  dzol(j,i)=czero
c	  ddzol(j,i)=czero
	  jm=jm+maxlog
	  if(j.eq.i) then
	    lip=lo(i)+1
	    sol(i,i)=rx**lip
	    dsol(i,i)=lip*rx**(lip-1)
c	    ddsol(i,i)=lip*(lip-1)*rx**(lip-2)
	  else
	    sol(j,i)=czero
	    dsol(j,i)=czero
c	    ddsol(j,i)=czero
	  endif
	  do is=1,ism
	    iy0=(jm+is-1)*maxpol
	    do ip=1,4
	      iy=iy0+ip
	      wy(iy)=czero
	      wyd(iy)=czero
	    enddo
          enddo
	enddo
304	format(3(2g11.4,3x))	
305	format(3(2g11.4,3x))	
	ipm=4
c    first iteration
	it=it+1
	jm=-maxlog
	do j=1,n
	  jm=jm+maxlog
	  ia=jm*maxpol
	  wa(ia+1)=vc(j,i)
	  wa(ia+2)=v0(j,i)
	  wa(ia+3)=v1(j,i)
	  wa(ia+4)=v2(j,i)
	  wad(ia+1)=czero
	  wad(ia+2)=czero
	  if(j.eq.i) wad(ia+2) = da
          wad(ia+3)=czero
	  wad(ia+4)=czero
	  do ip=1,4
	    lip=lo(i)+ip
	    ieta=lip*(lip+1)-lo(j)*(lo(j)+1)
	    ip1=jm*maxpol+ip
	    ip2=(jm+1)*maxpol+ip
	    if(ieta.eq.0) then
	      ih=2*lip+1
	      wy(ip2)=wa(ip1)/ih
	      wyd(ip2)=wad(ip1)/ih
	    else
              wy(ip1)=wa(ip1)/ieta
              wyd(ip1)=wad(ip1)/ieta
            endif
c    the values of sol,dsol,zol,dzol after the first iteration
            s1=czero
	    s2=czero
	    s3=czero
	    z1=czero
	    z2=czero
	    z3=czero
	    do is=1,ism
	      iy=(jm+is-1)*maxpol+ip
	      if(is.eq.1) then
	        clx1=cone
	        clx2=czero
	      else
	        clx1=clx
	        clx2=cone
	      endif
	      s1=s1+wy(iy)*clx1
	      c2=(lip+1)*clx1+(is-1)*clx2
	      s2=s2+wy(iy)*c2
	      c3=lip*(lip+1)*clx1+(is-1)*(2*lip+1)*clx2
	      s3=s3+wy(iy)*c3
	      z1=z1+wyd(iy)*clx1
	      z2=z2+wyd(iy)*c2
	      z3=z3+wyd(iy)*c3
	    enddo
	      rxp=rx**lip
	      rxp1=rxp*rx
	      rx1p=rxp/rx
              sol(j,i)=sol(j,i)+s1*rxp1
	      dsol(j,i)=dsol(j,i)+s2*rxp
c	      ddsol(j,i)=ddsol(j,i)+s3*rx1p
              zol(j,i)=zol(j,i)+z1*rxp1
	      dzol(j,i)=dzol(j,i)+z2*rxp
c	      ddzol(j,i)=ddzol(j,i)+z3*rx1p
	  enddo
        enddo
c	  do j=1,n
c       	     sum=czero
c	     zum=czero
c	  do l=1,n
c             sum=sum+(vc(j,l)/rx+v0(j,l)+v1(j,l)*rx+v2(j,l)*rx2)
c     &      *sol(l,i)
c	     del=czero
c	     if(j.eq.l) del=cone
c             zum=zum+(vc(j,l)/rx+v0(j,l)+v1(j,l)*rx+v2(j,l)*rx2)*
c     &       zol(l,i)+del*da*sol(l,i)
c   	 enddo
c             dif(j,i)=ddsol(j,i)-lo(j)*(lo(j)+1)*sol(j,i)/rx2-sum
c             diz(j,i)=ddzol(j,i)-lo(j)*(lo(j)+1)*zol(j,i)/rx2-zum
c   	 enddo
	ipmin=1
	ipmax=4
	ipm=ipmax-ipmin+1
	do while(iacc.eq.0.and.ierc.eq.0)
	  it=it+1
c	test for number of iterations
	  if(it.le.maxit) then
c	go on with iteration
	    jm=-maxlog
	    do j=1,n
	      jm=jm+maxlog
	      do is=1,ism
	        ia=(jm+is-1)*maxpol+1
	        wa(ia)=czero
	        wad(ia)=czero
	        lm=-maxlog
	        do l=1,n
	          lm=lm+maxlog
	          iy=(lm+is-1)*maxpol+1
	          wa(ia)=wa(ia)+vc(j,l)*wy(iy)
	          wad(ia)=wad(ia)+vc(j,l)*wyd(iy)
	        enddo
	        ia=ia+1
	        wa(ia)=czero
	        wad(ia)=czero
	        lm=-maxlog
                do l=1,n
	          lm=lm+maxlog
	          iy=(lm+is-1)*maxpol
	          iy1=iy+1
	          iy2=iy+2
	          wa(ia)=wa(ia)+vc(j,l)*wy(iy2)+v0(j,l)*wy(iy1)
	          del=czero
	          if(j.eq.l) del=cone
	          wad(ia)=wad(ia)+vc(j,l)*wyd(iy2)+v0(j,l)*wyd(iy1)+         
     !		  del*da*wy(iy1)
                enddo
	        ia=ia+1
                wa(ia)=czero
                wad(ia)=czero
                lm=-maxlog
                do l=1,n
                  lm=lm+maxlog
                  iy=(lm+is-1)*maxpol
                  iy1=iy+1
                  iy2=iy+2
                  iy3=iy+3
                  wa(ia)=wa(ia)+vc(j,l)*wy(iy3)+v0(j,l)*wy(iy2)+
     !		  v1(j,l)*wy(iy1)
                  del=czero
                  if(j.eq.l) del=cone
                  wad(ia)=wad(ia)+vc(j,l)*wyd(iy3)+
     !	          v0(j,l)*wyd(iy2)+v1(j,l)*wyd(iy1)+del*da*wy(iy2)
	        enddo
	        ia0=(jm+is-1)*maxpol
	        do ip=4,ipm
	          ia=ia0+ip
	          wa(ia)=czero
	          wad(ia)=czero
	          lm=-maxlog
	          do l=1,n
	            lm=lm+maxlog
	            iy0=(lm+is-1)*maxpol+ip
	            iy1=iy0-1
	            iy2=iy0-2
	            iy3=iy0-3
	            wa(ia)=wa(ia)+vc(j,l)*wy(iy0)+
     !              v0(j,l)*wy(iy1)+v1(j,l)*wy(iy2)+v2(j,l)*wy(iy3)       
                    del=czero
	            if(j.eq.l) del=cone
	            wad(ia)=wad(ia)+vc(j,l)*wyd(iy0)+v0(j,l)*wyd(iy1)+
     !  	    v1(j,l)*wyd(iy2)+v2(j,l)*wyd(iy3)+del*da*wy(iy1)
	          enddo
	        enddo
	        ip=ipm+1
	        ia=(jm+is-1)*maxpol+ip
	        wa(ia)=czero
	        wad(ia)=czero
	        lm=-maxlog
	        do l=1,n
	          lm=lm+maxlog
	          iy0=(lm+is-1)*maxpol+ipm
	          iy1=iy0-1
	          iy2=iy0-2
	          wa(ia)=wa(ia)+v0(j,l)*wy(iy0)+v1(j,l)*wy(iy1)+
     !		  v2(j,l)*wy(iy2)
                  del=czero
                  if(j.eq.l) del=cone
                  wad(ia)=wad(ia)+v0(j,l)*wyd(iy0)+v1(j,l)*wyd(iy1)+
     !	          v2(j,l)*wyd(iy2)+del*da*wy(iy0)
                enddo
	        ip=ipm+2
	        ia=(jm+is-1)*maxpol+ip
	        wa(ia)=czero
	        wad(ia)=czero
	        lm=-maxlog
	        do l=1,n
	          lm=lm+maxlog
	          iy0=(lm+is-1)*maxpol+ipm
	          iy1=iy0-1
	          wa(ia)=wa(ia)+v1(j,l)*wy(iy0)+v2(j,l)*wy(iy1)
	          wad(ia)=wad(ia)+v1(j,l)*wyd(iy0)+v2(j,l)*wyd(iy1)
                enddo
                ip=ipm+3
                ia=(jm+is-1)*maxpol+ip
                wa(ia)=czero
                wad(ia)=czero
                lm=-maxlog
                do l=1,n
                  lm=lm+maxlog
                  iy0=(lm+is-1)*maxpol+ipm
                  wa(ia)=wa(ia)+v2(j,l)*wy(iy0)
                  wad(ia)=wad(ia)+v2(j,l)*wyd(iy0)
                enddo
              enddo
	    enddo
	    ipm=ipm+3
	    ipmin=ipmin+1
	    ipmax=ipmax+4
	    ism1=ism
	    ism2=ism
	    do j=1,n
	      do ip=ipmin,ipmax
	        lip=lo(i)+ip
	        ieta=lip*(lip+1)-lo(j)*(lo(j)+1)
	        if(ieta.eq.0) then
	          ism1=ism+1
	        endif
	      enddo
	    enddo
	    ism=ism1
	    jm=-maxlog
	    do j=1,n
	      jm=jm+maxlog
	      do is=1,ism
	        iy0=(jm+is-1)*maxpol
	        do ip=1,ipm
	          iy=iy0+ip
	          wy(iy)=czero
	          wyd(iy)=czero
	        enddo
	      enddo
	    enddo
	    jm=-maxlog
	    do j=1,n
	      jm=jm+maxlog
	      wu(j)=czero
	      do ip=ipmin,ipmax
	        ipp=ip-ipmin+1
	        lip=lo(i)+ip
	        ieta=lip*(lip+1)-lo(j)*(lo(j)+1)
	        if(ieta.eq.0) then
	          ih=2*lip+1
	          iy1=(jm+ism2)*maxpol+ipp
	          iy0=(jm+ism2-1)*maxpol+ipp
	          wy(iy1)=wa(iy0)/(ih*ism2)
	          wyd(iy1)=wad(iy0)/(ih*ism2)
	          if(ism2.gt.1) then
	            k=ism2
	            do k1=ism2-1,1,-1
	              iy0=(jm+k-1)*maxpol+ipp
	              ia1=(jm+k1-1)*maxpol+ipp
	              iy1=(jm+k)*maxpol+ipp
	              wy(iy0)=(wa(ia1)-k*k1*wy(iy1))/(ih*k1)
	              wyd(iy0)=(wad(ia1)-k*k1*wyd(iy1))/(ih*k1)
	              k=k1
	            enddo
	          endif
	        else
	          ih1=lip*(lip+1)-lo(j)*(lo(j)+1)
	          iy=(jm+ism2-1)*maxpol+ipp
	          wy(iy)=wa(iy)/ih1
	          wyd(iy)=wad(iy)/ih1
	          if(ism2.ne.1) then
	            k=ism2-1
	            ih=2*lip+1
	            iy=(jm+ism2-1)*maxpol+ipp
	            ia=(jm+k-1)*maxpol+ipp
	            wy(ia)=(-ih*k*wy(iy)+wa(ia))/ih1
	            wyd(ia)=(-ih*k*wyd(iy)+wad(ia))/ih1
	          endif
	          if(ism2.gt.2) then
	            do k=ism2-2,1,-1
	              iy0=(jm+k-1)*maxpol+ipp
	              iy1=(jm+k)*maxpol+ipp
	              iy2=(jm+k+1)*maxpol+ipp
	              wy(iy0)=(wa(iy0)-k*(ih*wy(iy1)+(k+1)*wy(iy2)))/ih1      
	              wyd(iy0)=(wad(iy0)-k*(ih*wyd(iy1)+
     !	              (k+1)*wyd(iy2)))/ih1
                    enddo
                  endif
	        endif
c     the values of sol,dsol,zol,dzol after a new iteration
	        s1=czero
	        s2=czero
c	        s3=czero
	        z1=czero
	        z2=czero
c	        z3=czero
	        do is=1,ism
	          iy=(jm+is-1)*maxpol+ipp
	          if(is.eq.1) then
	            clx1=cone
	            clx2=czero
c	            clx3=czero
		  elseif(is.eq.2) then
		    clx1=clx
		    clx2=cone
c		    clx3=czero
		  elseif(is.eq.3) then
		    clx1=clx*clx
		    clx2=clx
c		    clx3=cone
	          else
	            clx1=clx**(is-1)
	            clx2=clx**(is-2)
c	            clx3=clx**(is-3)
	          endif
	          s1=s1+wy(iy)*clx1
	          c2=(lip+1)*clx1+(is-1)*clx2
	          s2=s2+wy(iy)*c2
c	          c3=lip*(lip+1)*clx1+(is-1)*(2*lip+1)*clx2+
c     !		  (is-1)*(is-2)*clx3
c                 s3=s3+wy(iy)*c3
	          z1=z1+wyd(iy)*clx1
	          z2=z2+wyd(iy)*c2
c	          z3=z3+wyd(iy)*c3
	        enddo
	      rxp=rx**lip
	      rxp1=rxp*rx
	      rx1p=rxp/rx
              sol(j,i)=sol(j,i)+s1*rxp1
	      dsol(j,i)=dsol(j,i)+s2*rxp
c	      ddsol(j,i)=ddsol(j,i)+s3*rx1p
              zol(j,i)=zol(j,i)+z1*rxp1
	      dzol(j,i)=dzol(j,i)+z2*rxp
c	      ddzol(j,i)=ddzol(j,i)+z3*rx1p
	        wu(j)=wu(j)+s1*rxp1
	      enddo
	    enddo
c	  do j=1,n
c	     sum=czero
c	     zum=czero
c	  do l=1,n
c             sum=sum+(vc(j,l)/rx+v0(j,l)+v1(j,l)*rx+v2(j,l)*rx2)
c     &      *sol(l,i)
c	     del=czero
c	     if(j.eq.l) del=cone
c             zum=zum+(vc(j,l)/rx+v0(j,l)+v1(j,l)*rx+v2(j,l)*rx2)*
c     &       zol(l,i)+del*da*sol(l,i)
c  	  enddo
c             dif(j,i)=ddsol(j,i)-lo(j)*(lo(j)+1)*sol(j,i)/rx2-sum
c             diz(j,i)=ddzol(j,i)-lo(j)*(lo(j)+1)*zol(j,i)/rx2-zum
c   	  enddo
c    test of accuracy
	    rolmax=rzero
	    do j=1,n
	      rabsol=abs(sol(j,i))
	      if(rabsol.gt.rolmax) rolmax=rabsol
	    enddo
	    reler=rzero
	    do j=1,n
	      rabsol=abs(wu(j))/rolmax
	      if(rabsol.gt.reler) reler=rabsol
	    enddo
            if(reler.le.reps) kont=kont+1
            if(reler.le.reps.and.kont.gt.1) then
              iter(i)=it
              iacc=1
            endif
	  else
	    ierc=1
	    iter(i)=-maxit
	  endif
	enddo
      enddo
c   the matrix  v0 is restored.
      do k=1,n
      ener=ak(k)*ak(k)
        v0(k,k) = v0(k,k) + ener
      enddo
      return
      end

c.................................................
      subroutine corr(rrv0,riv0,delta)
c     propagates solution and computes correction
c     for eigenvalue
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      include 'params'
      parameter(neq=2*ndim2,neq2=2*neq)
      parameter(maxlog=14,maxit=15,mxd=3*maxit+1,lwrk=mxd*maxlog*ndim)
      common /valp/ ck(ndim)
	common /eigend/ depth
      common /solsd/ sy(neq),syp(neq),dy(neq),dyp(neq),wght(ndd)
      common /solut/ sol(kmax,ndim2),dsol(kmax,ndim2)
      	common /pnuclear/ pnuc(kmax,ndim2)
      common /totsol/ tsol(kmax,ndim),tdsol(kmax,ndim)
      dimension a(ndd,ndd),b(ndd,ndd),c(ndd),d(ndd),
     *yy(ndd,ndd),zz(ndd,ndd),aa(ndd1,ndd1),indx(ndd),index(ndd),dd(ndd)
	dimension bl(ndim),cr(ndim)
      dimension cw1(ndim2),cw2(ndim2),cw3(ndim2),cw4(ndim2),back(ndim2)
      dimension indxr(ndim)
      external  fcn
      common /tol/ rtol
      common /pivot/ ipiv,jpiv
      common /radii/ rinit,rmatch,rfin,rdx
      common /matchn/ kmatc
      common /radial/ raz(0:kmax),kfin
	common /partit/ rpart(kmax),npart
	common/poti1/v1c(ndim,ndim),v10(ndim,ndim),v11(ndim,ndim),
     &  v12(ndim,ndim),lo(ndim)   
	common/worki1/aw(lwrk),yw(lwrk),adw(lwrk),ydw(lwrk),cw(ndim),
     &	reps,iter(ndim),ml
      common /dprec/ rprec
	data ci/(0.d0,1.d0)/,cunit/(1.d0,0.d0)/,czero/(0.d0,0.d0)/

	depth=dcmplx(rrv0,riv0)
	do i=1,ndd
	do j=1,ndd
	yy(i,j)=0.d0
	zz(i,j)=0.d0
	enddo
	enddo
*
*     initial conditions
*
	 call pi1(ndim,ndim,rpart(1),lo,v1c,v10,v11,v12,ml)
       call soli1(rpart(1),sy,syp,sy(ndim2+1),syp(ndim2+1)) 
       do i=1,ndim2
       sol(1,i)=sy(i)
       dsol(1,i)=syp(i)
	 dy(i)=sy(i)
	 dyp(i)=syp(i)
       enddo
	 if=1
c	 call dcunitm(back,ndim,ndim,ndim)
c      left-hand solution:
	 do kk=1,kmatc-1
	 r=rpart(kk)
	 rh=rpart(kk+1)-rpart(kk)
c	 if(kk.ne.1) then
c	 call renorm(ndim,sy,syp,sy(ndim2+1),syp(ndim2+1),cw1,cw2,cw3,
c     +   cw4,back,indxr)
c	 endif
       call prop(r,rh,kk,if,sy,syp)
       do i=1,ndim2
       sol(kk+1,i)=sy(i)
       dsol(kk+1,i)=syp(i)
       enddo
	 enddo
c         call dcgmpr(sy,back,cw3,ndim,ndim,ndim)
c         call dcgmpr(syp,back,cw4,ndim,ndim,ndim)
c	 do i=1,ndim2
c	 sy(i)=cw3(i)
c	 syp(i)=cw4(i)
c	 enddo
	 k=0
	 do j=1,ndim
	 do i=1,ndim
	 k=k+1
	 yy(i,j)=sy(k)
	 yy(i+ndim,j)=syp(k)
	 enddo
	 enddo

c      initial conditions for right-hand solution:
       call soli3(rpart(npart),dy,dyp,dy(ndim2+1),dyp(ndim2+1))
       do i=1,ndim2
       sol(npart,i)=dy(i)
       dsol(npart,i)=dyp(i)
 	 sy(i)=dy(i)
	 syp(i)=dyp(i)
       enddo
       if=2
	 do kk=npart-1,kmatc,-1
	 r=rpart(kk+1)
       rh = r - rpart(kk)
       call prop(r,rh,kk,if,dy,dyp)
       do i=1,ndim2
       sol(kk,i)=dy(i)
       dsol(kk,i)=dyp(i)
       enddo
	 enddo
	 k=0
	 do j=1,ndim
	 do i=1,ndim
	 k=k+1
	 yy(i,j+ndim)=-dy(k)
	 yy(i+ndim,j+ndim)=-dyp(k)
	 enddo
	 enddo

	do j=1,ndd
	do i=1,ndd
	a(i,j)=yy(i,j)
	enddo
	enddo


	ipiv=ndim+1
	jpiv=ndim+1

	ki=0
	do i=1,ndd
	if(i.ne.ipiv) then
	ki=ki+1
	kj=0
	do j=1,ndd
	if(j.ne.jpiv) then
	kj=kj+1
	aa(ki,kj)=a(i,j)
	endif
	enddo
	endif
	enddo
	ki=0
	do i=1,ndd
	if(i.ne.ipiv) then
	ki=ki+1
	d(ki)=-a(i,jpiv)
	endif
	enddo
      call ludcmp(aa,ndd1,indx,det,ierror)
      call lubksb(aa,ndd1,indx,d)
	do jj=1,ndim
	bl(jj)=d(jj)
	enddo
	cr(1)=cunit
	if(ndim.gt.1) then
	kj=0
	do jj=2,ndim
	kj=kj+1
	cr(jj)=d(ndim+kj)
	enddo
	endif

	phi=czero
	ppl=czero
	ppr=czero
	do j=1,ndim
	phi=phi+a(1,j)*bl(j)
	ppl=ppl+a(ipiv,j)*bl(j)
	ppr=ppr-a(ipiv,ndim+j)*cr(j)
	enddo
	do j=1,ndim
	wght(j)=bl(j)
	wght(ndim+j)=cr(j)
	enddo

	do k=1,kmatc-1
        kk=-ndim
        do j=1,ndim
        kk=kk+ndim
        do i=1,ndim
        sol(k,kk+i)=sol(k,kk+i)*wght(j)
        dsol(k,kk+i)=dsol(k,kk+i)*wght(j)
        enddo
        enddo
	enddo

         do k=npart,kmatc,-1
         kk=-ndim
         do j=1,ndim
         kk=kk+ndim
         do i=1,ndim
         sol(k,kk+i)=sol(k,kk+i)*wght(j+ndim)
         dsol(k,kk+i)=dsol(k,kk+i)*wght(j+ndim)
         enddo
         enddo
	 enddo
         do kk=1,npart
         do i=1,ndim
         sum=0.d0
         dsum=0.d0
         do j=1,ndim
         k=(j-1)*ndim+i
         sum=sum+sol(kk,k)
         dsum=dsum+dsol(kk,k)
         enddo
         tsol(kk,i)=sum
         tdsol(kk,i)=dsum
         enddo
         enddo

	sum=czero
	k=-ndim
	do j=1,ndim
        k=k+ndim
        do i=1,ndim
	quads=0.d0
	do kk=1,npart-1
	quadl=tsol(kk,i)*pnuc(kk,k+i)*tsol(kk,j)
	quadr=tsol(kk+1,i)*pnuc(kk+1,k+i)*tsol(kk+1,j)
	quadp=0.5d0*(quadl+quadr)*(rpart(kk+1)-rpart(kk))
	quads=quads+quadp
        enddo
	sum=sum+quads
	enddo
	enddo
c	write(*,*) 'ppl =',ppl
c	write(*,*) 'ppr =',ppr
	delta=phi*(ppr-ppl)/sum
        return
        end

      subroutine soli3(rx,dy,dyp,dyk,dykp)
c     right hand boundary condition
      include 'params'
      parameter (neq=2*ndim2)
      parameter(maxlog=14,maxit=15,mxd=3*maxit+1,lwrk=mxd*maxlog*ndim)
	implicit double complex (a-h,o-q,s-z)
	implicit double precision(r)
        common /valp/ ck(ndim)
        dimension dy(1),dyp(1),dyk(1),dykp(1)
c	dimension rfc(3),rfcp(3),rgc(3),rgcp(3)
	data czero/(0.d0,0.d0)/
	k=0
	do i=1,ndim
	do j=1,ndim
	k=k+1
	dy(k)=czero
	dyp(k)=czero
	dyk(k)=czero
	dykp(k)=czero
	enddo
	enddo
	ii=-ndim
	do 20 i=1,ndim
	call ident(i,l,m)
	coul=spot(i,i,rx)
	ii=ii+ndim
	call assolg(rx,ck(i),l,coul,dy(ii+i),dyp(ii+i),dyk(ii+i),
     &	dykp(ii+i))
c	alpha=0.5d0*coul
c	rck=ck(i) 
c	reta=alpha/rck
c	rho=rck*rx
c        call rcwfn(rho,reta,l,l,rfc,rfcp,rgc,rgcp,1.d-14,999.d0)
c	rgcp(l+1)=rgcp(l+1)*rck
c	write(45,*) 'i =',i
c	write(45,*) dreal(dy(ii+i)),dreal(dyp(ii+i)),
c     +  dreal(dyp(ii+i)/dy(ii+i))
c	write(45,*) rgc(l+1),rgcp(l+1),rgcp(l+1)/rgc(l+1)
20	continue
	k=0
	do i=1,ndim
	do j=1,ndim
	k=k+1
	dyk(k)=czero
	dykp(k)=czero
	enddo
	enddo
	return
	end

c******************************************************************
	subroutine assolg(rx,pe,l,coul,u,up,uk,ukp)
c*****************************************************************
c     this sbr computes the solution g_l of the coulomb equation *
c                                                                *
c     input parameters:                                          *
c     x = radius                                                 *
c     pe = k-value                                               *
c     l = angular momentum                                       *
c     coul = 2*k*eta, where eta is the sommerfeld parameter.     *
c     output parameters:                                         *
c     u = solution                                               *
c     up = d(u)/dr                                               *
c     uk = d(u)/dk                                               *
c     ukp= d(d(u)/dk)/dr                                         *
c*****************************************************************
      implicit double precision(r)
      implicit double complex(a-h,o-q,s-z)
        dimension cgc(0:2),cgcp(0:2),cfc(0:2),cfcp(0:2),sig(0:2)
	data ci/(0.d0,1.d0)/
            zlp=l
            nlp=1
            ifail=1
	    mode1=1
	    ipr=1
	    kfn=0
	alpha=0.5d0*coul
 	eta=alpha/pe
	xrho=pe*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	u=cgc(0)
	up=cgcp(0)*pe

	de=1.d-5*pe
	    mode1=1
	ckk=pe+2.d0*de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	u2=cgc(0)
	up2=cgcp(0)*ckk

	ckk=pe+de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	u1=cgc(0)
	up1=cgcp(0)*ckk

	ckk=pe-de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	v1=cgc(0)
	vp1=cgcp(0)*ckk

	ckk=pe-2.d0*de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	v2=cgc(0)
	vp2=cgcp(0)*ckk
	uk=(v2-8.d0*v1+8.d0*u1-u2)/(12.d0*de)
	ukp=(vp2-8.d0*vp1+8.d0*up1-up2)/(12.d0*de)
	return
	end

      subroutine rcwfn(rho,eta,minl,maxl,fc,fcp,gc,gcp,accur,step)
      double precision acc,accur,ai,ar,bi,br,d,del,di,dk,dp,dq,dr,eta,
     *etah,etar,eta2,f,fc,fcp,fp,g,gcp,gp,h,h2,h2ll,k,k1,k2,k3,k4,m1,
     *m2,m3,m4,p,pace,pl,pmx,q,r,rh,rho,rho2,rh2,r3,s,step,t,tf,tfp,
     *turn,w,wi,xll1,gc
      dimension fc(1),fcp(1),gc(1),gcp(1)
      pace=step
      acc=accur
      if(pace.lt.100.0d0) pace=100.0d0
      if(acc.lt.1.0d-15.or.acc.gt.1.0d-6) acc=1.0d-6
      r=rho
      ktr=1
      lmax=maxl
      lmin1=minl+1
      xll1=minl*lmin1
      eta2=eta*eta
      turn=eta+dsqrt(eta2+xll1)
      if(r.lt.turn.and.dabs(eta).ge.1.0d-6) ktr=-1
      ktrp=ktr
      go to 2
    1 r=turn
      tf=f
      tfp=fp
      lmax=minl
      ktrp=1
    2 etar=eta*r
      rho2=r*r
      pl=lmax+1
      pmx=pl+0.5d0
      fp=eta/pl+pl/r
      dk=etar*2.0d0
      del=0.0d0
      d=0.0d0
      f=1.0d0
      k=(pl*pl-pl+etar)*(2.0d0*pl-1.0d0)
      if(pl*pl+pl+etar.ne.0.0d0) go to 3
      r=r+1.0d-6
      go to 2
    3 h=(pl*pl+eta2)*(1.0d0-pl*pl)*rho2
      k=k+dk+pl*pl*6.0d0
      d=1.0d0/(d*h+k)
      del=del*(d*k-1.0d0)
      if(pl.lt.pmx) del=-r*(pl*pl+eta2)*(pl+1.0d0)*d/pl
      pl=pl+1.0d0
      fp=fp+del
      if(d.lt.0.0d0) f=-f
      if(pl.gt.20000.0d0) go to 11
      if(dabs(del/fp).ge.acc) go to 3
      fp=f*fp
      if(lmax.eq.minl) go to 5
      fc(lmax+1)=f
      fcp(lmax+1)=fp
      l=lmax
      do 4 lp=lmin1,lmax
      pl=l
      gc(l+1)=eta/pl+pl/r
      gcp(l+1)=dsqrt(eta2+pl*pl)/pl
      fc(l)=(gc(l+1)*fc(l+1)+fcp(l+1))/gcp(l+1)
      fcp(l)=gc(l+1)*fc(l)-gcp(l+1)*fc(l+1)
    4 l=l-1
      f=fc(lmin1)
      fp=fcp(lmin1)
    5 if(ktrp.eq.-1) go to 1
      p=0.0d0
      q=r-eta
      pl=0.0d0
      ar=-(eta2+xll1)
      ai=eta
      br=2.0d0*q
      bi=2.0d0
      wi=2.0d0*eta
      dr=br/(br*br+bi*bi)
      di=-bi/(br*br+bi*bi)
      dp=-(ar*di+ai*dr)
      dq=(ar*dr-ai*di)
    6 p=p+dp
      q=q+dq
      pl=pl+2.0d0
      ar=ar+pl
      ai=ai+wi
      bi=bi+2.0d0
      d=ar*dr-ai*di+br
      di=ai*dr+ar*di+bi
      t=1.0d0/(d*d+di*di)
      dr=t*d
      di=-t*di
      h=br*dr-bi*di-1.0d0
      k=bi*dr+br*di
      t=dp*h-dq*k
      dq=dp*k+dq*h
      dp=t
      if(pl.gt.46000.0d0) go to 11
      if(dabs(dp)+dabs(dq).ge.(dabs(p)+dabs(q))*acc) go to 6
      p=p/r
      q=q/r
      g=(fp-p*f)/q
      gp=p*g-q*f
      w=1.0d0/dsqrt(fp*g-f*gp)
      g=w*g
      gp=w*gp
      if(ktr.eq.1) go to 8
      f=tf
      fp=tfp
      lmax=maxl
      if(rho.lt.0.2d0*turn) pace=999.0d0
      r3=1.0d0/3.0d0
      h=(rho-turn)/(pace+1.0d0)
      h2=0.5d0*h
      i2=pace+0.001d0
      etah=eta*h
      h2ll=h2*xll1
      s=(etah+h2ll/r)/r-h2
    7 rh2=r+h2
      t=(etah+h2ll/rh2)/rh2-h2
      k1=h2*gp
      k2=h2*(gp+m1)
      m1=s*g
      m2=t*(g+k1)
      k3=h*(gp+m2)
      m3=t*(g+k2)
      m3=m3+m3
      k4=h2*(gp+m3)
      rh=r+h
      s=(etah+h2ll/rh)/rh-h2
      m4=s*(g+k3)
      g=g+(k1+k2+k2+k3+k4)*r3
      gp=gp+(m1+m2+m2+m3+m4)*r3
      r=rh
      i2=i2-1
      if(dabs(gp).gt.1.0d70) go to 11
      if(i2.ge.0) go to 7
      w=1.0d0/(fp*g-f*gp)
    8 gc(lmin1)=g
      gcp(lmin1)=gp
      if(lmax.eq.minl) go to 10
      do 9 l=lmin1,lmax
      t=gc(l+1)
      gc(l+1)=(gc(l)*gc(l+1)-gcp(l))/gcp(l+1)
      gcp(l+1)=gc(l)*gcp(l+1)-gc(l+1)*t
      fc(l+1)=w*fc(l+1)
    9 fcp(l+1)=w*fcp(l+1)
      fc(lmin1)=fc(lmin1)*w
      fcp(lmin1)=fcp(lmin1)*w
      return
   10 fc(lmin1)=w*f
      fcp(lmin1)=w*fp
      return
   11 w=0.0d0
      g=0.0d0
      gp=0.0d0
      go to 8
      end 

      subroutine solfi3(rx,dy,dyp,dyk,dykp)
      include 'params'
      parameter (neq=2*ndim2)
      parameter(maxlog=14,maxit=15,mxd=3*maxit+1,lwrk=mxd*maxlog*ndim)
	implicit double complex (a-h,o-q,s-z)
	implicit double precision(r)
        common /valp/ ck(ndim)
        dimension dy(1),dyp(1),dyk(1),dykp(1)
	data czero/(0.d0,0.d0)/
	k=0
	do i=1,ndim
	do j=1,ndim
	k=k+1
	dy(k)=czero
	dyp(k)=czero
	dyk(k)=czero
	dykp(k)=czero
	enddo
	enddo
	ii=-ndim
	do 20 i=1,ndim
	call ident(i,l,m)
	coul=spot(i,i,rx)
	ii=ii+ndim
	call assolf(rx,ck(i),l,coul,dy(ii+i),dyp(ii+i),dyk(ii+i),
     &	dykp(ii+i))
20	continue
	k=0
	do i=1,ndim
	do j=1,ndim
	k=k+1
	dyk(k)=czero
	dykp(k)=czero
	enddo
	enddo
	return
	end
c******************************************************************
	subroutine assolf(rx,pe,l,coul,u,up,uk,ukp)
c*****************************************************************
c     this sbr computes the solution f_l of the coulomb equation *

c     input parameters:                                          *
c     x = radius                                                 *
c     pe = k-value                                               *
c     l = angular momentum                                       *
c     coul = 2*k*eta, where eta is the sommerfeld parameter.     *
c     output parameters:                                         *
c     u = solution                                               *
c     up = d(u)/dr                                               *
c     uk = d(u)/dk                                               *
c     ukp= d(d(u)/dk)/dr                                         *
c*****************************************************************
      implicit double precision(r)
      implicit double complex(a-h,o-q,s-z)
        dimension cgc(0:2),cgcp(0:2),cfc(0:2),cfcp(0:2),sig(0:2)
	data ci/(0.d0,1.d0)/
            zlp=l
            nlp=1
            ifail=1
	    mode1=1
	    ipr=1
	    kfn=0
	alpha=0.5d0*coul
 	eta=alpha/pe
	xrho=pe*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	u=cfc(0)
	up=cfcp(0)*pe

	de=1.d-5*pe
	    mode1=1
	ckk=pe+2.d0*de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	u2=cfc(0)
	up2=cfcp(0)*ckk

	ckk=pe+de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	u1=cfc(0)
	up1=cfcp(0)*ckk

	ckk=pe-de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	v1=cfc(0)
	vp1=cfcp(0)*ckk

	ckk=pe-2.d0*de
 	eta=alpha/ckk
	xrho=ckk*rx
            call wclbes(xrho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	v2=cfc(0)
	vp2=cfcp(0)*ckk
	uk=(v2-8.d0*v1+8.d0*u1-u2)/(12.d0*de)
	ukp=(vp2-8.d0*vp1+8.d0*up1-up2)/(12.d0*de)
	return
	end

c******************************************************************
	subroutine assolh(rx,ck,l,coul,u,up,uk,ukp)
c*****************************************************************
c     this sbr computes the outgoing solution o_l = g_l + i f_l  *
c     of the coulomb equation using an asymptotic series.        *
c     the summation of the asymptotic series is done by a        *
c     summation method (levin) 
c     called routines: clevin
c
c     input parameters:                                          *
c     x = radius                                                 *
c     ck = k-value                                               *
c     l = angular momentum                                       *
c     coul = 2*k*eta, where eta is the sommerfeld parameter.     *
c     output parameters:                                         *
c     u = solution                                               *
c     up = d(u)/dr                                               *
c     uk = d(u)/dk                                               *
c     ukp= d(d(u)/dk)/dr                                         *
c*****************************************************************
      implicit double precision(r)
      implicit double complex(a-h,o-q,s-z)
	parameter (max=30,max1=max-1)
      	dimension arup(0:max),arlo(0:max)
	dimension p(max),q(max)
	dimension dif(max),summ(max)
	common /dsige/ sigmae
        dimension cgc(0:2),cgcp(0:2),cfc(0:2),cfcp(0:2),sig(0:2)
	data reps/1.d-13/,rsmall/1.d-1/
        or=dcmplx(1.d0,0.d0)
        oi=dcmplx(0.d0,1.d0)
	oz=dcmplx(0.d0,0.d0)
	x=rx
	rz=0.d0
	alpha=0.5d0*coul
 	eta=alpha/ck
        z=rx*ck
            zlp=l
            nlp=1
            ifail=1
	    mode1=11
	    ipr=1
	    kfn=0
            call wclbes(z,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     &      kfn,mode1,ifail,ipr)
	rpi2=2.d0*datan(1.d0)
	es=cdexp(oi*(sig(0)-l*rpi2))
	esk=-oi*alpha*es*sigmae/ck**2
	nmax=max
	nm1=nmax-1
	zlog2=cdlog(2.d0*z)
	se=cdexp(oi*(z-eta*zlog2))
	z2=2.d0*oi*z
	beta=or
	do 11 n=1,nmax
	p(n)=or
	do 11 j=1,n
	pj=(oi*eta+l+j)*(oi*eta-l+j-1)/j
	p(n)=p(n)*pj
11	continue
	n=1
	q(n)=oi*(2.d0*oi*eta+2*n-1)/n
	do 12 n=2,nmax
	q(n)=oz
	do j=1,n
	tj=or
	do i=1,n
	if(i.eq.j) then
	fi=oi*(2.d0*oi*eta+2*i-1)/i
	else
	fi=(oi*eta+l+i)*(oi*eta-l+i-1)/i
	endif
	tj=tj*fi
	enddo
	q(n)=q(n)+tj
	enddo
12	continue
	sofn=oz
	do 30 n=1,nm1
	term=p(n)/z2**n
	sofn=sofn+term
	if(cdabs(term).eq.rz) then
	sum1=sofn
	nused1=n
	go to 41
	endif
	tn=p(n+1)/z2**(n+1)
	zofn=tn
c  convergence acceleration by levin	
	inp=n-1
        call clevin(sofn,zofn,beta,inp,arup,arlo,max,estlim)
	if(inp.eq.0) then
	estlio=estlim
	go to 30
	endif
	rres=dreal(estlim)
	ries=dimag(estlim)
	rreso=dreal(estlio)
	rieso=dimag(estlio)
	rdifr =dabs(rres-rreso)
	if(dabs(rres).gt.rsmall) rdifr =dabs((rres-rreso)/rres)
	rdifi =dabs(ries-rieso)
	if(dabs(ries).gt.rsmall) rdifi =dabs((ries-rieso)/ries)
	dif(n)=dcmplx(rdifr,rdifi)
	summ(n)=estlim
        if(rdifr.le.reps.and.rdifi.le.reps) 
     +  then
	nused1=n
        go to 21     
	endif
	estlio=estlim
30 	continue
	write(*,*) 'tol not achieved'
	nused1=nm1
21	rdifmin=cdabs(dif(2))
	nmin=2
	do n=3,nused1
	if(cdabs(dif(n)).lt.rdifmin) then
	nmin=n
	rdifmin=cdabs(dif(n))
	endif
	enddo
	sum1=summ(nmin)
	nused1=nmin
41	u=se*(or+sum1)
*********************************************************
	sofn=oz
	do 31 n=1,nm1
	term=-n*p(n)/z2**(n+1)
	sofn=sofn+term
	if(cdabs(term).eq.rz) then
	sum2=sofn
	nused2=n
	go to 42
	endif
	tn=-(n+1)*p(n+1)/z2**(n+2)
	zofn=tn
	inp=n-1
c  convergence acceleration by levin	
        call clevin(sofn,zofn,beta,inp,arup,arlo,max,estlim)
	if(inp.eq.0) then
	estlio=estlim
	go to 31
	endif
	rres=dreal(estlim)
	ries=dimag(estlim)
	rreso=dreal(estlio)
	rieso=dimag(estlio)
	rdifr =dabs(rres-rreso)
	if(dabs(rres).gt.rsmall) rdifr =dabs((rres-rreso)/rres)
	rdifi =dabs(ries-rieso)
	if(dabs(ries).gt.rsmall) rdifi =dabs((ries-rieso)/ries)
	dif(n)=dcmplx(rdifr,rdifi)
	summ(n)=estlim
        if(rdifr.le.reps.and.rdifi.le.reps) 
     +  then
	nused2=n
c	sum2=estlim
        go to 22     
	endif
	estlio=estlim
31 	continue
	write(*,*) 'tol not achieved'
	nused2=nm1
c	sum2=estlim
22	rdifmin=cdabs(dif(2))
	nmin=2
	do n=3,nused2
	if(cdabs(dif(n)).lt.rdifmin) then
	nmin=n
	rdifmin=cdabs(dif(n))
	endif
	enddo
	sum2=summ(nmin)
	nused2=nmin
42	up=se*((or+sum1)*oi*(ck-eta/x)+2.d0*oi*ck*sum2)
*******************************************************
	sofn=oz
	do 32 n=1,nm1
	term=q(n)/z2**n
	sofn=sofn+term
	if(cdabs(term).eq.rz) then
	sum3=sofn
	nused3=n
	go to 43
	endif
	tn=q(n+1)/z2**(n+1)
	zofn=tn
c  convergence acceleration by levin	
	inp=n-1
        call clevin(sofn,zofn,beta,inp,arup,arlo,max,estlim)
	if(inp.eq.0) then
	estlio=estlim
	go to 32
	endif
	rres=dreal(estlim)
	ries=dimag(estlim)
	rreso=dreal(estlio)
	rieso=dimag(estlio)
	rdifr =dabs(rres-rreso)
	if(dabs(rres).gt.rsmall) rdifr =dabs((rres-rreso)/rres)
	rdifi =abs(ries-rieso)
	if(dabs(ries).gt.rsmall) rdifi =dabs((ries-rieso)/ries)
	dif(n)=dcmplx(rdifr,rdifi)
	summ(n)=estlim
        if(rdifr.le.reps.and.rdifi.le.reps) 
     +  then
	nused3=n
c	sum3=estlim
        go to 23     
	endif
	estlio=estlim
32 	continue
	write(*,*) 'tol not achieved'
	nused3=nm1
c	sum3=estlim
23	rdifmin=cdabs(dif(2))
	nmin=2
	do n=3,nused3
	if(cdabs(dif(n)).lt.rdifmin) then
	nmin=n
	rdifmin=cdabs(dif(n))
	endif
	enddo
	sum3=summ(nmin)
	nused3=nmin
43	ueta=-oi*zlog2*u+se*sum3
	uk=x*up/ck-alpha*ueta/ck**2
*********************************************************
	sofn=oz
	do 33 n=1,nm1
	term=-n*q(n)/z2**(n+1)
	sofn=sofn+term
	if(cdabs(term).eq.rz) then
	sum4=sofn
	nused4=n
	go to 44
	endif
	tn=-(n+1)*q(n+1)/z2**(n+2)
	zofn=tn
	inp=n-1
c  convergence acceleration by levin	
        call clevin(sofn,zofn,beta,inp,arup,arlo,max,estlim)
	if(inp.eq.0) then
	estlio=estlim
	go to 33
	endif
	rres=dreal(estlim)
	ries=dimag(estlim)
	rreso=dreal(estlio)
	rieso=dimag(estlio)
	rdifr =dabs(rres-rreso)
	if(dabs(rres).gt.rsmall) rdifr =dabs((rres-rreso)/rres)
	rdifi =abs(ries-rieso)
	if(dabs(ries).gt.rsmall) rdifi =dabs((ries-rieso)/ries)
	dif(n)=dcmplx(rdifr,rdifi)
	summ(n)=estlim
        if(rdifr.le.reps.and.rdifi.le.reps) 
     +  then
	nused4=n
c	sum4=estlim
        go to 24     
	endif
	estlio=estlim
33 	continue
	write(*,*) 'tol not achieved'
	nused4=nm1
c	sum4=estlim
24	rdifmin=cdabs(dif(2))
	nmin=2
	do n=3,nused4
	if(abs(dif(n)).lt.rdifmin) then
	nmin=n
	rdifmin=cdabs(dif(n))
	endif
	enddo
	sum4=summ(nmin)
	nused4=nmin
44      ukp=(oi*alpha/(ck*z)-z+2.d0*alpha/ck+l*(l+1)/z)*u+
     +        up/ck+oi*alpha*zlog2*up/ck**2-
     +        se*alpha*oi*(ck-eta/x)*sum3/ck**2- 
     +        se*2.d0*oi*alpha*sum4/ck
*******************************************************
	uk=es*uk+u*esk
	ukp=es*ukp+up*esk
	u=u*es
	up=up*es
	return
	end
c******************************************************************

************************************************************************
c   subroutine clevin(sofn,rofn,beta,n,arup,arlo,larray,estlim)
************************************************************************
*
*  description :
*  -------------
*
*  subroutine clevin computes levin's transformation.
*  (d. levin (j. comput. math., vol. b3, 1973, 371-388 ) using the
*  partial sums s(n) and terms a(n) of a given series as input data.
*
*  here, the definition of levin's sequence transformation according to
*  eq. (7.1-7) of e. j. weniger, comput. phys. rep., vol. 10 (1989),
*  189 - 371, is used.
*
*  in levin's sequence transformation it is assumed that the partial
*  sums s(n) of the series to be transformed can be written as
*
*  s(n) = s + r(n),
*
*  with s being the limit or antilimit of this series, and that the
*  exact remainder r(n) can be approximated by a remainder estimate
*  omega(n) in such a way that the ratio
*
*  [ s(n) - s ] / omega(n)
*
*  can be expressed as a by a poincare-type asymptotic expansion in
*  inverse powers of (n + beta).
*
*  levin suggested some simple remainder estimates omega(n) which can
*  be computed be determined from at most two terms a(n) of the series
*  to be transformed. with the help of these remainder estimates the
*  following variants of levin's sequence transformations result
*  (compare section 7.3 of e. j. weniger, comput. phys. rep.,
*  vol. 10 (1989), 189):
*
*  (n+beta) * a(n)                     : u transformation, eq. (7.3-5)
*  a(n)                                : t transformation, eq. (7.3-7)
*  a(n+1)                              : d transformation, eq. (7.3-7)
*  a(n + 1) * a(n) / [a(n) - a(n + 1)] : v transformation, eq. (7.3-11)
*
*  it is always assumed that the index n of the terms and the partial
*  sums s(n) satisfies n >= 0.
*
c  both numerator and denominator of this transformation are calculated 
c  with the help of a simple extension of a 2-dimensional nonlinear
c  3-term recurrence formula which was originally derived by 
c  t.fessler, w.f.ford and w.f.smith 
c  ( acm trans.math.software ,vol. 9 ( 1983 ), pp.346-354 ).           *
c                                                                      *
c  in this subroutine a variant of the so-called moving lozenge tech-  *
c  nique ( p.wynn, r.f.t.i.-chiffres, vol.9, pp.327-362 ( 1966 )) is   *
c  used, i.e., only one counterdiagonal of the numerator and denomina- *
c  tor tables have to be stored at a time. consequently, only two      *
c  1-dimensional arrays of appropriate length are needed in this       *
c  subroutine.                                                         *
c  in each call this subroutine calculates new counterdiagonals of the *
c  numerator and denominator tables and overwrites the previous entries.
c  in addition, a new estimate of the limit is calculated which is the *
c  in addition, a new estimate of the limit is calculated which is the *
c  ratio of the array elements arup(0) and arlo(0). hence, for every   *
c  new sequence element s-of-n this subroutine has to be called again.
c                                                                      *
c  input parameters :                                                  *
c  ----------------                                                    *
c                                                                      *
c  sofn   : element s(n) of the sequence to be accelerated             *
c                                                                      *
c  rofn   : estimate of the leading term of the poincare-type asympto- *
c           tic expansion of the remainder r(n) in inverse powers of 
c           (n+beta), i.e. omega(n).    
c           it is tacitly assumed that rofn is different from zero.    *
c                                                                      *
c  beta   : parameter contained in the inverse powers                  *
c           ( n + beta ) ** ( -j )                                     *
c           which describe the higher order contributions in the       *
c           poincare-type expansion of the remainder r(n).             *
c                                                                      *
c  n      : number of the last sequence element (n has to satisfy the  *
c           inequality n .ge. 0 )                                      *
c                                                                      *
c  arup   : 1-dimensional array to store the actual counterdiagonal of *
c           the numerator table.                                       *
c                                                                      *
c  arlo   : 1-dimensional array to store the actual counterdiagonal of *
c           the denominator table.                                     *
c                                                                      *
c  larray : length of the 1-dimensional arrays arup and arlo.          *
c                                                                      *
c  output parameter :                                                  *
c  ----------------                                                    *
c                                                                      *
c  estlim : approximation to the limit of the sequence s(n) which is   *
c           to be accelerated.                                         *
c           we have here estlim = arup(0)/arlo(0).                     *
c                                                                      *
c  machine-dependent parameters:                                       *
c  ----------------------------                                        *
c                                                                      *
c  huge,tiny : huge should be set close to but not identical with the  *
c              largest floating point number representable on the com- *
c              puter, and tiny should be set close to but not identical*
c              with the smallest floating point number representable   *
c              on the computer.                                        *
c                                                                      *
c***********************************************************************
      subroutine clevin(sofn,rofn,beta,n,arup,arlo,larray,estlim)
      implicit double complex(a-h,o-z)
      double precision huge, tiny, one,r00
      dimension arup(0:larray),arlo(0:larray)
      parameter ( huge = 1.d+30 , tiny = 1.d-30 , one = 1.d0 )
      data r00/0.d0/
      zone=dcmplx(one,r00)
c
c  check whether n is less than zero
c
      if (n.lt.0) then
        write(*,1000)
1000    format('*** error in subroutine clevin ***'/
     1         '=================================')
        write(6,1010) n
1010    format('illegal input parameter n which must be .ge. zero'/
     1         'here we have n = ',i4)
        stop
c
c  check whether the arrays arlo and arup are large enough to store
c  all the elements required
c
      else if (n.gt.larray) then
        write(*,1000)
        write(*,1020) n,larray
1020    format('illegal input parameters n and larray'/
     1  'we must have n .le. larray'/
     2  'here we have n = ',i4,' and larray = ',i4)
        return
      end if
c
c  this is reached if the general levin transformation has to be
c  calculated.
c
c  first it has to be checked whether n .ge. 1 holds, i.e., whether we
c  proceed with the calculation of the next counter diagonals of the
c  two arrays or whether we have to initialize a new calculation.
c  this is done if n = 0 holds.
c
      arup(n)   = sofn / rofn
      arlo(n)   = zone  / rofn
c
c  check whether n .gt. 0 holds. in that case the recursive computation
c  of the numerators and denominators of the general levin transformatio
c  is to be continued.
c
      if ( n .gt. 0 ) then
        arup(n-1) = arup(n) - arup(n-1)
        arlo(n-1) = arlo(n) - arlo(n-1)
c
c  check whether n .gt. 1 holds. in that case the actual recursion
c  of the numerators and denominators starts.
c
        if (n.gt.1) then
          bn1     = beta + float(n-1)
          bn2     = beta + float(n)
          coef    = bn1 / bn2
          do 10 j = 2,n
            fact  = (beta+float(n-j)) * coef**(j-2) / bn2
c
c  perform the calculation of the new elements of the two tables.
c
            arup(n-j)  = arup(n-j+1) - fact*arup(n-j)
            arlo(n-j)  = arlo(n-j+1) - fact*arlo(n-j)
10        continue
        end if
      end if
c
c  compute the estimate for the limit of the sequence.
c
        if (abs(arlo(0)).lt.tiny) then
          estlim = dcmplx(huge,huge)
        else
          estlim = arup(0)/arlo(0)
        end if
      return
      end
      subroutine ludcmp(a,n,indx,d,ierror)
      implicit double complex(a-h,o-q,s-z)
      implicit double precision (r)
      parameter (nmax=2000,tiny=(1.0d-30,1.0d-30))
      dimension a(n,n),indx(n),rvv(nmax)
      d=dcmplx(1.d0,0.d0)
	ierror=0
      do 12 i=1,n
        raamax=0.d0
        do 11 j=1,n
          if (abs(a(i,j)).gt.raamax) raamax=abs(a(i,j))
11      continue
        if (raamax.eq.0.d0)then
c	write (*,*) 'singular matrix in ludcmp'
	ierror=1
	return
	endif
        rvv(i)=1.d0/raamax
12    continue
      do 19 j=1,n
        if (j.gt.1) then
          do 14 i=1,j-1
            sum=a(i,j)
            if (i.gt.1)then
              do 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
13            continue
              a(i,j)=sum
            endif
14        continue
        endif
        raamax=0.d0
        do 16 i=j,n
          sum=a(i,j)
          if (j.gt.1)then
            do 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
15          continue
            a(i,j)=sum
          endif
          rdum=rvv(i)*abs(sum)
          if (rdum.ge.raamax) then
            imax=i
            raamax=rdum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          rvv(imax)=rvv(j)
        endif
        indx(j)=imax
        if(j.ne.n)then
          if(abs(a(j,j)).eq.0.d0) a(j,j)=tiny
          dum=1.d0/a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      if(abs(a(n,n)).eq.0.d0) a(n,n)=tiny
      return
      end

      subroutine lubksb(a,n,indx,b)
      implicit double complex(a-h,o-q,s-z)
      implicit double precision (r)
      dimension a(n,n),indx(n),b(n)
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.d0) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        if(i.lt.n)then
          do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
13        continue
        endif
        b(i)=sum/a(i,i)
14    continue
      return
      end

      function potnuc(i,j,r)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      include  'params'
        common /eigend/ depth 
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)

	data ci/(0.d0,1.d0)/
        if(i.eq.j) then
        pnsf=r2muhb2*vnd(r)
        else
        pnsf=0.d0
        endif
	call cmatd(r,cpot)
	potnuc=pnsf+r2muhb2*cpot(i,j)
      return
      end

      function pott(i,j,r)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      real*8 pcoul
      include 'params'
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
	common/poti1/v1c(ndim,ndim),v10(ndim,ndim),v11(ndim,ndim),
     &  v12(ndim,ndim),lo(ndim)   
	if(i.eq.j) then
	porb=lo(i)*(lo(i)+1)/r**2
	pnsf=r2muhb2*vn(r)
        pcsf=r2muhb2*pcoul(r)
	else
	porb=0.d0
	pnsf=0.d0
	pcsf=0.d0
	endif
	call cmat(r,cpot)
	pott=porb+pnsf+pcsf+r2muhb2*cpot(i,j)
      return
      end

      function potu(i,j,r)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      real*8 pcoul
      include 'params'
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
	common/poti1/v1c(ndim,ndim),v10(ndim,ndim),v11(ndim,ndim),
     &  v12(ndim,ndim),lo(ndim)   
        common /valene/ ce(ndim)
	if(i.eq.j) then
	porb=lo(i)*(lo(i)+1)/r**2
        pcsf=r2muhb2*pcoul(r)
	pene=-ce(i)
	else
	porb=0.d0
	pcsf=0.d0
	pene=0.d0
	endif
	call cmatu(r,cpot)
	potu=pene+porb+pcsf+r2muhb2*cpot(i,j)
        return
        end

      function pot(i,j,r)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      real*8 cpcoul
      include 'params'
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
	common/poti1/v1c(ndim,ndim),v10(ndim,ndim),v11(ndim,ndim),
     &  v12(ndim,ndim),lo(ndim)   
	if(i.eq.j) then
	pnsf=r2muhb2*vn(r)
        pcsf=r2muhb2*cpcoul(r)
	else
	pnsf=0.d0
	pcsf=0.d0
	endif
	call cmatc(r,cpot)
	pot=pnsf+pcsf+r2muhb2*cpot(i,j)
      return
      end

      function spot(i,j,r)
      implicit complex*16(a-h,o-q,s-z)
      implicit real*8 (r)
      real*8 spcoul
      include 'params'
      parameter (nlevelmax=25)
      dimension cpot(nlevelmax,nlevelmax)
        common /param/ r2muhb2,rv0,rw0,radp,rbeta,rdif,r0
	if(i.eq.j) then
        pcsf=r2muhb2*spcoul(r)
	else
	pcsf=0.d0
	endif
	call cmats(r,cpot)
	spot=pcsf+r2muhb2*cpot(i,j)
      return
      end

      subroutine ident(i,li,mi)
      include 'params'
	common /lorb/ lang(ndim)
      li=lang(i)
      return
      end

c********************************************************************
      subroutine wclbes(zz,eta1,zlmin,nl,fc,gc,fcp,gcp,sig,kfn,mode1,
     1                  ifail,ipr)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c  complex coulomb wavefunction program using steed's method           c
c  original title : coulcc                                             c
c                                                                      c
c  a. r. barnett           manchester  march   1981                    c
c  modified i.j. thompson  daresbury, sept. 1983 for complex functions c
c                                                                      c
c  the form (not the substance) of this program has been modified      c
c   by k.s. kolbig (cern)    december 1987                             c
c                                                                      c
c  original program  rcwfn       in    cpc  8 (1974) 377-395           c
c                 +  rcwff       in    cpc 11 (1976) 141-142           c
c                 +  coulfg      in    cpc 27 (1982) 147-166           c
c  description of real algorithm in    cpc 21 (1981) 297-314           c
c  description of complex algorithm    jcp 64 (1986) 490-509           c
c  this version written up       in    cpc 36 (1985) 363-372           c
c                                                                      c
c  wclbes returns f,g,g',g',sig for complex eta, zz, and zlmin,        c
c   for nl integer-spaced lambda values zlmin to zlmin+nl inclusive,   c
c   thus giving  complex-energy solutions to the coulomb schrodinger   c
c   equation,to the klein-gordon equation and to suitable forms of     c
c   the dirac equation ,also spherical & cylindrical bessel equations  c
c                                                                      c
c  if abs(mode1)                                                       c
c            = 1  get f,g,f',g'   for integer-spaced lambda values     c
c            = 2      f,g      unused arrays must be dimensioned in    c
c            = 3      f,  f'          call to at least length (0:1)    c
c            = 4      f                                                c
c            = 11 get f,h+,f',h+' ) if kfn=0, h+ = g + i.f        )    c
c            = 12     f,h+        )       >0, h+ = j + i.y = h(1) ) in c
c            = 21 get f,h-,f',h-' ) if kfn=0, h- = g - i.f        ) gc c
c            = 22     f,h-        )       >0, h- = j - i.y = h(2) )    c
c                                                                      c
c     if mode1 < 0 then the values returned are scaled by an exponen-  c
c                  tial factor (dependent only on zz) to bring nearer  c
c                  unity the functions for large abs(zz), small eta ,  c
c                  and abs(zl) < abs(zz).                              c
c        define scale = (  0        if mode1 > 0                       c
c                       (  imag(zz) if mode1 < 0  &  kfn < 3           c
c                       (  real(zz) if mode1 < 0  &  kfn = 3           c
c        then fc = exp(-abs(scale)) * ( f, j, j, or i)                 c
c         and gc = exp(-abs(scale)) * ( g, y, or y )                   c
c               or exp(scale)       * ( h+, h(1), or k)                c
c               or exp(-scale)      * ( h- or h(2) )                   c
c                                                                      c
c  if  kfn  =  0,-1  complex coulomb functions are returned   f & g    c
c           =  1   spherical bessel      "      "     "       j & y    c
c           =  2 cylindrical bessel      "      "     "       j & y    c
c           =  3 modified cyl. bessel    "      "     "       i & k    c
c                                                                      c
c          and where coulomb phase shifts put in sig if kfn=0 (not -1) c
c                                                                      c
c  the use of mode and kfn is independent                              c
c    (except that for kfn=3,  h(1) & h(2) are not given)               c
c                                                                      c
c  with negative orders lambda, wclbes can still be used but with      c
c    reduced accuracy as cf1 becomes unstable. the user is thus        c
c    strongly advised to use reflection formulae based on              c
c    h+-(zl,,) = h+-(-zl-1,,) * exp +-i(sig(zl)-sig(-zl-1)-(zl+1/2)pi) c
c                                                                      c
c  precision:  results to within 2-3 decimals of 'machine accuracy',   c
c              except in the following cases:                          c
c              (1) if cf1a fails because x too small or eta too large  c
c               the f solution  is less accurate if it decreases with  c
c               decreasing lambda (e.g. for lambda.le.-1 & eta.ne.0)   c
c              (2) if eta is large (e.g. >> 50) and x inside the       c
c                turning point, then progressively less accuracy       c
c              (3) if zlmin is around sqrt(accur) distance from an     c
c               integral order and abs(x) << 1, then errors present.   c
c               rerr traces the main roundoff errors.                  c
c                                                                      c
c   wclbes is coded for real*8 on ibm or equivalent  accur >= 10**-14  c
c          with a section of doubled real*16 for less roundoff errors. c
c          (if no doubled precision available, increase jmax to eg 100)c
c   use implicit complex*32 & real*16 on vs compiler accur >= 10**-32  c
c   for single precision cdc (48 bits) reassign                        c
c        double precision=real etc.                                    c
c                                                                      c
c   ipr    on input   = 0 : no printing of error messages              c
c                    ne 0 : print error messages on file 6             c
c   ifail  in output = -2 : argument out of range                      c
c                    = -1 : one of the continued fractions failed,     c
c                           or arithmetic check before final recursion c
c                    =  0 : all calculations satisfactory              c
c                    ge 0 : results available for orders up to & at    c
c                             position nl-ifail in the output arrays.  c
c                    = -3 : values at zlmin not found as over/underflowc
c                    = -4 : roundoff errors make results meaningless   c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c     machine dependent constants :                                    c
c                                                                      c
c     accu     target bound on relative error (except near 0 crossings)c
c               (accur should be at least 100 * acc8)                  c
c     acc8    smallest number with 1+acc8 .ne.1 in real*8  arithmetic  c
c     acc16    smallest number with 1+acc16.ne.1 in real*16 arithmetic c
c     fpmax    magnitude of largest floating point number * acc8       c
c     fpmin    magnitude of smallest floating point number / acc8      c
c                                                                      c
c     parameters determining region of calculations :                  c
c                                                                      c
c        r20      estimate of (2f0 iterations)/(cf2 iterations)        c
c        asym     minimum x/(eta**2+l) for cf1a to converge easily     c
c        xnear    minimum abs(x) for cf2 to converge accurately        c
c        limit    maximum no. iterations for cf1, cf2, and 1f1 series  c
c        jmax     size of work arrays for pade accelerations           c
c        ndrop    number of successive decrements to define instabilityc
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     c309r1 = cf2,   c309r2 = f11,    c309r3 = f20,
c     c309r4 = cf1r,  c309r5 = cf1c,   c309r6 = cf1a,
c     c309r7 = rcf,   c309r8 = ctidy.
c
      implicit complex*16 (a-h,o-z)
c
      common /coulc2/ nfp,n11,npq1,npq2,n20,kas1,kas2
      double precision zero,one,two,half
      double precision r20,asym,xnear
      double precision fpmax,fpmin,fplmax,fplmin
      double precision accu,acc8,acc16
      double precision hpi,tlog
      double precision err,rerr,absc,accur,acct,acch,accb,c309r4
      double precision paccq,eps,off,scale,sf,sfsh,ta,rk,omega,absx
 
      logical lpr,etane0,ifcp,rlel,donem,unstab,zlneg,axial,nocf2,npint
 
      parameter(zero = 0, one = 1, two = 2, half = one/two)
      parameter(ci = (0,1), cih = half*ci)
      parameter(r20 = 3, asym = 3, xnear = half)
      parameter(limit = 20000, ndrop = 5, jmax = 50)
 
      parameter(fpmax = 3d290, fpmin = 2d-300)
      parameter(accu = 1d-14, acc8 = 2d-16, acc16 = 2d-16)
 
      parameter(hpi  = 1.57079 63267 94896 619d0)
      parameter(tlog = 0.69314 71805 59945 309d0)
 
      dimension fc(0:*),gc(0:*),fcp(0:*),gcp(0:*),sig(0:*)
      dimension xrcf(jmax,4)
      integer npq(2),kas(2)
      equivalence (npq(1),npq1),(npq(2),npq2)
      equivalence (kas(1),kas1),(kas(2),kas2)
******************************************************************
	common /dsige/ sigmae
 
      nintc(w)=nint(dreal(w))
      absc(w)=abs(dreal(w))+abs(dimag(w))
      npint(w,accb)=absc(nintc(w)-w) .lt. accb .and. dreal(w) .lt. half
c
c	write(*,*) 'z :',zz
c	pause
      mode=mod(abs(mode1),10)
      ifcp=mod(mode,2) .eq. 1
      lpr=ipr .ne. 0
      ifail=-2
      n11=0
      nfp=0
      kas(1)=0
      kas(2)=0
      npq(1)=0
      npq(2)=0
      n20=0
 
      accur=max(accu,50*acc8)
      acct=half*accur
      acch=sqrt(accur)
      accb=sqrt(acch)
      rerr=acct
      fplmax=log(fpmax)
      fplmin=log(fpmin)
c
      cik=one
      if(kfn .ge. 3) cik=ci*sign(one,fpmin-dimag(zz))
      x=zz*cik
      eta=eta1
      if(kfn .gt. 0) eta=zero
      etane0=absc(eta) .gt. acc8
      etai=eta*ci
      dell=zero
      if(kfn .ge. 2) dell=half
      zm1=zlmin-dell
      scale=zero
      if(mode1 .lt. 0) scale=dimag(x)
c
      m1=1
      l1=m1+nl
      rlel=abs(dimag(eta))+abs(dimag(zm1)) .lt. acc8
      absx=abs(x)
      axial=rlel .and. abs(dimag(x)) .lt. acc8*absx
      if(mode .le. 2 .and. absx .lt. fpmin) go to 310
      xi=one/x
      xlog=log(x)
 
c       log with cut along the negative real axis, see also omega
 
      id=1
      donem=.false.
      unstab=.false.
      lf=m1
      ifail=-1
   10 zlm=zm1+(lf-m1)
      zll=zm1+(l1-m1)
c
c ***       zll  is final lambda value, or 0.5 smaller for j,y bessels
c
      z11=zll
      if(id .lt. 0) z11=zlm
      p11=ci*sign(one,acc8-dimag(eta))
      last=l1
c
c ***       find phase shifts and gamow factor at lambda = zll
c
      pk=zll+one
      aa=pk-etai
      ab=pk+etai
      bb=two*pk
      zlneg=npint(bb,accb)
      clgaa=wlogam(aa)
      cpsaa=-ci*wdigam(aa)
      clgab=clgaa
      cpsab=cpsaa
      if(etane0 .and. .not.rlel) then
      clgab=wlogam(ab)
      cpsab=ci*wdigam(ab)
      endif
      if(etane0 .and. rlel) then
      clgab=dconjg(clgaa)
      cpsab=dconjg(cpsaa)
      endif
      sigma=(clgaa-clgab)*cih
      sigmae=(cpsaa-cpsab)*cih
      if(kfn .eq. 0) sig(l1)=sigma
      if(.not.zlneg) cll=zll*tlog-hpi*eta-wlogam(bb)+(clgaa+clgab)*half
      theta=x-eta*(xlog+tlog)-zll*hpi+sigma
c
      ta=(dimag(aa)**2+dimag(ab)**2+abs(dreal(aa))+abs(dreal(ab)))*half
      if(id .gt. 0 .and. absx .lt. ta*asym .and. .not.zlneg) go to 20
c
c ***         use cf1 instead of cf1a, if predicted to converge faster,
c                 (otherwise using cf1a as it treats negative lambda &
c                  recurrence-unstable cases properly)
c
      rk=sign(one,dreal(x)+acc8)
      p=theta
      if(rk .lt. 0) p=-x+eta*(log(-x)+tlog)-zll*hpi-sigma
      xrcf(1,1)=pk
      f=rk*c309r6(x*rk,eta*rk,zll,p,acct,jmax,nfp,fest,err,fpmax,xrcf,
     1                                      xrcf(1,3),xrcf(1,4))
      fesl=log(fest)+abs(dimag(x))
      nfp=-nfp
      if(nfp .lt. 0 .or. unstab .and. err .lt. accb) go to 40
      if(.not.zlneg .or. unstab .and. err .gt. accb) go to 20
      if(lpr) write(6,1060) '-l',err
      if(err.gt.accb) go to 280
      go to 40
c
c ***    evaluate cf1  =  f   =  f'(zll,eta,x)/f(zll,eta,x)
c
   20 if(axial) then
	r_x=x
	r_eta=eta
	r_zll=zll
       f=c309r4(r_x,r_eta,r_zll,acc8,sf,rk,etane0,limit,err,nfp,
     1          fpmin,fpmax,lpr)
       fcl=sf
       tpk1=rk
      else
       f=c309r5(x,eta,zll,acc8,fcl,tpk1,etane0,limit,err,nfp,
     1          fpmin,fpmax,lpr)
      end if
      if(err .gt. one) then
       ifail=-1
       go to 290
      end if
c
c ***  make a simple check for cf1 being badly unstable:
c
      if(id .ge. 0) then
       unstab=dreal((one-eta*xi)*ci*dimag(theta)/f) .gt. zero
     1  .and. .not.axial .and. abs(dimag(theta)) .gt. -log(acc8)*half
     2  .and. absc(eta)+absc(zll) .lt. absc(x)
       if(unstab) then
        id=-1
        lf=l1
        l1=m1
        rerr=acct
        go to 10
       end if
      end if
c
c *** compare accumulated phase fcl with asymptotic phase for g(k+1) :
c     to determine estimate of f(zll) (with correct sign) to start recur
c
      w=x*x*(half/tpk1+one/tpk1**2)+eta*(eta-two*x)/tpk1
      fesl=(zll+one)*xlog+cll-w-log(fcl)
   40 fesl=fesl-abs(scale)
      rk=max(dreal(fesl),fplmin*half)
      fesl=dcmplx(min(rk,fplmax*half),dimag(fesl))
      fest=exp(fesl)
c
      rerr=max(rerr,err,acc8*abs(dreal(theta)))
c
      fcl=fest
      fpl=fcl*f
      if(ifcp) fcp(l1)=fpl
      fc(l1)=fcl
c
c *** downward recurrence to lambda = zlm. array gc,if present,stores rl
c
      i=max(-id,0)
      zl=zll+i
      mono=0
      off=abs(fcl)
      ta=absc(sigma)
 
c
c     corresponds to   do 70 l = l1-id,lf,-id
c
      l=l1-id
      lc70=(l1-lf)/id
   70 if(lc70 .le. 0) go to 80
      if(etane0) then
       if(rlel) then
        dsig=atan2(dreal(eta),dreal(zl))
	dsige=dreal(zl)/(dreal(zl)**2+dreal(eta)**2)
        rl=sqrt(dreal(zl)**2+dreal(eta)**2)
       else
        aa=zl-etai
        bb=zl+etai
        if(absc(aa) .lt. acch .or. absc(bb) .lt. acch) goto 50
        dsig=(log(aa)-log(bb))*cih
	dsige=half*(one/aa+one/bb)
        rl=aa*exp(ci*dsig)
       end if
       if(absc(sigma) .lt. ta*half) then
 
c               re-calculate sigma because of accumulating roundoffs:
 
        sl=(wlogam(zl+i-etai)-wlogam(zl+i+etai))*cih
	cpsaai=wdigam(zl+i-etai)
	cpsabi=wdigam(zl+i+etai)
	sle=half*(cpsaai+cpsabi)
        rl=(zl-etai)*exp(ci*id*(sigma-sl))
        sigma=sl
	sigmae=sle
        ta=zero
       else
        sigma=sigma-dsig*id
	sigmae=sigmae-dsige*id
       end if
       ta=max(ta,absc(sigma))
       sl=eta+zl*zl*xi
       pl=zero
       if(absc(zl) .gt. acch) pl=(sl*sl-rl*rl)/zl
       fcl1=(fcl*sl+id*zl*fpl)/rl
       sf=abs(fcl1)
       if(sf .gt. fpmax) go to 350
       fpl=(fpl*sl+id*pl*fcl)/rl
       if(mode .le. 1) gcp(l+id)=pl*id
      else
 
c                      eta = 0, including bessels.  nb rl==sl
 
       rl=zl*xi
       fcl1=fcl*rl+fpl*id
       sf=abs(fcl1)
       if(sf .gt. fpmax) go to 350
       fpl=(fcl1*rl-fcl)*id
      end if
      mono=mono+1
      if(sf .ge. off) mono=0
      fcl=fcl1
      off=sf
      fc(l)=fcl
      if(ifcp) fcp(l)=fpl
      if(kfn .eq. 0) sig(l)=sigma
      if(mode .le. 2) gc(l+id)=rl
      zl=zl-id
      if(mono .lt. ndrop .or. axial .or.
     1             dreal(zlm)*id .gt. -ndrop .and. .not.etane0) then
       l=l-id
       lc70=lc70-1
       go to 70
      end if
      unstab=.true.
c
c ***    take action if cannot or should not recur below this zl:
 
   50 zlm=zl
      lf=l
      if(id .lt. 0) go to 380
      if(.not.unstab) lf=l+1
      if(l+mono .lt. l1-2 .or. id .lt. 0 .or. .not.unstab) go to 80
 
c             otherwise, all l values (for stability) should be done
c                        in the reverse direction:
 
      id=-1
      lf=l1
      l1=m1
      rerr=acct
      go to 10
 
   80 if(fcl .eq. zero) fcl=acc8
      f=fpl/fcl
c
c *** check, if second time around, that the 'f' values agree!
c
      if(id .gt. 0) first=f
      if(donem) rerr=max(rerr,absc(f-first)/absc(f))
      if(donem) go to 90
c
      nocf2=.false.
      thetam=x-eta*(xlog+tlog)-zlm*hpi+sigma
c
c *** on left x-plane, determine omega by requiring cut on -x axis
c     on right x-plane, choose omega (using estimate based on thetam)
c       so h(omega) is smaller and recurs upwards accurately.
c     (x-plane boundary is shifted to give cf2(lh) a chance to converge)
c
      omega=sign(one,dimag(x)+acc8)
      if(dreal(x) .ge. xnear) omega=sign(one,dimag(thetam)+acc8)
      sfsh=exp(omega*scale-abs(scale))
      off=exp(min(two*max(abs(dimag(x)),abs(dimag(thetam)),
     1                         abs(dimag(zlm))*3),fplmax))
      eps=max(acc8,acct*half/off)
c
c ***    try first estimated omega, then its opposite,
c        to find the h(omega) linearly independent of f
c        i.e. maximise  cf1-cf2 = 1/(f h(omega)) , to minimise h(omega)
c
   90 do 100 l = 1,2
      lh=1
      if(omega .lt. zero) lh=2
      pm=ci*omega
      etap=eta*pm
      if(donem) go to 130
      pq1=zero
      paccq=one
      kase=0
c
c ***            check for small x, i.e. whether to avoid cf2 :
c
      if(mode .ge. 3 .and. absx .lt. one ) go to 190
      if(mode .lt. 3 .and. (nocf2 .or. absx .lt. xnear .and.
     1   absc(eta)*absx .lt. 5 .and. absc(zlm) .lt. 4)) then
       kase=5
       go to 120
      end if
c
c ***  evaluate   cf2 : pq1 = p + i.omega.q  at lambda = zlm
c
      pq1=c309r1(x,eta,zlm,pm,eps,limit,err,npq(lh),acc8,acch,
     1             lpr,accur,dell)
      err=err*max(one,absc(pq1)/max(absc(f-pq1),acc8))
c
c *** check if impossible to get f-pq accurately because of cancellation
c                original guess for omega (based on thetam) was wrong
c                use kase 5 or 6 if necessary if re(x) < xnear
      if(err .lt. acch) go to 110
      nocf2=dreal(x) .lt. xnear .and. abs(dimag(x)) .lt. -log(acc8)
  100 omega=-omega
      if(unstab) go to 360
      if(lpr .and. dreal(x) .lt. -xnear) write(6,1060) '-x',err
  110 rerr=max(rerr,err)
c
c ***  establish case of calculation required for irregular solution
c
  120 if(kase .ge. 5) go to 130
      if(dreal(x) .gt. xnear) then
 
c          estimate errors if kase 2 or 3 were to be used:
 
       paccq=eps*off*absc(pq1)/max(abs(dimag(pq1)),acc8)
      end if
      if(paccq .lt. accur) then
       kase=2
       if(axial) kase=3
      else
       kase=1
       if(npq(1)*r20 .lt. jmax) kase=4
c             i.e. change to kase=4 if the 2f0 predicted to converge
      end if
  130 go to (190,140,150,170,190,190), abs(kase)
  140 if(.not.donem)
c
c ***  evaluate   cf2 : pq2 = p - i.omega.q  at lambda = zlm   (kase 2)
c
     1  pq2=c309r1(x,eta,zlm,-pm,eps,limit,err,npq(3-lh),acc8,acch,
     2             lpr,accur,dell)
      p=(pq2+pq1)*half
      q=(pq2-pq1)*half*pm
      go to 160
  150 p=dreal(pq1)
      q=dimag(pq1)
c
c ***   with kase = 3 on the real axes, p and q are real & pq2 = pq1*
c
      pq2=dconjg(pq1)
c
c *** solve for fcm = f at lambda = zlm,then find norm factor w=fcm/fcl
c
  160 w=(pq1-f)*(pq2-f)
      sf=exp(-abs(scale))
      fcm=sqrt(q/w)*sf
 
c                  any sqrt given here is corrected by
c                  using sign for fcm nearest to phase of fcl
 
      if(dreal(fcm/fcl) .lt. zero) fcm=-fcm
      gam=(f-p)/q
      ta=absc(gam+pm)
      paccq=eps*max(ta,one/ta)
      hcl=fcm*(gam+pm)*(sfsh/(sf*sf))
 
c                            consider a kase = 1 calculation
 
      if(paccq .gt. accur .and. kase .gt. 0) then
       f11v=c309r2(x,eta,z11,p11,acct,limit,0,err,n11,fpmax,acc8,acc16)
       if(err .lt. paccq) go to 200
      end if
      rerr=max(rerr,paccq)
      go to 230
c
c *** arrive here if kase = 4
c     to evaluate the exponentially decreasing h(lh) directly.
c
  170 if(donem) go to 180
      aa=etap-zlm
      bb=etap+zlm+one
      f20v=c309r3(aa,bb,-half*pm*xi,acct,jmax,err,fpmax,n20,xrcf)
      if(n20 .le. 0) go to 190
      rerr=max(rerr,err)
      hcl=fpmin
      if(abs(dreal(pm*thetam)+omega*scale) .gt. fplmax) go to 330
  180 hcl=f20v*exp(pm*thetam+omega*scale)
      fcm=sfsh/((f-pq1)*hcl)
      go to 230
c
c *** arrive here if kase=1   (or if 2f0 tried mistakenly & failed)
c
c           for small values of x, calculate f(x,sl) directly from 1f1
c               using dreal*16 arithmetic if possible.
c           where z11 = zll if id>0, or = zlm if id<0
c
  190 f11v=c309r2(x,eta,z11,p11,acct,limit,0,err,n11,fpmax,acc8,acc16)
  200 if(n11 .lt. 0) then
 
c                               f11 failed from bb = negative integer
 
       if(lpr) write(6,1060) '-l',one
       ifail=-1
       go to 290
      end if
 
c                      consider a kase 2 or 3 calculation :
 
      if(err .gt. paccq .and. paccq .lt. accb) then
       kase=-2
       if(axial) kase=-3
       go to 130
      end if
      rerr=max(rerr,err)
      if(err .gt. fpmax) go to 370
      if(id .lt. 0) cll=z11*tlog-hpi*eta-wlogam(bb)
     1                        +wlogam(z11+one+p11*eta)-p11*sigma
      ek=(z11+one)*xlog-p11*x+cll-abs(scale)
      if(id .gt. 0) ek=ek-fesl+log(fcl)
      if(dreal(ek) .gt. fplmax) go to 350
      if(dreal(ek) .lt. fplmin) go to 340
      fcm=f11v*exp(ek)
      if(kase .ge. 5) then
       if(absc(zlm+zlm-nintc(zlm+zlm)) .lt. acch) kase=6
c
c ***  for abs(x) < xnear, then cf2 may not converge accurately, so
c ***      use an expansion for irregular soln from origin :
c
       sl=zlm
       zlneg=dreal(zlm) .lt. -one+accb
       if(kase .eq. 5 .or. zlneg) sl=-zlm-one
       pk=sl+one
       aa=pk-etap
       ab=pk+etap
       bb=two*pk
       clgaa=wlogam(aa)
       clgab=clgaa
       if(etane0) clgab=wlogam(ab)
       clgbb=wlogam(bb)
       if(kase .eq. 6 .and. .not.zlneg) then
        if(npint(aa,accur)) clgaa=clgab-two*pm*sigma
        if(npint(ab,accur)) clgab=clgaa+two*pm*sigma
       end if
       cll=sl*tlog-hpi*eta-clgbb+(clgaa+clgab)*half
       dsig=(clgaa-clgab)*pm*half
       if(kase .eq. 6) p11=-pm
       ek=pk*xlog-p11*x+cll-abs(scale)
       sf=exp(-abs(scale))
       chi=zero
       if(.not.(kase .eq. 5 .or. zlneg)) go to 210
c
c ***  use  g(l)  =  (cos(chi) * f(l) - f(-l-1)) /  sin(chi)
c
c      where chi = sig(l) - sig(-l-1) - (2l+1)*pi/2
c
       chi=sigma-dsig-(zlm-sl)*hpi
       f11v=c309r2(x,eta,sl,p11,acct,limit,0,err,npq(1),
     1             fpmax,acc8,acc16)
       rerr=max(rerr,err)
       if(kase .eq. 6) go to 210
       fesl=f11v*exp(ek)
       fcl1=exp(pm*chi)*fcm
       hcl=fcl1-fesl
       rerr=max(rerr,acct*max(absc(fcl1),absc(fesl))/absc(hcl))
       hcl=hcl/sin(chi)*(sfsh/(sf*sf))
       go to 220
c
c *** use the logarithmic expansion for the irregular solution (kase 6)
c        for the case that bb is integral so sin(chi) would be zero.
c
  210  rl=bb-one
       n=nintc(rl)
       zlog=xlog+tlog-pm*hpi
       chi=chi+pm*thetam+omega*scale+ab*zlog
       aa=one-aa
       if(npint(aa,accur)) then
        hcl=zero
       else
        if(id .gt. 0 .and. .not.zlneg) f11v=fcm*exp(-ek)
        hcl=exp(chi-clgbb-wlogam(aa))*(-1)**(n+1)*(f11v*zlog+
     1   c309r2(x,eta,sl,-pm,acct,limit,2,err,npq(2),fpmax,acc8,acc16))
        rerr=max(rerr,err)
       end if
       if(n .gt. 0) then
        ek=chi+wlogam(rl)-clgab-rl*zlog
        df=c309r2(x,eta,-sl-one,-pm,zero,n,0,err,l,fpmax,acc8,acc16)
        hcl=hcl+exp(ek)*df
       end if
       rerr=max(rerr,two*abs(bb-nintc(bb)))
  220  pq1=f-sfsh/(fcm*hcl)
      else
       if(mode .le. 2) hcl=sfsh/((f-pq1)*fcm)
       kase=1
      end if
c
c ***  now have absolute normalisations for coulomb functions
c          fcm & hcl  at lambda = zlm
c      so determine linear transformations for functions required :
c
  230 ih=abs(mode1)/10
      if(kfn .eq. 3) ih=(3-dimag(cik))/2+half
      p11=one
      if(ih .eq. 1) p11=ci
      if(ih .eq. 2) p11=-ci
      df=-pm
      if(ih .ge. 1) df=-pm+p11
      if(absc(df) .lt. acch) df=zero
c
c *** normalisations for spherical or cylindrical bessel functions
c
      if(kfn .le. 0) then
       alfa=zero
       beta=one
      else if(kfn .eq. 1) then
       alfa=xi
       beta=xi
      else
       alfa=xi*half
       beta=sqrt(xi/hpi)
       if(dreal(beta) .lt. zero) beta=-beta
      end if
      aa=one
      if(kfn .gt. 0) aa=-p11*beta
 
c                calculate rescaling factors for i & k output
 
      if(kfn .ge. 3) then
       p=exp((zlm+dell)*hpi*cik)
       aa=beta*hpi*p
       beta=beta/p
       q=cik*id
      end if
 
c                  calculate rescaling factors for gc output
 
      if(ih .eq. 0) then
       ta=abs(scale)+dimag(pm)*scale
       rk=zero
       if(ta .lt. fplmax) rk=exp(-ta)
      else
       ta=abs(scale)+dimag(p11)*scale
       if(absc(df) .gt. acch .and. ta .gt. fplmax) go to 320
       if(absc(df) .gt. acch) df=df*exp(ta)
       sf=two*(lh-ih)*scale
       rk=zero
       if(sf .gt. fplmax) go to 320
       if(sf .gt. fplmin) rk=exp(sf)
      end if
      kas((3-id)/2)=kase
      w=fcm/fcl
      if(log(absc(w))+log(absc(fc(lf))) .lt. fplmin) go to 340
      if(mode .ge. 3) go to 240
      if(lpr .and. absc(f-pq1) .lt. acch*absc(f))
     1                             write(6,1020) lh,zlm+dell
      hpl=hcl*pq1
      if(absc(hpl) .lt. fpmin .or. absc(hcl) .lt. fpmin) go to 330
c
c *** idward recurrence from hcl,hpl(lf) (stored gc(l) is rl if reqd)
c *** renormalise fc,fcp at each lambda
c ***    zl   = zlm - min(id,0) here
c
  240 do 270 l = lf,l1,id
      fcl=w*fc(l)
      if(absc(fcl) .lt. fpmin) go to 340
      if(ifcp) fpl=w*fcp(l)
      fc(l)=beta*fcl
      if(ifcp) fcp(l)=beta*(fpl-alfa*fcl)*cik
      fc(l)=c309r8(fc(l),accur)
      if(ifcp) fcp(l)=c309r8(fcp(l),accur)
      if(mode .ge. 3) go to 260
      if(l .eq. lf) go to 250
      zl=zl+id
      zid=zl*id
      rl=gc(l)
      if(etane0) then
       sl=eta+zl*zl*xi
       if(mode .eq. 1) then
        pl=gcp(l)
       else
        pl=zero
        if(absc(zl) .gt. acch) pl=(sl*sl-rl*rl)/zid
       end if
       hcl1=(sl*hcl-zid*hpl)/rl
       hpl=(sl*hpl-pl*hcl)/rl
      else
       hcl1=rl*hcl-hpl*id
       hpl=(hcl-rl*hcl1)*id
      end if
      hcl=hcl1
      if(absc(hcl) .gt. fpmax) go to 320
  250 gc(l)=aa*(rk*hcl+df*fcl)
      if(mode .eq. 1) gcp(l)=(aa*(rk*hpl+df*fpl)-alfa*gc(l))*cik
      gc(l)=c309r8(gc(l),accur)
      if(mode .eq. 1) gcp(l)=c309r8(gcp(l),accur)
      if(kfn .ge. 3) aa=aa*q
  260 if(kfn .ge. 3) beta=-beta*q
  270 last=min(last,(l1-l)*id)
      go to 280
c
c *** come here after all soft errors to determine how many l values ok
c
  310 if(lpr) write(6,1000) zz
      go to 999
  320 if(lpr) write(6,1010) zl+dell,'ir',hcl,'>',fpmax
      go to 280
  330 if(lpr) write(6,1010) zl+dell,'ir',hcl,'<',fpmin
      go to 280
  340 if(lpr) write(6,1010) zl+dell,'  ',fcl,'<',fpmin
      go to 280
  350 if(lpr) write(6,1010) zl+dell,'  ',fcl,'>',fpmax
      go to 280
  360 if(lpr) write(6,1030) zl+dell
      go to 280
  370 if(lpr) write(6,1040) z11,i
      ifail=-1
      go to 290
  380 if(lpr) write(6,1050) zlmin,zlm,zlm+one,zlmin+nl
      ifail=-1
      go to 290
  280 if(id .gt. 0 .or. last .eq. 0) ifail=last
      if(id .lt. 0 .and. last .ne. 0) ifail=-3
c
c *** come here after all errors for this l range (zlm,zll)
c
c *** so on first block, 'f' started decreasing monotonically,
c                        or hit bound states for low zl.
c     thus redo m1 to lf-1 in reverse direction, i.e. do
c      cf1a at zlmin & cf2 at zlm (midway between zlmin & zlmax)
c
  290 if(id .gt. 0 .and. lf .ne. m1) then
       id=-1
       if(.not.unstab) lf=lf-1
       donem=unstab
       lf=min(lf,l1)
       l1=m1
       go to 10
      end if
      if(ifail .lt. 0) go to 999
      if(lpr .and. rerr .gt. accb) write(6,1070) rerr
      if(rerr .gt. 0.1d0) ifail=-4
  999 do 998 l = 1,nl+1
      fc(l-1)=fc(l)
      gc(l-1)=gc(l)
      fcp(l-1)=fcp(l)
      gcp(l-1)=gcp(l)
  998 sig(l-1)=sig(l)
      return
c
 1000 format(1x,'***** cern c309 wclbes ... ',
     1 'cannot calculate irregular solutions for x =',
     2 1p,2d10.2,' abs(x) too small')
 1010 format(1x,'***** cern c309 wclbes ... ',
     1 'at zl =',2f8.3,' ',a2,'regular solution (',1p,2e10.1,') ',
     2 a1,e10.1)
 1020 format(1x,'***** cern c309 wclbes ... ',
     1 'warning: linear independence between ''f'' and ''h(',i1,
     2 ')'' is lost at zl = ',2f7.2/1x,'*****',22x,'(eg. coulomb ',
     3 'eigenstate or cf1 unstable)')
 1030 format(1x,'***** cern c309 wclbes ... ',
     2 '(eta & l)/x too large for cf1a, and cf1 unstable at l = ',2f8.2)
 1040 format(1x,'***** cern c309 wclbes ... ',
     1 'overflow in 1f1 series at zl = ',2f8.3,' at term',i5)
 1050 format(1x,'***** cern c309 wclbes ... ',
     1 'both bound-state poles and f-instabilities occur or multiple',
     2 ' instabilities present'/
     3   1x,'*****',22x,'try calling twice, 1st for zl from',2f8.3,
     4 ' to',2f8.3,' (incl)'/1x,'*****',41x,'2nd for zl from',2f8.3,
     5 ' to',2f8.3)
 1060 format(1x,'***** cern c309 wclbes ... ',
     1 'warning: as ''',a2,''' reflection rules not used ',
     2 'errors can be up to',1pd12.2)
 1070 format(1x,'***** cern c309 wclbes ... ',
     1 'warning: overall roundoff error approximately',1pe11.1)
      end

cdeck  id>, wdigam.
      function wdigam(u)
c*ns  complex*16 wdigam,w,u,v,h,r,p
c*ns  double precision c(6),x,y,a,pi,zero,one,half,c1,d
      complex*16 wdigam,  u,v,h,r,p
      double precision c(6),x,y,a,pi,zero,one,half,c1
 
      data pi /3.14159 26535 89793 24d0/
      data zero /0d0/, one /1d0/, half /0.5d0/, c1 /15d0/
 
      data c(1) / 8.33333 33333 33333 33d-2/
      data c(2) /-8.33333 33333 33333 33d-3/
      data c(3) / 3.96825 39682 53968 25d-3/
      data c(4) /-4.16666 66666 66666 67d-3/
      data c(5) / 7.57575 75757 57575 76d-3/
      data c(6) /-2.10927 96092 79609 28d-2/
 
c*ul    9 x=dreal(u)
      x=dreal(u)
      a=abs(x)
      if(dimag(u) .eq. zero .and. -a .eq. int(x)) then
       wdigam=zero
       write(6,100) x
       return
      end if
      if(x .lt. zero) u=-u
      v=u
      h=zero
      if(a .lt. c1) then
       h=one/v
       do 1 i = 1,14-int(a)
       v=v+one
    1  h=h+one/v
       v=v+one
      end if
      r=one/v**2
      p=r*c(1)
      do 2 i = 6,1,-1
    2 p=r*(c(i)+p)
      h=log(v)-half/v-p-h
      if(x .lt. zero) then
       v=pi*u
       x=dreal(v)
       y=tanh(dimag(v))
       a=sin(x)
       x=cos(x)
       h=h+one/u+pi*dcmplx(x,-a*y)/dcmplx(a,x*y)
      end if
      wdigam=h
 
  100 format(7x,'***** cern c307 wdigam ... argument equals ',
     1          'non-positive integer = ',e10.1)
      end
      complex function cdigam(z)
      complex z
      complex*16 wdigam,w,ww
      double precision d
 
      sround(d) = sngl(d+(d-dble(sngl(d))))
c*ia  d
      data d/0/
 
      w  = z
      ww = wdigam(w)
      cdigam = cmplx(sround(dreal(ww)),sround(dimag(ww)))
 
      end
cdeck  id>, wlogam.
      function wlogam(v)
c*ns  complex*16 wlogam,w,u,v,h,p,r
c*ns  double precision c(10),pi,t,x,y,ya,a,ur,ui,c1,c2,d
      complex*16 wlogam,  u,v,h,p,r
      double precision c(10),pi,t,x,y,ya,a,ur,ui,c1,c2
      double precision zero,one,seven,half
 
      data zero /0d0/, one /1d0/, seven /7d0/, half /0.5d0/
 
      data pi /3.14159 26535 89793 24d+0/
      data c1 /9.18938 53320 46727 42d-1/
      data c2 /1.14472 98858 49400 17d+0/
 
      data c( 1) / 8.33333 33333 33333 33d-2/
      data c( 2) /-2.77777 77777 77777 78d-3/
      data c( 3) / 7.93650 79365 07936 51d-4/
      data c( 4) /-5.95238 09523 80952 38d-4/
      data c( 5) / 8.41750 84175 08417 51d-4/
      data c( 6) /-1.91752 69175 26917 53d-3/
      data c( 7) / 6.41025 64102 56410 26d-3/
      data c( 8) /-2.95506 53594 77124 18d-2/
      data c( 9) / 1.79644 37236 88305 73d-1/
      data c(10) /-1.39243 22169 05901 12d+0/
 
c*ul9 x=dreal(v)
      x=dreal(v)
      y=dimag(v)
      if(y .eq. zero .and. -abs(x) .eq. int(x)) then
       wlogam=zero
       write(6,100) x
       return
      end if
      ya=abs(y)
      u=dcmplx(x,ya)
      if(x .lt. zero) u=one-u
      h=zero
      ur=dreal(u)
      if(ur .lt. seven) then
       ui=dimag(u)
       a=atan2(ui,ur)
       h=u
       do 1 i = 1,6-int(ur)
       ur=ur+one
       u=dcmplx(ur,ui)
       h=h*u
    1  a=a+atan2(ui,ur)
       h=dcmplx(half*log(dreal(h)**2+dimag(h)**2),a)
       u=u+one
      end if
      r=one/u**2
      p=r*c(10)
      do 2 i = 9,2,-1
    2 p=r*(c(i)+p)
      h=c1+(u-half)*log(u)-u+(c(1)+p)/u-h
      if(x .lt. zero) then
       ur=int(x)-one
       ui=pi*(x-ur)
       x=pi*ya
       t=exp(-x-x)
       a=sin(ui)
       t=x+half*log(t*a**2+(half*(one-t))**2)
       a=atan2(cos(ui)*tanh(x),a)-ur*pi
       h=c2-dcmplx(t,a)-h
      end if
      if(y .lt. zero) h=dconjg(h)
      wlogam=h
 
  100 format(7x,'***** cern c333 wlogam ... argument equals ',
     1          'non-positive integer = ',e10.1)
      end
      complex function clogam(z)
      complex z
      complex*16 wlogam,w,ww
      double precision d
 
      sround(d) = sngl(d+(d-dble(sngl(d))))
c*ia  d
      data d/0/
 
      w  = z
      ww = wlogam(w)
      clogam = cmplx(sround(dreal(ww)),sround(dimag(ww)))
 
      end
cdeck  id>, r1dp.
      function c309r1(x,eta,zl,pm,eps,limit,err,npq,acc8,acch,
     1             lpr,accur,dell)
c
c                                    (omega)        (omega)
c *** evaluate  cf2  = p + pm.q  =  h   (eta,x)' / h   (eta,x)
c                                    zl             zl
c     where pm = omega.i
c
      implicit complex*16(a-h,o-z)
      logical lpr
      double precision eps,err,acc8,acch,accur,ta,rk
      double precision absc,half
 
      parameter(half = 1d0/2d0)
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
 
      ta=limit+limit
      etap=eta*pm
      xi=1/x
      wi=etap+etap
      rk=0
      pq=(1-eta*xi)*pm
      aa=-(eta*eta+zl*zl+zl)+etap
      bb=2*(x-eta+pm)
      rl=xi*pm
      if(absc(bb) .lt. acch) then
       rl=rl*aa/(aa+rk+wi)
       pq=pq+rl*(bb+pm+pm)
       aa=aa+2*(rk+1+wi)
       bb=bb+4*pm
       rk=rk+4
      end if
      dd=1/bb
      dl=aa*dd*rl
   10 pq=pq+dl
      rk=rk+2
      aa=aa+rk+wi
      bb=bb+pm+pm
      dd=1/(aa*dd+bb)
      dl=dl*(bb*dd-1)
      err=absc(dl)/absc(pq)
      if(err .ge. max(eps,acc8*rk*half) .and. rk .le. ta) go to 10
c
      npq=half*rk
      c309r1=pq+dl
      if(lpr .and. npq .ge. limit-1 .and. err .gt. accur)
     1 write(6,1000) int(dimag(pm)),npq,err,zl+dell
      return
 1000 format(1x,'***** cern c309 wclbes ... ',
     2 'cf2(',i2,') not converged fully in ',i7,' iterations'/1x,27x,
     3 'error in irregular solution =',1p,d11.2,' at zl = ',2f8.3)
      end
cdeck  id>, r2dp.
      function c309r2(x,eta,zl,p,eps,limit,kind,err,nits,
     1                fpmax,acc8,acc16)
c
c *** evaluate the hypergeometric function
c                                        i
c            f (aa;bb; z) = sum  (aa)   z / ( (bb)  i! )
c           1 1              i       i            i
c
c     to accuracy eps with at most limit terms.
c  if kind = 0 : using extended precision but real arithmetic only,
c            1 : using normal precision in complex arithmetic,
c   or       2 : using normal complex arithmetic, but with wdigam factor
c
c  where aa, bb, and z are defined below
c
      implicit double precision(a-h,o-z)
      complex*16 x,eta,zl,p,aa,bb,z,c309r2,wdigam
      complex*16 dd,g,f,ai,bi,t
      double precision ar,br,gr,gi,dr,di,tr,ti,ur,ui,fi,fi1,den
 
      parameter(tbbb = 3d0/2d0)
 
      absc(aa)=abs(dreal(aa))+abs(dimag(aa))
      nintc(aa)=nint(dreal(aa))
 
      aa=zl+1-eta*p
      bb=2*zl+2
      z=2*p*x
      if(dreal(bb) .le. 0 .and. abs(bb-nintc(bb)) .lt.
     1 sqrt(sqrt(acc8))**3 .and. dreal(bb)+limit .ge. tbbb) then
       nits=-1
       return
      end if
      if(limit .le. 0) then
       c309r2=0
       err=0
       nits=1
       return
      end if
      ta=1
      rk=1
      if(kind .le. 0 .and. absc(z)*absc(aa) .gt. absc(bb)*1.0) then
       dr=1
       di=0
       gr=1
       gi=0
       ar=dreal(aa)
       br=dreal(bb)
       fi=0
       do 20 i = 2,limit
       fi1=fi+1
       tr=br*fi1
       ti=dimag(bb)*fi1
       den=1/(tr*tr+ti*ti)
       ur=(ar*tr+dimag(aa)*ti)*den
       ui=(dimag(aa)*tr-ar*ti)*den
       tr=ur*gr-ui*gi
       ti=ur*gi+ui*gr
       gr=dreal(z)*tr-dimag(z)*ti
       gi=dreal(z)*ti+dimag(z)*tr
       dr=dr+gr
       di=di+gi
       err=abs(gr)+abs(gi)
       if(err .gt. fpmax) go to 60
       rk=abs(dr)+abs(di)
       ta=max(ta,rk)
       if(err .lt. rk*eps .or. i .ge. 4 .and. err .lt. acc16) go to 30
       fi=fi1
       ar=ar+1
   20  br=br+1
c
   30  c309r2=dr+(0,1)*di
       err=acc16*ta/rk
      else
c
c*    if real*16 arithmetic is not available, (or already using it!),
c*    then use kind > 0
c
       g=1
       f=1
       if(kind .ge. 2) f=wdigam(aa)-wdigam(bb)-wdigam(g)
       dd=f
       do 40 i = 2,limit
       ai=aa+(i-2)
       bi=bb+(i-2)
       r=i-1
       g=g*z*ai/(bi*r)
c
c                       multiply by (psi(a+r)-psi(b+r)-psi(1+r))
c
       if(kind .eq. 2) f=f+1/ai-1/bi-1/r
       t=g*f
       dd=dd+t
       err=absc(t)
       if(err .gt. fpmax) go to 60
       rk=absc(dd)
       ta=max(ta,rk)
       if(err .lt. rk*eps .or. err .lt. acc8 .and. i .ge. 4) go to 50
   40  continue
 
   50  err=acc8*ta/rk
       c309r2=dd
      end if
   60 nits=i
      return
      end
cdeck  id>, r3dp.
      function c309r3(aa,bb,z,eps,jmax,re,fpmax,n,x)
c
c     evaluate the hypergeometric function
c                                             i
c            f (aa,bb;;z) = sum  (aa)  (bb)  z / i!
c           2 0              i       i     i
c
c     to accuracy eps with at most jmax terms.
c
c     if the terms start diverging,
c     the corresponding continued fraction is found by rcf
c     & evaluated progressively by steed's method to obtain convergence.
c
      implicit complex*16(a-h,o-z)
      dimension x(jmax,4)
      logical finite
      double precision ep,eps,at,atl,absc,re,fpmax
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
      nintc(w)=nint(dreal(w))
c
      re=0
      x(1,1)=1
      sum=x(1,1)
      atl=absc(x(1,1))
      f=sum
      d=1
      df=sum
      j=0
      ep=eps*(10*jmax)
      ma=-nintc(aa)
      mb=-nintc(bb)
      finite=abs(abs(dreal(aa))-ma) .lt. ep .and. abs(dimag(aa)) .lt. ep
     1  .or. abs(abs(dreal(bb))-mb) .lt. ep .and. abs(dimag(bb)) .lt. ep
      imax=jmax
      if(finite .and. ma .ge. 0) imax=min(ma+1,imax)
      if(finite .and. mb .ge. 0) imax=min(mb+1,imax)
      do 10 i = 2,imax
      x(i,1)=x(i-1,1)*z*(aa+i-2)*(bb+i-2)/(i-1)
      if(absc(x(i,1)) .gt. fpmax) then
       n=0
       c309r3=sum
       if(.not.finite) re=at/absc(sum)
       return
      end if
      at=absc(x(i,1))
      if(j .eq. 0) then
       sum=sum+x(i,1)
       if(at .lt. absc(sum)*eps) then
        n=i
        c309r3=sum
        if(.not.finite) re=at/absc(sum)
        return
       end if
      end if
      if(finite) go to 10
      if(j .gt. 0 .or. at .gt. atl .or. i .ge. jmax-2) j=j+1
      if(j .eq. 0) go to 10
      call c309r7(x(1,1),x(1,2),j,i,x(1,3),eps)
      if(i .lt. 0) then
       n=0
       c309r3=sum
       if(.not.finite) re=at/absc(sum)
       return
      end if
      do 50 k = max(j,2),i
      d=1/(d*x(k,2)+1)
      df=df*d-df
      f=f+df
      if(absc(df) .lt. absc(f)*eps .or.
     1   df .eq. 0 .and. f .eq. 0 .and. i .ge. 4) then
       n=k
       c309r3=f
       re=absc(df)/absc(f)
       return
      end if
   50 continue
      j=i
   10 atl=at
      if(.not.finite) i=-jmax
      n=i
      c309r3=sum
      if(.not.finite) re=at/absc(sum)
      return
      end
cdeck  id>, r4dp.
      function c309r4(x,eta,zl,eps,fcl,tpk1,etane0,limit,err,nfp,
     1              fpmin,fpmax,lpr)
c
c ***   evaluate cf1 = f'(zl,eta,x)/f(zl,eta,x)    (real)
c
      implicit double precision(a-h,o-z)
      logical lpr,etane0
 
      fcl=1
      xi=1/x
      pk=zl+1
      px=pk+limit
      ek=eta/pk
      f=ek+pk*xi
      if(abs(f) .lt. fpmin) f=fpmin
      d=0
      c=f
      small=sqrt(fpmin)
      rk2=1+ek*ek
c
c ***   begin cf1 loop on pk = k = lambda + 1
c
   10 pk1=pk+1
      tpk1=pk+pk1
      if(etane0) then
       ek=eta/pk
       rk2=1+ek*ek
       tk=tpk1*(xi+ek/pk1)
      else
       tk=tpk1*xi
      end if
      c=tk-rk2/c
      d=tk-rk2*d
      if(abs(c) .lt. fpmin) c=fpmin
      if(abs(d) .lt. fpmin) d=fpmin
      d=1/d
      df=d*c
      f=f*df
      fcl=fcl*d*tpk1*xi
      if(abs(fcl) .lt. small) fcl=fcl/small
      if(abs(fcl) .gt. fpmax) fcl=fcl*fpmin
      pk=pk1
      if(pk .le. px) then
       if(abs(df-1) .ge. eps) go to 10
       nfp=pk-zl-1
       err=eps*sqrt(real(nfp))
       c309r4=f
      else
       if(lpr) write (6,1000) limit,abs(x)
       err=2
      end if
      return
 1000 format(1x,'***** cern c309 wclbes ... cf1 (real) has failed ',
     1'to converge after',i10,' iterations as abs(x) =',f15.0)
      end
cdeck  id>, r5dp.
      function c309r5(x,eta,zl,eps,fcl,tpk1,etane0,limit,err,nfp,
     1              fpmin,fpmax,lpr)
c
c ***    evaluate cf1  =  f'(zl,eta,x)/f(zl,eta,x)  (complex)
c
      implicit complex*16(a-h,o-z)
      logical lpr,etane0
      double precision eps,err,fpmin,fpmax,absc,small,px
 
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
 
      fcl=1
      xi=1/x
      pk=zl+1
      px=pk+limit
      ek=eta/pk
      f=ek+pk*xi
      if(absc(f) .lt. fpmin) f=fpmin
      d=0
      c=f
      small=sqrt(fpmin)
      rk2=1+ek*ek
c
c ***   begin cf1 loop on pk = k = lambda + 1
c
   10 pk1=pk+1
      tpk1=pk+pk1
      if(etane0) then
       ek=eta/pk
       rk2=1+ek*ek
       tk=tpk1*(xi+ek/pk1)
      else
       tk=tpk1*xi
      end if
      c=tk-rk2/c
      d=tk-rk2*d
      if(absc(c) .lt. fpmin) c=fpmin
      if(absc(d) .lt. fpmin) d=fpmin
      d=1/d
      df=d*c
      f=f*df
      fcl=fcl*d*tpk1*xi
      if(absc(fcl) .lt. small) fcl=fcl/small
      if(absc(fcl) .gt. fpmax) fcl=fcl*fpmin
      pk=pk1
      if(dreal(pk) .le. px) then
       if(absc(df-1) .ge. eps) go to 10
       nfp=pk-zl-1
       err=eps*sqrt(real(nfp))
       c309r5=f
      else
       if(lpr) write (6,1000) limit,abs(x)
       err=2
      end if
      return
 1000 format(1x,'***** cern c309 wclbes ... cf1 (complex) has failed ',
     1'to converge after',i10,' iterations as abs(x) =',f15.0)
      end
cdeck  id>, r6dp.
      function c309r6(rho,eta,xl,psi,eps,nmax,nused,fcl,re,fpmax,xx,g,c)
c
c     evaluate the asymptotic expansion for the
c            logarithmic derivative of the regular solution
c
c ***        cf1a  =  f   =  f'(xl,eta,rho)/f(xl,eta,rho)
c
c      that is valid for dreal(rho)>0, and best for rho >> eta**2, xl,
c      and is derived from the 2f0 expansions for h+ and h-
c      e.g. by froeberg (rev. mod. physics vol 27, p399 , 1955)
c      some lines of this subprogram are for convenience copied from
c           takemasa, tamura & wolter cpc 17 (1979) 351.
c
c     evaluate to accuracy eps with at most nmax terms.
c
c     if the terms start diverging,
c     the corresponding continued fraction is found by rcf
c     & evaluated progressively by steed's method to obtain convergence.
c
      implicit complex*16(a-h,o-z)
      dimension xx(2,nmax),g(nmax),c(nmax)
      double precision re,eps,t1,t2,t3,at,atl,absc,fpmax
      double precision hpi
 
      parameter(hpi  = 1.57079 63267 94896 619d0)
 
      absc(w)=abs(dreal(w))+abs(dimag(w))
c
      t1=sin(dreal(psi))
      t2=cos(dreal(psi))
      atl=tanh(dimag(psi))
 
c         give cos(psi)/cosh(im(psi)), which always has correct sign
 
      cosl=dcmplx(t2,-t1*atl)
      tanl=dcmplx(t1,t2*atl)/cosl
      re=0
      xll1=xl*xl+xl
      etasq=eta*eta
      sl1=1
      sl=sl1
      sc1=0
      sc=sc1
      tl1=sc
      tl=tl1
      tc1=1-eta/rho
      tc=tc1
      fcl=tl+sl*tanl
      g(1)=(tc+sc*tanl)/fcl
      glast=g(1)
      atl=absc(glast)
      f=glast
      d=1
      df=glast
      j=0
      do 10 n = 2,nmax
      t1=n-1
      t2=2*t1-1
      t3=t1*t1-t1
      denom=2*rho*t1
      c1=(eta*t2)/denom
      c2=(etasq+xll1-t3)/denom
      sl2=c1*sl1-c2*tl1
      tl2=c1*tl1+c2*sl1
      sc2=c1*sc1-c2*tc1-sl2/rho
      tc2=c1*tc1+c2*sc1-tl2/rho
      sl=sl+sl2
      tl=tl+tl2
      sc=sc+sc2
      tc=tc+tc2
      sl1=sl2
      tl1=tl2
      sc1=sc2
      tc1=tc2
      fcl=tl+sl*tanl
      if(absc(fcl) .gt. fpmax .or. absc(fcl) .lt. 1./fpmax) then
       c309r6=g(1)
       fcl=1
       re=1
       nused=0
       return
      end if
      gsum=(tc+sc*tanl)/fcl
      g(n)=gsum-glast
      glast=gsum
      at=absc(g(n))
      if(at .lt. absc(gsum)*eps) then
       fcl=fcl*cosl
       c309r6=gsum
       re=at/absc(gsum)
       nused=n
       return
      end if
      if(j .gt. 0 .or. at .gt. atl .or. n .ge. nmax-2) j=j+1
      if(j .eq. 0) go to 10
      call c309r7(g,c,j,n,xx,eps)
      if(n .lt. 0) then
       c309r6=g(1)
       fcl=1
       re=1
       nused=0
       return
      end if
      do 60 k = max(j,2),n
      d=1/(d*c(k)+1)
      df=df*d-df
      f=f+df
      if(absc(df) .lt. absc(f)*eps .or.
     1   df .eq. 0 .and. f .eq. 0 .and. n .ge. 4) then
       c309r6=f
       fcl=fcl*cosl
       re=absc(df)/absc(f)
       nused=k
       return
      end if
   60 continue
      j=n
   10 atl=at
      c309r6=f
      fcl=fcl*cosl
      re=absc(df)/absc(f)
      nused=-nmax
      return
      end
cdeck  id>, r7dp.
      subroutine c309r7(a,b,ibeg,inum,xx,eps)
c
c*******************************************************************
c
c  rcf converts polynomial a to the corresponding continued
c         fraction, in 'normal'  form with coefficients b
c         by the 'p algorithm' of patry & gupta
c
c   a(z) = a1/z + a2/z**3 + a3/z**5 + ... + an/z**(2n-1)
c
c   b(z) = b1/z+ b2/z+ b3/z+ .../(z+ bn/z)
c
c  data:
c   a     vector a(k), k=1,inum         input
c   b     vector b(k), k=ibeg,inum      output
c   ibeg  order of first coef. calc.    input
c   inum  order of a, even or odd       input
c   xx    auxiliary vector of length .ge. length of vector b
c         caller provides space for a,b,xx
c   note that neither of the first two terms a(1) a(2) should be zero
c          & the user can start the calculation with any value of
c          ibeg provided the c.f. coefs have been already
c          calculated up to inum = ibeg-1
c          & the method breaks down as soon as the absolute value
c          of a c.f. coef. is less than eps.    at the time of the
c          break up  inum has been replaced by minus times the number
c          of this coefficient.
c   algorithm: j. patry & s. gupta, eir-bericht 247, november 1973
c              eidg. institut fur reaktorforschung
c              wuerenlingen, switzerland
c   see also:  haenggi, roesel & trautmann,
c              j. comput. phys., v. 137, (1980) 242-258
c   note:      restart procedure modified by i.j.thompson
c
c*******************************************************************
c
      implicit complex*16(a-h,o-z)
      dimension a(100),b(100),xx(2,100)
      logical even
      double precision eps
 
      ibn=inum
      if(ibeg .gt. 4) go to 50
      if(ibeg .eq. 4) go to 20
      b(1)=a(1)
      if(ibn .ge. 2) b(2)=-a(2)/a(1)
      if(ibn .lt. 3) return
      x0=a(3)/a(2)
      xx(2,1)=b(2)
      xx(1,1)=-x0
      xx(1,2)=0
      b(3)=-x0-b(2)
      x0=-b(3)*a(2)
      m=3
      mp12=2
      even=.true.
      if(ibn .le. 3) return
   20 if(abs(b(3)) .lt. eps*abs(x0)) then
       inum=-m
       return
      end if
      m=4
   30 x1=a(m)
      m2m1=mp12
      mp12=m2m1+1
      if(even) mp12=m2m1
      do 40 k = 2,mp12
   40 x1=x1+a(m-k+1)*xx(1,k-1)
      b(m)=-x1/x0
      if(m .ge. ibn) return
   50 if(abs(b(m)) .lt. eps*abs(x0)) then
       inum=-m
       return
      end if
      do 60 k = m2m1,2,-1
   60 xx(2,k)=xx(1,k)+b(m)*xx(2,k-1)
      xx(2,1)=xx(1,1)+b(m)
      do 70 k = 1,m2m1
      x0=xx(2,k)
      xx(2,k)=xx(1,k)
   70 xx(1,k)=x0
      x0=x1
      xx(1,m2m1+1)=0
      m=m+1
      even=.not.even
      go to 30
      end
cdeck  id>, r8dp.
      function c309r8(z,acc)
      complex*16 c309r8,z
      double precision acc
      double precision x,y,ax,ay,a
 
      x=dreal(z)
      y=dimag(z)
      ax=abs(x)
      ay=abs(y)
      a=5*acc*(ax+ay)
      if(ax .lt. a) x=0
      if(ay .lt. a) y=0
      c309r8=dcmplx(x,y)
      return
      end
    
    
