package SpaceGroupObject;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
@EXPORT = qw();

use strict;

#========================================================
# 静的メンバー関数
#========================================================
sub ConvertToFraction
{
	my ($val, $s) = @_;
	my $Threshold = 0.01;
	$s = '' if(!defined $s);

	return "" if(abs($val) < $Threshold);
	return "+$s" if(abs($val - 1) < $Threshold);
	return "-$s" if(abs($val + 1) < $Threshold);
	for(my $i = 1 ; $i <= 24 ; $i++) {
		my $v = $val*$i;
		if(abs($v - int($v)) < $Threshold) {
			$val = "$v/$i";
			$val = "+$val" if($val !~ /^-/); #if($val > 0);
			last;
		}
	}
	return "$val$s";
}

sub BuildSymmetryOperationFromMatrix
{
	my ($x1, $y1, $z1, $t1) = @_;
#print "($x1,$y1,$z1,$t1) => ";
	$x1 = &ConvertToFraction($x1, 'x');
	$y1 = &ConvertToFraction($y1, 'y');
	$z1 = &ConvertToFraction($z1, 'z');
	$t1 = &ConvertToFraction(Utils::Reduce01($t1));
	my $s = "$x1$y1$z1$t1";
	$s =~ s/^\+//;
#print "($x1,$y1,$z1,$t1):$s\n";
	return $s;
}

sub SymOpToMatrix
{
	my ($symopstr) = (@_);

	$symopstr = lc $symopstr;
	$symopstr =~ s/'//g;
	my ($x, $y, $z, $t);

	$symopstr =~ s/x/\$xx/gi;
	$symopstr =~ s/y/\$yy/gi;
	$symopstr =~ s/z/\$zz/gi;
	my $xx = 0;
	my $yy = 0;
	my $zz = 0;
	$t = eval($symopstr);
	$xx = 1;
	$x = eval($symopstr) - $t;
	$xx = 0;
	$yy = 1;
	$y = eval($symopstr) - $t;
	$yy = 0;
	$zz = 1;
	$z = eval($symopstr) - $t;
	$t = Utils::Reduce01($t);

	return ($x, $y, $z, $t);
}

#========================================================
# 変数取得関数
#========================================================
sub SPGName            { return shift->{'SPGName'}; }
sub SetSPGName         { my ($this,$name)=@_; return $this->{'SPGName'} = $name; }
sub iSPG               { return shift->{'iSPG'}; }
sub SetiSPG            { my ($this,$i)=@_; return $this->{'iSPG'} = $i; }
sub iSet               { return shift->{'iSet'}; }
sub SetiSet            { my ($this,$i)=@_; return $this->{'iSet'} = $i; }
sub Volume             { return shift->{'ITableVolume'}; }
sub nTranslation {
	my $this = shift;
	return $this->{'nTranslation'} if($this->{'nTranslation'});
	$this->{"TranslationVector[1][x]"} = 0;
	$this->{"TranslationVector[1][y]"} = 0;
	$this->{"TranslationVector[1][z]"} = 0;
	return $this->{'nTranslation'} = 1;
}
sub nSymmetryOperation { return shift->{'nSymmetryOperation'}; }
sub SymmetryOperation  { my ($this, $i)=@_; return $this->{"SymmetryOperation[$i]"}; }
sub LatticeSystem      { return shift->{'LatticeSystem'}; }

#========================================================
# コンストラクタ、デストラクタ
#========================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#========================================================
# 一般関数
#========================================================
sub SetP1
{
	my ($this) = @_;
	$this->ClearSymmetryOperation();
	$this->SetSPGName("P 1");
	$this->SetiSPG(1);
	$this->SetiSet(1);
	$this->AnalyzeTranslation();
	$this->AddSymmetryOperation('x,y,z');
	$this->SetLatticeSystem('triclinic');
}

sub DoSymmetryOperation
{
	my ($this, $iop, $x, $y, $z, $IsReduce01) = @_;
	my $nSymmetryOperation  = $this->nSymmetryOperation();
	my $iSymmetryOperation  = $iop % $nSymmetryOperation;
	my $iTranslateOperation = int($iop / $nSymmetryOperation);
	my ($tx, $ty, $tz)      = $this->TranslationVector($iTranslateOperation+1);
	my $symop = $this->SymmetryOperation($iSymmetryOperation+1);
	$symop = lc $symop;
	$symop =~ s/([xyz])/\$$1/g;
	my ($sx,$sy,$sz) = ($symop =~ /^(.*?),(.*?),(.*)$/);
#print "is=$iop: T=($tx, $ty, $tz), S=($sx, $sy, $sz)\n";
	my $xs = eval($sx)+$tx;
	my $ys = eval($sy)+$ty;
	my $zs = eval($sz)+$tz;
	if($IsReduce01) {
		$xs = Utils::Reduce01(AtomSiteObject::Round01($xs));
		$ys = Utils::Reduce01(AtomSiteObject::Round01($ys));
		$zs = Utils::Reduce01(AtomSiteObject::Round01($zs));
	}
	return ($xs,$ys,$zs);
}

