1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file maxwell3d.cpp
19 \author D. Martin, E. Lunéville
20 \since 09 oct 2010
21 \date 8 sept 2017
22
23 \brief Implementation of exact solutions of Maxwell 3d problems
24 */
25
26 #include "exactSolutions.hpp"
27 #include "../specialFunctions/specialFunctions.hpp"
28
29 namespace xlifepp
30 {
31
32 /*!
33 Exact solution to the Maxwell 3D problem outside a sphere :
34 Field scattered by solid sphere of radius Rad of an incoming vector field
35 Phi_w= (1,0,0)*exp^{(i*K*z}).
36 with electric field boundary condition Exn=0 on the sphere
37
38 Some computed values for K=1
39 X_1 X_2 X_3 Phi_1 Phi_2 Phi_3
40 1.0 0.0 0.0 | (1.75532926384 , 1.62107824182 ) (0 , 0 ) (-6.02796306309e-18,-1.58122117268e-18)
41 0.0 1.0 0.0 | (-1 , 4.81409564708e-17) (0 , 0 ) (0 , 0 )
42 0.0 0.0 1.0 | (-5.40302305868e-01,-8.41470984808e-01) (0 , 0 ) (0 , 0 )
43 1.0 1.0 1.0 | (-1.69184878965e-01, 1.61131192686e-01) (2.06589901137e-01, 2.52914553648e-01) (1.19848010392e-01, 3.28542217311e-01)
44 1.0 1.0 1.2 | (-1.86022424304e-01, 9.18508880491e-02) (1.49911272089e-01, 2.01342400004e-01) (9.37804387189e-02, 3.03355869671e-01)
45 1.0 1.0 1.4 | (-1.91775499206e-01, 3.34280295554e-02) (1.07142066815e-01, 1.59632591445e-01) (6.64556013889e-02, 2.71335010720e-01)
46 1.0 1.0 1.6 | (-1.88260942630e-01,-1.51180741005e-02) (7.54265673653e-02, 1.27220194457e-01) (4.11641109484e-02, 2.38176631400e-01)
47 1.0 1.0 1.8 | (-1.77261848423e-01,-5.51624601093e-02) (5.20175373132e-02, 1.02356587880e-01) (1.91974583103e-02, 2.06737287161e-01)
48 1.0 1.0 2.0 | (-1.60341338903e-01,-8.78244555191e-02) (3.47085585869e-02, 8.32071602821e-02) (8.24999561215e-04, 1.78106610047e-01)
49 1.2 1.4 1.0 | (-1.93790114095e-01, 7.12577504506e-02) (1.44884948487e-01, 2.06529528004e-01) (1.39532076671e-02, 1.97067304047e-01)
50 1.2 1.4 1.2 | (-1.90729686988e-01, 2.88712762800e-02) (1.12819966807e-01, 1.81439122245e-01) (9.10642428927e-03, 1.94073277948e-01)
51 1.2 1.4 1.4 | (-1.83137187727e-01,-1.04437244804e-02) (8.50160197427e-02, 1.57350273961e-01) (7.19355647585e-04, 1.84905785693e-01)
52 1.2 1.4 1.6 | (-1.70990951385e-01,-4.55518271208e-02) (6.17400731503e-02, 1.35416184077e-01) (-9.27523672141e-03, 1.71746716026e-01)
53 1.2 1.4 1.8 | (-1.54657449587e-01,-7.58575264387e-02) (4.27052811485e-02, 1.16028225273e-01) (-1.95184324637e-02, 1.56344795623e-01)
54 1.2 1.4 2.0 | (-1.34728108413e-01,-1.01074590219e-01) (2.73910156223e-02, 9.91350581457e-02) (-2.91430511034e-02, 1.39939464554e-01)
55 */
scatteredFieldMaxwellExn(const Point & p,Parameters & param)56 Vector<complex_t> scatteredFieldMaxwellExn(const Point& p, Parameters& param)
57 {
58 // We first recover the current value of the wave number K and radius of sphere
59 real_t k = 1., rad = 1.;
60 if (param.contains("k"))
61 {
62 k = real((param)("k"));
63 }
64 if (param.contains("radius"))
65 {
66 rad = real(param("radius"));
67 }
68 const real_t kr(k * rad);
69 const int nmax = 50;
70
71 // spherical Bessel functions up to order nmax for argument K*Rad
72 std::vector<real_t> js(nmax + 1), ys(nmax + 1);
73 js = sphericalbesselJ0N(kr, nmax);
74 ys = sphericalbesselY0N(kr, nmax);
75
76 // computation of coefficients an(K*Rad) and bn(K*Rad) until a minimum n is found such that
77 // | an(K*Rad) | <epsilon (the series is known to converge quickly)
78
79 // initialize n, 2n+1, i^n
80 const complex_t i1(complex_t(0., 1.));
81 int n = 0, nm1 = -1, np1 = 1;
82 real_t nx2p1 = 1.;
83 complex_t in = 1.;
84 std::vector<complex_t> a(nmax), b(nmax);
85 // H_{0}(K*Rad)
86 complex_t hn(js[0], ys[0]);
87 // Use recurrence relation to compute derivative : d/dz[H_{0}](z)=- H_{1}(z)
88 // d/dz[H_{0}](z)+H_{0}(z)=-H_{1}(z)+H_{0}(z)
89 complex_t hnpphn = -complex_t(js[1], ys[1]) + hn;
90 a[0] = -hnpphn.real() / hnpphn;
91 b[0] = -js[0] / hn;
92 while (std::abs(a[n]) > theZeroThreshold && n < nmax)
93 {
94 n++;
95 nx2p1 += 2.;
96 nm1++;
97 np1++;
98 hn = complex_t(js[n], ys[n]);
99 // Use recurrence relation to compute derivative d/dz[H_{n}](z) :
100 // (2n+1) d/dz[H_{n}](z)=n*H_{n-1}(z)-(n+1)*H_{n+1}(z)
101 // http://www.math.sfu.ca/~cbm/aands/page_439.htm (formula 10.1.20)
102 hnpphn = n * complex_t(js[nm1], ys[nm1]) - np1 * complex_t(js[np1], ys[np1]) + nx2p1 * hn;
103 a[n] = -in * nx2p1 * (hnpphn.real() / hnpphn);
104 b[n] = -in * nx2p1 * (js[n] / hn);
105 in *= i1;
106 }
107 int nf = n;
108 if (nf >= nmax)
109 {
110 where("scatteredFieldMaxwell_Exn");
111 error("conv_failed",nmax);
112 abort();
113 }
114 // spherical coordinates rho, theta and phi
115 real_t r(p[0]*p[0] + p[1]*p[1]), rho(std::sqrt(r + p[2]*p[2]));
116 r = std::sqrt(r);
117 real_t cost(p[2] / rho), sint(r / rho), cosp(1.), sinp(0.);
118 if (r > theEpsilon)
119 {
120 cosp = p[0] / r;
121 sinp = p[1] / r;
122 }
123
124 const dimen_t spaceDim(p.size());
125 Vector<complex_t> erho(spaceDim), etheta(spaceDim), ephi(spaceDim, 0.);
126 erho[0] = sint * cosp;
127 erho[1] = sint * sinp;
128 erho[2] = cost;
129 etheta[0] = cost * cosp;
130 etheta[1] = cost * sinp;
131 etheta[2] = -sint;
132 ephi[0] = -sinp;
133 ephi[1] = cosp;
134
135 // Compute partial sum of the series (including Legendre polynomials)
136 rho *= k;
137 // spherical Bessel functions up to order nf for argument K*rho
138 js = sphericalbesselJ0N(rho, nf);
139 ys = sphericalbesselY0N(rho, nf);
140
141 // initialize computation of Legendre polynomials and (2n+1)*cosp
142 real_t pcnm1 = 1., pcn = cost, pcnp1, n2p1cost = cost, twocost = 2 * cost, psnm1 = 0., psn = sint, psnp1, p0nm1 = 0., p0n = 1., p0np1, d0n = cost;
143 complex_t crho(0.), ctheta(0.), cphi(0.);
144 // initialize n-1, n+1, 2n+1, \sum_{1..n} (2i)
145 nm1 = 0;
146 np1 = 2;
147 nx2p1 = 3.;
148 real_t sumnx2 = 2.;
149 for (int n = 1; n < nf; nm1++, np1++, nx2p1 += 2, n++)
150 {
151 hn = complex_t(js[n], ys[n]);
152 // Use recurrence relation to compute derivative d/dz[H_{n}](z) :
153 // (2n+1) d/dz[H_{n}](z)=n*H_{n-1}(z)-(n+1)*H_{n+1}(z)
154 // to compute hnpphn=d/dz[H_{n}](K*rho)+H_{n}(K*rho)/(K*rho)
155 hnpphn = (n * complex_t(js[nm1], ys[nm1]) - np1 * complex_t(js[np1], ys[np1])) / nx2p1 + hn / rho;
156 crho += (a[n] * hn / rho) * psn;
157 hn *= i1 * b[n];
158 hnpphn *= a[n];
159 ctheta += (hnpphn * d0n + hn * p0n) / sumnx2;
160 cphi -= (hnpphn * p0n + hn * d0n) / sumnx2;
161
162 // Recurrence relation for Legendre polynomials
163 // (n+1) P_{n+1}(cosp)=(2n+1)*cosp*P_{n}(cosp)-n*P_{n-1}(cosp)
164 n2p1cost += twocost;
165 sumnx2 += 2 * np1;
166 pcnp1 = (n2p1cost * pcn - n * pcnm1) / np1;
167 pcnm1 = pcn;
168 pcn = pcnp1;
169 psnp1 = (n2p1cost * psn - np1 * psnm1) / n;
170 psnm1 = psn;
171 psn = psnp1;
172 p0np1 = (n2p1cost * p0n - np1 * p0nm1) / n;
173 p0nm1 = p0n;
174 p0n = p0np1;
175 d0n = sumnx2 * pcn - cost * p0n;
176 }
177 return (cosp * (crho * erho + ctheta * etheta) + sinp * cphi * ephi);
178 }
179
180 /*!
181 This function computes spherical harmonics Y_{l}^{m} up to order n on the unit sphere;
182 order n is computed by n=ylm.size()-1 and spherical harmonics ylm are computed for l=0,..,n and m=0,..,l.
183
184 Y_{l}^{m}=\f$(-1)^{(|m|-m) /2} \sqrt{ (2l+1)/4\pi*(l-|m|)!/(l+|m|)! } P_{l}^{|m|}(\cos(\theta)) e^{im \phi}\f$
185
186 which yields for m positive
187
188 Y_{l}^{m}=\f$\sqrt{ (2l+1)/4\pi*(l-m)!/(l+m)! } P_{l}^{m}(\cos(\theta)) e^{im \phi}\f$
189
190 Argument ylm is a "triangular 2d array" defined as a std::vector of vectors
191
192 the 1st std::vector carries the values : Y_{0}^{0}
193 the 2nd one carries the values : Y_{1}^{0} Y_{1}^{1}
194 and for l=2,..,n
195 the l-th std::vector carries the values: Y_{l}^{0} Y_{l}^{1} .. Y_{l}^{l}
196 if the size of std::vector ylm[l] is l+1.
197
198 For a given int n, container ylm can be declared (and defined) prior
199 to call to this function by
200 std::vector<std::vector<complex_t> > ylm(n+1);
201 for (int l=0; l<=n; l++) ylm[l]=std::vector<complex_t>(l+1);
202 */
sphericalHarmonics(const Point & x,std::vector<std::vector<complex_t>> & ylm)203 void sphericalHarmonics(const Point& x, std::vector<std::vector<complex_t> >& ylm)
204 {
205 // spherical coordinates rho, theta and phi
206 real_t rho(std::sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]));
207 real_t cost = x[2] / rho; //=cos(theta)
208 real_t phi = 0.;
209 if (std::abs(x[0]) > theEpsilon)
210 {
211 phi = std::atan(x[1] / x[0]);
212 }
213 //
214 // Compute Legendre functions: P_l^m (m=0,..,l; l=0,..,n)
215 //
216 int n = ylm.size() - 1;
217 // define 2d triangular array container for Legendre functions
218 std::vector<std::vector<real_t> > plm(n + 1);
219 for (int m = 0; m <= n; m++)
220 {
221 plm[m] = std::vector<real_t>(n + 1 - m);
222 }
223 // compute associated Legendre functions
224 legendreFunctions(cost, plm);
225
226 // Compute normalized spherical harmonics: Y_l^m (m=0,..,l; l=0,..,n)
227 std::vector<std::vector<complex_t> >::iterator it_yl(ylm.begin());
228 std::vector<complex_t>::iterator it_l, it_ylm;
229
230 int l = 0;
231 for (it_yl = ylm.begin(); it_yl != ylm.end(); it_yl++, l++)
232 {
233 real_t factRatio(1.), over4pi2lp1((2 * l + 1)*over4pi_);
234 // compute Y_{l}^{0}
235 *((*it_yl).begin()) = std::sqrt(over4pi2lp1) * plm[0][l];
236 int m(1);
237 if (l > 0)
238 {
239 factRatio = 1.;
240 int lpm(l + 1), lpm1(l);
241 for (it_ylm = (*it_yl).begin() + 1; it_ylm != (*it_yl).end() - 1; m++, lpm++, lpm1--, it_ylm++)
242 {
243 // compute ratio of factorials: (l-m)!/(l+m)! for m=1,..., l-1
244 factRatio /= lpm * lpm1;
245 // compute Y_{l}^{m} for m=1,..., l-1
246 *it_ylm = std::sqrt(over4pi2lp1 * factRatio) * plm[m][l - m] * std::exp(complex_t(0, m) * phi);
247 }
248 // compute Y_{l}^{l} for l > 0
249 factRatio /= 2 * l;
250 *((*it_yl).end() - 1) = std::sqrt(over4pi2lp1 * factRatio) * plm[l][0] * std::exp(complex_t(0, l) * phi);
251 }
252 }
253 }
254
255 /*!
256 This function computes spherical harmonics surface gradients gradY_{l}^{m} up to order n on the unit sphere;
257 order n is computed by n=ylm.size()-1 and spherical harmonics ylm are
258 computed for l=0,..,n and m=0,..,l.
259
260 Argument gradylm is a "triangular 2d array of 3d-vectors" defined as a vector of vectors of complex_t-Vectors.
261 the 1st std::vector carries the 3d-vectorial values : gradY_{0}^{0}
262 the 2nd one carries the 3d-vectorial values : gradY_{1}^{0} gradY_{1}^{1}
263 and for l=2,..,n
264 the l-th std::vector carries the 3d-vectorial values: gradY_{l}^{0} gradY_{l}^{1} .. gradY_{l}^{l}
265 where the size of std::vector gradylm[l] is l+1.
266
267 For a given int n, container gradylm can be declared (and defined) prior
268 to call to this function by
269 std::vector<std::vector<Vector<complex_t> > > gradylm(n+1);
270 for (int l=0; l<=n; l++) {
271 gradylm[l]=std::vector<Vector<complex_t> >(l+1);
272 for (int m=0; m<=l; m++) gradylm[l][m]=Vector<complex_t>(spaceDim);
273 }
274 where "spaceDim" is previously defined as the dimension space (necessarily 3 here).
275 TELL DANIEL : We may put a test that stops the program if spaceDim is not 3.
276 */
sphericalHarmonicsSurfaceGrad(const Point & x,std::vector<std::vector<Vector<complex_t>>> & gradylm)277 void sphericalHarmonicsSurfaceGrad(const Point& x, std::vector<std::vector<Vector<complex_t> > >& gradylm)
278 {
279 // spherical coordinates rho, theta and phi
280 real_t rho = std::sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
281 real_t cost = x[2] / rho; //=cos(theta)
282 real_t sint = std::sqrt(1. - cost * cost);
283 real_t phi(0.), cosp(1.), sinp(0.);
284 if (std::abs(x[0]) > theEpsilon)
285 {
286 phi = std::atan(x[1] / x[0]);
287 cosp = std::cos(phi);
288 sinp = std::sin(phi);
289 }
290
291 const dimen_t spaceDim(x.size());
292 Vector<complex_t> etheta(spaceDim), ephi(spaceDim, 0.);
293 etheta[0] = cost * cosp;
294 etheta[1] = cost * sinp;
295 etheta[2] = -sint;
296 ephi[0] = -sinp;
297 ephi[1] = cosp;
298
299 // Compute Legendre functions: P_l^m (m=0,..,l; l=0,..,n)
300 int n = gradylm.size() - 1;
301 // define 2d triangular array container for Legendre functions
302 std::vector<std::vector<real_t> > plm(n + 1);
303 for (int m = 0; m <= n; m++)
304 {
305 plm[m] = std::vector<real_t>(n + 1 - m);
306 }
307 // compute associated Legendre functions
308 legendreFunctions(cost, plm);
309
310 // Compute derivative of Legendre functions: P_l^m (m=0,..,l; l=0,..,n)
311 real_t factRatio(1.), over4pi2lp1;
312 complex_t coeffPhi;// coeff_theta, coeff
313 // define 2d triangular array container for Legendre functions derivatives
314 std::vector<std::vector<real_t> > dplm(n + 1);
315 for (int l = 0; l <= n; l++)
316 {
317 dplm[l] = std::vector<real_t>(l + 1);
318 }
319 // compute derivative of associated Legendre functions
320 legendreFunctionsDerivative(cost, plm, dplm);
321 std::vector<std::vector<real_t> >::iterator it_dpl(dplm.begin());
322 std::vector<real_t>::iterator it_dplm;
323
324 // Compute normalized spherical harmonics derivatives: gradY_l^m (m=0,..,l; l=0,..,n)
325 std::vector<std::vector<Vector<complex_t> > >::iterator it_gradyl(gradylm.begin());
326 std::vector<Vector<complex_t> >::iterator it_gradylm;
327
328 if (std::abs(sint) < theEpsilon)
329 {
330 *((*it_gradyl).begin()) = complex_t(0.); // Case l=0 where gradylm is equal to 0.
331 int l = 1;
332 real_t costl1 = cost;
333 it_dpl = dplm.begin() + 1;
334 for (it_gradyl = gradylm.begin() + 1; it_gradyl != gradylm.end(); it_gradyl++, l++, it_dpl++)
335 {
336 int llp1 = l * (l + 1);
337 factRatio = 1. / llp1;
338 over4pi2lp1 = (2 * l + 1) * over4pi_;
339 *((*it_gradyl).begin()) = complex_t(0.);
340 costl1 *= cost;
341
342 // compute gradY_{l}^{1} -- the other are 0
343 // gradY_{l}^{1}=\sqrt{ (2l+1)/4*pi*(l-1)!/(l+1)! }*e^{i phi} *
344 // (P'_{l}^{1}(\cos(theta)) e_theta-i (cos(theta))^{l+1} cos(phi) l(l+1)/2*e_phi)
345 // where e_theta and e_phi are the unit tangent std::vector on the unit sphere.
346 //
347 it_gradylm = (*it_gradyl).begin() + 1;
348 it_dplm = (*it_dpl).begin() + 1;
349 coeffPhi = - complex_t(0, 0.5 * costl1 * cosp * llp1);
350 *it_gradylm = std::sqrt(over4pi2lp1 * factRatio) * std::exp(complex_t(0, 1) * phi) * (*it_dplm * etheta + coeffPhi * ephi);
351 for (it_gradylm = (*it_gradyl).begin() + 2; it_gradylm != (*it_gradyl).end(); it_gradylm++)
352 {
353 *it_gradylm = complex_t(0.); // Case l>1 where gradylm is 0.
354 }
355 }
356 }
357 else
358 {
359 *((*it_gradyl).begin()) = complex_t(0.); // Case l=0 where gradylm is 0.
360 int l = 1;
361 it_dpl = dplm.begin() + 1;
362 for (it_gradyl = gradylm.begin() + 1; it_gradyl != gradylm.end(); it_gradyl++, l++, it_dpl++)
363 {
364 int lpm = l + 1, lpm1 = l, m = 0;
365 factRatio = 1.;
366 over4pi2lp1 = (2 * l + 1) * over4pi_;
367 it_dplm = (*it_dpl).begin();
368 for (it_gradylm = (*it_gradyl).begin(); it_gradylm != (*it_gradyl).end() - 1; m++, lpm++, lpm1--, it_gradylm++, it_dplm++)
369 {
370 // compute gradY_{l}^{m} for m=1,..., l
371 // gradY_{l}^{m}=sqrt{ (2l+1)/4*pi*(l-m)!/(l+m)! }*e^{i m phi} *
372 // (P'_{l}^{|m|}(cos(theta)) e_theta+im*P_{l}^{|m|}(cos(theta)) e_phi)
373 // where e_theta and e_phi are the unit tangent std::vector on the unit sphere.
374 //
375 // coeff=sqrt(over4pi2lp1*factRatio)*exp(complex_t(0,m)*phi);
376 // coeff_theta=*it_dplm;
377 coeffPhi = complex_t(0, m) * plm[m][l - m] / sint;
378 *it_gradylm = std::sqrt(over4pi2lp1 * factRatio) * std::exp(complex_t(0, m) * phi) * (*it_dplm * etheta + coeffPhi * ephi);
379 // compute ratio of factorials: (l-(m+1))!/(l+(m+1))! for m=O,..., l-1
380 factRatio /= lpm * lpm1;
381 }
382 it_gradylm = (*it_gradyl).end() - 1;
383 it_dplm = (*it_dpl).end() - 1;
384 coeffPhi = complex_t(0, l * plm[l][0] / sint);
385 *it_gradylm = std::sqrt(over4pi2lp1 * factRatio) * std::exp(complex_t(0, l) * phi) * (*it_dplm * etheta + coeffPhi * ephi);
386 }
387 }
388 }
389
390 // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
sphericalHarmonicsTest(const Point & x,const number_t n,std::ostream & out)391 void sphericalHarmonicsTest(const Point& x, const number_t n, std::ostream& out)
392 {
393 // Output function for Spherical Harmonics
394 std::vector<std::vector<complex_t> > yml(n + 1);
395 for (number_t l = 0; l <= yml.size() - 1; l++)
396 {
397 yml[l] = std::vector<complex_t>(l + 1);
398 }
399 sphericalHarmonics(x, yml);
400 out << " Spherical harmonics Y_l^m up to order " << yml.size() - 1 << " (m=0,..,l; l=0,..," << yml.size() - 1 << ")";
401 out << " for point(" << x[0] << ", " << x[1] << ", " << x[2] << ")" << std::endl;
402 out.setf(std::ios::scientific);
403 int l = 0;
404 for (std::vector<std::vector<complex_t> >::iterator it_yl = yml.begin(); it_yl != yml.end(); it_yl++, l++)
405 {
406 int m = 0;
407 out << " l=" << l << std::endl;
408 for (std::vector<complex_t>::iterator it_yml = (*it_yl).begin(); it_yml != (*it_yl).end(); it_yml++, m++)
409 {
410 out << " Y_" << l << "^" << m << "=" << std::setw(19) << std::setprecision(12) << std::endl;
411 }
412 }
413 out.unsetf(std::ios::scientific);
414 }
415
sphericalHarmonicsSurfaceGradTest(const Point & x,const number_t n,std::ostream & out)416 void sphericalHarmonicsSurfaceGradTest(const Point& x, const number_t n, std::ostream& out)
417 {
418 // Output function for Spherical Harmonics surface gradients
419 const dimen_t spaceDim(x.size());
420 std::vector<std::vector<Vector<complex_t> > > gradyml(n + 1);
421 for (number_t l = 0; l <= gradyml.size() - 1; l++)
422 {
423 gradyml[l] = std::vector<Vector<complex_t> >(l + 1);
424 for (number_t m = 0; m <= gradyml[l].size() - 1; m++)
425 {
426 gradyml[l][m] = Vector<complex_t>(spaceDim);
427 }
428 }
429 sphericalHarmonicsSurfaceGrad(x, gradyml);
430
431 out << " Spherical harmonics derivatives gradY_l^m up to order " << gradyml.size() - 1 << " (m=0,..,l; l=0,..," << gradyml.size() - 1 << ")";
432 out << " for point(" << x[0] << ", " << x[1] << ", " << x[2] << ")" << std::endl;
433 out.setf(std::ios::scientific);
434 int l = 0;
435 for (std::vector<std::vector<Vector<complex_t> > >::iterator it_gradyl = gradyml.begin(); it_gradyl != gradyml.end(); it_gradyl++, l++)
436 {
437 int m = 0;
438 out << " l=" << l << std::endl;
439 for (std::vector<Vector<complex_t> >::iterator grad_it_yml = (*it_gradyl).begin(); grad_it_yml != (*it_gradyl).end(); grad_it_yml++, m++)
440 {
441 out << " gradY_" << l << "^" << m << "=("
442 << std::setw(19) << std::setprecision(12) << (*grad_it_yml)[0] << " , " << (*grad_it_yml)[1] << " , " << (*grad_it_yml)[2] << "). " << std::endl;
443 }
444 }
445 out.unsetf(std::ios::scientific);
446 }
447
448 } // end of namespace xlifepp
449
450