source: util/src/NanoCreator.c@ e1a46d

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.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 betrag;
  int i;

  // first vector is untouched
  // second vector
  betrag = Projection(orthvector[axis[1]], orthvector[axis[0]]);
  fprintf(stdout,"%lg\t",betrag);
  for (i=0;i<NDIM;i++)
     orthvector[axis[1]][i] -= orthvector[axis[0]][i] * betrag;
  // third vector
  betrag = Projection(orthvector[axis[0]], orthvector[axis[2]]);
  fprintf(stdout,"%lg\t",betrag);
  for (i=0;i<NDIM;i++)
     orthvector[axis[2]][i] -= orthvector[axis[0]][i] * betrag;
  betrag = Projection(orthvector[axis[1]], orthvector[axis[2]]);
  fprintf(stdout,"%lg\n",betrag);
  for (i=0;i<NDIM;i++)
     orthvector[axis[2]][i] -= orthvector[axis[1]][i] * betrag;
}

/** 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;
}

/** Creates orthogonal vectors to directions axis[0] and axis[1], assuming that axis[2] is always orthogonal.
 * \param **OrthoVector contains vectors set on return with axis[2] equal to Vector[axis[2]]
 * \param **Vector vectors to orthogonalize against
 * \param *axis lookup for which direction is which.
 */
