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

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

use strict;
#use Math::Matrix;
use Math::MatrixReal;

use Sci qw($pi $a0 $torad $todeg);# acos asin);
use Sci::Algorism;
use Math::Complex;

use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::Rietan;

#==================================================================
# 静的メンバー関数
#==================================================================
sub RoundSymmetricPosition
{
	my ($x) = @_;
	return  1.0 / 3.0 if($x ==  0.3333 or $x ==  0.333);
	return -1.0 / 3.0 if($x == -0.3333 or $x == -0.333);
	return  2.0 / 3.0 if($x ==  0.6667 or $x ==  0.667);
	return -2.0 / 3.0 if($x == -0.6667 or $x == -0.667);
	return  1.0 / 6.0 if($x ==  0.1667 or $x ==  0.167);
	return -1.0 / 6.0 if($x == -0.1667 or $x == -0.167);
	return  5.0 / 6.0 if($x ==  0.8333 or $x ==  0.833);
	return -5.0 / 6.0 if($x == -0.8333 or $x == -0.833);
	return $x;
}

sub PolarizationFactor
{
	my ($Theta, $MonochromaterAlpha) = @_;
	my $cos2a = 1.0;
	if(defined $MonochromaterAlpha) {
		$cos2a = cos(2.0 * $MonochromaterAlpha * $torad);
		$cos2a = $cos2a * $cos2a;
	}
	my $cos2Q = cos($torad * 2.0 * $Theta);
	return 1.0 + $cos2a * $cos2Q * $cos2Q;
}

sub LorentzFactor
{
	my ($Theta) = @_;
	my $cosQ  = cos($torad * $Theta);
	my $sinQ = sin($torad * $Theta);
#	my $sinQ = sin($torad * 2.0 * $Theta);
	return 1.0 / $cosQ / $sinQ / $sinQ;

#	my $cosQ  = cos($torad * $Theta);
	my $sin2Q = sin($torad * 2.0 * $Theta);
	return $cosQ / $sin2Q / $sin2Q;
}

sub TemperatureFactor
{
# Biso in Angstrom, $s in nm^-1
	my ($Biso, $s) = @_;
	$s = 0.1 * $s;
	my $M = $Biso * $s * $s;
	return exp(-$M);
}

#==================================================================
# 変数取得関数
#==================================================================
sub SampleName     { return shift->{SampleName}; }
sub SetSampleName  { my ($this,$s)=@_; return $this->{SampleName} = $s; }
sub SetCrystalName { my ($this,$name) = @_; return $this->{CrystalName} = $name; }
sub CrystalName    { my ($this) = @_; return $this->{CrystalName}; }
sub SPGName        { my ($this) = @_; my $SPG = $this->GetCSpaceGroup(); return $SPG->SPGName(); }
sub iSPG           { my ($this) = @_; my $SPG = $this->GetCSpaceGroup(); return $SPG->iSPG(); }
sub iSet           { my ($this) = @_; my $SPG = $this->GetCSpaceGroup(); return $SPG->iSet(); }
sub FormulaUnit    { return shift->{'FormulaUnit'}; }
sub SetFormulaUnit { my($this,$Z)=@_; return $this->{'FormulaUnit'} = $Z; }
sub ChemicalComposition       { return shift->{'ChemicalComposition'}; }
sub SetChemicalComposition    { my($this,$c)=@_; return $this->{'ChemicalComposition'} = $c; }
sub SumChemicalComposition    { return shift->{'SumChemicalComposition'}; }
sub SetSumChemicalComposition { my($this,$c)=@_; return $this->{'SumChemicalComposition'} = $c; }
sub SetnAtomType              { my ($this,$n) = @_; return $this->{nAtomType} = $n; }
sub nAtomType                 { my ($this) = @_; return $this->{nAtomType}; }
sub SetnAsymmetricAtomSite    { my ($this,$n) = @_; return $this->{nAsymmetricAtomSite} = $n; }
sub nAsymmetricAtomSite       { my ($this) = @_; return $this->{nAsymmetricAtomSite}; }
sub SetnTotalExpandedAtomSite { my ($this,$n) = @_; return $this->{nTotalExpandedAtomSite} = $n; }
sub nTotalExpandedAtomSite { my ($this) = @_; return $this->{nTotalExpandedAtomSite}; }

sub SetCAsymmetricAtomSiteList
{
	my ($this, @atomlist) = @_;
	$this->{nAsymmetricAtomSite} = @atomlist;
	return $this->{RefAsymmetricAtomSiteList} = \@atomlist;
}
sub GetBravaisLattice { my ($this) = @_; my $SPG = $this->GetCSpaceGroup();
	return $SPG->GetBravaisLattice();
}

sub SetCAtomTypeList 
{ 
	my ($this, @atomlist) = @_; 
	$this->{nAtomType} = @atomlist; 
	return $this->{RefAtomTypeList} = \@atomlist;
}

sub LatticeSystem  {
	my ($this) = @_;
	return $this->{LatticeSystem} if($this->{LatticeSystem});
	my $SPG = $this->GetCSpaceGroup();
	return $SPG->LatticeSystem();

}

sub SetLatticeSystem
{
	my ($this, $latticesystem, $SearchByLatticeParameter,  $tollat, $tolangle) = @_;
	my $SPG = $this->GetCSpaceGroup();
	my $ls;
#	my $ls = $SPG->LatticeSystem();
#	if($ls eq '') {
		$SPG->SetLatticeSystem($latticesystem, $SearchByLatticeParameter, $tollat, $tolangle);
		$ls = $SPG->LatticeSystem();
#	}
	$this->{'LatticeSystem'} = $ls;
	return $ls;
}

sub GetCAtomTypeList
{
	my ($this) = @_;
	my $RefAtomList = $this->{"RefAtomTypeList"};
#print "Ref:$RefAtomList<br>\n";
	return () if($RefAtomList eq '');
	return @$RefAtomList;
}

sub GetCAtomType
{
	my ($this,$i) = @_;
	my $RefAtomList = $this->{"RefAtomTypeList"};
	my $atomtype =  $RefAtomList->[$i];
	if(!$atomtype) {
		return 0;
	}
	return $atomtype;
}

sub GetCAsymmetricAtomSiteList
{
	my ($this) = @_;
	my $RefAtomList = $this->{"RefAsymmetricAtomSiteList"};
	return () if($RefAtomList eq '');
	return @$RefAtomList;
}

sub Volume
{
	my ($this,$UseAtomicUnit) = @_;
	my $k = 1.0 / 0.529177;
	$k = 1.0 unless($UseAtomicUnit);
	$k = $k*$k*$k;

	return $this->{"Volume"} * $k;
}

sub SetSpaceGroup
{
	my ($this, $SPGName, $iSPG, $iSet) = @_;
	my $SPG = new SpaceGroup();
	$SPG->SetSPGName($SPGName);
	$SPG->SetiSPG($iSPG);
	$SPG->SetiSet($iSet) if(!defined $iSet);
	$SPG->AddSymmetryOperation("x,y,z");

	return $this->SetCSpaceGroup($SPG);
}

sub SetCSpaceGroup
{
	my ($this, $SPG) = @_;
	$SPG->SetLatticeParameters($this->LatticeParameters())
		if(!$SPG->{a});
	$this->{"SpaceGroup"} = $SPG;
	$this->SetLatticeSystem() if(!$SPG->LatticeSystem());
	return $SPG;
}

sub GetCSpaceGroup
{
	my ($this) = @_;
	if(not defined $this->{"SpaceGroup"} or $this->{"SpaceGroup"} eq '') {
		my $SPG = new SpaceGroup();
		$SPG->SetP1();
		$this->{"SpaceGroup"} = $SPG;
	}
	return $this->{"SpaceGroup"};
}

sub GetCAtomSite
{
	my ($this,$i) = @_;
	my $RefAtomList = $this->{RefAsymmetricAtomSiteList};
	my $atomsite =  $RefAtomList->[$i];
	if(!$atomsite) {
		return 0;
	}
	return $atomsite;
}

sub GetCExpandedAtomSiteList
{
	my ($this) = @_;
	my $RefAtomList = $this->{RefExpandedAtomSiteList};
	return () if($RefAtomList eq '');
	return @$RefAtomList;
}

sub GetCExpandedAtomSite
{
	my ($this,$i) = @_;
	my $RefAtomList = $this->{RefExpandedAtomSiteList};
	my $atomsite =  $RefAtomList->[$i];
	if(!$atomsite) {
		return 0;
	}
	return $atomsite;
}

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

	my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
	$alpha = 90 if(abs($alpha-90.0) < $tolangle);
	$beta  = 90 if(abs($beta -90.0) < $tolangle);
	$gamma = 90 if(abs($gamma-90.0) < $tolangle);
	$this->LatticeParameters($a, $b, $c, $alpha, $beta, $gamma);

	$this->SetLatticeSystem('', 1, $tollatt, $tolangle);
	my $SPG = $this->GetCSpaceGroup();
	my $ls  = $SPG->LatticeSystem();

	if($ls eq 'cubic') {
		my $ar = ($a+$b+$c) / 3.0;
		$this->SetLatticeParameters($ar,$ar,$ar,90,90,90);
	}
	elsif($ls eq 'rhombohedral' or $ls eq 'trigonal') {
		my $ar = ($a+$b+$c) / 3.0;
		my $al = ($alpha+$beta+$gamma) / 3.0;
		$this->SetLatticeParameters($ar,$ar,$ar,$alpha,$alpha,$alpha);
	}
	elsif($ls eq 'tetragonal') {
		my $ar = ($a+$b) / 2.0;
		$this->SetLatticeParameters($ar,$ar,$c,90,90,90);
	}
	elsif($ls eq 'hexagonal') {
		my $ar = ($a+$b) / 2.0;
		$this->SetLatticeParameters($ar,$ar,$c,90,90,120);
	}
	elsif($ls eq 'orthorhombic') {
		$this->SetLatticeParameters($a,$b,$c,90,90,90);
	}

	my @site = $this->GetCAsymmetricAtomSiteList();
	for(my $i = 0 ; $i < @site ; $i++) {
		my ($x, $y, $z) = $site[$i]->Position();
		$x = $this->RoundParameter($x, $tolpos);
		$y = $this->RoundParameter($y, $tolpos);
		$z = $this->RoundParameter($z, $tolpos);
		$site[$i]->SetPosition($x, $y, $z);
	}

	$this->ExpandCoordinates();
}

