diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml
index a742a9c41..c32013fb9 100644
--- a/.github/workflows/Downgrade.yml
+++ b/.github/workflows/Downgrade.yml
@@ -19,7 +19,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        version: ['1.10', '1']
+        version: ['1.10', '1.12']
     steps:
       - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
diff --git a/Project.toml b/Project.toml
index 9a7e511dd..31da95a0b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,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"
@@ -58,6 +59,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, 1.1"
 KernelAbstractions = "0.8, 0.9"
 LaTeXStrings = "1.1"
@@ -69,7 +71,7 @@ Parameters = "0.12"
 PolygonOps = "0.1.1"
 Polynomials = "2, 3, 4"
 ProgressMeter = "1.5"
-RadiationDetectorSignals = "0.3.5"
+RadiationDetectorSignals = "0.3.10"
 Random = "<0.0.1, 1"
 RecipesBase = "1"
 Requires = "1.1.3"
diff --git a/devenv/README.md b/devenv/README.md
index 10c37661c..ca081232d 100644
--- a/devenv/README.md
+++ b/devenv/README.md
@@ -3,7 +3,8 @@
 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.
+plus BenchmarkTools, Cthulhu and JLD2 (for snapshots and regression tests
+during development).
 
 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
diff --git a/devenv/cuda/Project.toml b/devenv/cuda/Project.toml
index c06810bdc..96f6a1032 100644
--- a/devenv/cuda/Project.toml
+++ b/devenv/cuda/Project.toml
@@ -18,6 +18,7 @@ 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"
diff --git a/docs/Project.toml b/docs/Project.toml
index 64d50ffb1..f1a4c97d4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
diff --git a/docs/src/tutorials/dead_layer_simulation_lit.jl b/docs/src/tutorials/dead_layer_simulation_lit.jl
index 99387fc95..2b6f2ee21 100644
--- a/docs/src/tutorials/dead_layer_simulation_lit.jl
+++ b/docs/src/tutorials/dead_layer_simulation_lit.jl
@@ -11,14 +11,14 @@ using SolidStateDetectors
 T = Float64
 
 # the geometry parameters of the model for following display
-mm=T(1/1000)
-det_rin=1*mm
-det_z=det_r=10*mm
-z_draw=det_z/2
-pn_r=8.957282*mm # this one was calculated by searching the zero impurity point (displayed in the following section)
+mm = T(1/1000)
+det_rin = 1*mm
+det_z = det_r = 10*mm
+z_draw = det_z/2
+pn_r = 8.957282*mm # this one was calculated by searching the zero impurity point (displayed in the following section)
 
 sim = Simulation{T}(SSD_examples[:TrueCoaxial])
-cfn=SSD_examples[:TrueCoaxial]
+cfn = SSD_examples[:TrueCoaxial]
 print(open(f -> read(f, String), cfn))
 plot(sim.detector, xunit = u"mm", yunit = u"mm", zunit = u"mm")
 #jl savefig("tutorial_det_dl.pdf") # hide
@@ -28,10 +28,10 @@ plot(sim.detector, xunit = u"mm", yunit = u"mm", zunit = u"mm")
 
 # ## Display the impurity profile defined
 # > Above position of the PN junction boundary `pn_r` is calculate by searching the point on the curve where the density is zero.
-r_list=(0*mm):(0.01*mm):det_r
+r_list = (0*mm):(0.01*mm):det_r
 imp_list = T[]
 for r in r_list
-    pt=CylindricalPoint{T}(r, 0, z_draw)
+    pt = CylindricalPoint{T}(r, 0, z_draw)
     push!(imp_list, SolidStateDetectors.get_impurity_density(sim.detector.semiconductor.impurity_density_model, pt))
 end
 imp_list = imp_list ./ 1e6 # in cm^-3
