reply-to-cpc/main.tex

\documentclass[a4paper]{article}



%% Language and font encodings
\usepackage[english]{babel}
\usepackage[utf8x]{inputenc}
\usepackage[T1]{fontenc}

%% Sets page size and margins
\usepackage[a4paper,top=3cm,bottom=2cm,left=3cm,right=3cm,marginparwidth=1.75cm]{geometry}

%% Useful packages
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage[colorinlistoftodos]{todonotes}
\usepackage[colorlinks=true, allcolors=blue]{hyperref}
\begin{document}

\begin{center}
  Response Letter
\end{center}
Dear Editor,
\vspace{10pt}

We thank the referees for the detailed reports and constructive comments for our manuscript
entitled ``Mie calculation of electromagnetic near-field for a multilayered sphere'' 
by K. Ladutenko, U. Pal, A. Rivera and O. Pe\~na-Rodr\'iguez (CPC manuscript CPC-D-15-00354). In order to recognize their important contribution for improving the manuscript we have added to the ``Acknowledgments'' section the following text ``The authors thank the anonymous reviewers for their valuable comments.''. We have tried to address all the comments and the detailed response to all the comments as well as the summary of the changes made to the manuscript are included below. All the changes in the manuscript are marked in green.
\vspace{10pt}
\\
Sincerely Yours,\\
The authors\\

\vspace{10pt}

