special.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 *                                                                           *
 *                       / /   / /    __    \ \  / /                         *
 *                      / /__ / /   / _ \    \ \/ /                          *
 *                     /  ___  /   | |/_/    / /\ \                          *
 *                    / /   / /    \_\      / /  \ \                         *
 *                                                                           *
 *                         Jakub Benda (c) 2014                              *
 *                     Charles University in Prague                          *
 *                                                                           *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef HEX_SPECF
#define HEX_SPECF

#include <cmath>

#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>

#define fac(n) gsl_sf_fact(n)

#define sqr(x) (gsl_sf_pow_int((x),2))

#include "arrays.h"
#include "complex.h"
#include "misc.h"

#define ALPHA_PLUS(n,l)  (sqrt(((n)-(l))*((n)+(l)+1)))
#define ALPHA_MINUS(n,l) (sqrt(((n)-(l)-1)*((n)+(l))))

#define DELTA(a,b)       ((a) == (b) ? 1 : 0)

namespace special
{

namespace constant
{

const double e              = M_E; // e

const double pi             = M_PI; // π
const double pi_half        = M_PI_2; // π/2
const double pi_quart       = M_PI_4; // π/4
const double pi_sqrt        = M_SQRTPI; // √π
const double pi_inv         = M_1_PI; // 1/π
const double pi_two_inv     = M_2_PI; // 2/π
const double pi_two_invsqrt = M_2_SQRTPI; // 2/√π

const double sqrt_two       = M_SQRT2; // √2
const double sqrt_three     = M_SQRT3; // √3
const double sqrt_half      = M_SQRT1_2; // 1/√2

const double log2e          = M_LOG2E;  // log₂ e
const double log10e         = M_LOG10E; // log e
const double ln2            = M_LN2; // ln 2
const double ln10           = M_LN10; // ln 10
const double lnpi           = M_LNPI; // ln π

const double euler          = M_EULER; // γ

// other special numbers
const double Inf = std::numeric_limits<double>::infinity();
const double Nan = std::numeric_limits<double>::quiet_NaN();

}; // end of namespace "special::constant"

namespace integral
{

template <class T> T trapz (double h, NumberArray<T> const & y)
{
    return h * (sum(y) - 0.5 * (y.front() + y.back()));
}

template <class T> T trapz (NumberArray<T> const & x, NumberArray<T> const & y)
{
    assert (x.size() == y.size());
    
    T sum = 0;
    for (unsigned i = 1; i < x.size(); i++)
        sum += 0.5 * (y[i] + y[i-1]) * (x[i] - x[i-1]);
    
    return sum;
}

template <class T> T simpson (double h, NumberArray<T> const & y)
{
    if (y.size() % 2 != 0)
        throw exception ("You need to use even number of grid points for Simpson integration.");
    
    T sum1 = 0, sum2 = 0;
    
    for (int i = 1; i < y.size(); i += 2)
        sum1 += y[i];
    for (int i = 2; i < y.size() - 1; i += 2)
        sum2 += y[i];
    
    return h * (y.front() + 4. * sum1 + 2. * sum2 + y.back()) / 3.;
}

template <class T> T pow_exp_hyperg1F1 (T a, T b, T c, T u, T v, T x, double epsrel = 1e-10, unsigned maxiter = 1000)
{
    // last added term and up-to-now sum of terms
    T term = 1. / (a + 1.);
    T bfactor = 1.;
    T suma = term;
    
    // auxiliary variables
    T c_over_b = c / b;
    
    // the outer sum
    for (unsigned n = 1; std::abs(term) > epsrel * std::abs(suma); n++)
    {
        // compute the inner sum
        T iterm = 1, isum = 0;
        for (unsigned k = 0; k <= n; k++)
        {
            // update the inner sum
            isum += iterm;
            
            // update term for the next k
            iterm *= (u + T(k)) / (v + T(k)) * ((n - k) / (k + 1.)) * c_over_b;
        }
        
        // update the term and the outer sum
        bfactor *= b * x / T(n);
        term = bfactor * isum / (a + T(n + 1));
        suma += term;
        
        // check if we run out of allowed iterations
        if (n == maxiter)
        {
            throw exception
            (
                "Maximal number of iterations (%d) reached in pow_exp_hyperg1F1.",
                maxiter
            );
        }
    }
    
    // return the result
    return suma;
}

template <class T> NumberArray<T> romberg (const ArrayView<T> y)
{
    NumberArray<T> z = y;
    unsigned N = y.size();
    
    T scale = 1;
    for (unsigned log4scale = 1; log4scale < N; log4scale++)
    {
        scale *= 4;
        for (unsigned j = log4scale; j < N; j++)
        {
            z[j] = (scale * z[j] - z[j-1]) / (scale - T(1));
        }
    }
    
    return z;
}

}; // end of namespace "special::integral"

int coulomb_zeros (double eta, int L, int nzeros, double * zeros, double epsrel = 1e-8);

std::vector<std::vector<int>> FdB_partition (int n);

template <class T, class OuterFunctionDerivative, class InnerFunctionDerivative> T chain_rule
(
    OuterFunctionDerivative Df,
    InnerFunctionDerivative Dg,
    int n,
    double x
)
{
    // no derivative : evaluate the function
    if (n == 0)
        return Df(0,Dg(0,x));
    
    // first and higher derivative : use the Faa di Bruno formula
    T suma = 0, term = 0;
    for (std::vector<int> & counts : FdB_partition(n))
    {
        // evaluate derivative of "f"
        if ((term = Df(std::accumulate(counts.begin(), counts.end(), 0),Dg(0,x))) == 0)
            continue;
        
        // evaluate all derivatives of "g"
        for (int j = 1; j <= n; j++)
        {
            term *= std::pow(Dg(j,x) / gsl_sf_fact(j), counts[j-1]) / gsl_sf_fact(counts[j-1]);
            
            if (term == 0)
                break;
        }
        
        // update the sum
        suma += term;
    }
    return gsl_sf_fact(n) * suma;
}

}; // end of namespace "special"