@@ -47,12 +47,11 @@ vline!([pn_r/mm], lw = 2, ls = :dash, color = :darkred, label = "PN junction bou
 using SolidStateDetectors: Electron, Hole
 cdm = sim.detector.semiconductor.charge_drift_model
 depth_list = 0:(0.01*mm):(det_r-pn_r)
-hole_mobility_list=[]
-electron_mobility_list=[]
+hole_mobility_list = []
+electron_mobility_list = []
 for depth in depth_list
-    r=det_r-depth
-    pt = CylindricalPoint{T}([r, 0, z_draw])
-    pt = CartesianPoint{T}(pt)
+    r = det_r-depth
+    pt = CartesianPoint{T}(r, 0, z_draw)
     push!(hole_mobility_list, SolidStateDetectors.calculate_mobility(cdm, pt, Hole))
     push!(electron_mobility_list, SolidStateDetectors.calculate_mobility(cdm, pt, Electron))
 end
@@ -70,9 +69,9 @@ plot!(ylabel = "Mobility [cm\$^2\$/V/s]", xlabel = "Depth to surface [mm]",
 calculate_electric_potential!(sim, max_n_iterations = 10, grid = Grid(sim), verbose = false, depletion_handling = true)
 g = sim.electric_potential.grid
 ax1, ax2, ax3 = g.axes
-bulk_tick_dis=0.05*mm
-dl_tick_dis=0.01*mm
-user_additional_ticks_ax1=sort(vcat(ax1.interval.left:bulk_tick_dis:pn_r, pn_r:dl_tick_dis:ax1.interval.right))
+bulk_tick_dis = 0.05*mm
+dl_tick_dis = 0.01*mm
+user_additional_ticks_ax1 = sort(vcat(ax1.interval.left:bulk_tick_dis:pn_r, pn_r:dl_tick_dis:ax1.interval.right))
 user_ax1 = typeof(ax1)(ax1.interval, user_additional_ticks_ax1)
 user_g = typeof(g)((user_ax1, ax2, ax3))
 calculate_electric_potential!(sim, refinement_limits = 0.1, use_nthreads = 8, grid = user_g, depletion_handling = true)
@@ -81,7 +80,7 @@ calculate_weighting_potential!(sim, 1, use_nthreads = 8, depletion_handling = tr
 calculate_weighting_potential!(sim, 2, use_nthreads = 8, depletion_handling = true);
 plot(
     begin
-        imp=plot(sim.imp_scale, φ = 0, xunit = u"mm", yunit = u"mm", title = "impurity scale")
+        imp = plot(sim.imp_scale, φ = 0, xunit = u"mm", yunit = u"mm", title = "impurity scale")
         vline!([pn_r/mm], lw = 2, ls = :dash, color = :darkred, label = "PN junction boundary", legendfontsize = 6)
     end,
     begin
@@ -109,25 +108,25 @@ pulse_list = []
 totTime = 5 # us
 totEnergy = 1 # 1 keV --> simulating ~339 carrier pairs
 max_nsteps = Int(totTime * 1000)
-N=Int(totEnergy*1000÷2.95)
+N = Int(totEnergy*1000÷2.95)
 for depth in depth_list
-    r=det_r-depth
+    r = det_r-depth
     energy_depos = fill(T(2.95/1000), N) * u"keV"
-    starting_positions = repeat([CylindricalPoint{T}(r, deg2rad(0), z_draw)], N)
+    starting_positions = repeat([CartesianPoint{T}(r, deg2rad(0), z_draw)], N)
     evt = Event(starting_positions, energy_depos);
     simulate!(evt, sim, Δt = 1u"ns", max_nsteps = max_nsteps, diffusion = true, end_drift_when_no_field = false, self_repulsion = false)
-    charge=ustrip(evt.waveforms[1].signal)
+    charge = ustrip(evt.waveforms[1].signal)
     push!(pulse_list, charge)
 end
-pulse_plot=plot()
+pulse_plot = plot()
 eff_list = []
 for (i, depth) in enumerate(depth_list)
-    depth=round(depth/mm, digits = 1)
+    depth = round(depth/mm, digits = 1)
     plot!(pulse_list[i], label = "Depth: $(depth) mm", lw = 2, xlabel = "Time [ns]", ylabel = "Amplitude [e]",
         legend = :topright, grid = :on, minorgrid = :on, frame = :box)
     push!(eff_list, maximum(pulse_list[i])/N)
 end
-cce_plot=plot(depth_list/mm, eff_list, xlabel = "Depth to surface [mm]", lw = 2, ylabel = "Charge collection efficiency",
+cce_plot = plot(depth_list/mm, eff_list, xlabel = "Depth to surface [mm]", lw = 2, ylabel = "Charge collection efficiency",
     frame = :box, grid = :on, minorgrid = :on, xticks = 0:0.2:1.2, yticks = 0:0.1:1, dpi = 500, color = :black, label = "")
 plot(pulse_plot, cce_plot, layout = (1, 2), size = (1000, 400), margin = 5Plots.mm)
 #jl savefig("tutorial_pulse_cce_dl.pdf") # hide
diff --git a/docs/src/tutorials/geant4_ssd_lit.jl b/docs/src/tutorials/geant4_ssd_lit.jl
index 9d9eb88de..54c2e6a60 100644
--- a/docs/src/tutorials/geant4_ssd_lit.jl
+++ b/docs/src/tutorials/geant4_ssd_lit.jl
@@ -12,7 +12,7 @@ using Geant4
 # The extension features a function that creates a `Geant4.G4JLApplcation` object from a SSD `Simulation` object and a particle source.
 
 using Plots
-using Unitful
+using ArraysOfArrays, Unitful
 
 # Two types of particle source are pre-defined in `SolidStateDetectors`:
 # 
@@ -89,7 +89,7 @@ events = run_geant4_simulation(app, N_events)
 
 plot(sim.detector, show_passives = false, size = (500,500), fmt = :png)
 plot!(source_1)
-plot!(CartesianPoint.(broadcast(p -> ustrip.(u"m", p), events[1:1000].pos.data)), ms = 0.5, msw = 0, color=:black, label = "")
+plot!(ustrip.(u"m", flatview(events[1:1000].pos)), ms = 0.5, msw = 0, color=:black, label = "")
 #jl savefig("events.pdf") # hide
 #md savefig("events.pdf") # hide
 #md savefig("events.svg"); nothing # hide
diff --git a/ext/Geant4/io_gdml.jl b/ext/Geant4/io_gdml.jl
index 3b48b5e61..24803fa8a 100644
--- a/ext/Geant4/io_gdml.jl
+++ b/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,
diff --git a/ext/SolidStateDetectorsGeant4Ext.jl b/ext/SolidStateDetectorsGeant4Ext.jl
index 019f5028b..feabfb2c1 100644
--- a/ext/SolidStateDetectorsGeant4Ext.jl
+++ b/ext/SolidStateDetectorsGeant4Ext.jl
@@ -219,7 +219,10 @@ function SolidStateDetectors.run_geant4_simulation(app::G4JLApplication, number_
         next!(p)
         evtno += 1
     end
-    return RadiationDetectorSignals.group_by_evtno(Table(evts))
+
+    evts_with_points = map(evt -> merge(evt, (pos = cartesian_zero + evt.pos,)), evts)
+
+    return RadiationDetectorSignals.group_by_evtno(Table(evts_with_points))
 end 
 
 
diff --git a/ext/SolidStateDetectorsLegendHDF5IOExt.jl b/ext/SolidStateDetectorsLegendHDF5IOExt.jl
index bd8c0559e..1f403da2d 100644
--- a/ext/SolidStateDetectorsLegendHDF5IOExt.jl
+++ b/ext/SolidStateDetectorsLegendHDF5IOExt.jl
@@ -6,6 +6,8 @@ import ..LegendHDF5IO
 
 using ..SolidStateDetectors
 using ..SolidStateDetectors: RealQuantity, SSDFloat, to_internal_units, chunked_ranges, LengthQuantity
+using ArraysOfArrays
+import Tables
 using TypedTables, Unitful
 using Format
 
@@ -24,7 +26,7 @@ end
 
 
 
-function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T},
+function SolidStateDetectors.simulate_waveforms( mcevents::AbstractVector{<:NamedTuple}, sim::Simulation{T},
                              output_dir::AbstractString, 
                              output_base_name::AbstractString = "generated_waveforms";
                              chunk_n_physics_events::Int = 1000, 
@@ -42,19 +44,22 @@ function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, si
     Δtime = T(to_internal_units(Δt)) 
     n_contacts = length(sim.detector.contacts)
     @info "Detector has $(n_contacts) contact(s)"
-    contacts = sim.detector.contacts
     if !ispath(output_dir) mkpath(output_dir) end
     nfmt(i::Int) = format(i, zeropadding = true, width = length(digits(n_total_physics_events)))
     evt_ranges = chunked_ranges(n_total_physics_events, chunk_n_physics_events)
-    @info "-> $(sum(length.(mcevents.edep))) energy depositions to simulate."
+    @info "-> $(length(mcevents)) events to simulate."
 
     for evtrange in evt_ranges
         ofn = joinpath(output_dir, "$(output_base_name)_evts_$(nfmt(first(evtrange)))-$(nfmt(last(evtrange))).h5")
         @info "Now simulating $(evtrange) and storing it in\n\t \"$ofn\""
         mcevents_sub = simulate_waveforms(mcevents[evtrange], sim; Δt, max_nsteps, diffusion, self_repulsion, number_of_carriers, number_of_shells, max_interaction_distance, end_drift_when_no_field, geometry_check, verbose)
 
+        # LH5 can't handle CartesianPoint, turn positions into CartesianVectors which will be saved as SVectors
+        pos_vec = VectorOfVectors([[CartesianPoint(p...) - cartesian_zero for p in ps] for ps in mcevents_sub.pos])
+        new_mcevents_sub = TypedTables.Table(merge(Tables.columns(mcevents_sub), (pos = pos_vec,)))
+
         LegendHDF5IO.lh5open(ofn, "w") do h5f
-            LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", mcevents_sub)
+            LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", new_mcevents_sub)
         end        
     end    
 end
diff --git a/src/ChargeDrift/ChargeDrift.jl b/src/ChargeDrift/ChargeDrift.jl
index b5b8084d9..e3e136309 100644
--- a/src/ChargeDrift/ChargeDrift.jl
+++ b/src/ChargeDrift/ChargeDrift.jl
@@ -137,7 +137,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
diff --git a/src/ChargeTrapping/ChargeTrapping.jl b/src/ChargeTrapping/ChargeTrapping.jl
index 7bd006a39..f749bca56 100644
--- a/src/ChargeTrapping/ChargeTrapping.jl
+++ b/src/ChargeTrapping/ChargeTrapping.jl
@@ -112,8 +112,7 @@ function _calculate_signal(
     nσ::T = ifelse(charge > 0, ctm.nσh, ctm.nσe)
     
     @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)
         w::T = i > 1 ? get_interpolation(wpot, path[i], S) : zero(T)
 
@@ -314,7 +313,7 @@ function _calculate_signal(
         Δt::T = (i > 1) ? (pathtimestamps[i] - pathtimestamps[i-1]) : zero(T)
         w::T = i > 1 ? get_interpolation(wpot, path[i], S) : zero(T)
 
-        Δ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)
         
         if in_inactive_region
diff --git a/src/ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl b/src/ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl
index 998dbe2e6..5c62b9054 100644
--- a/src/ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl
+++ b/src/ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl
@@ -19,8 +19,13 @@ module ConstructiveSolidGeometry
     using Unitful
     using YAML
 
+    import StatsBase
+
     using OrderedCollections: OrderedDict
 
+    import InverseFunctions
+    using InverseFunctions: inverse
+
     import Base: in, *, +, -, &, size, zero
 
     abstract type AbstractCoordinateSystem end
@@ -49,7 +54,16 @@ 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
@@ -57,6 +71,8 @@ module ConstructiveSolidGeometry
     _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 +82,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} = 
diff --git a/src/ConstructiveSolidGeometry/LinePrimitives/Edge.jl b/src/ConstructiveSolidGeometry/LinePrimitives/Edge.jl
index af245108d..69f6d4d8d 100644
--- a/src/ConstructiveSolidGeometry/LinePrimitives/Edge.jl
+++ b/src/ConstructiveSolidGeometry/LinePrimitives/Edge.jl
@@ -25,16 +25,13 @@ function Edge{T}(;
     Edge{T}(a, b)
 end
 
-Edge{T}(a::CylindricalPoint{T}, b::CylindricalPoint{T}) where {T} = Edge{T}(CartesianPoint(a), CartesianPoint(b))
-
-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))
@@ -44,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
diff --git a/src/ConstructiveSolidGeometry/LinePrimitives/Line.jl b/src/ConstructiveSolidGeometry/LinePrimitives/Line.jl
index 59957bf4d..52f498b3e 100644
--- a/src/ConstructiveSolidGeometry/LinePrimitives/Line.jl
+++ b/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)
diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/AffineFrames.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/AffineFrames.jl
new file mode 100644
index 000000000..14be8e504
--- /dev/null
+++ b/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
diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index 00cc5db8c..7426c7298 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -1,13 +1,50 @@
-# abstract type AbstractPlanarPoint{T} <: StaticArrays.FieldVector{2, T} end
-# abstract type AbstractPlanarVector{T} <: StaticArrays.FieldVector{2, T} end
-# 
-# struct PlanarPoint{T} <: AbstractPlanarPoint{T}
-#     u::T
-#     v::T
-# end
+"""
+    const AbstractCartesianPoint{T} = AbstractCoordinatePoint{T,Cartesian}
+
+Supertype for cartesian point types.
+"""
+const AbstractCartesianPoint{T} = AbstractCoordinatePoint{T,Cartesian}
+
+Base.:(==)(a::AbstractCartesianPoint, b::AbstractCartesianPoint) = get_x(a) == get_x(b) && get_y(a) == get_y(b) && get_z(a) == get_z(b)
+
+function Base.isapprox(a::AbstractCartesianPoint, b::AbstractCartesianPoint; kwargs...)
+    return isapprox(get_x(a), get_x(b); kwargs...) &&
+           isapprox(get_y(a), get_y(b); kwargs...) &&
+           isapprox(get_z(a), get_z(b); kwargs...)
+end
+
+function to_internal_units(x::AbstractCartesianPoint)
+    CartesianPoint(to_internal_units(get_x(x)), to_internal_units(get_y(x)), to_internal_units(get_z(x)))
+end
+
+
+
+# Unitful uses the `*` and `/` operators to combine values with units. But mathematically that's really not an algebraic
+# product (not defined for affine points), but a cartesian product, so supporting this should be fine:
+Base.:(*)(pt::AbstractCartesianPoint{<:Real}, u::Unitful.Units{<:Any,Unitful.𝐋}) = CartesianPoint(get_x(pt) * u, get_y(pt) * u, get_z(pt) * u)
+Base.:(/)(pt::AbstractCartesianPoint{<:Quantity{<:Real, Unitful.𝐋}}, u::Unitful.Units{<:Any,Unitful.𝐋}) = CartesianPoint(get_x(pt) / u, get_y(pt) / u, get_z(pt) / u)
+
+
+# In that spirit, and to enable constructs like `ustrip(u"mm", pt)` and `NoUnits(pt / u"mm")`, we'll support uconvert/ustrip
+# for CartesianPoint. Unitful doesn't encourage defining those for collections, but we'll view cartesian points as single
+# mathematical objects in regard to units:
+function Unitful.uconvert(u::Unitful.Units{<:Any,Unitful.𝐋}, pt::AbstractCartesianPoint{<:Quantity{<:Real, Unitful.𝐋}})
+    CartesianPoint(uconvert(u, get_x(pt)), uconvert(u, get_y(pt)), uconvert(u, get_z(pt)))
+end
+
+function Unitful.uconvert(u::Unitful.Units{<:Any,Unitful.NoDims}, pt::AbstractCartesianPoint{<:Quantity{<:Real, Unitful.NoDims}})
+    CartesianPoint(uconvert(u, get_x(pt)), uconvert(u, get_y(pt)), uconvert(u, get_z(pt)))
+end
+
+Unitful.ustrip(pt::AbstractCartesianPoint) = CartesianPoint(ustrip(get_x(pt)), ustrip(get_y(pt)), ustrip(get_z(pt)))
+
+
+
+
+
 
 """
-    struct CartesianPoint{T} <: AbstractCoordinatePoint{T, Cartesian}
+    struct CartesianPoint{T} <: AbstractCartesianPoint{T}
 
 Describes a three-dimensional point in Cartesian coordinates.
 
@@ -16,9 +53,12 @@ Describes a three-dimensional point in Cartesian coordinates.
 * `y`: y-coordinate (in m).
 * `z`: z-coordinate (in m).
 
+Given a point `pt = CartesianPoint(x, y, z)`, use `pt - cartesian_zero` to get
+a `CartesianVector` from the origin to point `pt`.
+
 See also [`CylindricalPoint`](@ref).
 """
-struct CartesianPoint{T} <: AbstractCoordinatePoint{T, Cartesian}
+struct CartesianPoint{T} <: AbstractCartesianPoint{T}
     x::T
     y::T
     z::T
@@ -26,31 +66,163 @@ end
 
 #Type promotion happens here
 function CartesianPoint(x::TX, y::TY, z::TZ) where {TX<:Real,TY<:Real,TZ<:Real}
+    # ToDo: Simplify this:
     eltypes = _csg_get_promoted_eltype.((TX,TY,TZ))
     T = float(promote_type(eltypes...))
     CartesianPoint{T}(T(x),T(y),T(z))
 end
 
-function CartesianPoint(;
-    x = 0,
-    y = 0,
-    z = 0
-)
-    CartesianPoint(x,y,z)
+CartesianPoint(; x = 0, y = 0, z = 0) = CartesianPoint(x, y, z)
+
+CartesianPoint{T}(;x = 0, y = 0, z = 0) where {T} = CartesianPoint{T}(T(x),T(y),T(z))
+
+# ToDo: Remove this, if possible
+CartesianPoint{T}(v::AbstractVector) where {T} = cartesian_zero + v
+
+CartesianPoint{T}(pt::CartesianPoint{T}) where {T} = pt
+
+function CartesianPoint{T}(pt::CartesianPoint{U}) where {T, U} 
+    return CartesianPoint{T}(convert(T, pt.x), convert(T, pt.y), convert(T, pt.z))
 end
 
-function CartesianPoint{T}(;
-    x = 0,
-    y = 0,
-    z = 0
-) where {T}
-    CartesianPoint{T}(T(x),T(y),T(z))
+function Base.convert(::Type{CartesianPoint{T}}, pt::CartesianPoint) where {T}
+    return CartesianPoint{T}(pt)
+end
+
+
+@inline get_x(pt::CartesianPoint) = pt.x
+@inline get_y(pt::CartesianPoint) = pt.y
+@inline get_z(pt::CartesianPoint) = pt.z
+
+
+Base.eltype(::Type{<:CartesianPoint{T}}) where {T} = T
+Base.keys(::CartesianPoint) = (:x, :y, :z)
+Base.getindex(p::CartesianPoint, k::Symbol) = getfield(p, k)
+
+AbstractCoordinatePoint{T, Cartesian}(x::Real, y::Real, z::Real) where T = CartesianPoint{T}(x, y, z)
+
+
+Base.zero(::CartesianPoint{T}) where {T} = CartesianPoint{T}(zero(T),zero(T),zero(T))
+# @inline Base.Tuple(pt::CartesianPoint) = (pt.x, pt.y, pt.z)
+@inline Base.copy(pt::CartesianPoint) = CartesianPoint(pt.x, pt.y, pt.z)
+
+@inline Base.:(+)(pt::CartesianPoint, v::CartesianVector) = CartesianPoint(pt.x + v.x, pt.y + v.y, pt.z + v.z)
+@inline Base.:(-)(pt::CartesianPoint, v::CartesianVector) = CartesianPoint(pt.x - v.x, pt.y - v.y, pt.z - v.z)
+
+@inline function Base.:(+)(pt::CartesianPoint, v::AbstractVector)
+    length(v) == 3 || throw(DimensionMismatch("Can only add vectors of length 3 to Cartesian points."))
+    CartesianPoint(pt.x + v[1], pt.y + v[2], pt.z + v[3])
+end
+
+@inline function Base.:(-)(pt::CartesianPoint, v::AbstractVector)
+    length(v) == 3 || throw(DimensionMismatch("Can only add vectors of length 3 to Cartesian points."))
+    CartesianPoint(pt.x - v[1], pt.y - v[2], pt.z - v[3])
 end
 
-zero(PT::Type{<:AbstractCoordinatePoint{T}}) where {T} = PT(zero(T),zero(T),zero(T))
+@inline Base.:(-)(a::CartesianPoint, b::CartesianPoint) = CartesianVector(a.x - b.x, a.y - b.y, a.z - b.z)
+
+@inline Base.transpose(pt::CartesianPoint) = CartesianPoint(transpose(pt.x), transpose(pt.y), transpose(pt.z))
+@inline Base.adjoint(pt::CartesianPoint) = CartesianPoint(adjoint(pt.x), adjoint(pt.y), adjoint(pt.z))
+
+# Barycentric combination
+function  Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}; dims = :) where T
+    # Equivalent to A .- Ref(cartesian_zero), but allocation-free:
+    B = reinterpret(CartesianVector{T}, A) # ToDo: Avoid this memory allocation.
+    return cartesian_zero + mean(B; dims = dims)
+end
+
+function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}, weights::StatsBase.AbstractWeights; dims = :) where T
+    # Equivalent to A .- Ref(cartesian_zero), but allocation-free:
+    B = reinterpret(CartesianVector{T}, A) # ToDo: Avoid this memory allocation.
+    return cartesian_zero + mean(B, weights; dims = dims)
+end
+
+# Base.broadcastable(p::CartesianPoint) = (p.x, p.y, p.z)
+
+
+
+"""
+    CartesianZero{T} <: AbstractCoordinatePoint{T, Cartesian}
+
+Represents origin of the Cartesian coordinate system.
+
+See also [`cartesian_zero`](@ref) and [`CartesianVector`](@ref).
+
+Constructors:
 
+```julia
+CartesianZero{T}()
+CartesianZero() == CartesianZero{Bool}()
+```
 """
