#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';
#use lib 'd:/MyWebs/cgi-bin/lib';

use lib '/home/tkamiya/bin/lib';
use lib '/home/tkamiya/bin';

use MyApplication;
use Crystal::CIF;
use Crystal::Crystal;

use strict;
#use warnings;

#===============================================
# Applicationクラス作製
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

my $LF = $App->LF();

# 実行プログラム名（デフォルトでは$0から取得）
	my $Program = $App->Program();
# アプリケーション名
	$App->SetAppName($Program);
# バージョン
	my $Version = $App->SetVersion("Ver 0.1");

#===============================================
# 環境設定
#===============================================
my $ProgramPath = $App->SpeculateProgramPath($0);
print "Program Path: $ProgramPath\n";

#最後の引数を0にすると、IniFileが存在しなくてもエラーがでない
my $IniFile = $App->OpenIniFile($ProgramPath, 0);

#===============================================
# 変数
#===============================================
my $CoordRMax = 2.7;
my $nMaxStep  = 6;

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--superlattice",     "--superlattice=[0|1]: Show this help");
$App->AddArgument("--help",     "--help    : Show this help");
exit 1 if($App->ReadArgs() != 1);
# Iniファイル読み込み
&ConfigureIniFileVariables();

my $CIFFile        = $App->GetGetArg(0);
$CIFFile = 'igzo.cif' if($CIFFile eq '');
my $DoSuperlattice = $App->GetGetArg('superlattice');

#CIFクラスを作成
print "key: $CIFFile\n";
unless($CIFFile) {
	$App->Usage();
	exit 1;
}
$App->{'CIF'} = new CIF;
#CIFクラスを取得
my $cif = $App->{'CIF'};
#CIFファイル名を引数に与えて読み込む
unless($cif->Read($CIFFile)) {
#unless($cif->Read($CIFFile)) {
	print "Error: Can not read [$CIFFile]\n\n";
	exit 2;
}
print "CIF: ", $cif->FileName(), "$LF";

my $Crystal = &ReadStructure($cif);
my $pSiteList = $Crystal->FindCoordinationsForAllAtoms(0, $CoordRMax, 0);
my $nSite     = @$pSiteList;

print "\n\n";
print "Find coordination structures...\n";
for(my $i = 0 ; $i < $nSite ; $i++) {
	my $pAtom0       = $pSiteList->[$i];
	my $label0       = $pAtom0->Label();
	my $type0        = $pAtom0->AtomName();
	my ($x0,$y0,$z0) = $pAtom0->Position(1);
	my $idSite0      = $pAtom0->{idSite};
	my $pA           = $pAtom0->{pCoordAtoms};

print "Coordination arround $i-th atom [$label0: $type0 (id=$idSite0)] ($x0, $y0, $z0)\n";
	for(my $i = 0 ; $i < @$pA ; $i++) {
		my $pCA            = $pA->[$i];
		my $iSite1         = $pCA->{iSite};
		my $idSite1        = $pCA->{idSite};
		my $label1         = $pCA->Label();
		my $type1          = $pCA->AtomName();
		my ($x1, $y1, $z1) = $pCA->Position();
		my ($X1, $Y1, $Z1) = @{$pCA->{pCoord}};
		my $r              = $pCA->{r};
printf "  %03d: [%3s: %4s (id = %03d)] (%8.4f, %8.4f, %8.4f): r = %8.4f  frac:(%6.3f, %6.3f, %6.3f)\n",
		$iSite1, $label1, $type1, $idSite1, $X1, $Y1, $Z1, $r, $x1, $y1, $z1;
	}
}

print "\n\n";
print "Find rings...\n";
my @AtomSiteList           = $Crystal->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite    = @AtomSiteList;
#my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
#my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();

my @CheckedId;
for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
	next if($CheckedId[$i]);
	$CheckedId[$i] = 1;
	my @link;
	my $ret = &FindRingsForAnAtom($Crystal, $pSiteList, \@CheckedId, $i, $CoordRMax, $nMaxStep+1, 0, \@link);
exit;
}

