added new internal field functions.

ChangeLog

@@ -1,11 +1,11 @@
 #---------------------------------------------------------------------
 # ChangeLog
-#
-# $Id$
 #---------------------------------------------------------------------
 
 changes since 0.11.0
 ===================================
+NEW 2014-02-07 added new internal field functions internFldGK and internFldLL to
+               PTheory. Online documentation still missing.
 NEW 2013-09-10 add a stretched Kubo-Toyabe function to the PTheory class, including 
                an update of the docu.
 NEW 2013-02-14 (i) adding t0 estimate for single histogram fits. The estimate

configure.ac

@@ -1,7 +1,7 @@
 AC_REVISION([$Id$])
 
-AC_PREREQ(2.59)
+AC_PREREQ(2.63)
-AC_INIT(musrfit, 0.11.0, andreas.suter@psi.ch)
+AC_INIT([musrfit],[0.11.0],[andreas.suter@psi.ch])
 AC_CONFIG_AUX_DIR(admin)
 AC_CANONICAL_HOST
 #AC_MSG_RESULT([${host} ${host_cpu} ${host_vendor} ${host_os}])
@@ -218,8 +218,8 @@ m4_ifdef([LT_INIT],
           LT_LANG([C++])
           LT_LANG([C])],
          [m4_defun([_LT_AC_LANG_F77_CONFIG], [:])
-          AM_DISABLE_STATIC
+          AC_DISABLE_STATIC([])
-          AM_PROG_LIBTOOL
+          LT_INIT
           AC_PROG_CXX
           AC_PROG_CC])
 