-    struct CylindricalPoint{T} <: AbstractCoordinatePoint{T, Cylindrical}
+struct CartesianZero{T} <: AbstractCoordinatePoint{T, Cartesian} end
+CartesianZero() = CartesianZero{Bool}()
+
+@inline get_x(pt::CartesianZero{T}) where T = zero(T)
+@inline get_y(pt::CartesianZero{T}) where T = zero(T)
+@inline get_z(pt::CartesianZero{T}) where T = zero(T)
+
+
+"""
+    const cartesian_zero = CartesianZero()
+
+Origin of the Cartesian coordinate system.
+
+See also [`CartesianZero`](@ref) and [`CartesianVector`](@ref).
+"""
+const cartesian_zero = CartesianZero()
+
+
+@inline Base.:(==)(::CartesianZero, ::CartesianZero) = true
+@inline Base.isapprox(::CartesianZero, ::CartesianZero; kwargs...) = true
+
+@inline Base.:(-)(pt::CartesianPoint, ::CartesianZero) = CartesianVector(pt.x, pt.y, pt.z)
+@inline Base.:(-)(::CartesianZero, pt::CartesianPoint) = CartesianVector(-pt.x, -pt.y, -pt.z)
+
+
+@inline Base.:(+)(::CartesianZero, v::CartesianVector) = CartesianPoint(v.x, v.y, v.z)
+@inline Base.:(-)(::CartesianZero, v::CartesianVector) = CartesianPoint(-v.x, -v.y, -v.z)
+
+@inline Base.:(+)(::CartesianZero, v::StaticVector{3}) = CartesianPoint(v[1], v[2], v[3])
+@inline Base.:(-)(::CartesianZero, v::StaticVector{3}) = CartesianPoint(-v[1], -v[2], -v[3])
+
+@inline function Base.:(+)(::CartesianZero, v::AbstractVector)
+    length(v) == 3 || throw(DimensionMismatch("Can only subtract vectors of length 3 from Cartesian points."))
+    CartesianPoint(v[1], v[2], v[3])
+end
+
+@inline function Base.:(-)(::CartesianZero, v::AbstractVector)
+    length(v) == 3 || throw(DimensionMismatch("Can only add vectors of length 3 to Cartesian points."))
+    CartesianPoint(-v[1], -v[2], -v[3])
+end
+
+@inline function Base.:(-)(::CartesianZero{T}, ::CartesianZero{U}) where {T,U}
+    R = promote_type(T, U)
+    return CartesianVector(zero(R), zero(R), zero(R))
+end
+
+@inline Base.zero(::CartesianZero{T}) where {T} = CartesianZero{T}()
+@inline Base.iszero(::CartesianZero) = true
+
+Base.:(*)(::CartesianZero{T}, u::Unitful.Units{<:Any,Unitful.𝐋}) where {T<:Real} = CartesianZero{Quantity{T, Unitful.𝐋, typeof(u)}}()
+
+function Unitful.uconvert(u::Unitful.Units{T,Unitful.NoDims}, pt::CartesianZero{<:Quantity{U, Unitful.NoDims}}) where {T,U}
+    CartesianZero{promote_type(T,U)}()
+end
+
+
+
+"""
+    const AbstractCylindricalPoint{T} = AbstractCoordinatePoint{T,Cylindrical}
+
+Supertype for cylindrical point types.
+"""
+const AbstractCylindricalPoint{T} = AbstractCoordinatePoint{T,Cylindrical}
+
+
+
+"""
+    struct CylindricalPoint{T} <: AbstractCylindricalPoint{T}
 
 Describes a three-dimensional point in cylindrical coordinates. 
 
@@ -64,7 +236,7 @@ Describes a three-dimensional point in cylindrical coordinates.
     
 See also [`CartesianPoint`](@ref).
 """
-struct CylindricalPoint{T} <: AbstractCoordinatePoint{T, Cylindrical}
+struct CylindricalPoint{T} <: AbstractCylindricalPoint{T}
     r::T
     φ::T
     z::T
@@ -73,26 +245,15 @@ struct CylindricalPoint{T} <: AbstractCoordinatePoint{T, Cylindrical}
 end
 
 function CylindricalPoint(r::TR, φ::TP, z::TZ) where {TR<:Real,TP<:Real,TZ<:Real}
+    # ToDo: Simplify this:
     eltypes = _csg_get_promoted_eltype.((TR,TP,TZ))
     T = float(promote_type(eltypes...))
     CylindricalPoint{T}(T(r),T(φ),T(z))
 end
 
-function CylindricalPoint(;
-    r = 0,
-    φ = 0,
-    z = 0
-)
-    CylindricalPoint(r,φ,z)
-end
+CylindricalPoint(; r = 0, φ = 0, z = 0) = CylindricalPoint(r,φ,z)
 
-function CylindricalPoint{T}(;
-    r = 0,
-    φ = 0,
-    z = 0
-) where {T}
-    CylindricalPoint{T}(T(r),T(φ),T(z))
-end
+CylindricalPoint{T}(; r = 0, φ = 0, z = 0) where {T} = CylindricalPoint{T}(T(r),T(φ),T(z))
 
 function CylindricalPoint(pt::CartesianPoint{T})::CylindricalPoint{T} where {T}
     return CylindricalPoint{T}(hypot(pt.x, pt.y), atan(pt.y, pt.x), pt.z)
@@ -109,6 +270,44 @@ end
 @inline _convert_point(pt::AbstractCoordinatePoint, ::Type{Cylindrical}) = CylindricalPoint(pt)
 @inline _convert_point(pt::AbstractCoordinatePoint, ::Type{Cartesian}) = CartesianPoint(pt)
 
+
+function Base.convert(::Type{CylindricalPoint{T}}, pt::CylindricalPoint{U}) where {T,U}
+    return CylindricalPoint{T}(convert(T, pt.r), convert(T, pt.φ), convert(T, pt.z))
+end
+
+AbstractCoordinatePoint{T, Cylindrical}(r::Real, φ::Real, z::Real) where T = CylindricalPoint{T}(r, φ, z)
+
+Base.:(==)(a::CylindricalPoint, b::CylindricalPoint) = a.r == b.r && a.φ == b.φ && a.z == b.z
+
+function Base.isapprox(a::CylindricalPoint, b::CylindricalPoint; kwargs...)
+    return isapprox(a.r, b.r; kwargs...) && isapprox(a.φ, b.φ; kwargs...) && isapprox(a.z, b.z; kwargs...)
+end
+
+Base.zero(::CylindricalPoint{T}) where {T} = CylindricalPoint{T}(zero(T), zero(T), zero(T))
+Base.iszero(pt::CylindricalPoint) = iszero(pt.r) && iszero(pt.z)
+
+@inline Base.copy(pt::CylindricalPoint) = CylindricalPoint(pt.r, pt.φ, pt.z)
+
+@inline Base.:(+)(pt::CylindricalPoint, v::CartesianVector) = CylindricalPoint(CartesianPoint(pt) + v)
+
+@inline Base.:(-)(pt::CylindricalPoint, v::CartesianVector) = CylindricalPoint(CartesianPoint(pt) - v)
+
+@inline Base.:(-)(a::CylindricalPoint, b::CylindricalPoint) = CartesianPoint(a) - CartesianPoint(b)
+
+@inline Base.transpose(pt::CylindricalPoint) = CylindricalPoint(transpose(pt.r), transpose(pt.φ), transpose(pt.z)) 
+@inline Base.adjoint(pt::CylindricalPoint) = CylindricalPoint(adjoint(pt.r), adjoint(pt.φ), adjoint(pt.z))
+
+# Barycentric combination
+function  Statistics.mean(A::AbstractArray{<:CylindricalPoint}; dims = :)
+    CylindricalPoint(mean(CartesianPoint, A; dims = dims))
+end
+
+function  Statistics.mean(A::AbstractArray{<:CylindricalPoint}, weights::StatsBase.AbstractWeights; dims = :)
+    B = CartesianPoint.(A)
+    CylindricalPoint(mean(B, weights; dims = dims))
+end
+
+
 # function _Δφ(φ1::T, φ2::T)::T where {T}
 #     δφ = mod(φ2 - φ1, T(2π))
 #     min(δφ, T(2π) - δφ)
diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl
index 6b9a176eb..6f57f1d2e 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl
@@ -1,9 +1,21 @@
-abstract type AbstractCoordinatePoint{T, S} <: StaticArrays.FieldVector{3, T} end
+abstract type AbstractCoordinatePoint{T, S} end
+
+zero(PT::Type{<:AbstractCoordinatePoint{T}}) where {T} = PT(zero(T),zero(T),zero(T))
+
+@inline StaticArrays.Size(::PT) where {PT<:AbstractCoordinatePoint} = Size(PT)
+@inline StaticArrays.Size(::Type{PT}) where {PT<:AbstractCoordinatePoint} = Size{(fieldcount(PT),)}()
+@inline Base.length(pt::AbstractCoordinatePoint) = prod(Size(pt))::Int
+@inline Base.size(pt::AbstractCoordinatePoint) = Tuple(Size(pt))
+
+Base.@propagate_inbounds Base.getindex(pt::AbstractCoordinatePoint, i::Int) = getfield(pt, i)
+
+
 abstract type AbstractCoordinateVector{T, S} <: StaticArrays.FieldVector{3, T} end
 
-include("Points.jl")
 include("Vectors.jl")
+include("Points.jl")
+include("AffineFrames.jl")
 
 @inline distance_squared(v::CartesianVector) = v.x * v.x + v.y * v.y + v.z * v.z
-@inline distance_squared(p1::CartesianPoint{T}, p2::CartesianPoint{T}) where {T <: Real} = distance_squared(CartesianVector(p1 .- p2))
-export distance_squared
\ No newline at end of file
+@inline distance_squared(p1::CartesianPoint{T}, p2::CartesianPoint{T}) where {T <: Real} = distance_squared(p1 - p2)
+export distance_squared
diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Vectors.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Vectors.jl
index 82ecd239a..f69f2501c 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Vectors.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Vectors.jl
@@ -13,7 +13,10 @@ Describes a three-dimensional vector in Cartesian coordinates.
 * `y`: y-coordinate (in m).
 * `z`: z-coordinate (in m).
 
