package ConductingMaterial;
use IniFile;
use Sci::Material;
@ISA = qw(IniFile Material);

#公開したいサブルーチン
@EXPORT = qw();

use strict;
#use warnings;

use CSV;
use Sci qw($pi $kB $e $e0 $me $h $hbar $KToeV $JToeV erfc);

use Sci::Material;
use Sci::Optimize;

#===============================================
# 大域変数
#===============================================
my $AddNDpToNe = 1;

#============================================================
#　変数等取得関数
#============================================================
sub Thickness        { return shift->{Thicnkess}; }
sub SetThickness     { my($this,$d)=@_; return $this->{Thicnkess} = $d; }
sub SetName          { my($this,$n)=@_; return $this->{Name} = $n; }
sub Name             { my($this)=@_; return $this->{Name}; }
sub Temperature      { return shift->{Temperature}; }
sub SetTemperature   { my($this,$T)=@_; return $this->{Temperature} = $T; }
sub SetEnergyStructure {
	my ($this, $EC, $EV, $EF, $me, $mh, $ND, $ED, $NA, $EA) = @_;
}
sub SetProperties {
	my ($this, $Eps, $mue, $taue, $mup, $taup) = @_;
}
#sub SetVariable
#{
#	my ($this, $Name, $value) = @_;
#	my ($ModelName, $VarName) = ($Name =~ /^(.*?)::(.*)$/);
#	if(defined $ModelName) {
#		my $n = $this->nDielectricModel();
#		my $p = $this->pDielectricModel();
#		if($ModelName =~ /^\d+$/) {
#			return $p->[$ModelName]->{$VarName} = $value;
#		}
#		for(my $i = 0 ; $i < $n ; $i++) {
#			if($p->[$i]->{name} eq $ModelName) {
#				return $p->[$i]->{$VarName} = $value;
#			}
#		}
#		return undef;
#	}
##print "Name=$Name\n";
#	return $this->{$Name} = $value;
#}

#sub pVariable
#{
#	my ($this, $Name) = @_;
#	my ($ModelName, $VarName) = ($Name =~ /^(.*?)::(.*)$/);
#	if(defined $ModelName) {
#		my $n = $this->nDielectricModel();
#		my $p = $this->pDielectricModel();
#		if($ModelName =~ /^\d+$/) {
#			return \$p->[$ModelName]->{$VarName};
#		}
#		for(my $i = 0 ; $i < $n ; $i++) {
#			if($p->[$i]->{name} eq $ModelName) {
##print "Hit: $ModelName=$VarName => ", $p->[$i]->{$VarName}, "\n";
#				return \$p->[$i]->{$VarName};
#			}
#		}
#		return undef;
#	}
##print "Name=$Name\n";
#	$this->{$Name} = '' if(!defined $this->{$Name});
#	return \$this->{$Name};
#}

#============================================================
#　コンストラクタ、デストラクタ
#============================================================
sub new
{
	my ($module, $material, $name, $app, %args) = @_;

	my $self = IniFile->new('', 0, 0);
	my $this = bless $self, ref($module) || $module;
	$this->SetApplication($app) if($app);
	Utils::MergeHash($this, \%args);
	$this->Initialize();
	$this->MakeMaterial($material, $name) if($material and $name);
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　一般関数
#============================================================
sub Initialize
{
	my ($this) = @_;
#	$this->{nDielectricModel}   = 0;
#	$this->{pDielectricModel}   = [];
#	$this->{idxDielectricModel} = {};
}

sub MakeMaterial
{
	my ($this, $material, $name) = @_;
	$this->SetName($material);
	if($material =~ /^specify$/i) {
		$this->SetName($name);
	}
	elsif($material =~ /^air$/i) {
	}
	elsif($material =~ /^c-Si$/i) {
	}
	elsif($material =~ /^a-Si:H$/i) {
	}
	elsif($material =~ /^ZnO$/i) {
	}
	elsif($material =~ /^a-InGaZnO4$/i) {
	}
}

sub PrepareForCalculation
{
	my ($this) = @_;

}

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

	my $Effme    = $this->{Effme};
	my $Em       = $this->{Em};
	my $pEnergyc = $this->{pEnergyc};
	my $pDc      = $this->{pDc};
	my $nE       = @$pEnergyc;
	my $tau0     = $this->{tau0};
	my $r        = $this->{r};
	my $pT       = $this->{PowerT};
	my $tau02    = $this->{tau02};
	my $r2       = $this->{r2};
	my $pT2      = $this->{PowerT2};
	my $Ecenter  = $this->{Ecenter};
	my $E0       = $this->{E0};

	my $K = -2.0 / 3.0 * $e * $e / $Effme / $me;
	my @f;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $pEnergyc->[$i];
		my ($dE, $t);
		if($E >= $Em) {
			$dE = $E-$Em;
			$t  = $this->tau($T, $dE, $tau0, $r, $pT, $tau02, $r2, $pT2);
		}
		else {
			$dE = 0.0;
			$t  = 0.0;
		}
		my $P = $this->BarrierProbability($E, $Ecenter, $E0);
#print "i=$i: $E,$Em: $K, $P, $t, $dE\n";
		$f[$i] = $K * $t * $dE * $this->dFEedE($E, $EF, $T) * $pDc->[$i] * $P;
	}
	my $sigma = Algorism::SinglePointIntegrateBySimpson($pEnergyc, \@f);