sub DoVelocitySymmetryOperation
{
	my ($this, $iop, $vx, $vy, $vz) = @_;
	return (0, 0, 0) if(not defined $vx);
	my $nSymmetryOperation = $this->nSymmetryOperation();
	my $iSymmetryOperation = $iop % $nSymmetryOperation;

	my $symop = $this->SymmetryOperation($iSymmetryOperation+1);
	my ($xo,$yo,$zo) = ($symop =~ /^(.*?),\s*(.*?),\s*(.*?)$/);
	my ($x1,$y1,$z1,$t1) = &SymOpToMatrix($xo);
	my ($x2,$y2,$z2,$t2) = &SymOpToMatrix($yo);
	my ($x3,$y3,$z3,$t3) = &SymOpToMatrix($zo);
	my $vxs = $x1*$vx + $y1*$vy + $z1*$vz;
	my $vys = $x2*$vx + $y2*$vy + $z2*$vz;
	my $vzs = $x3*$vx + $y3*$vy + $z3*$vz;
	return ($vxs,$vys,$vzs);
}

sub ClearSymmetryOperation
{
	my ($this) = @_;
	$this->{'nSymmetryOperation'} = 0;
#	$this->{"SymmetryOperation[$i]"} = $symop;
	return 0;
}

sub AddSymmetryOperation
{
	my ($this, $symop) = @_;
#print "symop: $symop\n";
	return unless($symop);
	my $n = $this->{'nSymmetryOperation'};
#print "n=$n: ";
	$n = 0 unless($this->{"SymmetryOperation[1]"});
#print "bef $n: $symop\n";
	for(my $i = 1 ; $i <= $n ; $i++) {
# suffixは1からはじまる
#print "i=$i: ", $this->{"SymmetryOperation[$i]"}, "\n";	
		return $n if(lc $symop eq lc $this->{"SymmetryOperation[$i]"});
	}
	$this->{'nSymmetryOperation'} = ++$n;
#print "nSy: ", $this->{'nSymmetryOperation'}, "\n";
	$this->{"SymmetryOperation[$n]"} = $symop;
#print "aft $n: $symop\n";
#	return $this->{'nSymmetryOperation'};
#print "n=$n\n";
	return $n;
}

sub SymmetryOperationByMatrix
{
	my ($this, $i) = @_;
	my $symop = $this->{"SymmetryOperation[$i]"};
	my ($xo,$yo,$zo) = ($symop =~ /^(.*?),\s*(.*?),\s*(.*?)$/);
	my ($x1,$y1,$z1,$t1) = &SymOpToMatrix($xo);
	my ($x2,$y2,$z2,$t2) = &SymOpToMatrix($yo);
	my ($x3,$y3,$z3,$t3) = &SymOpToMatrix($zo);
	return ($x1,$y1,$z1,$t1,
		$x2,$y2,$z2,$t2,
		$x3,$y3,$z3,$t3)
}

sub LatticeParameters {
	my ($this) = @_;
	return ($this->{'a'}, $this->{'b'}, $this->{'c'},
		$this->{'alpha'}, $this->{'beta'}, $this->{'gamma'});
}

sub SetLatticeParameters
{
	my ($this,$a,$b,$c,$alpha,$beta,$gamma) = @_;
	$a = Utils::Round($a, 6);
	$b = Utils::Round($b, 6);
	$c = Utils::Round($c, 6);
	$alpha = Utils::Round($alpha, 6);
	$beta  = Utils::Round($beta , 6);
	$gamma = Utils::Round($gamma, 6);

	$this->{'a'} = $a;
	$this->{'b'} = $b;
	$this->{'c'} = $c;
	$this->{'alpha'} = $alpha;
	$this->{'beta'}  = $beta;
	$this->{'gamma'} = $gamma;

	$this->SetLatticeSystem();

	return 1;
}

