#!/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 MyApplication;
use Sci qw($RyToeV $e $me);
use Crystal::WIEN2k;
use Crystal::XCrySDen;

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

#===============================================
# スクリプト大域変数
#===============================================
my $BoltzTraPPath        = "$ENV{HOME}/boltztrap/SRC_BoltzTraP110/x_trans";
my ($EF, $DeltaE, $Ecut) = (0.0, 0.0005, 0.4); # Ry
my $ne                   = 240;     # Number of Valence Electrons
my $lpfac                = 5;       # number of lattice points per k-point
my $efcut                = 0.15;    # Energy range of chemical potential arround EF used for integration
my ($Tmax, $Tgrid)       = (800.0, 50.0);
my $DOSMethod            = "TETRA"; # Scheme to obtain DOS. HISTO|TETRA

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",    "--Action=[all|ExtractData1|ExtractData2]", '');
$App->AddArgument("--DeltaE",    "--DeltaE=val [Ry]", '');
$App->AddArgument("--Ecut",      "--Ecut=val [Ry]", '');
$App->AddArgument("--ne",        "--ne=val: # of valence electrons", '');
$App->AddArgument("--lpfac",     "--lpfac=integer: # of lattice points per k-point", '');
$App->AddArgument("--efcut",     "--efcut=val: Energy range of chem. pot. arround EF for integration", '');
$App->AddArgument("--Tmax",      "--Tmax=val [K]", '');
$App->AddArgument("--Tgrid",     "--Tgrid=val [K]", '');
$App->AddArgument("--DOSMethod", "--DOSMethod=[TETRA|HISTO]", '');
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my %Args = $App->Args()->GetArgHash();
$App->SetDebug($Args{DebugMode}) if(defined $Args{DebugMode});

$DeltaE    = $Args{DeltaE}     if(defined $Args{DeltaE});
$Ecut      = $Args{Ecut}       if(defined $Args{EcutD});
$ne        = $Args{ne}         if(defined $Args{ne});
$lpfac     = $Args{lpfac}      if(defined $Args{lpfac});
$efcut     = $Args{efcut}      if(defined $Args{efcut});
$Tmax      = $Args{Tmax}       if(defined $Args{Tmax});
$Tgrid     = $Args{Tgrid}      if(defined $Args{Tgrid});
$DOSMethod = $Args{DOSMethod}  if(defined $Args{DOSMethod});

#==========================================
# メイン関数スタート
#==========================================
my $Action  = $Args{Action};
my $DataDir = $App->Args()->GetGetArg(0);
$DataDir = "." if($DataDir eq '');
$DataDir = Utils::ReduceDirectory($DataDir);
my $StructFile = &GetStructPath($DataDir);
if(!-e $StructFile) {
	$App->print("Error: StructFile [$StructFile] could not obtained.\n");
	exit 0;
}
my $Title        = &GetTitle($StructFile);
my $IntransPath  = Deps::ReplaceExtension($StructFile, ".intrans");
my $PrefPath = Deps::MakePath($DataDir, ".prefs-scf");

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

my $config = new IniFile();
$config->ReadAll($PrefPath);

$App->print("Struct: $StructFile\n");
$App->print("Lattice System: ", $Crystal->LatticeSystem(), "\n");
$App->print("Space Group   : ", $Crystal->SPGName(), "\n");
if($config->{spin} eq '') {
	$App->print("Spin polarized: no\n");
}
else {
	$App->print("Spin polarized: yes\n");
}

my $ret = 0;
if($Action eq 'all') {
	if($config->{spin} eq '') {
		&Execute($Crystal, $config, '', $Title, $IntransPath);
	}
	else {
		&Execute($Crystal, $config, 'up', $Title, $IntransPath);
		&Execute($Crystal, $config, 'dn', $Title, $IntransPath);
	}
}
elsif($Action eq 'ExtractData1') {
	my $EFWIEN2k = $WIEN2k->ReadFermiEnergy($StructFile);
	&ExtractData1($EFWIEN2k, $Crystal);
}
elsif($Action eq 'ExtractData2') {
	my $EFWIEN2k = $WIEN2k->ReadFermiEnergy($StructFile);
	&ExtractData2($EFWIEN2k, $Crystal);
}
else {
	$App->print("Error: Invald Action: [$Action]\n");
	$App->Usage();
}

