package Crystal;
use Crystal::CrystalObject;
@ISA = qw(CrystalObject);

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

use strict;

BEGIN {
}

sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub ExpandToSuperLattice
{
	my ($this,$nx,$ny,$nz,$crystal) = @_;

my $LF="<br>\n";
print "<H3>Expand to SuperLattice in Crystal.pm.</H3>\n";

	$nx = int($nx + 0.01);
	$ny = int($ny + 0.01);
	$nz = int($nz + 0.01);
	$nx = 1 if($nx < 1);
	$ny = 1 if($ny < 1);
	$nz = 1 if($nz < 1);
	my $nMult = $nx*$ny*$nz;

	$this->SetCrystalName($crystal->CrystalName());
	$this->SetCAtomTypeList($crystal->GetCAtomTypeList());
	$this->SetnAtomType($crystal->nAtomType());
	$this->SetCAsymmetricAtomSiteList($crystal->GetCAsymmetricAtomSiteList());
	$this->SetnAsymmetricAtomSite($crystal->nAsymmetricAtomSite());
	$this->SetDensity($crystal->Density());
	my ($a,$b,$c,$alpha,$beta,$gamma) = $crystal->LatticeParameters();
	$this->SetLatticeParameters($a*$nx, $b*$ny, $c*$nz, $alpha, $beta, $gamma);

	my @ExpandedAtomSiteList   = $crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $crystal->nTotalExpandedAtomSite();

#超格子へ展開する
	my @AtomSiteList;
	my $count = 0;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $label     = $ExpandedAtomSiteList[$i]->Label();
		my $type      = $ExpandedAtomSiteList[$i]->AtomName();
		my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
		my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
		my $id        = $ExpandedAtomSiteList[$i]->IdAsymmetricAtomSite();
		
		for(my $ix = 0 ; $ix < $nx ; $ix++) {
		for(my $iy = 0 ; $iy < $ny ; $iy++) {
		for(my $iz = 0 ; $iz < $nz ; $iz++) {
			my $x0 = ($x + $ix) / $nx;
			my $y0 = ($y + $iy) / $ny;
			my $z0 = ($z + $iz) / $nz;

			$AtomSiteList[$count] = new AtomSite();
			$AtomSiteList[$count]->SetLabel($label);
			$AtomSiteList[$count]->SetAtomName($type);
			$AtomSiteList[$count]->SetPosition($x0, $y0, $z0);
			$AtomSiteList[$count]->SetOccupancy($occupancy);
			$AtomSiteList[$count]->SetIdAsymmetricAtomSite($id);

			$count++;
		}
		}
		}
	}
#print "count: $count\n";
	$this->{"nTotalExpandedAtomSite"}  = $count;
	$this->{"RefExpandedAtomSiteList"} = \@AtomSiteList;

#多重度を超格子の値に更新
	my @nMultiplicityExpandedAtomSiteList
		 = $crystal->GetCnMultiplicityExpandedAtomSiteList();
	for(my $i = 0 ; $i < @nMultiplicityExpandedAtomSiteList ; $i++) {
		$nMultiplicityExpandedAtomSiteList[$i] *= $nMult;
	}
	$this->{"RefnExpandedAtomSite"} = \@nMultiplicityExpandedAtomSiteList;

#多重度を、全原子に設定
	for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
		my $atom = $AtomSiteList[$i];
		my $ia   = $atom->Id() - 1;
		$atom->SetMultiplicity($nMultiplicityExpandedAtomSiteList[$ia]);
		$this->SetIsSelected($i, 0);
	}

#IDをもとに部分占有原子を整理
#Occupancyに合わせて表示する原子数を選択する => IsSelectedに選択された原子を設定する。

#同じ位置にあるサイトの数
	my @nIdenticalSite;
#そのサイトを実際に占める原子の数
	my %nOccupiedAtom;
	my %iIdenticalPartialAtomWithSameId;
#ID=$idのサイトで、$j番目のサイト占める原子の名前　$AtomTypeInPartialAtom{"[$id][$j]"}
	my %AtomTypeInPartialAtom;
#ID=$idのサイトで、$j番目のサイトを占める原子のID　$AtomIdInPartialAtom{"[$id][$j]"}
	my %AtomIdInPartialAtom;
