math_functions.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2014  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2014  Peter Rucz <rucz@hit.bme.hu>
//
// This program 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.
//
// This program 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/>.
 
#include <boost/math/constants/constants.hpp>
 
#ifndef MATH_FUNCTIONS_HPP_INCLUDED
#define MATH_FUNCTIONS_HPP_INCLUDED
 
#include <cmath>
#include <complex>
#include <iostream>
#include <stdexcept>
 
namespace NiHu
{
 
template <typename T> 
int sgn(T val) 
{
    return (T(0) < val) - (val < T(0));
}
    
template <class T>
static T sinc(T const &x)
{
    if (std::abs(x) > 1e-3)
        return std::sin(x) / x;
    else
        return 1. - x*x/6. * (1. - x*x/20.);
}
 
namespace bessel
{
    template <class T>
    struct make_complex 
    { 
        typedef std::complex<T> type;
    };
 
    template <class T>
    struct make_complex<std::complex<T> > 
    {
        typedef std::complex<T> type;
    };
 
    double const large_lim(7.);
 
    template <class T>
    void mag_arg_large(int nu, T const &z, T &u, T &phi)
    {
        using namespace boost::math::double_constants;
 
        double mag[3][7] = {
            {
                1.,
                -6.25e-2,
                1.03515625e-1,
                -5.428466796875e-1,
                5.8486995697021484375,
                -1.06886793971061706543e2,
                2.96814293784275650978e3
            },
            {
                1.,
                1.875e-1,
                -1.93359375e-1,
                8.052978515625e-1,
                -7.7399539947509765625,
                1.32761824250221252441e2,
                -3.54330366536602377892e3
            },
            {
                1,
                9.3750000000e-01,
                7.9101562500e-01,
                -3.0487060547e+00,
                1.9199895859e+01,
                -2.5916166008e+02,
                6.0841126316e+03
            }
        };
        double arg[3][7] = {
            {
                -1.25e-1,
                6.51041666666666712926e-2,
                -2.09570312499999994449e-1, 
                1.63806588309151779370,
                -2.34751277499728736586e1,
                5.35640519510615945364e2,
                -1.78372796889474739146e4
            },
            {
                3.75e-1,
                -1.640625e-1,
                3.70898437500000011102e-1,
                -2.36939784458705338110,
                3.06240119934082031250e1,
                -6.59185221823778988437e2,
                2.11563140455278044101e4
            },
            {
                1.8750000000e0,
                -3.5156250000e-1,
                -1.4501953125e0,
                8.3074297224e0,
                -7.1739864349e1,
                1.2458673851e3,
                -3.5571449717e4             
            }
        };
 
        u = phi = T(0.);
        auto q(1./z/z);
        for (int m = sizeof(mag[0])/sizeof(double)-1; m >= 0; --m)
        {
            u = q * u + mag[nu][m];
            phi = q * phi + arg[nu][m];
        }
        phi = phi/z + z - (nu/2.+.25)*pi;
    }
 
    template <int nu, class T>
    T J_small(T const &z)
    {
        // upper limit for 1e-8 error
        int N = (int)(3+2.*std::abs(z));
 
        T q = z/2.;
 
        T res(1.), q2(q*q);
        for (int k = N; k > 0; --k)
            res = 1. - q2*(1./(k*(k+nu)))*res;
        for (int k = 1; k <= nu; ++k)
            res *= q/T(k);
        return res;
    }
 
    template <int nu, class T>
    T J_large(T const &z)
    {
        using boost::math::double_constants::pi;
 
        static_assert(nu >= 0 && nu <= 2 , "unimplemented Bessel J order");
 
        T mag, arg;
        if (std::real(z) < 0)
        {
            double const C(nu%2==0 ? 1. : -1);
            mag_arg_large(nu, -z, mag, arg);
            return std::sqrt(2./(pi * -z)) * mag * std::cos(arg) * C;
        }
        else
        {
            mag_arg_large(nu, z, mag, arg);
            return std::sqrt(2./(pi * z)) * mag * std::cos(arg);
        }
    }
 
