Accurate characterization of other work in JOSS paper

This issue is part of my portion of the JOSS review, https://github.com/openjournals/joss-reviews/issues/5478

The current draft of the JOSS paper discusses POPPY and prysm (caveat: latter is my software)

Commercial software for GO and PO calculations and analysis, such as OpticStudio (Zemax) and GRASP (TICRA), has already been developed, but to our knowledge PyPO is the first free open-source Python package that simulates planar and quadric (off-axis) reflector geometries using both GO and PO. Moreover, PyPO does not employ approximations often employed by other software packages such as POPPy (Perrin & others, 2012) and Prysm (Dube, 2019). Rather, PyPO directly solves the radiation integral, allowing for propagation between multiple (off-axis) reflector surfaces

At minimum, this paragraph should describe what approximation POPPY and prysm make that this software does not make.

If the "approximation" is a wave optics approach, the history for your repo begins October 27, 2022. prysm's raytracer dates to December 2021 and can model all geometries supported by PyPO.

If the "approximation" is the solution to the Rayleigh-Sommerfeld integral in either the Fraunhofer or Fresnel regime, I think the above paragraph misrepresents things in a meaningful way. POPPY and prysm use Fourier transforms to 'directly solve' the Rayleigh-Sommerfeld integral as it is written down in e.g. Born & Wolf or Goodman. Those integrals are themselves approximations and assume that the cosine obliquity of the exit pupil is "small." When it is not-small, one can include the cosine obliquity, which renders the integral exact once more.

If you want to make a claim that your software makes fewer approximations, I think you should include a demonstration that it is more accurate, as the purpose of making less approximations is surely higher accuracy. The ray-trace result on tutorial 1, finding the focus, is significantly less accurate that one would expect for the focus of a parabola for an on-axis object at infinity. That calculation should be of order ~ floating point epsilon.

With approximation, we meant the treatment of electromagnetic fields as vectorial entities, of which the polarisation state can change due to propagation through systems containing curved surfaces.

I have studied the documentation of Prysm, and I think it is a very nice software. I am truly sorry if the paragraph in the paper misrepresented things in this way and I totally agree that this paragraph needs to be sharpened and revised for the JOSS paper.

From what I understood in the docs of Prysm, there are two fundamental differences (and one big overlap) in how Prysm and PyPO operate. In the following discussion, I will stick to using Prysm as the example because it is mentioned that Prysm is a superset of POPPY and PROPER, and I thus assume (and after having a look at the docs of both these packages) that Prysm, POPPY and PROPER work in a similar fashion, physics-wise.

Difference #1: polarisation It seems that Prysm does plane-to-plane propagation using the Rayleigh-Sommerfeld integral. This integral is the Green's function based solution to the scalar Helmholtz equation on a plane. The solution proper is a convolution product of the Green's function with the scalar electric field distribution on the source plane. PyPO, on the other hand, uses the equivalent surface current approach, which is a solution to the vectorial Helmholtz equation using induced surface currents on scattering objects. This method involves using the boundary conditions of Maxwells equations at a perfectly reflecting surface to calculate the electric and magnetic surface currents that are induced. The solution to this vectorial integral is also in the form of a Green's function (the overlap between Prysm and PyPO!) convolution, but the amplitudes and phases of the resulting EM-field components are affected by these surface currents because of the convolution of the Green's function with the induced surface currents.

It seems that in Prysm, polarisation is only taken into account when calculating transmission/reflection coefficients for thin films (I claim this by a search for "polarization" in the Prysm docs). Here, the polarisation is used in the context of the local plane of incidence, judging from the two options "p" and "s". This is the "approximation" we mentioned in the paragraph, the capability of PyPO to calculate (induced) polarisation.

Thus, PyPO is capable of propagating polarised EM fields (see the tutorials in the docs where we actually calculate this in the form of the cross-polar efficiency of an optical system as a demonstration of this claim). This is interesting for applications where (induced) polarisation is of concern, such as design of radio telescopes and satellite dishes. This induced polarisation is brought forth by the curvature of optical elements which raises the next difference.

Difference #2: propagation and design In Prysm, propagation is done it seems by defining an aperture distribution on a plane. Then, aberrations are added to the model by expanding them into (orthogonal) polynomials acting on the shape of the wavefront. Then, using the Rayleigh-Sommerfeld integral with Fourier approach (with or without cosine obliquity), the source distribution is propagated through a system consisting of a number of planes, which I presume terminates at a plane as well. PyPO, however, propagates field distributions to (off-axis) quadric or planar surfaces, calculates the induced surface currents at these surfaces using the boundary conditions of Maxwell's equations at perfectly reflecting interfaces, and uses these calculated currents to calculate the EM-field incident on the next reflector, which can be curved as well. The placement of these curved reflectors in 3D space and the ability to visualise these curved surfaces also seems to be a difference between PyPO and Prysm. Therefore, I think that another differnce between PyPO and Prysm is in how curved surfaces are handled and how this influences the design process.

In summary, I think Prysm and PyPO are two different beasts, designed for different purposes, with different design principles.

Would it be acceptable for you if we replace the quoted paragraph in the JOSS paper by the following:

