#! @PERL@ -w
#
# This code helps in creating nanotube structures by specifying a unit cell and some
# multiplicative factors which are used to first generate a rectangular sheet, then
# rolling it up for the tube and again for the torus. All intermediate steps may be
# supplied by file and then only later geometries are created.
#
# The hows and whys are explained in the comments before each section.
# 
# \author Frederik Heber
# \date 15.06.07

use POSIX;

our $pi = 3.1415926535;			# no pi in perl
our $MYEPSILON = 1e-13;			# machine precision
	
our @vectorA1;
our @vectorA2;
our @vectorA3;
our @vector = ( \@vectorA1, \@vectorA2, \@vectorA3);
our @betrag;
our @RecivectorA1;
our @RecivectorA2;
our @RecivectorA3;
our @Recivector = ( \@RecivectorA1, \@RecivectorA2, \@RecivectorA3);
our @TubevectorA1;
our @TubevectorA2;
our @TubevectorA3;
our @Tubevector = ( \@TubevectorA1, \@TubevectorA2, \@TubevectorA3);
our	@Tubebetrag = (0,0,0);
our @TubevectorInverseA1;
our @TubevectorInverseA2;
our @TubevectorInverseA3;
our @TubevectorInverse = ( \@TubevectorInverseA1, \@TubevectorInverseA2, \@TubevectorInverseA3);
our	@TubebetragInverse = (0,0,0);
our @sheetnr;
our @buffer;
our @buffer2;
our @buffer3;

# ++++++++ Functions +++++++++++++++++++++++++++

# Convertes units using the unit programme
# \param $unit unit
# \param $SrcUnit source unit name
# \param $DstUnit destination unit name
# \return converted unit
sub UnitConversion
{
	local ($unit, $SrcUnit, $DstUnit) = ($_[0], $_[1], $_[2]);	# make local variables from given parameters
	$return = `units \'$unit, $SrcUnit,\' \'$DstUnit\' | head -n 1 | awk -F"* " {'print $2'}`;	
	return $return; # return value
}

# Adds commentary stuff (needed for further levels) to xyz files
# \param $filename  file name
# \param $atomicnumber number of atoms in cell
sub AddAtomicNumber
{
	local ($filename, $atomicnumber) = ($_[0], $_[1]);
	our @vector;
	our @Recivector;
	our @sheetnr;
	my $date = `date`;
	chomp($date);
	# fill buffer
	open(file1,$filename);
	@buffer = <file1>;
	close(file1);

	# open for writing and prepend
	open(file2,">$filename");
	print file2 "$atomicnumber\n";	# atomic number
	print file2 "\tgenerated with Nanotube creator on $date\n"; # ... and comment
	print file2 @buffer;	# append buffer

	# Add primitive vectors as comment
	print file2 "\n****************************************\n\n";
	print file2 "Primitive vectors\n";
	print file2 "a(1) = $vector[0][0]\t$vector[0][1]\t $vector[0][2]\n";
	print file2 "a(2) = $vector[1][0]\t$vector[1][1]\t $vector[1][2]\n";
	print file2 "a(3) = $vector[2][0]\t$vector[2][1]\t $vector[2][2]\n";
	print file2 "\nVolume = ".$volume."\n";
	print file2 "\nReciprocal Vectors\n";
	print file2 "b(1) = $Recivector[0][0]\t$Recivector[0][1]\t$Recivector[0][2]\t\n";
	print file2 "b(2) = $Recivector[1][0]\t$Recivector[1][1]\t$Recivector[1][2]\t\n";
	print file2 "b(3) = $Recivector[2][0]\t$Recivector[2][1]\t$Recivector[2][2]\t\n";

	close(file2);	# close file
}

# Transposes a 3x3 matrix
# \param **matrixref
sub Transpose
{
	local ($matrixref) = @_;
	for (my $i=0;$i<3;$i++) {
		for (my $j=0;$j<$i;$j++) {
			my $tmp = ${${$matrixref}[$i]}[$j];
			${${$matrixref}[$i]}[$j] = ${${$matrixref}[$j]}[$i];
			${${$matrixref}[$j]}[$i] = $tmp;
		}
	}
}

# Adds some more commentary stuff to xyz files
sub AddSheetInfo
{
	local ($filename, $majoraxis, $minoraxis, $noaxis, $n, $m, $radius, $length) = @_;
	# open for writing and prepend
	open(file2,">>$filename");
	
	# Add primitive vectors as comment
	print file2 "\n****************************************\n\n";
	print file2 "Axis $majoraxis\t$minoraxis\t$noaxis\n";
	print file2 "(n,m) $n $m\n";
	print file2 "factors $radius $length\n";
	close(file2);
}

