diff --git a/src/archetypes/spatial_dist.h b/src/archetypes/spatial_dist.h
index bf8abe9c..e0072526 100644
--- a/src/archetypes/spatial_dist.h
+++ b/src/archetypes/spatial_dist.h
@@ -43,8 +43,11 @@ namespace arch {
   struct Uniform : public SpatialDistribution<S, M> {
     Uniform(const M& metric) : SpatialDistribution<S, M> { metric } {}
 
-    Inline auto operator()(const coord_t<M::Dim>&) const -> real_t {
-      return ONE;
+    Inline void operator()(const coord_t<M::Dim>&,
+                           real_t& nppc_distribution,
+                           real_t& weight_distribution) const {
+      nppc_distribution = ONE;
+      weight_distribution = ONE;
     }
   };
 
@@ -68,7 +71,9 @@ namespace arch {
       , target_density { target_density }
       , target_max_density { target_max_density } {}
 
-    Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t {
+    Inline void operator()(const coord_t<M::Dim>& x_Ph,
+                           real_t& nppc_distribution,
+                           real_t& weight_distribution) const {
       coord_t<M::Dim> x_Cd { ZERO };
       metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, x_Cd);
       real_t dens { ZERO };
@@ -88,9 +93,11 @@ namespace arch {
       }
       const auto target = target_density(x_Ph);
       if (0.9 * target > dens) {
-        return (target - dens) / target_max_density;
+        nppc_distribution = (target - dens) / target_max_density;
+        weight_distribution = ONE;
       } else {
-        return ZERO;
+        nppc_distribution = ZERO;
+        weight_distribution = ONE;
       }
     }
   };
@@ -112,7 +119,9 @@ namespace arch {
       , idx { idx }
       , target_density { target_density } {}
 
-    Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t {
+    Inline void operator()(const coord_t<M::Dim>& x_Ph,
+                           real_t& nppc_distribution,
+                           real_t& weight_distribution) const {
       coord_t<M::Dim> x_Cd { ZERO };
       metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, x_Cd);
       real_t dens { ZERO };
@@ -131,9 +140,11 @@ namespace arch {
         raise::KernelError(HERE, "Invalid dimension");
       }
       if (0.9 * target_density > dens) {
-        return (target_density - dens) / target_density;
+        nppc_distribution = (target_density - dens) / target_density;
+        weight_distribution = ONE;
       } else {
-        return ZERO;
+        nppc_distribution = ZERO;
+        weight_distribution = ONE;
       }
     }
   };
diff --git a/src/archetypes/traits.h b/src/archetypes/traits.h
index a6a4ed15..f48d3e53 100644
--- a/src/archetypes/traits.h
+++ b/src/archetypes/traits.h
@@ -47,8 +47,9 @@ namespace arch {
     namespace spatialdist {
 
       template <class SD>
-      concept IsValid = requires(const SD& sdist, const coord_t<SD::D>& x_Ph) {
-        { sdist(x_Ph) } -> std::convertible_to<real_t>;
+      concept IsValid = requires(const SD& sdist, const coord_t<SD::D>& x_Ph,
+                                 real_t& ppc_dist, real_t& weight_dist) {
+        { sdist(x_Ph, ppc_dist, weight_dist) } -> std::same_as<void>;
       };
 
     } // namespace spatialdist
diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp
index 2f9e6a3e..9ed7c81c 100644
--- a/src/kernels/injectors.hpp
+++ b/src/kernels/injectors.hpp
@@ -625,8 +625,11 @@ namespace kernel {
       return idx_h();
     }
 
-    Inline auto injected_ppc(const coord_t<M::Dim>& x_Ph) const -> npart_t {
-      const auto ppc_real = ppc0 * spatial_dist(x_Ph);
+    Inline auto injected_ppc(const coord_t<M::Dim>& x_Ph,
+                             real_t& ppc_dist,
+                             real_t& weight_dist) const -> npart_t {
+      spatial_dist(x_Ph, ppc_dist, weight_dist);
+      const auto ppc_real = ppc0 * ppc_dist;
       auto       ppc      = static_cast<npart_t>(ppc_real);
       auto       rand_gen = random_pool.get_state();
       if (Random<real_t>(rand_gen) < (ppc_real - static_cast<real_t>(ppc))) {
@@ -692,15 +695,15 @@ namespace kernel {
         coord_t<Dim::_1D> x_Cd { i1_ + HALF };
         coord_t<Dim::_1D> x_Ph { ZERO };
         metric.template convert<Crd::Cd, Crd::Ph>(x_Cd, x_Ph);
-
-        const auto ppc = injected_ppc(x_Ph);
+        auto ppc_dist { ZERO }, weight_dist { ONE };
+        const auto ppc = injected_ppc(x_Ph, ppc_dist, weight_dist);
         if (ppc == 0) {
           return;
         }
 
-        auto weight = ONE;
+        auto weight = ONE * weight_dist;
         if constexpr (M::CoordType != Coord::Cart) {
-          weight = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0;
+          weight = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0 * weight_dist;
         }
         for (auto p { 0u }; p < ppc; ++p) {
           const auto index = Kokkos::atomic_fetch_add(&idx(), 1);
@@ -742,15 +745,16 @@ namespace kernel {
           x_Cd_[2] = ZERO;
         }
         metric.template convert<Crd::Cd, Crd::Ph>(x_Cd, x_Ph);
-
-        const auto ppc = injected_ppc(x_Ph);
+        
+        auto ppc_dist { ZERO }, weight_dist { ONE };
+        const auto ppc = injected_ppc(x_Ph, ppc_dist, weight_dist);
         if (ppc == 0) {
           return;
         }
 
-        auto weight = ONE;
+        auto weight = ONE * weight_dist;
         if constexpr (M::CoordType != Coord::Cart) {
-          weight = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * inv_V0;
+          weight = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * inv_V0 * weight_dist;
         }
         for (auto p { 0u }; p < ppc; ++p) {
           const auto index = Kokkos::atomic_fetch_add(&idx(), 1);
@@ -806,16 +810,16 @@ namespace kernel {
         coord_t<Dim::_3D> x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF };
         coord_t<Dim::_3D> x_Ph { ZERO };
         metric.template convert<Crd::Cd, Crd::Ph>(x_Cd, x_Ph);
-
-        const auto ppc = injected_ppc(x_Ph);
+        auto ppc_dist { ZERO }, weight_dist { ONE };
+        const auto ppc = injected_ppc(x_Ph, ppc_dist, weight_dist);
         if (ppc == 0) {
           return;
         }
 
-        auto weight = ONE;
+        auto weight = ONE * weight_dist;
         if constexpr (M::CoordType != Coord::Cart) {
           weight = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ + HALF }) *
-                   inv_V0;
+                   inv_V0 * weight_dist;
         }
         for (auto p { 0u }; p < ppc; ++p) {
           const auto index = Kokkos::atomic_fetch_add(&idx(), 1);