-See also [`CylindricalVector`](@ref).
+Given a vector `v = CartesianPoint(Δx, Δy, Δz)`, use `cartesian_zero + v` to
+get the `CartesianPoint(0 + Δx, 0 + Δy, 0 + Δz)`.
+
+See also [`CartesianPoint`](@ref) .
 """
 struct CartesianVector{T} <: AbstractCoordinateVector{T, Cartesian}
     x::T
@@ -23,28 +26,30 @@ end
 
 #Type promotion happens here
 function CartesianVector(x::TX, y::TY, z::TZ) where {TX<:Real,TY<:Real,TZ<:Real}
+    # ToDo: Simplify this:
     eltypes = _csg_get_promoted_eltype.((TX,TY,TZ))
     T = float(promote_type(eltypes...))
     CartesianVector{T}(T(x),T(y),T(z))
 end
 
-function CartesianVector(;
-    x = 0,
-    y = 0,
-    z = 0
-)
-    CartesianVector(x,y,z)
-end
+CartesianVector(; x = 0, y = 0, z = 0) = CartesianVector(x,y,z)
 
-function CartesianVector{T}(;
-    x = 0,
-    y = 0,
-    z = 0
-) where {T}
-    CartesianVector{T}(T(x),T(y),T(z))
-end
+CartesianVector{T}(; x = 0, y = 0, z = 0) where {T} = CartesianVector{T}(T(x),T(y),T(z))
+
+@inline Base.zero(::CartesianVector{T}) where {T} = CartesianVector{T}(zero(T),zero(T),zero(T))
+@inline Base.iszero(v::CartesianVector) = iszero(v.x) && iszero(v.y) && iszero(v.z)
+
+# Need to specialize multiplication and division with units for CartesianVector, otherwise result would just be an SVector:
+Base.:(*)(v::CartesianVector{<:Real}, u::Unitful.Units{<:Any,Unitful.𝐋}) = CartesianVector(v.x * u, v.y * u, v.z * u)
+Base.:(/)(v::CartesianVector{<:Quantity{<:Real, Unitful.𝐋}}, u::Unitful.Units) = CartesianVector(v.x / u, v.y / u, v.z / u)
+
+# For user convenience, to enable constructs like `ustrip(u"mm", v)` and `NoUnits(v / u"mm")`, we'll support uconvert/ustrip
+# for CartesianVector. Unitful doesn't encourage defining those for collections, but we'll view cartesian vectors as single
+# mathematical objects in regard to units:
+Unitful.uconvert(u::Unitful.Units{<:Any,Unitful.𝐋}, v::CartesianVector{<:Quantity{<:Real, Unitful.𝐋}}) = CartesianVector(uconvert(u, v.x), uconvert(u, v.y), uconvert(u, v.z))
+Unitful.uconvert(u::Unitful.Units{<:Any,Unitful.NoDims}, v::CartesianVector{<:Quantity{<:Real, Unitful.NoDims}}) = CartesianVector(uconvert(u, v.x), uconvert(u, v.y), uconvert(u, v.z))
+Unitful.ustrip(v::CartesianVector) = CartesianVector(ustrip(v.x), ustrip(v.y), ustrip(v.z))
 
-zero(VT::Type{<:AbstractCoordinateVector{T}}) where {T} = VT(zero(T),zero(T),zero(T))
 
 # @inline rotate(pt::CartesianPoint{T}, r::RotMatrix{3,T,TT}) where {T, TT} = r.mat * pt
 # @inline rotate(pt::CylindricalPoint{T}, r::RotMatrix{3,T,TT}) where {T, TT} = CylindricalPoint(rotate(CartesianPoint(pt), r))
@@ -60,48 +65,3 @@ zero(VT::Type{<:AbstractCoordinateVector{T}}) where {T} = VT(zero(T),zero(T),zer
 # @inline translate(pt::CylindricalPoint{T}, v::CartesianVector{T}) where {T} = CylindricalPoint(translate(CartesianPoint(pt), v))
 # @inline translate!(vpt::Vector{<:AbstractCoordinatePoint{T}}, v::CartesianVector{T}) where {T} = begin for i in eachindex(vpt) vpt[i] = translate(vpt[i], v) end; vpt end
 # @inline translate!(vvpt::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, v::CartesianVector{T}) where {T} =  begin for i in eachindex(vvpt) translate!(vvpt[i], v) end; vvpt end
-
-
-"""
-    struct CylindricalVector{T} <: AbstractCoordinateVector{T, Cylindrical}
-
-Describes a three-dimensional vector in cylindrical coordinates. 
-
-## Fields
-* `r`: Radius (in m).
-* `φ`: Polar angle (in rad).
-* `z`: `z`-coordinate (in m).
-
-!!! note 
-    `φ == 0` corresponds to the `x`-axis in the Cartesian coordinate system.
-    
-See also [`CartesianVector`](@ref).
-"""
-struct CylindricalVector{T} <: AbstractCoordinateVector{T, Cylindrical}
-    r::T
-    φ::T
-    z::T
-end
-
-#Type promotion happens here
-function CylindricalVector(r::TR, φ::TP, z::TZ) where {TR<:Real,TP<:Real,TZ<:Real}
-    eltypes = _csg_get_promoted_eltype.((TR,TP,TZ))
-    T = float(promote_type(eltypes...))
-    CylindricalVector{T}(T(r),T(φ),T(z))
-end
-
-function CylindricalVector(;
-    r = 0,
-    φ = 0,
-    z = 0
-)
-    CylindricalVector(r,φ,z)
-end
-
-function CylindricalVector{T}(;
-    r = 0,
-    φ = 0,
-    z = 0
-) where {T}
-    CylindricalVector{T}(T(r),T(φ),T(z))
-end
diff --git a/src/ConstructiveSolidGeometry/SurfacePrimitives/ConeMantle.jl b/src/ConstructiveSolidGeometry/SurfacePrimitives/ConeMantle.jl
index 868a6ca84..99f18352d 100644
--- a/src/ConstructiveSolidGeometry/SurfacePrimitives/ConeMantle.jl
+++ b/src/ConstructiveSolidGeometry/SurfacePrimitives/ConeMantle.jl
@@ -76,33 +76,73 @@ radius_at_z(cm::ConeMantle{T,Tuple{T,T}}, z::T) where {T} = radius_at_z(cm.hZ, c
 get_φ_limits(cm::ConeMantle{T,<:Any,T}) where {T} = T(0), cm.φ
 get_φ_limits(cm::ConeMantle{T,<:Any,Nothing}) where {T} = T(0), T(2π)
 
-function normal(cm::ConeMantle{T,T,<:Any,:inwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
-    pto = _transform_into_object_coordinate_system(pt, cm)
-    cyl = CylindricalPoint(pto)
-    return _transform_into_global_coordinate_system(
-            CartesianVector(CartesianPoint(CylindricalPoint{T}(-cyl.r, cyl.φ, zero(T)))), cm)
-end
+# function normal(cm::ConeMantle{T,T,<:Any,:inwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
+#     pto = _transform_into_object_coordinate_system(pt, cm)
+#     cyl = CylindricalPoint(pto)
+#     return _transform_into_global_coordinate_system(
+#             CartesianVector(CartesianPoint(CylindricalPoint{T}(-cyl.r, cyl.φ, zero(T)))), cm)
+# end
+# function normal(cm::ConeMantle{T,T,<:Any,:outwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
+#     pto = _transform_into_object_coordinate_system(pt, cm)
+#     cyl = CylindricalPoint(pto)
+#     return _transform_into_global_coordinate_system(
+#             CartesianVector(CartesianPoint(CylindricalPoint{T}(cyl.r, cyl.φ, zero(T)))), cm)
+# end
+
 function normal(cm::ConeMantle{T,T,<:Any,:outwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
-    pto = _transform_into_object_coordinate_system(pt, cm)
-    cyl = CylindricalPoint(pto)
-    return _transform_into_global_coordinate_system(
-            CartesianVector(CartesianPoint(CylindricalPoint{T}(cyl.r, cyl.φ, zero(T)))), cm)
+    p_local = _transform_into_object_coordinate_system(pt, cm)
+    x, y, _ = p_local
+    normal_local = normalize(CartesianVector(x, y, 0))
+    #FELIX z component ?
+    return _transform_into_global_coordinate_system(normal_local, cm)
 end
+
+function normal(cm::ConeMantle{T,T,<:Any,:inwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
+    p_local = _transform_into_object_coordinate_system(pt, cm)
+    x, y, _ = p_local
+    normal_local = -normalize(CartesianVector(x, y, 0))
+    return _transform_into_global_coordinate_system(normal_local, cm)
+end
+
+# function normal(cm::ConeMantle{T,Tuple{T,T},<:Any,:inwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
+#     pto = _transform_into_object_coordinate_system(pt, cm)
+#     cyl = CylindricalPoint(pto)
+#     Δr = cm.r[2] - cm.r[1]
+#     Δz = 2cm.hZ
+#     return _transform_into_global_coordinate_system(
+#             CartesianVector(CartesianPoint(CylindricalPoint{T}(-one(T), cyl.φ, Δr / Δz))), cm)
+# end
 function normal(cm::ConeMantle{T,Tuple{T,T},<:Any,:inwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
-    pto = _transform_into_object_coordinate_system(pt, cm)
-    cyl = CylindricalPoint(pto)
+    p_local = _transform_into_object_coordinate_system(pt, cm)
+    x = p_local.x
+    y = p_local.y
+    
     Δr = cm.r[2] - cm.r[1]
     Δz = 2cm.hZ
-    return _transform_into_global_coordinate_system(
-            CartesianVector(CartesianPoint(CylindricalPoint{T}(-one(T), cyl.φ, Δr / Δz))), cm)
+    norm_vec_local = normalize(CartesianVector(-x, -y, Δr / Δz))
+    return _transform_into_global_coordinate_system(norm_vec_local, cm)
+
 end
+
+# function normal(cm::ConeMantle{T,Tuple{T,T},<:Any,:outwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
+#     pto = _transform_into_object_coordinate_system(pt, cm)
+#     cyl = CylindricalPoint(pto)
+#     Δr = cm.r[2] - cm.r[1]
+#     Δz = 2cm.hZ
+#     return _transform_into_global_coordinate_system(
+#             CartesianVector(CartesianPoint(CylindricalPoint{T}( one(T), cyl.φ, -Δr / Δz))), cm)
+# end
+
 function normal(cm::ConeMantle{T,Tuple{T,T},<:Any,:outwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T}
-    pto = _transform_into_object_coordinate_system(pt, cm)
-    cyl = CylindricalPoint(pto)
+    p_local = _transform_into_object_coordinate_system(pt, cm)
+    x = p_local.x
+    y = p_local.y
+
     Δr = cm.r[2] - cm.r[1]
     Δz = 2cm.hZ
-    return _transform_into_global_coordinate_system(
-            CartesianVector(CartesianPoint(CylindricalPoint{T}( one(T), cyl.φ, -Δr / Δz))), cm)
+    norm_vec_local = normalize(CartesianVector(x, y, -Δr / Δz))
+    return _transform_into_global_coordinate_system(norm_vec_local, cm)
+
 end
 
 function vertices(cm::ConeMantle{T}, n_arc::Int)::Vector{CartesianPoint{T}} where {T}
diff --git a/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipsoidMantle.jl b/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipsoidMantle.jl
index c160a85d0..50919e1d0 100644
--- a/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipsoidMantle.jl
+++ b/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipsoidMantle.jl
@@ -97,15 +97,15 @@ function normal(em::EllipsoidMantle{T,NTuple{3,T},TP,TT,:outwards}, pt::Cartesia
     # not normalized, do we want this?
     # Or wrap this into somehting like `normal(em, pt) = normalize(direction(em, pt))` ?
     p = _transform_into_object_coordinate_system(pt, em)
-    obj_normal = CartesianPoint{T}(sign(p.x)*(p.x/em.r[1])^2, sign(p.y)*(p.y/em.r[2])^2, sign(p.z)*(p.z/em.r[3])^2) # We might want to store the inv(em.r) in the struct?
-    transform_into_global_coordinate_system(CartesianVector(_obj_normal, em))
+    obj_normal = CartesianVector{T}(sign(p.x)*(p.x/em.r[1])^2, sign(p.y)*(p.y/em.r[2])^2, sign(p.z)*(p.z/em.r[3])^2) # We might want to store the inv(em.r) in the struct?
+    transform_into_global_coordinate_system(obj_normal, em)
 end
 function normal(em::EllipsoidMantle{T,T,TP,TT,:outwards}, pt::CartesianPoint{T})::CartesianVector{T} where {T,TP,TT}
     # not normalized, do we want this?
     # Or wrap this into somehting like `normal(em, pt) = normalize(direction(em, pt))` ?
     p = _transform_into_object_coordinate_system(pt, em)
-    obj_normal = CartesianPoint{T}(sign(p.x)*(p.x/em.r)^2, sign(p.y)*(p.y/em.r)^2, sign(p.z)*(p.z/em.r)^2) # We might want to store the inv(em.r) in the struct?
-    _transform_into_global_coordinate_system(CartesianVector(obj_normal), em)
+    obj_normal = CartesianVector{T}(sign(p.x)*(p.x/em.r)^2, sign(p.y)*(p.y/em.r)^2, sign(p.z)*(p.z/em.r)^2) # We might want to store the inv(em.r) in the struct?
+    _transform_into_global_coordinate_system(obj_normal, em)
 end
 normal(em::EllipsoidMantle{T,TR,TP,TT,:inwards}, pt::CartesianPoint{T}) where {T,TR,TP,TT} = -normal(flip(em), pt)
 
diff --git a/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipticalSurface.jl b/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipticalSurface.jl
index 6e366d01e..2051418cf 100644
--- a/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipticalSurface.jl
+++ b/src/ConstructiveSolidGeometry/SurfacePrimitives/EllipticalSurface.jl
@@ -70,7 +70,7 @@ const Annulus{T} = EllipticalSurface{T,Tuple{T,T},Nothing}
 const PartialAnnulus{T} = EllipticalSurface{T,Tuple{T,T},T}
 
 Plane(es::EllipticalSurface{T}) where {T} = Plane{T}(es.origin, normal(es))
-
+# TODO is this correct?
 normal(es::EllipticalSurface{T}, ::CartesianPoint{T} = zero(CartesianPoint{T})) where {T} = 
     _transform_into_global_coordinate_system(CartesianVector{T}(zero(T), zero(T), one(T)), es)
 
diff --git a/src/ConstructiveSolidGeometry/SurfacePrimitives/Plane.jl b/src/ConstructiveSolidGeometry/SurfacePrimitives/Plane.jl
index 3dbe90466..850ed6eac 100644
--- a/src/ConstructiveSolidGeometry/SurfacePrimitives/Plane.jl
+++ b/src/ConstructiveSolidGeometry/SurfacePrimitives/Plane.jl
@@ -58,7 +58,7 @@ Calculates the intersections of a `Line` with a `Plane`.
 """
 function intersection(p::Plane{T}, line::Line{T}) where {T}
     ndir = normalize(line.direction)
