298 lines
4.4 KiB
Fortran
298 lines
4.4 KiB
Fortran
subroutine fit_min(narg)
|
|
! ------------------------
|
|
|
|
include 'fit.inc'
|
|
integer narg, nf
|
|
real earg
|
|
|
|
call fit_set_inthdl
|
|
call fit_check_range
|
|
|
|
nfcnmx=1000
|
|
if (narg .gt. 0) nfcnmx=narg
|
|
isw(1)=0
|
|
|
|
NF = NFCN
|
|
CALL SIMPLEX
|
|
if (ISW(1) .lt. 1) then
|
|
NFCNMX = NFCNMX + NF - NFCN
|
|
GO TO 360
|
|
endif
|
|
|
|
return
|
|
|
|
|
|
entry fit_fit(narg)
|
|
! -------------------
|
|
|
|
call fit_set_inthdl
|
|
call fit_check_range
|
|
|
|
nfcnmx=1000
|
|
if (narg .gt. 0) nfcnmx=narg
|
|
isw(1)=0
|
|
|
|
360 NF = NFCN
|
|
APSI = EPSI
|
|
CALL MIGRAD
|
|
IF (ISW(2) .LE. 2 .AND. ISW(1) .NE. 1) THEN
|
|
NFCNMX = NFCNMX + NF - NFCN
|
|
NF = NFCN
|
|
CALL SIMPLEX
|
|
ENDIF
|
|
return
|
|
|
|
|
|
entry fit_sim(narg)
|
|
! -------------------
|
|
|
|
call fit_set_inthdl
|
|
call fit_check_range
|
|
|
|
nfcnmx=1000
|
|
if (narg .gt. 0) nfcnmx=narg
|
|
isw(1)=0
|
|
|
|
call simplex
|
|
return
|
|
|
|
|
|
entry fit_mig(narg)
|
|
! -------------------
|
|
|
|
call fit_set_inthdl
|
|
call fit_check_range
|
|
|
|
nfcnmx=1000
|
|
if (narg .gt. 0) nfcnmx=narg
|
|
isw(1)=0
|
|
|
|
NF = NFCN
|
|
APSI = EPSI
|
|
CALL MIGRAD
|
|
return
|
|
|
|
|
|
entry fit_err(earg)
|
|
! -------------------
|
|
|
|
if (earg .gt. 0) then
|
|
epsi=epsi*earg/up
|
|
up=earg
|
|
else
|
|
epsi=epsi/up
|
|
UP=1.0
|
|
endif
|
|
IF (ISW(2).GE.1) CALL fit_print(1)
|
|
return
|
|
|
|
|
|
entry fit_pri(narg)
|
|
! -------------------
|
|
|
|
ISW(5) = narg
|
|
return
|
|
|
|
|
|
entry fit_epsi(earg)
|
|
! --------------------
|
|
|
|
if (earg .gt. 0) then
|
|
epsi=earg
|
|
else
|
|
epsi=0.1*up
|
|
endif
|
|
return
|
|
|
|
|
|
entry fit_vtest(earg)
|
|
! ---------------------
|
|
if (earg .gt. 0) then
|
|
vtest=earg
|
|
else
|
|
vtest=0.01
|
|
endif
|
|
return
|
|
|
|
end
|
|
|
|
|
|
subroutine fit_set_inthdl
|
|
|
|
include 'fit.inc'
|
|
external fit_inthdl
|
|
|
|
if (isyswr .ne. 0)
|
|
1 write(isyswr,'(t40,a)') 'Press Ctrl-C to abort fit-algorithm'
|
|
|
|
call sys_int_hdl(fit_inthdl)
|
|
end
|
|
|
|
subroutine fit_inthdl
|
|
include 'fit.inc'
|
|
nfcnmx=0
|
|
end
|
|
|
|
|
|
subroutine fit_reserr
|
|
|
|
implicit none
|
|
|
|
include 'fit.inc'
|
|
|
|
integer i
|
|
|
|
do i=1,nu
|
|
if (werr(i) .ne. 0) werr(i)=max(0.1,u(i)*0.1)
|
|
enddo
|
|
if (ififu .eq. 1) then
|
|
do i=3,nu,5
|
|
if (werr(i+1) .ne. 0) werr(i+1)=werr(i+3)+werr(i+4)
|
|
enddo
|
|
endif
|
|
call fit_set(0,0.0,0.0,0.0,0.0)
|
|
end
|
|
|
|
|
|
|
|
subroutine fit_chisq(chisq, istat)
|
|
|
|
real chisq
|
|
integer istat
|
|
|
|
include 'fit.inc'
|
|
|
|
chisq=amin
|
|
istat=isw(2)
|
|
end
|
|
|
|
|
|
|
|
subroutine fit_restore_wr
|
|
|
|
include 'fit.inc'
|
|
|
|
integer isave/0/
|
|
|
|
if (isave .ne. 0) then
|
|
isyswr=isave
|
|
isave=0
|
|
endif
|
|
return
|
|
|
|
entry fit_suspend_wr
|
|
if (isyswr .ne. 0) then
|
|
isave=isyswr
|
|
isyswr=0
|
|
endif
|
|
end
|
|
|
|
|
|
subroutine fit_covar
|
|
|
|
include 'fit.inc'
|
|
integer i,j,k
|
|
character line*130
|
|
|
|
if (isw(2) .lt. 3) then
|
|
print *,'covariance matrix undefined'
|
|
return
|
|
endif
|
|
do i=1,npar
|
|
do k=1,npar,16
|
|
write(line, '(16g8.2)') (v(i,j),j=k,min(k+15,npar))
|
|
! (100*v(i,j)/sqrt(abs(v(i,i)*v(j,j)))
|
|
! if (i .ge. k .and. i .le. k+15) then
|
|
! nam=' '
|
|
! do l=1,nu
|
|
! if (lcorsp(l) .eq. i) then
|
|
! nam=psho(l)
|
|
! goto 9
|
|
! endif
|
|
! enddo
|
|
!9 line((i-k)*8+1:(i-k+1)*8)=' '//nam
|
|
! endif
|
|
print *,line
|
|
enddo
|
|
if (i .gt. 26) print *
|
|
enddo
|
|
end
|
|
|
|
|
|
subroutine fit_true_err(ipar, fstep)
|
|
|
|
integer ipar
|
|
real fstep
|
|
|
|
include 'fit.inc'
|
|
|
|
real sums0, sums, uc, start, step, step0, sl0, sl1
|
|
real p0
|
|
integer i,idir,j
|
|
real usave(maxext), esave(maxext), trum, trup
|
|
|
|
external fit_restore_wr
|
|
|
|
step0=abs(werr(ipar))
|
|
if (step0 .eq. 0) then
|
|
print *,pnam(ipar),' is not free'
|
|
return
|
|
endif
|
|
call sys_err_hdl(fit_restore_wr)
|
|
call fit_suspend_wr
|
|
|
|
call fit_min(0)
|
|
call fit_fix(ipar)
|
|
call fit_min(0)
|
|
do i=1,nu
|
|
usave(i)=u(i)
|
|
esave(i)=werr(i)
|
|
enddo
|
|
! uc=up*max(1.0,amin)/nfree
|
|
uc=up/nfree
|
|
sums0=amin/uc
|
|
start=u(ipar)
|
|
print *,pnam(ipar),' chi^2'
|
|
print *,start,amin
|
|
|
|
do idir=-1,1,2
|
|
step=step0*idir
|
|
sums=sums0
|
|
do j=1,100
|
|
sl0=sqrt(max(0.0,sums-sums0))
|
|
do i=1,nu
|
|
werr(i)=esave(i)
|
|
enddo
|
|
p0=u(ipar)
|
|
call fit_set(ipar,start+step*j*fstep,0.0,0.0,-1.0)
|
|
call fit_min(0)
|
|
if (nfcnmx .eq. 0) goto 99
|
|
sums=amin/uc
|
|
print *,u(ipar),amin,sums-sums0
|
|
if (sums .gt. sums0+1) goto 50
|
|
enddo
|
|
50 sl1=sqrt(max(0.0,sums-sums0))
|
|
if (sl1 .le. sl0) then
|
|
print *,'error'
|
|
goto 99
|
|
endif
|
|
step=p0-start+(u(ipar)-p0)*(1-sl0)/(sl1-sl0)
|
|
if (idir .lt. 0) then
|
|
trum=step
|
|
else
|
|
trup=step
|
|
endif
|
|
do i=1,nu
|
|
u(i)=usave(i)
|
|
werr(i)=esave(i)
|
|
enddo
|
|
call fit_set(0,0.0,0.0,0.0,0.0)
|
|
call fit_min(0)
|
|
enddo
|
|
|
|
99 call fit_rel(ipar)
|
|
call fit_set(ipar, start, step0, 0.0, -1.0)
|
|
call fit_restore_wr
|
|
print *,'true error of ',pnam(ipar),':', trum, trup
|
|
end
|