sub SetLatticeSystem
{
	my ($this, $latticesystem, $SearchByLatticeParameter, $tollatt, $tolangle) = @_;
	$tollatt  = 1.0e-5   if(!defined $tollatt);
	$tolangle = $tollatt if(!defined $tolangle);

	my $tolangle2 = $tolangle * 1.01;

	if(defined $latticesystem and $latticesystem ne '') {
		return $this->{LatticeSystem} = $latticesystem;
	}
	my $SPGName = $this->{SPGName};
	my $a = $this->{a};
	my $b = $this->{b};
	my $c = $this->{c};
	my $alpha = $this->{alpha};
	my $beta  = $this->{beta};
	my $gamma = $this->{gamma};
#print "SetLatticeSystem: $a,$b,$c,$alpha,$beta,$gamma\n";
#print "SPG=$SPGName\n";
	return if($a == 0.0 or $b == 0.0 or $c == 0.0);

	if((!$SearchByLatticeParameter or $SearchByLatticeParameter == 2) and $SPGName =~ /4\s*3/i) {
		return $this->{LatticeSystem} = 'cubic';
	}

	if((!$SearchByLatticeParameter or $SearchByLatticeParameter == 2) and $SPGName =~ /4/i) {
		return $this->{LatticeSystem} = 'tetragonal';
	}

	if((!$SearchByLatticeParameter or $SearchByLatticeParameter == 2) and $SPGName =~ /^P\s*6/i) {
		return $this->{LatticeSystem} = 'hexagonal';
	}

	if((!$SearchByLatticeParameter or $SearchByLatticeParameter == 2) and $SPGName =~ /^R/i 
	   and abs($a-$b) < 1.0e-4 and abs($a-$c) < 1.0e-4
	   and abs($alpha-$beta) < 1.0e-4 and abs($alpha-$gamma) < 1.0e-4) {
#print "R\n";
#		return $this->{LatticeSystem} = 'rhombohedral';
		return $this->{LatticeSystem} = 'trigonal';
	}

#print "$a $b $c $alpha $beta $gamma: tollat=$tollatt, tolangle=$tolangle2\n";
	if(abs($a-$b) < $tollatt) {
#print "check: $a - $c < $tollatt: ", abs($a-$c), "\n";
		if(abs($a-$c) < $tollatt) {
#print "Yes\n";
			if(abs($alpha-90.0)<$tolangle2 and abs($beta-90.0)<$tolangle2 and abs($gamma-90.0)<$tolangle2) {
				return $this->{LatticeSystem} = 'cubic';
			}
			elsif(abs($alpha-$beta)<$tolangle2 and abs($beta-$gamma)<$tolangle2) {
#print "TRI\n";
#exit;
				return $this->{LatticeSystem} = 'trigonal';
			}
		}
#print "check 120\n";
		if(abs($alpha-90.0)<$tolangle2 and abs($beta-90.0)<$tolangle2 and abs($gamma-90.0)<$tolangle2) {
#print "Tet\n";
			return $this->{LatticeSystem} = 'tetragonal';
		}
		elsif(abs($alpha-90.0)<$tolangle2 and abs($beta-90.0)<$tolangle2 and abs($gamma-120.0)<$tolangle2) {
#print "H\n";
			return $this->{LatticeSystem} = 'hexagonal';
		}
	}

	if(abs($alpha-90.0)<$tolangle2 and abs($beta-90.0)<$tolangle2 and abs($gamma-90.0)<$tolangle2) {
#print "O\n";
		return $this->{LatticeSystem} = 'orthorhombic';
	}
	elsif(  (abs($alpha-90.0)<$tolangle2 and abs($beta-90.0) <$tolangle2) or
		(abs($alpha-90.0)<$tolangle2 and abs($gamma-90.0)<$tolangle2) or
		(abs($beta-90.0) <$tolangle2 and abs($gamma-90.0)<$tolangle2) ) {
#print "MC\n";
			return $this->{LatticeSystem} = 'monoclinic';
	}

#print "TRIC\n";
	return $this->{LatticeSystem} = 'triclinic';
}

