Code cleanup

src/specializations/Ray.jl

@@ -45,7 +45,8 @@ function integral(
 
     # Integrate f along the Ray
     domainunits = _units(ray(0))
-    return QuadGK.quadgk(t -> f(ray(t)) * domainunits, zero(FP), FP(Inf); rule.kwargs...)[1]
+    integrand(t) = f(ray(t)) * domainunits
+    return QuadGK.quadgk(integrand, zero(FP), FP(Inf); rule.kwargs...)[1]
 end
 
 function integral(
@@ -67,12 +68,12 @@ function integral(
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = FP[0.5]
+    testpoint_parametriccoord = zero(FP)
     integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, FP[0], FP[1]; rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, zero(FP), one(FP); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits

src/specializations/Tetrahedron.jl

@@ -24,11 +24,11 @@ function integral(
         rule::GaussKronrod,
         FP::Type{T} = Float64
 ) where {F <: Function, T <: AbstractFloat}
-    function inner∫₂(v, w)
-        QuadGK.quadgk(u -> f(tetrahedron(u, v, w)), FP(0), FP(1 - v - w); rule.kwargs...)[1]
-    end
-    inner∫₁(w) = QuadGK.quadgk(v -> inner∫₂(v, w), FP(0), FP(1 - w); rule.kwargs...)[1]
-    outer∫ = QuadGK.quadgk(w -> inner∫₁(w), FP(0), FP(1); rule.kwargs...)[1]
+    nil = zero(FP)
+    ∫uvw(u, v, w) = f(tetrahedron(u, v, w))
+    ∫vw(v, w) = QuadGK.quadgk(u -> ∫uvw(u, v, w), nil, FP(1 - v - w); rule.kwargs...)[1]
+    ∫w(w) = QuadGK.quadgk(v -> ∫vw(v, w), nil, FP(1 - w); rule.kwargs...)[1]
+    outer∫ = QuadGK.quadgk(inner∫₁, nil, one(FP); rule.kwargs...)[1]
 
     # Apply barycentric domain correction (volume: 1/6 → actual)
     return 6 * volume(tetrahedron) * outer∫