\newpage
\textbf{Reviewer \#1 comments}

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & There are several alternatives currently available, Mathematica GLMT Scripts, for example. It provides both near- and far-field results for multilayered sphere. Moreover, Mathematica can handle Bessel functions up to arbitrary precision without any difficulties, and the code is very transparent and compact. 
\end{tabular}
\vspace{0.5em}

We are aware of Mathematica GLMT Scripts and although it has interesting possibilities, we consider that both codes are not exactly comparable and, in any case, our code has some advantages:
\begin{itemize}
\item This script does not provide near-field evaluation inside the particle, at least none of the provided examples have it. Calculating the fields outside the particle is trivial and, indeed, it is possible with any Mie-type code because one only needs the scattering coefficients for this calculation.
\item Mathematica GLMT does not use Mathematica's ability to provide arbitrary precision for calculation of Mie coefficients for multilayer sphere. Actually GLMT Script itself references to our previous paper in CPC and uses the same evaluation of spherical functions via series with the same error due to $N_{\mathrm{stop}}$ selection. Equations needed to do the evaluation of the near-field distribution inside the multilayered spherical particle are provided in our present manuscript and, to the best of our knowledge, those optimized expressions are not currently available anywhere else.
\item You need to buy Mathematica license to use this script, which limits the availability of the code.
\item As it is published at \verb+http://photonicsdesign.jimdo.com/software/+ it does not provide any appropriate license conditions. Particularly, it is not clear if the code can be distributed,  modified (e.g. to put your own model parameters) or usable for academic or corporate research, etc. Our code Scattnlay is distributed under GPL (version 3) --- a well-known permissive license for open source software. 
\item Finally, it should be noted that the authors of Mathematica GLMT regularly use Scattnlay [M. Selmke, M. Braun, F. Cichos, arXiv:1105.3815v1 (2011); M. Selmke, M. Braun, F. Cichos, arXiv:1109.2772 (2011); M. Selmke, M. Braun, F. Cichos, ACS Nano 6 (2012) 2741–2749; M. Selmke, F. Cichos, J. Opt. Soc. Am. A 31 (2014) 2370–2384; M. Selmke, J. Quant. Spectrosc. Radiat. Transf. 162 (2015) 175–183] and their scripts (at least the initial versions) were heavily based on our code and, even now, most of their variables names are exactly the same that we use in Scattnlay. Hence, it is safe to say that Mathematica GLMT is a derivative work of our code, even if both projects have diverged over time.
\end{itemize}

The only code known to us that provides similar possibilities is MSTM by Mackowski and Mishchenko. However, it uses the T-matrix approach to do the evaluation, which is considerably less efficient for the case of concentric spheres and has no usage license defined. Actually, we  have verified our code against MSTM results too but choose not to include them in the manuscript due to abovementioned license restrictions of MSTM.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & The only advantage I see of the algorithm presented here is the speed, since it is implemented in the newest C++ language. But, to prove it benchmarking should be performed.
\end{tabular}
\vspace{0.5em}

We have dedicated a lot of effort in order to improve the speed and accuracy of our implementation. For instance, we have done an extensive optimization of the previous code that lead to more than a two-fold speed-up. However, this is by no means the main point of the manuscript. We thank the reviewer for the comment because it shows that our initial version of the manuscript does not stress the novelty enough. Some aspects of our work not available anywhere else are:
\begin{enumerate}
\item We have provided for the first time explicit expressions for Mie expansion coefficients inside the multilayer sphere.
\item We have proposed the calculation of the vector spherical harmonics in terms of the Ricatti-Bessel functions and their logarithmic derivatives and have checked the correctness of this approach.
\item We have verified the accuracy of the resulting implementation by comparing its results against those obtained with similar codes and full-wave 3D simulations.
\end{enumerate}
The manuscript has been modified in order to clarify the novelty.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & For currently presented test runs only core-shell structures are selected. The last example is, actually, from the same authors, which looks really confusing, since, I guess, they've  used exactly this code. So, they are test themselves with the same code - kind of weird. 
\end{tabular}
\vspace{0.5em}

Near-field calculations in our previous paper were obtained with full-wave commercial 3D electromagnetic simulation software CST MWS because this calculation was not available in Scattnlay at this time. Actually, this work inspired us to develop the near-field calculations in Scattnlay. However, the comment is very relevant because the lack of clarification in our manuscript can easily lead to misunderstanding. We have modified the manuscript after the corresponding reference ``reported recently [10]'' to include the following sentence ``, where the near-field distribution was calculated using full-wave commercial 3D electromagnetic simulation software CST MWS[reference to cst.com]''. Hopefully this will avoid confusions with this comparison.

On the other hand, we restricted ourselves mostly to compare against core-shell cases due to the lack of analytic software (including usage statements in the software license) to calculate the near-field of multilayered spheres but it should be noted that the results shown in Figure 4b correspond to a PEC sphere covered by 32 dielectric layers. There we verified the multilayered case against full-wave simulation done with finite-element method using CST MWS and finite-difference time-domain method using Lumerical FDTD. After doing this we are pretty sure that our code provides correct results (at least for the tested range of input parameters). As it was stated in ``Restrictions'' line of ``Program description'' section large complex part of layers refractive index leads to convergence problems.  Note that usage of full-wave 3D simulation required 4-5 orders of magnitude more computational resources (time) to get e.g. total scattering cross-section with a comparable accuracy. 

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & Also, in Fig. 2 they compare published and their results. But in all other figures they show only their results and refer to the literature for comparison. They should present all test and original results in the current manuscript.
\end{tabular}
\vspace{0.5em}

In Fig.2 we are not comparing with published results. Instead, we have downloaded the BHFIELD code developed by Suzuki and Lee and calculated the same example with our codes and with referenced software. To make it clearer we have modified the manuscript, replacing ``For this, we calculated electric and magnetic fields for the example provided with their computer code'' by ``For this, we used both codes to calculate electric and magnetic fields for the example provided with their computer code''. As for other figures, we only reference the relevant papers (see the corresponding references) because we do not have permission to reprint figures from those journals.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & On another note, since there are plenty of Mie based simulators which can be found at scattport.org I would encourage the authors to add GUI to their code to provide more flexibility and user-friendly environment. Otherwise, it's easier to code generalised Mie solution in Mathematica, which is perfectly suitable for it, and I don't see any advantage of the current program compare to others. Near field can now also be accessed in other programs I've mentioned and found at scattport.org
\end{tabular}
\vspace{0.5em}

We have checked each and every program listed at ``Mie type codes'' section (\href{http://www.scattport.org/index.php/light-scattering-software/mie-type-codes}{link}\footnote{http://www.scattport.org/index.php/light-scattering-software/mie-type-codes}) of Scattport before starting the development of our code (fall of 2014). While the software list is very large, most of the listed codes re-implement Mie solution as it was published in classical book of Bohren and Huffman ``Absorption and Scattering of Light by Small Particles'' (usually referenced as BHMIE) or re-implement MIEV0 code by Wiscombe. The few original approaches that we found were referenced in the manuscript as [11-18,22,23]. Note, that most of these implementations cannot evaluate field distribution inside the particle and those which can are restricted to the case of particles with one (bulk sphere) or two (core-shell) layers.

We totally agree with the reviewer in that it is a good idea to provide a GUI and such interface, named \href{http://www.scattport.org/index.php/programs-menu/mie-type-codes-menu/391-mielab}{MieLab}\footnote{http://www.scattport.org/index.php/programs-menu/mie-type-codes-menu/391-mielab}, exists for our previous code (restricted to far-field calculations). We have the intention of expanding MieLab to include near-field calculations but this will not be immediately ready because adding a GUI is very time consuming and we have no funding for this task.  Moreover, it should be noted that the paper describing our previous GUI-less code has more than 50 citations at the moment as it is listed in Scopus database. Hence, we believe that this code and the accompanying manuscript are valuable even without GUI. As an intermediate solution we provide bindings to Python language (and examples), which are the part of the released code. Python provides fairly good plotting possibilities (see figures 2,3, and 4 in the manuscript) via MatPlotLib module and its syntax is about the same complexity as used by Mathematica or Matlab. However, unlike those programs, Python can be used free of charge.

\newpage
\section{Reviewer \#2}

First of all we would like to thank this reviewer for his/her detailed review that has helped us to greatly improve the manuscript.


\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 
The paper does not make any claim to originality as far as the algorithm is concerned.
\end{tabular} 
\vspace{0.5em}

Indeed we believe that our algorithm for the near-field calculations is an improvement over the state of the art and that there are some original aspects to it. Next we copy the answer to a similar point raised by Reviewer 1:

We have dedicated a lot of effort in order to improve the speed and accuracy of our implementation. For instance, we have done an extensive optimization of the previous code that lead to more than a two-fold speed-up. However, this is by no means the main point of the manuscript. We thank the reviewer for the comment because it shows that our initial version of the manuscript does not stress the novelty enough. Some aspects of our work not available anywhere else are:
\begin{enumerate}
\item We have provided for the first time explicit expressions for Mie expansion coefficients inside the multilayer sphere.
\item We have proposed the calculation of the vector spherical harmonics in terms of the Ricatti-Bessel functions and their logarithmic derivatives and have checked the correctness of this approach.
\item We have verified the accuracy of the resulting implementation by comparing its results against those obtained with similar codes and full-wave 3D simulations.
\end{enumerate}
The manuscript has been modified in order to clarify the novelty.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & The authors’ citation list does acknowledge some of this earlier work.
The first recursive algorithm for the multi-layered sphere known to this
reviewer was published by Wait [2] in 1963, although this work does not
address the numerical issues tackled by the approach of Yang. The algorithm
is described in a 1970 textbook [3].
\end{tabular}
\vspace{0.5em}

We already provided Ref.~13 in the manuscript to Wait's paper. However, we found a  small typo in our manuscript that has been corrected. According to the journal`s web-site the  Wait paper was published in 1962 \href{http://link.springer.com/article/10.1007/BF02923455}{http://link.springer.com/article/10.1007/BF02923455}).

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 
The introductory section of the write-up includes an incidental citation of the
authors’ reference [18], which is a paper on the T-matrix method. This reviewer understands the T-matrix method as being applicable to a collection
of spheres side-by-side, not to a multi-layered sphere defined by concentric
spherical shells. If the authors agree, reference [18] should be omitted.
\end{tabular}
\vspace{0.5em}

The only restriction of T-matrix formulation as stated in Ref.~18 is that sphere surfaces should not intersect. This way, it is possible to define in MSTM several spheres embedded into each other to simulate a stratified sphere. Nevertheless, it should be noted that T-matrix is considerably less efficient than Mie theory for the case of concentric spheres. Moreover, MSTM does not provide any appropriate license condition, particularly it is not clear if it is valid or not to distribute, modify (e.g. to put your own model parameters) or use this code for academic or corporate research, etc. Scattnlay is distributed under GPL (version 3) --- a well-known permissive license for open source software. See the details at the MSTM web-page \href{http://www.eng.auburn.edu/~dmckwski/scatcodes/}{http://www.eng.auburn.edu/\char`\~dmckwski/scatcodes/} 

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 

This reviewer approves of the level of detail provided by the authors in
their derivation of the background theory and their algorithm. However, he
makes the following suggestion to improve the logical progression of the development and to avoid the near-repetition of the definition of the spherical
waves occurring in equations (2.1)–(2.4) and (9.1)–(9.4).
First, the reader could be given more guidance to the source of the
definitions by including the section in the citation of Bohren and Huffman
and the continue as follows: “Vector spherical harmonics can be calculated
using the expressions [12, Section 4.3], re-stated in terms of Ricatti-Bessel
functions for numerical reasons [17] as. . . ” followed by equations (9.1)–(9.4).

The definitions of $r_{n}$ , $\psi_{n}$ , $\zeta_{n}$ , $D^{(1)}_{n}$ and $D^{(3)}_{n}$ should be given at this early
point in the presentation so that they are ready to be used in equations
(5.1)–(5.4), (6.1)–(6.4) and (7.1)–(7.4), rather than being left until later.
If the authors choose to leave equations (2.1)–(2.4) as they occur in the
submitted version of the paper, they should correct the inconsistency in
notation: superscript (j) is used in (2.1) but (i) is used in (2.2)–(2.4).
\end{tabular}
\vspace{0.5em}

We have modified the paper according to the referee suggestion and further simplified eq. 6.1 - 7.4 (now eq. 7.1 - 8.5). We are very grateful for this suggestion since the explanation is now clearer and much more compact.


\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & This reviewer failed to reproduce the derivation of equations (6.1)–(6.4)
until he realised that there is a typographical error in equation (7.4), which
should read
\begin{equation*}
  T_4(m_{l+1}x_l) =  c^{(l+1)}_{n}  D^{(1)}_{n}(m_{l+1}x_l)  \psi_{n}(m_{l+1}x_l) 
          - b^{(l+1)}_{n}  D^{(3)}_{n}(m_{l+1}x_l)  \zeta_{n} (m_{l+1}x_l)\:,
\end{equation*}

The authors must correct this error and confirm that the correct T 4 has been
used in the code.
\end{tabular}
\vspace{0.5em}

Thank you for pointing this out, for sure this was a typo  that we have corrected in the updated manuscript. The code has the correct equation from \href{https://github.com/ovidiopr/scattnlay/blob/c72d4e94dee437cb6a4f170580dc13a4a4638f05/src/nmie-impl.hpp#L960}{the beginning}\footnote{https://github.com/ovidiopr/scattnlay/blob/c72d4e94dee437cb6a4f170580dc13a4a4638f05/src/nmie-impl.hpp\#L960 } (the link goes to the exact position in the code). Our derivation with the correct equation is also available \href{https://github.com/ovidiopr/scattnlay/blob/c72d4e94dee437cb6a4f170580dc13a4a4638f05/doc/EvalField.ipynb}{online}\footnote{https://github.com/ovidiopr/scattnlay/blob/c72d4e94dee437cb6a4f170580dc13a4a4638f05/doc/EvalField.ipynb}, see the output of the cell at the end of the file just after the word ``Finally''. The result was re-checked by substituting into initial equations and the error was corrected in the revised manuscript.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 
The authors have repeated the misunderstanding of Yang [1] regarding 
the justification of the boundary conditions at the centre of the sphere:
$a_n^{(1)} = 0$ and $b_n^{(1)} = 0$. The correct reason is given by Bohren and Huffman [4]
in their Sections 4.3 and 8.1. The point is that only the $j_n (kr)$ dependence
gives a finite, non-singular behaviour at the origin, and so can be the only
radial dependence occurring in the spherical wave expansion of the fields
in the core, region 1. A radial dependence of $j_n (kr)$ represents a standing
wave, that is, a superposition of both incoming and outgoing waves.
\end{tabular}
\vspace{0.5em}

We thank the reviewer for the provided explanation and we have included in the revised manuscript the discussion about the singularity of Bessel functions. However, in addition to this purely mathematical explanation we would like to have a more physical argument. It is clear for us now, that this part of the manuscript was badly written, a more correct way is  to say ```This can be understood in a more physical way by considering that there are no outward waves entering in the inner layer.''' 

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 
The authors mention that the code can accommodate a perfectly-conducting shell. Perfect conductor (PEC) is an idealisation that must be distinguished from real conductors. The only place a PEC layer should appear in
the model is at the core. It is futile to enclose other materials inside a PEC
shell since such a shell would exclude all fields from its interior.
There are several omissions from the reference list.
\end{tabular}
\vspace{0.5em}

Yes, we are aware of the fact that any layer below PEC will have no influence in the optical response. This has been discussed in our previous paper [J. Appl. Phys. 116 (2014) 184508] and for this reason we have not repeated the explanation here. However, from a computational point of view it might be interesting to be able to define any layer composed of PEC (for instance, to simplify de definition of some structures). In any case, the code handles this by simply ignoring any layer below the PEC. So, in this way the user has not to worry about this and the physical meaning is preserved.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 
No DOIs are given for [8], [9] and [22] and no ISBN is given for [12]. Several authors have second
initials which have been omitted: J.R. Wait for [13] and M.V. Bashevoy,
V.A. Fedotov and N.I. Zheludev for [20]. The final page number for [19]
should be given; this reviewer cannot ascertain whether it is 4735 or 4736
since he does not have ready access to the full paper. The full citation
for [21] is Nat. Commun. 5:3402(2014). [26] appears to be an unpublished
report, so readers should be told how to access it; a URL should be given.
\end{tabular}
\vspace{0.5em}

We have added DOIs for references [8] and [9] but reference [22] has no DOI. Moreover, all the other typos related to missing initials and incomplete citations have been fixed in the revised manuscript.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 
The Code Package
This reviewer has successfully compiled the standalone programs and run
the supplied tests. However, he has failed to compile the python and cython
extensions.
The authors should put much more effort into providing a greatly expanded README file, the first place any new user would look. Step-by-step
instructions should be provided in detail for the compilation of the various
programs in the package including any unusual dependencies on external libraries and commands. Step-by-step instructions for running the programs
should also be given, with an explanation of the command-line arguments
and expected outputs.
\end{tabular}
\vspace{0.5em}

We would be grateful if the reviewer can provide us more information about the error that he/she found when compiling the Python/Cython extensions. In addition to the README file we have provided a Makefile that should completely automate the process and gives a proper message on failure, including the name of the missing package (if any). Moreover, github provides several forums where we are always available to provide assistance with the compilation. We will be happy to improve the process when we have more information but this is more of a computational improvement rather than a problem with the manuscript.

\vspace{0.5em}
\begin{tabular}[!H]{l|p{0.9\textwidth}}
\quad & 
Recommendations
This paper and code should not be published in their present form.
They may be suitable for publication if the authors accept the suggestions and undertake the extensions and revisions described above.
The editors should refer any revisions to this reviewer to confirm that
the authors have acted on his suggestions and taken his recommendations
seriously.
The authors have chosen to use make their code publicly available on
the Web as well as submitting it to the CPC Library. The editors must
confirm that this action does not prevent the addition of the code to the
CPC Library.
\end{tabular}
\vspace{0.5em}

Of course this has to be checked by the editors but we would like to point out that the same was true for the first version of our code (available at this time in Sourceforge). We feel that this is a plus instead of a problem because the latest code is always available. For example, during last summer (after the initial submission to CPC) we have implemented multiple-precision in our code. It is a purely technical improvement that has no impact on the manuscript but it gives the possibility of simulating large lossy objects. This code is already available at master branch for usage. 


%\bibliographystyle{alpha}
%\bibliography{sample}

\end{document}