133 lines
2.1 KiB
Fortran
133 lines
2.1 KiB
Fortran
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
|