bessel.hh

Go to the documentation of this file.

 
#ifndef besselFormula_hh
#define besselFormula_hh
 
#include "toolbox/sequence.hh"
#include "formula.hh"
// Integrated from git clone https://github.com/valandil/complex_bessel.git
#ifdef HAS_COMPLEX_BESSEL
#include "basics/warnings/push.h"
#include "basics/warnings/ignore_warning_unused_variable.h"
#include "complex_bessel.h"
#include "basics/warnings/pop.h"
#endif
 
namespace concepts {
 
  // ************************************************************** besselJ0 **
 
  Real besselJ0(const Real x);
 
  // ************************************************************** besselY0 **
 
  Real besselY0(const Real x);
 
  // ************************************************************** besselJ1 **
 
  Real besselJ1(const Real x);
 
  // ************************************************************** besselY1 **
 
  Real besselY1(const Real x);
 
  // ************************************************************** besselJn **
 
  Real besselJn(const Real x, const int n);  
 
  Sequence<Real> besselJn(const Real x, const Sequence<int>& n);
 
  // ************************************************************** besselYn **
 
  Real besselYn(const Real x, const int n);
 
  Sequence<Real> besselJn(const Real x, const Sequence<int>& n);  
 
  // ************************************************************ hankel_1_n **
 
  Cmplx hankel_1_n(const Real x, const int n);
 
  Sequence<Cmplx> hankel_1_n(const Real x, const Sequence<int>& n);
 
  // ************************************************************ hankel_2_n **
 
  Cmplx hankel_2_n(const Real x, const int n);
 
  Sequence<Cmplx> hankel_2_n(const Real x, const Sequence<int>& n);
 
  // ****************************************************** hankel_1_deriv_n **
 
  Cmplx hankel_1_deriv_n(const Real x, const int n);
 
  Sequence<Cmplx> hankel_1_deriv_n(const Real x, const Sequence<int>& n);
 
  // ****************************************************** hankel_2_deriv_n **
 
  Cmplx hankel_2_deriv_n(const Real x, const int n);
 
  Sequence<Cmplx> hankel_2_deriv_n(const Real x, const Sequence<int>& n);
 
#ifdef HAS_COMPLEX_BESSEL
 
  // ************************************************************* besselJnu **
 
  Cmplx besselJnu(const Real x, const Real nu);  
 
  Cmplx besselJnu_deriv(const Real x, const Real nu);  
 
  // ************************************************************* besselYnu **
 
  Cmplx besselYnu(const Real x, const Real nu);  
 
  Cmplx besselYnu_deriv(const Real x, const Real nu);  
 
#endif
 
  // *************************************************************** BesselJ **
 
  template<int n>
  class BesselJ : public Formula<Real> {
  public:
    BesselJ(const Real m = 1.0, const Real r0 = 0.0) 
      : m_(m), r0_(r0, 0.0, 0.0) {}
    BesselJ(const Real m, const Real3d r0) 
      : m_(m), r0_(r0) {}
    virtual Real operator() (const Real p,    const Real t = 0.0) const;
    virtual Real operator() (const Real2d& p, const Real t = 0.0) const;
    virtual Real operator() (const Real3d& p, const Real t = 0.0) const;
    Real derivative(const Real p);
    Real derivative(const Real2d& p);
    Real derivative(const Real3d& p);
    
    virtual BesselJ<n>* clone() const {
      return new BesselJ<n>(m_, r0_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const Real m_;
    const Real3d r0_;
  };
  
  template<int n>
  Real BesselJ<n>::operator() (const Real p,    const Real t) const {
    return besselJn(m_*std::abs(p - r0_[0]),n);
  }
 
  template<int n>
  Real BesselJ<n>::operator() (const Real2d& p, const Real t) const {
    return besselJn(m_*(p - Real2d(r0_)).l2(),n);
  }
 
  template<int n>
  Real BesselJ<n>::operator() (const Real3d& p, const Real t) const {
    return besselJn(m_*(p - r0_).l2(),n);
  }
  
  template<int n>
  Real BesselJ<n>::derivative(const Real p) {
    Real x = m_*std::abs(p - r0_[0]);
    return -n/x*besselJn(x,n)+besselJn(x,n-1);
  }
 
  template<int n>
  Real BesselJ<n>::derivative(const Real2d& p) {
    Real x = m_*(p - Real2d(r0_)).l2();
    return -n/x*besselJn(x,n)+besselJn(x,n-1);
  }
 
  template<int n>
  Real BesselJ<n>::derivative(const Real3d& p) {
    Real x = m_*(p - r0_).l2();
    return -n/x*besselJn(x,n)+besselJn(x,n-1);
  }
  
  template<int n>
  std::ostream& BesselJ<n>::info(std::ostream& os) const {
    os << concepts::typeOf(*this)<< "(";
    if (m_ != 1) os << m_ << "*";
    if (r0_.l2() > 0)
      os << "(r - " << r0_ << ")";
    else
      os << "r";
    return os << ")";
  }
 
  // ************************************************************** BesselY **
 
  template<int n>
  class BesselY : public Formula<Real> {
  public:
    BesselY(const Real m = 1.0, const Real r0 = 0.0) 
      : m_(m), r0_(r0, 0.0, 0.0) {}
    BesselY(const Real m, const Real3d r0) 
      : m_(m), r0_(r0) {}
    virtual Real operator() (const Real p,    const Real t = 0.0) const;
    virtual Real operator() (const Real2d& p, const Real t = 0.0) const;
    virtual Real operator() (const Real3d& p, const Real t = 0.0) const;
    Real derivative(const Real p);
    Real derivative(const Real2d& p);
    Real derivative(const Real3d& p);
    
    virtual BesselY<n>* clone() const {
      return new BesselY<n>(m_, r0_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const Real m_;
    const Real3d r0_;
  };
  
  template<int n>
  Real BesselY<n>::operator() (const Real p,    const Real t) const {
    return besselYn(m_*std::abs(p - r0_[0]),n);
  }
 
  template<int n>
  Real BesselY<n>::operator() (const Real2d& p, const Real t) const {
    return besselYn(m_*(p - Real2d(r0_)).l2(),n);
  }
 
  template<int n>
  Real BesselY<n>::operator() (const Real3d& p, const Real t) const {
    return besselYn(m_*(p - r0_).l2(),n);
  }
  
  template<int n>
  Real BesselY<n>::derivative(const Real p) {
    Real x = m_*std::abs(p - r0_[0]);
    return -n/x*besselYn(x,n)+besselYn(x,n-1);
  }
 
  template<int n>
  Real BesselY<n>::derivative(const Real2d& p) {
    Real x = m_*(p - Real2d(r0_)).l2();
    return -n/x*besselYn(x,n)+besselYn(x,n-1);
  }
 
  template<int n>
  Real BesselY<n>::derivative(const Real3d& p) {
    Real x = m_*(p - r0_).l2();
    return -n/x*besselYn(x,n)+besselYn(x,n-1);
  }
  
  template<int n>
  std::ostream& BesselY<n>::info(std::ostream& os) const {
    os << concepts::typeOf(*this)<< "(";
    if (m_ != 1) os << m_ << "*";
    if (r0_.l2() > 0)
      os << "(r - " << r0_ << ")";
    else
      os << "r";
    return os << ")";
  }
 
} // namespace concepts
 
#endif // besselFormula_hh