Merge pull request JuliaPhysics#441 from JuliaPhysics/dev

Minor update in painting algorithm for `CylindricalGrid`

Author: fhagemann

src/ConstructiveSolidGeometry/SurfacePrimitives/ConeMantle.jl

@@ -230,6 +230,9 @@ function intersection(cm::ConeMantle{T,Tuple{T,T}}, l::Line{T}) where {T}
     end
     ints1 = obj_l.origin + λ1 * obj_l.direction 
     ints2 = obj_l.origin + λ2 * obj_l.direction 
+    # Could be reintroduced as sanity check
+    # ints1 = ints1 * ifelse(abs(ints1.z) <= cm.hZ, one(T), T(NaN))
+    # ints2 = ints2 * ifelse(abs(ints2.z) <= cm.hZ, one(T), T(NaN))
     return _transform_into_global_coordinate_system(ints1, cm), 
            _transform_into_global_coordinate_system(ints2, cm)
 end

src/ConstructiveSolidGeometry/SurfacePrimitives/TorusMantle.jl

@@ -214,33 +214,36 @@ Calculates the intersections of a `Line` with a `TorusMantle`.
 function intersection(tm::TorusMantle{T}, l::Line{T}) where {T}
     obj_l = _transform_into_object_coordinate_system(l, tm) # direction is not normalized
 
-    L1 = obj_l.origin.x
-    L2 = obj_l.origin.y
-    L3 = obj_l.origin.z
-    D1 = obj_l.direction.x
-    D2 = obj_l.direction.y
-    D3 = obj_l.direction.z
+    let T2 = Float64 # will this work on 32bit-machines?
 
-    R = tm.r_torus
-    r = tm.r_tube
+    L1::T2 = obj_l.origin.x |> T2
+    L2::T2 = obj_l.origin.y |> T2
+    L3::T2 = obj_l.origin.z |> T2
+    D1::T2 = obj_l.direction.x |> T2
+    D2::T2 = obj_l.direction.y |> T2
+    D3::T2 = obj_l.direction.z |> T2
 
-    A = L1^2 + L2^2 + L3^2 + R^2 - r^2
-    B = 2*(L1*D1 + L2*D2 + L3*D3)
-    C = D1^2 + D2^2 + D3^2
+    R::T2 = tm.r_torus |> T2
+    r::T2 = tm.r_tube  |> T2
 
-    a = (2*B*C) / C^2
-    b = (2*A*C + B^2 - 4*R^2*(D1^2 + D2^2)) / C^2
-    c = (2*A*B - 8*R^2*(L1*D1 + L2*D2)) / C^2
-    d = (A^2 - 4*R^2*(L1^2 + L2^2)) / C^2
+    A::T2 = L1^2 + L2^2 + L3^2 + R^2 - r^2
+    B::T2 = 2*(L1*D1 + L2*D2 + L3*D3)
+    C::T2 = D1^2 + D2^2 + D3^2
+
+    a::T2 = (2*B*C)
+    b::T2 = (2*A*C + B^2 - 4*R^2*(D1^2 + D2^2))
+    c::T2 = (2*A*B - 8*R^2*(L1*D1 + L2*D2))
+    d::T2 = (A^2 - 4*R^2*(L1^2 + L2^2))
 
     # Solve: `solve (sqrt((L1 + λ*D1)^2 + (L2 + λ*D2)^2)-R)^2 + (L3 + λ*D3)^2 = r^2 for λ`
 
 	# λ1, λ2, λ3, λ4 = roots_of_4th_order_polynomial(a, b, c, d) # That does not work for all combinations of a, b, c, d...
 	# fallback to Polynomials.jl, which is slower... We should improve `roots_of_4th_order_polynomial`... 
     return broadcast(λ -> _transform_into_global_coordinate_system(obj_l.origin + λ * obj_l.direction, tm), 
-        real.(Polynomials.roots(Polynomial((d, c, b, a, one(T))))))
+        T.(real.(Polynomials.roots(Polynomial((d, c, b, a, C^2))))))
     # pts[.!in.(pts, Ref(tm))] .= Ref(CartesianPoint{T}(NaN,NaN,NaN)) # for Partial Tori: discard all intersection points that are not part of the Torus
     # return pts
+    end
 end
 
 # These function compose cross sections of Tori at constant θ and ensure that the height is always positive

src/PotentialCalculation/Painting.jl

@@ -117,13 +117,13 @@ function paint!(point_types, potential, face::AbstractSurfacePrimitive{T}, geome
     eZ = CartesianVector{T}(0,0,1);
 
     widths_ax1 = diff(get_extended_ticks(grid[1]))
-    widths_ax2 = diff(get_extended_ticks(grid[2]))
+    # widths_ax2 = diff(get_extended_ticks(grid[2]))
     widths_ax3 = diff(get_extended_ticks(grid[3]))
     Δw_max_factor = T(1e-5)
     #= 
         Δw_max_factor is chosen by trying out different values for it. 
-        This value seems to be okay if the grid is not to unevenly spaced. 
-        But this can be secured via the `max_ratio`- and the `the max_tick_distance`-keywords.
+        This value seems to be okay if the grid is not too unevenly spaced. 
+        But this can be secured via the `max_ratio`- and the `max_tick_distance`-keywords.
         The critical parameter is `csgtol` inside the `in`-method, which is currently defined through 
         the widths of the voxel of the grid point next to the calculated intersection point.
         It would be probably better to turn this into an `NTuple{3,T}` and pass the widths of the voxel instead
@@ -176,7 +176,8 @@ function paint!(point_types, potential, face::AbstractSurfacePrimitive{T}, geome
                 i1 = searchsortednearest(ticks[1], pt_cyl[1])
                 csgtol = Δw_max_factor * max(
                     widths_ax1[i1], widths_ax1[i1+1],
-                    widths_ax2[i2], widths_ax2[i2+1],
+                    # skip φ (as done in the previous loop)
+                    # widths_ax2[i2], widths_ax2[i2+1]) 
                 ) 
                 if in(pt_cyl, grid) && abs(pt_cyl[2] - ticks[2][i2]) < T(0.1) && in(pt_car, geometry, csgtol)
                     point_types[i1, i2, i3] = zero(PointType)

test/test_real_detectors.jl

@@ -143,7 +143,6 @@ end
 end 
 @timed_testset "Simulate example detector: Toroidal" begin
     sim = Simulation{T}(SSD_examples[:CoaxialTorus])
-    SolidStateDetectors.apply_initial_state!(sim, ElectricPotential)
     timed_simulate!(sim, convergence_limit = 1e-5, device_array_type = device_array_type, refinement_limits = [0.2, 0.1, 0.05, 0.02, 0.01], 
         max_tick_distance = 0.5u"mm", verbose = false)
     evt = Event([CartesianPoint{T}(0.0075,0,0)])