; $Id: ref_bounded_jump_diffusion.pro,v 1.3 2006/07/26 14:41:35 azuah Exp $
;###############################################################################
;
; NAME:
;  REF_BOUNDED_JUMP_DIFFUSION
;
; PURPOSE:
;  Implementation of the bounded jump diffusion model for use in fitting
;  data in RAINS.
;
; CATEGORY:
;  DAVE, Data Analysis, RAINS, surface fitting
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ref_bounded_jump_diffusion,xvals,parms,				$
                    	parmnames = parmnames,					$
                    	init_guess = init_guess,				$
                    	qty = qty,								$
                    	step = step,							$
                    	single_parmnames = single_parmnames,	$
                    	multi_parmnames = multi_parmnames,		$
                    	res_ptr = res_ptr,						$
                    	all_y = all_y,							$
                    	xunits = xunits,						$
                    	_Extra = extra
; all_y keyword is relevant for this model
if n_elements(all_y) eq 0 then all_y = 1
nx = extra.nx & ny = extra.ny
x = xvals[0:nx-1] & y = xvals[nx:nx+ny-1]
ux = 1+bytarr(nx) & uy = 1+bytarr(ny)

parmnames =	['Confinement distance','Jump distance', $
			'HWHM: ',$
			'Intensity:'+strtrim(string(indgen(ny)),2), $
			'Center:'+strtrim(string(indgen(ny)),2)]

single_parmnames = parmnames[0:2]
multi_parmnames = ['Intensity','Center']
if arg_present(step) then begin
	step = dblarr(3+2*extra.ny)
endif
if (n_params() eq 0) or (n_params() eq 1) then begin
	if arg_present(init_guess) then begin
		intensity = max(qty,dimension = 1,index)
		qty_size = size(qty)
		for i = 0,qty_size[2]-1 do begin
			qty_max = max(qty[*,i],index)
			if i eq 0 then index_array = index else $
				index_array = [index_array,index]
		endfor
		center = x[index_array]
		; Guess the 1/2 width at the highest y-value
		; Get the low value first
		amp = max(qty[*,qty_size[2]-2])
		ht = 0.5*amp
		wh_lo = where(abs(amp-ht) le min(abs(amp-ht)) and (x lt 0.0),count_lo)
		x_lo_lim = x[wh_lo[0]]
 		wh_hi = where(abs(amp-ht) le min(abs(amp-ht)) and (x gt 0.0),count_hi)
		x_hi_lim = x[wh_hi[0]]
		hwhm = 0.5*(x_hi_lim - x_lo_lim)
		ro = 1.0
		l = 2.0
		init_guess = [l,ro,hwhm,intensity,center]
	endif
  return,-1
endif

l = parms[0]
ro = parms[1]
hwhm = parms[2]

intensity = ux#(parms[2:2+ny-1])
center = ux#(parms[2+ny:2+2*ny-1])

ao = ux#((sph_bessel(0.5*y*l,0))^2)
hbar = 0.6528	; ueV.ns or meV.ps
nterms = 40
sum = dblarr(nx,ny)
for i = 1,nterms-1 do begin
	ai = ux#(((2.0*y*l)^2/((y*l)^2-(i*!dpi)^2)^2)*(1.0-((-1.0)^i)*cos(y*l)))
	width = ux#(hwhm*(1.0-exp(-0.5*(i*!dpi*ro/l)^2)))#uy
	lor = width/(width^2+(x#uy-center)^2)
	sum = sum + ai*lor
endfor

xhi = max(x,min = xlo)
dx = x[1]-x[0] & dy = y[1]-y[0]
nshift = fix(0.5*((abs(xhi)-abs(xlo)))/dx)
resolution = *res_ptr
center = parms[2+ny:2+2*ny-1]   ; RTA
lineshape = ao*ref_delta(xvals,center,resolution,nx = nx,ny = ny)+$
            (1.0/!dpi)*scat_convolute(x,y,sum,resolution) ;RTA change scat_convolve->scat_convolute
yout = intensity*lineshape

return,yout
end
