Change AbstractCoordinatePoint to not be a FieldVector

Implement proper math for affine frames, Refactor code base where necessary.

Author: oschulz

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"

devenv/README.md

@@ -0,0 +1,11 @@
+# Julia development environments
+
+This directory contains the Julia project environments `cpu` and `cuda`
+that can be used when developing SolidStateDetectors with Julia >= v1.11.
+They contain all direct, test and doc-gen dependencies of SolidStateDetectors,
+plus BenchmarkTools and Cthulhu.
+
+Note: These environments can't be used with Julia versions <= v1.10, as they
+use a `[sources]` section in the `Project.toml` to ensure SolidStateDetectors
+is loaded from the local source directory, this Pkg feature was introduced
+in Julia v1.11.

devenv/cpu/Project.toml

@@ -0,0 +1,55 @@
+[deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+AffineMaps = "2c83c9a8-abf5-4329-a0d7-deffaf474661"
+ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
+Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
+Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44"
+HeterogeneousComputing = "2182be2a-124f-4a91-8389-f06db5907a21"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8"
+LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0"
+LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+ParallelProcessingTools = "8e8a01fc-6193-5ca1-a2f1-20776dae4199"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+RadiationDetectorSignals = "bf2c0563-65cf-5db2-a620-ceb7de82658c"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+SolidStateDetectors = "71e43887-2bd9-5f77-aebd-47f656f0a3f0"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
+YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+
+[sources]
+SolidStateDetectors = {path = "../.."}

devenv/cuda/Project.toml

@@ -0,0 +1,55 @@
+[deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+AffineMaps = "2c83c9a8-abf5-4329-a0d7-deffaf474661"
+ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
+Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
+Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44"
+HeterogeneousComputing = "2182be2a-124f-4a91-8389-f06db5907a21"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+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"
+LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8"
+LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0"
+LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+ParallelProcessingTools = "8e8a01fc-6193-5ca1-a2f1-20776dae4199"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+RadiationDetectorSignals = "bf2c0563-65cf-5db2-a620-ceb7de82658c"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+SolidStateDetectors = "71e43887-2bd9-5f77-aebd-47f656f0a3f0"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
+YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+
+[sources]
+SolidStateDetectors = {path = "../.."}

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)