Commercial software for GO and PO calculations and analysis, such as OpticStudio (Zemax) and GRASP (TICRA), has already been developed. Of these, GRASP also employs the method of equivalent surface currents. Open-source optical simulation software packages such as POPPy (Perrin & others, 2012) and Prysm (Dube, 2019) are also available. These open-source software packages use the Rayleigh-Sommerfeld integral equation, which is a solution to the scalar Helmholtz equation, and thus do not model polarisation of the propagated fields. PyPO aims to contribute to this open-source optical simulation ecosystem by offering a software package working on the principle of the equivalent surface currents, which allows for a fully vectorial propagation between multiple planar or (off-axis) quadric reflector surfaces. This is interesting in cases where polarisation is of concern, for example for the design of radio telescopes and satellite dishes, where receivers are highly sensitive to polarisation.

That point about the ray-trace focus finding being inaccurate was due to a default value set to 1e-2, a typo on our side. Now, in the current tutorial (not the GUI one I believe, this still needs to be fixed), the focus finding is more accurate (defaults to 1e-12).

It seems that Prysm does plane-to-plane propagation using the Rayleigh-Sommerfeld integral. This integral is the Green's function based solution to the scalar Helmholtz equation on a plane. The solution proper is a convolution product of the Green's function with the scalar electric field distribution on the source plane. PyPO, on the other hand, uses the equivalent surface current approach, which is a solution to the vectorial Helmholtz equation using induced surface currents on scattering objects. This method involves using the boundary conditions of Maxwells equations at a perfectly reflecting surface to calculate the electric and magnetic surface currents that are induced. The solution to this vectorial integral is also in the form of a Green's function (the overlap between Prysm and PyPO!) convolution, but the amplitudes and phases of the resulting EM-field components are affected by these surface currents because of the convolution of the Green's function with the induced surface currents.

I do not say this to split hairs, but the Rayleigh-Sommerfeld integral describes more than just plane-to-plane propagation, it also describes propagation between concentric spheres exactly. In the event that you are limited in a small regime on one of the spheres, it is also approximately the solution for sphere-to-plane or vice versa. This is the Fresnel approximation. In Hopkins' Image formation by a general optical system papers he covers how this is essential to being able to use the integral over the exit pupil exit sphere to represent a focusing propagation of arbitrary NA, provided the region of interest in the focal region is comparatively small, including arbitrary off-axisness and obliquity, provided the plane of interest is tangent to the chief ray.

prysm, POPPY, PROPER, ... do not deal with surfaces in the same manner as a raytracer; there is no concept of a generalized surface, just complex amplitude plates/sphereoids. In that sense, it is somewhat nonsensical to talk about paraboloids or hyperboloids, ... because that "primal" representation of a surface just is not how these programs work.

As a perhaps salient example, there is not really a way to actually represent a double gauss lens in any of the above programs, as concerns their physical optics capability. You would use a ray-tracer to trace rays to the exit pupil reference sphere for a given field point, and then you could propagate a field whose amplitude and phase were determined by the raytrace using wave optics. The non-wave-optics treatment of all of the lens elements means any wave behavior there is not in the simulation.

It seems that in Prysm, polarisation is only taken into account when calculating transmission/reflection coefficients for thin films (I claim this by a search for "polarization" in the Prysm docs). Here, the polarisation is used in the context of the local plane of incidence, judging from the two options "p" and "s". This is the "approximation" we mentioned in the paragraph, the capability of PyPO to calculate (induced) polarisation.

prysm contains an experimental polarization module, which includes an adapter that generalizes all of the propagation routines to operation on 2x2 Jones states. It is true that Ez is always assumed zero, the way it is written. This is not in the documentation, because it is not part of a release (and there has not been a new release in a while now).

That part of the tool is being used now by my intern to study polarization aberrations of large space telescopes and coronagraphs

In your edited paragraph, I would state (only) that the the other tools use Fourier methods to solve the Rayleigh-Sommerfeld integral, while your tool uses the surface currents approach. The version you have put above claims prysm does not do polarization, but the master branch right now contains a months-old polarization module. There are also polarization tools in various forks of POPPY, and in extensions of PROPER (such as the Roman-CGI models).

FYI, I see you have massively updated the docs (thank you!) -- but I went to look at the GUI tutorial to see the improved focus finding, and the images are all broken.

Ah, thank you for the detailed explanation. I will then update the paper with the edited paragraph and take your suggestion into account. I will let you know when its merged into main.

I will also have a look at the missing GUI images and fix that. The Jupyter notebooks do have the improved focus finder as well so there it can be checked.

And I am sorry about not mentioning Prysm's polarisation module. I could not find it in the docs so thought it was not there. But seeing now it is experimental, I fully understand that it is not included yet.

Could you let me know if you agree with the updated paragraph in the newly compiled JOSS paper as a more accurate state of the field? Then I will close this issue afterwards.

I believe the revised paragraph accurately captures the state of the field. As a grammar thing, at the moment it is now ... surface currents ... fourier methods ... surface currents. You may consider re-ordering the sentences to colocate all of the surface current related sentences.

PyPO-dev / PyPO

Accurate characterization of other work in JOSS paper #131