#include "boundary.hpp"

#define DEBUG 1
#define DoSingleStepOutput 0
#define DoTecplotOutput 1
#define DoRaster3DOutput 1
#define DoVRMLOutput 1
#define TecplotSuffix ".dat"
#define Raster3DSuffix ".r3d"
#define VRMLSUffix ".wrl"
#define HULLEPSILON MYEPSILON

// ======================================== Points on Boundary =================================

BoundaryPointSet::BoundaryPointSet()
{
	LinesCount = 0;
	Nr = -1;
}
;

BoundaryPointSet::BoundaryPointSet(atom *Walker)
{
	node = Walker;
	LinesCount = 0;
	Nr = Walker->nr;
}
;

BoundaryPointSet::~BoundaryPointSet()
{
	cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
	if (!lines.empty())
		cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
	node = NULL;
}
;

void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
{
	cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
			<< endl;
	if (line->endpoints[0] == this)
		{
			lines.insert(LinePair(line->endpoints[1]->Nr, line));
		}
	else
		{
			lines.insert(LinePair(line->endpoints[0]->Nr, line));
		}
	LinesCount++;
}
;

ostream &
operator <<(ostream &ost, BoundaryPointSet &a)
{
	ost << "[" << a.Nr << "|" << a.node->Name << "]";
	return ost;
}
;

// ======================================== Lines on Boundary =================================

BoundaryLineSet::BoundaryLineSet()
{
	for (int i = 0; i < 2; i++)
		endpoints[i] = NULL;
	TrianglesCount = 0;
	Nr = -1;
}
;

BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
{
	// set number
	Nr = number;
	// set endpoints in ascending order
	SetEndpointsOrdered(endpoints, Point[0], Point[1]);
	// add this line to the hash maps of both endpoints
	Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
	Point[1]->AddLine(this); //
	// clear triangles list
	TrianglesCount = 0;
	cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
}
;

BoundaryLineSet::~BoundaryLineSet()
{
	int Numbers[2];
	Numbers[0] = endpoints[1]->Nr;
	Numbers[1] = endpoints[0]->Nr;
	for (int i = 0; i < 2; i++) {
		cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
		endpoints[i]->lines.erase(Numbers[i]);
		if (endpoints[i]->lines.empty()) {
			cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
			if (endpoints[i] != NULL) {
				delete(endpoints[i]);
				endpoints[i] = NULL;
			} else
				cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
		} else
			cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
	}
	if (!triangles.empty())
		cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
}
;

void
BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
{
	cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
			<< endl;
	triangles.insert(TrianglePair(triangle->Nr, triangle));
	TrianglesCount++;
}
;

ostream &
operator <<(ostream &ost, BoundaryLineSet &a)
{
	ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
			<< a.endpoints[1]->node->Name << "]";
	return ost;
}
;

// ======================================== Triangles on Boundary =================================


BoundaryTriangleSet::BoundaryTriangleSet()
{
	for (int i = 0; i < 3; i++)
		{
			endpoints[i] = NULL;
			lines[i] = NULL;
		}
	Nr = -1;
}
;

BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3],
		int number)
{
	// set number
	Nr = number;
	// set lines
	cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
	for (int i = 0; i < 3; i++)
		{
			lines[i] = line[i];
			lines[i]->AddTriangle(this);
		}
	// get ascending order of endpoints
	map<int, class BoundaryPointSet *> OrderMap;
	for (int i = 0; i < 3; i++)
		// for all three lines
		for (int j = 0; j < 2; j++)
			{ // for both endpoints
				OrderMap.insert(pair<int, class BoundaryPointSet *> (
						line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
				// and we don't care whether insertion fails
			}
	// set endpoints
	int Counter = 0;
	cout << Verbose(6) << " with end points ";
	for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
			!= OrderMap.end(); runner++)
		{
			endpoints[Counter] = runner->second;
			cout << " " << *endpoints[Counter];
			Counter++;
		}
	if (Counter < 3)
		{
			cerr << "ERROR! We have a triangle with only two distinct endpoints!"
					<< endl;
			//exit(1);
		}
	cout << "." << endl;
}
;

BoundaryTriangleSet::~BoundaryTriangleSet()
{
	for (int i = 0; i < 3; i++) {
		cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
		lines[i]->triangles.erase(Nr);
		if (lines[i]->triangles.empty()) {
			cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
			if (lines[i] != NULL) {
				delete (lines[i]);
				lines[i] = NULL;
			} else
				cerr << "ERROR: This line " << i << " has already been free'd." << endl;
		} else
			cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
	}
}
;

void
BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
{
	// get normal vector
	NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
			&endpoints[2]->node->x);

	// make it always point inward (any offset vector onto plane projected onto normal vector suffices)
	if (NormalVector.Projection(&OtherVector) > 0)
		NormalVector.Scale(-1.);
}
;

ostream &
operator <<(ostream &ost, BoundaryTriangleSet &a)
{
	ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
			<< a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
	return ost;
}
;

// ========================================== F U N C T I O N S =================================

/** Finds the endpoint two lines are sharing.
 * \param *line1 first line
 * \param *line2 second line
 * \return point which is shared or NULL if none
 */
class BoundaryPointSet *
GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
{
	class BoundaryLineSet * lines[2] =
		{ line1, line2 };
	class BoundaryPointSet *node = NULL;
	map<int, class BoundaryPointSet *> OrderMap;
	pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
	for (int i = 0; i < 2; i++)
		// for both lines
		for (int j = 0; j < 2; j++)
			{ // for both endpoints
				OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
						lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
				if (!OrderTest.second)
					{ // if insertion fails, we have common endpoint
						node = OrderTest.first->second;
						cout << Verbose(5) << "Common endpoint of lines " << *line1
								<< " and " << *line2 << " is: " << *node << "." << endl;
						j = 2;
						i = 2;
						break;
					}
			}
	return node;
}
;

/** Determines the boundary points of a cluster.
 * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
 * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
 * center and first and last point in the triple, it is thrown out.
 * \param *out output stream for debugging
 * \param *mol molecule structure representing the cluster
 */
Boundaries *
GetBoundaryPoints(ofstream *out, molecule *mol)
{
	atom *Walker = NULL;
	PointMap PointsOnBoundary;
	LineMap LinesOnBoundary;
	TriangleMap TrianglesOnBoundary;

	*out << Verbose(1) << "Finding all boundary points." << endl;
	Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
	BoundariesTestPair BoundaryTestPair;
	Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
	double radius, angle;
	// 3a. Go through every axis
	for (int axis = 0; axis < NDIM; axis++)
		{
			AxisVector.Zero();
			AngleReferenceVector.Zero();
			AngleReferenceNormalVector.Zero();
			AxisVector.x[axis] = 1.;
			AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
			AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
			//		*out << Verbose(1) << "Axisvector is ";
			//		AxisVector.Output(out);
			//		*out << " and AngleReferenceVector is ";
			//		AngleReferenceVector.Output(out);
			//		*out << "." << endl;
			//		*out << " and AngleReferenceNormalVector is ";
			//		AngleReferenceNormalVector.Output(out);
			//		*out << "." << endl;
			// 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
			Walker = mol->start;
			while (Walker->next != mol->end)
				{
					Walker = Walker->next;
					Vector ProjectedVector;
					ProjectedVector.CopyVector(&Walker->x);
					ProjectedVector.ProjectOntoPlane(&AxisVector);
					// correct for negative side
					//if (Projection(y) < 0)
					//angle = 2.*M_PI - angle;
					radius = ProjectedVector.Norm();
					if (fabs(radius) > MYEPSILON)
						angle = ProjectedVector.Angle(&AngleReferenceVector);
					else
						angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues

					//*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
					if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
						{
							angle = 2. * M_PI - angle;
						}
					//*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
					//ProjectedVector.Output(out);
					//*out << endl;
					BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
							DistancePair (radius, Walker)));
					if (BoundaryTestPair.second)
						{ // successfully inserted
						}
					else
						{ // same point exists, check first r, then distance of original vectors to center of gravity
							*out << Verbose(2)
									<< "Encountered two vectors whose projection onto axis "
									<< axis << " is equal: " << endl;
							*out << Verbose(2) << "Present vector: ";
							BoundaryTestPair.first->second.second->x.Output(out);
							*out << endl;
							*out << Verbose(2) << "New vector: ";
							Walker->x.Output(out);
							*out << endl;
							double tmp = ProjectedVector.Norm();
							if (tmp > BoundaryTestPair.first->second.first)
								{
									BoundaryTestPair.first->second.first = tmp;
									BoundaryTestPair.first->second.second = Walker;
									*out << Verbose(2) << "Keeping new vector." << endl;
								}
							else if (tmp == BoundaryTestPair.first->second.first)
								{
									if (BoundaryTestPair.first->second.second->x.ScalarProduct(
											&BoundaryTestPair.first->second.second->x)
											< Walker->x.ScalarProduct(&Walker->x))
										{ // Norm() does a sqrt, which makes it a lot slower
											BoundaryTestPair.first->second.second = Walker;
											*out << Verbose(2) << "Keeping new vector." << endl;
										}
									else
										{
											*out << Verbose(2) << "Keeping present vector." << endl;
										}
								}
							else
								{
									*out << Verbose(2) << "Keeping present vector." << endl;
								}
						}
				}
			// printing all inserted for debugging
			//		{
			//			*out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
			//			int i=0;
			//			for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
			//				if (runner != BoundaryPoints[axis].begin())
			//					*out << ", " << i << ": " << *runner->second.second;
			//				else
			//					*out << i << ": " << *runner->second.second;
			//				i++;
			//			}
			//			*out << endl;
			//		}
			// 3c. throw out points whose distance is less than the mean of left and right neighbours
			bool flag = false;
			do
				{ // do as long as we still throw one out per round
					*out << Verbose(1)
							<< "Looking for candidates to kick out by convex condition ... "
							<< endl;
					flag = false;
					Boundaries::iterator left = BoundaryPoints[axis].end();
					Boundaries::iterator right = BoundaryPoints[axis].end();
					for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
							!= BoundaryPoints[axis].end(); runner++)
						{
							// set neighbours correctly
							if (runner == BoundaryPoints[axis].begin())
								{
									left = BoundaryPoints[axis].end();
								}
							else
								{
									left = runner;
								}
							left--;
							right = runner;
							right++;
							if (right == BoundaryPoints[axis].end())
								{
									right = BoundaryPoints[axis].begin();
								}
							// check distance

							// construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
								{
									Vector SideA, SideB, SideC, SideH;
									SideA.CopyVector(&left->second.second->x);
									SideA.ProjectOntoPlane(&AxisVector);
									//					*out << "SideA: ";
									//					SideA.Output(out);
									//					*out << endl;

									SideB.CopyVector(&right->second.second->x);
									SideB.ProjectOntoPlane(&AxisVector);
									//					*out << "SideB: ";
									//					SideB.Output(out);
									//					*out << endl;

									SideC.CopyVector(&left->second.second->x);
									SideC.SubtractVector(&right->second.second->x);
									SideC.ProjectOntoPlane(&AxisVector);
									//					*out << "SideC: ";
									//					SideC.Output(out);
									//					*out << endl;

									SideH.CopyVector(&runner->second.second->x);
									SideH.ProjectOntoPlane(&AxisVector);
									//					*out << "SideH: ";
									//					SideH.Output(out);
									//					*out << endl;

									// calculate each length
									double a = SideA.Norm();
									//double b = SideB.Norm();
									//double c = SideC.Norm();
									double h = SideH.Norm();
									// calculate the angles
									double alpha = SideA.Angle(&SideH);
									double beta = SideA.Angle(&SideC);
									double gamma = SideB.Angle(&SideH);
									double delta = SideC.Angle(&SideH);
									double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
											< M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
									//					*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
									//*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
									if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
											< MYEPSILON) && (h < MinDistance))
										{
											// throw out point
											//*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
											BoundaryPoints[axis].erase(runner);
											flag = true;
										}
								}
						}
				}
			while (flag);
		}
	return BoundaryPoints;
}
;