# Reads either value from stdin or recognizes if old value shall be used again by simple return key.
# \param $_[0] old value
# \return either old value or newly entered one
sub GetValue
{
	$input = <INPUT>;
	chomp $input;
	if ($input) {
		return $input;
	} else {
		print $_[0]."\n";
		return $_[0];
	}
}

# approximatively finds the greatest common divisor of two real numbers.
# \param @array 3x3 matrix with row unit cell vectors
# \param $axis1 major axis number
# \param $axis2 minor axis number
# \return approximative greatest common divisor
# taken from http://www.perlmonks.org/?node_id=56906
sub gcf {
  my ($x, $y) = @_;
  ($x, $y) = ($y, $x % $y) while $y;
  return $x;
}

# Evaluates scalar product
# \param vector array 1
# \param vector array 2
# \return skp
sub skp
{
	my ($vec1ref, $vec2ref) = (@_);
	my $result = 0;
	for (my $i=0;$i<3;$i++) {
		#print ${$vec1ref}[$i]." * ".${$vec2ref}[$i]."\n";
		$result += ${$vec1ref}[$i] * ${$vec2ref}[$i];
	}
	return $result;
}

# Calculate norm with help of skp()
sub norm
{
	return sqrt(skp(@_));
}

# Evaluates projection of one vector onto another
# \param $vec1ref reference to 3 vector upon which is projected
# \param $vec2ref reference to 3 vector which is projected
#	\return projection factor
sub projection
{
	my ($vec1ref, $vec2ref) = (@_);
	return skp(@_)/skp($vec1ref, $vec1ref);
}

# Matrix transformation
# \param $matrixref reference to 3x3 matrix A
# \param $vectorref reference to 3 vector b
# \return reference to resulting 3 vector Ab = x
sub MatrixTrafo
{
	my ($matrixref, $vectorref) = @_;
	@return = (0,0,0);
	for (my $i=0;$i<3;$i++) {
		for (my $j=0;$j<3;$j++) {
			#print ${${$matrixref}[$i]}[$j]." * ".${$vectorref}[$i]."\n";
			$return[$j] += ${${$matrixref}[$i]}[$j] * ${$vectorref}[$i];
		}
	}
	return @return;
}

# Fixed GramSchmidt-Orthogonalization for 3 vectors
# \param @orthvector reference to 3x3 matrix
# \param @orthbetrag reference to 3 vector with vector magnitudes
# \param @axis major-, minor- and noaxis for specific order for the GramSchmidt procedure
sub Orthogonalize
{
	local ($orthvector, $orthbetrag, $axis) = @_;
	local $betrag1;
	local $betrag2;
	
	# orthogonalize axis[0] vector
	#print "Orthvectors: ${${$orthvector}[0]}[0]\t${${$orthvector}[0]}[1]\t${${$orthvector}[0]}[2]\n";
	#print "Orthbetrag: ${$orthbetrag}[0]\n";
	#print "OrthAxis: ${$axis}[0]\t${$axis}[1]\t${$axis}[2]\n";
	$betrag1 = projection(${$orthvector}[${$axis}[1]],${$orthvector}[${$axis}[0]]);
	for (my $k=0;$k<3;$k++) {
		${${$orthvector}[$axis[0]]}[$k] -= ${${$orthvector}[$axis[1]]}[$k]*$betrag1;
	}
	# orthogonalize axis[2] vector
	$betrag1 = projection(${$orthvector}[${$axis}[0]],${$orthvector}[${$axis}[2]]);
	$betrag2 = projection(${$orthvector}[${$axis}[1]],${$orthvector}[${$axis}[2]]);
	for (my $k=0;$k<3;$k++) {
		${${$orthvector}[${$axis}[2]]}[$k] -= ${${$orthvector}[${$axis}[0]]}[$k]*$betrag1 + ${${$orthvector}[${$axis}[1]]}[$k]*$betrag2;
	}
}

# calculate scalar product of vectors and store.
# \param $vector reference to 3x3 matrix with row vectors
# \param $betrag reference to 3 vector for found magnitudes
sub DetermineVectorLengths
{
	local ($vector, $betrag) = @_;
	for (my $i=0; $i<3; $i++) {
		${$betrag}[$i] = skp(${$vector}[$i],${$vector}[$i]);
		print "\t$i: ${$betrag}[$i]";
	}
}

# Add two vectors
# \param $vec1ref reference to first vector
# \param $vec2ref reference to second vector
# \return reference to vector sum
sub VectorAdd
{
	my ($vec1ref, $vec2ref) = @_;
	@return = (0,0,0);
	for (my $k=0;$k<3;$k++) {
		$return[$k] = ${$vec1ref}[$k] + ${$vec2ref}[$k];
	}
	return @return;
}