-    λ = (p.normal ⋅ p.origin - p.normal ⋅ line.origin) / (p.normal ⋅ ndir)
+    λ = (p.normal ⋅ (p.origin - line.origin)) / (p.normal ⋅ ndir)
     (line.origin + λ * ndir,)
 end
 
diff --git a/src/ConstructiveSolidGeometry/Transformations.jl b/src/ConstructiveSolidGeometry/Transformations.jl
index 5bdc4376f..0bd059d5e 100644
--- a/src/ConstructiveSolidGeometry/Transformations.jl
+++ b/src/ConstructiveSolidGeometry/Transformations.jl
@@ -6,10 +6,17 @@ const Transformations{T} = NamedTuple{(:rotation, :translation), Tuple{SMatrix{3
 
 # Transformations{T}() where {T} = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(CartesianVector{T}))
 
-rotate(p::P, r::AbstractMatrix) where {P <: AbstractPrimitive} = P(p, origin = r * p.origin, rotation = r * p.rotation)
+# rotate(p::P, r::AbstractMatrix) where {P <: AbstractPrimitive} = P(p, origin = r * p.origin, rotation = r * p.rotation)
+
+function rotate(p::P, r::AbstractMatrix) where {P <: AbstractPrimitive}
+    # Compose rotation with the primitive's local frame:
+    P(p, origin = cartesian_zero + r * (p.origin - cartesian_zero), rotation = r * p.rotation)
+end
+# ToDo: Make this obsolete and then remove it:
 (*)(r::AbstractMatrix, p::AbstractPrimitive) = rotate(p, r)
 
 translate(p::P, v::CartesianVector) where {P <: AbstractPrimitive} = P(p, origin = p.origin + v, rotation = p.rotation)
+# ToDo: Make this obsolete and then remove it:
 (+)(p::AbstractPrimitive, v::CartesianVector) = translate(p, v)
 
 transform(g::AbstractPrimitive, t::Transformations) =
@@ -39,3 +46,5 @@ function Dictionary(m::SMatrix{3,3,T,9}) where {T}
     dict
 end
 
+Dictionary(pt::CartesianPoint) = [pt.x, pt.y, pt.z]
+
diff --git a/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl b/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl
index add510ef6..58a5b5feb 100644
--- a/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl
+++ b/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl
@@ -136,22 +136,22 @@ function Dictionary(b::Box{T})::OrderedDict{String, Any} where {T}
     dict["hX"] = b.hX
     dict["hY"] = b.hY
     dict["hZ"] = b.hZ
-    if b.origin != zero(CartesianVector{T}) dict["origin"] = b.origin end
-    if b.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(b.rotation) end
+    if !iszero(b.origin) dict["origin"] = Dictionary(b.origin) end
+    if !isone(b.rotation) dict["rotation"] = Dictionary(b.rotation) end
     OrderedDict{String,Any}("box" => dict)
 end
 
 function vertices(b::Box{T}) where {T}
-    return (
-        b.rotation * SVector{3,T}(-b.hX, -b.hY, -b.hZ) .+ b.origin,
-        b.rotation * SVector{3,T}(+b.hX, -b.hY, -b.hZ) .+ b.origin,
-        b.rotation * SVector{3,T}(+b.hX, +b.hY, -b.hZ) .+ b.origin,
-        b.rotation * SVector{3,T}(-b.hX, +b.hY, -b.hZ) .+ b.origin,
-        b.rotation * SVector{3,T}(-b.hX, -b.hY, +b.hZ) .+ b.origin,
-        b.rotation * SVector{3,T}(+b.hX, -b.hY, +b.hZ) .+ b.origin,
-        b.rotation * SVector{3,T}(+b.hX, +b.hY, +b.hZ) .+ b.origin,
-        b.rotation * SVector{3,T}(-b.hX, +b.hY, +b.hZ) .+ b.origin,
-    )
+    return _transform_into_global_coordinate_system.((
+        CartesianPoint{T}(-b.hX, -b.hY, -b.hZ),
+        CartesianPoint{T}(+b.hX, -b.hY, -b.hZ),
+        CartesianPoint{T}(-b.hX, +b.hY, -b.hZ),
+        CartesianPoint{T}(+b.hX, +b.hY, -b.hZ),
+        CartesianPoint{T}(-b.hX, -b.hY, +b.hZ),
+        CartesianPoint{T}(+b.hX, -b.hY, +b.hZ),
+        CartesianPoint{T}(-b.hX, +b.hY, +b.hZ),
+        CartesianPoint{T}(+b.hX, +b.hY, +b.hZ)
+    ), Ref(b))
 end
 
 function sample(b::Box{T})::Vector{CartesianPoint{T}} where {T} 
diff --git a/src/ConstructiveSolidGeometry/VolumePrimitives/Cone.jl b/src/ConstructiveSolidGeometry/VolumePrimitives/Cone.jl
index 0dccc5dce..a1884630f 100644
--- a/src/ConstructiveSolidGeometry/VolumePrimitives/Cone.jl
+++ b/src/ConstructiveSolidGeometry/VolumePrimitives/Cone.jl
@@ -716,8 +716,8 @@ function Dictionary(c::Cone{T, <:Any, TR})::OrderedDict{String, Any} where {T, T
     end
     if !isnothing(c.φ) dict["phi"] = OrderedDict("from" => "0°", "to" => string(rad2deg(c.φ))*"°") end
     dict["h"] = 2*c.hZ
-    if c.origin != zero(CartesianVector{T}) dict["origin"] = c.origin end
-    if c.rotation != one(SMatrix{3,3,T,9}) 
+    if !iszero(c.origin) dict["origin"] = Dictionary(c.origin) end
+    if !isone(c.rotation)
         d = Dictionary(c.rotation)
         if unique(keys(d)) == ["Z"]
             φ0 = mod2pi(_parse_value(T, d["Z"], internal_angle_unit))
diff --git a/src/ConstructiveSolidGeometry/VolumePrimitives/Ellipsoid.jl b/src/ConstructiveSolidGeometry/VolumePrimitives/Ellipsoid.jl
index ee0c0a481..715565d34 100644
--- a/src/ConstructiveSolidGeometry/VolumePrimitives/Ellipsoid.jl
+++ b/src/ConstructiveSolidGeometry/VolumePrimitives/Ellipsoid.jl
@@ -134,8 +134,8 @@ function Dictionary(e::Ellipsoid{T})::OrderedDict{String, Any} where {T}
     dict["r"] = e.r # always a Real 
     if !isnothing(e.φ) error("Partial Ellipsoid (`φ = φ`) is not yet supported.") end
     if !isnothing(e.θ) error("Partial Ellipsoid (`θ = θ`) is not yet supported.") end
-    if e.origin != zero(CartesianVector{T}) dict["origin"] = e.origin end
-    if e.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(e.rotation) end
+    if !iszero(e.origin) dict["origin"] = Dictionary(e.origin) end
+    if !isone(e.rotation) dict["rotation"] = Dictionary(e.rotation) end
     OrderedDict{String, Any}("sphere" => dict)
 end
 
diff --git a/src/ConstructiveSolidGeometry/VolumePrimitives/Polycone.jl b/src/ConstructiveSolidGeometry/VolumePrimitives/Polycone.jl
index 189a15cc1..2df022de9 100644
--- a/src/ConstructiveSolidGeometry/VolumePrimitives/Polycone.jl
+++ b/src/ConstructiveSolidGeometry/VolumePrimitives/Polycone.jl
@@ -176,8 +176,8 @@ function Dictionary(c::Polycone{T})::OrderedDict{String, Any} where {T}
     @assert isnothing(c.φ) "Polycone needs `φ` field to be nothing."
     dict["r"] = c.r
     dict["z"] = c.z
-    if c.origin != zero(CartesianVector{T}) dict["origin"] = c.origin end
-    if c.rotation != one(SMatrix{3,3,T,9})  dict["rotation"] = Dictionary(c.rotation) end 
+    if !iszero(c.origin) dict["origin"] = Dictionary(c.origin) end
+    if !isone(c.rotation) dict["rotation"] = Dictionary(c.rotation) end 
     OrderedDict{String, Any}("polycone" => dict)
 end
 
diff --git a/src/ConstructiveSolidGeometry/VolumePrimitives/RegularPrism.jl b/src/ConstructiveSolidGeometry/VolumePrimitives/RegularPrism.jl
index cca86df99..e98ed7342 100644
--- a/src/ConstructiveSolidGeometry/VolumePrimitives/RegularPrism.jl
+++ b/src/ConstructiveSolidGeometry/VolumePrimitives/RegularPrism.jl
@@ -124,8 +124,8 @@ function Dictionary(rp::RegularPrism{N,T, <:Any})::OrderedDict{String, Any} wher
     dict = OrderedDict{String, Any}()
     dict["r"] = rp.r # always a Real 
     dict["h"] = rp.hZ*2
-    if rp.origin != zero(CartesianVector{T}) dict["origin"] = rp.origin end
-    if rp.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(rp.rotation) end
+    if !iszero(rp.origin) dict["origin"] = Dictionary(rp.origin) end
+    if !isone(rp.rotation) dict["rotation"] = Dictionary(rp.rotation) end
     OrderedDict{String, Any}(PrismAliases[N] => dict)
 end
 
diff --git a/src/ConstructiveSolidGeometry/VolumePrimitives/Torus.jl b/src/ConstructiveSolidGeometry/VolumePrimitives/Torus.jl
index a6f400d83..679d1fe76 100644
--- a/src/ConstructiveSolidGeometry/VolumePrimitives/Torus.jl
+++ b/src/ConstructiveSolidGeometry/VolumePrimitives/Torus.jl
@@ -171,8 +171,8 @@ function Dictionary(t::Torus{T,<:Any,TR})::OrderedDict{String, Any} where {T,TR}
     dict["r_tube"] = TR <: Real ? t.r_tube : OrderedDict{String, Any}("from" => t.r_tube[1], "to" => t.r_tube[2]) 
     if !isnothing(t.φ) dict["phi"]   = OrderedDict("from" => "0°", "to" => string(rad2deg(t.φ))*"°") end
     if !isnothing(t.θ) dict["theta"] = OrderedDict("from" => string(rad2deg(t.θ[1]))*"°", "to" => string(rad2deg(t.θ[2]))*"°") end
-    if t.origin != zero(CartesianVector{T}) dict["origin"] = t.origin end
-    if t.rotation != one(SMatrix{3,3,T,9}) 
+    if !iszero(t.origin) dict["origin"] = Dictionary(t.origin) end
+    if !isone(t.rotation) 
         d = Dictionary(t.rotation)
         if unique(keys(d)) == ["Z"]
             φ0 = mod2pi(_parse_value(T, d["Z"], internal_angle_unit))
diff --git a/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl b/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl
index 547c72c76..206a9b20f 100644
--- a/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl
+++ b/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl
@@ -10,7 +10,8 @@ end
     if st == :samplesurface
         CSG_scale = ismissing(CSG_scale) ? get_scale(csg) : CSG_scale
         spacing = T(CSG_scale/n_samples)
-        filter(p -> in(p,csg), sample(primitive, spacing))
+        pts_nounits = filter(p -> in(p,csg), sample(primitive, spacing))
+        pts_nounits .* internal_length_unit
     elseif st == :slice
         isgr = occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
         CSG_scale = ismissing(CSG_scale) ? get_scale(csg) : CSG_scale
diff --git a/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/Points.jl
index 87f158fd3..77f255af3 100644
--- a/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/Points.jl
@@ -1,17 +1,13 @@
 @recipe function f(pt::CartesianPoint)
     xguide --> "x"
-    # xunit --> internal_length_unit
     yguide --> "y"
-    # yunit --> internal_length_unit
     zguide --> "z"
-    # zunit --> internal_length_unit
-    unitformat --> :slash
     if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
         aspect_ratio --> 1.0
-    end 
+    end
     @series begin
         seriestype --> :scatter
-        [pt.x]*internal_length_unit, [pt.y]*internal_length_unit, [pt.z]*internal_length_unit
+        [pt.x], [pt.y], [pt.z]
     end
 end
 
@@ -22,19 +18,24 @@ end
 end
 
 @recipe function f(v::AbstractVector{<:CartesianPoint})
+    default_markersize = 4.0
+    n_thresh = 10
+    n = length(v)
+    invmarkerscale = sqrt(max(1, 2 + (n - 1 - n_thresh) / n_thresh))
+    auto_markersize = max(default_markersize / invmarkerscale, 0.1)
+    auto_markerstrokewidth = auto_markersize < 1 ? zero(auto_markersize) : 0.1 * default_markersize
+
     xguide --> "x"
-    # xunit --> internal_length_unit
     yguide --> "y"
-    # yunit --> internal_length_unit
     zguide --> "z"
-    # zunit --> internal_length_unit
-    unitformat --> :slash
+    markersize --> auto_markersize
+    markerstrokewidth --> auto_markerstrokewidth
     if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
         aspect_ratio --> 1.0
     end 
     @series begin
         if !isempty(v) seriestype --> :scatter end
-        [v[i].x for i in eachindex(v)]*internal_length_unit, [v[i].y for i in eachindex(v)]*internal_length_unit, [v[i].z for i in eachindex(v)]*internal_length_unit
+        [v[i].x for i in eachindex(v)], [v[i].y for i in eachindex(v)], [v[i].z for i in eachindex(v)]
     end
 end
 
@@ -42,4 +43,4 @@ end
     @series begin
         map(x -> CartesianPoint(x), v)
     end
-end
\ No newline at end of file
+end
diff --git a/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/PointsAndVectors.jl b/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/PointsAndVectors.jl
index 798b8900a..3b4c21b2a 100644
--- a/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/PointsAndVectors.jl
+++ b/src/ConstructiveSolidGeometry/plotting/PointsAndVectors/PointsAndVectors.jl
@@ -1,4 +1,2 @@
 include("Points.jl")
 include("Vectors.jl")
-
-
diff --git a/src/Event/Event.jl b/src/Event/Event.jl
index c45d59573..d3943c9d8 100644
--- a/src/Event/Event.jl
+++ b/src/Event/Event.jl
@@ -139,10 +139,12 @@ end
 function move_charges_inside_semiconductor!(
         locations::AbstractVector{<:AbstractVector{CartesianPoint{T}}}, energies::AbstractVector{<:AbstractVector{T}},
         det::SolidStateDetector{T}; fraction::T = T(0.2), verbose::Bool = true) where {T <: SSDFloat}
+    global g_state = (;locations, energies, det)
     for n in eachindex(locations)
         idx_in = broadcast( pt -> pt in det.semiconductor, locations[n]);
         if !all(idx_in)
-            charge_center = sum(locations[n] .* energies[n]) / sum(energies[n])
+            edep_weights = ustrip.(energies[n])
+            charge_center = mean(locations[n], Weights(edep_weights))
             @assert charge_center in det.semiconductor "The center of the charge cloud ($(charge_center)) is not inside the semiconductor."
             surf = ConstructiveSolidGeometry.surfaces(det.semiconductor.geometry)
             for (k,m) in enumerate(findall(.!idx_in))
@@ -160,7 +162,7 @@ function move_charges_inside_semiconductor!(
                     end
                 end
             end
-            charge_center_new = sum(locations[n] .* energies[n]) / sum(energies[n])
+            charge_center_new = mean(locations[n], Weights(edep_weights))
             if verbose
                 @warn "$(sum(.!idx_in)) charges of the charge cloud at $(round.(charge_center, digits = (T == Float64 ? 12 : 6)))"*
                 " are outside. Moving them inside...\nThe new charge center is at $(round.(charge_center_new, digits = (T == Float64 ? 12 : 6))).\n"
diff --git a/src/Grids/Grids.jl b/src/Grids/Grids.jl
index 2db947fb4..a7cad9f71 100644
--- a/src/Grids/Grids.jl
+++ b/src/Grids/Grids.jl
@@ -36,9 +36,13 @@ CartesianGrid3D{T}(a) where {T} = Grid{T, 3, Cartesian, typeof(a)}(a)
 @inline size(grid::Grid{T, N, S}) where {T, N, S} = size.(grid.axes, 1)
 @inline length(grid::Grid{T, N, S}) where {T, N, S} = prod(size(grid))
 @inline getindex(grid::Grid{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = broadcast(getindex, grid.axes, I)
+# @inline getindex(grid::Grid{T, N, S}, I::NTuple{N, Int}) where {T, N, S} = ntuple(i -> grid.axes[i].ticks[inds[i]], Val(N))
 @inline getindex(grid::Grid{T, N, S}, i::Int) where {T, N, S} = getproperty(grid, :axes)[i]
 @inline getindex(grid::Grid{T, N, S}, s::Symbol) where {T, N, S} = getindex(grid, Val{s}())
 
+@inline getpoint(grid::Grid{T, N, S}, idxs::Vararg{Int, N}) where {T, N, S} = AbstractCoordinatePoint{T, S}(grid[idxs...]...)
+@inline getpoint(grid::Grid{T, N, S}, idxs::NTuple{N, Int}) where {T, N, S} = getpoint(grid, idxs...)
+
 @inline getproperty(grid::Grid{T, N, S}, s::Symbol) where {T, N, S} = getproperty(grid, Val{s}())
 @inline getproperty(grid::Grid{T}, ::Val{:axes}) where {T} = getfield(grid, :axes)
 
@@ -57,10 +61,6 @@ CartesianGrid3D{T}(a) where {T} = Grid{T, 3, Cartesian, typeof(a)}(a)
 @inline getindex(grid::CartesianGrid3D{T}, ::Val{:y}) where {T} = @inbounds grid.axes[2]
 @inline getindex(grid::CartesianGrid3D{T}, ::Val{:z}) where {T} = @inbounds grid.axes[3]
 
-@inline GridPoint(grid::Grid{T, 3, Cylindrical}, inds::NTuple{3, Int}) where {T} = 
-    CylindricalPoint{T}(broadcast(i -> grid.axes[i].ticks[inds[i]], (1, 2, 3)))
-@inline GridPoint(grid::Grid{T, 3, Cartesian}, inds::NTuple{3, Int}) where {T} = 
-    CartesianPoint{T}(broadcast(i -> grid.axes[i].ticks[inds[i]], (1, 2, 3)))
 
 function sizeof(grid::Grid{T, N, S}) where {T, N, S}
     return sum( sizeof.(grid.axes) )
diff --git a/src/MCEventsProcessing/MCEventsProcessing.jl b/src/MCEventsProcessing/MCEventsProcessing.jl
index b06205241..e4c89fcc1 100644
--- a/src/MCEventsProcessing/MCEventsProcessing.jl
+++ b/src/MCEventsProcessing/MCEventsProcessing.jl
@@ -1,8 +1,8 @@
 include("table_utils.jl")
 
 """
-    simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T}; kwargs...)
-    simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T}, output_dir::AbstractString, output_base_name::AbstractString; kwargs...)
+    simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::Simulation{T}; kwargs...)
+    simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::Simulation{T}, output_dir::AbstractString, output_base_name::AbstractString; kwargs...)
 
 Simulates the waveforms for all events defined in `mcevents` for a given [`Simulation`](@ref) by
 
@@ -49,7 +49,7 @@ simulate_waveforms(mcevents, sim, "output_dir", "my_basename", Δt = 1u"ns", ver
 !!! note
     Using values with units for `Δt` requires the package [Unitful.jl](https://github.com/PainterQubits/Unitful.jl).
 """
-function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
+function simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::Simulation{T};
                              Δt::RealQuantity = 4u"ns",
                              max_nsteps::Int = 1000,
                              diffusion::Bool = false,
@@ -60,6 +60,7 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
                              end_drift_when_no_field::Bool = true,
                              geometry_check::Bool = false,
                              verbose::Bool = false ) where {T <: SSDFloat}
+    mcevents = TypedTables.Table(mcevents_table)
     n_total_physics_events = length(mcevents)
     Δtime = T(to_internal_units(Δt)) 
     n_contacts = length(sim.detector.contacts)
@@ -100,25 +101,37 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
     return vcat(mcevents_chns...)  
 end
 
+
+# Support detectors hits with positions given as vectors:
+function _get_unitless_positions(point_vectors::AbstractVector{<:StaticVector{3}})
+    Ref(cartesian_zero) .+ to_internal_units.(point_vectors)
+end
+
+# Support detectors hits with positions given as points:
+_get_unitless_positions(points::AbstractVector{<:CartesianPoint}) = to_internal_units.(points)
+
 function _convertEnergyDepsToChargeDeps(
-        pos::AbstractVector{<:SVector{3}}, edep::AbstractVector{<:Quantity}, det::SolidStateDetector{T}; 
+        positions::AbstractVector, edep::AbstractVector{<:Quantity}, det::SolidStateDetector{T}; 
         particle_type::Type{PT} = Gamma, 
         radius::Vector{<:Union{<:Real, <:LengthQuantity}} = radius_guess.(T.(to_internal_units.(edep)), particle_type),
         max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN,
         kwargs...
-    ) where {T <: SSDFloat, PT <: ParticleType} 
+    ) where {T <: SSDFloat, PT <: ParticleType}
 
-    @assert eachindex(pos) == eachindex(edep) == eachindex(radius)
+    @assert eachindex(positions) == eachindex(edep) == eachindex(radius)
 
     d::T = T(to_internal_units(max_interaction_distance))
     @assert isnan(d) || d >= 0 "Max. interaction distance must be positive or NaN (no grouping), but $(max_interaction_distance) was given."
 
-    loc, ene, rad = group_points_by_distance(CartesianPoint.(map(x -> to_internal_units.(x), pos)), T.(to_internal_units.(edep)), T.(to_internal_units.(radius)), d)
+    unitless_pos =_get_unitless_positions(positions)
+    unitless_edep = T.(to_internal_units.(edep))
+    unitless_radius = T.(to_internal_units.(radius))
+    loc, ene, rad = group_points_by_distance(unitless_pos, unitless_edep, unitless_radius, d)
     _convertEnergyDepsToChargeDeps(loc, ene, det; radius = rad, kwargs...)
 end
 
 function _convertEnergyDepsToChargeDeps(
-        pos::AbstractVector{<:AbstractVector}, edep::AbstractVector{<:AbstractVector}, det::SolidStateDetector{T}; 
+        pos::AbstractVector{<:AbstractVector{<:CartesianPoint}}, edep::AbstractVector{<:AbstractVector}, det::SolidStateDetector{T}; 
         particle_type::Type{PT} = Gamma,
         radius::AbstractVector{<:AbstractVector{<:Union{<:Real, <:LengthQuantity}}} = map(e -> radius_guess.(to_internal_units.(e), particle_type), edep),
         number_of_carriers::Int = 1, number_of_shells::Int = 1, 
@@ -129,7 +142,7 @@ function _convertEnergyDepsToChargeDeps(
         iEdep_indep -> broadcast(
             i_together -> begin 
                 nbcc = NBodyChargeCloud(
-                    CartesianPoint(T.(to_internal_units.(pos[iEdep_indep][i_together]))), 
+                    CartesianPoint{T}(to_internal_units(pos[iEdep_indep][i_together])), 
                     T.(to_internal_units(edep[iEdep_indep][i_together])),
                     number_of_carriers,
                     number_of_shells = number_of_shells,
diff --git a/src/PlotRecipes/ChargeCloudModels.jl b/src/PlotRecipes/ChargeCloudModels.jl
index 08fa72ef0..ade6c88eb 100644
--- a/src/PlotRecipes/ChargeCloudModels.jl
+++ b/src/PlotRecipes/ChargeCloudModels.jl
@@ -290,7 +290,6 @@ end
 end
 
 @recipe function f(m::AbstractParticleSource; length = 0.01)
-    
     if hasproperty(m, :opening_angle)
         if iszero(m.opening_angle)
             v = m.position
@@ -307,7 +306,7 @@ end
             b = normalize(a × d)
             rot = hcat(a,b,d)
             cone = SolidStateDetectors.ConstructiveSolidGeometry.Cone(r = ((0,0),(0,length*sin(m.opening_angle))), hZ = length*cos(m.opening_angle)/2, 
-            origin = rot * [0,0,length*cos(m.opening_angle)/2] + m.position, 
+            origin = m.position + rot * CartesianVector(0,0,length*cos(m.opening_angle)/2), 
             rotation = rot)
 
             @series begin
@@ -318,7 +317,7 @@ end
             end
         end
     end
-    
+ 
     @series begin
         seriescolor := :gray
         markersize --> 5
@@ -327,4 +326,4 @@ end
         label --> string(typeof(m).name.name)
         m.position
     end
-end
\ No newline at end of file
+end
diff --git a/src/PlotRecipes/ElectricField.jl b/src/PlotRecipes/ElectricField.jl
index 6f11ea365..f9a268451 100644
--- a/src/PlotRecipes/ElectricField.jl
+++ b/src/PlotRecipes/ElectricField.jl
@@ -165,8 +165,8 @@ end
     spawn_positions = spawn_positions[findall(ipos -> ((spacing-1)+ipos)%spacing == 0, 1:length(spawn_positions))]
 
     for el_field in (el_field_itp, el_field_itp_inv)
-        paths::Array{CartesianPoint{T}, 2} = fill(zero(CartesianVector{T}), length(spawn_positions), max_nsteps)
-        last_step::Int = _drift_charge!(paths, Vector{T}(undef, max_nsteps), sim.detector, sim.point_types, sim.electric_potential.grid, CartesianPoint.(spawn_positions), ones(T, length(spawn_positions)), T(to_internal_units(Δt)), el_field, Electron, geometry_check = false, verbose = false)
+        paths = Array{CartesianPoint{T}, 2}(undef, length(spawn_positions), max_nsteps)
+        last_step::Int = _drift_charge!(paths, Vector{T}(undef, max_nsteps), sim.detector, sim.point_types, sim.electric_potential.grid, CartesianPoint.(spawn_positions), ones(T, length(spawn_positions)), T(to_internal_units(Δt)), el_field, Electron, geometry_check = false, verbose = false )
         for i in 1:size(paths, 1)
             path = @view paths[i, 1:last_step]
             @series begin
diff --git a/src/PotentialCalculation/Painting.jl b/src/PotentialCalculation/Painting.jl
index bca5f911c..6e7727a1c 100644
--- a/src/PotentialCalculation/Painting.jl
+++ b/src/PotentialCalculation/Painting.jl
@@ -167,10 +167,11 @@ function paint!(point_types, potential, face::AbstractSurfacePrimitive{T}, geome
     # end 
 
     for i3 in t_idx_r3 # z;   Maybe switch loops so that the direction of `l` has to be calculated less times..
-        o = CartesianPoint{T}(zero(T), zero(T), ticks[3][i3])
+        z = ticks[3][i3]
+        o = CartesianPoint{T}(zero(T), zero(T), z)
         for i2 in eachindex(ticks[2]) # φ
-            dir = CartesianVector(CartesianPoint(CylindricalPoint{T}(one(T), ticks[2][i2], zero(T)))) 
-            l = ConstructiveSolidGeometry.Line(o, dir) # dir should be normalized
+            radial_point = CartesianPoint(CylindricalPoint{T}(one(T), ticks[2][i2], z))
+            l = ConstructiveSolidGeometry.Line(o, radial_point)
             pts_car = ConstructiveSolidGeometry.intersection(face, l)
             for pt_car in pts_car
                 pt_cyl = CylindricalPoint(pt_car)
diff --git a/src/ScalarPotentials/ScalarPotential.jl b/src/ScalarPotentials/ScalarPotential.jl
index bbd970b1c..98e93d1d8 100644
--- a/src/ScalarPotentials/ScalarPotential.jl
+++ b/src/ScalarPotentials/ScalarPotential.jl
@@ -26,11 +26,10 @@ function getindex(ep::P, grid::Grid{T, N, S}) where {T, N, S, P <: ScalarPotenti
     gridsize::Tuple = size(grid)
     data::Array{T, N} = zeros(T, gridsize)
     ep_itp::Interpolations.Extrapolation{T, N} = interpolated_scalarfield(ep)
-    PT = (S == Cylindrical ? CylindricalPoint : CartesianPoint)
     for i1 in eachindex(grid[1])
         for i2 in eachindex(grid[2])
             for i3 in eachindex(grid[3])
-                data[i1, i2, i3] = get_interpolation(ep_itp, PT{T}(grid[i1, i2, i3]), S)
+                data[i1, i2, i3] = get_interpolation(ep_itp, getpoint(grid, i1, i2, i3), S)
             end
         end
     end
diff --git a/src/Simulation/Capacitance.jl b/src/Simulation/Capacitance.jl
index 0e849ffd2..33014e3d4 100644
--- a/src/Simulation/Capacitance.jl
+++ b/src/Simulation/Capacitance.jl
@@ -58,7 +58,7 @@ function _calculate_mutual_capacitance(grid::Grid{T, 3, CS}, grid_mps, int_ϵ_r,
                 w1, w2, w3 = voxel_widths(grid, i1, i2, i3)
                 dV = voxel_volume(grid, i1, i2, i3, w1, w2, w3)
 
-                pt_voxel_mid = GridPoint(grid_mps, (i1 + 1, i2 + 1, i3 + 1))
+                pt_voxel_mid = getpoint(grid_mps, (i1 + 1, i2 + 1, i3 + 1))
                 ϵ_r_voxel = get_interpolation(int_ϵ_r, pt_voxel_mid, CS)
 
                 efs_1 = _approximate_potential_gradient(int_p1, grid, i1, i2, i3, w1, w2, w3) 
@@ -72,14 +72,14 @@ function _calculate_mutual_capacitance(grid::Grid{T, 3, CS}, grid_mps, int_ϵ_r,
 end
 
 function _approximate_potential_gradient(int_p, grid::Grid{T, 3, CS}, i1, i2, i3, w1, w2, w3) where {T, CS}
-    p000 = get_interpolation(int_p, GridPoint(grid, (i1    , i2    , i3    )), CS)
-    p100 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2    , i3    )), CS)
-    p010 = get_interpolation(int_p, GridPoint(grid, (i1    , i2 + 1, i3    )), CS)
-    p110 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2 + 1, i3    )), CS)
-    p001 = get_interpolation(int_p, GridPoint(grid, (i1    , i2    , i3 + 1)), CS)
-    p101 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2    , i3 + 1)), CS)
-    p011 = get_interpolation(int_p, GridPoint(grid, (i1    , i2 + 1, i3 + 1)), CS)
-    p111 = get_interpolation(int_p, GridPoint(grid, (i1 + 1, i2 + 1, i3 + 1)), CS)
+    p000 = get_interpolation(int_p, getpoint(grid, (i1    , i2    , i3    )), CS)
+    p100 = get_interpolation(int_p, getpoint(grid, (i1 + 1, i2    , i3    )), CS)
+    p010 = get_interpolation(int_p, getpoint(grid, (i1    , i2 + 1, i3    )), CS)
+    p110 = get_interpolation(int_p, getpoint(grid, (i1 + 1, i2 + 1, i3    )), CS)
+    p001 = get_interpolation(int_p, getpoint(grid, (i1    , i2    , i3 + 1)), CS)
+    p101 = get_interpolation(int_p, getpoint(grid, (i1 + 1, i2    , i3 + 1)), CS)
+    p011 = get_interpolation(int_p, getpoint(grid, (i1    , i2 + 1, i3 + 1)), CS)
+    p111 = get_interpolation(int_p, getpoint(grid, (i1 + 1, i2 + 1, i3 + 1)), CS)
     efv1 = ( (p100 - p000) + (p110 - p010) + (p101 - p001) + (p111 - p011) ) / (4 * w1)
     efv2 = if CS == Cylindrical
         _w2 = (grid.axes[2].ticks[i2 + 1] - grid.axes[2].ticks[i2])
diff --git a/src/Simulation/Simulation.jl b/src/Simulation/Simulation.jl
index c590a7a86..1dee55e2b 100644
--- a/src/Simulation/Simulation.jl
+++ b/src/Simulation/Simulation.jl
@@ -849,7 +849,7 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
         else
             sor_consts = T.(sor_consts)
         end
-        min_tick_distance::NTuple{3, T} = if CS == Cylindrical
+        new_min_tick_distance::NTuple{3, T} = if CS == Cylindrical
             if !ismissing(min_tick_distance)
                 if min_tick_distance isa LengthQuantity
                     world_r_mid = (sim.world.intervals[1].right + sim.world.intervals[1].left)/2
@@ -962,7 +962,7 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                 else
                     abs.(ref_limits .* bias_voltage)
                 end
-                refine!(sim, ElectricPotential, max_diffs, min_tick_distance)
+                refine!(sim, ElectricPotential, max_diffs, new_min_tick_distance)
                 nt = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.electric_potential.grid), max_nthreads[iref+1], CS) : max_nthreads[iref+1]
                 verbose && println("Grid size: $(size(sim.electric_potential.data)) - $(onCPU ? "using $(nt) threads now" : "GPU")") 
                 update_till_convergence!( sim, potential_type, convergence_limit,
@@ -976,7 +976,7 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                                                 sor_consts = is_last_ref ? T(1) : sor_consts )
             else
                 max_diffs = abs.(ref_limits)
-                refine!(sim, WeightingPotential, contact_id, max_diffs, min_tick_distance)
+                refine!(sim, WeightingPotential, contact_id, max_diffs, new_min_tick_distance)
                 nt = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.weighting_potentials[contact_id].grid), max_nthreads[iref+1], CS) : max_nthreads[iref+1]
                 verbose && println("Grid size: $(size(sim.weighting_potentials[contact_id].data)) - $(onCPU ? "using $(nt) threads now" : "GPU")") 
                 update_till_convergence!( sim, potential_type, contact_id, convergence_limit,
diff --git a/src/SolidStateDetectors.jl b/src/SolidStateDetectors.jl
index 212800346..d5116866a 100644
--- a/src/SolidStateDetectors.jl
+++ b/src/SolidStateDetectors.jl
@@ -35,17 +35,19 @@ include("ka_compat.jl")
 include("ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl")
 using .ConstructiveSolidGeometry
 using .ConstructiveSolidGeometry:
-            CylindricalPoint, CartesianPoint, AbstractCoordinatePoint, _convert_point,
-            CartesianVector, CylindricalVector, AbstractCoordinateVector,
+            CartesianPoint, CartesianVector, CartesianZero, cartesian_zero, CylindricalPoint,
+            AbstractCoordinatePoint, _convert_point,
+            CartesianVector, AbstractCoordinateVector,
             Cartesian, Cylindrical, AbstractCoordinateSystem, CoordinateSystemType,
             Geometry, AbstractGeometry, AbstractSurfacePrimitive,
             parse_rotation_matrix, parse_translate_vector, parse_CSG_transformation,
             transform, CSG_dict, Transformations, combine_transformations,
             ConfigFileError, _parse_value,
-            LengthQuantity, AngleQuantity, get_scale
+            LengthQuantity, AngleQuantity, get_scale,
+            LocalAffineFrame, cartesian_zero, global_frame, frame_transformation
         
 import .ConstructiveSolidGeometry: sample, to_internal_units, from_internal_units
-export CartesianPoint, CartesianVector, CylindricalPoint
+export CartesianPoint, CartesianVector, CartesianZero, cartesian_zero, CylindricalPoint
 
 import Clustering
 import DataStructures
@@ -64,7 +66,7 @@ import Base: size, sizeof, length, getindex, setindex!, axes, getproperty, broad
 export SolidStateDetector
 export SSD_examples
 
-export Grid, GridPoint
+export Grid
 
 export ElectricPotential, PointTypes, EffectiveChargeDensity, DielectricDistribution, WeightingPotential, ElectricField
 export apply_initial_state!
@@ -128,6 +130,11 @@ include("PlotRecipes/PlotRecipes.jl")
 export @P_str # protected strings to overwrite plot labels with units
 
 
+export GridPoint
+@deprecate GridPoint(grid::Grid{T, N, S}, idxs::Vararg{Int, N}) where {T, N, S} getpoint(grid, idxs...)
+@deprecate GridPoint(grid::Grid{T, N, S}, idxs::NTuple{N, Int}) where {T, N, S} getpoint(grid, idxs)
+
+
 function __init__()
     @require LegendHDF5IO ="c9265ca6-b027-5446-b1a4-febfa8dd10b0" begin
         include("../ext/SolidStateDetectorsLegendHDF5IOExt.jl")
diff --git a/test/ConstructiveSolidGeometry/CSG_coordinates.jl b/test/ConstructiveSolidGeometry/CSG_coordinates.jl
new file mode 100644
index 000000000..cb8104953
--- /dev/null
+++ b/test/ConstructiveSolidGeometry/CSG_coordinates.jl
@@ -0,0 +1,72 @@
+# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).
+
+using Test
+
+using SolidStateDetectors.ConstructiveSolidGeometry: CartesianPoint, CartesianVector, 
+    CartesianZero, cartesian_zero, CylindricalPoint, LocalAffineFrame, global_frame, frame_transformation
+using StaticArrays: Size, SVector, SMatrix
+using InverseFunctions: inverse
+import Unitful
+
+@testset "points_and_vectors" begin
+    @testset "cartesian" begin
+        cart = @inferred CartesianVector(x=2f0,z=1f0)
+        @inferred CartesianVector{Float32}(x=2)
+        @test cart.x == Float32(2)
+        cart = @inferred CartesianPoint(x=2f0,z=1f0)
+        @inferred CartesianPoint{Float32}(x=2)
+
+        a = CartesianPoint(1.0, 2.0, 3.0)
+        b = CartesianPoint(3.0, 1.0, 2.0)
+        v = CartesianVector(0.1, 0.2, 0.3)
+        A = SMatrix{3,3}(1.0, 4.0, 7.0, 4.0, 5.0, 8.0, 7.0, 8.0, 9.0)
+
+        @test @inferred(a + v) == CartesianPoint(1.1, 2.2, 3.3)
+        @test @inferred(a - v) == CartesianPoint(0.9, 1.8, 2.7)
+        @test @inferred(a - v) == CartesianPoint(0.9, 1.8, 2.7)
+        @test @inferred(a - b) == CartesianVector(-2.0, 1.0, 1.0)
+
+        @test @inferred(zero(a) + (a - zero(a))) == a
+
+        @test @inferred(Size(a)) === Size(v)
+        @test @inferred(size(a)) === size(v)
+        @test @inferred(length(a)) === length(v)
+        @test @inferred(CartesianPoint(a[1], a[2], a[3])) === a
+        @test @inferred(CartesianPoint(a[1], a[2], a[3])) == a
+        @test @inferred(CartesianPoint(a[1], a[2], a[3])) ≈ a
+
+
+        frame = LocalAffineFrame(b, A)
+
+        f = frame_transformation(frame, global_frame)
+        @test @inferred(f(a)) == cartesian_zero + A * (a - cartesian_zero) + CartesianVector(b[1], b[2], b[3])
+        @test @inferred(inverse(f)(f(a))) ≈ a
+
+
+        @test @inferred(CartesianZero{Float32}() * Unitful.u"mm") === CartesianZero{typeof(zero(Float32) * Unitful.u"mm")}()
+        # ToDo: Add more tests!
+    end
+
+    @testset "cylindrical" begin
+        cyl = @inferred CylindricalPoint{Float32}(r=2.,z=1.)
+        cyl2 = @inferred CylindricalPoint(φ=3π)
+        @test CartesianPoint(cyl) == CartesianPoint(x=2f0,z=1f0)
+
+        a = CylindricalPoint(1.0, π/2, 3.0)
+        b = CylindricalPoint(3.0, π/2, 2.0)
+        v = CartesianVector(0.0, 0.2, 0.3)
+
+        @test  @inferred(a + v) ≈ CylindricalPoint(1.2, π/2, 3.3)
+        @test  @inferred(a - v) ≈ CylindricalPoint(0.8, π/2, 2.7)
+        @test  @inferred(a - b) ≈ CartesianVector(0, -2, 1)
+
+        @test @inferred(zero(a) + (a - zero(a))) == a
+
+        @test  @inferred(Size(a)) === Size(v)
+        @test  @inferred(size(a)) === size(v)
+        @test  @inferred(length(a)) === length(v)
+        @test  @inferred(CylindricalPoint(a[1], a[2], a[3])) === a
+        @test  @inferred(CylindricalPoint(a[1], a[2], a[3])) == a
+        @test  @inferred(CylindricalPoint(a[1], a[2], a[3])) ≈ a
+    end
+end
diff --git a/test/ConstructiveSolidGeometry/CSG_primitives.jl b/test/ConstructiveSolidGeometry/CSG_primitives.jl
index 005f297d0..90e75c6be 100644
--- a/test/ConstructiveSolidGeometry/CSG_primitives.jl
+++ b/test/ConstructiveSolidGeometry/CSG_primitives.jl
@@ -23,7 +23,7 @@ no_translations = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(Carte
                     "h"       => 1.0))
                 for i in 1:2
             ])
