program fit1
c  Finding delta that gives zero eigenvalue in the potential  V(r) = -1/r+lambda(1-exp(-delta r))/r
c     for ARBITRARY ANGULAR MOMENTUM
      implicit real*8 (a-h,o-z)
      common nq,en,enex,smax,alam,delta,v0,v1,lq,qll,nel
      common /y1y2/y1,y2
      external func
      print*,'Enter Lambda-guess, dlambda'
      read*,alam0,dd
      eps=1.d-8
      eta=0
      open(1,file='zero1.inp')
c      print*,'Enter Number of electrons, Quantum numbers'
      read(1,*)nel,nq,lq
c      print*,'Enter Delta0, Delta1 - for neutral atom and negative ion'
      read(1,*)y1,y2
c      
      qll=lq*(lq+1)/2.d0
      ifail=0
      call c05agf(alam0,dd,eps,eta,func,a,b,ifail)
      if(ifail.ne.0)print*,'IFAIL =',ifail
      print10,alam0,(nel-1)/alam0
  10  format(' -----',2f13.8)
  20  format(1x,2f13.8)
      stop      
      end
      
      real*8 function func(alam0)           
      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,nel
      common /y1y2/y1,y2
      external coeffn,bdyval,monit,report,fcn,fcn1,gcond
      pi=4*datan(1.d0)
      toldefa=1.d-8
      alam=alam0
c Define Delta as a function of Lambda!
c      if(nel.eq.2)delta=1.26422-0.408707*alam+0.0251437*alam**2
c      if(nel.eq.2)delta=1.25175-0.371094*alam
c      if(nel.eq.2)delta=1.26461-0.410073*alam+0.0261292*alam**2
c      if(nel.eq.4)delta=0.530339-0.215727*alam-0.0557214*alam**2
c      if(nel.eq.12)delta=0.486165
c     .    -0.870471*alam+1.10241*alam**2-0.588113*alam**3      
c      if(nel.eq.6)delta=0.402537-0.140678*alam-0.0441141*alam**2
c      if(nel.eq.7)delta=0.384966-0.139589*alam-0.0326446*alam**2
c      if(nel.eq.9)delta=0.396825-0.15399*alam-0.0273002*alam**2
c      if(nel.eq.10)delta=0.378044-0.131239*alam-0.0352188*alam**2
c      if(nel.eq.2 .and. nq.eq.2 .and. lq.eq.1)
c     .    delta=0.304889-0.08706*alam
c      if(nel.eq.14)delta=0.225918+0.015256*alam-0.129663*alam**2
c      if(nel.eq.15)delta=0.224682-0.0264879*alam-0.0885347*alam**2
c      if(nel.eq.16)delta=0.212692+0.00421499*alam-0.105775*alam**2
c      if(nel.eq.17)delta=0.208343-0.0158098*alam-0.0832303*alam**2
c      if(nel.eq.18)delta=0.20003-0.0201667*alam-0.0720394*alam**2
c      if(nel.eq.20)delta=0.372091-0.297303*alam
c      if(nel.eq.21)delta=0.273056-0.205269*alam
      x1=1-1.d0/nel
      x2=1
      delta=(y1*(alam-x2)-y2*(alam-x1))/(x1-x2)
      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=0
      en=en0
c      d=(1-alam)**2/2/nq**2
      d=1.d-4
  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
        print*,'xend1 ',xend
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      en=elam
c      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
        print*,'xend2 ',xend
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      en=elam
c      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)
      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
        print*,'xend3 ',xend
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      print10,'lam,del,xend,en',alam,delta,xend,elam
      en=elam
  40  format(f10.5,f16.8)
  10  format(1x,a15,3f12.6,f15.9)
      func=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,nel
      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,nel
      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,nel
      v0=-1
      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,nel
      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,nel
      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,nel
      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,nel
      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,nel
      a=veff(t)-en
      f(1)=0
      if(a.le.0)return
      f(1)=dsqrt(2*a)
      return
      end