sub ExpandSymmetryOperation
{
	my ($this,$iCenter) = @_;

	my $SPG = $this->SPGName();
	my $nSymmetryOperation = $this->nSymmetryOperation();
#print "nSymmetryOperation:$nSymmetryOperation\n";
	my $nTranslation = $this->nTranslation();
	my @list;

	for(my $is = 0 ; $is < $nSymmetryOperation ; $is++) {
		my $symop = $this->SymmetryOperation($is+1);
		my ($x1,$y1,$z1,$t1,$x2,$y2,$z2,$t2,$x3,$y3,$z3,$t3)
				= $this->SymmetryOperationByMatrix($is+1);
		my $sx = &BuildSymmetryOperationFromMatrix($x1, $y1, $z1, $t1);
		my $sy = &BuildSymmetryOperationFromMatrix($x2, $y2, $z2, $t2);
		my $sz = &BuildSymmetryOperationFromMatrix($x3, $y3, $z3, $t3);
		my $str = "$sx,$sy,$sz";
		my $IsPassed = 1;
		for(my $i = 0 ; $i < @list ; $i++) {
			if($list[$i] eq $str) {
				$IsPassed = 0;
				last;
			}
		}
		push(@list, $str) if($IsPassed);
	}

	if($iCenter) {
		my $ns = @list;
		for(my $is = 0 ; $is < $ns ; $is++) {
			my $symop = $this->SymmetryOperation($is+1);
			my ($x1,$y1,$z1,$t1,
			$x2,$y2,$z2,$t2,
			$x3,$y3,$z3,$t3)
				= $this->SymmetryOperationByMatrix($is+1);

			my $sx = &BuildSymmetryOperationFromMatrix(
					-$x1, -$y1, -$z1, -$t1);
			my $sy = &BuildSymmetryOperationFromMatrix(
					-$x2, -$y2, -$z2, -$t2);
			my $sz = &BuildSymmetryOperationFromMatrix(
					-$x3, -$y3, -$z3, -$t3);
			my $str = "$sx,$sy,$sz";
			my $IsPassed = 1;
			for(my $i = 0 ; $i < @list ; $i++) {
				if($list[$i] eq $str) {
					$IsPassed = 0;
					last;
				}
			}
			push(@list, $str) if($IsPassed);
		}
	}

	my $ns = @list;
	for(my $is = 0 ; $is < $ns ; $is++) {
			my $symop = $list[$is]; #$this->SymmetryOperation($is+1);
			my ($xo,$yo,$zo) = ($symop =~ /^(.*?),\s*(.*?),\s*(.*?)$/);
			my ($x1,$y1,$z1,$t1) = &SymOpToMatrix($xo);
			my ($x2,$y2,$z2,$t2) = &SymOpToMatrix($yo);
			my ($x3,$y3,$z3,$t3) = &SymOpToMatrix($zo);

		for(my $it = 0 ; $it < $nTranslation ; $it++) {
			my ($tx,$ty,$tz) = $this->TranslationVector($it+1);
			my $sx = &BuildSymmetryOperationFromMatrix(
					$x1, $y1, $z1, $t1+$tx);
			my $sy = &BuildSymmetryOperationFromMatrix(
					$x2, $y2, $z2, $t2+$ty);
			my $sz = &BuildSymmetryOperationFromMatrix(
					$x3, $y3, $z3, $t3+$tz);
			my $str = "$sx,$sy,$sz";
			my $IsPassed = 1;
			for(my $i = 0 ; $i < @list ; $i++) {
				if($list[$i] eq $str) {
					$IsPassed = 0;
					last;
				}
			}
			push(@list, $str) if($IsPassed);
		}
	}
	$this->{'nSymmetryOperation'} = @list;
	for(my $i = 0 ; $i < @list ; $i++) {
		my $i1 = $i+1;
		$this->{"SymmetryOperation[$i1]"} = $list[$i];
	}

	return (@list);
}

sub AnalyzeTranslation
{
	my ($this) = @_;
	
	my $SPG = &SPGName($this);
	$this->{"TranslationVector[1][x]"} = 0;
	$this->{"TranslationVector[1][y]"} = 0;
	$this->{"TranslationVector[1][z]"} = 0;

	if($SPG =~ /^'?I/i) {
		$this->{"nTranslation"} = 2;
		$this->{"TranslationVector[2][x]"} = 0.5;
		$this->{"TranslationVector[2][y]"} = 0.5;
		$this->{"TranslationVector[2][z]"} = 0.5;
		return $this->{"nTranslation"};
	}

	if($SPG =~ /^'?A/i) {
		$this->{"nTranslation"} = 2;
		$this->{"TranslationVector[2][x]"} = 0;
		$this->{"TranslationVector[2][y]"} = 0.5;
		$this->{"TranslationVector[2][z]"} = 0.5;
		return $this->{"nTranslation"};
	}

	if($SPG =~ /^'?B/i) {
		$this->{"nTranslation"} = 2;
		$this->{"TranslationVector[2][x]"} = 0.5;
		$this->{"TranslationVector[2][y]"} = 0;
		$this->{"TranslationVector[2][z]"} = 0.5;
		return $this->{"nTranslation"};
	}

	if($SPG =~ /^'?C/i) {
		$this->{"nTranslation"} = 2;
		$this->{"TranslationVector[2][x]"} = 0.5;
		$this->{"TranslationVector[2][y]"} = 0.5;
		$this->{"TranslationVector[2][z]"} = 0;
		return $this->{"nTranslation"};
	}

	if($SPG =~ /^'?F/i) {
		$this->{"nTranslation"} = 4;
		$this->{"TranslationVector[2][x]"} = 0;
		$this->{"TranslationVector[2][y]"} = 0.5;
		$this->{"TranslationVector[2][z]"} = 0.5;
		$this->{"TranslationVector[3][x]"} = 0.5;
		$this->{"TranslationVector[3][y]"} = 0;
		$this->{"TranslationVector[3][z]"} = 0.5;
		$this->{"TranslationVector[4][x]"} = 0.5;
		$this->{"TranslationVector[4][y]"} = 0.5;
		$this->{"TranslationVector[4][z]"} = 0;
		return $this->{"nTranslation"};
	}

	$this->{"nTranslation"} = 1;
	return $this->{"nTranslation"};
}