# Determines the inverse of a 3x3 matrix
# \param $matrix1ref matrix to be inverted
# \param $matrix2ref inverted matrix on exit
# formula taken from Wikipedia (Regulaere Matrix)
sub MatrixInversion
{
	local ($matrix1ref,$matrix2ref) = @_;
	my $det = Det($matrix1ref);
	#print "Determinant is $det\n";
	${${$matrix2ref}[0]}[0] = (${${$matrix1ref}[1]}[1]*${${$matrix1ref}[2]}[2] - ${${$matrix1ref}[1]}[2]*${${$matrix1ref}[2]}[1])/$det;
	${${$matrix2ref}[1]}[0] = (${${$matrix1ref}[0]}[2]*${${$matrix1ref}[2]}[1] - ${${$matrix1ref}[0]}[1]*${${$matrix1ref}[2]}[2])/$det;
	${${$matrix2ref}[2]}[0] = (${${$matrix1ref}[0]}[1]*${${$matrix1ref}[1]}[2] - ${${$matrix1ref}[0]}[2]*${${$matrix1ref}[1]}[1])/$det;
	${${$matrix2ref}[0]}[1] = (${${$matrix1ref}[1]}[2]*${${$matrix1ref}[2]}[0] - ${${$matrix1ref}[1]}[0]*${${$matrix1ref}[2]}[2])/$det;
	${${$matrix2ref}[1]}[1] = (${${$matrix1ref}[0]}[0]*${${$matrix1ref}[2]}[2] - ${${$matrix1ref}[0]}[2]*${${$matrix1ref}[2]}[0])/$det;
	${${$matrix2ref}[2]}[1] = (${${$matrix1ref}[0]}[2]*${${$matrix1ref}[1]}[0] - ${${$matrix1ref}[0]}[0]*${${$matrix1ref}[1]}[2])/$det;
	${${$matrix2ref}[0]}[2] = (${${$matrix1ref}[1]}[0]*${${$matrix1ref}[2]}[1] - ${${$matrix1ref}[1]}[1]*${${$matrix1ref}[2]}[0])/$det;
	${${$matrix2ref}[1]}[2] = (${${$matrix1ref}[0]}[1]*${${$matrix1ref}[2]}[0] - ${${$matrix1ref}[0]}[0]*${${$matrix1ref}[2]}[1])/$det;
	${${$matrix2ref}[2]}[2] = (${${$matrix1ref}[0]}[0]*${${$matrix1ref}[1]}[1] - ${${$matrix1ref}[0]}[1]*${${$matrix1ref}[1]}[0])/$det;

	print "Checking inverse ... ";
	my $flag = 0;
	for (my $i=0;$i<3;$i++) {
		for (my $j=0;$j<3;$j++) {
			my $tmp = 0;
			for (my $k=0;$k<3;$k++) {
				$tmp += ${${$matrix1ref}[$i]}[$k]*${${$matrix2ref}[$j]}[$k];
			}
			if ($flag == 0) {
				if ($i == $j) {
					$flag = fabs(1 - $tmp) > $MYEPSILON ? 1 : 0;
				} else {
					$flag = fabs($tmp) > $MYEPSILON ? 1 : 0;
				}
			}
			#print "$tmp";
			#if ($j != 2) { print ","; }
			#else { print "\n"; }
		}
		#if ($i != 2) { print ")\t("; }
		#else { print ")\n"; }
	}
	if ($flag == 0) {	print "ok\n"; } else { print "false\n"; }
}

# Evaluates determinant of 3x3 matrix
# \param $matrixref matrix whose determinant is to be calculated
sub Det
{
	local ($matrixref) = @_;
	my $return = 0;
	$return =   ${${$matrixref}[0]}[0] * (${${$matrixref}[1]}[1]*${${$matrixref}[2]}[2] - ${${$matrixref}[1]}[2]*${${$matrixref}[2]}[1])
						- ${${$matrixref}[1]}[1] * (${${$matrixref}[0]}[0]*${${$matrixref}[2]}[2] - ${${$matrixref}[0]}[2]*${${$matrixref}[2]}[0])
						+ ${${$matrixref}[2]}[2] * (${${$matrixref}[0]}[0]*${${$matrixref}[1]}[1] - ${${$matrixref}[0]}[1]*${${$matrixref}[1]}[0]);
	return $return;
}

# Parse for a pattern and return line split up by spaces
sub ParseFor
{
	my ($pattern, $arrayref) = @_;
	my $line;
	foreach $line (@{$arrayref}) {
		if ($line =~ /$pattern/)	{
			return (split(/\s+/,$line));
		} 
	}
}

# Reads contents of a file into an line-by-line array
# \param filename
# \return reference to array
sub ReadFile
{
	open(FILE,$_[0]);
	@storage = <FILE>;
	close(FILE);
	return @storage;
}