/** Determines greatest diameters of a cluster defined by its convex envelope.
 * Looks at lines parallel to one axis and where they intersect on the projected planes
 * \param *out output stream for debugging
 * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
 * \param *mol molecule structure representing the cluster
 * \param IsAngstroem whether we have angstroem or atomic units
 * \return NDIM array of the diameters
 */
double *
GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
		bool IsAngstroem)
{
	// get points on boundary of NULL was given as parameter
	bool BoundaryFreeFlag = false;
	Boundaries *BoundaryPoints = BoundaryPtr;
	if (BoundaryPoints == NULL)
		{
			BoundaryFreeFlag = true;
			BoundaryPoints = GetBoundaryPoints(out, mol);
		}
	else
		{
			*out << Verbose(1) << "Using given boundary points set." << endl;
		}
	// determine biggest "diameter" of cluster for each axis
	Boundaries::iterator Neighbour, OtherNeighbour;
	double *GreatestDiameter = new double[NDIM];
	for (int i = 0; i < NDIM; i++)
		GreatestDiameter[i] = 0.;
	double OldComponent, tmp, w1, w2;
	Vector DistanceVector, OtherVector;
	int component, Othercomponent;
	for (int axis = 0; axis < NDIM; axis++)
		{ // regard each projected plane
			//*out << Verbose(1) << "Current axis is " << axis << "." << endl;
			for (int j = 0; j < 2; j++)
				{ // and for both axis on the current plane
					component = (axis + j + 1) % NDIM;
					Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
					//*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
					for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
							!= BoundaryPoints[axis].end(); runner++)
						{
							//*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
							// seek for the neighbours pair where the Othercomponent sign flips
							Neighbour = runner;
							Neighbour++;
							if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
								Neighbour = BoundaryPoints[axis].begin();
							DistanceVector.CopyVector(&runner->second.second->x);
							DistanceVector.SubtractVector(&Neighbour->second.second->x);
							do
								{ // seek for neighbour pair where it flips
									OldComponent = DistanceVector.x[Othercomponent];
									Neighbour++;
									if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
										Neighbour = BoundaryPoints[axis].begin();
									DistanceVector.CopyVector(&runner->second.second->x);
									DistanceVector.SubtractVector(&Neighbour->second.second->x);
									//*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
								}
							while ((runner != Neighbour) && (fabs(OldComponent / fabs(
									OldComponent) - DistanceVector.x[Othercomponent] / fabs(
									DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
							if (runner != Neighbour)
								{
									OtherNeighbour = Neighbour;
									if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
										OtherNeighbour = BoundaryPoints[axis].end();
									OtherNeighbour--;
									//*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
									// now we have found the pair: Neighbour and OtherNeighbour
									OtherVector.CopyVector(&runner->second.second->x);
									OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
									//*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
									//*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
									// do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
									w1 = fabs(OtherVector.x[Othercomponent]);
									w2 = fabs(DistanceVector.x[Othercomponent]);
									tmp = fabs((w1 * DistanceVector.x[component] + w2
											* OtherVector.x[component]) / (w1 + w2));
									// mark if it has greater diameter
									//*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
									GreatestDiameter[component] = (GreatestDiameter[component]
											> tmp) ? GreatestDiameter[component] : tmp;
								} //else
							//*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
						}
				}
		}
	*out << Verbose(0) << "RESULT: The biggest diameters are "
			<< GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
			<< GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
			: "atomiclength") << "." << endl;

	// free reference lists
	if (BoundaryFreeFlag)
		delete[] (BoundaryPoints);

	return GreatestDiameter;
}
;

/** Creates the objects in a VRML file.
 * \param *out output stream for debugging
 * \param *vrmlfile output stream for tecplot data
 * \param *Tess Tesselation structure with constructed triangles
 * \param *mol molecule structure with atom positions
 */
void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
{
	atom *Walker = mol->start;
	bond *Binder = mol->first;
	int i;
	Vector *center = mol->DetermineCenterOfAll(out);
	if (vrmlfile != NULL) {
		//cout << Verbose(1) << "Writing Raster3D file ... ";
		*vrmlfile << "#VRML V2.0 utf8" << endl;
		*vrmlfile << "#Created by molecuilder" << endl;
		*vrmlfile << "#All atoms as spheres" << endl;
		while (Walker->next != mol->end) {
			Walker = Walker->next;
			*vrmlfile << "Sphere {" << endl << "	"; // 2 is sphere type
			for (i=0;i<NDIM;i++)
				*vrmlfile << Walker->x.x[i]+center->x[i] << " ";
			*vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
		}

		*vrmlfile << "# All bonds as vertices" << endl;
		while (Binder->next != mol->last) {
			Binder = Binder->next;
			*vrmlfile << "3" << endl << "	"; // 2 is round-ended cylinder type
			for (i=0;i<NDIM;i++)
				*vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
			*vrmlfile << "\t0.03\t";
			for (i=0;i<NDIM;i++)
				*vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
			*vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
		}

		*vrmlfile << "# All tesselation triangles" << endl;
		for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
			*vrmlfile << "1" << endl << "	"; // 1 is triangle type
			for (i=0;i<3;i++) { // print each node
				for (int j=0;j<NDIM;j++)	// and for each node all NDIM coordinates
					*vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
				*vrmlfile << "\t";
			}
			*vrmlfile << "1. 0. 0." << endl;	// red as colour
			*vrmlfile << "18" << endl << "	0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
		}
	} else {
		cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
	}
	delete(center);
};

/** Creates the objects in a raster3d file (renderable with a header.r3d).
 * \param *out output stream for debugging
 * \param *rasterfile output stream for tecplot data
 * \param *Tess Tesselation structure with constructed triangles
 * \param *mol molecule structure with atom positions
 */
void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
{
	atom *Walker = mol->start;
	bond *Binder = mol->first;
	int i;
	Vector *center = mol->DetermineCenterOfAll(out);
	if (rasterfile != NULL) {
		//cout << Verbose(1) << "Writing Raster3D file ... ";
		*rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
		*rasterfile << "@header.r3d" << endl;
		*rasterfile << "# All atoms as spheres" << endl;
		while (Walker->next != mol->end) {
			Walker = Walker->next;
			*rasterfile << "2" << endl << "	";	// 2 is sphere type
			for (i=0;i<NDIM;i++)
				*rasterfile << Walker->x.x[i]+center->x[i] << " ";
			*rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
		}

		*rasterfile << "# All bonds as vertices" << endl;
		while (Binder->next != mol->last) {
			Binder = Binder->next;
			*rasterfile << "3" << endl << "	";	// 2 is round-ended cylinder type
			for (i=0;i<NDIM;i++)
				*rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
			*rasterfile << "\t0.03\t";
			for (i=0;i<NDIM;i++)
				*rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
			*rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
		}

		*rasterfile << "# All tesselation triangles" << endl;
		*rasterfile << "8\n	25. -1.	 1. 1. 1.	 0.0		0 0 0 2\n	SOLID		 1.0 0.0 0.0\n	BACKFACE	0.3 0.3 1.0	 0 0\n";
		for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
			*rasterfile << "1" << endl << "	";	// 1 is triangle type
			for (i=0;i<3;i++) {	// print each node
				for (int j=0;j<NDIM;j++)	// and for each node all NDIM coordinates
					*rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
				*rasterfile << "\t";
			}
			*rasterfile << "1. 0. 0." << endl;	// red as colour
			//*rasterfile << "18" << endl << "	0.5 0.5 0.5" << endl;	// 18 is transparency type for previous object
		}
		*rasterfile << "9\n	terminating special property\n";
	} else {
		cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
	}
	delete(center);
};

/** This function creates the tecplot file, displaying the tesselation of the hull.
 * \param *out output stream for debugging
 * \param *tecplot output stream for tecplot data
 * \param N arbitrary number to differentiate various zones in the tecplot format
 */
