Change AbstractCoordinatePoint to not be a FieldVector

Project.toml

@@ -13,6 +13,7 @@ Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -57,6 +58,7 @@ GPUArrays = "8, 9, 10, 11"
 Geant4 = "0.1.13, 0.2"
 Interpolations = "0.14, 0.15, 0.16"
 IntervalSets = "0.6, 0.7"
+InverseFunctions = "0.1"
 JSON = "0.21.2"
 KernelAbstractions = "0.8, 0.9"
 LaTeXStrings = "1.1"

ext/Geant4/io_gdml.jl

@@ -26,7 +26,7 @@ using LightXML
 # Add <position> to <define> section, referenced in the geometry definition (in <solids>) via the name
 function create_position(e::AbstractConstructiveGeometry, x_define::XMLElement, name::String, v::Bool)
     # Calculate relative position of the volumes contained in the parent volume
-    x, y, z = parse_origin(e.b) .- parse_origin(e.a)
+    x, y, z = parse_origin(e.b) - parse_origin(e.a)
 
     xpos = new_child(x_define, "position")
     set_attributes(xpos, OrderedDict(
@@ -515,9 +515,10 @@ function add_to_world(x_world::XMLElement, x_define::XMLElement, e::AbstractGeom
 
     x_pos_ref = new_child(x_world, "positionref")
     set_attribute(x_pos_ref, "ref", "pos_" * name)
-        
+
     x_pos_ref = new_child(x_define, "position")
-    x, y, z = parse_origin(e)
+    parsed_e = parse_origin(e)
+    x, y, z = parsed_e.x, parsed_e.y, parsed_e.z
     set_attributes(x_pos_ref, OrderedDict(
         "name" => "pos_" * name,
         "x" => x,

src/ChargeDrift/ChargeDrift.jl

@@ -126,7 +126,7 @@ function _add_fieldvector_selfrepulsion!(step_vectors::Vector{CartesianVector{T}
         for m in eachindex(step_vectors)
             if done[m] continue end
             if m > n
-                direction::CartesianVector{T} = current_pos[n] .- current_pos[m]
+                direction::CartesianVector{T} = current_pos[n] - current_pos[m]
                 if iszero(direction) # if the two charges are at the exact same position
                     continue         # don't let them self-repel each other but treat them as same change
                 end                  # if diffusion is simulated, they will very likely be separated in the next step

src/ChargeTrapping/ChargeTrapping.jl

@@ -78,7 +78,7 @@ function _calculate_signal(
     running_sum::T = zero(T)
     q::T = charge
     @inbounds for i in eachindex(tmp_signal)
-        Δldrift::T = (i > 1) ? norm(path[i] .- path[i-1]) : zero(T)
+        Δldrift::T = (i > 1) ? norm(path[i] - path[i-1]) : zero(T)
         Δl::T = Δldrift > 0 ? hypot(Δldrift, vth * (pathtimestamps[i] - pathtimestamps[i-1])) : zero(T)
         Δq::T = q * nσ * Δl
         q -= Δq

src/ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl

@@ -21,6 +21,9 @@ module ConstructiveSolidGeometry
 
     using OrderedCollections: OrderedDict
 
+    import InverseFunctions
+    using InverseFunctions: inverse
+
     import Base: in, *, +, -, &, size, zero
 
     abstract type AbstractCoordinateSystem end
@@ -49,14 +52,25 @@ module ConstructiveSolidGeometry
     abstract type AbstractLinePrimitive{T} <: AbstractPrimitive{T} end
 
     abstract type AbstractConstructiveGeometry{T} <: AbstractGeometry{T} end
+
+    include("Units.jl")
+    include("PointsAndVectors/PointsAndVectors.jl")
+
+    affine_frame(p::AbstractPrimitive) = LocalAffineFrame(origin(p), rotation(p))
 
+    frame_transformation(a::AbstractAffineFrame, b::AbstractPrimitive) = frame_transformation(a, affine_frame(b))
+    frame_transformation(a::AbstractPrimitive, b::AbstractAffineFrame) = frame_transformation(affine_frame(a), b)
+    frame_transformation(a::AbstractPrimitive, b::AbstractPrimitive) = frame_transformation(affine_frame(a), affine_frame(b))
+
     _csg_convert_args(::Type{T}, r::Real) where T = convert(T, r) 
     _csg_convert_args(::Type{T}, r::Tuple) where T = broadcast(x -> _csg_convert_args(T, x), r)
     _csg_convert_args(::Type{T}, r::Nothing) where T = nothing
 
     _csg_get_promoted_eltype(::Type{T}) where {T <: AbstractArray} = eltype(T)
     _csg_get_promoted_eltype(::Type{T}) where {T <: Real} = T
     _csg_get_promoted_eltype(::Type{Nothing}) = Int
+    _csg_get_promoted_eltype(::Type{CartesianPoint{T}}) where {T<:Real} = T
+
     _csg_get_promoted_eltype(::Type{Tuple{T}}) where {T<:Real} = T
     _csg_get_promoted_eltype(::Type{Tuple{T1,T2}}) where {T1<:Real, T2<:Real} = promote_type(T1, T2)
     _csg_get_promoted_eltype(::Type{Tuple{T1,T2}}) where {T1<:Union{Real, Tuple}, T2<:Union{Real, Tuple}} = promote_type(_csg_get_promoted_eltype(T1), _csg_get_promoted_eltype(T2))
@@ -66,17 +80,17 @@ module ConstructiveSolidGeometry
     _handle_phi(φ, rotation) = (φ, rotation)
     _handle_phi(φ::Tuple, rotation) = (abs(φ[2]-φ[1]), rotation*RotZ(φ[1]))
 
-    include("Units.jl")
-    include("PointsAndVectors/PointsAndVectors.jl")
-
     rotation(p::AbstractPrimitive) = p.rotation
     origin(p::AbstractPrimitive) = p.origin
-    _transform_into_global_coordinate_system(pt::CartesianPoint, p::AbstractPrimitive) = (rotation(p) * pt) + origin(p)
-    _transform_into_global_coordinate_system(pt::CartesianVector, p::AbstractPrimitive) = rotation(p) * pt
+
+    # ToDo: Completely replace by frame_transformation and remove:
+    _transform_into_global_coordinate_system(pt::CartesianPoint, p::AbstractPrimitive) = frame_transformation(p, global_frame, pt)
+    _transform_into_global_coordinate_system(v::CartesianVector, p::AbstractPrimitive) = frame_transformation(p, global_frame, v)
     _transform_into_global_coordinate_system(pts::AbstractVector{<:CartesianPoint}, p::AbstractPrimitive) =
-        broadcast(pt -> _transform_into_global_coordinate_system(pt, p), pts)
-    _transform_into_object_coordinate_system(pt::CartesianPoint, p::AbstractPrimitive) = inv(rotation(p)) * (pt - origin(p)) 
-    _transform_into_object_coordinate_system(pt::CartesianVector, p::AbstractPrimitive) = inv(rotation(p)) * pt
+        frame_transformation(p, global_frame, pts)
+    _transform_into_object_coordinate_system(pt::CartesianPoint, p::AbstractPrimitive) = frame_transformation(global_frame, p, pt)
+    _transform_into_object_coordinate_system(v::CartesianVector, p::AbstractPrimitive) = frame_transformation(global_frame, p, v)
+
     in(pt::CartesianPoint{T}, p::AbstractPrimitive{T}, csgtol::T = csg_default_tol(T)) where {T} = 
         _in(_transform_into_object_coordinate_system(pt, p), p; csgtol = csgtol)
     in(pt::CylindricalPoint{T}, p::AbstractPrimitive{T}, csgtol::T = csg_default_tol(T)) where {T} = 

src/ConstructiveSolidGeometry/LinePrimitives/Edge.jl

@@ -25,14 +25,13 @@ function Edge{T}(;
     Edge{T}(a, b)
 end
 
-direction(e::Edge) = e.b - e.a
-
-Line(e::Edge{T}) where {T} = Line{T}(e.a, direction(e))
+Line(e::Edge) = Line(e.a, e.b)
 
 function distance(pt::CartesianPoint{T}, e::Edge{T}) where {T}
-    return if (pt - e.b) ⋅ direction(e) >= 0
+    v = e.b - e.a
+    return if (pt - e.b) ⋅ v >= 0
         norm(pt - e.b)
-    elseif (pt - e.a) ⋅ direction(e) <= 0 
+    elseif (pt - e.a) ⋅ v <= 0 
         norm(pt - e.a)
     else
         distance(pt, Line(e))
@@ -42,9 +41,9 @@ end
 function sample(e::Edge{T}; n = 2)::Vector{CartesianPoint{T}} where {T}
     xs = range(zero(T), stop = one(T), length = n)
     pts = Vector{CartesianPoint{T}}(undef, n)
-    dir = direction(e)
+    v = e.b - e.a
     for i in eachindex(pts)
-        pts[i] = e.a + xs[i] .* dir
+        pts[i] = e.a + xs[i] .* v
     end
     pts
 end

src/ConstructiveSolidGeometry/LinePrimitives/Line.jl

@@ -11,6 +11,13 @@ struct Line{T} <: AbstractLinePrimitive{T}
     end
 end
 
+# function Line(origin::CartesianPoint{T}, direction::CartesianVector{U}) where {T,U}
+#     R = promote_type(T,U)
+#     return Line(convert(CartesianPoint{R}, origin), convert(CartesianVector{R}, direction))
+# end
+
+Line(a::CartesianPoint, b::CartesianPoint) = Line(a, normalize(b - a))
+
 function Line(;
     origin = zero(CartesianPoint{Int}), 
     direction = CartesianVector{Int}(0,0,1)

src/ConstructiveSolidGeometry/PointsAndVectors/AffineFrames.jl

@@ -0,0 +1,66 @@
+frame_transformation(a, b, pt::AbstractCoordinatePoint) = frame_transformation(a, b)(pt)
+frame_transformation(a, b, pts::AbstractVector{<:AbstractCoordinatePoint}) = frame_transformation(a, b).(pts)
+
+frame_transformation(a, b, v::AbstractCoordinateVector) = frame_transformation(a, b)(v)
+frame_transformation(a, b, vs::AbstractVector{<:AbstractCoordinateVector}) = frame_transformation(a, b).(vs)
+
+
+
+abstract type AbstractAffineFrame end
+
+struct GlobalAffineFrame <: AbstractAffineFrame end
+
+const global_frame = GlobalAffineFrame()
+
+frame_transformation(::GlobalAffineFrame, ::GlobalAffineFrame) = identity
+
+
+
+struct LocalAffineFrame{PT<:AbstractCartesianPoint,LM} <: AbstractAffineFrame
+    # origin in the global frame
+    origin::PT
+
+    # linear operator (rotation matrix, etc.) in respect to the global frame
+    linop::LM
+end
+
+
+struct AFTransformLinOpFirst{V<:CartesianVector,LM}
+    linop::LM
+    offset::V
+end
+
+function (f::AFTransformLinOpFirst)(pt::AbstractCartesianPoint)
+    cartesian_zero + muladd(f.linop, pt - cartesian_zero, f.offset) 
+end
+(f::AFTransformLinOpFirst)(v::CartesianVector) = f.linop * v
+
+(f::AFTransformLinOpFirst)(pt::CylindricalPoint) = CylindricalPoint(f(CartesianPoint(pt)))
+
+InverseFunctions.inverse(f::AFTransformLinOpFirst) = AFTransformOffsetFirst(-f.offset, inv(f.linop))
+
+
+struct AFTransformOffsetFirst{V,M}
+    offset::V
+    linop::M
+end
+
+(f::AFTransformOffsetFirst)(pt::AbstractCartesianPoint) =  cartesian_zero + (f.linop * (pt - cartesian_zero + f.offset))
+(f::AFTransformOffsetFirst)(v::CartesianVector) = f.linop * v
+
+(f::AFTransformOffsetFirst)(pt::CylindricalPoint) = CylindricalPoint(f(CartesianPoint(pt)))
+
+InverseFunctions.inverse(f::AFTransformOffsetFirst) = AFTransformLinOpFirst(inv(f.linop), -f.offset)
+
+
+@inline function frame_transformation(frame::LocalAffineFrame, ::GlobalAffineFrame)
+    return AFTransformLinOpFirst(frame.linop, frame.origin - cartesian_zero)
+end
+
+@inline function frame_transformation(::GlobalAffineFrame, frame::LocalAffineFrame)
+    return AFTransformOffsetFirst(cartesian_zero - frame.origin, inv(frame.linop))
+end
+
+function frame_transformation(a::LocalAffineFrame, b::LocalAffineFrame)
+    return frame_transformation(a, global_frame) ∘ frame_transformation(global_frame, b)
+end