Point Spread Functions (#10)

* Implement plane wave interaction with Photodetector

* Fix review points

* add point source

* add CollimatedSource

* up vn to 0.8.5

* fix spot detector docs

* export reset_detector!

* add pd interact3d comments

* PSFDetector (warning: WIP)

* reset_photodetector! -> reset! and optimize PSF code for performance

* Fix syntax bugs

* spotdetector reset! adjust

* psf intensity bug fix

* opl docs fix

* rename reset! to empty! (dispatch base)

* up vn to 0.9 for breaking name change

* stricter refraction3d testset for (very) small angles

* add auto-limits (WIP) to PSFDetector

* Get rid of most allocations

* Use Auto-Diff for SDF gradients to fix issue #11

* Remove temporary issue fix

* add testset for jesus bug

* add BeamGroup testset

* rename Groups.jl to ObjectGroups.jl

* beam group docs

* fix typo

* fix docstring links

* rm create_spot_diagram

* Add primitive means of Strehl normalization

* Improve strehl ratio and centering

* Slay issue #11 again (happy Easter! It should have stayed dead) and fix specialization issues

* Remove photodetector PSF implementation

* Remove abs, match the sdf function behaviour

* Fix axis labeling and introduce a shift option for the user

* Adds the possibility to set a custom grid limit

* Remove Strehl, some docs

* beam group docs

* extended beam group tests

* typo

* minor changes 1

* minor changes 2

* import position

* export psfd (again)

* not a software

* fix some merge issues with regards to PR9

* Add a testset for PSFs and the UniformDiscSource

* Update UniformDiscSource doc strings

* Finish docs of UniformDiscSource

* Add docs

* rm BeamletOptics namespace when unnecessary

* rm solved field from psf detector (no purpose)

* typos

* UDSource docs update

* test psfd empty!

* fix num_rays kwarg in collimated sc

Author: TacHawkes

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.10.0"
+version = "0.10.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -44,4 +44,4 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test"]
+test = ["Aqua", "Test"]

docs/src/assets/beam_renders/collimated_sc.jl

