program BOSE
!	------------
!
! Simple user function example (straight line).
!
	implicit none
	external BOSE_FUN

	character str*32
	integer i,l

!---
! Welcome message

	call fit_vers(str)
	call str_trim(str, str, l)

	type '(X,2A)','Program FIT(BOSE)     Version ',str(1:l)
	do i=1,l
	  str(i:i)='-'
	enddo
	type '(X,2A/)','-----------------------------',str(1:l)

!---
! Function title and parameter names
!
	call fit_confun('lorenz/bose', bose_fun)	! function title, function
	call fit_userpar('BG:Bgr(0)')
	call fit_userpar('DB:dBgr/dX')
	call fit_userpar('EF:EF')
	call fit_userpar('T:Temp')
	do i=1,9
	  write(str, '(2(a,i1))') 'P',i,':Pos',i
	  call fit_userpar(str)
	  write(str, '(2(a,i1))') 'I',i,':Int',i
	  call fit_userpar(str)
	  write(str, '(2(a,i1))') 'W',i,':Wid',i
	  call fit_userpar(str)
	enddo
	call fit_main
	end



	real function bose_fun(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 i,i0,j,k
	parameter (i0=7)
	real x0,w0,y0,db,bg,kf,l0
	real voigt, bose_fact
	real xnew(9),ynew(9),wnew(9)

	if (mode .eq. 0) then

	  bose_fun=0
	  do i=i0,n-2,3
	    if (p(i+2) .ne. 0) then ! ignore delta functions (treated later)
	      bose_fun=bose_fun+p(i+1)*voigt(x-p(i), 0.0, p(i+2))
	      if (p(i+2) .lt. 0) then ! make a mirror peak for negative width
	        x0=x+p(i)
	        if (i .gt. i0 .and. p(i0+2) .eq. 0.0) x0=x0-2*p(i0) ! shift zero
	        bose_fun=bose_fun+p(i+1)*voigt(x0, 0.0, -p(i+2))
	      endif
	    endif
	  enddo
	  bose_fun=bose_fun*bose_fact(x/p(6))+p(3)+p(4)*x

	elseif (mode .eq. 1) then

! x-independent part

	  do i=i0,n-2,3
	    if (p(i+2) .eq. 0) then ! treat delta functions
	      call fit_delta(p(i),p(i+1))
	    endif
	  enddo
	  call fit_limit_xrange(-p(5),1e6)

	elseif (mode .eq. 2) then ! transform x -> t

	  x0=x+2*p(5)
	  if (x0 .ge. 0) then
	    bose_fun=sqrt(x0)
	  else
	    bose_fun=0
	  endif

	elseif (mode .eq. 3) then ! transform t -> x

	  bose_fun=x*x-2*p(5)

	else
	  if (nint(x) .eq. -1 .and. n .ge. 7) then ! convert from multi-voigt
	    print *
	    if (n .eq. 7) then
	      print *,'Convert from voigt'
	    else
	      print *,'Convert from multi-voigt'
	    endif
	    db=p(2)
	    bg=p(1)-p(3)*db ! different bg definition
	    j=3
	    x0=p(3)
	    y0=p(5)
	    w0=max(abs(p(6)), abs(p(7)))
	    l0=p(7)
	    do i=8,n,5
	      if (abs(p(i))  .lt. abs(x0)) then
	        j=i
	        x0=p(i)
	        y0=p(i+2)
	        w0=max(abs(p(i+3)), abs(p(i+4)))
	      endif
	    enddo
	    k=0
	    do i=3,n,5
	      if (p(i) .gt. 0 .and. i .ne. j) then
	        k=k+1
	        xnew(k)=p(i)
	        ynew(k)=p(i+2)
	        wnew(k)=-max(abs(p(i+3)), abs(p(i+4)))
	      endif
	    enddo
	    p(1)=w0
	    p(2)=w0*0.05
	    p(3)=bg
	    p(4)=db
	    kf=1.55
	    call fit_get_real('KFIX', kf)
	    p(5)=2.0723*kf*kf
	    p(6)=10   ! default Temp
	    call fit_get_real('Temp', p(6))
	    p(7)=x0
	    p(8)=y0
	    p(9)=l0
	    i=10
	    do j=1,k
	      p(i)=xnew(k)
	      p(i+1)=ynew(k)/bose_fact(p(i)/p(6))
	      p(i+2)=wnew(k)
	      i=i+3
	    enddo
	    x=i-1
	  else
	    print *
	    print *,'Up to 9 Lorenzians multiplied with bose_factor'
	1          ,', folded with gaussian'

	  endif
	  print *,'Negative Wid makes a mirror peak at -Pos'
	  print *
	endif
	end


        real function bose_fact(wot)

! argument: omega[meV] over T[K]
        real wot
	real K_meV, x
	parameter (K_meV=11.6048)

	x=wot*K_meV
        if (abs(x) .lt. 1e-5) then
          bose_fact=1.0
        else if (x .lt. -30) then
          bose_fact=0
        else
          bose_fact=x/(1-exp(-x))
        endif
        end