; $Id: ref_random_jump.pro,v 1.5 2007/03/20 16:03:38 azuah Exp $
;###############################################################################
;
; NAME:
;  REF_RANDOM_JUMP
;
; PURPOSE:
;  Implementation of the random 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_random_jump,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,                         $
                     lorentzian = lorentzian,               $
                     step = step,                           $
                     eisf = eisf,                           $
                     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 = [  'RMS vibrational displacement (A)',          $
               'bond length (A)',                           $
               'translational diffusivity (A^2/ps)',        $
               'site residence time (ps)',                  $
               'rotational correlation time (ps)',          $
               'Intensity:'+strtrim(string(indgen(ny)),2),  $
               'Center:'+strtrim(string(indgen(ny)),2),     $
               'EISF: '+strtrim(string(indgen(ny)),2)       ]
nsingle = 5
nmulti = 3
single_parmnames = parmnames[0:nsingle-1]
multi_parmnames = ['Intensity','Center','EISF']

if strupcase(xunits) eq 'MEV' then $
   hbar = 0.6594 else $
   hbar = 659.4

if arg_present(step) then begin
   step = dblarr(nsingle+3*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])
      center = x[index_array]
      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
      d = fwhm/(2.0*hbar*(y[ny-1]^2))
      tauo = 0.25 ; according to Sankar Nair, this is the approximate
                  ; residence time for bulk water at room temperature.
      u = 0.0
      tau = 0.1
      b = 1.0     ; initial guess for bond length
      center = x[index_array]
      eisf = 0.5+dblarr(ny)
      init_guess = [u,b,d,tauo,tau,intensity,center,eisf]
   endif
   return,-1
endif

u = parms[0]
b = parms[1]
d = parms[2]
tauo = parms[3]
tau = parms[4]
intensity = ux#(parms[nsingle:nsingle+ny-1])
center = ux#(parms[nsingle+ny:nsingle+2*ny-1])
center1 = parms[nsingle+ny:nsingle+2*ny-1] ; RTA
eisf = ux#(parms[nsingle+2*ny:nsingle+3*ny-1])

jo = ux#(sph_bessel(y*b,0))
j1 = ux#(sph_bessel(y*b,1))
j2 = ux#(sph_bessel(y*b,2))
gamma = d*y^2/(1.0+d*tauo*y^2)
gamma = ux#gamma
dwf_term = ux#(exp(-(y*u)^2/3.0))


x = x#uy
x = x-center

term1 = (jo^2)*hbar*gamma/((hbar*gamma)^2+x^2)
term2 = 3.0*(j1^2)*hbar*(gamma+1.0/(3.0*tau))/((hbar*(gamma+1.0/(3.0*tau)))^2+x^2)
term3 = 5.0*(j2^2)*hbar*(gamma+1.0/(1.0*tau))/((hbar*(gamma+1.0/(1.0*tau)))^2+x^2)
sqe_rjd = (1.0/!dpi)*(term1+term2+term3)

resolution = *res_ptr
delta_function = ref_delta(xvals,center1,resolution,nx = nx,ny = ny) ; RTA - center->center1

;lineshape = scat_convolute(x,y,sqe_rjd,resolution) ; RTA - ;scat_convolve->scat_convolute
lineshape = scat_convolute(xvals[0:nx-1],y,sqe_rjd,resolution) ; RTA - scat_convolve->scat_convolute
;yout = (intensity*dwf_term)*(eisf*lineshape+(1.0-eisf)*delta_function)
yout = (intensity*dwf_term)*(eisf*delta_function + (1.0-eisf)*lineshape) ; RTA - change based on PDF doc for func defn
return,yout
end
