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(T n, T a): m_n (n), m_a(a) {}
 
    DG_DEVICE
    T operator()(T x) const {
        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*x*m_a*m_a)*exp( x*m_a);
        return pow(-x*m_a,m_n)/tgamma(m_n+1)*exp(x*m_a);
    }
    DG_DEVICE
    T operator()(T x, T y) const { return this->operator()(x*y); }
 
    private:
    T m_n, m_a;
};
 
}//namespace mat
}//namespace dg