# EDが伝導帯に入っている場合、電子濃度 $Ne に非イオン化ドナー密度 $ND-$NDp を加える
if($AddNDpToNe) {
	if($this->{ED} >= $this->{Em}) {
		my $E   = $this->{ED};
		my $NDp = $this->CalNDp($EF, $T);
		my ($dE, $t);
		if($E >= $Em) {
			$dE = $E-$Em;
			$t  = $this->tau($T, $dE, $tau0, $r, $pT, $tau02, $r2, $pT2);
		}
		else {
			$dE = 0.0;
			$t  = 0.0;
		}
		my $P = $this->BarrierProbability($E, $Ecenter, $E0);
		$sigma += $K * $t * ($this->{ND} - $NDp) * $dE * $this->dFEedE($E, $EF, $T) * $P;
	}
}

	return $sigma * 1.0e4;
}

sub CalTauAveraged
{
	my ($this, $EF, $T, $n, $ConsiderBarrierProbability) = @_;
	$ConsiderBarrierProbability = 1 if(!defined $ConsiderBarrierProbability);

	my $Effme    = $this->{Effme};
	my $Em       = $this->{Em};
	my $pEnergyc = $this->{pEnergyc};
	my $pDc      = $this->{pDc};
	my $nE       = @$pEnergyc;
	my $tau0     = $this->{tau0};
	my $r        = $this->{r};
	my $pT       = $this->{PowerT};
	my $tau02    = $this->{tau02};
	my $r2       = $this->{r2};
	my $pT2      = $this->{PowerT2};
	my $Ecenter  = $this->{Ecenter};
	my $E0       = $this->{E0};

#print "t: $tau0 / $r : $tau02 / $r2\n";
	my @f;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $pEnergyc->[$i];
		my ($dE, $t);
		if($E >= $Em) {
			$dE = $E-$Em;
			$t  = $this->tau($T, $dE, $tau0, $r, $pT, $tau02, $r2, $pT2);
		}
		else {
			$dE = 0.0;
			$t  = 0.0;
		}
		$t = $t**$n if($t > 0.0);
		my $P;
		if($ConsiderBarrierProbability) {
			$P = $this->BarrierProbability($E, $Ecenter, $E0);
		}
		else {
			$P = 1.0;
		}
		$f[$i] = $t * $dE * $this->dFEedE($E, $EF, $T) * $pDc->[$i] * $P;
#my $d = $this->dFEedE($E, $EF, $T);
#print "f($i,$E,$EF)=$f[$i] (t=$t, P=$P, d=$d, Dc=$pDc->[$i]))\n";
	}
	my $tauaveraged = -2.0/3.0 * Algorism::SinglePointIntegrateBySimpson($pEnergyc, \@f);

	my $Ne = $this->CalNe($EF, $T);
