Merge pull request #236 from mantidproject/feature/9556_MuonDKT_function

Add DynamicKuboToyabe fitting function
mantidproject · Feb 20, 2015 · a2f3876 · a2f3876
2 parents c416586 + 9be647a
commit a2f3876
Show file tree

Hide file tree

Showing 5 changed files with 394 additions and 0 deletions.
diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/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

diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/DynamicKuboToyabe.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/DynamicKuboToyabe.h
@@ -0,0 +1,66 @@
+#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>
+
+namespace Mantid
+{
+  namespace CurveFitting
+  {
+    /** 
+     Provide Dynamic Kubo Toyabe function interface to IFunction1D for muon scientists.
+   
+     @author Raquel Alvarez, ISIS, RAL 
+     @date 18/02/2015 
+  
+     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_*/
diff --git a/Code/Mantid/Framework/CurveFitting/src/DynamicKuboToyabe.cpp b/Code/Mantid/Framework/CurveFitting/src/DynamicKuboToyabe.cpp
@@ -0,0 +1,151 @@
+//----------------------------------------------------------------------
+// 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)
+
+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));
+}
+
+// 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
+
+  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{
+      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
diff --git a/Code/Mantid/Framework/CurveFitting/test/DynamicKuboToyabeTest.h b/Code/Mantid/Framework/CurveFitting/test/DynamicKuboToyabeTest.h
@@ -0,0 +1,141 @@
+#ifndef DYNAMICKUBOTOYABETEST_H_
+#define DYNAMICKUBOTOYABETEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidCurveFitting/DynamicKuboToyabe.h"
+#include "MantidCurveFitting/StaticKuboToyabe.h"
+#include "MantidAPI/FunctionDomain1D.h"
+#include "MantidAPI/FunctionValues.h"
+
+using namespace Mantid::Kernel;
+using namespace Mantid::API;
+using namespace Mantid::CurveFitting;
+
+
+class DynamicKuboToyabeTest : public CxxTest::TestSuite
+{
+public:
+
+  void testZFZNDKTFunction()
+  {
+    // Test Dynamic Kubo Toyabe (DKT) for Zero Field (ZF) and Zero Nu (ZN)
+    // Function values must match exactly values from the Static Kubo Toyabe
+    const double asym = 1.0;
+    const double delta = 0.39;
+    const double field = 0;
+    const double nu = 0.0;
+
+    DynamicKuboToyabe dkt; 
+    dkt.initialize();
+    dkt.setParameter("Asym", asym);
+    dkt.setParameter("Delta",delta );
+    dkt.setParameter("Field",field);
+    dkt.setParameter("Nu",   nu);
+
+    StaticKuboToyabe skt;
+    skt.initialize();
+    skt.setParameter("A", asym);
+    skt.setParameter("Delta", delta);
+
+    // define 1d domain of 10 points in interval [0,10]
+    Mantid::API::FunctionDomain1DVector x(0,10,10);
+    Mantid::API::FunctionValues y1(x);
+    Mantid::API::FunctionValues y2(x);
+
+    TS_ASSERT_THROWS_NOTHING(dkt.function(x,y1));
+    TS_ASSERT_THROWS_NOTHING(skt.function(x,y2));
+
+    for(size_t i = 0; i < x.size(); ++i)
+    {
+      TS_ASSERT_DELTA( y1[i], y2[i], 1e-6 );
+    }
+  }
+
+  void testZFDKTFunction()
+  {
+    // Test Dynamic Kubo Toyabe (DKT) for Zero Field (ZF) (non-zero Nu)
+    const double asym = 1.0;
+    const double delta = 0.39;
+    const double field = 0;
+    const double nu = 1.0;
+
+    DynamicKuboToyabe dkt; 
+    dkt.initialize();
+    dkt.setParameter("Asym", asym);
+    dkt.setParameter("Delta",delta );
+    dkt.setParameter("Field",field);
+    dkt.setParameter("Nu",   nu);
+
+    // define 1d domain of 5 points in interval [0,5]
+    Mantid::API::FunctionDomain1DVector x(0,5,5);
+    Mantid::API::FunctionValues y(x);
+
+    TS_ASSERT_THROWS_NOTHING(dkt.function(x,y));
+    TS_ASSERT_DELTA( y[0], 1.000000, 0.000001);
+    TS_ASSERT_DELTA( y[1], 0.850107, 0.000001);
+    TS_ASSERT_DELTA( y[2], 0.625283, 0.000001);
+    TS_ASSERT_DELTA( y[3], 0.449064, 0.000001);
+    TS_ASSERT_DELTA( y[4], 0.323394, 0.000001);
+  }
+
+  void xtestZNDKTFunction()
+  {
+    // Test Dynamic Kubo Toyabe (DKT) for non-zero Field and Zero Nu (ZN)
+    const double asym = 1.0;
+    const double delta = 0.39;
+    const double field = 0.1;
+    const double nu = 0.0;
+
+    DynamicKuboToyabe dkt; 
+    dkt.initialize();
+    dkt.setParameter("Asym", asym);
+    dkt.setParameter("Delta",delta );
+    dkt.setParameter("Field",field);
+    dkt.setParameter("Nu",   nu);
+
+    // define 1d domain of 5 points in interval [0,5]
+    Mantid::API::FunctionDomain1DVector x(0,5,5);
+    Mantid::API::FunctionValues y(x);
+
+    TS_ASSERT_THROWS_NOTHING(dkt.function(x,y));
+
+    TS_ASSERT_DELTA( y[0], 1.000000, 0.000001);
+    TS_ASSERT_DELTA( y[1], 0.784636, 0.000001);
+    TS_ASSERT_DELTA( y[2], 0.353978, 0.000001);
+    TS_ASSERT_DELTA( y[3], 0.073286, 0.000001);
+    TS_ASSERT_DELTA( y[4], 0.055052, 0.000001);
+  }
+
+  void xtestDKTFunction()
+  {
+    // Test Dynamic Kubo Toyabe (DKT) (non-zero Field, non-zero Nu)
+    const double asym = 1.0;
+    const double delta = 0.39;
+    const double field = 0.1;
+    const double nu = 0.5;
+
+    DynamicKuboToyabe dkt; 
+    dkt.initialize();
+    dkt.setParameter("Asym", asym);
+    dkt.setParameter("Delta",delta );
+    dkt.setParameter("Field",field);
+    dkt.setParameter("Nu",   nu);
+
+    // define 1d domain of 5 points in interval [0,5]
+    Mantid::API::FunctionDomain1DVector x(0,5,5);
+    Mantid::API::FunctionValues y(x);
+
+    TS_ASSERT_THROWS_NOTHING(dkt.function(x,y));
+
+    TS_ASSERT_DELTA( y[0], 1.000000, 0.000001);
+    TS_ASSERT_DELTA( y[1], 0.822498, 0.000001);
+    TS_ASSERT_DELTA( y[2], 0.518536, 0.000001);
+    TS_ASSERT_DELTA( y[3], 0.295988, 0.000001);
+    TS_ASSERT_DELTA( y[4], 0.175489, 0.000001);
+  }
+
+
+};
+
+#endif /*DYNAMICKUBOTOYABETEST_H_*/