program fit1
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
c     for ARBITRARY ANGULAR MOMENTUM
      implicit real*8 (a-h,o-z)
      dimension tab(2,1000)
      common nq,en,enex,smax,alam,delta,v0,v1,lq,qll
      external func
      print*,'Enter Delta-guess, d-Delta'
      read*,delta0,dd
c      delta0=0.22
c      dd=0.02
      eps=1.d-9
      eta=0
      open(1,file='fit1-inp.dat')
      read(1,*)nel,nq,lq
      qll=lq*(lq+1)/2.d0
      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='fit1-out.dat')
c Input file has a table of (lambda Ei) for V = -1/r1-1/r2 + lambda/r12 found exactly
      alam0=0
      i=0
      do 3 n=ntab,1,-1
      i=i+1
      rz=tab(1,n)
      enz=tab(2,n)
      rlam=1/rz
      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,alam,delta1
  10  format(' -----',2f13.8)
      write(2,20)alam,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,lq,qll
      external coeffn,bdyval,monit,report,fcn,fcn1,gcond
      pi=4*datan(1.d0)
      toldefa=1.d-9
      delta=delta1
      tol0=toldefa
      tolg=tol0*10000
      toll=tol0/4
      x0=tol0**(1.d0/4)/3
  60  format(9x,a19,e10.2)
      ind=nq-lq-1
c      en0=-(1-alam)**2/2/nq**2
c      if(abs(en0).lt.1.d-3)en0=-1.d-3
      en0=-enex
      en=en0
      d=abs(en0)
      if(d.lt.1.d-3/nq**2)d=1.d-3/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
   2  if(veff(x2).gt.en)goto1
      x2=0.9*x2
      if(x2.lt.x0)print*,'Root of V(r) was not found!'
      goto2
   1  if(veff(x1).gt.en .or. veff(x2).lt.en)
     .   print*,'Root was not localised!'
      do 101 iter=1,100
      x=(x1+x2)/2
      vx=veff(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
  12  if(veff(x2).gt.en)goto11
      x2=0.9*x2
      if(x2.lt.x0)print*,'2 Root of V(r) was not found!'
      goto12
  11  if(veff(x1).gt.en .or. veff(x2).lt.en)
     .   print*,'2 Root was not localised!'
      do 102 iter=1,100
      x=(x1+x2)/2
      vx=veff(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
      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,'del,lam,en,ex',delta,alam,elam,enex
      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=veff(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,lq,qll
      h00=1
      h01=v0*h00/(lq+1)
      h02=(v0*h01+(v1-en)*h00)/(2*lq+3)
      fun=xl**(lq+1)*(h00+h01*xl+h02*xl**2)
      dfun=xl**lq*((lq+1)*h00
     .    +(lq+2)*h01*xl+(lq+3)*h02*xl**2)
      yl(1)=fun
      yl(2)=dfun/2
      yr(1)=1
      p2=2*(veff(xr)-elam)
      dp2=2*veffdif(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,lq,qll
      v0=-1
      v1=alam*delta
      v=-1/x+alam*(1-exp(-delta*x))/x
      return
      end
            real*8 function veff(x)
      implicit real*8 (a-h,o-z)
      common nq,en,enex,smax,alam,delta,v0,v1,lq,qll
c  !!!      
      v0=-1+alam
      v1=alam*delta
      veff=qll/x**2-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,lq,qll
      vdif=(1-alam)/x**2+alam*exp(-delta*x)/x**2*(1+delta*x)
      return
      end
            real*8 function veffdif(x)
      implicit real*8 (a-h,o-z)
      common nq,en,enex,smax,alam,delta,v0,v1,lq,qll
      veffdif=-2*qll/x**3+(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,lq,qll
      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,lq,qll
      f(1)=y(2)
      f(2)=2*(veff(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,lq,qll
      a=veff(t)-en
      f(1)=0
      if(a.le.0)return
      f(1)=dsqrt(2*a)
      return
      end