Astigmatic gausslets#54
Conversation
…nd detector integration - Added AstigmaticGaussianBeamlet with 9-ray structure (1 chief, 8 auxiliary). - Implemented synchronised tracing in System.jl to ensure all rays hit the same shapes. - Integrated AGB into the new Detector architecture (AstigmaticGaussianBeamletHit, field evaluation). - Added support for AGB in Beamsplitters (Thin, Plate, Cube) and PolarizationFilter. - Added Makie extension for AGB surface rendering. - Added comprehensive tests and modularized testing structure.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #54 +/- ##
==========================================
- Coverage 71.97% 71.63% -0.35%
==========================================
Files 57 61 +4
Lines 2901 3885 +984
==========================================
+ Hits 2088 2783 +695
- Misses 813 1102 +289 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: f154c5c9b8
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| # For auxiliary rays, we just need the new rays starting at npos | ||
| aux_ints = map(b -> begin | ||
| r = rays(b)[id] | ||
| BeamInteraction{T, Ray{T}}(nothing, Ray(npos, direction(r), wavelength(r))) |
There was a problem hiding this comment.
This rebuilds every auxiliary ray with the same origin npos from the chief ray, so after a polarization-filter hit all parabasal rays are forced to start at the chief intersection point instead of their own intersections. For astigmatic beamlets this collapses the transverse geometry and corrupts downstream waist/phase evolution whenever the auxiliary rays hit different points on the filter plane. Each auxiliary ray should be re-emitted from its own hit point (position(r) + length(r) * direction(r)).
Useful? React with 👍 / 👎.
| return PolarizedRay{F}( | ||
| Point3{F}(pos), | ||
| normalize(Point3{F}(dir)), | ||
| Point3{F}(dir), |
There was a problem hiding this comment.
The constructor now stores dir directly instead of normalizing it, which regresses prior behavior and makes PolarizedRay inconsistent with Ray. Many propagation/intersection calculations assume unit-length direction vectors; passing a non-unit dir now silently distorts geometric distances and optical-path-related phase calculations for polarized beams.
Useful? React with 👍 / 👎.
…e rotation correctly.
|
After some refactoring and a lot of headaches, the following script reproduces a Fraunhofer diffraction pattern like in this paper: # Reproduce fig. 12 of Harvey et. al. Optical Engineering 54(3), 035105 (March 2015)
using GLMakie, BeamletOptics
const BMO = BeamletOptics
const mm = 1e-3
const nm = 1e-9
## Figure 12: Fraunhofer diffraction of a square aperture
# D = 2 mm, λ = 632.8 nm, evaluated at 100 m
D = 2mm
λ = 632.8nm
L = 100.0 # meters
# N x N Gaussian beamlets
n_beams = 100
Δ = D / n_beams
w0s = Δ
xs = LinRange(-D / 2 + Δ / 2, D / 2 - Δ / 2, n_beams)
zs = LinRange(-D / 2 + Δ / 2, D / 2 - Δ / 2, n_beams)
beams = AstigmaticGaussianBeamlet[]
for x in xs
for z in zs
b = AstigmaticGaussianBeamlet([x, 0, z], [0, 1, 0], λ, w0s; E0 = [0, 0, 1.0])
push!(beams, b)
end
end
pd = Detector(3000mm)
translate_to3d!(pd, [0, L, 0])
system = System([pd])
# Trace all beams to the detector
solve_system!.(Ref(system), beams)
##
x, z, I = BMO.intensity(
pd, n = 300, x_min = -200mm, x_max = 200mm, z_min = -200mm, z_max = 200mm)
# Compare with analytical solution: I(x,z) ∝ sinc²(D*x/(λ*L)) * sinc²(D*z/(λ*L))
x_norm = x .* D ./ (λ * L)
z_norm = z .* D ./ (λ * L)
sinc_sq(u) = (u == 0) ? 1.0 : (sin(π * u) / (π * u))^2
analy_I = [sinc_sq(u_x) * sinc_sq(u_z) for u_x in x_norm, u_z in z_norm]
# Normalize numerical intensity to match analytical peak (1.0)
I_norm = I ./ maximum(I)
## 2D heatmap comparison
fig_2d = Figure(size = (800, 400))
ax_num = Axis(fig_2d[1, 1], title = "Numerical (AGB)", xlabel = "x (m)", ylabel = "z (m)")
heatmap!(ax_num, x, z, log10.(I_norm .+ 1e-15), colormap = :inferno, colorrange = (-4, 0))
ax_ana = Axis(fig_2d[1, 2], title = "Analytical", xlabel = "x (m)", ylabel = "z (m)")
heatmap!(ax_ana, x, z, log10.(analy_I .+ 1e-15), colormap = :inferno, colorrange = (-4, 0))
display(fig_2d)
##
# 1D comparison
# Take the slice at z = 0 (middle index)
mid_idx = argmin(abs.(z))
fig_1d = Figure()
ax_1d = Axis(fig_1d[1, 1], title = "Fraunhofer Diffraction Profile (z=0)",
xlabel = "x (m)", ylabel = "Normalized Irradiance", yscale = log10)
# The paper plots it on a log scale down to 1e-6
lines!(ax_1d, x, analy_I[:, mid_idx], label = "Analytical", color = :black, linewidth = 2)
lines!(ax_1d, x, I_norm[:, mid_idx], label = "Numerical (AGB)",
color = :red, linestyle = :dash, linewidth = 2)
ylims!(ax_1d, 1e-6, 1.5)
axislegend(ax_1d)
display(fig_1d)
|
|
Todos / Fixmes
|
…c field calculation
…/BeamletOptics.jl into AstigmaticGausslets
…/BeamletOptics.jl into AstigmaticGausslets
…' into AstigmaticGausslets
|
@codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 6e4db6eb41
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| p_proj = pt + dist * dir | ||
|
|
||
| v = p_proj - p0 | ||
| push!(pts_2D, Point2{T}(dot(v, ex), dot(v, ez))) |
There was a problem hiding this comment.
calc_local_pos for AstigmaticGaussianBeamletHit builds points from h1/h2, which are stored as Point3{Complex{T}}; this makes pt, p_proj, and the local dot results complex. Constructing Point2{T} from those values at this line will raise an InexactError whenever the imaginary part is nonzero (the normal astigmatic case), so spot_diagram/local-limit workflows can fail for astigmatic detector hits.
Useful? React with 👍 / 👎.
| # 3D waist point | ||
| pt = origin + waist_radius * (cos_t[k] * nd1 + sin_t[k] * nd2) | ||
| # Project onto detector plane along ray direction | ||
| dist = dot(p0 - pt, ey) / dot(dir, ey) |
There was a problem hiding this comment.
This projection now divides by dot(dir, ey) directly; for grazing/near-parallel incidence (dot(dir, ey) ≈ 0) it produces Inf/NaN projected points instead of a controlled failure path, which can corrupt spot bounds and field grids downstream. The previous implementation used line_plane_distance3d, which included a threshold check for this geometry edge case.
Useful? React with 👍 / 👎.
|
From my point of view this PR is good to go. The last remaining ToDo is:
|
TacHawkes
force-pushed
the
AstigmaticGausslets
branch
from
f8fdae6 to
0aea026
Compare
|
We might have to pre-render the double slit example. It seems to break doc CI. |
2 participants
PR Description - Astigmatic Gausslets (#54)
Summary
This PR introduces the
AstigmaticGaussianBeamlet(AGB) as a first-class beam type inBeamletOptics. AGBs extend the standard 3-rayGaussianBeamletto a 9-ray parabasal model, enabling analytical E-field computation for beams with elliptical cross-sections and independent phase curvature in two orthogonal transverse axes.This is required for physically correct simulation of:
Key Features & Improvements
Changes Since Initial Opening
z).