void
write_tecplot_file(ofstream *out, ofstream *tecplot,
		class Tesselation *TesselStruct, class molecule *mol, int N)
{
	if (tecplot != NULL)
		{
			*tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
			*tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
			*tecplot << "ZONE T=\"TRIANGLES" << N << "\", N="
					<< TesselStruct->PointsOnBoundaryCount << ", E="
					<< TesselStruct->TrianglesOnBoundaryCount
					<< ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
			int *LookupList = new int[mol->AtomCount];
			for (int i = 0; i < mol->AtomCount; i++)
				LookupList[i] = -1;

			// print atom coordinates
			*out << Verbose(2) << "The following triangles were created:";
			int Counter = 1;
			atom *Walker = NULL;
			for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target
					!= TesselStruct->PointsOnBoundary.end(); target++)
				{
					Walker = target->second->node;
					LookupList[Walker->nr] = Counter++;
					*tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " "
							<< Walker->x.x[2] << " " << endl;
				}
			*tecplot << endl;
			// print connectivity
			for (TriangleMap::iterator runner =
					TesselStruct->TrianglesOnBoundary.begin(); runner
					!= TesselStruct->TrianglesOnBoundary.end(); runner++)
				{
					*out << " " << runner->second->endpoints[0]->node->Name << "<->"
							<< runner->second->endpoints[1]->node->Name << "<->"
							<< runner->second->endpoints[2]->node->Name;
					*tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " "
							<< LookupList[runner->second->endpoints[1]->node->nr] << " "
							<< LookupList[runner->second->endpoints[2]->node->nr] << endl;
				}
			delete[] (LookupList);
			*out << endl;
		}
}

/** Determines the volume of a cluster.
 * Determines first the convex envelope, then tesselates it and calculates its volume.
 * \param *out output stream for debugging
 * \param *filename filename prefix for output of vertex data
 * \param *configuration needed for path to store convex envelope file
 * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
 * \param *mol molecule structure representing the cluster
 * \return determined volume of the cluster in cubed config:GetIsAngstroem()
 */
double
VolumeOfConvexEnvelope(ofstream *out, const char *filename, config *configuration,
		Boundaries *BoundaryPtr, molecule *mol)
{
	bool IsAngstroem = configuration->GetIsAngstroem();
	atom *Walker = NULL;
	struct Tesselation *TesselStruct = new Tesselation;
	bool BoundaryFreeFlag = false;
	Boundaries *BoundaryPoints = BoundaryPtr;
	double volume = 0.;
	double PyramidVolume = 0.;
	double G, h;
	Vector x, y;
	double a, b, c;

	//Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.

	// 1. calculate center of gravity
	*out << endl;
	Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);

	// 2. translate all points into CoG
	*out << Verbose(1) << "Translating system to Center of Gravity." << endl;
	Walker = mol->start;
	while (Walker->next != mol->end)
		{
			Walker = Walker->next;
			Walker->x.Translate(CenterOfGravity);
		}

	// 3. Find all points on the boundary
	if (BoundaryPoints == NULL)
		{
			BoundaryFreeFlag = true;
			BoundaryPoints = GetBoundaryPoints(out, mol);
		}
	else
		{
			*out << Verbose(1) << "Using given boundary points set." << endl;
		}

	// 4. fill the boundary point list
	for (int axis = 0; axis < NDIM; axis++)
		for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
				!= BoundaryPoints[axis].end(); runner++)
			{
				TesselStruct->AddPoint(runner->second.second);
			}

	*out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount
			<< " points on the convex boundary." << endl;
	// now we have the whole set of edge points in the BoundaryList

	// listing for debugging
	//	*out << Verbose(1) << "Listing PointsOnBoundary:";
	//	for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
	//		*out << " " << *runner->second;
	//	}
	//	*out << endl;

	// 5a. guess starting triangle
	TesselStruct->GuessStartingTriangle(out);

	// 5b. go through all lines, that are not yet part of two triangles (only of one so far)
	TesselStruct->TesselateOnBoundary(out, configuration, mol);

	*out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
			<< " triangles with " << TesselStruct->LinesOnBoundaryCount
			<< " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
			<< endl;

	// 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
	*out << Verbose(1)
			<< "Calculating the volume of the pyramids formed out of triangles and center of gravity."
			<< endl;
	for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
			!= TesselStruct->TrianglesOnBoundary.end(); runner++)
		{ // go through every triangle, calculate volume of its pyramid with CoG as peak
			x.CopyVector(&runner->second->endpoints[0]->node->x);
			x.SubtractVector(&runner->second->endpoints[1]->node->x);
			y.CopyVector(&runner->second->endpoints[0]->node->x);
			y.SubtractVector(&runner->second->endpoints[2]->node->x);
			a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
					&runner->second->endpoints[1]->node->x));
			b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
					&runner->second->endpoints[2]->node->x));
			c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
					&runner->second->endpoints[1]->node->x));
			G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
			x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
					&runner->second->endpoints[1]->node->x,
					&runner->second->endpoints[2]->node->x);
			x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
			h = x.Norm(); // distance of CoG to triangle
			PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
			*out << Verbose(2) << "Area of triangle is " << G << " "
					<< (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
					<< h << " and the volume is " << PyramidVolume << " "
					<< (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
			volume += PyramidVolume;
		}
	*out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)
			<< volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
			<< endl;

	// 7. translate all points back from CoG
	*out << Verbose(1) << "Translating system back from Center of Gravity."
			<< endl;
	CenterOfGravity->Scale(-1);
	Walker = mol->start;
	while (Walker->next != mol->end)
		{
			Walker = Walker->next;
			Walker->x.Translate(CenterOfGravity);
		}

	// 8. Store triangles in tecplot file
	string OutputName(filename);
	OutputName.append(TecplotSuffix);
	ofstream *tecplot = new ofstream(OutputName.c_str());
	write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
	tecplot->close();
	delete(tecplot);

	// free reference lists
	if (BoundaryFreeFlag)
		delete[] (BoundaryPoints);

	return volume;
}
;

/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
 * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
 * \param *out output stream for debugging
 * \param *configuration needed for path to store convex envelope file
 * \param *mol molecule structure representing the cluster
 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
 * \param celldensity desired average density in final cell
 */
void
PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
		double ClusterVolume, double celldensity)
{
	// transform to PAS
	mol->PrincipalAxisSystem(out, true);

	// some preparations beforehand
	bool IsAngstroem = configuration->GetIsAngstroem();
	Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
	double clustervolume;
	if (ClusterVolume == 0)
		clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
				BoundaryPoints, mol);
	else
		clustervolume = ClusterVolume;
	double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
			IsAngstroem);
	Vector BoxLengths;
	int repetition[NDIM] =
		{ 1, 1, 1 };
	int TotalNoClusters = 1;
	for (int i = 0; i < NDIM; i++)
		TotalNoClusters *= repetition[i];

	// sum up the atomic masses
	double totalmass = 0.;
	atom *Walker = mol->start;
	while (Walker->next != mol->end)
		{
			Walker = Walker->next;
			totalmass += Walker->type->mass;
		}
	*out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
			<< totalmass << " atomicmassunit." << endl;

	*out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
			<< totalmass / clustervolume << " atomicmassunit/"
			<< (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;

	// solve cubic polynomial
	*out << Verbose(1) << "Solving equidistant suspension in water problem ..."
			<< endl;
	double cellvolume;
	if (IsAngstroem)
		cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
				/ clustervolume)) / (celldensity - 1);
	else
		cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
				/ clustervolume)) / (celldensity - 1);
	*out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
			<< " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"
			: "atomiclength") << "^3." << endl;

	double minimumvolume = TotalNoClusters * (GreatestDiameter[0]
			* GreatestDiameter[1] * GreatestDiameter[2]);
	*out << Verbose(1)
			<< "Minimum volume of the convex envelope contained in a rectangular box is "
			<< minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"
			: "atomiclength") << "^3." << endl;
	if (minimumvolume > cellvolume)
		{
			cerr << Verbose(0)
					<< "ERROR: the containing box already has a greater volume than the envisaged cell volume!"
					<< endl;
			cout << Verbose(0)
					<< "Setting Box dimensions to minimum possible, the greatest diameters."
					<< endl;
			for (int i = 0; i < NDIM; i++)
				BoxLengths.x[i] = GreatestDiameter[i];
			mol->CenterEdge(out, &BoxLengths);
		}
	else
		{
			BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]
					* GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
			BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]
					* GreatestDiameter[1] + repetition[0] * repetition[2]
					* GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]
					* repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
			BoxLengths.x[2] = minimumvolume - cellvolume;
			double x0 = 0., x1 = 0., x2 = 0.;
			if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],
					BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
				*out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
						<< " ." << endl;
			else
				{
					*out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
							<< " and " << x1 << " and " << x2 << " ." << endl;
					x0 = x2; // sorted in ascending order
				}

			cellvolume = 1;
			for (int i = 0; i < NDIM; i++)
				{
					BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
					cellvolume *= BoxLengths.x[i];
				}

			// set new box dimensions
			*out << Verbose(0) << "Translating to box with these boundaries." << endl;
			mol->CenterInBox((ofstream *) &cout, &BoxLengths);
		}
	// update Box of atoms by boundary
	mol->SetBoxDimension(&BoxLengths);
	*out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
			<< BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
			<< BoxLengths.x[2] << " with total volume of " << cellvolume << " "
			<< (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
}
;

// =========================================================== class TESSELATION ===========================================

/** Constructor of class Tesselation.
 */
Tesselation::Tesselation()
{
	PointsOnBoundaryCount = 0;
	LinesOnBoundaryCount = 0;
	TrianglesOnBoundaryCount = 0;
	TriangleFilesWritten = 0;
}
;

/** Constructor of class Tesselation.
 * We have to free all points, lines and triangles.
 */
Tesselation::~Tesselation()
{
	cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
	for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
		if (runner->second != NULL) {
			delete (runner->second);
			runner->second = NULL;
		} else
			cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
	}
}
;

/** Gueses first starting triangle of the convex envelope.
 * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
 * \param *out output stream for debugging
 * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
 */
