Commit 4bd8f112 authored by Gonzalez, Miguel's avatar Gonzalez, Miguel
Browse files

New QENS model

parent 531a0737
......@@ -6168,6 +6168,94 @@ function QENS_Fixed_JumpsNsitesLogNormDist, X, P1, P2, P3, P4, P5, P6, P7, P8
end
function QENS_Mix1, X, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12
;+
;******************************************************************************
;** Delta +
;** Jump Translation Diffusion +
;** Localized Diffusion (Gaussian model from Volino) +
;** Isotropic Rotation
;**
;** This assumes 4 independent populations, each contributing to
;** a single type of motion. However there is no way to ensure
;** that P3+P5+P8 < 1, so the fit can go to non-physical values
;** such that the fraction for the last component 1-P3-P5-P8 < 0.
;**
;** As the position of the delta can be between two energy values,
;** the intensity of the delta is shared into 2 deltas placed at those two energies.
;**
;** Parameters:
;** P1 = Scaling factor
;** P2 = Center
;** P3 = Fraction of atoms contributing to delta
;** P4 = u^2 for DW term
;** P5 = Fraction of atoms contributing to 1st lorentzian (jump translational diff)
;** P6 = D (self-diffusion coefficient in A^2*meV) for 1st lorentzian
;** P7 = Residence time (meV^-1) for 1st lorentzian
;** P8 = Fraction of atoms contributing to 2nd lorentzian (Volino's model for localized diff)
;** P9 = D (self-diffusion coefficient in A^2*meV) for 2nd lorentzian
;** P10 = Variance <ux^2> in Volino's model for 2nd lorentzian
;** P11 = Radius of the sphere for the isotropic rotation for 3rd lorentzian
;** P12 = Rotational diffusion constant (meV) for 3rd lorentzian
;-
common c_str_fitqens, func_resol, Qcurrent1D, Qcurrent2D, mig2, mig3, mig4, mig5
nmax1 = 100 ;number of lorentzians to use in sum of localized diff part (theoretically is infinite)
nmax2 = 10 ;number of lorentzians to use in sum of rotational part (theoretically is infinite)
Y = QENS_delta(X, 1.0, P2, P3, 0.0) ;delta (area=1, position, fraction=P3, DW=0.0)
;first jump diffusion lorentzian
hwhm1 = P6 * Qcurrent1D^2 / (1.0 + P7 * P6 * Qcurrent1D^2)
Y += P5 * strnormLorentz(X, P2, hwhm1)
;second localized lorentzian
al = fltarr(nmax1+1) & hwhm2 = fltarr(nmax1+1)
arg = double(Qcurrent1D^2 * P10)
exparg = EXP((-arg)>(-70.))
;elastic term
if arg gt 0.0 then al[0] = exparg else al[0]=1.0
;quasielastic terms
for i = 1, nmax1 do begin
if arg gt 0.0 then begin
al[i] = exparg*(arg^i)/factorial(i)
hwhm2[i] = i*P9/P10
endif else al[i]=0.0
endfor
;elastic term (EISF)
Y += P8 * al[0] * QENS_delta(X, 1.0, P2, 1.0, 0.0) ;delta (area=1, position, fraction=1, DW=0.0)
;quasielastic terms
for i=1, nmax1 do begin
Y += P8 * al[i] * strnormLorentz(X, P2, hwhm2[i])
endfor
;third isotropic rotational lorentzian
jl = fltarr(nmax2+1) & hwhm3 = fltarr(nmax2+1)
arg = Qcurrent1D * P11
;elastic term
if arg gt 0.0 then jl[0] = sqrt(!pi/2/arg)*beselj(arg,0.5) else jl[0]=1.0
;quasielastic terms
for i = 1, nmax2 do begin
if arg gt 0.0 then jl[i] = sqrt(!pi/2/arg)*beselj(arg,i+0.5) else jl[i]=0.0
hwhm3[i] = i * (i+1) * P12
endfor
jl2 = jl^2
;elastic term (EISF)
Y += (1-P3-P5-P8) * jl2[0] * QENS_delta(X, 1.0, P2, 1.0, 0.0) ;delta (area=1, position, fraction=1, DW=0.0)
;quasielastic terms
for i=1, nmax2 do begin
Y += (1-P3-P5-P8) * (2*i+1) * jl2[i] * strnormLorentz(X, P2, hwhm3[i])
endfor
;scaling & DW
Y *= P1
if P4 gt 0. then Y *= EXP((-P4*Qcurrent1D^2/3.)>(-70.))
return, Y
end
;
;;End QENS functions
;******************************************************************************
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment