GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/util/Sphericals.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 318 390 81.5%
Functions: 12 13 92.3%
Branches: 298 634 47.0%
Line Branch Exec Source
1 // This file is part of Bembel, the higher order C++ boundary element library.
2 //
3 // Copyright (C) 2022 see <http://www.bembel.eu>
4 //
5 // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
6 // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
7 // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
8 // source code is subject to the GNU General Public License version 3 and
9 // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
10 // information.
11
12 #ifndef BEMBEL_SRC_UTIL_SPHERICALS_HPP_
13 #define BEMBEL_SRC_UTIL_SPHERICALS_HPP_
14 #include <Eigen/Dense>
15
16 #if !defined pi
17 #define pi BEMBEL_PI
18 #endif
19
20 namespace Bembel {
21
22 double evaluate_sphericals(Eigen::Vector3d x, Eigen::VectorXd cs,
23 unsigned int deg);
24
25 double evaluate_solid_sphericals(Eigen::Vector3d x, Eigen::VectorXd cs,
26 unsigned int deg, bool grad);
27
28 Eigen::Vector3d evaluate_dsphericals(Eigen::Vector3d x, Eigen::VectorXd cs,
29 unsigned int deg);
30
31 Eigen::Vector3d evaluate_dsolid_sphericals(Eigen::Vector3d x,
32 Eigen::VectorXd cs, unsigned int deg);
33
34 Eigen::Matrix<double, Eigen::Dynamic, 2> spherical_harmonics_full(
35 Eigen::Vector3d x, unsigned int N);
36
37 inline Eigen::Vector2d spherical_prev(Eigen::Vector3d x, int m, int n,
38 Eigen::Vector2d y1, Eigen::Vector2d y2);
39
40 Eigen::Matrix<double, 3, Eigen::Dynamic> Dsolid_harmonics_full(
41 Eigen::Vector3d x, unsigned int N, Eigen::MatrixXd spherical_val);
42
43 Eigen::Matrix<double, 3, 2> dsolid_spherical_prev(Eigen::Vector3d y, int m,
44 unsigned int n, Eigen::VectorXd L, double y_re, double y_im);
45
46 Eigen::VectorXd legendreFull(unsigned int N, double t);
47
48 double constant(int m, unsigned int n);
49
50 inline Eigen::Matrix3d functionalMatrix(Eigen::Vector3d z);
51
52 inline double pow_int(double x, int n);
53
54 /**
55 * \brief Evaluates the series \f$ \sum_{n = 0}^{\rm deg} \sum_{m = -n}^n c_m^n Y_m^n(x) \f$ for real coefficients,
56 * with the convenction that \f$ Y_{-m}^n := \overline{Y_m^n} \f$.
57 *
58 * \see K. Giebermann. _Schnelle Summationsverfahren zur numerischen Lösung von Integralgleichungen
59 * für Streuprobleme im_ \f$ \mathbb{R}^3 \f$. PhD thesis, Universität Karlsruhe (TH), Germany, 1997.
60 * \see R. von Rickenbach. _Boundary Element Methods for Shape Optimisation in Homogenisation_.
61 * MSc thesis, Universität Basel, Switzerland, 2021.
62 *
63 * @param x: The point of evaluation, a vector with length 1
64 * @param cs: The coefficients stored in the order \f$ [(0, 0), (1, -1), (1, 0), (1, 1)
65 * (2, -2), (2, -1), ... (n, n)] \f$
66 * @param deg: The degree
67 */
68 6615 double evaluate_sphericals(Eigen::Vector3d x, Eigen::VectorXd cs,
69 unsigned int deg) {
70 unsigned int m, n;
71 double z1[2], z2[2], z3[2], z1_start[2];
72 double r, fac, rootTimesZ, root_2, root_3;
73
74
2/4
✓ Branch 1 taken 6615 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6615 times.
 6615 assert(abs(x.norm() - 1) < Constants::generic_tolerance);
75
76 6615 r = z1[1] = 0;
77 6615 z1[0] = 0.5 / sqrt(pi);
78
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6615 times.
 6615 if (deg <= 1) {
79 ✗ return cs(0) * z1[0];
80 }
81
82
2/2
✓ Branch 0 taken 72765 times.
✓ Branch 1 taken 6615 times.
 79380 for (m = 0; m < deg - 1; m++) {
83
2/2
✓ Branch 0 taken 6615 times.
✓ Branch 1 taken 66150 times.
 72765 if (m == 0) {
84 6615 fac = 1.0;
85 } else {
86 66150 fac = 2.0;
87 }
88
89 72765 z1_start[0] = z1[0];
90 72765 z1_start[1] = z1[1];
91
1/2
✓ Branch 1 taken 72765 times.
✗ Branch 2 not taken.
 72765 rootTimesZ = sqrt(2 * m + 3) * x(2);
92 72765 z2[0] = rootTimesZ * z1[0];
93 72765 z2[1] = rootTimesZ * z1[1];
94
1/2
✓ Branch 1 taken 72765 times.
✗ Branch 2 not taken.
 72765 r += fac * cs(m * (m + 1) + m) * z1[0]; // + k[ m *(m+1)-m]*z1[1];
95
1/2
✓ Branch 1 taken 72765 times.
✗ Branch 2 not taken.
 72765 r += fac * cs((m + 1) * (m + 2) + m) * z2[0]; // + k[(m+1)*(m+2)-m]*z2[1];
96
2/2
✓ Branch 0 taken 363825 times.
✓ Branch 1 taken 72765 times.
 436590 for (n = m + 2; n < deg; n++) {
97 363825 root_2 = sqrt((2 * n + 1.0) / ((n - m) * (n + m)));
98
1/2
✓ Branch 1 taken 363825 times.
✗ Branch 2 not taken.
 363825 rootTimesZ = sqrt(2 * n - 1) * x(2);
99 363825 root_3 = sqrt(((n + m - 1.0) * (n - m - 1.0)) / (2 * n - 3));
100 363825 z3[0] = root_2 * (rootTimesZ * z2[0] - root_3 * z1[0]);
101 363825 z3[1] = root_2 * (rootTimesZ * z2[1] - root_3 * z1[1]);
102
1/2
✓ Branch 1 taken 363825 times.
✗ Branch 2 not taken.
 363825 r += fac * cs(n * (n + 1) + m) * z3[0]; // + k[n*(n+1)-m]*z3[1];
103 363825 z1[0] = z2[0];
104 363825 z1[1] = z2[1];
105 363825 z2[0] = z3[0];
106 363825 z2[1] = z3[1];
107 }
108 72765 root_2 = sqrt((2 * m + 3.0) / (2 * m + 2));
109
2/4
✓ Branch 1 taken 72765 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72765 times.
✗ Branch 5 not taken.
 72765 z1[0] = root_2 * (x(0) * z1_start[0] - x(1) * z1_start[1]);
110
2/4
✓ Branch 1 taken 72765 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72765 times.
✗ Branch 5 not taken.
 72765 z1[1] = root_2 * (x(0) * z1_start[1] + x(1) * z1_start[0]);
111 }
112
1/2
✓ Branch 1 taken 6615 times.
✗ Branch 2 not taken.
 6615 r += 2 * cs((deg - 1) * (deg + 1)) * z1[0]; // + k[(nk-1)*(nk-1)]*z1[1];
113 6615 return r;
114 }
115
116 /**
117 * \brief Evaluates the series \f$ \sum_{n = 0}^{\rm deg} \sum_{m = -n}^n
118 * |x|^n c_m^n Y_m^n( \frac{x}{|x|}) \f$ for real coefficients cs if grad is false, and
119 * \f$ \sum_{n = 0}^{\rm deg} \sum_{m = -n}^n n |x|^{n-1} c_m^n Y_m^n(\frac{x}{|x|}) \f$
120 * for real coefficients cs if grad is true
121 *
122 * @param x: The point of evaluation
123 * @param cs: The coefficients stored in the order \f$ [(0, 0), (1, -1), (1, 0), (1, 1)
124 * (2, -2), (2, -1), ... (n, n)] \f$
125 * @param deg: The degree
126 * @param grad: distinguish the two cases.
127 */
128 3279303 double evaluate_solid_sphericals(Eigen::Vector3d x, Eigen::VectorXd cs,
129 unsigned int deg, bool grad) {
130 unsigned int m, n;
131 double z1[2], z2[2], z3[2], z1_start[2];
132 double r, fac, rootTimesZ, root_2, root_3, norm, fac_tot;
133
1/2
✓ Branch 1 taken 3279303 times.
✗ Branch 2 not taken.
 3279303 Eigen::Vector3d y;
134
135 3279303 r = z1[1] = 0;
136 3279303 z1[0] = 0.5 / sqrt(pi);
137
3/4
✓ Branch 0 taken 61207 times.
✓ Branch 1 taken 3218096 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 61207 times.
 3279303 if (grad && deg <= 1) {
138 ✗ return 0.0;
139
3/4
✓ Branch 0 taken 3218096 times.
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3218096 times.
 3279303 } else if (!grad && deg <= 1) {
140 ✗ return cs(0) * z1[0];
141 }
142
143
1/2
✓ Branch 1 taken 3279303 times.
✗ Branch 2 not taken.
 3279303 norm = x.norm();
144
2/4
✓ Branch 1 taken 3279303 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3279303 times.
✗ Branch 5 not taken.
 3279303 y = x / norm;
145
146
2/2
✓ Branch 0 taken 36072333 times.
✓ Branch 1 taken 3279303 times.
 39351636 for (m = 0; m < deg - 1; m++) {
147
2/2
✓ Branch 0 taken 3279303 times.
✓ Branch 1 taken 32793030 times.
 36072333 if (m == 0) {
148 3279303 fac = 1.0;
149 } else {
150 32793030 fac = 2.0;
151 }
152
153
2/2
✓ Branch 0 taken 673277 times.
✓ Branch 1 taken 35399056 times.
 36072333 if (grad) {
154 673277 fac_tot = fac * m * pow_int(norm, m - 1);
155 } else {
156 35399056 fac_tot = fac * pow_int(norm, m);
157 }
158
159 36072333 z1_start[0] = z1[0];
160 36072333 z1_start[1] = z1[1];
161
1/2
✓ Branch 1 taken 36072333 times.
✗ Branch 2 not taken.
 36072333 rootTimesZ = sqrt(2 * m + 3) * y(2);
162 36072333 z2[0] = rootTimesZ * z1[0];
163 36072333 z2[1] = rootTimesZ * z1[1];
164
1/2
✓ Branch 1 taken 36072333 times.
✗ Branch 2 not taken.
 36072333 r += fac_tot * cs(m * (m + 1) + m) * z1[0];
165
166
2/2
✓ Branch 0 taken 673277 times.
✓ Branch 1 taken 35399056 times.
 36072333 if (grad) {
167 673277 fac_tot = fac * (m + 1) * pow_int(norm, m);
168 } else {
169 35399056 fac_tot *= norm;
170 }
171
1/2
✓ Branch 1 taken 36072333 times.
✗ Branch 2 not taken.
 36072333 r += fac_tot * cs((m + 1) * (m + 2) + m) * z2[0];
172
173
2/2
✓ Branch 0 taken 180361665 times.
✓ Branch 1 taken 36072333 times.
 216433998 for (n = m + 2; n < deg; n++) {
174
2/2
✓ Branch 0 taken 3366385 times.
✓ Branch 1 taken 176995280 times.
 180361665 if (grad) {
175 3366385 fac_tot = fac * n * pow_int(norm, n - 1);
176 } else {
177 176995280 fac_tot *= norm;
178 }
179 180361665 root_2 = sqrt((2 * n + 1.0) / ((n - m) * (n + m)));
180
1/2
✓ Branch 1 taken 180361665 times.
✗ Branch 2 not taken.
 180361665 rootTimesZ = sqrt(2 * n - 1) * y(2);
181 180361665 root_3 = sqrt(((n + m - 1.0) * (n - m - 1.0)) / (2 * n - 3));
182 180361665 z3[0] = root_2 * (rootTimesZ * z2[0] - root_3 * z1[0]);
183 180361665 z3[1] = root_2 * (rootTimesZ * z2[1] - root_3 * z1[1]);
184
1/2
✓ Branch 1 taken 180361665 times.
✗ Branch 2 not taken.
 180361665 r += fac_tot * cs(n * (n + 1) + m) * z3[0];
185 180361665 z1[0] = z2[0];
186 180361665 z1[1] = z2[1];
187 180361665 z2[0] = z3[0];
188 180361665 z2[1] = z3[1];
189 }
190 36072333 root_2 = sqrt((2 * m + 3.0) / (2 * m + 2));
191
2/4
✓ Branch 1 taken 36072333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36072333 times.
✗ Branch 5 not taken.
 36072333 z1[0] = root_2 * (y(0) * z1_start[0] - y(1) * z1_start[1]);
192
2/4
✓ Branch 1 taken 36072333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36072333 times.
✗ Branch 5 not taken.
 36072333 z1[1] = root_2 * (y(0) * z1_start[1] + y(1) * z1_start[0]);
193 }
194
195
2/2
✓ Branch 0 taken 61207 times.
✓ Branch 1 taken 3218096 times.
 3279303 if (grad) {
196 61207 fac_tot = fac * deg * pow_int(norm, deg - 1);
197 } else {
198 3218096 fac_tot *= norm;
199 }
200
1/2
✓ Branch 1 taken 3279303 times.
✗ Branch 2 not taken.
 3279303 r += fac_tot * cs((deg - 1) * (deg + 1)) * z1[0];
201 3279303 return r;
202 }
203
204 /**
205 * \brief Evaluates the series \f$ \sum_{n = 0}^{\rm deg} \sum_{m = -n}^n c_m^n
206 * \nabla Y_m^n(x) \f$ for real coefficients.
207 *
208 * @param x: The point of evaluation, a vector with length 1
209 * @param cs: The coefficients stored in the order \f$ [(0, 0), (1, -1), (1, 0), (1, 1)
210 * (2, -2), (2, -1), ... (n, n)] \f$
211 * @param deg: The degree
212 */
213 ✗ Eigen::Vector3d evaluate_dsphericals(Eigen::Vector3d x, Eigen::VectorXd cs,
214 unsigned int deg) {
215 unsigned int m, n;
216 double z1[2], z2[2], z3[2], z1_start[2];
217 ✗ Eigen::Vector3d dr;
218 double fac, rootTimesZ, root_2, root_3;
219
220 ✗ assert(abs(x.norm() - 1) < Constants::generic_tolerance);
221
222 ✗ z1[0] = sqrt(0.375 / pi);
223 ✗ z1[1] = 0;
224
225 ✗ dr = Eigen::Vector3d(0.0, 0.0, 0.0);
226
227 ✗ if (deg == 1) {
228 ✗ return dr;
229 }
230
231 ✗ for (m = 1; m < deg - 1; m++) {
232 ✗ z1_start[0] = z1[0];
233 ✗ z1_start[1] = z1[1];
234 ✗ rootTimesZ = sqrt(2 * m + 3) * x(2);
235 ✗ z2[0] = rootTimesZ * z1[0];
236 ✗ z2[1] = rootTimesZ * z1[1];
237 ✗ dr(0) += 2 * m * (cs(m * (m + 1) + m) * z1[0]);
238 // + a[ m *(m+1)-m]*z1[1]) for imaginary coefficients;
239 ✗ dr(0) += 2 * m * (cs((m + 1) * (m + 2) + m) * z2[0]);
240 // + a[(m+1)*(m+2)-m]*z2[1]) for imaginary coefficients;
241 ✗ dr(1) -= 2 * m * (cs(m * (m + 1) + m) * z1[1]);
242 // - a[ m *(m+1)-m]*z1[0]) for imaginary coefficients;
243 ✗ dr(1) -= 2 * m * (cs((m + 1) * (m + 2) + m) * z2[1]);
244 // - a[(m+1)*(m+2)-m]*z2[0]) for imaginary coefficients;
245
246 ✗ if (m == 1) {
247 ✗ fac = 1.0;
248 ✗ dr(2) += fac * sqrt(2 * m)
249 ✗ * (cs(m * (m + 1) + (m - 1)) * z1[0]
250 ✗ + cs(m * (m + 1) - (m - 1)) * z1[1]);
251 ✗ dr(2) += fac * sqrt(4 * m + 2)
252 ✗ * (cs((m + 1) * (m + 2) + (m - 1)) * z2[0]
253 ✗ + cs((m + 1) * (m + 2) - (m - 1)) * z2[1]);
254 } else {
255 ✗ fac = 2.0;
256 ✗ dr(2) += fac * sqrt(2 * m) * (cs(m * (m + 1) + (m - 1)) * z1[0]);
257 // + a[ m *(m+1)-(m-1)]*z1[1]) for imaginary coefficients;
258 ✗ dr(2) += fac * sqrt(4 * m + 2)
259 ✗ * (cs((m + 1) * (m + 2) + (m - 1)) * z2[0]);
260 // + a[(m+1)*(m+2)-(m-1)]*z2[1]) for imaginary coefficients;
261 }
262
263 ✗ for (n = m + 2; n < deg; n++) {
264 ✗ root_2 = sqrt((2 * n + 1.0) / ((n - m) * (n + m)));
265 ✗ rootTimesZ = sqrt(2 * n - 1) * x(2);
266 ✗ root_3 = sqrt(((n + m - 1.0) * (n - m - 1.0)) / (2 * n - 3));
267 ✗ z3[0] = root_2 * (rootTimesZ * z2[0] - root_3 * z1[0]);
268 ✗ z3[1] = root_2 * (rootTimesZ * z2[1] - root_3 * z1[1]);
269 ✗ dr(0) += 2 * m * (cs(n * (n + 1) + m) * z3[0]); // + a[n*(n+1)-m]*z3[1]);
270 ✗ dr(1) -= 2 * m * (cs(n * (n + 1) + m) * z3[1]); // - a[n*(n+1)-m]*z3[0]);
271 ✗ dr(2) += fac * sqrt((n + m) * (n - m + 1))
272 ✗ * (cs(n * (n + 1) + (m - 1)) * z3[0]); // + a[n*(n+1)-(m-1)]*z3[1]);
273 ✗ z1[0] = z2[0];
274 ✗ z1[1] = z2[1];
275 ✗ z2[0] = z3[0];
276 ✗ z2[1] = z3[1];
277 }
278 ✗ root_2 = sqrt((2 * m + 3.0) / (2 * m + 2));
279 ✗ z1[0] = root_2 * (x(0) * z1_start[0] - x(1) * z1_start[1]);
280 ✗ z1[1] = root_2 * (x(0) * z1_start[1] + x(1) * z1_start[0]);
281 }
282 ✗ dr(0) += (deg - 1) * (2 * cs((deg - 1) * (deg + 1)) * z1[0]);
283 // + a[(na-1)*(na-1)]*z1[1]) for imaginary coefficients;
284 ✗ dr(1) -= (deg - 1) * (2 * cs((deg - 1) * (deg + 1)) * z1[1]);
285 // - a[(na-1)*(na-1)]*z1[0]) for imaginary coefficients;
286 ✗ dr(2) += 2 * sqrt(2 * (deg - 1))
287 ✗ * (cs(deg * (deg - 1) + (deg - 2)) * z1[0]
288 ✗ + cs(deg * (deg - 1) - (deg - 2)) * z1[1]);
289
290 ✗ return dr;
291 }
292
293 /**
294 * \brief Evaluates the series \f$ \sum_{n = 0}^{\rm deg} \sum_{m = -n}^n
295 * |x|^{n-3} c_m^n (\nabla Y_m^n)(\frac{x}{|x|}) \f$ for real coefficients.
296 *
297 * @param x: The point of evaluation
298 * @param cs: The coefficients stored in the order \f$ [(0, 0), (1, -1), (1, 0), (1, 1)
299 * (2, -2), (2, -1), ... (n, n)] \f$
300 * @param deg: The degree
301 */
302 61207 Eigen::Vector3d evaluate_dsolid_sphericals(Eigen::Vector3d x,
303 Eigen::VectorXd cs, unsigned int deg) {
304 unsigned int m, n;
305 double z1[2], z2[2], z3[2], z1_start[2];
306
2/4
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61207 times.
✗ Branch 5 not taken.
 61207 Eigen::Vector3d dr, y;
307 double fac, rootTimesZ, root_2, root_3, norm, r_n3;
308
309 61207 z1[0] = sqrt(0.375 / pi);
310 61207 z1[1] = 0;
311
312
1/2
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
 61207 dr = Eigen::Vector3d(0.0, 0.0, 0.0);
313
1/2
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
 61207 norm = x.norm();
314
2/4
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61207 times.
✗ Branch 5 not taken.
 61207 y = x / norm;
315
316
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 61207 times.
 61207 if (deg == 1) {
317 ✗ return dr;
318 }
319
320
2/2
✓ Branch 0 taken 612070 times.
✓ Branch 1 taken 61207 times.
 673277 for (m = 1; m < deg - 1; m++) {
321 612070 r_n3 = pow_int(norm, m - 3);
322
323 612070 z1_start[0] = z1[0];
324 612070 z1_start[1] = z1[1];
325
1/2
✓ Branch 1 taken 612070 times.
✗ Branch 2 not taken.
 612070 rootTimesZ = sqrt(2 * m + 3) * y(2);
326 612070 z2[0] = rootTimesZ * z1[0];
327 612070 z2[1] = rootTimesZ * z1[1];
328
2/4
✓ Branch 1 taken 612070 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 612070 times.
✗ Branch 5 not taken.
 612070 dr(0) += 2 * m * r_n3 * (cs(m * (m + 1) + m) * z1[0]);
329
2/4
✓ Branch 1 taken 612070 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 612070 times.
✗ Branch 5 not taken.
 612070 dr(0) += 2 * m * r_n3 * norm * (cs((m + 1) * (m + 2) + m) * z2[0]);
330
2/4
✓ Branch 1 taken 612070 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 612070 times.
✗ Branch 5 not taken.
 612070 dr(1) -= 2 * m * r_n3 * (cs(m * (m + 1) + m) * z1[1]);
331
2/4
✓ Branch 1 taken 612070 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 612070 times.
✗ Branch 5 not taken.
 612070 dr(1) -= 2 * m * r_n3 * norm * (cs((m + 1) * (m + 2) + m) * z2[1]);
332
333
2/2
✓ Branch 0 taken 61207 times.
✓ Branch 1 taken 550863 times.
 612070 if (m == 1) {
334 61207 fac = 1.0;
335 183621 dr(2) += fac * r_n3 * sqrt(2 * m)
336
1/2
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
 61207 * (cs(m * (m + 1) + (m - 1)) * z1[0]
337
2/4
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61207 times.
✗ Branch 5 not taken.
 61207 + cs(m * (m + 1) - (m - 1)) * z1[1]);
338 61207 dr(2) += fac * r_n3 * norm * sqrt(4 * m + 2)
339
1/2
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
 61207 * (cs((m + 1) * (m + 2) + (m - 1)) * z2[0]
340
2/4
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61207 times.
✗ Branch 5 not taken.
 61207 + cs((m + 1) * (m + 2) - (m - 1)) * z2[1]);
341 } else {
342 550863 fac = 2.0;
343
2/4
✓ Branch 1 taken 550863 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 550863 times.
✗ Branch 5 not taken.
 550863 dr(2) += fac * r_n3 * sqrt(2 * m) * (cs(m * (m + 1) + (m - 1)) * z1[0]);
344 550863 dr(2) += fac * r_n3 * norm * sqrt(4 * m + 2)
345
2/4
✓ Branch 1 taken 550863 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 550863 times.
✗ Branch 5 not taken.
 550863 * (cs((m + 1) * (m + 2) + (m - 1)) * z2[0]);
346 }
347
348 612070 r_n3 *= norm;
349
350
2/2
✓ Branch 0 taken 2754315 times.
✓ Branch 1 taken 612070 times.
 3366385 for (n = m + 2; n < deg; n++) {
351 2754315 r_n3 *= norm;
352 2754315 root_2 = sqrt((2 * n + 1.0) / ((n - m) * (n + m)));
353
1/2
✓ Branch 1 taken 2754315 times.
✗ Branch 2 not taken.
 2754315 rootTimesZ = sqrt(2 * n - 1) * y(2);
354 2754315 root_3 = sqrt(((n + m - 1.0) * (n - m - 1.0)) / (2 * n - 3));
355 2754315 z3[0] = root_2 * (rootTimesZ * z2[0] - root_3 * z1[0]);
356 2754315 z3[1] = root_2 * (rootTimesZ * z2[1] - root_3 * z1[1]);
357
2/4
✓ Branch 1 taken 2754315 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2754315 times.
✗ Branch 5 not taken.
 2754315 dr(0) += 2 * m * r_n3 * (cs(n * (n + 1) + m) * z3[0]);
358
2/4
✓ Branch 1 taken 2754315 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2754315 times.
✗ Branch 5 not taken.
 2754315 dr(1) -= 2 * m * r_n3 * (cs(n * (n + 1) + m) * z3[1]);
359 8262945 dr(2) += fac * r_n3 * sqrt((n + m) * (n - m + 1))
360
2/4
✓ Branch 1 taken 2754315 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2754315 times.
✗ Branch 5 not taken.
 2754315 * (cs(n * (n + 1) + (m - 1)) * z3[0]);
361 2754315 z1[0] = z2[0];
362 2754315 z1[1] = z2[1];
363 2754315 z2[0] = z3[0];
364 2754315 z2[1] = z3[1];
365 }
366 612070 root_2 = sqrt((2 * m + 3.0) / (2 * m + 2));
367
2/4
✓ Branch 1 taken 612070 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 612070 times.
✗ Branch 5 not taken.
 612070 z1[0] = root_2 * (y(0) * z1_start[0] - y(1) * z1_start[1]);
368
2/4
✓ Branch 1 taken 612070 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 612070 times.
✗ Branch 5 not taken.
 612070 z1[1] = root_2 * (y(0) * z1_start[1] + y(1) * z1_start[0]);
369 }
370 61207 r_n3 *= norm;
371
2/4
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61207 times.
✗ Branch 5 not taken.
 61207 dr(0) += (deg - 1) * r_n3 * (2 * cs((deg - 1) * (deg + 1)) * z1[0]);
372
2/4
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61207 times.
✗ Branch 5 not taken.
 61207 dr(1) -= (deg - 1) * r_n3 * (2 * cs((deg - 1) * (deg + 1)) * z1[1]);
373 183621 dr(2) += 2 * r_n3 * sqrt(2 * (deg - 1))
374
1/2
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
 61207 * (cs(deg * (deg - 1) + (deg - 2)) * z1[0]
375
2/4
✓ Branch 1 taken 61207 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61207 times.
✗ Branch 5 not taken.
 61207 + cs(deg * (deg - 1) - (deg - 2)) * z1[1]);
376
377 61207 return dr;
378 }
379
380 /**
381 * \brief Calculates the the spherical harmonics \f$ Y_n^m(\frac{x}{|x|}) \f$,
382 * ordered by \f$ [Y_0^0, \, Y_1^{-1}, \, Y_1^0, \, Y_1^1, \,
383 * Y_2^{-2}, ..., Y_N^N] \f$
384 *
385 * @param x: the point of evaluation
386 * @param N: the maximal degree
387 */
388 13230 Eigen::Matrix<double, Eigen::Dynamic, 2> spherical_harmonics_full(
389 Eigen::Vector3d x, unsigned int N) {
390
391
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 Eigen::VectorXd real((N + 1) * (N + 1));
392
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 Eigen::VectorXd imag((N + 1) * (N + 1));
393
394
3/6
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13230 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13230 times.
✗ Branch 8 not taken.
 13230 Eigen::Vector3d y = x / x.norm();
395
396 int m, n;
397
398 /* handle n = 0 separately */
399
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 real(0) = 0.5 / sqrt(pi);
400
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 imag(0) = 0.0;
401
402 /* assemble two temporary vectors */
403
3/6
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13230 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13230 times.
✗ Branch 8 not taken.
 13230 Eigen::Vector2d z1, z2, tmp;
404
405
2/2
✓ Branch 0 taken 158760 times.
✓ Branch 1 taken 13230 times.
 171990 for (n = 1; n <= N; n++) {
406
2/2
✓ Branch 0 taken 873180 times.
✓ Branch 1 taken 158760 times.
 1031940 for (m = 0; m < n - 1; m++) {
407
2/4
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
 873180 z1(0) = real((n - 1) * (n - 1) + n - 1 + m);
408
2/4
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
 873180 z1(1) = real((n - 2) * (n - 2) + n - 2 + m);
409
2/4
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
 873180 z2(0) = imag((n - 1) * (n - 1) + n - 1 + m);
410
2/4
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
 873180 z2(1) = imag((n - 2) * (n - 2) + n - 2 + m);
411
412
4/8
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 873180 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 873180 times.
✗ Branch 11 not taken.
 873180 tmp = spherical_prev(y, m, n, z1, z2);
413
414
2/4
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
 873180 real(n * n + n + m) = tmp(0);
415
2/4
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
 873180 imag(n * n + n + m) = tmp(1);
416
417 // if m > 0, copy the value with the prefactor
418
2/2
✓ Branch 0 taken 727650 times.
✓ Branch 1 taken 145530 times.
 873180 if (m > 0) {
419
2/4
✓ Branch 1 taken 727650 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 727650 times.
✗ Branch 5 not taken.
 727650 real(n * n + n - m) = tmp(0);
420
2/4
✓ Branch 1 taken 727650 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 727650 times.
✗ Branch 5 not taken.
 727650 imag(n * n + n - m) = -tmp(1);
421 }
422 }
423
424 /* for m = n-1 */
425 158760 m = n - 1;
426
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 z1(0) = real(m * m + m + m);
427
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 z2(0) = imag(m * m + m + m);
428
429
4/8
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 158760 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 158760 times.
✗ Branch 11 not taken.
 158760 tmp = spherical_prev(y, m, n, z1, z2);
430
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 real(n * n + n + m) = tmp(0);
431
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 imag(n * n + n + m) = tmp(1);
432
433
2/2
✓ Branch 0 taken 145530 times.
✓ Branch 1 taken 13230 times.
 158760 if (m > 0) {
434
2/4
✓ Branch 1 taken 145530 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 145530 times.
✗ Branch 5 not taken.
 145530 real(n * n + n - m) = tmp(0);
435
2/4
✓ Branch 1 taken 145530 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 145530 times.
✗ Branch 5 not taken.
 145530 imag(n * n + n - m) = -tmp(1);
436 }
437
438 /* for m = n */
439 158760 m = n;
440
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 z1(0) = real((n - 1) * (n - 1) + n - 1 + n - 1);
441
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 z2(0) = imag((n - 1) * (n - 1) + n - 1 + n - 1);
442
443
4/8
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 158760 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 158760 times.
✗ Branch 11 not taken.
 158760 tmp = spherical_prev(y, m, n, z1, z2);
444
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 real(n * n + n + m) = tmp(0);
445
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 imag(n * n + n + m) = tmp(1);
446
447
1/2
✓ Branch 0 taken 158760 times.
✗ Branch 1 not taken.
 158760 if (m > 0) {
448
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 real(n * n + n - m) = tmp(0);
449
2/4
✓ Branch 1 taken 158760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158760 times.
✗ Branch 5 not taken.
 158760 imag(n * n + n - m) = -tmp(1);
450 }
451 }
452
453
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 Eigen::MatrixXd res((N + 1) * (N + 1), 2);
454
2/4
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13230 times.
✗ Branch 5 not taken.
 13230 res.col(0) = real;
455
2/4
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13230 times.
✗ Branch 5 not taken.
 13230 res.col(1) = imag;
456
457
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 26460 return res;
458 13230 }
459
460 /**
461 * \brief Calculates the spherical harmonic \f$ Y_n^m(x) \f$ based on previous values
462 */
463 1190700 inline Eigen::Vector2d spherical_prev(Eigen::Vector3d x, int m, int n,
464 Eigen::Vector2d y1, Eigen::Vector2d y2) {
465 1190700 Eigen::Vector2d z;
466
467
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1190700 times.
 1190700 assert(abs(x.norm() - 1) < 1e-14);
468
469
3/4
✓ Branch 0 taken 158760 times.
✓ Branch 1 taken 1031940 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 158760 times.
 1190700 if ((m == 0) && (n == 0)) {
470 ✗ z(0) = 0.5 / sqrt(pi);
471 ✗ z(1) = 0;
472
2/2
✓ Branch 0 taken 158760 times.
✓ Branch 1 taken 1031940 times.
 1190700 } else if (m == n) {
473 158760 z(0) = sqrt((2 * m + 1.0) / (2 * m)) * (x(0) * y1(0) - x(1) * y2(0));
474 158760 z(1) = sqrt((2 * m + 1.0) / (2 * m)) * (x(0) * y2(0) + x(1) * y1(0));
475
2/2
✓ Branch 0 taken 158760 times.
✓ Branch 1 taken 873180 times.
 1031940 } else if (m + 1 == n) {
476 158760 z(0) = y1(0) * sqrt(2 * m + 3) * x(2);
477 158760 z(1) = y2(0) * sqrt(2 * m + 3) * x(2);
478 } else {
479 2619540 z(0) = sqrt((2 * n + 1.0) / ((n - m) * (n + m)))
480 873180 * (sqrt(2 * n - 1) * x(2) * y1(0)
481 873180 - sqrt(((n + m - 1.0) * (n - m - 1.0)) / (2 * n - 3)) * y1(1));
482 2619540 z(1) = sqrt((2 * n + 1.0) / ((n - m) * (n + m)))
483 1746360 * (sqrt(2 * n - 1) * x(2) * y2(0)
484 873180 - sqrt(((n + m - 1.0) * (n - m - 1.0)) / (2 * n - 3)) * y2(1));
485 }
486 1190700 return z;
487 }
488
489 /**
490 * \brief Calculates all gradients of the solid harmonics, given the Legendre Coefficients L.
491 */
492 13230 Eigen::Matrix<double, 3, Eigen::Dynamic> Dsolid_harmonics_full(
493 Eigen::Vector3d x, unsigned int N, Eigen::MatrixXd spherical_val) {
494
495
3/6
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13230 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13230 times.
✗ Branch 8 not taken.
 13230 Eigen::VectorXd L = legendreFull(N, x(2) / x.norm());
496
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 Eigen::Matrix<double, 3, 2> z;
497
498
2/4
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13230 times.
✗ Branch 5 not taken.
 13230 Eigen::MatrixXd reals(3, (N + 1) * (N + 1)), imags(3, (N + 1) * (N + 1));
499
500 int m, n;
501
502
2/2
✓ Branch 0 taken 171990 times.
✓ Branch 1 taken 13230 times.
 185220 for (n = 0; n <= N; n++) {
503
2/2
✓ Branch 0 taken 1203930 times.
✓ Branch 1 taken 171990 times.
 1375920 for (m = 0; m <= n; m++) {
504
3/6
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1203930 times.
✗ Branch 8 not taken.
 2407860 z = dsolid_spherical_prev(x, m, n, L, spherical_val(n * n + n + m, 0),
505
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 2407860 spherical_val(n * n + n + m, 1));
506
507
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 reals(0, n * n + n + m) = z(0, 0);
508
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 reals(1, n * n + n + m) = z(1, 0);
509
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 reals(2, n * n + n + m) = z(2, 0);
510
511
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 imags(0, n * n + n + m) = z(0, 1);
512
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 imags(1, n * n + n + m) = z(1, 1);
513
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 imags(2, n * n + n + m) = z(2, 1);
514
515
2/2
✓ Branch 0 taken 1031940 times.
✓ Branch 1 taken 171990 times.
 1203930 if (m > 0) {
516
2/4
✓ Branch 1 taken 1031940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1031940 times.
✗ Branch 5 not taken.
 1031940 reals(0, n * n + n - m) = z(0, 0);
517
2/4
✓ Branch 1 taken 1031940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1031940 times.
✗ Branch 5 not taken.
 1031940 reals(1, n * n + n - m) = z(1, 0);
518
2/4
✓ Branch 1 taken 1031940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1031940 times.
✗ Branch 5 not taken.
 1031940 reals(2, n * n + n - m) = z(2, 0);
519
520
2/4
✓ Branch 1 taken 1031940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1031940 times.
✗ Branch 5 not taken.
 1031940 imags(0, n * n + n - m) = -z(0, 1);
521
2/4
✓ Branch 1 taken 1031940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1031940 times.
✗ Branch 5 not taken.
 1031940 imags(1, n * n + n - m) = -z(1, 1);
522
2/4
✓ Branch 1 taken 1031940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1031940 times.
✗ Branch 5 not taken.
 1031940 imags(2, n * n + n - m) = -z(2, 1);
523 }
524 }
525 }
526
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 26460 return reals;
527 13230 }
528
529 /**
530 * \brief Calculates \f$ \nabla Y_m^n(x) \f$ based on the Legendre Coefficients L.
531 */
532 1203930 inline Eigen::Matrix<double, 3, 2> dspherical_prev(Eigen::Vector3d x, int m,
533 unsigned int n, Eigen::VectorXd L) {
534
535 double c;
536 unsigned int i;
537
538
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1203930 times.
 1203930 assert(abs(x.norm() - 1) < 1e-14);
539
540 1203930 Eigen::Matrix<double, 3, 2> z;
541
542
2/2
✓ Branch 0 taken 171990 times.
✓ Branch 1 taken 1031940 times.
 1203930 if (m == 0) {
543 171990 z(0, 0) = z(0, 1) = z(2, 1) = 0;
544
2/2
✓ Branch 0 taken 13230 times.
✓ Branch 1 taken 158760 times.
 171990 if (1 > n) {
545 13230 z(2, 0) = 0;
546 } else {
547 158760 z(2, 0) = L((n * (n + 1)) / 2 + 1);
548 }
549 } else {
550
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1031940 times.
 1031940 if (m > n) {
551 ✗ z(0, 0) = 0;
552 } else {
553 1031940 z(0, 0) = L((n * (n + 1)) / 2 + m);
554 }
555
556
2/2
✓ Branch 0 taken 158760 times.
✓ Branch 1 taken 873180 times.
 1031940 if (m + 1 > n) {
557 158760 z(2, 0) = 0;
558 } else {
559 873180 z(2, 0) = L((n * (n + 1)) / 2 + m + 1);
560 }
561 1031940 z(0, 1) = z(2, 1) = 0;
562
563 1031940 c = x(0) * z(2, 0) - x(1) * z(2, 1);
564 1031940 z(2, 1) = x(0) * z(2, 1) + x(1) * z(2, 0);
565 1031940 z(2, 0) = c;
566
567
2/2
✓ Branch 0 taken 3783780 times.
✓ Branch 1 taken 1031940 times.
 4815720 for (i = 1; i < m; i++) {
568 3783780 c = x(0) * z(0, 0) - x(1) * z(0, 1);
569 3783780 z(0, 1) = x(0) * z(0, 1) + x(1) * z(0, 0);
570 3783780 z(0, 0) = c;
571
572 3783780 c = x(0) * z(2, 0) - x(1) * z(2, 1);
573 3783780 z(2, 1) = x(0) * z(2, 1) + x(1) * z(2, 0);
574 3783780 z(2, 0) = c;
575 }
576 }
577
578 1203930 c = constant(m, n);
579 1203930 z(0, 0) *= m * c;
580 1203930 z(0, 1) *= m * c;
581 1203930 z(1, 0) = -z(0, 1);
582 1203930 z(1, 1) = +z(0, 0);
583 1203930 z(2, 0) *= c;
584 1203930 z(2, 1) *= c;
585
586 1203930 return z;
587 }
588
589 /**
590 * \brief Calculates the derivative of the solid harmonics function \f$ \phi_n^m \f$ at
591 * the point \f$ y \f$ with the spherical values.
592 */
593 1203930 Eigen::Matrix<double, 3, 2> dsolid_spherical_prev(Eigen::Vector3d y, int m,
594 unsigned int n, Eigen::VectorXd L, double y_re, double y_im) {
595
596 // normalise y
597 // x denotes the normalised y and is passed to the original functions
598
1/2
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
 1203930 double r = y.norm();
599
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 Eigen::Vector3d x = y / r;
600
1/2
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
 1203930 Eigen::Matrix<double, 3, 2> z_harm;
601
1/2
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
 1203930 Eigen::Matrix3d A;
602 unsigned int m_tilde;
603
604 // part with the gradient of Y
605
1/2
✓ Branch 0 taken 1203930 times.
✗ Branch 1 not taken.
 1203930 if (m >= 0) {
606 1203930 m_tilde = m;
607
3/6
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1203930 times.
✗ Branch 8 not taken.
 1203930 z_harm = dspherical_prev(x, m_tilde, n, L);
608
609 // multiply the gradient of Y by r^(n-3) and the matrix A
610 double r_n3;
611
2/2
✓ Branch 0 taken 1071630 times.
✓ Branch 1 taken 132300 times.
 1203930 if (n > 3) {
612 1071630 r_n3 = pow(r, n - 3);
613
2/2
✓ Branch 0 taken 52920 times.
✓ Branch 1 taken 79380 times.
 132300 } else if (n == 3) {
614 52920 r_n3 = 1.0;
615 } else {
616 79380 r_n3 = pow(1.0 / r, 3 - n);
617 }
618
619
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 1203930 A = functionalMatrix(y);
620
3/6
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1203930 times.
✗ Branch 8 not taken.
 1203930 z_harm = r_n3 * (A * z_harm);
621
622 } else {
623 ✗ m_tilde = -m;
624 ✗ z_harm = dspherical_prev(x, m_tilde, n, L);
625
626 // define r^(n-3) as above
627 double r_n3;
628 ✗ if (n > 3) {
629 ✗ r_n3 = pow(r, n - 3);
630 ✗ } else if (n == 3) {
631 ✗ r_n3 = 1.0;
632 } else {
633 ✗ r_n3 = pow(1 / r, 3 - n);
634 }
635
636 // get the part of the gradient of Y,
637 // conjugate and multiply by the functional matrix
638 ✗ A = functionalMatrix(y);
639 ✗ z_harm = r_n3 * (A * z_harm);
640
641 // get the complex conjugation
642 ✗ z_harm.col(1) *= -1.0;
643 }
644
645 // calculate the radial part and multiply with n*r^{n-1}
646
1/2
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
 1203930 Eigen::Matrix<double, 3, 2> z_rad;
647
648 // case n = 0: r^n = 1 is constant
649
2/2
✓ Branch 0 taken 13230 times.
✓ Branch 1 taken 1190700 times.
 1203930 if (n == 0) {
650
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 z_rad.setZero();
651 } else {
652 // the gradient of r^n is n*r^(n-2)*y = nr^(n-1)*x
653 // again adapt the power
654 double nr_n1;
655
2/2
✓ Branch 0 taken 26460 times.
✓ Branch 1 taken 1164240 times.
 1190700 if (n == 1) {
656 26460 nr_n1 = 1.0;
657 } else {
658 1164240 nr_n1 = n * pow(r, n - 1);
659 }
660
661 double y_imt;
662
1/2
✓ Branch 0 taken 1190700 times.
✗ Branch 1 not taken.
 1190700 if (m >= 0) {
663 1190700 y_imt = y_im;
664 } else {
665 ✗ y_imt = -y_im;
666 }
667
668
3/6
✓ Branch 1 taken 1190700 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1190700 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1190700 times.
✗ Branch 8 not taken.
 1190700 z_rad.col(0) = y_re * nr_n1 * x;
669
3/6
✓ Branch 1 taken 1190700 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1190700 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1190700 times.
✗ Branch 8 not taken.
 1190700 z_rad.col(1) = y_imt * nr_n1 * x;
670 }
671
2/4
✓ Branch 1 taken 1203930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1203930 times.
✗ Branch 5 not taken.
 2407860 return z_harm + z_rad;
672 }
673
674 /**
675 * \brief Returns the values of the spherical polynomials \f$ P_n^m(t) \f$
676 */
677 13230 Eigen::VectorXd legendreFull(unsigned int N, double t) {
678
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 Eigen::VectorXd L(((N + 1) * (N + 2)) / 2);
679 int n, s, m;
680
681 /* n = 0 */
682
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 L(0) = 1;
683
684 /* n = 1 */
685
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 L(1) = t;
686
1/2
✓ Branch 1 taken 13230 times.
✗ Branch 2 not taken.
 13230 L(2) = 1;
687
688
2/2
✓ Branch 0 taken 145530 times.
✓ Branch 1 taken 13230 times.
 158760 for (n = 2; n <= N; n++) {
689 145530 s = (n * (n + 1)) / 2;
690
2/2
✓ Branch 0 taken 873180 times.
✓ Branch 1 taken 145530 times.
 1018710 for (m = 0; m < n - 1; m++) {
691
1/2
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
 873180 L(s + m) = ((2 * n - 1) * t * L((n * (n - 1)) / 2 + m)
692
2/4
✓ Branch 1 taken 873180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 873180 times.
✗ Branch 5 not taken.
 873180 - (n + m - 1) * L(((n - 1) * (n - 2)) / 2 + m)) / (n - m);
693 }
694
695 /* m = n-1 */
696 145530 m = n - 1;
697
2/4
✓ Branch 1 taken 145530 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 145530 times.
✗ Branch 5 not taken.
 145530 L(s + m) = (2 * m + 1) * t * L((m * (m + 1)) / 2 + m);
698
699 /* m = n */
700 145530 m = n;
701
2/4
✓ Branch 1 taken 145530 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 145530 times.
✗ Branch 5 not taken.
 145530 L(s + m) = (2 * m - 1) * L((n * (n - 1)) / 2 + n - 1);
702 }
703 13230 return L;
704 ✗ }
705
706 /**
707 * \brief gives the norming factor of the spherical polynomials
708 */
709 1203930 double constant(int m, unsigned int n) {
710 unsigned int i;
711 1203930 double c = sqrt((2 * n + 1) / (4 * pi));
712
713
2/2
✓ Branch 0 taken 4815720 times.
✓ Branch 1 taken 1203930 times.
 6019650 for (i = 0; i < m; i++) {
714 4815720 c *= 1 / sqrt((n + i + 1) * (n - i));
715 }
716
717 1203930 return c;
718 }
719
720 /**
721 * \brief gives \f$ |z|^3 \f$ times the Jacobi Matrix of the transformation \f$ z \mapsto \frac{z}{|z|} \f$
722 */
723 1265137 inline Eigen::Matrix3d functionalMatrix(Eigen::Vector3d z) {
724 1265137 Eigen::Matrix3d M;
725
726 1265137 double r = z.norm();
727
728 1265137 M(0, 0) = r * r - z(0) * z(0);
729 1265137 M(1, 1) = r * r - z(1) * z(1);
730 1265137 M(2, 2) = r * r - z(2) * z(2);
731
732 1265137 M(1, 0) = M(0, 1) = -z(0) * z(1);
733 1265137 M(2, 0) = M(0, 2) = -z(0) * z(2);
734 1265137 M(2, 1) = M(1, 2) = -z(1) * z(2);
735
736 1265137 return M;
737 }
738
739 /**
740 * \brief returns the \f$ n \f$-th power of \f$ x \f$ without using pow.
741 */
742 40785272 inline double pow_int(double x, int n) {
743
2/2
✓ Branch 0 taken 3401717 times.
✓ Branch 1 taken 37383555 times.
 40785272 if (n == 0) {
744 3401717 return 1.0;
745 }
746
747 unsigned int ul, k;
748
2/2
✓ Branch 0 taken 37199934 times.
✓ Branch 1 taken 183621 times.
 37383555 if (n > 0) {
749 37199934 ul = n;
750 } else {
751 183621 ul = -n;
752 }
753
754 37383555 double res = x;
755
2/2
✓ Branch 0 taken 171929021 times.
✓ Branch 1 taken 37383555 times.
 209312576 for (k = 1; k < ul; k++) {
756 171929021 res *= x;
757 }
758
759
2/2
✓ Branch 0 taken 37199934 times.
✓ Branch 1 taken 183621 times.
 37383555 if (n > 0) {
760 37199934 return res;
761 } else {
762 183621 return 1.0 / res;
763 }
764 }
765
766 } // namespace Bembel
767
768 #endif // BEMBEL_SRC_UTIL_SPHERICALS_HPP_
769