Index: /cpp/frams/genetics/fS/fS_conv.cpp
===================================================================
--- /cpp/frams/genetics/fS/fS_conv.cpp	(revision 1005)
+++ /cpp/frams/genetics/fS/fS_conv.cpp	(revision 1006)
@@ -18,8 +18,5 @@
 	}
 
-	Model model;
-	model.open(using_checkpoints);
-	genotype->buildModel(model);
-	model.close();
+	Model model = genotype->buildModel(using_checkpoints);
 	delete genotype;
 
Index: /cpp/frams/genetics/fS/fS_conv.h
===================================================================
--- /cpp/frams/genetics/fS/fS_conv.h	(revision 1005)
+++ /cpp/frams/genetics/fS/fS_conv.h	(revision 1006)
@@ -10,5 +10,5 @@
 
 /**
- * Genotype converter from fS to f0.
+ * Genotype converter from fS to f0s.
  */
 class GenoConv_fS0 : public GenoConverter
Index: /cpp/frams/genetics/fS/fS_general.cpp
===================================================================
--- /cpp/frams/genetics/fS/fS_general.cpp	(revision 1005)
+++ /cpp/frams/genetics/fS/fS_general.cpp	(revision 1006)
@@ -17,19 +17,22 @@
 bool fS_Genotype::TURN_WITH_ROTATION = false;
 
-const std::map<string, double> defaultValues = {
-		{INGESTION,      0.5},//Model::getDefPart().ingest },
-		{FRICTION,       0.5},//Model::getDefPart().friction},
-		{STIFFNESS,	 	 0.5},//Model::getDefJoint().stif},
-		{ROT_X,          0.0},
-		{ROT_Y,          0.0},
-		{ROT_Z,          0.0},
-		{RX,             0.0},
-		{RY,             0.0},
-		{RZ,             0.0},
-		{SIZE,           1.0},
-		{SIZE_X,         0.5},//Model::getDefPart().scale.x},
-		{SIZE_Y,         0.5},//Model::getDefPart().scale.y},
-		{SIZE_Z,         0.5}//Model::getDefPart().scale.z}
-};
+void Node::prepareParams()
+{
+	defaultValues = {
+			{INGESTION, Model::getDefPart().ingest},
+			{FRICTION,  Model::getDefPart().friction},
+			{STIFFNESS, Model::getDefJoint().stif},
+			{ROT_X,     0.0},
+			{ROT_Y,     0.0},
+			{ROT_Z,     0.0},
+			{RX,        0.0},
+			{RY,        0.0},
+			{RZ,        0.0},
+			{SIZE,      1.0},
+			{SIZE_X,    Model::getDefPart().scale.x},
+			{SIZE_Y,    Model::getDefPart().scale.y},
+			{SIZE_Z,    Model::getDefPart().scale.z}
+	};
+}
 
 double fS_stod(const string&  str, int start, size_t* size)
