#!/usr/bin/perl -w

use lib 'd:/git/tkProg/tklib/Perl/lib';
use lib '/home/share/tkProg/tklib/Perl/lib';

use strict;

use Jcode;
#use encoding 'sjis';

use Sci qw($pi $kB $e $e0 $me $h $hbar $KToeV erfc);
use Sci::Material;
use Sci::Optimize;

my $material = new Material;
my $opt      = new Optimize;

# EF計算収束条件
my $EFEPS    = 1.0e-6;
my $EFSEPS   = 1.0e-12;
my $nMaxIter = 1000;
my $IsPrint  = 0;

# 温度依存性計算範囲
my $T      = 300.0;
my $Tstart =  40.0;
my $Tend   = 320.1;
my $Tstep  =  40.0;

# バンド端構造
my $Effme = 0.35;
my $Effmh = 1.5;
my $Ec    = 3.0;
my $Ev    = 0.0;
my $CBMultiplicity = 1.0;
my $VBMultiplicity = 1.0;

# ドーパント
my $ND = 1.9e19;
my $ED = $Ec + 1.000; # eV
my $NA = 0.0e16;
my $EA = $Ev + 0.045; # eV

# 伝導帯裾状態
my $Eu    = 0.01;   # eV
my $Ntail = 4.0e19; # cm-3
my $Nt0   = 0.0; # $Ntail / $Eu;
my $Eth   = $Ec + 1;
my $Ecmin = $Ec - 0.2; # Minimum energy for CB integration
#printf "Ntail=%10.4g Nt0=%10.4g cm-3 E0=%6.3f eV\n", $Ntail, $Nt0, $Eu;

# 散乱
my $Em   = $Ec;     # Mobility edge
my $tau0 = 1.0e-12;
$tau0 = 2.0 * 14.2 / 13.1 * $tau0;
my $r    = 0.5; # 散乱因子

# パーコレーション
my $E0      = 0.05;      # エネルギー幅 (Gaussian: exp( -[($E-$Ecenter)/$E0]^2 ))
my $Ecenter = $Ec + 0.1; # エネルギー中央値
if(0) {
	for(my $x = -3 ; $x <= 3 ; $x += 0.1) {
		my $P = &BarrierProbability($x, 0.0, 1.0);
		print "$x: $P\n";
	}
	exit;
}

# parameters: $tau0, $ND, $Ecenter, $E0
my $iSample = 5;
# For Thermally activated
if($iSample == -1) {
	$ED      = $Ec - 0.080; # eV
	$ND      = 7.3e15;
	$E0      = 0.0;
	$Ecenter = $Ec + 0.073;
}

# For #4
if($iSample == 4) {
	$ED      = $Ec - 0.080; # eV
	$ND      = 7.3e15;
	$E0      = 0.018;
	$Ecenter = $Ec + 0.073;
}

# For #3
if($iSample == 3) {
	$ED      = $Ec + 1.000; # eV
	$ND      = 1.3e17;
	$E0      = 0.018;
	$Ecenter = $Ec + 0.05;
}

# For #2
if($iSample == 2) {
	$ED      = $Ec + 0.013; # eV
	$ND      = 4.5e18;
	$E0      = 0.022;
	$Ecenter = $Ec + 0.036;
}

# For #6
if($iSample == 6) {
	$ED      = $Ec + 1.000; # eV
	$ND      = 1.9e19;
	$E0      = 0.008;
	$Ecenter = $Ec + 0.078;
}

# For #1
if($iSample == 1) {
	$ED      = $Ec + 1.000; # eV
	$ND      = 8.3e19;
	$E0      = 0.008;
	$Ecenter = $Ec + 0.078;
}

# For #5
if($iSample == 5) {
	$ED      = $Ec + 1.000; # eV
	$ND      = 9.3e19;
	$E0      = 0.008;
	$Ecenter = $Ec + 0.078;
}

print "\nND=$ND cm-3  ED=$ED eV  E0=$E0 Ecntr=$Ecenter tau0=$tau0\n\n";

