HomerReid / scuff-em

A comprehensive and full-featured computational physics suite for boundary-element analysis of electromagnetic scattering, fluctuation-induced phenomena (Casimir forces and radiative heat transfer), nanophotonics, RF device engineering, electrostatics, and more. Includes a core library with C++ and python APIs as well as many command-line applications.
http://www.homerreid.com/scuff-em
GNU General Public License v2.0
125 stars 50 forks source link

Constant factor discrepancy in scuff-tmatrix? #141

Open odmiller opened 7 years ago

odmiller commented 7 years ago

Hi Homer,

I am using the delightful scuff-tmatrix tool. I am puzzled by what appears to be a constant multiplicative-factor discrepancy in the expected T-Matrix values and the known T-Matrix values for a sphere. For a sphere of radius r=1, constant permittivity epsilon=3, and a frequency k=6, I can use Equations 12(c) and 12(d) in your VSW technical memo to compute the non-zero diagonal elements of the T matrix. I find:

T_M = -0.96524023228608 + 0.183170757115865i T_E = -0.642139756924389 + 0.47937072240749i

for the \ell = 1 components. I have attached "scuff_tmatrix.m," which I used to compute these values. These values exactly agree with an independent Mie-theory T-Matrix code that I have. So I take these values to be correct. When I run

scuff-tmatrix --geometry sphere.scuffgeo --lMax 3 --Omega 6 --FileBase sphere_tmat

I get different values for the l=1 magnetic and electric T-Matrix elements:

T_M = 1.45 - 0.27i T_E = 0.94 - 0.73i

which are rough values that vary slightly for m=-1,0,1. I've investigated this discrepancy for a variety of sphere permittivities and frequencies, and I find that the values produced by scuff-em are consistently different from the analytical values by the factor

-k/4

where k is the frequency. Is there a reason for this discrepancy, and does it generalize to all structures? (Another possibility is that I've gotten the results of Eqs. (12c,12d) of your technical memo wrong, but the fact that they agree with an independent code suggests that they are correct.)

scuff_tmatrix_question.zip

[Note: I got the same results with a version of scuff-em from last October as well as the updated version as of today]

HomerReid commented 7 years ago

Thanks for pointing this out, and thanks especially for spending time and effort on detective work to diagnose the issue and compile a detailed bug report.

This is an old calculation that was never vetted particularly thoroughly. To refresh my memory, I just updated the spherical-wave memo with a new section (Section 10) outlining my strategy for computing T-matrix elements from surface currents. Equation (52) in the latest version of this memo corresponds to the calculation implemented in src/libscuff/GetSphericalMoments.cc, particularly equations 253-261. But it looks like the calculation as implemented in the code might be overcounting contributions---specifically, double-counting---compared to what is in the memo. So that would explain one factor of 2.

As for the spurious frequency factor, sure enough line 173 of src/applications/scuff-tmatrix/scuff-tmatrix.cc seems, for no apparent reason, to be multiplying the T-matrix expressions by a factor of -1.0Omega. (I see that this code comes with an archaic comment complaining that I don't understand the minus sign, and now I can conveniently augment that to say I don't understand the Omega factor either*, which feels at least like progress in the direction of consistency over the past 6 years.) At some point I might have gotten confused over the distinction between spherical multiple moments, which have units, and T-matrix elements, which are dimensionless. I might have inserted one layer too many of normalization factors.

At any rate, it certainly looks like there is a spurious factor of -1.0*Omega and possibly another factor of 2 here, so that gets us almost to the factor you identified. With a little more digging presumably the full issue can be resolved, at which point I'll update the code. In the meantime, I think the answers to your questions are: yes, there is a reason for this discrepancy (namely: sheer, unmitigated incompetence) and yes, it should generalize to all structures.

odmiller commented 7 years ago

Thanks for the quick response! I'll proceed cautiously but optimistically, and am happy to close the issue once it is fully resolved.

UmrathSt commented 6 years ago

I can also confirm this factor -omega/4 for a perfectly conducting sphere where for example the magnetic mie-coefficient b_l(kR) should read: b_l(kR) = - j(l+0.5, kR)/h1(l+0.5, kR) (Magnetic-Magnetic T-matrix entry for a PEC sphere), but comes out as: b_l(kR) * 4 / omega, with omega = k/c

HomerReid commented 6 years ago

In 4827a293510e55d1ad959ede5b0ed0eecddbd9ba I have I have overhauled the calculation of spherical multipole moments and revised the SCUFF-TMATRIX application code. I have also revamped the online documentation:

http://homerreid.github.io/scuff-em-documentation/applications/scuff-tmatrix/scuff-tmatrix/

I think all spurious factors have been eliminated and it is now doing the right thing, but it could always use another pair of eyes. Let me know there are any further issues, questions, comments, etc.

UmrathSt commented 6 years ago

Thank you, Homer. I can confirm that the TMatrix results for a PEC sphere now coincide with Mie-theory. Best regards

texnokrates commented 5 years ago

Hi, I think there is still something wrong. For a sphere with radius R > 1 I get T-matrix elements with absolute values greater than one, which is unphysical. The T-matrix values obtained by scuff-tmatrix seem to be R times higher than the values predicted by the Lorenz-Mie formulae. Here is an example with a sphere made of material with ε = 10 (background ε = 1), R = 1.5; dots are the scuff-tmatrix results, lines are the Lorenz-Mie predictions: e10sphere1 5 e_1 e10sphere1 5 m_1

scuff-tmatrix_inconsistency.tar.gz