[ANN] TransferMatrix.jl: a general 4 x 4 transfer matrix for optics of multilayer media

Package Announcements
, ,
1

I’m announcing my first Julia package, TransferMatrix.jl, a general 4 x 4 transfer matrix implementation in Julia.

The transfer matrix method propagates an electric field through a multi-layered medium. The results are the transmitted and reflected light for a given wavelength, and the electric field profile may also be calculated.

TransferMatrix.jl-quarter-wave
TransferMatrix.jl-quarter-wave800×600 49.3 KB

There are a lot of transfer matrix implementations on the web. It seems any graduate student in thin-films or optics writes a transfer matrix program. Why make another one? Many of these programs are one-offs with little documentation and are not very reusable. It is often difficult to get the program up and running or extend the code for a new research idea. Further, most use the 2 x 2 method and do not strive for full generality.

This implementation includes several key features:

Modular

The code is broken into small functions so that the user can easily implement their own solution at any step in the transfer matrix calculation by swapping out one or more functions.

General

This is a 4 x 4 transfer matrix based on the latest research trying to create a general implementation with broad applicability that is numerically stable and avoids singularities, a well-known problem for transfer matrix algorithms historically (a reference list is on the package website).

Documentation

Everything is fully documented and each piece of code that comes from a paper includes the reference and DOI at the function level. Furthermore, step-by-step tutorials are given to get the user running right away. All of the code on the tutorial website is runnable as-is.

User-friendly

It is very easy to build a multi-layered structure layer by layer. I have tried to make the boiler-plate code very low. This includes functionality to read refractive index files downloaded from refractiveindex.info. Complicated simulations are also easy to build and share with a YAML file. And together with Julia’s speed, it is fast and easy to create a transmittance / reflectance spectrum as a function of both wavelength and angle with respect to the incident light to produce 2D contour plots.

Reproducible

It is easy to make a complicated structure and share it with others with a YAML file. You can specify all of the simulation parameters (wavelength range, refractive indices, files, device angle range, etc.). It is simple to set up a periodic structure with an arbitrary number of periods in the file. This facilitates reproducibility and also makes it easy to try several variations and save them as separate config files.

Goals

My goal is for this transfer matrix program is to be readily reusable and extensible. While it may be instructive to create one’s own implementation, there is a lot of redundancy out there. I would hope that this package becomes at least a first stop for someone looking for a transfer matrix algorithm. If the Julia community thinks it is robust and useful, then I would like to see it become part of a set of synergistic science, physics, or optics packages in the ecosystem.

This is my first package written with a broad audience in mind. Everything else I’ve written has basically been just for me. I want as much feedback as I can get. In particular:

  1. I’m not certain that I’m taking full advantage of Julia’s multiple dispatch
  2. I’m not sure if the current way I’m using types is appropriate
  3. Can this code be made even more reusable by other researchers? What have I missed?

I would be very thankful to anyone generous enough to give their time to provide feedback. Check out the tutorials and documentation on the package website, built with Documenter.jl.

Similar packages

I could not find similar packages in Julia, but there are two that seem popular in Python:

12 Likes
2

You might consider computing scattering matrices (built-up layer-by-layer via Redheffer products) rather than transfer matrices. Unfortunately, getting the scattering matrix from the transfer matrix becomes exponentially ill-conditioned if there are evanescent waves, e.g. in a metal or in a photonic bandgap.

More generally, the 4x4 scattering matrix for multilayered media is a special case of the RCWA method (generalizing from translation-invariant to periodic layers), for which there are nowadays several implementations in C++ and Python, some of which are compatible with AD … would be great to see this in Julia.

8 Likes
3

Could you point me to a couple of popular scattering matrix packages? My main research is experimental ultrafast IR spectroscopy, so my knowledge of the difference between transfer/scattering matrices is limited. I’m aware of calculating scattering matrices from the transfer matrix elements, but I am not aware of the problems when doing this. I would want to take these into account if I write a scattering matrix program.

Is there a generally accepted method for computing the scattering matrix (more than just starting from the transfer matrix components)? Or would a sort of basic implementation of the scattering matrix still be useful? I would want to follow a standard procedure from a research paper or textbook. That said, I can imagine doing a WavePropagation.jl package and letting the user choose which computation to perform.

