diff --git a/src/MultipleScattering.jl b/src/MultipleScattering.jl index e845db00..0bb2de7a 100644 --- a/src/MultipleScattering.jl +++ b/src/MultipleScattering.jl @@ -3,7 +3,7 @@ __precompile__() module MultipleScattering ## Shapes -export Shape, Circle, Rectangle, Box, EmptyShape, Halfspace, TimeOfFlight, TimeOfFlightFromPoint +export Shape, Circle, Rectangle, Box, EmptyShape, Halfspace, Plate, TimeOfFlight, TimeOfFlightFromPoint export outer_radius, volume, name, iscongruent, (≅), congruent, in, issubset, origin, shape, Sphere, (==), isequal, show export boundary_functions, boundary_points, boundary_data, bounding_box, box_corners diff --git a/src/physics/special_functions.jl b/src/physics/special_functions.jl index 96a71e0d..029ef810 100644 --- a/src/physics/special_functions.jl +++ b/src/physics/special_functions.jl @@ -190,7 +190,7 @@ cartesian_to_radial_coordinates(x::Vector) = cartesian_to_radial_coordinates(SVe radial_to_cartesian_coordinates(θ::Vector) = radial_to_cartesian_coordinates(SVector(θ...)) function cartesian_to_radial_coordinates(x::SVector{3,CT}) where CT - r = sqrt(sum(x .^2)) # note this should be complex if x is complex + r = sqrt(sum(x .^2)) # note this is, and should be, a complex number when x is a complex vector θ = atan(sqrt(x[1]^2+x[2]^2),x[3]) φ = atan(x[2], x[1]) return [r,θ,φ] diff --git a/src/shapes/halfspace.jl b/src/shapes/halfspace.jl index 7419ea88..bb8e96cd 100644 --- a/src/shapes/halfspace.jl +++ b/src/shapes/halfspace.jl @@ -1,7 +1,7 @@ """ Halfspace(normal::AbstractVector{T} [, origin::AbstractVector{T}=zeros()]) -N-dimensional [`Shape`](@ref) defined by all points x = origin - λ * normal, for any λ > 0. +A halfspace defined by all the points ``\\mathbf x`` that satify ``(\\mathbf x - \\mathbf o) \\cdot \\mathbf n < 0`` where ``\\mathbf n`` is the unit normal and ``\\mathbf o`` is the origin. """ struct Halfspace{T,Dim} <: Shape{T,Dim} normal::SVector{Dim,T} # outward pointing normal vector diff --git a/src/shapes/plate.jl b/src/shapes/plate.jl new file mode 100644 index 00000000..695fda93 --- /dev/null +++ b/src/shapes/plate.jl @@ -0,0 +1,69 @@ +""" + Plate(normal::AbstractVector{T}, width::T, [, origin::AbstractVector{T}=zeros()]) + +A plate defined by all the points ``\\mathbf x`` that satify ``|(\\mathbf x - \\mathbf o) \\cdot \\mathbf n| < w /2`` where ``\\mathbf n`` is the unit normal, ``\\mathbf o`` is the origin, and ``w`` is the width. +""" +struct Plate{T,Dim} <: Shape{T,Dim} + normal::SVector{Dim,T} # outward pointing normal vector + width::T # the width + origin::SVector{Dim,T} +end + +# Constructors +Plate(normal::AbstractVector{T}, width::T, origin::AbstractVector{T} = zeros(T,length(normal))) where {T} = Plate{T,length(normal)}(normal ./ norm(normal), width, origin) + +Plate(normal::NTuple{Dim,T}, width::T, origin::AbstractVector{T} = zeros(T,Dim)) where {T,Dim} = Plate{T,Dim}(normal ./ norm(normal), width, origin) + +name(shape::Plate) = "Plate" + +volume(shape::Plate) = Inf +outer_radius(hs::Plate{T}) where T = T(Inf) + +# import Base.issubset +# function issubset(inner_rect::Plate{T}, outer_rect::Plate{T}) where T +# all(topright(inner_rect) .<= topright(outer_rect)) && +# all(bottomleft(inner_rect) .>= bottomleft(outer_rect)) +# end + +import Base.in +function in(x::AbstractVector{T}, p::Plate)::Bool where T + abs(dot(x - p.origin, p.normal)) < p.width / T(2) +end + +import Base.(==) +function ==(h1::Plate{T}, h2::Plate{T}) where T + h1.origin == h2.origin && + h1.width == h2.width && + h1.normal == h2.normal +end + +import Base.isequal +function isequal(h1::Plate{T}, h2::Plate{T}) where T + isequal(h1.origin, h2.origin) && + isequal(h1.width, h2.width) && + isequal(h1.normal, h2.normal) +end + +function iscongruent(h1::Plate{T}, h2::Plate{T}) where T + (h1.normal == h2.normal) && (h1.width == h2.width) +end + +function congruent(h::Plate{T}, x) where T + Plate(h.normal, h.width, x) +end + +# function boundary_functions(h::Plate{T,2}, scale::T = 10.0) where T +# surface_vec = [-h.normal[2], h.normal[1]] +# +# function x(t) +# check_boundary_coord_range(t) +# h.origin[1] + surface_vec[1] * width/2 + t * scale * surface_vec[1] +# end +# +# function y(t) +# check_boundary_coord_range(t) +# h.origin[2] + t * scale * surface_vec[2] +# end +# +# return x, y +# end diff --git a/src/shapes/rectangle.jl b/src/shapes/rectangle.jl deleted file mode 100644 index d9134e5a..00000000 --- a/src/shapes/rectangle.jl +++ /dev/null @@ -1,100 +0,0 @@ -""" - Rectangle([origin::AbstractVector{T}=zeros(),] width::T, Height::T) - Rectangle(bottomleft::AbstractVector{T}, topright::AbstractVector{T}) - -2D [`Shape`](@ref) with axis aligned sides, defined by width, height and origin (at the center). -""" -struct Rectangle{T} <: Shape{T,2} - origin::SVector{2,T} - width::T - height::T -end - -# Alternative constructor, where bottomleft and topright corners are specified -function Rectangle(bottomleft::Union{AbstractVector{T},NTuple{2,T}}, topright::Union{AbstractVector{T},NTuple{2,T}}) where {T} - origin = (bottomleft .+ topright)./2 - width = topright[1] - bottomleft[1] - height = topright[2] - bottomleft[2] - if sum(bottomleft .< topright) < 2 - error("bottomleft: $bottomleft should be below and left of topright: $topright.") - end - Rectangle{T}(origin, width, height) -end - -# Alternate constructors, where type is inferred naturally -Rectangle(origin::Tuple{T,T}, width::T, height::T) where {T} = Rectangle{T}(origin, width, height) -Rectangle(origin::Vector{T}, width::T, height::T) where {T} = Rectangle{T}(origin, width, height) -# If no position is given, assume origin is at zero -Rectangle(width::T, height::T) where {T} = Rectangle{T}(SVector(zero(T),zero(T)), width, height) - -name(shape::Rectangle) = "Rectangle" - -outer_radius(r::Rectangle) = sqrt((r.width/2)^2 + (r.height/2)^2) -volume(rectangle::Rectangle) = rectangle.width*rectangle.height - -import Base.issubset -function issubset(inner_rect::Rectangle{T}, outer_rect::Rectangle{T}) where T - all(topright(inner_rect) .<= topright(outer_rect)) && - all(bottomleft(inner_rect) .>= bottomleft(outer_rect)) -end - -import Base.in -function in(x::AbstractVector, r::Rectangle)::Bool - all(abs.(x .- r.origin) .<= SVector(r.width, r.height)) -end - -import Base.(==) -function ==(r1::Rectangle{T}, r2::Rectangle{T}) where T - r1.origin == r2.origin && - r1.width == r2.width && - r1.height == r2.height -end - -import Base.isequal -function isequal(r1::Rectangle{T}, r2::Rectangle{T}) where T - isequal(r1.origin, r2.origin) && - isequal(r1.width , r2.width ) && - isequal(r1.height, r2.height) -end - -function iscongruent(r1::Rectangle{T}, r2::Rectangle{T}) where T - r1.width == r2.width && - r1.height == r2.height -end - -function congruent(r::Rectangle{T}, x) where T - Rectangle{T}(x, r.width, r.height) -end - -# Rectangle bounds itself -bounding_box(rect::Rectangle) = rect - -"Return SVector with the coordinates of the bottom left of a rectangle" -bottomleft(rect::Rectangle) = origin(rect) .- SVector(rect.width/2, rect.height/2) - -"Return SVector with the coordinates of the top right of a rectangle" -topright(rect::Rectangle) = origin(rect) .+ SVector(rect.width/2, rect.height/2) - - -function boundary_functions(rect::Rectangle) - - function x(t) - check_boundary_coord_range(t) - if t <= 1//4 x = bottomleft(rect)[1] + 4t*rect.width - elseif t <= 2//4 x = topright(rect)[1] - elseif t <= 3//4 x = topright(rect)[1] - 4*(t-2//4)*rect.width - else x = bottomleft(rect)[1] - end - end - - function y(t) - check_boundary_coord_range(t) - if t <= 1//4 x = bottomleft(rect)[2] - elseif t <= 2//4 x = bottomleft(rect)[2] + 4*(t-1//4)*rect.height - elseif t <= 3//4 x = topright(rect)[2] - else x = topright(rect)[2] - 4*(t-3//4)*rect.height - end - end - - return x, y -end diff --git a/src/shapes/shapes.jl b/src/shapes/shapes.jl index b932beea..da1074e3 100644 --- a/src/shapes/shapes.jl +++ b/src/shapes/shapes.jl @@ -41,9 +41,11 @@ end include("box.jl") include("sphere.jl") include("halfspace.jl") +include("plate.jl") include("time_of_flight.jl") include("time_of_flight_from_point.jl") include("empty_shape.jl") + """ points_in_shape(Shape; res = 20, xres = res, yres = res, exclude_region = EmptyShape(region), kws...) diff --git a/test/shapetests.jl b/test/shapetests.jl index 6d6cfd03..4f492e1a 100644 --- a/test/shapetests.jl +++ b/test/shapetests.jl @@ -142,4 +142,20 @@ @test true end + + @testset "Infinite volume shapes" begin + normal = rand(3); + origin = rand(3); + width = 1.0 + rand(); + p = Plate(normal, width, origin) + h = Halfspace(normal, origin) + + shapes = [p,h]; + + @test isempty(findall( (volume.(shapes) .== Inf) .== 0)) + + xs = [rand(3) for i = 1:10]; + [[x ∈ s for x in xs] for s in shapes] + + end end