exit $ret;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================
sub GetTitle
{
	my ($path) = @_;

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($path);
	my $Title   = $filebody;
	$Title =~ s/\s/_/g;
	return $Title;
}

sub GetStructPath
{
	my ($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 from Data Directory [$DataDir].\n");
			return undef;
		}
	}
	return $StructFile;
}

sub Execute
{
	my ($Crystal, $config, $spin, $Title,$IntransPath) = @_;

	$App->print("\nCalculate Transport Properties by BoltzTraP for spin=[$spin]\n");
	$App->print("BoltzTraP Path: $BoltzTraPPath\n");
	$App->print(".struct Path  : $StructFile\n");
	$App->print(".intrans Path : $IntransPath\n");
	$EF = $WIEN2k->ReadFermiEnergy($StructFile);
	$App->print("EF: $EF Ry\n");

	open(INT, ">$IntransPath") or die "$!: Can not write to [$IntransPath].\n";
	print INT<<EOT;
WIEN                      # Format of DOS
0 1 0 0.0                 # iskip (not presently used) idebug setgap shiftgap
$EF $DeltaE $Ecut $ne  # Fermilevel (Ry), energygrid, energy span around Fermilevel, number of electrons
CALC                      # CALC (calculate expansion coeff), NOCALC read from file
$lpfac                         # lpfac, number of latt-points per k-point
BOLTZ                     # run mode (only BOLTZ is supported)
$efcut                     # (efcut) energy range of chemical potential
$Tmax $Tgrid                # Tmax, temperature grid
-1.                        # energyrange of bands given individual DOS output sig_xxx and dos_xxx (xxx is band number)
$DOSMethod
EOT
	close(INT);

	my $opt = "";
	if($spin ne '') {
		$opt .= " -$spin";
	}
	if($config->{so} eq 'on') {
		$opt .= " -so";
	}
	$opt = "BoltzTraP" if($opt eq '');
	my $cmd = "$BoltzTraPPath BoltzTraP $opt";
	if(Utils::Execute($cmd)) {
		$App->print("Error: Can not execute [$cmd].\n");
		exit 0;
	}
	if($spin ne '') {
		foreach my $key qw(condtens halltens outputtrans sigxx sigxxx transdos engre) {
			my $OldPath = Deps::ReplaceExtension($StructFile, ".$key");
			my $NewPath = Deps::ReplaceExtension($StructFile, ".${key}$spin");
			Utils::Execute("mv $OldPath $NewPath");
		}
	}
}