Related, is the scientific community moving away from transfer matrices and towards scattering matrices? Like transfer matrices are becoming irrelevant? Or is the situation more complicated than that?

References I know of:
I have the textbook Wave Propagation by Peter Marks and Costas Soukoulis, which has sections on scattering matrices, and I know of Raymond Rumpf’s lectures on scattering matrices (YouTube lectures here and I think he has lecture notes online somewhere).

1 Like
4

The Stanford S4 code is a good example of a free/open-source C++ RCWA scattering-matrix package. Its author wrote a short review of S-matrix calculations via Redheffer products, and one of my students put some more references in the Redheffer-product Wikipedia article — the technique is totally standard and accepted.

grcwa is a nice AD-compatible RCWA package in Python.

The advantages of Redheffer scattering-matrix methods over transfer matrices have been known for many decades now. S-matrix formulations have always been essential for RCWA methods, because as soon as you have finite-period unit cells you have lots of evanescent modes and transfer matrices quickly become a disaster. You can get away with transfer matrices for the simpler case of planar multilayer-film systems (or for more complicated structures if they are thin), and they are much simpler to explain and popular in textbooks, so you still see people using them — but they are terrible as soon as you start going to lots of layers, especially if you study photonic-bandgap systems. (I’ve seen people who didn’t know better resorting to arbitrary-precision transfer-matrix calculations in Mathematica to get around the ill-conditioning!) I think it’s pretty well known among computational-photonics people that S matrices fix these problems, but it hasn’t penetrated into the introductory pedagogy.

(IIRC, transfer matrices are okay as long as you don’t have any fields that decay by more than 1e-8 or so across your whole structure. After that point, the ill-conditioning starts to get so bad that you get noticably inaccurate results even in double precision, and eventually you get complete garbage. Unfortunately, for thick multilayer structures there are lots of situations in which you get enormous exponential attenuation.)

Physically, scattering matrices (which give you transmission and reflection coefficients), are usually what people want to compute at the end of the day. You can get them from the transfer matrices, but the problem is that if the transfer matrix is ill-conditioned then the resulting scattering matrix can be very inaccurate. It’s much better to build up the scattering matrix layer by layer using Redheffer products (although authors don’t always call it that — the same machinery seems to have been rediscovered independently multiple times since Redheffer in 1959).

10 Likes
5

In fact, I’m guilty of this myself on more than one occasion — I only learned about the Redheffer origins from Victor Liu of the Stanford S4 code. I once did a bunch of work with cylindrical multilayer films for Bragg fiber calculations (see also chapter 9 of our textbook), and initially I used 4x4 “Bessel” transfer matrices following the original 1978 paper by Yeh. But when we went to lots of layers, we eventually ran into numerical problems and I ended up re-inventing the S-matrix technique in our code (without knowing it).

A long time ago, I also implemented a cylindrical S-matrix scattering code in Julia for scattering off of 2d multilayer cylinders (1x1 matrices since it is 2d with one polarization, and there is only a single incident wave and a single reflected wave for a cylinder at a given angular order e^{im\phi}). The code is available as a Julia notebook, but beware that this is from 2014 so it won’t work in Julia 1.x without some updating.

5 Likes
6

As an experimental fluids guy who needs nothing more than 2x2 ray transfer matrices, I do wish someone would publish a package for that simple purpose. I’ve plugged the definitions in by hand each time I’ve needed them, but a package would be nice. I should probably just take it up myself.

Separately, note my interface package to refractiveindex.info: GitHub - stillyslalom/RefractiveIndex.jl: Interface to https://refractiveindex.info/. Please let me know if there are missing features or anything that’d make it more useful as a potential dependency instead of loading refractive index data ad hoc.

2 Likes
7

Rather than physics or optics, I came up through the microwave/antenna engineering community. In my discipline, the Redheffer star product has also been forgotten and rediscovered countless times. It has applications for microwave circuit interconnections and for stratified dielectric layers that occur in radomes, frequency selective surfaces, polarization selective surfaces, and so-called “meta surfaces”. I use this technique in my Julia package PSSFSS - analysis of polarization and frequency selective surfaces in Julia. The main outputs from PSSFSS are either the composite scattering matrix entries, or quantities directly derived from them. Although I do reference Redheffer, I also provide a derivation of the star product from first principles for the reader’s convenience in chapter 3 of the theory documentation.