exit 0;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# サブルーチン
#==========================================
sub FindRingsForAnAtom
{
	my ($Crystal, $pSiteList, $pCheckedId, $iSite, $CoordRMax, $nMaxStep, $iStep, $plink, $xi, $yi, $zi) = @_;
#ステップ数 (リング構成原子数) を超えた
	return undef if($iStep > $nMaxStep);

	my $pAtom0         = $pSiteList->[$iSite];
	my $label0         = $pAtom0->Label();
#	my $type0          = $pAtom0->AtomName();
	my $idSite0        = $pAtom0->{idSite};
	my ($x0, $y0, $z0) = $pAtom0->Position(0);
	if(!defined $xi) {
		($xi, $yi, $zi) = ($x0, $y0, $z0);
	}
	my $pCoord         = $pAtom0->{pCoordAtoms};

	$plink->[$iStep] = {
		pAtom     => $pAtom0,
		Label     => $label0,
		Name      => $pAtom0->AtomName(),
		pPosition => [$x0, $y0, $z0],
		};

# 配位子リストのあるAtomSiteをさがす
	my $pAtomSite      = $pSiteList->[$idSite0];
# Base座標: 配位子リストのあるAtomSiteの座標
#           ($x0, $y0, $z0)とは違う場合がある
	my ($xBase, $yBase, $zBase) = $pAtomSite->Position(1);
#print "Base ($xBase, $yBase, $zBase) for $label0 iSite=$iSite, id=$idSite0\n";
print "iStep=$iStep (< $nMaxStep): Center atom $label0 (iSite=$iSite) ($x0, $y0, $z0)\n";

	for(my $i = 0 ; $i < @$pCoord ; $i++) {
		my $pCA            = $pCoord->[$i];
		my $label1         = $pCA->Label();
		my $iSite1         = $pCA->{iSite};
		my $idSite1        = $pCA->{idSite};
# 配位子の座標はBase座標 ($xBase, $yBase, $zBase) を基準
		my ($x1, $y1, $z1) = $pCA->Position(0);
		my ($X1, $Y1, $Z1) = @{$pCA->{pCoord}};
#print "($x1, $y1, $z1)-($X1, $Y1, $Z1)  Base ($xBase, $yBase, $zBase)\n";
print "    Center $label0($iSite): i=$i Coordinated atom $label1 ($x1, $y1, $z1): idSite1=$idSite1  CheckedId=$pCheckedId->[$idSite1]\n";

# うまく機能しないので、とりあえずはずす
#		next if($pCheckedId->[$idSite1]);
		$pCheckedId->[$idSite1] = 1;

#		my $r              = $pCA->{r};
# 配位子の座標は、中心原子のBase座標 ($xBase, $yBase, $zBase) を基準
# 中心原子の座標 ($x0, $y0, $z0) とBase座標 ($xBase, $yBase, $zBase) の平行移動分を考慮して
# 初期原子座標からの距離を計算
		my $X = $X1 - ($xBase - $x0);
		my $Y = $Y1 - ($yBase - $y0);
		my $Z = $Z1 - ($zBase - $z0);
#print "X = $X1 - ($xBase - $x0) - $xi, Y = $Y1 - ($yBase - $y0) - $yi, Z = $Z1 - ($zBase - $z0) - $zi\n";
		my $r = $Crystal->GetInterAtomicDistance($xi, $yi, $zi, $X, $Y, $Z);

#ステップ数は3以上でないとリングにならない
		if($iStep > 2 and $r < 1.0e-3) {
print "*******Returned    : Level $iStep (<= $nMaxStep): $label1 ($X1, $Y1, $Z1)=>($X, $Y, $Z): r=$r (<$CoordRMax)\n";
			$plink->[$iStep+1] = {
				pAtom     => $pCA,
				Label     => $label1,
				Name      => $pCA->AtomName(),
				pPosition => [$X1, $Y1, $Z1],
				};
			for(my $l = 0 ; $l <= $iStep+1 ; $l++) {
				my $p      = $plink->[$l];
				my $pAtoml = $p->{pAtom};
				my ($xl, $yl, $zl) = $pAtoml->Position();
				print "   $l: $p->{Label}:$p->{Name} ($xl, $yl, $zl)\n";
			}
print "\n";
			return ($iStep+1, $plink);
		}
#残りステップでは原点に戻れない
		elsif($r > $CoordRMax * ($nMaxStep - $iStep)) {
print "       Not returned, terminated: Level $iStep (<= $nMaxStep): $label1 ($X1, $Y1, $Z1)=>($X, $Y, $Z): r=$r (<$CoordRMax)\n";
			next;
		}
		else {
print "       Not returned yet: Level $iStep (<= $nMaxStep): $label1 ($X, $Y, $Z)=>($X, $Y, $Z): r=$r (<$CoordRMax)\n";
		}

		&FindRingsForAnAtom($Crystal, $pSiteList, $pCheckedId, $iSite1, $CoordRMax, $nMaxStep, $iStep+1, $plink, $xi, $yi, $zi);
#print ">>Push any key\n";
#<>;
	}
}

sub ConfigureIniFileVariables
{
	$App->AddIniFileVariable("\\Preferences\\WorkDir",        "WorkDir");
	$App->ReadSetting();
}