@@ -25,3 +25,23 @@ cs_view = [
 
 set_view(cs_ax, cs_view)
 save("collimated_beam_source.png", cs_fig; px_per_unit=8, update = false)
+
+## uniform disc source
+cs = UniformDiscSource([0,0,0], [0,1,0], 15e-3; num_rays=100)
+
+cs_fig = Figure(; size=(600,200))
+display(cs_fig)
+cs_ax = LScene(cs_fig[1,1])
+hide_axis(cs_ax)
+
+render!(cs_ax, cs, show_pos=true, flen=0.05, color=:red, render_every=1)
+
+cs_view = [
+    0.26513     0.964213    7.68482e-16  -0.0180273
+    0.0488896  -0.0134432   0.998714      0.00108822
+    0.962972   -0.264789   -0.0507042    -0.025119
+    0.0         0.0         0.0           1.0
+]
+
+set_view(cs_ax, cs_view)
+save("collimated_uniform_beam_source.png", cs_fig; px_per_unit=8, update = false)

docs/src/assets/mi_assets/michelson_plots.jl

@@ -127,13 +127,13 @@ heat2 = Axis(fringes_fig[1, 2], xlabel="x [mm]", ylabel="y [mm]", title="After r
 hidedecorations!(heat1)
 hidedecorations!(heat2)
 
-hm = heatmap!(heat1, pd.x*1e3, pd.y*1e3, BeamletOptics.intensity(pd), colormap=:viridis)
+hm = heatmap!(heat1, pd.x*1e3, pd.y*1e3, intensity(pd), colormap=:viridis)
 
 zrotate3d!(m1, 1e-3)
 empty!(pd)
 solve_system!(system, beam)
 
-hm = heatmap!(heat2, pd.x*1e3, pd.y*1e3, BeamletOptics.intensity(pd), colormap=:viridis)
+hm = heatmap!(heat2, pd.x*1e3, pd.y*1e3, intensity(pd), colormap=:viridis)
 
 save("mi_fringes.png", fringes_fig, px_per_unit=4)
 

docs/src/assets/photodetector_showcase.jl

@@ -21,7 +21,7 @@ solve_system!(system, g2)
 ## render fringes
 fringes_fig = Figure()
 heat = Axis(fringes_fig[1, 1], xlabel="x [mm]", ylabel="y [mm]", aspect=1)
-hm = heatmap!(heat, pd.x*1e3, pd.y*1e3, BeamletOptics.intensity(pd), colormap=:viridis)
+hm = heatmap!(heat, pd.x*1e3, pd.y*1e3, intensity(pd), colormap=:viridis)
 cb = Colorbar(fringes_fig[1, 2], hm, label="Intensity [W/m²]")
 
 ## render system

docs/src/assets/psfdetector_showcase.jl

@@ -0,0 +1,59 @@
+using GLMakie, BeamletOptics
+
+const mm = 1e-3
+
+l = 1e-3
+R1 = 100e-3
+R2 = Inf
+d = 25.4e-3
+n = 1.5
+λ = 1e-6
+
+cs = UniformDiscSource([0, -10mm, 0], [0, 1, 0], 15e-3, λ)
+lens = SphericalLens(R1, R2, l, d, x -> n)
+
+psfd = PSFDetector(10e-3)
+translate3d!(psfd, [0, 200mm + 0.13mm, 0])
+
+sys = System([lens, psfd])
+
+solve_system!(sys, cs)
+
+x, z, I_num = intensity(psfd; n=500, crop_factor=10)
+
+fringes_fig = Figure()
+srf = Axis3(fringes_fig[1, 1], xlabel="x [µm]", ylabel="z [µm]", elevation=0.9, azimuth=2.372)
+hm = surface!(srf, x*1e6, z*1e6, I_num, colormap=:viridis)
+
+hidezdecorations!(srf)
+
+##
+k = -0.675
+d = 75.0e-3
+l = 15e-3
+radius = 76.68e-3
+A = [0*(1e3)^1, 2.7709219e-8*(1e3)^3, 6.418186e-13*(1e3)^5, -1.5724014e-17*(1e3)^7, -2.7768768e-21*(1e3)^9, -2.590162e-25*(1e3)^11]
+AL75150 = Lens(
+    EvenAsphericalSurface(radius, d, k, A),
+    l,
+    n -> 1.5006520430
+)
+
+xrotate3d!(AL75150, deg2rad(-0.5))
+
+pd = PSFDetector(15e-3)
+
+translate3d!(pd, [0, 158.1779e-3, 0.0])
+system = System([AL75150, pd])
+
+ps = UniformDiscSource([0, -0.1, 0], [0,1,0], 0.8*d, 1550e-9)
+
+solve_system!(system, ps)
+
+asph_fig = Figure(size=(800, 500))
+
+xs, ys, _intensity = intensity(pd, n=801, x_min=-25e-6, x_max=25e-6, z_min=-25e-6, z_max=25e-6, z0_shift=-50e-6)
+
+ax = Axis(asph_fig[1,1], xlabel="x [µm]", ylabel="z [µm]", aspect=1)
+
+hm = heatmap!(xs*1e6, ys*1e6, _intensity)

docs/src/assets/spotdetector_showcase.jl

@@ -31,9 +31,7 @@ cs = CollimatedSource([0,-50mm,0], [0,1,0], aperture, 1e-6; num_rings, num_rays)
 
 t1 = @timed solve_system!(system, cs)
 
-for i = 1:50:length(BeamletOptics.beams(cs))
-    render!(system_ax, BeamletOptics.beams(cs)[i], color=:blue, show_pos=true)
-end
+render!(system_ax, cs, color=:blue, show_pos=true, render_every=50)
 
 ## render diagram
 spot_fig = Figure(size=(600,400))

docs/src/basics/beams.md

@@ -46,6 +46,16 @@ CollimatedSource(::AbstractArray{<:Real}, ::AbstractArray{<:Real}, ::Real, ::Rea
 
 ![Collimated group of beams](collimated_beam_source.png)
 
+A special constructor called [`UniformDiscSource`](@ref) is available, which offers an equal-area
+sampling (Fibonnaci-pattern) sampling and is thus favorable in situations where the weighting of the
+individual beams becomes important, e.g. for calculating a point spread function using [`PSFDetector`](@ref).
+
+```@docs; canonical=false
+UniformDiscSource
+```
+
+![Collimated uniform group of beams](collimated_uniform_beam_source.png)
+
 ### Point beam source
 
 The `PointSource` type is used to model emission from a spatially localized source that radiates [`Beam`](@ref)s in a range of directions. This is commonly used to simulate conical emission patterns, such as light emerging from a fiber tip or a light source for a lens objective with a known focal distance. You can specify the origin and a propagation direction, which are then used to construct the monochromatic `PointSource`.

docs/src/components/detectors.md

@@ -26,7 +26,7 @@ save("pd_showcase.png", detector_fig, px_per_unit=4)
 nothing
 ```
 
-The `interact3d` model of the [`Photodetector`](@ref) can store complex electric field (E-field) values from intersecting [`GaussianBeamlet`](@ref)s, enabling the reconstruction of spatial intensity distribution across its active surface. This data can be used to calculate e.g. beam interference patterns via the [`BeamletOptics.intensity`](@ref) function. The [`BeamletOptics.optical_power`](@ref) method can be used in order to obtain the total optical power at the detector. Below a rendered example of a detector model ([FDS010](https://www.thorlabs.com/thorproduct.cfm?partnumber=FDS010)) can be seen. The detector active area is marked in blue (1x1 mm²). 
+The `interact3d` model of the [`Photodetector`](@ref) can store complex electric field (E-field) values from intersecting [`GaussianBeamlet`](@ref)s, enabling the reconstruction of spatial intensity distribution across its active surface. This data can be used to calculate e.g. beam interference patterns via the [`intensity`](@ref) function. The [`BeamletOptics.optical_power`](@ref) method can be used in order to obtain the total optical power at the detector. Below a rendered example of a detector model ([FDS010](https://www.thorlabs.com/thorproduct.cfm?partnumber=FDS010)) can be seen. The detector active area is marked in blue (1x1 mm²). 
 
 ![Photodetector showcase](pd_showcase.png)
 
@@ -62,4 +62,116 @@ Below an optical system consisting of a collection of collimated [`Beam`](@ref)s
 
 The beam bundle used to generate the spot diagram was created via the [`CollimatedSource`](@ref) constructor. The resulting spot diagram of the lens shown above is visualized below.
 
-![Spot diagram showcase](spot_diagram_showcase.png)
+![Spot diagram showcase](spot_diagram_showcase.png)
+
+## Point-spread-function detector type
+
+!!! warning "Experimental feature"
+    The point spread function estimation is a highly experimental feature. It does not use
+    pupils (yet) but merely uses superposition of the ray-attached plane-waves. While this
+    gives qualitatively sound results, it requires good sampling of the problem to obtain
+    quantitatively good results. Currently no Strehl-ratio is calculated due to that.
+
+The package offers a simple method to estimate the point spread function of a system. It is 
+currently limited and requires careful assessment by the user, if the results are to be trusted.
+
+To analyze the PSF of a imaging system a [`PSFDetector`](@ref) is added to the system at the plane
+and orientation where the PSF is requested. This is the same approach as for the other detector types.
+
+```@docs; canonical=false
+PSFDetector(::Real)
+```
+
+The intensity map together with the coordinate system of the detector can be retrieved after solving the system by calling the [`intensity`](@ref) function.
+
+```@docs; canonical=false
+intensity(::PSFDetector)
+```
+
+!!! warning
+    When dealing with a collimated source as the input to your optical system, where you want to calculate the PSF, **DO NOT** use the [`CollimatedSource`](@ref) beam group directly but instead use the [`UniformDiscSource`](@ref) constructor. This function returns a `CollimatedSource` with an equal-area sampling, which correctly weights the outer beams in relation to the inner beams. Otherwise the results might be wrong.
+
+
+```@eval
+asset_dir = joinpath(@__DIR__, "..", "assets")
+
+Base.include(@__MODULE__, joinpath(asset_dir, "psfdetector_showcase.jl"))
+
+save("psf_airy_showcase.png", fringes_fig, px_per_unit=4)
+save("psf_tilted_showcase.png", asph_fig, px_per_unit=4)
+
+nothing
+```
+
+### Airy-Disc Example
+
+This is a classic example where a collimated circular beam is imaged onto a point by a singlet lens.
+Due to the finite size of the aperture stop (in this case given by the 15 mm size of the beam), the diffraction
+limited intensity pattern is given by the Airy-disc:
+
+```math
+I(r)=I_0\!\left[\frac{2J_1\!\bigl(\pi D r/(\lambda f)\bigr)}{\pi D r/(\lambda f)}\right]^2
+```
+
+With ``r`` the radius from the origin, ``I_0`` the maximum intensity, ``J_1`` the Bessel function of the first kind of order one, ``D`` the aperture width, ``\lambda`` the wavelength and the focal length ``f``.
+
+```julia
+# example parameters
+l = 1e-3
+R1 = 100e-3
+R2 = Inf
+d = 25.4e-3
+n = 1.5
+λ = 1e-6
+
+# generate uniform source, lens and PSF detector
+cs = UniformDiscSource([0, -10mm, 0], [0, 1, 0], 15e-3, λ)
+lens = SphericalLens(R1, R2, l, d, x -> n)
+psfd = PSFDetector(10e-3)
+
+# shift detector into focus
+translate3d!(psfd, [0, 200mm + 0.13mm, 0])
+
+# build system
+sys = System([lens, psfd])
+
+solve_system!(sys, cs)
+
+# retrieve intensity
+x, z, I_num = intensity(psfd; n=500, crop_factor=5)
+```
+
+Visualizing the result yields the expected Airy-disk pattern.
+
+![Airy disc PSF](psf_airy_showcase.png)
+
+### Coma and Astigmatism Example
+
+In this example, an aspheric lens images the collimated source onto a point but is tilted around the x-axis by 0.5 degrees.
+This results in aberrations distorting the stigmatic imaging and leading to coma and astigmatism.
+
+```julia
+k = -0.675
+d = 75.0e-3
+l = 15e-3
+radius = 76.68e-3
+A = [0*(1e3)^1, 2.7709219e-8*(1e3)^3, 6.418186e-13*(1e3)^5, -1.5724014e-17*(1e3)^7, -2.7768768e-21*(1e3)^9, -2.590162e-25*(1e3)^11]
+AL75150 = Lens(
+    EvenAsphericalSurface(radius, d, k, A),
+    l,
+    n -> 1.5006520430
+)
+
+xrotate3d!(AL75150, deg2rad(-0.5))
+
+pd = PSFDetector(15e-3)
+
+translate3d!(pd, [0, 158.1779e-3, 0.0])
+system = System([AL75150, pd])
+
+ps = UniformDiscSource([0, -0.1, 0], [0,1,0], 0.8*d, 1550e-9)
+
+solve_system!(system, ps)
+```
+
+![Tilted asphere PSF](psf_tilted_showcase.png)

docs/src/tutorials/michelson.md

@@ -216,14 +216,14 @@ fringes_fig = Figure()
 heat1 = Axis(fringes_fig[1, 1], aspect=1)
 heat2 = Axis(fringes_fig[1, 2], aspect=1)
 
-hm = heatmap!(heat1, pd.x, pd.y, BeamletOptics.intensity(pd), colormap=:viridis)
+hm = heatmap!(heat1, pd.x, pd.y, intensity(pd), colormap=:viridis)
 
 # rotate m1, reset pd field data, resolve system
 zrotate3d!(m1, 1e-3)
 empty!(pd)
 solve_system!(system, beam)
 
-hm = heatmap!(heat2, pd.x, pd.y, BeamletOptics.intensity(pd), colormap=:viridis)
+hm = heatmap!(heat2, pd.x, pd.y, intensity(pd), colormap=:viridis)
 ```
 
 By experimenting with different mirror angles, arm lengths, or beamsplitter properties, you can observe how interference fringes evolve and gain insights into the stability and sensitivity of the interferometric setup. This can be important to optimize alignment and achieve high contrast fringes.

src/BeamGroups.jl

@@ -45,7 +45,7 @@ The following inputs and arguments can be used to configure the [`PointSource`](
 - `num_rings`: number of concentric beam rings, default is 10
 - `num_rays`: total number of rays in the source, default is 100x num_rings
 
-!!! warning 
+!!! warning
     The orthogonal basis vectors for the beam generation are generated randomly.
 """
 function PointSource(
@@ -138,15 +138,15 @@ The following inputs and arguments can be used to configure the [`CollimatedSour
 
 - `pos`: center beam starting position
 - `dir`: center beam starting direction
-- `diameter`: outer beam bundle diameter in [m] 
+- `diameter`: outer beam bundle diameter in [m]
 - `λ = 1e-6`: wavelength in [m], default val. is 1000 nm
 
 ## Keyword Arguments
 
 - `num_rings`: number of concentric beam rings, default is 10
 - `num_rays`: total number of rays in the source, default is 100x num_rings
 
-!!! warning 
+!!! warning
     The orthogonal basis vectors for the beam generation are generated randomly.
 """
 function CollimatedSource(
@@ -192,4 +192,54 @@ function CollimatedSource(
         end
     end
     return CollimatedSource(beams, T(diameter))
-end
+end
+
+"""
+    UniformDiscSource(pos, dir, diameter, λ; num_rays=1_000)
+
+Generates a ray fan with *equal area per ray* across a circular pupil
+using the deterministic sunflower (Fibonacci) pattern.
+
+!!! note
+    This is merely a [`CollimatedSource`](@ref) constructor which uses Fibonacci sampling
+    instead of a linear grid.
+
+# Arguments
+
+The following inputs and arguments can be used to configure the underlying [`CollimatedSource`](@ref):
+
+## Inputs
+
+- `pos`: center beam starting position
+- `dir`: center beam starting direction
+- `diameter`: outer beam bundle diameter in [m]
+- `λ = 1e-6`: wavelength in [m]
+
+## Keyword Arguments
+
+- `num_rays=1000`: total number of rays in the source
+"""
+function UniformDiscSource(
+        pos::AbstractArray{P},
+        dir::AbstractArray{D1},
+        diameter::D2,
+        λ::L=1e-6;
+        # kwargs
+        num_rays::Int=1_000
+    ) where {P<:Real, D1<:Real, D2<:Real, L<:Real}
+    T  = promote_type(P, D1, D2, L)
+    R  = diameter / 2
+    φ0 = 2π / (1 + √5)         # golden angle
+    beams = Vector{Beam{T,Ray{T}}}(undef, num_rays)
+    # orthogonal basis in the pupil plane
+    e1 = normal3d(dir)
+    e2 = normalize(cross(dir, e1))
+    for k in 0:num_rays-1
+        ρ  = √((k + 0.5)/num_rays)     # equal-area radius
+        φ  = k * φ0
+        r  = R * ρ
+        x  =  r * cos(φ) * e1 + r * sin(φ) * e2
+        beams[k+1] = Beam(pos + x, dir, λ)
+    end
+    return CollimatedSource(beams, T(diameter))
+end