sub SymmetrizeParameter
{
	my ($this, $ls, $x, $y, $z, $tol) = @_;

	$ls = $this->LatticeSystem() if(!defined $ls or $ls eq '');
	if($ls eq 'cubic' or $ls eq 'rhombohedral' or $ls eq 'trigonal') {
		$x = $y = $z = ($x+$y+$z) / 3.0;
	}
	elsif($ls eq 'tetragonal' or $ls eq 'hexagonal') {
		$x = $y = ($x+$y) / 2.0;
	}
	$x = $this->RoundParameter($x, $tol);
	$y = $this->RoundParameter($y, $tol);
	$z = $this->RoundParameter($z, $tol);

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

sub RoundParameter
{
	my ($this, $x, $tol) = @_;
	return $tol * int( ($x+0.1*$tol) / $tol );
}

sub GetCnMultiplicityExpandedAtomSiteList
{
	my ($this) = @_;
	my $RefnAtomList = $this->{"RefnExpandedAtomSite"};
	return () if($RefnAtomList eq '');
	return @$RefnAtomList;
}

#How to build Crystal (Ex. from CIF.pm)
#	my $crystal = new Crystal();
#	$crystal->SetLatticeParameters(
#		$this->LatticeParameters() );
#	my $SPG = $this->GetCSpaceGroup();
#	$crystal->SetCSpaceGroup($SPG);
#	$crystal->SetCAtomTypeList(
#		$this->GetCAtomTypeList() );
#	$crystal->SetCAsymmetricAtomSiteList(
#		$this->GetCAsymmetricAtomSiteList() );

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

require Math::Matrix;
require Math::MatrixReal;

	my $SPG = new SpaceGroup();
	$SPG->SetP1();
	$this->{"SpaceGroup"} = $SPG;

	my @AtomTypeList;
#	my $at1 = new AtomType();
#	push(@AtomTypeList, $at1);
#	pop(@AtomTypeList);
	$this->{"RefAtomTypeList"} = \@AtomTypeList;;
	$this->{"nAtomType"} = 0;

	my @AsAtomSiteList;
	my $at2 = new AtomSite();
	push(@AsAtomSiteList, $at2);
	pop(@AsAtomSiteList);
	$this->{"RefAsymmetricAtomSiteList"} = \@AsAtomSiteList;
	$this->{"nAsymmetricAtomSite"} = 0;

	my @ExAtomSiteList;
	my $at3 = new AtomSite();
	push(@ExAtomSiteList, $at3);
	pop(@ExAtomSiteList);
	$this->{"RefExpandedAtomSiteList"} = \@ExAtomSiteList;

	my @nAtomList;
	$this->{"RefnExpandedAtomSite"} = \@nAtomList;

	$this->{"nTotalExpandedAtomSite"}  = 0;

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#==================================================================
# 一般メンバー関数
#==================================================================
sub PreferentialOrientation
{
	my ($this, $IsSheet, $ho, $ko, $lo, $p1, $p2, $h, $k, $l) = @_;
	return 1.0 if($p1 == 1.0 or $p2 == 0.0);
	my @VO    = $this->GetVectorFromHKL($ho, $ko, $lo);
	my @Vhkl  = $this->GetVectorFromHKL($h, $k, $l);
	my $eVO   = Sci::NormalizeVector(\@VO);
	my $eVhkl = Sci::NormalizeVector(\@Vhkl);
	my $c = Sci::InnerProduct($eVO, $eVhkl);
	my $Q = acos(abs($c));
	$Q = $pi / 2.0 - $Q if(!$IsSheet);
#print "PO [$IsSheet]: $ho,$ko,$lo - $h,$k,$l => $Q\n";
	return $p1 + (1.0 - $p1) * exp(-$p2 * $Q * $Q);
}

sub IsVariablePosition
{
	my ($this, $x, $y, $z, $pos)= @_;

	my @x1 = ($x, $y, $z);
	my $c0 = new Crystal();
	$c0->SetLatticeParameters($this->LatticeParameters());
	$c0->SetCSpaceGroup($this->GetCSpaceGroup());
	$c0->AddAtomType('H');
	$c0->AddAtomSite('H1', 'H', $x1[0], $x1[1], $x1[2]);
	$c0->ExpandCoordinates();
	my @sites = $c0->GetCExpandedAtomSiteList();
	my $nSite = @sites;

	my $VariablePos = 1.0;
	$x1[$pos] += 0.14142;
	my $c1 = new Crystal();
	$c1->SetLatticeParameters($this->LatticeParameters());
	$c1->SetCSpaceGroup($this->GetCSpaceGroup());
	$c1->AddAtomType('H');
	$c1->AddAtomSite('H1', 'H', $x1[0], $x1[1], $x1[2]);
	$c1->ExpandCoordinates();
	@sites = $c1->GetCExpandedAtomSiteList();
	$VariablePos = 0.0 if($nSite < @sites);

	return $VariablePos;
}

sub SetMinimalLatticeParameters
{
	my ($this, @a) = @_;
	my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
	if($ls eq 'cubic') {
		$this->SetLatticeParameters($a[0], $a[0], $a[0], 90.0, 90.0, 90.0);
	}
	elsif($ls =~ /tetra/) {
		$this->SetLatticeParameters($a[0], $a[0], $a[1], 90.0, 90.0, 90.0);
	}
	elsif($ls =~ /ortho/) {
		$this->SetLatticeParameters($a[0], $a[1], $a[2], 90.0, 90.0, 90.0);
	}
	elsif($ls =~ /(rhomb|trig)/) {
		$this->SetLatticeParameters($a[0], $a[0], $a[0], $a[1], $a[1], $a[1]);
	}
	elsif($ls =~ /hex/) {
		$this->SetLatticeParameters($a[0], $a[0], $a[1], 90.0, 90.0, 120.0);
	}
	elsif($ls =~ /mono/) {
		$this->SetLatticeParameters($a[0], $a[1], $a[2], 90.0, $a[3], 90.0);
	}
	else {
		$this->SetLatticeParameters($a[0], $a[1], $a[2], $a[3], $a[4], $a[5]);
	}
	my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
#print "latt: $a $b $c $alpha $beta $gamma\n";
	return ($a,$b,$c,$alpha,$beta,$gamma);
}

sub SetLatticeParameters
{
	my ($this, $a, $b, $c, $alpha, $beta, $gamma) = @_;
#	return if(!defined $a or $a == 0);

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

	$this->CalcMetrics();

	return ($a, $b, $c, $alpha, $beta, $gamma);
}

sub LatticeParameters
{
	my ($this,$UseAtomicUnit) = @_;

	my $a = $this->{'a'};
#$a = 1.0 if(not defined $a);
	my $b = $this->{'b'};
#$b = $a if(not defined $b);
	my $c = $this->{'c'};
#$c = $b if(not defined $c);
	my $alpha = $this->{'alpha'};
#$alpha = 90.0 if(not defined $alpha);
	my $beta  = $this->{'beta'};
#$beta  = 90.0 if(not defined $beta);
	my $gamma = $this->{'gamma'};
#$gamma = 90.0 if(not defined $gamma);

	my $k = 1.0 / 0.529177;
	$k = 1.0 unless($UseAtomicUnit);

	return ($k*$a, $k*$b, $k*$c, $alpha, $beta, $gamma);
}

sub GetReciprocalDistanceFromK
{
	my ($this, $kx0, $ky0, $kz0, $kx1, $ky1, $kz1) = @_;
#printf "k: (%6.3g,%6.3g,%6.3g) - (%6.3g,%6.3g,%6.3g): ",
#	$kx0,$ky0,$kz0,$kx1,$ky1,$kz1;

	$kx0 = $kx1 - $kx0 if(defined $kx1);
	$ky0 = $ky1 - $ky0 if(defined $ky1);
	$kz0 = $kz1 - $kz0 if(defined $kz1);
#printf "(%6.3g,%6.3g,%6.3g)\n", $kx0, $ky0, $kz0;

if(0) {
	my $Ra = $this->{'Ra'};
	my $Rb = $this->{'Rb'};
	my $Rc = $this->{'Rc'};
#print "Rec lat: $Ra $Rb $Rc\n";
	$kx0 = $kx0 / $Ra;
	$ky0 = $ky0 / $Rb;
	$kz0 = $kz0 / $Rc;
}

	my $Rg11 = $this->{"Rg11"};
	my $Rg12 = $this->{"Rg12"};
	my $Rg13 = $this->{"Rg13"};
	my $Rg21 = $this->{"Rg21"};
	my $Rg22 = $this->{"Rg22"};
	my $Rg23 = $this->{"Rg23"};
	my $Rg31 = $this->{"Rg31"};
	my $Rg32 = $this->{"Rg32"};
	my $Rg33 = $this->{"Rg33"};
#print "Rg: ($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)<br>\n";
#print "($h,$k,$l)-($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)<br>\n";
	my $d = $Rg11*$kx0*$kx0 + $Rg22*$ky0*$ky0 + $Rg33*$kz0*$kz0
	      + 2.0*$Rg12*$kx0*$ky0 + 2.0*$Rg23*$ky0*$kz0 + 2.0*$Rg13*$kx0*$kz0;
	$d = 0.0 if($d < 0.0);
	$d = sqrt($d);
#print "d=$d<br>\n";
#print "enter: $d ($h,$k,$l),($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)<br>";
	return $d;
}

sub CalculateHKLInterplanarSpacing
{
	my ($this, $h, $k, $l) = @_;
	my $Rg11 = $this->{"Rg11"};
	my $Rg12 = $this->{"Rg12"};
	my $Rg13 = $this->{"Rg13"};
	my $Rg21 = $this->{"Rg21"};
	my $Rg22 = $this->{"Rg22"};
	my $Rg23 = $this->{"Rg23"};
	my $Rg31 = $this->{"Rg31"};
	my $Rg32 = $this->{"Rg32"};
	my $Rg33 = $this->{"Rg33"};
#print "($h,$k,$l)-($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)<br>\n";
	my $d = $Rg11*$h*$h + $Rg22*$k*$k + $Rg33*$l*$l
	      + 2.0*$Rg12*$h*$k + 2.0*$Rg23*$k*$l + 2.0*$Rg13*$l*$h;
#print "d=$d<br>";
#print "enter: $d ($h,$k,$l),($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)<br>";
	return 0.0 if($d == 0.0);
	$d = sqrt(1.0 / $d);
#print "d=$d<br>";
	return $d;
}

sub CalculateDiffractionAngleFromInterplanarSpacing
{
	my ($this, $wl, $d) = @_;
	return 0.0 if($d == 0.0);
	my $sq = $wl * 0.5 / $d;
#print "sq: $sq<br>";
	return 0.0 if(abs($sq) > 1.0);
	my $Q = Sci::dasin($sq);
	return $Q;
}

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

	my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
if($a == 0.0) {
	print "Error: CrystalObject.pm::CalcMetrics: Invalid Lattice Parameters.<br>\n";
	print "Lattice: $a $b $c $alpha $beta $gamma<br>\n";
}
	my $torad = Sci::torad(); #3.1415926535 / 180.0;
	$this->{g11} = $a * $a;
	$this->{g12} = $this->{g21} = $a * $b * cos($torad * $gamma);
	$this->{g12} = $this->{g21} = 0.0 if(abs($this->{"g12"}) < 1.0e-7);
	$this->{g13} = $this->{g31} = $a * $c * cos($torad * $beta);
	$this->{g13} = $this->{g31} = 0.0 if(abs($this->{"g13"}) < 1.0e-7);
	$this->{g22} = $b * $b;
	$this->{g23} = $this->{g32} = $b * $c * cos($torad * $alpha);
	$this->{g23} = $this->{g32} = 0.0 if(abs($this->{"g23"}) < 1.0e-7);
	$this->{g33} = $c * $c;
	my $V = $this->CalculateVolume();
#print "V=$V<BR>";
#print "a=$a\n";
	my $ax = $a;
	my $ay = 0.0;
	my $az = 0.0;
	my $bx = $b * Sci::dcos($gamma);
	my $by = $b * Sci::dsin($gamma);
	my $bz = 0.0;
	my $cx = $c * Sci::dcos($beta);
	my $cy = ($c * Sci::dcos($alpha) - $cx * Sci::dcos($gamma)) / Sci::dsin($gamma);
	my $cz = sqrt($c*$c - $cx*$cx - $cy*$cy);
	$this->{a11} = $ax;
	$this->{a12} = $ay;
	$this->{a13} = $az;
	$this->{a21} = $bx;
	$this->{a22} = $by;
	$this->{a23} = $bz;
	$this->{a31} = $cx;
	$this->{a32} = $cy;
	$this->{a33} = $cz;

	my $sina = Sci::dsin($alpha);
	my $sinb = Sci::dsin($beta);
	my $sinc = Sci::dsin($gamma);
	my $cosa = Sci::dcos($alpha);
	my $cosb = Sci::dcos($beta);
	my $cosc = Sci::dcos($gamma);
	$this->{Rg11} = $b*$b*$c*$c*$sina*$sina / $V / $V;
	$this->{Rg12} = $this->{Rg21} = $a*$b*$c*$c * ($cosa*$cosb - $cosc) / $V / $V;
	$this->{Rg12} = $this->{Rg21} = 0.0 if(abs($this->{"g12"}) < 1.0e-7);
	$this->{Rg13} = $this->{Rg31} = $a*$b*$b*$c * ($cosc*$cosa - $cosb) / $V / $V;
	$this->{Rg13} = $this->{Rg31} = 0.0 if(abs($this->{"g13"}) < 1.0e-7);
	$this->{Rg22} = $a*$a*$c*$c*$sinb*$sinb / $V / $V;
	$this->{Rg23} = $this->{Rg32} = $a*$a*$b*$c * ($cosb*$cosc - $cosa) / $V / $V;
	$this->{Rg23} = $this->{Rg32} = 0.0 if(abs($this->{"g23"}) < 1.0e-7);
	$this->{Rg33} = $a*$a*$b*$b*$sinc*$sinc / $V / $V;

	my $Ra = $this->{Ra} = sqrt($this->{Rg11});
	my $Rb = $this->{Rb} = sqrt($this->{Rg22});
	my $Rc = $this->{Rc} = sqrt($this->{Rg33});
	my $Ralpha = $this->{Ralpha} = $todeg * acos($this->{Rg23} / $this->{Rb} / $this->{Rc});
	my $Rbeta  = $this->{Rbeta}  = $todeg * acos($this->{Rg31} / $this->{Rc} / $this->{Ra});
	my $Rgamma = $this->{Rgamma} = $todeg * acos($this->{Rg12} / $this->{Ra} / $this->{Rb});
#print "RLatt: $Ra, $Rb, $Rc,  $Ralpha, $Rbeta, $Rgamma\n";
	my $kV = 1.0 / $V;
	my ($Rax, $Ray, $Raz) = Sci::VectorProduct([$bx, $by, $bz], [$cx, $cy, $cz], $kV);
	my ($Rbx, $Rby, $Rbz) = Sci::VectorProduct([$cx, $cy, $cz], [$ax, $ay, $az], $kV);
	my ($Rcx, $Rcy, $Rcz) = Sci::VectorProduct([$ax, $ay, $az], [$bx, $by, $bz], $kV);
#my $cc = Sci::ScalarProduct([$Rax, $Ray, $Raz], [$Rax, $Ray, $Raz]);
#print "cc=$cc, $this->{Rg11}\n";
#$cc = Sci::ScalarProduct([$Rax, $Ray, $Raz], [$Rbx, $Rby, $Rbz]);
#print "cc=$cc, $this->{Rg12}\n";
#printf "Rx*1=(%7.3f, %7.3f, %7.3f)\n", $Rax, $Ray, $Raz;
#printf "Ry*1=(%7.3f, %7.3f, %7.3f)\n", $Rbx, $Rby, $Rbz;
#printf "Rz*1=(%7.3f, %7.3f, %7.3f)\n", $Rcx, $Rcy, $Rcz;
if(0) {
	my $Rax = $Ra;
	my $Ray = 0.0;
	my $Raz = 0.0;
	my $Rbx = $Rb * Sci::dcos($Rgamma);
	my $Rby = $Rb * Sci::dsin($Rgamma);
	my $Rbz = 0.0;
	my $Rcx = $Rc * Sci::dcos($Rbeta);
	my $Rcy = ($Rc * Sci::dcos($Ralpha) - $Rcx * Sci::dcos($Rgamma)) / Sci::dsin($Rgamma);
	my $Rcz = sqrt($Rc*$Rc - $Rcx*$Rcx - $Rcy*$Rcy);
}
	$this->{Ra11} = $Rax;
	$this->{Ra12} = $Ray;
	$this->{Ra13} = $Raz;
	$this->{Ra21} = $Rbx;
	$this->{Ra22} = $Rby;
	$this->{Ra23} = $Rbz;
	$this->{Ra31} = $Rcx;
	$this->{Ra32} = $Rcy;
	$this->{Ra33} = $Rcz;
#printf "Rx*2=(%7.3f, %7.3f, %7.3f)\n", $Rax, $Ray, $Raz;
#printf "Ry*2=(%7.3f, %7.3f, %7.3f)\n", $Rbx, $Rby, $Rbz;
#printf "Rz*2=(%7.3f, %7.3f, %7.3f)\n", $Rcx, $Rcy, $Rcz;
	return 1;
}

