From 4305c771c71cba7592d1a417ba8d7fa8d4b6948f Mon Sep 17 00:00:00 2001 From: mattderstine Date: Sun, 10 Oct 2021 15:44:18 -0700 Subject: [PATCH 01/19] first step of changes to move Aspheric function to it's own Surface type --- src/Geometry/Primitives/AsphericSurface.jl | 220 +++++++++++++++++++++ src/Geometry/Primitives/Zernike.jl | 105 +++------- 2 files changed, 246 insertions(+), 79 deletions(-) create mode 100644 src/Geometry/Primitives/AsphericSurface.jl diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl new file mode 100644 index 000000000..a105581b1 --- /dev/null +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -0,0 +1,220 @@ +# MIT license +# Copyright (c) Microsoft Corporation. All rights reserved. +# See LICENSE in the project root for full license information. + + + +######################################################################################################### + + + +""" + AsphericSurface{T,N,Q} <: ParametricSurface{T,N} + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter, Zernike terms are normalized according to the `normradius` parameter. +`T` is the datatype, `N` is the dimensionality, and `Q` is the number of aspheric terms. + +The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space. +Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned. + + +```julia +AsphericSurface(semidiameter, radius, conic, aspherics=nothing, normradius = semidiameter,) +``` + +`aspherics` should be a vector containing the aspheric polynomial coefficients starting at A3 + +The sag is defined by the equation + +```math +z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir^{i} +``` + +where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` . +""" +struct AsphericSurface{T,N,Q} <: ParametricSurface{T,N} + semidiameter::T + aspherics::SVector{Q,T} + curvature::T + conic::T + normradius::T + boundingcylinder::Cylinder{T,N} + + function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T, aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} + @assert semidiameter > 0 + @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) + @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" + + acs = [] + if aspherics !== nothing + asphericTerms = aspherics[:][1] + minAsphericTerm = minimum(asphericTerms) + maxAsphericTerm = maximum(asphericTerms) + @assert minAsphericTerm > 2 "Aspheric Orders must be A3 or higher (A$minAsphericTerm)" + acs = zeros(T, maxAsphericTerm -2 ) + for (i, k) in aspherics + acs[i-2] = k + end + end + Q = length(acs) + new{T,3,P,Q}(semidiameter, SVector{P,Tuple{Int,Int,T}}(zcs), SVector{Q,T}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + end + +end +export AsphericSurface + +uvrange(::Type{AsphericSurface{T,N,Q}}) where {T<:Real,N,Q} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ + +semidiameter(z::AsphericSurface{T}) where {T<:Real} = z.semidiameter +halfsizeu(z::AsphericSurface{T}) where {T<:Real} = semidiameter(z) +halfsizev(z::AsphericSurface{T}) where {T<:Real} = semidiameter(z) + +boundingobj(z::AsphericSurface{T}) where {T<:Real} = z.boundingcylinder + +function point(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} + rad = z.semidiameter + r = ρ * rad + r2 = r^2 + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN) + end + h = z.curvature * r2 / (one(T) + sqrt(t)) + # sum aspheric + prod = r*r + for asp in aspherics + prod *= r + h += asp * prod + end + # sum zernike + return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) +end + +function partials(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} + rad = z.semidiameter + r = ρ * rad + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) + end + dhdρ = rad * z.curvature * r * sqrt(t) / t + # sum aspherics partial + prod = r + for (m,asp) in enumerate(aspherics) + prod *= r + dhdρ += rad * (m+2) * asp * prod #aspherics has order 3 and above + end + dhdϕ = zero(T) + cosϕ = cos(ϕ) + sinϕ = sin(ϕ) + pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) + pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) + return pu, pv +end + +function normal(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} + du, dv = partials(z, ρ, ϕ) + if ρ == zero(T) && norm(dv) == zero(T) + # in cases where there is no δϕ at ρ = 0 (i.e. anything which is rotationally symetric) + # then we get some big problems, hardcoding this case solves the problems + return SVector{3,T}(0, 0, 1) + end + return normalize(cross(du, dv)) +end + +function uv(z::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real,Q} + # avoid divide by zero for ForwardDiff + ϕ = NaNsafeatan(p[2], p[1]) + if p[1] == zero(T) && p[2] == zero(T) + ρ = zero(T) + else + ρ = sqrt(p[1]^2 + p[2]^2) / semidiameter(z) + end + return SVector{2,T}(ρ, ϕ) +end + +function onsurface(surf::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real,Q} + ρ, ϕ = uv(surf, p) + if ρ > one(T) + return false + else + surfpoint = point(surf, ρ, ϕ) + return samepoint(p[3], surfpoint[3]) + end +end + +function inside(surf::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real,Q} + ρ, ϕ = uv(surf, p) + if ρ > one(T) + return false + else + surfpoint = point(surf, ρ, ϕ) + return p[3] < surfpoint[3] + end +end + +######################################################################################################### + +# Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis. +function surfaceintersection(surf::AcceleratedParametricSurface{T,3,AsphericSurface{T,3,Q}}, r::AbstractRay{T,3}) where {T<:Real,Q} + cylint = surfaceintersection(surf.surface.boundingcylinder, r) + if cylint isa EmptyInterval{T} + return EmptyInterval(T) + else + if doesintersect(surf.triangles_bbox, r) || inside(surf.triangles_bbox, origin(r)) + surfint = triangulatedintersection(surf, r) + if !(surfint isa EmptyInterval{T}) + return intervalintersection(cylint, surfint) + end + end + # hasn't hit the surface + if lower(cylint) isa RayOrigin{T} && upper(cylint) isa Infinity{T} + if inside(surf.surface, origin(r)) + return Interval(RayOrigin(T), Infinity(T)) + else + return EmptyInterval(T) + end + # otherwise check that the intersection is underneath the surface + else + p = point(closestintersection(cylint, false)) + ρ, ϕ = uv(surf, p) + surfpoint = point(surf.surface, ρ, ϕ) + if p[3] < surfpoint[3] + return cylint # TODO!! UV (and interface) issues? + else + return EmptyInterval(T) + end + end + end +end + +function AcceleratedParametricSurface(surf::T, numsamples::Int = 17; interface::NullOrFresnel{S} = NullInterface(S)) where {S<:Real,N,T<:AsphericSurface{S,N}} + # Zernike uses ρ, ϕ uv space so need to modify extension of triangulation + a = AcceleratedParametricSurface(surf, triangulate(surf, numsamples, true, false, true, false), interface = interface) + emptytrianglepool!(S) + return a +end + +function BoundingBox(surf::AsphericSurface{T,3,Q}) where {T<:Real,Q} + xmin = -semidiameter(surf) + xmax = semidiameter(surf) + ymin = -semidiameter(surf) + ymax = semidiameter(surf) + # aspherics can be more comlpicated, so just take sum of all negative and all positives + amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in surf.aspherics) : zero(T) + amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in surf.aspherics) : zero(T) + zmin = amin + zmax = amax + # curvature only goes one way + q = one(T) - (one(T) + surf.conic) * surf.curvature^2 * surf.semidiameter^2 + if q < zero(T) + throw(ErrorException("The surface is invalid, no bounding box can be constructed")) + end + hmax = surf.curvature * surf.semidiameter^2 / (one(T) + sqrt(q)) + if hmax > zero(T) + zmax += hmax + else + zmin += hmax + end + return BoundingBox(xmin, xmax, ymin, ymax, zmin, zmax) +end diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 3bb0180ff..b63f5915e 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -130,7 +130,7 @@ For convenience the input `zcoeff` can be indexed using either OSA or Noll conve ZernikeSurface(semidiameter, radius = Inf, conic = 0, zcoeff = nothing, aspherics = nothing, normradius = semidiameter, indexing = ZernikeIndexingOSA) ``` -`zcoeff` and `aspherics` should be vectors containing tuples of the form `(i, v)` where `i` is either the index of the Zernike term for the corresponding `indexing`, or the polynomial power of the aspheric term (must be even) and `v` is the corresponding coefficient ``A_i`` or ``\\alpha_i`` respectively.. +`zcoeff` and `aspherics` should be vectors containing tuples of the form `(i, v)` where `i` is either the index of the Zernike term for the corresponding `indexing`, or the polynomial power of the aspheric term (may be even or odd) and `v` is the corresponding coefficient ``A_i`` or ``\\alpha_i`` respectively.. The sag is defined by the equation @@ -141,18 +141,20 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Z_n`` is the nᵗʰ Zernike polynomial. """ struct ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} - semidiameter::T + asp::AsphericSurface{T,N,Q} + #semidiameter::T coeffs::SVector{P,Tuple{Int,Int,T}} - aspherics::SVector{Q,Tuple{Int,T}} - curvature::T - conic::T - normradius::T + #aspherics::SVector{Q,Tuple{Int,T}} + #curvature::T + #conic::T + #normradius::T boundingcylinder::Cylinder{T,N} function ZernikeSurface(semidiameter::T; radius::T = typemax(T), conic::T = zero(T), zcoeff::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter, indexing::ZernikeIndexType = ZernikeIndexingOSA) where {T<:Real} - @assert semidiameter > 0 - @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) - @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" + asp = AsphericSurface(semidiameter; radius, conic, aspherics, normradius) + #@assert semidiameter > 0 + #@assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) + #@assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" zcs = [] if zcoeff !== nothing for (i, k) in zcoeff @@ -167,17 +169,8 @@ struct ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} end end P = length(zcs) - acs = [] - if aspherics !== nothing - for (i, k) in aspherics - if abs(k) > zero(T) - @assert i % 2 == 0 "Only even aspheric terms are supported" - push!(acs, (i ÷ 2, k)) - end - end - end - Q = length(acs) - new{T,3,P,Q}(semidiameter, SVector{P,Tuple{Int,Int,T}}(zcs), SVector{Q,Tuple{Int,T}}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + + new{T,3,P,Q}(asp, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end end @@ -185,62 +178,35 @@ export ZernikeSurface uvrange(::Type{ZernikeSurface{T,N,P,Q}}) where {T<:Real,N,P,Q} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ -semidiameter(z::ZernikeSurface{T}) where {T<:Real} = z.semidiameter +semidiameter(z::ZernikeSurface{T}) where {T<:Real} = z.asp.semidiameter halfsizeu(z::ZernikeSurface{T}) where {T<:Real} = semidiameter(z) halfsizev(z::ZernikeSurface{T}) where {T<:Real} = semidiameter(z) boundingobj(z::ZernikeSurface{T}) where {T<:Real} = z.boundingcylinder function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q} - rad = z.semidiameter - r = ρ * rad - r2 = r^2 - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN) - end - h = z.curvature * r2 / (one(T) + sqrt(t)) - # sum aspheric - @inbounds @simd for i in 1:Q - (m, k) = z.aspherics[i] - h += k * r2^m - end + pnt = point(z.asp, ρ, ϕ) + # sum zernike u = r / z.normradius @inbounds @simd for m in 1:P (R, S, k) = z.coeffs[m] - h += k * Zernike.ζ(R, S, u, ϕ) + pnt[3] += k * Zernike.ζ(R, S, u, ϕ) end - return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) + return pnt #pnt is already the right type and length end function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q} - rad = z.semidiameter - r = ρ * rad - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) - end - dhdρ = rad * z.curvature * r * sqrt(t) / t - # sum aspherics partial - @inbounds @simd for i in 1:Q - (m, k) = z.aspherics[i] - dhdρ += rad * (2 * m) * k * r^(2 * m - 1) - end + pu,pv = partials(z.asp, ρ, ϕ) # sum zernike partials - dhdϕ = zero(T) - n = rad / z.normradius + n = rad / z.asp.normradius u = ρ * n @inbounds @simd for m in 1:P (R, S, k) = z.coeffs[m] du, dϕ = Zernike.δζ(R, S, u, ϕ) - dhdρ += k * du * n # want the derivative wrt ρ, not u - dhdϕ += k * dϕ + pu[3] += k * du * n # want the derivative wrt ρ, not u + pv[3] += k * dϕ end - cosϕ = cos(ϕ) - sinϕ = sin(ϕ) - pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) - pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) return pu, pv end @@ -328,30 +294,11 @@ function AcceleratedParametricSurface(surf::T, numsamples::Int = 17; interface:: end function BoundingBox(surf::ZernikeSurface{T,3,P,Q}) where {T<:Real,P,Q} - xmin = -semidiameter(surf) - xmax = semidiameter(surf) - ymin = -semidiameter(surf) - ymax = semidiameter(surf) + bb = BoundingBox(z.asp) # zernike terms have condition than |Zᵢ| <= 1 # so this gives us a (loose) bounding box ak = P > 0 ? sum(abs.(Zernike.normalisation(T, N, M) * k for (N, M, k) in surf.coeffs)) : zero(T) - zmin = -ak - zmax = ak - # aspherics can be more comlpicated, so just take sum of all negative and all positives - amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in surf.aspherics) : zero(T) - amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in surf.aspherics) : zero(T) - zmin += amin - zmax += amax - # curvature only goes one way - q = one(T) - (one(T) + surf.conic) * surf.curvature^2 * surf.semidiameter^2 - if q < zero(T) - throw(ErrorException("The surface is invalid, no bounding box can be constructed")) - end - hmax = surf.curvature * surf.semidiameter^2 / (one(T) + sqrt(q)) - if hmax > zero(T) - zmax += hmax - else - zmin += hmax - end - return BoundingBox(xmin, xmax, ymin, ymax, zmin, zmax) + bb.zmin -= ak + bb.zmax += ak + return BoundingBox(bb.xmin, bb.xmax, bb.ymin, bb.ymax, bb.zmin, bb.zmax) #could just return bb, but this way is safer end From 8007032cc99a94b357b9985922893c757772f3da Mon Sep 17 00:00:00 2001 From: mattderstine Date: Sun, 10 Oct 2021 15:58:07 -0700 Subject: [PATCH 02/19] Fix Geometry.jl to include AsphericSurface.jl --- src/Geometry/Geometry.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl index 011158618..9f7dfab68 100644 --- a/src/Geometry/Geometry.jl +++ b/src/Geometry/Geometry.jl @@ -51,6 +51,7 @@ include("Primitives/Curves/PowerBasis.jl") include("CSG/CSG.jl") +include("Primitives/AsphericSurface.jl") include("Primitives/Zernike.jl") include("Primitives/Qtype.jl") include("Primitives/Chebyshev.jl") From 9d5ab4f204cbd2d50036e88e506e7a11a01e4d12 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Sun, 10 Oct 2021 21:35:56 -0700 Subject: [PATCH 03/19] change to include A1 and up --- src/Geometry/Primitives/AsphericSurface.jl | 23 +++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index a105581b1..49429b40e 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -47,17 +47,17 @@ struct AsphericSurface{T,N,Q} <: ParametricSurface{T,N} acs = [] if aspherics !== nothing - asphericTerms = aspherics[:][1] + asphericTerms = [i for (i, ) in aspherics] minAsphericTerm = minimum(asphericTerms) maxAsphericTerm = maximum(asphericTerms) - @assert minAsphericTerm > 2 "Aspheric Orders must be A3 or higher (A$minAsphericTerm)" - acs = zeros(T, maxAsphericTerm -2 ) + @assert minAsphericTerm > 0 "Aspheric Terms must be A1 or higher (A$minAsphericTerm)" + acs = zeros(T, maxAsphericTerm ) for (i, k) in aspherics - acs[i-2] = k + acs[i] = k end end Q = length(acs) - new{T,3,P,Q}(semidiameter, SVector{P,Tuple{Int,Int,T}}(zcs), SVector{Q,T}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + new{T,3,Q}(semidiameter, SVector{P,Tuple{Int,Int,T}}(zcs), SVector{Q,T}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end end @@ -81,12 +81,11 @@ function point(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<: end h = z.curvature * r2 / (one(T) + sqrt(t)) # sum aspheric - prod = r*r + prod = one(T) for asp in aspherics prod *= r h += asp * prod end - # sum zernike return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) end @@ -99,10 +98,12 @@ function partials(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},S end dhdρ = rad * z.curvature * r * sqrt(t) / t # sum aspherics partial - prod = r - for (m,asp) in enumerate(aspherics) - prod *= r - dhdρ += rad * (m+2) * asp * prod #aspherics has order 3 and above + ((m, asp), rest) = peel(enumerate(aspherics)) + dhdρ += rad * asp #first term m=1 and prod = one(T) + prod = one(T) + for (m,asp) in rest + prod *= r + dhdρ += rad * m * asp * prod end dhdϕ = zero(T) cosϕ = cos(ϕ) From cfb7bb76fa0f9090cfe30524731d62c0662577f9 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Wed, 13 Oct 2021 16:00:05 -0700 Subject: [PATCH 04/19] Working AsphericSurface type, tests etc --- src/Geometry/Primitives/AsphericSurface.jl | 26 ++++++------ src/Geometry/Primitives/Zernike.jl | 27 +++++++------ test/TestData/lenses.jl | 20 ++++++++-- test/TestData/systems.jl | 4 +- test/testsets/Comparison.jl | 46 +++++++++++++++++++++- test/testsets/Intersection.jl | 1 - 6 files changed, 95 insertions(+), 29 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 49429b40e..38e2c004d 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -11,7 +11,7 @@ """ AsphericSurface{T,N,Q} <: ParametricSurface{T,N} -Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter, Zernike terms are normalized according to the `normradius` parameter. +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter,. `T` is the datatype, `N` is the dimensionality, and `Q` is the number of aspheric terms. The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space. @@ -19,10 +19,10 @@ Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defi ```julia -AsphericSurface(semidiameter, radius, conic, aspherics=nothing, normradius = semidiameter,) +AsphericSurface(semidiameter; radius, conic, aspherics=nothing, normradius = semidiameter) ``` -`aspherics` should be a vector containing the aspheric polynomial coefficients starting at A3 +`aspherics` should be a vector containing the aspheric polynomial coefficients starting at A1 The sag is defined by the equation @@ -40,7 +40,7 @@ struct AsphericSurface{T,N,Q} <: ParametricSurface{T,N} normradius::T boundingcylinder::Cylinder{T,N} - function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T, aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} + function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T= zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} @assert semidiameter > 0 @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" @@ -57,7 +57,7 @@ struct AsphericSurface{T,N,Q} <: ParametricSurface{T,N} end end Q = length(acs) - new{T,3,Q}(semidiameter, SVector{P,Tuple{Int,Int,T}}(zcs), SVector{Q,T}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + new{T,3,Q}(semidiameter, SVector{Q,T}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end end @@ -82,7 +82,7 @@ function point(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<: h = z.curvature * r2 / (one(T) + sqrt(t)) # sum aspheric prod = one(T) - for asp in aspherics + for asp in z.aspherics prod *= r h += asp * prod end @@ -98,12 +98,14 @@ function partials(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},S end dhdρ = rad * z.curvature * r * sqrt(t) / t # sum aspherics partial - ((m, asp), rest) = peel(enumerate(aspherics)) - dhdρ += rad * asp #first term m=1 and prod = one(T) - prod = one(T) - for (m,asp) in rest - prod *= r - dhdρ += rad * m * asp * prod + if length(z.aspherics) != 0 + ((m, asp), rest) = Iterators.peel(enumerate(z.aspherics)) + dhdρ += rad * asp #first term m=1 and prod = one(T) + prod = one(T) + for (m,asp) in rest + prod *= r + dhdρ += rad * m * asp * prod + end end dhdϕ = zero(T) cosϕ = cos(ϕ) diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index b63f5915e..658c335ba 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -140,7 +140,7 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Z_n`` is the nᵗʰ Zernike polynomial. """ -struct ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} +struct ZernikeSurface{T,N,P, Q} <: ParametricSurface{T,N} asp::AsphericSurface{T,N,Q} #semidiameter::T coeffs::SVector{P,Tuple{Int,Int,T}} @@ -152,9 +152,8 @@ struct ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} function ZernikeSurface(semidiameter::T; radius::T = typemax(T), conic::T = zero(T), zcoeff::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter, indexing::ZernikeIndexType = ZernikeIndexingOSA) where {T<:Real} asp = AsphericSurface(semidiameter; radius, conic, aspherics, normradius) - #@assert semidiameter > 0 - #@assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) - #@assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" + Q = length(asp.aspherics) #this is not the same as the aspherics variable passed to the function! + zcs = [] if zcoeff !== nothing for (i, k) in zcoeff @@ -170,7 +169,7 @@ struct ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} end P = length(zcs) - new{T,3,P,Q}(asp, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + new{T,3,P, Q}(asp, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end end @@ -188,26 +187,32 @@ function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T< pnt = point(z.asp, ρ, ϕ) # sum zernike - u = r / z.normradius + rad=z.asp.semidiameter + r=ρ * rad + u = r / z.asp.normradius + h = zero(T) @inbounds @simd for m in 1:P (R, S, k) = z.coeffs[m] - pnt[3] += k * Zernike.ζ(R, S, u, ϕ) + h += k * Zernike.ζ(R, S, u, ϕ) end - return pnt #pnt is already the right type and length + return SVector{3,T}(pnt[1], pnt[2], pnt[3] + h) end function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q} pu,pv = partials(z.asp, ρ, ϕ) # sum zernike partials + rad=z.asp.semidiameter n = rad / z.asp.normradius u = ρ * n + dpu = zero(T) + dpv = zero(T) @inbounds @simd for m in 1:P (R, S, k) = z.coeffs[m] du, dϕ = Zernike.δζ(R, S, u, ϕ) - pu[3] += k * du * n # want the derivative wrt ρ, not u - pv[3] += k * dϕ + dpu += k * du * n # want the derivative wrt ρ, not u + dpv += k * dϕ end - return pu, pv + return SVector{3,T}(pu[1], pu[2], pu[3] + dpu), SVector{3,T}(pv[1], pv[2], pv[3] + dpv) end function normal(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q} diff --git a/test/TestData/lenses.jl b/test/TestData/lenses.jl index 31d4d2f30..c746958ed 100644 --- a/test/TestData/lenses.jl +++ b/test/TestData/lenses.jl @@ -9,15 +9,29 @@ function zernikelens() return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end +function evensphericlens() + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, radius = -50.0,aspherics = [(2, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), Geometry.translation(0.0, 0.0, 5.0)) + botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) + barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) + return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) +end + +function oddasphericlens() + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, radius = -50.0,aspherics = [(3, 0.001), (5, -0.0001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), Geometry.translation(0.0, 0.0, 5.0)) + botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) + barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) + return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) +end + function asphericlens() - topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, aspherics = [(2, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, aspherics = [(3, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end -function coniclensZ() - topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) +function coniclensA() + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) diff --git a/test/TestData/systems.jl b/test/TestData/systems.jl index b8d37bfeb..30ea7f32c 100644 --- a/test/TestData/systems.jl +++ b/test/TestData/systems.jl @@ -148,7 +148,9 @@ end zernikesystem() = CSGOpticalSystem(LensAssembly(zernikelens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) asphericsystem() = CSGOpticalSystem(LensAssembly(asphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) -conicsystemZ() = CSGOpticalSystem(LensAssembly(coniclensZ()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +evenapshericsystem()=CSGOpticalSystem(LensAssembly(evensphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +oddapshericsystem()=CSGOpticalSystem(LensAssembly(oddsphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +conicsystemA() = CSGOpticalSystem(LensAssembly(coniclensA()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) conicsystemQ() = CSGOpticalSystem(LensAssembly(coniclensQ()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) qtypesystem() = CSGOpticalSystem(LensAssembly(qtypelens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) chebyshevsystem() = CSGOpticalSystem(LensAssembly(chebyshevlens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) diff --git a/test/testsets/Comparison.jl b/test/testsets/Comparison.jl index 00a858638..f9f589f85 100644 --- a/test/testsets/Comparison.jl +++ b/test/testsets/Comparison.jl @@ -157,7 +157,7 @@ r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ) r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ) - a = TestData.conicsystemZ() + a = TestData.conicsystemA() res = trace(a, r1, test = true) @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) @@ -195,6 +195,50 @@ @test isapprox(point(track[end]), [5.49905367197174, 5.66882664623822, 0.0], rtol = COMP_TOLERANCE) @test isapprox(pathlength(track[end]), 47.6825931025333, rtol = COMP_TOLERANCE) + a = evenasphericsystem() + res = trace(a, r1, test = true) + @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) + res = trace(a, r2, test = true) + @test isapprox(point(res), [-1.1537225076569997, -1.153722507656998, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 74.08191763367685, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r3, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [-4.684043128742365, -4.684043128742363, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 45.610874251552275, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r4, test = true, trackrays = track) + @test isapprox(point(res), [0.0, 15.887441944041871, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 77.14042689028618, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r5, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [5.095756154983988, 5.24166010736873, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 45.618062468023375, rtol = COMP_TOLERANCE) + + a = oddasphericsystem() + res = trace(a, r1, test = true) + @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) + res = trace(a, r2, test = true) + @test isapprox(point(res), [0.6256581890889898, 0.6256581890889894, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 73.97868203031574, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r3, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [-6.097868437446996, -6.097868437446997, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 48.34832912579036, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r4, test = true, trackrays = track) + @test isapprox(point(res), [0.0, 7.485281534548376, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 75.07125754375438, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r5, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [6.197436823172213, 6.440993733499131, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 48.292537161700146, rtol = COMP_TOLERANCE) + a = TestData.zernikesystem() res = trace(a, r1, test = true) @test isapprox(point(res), [0.0, -8.9787010034042, -67.8], rtol = COMP_TOLERANCE) diff --git a/test/testsets/Intersection.jl b/test/testsets/Intersection.jl index 96682cb6c..f4dc12d9d 100644 --- a/test/testsets/Intersection.jl +++ b/test/testsets/Intersection.jl @@ -658,7 +658,6 @@ @testset "Zernike" begin @test_nowarn ZernikeSurface(1.5) @test_throws AssertionError ZernikeSurface(0.0) - @test_throws AssertionError ZernikeSurface(1.5, aspherics = [(1, 0.1)]) z1 = TestData.zernikesurface1() az1 = AcceleratedParametricSurface(z1, 20) From e12c8d8982abf4ec753d12ba5f25e3ec04c8fdc9 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Wed, 13 Oct 2021 16:04:50 -0700 Subject: [PATCH 05/19] fix documentation for AsphericSurface --- src/Geometry/Primitives/AsphericSurface.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 38e2c004d..73f0d1f02 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -22,7 +22,7 @@ Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defi AsphericSurface(semidiameter; radius, conic, aspherics=nothing, normradius = semidiameter) ``` -`aspherics` should be a vector containing the aspheric polynomial coefficients starting at A1 +`aspherics` aspherics should be vectors containing tuples of the form (i, v) where i is the polynomial power of the aspheric term The sag is defined by the equation @@ -30,7 +30,7 @@ The sag is defined by the equation z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir^{i} ``` -where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` . +where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, and ``k = \\texttt{conic}`` . """ struct AsphericSurface{T,N,Q} <: ParametricSurface{T,N} semidiameter::T From a6599ed99c657be2ecaf014807a220ea927310b9 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Fri, 15 Oct 2021 19:40:07 -0700 Subject: [PATCH 06/19] all working, but allocation problems --- src/Geometry/Primitives/AsphericSurface.jl | 295 ++++++++++++++++++--- src/Geometry/Primitives/Zernike.jl | 4 +- src/Optical/Lenses.jl | 2 +- test/TestData/lenses.jl | 2 +- test/TestData/systems.jl | 4 +- test/testsets/Comparison.jl | 4 +- 6 files changed, 260 insertions(+), 51 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 73f0d1f02..13c897ead 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -9,7 +9,10 @@ """ - AsphericSurface{T,N,Q} <: ParametricSurface{T,N} + +```julia +AsphericSurface(semidiameter; radius, conic, aspherics=nothing, normradius = semidiameter) +``` Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter,. `T` is the datatype, `N` is the dimensionality, and `Q` is the number of aspheric terms. @@ -17,11 +20,6 @@ Surface incorporating an aspheric polynomial - radius, conic and aspherics are d The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space. Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned. - -```julia -AsphericSurface(semidiameter; radius, conic, aspherics=nothing, normradius = semidiameter) -``` - `aspherics` aspherics should be vectors containing tuples of the form (i, v) where i is the polynomial power of the aspheric term The sag is defined by the equation @@ -31,39 +29,131 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir ``` where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, and ``k = \\texttt{conic}`` . + +The function checks if the aspheric terms are even, odd or both and uses `EvenAsphericSurface()`, +`OddAsphericSurface()`, or `OddEvenAsphericSurface()` as appropriate. + +""" +abstract type AsphericSurface{T,N} <: ParametricSurface{T,N} end + +function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T= zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} + @assert semidiameter > 0 + @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) + @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" + + acs = [] + if isnothing(aspherics) + surf = EvenAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(acs); normradius) + else + asphericTerms = [i for (i, ) in aspherics] + minAsphericTerm = minimum(asphericTerms) + maxAsphericTerm = maximum(asphericTerms) + @assert minAsphericTerm > 0 "Aspheric Terms must be Order 1 or higher ($minAsphericTerm)" + acs = zeros(T, maxAsphericTerm ) + for (i, k) in aspherics + acs[i] = k + end + odd = any([isodd(i) && a != zero(T) for (i, a) in enumerate(acs)]) + even = any([iseven(i) && a != zero(T) for (i, a) in enumerate(acs)]) + if odd && even + Q = maxAsphericTerm + surf = OddEvenAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(acs); normradius) + elseif even + eacs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)] + #Q = length(eacs) + surf = EvenAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(eacs); normradius) + else #odd + oacs = [acs[2i-1] for i in 1:((maxAsphericTerm+1) ÷ 2)] + Q = length(oacs) + surf = OddAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(oacs); normradius) + end + end + + + return surf #the returned value is an AspericSurface but is also either EvenAphericSurface, OddAphericSurface or OddEvenAsphericSurface +end + +# These methods might be more nicely done with macros to avoid code duplication + """ -struct AsphericSurface{T,N,Q} <: ParametricSurface{T,N} + +```julia +EvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) +``` + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. + +`aspherics` should be an array of the even coefficients of the aspheric polynomial + + +""" +struct EvenAsphericSurface{T,N, Q} <: AsphericSurface{T,N} semidiameter::T - aspherics::SVector{Q,T} curvature::T - conic::T + conic::T + aspherics::SVector{Q,T} #not static so can be optimized normradius::T boundingcylinder::Cylinder{T,N} + function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + new{T,3, Q}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + end +end - function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T= zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} - @assert semidiameter > 0 - @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) - @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" - - acs = [] - if aspherics !== nothing - asphericTerms = [i for (i, ) in aspherics] - minAsphericTerm = minimum(asphericTerms) - maxAsphericTerm = maximum(asphericTerms) - @assert minAsphericTerm > 0 "Aspheric Terms must be A1 or higher (A$minAsphericTerm)" - acs = zeros(T, maxAsphericTerm ) - for (i, k) in aspherics - acs[i] = k - end - end - Q = length(acs) - new{T,3,Q}(semidiameter, SVector{Q,T}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder +""" + +```julia +OddAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) +``` + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. + +`aspherics` should be an array of the odd coefficients of the aspheric polynomial starting with A1 + +""" + +struct OddAsphericSurface{T,N, Q} <: AsphericSurface{T,N} + semidiameter::T + curvature::T + conic::T + aspherics::SVector{Q,T} + normradius::T + boundingcylinder::Cylinder{T,N} + function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + new{T,3,Q}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) end +end +""" + +```julia +OddEvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) +``` + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. + +`aspherics` should be an array of the both odd and even coefficients of the aspheric polynomial starting with A1 + +""" +struct OddEvenAsphericSurface{T,N, Q} <: AsphericSurface{T,N} + semidiameter::T + curvature::T + conic::T + aspherics::SVector{Q,T} + normradius::T + boundingcylinder::Cylinder{T,N} + function OddEvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + new{T,3, Q}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + end end -export AsphericSurface -uvrange(::Type{AsphericSurface{T,N,Q}}) where {T<:Real,N,Q} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ + +export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface + + +uvrange(::AsphericSurface{T,N}) where {T<:Real,N} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ semidiameter(z::AsphericSurface{T}) where {T<:Real} = z.semidiameter halfsizeu(z::AsphericSurface{T}) where {T<:Real} = semidiameter(z) @@ -71,7 +161,51 @@ halfsizev(z::AsphericSurface{T}) where {T<:Real} = semidiameter(z) boundingobj(z::AsphericSurface{T}) where {T<:Real} = z.boundingcylinder -function point(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} +function point(z::OddEvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} + rad = z.semidiameter + r = ρ * rad + r2 = r^2 + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN) + end + h = z.curvature * r2 / (one(T) + sqrt(t)) + # sum aspheric + if Q != 0 + asp,rest = Iterators.peel(z.aspherics) + prod = r + h += asp * prod + for asp in rest + prod *= r + h += asp * prod + end + end + return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) +end + +function point(z::EvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real, Q} + rad = z.semidiameter + r = ρ * rad + r2 = r^2 + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN) + end + h = z.curvature * r2 / (one(T) + sqrt(t)) + # sum aspheric + if Q != 0 + asp,rest = Iterators.peel(z.aspherics) + prod = r2 + h += asp * prod + for asp in rest + prod *= r2 + h += asp * prod + end + end + return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) +end + +function point(z::OddAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} rad = z.semidiameter r = ρ * rad r2 = r^2 @@ -81,15 +215,19 @@ function point(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<: end h = z.curvature * r2 / (one(T) + sqrt(t)) # sum aspheric - prod = one(T) - for asp in z.aspherics - prod *= r + if Q != 0 + asp,rest = Iterators.peel(z.aspherics) + prod = r h += asp * prod + for asp in rest + prod *= r2 #one extra mulitply on last iteration + h += asp * prod + end end return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) end -function partials(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} +function partials(z::OddEvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} rad = z.semidiameter r = ρ * rad t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 @@ -98,7 +236,7 @@ function partials(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},S end dhdρ = rad * z.curvature * r * sqrt(t) / t # sum aspherics partial - if length(z.aspherics) != 0 + if Q != 0 ((m, asp), rest) = Iterators.peel(enumerate(z.aspherics)) dhdρ += rad * asp #first term m=1 and prod = one(T) prod = one(T) @@ -115,7 +253,61 @@ function partials(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},S return pu, pv end -function normal(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} +function partials(z::EvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} + rad = z.semidiameter + r = ρ * rad + r2 = r*r + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) + end + dhdρ = rad * z.curvature * r * sqrt(t) / t + # sum aspherics partial + if Q != 0 + ((m, asp), rest) = Iterators.peel(enumerate(z.aspherics)) + prod = r + dhdρ += rad * asp * 2 * prod #first term m=1*2 and prod = r + for (m,asp) in rest + prod *= r2 + dhdρ += rad * 2m * asp * prod + end + end + dhdϕ = zero(T) + cosϕ = cos(ϕ) + sinϕ = sin(ϕ) + pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) + pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) + return pu, pv +end + +function partials(z::OddAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} + rad = z.semidiameter + r = ρ * rad + r2 = r*r + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) + end + dhdρ = rad * z.curvature * r * sqrt(t) / t + # sum aspherics partial + if Q != 0 + ((m, asp), rest) = Iterators.peel(enumerate(z.aspherics)) + dhdρ += rad * asp #first term m=1 and prod = one(T) + prod = one(T) + for (m,asp) in rest + prod *= r2 + dhdρ += rad * (2m-1) * asp * prod + end + end + dhdϕ = zero(T) + cosϕ = cos(ϕ) + sinϕ = sin(ϕ) + pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) + pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) + return pu, pv +end + +function normal(z::AsphericSurface{T,3}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real} du, dv = partials(z, ρ, ϕ) if ρ == zero(T) && norm(dv) == zero(T) # in cases where there is no δϕ at ρ = 0 (i.e. anything which is rotationally symetric) @@ -125,7 +317,7 @@ function normal(z::AsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T< return normalize(cross(du, dv)) end -function uv(z::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real,Q} +function uv(z::AsphericSurface{T,3}, p::SVector{3,T}) where {T<:Real} # avoid divide by zero for ForwardDiff ϕ = NaNsafeatan(p[2], p[1]) if p[1] == zero(T) && p[2] == zero(T) @@ -136,7 +328,7 @@ function uv(z::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real,Q} return SVector{2,T}(ρ, ϕ) end -function onsurface(surf::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real,Q} +function onsurface(surf::AsphericSurface{T,3}, p::SVector{3,T}) where {T<:Real} ρ, ϕ = uv(surf, p) if ρ > one(T) return false @@ -146,7 +338,7 @@ function onsurface(surf::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real end end -function inside(surf::AsphericSurface{T,3,Q}, p::SVector{3,T}) where {T<:Real,Q} +function inside(surf::AsphericSurface{T,3}, p::SVector{3,T}) where {T<:Real} ρ, ϕ = uv(surf, p) if ρ > one(T) return false @@ -159,7 +351,7 @@ end ######################################################################################################### # Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis. -function surfaceintersection(surf::AcceleratedParametricSurface{T,3,AsphericSurface{T,3,Q}}, r::AbstractRay{T,3}) where {T<:Real,Q} +function surfaceintersection(surf::AcceleratedParametricSurface{T,3,AsphericSurface{T,3}}, r::AbstractRay{T,3}) where {T<:Real} cylint = surfaceintersection(surf.surface.boundingcylinder, r) if cylint isa EmptyInterval{T} return EmptyInterval(T) @@ -198,14 +390,31 @@ function AcceleratedParametricSurface(surf::T, numsamples::Int = 17; interface:: return a end -function BoundingBox(surf::AsphericSurface{T,3,Q}) where {T<:Real,Q} +function asphericSag(surf::EvenAsphericSurface{T,3,Q}) where {T<:Real,Q} + amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) + amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) + return amin, amax +end + +function asphericSag(surf::OddAsphericSurface{T,3,Q}) where {T<:Real,Q} + amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) + amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) + return amin, amax +end + +function asphericSag(surf::OddEvenAsphericSurface{T,3,Q}) where {T<:Real,Q} + amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) + amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) + return amin, amax +end + +function BoundingBox(surf::AsphericSurface{T,3}) where {T<:Real} xmin = -semidiameter(surf) xmax = semidiameter(surf) ymin = -semidiameter(surf) ymax = semidiameter(surf) # aspherics can be more comlpicated, so just take sum of all negative and all positives - amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in surf.aspherics) : zero(T) - amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in surf.aspherics) : zero(T) + amin, amax = asphericSag(surf) zmin = amin zmax = amax # curvature only goes one way diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 658c335ba..d6bb9e9f6 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -141,7 +141,7 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Z_n`` is the nᵗʰ Zernike polynomial. """ struct ZernikeSurface{T,N,P, Q} <: ParametricSurface{T,N} - asp::AsphericSurface{T,N,Q} + asp::AsphericSurface{T,N} #semidiameter::T coeffs::SVector{P,Tuple{Int,Int,T}} #aspherics::SVector{Q,Tuple{Int,T}} @@ -195,7 +195,7 @@ function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T< (R, S, k) = z.coeffs[m] h += k * Zernike.ζ(R, S, u, ϕ) end - return SVector{3,T}(pnt[1], pnt[2], pnt[3] + h) + return SVector{3,T}(pnt[1], pnt[2], pnt[3] + h) end function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q} diff --git a/src/Optical/Lenses.jl b/src/Optical/Lenses.jl index 85396cac4..5ea3cdc42 100644 --- a/src/Optical/Lenses.jl +++ b/src/Optical/Lenses.jl @@ -69,7 +69,7 @@ function AsphericLens(insidematerial::T, frontvertex::S, frontradius::S, frontco if frontaspherics !== nothing frontaspherics = [(i, -k) for (i, k) in frontaspherics] end - surf = AcceleratedParametricSurface(ZernikeSurface(semidiameter + frontdecenter_l + S(0.01), radius = -frontradius, conic = frontconic, aspherics = frontaspherics), nsamples, interface = opticinterface(S, insidematerial, lastmaterial, frontsurfacereflectance, interfacemode)) + surf = AcceleratedParametricSurface(AsphericSurface(semidiameter + frontdecenter_l + S(0.01), radius = -frontradius, conic = frontconic, aspherics = frontaspherics), nsamples, interface = opticinterface(S, insidematerial, lastmaterial, frontsurfacereflectance, interfacemode)) lens_front = leaf(surf, translation(S, frontdecenter[1], frontdecenter[2], frontvertex)) end if isinf(backradius) && (backaspherics === nothing) diff --git a/test/TestData/lenses.jl b/test/TestData/lenses.jl index c746958ed..b9341fdf7 100644 --- a/test/TestData/lenses.jl +++ b/test/TestData/lenses.jl @@ -9,7 +9,7 @@ function zernikelens() return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end -function evensphericlens() +function evenasphericlens() topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, radius = -50.0,aspherics = [(2, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), Geometry.translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) diff --git a/test/TestData/systems.jl b/test/TestData/systems.jl index 30ea7f32c..d5fe9d767 100644 --- a/test/TestData/systems.jl +++ b/test/TestData/systems.jl @@ -148,8 +148,8 @@ end zernikesystem() = CSGOpticalSystem(LensAssembly(zernikelens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) asphericsystem() = CSGOpticalSystem(LensAssembly(asphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) -evenapshericsystem()=CSGOpticalSystem(LensAssembly(evensphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) -oddapshericsystem()=CSGOpticalSystem(LensAssembly(oddsphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +evenasphericsystem()=CSGOpticalSystem(LensAssembly(evenasphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +oddasphericsystem()=CSGOpticalSystem(LensAssembly(oddasphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) conicsystemA() = CSGOpticalSystem(LensAssembly(coniclensA()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) conicsystemQ() = CSGOpticalSystem(LensAssembly(coniclensQ()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) qtypesystem() = CSGOpticalSystem(LensAssembly(qtypelens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) diff --git a/test/testsets/Comparison.jl b/test/testsets/Comparison.jl index f9f589f85..5eb392fd8 100644 --- a/test/testsets/Comparison.jl +++ b/test/testsets/Comparison.jl @@ -195,7 +195,7 @@ @test isapprox(point(track[end]), [5.49905367197174, 5.66882664623822, 0.0], rtol = COMP_TOLERANCE) @test isapprox(pathlength(track[end]), 47.6825931025333, rtol = COMP_TOLERANCE) - a = evenasphericsystem() + a = TestData.evenasphericsystem() res = trace(a, r1, test = true) @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) @@ -217,7 +217,7 @@ @test isapprox(point(track[end]), [5.095756154983988, 5.24166010736873, 0.0], rtol = COMP_TOLERANCE) @test isapprox(pathlength(track[end]), 45.618062468023375, rtol = COMP_TOLERANCE) - a = oddasphericsystem() + a = TestData.oddasphericsystem() res = trace(a, r1, test = true) @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) From 1e37fb4e3478800a7467c0a780eb61cc12b14b41 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Fri, 15 Oct 2021 20:54:26 -0700 Subject: [PATCH 07/19] Simple change to test --- src/Geometry/Primitives/AsphericSurface.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 13c897ead..f1799fffe 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -7,7 +7,6 @@ ######################################################################################################### - """ ```julia From 7b7cb4f4c59a585d3724316a6bac323e21dcdecc Mon Sep 17 00:00:00 2001 From: mattderstine Date: Fri, 22 Oct 2021 12:04:37 -0700 Subject: [PATCH 08/19] Fix variable names of ZernikeSurface --- src/Geometry/Primitives/Zernike.jl | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index d6bb9e9f6..287f4f7b5 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -141,13 +141,8 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Z_n`` is the nᵗʰ Zernike polynomial. """ struct ZernikeSurface{T,N,P, Q} <: ParametricSurface{T,N} - asp::AsphericSurface{T,N} - #semidiameter::T + asp::AsphericSurface{T,N,Q, R} coeffs::SVector{P,Tuple{Int,Int,T}} - #aspherics::SVector{Q,Tuple{Int,T}} - #curvature::T - #conic::T - #normradius::T boundingcylinder::Cylinder{T,N} function ZernikeSurface(semidiameter::T; radius::T = typemax(T), conic::T = zero(T), zcoeff::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter, indexing::ZernikeIndexType = ZernikeIndexingOSA) where {T<:Real} @@ -199,20 +194,20 @@ function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T< end function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q} - pu,pv = partials(z.asp, ρ, ϕ) + pρ,pϕ = partials(z.asp, ρ, ϕ) # sum zernike partials rad=z.asp.semidiameter n = rad / z.asp.normradius u = ρ * n - dpu = zero(T) - dpv = zero(T) + dhdρ = zero(T) + dhdϕ = zero(T) @inbounds @simd for m in 1:P (R, S, k) = z.coeffs[m] du, dϕ = Zernike.δζ(R, S, u, ϕ) - dpu += k * du * n # want the derivative wrt ρ, not u - dpv += k * dϕ + dhdρ += k * du * n # want the derivative wrt ρ, not u + dhdϕ += k * dϕ end - return SVector{3,T}(pu[1], pu[2], pu[3] + dpu), SVector{3,T}(pv[1], pv[2], pv[3] + dpv) + return SVector{3,T}(pρ[1], pρ[2], pρ[3] + dhdρ), SVector{3,T}(pϕ[1], pϕ[2], pϕ[3] + dhdϕ) end function normal(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q} From 14f2c87289cff7d37fe5b0effdc5b517b6c442e4 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Fri, 29 Oct 2021 17:05:53 -0700 Subject: [PATCH 09/19] evaluate based on on aspheric pseudo-types --- src/Geometry/Primitives/AsphericSurface.jl | 286 +++++++-------------- src/Geometry/Primitives/Zernike.jl | 12 +- test/testsets/JuliaLang.jl | 1 + 3 files changed, 105 insertions(+), 194 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index f1799fffe..a59923564 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -33,46 +33,77 @@ The function checks if the aspheric terms are even, odd or both and uses `EvenAs `OddAsphericSurface()`, or `OddEvenAsphericSurface()` as appropriate. """ -abstract type AsphericSurface{T,N} <: ParametricSurface{T,N} end +#"pseudo-types" of aspheres +const CONIC = 0 +const ODD = 1 +const EVEN = 2 +const ODDEVEN = 3 -function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T= zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} - @assert semidiameter > 0 - @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) - @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" +struct AsphericSurface{T,N, Q, M} <: ParametricSurface{T,N} + semidiameter::T + curvature::T + conic::T + aspherics::SVector{Q,T} + normradius::T + boundingcylinder::Cylinder{T,N} - acs = [] - if isnothing(aspherics) - surf = EvenAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(acs); normradius) - else - asphericTerms = [i for (i, ) in aspherics] - minAsphericTerm = minimum(asphericTerms) - maxAsphericTerm = maximum(asphericTerms) - @assert minAsphericTerm > 0 "Aspheric Terms must be Order 1 or higher ($minAsphericTerm)" - acs = zeros(T, maxAsphericTerm ) - for (i, k) in aspherics - acs[i] = k - end - odd = any([isodd(i) && a != zero(T) for (i, a) in enumerate(acs)]) - even = any([iseven(i) && a != zero(T) for (i, a) in enumerate(acs)]) - if odd && even - Q = maxAsphericTerm - surf = OddEvenAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(acs); normradius) - elseif even - eacs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)] - #Q = length(eacs) - surf = EvenAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(eacs); normradius) - else #odd - oacs = [acs[2i-1] for i in 1:((maxAsphericTerm+1) ÷ 2)] - Q = length(oacs) - surf = OddAsphericSurface(semidiameter, 1 / radius, conic, Vector{T}(oacs); normradius) + function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T= zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} + @assert semidiameter > 0 + @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) + @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" + + acs = [] + if aspherics===nothing + surf = new{T,3, 0, CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder{T,3}(semidiameter, interface = opaqueinterface(T))) + else + asphericTerms = [i for (i, ) in aspherics] + minAsphericTerm = minimum(asphericTerms) + maxAsphericTerm = maximum(asphericTerms) + @assert minAsphericTerm > 0 "Aspheric Terms must be Order 1 or higher ($minAsphericTerm)" + acs = zeros(T, maxAsphericTerm ) + for (i, k) in aspherics + acs[i] = k + end + odd = any([isodd(i) && a != zero(T) for (i, a) in enumerate(acs)]) + even = any([iseven(i) && a != zero(T) for (i, a) in enumerate(acs)]) + if odd && even + Q = maxAsphericTerm + surf = new{T, 3, Q, ODDEVEN}(semidiameter, 1 / radius, conic, Vector{T}(acs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + elseif even + eacs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)] + Q = length(eacs) + surf = new{T, 3, Q, EVEN}(semidiameter, 1 / radius, conic, Vector{T}(eacs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + else #odd + oacs = [acs[2i-1] for i in 1:((maxAsphericTerm+1) ÷ 2)] + Q = length(oacs) + surf = new{T, 3, Q, ODD}(semidiameter, 1 / radius, conic, Vector{T}(oacs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + end end + return surf + end + function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + new{T,3, Q, EVEN}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + end + function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + new{T,3,Q, ODD}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + end + function OddEvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + new{T,3,Q,CONIC}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) end - - return surf #the returned value is an AspericSurface but is also either EvenAphericSurface, OddAphericSurface or OddEvenAsphericSurface end -# These methods might be more nicely done with macros to avoid code duplication + +asphericType(z::AsphericSurface{T,3,Q,EVEN}) where {T<:Real,Q} = EVEN + +asphericType(z::AsphericSurface{T,3,Q,CONIC}) where {T<:Real,Q} = CONIC + +asphericType(z::AsphericSurface{T,3,Q,ODD}) where {T<:Real,Q} = ODD + +asphericType(z::AsphericSurface{T,3,Q,ODDEVEN}) where {T<:Real,Q} = ODDEVEN """ @@ -86,18 +117,6 @@ Surface incorporating an aspheric polynomial - radius, conic and aspherics are d """ -struct EvenAsphericSurface{T,N, Q} <: AsphericSurface{T,N} - semidiameter::T - curvature::T - conic::T - aspherics::SVector{Q,T} #not static so can be optimized - normradius::T - boundingcylinder::Cylinder{T,N} - function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real - Q=length(aspherics) - new{T,3, Q}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - end -end """ @@ -111,19 +130,6 @@ Surface incorporating an aspheric polynomial - radius, conic and aspherics are d """ -struct OddAsphericSurface{T,N, Q} <: AsphericSurface{T,N} - semidiameter::T - curvature::T - conic::T - aspherics::SVector{Q,T} - normradius::T - boundingcylinder::Cylinder{T,N} - function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real - Q=length(aspherics) - new{T,3,Q}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - end -end - """ ```julia @@ -135,18 +141,6 @@ Surface incorporating an aspheric polynomial - radius, conic and aspherics are d `aspherics` should be an array of the both odd and even coefficients of the aspheric polynomial starting with A1 """ -struct OddEvenAsphericSurface{T,N, Q} <: AsphericSurface{T,N} - semidiameter::T - curvature::T - conic::T - aspherics::SVector{Q,T} - normradius::T - boundingcylinder::Cylinder{T,N} - function OddEvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real - Q=length(aspherics) - new{T,3, Q}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - end -end export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface @@ -160,51 +154,11 @@ halfsizev(z::AsphericSurface{T}) where {T<:Real} = semidiameter(z) boundingobj(z::AsphericSurface{T}) where {T<:Real} = z.boundingcylinder -function point(z::OddEvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} - rad = z.semidiameter - r = ρ * rad - r2 = r^2 - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN) - end - h = z.curvature * r2 / (one(T) + sqrt(t)) - # sum aspheric - if Q != 0 - asp,rest = Iterators.peel(z.aspherics) - prod = r - h += asp * prod - for asp in rest - prod *= r - h += asp * prod - end - end - return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) -end - -function point(z::EvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real, Q} - rad = z.semidiameter - r = ρ * rad - r2 = r^2 - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN) - end - h = z.curvature * r2 / (one(T) + sqrt(t)) - # sum aspheric - if Q != 0 - asp,rest = Iterators.peel(z.aspherics) - prod = r2 - h += asp * prod - for asp in rest - prod *= r2 - h += asp * prod - end - end - return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) -end +prod_step(z::AsphericSurface{T,N,Q,ODDEVEN}, r, r2) where {T<:Real,N,Q} = r, r +prod_step(z::AsphericSurface{T,N,Q,ODD}, r, r2) where {T<:Real,N,Q} = r, r2 +prod_step(z::AsphericSurface{T,N,Q,EVEN}, r, r2) where {T<:Real,N,Q} = r2, r2 -function point(z::OddAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q} +function point(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q,M} rad = z.semidiameter r = ρ * rad r2 = r^2 @@ -213,46 +167,24 @@ function point(z::OddAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::SVector{3,T} where { return SVector{3,T}(NaN, NaN, NaN) end h = z.curvature * r2 / (one(T) + sqrt(t)) - # sum aspheric - if Q != 0 + # sum aspheric + if M != CONIC + prod, step = prod_step(z, r, r2) #multiple dispatch on R asp,rest = Iterators.peel(z.aspherics) - prod = r h += asp * prod for asp in rest - prod *= r2 #one extra mulitply on last iteration + prod *= step h += asp * prod end end return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) end -function partials(z::OddEvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} - rad = z.semidiameter - r = ρ * rad - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) - end - dhdρ = rad * z.curvature * r * sqrt(t) / t - # sum aspherics partial - if Q != 0 - ((m, asp), rest) = Iterators.peel(enumerate(z.aspherics)) - dhdρ += rad * asp #first term m=1 and prod = one(T) - prod = one(T) - for (m,asp) in rest - prod *= r - dhdρ += rad * m * asp * prod - end - end - dhdϕ = zero(T) - cosϕ = cos(ϕ) - sinϕ = sin(ϕ) - pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) - pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) - return pu, pv -end +partial_prod_step(z::AsphericSurface{T, 3, Q, EVEN}) where {T<:Real,Q} = r, r2, 2:2:2Q +partial_prod_step(z::AsphericSurface{T, 3, Q, ODD}) where {T<:Real,Q} = one(T), r2, 1:2:(2Q-1) +partial_prod_step(z::AsphericSurface{T, 3, Q, ODDEVEN}) where {T<:Real,Q} = one(T), r, 1:1:Q -function partials(z::EvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} +function partials(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q,M} rad = z.semidiameter r = ρ * rad r2 = r*r @@ -262,48 +194,20 @@ function partials(z::EvenAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3, end dhdρ = rad * z.curvature * r * sqrt(t) / t # sum aspherics partial - if Q != 0 - ((m, asp), rest) = Iterators.peel(enumerate(z.aspherics)) - prod = r + if M != CONIC + prod, step, mIter = partial_prod_step(z, r, r2) + + ((m, asp), rest) = Iterators.peel(zip(mIter,z.aspherics)) dhdρ += rad * asp * 2 * prod #first term m=1*2 and prod = r for (m,asp) in rest - prod *= r2 - dhdρ += rad * 2m * asp * prod - end - end - dhdϕ = zero(T) - cosϕ = cos(ϕ) - sinϕ = sin(ϕ) - pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) - pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) - return pu, pv -end - -function partials(z::OddAsphericSurface{T,3,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q} - rad = z.semidiameter - r = ρ * rad - r2 = r*r - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) - end - dhdρ = rad * z.curvature * r * sqrt(t) / t - # sum aspherics partial - if Q != 0 - ((m, asp), rest) = Iterators.peel(enumerate(z.aspherics)) - dhdρ += rad * asp #first term m=1 and prod = one(T) - prod = one(T) - for (m,asp) in rest - prod *= r2 - dhdρ += rad * (2m-1) * asp * prod + prod *= step + dhdρ += rad * m * asp * prod end end dhdϕ = zero(T) cosϕ = cos(ϕ) sinϕ = sin(ϕ) - pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) - pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) - return pu, pv + return SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ), SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) end function normal(z::AsphericSurface{T,3}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real} @@ -389,21 +293,27 @@ function AcceleratedParametricSurface(surf::T, numsamples::Int = 17; interface:: return a end -function asphericSag(surf::EvenAsphericSurface{T,3,Q}) where {T<:Real,Q} - amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) - amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) +function asphericSag(surf::AsphericSurface{T,3,Q, EVEN}) where {T<:Real,Q} + amin = sum(k < zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) + amax = sum(k > zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) return amin, amax end -function asphericSag(surf::OddAsphericSurface{T,3,Q}) where {T<:Real,Q} - amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) - amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) +function asphericSag(surf::AsphericSurface{T,3,Q, ODD}) where {T<:Real,Q} + amin = sum(k < zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) + amax = sum(k > zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) return amin, amax end -function asphericSag(surf::OddEvenAsphericSurface{T,3,Q}) where {T<:Real,Q} - amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) - amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) : zero(T) +function asphericSag(surf::AsphericSurface{T,3,Q, ODDEVEN}) where {T<:Real,Q} + amin = sum(k < zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) + amax = sum(k > zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) + return amin, amax +end + +function asphericSag(surf::AsphericSurface{T,3,Q, CONIC}) where {T<:Real,Q} + amin = zero(T) + amax = zero(T) return amin, amax end diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 287f4f7b5..4b035a892 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -140,8 +140,8 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Z_n`` is the nᵗʰ Zernike polynomial. """ -struct ZernikeSurface{T,N,P, Q} <: ParametricSurface{T,N} - asp::AsphericSurface{T,N,Q, R} +struct ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N} + asp::AsphericSurface{T,N,Q,M} coeffs::SVector{P,Tuple{Int,Int,T}} boundingcylinder::Cylinder{T,N} @@ -163,8 +163,8 @@ struct ZernikeSurface{T,N,P, Q} <: ParametricSurface{T,N} end end P = length(zcs) - - new{T,3,P, Q}(asp, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + M=asphericType(asp) + new{T,3,P,Q,M}(asp, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end end @@ -182,7 +182,7 @@ function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T< pnt = point(z.asp, ρ, ϕ) # sum zernike - rad=z.asp.semidiameter + rad=semidiameter(z.asp) r=ρ * rad u = r / z.asp.normradius h = zero(T) @@ -197,7 +197,7 @@ function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T}, pρ,pϕ = partials(z.asp, ρ, ϕ) # sum zernike partials rad=z.asp.semidiameter - n = rad / z.asp.normradius + n = rad / semidiameter(z.asp) u = ρ * n dhdρ = zero(T) dhdϕ = zero(T) diff --git a/test/testsets/JuliaLang.jl b/test/testsets/JuliaLang.jl index 75f8b4c79..268b75916 100644 --- a/test/testsets/JuliaLang.jl +++ b/test/testsets/JuliaLang.jl @@ -24,6 +24,7 @@ and pop it from the set (OpticSim.MultiHologramInterface, Tuple{Vararg{HologramInterface{T},N}} where {N} where {T<:Real}), ] : [], OpticSim.Zernike => [], + #OpticSim.AsphericSurface => [], OpticSim.QType => [], OpticSim.Vis => [ (OpticSim.Vis.drawcurves, Tuple{Vararg{Spline{P,S,N,M},N1} where N1} where {M} where {N} where {S} where {P}), From 456d57bff1e4b926837d5f125f6d247e25d69d38 Mon Sep 17 00:00:00 2001 From: mattderstine <59268355+mattderstine@users.noreply.github.com> Date: Sat, 30 Oct 2021 11:18:07 -0700 Subject: [PATCH 10/19] Update src/Geometry/Primitives/AsphericSurface.jl Co-authored-by: Charlie Hewitt --- src/Geometry/Primitives/AsphericSurface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index a59923564..19881780f 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -180,7 +180,7 @@ function point(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) end -partial_prod_step(z::AsphericSurface{T, 3, Q, EVEN}) where {T<:Real,Q} = r, r2, 2:2:2Q +partial_prod_step(z::AsphericSurface{T,3,Q,EVEN}, r: T, r2: T) where {T<:Real,Q} = r, r2, 2:2:2Q partial_prod_step(z::AsphericSurface{T, 3, Q, ODD}) where {T<:Real,Q} = one(T), r2, 1:2:(2Q-1) partial_prod_step(z::AsphericSurface{T, 3, Q, ODDEVEN}) where {T<:Real,Q} = one(T), r, 1:1:Q From 4850084347078a0043b4b60381bdba7fb9f1c67f Mon Sep 17 00:00:00 2001 From: mattderstine <59268355+mattderstine@users.noreply.github.com> Date: Sat, 30 Oct 2021 11:18:14 -0700 Subject: [PATCH 11/19] Update src/Geometry/Primitives/Zernike.jl Co-authored-by: Charlie Hewitt --- src/Geometry/Primitives/Zernike.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 4b035a892..6dac9bda7 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -163,7 +163,7 @@ struct ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N} end end P = length(zcs) - M=asphericType(asp) + M = asphericType(asp) new{T,3,P,Q,M}(asp, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end From 9d1513a3f86145f17ea6e47986de6cad59b754f9 Mon Sep 17 00:00:00 2001 From: mattderstine <59268355+mattderstine@users.noreply.github.com> Date: Sat, 30 Oct 2021 11:18:21 -0700 Subject: [PATCH 12/19] Update src/Geometry/Primitives/Zernike.jl Co-authored-by: Charlie Hewitt --- src/Geometry/Primitives/Zernike.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 6dac9bda7..2a2618860 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -182,7 +182,7 @@ function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T< pnt = point(z.asp, ρ, ϕ) # sum zernike - rad=semidiameter(z.asp) + rad = semidiameter(z.asp) r=ρ * rad u = r / z.asp.normradius h = zero(T) From 8e620330417d03c78db842301fa99dbf8350c21d Mon Sep 17 00:00:00 2001 From: mattderstine <59268355+mattderstine@users.noreply.github.com> Date: Sat, 30 Oct 2021 11:19:15 -0700 Subject: [PATCH 13/19] Update src/Geometry/Primitives/AsphericSurface.jl Co-authored-by: Charlie Hewitt --- src/Geometry/Primitives/AsphericSurface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 19881780f..849a7ed5c 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -83,7 +83,7 @@ struct AsphericSurface{T,N, Q, M} <: ParametricSurface{T,N} end function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real Q=length(aspherics) - new{T,3, Q, EVEN}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + new{T,3,Q,EVEN}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) end function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real Q=length(aspherics) From 438536057de1c56eb1ca22a890a174db6bd4abcc Mon Sep 17 00:00:00 2001 From: mattderstine <59268355+mattderstine@users.noreply.github.com> Date: Sat, 30 Oct 2021 11:19:36 -0700 Subject: [PATCH 14/19] Update src/Geometry/Primitives/AsphericSurface.jl Co-authored-by: Charlie Hewitt --- src/Geometry/Primitives/AsphericSurface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 849a7ed5c..dbc93d7ca 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -54,7 +54,7 @@ struct AsphericSurface{T,N, Q, M} <: ParametricSurface{T,N} acs = [] if aspherics===nothing - surf = new{T,3, 0, CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder{T,3}(semidiameter, interface = opaqueinterface(T))) + surf = new{T,3,0,CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder{T,3}(semidiameter, interface = opaqueinterface(T))) else asphericTerms = [i for (i, ) in aspherics] minAsphericTerm = minimum(asphericTerms) From 3700966d90197e4c45396d707bb93c3d06de4944 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Sat, 30 Oct 2021 11:21:24 -0700 Subject: [PATCH 15/19] fix argument lists, other minor mistakes --- src/Geometry/Primitives/AsphericSurface.jl | 14 ++++++++------ src/Geometry/Primitives/Zernike.jl | 2 +- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index a59923564..4316d01b5 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -54,7 +54,7 @@ struct AsphericSurface{T,N, Q, M} <: ParametricSurface{T,N} acs = [] if aspherics===nothing - surf = new{T,3, 0, CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder{T,3}(semidiameter, interface = opaqueinterface(T))) + surf = new{T,3, 0, CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) else asphericTerms = [i for (i, ) in aspherics] minAsphericTerm = minimum(asphericTerms) @@ -73,10 +73,12 @@ struct AsphericSurface{T,N, Q, M} <: ParametricSurface{T,N} eacs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)] Q = length(eacs) surf = new{T, 3, Q, EVEN}(semidiameter, 1 / radius, conic, Vector{T}(eacs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - else #odd + elseif odd oacs = [acs[2i-1] for i in 1:((maxAsphericTerm+1) ÷ 2)] Q = length(oacs) surf = new{T, 3, Q, ODD}(semidiameter, 1 / radius, conic, Vector{T}(oacs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + else #CONIC + surf = new{T,3, 0, CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) end end return surf @@ -143,7 +145,7 @@ Surface incorporating an aspheric polynomial - radius, conic and aspherics are d """ -export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface +export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface, asphericType uvrange(::AsphericSurface{T,N}) where {T<:Real,N} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ @@ -180,9 +182,9 @@ function point(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) end -partial_prod_step(z::AsphericSurface{T, 3, Q, EVEN}) where {T<:Real,Q} = r, r2, 2:2:2Q -partial_prod_step(z::AsphericSurface{T, 3, Q, ODD}) where {T<:Real,Q} = one(T), r2, 1:2:(2Q-1) -partial_prod_step(z::AsphericSurface{T, 3, Q, ODDEVEN}) where {T<:Real,Q} = one(T), r, 1:1:Q +partial_prod_step(z::AsphericSurface{T, 3, Q, EVEN}, r, r2) where {T<:Real,Q} = r, r2, 2:2:2Q +partial_prod_step(z::AsphericSurface{T, 3, Q, ODD}, r, r2) where {T<:Real,Q} = one(T), r2, 1:2:(2Q-1) +partial_prod_step(z::AsphericSurface{T, 3, Q, ODDEVEN}, r, r2) where {T<:Real,Q} = one(T), r, 1:1:Q function partials(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q,M} rad = z.semidiameter diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 4b035a892..c22b98fb0 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -170,7 +170,7 @@ struct ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N} end export ZernikeSurface -uvrange(::Type{ZernikeSurface{T,N,P,Q}}) where {T<:Real,N,P,Q} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ +uvrange(::Type{ZernikeSurface{T,N}}) where {T<:Real,N} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ semidiameter(z::ZernikeSurface{T}) where {T<:Real} = z.asp.semidiameter halfsizeu(z::ZernikeSurface{T}) where {T<:Real} = semidiameter(z) From ce1321ced04914aece6ecc3a60c3f1080c330423 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Tue, 9 Nov 2021 11:32:20 -0800 Subject: [PATCH 16/19] 1st version to pass all tests. Still need benchmark tests and documentation --- src/Geometry/Primitives/AsphericSurface.jl | 15 +++++-------- src/Geometry/Primitives/Zernike.jl | 26 +++++++++++----------- test/TestData/surfaces.jl | 3 +++ 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 4f2bf5334..02af5f54a 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -34,12 +34,9 @@ The function checks if the aspheric terms are even, odd or both and uses `EvenAs """ #"pseudo-types" of aspheres -const CONIC = 0 -const ODD = 1 -const EVEN = 2 -const ODDEVEN = 3 +@enum AsphSurfaceType CONIC ODD EVEN ODDEVEN -struct AsphericSurface{T,N, Q, M} <: ParametricSurface{T,N} +struct AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} semidiameter::T curvature::T conic::T @@ -171,7 +168,7 @@ function point(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T h = z.curvature * r2 / (one(T) + sqrt(t)) # sum aspheric if M != CONIC - prod, step = prod_step(z, r, r2) #multiple dispatch on R + prod, step = prod_step(z, r, r2) #multiple dispatch on M asp,rest = Iterators.peel(z.aspherics) h += asp * prod for asp in rest @@ -182,9 +179,9 @@ function point(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) end -partial_prod_step(z::AsphericSurface{T,3,Q,EVEN}, r: T, r2: T) where {T<:Real,Q} = r, r2, 2:2:2Q -partial_prod_step(z::AsphericSurface{T,3,Q,ODD}) where {T<:Real,Q} = one(T), r2, 1:2:(2Q-1) -partial_prod_step(z::AsphericSurface{T,3,Q,ODDEVEN}) where {T<:Real,Q} = one(T), r, 1:1:Q +partial_prod_step(z::AsphericSurface{T,3,Q,EVEN}, r::T, r2::T) where {T<:Real,Q} = r, r2, 2:2:2Q +partial_prod_step(z::AsphericSurface{T,3,Q,ODD}, r::T, r2::T) where {T<:Real,Q} = one(T), r2, 1:2:(2Q-1) +partial_prod_step(z::AsphericSurface{T,3,Q,ODDEVEN}, r::T, r2::T) where {T<:Real,Q} = one(T), r, 1:1:Q function partials(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q,M} rad = z.semidiameter diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 7523c40cf..39706bf9d 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -164,13 +164,13 @@ struct ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N} end P = length(zcs) M = asphericType(asp) - new{T,3,P,Q,M}(asp, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + new{T,3,P,Q,M}(asp::AsphericSurface{T,3,Q,M}, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end end export ZernikeSurface -uvrange(::Type{ZernikeSurface{T,N}}) where {T<:Real,N} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ +uvrange(::Type{ZernikeSurface{T,N,P,Q,M}}) where {T<:Real,N,P,Q,M} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ semidiameter(z::ZernikeSurface{T}) where {T<:Real} = z.asp.semidiameter halfsizeu(z::ZernikeSurface{T}) where {T<:Real} = semidiameter(z) @@ -178,12 +178,12 @@ halfsizev(z::ZernikeSurface{T}) where {T<:Real} = semidiameter(z) boundingobj(z::ZernikeSurface{T}) where {T<:Real} = z.boundingcylinder -function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q} +function point(z::ZernikeSurface{T,3,P,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q,M} pnt = point(z.asp, ρ, ϕ) # sum zernike rad = semidiameter(z.asp) - r=ρ * rad + r = ρ * rad u = r / z.asp.normradius h = zero(T) @inbounds @simd for m in 1:P @@ -193,11 +193,11 @@ function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T< return SVector{3,T}(pnt[1], pnt[2], pnt[3] + h) end -function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q} +function partials(z::ZernikeSurface{T,3,P,Q,M}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q,M} pρ,pϕ = partials(z.asp, ρ, ϕ) # sum zernike partials rad=z.asp.semidiameter - n = rad / semidiameter(z.asp) + n = rad / z.asp.normradius u = ρ * n dhdρ = zero(T) dhdϕ = zero(T) @@ -210,7 +210,7 @@ function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T}, return SVector{3,T}(pρ[1], pρ[2], pρ[3] + dhdρ), SVector{3,T}(pϕ[1], pϕ[2], pϕ[3] + dhdϕ) end -function normal(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q} +function normal(z::ZernikeSurface{T,3,P,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q,M} du, dv = partials(z, ρ, ϕ) if ρ == zero(T) && norm(dv) == zero(T) # in cases where there is no δϕ at ρ = 0 (i.e. anything which is rotationally symetric) @@ -220,7 +220,7 @@ function normal(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T return normalize(cross(du, dv)) end -function uv(z::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} +function uv(z::ZernikeSurface{T,3,P,Q,M}, p::SVector{3,T}) where {T<:Real,P,Q,M} # avoid divide by zero for ForwardDiff ϕ = NaNsafeatan(p[2], p[1]) if p[1] == zero(T) && p[2] == zero(T) @@ -231,7 +231,7 @@ function uv(z::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} return SVector{2,T}(ρ, ϕ) end -function onsurface(surf::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} +function onsurface(surf::ZernikeSurface{T,3,P,Q,M}, p::SVector{3,T}) where {T<:Real,P,Q,M} ρ, ϕ = uv(surf, p) if ρ > one(T) return false @@ -241,7 +241,7 @@ function onsurface(surf::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Rea end end -function inside(surf::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} +function inside(surf::ZernikeSurface{T,3,P,Q,M}, p::SVector{3,T}) where {T<:Real,P,Q,M} ρ, ϕ = uv(surf, p) if ρ > one(T) return false @@ -254,7 +254,7 @@ end ######################################################################################################### # Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis. -function surfaceintersection(surf::AcceleratedParametricSurface{T,3,ZernikeSurface{T,3,P,Q}}, r::AbstractRay{T,3}) where {T<:Real,P,Q} +function surfaceintersection(surf::AcceleratedParametricSurface{T,3,ZernikeSurface{T,3,P,Q,M}}, r::AbstractRay{T,3}) where {T<:Real,P,Q,M} cylint = surfaceintersection(surf.surface.boundingcylinder, r) if cylint isa EmptyInterval{T} return EmptyInterval(T) @@ -293,11 +293,11 @@ function AcceleratedParametricSurface(surf::T, numsamples::Int = 17; interface:: return a end -function BoundingBox(surf::ZernikeSurface{T,3,P,Q}) where {T<:Real,P,Q} +function BoundingBox(surf::ZernikeSurface{T,3,P,Q,M}) where {T<:Real,P,Q,M} bb = BoundingBox(z.asp) # zernike terms have condition than |Zᵢ| <= 1 # so this gives us a (loose) bounding box - ak = P > 0 ? sum(abs.(Zernike.normalisation(T, N, M) * k for (N, M, k) in surf.coeffs)) : zero(T) + ak = P > 0 ? sum(abs.(Zernike.normalisation(T, n, m) * k for (n, m, k) in surf.coeffs)) : zero(T) bb.zmin -= ak bb.zmax += ak return BoundingBox(bb.xmin, bb.xmax, bb.ymin, bb.ymax, bb.zmin, bb.zmax) #could just return bb, but this way is safer diff --git a/test/TestData/surfaces.jl b/test/TestData/surfaces.jl index 4b24520a8..de6f8b2d7 100644 --- a/test/TestData/surfaces.jl +++ b/test/TestData/surfaces.jl @@ -84,6 +84,9 @@ zernikesurface2() = ZernikeSurface(1.5, radius = 5.0, conic = 1.0, zcoeff = [(7, zernikesurface3() = ZernikeSurface(9.5, zcoeff = [(1, 0.0), (4, -1.0), (7, 0.5)]) conicsurface() = AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0)) +# +# TODO need to add surfaces to test all types of Aspheric Surfaces +#conicsurface() = AcceleratedParametricSurface(AsphericSurface(9.5, radius = -15.0, conic = -5.0)) function simplesaggrid() points = map( From abd7776f339626d46283602fff8382857285f510 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Wed, 10 Nov 2021 10:23:24 -0800 Subject: [PATCH 17/19] Ready for final check --- docs/src/primitives.md | 1 + docs/src/ref.md | 6 ++ src/Geometry/Primitives/AsphericSurface.jl | 84 +++++++++++----------- src/Geometry/Primitives/Zernike.jl | 9 ++- test/Benchmarks/Benchmarks.jl | 5 +- test/TestData/surfaces.jl | 9 +-- 6 files changed, 66 insertions(+), 48 deletions(-) diff --git a/docs/src/primitives.md b/docs/src/primitives.md index 6b3058608..12546a507 100644 --- a/docs/src/primitives.md +++ b/docs/src/primitives.md @@ -60,6 +60,7 @@ Cylinder Plane Sphere SphericalCap +AsphericSurface ZernikeSurface BezierSurface BSplineSurface diff --git a/docs/src/ref.md b/docs/src/ref.md index 028d09221..9d947cccf 100644 --- a/docs/src/ref.md +++ b/docs/src/ref.md @@ -19,6 +19,12 @@ Modules = [OpticSim] Modules = [OpticSim.Geometry] ``` +## Aspheric + +```@autodocs +Modules = [OpticSim.AsphericSurface] +``` + ## Zernike ```@autodocs diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 02af5f54a..05c48132e 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -9,13 +9,16 @@ """ +AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter,. +`T` is the datatype, `N` is the dimensionality, `Q` is the number of aspheric terms, and `M` is the type of aspheric polynomial. + + ```julia AsphericSurface(semidiameter; radius, conic, aspherics=nothing, normradius = semidiameter) ``` -Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter,. -`T` is the datatype, `N` is the dimensionality, and `Q` is the number of aspheric terms. - The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space. Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned. @@ -29,8 +32,7 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, and ``k = \\texttt{conic}`` . -The function checks if the aspheric terms are even, odd or both and uses `EvenAsphericSurface()`, -`OddAsphericSurface()`, or `OddEvenAsphericSurface()` as appropriate. +The function checks if the aspheric terms are missing, even, odd or both and uses the appropriate polynomial evaluation strategy. """ #"pseudo-types" of aspheres @@ -44,14 +46,18 @@ struct AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} normradius::T boundingcylinder::Cylinder{T,N} - function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T= zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} + function AsphericSurface(M::AsphSurfaceType, semidiameter::T, curvature::T, conic::T, aspherics::SVector{Q,T}, normradius::T, boundingcylinder) where {T<:Real,Q} + new{T,3,Q,M}(semidiameter, curvature, conic, aspherics, normradius, boundingcylinder) + end + + function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T = zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} @assert semidiameter > 0 @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" acs = [] if aspherics===nothing - surf = new{T,3,0,CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + M = CONIC else asphericTerms = [i for (i, ) in aspherics] minAsphericTerm = minimum(asphericTerms) @@ -64,46 +70,25 @@ struct AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} odd = any([isodd(i) && a != zero(T) for (i, a) in enumerate(acs)]) even = any([iseven(i) && a != zero(T) for (i, a) in enumerate(acs)]) if odd && even - Q = maxAsphericTerm - surf = new{T, 3, Q, ODDEVEN}(semidiameter, 1 / radius, conic, Vector{T}(acs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + M = ODDEVEN elseif even - eacs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)] - Q = length(eacs) - surf = new{T, 3, Q, EVEN}(semidiameter, 1 / radius, conic, Vector{T}(eacs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + M = EVEN + acs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)] elseif odd - oacs = [acs[2i-1] for i in 1:((maxAsphericTerm+1) ÷ 2)] - Q = length(oacs) - surf = new{T, 3, Q, ODD}(semidiameter, 1 / radius, conic, Vector{T}(oacs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - else #CONIC - surf = new{T,3, 0, CONIC}(semidiameter, 1/radius, conic, Vector{T}(acs), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + M = ODD + acs = [acs[2i-1] for i in 1:((maxAsphericTerm+1) ÷ 2)] + else #there are no nonzero aspherics terms in the list + M = CONIC + acs = [] end end + Q = length(acs) + surf = new{T,3, Q, M}(semidiameter, 1/radius, conic, acs, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) return surf end - function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real - Q=length(aspherics) - new{T,3,Q,EVEN}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - end - function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real - Q=length(aspherics) - new{T,3,Q, ODD}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - end - function OddEvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real - Q=length(aspherics) - new{T,3,Q,CONIC}(semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) - end end - -asphericType(z::AsphericSurface{T,3,Q,EVEN}) where {T<:Real,Q} = EVEN - -asphericType(z::AsphericSurface{T,3,Q,CONIC}) where {T<:Real,Q} = CONIC - -asphericType(z::AsphericSurface{T,3,Q,ODD}) where {T<:Real,Q} = ODD - -asphericType(z::AsphericSurface{T,3,Q,ODDEVEN}) where {T<:Real,Q} = ODDEVEN - """ ```julia @@ -112,10 +97,14 @@ EvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. -`aspherics` should be an array of the even coefficients of the aspheric polynomial +`aspherics` should be an array of the even coefficients of the aspheric polynomial starting with A2 """ +function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + AsphericSurface(EVEN, semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) +end """ @@ -128,6 +117,10 @@ Surface incorporating an aspheric polynomial - radius, conic and aspherics are d `aspherics` should be an array of the odd coefficients of the aspheric polynomial starting with A1 """ +function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + AsphericSurface(ODD, semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) +end """ @@ -140,9 +133,20 @@ Surface incorporating an aspheric polynomial - radius, conic and aspherics are d `aspherics` should be an array of the both odd and even coefficients of the aspheric polynomial starting with A1 """ +function OddEvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + AsphericSurface(ODDEVEN, semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) +end + +#don't have a function for CONIC as there are better ways to define a CONIC +asphericType(z::AsphericSurface{T,3,Q,M}) where {T<:Real,Q,M} = M # does not seem to work +#asphericType(z::AsphericSurface{T,3,Q,EVEN}) where {T<:Real,Q} = EVEN +#asphericType(z::AsphericSurface{T,3,Q,CONIC}) where {T<:Real,Q} = CONIC +#asphericType(z::AsphericSurface{T,3,Q,ODD}) where {T<:Real,Q} = ODD +#asphericType(z::AsphericSurface{T,3,Q,ODDEVEN}) where {T<:Real,Q} = ODDEVEN -export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface, asphericType +export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface, asphericType, EVEN, CONIC, ODD, ODDEVEN uvrange(::AsphericSurface{T,N}) where {T<:Real,N} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 39706bf9d..b3df4a5d5 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -116,10 +116,10 @@ Either `ZernikeIndexingOSA` or `ZernikeIndexingNoll`, see [Zernike polynomials w export ZernikeIndexType, ZernikeIndexingOSA, ZernikeIndexingNoll """ - ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} + ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N} Surface incorporating the Zernike polynomials - radius, conic and aspherics are defined relative to absolute semi-diameter, Zernike terms are normalized according to the `normradius` parameter. -`T` is the datatype, `N` is the dimensionality, `P` is the number of Zernike terms and `Q` is the number of aspheric terms. Only even aspheric terms are supported. +`T` is the datatype, `N` is the dimensionality, `P` is the number of Zernike terms, `Q` is the number of aspheric terms and `M` is the Aspheric Type. The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space. Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned. @@ -130,7 +130,10 @@ For convenience the input `zcoeff` can be indexed using either OSA or Noll conve ZernikeSurface(semidiameter, radius = Inf, conic = 0, zcoeff = nothing, aspherics = nothing, normradius = semidiameter, indexing = ZernikeIndexingOSA) ``` -`zcoeff` and `aspherics` should be vectors containing tuples of the form `(i, v)` where `i` is either the index of the Zernike term for the corresponding `indexing`, or the polynomial power of the aspheric term (may be even or odd) and `v` is the corresponding coefficient ``A_i`` or ``\\alpha_i`` respectively.. +`zcoeff` and `aspherics` should be vectors containing tuples of the form `(i, v)` where `i` is either the index of the Zernike term +for the corresponding `indexing`, or the polynomial power of the aspheric term (may be even or odd) and +`v` is the corresponding coefficient ``A_i`` or ``\\alpha_i`` respectively.. `M` will be determined from the terms entered to optimize +the evaluation of the aspheric polynomial. The sag is defined by the equation diff --git a/test/Benchmarks/Benchmarks.jl b/test/Benchmarks/Benchmarks.jl index 014222958..d9713a451 100644 --- a/test/Benchmarks/Benchmarks.jl +++ b/test/Benchmarks/Benchmarks.jl @@ -56,13 +56,16 @@ simple_surface_benchmarks() = plane, rectangle, ellipse, triangle, sphere_outsid zernike_surface1() = surfaceintersection, (AcceleratedParametricSurface(TestData.zernikesurface1(), 30), perturbrayz()) zernike_surface2() = surfaceintersection, (AcceleratedParametricSurface(TestData.zernikesurface3(), 30), perturbrayz()) +evenaspheric_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.evenasphericsurface(), 30), perturbrayz()) +oddaspheric_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.oddasphericsurface(), 30), perturbrayz()) +oddevenaspheric_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.oddevenasphericsurface(), 30), perturbrayz()) qtype_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.qtypesurface1(), 30), perturbrayz()) bezier_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.beziersurface(), 30), bezierray()) chebyshev_surface1() = surfaceintersection, (AcceleratedParametricSurface(TestData.chebyshevsurface1(), 30), perturbrayz()) chebyshev_surface2() = surfaceintersection, (AcceleratedParametricSurface(TestData.chebyshevsurface2(), 30), perturbrayz()) gridsag_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.gridsagsurfacebicubic(), 30), perturbrayz()) -accel_surface_benchmarks() = zernike_surface1, zernike_surface2, qtype_surface, bezier_surface, chebyshev_surface1, chebyshev_surface2, gridsag_surface +accel_surface_benchmarks() = zernike_surface1, zernike_surface2, evenaspheric_surface, oddaspheric_surface, oddevenaspheric_surface, qtype_surface, bezier_surface, chebyshev_surface1, chebyshev_surface2, gridsag_surface all_benchmarks() = (simple_surface_benchmarks()..., accel_surface_benchmarks()..., system_benchmarks()...) diff --git a/test/TestData/surfaces.jl b/test/TestData/surfaces.jl index de6f8b2d7..801c3dc13 100644 --- a/test/TestData/surfaces.jl +++ b/test/TestData/surfaces.jl @@ -83,10 +83,11 @@ zernikesurface2() = ZernikeSurface(1.5, radius = 5.0, conic = 1.0, zcoeff = [(7, zernikesurface3() = ZernikeSurface(9.5, zcoeff = [(1, 0.0), (4, -1.0), (7, 0.5)]) -conicsurface() = AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0)) -# -# TODO need to add surfaces to test all types of Aspheric Surfaces -#conicsurface() = AcceleratedParametricSurface(AsphericSurface(9.5, radius = -15.0, conic = -5.0)) +conicsurface() = AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0)) #this also tests the CONIC parameter of AsphericSurface + +oddasphericsurface() = OddAsphericSurface(9.5, 1.0/-50.0, 1.0, [0.0, 0.001, -0.0001]) +evenasphericsurface() = EvenAsphericSurface(9.5, 1.0/-50.0, 1.0, [0.001, -0.0001]) +oddevenasphericsurface() = OddEvenAsphericSurface(9.5, 1.0/-50.0, 1.0, [0.0, 0.001, -0.0001]) function simplesaggrid() points = map( From f4521cd0ffde0b766492bcbb8181a91d8e1e79b9 Mon Sep 17 00:00:00 2001 From: mattderstine Date: Thu, 11 Nov 2021 12:23:28 -0800 Subject: [PATCH 18/19] fixed: documentation, last tests, axisym optisur --- src/Geometry/Primitives/AsphericSurface.jl | 47 ++++++++++------------ src/Optical/Lenses.jl | 2 +- test/testsets/JuliaLang.jl | 1 - 3 files changed, 23 insertions(+), 27 deletions(-) diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl index 05c48132e..ec19bd92d 100644 --- a/src/Geometry/Primitives/AsphericSurface.jl +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -6,10 +6,15 @@ ######################################################################################################### - +#"pseudo-types" of aspheres +""" +AphericSurfaces polynomial evaluation is optimized for no terms `CONIC`, odd terms `ODD`, even terms `EVEN` or any `ODDEVEN` """ +@enum AsphericSurfaceType CONIC ODD EVEN ODDEVEN -AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} + +""" + AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter,. `T` is the datatype, `N` is the dimensionality, `Q` is the number of aspheric terms, and `M` is the type of aspheric polynomial. @@ -32,12 +37,9 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, and ``k = \\texttt{conic}`` . -The function checks if the aspheric terms are missing, even, odd or both and uses the appropriate polynomial evaluation strategy. +The function checks if the aspheric terms are missing, even, odd or both and uses an efficient polynomial evaluation strategy. """ -#"pseudo-types" of aspheres -@enum AsphSurfaceType CONIC ODD EVEN ODDEVEN - struct AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} semidiameter::T curvature::T @@ -46,7 +48,7 @@ struct AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} normradius::T boundingcylinder::Cylinder{T,N} - function AsphericSurface(M::AsphSurfaceType, semidiameter::T, curvature::T, conic::T, aspherics::SVector{Q,T}, normradius::T, boundingcylinder) where {T<:Real,Q} + function AsphericSurface(M::AsphericSurfaceType, semidiameter::T, curvature::T, conic::T, aspherics::SVector{Q,T}, normradius::T, boundingcylinder) where {T<:Real,Q} new{T,3,Q,M}(semidiameter, curvature, conic, aspherics, normradius, boundingcylinder) end @@ -90,10 +92,7 @@ struct AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} end """ - -```julia -EvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) -``` + EvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. @@ -107,10 +106,7 @@ function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics: end """ - -```julia -OddAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) -``` + OddAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. @@ -123,10 +119,7 @@ function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics:: end """ - -```julia -OddEvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) -``` + OddEvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. @@ -138,13 +131,17 @@ function OddEvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspheri AsphericSurface(ODDEVEN, semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) end -#don't have a function for CONIC as there are better ways to define a CONIC +#don't have a function for CONIC as it would better to directly solve for the surface intersection instead of the rootfinding of a ParametricSurface + +""" + asphericType(surf::AsphericSurface) -asphericType(z::AsphericSurface{T,3,Q,M}) where {T<:Real,Q,M} = M # does not seem to work -#asphericType(z::AsphericSurface{T,3,Q,EVEN}) where {T<:Real,Q} = EVEN -#asphericType(z::AsphericSurface{T,3,Q,CONIC}) where {T<:Real,Q} = CONIC -#asphericType(z::AsphericSurface{T,3,Q,ODD}) where {T<:Real,Q} = ODD -#asphericType(z::AsphericSurface{T,3,Q,ODDEVEN}) where {T<:Real,Q} = ODDEVEN +Query the polynomial type of `asp. Returns CONIC, ODD, EVEN, or ODDEVEN. CONIC corresponds to no aspheric terms, ODD +means it only has odd aspheric terms, EVEN means only even aspheric terms and ODDEVEN means both even and odd terms. + +This function is to enable proper interpretation of `surf.aspherics` by any optimization routines that directly query the aspheric coefficients. +""" +asphericType(z::AsphericSurface{T,3,Q,M}) where {T<:Real,Q,M} = M export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface, asphericType, EVEN, CONIC, ODD, ODDEVEN diff --git a/src/Optical/Lenses.jl b/src/Optical/Lenses.jl index 5ea3cdc42..eb81ae196 100644 --- a/src/Optical/Lenses.jl +++ b/src/Optical/Lenses.jl @@ -103,7 +103,7 @@ function AsphericLens(insidematerial::T, frontvertex::S, frontradius::S, frontco if backaspherics !== nothing backaspherics = Tuple{Int,S}.(backaspherics) end - surf = AcceleratedParametricSurface(ZernikeSurface(semidiameter + backdecenter_l + S(0.01), radius = backradius, conic = backconic, aspherics = backaspherics), nsamples, interface = opticinterface(S, insidematerial, nextmaterial, backsurfacereflectance, interfacemode)) + surf = AcceleratedParametricSurface(AsphericSurface(semidiameter + backdecenter_l + S(0.01), radius = backradius, conic = backconic, aspherics = backaspherics), nsamples, interface = opticinterface(S, insidematerial, nextmaterial, backsurfacereflectance, interfacemode)) lens_rear = leaf(surf, Transform{S}(zero(S), S(π), zero(S), backdecenter[1], backdecenter[2], frontvertex - thickness)) end extra_front = frontradius >= zero(S) || isinf(frontradius) ? zero(S) : abs(frontradius) - sqrt(frontradius^2 - semidiameter^2) diff --git a/test/testsets/JuliaLang.jl b/test/testsets/JuliaLang.jl index 268b75916..75f8b4c79 100644 --- a/test/testsets/JuliaLang.jl +++ b/test/testsets/JuliaLang.jl @@ -24,7 +24,6 @@ and pop it from the set (OpticSim.MultiHologramInterface, Tuple{Vararg{HologramInterface{T},N}} where {N} where {T<:Real}), ] : [], OpticSim.Zernike => [], - #OpticSim.AsphericSurface => [], OpticSim.QType => [], OpticSim.Vis => [ (OpticSim.Vis.drawcurves, Tuple{Vararg{Spline{P,S,N,M},N1} where N1} where {M} where {N} where {S} where {P}), From b2bdaa1fde4b3978af6de9ebaea492f3386d8fef Mon Sep 17 00:00:00 2001 From: Charlie Hewitt Date: Thu, 11 Nov 2021 21:32:36 +0000 Subject: [PATCH 19/19] Update ref.md --- docs/src/ref.md | 6 ------ 1 file changed, 6 deletions(-) diff --git a/docs/src/ref.md b/docs/src/ref.md index 9d947cccf..028d09221 100644 --- a/docs/src/ref.md +++ b/docs/src/ref.md @@ -19,12 +19,6 @@ Modules = [OpticSim] Modules = [OpticSim.Geometry] ``` -## Aspheric - -```@autodocs -Modules = [OpticSim.AsphericSurface] -``` - ## Zernike ```@autodocs