functors.h

Go to the documentation of this file.

#pragma once
#include <boost/math/special_functions.hpp>
#include "../dg/backend/config.h" // DG_DEVICE
namespace dg {
namespace mat {
template < class T = double>
struct BESSELI0
{
    BESSELI0( ) {}
    T operator() ( T x) const
    {
        return boost::math::cyl_bessel_i(0, x);
    }
};
 
template< class T = double >
struct BesselJ
{
    BesselJ(unsigned n): m_n(n) {}
 
    T operator() ( T x) const
    {
        return boost::math::cyl_bessel_j(m_n, x);
    }
    private:
        unsigned m_n;
};
 
template< class T = double >
struct LaguerreL
{
    LaguerreL(unsigned n): m_n(n) {}
 
    T operator() ( T x) const
    {
        return boost::math::laguerre(m_n, x);
    }
    private:
        unsigned m_n;
};
 
 
template < class T = double>
struct GAMMA0
{
    GAMMA0( ) {}
    T operator() ( T x) const
    {
        return exp(x)*boost::math::cyl_bessel_i(0, x);
    }
};
 
 
template<class T>
T phi1( T x){
    if ( fabs( x) < 1e-16 )
        return 1.;
    if ( fabs(x) < 1)
        return expm1( x)/x;
    return (exp( x) - 1)/x;
}
 
template<class T>
T phi2( T x){
    if ( fabs( x) < 1e-16 )
        return 1./2;
    if ( fabs(x) < 1e-2)
        return 1./2.+x*(1./6.+x*(1./24. + x/120.));
    return ((exp( x) - 1)/x - 1)/x;
}
 
template<class T>
T phi3( T x){
    if ( fabs( x) < 1e-16 )
        return 1./6.;
    if ( fabs(x) < 1e-2)
        return 1./6. + x*(1./24.+x*(1./120. +x/720.));
    return (((exp( x) - 1)/x - 1)/x-1./2.)/x;
}
template<class T>
T phi4( T x){
    if ( fabs( x) < 1e-16 )
        return 1./24.;
    if ( fabs(x) < 1e-2)
        return 1./24. + x*(1./120.+x*(1./720. + x/5040));
    return ((((exp( x) - 1)/x - 1)/x-1./2.)/x-1./6.)/x;
}
 
template<class T = double>
struct GyrolagK
{
    GyrolagK(unsigned n, T a = T(1)): m_n (n), m_a(a) {}
 
    DG_DEVICE
    T operator()(T x) const {
        // Factorials overflow int32 for n>12
        static constexpr unsigned fact[] = {1, 1, 2, 6, 24, 120, 720,
            5040, 40320, 362880, 3628800, 39916800, 479001600};
 
        if( m_n == 0)
            return exp(-x*m_a); // faster to evaluate than tgamma and pow ...
        if( m_n == 1)
            return (x*m_a)*exp( -x*m_a);
        if( m_n == 2)
            return 0.5*(x*m_a)*(x*m_a)*exp( -x*m_a);
        if( m_n == 3)
            return (x*m_a)*(x*m_a)*(x*m_a)*exp( -x*m_a)/6.0;
        return pow(x*m_a,m_n)/fact[m_n]*exp(-x*m_a);
    }
    DG_DEVICE
    T operator()(T x, T y) const { return this->operator()(x*y); }
 
    private:
    unsigned m_n;
    T m_a;
};
 
template<class T = double>
struct DGyrolagK
{
    DGyrolagK(unsigned n, T a = T(1)): m_n (n), m_a(a), m_kn(n,a), m_knm(n-1,a) {}
    T  DG_DEVICE operator()(T x) const {
        if( m_n == 0)
            return-m_a*exp(-x*m_a);
        return m_a*(m_knm(x) - m_kn(x));
    }
    T  DG_DEVICE operator()(T x, T y) const { return this->operator()(x*y); }
 
    private:
    unsigned m_n;
    T m_a;
    GyrolagK<T> m_kn, m_knm;
};
 
}//namespace mat
}//namespace dg