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