    template <int nu, class T>
    T J(T const &z)
    {
         return std::abs(z) < large_lim ? J_small<nu>(z) : J_large<nu>(z);
    }
 
 
    template <int nu>
    std::complex<double> Y_small(std::complex<double> const &z)
    {
        using boost::math::double_constants::pi;
        using boost::math::double_constants::euler;
        
        // upper limit for 1e-8 error
        int const N = (int)(4+2.*std::abs(z));
 
        std::complex<double> q(z/2.0), q2(q*q);
        std::complex<double> first(2.0*J_small<nu>(z)*(std::log(q) + euler));
        std::complex<double> second;
        switch (nu)
        {
            case 0: second = 0.; break;
            case 1: second = -1./q; break;
            case 2: second = -1. -1./q2; break;
        }
 
        std::complex<double> third(0.0);
        std::complex<double> q2pow(-std::pow(q, nu));
        
        double a(0.);
        for (int k = 1; k <= nu; ++k)
            a += 1./k;
        
        double div = 1.;
        for (int k = 1; k <= nu; ++k)
            div /= k;
        
        for (int k = 0; k < N; ++k)
        {
            third += div*q2pow*a;
 
            div /= -(k+1)*(k+nu+1);
            q2pow *= q2;
            a += 1./(k+1.) + 1./(k+nu+1.);
        }
 
        return 1./pi * (first + second + third);
    }
 
    template <int nu>
    std::complex<double> Y_large(std::complex<double> const &z)
    {
        using boost::math::double_constants::pi;
        
        static_assert(nu >= 0 || nu <= 2, "unimplemented Bessel Y order");
        
        std::complex<double> mag, arg;
        if (std::real(z) < 0)
        {
            double const C1(nu%2 == 0 ? 1. : -1.);
            std::complex<double> const C2(0., std::imag(z) < 0 ? -2. : 2.);
            std::complex<double> const MAG(std::sqrt(C2*C2+1.));
            std::complex<double> const ARG(std::atan(1./C2));
            mag_arg_large(nu, -z, mag, arg);
            return std::sqrt(2./(pi * -z)) * mag * MAG * std::cos(arg-ARG) * C1;
        }
        else
        {
            mag_arg_large(nu, z, mag, arg);
            return std::sqrt(2./(pi * z)) * mag * std::sin(arg);
        }
    }
 
 
    template <int nu>
    std::complex<double> Y(std::complex<double> const &z)
    {
        return std::abs(z) < large_lim ? Y_small<nu>(z) : Y_large<nu>(z);
    }
 
    template <int nu, int kind, class T>
    typename make_complex<T>::type H_large(T const &z)
    {
        using namespace std::complex_literals;
        using boost::math::double_constants::pi;
 
        static_assert((kind == 1) || (kind == 2), "invalid kind argument of bessel::H");
        double const C = (kind == 2 ? -1. : 1.);
 
        typename make_complex<T>::type mag, arg;
 
        if (std::real(z) < 0 && std::abs(std::imag(z)) < 8.)
        {
            double const C1(nu%2 == 0 ? 1. : -1.);
 
            mag_arg_large(nu, -z, mag, arg);
 
            double sgn = std::imag(z) < 0 ? -1. : 1.;
            double MAG(-sgn*std::sqrt(3.));
            std::complex<double> const ARG(std::atan(-sgn*.5i));
 
            return std::sqrt(2./(pi * -z)) * mag * (std::cos(arg) + C * MAG * std::cos(arg-ARG)) * C1;
        }
 
        mag_arg_large(nu, z, mag, arg);
        
        return std::sqrt(2./pi/z) * mag * std::exp(C*1.i*arg);
    }
 
    template <int nu, int kind, class T>
    typename make_complex<T>::type H(T const &z)
    {
        using namespace std::complex_literals;
        static_assert(kind == 1 || kind == 2, "invalid kind argument of bessel::H");
        double const C = (kind == 2 ? -1. : 1.);
        if (std::abs(z) < large_lim)
            return J_small<nu>(z) + C * 1.i * Y_small<nu>(z);
        else
            return H_large<nu, kind>(z);
    }
 
    template <int nu, class T>
    typename make_complex<T>::type K(T const &z);
} // end of namespace bessel
 
 
template <class T = double>
T chebroot(size_t n, size_t i)
{
    using boost::math::double_constants::pi;
    return std::cos(pi / 2. * (2 * (n - i) - 1) / T(n));
}
 
 
template <class I>
I Ipow(I base, I exp)
{
    I res = 1;
    for (unsigned i = 0; i < exp; ++i)
        res *= base;
    return res;
}
 
 
template < unsigned Base, unsigned Exp >
struct IpowC
{
    static unsigned const value = IpowC<Base, Exp - 1>::value * Base;
};
 
template <unsigned Base>
struct IpowC < Base, 0 >
{
    static unsigned const value = 1U;
};
 
template <unsigned Exp>
struct IpowC < 0, Exp >
{
    static unsigned const value = 0U;
};
 
template <>
struct IpowC < 0, 0 >
{
    // undefined
};
 
 
} // end of namespace NiHu
 
#endif /* MATH_FUNCTIONS_HPP_INCLUDED */