Optimize calculation of reflected Green's tensor

[GoogleCodeExporter]

Currently, somnec.c is used to compute a table of Sommerfeld integrals. There 
are two problems.

1) It is used as black box, without sufficient understanding of its accuracy. 
For some applications, the accuracy is probably not sufficient, for others - 
some of it can be sacrificed to get computational speed. This should be studied 
in more details (see comments in somnec.c). 

2) It takes considerable time, currently approximately equivalent to 200 
iterations. This can be prohibitive, e.g. for biological applications (large 
grid, but small number of iterations). So studying other alternatives 
(described in the beginning of issue 101) is still very relevant. However, the 
issue is not easy (preliminary work showed that it is not that easy to make a 
general "black box" routine with at least any automatic convergence control).

Alternative approach is to use 2D interpolating table to calculate Sommerfeld 
integrals - but that adds another free parameter, effecting the final accuracy.

Also, there is a connected (more specific) issue 175.

#101, #175

Original issue reported on code.google.com by yurkin on 25 Sep 2013 at 8:28

[GoogleCodeExporter]

We used the code from NEC2C from 
http://sharon.esrac.ele.tue.nl/users/pe1rxq/nec2c.rxq/nec2c.rxq-0.2.tar.gz. But 
later versions (1.3) are available at http://www.qsl.net/5b4az/ . It is not 
clear, however, if any changes to the Sommerfeld integrals part has been made.

Original comment by yurkin on 2 Feb 2014 at 9:33

[myurkin]

#175

[myurkin]

Concerning the optimization and accuracy, there are a lot of ideas in Schmehl's thesis (1994) - https://www.researchgate.net/publication/260426740_The_Coupled-Dipole_Method_for_Light_Scattering_from_Particles_on_Plane_Surfaces . It includes building a table and doing 2D interpolation on it.

The major question is what is the proper way to control accuracy. Constant relative error (as it seems to be now) may be not appropriate, since at large distances it is difficult to obtain and not really needed (since the values are then added to a much larger image-dipole term).

Probably, the best way would be to make an accuracy parameters (currently defines MAXH and CRIT in somnec.c) into function arguments to be adjusted at each call. Then a smaller CRIT can be used at larger distances.

This may help with problems as in #215

[myurkin]

Another relevant ideas for improvement (bug fixes) were proposed by Sheila Edalatpour in May 2014. I have not been able to look into it yet, so I am putting it here for completeness.

Currently, I am adding the surface interactions to the T-DDA, and I am using the Fortran version of the NEC for calculating the Sommerfeld integrals. I think you are using the C version of the same code in ADDA. I noticed that there is a bug in the code, which arises in some special cases. The bug causes the integration path to cross the branch cuts, and therefore the integration does not converge. I thought that you will be interested in knowing more about this mistake.

The problem is in evlua function in the somnec.c file. I attached two figures to the email. Figure 1 shows the general path of the integration for Hankel form of the Sommerfeld integrals and Fig. 2 shows the integration path when this mistake happens. If the condition ((ck2+1.)> creal(ck1)*1.01) in line 357 of the somnec.c is satisfied (which is the case when the real part of the refractive index of the surface is less than 1), the real part of the breakpoint (BK) is set to ck2+1.. The real part of the point CP3 is set to 1.02*ck2 in line 293. Therefore, if 0.02*ck2 is greater than 1 (which usually is), the real part of the breakpoint will be smaller than the real part of CP3 (please see Fig. 2). In this case the integration path goes toward the negative real axis (instead of going to + infinity) and it crosses the branch cuts. For fixing this problem, I simply changed the real part of the breakpoint to 1.03*ck2 (by modifying line 358 from rmis=ck2+1. to rmis=ck2*1.03) to make sure that the breakpoint is always on the right hand side of the CP3.

I also noticed that there is an error message in the Fortran version of the function rom1 (it alerts about the divergence of the integral), which does not exist in somnec.c. I think this message is very important. From this message I was able to find this error.

If you run the evlua function with the parameters listed below, you will be able to observe the bug.
Frequency = 9e13 (rad/sec)
Refractive index of the surface = 0.96+2.16*i
The position of the dipole i = (7e-8,2e-7,2.5e-7) (meter)
The position of the dipole j = (2e-7,6e-7,5e-7) (meter)

[myurkin]

Similar ideas for 2D interpolation are described in the following paper for another method. But is completely applicable to the DDA.
J. Waxenegger, A. Trügler, and U. Hohenester, “Plasmonics simulations with the MNPBEM toolbox: Consideration of substrates and layer structures,” Comput. Phys. Commun. 193, 138–150 (2015).

[myurkin]

/cc @anna-ae

[myurkin]

This will also allow using rectangular dipoles with unequal dX and dY in the presence of substrate. See #245

[anna-ae]

I added a preliminary implementation of the new method of calculating Sommerfeld integrals in my fork. The method includes applying a transformation to the integrals, calculating a 2D table and finding the values for all the necessary points by interpolation.

[myurkin]

A nice description of various asymptotics and singular points in the complex plane, relevant for evaluation of Sommerfeld integrals, is given in Section 4.4 of Osipov A.V. and Tretyakov S.A. Modern Electromagnetic Scattering Theory with Applications, John Wiley & Sons (2017).