#!/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 Sci qw($hbar $me $e);# acos asin);
use JFile;
use Crystal::TightBinding;

#==============================================================
# General parameters
#==============================================================
my $CEk = $hbar * $hbar / 2.0 / $me * 1.0e20 / $e;

my $InputFile   = ($ARGV[0])? $ARGV[0] : 'TB.in';

#==============================================================
# Read files
#==============================================================
my $tb      = new TightBinding;
my $pParams = $tb->ReadInput($InputFile);
if($pParams->{Diagonalization} ne 'Math::MatrixReal' and !-f $pParams->{Diagonalization}) {
	print "Error: Can not find the diagonalization program [$pParams->{Diagonalization}].\n";
	exit;
}
if(!-f $pParams->{CIFFile}) {
	print "Error: Can not find the CIF file [$pParams->{CIFFile}].\n";
	exit;
}
if(!-f $pParams->{KPOINTSFile}) {
	print "Error: Can not find the KPOINTS file [$pParams->{KPOINTSFile}].\n";
	exit;
}

my ($n, $pKPoints) = $tb->ReadKPOINTS($pParams->{KPOINTSFile});
my $Crystal        = $tb->ReadCIF($pParams->{CIFFile});

#==============================================================
# Main routine
#==============================================================
print "#==========================================#\n";
print "# Free electron model                      #\n";
print "#==========================================#\n";
print "\n";

$tb->PrintInputInf();
my $krange = $tb->{pInput}{kRange};
my $me     = $tb->{pInput}{me};
my $UG     = $tb->{pInput}{UG};
$CEk = $CEk / $me;

print "# Crystal structure\n";
print "Crystal structure read from [$pParams->{CIFFile}]\n";
$tb->PrintCrystalInf();

#==============================================================
# Build Orbitals
#==============================================================
my @Orbitals;
my $c = 0;
print "# Build orbitals\n";
for(my $ikz = -$krange ; $ikz <= $krange ; $ikz++) {
for(my $iky = -$krange ; $iky <= $krange ; $iky++) {
for(my $ikx = -$krange ; $ikx <= $krange ; $ikx++) {
	$Orbitals[$c] = [ $ikx, $iky, $ikz, "$ikx $iky $ikz" ];
	$c++;
}
}
}
my $nOrb = $c;
print "nOrbitals: $nOrb\n";
for(my $io = 0 ; $io < $nOrb ; $io++) {
	printf "  %2d: ik=%2d %2d %2d\n", $io, $Orbitals[$io][0], $Orbitals[$io][1], $Orbitals[$io][2];
}
print "\n";
#exit;

#==============================================================
# Build k list
#==============================================================
my $pCalKList = $tb->BuildKList();
my @k         = @$pCalKList;
my $nCalK     = @k;
#exit;

print  "# k points\n";
print  "nK(target)=$n\n";
printf "nK(adjusted)=$nCalK\n";
for(my $i = 0 ; $i < $nCalK ; $i++) {
	my ($kx0, $ky0, $kz0, $name0) = @{$k[$i]};
	printf "k[$i]: (%8.4f, %8.4f, %8.4f, %s)\n", $kx0, $ky0, $kz0, $name0;
}
print "\n";
#exit;

my $out = JFile->new($pParams->{BandFile}, 'w') or die "$!: Can not write to [$pParams->{BandFile}].\n";
$out->print("name,kx,ky,kz,k");
for(my $io = 0 ; $io < $nOrb ; $io++) {
	$out->print(",E($Orbitals[$io][3])");
}
$out->print("\n");

for(my $i = 0 ; $i < @k ; $i++) {
	my $p = $k[$i];
	my ($kx, $ky, $kz, $ktot, $name) = @$p;

	print "#================================================\n";
	print "# k = ($kx, $ky, $kz): $name\n";
	print "#================================================\n";
	print "\n";

	$out->print("$name,$kx,$ky,$kz,$ktot");
	for(my $io = 0 ; $io < $nOrb ; $io++) {
		my ($ikx, $iky, $ikz, $name) = @{$Orbitals[$io]};

		my $k = $Crystal->GetReciprocalDistanceFromK($kx + $ikx, $ky + $iky, $kz + $ikz, 0, 0, 0);
		my $E = $CEk * $k * $k;
		$out->print(",$E");
	}
	$out->print("\n");
}
$out->Close();

exit;