subroutine fit_userfun(name, ufun)
!	----------------------------------

	include 'fit.inc'

	character*(*) name
	external ufun
	real ufun
	integer sys_adrr_rriii
	
	usernp=0
	npd=0
	usertit=name
	userfun = sys_adrr_rriii(ufun)
	end


	subroutine fit_userpar(name)
!	----------------------------

	include 'fit.inc'

	character*(*) name

	integer i

	if (npd .ne. 0)
	1  stop 'per-dataset parameters must be defined at end!'
	usernp=usernp+1
	if (usernp .gt. maxext) stop 'too many user parameters'
	i=index(name,':')
	userpar(usernp)=name(i+1:)
	if (i .gt. 1) then
	  call str_upcase(usersho(usernp),name(1:i-1))
	else
	  usersho(usernp)=' '
	endif
	usercct(usernp)=cinfo
	end


	subroutine fit_userpdp(name)
!	----------------------------

	include 'fit.inc'

	character*(*) name

	integer i

	npd=npd+1
	if (npd .gt. maxpd .or. npd+usernp .gt. maxext)
	1 stop 'too many per-dataset parameters'
	userpar(usernp+npd)=name
	do i=1,maxset
	  pdpar(npd, i)=none
	enddo
c	usercct(usernp+npd)=cinfo ! not yet implemented
	end


	subroutine fit_usercinfo(calcinfo)
!	----------------------------------

	include 'fit.inc'

	integer calcinfo

	cinfo=calcinfo
	end



	subroutine fit_ncurves(ncur)
!	----------------------------

	include 'fit.inc'

	integer ncur
	integer i

	nset=ncur
	if (nset .le. 0) then
	  do i=1,npkt
	    nset=max(nset,iset(i))
	  enddo
	endif
	end



	subroutine fit_user1st

	include 'fit.inc'

	real dum
	real sys_rfun_rriii

	integer calcinfo,i

	calcinfo=0
	do i=1,nu
	  if (u(i) .ne. pold(i)) then
	    calcinfo=ior(calcinfo, usercct(i))
	  endif
	enddo

	dum=sys_rfun_rriii(userfun, 0.0, u, nu, 1, calcinfo)
	end


	integer function fit_userini(jfifu)

	include 'fit.inc'

	integer i,jfifu
	real dum, xx

	real sys_rfun_rriii

	do i=1,nu
	  pold(i)=1.234e-35
	enddo
	xx=-float(jfifu)
	dum=sys_rfun_rriii(userfun, xx, u, nu, -1, 0)
	fit_userini=nint(xx)
	if (fit_userini .lt. 0) fit_userini=0
	end


	subroutine fit_check_range

	include 'fit.inc'

	real xmin, xmax
	real xmin0/-1./,xmax0/1.0/
	integer i

	if (ififu .ne. 7) RETURN
	xmin0=xval(nxmin)
	xmax0=xmin0
	do i=nxmin,nxmax
	  if (xval(i) .lt. xmin0) xmin0=xval(i)
	  if (xval(i) .gt. xmax0) xmax0=xval(i)
	enddo
	return

	entry fit_put_xrange(xmin, xmax)

	xmin0=min(xmin,xmax)
	xmax0=max(xmin,xmax)
	RETURN

	entry fit_limit_xrange(xmin, xmax)

	xmin0=max(xmin0,min(xmin,xmax))
	xmax0=min(xmax0,max(xmin,xmax))
	RETURN

	entry fit_get_xrange(xmin, xmax)

	xmin=xmin0
	xmax=xmax0
	end


	subroutine fit_fixit(i)

	include 'fit.inc'

	integer i

	lcode(i)=0
	end


	subroutine fit_delta(xd,sd)

	include 'fit.inc'

	real sd,xd
	logical read

	integer nd,nr
	real xdi(maxpeak), sdi(maxpeak)

	logical disable/.true./
	save nd,nr,xdi,sdi,disable

	if (sd .eq. 0 .or. disable) RETURN
	if (nd .ge. maxpeak) then
	  print *,'too many delta peaks'
	  RETURN
	endif
	nd=nd+1
	xdi(nd)=xd
	sdi(nd)=sd
	return

	entry fit_init_delta(read)
	if (read) then
	  nr=0
	else
	  nd=0
	endif
	disable=read
	return

	entry fit_get_delta(xd,sd)
	if (nr .ge. nd) then
	  sd=0
	  xd=0
	else
	  nr=nr+1
	  sd=sdi(nr)
	  xd=xdi(nr)
	endif
	end


	subroutine fit_confun(name, ufun)
!	---------------------------------

	include 'fit.inc'
	include 'fit_user.inc'
	
	character*(*) name
	external ufun, fit_cfun
	real ufun
	real fit_cfun
	integer sys_adrr_rriii
	
	data xbase/0.0/
	data show/.false./
	data ecnt/0/ 

	call fit_userfun(name, ufun)
	confun=userfun
	userfun = sys_adrr_rriii(fit_cfun)
	call fit_userpar('GW:GaussWid')
	call fit_userpar('XS:CalcStep')
	end


	subroutine fit_xbase(x_base)
