greens_functions/SphericalBesselGenerator.hpp at master · ecell/greens_functions · GitHub

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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
#ifndef GREENS_FUNCTIONS_SPHERICALBESSELGENERATOR_HPP
#define GREENS_FUNCTIONS_SPHERICALBESSELGENERATOR_HPP

#include <cmath>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_bessel.h>

#include "Defs.hpp"

#ifdef NO_BESSEL_TABLE
#include "tablegen/sjy_table.hpp"

namespace sb_table
{

struct Table
{
    unsigned int N;
    double x_start;
    double delta_x;
    std::vector<double> y;
};

const unsigned int sj_table_min = 4;
const unsigned int sj_table_max = 51;
const unsigned int sy_table_min = 3;
const unsigned int sy_table_max = 40;
const unsigned int sjy_table_resolution = 35;

} // sb_table
#endif

namespace greens_functions
{

class SphericalBesselGenerator
{
public:

    SphericalBesselGenerator()
    {
#ifdef NO_BESSEL_TABLE
        sjy_table table = jnyn(
            std::max(sb_table::sj_table_max, sb_table::sy_table_max), sb_table::sjy_table_resolution);

        sj_table_.resize(sb_table::sj_table_max + 1);
        for (unsigned int n(sb_table::sj_table_min); n<= sb_table::sj_table_max; ++n)
        {
            const int start(searchsorted(table.z, minz_j(n)));
            const double z_start(table.z.at(start));
            const int end(searchsorted(table.z, maxz_j(n)));
            const std::vector<double> js(get_sub_sequence_from_matrix2(table.j, table.jdot, n, start, end));

            const sb_table::Table sj_table_n = {end - start, z_start, table.delta, js};
            sj_table_[n] = sj_table_n;
        }

        sy_table_.resize(sb_table::sy_table_max + 1);
        for (unsigned int n(sb_table::sy_table_min); n<= sb_table::sy_table_max; ++n)
        {
            const int start(searchsorted(table.z, minz_y(n)));
            const double z_start(table.z.at(start));
            const int end(searchsorted(table.z, maxz_y(n)));
            const std::vector<double> ys(get_sub_sequence_from_matrix2(table.y, table.ydot, n, start, end));

            const sb_table::Table sy_table_n = {end - start, z_start, table.delta, ys};
            sy_table_[n] = sy_table_n;
        }
#endif
    }

    ~SphericalBesselGenerator()
    {
        ; // do nothing
    }

    Real j(UnsignedInteger n, Real z) const;

    Real y(UnsignedInteger n, Real z) const;

    static UnsignedInteger getMinNJ();
    static UnsignedInteger getMinNY();
    static UnsignedInteger getMaxNJ();
    static UnsignedInteger getMaxNY();

    static SphericalBesselGenerator const& instance();

#ifdef NO_BESSEL_TABLE
    Real _j_table(UnsignedInteger n, Real z) const;
    Real _y_table(UnsignedInteger n, Real z) const;
    sb_table::Table const* getSJTable(UnsignedInteger n) const;
    sb_table::Table const* getSYTable(UnsignedInteger n) const;

private:

    std::vector<sb_table::Table> sj_table_;
    std::vector<sb_table::Table> sy_table_;
#endif
};

}


#endif /* __SPHERICALBESSELGENERATOR_HPP */