Is this package alive?

[KronosTheLate]

Is anyone maintining this package? I am concidering implementing functionality found here, with some extensions, and I was wondering if forking this repo, making changes, and registering the package was a good wa to go about it, so as to not recreate functionality already found here.

[roflmaostc]

Would be cool to bring this back again.

Can be super illustrative with some fancy Pluto notebooks too!

@ngedwin98 what do you think? To make collaboration and visibility greater, we could move it to @JuliaPhysics

[ngedwin98]

Apologies for the very delayed reply. While I am not actively maintaining this project at the moment, I would love to see it registered and be used as a foundation for something useful. I'm also happy to discuss what needs to be done to get it there.

[KronosTheLate]

One of my main problems was the need to extract the matrices and multiply them manually, still having to reverse the order. I was hoping/expecting to be able to define a system in terms of components such as FreeSpace, in the order that the light sees them (e.g. in a vector, such as sys = [FreeSpace(1), OtherComponentsHere]]), and just apply that system to a ray state (propegated_beam = sys * initial_beam_statevector). That is the convenience/sugar I feel like one has the opportunity to add in such a package.

[KronosTheLate]

Perhaps it would make sense to make another github user a contributor, allowing them to register and maintain the package? Clearly things are moving so slow here that there is a need for more hands on deck xD (With all respect, we are all busy and doing this voluntarily)

[roflmaostc]

@KronosTheLate We can create under this? I've got permissions for it. https://github.com/JuliaPhysics/

But I'm afraid that we would need to adapt the code quite a bit since abstractly typed structs are not recommended See https://github.com/ngedwin98/ABCDBeamTrace.jl/blob/f7b3ab35586bdb48d52e17cd9355376bacf63a3f/src/elements.jl#L2

[KronosTheLate]

That would make a lot of sense. I have a few other projects going on, so I will not have time to contribute to the transition. But I am all for it, as this is a useful package and it would be great to see it improved and maintained!

[roflmaostc]

All right, I might give it a shot and get inspired by the source code here. The mathematics and physics behind it should be fully covered by Wikipedia, I guess.

[KronosTheLate]

I would believe so. All I know is that you have the state of a beam represented by 2 "coordinates": The position above the optical axis (can be negative), and the angle out from the optical axis. To see the position and angle after some transformation (free space propegation, crossing a material interface etc) you multiply a 2x2 matrix from the left with the 2x1 vector that represents the state.

Arbitrarily many transformations can be represented by a single 2x2 matrix, as many matrices can be multiplied before any state is actually known. But this package seems to mainly be concerned with easy construction of individual 2x2 matrices.

Important detail to get right: The order of matrix application. If the light undergoes transformations represented by matrices A, then B, then C (e.g. free space propegation, interface change, glass propegation), then the final state s_f can be found in terms of the initial state s_i as s_f = C B A s_i (multiplication is implicit). This is clear when remembering that matrix multiplication is associative, meaning that we can draw parenthesis as we like, so that we are computing C (B (A s_i)). First the "operator" A is "applied", then B, then C. The complete "system matric" would then be given by C B A.

[roflmaostc]

Just to catch up, I started to port some parts: https://github.com/JuliaPhysics/ABCDMatrixOptics.jl

Gaussian is not ported yet, the plotting mechanism isn't in.

I got inspired by this PR to include the GeometricBeam

[ngedwin98]

Sorry for the lack of activity on this; indeed, I haven't worked much with free-space optical systems recently. For what it's worth, I'm definitely in support of forking this to JuliaPhysics for visibility and collaboration. I'll probably be more helpful in terms of suggesting new features than direct code contributions (aside from the original code commits here), but happy to be plugged into any discussions moving forward. If a more formal fork into JuliaPhysics would be useful, I can facilitate that as well.

One thing I wanted to note just glancing at the early port by @roflmaostc (particularly the implementation of the RTM function): In this package, I was using a convention in which the transfer matrices are independent of the refractive index, and the refractive index is instead a property of the state (i.e., the ray or the Gaussian beam). As I understand it, this can be done because the refractive index is just a rescaling of the length scale, which can be absorbed into the propagation constant (wavelength or k-vector). Furthermore, its algebraic structure is just multiplication by scalars, so that's also why the Interface elements are defined by the ratio of input and output refractive indices, rather than the refractive index itself (as suggested by the notations η rather than n, respectively). I'm not sure if we want to keep doing this moving forward, but it led to a fairly clean workflow for the use cases I had back then. Also would be interested to know if there might be some use cases where this rescaling can't be done.