# A sheet, specified by axis[0] and axis[1] with vectors in $matrixref, is diametered
# \param $matrixref reference to 3x3 matrix with row vectors
# \param $axis reference to 3 vector with permutation of axis indices [0,1,2]
# \param $lengthfactor factor for axis[0]
# \param $radiusfactor factor for axis[1]
#	\return biggest diameter of sheet 
sub DetermineBiggestDiameter
{
	local ($matrixref, $axis, $lengthfactor, $radiusfactor) = @_;
	# get greatest diameter of new sheet
	my @diameter = (0,0);
	for(my $i=0;$i<3;$i++) {
		$diameter[0] += (${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}-${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor})
									 *(${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}-${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor});
		$diameter[1] += (${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}+${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor})
									 *(${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}+${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor});
	}
	if (($diameter[0] - $diameter[1]) > $MYEPSILON) {
		$biggest = 0;
	} else {
		$biggest = 1;
	}	
	@diameter = (sqrt($diameter[0]),sqrt($diameter[1]));
	printf("\n\nMajor diameter of the sheet is %5.5f, minor diameter is %5.5f.\n",$diameter[$biggest],$diameter[($biggest+1)%2]);
	
	return $diameter[$biggest];
}	


# ++++++++ Beginning +++++++++++++++++++++++++++

print "Nanotube Creator.\n=================\n\n";	# Welcome msg

my $file;

# Read command line arguments
if (!$ARGV[1]) {
	print "Nanotubes <file> <status>\n\t<file> is the basis used for the various xyz files that are created.\n\t<status> is either None, Cell, Sheet or Tube,\n\t\t whether the given file does not yet exist, specifies a unit cell, a sheet or a tube.\n";
	exit 255;
} else {
	$file = $ARGV[0];
	$stage = $ARGV[1];
	$file =~ s/(.*)\..*/$1/;
 	print "I will use \'$file' as base for the filenames.\n\n";
	
}

# for text input
open(INPUT,'-');

# Get primitive vectors, either ...
if ($stage !~ /[Nn]on/) {
	# ... from end of cell file
	print "Opening ".$file.".Cell.xyz to read primitive vectors in lines beginning with a(i) =  ...\n";
	open(CELL, $file.".Cell.xyz");
	@buffer = <CELL>;
	local $flag = 3;	# Marks whether vectors were found or not
	for($i=0; $i<@buffer; $i++) {
		#print "SEEK: $buffer[$i]";
		if ($buffer[$i] =~ /a\(1\) /) {			
			($name, $equal, @vectorA1) = split(/\s+/,$buffer[$i]);
			print "$name was found: ($vector[0][0], $vector[0][1], $vector[0][2])\n";
			$flag--;
		} elsif ($buffer[$i] =~ /a\(2\) /) {			
			($name, $equal, @vectorA2) = split(/\s+/,$buffer[$i]);
			print "$name was found: ($vector[1][0], $vector[1][1], $vector[1][2])\n";
			$flag--;
		} elsif ($buffer[$i] =~ /a\(3\) /) {			
			($name, $equal, @vectorA3) = split(/\s+/,$buffer[$i]);
			print "$name was found: ($vector[2][0], $vector[2][1], $vector[2][2])\n";
			$flag--;
		}
	}
	# parse number of atems in unit cell
	$numbercell = $buffer[0]; # number of atoms is first number in first line
	chomp $numbercell;
	print "\nThere are $numbercell atoms in the unit cell.\n";

	if ($flag != 0) { $stage = "None";	} # Set to None if vectors could not be found 
}

if ($stage =~ /[Nn]on/) {
	# ... or from INPUT
	# give some explanation to help imagination
	print "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";
	
	print "Enter 1st primitive vector (x z y): ";
	$line = <INPUT>;
	chomp($line);
	@vectorA1 = split(/\s+/,$line);
	print "Enter 2nd primitive vector (x z y): ";
	$line = <INPUT>;
	chomp($line);
	@vectorA2 = split(/\s+/,$line);
	print "Enter 3rd primitive vector (x z y): ";
	$line = <INPUT>;
	chomp($line);
	@vectorA3 = split(/\s+/,$line);
}	

print "Cell vectors: (@{$vector[0]})\t(@{$vector[1]})\t(@{$vector[2]})\n";

# Calculate some remaining stuff which comes at the end of the xyz file
# volume is determinant of primitive vectors seen as a matrix (column vectors)
$volume = Det(\@vector);
print "Volume is $volume\n";
MatrixInversion(\@vector, \@Recivector);
print "Reciprocal cell vectors: (@{$Recivector[0]})\t(@{$Recivector[1]})\t(@{$Recivector[2]})\n";
$Recivolume = Det(\@Recivector);


# ============ CELL ===========================
# The cell is simply created by transforming relative coordinates within the cell
# into cartesian ones using the unit cell vectors.