# EDが伝導帯に入っている場合、電子濃度 $Ne に非イオン化ドナー密度 $ND-$NDp を加える
if($AddNDpToNe) {
	if($this->{ED} >= $this->{Em}) {
		my $E   = $this->{ED};
		my $NDp = $this->CalNDp($EF, $T);
		my ($dE, $t);
		if($E >= $Em) {
			$dE = $E-$Em;
			$t  = $this->tau($T, $dE, $tau0, $r, $pT, $tau02, $r2, $pT2);
		}
		else {
			$dE = 0.0;
			$t  = 0.0;
		}
		$t = $t**$n if($t > 0.0);
		my $P;
		if($ConsiderBarrierProbability) {
			$P = $this->BarrierProbability($E, $Ecenter, $E0);
		}
		else {
			$P = 1.0;
		}
#print "Ne: $Ne ($tauaveraged) => ";
		$tauaveraged += -2.0/3.0 * $t * ($this->{ND} - $NDp) * $this->dFEedE($E, $EF, $T) * $P;
		$Ne += $this->{ND} - $NDp;
#print "$Ne ($tauaveraged)\n";

	}
}

	return $tauaveraged / $Ne;
}

sub CalDriftMobility
{
	my ($this, $EF, $T) = @_;

	my $Effme       = $this->{Effme};
	my $tauaveraged = $this->CalTauAveraged($EF, $T, 1);
	return $e * $tauaveraged / $Effme / $me * 1.0e4;
}

sub CalHallMobility
{
	my ($this, $EF, $T) = @_;

	my $Effme        = $this->{Effme};
	my $tauaveraged  = $this->CalTauAveraged($EF, $T, 1);
	my $tau2averaged = $this->CalTauAveraged($EF, $T, 2);
	my $R = 0.0;
	$R = $tau2averaged / $tauaveraged if($tauaveraged > 0.0);

	return $e * $R / $Effme / $me * 1.0e4;
}

sub CalSeebeck
{
	my ($this, $EF, $T, $r, $iChargePolarity) = @_;
	$r = $this->{r} if(!defined $r);
	$iChargePolarity = -1 if(!defined $iChargePolarity);

	my $Effme    = $this->{Effme};
	my $Em       = $this->{Em};
	my $pEnergyc = $this->{pEnergyc};
	my $pDc      = $this->{pDc};
	my $nE       = @$pEnergyc;
	my $tau0     = $this->{tau0};
	my $pT       = $this->{PowerT};
	my $tau02    = $this->{tau02};
	my $r2       = $this->{r2};
	my $pT2      = $this->{PowerT2};
	my $Ecenter  = $this->{Ecenter};
	my $E0       = $this->{E0};

	my @f;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $pEnergyc->[$i];
		my ($dE, $t);
		if($E >= $Em) {
			$dE = $E-$Em;
			$t  = $this->tau($T, $dE, $tau0, $r, $pT, $tau02, $r2, $pT2);
		}
		else {
			$dE = 0.0;
			$t  = 0.0;
		}
		my $P = $this->BarrierProbability($E, $Ecenter, $E0);
#print "i=$i: $E,$Em: $K, $P, $t, $dE\n";
		$f[$i] = $t * $dE * $dE * $this->dFEedE($E, $EF, $T) * $pDc->[$i] * $P;
	}
	my $S  = -2.0/3.0 * Algorism::SinglePointIntegrateBySimpson($pEnergyc, \@f);
	my $Ne = $this->CalNe($EF, $T);

# EDが伝導帯に入っている場合、電子濃度 $Ne に非イオン化ドナー密度 $ND-$NDp を加える
if($AddNDpToNe) {
	if($this->{ED} >= $this->{Em}) {
		my $E   = $this->{ED};
		my $NDp = $this->CalNDp($EF, $T);
		my ($dE, $t);
		if($E >= $Em) {
			$dE = $E-$Em;
			$t  = $this->tau($T, $dE, $tau0, $r, $pT, $tau02, $r2, $pT2);
		}
		else {
			$dE = 0.0;
			$t  = 0.0;
		}
		my $P = $this->BarrierProbability($E, $Ecenter, $E0);
		$S += -2.0/3.0 * $t * ($this->{ND} - $NDp) * $dE * $dE * $this->dFEedE($E, $EF, $T) * $P;
		$Ne += $this->{ND} - $NDp;
	}
}

	my $tauaveraged  = $this->CalTauAveraged($EF, $T, 1);
	my $R = 0.0;
	$R = $S / $Ne / $tauaveraged if($Ne > 0.0 and $tauaveraged > 0.0);

#	$S = (1.0 / $T) * $R; #  - ($EF - $Em) + $T * $dEFdT
	$S = (1.0 / $T) * ($R - ($EF - $Em)); # + $T * $dEFdT
	return -$S if($iChargePolarity <= 0);
	return $S;
}

