#!/usr/bin/perl

use lib 'c:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';

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

use strict;

use MyApplication;
use Sci qw($pi);
use Crystal::CIF;
use Crystal::Crystal;
use Crystal::PointGroup;
use Crystal::SpaceGroup;

#use warnings;

#===============================================
# Global変数
#===============================================

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

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--TolD",                "--TolD=val: Torelance for atomic position (A) [Def 0.01 A]", "0.01");
$App->AddArgument("--xc0",                 "--xc0=val [Def 0.0]",                 "0.0");
$App->AddArgument("--yc0",                 "--yc0=val [Def 0.0]",                 "0.0");
$App->AddArgument("--zc0",                 "--zc0=val [Def 0.0]",                 "0.0");
$App->AddArgument("--ScanAll",             "--ScanAll=[0|1] [Def 0]",             "0");
$App->AddArgument("--FindSiteSymmetry",    "--FindSiteSymmetry=[0|1] [Def 0]",    "0");
$App->AddArgument("--ShowCrystalInf",      "--ShowCrystalInf=[0|1] [Def 0]",      "0");
$App->AddArgument("--ShowAllCoordination", "--ShowAllCoordination=[0|1] [Def 0]", "0");
$App->AddArgument("--RCoordMax",           "--RCoordMax=val [Def 7.0]",           "7.0");
$App->AddArgument("--ShowFoundStructure",  "--ShowFoundStructure=[0|1] [Def 0]",  "0");
$App->AddArgument("--help",                "--help: Show this help");
exit 1 if($App->ReadArgs() != 1);

my $OverlapCriteriaRadius = $App->GetGetArg('TolD');
   $OverlapCriteriaRadius = 0.01 if(!defined $OverlapCriteriaRadius);
my $xc0                   = $App->GetGetArg('xc0');
   $xc0                   = 0.0 if(!defined $xc0);
my $yc0                   = $App->GetGetArg('zc0');
   $yc0                   = 0.0 if(!defined $yc0);
my $zc0                   = $App->GetGetArg('zc0');
   $zc0                   = 0.0 if(!defined $zc0);
my $ScanAll               = $App->GetGetArg('ScanAll');
   $ScanAll               = 0 if(!defined $ScanAll);
my $FindSiteSymmetry      = $App->GetGetArg('FindSiteSymmetry');
   $FindSiteSymmetry      = 0 if(!defined $FindSiteSymmetry);
my $ShowCrystalInf        = $App->GetGetArg('ShowCrystalInf');
   $ShowCrystalInf        = 0 if(!defined $ShowCrystalInf);
my $ShowAllCoordination   = $App->GetGetArg('ShowAllCoordination');
   $ShowAllCoordination   = 0 if(!defined $ShowAllCoordination);
my $RCoordMax             = $App->GetGetArg('RCoordMax');
   $RCoordMax             = 7.0 if(!defined $RCoordMax);
my $ShowFoundStructure    = $App->GetGetArg('ShowFoundStructure');
   $ShowFoundStructure    = 0 if(!defined $ShowFoundStructure);
my $CIFFile               = $App->GetGetArg(0);
   $CIFFile               = "NaCl-converted.cif" if(!defined $CIFFile);

#CIFクラスを取得
my $cif = new CIF;
#CIFファイル名を引数に与えて読み込む
unless($cif->Read($CIFFile)) {
	print "Error: Can not read [$CIFFile]\n\n";
	exit 2;
}
print "CIF                  : ", $cif->FileName(), "\n";
print "Initial center       : ($xc0, $yc0, $zc0)\n";
print "ShowCrystalInf       : $ShowCrystalInf\n";
print "ShowAllCoordination  : $ShowAllCoordination\n";
print "ShowFoundStructure   : $ShowFoundStructure\n";
print "ScanAll              : $ScanAll\n";
print "FindSiteSymmetry     : $FindSiteSymmetry\n";
print "OverlapCriteriaRadius: $OverlapCriteriaRadius A\n";
print "RCoordMax            : $RCoordMax A\n";

my $IsPrint = 0;
my $Crystal = $cif->GetCCrystal();
&MyFindPossibleSPGs($Crystal, $CIFFile,
	$ShowCrystalInf, $ShowAllCoordination, $ShowFoundStructure, $ScanAll,
	$FindSiteSymmetry, $xc0, $yc0, $zc0, $OverlapCriteriaRadius, $ShowAllCoordination, $RCoordMax, $IsPrint);

exit 0;

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

#==========================================
# サブルーチン
#==========================================
sub MyFindPossibleSPGs
{
	my ($Crystal, $CIFFile, $ShowCrystalInf, $ShowAllCoordination, $ShowFoundStructure, $ScanAll,
		$FindSiteSymmetry, $xc0, $yc0, $zc0, 
		$OverlapCriteriaRadius, $ShowAllCoordination, $RCoordMax, $IsPrint) = @_;

#CIFクラスの内容から、Crystalクラスを作成
#作成した時点で、Metricsの計算、イオン情報(原子量など)の読み込み、
#格子体積の計算は終わっている

	print "\n";
	my $CrystalName = $Crystal->CrystalName();
	print "CrystalName: $CrystalName\n";
	my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	print "Lattice paramters: $a $b $c  $alpha $beta $gamma\n";
#	my ($g11,$g12,$g13,$g21,$g22,$g23,$g31,$g32,$g33) = $Crystal->Metrics();
#	print "Metrics: $g11 $g22 $g33  $g12 $g23 $g31\n";
#	my $vol = $Crystal->Volume();
#	print "Volume: $vol A^3\n";
#対称要素を使い、非対称単位中の原子を展開する
#この際、Densityも計算される
	$Crystal->ExpandCoordinates();
#	my $Density = $Crystal->Density();
#	print "Density: $Density g/cm^3\n";

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

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

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

#Crystalクラス中の、原子の種類 AtomTypeクラスのリストをとる
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	print "nAtomType: $nAtomType\n";
	if($ShowCrystalInf) {
		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)\n";
		}
		print "\n";
	}
	
