Re #10499 RRFMuon base implementation

Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RRFMuon.h

@@ -0,0 +1,77 @@
+#ifndef MANTID_ALGORITHM_RRFMUON_H_
+#define MANTID_ALGORITHM_RRFMUON_H_
+
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidAPI/Algorithm.h"
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_fft_complex.h>
+
+namespace Mantid
+{
+  namespace Algorithms
+  {
+    /**Algorithm for calculating Muon decay envelope.
+
+    @author Raquel Alvarez, ISIS, RAL
+    @date 5/12/2014
+
+    Copyright &copy; 2014-12 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 RRFMuon : public API::Algorithm
+    {
+    public:
+      /// Default constructor
+      RRFMuon() {};
+      /// Destructor
+      virtual ~RRFMuon() {};
+      /// Algorithm's name for identification overriding a virtual method
+      virtual const std::string name() const { return "RRFMuon";}
+      ///Summary of algorithm's purpose
+      virtual const std::string summary() const {return "Calculate Muon decay envelop from squashograms and input frequency.";}
+
+      /// Algorithm's version for identification overriding a virtual method
+      virtual int version() const { return 1;}
+      /// Algorithm's category for identification overriding a virtual method
+      virtual const std::string category() const { return "Muon";}
+
+    private:
+
+      /// Initialise the properties
+      void init();
+      /// Run the algorithm
+      void exec();
+      /// Compute the Hilbert transform
+      void hilbert (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs);
+
+      /// Frequency of the oscillations
+      double m_freq;
+      /// Number of data points per histogram
+      size_t m_nData;
+      /// Number of histograms
+      size_t m_nHisto;
+    };
+
+  } // namespace Algorithms
+} // namespace Mantid
+
+#endif /*MANTID_ALGORITHM_RRFMUON_H_*/

Code/Mantid/Framework/Algorithms/src/RRFMuon.cpp

@@ -0,0 +1,159 @@
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidAlgorithms/RRFMuon.h"
+#include <boost/shared_array.hpp>
+
+#define REAL(z,i) ((z)[2*(i)])
+#define IMAG(z,i) ((z)[2*(i)+1])
+
+namespace Mantid
+{
+namespace Algorithms
+{
+
+using namespace Kernel;
+using API::Progress;
+
+// Register the class into the algorithm factory
+DECLARE_ALGORITHM(RRFMuon)
+
+/** Initialisation method. Declares properties to be used in algorithm.
+ *
+ */
+void RRFMuon::init()
+{
+
+  declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("InputWorkspace", "", Direction::Input), 
+    "Name of the input workspace containing the spectra");
+
+  declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace", "", Direction::Output), 
+    "Name of the output workspace containing envelope of the decay" );
+
+  declareProperty(new PropertyWithValue<double>("Frequency", 0, Direction::Input), 
+    "Frequency of the oscillations");
+
+}
+
+/** Executes the algorithm
+ *
+ */
+void RRFMuon::exec()
+{
+	// Get input workspace
+  API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
+
+  // Get number of histograms
+  m_nHisto = inputWs->getNumberHistograms();
+  // Get number of data points
+  m_nData = inputWs->blocksize();
+
+  //if (m_nHisto!=2)
+  //{
+  //  throw std::runtime_error("RRFMuon: InputWorkspace must contain two spectra");
+  //}
+
+
+  // Create output workspace with two spectra
+  API::MatrixWorkspace_sptr outputWs = boost::dynamic_pointer_cast<API::MatrixWorkspace>(
+    API::WorkspaceFactory::Instance().create("Workspace2D", m_nHisto, m_nData+1, m_nData));
+  outputWs->getAxis(0)->unit() = inputWs->getAxis(0)->unit();
+
+  hilbert(inputWs,outputWs);
+
+  // Get input frequency
+  m_freq = getProperty("Frequency");
+
+
+  // Set output workspace
+  setProperty("OutputWorkspace", outputWs);
+
+}
+
+void RRFMuon::hilbert (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs)
+{
+  const int ySize = static_cast<int>(inputWs->blocksize());
+
+  // Alloc memory
+  gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc(ySize);
+  gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc(ySize);
+  boost::shared_array<double> data(new double[2*ySize]);
+
+
+  // Data initialisation
+  for(int i=0; i<ySize; i++)
+  {
+    data[2*i  ] = inputWs->dataY(0)[i]; // even indexes filled with the real part
+    data[2*i+1] = 0; // odd indexes filled with the imaginary part
+  }
+
+  // Forward FFT
+  gsl_fft_complex_forward (data.get(), 1, ySize, wavetable, workspace);
+
+  // Compute Hilbert transform
+  boost::shared_array<double> hilbertData(new double[2*ySize]);
+  for(int i=0; i<ySize; i++)
+  {
+    if ( i < ySize/2 )
+    {
+      hilbertData[2*i  ] =  data[2*i+1];
+      hilbertData[2*i+1] = -data[2*i];
+    }
+    else if ( i > ySize/2 )
+    {
+      hilbertData[2*i  ] = -data[2*i+1];
+      hilbertData[2*i+1] =  data[2*i];
+    }
+    else
+    {
+      hilbertData[2*i  ] = 0;
+      hilbertData[2*i+1] = 0;
+    }
+  }
+  gsl_fft_complex_inverse (hilbertData.get(), 1, ySize, wavetable, workspace);
+
+
+  // Write results into output workspace
+  outputWs->dataX(0) = inputWs->dataX(0);
+  auto oldData = inputWs->dataY(0);
+  for (int i=0; i<ySize; i++)
+  {
+    outputWs->dataY(0)[i] = hilbertData[2*i]; // uncommet to print hilbert transform
+    outputWs->dataY(0)[i] = sqrt( oldData[i]*oldData[i] + hilbertData[2*i]*hilbertData[2*i]);
+  }
+
+  // Free memory
+  gsl_fft_complex_wavetable_free (wavetable);
+  gsl_fft_complex_workspace_free (workspace);
+}
+
+
+//void hilbert(double *source, double *target, double *filt, int npt, int lfilt)
+//{
+//  int i=1;
+//  int l=1;
+//  double yt =0.0;
+//  for (l=1; l<=npt-lfilt+1; l++)
+//  {
+//    yt = 0.0;
+//    for (i=1; i<=lfilt; i++)
+//      yt = yt+source[l+i-1]*filt[lfilt+1-i];
+//    target[l] = yt;
+//  }
+//  /* shifting points */
+//  for (i=1; i<=npt-lfilt; i++)
+//  {
+//    target[i] = 0.5*(target[i]+target[i+1]);
+//  }
+//  for (i=npt-lfilt; i>=1; i--)
+//    target[i+lfilt/2]=target[i];
+//  /* adding zeros */
+//  for (i=1; i<=lfilt/2; i++)
+//  {
+//    target[i] = 0.0;
+//    target[npt+1-i] = 0.0;
+//  }
+//}
+
+}
+}

