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