program ufit				! change FIT to your own program name
!	------------
!
! Simple user function example (straight line).
!
	implicit none
	external FIT_ufun			! change FIT_ufun to your own function name

!---
! Welcome message

	print '(5(/X,A))'
	1,'Program UFIT'
	1,'------------'
	1,'User function: sum of lorentzian folded with meas. resolution'

!---
! Function title and parameter names
!
	call fit_userfun('Quasielastic', fit_ufun)	! function title, function
	call fit_userpar('B:Bg(0)')     ! first parameter: background at zero
	call fit_userpar('D:dBg/dX')    ! second parameter: slope
	call fit_userpar('S:bg.scale')  ! background slope
	call fit_userpar('G:fwhm gaussian')
	call fit_userpar('P:Pos')       ! position
	call fit_userpar('I1:IntInt 1') ! 1st lorentzian intensity
	call fit_userpar('L1:fwhm L 1') ! 1st lorentzian width
	call fit_userpar('I2:IntInt 2') ! 2nd lorentzian intensity
	call fit_userpar('L2:fwhm L 2') ! 2nd lorentzian width
	call fit_userpar('I3:IntInt 3') ! 3rd lorentzian intensity
	call fit_userpar('L3:fwhm L 3') ! 3rd lorentzian width
	call fit_main
	end


	real function fit_ufun(x,p,n,mode,cinfo)
!	-------------------------------------------

	implicit none

	real x		! x-value
	integer n	! number of parameters
	real p(n)	! parameters
	integer mode	! mode
	integer cinfo	! calculation information (see below)

	integer npnt
	parameter (npnt=10000)
	real xx(npnt), yy(npnt)
	real gg,xp,b,q
	integer idx/1/, nb/0/

	real voigt

	if (mode .eq. 0) then

! Define here your own function

	  xp=x-p(5)
	  gg=p(4)
	  fit_ufun=p(1)+xp*p(2)
	1         +p(6)*voigt(xp, gg, p(7))
	1         +p(8)*voigt(xp, gg, p(9))
	1         +p(10)*voigt(xp, gg, p(11))
	  if (idx .le. 0 .or. idx .ge. nb) stop 'FIT_UFUN: illegal IDX'
10	  if (x .gt. xx(idx+1)) then
	    if (idx .lt. nb-1) then
	      idx=idx+1
	      goto 10
	    endif
	  else
20	    if (x .lt. xx(idx)) then
	      if (idx .gt. 1) then
	        idx=idx-1
	        goto 20
	      endif
	    endif
	  endif
	  q=(x-xx(idx))/(xx(idx+1)-xx(idx))
	  b=yy(idx)*(1-q)+q*yy(idx+1)
	  fit_ufun=fit_ufun+p(3)*b

	elseif (mode .lt. 0) then

	  call fit_sort(0,0) ! sort data
	  call fit_get_array('X', xx, npnt, nb)
	  if (nb .ge. npnt) print *,'background points limit reached:',npnt
	  idx=1
	  call fit_get_array('Y', yy, npnt, nb)

	endif
	end