#!/usr/bin/perl

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;

use Deps;
use Utils;
use JFile;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::WIEN2k;

my $SGROUPPath   = ($^O =~ /^MSWin/)? "D:/Programs/Perl/WIEN2k/sgroup.exe" : "sgroup";

my $InputFile         = "Si-final.cif";
my ($drive, $dir, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($InputFile);
my $SampleName    = $filebody;

my $SimplifiedCIFFile = $InputFile;
$SimplifiedCIFFile =~ s/\.cif/-simple\.cif/i;
my $OutputFile    = $InputFile;
$OutputFile =~ s/\.cif/-converted\.cif/i;
my $StructPath    = $InputFile;
$StructPath       =~ s/\.cif/\.struct/i;
my $SymStructPath = $InputFile;
$SymStructPath    =~ s/\.cif/\-Symmetrized.struct/i;
my $SymCIFPath    = $InputFile;
$SymCIFPath       =~ s/\.cif/\-Symmetrized.cif/i;

#基本ベクトルの変換行列
my $T    = new Math::MatrixReal(3, 3);
my $sT   = new Math::MatrixReal(3, 3);
my $sxyz = new Math::MatrixReal(3, 1);
my $shkl = new Math::MatrixReal(3, 1);

print "<H1>Lattice conversion</H1>\n";

my $CheckConsistency = 1; #&GetArg('CheckConsistency');
my $ParametersSelect = "ai"; #ai, xyz, a*i, hkl
my $DirectionSelect  = "OriginalToConverted"; #OriginalToConverted, ConvertedlToOriginal
my $PresetRule = "PrimFCC"; #RhombHex, HexRhomb, HexOrtho, FCCPrim, BCCPrim, ACenterPrim, BCenterPrim, CCenterPrim
print "<b>PresetRule:</b> $PresetRule\n";
print "<b>Direction:</b> $DirectionSelect\n";
print "<b>Check consistency:</b>$CheckConsistency\n";

	my $a     = 0; #&GetArg('a'); $a = eval($a);
	my $b     = 0; #&GetArg('b'); $b = eval($b);
	my $c     = 0; #&GetArg('c'); $c = eval($c);
	my $alpha = 0; #&GetArg('alpha'); $alpha = eval($alpha);
	my $beta  = 0; #&GetArg('beta');  $beta  = eval($beta);
	my $gamma = 0; #&GetArg('gamma'); $gamma = eval($gamma);

	my $v = 0; #&GetArg('sx'); $v = eval($v);
	$sxyz->assign(1, 1, $v);
	$v = 0; #&GetArg('sy'); $v = eval($v);
	$sxyz->assign(2, 1, $v);
	$v = 0; #&GetArg('sz'); $v = eval($v);
	$sxyz->assign(3, 1, $v);

	$v = 0; #&GetArg('sh'); $v = eval($v);
	$shkl->assign(1, 1, $v);
	$v = 0; #&GetArg('sk'); $v = eval($v);
	$shkl->assign(2, 1, $v);
	$v = 0; #&GetArg('sl'); $v = eval($v);
	$shkl->assign(3, 1, $v);

my $CIF = new CIF();
unless($CIF->Read($InputFile)) {
	print "<H2>Error: Can not read $InputFile.</H2>\n";
	return -2;
}
if($CIF->WriteSimpleCIFFile($SimplifiedCIFFile)) {
}

my $Crystal = $CIF->GetCCrystal();
$Crystal->ExpandCoordinates();
my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "Lattice parameters:\n";
print "a=$a b=$b c=$c  alpha=$alpha beta=$beta gamma=$gamma\n";
print "Volume of original lattice: ", $Crystal->Volume(), " A<sup>3</sup>\n";
print "Density of original lattice: ", $Crystal->Density(), " g/cm<sup>3</sup>\n";
print "\n";

my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();
print "<b>Real space metrics conversion:</b>\n";
printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;

&LatticeConversion();

exit;


sub LatticeConversion
{
	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') {
		$sT = Math::MatrixReal->new_from_cols( 
				[ [-0.5, 0.5, 0.5], 
				  [ 0.5,-0.5, 0.5],
				  [ 0.5, 0.5,-0.5] ] );
		$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 {
		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 "<h2>Conversion rule:</h2>\n";
	my $IsChosen = 0;
	if($ParametersSelect eq 'ai') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "<b>Real space vectors:</b> (<b>a'<sub>i</sub></b>) = "
				."(t<sub>ij</sub>)(<b>a<sub>i</sub></b>)\n";
			$T = $sT;
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "<b>Real space vectors:</b> (<b>a<sub>i</sub></b>) = "
				."(t<sub>ij</sub>)(<b>a'<sub>i</sub></b>)\n";
			$T = $sT->inverse();
			$IsChosen = 1;
		}
	}
	elsif($ParametersSelect eq 'xyz') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "<b>Real space coordinates:</b> (x'<sub>i</sub>) = "
				."(t<sub>ij</sub>)(x<sub>i</sub>)\n";
			$T->transpose($sT);
			$T = $T->inverse();
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "<b>Real space vectors:</b> (x<sub>i</sub>) = "
				."(t<sub>ij</sub>)(x'<sub>i</sub>)\n";
			$T->transpose($sT);
			$IsChosen = 1;
		}
	}
	elsif($ParametersSelect eq 'a*i') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "<b>Reciprocal space vectors:</b> (<b>a'*<sub>i</sub></b>) = "
				."(t<sub>ij</sub>)(<b>a*<sub>i</sub></b>)\n";
			$T->transpose($sT);
			$T = $T->inverse();
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "<b>Reciprocal space vectors:</b> (<b>a*<sub>i</sub></b>) = "
				."(t<sub>ij</sub>)(<b>a'*<sub>i</sub></b>)\n";
			$T->transpose($sT);
			$IsChosen = 1;
		}
	}
	elsif($ParametersSelect eq 'hkl') {
		if($DirectionSelect eq 'OriginalToConverted') {
			print "<b>Reciprocal space coordinates:</b> (h'<sub>i</sub>) = "
				."(t<sub>ij</sub>)(h<sub>i</sub>)\n";
			$T = $sT;
			$IsChosen = 1;
		}
		elsif($DirectionSelect eq 'ConvertedlToOriginal') {
			print "<b>Reciprocal space coordinates:</b> (h<sub>i</sub>) = "
				."(t<sub>ij</sub>)(h'i</sub>)\n";
			$T = $sT->inverse();
			$IsChosen = 1;
		}
	}
	unless($IsChosen) {
		print "Selected option {Parameters: $ParametersSelect} in {$DirectionSelect}"
			." direction is not supported.\n";
		return 1;
	}
	
	print "<b>Input matrix:</b> (t<sub>ij</sub>)\n";
	print "<pre>\n";
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$sT->element(1,1), $sT->element(1,2), $sT->element(1,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$sT->element(2,1), $sT->element(2,2), $sT->element(2,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$sT->element(3,1), $sT->element(3,2), $sT->element(3,3);
	print "</pre>\n";
	print "<hr>\n";

#実空間べクトル(ai)の変換行列
	print "<b>Real space vector:</b> (<b>a'<sub>i</sub></b>) = "
		."(T<sub>ij</sub>)(<b>a<sub>i</sub></b>)\n";
	print "<pre>\n";
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$T->element(1,1), $T->element(1,2), $T->element(1,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$T->element(2,1), $T->element(2,2), $T->element(2,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$T->element(3,1), $T->element(3,2), $T->element(3,3);
	print "</pre>\n";

#実空間座標(x,y,z)の変換行列
	my $tRT = new Math::MatrixReal(3, 3);
	$tRT->transpose($T);
	$tRT = $tRT->inverse();
	print "<b>Real space coordinate conversion matrix</b> (x') = "
		."((T<sub>ij</sub>)<sup>t</sup>)<sup>-1</sup>(x)\n";
	print "<pre>\n";
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$tRT->element(1,1), $tRT->element(1,2), $tRT->element(1,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$tRT->element(2,1), $tRT->element(2,2), $tRT->element(2,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$tRT->element(3,1), $tRT->element(3,2), $tRT->element(3,3);
	print "</pre>\n";

#逆空間ベクトル(a*i)の変換行列: $tRT
	my $RT = new Math::MatrixReal(3, 3);
	$RT->transpose($T);
	$RT = $RT->inverse();
	print "<b>Reciprocal space vector conversion matrix</b> (<b>a<sup>*</sup>'</b>) = "
		."((T<sub>ij</sub>)<sup>t</sup>)<sup>-1</sup>(<b>a<sup>*</sup></b>)\n";
	print "<pre>\n";
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$RT->element(1,1), $RT->element(1,2), $RT->element(1,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$RT->element(2,1), $RT->element(2,2), $RT->element(2,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$RT->element(3,1), $RT->element(3,2), $RT->element(3,3);
	print "</pre>\n";
	print "\n";

#逆空間座標(h k l)の変換行列: $T
	my $tR = new Math::MatrixReal(3, 3);
	$tR = $T;
	print "<b>Reciprocal space coordinate conversion matrix</b> (h') = (T<sub>ij</sub>)(h)\n";
	print "<pre>\n";
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$tR->element(1,1), $tR->element(1,2), $tR->element(1,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$tR->element(2,1), $tR->element(2,2), $tR->element(2,3);
	printf "     |%10.4f %10.4f %10.4f|\n", 
		$tR->element(3,1), $tR->element(3,2), $tR->element(3,3);
	print "</pre>\n";
	print "\n";
	print "<hr>\n";

#座標(sxyz)の変換: tRT
	if($sxyz->element(1,1) ne '') {
		my $a1 = 
		my $txyz = $tRT * $sxyz;
		print "<b>Real space coordinate conversion:</b>\n";
		printf "(%10.4f, %10.4f, %10.4f) => (%10.4f, %10.4f, %10.4f)\n",
			$sxyz->element(1,1), $sxyz->element(2,1), $sxyz->element(3,1),
			$txyz->element(1,1), $txyz->element(2,1), $txyz->element(3,1);
		print "\n";
	}

#逆空間座標(shkl)の変換: tR
	if($shkl->element(1,1) ne '') {
		my $thkl = $tR * $shkl;
		print "<b>Reciprocal space coordinate conversion:</b>\n";
		printf "(%10.4f, %10.4f, %10.4f) => (%10.4f, %10.4f, %10.4f)\n",
			$shkl->element(1,1), $shkl->element(2,1), $shkl->element(3,1),
			$thkl->element(1,1), $thkl->element(2,1), $thkl->element(3,1);
		print "\n";
	}

#格子定数の変換
	if($a ne '') {
		print "<b>Original lattice parameters:</b>\n";
		printf "a=%12.6f b=%12.6f c=%12.6f  alpha=%8.4f beta=%8.4f gamma=%8.4f\n",
			$a, $b, $c, $alpha, $beta, $gamma;
		my $Crystal = new Crystal;
		$Crystal->SetLatticeParameters($a,$b,$c,$alpha,$beta,$gamma);

		my $NewCrystal = $Crystal->ConvertLattice($T, 1);

		my ($a2,$b2,$c2,$alpha2,$beta2,$gamma2) = $NewCrystal->LatticeParameters();
		print "<b>Converted lattice parameters:</b>\n";
		printf "a'=%12.6f b'=%12.6f c'=%12.6f  alpha'=%8.4f beta'=%8.4f gamma'=%8.4f\n",
			$a2, $b2, $c2, $alpha2, $beta2, $gamma2;
		print "\n";

		my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33)
			= $Crystal->Metrics();

		print "<b>Real space metrics conversion:</b>\n";
		print "<pre>\n";
		printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
		printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
		printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;
		print "</pre>\n";
		print "=>\n";

		($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33)
			= $NewCrystal->Metrics();
		print "<pre>\n";
		printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
		printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
		printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;
		print "</pre>\n";
	}

	{
		print "<hr>\n";
		print "Convert CIF File ($InputFile).\n";
		print "<b>Links:</b>\n";

		if($InputFile !~ /\.cif$/i) {
			print "<H2>Error: ($InputFile) is not a CIF file.</H2>\n";
			return;
		}



my $NewCrystal = $Crystal->ConvertLattice($T, 0, $CheckConsistency);
my ($a2,$b2,$c2,$alpha2,$beta2,$gamma2) = $NewCrystal->LatticeParameters();
($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $NewCrystal->Metrics();
printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;

		print "<b>Converted Lattice parameters:</b>\n";
		print "a'=$a2 b'=$b2 c'=$c2  alpha'=$alpha2 beta'=$beta2 gamma'=$gamma2\n";
		print "<b>Volume of converted lattice:</b> ", $NewCrystal->Volume(), " A<sup>3</sup>\n";
		print "<b>Density of converted lattice:</b> ", $NewCrystal->Density(), 
					" g/cm<sup>3</sup>\n";

		my $NewCIF = new CIF();
		$NewCIF->SetFileName($OutputFile);
		$NewCIF->SetCrystalName($Crystal->CrystalName());
		$NewCIF->SetLatticeParameters($a2, $b2, $c2, $alpha2, $beta2, $gamma2);
		$NewCIF->SetSpaceGroup("P 1", 1);
		$NewCIF->SetVolume($NewCrystal->Volume());
		$NewCIF->AddSymmetryOperation("x,y,z");
		my @ExpandedAtomSiteList = $NewCrystal->GetCExpandedAtomSiteList();
		my $nTotalExpandedAtomSite = $NewCrystal->nTotalExpandedAtomSite();
		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();

			$NewCIF->AddAsymmetricAtomSite($label, $atomname, $x, $y, $z,$occ);
		}
		$NewCIF->FillCIFData();

		unless($NewCIF->WriteSimpleCIFFile($OutputFile)) {
			print "<H2>Error: Could not write to ($OutputFile).</H2>\n";
			return -1;
		}

		if(abs($Crystal->Density() - $NewCrystal->Density()) / $Crystal->Density()
				> 0.05) {
			print "<H1>Causion!!.</H1>\n";
			print "<H2>Densities of the original and converted cells are "
				."different.</H2>\n";
			print "<H2>The conversion rule employed would not be compatible with "
				."the original cell.</H2>\n";
		}

		if($PresetRule eq 'HexOrtho' or 
		   $PresetRule eq 'FCCPrim' or
		   $PresetRule eq 'BCCPrim' or
		   $PresetRule eq 'PrimFCC' or
		   $PresetRule eq 'PrimBCC'
		   ) {

			print "<H2>Find symmetry for the converted lattice</H2>";

			my ($a0,$b0,$c0,$alpha0,$beta0,$gamma0) = $NewCrystal->LatticeParameters();
			if($PresetRule eq 'HexOrtho') {
				$NewCrystal->SetLatticeParameters($a0, $a0*1.3, $c0,$alpha0,$beta0,$gamma0);
			}
			elsif($PresetRule eq 'FCCPrim' or $PresetRule eq 'BCCPrim') {
				$NewCrystal->SetLatticeParameters($a0,$b0, $c0, 65,65,65);
			}

			print "<H3>MakeStruct: Save to $StructPath</H3>";
			my $WIEN = new WIEN2k;
			$WIEN->SetSampleName($SampleName);
			$WIEN->SaveStructFile($NewCrystal, $StructPath, 0, 0, 0, 0);

			print "<H3>Execute sgroup: Save to $SymStructPath</H3>";
			print("$SGROUPPath -wi $StructPath\n");
			open(PIPEIN, "$SGROUPPath -wi $StructPath |");
			open(OUT, ">$StructPath.out");
			while(<PIPEIN>) {
				print OUT $_;
			}
			close(OUT);
			close(PIPEIN);

			print("$SGROUPPath -wi -wo $StructPath $SymStructPath\n");
			system("$SGROUPPath -wi -wo $StructPath $SymStructPath");

			my $WIEN2k2 = new WIEN2k;
			my $SymCrystal = $WIEN2k2->ReadStructFile($SymStructPath, 0);
			unless($SymCrystal) {
				print("Error: Can not read [$SymStructPath].\n");
				return 0;
			}

			my $SPG = $SymCrystal->GetCSpaceGroup();
			my $SPGName = $SPG->SPGName();
			if($SPGName =~ /^R/i or $SPGName =~ /^P\s*\-?3/i) {
				$SPG->SetSPGName("$SPGName R");
				my ($a,$b,$c,$alpha,$beta,$gamma) = $SymCrystal->LatticeParameters();
				print("  Hexagonal axis: a=$a   c=$c\n");
				my $ar = sqrt(3.0*$a*$a + $c*$c) / 3.0;
				my $alphah = 1.5 / sqrt(3.0 + ($c/$a)*($c/$a));
				$alphah = 2.0 * Sci::asin($alphah) * Sci::todeg();
				$alphah = $alpha0 if(abs($alphah -65.0) < 1.0e-3);
				print("  Rhombohedral axis: a=$ar   alpha=$alphah\n");
				$SymCrystal->SetLatticeParameters($ar,$ar,$ar,$alphah,$alphah,$alphah);
			}
			my $CIF2 = new CIF;
			if(!$CIF2->CreateCIFFileFromCCrystal($SymCrystal, $SymCIFPath, 0, 'unix')) {
				print("Error: Can not create [$SymCIFPath].\n");
				return 0;
			}
		}
	}

	return 1;
}