subroutine fit_merge(step)

	implicit none

	include 'fit.inc'

	real step
	real x0, xsum, ysum, ssum, rsum, s, r
	integer i,j

	if (nxmin .ge. nxmax) return
	nset=0
	do i=1,nxmin-1
	  nset=max(nset,iset(i))
	enddo
	do i=nxmax+1,npkt
	  nset=max(nset,iset(i))
	enddo
	do i=nxmin,nxmax
	  if (iset(i) .le. nset) goto 10
	enddo
	nset=nset+1
10	call fit_sort(nxmin, nxmax)
	if (step .eq. 0) then
	  s=(xval(nxmax)-xval(nxmin))/(nxmax-nxmin)/2
	  print *,'merge tolerance ',s
	else
	  s=step
	endif
	i=nxmin
	j=nxmin
	do while (i .le. nxmax)
	  rsum=rmon(i)			! sum of weights
	  xsum=xval(i)*rsum		! sum of xval
	  ysum=YVAL(i)*rsum		! sum of weightet YVAL
	  ssum=sig(i)*sig(i)*rsum*rsum	! sum of weightet sig^2
	  x0=0.75*(xval(i)+s)
	  i=i+1
	  do while (i .le. nxmax .and. xval(i) .lt. x0+xval(i-1)/4)
	    r=rmon(i)
	    rsum=rsum+r
	    xsum=xsum+xval(i)*r
	    ysum=ysum+YVAL(i)*r
	    ssum=ssum+sig(i)*sig(i)*r*r
	    i=i+1
	  enddo
	  if (rsum .gt. 0) then
	    xval(j)=xsum/rsum
	    YVAL(j)=ysum/rsum
	    sig(j)=sqrt(ssum)/rsum
	    rmon(j)=rsum
	    iset(j)=nset
	    j=j+1
	  endif
	enddo
	nxmax=j-1

! shift points back after fit window
	do while (i .le. npkt)
	  xval(j)=xval(i)
	  YVAL(j)=YVAL(i)
	  sig(j)=sig(i)
	  rmon(j)=rmon(i)
	  iset(j)=iset(i)
	  i=i+1
	  j=j+1
	enddo	  
	npkt=j-1
	end


	subroutine fit_sort_load(i)

	include 'fit.inc'

	integer i,j,tst

	integer ibuf
	real buf,el/0/

	el=xval(i)
	return

	entry fit_sort_exch(i,j)

	buf=xval(i)
	xval(i)=xval(j)
	xval(j)=buf
	buf=YVAL(i)
	YVAL(i)=YVAL(j)
	YVAL(j)=buf
	buf=sig(i)
	sig(i)=sig(j)
	sig(j)=buf
	buf=rmon(i)
	rmon(i)=rmon(j)
	rmon(j)=buf
	ibuf=iset(i)
	iset(i)=iset(j)
	iset(j)=ibuf
	return

	entry fit_sort_comp(i,tst)

	if (xval(i) .gt. el) then
	  tst=1
	elseif (xval(i) .lt. el) then
	  tst=-1
	else
	  tst=0
	endif
	end


	subroutine fit_sort(from, to)

	implicit none

	integer from, to

	include 'fit.inc'

	external fit_sort_load, fit_sort_exch, fit_sort_comp

	if (from .eq. 0 .and. to .eq. 0) then
	  call quick_sort(nxmin, nxmax
	1, fit_sort_load, fit_sort_exch, fit_sort_comp, 0)
	else
	  call quick_sort(from, to
	1, fit_sort_load, fit_sort_exch, fit_sort_comp, 0)
	endif
	end