From 72397e86f117b57fe223ffac12886a064e32e46a Mon Sep 17 00:00:00 2001 From: Maxim Yurkin Date: Sat, 10 Sep 2022 16:10:26 +0700 Subject: [PATCH] Lays the foundation for further improvement (and corrections) in calculation of Sommerfeld integrals - updates somnec.c according to the latest version of NEC2C (1.3) found on the Internet. This includes cosmetic changes, zero-initialization of several variable, moving declaration of several variables inside the cycles (where they are used), and removal of several code block that were never reached. The choice of tabs and indents is now uniform over the whole file. - there is a change in rom1() in somnec.c (removal of test for `nt`), which may potentially change the results, but existing tests does not indicate that (requires further study). - adds a lot of comments to somnec.c, mostly explaining how the integration contours and branch cuts for functions involving square roots are chosen in the complex plane. Also mentions several ideas for further improvements. - evlua() now accepts additional argument (mode) to enable testing different equivalent ways to compute the integral. The main ADDA code uses mode=0 (automatic choice) leading to no changes in operation. - creates header file somnec.h (to be able to use somnec uniformly from several places) - creates a standalone test module in somnec_test.c - the new files fully adhere to ADDA code style, while somnec.c is only half way there --- src/README.md | 3 +- src/interaction.c | 6 +- src/somnec.c | 1800 +++++++++++++++++++++++---------------------- src/somnec.h | 28 + src/somnec_test.c | 77 ++ 5 files changed, 1038 insertions(+), 876 deletions(-) create mode 100644 src/somnec.h create mode 100644 src/somnec_test.c diff --git a/src/README.md b/src/README.md index 2346ee4b..18b41ada 100644 --- a/src/README.md +++ b/src/README.md @@ -13,5 +13,6 @@ Main source of ADDA, including Makefiles and other scripts. C++ and Fortran sour * `iw_compile.bat` - batch script to compile ADDA with Intel compilers on Windows * `mt19937ar.c/h` – source and header files for Mersenne twister random generator * `oclkernels.cl` - OpenCL kernels -* `somnec.c` - routines to compute Sommerfeld integrals from NEC2C package +* `somnec.c/h` - source and header files for computation of Sommerfeld integrals (based on the NEC2C package) +* `somnec_test.c` - a standalone routine to test the calculation of Sommerfeld integrals * `updgithash.sh` - script to obtain git hash of the source and store it into generated `githash.h` \ No newline at end of file diff --git a/src/interaction.c b/src/interaction.c index 9b822158..ec521176 100644 --- a/src/interaction.c +++ b/src/interaction.c @@ -20,6 +20,7 @@ #include "igt_so.h" #include "io.h" #include "memory.h" +#include "somnec.h" #include "vars.h" // system headers #include // for DBL_EPSILON @@ -93,9 +94,6 @@ void propaespacelibreintadda_(const double *Rij,const double *ka,const double *g #endif // sinint.c void cisi(double x,double *ci,double *si); -// somnec.c -void som_init(complex double epscf); -void evlua(double zphIn,double rhoIn,complex double *erv, complex double *ezv,complex double *erh, complex double *eph); // this is used for debugging, should be empty define, when not required #define PRINT_GVAL /*printf("%s: %d,%d,%d: %g%+gi, %g%+gi, %g%+gi,\n%g%+gi, %g%+gi, %g%+gi\n",__func__,i,j,k,\ @@ -705,7 +703,7 @@ static inline void SingleSomIntegral(double rho,const double z,doublecomplex val const double isc=pow(scale,3); // this is subject to under/overflow if (rho==0) rho=z*0.00000001; // a hack to overcome the poor precision of somnec for rho=0; - evlua(z*scale,rho*scale,vals,vals+1,vals+2,vals+3); + evlua(z*scale,rho*scale,vals,vals+1,vals+2,vals+3,0); vals[0]*=isc; vals[1]*=isc; vals[2]*=isc; diff --git a/src/somnec.c b/src/somnec.c index 4339179b..4b8b60c3 100644 --- a/src/somnec.c +++ b/src/somnec.c @@ -1,13 +1,12 @@ -/* last change: pgm 8 nov 2000 1:04 pm */ -/* program somnec(input,output,tape21) */ - -/* Original code is NEC2C (in public domain) obtained from http://www.si-list.net/swindex.html, - * http://sharon.esrac.ele.tue.nl/users/pe1rxq/nec2c.rxq/nec2c.rxq-0.2.tar.gz +/* Original code (somnec.c) is part of NEC2C v.1.3 (in public domain) obtained + * from https://www.qsl.net/5b4az/pkg/nec2/nec2c/nec2c-1.3.tar.bz2 * This code is incorporated into ADDA with several changes, which are explicitly indicated by comments. Those include: * - no generation of interpolation grid, only single run - * - numerical precision was changed to double * - conjugation (that was in place to couple with other parts of nec2 code) was removed * - a few cosmetic changes to remove compiler warnings (with -pedantic flag) + * - replacement of 20 to MAXH in couple of places + * + * Further updates */ /* TODO: Systematic accuracy study of this code is required. At least 7 digits of precision are desired (for test runs) @@ -16,24 +15,25 @@ */ // ADDA: the cycling to generate interpolation grid was removed, now it is a single-run routine -/* program to generate nec interpolation grids for fields due to */ -/* ground. field components are computed by numerical evaluation */ -/* of modified sommerfeld integrals. */ -// ADDA: double precision is now used throughout -/* somnec2d is a double precision version of somnec for use with */ -/* nec2d. an alternate version (somnec2sd) is also provided in which */ -/* computation is in single precision but the output file is written */ -/* in double precision for use with nec2d. somnec2sd runs about twic */ -/* as fast as the full double precision somnec2d. the difference */ -/* between nec2d results using a for021 file from this code rather */ -/* than from somnec2sd was insignficant in the cases tested. */ - -/* changes made by j bergervoet, 31-5-95: */ -/* parameter 0. --> 0.d0 in calling of routine test */ -/* status of output files set to 'unknown' */ +/* program to generate nec interpolation grids for fields due to + ground. field components are computed by numerical evaluation + of modified sommerfeld integrals. + + somnec2d is a double precision version of somnec for use with + nec2d. an alternate version (somnec2sd) is also provided in which + computation is in single precision but the output file is written + in double precision for use with nec2d. somnec2sd runs about twice + as fast as the full double precision somnec2d. the difference + between nec2d results using a for021 file from this code rather + than from somnec2sd was insignificant in the cases tested. + + changes made by j bergervoet, 31-5-95: + parameter 0. --> 0.d0 in calling of routine test + status of output files set to 'unknown' */ + +#include "somnec.h" // corresponding header // ADDA: The following definitions are extracted from the file nec2c.h (only those that are used here) -#include #include #include #include @@ -46,23 +46,23 @@ #endif /* common constants */ -#define PI 3.141592654 -#define TP 6.283185308 -#define PTP .6283185308 -#define PI10 31.41592654 -#define GAMMA .5772156649 -#define C1 -.02457850915 -#define C2 .3674669052 -#define C3 .7978845608 -#define P10 .0703125 -#define P20 .1121520996 -#define Q10 .125 -#define Q20 .0732421875 -#define P11 .1171875 -#define P21 .1441955566 -#define Q11 .375 -#define Q21 .1025390625 -#define POF .7853981635 +#define PI 3.141592654 +#define TP 6.283185308 +#define PTP 0.6283185308 +#define PI10 31.41592654 +#define GAMMA 0.5772156649 +#define C1 -0.02457850915 +#define C2 0.3674669052 +#define C3 0.7978845608 +#define P10 0.0703125 +#define P20 0.1121520996 +#define Q10 0.125 +#define Q20 0.0732421875 +#define P11 0.1171875 +#define P21 0.1441955566 +#define Q11 0.375 +#define Q21 0.1025390625 +#define POF 0.7853981635 /* MAXH=100 seems sufficient to obtain convergence in all cases; larger values doesn't seem to improve the accuracy * even when smaller CRIT is used. Probably other approximations are employed somewhere. * However, larger MAXH is required for distances larger than several wavelengths to reach specified CRIT. In practice, @@ -70,28 +70,27 @@ * reflection term. Thus, we increased MAXH to 500 for now, but final solution is left for * https://github.com/adda-team/adda/issues/176 * If you get convergence errors, increase MAXH further + * Also, tsmag below may affect the accuracy. Currently, the relative errors due to tsmag is <= 3e-7 + * TODO: increase accuracy of constants above to double */ -#define MAXH 500 -#define CRIT 1.0E-4 -#define NM 131072 -#define NTS 4 +#define MAXH 500 +#define CRIT 1.0E-4 +#define NM 131072 +#define NTS 4 #define cmplx(r, i) ((r)+(i)*I) -void som_init(complex double epscf); -static void bessel(complex double z, complex double *j0, complex double *j0p); -void evlua(double zphIn,double rhoIn,complex double *erv, complex double *ezv, - complex double *erh, complex double *eph); -static void gshank(complex double start, complex double dela, complex double *sum, +static void bessel(complex double z, complex double *j0, complex double *j0p); +static void gshank(complex double start, complex double dela, complex double *sum, int nans, complex double *seed, int ibk, complex double bk, complex double delb); -static void hankel(complex double z, complex double *h0, complex double *h0p); -static void lambda(double t, complex double *xlam, complex double *dxlam); -static void rom1(int n, complex double *sum, int nx); -static void saoa( double t, complex double *ans); +static void hankel(complex double z, complex double *h0, complex double *h0p); +static void lambda(double t, complex double *xlam, complex double *dxlam); +static void rom1(int n, complex double *sum, int nx); +static void saoa( double t, complex double *ans); -static void test(double f1r, double f2r, double *tr, double f1i, +static void test(double f1r, double f2r, double *tr, double f1i, double f2i, double *ti, double dmin); -static void abort_on_error(int why); +static void abort_on_error(int why); /* common /evlcom/ */ static int jh; @@ -114,274 +113,318 @@ static inline double cAbs2(const complex double z) void som_init(complex double epscf) { - complex double erv, ezv; - - ck2=TP; - ck2sq=ck2*ck2; - - /* sommerfeld integral evaluation uses exp(-jwt) */ - // ADDA: removed conjugation of epscf and in evaluation of integrals below - - ck1sq=ck2sq*epscf; - ck1=csqrt(ck1sq); - ck1r=creal(ck1); - tkmag=100.*cabs(ck1); - tsmag=100.*cAbs2(ck1); - cksm=ck2sq/(ck1sq+ck2sq); - ct1=.5*(ck1sq-ck2sq); - erv=ck1sq*ck1sq; - ezv=ck2sq*ck2sq; - ct2=.125*(erv-ezv); - erv *= ck1sq; - ezv *= ck2sq; - ct3=.0625*(erv-ezv); - - return; + complex double erv, ezv; + + ck2=TP; + ck2sq=ck2*ck2; + + /* sommerfeld integral evaluation uses exp(-jwt) */ + // ADDA: removed conjugation of epscf and in evaluation of integrals below + + ck1sq=ck2sq*epscf; + ck1=csqrt(ck1sq); + ck1r=creal(ck1); + tkmag=100.*cabs(ck1); + // TODO: the largest of |ck1|, |ck2| should be used to define tsmag + tsmag=100.*cAbs2(ck1); + cksm=ck2sq/(ck1sq+ck2sq); + ct1=.5*(ck1sq-ck2sq); + erv=ck1sq*ck1sq; + ezv=ck2sq*ck2sq; + ct2=.125*(erv-ezv); + erv *= ck1sq; + ezv *= ck2sq; + ct3=.0625*(erv-ezv); + + return; } /*-----------------------------------------------------------------------*/ /* bessel evaluates the zero-order bessel function */ /* and its derivative for complex argument z. */ -static void bessel(complex double z, complex double *j0, complex double *j0p ) +static void bessel( complex double z, complex double *j0, complex double *j0p ) { - int k, i, ib, iz, miz; - static int m[101], init = FALSE; - static double a1[25], a2[25]; - double tst, zms; - complex double p0z, p1z, q0z, q1z, zi, zi2, zk, cz, sz, j0x, j0px; - - /* initialization of constants */ - if( ! init ) - { - for( k = 1; k <= 25; k++ ) - { - i = k-1; - a1[i]=-.25/(k*k); - a2[i]=1.0/(k+1.0); - } - - for( i = 1; i <= 101; i++ ) - { - tst=1.0; - for( k = 0; k < 24; k++ ) - { - init = k; - tst *= -i*a1[k]; - if( tst < 1.0e-6 ) - break; - } - - m[i-1] = init+1; - } /* for( i = 1; i<= 101; i++ ) */ - - init = TRUE; - } /* if(init == 0) */ - - zms=cAbs2(z); - if(zms <= 1.e-12) - { - *j0=1; - *j0p=-.5*z; - return; - } - - ib=0; - if(zms <= 37.21) - { - if(zms > 36.) - ib=1; - - /* series expansion */ - iz=zms; - miz=m[iz]; - *j0=1; - *j0p=*j0; - zk=*j0; - zi=z*z; - - for( k = 0; k < miz; k++ ) - { - zk *= a1[k]*zi; - *j0 += zk; - *j0p += a2[k]*zk; - } - *j0p *= -.5*z; - - if(ib == 0) - return; - - j0x=*j0; - j0px=*j0p; - } - - /* asymptotic expansion */ - zi=1./z; - zi2=zi*zi; - p0z=1.+(P20*zi2-P10)*zi2; - p1z=1.+(P11-P21*zi2)*zi2; - q0z=(Q20*zi2-Q10)*zi; - q1z=(Q11-Q21*zi2)*zi; - zk=cexp(I*(z-POF)); - zi2=1./zk; - cz=.5*(zk+zi2); - sz=I*.5*(zi2-zk); - zk=C3*csqrt(zi); - *j0=zk*(p0z*cz-q0z*sz); - *j0p=-zk*(p1z*sz+q1z*cz); - - if(ib == 0) - return; - - zms=cos((sqrt(zms)-6.)*PI10); - *j0=.5*(j0x*(1.+zms)+ *j0*(1.-zms)); - *j0p=.5*(j0px*(1.+zms)+ *j0p*(1.-zms)); - - return; + int k, ib; + static int m[101], init = FALSE; + static double a1[25], a2[25]; + double zms; + complex double p0z, p1z, q0z, q1z, zi, zi2, zk, cz, sz, j0x=0, j0px=0; + + /* initialization of constants */ + if( ! init ) + { + int i; + for( k = 1; k <= 25; k++ ) + { + i = k-1; + a1[i]=-.25/(k*k); + a2[i]=1.0/(k+1.0); + } + + for( i = 1; i <= 101; i++ ) + { + double tst = 1.0; + for( k = 0; k < 24; k++ ) + { + init = k; + tst *= -i*a1[k]; + if( tst < 1.0e-6 ) + break; + } + + m[i-1] = init+1; + } /* for( i = 1; i<= 101; i++ ) */ + + init = TRUE; + } /* if(init == 0) */ + + zms=cAbs2(z); + if(zms <= 1.e-12) + { + *j0=1; + *j0p=-.5*z; + return; + } + + ib=0; + if(zms <= 37.21) + { + if(zms > 36.) + ib=1; + + /* series expansion */ + int iz=(int)zms; + int miz=m[iz]; + *j0=1; + *j0p=*j0; + zk=*j0; + zi=z*z; + + for( k = 0; k < miz; k++ ) + { + zk *= a1[k]*zi; + *j0 += zk; + *j0p += a2[k]*zk; + } + *j0p *= -.5*z; + + if(ib == 0) + return; + + j0x=*j0; + j0px=*j0p; + } + + /* asymptotic expansion */ + zi=1./z; + zi2=zi*zi; + p0z=1.+(P20*zi2-P10)*zi2; + p1z=1.+(P11-P21*zi2)*zi2; + q0z=(Q20*zi2-Q10)*zi; + q1z=(Q11-Q21*zi2)*zi; + zk=cexp(I*(z-POF)); + zi2=1./zk; + cz=.5*(zk+zi2); + sz=I*.5*(zi2-zk); + zk=C3*csqrt(zi); + *j0=zk*(p0z*cz-q0z*sz); + *j0p=-zk*(p1z*sz+q1z*cz); + + if(ib == 0) + return; + + zms=cos((sqrt(zms)-6.)*PI10); + *j0=.5*(j0x*(1.+zms)+ *j0*(1.-zms)); + *j0p=.5*(j0px*(1.+zms)+ *j0p*(1.-zms)); + + return; } /*-----------------------------------------------------------------------*/ /* evlua controls the integration contour in the complex */ /* lambda plane for evaluation of the sommerfeld integrals */ +/* the logic for choosing integration contours follows Section 4.2 "Numerical Evaluation of the Sommerfeld Integrals" of + * http://users.tpg.com.au/micksw012/files/gnec-theory.pdf (gNEC Theory v. 0.9.15, 02.02.2018) with the only change of + * conjugating all complex variables due to different time-harmonic convention. In the following "on the figure" refers to + * contour descriptions in this section. + * Improvements upon this logic (if any) will be noted explicitly. + */ void evlua(double zphIn,double rhoIn, complex double *erv, complex double *ezv, - complex double *erh, complex double *eph ) + complex double *erh, complex double *eph, int mode) +/* mode determines the choice of contour: automatic (0), Bessel (1), shorter Hankel (2), longer Hankel (3) + * mode=2 is effectively used is some other value is specified (TODO: add a meaningful test for that) + */ { - int i, jump; - static double del, slope, rmis; - static complex double cp1, cp2, cp3, bk, delta, delta2, sum[6], ans[6]; - - zph=zphIn; - rho=rhoIn; - - del=zph; - if( rho > del ) - del=rho; - - if(zph >= 2.*rho) - { - /* bessel function form of sommerfeld integrals */ - jh=0; - a=0; - del=1./del; - - if( del > tkmag) - { - b=cmplx(.1*tkmag,-.1*tkmag); - rom1(6,sum,2); - a=b; - b=cmplx(del,-del); - rom1 (6,ans,2); - for( i = 0; i < 6; i++ ) - sum[i] += ans[i]; - } - else - { - b=cmplx(del,-del); - rom1(6,sum,2); - } - - delta=PTP*del; - gshank(b,delta,ans,6,sum,0,b,b); - ans[5] *= ck1; - - /* ADDA: conjugate was removed */ - *erv=ck1sq*ans[2]; - *ezv=ck1sq*(ans[1]+ck2sq*ans[4]); - *erh=ck2sq*(ans[0]+ans[5]); - *eph=-ck2sq*(ans[3]+ans[5]); - - return; - - } /* if(zph >= 2.*rho) */ - - /* hankel function form of sommerfeld integrals */ - jh=1; - cp1=cmplx(0.0,.4*ck2); - cp2=cmplx(.6*ck2,-.2*ck2); - cp3=cmplx(1.02*ck2,-.2*ck2); - a=cp1; - b=cp2; - rom1(6,sum,2); - a=cp2; - b=cp3; - rom1(6,ans,2); - - for( i = 0; i < 6; i++ ) - sum[i]=-(sum[i]+ans[i]); - - /* path from imaginary axis to -infinity */ - if(zph > .001*rho) - slope=rho/zph; - else - slope=1000.; - - del=PTP/del; - delta=cmplx(-1.0,slope)*del/sqrt(1.+slope*slope); - delta2=-conj(delta); - gshank(cp1,delta,ans,6,sum,0,bk,bk); - rmis=rho*(creal(ck1)-ck2); - - jump = FALSE; - if( (rmis >= 2.*ck2) && (rho >= 1.e-10) ) - { - if(zph >= 1.e-10) - { - bk=cmplx(-zph,rho)*(ck1-cp3); - rmis=-creal(bk)/fabs(cimag(bk)); - if(rmis > 4.*rho/zph) - jump = TRUE; - } - - if( ! jump ) - { - /* integrate up between branch cuts, then to + infinity */ - cp1=ck1-cmplx(.1,.2); - cp2=cp1+.2; - bk=cmplx(0.,del); - gshank(cp1,bk,sum,6,ans,0,bk,bk); - a=cp1; - b=cp2; - rom1(6,ans,1); - for( i = 0; i < 6; i++ ) - ans[i] -= sum[i]; - - gshank(cp3,bk,sum,6,ans,0,bk,bk); - gshank(cp2,delta2,ans,6,sum,0,bk,bk); - } - - jump = TRUE; - - } /* if( (rmis >= 2.*ck2) || (rho >= 1.e-10) ) */ - else - jump = FALSE; - - if( ! jump ) - { - /* integrate below branch points, then to + infinity */ - for( i = 0; i < 6; i++ ) - sum[i]=-ans[i]; - - rmis=creal(ck1)*1.01; - if( (ck2+1.) > rmis ) - rmis=ck2+1.; - - bk=cmplx(rmis,.99*cimag(ck1)); - delta=bk-cp3; - delta *= del/cabs(delta); - gshank(cp3,delta,ans,6,sum,1,bk,delta2); - - } /* if( ! jump ) */ - - ans[5] *= ck1; - - /* ADDA: conjugate was removed */ - *erv=ck1sq*ans[2]; - *ezv=ck1sq*(ans[1]+ck2sq*ans[4]); - *erh=ck2sq*(ans[0]+ans[5]); - *eph=-ck2sq*(ans[3]+ans[5]); - - return; + int i, jump; + static double del, slope, rmis; + static complex double cp1, cp2, cp3, bk, delta, delta2, sum[6], ans[6]; + + zph=zphIn; + rho=rhoIn; + + del=zph; + if( rho > del ) + del=rho; + + if (mode==1 || (mode==0 && zph >= 2.*rho)) + { + /* Bessel-function form of Sommerfeld integrals */ + jh=0; + a=0; + del=1./del; + + if( del > tkmag) + { + /* when both rho and Z are small, the contour from zero to inflection point is divided into two parts, + * probably to better control the accuracy under the large variation of the integrand over the contour + * (at least, of the J0'). This is not explained in the docs, and it is not clear, why the k1 (not k2) is + * used for the boundary value tkmag. Still, any changes will require extensive testing to justify. + */ + b=cmplx(.1*tkmag,-.1*tkmag); + rom1(6,sum,2); // from zero to intermediate point + a=b; + b=cmplx(del,-del); + rom1 (6,ans,2); // from intermediate point to the inflection one on the figure + for( i = 0; i < 6; i++ ) + sum[i] += ans[i]; + } + else + { + b=cmplx(del,-del); + rom1(6,sum,2); // from zero to inflection point on the figure + } + + delta=PTP*del; + gshank(b,delta,ans,6,sum,0,b,b); // horizontal line to infinity on the figure + ans[5] *= ck1; + + /* ADDA: conjugate was removed */ + *erv=ck1sq*ans[2]; + *ezv=ck1sq*(ans[1]+ck2sq*ans[4]); + *erh=ck2sq*(ans[0]+ans[5]); + *eph=-ck2sq*(ans[3]+ans[5]); + + return; + + } /* if(zph >= 2.*rho) */ + + /* Hankel-function form of Sommerfeld integrals, based on H0^(1)(rho*xl) */ + jh=1; + cp1=cmplx(0.0,.4*ck2); + cp2=cmplx(.6*ck2,-.2*ck2); + cp3=cmplx(1.02*ck2,-.2*ck2); + a=cp1; + b=cp2; + rom1(6,sum,2); // from a to b on the figure + a=cp2; + b=cp3; + rom1(6,ans,2); // from b to c on the figure + + /* minus signs are used somewhat complicatedly in the following, motivated by the fact that some integrals are + * computed from infinity to a point, i.e., in reverse direction to that assumed by gshank + */ + for( i = 0; i < 6; i++ ) + sum[i]=-(sum[i]+ans[i]); + + /* path from imaginary axis to -infinity */ + /* Overall, the paths to infinity are chosen so that the exp(-gamma2*Z)*H0(1)(rho*xl) asymptotically decays as + * exp(-a*|xl|) with real a>0 (for large |xl|). Below, the difference between delta and delta2 is due to their usage + * in the second and first quadrant, respectively, where gamma2 ~ -z and z. + */ + if(zph > .001*rho) + slope=rho/zph; + // this avoids very large numbers below, but convergence is guaranteed to be excellent along this direction + else + slope=1000.; + + del=PTP/del; + delta=cmplx(-1.0,slope)*del/sqrt(1.+slope*slope); + delta2=-conj(delta); + //TODO: replace bk with 0 in last two arguments to gshank, whenever they are not used (i.e., when ibk=0) + gshank(cp1,delta,ans,6,sum,0,bk,bk); // from infinity to a on the figure + // by this point ans hold minus integral from infinity to c through a and b + rmis=rho*(creal(ck1)-ck2); + + jump = FALSE; + if (mode==3 || (mode==0 && rmis >= 2.*ck2 && rho >= 1.e-10)) + { + if(zph >= 1.e-10) + { + /* Here we choose between two alternatives: (1) integrate directly between two branch points (c to d on + * figure) or (2) from c to infinity and then back to d along two vertical lines near the corresponding + * branch cuts. The choice of one over another depends on the |arg(a)| in asymptotic expression exp(-a*|xl|) + * for the integrand. We want to minimize this argument. + * In these cases, |tg[arg(a)]| is approximately tg[arg((Z-i*rho)*(d-c))] and Z/rho, respectively. Inverse + * of these quantities are compared in the following. + * Factor 4 seems to be an empirical one to balance the two options (either one finite integral or two + * infinite ones). However, this specific values is not discussed in any docs, moreover, based on the + * balance arguments, it is more logical for it to be 1/4 or 1/2. The only possible argument in favor of 4 + * is that a finite integral does not necessarily reaches the asymptotic behavior, and thus is more + * complicated for adaptive routines. + * Anyway, changing this factor can only be done after careful testing of computational times to reach the + * prescribed accuracy. + */ + bk=cmplx(-zph,rho)*(ck1-cp3); + rmis=-creal(bk)/fabs(cimag(bk)); + // TODO: better to use inverse comparison, since Z can be 0 + // TODO: the following should be corrected not to skip the case below (from c to d) + if(rmis > 4.*rho/zph) + jump = TRUE; + } + + if( ! jump ) + { + /* integrate up to infinity and back between branch cuts, then to + infinity on the right side */ + cp1=ck1-cmplx(.1,.2); + cp2=cp1+.2; + bk=cmplx(0.,del); + gshank(cp1,bk,sum,6,ans,0,bk,bk); // from e to infinity on the figure + a=cp1; + b=cp2; + rom1(6,ans,1); + // invert sign from previous parts of contour + for( i = 0; i < 6; i++ ) + ans[i] -= sum[i]; + + gshank(cp3,bk,sum,6,ans,0,bk,bk); // from c to infinity on the figure + gshank(cp2,delta2,ans,6,sum,0,bk,bk); // from f to infinity on the figure + } + + jump = TRUE; + + } /* if( (rmis >= 2.*ck2) && (rho >= 1.e-10) ) */ + else + jump = FALSE; + + if( ! jump ) // this is always enabled when mode==2 (or any other unrecognized value) + { + /* integrate below branch points, then to + infinity */ + // invert sign from previous parts of contour + for( i = 0; i < 6; i++ ) + sum[i]=-ans[i]; + + rmis=creal(ck1)*1.01; + if( (ck2+1.) > rmis ) + rmis=ck2+1.; + + bk=cmplx(rmis,.99*cimag(ck1)); + delta=bk-cp3; + delta *= del/cabs(delta); + gshank(cp3,delta,ans,6,sum,1,bk,delta2); // from c to d and then to infinity on the figure + + } /* if( ! jump ) */ + + ans[5] *= ck1; + + /* ADDA: conjugate was removed */ + *erv=ck1sq*ans[2]; + *ezv=ck1sq*(ans[1]+ck2sq*ans[4]); + *erh=ck2sq*(ans[0]+ans[5]); + *eph=-ck2sq*(ans[3]+ans[5]); + + return; } /*-----------------------------------------------------------------------*/ @@ -391,170 +434,151 @@ void evlua(double zphIn,double rhoIn, complex double *erv, complex double *ezv, /* the step increment may be changed from dela to delb. shank's */ /* algorithm to accelerate convergence of a slowly converging series */ /* is used */ -static void gshank( complex double start, complex double dela, complex double *sum, - int nans, complex double *seed, int ibk, complex double bk, complex double delb ) +static void gshank( complex double start, complex double dela, + complex double *sum, int nans, complex double *seed, + int ibk, complex double bk, complex double delb ) { - int ibx, j, i, jm, intx, inx, brk=0; - static double rbk, amg, den, denm; - static complex double a1, a2, as1, as2, del, aa; - static complex double q1[6][MAXH], q2[6][MAXH], ans1[6], ans2[6]; - - rbk=creal(bk); - del=dela; - if(ibk == 0) - ibx=1; - else - ibx=0; - - for( i = 0; i < nans; i++ ) - ans2[i]=seed[i]; - - b=start; - for( intx = 1; intx <= MAXH; intx++ ) - { - inx=intx-1; - a=b; - b += del; - - if( (ibx == 0) && (creal(b) >= rbk) ) - { - /* hit break point. reset seed and start over. */ - ibx=1; - b=bk; - del=delb; - rom1(nans,sum,2); - if( ibx != 2 ) - { - for( i = 0; i < nans; i++ ) - ans2[i] += sum[i]; - intx = 0; - continue; - } - - for( i = 0; i < nans; i++ ) - ans2[i]=ans1[i]+sum[i]; - intx = 0; - continue; - - } /* if( (ibx == 0) && (creal(b) >= rbk) ) */ - - rom1(nans,sum,2); - for( i = 0; i < nans; i++ ) - ans1[i] = ans2[i]+sum[i]; - a=b; - b += del; - - if( (ibx == 0) && (creal(b) >= rbk) ) - { - /* hit break point. reset seed and start over. */ - ibx=2; - b=bk; - del=delb; - rom1(nans,sum,2); - if( ibx != 2 ) - { - for( i = 0; i < nans; i++ ) - ans2[i] += sum[i]; - intx = 0; - continue; - } - - for( i = 0; i < nans; i++ ) - ans2[i] = ans1[i]+sum[i]; - intx = 0; - continue; - - } /* if( (ibx == 0) && (creal(b) >= rbk) ) */ - - rom1(nans,sum,2); - for( i = 0; i < nans; i++ ) - ans2[i]=ans1[i]+sum[i]; - - den=0.; - for( i = 0; i < nans; i++ ) - { - as1=ans1[i]; - as2=ans2[i]; - - if(intx >= 2) - { - for( j = 1; j < intx; j++ ) - { - jm=j-1; - aa=q2[i][jm]; - a1=q1[i][jm]+as1-2.*aa; - - if( (creal(a1) != 0.) || (cimag(a1) != 0.) ) - { - a2=aa-q1[i][jm]; - a1=q1[i][jm]-a2*a2/a1; - } - else - a1=q1[i][jm]; - - a2=aa+as2-2.*as1; - if( (creal(a2) != 0.) || (cimag(a2) != 0.) ) - a2=aa-(as1-aa)*(as1-aa)/a2; - else - a2=aa; - - q1[i][jm]=as1; - q2[i][jm]=as2; - as1=a1; - as2=a2; - - } /* for( j = 1; i < intx; i++ ) */ - - } /* if(intx >= 2) */ - - q1[i][intx-1]=as1; - q2[i][intx-1]=as2; - amg=fabs(creal(as2))+fabs(cimag(as2)); - if(amg > den) - den=amg; - - } /* for( i = 0; i < nans; i++ ) */ - - denm=1.e-3*den*CRIT; - jm=intx-3; - if(jm < 1) - jm=1; - - for( j = jm-1; j < intx; j++ ) - { - brk = FALSE; - for( i = 0; i < nans; i++ ) - { - a1=q2[i][j]; - den=(fabs(creal(a1))+fabs(cimag(a1)))*CRIT; - if(den < denm) - den=denm; - a1=q1[i][j]-a1; - amg=fabs(creal(a1)+fabs(cimag(a1))); - if(amg > den) - { - brk = TRUE; - break; - } - - } /* for( i = 0; i < nans; i++ ) */ - - if( brk ) break; - - } /* for( j = jm-1; j < intx; j++ ) */ - - if( ! brk ) - { - for( i = 0; i < nans; i++ ) - sum[i]=.5*(q1[i][inx]+q2[i][inx]); - return; - } + int ibx, j, i, jm, intx, inx, brk=0; + static double rbk, amg, den, denm; + static complex double a1, a2, as1, as2, del, aa; + static complex double q1[6][MAXH], q2[6][MAXH], ans1[6], ans2[6]; + + rbk=creal(bk); + del=dela; + if(ibk == 0) + ibx=1; + else + ibx=0; - } /* for( intx = 1; intx <= maxh; intx++ ) */ + for( i = 0; i < nans; i++ ) + ans2[i]=seed[i]; + b=start; - /* No convergence */ - printf("z=%g, rho=%g\n",zph,rho); - abort_on_error(-6); + for( intx = 1; intx <= MAXH; intx++ ) + { + inx=intx-1; + a=b; + b += del; + + if( (ibx == 0) && (creal(b) >= rbk) ) + { + /* hit break point. reset seed and start over. */ + ibx=1; + b=bk; + del=delb; + rom1(nans,sum,2); + for( i = 0; i < nans; i++ ) + ans2[i] += sum[i]; + intx = 0; + continue; + } /* if( (ibx == 0) && (creal(b) >= rbk) ) */ + + rom1(nans,sum,2); + for( i = 0; i < nans; i++ ) + ans1[i] = ans2[i]+sum[i]; + a=b; + b += del; + if( (ibx == 0) && (creal(b) >= rbk) ) + { + /* hit break point. reset seed and start over. */ + ibx=2; + b=bk; + del=delb; + rom1(nans,sum,2); + for( i = 0; i < nans; i++ ) + ans2[i] = ans1[i]+sum[i]; + intx = 0; + continue; + } /* if( (ibx == 0) && (creal(b) >= rbk) ) */ + + rom1(nans,sum,2); + for( i = 0; i < nans; i++ ) + ans2[i]=ans1[i]+sum[i]; + + den=0.; + for( i = 0; i < nans; i++ ) + { + as1=ans1[i]; + as2=ans2[i]; + + if(intx >= 2) + { + for( j = 1; j < intx; j++ ) + { + jm=j-1; + aa=q2[i][jm]; + a1=q1[i][jm]+as1-2.*aa; + + if( (creal(a1) != 0.) || (cimag(a1) != 0.) ) + { + a2=aa-q1[i][jm]; + a1=q1[i][jm]-a2*a2/a1; + } + else + a1=q1[i][jm]; + + a2=aa+as2-2.*as1; + if( (creal(a2) != 0.) || (cimag(a2) != 0.) ) + a2=aa-(as1-aa)*(as1-aa)/a2; + else + a2=aa; + + q1[i][jm]=as1; + q2[i][jm]=as2; + as1=a1; + as2=a2; + } /* for( j = 1; i < intx; i++ ) */ + + } /* if(intx >= 2) */ + + q1[i][intx-1]=as1; + q2[i][intx-1]=as2; + amg=fabs(creal(as2))+fabs(cimag(as2)); + if(amg > den) + den=amg; + + } /* for( i = 0; i < nans; i++ ) */ + + denm=1.e-3*den*CRIT; + jm=intx-3; + if(jm < 1) + jm=1; + + for( j = jm-1; j < intx; j++ ) + { + brk = FALSE; + for( i = 0; i < nans; i++ ) + { + a1=q2[i][j]; + den=(fabs(creal(a1))+fabs(cimag(a1)))*CRIT; + if(den < denm) + den=denm; + a1=q1[i][j]-a1; + amg=fabs(creal(a1)+fabs(cimag(a1))); + if(amg > den) + { + brk = TRUE; + break; + } + + } /* for( i = 0; i < nans; i++ ) */ + + if( brk ) break; + + } /* for( j = jm-1; j < intx; j++ ) */ + + if( ! brk ) + { + for( i = 0; i < nans; i++ ) + sum[i]=.5*(q1[i][inx]+q2[i][inx]); + return; + } + + } /* for( intx = 1; intx <= maxh; intx++ ) */ + + /* No convergence */ + printf("z=%g, rho=%g\n",zph,rho); + abort_on_error(-6); } /*-----------------------------------------------------------------------*/ @@ -563,105 +587,106 @@ static void gshank( complex double start, complex double dela, complex double *s /* order zero, and its derivative for complex argument z */ static void hankel( complex double z, complex double *h0, complex double *h0p ) { - int i, k, ib, iz, miz; - static int m[101], init = FALSE; - static double a1[25], a2[25], a3[25], a4[25], psi, tst, zms; - complex double clogz, j0, j0p, p0z, p1z, q0z, q1z, y0, y0p, zi, zi2, zk; - - /* initialization of constants */ - if( ! init ) - { - psi=-GAMMA; - for( k = 1; k <= 25; k++ ) - { - i = k-1; - a1[i]=-.25/(k*k); - a2[i]=1.0/(k+1.0); - psi += 1.0/k; - a3[i]=psi+psi; - a4[i]=(psi+psi+1.0/(k+1.0))/(k+1.0); - } - - for( i = 1; i <= 101; i++ ) - { - tst=1.0; - for( k = 0; k < 24; k++ ) - { - init = k; - tst *= -i*a1[k]; - if(tst*a3[k] < 1.e-6) - break; - } - m[i-1]=init+1; - } - - init = TRUE; - - } /* if( ! init ) */ - - zms=cAbs2(z); - if(zms == 0.) - abort_on_error(-7); - - ib=0; - if(zms <= 16.81) - { - if(zms > 16.) - ib=1; - - /* series expansion */ - iz=zms; - miz=m[iz]; - j0=1; - j0p=j0; - y0=0; - y0p=y0; - zk=j0; - zi=z*z; - - for( k = 0; k < miz; k++ ) - { - zk *= a1[k]*zi; - j0 += zk; - j0p += a2[k]*zk; - y0 += a3[k]*zk; - y0p += a4[k]*zk; - } - - j0p *= -.5*z; - clogz=clog(.5*z); - y0=(2.*j0*clogz-y0)/PI+C2; - y0p=(2./z+2.*j0p*clogz+.5*y0p*z)/PI+C1*z; - *h0=j0+I*y0; - *h0p=j0p+I*y0p; - - if(ib == 0) - return; - - y0=*h0; - y0p=*h0p; - - } /* if(zms <= 16.81) */ - - /* asymptotic expansion */ - zi=1./z; - zi2=zi*zi; - p0z=1.+(P20*zi2-P10)*zi2; - p1z=1.+(P11-P21*zi2)*zi2; - q0z=(Q20*zi2-Q10)*zi; - q1z=(Q11-Q21*zi2)*zi; - zk=cexp(I*(z-POF))*csqrt(zi)*C3; - *h0=zk*(p0z+I*q0z); - *h0p=I*zk*(p1z+I*q1z); - - if(ib == 0) - return; - - zms=cos((sqrt(zms)-4.)*31.41592654); - *h0=.5*(y0*(1.+zms)+ *h0*(1.-zms)); - *h0p=.5*(y0p*(1.+zms)+ *h0p*(1.-zms)); - - return; + int k, ib; + static int m[101], init = FALSE; + static double a1[25], a2[25], a3[25], a4[25], psi, tst, zms; + complex double clogz, j0, j0p, p0z, p1z, q0z, q1z, y0=0, y0p=0, zi, zi2, zk; + + /* initialization of constants */ + if( ! init ) + { + int i; + psi=-GAMMA; + for( k = 1; k <= 25; k++ ) + { + i = k-1; + a1[i]=-.25/(k*k); + a2[i]=1.0/(k+1.0); + psi += 1.0/k; + a3[i]=psi+psi; + a4[i]=(psi+psi+1.0/(k+1.0))/(k+1.0); + } + + for( i = 1; i <= 101; i++ ) + { + tst=1.0; + for( k = 0; k < 24; k++ ) + { + init = k; + tst *= -i*a1[k]; + if(tst*a3[k] < 1.e-6) + break; + } + m[i-1]=init+1; + } + + init = TRUE; + + } /* if( ! init ) */ + + zms=cAbs2(z); + if(zms == 0.) + abort_on_error(-7); + + ib=0; + if(zms <= 16.81) + { + if(zms > 16.) + ib=1; + + /* series expansion */ + int iz=(int)zms; + int miz=m[iz]; + j0=1; + j0p=j0; + y0=0; + y0p=y0; + zk=j0; + zi=z*z; + + for( k = 0; k < miz; k++ ) + { + zk *= a1[k]*zi; + j0 += zk; + j0p += a2[k]*zk; + y0 += a3[k]*zk; + y0p += a4[k]*zk; + } + + j0p *= -.5*z; + clogz=clog(.5*z); + y0=(2.*j0*clogz-y0)/PI+C2; + y0p=(2./z+2.*j0p*clogz+.5*y0p*z)/PI+C1*z; + *h0=j0+I*y0; + *h0p=j0p+I*y0p; + + if(ib == 0) + return; + + y0=*h0; + y0p=*h0p; + + } /* if(zms <= 16.81) */ + + /* asymptotic expansion */ + zi=1./z; + zi2=zi*zi; + p0z=1.+(P20*zi2-P10)*zi2; + p1z=1.+(P11-P21*zi2)*zi2; + q0z=(Q20*zi2-Q10)*zi; + q1z=(Q11-Q21*zi2)*zi; + zk=cexp(I*(z-POF))*csqrt(zi)*C3; + *h0=zk*(p0z+I*q0z); + *h0p=I*zk*(p1z+I*q1z); + + if(ib == 0) + return; + + zms=cos((sqrt(zms)-4.)*31.41592654); + *h0=.5*(y0*(1.+zms)+ *h0*(1.-zms)); + *h0p=.5*(y0p*(1.+zms)+ *h0p*(1.-zms)); + + return; } /*-----------------------------------------------------------------------*/ @@ -669,9 +694,9 @@ static void hankel( complex double z, complex double *h0, complex double *h0p ) /* compute integration parameter xlam=lambda from parameter t. */ static void lambda( double t, complex double *xlam, complex double *dxlam ) { - *dxlam=b-a; - *xlam=a+*dxlam*t; - return; + *dxlam=b-a; + *xlam=a+*dxlam*t; + return; } /*-----------------------------------------------------------------------*/ @@ -680,162 +705,163 @@ static void lambda( double t, complex double *xlam, complex double *dxlam ) /* the method of variable interval width romberg integration is used. */ static void rom1( int n, complex double *sum, int nx ) { - int jump, lstep, nogo, i, ns, nt; - static double z, ze, s, ep, zend, dz=0., dzot=0., tr, ti; - static complex double t00, t11, t02; - static complex double g1[6], g2[6], g3[6], g4[6], g5[6], t01[6], t10[6], t20[6]; - - lstep=0; - z=0.; - ze=1.; - s=1.; - ep=s/(1.e4*NM); - zend=ze-ep; - for( i = 0; i < n; i++ ) - sum[i]=0; - ns=nx; - nt=0; - saoa(z,g1); - - jump = FALSE; - while( TRUE ) - { - if( ! jump ) - { - dz=s/ns; - if( (z+dz) > ze ) - { - dz=ze-z; - if( dz <= ep ) - return; - } - - dzot=dz*.5; - saoa(z+dzot,g3); - saoa(z+dz,g5); - - } /* if( ! jump ) */ - - nogo=FALSE; - for( i = 0; i < n; i++ ) - { - t00=(g1[i]+g5[i])*dzot; - t01[i]=(t00+dz*g3[i])*.5; - t10[i]=(4.*t01[i]-t00)/3.; - - /* test convergence of 3 point romberg result */ - test( creal(t01[i]), creal(t10[i]), &tr, cimag(t01[i]), cimag(t10[i]), &ti, 0. ); - if( (tr > CRIT) || (ti > CRIT) ) - nogo = TRUE; - } - - if( ! nogo ) - { - for( i = 0; i < n; i++ ) - sum[i] += t10[i]; - - nt += 2; - z += dz; - if(z > zend) - return; - - for( i = 0; i < n; i++ ) - g1[i]=g5[i]; - - if( (nt >= NTS) && (ns > nx) ) - { - ns=ns/2; - nt=1; - } - - jump = FALSE; - continue; - - } /* if( ! nogo ) */ - - saoa(z+dz*.25,g2); - saoa(z+dz*.75,g4); - nogo=FALSE; - for( i = 0; i < n; i++ ) - { - t02=(t01[i]+dzot*(g2[i]+g4[i]))*.5; - t11=(4.*t02-t01[i])/3.; - t20[i]=(16.*t11-t10[i])/15.; - - /* test convergence of 5 point romberg result */ - test( creal(t11), creal(t20[i]), &tr, cimag(t11), cimag(t20[i]), &ti, 0. ); - if( (tr > CRIT) || (ti > CRIT) ) - nogo = TRUE; - } - - if( ! nogo ) - { - for( i = 0; i < n; i++ ) - sum[i] += t20[i]; - - nt++; - z += dz; - if(z > zend) - return; - - for( i = 0; i < n; i++ ) - g1[i]=g5[i]; - - if( (nt >= NTS) && (ns > nx) ) - { - ns=ns/2; - nt=1; - } - - jump = FALSE; - continue; - - } /* if( ! nogo ) */ - - nt=0; - if(ns < NM) - { - ns *= 2; - dz=s/ns; - dzot=dz*.5; - - for( i = 0; i < n; i++ ) - { - g5[i]=g3[i]; - g3[i]=g2[i]; - } - - jump = TRUE; - continue; - - } /* if(ns < nm) */ - - if( ! lstep ) - { - lstep = TRUE; - lambda( z, &t00, &t11 ); - } - - for( i = 0; i < n; i++ ) - sum[i] += t20[i]; - - nt++; - z += dz; - if(z > zend) - return; - - for( i = 0; i < n; i++ ) - g1[i]=g5[i]; - - if( (nt >= NTS) && (ns > nx) ) - { - ns /= 2; - nt=1; - } - - jump = FALSE; - - } /* while( TRUE ) */ + int jump, lstep, nogo, i, ns, nt; + static double z, ze, s, ep, zend, dz=0., dzot=0., tr, ti; + static complex double t00, t11, t02; + static complex double g1[6], g2[6], g3[6], g4[6], g5[6], t01[6], t10[6], t20[6]; + + lstep=0; + z=0.; + ze=1.; + s=1.; + ep=s/(1.e4*NM); + zend=ze-ep; + for( i = 0; i < n; i++ ) + sum[i]=0; + ns=nx; + nt=0; + saoa(z,g1); + + jump = FALSE; + while( TRUE ) + { + if( ! jump ) + { + dz=s/ns; + if( (z+dz) > ze ) + { + dz=ze-z; + if( dz <= ep ) + return; + } + + dzot=dz*.5; + saoa(z+dzot,g3); + saoa(z+dz,g5); + + } /* if( ! jump ) */ + + nogo=FALSE; + for( i = 0; i < n; i++ ) + { + t00=(g1[i]+g5[i])*dzot; + t01[i]=(t00+dz*g3[i])*.5; + t10[i]=(4.*t01[i]-t00)/3.; + + /* test convergence of 3 point romberg result */ + test( creal(t01[i]), creal(t10[i]), &tr, cimag(t01[i]), cimag(t10[i]), &ti, 0. ); + if( (tr > CRIT) || (ti > CRIT) ) + nogo = TRUE; + } + + if( ! nogo ) + { + for( i = 0; i < n; i++ ) + sum[i] += t10[i]; + nt += 2; + + z += dz; + if(z > zend) + return; + + for( i = 0; i < n; i++ ) + g1[i]=g5[i]; + + if( (nt >= NTS) && (ns > nx) ) + { + ns=ns/2; + nt=1; + } + + jump = FALSE; + continue; + + } /* if( ! nogo ) */ + + saoa(z+dz*.25,g2); + saoa(z+dz*.75,g4); + nogo=FALSE; + for( i = 0; i < n; i++ ) + { + t02=(t01[i]+dzot*(g2[i]+g4[i]))*.5; + t11=(4.*t02-t01[i])/3.; + t20[i]=(16.*t11-t10[i])/15.; + + /* test convergence of 5 point romberg result */ + test( creal(t11), creal(t20[i]), &tr, cimag(t11), cimag(t20[i]), &ti, 0. ); + if( (tr > CRIT) || (ti > CRIT) ) + nogo = TRUE; + } + + if( ! nogo ) + { + for( i = 0; i < n; i++ ) + sum[i] += t20[i]; + + nt++; + z += dz; + if(z > zend) + return; + + for( i = 0; i < n; i++ ) + g1[i]=g5[i]; + + if( (nt >= NTS) && (ns > nx) ) + { + ns=ns/2; + nt=1; + } + + jump = FALSE; + continue; + + } /* if( ! nogo ) */ + + nt=0; + if(ns < NM) + { + ns *= 2; + dz=s/ns; + dzot=dz*.5; + + for( i = 0; i < n; i++ ) + { + g5[i]=g3[i]; + g3[i]=g2[i]; + } + + jump = TRUE; + continue; + + } /* if(ns < nm) */ + + if( ! lstep ) + { + lstep = TRUE; + lambda( z, &t00, &t11 ); + } + + for( i = 0; i < n; i++ ) + sum[i] += t20[i]; + + nt++; + z += dz; + if(z > zend) + return; + + for( i = 0; i < n; i++ ) + g1[i]=g5[i]; + +// ADDA: the following was present in NEC2C v.0.2 but disappeared in v.1.3. Not clear, what is the implication +// if( (nt >= NTS) && (ns > nx) ) +// { +// ns /= 2; +// nt=1; +// } + + jump = FALSE; + + } /* while( TRUE ) */ } @@ -845,126 +871,157 @@ static void rom1( int n, complex double *sum, int nx ) /* integrals for source and observer above ground */ static void saoa( double t, complex double *ans) { - double xlr, sign; - static complex double xl, dxl, cgam1, cgam2, b0, b0p, com, dgam, den1, den2; - - lambda(t, &xl, &dxl); - if( jh == 0 ) - { - /* bessel function form */ - bessel(xl*rho, &b0, &b0p); - b0 *=2.; - b0p *=2.; - cgam1=csqrt(xl*xl-ck1sq); - cgam2=csqrt(xl*xl-ck2sq); - if(creal(cgam1) == 0.) - cgam1=cmplx(0.,-fabs(cimag(cgam1))); - if(creal(cgam2) == 0.) - cgam2=cmplx(0.,-fabs(cimag(cgam2))); - } - else - { - /* hankel function form */ - hankel(xl*rho, &b0, &b0p); - com=xl-ck1; - cgam1=csqrt(xl+ck1)*csqrt(com); - if(creal(com) < 0. && cimag(com) >= 0.) - cgam1=-cgam1; - com=xl-ck2; - cgam2=csqrt(xl+ck2)*csqrt(com); - if(creal(com) < 0. && cimag(com) >= 0.) - cgam2=-cgam2; - } - - xlr=cAbs2(xl); - if(xlr >= tsmag) - { - if(cimag(xl) >= 0.) - { - xlr=creal(xl); - if(xlr >= ck2) - { - if(xlr <= ck1r) - dgam=cgam2-cgam1; + double xlr; + static complex double xl, dxl, cgam1, cgam2, b0, b0p, com, dgam, den1, den2; + + lambda(t, &xl, &dxl); + // evaluate gamma1 and gamma2 (cgam1, cgam2) + if( jh == 0 ) + { + /* bessel function form */ + /* Assuming k1,2 to belong to the first quadrant of the complex plane, the following expressions have branch + * cuts from k to i*inf in an arc that stays within the first quadrant with Im(z)>=Im(k) and symmetrically from + * -k to -i*inf. + * For Re(k)=0, this is exactly the vertical line from k to k+i*inf (as assumed in the theory). + * For Im(k)=0, it is a corner from k to 0 and further to i*inf. + * This is sufficient to avoid integration contour, used for the Bessel form (concentrated fully in the fourth + * quadrant). The exact zero falls on the branch cut for Im(k)=0, but the function value is zero in this point + * anyway. The last part of the code keeps the functions continuous at the boundary of the fourth quadrant, when + * Im(k)=0, although that does not seem strictly necessary for the employed contours. + */ + bessel(xl*rho, &b0, &b0p); + b0 *=2.; + b0p *=2.; + cgam1=csqrt(xl*xl-ck1sq); + cgam2=csqrt(xl*xl-ck2sq); + if(creal(cgam1) == 0.) + cgam1=cmplx(0.,-fabs(cimag(cgam1))); + if(creal(cgam2) == 0.) + cgam2=cmplx(0.,-fabs(cimag(cgam2))); + } + else + { + /* hankel function form */ + /* Assuming k1,2 to belong to the first quadrant of the complex plane, the following expressions have branch + * cuts from k vertically to k+i*inf and from -k horizontally to -k-inf. The first one is the same as in the + * manuals of NEC code (where the integration contours are discussed), while the second one is different. But + * it still stays in the third quadrant, which is sufficient for the used integration contours. + */ + hankel(xl*rho, &b0, &b0p); + com=xl-ck1; + cgam1=csqrt(xl+ck1)*csqrt(com); + if(creal(com) < 0. && cimag(com) >= 0.) + cgam1=-cgam1; + com=xl-ck2; + cgam2=csqrt(xl+ck2)*csqrt(com); + if(creal(com) < 0. && cimag(com) >= 0.) + cgam2=-cgam2; + } + + /* evaluate dgam=gamma2-gamma1 avoiding loss of precision for large |xl| + * It uses the expansion of gamma=sqrt(z^2-k^2) with error of O(k^8/z^7), which is valid (with sign=1) for |z|>>|k| + * and Re(z)>Re(k)sign(Im(z)), assuming k is in the first quadrant. + * If sign=-1, then the expansion works for |z|>>|k| and Re(z)= tsmag) + { + double sign; + if(cimag(xl) >= 0.) + { + xlr=creal(xl); + if(xlr >= ck2) + { + /* the original expression is fine between the two branch cuts, since then gamma1 is close to -gamma2 + * instead of gamma2 (for large |z|), i.e., they do not cancel each other + */ + if(xlr <= ck1r) + dgam=cgam2-cgam1; + else + { + sign=1.; + dgam=1./(xl*xl); + dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl; + } + } + else + /* this covers the whole second quadrant (which seems sufficient for all used contours), + * but the values for 0=Re(k2). Otherwise, + * additionally Re(xl)= ck2) */ + + } /* if(cimag(xl) >= 0.) */ + else + // this simple condition is sufficient for Im(xl)<0 if we assume that third quadrant is never reached + { + sign=1.; + dgam=1./(xl*xl); + dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl; + } + + } /* if(xlr < tsmag) */ else + dgam=cgam2-cgam1; + + // calculate integrand multipliers D1 and D2 (up to a factor of 2). D2 is rewritten to avoid loss of precision + den2=cksm*dgam/(cgam2*(ck1sq*cgam2+ck2sq*cgam1)); + den1=1./(cgam1+cgam2)-cksm/cgam2; + com=dxl*xl*cexp(-cgam2*zph); + ans[5]=com*b0*den1/ck1; + com *= den2; + + if(rho != 0.) { - sign=1.; - dgam=1./(xl*xl); - dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl; + b0p=b0p/rho; + ans[0]=-com*xl*(b0p+b0*xl); + ans[3]=com*xl*b0p; } - } - else - { - sign=-1.; - dgam=1./(xl*xl); - dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl; - } /* if(xlr >= ck2) */ - - } /* if(cimag(xl) >= 0.) */ - else - { - sign=1.; - dgam=1./(xl*xl); - dgam=sign*((ct3*dgam+ct2)*dgam+ct1)/xl; - } - - } /* if(xlr < tsmag) */ - else - dgam=cgam2-cgam1; - - den2=cksm*dgam/(cgam2*(ck1sq*cgam2+ck2sq*cgam1)); - den1=1./(cgam1+cgam2)-cksm/cgam2; - com=dxl*xl*cexp(-cgam2*zph); - ans[5]=com*b0*den1/ck1; - com *= den2; - - if(rho != 0.) - { - b0p=b0p/rho; - ans[0]=-com*xl*(b0p+b0*xl); - ans[3]=com*xl*b0p; - } - else - { - ans[0]=-com*xl*xl*.5; - ans[3]=ans[0]; - } - - ans[1]=com*cgam2*cgam2*b0; - ans[2]=-ans[3]*cgam2*rho; - ans[4]=com*b0; - - return; + else + { + ans[0]=-com*xl*xl*.5; + ans[3]=ans[0]; + } + + ans[1]=com*cgam2*cgam2*b0; + ans[2]=-ans[3]*cgam2*rho; + ans[4]=com*b0; + + return; } /*-----------------------------------------------------------------------*/ /* test for convergence in numerical integration */ -// ADDA: copied from nec2c.c +// ADDA: copied from calculation.c static void test( double f1r, double f2r, double *tr, - double f1i, double f2i, double *ti, double dmin ) + double f1i, double f2i, double *ti, double dmin ) { - double den; + double den; - den= fabs( f2r); - *tr= fabs( f2i); + den= fabs( f2r); + *tr= fabs( f2i); - if( den < *tr) - den= *tr; - if( den < dmin) - den= dmin; + if( den < *tr) + den= *tr; + if( den < dmin) + den= dmin; - if( den < 1.0e-37) - { - *tr=0.; - *ti=0.; - return; - } + if( den < 1.0e-37) + { + *tr=0.; + *ti=0.; + return; + } - *tr= fabs(( f1r- f2r)/ den); - *ti= fabs(( f1i- f2i)/ den); + *tr= fabs(( f1r- f2r)/ den); + *ti= fabs(( f1i- f2i)/ den); - return; + return; } /*------------------------------------------------------------------------*/ @@ -973,47 +1030,48 @@ static void test( double f1r, double f2r, double *tr, * * prints an error message and exits */ - +// ADDA: copied from misc.c, first 5 cases are redundant static void abort_on_error( int why ) { - switch( why ) - { - case -1 : /* abort if input file name too long */ - fprintf(stderr, "%s\n", - "nec2c: input file name too long - aborting"); - break; - - case -2 : /* abort if output file name too long */ - fprintf(stderr, "%s\n", - "nec2c: output file name too long - aborting"); - break; - - case -3 : /* abort on input file read error */ - fprintf(stderr, "%s\n", - "nec2c: error reading input file - aborting"); - break; - - case -4 : /* Abort on malloc failure */ - fprintf(stderr, "%s\n", - "nec2c: A memory allocation request has failed - aborting"); - break; - - case -5 : /* Abort if a GF card is read */ - fprintf(stderr, "%s\n", - "nec2c: NGF solution option not supported - aborting"); - break; - - case -6: /* No convergence in gshank() */ - fprintf(stderr,"No convergence in gshank() - aborting. Try to increase MAXH in somnec.c and recompile\n"); - break; - - case -7: /* Error in hankel() */ - fprintf(stderr, "%s\n", - "nec2c: hankel not valid for z=0. - aborting"); - break; - } /* switch( why ) */ - - /* clean up and quit */ - exit( why ); + switch( why ) + { + case -1 : /* abort if input file name too long */ + fprintf(stderr, "%s\n", + "nec2c: input file name too long - aborting"); + break; + + case -2 : /* abort if output file name too long */ + fprintf(stderr, "%s\n", + "nec2c: output file name too long - aborting"); + break; + + case -3 : /* abort on input file read error */ + fprintf(stderr, "%s\n", + "nec2c: error reading input file - aborting"); + break; + + case -4 : /* Abort on malloc failure */ + fprintf(stderr, "%s\n", + "nec2c: A memory allocation request has failed - aborting"); + break; + + case -5 : /* Abort if a GF card is read */ + fprintf(stderr, "%s\n", + "nec2c: NGF solution option not supported - aborting"); + break; + + // ADDA: updated comment to be specific to ADDA + case -6: /* No convergence in gshank() */ + fprintf(stderr,"No convergence in gshank() - aborting. Try to increase MAXH in somnec.c and recompile\n"); + break; + + case -7: /* Error in hankel() */ + fprintf(stderr, "%s\n", + "nec2c: hankel not valid for z=0. - aborting"); + break; + } /* switch( why ) */ + + /* clean up and quit */ + exit( why ); } /* end of abort_on_error() */ diff --git a/src/somnec.h b/src/somnec.h new file mode 100644 index 00000000..8232907f --- /dev/null +++ b/src/somnec.h @@ -0,0 +1,28 @@ +/* Header for calculation of Sommerfeld integrals + * + * Copyright (C) ADDA contributors + * This file is part of ADDA. + * + * ADDA is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. + * + * ADDA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License along with ADDA. If not, see + * . + */ +#ifndef __somnec_h +#define __somnec_h + +// system headers +/* ADDA uses doublecomplex macro instead of standard 'complex double', but we use the latter in somnec.c/h to + * keep them standalone + */ +#include + +void som_init(complex double epscf); +void evlua(double zphIn,double rhoIn,complex double *erv,complex double *ezv,complex double *erh,complex double *eph, + int mode); + +#endif // __somnec_h diff --git a/src/somnec_test.c b/src/somnec_test.c new file mode 100644 index 00000000..3dd0a742 --- /dev/null +++ b/src/somnec_test.c @@ -0,0 +1,77 @@ +/* This is a standalone module for testing the calculation of Sommerfeld integrals (in somnec.c/h). + * It is not yet included as a target in Makefile, so need to be compiled manually by a command like: + * gcc -O3 -Wall -o somnec_test somnec*.c + * + * Copyright (C) ADDA contributors + * This file is part of ADDA. + * + * ADDA is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. + * + * ADDA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License along with ADDA. If not, see + * . + */ +// project headers +#include "somnec.h" +// system headers +#include + +// useful macros for printing complex numbers +#define CFORM "%.10g%+.10gi" +#define REIM(a) creal(a),cimag(a) + +//====================================================================================================================== + +int main(int argc,char **argv) +{ + double mre,mim; + double rho,Z; + complex double eps,erv0,ezv0,erh0,eph0,erv1,ezv1,erh1,eph1,erv2,ezv2,erh2,eph2,erv3,ezv3,erh3,eph3; + // TODO: add meaningful error messages + if (argc != 5) return 1; + if (sscanf(argv[1],"%lf",&mre)!=1) return 2; + if (sscanf(argv[2],"%lf",&mim)!=1) return 2; + if (sscanf(argv[3],"%lf",&rho)!=1) return 2; + if (sscanf(argv[4],"%lf",&Z)!=1) return 2; + + eps=mre+I*mim; + eps=eps*eps; + som_init(eps); + + evlua(Z,rho,&erv0,&ezv0,&erh0,&eph0,0); + printf("Mode 0:\n"); + printf("Irv="CFORM"\n",REIM(erv0)); + printf("Izv="CFORM"\n",REIM(ezv0)); + printf("Irh="CFORM"\n",REIM(erh0)); + printf("Iph="CFORM"\n",REIM(eph0)); + + evlua(Z,rho,&erv1,&ezv1,&erh1,&eph1,1); + printf("\nMode 1:\n"); + printf("Irv="CFORM"\n",REIM(erv1)); + printf("Izv="CFORM"\n",REIM(ezv1)); + printf("Irh="CFORM"\n",REIM(erh1)); + printf("Iph="CFORM"\n",REIM(eph1)); + + evlua(Z,rho,&erv2,&ezv2,&erh2,&eph2,2); + printf("\nMode 2:\n"); + printf("Irv="CFORM"\n",REIM(erv2)); + printf("Izv="CFORM"\n",REIM(ezv2)); + printf("Irh="CFORM"\n",REIM(erh2)); + printf("Iph="CFORM"\n",REIM(eph2)); + + evlua(Z,rho,&erv3,&ezv3,&erh3,&eph3,3); + printf("\nMode 3:\n"); + printf("Irv="CFORM"\n",REIM(erv3)); + printf("Izv="CFORM"\n",REIM(ezv3)); + printf("Irh="CFORM"\n",REIM(erh3)); + printf("Iph="CFORM"\n",REIM(eph3)); + + printf("\nSum of absolute errors 0vs1: %g\n",cabs(erv0-erv1)+cabs(ezv0-ezv1)+cabs(erh0-erh1)+cabs(eph0-eph1)); + printf("\nSum of absolute errors 2vs1: %g\n",cabs(erv2-erv1)+cabs(ezv2-ezv1)+cabs(erh2-erh1)+cabs(eph2-eph1)); + printf("\nSum of absolute errors 3vs1: %g\n",cabs(erv3-erv1)+cabs(ezv3-ezv1)+cabs(erh3-erh1)+cabs(eph3-eph1)); + + return 0; +}