Handle two ring lenses (#7)

* Handle two ring lenses

* Improve ring code

* Implement a bunch of fixes to lens rings...

* Remove forgotten debug statements

* Add ring test cases to avoid regression

* up vn to 0.8.3

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.8.2"
+version = "0.8.3"
 
 
 [deps]

src/OpticalComponents/Lenses.jl

@@ -184,15 +184,12 @@ function Lens(
 
     # Initialize remaining cylindrical section length.
     l0 = center_thickness
-
     # Front Surface
     front = sdf(front_surface, ForwardOrientation())
     l0 -= isnothing(front) ? zero(l0) : thickness(front)
-
     # Back Surface
     back = sdf(back_surface, BackwardOrientation())
     l0 -= isnothing(back) ? zero(l0) : thickness(back)
-
     # Use MeniscusLensSDF if cylinder length is non-positive
     if l0 ≤ 0
         if sign(radius(front_surface)) == sign(radius(back_surface))
@@ -222,39 +219,89 @@ function Lens(
         end
         shape = mid
 
-        # Add mechanical ring if md_mid > d_mid
-        if md_mid > d_mid
-            ring_thickness = thickness(mid)
-            ring_center = position(mid)[2] + ring_thickness / 2
-            if front !== nothing
-                s = edge_sag(front_surface, front)
-                ring_thickness -= s
-                ring_center += s / 2
+        d_front = diameter(front_surface)
+        d_back = diameter(back_surface)
+        d_min = min(d_front, d_back)
+        d_max = max(d_front, d_back)
+
+        if md_mid < d_min
+            @warn "Mechanical diameter is less than clear aperture; parameter md has been ignored."
+        else
+            # add leveling ring if there is a step between front and back clear apertures
+            if d_front != d_back
+                if d_back > d_front
+                    # Step exists on the front side: level front to match back.
+                    leveling_thickness = l0
+                    if front !== nothing
+                        s_front = edge_sag(front_surface, front)
+                        if s_front < 0
+                            # edge curves towards negative so it is a concave type shape,
+                            # which has to be covered by the ring
+                            leveling_thickness += abs(s_front) + thickness(front)
+                        end
+                    end
+
+                    ring = RingSDF(d_front / 2, (d_back - d_front) / 2, leveling_thickness)
+                    translate3d!(ring, [0, edge_sag(front_surface, front) + leveling_thickness / 2, 0])
+                    shape += ring
+                else  # d_front > d_back
+                    # Step exists on the back side: level back to match front.
+                    leveling_thickness = l0
+                    if back !== nothing
+                        s_back = edge_sag(back_surface, back)
+                        if (s_back - thickness(back)) > 0
+                            # edge curves towards positive so it is a concave type shape,
+                            # which has to be covered by the ring
+                            leveling_thickness += abs(s_back) + thickness(back)
+                        end
+                    end
+                    leveling_center = position(mid)[2] +
+                                      (l0/2 + (back !== nothing ? edge_sag(back_surface, back) : 0))
+                    if back !== nothing
+                        leveling_center += s_back / 2
+                    end
+                    ring = RingSDF(d_back / 2, (d_front - d_back) / 2, leveling_thickness)
+                    translate3d!(ring, [0, thickness(front) + leveling_thickness/2, 0])
+                    shape += ring
+                end
             end
-            if back !== nothing
-                s = edge_sag(back_surface, back)
-                ring_thickness += s
-                ring_center += s / 2
+
+            # add outer ring if the mechanical diameter exceeds the larger clear aperture.
+            if md_mid > d_max
+                outer_thickness = thickness(mid)
+                outer_center = position(mid)[2] + outer_thickness / 2
+                if front !== nothing
+                    s_front = edge_sag(front_surface, front)
+                    outer_thickness -= s_front
+                    outer_center += s_front / 2
+                end
+                if back !== nothing
+                    s_back = edge_sag(back_surface, back)
+                    outer_thickness += s_back
+                    outer_center += s_back / 2
+                end
+                ring = RingSDF(d_max / 2, (md_mid - d_max) / 2, outer_thickness)
+                translate3d!(ring, [0, outer_center, 0])
+                shape += ring
             end
-            ring = RingSDF(d_mid / 2, (md_mid - d_mid) / 2, ring_thickness)
-            translate3d!(ring, [0, ring_center, 0])
-            shape += ring
-        elseif md_mid < d_mid
-            @warn "Mechanical diameter is less than clear aperture; parameter md has been ignored."
         end
     end
 
     return Lens(shape, n)
 end
 
-Lens(front_surface::AbstractRotationallySymmetricSurface, center_thickness::Real, n::RefractiveIndex) = Lens(
-    front_surface,
-    CircularFlatSurface(diameter(front_surface)),
-    center_thickness,
-    n
-)
+function Lens(front_surface::AbstractRotationallySymmetricSurface,
+        center_thickness::Real, n::RefractiveIndex)
+    Lens(
+        front_surface,
+        CircularFlatSurface(diameter(front_surface)),
+        center_thickness,
+        n
+    )
+end
 
-function Lens(front_surface::CircularFlatSurface, back_surface::CircularFlatSurface, center_thickness::Real, n::RefractiveIndex)
+function Lens(front_surface::CircularFlatSurface, back_surface::CircularFlatSurface,
+        center_thickness::Real, n::RefractiveIndex)
     d_mid = min(diameter(front_surface), diameter(back_surface))
 
     return Lens(
@@ -306,7 +353,7 @@ function Lens(
     else
         # Construct central box and add front/back surfaces
         mid = BoxSDF(h, l0, d_mid)
-        translate3d!(mid, [0, l0/2, 0])
+        translate3d!(mid, [0, l0 / 2, 0])
         if front !== nothing
             translate3d!(mid, [0, thickness(front), 0])
             mid += front
@@ -342,15 +389,20 @@ function Lens(
     return Lens(shape, n)
 end
 
-Lens(front_surface::AbstractCylindricalSurface, center_thickness::Real, n::RefractiveIndex) = Lens(
-    front_surface,
-    RectangularFlatSurface(diameter(front_surface)),
-    center_thickness,
-    n
-)
+function Lens(front_surface::AbstractCylindricalSurface,
+        center_thickness::Real, n::RefractiveIndex)
+    Lens(
+        front_surface,
+        RectangularFlatSurface(diameter(front_surface)),
+        center_thickness,
+        n
+    )
+end
 
-function cylindric_lens_outer_parameters(f::AbstractCylindricalSurface, b::AbstractCylindricalSurface)
-    height(f) != height(b) && throw(ArgumentError("height of front and back surface have to match for cylindric lenses"))
+function cylindric_lens_outer_parameters(
+        f::AbstractCylindricalSurface, b::AbstractCylindricalSurface)
+    height(f) != height(b) &&
+        throw(ArgumentError("height of front and back surface have to match for cylindric lenses"))
 
     d_mid = min(diameter(f), diameter(b))
     md_mid = max(mechanical_diameter(f), mechanical_diameter(b))
@@ -359,27 +411,33 @@ function cylindric_lens_outer_parameters(f::AbstractCylindricalSurface, b::Abstr
     return d_mid, md_mid, h
 end
 
-function cylindric_lens_outer_parameters(f::AbstractCylindricalSurface, ::RectangularFlatSurface)
+function cylindric_lens_outer_parameters(
+        f::AbstractCylindricalSurface, ::RectangularFlatSurface)
     d_mid = diameter(f)
     md_mid = mechanical_diameter(f)
     h = height(f)
 
     return d_mid, md_mid, h
 end
 
-cylindric_lens_outer_parameters(f::RectangularFlatSurface, b::AbstractCylindricalSurface) = cylindric_lens_outer_parameters(b, f)
+function cylindric_lens_outer_parameters(
+        f::RectangularFlatSurface, b::AbstractCylindricalSurface)
+    cylindric_lens_outer_parameters(b, f)
+end
 
-function cylindric_lens_outer_parameters(f::RectangularFlatSurface, b::RectangularFlatSurface)
+function cylindric_lens_outer_parameters(
+        f::RectangularFlatSurface, b::RectangularFlatSurface)
     d_mid = diameter(f)
     md_mid = mechanical_diameter(f)
 
     return d_mid, md_mid, d_mid
 end
 
-function Lens(front_surface::RectangularFlatSurface, back_surface::RectangularFlatSurface, center_thickness::Real, n::RefractiveIndex)
+function Lens(front_surface::RectangularFlatSurface, back_surface::RectangularFlatSurface,
+        center_thickness::Real, n::RefractiveIndex)
     d_mid = min(diameter(front_surface), diameter(back_surface))
     mid = BoxSDF(d_mid, center_thickness, d_mid)
-    translate3d!(mid, [0, center_thickness/2, 0])
+    translate3d!(mid, [0, center_thickness / 2, 0])
 
     return Lens(
         mid,

src/SDFs/AsphericalLensSDF.jl

@@ -27,7 +27,8 @@ mutable struct ConvexAsphericalSurfaceSDF{T} <: AbstractAsphericalSurfaceSDF{T}
 end
 
 function thickness(s::ConvexAsphericalSurfaceSDF)
-    abs(aspheric_equation(s.diameter / 2, 1 / s.radius, s.conic_constant, s.coefficients))
+    sag = aspheric_equation(s.diameter / 2, 1 / s.radius, s.conic_constant, s.coefficients)
+    return s.max_sag[1] > 0 && sag < 0 ? s.max_sag[1] : abs(sag)
 end
 
 # Constructor for ConvexAsphericalSurfaceSDF

src/SDFs/SphericalLensSDF.jl

@@ -465,7 +465,7 @@ end
 
 mechanical_diameter(s::SphericalSurface) = s.mechanical_diameter
 
-edge_sag(::SphericalSurface, sd::AbstractSphericalSurfaceSDF) = abs(sag(sd))
+edge_sag(::SphericalSurface, sd::AbstractSphericalSurfaceSDF) = sag(sd)
 
 function sdf(s::SphericalSurface, ot::AbstractOrientationType)
     isinf(radius(s)) && return nothing

test/runtests.jl

@@ -1316,6 +1316,39 @@ end
 
         # test lens thickness
         @test BMO.thickness(lens) ≈ l
+
+        # test ring generation
+        NBK7 = DiscreteRefractiveIndex([532e-9, 1064e-9], [1.5195, 1.5066])
+        mm = 1e-3
+
+        s1 = Lens(
+            SphericalSurface(38.184mm, 2*1.840mm, 2*2.380mm),
+            SphericalSurface(3.467mm, 2*2.060mm, 2*2.380mm),
+            0.5mm,
+            NBK7
+        )
+
+        @test 2*s1.shape.sdfs[5].hthickness ≈ 0.001134 atol=1e-6
+
+        s2 = Lens(
+            SphericalSurface(3.467mm, 2*2.060mm, 2*2.380mm),
+            SphericalSurface(-5.020mm, 2*2.380mm, 2*2.380mm),
+            2.5mm,
+            NBK7
+        )
+
+        @test 2*s2.shape.sdfs[4].hthickness ≈ 0.001221590 atol=1e-8
+
+        # doublet test case
+        s1 = SphericalSurface(7.744mm, 2*2.812mm, 2*3mm)
+        s2 = SphericalSurface(-3.642mm, 2*3mm)
+        s3 = SphericalSurface(-14.413mm, 2*2.812mm, 2*3mm)
+
+        dl21 = Lens(s1, s2, 3.4mm, NBK7)
+        dl22 = Lens(s2, s3, 1.0mm, NBK7)
+
+        @test 2*dl21.shape.sdfs[4].hthickness ≈ 0.001294398 atol=1e-6
+        @test 2*dl22.shape.sdfs[4].hthickness ≈ 0.000723025 atol=1e-6
     end
 end
 
@@ -1486,6 +1519,11 @@ end
             # test if the beam is correctly focussed
             @test f_pos[3]≈0 atol=1e-7
         end
+
+        # test rings
+        @test 2*L1.shape.sdfs[4].hthickness ≈ 0.00060839 atol=1e-6
+        @test 2*L2.shape.sdfs[4].hthickness ≈ 0.00057497 atol=1e-6
+        @test 2*L3.shape.sdfs[4].hthickness ≈ 0.00048395 atol=1e-6
     end
 end