void
Tesselation::GuessStartingTriangle(ofstream *out)
{
	// 4b. create a starting triangle
	// 4b1. create all distances
	DistanceMultiMap DistanceMMap;
	double distance, tmp;
	Vector PlaneVector, TrialVector;
	PointMap::iterator A, B, C; // three nodes of the first triangle
	A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily

	// with A chosen, take each pair B,C and sort
	if (A != PointsOnBoundary.end())
		{
			B = A;
			B++;
			for (; B != PointsOnBoundary.end(); B++)
				{
					C = B;
					C++;
					for (; C != PointsOnBoundary.end(); C++)
						{
							tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
							distance = tmp * tmp;
							tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
							distance += tmp * tmp;
							tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
							distance += tmp * tmp;
							DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
									PointMap::iterator, PointMap::iterator> (B, C)));
						}
				}
		}
	//		// listing distances
	//		*out << Verbose(1) << "Listing DistanceMMap:";
	//		for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
	//			*out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
	//		}
	//		*out << endl;
	// 4b2. pick three baselines forming a triangle
	// 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
	DistanceMultiMap::iterator baseline = DistanceMMap.begin();
	for (; baseline != DistanceMMap.end(); baseline++)
		{
			// we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
			// 2. next, we have to check whether all points reside on only one side of the triangle
			// 3. construct plane vector
			PlaneVector.MakeNormalVector(&A->second->node->x,
					&baseline->second.first->second->node->x,
					&baseline->second.second->second->node->x);
			*out << Verbose(2) << "Plane vector of candidate triangle is ";
			PlaneVector.Output(out);
			*out << endl;
			// 4. loop over all points
			double sign = 0.;
			PointMap::iterator checker = PointsOnBoundary.begin();
			for (; checker != PointsOnBoundary.end(); checker++)
				{
					// (neglecting A,B,C)
					if ((checker == A) || (checker == baseline->second.first) || (checker
							== baseline->second.second))
						continue;
					// 4a. project onto plane vector
					TrialVector.CopyVector(&checker->second->node->x);
					TrialVector.SubtractVector(&A->second->node->x);
					distance = TrialVector.Projection(&PlaneVector);
					if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
						continue;
					*out << Verbose(3) << "Projection of " << checker->second->node->Name
							<< " yields distance of " << distance << "." << endl;
					tmp = distance / fabs(distance);
					// 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
					if ((sign != 0) && (tmp != sign))
						{
							// 4c. If so, break 4. loop and continue with next candidate in 1. loop
							*out << Verbose(2) << "Current candidates: "
									<< A->second->node->Name << ","
									<< baseline->second.first->second->node->Name << ","
									<< baseline->second.second->second->node->Name << " leave "
									<< checker->second->node->Name << " outside the convex hull."
									<< endl;
							break;
						}
					else
						{ // note the sign for later
							*out << Verbose(2) << "Current candidates: "
									<< A->second->node->Name << ","
									<< baseline->second.first->second->node->Name << ","
									<< baseline->second.second->second->node->Name << " leave "
									<< checker->second->node->Name << " inside the convex hull."
									<< endl;
							sign = tmp;
						}
					// 4d. Check whether the point is inside the triangle (check distance to each node
					tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
					int innerpoint = 0;
					if ((tmp < A->second->node->x.DistanceSquared(
							&baseline->second.first->second->node->x)) && (tmp
							< A->second->node->x.DistanceSquared(
									&baseline->second.second->second->node->x)))
						innerpoint++;
					tmp = checker->second->node->x.DistanceSquared(
							&baseline->second.first->second->node->x);
					if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
							&A->second->node->x)) && (tmp
							< baseline->second.first->second->node->x.DistanceSquared(
									&baseline->second.second->second->node->x)))
						innerpoint++;
					tmp = checker->second->node->x.DistanceSquared(
							&baseline->second.second->second->node->x);
					if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
							&baseline->second.first->second->node->x)) && (tmp
							< baseline->second.second->second->node->x.DistanceSquared(
									&A->second->node->x)))
						innerpoint++;
					// 4e. If so, break 4. loop and continue with next candidate in 1. loop
					if (innerpoint == 3)
						break;
				}
			// 5. come this far, all on same side? Then break 1. loop and construct triangle
			if (checker == PointsOnBoundary.end())
				{
					*out << "Looks like we have a candidate!" << endl;
					break;
				}
		}
	if (baseline != DistanceMMap.end())
		{
			BPS[0] = baseline->second.first->second;
			BPS[1] = baseline->second.second->second;
			BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
			BPS[0] = A->second;
			BPS[1] = baseline->second.second->second;
			BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
			BPS[0] = baseline->second.first->second;
			BPS[1] = A->second;
			BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);

			// 4b3. insert created triangle
			BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
			TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
			TrianglesOnBoundaryCount++;
			for (int i = 0; i < NDIM; i++)
				{
					LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
					LinesOnBoundaryCount++;
				}

			*out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
		}
	else
		{
			*out << Verbose(1) << "No starting triangle found." << endl;
			exit(255);
		}
}
;

/** Tesselates the convex envelope of a cluster from a single starting triangle.
 * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
 * 2 triangles. Hence, we go through all current lines:
 * -# if the lines contains to only one triangle
 * -# We search all points in the boundary
 *		-# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors
 *		-# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
 *			 baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
 *		-# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
 * \param *out output stream for debugging
 * \param *configuration for IsAngstroem
 * \param *mol the cluster as a molecule structure
 */
void
Tesselation::TesselateOnBoundary(ofstream *out, config *configuration,
		molecule *mol)
{
	bool flag;
	PointMap::iterator winner;
	class BoundaryPointSet *peak = NULL;
	double SmallestAngle, TempAngle;
	Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
			PropagationVector;
	LineMap::iterator LineChecker[2];
	do
		{
			flag = false;
			for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
					!= LinesOnBoundary.end(); baseline++)
				if (baseline->second->TrianglesCount == 1)
					{
						*out << Verbose(2) << "Current baseline is between "
								<< *(baseline->second) << "." << endl;
						// 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
						SmallestAngle = M_PI;
						BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
						// get peak point with respect to this base line's only triangle
						for (int i = 0; i < 3; i++)
							if ((BTS->endpoints[i] != baseline->second->endpoints[0])
									&& (BTS->endpoints[i] != baseline->second->endpoints[1]))
								peak = BTS->endpoints[i];
						*out << Verbose(3) << " and has peak " << *peak << "." << endl;
						// normal vector of triangle
						BTS->GetNormalVector(NormalVector);
						*out << Verbose(4) << "NormalVector of base triangle is ";
						NormalVector.Output(out);
						*out << endl;
						// offset to center of triangle
						CenterVector.Zero();
						for (int i = 0; i < 3; i++)
							CenterVector.AddVector(&BTS->endpoints[i]->node->x);
						CenterVector.Scale(1. / 3.);
						*out << Verbose(4) << "CenterVector of base triangle is ";
						CenterVector.Output(out);
						*out << endl;
						// vector in propagation direction (out of triangle)
						// project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
						TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
						TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
						PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
						TempVector.CopyVector(&CenterVector);
						TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
						//*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
						if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
							PropagationVector.Scale(-1.);
						*out << Verbose(4) << "PropagationVector of base triangle is ";
						PropagationVector.Output(out);
						*out << endl;
						winner = PointsOnBoundary.end();
						for (PointMap::iterator target = PointsOnBoundary.begin(); target
								!= PointsOnBoundary.end(); target++)
							if ((target->second != baseline->second->endpoints[0])
									&& (target->second != baseline->second->endpoints[1]))
								{ // don't take the same endpoints
									*out << Verbose(3) << "Target point is " << *(target->second)
											<< ":";
									bool continueflag = true;

									VirtualNormalVector.CopyVector(
											&baseline->second->endpoints[0]->node->x);
									VirtualNormalVector.AddVector(
											&baseline->second->endpoints[0]->node->x);
									VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
									VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
									TempAngle = VirtualNormalVector.Angle(&PropagationVector);
									continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
									if (!continueflag)
										{
											*out << Verbose(4)
													<< "Angle between propagation direction and base line to "
													<< *(target->second) << " is " << TempAngle
													<< ", bad direction!" << endl;
											continue;
										}
									else
										*out << Verbose(4)
												<< "Angle between propagation direction and base line to "
												<< *(target->second) << " is " << TempAngle
												<< ", good direction!" << endl;
									LineChecker[0] = baseline->second->endpoints[0]->lines.find(
											target->first);
									LineChecker[1] = baseline->second->endpoints[1]->lines.find(
											target->first);
									//						if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
									//							*out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
									//						else
									//							*out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
									//						if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
									//							*out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
									//						else
									//							*out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
									// check first endpoint (if any connecting line goes to target or at least not more than 1)
									continueflag = continueflag && (((LineChecker[0]
											== baseline->second->endpoints[0]->lines.end())
											|| (LineChecker[0]->second->TrianglesCount == 1)));
									if (!continueflag)
										{
											*out << Verbose(4) << *(baseline->second->endpoints[0])
													<< " has line " << *(LineChecker[0]->second)
													<< " to " << *(target->second)
													<< " as endpoint with "
													<< LineChecker[0]->second->TrianglesCount
													<< " triangles." << endl;
											continue;
										}
									// check second endpoint (if any connecting line goes to target or at least not more than 1)
									continueflag = continueflag && (((LineChecker[1]
											== baseline->second->endpoints[1]->lines.end())
											|| (LineChecker[1]->second->TrianglesCount == 1)));
									if (!continueflag)
										{
											*out << Verbose(4) << *(baseline->second->endpoints[1])
													<< " has line " << *(LineChecker[1]->second)
													<< " to " << *(target->second)
													<< " as endpoint with "
													<< LineChecker[1]->second->TrianglesCount
													<< " triangles." << endl;
											continue;
										}
									// check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
									continueflag = continueflag && (!(((LineChecker[0]
											!= baseline->second->endpoints[0]->lines.end())
											&& (LineChecker[1]
													!= baseline->second->endpoints[1]->lines.end())
											&& (GetCommonEndpoint(LineChecker[0]->second,
													LineChecker[1]->second) == peak))));
									if (!continueflag)
										{
											*out << Verbose(4) << "Current target is peak!" << endl;
											continue;
										}
									// in case NOT both were found
									if (continueflag)
										{ // create virtually this triangle, get its normal vector, calculate angle
											flag = true;
											VirtualNormalVector.MakeNormalVector(
													&baseline->second->endpoints[0]->node->x,
													&baseline->second->endpoints[1]->node->x,
													&target->second->node->x);
											// make it always point inward
											if (baseline->second->endpoints[0]->node->x.Projection(
													&VirtualNormalVector) > 0)
												VirtualNormalVector.Scale(-1.);
											// calculate angle
											TempAngle = NormalVector.Angle(&VirtualNormalVector);
											*out << Verbose(4) << "NormalVector is ";
											VirtualNormalVector.Output(out);
											*out << " and the angle is " << TempAngle << "." << endl;
											if (SmallestAngle > TempAngle)
												{ // set to new possible winner
													SmallestAngle = TempAngle;
													winner = target;
												}
										}
								}
						// 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
						if (winner != PointsOnBoundary.end())
							{
								*out << Verbose(2) << "Winning target point is "
										<< *(winner->second) << " with angle " << SmallestAngle
										<< "." << endl;
								// create the lins of not yet present
								BLS[0] = baseline->second;
								// 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
								LineChecker[0] = baseline->second->endpoints[0]->lines.find(
										winner->first);
								LineChecker[1] = baseline->second->endpoints[1]->lines.find(
										winner->first);
								if (LineChecker[0]
										== baseline->second->endpoints[0]->lines.end())
									{ // create
										BPS[0] = baseline->second->endpoints[0];
										BPS[1] = winner->second;
										BLS[1] = new class BoundaryLineSet(BPS,
												LinesOnBoundaryCount);
										LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
												BLS[1]));
										LinesOnBoundaryCount++;
									}
								else
									BLS[1] = LineChecker[0]->second;
								if (LineChecker[1]
										== baseline->second->endpoints[1]->lines.end())
									{ // create
										BPS[0] = baseline->second->endpoints[1];
										BPS[1] = winner->second;
										BLS[2] = new class BoundaryLineSet(BPS,
												LinesOnBoundaryCount);
										LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
												BLS[2]));
										LinesOnBoundaryCount++;
									}
								else
									BLS[2] = LineChecker[1]->second;
								BTS = new class BoundaryTriangleSet(BLS,
										TrianglesOnBoundaryCount);
								TrianglesOnBoundary.insert(TrianglePair(
										TrianglesOnBoundaryCount, BTS));
								TrianglesOnBoundaryCount++;
							}
						else
							{
								*out << Verbose(1)
										<< "I could not determine a winner for this baseline "
										<< *(baseline->second) << "." << endl;
							}

						// 5d. If the set of lines is not yet empty, go to 5. and continue
					}
				else
					*out << Verbose(2) << "Baseline candidate " << *(baseline->second)
							<< " has a triangle count of "
							<< baseline->second->TrianglesCount << "." << endl;
		}
	while (flag);

}
;