if ($stage =~ /[Nn]on/) {
	# enter number of atoms in unit cell
	print "\nHow many atoms are in the unit cell: ";
	$numbercell = <INPUT>;
	chomp($numbercell);
	
	# open unit cell file
	open(XYZ,">$file".".Cell.xyz");	# for xyz output
	
	print "\nNext, you have to enter each atom in the cell as follows, e.g.\n";
	print "Enter \'ChemicalSymbol X Y Z\' (relative to primitive vectors): C 0.5 0.25 0.5\n\n";
	# Enter element and coordinates of each
	for ($nr = 0; $nr < $numbercell; $nr++) {
		do {
			print "Enter for ".($nr+1)." atom \'ChemicalSymbol X Y Z\' (relative to primitive vectors): ";
			($name, @atom) = split(/\s+/,<INPUT>);
		} while (@atom < 3);
		my @x = MatrixTrafo(\@vector, \@atom);
		print XYZ $name."\t".($x[0])."\t".($x[1])."\t".($x[2])."\n";
	}
	
	close(XYZ);
	AddAtomicNumber($file.".Cell.xyz", $numbercell); # prepend atomic number and comment
	print "\nThere are $numbercell atoms in the created cell.\n"; 
	
	# Read cell file into buffer for later geometries
	@buffer = ReadFile("<$file".".Cell.xyz");

	$stage = "Cell"
} 

# evaluate length of each unit cell vector
print "\nVektorbetraege: ";
DetermineVectorLengths(\@vector,\@betrag);

# ============ 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.

if ($stage =~ /[Cc]ell/) {
	# get radial axis of tube
	print "\nSpecify the axis permutation that is going to be perpendicular to the sheet [tubeaxis, torusaxis, noaxis]: ";
	$line = <INPUT>;
	chomp($line);
	if (!$line) {
		$majoraxis = 0;
		$minoraxis = 1;
		$noaxis = 2;
	} else {
		($majoraxis, $minoraxis, $noaxis) = split(/\s+/,$line);
	}

	# find GCD for the two "longest" vector magnitudes
	#if ($betrag[($noaxis+1)%3] < $betrag[($noaxis+2)%3]) {
	#	$majoraxis = ($noaxis+2)%3;
	#	$minoraxis = ($noaxis+1)%3;
	#} else {
	#	$majoraxis = ($noaxis+1)%3;
	#	$minoraxis = ($noaxis+2)%3;
	#}
	#$gcf = ceil(1000*$betrag[$majoraxis])/gcf(ceil(1000*$betrag[$majoraxis]), ceil(1000*$betrag[$minoraxis]));	# factor to a desired precision and round as they may not always match
	#print "Major axis is $majoraxis, minor axis is $minoraxis and greatest common factor is $gcf: ".$betrag[$majoraxis]." == ".($betrag[$minoraxis]*$gcf)."\n";
	
	my $flag = 0;
	do {	# the discrete angle
		print "\nNow specify the two natural numbers (n m) defining the chiral angle: ";
		($n,$m) = split(/\s+/,<INPUT>);
		chomp($m);
	
		# Compute rotated cell vectors	
		for ($i=0; $i<3; $i++) {
			$Tubevector[$majoraxis][$i] =  $n*$vector[$majoraxis][$i] + $m*$vector[$minoraxis][$i];
			$Tubevector[$minoraxis][$i] = -$m*$vector[$majoraxis][$i] + $n*$vector[$minoraxis][$i];
			$Tubevector[$noaxis][$i] = $vector[$noaxis][$i];
		}
		print "\nTubevektorbetraege: ";
		DetermineVectorLengths(\@Tubevector,\@Tubebetrag);

		# calculate angle and radius		
		$tubecircum = sqrt($Tubebetrag[$minoraxis]);
		$tubelength = sqrt($Tubebetrag[$majoraxis]);
		$skp = skp($vector[$majoraxis],$Tubevector[$majoraxis]);
		$chiralangle = acos($skp/sqrt($betrag[$majoraxis]*$Tubebetrag[$majoraxis]))/$pi*180;

		print "\nGive integer factors for radius and length of tube: ";
		$line = <INPUT>;
		chomp($line);
		($radiusfactor, $lengthfactor) = split(/\s+/, $line);
	
		printf("\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", $chiralangle, $tubecircum*$radiusfactor/$pi, $tubelength*$lengthfactor, $tubelength*$lengthfactor/$pi);
		
		print "Satisfied? ";
		$flag = <INPUT>;
	} while ($flag !~ /[yY]/);
} else {
	# Read sheetnumbers from file
	@buffer2 = ReadFile($file.".Sheet.xyz");
	($pattern, $majoraxis, $minoraxis, $noaxis) = ParseFor("Axis", \@buffer2);
	($pattern, $n, $m) = ParseFor("(n,m)", \@buffer2);
	($pattern, $radiusfactor, $lengthfactor) = ParseFor("factors", \@buffer2);
		
	# Compute rotated cell vectors	
	for ($i=0; $i<3; $i++) {
		$Tubevector[$majoraxis][$i] =  $n*$vector[$majoraxis][$i] + $m*$vector[$minoraxis][$i];
		$Tubevector[$minoraxis][$i] = -$m*$vector[$majoraxis][$i] + $n*$vector[$minoraxis][$i];
		$Tubevector[$noaxis][$i] = $vector[$noaxis][$i];
	}
	print "\nTubevektorbetraege: ";
	DetermineVectorLengths(\@Tubevector,\@Tubebetrag);

	# calculate angle and radius		
	$tubecircum = sqrt($Tubebetrag[$minoraxis]);
	$tubelength = sqrt($Tubebetrag[$majoraxis]);
	$skp = skp($vector[$majoraxis],$Tubevector[$majoraxis]);
	$chiralangle = acos($skp/sqrt($betrag[$majoraxis]*$Tubebetrag[$majoraxis]))/$pi*180;
	printf("\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", $chiralangle, $tubecircum*$radiusfactor/$pi, $tubelength*$lengthfactor, $tubelength*$lengthfactor/$pi);
}

