#!/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;
use Crystal::VASP;
use Sci::Algorism;

my ($x, $y, $z) = (0.083333, 0.083333, 0.663406);
my $OutputCSV = ($ARGV[0])? $ARGV[0] : "LDOS.csv";
my ($Emin, $Emax) = (defined $ARGV[2])? ($ARGV[1], $ARGV[2]) : (-1.5, 2);

my $OUTCAR  = "../SCF/OUTCAR";
my $EIGENVAL = "EIGENVAL";

print "Output: [$OutputCSV]\n";
print "E range: $Emin - $Emax eV\n";

my $out = new JFile($OutputCSV, "w") or die "Error: $!: Can not write to [$OutputCSV].\n";
$out->print("i,iLevel,iK,weight(K),Energy(eV),Density,Density(intg),MinV,MaxV,S,S2\n");

my $vasp = new VASP;

my $EF = $vasp->ReadFermiEnergy($OUTCAR);
print "EF = $EF eV\n";

my ($nn, $nK, $nLevels, $pEnergyBands, $pKPointInf) = $vasp->ReadEIGENVALtoArray($EIGENVAL, $EF);
print "nK = $nK   nLevels = $nLevels\n";

my ($LevelMin, $LevelMax) = (-1, -1);
for(my $i = 0 ; $i < @$pEnergyBands ; $i++) {
	for(my $iK = 0 ; $iK < $nK ; $iK++) {
		for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
			my $eup = $pEnergyBands->[$iL][$iK][0];
			my $edn = $pEnergyBands->[$iL][$iK][1];
#print "$i,$iK,$iL: $eup eV\n";

			if(defined $eup and $Emin <= $eup and $eup < $Emax) {
#print "$i,$iK,$iL: $eup eV\n";
				if($LevelMin == -1 or $iL < $LevelMin) {
					$LevelMin = $iL;
				}
				if($LevelMax == -1 or $LevelMax < $iL) {
#print "up: $i,$iK,$iL: $eup eV\n";
					$LevelMax = $iL;
				}
			}
			if(defined $edn and $Emin <= $edn and $edn < $Emax) {
				if($LevelMin == -1 or $iL < $LevelMin) {
					$LevelMin = $iL;
				}
				if($LevelMax == -1 or $LevelMax < $iL) {
#print "dn: $i,$iK,$iL: $eup eV\n";
					$LevelMax = $iL;
				}
			}
		}
	}
}
print "iLevel range: $LevelMin - $LevelMax\n";

my $dx = 16 * 0.3;
my $dy = 16 * 0.3;
my $dz = 22 * 0.3;
my $count = 0;
for(my $iL = $LevelMin ; $iL <= $LevelMax ; $iL++) {
	my $iL1 = $iL+1;
	for(my $iK = 0 ; $iK < $nK ; $iK++) {
		my $iK1 = $iK+1;
		my $eup = $pEnergyBands->[$iL][$iK][0];
		my $edn = $pEnergyBands->[$iL][$iK][1];

		my $fmask = sprintf("PARCHG.%04d.%04d*", $iL1, $iK1);
		my @f = glob($fmask);
		my $path = $f[0];

		next if(!defined $path or !-f $path);
		my $pDensity = $vasp->ReadDensityFileToHash($path);
		my $pValues = $pDensity->{pValues};
		my $nx = $pDensity->{nx};
		my $ny = $pDensity->{ny};
		my $nz = $pDensity->{nz};
		my $ix = int($x * $nx);
		my $iy = int($y * $ny);
		my $iz = int($z * $nz);
my ($MinV, $MaxV) = (1.0e100, -1.0e100);
my $S   = 0.0;
my $S2 = 0.0;
my $V   = 0.0;
for(my $jx = 0 ; $jx < $nx ; $jx++) {
	for(my $jy = 0 ; $jy < $ny ; $jy++) {
		for(my $jz = 0 ; $jz < $nz ; $jz++) {
			my $val = $pValues->[$jx][$jy][$jz];
			$S += $val;
			$S2 += $val*$val;
			$MinV = $val if($MinV > $val);
			$MaxV = $val if($MaxV < $val);
			if($ix-$dx <= $jx and $jx <= $ix+$dx and
			   $iy-$dy <= $jy and $jy <= $iy+$dy and
			   $iz-$dz <= $jz and $jz <= $iz+$dz) {
			   	$V += $val;
			}
		}
	}
}
$S   /= $nx*$ny*$nz;
$S2 /= $nx*$ny*$nz;
printf "$count: $path: val=($MinV, $MaxV)  S=%3.2f  S2=%6.2f\n", $S, $S2;

		my $v = Algorism::InterpolateValuesInCube(
					$nx, $ny, $nz, $pValues, $ix, $iy, $iz);
		my $pKInfArray = $pKPointInf->[$iK];

		$count++;
printf "  w=%6.4f ($ix,$iy,$iz): %6.3f eV: v=%6.4f  V=%6.4f\n", $pKInfArray->[3], $eup, $v, $V;
$out->print("$count,$iL1,$iK1,$pKInfArray->[3],$eup,$v,$V,$MinV,$MaxV,$S,$S2\n");

		delete $pDensity->{pValues};
		undef $pDensity;
	}
}

$out->Close();
exit;