/** Adds an atom to the tesselation::PointsOnBoundary list.
 * \param *Walker atom to add
 */
void
Tesselation::AddPoint(atom *Walker)
{
	PointTestPair InsertUnique;
	BPS[0] = new class BoundaryPointSet(Walker);
	InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
	if (InsertUnique.second) // if new point was not present before, increase counter
		PointsOnBoundaryCount++;
}
;

/** Adds point to Tesselation::PointsOnBoundary if not yet present.
 * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
 * @param Candidate point to add
 * @param n index for this point in Tesselation::TPS array
 */
void
Tesselation::AddTrianglePoint(atom* Candidate, int n)
{
	PointTestPair InsertUnique;
	TPS[n] = new class BoundaryPointSet(Candidate);
	InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
	if (InsertUnique.second) // if new point was not present before, increase counter
		{
			PointsOnBoundaryCount++;
		}
	else
		{
			delete TPS[n];
			cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node)
					<< " gibt's schon in der PointMap." << endl;
			TPS[n] = (InsertUnique.first)->second;
		}
}
;

/** Function tries to add line from current Points in BPS to BoundaryLineSet.
 * If succesful it raises the line count and inserts the new line into the BLS,
 * if unsuccesful, it writes the line which had been present into the BLS, deleting the new constructed one.
 * @param *a first endpoint
 * @param *b second endpoint
 * @param n index of Tesselation::BLS giving the line with both endpoints
 */
void
Tesselation::AddTriangleLine(class BoundaryPointSet *a,
		class BoundaryPointSet *b, int n)
{
	LineMap::iterator LineWalker;
	//cout << "Manually checking endpoints for line." << endl;
	if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr)
	//If a line is there, how do I recognize that beyond a shadow of a doubt?
		{
			//cout << Verbose(2) << "Line exists already, retrieving it from LinesOnBoundarySet" << endl;

			LineWalker = LinesOnBoundary.end();
			LineWalker--;

			while (LineWalker->second->endpoints[0]->node->nr != min(a->node->nr,
					b->node->nr) or LineWalker->second->endpoints[1]->node->nr != max(
					a->node->nr, b->node->nr))
				{
					//cout << Verbose(1) << "Looking for line which already exists"<< endl;
					LineWalker--;
				}
			BPS[0] = LineWalker->second->endpoints[0];
			BPS[1] = LineWalker->second->endpoints[1];
			BLS[n] = LineWalker->second;

		}
	else
		{
			cout << Verbose(2)
					<< "Adding line which has not been used before between "
					<< *(a->node) << " and " << *(b->node) << "." << endl;
			BPS[0] = a;
			BPS[1] = b;
			BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);

			LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
			LinesOnBoundaryCount++;

		}
}
;

/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
 * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
 */
void
Tesselation::AddTriangleToLines()
{

	cout << Verbose(1) << "Adding triangle to its lines" << endl;
	int i = 0;
	TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
	TrianglesOnBoundaryCount++;

	/*
	 * this is apparently done when constructing triangle

	 for (i=0; i<3; i++)
	 {
	 BLS[i]->AddTriangle(BTS);
	 }
	 */
}
;


double det_get(gsl_matrix *A, int inPlace) {
	/*
	inPlace = 1 => A is replaced with the LU decomposed copy.
	inPlace = 0 => A is retained, and a copy is used for LU.
	*/

	double det;
	int signum;
	gsl_permutation *p = gsl_permutation_alloc(A->size1);
	gsl_matrix *tmpA;

	if (inPlace)
	tmpA = A;
	else {
	gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
	gsl_matrix_memcpy(tmpA , A);
	}


	gsl_linalg_LU_decomp(tmpA , p , &signum);
	det = gsl_linalg_LU_det(tmpA , signum);
	gsl_permutation_free(p);
	if (! inPlace)
	gsl_matrix_free(tmpA);

	return det;
};

void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
{
	gsl_matrix *A = gsl_matrix_calloc(3,3);
	double m11, m12, m13, m14;

	for(int i=0;i<3;i++) {
		gsl_matrix_set(A, i, 0, a.x[i]);
		gsl_matrix_set(A, i, 1, b.x[i]);
		gsl_matrix_set(A, i, 2, c.x[i]);
	}
	m11 = det_get(A, 1);

	for(int i=0;i<3;i++) {
		gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
		gsl_matrix_set(A, i, 1, b.x[i]);
		gsl_matrix_set(A, i, 2, c.x[i]);
	}
	m12 = det_get(A, 1);

	for(int i=0;i<3;i++) {
		gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
		gsl_matrix_set(A, i, 1, a.x[i]);
		gsl_matrix_set(A, i, 2, c.x[i]);
	}
	m13 = det_get(A, 1);

	for(int i=0;i<3;i++) {
		gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
		gsl_matrix_set(A, i, 1, a.x[i]);
		gsl_matrix_set(A, i, 2, b.x[i]);
	}
	m14 = det_get(A, 1);

	if (fabs(m11) < MYEPSILON)
		cerr << "ERROR: three points are colinear." << endl;

	center->x[0] =	0.5 * m12/ m11;
	center->x[1] = -0.5 * m13/ m11;
	center->x[2] =	0.5 * m14/ m11;

	if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
		cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;

	gsl_matrix_free(A);
};



/**
 * Function returns center of sphere with RADIUS, which rests on points a, b, c
 * @param Center this vector will be used for return
 * @param a vector first point of triangle
 * @param b vector second point of triangle
 * @param c vector third point of triangle
 * @param Direction vector indicates up/down
 * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
 * @param Halfplaneindicator double indicates whether Direction is up or down
 * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
 * @param alpha double angle at a
 * @param beta double, angle at b
 * @param gamma, double, angle at c
 * @param Radius, double
 * @param Umkreisradius double radius of circumscribing circle
 */
void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
		double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
{
	Vector TempNormal, helper;
	double Restradius;
	Vector OtherCenter;
	cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
	Center->Zero();
	helper.CopyVector(&a);
	helper.Scale(sin(2.*alpha));
	Center->AddVector(&helper);
	helper.CopyVector(&b);
	helper.Scale(sin(2.*beta));
	Center->AddVector(&helper);
	helper.CopyVector(&c);
	helper.Scale(sin(2.*gamma));
	Center->AddVector(&helper);
	//*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
	Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
	NewUmkreismittelpunkt->CopyVector(Center);
	cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
	// Here we calculated center of circumscribing circle, using barycentric coordinates
	cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";

	TempNormal.CopyVector(&a);
	TempNormal.SubtractVector(&b);
	helper.CopyVector(&a);
	helper.SubtractVector(&c);
	TempNormal.VectorProduct(&helper);
	if (fabs(HalfplaneIndicator) < MYEPSILON)
		{
			if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
				{
					TempNormal.Scale(-1);
				}
		}
	else
		{
			if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
				{
					TempNormal.Scale(-1);
				}
		}

	TempNormal.Normalize();
	Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
	cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
	TempNormal.Scale(Restradius);
	cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";

	Center->AddVector(&TempNormal);
	cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
	get_sphere(&OtherCenter, a, b, c, RADIUS);
	cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
	cout << Verbose(3) << "End of Get_center_of_sphere.\n";
};

/** This recursive function finds a third point, to form a triangle with two given ones.
 * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
 *	supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
 *	upon which we operate.
 *	If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
 *	direction and angle into Storage.
 *	We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
 *	with all neighbours of the candidate.
 * @param a first point
 * @param b second point
 * *param c atom old third point of old triangle
 * @param Candidate base point along whose bonds to start looking from
 * @param Parent point to avoid during search as its wrong direction
 * @param RecursionLevel contains current recursion depth
 * @param Chord baseline vector of first and second point
 * @param direction1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
 * @param OldNormal normal of the triangle which the baseline belongs to
 * @param ReferencePoint Vector of center of circumscribing circle from which we look towards center of sphere
 * @param Opt_Candidate candidate reference to return
 * @param Storage array containing two angles of current Opt_Candidate
 * @param RADIUS radius of ball
 * @param mol molecule structure with atoms and bonds
 */
