; $Id: ref_diffusion_in_a_sphere.pro,v 1.8 2006/07/26 14:41:35 azuah Exp $
;###############################################################################
;
; NAME:
;  REF_DIFFUSION_IN_A_SPHERE
;
; PURPOSE:
;  Implementation of a diffusion-in-a-sphere 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_diffusion_in_a_sphere,xvals,parms,					$
                    	parmnames = parmnames,					$
                    	init_guess = init_guess,				$
                    	qty = qty,								$
                    	single_parmnames = single_parmnames,	$
                    	multi_parmnames = multi_parmnames,		$
                    	res_ptr = res_ptr,						$
                    	all_y = all_y,							$
                    	step = step,							$
                    	xunits = xunits,						$
                    	_Extra = extra
; all_y keyword is not 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 =	['Sphere radius','Diffusion constant', $
			'Intensity:'+strtrim(string(indgen(ny)),2), $
			'Center:'+strtrim(string(indgen(ny)),2)]

single_parmnames = parmnames[0:1]
multi_parmnames = ['Intensity','Center']
if arg_present(step) then begin
	step = dblarr(2+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
		; Guess the 1/2 width at the highest y-value
		; Get the low value first
		amp = max(qty[*,qty_size[2]-1])
		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]]
		fwhm = x_hi_lim - x_lo_lim
		q_value = y[qty_size[2]-2]
		hbar = 0.6528
		if strupcase(xunits) eq 'MEV' then hbar = hbar*1.d4 else hbar = hbar*1.d7
		d = 0.5*fwhm/(hbar*q_value^2)
		a = 2.0 & center = x[index_array]
		init_guess = [a,d,intensity,center]
	endif
	return,-1
endif
ret = extra.otable->get_xnl(xnl = xnl, l = l)

a = parms[0]
D = parms[1]
hbar = 0.6528
if strupcase(xunits) eq 'MEV' then hbar = hbar*1.d4 else hbar = hbar*1.d7;D = D*1.d4 else D = D*1.d7
ntable = n_elements(xnl)
intensity = ux#(parms[2:2+ny-1])
center = ux#(parms[2+ny:2+2*ny-1])
center1 = parms[2+ny:2+2*ny-1]  ; RTA
resolution = *res_ptr

; Now calculate the EISF
ao = ux#(((3.0*sph_bessel(y*a,1))/y*a)^2)
sum = dblarr(nx,ny)
for i = 1,ntable-1 do begin

	anl = 1.5*uy##(((sph_bessel(xnl[i],1))^2)*((xnl[i])^2-l[i]*(l[i]+1.0))/(xnl[i])^2)
	ok = where(y*a ne xnl[i],count_ok)

	if count_ok gt 0 then begin
		first_term = (6.0*xnl[i]^2)/((xnl[i])^2-l[i]*(l[i]+1.0))
		second_term_num = (y[ok]*a*sph_bessel(y[ok]*a,l[i]+1)-l[i]*sph_bessel(y[ok]*a,l[i]))
		second_term_den = (y[ok]*a)^2-(xnl[i])^2
		second_term = second_term_num/second_term_den
		anl[ok] = (second_term)
	endif

	num_arg = ux#((xnl[i]^2)*hbar*D/a^2)#uy
	sum = sum + (2.0*l[i]+1)*(ux#anl)*num_arg/$
		(num_arg^2+(x#uy-center)^2)
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)
lineshape = ao*ref_delta(xvals,center1,resolution,nx = nx,ny = ny)+$ ; RTA - center->center1
            (1.0/!dpi)*scat_convolute(x,y,sum,resolution) ;RTA - scat_convolve->scat_convolute
yout = intensity*lineshape
return,yout
end