sub ExtractData1
{
	my ($EFWien2k, $Crystal) = @_;

my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($DataDir);
my $FileHeader = $filename;#"LaCuOSe";
my $Volume = $Crystal->Volume(); #4.067 * 4.067 * 8.7975 * 1.0e-24; # cm3
my $tau    = 2.7e-15; #1.5e-15; #s

my $T = 300;
my ($EFStart, $EFEnd, $EFStep) = (-1.0, 2.0, 0.02);
my $EFOutCSV = "${T}K.csv";

my $EF = -0.1933;
my $TOutCSV = "${EF}eV.csv";

my $FileFormat = 'Trace';
#my $FileFormat = 'CondTens';
#my $FileFormat = 'HallTens';

my @TraceFileFields;
my $InFile;

if($FileFormat eq 'Trace') {
	$InFile = "$FileHeader.trace";
	@TraceFileFields = ("EF(eV)", "T(K)", "N(cm-3)", "D(EF)",  "S(uV/K)", "sigma(S/cm)", "RH", "kappa0", "c", "chi");
}
elsif($FileFormat eq 'CondTens') {
	$InFile = "$FileHeader.condtens";
	@TraceFileFields = ("EF(eV)", "T(K)", "N(cm-3)", 
						"sigma_xx", "sigma_xy", "sigma_xz", "sigma_yx", "sigma_yy", "sigma_yz", "sigma_zx", "sigma_zy", "sigma_zz", 
						"S_xx", "S_xy", "S_xz", "S_yx", "S_yy", "S_yz", "S_zx", "S_zy", "S_zz",
						"k_xx", "k_xy", "k_xz", "k_yx", "k_yy", "k_yz", "k_zx", "k_zy", "k_zz");
}
elsif($FileFormat eq 'HallTens') {
	$InFile = "$FileHeader.halltens";
	@TraceFileFields = ("EF(eV)", "T(K)", "N(cm-3)");
	my $c = 3;
	for(my $i = 1 ; $i <= 3 ; $i++) {
		for(my $j = 1 ; $j <= 3 ; $j++) {
			for(my $k = 1 ; $k <= 3 ; $k++) {
				$TraceFileFields[$c++] = "R$i$j$k";
			}
		}
	}
}
else {
	print("Invalid File Format [$FileFormat]\n");
	exit;
}


my $in = new JFile;
$in->Open($InFile, "r") or die "Error: can not read [$InFile].\n";
my $outEF = new JFile;
$outEF->Open($EFOutCSV, "w") or die "Error: can not write to [$EFOutCSV].\n";
my $outT = new JFile;
$outT->Open($TOutCSV, "w") or die "Error: can not write to [$TOutCSV].\n";

my $line = $in->ReadLine();
$outEF->print(join(',', @TraceFileFields));
$outEF->print("\n");
$outT->print(join(',', @TraceFileFields));
$outT->print("\n");

my @lines;
my $c = 0;
my $EFatDmin = -1e10;
my $DEFmin   = 1.0e100;
while(1) {
	$lines[$c] = $in->ReadLine();
	last if(!defined $lines[$c]);
	Utils::DelSpace($lines[$c]);
	last if($lines[$c] eq '');

	if($FileFormat eq 'Trace') {
		my ($EFi, $Ti, $Ni, @others) = Utils::Split("\\s+", $lines[$c]);
		$EFi     = ($EFi - $EFWien2k) * $RyToeV;
#print "EF=$EFi\n";
		if($EFi < 0) {
			$c++;
			next;
		}

		my ($DEFi, $Si, $sigmai, $RHi, $kappa0i, $ci, $chii) = @others;
#print "D(EF)=$DEFi\n";
		if($DEFi < $DEFmin) {
			$DEFmin = $DEFi;
			$EFatDmin = $EFi;
		}
	}

	$c++;
}
my $nData = $c;
#DEFはRy単位で２倍になっているので、その分を修正している
$DEFmin /= $Volume * $RyToeV * 2; # n(u) (e/uc/eV)
print "nData=$nData\n";
print "EF at minimum D(EF): $EFatDmin eV, D(EF)=$DEFmin e/cm3/Ry\n";
#exit;

my $EFnow = $EFStart;
for(my $i = 0 ; $i < $nData ; $i++) {
	$line = $lines[$i];
	Utils::DelSpace($line);

	my ($EFi, $Ti, $Ni, @others) = Utils::Split("\\s+", $line);
	$EFi     = ($EFi - $EFWien2k) * $RyToeV;
	$Ni     /= $Volume; # N (e/uc)

	if($FileFormat eq 'Trace') {
		my ($DEFi, $Si, $sigmai, $RHi, $kappa0i, $ci, $chii) = @others;
#DEFはRy単位で２倍になっているので、その分を修正している
		$DEFi   /= $Volume * $RyToeV * 2.0; # n(u) (e/uc/eV)
		$Si     *= 1.0e6;   # S (uV/K)
		$sigmai *= $tau * 1.0e-2;  # sigma (S/cm)
		if($EFnow <= $EFi and $EFnow <= $EFEnd and abs($Ti - $T) < 1.0e-3) {
			print "$Ti K: EF = $EFi eV\n";
			$outEF->print(join(',', $EFi, $Ti, $Ni, $DEFi, $Si, $sigmai, $RHi, $kappa0i, $ci, $chii));
			$outEF->print("\n");
			while($EFnow <= $EFi) {
				$EFnow += $EFStep;
			}
		}
		if(abs($EFi - $EF) < 1.0e-3) {
			print "$Ti K: EF = $EFi eV\n";
			$outT->print(join(',', $EFi, $Ti, $Ni, $DEFi, $Si, $sigmai, $RHi, $kappa0i, $ci, $chii));
			$outT->print("\n");
		}
#print "EF now: $EFnow eV\n";
#exit;
	}
	elsif($FileFormat eq 'CondTens') {
		for(my $i = 0 ; $i < 9 ; $i++) {
			$others[$i]   *= $tau * 1.0e-2; # sigma (S/cm)
			$others[$i+9] *= 1.0e6;         # S (uV/K)
		}
		if($EFnow <= $EFi and $EFnow <= $EFEnd and abs($Ti - $T) < 1.0e-3) {
			print "$Ti K: EF = $EFi eV\n";
			$outEF->print(join(',', $EFi, $Ti, $Ni, @others));
			$outEF->print("\n");
			$EFnow += $EFStep;
		}
		if(abs($EFi - $EF) < 1.0e-3) {
			print "$Ti K: EF = $EFi eV\n";
			$outT->print(join(',', $EFi, $Ti, $Ni, @others));
			$outT->print("\n");
		}
	}
	elsif($FileFormat eq 'HallTens') {
		if($EFnow <= $EFi and $EFnow <= $EFEnd and abs($Ti - $T) < 1.0e-3) {
			print "$Ti K: EF = $EFi eV\n";
			$outEF->print(join(',', $EFi, $Ti, $Ni, @others));
			$outEF->print("\n");
			$EFnow += $EFStep;
		}
		if(abs($EFi - $EF) < 1.0e-3) {
			print "$Ti K: EF = $EFi eV\n";
			$outT->print(join(',', $EFi, $Ti, $Ni, @others));
			$outT->print("\n");
		}
	}
}

$outT->Close();
$outEF->Close();
$in->Close();

}

