#!/usr/bin/perl

use strict;

use lib 'c:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';
use lib "$ENV{TkPerlDir}/lib";

use Deps;
use Utils;
use JFile;
use CSV;

use MyApplication;

use Sci qw($e $me $pi $h $hbar);
use Sci::Algorism;
use Crystal::VASP;

BEGIN {
#	$| = 1;
	open (STDERR, ">&STDOUT");
}

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#==========================================
# コマンドラインオプション読み込み
#==========================================
#$App->AddArgument("--Action", "--Action=[]", '');
$App->AddArgument("--EPSdir",  "--EPSdir=val (Def:0.1)", 0.1);
$App->AddArgument("--nSmooth", "--nSmooth=val (Def:11)", 11);
$App->AddArgument("--EF",      "--EF=val (Def:)",   "");
$App->AddArgument("--iVBM",    "--iVBM=val (Def:)", "");
exit 1 if($App->ReadArgs(1, "sjis", 0) != 1);
my $Args = $App->Args();

#===============================================
# 大域変数
#===============================================
#my $CARDir = "PbSnS2";
#$CARDir = $ARGV[0] if(defined $ARGV[0]);
my $CARDir    = $Args->GetGetArg(0);
my $EFShift   = $Args->GetGetArg('EF');
my $iVBMShift = $Args->GetGetArg('iVBM');
#$iVBMShift = 0 if(!defined $iVBMShift);
my $nSmooth   = $Args->GetGetArg('nSmooth');
my $EPSdir    = $Args->GetGetArg('EPSdir');
$EPSdir = 0.1 if(!defined $EPSdir);

my $vasp = new VASP;
my ($INCAR, $POSCAR, $KPOINTS, $OUTCAR);
my $BandDir = Utils::MakePath($CARDir, ["BandX"], '/', 0);
#print "BD[$BandDir]\n";
if(-d $BandDir) {
	$INCAR     = $vasp->GetINCARFileName($CARDir);
	$POSCAR    = Utils::MakePath($BandDir, ["POSCAR"], '/', 0); #$vasp->GetPOSCARFileName($CARDir);
	$KPOINTS   = Utils::MakePath($BandDir, ["KPOINTS"], '/', 0);
	$OUTCAR    = $vasp->GetOUTCARFileName($POSCAR);
}
else {
	$BandDir   = $CARDir;
	$INCAR     = $vasp->GetINCARFileName($CARDir);
	$POSCAR    = $vasp->GetPOSCARFileName($CARDir);
	$KPOINTS   = Utils::MakePath($CARDir, ["KPOINTS"], '/', 0);
	$OUTCAR    = Utils::MakePath($CARDir, ["OUTCAR"], '/', 0);
}
my $InputFile = Utils::MakePath($BandDir, ["EIGENVAL.band"], '/', 0); #$vasp->GetPOSCARFileName($CARDir);
   $InputFile = Utils::MakePath($BandDir, ["EIGENVAL"], '/', 0) if(!-f $InputFile);

my $EIGENVAL  = $InputFile; #$vasp->GetEIGENVALFileName($POSCAR);
my $OutputFile = "EffectiveMassBand.csv";

print "INCAR     : [$INCAR]\n";
print "POSCAR    : [$POSCAR]\n";
print "KPOINTS   : [$KPOINTS]\n";
print "OUTCAR    : [$OUTCAR]\n";
print "InputFile : [$InputFile]\n";
print "OutputFile: [$OutputFile]\n";
if(!-f $InputFile) {
	die "Error: Can not read [$InputFile].\n";
}
if(!-f $INCAR) {
	die "Error: Can not read [$INCAR].\n";
}
if(!-f $POSCAR) {
	die "Error: Can not read [$POSCAR].\n";
}
if(!-f $KPOINTS) {
	die "Error: Can not read [$KPOINTS].\n";
}
if(!-f $OUTCAR) {
	die "Error: Can not read [$OUTCAR].\n";
}

my $out = new JFile;
$out->Open($OutputFile, "w") or die "Error: Can not write to [$OutputFile].\n";

my $Crystal = $vasp->ReadStructureFromCARFiles($POSCAR, 1) or die "Error: Can not read [$POSCAR].\n";
my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
print "Latt: $a, $b, $c angstrom    $alpha, $beta, $gamma\n";
my $Volume = $Crystal->Volume();
print "Volume: $Volume A^3\n";