#Crystalクラス中の、非対称単位中の原子 AsymmetricAtomSiteクラスのリストをとる
	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
	print "nAsymmetricAtomSite: $nAsymmetricAtomSite\n";
	if($ShowCrystalInf or $FindSiteSymmetry) {
		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]\n";

			my ($ret, $pMatchedSyms) = $Crystal->FindSiteSymmetries($x, $y, $z, $OverlapCriteriaRadius, $IsPrint) 
							if($FindSiteSymmetry);
		}
	}

	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
	print "nExpanededAtomSite: $nTotalExpandedAtomSite\n";
	if($ShowCrystalInf) {
		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 $iAtomType = $Crystal->FindiAtomType($type);

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

	my ($pMostPossibleSPG, $pPossibleSPGs) =
		$Crystal->FindPossibleSPGs($ShowCrystalInf, $ShowAllCoordination, $ShowFoundStructure, $ScanAll,
			$FindSiteSymmetry, $xc0, $yc0, $zc0, 
			$OverlapCriteriaRadius, $ShowAllCoordination, $RCoordMax, $IsPrint);

	my ($drive, $dir, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($CIFFile);
	my $TranslatedCIFFile = "$filebody-Translated.cif";
	if(!CIF->new()->CreateCIFFileFromCCrystal($Crystal, $TranslatedCIFFile, 0, 'unix')) {
		print("Error: Can not create [$TranslatedCIFFile].\n");
		return 0;
	}

	my $MostProbableiSPG   = $pMostPossibleSPG->{iSPG};
	my $MostProbableiSet   = $pMostPossibleSPG->{iSet};
	my $SPGName            = $pMostPossibleSPG->{SPGName};
	my $LatticeSystem      = $pMostPossibleSPG->{LatticeSystem};
	my $nTranslation       = $pMostPossibleSPG->{nTranslation};
	my $nSymmetryOperation = $pMostPossibleSPG->{nSymmetryOperation};

	print "\n";
	print "Most matched SPG $SPGName (#$MostProbableiSPG-$MostProbableiSet LS=$LatticeSystem)\n";
	print "   nSym=$nSymmetryOperation  nTrans=$nTranslation\n";
	print "\n";

	for(my $i = 0 ; $i < @$pPossibleSPGs ; $i++) {
		my $SPG                = $pPossibleSPGs->[$i]{pSPG};
		my $iSPG               = $pPossibleSPGs->[$i]{iSPG};
		my $iSet               = $pPossibleSPGs->[$i]{iSet};
		my $SPGName            = $pPossibleSPGs->[$i]{SPGName};
		my $LatticeSystem      = $pPossibleSPGs->[$i]{LatticeSystem};
		my $nTranslation       = $pPossibleSPGs->[$i]{nTranslation};
		my $nSymmetryOperation = $pPossibleSPGs->[$i]{nSymmetryOperation};
		my $pAsymSites         = $pPossibleSPGs->[$i]{pAsymSites};
		my $nASSite            = $pPossibleSPGs->[$i]{nASSite};

		print "SPG $SPGName (#$iSPG-$iSet LS=$LatticeSystem) is matched\n";
		print "   nSym=$nSymmetryOperation  nTrans=$nTranslation\n";
		print "   $nASSite asymmetric sites are found.\n";

		my $NewCrystal = new Crystal();
		$NewCrystal->SetCrystalName($CrystalName);
		$NewCrystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
		$NewCrystal->SetCSpaceGroup($SPG);

		for(my $i = 0 ; $i < @$pAsymSites ; $i++) {
#print "SiteId=$pSiteId->[$i]\n";
#					next if(defined $pSiteId->[$i]);

			my $site      = $pAsymSites->[$i];
			my $label     = $site->Label();
			my $atomname  = $site->AtomName();
			my ($x,$y,$z) = $site->Position(1);
			my $occupancy = $site->Occupancy();
			print "   $i: $atomname: ($x, $y, $z) occ=$occupancy\n" if($ShowFoundStructure);

			$NewCrystal->AddAtomSite($label, $atomname, $x, $y, $z, $occupancy);
		}

		my $NewCIFName = $CIFFile;
		$NewCIFName    =~ s/\.cif//i;
		$SPGName       =~ s/\s+//g;
		$NewCIFName   .= "-$SPGName.cif";
		print "** Save CIF file to [$NewCIFName].\n";
		CIF->new()->CreateCIFFileFromCCrystal($NewCrystal, $NewCIFName, undef, undef, $ShowFoundStructure);
	}

	return;
}

1;