Changeset 776b64 for src

Timestamp:

Oct 27, 2009, 4:11:22 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

fb73b8

Parents:

ad37ab

Message:

Huge refactoring to make const what is const (ticket #38), continued.

• too many changes because of too many cross-references to be able to list them up here.
• NOTE that "make check" runs fine and did catch several error.
• note that we had to use const_iterator several times when the map, ... was declared const.
• at times we changed an allocated LinkedCell LCList(...) into

const LinkedCell *LCList;
LCList = new LinkedCell(...);

• also mutable (see ticket #5) was used, e.g. for molecule::InternalPointer (PointCloud changes are allowed, because they are just accounting).

Signed-off-by: Frederik Heber <heber@…>

Location:

src

Files:

• analysis_correlation.cpp (modified) (9 diffs)
• analysis_correlation.hpp (modified) (3 diffs)
• boundary.cpp (modified) (45 diffs)
• boundary.hpp (modified) (1 diff)
• builder.cpp (modified) (8 diffs)
• ellipsoid.cpp (modified) (5 diffs)
• leastsquaremin.cpp (modified) (1 diff)
• leastsquaremin.hpp (modified) (1 diff)
• linkedcell.cpp (modified) (9 diffs)
• linkedcell.hpp (modified) (1 diff)
• molecule.cpp (modified) (2 diffs)
• molecule.hpp (modified) (4 diffs)
• molecule_geometry.cpp (modified) (2 diffs)
• molecule_graph.cpp (modified) (3 diffs)
• molecule_pointcloud.cpp (modified) (9 diffs)
• moleculelist.cpp (modified) (3 diffs)
• tesselation.cpp (modified) (52 diffs)
• tesselation.hpp (modified) (10 diffs)
• tesselationhelpers.cpp (modified) (20 diffs)
• tesselationhelpers.hpp (modified) (1 diff)
• unittests/ActOnAllUnitTest.cpp (modified) (2 diffs)
• unittests/AnalysisCorrelationToPointUnitTest.cpp (modified) (1 diff)
• unittests/AnalysisCorrelationToSurfaceUnitTest.cpp (modified) (2 diffs)
• vector.cpp (modified) (43 diffs)
• vector.hpp (modified) (1 diff)

src/analysis_correlation.cpp

@@ -24 +24 @@
  * \return Map of doubles with values the pair of the two atoms.
  */
-PairCorrelationMap *PairCorrelation( ofstream *out, molecule *mol, element *type1, element *type2 )
+PairCorrelationMap *PairCorrelation( const ofstream *out, const molecule * const mol, const element * const type1, const element * const type2 )
 {
   PairCorrelationMap *outmap = NULL;
@@ -61 +61 @@
  * \return Map of dobules with values as pairs of atom and the vector
  */
-CorrelationToPointMap *CorrelationToPoint( ofstream *out, molecule *mol, element *type, Vector *point )
+CorrelationToPointMap *CorrelationToPoint( const ofstream *out, const molecule * const mol, const element * const type, const Vector *point )
 {
   CorrelationToPointMap *outmap = NULL;
@@ -76 +76 @@
     if ((type == NULL) || (Walker->type == type)) {
       distance = Walker->node->Distance(point);
-      outmap->insert ( pair<double, pair<atom *, Vector*> >(distance, pair<atom *, Vector*> (Walker, point) ) );
+      outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
     }
   }
@@ -91 +91 @@
  * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
  */
-CorrelationToSurfaceMap *CorrelationToSurface( ofstream *out, molecule *mol, element *type, Tesselation *Surface, LinkedCell *LC )
+CorrelationToSurfaceMap *CorrelationToSurface( ofstream *out, const molecule * const mol, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
 {
 
@@ -124 +124 @@
  * \param BinStart first bin
  */
-double GetBin ( double value, double BinWidth, double BinStart )
+double GetBin ( const double value, const double BinWidth, const double BinStart )
 {
   double bin =(double) (floor((value - BinStart)/BinWidth));
@@ -135 +135 @@
  * \param *map map to write
  */
-void OutputCorrelation( ofstream *file, BinPairMap *map )
+void OutputCorrelation( ofstream *file, const BinPairMap * const map )
 {
   *file << "# BinStart\tCount" << endl;
-  for (BinPairMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
+  for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << runner->second << endl;
   }
@@ -147 +147 @@
  * \param *map map to write
  */
-void OutputPairCorrelation( ofstream *file, PairCorrelationMap *map )
+void OutputPairCorrelation( ofstream *file, const PairCorrelationMap * const map )
 {
   *file << "# BinStart\tAtom1\tAtom2" << endl;
-  for (PairCorrelationMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
+  for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
   }
@@ -159 +159 @@
  * \param *map map to write
  */
-void OutputCorrelationToPoint( ofstream *file, CorrelationToPointMap *map )
+void OutputCorrelationToPoint( ofstream *file, const CorrelationToPointMap * const map )
 {
   *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
-  for (CorrelationToPointMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
+  for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first;
     for (int i=0;i<NDIM;i++)
@@ -174 +174 @@
  * \param *map map to write
  */
-void OutputCorrelationToSurface( ofstream *file, CorrelationToSurfaceMap *map )
+void OutputCorrelationToSurface( ofstream *file, const CorrelationToSurfaceMap * const map )
 {
   *file << "# BinStart\tTriangle" << endl;
-  for (CorrelationToSurfaceMap::iterator runner = map->begin(); runner != map->end(); ++runner) {
+  for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first << "\t" << *(runner->second.second) << endl;
   }

src/analysis_correlation.hpp

@@ -39 +39 @@
 
 typedef multimap<double, pair<atom *, atom *> > PairCorrelationMap;
-typedef multimap<double, pair<atom *, Vector *> > CorrelationToPointMap;
+typedef multimap<double, pair<atom *, const Vector *> > CorrelationToPointMap;
 typedef multimap<double, pair<atom *, BoundaryTriangleSet *> > CorrelationToSurfaceMap;
 typedef map<double, int> BinPairMap;
@@ -45 +45 @@
 /********************************************** declarations *******************************/
 
-PairCorrelationMap *PairCorrelation( ofstream *out, molecule *mol, element *type1, element *type2 );
-CorrelationToPointMap *CorrelationToPoint( ofstream *out, molecule *mol, element *type, Vector *point );
-CorrelationToSurfaceMap *CorrelationToSurface( ofstream *out, molecule *mol, element *type, Tesselation *Surface, LinkedCell *LC );
-double GetBin ( double value, double BinWidth, double BinStart );
-void OutputCorrelation( ofstream *file, BinPairMap *map );
-void OutputPairCorrelation( ofstream *file, BinPairMap *map );
-void OutputCorrelationToPoint( ofstream *file, BinPairMap *map );
-void OutputCorrelationToSurface( ofstream *file, BinPairMap *map );
+PairCorrelationMap *PairCorrelation( const ofstream *out, const molecule * const mol, const element * const type1, const element * const  type2 );
+CorrelationToPointMap *CorrelationToPoint( const ofstream *out, const molecule * const mol, const element * const type, const Vector *point );
+CorrelationToSurfaceMap *CorrelationToSurface( ofstream *out, const molecule * const mol, const element * const type, const Tesselation * const Surface, const LinkedCell *LC );
+double GetBin ( const double value, const double BinWidth, const double BinStart );
+void OutputCorrelation( ofstream *file, const BinPairMap * const map );
+void OutputPairCorrelation( ofstream *file, const BinPairMap * const map );
+void OutputCorrelationToPoint( ofstream *file, const BinPairMap * const map );
+void OutputCorrelationToSurface( ofstream *file, const BinPairMap * const map );
 
 
@@ -95 +95 @@
  * \return Map of doubles (the bins) with counts as values
  */
-template <typename T> BinPairMap *BinData ( ofstream *out,  T *map, double BinWidth, double BinStart, double BinEnd )
+template <typename T> BinPairMap *BinData ( ofstream *out,  T *map, const double BinWidth, const double BinStart, const double BinEnd )
 {
   BinPairMap *outmap = new BinPairMap;

src/boundary.cpp

@@ -26 +26 @@
  * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
  * \param *mol molecule structure representing the cluster
+ * \param *&TesselStruct Tesselation structure with triangles
  * \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)
+double *GetDiametersOfCluster(ofstream *out, const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
 {
   // get points on boundary of NULL was given as parameter
   bool BoundaryFreeFlag = false;
-  Boundaries *BoundaryPoints = BoundaryPtr;
   double OldComponent = 0.;
   double tmp = 0.;
@@ -42 +42 @@
   int component = 0;
   int Othercomponent = 0;
-  Boundaries::iterator Neighbour;
-  Boundaries::iterator OtherNeighbour;
+  Boundaries::const_iterator Neighbour;
+  Boundaries::const_iterator OtherNeighbour;
   double *GreatestDiameter = new double[NDIM];
 
-  if (BoundaryPoints == NULL) {
+  const Boundaries *BoundaryPoints;
+  if (BoundaryPtr == NULL) {
     BoundaryFreeFlag = true;
-    BoundaryPoints = GetBoundaryPoints(out, mol);
+    BoundaryPoints = GetBoundaryPoints(out, mol, TesselStruct);
   } else {
+    BoundaryPoints = BoundaryPtr;
     *out << Verbose(1) << "Using given boundary points set." << endl;
   }
@@ -63 +65 @@
           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++)
-            {
+          for (Boundaries::const_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
@@ -74 +74 @@
               DistanceVector.CopyVector(&runner->second.second->x);
               DistanceVector.SubtractVector(&Neighbour->second.second->x);
-              do
-                { // seek for neighbour pair where it flips
+              do { // seek for neighbour pair where it flips
                   OldComponent = DistanceVector.x[Othercomponent];
                   Neighbour++;
@@ -83 +82 @@
                   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(
+                } 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)
-                {
+              if (runner != Neighbour) {
                   OtherNeighbour = Neighbour;
                   if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
@@ -133 +130 @@
  * \param *out output stream for debugging
  * \param *mol molecule structure representing the cluster
- */
-Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol)
+ * \param *&TesselStruct pointer to Tesselation structure
+ */
+Boundaries *GetBoundaryPoints(ofstream *out, const molecule *mol, Tesselation *&TesselStruct)
 {
   atom *Walker = NULL;
@@ -148 +146 @@
   Vector ProjectedVector;
   Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
-  double radius = 0.;
   double angle = 0.;
 
@@ -172 +169 @@
 
       // correct for negative side
-      radius = ProjectedVector.NormSquared();
+      const double radius = ProjectedVector.NormSquared();
       if (fabs(radius) > MYEPSILON)
         angle = ProjectedVector.Angle(&AngleReferenceVector);
@@ -188 +185 @@
         *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
         *out << Verbose(2) << "New vector: " << *Walker << endl;
-        double tmp = ProjectedVector.NormSquared();
-        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
-          BoundaryTestPair.first->second.first = tmp;
+        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
+        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
+          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
           BoundaryTestPair.first->second.second = Walker;
-          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
-        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
+          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
+        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
           helper.CopyVector(&Walker->x);
           helper.SubtractVector(MolCenter);
-          tmp = helper.NormSquared();
+          const double oldhelperNorm = helper.NormSquared();
           helper.CopyVector(&BoundaryTestPair.first->second.second->x);
           helper.SubtractVector(MolCenter);
-          if (helper.NormSquared() < tmp) {
+          if (helper.NormSquared() < oldhelperNorm) {
             BoundaryTestPair.first->second.second = Walker;
             *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
           } else {
-            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
+            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl;
           }
         } else {
-          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
+          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
         }
       }
@@ -305 +302 @@
 /** Tesselates the convex boundary by finding all boundary points.
  * \param *out output stream for debugging
- * \param *mol molecule structure with Atom's and Bond's
+ * \param *mol molecule structure with Atom's and Bond's.
  * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
  * \param *LCList atoms in LinkedCell list
@@ -311 +308 @@
  * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
  */
-void FindConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename)
+void FindConvexBorder(ofstream *out, const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
 {
   bool BoundaryFreeFlag = false;
@@ -318 +315 @@
   cout << Verbose(1) << "Begin of FindConvexBorder" << endl;
 
-  if (mol->TesselStruct != NULL) // free if allocated
-    delete(mol->TesselStruct);
-  mol->TesselStruct = new class Tesselation;
+  if (TesselStruct != NULL) // free if allocated
+    delete(TesselStruct);
+  TesselStruct = new class Tesselation;
 
   // 1. Find all points on the boundary
   if (BoundaryPoints == NULL) {
       BoundaryFreeFlag = true;
-      BoundaryPoints = GetBoundaryPoints(out, mol);
+      BoundaryPoints = GetBoundaryPoints(out, mol, TesselStruct);
   } else {
       *out << Verbose(1) << "Using given boundary points set." << endl;
@@ -348 +345 @@
   for (int axis = 0; axis < NDIM; axis++)
     for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
-        if (!mol->TesselStruct->AddBoundaryPoint(runner->second.second, 0))
+        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
           *out << Verbose(3) << "WARNING: Point " << *(runner->second.second) << " is already present!" << endl;
 
-  *out << Verbose(2) << "I found " << mol->TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
+  *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
 
@@ -362 +359 @@
 
   // 3a. guess starting triangle
-  mol->TesselStruct->GuessStartingTriangle(out);
+  TesselStruct->GuessStartingTriangle(out);
 
   // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
-  mol->TesselStruct->TesselateOnBoundary(out, mol);
+  TesselStruct->TesselateOnBoundary(out, mol);
 
   // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
-  if (!mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList))
+  if (!TesselStruct->InsertStraddlingPoints(out, mol, LCList))
     *out << Verbose(1) << "Insertion of straddling points failed!" << endl;
 
-  *out << Verbose(2) << "I created " << mol->TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << mol->TesselStruct->LinesOnBoundary.size() << " lines and " << mol->TesselStruct->PointsOnBoundary.size() << " points." << endl;
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
 
   // 4. Store triangles in tecplot file
@@ -380 +377 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
+      WriteTecplotFile(out, tecplot, TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -389 +386 @@
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
+      WriteRaster3dFile(out, rasterplot, TesselStruct, mol);
       rasterplot->close();
       delete(rasterplot);
@@ -400 +397 @@
   do {
     AllConvex = true;
-    for (LineMap::iterator LineRunner = mol->TesselStruct->LinesOnBoundary.begin(); LineRunner != mol->TesselStruct->LinesOnBoundary.end(); LineRunner++) {
+    for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
       line = LineRunner->second;
       *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
@@ -407 +404 @@
 
         // flip the line
-        if (mol->TesselStruct->PickFarthestofTwoBaselines(out, line) == 0.)
+        if (TesselStruct->PickFarthestofTwoBaselines(out, line) == 0.)
           *out << Verbose(1) << "ERROR: Correction of concave baselines failed!" << endl;
         else {
-          mol->TesselStruct->FlipBaseline(out, line);
+          TesselStruct->FlipBaseline(out, line);
           *out << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
         }
@@ -418 +415 @@
 
   // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
-//  if (!mol->TesselStruct->CorrectConcaveTesselPoints(out))
+//  if (!TesselStruct->CorrectConcaveTesselPoints(out))
 //    *out << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
 
-  *out << Verbose(2) << "I created " << mol->TesselStruct->TrianglesOnBoundary.size() << " triangles with " << mol->TesselStruct->LinesOnBoundary.size() << " lines and " << mol->TesselStruct->PointsOnBoundary.size() << " points." << endl;
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
 
   // 4. Store triangles in tecplot file
@@ -429 +426 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
+      WriteTecplotFile(out, tecplot, TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -437 +434 @@
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
+      WriteRaster3dFile(out, rasterplot, TesselStruct, mol);
       rasterplot->close();
       delete(rasterplot);
@@ -458 +455 @@
  * \return true - all removed, false - something went wrong
  */
-bool RemoveAllBoundaryPoints(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename)
+bool RemoveAllBoundaryPoints(ofstream *out, class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
 {
   int i=0;
@@ -481 +478 @@
     // store envelope
     sprintf(number, "-%04d", i++);
-    StoreTrianglesinFile(out, mol, filename, number);
+    StoreTrianglesinFile(out, mol, (const Tesselation *&)TesselStruct, filename, number);
   }
 
@@ -511 +508 @@
  * \return volume difference between the non- and the created convex envelope
  */
-double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename)
+double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
 {
   double volume = 0;
@@ -539 +536 @@
     sprintf(dummy, "-first-%d", run);
     //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-    StoreTrianglesinFile(out, mol, filename, dummy);
+    StoreTrianglesinFile(out, mol, (const Tesselation *&)TesselStruct, filename, dummy);
 
     PointRunner = TesselStruct->PointsOnBoundary.begin();
@@ -555 +552 @@
           volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
           sprintf(dummy, "-first-%d", ++run);
-          StoreTrianglesinFile(out, mol, filename, dummy);
+          StoreTrianglesinFile(out, mol, (const Tesselation *&)TesselStruct, filename, dummy);
           Concavity = true;
           break;
@@ -565 +562 @@
     sprintf(dummy, "-second-%d", run);
     //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-    StoreTrianglesinFile(out, mol, filename, dummy);
+    StoreTrianglesinFile(out, mol, (const Tesselation *&)TesselStruct, filename, dummy);
 
     // second step: PickFarthestofTwoBaselines
@@ -579 +576 @@
         volume += tmp;
         if (tmp != 0.) {
-          mol->TesselStruct->FlipBaseline(out, line);
+          TesselStruct->FlipBaseline(out, line);
           Concavity = true;
         }
@@ -611 +608 @@
 
   CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-  StoreTrianglesinFile(out, mol, filename, "");
+  StoreTrianglesinFile(out, mol, (const Tesselation *&)TesselStruct, filename, "");
 
   // end
@@ -668 +665 @@
  * \param *out output stream for debugging
  * \param *mol molecule with atoms and bonds
+ * \param *&TesselStruct Tesselation with boundary triangles
  * \param *filename prefix of filename
  * \param *extraSuffix intermediate suffix
  */
-void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix)
+void StoreTrianglesinFile(ofstream *out, const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
 {
   // 4. Store triangles in tecplot file
@@ -680 +678 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
+      WriteTecplotFile(out, tecplot, TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -689 +687 @@
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
+      WriteRaster3dFile(out, rasterplot, TesselStruct, mol);
       rasterplot->close();
       delete(rasterplot);
@@ -701 +699 @@
  * \param *configuration needed for path to store convex envelope file
  * \param *mol molecule structure representing the cluster
+ * \param *&TesselStruct Tesselation structure with triangles on return
  * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
  * \param celldensity desired average density in final cell
@@ -722 +721 @@
 
   IsAngstroem = configuration->GetIsAngstroem();
-  GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
-  BoundaryPoints = GetBoundaryPoints(out, mol);
+  GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, TesselStruct, IsAngstroem);
+  BoundaryPoints = GetBoundaryPoints(out, mol, TesselStruct);
   LinkedCell LCList(mol, 10.);
-  FindConvexBorder(out, mol, &LCList, NULL);
+  FindConvexBorder(out, mol, TesselStruct, &LCList, NULL);
 
   // some preparations beforehand
@@ -821 +820 @@
   LinkedCell *LCList[List->ListOfMolecules.size()];
   double phi[NDIM];
+  class Tesselation *TesselStruct[List->ListOfMolecules.size()];
 
   *out << Verbose(0) << "Begin of FillBoxWithMolecule" << endl;
@@ -828 +828 @@
     *out << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
     LCList[i] = new LinkedCell((*ListRunner), 5.); // get linked cell list
-    if ((*ListRunner)->TesselStruct == NULL) {
+    if (TesselStruct[i] == NULL) {
       *out << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
-      FindNonConvexBorder((ofstream *)&cout, (*ListRunner), LCList[i], 5., NULL);
+      FindNonConvexBorder((ofstream *)&cout, (*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
     }
     i++;
@@ -868 +868 @@
         for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
           // get linked cell list
-          if ((*ListRunner)->TesselStruct == NULL) {
+          if (TesselStruct[i] == NULL) {
             *out << Verbose(1) << "ERROR: TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
             FillIt = false;
-          } else
-            FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInnerPoint(out, CurrentPosition, LCList[i++]));
+          } else {
+            FillIt = FillIt && (!TesselStruct[i]->IsInnerPoint(out, CurrentPosition, LCList[i]));
+            i++;
+          }
         }
 
@@ -947 +949 @@
  * \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 atoms in LinkedCell list
+ * \param *&TesselStruct Tesselation filled with points, lines and triangles on boundary on return
+ * \param *&LCList atoms in LinkedCell list
  * \param RADIUS radius of the virtual sphere
  * \param *filename filename prefix for output of vertex data
  */
-void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const double RADIUS, const char *filename = NULL)
+void FindNonConvexBorder(ofstream *out, const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
 {
   bool freeLC = false;
@@ -961 +963 @@
 
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
-  if (mol->TesselStruct == NULL) {
+  if (TesselStruct == NULL) {
     *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
-    mol->TesselStruct = new Tesselation;
+    TesselStruct= new Tesselation;
   } else {
-    delete(mol->TesselStruct);
+    delete(TesselStruct);
     *out << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
-    mol->TesselStruct = new Tesselation;
+    TesselStruct = new Tesselation;
   }
 
@@ -979 +981 @@
 
   // 1. get starting triangle
-  mol->TesselStruct->FindStartingTriangle(out, RADIUS, LCList);
+  TesselStruct->FindStartingTriangle(out, RADIUS, LCList);
 
   // 2. expand from there
-  baseline = mol->TesselStruct->LinesOnBoundary.begin();
+  baseline = TesselStruct->LinesOnBoundary.begin();
   baseline++; // skip first line
-  while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (OneLoopWithoutSuccessFlag)) {
+  while ((baseline != TesselStruct->LinesOnBoundary.end()) || (OneLoopWithoutSuccessFlag)) {
     if (baseline->second->triangles.size() == 1) {
       // 3. find next triangle
-      TesselationFailFlag = mol->TesselStruct->FindNextSuitableTriangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one.
+      TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one.
       OneLoopWithoutSuccessFlag = OneLoopWithoutSuccessFlag || TesselationFailFlag;
       if (!TesselationFailFlag)
@@ -994 +996 @@
       // write temporary envelope
       if (filename != NULL) {
-        if ((DoSingleStepOutput && ((mol->TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
-          mol->TesselStruct->Output(out, filename, mol);
+        if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
+          TesselStruct->Output(out, filename, mol);
         }
       }
-      baseline = mol->TesselStruct->LinesOnBoundary.end();
+      baseline = TesselStruct->LinesOnBoundary.end();
       *out << Verbose(2) << "Baseline set to end." << endl;
     } else {
@@ -1006 +1008 @@
     }
 
-    if ((baseline == mol->TesselStruct->LinesOnBoundary.end()) && (OneLoopWithoutSuccessFlag)) {
-      baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+    if ((baseline == TesselStruct->LinesOnBoundary.end()) && (OneLoopWithoutSuccessFlag)) {
+      baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
       OneLoopWithoutSuccessFlag = false;
     }
@@ -1013 +1015 @@
   }
   // check envelope for consistency
-  CheckListOfBaselines(out, mol->TesselStruct);
+  CheckListOfBaselines(out, TesselStruct);
 
   // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
-  //mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
+  //->InsertStraddlingPoints(out, mol, LCList);
 //  mol->GoToFirst();
 //  class TesselPoint *Runner = NULL;
@@ -1022 +1024 @@
 //    Runner = mol->GetPoint();
 //    *out << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
-//    if (!mol->TesselStruct->IsInnerPoint(out, Runner, LCList)) {
+//    if (!->IsInnerPoint(out, Runner, LCList)) {
 //      *out << Verbose(2) << Runner->Name << " is outside of envelope, adding via degenerated triangles." << endl;
-//      mol->TesselStruct->AddBoundaryPointByDegeneratedTriangle(out, Runner, LCList);
+//      ->AddBoundaryPointByDegeneratedTriangle(out, Runner, LCList);
 //    } else {
 //      *out << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
@@ -1032 +1034 @@
 
   // Purges surplus triangles.
-  mol->TesselStruct->RemoveDegeneratedTriangles();
+  TesselStruct->RemoveDegeneratedTriangles();
 
   // check envelope for consistency
-  CheckListOfBaselines(out, mol->TesselStruct);
+  CheckListOfBaselines(out, TesselStruct);
 
   // write final envelope
-  CalculateConcavityPerBoundaryPoint(out, mol->TesselStruct);
-  StoreTrianglesinFile(out, mol, filename, "");
+  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, (const Tesselation *&)TesselStruct, filename, "");
 
   if (freeLC)

src/boundary.hpp

@@ -48 +48 @@
 /********************************************** declarations *******************************/
 
+double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
+molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
+void FindConvexBorder(ofstream *out, const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
+Vector* FindEmbeddingHole(ofstream *out, MoleculeListClass *mols, molecule *srcmol);
+void FindNextSuitablePoint(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
+void FindNonConvexBorder(ofstream *out, const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LC, const double RADIUS, const char *tempbasename);
+Boundaries *GetBoundaryPoints(ofstream *out, const molecule *mol, Tesselation *&TesselStruct);
+double * GetDiametersOfCluster(ofstream *out, const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem);
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
+bool RemoveAllBoundaryPoints(ofstream *out, class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
+void StoreTrianglesinFile(ofstream *out, const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix);
 double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration);
-double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
-molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
-void FindConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename);
-void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LC, const double RADIUS, const char *tempbasename);
-double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename);
-void FindNextSuitablePoint(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
-Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);
-void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix);
-bool RemoveAllBoundaryPoints(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename);
-Vector* FindEmbeddingHole(ofstream *out, MoleculeListClass *mols, molecule *srcmol);
 
 

src/builder.cpp

@@ -254 +254 @@
         } while ((j != -1) && (i<128));
         if (i >= 2) {
-          first->x.LSQdistance(atoms, i);
+          first->x.LSQdistance((const Vector **)atoms, i);
 
           first->x.Output((ofstream *)&cout);
@@ -592 +592 @@
       {
         cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-        LinkedCell LCList(mol, 10.);
         class Tesselation *TesselStruct = NULL;
-        FindConvexBorder((ofstream *)&cout, mol, &LCList, NULL);
+        const LinkedCell *LCList = NULL;
+        LCList = new LinkedCell(mol, 10.);
+        FindConvexBorder((ofstream *)&cout, mol, TesselStruct, LCList, NULL);
         double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, TesselStruct, configuration);
-        cout << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;
+        cout << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;\
+        delete(LCList);
         delete(TesselStruct);
       }
@@ -720 +722 @@
        cin >> factor[2];
        valid = true;
-       mol->Scale(&factor);
+       mol->Scale((const double ** const)&factor);
        delete[](factor);
       }
@@ -1588 +1590 @@
                 // get the boundary
                 class molecule *Boundary = new molecule(periode);
+                class Tesselation *TesselStruct = NULL;
                 struct ConfigFileBuffer *FileBuffer = NULL;
                 PrepareFileBuffer(argv[argptr], FileBuffer);
                 LoadMolecule(Boundary, FileBuffer, periode, false);
-                LinkedCell LCList(Boundary, 2.*radius);
+                const LinkedCell *LCList = NULL;
+                LCList = new LinkedCell(Boundary, 2.*radius);
                 element *oxygen = periode->FindElement(8);
-                FindNonConvexBorder((ofstream *)&cout, Boundary, &LCList, radius, NULL);
-
-                CorrelationToSurfaceMap *surfacemap = CorrelationToSurface( (ofstream *)&cout, mol, oxygen, Boundary->TesselStruct, &LCList );
+                FindNonConvexBorder((ofstream *)&cout, Boundary, TesselStruct, LCList, radius, NULL);
+                CorrelationToSurfaceMap *surfacemap = CorrelationToSurface( (ofstream *)&cout, mol, oxygen, TesselStruct, LCList );
                 BinPairMap *binmap = BinData( (ofstream *)&cout, surfacemap, 0.5, 0., 0. );
                 OutputCorrelation ( &binoutput, binmap );
                 output.close();
                 binoutput.close();
+                delete(LCList);
                 delete(FileBuffer);
+                delete(TesselStruct);
                 delete(Boundary);
                 argptr+=3;
@@ -1675 +1680 @@
                 cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
               } else {
-                class Tesselation T;
+                class Tesselation *T = NULL;
+                const LinkedCell *LCList = NULL;
                 string filename(argv[argptr+1]);
                 filename.append(".csv");
@@ -1681 +1687 @@
                 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
                 start = clock();
-                LinkedCell LCList(mol, atof(argv[argptr])*2.);
-                FindNonConvexBorder((ofstream *)&cout, mol, &LCList, atof(argv[argptr]), argv[argptr+1]);
-                //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
+                LCList = new LinkedCell(mol, atof(argv[argptr])*2.);
+                FindNonConvexBorder((ofstream *)&cout, mol, T, LCList, atof(argv[argptr]), argv[argptr+1]);
+                //FindDistributionOfEllipsoids((ofstream *)&cout, T, &LCList, N, number, filename.c_str());
                 end = clock();
                 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+                delete(LCList);
                 argptr+=2;
               }
@@ -1806 +1813 @@
                 factor[1] = atof(argv[argptr+1]);
                 factor[2] = atof(argv[argptr+2]);
-                mol->Scale(&factor);
+                mol->Scale((const double ** const)&factor);
                 for (int i=0;i<NDIM;i++) {
                   j += i+1;
@@ -1935 +1942 @@
                 cerr << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl;
               } else {
+                class Tesselation *TesselStruct = NULL;
+                const LinkedCell *LCList = NULL;
                 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
                 cout << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;
                 cout << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl;
-                LinkedCell LCList(mol, 10.);
-                //FindConvexBorder((ofstream *)&cout, mol, &LCList, argv[argptr]);
-                FindNonConvexBorder((ofstream *)&cout, mol, &LCList, 5., argv[argptr+1]);
-//                RemoveAllBoundaryPoints((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
-                double volumedifference = ConvexizeNonconvexEnvelope((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
-                double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, mol->TesselStruct, &configuration);
+                LCList = new LinkedCell(mol, 10.);
+                //FindConvexBorder((ofstream *)&cout, mol, LCList, argv[argptr]);
+                FindNonConvexBorder((ofstream *)&cout, mol, TesselStruct, LCList, 5., argv[argptr+1]);
+//                RemoveAllBoundaryPoints((ofstream *)&cout, TesselStruct, mol, argv[argptr]);
+                double volumedifference = ConvexizeNonconvexEnvelope((ofstream *)&cout, TesselStruct, mol, argv[argptr]);
+                double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, TesselStruct, &configuration);
                 cout << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                 cout << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
+                delete(TesselStruct);
+                delete(LCList);
                 argptr+=2;
               }

src/ellipsoid.cpp

@@ -234 +234 @@
   int index;
   TesselPoint *Candidate = NULL;
-  LinkedNodes *List = NULL;
   *out << Verbose(2) << "Begin of PickRandomPointSet" << endl;
 
@@ -249 +248 @@
     *out << Verbose(2) << "INFO: Center cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " ... ";
     // get random cell
-    List = LC->GetCurrentCell();
+    const LinkedNodes *List = LC->GetCurrentCell();
     if (List == NULL) {  // set index to it
       continue;
@@ -268 +267 @@
       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();
+          const LinkedNodes *List = LC->GetCurrentCell();
           PointsLeft += List->size();
         }
@@ -293 +292 @@
       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();
+          const LinkedNodes *List = LC->GetCurrentCell();
 //          *out << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
           if (List != NULL) {
@@ -300 +299 @@
 //            else
 //              *out << Verbose(2) << "Cell is empty ... " << endl;
-            for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
               if ((current != PickedAtomNrs.end()) && (*current == index)) {
                 Candidate = (*Runner);

src/leastsquaremin.cpp

@@ -16 +16 @@
  * \return sum of square distances
  */
-double LSQ (const gsl_vector * x, void * params)
+double LSQ (const gsl_vector * const x, void * params)
 {
   double sum = 0.;
   struct LSQ_params *par = (struct LSQ_params *)params;
-  Vector **vectors = par->vectors;
+  const Vector ** vectors = par->vectors;
   int num = par->num;
 

src/leastsquaremin.hpp

@@ -31 +31 @@
  */
 struct LSQ_params {
-  Vector **vectors;
+  const Vector **vectors;
   int num;
 };
 
-double LSQ(const gsl_vector * x, void * params);
+double LSQ(const gsl_vector * const x, void * params);
 
 /** Parameter structure for least square minimsation.

src/linkedcell.cpp

@@ -33 +33 @@
  * \param RADIUS edge length of cells
  */
-LinkedCell::LinkedCell(PointCloud *set, double radius)
+LinkedCell::LinkedCell(const PointCloud * const set, const double radius)
 {
   TesselPoint *Walker = NULL;
@@ -109 +109 @@
  * \param RADIUS edge length of cells
  */
-LinkedCell::LinkedCell(LinkedNodes *set, double radius)
+LinkedCell::LinkedCell(LinkedNodes *set, const double radius)
 {
   class TesselPoint *Walker = NULL;
@@ -192 +192 @@
  * \return if all in intervals - true, else -false
  */
-bool LinkedCell::CheckBounds()
+bool LinkedCell::CheckBounds() const
 {
   bool status = true;
@@ -207 +207 @@
  * \return if all in intervals - true, else -false
  */
-bool LinkedCell::CheckBounds(int relative[NDIM])
+bool LinkedCell::CheckBounds(const int relative[NDIM]) const
 {
   bool status = true;
@@ -219 +219 @@
  * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
  */
-LinkedNodes* LinkedCell::GetCurrentCell()
+const LinkedNodes* LinkedCell::GetCurrentCell() const
 {
   if (CheckBounds()) {
@@ -233 +233 @@
  * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[]+relative[] are out of bounds.
  */
-LinkedNodes* LinkedCell::GetRelativeToCurrentCell(int relative[NDIM])
+const LinkedNodes* LinkedCell::GetRelativeToCurrentCell(const int relative[NDIM]) const
 {
   if (CheckBounds(relative)) {
@@ -247 +247 @@
  * \return if the atom is also found in this cell - true, else - false
  */
-bool LinkedCell::SetIndexToNode(const TesselPoint *Walker)
+bool LinkedCell::SetIndexToNode(const TesselPoint * const Walker) const
 {
   bool status = false;
@@ -268 +268 @@
  * \param *upper upper bounds
  */
-void LinkedCell::GetNeighbourBounds(int lower[NDIM], int upper[NDIM])
+void LinkedCell::GetNeighbourBounds(int lower[NDIM], int upper[NDIM]) const
 {
   for (int i=0;i<NDIM;i++) {
@@ -288 +288 @@
  * \return Vector is inside bounding box - true, else - false
  */
-bool LinkedCell::SetIndexToVector(const Vector *x)
+bool LinkedCell::SetIndexToVector(const Vector * const x) const
 {
   bool status = true;

src/linkedcell.hpp

@@ -46 +46 @@
     double RADIUS;    // cell edge length
     int N[NDIM];      // number of cells per axis
-    int n[NDIM];      // temporary variable for current cell per axis
-    int index;        // temporary index variable , access by index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
+    mutable int n[NDIM];      // temporary variable for current cell per axis
+    mutable int index;        // temporary index variable , access by index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
 
     LinkedCell();
-    LinkedCell(PointCloud *set, double RADIUS);
-    LinkedCell(LinkedNodes *set, double radius);
+    LinkedCell(const  PointCloud * const set, const double RADIUS);
+    LinkedCell(LinkedNodes *set, const double radius);
     ~LinkedCell();
-    LinkedNodes* GetCurrentCell();
-    LinkedNodes* GetRelativeToCurrentCell(int relative[NDIM]);
-    bool SetIndexToNode(const TesselPoint *Walker);
-    bool SetIndexToVector(const Vector *x);
-    bool CheckBounds();
-    bool CheckBounds(int relative[NDIM]);
-    void GetNeighbourBounds(int lower[NDIM], int upper[NDIM]);
+    const LinkedNodes* GetCurrentCell()const ;
+    const LinkedNodes* GetRelativeToCurrentCell(const int relative[NDIM])const ;
+    bool SetIndexToNode(const TesselPoint * const Walker)const ;
+    bool SetIndexToVector(const Vector * const x)const ;
+    bool CheckBounds()const ;
+    bool CheckBounds(const int relative[NDIM])const ;
+    void GetNeighbourBounds(int lower[NDIM], int upper[NDIM])const ;
 
     // not implemented yet

src/molecule.cpp

@@ -56 +56 @@
   IndexNr  = -1;
   ActiveFlag = false;
-  TesselStruct = NULL;
 };
 
@@ -64 +63 @@
 molecule::~molecule()
 {
-  if (TesselStruct != NULL)
-    delete(TesselStruct);
   CleanupMolecule();
   delete(first);

src/molecule.hpp

@@ -99 +99 @@
     char name[MAXSTRINGSIZE];         //!< arbitrary name
     int IndexNr;        //!< index of molecule in a MoleculeListClass
-    class Tesselation *TesselStruct;
 
   molecule(periodentafel *teil);
@@ -105 +104 @@
 
   // re-definition of virtual functions from PointCloud
-  Vector *GetCenter(ofstream *out);
-  TesselPoint *GetPoint();
-  TesselPoint *GetTerminalPoint();
-  void GoToNext();
-  void GoToPrevious();
-  void GoToFirst();
-  void GoToLast();
-  bool IsEmpty();
-  bool IsEnd();
+  Vector *GetCenter(ofstream *out) const ;
+  TesselPoint *GetPoint() const ;
+  TesselPoint *GetTerminalPoint() const ;
+  void GoToNext() const ;
+  void GoToPrevious() const ;
+  void GoToFirst() const ;
+  void GoToLast() const ;
+  bool IsEmpty() const ;
+  bool IsEnd() const ;
 
   // templates for allowing global manipulation of all vectors
@@ -217 +216 @@
   void Mirror(const Vector *x);
   void Align(Vector *n);
-  void Scale(double **factor);
+  void Scale(const double ** const factor);
   void DeterminePeriodicCenter(Vector &center);
   Vector * DetermineCenterOfGravity(ofstream *out);
-  Vector * DetermineCenterOfAll(ofstream *out);
+  Vector * DetermineCenterOfAll(ofstream *out) const;
   void SetNameFromFilename(const char *filename);
   void SetBoxDimension(Vector *dim);
@@ -298 +297 @@
   private:
   int last_atom;      //!< number given to last atom
-  atom *InternalPointer;  //!< internal pointer for PointCloud
+  mutable atom *InternalPointer;  //!< internal pointer for PointCloud
 };
 

src/molecule_geometry.cpp

@@ -122 +122 @@
  * \return pointer to center of all vector
  */
-Vector * molecule::DetermineCenterOfAll(ofstream *out)
+Vector * molecule::DetermineCenterOfAll(ofstream *out) const
 {
   atom *ptr = start->next;  // start at first in list
@@ -198 +198 @@
  * \param *factor pointer to scaling factor
  */
-void molecule::Scale(double **factor)
+void molecule::Scale(const double ** const factor)
 {
   atom *ptr = start;

src/molecule_graph.cpp

@@ -102 +102 @@
   double distance, MinDistance, MaxDistance;
   LinkedCell *LC = NULL;
-  LinkedNodes *List = NULL;
-  LinkedNodes *OtherList = NULL;
 
   BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
@@ -131 +129 @@
       for (LC->n[1] = 0; LC->n[1] < LC->N[1]; LC->n[1]++)
         for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
-          List = LC->GetCurrentCell();
+          const LinkedNodes *List = LC->GetCurrentCell();
           //*out << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
           if (List != NULL) {
-            for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
               Walker = AtomMap[(*Runner)->nr];
               //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
@@ -141 +139 @@
                 for (n[1] = -1; n[1] <= 1; n[1]++)
                   for (n[2] = -1; n[2] <= 1; n[2]++) {
-                    OtherList = LC->GetRelativeToCurrentCell(n);
+                    const LinkedNodes *OtherList = LC->GetRelativeToCurrentCell(n);
                     //*out << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
                     if (OtherList != NULL) {
-                      for (LinkedNodes::iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
+                      for (LinkedNodes::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
                         if ((*OtherRunner)->nr > Walker->nr) {
                           OtherWalker = AtomMap[(*OtherRunner)->nr];

src/molecule_pointcloud.cpp

@@ -18 +18 @@
  * \return pointer to allocated with central coordinates
  */
-Vector *molecule::GetCenter(ofstream *out)
+Vector *molecule::GetCenter(ofstream *out) const
 {
   Vector *center = DetermineCenterOfAll(out);
@@ -27 +27 @@
  * \return pointer to atom or NULL if none present
  */
-TesselPoint *molecule::GetPoint()
+TesselPoint *molecule::GetPoint() const
 {
   if ((InternalPointer != start) && (InternalPointer != end))
@@ -38 +38 @@
  * \return pointer to end marker
  */
-TesselPoint *molecule::GetTerminalPoint()
+TesselPoint *molecule::GetTerminalPoint() const
 {
   return end;
@@ -46 +46 @@
  * Stops at last one.
  */
-void molecule::GoToNext()
+void molecule::GoToNext() const
 {
   if (InternalPointer != end)
@@ -55 +55 @@
  * Stops at first one.
  */
-void molecule::GoToPrevious()
+void molecule::GoToPrevious() const
 {
   if (InternalPointer->previous != start)
@@ -63 +63 @@
 /** Goes to first atom.
  */
-void molecule::GoToFirst()
+void molecule::GoToFirst() const
 {
   InternalPointer = start->next;
@@ -70 +70 @@
 /** Goes to last atom.
  */
-void molecule::GoToLast()
+void molecule::GoToLast() const
 {
   InternalPointer = end->previous;
@@ -78 +78 @@
  * \return true - no atoms, false - not empty
  */
-bool molecule::IsEmpty()
+bool molecule::IsEmpty() const
 {
   return (start->next == end);
@@ -86 +86 @@
  * \return true - current atom is last one, false - is not last one
  */
-bool molecule::IsEnd()
+bool molecule::IsEnd() const
 {
   return (InternalPointer == end);

src/moleculelist.cpp

@@ -306 +306 @@
 bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
 {
+  LinkedCell *LCList = NULL;
+  Tesselation *TesselStruct = NULL;
   if ((srcmol == NULL) || (mol == NULL)) {
     cout << Verbose(1) << "ERROR: Either fixed or variable molecule is given as NULL." << endl;
@@ -312 +314 @@
 
   // calculate envelope for *mol
-  LinkedCell *LCList = new LinkedCell(mol, 8.);
-  FindNonConvexBorder((ofstream *)&cout, mol, LCList, 4., NULL);
-  if (mol->TesselStruct == NULL) {
+  LCList = new LinkedCell(mol, 8.);
+  FindNonConvexBorder((ofstream *)&cout, mol, TesselStruct, (const LinkedCell *&)LCList, 4., NULL);
+  if (TesselStruct == NULL) {
     cout << Verbose(1) << "ERROR: Could not tesselate the fixed molecule." << endl;
     return false;
   }
   delete(LCList);
-  LCList = new LinkedCell(mol->TesselStruct, 8.);  // re-create with boundary points only!
+  LCList = new LinkedCell(TesselStruct, 8.);  // re-create with boundary points only!
 
   // prepare index list for bonds
@@ -333 +335 @@
     Walker = Walker->next;
     cout << Verbose(2) << "INFO: Current Walker is " << *Walker << "." << endl;
-    if (!mol->TesselStruct->IsInnerPoint((ofstream *)&cout, Walker->x, LCList)) {
+    if (!TesselStruct->IsInnerPoint((ofstream *)&cout, Walker->x, LCList)) {
       CopyAtoms[Walker->nr] = new atom(Walker);
       mol->AddAtom(CopyAtoms[Walker->nr]);

src/tesselation.cpp

@@ -31 +31 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint *Walker)
+BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker)
 {
   node = Walker;
@@ -73 +73 @@
  * \param &a boundary point
  */
-ostream & operator <<(ostream &ost, BoundaryPointSet &a)
+ostream & operator <<(ostream &ost, const BoundaryPointSet &a)
 {
   ost << "[" << a.Nr << "|" << a.node->Name << " at " << *a.node->node << "]";
@@ -96 +96 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
 {
   // set number
@@ -275 +275 @@
  * \param &a boundary line
  */
-ostream & operator <<(ostream &ost, BoundaryLineSet &a)
+ostream & operator <<(ostream &ost, const  BoundaryLineSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]";
@@ -532 +532 @@
  * \param &a boundary triangle
  */
-ostream &operator <<(ostream &ost, BoundaryTriangleSet &a)
+ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
@@ -642 +642 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-Vector * Tesselation::GetCenter(ofstream *out)
+Vector * Tesselation::GetCenter(ofstream *out) const
 {
   Vector *Center = new Vector(0.,0.,0.);
@@ -657 +657 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-TesselPoint * Tesselation::GetPoint()
+TesselPoint * Tesselation::GetPoint() const
 {
   return (InternalPointer->second->node);
@@ -665 +665 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-TesselPoint * Tesselation::GetTerminalPoint()
-{
-  PointMap::iterator Runner = PointsOnBoundary.end();
+TesselPoint * Tesselation::GetTerminalPoint() const
+{
+  PointMap::const_iterator Runner = PointsOnBoundary.end();
   Runner--;
   return (Runner->second->node);
@@ -675 +675 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-void Tesselation::GoToNext()
+void Tesselation::GoToNext() const
 {
   if (InternalPointer != PointsOnBoundary.end())
@@ -684 +684 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-void Tesselation::GoToPrevious()
+void Tesselation::GoToPrevious() const
 {
   if (InternalPointer != PointsOnBoundary.begin())
@@ -693 +693 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-void Tesselation::GoToFirst()
+void Tesselation::GoToFirst() const
 {
   InternalPointer = PointsOnBoundary.begin();
@@ -700 +700 @@
 /** PointCloud implementation of GoToLast.
  * Uses PointsOnBoundary and STL stuff.
- */   
-void Tesselation::GoToLast()
+ */
+void Tesselation::GoToLast() const
 {
   InternalPointer = PointsOnBoundary.end();
@@ -710 +710 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-bool Tesselation::IsEmpty()
+bool Tesselation::IsEmpty() const
 {
   return (PointsOnBoundary.empty());
@@ -718 +718 @@
  * Uses PointsOnBoundary and STL stuff.
  */   
-bool Tesselation::IsEnd()
+bool Tesselation::IsEnd() const
 {
   return (InternalPointer == PointsOnBoundary.end());
@@ -899 +899 @@
  * \param *cloud cluster of points
  */
-void Tesselation::TesselateOnBoundary(ofstream *out, PointCloud *cloud)
+void Tesselation::TesselateOnBoundary(ofstream *out, const PointCloud * const cloud)
 {
   bool flag;
@@ -1097 +1097 @@
  * \return true - all straddling points insert, false - something went wrong
  */
-bool Tesselation::InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC)
+bool Tesselation::InsertStraddlingPoints(ofstream *out, const PointCloud *cloud, const LinkedCell *LC)
 {
   Vector Intersection, Normal;
@@ -1207 +1207 @@
  * \return true - new point was added, false - point already present
  */
-bool
-Tesselation::AddBoundaryPoint(TesselPoint *Walker, int n)
+bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n)
 {
   PointTestPair InsertUnique;
@@ -1229 +1228 @@
  * @param n index for this point in Tesselation::TPS array
  */
-void
-Tesselation::AddTesselationPoint(TesselPoint* Candidate, int n)
+void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n)
 {
   PointTestPair InsertUnique;
@@ -1252 +1250 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
+void Tesselation::AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) {
   bool insertNewLine = true;
 
@@ -1290 +1288 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
+void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
 {
   cout << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
@@ -1323 +1321 @@
  * \param nr triangle number
  */
-void Tesselation::AddTesselationTriangle(int nr)
+void Tesselation::AddTesselationTriangle(const int nr)
 {
   cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
@@ -1554 +1552 @@
  * \param *LC LinkedCell structure with neighbouring TesselPoint's
  */
-void Tesselation::FindStartingTriangle(ofstream *out, const double RADIUS, LinkedCell *LC)
+void Tesselation::FindStartingTriangle(ofstream *out, const double RADIUS, const LinkedCell *LC)
 {
   cout << Verbose(1) << "Begin of FindStartingTriangle\n";
   int i = 0;
-  LinkedNodes *List = NULL;
   TesselPoint* FirstPoint = NULL;
   TesselPoint* SecondPoint = NULL;
@@ -1580 +1577 @@
     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();
+        const LinkedNodes *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 (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
+          for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) {
             if ((*Runner)->node->x[i] > maxCoordinate[i]) {
               cout << Verbose(2) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
@@ -1645 +1642 @@
 
     //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    FindThirdPointForTesselation(
-      Oben, SearchDirection, helper, BLS[0], NULL, *&OptCandidates, &ShortestAngle, RADIUS, LC
-    );
+    FindThirdPointForTesselation(Oben, SearchDirection, helper, BLS[0], NULL, *&OptCandidates, &ShortestAngle, RADIUS, LC);
     cout << Verbose(1) << "List of third Points is ";
     for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
@@ -1713 +1708 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
 {
   cout << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
@@ -1774 +1769 @@
 
     // add third point
-    FindThirdPointForTesselation(
-      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, OptCandidates,
-      &ShortestAngle, RADIUS, LC
-    );
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, OptCandidates, &ShortestAngle, RADIUS, LC);
 
   } else {
@@ -1813 +1805 @@
       AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
 
-      if (CheckLineCriteriaForDegeneratedTriangle(TPS)) {
+      if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet **)TPS)) {
         AddTesselationLine(TPS[0], TPS[1], 0);
         AddTesselationLine(TPS[0], TPS[2], 1);
@@ -1842 +1834 @@
         // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
         // i.e. at least one of the three lines must be present with TriangleCount <= 1
-        if (CheckLineCriteriaForDegeneratedTriangle(TPS)) {
+        if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet **)TPS)) {
           AddTesselationLine(TPS[0], TPS[1], 0);
           AddTesselationLine(TPS[0], TPS[2], 1);
@@ -1949 +1941 @@
 };
 
-void Tesselation::PrintAllBoundaryPoints(ofstream *out)
+void Tesselation::PrintAllBoundaryPoints(ofstream *out) const
 {
   // print all lines
   *out << Verbose(1) << "Printing all boundary points for debugging:" << endl;
-  for (PointMap::iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++)
+  for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++)
     *out << Verbose(2) << *(PointRunner->second) << endl;
 };
 
-void Tesselation::PrintAllBoundaryLines(ofstream *out)
+void Tesselation::PrintAllBoundaryLines(ofstream *out) const
 {
   // print all lines
   *out << Verbose(1) << "Printing all boundary lines for debugging:" << endl;
-  for (LineMap::iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
+  for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
     *out << Verbose(2) << *(LineRunner->second) << endl;
 };
 
-void Tesselation::PrintAllBoundaryTriangles(ofstream *out)
+void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const
 {
   // print all triangles
   *out << Verbose(1) << "Printing all boundary triangles for debugging:" << endl;
-  for (TriangleMap::iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
+  for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
     *out << Verbose(2) << *(TriangleRunner->second) << endl;
 };
@@ -2165 +2157 @@
  * \param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
+void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC)
 {
   cout << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl;
   Vector AngleCheck;
   class TesselPoint* Candidate = NULL;
-  double norm = -1., angle;
-  LinkedNodes *List = NULL;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  double norm = -1.;
+  double angle = 0.;
+  int N[NDIM];
+  int Nlower[NDIM];
+  int Nupper[NDIM];
 
   if (LC->SetIndexToNode(a)) {  // get cell for the starting point
@@ -2198 +2192 @@
     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();
+        const LinkedNodes *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 (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             Candidate = (*Runner);
             // check if we only have one unique point yet ...
@@ -2285 +2279 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint  *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
+void Tesselation::FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint  *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC)
 {
   Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
@@ -2294 +2288 @@
   Vector NewNormalVector;   // normal vector of the Candidate's triangle
   Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
-  LinkedNodes *List = NULL;
   double CircleRadius; // radius of this circle
   double radius;
@@ -2357 +2350 @@
         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();
+            const LinkedNodes *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 (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+              for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
                 Candidate = (*Runner);
 
@@ -2471 +2464 @@
  * \return point which is shared or NULL if none
  */
-class BoundaryPointSet *Tesselation::GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
-{
-  class BoundaryLineSet * lines[2] =
-    { line1, line2 };
+class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const
+{
+  const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
   map<int, class BoundaryPointSet *> OrderMap;
@@ -2502 +2494 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
  */
-list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC)
+list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(ofstream *out, const Vector *x, const LinkedCell* LC) const
 {
   TesselPoint *trianglePoints[3];
@@ -2522 +2514 @@
   if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
     *out << Verbose(3) << "Point is right on a tesselation point, no nearest triangle." << endl;
-    PointMap::iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
+    PointMap::const_iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
     triangles = new list<BoundaryTriangleSet*>;
     if (PointRunner != PointsOnBoundary.end()) {
@@ -2578 +2570 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x, LinkedCell* LC)
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, const Vector *x, const LinkedCell* LC) const
 {
   class BoundaryTriangleSet *result = NULL;
@@ -2614 +2606 @@
  * @return true if the point is inside the tesselation structure, false otherwise
  */
-bool Tesselation::IsInnerPoint(ofstream *out, Vector Point, LinkedCell* LC)
+bool Tesselation::IsInnerPoint(ofstream *out, const Vector &Point, const LinkedCell* const LC) const
 {
   class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, &Point, LC);
@@ -2644 +2636 @@
  * @return true if the point is inside the tesselation structure, false otherwise
  */
-bool Tesselation::IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC)
+bool Tesselation::IsInnerPoint(ofstream *out, const TesselPoint * const Point, const LinkedCell* const LC) const
 {
   return IsInnerPoint(out, *(Point->node), LC);
@@ -2655 +2647 @@
  * @return set of the all points linked to the provided one
  */
-set<TesselPoint*> * Tesselation::GetAllConnectedPoints(ofstream *out, TesselPoint* Point)
+set<TesselPoint*> * Tesselation::GetAllConnectedPoints(ofstream *out, const TesselPoint* const Point) const
 {
   set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
@@ -2665 +2657 @@
 
   // find the respective boundary point
-  PointMap::iterator PointRunner = PointsOnBoundary.find(Point->nr);
+  PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
   if (PointRunner != PointsOnBoundary.end()) {
     ReferencePoint = PointRunner->second;
@@ -2675 +2667 @@
   // little trick so that we look just through lines connect to the BoundaryPoint
   // OR fall-back to look through all lines if there is no such BoundaryPoint
-  LineMap *Lines = &LinesOnBoundary;
+  const LineMap *Lines;;
   if (ReferencePoint != NULL)
     Lines = &(ReferencePoint->lines);
-  LineMap::iterator findLines = Lines->begin();
+  else
+    Lines = &LinesOnBoundary;
+  LineMap::const_iterator findLines = Lines->begin();
   while (findLines != Lines->end()) {
    takePoint = false;
@@ -2719 +2713 @@
  * @return list of the all points linked to the provided one
  */
-list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector *Reference)
+list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, const TesselPoint* const Point, const Vector * const Reference) const
 {
   map<double, TesselPoint*> anglesOfPoints;
@@ -2733 +2727 @@
 
   // calculate central point
-  for (set<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
+  for (set<TesselPoint*>::const_iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //*out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -2796 +2790 @@
  * @return list of the all points linked to the provided one
  */
-list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(ofstream *out, const TesselPoint* const Point) const
 {
   map<double, TesselPoint*> anglesOfPoints;
@@ -2813 +2807 @@
 
   // find the respective boundary point
-  PointMap::iterator PointRunner = PointsOnBoundary.find(Point->nr);
+  PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
   if (PointRunner != PointsOnBoundary.end()) {
     ReferencePoint = PointRunner->second;
@@ -2911 +2905 @@
  * @return list of the closed paths
  */
-list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(ofstream *out, const TesselPoint* const Point) const
 {
   list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(out, Point);
@@ -2972 +2966 @@
  * \return pointer to allocated list of triangles
  */
-set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(ofstream *out, class BoundaryPointSet *Point)
+set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(ofstream *out, const BoundaryPointSet * const Point) const
 {
   set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
@@ -2980 +2974 @@
   } else {
     // go through its lines and insert all triangles
-    for (LineMap::iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
+    for (LineMap::const_iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
       for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
       connectedTriangles->insert(TriangleRunner->second);
@@ -3224 +3218 @@
  *         will usually be one, in case of degeneration, there will be two
  */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(TesselPoint* Points[3]) 
+list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
 {
   list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
-  LineMap::iterator FindLine;
-  PointMap::iterator FindPoint;
-  TriangleMap::iterator FindTriangle;
+  LineMap::const_iterator FindLine;
+  TriangleMap::const_iterator FindTriangle;
   class BoundaryPointSet *TrianglePoints[3];
 
   for (int i = 0; i < 3; i++) {
-    FindPoint = PointsOnBoundary.find(Points[i]->nr);
+    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
     if (FindPoint != PointsOnBoundary.end()) {
       TrianglePoints[i] = FindPoint->second;
@@ -3544 +3537 @@
  * \param *cloud PointCloud structure with all nodes
  */
-void Tesselation::Output(ofstream *out, const char *filename, PointCloud *cloud)
+void Tesselation::Output(ofstream *out, const char *filename, const PointCloud * const cloud)
 {
   ofstream *tempstream = NULL;

src/tesselation.hpp

@@ -83 +83 @@
   public:
     BoundaryPointSet();
-    BoundaryPointSet(TesselPoint *Walker);
+    BoundaryPointSet(TesselPoint * Walker);
     ~BoundaryPointSet();
 
@@ -95 +95 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryPointSet &a);
+ostream & operator << (ostream &ost, const BoundaryPointSet &a);
 
 // ======================================================== class BoundaryLineSet ==========================================
@@ -102 +102 @@
   public:
     BoundaryLineSet();
-    BoundaryLineSet(class BoundaryPointSet *Point[2], int number);
+    BoundaryLineSet(class BoundaryPointSet *Point[2], const int number);
     ~BoundaryLineSet();
 
@@ -116 +116 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryLineSet &a);
+ostream & operator << (ostream &ost, const BoundaryLineSet &a);
 
 // ======================================================== class BoundaryTriangleSet =======================================
@@ -142 +142 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryTriangleSet &a);
+ostream & operator << (ostream &ost, const BoundaryTriangleSet &a);
 
 // =========================================================== class TESSELPOINT ===========================================
@@ -170 +170 @@
   virtual ~PointCloud();
 
-  virtual Vector *GetCenter(ofstream *out) { return NULL; };
-  virtual TesselPoint *GetPoint() { return NULL; };
-  virtual TesselPoint *GetTerminalPoint() { return NULL; };
-  virtual void GoToNext() {};
-  virtual void GoToPrevious() {};
-  virtual void GoToFirst() {};
-  virtual void GoToLast() {};
-  virtual bool IsEmpty() { return false; };
-  virtual bool IsEnd() { return false; };
+  virtual Vector *GetCenter(ofstream *out) const { return NULL; };
+  virtual TesselPoint *GetPoint() const { return NULL; };
+  virtual TesselPoint *GetTerminalPoint() const { return NULL; };
+  virtual void GoToNext() const {};
+  virtual void GoToPrevious() const {};
+  virtual void GoToFirst() const {};
+  virtual void GoToLast() const {};
+  virtual bool IsEmpty() const { return false; };
+  virtual bool IsEnd() const { return false; };
 };
 
@@ -204 +204 @@
     virtual ~Tesselation();
 
-    void AddTesselationPoint(TesselPoint* Candidate, int n);
-    void AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
+    void AddTesselationPoint(TesselPoint* Candidate, const int n);
+    void AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
+    void AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddTesselationTriangle();
-    void AddTesselationTriangle(int nr);
+    void AddTesselationTriangle(const int nr);
     void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle);
     void RemoveTesselationLine(class BoundaryLineSet *line);
     void RemoveTesselationPoint(class BoundaryPointSet *point);
 
-    class BoundaryPointSet *GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2);
 
     // concave envelope
-    void FindStartingTriangle(ofstream *out, const double RADIUS, class LinkedCell *LC);
-    void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, class LinkedCell *LC);
-    void FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, class LinkedCell *LC);
-    bool FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC);
+    void FindStartingTriangle(ofstream *out, const double RADIUS, const LinkedCell *LC);
+    void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC);
+    void FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC);
+    bool FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
     int CheckPresenceOfTriangle(ofstream *out, class TesselPoint *Candidates[3]);
     class BoundaryTriangleSet * GetPresentTriangle(ofstream *out, TesselPoint *Candidates[3]);
 
     // convex envelope
-    void TesselateOnBoundary(ofstream *out, PointCloud *cloud);
+    void TesselateOnBoundary(ofstream *out, const PointCloud * const cloud);
     void GuessStartingTriangle(ofstream *out);
-    bool InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC);
+    bool InsertStraddlingPoints(ofstream *out, const PointCloud *cloud, const LinkedCell *LC);
     double RemovePointFromTesselatedSurface(ofstream *out, class BoundaryPointSet *point);
     class BoundaryLineSet * FlipBaseline(ofstream *out, class BoundaryLineSet *Base);
@@ -236 +235 @@
     void AddBoundaryPointByDegeneratedTriangle(ofstream *out, class TesselPoint *point, LinkedCell *LC);
 
-    set<TesselPoint*> * GetAllConnectedPoints(ofstream *out, TesselPoint* Point);
-    set<BoundaryTriangleSet*> *GetAllTriangles(ofstream *out, class BoundaryPointSet *Point);
-    list<list<TesselPoint*> *> * GetPathsOfConnectedPoints(ofstream *out, TesselPoint* Point);
-    list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(ofstream *out, TesselPoint* Point);
-    list<TesselPoint*> * GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector *Reference = NULL);
-    list<BoundaryTriangleSet*> *FindTriangles(TesselPoint* Points[3]);
-    list<BoundaryTriangleSet*> * FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC);
-    class BoundaryTriangleSet * FindClosestTriangleToPoint(ofstream *out, Vector *x, LinkedCell* LC);
-    bool IsInnerPoint(ofstream *out, Vector Point, LinkedCell* LC);
-    bool IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC);
-    bool AddBoundaryPoint(TesselPoint *Walker, int n);
+    set<TesselPoint*> * GetAllConnectedPoints(ofstream *out, const TesselPoint* const Point) const;
+    set<BoundaryTriangleSet*> *GetAllTriangles(ofstream *out, const BoundaryPointSet * const Point) const;
+    list<list<TesselPoint*> *> * GetPathsOfConnectedPoints(ofstream *out, const TesselPoint* const Point) const;
+    list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(ofstream *out, const TesselPoint* const Point) const;
+    list<TesselPoint*> * GetCircleOfConnectedPoints(ofstream *out, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
+    class BoundaryPointSet *GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
+    list<BoundaryTriangleSet*> *FindTriangles(const TesselPoint* const Points[3]) const;
+    list<BoundaryTriangleSet*> * FindClosestTrianglesToPoint(ofstream *out, const Vector *x, const LinkedCell* LC) const;
+    class BoundaryTriangleSet * FindClosestTriangleToPoint(ofstream *out, const Vector *x, const LinkedCell* LC) const;
+    bool IsInnerPoint(ofstream *out, const Vector &Point, const LinkedCell* const LC) const;
+    bool IsInnerPoint(ofstream *out, const TesselPoint * const Point, const LinkedCell* const LC) const;
+    bool AddBoundaryPoint(TesselPoint * Walker, const int n);
 
     // print for debugging
-    void PrintAllBoundaryPoints(ofstream *out);
-    void PrintAllBoundaryLines(ofstream *out);
-    void PrintAllBoundaryTriangles(ofstream *out);
+    void PrintAllBoundaryPoints(ofstream *out) const;
+    void PrintAllBoundaryLines(ofstream *out) const;
+    void PrintAllBoundaryTriangles(ofstream *out) const;
 
     // store envelope in file
-    void Output(ofstream *out, const char *filename, PointCloud *cloud);
+    void Output(ofstream *out, const char *filename, const PointCloud * const cloud);
 
     PointMap PointsOnBoundary;
@@ -264 +264 @@
 
     // PointCloud implementation for PointsOnBoundary
-    virtual Vector *GetCenter(ofstream *out);
-    virtual TesselPoint *GetPoint();
-    virtual TesselPoint *GetTerminalPoint();
-    virtual void GoToNext();
-    virtual void GoToPrevious();
-    virtual void GoToFirst();
-    virtual void GoToLast();
-    virtual bool IsEmpty();
-    virtual bool IsEnd();
+    virtual Vector *GetCenter(ofstream *out) const;
+    virtual TesselPoint *GetPoint() const;
+    virtual TesselPoint *GetTerminalPoint() const;
+    virtual void GoToNext() const;
+    virtual void GoToPrevious() const;
+    virtual void GoToFirst() const;
+    virtual void GoToLast() const;
+    virtual bool IsEmpty() const;
+    virtual bool IsEnd() const;
 
     class BoundaryPointSet *BPS[2];
@@ -283 +283 @@
     class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
 
-    PointMap::iterator InternalPointer;
+    mutable PointMap::const_iterator InternalPointer;
 };
 

src/tesselationhelpers.cpp

@@ -452 +452 @@
  * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
  */
-bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3])
+bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet *nodes[3])
 {
   bool result = false;
@@ -461 +461 @@
     for (int j=i+1; j<3; j++) {
       if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
-        LineMap::iterator FindLine;
-        pair<LineMap::iterator,LineMap::iterator> FindPair;
+        LineMap::const_iterator FindLine;
+        pair<LineMap::const_iterator,LineMap::const_iterator> FindPair;
         FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
         for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
@@ -486 +486 @@
 /** Sort function for the candidate list.
  */
-bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
+bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
 {
   Vector BaseLineVector, OrthogonalVector, helper;
@@ -538 +538 @@
 TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
 {
-  LinkedNodes *List = NULL;
   TesselPoint* closestPoint = NULL;
   TesselPoint* secondClosestPoint = NULL;
@@ -556 +555 @@
     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();
+        const LinkedNodes *List = LC->GetCurrentCell();
         //cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
         if (List != NULL) {
-          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             helper.CopyVector(Point);
             helper.SubtractVector((*Runner)->node);
@@ -591 +590 @@
  * @return point which is closest to the provided one, NULL if none found
  */
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
-{
-  LinkedNodes *List = NULL;
+TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
+{
   TesselPoint* closestPoint = NULL;
   SecondPoint = NULL;
@@ -611 +609 @@
     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();
+        const LinkedNodes *List = LC->GetCurrentCell();
         //cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
         if (List != NULL) {
-          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             helper.CopyVector(Point);
             helper.SubtractVector((*Runner)->node);
@@ -651 +649 @@
  * \return Vector on reference line that has closest distance
  */
-Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
+Vector * GetClosestPointBetweenLine(ofstream *out, const BoundaryLineSet *Base, const BoundaryLineSet *OtherBase)
 {
   // construct the plane of the two baselines (i.e. take both their directional vectors)
@@ -691 +689 @@
  * \return distance between \a *x and plane defined by \a *triangle, -1 - if something went wrong
  */
-double DistanceToTrianglePlane(ofstream *out, Vector *x, BoundaryTriangleSet *triangle)
+double DistanceToTrianglePlane(ofstream *out, const Vector *x, const BoundaryTriangleSet *triangle)
 {
   double distance = 0.;
@@ -707 +705 @@
  * \param *mol molecule structure with atom positions
  */
-void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, PointCloud *cloud)
+void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, const Tesselation *Tess, const PointCloud *cloud)
 {
   TesselPoint *Walker = NULL;
@@ -728 +726 @@
 
     *vrmlfile << "# All tesselation triangles" << endl;
-    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+    for (TriangleMap::const_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
@@ -750 +748 @@
  * \param *mol molecule structure with atom positions
  */
-void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
+void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, const Tesselation *Tess, const PointCloud *cloud)
 {
   Vector helper;
@@ -775 +773 @@
  * \param *mol molecule structure with atom positions
  */
-void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
+void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, const Tesselation *Tess, const PointCloud *cloud)
 {
   TesselPoint *Walker = NULL;
@@ -797 +795 @@
     *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++) {
+    for (TriangleMap::const_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
@@ -820 +818 @@
  * \param N arbitrary number to differentiate various zones in the tecplot format
  */
-void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, PointCloud *cloud, int N)
+void WriteTecplotFile(ofstream *out, ofstream *tecplot, const Tesselation *TesselStruct, const PointCloud *cloud, const int N)
 {
   if ((tecplot != NULL) && (TesselStruct != NULL)) {
@@ -840 +838 @@
     int Counter = 1;
     TesselPoint *Walker = NULL;
-    for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
+    for (PointMap::const_iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
       Walker = target->second->node;
       LookupList[Walker->nr] = Counter++;
@@ -847 +845 @@
     *tecplot << endl;
     // print connectivity
-    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
+    for (TriangleMap::const_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;
@@ -861 +859 @@
  * \param *TesselStruct pointer to Tesselation structure
  */
-void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct)
+void CalculateConcavityPerBoundaryPoint(ofstream *out, const Tesselation *TesselStruct)
 {
   class BoundaryPointSet *point = NULL;
@@ -868 +866 @@
   //*out << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
   // calculate remaining concavity
-  for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     point = PointRunner->second;
     *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
@@ -888 +886 @@
  * \return true - all have exactly two triangles, false - some not, list is printed to screen
  */
-bool CheckListOfBaselines(ofstream *out, Tesselation *TesselStruct)
-{
-  LineMap::iterator testline;
+bool CheckListOfBaselines(ofstream *out, const Tesselation *TesselStruct)
+{
+  LineMap::const_iterator testline;
   bool result = false;
   int counter = 0;

src/tesselationhelpers.hpp

@@ -57 +57 @@
 double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector);
 
-bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3]);
-bool SortCandidates(class CandidateForTesselation* candidate1, class CandidateForTesselation* candidate2);
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC);
-TesselPoint* FindSecondClosestPoint(const Vector*, LinkedCell*);
-Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase);
+bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet *nodes[3]);
+bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2);
+TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC);
+TesselPoint* FindSecondClosestPoint(const Vector*, const LinkedCell*);
+Vector * GetClosestPointBetweenLine(ofstream *out, const BoundaryLineSet *Base, const BoundaryLineSet *OtherBase);
 
-void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, PointCloud *cloud, int N);
-void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud);
-void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud);
-void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, PointCloud *cloud);
-void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct);
-double DistanceToTrianglePlane(ofstream *out, Vector *x, BoundaryTriangleSet *triangle);
+void WriteTecplotFile(ofstream *out, ofstream *tecplot, const Tesselation *TesselStruct, const PointCloud *cloud, const int N);
+void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, const Tesselation *Tess, const PointCloud *cloud);
+void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, const Tesselation *Tess, const PointCloud *cloud);
+void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, const Tesselation *Tess, const PointCloud *cloud);
+void CalculateConcavityPerBoundaryPoint(ofstream *out, const Tesselation *TesselStruct);
+double DistanceToTrianglePlane(ofstream *out, const Vector *x, const BoundaryTriangleSet *triangle);
 
-bool CheckListOfBaselines(ofstream *out, Tesselation *TesselStruct);
+bool CheckListOfBaselines(ofstream *out, const Tesselation *TesselStruct);
 
 

src/unittests/ActOnAllUnitTest.cpp

@@ -62 +62 @@
 
   // scaling by value
-  VL.ActOnAll( (void (Vector::*)(double)) &Vector::Scale, 2. );
+  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, 2. );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
 
-  VL.ActOnAll( (void (Vector::*)(double)) &Vector::Scale, 0.5 );
+  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, 0.5 );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
 
   // scaling by ref
-  VL.ActOnAll( (void (Vector::*)(double *)) &Vector::Scale, &factor );
+  VL.ActOnAll( (void (Vector::*)(const double * const)) &Vector::Scale, (const double * const)&factor );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
 
-  VL.ActOnAll( (void (Vector::*)(double *)) &Vector::Scale, &inverse );
+  VL.ActOnAll( (void (Vector::*)(const double * const)) &Vector::Scale, (const double * const)&inverse );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
 
@@ -82 +82 @@
     inverses[i] = 1./factors[i];
   }
-  VL.ActOnAll( (void (Vector::*)(double **)) &Vector::Scale, &factors );
+  VL.ActOnAll( (void (Vector::*)(const double ** const)) &Vector::Scale, (const double ** const)&factors );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
 
-  VL.ActOnAll( (void (Vector::*)(double **)) &Vector::Scale, &inverses );
+  VL.ActOnAll( (void (Vector::*)(const double ** const)) &Vector::Scale, (const double ** const)&inverses );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
 };

src/unittests/AnalysisCorrelationToPointUnitTest.cpp

@@ -77 +77 @@
 
   // init maps
-  pointmap = CorrelationToPoint( (ofstream *)&cout, TestMolecule, hydrogen, point );
+  pointmap = CorrelationToPoint( (ofstream *)&cout, (const molecule * const)TestMolecule, (const element * const)hydrogen, (const Vector *)point );
   binmap = NULL;
 

src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

@@ -81 +81 @@
   // init tesselation and linked cell
   Surface = new Tesselation;
-  TestMolecule->TesselStruct = Surface;
-  FindNonConvexBorder((ofstream *)&cerr, TestMolecule, LC, 2.5, NULL);
+  FindNonConvexBorder((ofstream *)&cerr, TestMolecule, Surface, (const LinkedCell *&)LC, 2.5, NULL);
   LC = new LinkedCell(TestMolecule, 5.);
   CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() );
@@ -123 +122 @@
   // remove
   delete(TestMolecule);
-  // note that Surface and all the atoms are cleaned by TestMolecule
+  delete(Surface);
+  // note that all the atoms are cleaned by TestMolecule
   delete(LC);
   delete(tafel);

src/vector.cpp

@@ -21 +21 @@
 /** Constructor of class vector.
  */
-Vector::Vector(double x1, double x2, double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
+Vector::Vector(const double x1, const double x2, const double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
 
 /** Desctructor of class vector.
@@ -31 +31 @@
  * \return \f$| x - y |^2\f$
  */
-double Vector::DistanceSquared(const Vector *y) const
+double Vector::DistanceSquared(const Vector * const y) const
 {
   double res = 0.;
@@ -43 +43 @@
  * \return \f$| x - y |\f$
  */
-double Vector::Distance(const Vector *y) const
+double Vector::Distance(const Vector * const y) const
 {
   double res = 0.;
@@ -56 +56 @@
  * \return \f$| x - y |\f$
  */
-double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const
+double Vector::PeriodicDistance(const Vector * const y, const double * const cell_size) const
 {
   double res = Distance(y), tmp, matrix[NDIM*NDIM];
@@ -94 +94 @@
  * \return \f$| x - y |^2\f$
  */
-double Vector::PeriodicDistanceSquared(const Vector *y, const double *cell_size) const
+double Vector::PeriodicDistanceSquared(const Vector * const y, const double * const cell_size) const
 {
   double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
@@ -131 +131 @@
  * Tries to translate a vector into each adjacent neighbouring cell.
  */
-void Vector::KeepPeriodic(ofstream *out, double *matrix)
+void Vector::KeepPeriodic(ofstream *out, const double * const matrix)
 {
 //  int N[NDIM];
@@ -162 +162 @@
  * \return \f$\langle x, y \rangle\f$
  */
-double Vector::ScalarProduct(const Vector *y) const
+double Vector::ScalarProduct(const Vector * const y) const
 {
   double res = 0.;
@@ -177 +177 @@
  *  \return \f$ x \times y \f&
  */
-void Vector::VectorProduct(const Vector *y)
+void Vector::VectorProduct(const Vector * const y)
 {
   Vector tmp;
@@ -184 +184 @@
   tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
   this->CopyVector(&tmp);
-
 };
 
@@ -192 +191 @@
  * \return \f$\langle x, y \rangle\f$
  */
-void Vector::ProjectOntoPlane(const Vector *y)
+void Vector::ProjectOntoPlane(const Vector * const y)
 {
   Vector tmp;
@@ -217 +216 @@
  * \return true -  \a this contains intersection point on return, false - line is parallel to plane
  */
-bool Vector::GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *Origin, Vector *LineVector)
+bool Vector::GetIntersectionWithPlane(ofstream *out, const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
 {
   double factor;
@@ -264 +263 @@
  * \return distance to plane
  */
-double Vector::DistanceToPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset)
+double Vector::DistanceToPlane(ofstream *out, const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
 {
   Vector temp;
@@ -292 +291 @@
  * \return true - \a this will contain the intersection on return, false - lines are parallel
  */
-bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b, const Vector *PlaneNormal)
+bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
 {
   bool result = true;
@@ -375 +374 @@
  * \param *y array to second vector
  */
-void Vector::ProjectIt(const Vector *y)
+void Vector::ProjectIt(const Vector * const y)
 {
   Vector helper(*y);
@@ -386 +385 @@
  * \return Vector
  */
-Vector Vector::Projection(const Vector *y) const
+Vector Vector::Projection(const Vector * const y) const
 {
   Vector helper(*y);
@@ -435 +434 @@
 /** Zeros all components of this vector.
  */
-void Vector::One(double one)
+void Vector::One(const double one)
 {
   for (int i=NDIM;i--;)
@@ -443 +442 @@
 /** Initialises all components of this vector.
  */
-void Vector::Init(double x1, double x2, double x3)
+void Vector::Init(const double x1, const double x2, const double x3)
 {
   x[0] = x1;
@@ -469 +468 @@
  * @return true - vector is normalized, false - vector is not
  */
-bool Vector::IsNormalTo(const Vector *normal) const
+bool Vector::IsNormalTo(const Vector * const normal) const
 {
   if (ScalarProduct(normal) < MYEPSILON)
@@ -481 +480 @@
  * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
  */
-double Vector::Angle(const Vector *y) const
+double Vector::Angle(const Vector * const y) const
 {
   double norm1 = Norm(), norm2 = y->Norm();
@@ -500 +499 @@
  * \param alpha rotation angle in radian
  */
-void Vector::RotateVector(const Vector *axis, const double alpha)
+void Vector::RotateVector(const Vector * const axis, const double alpha)
 {
   Vector a,y;
@@ -656 +655 @@
  * \param *factor pointer to scaling factor
  */
-void Vector::Scale(double **factor)
+void Vector::Scale(const double ** const factor)
 {
   for (int i=NDIM;i--;)
@@ -662 +661 @@
 };
 
-void Vector::Scale(double *factor)
+void Vector::Scale(const double * const factor)
 {
   for (int i=NDIM;i--;)
@@ -668 +667 @@
 };
 
-void Vector::Scale(double factor)
+void Vector::Scale(const double factor)
 {
   for (int i=NDIM;i--;)
@@ -677 +676 @@
  * \param trans[] translation vector.
  */
-void Vector::Translate(const Vector *trans)
+void Vector::Translate(const Vector * const trans)
 {
   for (int i=NDIM;i--;)
@@ -687 +686 @@
  * \param *Minv inverse matrix
  */
-void Vector::WrapPeriodically(const double *M, const double *Minv)
+void Vector::WrapPeriodically(const double * const M, const double * const Minv)
 {
   MatrixMultiplication(Minv);
@@ -704 +703 @@
  * \param *matrix NDIM_NDIM array
  */
-void Vector::MatrixMultiplication(const double *M)
+void Vector::MatrixMultiplication(const double * const M)
 {
   Vector C;
@@ -719 +718 @@
  * \param *matrix NDIM_NDIM array
  */
-double * Vector::InverseMatrix(double *A)
+double * Vector::InverseMatrix( const double * const A)
 {
   double *B = Malloc<double>(NDIM * NDIM, "Vector::InverseMatrix: *B");
@@ -746 +745 @@
  * \param *matrix NDIM_NDIM array
  */
-void Vector::InverseMatrixMultiplication(const double *A)
+void Vector::InverseMatrixMultiplication(const double * const A)
 {
   Vector C;
@@ -786 +785 @@
  * \param *factors three-component vector with the factor for each given vector
  */
-void Vector::LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors)
+void Vector::LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors)
 {
   for(int i=NDIM;i--;)
@@ -795 +794 @@
  * \param n[] normal vector of mirror plane.
  */
-void Vector::Mirror(const Vector *n)
+void Vector::Mirror(const Vector * const n)
 {
   double projection;
@@ -817 +816 @@
  * \return true - success, vectors are linear independent, false - failure due to linear dependency
  */
-bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2, const Vector *y3)
+bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2, const Vector * const y3)
 {
   Vector x1, x2;
@@ -853 +852 @@
  * \return true - success, vectors are linear independent, false - failure due to linear dependency
  */
-bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2)
+bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2)
 {
   Vector x1,x2;
@@ -884 +883 @@
  * \return true - success, false - vector is zero
  */
-bool Vector::MakeNormalVector(const Vector *y1)
+bool Vector::MakeNormalVector(const Vector * const y1)
 {
   bool result = false;
@@ -904 +903 @@
  * \return true - success, false - failure (null vector given)
  */
-bool Vector::GetOneNormalVector(const Vector *GivenVector)
+bool Vector::GetOneNormalVector(const Vector * const GivenVector)
 {
   int Components[NDIM]; // contains indices of non-zero components
@@ -950 +949 @@
  * \return scaling parameter for this vector
  */
-double Vector::CutsPlaneAt(Vector *A, Vector *B, Vector *C)
+double Vector::CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const
 {
 //  cout << Verbose(3) << "For comparison: ";
@@ -965 +964 @@
  * \return true if success, false if failed due to linear dependency
  */
-bool Vector::LSQdistance(Vector **vectors, int num)
+bool Vector::LSQdistance(const Vector **vectors, int num)
 {
   int j;
@@ -1047 +1046 @@
  * \param *y vector
  */
-void Vector::AddVector(const Vector *y)
+void Vector::AddVector(const Vector * const y)
 {
   for (int i=NDIM;i--;)
@@ -1056 +1055 @@
  * \param *y vector
  */
-void Vector::SubtractVector(const Vector *y)
+void Vector::SubtractVector(const Vector * const y)
 {
   for (int i=NDIM;i--;)
@@ -1065 +1064 @@
  * \param *y vector
  */
-void Vector::CopyVector(const Vector *y)
+void Vector::CopyVector(const Vector * const y)
 {
   for (int i=NDIM;i--;)
@@ -1074 +1073 @@
  * \param y vector
  */
-void Vector::CopyVector(const Vector y)
+void Vector::CopyVector(const Vector &y)
 {
   for (int i=NDIM;i--;)
@@ -1085 +1084 @@
  * \param check whether bounds shall be checked (true) or not (false)
  */
-void Vector::AskPosition(double *cell_size, bool check)
+void Vector::AskPosition(const double * const cell_size, const bool check)
 {
   char coords[3] = {'x','y','z'};
@@ -1113 +1112 @@
  * \bug this is not yet working properly
  */
-bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
+bool Vector::SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c)
 {
   double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
@@ -1276 +1275 @@
  * @param three vectors forming the matrix that defines the shape of the parallelpiped
  */
-bool Vector::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
+bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
 {
   Vector a;

src/vector.hpp

@@ -30 +30 @@
   ~Vector();
 
-  double Distance(const Vector *y) const;
-  double DistanceSquared(const Vector *y) const;
-  double DistanceToPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset);
-  double PeriodicDistance(const Vector *y, const double *cell_size) const;
-  double PeriodicDistanceSquared(const Vector *y, const double *cell_size) const;
-  double ScalarProduct(const Vector *y) const;
+  double Distance(const Vector * const y) const;
+  double DistanceSquared(const Vector * const y) const;
+  double DistanceToPlane(ofstream *out, const Vector * const PlaneNormal, const Vector * const PlaneOffset) const;
+  double PeriodicDistance(const Vector * const y, const double * const cell_size) const;
+  double PeriodicDistanceSquared(const Vector * const y, const double * const cell_size) const;
+  double ScalarProduct(const Vector * const y) const;
   double Norm() const;
   double NormSquared() const;
-  double Angle(const Vector *y) const;
+  double Angle(const Vector * const y) const;
   bool IsZero() const;
   bool IsOne() const;
-  bool IsNormalTo(const Vector *normal) const;
+  bool IsNormalTo(const Vector * const normal) const;
 
-  void AddVector(const Vector *y);
-  void SubtractVector(const Vector *y);
-  void CopyVector(const Vector *y);
-  void CopyVector(const Vector y);
-  void RotateVector(const Vector *y, const double alpha);
-  void VectorProduct(const Vector *y);
-  void ProjectOntoPlane(const Vector *y);
-  void ProjectIt(const Vector *y);
-  Vector Projection(const Vector *y) const;
+  void AddVector(const Vector * const y);
+  void SubtractVector(const Vector * const y);
+  void CopyVector(const Vector * const y);
+  void CopyVector(const Vector &y);
+  void RotateVector(const Vector * const y, const double alpha);
+  void VectorProduct(const Vector * const y);
+  void ProjectOntoPlane(const Vector * const y);
+  void ProjectIt(const Vector * const y);
+  Vector Projection(const Vector * const y) const;
   void Zero();
-  void One(double one);
-  void Init(double x1, double x2, double x3);
+  void One(const double one);
+  void Init(const double x1, const double x2, const double x3);
   void Normalize();
-  void Translate(const Vector *x);
-  void Mirror(const Vector *x);
-  void Scale(double **factor);
-  void Scale(double *factor);
-  void Scale(double factor);
-  void MatrixMultiplication(const double *M);
-  double * InverseMatrix(double *A);
-  void InverseMatrixMultiplication(const double *M);
-  void KeepPeriodic(ofstream *out, double *matrix);
-  void LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors);
-  double CutsPlaneAt(Vector *A, Vector *B, Vector *C);
-  bool GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *Origin, Vector *LineVector);
-  bool GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b, const Vector *Normal = NULL);
-  bool GetOneNormalVector(const Vector *x1);
-  bool MakeNormalVector(const Vector *y1);
-  bool MakeNormalVector(const Vector *y1, const Vector *y2);
-  bool MakeNormalVector(const Vector *x1, const Vector *x2, const Vector *x3);
-  bool SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c);
-  bool LSQdistance(Vector **vectors, int dim);
-  void AskPosition(double *cell_size, bool check);
+  void Translate(const Vector * const x);
+  void Mirror(const Vector * const x);
+  void Scale(const double ** const factor);
+  void Scale(const double * const factor);
+  void Scale(const double factor);
+  void MatrixMultiplication(const double * const M);
+  double * InverseMatrix(const double * const A);
+  void InverseMatrixMultiplication(const double * const M);
+  void KeepPeriodic(ofstream *out, const double * const matrix);
+  void LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors);
+  double CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const;
+  bool GetIntersectionWithPlane(ofstream *out, const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector);
+  bool GetIntersectionOfTwoLinesOnPlane(ofstream *out, const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *Normal = NULL);
+  bool GetOneNormalVector(const Vector * const x1);
+  bool MakeNormalVector(const Vector * const y1);
+  bool MakeNormalVector(const Vector * const y1, const Vector * const y2);
+  bool MakeNormalVector(const Vector * const x1, const Vector * const x2, const Vector * const x3);
+  bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);
+  bool LSQdistance(const Vector ** vectors, int dim);
+  void AskPosition(const double * const cell_size, const bool check);
   bool Output(ofstream *out) const;
-  bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
-  void WrapPeriodically(const double *M, const double *Minv);
+  bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
+  void WrapPeriodically(const double * const M, const double * const Minv);
 };