sub GetTranslationArrays
{
	my @TranslationArrays = (
		{ 'Matrix' => [0.5, 0.5, 0.0], 'Name' => 'C' },
		{ 'Matrix' => [0.5, 0.0, 0.5], 'Name' => 'B' },
		{ 'Matrix' => [0.0, 0.5, 0.5], 'Name' => 'A' },
		{ 'Matrix' => [0.5, 0.5, 0.5], 'Name' => 'I' },
		);
	
	return @TranslationArrays;
}

sub GetTranslationalVectors
{
	my ($this) = @_;

	my $SPG        = $this->GetCSpaceGroup();
	my $SPGName    = $SPG->SPGName();

	my @VT = ([  0,   0,   0]);
	if($SPGName =~ /^\s*F/i) {
		@VT = (
			[  0,   0,   0],
			[  0, 1/2, 1/2],
			[1/2,   0, 1/2],
			[1/2, 1/2,   0],
			);
	}
	elsif($SPGName =~ /^\s*I/i) {
		@VT = (
			[  0,   0,   0],
			[1/2, 1/2, 1/2],
			);
	}
	elsif($SPGName =~ /^\s*A/i) {
		@VT = (
			[  0,   0,   0],
			[  0, 1/2, 1/2],
			);
	}
	elsif($SPGName =~ /^\s*B/i) {
		@VT = (
			[  0,   0,   0],
			[1/2,   0, 1/2],
			);
	}
	elsif($SPGName =~ /^\s*C/i) {
		@VT = (
			[  0,   0,   0],
			[1/2, 1/2,   0],
			);
	}
	return \@VT;
}

sub TranslationVector
{
	my ($this, $i) = @_;
	
	my $x = $this->{"TranslationVector[$i][x]"};
	my $y = $this->{"TranslationVector[$i][y]"};
	my $z = $this->{"TranslationVector[$i][z]"};
	return ($x, $y, $z);
}

sub GetMatrixConventionalToPrimitiveCell
{
	my ($this, $ReferToConventionalCell) = @_;
	my $SPGName = ($ReferToConventionalCell and $this->{ConventionalSPGName})? 
				$this->{ConventionalSPGName} : $this->SPGName();
#print "2********[$this->{ConventionalSPGName}][$SPGName]\n";

	if($SPGName =~ /^\s*F/i) {
		return [ [ 0.5, 0.5,   0], 
			 [ 0,   0.5, 0.5],
			 [ 0.5,   0, 0.5] ]
	}
	elsif($SPGName =~ /^\s*I/i) {
		return [ [-0.5, 0.5, 0.5], 
			 [ 0.5,-0.5, 0.5],
			 [ 0.5, 0.5,-0.5] ];
	}
	elsif($SPGName =~ /^\s*A/i) {
		return [ [ 1.0, 0.0, 0.0], 
			 [ 0.0, 0.5, 0.5],
			 [ 0.0,-0.5, 0.5] ];
	}
	elsif($SPGName =~ /^\s*B/i) {
		return [ [ 0.5, 0.0, 0.5], 
			 [ 0.0, 1.0, 0.0],
			 [-0.5, 0.0, 0.5] ];
	}
	elsif($SPGName =~ /^\s*C/i) {
		return [ [ 0.5, 0.5, 0.0], 
			 [-0.5, 0.5, 0.0],
			 [ 0.0, 0.0, 1.0] ];
	}
	return [ [ 1.0, 0.0, 0.0], 
		 [ 0.0, 1.0, 0.0],
		 [ 0.0, 0.0, 1.0] ];
}

sub GetMatrixConventionalToPrimitiveReciprocalCell
{
	my ($this) = @_;
	my $SPGName = $this->SPGName();

	if($SPGName =~ /^\s*F/i) {
		return [ [-1.0, 1.0, 1.0], 
			 [ 1.0,-1.0, 1.0],
			 [ 1.0, 1.0,-1.0] ];
	}
	elsif($SPGName =~ /^\s*I/i) {
		return [ [ 1.0, 1.0,   0], 
			 [ 0,   1.0, 1.0],
			 [ 1.0,   0, 1.0] ]
	}
	elsif($SPGName =~ /^\s*A/i) {
		return [ [ 1.0, 0.0, 0.0], 
			 [ 0.0, 1.0, 1.0],
			 [ 0.0,-1.0, 1.0] ];
	}
	elsif($SPGName =~ /^\s*B/i) {
		return [ [ 1.0, 0.0, 1.0], 
			 [ 0.0, 1.0, 0.0],
			 [-1.0, 0.0, 1.0] ];
	}
	elsif($SPGName =~ /^\s*C/i) {
		return [ [ 1.0, 1.0, 0.0], 
			 [-1.0, 1.0, 0.0],
			 [ 0.0, 0.0, 1.0] ];
	}
	return [ [ 1.0, 0.0, 0.0], 
		 [ 0.0, 1.0, 0.0],
		 [ 0.0, 0.0, 1.0] ];
}

sub IsFaceCenteredCell
{
	my ($this, $pAS, $eps) = @_;

	my @m = (
		[0.0, 0.5, 0.5],
		[0.5, 0.0, 0.5],
		[0.5, 0.5, 0.0],
		);

	return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}