The first reference to Redheffer’s formula that I can find in the antenna/microwave literature is in 1961 [1]. Years later, in 1986, the formula is re-derived without attribution by Chu and Ito [2] in a widely cited paper. The formula is derived yet again in 1987, this time for frequency selective surfaces, by Cwik and Mittra [3]. (IMO the notation used in [3] is poorly chosen and much more opaque than necessary). The book [4] makes explicit use of the Refheffer star notation in Section 6.3 and discusses how this formulation is superior to the technique of multiplying wave transmission matrices.

Considering its importance, it’s surprising that the Redheffer product isn’t more well known (at least in my area). I’m grateful for the Wikipedia article mentioned by @stevengj .

References:
[1] @ ARTICLE{stka:61,
author={Stock, D.J.R. and Kaplan, L.J.},
journal={IRE Transactions on Microwave Theory and Techniques},
title={A Comment on the Scattering Matrix of Cascaded 2n-Ports (Correspondence)},
year={1961},
volume={9},
number={5},
pages={454-454},
doi={10.1109/TMTT.1961.1125369}}

[2] @ ARTICLE{chit:86,
author={Tak Sum Chu and Itoh, T.},
journal={IEEE Transactions on Microwave Theory and Techniques},
title={Generalized Scattering Matrix Method for Analysis of Cascaded and Offset Microstrip Step Discontinuities},
year={1986},
volume={34},
number={2},
pages={280-284},
doi={10.1109/TMTT.1986.1133323}}

[3] @ article{cwmi:87,
author={Cwik, T. and Mittra, R.},
journal={IEEE Transactions on Antennas and Propagation},
title={The cascade connection of planar periodic surfaces and lossy dielectric layers to form an arbitrary periodic screen},
year={1987},
volume={35},
number={12},
pages={1397-1405},
doi={10.1109/TAP.1987.1144054}}

[4] @ book{vardaxoglou1997frequency,
title={Frequency selective surfaces: analysis and design},
author={Vardaxoglou, John C},
year={1997},
publisher={John Wiley and Sons}
}

4 Likes
8

Thank you for taking the time to explain this in such depth! You’ve convinced me that a scattering matrix implementation would be beneficial to have in Julia, rather than the transfer matrix. I’m surprised that the Redheffer star product has been rediscovered over and over (and in different sub-disciplines @PeterSimon!) and that nobody has gotten the message. It seems to me now that a robust and easy-to-use S-matrix program would be important not only for the community, but also to point to for pedagogical purposes.

I’m going to look over the resources you linked to (and others) and make a plan. I’m thinking I’ll work on it more publicly than I did for TransferMatrix.jl, since it seems from the comments here that there is interest in electrodynamics simulations and I would appreciate feedback along the way. I think I won’t have to modify the convenience features in TransferMatrix.jl – I will just need to write the S-matrix code.

Long-term, I’d like to see more physics packages in pure Julia instead of the split Python interface and C/C++ backends that is the current state of things. I think I can contribute from my corner of spectroscopy, so if a scattering matrix implementation would have broad appeal then it would be nice to have one as a part of the Julia physics ecosystem.

6 Likes
9

Your interface package does look interesting. I’ll take a look at it.

As for ray optics, does PhysicalOptics.jl have what you’re looking for?

1 Like
10

Do you think it would make sense to have a separate package just for the Redheffer star product? I don’t think it’s general enough to be part of the standard linear algebra package, but it seems like it’s used in more contexts than just scattering matrices.

11

I would say that it’s not worth a package, since Julia’s powerful and concise linear algebra syntax allows the star product to be implemented in just a few lines of code.

I took a quick look at your code and it appears that you are dealing with matrices of order 6 or less. In that case I suggest that you use StaticArrays instead of native Julia arrays. This will eliminate many allocations and result in much faster linear algebra operations.

3 Likes
12

