diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp
index a9c1c05c..d1be50fd 100644
--- a/src/classes/PTheory.cpp
+++ b/src/classes/PTheory.cpp
@@ -471,16 +471,19 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
           return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
         case THEORY_STATIC_ZF_NK:
-          return StaticNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+          return StaticNKZF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);          
         case THEORY_STATIC_TF_NK:
-          return StaticNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+          return StaticNKTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
         case THEORY_DYNAMIC_ZF_NK:
-          return DynamicNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+          return DynamicNKZF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
         case THEORY_DYNAMIC_TF_NK:
-          return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+          return DynamicNKTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+                 fAdd->Func(t, paramValues, funcValues);
+        case THEORY_F_MU_F:
+          return FmuF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
         case THEORY_POLYNOM:
           return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
@@ -560,6 +563,8 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
           return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
         case THEORY_MU_MINUS_EXP:
           return MuMinusExpTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+        case THEORY_F_MU_F:
+          return FmuF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
         case THEORY_POLYNOM:
           return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
         case THEORY_USER_FCN:
@@ -635,6 +640,8 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
           return DynamicNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
         case THEORY_MU_MINUS_EXP:
           return MuMinusExpTF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+        case THEORY_F_MU_F:
+            return FmuF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
         case THEORY_POLYNOM:
           return Polynom(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
         case THEORY_USER_FCN:
@@ -708,6 +715,8 @@ Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoub
           return DynamicNKTF(t, paramValues, funcValues);
         case THEORY_MU_MINUS_EXP:
           return MuMinusExpTF(t, paramValues, funcValues);
+        case THEORY_F_MU_F:
+          return FmuF(t, paramValues, funcValues);
         case THEORY_POLYNOM:
           return Polynom(t, paramValues, funcValues);
         case THEORY_USER_FCN:
@@ -2590,6 +2599,49 @@ Double_t PTheory::DynamicNKTF(Double_t t, const PDoubleVector& paramValues, cons
   return result;
 }
 
+//--------------------------------------------------------------------------
+/**
+ * <p>F-\f$\mu\f-F polaritation function.
+ * For details see e.g. "Muon Spectroscopy - An Introduction", S. Blundell, etal.
+ *
+ * @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::FmuF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+{
+    // expected paramters: w_d [0], [tshift [1]]
+
+    Double_t val[2];
+
+    assert(fParamNo.size() <= 2);
+
+    if (t < 0.0)
+        return 1.0;
+
+    // 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() == 1) // no tshift
+        tt = t;
+    else // tshift present
+        tt = t-val[1];
+
+    const Double_t sqrt3 = sqrt(3.0);
+    const Double_t wd_t = val[0]*tt;
+
+    return (3.0+cos(sqrt3*wd_t)+(1.0-1.0/sqrt3)*cos(((3.0-sqrt3)/2.0)*wd_t)+(1.0+1.0/sqrt3)*cos(((3.0 + sqrt3)/2.0)*wd_t))/6.0;
+}
+
 //--------------------------------------------------------------------------
 /**
  * <p> theory function: polynom
diff --git a/src/include/PTheory.h b/src/include/PTheory.h
index cb01317e..ddbccf9c 100644
--- a/src/include/PTheory.h
+++ b/src/include/PTheory.h
@@ -73,9 +73,10 @@
 #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
+#define THEORY_F_MU_F                       30
+#define THEORY_MU_MINUS_EXP                 31
+#define THEORY_POLYNOM                      32
+#define THEORY_USER_FCN                     33
 
 // 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
@@ -109,10 +110,11 @@
 #define THEORY_PARAM_STATIC_TF_NK                 4 // phase, frequency, damping D0, R_b=DGbG/D0 (tshift)
 #define THEORY_PARAM_DYNAMIC_ZF_NK                3 // damping D0, R_b=DGbG/D0, nu_c (tshift)
 #define THEORY_PARAM_DYNAMIC_TF_NK                5 // phase, frequency, damping D0, R_b=DGbG/D0, nu_c (tshift)
+#define THEORY_PARAM_F_MU_F                       1 // frequency (tshift)
 #define THEORY_PARAM_MU_MINUS_EXP                 6 // N0, tau, A, damping, phase, frequency (tshift)
 
 // number of available user functions
-#define THEORY_MAX 33
+#define THEORY_MAX 34
 
 // maximal number of parameters. Needed in the contents of LF
 #define THEORY_MAX_PARAM 10
@@ -234,6 +236,9 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = {
   {THEORY_DYNAMIC_TF_NK, THEORY_PARAM_DYNAMIC_TF_NK, false,
    "dynamicNKTF", "dnktf", "(phase frequency damping_D0 R_b nu_c)", "(phase frequency damping_D0 R_b nu_c tshift)"},
 
+  {THEORY_F_MU_F, THEORY_PARAM_F_MU_F, false,
+   "F_mu_F", "fmuf", "(frequency)", "(frequency tshift)"},
+
   {THEORY_MU_MINUS_EXP, THEORY_PARAM_MU_MINUS_EXP, false,
    "muMinusExpTF", "mmsetf", "(N0 tau A lambda phase nu)", "(N0 tau A lambda phase nu tshift)"},
 
@@ -295,6 +300,7 @@ class PTheory
     virtual Double_t StaticNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t DynamicNKZF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t DynamicNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t FmuF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t MuMinusExpTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t Polynom(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t UserFcn(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;