#!/usr/bin/perl

use lib 'd:/tkProg.main/tklib/Perl/lib';
#use lib 'd:/Programs/Perl/lib';
#use lib 'c:/Programs/Perl/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);

#===============================================
# コマンドラインオプション読み込み
#===============================================
$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);
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";

&Execute();

print("\nPress ENTER to terminate>>");
<STDIN>;
exit 0;

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

#==========================================
# サブルーチン
#==========================================
sub ConfigureIniFileVariables
{
	$App->AddIniFileVariable("\\Preferences\\WorkDir",        "WorkDir");
	$App->ReadSetting();
}

sub Execute
{
#CIFクラスの内容から、Crystalクラスを作成
#作成した時点で、Metricsの計算、イオン情報(原子量など)の読み込み、
#格子体積の計算は終わっている
	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";

	if(!$DoSuperlattice) {
		print("\nPress ENTER to terminate>>");
		<STDIN>;
	}


#$Crystalをもとに、超格子の結晶格子$SuperLatticeを作る
	my ($nx,$ny,$nz) = (2,2,2);
	$Crystal->CreateSuperLattice($nx,$ny,$nz);
	my $SuperLattice = $Crystal->GetCSuperLattice();

	print "$LF";
	print "SuplaerLattice $nx x $ny x $nz$LF";
	my ($a,$b,$c,$alpha,$beta,$gamma) = $SuperLattice->LatticeParameters();
	print "Cell: $a $b $c  $alpha $beta $gamma$LF";

	my @ExpandedAtomSiteList = $SuperLattice->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $SuperLattice->nTotalExpandedAtomSite();
	my @nMultiplicityExpandedAtomSiteList 
		= $SuperLattice->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";
	}
}