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