#!/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 strict;
#use warnings;
use Math::MatrixReal;
use Math::Complex;

use Sci qw($pi $a0 $torad $todeg);# acos asin);
use Sci::Algorism;
use Sci::MyMatrixReal;

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

my $InFile         = ($ARGV[0] eq '')? 'FCC-PrimitiveCell.Cartesian.KPOINTS' : $ARGV[0];
#my $InFile        = ($ARGV[0] eq '')? 'FCC-PrimitiveCell.Reciprocal.KPOINTS' : $ARGV[0];
my $OutFile        = ($ARGV[1] eq '')? "$InFile.conved" : $ARGV[1];
#No P1 RhombHex HexRhomb HexOrtho FCCPrim BCCPrim A/B/CCenterPrim
my $BravaisLattice = ($ARGV[2] eq '')? "FCCPrim" : $ARGV[2];

if($ARGV[0] eq '--help') {
	&usage();
	exit;
}

my ($InputMode, $pContent) = &ReadKPOINTS($InFile);

my $Direction;
if($InputMode eq 'Cartesian') {
	$Direction = 'OriginalToConverted';
}
else {
	$Direction = 'ConvertedlToOriginal';
}

# 1st: No P1 RhombHex HexRhomb HexOrtho FCCPrim BCCPrim A/B/CCenterPrim
# 2nd: OriginalToConverted ConvertedlToOriginal
# 3rd: ai xyz a*i hkl
my ($sT, $T, $tRT, $RT, $tR) = Crystal->GetLatticeConversionMatrix($BravaisLattice, $Direction, 'ai', 1);
#1st argで指定された変換行列  : $sT
#2nd argで選択された変換行列  : $T
#実空間座標(x,y,z)の変換行列  : $tRT = transpose(T)
#逆空間ベクトル(a*i)の変換行列: $RT  = transpose(T)
#逆空間座標(h k l)の変換行列  : $tR  = $T
print "Conv M specified by 1st arg: sT=\n", $sT;
print "Conv M selected by 2nd arg :  T=\n", $T;
#print "Conv M for (x,y,z)         : RT=\n", $RT;
#print "Conv M for (a*_i)          :tRT=\n", $tRT;
#print "Conv M for (h k l)         : tT=\n", $tR;

#my $GT = $T; # equal to $RT?
#$GT->transpose($GT);
#$GT = $GT->inverse();
#print "GT=\n";
#print $GT;

print "\n";
print "Conversion for [$BravaisLattice]\n";
print "Input data is for [$InputMode]\n";

my $out = JFile->new($OutFile, 'w') or die "$!: Can not write to [$OutFile]";
$out->print($pContent->[0]);
$out->print($pContent->[1]);
$out->print($pContent->[2]);
if($InputMode eq 'Cartesian') {
	$out->print("Reciprocal\n");
}
else {
	$out->print("Cartesian\n");
}

for(my $i = 4 ; $i < @$pContent ; $i++) {
#print "$i: $pContent->[$i]";
	my ($kx, $ky, $kz, $name) = Utils::Split("\\s+", $pContent->[$i]);
	if($kz eq '') {
		$out->print($pContent->[$i]);
		next;
	}

	print " k[$i] = ($kx, $ky, $kz) [$name]\n";

	my $hkl = new Math::MatrixReal(3, 1);
	$hkl = Math::MatrixReal->new_from_cols([ [$kx, $ky, $kz] ]);

	my $chkl = $T * $hkl;
#	my $chkl = $RT * $hkl;
	$out->printf("%10.6f %10.6f %10.6f  %s\n", $chkl->element(1,1), $chkl->element(2,1), $chkl->element(3,1), $name);
print $chkl;
}

$out->Close();

exit;


#======================================================================
#
#======================================================================
sub ReadKPOINTS
{
	my ($InFile) = @_;

print "Read [$InFile]\n";
	my $in = JFile->new($InFile, 'r') or die "$!: Can not read [$InFile].\n";
	my @content;
	for(my $i = 0 ; ; $i++) {
		$content[$i] = $in->ReadLine();
#print "$i:$content[$i]";
		last if($in->eof());

	}
	$in->Close();
	my $mode = ($content[3] =~ /^C/i)? 'Cartesian' : 'Reciprocal';

	return ($mode, \@content);
}

sub usage
{
print "\n";
print "perl ConvertKPOINTS.pl InFile OutFile InputMode\n";
print "   InputMode: [No P1 RhombHex HexRhomb HexOrtho FCCPrim BCCPrim A/B/CCenterPrim]\n";
print "\n";
}