Changes in str_fit to address issues #28 and #29 related to QENS models

lamp_mac/str_fit.pro

@@ -5995,6 +5995,68 @@ function QENS_IL2, X, P1, P2, P3, P4, P5, P6, P7, P8
 end
 
 
+function QENS_IL3, X, P1, P2, P3, P4, P5, P6, P7, P8, P9
+  ;+
+  ;******************************************************************************
+  ;** Jump-diffusion translation x Confined translation (Volino) x Localized motion
+  ;**
+  ;** Parameters:
+  ;**   P1  = Scaling factor
+  ;**   P2  = Center
+  ;**   P3  = u^2 for DW term
+  ;**   P4  = D for long-range translation
+  ;**   P5  = Residence time for long-range translation
+  ;**   P6  = D for confined translation
+  ;**   P7  = u^2 for confined translation
+  ;**   P8  = EISF of localized lorentzian
+  ;**   P9  = HWHM of localized lorentzian
+  ;-
+
+  common c_str_fitqens, func_resol, Qcurrent1D, Qcurrent2D, mig2, mig3, mig4, mig5
+
+  ;width of long-range translational lorentzian
+  hwhm_t = P4 * Qcurrent1D^2 / (1.0 + P5 * P4 * Qcurrent1D^2)
+  
+  ;A_i and Gamma_i terms in Volino's model 
+  nmax = 100 ;number of lorentzians to use in sum (theoretically is infinite)
+  al = fltarr(nmax+1) & hwhm_v = fltarr(nmax+1)
+  arg = double(Qcurrent1D^2 * P7)
+  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, nmax do begin
+    if arg gt 0.0 then begin
+      al[i] = exparg*(arg^i)/factorial(i)
+      hwhm_v[i] = i*P6/P7
+    endif else al[i]=0.0
+  endfor
+
+  ;First term = A0_V * A0_L * Lorentz(Gamma_T)
+  Y = al[0] * P8 * strnormLorentz(X, P2, hwhm_t)
+  
+  ;Second term = A0_L * Sum(AV_i * Lorentz(Gamma_T+GammaV_i)
+  for i = 1, nmax do begin
+    Y += P8 * al[i] * strnormLorentz(X, P2, hwhm_t+hwhm_v[i])
+  endfor
+  
+  ;Third term = A0_V * (1-A0_L) * Lorentz(Gamma_T + Gamma_L)
+  Y += al[0] * (1-P8) * strnormLorentz(X, P2, hwhm_t+P9)
+  
+  ;Fourth term = (1-A0_L) * Sum(AV_i * * Lorentz(Gamma_T+GammaV_i+Gamma_L)
+  for i = 1, nmax do begin
+    Y += (1-P8) * al[i] * strnormLorentz(X, P2, hwhm_t+hwhm_v[i]+P9)
+  endfor
+
+  ;scaling & DW
+  Y *= P1
+  if P3 gt 0. then Y *= EXP((-P3*Qcurrent1D^2/3.)>(-70.))
+
+  return, Y
+
+end
+
+
 function QENS_Fixed_JumpsNsites, X, P1, P2, P3, P4, P5, P6, P7
   ;+
   ;********************************************************
@@ -6161,8 +6223,8 @@ function QENS_Fixed_JumpsNsitesLogNormDist, X, P1, P2, P3, P4, P5, P6, P7, P8
   endfor
 
   ;scaling x participation ratio & DW
-  Y *= P1 * P3
-  if P4 gt 0. then Y *= EXP((-P4*Qcurrent1D^2/3.)>(-70.))
+  Y *= P1
+  if P3 gt 0. then Y *= EXP((-P3*Qcurrent1D^2/3.)>(-70.))
 
   return, Y
 
@@ -6178,8 +6240,8 @@ function QENS_Mix1, X, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12
   ;**
   ;**   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.
+  ;**   that P4+P5+P8 < 1, so the fit can go to non-physical values
+  ;**   such that the fraction for the last component 1-P4-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.
@@ -6187,8 +6249,8 @@ function QENS_Mix1, X, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12
   ;** Parameters:
   ;**   P1  = Scaling factor
   ;**   P2  = Center
-  ;**   P3  = Fraction of atoms contributing to delta
-  ;**   P4  = u^2 for DW term
+  ;**   P3  = u^2 for DW term
+  ;**   P4  = Fraction of atoms contributing to delta
   ;**   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
@@ -6204,7 +6266,7 @@ function QENS_Mix1, X, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12
   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)
+  Y = QENS_delta(X, 1.0, P2, P4, 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)
@@ -6242,15 +6304,15 @@ function QENS_Mix1, X, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, 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)
+  Y += (1-P4-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])
+    Y += (1-P4-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.))
+  if P3 gt 0. then Y *= EXP((-P3*Qcurrent1D^2/3.)>(-70.))
 
   return, Y