@@ -229,7 +229,7 @@ dnl -----------------------------------------------
 
 FFTW3_FOUND=0
 AC_ARG_WITH([fftw3],
-    [AC_HELP_STRING([--with-fftw3],[prefix of FFTW3 installation, e.g. /usr/local or /usr])],
+    [AS_HELP_STRING([--with-fftw3],[prefix of FFTW3 installation, e.g. /usr/local or /usr])],
     [FFTW3_PREFIX=$with_fftw3
      AC_MSG_CHECKING([whether FFTW3 can be found at the specified location])
      if !(test -r ${FFTW3_PREFIX}/include/fftw3.h); then
@@ -275,7 +275,7 @@ dnl -----------------------------------------------
 
 GSL_FOUND=0
 AC_ARG_WITH([gsl],
-    [AC_HELP_STRING([--with-gsl],[prefix of GSL installation. e.g. /usr/local or /usr])],
+    [AS_HELP_STRING([--with-gsl],[prefix of GSL installation. e.g. /usr/local or /usr])],
     [GSL_PREFIX=$with_gsl
      AC_MSG_CHECKING([whether GSL can be found at the specified location])
      if !(test -r ${GSL_PREFIX}/include/gsl/gsl_sf_hyperg.h); then
@@ -323,7 +323,7 @@ dnl -----------------------------------------------
 
 BOOST_FOUND=0
 AC_ARG_WITH([boostinc],
-    [AC_HELP_STRING([--with-boostinc],[path to the header files of the BOOST installation, e.g. /usr/local/include or /usr/include])],
+    [AS_HELP_STRING([--with-boostinc],[path to the header files of the BOOST installation, e.g. /usr/local/include or /usr/include])],
     [BOOST_INCLUDE=$with_boostinc
      AC_MSG_CHECKING([whether BOOST can be found at the specified location])
      if !(test -r ${BOOST_INCLUDE}/boost/spirit.hpp) && !(test -r ${BOOST_INCLUDE}/boost/spirit/include/classic_spirit.hpp); then
@@ -371,7 +371,7 @@ dnl Ask user for path to libxml2
 dnl -----------------------------------------------
 LIBXML2_FOUND=0
 AC_ARG_WITH([libxml2],
-    [AC_HELP_STRING([--with-libxml2],[prefix of the libxml2 installation, e.g. /usr/local/include])],
+    [AS_HELP_STRING([--with-libxml2],[prefix of the libxml2 installation, e.g. /usr/local/include])],
     [LIBXML2_PREFIX=$with_libxml2
     AC_MSG_CHECKING([whether libxml2 can be found at the specified location])
      if !(test -r ${LIBXML2_PREFIX}/libxml2/libxml/xmlreader.h) && !(test -r ${LIBXML2_PREFIX}/libxml2/libxml/parser.h) && !(test -r ${LIBXML2_PREFIX}/libxml2/libxml/xmlschemas.h); then
@@ -466,11 +466,11 @@ AC_SUBST(ROOT_CFLAGS)
 dnl -----------------------------------------------
 dnl Check for NeXus, HDF5, etc.
 dnl -----------------------------------------------
-AC_ARG_ENABLE([NeXus], [AC_HELP_STRING([--enable-NeXus], [build optional NeXus support [default=no]])],
+AC_ARG_ENABLE([NeXus], [AS_HELP_STRING([--enable-NeXus],[build optional NeXus support [default=no]])],
     [
     # checking for hdf4
     AC_ARG_WITH([hdf4],
-        [AC_HELP_STRING([--with-hdf4], [prefix of the HDF4 installation (needed for NeXus), e.g. /usr/local/hdf4 or /opt/hdf4])],
+        [AS_HELP_STRING([--with-hdf4],[prefix of the HDF4 installation (needed for NeXus), e.g. /usr/local/hdf4 or /opt/hdf4])],
         [HDF4_PREFIX=$with_hdf4
          AC_MSG_CHECKING([whether HDF4 can be found at the specified location])
          if !(test -r ${HDF4_PREFIX}/include/hdf.h); then
@@ -510,7 +510,7 @@ AC_ARG_ENABLE([NeXus], [AC_HELP_STRING([--enable-NeXus], [build optional NeXus s
 
     # checking for hdf5
     AC_ARG_WITH([hdf5],
-        [AC_HELP_STRING([--with-hdf5], [prefix of the HDF5 installation (needed for NeXus), e.g. /usr/local/hdf5 or /opt/hdf5])],
+        [AS_HELP_STRING([--with-hdf5],[prefix of the HDF5 installation (needed for NeXus), e.g. /usr/local/hdf5 or /opt/hdf5])],
         [HDF5_PREFIX=$with_hdf5
          AC_MSG_CHECKING([whether HDF5 can be found at the specified location])
          if !(test -r ${HDF5_PREFIX}/include/hdf5.h); then
@@ -548,7 +548,7 @@ AC_ARG_ENABLE([NeXus], [AC_HELP_STRING([--enable-NeXus], [build optional NeXus s
     # checking for NeXus
     NEXUS_FOUND=0
     AC_ARG_WITH([nexus],
-        [AC_HELP_STRING([--with-nexus], [prefix of the NeXus installation, e.g. /usr/local])],
+        [AS_HELP_STRING([--with-nexus],[prefix of the NeXus installation, e.g. /usr/local])],
         [NEXUS_PREFIX=$with_nexus
          AC_MSG_CHECKING([whether NeXus can be found at the specified location])
          if !(test -r ${NEXUS_PREFIX}/include/napi.h); then
@@ -601,11 +601,11 @@ dnl -----------------------------------------------
 dnl Ask user if BMWlibs should be built and where to find the CUBA installation
 dnl -----------------------------------------------
 
-AC_ARG_ENABLE([BMWlibs], [AC_HELP_STRING([--enable-BMWlibs],[build optional BMW plug-ins [default=no]])],
+AC_ARG_ENABLE([BMWlibs], [AS_HELP_STRING([--enable-BMWlibs],[build optional BMW plug-ins [default=no]])],
     [CUBA_FOUND=0
     BUILD_CUBA=0
     AC_ARG_WITH([cuba],
-        [AC_HELP_STRING([--with-cuba],[prefix of CUBA installation, e.g. /usr/local or /usr])],
+        [AS_HELP_STRING([--with-cuba],[prefix of CUBA installation, e.g. /usr/local or /usr])],
         [CUBA_PREFIX=$with_cuba
         AC_MSG_CHECKING([whether CUBA can be found at the specified location])
         if !(test -r ${CUBA_PREFIX}/include/cuba.h); then
@@ -653,7 +653,17 @@ AC_ARG_ENABLE([BMWlibs], [AC_HELP_STRING([--enable-BMWlibs],[build optional BMW
 
       AC_C_CONST
       AC_C_INLINE
-      AC_C_LONG_DOUBLE
+      AC_DIAGNOSE([obsolete],[The macro `AC_C_LONG_DOUBLE' is obsolete.
+You should use `AC_TYPE_LONG_DOUBLE' or `AC_TYPE_LONG_DOUBLE_WIDER' instead.
+])dnl
+
+    AC_TYPE_LONG_DOUBLE_WIDER
+    ac_cv_c_long_double=$ac_cv_type_long_double_wider
+    if test $ac_cv_c_long_double = yes; then
+      AC_DEFINE([HAVE_LONG_DOUBLE],[1],[Define to 1 if the type `long double' works and has more range or
+	 precision than `double'.])
+    fi
+  
 
       AC_CHECK_FUNCS([powl])
       AC_CHECK_FUNCS([erf])
@@ -721,7 +731,7 @@ dnl -----------------------------------------------
 dnl Ask user if ASlibs should be built
 dnl -----------------------------------------------
 
-AC_ARG_ENABLE([ASlibs], [AC_HELP_STRING([--enable-ASlibs],[build optional AS plug-ins [default=no]])],[BUILD_AS_LIBS=1], [BUILD_AS_LIBS=0])
+AC_ARG_ENABLE([ASlibs], [AS_HELP_STRING([--enable-ASlibs],[build optional AS plug-ins [default=no]])],[BUILD_AS_LIBS=1], [BUILD_AS_LIBS=0])
 
 dnl -----------------------------------------------
 dnl Set some paths and flags for PMusr, MusrRoot, TLemRunHeader, Class_MuSR_PSI, mud, NeXus (if enabled), etc.
@@ -810,7 +820,7 @@ dnl -----------------------------------------------
 dnl Ask user if OpenMP support should be disabled (used for parallel chisq calculation and in libFitPofB)
 dnl -----------------------------------------------
 
-AC_ARG_ENABLE([omp], [AC_HELP_STRING([--enable-omp],[build musrfit with OpenMP support [default=yes]])])
+AC_ARG_ENABLE([omp], [AS_HELP_STRING([--enable-omp],[build musrfit with OpenMP support [default=yes]])])
 
 if test "x$enable_omp" != "xno"; then
   SAVED_CXXFLAGS="$CXXFLAGS"
@@ -825,7 +835,7 @@ dnl -----------------------------------------------
 dnl Ask user if the building of musredit/musrgui should be disabled
 dnl -----------------------------------------------
 
-AC_ARG_ENABLE([editor], [AC_HELP_STRING([--enable-editor],[check for Qt and build musredit or musrgui [default=yes]])])
+AC_ARG_ENABLE([editor], [AS_HELP_STRING([--enable-editor],[check for Qt and build musredit or musrgui [default=yes]])])
 
 if test "x$enable_editor" != "xno"; then
 
@@ -847,7 +857,7 @@ if test "x$enable_editor" != "xno"; then
   QTEDITOR=none
 
   AC_ARG_WITH([qt3],
-        [AC_HELP_STRING([--with-qt3],[prefix of the Qt3 installation, e.g. /usr/lib/qt-3.3])],
+        [AS_HELP_STRING([--with-qt3],[prefix of the Qt3 installation, e.g. /usr/lib/qt-3.3])],
         [QT3_PREFIX=$with_qt3
         AC_MSG_CHECKING([whether qmake for Qt3 can be found at the specified location])
         for QMAKE in qmake qmake-qt3
@@ -876,7 +886,7 @@ if test "x$enable_editor" != "xno"; then
   )
 
   AC_ARG_WITH([qt4],
-        [AC_HELP_STRING([--with-qt4],[prefix of the Qt4 installation, e.g. /opt/qt-4.7])],
+        [AS_HELP_STRING([--with-qt4],[prefix of the Qt4 installation, e.g. /opt/qt-4.7])],
         [QT4_PREFIX=$with_qt4
         AC_MSG_CHECKING([whether qmake for Qt4 can be found at the specified location])
         for QMAKE in qmake qmake-qt4

src/classes/PTheory.cpp

@@ -5,12 +5,10 @@
   Author: Andreas Suter
   e-mail: andreas.suter@psi.ch
 
-  $Id$
-
 ***************************************************************************/
 
 /***************************************************************************
- *   Copyright (C) 2007-2013 by Andreas Suter                              *
+ *   Copyright (C) 2007-2014 by Andreas Suter                              *
  *   andreas.suter@psi.ch                                                  *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
@@ -460,13 +458,21 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
           return Abragam(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + 
                  fAdd->Func(t, paramValues, funcValues);
           break;
+        case THEORY_TF_COS:
+          return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+                 fAdd->Func(t, paramValues, funcValues);
+          break;
         case THEORY_INTERNAL_FIELD:
           return InternalField(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
           break;
-        case THEORY_TF_COS:
+        case THEORY_INTERNAL_FIELD_KORNILOV:
-          return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + 
+          return InternalFieldGK(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
-                 fAdd->Func(t, paramValues, funcValues);
+              fAdd->Func(t, paramValues, funcValues);
+          break;
+        case THEORY_INTERNAL_FIELD_LARKIN:
+          return InternalFieldLL(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+              fAdd->Func(t, paramValues, funcValues);
           break;
         case THEORY_BESSEL:
           return Bessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + 
@@ -559,11 +565,17 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
         case THEORY_ABRAGAM:
           return Abragam(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
           break;
+        case THEORY_TF_COS:
+          return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+          break;
         case THEORY_INTERNAL_FIELD:
           return InternalField(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
           break;
-        case THEORY_TF_COS:
+        case THEORY_INTERNAL_FIELD_KORNILOV:
-          return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+          return InternalFieldGK(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+          break;
+        case THEORY_INTERNAL_FIELD_LARKIN:
+          return InternalFieldLL(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
           break;
         case THEORY_BESSEL:
           return Bessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
@@ -649,11 +661,17 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
         case THEORY_ABRAGAM:
           return Abragam(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
           break;
+        case THEORY_TF_COS:
+          return TFCos(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+          break;
         case THEORY_INTERNAL_FIELD:
           return InternalField(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
           break;
-        case THEORY_TF_COS:
+        case THEORY_INTERNAL_FIELD_KORNILOV:
-          return TFCos(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+          return InternalFieldGK(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+          break;
+        case THEORY_INTERNAL_FIELD_LARKIN:
+          return InternalFieldLL(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
           break;
         case THEORY_BESSEL:
           return Bessel(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
@@ -737,11 +755,17 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co
         case THEORY_ABRAGAM:
           return Abragam(t, paramValues, funcValues);
           break;
+        case THEORY_TF_COS:
+          return TFCos(t, paramValues, funcValues);
+          break;
         case THEORY_INTERNAL_FIELD:
           return InternalField(t, paramValues, funcValues);
           break;
-        case THEORY_TF_COS:
+        case THEORY_INTERNAL_FIELD_KORNILOV:
-          return TFCos(t, paramValues, funcValues);
+          return InternalFieldGK(t, paramValues, funcValues);
+          break;
+        case THEORY_INTERNAL_FIELD_LARKIN:
+          return InternalFieldLL(t, paramValues, funcValues);
           break;
         case THEORY_BESSEL:
           return Bessel(t, paramValues, funcValues);
@@ -1955,6 +1979,47 @@ Double_t PTheory::Abragam(register Double_t t, const PDoubleVector& paramValues,
                     (TMath::Exp(-gamma_t)-1.0+gamma_t));
 }
 
+//--------------------------------------------------------------------------
+/**
+ * <p> theory function: cosine including phase
+ *
+ * \f[ = \cos(2\pi\nu t + \varphi) \f]
+ * or the time shifted version of it.
+ *
+ * <b>meaning of paramValues:</b> \f$\nu\f$, \f$\varphi\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::TFCos(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+{
+  // expected parameters: phase frequency [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];
+
+  return TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
+}
+
 //--------------------------------------------------------------------------
 /**
  * <p> theory function: internal field function
@@ -1999,12 +2064,12 @@ Double_t PTheory::InternalField(register Double_t t, const PDoubleVector& paramV
 
 //--------------------------------------------------------------------------
 /**
- * <p> theory function: cosine including phase
+ * <p> theory function: internal field function Gauss-Kornilov (see Physics Letters A <b>153</b>, 364 (1991)).
  *
- * \f[ = \cos(2\pi\nu t + \varphi) \f]
+ * \f[ = \alpha\,\left[\cos(2\pi\nu t)-\frac{\sigma^2 t}{2\pi\nu}\sin(2\pi\nu t)\right]\exp(-[\sigma t]^2/2)+(1-\alpha)\,\exp(-(\lambda t)^\beta)\f]
- * or the time shifted version of it.
+ * or the time shifted version of it. For the powder averaged case, \f$\alpha=2/3\f$.
  *
- * <b>meaning of paramValues:</b> \f$\nu\f$, \f$\varphi\f$ [, \f$t_{\rm shift}\f$]
+ * <b>meaning of paramValues:</b> \f$\alpha\f$, \f$\nu\f$, \f$\sigma\f$, \f$\lambda\f$, \f$\beta\f$ [, \f$t_{\rm shift}\f$]
  *
  * <b>return:</b> function value
  *
@@ -2012,13 +2077,13 @@ Double_t PTheory::InternalField(register Double_t t, const PDoubleVector& paramV
  * \param paramValues parameter values
  * \param funcValues vector with the functions (i.e. functions of the parameters)
  */
-Double_t PTheory::TFCos(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+Double_t PTheory::InternalFieldGK(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
 {
-  // expected parameters: phase frequency [tshift]
+  // expected parameters: [0]:fraction [1]:frequency [2]:sigma [3]:lambda [4]:beta [[5]:tshift]
 
-  Double_t val[3];
+  Double_t val[6];
 
-  assert(fParamNo.size() <= 3);
+  assert(fParamNo.size() <= 6);
 
   // check if FUNCTIONS are used
   for (UInt_t i=0; i<fParamNo.size(); i++) {
@@ -2030,12 +2095,75 @@ Double_t PTheory::TFCos(register Double_t t, const PDoubleVector& paramValues, c
   }
 
   Double_t tt;
-  if (fParamNo.size() == 2) // no tshift
+  if (fParamNo.size() == 5) // no tshift
     tt = t;
   else // tshift present
-    tt = t-val[2];
+    tt = t-val[5];
 
-  return TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
+  Double_t result = 0.0;
+  Double_t nu_t = val[1]*tt;
+  Double_t rateLF = TMath::Power(val[3]*tt, val[4]);
+  Double_t rate2 = val[2]*val[2]*tt*tt; // (sigma t)^2
+
+  if (val[1] < 0.01) { // internal field frequency is approaching zero
+    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(1.0-rate2)*TMath::Exp(-0.5*rate2);
+  } else {
+    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(TMath::Cos(nu_t)-val[2]*val[2]*tt/val[1]*TMath::Sin(nu_t))*TMath::Exp(-0.5*rate2);
+  }
+
+  return result;
+}
+
+//--------------------------------------------------------------------------
+/**
+ * <p> theory function: internal field function Lorentz-Larkin (see Physica B: Condensed Matter <b>289-290</b>, 153 (2000)).
+ *
+ * \f[ = \alpha\,\left[\cos(2\pi\nu t)-\frac{a}{2\pi\nu}\sin(2\pi\nu t)\right]\exp(-a t)+(1-\alpha)\,\exp(-(\lambda t)^\beta)\f]
+ * or the time shifted version of it. For the powder averaged case, \f$\alpha=2/3\f$.
+ *
+ * <b>meaning of paramValues:</b> \f$\alpha\f$, \f$\nu\f$, \f$a\f$, \f$\lambda\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::InternalFieldLL(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
+{
+  // expected parameters: [0]:fraction [1]:frequency [2]:a [3]:lambda [4]:beta [[5]:tshift]
+
+  Double_t val[6];
+
+  assert(fParamNo.size() <= 6);
+
+  // 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() == 5) // no tshift
+    tt = t;
+  else // tshift present
+    tt = t-val[5];
+
+  Double_t result = 0.0;
+  Double_t nu_t = val[1]*tt;
+  Double_t rateLF = TMath::Power(val[3]*tt, val[4]);
+  Double_t a_t = val[2]*tt; // a t
+
+  if (val[1] < 0.01) { // internal field frequency is approaching zero
+    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(1.0-a_t)*TMath::Exp(-a_t);
+  } else {
+    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(TMath::Cos(nu_t)-val[3]/val[1]*TMath::Sin(nu_t))*TMath::Exp(-a_t);
+  }
+
+  return result;
 }
 
 //--------------------------------------------------------------------------

src/include/PTheory.h

@@ -5,12 +5,10 @@
   Author: Andreas Suter
   e-mail: andreas.suter@psi.ch
 
-  $Id$
-
 ***************************************************************************/
 
 /***************************************************************************
- *   Copyright (C) 2007-2013 by Andreas Suter                              *
+ *   Copyright (C) 2007-2014 by Andreas Suter                              *
  *   andreas.suter@psi.ch                                                  *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
@@ -61,17 +59,19 @@
 #define THEORY_SPIN_GLASS                   13
 #define THEORY_RANDOM_ANISOTROPIC_HYPERFINE 14
 #define THEORY_ABRAGAM                      15
-#define THEORY_INTERNAL_FIELD               16
+#define THEORY_TF_COS                       16
-#define THEORY_TF_COS                       17
+#define THEORY_INTERNAL_FIELD               17
-#define THEORY_BESSEL                       18
+#define THEORY_INTERNAL_FIELD_KORNILOV      18
-#define THEORY_INTERNAL_BESSEL              19
+#define THEORY_INTERNAL_FIELD_LARKIN        19
-#define THEORY_SKEWED_GAUSS                 20
+#define THEORY_BESSEL                       20
-#define THEORY_STATIC_ZF_NK                 21
+#define THEORY_INTERNAL_BESSEL              21
-#define THEORY_STATIC_TF_NK                 22
+#define THEORY_SKEWED_GAUSS                 22
-#define THEORY_DYNAMIC_ZF_NK                23
+#define THEORY_STATIC_ZF_NK                 23
-#define THEORY_DYNAMIC_TF_NK                24
+#define THEORY_STATIC_TF_NK                 24
-#define THEORY_POLYNOM                      25
+#define THEORY_DYNAMIC_ZF_NK                25
-#define THEORY_USER_FCN                     26
+#define THEORY_DYNAMIC_TF_NK                26
+#define THEORY_POLYNOM                      27
+#define THEORY_USER_FCN                     28
 
 // 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
@@ -91,8 +91,10 @@
 #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)
-#define THEORY_PARAM_INTERNAL_FIELD               5 // fraction, phase, frequency, TF damping, damping (tshift)
 #define THEORY_PARAM_TF_COS                       2 // phase, frequency (tshift)
+#define THEORY_PARAM_INTERNAL_FIELD               5 // fraction, phase, frequency, TF damping, damping (tshift)
+#define THEORY_PARAM_INTERNAL_FIELD_KORNILOV      5 // fraction, frequency, TF damping, damping, beta (tshift)
+#define THEORY_PARAM_INTERNAL_FIELD_LARKIN        4 // fraction, frequency, TF damping, damping (tshift)
 #define THEORY_PARAM_BESSEL                       2 // phase, frequency (tshift)
 #define THEORY_PARAM_INTERNAL_BESSEL              5 // fraction, phase, frequency, TF damping, LF damping (tshift)
 #define THEORY_PARAM_SKEWED_GAUSS                 4 // phase, frequency, rate minus, rate plus (tshift)
@@ -102,7 +104,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 27
+#define THEORY_MAX 29
 
 // maximal number of parameters. Needed in the contents of LF
 #define THEORY_MAX_PARAM 10
@@ -182,11 +184,17 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = {
   {THEORY_ABRAGAM, THEORY_PARAM_ABRAGAM, false,
    "abragam", "ab", "(rate hopprate)", "(rate hopprate tshift)"},
 
+  {THEORY_TF_COS, THEORY_PARAM_TF_COS, false,
+   "TFieldCos", "tf", "(phase frequency)", "(phase frequency tshift)"},
+
   {THEORY_INTERNAL_FIELD, THEORY_PARAM_INTERNAL_FIELD, false,
    "internFld", "if", "(fraction phase frequency Trate Lrate)", "(fraction phase frequency Trate Lrate tshift)"},
 
-  {THEORY_TF_COS, THEORY_PARAM_TF_COS, false,
+  {THEORY_INTERNAL_FIELD_KORNILOV, THEORY_PARAM_INTERNAL_FIELD_KORNILOV, false,
-   "TFieldCos", "tf", "(phase frequency)", "(phase frequency tshift)"},
+   "internFldGK", "ifgk", "(fraction phase frequency Trate Lrate)", "(fraction phase frequency Trate Lrate tshift)"},
+
+  {THEORY_INTERNAL_FIELD_LARKIN, THEORY_PARAM_INTERNAL_FIELD_LARKIN, false,
+   "internFldLL", "ifll", "(fraction phase frequency Trate Lrate)", "(fraction phase frequency Trate Lrate tshift)"},
 
   {THEORY_BESSEL, THEORY_PARAM_BESSEL, false,
    "bessel", "b", "(phase frequency)", "(phase frequency tshift)"},
@@ -253,8 +261,10 @@ class PTheory
     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;
-    virtual Double_t InternalField(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t TFCos(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t InternalField(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t InternalFieldGK(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
+    virtual Double_t InternalFieldLL(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t Bessel(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t InternalBessel(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual Double_t SkewedGauss(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;