my $Eg   = $Ec - $Ev;
my $Nc   = $material->Nc($Effme*$me,  $CBMultiplicity, $T, 0) * 1.0e-6;
my $Dc0  = $material->Dc0($Effme*$me, $CBMultiplicity, 0) * $e**1.5 * 1.0e-6; # cm^-3/eV^1.5
my $Nv   = $material->Nv($Effmh*$me,  $VBMultiplicity, $T, 0) * 1.0e-6;
my $Dv0  = $material->Dv0($Effmh*$me, $VBMultiplicity, 0) * $e**1.5 * 1.0e-6; # cm^-3/eV^1.5
my $Ni   = $material->Ni($T, $Eg);

printf "Nc=%10.4g cm-3  Dc0=%10.4g cm-3/eV1.5\n", $Nc, $Dc0;
printf "Nv=%10.4g cm-3  Dv0=%10.4g cm-3/eV1.5\n", $Nv, $Dv0;
printf "Ni=%10.4g cm-3\n", $Ni;

my $EF     = 1.5;

my $dE     = $T * $KToeV * 0.02;
my $Erange = $T * $KToeV * 20.0 + $Ec - $Ecmin;
my $nE     = int($Erange / $dE + 1.001);
$nE = int($nE / 2) * 2 + 1;
printf "Integration range: %6.3f, %10.3g eV step [%d]\n", $Erange, $dE, $nE;

my (@Energyc, @Dc);
&CalDOSc();
my (@Energyv, @Dv);
&CalDOSv();

my ($Ne, $Nh);
if(0) {
	$Ne = &CalNe($EF, $T);
	$Nh = &CalNh($EF, $T);
	printf "EF=%6.3f: Ne=%10.4g Nh=%10.4g cm-3\n", $EF, $Ne, $Nh;
}
if(0) {
	print "\n";
	my ($EF0, $S0) = &CalEF($T, $EFEPS, $EFSEPS, 0);
	$Ne  = &CalNe($EF0, $T);
	$Nh  = &CalNh($EF0, $T);
	printf  "T=%5.2f:  EF=%6.3f eV: Ne=%10.3g Nh=%10.3g cm-3 S=%8.3g\n",
					$T, $EF0, $Ne, $Nh, $S0;
}

my (@T, @Ne, @Nh, @ue, @sigma, @S);
for(my $i = 0 ; ; $i++) {
	my $Tcal = $Tstart + $i * $Tstep;
	last if($Tcal > $Tend);

	my ($EFopt, $S) = &CalEF($Tcal, $EFEPS, $EFSEPS, 0);
	$Ne  = &CalNe($EFopt, $Tcal);
	$Nh  = &CalNh($EFopt, $Tcal);
	my $NDp = &CalNDp($EFopt, $Tcal);
	my $NAm = &CalNAm($EFopt, $Tcal);
	my $sigma = &CalSigma($EFopt, $Tcal);
# EDが伝導帯に入っている場合、電子濃度 $Ne に非イオン化ドナー密度 $ND-$NDp を加える
	$Ne += $ND - $NDp if($ED >= $Em);
printf  "%5.2f: EF=%5.3feV: Ne=%10.3g NDp=%10.3g s=%10.3g S/cm\n",
				$Tcal, $EFopt, $Ne, $NDp, $sigma;

	$T[$i]     = $Tcal;
	$Ne[$i]    = $Ne;
	$ue[$i]    = $sigma / $e / $Ne;
	$Nh[$i]    = $Nh;
	$sigma[$i] = $sigma;
	$S[$i]     = $S;

	$EF = $EFopt;
}

# 温度依存性データをCSVファイルに保存
my $f = "Transport.csv";
my $out = new JFile($f, "w");
$out->print("T(K),1000/T(K-1),Ne(cm-3),ue(cm2/Vs),Nh(cm-3),sigma(S/cm),S(cm-3)\n");
for(my $i = 0 ; $i < @T ; $i++) {
	my $InvT = 1000.0 / $T[$i];
	$out->print("$T[$i],$InvT,$Ne[$i],$ue[$i],$Nh[$i],$sigma[$i],$S[$i]\n");
}
$out->Close();