#ID=$idのサイトで、$j番目のサイトを占める原子のIndex　$AtomIndexInPartialAtom{"[$id][$j]"}
	my %AtomIndexInPartialAtom;
#多重度
	my @nSameIdSite;
	my %SameIdSiteIndex;


#同じIDをもつサイトのindexを$SameIdSiteIndex{"[$idm1][$sid]"} に設定する
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++){
		my $atom      = $ExpandedAtomSiteList[$i];
		my $type      = $atom->AtomName();
		my ($x,$y,$z) = $atom->Position(1);
		my $occ       = $atom->Occupancy();
		my $mult      = $atom->Multiplicity();
		my $id        = $atom->Id();
		my $idm1      = $id-1;
		$nSameIdSite[$idm1]++;
		my $sid = $nSameIdSite[$idm1]-1;
		$SameIdSiteIndex{"[$idm1][$sid]"} = $i;
#printf "#++%4d[%4s]: (%7.4f, %7.4f, %7.4f) [%6.4f] (mult=%3d: id=%2d)$LF",
#	$i+1, $type, $x, $y, $z, $occ, $mult, $id;
	}

print "<b>#Check number of atoms occupying the sites.</b>$LF";
	my @IsPrintedById;
	my @IsPrinted;
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $PrintCount = 0;
		my $atom = $ExpandedAtomSiteList[$i];
		my $type = $atom->AtomName();
		my ($x,$y,$z) = $atom->Position(1);
		my $occ  = $atom->Occupancy();
		my $id   = $atom->Id();
		my $mult = $atom->Multiplicity();
		my $idm1 = $id-1;
#該当する等価な全てサイトのうち、何個の原子が存在するか
		$nOccupiedAtom{"[$i][$id]"} = int($mult * $occ + 0.01);
		my $nc = $nOccupiedAtom{"[$i][$id]"};
		next if($IsPrintedById[$id]);
		next if($occ >= 0.999999);

#print "+++A(i:$i): $type: occ=$occ: nMult=$mult: nAtom=$nc (id=$id)$LF";

#同じ位置にある原子のIDのリストを作る
		$iIdenticalPartialAtomWithSameId{"[$i][0]"} = $i;

		$IsPrinted[$i] = 1;
		$IsPrintedById[$id] = 1;

#同じ位置にあるサイトの数を$cにいれる
		my $c = 1;
		for(my $j = $i+1 ; $j < @ExpandedAtomSiteList ; $j++) {
			next if($IsPrinted[$j]);
			my $atom2 = $ExpandedAtomSiteList[$j];
			my $type2 = $atom2->AtomName();;
			my ($x2,$y2,$z2) = $atom2->Position(1);
			my $occ2  = $atom2->Occupancy();
			my $id2   = $atom2->Id();
			my $mult2 = $atom2->Multiplicity();

			next if($occ2 >= 0.999999);
			my $dis = $this->GetNearestInterAtomicDistance(
					$x, $y, $z, $x2, $y2, $z2);
			next if($dis > 0.001);

			$iIdenticalPartialAtomWithSameId{"[$i][$c]"} = $j;

#$i番目のサイトと同じ位置にあるID=$id2のサイト、実際にいくつの原子があるか
			$nOccupiedAtom{"[$i][$id2]"} = int($mult2 * $occ2 + 0.01);

			$IsPrinted[$j] = 1;
			$IsPrintedById[$id2] = 1;
			$c++;

my $nc2 = $nOccupiedAtom{"[$i][$id2]"};
#print "+++***(j:$j)$type2: occ=$occ2: nMult=$mult2: nAtom=$nc2 (id=$id2)$LF";
		}

#$i番と同じ位置にあるサイトサイトの数を$cにいれる
#$iは、0から始まる（ID-1に等しい）
		$nIdenticalSite[$i] = $c;
	}