sub IsBodyCenteredCell
{
	my ($this, $pAS, $eps) = @_;

	my @m = (
		[0.5, 0.5, 0.5]
		);

	return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}

sub IsACenteredCell
{
	my ($this, $pAS, $eps) = @_;

	my @m = (
		[0.0, 0.5, 0.5]
		);

	return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}

sub IsBCenteredCell
{
	my ($this, $pAS, $eps) = @_;

	my @m = (
		[0.5, 0.0, 0.5]
		);

	return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}

sub IsCCenteredCell
{
	my ($this, $pAS, $eps) = @_;

	my @m = (
		[0.5, 0.5, 0.0]
		);

	return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}

sub GuessBravaisLattice
{
	my ($this, $tollatt, $tolangle, $EPSCoord) = @_;

	my ($a, $b, $c, $alpha, $beta, $gamma) = $this->LatticeParameters();
	my $SPG     = $this->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	my $ls      = lc $this->LatticeSystem();

	my $LatticeSystem = ($ls eq '')? $SPG->SetLatticeSystem(undef, 1, $tollatt, $tolangle) : $ls;

	my $BravaisLattice;
	if($SPGName =~ /^\s*F/i) {
		$BravaisLattice = 'F';
	}
	elsif($SPGName =~ /^\s*I/i) {
		$BravaisLattice = 'I';
	}
	elsif($SPGName =~ /^\s*A/i) {
		$BravaisLattice = 'A';
	}
	elsif($SPGName =~ /^\s*B/i) {
		$BravaisLattice = 'B';
	}
	elsif($SPGName =~ /^\s*C/i) {
		$BravaisLattice = 'C';
	}
#print "SPG=$SPGName  BL=$BravaisLattice\n";

	if($BravaisLattice eq '') {
		my @as = $this->GetCExpandedAtomSiteListByOutputMode();
		if($this->IsFaceCenteredCell(\@as, $EPSCoord)) {
			$BravaisLattice = 'F';
		}
		elsif($this->IsBodyCenteredCell(\@as, $EPSCoord)) {
			$BravaisLattice = 'I';
		}
		elsif($this->IsACenteredCell(\@as, $EPSCoord)) {
			$BravaisLattice = 'A';
		}
		elsif($this->IsBCenteredCell(\@as, $EPSCoord)) {
			$BravaisLattice = 'B';
		}
		elsif($this->IsCCenteredCell(\@as, $EPSCoord)) {
			$BravaisLattice = 'C';
		}
		else {
			$BravaisLattice = 'P';
		}
	}

	return ($LatticeSystem, $BravaisLattice);
}

sub LatticeSystemFromSPGName
{
	my ($this, $SPGName) = @_;
	return 'cubic'      if($SPGName =~ /4/ and $SPGName =~ /3/);
	return 'cubic'      if($SPGName =~ /[FI]/ and $SPGName =~ /3/);
	return 'tetragonal' if($SPGName =~ /4/);
	return 'hexagonal'  if($SPGName =~ /6/);

	return 'trigonal'   if($SPGName =~ /3/);
	return 'orthogonal' if($SPGName =~ /2.*2/);
	return 'orthogonal' if($SPGName =~ /[mnabc].*[mnabc]/);
	return 'monoclinic' if($SPGName =~ /2/);
	return 'monoclinic' if($SPGName =~ /^\s*\S\s*[mnabc]/);
	return 'triclinic'  if($SPGName =~ /1/);
}

sub GetBravaisLattice
{
	my ($this) = @_;
	my $SPGName = $this->SPGName();
	my $LaueGroup = $this->LatticeSystem();
	my ($BravaisLattice) = ($SPGName =~ /^(\S)+/);
	if($BravaisLattice =~ /P/) {
		$BravaisLattice = "P($LaueGroup)";
	}
	elsif($BravaisLattice =~ /F/) {
		$BravaisLattice = "FC($LaueGroup)";
	}
	elsif($BravaisLattice =~ /I/) {
		$BravaisLattice = "BC($LaueGroup)";
	}
	elsif($BravaisLattice =~ /A/) {
		$BravaisLattice = "A($LaueGroup)";
	}
	elsif($BravaisLattice =~ /B/) {
		$BravaisLattice = "B($LaueGroup)";
	}
	elsif($BravaisLattice =~ /C/) {
		$BravaisLattice = "C($LaueGroup)";
	}
	$BravaisLattice = 'R(Rhombohedral)' if($BravaisLattice eq 'R');
#	$BravaisLattice = 'R(Rhombohedral)' if($BravaisLattice eq 'P(Trigonal)' or $BravaisLattice eq 'R');
	
	return $BravaisLattice;
}