sub ReadStructure
{
	my ($cif) = @_;

	my $Crystal = $cif->GetCCrystal();
$Crystal->SplitPartialSites();

	print "$LF";
	my $CrystalName = $Crystal->CrystalName();
	print "CrystalName: $CrystalName\n";
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	printf "cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a, $b, $c, $alpha, $beta, $gamma;
	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
	printf "  Vectors: (%9.6f %9.6f %9.6f)\n", $a11, $a12, $a13;
	printf "           (%9.6f %9.6f %9.6f)\n", $a21, $a22, $a23;
	printf "           (%9.6f %9.6f %9.6f)\n", $a31, $a32, $a33;
	my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();
	printf "  Metrics: |%12.6f %12.6f %12.6f|\n", $g11, $g12, $g13;
	printf "           |%12.6f %12.6f %12.6f|\n", $g21, $g22, $g23;
	printf "           |%12.6f %12.6f %12.6f|\n", $g31, $g32, $g33;
	my $vol = $Crystal->Volume();
	print "  Volume: $vol A^3$LF";
#対称要素を使い、非対称単位中の原子を展開する
#この際、Densityも計算される
	$Crystal->ExpandCoordinates();
	my $Density = $Crystal->Density();
	print "  Density: $Density g/cm^3$LF";

	print "$LF";
	my ($Ra,$Rb,$Rc,$Ralpha,$Rbeta,$Rgamma) = $Crystal->ReciprocalLatticeParameters();
	printf "Reciprocal cell (w/o 2pi): %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $Ra, $Rb, $Rc, $Ralpha, $Rbeta, $Rgamma;
	my ($Ra11, $Ra12, $Ra13, $Ra21, $Ra22, $Ra23, $Ra31, $Ra32, $Ra33) = $Crystal->ReciprocalLatticeVectors();
	printf "  Vectors: (%9.6f %9.6f %9.6f)\n", $Ra11, $Ra12, $Ra13;
	printf "           (%9.6f %9.6f %9.6f)\n", $Ra21, $Ra22, $Ra23;
	printf "           (%9.6f %9.6f %9.6f)\n", $Ra31, $Ra32, $Ra33;
	my ($Rg11, $Rg12, $Rg13, $Rg21, $Rg22, $Rg23, $Rg31, $Rg32, $Rg33) = $Crystal->RMetrics();
	printf "  Metrics: |%9.6f %9.6f %9.6f|\n", $Rg11, $Rg12, $Rg13;
	printf "           |%9.6f %9.6f %9.6f|\n", $Rg21, $Rg22, $Rg23;
	printf "           |%9.6f %9.6f %9.6f|\n", $Rg31 ,$Rg32, $Rg33;
	my $Rvol = $Crystal->CalculateReciprocalVolume();
	print "  Volume: $Rvol A^-3$LF";

	print "$LF";
#Crystalクラスから、SpaceGroupクラスを取り出す
#この時点で、並進対称要素の抽出、LatticeSystemの抽出は完了している
	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	my $iSPG    = $SPG->iSPG();
	my $LatticeSystem = $SPG->LatticeSystem();
	print "Space Group: $SPGName ($iSPG)$LF";
	print "Lattice System: $LatticeSystem$LF";

	print "$LF";
	my $nTranslation = $SPG->nTranslation();
	print "nTranslation: $nTranslation$LF";
	for(my $i = 0 ; $i < $nTranslation ; $i++) {
		my ($x,$y,$z) = $SPG->TranslationVector($i+1);
		print "   ($x, $y, $z)$LF";
	}

	print "$LF";
	my $nSymmetryOperation = $SPG->nSymmetryOperation();
	print "nSymmetryOperation: $nSymmetryOperation$LF";
	for(my $i = 0 ; $i < $nSymmetryOperation ; $i++) {
		my $symop = $SPG->SymmetryOperation($i+1);
		print "   $symop$LF";
	
		my ($x1,$y1,$z1,$t1,
			$x2,$y2,$z2,$t2,
			$x3,$y3,$z3,$t3)
			= $SPG->SymmetryOperationByMatrix($i+1);
		print "      $x1 $y1 $z1  $t1$LF";
		print "      $x2 $y2 $z2  $t2$LF";
		print "      $x3 $y3 $z3  $t3$LF";
	}

	print "$LF";
#Crystalクラス中の、原子の種類 AtomTypeクラスのリストをとる
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	print "nAtomType: $nAtomType$LF";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $label    = $AtomTypeList[$i]->Label();
		my $atomname = $AtomTypeList[$i]->AtomName();
		my $charge   = $AtomTypeList[$i]->Charge();
		my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
		my $AtomicMass   = $AtomTypeList[$i]->AtomicMass();
		my $i1 = $i+1;
		print "   #$i1: $label : $atomname [$AtomicNumber] ($charge) "
			."($AtomicMass)$LF";
	}

	print "$LF";
#Crystalクラス中の、非対称単位中の原子 AsymmetricAtomSiteクラスのリストをとる
	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
	print "nAsymmetricAtomSite: $nAsymmetricAtomSite$LF";
	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
		my $label     = $AtomSiteList[$i]->Label();
		my $type      = $AtomSiteList[$i]->AtomName();
		my ($x,$y,$z) = $AtomSiteList[$i]->Position(1);
		my $occupancy = $AtomSiteList[$i]->Occupancy();
		my $iAtomType = $Crystal->FindiAtomType($type);

		print "   $label ($type)[#$iAtomType]: ($x, $y, $z) [$occupancy]$LF";
	}

	print "$LF";
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
	my @nMultiplicityExpandedAtomSiteList 
		= $Crystal->GetCnMultiplicityExpandedAtomSiteList();
	print "nTotalExpandedAtomSite: $nTotalExpandedAtomSite$LF";
	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();
		my $mult = $nMultiplicityExpandedAtomSiteList[$id-1];
		my $i1 = $i+1;
		print "   $i1: [$id]$label ($type): ($x, $y, $z) [$occupancy][$mult]$LF";
	}
	print "$LF";

	return $Crystal;
}