add a stretched Kubo-Toyabe function to the PTheory class, including an update of the docu. Motivated by Z. Salman, which also provided a patch, which needed some improvements.

2013-09-10 09:03:21 +00:00
parent e1d99a07cb
commit 3443820e83
18 changed files with 311 additions and 182 deletions
--- a/src/classes/PTheory.cpp
+++ b/src/classes/PTheory.cpp
@ -444,6 +444,10 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
          return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
          break;
+        case THEORY_STR_KT:
+          return StrKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+                 fAdd->Func(t, paramValues, funcValues);
+          break;
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
@ -543,6 +547,9 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
        case THEORY_COMBI_LGKT:
          return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
          break;
+        case THEORY_STR_KT:
+          return StrKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+          break;
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
          break;
@ -630,6 +637,9 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
        case THEORY_COMBI_LGKT:
          return CombiLGKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
          break;
+        case THEORY_STR_KT:
+          return StrKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+          break;
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
          break;
@ -715,6 +725,9 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
        case THEORY_COMBI_LGKT:
          return CombiLGKT(t, paramValues, funcValues);
          break;
+        case THEORY_STR_KT:
+          return StrKT(t, paramValues, funcValues);
+          break;
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues);
          break;
@ -1752,6 +1765,56 @@ Double_t PTheory::CombiLGKT(register Double_t t, const PDoubleVector& paramValue
         (1.0 + 2.0*(1.0-lambdaL_t-lambdaG_t_2)*TMath::Exp(-(lambdaL_t+0.5*lambdaG_t_2)));
 }

+//--------------------------------------------------------------------------
+/**
+ * <p> theory function: stretched Kubo-Toyabe in zero field.
+ * Ref: Crook M. R. and Cywinski R., J. Phys. Condens. Matter, 9 (1997) 1149.
+ *
+ * \f[ =  1/3 + 2/3 \left(1-(\sigma t)^\beta \right) \exp\left(-(\sigma t)^\beta / \beta\right)\f]
+ * or the time shifted version of it.
+ *
+ * <b>meaning of paramValues:</b> \f$\sigma\f$, \f$\beta\f$ [, \f$t_{\rm shift}\f$]
+ *
+ * <b>return:</b> function value
+ *
+ * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit
+ * \param paramValues parameter values
+ * \param funcValues vector with the functions (i.e. functions of the parameters)
+ */
+Double_t PTheory::StrKT(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+{
+  // expected parameters: sigma beta [tshift]
+
+  Double_t val[3];
+
+  assert(fParamNo.size() <= 3);
+
+  // check if FUNCTIONS are used
+  for (UInt_t i=0; i<fParamNo.size(); i++) {
+    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
+      val[i] = paramValues[fParamNo[i]];
+    } else { // function
+      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
+    }
+  }
+
+  // check for beta too small (beta < 0.1) in which case numerical problems could arise and the function is anyhow
+  // almost identical to a constant of 1/3.
+  if (val[1] < 0.1)
+    return 0.333333333333333;
+
+  Double_t tt;
+  if (fParamNo.size() == 2) // no tshift
+    tt = t;
+  else // tshift present
+    tt = t-val[2];
+
+  Double_t sigma_t   = TMath::Power(tt*val[0],val[1]);
+
+  return 0.333333333333333 * 
+    (1.0 + 2.0*(1.0-sigma_t)*TMath::Exp(-sigma_t/val[1]));
+}
+
 //--------------------------------------------------------------------------
 /**
 * <p> theory function: spin glass function
--- a/src/include/PTheory.h
+++ b/src/include/PTheory.h
@ -57,20 +57,21 @@
 #define THEORY_STATIC_LORENTZ_KT_LF          9
 #define THEORY_DYNAMIC_LORENTZ_KT_LF        10
 #define THEORY_COMBI_LGKT                   11
-#define THEORY_SPIN_GLASS                   12
-#define THEORY_RANDOM_ANISOTROPIC_HYPERFINE 13
-#define THEORY_ABRAGAM                      14
-#define THEORY_INTERNAL_FIELD               15
-#define THEORY_TF_COS                       16
-#define THEORY_BESSEL                       17
-#define THEORY_INTERNAL_BESSEL              18
-#define THEORY_SKEWED_GAUSS                 19
-#define THEORY_STATIC_ZF_NK                 20
-#define THEORY_STATIC_TF_NK                 21
-#define THEORY_DYNAMIC_ZF_NK                22
-#define THEORY_DYNAMIC_TF_NK                23
-#define THEORY_POLYNOM                      24
-#define THEORY_USER_FCN                     25
+#define THEORY_STR_KT                       12
+#define THEORY_SPIN_GLASS                   13
+#define THEORY_RANDOM_ANISOTROPIC_HYPERFINE 14
+#define THEORY_ABRAGAM                      15
+#define THEORY_INTERNAL_FIELD               16
+#define THEORY_TF_COS                       17
+#define THEORY_BESSEL                       18
+#define THEORY_INTERNAL_BESSEL              19
+#define THEORY_SKEWED_GAUSS                 20
+#define THEORY_STATIC_ZF_NK                 21
+#define THEORY_STATIC_TF_NK                 22
+#define THEORY_DYNAMIC_ZF_NK                23
+#define THEORY_DYNAMIC_TF_NK                24
+#define THEORY_POLYNOM                      25
+#define THEORY_USER_FCN                     26

 // function parameter tags, i.e. how many parameters has a specific function
 // if there is a comment with a (tshift), the number of parameters is increased by one
@ -86,6 +87,7 @@
 #define THEORY_PARAM_STATIC_LORENTZ_KT_LF         2 // frequency, damping (tshift)
 #define THEORY_PARAM_DYNAMIC_LORENTZ_KT_LF        3 // frequency, damping, hop-rate (tshift)
 #define THEORY_PARAM_COMBI_LGKT                   2 // Lorentz rate, Gauss rate (tshift)
+#define THEORY_PARAM_STR_KT                       2 // rate, exponent (tshift)
 #define THEORY_PARAM_SPIN_GLASS                   3 // rate, hop-rate, order parameter (tshift)
 #define THEORY_PARAM_RANDOM_ANISOTROPIC_HYPERFINE 2 // frequency, rate (tshift)
 #define THEORY_PARAM_ABRAGAM                      2 // rate, hop-rate (tshift)
@ -100,7 +102,7 @@
 #define THEORY_PARAM_DYNAMIC_TF_NK                5 // phase, frequency, damping D0, R_b=DGbG/D0, nu_c (tshift)

 // number of available user functions
-#define THEORY_MAX 26
+#define THEORY_MAX 27

 // maximal number of parameters. Needed in the contents of LF
 #define THEORY_MAX_PARAM 10
@ -166,7 +168,10 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = {
   "dynExpKTLF", "dektlf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"},

  {THEORY_COMBI_LGKT, THEORY_PARAM_COMBI_LGKT, false,
-   "combiLGKT", "lgkt", "(LorentzRate GaussRate)", "(LorentzRate GaussRate tshift)"},
+   "combiLGKT", "lgkt", "(lorentzRate gaussRate)", "(lorentzRate gaussRate tshift)"},
+
+  {THEORY_STR_KT, THEORY_PARAM_STR_KT, false,
+   "strKT", "skt", "(rate beta)", "(rate beta tshift)"},

  {THEORY_SPIN_GLASS, THEORY_PARAM_SPIN_GLASS, false,
   "spinGlass",  "spg", "(rate hopprate order)", "(rate hopprate order tshift)"},
@ -244,6 +249,7 @@ class PTheory
    virtual Double_t StaticLorentzKTLF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
    virtual Double_t DynamicLorentzKTLF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
    virtual Double_t CombiLGKT(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t StrKT(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
    virtual Double_t SpinGlass(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
    virtual Double_t RandomAnisotropicHyperfine(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
    virtual Double_t Abragam(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
--- a/src/musredit/musredit_startup.xml.in
+++ b/src/musredit/musredit_startup.xml.in
@ -130,11 +130,18 @@
    </func>
    <func>
      <name>combiLGKT</name>
-      <comment>(LorentzRate GaussRate)</comment>
+      <comment>(lorentzRate gaussRate)</comment>
      <label>combined Lorentz-Gauss KT</label>
      <pixmap>combiLGKT.png</pixmap>
      <params>2</params>
    </func>
+    <func>
+      <name>strKT</name>
+      <comment>(rate beta)</comment>
+      <label>stretched Kubo-Toyabe</label>
+      <pixmap></pixmap>
+      <params>2</params>
+    </func>
    <func>
      <name>spinGlass</name>
      <comment>(rate hopping-rate order)</comment>
--- a/src/musrgui/musrgui_startup.xml.in
+++ b/src/musrgui/musrgui_startup.xml.in
@ -115,11 +115,18 @@
    </func>
    <func>
      <name>combiLGKT</name>
-      <comment>(LorentzRate GaussRate)</comment>
+      <comment>(lorentzRate gaussRate)</comment>
      <label>combined Lorentz-Gauss KT</label>
      <pixmap>combiLGKT.png</pixmap>
      <params>2</params>
    </func>
+    <func>
+      <name>strKT</name>
+      <comment>(rate beta)</comment>
+      <label>stretched Kubo-Toyabe</label>
+      <pixmap></pixmap>
+      <params>2</params>
+    </func>
    <func>
      <name>spinGlass</name>
      <comment>(rate hopping-rate order)</comment>