sub LatticeVectors
{
	my ($this,$UseAtomicUnit) = @_;
	my $k = 1.0 / 0.529177;
	$k = 1.0 unless($UseAtomicUnit);

	return ($k*$this->{a11}, $k*$this->{a12}, $k*$this->{a13},
			$k*$this->{a21}, $k*$this->{a22}, $k*$this->{a23},
			$k*$this->{a31}, $k*$this->{a32}, $k*$this->{a33});
}

sub ReciprocalLatticeVectors
{
	my ($this,$UseAtomicUnit) = @_;
	my $k = 0.529177;
	$k = 1.0 unless($UseAtomicUnit);

	return ($k*$this->{Ra11}, $k*$this->{Ra12}, $k*$this->{Ra13},
			$k*$this->{Ra21}, $k*$this->{Ra22}, $k*$this->{Ra23},
			$k*$this->{Ra31}, $k*$this->{Ra32}, $k*$this->{Ra33});
}

sub GetPrimitiveCrystal
{
	my ($this, $LatticeParameterOnly, $CheckConsistency, $IsPrint) = @_;

	my $T = new Math::MatrixReal(3, 3);
	my $matrix = $this->GetMatrixConventionalToPrimitiveCell();
	$T = Math::MatrixReal->new_from_cols($matrix);
	$T->transpose($T);
	return $this->ConvertLattice($T, $LatticeParameterOnly, $CheckConsistency, $IsPrint);
}

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

	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 GetYRangeForIteration
{
	my ($this, $ix) = @_;

	my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
	if($ls eq 'cubic' or $ls =~ /(rhomb|trig)/) {
		return $ix;
	}
	elsif($ls =~ /(rhomb|trig)/) {
		return $ix;
	}
	elsif($ls =~ /tetra/ or $ls =~ /hex/) {
		return $ix;
	}
	elsif($ls =~ /ortho/) {
		return 0;
	}
	elsif($ls =~ /mono/) {
		return 0;
	}
	else {
		return 0;
	}
}

sub GetZRangeForIteration
{
	my ($this, $ix, $iy) = @_;

	my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
	if($ls eq 'cubic' or $ls =~ /(rhomb|trig)/) {
		return $iy;
	}
	elsif($ls =~ /tetra/ or $ls =~ /hex/) {
		return 0;
	}
	elsif($ls =~ /ortho/) {
		return 0;
	}
	elsif($ls =~ /mono/) {
		return 0;
	}
	else {
		return 0;
	}
}