my $pINCARHash = $vasp->ReadINCARtoHash($INCAR);
my $ISPIN      = $pINCARHash->{ISPIN};
my $IsHF       = ($pINCARHash->{LHFCALC} =~ /T/i)? 1 : 0;
my $EF         = $vasp->ReadFermiEnergy($OUTCAR);
if(defined $EFShift) {
	if($EFShift =~ /^\+(.*)/) {
		$EF += $1;
	}
	else {
		$EF = $EFShift;
	}
}
print "ISPIN  : $ISPIN\n";
print "IsHF   : $IsHF\n";
print "nSmooth: $nSmooth\n";
print "EF     : $EF (EFShift: $EFShift)\n";

my $pKPOINTSHash = $vasp->ReadKPOINTStoHash($KPOINTS);
if($IsHF) {
}
else {
	print "Lattice mode in KPOINTS: $pKPOINTSHash->{LatticeMode}\n";
	if($pKPOINTSHash->{LatticeMode} !~ /^R/i) {
		print "Error: Lattice mode [$pKPOINTSHash->{LatticeMode}] is not supported.\n";
		exit;
	}
}

print "Read k points from [$KPOINTS].\n";
my ($nKPoints, $pKPOINTSLines, $pKPOINTSList, $iK0) = $vasp->ReadKListFromKPOINTS($KPOINTS, 0);
#my ($nKPoints, $pKPOINTSLines, $pKPOINTSList, $iK0) = $vasp->ReadKListFromKPOINTS($KPOINTS, 1);
print "  iK0(skip) = $iK0\n";
print "  nKPoints  = $nKPoints\n";
print "\n";

print "Read eigenvalues from [$EIGENVAL].\n";
$vasp->ReadEIGENVAL($EIGENVAL, undef, 0) or die "Error: Can not read [$EIGENVAL].\n";
#$vasp->ReadEIGENVAL($EIGENVAL, undef, 1) or die "Error: Can not read [$EIGENVAL].\n";
print "  iK0=$vasp->{iKPoint0}\n";
print "  nKPoints=$vasp->{iKPoint1}\n";

my $EnergyBandupArray = $vasp->{EnergyBandupArray};
my $nEnergyBand = $EnergyBandupArray->nEnergyBand();
print "nEnergyBand: $nEnergyBand\n";

my $iVBM;
my $nk;
for(my $il = 0 ; $il < $nEnergyBand ; $il++) {
	my $EB = $EnergyBandupArray->GetEnergyBand($il);
	$nk = $EB->nData();
	my $ik = 0;
	my ($kx, $ky, $kz, $dk, $e, $n) = $EB->GetData($ik);

	if(defined $EFShift) {
		if($e > $EF) {
			$iVBM = $il-1;
			last;
		}
	}
	else {
		if($e > 0.0) {
#		if($n < 0.01) {
			$iVBM = $il-1;
			last;
		}
	}
}
if(defined $iVBMShift) {
	if($iVBMShift =~ /^\+(.*)/) {
		$iVBM += $1;
	}
	elsif($iVBMShift ne '') {
		$iVBM = $iVBMShift 
	}
}
print "iVBM: $iVBM (shift: $iVBMShift)\n";

my $iCBM = $iVBM + 1;
#my $VBM3 = $EnergyBandupArray->GetEnergyBand($iVBM-4);
#my $VBM2 = $EnergyBandupArray->GetEnergyBand($iVBM-2);
my $VBM3 = $EnergyBandupArray->GetEnergyBand($iVBM-2);
my $VBM2 = $EnergyBandupArray->GetEnergyBand($iVBM-1);
my $VBM1 = $EnergyBandupArray->GetEnergyBand($iVBM);
my $CBM1 = $EnergyBandupArray->GetEnergyBand($iCBM);
my $CBM2 = $EnergyBandupArray->GetEnergyBand($iCBM+1);

my $pBoundaryDistances = $VBM1->{BoundaryDistances};
my $pBoundaryPositions = $VBM1->{BoundaryPositions};
print "k boundaries:\n";
for(my $i = 0 ; $i < @$pBoundaryDistances ; $i++) {
	my @b = Utils::Split("\\s+", $pBoundaryPositions->[$i]);
	printf "  $i: (%6.4f %6.4f %6.4f)  |k|=%12.4g\n", $b[0], $b[1], $b[2], $pBoundaryDistances->[$i];
}

$out->print("kx,ky,kz,k,EVBM-2(eV),EVBM-1(eV),EVBM(eV),ECBM(eV),ECBM+1(eV),"
		   ."d2EVBM(-2)/dE2,d2EVBM(-1)/dE2,d2EVBM/dE2,d2ECBM/dE2,d2ECBM(+1)/dE2,"
		   ."mh*(-2),mh*(-1),mh*,me*,me*(+1)");