#結果の表示
	for(my $id0 = 0 ; $id0 < @ExpandedAtomSiteList ; $id0++) {
#$nIdenticalSiteは部分共有原子の$iのみで定義される。
#同じ位置にある原子の値は0なので、この条件文でスキップできる
		next if(!defined $nIdenticalSite[$id0]or $nIdenticalSite[$id0] <= 0);

		my $i  = $iIdenticalPartialAtomWithSameId{"[$id0][0]"};
		my $atom = $ExpandedAtomSiteList[$i];
		my ($x0,$y0,$z0) = $atom->Position(1);
		my $Id           = $atom->Id($i);
		my $occupancy    = $atom->Occupancy($i);
		my $mult         = $atom->Multiplicity($i);
		my $Idm1 = $Id-1;

#$nIdenticalSiteは0から始まる配列
print "<b>ID=</b>$Id: nIdenticalSite = ", $nIdenticalSite[$id0], "  at ($x0, $y0, $z0)$LF";
print "+++ i=";
		for(my $sn = 0 ; $sn < $nIdenticalSite[$id0] ; $sn++) {
			my $j = $iIdenticalPartialAtomWithSameId{"[$id0][$sn]"};
print "$j ";
		}
print "$LF";

		for(my $sn = 0 ; $sn < $nIdenticalSite[$id0] ; $sn++) {
			my $j = $iIdenticalPartialAtomWithSameId{"[$i][$sn]"};
			next if(!defined($j) or $j eq '');
			my $atom  = $ExpandedAtomSiteList[$j];
			my $Id2          = $atom->Id();
			my $Id2m1        = $Id2-1;
			my $type2        = $atom->AtomName();
			my ($x2,$y2,$z2) = $atom->Position(1);
			my $occ2         = $atom->Occupancy();
			my $mult2        = $atom->Multiplicity();
			my $nAtom2 = $nOccupiedAtom{"[$i][$Id2]"};

print "+++ i=$j(Id=$Id2): $type2 ($x2, $y2, $z2)"
	. " [occ=$occ2] [nMult=$mult2] [N=$nAtom2]$LF";
print "+++*** SameIdSiteIndex";
			for(my $j = 0 ; $j < $mult ; $j++) {
#				my $idx = $this->SameIdSiteIndex($Id2m1, $j);
				my $idx = $SameIdSiteIndex{"[$Id2m1][$j]"};
printf "%3d ", $idx;
			}
print "$LF";
		}


#$SameIdSiteIndex{"[$Id2][$j]"}から、実際に割り当てる原子を決める
		my $nSiteWithSameId0 = $atom->Multiplicity();
		my $ResidualNum      = $nSiteWithSameId0;
		my $nIdenticalSite0  = $nIdenticalSite[$i];
#print "nSiteWithSameId0: $nSiteWithSameId0, nIdenticalSite0: $nIdenticalSite0$LF";

		for(my $sn = 0 ; $sn < $nIdenticalSite0 ; $sn++) {
			my $jatom = $iIdenticalPartialAtomWithSameId{"[$i][$sn]"};
#print "jatom=$jatom$LF";
			my $atom0 = $ExpandedAtomSiteList[$jatom];
			my $Id0   = $atom0->Id();
			my $type  = $atom0->AtomName();
			my $nAtom = $nOccupiedAtom{"[$i][$Id0]"};
#print "j=$jatom: Id=$Id0: $type: nOccupiedAtom=$nAtom$LF";
			for(my $is = 0 ; $is < $nAtom ; $is++) {
				my $r = int(rand($ResidualNum));
				for(my $j = 0 ; $j < $nSiteWithSameId0 ; $j++) {
					next if(defined($AtomTypeInPartialAtom{"[$Id][$j]"}));
					$r--;
					if($r < 0) {
#ID=$idの原子と同一のサイトで、$j番目のサイトをしめる原子 $AtomTypeInPartialAtom{"[$id][$j]"}
						$AtomTypeInPartialAtom{"[$Id][$j]"} = $type;
						$AtomIdInPartialAtom{"[$Id][$j]"}   = $Id0;
						my $Id0m1 = $Id0 - 1;
						my $Index = $SameIdSiteIndex{"[$Id0m1][$j]"}; 
#表示する原子を設定する

						$this->SetIsSelected($Index, 1);
						$AtomIndexInPartialAtom{"[$Id][$j]"} = $Index;
#print "Id0=$Id0, sn=$j, Idx=$Index$LF";
						$ResidualNum--;
						last;
					}
				}
			}
		}
	}