sub ConvertLatticeIndexByLatticeSystem
{
	my ($this, $ix, $iy, $iz, $SortDirection) = @_;
	$SortDirection = 1 if(!defined $SortDirection);

	if($SortDirection == 0) {
		return (0, $ix, $iy, $iz);
	}
	
	my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
	my $IsConverted = 0;
	if($ls eq 'cubic' or $ls =~ /(rhomb|trig)/) {
		$IsConverted = 1 if($ix > $iy or $iy > $iz or $ix > $iz);
		if($SortDirection > 0) {
			return ($IsConverted, sort { return $a <=> $b; } ($ix, $iy, $iz) );
		}
		return ($IsConverted, sort { return $b <=> $a; } ($ix, $iy, $iz) );
	}
	elsif($ls =~ /tetra/ or $ls =~ /hex/) {
		$IsConverted = 1 if($ix > $iy);
		my @a;
		if($SortDirection > 0) {
			@a =  sort { return $a <=> $b; } ($ix, $iy);
		}
		else {
			@a =  sort { return $b <=> $a; } ($ix, $iy);
		}
		return ($IsConverted, @a, $iz);
	}
	else {
		die "\nCrystaObject::ConvertLatticeIndexByLatticeSystem: Error: Lattice system [$ls] is not supported.\n";
	}
	return ($IsConverted, $ix, $iy, $iz);
}

sub Metrics
{
	my ($this) = @_;
	my $g11 = $this->{g11};
	my $g12 = $this->{g12};
	my $g13 = $this->{g13};
	my $g21 = $this->{g21};
	my $g22 = $this->{g22};
	my $g23 = $this->{g23};
	my $g31 = $this->{g31};
	my $g32 = $this->{g32};
	my $g33 = $this->{g33};
	return ($g11, $g12, $g13,
			$g21, $g22, $g23,
			$g31, $g32, $g33)
}

sub GetPerpendicularHKL
{
	my ($this, $h1, $k1, $l1, $h2, $k2, $l2) = @_;

	my @v1 = $this->GetVectorFromHKL($h1, $k1, $l1);
	my @v2 = $this->GetVectorFromHKL($h2, $k2, $l2);
	my @v3 = Sci::VectorProduct(\@v1, \@v2);
#	my @v3 = ($v1[1]*$v2[2] - $v1[2]*$v2[1],
#			  $v1[2]*$v2[0] - $v1[0]*$v2[2],
#			  $v1[0]*$v2[1] - $v1[1]*$v2[0]);
#printf "v1=(%7.5f, %7.5f, %7.5f)\n", @v1;
#printf "v2=(%7.5f, %7.5f, %7.5f)\n", @v2;
#printf "v3=(%7.5f, %7.5f, %7.5f)\n", @v3;

	my $R = new Math::MatrixReal(3, 3);
	$R->assign(1, 1, $this->{Ra11});
	$R->assign(1, 2, $this->{Ra21});
	$R->assign(1, 3, $this->{Ra31});
	$R->assign(2, 1, $this->{Ra12});
	$R->assign(2, 2, $this->{Ra22});
	$R->assign(2, 3, $this->{Ra32});
	$R->assign(3, 1, $this->{Ra13});
	$R->assign(3, 2, $this->{Ra23});
	$R->assign(3, 3, $this->{Ra33});
#print "R=$R\n";
	my $RR = $R->inverse();
#print "RR=$RR\n";
	my @hkl = ($RR->element(1,1)*$v3[0] + $RR->element(1,2)*$v3[1] + $RR->element(1,3)*$v3[2],
			   $RR->element(2,1)*$v3[0] + $RR->element(2,2)*$v3[1] + $RR->element(2,3)*$v3[2],
			   $RR->element(3,1)*$v3[0] + $RR->element(3,2)*$v3[1] + $RR->element(3,3)*$v3[2]);
	my ($min, $max) = Utils::CalMinMax(\@hkl, 1);
#printf "hkl=%6.3f %6.3f %6.3f\n", @hkl;
	$min = 1.0e30;
	for(my $i = 0 ; $i < 3 ; $i++) {
		$hkl[$i] = 0.0 if(abs($hkl[$i]) < $max * 1.0e-15);
		$min = $hkl[$i] if($hkl[$i] != 0 and $min > $hkl[$i]);
	}
	if($min < 1.0e30) {
		for(my $i = 0 ; $i < 3 ; $i++) {
			$hkl[$i] /= $min;
		}
	}
#printf "hkl=%6.3f %6.3f %6.3f\n", @hkl;
#my @v32 =  $this->GetVectorFromHKL(@hkl);
#printf "v3=(%7.5f, %7.5f, %7.5f)\n", @v32;
	return @hkl;
}

sub GetVectorFromHKL
{
	my ($this, $h, $k, $l) = @_;
	my @V = ($this->{Ra11} * $h + $this->{Ra21} * $k + $this->{Ra31} * $l,
		     $this->{Ra12} * $h + $this->{Ra22} * $k + $this->{Ra32} * $l,
		     $this->{Ra13} * $h + $this->{Ra23} * $k + $this->{Ra33} * $l);
#print "hkl=$h,$k,$l: ", join(',', @V), ")\n";
	return @V;
}

sub RMetrics
{
	my ($this) = @_;
	my $Rg11 = $this->{Rg11};
	my $Rg12 = $this->{Rg12};
	my $Rg13 = $this->{Rg13};
	my $Rg21 = $this->{Rg21};
	my $Rg22 = $this->{Rg22};
	my $Rg23 = $this->{Rg23};
	my $Rg31 = $this->{Rg31};
	my $Rg32 = $this->{Rg32};
	my $Rg33 = $this->{Rg33};
	return ($Rg11, $Rg12, $Rg13,
		$Rg21, $Rg22, $Rg23,
		$Rg31, $Rg32, $Rg33)
}

sub SetLatticeVector
{
	my ($this, $a11, $a12, $a13,
	    $a21, $a22, $a23, $a31, $a32, $a33) = @_;
	$this->{a11} = $a11;
	$this->{a12} = $a12;
	$this->{a13} = $a13;
	$this->{a21} = $a21;
	$this->{a22} = $a22;
	$this->{a23} = $a23;
	$this->{a31} = $a31;
	$this->{a32} = $a32;
	$this->{a33} = $a33;
	$this->CalMetricsFromLatticeVector();
	$this->CalLatticeParametersFromMetrics();
	return;
}

sub CalMetricsFromLatticeVector
{
	my ($this) = @_;
	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = (
		$this->{a11}, $this->{a12}, $this->{a13},
		$this->{a21}, $this->{a22}, $this->{a23},
		$this->{a31}, $this->{a32}, $this->{a33});

	$this->{g11} = $a11*$a11 + $a12*$a12 + $a13*$a13;
	$this->{g22} = $a21*$a21 + $a22*$a22 + $a23*$a23;
	$this->{g33} = $a31*$a31 + $a32*$a32 + $a33*$a33;
	$this->{g12} = $this->{g21} = $a11*$a21 + $a12*$a22 + $a13*$a33;
	$this->{g23} = $this->{g32} = $a21*$a31 + $a22*$a32 + $a23*$a33;
	$this->{g31} = $this->{g13} = $a31*$a11 + $a32*$a12 + $a33*$a13;
	return ($this->{g11}, $this->{g12}, $this->{g13},
			$this->{g21}, $this->{g22}, $this->{g23},
			$this->{g31}, $this->{g32}, $this->{g33})
}

sub CalculateLatticeConstantFromVector
{
	my ($this) = @_;
	my $a11 = $this->{a11};
	my $a12 = $this->{a12};
	my $a13 = $this->{a13};
	my $a21 = $this->{a21};
	my $a22 = $this->{a22};
	my $a23 = $this->{a23};
	my $a31 = $this->{a31};
	my $a32 = $this->{a32};
	my $a33 = $this->{a33};
	
	my $a = sqrt($a11*$a11 + $a12*$a12 + $a13*$a13);
	my $b = sqrt($a21*$a21 + $a22*$a22 + $a23*$a23);
	my $c = sqrt($a31*$a31 + $a32*$a32 + $a33*$a33);
	my $ab = $a11*$a21 + $a12*$a22 + $a13*$a23;
	my $bc = $a21*$a31 + $a22*$a32 + $a23*$a33;
	my $ca = $a31*$a11 + $a32*$a12 + $a33*$a13;
	my $cosg = $ab / $a / $b;
	my $cosa = $bc / $b / $c;
	my $cosb = $ca / $c / $a;
	my $alpha = Sci::dacos($cosa);
	my $beta  = Sci::dacos($cosb);
	my $gamma = Sci::dacos($cosg);
	$a = Utils::Round($a+0.0, 6); #sprintf("%11.6f", $a) + 0.0;
	$b = Utils::Round($b+0.0, 6); #sprintf("%11.6f", $b) + 0.0;
	$c = Utils::Round($c+0.0, 6); #sprintf("%11.6f", $c) + 0.0;
	$alpha = Utils::Round($alpha+0.0, 6); #sprintf("%12.6f", $alpha) + 0.0;
	$beta  = Utils::Round($beta+0.0, 6); #sprintf("%12.6f", $beta) + 0.0;
	$gamma = Utils::Round($gamma+0.0, 6); #sprintf("%12.6f", $gamma) + 0.0;
#print "a=$a $b $c $alpha $beta $gamma\n";
	return ($a,$b,$c,$alpha,$beta,$gamma);
}

sub FindNearestEquivalentFractionCoordinate
{
	my ($this, $x2, $x) = @_;
	my $xret = $x2;
	my $mind = abs($x2 - $x);
	if($mind > abs($x2+1.0 - $x)) {
		$xret = $x2+1.0;
		$mind = abs($x2+1.0 - $x);
	}
	if($mind > abs($x2-1.0 - $x)) {
		$xret = $x2-1.0;
		$mind = abs($x2-1.0 - $x);
	}
	return $xret;
}

sub FractionalToCartesian
{
	my ($this, $x,$y,$z) = @_;
	
	my $a11 = $this->{a11};
	my $a12 = $this->{a12};
	my $a13 = $this->{a13};
	my $a21 = $this->{a21};
	my $a22 = $this->{a22};
	my $a23 = $this->{a23};
	my $a31 = $this->{a31};
	my $a32 = $this->{a32};
	my $a33 = $this->{a33};
	my $xc = $a11*$x + $a21*$y + $a31*$z;
	my $yc = $a12*$x + $a22*$y + $a32*$z;
	my $zc = $a13*$x + $a23*$y + $a33*$z;
	return ($xc,$yc,$zc);
}

