1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5 //
6 // Contributor(s):
7 // Ruth Sabariego
8 // Xavier Antoine
9 //
10
11 // g++ -std=c++11 on mingw does not define bessel functions
12 #if defined (WIN32) && !defined(__CYGWIN__)
13 #undef __STRICT_ANSI__
14 #endif
15
16 #include <math.h>
17 #include <stdlib.h>
18 #include "ProData.h"
19 #include "F.h"
20 #include "Legendre.h"
21 #include "Bessel.h"
22 #include "MallocUtils.h"
23 #include "Message.h"
24
25 #define SQU(a) ((a)*(a))
26
27 /* some utility functions to deal with complex numbers */
28
29 typedef struct {
30 double r;
31 double i;
32 } cplx;
33
Csum(cplx a,cplx b)34 static cplx Csum(cplx a, cplx b)
35 {
36 cplx s;
37 s.r = a.r + b.r;
38 s.i = a.i + b.i;
39 return(s);
40 }
41
Csub(cplx a,cplx b)42 static cplx Csub(cplx a, cplx b)
43 {
44 cplx s;
45 s.r = a.r - b.r;
46 s.i = a.i - b.i;
47 return(s);
48 }
49
Csubr(double a,cplx b)50 static cplx Csubr(double a, cplx b)
51 {
52 cplx s;
53 s.r = a - b.r;
54 s.i = - b.i;
55 return(s);
56 }
57
Cprod(cplx a,cplx b)58 static cplx Cprod(cplx a, cplx b)
59 {
60 cplx s;
61 s.r = a.r * b.r - a.i * b.i;
62 s.i = a.r * b.i + a.i * b.r;
63 return(s);
64 }
65
Cdiv(cplx a,cplx b)66 static cplx Cdiv(cplx a, cplx b)
67 {
68 cplx s;
69 double den;
70 den = b.r * b.r + b.i * b.i;
71 s.r = (a.r * b.r + a.i * b.i) / den;
72 s.i = (a.i * b.r - a.r * b.i) / den;
73 return(s);
74 }
75
Cdivr(double a,cplx b)76 static cplx Cdivr(double a, cplx b)
77 {
78 cplx s;
79 double den;
80 den = b.r * b.r + b.i * b.i;
81 s.r = (a * b.r) / den;
82 s.i = (- a * b.i) / den;
83 return(s);
84 }
85
Cconj(cplx a)86 static cplx Cconj(cplx a)
87 {
88 cplx s;
89 s.r = a.r;
90 s.i = -a.i;
91 return(s);
92 }
93
Cneg(cplx a)94 static cplx Cneg(cplx a)
95 {
96 cplx s;
97 s.r = -a.r;
98 s.i = -a.i;
99 return(s);
100 }
101
Cmodu(cplx a)102 static double Cmodu(cplx a)
103 {
104 return(sqrt(a.r * a.r + a.i * a.i));
105 }
106
Cpow(cplx a,double b)107 static cplx Cpow(cplx a, double b)
108 {
109 cplx s;
110 double mod, arg;
111 mod = a.r * a.r + a.i * a.i;
112 arg = atan2(a.i,a.r);
113 mod = pow(mod,0.5*b);
114 arg *= b;
115 s.r = mod * cos(arg);
116 s.i = mod * sin(arg);
117
118 return(s);
119 }
120
Cprodr(double a,cplx b)121 static cplx Cprodr(double a, cplx b)
122 {
123 cplx s;
124 s.r = a * b.r;
125 s.i = a * b.i;
126 return(s);
127 }
128
129 /* ------------------------------------------------------------------------ */
130 /* Exact solutions for spheres */
131 /* ------------------------------------------------------------------------ */
132
133 /* Scattering by solid PEC sphere. Returns theta-component of surface
134 current */
135
F_JFIE_SphTheta(F_ARG)136 void F_JFIE_SphTheta(F_ARG)
137 {
138 double k0, r, kr, e0, eta, theta, phi, a1, b1, c1, d1, den1, P, P0, dP ;
139 double ctheta, stheta, cteRe1, cteRe2, a2, b2, c2, d2, den2 ;
140 int i, n ;
141
142 theta = atan2(sqrt( A->Val[0]* A->Val[0] + A->Val[1]*A->Val[1] ), A->Val[2]);
143 phi = atan2( A->Val[1], A->Val[0] ) ;
144
145 k0 = Fct->Para[0] ;
146 eta = Fct->Para[1] ;
147 e0 = Fct->Para[2] ;
148 r = Fct->Para[3] ;
149
150 kr = k0*r ;
151 n = 50 ;
152
153 V->Val[0] = 0.;
154 V->Val[MAX_DIM] = 0. ;
155
156 if ( theta == 0. ) theta += 1e-7; /* Warning! This is an approximation. */
157 if ( theta == M_PI || theta == -M_PI ) theta -= 1e-7;
158
159 for (i = 1 ; i <= n ; i++ ){
160 ctheta = cos(theta);
161 stheta = sin(theta);
162
163 P = Legendre(i,1,ctheta);
164 P0 = Legendre(i,0,ctheta);
165
166 dP = (i+1)*i* P0/stheta-(ctheta/(ctheta*ctheta-1))* P;
167
168 cteRe1 = (2*i+1) * stheta * dP/i/(i+1);
169 cteRe2 = (2*i+1) * P/stheta/i/(i+1);
170
171 a1 = cos((1-i)*M_PI/2) ;
172 b1 = sin((1-i)*M_PI/2) ;
173 c1 = -AltSpherical_j_n(i+1, kr) + (i+1) * AltSpherical_j_n(i, kr)/kr ; /* Derivative */
174 d1 = -(-AltSpherical_y_n(i+1, kr) + (i+1) * AltSpherical_y_n(i, kr)/kr) ;
175
176 a2 = cos((2-i)*M_PI/2) ;
177 b2 = sin((2-i)*M_PI/2) ;
178 c2 = AltSpherical_j_n(i, kr) ;
179 d2 = -AltSpherical_y_n(i, kr) ;
180
181 den1 = c1*c1+d1*d1 ;
182 den2 = c2*c2+d2*d2 ;
183
184 V->Val[0] += cteRe1*(a1*c1+b1*d1)/den1 + cteRe2*(a2*c2+b2*d2)/den2 ;
185 V->Val[MAX_DIM] += cteRe1*(b1*c1-a1*d1)/den1 + cteRe2*(b2*c2-a2*d2)/den2 ;
186 }
187
188 V->Val[0] *= e0*cos(phi)/eta/kr ;
189 V->Val[MAX_DIM] *= e0*cos(phi)/eta/kr ;
190
191 V->Type = SCALAR ;
192 }
193
194 /* Scattering by solid PEC sphere. Returns theta-component of RCS */
195
F_RCS_SphTheta(F_ARG)196 void F_RCS_SphTheta(F_ARG)
197 {
198 double k0, r, kr, e0, rinf, krinf, theta, phi, a1 =0., b1=0., d1, den1, P, P0, dP ;
199 double J, J_1, dJ, ctheta, stheta, cteRe1, cteRe2, a2, b2, d2, den2, lambda ;
200 int i, n ;
201
202 theta = atan2(sqrt( A->Val[0]* A->Val[0] + A->Val[1]*A->Val[1] ), A->Val[2]);
203 phi = atan2( A->Val[1], A->Val[0] ) ;
204
205 k0 = Fct->Para[0] ;
206 e0 = Fct->Para[1] ;
207 r = Fct->Para[2] ;
208 rinf = Fct->Para[3] ;
209
210 kr = k0*r ;
211 krinf = k0*rinf ;
212 lambda = 2*M_PI/k0 ;
213
214 n = 50 ;
215
216 if ( theta == 0. ) theta += 1e-7; /* Warning! This is an approximation. */
217 if ( theta == M_PI || theta == -M_PI ) theta -= 1e-7;
218
219 for (i = 1 ; i <= n ; i++ ){
220 ctheta = cos(theta);
221 stheta = sin(theta);
222
223 P = Legendre(i,1,ctheta);
224 P0 = Legendre(i,0,ctheta);
225 dP = (i+1)*i* P0/stheta-(ctheta/(ctheta*ctheta-1))* P;
226
227 J = AltSpherical_j_n(i, kr) ;
228 J_1 = AltSpherical_j_n(i+1, kr) ;
229 dJ = -J_1 + (i + 1) * J/kr ;
230
231 cteRe1 = -(2*i+1) * stheta * dP * dJ /i/(i+1);
232 cteRe2 = (2*i+1) * P * J /stheta/i/(i+1);
233
234 d1 = -(-AltSpherical_y_n(i+1, kr) + (i+1) * AltSpherical_y_n(i, kr)/kr) ;
235
236 d2 = -AltSpherical_y_n(i, kr) ;
237
238 den1 = dJ*dJ+d1*d1 ;
239 den2 = J*J+d2*d2 ;
240
241 a1 += cteRe1 * dJ /den1 + cteRe2 * J /den2 ;
242 b1 += cteRe1*(-d1) /den1 + cteRe2*(-d2) /den2 ;
243 }
244
245 a2 = e0*cos(phi)*sin(krinf)/krinf ;
246 b2 = e0*cos(phi)*cos(krinf)/krinf ;
247
248 V->Val[0] = 10*log10( 4*M_PI*SQU(rinf/lambda)*(SQU(a1*a2-b1*b2) + SQU(a1*b2+a2*b1)) );
249 V->Val[MAX_DIM] = 0. ;
250
251 V->Type = SCALAR ;
252 }
253
254 /* Scattering by solid PEC sphere. Returns phi-component of surface current */
255
F_JFIE_SphPhi(F_ARG)256 void F_JFIE_SphPhi(F_ARG)
257 {
258 double k0, r, kr, e0, eta, theta, phi, a1, b1, c1, d1, den1, P, P0, dP ;
259 double ctheta, stheta, cteRe1, cteRe2, a2, b2, c2, d2, den2 ;
260 int i, n ;
261
262 theta = atan2( sqrt( A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1] ), A->Val[2]);
263 phi = atan2( A->Val[1], A->Val[0] ) ;
264
265 k0 = Fct->Para[0] ;
266 eta = Fct->Para[1] ;
267 e0 = Fct->Para[2] ;
268 r = Fct->Para[3] ;
269
270 kr = k0*r ;
271 n = 50 ;
272
273 V->Val[0] = 0.;
274 V->Val[MAX_DIM] = 0. ;
275
276 if ( theta == 0. ) theta += 1e-7; /* Warning! This is an approximation. */
277 if ( theta == M_PI || theta == -M_PI ) theta -= 1e-7;
278
279 for (i = 1 ; i <= n ; i++ ){
280 ctheta = cos(theta);
281 stheta = sin(theta);
282
283 P = Legendre(i,1,ctheta);
284 P0 = Legendre(i,0,ctheta);
285
286 dP = (i+1)*i* P0/stheta - ctheta/(ctheta*ctheta-1)*P; /* Derivative */
287
288 cteRe1 = (2*i+1) * P /i/(i+1)/stheta;
289 cteRe2 = (2*i+1) * stheta * dP/i/(i+1);
290
291 a1 = cos((1-i)*M_PI/2) ;
292 b1 = sin((1-i)*M_PI/2) ;
293 c1 = -AltSpherical_j_n(i+1, kr) + (i+1)*AltSpherical_j_n(i, kr)/kr ; /* Derivative */
294 d1 = -(-AltSpherical_y_n(i+1, kr) + (i+1)*AltSpherical_y_n(i, kr)/kr) ;
295
296 a2 = cos((2-i)*M_PI/2) ;
297 b2 = sin((2-i)*M_PI/2) ;
298 c2 = AltSpherical_j_n(i, kr) ;
299 d2 = -AltSpherical_y_n(i, kr) ;
300
301 den1 = c1*c1+d1*d1 ;
302 den2 = c2*c2+d2*d2 ;
303
304 V->Val[0] += cteRe1*(a1*c1+b1*d1)/den1 + cteRe2*(a2*c2+b2*d2)/den2 ;
305 V->Val[MAX_DIM] += cteRe1*(b1*c1-a1*d1)/den1 + cteRe2*(b2*c2-a2*d2)/den2 ;
306 }
307
308 V->Val[0] *= e0*sin(phi)/eta/kr ;
309 V->Val[MAX_DIM] *= e0*sin(phi)/eta/kr ;
310
311 V->Type = SCALAR ;
312 }
313
314 /* Scattering by solid PEC sphere. Returns phi-component of RCS */
315
F_RCS_SphPhi(F_ARG)316 void F_RCS_SphPhi(F_ARG)
317 {
318 double k0, r, kr, e0, rinf, krinf, theta, phi, a1 =0., b1=0., d1, den1, P, P0, dP ;
319 double J, J_1, dJ, ctheta, stheta, cteRe1, cteRe2, a2, b2, d2, den2, lambda ;
320 int i, n ;
321
322 theta = atan2(sqrt( A->Val[0]* A->Val[0] + A->Val[1]*A->Val[1] ), A->Val[2]);
323 phi = M_PI/2 ;
324
325 k0 = Fct->Para[0] ;
326 e0 = Fct->Para[1] ;
327 r = Fct->Para[2] ;
328 rinf = Fct->Para[3] ;
329
330 kr = k0*r ;
331 krinf = k0*rinf ;
332 lambda = 2*M_PI/k0 ;
333
334 n = 50 ;
335
336 if ( theta == 0. ) theta += 1e-7; /* Warning! This is an approximation. */
337 if ( theta == M_PI || theta == -M_PI ) theta -= 1e-7;
338
339 for (i = 1 ; i <= n ; i++ ){
340 ctheta = cos(theta);
341 stheta = sin(theta);
342
343 P = Legendre(i,1,ctheta);
344 P0 = Legendre(i,0,ctheta);
345 dP = (i+1)*i* P0/stheta-(ctheta/(ctheta*ctheta-1))* P;
346
347 J = AltSpherical_j_n(i, kr) ;
348 J_1 = AltSpherical_j_n(i+1, kr) ;
349 dJ = -J_1 + (i + 1) * J/kr ;
350
351 cteRe1 = -(2*i+1) * P * dJ /stheta/i/(i+1);
352 cteRe2 = (2*i+1) * stheta * dP * J/i/(i+1);
353
354 d1 = -(-AltSpherical_y_n(i+1, kr) + (i+1) * AltSpherical_y_n(i, kr)/kr) ;
355 d2 = -AltSpherical_y_n(i, kr) ;
356
357 den1 = dJ*dJ+d1*d1 ;
358 den2 = J*J+d2*d2 ;
359
360 a1 += cteRe1 * dJ /den1 + cteRe2 * J /den2 ;
361 b1 += cteRe1*(-d1) /den1 + cteRe2*(-d2) /den2 ;
362 }
363
364 a2 = e0*sin(phi)*sin(krinf)/krinf ;
365 b2 = e0*sin(phi)*cos(krinf)/krinf ;
366
367 V->Val[0] = 10*log10( 4*M_PI*SQU(rinf/lambda)*(SQU(a1*a2-b1*b2) + SQU(a1*b2+a2*b1)) );
368 V->Val[MAX_DIM] = 0. ;
369
370 V->Type = SCALAR ;
371 }
372
373 /* Scattering by a perfectly conducting sphere of radius a, under plane wave
374 incidence pol*e^{ik \alpha\dot\r}, with alpha = (0,0,-1) and pol =
375 (1,0,0). Returns the scattered electric field anywhere outside the sphere
376 (From Balanis, Advanced Engineering Electromagnetics, sec 11.8, p. 660)
377
378 Warning This is probably wring :-)
379
380 */
381
382
383 // diffraction onde plane par une sphere diectrique en -iwt
F_ElectricFieldDielectricSphereMwt(F_ARG)384 void F_ElectricFieldDielectricSphereMwt(F_ARG)
385 {
386 double x = A->Val[0];
387 double y = A->Val[1];
388 double z = A->Val[2];
389 double r = sqrt(x * x + y * y + z * z);
390 double theta = atan2(sqrt(x * x + y * y), z);
391 double phi = atan2(y, x);
392
393
394 double k = Fct->Para[0] ;
395 double a = Fct->Para[1] ;
396 double kr = k * r;
397 double ka = k * a;
398
399 double epsi = 2.25;
400 double ki = k*sqrt(epsi);
401 double Zi = 1.0 / sqrt(epsi);
402 double kia = ki * a;
403 double kir = ki * r;
404
405 int ns = (int)k + 12;
406
407 std::vector<std::complex<double> > Hnkr(ns + 1), Hnka(ns + 1), Hnkir(ns + 1), Hnkia(ns + 1);
408 for (int n = 1 ; n <= ns ; n++){
409 Hnkr[n] = std::complex<double>(AltSpherical_j_n(n, kr), AltSpherical_y_n(n, kr));
410 Hnka[n] = std::complex<double>(AltSpherical_j_n(n, ka), AltSpherical_y_n(n, ka));
411 Hnkia[n] = std::complex<double>(AltSpherical_j_n(n, kia), AltSpherical_y_n(n, kia));
412 Hnkir[n] = std::complex<double>(AltSpherical_j_n(n, kir), AltSpherical_y_n(n, kir));
413 }
414 double ctheta = cos(theta);
415 double stheta = sin(theta);
416
417 std::complex<double> Er(0.,0), Etheta(0.,0), Ephi(0.,0), I(0, 1.);
418
419
420
421
422 if( theta == 0.) {
423
424 if (r >= a ) {
425 for (int n = 1 ; n < ns ; n++){
426 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
427
428 double A1 = -n * (n + 1.) / 2.;
429
430 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
431 std::complex<double> dHnkia = - Hnkia[n + 1] + (n + 1.) * Hnkia[n] / kia;
432
433 std::complex<double> aln = (dHnka.real() * Hnkia[n].real() - Zi * dHnkia.real() * Hnka[n].real() ) /
434 (Zi*dHnkia.real() * Hnka[n] - dHnka * Hnkia[n].real()) ;
435 std::complex<double> ben = (dHnkia.real() * Hnka[n].real() - Zi * dHnka.real() * Hnkia[n].real() ) /
436 (Zi*Hnkia[n].real() * dHnka - dHnkia.real() * Hnka[n]) ;
437
438 std::complex<double> dn = aln*an;
439
440 std::complex<double> dHnkr = - Hnkr[n + 1] + (n + 1.) * Hnkr[n] / kr;
441 std::complex<double> d2Hnkr = Hnkr[n] / kr;
442 std::complex<double> jnkr = Hnkr[n].real() / kr;
443
444 double Pn1 = Legendre(n, 1, ctheta);
445
446 Er += (n * (n + 1.) * d2Hnkr * dn + n * (n + 1.) * jnkr * an ) * Pn1;
447 Etheta += an * (I * aln * dHnkr * A1 + ben * Hnkr[n] * A1 + I * dHnkr.real() * A1 + Hnkr[n].real() * A1 );
448 Ephi += an * (I * aln * dHnkr * A1 + ben * Hnkr[n] * A1 + I * dHnkr.real() * A1 + Hnkr[n].real() * A1 );
449 }
450
451 Er *= I * cos(phi) / kr;
452 Etheta *= (1. / kr) * (cos(phi));
453 Ephi *= (-1. / kr) * (sin(phi));
454 }
455 else {
456 for (int n = 1 ; n < ns ; n++){
457 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
458
459 double A1 = -n * (n + 1.) / 2.;
460
461 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
462 std::complex<double> dHnkia = - Hnkia[n + 1] + (n + 1.) * Hnkia[n] / kia;
463 std::complex<double> dHnkir = - Hnkir[n + 1] + (n + 1.) * Hnkir[n] / kir;
464
465 std::complex<double> den = (ki * Zi / k) * (dHnka * Hnka[n].real() - dHnka.real() * Hnka[n] ) /
466 (Zi*Hnkia[n].real() * dHnka - dHnkia.real() * Hnka[n]) ;
467 std::complex<double> gan = (ki * Zi / k) * (dHnka.real() * Hnka[n] - dHnka * Hnka[n].real() ) /
468 (Zi*dHnkia.real() * Hnka[n] - dHnka * Hnkia[n].real()) ;
469
470 std::complex<double> dn = gan*an;
471
472 std::complex<double> jnkir = Hnkir[n].real() / kir;
473
474 double Pn1 = Legendre(n, 1, ctheta);
475
476 Er += (n * (n + 1.) * jnkir * dn ) * Pn1;
477 Etheta += an * (I * gan * dHnkir.real() * A1 + den * Hnkir[n].real() * A1 );
478 Ephi += an * (I * gan * dHnkir.real() * A1 + den * Hnkir[n].real() * A1 );
479 }
480
481 Er *= I * cos(phi) / kir;
482 Etheta *= (1. / kir) * (cos(phi));
483 Ephi *= (-1. / kir) * (sin(phi));
484
485 }
486 }
487
488 else if( theta == M_PI ) {
489
490 if (r >= a ) {
491 for (int n = 1 ; n < ns ; n++){
492 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
493
494 double A2 = -pow(-1, n + 1) * n * (n + 1.) / 2.;
495 double A3 = -pow(-1, n ) * n * (n + 1.) / 2.;
496
497 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
498 std::complex<double> dHnkia = - Hnkia[n + 1] + (n + 1.) * Hnkia[n] / kia;
499
500 std::complex<double> aln = (dHnka.real() * Hnkia[n].real() - Zi * dHnkia.real() * Hnka[n].real() ) /
501 (Zi*dHnkia.real() * Hnka[n] - dHnka * Hnkia[n].real()) ;
502 std::complex<double> ben = (dHnkia.real() * Hnka[n].real() - Zi * dHnka.real() * Hnkia[n].real() ) /
503 (Zi*Hnkia[n].real() * dHnka - dHnkia.real() * Hnka[n]) ;
504
505 std::complex<double> dn = aln*an;
506
507 std::complex<double> dHnkr = - Hnkr[n + 1] + (n + 1.) * Hnkr[n] / kr;
508 std::complex<double> d2Hnkr = Hnkr[n] / kr;
509 std::complex<double> jnkr = Hnkr[n].real() / kr;
510
511 double Pn1 = Legendre(n, 1, ctheta);
512
513 Er += (n * (n + 1.) * d2Hnkr * dn + n * (n + 1.) * jnkr * an ) * Pn1;
514 Etheta += an * (I * aln * dHnkr * A3 + ben * Hnkr[n] * A2 + I * dHnkr.real() * A3 + Hnkr[n].real() * A2 );
515 Ephi += an * (I * aln * dHnkr * A2 + ben * Hnkr[n] * A3 + I * dHnkr.real() * A2 + Hnkr[n].real() * A3 );
516 }
517
518 Er *= I * cos(phi) / kr;
519 Etheta *= (1. / kr) * (cos(phi));
520 Ephi *= (-1. / kr) * (sin(phi));
521 }
522 else {
523 for (int n = 1 ; n < ns ; n++){
524 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
525
526 double A2 = -pow(-1, n + 1) * n * (n + 1.) / 2.;
527 double A3 = -pow(-1, n ) * n * (n + 1.) / 2.;
528
529 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
530 std::complex<double> dHnkia = - Hnkia[n + 1] + (n + 1.) * Hnkia[n] / kia;
531 std::complex<double> dHnkir = - Hnkir[n + 1] + (n + 1.) * Hnkir[n] / kir;
532
533 std::complex<double> den = (ki * Zi / k) * (dHnka * Hnka[n].real() - dHnka.real() * Hnka[n] ) /
534 (Zi*Hnkia[n].real() * dHnka - dHnkia.real() * Hnka[n]) ;
535 std::complex<double> gan = (ki * Zi / k) * (dHnka.real() * Hnka[n] - dHnka * Hnka[n].real() ) /
536 (Zi*dHnkia.real() * Hnka[n] - dHnka * Hnkia[n].real()) ;
537
538 std::complex<double> dn = gan*an;
539
540 std::complex<double> jnkir = Hnkir[n].real() / kir;
541
542 double Pn1 = Legendre(n, 1, ctheta);
543
544 Er += (n * (n + 1.) * jnkir * dn ) * Pn1;
545 Etheta += an * (I * gan * dHnkir.real() * A3 + den * Hnkir[n].real() * A2 );
546 Ephi += an * (I * gan * dHnkir.real() * A2 + den * Hnkir[n].real() * A3 );
547 }
548
549 Er *= I * cos(phi) / kir;
550 Etheta *= (1. / kir) * (cos(phi));
551 Ephi *= (-1. / kir) * (sin(phi));
552
553 }
554 }
555
556 else {
557 if (r >= a ) {
558 for (int n = 1 ; n < ns ; n++){
559 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
560
561
562 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
563 std::complex<double> dHnkia = - Hnkia[n + 1] + (n + 1.) * Hnkia[n] / kia;
564
565 std::complex<double> aln = (dHnka.real() * Hnkia[n].real() - Zi * dHnkia.real() * Hnka[n].real() ) /
566 (Zi*dHnkia.real() * Hnka[n] - dHnka * Hnkia[n].real()) ;
567 std::complex<double> ben = (dHnkia.real() * Hnka[n].real() - Zi * dHnka.real() * Hnkia[n].real() ) /
568 (Zi*Hnkia[n].real() * dHnka - dHnkia.real() * Hnka[n]) ;
569
570 std::complex<double> dn = aln*an;
571
572 std::complex<double> dHnkr = - Hnkr[n + 1] + (n + 1.) * Hnkr[n] / kr;
573 std::complex<double> d2Hnkr = Hnkr[n] / kr;
574 std::complex<double> jnkr = Hnkr[n].real() / kr;
575
576 double Pn1 = Legendre(n, 1, ctheta);
577 double Pn11 = Legendre(n+1, 1, ctheta);
578 double dPn1 = n * Pn11 - (n + 1) * ctheta * Pn1 ;
579
580
581 Er += (n * (n + 1.) * d2Hnkr * dn + n * (n + 1.) * jnkr * an ) * Pn1;
582 Etheta += an * (I * aln * dHnkr * dPn1 + ben * Hnkr[n] * Pn1 + I * dHnkr.real() * dPn1 + Hnkr[n].real() * Pn1 );
583 Ephi += an * (I * aln * dHnkr * Pn1 + ben * Hnkr[n] * dPn1 + I * dHnkr.real() * Pn1 + Hnkr[n].real() * dPn1 );
584 }
585
586 Er *= I * cos(phi) / kr;
587 Etheta *= (1. / kr) * (cos(phi)/stheta);
588 Ephi *= (-1. / kr) * (sin(phi)/stheta);
589
590 }
591 else {
592 for (int n = 1 ; n < ns ; n++){
593 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
594
595 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
596 std::complex<double> dHnkia = - Hnkia[n + 1] + (n + 1.) * Hnkia[n] / kia;
597 std::complex<double> dHnkir = - Hnkir[n + 1] + (n + 1.) * Hnkir[n] / kir;
598
599 std::complex<double> den = (ki * Zi / k) * (dHnka * Hnka[n].real() - dHnka.real() * Hnka[n] ) /
600 (Zi*Hnkia[n].real() * dHnka - dHnkia.real() * Hnka[n]) ;
601 std::complex<double> gan = (ki * Zi / k) * (dHnka.real() * Hnka[n] - dHnka * Hnka[n].real() ) /
602 (Zi*dHnkia.real() * Hnka[n] - dHnka * Hnkia[n].real()) ;
603
604 std::complex<double> dn = gan*an;
605
606 std::complex<double> jnkir = Hnkir[n].real() / kir;
607
608 double Pn1 = Legendre(n, 1, ctheta);
609 double Pn11 = Legendre(n+1, 1, ctheta);
610 double dPn1 = n * Pn11 - (n + 1) * ctheta * Pn1 ;
611
612 Er += (n * (n + 1.) * jnkir * dn ) * Pn1;
613 Etheta += an * (I * gan * dHnkir.real() * dPn1 + den * Hnkir[n].real() * Pn1 );
614 Ephi += an * (I * gan * dHnkir.real() * Pn1 + den * Hnkir[n].real() * dPn1 );
615 }
616
617 Er *= I * cos(phi) / kir;
618 Etheta *= (1. / kir) * (cos(phi)/stheta);
619 Ephi *= (-1. / kir) * (sin(phi)/stheta);
620
621 }
622 }
623
624
625
626 // r, theta, phi components
627 std::complex<double> rtp[3] = {Er, Etheta, Ephi};
628
629 double mat[3][3];
630 // r basis vector
631 mat[0][0] = sin(theta) * cos(phi) ;
632 mat[0][1] = sin(theta) * sin(phi) ;
633 mat[0][2] = cos(theta) ;
634 // theta basis vector
635 mat[1][0] = cos(theta) * cos(phi) ;
636 mat[1][1] = cos(theta) * sin(phi) ;
637 mat[1][2] = - sin(theta) ;
638 // phi basis vector
639 mat[2][0] = - sin(phi) ;
640 mat[2][1] = cos(phi);
641 mat[2][2] = 0.;
642
643 // x, y, z components
644 std::complex<double> xyz[3] = {0., 0., 0.};
645 for(int i = 0; i < 3; i++)
646 for(int j = 0; j < 3; j++)
647 xyz[i] = xyz[i] + mat[j][i] * rtp[j];
648
649 V->Val[0] = xyz[0].real();
650 V->Val[1] = xyz[1].real();
651 V->Val[2] = xyz[2].real();
652 V->Val[MAX_DIM] = xyz[0].imag();
653 V->Val[MAX_DIM+1] = xyz[1].imag();
654 V->Val[MAX_DIM+2] = xyz[2].imag();
655
656 V->Type = VECTOR ;
657 }
658
659 // diffraction onde plane par une sphere PEC en -iwt
F_ElectricFieldPerfectlyConductingSphereMwt(F_ARG)660 void F_ElectricFieldPerfectlyConductingSphereMwt(F_ARG)
661 {
662 double x = A->Val[0];
663 double y = A->Val[1];
664 double z = A->Val[2];
665 double r = sqrt(x * x + y * y + z * z);
666 double theta = atan2(sqrt(x * x + y * y), z);
667 double phi = atan2(y, x);
668
669
670 double k = Fct->Para[0] ;
671 double a = Fct->Para[1] ;
672 double kr = k * r;
673 double ka = k * a;
674
675 int ns = (int)k + 12;
676
677 std::vector<std::complex<double> > Hnkr(ns + 1), Hnka(ns + 1);
678 for (int n = 1 ; n <= ns ; n++){
679 Hnkr[n] = std::complex<double>(AltSpherical_j_n(n, kr), AltSpherical_y_n(n, kr));
680 Hnka[n] = std::complex<double>(AltSpherical_j_n(n, ka), AltSpherical_y_n(n, ka));
681
682 }
683 double ctheta = cos(theta);
684 double stheta = sin(theta);
685
686 std::complex<double> Er(0.,0), Etheta(0.,0), Ephi(0.,0), I(0, 1.);
687
688 if( theta == 0.) {
689 for (int n = 1 ; n < ns ; n++){
690 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
691
692 double A1 = n * (n + 1.) / 2.;
693 // PS: the following is correct (Hn(z) = z hn(z) is not a standard spherical
694 // bessel function!)
695 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
696 std::complex<double> bn = - dHnka.real() / dHnka;
697 std::complex<double> cn = - Hnka[n].real() / Hnka[n];
698 std::complex<double> dn = bn*an;
699
700 std::complex<double> dHnkr = - Hnkr[n + 1] + (n + 1.) * Hnkr[n] / kr;
701 std::complex<double> d2Hnkr = Hnkr[n] / kr;
702
703 double Pn1 = Legendre(n, 1, ctheta);
704 Er += n * (n + 1.) * d2Hnkr * Pn1 * dn;
705 Etheta += an * (I * bn * dHnkr * A1 + cn * Hnkr[n] * A1);
706 Ephi += an * (I * bn * dHnkr * A1 + cn * Hnkr[n] * A1);
707 }
708
709 Er *= I * cos(phi) / kr;
710 Etheta *= (1. / kr) * (cos(phi));
711 Ephi *= (-1. / kr) * (sin(phi));
712 }
713 else if( theta == M_PI ) {
714 for (int n = 1 ; n < ns ; n++){
715 std::complex<double> an = std::pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
716
717 double A2 = pow(-1, n + 1) * n * (n + 1.) / 2.;
718 double A3 = pow(-1, n ) * n * (n + 1.) / 2.;
719
720 // PS: the following is correct (Hn(z) = z hn(z) is not a standard spherical
721 // bessel function!)
722 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
723 std::complex<double> bn = - dHnka.real() / dHnka;
724 std::complex<double> cn = - Hnka[n].real() / Hnka[n];
725 std::complex<double> dn = bn*an;
726
727 std::complex<double> dHnkr = - Hnkr[n + 1] + (n + 1.) * Hnkr[n] / kr;
728 std::complex<double> d2Hnkr = Hnkr[n] / kr;
729
730 double Pn1 = Legendre(n, 1, ctheta);
731 Er += n * (n + 1.) * d2Hnkr * Pn1 * dn;
732 Etheta += an * (I * bn * dHnkr * A3 + cn * Hnkr[n] * A2 );
733 Ephi += an * (I * bn * dHnkr * A2 + cn * Hnkr[n] * A3);
734 }
735
736 Er *= I * cos(phi) / kr;
737 Etheta *= (1.0 / kr) * cos(phi);
738 Ephi *= (-1.0 / kr) * sin(phi);
739 }
740 else {
741
742 for (int n = 1 ; n < ns ; n++){
743 std::complex<double> an = std::pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
744
745 // PS: the following is correct (Hn(z) = z hn(z) is not a standard spherical
746 // bessel function!)
747 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
748 std::complex<double> bn = - dHnka.real() / dHnka;
749 std::complex<double> cn = - Hnka[n].real() / Hnka[n];
750 std::complex<double> dn = bn * an;
751
752 std::complex<double> dHnkr = - Hnkr[n + 1] + (n + 1.) * Hnkr[n] / kr;
753 std::complex<double> d2Hnkr = Hnkr[n] / kr;
754
755 double Pn1 = Legendre(n, 1, ctheta);
756 double Pn11 = Legendre(n+1, 1, ctheta);
757 double dPn1 = n * Pn11 - (n + 1) * ctheta * Pn1 ;
758 Er += n * (n + 1.) * d2Hnkr * Pn1 * dn;
759 Etheta += an * (I * bn * dHnkr * dPn1 + cn * Hnkr[n] * Pn1 );
760 Ephi += an * (I * bn * dHnkr * Pn1 + cn * Hnkr[n] * dPn1);
761 }
762
763 Er *= I * cos(phi) / kr;
764 Etheta *= (1. / kr) * (cos(phi)/stheta);
765 Ephi *= (-1. / kr) * (sin(phi)/stheta);
766 }
767
768 // r, theta, phi components
769 std::complex<double> rtp[3] = {Er, Etheta, Ephi};
770
771 double mat[3][3];
772 // r basis vector
773 mat[0][0] = sin(theta) * cos(phi) ;
774 mat[0][1] = sin(theta) * sin(phi) ;
775 mat[0][2] = cos(theta) ;
776 // theta basis vector
777 mat[1][0] = cos(theta) * cos(phi) ;
778 mat[1][1] = cos(theta) * sin(phi) ;
779 mat[1][2] = - sin(theta) ;
780 // phi basis vector
781 mat[2][0] = - sin(phi) ;
782 mat[2][1] = cos(phi);
783 mat[2][2] = 0.;
784
785 // x, y, z components
786 std::complex<double> xyz[3] = {0., 0., 0.};
787 for(int i = 0; i < 3; i++)
788 for(int j = 0; j < 3; j++)
789 xyz[i] = xyz[i] + mat[j][i] * rtp[j];
790
791 V->Val[0] = xyz[0].real();
792 V->Val[1] = xyz[1].real();
793 V->Val[2] = xyz[2].real();
794 V->Val[MAX_DIM] = xyz[0].imag();
795 V->Val[MAX_DIM+1] = xyz[1].imag();
796 V->Val[MAX_DIM+2] = xyz[2].imag();
797
798 V->Type = VECTOR ;
799 }
800
801 // calcul la solution exact de OSRC sur la sphere
F_ExactOsrcSolutionPerfectlyConductingSphereMwt(F_ARG)802 void F_ExactOsrcSolutionPerfectlyConductingSphereMwt(F_ARG)
803 {
804 double x = A->Val[0];
805 double y = A->Val[1];
806 double z = A->Val[2];
807 double theta = atan2(sqrt(x * x + y * y), z);
808 double phi = atan2(y, x);
809
810 double k = Fct->Para[0] ;
811 double a = Fct->Para[1] ;
812
813 double ka = k * a;
814
815 int ns = (int)k + 10;
816
817 std::vector<std::complex<double> > Hnka(ns + 1);
818 for (int n = 1 ; n <= ns ; n++){
819 Hnka[n] = std::complex<double>(AltSpherical_j_n(n, ka), AltSpherical_y_n(n, ka));
820
821 }
822 double ctheta = cos(theta);
823 double stheta = sin(theta);
824
825 std::complex<double> Er(0.,0), Etheta(0.,0), Ephi(0.,0), I(0, 1.);
826
827 if( theta == 0.) {
828 for (int n = 1 ; n < ns ; n++){
829 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
830
831 double A1 = n * (n + 1.0) / 2.;
832 double mu_n = 1 - n * (n + 1.0) / (k * k);
833
834 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
835 std::complex<double> bn = dHnka.real() ;
836 std::complex<double> cn = Hnka[n].real();
837
838
839 if ( k * k >= n * (n + 1) ) {
840 Etheta += an * (-I * cn * A1 * sqrt(mu_n) + bn * A1 / sqrt(mu_n));
841 Ephi += an * ( I * cn * A1 * sqrt(mu_n) - bn * A1 / sqrt(mu_n)); }
842 else{
843 Etheta += an * (-I * cn * A1 * I * sqrt(-mu_n) - I * bn * A1 / sqrt(-mu_n));
844 Ephi += an * ( I * cn * A1 * I * sqrt(-mu_n) + I * bn * A1 / sqrt(-mu_n));
845 }
846 }
847
848 Etheta *= cos(phi);
849 Ephi *= sin(phi);
850 }
851 else if( theta == M_PI) {
852 for (int n = 1 ; n < ns ; n++){
853 std::complex<double> an = std::pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
854 double A2 = pow(-1.0, n + 1.0) * n * (n + 1.) / 2.;
855 double A3 = pow(-1.0, n + 0.0) * n * (n + 1.) / 2.;
856 double mu_n = 1 - n * (n + 1.0) / (k * k);
857
858 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
859 std::complex<double> bn = dHnka.real();
860 std::complex<double> cn = Hnka[n].real();
861
862 if ( k * k >= n * (n+1)) {
863 Etheta += an * (-I * cn * A2 * sqrt(mu_n) + bn * A3 / sqrt(mu_n));
864 Ephi += an * ( I * cn * A3 * sqrt(mu_n) - bn * A2 / sqrt(mu_n));}
865 else {
866 Etheta += an * (-I * cn * A2 * I * sqrt(-mu_n) - I * bn * A3 / sqrt(-mu_n));
867 Ephi += an * ( I * cn * A3 * I * sqrt(-mu_n) + I * bn * A2 / sqrt(-mu_n));
868 }
869
870 }
871
872 Etheta *= cos(phi);
873 Ephi *= sin(phi);
874 }
875 else{
876 for (int n = 1 ; n < ns ; n++){
877 std::complex<double> an = std::pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
878
879 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
880 std::complex<double> bn = dHnka.real();
881 std::complex<double> cn = Hnka[n].real();
882
883 double mu_n = 1 - n * (n + 1.0) / (k * k);
884 double Pn1 = Legendre(n, 1, ctheta);
885 double Pn11 = Legendre(n+1, 1, ctheta);
886 double dPn1 = n * Pn11 - (n + 1) * ctheta * Pn1 ;
887
888 if ( k * k >= n * (n+1)) {
889 Etheta += an * (-I * cn * Pn1 * sqrt(mu_n) + bn * dPn1 / sqrt(mu_n));
890 Ephi += an * ( I * cn * dPn1 * sqrt(mu_n) - bn * Pn1 / sqrt(mu_n));}
891 else {
892 Etheta += an * (-I * cn * Pn1 * I * sqrt(-mu_n) -I * bn * dPn1 / sqrt(-mu_n));
893 Ephi += an * ( I * cn * dPn1 * I * sqrt(-mu_n) + I * bn * Pn1 / sqrt(-mu_n));
894 }
895
896 }
897
898 Etheta *= cos(phi)/stheta;
899 Ephi *= sin(phi) /stheta;
900 }
901
902 Etheta *= -I / k;
903 Ephi *= -I / k;
904
905 // r, theta, phi components
906 std::complex<double> rtp[3] = {Er, Etheta, Ephi};
907
908 double mat[3][3];
909 // r basis vector
910 mat[0][0] = sin(theta) * cos(phi) ;
911 mat[0][1] = sin(theta) * sin(phi) ;
912 mat[0][2] = cos(theta) ;
913 // theta basis vector
914 mat[1][0] = cos(theta) * cos(phi) ;
915 mat[1][1] = cos(theta) * sin(phi) ;
916 mat[1][2] = - sin(theta) ;
917 // phi basis vector
918 mat[2][0] = - sin(phi) ;
919 mat[2][1] = cos(phi);
920 mat[2][2] = 0.;
921
922 // x, y, z components
923 std::complex<double> xyz[3] = {0., 0., 0.};
924 for(int i = 0; i < 3; i++)
925 for(int j = 0; j < 3; j++)
926 xyz[i] = xyz[i] + mat[j][i] * rtp[j];
927
928 V->Val[0] = xyz[0].real();
929 V->Val[1] = xyz[1].real();
930 V->Val[2] = xyz[2].real();
931 V->Val[MAX_DIM] = xyz[0].imag();
932 V->Val[MAX_DIM+1] = xyz[1].imag();
933 V->Val[MAX_DIM+2] = xyz[2].imag();
934
935 V->Type = VECTOR ;
936 }
937
938 // returne n /\ H en -iwt
F_CurrentPerfectlyConductingSphereMwt(F_ARG)939 void F_CurrentPerfectlyConductingSphereMwt(F_ARG)
940 {
941 double x = A->Val[0];
942 double y = A->Val[1];
943 double z = A->Val[2];
944
945 double theta = atan2(sqrt(x * x + y * y), z);
946 double phi = atan2(y, x);
947
948 double k = Fct->Para[0] ;
949 double a = Fct->Para[1] ;
950 double Z0 = Fct->Para[2] ;
951
952 double ka = k * a;
953
954 int ns = (int)k + 10;
955
956 std::vector<std::complex<double> > Hnka(ns + 1);
957 for (int n = 1 ; n <= ns ; n++){
958 Hnka[n] = std::complex<double>(AltSpherical_j_n(n, ka), AltSpherical_y_n(n, ka));
959
960 }
961 double ctheta = cos(theta);
962 double stheta = sin(theta);
963
964 std::complex<double> Er(0.,0), Etheta(0.,0), Ephi(0.,0), I(0, 1.);
965
966 if( theta == 0.) {
967 for (int n = 1 ; n < ns ; n++){
968 std::complex<double> an = pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
969
970 double A1 = n * (n + 1.0) / 2.;
971
972 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
973 std::complex<double> bn = - dHnka.real() / dHnka;
974 std::complex<double> cn = - Hnka[n].real() / Hnka[n];
975
976
977 Etheta += an * (I * cn * dHnka * A1 + bn * Hnka[n] * A1);
978 Ephi += an * (I * cn * dHnka * A1 + bn * Hnka[n] * A1);
979 }
980
981 Etheta *= (-1. / (Z0*ka)) * (sin(phi));
982 Ephi *= (-1. / (Z0*ka)) * (cos(phi));
983 }
984 else if( theta == M_PI) {
985 for (int n = 1 ; n < ns ; n++){
986 std::complex<double> an = std::pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
987 double A2 = pow(-1.0, n + 1.0) * n * (n + 1.) / 2.;
988 double A3 = pow(-1.0, n + 0.0) * n * (n + 1.) / 2.;
989
990 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
991 std::complex<double> bn = - dHnka.real() / dHnka;
992 std::complex<double> cn = - Hnka[n].real() / Hnka[n];
993
994
995 Etheta += an * (I * cn * dHnka * A3 + bn * Hnka[n] * A2);
996 Ephi += an * (I * cn * dHnka * A2 + bn * Hnka[n] * A3);
997 }
998
999 Etheta *= (-1.0 / (Z0*ka)) * sin(phi);
1000 Ephi *= (-1.0 / (Z0*ka)) * cos(phi);
1001 }
1002 else{
1003 for (int n = 1 ; n < ns ; n++){
1004 std::complex<double> an = std::pow(-I, n) * (2. * n + 1.) / (n * (n + 1.));
1005
1006 std::complex<double> dHnka = - Hnka[n + 1] + (n + 1.) * Hnka[n] / ka;
1007 std::complex<double> bn = - dHnka.real() / dHnka;
1008 std::complex<double> cn = - Hnka[n].real() / Hnka[n];
1009
1010 double Pn1 = Legendre(n, 1, ctheta);
1011 double Pn11 = Legendre(n+1, 1, ctheta);
1012 double dPn1 = n * Pn11 - (n + 1) * ctheta * Pn1 ;
1013
1014 Etheta += an * (I * cn * dHnka * dPn1 + bn * Hnka[n] * Pn1 );
1015 Ephi += an * (I * cn * dHnka * Pn1 + bn * Hnka[n] * dPn1);
1016 }
1017
1018 Etheta *= (-1. / (Z0*ka)) * (sin(phi)/stheta);
1019 Ephi *= (-1. / (Z0*ka)) * (cos(phi) /stheta);
1020 }
1021
1022 // r, theta, phi components
1023 std::complex<double> rtp[3] = {Er, -Ephi, Etheta};
1024
1025 double mat[3][3];
1026 // r basis vector
1027 mat[0][0] = sin(theta) * cos(phi) ;
1028 mat[0][1] = sin(theta) * sin(phi) ;
1029 mat[0][2] = cos(theta) ;
1030 // theta basis vector
1031 mat[1][0] = cos(theta) * cos(phi) ;
1032 mat[1][1] = cos(theta) * sin(phi) ;
1033 mat[1][2] = - sin(theta) ;
1034 // phi basis vector
1035 mat[2][0] = - sin(phi) ;
1036 mat[2][1] = cos(phi);
1037 mat[2][2] = 0.;
1038
1039 // x, y, z components
1040 std::complex<double> xyz[3] = {0., 0., 0.};
1041 for(int i = 0; i < 3; i++)
1042 for(int j = 0; j < 3; j++)
1043 xyz[i] = xyz[i] + mat[j][i] * rtp[j];
1044
1045 V->Val[0] = xyz[0].real();
1046 V->Val[1] = xyz[1].real();
1047 V->Val[2] = xyz[2].real();
1048 V->Val[MAX_DIM] = xyz[0].imag();
1049 V->Val[MAX_DIM+1] = xyz[1].imag();
1050 V->Val[MAX_DIM+2] = xyz[2].imag();
1051
1052 V->Type = VECTOR ;
1053 }
1054
1055 // version avec +iwt
1056 /* Scattering by a perfectly conducting sphere of radius R, under plane wave
1057 incidence pol*e^{ik \alpha\dot\r}, with alpha = (0,0,-1) and pol =
1058 (1,0,0). Returns surface current (From Harrington, Time-harmonic
1059 electromagnetic fields, p. 294) */
1060
F_CurrentPerfectlyConductingSphere(F_ARG)1061 void F_CurrentPerfectlyConductingSphere(F_ARG)
1062 {
1063 cplx I = {0., 1.}, tmp, *hn, coef1, coef2, an, jtheta, jphi, rtp[3], xyz[3];
1064 double k, R, r, kR, theta, phi, Z0, ctheta, stheta, Pn0, Pn1, dPn1, mat[3][3], x, y, z ;
1065 int n, ns, i, j ;
1066
1067 x = A->Val[0];
1068 y = A->Val[1];
1069 z = A->Val[2];
1070 r = sqrt(x*x+y*y+z*z);
1071 theta = atan2(sqrt(x*x+y*y), z);
1072 phi = atan2(y,x);
1073
1074 // warning: approximation
1075 if (theta == 0. ) theta += 1e-6;
1076 if (theta == M_PI || theta == -M_PI) theta -= 1e-6;
1077
1078 k = Fct->Para[0] ;
1079 R = Fct->Para[1] ;
1080 Z0 = Fct->Para[2] ; // impedance of vacuum = sqrt(mu_0/eps_0) \approx 120*pi
1081 kR = k * R;
1082
1083 // test position to check if on sphere
1084 if(fabs(r - R) > 1.e-3) Message::Error("Evaluation point not on sphere");
1085
1086 V->Val[0] = 0.;
1087 V->Val[MAX_DIM] = 0. ;
1088
1089 ns = (int)k + 10;
1090
1091 hn = (cplx*)Malloc((ns + 1)*sizeof(cplx));
1092
1093 for (n = 0 ; n < ns + 1 ; n++){
1094 hn[n].r = AltSpherical_j_n(n, kR);
1095 hn[n].i = - AltSpherical_y_n(n, kR);
1096 }
1097
1098 ctheta = cos(theta);
1099 stheta = sin(theta);
1100
1101 jtheta.r = 0;
1102 jtheta.i = 0;
1103 jphi.r = 0;
1104 jphi.i = 0;
1105
1106 for (n = 1 ; n < ns ; n++){
1107 // 1 / \hat{H}_n^2 (ka)
1108 coef1 = Cdivr( 1.0 , hn[n] );
1109 // 1 / \hat{H}_n^2' (ka)
1110 coef2 = Cdivr( 1.0 , Csub( Cprodr( (double)(n+1) / kR , hn[n]) , hn[n+1] ) );
1111
1112 Pn0 = Legendre(n, 0, ctheta);
1113 Pn1 = Legendre(n, 1, ctheta);
1114
1115 dPn1 = (n+1)*n* Pn0/stheta - (ctheta/(ctheta*ctheta-1))* Pn1;
1116 an = Cprodr( (2.*n+1.) / (double)(n * (n+1.)) , Cpow(I, -n) );
1117
1118 tmp = Cprod( an , Csum( Cprodr( stheta * dPn1 , coef2 ) ,
1119 Cprodr( Pn1 / stheta , Cprod( I , coef1 )) ) );
1120 jtheta = Csum( jtheta, tmp );
1121
1122 tmp = Cprod( an , Csub( Cprodr( Pn1 / stheta , coef2 ) ,
1123 Cprodr( dPn1 * stheta , Cdiv(coef1 , I)) ) );
1124 jphi = Csum( jphi, tmp );
1125 }
1126
1127 Free(hn);
1128
1129 tmp = Cprodr( cos(phi)/(kR*Z0), I);
1130 jtheta = Cprod( jtheta, tmp );
1131
1132 tmp = Cprodr( sin(phi)/(kR*Z0), I );
1133 jphi = Cprod( jphi, tmp );
1134
1135 // r, theta, phi components
1136 rtp[0].r = 0;
1137 rtp[0].i = 0;
1138 rtp[1] = jtheta;
1139 rtp[2] = jphi;
1140
1141 // r basis vector
1142 mat[0][0] = sin(theta) * cos(phi) ;
1143 mat[0][1] = sin(theta) * sin(phi) ;
1144 mat[0][2] = cos(theta) ;
1145 // theta basis vector
1146 mat[1][0] = cos(theta) * cos(phi) ;
1147 mat[1][1] = cos(theta) * sin(phi) ;
1148 mat[1][2] = - sin(theta) ;
1149 // phi basis vector
1150 mat[2][0] = - sin(phi) ;
1151 mat[2][1] = cos(phi);
1152 mat[2][2] = 0.;
1153
1154 // x, y, z components
1155 for(i = 0; i < 3; i++){
1156 xyz[i].r = 0;
1157 xyz[i].i = 0;
1158 for(j = 0; j < 3; j++)
1159 xyz[i] = Csum( xyz[i] , Cprodr(mat[j][i] , rtp[j]) );
1160 }
1161
1162 V->Val[0] = xyz[0].r;
1163 V->Val[1] = xyz[1].r;
1164 V->Val[2] = xyz[2].r;
1165 V->Val[MAX_DIM] = xyz[0].i;
1166 V->Val[MAX_DIM+1] = xyz[1].i;
1167 V->Val[MAX_DIM+2] = xyz[2].i;
1168
1169 V->Type = VECTOR ;
1170 }
1171
1172 /* Scattering by an acoustically soft sphere (exterior Dirichlet problem) of
1173 radius R, under plane wave incidence e^{ikx}. Returns scattered field
1174 outside. (Colton and Kress, Inverse Acoustic..., p 51, eq. 3.29) */
1175
F_AcousticFieldSoftSphere(F_ARG)1176 void F_AcousticFieldSoftSphere(F_ARG)
1177 {
1178 double x = A->Val[0];
1179 double y = A->Val[1];
1180 double z = A->Val[2];
1181 double r = sqrt(x*x + y*y + z*z);
1182 double k = Fct->Para[0];
1183 double R = Fct->Para[1];
1184 double kr = k*r;
1185 double kR = k*R;
1186
1187 // 3rd/4th/5th parameters: incidence direction
1188 double XdotK;
1189 if(Fct->NbrParameters > 4){
1190 double dx = Fct->Para[2];
1191 double dy = Fct->Para[3];
1192 double dz = Fct->Para[4];
1193 double dr = sqrt(dx*dx + dy*dy + dz*dz);
1194 XdotK = (x*dx + y*dy + z*dz)/(r*dr);
1195 }
1196 else{
1197 XdotK = x/r;
1198 }
1199 XdotK = (XdotK > 1) ? 1 : XdotK;
1200 XdotK = (XdotK < -1) ? -1 : XdotK;
1201
1202 // 6th/7th parameters: range of modes
1203 int nStart = 0;
1204 int nEnd = (int)(kR) + 10;
1205 if(Fct->NbrParameters > 5){
1206 nStart = Fct->Para[5];
1207 nEnd = (Fct->NbrParameters > 6) ? Fct->Para[6] : nStart+1;
1208 }
1209
1210 std::complex<double> I(0,1);
1211 double vr=0, vi=0;
1212 #if defined(_OPENMP)
1213 #pragma omp parallel for reduction(+: vr,vi)
1214 #endif
1215 for (int n=nStart ; n<nEnd; n++){
1216 std::complex<double> hnkR( Spherical_j_n(n,kR), Spherical_y_n(n,kR) );
1217 std::complex<double> hnkr( Spherical_j_n(n,kr), Spherical_y_n(n,kr) );
1218 std::complex<double> tmp1 = std::pow(I,n) * hnkr / hnkR;
1219 double tmp2 = -(2*n+1) * std::real(hnkR) * Legendre(n, 0, XdotK);
1220 vr += tmp2 * std::real(tmp1);
1221 vi += tmp2 * std::imag(tmp1);
1222 }
1223 V->Val[0] = vr;
1224 V->Val[MAX_DIM] = vi;
1225 V->Type = SCALAR ;
1226 }
1227
Dhn_Spherical(cplx * hntab,int n,double x)1228 cplx Dhn_Spherical(cplx *hntab, int n, double x)
1229 {
1230 return Csub( Cprodr( (double)n/x, hntab[n] ) , hntab[n+1] );
1231 }
1232
1233 /* Scattering by acoustically soft circular sphere of radius R0,
1234 under plane wave incidence e^{ikx}, with artificial boundary
1235 condition at R1. Returns exact solution of the (interior!) problem
1236 between R0 and R1. */
1237
F_AcousticFieldSoftSphereABC(F_ARG)1238 void F_AcousticFieldSoftSphereABC(F_ARG)
1239 {
1240 cplx I = {0.,1.}, tmp, alpha1, alpha2, delta, am, bm, lambda;
1241 cplx h1nkR0, *h1nkR1tab, *h2nkR1tab, h1nkr;
1242
1243 double k, R0, R1, r, kr, kR0, kR1, theta, cosfact, sinfact, fact;
1244 int n, ns, ABCtype, SingleMode ;
1245
1246 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1] + A->Val[2]*A->Val[2]) ;
1247 theta = acos(A->Val[0] / r); // angle between position vector and (1,0,0)
1248
1249 k = Fct->Para[0] ;
1250 R0 = Fct->Para[1] ;
1251 R1 = Fct->Para[2] ;
1252 ABCtype = (int)Fct->Para[3] ;
1253 SingleMode = (int)Fct->Para[4] ;
1254 kr = k * r;
1255 kR0 = k * R0;
1256 kR1 = k * R1;
1257
1258 if(ABCtype == 1){ /* Sommerfeld */
1259 lambda = Cprodr(-k, I);
1260 }
1261 else{
1262 Message::Error("ABC type not yet implemented");
1263 }
1264
1265 V->Val[0] = 0.;
1266 V->Val[MAX_DIM] = 0. ;
1267
1268 ns = (int)k + 11;
1269
1270 h1nkR1tab = (cplx*)Malloc(ns * sizeof(cplx));
1271 for (n = 0 ; n < ns ; n++){
1272 h1nkR1tab[n].r = Spherical_j_n(n, kR1);
1273 h1nkR1tab[n].i = Spherical_y_n(n, kR1);
1274 }
1275
1276 h2nkR1tab = (cplx*)Malloc(ns * sizeof(cplx));
1277 for (n = 0 ; n < ns ; n++){
1278 h2nkR1tab[n] = Cconj(h1nkR1tab[n]);
1279 }
1280
1281 for (n = 0 ; n < ns-1 ; n++){
1282 if(SingleMode >= 0 && SingleMode != n) continue;
1283
1284 h1nkR0.r = Spherical_j_n(n, kR0);
1285 h1nkR0.i = Spherical_y_n(n, kR0);
1286
1287 h1nkr.r = Spherical_j_n(n,kr);
1288 h1nkr.i = Spherical_y_n(n,kr);
1289
1290 alpha1 = Csum( Cprodr(k, Dhn_Spherical(h1nkR1tab, n, kR1)) ,
1291 Cprod(lambda, h1nkR1tab[n]) );
1292 alpha2 = Csum( Cprodr(k, Dhn_Spherical(h2nkR1tab, n, kR1)) ,
1293 Cprod(lambda, h2nkR1tab[n]) );
1294 delta = Csub( Cprod( alpha1 , Cconj(h1nkR0) ) ,
1295 Cprod( alpha2 , h1nkR0 ) );
1296
1297 if(Cmodu(delta) < 1.e-6) break;
1298
1299 am = Cdiv( Cprodr(h1nkR0.r, alpha2) ,
1300 delta );
1301 bm = Cdiv( Cprodr(-h1nkR0.r, alpha1) ,
1302 delta );
1303
1304 if(SingleMode >= 0 && SingleMode == n){
1305 tmp = Csum( Cprod( am , h1nkr ) , Cprod( bm , Cconj(h1nkr) ) );
1306 cosfact = (2*n+1) * Legendre(n, 0, cos(theta));
1307 sinfact = (2*n+1) * Legendre(n, 0, sin(theta));
1308 V->Val[0] += cosfact * tmp.r - sinfact * tmp.i;
1309 V->Val[MAX_DIM] += cosfact * tmp.i + sinfact * tmp.r;
1310 }
1311 else{
1312 tmp = Cprod( Cpow(I,n) , Csum( Cprod( am , h1nkr ) ,
1313 Cprod( bm , Cconj(h1nkr) ) ) );
1314 fact = (2*n+1) * Legendre(n, 0, cos(theta));
1315 V->Val[0] += fact * tmp.r;
1316 V->Val[MAX_DIM] += fact * tmp.i;
1317 }
1318 }
1319
1320 Free(h1nkR1tab);
1321 Free(h2nkR1tab);
1322
1323 if(SingleMode < 0){
1324 V->Val[0] *= 1;
1325 V->Val[MAX_DIM] *= 1;
1326 }
1327
1328 V->Type = SCALAR ;
1329 }
1330
1331 /* Scattering by an acoustically soft sphere (exterior Dirichlet problem) of
1332 radius R, under plane wave incidence e^{ikx}. Returns radial derivative of
1333 scattered field outside */
1334
F_DrAcousticFieldSoftSphere(F_ARG)1335 void F_DrAcousticFieldSoftSphere(F_ARG)
1336 {
1337 cplx I = {0.,1.}, hnkR, tmp, *hnkrtab;
1338 double k, R, r, kr, kR, theta, fact ;
1339 int n, ns ;
1340
1341 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1] + A->Val[2]*A->Val[2]) ;
1342 theta = acos(A->Val[0] / r); // angle between position vector and (1,0,0)
1343
1344 k = Fct->Para[0] ;
1345 R = Fct->Para[1] ;
1346 kr = k*r;
1347 kR = k*R;
1348
1349 V->Val[0] = 0.;
1350 V->Val[MAX_DIM] = 0. ;
1351
1352 ns = (int)k + 10;
1353
1354 hnkrtab = (cplx*)Malloc((ns + 1)*sizeof(cplx));
1355
1356 for (n = 0 ; n < ns + 1 ; n++){
1357 hnkrtab[n].r = Spherical_j_n(n, kr);
1358 hnkrtab[n].i = Spherical_y_n(n, kr);
1359 }
1360
1361 for (n = 0 ; n < ns ; n++){
1362 hnkR.r = Spherical_j_n(n, kR);
1363 hnkR.i = Spherical_y_n(n, kR);
1364
1365 tmp = Cdiv( Cprod( Cpow(I,n) , Cprodr(hnkR.r * k, Dhn_Spherical(hnkrtab, n, kr)) ) ,
1366 hnkR );
1367
1368 fact = (2*n+1) * Legendre(n, 0, cos(theta));
1369
1370 V->Val[0] += fact * tmp.r;
1371 V->Val[MAX_DIM] += fact * tmp.i;
1372 }
1373
1374 Free(hnkrtab);
1375
1376 V->Val[0] *= -1;
1377 V->Val[MAX_DIM] *= -1;
1378
1379 V->Type = SCALAR ;
1380 }
1381
1382 /* Scattering by an acoustically soft sphere (exterior Dirichlet problem) of
1383 radius R, under plane wave incidence e^{ikx}. Returns RCS. (Colton and
1384 Kress, Inverse Acoustic..., p 52, eq. 3.30) */
1385
F_RCSSoftSphere(F_ARG)1386 void F_RCSSoftSphere(F_ARG)
1387 {
1388 cplx I = {0.,1.}, hnkR, tmp, res;
1389 double k, R, r, kR, theta, fact, val ;
1390 int n, ns ;
1391
1392 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1] + A->Val[2]*A->Val[2]) ;
1393 theta = acos(A->Val[0] / r); // angle between position vector and (1,0,0)
1394
1395 k = Fct->Para[0] ;
1396 R = Fct->Para[1] ;
1397 kR = k*R;
1398
1399 res.r = 0.;
1400 res.i = 0. ;
1401
1402 ns = (int)k + 10;
1403
1404 for (n = 0 ; n < ns ; n++){
1405 hnkR.r = Spherical_j_n(n, kR);
1406 hnkR.i = Spherical_y_n(n, kR);
1407
1408 tmp = Cdivr( hnkR.r, hnkR );
1409
1410 fact = (2*n+1) * Legendre(n, 0, cos(theta));
1411
1412 res.r += fact * tmp.r;
1413 res.i += fact * tmp.i;
1414 }
1415
1416 val = Cmodu( Cprodr( 1./k , Cprod(res, I) ) );
1417 val *= val;
1418 val *= 4. * M_PI;
1419 val = 10. * log10(val);
1420
1421 V->Val[0] = val;
1422 V->Val[MAX_DIM] = 0.;
1423
1424 V->Type = SCALAR ;
1425 }
1426
1427 /* Scattering by an acoustically hard sphere (exterior Neumann problem) of
1428 radius R, under plane wave incidence e^{ikx}. Returns scattered field
1429 outside */
1430
F_AcousticFieldHardSphere(F_ARG)1431 void F_AcousticFieldHardSphere(F_ARG)
1432 {
1433 double x = A->Val[0];
1434 double y = A->Val[1];
1435 double z = A->Val[2];
1436 double r = sqrt(x*x + y*y + z*z);
1437 double k = Fct->Para[0];
1438 double R = Fct->Para[1];
1439 double kr = k*r;
1440 double kR = k*R;
1441
1442 // 3rd/4th/5th parameters: incidence direction
1443 double XdotK;
1444 if(Fct->NbrParameters > 4){
1445 double dx = Fct->Para[2];
1446 double dy = Fct->Para[3];
1447 double dz = Fct->Para[4];
1448 double dr = sqrt(dx*dx + dy*dy + dz*dz);
1449 XdotK = (x*dx + y*dy + z*dz)/(r*dr);
1450 }
1451 else{
1452 XdotK = x/r;
1453 }
1454 XdotK = (XdotK > 1) ? 1 : XdotK;
1455 XdotK = (XdotK < -1) ? -1 : XdotK;
1456
1457 // 6th/7th parameters: range of modes
1458 int nStart = 0;
1459 int nEnd = (int)(kR) + 10;
1460 if(Fct->NbrParameters > 5){
1461 nStart = Fct->Para[5];
1462 nEnd = (Fct->NbrParameters > 6) ? Fct->Para[6] : nStart+1;
1463 }
1464
1465 std::complex<double> *hnkRtab;
1466 hnkRtab = new std::complex<double>[nEnd+1];
1467 #if defined(_OPENMP)
1468 #pragma omp parallel for
1469 #endif
1470 for (int n=nStart; n<nEnd+1; n++){
1471 hnkRtab[n] = std::complex<double>(Spherical_j_n(n,kR), Spherical_y_n(n,kR));
1472 }
1473
1474 std::complex<double> I(0,1);
1475 double vr=0, vi=0;
1476 #if defined(_OPENMP)
1477 #pragma omp parallel for reduction(+: vr,vi)
1478 #endif
1479 for (int n=nStart ; n<nEnd; n++){
1480 std::complex<double> hnkr( Spherical_j_n(n,kr), Spherical_y_n(n,kr) );
1481 std::complex<double> DhnkR = ((double)n/kR) * hnkRtab[n] - hnkRtab[n+1];
1482 std::complex<double> tmp1 = std::pow(I,n) * hnkr / DhnkR;
1483 double tmp2 = -(2*n+1) * std::real(DhnkR) * Legendre(n, 0, XdotK);
1484 vr += tmp2 * std::real(tmp1);
1485 vi += tmp2 * std::imag(tmp1);
1486 }
1487 V->Val[0] = vr;
1488 V->Val[MAX_DIM] = vi;
1489 V->Type = SCALAR ;
1490
1491 delete hnkRtab;
1492 }
1493
1494 /* Scattering by an acoustically hard sphere (exterior Dirichlet problem) of
1495 radius R, under plane wave incidence e^{ikx}. Returns RCS */
1496
F_RCSHardSphere(F_ARG)1497 void F_RCSHardSphere(F_ARG)
1498 {
1499 cplx I = {0.,1.}, DhnkR, tmp, res, *hnkRtab;
1500 double k, R, r, kR, theta, fact, val ;
1501 int n, ns ;
1502
1503 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1] + A->Val[2]*A->Val[2]) ;
1504 theta = acos(A->Val[0] / r); // angle between position vector and (1,0,0)
1505
1506 k = Fct->Para[0] ;
1507 R = Fct->Para[1] ;
1508 kR = k*R;
1509
1510 res.r = 0.;
1511 res.i = 0. ;
1512
1513 ns = (int)k + 10;
1514
1515 hnkRtab = (cplx*)Malloc((ns + 1)*sizeof(cplx));
1516
1517 for (n = 0 ; n < ns + 1 ; n++){
1518 hnkRtab[n].r = Spherical_j_n(n, kR);
1519 hnkRtab[n].i = Spherical_y_n(n, kR);
1520 }
1521
1522 for (n = 0 ; n < ns ; n++){
1523 DhnkR = Dhn_Spherical(hnkRtab, n, kR);
1524
1525 tmp = Cdivr( DhnkR.r, DhnkR );
1526
1527 fact = (2*n+1) * Legendre(n, 0, cos(theta));
1528
1529 res.r += fact * tmp.r;
1530 res.i += fact * tmp.i;
1531 }
1532
1533 Free(hnkRtab);
1534
1535 val = Cmodu( Cprodr( 1./k , Cprod(res, I) ) );
1536 val *= val;
1537 val *= 4. * M_PI;
1538 val = 10. * log10(val);
1539
1540 V->Val[0] = val;
1541 V->Val[MAX_DIM] = 0.;
1542
1543 V->Type = SCALAR ;
1544 }
1545
1546 /* ------------------------------------------------------------------------ */
1547 /* Exact solutions for cylinders */
1548 /* ------------------------------------------------------------------------ */
1549
1550 /* Scattering by solid PEC cylinder, incident wave z-polarized.
1551 Returns current on cylinder surface */
1552
F_JFIE_ZPolCyl(F_ARG)1553 void F_JFIE_ZPolCyl(F_ARG)
1554 {
1555 double k0, r, kr, e0, eta, phi, a, b, c, d, den ;
1556 int i, ns ;
1557
1558 phi = atan2( A->Val[1], A->Val[0] ) ;
1559
1560 k0 = Fct->Para[0] ;
1561 eta = Fct->Para[1] ;
1562 e0 = Fct->Para[2] ;
1563 r = Fct->Para[3] ;
1564
1565 kr = k0*r ;
1566 ns = 100 ;
1567
1568
1569 V->Val[0] = 0.;
1570 V->Val[MAX_DIM] = 0. ;
1571
1572 for (i = -ns ; i <= ns ; i++ ){
1573 a = cos(i*(phi-(M_PI/2))) ;
1574 b = sin(i*(phi-(M_PI/2))) ;
1575 c = jn(i,kr) ;
1576 d = -yn(i,kr) ;
1577
1578 den = c*c+d*d ;
1579
1580 V->Val[0] += (a*c+b*d)/den ;
1581 V->Val[MAX_DIM] += (b*c-a*d)/den ;
1582 }
1583
1584 V->Val[0] *= -2*e0/kr/eta/M_PI ;
1585 V->Val[MAX_DIM] *= -2*e0/kr/eta/M_PI ;
1586
1587 V->Type = SCALAR ;
1588 }
1589
1590 /* Scattering by solid PEC cylinder, incident wave z-polarized.
1591 Returns RCS */
1592
F_RCS_ZPolCyl(F_ARG)1593 void F_RCS_ZPolCyl(F_ARG)
1594 {
1595 double k0, r, kr, rinf, krinf, phi, a, b, d,den ;
1596 double lambda, bjn, rr = 0., ri = 0. ;
1597 int i, ns ;
1598
1599 phi = atan2( A->Val[1], A->Val[0] ) ;
1600
1601 k0 = Fct->Para[0] ;
1602 r = Fct->Para[1] ;
1603 rinf = Fct->Para[2] ;
1604
1605 kr = k0*r ;
1606 krinf = k0*rinf ;
1607 lambda = 2*M_PI/k0 ;
1608
1609 ns = 100 ;
1610
1611 for (i = -ns ; i <= ns ; i++ ){
1612 bjn = jn(i,kr) ;
1613 a = bjn * cos(i*phi) ;
1614 b = bjn * sin(i*phi) ;
1615 d = -yn(i,kr) ;
1616
1617 den = bjn*bjn+d*d ;
1618
1619 rr += (a*bjn+b*d)/den ;
1620 ri += (b*bjn-a*d)/den ;
1621 }
1622
1623 V->Val[0] = 10*log10( 4*M_PI*SQU(rinf/lambda) * 2/krinf/M_PI *(SQU(rr)+SQU(ri)) ) ;
1624 V->Val[MAX_DIM] = 0. ;
1625
1626 V->Type = SCALAR ;
1627 }
1628
1629 /* Scattering by solid PEC cylinder, incident wave polarized
1630 transversely to z. Returns current on cylinder surface */
1631
F_JFIE_TransZPolCyl(F_ARG)1632 void F_JFIE_TransZPolCyl(F_ARG)
1633 {
1634 double k0, r, kr, h0, phi, a, b, c, d, den ;
1635 int i, ns ;
1636
1637 phi = atan2( A->Val[1], A->Val[0] ) ;
1638
1639 k0 = Fct->Para[0] ;
1640 h0 = Fct->Para[1] ;
1641 r = Fct->Para[2] ;
1642
1643 kr = k0*r ;
1644 ns = 100 ;
1645
1646 V->Val[0] = 0.;
1647 V->Val[MAX_DIM] = 0. ;
1648
1649 for (i = -ns ; i <= ns ; i++ ){
1650 a = cos(M_PI/2 +i*(phi-(M_PI/2))) ;
1651 b = sin(M_PI/2 +i*(phi-(M_PI/2))) ;
1652 c = -jn(i+1,kr) + (i/kr)*jn(i,kr) ;
1653 d = yn(i+1,kr) - (i/kr)*yn(i,kr) ;
1654
1655 den = c*c+d*d ;
1656
1657 V->Val[0] += (a*c+b*d)/den ;
1658 V->Val[MAX_DIM] += (b*c-a*d)/den ;
1659 }
1660
1661 V->Val[0] *= 2*h0/kr/M_PI ;
1662 V->Val[MAX_DIM] *= 2*h0/kr/M_PI ;
1663
1664 V->Type = SCALAR ;
1665 }
1666
1667 /* Scattering by acoustically soft circular cylinder of radius R,
1668 under plane wave incidence e^{ikx}. Returns scatterered field
1669 outside */
1670
F_AcousticFieldSoftCylinder(F_ARG)1671 void F_AcousticFieldSoftCylinder(F_ARG)
1672 {
1673 double theta = atan2(A->Val[1], A->Val[0]);
1674 double r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]);
1675 double k = Fct->Para[0];
1676 double R = Fct->Para[1];
1677 double kr = k*r;
1678 double kR = k*R;
1679
1680 // 3rd parameter: incidence direction
1681 if(Fct->NbrParameters > 2){
1682 double thetaInc = Fct->Para[2];
1683 theta += thetaInc;
1684 }
1685
1686 // 4th/5th parameters: range of modes
1687 int nStart = 0;
1688 int nEnd = (int)(kR) + 10;
1689 if(Fct->NbrParameters > 3){
1690 nStart = Fct->Para[3];
1691 nEnd = (Fct->NbrParameters > 4) ? Fct->Para[4] : nStart+1;
1692 }
1693
1694 std::complex<double> I(0,1);
1695 double vr=0, vi=0;
1696 #if defined(_OPENMP)
1697 #pragma omp parallel for reduction(+: vr,vi)
1698 #endif
1699 for(int n = nStart ; n < nEnd ; n++){
1700 std::complex<double> HnkR( jn(n,kR), yn(n,kR) );
1701 std::complex<double> Hnkr( jn(n,kr), yn(n,kr) );
1702 std::complex<double> tmp1 = std::pow(I,n) * Hnkr/HnkR;
1703 double tmp2 = - (!n ? 1. : 2.) * cos(n*theta) * std::real(HnkR);
1704 std::complex<double> val = tmp1 * tmp2;
1705 vr += std::real(val);
1706 vi += std::imag(val);
1707 }
1708 V->Val[0] = vr;
1709 V->Val[MAX_DIM] = vi;
1710 V->Type = SCALAR;
1711 }
1712
DHn(cplx * Hnkrtab,int n,double x)1713 cplx DHn(cplx *Hnkrtab, int n, double x)
1714 {
1715 if(n == 0){
1716 return Cneg(Hnkrtab[1]);
1717 }
1718 else{
1719 return Csub( Hnkrtab[n-1] , Cprodr((double)n/x, Hnkrtab[n]) );
1720 }
1721 }
1722
1723 /* Scattering by acoustically soft circular cylinder of radius R0,
1724 under plane wave incidence e^{ikx}, with artificial boundary
1725 condition at R1. Returns exact solution of the (interior!) problem
1726 between R0 and R1. */
1727
F_AcousticFieldSoftCylinderABC(F_ARG)1728 void F_AcousticFieldSoftCylinderABC(F_ARG)
1729 {
1730 cplx I = {0.,1.}, tmp, alpha1, alpha2, delta, am, bm, lambda, coef;
1731 cplx H1nkR0, *H1nkR1tab, *H2nkR1tab, H1nkr, alphaBT, betaBT, keps = {0., 0.};
1732
1733 double k, R0, R1, r, kr, kR0, kR1, theta, cost, sint, kappa ;
1734 int n, ns, ABCtype, SingleMode ;
1735
1736 theta = atan2(A->Val[1], A->Val[0]) ;
1737 r = sqrt(A->Val[0] * A->Val[0] + A->Val[1] * A->Val[1]) ;
1738
1739 k = Fct->Para[0] ;
1740 R0 = Fct->Para[1] ;
1741 R1 = Fct->Para[2] ;
1742 ABCtype = (int)Fct->Para[3] ;
1743 SingleMode = (int)Fct->Para[4] ;
1744 kr = k * r;
1745 kR0 = k * R0;
1746 kR1 = k * R1;
1747
1748 if(ABCtype == 1){ /* Sommerfeld */
1749 lambda = Cprodr(-k, I);
1750 }
1751 else if(ABCtype == 2){ /* Bayliss-Turkel */
1752 /*
1753 alphaBT[] = 1/(2*R1) - I[]/(8*k*R1^2*(1+I[]/(k*R1)));
1754 betaBT[] = - 1/(2*I[]*k*(1+I[]/(k*R1)));
1755 */
1756 coef.r = 2*k;
1757 coef.i = 2/R1;
1758 alphaBT = Csubr( 1/(2*R1) , Cdiv(I , Cprodr(4*R1*R1 , coef) ) );
1759 betaBT = Cdiv(I , coef);
1760 }
1761 else if(ABCtype == 3){ /* Pade */
1762 kappa = 1./R1; /* for circular boundary only! */
1763 keps.r = k;
1764 keps.i = 0.4 * pow(k, 1./3.) * pow(kappa, 2./3.);
1765 }
1766 else{
1767 Message::Error("Unknown ABC type");
1768 }
1769
1770 V->Val[0] = 0.;
1771 V->Val[MAX_DIM] = 0. ;
1772
1773 ns = (int)k + 10;
1774
1775 H1nkR1tab = (cplx*)Malloc(ns * sizeof(cplx));
1776 for (n = 0 ; n < ns ; n++){
1777 H1nkR1tab[n].r = jn(n, kR1);
1778 H1nkR1tab[n].i = yn(n, kR1);
1779 }
1780
1781 H2nkR1tab = (cplx*)Malloc(ns * sizeof(cplx));
1782 for (n = 0 ; n < ns ; n++){
1783 H2nkR1tab[n] = Cconj(H1nkR1tab[n]);
1784 }
1785
1786 for (n = 0 ; n < ns ; n++){
1787 if(SingleMode >= 0 && SingleMode != n) continue;
1788
1789 H1nkR0.r = jn(n, kR0);
1790 H1nkR0.i = yn(n, kR0);
1791
1792 H1nkr.r = jn(n,kr);
1793 H1nkr.i = yn(n,kr);
1794
1795 if(ABCtype == 2){
1796 lambda = Csum( Csum( Cprodr(-k, I) , alphaBT ) ,
1797 Cprodr( n*n/(R1*R1) , betaBT ) );
1798 }
1799 else if(ABCtype == 3){
1800 lambda = Cprod( Cprodr(-k, I) ,
1801 Cpow( Csubr(1 , Cdivr(n*n/(R1*R1) , Cprod(keps , keps))) , 0.5));
1802
1803 }
1804
1805 alpha1 = Csum( Cprodr(k, DHn(H1nkR1tab, n, kR1)) ,
1806 Cprod(lambda, H1nkR1tab[n]) );
1807 alpha2 = Csum( Cprodr(k, DHn(H2nkR1tab, n, kR1)) ,
1808 Cprod(lambda, H2nkR1tab[n]) );
1809 delta = Csub( Cprod( alpha1 , Cconj(H1nkR0) ) ,
1810 Cprod( alpha2 , H1nkR0 ) );
1811
1812 if(Cmodu(delta) < 1.e-6) break;
1813
1814 am = Cdiv( Cprodr(H1nkR0.r, alpha2) ,
1815 delta );
1816 bm = Cdiv( Cprodr(-H1nkR0.r, alpha1) ,
1817 delta );
1818
1819 if(SingleMode >= 0 && SingleMode == n){
1820 tmp = Csum( Cprod( am , H1nkr ) , Cprod( bm , Cconj(H1nkr) ) );
1821 cost = cos(n * theta);
1822 sint = sin(n * theta);
1823 V->Val[0] += cost * tmp.r - sint * tmp.i;
1824 V->Val[MAX_DIM] += cost * tmp.i + sint * tmp.r;
1825 }
1826 else{
1827 tmp = Cprod( Cpow(I,n) , Csum( Cprod( am , H1nkr ) ,
1828 Cprod( bm , Cconj(H1nkr) ) ) );
1829 cost = cos(n * theta);
1830 V->Val[0] += cost * tmp.r * (!n ? 0.5 : 1.);
1831 V->Val[MAX_DIM] += cost * tmp.i * (!n ? 0.5 : 1.);
1832 }
1833 }
1834
1835 Free(H1nkR1tab);
1836 Free(H2nkR1tab);
1837
1838 if(SingleMode < 0){
1839 V->Val[0] *= 2;
1840 V->Val[MAX_DIM] *= 2;
1841 }
1842
1843 V->Type = SCALAR ;
1844 }
1845
1846 /* Scattering by acoustically soft circular cylinder of radius R,
1847 under plane wave incidence e^{ikx}. Returns radial derivative of
1848 the solution of the Helmholtz equation outside */
1849
F_DrAcousticFieldSoftCylinder(F_ARG)1850 void F_DrAcousticFieldSoftCylinder(F_ARG)
1851 {
1852 cplx I = {0.,1.}, HnkR, tmp, *Hnkrtab;
1853 double k, R, r, kr, kR, theta, cost ;
1854 int n, ns ;
1855
1856 theta = atan2(A->Val[1], A->Val[0]) ;
1857 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]) ;
1858
1859 k = Fct->Para[0] ;
1860 R = Fct->Para[1] ;
1861 kr = k*r;
1862 kR = k*R;
1863
1864 V->Val[0] = 0.;
1865 V->Val[MAX_DIM] = 0. ;
1866
1867 ns = (int)k + 10;
1868
1869 Hnkrtab = (cplx*)Malloc(ns*sizeof(cplx));
1870
1871 for (n = 0 ; n < ns ; n++){
1872 Hnkrtab[n].r = jn(n,kr);
1873 Hnkrtab[n].i = yn(n,kr);
1874 }
1875
1876 for (n = 0 ; n < ns ; n++){
1877 HnkR.r = jn(n,kR);
1878 HnkR.i = yn(n,kR);
1879
1880 tmp = Cdiv( Cprod( Cpow(I,n) , Cprodr( HnkR.r, DHn(Hnkrtab, n, kr) ) ) , HnkR );
1881
1882 cost = cos(n*theta);
1883
1884 V->Val[0] += cost * tmp.r * (!n ? 0.5 : 1.);
1885 V->Val[MAX_DIM] += cost * tmp.i * (!n ? 0.5 : 1.);
1886 }
1887
1888 Free(Hnkrtab);
1889
1890 V->Val[0] *= -2 * k;
1891 V->Val[MAX_DIM] *= -2 * k;
1892
1893 V->Type = SCALAR ;
1894 }
1895
1896 /* Scattering by acoustically soft circular cylinder of radius R,
1897 under plane wave incidence e^{ikx}. Returns RCS */
1898
F_RCSSoftCylinder(F_ARG)1899 void F_RCSSoftCylinder(F_ARG)
1900 {
1901 cplx I = {0.,1.}, HnkR, Hnkr, res, tmp;
1902 double k, R, r, kR, theta, cost, val ;
1903 int n, ns ;
1904
1905 theta = atan2(A->Val[1], A->Val[0]) ;
1906 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]) ;
1907
1908 k = Fct->Para[0] ;
1909 R = Fct->Para[1] ;
1910 kR = k*R;
1911
1912 res.r = 0.;
1913 res.i = 0. ;
1914
1915 ns = (int)k + 10;
1916
1917 for (n = 0 ; n < ns ; n++){
1918
1919 HnkR.r = jn(n,kR);
1920 HnkR.i = yn(n,kR);
1921
1922 /* leaving r in following asymptotic formula for clarity (see
1923 Colton and Kress, Inverse Acoustic..., p. 65, eq. 3.59) */
1924 Hnkr.r = cos(k*r - n*M_PI/2. - M_PI/4.) / sqrt(k*r) * sqrt(2./M_PI);
1925 Hnkr.i = sin(k*r - n*M_PI/2. - M_PI/4.) / sqrt(k*r) * sqrt(2./M_PI);
1926
1927 tmp = Cdiv( Cprod( Cpow(I,n) , Cprodr( HnkR.r, Hnkr) ) , HnkR );
1928
1929 cost = cos(n*theta);
1930
1931 res.r += cost * tmp.r * (!n ? 0.5 : 1.);
1932 res.i += cost * tmp.i * (!n ? 0.5 : 1.);
1933 }
1934
1935 res.r *= -2;
1936 res.i *= -2;
1937
1938 val = Cmodu(res);
1939 val *= val;
1940 val *= 2. * M_PI * r;
1941 val = 10. * log10(val);
1942
1943 V->Val[0] = val;
1944 V->Val[MAX_DIM] = 0.;
1945
1946 V->Type = SCALAR ;
1947 }
1948
1949 /* Scattering by acoustically hard circular cylinder of radius R,
1950 under plane wave incidence e^{ikx}. Returns scatterered field
1951 outside */
1952
F_AcousticFieldHardCylinder(F_ARG)1953 void F_AcousticFieldHardCylinder(F_ARG)
1954 {
1955 double theta = atan2(A->Val[1], A->Val[0]);
1956 double r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]);
1957 double k = Fct->Para[0];
1958 double R = Fct->Para[1];
1959 double kr = k*r;
1960 double kR = k*R;
1961
1962 // 3rd parameter: incidence direction
1963 if(Fct->NbrParameters > 2){
1964 double thetaInc = Fct->Para[2];
1965 theta += thetaInc;
1966 }
1967
1968 // 4th/5th parameters: range of modes
1969 int nStart = 0;
1970 int nEnd = (int)(kR) + 10;
1971 if(Fct->NbrParameters > 3){
1972 nStart = Fct->Para[3];
1973 nEnd = (Fct->NbrParameters > 4) ? Fct->Para[4] : nStart+1;
1974 }
1975
1976 std::complex<double> *HnkRtab;
1977 HnkRtab = new std::complex<double>[nEnd];
1978 #if defined(_OPENMP)
1979 #pragma omp parallel for
1980 #endif
1981 for (int n=nStart; n<nEnd; n++){
1982 HnkRtab[n] = std::complex<double>(jn(n,kR),yn(n,kR));
1983 }
1984
1985 std::complex<double> I(0,1);
1986 double vr=0, vi=0;
1987 #if defined(_OPENMP)
1988 #pragma omp parallel for reduction(+: vr,vi)
1989 #endif
1990 for (int n=nStart; n<nEnd; n++){
1991 std::complex<double> Hnkr( jn(n,kr), yn(n,kr) );
1992 std::complex<double> dHnkR = (!n ? -HnkRtab[1] : HnkRtab[n-1] - (double)n/kR * HnkRtab[n]);
1993 std::complex<double> tmp1 = std::pow(I,n) * Hnkr/dHnkR;
1994 double tmp2 = - (!n ? 1. : 2.) * cos(n*theta) * std::real(dHnkR);
1995 std::complex<double> val = tmp1 * tmp2;
1996 vr += std::real(val);
1997 vi += std::imag(val);
1998 }
1999
2000 delete HnkRtab;
2001
2002 V->Val[0] = vr;
2003 V->Val[MAX_DIM] = vi;
2004 V->Type = SCALAR ;
2005 }
2006
2007 /* Scattering by acoustically hard circular cylinder of radius R,
2008 under plane wave incidence e^{ikx}. Returns the angular derivative
2009 of the solution outside */
2010
F_DthetaAcousticFieldHardCylinder(F_ARG)2011 void F_DthetaAcousticFieldHardCylinder(F_ARG)
2012 {
2013 cplx I = {0.,1.}, Hnkr, dHnkR, tmp, *HnkRtab;
2014 double k, R, r, kr, kR, theta, sint ;
2015 int n, ns ;
2016
2017 theta = atan2(A->Val[1], A->Val[0]) ;
2018 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]) ;
2019
2020 k = Fct->Para[0] ;
2021 R = Fct->Para[1] ;
2022 kr = k*r;
2023 kR = k*R;
2024
2025 V->Val[0] = 0.;
2026 V->Val[MAX_DIM] = 0. ;
2027
2028 ns = (int)k + 10;
2029
2030 HnkRtab = (cplx*)Malloc(ns*sizeof(cplx));
2031
2032 for (n = 0 ; n < ns ; n++){
2033 HnkRtab[n].r = jn(n,kR);
2034 HnkRtab[n].i = yn(n,kR);
2035 }
2036
2037 for (n = 0 ; n < ns ; n++){
2038 Hnkr.r = jn(n,kr);
2039 Hnkr.i = yn(n,kr);
2040
2041 dHnkR = DHn(HnkRtab, n, kR);
2042
2043 tmp = Cdiv( Cprod( Cpow(I,n) , Cprodr( dHnkR.r, Hnkr) ) , dHnkR );
2044
2045 sint = sin(n*theta);
2046
2047 V->Val[0] += - n * sint * tmp.r * (!n ? 0.5 : 1.);
2048 V->Val[MAX_DIM] += - n * sint * tmp.i * (!n ? 0.5 : 1.);
2049 }
2050
2051 Free(HnkRtab);
2052
2053 V->Val[0] *= -2 ;
2054 V->Val[MAX_DIM] *= -2 ;
2055
2056 V->Type = SCALAR ;
2057 }
2058
2059 /* Scattering by acoustically hard circular cylinder of radius R0,
2060 under plane wave incidence e^{ikx}, with artificial boundary
2061 condition at R1. Returns exact solution of the (interior!) problem
2062 between R0 and R1. */
2063
F_AcousticFieldHardCylinderABC(F_ARG)2064 void F_AcousticFieldHardCylinderABC(F_ARG)
2065 {
2066 cplx I = {0.,1.}, tmp, alpha1, alpha2, delta, am, bm, lambda, coef;
2067 cplx *H1nkR0tab, *H2nkR0tab, *H1nkR1tab, *H2nkR1tab, H1nkr, alphaBT, betaBT;
2068
2069 double k, R0, R1, r, kr, kR0, kR1, theta, cost, sint ;
2070 int n, ns, ABCtype, SingleMode ;
2071
2072 theta = atan2(A->Val[1], A->Val[0]) ;
2073 r = sqrt(A->Val[0] * A->Val[0] + A->Val[1] * A->Val[1]) ;
2074
2075 k = Fct->Para[0] ;
2076 R0 = Fct->Para[1] ;
2077 R1 = Fct->Para[2] ;
2078 ABCtype = (int)Fct->Para[3] ;
2079 SingleMode = (int)Fct->Para[4] ;
2080 kr = k * r;
2081 kR0 = k * R0;
2082 kR1 = k * R1;
2083
2084 if(ABCtype == 1){ /* Sommerfeld */
2085 lambda = Cprodr(-k, I);
2086 }
2087 else if(ABCtype == 2){ /* Bayliss-Turkel */
2088 /*
2089 alphaBT[] = 1/(2*R1) - I[]/(8*k*R1^2*(1+I[]/(k*R1)));
2090 betaBT[] = - 1/(2*I[]*k*(1+I[]/(k*R1)));
2091 */
2092 coef.r = 2*k;
2093 coef.i = 2/R1;
2094 alphaBT = Csubr( 1/(2*R1) , Cdiv(I , Cprodr(4*R1*R1 , coef) ) );
2095 betaBT = Cdiv(I , coef);
2096 }
2097 else{
2098 Message::Error("Unknown ABC type");
2099 }
2100
2101 V->Val[0] = 0.;
2102 V->Val[MAX_DIM] = 0. ;
2103
2104 ns = (int)k + 10;
2105
2106 H1nkR0tab = (cplx*)Malloc(ns * sizeof(cplx));
2107 for (n = 0 ; n < ns ; n++){
2108 H1nkR0tab[n].r = jn(n, kR0);
2109 H1nkR0tab[n].i = yn(n, kR0);
2110 }
2111
2112 H2nkR0tab = (cplx*)Malloc(ns * sizeof(cplx));
2113 for (n = 0 ; n < ns ; n++){
2114 H2nkR0tab[n] = Cconj(H1nkR0tab[n]);
2115 }
2116
2117 H1nkR1tab = (cplx*)Malloc(ns * sizeof(cplx));
2118 for (n = 0 ; n < ns ; n++){
2119 H1nkR1tab[n].r = jn(n, kR1);
2120 H1nkR1tab[n].i = yn(n, kR1);
2121 }
2122
2123 H2nkR1tab = (cplx*)Malloc(ns * sizeof(cplx));
2124 for (n = 0 ; n < ns ; n++){
2125 H2nkR1tab[n] = Cconj(H1nkR1tab[n]);
2126 }
2127
2128 for (n = 0 ; n < ns ; n++){
2129 if(SingleMode >= 0 && SingleMode != n) continue;
2130
2131 H1nkr.r = jn(n,kr);
2132 H1nkr.i = yn(n,kr);
2133
2134 if(ABCtype == 2){
2135 lambda = Csum( Csum( Cprodr(-k, I) , alphaBT ) ,
2136 Cprodr( n*n/(R1*R1) , betaBT ) );
2137 }
2138
2139 alpha1 = Csum( Cprodr(k, DHn(H1nkR1tab, n, kR1)) ,
2140 Cprod(lambda, H1nkR1tab[n]) );
2141 alpha2 = Csum( Cprodr(k, DHn(H2nkR1tab, n, kR1)) ,
2142 Cprod(lambda, H2nkR1tab[n]) );
2143 delta = Cprodr( k , Csub( Cprod( alpha1 , DHn(H2nkR0tab, n, kR0) ) ,
2144 Cprod( alpha2 , DHn(H1nkR0tab, n, kR0) ) ) );
2145
2146 if(Cmodu(delta) < 1.e-6) break;
2147
2148 am = Cdiv( Cprodr(k * DHn(H1nkR0tab, n, kR0).r, alpha2) ,
2149 delta );
2150 bm = Cdiv( Cprodr(-k * DHn(H1nkR0tab, n, kR0).r, alpha1) ,
2151 delta );
2152
2153 if(SingleMode >= 0 && SingleMode == n){
2154 tmp = Csum( Cprod( am , H1nkr ) , Cprod( bm , Cconj(H1nkr) ) ) ;
2155 cost = cos(n * theta);
2156 sint = sin(n * theta);
2157 V->Val[0] += cost * tmp.r - sint * tmp.i;
2158 V->Val[MAX_DIM] += cost * tmp.i + sint * tmp.r;
2159 }
2160 else{
2161 tmp = Cprod( Cpow(I,n) , Csum( Cprod( am , H1nkr ) ,
2162 Cprod( bm , Cconj(H1nkr) ) ) );
2163 cost = cos(n * theta);
2164 V->Val[0] += cost * tmp.r * (!n ? 0.5 : 1.);
2165 V->Val[MAX_DIM] += cost * tmp.i * (!n ? 0.5 : 1.);
2166 }
2167 }
2168
2169 Free(H1nkR0tab);
2170 Free(H2nkR0tab);
2171 Free(H1nkR1tab);
2172 Free(H2nkR1tab);
2173
2174 if(SingleMode < 0){
2175 V->Val[0] *= 2;
2176 V->Val[MAX_DIM] *= 2;
2177 }
2178
2179 V->Type = SCALAR ;
2180 }
2181
2182 /* Scattering by acoustically hard circular cylinder of radius R,
2183 under plane wave incidence e^{ikx}. Returns RCS. */
2184
F_RCSHardCylinder(F_ARG)2185 void F_RCSHardCylinder(F_ARG)
2186 {
2187 cplx I = {0.,1.}, Hnkr, dHnkR, res, tmp, *HnkRtab;
2188 double k, R, r, kR, theta, cost, val ;
2189 int n, ns ;
2190
2191 theta = atan2(A->Val[1], A->Val[0]) ;
2192 r = sqrt(A->Val[0]*A->Val[0] + A->Val[1]*A->Val[1]) ;
2193
2194 k = Fct->Para[0] ;
2195 R = Fct->Para[1] ;
2196 kR = k*R;
2197
2198 res.r = 0.;
2199 res.i = 0. ;
2200
2201 ns = (int)k + 10;
2202
2203 HnkRtab = (cplx*)Malloc(ns*sizeof(cplx));
2204
2205 for (n = 0 ; n < ns ; n++){
2206 HnkRtab[n].r = jn(n,kR);
2207 HnkRtab[n].i = yn(n,kR);
2208 }
2209
2210 for (n = 0 ; n < ns ; n++){
2211 /* leaving r in following asymptotic formula for clarity (see
2212 Colton and Kress, Inverse Acoustic..., p. 65, eq. 3.59) */
2213 Hnkr.r = cos(k*r - n*M_PI/2. - M_PI/4.) / sqrt(k*r) * sqrt(2./M_PI);
2214 Hnkr.i = sin(k*r - n*M_PI/2. - M_PI/4.) / sqrt(k*r) * sqrt(2./M_PI);
2215
2216 dHnkR = DHn(HnkRtab, n, kR);
2217
2218 tmp = Cdiv( Cprod( Cpow(I,n) , Cprodr( dHnkR.r, Hnkr) ) , dHnkR );
2219
2220 cost = cos(n*theta);
2221
2222 res.r += cost * tmp.r * (!n ? 0.5 : 1.);
2223 res.i += cost * tmp.i * (!n ? 0.5 : 1.);
2224 }
2225
2226 Free(HnkRtab);
2227
2228 res.r *= -2;
2229 res.i *= -2;
2230
2231 val = Cmodu(res);
2232 val *= val;
2233 val *= 2. * M_PI * r;
2234 val = 10. * log10(val);
2235
2236 V->Val[0] = val;
2237 V->Val[MAX_DIM] = 0.;
2238
2239 V->Type = SCALAR ;
2240 }
2241
2242 /* ------------------------------------------------------------------------ */
2243 /* On Surface Radiation Conditions (OSRC) */
2244 /* ------------------------------------------------------------------------ */
2245
2246 /* Coefficients C0, Aj and Bj: see papers
2247 1) Kechroud, Antoine & Soulaimani, Nuemrical accuracy of a
2248 Pade-type non-reflecting..., IJNME 2005
2249 2) Antoine, Darbas & Lu, An improved surface radiation condition...
2250 CMAME, 2006(?) */
2251
aj(int j,int N)2252 static double aj(int j, int N)
2253 {
2254 return 2./(2.*N + 1.) * SQU(sin((double)j * M_PI/(2.*N + 1.))) ;
2255 }
2256
bj(int j,int N)2257 static double bj(int j, int N)
2258 {
2259 return SQU(cos((double)j * M_PI/(2.*N + 1.))) ;
2260 }
2261
padeC0(int N,double theta)2262 static std::complex<double> padeC0(int N, double theta)
2263 {
2264 std::complex<double> sum = std::complex<double>(1, 0);
2265 std::complex<double> one = std::complex<double>(1, 0);
2266 std::complex<double> z = std::complex<double>(cos(-theta) - 1, sin(-theta));
2267
2268 for(int j = 1; j <= N; j++)
2269 sum += (z * aj(j, N)) / (one + z * bj(j, N));
2270
2271 z = std::complex<double>(cos(theta / 2.), sin(theta / 2.));
2272
2273 return sum * z;
2274 }
2275
padeA(int j,int N,double theta)2276 static std::complex<double> padeA(int j, int N, double theta)
2277 {
2278 std::complex<double> one = std::complex<double>(1, 0);
2279 std::complex<double> res;
2280 std::complex<double> z;
2281
2282 z = std::complex<double>(cos(-theta / 2.), sin(-theta / 2.));
2283 res = z * aj(j, N);
2284
2285 z = std::complex<double>(cos(-theta) - 1., sin(-theta));
2286 res = res / ((one + z * bj(j, N)) * (one + z * bj(j, N)));
2287
2288 return res;
2289 }
2290
padeB(int j,int N,double theta)2291 static std::complex<double> padeB(int j, int N, double theta)
2292 {
2293 std::complex<double> one = std::complex<double>(1, 0);
2294 std::complex<double> res;
2295 std::complex<double> z;
2296
2297 z = std::complex<double>(cos(-theta), sin(-theta));
2298 res = z * bj(j, N);
2299
2300 z = std::complex<double>(cos(-theta) - 1., sin(-theta));
2301 res = res / (one + z * bj(j, N));
2302
2303 return res;
2304 }
2305
padeR0(int N,double theta)2306 static std::complex<double> padeR0(int N, double theta)
2307 {
2308 std::complex<double> sum = padeC0(N, theta);
2309
2310 for(int j = 1; j <= N; j++)
2311 sum += padeA(j, N, theta) / padeB(j, N, theta);
2312
2313 return sum;
2314 }
2315
F_OSRC_C0(F_ARG)2316 void F_OSRC_C0(F_ARG)
2317 {
2318 int N;
2319 double theta;
2320
2321 N = (int)Fct->Para[0];
2322 theta = Fct->Para[1];
2323
2324 std::complex<double> C0 = padeC0(N, theta);
2325
2326 V->Val[0] = C0.real();
2327 V->Val[MAX_DIM] = C0.imag();
2328 V->Type = SCALAR;
2329 }
2330
F_OSRC_R0(F_ARG)2331 void F_OSRC_R0(F_ARG)
2332 {
2333 int N;
2334 double theta;
2335
2336 N = (int)Fct->Para[0];
2337 theta = Fct->Para[1];
2338
2339 std::complex<double> C0 = padeR0(N, theta);
2340
2341 V->Val[0] = C0.real();
2342 V->Val[MAX_DIM] = C0.imag();
2343 V->Type = SCALAR;
2344 }
2345
F_OSRC_Aj(F_ARG)2346 void F_OSRC_Aj(F_ARG)
2347 {
2348 int j, N;
2349 double theta;
2350
2351 j = (int)Fct->Para[0];
2352 N = (int)Fct->Para[1];
2353 theta = Fct->Para[2];
2354
2355 std::complex<double> Aj = padeA(j, N, theta);
2356
2357 V->Val[0] = Aj.real();
2358 V->Val[MAX_DIM] = Aj.imag();
2359 V->Type = SCALAR;
2360 }
2361
F_OSRC_Bj(F_ARG)2362 void F_OSRC_Bj(F_ARG)
2363 {
2364 int j, N;
2365 double theta;
2366
2367 j = (int)Fct->Para[0];
2368 N = (int)Fct->Para[1];
2369 theta = Fct->Para[2];
2370
2371 std::complex<double> Bj = padeB(j, N, theta);
2372
2373 V->Val[0] = Bj.real();
2374 V->Val[MAX_DIM] = Bj.imag();
2375 V->Type = SCALAR;
2376 }
2377
2378 /* --------------------------------------------------------------------- */
2379 /* Vector spherical harmonics (aka Vector Partial Waves) */
2380 /* time sign : beware! */
2381 /* Implemented following notations and conventions in : */
2382 /* B. Stout, "Spherical harmonic Lattice Sums for Gratings", */
2383 /* Ch. 6 of "Gratings: Theory and Numeric Applications", */
2384 /* (see Eq. (6.200), p.6.37) in Chap. 6 of */
2385 /* "Gratings: Theory and Numeric Applications", */
2386 /* Ed. E. Popov (Institut Fresnel, CNRS, AMU, 2012) */
2387 /* See also Jackson, 3rd Ed., Chap. 9, p. 431... */
2388 /* See also Chew, p. 185... */
2389 /* --------------------------------------------------------------------- */
2390
sph_pnm(int n,int m,double u)2391 double sph_pnm(int n, int m, double u)
2392 {
2393 int k, kn, km;
2394 double knf, kmf, **pmntab;
2395
2396 if(abs(m) > n)
2397 Message::Error("|m|<=n for the normalized pnm's");
2398
2399 pmntab = (double**)Malloc( (2*n+1) * sizeof(double *));
2400 for (k = 0; k <= 2*n ; k++){
2401 pmntab[k] = (double*)Malloc((n+1) * sizeof(double));
2402 }
2403 for (km = 0; km <= 2*n ; km++){
2404 for (kn = 0; kn <= n ; kn++){
2405 pmntab[km][kn]=0.;
2406 }
2407 }
2408
2409 // initialize recur.
2410 pmntab[n][0]=1./sqrt(4.*M_PI);
2411
2412 // fill diag and first diag
2413 for (kn = 1; kn <= n ; kn++){
2414 knf = (double)kn;
2415 pmntab[n+kn][kn] = -sqrt((2.*knf+1.)/(2.*knf)) * sqrt(1.-pow(u,2))
2416 * pmntab[n+kn-1][kn-1];
2417 pmntab[n+kn-1][kn] = u*sqrt(2.*knf+1.) * pmntab[n+kn-1][kn-1];
2418 }
2419 for (kn = 2; kn <= n ; kn++){
2420 for (km = 0; km <= kn-2 ; km++){
2421 knf = (double)kn;
2422 kmf = (double)km;
2423 pmntab[n+km][kn] = sqrt((4.*knf*knf-1.)/(pow(knf,2)-pow(kmf,2)))
2424 * pmntab[n+km][kn-1] * u
2425 - sqrt(((2.*knf+1.)*(pow(knf-1.,2)-pow(kmf,2)))
2426 / ((2.*knf-3.)*(pow(knf,2)-pow(kmf,2))))
2427 * pmntab[n+km][kn-2];
2428 }
2429 }
2430 for (kn = 0; kn <= n ; kn++){
2431 for (km = n-1; km >= 0 ; km--){
2432 pmntab[km][kn] = pow(-1,n-km)*pmntab[2*n-km][kn];
2433 }
2434 }
2435 double res = pmntab[n+m][n];
2436 for (k = 0; k <= 2*n ; k++){
2437 Free(pmntab[k]);
2438 }
2439 Free(pmntab);
2440 return res;
2441 }
2442
sph_unm(int n,int m,double u)2443 double sph_unm(int n, int m, double u)
2444 {
2445 int k, kn, km;
2446 double knf, kmf, **umntab;
2447
2448 if(abs(m) > n)
2449 Message::Error("|m|<=n for the normalized unm's");
2450 umntab = (double**)Malloc( (2*n+1) * sizeof(double *));
2451 for (k = 0; k <= 2*n ; k++){
2452 umntab[k] = (double*)Malloc((n+1) * sizeof(double));
2453 }
2454
2455 for (km = 0; km <= 2*n ; km++){
2456 for (kn = 0; kn <= n ; kn++){
2457 umntab[km][kn]=0.;
2458 }
2459 }
2460 // initialize recur.
2461 umntab[n][0]=0.;
2462 umntab[n+1][1]=-0.25*sqrt(3./M_PI);
2463 // fill diag and first diag
2464 for (kn = 2; kn <= n ; kn++){
2465 knf = (double)kn;
2466 umntab[n+kn][kn] = -sqrt((knf*(2.*knf+1.))
2467 /(2.*(knf+1.)*(knf-1.)))
2468 *sqrt(1.-pow(u,2))*umntab[n+kn-1][kn-1];
2469 umntab[n+kn-1][kn] = sqrt((2.*knf+1.)*(knf-1.)/(knf+1.))
2470 *u* umntab[n+kn-1][kn-1];
2471 }
2472 for (kn = 2; kn <= n ; kn++){
2473 for (km = 0; km <= kn-2 ; km++){
2474 knf = (double)kn;
2475 kmf = (double)km;
2476 umntab[n+km][kn] = sqrt(((4.*pow(knf,2)-1.)*(knf-1.))
2477 /((pow(knf,2)-pow(kmf,2))*(knf+1.)))
2478 *umntab[n+km][kn-1]*u
2479 -sqrt(((2.*knf+1.)*(knf-1.)*(knf-2.)
2480 *(knf-kmf-1.)*(knf+kmf-1.))
2481 /((2.*knf-3.)*(pow(knf,2)-pow(kmf,2))*knf*(knf+1.)))
2482 *umntab[n+km][kn-2];
2483 }
2484 }
2485 for (kn = 0; kn <= n ; kn++){
2486 for (km = n-1; km >= 0 ; km--){
2487 umntab[km][kn] = pow(-1,n-km+1)*umntab[2*n-km][kn];
2488 }
2489 }
2490 double res = umntab[n+m][n];
2491 for (k = 0; k <= 2*n ; k++){
2492 Free(umntab[k]);
2493 }
2494 Free(umntab);
2495 return res;
2496 }
2497
sph_snm(int n,int m,double u)2498 double sph_snm(int n, int m, double u)
2499 {
2500 int k, kn, km;
2501 double knf, kmf, **umntab, **smntab;
2502 if(abs(m) > n)
2503 Message::Error("|m|<=n for the normalized snm's");
2504 umntab = (double**)Malloc( (2*n+1) * sizeof(double *));
2505 smntab = (double**)Malloc( (2*n+1) * sizeof(double *));
2506 for (k = 0; k <= 2*n ; k++){
2507 umntab[k] = (double*)Malloc((n+1) * sizeof(double));
2508 smntab[k] = (double*)Malloc((n+1) * sizeof(double));
2509 }
2510 for (km = 0; km <= 2*n ; km++){
2511 for (kn = 0; kn <= n ; kn++){
2512 umntab[km][kn]=0.;
2513 smntab[km][kn]=0.;
2514 }
2515 }
2516
2517 // initialize recur.
2518 umntab[n][0]=0.;
2519 umntab[n+1][1]=-0.25*sqrt(3./M_PI);
2520 // fill diag and first diag
2521 for (kn = 2; kn <= n ; kn++){
2522 knf = (double)kn;
2523 umntab[n+kn][kn] =-sqrt((knf*(2.*knf+1.))/(2.*(knf+1.)*(knf-1.)))
2524 *sqrt(1.-pow(u,2))*umntab[n+kn-1][kn-1];
2525 umntab[n+kn-1][kn]= sqrt((2.*knf+1.)*(knf-1.)/(knf+1.))
2526 *u*umntab[n+kn-1][kn-1];
2527 }
2528 for (kn = 2; kn <= n ; kn++){
2529 for (km = 0; km <= kn-2 ; km++){
2530 knf = (double)kn;
2531 kmf = (double)km;
2532 umntab[n+km][kn] = sqrt(((4.*pow(knf,2)-1.)*(knf-1.))
2533 /((pow(knf,2)-pow(kmf,2))*(knf+1.)))
2534 *umntab[n+km][kn-1]*u
2535 -sqrt(((2.*knf+1.)*(knf-1.)*(knf-2.)
2536 *(knf-kmf-1.)*(knf+kmf-1.))
2537 /((2.*knf-3.)*(pow(knf,2)
2538 -pow(kmf,2))*knf*(knf+1.)))
2539 *umntab[n+km][kn-2];
2540 }
2541 }
2542 for (kn = 0; kn <= n ; kn++){
2543 for (km = n-1; km >= 0 ; km--){
2544 umntab[km][kn] = pow(-1,n-km+1)*umntab[2*n-km][kn];
2545 }
2546 }
2547 for (kn = 0; kn <= n ; kn++){
2548 km = 0;
2549 kmf = 0.;
2550 knf = (double)kn;
2551 smntab[n+km][kn]=1./(kmf+1)*sqrt((knf+kmf+1.)*(knf-kmf))
2552 *sqrt(1.-pow(u,2))* umntab[n+km+1][kn]
2553 +u*umntab[n+km][kn];
2554 }
2555 for (kn = 1; kn <= n ; kn++){
2556 for (km = 1; km <= kn ; km++){
2557 knf = (double)kn;
2558 kmf = (double)km;
2559 smntab[n+km][kn] = knf/kmf*u*umntab[n+km][kn]-(kmf+knf)/kmf*
2560 sqrt(((2.*knf+1.)*(knf-kmf)*(knf-1.))
2561 /((2.*knf-1.)*(knf+kmf)*(knf+1.)))
2562 *umntab[n+km][kn-1];
2563 }
2564 }
2565 for (kn = 0; kn <= n ; kn++){
2566 for (km = n-1; km >= 0 ; km--){
2567 smntab[km][kn] = pow(-1,n-km)*smntab[2*n-km][kn];
2568 }
2569 }
2570 double res = smntab[n+m][n];
2571 for (k = 0; k <= 2*n ; k++){
2572 Free(umntab[k]);
2573 Free(smntab[k]);
2574 }
2575 Free(smntab);
2576 return res;
2577 }
2578
F_pnm(F_ARG)2579 void F_pnm(F_ARG)
2580 {
2581 int n, m;
2582 double u;
2583 if(A->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
2584 Message::Error("Non scalar argument(s) for the normalized pnm's");
2585 n = (int)A->Val[0];
2586 m = (int)(A+1)->Val[0];
2587 u = (A+2)->Val[0];
2588 V->Val[0] = sph_pnm(n,m,u);
2589 V->Type = SCALAR;
2590 }
F_unm(F_ARG)2591 void F_unm(F_ARG)
2592 {
2593 int n, m;
2594 double u;
2595 if(A->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
2596 Message::Error("Non scalar argument(s) for the normalized unm's");
2597 n = (int)A->Val[0];
2598 m = (int)(A+1)->Val[0];
2599 u = (A+2)->Val[0];
2600 V->Val[0] = sph_unm(n,m,u);
2601 V->Type = SCALAR;
2602 }
2603
F_snm(F_ARG)2604 void F_snm(F_ARG)
2605 {
2606 int n, m;
2607 double u;
2608 if(A->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
2609 Message::Error("Non scalar argument(s) for the normalized snm's");
2610 n = (int)A->Val[0];
2611 m = (int)(A+1)->Val[0];
2612 u = (A+2)->Val[0];
2613 V->Val[0] = sph_snm(n,m,u);
2614 V->Type = SCALAR;
2615 }
2616
F_Xnm(F_ARG)2617 void F_Xnm(F_ARG)
2618 {
2619 int n, m;
2620 double x, y, z;
2621 double r, theta, phi;
2622 double Xnm_r_re,Xnm_r_im,Xnm_t_re,Xnm_t_im,Xnm_p_re,Xnm_p_im;
2623 double Xnm_x_re,Xnm_x_im,Xnm_y_re,Xnm_y_im,Xnm_z_re,Xnm_z_im;
2624 double costheta, unm_costheta, snm_costheta;
2625 double exp_jm_phi_re, exp_jm_phi_im;
2626 double avoid_r_singularity = 1.E-13;
2627
2628 if( A->Type != SCALAR
2629 || (A+1)->Type != SCALAR
2630 || (A+2)->Type != SCALAR
2631 || (A+3)->Type != SCALAR
2632 || (A+4)->Type != SCALAR)
2633 Message::Error("Non scalar argument(s) for the Mnm's");
2634 n = (int)A->Val[0];
2635 m = (int)(A+1)->Val[0];
2636 x = (A+2)->Val[0];
2637 y = (A+3)->Val[0];
2638 z = (A+4)->Val[0];
2639 // k0 = (A+5)->Val[0];
2640
2641 if(n < 0) Message::Error("n should be a positive integer for the Xnm's");
2642 if(abs(m) > n) Message::Error("|m|<=n is required for the Xnm's");
2643
2644 r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
2645 if (r<avoid_r_singularity) r=avoid_r_singularity;
2646 costheta = z/r;
2647 theta = acos(costheta);
2648 phi = atan2(-y,-x)+M_PI;
2649
2650 unm_costheta = sph_unm(n,m,costheta);
2651 snm_costheta = sph_snm(n,m,costheta);
2652 exp_jm_phi_re = cos( ((double)m) *phi);
2653 exp_jm_phi_im = sin( ((double)m) *phi);
2654
2655 // in spherical coord
2656 Xnm_r_re = 0.;
2657 Xnm_r_im = 0.;
2658 Xnm_t_re = -1.*unm_costheta * exp_jm_phi_im;
2659 Xnm_t_im = unm_costheta * exp_jm_phi_re;
2660 Xnm_p_re = -1.*snm_costheta * exp_jm_phi_re;
2661 Xnm_p_im = -1.*snm_costheta * exp_jm_phi_im;
2662
2663 // in cart coord
2664 Xnm_x_re = sin(theta)*cos(phi)*Xnm_r_re
2665 + cos(theta)*cos(phi)*Xnm_t_re
2666 - sin(phi)*Xnm_p_re ;
2667 Xnm_y_re = sin(theta)*sin(phi)*Xnm_r_re
2668 + cos(theta)*sin(phi)*Xnm_t_re
2669 + cos(phi)*Xnm_p_re ;
2670 Xnm_z_re = cos(theta)*Xnm_r_re
2671 - sin(theta)*Xnm_t_re;
2672 Xnm_x_im = sin(theta)*cos(phi)*Xnm_r_im
2673 + cos(theta)*cos(phi)*Xnm_t_im
2674 - sin(phi)*Xnm_p_im ;
2675 Xnm_y_im = sin(theta)*sin(phi)*Xnm_r_im
2676 + cos(theta)*sin(phi)*Xnm_t_im
2677 + cos(phi)*Xnm_p_im ;
2678 Xnm_z_im = cos(theta)*Xnm_r_im
2679 - sin(theta)*Xnm_t_im;
2680
2681 V->Type = VECTOR;
2682 V->Val[0] = Xnm_x_re ; V->Val[MAX_DIM ] = Xnm_x_im ;
2683 V->Val[1] = Xnm_y_re ; V->Val[MAX_DIM+1] = Xnm_y_im ;
2684 V->Val[2] = Xnm_z_re ; V->Val[MAX_DIM+2] = Xnm_z_im ;
2685 }
F_Ynm(F_ARG)2686 void F_Ynm(F_ARG)
2687 {
2688 int n, m;
2689 double x, y, z;
2690 double r, theta, phi;
2691 double Ynm_r_re,Ynm_r_im,Ynm_t_re,Ynm_t_im,Ynm_p_re,Ynm_p_im;
2692 double Ynm_x_re,Ynm_x_im,Ynm_y_re,Ynm_y_im,Ynm_z_re,Ynm_z_im;
2693 double costheta, pnm_costheta;
2694 double exp_jm_phi_re, exp_jm_phi_im;
2695 double avoid_r_singularity = 1.E-13;
2696
2697 if( A->Type != SCALAR
2698 || (A+1)->Type != SCALAR
2699 || (A+2)->Type != SCALAR
2700 || (A+3)->Type != SCALAR
2701 || (A+4)->Type != SCALAR)
2702 Message::Error("Non scalar argument(s) for the Mnm's");
2703 n = (int)A->Val[0];
2704 m = (int)(A+1)->Val[0];
2705 x = (A+2)->Val[0];
2706 y = (A+3)->Val[0];
2707 z = (A+4)->Val[0];
2708 // k0 = (A+5)->Val[0];
2709
2710 if(n < 0) Message::Error("n should be a positive integer for the Ynm's");
2711 if(abs(m) > n) Message::Error("|m|<=n is required for the Ynm's");
2712
2713 r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
2714 if (r<avoid_r_singularity) r=avoid_r_singularity;
2715 costheta = z/r;
2716 theta = acos(costheta);
2717 phi = atan2(-y,-x)+M_PI;
2718
2719 pnm_costheta = sph_pnm(n,m,costheta);
2720 // unm_costheta = sph_unm(n,m,costheta);
2721 // snm_costheta = sph_snm(n,m,costheta);
2722 exp_jm_phi_re = cos( ((double)m) *phi);
2723 exp_jm_phi_im = sin( ((double)m) *phi);
2724
2725 Ynm_r_re = pnm_costheta * exp_jm_phi_re;
2726 Ynm_r_im = pnm_costheta * exp_jm_phi_im;
2727 Ynm_t_re = 0.;
2728 Ynm_t_im = 0.;
2729 Ynm_p_re = 0.;
2730 Ynm_p_im = 0.;
2731
2732 Ynm_x_re = sin(theta)*cos(phi)*Ynm_r_re
2733 + cos(theta)*cos(phi)*Ynm_t_re
2734 - sin(phi)*Ynm_p_re ;
2735 Ynm_y_re = sin(theta)*sin(phi)*Ynm_r_re
2736 + cos(theta)*sin(phi)*Ynm_t_re
2737 + cos(phi)*Ynm_p_re ;
2738 Ynm_z_re = cos(theta) *Ynm_r_re
2739 - sin(theta) *Ynm_t_re;
2740 Ynm_x_im = sin(theta)*cos(phi)*Ynm_r_im
2741 + cos(theta)*cos(phi)*Ynm_t_im
2742 - sin(phi)*Ynm_p_im ;
2743 Ynm_y_im = sin(theta)*sin(phi)*Ynm_r_im
2744 + cos(theta)*sin(phi)*Ynm_t_im
2745 + cos(phi)*Ynm_p_im ;
2746 Ynm_z_im = cos(theta) *Ynm_r_im
2747 - sin(theta) *Ynm_t_im;
2748
2749 V->Type = VECTOR;
2750 V->Val[0] = Ynm_x_re ; V->Val[MAX_DIM ] = Ynm_x_im ;
2751 V->Val[1] = Ynm_y_re ; V->Val[MAX_DIM+1] = Ynm_y_im ;
2752 V->Val[2] = Ynm_z_re ; V->Val[MAX_DIM+2] = Ynm_z_im ;
2753 }
2754
F_Znm(F_ARG)2755 void F_Znm(F_ARG)
2756 {
2757 int n, m;
2758 double x, y, z;
2759 double r, theta, phi;
2760 double Znm_r_re,Znm_r_im,Znm_t_re,Znm_t_im,Znm_p_re,Znm_p_im;
2761 double Znm_x_re,Znm_x_im,Znm_y_re,Znm_y_im,Znm_z_re,Znm_z_im;
2762 double costheta, unm_costheta, snm_costheta;
2763 double exp_jm_phi_re, exp_jm_phi_im;
2764 double avoid_r_singularity = 1.E-13;
2765
2766 if( A->Type != SCALAR
2767 || (A+1)->Type != SCALAR
2768 || (A+2)->Type != SCALAR
2769 || (A+3)->Type != SCALAR
2770 || (A+4)->Type != SCALAR)
2771 Message::Error("Non scalar argument(s) for the Mnm's");
2772 n = (int)A->Val[0];
2773 m = (int)(A+1)->Val[0];
2774 x = (A+2)->Val[0];
2775 y = (A+3)->Val[0];
2776 z = (A+4)->Val[0];
2777 // k0 = (A+5)->Val[0];
2778
2779 if(n < 0)
2780 Message::Error("n should be a positive integer for the Znm's");
2781 if(abs(m) > n)
2782 Message::Error("|m|<=n is required for the Znm's");
2783
2784 r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
2785 if (r<avoid_r_singularity) r=avoid_r_singularity;
2786 costheta = z/r;
2787 theta = acos(costheta);
2788 phi = atan2(-y,-x)+M_PI;
2789
2790 unm_costheta = sph_unm(n,m,costheta);
2791 snm_costheta = sph_snm(n,m,costheta);
2792 exp_jm_phi_re = cos( ((double)m) *phi);
2793 exp_jm_phi_im = sin( ((double)m) *phi);
2794
2795 Znm_r_re = 0.;
2796 Znm_r_im = 0.;
2797 Znm_t_re = snm_costheta * exp_jm_phi_re;
2798 Znm_t_im = snm_costheta * exp_jm_phi_im;
2799 Znm_p_re = -1.*unm_costheta * exp_jm_phi_im;
2800 Znm_p_im = unm_costheta * exp_jm_phi_re;
2801
2802 Znm_x_re = sin(theta)*cos(phi)*Znm_r_re
2803 + cos(theta)*cos(phi)*Znm_t_re
2804 - sin(phi)*Znm_p_re ;
2805 Znm_y_re = sin(theta)*sin(phi)*Znm_r_re
2806 + cos(theta)*sin(phi)*Znm_t_re
2807 + cos(phi)*Znm_p_re ;
2808 Znm_z_re = cos(theta)*Znm_r_re
2809 - sin(theta)*Znm_t_re;
2810 Znm_x_im = sin(theta)*cos(phi)*Znm_r_im
2811 + cos(theta)*cos(phi)*Znm_t_im
2812 - sin(phi)*Znm_p_im ;
2813 Znm_y_im = sin(theta)*sin(phi)*Znm_r_im
2814 + cos(theta)*sin(phi)*Znm_t_im
2815 + cos(phi)*Znm_p_im ;
2816 Znm_z_im = cos(theta)*Znm_r_im
2817 - sin(theta)*Znm_t_im;
2818
2819 V->Type = VECTOR;
2820 V->Val[0] = Znm_x_re ; V->Val[MAX_DIM ] = Znm_x_im ;
2821 V->Val[1] = Znm_y_re ; V->Val[MAX_DIM+1] = Znm_y_im ;
2822 V->Val[2] = Znm_z_re ; V->Val[MAX_DIM+2] = Znm_z_im ;
2823 }
2824
F_Mnm(F_ARG)2825 void F_Mnm(F_ARG)
2826 {
2827 int Mtype, n, m;
2828 double x, y, z, k0;
2829 double r=0., theta=0., phi=0.;
2830 double Xnm_t_re,Xnm_t_im,Xnm_p_re,Xnm_p_im;
2831 double Mnm_r_re,Mnm_r_im,Mnm_t_re,Mnm_t_im,Mnm_p_re,Mnm_p_im;
2832 double Mnm_x_re,Mnm_x_im,Mnm_y_re,Mnm_y_im,Mnm_z_re,Mnm_z_im;
2833 double costheta, unm_costheta, snm_costheta;
2834 double exp_jm_phi_re, exp_jm_phi_im;
2835 double sph_bessel_n_ofkr_re, sph_bessel_n_ofkr_im;
2836 double avoid_r_singularity = 1.E-13;
2837
2838 if( A->Type != SCALAR
2839 || (A+1)->Type != SCALAR
2840 || (A+2)->Type != SCALAR
2841 || (A+3)->Type != VECTOR
2842 || (A+4)->Type != SCALAR)
2843 Message::Error("Check types of arguments for the Mnm's");
2844 Mtype = (int)A->Val[0];
2845 n = (int)(A+1)->Val[0];
2846 m = (int)(A+2)->Val[0];
2847 x = (A+3)->Val[0];
2848 y = (A+3)->Val[1];
2849 z = (A+3)->Val[2];
2850 k0 = (A+4)->Val[0];
2851
2852 if(n < 0) Message::Error("n should be a positive integer for the Mnm's");
2853 if(abs(m) > n) Message::Error("|m|<=n is required for the Mnm's");
2854
2855 r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
2856 if (r<avoid_r_singularity) r=avoid_r_singularity;
2857 costheta = z/r;
2858 theta = acos(costheta);
2859 phi = atan2(-y,-x)+M_PI;
2860
2861 // pnm_costheta = sph_pnm(n,m,costheta);
2862 unm_costheta = sph_unm(n,m,costheta);
2863 snm_costheta = sph_snm(n,m,costheta);
2864 exp_jm_phi_re = cos( ((double)m) *phi);
2865 exp_jm_phi_im = sin( ((double)m) *phi);
2866
2867 Xnm_t_re = -1.*unm_costheta * exp_jm_phi_im;
2868 Xnm_t_im = unm_costheta * exp_jm_phi_re;
2869 Xnm_p_re = -1.*snm_costheta * exp_jm_phi_re;
2870 Xnm_p_im = -1.*snm_costheta * exp_jm_phi_im;
2871
2872 switch (Mtype) {
2873 case 1: // Spherical Bessel function of the first kind j_n
2874 sph_bessel_n_ofkr_re = Spherical_j_n(n, k0*r);
2875 sph_bessel_n_ofkr_im = 0;
2876 break;
2877 case 2: // Spherical Bessel function of the second kind y_n
2878 sph_bessel_n_ofkr_re = Spherical_y_n(n, k0*r);
2879 sph_bessel_n_ofkr_im = 0;
2880 break;
2881 case 3: // Spherical Hankel function of the first kind h^1_n
2882 sph_bessel_n_ofkr_re = Spherical_j_n(n, k0*r);
2883 sph_bessel_n_ofkr_im = Spherical_y_n(n, k0*r);
2884 break;
2885 case 4: // Spherical Hankel function of the second kind h^2_n
2886 sph_bessel_n_ofkr_re = Spherical_j_n(n, k0*r);
2887 sph_bessel_n_ofkr_im =-Spherical_y_n(n, k0*r);
2888 break;
2889 default:
2890 sph_bessel_n_ofkr_re =0.;
2891 sph_bessel_n_ofkr_im =0.;
2892 Message::Error("First argument for Nnm's should be 1,2,3 or 4");
2893 break;
2894 }
2895
2896 Mnm_r_re = 0.;
2897 Mnm_r_im = 0.;
2898 Mnm_t_re = sph_bessel_n_ofkr_re*Xnm_t_re - sph_bessel_n_ofkr_im*Xnm_t_im;
2899 Mnm_t_im = sph_bessel_n_ofkr_re*Xnm_t_im + sph_bessel_n_ofkr_im*Xnm_t_re;
2900 Mnm_p_re = sph_bessel_n_ofkr_re*Xnm_p_re - sph_bessel_n_ofkr_im*Xnm_p_im;
2901 Mnm_p_im = sph_bessel_n_ofkr_re*Xnm_p_im + sph_bessel_n_ofkr_im*Xnm_p_re;
2902
2903 Mnm_x_re = sin(theta)*cos(phi)*Mnm_r_re
2904 + cos(theta)*cos(phi)*Mnm_t_re
2905 - sin(phi)*Mnm_p_re ;
2906 Mnm_y_re = sin(theta)*sin(phi)*Mnm_r_re
2907 + cos(theta)*sin(phi)*Mnm_t_re
2908 + cos(phi)*Mnm_p_re ;
2909 Mnm_z_re = cos(theta)*Mnm_r_re
2910 - sin(theta)*Mnm_t_re;
2911 Mnm_x_im = sin(theta)*cos(phi)*Mnm_r_im
2912 + cos(theta)*cos(phi)*Mnm_t_im
2913 - sin(phi)*Mnm_p_im ;
2914 Mnm_y_im = sin(theta)*sin(phi)*Mnm_r_im
2915 + cos(theta)*sin(phi)*Mnm_t_im
2916 + cos(phi)*Mnm_p_im ;
2917 Mnm_z_im = cos(theta)*Mnm_r_im
2918 - sin(theta)*Mnm_t_im;
2919
2920 V->Type = VECTOR;
2921 V->Val[0] = Mnm_x_re ; V->Val[MAX_DIM ] = Mnm_x_im ;
2922 V->Val[1] = Mnm_y_re ; V->Val[MAX_DIM+1] = Mnm_y_im ;
2923 V->Val[2] = Mnm_z_re ; V->Val[MAX_DIM+2] = Mnm_z_im ;
2924 }
2925
F_Nnm(F_ARG)2926 void F_Nnm(F_ARG)
2927 {
2928 int Ntype, n, m;
2929 double x, y, z, k0;
2930 double r, theta, phi;
2931 double Ynm_r_re,Ynm_r_im;
2932 double Znm_t_re,Znm_t_im,Znm_p_re,Znm_p_im;
2933 double Nnm_r_re,Nnm_t_re,Nnm_p_re,Nnm_r_im,Nnm_t_im,Nnm_p_im;
2934 double Nnm_x_re,Nnm_y_re,Nnm_z_re,Nnm_x_im,Nnm_y_im,Nnm_z_im;
2935 double costheta, pnm_costheta, unm_costheta, snm_costheta;
2936 double exp_jm_phi_re, exp_jm_phi_im;
2937 double sph_bessel_n_ofkr_re, sph_bessel_nminus1_ofkr_re, dRicatti_dx_ofkr_re;
2938 double sph_bessel_n_ofkr_im, sph_bessel_nminus1_ofkr_im, dRicatti_dx_ofkr_im;
2939 double avoid_r_singularity = 1.E-13;
2940
2941 if( A->Type != SCALAR
2942 || (A+1)->Type != SCALAR
2943 || (A+2)->Type != SCALAR
2944 || (A+3)->Type != VECTOR
2945 || (A+4)->Type != SCALAR)
2946 Message::Error("Check types of arguments for the Mnm's");
2947 Ntype = (int)A->Val[0];
2948 n = (int)(A+1)->Val[0];
2949 m = (int)(A+2)->Val[0];
2950 x = (A+3)->Val[0];
2951 y = (A+3)->Val[1];
2952 z = (A+3)->Val[2];
2953 k0 = (A+4)->Val[0];
2954
2955 if(n < 0) Message::Error("n should be a positive integer for the Nnm's");
2956 if(abs(m) > n) Message::Error("|m|<=n is required for the Nnm's");
2957
2958 r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
2959 if (r<avoid_r_singularity) r=avoid_r_singularity;
2960 costheta = z/r;
2961 theta = acos(costheta);
2962 phi = atan2(-y,-x)+M_PI;
2963
2964 pnm_costheta = sph_pnm(n,m,costheta);
2965 unm_costheta = sph_unm(n,m,costheta);
2966 snm_costheta = sph_snm(n,m,costheta);
2967 exp_jm_phi_re = cos((double)m *phi);
2968 exp_jm_phi_im = sin((double)m *phi);
2969
2970 Ynm_r_re = pnm_costheta * exp_jm_phi_re;
2971 Ynm_r_im = pnm_costheta * exp_jm_phi_im;
2972
2973 Znm_t_re = snm_costheta * exp_jm_phi_re;
2974 Znm_t_im = snm_costheta * exp_jm_phi_im;
2975 Znm_p_re = -1.*unm_costheta * exp_jm_phi_im;
2976 Znm_p_im = unm_costheta * exp_jm_phi_re;
2977 switch (Ntype) {
2978 case 1: // Spherical Bessel function of the first kind j_n
2979 sph_bessel_n_ofkr_re = Spherical_j_n(n ,k0*r);
2980 sph_bessel_nminus1_ofkr_re = Spherical_j_n(n-1,k0*r);
2981 sph_bessel_n_ofkr_im = 0;
2982 sph_bessel_nminus1_ofkr_im = 0;
2983 break;
2984 case 2: // Spherical Bessel function of the second kind y_n
2985 sph_bessel_n_ofkr_re = Spherical_y_n(n ,k0*r);
2986 sph_bessel_nminus1_ofkr_re = Spherical_y_n(n-1,k0*r);
2987 sph_bessel_n_ofkr_im = 0;
2988 sph_bessel_nminus1_ofkr_im = 0;
2989 break;
2990 case 3: // Spherical Hankel function of the first kind h^1_n
2991 sph_bessel_n_ofkr_re = Spherical_j_n(n ,k0*r);
2992 sph_bessel_nminus1_ofkr_re = Spherical_j_n(n-1,k0*r);
2993 sph_bessel_n_ofkr_im = Spherical_y_n(n ,k0*r);
2994 sph_bessel_nminus1_ofkr_im = Spherical_y_n(n-1,k0*r);
2995 break;
2996 case 4: // Spherical Hankel function of the second kind h^2_n
2997 sph_bessel_n_ofkr_re = Spherical_j_n(n ,k0*r);
2998 sph_bessel_nminus1_ofkr_re = Spherical_j_n(n-1,k0*r);
2999 sph_bessel_n_ofkr_im =-Spherical_y_n(n ,k0*r);
3000 sph_bessel_nminus1_ofkr_im =-Spherical_y_n(n-1,k0*r);
3001 break;
3002 default:
3003 sph_bessel_n_ofkr_re = 0.;
3004 sph_bessel_nminus1_ofkr_re = 0.;
3005 sph_bessel_n_ofkr_im = 0.;
3006 sph_bessel_nminus1_ofkr_im = 0.;
3007 Message::Error("First argument for Nnm's should be 1,2,3 or 4");
3008 break;
3009 }
3010
3011 dRicatti_dx_ofkr_re = k0*r * (sph_bessel_nminus1_ofkr_re-
3012 (((double)n+1.)/(k0*r)) * sph_bessel_n_ofkr_re)
3013 + sph_bessel_n_ofkr_re ;
3014 dRicatti_dx_ofkr_im = k0*r * (sph_bessel_nminus1_ofkr_im-
3015 (((double)n+1.)/(k0*r)) * sph_bessel_n_ofkr_im)
3016 + sph_bessel_n_ofkr_im ;
3017
3018 Nnm_r_re = 1./(k0*r) * sqrt((((double)n)*((double)n+1.)))
3019 *(sph_bessel_n_ofkr_re*Ynm_r_re - sph_bessel_n_ofkr_im*Ynm_r_im);
3020 Nnm_r_im = 1./(k0*r) * sqrt( (((double)n)*((double)n+1.)) )
3021 *(sph_bessel_n_ofkr_re*Ynm_r_im + sph_bessel_n_ofkr_im*Ynm_r_re);
3022 Nnm_t_re = 1./(k0*r)
3023 *(dRicatti_dx_ofkr_re*Znm_t_re - dRicatti_dx_ofkr_im*Znm_t_im);
3024 Nnm_t_im = 1./(k0*r)
3025 *(dRicatti_dx_ofkr_re*Znm_t_im + dRicatti_dx_ofkr_im*Znm_t_re);
3026 Nnm_p_re = 1./(k0*r)
3027 *(dRicatti_dx_ofkr_re*Znm_p_re - dRicatti_dx_ofkr_im*Znm_p_im);
3028 Nnm_p_im = 1./(k0*r)
3029 *(dRicatti_dx_ofkr_re*Znm_p_im + dRicatti_dx_ofkr_im*Znm_p_re);
3030
3031 Nnm_x_re = sin(theta)*cos(phi)*Nnm_r_re
3032 + cos(theta)*cos(phi)*Nnm_t_re
3033 - sin(phi)*Nnm_p_re ;
3034 Nnm_y_re = sin(theta)*sin(phi)*Nnm_r_re
3035 + cos(theta)*sin(phi)*Nnm_t_re
3036 + cos(phi)*Nnm_p_re ;
3037 Nnm_z_re = cos(theta)*Nnm_r_re
3038 - sin(theta)*Nnm_t_re;
3039 Nnm_x_im = sin(theta)*cos(phi)*Nnm_r_im
3040 + cos(theta)*cos(phi)*Nnm_t_im
3041 - sin(phi)*Nnm_p_im ;
3042 Nnm_y_im = sin(theta)*sin(phi)*Nnm_r_im
3043 + cos(theta)*sin(phi)*Nnm_t_im
3044 + cos(phi)*Nnm_p_im ;
3045 Nnm_z_im = cos(theta)*Nnm_r_im
3046 - sin(theta)*Nnm_t_im;
3047
3048 V->Type = VECTOR;
3049 V->Val[0] = Nnm_x_re ; V->Val[MAX_DIM ] = Nnm_x_im ;
3050 V->Val[1] = Nnm_y_re ; V->Val[MAX_DIM+1] = Nnm_y_im ;
3051 V->Val[2] = Nnm_z_re ; V->Val[MAX_DIM+2] = Nnm_z_im ;
3052 }
3053
3054 /* ------------------------------------------------------------ */
3055 /* Dyadic Green's Function for homogeneous lossless media) */
3056 /* See e.g. Chew, Chap. 7, p. 375 */
3057 /* Basic usage : */
3058 /* Dyad_Green[] = DyadGreenHom[x_p,y_p,z_p,XYZ[],kb,siwt]; */
3059 /* ------------------------------------------------------------ */
F_DyadGreenHom(F_ARG)3060 void F_DyadGreenHom(F_ARG)
3061 {
3062 double x, y, z;
3063 double x_p, y_p, z_p;
3064 double kb, siwt;
3065 double normrr_p;
3066 double RR_xx,RR_xy,RR_xz,RR_yx,RR_yy,RR_yz,RR_zx,RR_zy,RR_zz;
3067 std::complex<double> Gsca,fact_diag,fact_glob;
3068 std::complex<double> G_xx,G_xy,G_xz,G_yx,G_yy,G_yz,G_zx,G_zy,G_zz;
3069 std::complex<double> I = std::complex<double>(0., 1.);
3070
3071 if( A->Type != SCALAR
3072 || (A+1)->Type != SCALAR
3073 || (A+2)->Type != SCALAR
3074 || (A+3)->Type != VECTOR
3075 || (A+4)->Type != SCALAR
3076 || (A+5)->Type != SCALAR)
3077 Message::Error("Non scalar argument(s) for the GreenHom");
3078 x_p = A->Val[0];
3079 y_p = (A+1)->Val[0];
3080 z_p = (A+2)->Val[0];
3081 x = (A+3)->Val[0];
3082 y = (A+3)->Val[1];
3083 z = (A+3)->Val[2];
3084 kb = (A+4)->Val[0];
3085 siwt = (A+5)->Val[0];
3086
3087 normrr_p = std::sqrt(std::pow((x-x_p),2)+std::pow((y-y_p),2)+std::pow((z-z_p),2));
3088
3089 Gsca = std::exp(-1.*siwt*I*kb*normrr_p)/(4.*M_PI*normrr_p);
3090 fact_diag = 1. + (I*kb*normrr_p-1.)/std::pow((kb*normrr_p),2);
3091 fact_glob = (3.-3.*I*kb*normrr_p-std::pow((kb*normrr_p),2))/
3092 std::pow((kb*std::pow(normrr_p,2)),2);
3093
3094 RR_xx = (x-x_p)*(x-x_p);
3095 RR_xy = (x-x_p)*(y-y_p);
3096 RR_xz = (x-x_p)*(z-z_p);
3097 RR_yx = (y-y_p)*(x-x_p);
3098 RR_yy = (y-y_p)*(y-y_p);
3099 RR_yz = (y-y_p)*(z-z_p);
3100 RR_zx = (z-z_p)*(x-x_p);
3101 RR_zy = (z-z_p)*(y-y_p);
3102 RR_zz = (z-z_p)*(z-z_p);
3103
3104 G_xx = Gsca*(fact_glob*RR_xx+fact_diag);
3105 G_xy = Gsca*(fact_glob*RR_xy);
3106 G_xz = Gsca*(fact_glob*RR_xz);
3107 G_yx = Gsca*(fact_glob*RR_yx);
3108 G_yy = Gsca*(fact_glob*RR_yy+fact_diag);
3109 G_yz = Gsca*(fact_glob*RR_yz);
3110 G_zx = Gsca*(fact_glob*RR_zx);
3111 G_zy = Gsca*(fact_glob*RR_zy);
3112 G_zz = Gsca*(fact_glob*RR_zz+fact_diag);
3113
3114 V->Type = TENSOR;
3115 V->Val[0] = G_xx.real() ; V->Val[MAX_DIM+0] = G_xx.imag() ;
3116 V->Val[1] = G_xy.real() ; V->Val[MAX_DIM+1] = G_xy.imag() ;
3117 V->Val[2] = G_xz.real() ; V->Val[MAX_DIM+2] = G_xz.imag() ;
3118 V->Val[3] = G_yx.real() ; V->Val[MAX_DIM+3] = G_yx.imag() ;
3119 V->Val[4] = G_yy.real() ; V->Val[MAX_DIM+4] = G_yy.imag() ;
3120 V->Val[5] = G_yz.real() ; V->Val[MAX_DIM+5] = G_yz.imag() ;
3121 V->Val[6] = G_zx.real() ; V->Val[MAX_DIM+6] = G_zx.imag() ;
3122 V->Val[7] = G_zy.real() ; V->Val[MAX_DIM+7] = G_zy.imag() ;
3123 V->Val[8] = G_zz.real() ; V->Val[MAX_DIM+8] = G_zz.imag() ;
3124 }
F_CurlDyadGreenHom(F_ARG)3125 void F_CurlDyadGreenHom(F_ARG)
3126 {
3127 double x, y, z;
3128 double x_p, y_p, z_p;
3129 double kb, siwt;
3130 double normrr_p;
3131 std::complex<double> Gsca,dx_Gsca,dy_Gsca,dz_Gsca;
3132 std::complex<double> curlG_xx,curlG_xy,curlG_xz,curlG_yx,curlG_yy,curlG_yz,curlG_zx,curlG_zy,curlG_zz;
3133 std::complex<double> I = std::complex<double>(0., 1.);
3134
3135 if( A->Type != SCALAR
3136 || (A+1)->Type != SCALAR
3137 || (A+2)->Type != SCALAR
3138 || (A+3)->Type != VECTOR
3139 || (A+4)->Type != SCALAR
3140 || (A+5)->Type != SCALAR)
3141 Message::Error("Non scalar argument(s) for the GreenHom");
3142 x_p = A->Val[0];
3143 y_p = (A+1)->Val[0];
3144 z_p = (A+2)->Val[0];
3145 x = (A+3)->Val[0];
3146 y = (A+3)->Val[1];
3147 z = (A+3)->Val[2];
3148 kb = (A+4)->Val[0];
3149 siwt = (A+5)->Val[0];
3150
3151 normrr_p = std::sqrt(std::pow((x-x_p),2)+std::pow((y-y_p),2)+std::pow((z-z_p),2));
3152
3153 Gsca = std::exp(-1.*siwt*I*kb*normrr_p)/(4.*M_PI*normrr_p);
3154
3155 dx_Gsca = (I*kb-1/normrr_p)*Gsca*(x-x_p);
3156 dy_Gsca = (I*kb-1/normrr_p)*Gsca*(y-y_p);
3157 dz_Gsca = (I*kb-1/normrr_p)*Gsca*(z-z_p);
3158
3159 curlG_xx = 0.;
3160 curlG_xy = -dz_Gsca;
3161 curlG_xz = dy_Gsca;
3162 curlG_yx = dz_Gsca;
3163 curlG_yy = 0.;
3164 curlG_yz = -dx_Gsca;
3165 curlG_zx = -dy_Gsca;
3166 curlG_zy = dx_Gsca;
3167 curlG_zz = 0.;
3168
3169 V->Type = TENSOR;
3170 V->Val[0] = curlG_xx.real() ; V->Val[MAX_DIM+0] = curlG_xx.imag() ;
3171 V->Val[1] = curlG_xy.real() ; V->Val[MAX_DIM+1] = curlG_xy.imag() ;
3172 V->Val[2] = curlG_xz.real() ; V->Val[MAX_DIM+2] = curlG_xz.imag() ;
3173 V->Val[3] = curlG_yx.real() ; V->Val[MAX_DIM+3] = curlG_yx.imag() ;
3174 V->Val[4] = curlG_yy.real() ; V->Val[MAX_DIM+4] = curlG_yy.imag() ;
3175 V->Val[5] = curlG_yz.real() ; V->Val[MAX_DIM+5] = curlG_yz.imag() ;
3176 V->Val[6] = curlG_zx.real() ; V->Val[MAX_DIM+6] = curlG_zx.imag() ;
3177 V->Val[7] = curlG_zy.real() ; V->Val[MAX_DIM+7] = curlG_zy.imag() ;
3178 V->Val[8] = curlG_zz.real() ; V->Val[MAX_DIM+8] = curlG_zz.imag() ;
3179 }
3180 #undef F_ARG
3181