my ($kx, $ky, $kz, $dk, $EVBM3, $EVBM2, $EVBM1, $ECBM1, $ECBM2, $n);
my $iBoundary = 1;
my $NextBoundary = $pBoundaryDistances->[$iBoundary];
my (@line, @d, @ev3, @ev2, @ev1, @ec1, @ec2);
my $ic = 0;
my ($d2EVBM3dk2, $d2EVBM2dk2, $d2EVBM1dk2, $d2ECBM1dk2, $d2ECBM2dk2);
my $Kk = 2.0 * $pi / 1.0e-10;
#for(my $ik = $iK0 ; $ik < $nk ; $ik++) {
for(my $ik = 0 ; $ik < $nk ; $ik++) {
#print "ik=$ik/$nk\n";
	($kx, $ky, $kz, $dk, $EVBM3, $n) = $VBM3->GetData($ik);
	($kx, $ky, $kz, $dk, $EVBM2, $n) = $VBM2->GetData($ik);
	($kx, $ky, $kz, $dk, $EVBM1, $n) = $VBM1->GetData($ik);
	($kx, $ky, $kz, $dk, $ECBM1, $n) = $CBM1->GetData($ik);
	($kx, $ky, $kz, $dk, $ECBM2, $n) = $CBM2->GetData($ik);
	$kx += 0.0;
	$ky += 0.0;
	$kz += 0.0;
#print "ic=$ic  ik=$ik: d[$ic-1]=", $d[$ic-1], "  dk=$dk   Next b=$NextBoundary\n";
	if($ik == 0 or $d[$ic-1] != $dk) {
		$line[$ic] = "$kx,$ky,$kz,$dk,$EVBM3,$EVBM2,$EVBM1,$ECBM1,$ECBM2";
		$d[$ic]    = $dk;
		$ev3[$ic]  = $EVBM3;
		$ev2[$ic]  = $EVBM2;
		$ev1[$ic]  = $EVBM1;
		$ec1[$ic]  = $ECBM1;
		$ec2[$ic]  = $ECBM2;
		$ic++;
	}

#	if($dk > $NextBoundary - 1.0e-3) {
	if($dk > $NextBoundary - 1.0e-5) {
	$out->print("\n");
		for(my $i = 0 ; $i < $ic ; $i++) {
#print "i=$i/$ic: $line[$i]\n";
		}
		for(my $i = 0 ; $i < $ic ; $i++) {
			$d2EVBM3dk2 = &Cald2Edk2($i, \@d, \@ev3);
			$d2EVBM2dk2 = &Cald2Edk2($i, \@d, \@ev2);
			$d2EVBM1dk2 = &Cald2Edk2($i, \@d, \@ev1);
			$d2ECBM1dk2 = &Cald2Edk2($i, \@d, \@ec1);
			$d2ECBM2dk2 = &Cald2Edk2($i, \@d, \@ec2);
			my $mhs3 = ($d2EVBM3dk2 == 0)? 0 : $hbar * $hbar / ($d2EVBM3dk2 * $e / $Kk / $Kk) / $me;
			my $mhs2 = ($d2EVBM2dk2 == 0)? 0 : $hbar * $hbar / ($d2EVBM2dk2 * $e / $Kk / $Kk) / $me;
			my $mhs1 = ($d2EVBM1dk2 == 0)? 0 : $hbar * $hbar / ($d2EVBM1dk2 * $e / $Kk / $Kk) / $me;
			my $mes1 = ($d2ECBM1dk2 == 0)? 0 : $hbar * $hbar / ($d2ECBM1dk2 * $e / $Kk / $Kk) / $me;
			my $mes2 = ($d2ECBM2dk2 == 0)? 0 : $hbar * $hbar / ($d2ECBM2dk2 * $e / $Kk / $Kk) / $me;
#print "ic=$i: $d[$i]  $ev1[$i]  $ec1[$i]  $d2EVBM1dk2  $d2ECBM1dk2\n";
#print "      m=$mhs1, $mes1\n";
			$out->print("$line[$i],$d2EVBM3dk2,$d2EVBM2dk2,$d2EVBM1dk2,$d2ECBM1dk2,$d2ECBM2dk2,$mhs3,$mhs2,$mhs1,$mes1,$mes2\n");
		}

		@d   = ();
		@ev3 = ();
		@ev2 = ();
		@ev1 = ();
		@ec1 = ();
		@ec2 = ();
		$ic  = 0;
		$line[$ic] = "$kx,$ky,$kz,$dk,$EVBM3,$EVBM2,$EVBM1,$ECBM1,$ECBM2";
		$d[$ic]  = $dk;
		$ev3[$ic] = $EVBM3;
		$ev2[$ic] = $EVBM2;
		$ev1[$ic] = $EVBM1;
		$ec1[$ic] = $ECBM1;
		$ec2[$ic] = $ECBM2;
		$ic++;
		$NextBoundary = $pBoundaryDistances->[++$iBoundary];
	}
}
$out->Close();

