remove nearestEntity and use getEntities instead

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));

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. */

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));

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

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);

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

SKIRT/core/Snapshot.cpp

@@ -4,6 +4,7 @@
 ///////////////////////////////////////////////////////////////// */
 
 #include "Snapshot.hpp"
+#include "EntityCollection.hpp"
 #include "Log.hpp"
 #include "Random.hpp"
 #include "StringUtils.hpp"
@@ -242,29 +243,13 @@ 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()];
 }
 
 ////////////////////////////////////////////////////////////////////
 
-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); });
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -335,8 +315,10 @@ double Snapshot::velocityDispersion(int m) const
 
 double Snapshot::velocityDispersion(Position bfr) const
 {
-    int m = nearestEntity(bfr);
-    return m >= 0 ? velocityDispersion(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 velocityDispersion(m); },
+                                 [this](int m) { return currentMass(m); });
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -351,8 +333,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 +352,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

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,7 +464,7 @@ 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
+    /** This function returns the velocity dispersion 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. */
@@ -500,7 +478,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 +493,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. */

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));