diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp
index 2de4dc7c..a9c1c05c 100644
--- a/src/classes/PTheory.cpp
+++ b/src/classes/PTheory.cpp
@@ -425,6 +425,15 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
         case THEORY_DYNAMIC_LORENTZ_KT_LF:
           return DynamicLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
+          return DynamicGauLorKTZFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+                 fAdd->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
+          return DynamicGauLorKTLFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+                 fAdd->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_KT_LF:
+          return DynamicGauLorKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+                 fAdd->Func(t, paramValues, funcValues);
         case THEORY_COMBI_LGKT:
           return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
@@ -511,6 +520,12 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
           return StaticLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
         case THEORY_DYNAMIC_LORENTZ_KT_LF:
           return DynamicLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
+          return DynamicGauLorKTZFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
+          return DynamicGauLorKTLFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_KT_LF:
+          return DynamicGauLorKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
         case THEORY_COMBI_LGKT:
           return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
         case THEORY_STR_KT:
@@ -580,6 +595,12 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
           return StaticLorentzKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
         case THEORY_DYNAMIC_LORENTZ_KT_LF:
           return DynamicLorentzKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
+          return DynamicGauLorKTZFFast(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
+          return DynamicGauLorKTLFFast(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_KT_LF:
+          return DynamicGauLorKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
         case THEORY_COMBI_LGKT:
           return CombiLGKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
         case THEORY_STR_KT:
@@ -647,6 +668,12 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
           return StaticLorentzKTLF(t, paramValues, funcValues);
         case THEORY_DYNAMIC_LORENTZ_KT_LF:
           return DynamicLorentzKTLF(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
+          return DynamicGauLorKTZFFast(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
+          return DynamicGauLorKTLFFast(t, paramValues, funcValues);
+        case THEORY_DYNAMIC_GAULOR_KT_LF:
+          return DynamicGauLorKTLF(t, paramValues, funcValues);
         case THEORY_COMBI_LGKT:
           return CombiLGKT(t, paramValues, funcValues);
         case THEORY_STR_KT:
@@ -1633,6 +1660,108 @@ Double_t PTheory::DynamicLorentzKTLF(Double_t t, const PDoubleVector& paramValue
 
 }
 
+//--------------------------------------------------------------------------
+/**
+ * <p>Local Gaussian, global Lorentzian approximation in the limit
+ * \f[ \nu_c \gg \gamma_\mu \Delta_{\rm L} \f] in ZF.
+ * For details see "Muon Spin Rotation, Relaxation, and Resonance",
+ * A. Yaouanc and P. Dalmas Sec. 6.4, Eq.(6.89).
+ *
+ * @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)
+ *
+ * @return Polarization value of this function
+ */
+Double_t PTheory::DynamicGauLorKTZFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+{
+  // expected parameters: damping hopping [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];
+    }
+  }
+
+  Double_t tt;
+  if (fParamNo.size() == 2) // no tshift
+    tt = t;
+  else // tshift present
+    tt = t-val[2];
+
+  Double_t nut = val[1]*tt;
+  return exp(-sqrt(4.0*pow(val[0]/val[1], 2.0)*(exp(-nut)-1.0+nut)));
+}
+
+//--------------------------------------------------------------------------
+/**
+ * <p>Local Gaussian, global Lorentzian approximation in the limit
+ * \f[ \nu_c \gg \gamma_\mu \Delta_{\rm L} \f] in LF.
+ * For details see "Muon Spin Rotation, Relaxation, and Resonance",
+ * A. Yaouanc and P. Dalmas Sec. 6.4, Eq.(6.93).
+ *
+ * @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)
+ *
+ * @return Polarization value of this function
+ */
+Double_t PTheory::DynamicGauLorKTLFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+{
+  // expected parameters: frequency damping hopping [tshift]
+
+  Double_t val[4];
+
+  assert(fParamNo.size() <= 4);
+
+  // 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];
+    }
+  }
+
+  Double_t tt;
+  if (fParamNo.size() == 3) // no tshift
+    tt = t;
+  else // tshift present
+    tt = t-val[3];
+
+  Double_t w0 = TMath::TwoPi()*val[0];
+  Double_t w0_2 = w0*w0;
+  Double_t nu_2 = val[2]*val[2];
+  Double_t nu_t = val[2]*tt;
+  Double_t w0_t = w0*tt;
+  Double_t Gamma_t = ((w0_2+nu_2)*nu_t+(w0_2-nu_2)*(1.0-exp(-nu_t)*cos(w0_t))-2.0*val[2]*w0*exp(-nu_t)*sin(w0_t))/pow(w0_2+nu_2,2.0);
+  if (Gamma_t < 0.0)
+    Gamma_t = 0.0;
+
+  return exp(-sqrt(4.0*val[1]*val[1]*Gamma_t));
+}
+
+//--------------------------------------------------------------------------
+/**
+ * @brief PTheory::DynamicGauLorKTLF
+ * @param t
+ * @param paramValues
+ * @param funcValues
+ * @return
+ */
+Double_t PTheory::DynamicGauLorKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+{
+    // NOT YET IMPLEMENTED. Will be Eq.(6.86)
+    return 0.0;
+}
+
 //--------------------------------------------------------------------------
 /**
  * <p> theory function: dynamic Lorentzain Kubo-Toyabe in longitudinal applied field
diff --git a/src/include/PTheory.h b/src/include/PTheory.h
index 06f94f29..cb01317e 100644
--- a/src/include/PTheory.h
+++ b/src/include/PTheory.h
@@ -54,25 +54,28 @@
 #define THEORY_STATIC_LORENTZ_KT             8
 #define THEORY_STATIC_LORENTZ_KT_LF          9
 #define THEORY_DYNAMIC_LORENTZ_KT_LF        10
-#define THEORY_COMBI_LGKT                   11
-#define THEORY_STR_KT                       12
-#define THEORY_SPIN_GLASS                   13
-#define THEORY_RANDOM_ANISOTROPIC_HYPERFINE 14
-#define THEORY_ABRAGAM                      15
-#define THEORY_TF_COS                       16
-#define THEORY_INTERNAL_FIELD               17
-#define THEORY_INTERNAL_FIELD_KORNILOV      18
-#define THEORY_INTERNAL_FIELD_LARKIN        19
-#define THEORY_BESSEL                       20
-#define THEORY_INTERNAL_BESSEL              21
-#define THEORY_SKEWED_GAUSS                 22
-#define THEORY_STATIC_ZF_NK                 23
-#define THEORY_STATIC_TF_NK                 24
-#define THEORY_DYNAMIC_ZF_NK                25
-#define THEORY_DYNAMIC_TF_NK                26
-#define THEORY_MU_MINUS_EXP                 27
-#define THEORY_POLYNOM                      28
-#define THEORY_USER_FCN                     29
+#define THEORY_DYNAMIC_GAULOR_FAST_KT_ZF    11
+#define THEORY_DYNAMIC_GAULOR_FAST_KT_LF    12
+#define THEORY_DYNAMIC_GAULOR_KT_LF         13
+#define THEORY_COMBI_LGKT                   14
+#define THEORY_STR_KT                       15
+#define THEORY_SPIN_GLASS                   16
+#define THEORY_RANDOM_ANISOTROPIC_HYPERFINE 17
+#define THEORY_ABRAGAM                      18
+#define THEORY_TF_COS                       19
+#define THEORY_INTERNAL_FIELD               20
+#define THEORY_INTERNAL_FIELD_KORNILOV      21
+#define THEORY_INTERNAL_FIELD_LARKIN        22
+#define THEORY_BESSEL                       23
+#define THEORY_INTERNAL_BESSEL              24
+#define THEORY_SKEWED_GAUSS                 25
+#define THEORY_STATIC_ZF_NK                 26
+#define THEORY_STATIC_TF_NK                 27
+#define THEORY_DYNAMIC_ZF_NK                28
+#define THEORY_DYNAMIC_TF_NK                29
+#define THEORY_MU_MINUS_EXP                 30
+#define THEORY_POLYNOM                      31
+#define THEORY_USER_FCN                     32
 
 // 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
@@ -87,6 +90,9 @@
 #define THEORY_PARAM_STATIC_LORENTZ_KT            1 // damping (tshift)
 #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_DYNAMIC_GAULOR_FAST_KT_ZF    2 // damping, hop-rate (tshift)
+#define THEORY_PARAM_DYNAMIC_GAULOR_FAST_KT_LF    3 // frequency, damping, hop-rate (tshift)
+#define THEORY_PARAM_DYNAMIC_GAULOR_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)
@@ -106,7 +112,7 @@
 #define THEORY_PARAM_MU_MINUS_EXP                 6 // N0, tau, A, damping, phase, frequency (tshift)
 
 // number of available user functions
-#define THEORY_MAX 30
+#define THEORY_MAX 33
 
 // maximal number of parameters. Needed in the contents of LF
 #define THEORY_MAX_PARAM 10
@@ -171,6 +177,15 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = {
   {THEORY_DYNAMIC_LORENTZ_KT_LF, THEORY_PARAM_DYNAMIC_LORENTZ_KT_LF, true,
    "dynExpKTLF", "dektlf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"},
 
+  {THEORY_DYNAMIC_GAULOR_FAST_KT_ZF, THEORY_PARAM_DYNAMIC_GAULOR_FAST_KT_ZF, true,
+   "dynGLKT_F_ZF", "dglktfzf", "(damping hopping-rate)", "(damping hopping-rate tshift)"},
+
+  {THEORY_DYNAMIC_GAULOR_FAST_KT_LF, THEORY_PARAM_DYNAMIC_GAULOR_FAST_KT_LF, true,
+   "dynGLKT_F_LF", "dglktflf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"},
+
+  {THEORY_DYNAMIC_GAULOR_KT_LF, THEORY_PARAM_DYNAMIC_GAULOR_KT_LF, true,
+   "dynGLKT_LF", "dglktlf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"},
+
   {THEORY_COMBI_LGKT, THEORY_PARAM_COMBI_LGKT, false,
    "combiLGKT", "lgkt", "(lorentzRate gaussRate)", "(lorentzRate gaussRate tshift)"},
 
@@ -261,6 +276,9 @@ class PTheory
     virtual Double_t StaticLorentzKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t StaticLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t DynamicLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t DynamicGauLorKTZFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t DynamicGauLorKTLFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t DynamicGauLorKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t CombiLGKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t StrKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t SpinGlass(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;