CLASSpp Manual
Cosmology reference and developer manual
Loading...
Searching...
No Matches
hyperspherical.h
1
5#ifndef __HYPERSPHERICAL__
6#define __HYPERSPHERICAL__
7
8#include <algorithm>
9#include <vector>
10
11#include "common.h"
12
13#define _HYPER_OVERFLOW_ 1e200
14#define _ONE_OVER_HYPER_OVERFLOW_ 1e-200
15#define _HYPER_SAFETY_ 1e-5
16#define _TRIG_PRECISSION_ 1e-7
17#define _HYPER_BLOCK_ 8
18#define _HYPER_CHUNK_ 16
19#define _TWO_OVER_THREE_ 0.666666666666666666666666666667e0
20#define _HIS_BYTE_ALIGNMENT_ 16
21
22typedef struct HypersphericalInterpolationStructure {
23 int K; //Sign of the curvature, (0,-1,1)
24 double beta;
25 double delta_x; //x-spacing. (xvec is uniformly spaced)
26 int trig_order; //Order of the interpolation formula for SinK and CosK.
27 int l_size; //Number of l values
28 std::vector<int> l; //Vector of l values stored
29 std::vector<double> chi_at_phimin; // vector x_min[index-l] below which neglect Bessels
30 int x_size; //Number of x-values
31 std::vector<double> x; //Pointer to x-values
32 std::vector<double> sinK; //Vector of sin_K(xvec)
33 std::vector<double> cotK; //Vector of cot_K(xvec)
34 std::vector<double> phi; //array of size nl*nx. [y_{l1}(x1) t_{l1}(x2)...]
35 std::vector<double> dphi; //Same as phivec, but containing derivatives.
36
37 HypersphericalInterpolationStructure() = default;
38 HypersphericalInterpolationStructure(int K,
39 double beta,
40 int nl,
41 const int* lvec,
42 double xmin,
43 double xmax,
44 double sampling,
45 int l_WKB,
46 double phiminabs);
47} HyperInterpStruct;
48
49struct WKB_parameters {
50 int K;
51 int l;
52 double beta;
53 double phiminabs;
54};
55
56void hyperspherical_forwards_recurrence(int K,
57 int lmax,
58 double beta,
59 double x,
60 double sinK,
61 double cotK,
62 double* sqrtK,
63 double* one_over_sqrtK,
64 double* PhiL);
65void hyperspherical_forwards_recurrence_chunk(int K,
66 int lmax,
67 double beta,
68 double* x,
69 double* sinK,
70 double* cotK,
71 int chunk,
72 double* sqrtK,
73 double* one_over_sqrtK,
74 double* PhiL);
75void hyperspherical_backwards_recurrence(int K,
76 int lmax,
77 double beta,
78 double x,
79 double sinK,
80 double cotK,
81 double* sqrtK,
82 double* one_over_sqrtK,
83 double* PhiL);
84
85void hyperspherical_backwards_recurrence_chunk(int K,
86 int lmax,
87 double beta,
88 double* x,
89 double* sinK,
90 double* cotK,
91 int chunk,
92 double* sqrtK,
93 double* one_over_sqrtK,
94 double* PhiL);
95
96bool hyperspherical_WKB(int K, int l, double beta, double y, double* Phi);
97void hyperspherical_bessel_direct_vector(
98 int K, double beta, int* lvec, int nl, double* xvec, int nx, double* Phi);
99void ClosedModY(int l, int beta, double* y, int* phisign, int* dphisign);
100bool get_CF1(int K, int l, double beta, double cotK, double* CF, int* isign);
101bool CF1_from_Gegenbauer(int l, int beta, double sinK, double cotK, double* CF);
102double airy_cheb_approx(double z);
103double coef1(double z);
104double coef2(double z);
105double coef3(double z);
106double coef4(double z);
107double cheb(double x, int n, const double A[]);
108
109double PhiWKB_minus_phiminabs(double x, void* param);
110
111void hyperspherical_get_xmin_from_Airy(
112 int K, int l, double beta, double xtol, double phiminabs, double* xmin, int* fevals);
113
114bool fzero_ridder(double (*func)(double, void*),
115 double x1,
116 double x2,
117 double xtol,
118 void* param,
119 double* Fx1,
120 double* Fx2,
121 double* xzero,
122 int* fevals);
123
124void hyperspherical_get_xmin_from_approx(
125 int K, int l, double nu, double ignore1, double phiminabs, double* xmin, int* ignore2);
126
143template <int Order, bool DoPhi, bool DoDPhi, bool DoD2Phi>
144void hyperspherical_Hermite_interpolation(const HyperInterpStruct* pHIS,
145 int nxi,
146 int lnum,
147 const double* xinterp,
148 double* Phi,
149 double* dPhi,
150 double* d2Phi);
151
152// The definition lives in hyperspherical.cpp (same translation unit as
153// ClosedModY, which must inline into the closed-case loop); these are the
154// combinations the transfer module dispatches to.
155extern template void hyperspherical_Hermite_interpolation<4, true, false, false>(
156 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
157extern template void hyperspherical_Hermite_interpolation<4, false, true, false>(
158 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
159extern template void hyperspherical_Hermite_interpolation<4, true, true, false>(
160 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
161extern template void hyperspherical_Hermite_interpolation<4, true, false, true>(
162 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
163extern template void hyperspherical_Hermite_interpolation<4, true, true, true>(
164 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
165extern template void hyperspherical_Hermite_interpolation<6, true, false, false>(
166 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
167extern template void hyperspherical_Hermite_interpolation<6, false, true, false>(
168 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
169extern template void hyperspherical_Hermite_interpolation<6, true, true, false>(
170 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
171extern template void hyperspherical_Hermite_interpolation<6, true, false, true>(
172 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
173extern template void hyperspherical_Hermite_interpolation<6, true, true, true>(
174 const HyperInterpStruct*, int, int, const double*, double*, double*, double*);
175
176#endif