Lays the foundation for further improvement (and corrections) in calc…

…ulation 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

@@ -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`

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 <float.h> // 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;