471 lines
12 KiB
Fortran
471 lines
12 KiB
Fortran
SUBROUTINE fit_fun(funarg, nargs, parg, earg)
|
|
! --------------------------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
integer funarg ! negative: ask, other values: see format 1100 at bottom
|
|
integer nargs ! functions 2,3,4,7:
|
|
! nargs=0: ask for parameters
|
|
! nargs>0: preset parameters and errors
|
|
! functions 0,1:
|
|
! nargs=0: auto parameters
|
|
! nargs=1: preset position, auto width/intensity
|
|
! function 5: narg unused
|
|
! function 6: auto parameters
|
|
! nargs=0: normal, nargs=1: multiply with with parg(1)
|
|
real parg(*) ! preset parameter
|
|
real earg(*) ! preset error
|
|
logical fixed, strange
|
|
|
|
integer npfix,i,iexp,nsum,npuser,k,npeak
|
|
integer jfifu, nu_user, narg
|
|
real a,b,ukp,ukm,sumx,sumy,sumxx,sumxy,xr,yr,xl,yl,uk
|
|
real wk, YVALmax, wmult
|
|
|
|
real voigt ! function
|
|
integer fit_userini
|
|
|
|
jfifu=ififu
|
|
nu_user=0
|
|
if (up .eq. 0.0)
|
|
1 stop 'error in FIT_FUN: FIT_INIT must be called first'
|
|
13 if (funarg .lt. 0) then
|
|
WRITE (ISYSWR,1100) usertit
|
|
READ (ISYSRD,1101,END=999,err=999) i
|
|
IF (i.LT.0 .OR. i.GT.8) GO TO 13
|
|
ififu=i
|
|
narg=0
|
|
else
|
|
ififu=funarg
|
|
narg=nargs
|
|
endif
|
|
|
|
wmult=1.0
|
|
if (ififu .eq. 6) then ! for strange peak use gaussian to determine start values
|
|
if (narg .eq. 1) wmult=parg(1)
|
|
narg=0
|
|
ififu=0
|
|
strange=.true.
|
|
else
|
|
strange=.false.
|
|
endif
|
|
|
|
10 NFCN=1
|
|
ISW(1)=0
|
|
ISW(2)=0
|
|
ISW(3)=0
|
|
SIGMA = 0.
|
|
NPFIX = 0
|
|
NPAR = 0
|
|
DO 11 I= 1, MAXEXT
|
|
if (ififu .ne. 6) then
|
|
! U(I) = 0.0
|
|
WERRS(I) = 0.
|
|
endif
|
|
ICSW(I) = 0
|
|
ICTO(I) = 0
|
|
CFAC(I) = 0.
|
|
coff(I) = 0.
|
|
LCODE(I) = 0
|
|
11 LCORSP (I) = 0
|
|
|
|
IF (IFIFU.LE.1) THEN
|
|
|
|
C . . . . . SINGLE PEAK PARAMETERS . .
|
|
|
|
pnam(1)= 'Bg(Pos) '
|
|
psho(1)= 'B'
|
|
pnam(2)= 'dBg/dx '
|
|
psho(2)= 'D'
|
|
pnam(3)= 'Posi. 1'
|
|
psho(3)= 'P1'
|
|
pnam(4)= 'MaxInt 1'
|
|
psho(4)= 'M1'
|
|
pnam(5)= 'IntInt 1'
|
|
psho(5)= 'I1'
|
|
pnam(6)= 'Fwhm G 1'
|
|
psho(6)= 'G1'
|
|
pnam(7)= 'Fwhm L 1'
|
|
psho(7)= 'L1'
|
|
do i=1,7
|
|
lcode(i) = 1
|
|
enddo
|
|
if (ififu .eq. 0) then ! pure gaussian
|
|
lcode(7)=0
|
|
u(7)=0
|
|
werr(7)=0
|
|
endif
|
|
ififu=1
|
|
if (narg .eq. 1) then ! preselected position
|
|
u(3)=parg(1)
|
|
lcode(3)=0 ! fix pos for findauto
|
|
endif
|
|
nu=2
|
|
call fit_findauto
|
|
nu=7
|
|
if (narg .eq. 1) lcode(3)=1 ! release position
|
|
|
|
elseif (ififu .eq. 6) then
|
|
|
|
C . . . . . . STRANGE PEAK . . . . .
|
|
|
|
PNAM(3)='Position'
|
|
psho(3)='P'
|
|
PNAM(4)='fw(Bg) '
|
|
psho(4)='W'
|
|
NU=4
|
|
WERR(1)=1
|
|
LCODE(1)=1
|
|
WERR(2)=1
|
|
LCODE(2)=1
|
|
WERR(3)=0
|
|
LCODE(3)=0
|
|
WERR(4)=0
|
|
LCODE(4)=0
|
|
if (u(4) .gt. 0) then
|
|
U(4)=u(6)*wmult*sqrt(max(1.0,log(2*u(4)/werr(1))/1.25))
|
|
else
|
|
u(4)=u(6)
|
|
endif
|
|
call extoin
|
|
call fit_fit(0)
|
|
|
|
ELSEIF (IFIFU.EQ.5) THEN
|
|
|
|
C. . . . . CRITCAL EXPONENT PARAMETERS . .
|
|
|
|
do i=1,7
|
|
psho(i)=' '
|
|
enddo
|
|
PNAM(1)= 'TC '
|
|
PNAM(2)= 'BETA '
|
|
PNAM(3)= 'CONST '
|
|
PNAM(4)= 'BG '
|
|
201 WRITE (ISYSWR,1200)
|
|
READ(ISYSRD,1101,ERR=201,END=999) IEXP
|
|
202 WRITE (ISYSWR,1201)
|
|
READ(ISYSRD,1107,ERR=202,END=999) U(1),WERR(1),A,B
|
|
IF (A) 212,210,212
|
|
210 IF (B) 212,211,212
|
|
211 LCODE(1) = 1
|
|
GO TO 220
|
|
212 UKP = U(1)+ABS(WERR(1))
|
|
UKM = U(1)-ABS(WERR(1))
|
|
IF (A.GE.UKM) GOTO 213
|
|
IF (B.LE.UKP) GOTO 213
|
|
LCODE(1) = 4
|
|
ALIM(1) = A
|
|
BLIM(1) = B
|
|
GOTO 220
|
|
213 WRITE (ISYSWR,1111)
|
|
GOTO 202
|
|
220 IF (WERR(1).EQ.0.) LCODE(1) = 0
|
|
NU = 3
|
|
U(6) = -1.
|
|
IF (IEXP.EQ.2) U(6) = 1.
|
|
IF (IEXP.EQ.4) U(6) = 1.
|
|
IF (IEXP.GE.2) PNAM(2)= 'GAMMA '
|
|
IF (IEXP.GE.4) PNAM(2)= 'NY '
|
|
IF (IEXP.GE.4) GOTO 230
|
|
NU = 4
|
|
221 WRITE (ISYSWR,1202)
|
|
READ(ISYSRD,1107,ERR=221,END=999) U(4),WERR(4)
|
|
IF (WERR(4).NE.0.) LCODE(4) = 1
|
|
230 NSUM = 0
|
|
SUMX = 0.
|
|
SUMY = 0.
|
|
SUMXX = 0.
|
|
SUMXY = 0.
|
|
DO 231 I = NXMIN,NXMAX
|
|
XR = U(6)*(XVAL(I)-U(1))/U(1)
|
|
IF (XR.LT.1E-5) GOTO 231
|
|
YR = YVAL(I) - U(4)
|
|
IF (YR.LT.1E-10) GOTO 231
|
|
XL = ALOG(XR)
|
|
YL = ALOG(YR)
|
|
SUMX = SUMX + XL
|
|
SUMY = SUMY + YL
|
|
SUMXX = SUMXX + XL*XL
|
|
SUMXY = SUMXY + XL*YL
|
|
NSUM = NSUM + 1
|
|
231 CONTINUE
|
|
IF (NSUM.GE.2) GOTO 232
|
|
WRITE (ISYSWR,1203)
|
|
GOTO 10
|
|
232 U(2) = (SUMXY-SUMX*SUMY/NSUM)/(SUMXX-SUMX*SUMX/NSUM)
|
|
WERR(2) = ABS(0.1*U(2))
|
|
LCODE(2) = 1
|
|
U(3) = EXP(SUMY/NSUM-U(2)*SUMX/NSUM)
|
|
WERR(3) = ABS(0.1*U(3))
|
|
LCODE(3) = 1
|
|
u(5)=0.
|
|
IF (IEXP.GT.0) GOTO 259
|
|
250 NU = 7
|
|
PNAM(4)= 'CRI.A '
|
|
PNAM(5)= 'CRI.C '
|
|
PNAM(6)= 'CRI.S '
|
|
PNAM(7)= 'BG '
|
|
U(7) = U(4)
|
|
WERR(7) = WERR(4)
|
|
LCODE(7) = LCODE(4)
|
|
U(6) = 3.
|
|
WERR(6) = 0.
|
|
LCODE(6) = 0
|
|
LCODE(5) = 1
|
|
LCODE(4) = 1
|
|
NSUM = 0
|
|
SUMX = 0.
|
|
SUMY = 0.
|
|
SUMXX = 0.
|
|
SUMXY = 0.
|
|
DO 255 I = NXMIN,NXMAX
|
|
XR = (XVAL(I)-U(1))/U(1)
|
|
IF (XR.GT.-1.E-4) GOTO 251
|
|
YR = YVAL(I) - U(3)*((-XR)**U(2)) - U(7)
|
|
XR = -XR*U(6)
|
|
GOTO 252
|
|
251 IF (XR.LT.1.E-4) GOTO 255
|
|
YR = YVAL(I) - U(7)
|
|
252 IF (YR.LT.1E-10) GOTO 255
|
|
XL = XR
|
|
YL = ALOG(YR)
|
|
SUMX = SUMX + XL
|
|
SUMY = SUMY + YL
|
|
SUMXX = SUMXX + XL*XL
|
|
SUMXY = SUMXY + XL*YL
|
|
NSUM = NSUM + 1
|
|
255 CONTINUE
|
|
IF (NSUM.GE.2) GOTO 256
|
|
PNAM(2)= 'CR.SCAT.'
|
|
WRITE (ISYSWR,1203) PNAM(2)
|
|
GOTO 10
|
|
256 U(4) = (SUMXY-SUMX*SUMY/NSUM)/(SUMXX-SUMX*SUMX/NSUM)
|
|
WERR(4) = ABS(0.1*U(4))
|
|
U(5) = EXP(SUMY/NSUM-U(4)*SUMX/NSUM)
|
|
WERR(5) = ABS(0.1*U(5))
|
|
259 CONTINUE
|
|
|
|
ELSEIF (IFIFU .EQ. 8) THEN
|
|
|
|
C . . . . . . PLOT ONLY
|
|
|
|
nu=0
|
|
RETURN
|
|
|
|
ELSE ! standard input
|
|
if (ififu.eq.7) then
|
|
if (usernp .eq. 0 .or. userfun .eq. 0) then
|
|
print *,'No user function implemented'
|
|
goto 999
|
|
endif
|
|
npuser=usernp
|
|
do k=1,npuser
|
|
pnam(k)=userpar(k)
|
|
psho(k)=usersho(k)
|
|
lcode(k)=1
|
|
enddo
|
|
if (nu .eq. 0) nu=1
|
|
nu_user=fit_userini(jfifu)
|
|
if (nu_user .gt. 0) npuser=nu_user
|
|
nu=0 ! nu is incremented for each check parameter below
|
|
else ! multi-voigt based functions
|
|
nu=0
|
|
if (narg .eq. 0) then
|
|
print *
|
|
print *,' Hint: In general, it is easier to use Single',
|
|
1 ' Gaussian'
|
|
print *,' for start and then add new peaks graphically or'
|
|
print *,' with command NEWPEAK'
|
|
endif
|
|
npeak=maxpeak
|
|
if (ififu .eq. 4) npeak=2
|
|
npuser=2+npeak*5
|
|
write(pnam, 291) 'Bg(Pos1)','dBg/dx ',
|
|
1 ('Posi. ',i
|
|
1 ,'MaxInt',i
|
|
1 ,'IntInt',i
|
|
1 ,'Fwhm G',i
|
|
1 ,'Fwhm L',i, i=1,npeak)
|
|
291 format (a8/a8,999(/a6,i2))
|
|
write(psho, 292) 'B','D',
|
|
1 ('P',i,'M',i,'I',i,'G',i,'L',i, i=1,npeak)
|
|
292 format (a1/a1,45(/a1,i1),999(5(/a1,i2)))
|
|
do k=1,npuser
|
|
lcode(k)=1
|
|
enddo
|
|
endif
|
|
C. . . . . . . . . . . . . . . . . . .
|
|
C. . . . . LOOP TO SET PARAMETERS . .
|
|
DO 590 K=1,npuser
|
|
521 if (nu_user .gt. 0) then
|
|
uk=u(k)
|
|
wk=0
|
|
fixed=.true.
|
|
a=0
|
|
b=0
|
|
elseif (funarg .lt. 0 .or. narg .eq. 0) then
|
|
if (k .eq. 1) WRITE (ISYSWR,1102)
|
|
if (ififu .eq. 2 .and. k .gt. 2 .and. mod(k,5) .eq. 2 .or.
|
|
1 ififu .eq. 4 .and. k .eq. 7) then ! do not ask for lorentzian HW
|
|
u(k)=0
|
|
werr(k)=0
|
|
lcode(k)=0
|
|
nu=nu+1
|
|
goto 590
|
|
endif
|
|
if (ififu .le. 4 .and. mod(k,5) .eq. 0) then ! Int.Int.
|
|
lcode(k)=0
|
|
nu=nu+1
|
|
goto 590
|
|
endif
|
|
522 WRITE (ISYSWR,1106) K,PNAM(K),': '
|
|
READ (ISYSRD,'(4F20.0)',ERR=522,END=529)
|
|
1 UK,WK,A,B
|
|
fixed=wk .eq. 0
|
|
IF (UK.EQ.999.) goto 529
|
|
else
|
|
if (k .gt. narg) goto 529
|
|
uk=parg(k)
|
|
wk=earg(k)
|
|
fixed=wk .eq. 0.0
|
|
a=0
|
|
b=0
|
|
endif
|
|
C. . . . . LIMIT CHECK . . . . . . .
|
|
IF (A .EQ. 0.0 .AND. B .EQ. 0.0) GOTO 540
|
|
UKP = UK+ABS(WK)
|
|
UKM = UK-ABS(WK)
|
|
if (a .ge. ukm .or. b .le. ukp) then
|
|
WRITE (ISYSWR,1111)
|
|
if (funarg .lt. 0 .or. narg .eq. 0) goto 521
|
|
a=ukm
|
|
b=ukp
|
|
endif
|
|
LCODE(K) = 4
|
|
ALIM(K) = A
|
|
BLIM(K) = B
|
|
C. . . . . FWHM 1 EQUAL TO ZERO . . . .
|
|
540 if (ififu .le. 4 .and. uk .eq. 0 .and.
|
|
1 k .gt. 1 .and. mod(k,5) .eq. 1) then ! gaussian hw
|
|
if (k .eq. 5) then
|
|
if (ififu .eq. 2) then
|
|
WRITE (ISYSWR,1112)
|
|
if (funarg .lt. 0 .or. narg .eq. 0) goto 521
|
|
uk=1
|
|
endif
|
|
goto 570
|
|
endif
|
|
C. . . . . FWHM CORRELATION . . .
|
|
560 ICSW(K) = 1
|
|
ICTO(K) = 5
|
|
CFAC(K) = 1.
|
|
WERRS(K) = WK
|
|
UK = U(5)
|
|
WK = 0.
|
|
WRITE(ISYSWR,1110) K,PNAM(K)
|
|
C. . . . . END FWHM COR . . . . .
|
|
endif
|
|
570 NU=NU+1
|
|
U(K) = UK
|
|
WERR(K) = WK
|
|
if (fixed) then
|
|
lcode(k)=0
|
|
werr(k)=0
|
|
endif
|
|
590 CONTINUE
|
|
goto 599
|
|
|
|
529 if (ififu .le. 4 .and. mod(nu,5) .ne. 2) then
|
|
nu=nu-mod(nu+2,5)
|
|
WRITE (ISYSWR, '(/X,A,/)')
|
|
1 'Not all parameters are defined -> last peak removed'
|
|
endif
|
|
|
|
599 continue
|
|
ENDIF
|
|
if (ififu .le. 4) then ! calculate int.intensities
|
|
ififu=1
|
|
do k=4,nu,5
|
|
YVALmax=voigt(0.0, u(k+2), u(k+3))
|
|
u(k+1)=u(k)/YVALmax
|
|
if (lcode(k) .ne. 1) then ! fixed or limited Max.Int
|
|
werr(k)=0
|
|
icsw(k+1)=-1
|
|
icto(k+1)=k
|
|
else
|
|
werr(k+1)=werr(k)/YVALmax
|
|
lcode(k)=0
|
|
lcode(k+1)=1
|
|
icsw(k)=-1
|
|
icto(k)=k+1
|
|
endif
|
|
enddo
|
|
endif
|
|
|
|
call extoin
|
|
call fit_check_range
|
|
|
|
if (ififu .eq. 8) return
|
|
if (isyswr .ne. 0) WRITE (ISYSWR,2000)
|
|
2000 FORMAT (/' First entry to FCN')
|
|
call fnctn(x, amin)
|
|
if (strange) then
|
|
strange=.false.
|
|
write(isyswr,'(/x,a/)')
|
|
1 '--- STRANGE peak: fit one gaussian'
|
|
call fit_fit(0)
|
|
write(isyswr,'(/x,a/)')
|
|
1 '--- STRANGE peak: calculate P and W from gaussian'
|
|
ififu=6
|
|
goto 10
|
|
endif
|
|
if (ififu .eq. 6) then
|
|
write(isyswr,'(/x,a/)')
|
|
1 '--- STRANGE peak: fit background'
|
|
call fit_fit(0)
|
|
return
|
|
endif
|
|
|
|
CALL FIT_PRINT (1)
|
|
WK=AMIN
|
|
NFCN=1
|
|
call fnctn(x, wk)
|
|
IF (WK.EQ.AMIN) RETURN
|
|
WRITE (ISYSWR,2001) AMIN,WK
|
|
2001 FORMAT (/' FCN is time-dependent: 1. Call F=',E16.7
|
|
1 /24X, '2. Call F=',E16.7)
|
|
999 RETURN
|
|
|
|
1100 FORMAT(
|
|
1/8X,'Single Gaussian = 0'
|
|
1/8X,'Single Voigtian = 1 (Lorentzian folded with Gaussian)',
|
|
1/8X,'Multi Gaussian = 2'
|
|
1/8X,'Multi Voigtian = 3'
|
|
1/8X,'Gaussian + Voigtian = 4'
|
|
1/8X,'Crit.Exponent = 5'
|
|
1/8X,'Strange = 6'
|
|
1/8X,A, ' = 7 (User function)'
|
|
1/8X,'Plot only = 8'//
|
|
1/8X,'Select function [0] : ',$)
|
|
|
|
1101 FORMAT (I1)
|
|
1102 FORMAT(/,' Parameter: Value,Error,(Lower-,Upper-Lim'
|
|
1,'it) end = 999',/)
|
|
1103 FORMAT (/,' Only ',I2,' variable parameters allowed',/)
|
|
1106 FORMAT (I6,4X,A8,A,$)
|
|
1107 FORMAT (4F10.0)
|
|
1110 FORMAT ('+',I5,4X,A8,' = Fwhm G 1')
|
|
1111 FORMAT (8X,'Input error: Parameter out of limits')
|
|
1112 FORMAT (8X,'Input error: Fwhm G 1 = 0.')
|
|
1200 FORMAT(/,' EXPONENT: BRAGG+CRIT.SCATTERING: = 0 ',
|
|
,/, ' ORD.PARAMETER: BETA,T<TC = 1 ',
|
|
,/, ' FLUCTUATIONS: GAMMA,T>TC = 2 ',
|
|
,/, ' GAMMA,T<TC = 3 ',
|
|
,/, ' CORREL.LENGTH: NY ,T>TC = 4 ',
|
|
,/, ' NY ,T<TC = 5 ',$)
|
|
1201 FORMAT(/,' ORDERING TEMPERATURE ',
|
|
, ' VALUE,ERROR,(LOWER-,UPPER-LIMIT) : ',$)
|
|
1202 FORMAT(/,' TEMP.INDEP.BACKGROUND VALUE,ERROR : ',$)
|
|
1203 FORMAT(/,' NOT ENOUGH DATAPOINTS TO DETERMINE ',A8)
|
|
END
|