Regarding @KronosTheLate's point on multiplying RTMs and elements, perhaps the following might work?

1. Introduce a new "generic" Element whose action is defined solely by its RTM, rather than having its RTM analytically inferred from a single parameter (like a lens would be).
2. Define a reduce method for Element sequences that just returns a single Element of the sort above, defined via an RTM equal to the matrix product. This would include any ordering reversal that is relevant.
3. To access the new RTM, just access it from the new Element.

In total, it would be something like RTM(reduce(*, vector_of_elements...))or even easier RTM(*(vector_of_elements...)). This seems better than the less transparent approach of defining RTM on Vector{<:Element} manually. We would basically enforce that the algebra of Elements composes left to right despite using matrices (which compose right to left) as the underlying mathematical backend.

Incidentally, as I write this, any thoughts on RTM vs (probably the more standard) rtm, or abcd, etc.?

[roflmaostc]

Ah now I understand your \eta convention, I blindly translated it and was confused by it.

Yes, I didn't like this manual fixing of the interface refractive index. Because the interface matrix somehow depends on what you multiply it with (the refractive index of the beam slightly before).

This issue is fixed by \eta? Never saw that before, but could work.

[KronosTheLate]

I have few inntuitions on naming and syntax, as I have little experience in this field. I can mainly say that I agree with everything written so far.

Just to spitball, I have in my head an idea along the lines of

elements = [FreeSpace(1), ThinLens(1)] # in order encountered by ray, i.e. along positive z-axis

initial_state = [initial_r, initial_theta]

final_state = apply(elements, initial_state)

I am not sure if that is the most elegant way to do it, but it aligns well with my mental model of what happens in what order between which objects.

[roflmaostc]

Ok regarding the convention: M = Lens FreeSpace vector

We can it realize either as: [Lens, FreeSpace] or [Free Space, Lens]

One would be more aligned with optical convention left -> right. The other one would be more consistent with the matrix description

What's your favorite?

[ngedwin98]

@KronosTheLate To get that functionality, we could first implement the product operation on Element as discussed above, then just use the transform function already defined, which takes an optical element and a state and produces the transformed state:

elements = [FreeSpace(1), ThinLens(1)] # In order encountered by ray, i.e. along positive z-axis
initial_state = Beam(λ, z0, n0, x0, k0, q0) # Convert r and theta into beam parameters
final_state = transform(prod(elements), initial_state)

I think we want to make a distinction between the optical elements and the transfer matrices. As I believe @KronosTheLate suggests, it's much more intuitive for us to write down an optical system by describing the optical elements in the order that the light beam encounters them, so left to right. Thus, the way the product operates on Element sequences should be left to right, as in the above example. But to implement that using matrices, we need to multiply the corresponding matrices right to left, which would be an implementation detail for the product operator on elements. I think it should be possible to develop a nice interface for working with optical elements without ever touching the matrices, so access to that lower layer via RTM can become strictly optional for those who want to do something with the matrices for whatever reason.

[ngedwin98]

@roflmaostc Regarding η, yes, the line in which it's "handled" is here -- https://github.com/ngedwin98/ABCDBeamTrace.jl/blob/f7b3ab35586bdb48d52e17cd9355376bacf63a3f/src/beamtrace.jl#L17 Assuming we clearly define η as the ratio of the new refractive index to the old across the interface, then n/η produces the new refractive index. If there is a chain of such ratios in an optical system, then simply multiplying the η factors suffice to capture them all, without needing immediate access to any "current" refractive index.

[roflmaostc]

But would that cover tee use case if someone does an optical system like

[Interface(1.3), FreeSpace(100), Lens(100), Interface(1.0), FreeSpace(100)]?

So a system where I also use the Interface to "turn on" a certain material. This would be a water dipping microscope for example.

In your approach, the user would need to manually take care of the previous refractive index, correct?

I think we could do that but since we anyway do not really work with raw matrices, I like the approach more where the ray carries the current refractive index and plugs the right refractive index into the next element, as in the case for the interface.

But it also breaks the fact that the matrices are commutative since in that case one would need to do it really from right to left to carry the right n

[ngedwin98]

Regarding that example, if you want to move from vacuum into a medium with index 1.3, propagate for 100 mm, pass through a 100-mm Lens, and then move back into vacuum to propagate for 100 mm, I believe that would correspond to the system