-    
+            
             c = Geometry(T, dict, default_units, no_translations)
     
             # Conversion from Geometry -> Dict and Dict -> Geometry should result in the same geometry
@@ -195,8 +195,8 @@ no_translations = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(Carte
         @test in(CartesianPoint{Float32}(1,0,0),ellip1)
         @test !in(CartesianPoint{Float32}(1,0,0),ellip_open)
     
-        ellip_closed_trafo = @inferred CSG.Ellipsoid(CSG.ClosedPrimitive,r=1.0, origin=1/sqrt(3)*CartesianPoint{Float32}(1,1,1),rotation=SMatrix{3}(0.5,sqrt(3)/2,0,-sqrt(3)/2,0.5,0,0,0,1))
-        ellip_open_trafo = @inferred CSG.Ellipsoid(CSG.OpenPrimitive,r=1.0, origin=1/sqrt(3)*CartesianPoint{Float32}(1,1,1),rotation=SMatrix{3}(0.5,sqrt(3)/2,0,-sqrt(3)/2,0.5,0,0,0,1))
+        ellip_closed_trafo = @inferred CSG.Ellipsoid(CSG.ClosedPrimitive,r=1.0, origin=CartesianPoint{Float64}(1/sqrt(3),1/sqrt(3),1/sqrt(3)),rotation=SMatrix{3}(0.5,sqrt(3)/2,0,-sqrt(3)/2,0.5,0,0,0,1))
+        ellip_open_trafo = @inferred CSG.Ellipsoid(CSG.OpenPrimitive,r=1.0, origin=CartesianPoint{Float64}(1/sqrt(3),1/sqrt(3),1/sqrt(3)),rotation=SMatrix{3}(0.5,sqrt(3)/2,0,-sqrt(3)/2,0.5,0,0,0,1))
         @test !in(CartesianPoint{Float64}(-1e-8,0,0),ellip_closed_trafo)
         @test in(CartesianPoint{Float64}(1e-8,0,0),ellip_open_trafo)
     end
@@ -299,22 +299,6 @@ no_translations = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(Carte
         edge2 = @inferred CSG.Edge(a = CartesianPoint{Float32}(0,0,0), b = CartesianPoint{Float16}(0,0,1)) 
         @test edge1 == edge2
     end
-    @testset "Vector" begin
-        cart = @inferred CSG.CartesianVector(x=2f0,z=1f0)
-        @inferred CSG.CartesianVector{Float32}(x=2)
-        cyl = @inferred CSG.CylindricalVector{Float32}(r=2.,z=1.)
-        cyl2 = @inferred CSG.CylindricalVector(φ=3π)
-        @test CartesianVector(cyl) == cart
-        @test cart.x == Float32(2)
-    end
-    @testset "Point" begin
-        cart = @inferred CSG.CartesianPoint(x=2f0,z=1f0)
-        @inferred CSG.CartesianPoint{Float32}(x=2)
-        cyl = @inferred CSG.CylindricalPoint{Float32}(r=2.,z=1.)
-        cyl2 = @inferred CSG.CylindricalPoint(φ=3π)
-        @test CartesianPoint(cyl) == cart
-        @test cyl2.φ == T(π)
-    end
 end
 
 @testset "Test geometry parsing" begin
diff --git a/test/Project.toml b/test/Project.toml
index b9649e3af..9b4965708 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,6 +2,7 @@
 ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8"
 LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/test/runtests.jl b/test/runtests.jl
index 08afeeac2..4f6e581f7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -17,6 +17,7 @@ end
 
 @timed_testset "ConstructiveSolidGeometry" begin
     include("ConstructiveSolidGeometry/CSG_IO.jl")
+    include("ConstructiveSolidGeometry/CSG_coordinates.jl")
     include("ConstructiveSolidGeometry/CSG_primitives.jl")
 end