Complex hydro_P(unsigned n, unsigned l, Complex z);

Complex dhydro_P(unsigned n, unsigned l, Complex z);

//
// Ricatti-Bessel functions.
//

Complex ric_j(int n, Complex z);

cArray ric_jv(int lmax, Complex z);

Complex dric_j(int n, Complex z);

cArray dric_jv(int lmax, Complex z);

inline double ric_j (int n, double x)
{
    return x * gsl_sf_bessel_jl(n,x);
}

inline double ric_n (int n, double x)
{
    return x * gsl_sf_bessel_yl(n,x);
}

//
// Modified Bessel function of the first kind.
//

inline double sph_i_scaled (int n, double x)
{
    return gsl_sf_bessel_il_scaled(n,x);
}

inline double ric_i_scaled (int n, double x)
{
    return x * gsl_sf_bessel_il_scaled(n,x);
}

inline double sph_i (int n, double x)
{
    return exp(x) * sph_i_scaled(n,x);
}

inline double ric_i (int n, double x)
{
    return exp(x)*ric_i_scaled(n,x);
}


//
// Modified Bessel function of the second kind.
//

inline double sph_k_scaled (int n, double x)
{
    return gsl_sf_bessel_kl_scaled(n,x) * M_2_PI;
}

inline double ric_k_scaled (int n, double x)
{
    return x * sph_k_scaled(n,x);
}

inline double sph_k (int n, double x)
{
    return exp(-x) * sph_k_scaled(n,x);
}

inline double ric_k (int n, double x)
{
    return exp(-x) * ric_k_scaled(n,x);
}

//
// Ricatti-Hankel functions.
//

inline Complex ric_h_plus (int n, double x)
{
    return Complex(ric_j(n,x), ric_n(n,x));
}

//
// Other special functions.
//

Complex sphY(int l, int m, double theta, double phi);

Complex sphBiY(int l1, int l2, int L, int M, double theta1, double phi1, double theta2, double phi2);

template <typename T> class RadialFunction
{
public:
    
    virtual T operator() (double x) const = 0;
};

inline double dric_j(int n, double x)
{
    if (n == 0)
        return cos(x);
    
    return gsl_sf_bessel_jl(n,x) + (n * ric_j(n-1,x) - (n+1) * ric_j(n+1,x)) / (2*n+1);
}

inline double dric_n(int n, double x)
{
    if (n == 0)
        return sin(x);
    
    return gsl_sf_bessel_yl(n,x) + (n * ric_n(n-1,x) - (n+1) * ric_n(n+1,x)) / (2*n+1);
}

inline double dric_i(int n, double x)
{
    if (n == 0)
        return cosh(x);
    
    return exp(x) * (sph_i_scaled(n,x) + (n * ric_i_scaled(n-1,x) + (n+1) * ric_i_scaled(n+1,x)) / (2*n+1));
}

inline double dric_k(int n, double x)
{
    if (n == 0)
        return -exp(-x);
    
    return exp(-x) * (sph_k_scaled(n,x) - (n * ric_k_scaled(n-1,x) + (n+1) * ric_k_scaled(n+1,x)) / (2*n+1));
}

inline double dric_i_scaled(int n, double x)
{
    if (n == 0)
        return exp(-2*x);
    
    return sph_i_scaled(n,x) - ric_i_scaled(n,x) + (n * ric_i_scaled(n-1,x) + (n+1) * ric_i_scaled(n+1,x)) / (2*n+1);
}

inline double dric_k_scaled(int n, double x)
{
    if (n == 0)
        return 0;
    
    return sph_k_scaled(n,x) + ric_k_scaled(n,x) - (n * ric_k_scaled(n-1,x) + (n+1) * ric_k_scaled(n+1,x)) / (2*n+1);
}

inline Complex dric_h_plus (int n, double x)
{
    return Complex(dric_j(n,x),dric_n(n,x));
}

int coul_F_michel(int l, double k, double r, double& F, double& Fp);

int coul_F (int l, double k, double r, double & F, double & Fp);

double coul_F_asy(int l, double k, double r, double sigma = special::constant::Nan);

double coul_F_sigma(int l, double k);

bool makes_triangle (int two_j1, int two_j2, int two_j3);

double logdelta (int two_j1, int two_j2, int two_j3);

double Wigner3j_2 (int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3);
#define Wigner3j(a,b,c,d,e,f) Wigner3j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))


double Wigner6j_2 (int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6);
#define Wigner6j(a,b,c,d,e,f) Wigner6j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))


double computef (int lambda, int l1, int l2, int l1p, int l2p, int L);

double ClebschGordan(int l1, int m1, int l2, int m2, int L, int M);

double Gaunt(int l1, int m1, int l2, int m2, int l, int m);

int triangle_count(int L, int maxell);

double Hyper2F1 (double a, double b, double c, double x);

#endif