Code/Mantid/Framework/Algorithms/test/RRFMuonTest.h

@@ -0,0 +1,37 @@
+#ifndef MANTID_ALGORITHMS_RRFMUON_H_
+#define MANTID_ALGORITHMS_RRFMUON_H_
+
+#include <cxxtest/TestSuite.h>
+#include "MantidAlgorithms/RRFMuon.h"
+#include <Poco/File.h>
+#include <stdexcept>
+
+using namespace Mantid::Algorithms;
+using namespace Mantid::API;
+
+class RRFMuonTest : public CxxTest::TestSuite
+{
+public:
+
+  void testName()
+  {
+    TS_ASSERT_EQUALS( rrfMuon.name(), "RRFMuon" );
+  }
+
+  void testCategory()
+  {
+    TS_ASSERT_EQUALS( rrfMuon.category(), "Muon" )
+  }
+
+  void testInit()
+  {
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.initialize() );
+    TS_ASSERT( rrfMuon.isInitialized() )
+  }
+
+
+private:
+  RRFMuon rrfMuon;
+};
+
+#endif /* MANTID_ALGORITHMS_RRFMUON_H_ */

Code/Mantid/docs/source/algorithms/RRFMuon-v1.rst

@@ -0,0 +1,15 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Assuming that the *InputWorkspace* contains two squashograms, the algorithm returns a workspace
+containing the envelope of the muon decay, for a given frequency.
+
+.. categories::