@@ -76,6 +79,8 @@
 
 
-fS_Neuron::fS_Neuron(const char *str, int start, int length)
-{
+fS_Neuron::fS_Neuron(const char *str, int _start, int length)
+{
+	start = _start + 1;
+	end = start + length;
 	if (length == 0)
 		return;
@@ -123,4 +128,5 @@
 Node::Node(Substring &restOfGeno, Node *_parent, GenotypeParams _genotypeParams)
 {
+	prepareParams();
 	partDescription = new Substring(restOfGeno);
 	genotypeParams = _genotypeParams;
@@ -165,5 +171,5 @@
 	for (int i = 0; i < restOfGenotype.len; i++)
 	{
-		if (GENE_TO_SHAPETYPE.find(restOfGenotype.at(i)) != GENE_TO_SHAPETYPE.end())
+		if (GENE_TO_SHAPE.find(restOfGenotype.at(i)) != GENE_TO_SHAPE.end())
 			return i;
 	}
@@ -193,6 +199,6 @@
 void Node::extractPartType(Substring &restOfGenotype)
 {
-	auto itr = GENE_TO_SHAPETYPE.find(restOfGenotype.at(0));
-	if (itr == GENE_TO_SHAPETYPE.end())
+	auto itr = GENE_TO_SHAPE.find(restOfGenotype.at(0));
+	if (itr == GENE_TO_SHAPE.end())
 		throw fS_Exception("Invalid part type", restOfGenotype.start);
 
@@ -465,4 +471,5 @@
 		Neuro *neuro = new Neuro(*neurons[i]);
 		model.addNeuro(neuro);
+		neuro->addMapping(MultiRange(IRange(neurons[i]->start, neurons[i]->end)));
 		if (neuro->getClass()->preflocation == NeuroClass::PREFER_JOINT && parent != nullptr)
 		{
@@ -536,5 +543,5 @@
 		result += std::string(count, mod).c_str();
 	}
-	result += SHAPETYPE_TO_GENE.at(partType);
+	result += SHAPE_TO_GENE.at(partType);
 
 	if (!neurons.empty())
@@ -649,9 +656,16 @@
 }
 
-void fS_Genotype::buildModel(Model &model)
-{
+Model fS_Genotype::buildModel(bool using_checkpoints)
+{
+
+	Model model;
+	model.open(using_checkpoints);
+
 	getState();
 	startNode->buildModel(model, nullptr);
 	buildNeuroConnections(model);
+
+	model.close();
+	return model;
 }
 
Index: /cpp/frams/genetics/fS/fS_general.h
===================================================================
--- /cpp/frams/genetics/fS/fS_general.h	(revision 1005)
+++ /cpp/frams/genetics/fS/fS_general.h	(revision 1006)
@@ -27,5 +27,5 @@
 const char NEURON_END = ']';
 const char NEURON_SEPARATOR = ';';
-const SString NEURON_INTERNAL_SEPARATOR("_");
+const SString NEURON_INTERNAL_SEPARATOR("'");
 #define NEURON_I_W_SEPARATOR ':'
 //@}
@@ -69,5 +69,5 @@
 const char CUBOID = 'C';
 const char CYLINDER = 'R';
-const std::unordered_map<Part::Shape, char> SHAPETYPE_TO_GENE = {
+const std::unordered_map<Part::Shape, char> SHAPE_TO_GENE = {
 		{Part::Shape::SHAPE_ELLIPSOID, ELLIPSOID},
 		{Part::Shape::SHAPE_CUBOID,    CUBOID},
@@ -76,10 +76,10 @@
 
 // This map is inverse to SHAPE_TO_SYMBOL. Those two should be compatible
-const std::unordered_map<char, Part::Shape> GENE_TO_SHAPETYPE = {
+const std::unordered_map<char, Part::Shape> GENE_TO_SHAPE = {
 		{ELLIPSOID, Part::Shape::SHAPE_ELLIPSOID},
 		{CUBOID,    Part::Shape::SHAPE_CUBOID},
 		{CYLINDER,  Part::Shape::SHAPE_CYLINDER},
 };
-const int SHAPE_COUNT = 3;    // This should be the count of SHAPETYPE_TO_GENE and GENE_TO_SHAPETYPE
+const int SHAPE_COUNT = 3;    // This should be the count of SHAPE_TO_GENE and GENE_TO_SHAPE
 
 const char DEFAULT_JOINT = 'a';
@@ -100,6 +100,4 @@
 };
 
-extern const std::map<string, double> defaultValues;
-
 /** @name Number of tries of performing a mutation before GENOPER_FAIL is returned */
 #define mutationTries  20
@@ -198,7 +196,6 @@
 	MultiRange toMultiRange()
 	{
-		MultiRange range;
-		range.add(start, start + len - 1);
-		return range;
+		int end = start + len - 1;
+		return MultiRange(IRange(start, end));
 	}
 };
@@ -244,4 +241,5 @@
 {
 public:
+	int start, end;
 	std::map<int, double> inputs;
 
@@ -274,10 +272,12 @@
 	Part *part;     /// A part object built from node. Used in building the Model
 	int partCodeLen; /// The length of substring that directly describes the corresponding part
+	std::map<string, double> defaultValues;
 	GenotypeParams genotypeParams;
-
 
 	vector<Node *> children;    /// Vector of all direct children
 	std::map<char, int> modifiers;     /// Vector of all modifiers
 	vector<fS_Neuron *> neurons;    /// Vector of all the neurons
+
+	void prepareParams();
 
 	double getDistance();
@@ -398,14 +398,4 @@
 
 	/**
-	 * Change the value of the size parameter by given multiplier
-	 * Do not change the value if any of the size restrictions is not satisfied
-	 * @param paramKey
-	 * @param multiplier
-	 * @param ensureCircleSection
-	 * @return True if the parameter value was change, false otherwise
-	 */
-	bool mutateSizeParam(string paramKey,  bool ensureCircleSection);
-
-	/**
 	 * Counts all the nodes in subtree
 	 * @return node count
@@ -521,5 +511,5 @@
 	 * @param a reference to a model that will contain a built model
 	 */
-	void buildModel(Model &model);
+	Model buildModel(bool using_checkpoints);
 
 	/**
Index: /cpp/frams/genetics/fS/fS_oper.cpp
===================================================================
--- /cpp/frams/genetics/fS/fS_oper.cpp	(revision 1005)
+++ /cpp/frams/genetics/fS/fS_oper.cpp	(revision 1006)
@@ -36,40 +36,42 @@
 
 
-
-const std::map<string, double> minValues = {
-		{INGESTION,      0},//Model::getMinPart().ingest},
-		{FRICTION,       0},//Model::getMinPart().friction},
-		{STIFFNESS,	 	 0.1},
-		{ROT_X,          -M_PI},
-		{ROT_Y,          -M_PI},
-		{ROT_Z,          -M_PI},
-		{RX,             -M_PI},
-		{RY,             -M_PI},
-		{RZ,             -M_PI},
-		{SIZE,           0.01},
-		{SIZE_X,         0},//Model::getMinPart().scale.x},
-		{SIZE_Y,         0},//Model::getMinPart().scale.y},
-		{SIZE_Z,         0}//Model::getMinPart().scale.z}
-};
-
-const std::map<string, double> maxValues = {
-		{INGESTION,      1},//Model::getMaxPart().ingest},
-		{FRICTION,       1},//Model::getMaxPart().friction},
-		{STIFFNESS,	 	 0.5},
-		{ROT_X,          M_PI},
-		{ROT_Y,          M_PI},
-		{ROT_Z,          M_PI},
-		{RX,             M_PI},
-		{RY,             M_PI},
-		{RZ,             M_PI},
-		{SIZE,           100.0},
-		{SIZE_X,         1},//Model::getMaxPart().scale.x},
-		{SIZE_Y,         1},//Model::getMaxPart().scale.y},
-		{SIZE_Z,         1}//Model::getMaxPart().scale.z}
-};
-
+void GenoOper_fS::prepareParams()
+{
+	minValues = {
+			{INGESTION, Model::getMinPart().ingest},
+			{FRICTION,  Model::getMinPart().friction},
+			{STIFFNESS, 0.1},
+			{ROT_X,     -M_PI},
+			{ROT_Y,     -M_PI},
+			{ROT_Z,     -M_PI},
+			{RX,        -M_PI},
+			{RY,        -M_PI},
+			{RZ,        -M_PI},
+			{SIZE,      0.01},
+			{SIZE_X,    Model::getMinPart().scale.x},
+			{SIZE_Y,    Model::getMinPart().scale.y},
+			{SIZE_Z,    Model::getMinPart().scale.z}
+	};
+
+	maxValues = {
+			{INGESTION, Model::getMaxPart().ingest},
+			{FRICTION,  Model::getMaxPart().friction},
+			{STIFFNESS, 0.5},
+			{ROT_X,     M_PI},
+			{ROT_Y,     M_PI},
+			{ROT_Z,     M_PI},
+			{RX,        M_PI},
+			{RY,        M_PI},
+			{RZ,        M_PI},
+			{SIZE,      100.0},
+			{SIZE_X,    Model::getMaxPart().scale.x},
+			{SIZE_Y,    Model::getMaxPart().scale.y},
+			{SIZE_Z,    Model::getMaxPart().scale.z}
+	};
+}
 
 GenoOper_fS::GenoOper_fS()
 {
+	prepareParams();
 	par.setParamTab(genooper_fS_paramtab);
 	par.select(this);
@@ -332,5 +334,5 @@
 	geno.getState();
 	Node *node = geno.chooseNode();
-	char partType = SHAPETYPE_TO_GENE.at(availablePartShapes[rndUint(availablePartShapes.size())]);
+	char partType = SHAPE_TO_GENE.at(availablePartShapes[rndUint(availablePartShapes.size())]);
 
 	Substring substring(&partType, 0, 1);
@@ -374,7 +376,7 @@
 	{
 		geno.getState();
-		newNode->mutateSizeParam(SIZE_X, true);
-		newNode->mutateSizeParam(SIZE_Y, true);
-		newNode->mutateSizeParam(SIZE_Z, true);
+		mutateSizeParam(newNode, SIZE_X, true);
+		mutateSizeParam(newNode, SIZE_Y, true);
+		mutateSizeParam(newNode, SIZE_Z, true);
 	}
 	return true;
@@ -483,4 +485,5 @@
 	bool isRadius = isRadiusOfBase || key == SIZE_X;
 	if (ensureCircleSection && isRadius)
+	if (ensureCircleSection && isRadius)
 	{
 		if (randomNode->partType == Part::Shape::SHAPE_ELLIPSOID)
@@ -490,5 +493,5 @@
 	}
 	// Add modified default value for param
-	randomNode->params[key] = mutateCreep('f', defaultValues.at(key), minValues.at(key), maxValues.at(key), true);
+	randomNode->params[key] = mutateCreep('f', randomNode->defaultValues.at(key), minValues.at(key), maxValues.at(key), true);
 	return true;
 }
@@ -530,5 +533,5 @@
 				return true;
 			} else
-				return randomNode->mutateSizeParam(it->first, ensureCircleSection);
+				return mutateSizeParam(randomNode, it->first, ensureCircleSection);
 		}
 	}
@@ -696,8 +699,8 @@
 }
 
-bool Node::mutateSizeParam(string key, bool ensureCircleSection)
-{
-	double oldValue = getParam(key);
-	double volume = calculateVolume();
+bool GenoOper_fS::mutateSizeParam(Node *node, string key, bool ensureCircleSection)
+{
+	double oldValue = node->getParam(key);
+	double volume = node->calculateVolume();
 	double valueAtMinVolume, valueAtMaxVolume;
 	if(key == SIZE)
@@ -715,12 +718,12 @@
 	double max = std::min(maxValues.at(key), valueAtMaxVolume);
 
-	params[key] = GenoOperators::mutateCreep('f', getParam(key), min, max, true);
-
-	if (!ensureCircleSection || isPartSizeValid())
+	node->params[key] = GenoOperators::mutateCreep('f', node->getParam(key), min, max, true);
+
+	if (!ensureCircleSection || node->isPartSizeValid())
 		return true;
 	else
 	{
-		params[key] = oldValue;
-		return false;
-	}
-}
+		node->params[key] = oldValue;
+		return false;
+	}
+}
Index: /cpp/frams/genetics/fS/fS_oper.h
===================================================================
--- /cpp/frams/genetics/fS/fS_oper.h	(revision 1005)
+++ /cpp/frams/genetics/fS/fS_oper.h	(revision 1006)
@@ -32,4 +32,6 @@
 const int PARENT_COUNT = 2;
 
+
+
 class GenoOper_fS : public GenoOperators
 {
@@ -40,4 +42,7 @@
 	paInt useElli, useCub,  useCyl;
 	paInt strongAddPart;
+
+	std::map<string, double> minValues;
+	std::map<string, double> maxValues;
 
 	GenoOper_fS();
@@ -129,4 +134,16 @@
 
 	bool changeNeuroParam(fS_Genotype &geno);
+
+	/**
+	 * Change the value of the size parameter by given multiplier
+	 * Do not change the value if any of the size restrictions is not satisfied
+	 * @param paramKey
+	 * @param multiplier
+	 * @param ensureCircleSection
+	 * @return True if the parameter value was change, false otherwise
+	 */
+	bool mutateSizeParam(Node *node, string key, bool ensureCircleSection);
+
+	void prepareParams();
 };
 
Index: /cpp/frams/genetics/fS/part_distance_estimator.h
===================================================================
--- /cpp/frams/genetics/fS/part_distance_estimator.h	(revision 1006)
+++ /cpp/frams/genetics/fS/part_distance_estimator.h	(revision 1006)
@@ -0,0 +1,235 @@
+// This file is a part of Framsticks SDK.  http://www.framsticks.com/
+// Copyright (C) 2019-2020  Maciej Komosinski and Szymon Ulatowski.
+// See LICENSE.txt for details.
+
+#ifndef _PART_DISTANCE_ESTIMATOR_H_
+#define _PART_DISTANCE_ESTIMATOR_H_
+
+class fS_Utils
+{
+public:
+	static void rotateVector(Pt3D &vector, const Pt3D &rotation)
+	{
+		Orient rotmatrix = Orient_1;
+		rotmatrix.rotate(rotation);
+		vector = rotmatrix.transform(vector);
+	}
+
+	static double avg(double a, double b)
+	{
+		return 0.5 * (a + b);
+	}
+
+	static double min3(const Pt3D &p)
+	{
+		double tmp = p.x;
+		if (p.y < tmp)
+			tmp = p.y;
+		if (p.z < tmp)
+			tmp = p.z;
+		return tmp;
+	}
+
+	static double max3(const Pt3D &p)
+	{
+		double tmp = p.x;
+		if (p.y > tmp)
+			tmp = p.y;
+		if (p.z > tmp)
+			tmp = p.z;
+		return tmp;
+	}
+};
+
+class PartDistanceEstimator
+{
+	/**
+	 * Used in finding the proper distance between the parts
+	 * distance between spheres / sphere radius
+	 * That default value can be changed in certain cases
+	 * */
+	static constexpr float SPHERE_RELATIVE_DISTANCE = 0.5;
+	/**
+	 * Used in finding the proper distance between the parts
+	 * The maximal allowed value for
+	 * maximal radius of the node / sphere radius
+	 */
+	static const int MAX_DIAMETER_QUOTIENT = 30;
+	/**
+	 * The tolerance of the value of distance between parts
+	 */
+	static constexpr double SPHERE_DISTANCE_TOLERANCE = 0.96;
+
+	static constexpr double SPHERE_SIZE_QUOTIENT = 0.25;
+
+
+	static bool isInsidePart(Part::Shape shape, const Pt3D &partRadii, const Pt3D &center, double sphereRadius)
+	{
+		if(sphereRadius >= fS_Utils::min3(partRadii))
+			return true; // TODO Special case
+		Pt3D mostRemote = Pt3D(fabs(center.x), fabs(center.y), fabs(center.z)) + Pt3D(sphereRadius);
+		bool isInEllipsis;
+
+		bool result = false;
+
+		switch (shape)
+		{
+			case Part::Shape::SHAPE_CUBOID:
+				if(mostRemote.x <= partRadii.x && mostRemote.y <= partRadii.y && mostRemote.z <= partRadii.z)
+					result = true;
+				break;
+			case Part::Shape::SHAPE_CYLINDER:
+				isInEllipsis = pow(center.y / (partRadii.y - sphereRadius), 2) + pow(center.z / (partRadii.z - sphereRadius), 2) <= 1.0;
+				if (mostRemote.x <= partRadii.x && isInEllipsis)
+					result = true;
+				break;
+			case Part::Shape::SHAPE_ELLIPSOID:
+				if(pow(center.x / (partRadii.x - sphereRadius), 2) + pow(center.y / (partRadii.y - sphereRadius), 2) + pow(center.z / (partRadii.z - sphereRadius), 2) <= 1.0)
+					result = true;
+				break;
+			default:
+				logMessage("fS", "calculateVolume", LOG_ERROR, "Invalid part type");
+		}
+		return result;
+	}
+
+public:
+	static vector<Pt3D> findSphereCenters(Part::Shape shape, double &sphereRadius, const Pt3D &radii, const Pt3D &rotations)
+	{
+		double sphereRelativeDistance = SPHERE_RELATIVE_DISTANCE;
+		double minRadius = fS_Utils::min3(radii);
+		if(minRadius <= 0)
+			throw fS_Exception("Invalid part size in PartDistanceEstimator", 0);
+		double maxRadius = fS_Utils::max3(radii);
+		sphereRadius = SPHERE_SIZE_QUOTIENT * minRadius;
+		if (MAX_DIAMETER_QUOTIENT < maxRadius / sphereRadius)
+		{
+			// When max radius is much bigger than sphereRadius and there are to many spheresZ
+			sphereRelativeDistance = 1.0;   // Make the spheres adjacent to speed up the computation
+			sphereRadius = maxRadius / MAX_DIAMETER_QUOTIENT;
+		}
+		else if(shape == Part::Shape::SHAPE_ELLIPSOID)
+			sphereRadius = minRadius;
+
+		double sphereDiameter = 2 * sphereRadius;
+
+		double radiiArr[3]{radii.x, radii.y, radii.z};
+		vector<double> coordinates[3];
+		for (int dim = 0; dim < 3; dim++)
+		{
+			double spaceForSphereCenters = 2 * radiiArr[dim] - sphereDiameter;
+			if (spaceForSphereCenters > 0)
+			{
+				int lastIndex = ceil(spaceForSphereCenters / (sphereDiameter * sphereRelativeDistance));
+				for (int i = 0; i <= lastIndex; i++)
+				{
+					coordinates[dim].push_back(spaceForSphereCenters * (double(i) / lastIndex - 0.5));
+				}
+			}
+			else
+				coordinates[dim].push_back(0.0);
+		}
+
+		vector<Pt3D> centers;
+		for (int xi = 0; xi < int(coordinates[0].size()); xi++)
+		{
+			for (int yi = 0; yi < int(coordinates[1].size()); yi++)
+			{
+				for (int zi = 0; zi < int(coordinates[2].size()); zi++)
+				{
+						Pt3D center = Pt3D(coordinates[0][xi], coordinates[1][yi], coordinates[2][zi]);
+						if (isInsidePart(shape, radii, center, sphereRadius))
+						{
+							fS_Utils::rotateVector(center, rotations);
+							centers.push_back(center);
+						}
+				}
+			}
+		}
+		return centers;
+	}
+
+	static int isCollision(vector<Pt3D> &centersParent, vector<Pt3D> &centers, Pt3D &vectorBetweenParts,
+					double distanceThreshold)
+	{
+		double upperThresholdSq = distanceThreshold * distanceThreshold;
+		double lowerThresholdSq = pow(SPHERE_DISTANCE_TOLERANCE * distanceThreshold, 2);
+		double distanceSq;
+		double dx, dy, dz;
+		bool existsAdjacent = false;
+		Pt3D *tmpPoint;
+		for (int sc = 0; sc < int(centers.size()); sc++)
+		{
+			Pt3D shiftedSphere = Pt3D(centers[sc]);
+			shiftedSphere += vectorBetweenParts;
+			for (int psc = 0; psc < int(centersParent.size()); psc++)
+			{
+				tmpPoint = &centersParent[psc];
+				dx = shiftedSphere.x - tmpPoint->x;
+				dy = shiftedSphere.y - tmpPoint->y;
+				dz = shiftedSphere.z - tmpPoint->z;
+				distanceSq = dx * dx + dy * dy + dz * dz;
+
+				if (distanceSq <= upperThresholdSq)
+				{
+					if (distanceSq >= lowerThresholdSq)
+						existsAdjacent = true;
+					else
+					{
+						return COLLISION;
+					}
+				}
+			}
+		}
+		if (existsAdjacent)
+			return ADJACENT;
+		else
+			return DISJOINT;
+	}
+};
+
+double Node::getDistance()
+{
+	Pt3D size = calculateSize();
+	Pt3D parentSize = parent->calculateSize();	// Here we are sure that parent is not nullptr
+	double parentSphereRadius, sphereRadius;
+	vector<Pt3D> centersParent = PartDistanceEstimator::findSphereCenters(parent->partType, parentSphereRadius, parentSize, parent->getRotation());
+	vector<Pt3D> centers = PartDistanceEstimator::findSphereCenters(partType, sphereRadius, size, getRotation());
+
+	double distanceThreshold = sphereRadius + parentSphereRadius;
+	double minDistance = 0.0;
+	double maxDistance = 2 * (fS_Utils::max3(parentSize) + fS_Utils::max3(size));
+	double currentDistance = fS_Utils::avg(maxDistance, minDistance);
+	int result = -1;
+	int iterationNo = 0;
+	while (result != ADJACENT)
+	{
+		iterationNo++;
+		Pt3D vectorBetweenParts = state->v * currentDistance;
+		result = PartDistanceEstimator::isCollision(centersParent, centers, vectorBetweenParts, distanceThreshold);
+
+		if (result == DISJOINT)
+		{
+			maxDistance = currentDistance;
+			currentDistance = fS_Utils::avg(currentDistance, minDistance);
+		} else if (result == COLLISION)
+		{
+			minDistance = currentDistance;
+			currentDistance = fS_Utils::avg(maxDistance, currentDistance);
+		}
+
+		if(maxDistance <= 0 || iterationNo > 1000)
+			throw fS_Exception("Computing of distances between parts failed", 0);
+		if (currentDistance > maxDistance)
+		{
+			throw fS_Exception("Internal error; the maximal distance between parts exceeded.", 0);
+		}
+		if (currentDistance < minDistance)
+			throw fS_Exception("Internal error; the minimal distance between parts exceeded.", 0);
+
+	}
+
+	return currentDistance;
+}
+
+#endif //_PART_DISTANCE_ESTIMATOR_H_
