#!/usr/bin/perl -w

use lib 'D:/Programs/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 $NeBase      = 1.0e18;
my $NeBasestart = 1.0e10;
my $NeBaseend   = 1.0e21;
my $nNeBase     =  12;
my $kNeBase     =  (log($NeBaseend / $NeBasestart) / log(10) / ($nNeBase-1));
#print "k0=$kNeBase\n";
$kNeBase = 10**$kNeBase;
#print "k0=$kNeBase\n";

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

# ドーパント
my $ND = 0.0e15;
my $ED = $Ec - 0.38; # eV
my $NA = 1.0e18;
my $EA = $Ev + 0.55; # 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.00;      # エネルギー幅 (Gaussian: exp( -[($E-$Ecenter)/$E0]^2 ))
my $Ecenter = $Ec + 0.0; # エネルギー中央値

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

print  "\n";
print  "ND=$ND cm-3  ED=$ED eV  E0=$E0 Ecntr=$Ecenter tau0=$tau0\n\n";
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);
my (@T, @EF, @Ne0, @Ne, @Nh, @ue, @sigma, @S);

#print "k=$kNeBase\n";
my $Ne0 = -$NeBasestart;
for(my $i = 0 ; $i < $nNeBase ; $i++) {
#	print "i=$i: Ne=$Ne0\n";
	my ($EFopt, $S) = &CalEF($T, $Ne0, $EFEPS, $EFSEPS, 0);
	$Ne  = &CalNe($EFopt, $T);
	$Nh  = &CalNh($EFopt, $T);
	printf "Ne=%g cm-3: EF=%f eV: Ne=%7.4g Nh=%7.4g  (S=%5.3g)\n", $Ne0, $EFopt, $Ne, $Nh, $S;
	$T[$i]   = $T;
	$EF[$i]  = $EFopt;
	$Ne0[$i] = $Ne0;
	$Ne[$i]  = $Ne;
	$Nh[$i]  = $Nh;
	$S[$i]   = $S;
	$Ne0 *= $kNeBase;
}
# データをCSVファイルに保存
my $f = "Transport.csv";
my $out = new JFile($f, "w");
$out->print("T(K),1000/T(K-1),Ne0(cm-3),EF(eV),Ne(cm-3),Nh(cm-3),S(cm-3)\n");
for(my $i = 0 ; $i < @T ; $i++) {
	my $InvT = 1000.0 / $T[$i];
	$out->print("$T[$i],$InvT,$Ne0[$i],$EF[$i],$Ne[$i],$Nh[$i],$S[$i]\n");
}
$out->Close();
exit;

for(my $i = 0 ; ; $i++) {
	my $Tcal = $Tstart + $i * $Tstep;
	last if($Tcal > $Tend);

	my ($EFopt, $S) = &CalEF($Tcal, $NeBase, $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, $NeBase, $xEPS, $SEPS, $IsPrint) = @_;

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

# 二分法でフェルミ準位を計算 (二分法ルーチン内蔵)
sub CalEF0
{
	my ($T, $NeBase, $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, $NeBase, $Emin, $IsPrint);
	my $Smax = &CalS($T, $NeBase, $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, $NeBase, $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, $NeBase, $Emax, $IsPrint);
		}
		else {
			$Emin = $Ehalf;
			$Smin = &CalS($T, $NeBase, $Emin, $IsPrint);
		}
	}
	return (($Emin+$Emax)/2.0, ($Smin+$Smax)/2.0);
}

# 過剰電荷の計算。過剰電荷＝＝0のEFを探すために使う
sub CalS 
{
	my ($T, $NeBase, $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) - $NeBase;
	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) );
}