#!/usr/bin/perl

use strict;
use lib 'd:/Programs/Perl/lib';
use lib "$ENV{TkPerlDir}/lib";
use Utils;
use JFile;
use CSV;
use Sci qw($RyToeV $e $me);

if(!$ARGV[0]) {
	print "Error: usage: perl ConvertBoltzTraPFiles.pl FileHeader\n\n";
	exit;
}
my $FileHeader = $ARGV[0];

my $EFWien2k = 0.57744; # Ry
my $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();
exit;