I see, that is pretty simple. Thanks for linking to your code. Is your scattering matrix calculation specific to antenna design and electrical engineering? I ask because if you already have an efficient scattering matrix that can be used more broadly, I wouldn’t want to write redundant code. But it seems like your S-matrix calculation is pretty deeply embedded in the rest of the package and wouldn’t necessarily lend itself to stand-alone electrodynamics simulations?

Thank you for this suggestion. I didn’t know if it would be worth it to use StaticArrays and I was still learning so much about Julia while writing TransferMatrix.jl. I will make this a priority now, and will also use it in future similar projects.

2 Likes
13

No, but my situation is different than yours. As mentioned by @stevengl, my code deals with periodic unit cells, so the scattering matrix includes both propagating and evanescent modes. Therefore its order can be as large as 1000 or so; it needs to be realized as an ordinary Julia allocated array. On the other hand, for your TransferMatrix package at least, once you switch from transfer matrices to an S-matrix formulation, I believe your matrices S_{11}, S_{12}, S_{21}, and S_{22} will all be of order 4 (one TE and one TM mode at each side of a slab) which is well within the useful size of a StaticArray. So I think you want to write your own version that uses StaticArray matrices. But feel free to study my code if you wish and copy/use anything that looks useful. That’s why I gave it an MIT license–I hope it can be of use to others :smile:.

2 Likes
14

Other implementations of note, which haven’t been linked yet:

  1. Julia: GitHub - jonschlipf/RigorousCoupledWaveAnalysis.jl: Rigorous Coupled-Wave Analysis (RCWA) for nanophotonics simulations
  2. Julia: GitHub - lxvm/DeltaRCWA.jl: An nanophotonics solver for inverse design of metamaterials
    i. I believe the use-case is narrowly tailored here
  3. Autodiff: GitHub - Columbia-Computational-X-Lab/DiffSMat: Differentiable scattering matrix computation for designing photonic devices
    i. In contrast to grcwa above, this approach derives the forward differential of the exponential by hand, which the authors claim is more numerically robust than a pure autodiff approach because it stably resolves the issues around differentiating the eigendecomposition

Overall, RCWA (also known as the Fourier Modal Method) can be extended quite a bit from what you find off-the-shelf, including to perfectly matched layers, non-Fourier basis functions, etc.

5 Likes
15

That makes sense.

Thank you! I will absolutely be using it as a reference.

16

Since you mentioned it, PhysicalOptics is not maintained and not really useful currently.

I somehow struggled to get everything working with CUDA+Zygote in certain cases. Also I’m not quite sure which design (pure arrays, or own structs) is the best.

But what I need for my work is something like wave propagation and simulation of diffraction. And everything differentiable.

1 Like
17

Thanks for letting me know. I think it’s useful to know what physics packages are wanting in Julia. I don’t have a good sense of what to work on outside of my own research, but it would be nice to have a list and that way community efforts might be more focused.

18

When I created my package, I thought about parameterizing the element types of the scattering matrices (rather than hard-wiring them to ComplexF64) with a view towards making them amenable to auto differentiation (AD). But I decided not to for two reasons:

  1. As I understand it, AD doesn’t work for complex numbers.
  2. The matrices are fairly large (orders up to, say, 1000) and dense. In order to get maximal efficiency for the linear algebra operations I need, I want to use lapack/BLAS routines that can only work with native floating point types like ComplexF64. Even if AD supported complex types, I expect that with fancy dual numbers the linear algebra calculations would have to fall back to generic methods and be unbearably slow.

If my understanding is incorrect I’d like to hear about it and would be willing to modify my approach.

19

Yes those concerns were also regarding PhysicalOptics.jl :smiley:

  1. In principle e.g. Zygote supports complex numbers. There are only some issues with CUDA sometimes.

  2. BLAS routines also work with ComplexF32, don’t they? That should yield a speed-up in many cases where accuracy might not be an issue.

1 Like
20

Regarding point 1, thanks, that’s good to know. And wrt point 2, I hadn’t considered 32-bit. I’ve been using nothing but “double precision” since my Fortran days in the 80’s! I’ll have to give that some consideration.

TBH I have no experience coding or even using AD. When optimizing designs using my existing code, I use black box optimizers such as CMAEvolutionStrategy.jl. The nice thing about those optimizers is that they don’t require a smooth objective function. So I can add, say, discontinuous constraints without any problems.

1 Like