Index: util/src/NanoCreator.c
===================================================================
--- util/src/NanoCreator.c	(revision 1f42095177a4310820f5d7500d096a9c3a93b46d)
+++ util/src/NanoCreator.c	(revision 1f42095177a4310820f5d7500d096a9c3a93b46d)
@@ -0,0 +1,1290 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+
+#include "NanoCreator.h"
+
+// ---------------------------------- F U N C T I O N S ----------------------------------------------
+
+// ================================= File functions ==============================
+
+#define LINE_SIZE 80
+#define NDIM  3
+
+/** Allocs memory and prints a message on fail.
+ * \param *size size of alloc in bytes
+ * \param *msg error msg
+ * \return pointer to allocated memory
+ */
+void * Malloc (size_t size, const char *msg) 
+{
+  void *ptr = malloc(size);
+  if (ptr == NULL) {
+    if (msg == NULL)
+      fprintf(stdout, "ERROR: Malloc\n");
+    else
+      fprintf(stdout, "ERROR: Malloc - %s\n", msg);
+    return NULL;
+  } else {
+    return ptr;
+  }
+}
+
+/** Callocs memory and prints a message on fail.
+ * \param *size size of alloc in bytes
+ * \param *value pointer to initial value
+ * \param *msg error msg
+ * \return pointer to allocated memory
+ */
+void * Calloc (size_t size, double value, const char *msg) 
+{
+  void *ptr = calloc(size, value);
+  if (ptr == NULL) {
+    if (msg == NULL)
+      fprintf(stdout, "ERROR: Calloc\n");
+    else
+      fprintf(stdout, "ERROR: Calloc - %s\n", msg);
+    return NULL;
+  } else {
+    return ptr;
+  }
+}
+
+/** Frees memory only if ptr != NULL.
+ * \param *ptr pointer to array
+ * \param *msg
+ */
+void Free(void * ptr, const char *msg)
+{
+  if (ptr != NULL)
+    free(ptr);
+  else {
+    if (msg == NULL)
+      fprintf(stdout, "ERROR: Free\n");
+    else
+      fprintf(stdout, "ERROR: Free - %s\n", msg);
+  }
+}
+
+/** Allocs memory and prints a message on fail.
+ * \param **old_ptr pointer to old memory range 
+ * \param *newsize new size of alloc in bytes
+ * \param *msg error msg
+ * \return pointer to allocated memory
+ */
+void * Realloc (void *old_ptr, size_t newsize, const char *msg) 
+{
+  if (old_ptr == NULL) {
+    fprintf(stdout, "ERROR: Realloc -  old_ptr NULL\n");
+    exit(255);
+  }
+  void *ptr = realloc(old_ptr, newsize);
+  if (ptr == NULL) {
+    if (msg == NULL)
+      fprintf(stdout, "ERROR: Realloc\n");
+    else
+      fprintf(stdout, "ERROR: Realloc - %s\n", msg);
+    return NULL;
+  } else {
+    return ptr;
+  }
+}
+
+/** Reads a file's contents into a char buffer of appropiate size.
+ * \param *filename name of file
+ * \param pointer to integer holding read/allocated buffer length
+ * \return pointer to allocated segment with contents
+ */
+char * ReadBuffer (char *filename, int *bufferlength) 
+{
+  if ((filename == NULL) || (bufferlength == NULL)) {
+    fprintf(stderr, "ERROR: ReadBuffer - ptr to filename or bufferlength NULL\n");
+    exit(255);
+  }
+  char *buffer = Malloc(sizeof(char)*LINE_SIZE, "ReadBuffer: buffer");
+  int i = 0, number;
+  FILE *file = fopen(filename, "r");
+  if (file == NULL) {
+    fprintf(stderr, "File open error: %s", filename);
+    exit(255);
+  }
+  while ((number = fread(&buffer[i], sizeof(char), LINE_SIZE, file)) == LINE_SIZE) {
+    //fprintf(stdout, "%s", &buffer[i]);
+    i+= LINE_SIZE;
+    buffer = (char *) Realloc(buffer, i+LINE_SIZE, "ReadBuffer: buffer");
+  }
+  fclose(file);
+  fprintf(stdout, "Buffer length is %i\n", i);
+  *bufferlength = i+(number);
+  return buffer;
+}
+
+/** Extracts next line out of a buffer.
+ * \param *buffer buffer to parse for newline
+ * \param *line contains complete line on return
+ * \return length of line, 0 if end of file
+ */
+int GetNextline(char *buffer, char *line)
+{
+  if ((buffer == NULL) || (line == NULL)) {
+    fprintf(stderr, "ERROR: GetNextline - ptr to buffer or line NULL\n");
+    exit(255);
+  }
+  int length;
+  char *ptr = strchr(buffer, '\n');
+  //fprintf(stdout, "Newline at %p from %p\n", ptr, buffer);
+  if (ptr == NULL) { // buffer ends
+    return 0;
+  } else {
+    //fprintf(stdout, "length of line is %d\n", length);
+    length = (int)(ptr - buffer)/sizeof(char);
+    strncpy(line, buffer, length);
+    line[length]='\0';
+    return length+1;
+  }
+}
+
+/** Adds commentary stuff (needed for further stages) to Cell xyz files.
+ * \param *filename  file name
+ * \param atomicnumner number of atoms in xyz
+ * \param **Vector list of three unit cell vectors
+ * \param **Recivector list of three reciprocal unit cell vectors
+ * \param atomicnumber number of atoms in cell
+ */
+void AddAtomicNumber(char *filename, int atomicnumber, double **Vector, double **Recivector) 
+{
+  if ((filename == NULL) || (Vector == NULL) || (Recivector == NULL)) {
+    fprintf(stdout, "ERROR: AddAtomicNumber - ptr to filename, Vector or Recivector NULL\n");
+    exit(255);
+  }
+  int bufferlength;
+  char *buffer = ReadBuffer(filename, &bufferlength);
+  FILE *file2 = fopen(filename, "w+");
+  if (file2 == NULL) {
+    fprintf(stdout, "ERROR: AddAtomicNumber: %s can't open for writing\n", filename);
+    exit(255);
+  }
+  double volume = Determinant(Vector);
+  time_t now;
+  
+  now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
+  // open for writing and prepend
+  fprintf(file2,"%d\n", atomicnumber); // 2
+  fprintf(file2,"\tgenerated with Nanotube creator on %s", ctime(&now));
+  fwrite(buffer, sizeof(char), bufferlength, file2); // append buffer
+
+  // Add primitive vectors as comment
+  fprintf(file2,"\n****************************************\n\n");
+  fprintf(file2,"Primitive vectors\n");
+  fprintf(file2,"a(1) = %f\t%f\t%f\n", Vector[0][0], Vector[0][1], Vector[0][2]);
+  fprintf(file2,"a(2) = %f\t%f\t%f\n", Vector[1][0], Vector[1][1], Vector[1][2]);
+  fprintf(file2,"a(3) = %f\t%f\t%f\n", Vector[2][0], Vector[2][1], Vector[2][2]);
+  fprintf(file2,"\nVolume = %f", volume);
+  fprintf(file2,"\nReciprocal Vectors\n");
+  fprintf(file2,"b(1) = %f\t%f\t%f\n", Recivector[0][0], Recivector[0][1], Recivector[0][2]);
+  fprintf(file2,"b(2) = %f\t%f\t%f\n", Recivector[1][0], Recivector[1][1], Recivector[1][2]);
+  fprintf(file2,"b(3) = %f\t%f\t%f\n", Recivector[2][0], Recivector[2][1], Recivector[2][2]);
+
+  fclose(file2); // close file
+  Free(buffer, "AddAtomicNumber: buffer");
+}
+
+/** Adds commentary stuff (needed for further stages) to Sheet xyz files.
+ * \param *filename file name
+ * \param *axis array with major, minor and no axis
+ * \param *chiral pointer to array with both chiral values
+ * \param *factors pointer to array with length and radius factor
+ * \param seed random number seed
+ * \param numbercell number of atoms in unit cell, needed as length of \a *randomness 
+ * \param *randomness for each atom in unit cell a factor between 0..1, giving its probability of appearance
+ */
+void AddSheetInfo(char *filename, int *axis, int *chiral, int *factors, int seed, int numbercell, double *randomness)
+{
+  int i;
+  if ((filename == NULL) || (axis == NULL) || (chiral == NULL) || (factors == NULL)) {
+    fprintf(stdout, "ERROR: AddSheetInfo - ptr to filename, axis, chiral or factors NULL\n");
+    exit(255);
+  }
+  // open for writing and append
+  FILE *file2 = fopen(filename,"a");
+  if (file2 == NULL) {
+    fprintf(stderr, "ERROR: AddSheetInfo - can't open %s for appending\n", filename);
+    exit(255);
+  }
+  // Add primitive vectors as comment
+  fprintf(file2,"\n****************************************\n\n");
+  fprintf(file2,"Axis %d\t%d\t%d\n", axis[0], axis[1], axis[2]);
+  fprintf(file2,"(n,m) %d\t%d\n", chiral[0], chiral[1]);
+  fprintf(file2,"factors %d\t%d\n", factors[0], factors[1]);
+  fprintf(file2,"seed %d\n", seed);
+  fprintf(file2,"\nRandomness\n");
+  for (i=0; i<numbercell; i++) {
+    fprintf(file2,"%d %g\n", i, randomness[i]);
+  }
+  fclose(file2);
+}
+
+
+// ================================= Vector functions ==============================
+
+/** Transforms a vector b with a matrix A: Ab = x.
+ * \param **matrixref reference to NDIMxNDIM matrix A
+ * \param *vectorref reference to NDIM vector b
+ * \return reference to resulting NDIM vector Ab = x
+ */
+double *MatrixTrafo(double **matrixref, double *vectorref)
+{
+  if ((matrixref == NULL) || (vectorref == NULL)) {
+    fprintf(stderr, "ERROR: MatrixTrafo: ptr to matrixref or vectorref NULL\n");
+    exit(255);
+  }
+  //double *returnvector = Calloc(sizeof(double)*NDIM, 0., "MatrixTrafo: returnvector");
+  double *returnvector = calloc(sizeof(double)*NDIM, 0.);
+  if (returnvector == NULL) {
+    fprintf(stderr, "ERROR: MatrixTrafo - returnvector\n");
+    exit(255);
+  }
+  int i,j;
+  
+  for (i=0;i<NDIM;i++)
+    for (j=0;j<NDIM;j++)
+      returnvector[j] += matrixref[i][j] * vectorref[i];
+      
+  return returnvector;  
+}
+double *MatrixTrafoInverse(double *vectorref, double **matrixref)
+{
+  if ((matrixref == NULL) || (vectorref == NULL)) {
+    fprintf(stderr, "ERROR: MatrixTrafo: ptr to matrixref or vectorref NULL\n");
+    exit(255);
+  }
+  //double *returnvector = Calloc(sizeof(double)*NDIM, 0., "MatrixTrafo: returnvector");
+  double *returnvector = calloc(sizeof(double)*NDIM, 0.);
+  if (returnvector == NULL) {
+    fprintf(stderr, "ERROR: MatrixTrafo - returnvector\n");
+    exit(255);
+  }
+  int i,j;
+  
+  for (i=0;i<NDIM;i++)
+    for (j=0;j<NDIM;j++)
+      returnvector[i] += matrixref[i][j] * vectorref[j];
+      
+  return returnvector;  
+}
+
+/** Inverts a NDIMxNDIM matrix.
+ * \param **matrix to be inverted
+ * \param **inverse allocated space for inverse of \a **matrix
+ */
+void MatrixInversion(double **matrix, double **inverse) 
+{
+  if ((matrix == NULL) || (inverse == NULL)) {
+    fprintf(stderr, "ERROR: MatrixInversion: ptr to matrix or inverse NULL\n");
+    exit(255);
+  }
+  // determine inverse
+  double det = Determinant(matrix);
+  inverse[0][0] = (matrix[1][1]*matrix[2][2] - matrix[1][2]*matrix[2][1])/det;  
+  inverse[1][0] = (matrix[0][2]*matrix[2][1] - matrix[0][1]*matrix[2][2])/det;  
+  inverse[2][0] = (matrix[0][1]*matrix[1][2] - matrix[0][2]*matrix[1][1])/det;  
+  inverse[0][1] = (matrix[1][2]*matrix[2][0] - matrix[1][0]*matrix[2][2])/det;  
+  inverse[1][1] = (matrix[0][0]*matrix[2][2] - matrix[0][2]*matrix[2][0])/det;  
+  inverse[2][1] = (matrix[0][2]*matrix[1][0] - matrix[0][0]*matrix[1][2])/det;  
+  inverse[0][2] = (matrix[1][0]*matrix[2][1] - matrix[1][1]*matrix[2][0])/det;  
+  inverse[1][2] = (matrix[0][1]*matrix[2][0] - matrix[0][0]*matrix[2][1])/det;  
+  inverse[2][2] = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0])/det;
+   
+  // check inverse
+  int flag = 0; 
+  int i,j,k;
+  double tmp;   
+  fprintf(stdout, "Checking inverse ... ");
+  for (i=0;i<NDIM;i++)
+    for (j=0;j<NDIM;j++) {
+      tmp = 0.;
+      for (k=0;k<NDIM;k++)
+        tmp += matrix[i][k]*inverse[j][k];
+      if (!flag) {
+        if (i == j) {
+          flag = (fabs(1.-tmp) > MYEPSILON) ? 1 : 0;
+        } else {
+          flag = (fabs(tmp) > MYEPSILON) ? 1 : 0;
+        }        
+      }
+    }
+  if (!flag)
+    fprintf(stdout, "ok\n");
+  else
+    fprintf(stdout, "false\n");
+}
+
+/** Flips to double numbers in place.
+ * \param *number1 pointer to first double
+ * \param *number2 pointer to second double
+ */
+void flip(double *number1, double *number2)
+{
+  if ((number1 == NULL) || (number2 == NULL)) {
+    fprintf(stderr, "ERROR: flip: ptr to number1 or number2 NULL\n");
+    exit(255);
+  }
+  double tmp = *number1;
+  *number1 = *number2;
+  *number2 = tmp;
+}
+
+/** Transposes a matrix.
+ * \param **matrix pointer to NDIMxNDIM-matrix array
+ */
+void Transpose(double **matrix)
+{
+  if (matrix == NULL) {
+    fprintf(stderr, "ERROR: Transpose: ptr to matrix NULL\n");
+    exit(255);
+  }
+  int i,j;
+  for (i=0;i<NDIM;i++)
+    for (j=0;j<i;j++)
+      flip(&matrix[i][j],&matrix[j][i]);
+}
+ 
+
+/** Computes the determinant of a NDIMxNDIM matrix.
+ * \param **matrix pointer to matrix array
+ * \return det(matrix)
+ */
+double Determinant(double **matrix) {
+  if (matrix == NULL) {
+    fprintf(stderr, "ERROR: Determinant: ptr to Determinant NULL\n");
+    exit(255);
+  }
+  double det =   matrix[0][0] * (matrix[1][1]*matrix[2][2] - matrix[1][2]*matrix[2][1])
+               - matrix[1][1] * (matrix[0][0]*matrix[2][2] - matrix[0][2]*matrix[2][0])
+               + matrix[2][2] * (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);  
+  return det;
+}
+
+/** Adds \a *vector1 onto \a *vector2 coefficient-wise.
+ * \param *vector1 first vector, on return contains sum
+ * \param *vector2 vector which is projected
+ * \return sum of the two vectors
+ */
+double * VectorAdd(double *vector1, double *vector2)
+{
+  if ((vector1 == NULL) || (vector2 == NULL)) {
+    fprintf(stderr, "ERROR: VectorAdd: ptr to vector1 or vector2 NULL\n");
+    exit(255);
+  }
+  //double *returnvector = Calloc(sizeof(double)*NDIM, 0., "VectorAdd: returnvector");
+  double *returnvector = calloc(sizeof(double)*NDIM, 0.);
+  if (returnvector == NULL) {
+    fprintf(stderr, "ERROR: VectorAdd - returnvector\n");
+    exit(255);
+  }
+  int i;
+  
+  for (i=0;i<NDIM;i++)
+    returnvector[i] = vector1[i] + vector2[i];
+    
+  return returnvector;
+}
+
+/**  Fixed GramSchmidt-Orthogonalization for NDIM vectors
+ * \param @orthvector reference to NDIMxNDIM matrix
+ * \param @orthbetrag reference to NDIM vector with vector magnitudes
+ * \param @axis major-, minor- and noaxis for specific order for the GramSchmidt procedure
+ */
+void Orthogonalize(double **orthvector, int *axis)
+{
+  if ((orthvector == NULL) || (axis == NULL)) {
+    fprintf(stderr, "ERROR: Orthogonalize: ptr to orthvector or axis NULL\n");
+    exit(255);
+  }
+  double betrag1, betrag2;
+  int i;
+  
+  // first vector is untouched
+  // second vector
+  betrag1 = Projection(orthvector[axis[1]], orthvector[axis[0]]); 
+  for (i=0;i<NDIM;i++)
+     orthvector[axis[0]][i] -= orthvector[axis[1]][i] * betrag1;
+  // third vector
+  betrag1 = Projection(orthvector[axis[0]], orthvector[axis[2]]); 
+  betrag2 = Projection(orthvector[axis[1]], orthvector[axis[2]]); 
+  for (i=0;i<NDIM;i++)
+     orthvector[axis[2]][i] -= orthvector[axis[0]][i] * betrag1 + orthvector[axis[1]][i] * betrag2;
+}
+
+/** Computes projection of \a *vector2 onto \a *vector1.
+ * \param *vector1 reference vector
+ * \param *vector2 vector which is projected
+ * \return projection 
+ */
+double Projection(double *vector1, double *vector2)
+{
+  if ((vector1 == NULL) || (vector2 == NULL)) {
+    fprintf(stderr, "ERROR: Projection: ptr to vector1 or vector2 NULL\n");
+    exit(255);
+  }
+  return (ScalarProduct(vector1, vector2)/Norm(vector1)/Norm(vector2));
+}
+
+/** Determine scalar product between two vectors.
+ * \param *vector1 first vector
+ * \param *vector2 second vector
+ * \return scalar product
+ */
+double ScalarProduct(double *vector1, double *vector2)
+{
+  if ((vector1 == NULL) || (vector2 == NULL)) {
+    fprintf(stderr, "ERROR: ScalarProduct: ptr to vector1 or vector2 NULL\n");
+    exit(255);
+  }
+  double scp = 0.;
+  int i;
+  
+  for (i=0;i<NDIM;i++)
+    scp += vector1[i] * vector2[i];
+    
+  return scp;
+}
+
+/** Computes norm of \a *vector.
+ * \param *vector pointer to NDIM vector
+ * \return norm of \a *vector
+ */
+double Norm(double *vector) 
+{
+  if (vector == NULL) {
+    fprintf(stderr, "ERROR: Norm: ptr to vector NULL\n");
+    exit(255);
+  }
+  return sqrt(ScalarProduct(vector, vector));
+}
+
+/** Prints vector to \a *file.
+ * \param *file file or e.g. stdout
+ * \param *vector vector to be printed
+ */
+void PrintVector(FILE *file, double *vector)
+{
+  if ((file == NULL) || (vector == NULL)) {
+    fprintf(stderr, "ERROR: PrintVector: ptr to file or vector NULL\n");
+    exit(255);
+  }
+  int i;
+  for (i=0;i<NDIM;i++)
+    fprintf(file, "%5.5f\t", vector[i]);
+  fprintf(file, "\n");
+}
+
+/** Prints matrix to \a *file.
+ * \param *file file or e.g. stdout
+ * \param **matrix matrix to be printed
+ */
+void PrintMatrix(FILE *file, double **matrix)
+{
+  if ((file == NULL) || (matrix == NULL)) {
+    fprintf(stderr, "ERROR: PrintMatrix: ptr to file or matrix NULL\n");
+    exit(255);
+  }
+  int i,j;
+  for (i=0;i<NDIM;i++) {
+    for (j=0;j<NDIM;j++)
+      fprintf(file, "%5.5f\t", matrix[i][j]);
+    fprintf(file, "\n");
+  }
+}
+
+/** Returns greatest common denominator.
+ * \param a first integer
+ * \param b second integer
+ * \return GCD of a and b
+ */
+int GCD(int a, int b)
+{
+  int c;
+  do {
+    c = a % b;        /* Rest of integer divison */
+    a = b; b = c;              /* flip the two values */
+  } while( c != 0);
+  return a;
+}
+
+/** Determines the biggest diameter of a sheet.
+ * \param **matrix reference to NDIMxNDIM matrix with row vectors
+ * \param *axis reference to NDIM vector with permutation of axis indices [0,1,2]
+ * \param *factors factorsfor axis[0] and axis[1]
+ * \return biggest diameter of sheet
+*/
+double DetermineBiggestDiameter(double **matrix, int *axis, int *factors)
+{
+  if ((axis == NULL) || (factors == NULL) || (matrix == NULL)) {
+    fprintf(stderr, "ERROR: DetermineBiggestDiameter: ptr to factors, axis or matrix NULL\n");
+    exit(255);
+  }
+  double diameter[2] = {0.,0.};
+  int i, biggest;
+  
+  for (i=0;i<NDIM;i++) {
+    diameter[0] += (matrix[axis[0]][i]*factors[0] - matrix[axis[1]][i]*factors[1]) * (matrix[axis[0]][i]*factors[0] - matrix[axis[1]][i]*factors[1]);
+    diameter[1] += (matrix[axis[0]][i]*factors[0] + matrix[axis[1]][i]*factors[1]) * (matrix[axis[0]][i]*factors[0] + matrix[axis[1]][i]*factors[1]);
+  }
+  if ((diameter[0] - diameter[1]) > MYEPSILON) {
+    biggest = 0;
+  } else {
+    biggest = 1;
+  }
+  diameter[0] = sqrt(diameter[0]);
+  diameter[1] = sqrt(diameter[1]);
+  fprintf(stdout, "\n\nMajor diameter of the sheet is %5.5f, minor diameter is %5.5f.\n",diameter[biggest],diameter[(biggest+1)%2]);
+
+  return diameter[biggest];
+} 
+
+/** Determines the center of gravity of atoms in a buffer \a bufptr with given \a number
+ * \param *bufptr pointer to char buffer with atoms in (name x y z)-manner
+ * \param number number of atoms/lines to scan
+ * \return NDIM vector (allocated doubles) pointing back to center of gravity
+ */
+double * CenterOfGravity(char *bufptr, int number)
+{
+  if (bufptr == NULL) {
+    fprintf(stderr, "ERROR: CenterOfGravity - bufptr NULL\n");
+    exit(255);
+  }
+  double *cog = calloc(sizeof(double)*NDIM, 0.);
+  if (cog == NULL) {
+    fprintf(stderr, "ERROR: CenterOfGravity - cog\n");
+    exit(255);
+  }
+  double *atom = Malloc(sizeof(double)*NDIM, "CenterOfGravity: atom");
+  char name[255], line[255];
+  int i,j;
+  
+  // Determine center of gravity
+  for (i=0;i<number;i++) {
+    bufptr += GetNextline(bufptr, line);
+    sscanf(line, "%s %lg %lg %lg", name, &atom[0], &atom[1], &atom[2]);
+    //fprintf(stdout, "Read Atom %s %lg %lg %lg\n", name, atom[0], atom[1], atom[2]);
+    for (j=0;j<NDIM;j++)
+      cog[j] += atom[j];
+  }
+  for (i=0;i<NDIM;i++)
+    cog[i] /= -number;
+    
+  Free(atom, "CenterOfGravity: atom");
+  return cog;
+}
+
+// ================================= other functions ==============================
+
+void Debug(char *msg)
+{
+  if (msg == NULL) {
+    fprintf(stderr, "ERROR: Debug: ptr to msg NULL\n");
+    exit(255);
+  }
+  fprintf(stdout, "%s", msg);
+}
+
+// --------------------------------------- M A I N ---------------------------------------------------
+int main(int argc, char **argv) {
+	char *filename = NULL; 
+  char *CellFilename = NULL, *SheetFilename = NULL, *TubeFilename = NULL, *TorusFilename = NULL;
+  char *SheetFilenameAligned = NULL, *TubeFilenameAligned = NULL;
+	double **Vector, **Recivector = NULL, **Tubevector = NULL, **TubevectorInverse = NULL;
+  double *atom = NULL, *atom_transformed = NULL;
+  double *x = NULL, *coord = NULL, *tempvector = NULL, *offset = NULL;
+  //double *cog = NULL, *cog_projected = NULL;
+  char stage[6];
+  int axis[NDIM] = {0,1,2}, chiral[2] = {1,1}, factors[NDIM] = {1,1,1};
+  char name[255], line[255], input = 'n';
+  char *CellBuffer = NULL, *SheetBuffer = NULL, *TubeBuffer = NULL, *bufptr = NULL;
+  double *randomness = NULL, percentage; // array with percentages for presence in sheet and beyond
+  int i,j, ggT;
+  int length;
+  
+	// Read command line arguments	
+	if (argc <= 2) {
+		fprintf(stdout, "Usage: %s <file> <stage>\n\tWhere <file> specifies a file to start from <stage> or a basename\n\t<stage> is either None, Cell, Sheet, Tube, Torus and specifies where to start the rolling up from.\n\tNote: The .Aligned. files can't be used (rotation is essential).\n", argv[0]);
+		exit(255);
+	} else {
+		// store variables
+		filename = argv[1];
+		strncpy(stage, argv[2], 5);
+    fprintf(stdout, "I will begin at stage %s.\n", stage);
+		
+		// build filenames
+    char *ptr = strrchr(filename, '.');
+    length = strlen(filename);
+    if (ptr != NULL) {
+      length = ((int)(ptr-filename) >= length-5) ? (int)(ptr-filename) : length;  
+      filename[length] = '\0'; // eventueller
+    }
+    fprintf(stdout, "I will use \'%s' as base for the filenames.\n\n", filename);
+    CellFilename = Malloc(sizeof(char)*(length+10), "Main: CellFilename");
+    SheetFilename = Malloc(sizeof(char)*(length+11), "Main: SheetFilename");
+    TubeFilename = Malloc(sizeof(char)*(length+10), "Main: TubeFilename");
+    TorusFilename = Malloc(sizeof(char)*(length+11), "Main: TorusFilename");
+    SheetFilenameAligned = Malloc(sizeof(char)*(length+20), "Main: SheetFilenameAligned");
+    TubeFilenameAligned = Malloc(sizeof(char)*(length+19), "Main: TubeFilenameAligned");
+		sprintf(CellFilename, "%s.Cell.xyz", filename);
+    sprintf(SheetFilename, "%s.Sheet.xyz", filename);
+    sprintf(TubeFilename, "%s.Tube.xyz", filename);
+    sprintf(TorusFilename, "%s.Torus.xyz", filename);
+    sprintf(SheetFilenameAligned, "%s.Sheet.Aligned.xyz", filename);
+    sprintf(TubeFilenameAligned, "%s.Tube.Aligned.xyz", filename);
+	}
+	
+  // Allocating memory
+  Debug ("Allocating memory\n");
+  atom = (double *) Malloc(sizeof(double)*NDIM, "Main: atom");
+  Vector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Vector");
+  Recivector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Recivector");
+  Tubevector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Tubevector");
+  TubevectorInverse = (double **) Malloc(sizeof(double *)*NDIM, "Main: *TubevectorInverse");
+  for (i=0; i<NDIM; i++ )  {
+    Vector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Vector");
+    Recivector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Recivector");
+    Tubevector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Tubevector");
+    TubevectorInverse[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: TubevectorInverse");
+  }
+  
+  // ======================== STAGE: Cell ==============================
+  // The cell is simply created by transforming relative coordinates within the cell
+  // into cartesian ones using the unit cell vectors.
+    
+  double volume;
+  int numbercell;
+  FILE *CellFile;
+
+  Debug ("STAGE: None\n");
+  // Read cell vectors from stdin or from file
+	if (!strncmp(stage, "Non", 3)) {
+  	fprintf(stdout, "You will give the unit cell of the given substance.\nAfterwards, the programme will create a Sheet, a Tube and a Torus, each with their own xyz-file named accordingly.\n\n");
+    fprintf(stdout, "Enter 1st primitive vector: ");
+    fscanf(stdin, "%lg %lg %lg", &Vector[0][0], &Vector[0][1], &Vector[0][2]);	
+    fprintf(stdout, "Enter 2nd primitive vector: ");
+    fscanf(stdin, "%lg %lg %lg", &Vector[1][0], &Vector[1][1], &Vector[1][2]);  
+    fprintf(stdout, "Enter 3rd primitive vector: ");
+    fscanf(stdin, "%lg %lg %lg", &Vector[2][0], &Vector[2][1], &Vector[2][2]);
+    fprintf(stdout,"Unit vectors are\n");
+    PrintMatrix(stdout, Vector);  
+	} else {
+    char *ptr = NULL;
+    char dummy[10];
+    CellBuffer = bufptr = ReadBuffer(CellFilename, &length);
+    for (i=0;i<NDIM;i++) {
+      sprintf(dummy, "a(%i) = ", i+1);
+      fprintf(stdout, "%s", dummy);
+      while ((length = GetNextline(bufptr, line)) != -1) {
+        bufptr += (length)*sizeof(char);
+        //fprintf(stdout, "LINE at %p: %s\n", bufptr, line);
+        if ((ptr = strstr(line, dummy)) != NULL)
+          break;
+      }
+      ptr += strlen(dummy);
+      sscanf(ptr, "%lg %lg %lg", &Vector[i][0], &Vector[i][1], &Vector[i][2]);
+      fprintf(stdout, "%5.5lg %5.5lg %5.5lg\n", Vector[i][0], Vector[i][1], Vector[i][2]); 
+    }
+	}
+
+  volume = Determinant(Vector);
+  fprintf(stdout,"Volume is %lg\n", volume);
+  MatrixInversion(Vector, Recivector);  
+  //Transpose(Recivector);   // inverse's got row vectors if normal matrix' got column ones
+  fprintf(stdout, "Reciprocal vector is ");
+  PrintMatrix(stdout, Recivector);
+  fprintf(stdout, "Reciprocal volume is %lg\n", Determinant(Recivector));
+
+  fprintf(stdout, "Vector magnitudes: %5.5lg %5.5lg %5.5lg\n", Norm(Vector[0]), Norm(Vector[1]), Norm(Vector[2]));
+   
+  Debug ("STAGE: Cell\n");
+  if (!strncmp(stage, "Non", 3)) {
+    fprintf(stdout, "\nHow many atoms are in the unit cell: ");
+    fscanf(stdin, "%i", &numbercell);
+    CellFile = fopen(CellFilename, "w");
+    if (CellFile == NULL) {
+      fprintf(stderr, "ERROR: main - can't open %s for writing\n", CellFilename);
+      exit(255);
+    }
+    fprintf(stdout, "\nNext, you have to enter each atom in the cell as follows, e.g.\n");
+    fprintf(stdout, "Enter \'ChemicalSymbol X Y Z\' (relative to primitive vectors): C 0.5 0.25 0.5\n\n");
+    for (i = 0; i < numbercell; i++) {
+      fprintf(stdout, "Enter for atom %i \'ChemicalSymbol X Y Z\': ", i+1);
+      fscanf(stdin, "%s %lg %lg %lg", name, &atom[0], &atom[1], &atom[2]);
+      tempvector = MatrixTrafo(Vector, atom);
+      fprintf(stdout, "Atom %i: %s %5.5lg %5.5lg %5.5lg\n", i, name, tempvector[0], tempvector[1], tempvector[2]);
+      fprintf(stdout, "Probe: %s %5.5lg %5.5lg %5.5lg\n", name,
+        atom[0]*Vector[0][0]+atom[1]*Vector[1][0]+atom[2]*Vector[2][0],      
+        atom[0]*Vector[0][1]+atom[1]*Vector[1][1]+atom[2]*Vector[2][1],      
+        atom[0]*Vector[0][2]+atom[1]*Vector[1][2]+atom[2]*Vector[2][2]    
+        );
+      fprintf(CellFile, "%s %lg %lg %lg\n", name, tempvector[0], tempvector[1], tempvector[2]);
+      Free(tempvector, "Main: At stage Cell - tempvector");
+    }
+    fflush(CellFile);
+    fclose(CellFile);
+    AddAtomicNumber(CellFilename, numbercell, Vector, Recivector);
+
+    CellBuffer = ReadBuffer(CellFilename, &length);
+        
+    sprintf(stage, "Cell");
+  } else {
+    bufptr = CellBuffer;
+    GetNextline(bufptr, line);
+    sscanf(line, "%i", &numbercell);
+  }
+ 
+  fprintf(stdout, "\nThere are %i atoms in the cell.\n", numbercell);
+    
+  // ======================== STAGE: Sheet =============================
+  // The sheet is a bit more complex. We read the cell in cartesian coordinates
+  // from the file. Next, we have to rotate the unit cell vectors by the so called
+  // chiral angle. This gives a slanted and elongated section upon the sheet of
+  // periodically repeated original unit cells. It only matches up if the factors
+  // were all integer! (That's why the rotation is discrete and the chiral angle
+  // specified not as (cos alpha, sin alpha) but as (n,m)) Also, we want this
+  // section to be rectangular, thus we orthogonalize the original unit vectors
+  // to gain our (later) tube axis.
+  // By looking at the biggest possible diameter we know whose original cells to
+  // look at and check if their respective compounds (contained atoms) still reside
+  // in the rotated, elongated section we need for the later tube.
+  // Then in a for loop we go through every cell. By projecting the vector leading
+  // from the origin to the specific atom down onto the major and minor axis we
+  // know if it's still within the boundaries spanned by these rotated and elongated
+  // (radius-, length factor) unit vectors or not. If yes, its coordinates are
+  // written to file.
+
+  int numbersheet, biggestdiameter, sheetnr[NDIM], tmp, seed;
+  double x1,x2,x3, angle;
+  char flag = 'n';
+  FILE *SheetFile = NULL;
+  FILE *SheetFileAligned = NULL;
+  
+  Debug ("STAGE: Sheet\n");
+  if (!strncmp(stage, "Cell", 4)) {
+    fprintf(stdout, "\nEnter seed unequal 0 if any of the bonds shall have a randomness in their being: ");
+    fscanf(stdin, "%d", &seed);
+    if (seed != 0) 
+      input = 'y';
+    randomness = (double *) Malloc(sizeof(double)*numbercell, "Main: at sheet - randomness");
+    for(i=0;i<numbercell;i++)
+      randomness[i] = 0.;
+    i = 0;
+    fprintf(stdout, "\n");
+    while (input == 'y') {
+      fprintf(stdout, "Enter atom number (-1  0 to end) and percentage (0.0...1.0): ");
+      fscanf(stdin, "%d %lg", &i, &percentage);
+      if (i == -1) { input = 'n'; fprintf(stdout, "Breaking\n"); }
+      else { randomness[i] = 1.-percentage; }
+    };
+    
+    fprintf(stdout, "\nSpecify the axis permutation that is going to be perpendicular to the sheet [tubeaxis, torusaxis, noaxis]: ");
+    fscanf(stdin, "%d %d %d", &axis[0], &axis[1], &axis[2]);
+    fprintf(stdout, "axis: %d %d %d\n", axis[0], axis[1], axis[2]);
+    do {      
+      fprintf(stdout, "\nNow specify the two natural numbers (m n) defining the chiral angle, \nif the result is crap, try flipping to (m,n): ");
+      fscanf(stdin, "%d %d", &chiral[0], &chiral[1]);
+      ggT = GCD(2*chiral[1]+chiral[0],2*chiral[0]+chiral[1]);
+      fprintf(stdout, "Greatest Common Denominator of (2n+m, 2m+n) is %d\n", ggT);
+      fprintf(stdout, "chiral0: %d\tchiral1: %d\n", chiral[0], chiral[1]);
+      for (i=0;i<NDIM;i++) {
+        //Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i];
+        Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
+        //Tubevector[axis[0]][i] = (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[0]][i] + (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[1]][i];
+        Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
+        //Tubevector[axis[1]][i] = (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[0]][i] + (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[1]][i];
+//        fprintf(stderr, "Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i]\n = %lg * %lg + %lg * %lg = %lg + %lg = %lg\n\n", 
+//          (double)chiral[0], Vector[axis[0]][i], (double)chiral[1], Vector[axis[1]][i], 
+//          (double)chiral[0] * Vector[axis[0]][i], (double)chiral[1] * Vector[axis[1]][i], 
+//          Tubevector[axis[0]][i]);
+        Tubevector[axis[2]][i] = Vector[axis[2]][i];
+      }
+/*      fprintf(stdout, "Vector\n");
+      PrintMatrix(stdout, Vector);
+      fprintf(stdout, "Tubevector\n");
+      PrintMatrix(stdout, Tubevector);
+      for (i=0;i<NDIM;i++) {
+        fprintf(stdout, "Tubevector %d in Unit cell vectors:\t", axis[i]);
+        tempvector = MatrixTrafoInverse(Tubevector[axis[i]], Recivector);
+        PrintVector(stdout, tempvector);
+        Free(tempvector, "Main:tempvector");
+      }*/
+
+      // Give info for length and radius factors
+      fprintf(stdout, "\nThe chiral angle then is %lg degrees with tube radius %5.5f A and length %5.5f A, i.e. torus radius of %5.5f A.\n",
+        acos(Projection(Vector[axis[0]], Tubevector[axis[0]]))/M_PI*180.,
+        Norm(Tubevector[axis[0]])/(2.*M_PI),
+        Norm(Tubevector[axis[1]]),
+        Norm(Tubevector[axis[1]])/(2.*M_PI)
+        );
+      fprintf(stdout, "\nGive integer factors for length and radius of tube (multiple of %d suggested) : ", ggT);
+      fscanf(stdin, "%d %d", &factors[1], &factors[0]);
+      fprintf(stdout, "\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. torus radius of %5.5f A.\n",
+        acos(Projection(Vector[axis[0]], Tubevector[axis[0]]))/M_PI*180.,
+        (double)factors[0]*Norm(Tubevector[axis[0]])/(2.*M_PI),
+        (double)factors[1]*Norm(Tubevector[axis[1]]),
+        (double)factors[1]*Norm(Tubevector[axis[1]])/(2.*M_PI)
+        );
+      fprintf(stdout, "Satisfied? [yn] ");
+      fscanf(stdin, "%c", &flag);
+      fscanf(stdin, "%c", &flag);
+    } while (flag != 'y');
+  } else {
+    char *ptr = NULL;
+    char dummy[10];
+    double dummydouble;
+
+    SheetBuffer = bufptr = ReadBuffer(SheetFilename, &length);
+    bufptr += (GetNextline(bufptr, line))*sizeof(char);
+    sscanf(line, "%d", &numbersheet);
+    
+    // retrieve axis permutation
+    sprintf(dummy,  "Axis");
+    fprintf(stdout, "%s ", dummy);
+    while ((length = GetNextline(bufptr, line)) != 0) {
+      bufptr += (length)*sizeof(char);
+      if ((ptr = strstr(line, dummy)) != NULL)
+        break;
+    }
+    if (length == 0) {
+      fprintf(stderr, "ERROR: Main at stage Sheet - could not find %s in %s\n", dummy, SheetFilename);
+      exit(255);
+    }
+    ptr += strlen(dummy);
+    sscanf(ptr, "%d %d %d", &axis[0], &axis[1], &axis[2]);
+    fprintf(stdout, "%d %d %d\n", axis[0], axis[1], axis[2]); 
+    
+    // retrieve chiral numbers
+    sprintf(dummy,  "(n,m)");
+    fprintf(stdout, "%s ", dummy);
+    while ((length = GetNextline(bufptr, line)) != 0) {
+      bufptr += (length)*sizeof(char);
+      if ((ptr = strstr(line, dummy)) != NULL)
+        break;
+    }
+    if (length == 0) {
+      fprintf(stderr, "ERROR: Main at stage Sheet - could not find %s in %s\n", dummy, SheetFilename);
+      exit(255);
+    }      
+    ptr += strlen(dummy);
+    sscanf(ptr, "%d %d", &chiral[0], &chiral[1]);
+    fprintf(stdout, "%d %d\n", chiral[0], chiral[1]); 
+    ggT = GCD(2*chiral[1]+chiral[0],2*chiral[0]+chiral[1]);
+    fprintf(stdout, "Greatest Common Denominator of (2n+m, 2m+n) is %d\n", ggT);
+    
+    // retrieve length and radius factors
+    sprintf(dummy,  "factors");
+    fprintf(stdout, "%s ", dummy);
+    while ((length = GetNextline(bufptr, line)) != 0) {
+      bufptr += (length)*sizeof(char);
+      if ((ptr = strstr(line, dummy)) != NULL)
+        break;
+    }
+    if (length == 0) {
+      fprintf(stderr, "ERROR: Main at stage Sheet - could not find %s in %s\n", dummy, SheetFilename);
+      exit(255);
+    }
+    ptr += strlen(dummy);
+    sscanf(ptr, "%d %d %d", &factors[0], &factors[1], &factors[2]);
+    fprintf(stdout, "%d %d %d\n", factors[0], factors[1], factors[2]); 
+
+    // create Tubevectors
+    for (i=0;i<NDIM;i++) {
+      Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
+      //Tubevector[axis[0]][i] = (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[0]][i] + (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[1]][i];
+      Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
+      //Tubevector[axis[1]][i] = (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[0]][i] + (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[1]][i];
+      Tubevector[axis[2]][i] = Vector[axis[2]][i];
+    }
+
+    // retrieve seed ... 
+    randomness = (double *) Calloc(sizeof(double)*numbercell, 0., "Main: at sheet - randomness");
+    sprintf(dummy,  "seed");
+    fprintf(stdout, "%s ", dummy);
+    while ((length = GetNextline(bufptr, line)) != 0) {
+      bufptr += (length)*sizeof(char);
+      if ((ptr = strstr(line, dummy)) != NULL)
+        break;
+    }
+    if (length == 0) {
+      fprintf(stderr, "ERROR: Main at stage Sheet - could not find %s in %s\n", dummy, SheetFilename);
+      exit(255);
+    }      
+    ptr += strlen(dummy);
+    sscanf(ptr, "%d", &seed);
+    fprintf(stdout, "%d\n", seed); 
+
+    // ... and randomness
+    if (seed != 0) {  // only parse for values if a seed, i.e. randomness wanted, was specified
+      sprintf(dummy,  "Randomness");
+      fprintf(stdout, "%s\n", dummy);
+      while ((length = GetNextline(bufptr, line)) != 0) {
+        bufptr += (length)*sizeof(char);
+        if ((ptr = strstr(line, dummy)) != NULL)
+          break;
+      }
+      if (length == 0) {
+        fprintf(stderr, "ERROR: Main at stage Sheet - could not find %s in %s\n", dummy, SheetFilename);
+        exit(255);
+      }
+      sprintf(dummy,  "probability values");
+      for (i=0;i<numbercell;i++) {
+        length = GetNextline(bufptr, line);
+        if (length == 0) {
+          fprintf(stderr, "ERROR: Main at stage Sheet - could not find %s in %s\n", dummy, SheetFilename);
+          exit(255);
+        }
+        bufptr += (length)*sizeof(char);
+        sscanf(line, "%d %lg", &j, &dummydouble);
+        randomness[j] = dummydouble;
+        fprintf(stdout, "%d %g\n", j, randomness[j]);      
+      }
+    }
+  }
+
+  //Orthogonalize(Tubevector,axis);
+  angle = acos(Projection(Vector[axis[0]], Vector[axis[1]])); // calcs angle between shanks in unit cell
+  fprintf(stdout, "The basic angle between the two shanks of the unit cell is %lg\n", angle/M_PI*180.);
+  angle = acos(Projection(Tubevector[axis[0]], Tubevector[axis[1]])); // calcs angle between shanks in unit cell
+  fprintf(stdout, "The basic angle between the two shanks of the tube unit cell is %lg\n", angle/M_PI*180.);
+  //angle -= acos(Projection(Vector[axis[0]], Tubevector[axis[0]]));
+  //angle = 30./180.*M_PI - acos(Projection(Vector[axis[0]], Tubevector[axis[0]]));
+  angle = - acos(Projection(Tubevector[axis[0]], Vector[axis[0]])); 
+  fprintf(stdout, "The relative alignment rotation angle then is %lg\n", angle/M_PI*180.);
+  angle += acos(Projection(Vector[axis[0]], Vector[axis[1]]))/2.; 
+  fprintf(stdout, "The absolute alignment rotation angle then is %lg\n", angle/M_PI*180.);
+  fprintf(stdout, "\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. torus radius of %5.5f A.\n",
+    acos(Projection(Vector[axis[0]], Tubevector[axis[0]]))/M_PI*180.,
+    (double)factors[0]*Norm(Tubevector[axis[0]])/(2.*M_PI),
+    (double)factors[1]*Norm(Tubevector[axis[1]]),
+    (double)factors[1]*Norm(Tubevector[axis[1]])/(2.*M_PI)
+    );
+  Orthogonalize(Tubevector, axis);   // with correct translational vector, not needed anymore (? what's been done here. Hence, re-inserted)
+  fprintf(stdout, "Tubevector magnitudes: %5.5lg %5.5lg %5.5lg\n", Norm(Tubevector[0]), Norm(Tubevector[1]), Norm(Tubevector[2]));
+  fprintf(stdout, "Tubebvectors are \n");
+  PrintMatrix(stdout, Tubevector);
+  MatrixInversion(Tubevector, TubevectorInverse);
+  //Transpose(TubevectorInverse);
+      fprintf(stdout, "Vector\n");
+      PrintMatrix(stdout, Vector);
+      fprintf(stdout, "TubevectorInverse\n");
+      PrintMatrix(stdout, TubevectorInverse);
+      for (i=0;i<NDIM;i++) {
+        fprintf(stdout, "Vector %d in TubeVectorInverse vectors:\t", axis[i]);
+        tempvector = MatrixTrafoInverse(Vector[axis[i]], TubevectorInverse);
+        PrintVector(stdout, tempvector);
+        Free(tempvector, "Main:tempvector");
+      }
+  fprintf(stdout, "Reciprocal Tubebvectors are \n");
+  PrintMatrix(stdout, TubevectorInverse);
+  fprintf(stdout, "Tubevector magnitudes: %5.5lg %5.5lg %5.5lg\n", Norm(Tubevector[0]), Norm(Tubevector[1]), Norm(Tubevector[2]));
+
+  biggestdiameter = DetermineBiggestDiameter(Tubevector, axis, factors);
+  for (i=0;i<NDIM;i++) {
+    sheetnr[i] = 0;
+  }
+  for (i=0;i<NDIM;i++) {
+    for (j=0;j<NDIM;j++) {
+//      sheetnr[j] =  ceil(biggestdiameter/Norm(Vector[j]));
+      if (fabs(Vector[i][j]) > MYEPSILON) {
+        tmp = ceil(biggestdiameter/fabs(Vector[i][j]));
+      } else {
+        tmp = 0;
+      }
+      sheetnr[j] = sheetnr[j] > tmp ? sheetnr[j] : tmp;
+    }
+  } 
+  fprintf(stdout, "Maximum indices to regard: %d %d %d\n", sheetnr[0], sheetnr[1], sheetnr[2]);
+  for (i=0;i<NDIM;i++) {
+    fprintf(stdout, "For axis %d: (%5.5lg\t%5.5lg\t%5.5lg) with %5.5lg\n", i, (Vector[i][0]*sheetnr[i]), (Vector[i][1]*sheetnr[i]), (Vector[i][2]*sheetnr[i]), Norm(Vector[i]));
+  }
+
+  //if (!strncmp(stage, "Cell", 4)) {
+    // parse in atoms for quicker processing
+    struct Atoms *atombuffer = malloc(sizeof(struct Atoms)*numbercell);
+    bufptr = CellBuffer; 
+    bufptr += GetNextline(bufptr, line)*sizeof(char);
+    bufptr += GetNextline(bufptr, line)*sizeof(char);
+    for (i=0;i<numbercell;i++) {
+      if ((length = GetNextline(bufptr, line)) != 0) {
+        bufptr += length*sizeof(char);
+        sscanf(line, "%s %lg %lg %lg", atombuffer[i].name, &(atombuffer[i].x[0]), &(atombuffer[i].x[1]), &(atombuffer[i].x[2]));
+        fprintf(stdout, "Read Atombuffer Nr %i: %s %5.5lg %5.5lg %5.5lg\n", i+1, atombuffer[i].name, atombuffer[i].x[0], atombuffer[i].x[1], atombuffer[i].x[2]);
+      } else {
+        fprintf(stdout, "Error reading Atom Nr. %i\n", i+1);
+        break;
+      }
+    }
+    SheetFile = fopen(SheetFilename, "w");
+    if (SheetFile == NULL) {
+      fprintf(stderr, "ERROR: main - can't open %s for writing\n", SheetFilename);
+      exit(255);
+    }
+    SheetFileAligned = fopen(SheetFilenameAligned, "w");
+    if (SheetFile == NULL) {
+      fprintf(stderr, "ERROR: main - can't open %s for writing\n", SheetFilenameAligned);
+      exit(255);
+    }
+    // Now create the sheet
+    double index[NDIM];
+    int nr;//, nummer = 0;
+    numbersheet = 0;
+    index[axis[2]] = 0;
+    // initialise pseudo random number generator with given seed
+    fprintf(stdout, "Initialising pseudo random number generator with given seed %d.\n", seed);
+    srand(seed);
+    //for (index[axis[0]] = 0; index[axis[0]] <= sheetnr[axis[0]]; index[axis[0]]++) { // NOTE: minor axis may start from 0! Check on this later ...
+    for (index[axis[0]] = -sheetnr[axis[0]]+1; index[axis[0]] < sheetnr[axis[0]]; index[axis[0]]++) { // NOTE: minor axis may start from 0! Check on this later ...
+      //for (index[axis[1]] = 0; index[axis[1]] <= sheetnr[axis[1]]; index[axis[1]]++) {  // These are all the cells that need be checked on
+      for (index[axis[1]] = -sheetnr[axis[1]]+1; index[axis[1]] < sheetnr[axis[1]]; index[axis[1]]++) {  // These are all the cells that need be checked on
+        // Calculate offset in cartesian coordinates
+        offset = MatrixTrafo(Vector, index);
+
+        //fprintf(stdout, "Now dealing with numbercell atoms in unit cell at R = (%lg,%lg,%lg)\n", offset[0], offset[1], offset[2]);
+        for (nr = 0; nr < numbercell; nr++) {
+          percentage = rand()/(RAND_MAX+1.0);
+          //fprintf(stdout, "Lucky number for %d is %lg >? %lg\n", nr, percentage, randomness[nr]); 
+          if (percentage >= randomness[nr]) {
+            // Create coordinates at atom site
+            coord = VectorAdd(atombuffer[nr].x, offset);
+            //fprintf(stdout, "Atom Nr. %i: ", (numbersheet+1));         
+            //PrintVector(stdout, coord);
+            // project down on major and minor Tubevectors and check for length if out of sheet
+            tempvector = MatrixTrafoInverse(coord, TubevectorInverse);
+            if (((tempvector[axis[0]] + MYEPSILON) >= 0) && ((factors[0] - tempvector[axis[0]]) > MYEPSILON) &&
+                ((tempvector[axis[1]] + MYEPSILON) >= 0) && ((factors[1] - tempvector[axis[1]]) > MYEPSILON) &&
+                ((tempvector[axis[2]] + MYEPSILON) >= 0) && ((factors[2] - tempvector[axis[2]]) > MYEPSILON)) { // check if within rotated cell          numbersheet++;
+              //if (nummer >= 2)  strcpy(atombuffer[nr].name, "O");
+              //nummer++;
+              fprintf(SheetFile, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", atombuffer[nr].name, coord[0], coord[1], coord[2]);
+              // rotate to align the sheet in xy plane
+              x1 = coord[0]*cos(-angle) + coord[1] * sin(-angle);
+              x2 = coord[0]*(-sin(-angle)) + coord[1] * cos(-angle);
+              x3 = coord[2];
+              fprintf(SheetFileAligned, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", atombuffer[nr].name, x1, x2, x3);
+              //fprintf(SheetFile, "O\t%5.5lg\t%5.5lg\t%5.5lg\n", coord[0], coord[1], coord[2]);
+              //fprintf(stdout, "%s/%d\t(%lg\t%lg\t%lg)\t", atombuffer[nr].name, numbersheet+1, coord[0], coord[1], coord[2]);
+              //PrintVector(stdout, tempvector);
+              numbersheet++;
+              //fprintf(stdout, "%i,", nr);
+            } //else {
+              //numbersheet++;
+              //fprintf(SheetFile, "B\t%lg\t%lg\t%lg\n", coord[0], coord[1], coord[2]);
+              //fprintf(stdout, "O \t(%lg\t%lg\t%lg)\n", coord[0], coord[1], coord[2]);
+              //fprintf(stdout, "!!%i!!, ", nr);
+            //}
+            Free(tempvector, "Main: At stage Sheet - tempvector");
+            Free(coord, "Main: At stage Sheet - coord");
+          }
+        }
+        Free(offset, "Main: At stage Sheet - offset");
+      }
+      //fprintf(stdout, "\n";
+    }
+
+    fclose(SheetFile);
+    fclose(SheetFileAligned);
+    AddAtomicNumber(SheetFilename,numbersheet, Vector, Recivector); // prepend atomic number and comment
+    AddAtomicNumber(SheetFilenameAligned,numbersheet, Vector, Recivector); // prepend atomic number and comment
+    AddSheetInfo(SheetFilename,axis,chiral, factors, seed, numbercell, randomness);
+    fprintf(stdout, "\nThere are %i atoms in the created sheet.\n", numbersheet);
+
+    strncpy(stage, "Sheet", 5);
+  //}
+  SheetBuffer = ReadBuffer(SheetFilename, &length);
+  
+  
+  // ======================== STAGE: Tube ==============================
+  // The tube starts with the rectangular (due to the orthogonalization) sheet
+  // just created (or read). Along the minor axis it is rolled up, i.e. projected
+  // from a 2d surface onto a cylindrical surface (x,y,z <-> r,alpha,z). The only
+  // thing that's a bit complex is that the sheet it not aligned along the cartesian
+  // axis but along major and minor. That's why we have to transform the atomic
+  // cartesian coordinates into the orthogonal tubevector base, do the rolling up
+  // there (and regard that minor and major axis must not necessarily be of equal
+  // length) and afterwards transform back again (where we need the $halfaxis due to
+  // the above possible inequality).
+
+  FILE *TubeFile = NULL;
+  FILE *TubeFileAligned = NULL;
+  
+  Debug ("STAGE: Tube\n");
+  if (!strncmp(stage, "Sheet", 4)) {
+    TubeFile = fopen(TubeFilename, "w");
+    if (TubeFile == NULL) {
+      fprintf(stderr, "ERROR: Main - can't open %s for writing\n", TubeFilename);
+      exit(255);
+    }
+    TubeFileAligned = fopen(TubeFilenameAligned, "w");
+    if (TubeFile == NULL) {
+      fprintf(stderr, "ERROR: Main - can't open %s for writing\n", TubeFilenameAligned);
+      exit(255);
+    }
+    bufptr = SheetBuffer;
+    bufptr += GetNextline(bufptr, line);  // write numbers to file
+    bufptr += GetNextline(bufptr, line);  // write comment to file
+
+    //cog = CenterOfGravity(bufptr, numbersheet);
+    //cog_projected = MatrixTrafoInverse(cog, TubevectorInverse);
+    //fprintf(stdout, "\nCenter of Gravity at (%5.5lg\t%5.5lg\t%5.5lg) and projected at (%5.5lg\t%5.5lg\t%5.5lg)\n", cog[0], cog[1], cog[2], cog_projected[0], cog_projected[1], cog_projected[2]);
+    
+    // restart
+    bufptr = SheetBuffer;
+    bufptr += GetNextline(bufptr, line);  // write numbers to file
+    bufptr += GetNextline(bufptr, line);  // write numbers to file
+
+    // determine half axis as tube vector not necessarily have the same length
+    double halfaxis[NDIM];
+    for (i=0;i<NDIM;i++)
+      halfaxis[i] = factors[0]*Norm(Tubevector[axis[0]])/Norm(Tubevector[i]);
+      
+    double arg, radius;  
+    for (i=0;i<numbersheet;i++) {
+      // scan next atom 
+      bufptr += GetNextline(bufptr, line);
+      sscanf(line, "%s %lg %lg %lg", name, &atom[0], &atom[1], &atom[2]);
+
+      //  transform atom coordinates in cartesian system to the axis eigensystem
+      x =  MatrixTrafoInverse(atom, TubevectorInverse);
+      //x = VectorAdd(y, cog_projected);
+      //free(y);
+      
+      // roll up (project (x,y,z) on cylindrical coordinates (radius,arg,z))
+      arg = 2.*M_PI*x[axis[0]]/(factors[0]) - M_PI;  // is angle
+      radius = 1./(2.*M_PI);  // is length of sheet in units of axis vector, divide by pi to get radius (from circumference)
+      // fprintf(stdout, "arg: %5.2f (c%2.2f,s%2.2f)\t",$arg, cos($arg), sin($arg));
+      x[axis[0]] = cos(arg)*halfaxis[axis[0]]*radius; //(radius+x[axis[2]]/halfaxis[axis[2]]); // as both vectors are not normalized additional betrag has to be taken into account!
+      x[axis[2]] = sin(arg)*halfaxis[axis[2]]*radius; //(radius+x[axis[2]]/halfaxis[axis[2]]); // due to the back-transformation from eigensystem to cartesian one
+      //fprintf(stdout, "rotated: (%5.2f,%5.2f,%5.2f)\n",x[0],x[1],x[2]);
+      atom_transformed = MatrixTrafo(Tubevector, x);
+      fprintf(TubeFile, "%s\t%lg\t%lg\t%lg\n", name, atom_transformed[0], atom_transformed[1], atom_transformed[2]);
+      // rotate and flip to align tube in z-direction
+      x1 = atom_transformed[0]*cos(angle) + atom_transformed[1] * sin(angle);
+      x2 = atom_transformed[0]*(-sin(angle)) + atom_transformed[1] * cos(angle);
+      x3 = atom_transformed[2];
+      fprintf(TubeFileAligned, "%s\t%lg\t%lg\t%lg\n", name, x1, x3, x2);
+      //fprintf(stdout, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", name, atom_transformed[0], atom_transformed[1] ,atom_transformed[2]);
+
+      Free(x, "Main: at stage Tube - x"); 
+      Free(atom_transformed, "Main: at stage Tube - atom_transformed");
+    }
+
+    
+    fclose(TubeFile);
+    fclose(TubeFileAligned);
+    //free(cog);
+    //free(cog_projected);
+    AddAtomicNumber(TubeFilename,numbersheet, Vector, Recivector); // prepend atomic number and comment
+    AddAtomicNumber(TubeFilenameAligned,numbersheet, Vector, Recivector); // prepend atomic number and comment
+    AddSheetInfo(TubeFilename,axis,chiral, factors, seed, numbercell, randomness);
+    fprintf(stdout, "\nThere are %i atoms in the created tube.\n", numbersheet);
+    
+    strncpy(stage, "Tube", 4);
+  } else {
+  }
+
+  TubeBuffer = ReadBuffer(TubeFilename, &length);
+
+  // ======================== STAGE: Torus =============================
+  // The procedure for the torus is very much alike to the one used to make the
+  // tube. Only the projection is not from 2d surface onto a cylindrical one but
+  // from a cylindrial onto a torus surface
+  // (x,y,z) <-> (cos(s)*(R+r*cos(t)), sin(s)*(R+rcos(t)), r*sin(t)).
+  // Here t is the angle within the tube with radius r, s is the torus angle with
+  // radius R. We get R from the tubelength (that's why we need lengthfactor to
+  // make it long enough). And due to fact that we have it already upon a cylindrical
+  // surface, r*cos(t) and r*sin(t) already reside in $minoraxis and $noaxis.
+
+  FILE *TorusFile;
+  
+  Debug ("STAGE: Torus\n");
+  if (!strncmp(stage, "Tube", 4)) {
+    TorusFile = fopen(TorusFilename, "w");
+    if (TorusFile == NULL) {
+      fprintf(stderr, "ERROR: main - can't open %s for writing\n", TorusFilename);
+      exit(255);
+    }
+    bufptr = TubeBuffer;
+    bufptr += GetNextline(bufptr, line);  // write numbers to file
+    bufptr += GetNextline(bufptr, line);  // write comment to file
+    
+    //cog = CenterOfGravity(bufptr, numbersheet);
+    //cog_projected = MatrixTrafoInverse(cog, TubevectorInverse);
+    //fprintf(stdout, "\nCenter of Gravity at (%5.5lg\t%5.5lg\t%5.5lg) and projected at (%5.5lg\t%5.5lg\t%5.5lg)\n", cog[0], cog[1], cog[2], cog_projected[0], cog_projected[1], cog_projected[2]);
+    
+    // determine half axis as tube vectors not necessarily have same length
+    double halfaxis[NDIM];
+    for (i=0;i<NDIM;i++)
+      halfaxis[i] = Norm(Tubevector[axis[1]])/Norm(Tubevector[i]);
+      
+    double arg, radius;  
+    for (i=0;i<numbersheet;i++) {
+      // scan next atom 
+      bufptr += GetNextline(bufptr, line);
+      sscanf(line, "%s %lg %lg %lg", name, &atom[0], &atom[1], &atom[2]);
+
+      //  transform atom coordinates in cartesian system to the axis eigensystem
+      x = MatrixTrafoInverse(atom, TubevectorInverse);
+      //x = VectorAdd(y, cog_projected);
+      //free(y);
+      
+      // roll up (project (x,y,z) on cylindrical coordinates (radius,arg,z))
+      arg = 2.*M_PI*x[axis[1]]/(factors[1]) - M_PI;  // is angle
+      radius = (factors[1])/(2.*M_PI) + x[axis[0]]/halfaxis[axis[0]];  // is length of sheet in units of axis vector, divide by pi to get radius (from circumference)
+      // fprintf(stdout, "arg: %5.2f (c%2.2f,s%2.2f)\t",$arg, cos($arg), sin($arg));
+      x[axis[0]] = cos(arg)*halfaxis[axis[0]]*radius; // as both vectors are not normalized additional betrag has to be taken into account!
+      x[axis[1]] = sin(arg)*halfaxis[axis[1]]*radius; // due to the back-transformation from eigensystem to cartesian one
+      //fprintf(stdout, "rotated: (%5.2f,%5.2f,%5.2f)\n",x[0],x[1],x[2]);
+      atom_transformed = MatrixTrafo(Tubevector, x);
+      fprintf(TorusFile, "%s\t%lg\t%lg\t%lg\n", name, atom_transformed[0], atom_transformed[1] ,atom_transformed[2]);
+      //fprintf(stdout, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", name, atom_transformed[0], atom_transformed[1] ,atom_transformed[2]);
+
+      Free(x, "Main: at stage Torus - x"); 
+      Free(atom_transformed, "Main: at stage Torus - atom_transformed");
+    }
+    
+    fclose(TorusFile);
+    //free(cog);
+    //free(cog_projected);
+    AddAtomicNumber(TorusFilename,numbersheet, Vector, Recivector); // prepend atomic number and comment
+    AddSheetInfo(TorusFilename,axis,chiral, factors, seed, numbercell, randomness);
+    fprintf(stdout, "\nThere are %i atoms in the created torus.\n", numbersheet);
+
+    strncpy(stage, "Torus", 5);
+  } else {
+  }
+
+  // Free memory
+  for (i=0; i<NDIM; i++ )  {
+    Free(Vector[i], "Main: end of stages - *Vector");
+    Free(Recivector[i], "Main: end of stages - *Recivector");
+    Free(Tubevector[i], "Main: end of stages - *Tubevector");
+    Free(TubevectorInverse[i], "Main: end of stages - *TubevectorInverse");
+  }
+  Free(atom, "Main: end of stages - atom");
+  Free(Vector, "Main: end of stages - Vector");
+  Free(Recivector, "Main: end of stages - Recivector");
+  Free(Tubevector, "Main: end of stages - Tubevector");
+  Free(TubevectorInverse, "Main: end of stages - TubevectorInverse");
+  Free(randomness, "Main: at stage Sheet - randomness");
+
+  if (CellBuffer != NULL) Free(CellBuffer, "Main: end of stages - CellBuffer");
+  if (SheetBuffer != NULL) Free(SheetBuffer, "Main: end of stages - SheetBuffer");
+  if (TubeBuffer != NULL) Free(TubeBuffer, "Main: end of stages - TubeBuffer");
+
+  Free(CellFilename, "Main: end of stafes - CellFilename");
+  Free(SheetFilename, "Main: end of stafes - CellFilename");
+  Free(TubeFilename, "Main: end of stafes - CellFilename");
+  Free(TorusFilename, "Main: end of stafes - CellFilename");
+  Free(SheetFilenameAligned, "Main: end of stafes - CellFilename");
+  Free(TubeFilenameAligned, "Main: end of stafes - CellFilename");
+  
+  // exit
+  exit(0);
+}
