Add DynamicKuboToyabe fitting function

Code/Mantid/Framework/CurveFitting/CMakeLists.txt

@@ -29,6 +29,7 @@ set ( SRC_FILES
 	src/DerivMinimizer.cpp
 	src/DiffRotDiscreteCircle.cpp
 	src/DiffSphere.cpp
+	src/DynamicKuboToyabe.cpp
 	src/EndErfc.cpp
 	src/ExpDecay.cpp
 	src/ExpDecayMuon.cpp
@@ -135,6 +136,7 @@ set ( INC_FILES
 	inc/MantidCurveFitting/DiffRotDiscreteCircle.h
 	inc/MantidCurveFitting/DiffSphere.h
 	inc/MantidCurveFitting/DllConfig.h
+	inc/MantidCurveFitting/DynamicKuboToyabe.h
 	inc/MantidCurveFitting/EndErfc.h
 	inc/MantidCurveFitting/ExpDecay.h
 	inc/MantidCurveFitting/ExpDecayMuon.h
@@ -237,6 +239,7 @@ set ( TEST_FILES
 	DeltaFunctionTest.h
 	DiffRotDiscreteCircleTest.h
 	DiffSphereTest.h
+	DynamicKuboToyabeTest.h
 	EndErfcTest.h
 	ExpDecayMuonTest.h
 	ExpDecayOscTest.h

Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/DynamicKuboToyabe.h

@@ -0,0 +1,72 @@
+#ifndef MANTID_CURVEFITTING_DYNAMICKUBOTOYABE_H_
+#define MANTID_CURVEFITTING_DYNAMICKUBOTOYABE_H_
+
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidAPI/IPeakFunction.h"
+#include "MantidAPI/IFunctionMW.h"
+#include "MantidAPI/IFunctionWithLocation.h"
+#include <cmath>
+
+//#include "MantidAPI/ParamFunction.h"
+//#include "MantidAPI/IPeakFunction.h"
+
+//#include "MantidAPI/IFunctionMW.h"
+//#include "MantidAPI/IFunction1D.h"
+
+namespace Mantid
+{
+  namespace CurveFitting
+  {
+    /** 
+     Provide Dynamic Kubo Toyabe function interface to IPeakFunction for muon scientists.
+
+     @author Karl Palmen, ISIS, RAL 
+     @date 21/03/2012 
+
+     Copyright &copy; 2007-2012 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory 
+
+     This file is part of Mantid. 
+
+     Mantid is free software; you can redistribute it and/or modify 
+     it under the terms of the GNU General Public License as published by 
+     the Free Software Foundation; either version 3 of the License, or 
+     (at your option) any later version. 
+
+     Mantid is distributed in the hope that it will be useful, 
+     but WITHOUT ANY WARRANTY; without even the implied warranty of 
+     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
+     GNU General Public License for more details. 
+
+     You should have received a copy of the GNU General Public License 
+     along with this program.  If not, see <http://www.gnu.org/licenses/>. 
+
+     File change history is stored at: <https://github.com/mantidproject/mantid> 
+     Code Documentation is available at: <http://doxygen.mantidproject.org> 
+     */ 
+
+    class DLLExport DynamicKuboToyabe :  public API::ParamFunction, public API::IFunction1D
+    {
+    public:
+
+      /// Destructor
+      virtual ~DynamicKuboToyabe() {}
+
+      /// overwrite base class methods
+      std::string name()const{return "DynamicKuboToyabe";}
+      virtual const std::string category() const { return "Muon";}
+
+    protected:
+      virtual void function1D(double* out, const double* xValues, const size_t nData)const;
+      virtual void functionDeriv1D(API::Jacobian* out, const double* xValues, const size_t nData);
+      virtual void functionDeriv(const API::FunctionDomain& domain, API::Jacobian& jacobian);
+      virtual void init();
+      virtual void setActiveParameter(size_t i, double value);
+
+    };
+
+  } // namespace CurveFitting
+} // namespace Mantid
+
+#endif /*MANTID_CURVEFITTING_DYNAMICKUBOTOYABE_H_*/

Code/Mantid/Framework/CurveFitting/src/DynamicKuboToyabe.cpp

@@ -0,0 +1,170 @@
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidCurveFitting/DynamicKuboToyabe.h"
+#include "MantidAPI/Jacobian.h"
+#include "MantidAPI/FunctionFactory.h"
+#include <vector>
+
+namespace Mantid
+{
+namespace CurveFitting
+{
+
+using namespace Kernel;
+using namespace API;
+
+DECLARE_FUNCTION(DynamicKuboToyabe)
+
+// ** MODIFY THIS **
+// Here specify/declare the parameters of your Fit Function
+// 
+// declareParameter takes three arguments:
+//
+//   1st: The name of the parameter
+//   2nd: The default (initial) value of the parameter
+//   3rd: A description of the parameter (optional)
+//
+void DynamicKuboToyabe::init()
+{
+  declareParameter("Asym",  0.2, "Amplitude at time 0");
+  declareParameter("Delta", 0.2, "Local field");
+  declareParameter("Field", 0.0, "External field");
+  declareParameter("Nu",    0.0, "Hopping rate");
+}
+
+// Static Zero Field Kubo Toyabe relaxation function
+double ZFKT (const double x, const double G){
+
+  const double q = G*G*x*x;
+  return (0.3333333333 + 0.6666666667*exp(-0.5*q)*(1-q));
+}
+
+// Static Non-Zero field Kubo Toyabe relaxation function
+//double HKT (const double x, const double G, const double F) 
+//{
+//  throw std::runtime_error("HKT not implemented yet");
+//}
+
+// Dynamic Kubo-Toyabe
+double getDKT (double t, double G, double v){
+
+  const int tsmax = 656; // Length of the time axis, 32 us of valid data
+  const double eps = 0.05; // Bin width for calculations
+  //  const int stk = 25; // Not used for the moment
+
+  static double oldG=-1., oldV=-1.;
+  static std::vector<double> gStat(tsmax), gDyn(tsmax);
+
+
+  if ( (G != oldG) || (v != oldV) ){
+
+    // If G or v have changed with respect to the 
+    // previous call, we need to re-do the computations
+
+
+    if ( G != oldG ){
+
+      // But we only need to
+      // re-compute gStat if G has changed
+
+      // Generate static Kubo-Toyabe
+      for (int k=0; k<tsmax; k++){
+        gStat[k]= ZFKT(k*eps,G);
+      }
+      // Store new G value
+      oldG =G;
+    }
+
+    // Store new v value
+    oldV =v;
+
+    double hop = v*eps;
+
+    // Generate dynamic Kubo Toyabe
+    for (int k=0; k<tsmax; k++){
+      double y=gStat[k];
+      for (int j=k-1; j>0; j--){
+        y=y*(1-hop)+hop*gDyn[k-j]*gStat[j];
+      }
+      gDyn[k]=y;
+    }
+  }
+
+  // Interpolate table 
+  // If beyond end, extrapolate
+  int x=int(fabs(t)/eps);
+  if (x>tsmax-2)
+    x = tsmax-2;
+  double xe=(fabs(t)/eps)-x;
+  return gDyn[x]*(1-xe)+xe*gDyn[x+1];
+
+}
+
+// Dynamic Kubo Toyabe function
+void DynamicKuboToyabe::function1D(double* out, const double* xValues, const size_t nData)const
+{
+  const double& A = getParameter("Asym");
+  const double& G = fabs(getParameter("Delta"));
+  const double& F = fabs(getParameter("Field"));
+  const double& v = fabs(getParameter("Nu"));
+
+
+    // Zero hopping rate
+	if (v == 0.0) {
+
+    // Zero external field
+    if ( F == 0.0 ){
+      for (size_t i = 0; i < nData; i++) {
+        out[i] = A*ZFKT(xValues[i],G);
+      }
+    }
+    // Non-zero external field
+    else{
+      //for (size_t i = 0; i < nData; i++) {
+      //  out[i] = A*HKT(xValues[i],G,F);
+      //}
+      throw std::runtime_error("HKT() not implemented yet");
+    }
+	} 
+
+  // Non-zero hopping rate
+	else {
+
+    if ( F==0.0 ) {
+
+      for (size_t i = 0; i<nData; i++){
+        out[i] = A*getDKT(xValues[i],G,v);
+      }
+
+	  } else {
+
+	    // Non-zero field
+      throw std::runtime_error("HKT() not implemented yet");
+    }
+
+	} // else hopping rate != 0
+
+
+}
+
+
+
+void DynamicKuboToyabe::functionDeriv(const API::FunctionDomain& domain, API::Jacobian& jacobian)
+{
+  calNumericalDeriv(domain, jacobian);
+}
+
+void DynamicKuboToyabe::functionDeriv1D(API::Jacobian* , const double* , const size_t )
+{
+  throw Mantid::Kernel::Exception::NotImplementedError("functionDerivLocal is not implemented for DynamicKuboToyabe.");
+}
+
+void DynamicKuboToyabe::setActiveParameter(size_t i, double value) {
+
+  setParameter( i, fabs(value), false);
+
+}
+
+} // namespace CurveFitting
+} // namespace Mantid