233 lines
4.1 KiB
Fortran
233 lines
4.1 KiB
Fortran
real function FIFUN(XX)
|
|
C -----------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
real xx
|
|
|
|
real fifunp
|
|
|
|
fifun=fifunp(xx, u)
|
|
end
|
|
|
|
|
|
subroutine FNCTN(PINT,F)
|
|
c ------------------------
|
|
include 'fit.inc'
|
|
c
|
|
real pint(*),f
|
|
|
|
integer i,j
|
|
real sum
|
|
|
|
real fifun ! function
|
|
|
|
call intoex(pint)
|
|
nfcn=nfcn+1
|
|
|
|
if (ififu .eq. 6) then
|
|
j=0
|
|
sum = 0.
|
|
DO I=NXMIN,NXMAX
|
|
if (abs(xval(i)-u(3)) .gt. u(4)/2) then
|
|
sum = sum + ((YVAL(I)-FIFUN(XVAL(I)))/SIG(I))**2
|
|
j=j+1
|
|
endif
|
|
ENDDO
|
|
if (j .lt. 2) then
|
|
sum = sum + ((YVAL(nxmin)-FIFUN(XVAL(nxmin)))/SIG(nxmin))**2
|
|
1 + ((YVAL(nxmax)-FIFUN(XVAL(nxmax)))/SIG(nxmax))**2
|
|
j=j+2
|
|
endif
|
|
nfree=max(1,j-npar)
|
|
F = sum/nfree
|
|
else
|
|
if (ififu .eq. 7) call fit_user1st
|
|
sum = 0.
|
|
DO I=NXMIN,NXMAX
|
|
actset=iset(i)
|
|
sum = sum + ((YVAL(I)-FIFUN(XVAL(I)))/SIG(I))**2
|
|
ENDDO
|
|
F = sum/nfree
|
|
endif
|
|
END
|
|
|
|
|
|
real FUNCTION FIFUNP(XX, P)
|
|
C ---------------------------
|
|
|
|
include 'fit.inc'
|
|
|
|
real xx, p(maxext)
|
|
|
|
real xr
|
|
integer i,j
|
|
real sys_rfun_rriii, voigt
|
|
|
|
goto (100,100,100,100,500,600,700) ififu
|
|
fifunp=0
|
|
return
|
|
|
|
C. . . . . VOIGT . . .
|
|
|
|
100 continue
|
|
if (nu .eq. 2) then
|
|
fifunp = p(1)+p(2)*xx
|
|
return
|
|
endif
|
|
fifunp = p(1)+p(2)*(xx-p(3))
|
|
do i=3,nu,5
|
|
fifunp=fifunp+p(i+2)*voigt(xx-p(i), p(i+3), p(i+4))
|
|
enddo
|
|
return
|
|
|
|
C. . . . . CRITICAL EXPONENT. . .
|
|
|
|
500 IF (P(5).EQ.0.) THEN
|
|
XR = P(6)*(XX-P(1))/P(1)
|
|
IF (XR.LT.1E-6) then
|
|
fifunp=p(4)
|
|
else
|
|
FIFUNP = P(3)*(XR**P(2)) + P(4)
|
|
endif
|
|
RETURN
|
|
ENDIF
|
|
XR = (XX-P(1))/P(1)
|
|
IF (XR.LE.-1.E-10) THEN
|
|
FIFUNP = P(3)*((-XR)**P(2)) + P(5)*EXP((-XR*P(6))*P(4)) + P(7)
|
|
RETURN
|
|
ENDIF
|
|
FIFUNP = P(5)*EXP(XR*P(4)) + P(7)
|
|
RETURN
|
|
|
|
C. . . . . STRANGE . . . . . . . .
|
|
|
|
600 fifunp = p(1)+p(2)*(xx-p(3))
|
|
return
|
|
|
|
C. . . . . SPECIAL FUNCTION . . .
|
|
|
|
700 if (nu .gt. 0) then
|
|
if (npd .eq. 0) then
|
|
fifunp=sys_rfun_rriii(userfun, xx, p, nu, 0, actset)
|
|
else
|
|
j=nu
|
|
do i=1,npd
|
|
j=j+1
|
|
p(j)=pdpar(i,min(actset,maxset))
|
|
enddo
|
|
fifunp=sys_rfun_rriii(userfun, xx, p, nu+npd, 0, actset)
|
|
endif
|
|
else
|
|
fifunp=0
|
|
endif
|
|
return
|
|
end
|
|
|
|
|
|
real function voigt(x, gg, gl)
|
|
|
|
implicit none
|
|
|
|
real x,gg,gl
|
|
|
|
real vint, sqrtln2, pi
|
|
parameter (vint=0.939437278699651, sqrtln2=0.832554611157697)
|
|
parameter (pi=3.14159265)
|
|
|
|
real xg, c
|
|
|
|
complex cwerf
|
|
|
|
if (gl .eq. 0) then ! pure gaussian
|
|
if (gg .eq. 0) then
|
|
if (x .eq. 0) then
|
|
voigt=1e18
|
|
else
|
|
voigt=0
|
|
endif
|
|
else
|
|
xg=(x/gg*(sqrtln2*2))**2
|
|
if (xg .lt. 50.) then
|
|
voigt=exp(-xg)*vint/abs(gg)
|
|
else
|
|
voigt=0
|
|
endif
|
|
endif
|
|
elseif (abs(gl) .gt. abs(gg)*100) then ! pure lorentzian
|
|
voigt=abs(gl)/(gl*gl/4+x*x)/(2*pi)
|
|
elseif (abs(x) .gt. 12*abs(gg) .and.
|
|
1 (abs(gg) .gt. abs(gl) .or.
|
|
1 abs(x/gl) .gt. 12*sqrt(abs(gg/gl)))) then ! approx. for high x
|
|
voigt=abs(gl)/(gl*gl/4+x*x)/(2*pi)/(1-0.541*(gg/x)**2)
|
|
else
|
|
c=sqrtln2/abs(gg)
|
|
voigt=real(cwerf(2*x*c, abs(gl)*c))*vint/abs(gg)
|
|
endif
|
|
end
|
|
|
|
|
|
COMPLEX FUNCTION CWERF(xx,yy)
|
|
|
|
COMPLEX VALUE
|
|
REAL LAMBDA
|
|
real const
|
|
parameter (CONST=1.12837916709551)
|
|
|
|
real xx,yy,x,y,h,s,h2,r1,r2,s1,s2,fn,c,t1,t2
|
|
integer nc,nu,n
|
|
|
|
X=ABS(XX)
|
|
Y=ABS(YY)
|
|
IF (Y .GE. 4.29 .OR. X .GE. 5.33) THEN
|
|
H=0.
|
|
NC=0
|
|
NU=9
|
|
LAMBDA=0.
|
|
ELSE
|
|
S=(1.0-Y/4.29)*SQRT(1.0-X*X/28.41)
|
|
H=1.6*S
|
|
H2=2.0*H
|
|
NC=6+INT(23.0*S)
|
|
NU=10+INT(21.0*S)
|
|
LAMBDA=H2**NC
|
|
NC=NC+1
|
|
IF (LAMBDA .EQ. 0.0) NC=0
|
|
END IF
|
|
R1=0.
|
|
R2=0.
|
|
S1=0.
|
|
S2=0.
|
|
DO N=NU,1,-1
|
|
FN=N
|
|
T1=Y+H+FN*R1
|
|
T2=X-FN*R2
|
|
C=0.5/(T1*T1+T2*T2)
|
|
R1=C*T1
|
|
R2=C*T2
|
|
IF (N .LE. NC) THEN
|
|
T1=LAMBDA+S1
|
|
S1=R1*T1-R2*S2
|
|
S2=R2*T1+R1*S2
|
|
LAMBDA=LAMBDA/H2
|
|
ENDIF
|
|
ENDDO
|
|
IF (NC .EQ. 0) THEN
|
|
S1=R1
|
|
S2=R2
|
|
ENDIF
|
|
S1=CONST*S1
|
|
c IF (Y .EQ. 0.0) RS1=EXP(-X*X)
|
|
VALUE=CMPLX(S1,CONST *S2)
|
|
IF(YY .GE. 0.0) THEN
|
|
IF(XX .LT. 0.0) VALUE=CONJG(VALUE)
|
|
ELSE
|
|
value=-value
|
|
VALUE=2.0*CEXP(-CMPLX(X,Y)**2)-VALUE
|
|
ccc print *,xx,yy,value
|
|
IF(XX .GT. 0.0) VALUE=CONJG(VALUE)
|
|
ENDIF
|
|
CWERF=VALUE
|
|
RETURN
|
|
END
|