sub CartesianToFractional
{
	my ($this, $xc,$yc,$zc) = @_;
	
	my $A = new Math::MatrixReal(3, 3);
	$A->assign(1, 1, $this->{"a11"});
	$A->assign(1, 2, $this->{"a21"});
	$A->assign(1, 3, $this->{"a31"});
#	$A->assign(1, 2, $this->{"a12"});
#	$A->assign(1, 3, $this->{"a13"});
	$A->assign(2, 1, $this->{"a12"});
#	$A->assign(2, 1, $this->{"a21"});
	$A->assign(2, 2, $this->{"a22"});
	$A->assign(2, 3, $this->{"a32"});
#	$A->assign(2, 3, $this->{"a23"});
	$A->assign(3, 1, $this->{"a13"});
	$A->assign(3, 2, $this->{"a23"});
#	$A->assign(3, 1, $this->{"a31"});
#	$A->assign(3, 2, $this->{"a32"});
	$A->assign(3, 3, $this->{"a33"});
#	print "A:<br>\n";print $A;print "<br>\n";

	my $RA = $A->inverse();
#	print "RA:<br>\n";print $RA;print "<br>\n";

	my $X = new Math::MatrixReal(3,1);
	$X->assign(1, 1, $xc);
	$X->assign(2, 1, $yc);
	$X->assign(3, 1, $zc);
#	print "X:<br>\n";print $X;print "<br>\n";
	
	my $XX = $RA * $X;
#	print "XX:<br>\n";print $XX;print "<br>\n";

	return ($XX->element(1,1), 
		$XX->element(2,1), 
		$XX->element(3,1));
}

sub CalculateVolume
{
	my ($this) = @_;
	my $a = $this->{"a"};
	my $b = $this->{"b"};
	my $c = $this->{"c"};
	my $alpha = $this->{"alpha"};
	my $beta  = $this->{"beta"};
	my $gamma = $this->{"gamma"};
	my $cosa = Sci::dcos($alpha);
	my $cosb = Sci::dcos($beta);
	my $cosg = Sci::dcos($gamma);
	my $vol = 1.0 - $cosa*$cosa - $cosb*$cosb - $cosg*$cosg 
		+ 2.0 * $cosa*$cosb*$cosg;
	$vol = $a*$b*$c*sqrt($vol);
	$vol = sprintf("%14.6g", $vol) + 0.0;
	$this->{"Volume"} = Sci::Round($vol, 4);
	return $vol;
}

sub FindAtomTypeByName
{
	my ($this, $atomname) = @_;
	$atomname =~ s/([\d+-]+)//;

	my $atomlistref = $this->{"RefAtomTypeList"};
	for(my $i = 0 ; $i < @$atomlistref ; $i++) {
		my $type = $atomlistref->[$i];
		my $name = $type->AtomNameOnly();
#print "name: $atomname <=> $name\n";
		return $type if(uc $name eq uc $atomname);
	}
	return undef;
}

sub AddAtomType
{
	my ($this, $atomname, $CheckRegistered) = @_;
	$CheckRegistered = 1 if(!defined $CheckRegistered);

	my $atomlistref = $this->{"RefAtomTypeList"};
	if($atomlistref) {
		if($CheckRegistered and $this->FindAtomTypeByName($atomname)) {
			return;
		}
		my $at = new AtomType();
		$at->SetAtomName($atomname);
		push(@$atomlistref, $at);

		$this->{"nAtomType"} = scalar @$atomlistref;
	}
	else {
		my @atomlist;
		my $at = new AtomType();
		$at->SetAtomName($atomname);
		push(@atomlist, $at);

		$this->{"nAtomType"} = @atomlist;
		$this->{"RefAtomTypeList"} = \@atomlist;
	}

#my $ref = $this->{"RefAtomTypeList"} ;
#print "ref: ", $ref->[$this->{"nAtomType"}-1]->AtomName(), "<BR>";
	return $this->{"nAtomType"};
}

sub GetCBravaisLatticeAsymmetricAtomSiteList
{
	my ($this) = @_;
	$this->ExpandCoordinates();
	my $SPG = $this->GetCSpaceGroup();
	my @ExpandedAtomSiteList   = $this->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();

	my (@Tx, @Ty, @Tz);
	my $nT = $SPG->nTranslation();
	for(my $it = 1 ; $it <= $nT ; $it++) {
		($Tx[$it-1], $Ty[$it-1], $Tz[$it-1]) = $SPG->TranslationVector($it);
	}

	my @AtomSite;
	my $c = 0;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $label     = $atom->Label();
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position();
		my $occ       = $atom->Occupancy();
#print "$i: $atomname $charge [$label]: ($x, $y, $z) [occ=$occ]\n";

		my $IsExist = 0;
		for(my $it = 1 ; $it < $nT ; $it++) {
			my ($tx, $ty, $tz) = ($Tx[$it], $Ty[$it], $Tz[$it]);
			my $x1 = Utils::Reduce01($x + $tx);
			my $y1 = Utils::Reduce01($y + $ty);
			my $z1 = Utils::Reduce01($z + $tz);
			for(my $j = 0 ; $j < $c ; $j++) {
				my $atom2 = $AtomSite[$j];
				my ($x2, $y2, $z2) = $atom2->Position();
				my $dis = $this->GetNearestInterAtomicDistance($x1, $y1, $z1, $x2, $y2, $z2);
#print "$i-$it-$j: ($x,$y,$z): ($x1,$y1,$z1)-($x2,$y2,$z2)=$dis\n";
				if($dis < 1.0e-6) {
					$IsExist = 1;
					last;
				}
			}
			last if($IsExist);
		}
		if(!$IsExist) {
			$AtomSite[$c] = $atom;
			$c++;
		}
	}

	return @AtomSite;
}

sub AddAtomSite
{
	my ($this, $label, $atomname, $x, $y, $z, $occ, $fx, $fy, $fz) = @_;
#print "AddAtomSite: F=$fx, $y, $fz\n";
#	return $this->{nAsymmetricAtomSite} if($atomname eq '');

	my $atomlistref = $this->{"RefAsymmetricAtomSiteList"};
	if($atomlistref) {
		my $at = new AtomSite();
		$at->SetLabel($label);
		$at->SetAtomName($atomname);
		$at->SetPosition($x, $y, $z);
		$at->SetForce($fx, $fy, $fz) if(defined $fz);
		$at->SetOccupancy($occ);
		$at->SetId(@$atomlistref + 1);
		push(@$atomlistref, $at);

		$this->{nAsymmetricAtomSite} = scalar @$atomlistref;
	}
	else {
		my @atomlist;
		my $at = new AtomSite();
		$at->SetLabel($label);
		$at->SetAtomName($atomname);
		$at->SetPosition($x, $y, $z);
		$at->SetForce($fx, $fy, $fz) if(defined $fz);
		$at->SetOccupancy($occ);
		$at->SetId(1);
		push(@atomlist, $at);

		$this->{nAsymmetricAtomSite}       = @atomlist;
		$this->{RefAsymmetricAtomSiteList} = \@atomlist;
	}
	$this->AddAtomType($atomname) if($atomname ne '');

	return $this->{nAsymmetricAtomSite};
}

sub GetNearestInterAtomicDistance
{
	my ($this, $x0, $y0, $z0, $x1, $y1, $z1) = @_;;
	$x0 = $this->FindNearestEquivalentFractionCoordinate($x0, $x1);
	$y0 = $this->FindNearestEquivalentFractionCoordinate($y0, $y1);
	$z0 = $this->FindNearestEquivalentFractionCoordinate($z0, $z1);
	return $this->GetInterAtomicDistance($x0, $y0, $z0, $x1, $y1, $z1);
}

sub GetInterAtomicDistance
{
	my ($this, $x0, $y0, $z0, $x1, $y1, $z1) = @_;;
	my $dx = $x1 - $x0;
	my $dy = $y1 - $y0;
	my $dz = $z1 - $z0;
	my $g11 = $this->{"g11"};
	my $g22 = $this->{"g22"};
	my $g33 = $this->{"g33"};
	my $g12 = 2.0* $this->{"g12"};
	my $g13 = 2.0* $this->{"g13"};
	my $g23 = 2.0* $this->{"g23"};
	my $r2 = $g11*$dx*$dx + $g22*$dy*$dy + $g33*$dz*$dz
	       + $g12*$dx*$dy + $g13*$dx*$dz + $g23*$dy*$dz;
	return sqrt($r2);
}

sub GetInterAtomicAngle
{
	my ($this, $x0, $y0, $z0, $x1, $y1, $z1,
		   $x2, $y2, $z2) = @_;;

	my $dis01 = $this->GetInterAtomicDistance($x0,$y0,$z0,$x1,$y1,$z1);
	return 0.0 if($dis01 == 0.0);
	my $dis02 = $this->GetInterAtomicDistance($x0,$y0,$z0,$x2,$y2,$z2);
	return 0.0 if($dis02 == 0.0);
#print "dis: $dis01, $dis02\n";
	my $dx01 = $x1 - $x0;
	my $dy01 = $y1 - $y0;
	my $dz01 = $z1 - $z0;
	my $dx02 = $x2 - $x0;
	my $dy02 = $y2 - $y0;
	my $dz02 = $z2 - $z0;
	my $g11 = $this->{"g11"};
	my $g22 = $this->{"g22"};
	my $g33 = $this->{"g33"};
	my $g12 = 2.0* $this->{"g12"};
	my $g13 = 2.0* $this->{"g13"};
	my $g23 = 2.0* $this->{"g23"};
	my $ip  = $g11*$dx01*$dx02 + $g22*$dy01*$dy02 + $g33*$dz01*$dz02
		+ $g12*$dx01*$dy02 + $g13*$dx01*$dz02
		+ $g12*$dy01*$dx02 + $g23*$dy01*$dz02
		+ $g13*$dz01*$dx02 + $g23*$dz01*$dy02;

	my $cosa = $ip / $dis01 / $dis02;
#print "cos: $cosa<br>";
	my $angle = Sci::dacos($cosa);
#print "cos: $cosa $angle<br>";
	$angle = 360.0 - $angle if($angle > 180.0001);

	return $angle;
}

