#!/usr/bin/perl

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


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    = 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();