# Now for the rotated sheet we need orthogonal major and minor axis, but with major unchanged
@axis = ($majoraxis,$minoraxis,$noaxis);
Orthogonalize(\@Tubevector,\@Tubebetrag,\@axis);
# print new Tubebetraege
print "\nTubevektorbetraege after Orthogonalization: ";
DetermineVectorLengths(\@Tubevector,\@Tubebetrag);

print "Tube vectors: (@{$Tubevector[0]})\t(@{$Tubevector[1]})\t(@{$Tubevector[2]})\n";
MatrixInversion(\@Tubevector, \@TubevectorInverse);
Transpose(\@TubevectorInverse);
Transpose(\@Recivector);
print "Inverted tube vectors: (@{$TubevectorInverse[0]})\t(@{$TubevectorInverse[1]})\t(@{$TubevectorInverse[2]})\n";

# check touched original(!) cells in each direction (in a circle with biggest sheet diameter as radius)
$biggestdiameter = DetermineBiggestDiameter(\@Tubevector, \@axis, \$lengthfactor, \$radiusfactor);
for ($i=0;$i<3;$i++) {
	$sheetnr[$i] = 0;
}
for ($i=0;$i<3;$i++) {
	for ($j=0;$j<3;$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;
	}
}
print "Maximum indices to regard: @sheetnr\n";
for ($i=0;$i<3;$i++) {
	print "For axis $i: (".($vector[$i][0]*$sheetnr[$i])."\t".($vector[$i][1]*$sheetnr[$i])."\t".($vector[$i][2]*$sheetnr[$i]).") with ".sqrt($betrag[$i])."\n";
}