#iAtomTypeは1から始まる
sub FindiAtomType
{
	my ($this, $atomname) = @_;

	$atomname =~ s/^(\D+)(\d.*)$/$1/;
	$atomname = uc $atomname;
	my @AtomTypeList = $this->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name = $AtomTypeList[$i]->AtomName();
		$name =~ s/^(\D+)(\d.*)$/$1/;
		return $i+1 if($atomname eq uc $name);
	}
	return -1;
}

sub FindIdenticalAsymmetricSite
{
	my ($this, $atomname, $x, $y, $z) = @_;
	$atomname =~ s/\d[\+\-]?//;
	my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
	for(my $ia = 0 ; $ia < $nAsymmetricAtomSite ; $ia++) {
		my $site = $AtomSiteList[$ia];
		my $name = $site->AtomNameOnly();
		next if(uc $name ne uc $atomname);
		my ($x1, $y1, $z1) = $site->Position();

		my $dis = $this->GetNearestInterAtomicDistance(
				$x, $y, $z, $x1, $y1, $z1);
		return $ia if($dis < 0.01);
	}
	return undef;
}

sub ExpandCoordinates
{
	my ($this, $DoTranslation) = @_;
	$DoTranslation = 1 if(!defined $DoTranslation);

#初期化
	$this->{nTotalExpandedAtomSite}  = 0;
	$this->{Density}  = 0.0;

#同一原子の判断半径を決める
	my ($a, $b, $c) = $this->LatticeParameters();
	$a = $b if($a < $b);
	$a = $c if($a < $c);
	my $OverlapCriteriaRadius = 0.001 * $a;

#非対称構造中のパラメータ
	my @AtomTypeList = $this->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
	my $SPG = $this->GetCSpaceGroup();
	my $nTranslation = $SPG->nTranslation();
	$nTranslation = 1 if(!$DoTranslation);
	my $nSymmetryOperation = $SPG->nSymmetryOperation();
	my $nTotalSymmetry = $nTranslation * $nSymmetryOperation;
#print "nTotalSymmetry: $nTotalSymmetry\n";
#座標展開を開始
	my @ExpandedAtomSiteList;
	my @nExpandedAtomSite;
	my $count = 0;
	for(my $ia = 0 ; $ia < $nAsymmetricAtomSite ; $ia++) {
		$nExpandedAtomSite[$ia] = 0;

		my $atom         = $AtomSiteList[$ia];
		my $label        = $atom->Label();
		my $AtomName     = $atom->AtomName();
		my $iAtomType    = $this->FindiAtomType($AtomName);
		my ($x,$y,$z)    = $atom->Position(1);
		my ($fx,$fy,$fz) = $atom->Force();
		$x = CrystalObject::RoundSymmetricPosition($x);
		$y = CrystalObject::RoundSymmetricPosition($y);
		$z = CrystalObject::RoundSymmetricPosition($z);
		my $Occupancy    = $atom->Occupancy();
		my ($vx,$vy,$vz) = $atom->Velocity();
		$atom->SetId($ia+1);

		$ExpandedAtomSiteList[$count] = new AtomSite();
		$ExpandedAtomSiteList[$count]->SetIdAsymmetricAtomSite($ia+1);
		$ExpandedAtomSiteList[$count]->SetiAtomType($iAtomType);
		$ExpandedAtomSiteList[$count]->SetLabel($label);
		$ExpandedAtomSiteList[$count]->SetAtomName($AtomName);
		$ExpandedAtomSiteList[$count]->SetPosition($x,$y,$z);
		$ExpandedAtomSiteList[$count]->SetForce($fx,$fy,$fz) if(defined $fx);
		$ExpandedAtomSiteList[$count]->SetOccupancy($Occupancy);
		$ExpandedAtomSiteList[$count]->SetVelocity($vx,$vy,$vz) if(defined $vx);

		$nExpandedAtomSite[$ia] += 1;
		$count++;
#print "($ia) x,y,z: $x,$y,$z<br>\n";
		for(my $is = 0 ; $is < $nTotalSymmetry ; $is++) {
			my ($xs, $ys, $zs) = $SPG->DoSymmetryOperation($is, $x, $y, $z, 1);
#print "        ($is) xs,ys,zs: $xs,$ys,$zs<br>\n";
			my ($fxs, $fys, $fzs) = $SPG->DoVelocitySymmetryOperation($is, $fx, $fy, $fz, 1);
			my ($vxs, $vys, $vzs) = $SPG->DoVelocitySymmetryOperation($is, $vx, $vy, $vz, 1);

			my $IsPassed = 1;
			for(my $it = 0 ; $it < $count ; $it++) {
				my $iat = $ExpandedAtomSiteList[$it]->iAtomType();
				next if($iat != $iAtomType);
				my ($x0,$y0,$z0) = $ExpandedAtomSiteList[$it]->Position(1);

				my $dis = $this->GetNearestInterAtomicDistance(
							$x0, $y0, $z0, $xs, $ys, $zs);
#print "$label(ia=$ia:is=$is:it=$it): ";
#printf "(%7.4f,%7.4f,%7.4f)-(%7.4f,%7.4f,%7.4f): dis: %f (%f)<br>\n",
#	$x0, $y0, $z0, $xs, $ys, $zs, $dis, $OverlapCriteriaRadius;
				if($dis < $OverlapCriteriaRadius) {
					$IsPassed = 0;
					last;
				}
			}
			next unless($IsPassed);
			$ExpandedAtomSiteList[$count] = new AtomSite();
			$ExpandedAtomSiteList[$count]->SetIdAsymmetricAtomSite($ia+1);
			$ExpandedAtomSiteList[$count]->SetiAtomType($iAtomType);
			$ExpandedAtomSiteList[$count]->SetLabel($label);
			$ExpandedAtomSiteList[$count]->SetAtomName($AtomName);
			$ExpandedAtomSiteList[$count]->SetPosition($xs,$ys,$zs);
			$ExpandedAtomSiteList[$count]->SetForce($fxs,$fys,$fzs);
			$ExpandedAtomSiteList[$count]->SetOccupancy($Occupancy);
			$ExpandedAtomSiteList[$count]->SetVelocity($vxs,$vys,$vzs);
#Idは1からはじまる
			$ExpandedAtomSiteList[$count]->SetId($ia+1);
#print "v: ($vx,$vy,$vz) ($vxs, $vys, $vzs)<br>\n";
			$nExpandedAtomSite[$ia] += 1;
			$count++;
		}
	}
	$this->{nTotalExpandedAtomSite}  = $count;
	$this->{RefExpandedAtomSiteList} = \@ExpandedAtomSiteList;
	$this->{RefnExpandedAtomSite}    = \@nExpandedAtomSite;

#多重度を、非対称原子に設定
	for(my $ia = 0 ; $ia < @AtomSiteList ; $ia++) {
		my $atom = $AtomSiteList[$ia];
#print "atom: ($ia): m=$nExpandedAtomSite[$ia]<br>\n";
		$atom->SetMultiplicity($nExpandedAtomSite[$ia]);
	}
#多重度を、全原子に設定
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom = $ExpandedAtomSiteList[$i];
		my $ia   = $atom->Id() - 1;
#print "atom: ($i,$ia): m=$nExpandedAtomSite[$ia]<br>\n";
		$atom->SetMultiplicity($nExpandedAtomSite[$ia]);
	}

	my $Density = $this->CalculateDensity();
	$this->{Density} = $Density;

	$this->AnalyzeChemicalComposition();

	return;
}

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

	my @ExpandedAtomSiteList   = $this->GetCExpandedAtomSiteList();
	my %AtomCount;
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $site = $ExpandedAtomSiteList[$i];
		my $name = $site->AtomNameOnly(1);
#		$name =~ s/\{[^\{\}]*\}//g;
		$AtomCount{$name}++;
#print "$i: $name [$AtomCount{$name}]\n";
	}

	my @Names = keys %AtomCount;
	my $Composition = '';
	my $MinN = 1000000;
	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next if($name =~ /^(O|S|Se|Te|F|Cl|I|N|As)$/i);
		my $n = $AtomCount{$name};
		$MinN = $n if($MinN > $n);
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
#print "i=$i: [$n] [$MinN] [$name]\n";
	}
	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next unless($name =~ /^(O|S|Se|Te|F|Cl|I|N|As)$/i);
		my $n = $AtomCount{$name};
		$MinN = $n if($MinN > $n);
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
	}
	$Composition =~ s/^\s*//;
	$this->SetSumChemicalComposition($Composition);

	my @elements = Algorism::Factorize($MinN);
	my $fac = 1;
#print "$MinN = ";
	for(my $i = 0 ; $i < @elements ; $i++) {
		my $num = $elements[$i];
#print " * "  if($i > 0);
#print "$num";

		my $IsCommonFactor = 1;
		for(my $j = 0 ; $j < @Names ; $j++) {
			my $name = $Names[$j];
			my $n = $AtomCount{$name};
			unless(Algorism::IsFactor($n, $num*$fac)) {
				$IsCommonFactor = 0;
			}
		}
		$fac *= $num if($IsCommonFactor);
#print "fac[$i]=$fac\n";
	}
#print "\n";
	$this->SetFormulaUnit($fac);
#print "Z = $fac\n";

	$Composition = '';
	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next if($name =~ /^(O|S|Se|Te|F|Cl|I|N|As|)$/i);
		my $n = $AtomCount{$name} / $fac;
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
	}
	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next unless($name =~ /^(O|S|Se|Te|F|Cl|I|N|As|)$/i);
		my $n = $AtomCount{$name} / $fac;
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
	}
	$Composition =~ s/^\s*//;
	$this->SetChemicalComposition($Composition);
}

