diff --git a/SKIRT/core/AdaptiveMeshSnapshot.cpp b/SKIRT/core/AdaptiveMeshSnapshot.cpp
index 324a74b5..74b31869 100644
--- a/SKIRT/core/AdaptiveMeshSnapshot.cpp
+++ b/SKIRT/core/AdaptiveMeshSnapshot.cpp
@@ -11,9 +11,7 @@
 #include "PathSegmentGenerator.hpp"
 #include "Random.hpp"
 #include "SpatialGridPath.hpp"
-#include "StringUtils.hpp"
 #include "TextInFile.hpp"
-#include "Units.hpp"
 
 ////////////////////////////////////////////////////////////////////
 
@@ -373,13 +371,6 @@ const Array& AdaptiveMeshSnapshot::properties(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-int AdaptiveMeshSnapshot::nearestEntity(Position bfr) const
-{
-    return cellIndex(bfr);
-}
-
-////////////////////////////////////////////////////////////////////
-
 void AdaptiveMeshSnapshot::getEntities(EntityCollection& entities, Position bfr) const
 {
     entities.addSingle(cellIndex(bfr));
diff --git a/SKIRT/core/AdaptiveMeshSnapshot.hpp b/SKIRT/core/AdaptiveMeshSnapshot.hpp
index fe637b60..44ed35ae 100644
--- a/SKIRT/core/AdaptiveMeshSnapshot.hpp
+++ b/SKIRT/core/AdaptiveMeshSnapshot.hpp
@@ -228,10 +228,6 @@ class AdaptiveMeshSnapshot : public Snapshot
         range, the behavior is undefined. */
     const Array& properties(int m) const override;
 
-    /** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ of the cell containing
-        the specified point \f${\bf{r}}\f$, or -1 if the point is outside the domain. */
-    int nearestEntity(Position bfr) const override;
-
 public:
     /** This function sets the specified entity collection to the cell containing the specified
         point \f${\bf{r}}\f$, or to the empty collection if the point is outside the domain. */
diff --git a/SKIRT/core/CellSnapshot.cpp b/SKIRT/core/CellSnapshot.cpp
index 87614ad8..67b319dd 100644
--- a/SKIRT/core/CellSnapshot.cpp
+++ b/SKIRT/core/CellSnapshot.cpp
@@ -10,7 +10,6 @@
 #include "Random.hpp"
 #include "StringUtils.hpp"
 #include "TextInFile.hpp"
-#include "Units.hpp"
 
 ////////////////////////////////////////////////////////////////////
 
@@ -404,13 +403,6 @@ const Array& CellSnapshot::properties(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-int CellSnapshot::nearestEntity(Position bfr) const
-{
-    return _grid ? _grid->cellIndexFor(bfr) : -1;
-}
-
-////////////////////////////////////////////////////////////////////
-
 void CellSnapshot::getEntities(EntityCollection& entities, Position bfr) const
 {
     entities.addSingle(_grid->cellIndexFor(bfr));
diff --git a/SKIRT/core/CellSnapshot.hpp b/SKIRT/core/CellSnapshot.hpp
index 23401082..304f00b7 100644
--- a/SKIRT/core/CellSnapshot.hpp
+++ b/SKIRT/core/CellSnapshot.hpp
@@ -109,11 +109,6 @@ class CellSnapshot : public Snapshot
         range, the behavior is undefined. */
     const Array& properties(int m) const override;
 
-    /** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ of the cell containing
-        the specified point \f${\bf{r}}\f$, or -1 if the point is outside the domain, if there are
-        no cells in the snapshot, or if the search data structures were not created. */
-    int nearestEntity(Position bfr) const override;
-
 public:
     /** This function sets the specified entity collection to the cell containing the specified
         point \f${\bf{r}}\f$, or to the empty collection if the point is outside the domain or if
diff --git a/SKIRT/core/ParticleSnapshot.cpp b/SKIRT/core/ParticleSnapshot.cpp
index b691bbc9..f2ee7108 100644
--- a/SKIRT/core/ParticleSnapshot.cpp
+++ b/SKIRT/core/ParticleSnapshot.cpp
@@ -11,7 +11,6 @@
 #include "SmoothingKernel.hpp"
 #include "StringUtils.hpp"
 #include "TextInFile.hpp"
-#include "Units.hpp"
 
 ////////////////////////////////////////////////////////////////////
 
@@ -212,24 +211,6 @@ class ParticleSnapshot::ParticleGrid : public Box
         return _listv[index(_m, i, j, k)];
     }
 
-    /** This function returns a pointer to the particle centered nearest to the specified position,
-        or the null pointer if the specified position is outside of the grid. */
-    const Particle* nearestParticle(Vec r) const
-    {
-        const Particle* nearestParticle = nullptr;
-        double nearestSquaredDistance = std::numeric_limits<double>::infinity();
-        for (const Particle* particle : particlesFor(r))
-        {
-            double d2 = (r - particle->center()).norm2();
-            if (d2 < nearestSquaredDistance)
-            {
-                nearestParticle = particle;
-                nearestSquaredDistance = d2;
-            }
-        }
-        return nearestParticle;
-    }
-
     /** This function replaces the contents of the specified entity collection by the set of
         particles with a smoothing kernel that overlaps the specified point \f${\bf{r}}\f$. The
         weight corresponding to each particle is set to the particle's smoothing kernel value at
@@ -555,14 +536,6 @@ const Array& ParticleSnapshot::properties(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-int ParticleSnapshot::nearestEntity(Position bfr) const
-{
-    const Particle* nearestParticle = _grid ? _grid->nearestParticle(bfr) : nullptr;
-    return nearestParticle ? nearestParticle->index() : -1;
-}
-
-////////////////////////////////////////////////////////////////////
-
 void ParticleSnapshot::getEntities(EntityCollection& entities, Position bfr) const
 {
     _grid->getEntities(entities, bfr, _kernel);
diff --git a/SKIRT/core/ParticleSnapshot.hpp b/SKIRT/core/ParticleSnapshot.hpp
index 867ab7e3..56c5b4d6 100644
--- a/SKIRT/core/ParticleSnapshot.hpp
+++ b/SKIRT/core/ParticleSnapshot.hpp
@@ -108,12 +108,6 @@ class ParticleSnapshot : public Snapshot
         of range, the behavior is undefined. */
     const Array& properties(int m) const override;
 
-    /** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ of the particle whose
-        center is nearest to the specified point \f${\bf{r}}\f$, or -1 if the point is outside the
-        domain, if there are no particles in the snapshot, or if the search data structures were
-        not created. */
-    int nearestEntity(Position bfr) const override;
-
 public:
     /** This function replaces the contents of the specified entity collection by the set of
         particles with a smoothing kernel that overlaps the specified point \f${\bf{r}}\f$. The
diff --git a/SKIRT/core/Snapshot.cpp b/SKIRT/core/Snapshot.cpp
index 23d4f0c8..7ed05fa7 100644
--- a/SKIRT/core/Snapshot.cpp
+++ b/SKIRT/core/Snapshot.cpp
@@ -4,6 +4,7 @@
 ///////////////////////////////////////////////////////////////// */
 
 #include "Snapshot.hpp"
+#include "EntityCollection.hpp"
 #include "Log.hpp"
 #include "Random.hpp"
 #include "StringUtils.hpp"
@@ -242,14 +243,6 @@ double Snapshot::initialMass(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-double Snapshot::initialMass(Position bfr) const
-{
-    int m = nearestEntity(bfr);
-    return m >= 0 ? initialMass(m) : 0.;
-}
-
-////////////////////////////////////////////////////////////////////
-
 double Snapshot::currentMass(int m) const
 {
     return properties(m)[currentMassIndex()];
@@ -257,14 +250,6 @@ double Snapshot::currentMass(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-double Snapshot::currentMass(Position bfr) const
-{
-    int m = nearestEntity(bfr);
-    return m >= 0 ? currentMass(m) : 0.;
-}
-
-////////////////////////////////////////////////////////////////////
-
 double Snapshot::metallicity(int m) const
 {
     return properties(m)[metallicityIndex()];
@@ -274,8 +259,9 @@ double Snapshot::metallicity(int m) const
 
 double Snapshot::metallicity(Position bfr) const
 {
-    int m = nearestEntity(bfr);
-    return m >= 0 ? metallicity(m) : 0.;
+    thread_local EntityCollection entities;  // can be reused for all queries in a given execution thread
+    getEntities(entities, bfr);
+    return entities.averageValue([this](int m) { return metallicity(m); }, [this](int m) { return currentMass(m); });
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -287,14 +273,6 @@ double Snapshot::age(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-double Snapshot::age(Position bfr) const
-{
-    int m = nearestEntity(bfr);
-    return m >= 0 ? age(m) : 0.;
-}
-
-////////////////////////////////////////////////////////////////////
-
 double Snapshot::temperature(int m) const
 {
     return properties(m)[temperatureIndex()];
@@ -304,8 +282,9 @@ double Snapshot::temperature(int m) const
 
 double Snapshot::temperature(Position bfr) const
 {
-    int m = nearestEntity(bfr);
-    return m >= 0 ? temperature(m) : 0.;
+    thread_local EntityCollection entities;  // can be reused for all queries in a given execution thread
+    getEntities(entities, bfr);
+    return entities.averageValue([this](int m) { return temperature(m); }, [this](int m) { return currentMass(m); });
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -320,8 +299,9 @@ Vec Snapshot::velocity(int m) const
 
 Vec Snapshot::velocity(Position bfr) const
 {
-    int m = nearestEntity(bfr);
-    return m >= 0 ? velocity(m) : Vec();
+    thread_local EntityCollection entities;  // can be reused for all queries in a given execution thread
+    getEntities(entities, bfr);
+    return entities.averageValue([this](int m) { return velocity(m); }, [this](int m) { return currentMass(m); });
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -333,14 +313,6 @@ double Snapshot::velocityDispersion(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-double Snapshot::velocityDispersion(Position bfr) const
-{
-    int m = nearestEntity(bfr);
-    return m >= 0 ? velocityDispersion(m) : 0.;
-}
-
-////////////////////////////////////////////////////////////////////
-
 Vec Snapshot::magneticField(int m) const
 {
     const auto& propv = properties(m);
@@ -351,8 +323,9 @@ Vec Snapshot::magneticField(int m) const
 
 Vec Snapshot::magneticField(Position bfr) const
 {
-    int m = nearestEntity(bfr);
-    return m >= 0 ? magneticField(m) : Vec();
+    thread_local EntityCollection entities;  // can be reused for all queries in a given execution thread
+    getEntities(entities, bfr);
+    return entities.averageValue([this](int m) { return magneticField(m); }, [this](int m) { return currentMass(m); });
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -369,7 +342,22 @@ void Snapshot::parameters(int m, Array& params) const
 
 void Snapshot::parameters(Position bfr, Array& params) const
 {
-    int m = nearestEntity(bfr);
+    thread_local EntityCollection entities;  // can be reused for all queries in a given execution thread
+    getEntities(entities, bfr);
+
+    // look for the entity with the highest weight
+    double wmax = 0.;
+    int m = -1;
+    for (const auto& entity : entities)
+    {
+        if (entity.second > wmax)
+        {
+            wmax = entity.second;
+            m = entity.first;
+        }
+    }
+
+    // if found, use it
     if (m >= 0)
         parameters(m, params);
     else
diff --git a/SKIRT/core/Snapshot.hpp b/SKIRT/core/Snapshot.hpp
index 946be3f9..edb4a14c 100644
--- a/SKIRT/core/Snapshot.hpp
+++ b/SKIRT/core/Snapshot.hpp
@@ -359,13 +359,6 @@ class Snapshot
         range, the behavior is undefined. */
     virtual const Array& properties(int m) const = 0;
 
-    /** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ for the entity at or
-        nearest to the specified point \f${\bf{r}}\f$, or -1 if the point is outside the domain or
-        if there are no entities in the snapshot. For a cell-based snapshot, the function returns
-        the index of the cell containing the given point. For a particle-based snapshot, the
-        function returns the index of the particle whose center is nearest to the given point. */
-    virtual int nearestEntity(Position bfr) const = 0;
-
 public:
     /** This function replaces the contents of the specified entity collection by the set of
         entities that overlap the specified point \f${\bf{r}}\f$, with their corresponding weights.
@@ -407,11 +400,6 @@ class Snapshot
         range, the behavior is undefined. */
     double initialMass(int m) const;
 
-    /** This function returns the initial mass of the entity nearest to (or at) the specified point
-        \f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
-        initial mass is not being imported, the behavior is undefined. */
-    double initialMass(Position bfr) const;
-
     /** This function returns true if the current mass is being imported, and false otherwise. */
     bool hasCurrentMass() const { return _currentMassIndex >= 0; }
 
@@ -420,11 +408,6 @@ class Snapshot
         range, the behavior is undefined. */
     double currentMass(int m) const;
 
-    /** This function returns the current mass of the entity nearest to (or at) the specified point
-        \f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
-        current mass is not being imported, the behavior is undefined. */
-    double currentMass(Position bfr) const;
-
     /** This function returns true if the metallicity is being imported, and false otherwise. */
     bool hasMetallicity() const { return _metallicityIndex >= 0; }
 
@@ -433,7 +416,7 @@ class Snapshot
         range, the behavior is undefined. */
     double metallicity(int m) const;
 
-    /** This function returns the metallicity of the entity nearest to (or at) the specified point
+    /** This function returns the metallicity at the specified point
         \f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
         metallicity is not being imported, the behavior is undefined. */
     double metallicity(Position bfr) const;
@@ -446,11 +429,6 @@ class Snapshot
         */
     double age(int m) const;
 
-    /** This function returns the age of the entity nearest to (or at) the specified point
-        \f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the age
-        is not being imported, the behavior is undefined. */
-    double age(Position bfr) const;
-
     /** This function returns true if the temperature is being imported, and false otherwise. */
     bool hasTemperature() const { return _temperatureIndex >= 0; }
 
@@ -459,7 +437,7 @@ class Snapshot
         range, the behavior is undefined. */
     double temperature(int m) const;
 
-    /** This function returns the temperature of the entity nearest to (or at) the specified point
+    /** This function returns the temperature at the specified point
         \f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
         temperature is not being imported, the behavior is undefined. */
     double temperature(Position bfr) const;
@@ -472,7 +450,7 @@ class Snapshot
         the behavior is undefined. */
     Vec velocity(int m) const;
 
-    /** This function returns the velocity of the entity nearest to (or at) the specified point
+    /** This function returns the velocity at the specified point
         \f${\bf{r}}\f$. If the point is outside the domain, the function returns zero velocity. If
         the velocity is not being imported, the behavior is undefined. */
     Vec velocity(Position bfr) const;
@@ -486,12 +464,6 @@ class Snapshot
         of range, the behavior is undefined. */
     double velocityDispersion(int m) const;
 
-    /** This function returns the velocity dispersion of the entity nearest to (or at) the
-        specified point \f${\bf{r}}\f$. If the point is outside the domain, the function returns
-        zero dispersion. If the velocity dispersion is not being imported, the behavior is
-        undefined. */
-    double velocityDispersion(Position bfr) const;
-
     /** This function returns true if the magnetic field is being imported, and false otherwise. */
     bool hasMagneticField() const { return _magneticFieldIndex >= 0; }
 
@@ -500,7 +472,7 @@ class Snapshot
         range, the behavior is undefined. */
     Vec magneticField(int m) const;
 
-    /** This function returns the magnetic field vector of the entity nearest to (or at) the
+    /** This function returns the magnetic field vector at the
         specified point \f${\bf{r}}\f$. If the point is outside the domain, the function returns a
         zero magnetic field. If the magnetic field is not being imported, the behavior is
         undefined. */
@@ -515,7 +487,7 @@ class Snapshot
         index is out of range, the behavior is undefined. */
     void parameters(int m, Array& params) const;
 
-    /** This function stores the parameters of the entity nearest to (or at) the specified point
+    /** This function stores the parameters at the specified point
         \f${\bf{r}}\f$ into the given array. If the point is outside the domain, the function
         returns the appropriate number of zero parameter values. If parameters are not being
         imported, the behavior is undefined. */
diff --git a/SKIRT/core/VoronoiMeshSnapshot.cpp b/SKIRT/core/VoronoiMeshSnapshot.cpp
index ffd925d6..3f3b683b 100644
--- a/SKIRT/core/VoronoiMeshSnapshot.cpp
+++ b/SKIRT/core/VoronoiMeshSnapshot.cpp
@@ -19,7 +19,6 @@
 #include "StringUtils.hpp"
 #include "Table.hpp"
 #include "TextInFile.hpp"
-#include "Units.hpp"
 #include <set>
 #include "container.hh"
 
@@ -1108,13 +1107,6 @@ const Array& VoronoiMeshSnapshot::properties(int m) const
 
 ////////////////////////////////////////////////////////////////////
 
-int VoronoiMeshSnapshot::nearestEntity(Position bfr) const
-{
-    return _blocktrees.size() ? cellIndex(bfr) : -1;
-}
-
-////////////////////////////////////////////////////////////////////
-
 void VoronoiMeshSnapshot::getEntities(EntityCollection& entities, Position bfr) const
 {
     entities.addSingle(cellIndex(bfr));
diff --git a/SKIRT/core/VoronoiMeshSnapshot.hpp b/SKIRT/core/VoronoiMeshSnapshot.hpp
index ca3def78..2629ace4 100644
--- a/SKIRT/core/VoronoiMeshSnapshot.hpp
+++ b/SKIRT/core/VoronoiMeshSnapshot.hpp
@@ -346,11 +346,6 @@ class VoronoiMeshSnapshot : public Snapshot
         range, the behavior is undefined. */
     const Array& properties(int m) const override;
 
-    /** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ of the cell containing
-        the specified point \f${\bf{r}}\f$, or -1 if the point is outside the domain, if there
-        are no cells in the snapshot, or if the search data structures were not created. */
-    int nearestEntity(Position bfr) const override;
-
 public:
     /** This function sets the specified entity collection to the cell containing the specified
         point \f${\bf{r}}\f$, or to the empty collection if the point is outside the domain or if
diff --git a/SKIRT/utils/EntityCollection.cpp b/SKIRT/utils/EntityCollection.cpp
index d18acb63..25234ffe 100644
--- a/SKIRT/utils/EntityCollection.cpp
+++ b/SKIRT/utils/EntityCollection.cpp
@@ -63,6 +63,19 @@ std::pair<double, double> EntityCollection::average(std::function<double(int m)>
 
 ////////////////////////////////////////////////////////////////////
 
+double EntityCollection::averageValue(std::function<double(int)> value, std::function<double(int)> weight)
+{
+    auto numEntities = _entities.size();
+    if (numEntities == 0) return 0.;
+    if (numEntities == 1) return value(_entities.cbegin()->first);
+
+    double sumvw, sumw;
+    std::tie(sumvw, sumw) = average(value, weight);
+    return sumw > 0. ? sumvw / sumw : 0.;
+}
+
+////////////////////////////////////////////////////////////////////
+
 std::pair<Vec, double> EntityCollection::average(std::function<Vec(int m)> value, std::function<double(int m)> weight)
 {
     Vec sumvw;
@@ -78,3 +91,17 @@ std::pair<Vec, double> EntityCollection::average(std::function<Vec(int m)> value
 }
 
 ////////////////////////////////////////////////////////////////////
+
+Vec EntityCollection::averageValue(std::function<Vec(int)> value, std::function<double(int)> weight)
+{
+    auto numEntities = _entities.size();
+    if (numEntities == 0) return Vec();
+    if (numEntities == 1) return value(_entities.cbegin()->first);
+
+    Vec sumvw;
+    double sumw;
+    std::tie(sumvw, sumw) = average(value, weight);
+    return sumw > 0. ? sumvw / sumw : Vec();
+}
+
+////////////////////////////////////////////////////////////////////
diff --git a/SKIRT/utils/EntityCollection.hpp b/SKIRT/utils/EntityCollection.hpp
index d74703f6..3888ae49 100644
--- a/SKIRT/utils/EntityCollection.hpp
+++ b/SKIRT/utils/EntityCollection.hpp
@@ -83,6 +83,17 @@ class EntityCollection
         \f$\sum_m \omega(m)\,w_m\f$. */
     std::pair<double, double> average(std::function<double(int m)> value, std::function<double(int m)> weight);
 
+    /** This function returns the weighted average of a given scalar field \f$f(m)\f$ with given
+        external weight \f$\omega(m)\f$ over all entities in the collection. The arguments
+        respectively specify the scalar field \f$f(m)\f$ and the corresponding external weight
+        \f$\omega(m)\f$. These functions should return the field value respectively the external
+        weight corresponding to a given entity index.
+
+        The function combines the external weights with the weights stored internally for each
+        entity in the collection. Specifically, it calculates \f$\sum_m f(m)\,\omega(m)\,w_m /
+        \sum_m \omega(m)\,w_m\f$. */
+    double averageValue(std::function<double(int m)> value, std::function<double(int m)> weight);
+
     /** This function returns the nominator and denominator for the weighted average of a given
         vector field \f${\bf{f}}(m)\f$ with given external weight \f$\omega(m)\f$ over all entities
         in the collection. The arguments respectively specify the vector field \f${\bf{f}}(m)\f$
@@ -94,6 +105,17 @@ class EntityCollection
         {\bf{f}}(m)\,\omega(m)\,w_m\f$ and \f$\sum_m \omega(m)\,w_m\f$. */
     std::pair<Vec, double> average(std::function<Vec(int m)> value, std::function<double(int m)> weight);
 
+    /** This function returns the weighted average of a given vector field \f${\bf{f}}(m)\f$ with
+        given external weight \f$\omega(m)\f$ over all entities in the collection. The arguments
+        respectively specify the vector field \f${\bf{f}}(m)\f$ and the corresponding external
+        weight \f$\omega(m)\f$. These functions should return the field value respectively the
+        external weight corresponding to a given entity index.
+
+        The function combines the external weights with the weights stored internally for each
+        entity in the collection. Specifically, it calculates \f$\sum_m {\bf{f}}(m)\,\omega(m)\,w_m
+        / \sum_m \omega(m)\,w_m\f$. */
+    Vec averageValue(std::function<Vec(int m)> value, std::function<double(int m)> weight);
+
     // ------- Data members -------
 
 private: