Interpolator.h

// -*- C++ -*-
//
// Interpolator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_Interpolator_H
#define HERWIG_Interpolator_H
//
// This is the declaration of the Interpolator class.
//
 
#include "ThePEG/Interface/Interfaced.h"
#include <cassert>
 
namespace Herwig {
 
using namespace ThePEG;
 
template <typename ValT, typename ArgT>
class Interpolator: public Interfaced {
 
public:
 
  typedef typename Ptr<Interpolator<ValT,ArgT> >::pointer Ptr;
 
  Interpolator() : _order(3), _copyx(5),_copyfun(5) {}
 
  Interpolator(const vector<ValT> & f, 
         const vector<ArgT> & x, 
         unsigned int order) 
    : _fun(f.size(),0.0),_xval(x.size(),0.0),_order(order),
      _funit(TypeTraits<ValT>::baseunit()), 
      _xunit(TypeTraits<ArgT>::baseunit()),
      _copyx(order+2),_copyfun(order+2) {
    assert(_order>0);
    assert(x.size() == f.size());
    for (size_t i = 0; i < f.size(); ++i) {
      _fun [i] = f[i] / _funit;
      _xval[i] = x[i] / _xunit;
    }
  }
 
  template <size_t N>
  Interpolator(const array<ValT,N> & f, 
         const array<ArgT,N> & x, 
         unsigned int order) 
    : _fun(N,0.0),_xval(N,0.0),_order(order),
      _funit(TypeTraits<ValT>::baseunit()), 
      _xunit(TypeTraits<ArgT>::baseunit()),
      _copyx(order+2),_copyfun(order+2) {
    assert(_order>0);
    for (size_t i = 0; i < N; ++i) {
      _fun [i] = f[i] / _funit;
      _xval[i] = x[i] / _xunit;
    }
  }
 
  Interpolator(size_t size, 
               const double f[], ValT funit,
               const double x[], ArgT xunit,
               unsigned int order) 
    : _fun(size,0.0),_xval(size,0.0),_order(order),
      _funit(funit),_xunit(xunit), _copyx(order+2),_copyfun(order+2) {
    assert(_order>0);
    for (size_t i = 0; i < size; ++i) {
      _fun [i] = f[i];
      _xval[i] = x[i];
    }
  }
 
  ValT operator () (ArgT) const;
  typedef ValT ValType;
  typedef ArgT ArgType;
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const { return new_ptr(*this); }
 
  virtual IBPtr fullclone() const { return new_ptr(*this); }
 
private:
 
  Interpolator & operator=(const Interpolator &) = delete;
  
private:
  
  vector<double> _fun;
 
  vector<double> _xval;
 
  unsigned int _order;
 
  ValT _funit;
 
  ArgT _xunit;
 
  mutable vector<double> _copyx;
 
  mutable vector<double> _copyfun;
 
};
 
template <typename ValT, typename ArgT>
inline typename Interpolator<ValT,ArgT>::Ptr
make_InterpolatorPtr(size_t size, 
         const double f[], ValT funit,
         const double x[], ArgT xunit,
         unsigned int order)
{
  return new_ptr(Interpolator<ValT,ArgT>(size,
           f,funit,
           x,xunit,
           order));
}
 
template <typename ValT, typename ArgT>
inline typename Interpolator<ValT,ArgT>::Ptr
make_InterpolatorPtr(const typename std::vector<ValT> & f, 
         const typename std::vector<ArgT> & x, 
         unsigned int order)
{
  return new_ptr(Interpolator<ValT,ArgT>(f,x,order));
}
 
template <typename ValT, typename ArgT, size_t N>
inline typename Interpolator<ValT,ArgT>::Ptr
make_InterpolatorPtr(const typename std::array<ValT,N> & f, 
         const typename std::array<ArgT,N> & x, 
         unsigned int order)
{
  return new_ptr(Interpolator<ValT,ArgT>(f,x,order));
}
 
}
 
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "Interpolator.tcc"
#endif
 
#endif /* HERWIG_Interpolator_H */