Astigmatic gausslets

.github/workflows/CI.yml

@@ -23,8 +23,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
+          - '1.12'
           - '1.10'
-          - '1.9'
           - 'pre'
         os:
           - ubuntu-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "BeamletOptics"
 uuid = "c387492b-ffec-4e51-9a46-bd230226031c"
-authors = ["Hugo Uittenbosch <hugo.uittenbosch@dlr.de>, Oliver Kliebisch <oliver.kliebisch@dlr.de> and Thomas Dekorsy <thomas.dekorsy@dlr.de>"]
 version = "0.12.1"
+authors = ["Hugo Uittenbosch <hugo.uittenbosch@dlr.de>, Oliver Kliebisch <oliver.kliebisch@dlr.de> and Thomas Dekorsy <thomas.dekorsy@dlr.de>"]
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -13,6 +13,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MarchingCubes = "299715c1-40a9-479a-aaf9-4a633d36f717"
 MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1"
 
@@ -30,18 +32,22 @@ ForwardDiff = "0.10, 1.0.1"
 GeometryBasics = "0.5"
 InteractiveUtils = "1"
 LinearAlgebra = "1"
+Logging = "1"
 Makie = "0.22, 0.23, 0.24"
 MarchingCubes = "0.1"
 MeshIO = "0.5"
 PrecompileTools = "1"
+Preferences = "1.5.2"
+Random = "1"
 StaticArrays = "1"
 Test = "1"
 Trapz = "2"
-julia = "1.9.4"
+julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test"]
+test = ["Aqua", "Test", "Logging"]

docs/make.jl