sub IsCompatibleLatticeSystem
{
	my ($this, $CrystalLS, $SPGLS) = @_;

	if($CrystalLS =~ /cub/i) {
		return 1;
	}
	if($CrystalLS =~ /tet/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 1;
	}
	if($CrystalLS =~ /ort/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 1;
	}
	if($CrystalLS =~ /hex/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 1;
	}
	if($CrystalLS =~ /trig/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 1;
	}
	if($CrystalLS =~ /rho/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 1;
	}
	if($CrystalLS =~ /mon/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 0 if($SPGLS =~ /hex/i);
		return 0 if($SPGLS =~ /trig/i);
		return 0 if($SPGLS =~ /rho/i);
		return 1;
	}
	return 1 if($SPGLS =~ /trig/i or $SPGLS eq '');
	return 0;
}

sub GetLatticeConversionMatrix
{
	my ($this,  $PresetRule, $DirectionSelect, $ParametersSelect) = @_;

	my $sT  = new Math::MatrixReal(3, 3);
	my $T   = new Math::MatrixReal(3, 3);
	my $tRT = new Math::MatrixReal(3, 3);
	my $RT  = new Math::MatrixReal(3, 3);
	my $tR  = new Math::MatrixReal(3, 3);

	if($PresetRule eq 'No' or $PresetRule eq 'P1') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 1,  0,  0], 
				  [ 0,  1,  0],
				  [ 0,  0,  1] ] );
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'RhombHex') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 1, -1,  0], 
				  [ 0,  1, -1],
				  [ 1,  1,  1] ] );
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'HexRhomb') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 2/3, 1/3, 1/3],
				  [-1/3, 1/3, 1/3],
				  [-1/3,-2/3, 1/3] ]);
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'HexOrtho') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 1, 0, 0], 
				  [ 1, 2, 0],
				  [ 0, 0, 1] ] );
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'FCCPrim' or $PresetRule eq 'PrimFCC') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 0.5, 0.5,   0], 
				  [ 0,   0.5, 0.5],
				  [ 0.5,   0, 0.5] ] );
		if($PresetRule eq 'PrimFCC') {
			$sT = $sT->inverse();
		}
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'BCCPrim' or $PresetRule eq 'PrimBCC') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [-0.5, 0.5, 0.5], 
				  [ 0.5,-0.5, 0.5],
				  [ 0.5, 0.5,-0.5] ] );
		if($PresetRule eq 'PrimBCC') {
			$sT = $sT->inverse();
		}
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'ACenterPrim') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 1.0, 0.0, 0.0], 
				  [ 0.0, 0.5, 0.5],
				  [ 0.0,-0.5, 0.5] ] );
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'BCenterPrim') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 0.5, 0.0, 0.5], 
				  [ 0.0, 1.0, 0.0],
				  [-0.5, 0.0, 0.5] ] );
		$sT->transpose($sT);
	}
	elsif($PresetRule eq 'CCenterPrim') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [ 0.5, 0.5, 0.0], 
				  [-0.5, 0.5, 0.0],
				  [ 0.0, 0.0, 1.0] ] );
		$sT->transpose($sT);
	}
	else {
die "Error: Invalid conversion rule [$PresetRule].\n";

		print "  Rule is specified by Matrix:\n";
		$sT->assign(1, 1, eval(&GetArg('t11')));
		$sT->assign(1, 2, eval(&GetArg('t12')));
		$sT->assign(1, 3, eval(&GetArg('t13')));
		$sT->assign(2, 1, eval(&GetArg('t21')));
		$sT->assign(2, 2, eval(&GetArg('t22')));
		$sT->assign(2, 3, eval(&GetArg('t23')));
		$sT->assign(3, 1, eval(&GetArg('t31')));
		$sT->assign(3, 2, eval(&GetArg('t32')));
		$sT->assign(3, 3, eval(&GetArg('t33')));
	}

print "Conversion rule:\n";
	my $IsChosen = 0;
	if($ParametersSelect eq 'ai') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "Real space vectors: (a'i) = (tij)(ai)\n";
			$T = $sT;
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "Real space vectors: (ai) = (tij)(a'i)\n";
			$T = $sT->inverse();
			$IsChosen = 1;
		}
	}
	elsif($ParametersSelect eq 'xyz') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "Real space coordinates: (x'i) = (tij)(xi)\n";
			$T->transpose($sT);
			$T = $T->inverse();
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "Real space vectors: (xi) = (tij)(x'i)\n";
			$T->transpose($sT);
			$IsChosen = 1;
		}
	}
	elsif($ParametersSelect eq 'a*i') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "Reciprocal space vectors: (a'*i) = (tij)(a*i)\n";
			$T->transpose($sT);
			$T = $T->inverse();
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "Reciprocal space vectors: (a*i) = (tij)(a'*i)\n";
			$T->transpose($sT);
			$IsChosen = 1;
		}
	}
	elsif($ParametersSelect eq 'hkl') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "Reciprocal space coordinates: (h'i) = (tij)(hi)\n";
			$T = $sT;
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "Reciprocal space coordinates: (hi) = (tij)(h'i)\n";
			$T = $sT->inverse();
			$IsChosen = 1;
		}
	}
	unless($IsChosen) {
print "Selected option {Parameters: $ParametersSelect} in {$DirectionSelect} direction is not supported.\n";
		return 1;
	}

