program fit
c  Fitting delta in the potential  V(r) = -1/r+lambda(1-exp(-delta r))/r
c     in order to obtain exact energy for scaled two-electron atoms
      implicit real*8 (a-h,o-z)
      dimension tab(2,1000)
      common nq,en,enex,smax,alam,delta,v0,v1
      external func
      delta0=0.125
      dd=0.005
      eps=1.d-8
      eta=0
      open(1,file='fit-inp.dat')
      read(1,*)nel,nq
      ntab=0
   1  if(eof(1))goto2
      ntab=ntab+1
      read(1,*)alam,enex
      tab(1,ntab)=alam
      tab(2,ntab)=enex
      goto1
   2  print*,'Ntab =',ntab      
      open(2,file='fit-out.dat')
c Input file has a table of (lambda Ei) for V = -1/r1-1/r2 + lambda/r12 found from 1/Z-expansion
      alam0=0
      i=0
      do 3 n=ntab,1,-1
      i=i+1
      rz=tab(1,n)
      enz=tab(2,n)
      rlam=1/rz
c     For lithium, berillium, Na, Mg *2 or *3, *10, *11
      alam=rlam *(nel-1)
      enex=enz/rz**2
      dlam1=alam-alam0      
      delta1=delta0
      if(i.le.2)goto4
      delta1=delta1+dd*dlam1/dlam0
      dd=dd/10
   4  ifail=0
      call c05agf(delta1,dd,eps,eta,func,a,b,ifail)
      if(ifail.ne.0)print*,'IFAIL =',ifail
      print10,rlam,delta1
  10  format(' -----',2f13.8)
      write(2,20)rlam,delta1
  20  format(1x,2f13.8)      
      dd=delta1-delta0
      if(abs(dd).lt.1.d-4)dd=1.d-4
      delta0=delta1
      alam0=alam
      dlam0=dlam1
   3  continue      
      end
      
      real*8 function func(delta1)           
      implicit real*8 (a-h,o-z)
      dimension xpoint(4),hmax(2,4),y(9),w(9,9)
      common nq,en,enex,smax,alam,delta,v0,v1
      external coeffn,bdyval,monit,report,fcn,fcn1,gcond
      pi=4*datan(1.d0)
      toldefa=1.d-9
c      nq=1
      delta=delta1
      tol0=toldefa
      tolg=tol0*10000
      toll=tol0/4
      x0=tol0**(1.d0/4)/3
  60  format(9x,a19,e10.2)
      ind=nq-1
      en0=-(1-alam)**2/2/nq**2
      if(abs(en0).lt.1.d-3)en0=-1.d-3
      en=en0
      d=(1-alam)**2/2/nq**2
  30  format(' lambda, q.number, Ecl, dE',f6.3,i4,2f20.12)
c   Xend is determined from eq. exp(-q.cl. int.)=sqrt(tolg)
      smax=-dlog(tolg)/2
      x1=x0
      x2=1.d9
      if(v(x1).gt.en .or. v(x2).lt.en)print*,'Root was not localised!'
      do 101 iter=1,100
      x=(x1+x2)/2
      vx=v(x)
      if(vx.lt.en)x1=x
 101  if(vx.ge.en)x2=x
      y(1)=0
      tol=tolg
      ifail=0
      xend=10000
      call d02bhf(x,xend,1,y,tol,0,0.d0,fcn1,gcond,w,ifail)
      xend=x
      xpoint(1)=x0
      xpoint(2)=x0
      xpoint(3)=xend
      xpoint(4)=xend
      match=0
      tol=tolg
      elam=en-d
      delam=2*d
      hmax(1,1)=0
      maxit=100
      maxfun=100 000
      ifail=0
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      en=elam
      d=delam*2
c   Xend is determined from eq. exp(-q.cl. int.)=sqrt(tol0)
      smax=-dlog(tol0)/2
      x1=x0
      x2=1.d9
      if(v(x1).gt.en .or. v(x2).lt.en)print*,'Root was not localised!'
      do 102 iter=1,100
      x=(x1+x2)/2
      vx=v(x)
      if(vx.lt.en)x1=x
 102  if(vx.ge.en)x2=x
      y(1)=0
      tol=tolg
      ifail=0
      xend=10000
      call d02bhf(x,xend,1,y,tol,0,0.d0,fcn1,gcond,w,ifail)
      xend=x
      xpoint(1)=x0
      xpoint(2)=x0
      xpoint(3)=xend
      xpoint(4)=xend
      match=0
      tol=tol0
      elam=en-d
      delam=2*d
      hmax(1,1)=0
      maxit=100
      maxfun=100 000
      ifail=0
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      en=elam
      d=delam*20
