Add a BoseEinsteinDistribution class. Refs #5397

mantidproject · Jul 12, 2012 · 7ed5a88 · 7ed5a88
1 parent 4f30e1e
commit 7ed5a88
Show file tree

Hide file tree

Showing 4 changed files with 257 additions and 0 deletions.
diff --git a/Code/Mantid/Framework/Kernel/CMakeLists.txt b/Code/Mantid/Framework/Kernel/CMakeLists.txt
@@ -35,6 +35,7 @@ set ( SRC_FILES
 	src/LogFilter.cpp
 	src/LogParser.cpp
 	src/Logger.cpp
+	src/Math/Distributions/BoseEinsteinDistribution.cpp
 	src/MagneticIon.cpp
 	src/MandatoryValidator.cpp
 	src/MantidVersion.cpp
@@ -146,6 +147,7 @@ set ( INC_FILES
 	inc/MantidKernel/LogFilter.h
 	inc/MantidKernel/LogParser.h
 	inc/MantidKernel/Logger.h
+	inc/MantidKernel/Math/Distributions/BoseEinsteinDistribution.h
 	inc/MantidKernel/MRUList.h
 	inc/MantidKernel/MagneticIon.h
 	inc/MantidKernel/MandatoryValidator.h
@@ -217,6 +219,7 @@ set ( TEST_FILES
 	test/AtomTest.h
 	test/BinFinderTest.h
 	test/BinaryFileTest.h
+	test/BoseEinsteinDistributionTest.h
 	test/BoundedValidatorTest.h
 	test/CPUTimerTest.h
 	test/CacheTest.h

diff --git a/Code/Mantid/Framework/Kernel/inc/MantidKernel/Math/Distributions/BoseEinsteinDistribution.h b/Code/Mantid/Framework/Kernel/inc/MantidKernel/Math/Distributions/BoseEinsteinDistribution.h
@@ -0,0 +1,51 @@
+#ifndef MANTID_KERNEL_BOSEEINSTEINDISTRIBUTION_H_
+#define MANTID_KERNEL_BOSEEINSTEINDISTRIBUTION_H_
+/**
+    Copyright &copy; 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://svn.mantidproject.org/mantid/trunk/Code/Mantid>.
+    Code Documentation is available at: <http://doxygen.mantidproject.org>
+*/
+#include "MantidKernel/DllConfig.h"
+
+namespace Mantid
+{
+  namespace Kernel
+  {
+    namespace Math
+    {
+      /**
+       * Defines a static class for computing the coefficient from a
+       * Bose-Einstein distribution for a given energy in meV & temperature in Kelvin
+       */
+      class MANTID_KERNEL_DLL BoseEinsteinDistribution
+      {
+      public:
+        /// Calculate the expected number of particles in an energy state at a given temperature
+        /// for a degenerate distribution with zero chemical potential
+        static double n(const double energy, const double temperature);
+        /// Calculate the \f$(n+1)\epsilon\f$ for a degenerate distribution with zero chemical potential
+        /// where n is the Bose-Einstein distribution
+        static double np1Eps(const double energy, const double temperature);
+      };
+
+    }
+  }
+}
+
+#endif /* MANTID_KERNEL_BOSEEINSTEINDISTRIBUTION_H_ */
diff --git a/Code/Mantid/Framework/Kernel/src/Math/Distributions/BoseEinsteinDistribution.cpp b/Code/Mantid/Framework/Kernel/src/Math/Distributions/BoseEinsteinDistribution.cpp
@@ -0,0 +1,113 @@
+#include "MantidKernel/Math/Distributions/BoseEinsteinDistribution.h"
+#include "MantidKernel/Exception.h"
+#include "MantidKernel/PhysicalConstants.h"
+
+#include <boost/lexical_cast.hpp>
+
+namespace Mantid {  namespace Kernel
+{
+  namespace Math
+  {
+    namespace
+    {
+      /// Tolerance to consider value zero
+      double ZERO_EPS = 1e-12;
+      /// Forward declaration of helper to raise domain error
+      void throwDomainError(const std::string & msg, const double value);
+      /// Forward declaration of helper to calculate distribution
+      double yOver1MinusExpY(const double);
+    }
+
+    /**
+     * Calculate the expected number of particles in an energy state at a given temperature
+     * for a Bose-Einstein distribution: \f$1/[exp(\epsilon/k_BT) - 1]\f$. Note this function
+     * is not well behaved for y=0. @see np1Eps
+     * @param energy :: The value of the energy of the state in meV
+     * @param temperature :: The temperature in Kelvin
+     * @return The value of the distribution
+     */
+    double BoseEinsteinDistribution::n(const double energy, const double temperature)
+    {
+      static const double kbT = PhysicalConstants::BoltzmannConstant*temperature;
+      if(std::abs(temperature) < ZERO_EPS)
+        throwDomainError("BoseEinsteinDistribution::n - Temperature very small, function not well behaved", temperature);
+
+      const double beta = energy/kbT;
+      if(std::abs(beta) < ZERO_EPS)
+        throwDomainError("BoseEinsteinDistribution::n - Exponent very small, function not well-behaved", beta);
+
+      return 1.0/(std::exp(beta) - 1.0);
+    }
+
+    /**
+     * Calculate the expected number of particles in an energy state at a given temperature
+     * for a Bose-Einstein distribution: \f$1/[exp(\epsilon/k_BT) - 1]\f$. Deals with edge cases where
+     *    exponent = 0 using a taylor series
+     *    large & negative exponent by multiplying by \f$exp(-y)\f$
+     * @param energy :: The value of the energy of the state in meV
+     * @param temperature :: The temperature in Kelvin
+     * @return The value of the distribution
+     */
+    double BoseEinsteinDistribution::np1Eps(const double energy, const double temperature)
+    {
+      if(temperature < ZERO_EPS)
+      {
+        if(energy < 0.0) return 0.0;
+        else return energy;
+      }
+      else
+      {
+        const double kBT = (PhysicalConstants::BoltzmannConstant*temperature);
+        return  kBT * yOver1MinusExpY(energy/kBT);
+      }
+    }
+
+    //---------------------------------------------------------------------
+    // Non-member functions
+    //---------------------------------------------------------------------
+    namespace
+    {
+      /**
+       * Raises a domain error with the given error message
+       * @param msg :: Error message
+       * @param value :: Appends string Value=value after
+       */
+      void throwDomainError(const std::string & msg, const double value)
+      {
+        throw std::domain_error(msg + " Value=" + boost::lexical_cast<std::string>(value));
+      }
+
+      /**
+       * Calculate y/(1-exp(-y)) dealing with edges cases at
+       *   y = 0 : using Taylor series
+       *   large negative y: multiplying top & bottom by exp(-y)
+       * @param y The value of y
+       * @return The value of y/(1-exp(-y))
+       */
+      double yOver1MinusExpY(const double y)
+      {
+        double magnitudeY = std::abs(y);
+        if(magnitudeY > 0.1 )
+        {
+          const double expMinusY = std::exp(-magnitudeY);
+          return magnitudeY / (1.0 - expMinusY);
+        }
+        else
+        {
+          // Taylor series coefficients
+          static const double by2(0.5), by6(1.0/6.0), by60(1.0/60.0),
+              by42(1.0/42.0), by40(1.0/40.0);
+          const double ysqr = y*y;
+          return 1.0 + by2*y*( 1.0 + by6*y*( 1.0 - by60*ysqr*(1.0-by42*ysqr*(1.0 - by40*ysqr ))));
+        }
+      }
+    } //end anonyomous
+
+  } // end Math
+}} // end Mantid::Kernel
+
+
+
+
+
+
diff --git a/Code/Mantid/Framework/Kernel/test/BoseEinsteinDistributionTest.h b/Code/Mantid/Framework/Kernel/test/BoseEinsteinDistributionTest.h
@@ -0,0 +1,90 @@
+#ifndef MANTID_KERNEL_BOSEEINSTEINDISTRIBUTIONTEST_H_
+#define MANTID_KERNEL_BOSEEINSTEINDISTRIBUTIONTEST_H_
+
+#include "MantidKernel/Math/Distributions/BoseEinsteinDistribution.h"
+#include "MantidKernel/PhysicalConstants.h"
+
+#include <cxxtest/TestSuite.h>
+
+#include <iostream>
+#include <iomanip>
+
+class BoseEinsteinDistributionTest : public CxxTest::TestSuite
+{
+public:
+
+  void test_Standard_Distribution_Gives_Correct_Value_Away_From_Edge()
+  {
+    using namespace Mantid::Kernel::Math;
+    const double energy = 30.0;
+    const double temperature = 35.0;
+
+    TS_ASSERT_DELTA(BoseEinsteinDistribution::n(energy, temperature), 0.000047886213, 1e-12);
+  }
+
+  void test_Standard_Distribution_Throws_When_Energy_Or_Temperature_Is_Zero()
+  {
+    using namespace Mantid::Kernel::Math;
+    const double energy = 0.0;
+    const double temperature = 35.0;
+
+    TS_ASSERT_THROWS(BoseEinsteinDistribution::n(energy, temperature), std::domain_error);
+    TS_ASSERT_THROWS(BoseEinsteinDistribution::n(temperature, energy), std::domain_error);
+  }
+
+  void test_np1Eps_Is_Returns_Energy_When_Temp_Is_Negative_And_Energy_Positive()
+  {
+    using namespace Mantid::Kernel::Math;
+    const double energy = 200;
+    const double temperature = -35.0;
+
+    const double expected = energy;
+    TS_ASSERT_DELTA(BoseEinsteinDistribution::np1Eps(energy, temperature), expected, 1e-12);
+  }
+
+  void test_np1Eps_Returns_kbT_When_Answer_When_Exponent_Is_Zero()
+  {
+    using namespace Mantid::Kernel::Math;
+    const double energy = 0.0;
+    const double temperature = 35.0;
+
+    const double expected = Mantid::PhysicalConstants::BoltzmannConstant*temperature;
+    TS_ASSERT_DELTA(BoseEinsteinDistribution::np1Eps(energy, temperature), expected, 1e-12);
+  }
+
+  void test_np1Eps_Is_Returns_Zero_When_Temp_Is_Negative_And_Energy_Negative()
+  {
+    using namespace Mantid::Kernel::Math;
+    const double energy = -200;
+    const double temperature = -35.0;
+
+    const double expected = 0.0;
+    TS_ASSERT_DELTA(BoseEinsteinDistribution::np1Eps(energy, temperature), expected, 1e-12);
+  }
+
+  void test_np1Eps_Is_Well_Behaved_When_Exponent_Is_Larger_Than_Point1()
+  {
+    using namespace Mantid::Kernel::Math;
+    const double energy = 20;
+    const double temperature = 29.0;
+
+    const double expected = 20.006690611537;
+
+    TS_ASSERT_DELTA(BoseEinsteinDistribution::np1Eps(energy, temperature), expected, 1e-12);
+  }
+
+  void test_np1Eps_Is_Well_Behaved_When_Abs_Exponent_Is_Larger_Than_Point1_But_Large_And_Negative()
+  {
+    using namespace Mantid::Kernel::Math;
+    const double energy = -20;
+    const double temperature = 35.0;
+
+    const double expected = 20.026407635389;
+
+
+    TS_ASSERT_DELTA(BoseEinsteinDistribution::np1Eps(energy, temperature), expected, 1e-12);
+  }
+
+};
+
+#endif /* MANTID_KERNEL_BOSEEINSTEINDISTRIBUTIONTEST_H_ */