@@ -34,12 +34,18 @@ makedocs(;
                 "Aspherical lenses"         => joinpath("examples", "aspherical_lenses.md"),
                 "Double Gauss lens"         => joinpath("examples", "double_gauss.md"),
                 "Lens groups"               => joinpath("examples", "lens_groups.md"),
+                "Double slit"               => joinpath("examples", "double_slit.md"),
             ],
         ],
         "Basics" => Any[
             "Introduction"                  => joinpath("basics", "intro.md"),
             "Rays"                          => joinpath("basics", "rays.md"),
-            "Beams"                         => joinpath("basics", "beams.md"),
+            "Beams" => Any[
+                "Basic beam"                => joinpath("basics", "beams", "beams.md"),
+                "Stigmatic Gaussian"        => joinpath("basics", "beams", "stigmatic_beam.md"),
+                "Astigmatic Gaussian"       => joinpath("basics", "beams", "astigmatic_beam.md"),
+                "Beam groups"               => joinpath("basics", "beams", "beam_groups.md"),
+            ],
             "Optical components" => Any[
                 "Overview"                  => joinpath("basics", "components", "components.md"),
                 "Mirrors"                   => joinpath("basics", "components", "mirrors.md"),

docs/src/assets/beam_renders/agb_showcase.jl

@@ -0,0 +1,98 @@
+using GLMakie, BeamletOptics
+
+GLMakie.activate!(; ssao=true)
+
+const BMO = BeamletOptics
+
+const mm = 1e-3
+
+include(joinpath(@__DIR__, "..", "render_utils.jl"))
+
+##
+c_view = [
+  0.605459    0.783477  -0.139944  -0.0298959
+ -0.0484677   0.211806   0.976109  -0.0106539
+  0.7944     -0.584211   0.166213  -1.69474
+  0.0         0.0        0.0        1.0
+]
+
+agb = AstigmaticGaussianBeamlet([0,0,0], [0,1,0], 9e-4, 2e-3; support=[0,0,1])
+
+l = ThinLens(20e-3, 20e-3, BMO.inch, 1.5)
+s = System([l])
+
+translate3d!(l, [0, 30e-3, 0])
+
+fig1 = Figure(size=(600,300))
+display(fig1)
+ax = LScene(fig1[1,1])
+hide_axis(ax)
+
+render!(ax, agb; flen=0.1-0.03, show_beams=true, show_pos=true, color=RGBAf(1,0,0,0.1))
+
+set_orthographic(ax)
+set_view(ax, c_view)
+
+save("agbtest1.png", fig1; px_per_unit=8, update = false)
+
+##
+solve_system!(s, agb; check_invariant=false)
+
+fig = Figure(size=(600,300))
+display(fig)
+ax = LScene(fig[1,1])
+hide_axis(ax)
+render!(ax, agb; flen=0.1-0.03, show_beams=true, show_pos=true, color=RGBAf(1,0,0,0.1))
+render!(ax, l; color=RGBAf(1, 1, 1, .1))
+
+c_view = [
+  0.653709   0.756746  -5.27356e-16  -0.0173098
+ -0.352701   0.304679   0.884744     -0.00702154
+  0.669526  -0.578366   0.466077     -0.0223342
+  0.0        0.0        0.0           1.0
+]
+
+set_view(ax, c_view)
+save("agbtest2.png", fig; px_per_unit=8, update = false)
+
+##
+s1 = BMO.CylinderSDF(BMO.inch/2, BMO.inch)
+s2 = BMO.CylinderSDF(BMO.inch/2, BMO.inch)
+l1 = BMO.Lens(s1, n->1.5)
+l2 = BMO.Lens(s2, n->1.5)
+translate3d!(l1, [0, 50mm,  0])
+translate3d!(l2, [0, 100mm, 0])
+xrotate3d!(l1, deg2rad(90))
+xrotate3d!(l2, deg2rad(90))
+yrotate3d!(l1, deg2rad(90))
+translate3d!(l2, [5mm, 0, 0])
+
+system = System([l1, l2])
+
+z_os = 0e-3
+w0 = 2mm
+λ0 = 1000e-9
+dir = [0,1,0]
+E0 = [0,0,1]
+beam = AstigmaticGaussianBeamlet([0., 0, z_os], dir, λ0, w0; E0, support=[0,0,1])
+
+solve_system!(system, beam; check_invariant=false)
+
+##
+fig = Figure(size=(600,400))
+display(fig)
+ax = LScene(fig[1,1])
+hide_axis(ax)
+render!(ax, l1; color=RGBAf(1, 1, 1, .2))
+render!(ax, l2; color=RGBAf(1, 1, 1, .2))
+render!(ax, beam; color=RGBAf(1,0,0,0.2), z_res=4, r_res=20, show_waist=true, markersize=4, flen=0.1)
+
+view = [
+  0.763832   0.645415  3.88578e-16  -0.0541476
+ -0.444252   0.525762  0.725407     -0.0339012
+  0.468188  -0.554089  0.68832      -0.0967937
+  0.0        0.0       0.0           1.0
+]
+
+set_view(ax, view)
+save("agbtest3.png", fig; px_per_unit=8, update = false)

docs/src/assets/cond_save.jl

@@ -21,14 +21,37 @@ function conditional_include(fname::String; use_placeholder::Bool=!haskey(ENV, "
     # Check if to run script
     if use_placeholder        
         content = read(fname, String)
-        # match save("filename")
-        pattern = r"""save\(\s*"([^"]+)"""
-        filenames = [m.captures[1] for m in eachmatch(pattern, content)]
-        for f in filenames
-            fig = Figure(size=(600,300))
-            Box(fig[1, 1], color = :gray)
-            save(f, fig)
-            @info "Saving placeholder for $f"
+        # match save("filename") but ignore comments, i.e. # save(...)
+        pattern = Regex("^[ \\t]*save\\(\\s*\"([^\"]+)\"\\s*,\\s*([a-zA-Z_][a-zA-Z0-9_]*)", "m")
+        for match in eachmatch(pattern, content)
+            save_name = match.captures[1]
+            fig_name = match.captures[2]
+            offset = match.offset
+            prior_content = content[1:offset]
+            # find last missing fig_name variable match
+            npattern = Regex("\\Q$fig_name\\E\\s*=\\s*Figure\\(.*?size\\s*=\\s*\\(([\\d]+)\\s*,\\s*([\\d]+)\\)")
+            matches = collect(eachmatch(npattern, prior_content))
+            if !isempty(matches)
+                nmatch = last(matches)
+                _width =  parse(Int, nmatch.captures[1])
+                _height = parse(Int, nmatch.captures[2])
+            else
+                @info "Using default height for $fig_name"
+                _width =  600
+                _height = 300
+            end
+            # save placeholder
+            fig = Figure(size=(_width, _height))
+            ax = Axis(fig[1,1], backgroundcolor=:gray)
+            hidedecorations!(ax)
+            text!(ax, 0, 0;
+                text="$save_name\n($_width x $_height)",
+                align=(:center, :center),
+                color=:white,
+                fontsize=20
+            )
+            save(save_name, fig)
+            @info "Saving placeholder for $save_name (Size: $_width x $_height, fig_var_name=$fig_name)"
         end
     else
         @info "Running script for $fname"

docs/src/assets/examples/double_slit.jl

@@ -0,0 +1,131 @@
+using BeamletOptics, GLMakie
+
+function get_view(ls::LScene)
+    cam = ls.scene.camera_controls
+    eye = cam.eyeposition[]
+    lookat = cam.lookat[]
+    up = cam.upvector[]
+    return eye, lookat, up
+end
+
+function set_view(ls::LScene, eye, lookat, up)
+    cam = ls.scene.camera_controls
+    cam.eyeposition[] = Vec3f(eye...)
+    cam.lookat[] = Vec3f(lookat...)
+    cam.upvector[] = Vec3f(up...)
+    update_cam!(ls.scene, cam)
+end
+
+# Parameters
+λ = 633e-9          # 633 nm (HeNe Laser)
+a = 10e-6           # Slit width (10 µm)
+d = 100e-6          # Slit separation (100 µm)
+L = 0.2             # Propagation distance (0.2 m)
+W = 0.5e-3          # calculation window (0.5 mm)
+nx = 251            # resolution
+
+# Define the Aperture Mask (Double Slit)
+x = LinRange(-W / 2, W / 2, nx)
+z = LinRange(-W / 2, W / 2, nx)
+amplitude = zeros(nx, nx)
+
+for i in 1:nx, j in 1:nx
+    # Slits are long along z, narrow along x
+    in_slit1 = abs(x[i] - d / 2) < a / 2
+    in_slit2 = abs(x[i] + d / 2) < a / 2
+
+    if (in_slit1 || in_slit2)
+        amplitude[i, j] = 1.0
+    end
+end
+
+phase = zeros(nx, nx)
+dir = [0.0, 1.0, 0.0] # Propagation along Y
+##
+# Decompose into Astigmatic Beamlets
+@info "Decomposing Double Slit into AGBs..."
+# We explicitly define the basis to ensure Grid X = World X
+e1 = [1.0, 0.0, 0.0]
+e2 = [0.0, 0.0, 1.0]
+beams = WavefrontBeamletDecomposition(x, z, amplitude, phase, dir, λ;
+    threshold = 1e-3, overlap = 3.0, basis = (e1, e2))
+
+@info "Number of beamlets generated: $(length(beams))"
+
+# Define the System and Propagate
+target_pd = Detector(30e-3)
+translate_to3d!(target_pd, [0.0, L, 0.0])
+system = System([target_pd])
+
+@info "Propagating to screen at L = $L m..."
+solve_system!(system, beams)
+
+# Extract Results
+@info "Calculating Intensity Pattern..."
+n_eval = 500
+# Evaluate over a 30mm window to see many fringes
+xs_eval, zs_eval, I = intensity(target_pd; n = n_eval,
+    x_min = -15e-3, x_max = 15e-3, z_min = -15e-3, z_max = 15e-3)
+
+##
+# Visualization
+fig = Figure(size = (800, 600), fontsize = 18)
+display(fig)
+
+# Axis 1: 3D Scene (The Setup)
+ax3d = LScene(fig[1:2, 1], show_axis = false)
+
+# We create a mesh grid for the surface to ensure correct orientation in 3D
+X_mask = [xv for xv in x, zv in z]
+Z_mask = [zv for xv in x, zv in z]
+Y_mask = fill(0.0, size(X_mask))
+surface!(ax3d, X_mask, Y_mask, Z_mask,
+    color = amplitude, colormap = :binary, transparency = true)
+
+# Render Intensity Results (at y = L)
+X_res = [xv for xv in xs_eval, zv in zs_eval]
+Z_res = [zv for xv in xs_eval, zv in zs_eval]
+Y_res = fill(L, size(X_res))
+I_norm = I ./ (maximum(I) + 1e-12) # Avoid division by zero
+plt_res = surface!(ax3d, X_res, Y_res, Z_res,
+    color = I_norm, colormap = :magma,
+    transparency = false,
+    shading = NoShading,
+    depth_shift = -1e-3) # Subtle shift to win Z-fighting without clipping
+
+# Render a few representative beamlets
+step = max(1, length(beams) ÷ 15)
+for i in 1:step:length(beams)
+    render!(ax3d, beams[i], color = (:cyan, 0.05), show_beams = false, flen = 0.0)
+end
+
+# Render the system components (Detector frame)
+#render!(ax3d, system)
+
+# Axis 2: 2D Intensity Map (Zoomed)
+ax2d = Axis(fig[1, 2],
+    title = "Diffraction Pattern (Detector Plane)",
+    xlabel = "x [mm]", ylabel = "z [mm]", yticks = -1:2:1,
+    aspect = DataAspect())
+hm = heatmap!(ax2d, xs_eval .* 1000, zs_eval .* 1000, I, colormap = :magma)
+ylims!(-1.1, 1.1)
+
+# Axis 3: 1D Cross-section
+ax1d = Axis(fig[2, 2],
+    title = "Interference Fringes (x-slice)",
+    xlabel = "x Position [mm]", ylabel = "Intensity",
+    yticksvisible = false, yticklabelsvisible = false)
+lines!(ax1d, xs_eval .* 1000, I[:, size(I, 2) ÷ 2], color = :red, linewidth = 2)
+
+# This makes the 1.2m distance look shorter so we can zoom in on X/Z details
+scale!(ax3d.scene, 1.0, 0.15, 1.0)
+
+# Set a nice default camera view
+cl = ([0.005679703099338854, -0.010364784764641269, 0.002815435213044604],
+    [-0.003708494535423113, 0.026029085430702627, 0.0009259786679933707],
+    [-0.021598444746385476, 0.046277584411491254, 0.998695114799569])
+
+rowsize!(fig.layout, 1, Relative(1 // 4))
+
+set_view(ax3d, cl...)
+save("agb_doubleslit_experiment.png", fig; px_per_unit = 8, update = false)