subroutine bend_fit(nl, za,bend, nummin,nummax,zbeg,zend,bendfit,
     &                  ScaleH)
c   
c
      implicit none
      integer nl,i,jj,jj1,jj2,Ibeg,Iend,idelta,numdelta
      integer nummin,nummax,num,nhalf,ndiff
      real*4 ZA(nl),BEND(nl),ScaleH(nl),bendfit(nl),zbeg,zend,slope
      REAL*8 DeltaZ,anum,sumx,sumxsq,sumlna,sumxlna,suma,sumxa,
     &sumxsqa,sumalna,sumxalna
      do i=1,nl
        if(za(i).ge.zbeg) Ibeg=i
        if(za(i).ge.zend) Iend=i
      enddo
      idelta=Iend-Ibeg
      numdelta=nummax-nummin
      slope=real(numdelta)/real(idelta)
c      write(*,*) Ibeg,Iend,idelta,numdelta,slope,nummin,nummax
c       stop
      do i=2,nl-1
        sumx=0.0D0
        sumxsq=0.0D0
        sumlna=0.0D0
        sumxlna=0.0D0
        suma=0.0D0
        sumxa=0.0D0
        sumxsqa=0.0D0
        sumalna=0.0D0
        sumxalna=0.0D0
        IF(numdelta.gt.0) then
          if(i.le.ibeg) THEN
            nhalf=nummax
          elseif(i.ge.iend) THEN
            nhalf=nummin
          else
            nhalf=nummax-INT(slope*(i-Ibeg))
          endif
        ELSE
          nhalf=nummax
        ENDIF
        ndiff=min(i-1,nl-i)
        if(nhalf.gt.ndiff) nhalf=ndiff
        if(nhalf.lt.1) nhalf=1
        num=1+2*nhalf
        anum=dble(num)
        jj1=i-nhalf
        jj2=i+nhalf
        DO jj=jj1,jj2
          DeltaZ=za(jj)-za(i)
          sumx=sumx+DeltaZ
          sumxsq=sumxsq+Deltaz**2
          sumlna=sumlna+alog(bend(jj))
          sumxlna=sumxlna+DeltaZ*alog(bend(jj))
          suma=suma+bend(jj)
          sumxa=sumxa+DeltaZ*bend(jj)
          sumxsqa=sumxsqa+bend(jj)*DeltaZ**2
          sumalna=sumalna+bend(jj)*alog(bend(jj))
          sumxalna=sumxalna+DeltaZ*bend(jj)*alog(bend(jj))
        enddo  
C        ScaleH(i)=-REAL(((anum*sumxsq-sumx**2)/
C     &             (anum*sumxlna-sumx*sumlna)))
C        bendfit(i)=exp(REAL((sumlna*sumxsq-sumx*sumxlna)/
C     &             (anum*sumxsq-sumx**2)))
        ScaleH(i)=-REAL((suma*sumxsqa-sumxa**2)/
     &  (suma*sumxalna-sumxa*sumalna))
        bendfit(i)=exp(REAL((sumxsqa*sumalna-sumxa*sumxalna)/
     &             (suma*sumxsqa-sumxa**2)))
C        write(*,*) ZA(I),bendfit(i),ScaleH(I),anum
      ENDDO
      bendfit(1)=bend(1)
      ScaleH(1)=ScaleH(2)
      bendfit(nl)=bend(nl-1)
      ScaleH(nl)=ScaleH(nl-1)
      RETURN
      END