void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
		int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
		atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
{
	cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
	cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
	cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
	cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
	cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
	cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
	/* OldNormal is normal vector on the old triangle
	 * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
	 * Chord points from b to a!!!
	 */
	Vector dif_a; //Vector from a to candidate
	Vector dif_b; //Vector from b to candidate
	Vector AngleCheck;
	Vector TempNormal, Umkreismittelpunkt;
	Vector Mittelpunkt;

	double CurrentEpsilon = 0.1;
	double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
	double BallAngle, AlternativeSign;
	atom *Walker; // variable atom point

	Vector NewUmkreismittelpunkt;

	if (a != Candidate and b != Candidate and c != Candidate) {
		cout << Verbose(3) << "We have a unique candidate!" << endl;
		dif_a.CopyVector(&(a->x));
		dif_a.SubtractVector(&(Candidate->x));
		dif_b.CopyVector(&(b->x));
		dif_b.SubtractVector(&(Candidate->x));
		AngleCheck.CopyVector(&(Candidate->x));
		AngleCheck.SubtractVector(&(a->x));
		AngleCheck.ProjectOntoPlane(Chord);

		SideA = dif_b.Norm();
		SideB = dif_a.Norm();
		SideC = Chord->Norm();
		//Chord->Scale(-1);

		alpha = Chord->Angle(&dif_a);
		beta = M_PI - Chord->Angle(&dif_b);
		gamma = dif_a.Angle(&dif_b);

		cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;

		if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
			cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
			cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
		}

		if ((M_PI*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {
			Umkreisradius = SideA / 2.0 / sin(alpha);
			//cout << Umkreisradius << endl;
			//cout << SideB / 2.0 / sin(beta) << endl;
			//cout << SideC / 2.0 / sin(gamma) << endl;

			if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
				cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
				cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
				sign = AngleCheck.ScalarProduct(direction1);
				if (fabs(sign)<MYEPSILON) {
					if (AngleCheck.ScalarProduct(OldNormal)<0) {
						sign =0;
						AlternativeSign=1;
					} else {
						sign =0;
						AlternativeSign=-1;
					}
				} else {
					sign /= fabs(sign);
				}
				if (sign >= 0) {
					cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
					Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
					Mittelpunkt.SubtractVector(&ReferencePoint);
					cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;
					BallAngle = Mittelpunkt.Angle(OldNormal);
					cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;

					//cout << "direction1 is " << *direction1 << "." << endl;
					//cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl;
					//cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;

					NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);

					if ((Mittelpunkt.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
						if (Storage[0]< -1.5) { // first Candidate at all
							if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
								cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
								Opt_Candidate = Candidate;
								Storage[0] = sign;
								Storage[1] = AlternativeSign;
								Storage[2] = BallAngle;
								cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;
							} else
								cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
						} else {
							if ( Storage[2] > BallAngle) {
								if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
									cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
									Opt_Candidate = Candidate;
									Storage[0] = sign;
									Storage[1] = AlternativeSign;
									Storage[2] = BallAngle;
									cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;
								} else
									cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
							} else {
								if (DEBUG) {
									cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
								}
							}
						}
					} else {
						if (DEBUG) {
							cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
						}
					}
				} else {
					if (DEBUG) {
						cout << Verbose(3) << *Candidate << " is not in search direction." << endl;
					}
				}
			} else {
				if (DEBUG) {
					cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
				}
			}
		} else {
			if (DEBUG) {
				cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
			}
		}
	} else {
		if (DEBUG) {
			cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
		}
	}

	if (RecursionLevel < 5) { // Seven is the recursion level threshold.
		for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
			Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
			if (Walker == Parent) { // don't go back the same bond
				continue;
			} else {
				Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
			}
		}
	}
	cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
}
;


/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
 * \param *Center new center on return
 * \param *a first point
 * \param *b second point
 * \param *c third point
 */
void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
{
	Vector helper;
	double alpha, beta, gamma;
	Vector SideA, SideB, SideC;
	SideA.CopyVector(b);
	SideA.SubtractVector(c);
	SideB.CopyVector(c);
	SideB.SubtractVector(a);
	SideC.CopyVector(a);
	SideC.SubtractVector(b);
	alpha = M_PI - SideB.Angle(&SideC);
	beta = M_PI - SideC.Angle(&SideA);
	gamma = M_PI - SideA.Angle(&SideB);
	cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
	if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
		cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;

	Center->Zero();
	helper.CopyVector(a);
	helper.Scale(sin(2.*alpha));
	Center->AddVector(&helper);
	helper.CopyVector(b);
	helper.Scale(sin(2.*beta));
	Center->AddVector(&helper);
	helper.CopyVector(c);
	helper.Scale(sin(2.*gamma));
	Center->AddVector(&helper);
	Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
};

/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
 * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
 * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
 * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
 * \param CircleCenter Center of the parameter circle
 * \param CirclePlaneNormal normal vector to plane of the parameter circle
 * \param CircleRadius radius of the parameter circle
 * \param NewSphereCenter new center of a circumcircle
 * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
 * \param NormalVector normal vector
 * \param SearchDirection search direction to make angle unique on return.
 * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
 */
double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
{
	Vector helper;
	double radius, alpha;

	helper.CopyVector(&NewSphereCenter);
	// test whether new center is on the parameter circle's plane
	if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
		cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))	<< "!" << endl;
		helper.ProjectOntoPlane(&CirclePlaneNormal);
	}
	radius = helper.ScalarProduct(&helper);
	// test whether the new center vector has length of CircleRadius
	if (fabs(radius - CircleRadius) > HULLEPSILON)
		cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
	alpha = helper.Angle(&OldSphereCenter);
	// make the angle unique by checking the halfplanes/search direction
	if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)	// acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
		alpha = 2.*M_PI - alpha;
	cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
	radius = helper.Distance(&OldSphereCenter);
	helper.ProjectOntoPlane(&NormalVector);
	// check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
	if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
		cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
		return alpha;
	} else {
		cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
		return 2.*M_PI;
	}
};


/** This recursive function finds a third point, to form a triangle with two given ones.
 * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
 * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
 * the center of the sphere is still fixed up to a single parameter. The band of possible values
 * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
 * us the "null" on this circle, the new center of the candidate point will be some way along this
 * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
 * by the normal vector of the base triangle that always points outwards by construction.
 * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
 * We construct the normal vector that defines the plane this circle lies in, it is just in the
 * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
 * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
 * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
 * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
 * both.
 * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
 * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
 * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
 * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
 * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
 * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
 * @param BaseTriangle BoundaryTriangleSet of the current base triangle with all three points
 * @param BaseLine BoundaryLineSet of BaseTriangle with the current base line
 * @param OptCandidate candidate reference on return
 * @param OptCandidateCenter candidate's sphere center on return
 * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
 * @param RADIUS radius of sphere
 * @param *LC LinkedCell structure with neighbouring atoms
 */
// void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
// {
//	 atom *Walker = NULL;
//	 Vector CircleCenter;	// center of the circle, i.e. of the band of sphere's centers
//	 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
//	 Vector OldSphereCenter;	 // center of the sphere defined by the three points of BaseTriangle
//	 Vector NewSphereCenter;	 // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
//	 Vector OtherNewSphereCenter;	 // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
//	 Vector NewNormalVector;	 // normal vector of the Candidate's triangle
//	 Vector SearchDirection;	 // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)
//	 Vector helper;
//	 LinkedAtoms *List = NULL;
//	 double CircleRadius; // radius of this circle
//	 double radius;
//	 double alpha, Otheralpha; // angles (i.e. parameter for the circle).
//	 double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
//	 int N[NDIM], Nlower[NDIM], Nupper[NDIM];
//	 atom *Candidate = NULL;
//
//	 cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;
//
//	 cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;
//
//	 // construct center of circle
//	 CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
//	 CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
//	 CircleCenter.Scale(0.5);
//
//	 // construct normal vector of circle
//	 CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
//	 CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
//
//	 // calculate squared radius of circle
//	 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
//	 if (radius/4. < RADIUS*RADIUS) {
//		 CircleRadius = RADIUS*RADIUS - radius/4.;
//		 CirclePlaneNormal.Normalize();
//		 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
//
//		 // construct old center
//		 GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));
//		 helper.CopyVector(&BaseTriangle->NormalVector);	// normal vector ensures that this is correct center of the two possible ones
//		 radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
//		 helper.Scale(sqrt(RADIUS*RADIUS - radius));
//		 OldSphereCenter.AddVector(&helper);
//		 OldSphereCenter.SubtractVector(&CircleCenter);
//		 cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
//
//		 // test whether old center is on the band's plane
//		 if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
//			 cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
//			 OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
//		 }
//		 radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
//		 if (fabs(radius - CircleRadius) < HULLEPSILON) {
//
//			 // construct SearchDirection
//			 SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);
//			 helper.CopyVector(&BaseLine->endpoints[0]->node->x);
//			 for(int i=0;i<3;i++)	// just take next different endpoint
//				 if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {
//					 helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);
//				 }
//			 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)	// ohoh, SearchDirection points inwards!
//				 SearchDirection.Scale(-1.);
//			 SearchDirection.ProjectOntoPlane(&OldSphereCenter);
//			 SearchDirection.Normalize();
//			 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
//			 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {	// rotated the wrong way!
//				 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
//			 }
//
//			 if (LC->SetIndexToVector(&CircleCenter)) {	// get cell for the starting atom
//				 for(int i=0;i<NDIM;i++) // store indices of this cell
//					 N[i] = LC->n[i];
//				 cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
//			 } else {
//				 cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
//				 return;
//			 }
//			 // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
//			 cout << Verbose(2) << "LC Intervals:";
//			 for (int i=0;i<NDIM;i++) {
//				 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
//				 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
//				 cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
//			 }
//			 cout << endl;
//			 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
//				 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
//					 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
//						 List = LC->GetCurrentCell();
//						 cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
//						 if (List != NULL) {
//							 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
//								 Candidate = (*Runner);
//
//								 // check for three unique points
//								 if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {
//									 cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
//
//									 // construct both new centers
//									 GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
//									 OtherNewSphereCenter.CopyVector(&NewSphereCenter);
//
//									 if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
//										 helper.CopyVector(&NewNormalVector);
//										 cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
//										 radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
//										 if (radius < RADIUS*RADIUS) {
//											 helper.Scale(sqrt(RADIUS*RADIUS - radius));
//											 cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
//											 NewSphereCenter.AddVector(&helper);
//											 NewSphereCenter.SubtractVector(&CircleCenter);
//											 cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
//
//											 helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
//											 OtherNewSphereCenter.AddVector(&helper);
//											 OtherNewSphereCenter.SubtractVector(&CircleCenter);
//											 cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
//
//											 // check both possible centers
//											 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
//											 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
//											 alpha = min(alpha, Otheralpha);
//											 if (*ShortestAngle > alpha) {
//													 OptCandidate = Candidate;
//													 *ShortestAngle = alpha;
//													 if (alpha != Otheralpha)
//														 OptCandidateCenter->CopyVector(&NewSphereCenter);
//													 else
//														 OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
//													 cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
//											 } else {
//												 if (OptCandidate != NULL)
//													 cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
//												 else
//													 cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
//											 }
//
//										 } else {
//											 cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
//										 }
//									 } else {
//										 cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
//									 }
//								 } else {
//									 cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;
//								 }
//							 }
//						 }
//					 }
//		 } else {
//			 cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
//		 }
//	 } else {
//		 cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;
//	 }
//
//	 cout << Verbose(1) << "End of Find_next_suitable_point" << endl;
// };