sub ExtractData2
{
	my ($EFWien2k, $Crystal) = @_;

my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($DataDir);
my $FileHeader = $filename;#"LaCuOSe";

my $Volume = $Crystal->Volume(); #4.067 * 4.067 * 8.7975 * 1.0e-24; # cm3
my $tau    = 3.4e-15; #1.5e-15; #s

my $T = 300;
my ($EFStart, $EFEnd, $EFStep) = (-1.0, 2.0, 0.02);
my $EF   = -0.1933;
my $EVBM = 0.0;

my $EFOutCSV = "${T}K-All.csv";
my $TOutCSV  = "${EF}eV-All.csv";

my $TraceFile = "$FileHeader.trace";
my @TraceFileFields = ("EF(eV)", "T(K)", "N(cm-3)", "D(EF)",  "S(uV/K)", "sigma(S/cm)", "RH", "kappa0", "c", "chi");
my $CondTensFile = "$FileHeader.condtens";
my @CondTensFileFields = ("EF(eV)", "T(K)", "N(cm-3)", 
						"sigma_xx", "sigma_xy", "sigma_xz", "sigma_yx", "sigma_yy", "sigma_yz", "sigma_zx", "sigma_zy", "sigma_zz", 
						"S_xx", "S_xy", "S_xz", "S_yx", "S_yy", "S_yz", "S_zx", "S_zy", "S_zz",
						"k_xx", "k_xy", "k_xz", "k_yx", "k_yy", "k_yz", "k_zx", "k_zy", "k_zz");
my $HallTensFile = "$FileHeader.halltens";
my @HalLTensFileFields = ("EF(eV)", "T(K)", "N(cm-3)");
my $c = 3;
my @xyz = qw(x y z);
for(my $i = 0 ; $i < 3 ; $i++) {
	for(my $j = 0 ; $j < 3 ; $j++) {
		for(my $k = 0 ; $k < 3 ; $k++) {
			$HalLTensFileFields[$c++] = "RH_$xyz[$i]$xyz[$j]$xyz[$k]";
		}
	}
}

my $TraceCSV = new CSV;
$TraceCSV->ReadFreeFormat($TraceFile, "\\s+", \@TraceFileFields, 0);
my $pEF0 = $TraceCSV->GetXData("EF.*");
my $pT0  = $TraceCSV->GetXData("T.*");
my $pN0  = $TraceCSV->GetXData("N.*");
my $pDEF = $TraceCSV->GetXData("D\\(EF\\).*");
#print "p=$pEF0,$pT0,$pN0,$pDEF\n";

my $CondTensCSV = new CSV;
$CondTensCSV->ReadFreeFormat($CondTensFile, "\\s+", \@CondTensFileFields, 0);
my $psxx = $CondTensCSV->GetXData("sigma_xx.*");
my $pszz = $CondTensCSV->GetXData("sigma_zz.*");
my $pSxx = $CondTensCSV->GetXData("S_xx.*");
my $pSzz = $CondTensCSV->GetXData("S_zz.*");
#print "p=$psxx,$pszz,$pSxx,$pSzz\n";

my $HallTensCSV = new CSV;
$HallTensCSV->ReadFreeFormat($HallTensFile, "\\s+", \@HalLTensFileFields, 0);
#my @pData = $HallTensCSV->LabelArray();
#foreach my $s (@pData) {
#	print "s=$s, ", join(', ', @$s), "\n";
#}
my $pRHyxz = $HallTensCSV->GetXData("RH_yxz.*");
my $pRHzyx = $HallTensCSV->GetXData("RH_zyx.*");
#print "p=$pRHyxz, $pRHzyx\n";
my $nData = @$pEF0;
print "nData=$nData\n";

my $outEF = new JFile;
$outEF->Open($EFOutCSV, "w") or die "Error: can not write to [$EFOutCSV].\n";
my $outT = new JFile;
$outT->Open($TOutCSV, "w") or die "Error: can not write to [$TOutCSV].\n";

my @Nint;
my @OutputFields = ("EF(eV)", "T(K)", "N(cm-3)", "Nint(cm-3)", "DEF(cm-3eV-1)", "sigma_xx(S/cm)", "sigma_zz(S/cm)", 
					"RH_yxz(cm3/C)", "RH_yzx(cm3/C)", "nHall_xx(cm-3)", "nHall_zz(cm-3)", "uHall_xx(cm2/Vs)", "uHall_zz(cm2/Vs)", 
					"m_xx(me)", "m_zz(me)", "S_xx(uV/K)", "S_zz(uV/K)");
$outEF->print(join(',', @OutputFields), "\n");
$outT->print(join(',', @OutputFields), "\n");

my $EFnow = $EFStart;
my $iDEFMin = 0;
my $DEFMin = 1.0e100;
for(my $i = 0 ; $i < $nData ; $i++) {
	my ($EFi, $Ti, $DEFi) = ($pEF0->[$i], $pT0->[$i], $pDEF->[$i]);
	$EFi = ($EFi - $EFWien2k) * $RyToeV;
	if($EVBM <= $EFi and $EFnow <= $EFi and $EFnow <= $EFEnd and abs($Ti - $T) < 1.0e-3) {
		if($DEFi < $DEFMin) {
			$DEFMin = $DEFi;
			$iDEFMin = $i;
		}
		while($EFnow <= $EFi) {
			$EFnow += $EFStep;
		}
	}
}
print "D(EF) min = $pDEF->[$iDEFMin] at i=$iDEFMin, E=$pEF0->[$iDEFMin] Ry\n";
if($iDEFMin > 0) {
	$Nint[$iDEFMin] = 0;
	for(my $i = $iDEFMin-1 ; $i >= 0 ; $i--) {
		$Nint[$i] = $Nint[$i+1] + ($pDEF->[$i+1] + $pDEF->[$i]) / 2.0 * ($pEF0->[$i+1] - $pEF0->[$i]) / $Volume / 2.0;
	}
	for(my $i = $iDEFMin+1 ; $i < $nData ; $i++) {
		$Nint[$i] = $Nint[$i-1] - ($pDEF->[$i] + $pDEF->[$i-1]) / 2.0 * ($pEF0->[$i] - $pEF0->[$i-1]) / $Volume / 2.0;
	}
}

my $EFnow = $EFStart;
for(my $i = 0 ; $i < $nData ; $i++) {
	my ($EFi, $Ti, $Ni) = ($pEF0->[$i], $pT0->[$i], $pN0->[$i]);
	my ($DEFi, $sxx, $szz, $Sxx, $Szz) = ($pDEF->[$i], $psxx->[$i], $pszz->[$i], $pSxx->[$i], $pSzz->[$i]);
	my ($RHyxz, $RHzyx) = ($pRHyxz->[$i], $pRHzyx->[$i]);

if(1) {
	$EFi = ($EFi - $EFWien2k) * $RyToeV;
	$Ni /= $Volume; # N (e/uc)
	$DEFi /= $Volume * $RyToeV * 2.0; # (e/cm3/eV) #DEFはRy単位で２倍になっているので、その分を修正している
	$sxx  *= $tau * 1.0e-2;  # sigma (S/cm)
	$szz  *= $tau * 1.0e-2;  # sigma (S/cm)
	$Sxx  *= 1.0e6;   # S (uV/K)
	$Szz  *= 1.0e6;   # S (uV/K)
	$RHyxz *= 1.0e6; # cm3/C
	$RHzyx *= 1.0e6; # cm3/C
}
	my $nHallxx = ($RHyxz == 0.0)? 1.0e100 : 1.0 / $e / $RHyxz;
	my $nHallzz = ($RHzyx == 0.0)? 1.0e100 : 1.0 / $e / $RHzyx;
	my $uHallxx = ($nHallxx == 0.0)? 1.0e100 : $sxx / $e / $nHallxx;
	my $uHallzz = ($nHallzz == 0.0)? 1.0e100 : $szz / $e / $nHallzz;
	my $mxx     = ($uHallxx == 0.0)? 1.0e100 : $e * $tau / ($uHallxx * 1.0e-4) / $me;
	my $mzz     = ($uHallzz == 0.0)? 1.0e100 : $e * $tau / ($uHallzz * 1.0e-4) / $me;

	if($EFnow <= $EFi and $EFnow <= $EFEnd and abs($Ti - $T) < 1.0e-3) {
		print "$Ti K: EF = $EFi eV\n";
		$outEF->print(join(',', $EFi, $Ti, $Ni, $Nint[$i], $DEFi, $sxx, $szz, $RHyxz, $RHzyx, $nHallxx, $nHallzz, $uHallxx, $uHallzz, $mxx, $mzz, $Sxx, $Szz,), "\n");
		while($EFnow <= $EFi) {
			$EFnow += $EFStep;
		}
	}
	if(abs($EFi - $EF) < 1.0e-3) {
		print "$Ti K: EF = $EFi eV\n";
		$outT->print(join(',', $EFi, $Ti, $Ni, $DEFi, $sxx, $szz, $RHyxz, $RHzyx, $nHallxx, $nHallzz, $uHallxx, $uHallzz, $mxx, $mzz, $Sxx, $Szz,), "\n");
	}
}

$outT->Close();
$outEF->Close();
}