[Interface(1.3/1), FreeSpace(100), ThinLens(100), Interface(1/1.3), FreeSpace(100)]

To get the full simulation, we would initialize the state as, say, state = Beam(1e-3, 1, 1) (corresponding to a 1-mm radius waist right before the interface for a wavelength of 1 micron and an initial refractive index of vacuum).

As far as whether the user has to track the refractive index, it would depend on what they want to do. If the actual refractive index at the end of propagation (or anywhere intermediate) is desired, then yes, the user needs to update the refractive index according to n *= η. In fact, this is what is done in the transform function to describe what happens to a beam. However, if the user simply wants to have an Element or an RTM corresponding to the transfer function through the system, then they only need to know the refractive index at the end relative to that in the beginning. In this example, that relative factor is 1. I sort of think of these two pictures as "Schroedinger-like" vs. "Heisenberg-like" -- the former picture tracks the current state of the beam or ray as well as the current refractive index, while the latter only tracks changes in these quantities, and are in fact linear operators describing such changes. Hopefully we can have this package seamlessly support both pictures.

Does that help? I can also try cooking up another more complicated example to see if we can get the bottom of it.

[roflmaostc]

Yeah, what I mean is that the user effectively has to define the ratio themselves in your picture. The ration 1/1.3 results from the fact that the present media is 1.3 and you go back to n=1. But it's manual step.

In my formalism [Interface(1.3), FreeSpace(100), Lens(100), Interface(1.0), FreeSpace(100)]? the code could handle tracking which media it is in, and the user would only define different interface (as in a ray tracing tool).

Does that make sense?

[roflmaostc]

But maybe we should go the simple approach and define the matrices exactly as they are defined on Wikipedia.

So a parallel glass plate with diameter d would get defined by:

GlassPlate(n1, n2, d) = [Interface(n1=n1, n2=n1), FreeSpace(d), Interface(n1=n2, n2=n1)].

Slightly verbose, but easy to understand.

[ngedwin98]

Ahh, I see what you mean now! Indeed, any time you have a pair of interfaces, you ultimately need to specify the "inner" index twice. I suppose this is due to the "local" nature of each RTM, and "fusing" through the index changes requires global information. I certainly wouldn't complain if I got some syntactic sugar for this, but I suspect it would make the formalism a little more complicated than the usual version (as defined on Wikipedia).

Your suggested use of a two-argument Interface is brilliant. No more digging through the source code for how I defined the ratio. :)

[roflmaostc]

Yeah maybe we can come up with a macro which would do the manual setting of the refractive indices:

Something like: @trace_n [Interface(n2), FreeSpace(100), Interface(n1)] which would convert it to the above expression.

Great, I think we converged and I'll revert my code to the cleaner solution haha

[roflmaostc]

So implemented it.

Next big feature would be inclusion of GaussianBeam since I only cover GeometricBeam effectively. Then we could also think about including DynamicQuantities to make it more user friendly.

Not sure how much work it is to get it done with the current types.

Plotting does not work yet, never done that. Soo if you're motivated to do it :laughing:

[roflmaostc]

Plotting seems to work for GaussianBeams

[KronosTheLate]

Wow, loving the activity on here! The plot looks great <3

[roflmaostc]

Copy pasted it from @ngedwin98 :smiley:

Now trying to figure out, how the new package extensions scheme fits in there for the Plotting.

[roflmaostc]

Figure it out. There is now a weakdeps on Plots.

Whenever Plots is loaded, it automatically loads ABCDMatrixOpticsPlotExt which pulls in the Colors and Interpolations dependency. And it defines the plotting mechanism.

So loading time of the package should be fine. Plotting only works with GaussianBeam currently.

[roflmaostc]

I assembled also a small ToDo list.

So if you're keen @KronosTheLate you can help out! I might have some spare time to help out the next days. I could take care of docstrings, more testing, mirrors.

Nice (static) Pluto examples which can be accessed over the browser would be cool to have. Testing we can just check with codecov and add the missing lines.

Once we have the essentials, we can push it over the line and register it :smile:

[KronosTheLate]

I will be busy with other projects unfortunately. But thank you so much for your work!

[roflmaostc]

I think it's in a good shape to get registered. Feel free to play around with it, then I can trigger it next week https://github.com/JuliaPhysics/ABCDMatrixOptics.jl