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

	my $DataDir = ($ARGV[0] eq '')? "." : $ARGV[0];
	$DataDir = Utils::ReduceDirectory($DataDir);
#	my ($a, $b, $c) = (0.135, 0.135, 0.0625);
	my ($EFmin, $EFmax) = (-1.0, 1.5);
	my $Eoffset = 0.0;# 100;
	my $div     = 60;
	my ($nx, $ny, $nz) = (10, 10, 10);
	my $OneByOne = 1;

	my $EF = ($EFmin + $EFmax) / 2.0;

	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 $OutFile = "$filebody.bxsf";
	my $Title   = $filebody;
	$Title =~ s/\s/_/g;

	print("  Read $StructFile.\n");

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

	my $PrefPath = Deps::MakePath($DataDir, ".prefs-scf");
	my $ini = new IniFile();
	$ini->ReadAll($PrefPath);
	my @InFile;
	if($ini->{spin} eq '') {
		$InFile[0]  = "$filebody.spaghetti_ene";
	}
	else {
		$InFile[0]  = "$filebody.spaghettiup_ene";
		$InFile[1]  = "$filebody.spaghettiup_ene";
	}
	if(!-e $InFile[0]) {
		print("Error: Can not read [$InFile[0]]\n");
		exit;
	}

	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
print "SPG: $SPGName\n";

	my $T = new Math::MatrixReal(3, 3);
	my $matrix = $Crystal->GetMatrixConventionalToPrimitiveCell();
	$T = Math::MatrixReal->new_from_cols($matrix);
	$T->transpose($T);
	my $PrimCrystal = $Crystal->ConvertLattice($T, 0, 1);

#実空間べクトル(ai)の変換行列
if(0) {
	print "Conversion matrix for real space vector: (a'i) = (Tij)(ai)\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);
}
if(0) {
	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 ($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";
}

my $iBand = 0;
my $nBand = 0;
my @BandEnergy;
my @BandIndex;
foreach my $InFile (@InFile) {
	open(IN, "$InFile")    or die "$!: Can not read [$InFile]\n";

	while(1) {
		$iBand++;
		my $line = <IN>;
		last if(!defined $line);
		$line = Utils::DelSpace($line);
		last if($line eq '');
print "Index line: $line\n";

		my ($emin, $emax) = (1.0e10, -1.0e10);
		my @Energy = ([], [], []);
		for(my $ix = 0 ; $ix <= $nx ; $ix++) {
			for(my $iy = 0 ; $iy <= $ny ; $iy++) {
				for(my $iz = 0 ; $iz <= $nz ; $iz++) {
					my ($IsConverted, $cix, $ciy, $ciz) = $Crystal->ConvertLatticeIndexByLatticeSystem();
					if($IsConverted) {
						$Energy[$ix][$iy][$iz] = $Energy[$cix][$ciy][$ciz];
						next;
					}

					$line = Utils::DelSpace(<IN>);
					my ($kx, $ky, $kz, $dk, $e) = Utils::Split("\\s+", $line);
					if(!defined $e) {
						print("Error: Wrong format\n");
						exit;
					}
					$Energy[$ix][$iy][$iz] = $e;
					$emin = $e if($e < $emin);
					$emax = $e if($e > $emax);
				}
			}
		}
		next if($emin > $EFmax or $emax < $EFmin);

print "  To be recorded: e=($emin - $emax)\n";
		$BandEnergy[$nBand] = \@Energy;
		$BandIndex[$nBand]  = $iBand;
		$nBand++;
	}
	close(IN);
}

print("Save to [$OutFile]\n");
	open(OUT, ">$OutFile") or die "$!: Can not write to [$OutFile]\n";
	my ($nx1, $ny1, $nz1) = ($nx+1, $ny+1, $nz+1);
	print OUT <<EOT;
 BEGIN_INFO
   #
   # This is a Band-XCRYSDEN-Structure-File
   # Produced by ene2bxsf.pl
   #
   # Launch as: xcrysden --bxsf $OutFile
   #
   Fermi Energy: $EF
 END_INFO

 BEGIN_BLOCK_BANDGRID_3D
   $Title
   BEGIN_BANDGRID_3D_fermi
     $nBand
     $nx1 $ny1 $nz1
EOT
	printf(OUT "    %12.4g %12.4g %12.4g\n", 0.0, 0.0, 0.0);
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra11}, $PrimCrystal->{Ra12}, $PrimCrystal->{Ra13});
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra21}, $PrimCrystal->{Ra22}, $PrimCrystal->{Ra23});
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra31}, $PrimCrystal->{Ra32}, $PrimCrystal->{Ra33});

	for(my $ib = 0 ; $ib < $nBand ; $ib++) {
		printf(OUT "   BAND: %3d\n", $BandIndex[$ib]);
		my $ie = 0;
		for(my $ix = 0 ; $ix <= $nx ; $ix++) {
			for(my $iy = 0 ; $iy <= $ny ; $iy++) {
				for(my $iz = 0 ; $iz <= $nz ; $iz++) {
					printf(OUT "    %8.5f", $BandEnergy[$ib][$ix][$iy][$iz]);
					print(OUT "\n") if($OneByOne);
					$ie++;
					if($ie % 10 == 0) {
#					if($ie % 10 == 0 and $iz != $nz) {
						print(OUT "\n") if(!$OneByOne);
					}
				}
#				print(OUT "\n");
			}
			if($OneByOne) {
			}
			else {
				print(OUT "\n\n");
			}
			$ie = 0;
		}
	}
	print OUT <<EOT;
   END_BANDGRID_3D
 END_BLOCK_BANDGRID_3D
EOT

	close(OUT);

	exit;