sub CalSeebeckWithNonDegenerateApproximation
{
	my ($this, $EF, $T, $r, $iChargePolarity) = @_;
	$r = $this->{r} if(!defined $r);
	$iChargePolarity = -1 if(!defined $iChargePolarity);

	return -($kB/$e) * ($r+2.0 - ($EF-$this->{Ec})*$e/$kB/$T) if($iChargePolarity <= 0);
	return  ($kB/$e) * ($r+2.0 - ($this->{Ev}-$EF)*$e/$kB/$T);
}

sub CalSeebeckWithDegenerateApproximation
{
	my ($this, $EF, $T, $r, $iChargePolarity) = @_;
	$r = $this->{r} if(!defined $r);
	$iChargePolarity = -1 if(!defined $iChargePolarity);

	return -($kB/$e) * ($r+1.0) * ($pi*$pi/3.0) / (($EF-$this->{Ec})*$e/$kB/$T) if($iChargePolarity <= 0);
	return  ($kB/$e) * ($r+1.0) * ($pi*$pi/3.0) / (($this->{Ev}-$EF)*$e/$kB/$T);
}

# エネルギー依存緩和時間を計算
sub tau
{
	my ($this, $T, $E, @taur) = @_;
	return 0.0 if($E < 0.0);
	my $tinv = 0.0;
	my $n = @taur;
	for(my $i = 0 ; $i < $n ; $i += 3) {
		my $tau0 = $taur[$i];
		next if($tau0 == 0.0);
		my $r    = $taur[$i+1];
		return 0.0 if($E == 0.0 and $r != 0.5);
		my $pT   = $taur[$i+2];
		my $kT = ($pT == 0.0)? 1.0 : $T**$pT;

		my $tau  = $tau0 * $kT * $E**($r-0.5);
		return 0.0 if($tau == 0.0);

		$tinv += 1.0 / $tau;
#my $t = 1.0 / $tinv;
#print "i=$i: $tau / $t\n";
	}
	return 1.0 / $tinv;
}

# エネルギー$Eを持つ電子がGaussian分布をする障壁を越えられる確率を計算
# 誤差関数erfcを利用して計算している
sub BarrierProbability
{
	my ($this, $E, $Ec, $E0) = @_;
#print "E=$E  Ec=$Ec\n";
	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 ($this, $Dc0, $Ec, $Nt0, $Eth, $Eu, $nE, $Estart, $Estep) = @_;
#print "Dc0= $Dc0, $Ec, $Nt0, Eth=$Eth, $Eu, nE=$nE, $Estart, $Estep\n";
	my (@Energyc, @Dc);
	my (@Dtail);
	my $IsTail = 1;
	$IsTail = 0 if($Nt0 == 0.0);
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Estart + $i * $Estep;
		$Energyc[$i] = $E;
		$Dc[$i] = ($E > $Ec)? $Dc0 * sqrt($E - $Ec) : 0.0;
		$Dtail[$i] = ($Eu > 0.0)? $Nt0 * exp(($E-$Eth) / $Eu) : 0.0;
		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;
		}
	}
#print "nE=$nE\n";
	$this->{pEnergyc} = \@Energyc;
	$this->{pDc}      = \@Dc;
	return (\@Energyc, \@Dc);
}

# パラボリック価電子帯の状態密度
sub CalDOSv
{
	my ($this, $Dv0, $Ev, $Nt0, $Eth, $Eu, $nE, $Estart, $Estep) = @_;

	my (@Energyv, @Dv);
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Estart + $i * $Estep;
		$Energyv[$i] = $E;
		$Dv[$i] = ($E < $Ev)? $Dv0 * sqrt($Ev - $E) : 0.0;
	}

	$this->{pEnergyv} = \@Energyv;
	$this->{pDv}      = \@Dv;
	return (\@Energyv, \@Dv);
}

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

	my $Emax = 30.0 * $T * $KToeV;
	if($this->{NA} < 1.0e12 and $this->{ND} < 1.0e12 and $this->{Ec}-$this->{Ev} > $Emax) {
		my $EF = $this->CalEFWithApproximation($T);
		return ($EF, -1.0);
	}
	my $opt = new Optimize;