sub CalculateDensity
{
	my ($this) = @_;
	my $vol = $this->Volume();
#	return $this->{Density} if($vol == 0.0);

	my $Density = 0.0;
	my @AtomTypeList = $this->GetCAtomTypeList();
	my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = $this->{'nAsymmetricAtomSite'};
	my @nMultiplicityExpandedAtomSiteList 
		= $this->GetCnMultiplicityExpandedAtomSiteList();
	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
		my $AtomSite = $AtomSiteList[$i];
		my $AtomName = $AtomSite->AtomNameOnly(1);
		my $iAtomType = $this->FindiAtomType($AtomName);
		my $AtomType = $AtomTypeList[$iAtomType-1];
		my $occupancy = $AtomSite->Occupancy();
#print "i=$i: $AtomSite, $AtomName, $iAtomType, $AtomType, $occupancy\n";
		my $AtomicMass = $AtomType->AtomicMass();
		my $mult = $nMultiplicityExpandedAtomSiteList[$i];
#print "$i: $AtomName ($iAtomType): $AtomicMass\n";
#print "$i: $AtomName ($AtomicMass) [$mult]\n";

		$AtomicMass = 0.0 if(not defined $AtomicMass);
		$occupancy = 1.0 if(not defined $occupancy);
		$mult = 1.0 if(not defined $mult);
		$Density += $AtomicMass*$mult*$occupancy;
	}

	my $NA = Sci::NA(); #6.0221367e23;
#print "vol: $vol  NA: $NA\n";
	return $this->{"Density"} = Sci::Round($Density / $NA / $vol * 1.0e24, 4);
}

sub SetDensity
{
	my ($this, $density) = @_;
	return $this->{"Density"} = $density;
}

sub Density
{
	my ($this) = @_;
	return $this->{"Density"};
}

sub ConvertLattice
{
	my ($this, $T, $LatticeParameterOnly, $CheckConsistency, $IsPrint) = @_;
	my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();

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

	my $g = new Math::MatrixReal(3, 3);
	$g->assign(1, 1, $this->{'g11'});
	$g->assign(1, 2, $this->{'g12'});
	$g->assign(1, 3, $this->{'g13'});
	$g->assign(2, 1, $this->{'g21'});
	$g->assign(2, 2, $this->{'g22'});
	$g->assign(2, 3, $this->{'g23'});
	$g->assign(3, 1, $this->{'g31'});
	$g->assign(3, 2, $this->{'g32'});
	$g->assign(3, 3, $this->{'g33'});
	my $g2 = new Math::MatrixReal(3, 3);
	for(my $i1 = 1 ; $i1 <= 3 ; $i1++) {
	for(my $i2 = 1 ; $i2 <= 3 ; $i2++) {
		my $val = 0.0;
		for(my $j1 = 1 ; $j1 <= 3 ; $j1++) {
		for(my $j2 = 1 ; $j2 <= 3 ; $j2++) {
			$val += $g->element($j1,$j2) 
			      * $T->element($i1,$j1)
			      * $T->element($i2,$j2);
		}
		}
		$g2->assign($i1, $i2, $val);
	}
	}
	my $NewCrystal = new Crystal();
	$NewCrystal->SetSpaceGroup("P 1", 1);
	$NewCrystal->SetMetrics(
		$g2->element(1,1), $g2->element(1,2), $g2->element(1,3),
		$g2->element(2,1), $g2->element(2,2), $g2->element(2,3),
		$g2->element(3,1), $g2->element(3,2), $g2->element(3,3));
	$NewCrystal->CalLatticeParametersFromMetrics();

	return $NewCrystal if($LatticeParameterOnly);

	$this->ExpandCoordinates();

	my $amin = 0;
	my $amax = 0;
	for(my $i = 1 ; $i <= 3 ; $i++) {
		my $e = $T->element($i,1);
		$amin = int($e-0.999) if($e < $amin);
		$amax = int($e-0.001) if($e > $amax);
	}
	my $bmin = 0;
	my $bmax = 0;
	for(my $i = 1 ; $i <= 3 ; $i++) {
		my $e = $T->element($i,2);
		$bmin = int($e-0.999) if($e < $bmin);
		$bmax = int($e-0.001) if($e > $bmax);
	}
	my $cmin = 0;
	my $cmax = 0;
	for(my $i = 1 ; $i <= 3 ; $i++) {
		my $e = $T->element($i,3);
		$cmin = int($e-0.999) if($e < $cmin);
		$cmax = int($e-0.001) if($e > $cmax);
	}
	if($CheckConsistency) {
		$amax++; $bmax++; $cmax++;
	}
#$bmin = 0;
print "Expansion range: ($amin,$bmin,$cmin) - ($amax,$bmax,$cmax)<br>\n" if($IsPrint);
	my @ExpandedAtomSiteList   = $this->GetCExpandedAtomSiteList();
	$nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
	my $count = 0;
	my %AtomTypeCount;

	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atom = $ExpandedAtomSiteList[$i];
		my $atomname = $atom->AtomNameOnly();
		my $charge   = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occ = $atom->Occupancy();

		for(my $ia = $amin ; $ia <= $amax ; $ia++) {
		for(my $ib = $bmin ; $ib <= $bmax ; $ib++) {
		for(my $ic = $cmin ; $ic <= $cmax ; $ic++) {
			my $xyz = new Math::MatrixReal(3, 1);
			$xyz->assign(1,1, $x+$ia);
			$xyz->assign(2,1, $y+$ib);
			$xyz->assign(3,1, $z+$ic);
			$xyz = $tRT * $xyz;
			my $x1 = Utils::Reduce01($xyz->element(1,1));
			my $y1 = Utils::Reduce01($xyz->element(2,1));
			my $z1 = Utils::Reduce01($xyz->element(3,1));
#print "(",$x+$ia,",",$y+$ib,",",$z+$ic,",) - ($x1,$y1,$z1)<br>";
			my $WillAdd = 1;
			my $MinimumDistance = 0.1;
			if($NewCrystal->nAsymmetricAtomSite() == 0) {
				$AtomTypeCount{$atomname}++;
				my $label = $atomname . $AtomTypeCount{$atomname};
				$NewCrystal->AddAtomSite($label, $atomname, $x1, $y1, $z1, $occ);
#print "Passed0 ($x1,$y1,$z1)<br>";
				next;
			}
			my @AtomList = $NewCrystal->GetCAsymmetricAtomSiteList();
			my $nAtom = @AtomList;
			for(my $j = 0 ; $j < @AtomList ; $j++) {
				my $atom0         = $AtomList[$j];
#				my $atom0         = $NewCrystal->GetCAtomSite($j);
				my $atomname0     = $atom0->AtomNameOnly();
				my ($x0,$y0,$z0)  = $atom0->Position(1);
				next if($atomname ne $atomname0);

				my $dis = $NewCrystal->GetNearestInterAtomicDistance(
					$x0, $y0, $z0, $x1, $y1, $z1);
#print "$nAtom: ($i,$j): ($ia,$ib,$ic): ($x0,$y0,$z0)-($x1,$y1,$z1) ";
#print "dis: $dis<br>";
				if($dis < $MinimumDistance) {
#print "Dropped<br>";
					$WillAdd = 0;
					last;
				}
			}
			if($WillAdd) {
#print "Passed ($x1,$y1,$z1)<br>";
				$AtomTypeCount{$atomname}++;
				my $label = $atomname . $AtomTypeCount{$atomname};
				$NewCrystal->AddAtomSite($label, $atomname, $x1, $y1, $z1, $occ);
			}
		}
		}
		}
		$count++;
	}
	$NewCrystal->ExpandCoordinates();

	return $NewCrystal;
}

sub SetMetrics
{
	my ($this,$g11,$g12,$g13,$g21,$g22,$g23,$g31,$g32,$g33)
		= @_;
	$this->{'g11'} = $g11;
	$this->{'g12'} = $g12;
	$this->{'g13'} = $g13;
	$this->{'g21'} = $g21;
	$this->{'g22'} = $g22;
	$this->{'g23'} = $g23;
	$this->{'g31'} = $g31;
	$this->{'g32'} = $g32;
	$this->{'g33'} = $g33;
	return ($this,$g11,$g12,$g13,$g21,$g22,$g23,$g31,$g32,$g33);
}

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

	my $a = sqrt($this->{'g11'});
	my $b = sqrt($this->{'g22'});
	my $c = sqrt($this->{'g33'});
	my $alpha = $this->{'g23'} / $b / $c;
	$alpha = Sci::dacos($alpha);
	my $beta  = $this->{'g13'} / $a / $c;
	$beta  = Sci::dacos($beta);
	my $gamma = $this->{'g12'} / $a / $b;
	$gamma = Sci::dacos($gamma);
	$a = Sci::Round($a, 6);
	$b = Sci::Round($b, 6);
	$c = Sci::Round($c, 6);
	$alpha = Sci::Round($alpha, 4);
	$beta  = Sci::Round($beta , 4);
	$gamma = Sci::Round($gamma, 4);
	$this->SetLatticeParameters($a,$b,$c,$alpha,$beta,$gamma);
}

sub SetIsSelected
{
	my($this, $index, $IsSelected) = @_;
	my $RefAtomList = $this->{"RefExpandedAtomSiteList"};
	my $atom = $RefAtomList->[$index];
	if($atom ne '') {
		return $atom->SetIsSelected($IsSelected);
	}
	return 0;
}

sub IsSelected
{
	my($this, $index) = @_;
	my $RefAtomList = $this->{"RefExpandedAtomSiteList"};
	my $atom = $RefAtomList->[$index];
	if($atom ne '') {
		return $atom->IsSelected();
	}
	return 0;
}

sub SameIdSiteIndex
{
	my ($this, $idm1, $sid) = @_;
	my $SameIdSiteIndex = $this->{'SameIdSiteIndex'};
	return $SameIdSiteIndex->{"[$idm1][$sid]"};
}

sub SortAtomTypeOrder
{
	my ($this, $sorttype) = @_;
#print "<H1>SORT: $sorttype</H1>";
	if($sorttype =~ /ChargeAscend/i) {
		my $ref = $this->{"RefAtomTypeList"};
		@$ref = sort { $a->Charge() <=> $b->Charge() } @$ref;
	}
	elsif($sorttype =~ /ChargeDescend/i) {
		my $ref = $this->{"RefAtomTypeList"};
		@$ref = sort { -$a->Charge() <=> -$b->Charge() } @$ref;
	}
	return 1;
}