/** Checks whether the triangle consisting of the three atoms is already present.
 * Searches for the points in Tesselation::PointsOnBoundary and checks their
 * lines. If any of the three edges already has two triangles attached, false is
 * returned.
 * \param *out output stream for debugging
 * \param *Candidates endpoints of the triangle candidate
 * \return false - triangle invalid due to edge criteria, true - triangle may be added.
 */
bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
	LineMap::iterator FindLine;
	PointMap::iterator FindPoint;
	bool Present[3];

	*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
	for (int i=0;i<3;i++) { // check through all endpoints
		FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
		if (FindPoint != PointsOnBoundary.end())
			TPS[i] = FindPoint->second;
		else
			TPS[i] = NULL;
	}

	// check lines
	for (int i=0;i<3;i++)
		if (TPS[i] != NULL)
			for (int j=i;j<3;j++)
				if (TPS[j] != NULL) {
					FindLine = TPS[i]->lines.find(TPS[j]->node->nr);
					if ((FindLine != TPS[i]->lines.end()) && (FindLine->second->TrianglesCount > 1)) {
						*out << "WARNING: Line " << *FindLine->second << " already present with " << FindLine->second->TrianglesCount << " triangles attached." << endl;
						*out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
						return false;
					}
				}
	*out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
	return true;
};

/** This recursive function finds a third point, to form a triangle with two given ones.
 * Note that this function is for the starting triangle.
 * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
 * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
 * the center of the sphere is still fixed up to a single parameter. The band of possible values
 * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
 * us the "null" on this circle, the new center of the candidate point will be some way along this
 * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
 * by the normal vector of the base triangle that always points outwards by construction.
 * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
 * We construct the normal vector that defines the plane this circle lies in, it is just in the
 * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
 * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
 * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
 * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
 * both.
 * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
 * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
 * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
 * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
 * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
 * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
 * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
 * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
 * @param BaseLine BoundaryLineSet with the current base line
 * @param ThirdNode third atom to avoid in search
 * @param OptCandidate candidate reference on return
 * @param OptCandidateCenter candidate's sphere center on return
 * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
 * @param RADIUS radius of sphere
 * @param *LC LinkedCell structure with neighbouring atoms
 */
void Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, atom *ThirdNode, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
{
	atom *Walker = NULL;
	Vector CircleCenter;	// center of the circle, i.e. of the band of sphere's centers
	Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
	Vector NewSphereCenter;	 // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
	Vector OtherNewSphereCenter;	 // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
	Vector NewNormalVector;	 // normal vector of the Candidate's triangle
	Vector helper;
	LinkedAtoms *List = NULL;
	double CircleRadius; // radius of this circle
	double radius;
	double alpha, Otheralpha; // angles (i.e. parameter for the circle).
	double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
	int N[NDIM], Nlower[NDIM], Nupper[NDIM];
	atom *Candidate = NULL;

	cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;

	cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;

	// construct center of circle
	CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
	CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
	CircleCenter.Scale(0.5);

	// construct normal vector of circle
	CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
	CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);

	// calculate squared radius of circle
	radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
	if (radius/4. < RADIUS*RADIUS) {
		CircleRadius = RADIUS*RADIUS - radius/4.;
		CirclePlaneNormal.Normalize();
		cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;

		// test whether old center is on the band's plane
		if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
			cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
			OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
		}
		radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
		if (fabs(radius - CircleRadius) < HULLEPSILON) {

			// check SearchDirection
			cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
			if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {	// rotated the wrong way!
				cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
			}
			// get cell for the starting atom
			if (LC->SetIndexToVector(&CircleCenter)) {
					for(int i=0;i<NDIM;i++) // store indices of this cell
					N[i] = LC->n[i];
				cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
			} else {
				cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
				return;
			}
			// then go through the current and all neighbouring cells and check the contained atoms for possible candidates
			cout << Verbose(2) << "LC Intervals:";
			for (int i=0;i<NDIM;i++) {
				Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
				Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
				cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
			}
			cout << endl;
			for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
				for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
					for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
						List = LC->GetCurrentCell();
						//cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
						if (List != NULL) {
							for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
								Candidate = (*Runner);

								// check for three unique points
								if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) {
									cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;

									// construct both new centers
									GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
									OtherNewSphereCenter.CopyVector(&NewSphereCenter);

									if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
										helper.CopyVector(&NewNormalVector);
										cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
										radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
										if (radius < RADIUS*RADIUS) {
											helper.Scale(sqrt(RADIUS*RADIUS - radius));
											cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
											NewSphereCenter.AddVector(&helper);
											NewSphereCenter.SubtractVector(&CircleCenter);
											cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;

											helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
											OtherNewSphereCenter.AddVector(&helper);
											OtherNewSphereCenter.SubtractVector(&CircleCenter);
											cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;

											// check both possible centers
											alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
											Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
											alpha = min(alpha, Otheralpha);
											if (*ShortestAngle > alpha) {
													OptCandidate = Candidate;
													*ShortestAngle = alpha;
													if (alpha != Otheralpha)
														OptCandidateCenter->CopyVector(&NewSphereCenter);
													else
														OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
													cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
											} else {
												if (OptCandidate != NULL)
													cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
												else
													cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
											}

										} else {
											cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
										}
									} else {
										cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
									}
								} else {
									if (ThirdNode != NULL)
										cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
									else
										cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
								}
							}
						}
					}
		} else {
			cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
		}
	} else {
		if (ThirdNode != NULL)
			cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
		else
			cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
	}

	cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
};

/** Finds the second point of starting triangle.
 * \param *a first atom
 * \param *Candidate pointer to candidate atom on return
 * \param Oben vector indicating the outside
 * \param Opt_Candidate reference to recommended candidate on return
 * \param Storage[3] array storing angles and other candidate information
 * \param RADIUS radius of virtual sphere
 * \param *LC LinkedCell structure with neighbouring atoms
 */
void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
{
	cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
	int i;
	Vector AngleCheck;
	atom* Walker;
	double norm = -1., angle;
	LinkedAtoms *List = NULL;
	int N[NDIM], Nlower[NDIM], Nupper[NDIM];

	if (LC->SetIndexToAtom(a)) {	// get cell for the starting atom
		for(int i=0;i<NDIM;i++) // store indices of this cell
			N[i] = LC->n[i];
	} else {
		cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
		return;
	}
	// then go through the current and all neighbouring cells and check the contained atoms for possible candidates
	cout << Verbose(2) << "LC Intervals:";
	for (int i=0;i<NDIM;i++) {
		Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
		Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
		cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
	}
	cout << endl;
	for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
		for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
			for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
				List = LC->GetCurrentCell();
				//cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
				if (List != NULL) {
					for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
						Candidate = (*Runner);
								// check if we only have one unique point yet ...
								if (a != Candidate) {
			cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
			AngleCheck.CopyVector(&(Candidate->x));
			AngleCheck.SubtractVector(&(a->x));
			norm = AngleCheck.Norm();
							// second point shall have smallest angle with respect to Oben vector
							if (norm < RADIUS) {
								angle = AngleCheck.Angle(&Oben);
								if (angle < Storage[0]) {
									//cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
									cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
									Opt_Candidate = Candidate;
									Storage[0] = AngleCheck.Angle(&Oben);
									//cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
								} else {
									cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
								}
							} else {
								cout << "Refused due to Radius " << norm << endl;
							}
						}
					}
				}
			}
	cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
};

/** Finds the starting triangle for find_non_convex_border().
 * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
 * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
 * point are called.
 * \param RADIUS radius of virtual rolling sphere
 * \param *LC LinkedCell structure with neighbouring atoms
 */