if(0) {
print "<b>Summary:</b>$LF";
#	my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom         = $ExpandedAtomSiteList[$i];
		my $name         = $atom->AtomName();
		my $occupancy    = $atom->Occupancy();
#表示する原子を設定する
		$this->SetIsSelected($i, 1) if($occupancy >= 1);
		next if(!$this->IsSelected($i));

		my ($x,$y,$z)    = $atom->Position(1);
		my $Id           = $atom->Id();
print "+++$i: $name (id=$Id): ($x, $y, $z) [$occupancy]$LF";
	}
}

print "<H3>Expand to Superlattice in Crystal.pm: End</H3>\n";
print "$LF";

	return $count;
}

sub SetOutputMode
{
	my ($this, $mode) = @_;
	return $this->{'OutputMode'} = $mode;
}

sub OutputMode
{
	my ($this) = @_;
	return $this->{'OutputMode'};
}

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

	if($this->{OutputMode} eq 'expanded' or
	   $this->{OutputMode} eq 'choose' ) {
		return "P 1";
	}
	return $this->SPGName();
}

sub iSPGByOutputMode
{
	my ($this) = @_;
	if($this->{'OutputMode'} eq 'expanded' or
	   $this->{'OutputMode'} eq 'choose' ) {
		return 1;
	}
	return $this->iSPG();
}

sub iSetByOutputMode
{
	my ($this) = @_;
	if($this->{'OutputMode'} eq 'expanded' or
	   $this->{'OutputMode'} eq 'choose' ) {
		return 1;
	}
	return $this->iSet();
}

sub LatticeParametersByOutputMode
{
	my ($this, $UseAtomicUnit) = @_;
	if($this->{'OutputMode'} eq 'expanded' or
	   $this->{'OutputMode'} eq 'choose' ) {
		my $SuperLattice = $this->{'CSuperLattice'};
		return $SuperLattice->LatticeParameters($UseAtomicUnit)
			if($SuperLattice);
		return $this->LatticeParameters($UseAtomicUnit);
	}
	return $this->LatticeParameters($UseAtomicUnit);
}

sub LatticeVectorsByOutputMode
{
	my ($this, $UseAtomicUnit) = @_;
	my $k = ($UseAtomicUnit)? 1.0 / 0.529177 : 1.0;
#my $SPG     = $this->GetCSpaceGroup();
#my $SPGName = $SPG->SPGName();
#print "SPG in Crystal.pm:[$SPGName]\n";
#exit;
	if($this->{'OutputMode'} eq 'expanded' or
	   $this->{'OutputMode'} eq 'choose' ) {
		my $SuperLattice = $this->{'CSuperLattice'};
		return $SuperLattice->LatticeVectors($UseAtomicUnit)
			if($SuperLattice);
		return $this->LatticeVectors($UseAtomicUnit);
	}
	return $this->LatticeVectors($UseAtomicUnit);
}

sub GetCExpandedAtomSiteListByOutputMode
{
	my ($this,$UseAtomicUnit) = @_;
	if(defined $this->{'OutputMode'} and $this->{'OutputMode'} eq 'expanded') {
		my $SuperLattice = $this->{'CSuperLattice'};
		return $SuperLattice->GetCExpandedAtomSiteList()
			if($SuperLattice);
		return $this->GetCExpandedAtomSiteList();
	}
	elsif(defined $this->{'OutputMode'} and $this->{'OutputMode'} eq 'choose' ) {
		my $SuperLattice = $this->{'CSuperLattice'};
		if($SuperLattice) {
			my @slist = $SuperLattice->GetCExpandedAtomSiteList();
			my @tlist;
			for(my $i = 0 ; $i < @slist ; $i++) {
				my $atom = $slist[$i];
				next if($atom->Occupancy() < 1.0 and
					!$atom->IsSelected() );
				push(@tlist, $atom);
			}
			return @tlist;
		}
		return $this->GetCExpandedAtomSiteList();
	}
	return $this->GetCAsymmetricAtomSiteList();
}

sub CreateSuperLattice
{
	my ($this, $nx, $ny, $nz) = @_;

	my $SuperLattice = new Crystal();
	$SuperLattice->ExpandToSuperLattice($nx,$ny,$nz,$this);
	$this->{'CSuperLattice'} = $SuperLattice;
	return 1;
}

sub GetCSuperLattice
{
	my ($this) = @_;
	return $this->{'CSuperLattice'};
}
1;