exit;

# Percolationを考慮して電気伝導度を計算
sub CalSigma
{
	my ($EF, $T) = @_;

	my $K = -2.0 / 3.0 * $e * $e / $Effme / $me;
	my @f;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Energyc[$i];
		my ($dE, $t);
		if($E >= $Em) {
			$dE = $E-$Em;
			$t  = &tau($dE, $tau0, $r);
		}
		else {
			$dE = 0.0;
			$t  = 0.0;
		}
		my $P = &BarrierProbability($E, $Ecenter, $E0);
		$f[$i] = $K * $t * $dE * &dFEedE($E, $EF, $T) * $Dc[$i] * $P;
	}
	return  Algorism::SinglePointIntegrateBySimpson(\@Energyc, \@f);
}

# エネルギー依存緩和時間を計算
sub tau
{
	my ($E, $tau0, $r) = @_;
	return $tau0 * $E**($r-0.5);
}

# エネルギー$Eを持つ電子がGaussian分布をする障壁を越えられる確率を計算
# 誤差関数erfcを利用して計算している
sub BarrierProbability
{
	my ($E, $Ec, $E0) = @_;

	if($E0 == 0.0) {
		return 1.0 if($E >= $Ec);
		return 0.0;
	}

	my $x = ($E - $Ec) / $E0;
	if($x < 0.0) {
		my $r = erfc(-$x);
		return 0.5 * (1.0 - $r);
	}
	my $ret = erfc($x);
	return 0.5 * (1.0 + $ret);
}

# 裾状態を含んだパラボリック伝導帯の状態密度
sub CalDOSc
{
	my (@Dtail);

	my $IsTail = 1;
	$IsTail = 0 if($Nt0 == 0.0);
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Ecmin + $i * $dE;
		$Energyc[$i] = $E;
		$Dc[$i] = ($E > $Ec)? $Dc0 * sqrt($E - $Ec) : 0.0;
		$Dtail[$i] = $Nt0 * exp(($E-$Eth) / $Eu);
		if($IsTail and $Dtail[$i] > $Dc[$i]) {
			$Dc[$i] = $Dtail[$i];
		}
		elsif($IsTail) {
			my $Ntrap = $Nt0 * $Eu * exp(($E-$Eth) / $Eu);
print "Eborder=$E  Eth=$Eth eV  Ntrap=$Ntrap Nt0=$Nt0 cm-3\n";
			$IsTail = 0;
		}
	}

	my $f = "Dc.csv";
	my $out = new JFile($f, "w");
	$out->print("E,Dc,Dtail\n");
	for(my $i = 0 ; $i < $nE ; $i++) {
		$out->print("$Energyc[$i],$Dc[$i],$Dtail[$i]\n");
	}
	$out->Close();
}

# パラボリック価電子帯の状態密度
sub CalDOSv
{
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Ev - $i * $dE;
		$Energyv[$i] = $E;
		$Dv[$i] = ($E < $Ev)? $Dv0 * sqrt($Ev - $E) : 0.0;
	}
}

# 二分法でフェルミ準位を計算
sub CalEF
{
	my ($T, $xEPS, $SEPS, $IsPrint) = @_;

	my ($pVars, $S) = $opt->BinarySearch(
							$Ev-1.0, $Ec+1.0, sub { &CalS($T, @_); },
							$nMaxIter, $xEPS, $SEPS, $IsPrint
						);
	return ($pVars->[0], $S);
}