#	$xEPS = $Emax;
	my ($pVars, $S) = $opt->BinarySearch(
					$this->{Ev}-1.0, $this->{Ec}+5.0, sub { $this->CalS($T, @_); },
					$nMaxIter, $xEPS, $SEPS, $IsPrint
				);
	my $EF = $pVars->[0];
if($AddNDpToNe) {
	$opt->SetDumpingFactor(0.0, 0.1, 5);
	($pVars, $S) = $opt->Optimize("ModifiedNewton", [$EF], [1], [0.001], 
				$SEPS, $nMaxIter, $IsPrint,
				sub { $this->CalNewtonS2($T, @_); },
				sub {},
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
}

	$EF = $pVars->[0];
	return ($EF, $S);
}

sub EiFromEffectiveMass
{
	my ($this, $T) = @_;
	my $EF = 0.5*($this->{Ec}+$this->{Ev}) + 0.75*$kB*$T/$e * log($this->{Effmh}/$this->{Effme});
	return $EF;
}

sub CalEFWithApproximation
{
	my ($this, $T) = @_;

	if($this->{ND} > 0.0 and $this->{ND} > $this->{NA}) {
		return $this->{Ec} + $kB*$T * log($this->{ND}/$this->{Nc}) * $JToeV;
	}
	if($this->{NA} > 0.0) {
		return $this->{Ev} - $kB*$T * log($this->{NA}/$this->{Nv}) * $JToeV;
	}
	return $this->EiFromEffectiveMass($T);
}


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

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

sub CalNewtonS2
{
	my ($this, $T, $pVars, $IsPrint) = @_;

	my $EF = $pVars->[0];
	my $Ne  = $this->CalNe($EF, $T);
	my $Nh  = $this->CalNh($EF, $T);
	my $NDp = $this->CalNDp($EF, $T);
	my $NAm = $this->CalNAm($EF, $T);
#print "EF=$EF  Ne=$Ne ND=$ND, $NDp  NA=$NA, $NAm\n";
	my $dN = $Ne - $NDp - ($Nh - $NAm);
	my $S2 = $dN*$dN;
printf("EF=%6.3f: Ne=%10.3g NDp=%10.3g Nh=%10.3g (%g)\n", 
				$EF, $Ne, $NDp, $Nh, $S2) if($IsPrint);
#printf("EF=%6.3f: Ne=%10.3g NDp=%10.3g Nh=%10.3g NAm=%10.3g (%g)\n", 
#				$EF, $Ne, $NDp, $Nh, $NAm, $S) if($IsPrint);
	return $S2;
}

# 伝導帯中の占有電子数
# EDが伝導帯中にある場合でも、ドナー準位を占める電子数は数えていない
sub CalNe
{
	my ($this, $EF, $T) = @_;

	my $pEnergyc = $this->{pEnergyc};
	my $pDc      = $this->{pDc};
	my $nE       = @$pEnergyc;
	my @Dfc;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $pEnergyc->[$i];
		$Dfc[$i] = $pDc->[$i] * $this->FEe($E, $EF, $T);
	}
	return Algorism::SinglePointIntegrateBySimpson($pEnergyc, \@Dfc);
}

# 価電子中の占有正孔数
# EAが価電子帯中にある場合でも、アクセプター準位を占める正孔数は数えていない
sub CalNh
{
	my ($this, $EF, $T) = @_;

	my $pEnergyv = $this->{pEnergyv};
	my $pDv      = $this->{pDv};
	my $nE       = @$pEnergyv;
	my @Dfv;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $pEnergyv->[$i];
		$Dfv[$i] = $pDv->[$i] * $this->FEh($E, $EF, $T);
	}
#	return Algorism::SinglePointIntegrateBySimpson($pEnergyv, \@Dfv);
	return -Algorism::SinglePointIntegrateBySimpson($pEnergyv, \@Dfv);
}

# イオン化ドナー密度
sub CalNDp
{
	my ($this, $EF, $T) = @_;
	my $FE = $this->FEh($this->{ED}, $EF, $T);
#print "$this->{ND}, $this->{ED}, $EF, $T, $FE\n";
	return $this->{ND} * $FE;
}

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

# 電子のFermi-Dirac分布のエネルギー微分
sub dFEedE
{
	my ($this, $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 ($this, $E, $EF, $T) = @_;
#print "e: $this,$E,$EF,$T\n";
	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 ($this, $E, $EF, $T) = @_;
#print "h: $this,$E,$EF,$T\n";
	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) );
}

1;