499 lines
9.9 KiB
Fortran
499 lines
9.9 KiB
Fortran
subroutine fit_userfun(name, ufun)
|
|
! ----------------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
character*(*) name
|
|
external ufun
|
|
real ufun
|
|
integer sys_adrr_rriii
|
|
|
|
usernp=0
|
|
npd=0
|
|
usertit=name
|
|
userfun = sys_adrr_rriii(ufun)
|
|
end
|
|
|
|
|
|
subroutine fit_userpar(name)
|
|
! ----------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
character*(*) name
|
|
|
|
integer i
|
|
|
|
if (npd .ne. 0)
|
|
1 stop 'per-dataset parameters must be defined at end!'
|
|
usernp=usernp+1
|
|
if (usernp .gt. maxext) stop 'too many user parameters'
|
|
i=index(name,':')
|
|
userpar(usernp)=name(i+1:)
|
|
if (i .gt. 1) then
|
|
call str_upcase(usersho(usernp),name(1:i-1))
|
|
else
|
|
usersho(usernp)=' '
|
|
endif
|
|
usercct(usernp)=cinfo
|
|
end
|
|
|
|
|
|
subroutine fit_userpdp(name)
|
|
! ----------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
character*(*) name
|
|
|
|
integer i
|
|
|
|
npd=npd+1
|
|
if (npd .gt. maxpd .or. npd+usernp .gt. maxext)
|
|
1 stop 'too many per-dataset parameters'
|
|
userpar(usernp+npd)=name
|
|
do i=1,maxset
|
|
pdpar(npd, i)=none
|
|
enddo
|
|
c usercct(usernp+npd)=cinfo ! not yet implemented
|
|
end
|
|
|
|
|
|
subroutine fit_usercinfo(calcinfo)
|
|
! ----------------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
integer calcinfo
|
|
|
|
cinfo=calcinfo
|
|
end
|
|
|
|
|
|
|
|
subroutine fit_ncurves(ncur)
|
|
! ----------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
integer ncur
|
|
integer i
|
|
|
|
nset=ncur
|
|
if (nset .le. 0) then
|
|
do i=1,npkt
|
|
nset=max(nset,iset(i))
|
|
enddo
|
|
endif
|
|
end
|
|
|
|
|
|
|
|
subroutine fit_user1st
|
|
|
|
include 'fit.inc'
|
|
|
|
real dum
|
|
real sys_rfun_rriii
|
|
|
|
integer calcinfo,i
|
|
|
|
calcinfo=0
|
|
do i=1,nu
|
|
if (u(i) .ne. pold(i)) then
|
|
calcinfo=ior(calcinfo, usercct(i))
|
|
endif
|
|
enddo
|
|
|
|
dum=sys_rfun_rriii(userfun, 0.0, u, nu, 1, calcinfo)
|
|
end
|
|
|
|
|
|
integer function fit_userini(jfifu)
|
|
|
|
include 'fit.inc'
|
|
|
|
integer i,jfifu
|
|
real dum, xx
|
|
|
|
real sys_rfun_rriii
|
|
|
|
do i=1,nu
|
|
pold(i)=1.234e-35
|
|
enddo
|
|
xx=-float(jfifu)
|
|
dum=sys_rfun_rriii(userfun, xx, u, nu, -1, 0)
|
|
fit_userini=nint(xx)
|
|
if (fit_userini .lt. 0) fit_userini=0
|
|
end
|
|
|
|
|
|
subroutine fit_check_range
|
|
|
|
include 'fit.inc'
|
|
|
|
real xmin, xmax
|
|
real xmin0/-1./,xmax0/1.0/
|
|
integer i
|
|
|
|
if (ififu .ne. 7) RETURN
|
|
xmin0=xval(nxmin)
|
|
xmax0=xmin0
|
|
do i=nxmin,nxmax
|
|
if (xval(i) .lt. xmin0) xmin0=xval(i)
|
|
if (xval(i) .gt. xmax0) xmax0=xval(i)
|
|
enddo
|
|
return
|
|
|
|
entry fit_put_xrange(xmin, xmax)
|
|
|
|
xmin0=min(xmin,xmax)
|
|
xmax0=max(xmin,xmax)
|
|
RETURN
|
|
|
|
entry fit_limit_xrange(xmin, xmax)
|
|
|
|
xmin0=max(xmin0,min(xmin,xmax))
|
|
xmax0=min(xmax0,max(xmin,xmax))
|
|
RETURN
|
|
|
|
entry fit_get_xrange(xmin, xmax)
|
|
|
|
xmin=xmin0
|
|
xmax=xmax0
|
|
end
|
|
|
|
|
|
subroutine fit_fixit(i)
|
|
|
|
include 'fit.inc'
|
|
|
|
integer i
|
|
|
|
lcode(i)=0
|
|
end
|
|
|
|
|
|
subroutine fit_delta(xd,sd)
|
|
|
|
include 'fit.inc'
|
|
|
|
real sd,xd
|
|
logical read
|
|
|
|
integer nd,nr
|
|
real xdi(maxpeak), sdi(maxpeak)
|
|
|
|
logical disable/.true./
|
|
save nd,nr,xdi,sdi,disable
|
|
|
|
if (sd .eq. 0 .or. disable) RETURN
|
|
if (nd .ge. maxpeak) then
|
|
print *,'too many delta peaks'
|
|
RETURN
|
|
endif
|
|
nd=nd+1
|
|
xdi(nd)=xd
|
|
sdi(nd)=sd
|
|
return
|
|
|
|
entry fit_init_delta(read)
|
|
if (read) then
|
|
nr=0
|
|
else
|
|
nd=0
|
|
endif
|
|
disable=read
|
|
return
|
|
|
|
entry fit_get_delta(xd,sd)
|
|
if (nr .ge. nd) then
|
|
sd=0
|
|
xd=0
|
|
else
|
|
nr=nr+1
|
|
sd=sdi(nr)
|
|
xd=xdi(nr)
|
|
endif
|
|
end
|
|
|
|
|
|
subroutine fit_confun(name, ufun)
|
|
! ---------------------------------
|
|
|
|
include 'fit.inc'
|
|
include 'fit_user.inc'
|
|
|
|
character*(*) name
|
|
external ufun, fit_cfun
|
|
real ufun
|
|
real fit_cfun
|
|
integer sys_adrr_rriii
|
|
|
|
data xbase/0.0/
|
|
data show/.false./
|
|
data ecnt/0/
|
|
|
|
call fit_userfun(name, ufun)
|
|
confun=userfun
|
|
userfun = sys_adrr_rriii(fit_cfun)
|
|
call fit_userpar('GW:GaussWid')
|
|
call fit_userpar('XS:CalcStep')
|
|
end
|
|
|
|
|
|
subroutine fit_xbase(x_base)
|
|
! ----------------------------
|
|
|
|
real x_base
|
|
|
|
include 'fit_user.inc'
|
|
|
|
xbase=x_base
|
|
end
|
|
|
|
|
|
real function fit_cfun(x,p,n,mode,cinfo)
|
|
! ----------------------------------------
|
|
|
|
implicit none
|
|
|
|
real x ! x-value
|
|
integer n ! number of parameters
|
|
real p(*) ! parameters
|
|
integer mode ! mode ! 0: x-dependent part, 1: x-independent part, -1: par. independent part
|
|
integer cinfo ! calculation information (not used here)
|
|
|
|
include 'fit_user.inc'
|
|
|
|
integer maxt, maxg, ng, imaxg, m_f, m_i, m_t, m_x
|
|
real pn, gn
|
|
parameter (m_f=0, m_i=1, m_t=2, m_x=3)
|
|
! m_f: call function directly
|
|
! m_i: initialize function
|
|
! m_t: transform x -> t
|
|
! m_x: transform t -> x
|
|
parameter (pn=2.0) ! step size of larger gaussian = fwhm/pn
|
|
parameter (gn=2.0) ! total width = gn*fwhm
|
|
parameter (imaxg=128) ! maximal fwhm (in steps)
|
|
integer i,j,k,m
|
|
|
|
real ts ! t-step
|
|
real t ! transformed x
|
|
real ta, tb ! limits of t-space
|
|
|
|
real gw ! total gaussian width (in t-space)
|
|
real wg ! total peak weight
|
|
real sum ! summing variable
|
|
integer kstp ! step for folding
|
|
integer mfold ! number of folding steps
|
|
real wd ! individual peak width
|
|
real t0 ! base for integration
|
|
real lasty ! last value for integr.
|
|
|
|
real xa, xb ! x-range
|
|
real xs ! x-step at xbase
|
|
real q ! ratio between fwhm's of peaks for folding steps
|
|
real stp ! real step
|
|
real xd, sd ! delta peak parameters
|
|
|
|
real x0, x1, x2, qq, td, xsi
|
|
|
|
parameter (maxt=10000) ! maximal number of points
|
|
real y(0:maxt),yf(0:maxt) ! array before and after folding
|
|
|
|
parameter (maxg=100) ! maximal number of points per gaussian (> gn*5)
|
|
real gs(0:maxg) ! precalculated gaussian
|
|
|
|
save y,ta,tb,xa,xb,gw,kstp
|
|
|
|
external sys_rfun_rriii, voigt, spline3, spline4
|
|
real sys_rfun_rriii, voigt, spline3, spline4
|
|
|
|
if (mode .eq. 0) then
|
|
if (gw .eq. 0) then
|
|
fit_cfun=sys_rfun_rriii(confun,x,p,n,m_f,cinfo)
|
|
RETURN
|
|
endif
|
|
if (x .lt. xa .or. x .gt. xb) then
|
|
fit_cfun=sys_rfun_rriii(confun,x,p,n,m_f,cinfo)
|
|
ecnt=ecnt+1
|
|
show=.true.
|
|
else
|
|
t=sys_rfun_rriii(confun,x,p,n,m_t,cinfo)/ts
|
|
t=(t-ta)/istp
|
|
i=nint(t-0.5)
|
|
if (i .le. 1) then
|
|
fit_cfun=spline3(t, y(0), y(1), y(2))
|
|
elseif (i .ge. npy-1) then
|
|
fit_cfun=spline3(npy-t, y(npy), y(npy-1), y(npy-2))
|
|
else
|
|
fit_cfun=spline4(t-i, y(i-1), y(i), y(i+1), y(i+2))
|
|
endif
|
|
endif
|
|
else if (mode .gt. 0) then
|
|
show=.true.
|
|
call fit_init_delta(.false.)
|
|
if (p(1) .eq. 0.0) then
|
|
fit_cfun=sys_rfun_rriii(confun,x,p,n,m_f,cinfo)
|
|
RETURN
|
|
endif
|
|
call fit_get_xrange(xa,xb)
|
|
fit_cfun=sys_rfun_rriii(confun,x,p,n,m_i,cinfo)
|
|
call fit_init_delta(.true.)
|
|
ta=sys_rfun_rriii(confun,xa,p,n,m_t,cinfo)
|
|
tb=sys_rfun_rriii(confun,xb,p,n,m_t,cinfo)
|
|
xs=p(2)
|
|
if (xs .eq. 0) xs=p(1)/10
|
|
if (xbase-max(p(1),xs) .lt. xa .or.
|
|
1 xbase+max(p(1),xs) .gt. xb) then
|
|
qq=(tb-ta)/(xb-xa)
|
|
ts=abs(xs*qq)
|
|
gw=abs(p(1)*qq)/ts
|
|
else
|
|
ts=abs(sys_rfun_rriii(confun,xbase+xs,p,n,m_t,cinfo)
|
|
1 -sys_rfun_rriii(confun,xbase-xs,p,n,m_t,cinfo))
|
|
if (ts .eq. 0) ts=1.0
|
|
gw=abs(sys_rfun_rriii(confun,xbase+p(1)/2,p,n,m_t,cinfo)
|
|
1 -sys_rfun_rriii(confun,xbase-p(1)/2,p,n,m_t,cinfo))/ts
|
|
endif
|
|
ta=int((ta/ts-gn*gw))
|
|
tb=int((tb/ts+gn*gw)+1)
|
|
npy=nint(tb-ta)
|
|
istp=max(int(gw/imaxg+1), (npy+maxt-1)/maxt)
|
|
npy=npy/istp
|
|
|
|
t0=ta-0.5*istp
|
|
x0=sys_rfun_rriii(confun,t0*ts,p,n,m_x,cinfo)
|
|
lasty=sys_rfun_rriii(confun,x0,p,n,m_f,cinfo)
|
|
k=0
|
|
do i=0,npy
|
|
x1=sys_rfun_rriii(confun,(t0+i*istp)*ts,p,n,m_x,cinfo)
|
|
! simpson integral over interval (x0...x1)
|
|
xsi=(x1-x0)/istp
|
|
x2=x1-xsi/2
|
|
sum=lasty
|
|
do j=1,istp
|
|
lasty=sys_rfun_rriii(confun,x0+j*xsi,p,n,m_f,cinfo)
|
|
sum=sum+4*sys_rfun_rriii(confun,x2+j*xsi,p,n,m_f,cinfo)
|
|
1 +2*lasty
|
|
enddo
|
|
x0=x1
|
|
y(i)=(sum-lasty)/6/istp
|
|
k=k+istp
|
|
enddo
|
|
call fit_get_delta(xd, sd)
|
|
do while (sd .ne. 0)
|
|
td=sys_rfun_rriii(confun,xd,p,n,m_t,cinfo)/ts
|
|
sd=sd/abs(
|
|
1 sys_rfun_rriii(confun,(td+0.5*istp)*ts,p,n,m_x,cinfo)
|
|
1 -sys_rfun_rriii(confun,(td-0.5*istp)*ts,p,n,m_x,cinfo) )
|
|
if (td .ge. ta .and. td .le. tb) then
|
|
t=(td-ta)/istp
|
|
i=nint(t-0.4999)
|
|
if (i .ge. npy) i=npy-1
|
|
y(i)=y(i)+sd*(i+1-t)
|
|
y(i+1)=y(i+1)+sd*(t-i)
|
|
endif
|
|
call fit_get_delta(xd, sd)
|
|
enddo
|
|
|
|
if (gw .le. istp) RETURN
|
|
|
|
mfold=max(1,nint(log(gw/pn/istp)))
|
|
sum=0
|
|
q=max(1.1,(gw/pn/istp)**(1.0/mfold))
|
|
qq=1/(q*q)
|
|
qq=sqrt((1-qq**mfold)/(1-qq)) ! sqrt of sum of squared widths divided by gw
|
|
wd=gw/istp/q**(mfold-1)/qq
|
|
|
|
icnt=0
|
|
stp=1.0
|
|
do m=1,mfold
|
|
kstp=int(stp)
|
|
stp=stp*q
|
|
ng=min(maxg,int(gn*wd/kstp+0.9))
|
|
icnt=icnt+(2*ng+1)*(npy+1)-(ng-1)*ng+1
|
|
do i=0,ng
|
|
gs(i)=voigt(float(i), wd/kstp, 0.0)
|
|
enddo
|
|
do i=0,npy
|
|
wg=0
|
|
sum=0
|
|
k=i
|
|
do j=0,min(ng,i/kstp)
|
|
sum=sum+gs(j)*y(k)
|
|
wg=wg+gs(j)
|
|
k=k-kstp
|
|
enddo
|
|
k=i+kstp
|
|
do j=1,min(ng,(npy-i)/kstp)
|
|
sum=sum+gs(j)*y(k)
|
|
wg=wg+gs(j)
|
|
k=k+kstp
|
|
enddo
|
|
yf(i)=sum/wg
|
|
enddo
|
|
do i=ng*kstp,npy-ng*kstp
|
|
y(i)=yf(i)
|
|
enddo
|
|
wd=wd*q
|
|
enddo
|
|
else
|
|
fit_cfun=sys_rfun_rriii(confun,x,p,n,mode,cinfo)
|
|
call fit_fixit(2)
|
|
endif
|
|
end
|
|
|
|
|
|
subroutine fit_confun_print
|
|
|
|
include 'fit_user.inc'
|
|
|
|
if (show) then
|
|
print *
|
|
print '(i10,a,i9,a)'
|
|
1 , npy*istp*2+1, ' fn. calls, fold. multiplications: ',icnt
|
|
if (ecnt .gt. 0) then
|
|
print *,ecnt,' outside range calls counted'
|
|
ecnt=0
|
|
endif
|
|
show=.false.
|
|
endif
|
|
end
|
|
|
|
|
|
real function spline4(x, y1, y2, y3, y4)
|
|
|
|
! spline interpolation.
|
|
! the result is valid for values x between 0 and 1 (x-coord of y2 and y3)
|
|
|
|
real x, y1, y2, y3, y4
|
|
real ss, sl, sr
|
|
|
|
ss=y3-y2 ! median slope
|
|
sl=(y3-y1)/2 ! slope at point 2
|
|
sr=(y4-y2)/2 ! slope at point 3
|
|
spline4=y2+x*ss ! linear interpol.
|
|
1 + (sl-ss)*x*(x-1)**2 ! correction for left slope
|
|
1 - (sr-ss)*(1-x)*x**2 ! correction for right slope
|
|
end
|
|
|
|
|
|
real function spline3(x, y1, y2, y3)
|
|
|
|
! spline interpolation.
|
|
! the result is valid for values x between 0 and 1 (x-coord of y1 and y2)
|
|
! or, as an extrapolation to the left of x1
|
|
|
|
real x, y1, y2, y3
|
|
real sl, sr, s1
|
|
|
|
s1=y2-y1
|
|
sr=(y3-y1)/2 ! interpolate slope at point 2
|
|
sl=2*sr-s1 ! extrapolate slope at point 1
|
|
spline3=y1+x*s1 ! linear interpol.
|
|
1 + (sl-s1)*x*(x-1)**2 ! correction for left slope
|
|
1 - (sr-s1)*(1-x)*x**2 ! correction for right slope
|
|
end
|