#!/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 Utils;
use Crystal::WIEN2k;

	my $DataDir = ($ARGV[0] eq '')? "." : $ARGV[0];
	$DataDir = Utils::ReduceDirectory($DataDir);
	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($DataDir);
	my $StructFile = Deps::MakePath($DataDir, "$filename.struct", 0);
	if(!-e $StructFile) {
		my $fmask = Deps::MakePath($DataDir, "*.struct");
		my @files = glob($fmask);
		if($files[0] =~ /setrmt\.struct$/) {
			$StructFile = $files[1];
		}
		else {
			$StructFile = $files[0];
		}
		unless($StructFile) {
			print("Struct file is not found.\n");
			return 0;
		}
	}
	($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($StructFile);
	my $Title   = $filebody;
	$Title =~ s/\s/_/g;

	my $KListPath = "$filebody.klist";
	my $div = 60;
	my $K = 4;
	my ($nx, $ny, $nz) = (10, 10, 10);
	my ($emin, $emax) = (-2.0, 1.5);

	my $WIEN2k = new WIEN2k;
	my $Crystal = $WIEN2k->ReadStructFile($StructFile, 0);
	unless($Crystal) {
		print("Error: Can not read [$StructFile].\n");
		return 0;
	}
	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
print "SPG: $SPGName\n";

	my $T = new Math::MatrixReal(3, 3);
	if($SPGName =~ /^\s*F/i) {
		$T = Math::MatrixReal->new_from_cols( 
			[ [ 0.5, 0.5,   0], 
			  [ 0,   0.5, 0.5],
			  [ 0.5,   0, 0.5] ] );
		$T->transpose($T);
	}
	elsif($SPGName =~ /^\s*I/i) {
		$T = Math::MatrixReal->new_from_cols( 
			[ [-0.5, 0.5, 0.5], 
			  [ 0.5,-0.5, 0.5],
			  [ 0.5, 0.5,-0.5] ] );
		$T->transpose($T);
	}
	elsif($SPGName =~ /^\s*A/i) {
		$T = Math::MatrixReal->new_from_cols( 
			[ [ 1.0, 0.0, 0.0], 
			  [ 0.0, 0.5, 0.5],
			  [ 0.0,-0.5, 0.5] ] );
		$T->transpose($T);
	}
	elsif($SPGName =~ /^\s*B/i) {
		$T = Math::MatrixReal->new_from_cols( 
			[ [ 0.5, 0.0, 0.5], 
			  [ 0.0, 1.0, 0.0],
			  [-0.5, 0.0, 0.5] ] );
		$T->transpose($T);
	}
	elsif($SPGName =~ /^\s*C/i) {
		$T = Math::MatrixReal->new_from_cols( 
			[ [ 0.5, 0.5, 0.0], 
			  [-0.5, 0.5, 0.0],
			  [ 0.0, 0.0, 1.0] ] );
		$T->transpose($T);
	}

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

	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
	print "Lattice vector: ($a11, $a12, $a13)\n";
	print "                ($a21, $a22, $a23)\n";
	print "                ($a31, $a32, $a33)\n";
	($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->ReciprocalLatticeVectors();
	print "RLattice vector: ($a11, $a12, $a13)\n";
	print "                 ($a21, $a22, $a23)\n";
	print "                 ($a31, $a32, $a33)\n";

	my $PrimCrystal = $Crystal->ConvertLattice($T, 0, 1);
	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $PrimCrystal->LatticeVectors();
	print "Lattice vector: ($a11, $a12, $a13)\n";
	print "                ($a21, $a22, $a23)\n";
	print "                ($a31, $a32, $a33)\n";
	($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $PrimCrystal->ReciprocalLatticeVectors();
	print "RLattice vector: ($a11, $a12, $a13)\n";
	print "                 ($a21, $a22, $a23)\n";
	print "                 ($a31, $a32, $a33)\n";
exit;

	open(OUT, ">$KListPath") or die "$!: Can not write to [$KListPath]\n";

	my $count = 1;
	my $BandName = "";
	for(my $ix = 0 ; $ix <= $nx ; $ix++) {
		for(my $iy = 0 ; $iy <= $ny ; $iy++) {
			for(my $iz = 0 ; $iz <= $nz ; $iz++) {
				if($iz == 0) {
					$BandName = sprintf("B%02d%02d", $ix+1, $iy+1);
				}
				else {
					$BandName = "";
				}

				my $x = $K * $div * $ix / 2 / $nx;
				my $y = $K * $div * $iy / 2 / $ny;
				my $z = $K * $div * $iz / 2 / $nz;
				if($count == 1) {
					printf(OUT "%5s       %3d  %3d  %3d  %3d%4.1f %4.1f %4.1f\n",
						$BandName, $x, $y, $z, $div, 1.0, $emin, $emax);
				}
				else {
					printf(OUT "%5s       %3d  %3d  %3d  %3d%4.1f\n",
						$BandName, $x, $y, $z, $div, 1.0);
				}
				$count++;
			}
		}
	}
	printf(OUT "END\n\n");
	close(OUT);
	exit;