if ($stage =~ /[Cc]ell/) {
	# open sheet file
	open(XYZ,">$file".".Sheet.xyz");	# for xyz output
		
	# Quicker: scan all atoms beforehand
	for (my $nr = 0; $nr < $numbercell; $nr++) {
		push (@names, "none");		# These two are necessary as with ($name, @atom) = split(..); the same variable would be used and thus overwritten!
		push (@atoms, [0,0,0]);		# They sort of allocate the necessary memory
		($names[$nr],$atoms[$nr][0],$atoms[$nr][1],$atoms[$nr][2]) = split(/\s+/, $buffer[2+$nr]);	# take line by line
		print "Prereading $names[$nr]_$nr\t@{$atoms[$nr]}\n";
	}
	
	# Now create the sheet	
	$numbersheet = 0;
	@index = (0,0,0);
	@projection = (0,0,0);
	for ($index[$majoraxis] = -$sheetnr[$majoraxis]; $index[$majoraxis] < $sheetnr[$majoraxis]; $index[$majoraxis]++) {	# NOTE: minor axis may start from 0! Check on this later ...
		for ($index[$minoraxis] = 0; $index[$minoraxis] < $sheetnr[$minoraxis]; $index[$minoraxis]++) {	# These are all the cells that need be checked on
		#for ($index[$minoraxis] = -$sheetnr[$minoraxis]; $index[$minoraxis] <= $sheetnr[$minoraxis]; $index[$minoraxis]++) {	# These are all the cells that need be checked on
			# Calculate offset
			@offset = MatrixTrafo(\@vector, \@index);
			
			# print "Now dealing with $numbercell atoms in unit cell at R = (@offset): ";
			my @coord;
			for (my $nr = 0; $nr < $numbercell; $nr++) {
				# Create coordinates at atom site
				@coord = VectorAdd($atoms[$nr], \@offset);
				
				# project down on major and minor Tubevectors and check for length if out of sheet
				#@projection = MatrixTrafo(\@Recivector, \@coord);
				#print "length: $lengthfactor\tradius: $radiusfactor\tdepth: 1\taxis: $majoraxis $minoraxis $noaxis\tEps: $MYEPSILON ";
				#if ($numbersheet < 2) {
				#if ((($projection[$majoraxis] + $MYEPSILON) >= 0) && ((1 - $projection[$majoraxis]) > $MYEPSILON) &&
				#		(($projection[$minoraxis] + $MYEPSILON) >= 0) && ((1 - $projection[$minoraxis]) > $MYEPSILON) &&
			  #		(($projection[$noaxis] + $MYEPSILON) >= 0) && ((1 - $projection[$noaxis]) > $MYEPSILON)) { # check if within rotated cell
				#	$char = 'C';
				#} else {
				#	$char = 'O';
				#}
				@projection = MatrixTrafo(\@TubevectorInverse, \@coord);
				#print "projection: $projection[$majoraxis] $projection[$minoraxis] $projection[$noaxis]\n";
				#print "Atom Nr. ".($numbersheet+1).": @projection\n";
				if ((($projection[$majoraxis] + $MYEPSILON) >= 0) && (($lengthfactor - $projection[$majoraxis]) > $MYEPSILON) &&
						(($projection[$minoraxis] + $MYEPSILON) >= 0) && (($radiusfactor - $projection[$minoraxis]) > $MYEPSILON) &&
						(($projection[$noaxis] + $MYEPSILON) >= 0) && ((1 - $projection[$noaxis]) > $MYEPSILON)) { # check if within rotated cell
					$numbersheet++;
					print XYZ "$names[$nr]\t".($coord[0])."\t".($coord[1])."\t".($coord[2])."\n";
					#print XYZ "$char\t".($coord[0])."\t".($coord[1])."\t".($coord[2])."\n";
					#print "$nr,";
				}	#else {
					#$numbersheet++;
					#print XYZ "$char\t".($coord[0])."\t".($coord[1])."\t".($coord[2])."\n";
					#print "#$nr#, ";
				#}
			}
			#print "\n";
		}
	}
	
	close(XYZ);
	AddAtomicNumber($file.".Sheet.xyz",$numbersheet); # prepend atomic number and comment
	AddSheetInfo($file.".Sheet.xyz",$majoraxis,$minoraxis,$noaxis,$n,$m,$radiusfactor,$lengthfactor);
	print "\nThere are $numbersheet atoms in the created sheet.\n"; 

	$stage = "Sheet";
}

open (file3, $file.".Sheet.xyz");
@buffer2 = <file3>;
close(file3);
$numbersheet = $buffer2[0];
chomp($numbersheet);

# ============ 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).

if ($stage =~ /[Ss]heet/) {
	# open tube file
	open(XYZ,">$file".".Tube.xyz");	# for xyz output
	
	# Prepend begin
	print XYZ "$buffer2[0]";	# write number of atoms ...
	print XYZ "$buffer2[1]";	# ... and comment
	
	# determine center of gravity
#	@cog = (0,0,0);
#	for ($nr=0;$nr<$buffer2[0];$nr++) {		# for each atom
#		($name, @atom) = split(/\s+/, $buffer2[2+$nr]);
#		if (@atom >= 3) {
#			for (my $i=0;$i<3;$i++) {
#				$cog[$i] += $atom[$i];
#			}
#		}
#	}
#	for (my $i=0;$i<3;$i++) {
#		$cog[$i] /= -$buffer2[0];
#	}
#	@cog_projected = MatrixTrafo(\@TubevectorInverse, \@cog);

#	print "\nCenter of Gravity at (@cog) and projected at (@cog_projected)\n";

	# determine half axis as tube vector not necessarily half same length
	for (my $i=0;$i<3;$i++) {
		$halfaxis[$i] = $radiusfactor*($tubecircum)/sqrt($Tubebetrag[$i]);
	}
	# Rolling up the sheet
	my @x = (0,0,0);
	print "\nRadiusfactor: $radiusfactor\t pi: $pi\taxis ($majoraxis,$minoraxis,$noaxis)\n";
	for ($nr=0;$nr<$buffer2[0];$nr++) {		# for each atom
		($name, @atom) = split(/\s+/, $buffer2[2+$nr]);
		if (@atom >= 3) {
			#printf("Atom $nr: (%5.2f,%5.2f,%5.2f)\t",$atom[0],$atom[1],$atom[2]);
			# transform atom coordinates in cartesian system to the axis eigensystem
			@x = MatrixTrafo(\@TubevectorInverse, \@atom);
			#@x = VectorAdd(\@y, \@cog_projected);
			#printf("projected: (%5.2f,%5.2f,%5.2f)\t",$x[0],$x[1],$x[2]);

			# roll up (project (x,y,z) on cylindrical coordinates (radius,arg,z))
			$arg = 2*$pi*$x[$minoraxis]/($radiusfactor) - $pi;	# is angle
			$radius = 1/(2*$pi);	# is length of sheet in units of axis vector, divide by pi to get radius (from circumference)
			#printf("arg: %5.2f (c%2.2f,s%2.2f)\t",$arg, cos($arg), sin($arg));
			$x[$minoraxis] = cos($arg)*$halfaxis[$minoraxis]*$radius;	# due to the back-transformation from eigensystem to cartesian one
			$x[$noaxis] = sin($arg)*$halfaxis[$noaxis]*$radius;	# as both vectors are not normalized additional betrag has to be taken into account!
			#printf("rotated: (%5.2f,%5.2f,%5.2f)\n",$x[0],$x[1],$x[2]);
			@atom = MatrixTrafo(\@Tubevector, \@x);
			print XYZ $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
			#print $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
		} else {
			print "Error reading atom Nr. $nr\n";
		}
	}
	# append rest
	for ($i=$nr+2;$i<@buffer2;$i++) {
		print XYZ "$buffer2[$i]";
	}
	print XYZ "tubecircum $tubecircum\n";
	close(XYZ);
	print "\nThere are $numbersheet atoms in the created tube.\n"; 
	
	$stage = "Tube";
} 

