eko 0.15.3

Evolution Kernel Operators
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

//! Tools for Mellin inversion.
//!
//! We provide all necessary toold to deal with the
//! [inverse Mellin transformation](https://en.wikipedia.org/wiki/Mellin_inversion_theorem).

use num::complex::Complex;
use std::f64::consts::PI;

/// Talbot inversion path.
///
/// Implements the algorithm presented in [\[Abate\]](crate::bib::Abate).
///
/// $$p_{\text{Talbot}}(t) =  o + r \cdot ( \theta \cot(\theta) + i\theta) ~ \text{with}~ \theta = \pi(2t-1)$$
///
/// The default values for the parameters $r,o$ are given by
/// $r = \frac{2}{5} \frac{16}{0.1 - \ln(x)}$ and $o = 0$ for
/// the non-singlet integrals and $o = 1$ for the singlet sector.
/// Note that the non-singlet kernels evolve poles only up to
/// $N=0$ whereas the singlet kernels have poles up to $N=1$.
pub struct TalbotPath {
    /// integration variable
    t: f64,

    /// bending variable
    r: f64,

    /// real offset
    o: f64,
}

impl TalbotPath {
    /// Auxiliary angle.
    fn theta(&self) -> f64 {
        PI * (2.0 * self.t - 1.0)
    }

    /// Constructor from parameters.
    pub fn new(t: f64, logx: f64, is_singlet: bool) -> Self {
        // The prescription suggested by Abate for r is 0.4 * M / ( - logx)
        // with M the number of accurate digits; Maria Ubiali suggested in her thesis M = 16.
        // However, this seems to yield unstable results for the OME in the large x region
        // so we add an additional regularization, which makes the path less "edgy" there.
        Self {
            t,
            r: 0.4 * 16.0 / (0.1 - logx),
            o: if is_singlet { 1. } else { 0. },
        }
    }

    /// Mellin N.
    pub fn n(&self) -> Complex<f64> {
        let theta = self.theta();
        // treat singular point separately
        let re = if 0.5 == self.t {
            1.0
        } else {
            theta / theta.tan()
        };
        let im = theta;
        self.o + self.r * Complex::<f64> { re, im }
    }

    /// Transformation jacobian.
    pub fn jac(&self) -> Complex<f64> {
        let theta = self.theta();
        // treat singular point separately
        let re = if 0.5 == self.t {
            0.0
        } else {
            1.0 / theta.tan() - theta / theta.sin().powf(2.)
        };
        let im = 1.0;
        self.r * PI * 2.0 * Complex::<f64> { re, im }
    }

    /// Mellin inversion prefactor.
    pub fn prefactor(&self) -> Complex<f64> {
        Complex::<f64> {
            re: 0.,
            im: -1. / PI,
        }
    }
}