void CreateOrthogonalAxisVectors(double **OrthoVector, double **Vector, int *axis) {
  int i,j;
  double factor;
  // allocate memory
  int *TempAxis = (int *) Malloc(sizeof(int)*NDIM, "Main: *TempAxis");
  double **TempVectors = (double **) Malloc(sizeof(double *)*NDIM, "Main: *TempVectors");
  for (i=0; i<NDIM; i++)
    TempVectors[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: TempVectors");

  for (i=0; i<NDIM; i++)
    for (j=0; j<NDIM; j++)
      TempVectors[i][j] = Vector[i][j];
  // GramSchmidt generates Vector[1] orthogonal to Vector[0] and Vector[2] ortho. to Vector[1] and Vector[0]!
  TempAxis[0] = axis[2];  // (axis 2 is the orthogonal plane axis)
  TempAxis[1] = axis[0];
  TempAxis[2] = axis[1];
  Orthogonalize(TempVectors, TempAxis);
  factor = Norm(Vector[axis[0]])/Norm(TempVectors[TempAxis[2]]);
  factor *= (Projection(TempVectors[TempAxis[2]], Vector[axis[0]]) > 0) ? 1. : -1.;
  for (i=0; i<NDIM; i++)
    OrthoVector[axis[0]][i] = TempVectors[TempAxis[2]][i]*factor;

  TempAxis[1] = axis[1];
  TempAxis[2] = axis[0];
  for (i=0; i<NDIM; i++)
    for (j=0; j<NDIM; j++)
      TempVectors[i][j] = Vector[i][j];
  Orthogonalize(TempVectors, TempAxis);
  factor = Norm(Vector[axis[1]])/Norm(TempVectors[TempAxis[2]]);
  factor *= (Projection(TempVectors[TempAxis[2]], Vector[axis[1]]) > 0) ? 1. : -1.;
  for (i=0; i<NDIM; i++)
    OrthoVector[axis[1]][i] = TempVectors[TempAxis[2]][i]*factor;

  for (i=0; i<NDIM; i++)
    OrthoVector[axis[2]][i] = Vector[axis[2]][i];
  //print vectors
  fprintf(stdout, "Orthogonal vectors are: \n");
  for (i=0; i<NDIM; i++) {
    for (j=0; j<NDIM; j++)
      fprintf(stdout, "%lg\t", OrthoVector[axis[i]][j]);
    fprintf(stdout, "\n");
  }

  // free memory
  Free(TempAxis, "CreateOrthogonalAxisVectors: *TempAxis");
  for (i=0; i<NDIM; i++ )
    Free(TempVectors[i], "CreateOrthogonalAxisVectors: TempVectors");
  Free(TempVectors, "CreateOrthogonalAxisVectors: *TempVectors");
};

// ================================= other functions ==============================

/** Prints a debug message.
 * \param *msg debug message
 */
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, **OrthoVector, **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");
  OrthoVector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *OrthoVector");
  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");
    OrthoVector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: OrthoVector");
    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]);

    // create orthogonal vectors individually for each unit cell vector
    CreateOrthogonalAxisVectors(OrthoVector, Vector, axis);
    fprintf(stdout, "Orthogonal vector axis[0] %lg\n", Projection(Vector[axis[0]], OrthoVector[axis[0]]));
    fprintf(stdout, "Orthogonal vector axis[1] %lg\n", Projection(Vector[axis[1]], OrthoVector[axis[1]]));

    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] = (double)chiral[0] * OrthoVector[axis[0]][i] - (double)chiral[1] * OrthoVector[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];
      }
      // here we assume, that Vector[axis[2]] is along z direction!
      gsl_matrix *M = gsl_matrix_alloc(2,2);
      gsl_matrix *C = gsl_matrix_alloc(2,2);
      gsl_matrix *evec = gsl_matrix_alloc(2,2);
      gsl_vector *eval = gsl_vector_alloc(2);
      gsl_vector *v = gsl_vector_alloc(2);
      gsl_vector *u = gsl_vector_alloc(2);
      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(2);
      gsl_matrix_set(C, 0,0, Vector[axis[0]][0]);
      gsl_matrix_set(C, 1,0, Vector[axis[0]][1]);
      gsl_matrix_set(C, 0,1, Vector[axis[1]][0]);
      gsl_matrix_set(C, 1,1, Vector[axis[1]][1]);
      gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0, C, C, 0.0, M);
      fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
      gsl_eigen_symmv(M, eval, evec, w);
      gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_ABS_DESC);
      fprintf(stdout, "Eigenvalues: %lg\t%lg\n", gsl_vector_get(eval,0), gsl_vector_get(eval,1));
      fprintf(stdout, "Eigenvectors: \t%lg\t%lg\n\t\t%lg\t%lg\n", gsl_matrix_get(evec,0,0), gsl_matrix_get(evec,0,1), gsl_matrix_get(evec,1,0), gsl_matrix_get(evec,1,1));
      gsl_matrix_set(M, 0,0, 0.);
      gsl_matrix_set(M, 1,0, 1.);
      gsl_matrix_set(M, 0,1, -gsl_vector_get(eval,1)/gsl_vector_get(eval,0));
      gsl_matrix_set(M, 1,1, 0.);
      gsl_vector_set(v,0,(double)chiral[0]);
      gsl_vector_set(v,1,(double)chiral[1]);
      gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0, evec, M, 0.0, C);
      gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0, C, evec, 0.0, M);
      fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
      gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
      fprintf(stdout, "Looking for factor to integer...\n");
      for(i=1;i<(chiral[0]+chiral[1])*(chiral[0]+chiral[1]);i++) {
        x1 = gsl_vector_get(u,0)*(double)i;
        x2 = gsl_vector_get(u,1)*(double)i;
        x3 =
        fprintf(stdout, "%d: %d\t%d vs. %lg\t%lg\n",i, ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), (x1), (x2));
        if (( fabs( ((int)(x1+x1/fabs(x1)*.5)) - (x1) ) < 1e-6) && ( fabs( ((int)(x2+x2/fabs(x2)*.5)) - (x2) ) < 1e-6 )) {
          gsl_blas_dscal((double)i, u);
          break;
        }
      }
      fprintf(stdout, "(c,d) = (%lg,%lg)\n",gsl_vector_get(u,0), gsl_vector_get(u,1));

      // get length
      double x[NDIM];
      for (i=0;i<NDIM;i++)
        x[i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
      angle = Norm(x)/Norm(Tubevector[axis[1]]) ;//ScalarProduct(x,Tubevector[axis[1]])/Norm(Tubevector[axis[1]]);
      fprintf(stdout, "angle is %lg\n", angle);
      for (i=0;i<NDIM;i++) {
        Tubevector[axis[1]][i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
      }

      // Probe
      gsl_matrix_set(M, 0,0, Vector[axis[0]][0]);
      gsl_matrix_set(M, 1,0, Vector[axis[0]][1]);
      gsl_matrix_set(M, 0,1, Vector[axis[1]][0]);
      gsl_matrix_set(M, 1,1, Vector[axis[1]][1]);
      gsl_vector_set(v,0,(double)chiral[0]);
      gsl_vector_set(v,1,(double)chiral[1]);
      gsl_blas_dgemv(CblasNoTrans, 1.0, M, u, 0.0, eval);
      gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
      x1=1.;
      gsl_blas_ddot(u,eval,&x1);
      fprintf(stdout, "Testing (c,d): (a,b) M^t M (c,d)^t = 0 ? : %lg\n", x1);

      gsl_matrix_free(M);
      gsl_matrix_free(C);
      gsl_matrix_free(evec);
      gsl_vector_free(eval);
      gsl_vector_free(v);
      gsl_vector_free(u);
      gsl_eigen_symmv_free(w);

      if (fabs(x1) > 1e-6) {
        fprintf(stderr,"Resulting TubeVectors of axis %d and %d and not orthogonal, aborting.\n", axis[0], axis[1]);
        return(128);
      }


      angle = Projection(Tubevector[axis[1]], Vector[axis[0]]);
      fprintf(stdout, "Projection Tubevector1 axis[0] %lg %lg\n", angle, 1./angle);
      angle = Projection(Tubevector[axis[1]], Vector[axis[1]]);
      fprintf(stdout, "Projection Tubevector1 axis[1] %lg %lg\n", angle, 1./angle);

/*      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 orthogonal vectors individually for each unit cell vector
    CreateOrthogonalAxisVectors(OrthoVector, Vector, axis);
    fprintf(stdout, "Orthogonal vector axis[0] %lg", Projection(Vector[axis[0]], OrthoVector[axis[0]]));
    fprintf(stdout, "Orthogonal vector axis[1] %lg", Projection(Vector[axis[1]], OrthoVector[axis[1]]));
    // 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] * OrthoVector[axis[0]][i] + chiral[1] * OrthoVector[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];
    }
    // here we assume, that Vector[axis[2]] is along z direction!
    gsl_matrix *M = gsl_matrix_alloc(2,2);
    gsl_matrix *C = gsl_matrix_alloc(2,2);
    gsl_matrix *evec = gsl_matrix_alloc(2,2);
    gsl_vector *eval = gsl_vector_alloc(2);
    gsl_vector *v = gsl_vector_alloc(2);
    gsl_vector *u = gsl_vector_alloc(2);
    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(2);
    gsl_matrix_set(C, 0,0, Vector[axis[0]][0]);
    gsl_matrix_set(C, 1,0, Vector[axis[0]][1]);
    gsl_matrix_set(C, 0,1, Vector[axis[1]][0]);
    gsl_matrix_set(C, 1,1, Vector[axis[1]][1]);
    gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0, C, C, 0.0, M);
    fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
    gsl_eigen_symmv(M, eval, evec, w);
    gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_ABS_DESC);
    fprintf(stdout, "Eigenvalues: %lg\t%lg\n", gsl_vector_get(eval,0), gsl_vector_get(eval,1));
    fprintf(stdout, "Eigenvectors: \t%lg\t%lg\n\t\t%lg\t%lg\n", gsl_matrix_get(evec,0,0), gsl_matrix_get(evec,0,1), gsl_matrix_get(evec,1,0), gsl_matrix_get(evec,1,1));
    gsl_matrix_set(M, 0,0, 0.);
    gsl_matrix_set(M, 1,0, 1.);
    gsl_matrix_set(M, 0,1, -gsl_vector_get(eval,1)/gsl_vector_get(eval,0));
    gsl_matrix_set(M, 1,1, 0.);
    gsl_vector_set(v,0,(double)chiral[0]);
    gsl_vector_set(v,1,(double)chiral[1]);
    gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0, evec, M, 0.0, C);
    gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0, C, evec, 0.0, M);
    fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
    gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
    fprintf(stdout, "Looking for factor to integer...\n");
    for(i=1;i<(chiral[0]+chiral[1])*(chiral[0]+chiral[1]);i++) {
      x1 = gsl_vector_get(u,0)*(double)i;
      x2 = gsl_vector_get(u,1)*(double)i;
      x3 =
      fprintf(stdout, "%d: %d\t%d vs. %lg\t%lg\n",i, ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), (x1), (x2));
      if (( fabs( ((int)(x1+x1/fabs(x1)*.5)) - (x1) ) < 1e-6) && ( fabs( ((int)(x2+x2/fabs(x2)*.5)) - (x2) ) < 1e-6 )) {
        gsl_blas_dscal((double)i, u);
        break;
      }
    }
    fprintf(stdout, "(c,d) = (%lg,%lg)\n",gsl_vector_get(u,0), gsl_vector_get(u,1));

    // get length
    double x[NDIM];
    for (i=0;i<NDIM;i++)
      x[i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
    angle = Norm(x)/Norm(Tubevector[axis[1]]) ;//ScalarProduct(x,Tubevector[axis[1]])/Norm(Tubevector[axis[1]]);
    fprintf(stdout, "angle is %lg\n", angle);
    for (i=0;i<NDIM;i++) {
      Tubevector[axis[1]][i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
    }

    // Probe
    gsl_matrix_set(M, 0,0, Vector[axis[0]][0]);
    gsl_matrix_set(M, 1,0, Vector[axis[0]][1]);
    gsl_matrix_set(M, 0,1, Vector[axis[1]][0]);
    gsl_matrix_set(M, 1,1, Vector[axis[1]][1]);
    gsl_vector_set(v,0,(double)chiral[0]);
    gsl_vector_set(v,1,(double)chiral[1]);
    gsl_blas_dgemv(CblasNoTrans, 1.0, M, u, 0.0, eval);
    gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
    x1=1.;
    gsl_blas_ddot(u,eval,&x1);
    fprintf(stdout, "Testing (c,d): (a,b) M^t M (c,d)^t = 0 ? : %lg\n", x1);

    gsl_matrix_free(M);
    gsl_matrix_free(C);
    gsl_matrix_free(evec);
    gsl_vector_free(eval);
    gsl_vector_free(v);
    gsl_vector_free(u);
    gsl_eigen_symmv_free(w);

    if (fabs(x1) > 1e-6) {
      fprintf(stderr,"Resulting TubeVectors of axis %d and %d and not orthogonal, aborting.\n", axis[0], axis[1]);
      return(128);
    }

    // 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]);
      }
    }
  }

  //int OrthOrder[NDIM] = { axis[2], axis[0], axis[1] };
  //Orthogonalize(Tubevector,OrthOrder);
  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 %lg\n", angle/M_PI*180., Projection(Vector[axis[0]], Vector[axis[1]]));
  if ( angle/M_PI*180. > 90 ) {
    fprintf(stderr, "There seems to be something wrong with the unit cell! for nanotube the angle should be 60 degrees for example!\n");
    return 1;
  }
  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 %lg\n", angle/M_PI*180., Projection(Tubevector[axis[0]], Tubevector[axis[1]]));
  //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.);
  if (fabs(Tubevector[axis[0]][0]) > MYEPSILON)
    angle = -M_PI/2. + acos(Tubevector[axis[0]][0]/Norm(Tubevector[axis[0]]));
  else
    angle = 0.;
  fprintf(stdout, "The absolute alignment rotation angle then is %lg %lg\n", angle/M_PI*180., Tubevector[axis[0]][0]/Norm(Tubevector[axis[0]]));
  fprintf(stdout, "\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. final 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, "Tubevectors 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+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+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, x3, x2, x1);  // order so that symmetry is along z axis
      //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);
}