# Read tubecircum from file
open (file3, "<$file".".Tube.xyz");
@buffer3 = <file3>;
close(file3);
($pattern, $tubecircum) = ParseFor("tubecircum", \@buffer3);	


# ============ 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.

if ($stage =~ /[Tt]ube/) {

	# open torus file
	open(XYZ,">$file".".Torus.xyz");	# for xyz output
	
	# Prepend begin
	print XYZ "$buffer3[0]";	# write number of atoms ...
	print XYZ "$buffer3[1]";	# ... and comment

	# determine center of gravity
#	@cog = (0,0,0);
#	for ($nr=0;$nr<$buffer3[0];$nr++) {		# for each atom
#		($name, @atom) = split(/\s+/, $buffer3[2+$nr]);
#		if (@atom >= 3) {
#			for (my $i=0;$i<3;$i++) {
#				$cog[$i] += $atom[$i];
#			}
#		}
#	}
#	for (my $i=0;$i<3;$i++) {
#		$cog[$i] /= -$buffer3[0];
#	}
#	@ocog_projected = MatrixTrafo(\@TubevectorInverse, \@cog);
	
#	print "\nCenter of Gravity at (@cog) and projected at (@cog_projected)\n";

	# Determine	half axis
	for (my $i=0;$i<3;$i++) {
		$halfaxis[$i] = ($tubelength)/sqrt($Tubebetrag[$i]);
	}
	# Rolling up the sheet
	my @x = (0,0,0);
	#print "\nRadiusfactor: $radiusfactor\t pi: $pi\taxis ($majoraxis,$minoraxis,$noaxis)\n";
	for ($nr=0;$nr<$buffer3[0];$nr++) {		# for each atom
		($name, @atom) = split(/\s+/, $buffer3[2+$nr]);
		if (@atom >= 3) {
			#printf("Atom $nr: (%5.2f,%5.2f,%5.2f)\t",$atom[0],$atom[1],$atom[2]);

			# transform atom coordinates in cartesian system to the axis eigensystem and shift into center of gravity
			@x = MatrixTrafo(\@TubevectorInverse, \@atom);
			#@x = VectorAdd(\@y, \@cog_projected);
			#printf("projected: (%5.2f,%5.2f,%5.2f)\t",$x[0],$x[1],$x[2]);

			# roll up (project (x,z,y) on cylindrical coordinates (radius,arg,z))
			$arg = 2*$pi*$x[$majoraxis]/($lengthfactor)-$pi;	# is angle
			$radius = ($lengthfactor/(2*$pi) + $x[$minoraxis]/$halfaxis[$minoraxis]);	# is R + r*cos(t)
			#printf("arg: %5.2f (c%2.2f,s%2.2f)\t",$arg, cos($arg), sin($arg));
			$x[$majoraxis] = cos($arg)*$halfaxis[$majoraxis]*($radius);	# as both vectors are not normalized additional betrag has to be taken into account!
			$x[$minoraxis] = sin($arg)*$halfaxis[$minoraxis]*($radius);	# due to the back-transformation from eigensystem to cartesian one
			#printf("rotated: (%5.2f,%5.2f,%5.2f)\n",$x[0],$x[1],$x[2]);
			@atom = MatrixTrafo(\@Tubevector, \@x);

			print XYZ $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
			#print $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
		} else {
			print "Error reading atom Nr. $nr\n";
		}
	}

	# append rest
	for ($i=$nr+2;$i<@buffer3;$i++) {
		print XYZ "$buffer3[$i]";
	}
	print XYZ "tubecircum $tubecircum\n";
	close(XYZ);

	$stage = "Torus";
}

print "\nFinished with geometry $stage.\n";

# close stdin and exit
close(INPUT);
exit 0;

# ++++++++ This is the end +++++++++++++++++++++