exit;

sub Cald2Edk2
{
	my ($i, $pd, $pe) = @_;

	my $n = @$pe;
	return undef if($n <= 2);

	if($i == 0) {
#		my $dedk1 = ($pe->[1] - $pe->[0])/($pd->[1] - $pd->[0]);
#		my $dedk2 = ($pe->[0] - $pe->[1])/($pd->[1] - $pd->[0]);
		my $dedk1 = ($pe->[2] - $pe->[1])/($pd->[2] - $pd->[1]);
		my $dedk2 = ($pe->[1] - $pe->[0])/($pd->[1] - $pd->[0]);
		return ($dedk1 - $dedk2) / ($pd->[2] - $pd->[0]) * 2;
	}
	elsif($i == $n-1) {
#		my $dedk1 = ($pe->[$i]   - $pe->[$i-1])/($pd->[$i] - $pd->[$i-1]);
#		my $dedk2 = ($pe->[$i-1] - $pe->[$i])  /($pd->[$i] - $pd->[$i-1]);
		my $dedk1 = ($pe->[$i]   - $pe->[$i-1])/($pd->[$i]   - $pd->[$i-1]);
		my $dedk2 = ($pe->[$i-1] - $pe->[$i-2])/($pd->[$i-1] - $pd->[$i-2]);
		return ($dedk1 - $dedk2) / ($pd->[$i] - $pd->[$i-2]) * 2;
	}
	my $dedk1 = ($pe->[$i+1] - $pe->[$i])/($pd->[$i+1] - $pd->[$i]);
	my $dedk2 = ($pe->[$i] - $pe->[$i-1])/($pd->[$i] - $pd->[$i-1]);
	return ($dedk1 - $dedk2) / ($pd->[$i+1] - $pd->[$i-1]) * 2;
}

# not correct
sub Cald2Edk2_notcorrect
{
	my ($i, $pd, $pe) = @_;

	my $n = @$pe;
#print"*n=$n\n";
	return undef if($n <= 2);

	my (@x, @y);
	if($i == 0) {
#		my $nn = 3;
		my $nn = ($n == 3)? 3 : 4;
		for(my $j = 0 ; $j < $nn ; $j++) {
			$x[$j] = $pd->[$j];
			$y[$j] = $pe->[$j];
		}
#print "n(1st)=", scalar @x, "\n";
	}
	elsif($i == $n-1) {
#		my $nn = 3;
		my $nn = ($n == 3)? 3 : 4;
		for(my $j = 0 ; $j < $nn ; $j++) {
			my $idx = $n-1 - ($nn-1) + $j;
#print "j=$j  idx=$idx\n";
			$x[$j] = $pd->[$idx];
			$y[$j] = $pe->[$idx];
		}
#print "n(last)=", scalar @x, "\n";
	}
	else {
		my $nn = 3;
		for(my $j = 0 ; $j < $nn ; $j++) {
			my $idx = $i - 1 + $j;
#print "j=$j  idx=$idx\n";
			$x[$j] = $pd->[$idx];
			$y[$j] = $pe->[$idx];
		}
#print "n($i)=", scalar @x, "\n";
	}

	my @f1;
	for(my $j = 0 ; $j < @x ; $j++) {
		$f1[$j] = Algorism::DifferentiateByLagrange(\@x, \@y, $x[$j]);
	}
	my $d2Edk2 = Algorism::DifferentiateByLagrange(\@x, \@f1, $pd->[$i]);

#	for(my $j = 0 ; $j < @$pd ; $j++) {
#		$f1[$j] = Algorism::DifferentiateByCubicSpline($pd, $pe, $pd->[$j]);
#	}
#	my $d2Edk2 = Algorism::DifferentiateByCubicSpline($pd, \@f1, $pd->[$i]);

	return $d2Edk2;
}