c      d=delam*2
c   Repeating the calculation with accelerating accuracy
c   Xend is determined from eq. exp(-q.cl. int.)=sqrt(toll)
      smax=-dlog(toll)/2
      x1=x0
      x2=1.d9
      if(v(x1).gt.en .or. v(x2).lt.en)print*,'Root was not localised!'
      do 103 iter=1,100
      x=(x1+x2)/2
      vx=v(x)
      if(vx.lt.en)x1=x
 103  if(vx.ge.en)x2=x
      y(1)=0
      tol=tolg
      ifail=0
      xend=10000
      call d02bhf(x,xend,1,y,tol,0,0.d0,fcn1,gcond,w,ifail)
      xend=x
      xpoint(1)=x0
      xpoint(2)=x0
      xpoint(3)=xend
      xpoint(4)=xend
      match=0
      tol=toll
      elam=en-d
      delam=2*d
      hmax(1,1)=0
      maxit=100
      maxfun=100 000
      ifail=0
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      print10,'delt,xend,en,den',delta,xend,elam,delam
      en=elam
  40  format(f10.5,f16.8)
  10  format(1x,a15,4f12.6)
      func=enex+en
      return
      end
            subroutine coeffn(p,q,dqdl,x,elam,jint)
      implicit real*8 (a-h,o-z)
      p=-.5d0
      dqdl=-1.d0
      q=v(x)-elam
      return
      end
            subroutine bdyval(xl,xr,elam,yl,yr)
      implicit real*8 (a-h,o-z)
      dimension yl(3),yr(3)
      common nq,en,enex,smax,alam,delta,v0,v1
      h00=1
      h01=v0*h00
      h02=(v0*h01+(v1-en)*h00)/3
      fun=xl*(h00+h01*xl+h02*xl**2)
      dfun=h00+2*h01*xl+3*h02*xl**2
      yl(1)=fun
      yl(2)=dfun/2
      yr(1)=1
      p2=2*(v(xr)-elam)
      dp2=2*vdif(xr)
      y=p2
      if(y.lt.0)y=-y
      p=dsqrt(y)
      yr(2)=-(dp2/p2/4+p)
      return
      end
            subroutine monit(maxit,iflag,elam,finfo)
      implicit real*8 (a-h,o-z)
      dimension finfo(15)
      nfun=finfo(4)+.1
c      print10,maxit,nfun,elam,iflag
  10  format('+',29x,'it,nfun,elam,ifl=',i3,i7,f20.12,i3)
      return
      end
            subroutine report(x,vv,jint)
      implicit real*8 (a-h,o-z)
      dimension vv(3)
      return
      end
            real*8 function v(x)
      implicit real*8 (a-h,o-z)
      common nq,en,enex,smax,alam,delta,v0,v1
      v0=-1
      v1=alam*delta
      v=-1/x+alam*(1-exp(-delta*x))/x
      return
      end
            real*8 function vdif(x)
      implicit real*8 (a-h,o-z)
      common nq,en,enex,smax,alam,delta,v0,v1
      vdif=(1-alam)/x**2+alam*exp(-delta*x)/x**2*(1+delta*x)
      return
      end
            real*8 function gcond(t,y)
      implicit real*8 (a-h,o-z)
      dimension y(1)
      common nq,en,enex,smax,alam,delta,v0,v1
      gcond=y(1)-smax
      return
      end
            subroutine fcn(t,y,f)
      implicit real*8 (a-h,o-z)
      dimension y(3),f(3)
      common nq,en,enex,smax,alam,delta,v0,v1
	f(1)=y(2)
      f(2)=2*(v(t)-en)*y(1)
      f(3)=y(1)**2
      return
      end
            subroutine fcn1(t,y,f)
      implicit real*8 (a-h,o-z)
      dimension y(1),f(1)
      common nq,en,enex,smax,alam,delta,v0,v1
      a=v(t)-en
      f(1)=0
      if(a.le.0)return
      f(1)=dsqrt(2*a)
      return
      end