Working on loopless beamlet calculation in restructure branch - field isn't consistent

[Jashcraf]

Simulation parameters

• nrays = 24
• npix = 64
• wavelength = 1.65e-6

The field from the loopless calculation looks like a gaussian that drops to 0

Whereas the field from the current main branch (to be named legacy) Looks like what we expect

This issue is to track progress on fixing the loopless calculation

[Jashcraf]

I feel quite good about the code up to the construction of the orthogonal transformation matrices because the z-components that are supposed to be on the transversal plane are correctly zero-valued. It follows from this that the Delta calculation should be good, but maybe there's a rotation when the matrix multiplication happens that I didn't account for.

[Jashcraf]

Test the OPD

Restructure Branch

Main Branch

Apologies for the difference in ray distribution and colormap, but it appears we are loading the OPD correctly. Now I wonder if we can construct a similar plot of the Deltas?

EDIT: Turns out I didn't even implement fibonacci sampling, we should switch to cartesian on both for future tests.

[Jashcraf]

Testing Delta v.s. Detector Coordinate

Getting these to operate on the same beamlet might be a bit challenging, but if we stick to the same sample scheme I think it should work.

Restructure Branch

Main Branch

IT LOOKS LIKE THERE'S A MINUS SIGN. Turns out I flipped the LHS and RHS when defining the variables. If we correct that error it looks like this

So the Delta computation looks about right as well. In all cases all 4 of the beamlets look about the same, which is interesting. Though if they all come from a corner (which, they are in fig below) then it makes sense that the rotated detector would look about the same.

[Jashcraf]

Testing the ABCD Matrix computation

Restructure Branch

Main Branch

Now this is considerably different, so we have a place to start debugging the code.

[Jashcraf]

There was a slight discrepancy in the differentials because in the restructure branch I was using the normalized differentials and in the main I was using the un-normalized, which is correct. Change was made but it only had a small bearing on the final result

[Jashcraf]

Position and Directions of Base Ray on Transversal Plane

Restructure

keep in mind that there's more data here bc we aren't looping over nbeamlets

Main

The numbers aren't exactly the same but we see the same trend which (could) be attributed to rounding error at <1e-16

• x is consistent but slightly different
• y is the same between the two
• z looks different but contains a lot of the same values at numbers > ~1e-20. It might just be an ordering thing.

It might be worth investigating the actual rays, but first let's look at the same data for the differential rays.

[Jashcraf]

dPx - Looks good!

Restructure

Main

[Jashcraf]

dPy - looks good!

Restructure

Main

[Jashcraf]

dHx - W R O N G

already suspect of this because of the k vector

Restructure

Main

[Jashcraf]

dHy - W R O N G

It actually looks wrong by the same exact factor for both of them, what's going on?

Restructure

Main

[Jashcraf]

Let's start at the top (cue spiderverse) with the original ray data from the Zemax ray trace

Restructure

Main

So it actually looks like the positions and directions of the dL and dM rays are suspect. The positions are different by like a factor of 50?

[Jashcraf]

It appears that we have the same differential but a different initial rayset, what's causing this error?

[Jashcraf]

Update, wasn't converting the divergence angle into degrees. This was changed except for the differential that is ultimately passed to the abcd matrix computation. Now we are close, but not quite there.

Restructure

Main

What could it be? The ray transfer matrix calculation looks virtually identical now.

Restructure

Main

So everything before that must be basically fine. I think the remainder is to look for typos in the actual beam~

[Jashcraf]

Got part of the way by fixing the Qinv calculation, units were wrong so now that's fixed

[Jashcraf]

It looks like Qpinv is different between the two methods. Qinv isn't different so this is the problem line Qpinv = (C + D @ Qinv) @ np.linalg.inv(A + B @ Qinv)

All of these operations should broadcast because they are in the right order. But it appears that maybe I'm misunderstanding something.

[Jashcraf]

The broadcasted computations work exactly how I think they work. I guess the only thing remaining is to literally print the ray transfer matrices and see if those match up?

[Jashcraf]

It looks like there are very small differences at the 10^-10 level in the A ray transfer matrix, 10^-5 in the B ray transfer matrix, 10^-6 in the C ray transfer matrix (with a negative sign?).

[Jashcraf]

Okay here's something weird

Got a PSF that looks meaningful, but it was because of a typo in the amplitude Amplitude = 1/(np.sqrt(np.linalg.det(A + B @ Qpinv)))

In the main branch it's Qinv, which should be correct

So it's like this is cancelling with what the error is. For reference, when this is Qinv this is what we get

[Jashcraf]

For the time being I think I'm going to call this good. We have an accurate calculation using the matrix method of beamlet decomposition - so we can start profiling the code knowing that it produces physical results.

The "accelerated computing" paper can start going, commited this version of the code.

[Jashcraf]

Closing this for now because this doesn't seem to complicate the physics. Though it may be worth re-writing beamlets.py in the future. This is also a great example of a debugging journal.