#実空間座標(x,y,z)の変換行列
	$tRT->transpose($T);
	$tRT = $tRT->inverse();

#逆空間ベクトル(a*i)の変換行列: $tRT
	$RT->transpose($T);
	$RT = $RT->inverse();

#逆空間座標(h k l)の変換行列: $T
	$tR = $T;

	return ($sT, $T, $tRT, $RT, $tR);
}

sub GetBravaisToPrimitiveCellMatrix
{
	my ($this) = @_;

	my $SPG        = $this->GetCSpaceGroup();
	my $SPGName    = $SPG->SPGName();

	my $sT = new Math::MatrixReal(3, 3);
	my $T;
	if($SPGName =~ /^\s*F/i) {
		$sT = Math::MatrixReal->new_from_cols( 
			[ [ 0.5, 0.5,   0],
			  [ 0,   0.5, 0.5],
			  [ 0.5,   0, 0.5] ] );
	}
	elsif($SPGName =~ /^\s*I/i) {
		$sT = Math::MatrixReal->new_from_cols( 
			[ [-0.5, 0.5, 0.5],
			  [ 0.5,-0.5, 0.5],
			  [ 0.5, 0.5,-0.5] ] );
	}
	elsif($SPGName =~ /^\s*A/i) {
		$sT = Math::MatrixReal->new_from_cols( 
			[ [ 1.0, 0.0, 0.0],
			  [ 0.0, 0.5, 0.5],
			  [ 0.0,-0.5, 0.5] ] );
	}
	elsif($SPGName =~ /^\s*B/i) {
		$sT = Math::MatrixReal->new_from_cols( 
			[ [ 0.5, 0.0, 0.5],
			  [ 0.0, 1.0, 0.0],
			  [-0.5, 0.0, 0.5] ] );
	}
	elsif($SPGName =~ /^\s*C/i) {
		$sT = Math::MatrixReal->new_from_cols( 
			[ [ 0.5, 0.5, 0.0],
			  [-0.5, 0.5, 0.0],
			  [ 0.0, 0.0, 1.0] ] );
	}
	$T = $sT;
	$T->transpose($T);
	
	return ($sT, $T);
}

sub GetPrimitiveCell
{
	my ($this, $IsPrint) = @_;

#実空間べクトル(ai)の変換行列
	my ($sT, $T) = $this->GetBravaisToPrimitiveCellMatrix();
if($IsPrint) {
	print "Conversion matrix for real space vector: (a'i) = (Tij)(ai)\n";
	printf "   |%10.4f %10.4f %10.4f|\n", $T->element(1,1), $T->element(1,2), $T->element(1,3);
	printf "   |%10.4f %10.4f %10.4f|\n", $T->element(2,1), $T->element(2,2), $T->element(2,3);
	printf "   |%10.4f %10.4f %10.4f|\n", $T->element(3,1), $T->element(3,2), $T->element(3,3);
}
	return $this->ConvertLattice($T, 0, 1);
}

sub ConvertByMatrix
{
	my ($this, $x, $y, $z, $T) = @_;
	
	my $x1 = $T->element(1,1)*$x + $T->element(1,2)*$y + $T->element(1,3)*$z;
	my $y1 = $T->element(2,1)*$x + $T->element(2,2)*$y + $T->element(2,3)*$z;
	my $z1 = $T->element(3,1)*$x + $T->element(3,2)*$y + $T->element(3,3)*$z;
	return ($x1, $y1, $z1);
}

sub PrimitiveLatticeVectors
{
	my ($this, $UseAtomicUnit) = @_;
	my $k = ($UseAtomicUnit)? 1.0 / 0.529177 : 1.0;

	my ($ax,$ay,$az,$bx,$by,$bz,$cx,$cy,$cz) = $this->LatticeVectors($UseAtomicUnit);
	my ($sT, $T) = $this->GetBravaisToPrimitiveCellMatrix();
	my ($ax1, $ay1, $az1) = $this->ConvertByMatrix($ax, $ay, $az, $sT);
	my ($bx1, $by1, $bz1) = $this->ConvertByMatrix($bx, $by, $bz, $sT);
	my ($cx1, $cy1, $cz1) = $this->ConvertByMatrix($cx, $cy, $cz, $sT);
#print "vector ($ax1, $ay1, $az1, $bx1, $by1, $bz1, $cx1, $cy1, $cz1)\n";
#exit;
	return ($ax1, $ay1, $az1, $bx1, $by1, $bz1, $cx1, $cy1, $cz1);
}

1;