sub FillAtomTypeData
{
	my ($this) = @_;
	my $ref = $this->{"RefAtomTypeList"};
	for(my $i = 0 ; $i < @$ref ; $i++) {
		$ref->[$i]->FillAtomTypeData();
	}
}


# $s = sin Q / lambda (in nm)
sub AtomicScatteringFactor
{
	my ($this, $s) = @_;
	return $this->asf($s);
}

# $s = sin Q / lambda (in nm)
sub asfElectron
{
	my ($this, $iAtomType, $s) = @_;
	$s = 1.0e-5 if($s <= 0.0);
	my $s2 = $s * $s;
	my $k  = 1.0e-2 * 8.0 * $pi * $pi / ($a0 * 1.0e9);
	my $Z  = $this->{"ASFDCAtomNameZ$iAtomType"};

	my $asfXray = $this->asf($iAtomType, $s);
#print "s=$s, i=$iAtomType, Z=$Z, a=$asfXray\n";
#exit;
	my $val = $k * ($Z - $asfXray) / $s2;
#print "va [$iAtomType]=$val = $k * ($Z - $asfXray) / $s2\n";
	return $val;
}

sub ReadASF
{
	my ($this, $Source) = @_;

	my @AtomTypeList = $this->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $Rietan = new Rietan;
		my $AtomName;
		my $AtomNameOnly = $AtomTypeList[$i]->AtomNameOnly(1);
		my $Charge   = $AtomTypeList[$i]->Charge();
		if($Charge == 1) {
			$AtomName = "$AtomNameOnly+";
		}
		elsif($Charge == -1) {
			$AtomName = "$AtomNameOnly-";
		}
		elsif($Charge > 1) {
			$AtomName = "$AtomNameOnly+";
		}
		elsif($Charge < -1) {
			$Charge = -$Charge;
			$AtomName = "$AtomNameOnly-";
		}
#print "a1: $AtomName\n";
		my @a = $Rietan->ReadASFParameters($AtomName);
		if(@a < 9) {
			$AtomName = $AtomTypeList[$i]->AtomNameOnly(1);
#print "a2: $AtomName\n";
			@a = $Rietan->ReadASFParameters($AtomName);
		}
		if(@a < 9) {
			print "Error: Can not read ASFDC for [$AtomName].\n";
			next;
		}
#print "array: ", join(',', @a), "\n";
#print "i=$i: $AtomName\n";
		$this->{"ASFDCRed$i"}              = $Rietan;
		$this->{"ASFDCRed$AtomName"}       = $Rietan;
		$this->{"ASFDCAtomName$i"}         = $Rietan->{ASFDCAtomName};
		$this->{"ASFDCAtomName$AtomName"}  = $Rietan->{ASFDCAtomName};
		my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
		$this->{"ASFDCAtomNameZ$i"}            = $AtomicNumber;
		$this->{"ASFDCAtomNameZ$AtomNameOnly"} = $AtomicNumber;
print "ASFDC AtomName [$i]: $Rietan->{ASFDCAtomName} [Zi=$AtomicNumber]\n";
#print "   ao:$AtomNameOnly: $i\n";
	}
	$this->{ASFDCRed} = 1;
}

# $s = sin Q / lambda (in nm)
sub asf
{
	my ($this, $iAtomType, $s, $Source) = @_;
	$Source = 'Cu' if(!defined $Source);
	$this->ReadASF($Source) 	if(!defined $this->{ASFDCRed});

	my $Rietan = $this->{"ASFDCRed$iAtomType"};
	my $pP     = $Rietan->{pASFDCParameters};
	my ($A1,$B1,$A2,$B2,$A3,$B3,$A4,$B4,$C, $b,
		$Crdf1,$Crdf2, $Fedf1,$Fedf2, $Codf1,$Codf2, $Cudf1,$Cudf2, 
		$Modf1, $Modf2, $Agdf1,$Agdf2, $Kbdf1, $Kbdf2) = @$pP;
#	if($Source eq 'n') {
#		return $b;
#	}

	my $s2 = 0.01 * $s * $s;
	my $asf = $A1*exp(-$B1*$s2)+$A2*exp(-$B2*$s2)+$A3*exp(-$B3*$s2)+$A4*exp(-$B4*$s2)+$C;
	if($Source eq 'Cr' and defined $Crdf2) {
		$asf = cplx($asf + $Crdf1, $Crdf2);
	}
	elsif($Source eq 'Fe' and defined $Fedf2) {
		$asf = cplx($asf + $Fedf1, $Fedf2);
	}
	elsif($Source eq 'Co' and defined $Codf2) {
		$asf = cplx($asf + $Codf1, $Codf2);
	}
	elsif($Source eq 'Cu' and defined $Cudf2) {
#print "d=$Cudf1,$Cudf2\n";
		$asf = cplx($asf + $Cudf1, $Cudf2);
	}
	elsif($Source eq 'Mo' and defined $Modf2) {
		$asf = cplx($asf + $Modf1, $Modf2);
	}
	elsif($Source eq 'Ag' and defined $Agdf2) {
		$asf = cplx($asf + $Agdf1, $Agdf2);
	}
	elsif($Source eq 'Kb' and defined $Kbdf2) {
		$asf = cplx($asf + $Kbdf1, $Kbdf2);
	}

	return $asf;
}

sub Fhkl
{
	my ($this, $h, $k, $l, $IsElectron) = @_;
	$IsElectron = 0 if(!defined $IsElectron);

	my $F  = cplx(0, 0);
	my @ExpandedAtomSiteList   = $this->GetCExpandedAtomSiteList();
	if(@ExpandedAtomSiteList == 0) {
		$this->ExpandCoordinates();
		@ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
	}

	my $dhkl = $this->CalculateHKLInterplanarSpacing($h, $k, $l);
	my $s    = 1.0 / 2.0 / ($dhkl * 0.1); # in nm
#print "d=$dhkl  s=$s\n";
	my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $AtomName  = $ExpandedAtomSiteList[$i]->AtomNameOnly(1);
		my $iAtomType = $ExpandedAtomSiteList[$i]->iAtomType() - 1;
		my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
		my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
		my $asf;
		if($IsElectron) {
			$asf = $occupancy * $this->asfElectron($iAtomType, $s);
		}
		else {
			$asf = $occupancy * $this->asf($iAtomType, $s);
		}
#print "i=$i [$h,$k,$l]: pos($x,$y,$z) o=$occupancy  asf=$asf ($AtomName:", $this->{"ASFDCAtomName$iAtomType"}, ")\n";
		my $phi = 2.0 * $pi * ($h * $x + $k * $y + $l * $z);
		$F += $asf * exp(i * $phi);
	}
	return (Re($F), Im($F));
#	return ($Fr, $Fi);
}

sub OrthorhombicMultiplicity
{
	my ($this, $m, $h, $k, $l) = @_;
	$m *= 2 if($h != 0);
	$m *= 2 if($k != 0);
	$m *= 2 if($l != 0);
	return $m;
}

sub TetragonalMultiplicity
{
	my ($this, $m, $h, $k, $l) = @_;
	$h = abs($h);
	$k = abs($k);
	$m *= 2 if($h != 0);
	$m *= 2 if($k != 0);
	$m *= 2 if($l != 0);
	if($h != $k) {
		$m *= 2;
	}
	return $m;
}

sub HexagonalMultiplicity
{
	my ($this, $m, $h, $k, $l) = @_;
	return 16 if($h != 0 and $k != 0 and $l != 0);
	return  8 if($h == $k and $l != 0);
	return 12 if($h == 0 and $k != 0 and $l != 0);
	return 12 if($h != 0 and $k == 0 and $l != 0);
	return 12 if($h != 0 and $k != 0 and $l == 0);
	return  6 if($h == $k and $l == 0);
	return  6 if($h == 0 and $k != 0 and $l == 0);
	return  6 if($h != 0 and $k == 0 and $l == 0);
	return  2 if($h == 0 and $k == 0 and $l != 0);
print "Error in CrystalObject::HexagonalMultiplicity: Invalid hkl [$h,$k,$l]\n";
	return undef;
}

sub TrigonalMultiplicity
{
	my ($this, $m, $h, $k, $l) = @_;
	return $this->HexagonalMultiplicity($m, $h, $k, $l);
}

sub RhombohedralMultiplicity
{
	my ($this, $m, $h, $k, $l) = @_;
	return $this->HexagonalMultiplicity($m, $h, $k, $l);
}

sub CubicMultiplicity
{
	my ($this, $m, $h, $k, $l) = @_;
	$h = abs($h);
	$k = abs($k);
	$l = abs($l);
	$m *= 2 if($h != 0);
	$m *= 2 if($k != 0);
	$m *= 2 if($l != 0);
	if($h != $k and $k != $l and $l != $h) {
		$m *= 6;
	}
	elsif($h == $k and $k != $l) {
		$m *= 3;
	}
	elsif($k == $l and $l != $h) {
		$m *= 3;
	}
	elsif($l == $h and $h != $k) {
		$m *= 3;
	}
	return $m;
}

sub Multiplicity
{
	my ($this, $h, $k, $l) = @_;
	my $ls = $this->LatticeSystem();
#print "ls: $ls\n";
	my $m = 1;
	if($ls eq 'cubic') {
		return $this->CubicMultiplicity($m, $h, $k, $l);
	}
	elsif($ls eq 'tetragonal') {
		return $this->TetragonalMultiplicity($m, $h, $k, $l);
	}
	elsif($ls eq 'trigonal') {
		return $this->TrigonalMultiplicity($m, $h, $k, $l);
	}
	elsif($ls eq 'hexagonal') {
		return $this->HexagonalMultiplicity($m, $h, $k, $l);
	}
	elsif($ls eq 'orthorhombic') {
		return $this->OrthorhombicMultiplicity($m, $h, $k, $l);
	}
	elsif($ls eq 'triclinic') {
		return 2;
	}
	else {
		print "Error in CrystalObject::Multiplicity: Lattice system [$ls] is not supported.\n";
		return undef;
	}
}

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 '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);
}
1;