"Fossies" - the Fresh Open Source Software Archive

Member "getdp-3.1.0-source/Functions/F_Analytic.cpp" (4 Apr 2019, 93072 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "F_Analytic.cpp" see the Fossies "Dox" file reference documentation.

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