!	----------------------------

	real x_base
	
	include 'fit_user.inc'
	
	xbase=x_base
	end


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

	implicit none

	real x		! x-value
	integer n	! number of parameters
	real p(*)	! parameters
	integer mode	! mode ! 0: x-dependent part, 1: x-independent part, -1: par. independent part
	integer cinfo	! calculation information (not used here)

	include 'fit_user.inc'

	integer maxt, maxg, ng, imaxg, m_f, m_i, m_t, m_x
	real pn, gn
	parameter (m_f=0, m_i=1, m_t=2, m_x=3)
! m_f:  call function directly
! m_i:  initialize function
! m_t:  transform x -> t
! m_x:  transform t -> x
	parameter (pn=2.0) ! step size of larger gaussian = fwhm/pn
	parameter (gn=2.0) ! total width = gn*fwhm
	parameter (imaxg=128) ! maximal fwhm (in steps)
	integer i,j,k,m

	real ts ! t-step
	real t ! transformed x
	real ta, tb  ! limits of t-space

	real gw ! total gaussian width (in t-space)
	real wg ! total peak weight
	real sum ! summing variable
	integer kstp ! step for folding
	integer mfold ! number of folding steps
	real wd ! individual peak width
	real t0 ! base for integration
	real lasty ! last value for integr.

	real xa, xb ! x-range
	real xs ! x-step at xbase
	real q ! ratio between fwhm's of peaks for folding steps
	real stp ! real step
	real xd, sd ! delta peak parameters

	real x0, x1, x2, qq, td, xsi

	parameter (maxt=10000) ! maximal number of points
	real y(0:maxt),yf(0:maxt) ! array before and after folding

	parameter (maxg=100) ! maximal number of points per gaussian (> gn*5)
	real gs(0:maxg) ! precalculated gaussian

	save y,ta,tb,xa,xb,gw,kstp

	external sys_rfun_rriii, voigt, spline3, spline4
	real     sys_rfun_rriii, voigt, spline3, spline4

	if (mode .eq. 0) then
	  if (gw .eq. 0) then
	    fit_cfun=sys_rfun_rriii(confun,x,p,n,m_f,cinfo)
	    RETURN
	  endif
	  if (x .lt. xa .or. x .gt. xb) then
	    fit_cfun=sys_rfun_rriii(confun,x,p,n,m_f,cinfo)
	    ecnt=ecnt+1
	    show=.true.
	  else
	    t=sys_rfun_rriii(confun,x,p,n,m_t,cinfo)/ts
	    t=(t-ta)/istp
	    i=nint(t-0.5)
	    if (i .le. 1) then
	      fit_cfun=spline3(t, y(0), y(1), y(2))
	    elseif (i .ge. npy-1) then
	      fit_cfun=spline3(npy-t, y(npy), y(npy-1), y(npy-2))
	    else
	      fit_cfun=spline4(t-i, y(i-1), y(i), y(i+1), y(i+2))
	    endif
	  endif
	else if (mode .gt. 0) then
	  show=.true.
	  call fit_init_delta(.false.)
	  if (p(1) .eq. 0.0) then
	    fit_cfun=sys_rfun_rriii(confun,x,p,n,m_f,cinfo)
	    RETURN
	  endif
	  call fit_get_xrange(xa,xb)
	  fit_cfun=sys_rfun_rriii(confun,x,p,n,m_i,cinfo)
	  call fit_init_delta(.true.)
	  ta=sys_rfun_rriii(confun,xa,p,n,m_t,cinfo)
	  tb=sys_rfun_rriii(confun,xb,p,n,m_t,cinfo)
	  xs=p(2)
	  if (xs .eq. 0) xs=p(1)/10
	  if (xbase-max(p(1),xs) .lt. xa .or.
	1     xbase+max(p(1),xs) .gt. xb) then
	    qq=(tb-ta)/(xb-xa)
	    ts=abs(xs*qq)
	    gw=abs(p(1)*qq)/ts
	  else
	    ts=abs(sys_rfun_rriii(confun,xbase+xs,p,n,m_t,cinfo)
	1         -sys_rfun_rriii(confun,xbase-xs,p,n,m_t,cinfo))
	    if (ts .eq. 0) ts=1.0
	    gw=abs(sys_rfun_rriii(confun,xbase+p(1)/2,p,n,m_t,cinfo)
	1         -sys_rfun_rriii(confun,xbase-p(1)/2,p,n,m_t,cinfo))/ts
	  endif
	  ta=int((ta/ts-gn*gw))
	  tb=int((tb/ts+gn*gw)+1)
	  npy=nint(tb-ta)
	  istp=max(int(gw/imaxg+1), (npy+maxt-1)/maxt)
	  npy=npy/istp

	  t0=ta-0.5*istp
	  x0=sys_rfun_rriii(confun,t0*ts,p,n,m_x,cinfo)
	  lasty=sys_rfun_rriii(confun,x0,p,n,m_f,cinfo)
	  k=0
	  do i=0,npy
	    x1=sys_rfun_rriii(confun,(t0+i*istp)*ts,p,n,m_x,cinfo)
	    ! simpson integral over interval (x0...x1)
	    xsi=(x1-x0)/istp
	    x2=x1-xsi/2
	    sum=lasty
	    do j=1,istp
	      lasty=sys_rfun_rriii(confun,x0+j*xsi,p,n,m_f,cinfo)
	      sum=sum+4*sys_rfun_rriii(confun,x2+j*xsi,p,n,m_f,cinfo)
	1            +2*lasty
	    enddo
	    x0=x1
	    y(i)=(sum-lasty)/6/istp
	    k=k+istp
	  enddo
	  call fit_get_delta(xd, sd)	  
	  do while (sd .ne. 0)
	    td=sys_rfun_rriii(confun,xd,p,n,m_t,cinfo)/ts
	    sd=sd/abs(
	1       sys_rfun_rriii(confun,(td+0.5*istp)*ts,p,n,m_x,cinfo)
	1      -sys_rfun_rriii(confun,(td-0.5*istp)*ts,p,n,m_x,cinfo) )
	    if (td .ge. ta .and. td .le. tb) then
	      t=(td-ta)/istp
	      i=nint(t-0.4999)
	      if (i .ge. npy) i=npy-1
	      y(i)=y(i)+sd*(i+1-t)
	      y(i+1)=y(i+1)+sd*(t-i)
	    endif
	    call fit_get_delta(xd, sd)
	  enddo
	  
	  if (gw .le. istp) RETURN

	  mfold=max(1,nint(log(gw/pn/istp)))
	  sum=0
	  q=max(1.1,(gw/pn/istp)**(1.0/mfold))
	  qq=1/(q*q)
	  qq=sqrt((1-qq**mfold)/(1-qq)) ! sqrt of sum of squared widths divided by gw
	  wd=gw/istp/q**(mfold-1)/qq

	  icnt=0
	  stp=1.0
	  do m=1,mfold
	    kstp=int(stp)
	    stp=stp*q
	    ng=min(maxg,int(gn*wd/kstp+0.9))
	    icnt=icnt+(2*ng+1)*(npy+1)-(ng-1)*ng+1
	    do i=0,ng
	      gs(i)=voigt(float(i), wd/kstp, 0.0)
	    enddo
	    do i=0,npy
	      wg=0
	      sum=0
	      k=i
	      do j=0,min(ng,i/kstp)
	        sum=sum+gs(j)*y(k)
	        wg=wg+gs(j)
	        k=k-kstp
	      enddo
	      k=i+kstp
	      do j=1,min(ng,(npy-i)/kstp)
	        sum=sum+gs(j)*y(k)
	        wg=wg+gs(j)
	        k=k+kstp
	      enddo
	      yf(i)=sum/wg
	    enddo
	    do i=ng*kstp,npy-ng*kstp
	      y(i)=yf(i)
	    enddo
	    wd=wd*q
	  enddo
	else
	  fit_cfun=sys_rfun_rriii(confun,x,p,n,mode,cinfo)
	  call fit_fixit(2)
	endif
	end


	subroutine fit_confun_print

	include 'fit_user.inc'

	if (show) then
	  print *
	  print '(i10,a,i9,a)'
	1   , npy*istp*2+1, ' fn. calls, fold. multiplications: ',icnt
	  if (ecnt .gt. 0) then
	    print *,ecnt,' outside range calls counted'
	    ecnt=0
	  endif
	  show=.false.
	endif
	end


	real function spline4(x, y1, y2, y3, y4)

! spline interpolation.
! the result is valid for values x between 0 and 1 (x-coord of y2 and y3)

	real x, y1, y2, y3, y4
	real ss, sl, sr

	ss=y3-y2                 ! median slope
	sl=(y3-y1)/2             ! slope at point 2
	sr=(y4-y2)/2             ! slope at point 3
	spline4=y2+x*ss          ! linear interpol.
	1 + (sl-ss)*x*(x-1)**2   ! correction for left slope
	1 - (sr-ss)*(1-x)*x**2   ! correction for right slope
	end


	real function spline3(x, y1, y2, y3)

! spline interpolation.
! the result is valid for values x between 0 and 1 (x-coord of y1 and y2)
! or, as an extrapolation to the left of x1

	real x, y1, y2, y3
	real sl, sr, s1

	s1=y2-y1
	sr=(y3-y1)/2                 ! interpolate slope at point 2
	sl=2*sr-s1                   ! extrapolate slope at point 1
	spline3=y1+x*s1              ! linear interpol.
	1 + (sl-s1)*x*(x-1)**2       ! correction for left slope
	1 - (sr-s1)*(1-x)*x**2       ! correction for right slope
	end