# 二分法でフェルミ準位を計算 (二分法ルーチン内蔵)
sub CalEF0
{
	my ($T, $xEPS, $SEPS, $IsPrint) = @_;
#print "\nSearch EF...\n";
#print "EPS: $EPS\n";

	my $Emin = $Ev-1.0;
	my $Emax = $Ec+1.0;
	my $Smin = &CalS($T, $Emin, $IsPrint);
	my $Smax = &CalS($T, $Emax, $IsPrint);
	if($Smin * $Smax > 0.0) {
		print "S($Smin)@$Emin and S($Smax)@$Emax must be negative and positive.\n";
		return (0.0, 0.0);
	}

	for(my $i = 0 ; $i < $nMaxIter ; $i++) {
		return (($Emin+$Emax)/2.0, ($Smin+$Smax)/2.0) if(abs($Emax - $Emin) < $xEPS);

		my $Ehalf  = ($Emin + $Emax) / 2.0;
		my $Shalf = &CalS($T, $Ehalf,      $IsPrint);
#print "$i: E: $Emin - $Ehalf - $Emax: $Smin - $Shalf - $Smax\n";
		return ($Ehalf, $Shalf) if(abs($Shalf) < $SEPS);

		if($Smin * $Shalf < 0.0) {
			$Emax = $Ehalf;
			$Smax = &CalS($T, $Emax, $IsPrint);
		}
		else {
			$Emin = $Ehalf;
			$Smin = &CalS($T, $Emin, $IsPrint);
		}
	}
	return (($Emin+$Emax)/2.0, ($Smin+$Smax)/2.0);
}

# 過剰電荷の計算。過剰電荷＝＝0のEFを探すために使う
sub CalS 
{
	my ($T, $EF, $iPrintLevel) = @_;

	my $Ne  = &CalNe($EF, $T);
	my $Nh  = &CalNh($EF, $T);
	my $NDp = &CalNDp($EF, $T);
	my $NAm = &CalNAm($EF, $T);
#print "EF=$EF  Ne=$Ne ND=$ND, $NDp  NA=$NA, $NAm\n";
	my $dN = $Ne - $NDp - ($Nh - $NAm);
	my $S = $dN;
printf("EF=%6.3f: Ne=$Ne Nh=$Nh (%g)\n", $EF, $S) if($iPrintLevel >= 2);
	return $S;
}

# 伝導帯中の占有電子数
# EDが伝導帯中にある場合でも、ドナー準位を占める電子数は数えていない
sub CalNe
{
	my ($EF, $T) = @_;
	my @Dfc;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Energyc[$i];
		$Dfc[$i] = $Dc[$i] * &FEe($E, $EF, $T);
	}
	return  Algorism::SinglePointIntegrateBySimpson(\@Energyc, \@Dfc);
}

# 価電子中の占有正孔数
# EAが価電子帯中にある場合でも、アクセプター準位を占める正孔数は数えていない
sub CalNh
{
	my ($EF, $T) = @_;
	my @Dfv;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Energyv[$i];
		$Dfv[$i] = $Dv[$i] * &FEh($E, $EF, $T);
	}
	return -Algorism::SinglePointIntegrateBySimpson(\@Energyv, \@Dfv);
}

# イオン化ドナー密度
sub CalNDp
{
	my ($EF, $T) = @_;
	return $ND * &FEh($ED, $EF, $T);
}

# イオン化アクセプター密度
sub CalNAm
{
	my ($EF, $T) = @_;
	return $NA * &FEe($EA, $EF, $T);
}

# 電子のFermi-Dirac分布のエネルギー微分
sub dFEedE
{
	my ($E, $EF, $T) = @_;
	return 0.0 if($T == 0.0);

	my $ex  = exp(($E-$EF) / $T / $KToeV);
	my $ex1 = 1.0 + $ex;
	return -1.0 / $T / $KToeV * $ex / ($ex1*$ex1);
}

# 電子のFermi-Dirac分布
sub FEe
{
	my ($E, $EF, $T) = @_;
	if($T == 0.0) {
		return 0.0 if($E > $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($E-$EF) / $T / $KToeV) );
}

# 正孔のFermi-Dirac分布
sub FEh
{
	my ($E, $EF, $T) = @_;
	if($T == 0.0) {
		return 0.0 if($E < $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($EF-$E) / $T / $KToeV) );
}