void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
{
	cout << Verbose(1) << "Begin of Find_starting_triangle\n";
	int i = 0;
	LinkedAtoms *List = NULL;
	atom* Walker;
	atom* FirstPoint;
	atom* SecondPoint;
	atom* MaxAtom[NDIM];
	double max_coordinate[NDIM];
	Vector Oben;
	Vector helper;
	Vector Chord;
	Vector SearchDirection;
	Vector OptCandidateCenter;

	Oben.Zero();

	for (i = 0; i < 3; i++) {
		MaxAtom[i] = NULL;
		max_coordinate[i] = -1;
	}

	// 1. searching topmost atom with respect to each axis
	for (int i=0;i<NDIM;i++) { // each axis
		LC->n[i] = LC->N[i]-1; // current axis is topmost cell
		for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
			for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
				List = LC->GetCurrentCell();
				//cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
				if (List != NULL) {
					for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
						cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl;
						if ((*Runner)->x.x[i] > max_coordinate[i]) {
							max_coordinate[i] = (*Runner)->x.x[i];
							MaxAtom[i] = (*Runner);
						}
					}
				} else {
					cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
				}
			}
	}

	cout << Verbose(2) << "Found maximum coordinates: ";
	for (int i=0;i<NDIM;i++)
		cout << i << ": " << *MaxAtom[i] << "\t";
	cout << endl;
	const int k = 1;		// arbitrary choice
	Oben.x[k] = 1.;
	FirstPoint = MaxAtom[k];
	cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;

	// add first point
	AddTrianglePoint(FirstPoint, 0);

	double ShortestAngle;
	atom* Opt_Candidate = NULL;
	ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.

	Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
	SecondPoint = Opt_Candidate;
	cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";

	// add second point and first baseline
	AddTrianglePoint(SecondPoint, 1);
	AddTriangleLine(TPS[0], TPS[1], 0);

	helper.CopyVector(&(FirstPoint->x));
	helper.SubtractVector(&(SecondPoint->x));
	helper.Normalize();
	Oben.ProjectOntoPlane(&helper);
	Oben.Normalize();
	helper.VectorProduct(&Oben);
	ShortestAngle = 2.*M_PI; // This will indicate the quadrant.

	Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
	Chord.SubtractVector(&(SecondPoint->x));
	double radius = Chord.ScalarProduct(&Chord);
	double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
	helper.CopyVector(&Oben);
	helper.Scale(CircleRadius);
	// Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)

	cout << Verbose(2) << "Looking for third point candidates \n";
	// look in one direction of baseline for initial candidate
	Opt_Candidate = NULL;
	SearchDirection.MakeNormalVector(&Chord, &Oben);	// whether we look "left" first or "right" first is not important ...

	cout << Verbose(1) << "Looking for third point candidates ...\n";
	Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
	cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;

	// add third point
	AddTrianglePoint(Opt_Candidate, 2);

	// FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate

	// Finally, we only have to add the found further lines
	AddTriangleLine(TPS[1], TPS[2], 1);
	AddTriangleLine(TPS[0], TPS[2], 2);
	// ... and triangles to the Maps of the Tesselation class
	BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
	AddTriangleToLines();
	// ... and calculate its normal vector (with correct orientation)
	OptCandidateCenter.Scale(-1.);
	cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl;
	BTS->GetNormalVector(OptCandidateCenter);
	cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
	cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
	cout << Verbose(1) << "End of Find_starting_triangle\n";
};

/** This function finds a triangle to a line, adjacent to an existing one.
 * @param out output stream for debugging
 * @param *mol molecule with Atom's and Bond's
 * @param Line current baseline to search from
 * @param T current triangle which \a Line is edge of
 * @param RADIUS radius of the rolling ball
 * @param N number of found triangles
 * @param *filename filename base for intermediate envelopes
 * @param *LC LinkedCell structure with neighbouring atoms
 */
bool Tesselation::Find_next_suitable_triangle(ofstream *out,
		molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
		const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
{
	cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
	ofstream *tempstream = NULL;
	char NumberName[255];
	double tmp;

	atom* Opt_Candidate = NULL;
	Vector OptCandidateCenter;

	Vector CircleCenter;
	Vector CirclePlaneNormal;
	Vector OldSphereCenter;
	Vector SearchDirection;
	Vector helper;
	atom *ThirdNode = NULL;
	double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
	double radius, CircleRadius;

	cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
	for (int i=0;i<3;i++)
		if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
			ThirdNode = T.endpoints[i]->node;

	// construct center of circle
	CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
	CircleCenter.AddVector(&Line.endpoints[1]->node->x);
	CircleCenter.Scale(0.5);

	// construct normal vector of circle
	CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
	CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);

	// calculate squared radius of circle
	radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
	if (radius/4. < RADIUS*RADIUS) {
		CircleRadius = RADIUS*RADIUS - radius/4.;
		CirclePlaneNormal.Normalize();
		cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;

		// construct old center
		GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
		helper.CopyVector(&T.NormalVector);	// normal vector ensures that this is correct center of the two possible ones
		radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
		helper.Scale(sqrt(RADIUS*RADIUS - radius));
		OldSphereCenter.AddVector(&helper);
		OldSphereCenter.SubtractVector(&CircleCenter);
		cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;

		// construct SearchDirection
		SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
		helper.CopyVector(&Line.endpoints[0]->node->x);
		helper.SubtractVector(&ThirdNode->x);
		if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)	// ohoh, SearchDirection points inwards!
			SearchDirection.Scale(-1.);
		SearchDirection.ProjectOntoPlane(&OldSphereCenter);
		SearchDirection.Normalize();
		cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
		if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {	// rotated the wrong way!
			cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
		}

		// add third point
		cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
		Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);

	} else {
		cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
	}

	if (Opt_Candidate == NULL) {
		cerr << "WARNING: Could not find a suitable candidate." << endl;
		return false;
	}
	cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << endl;

	// check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
	atom *AtomCandidates[3];
	AtomCandidates[0] = Opt_Candidate;
	AtomCandidates[1] = Line.endpoints[0]->node;
	AtomCandidates[2] = Line.endpoints[1]->node;
	bool flag = CheckPresenceOfTriangle(out, AtomCandidates);

	if (flag) { // if so, add
		AddTrianglePoint(Opt_Candidate, 0);
		AddTrianglePoint(Line.endpoints[0]->node, 1);
		AddTrianglePoint(Line.endpoints[1]->node, 2);

		AddTriangleLine(TPS[0], TPS[1], 0);
		AddTriangleLine(TPS[0], TPS[2], 1);
		AddTriangleLine(TPS[1], TPS[2], 2);

		BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
		AddTriangleToLines();

		OptCandidateCenter.Scale(-1.);
		BTS->GetNormalVector(OptCandidateCenter);
		OptCandidateCenter.Scale(-1.);

		cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
		cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
	} else { // else, yell and do nothing
		cout << Verbose(1) << "This triangle consisting of ";
		cout << *Opt_Candidate << ", ";
		cout << *Line.endpoints[0]->node << " and ";
		cout << *Line.endpoints[1]->node << " ";
		cout << "is invalid!" << endl;
		return false;
	}

	if (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration
		sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
		if (DoTecplotOutput) {
			string NameofTempFile(tempbasename);
			NameofTempFile.append(NumberName);
			for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos))
				NameofTempFile.erase(npos, 1);
			NameofTempFile.append(TecplotSuffix);
			cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
			tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
			write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
			tempstream->close();
			tempstream->flush();
			delete(tempstream);
		}

		if (DoRaster3DOutput) {
			string NameofTempFile(tempbasename);
			NameofTempFile.append(NumberName);
			for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos))
				NameofTempFile.erase(npos, 1);
			NameofTempFile.append(Raster3DSuffix);
			cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
			tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
			write_raster3d_file(out, tempstream, this, mol);
			// include the current position of the virtual sphere in the temporary raster3d file
			// make the circumsphere's center absolute again
			helper.CopyVector(&Line.endpoints[0]->node->x);
			helper.AddVector(&Line.endpoints[1]->node->x);
			helper.Scale(0.5);
			OptCandidateCenter.AddVector(&helper);
			Vector *center = mol->DetermineCenterOfAll(out);
			OptCandidateCenter.AddVector(center);
			delete(center);
			// and add to file plus translucency object
			*tempstream << "# current virtual sphere\n";
			*tempstream << "8\n	25.0		0.6		 -1.0 -1.0 -1.0		 0.2				0 0 0 0\n";
			*tempstream << "2\n	" << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n";
			*tempstream << "9\n	terminating special property\n";
			tempstream->close();
			tempstream->flush();
			delete(tempstream);
		}
		if (DoTecplotOutput || DoRaster3DOutput)
			TriangleFilesWritten++;
	}

	cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
	return true;
};

/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
 * \param *out output stream for debugging
 * \param *mol molecule structure with Atom's and Bond's
 * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
 * \param *LCList linked cell list of all atoms
 * \param *filename filename prefix for output of vertex data
 * \para RADIUS radius of the virtual sphere
 */
void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
{
	int N = 0;
	bool freeTess = false;
	*out << Verbose(1) << "Entering search for non convex hull. " << endl;
	if (Tess == NULL) {
		*out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
		Tess = new Tesselation;
		freeTess = true;
	}
	bool freeLC = false;
	LineMap::iterator baseline;
	*out << Verbose(0) << "Begin of Find_non_convex_border\n";
	bool flag = false;	// marks whether we went once through all baselines without finding any without two triangles
	bool failflag = false;

	if (LCList == NULL) {
		LCList = new LinkedCell(mol, 2.*RADIUS);
		freeLC = true;
	}

	Tess->Find_starting_triangle(out, mol, RADIUS, LCList);

	baseline = Tess->LinesOnBoundary.begin();
	while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
		if (baseline->second->TrianglesCount == 1) {
			failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
			flag = flag || failflag;
			if (!failflag)
				cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
		} else {
			cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
		}
		N++;
		baseline++;
		if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
			baseline = Tess->LinesOnBoundary.begin();	 // restart if we reach end due to newly inserted lines
			flag = false;
		}
	}
	if (1) { //failflag) {
		*out << Verbose(1) << "Writing final tecplot file\n";
		if (DoTecplotOutput) {
			string OutputName(filename);
			OutputName.append(TecplotSuffix);
			ofstream *tecplot = new ofstream(OutputName.c_str());
			write_tecplot_file(out, tecplot, Tess, mol, -1);
			tecplot->close();
			delete(tecplot);
		}
		if (DoRaster3DOutput) {
			string OutputName(filename);
			OutputName.append(Raster3DSuffix);
			ofstream *raster = new ofstream(OutputName.c_str());
			write_raster3d_file(out, raster, Tess, mol);
			raster->close();
			delete(raster);
		}
	} else {
		cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
	}
	if (freeTess)
		delete(Tess);
	if (freeLC)
		delete(LCList);
	*out << Verbose(0) << "End of Find_non_convex_border\n";
};