#===============================================
# Material
#===============================================
package Material;
use Common;
@ISA = qw(Common);

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

use strict;
use Math::Complex;

use Utils;
use JFile;
use CSV;
use GraphData;
use Sci::GeneralFileFormat;
use Sci qw($pi $e0 $e $c $h $hbar $me $kB $JToeV);

my $pi    = Sci::pi();
my $e0    = Sci::e0();
my $e     = Sci::e();
my $c     = Sci::c();
my $h     = Sci::h();
my $hbar  = Sci::hbar();
my $me    = Sci::me();
my $kB    = Sci::kB();
my $JToeV = Sci::JToeV();

my @y = (0, 0.05,   0.1,    0.15,   0.2,    0.25,   0.3,    0.35,   0.4,
	     0.45,   0.5,    0.55,   0.6,    0.65,   0.7,    0.75,   0.8,    0.85,   0.9,    0.95,   1);
my @v = (1, 0.9984, 0.9817, 0.9622, 0.937,  0.9068, 0.8718, 0.8323, 0.7888, 
		 0.7413, 0.69,   0.6351, 0.5768, 0.5152, 0.4504, 0.3825, 0.3117, 0.2379, 0.1613, 0.082,  0);
my @t = (1, 1.0011, 1.0036, 1.007,  1.0111, 1.0157, 1.0207, 1.0262, 1.0319, 
		 1.0378, 1.0439, 1.0502, 1.0565, 1.0631, 1.0697, 1.0765, 1.0832, 1.09,   1.0969, 1.1037, 1.1107);
my $vArray = new GraphData;
$vArray->SetXDataArray(0, \@y);
$vArray->SetYDataArray(0, \@v);
$vArray->CalMinMax();
my $tArray = new GraphData;
$tArray->SetXDataArray(0, \@y);
$tArray->SetYDataArray(0, \@t);
$tArray->CalMinMax();


sub new
{
	my ($module, %arg) = @_;
	my $this = {};
	Utils::MergeHash($this, \%arg);
	bless $this;

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub PNJunctionParameters
{
	my ($this, $Vbi,
			   $ND, $nDielectricConstant,
			   $NA, $pDielectricConstant,
			   $IsPrint) = @_;
	my $Ntot = $ND + $NA;
	my $en = $nDielectricConstant;
	my $ep = $pDielectricConstant;
	my $W = sqrt( 2.0/$e/$NA/$ND * $Ntot*$Ntot / ($ep*$NA + $en*$ND) * $en*$ep * $Vbi );
	my $xn = $W / $Ntot * $NA;
	my $xp = $W / $Ntot * $ND;
	my $Enmax = $e * $ND / $nDielectricConstant * $xn;
	my $Epmax = $e * $NA / $pDielectricConstant * $xp;
	if($IsPrint) {
		$this->printf("Depletion layer thickness: %12.6g nm\n", $W * 1e9);
		$this->printf("  n: %12.6g nm\n", $xn * 1e9);
		$this->printf("  p: %12.6g nm\n", $xp * 1e9);
		$this->printf("Maxium electric field:\n");
		$this->printf("  n: %12.6g V/cm\n", $Enmax * 1.0e-2);
		$this->printf("  p: %12.6g V/cm\n", $Epmax * 1.0e-2);
	}
}

sub EiFromDOS
{
	my ($this, $Temperature, $IsPrint) = @_;
	my $EF = 0.5*($this->{Ec}+$this->{Ev}) + 0.5*$kB*$Temperature/$e * log($this->{Nv}/$this->{Nc});
	$this->print("EF=%g eV\n", $this->{EF}) if($IsPrint);
	return $EF;
}

sub EF
{
	my ($this, $Temperature, $IsPrint) = @_;
	my $EF = 0.5*($this->{Ec}+$this->{Ev}) + 0.5*$kB*$Temperature * log($this->{Nv}/$this->{Nv});
	$this->print("EF=%g eV\n", $this->{EF}) if($IsPrint);
	return $EF;
}

sub Eg
{
	my ($this, $Ec, $Ev, $IsPrint) = @_;
	$this->{Ec} = $Ec;
	$this->{Ev} = $Ev;
	$this->{Eg} = $Ec - $Ev;
	$this->print("Eg=%g eV\n", $this->{Eg}) if($IsPrint);
	return $this->{Eg};
}

sub nh
{
	my ($this, $Temperature) = @_;
	my $ne = $this->{Nc} * exp( -($this->{Ec}-$this->{EF}) / $kB / $Temperature );
	return $ne;
}

sub ne
{
	my ($this, $Temperature) = @_;
	my $ne = $this->{Nc} * exp( -($this->{Ec}-$this->{EF}) / $kB / $Temperature );
	return $ne;
}

sub Ni
{
	my ($this, $Temperature, $Eg, $IsPrint) = @_;
	$this->{Ni} = sqrt($this->{Nc} * $this->{Nv} * exp(-$Eg * $e / $kB / $Temperature));
	$this->printf("  Ni = %12.6g cm-3\n", $this->{Ni}*1.0e-6) if($IsPrint);
	return $this->{Ni};
}

sub Nv
{
	my ($this, $HoleEffectiveMass, $VBMultiplicity, $Temperature, $IsPrint) = @_;
	$this->{Nv} = 2.0 * $VBMultiplicity * (2.0*$pi * $HoleEffectiveMass * $kB * $Temperature / $h / $h)**1.5;
	$this->printf("  Nv = %12.6g cm-3\n", $this->{Nv}*1.0e-6) if($IsPrint);
	return $this->{Nv};
}

sub Nc
{
	my ($this, $ElectronEffectiveMass, $CBMultiplicity, $Temperature, $IsPrint) = @_;
	$this->{Nc} = 2.0 * $CBMultiplicity * (2.0*$pi * $ElectronEffectiveMass * $kB * $Temperature / $h / $h)**1.5;
	$this->printf("  Nc = %12.6g cm-3\n", $this->{Nc}*1.0e-6) if($IsPrint);
	return $this->{Nc};
}

sub Dv0
{
	my ($this, $HoleEffectiveMass, $VBMultiplicity, $IsPrint) = @_;
	$this->{Dv0} = $VBMultiplicity / 2.0 / $pi / $pi * (2.0 * $HoleEffectiveMass / $hbar / $hbar)**(1.5);
	$this->printf("  Dv(Ev-E) = %12.6g*sqrt(Ev-E)\n", $this->{Dv0}) if($IsPrint);
	return $this->{Dv0};
}

sub Dc0
{
	my ($this, $ElectronEffectiveMass, $CBMultiplicity, $IsPrint) = @_;
	$this->{Dc0} = $CBMultiplicity / 2.0 / $pi / $pi * (2.0 * $ElectronEffectiveMass / $hbar / $hbar)**(1.5);
	$this->printf("  Dc(Ec-E) = %12.6g*sqrt(Ec-E)\n", $this->{Dc0}) if($IsPrint);
	return $this->{Dc0};
}

sub CalEFFromNe
{
	my ($this, $Ne, $Temperature, $IsPrint) = @_;
	my $EF = $this->{Ec} + $kB * $Temperature / $e * log($Ne / $this->{Nc});
	$this->printf("  EF at Ne=%g cm-3 = %12.6g eV\n", $Ne*1.0e-6, $EF) if($IsPrint);
	return $EF;
}

sub CalEFFromNh
{
	my ($this, $Nh, $Temperature, $IsPrint) = @_;
	my $EF = $this->{Ev} - $kB * $Temperature / $e * log($Nh / $this->{Nv});
	$this->printf("  EF at Nh=%g cm-3 = %12.6g eV\n", $Nh*1.0e-6, $EF) if($IsPrint);
	return $EF;
}

sub DOS
{
	my ($this, $E) = @_;
	if($E >= $this->{Ec}) {
		return $this->{Dc0} * sqrt($E - $this->{Ec});
	}
	if($E <= $this->{Ev}) {
		return $this->{Dv0} * sqrt($this->{Ev} - $E);
	}
	return 0.0;
}

sub BMShift
{
	my ($this, $Ne, $D0, $IsPrint) = @_;
#print "	 ($Ne / $D0 * 1.5)**(2.0/3.0);\n";
	my $EF = ($Ne / $D0 * 1.5)**(2.0/3.0);
	$this->printf("EBM=EF(metal): %6.3f eV\n", $EF / $e) if($IsPrint);
	return $EF;
}

sub ElectronDiffusionLength
{
	my ($this, $ElectronDiffusionConstant, $ElectronLifeTime, $IsPrint) = @_;
	my $ElectronDiffusionLength = sqrt($ElectronDiffusionConstant * $ElectronLifeTime);
	$this->printf("  Diffusion length  : %g nm for lifetime of $ElectronLifeTime s\n", 
				$ElectronDiffusionLength * 1.0e9, $ElectronLifeTime) if($IsPrint);
	return $ElectronDiffusionLength;
}

sub ElectronDiffusionConstant
{
	my ($this, $Temperature, $ElectronMobility, $IsPrint) = @_;
	my $ElectronDiffusionConstant = $ElectronMobility / $e * $kB * $Temperature;
	$this->printf("  Diffusion constant: %g cm2/s\n", $ElectronDiffusionConstant * 1.0e4) if($IsPrint);
	return $ElectronDiffusionConstant;
}

sub ElectronDriftVelocity
{
	my ($this, $V, $Thickness, $ElectronMobility, $IsPrint) = @_;
	my $E = $V / $Thickness;
	my $vdrift = $ElectronMobility * $E;
	$this->printf("  Drift velocity(e): %e m/s\n", $vdrift) if($IsPrint);
	return $vdrift;
}

sub ElectronDriftMeanFreePath
{
	my ($this, $V, $Thickness, $ElectronMobility, $MomentumRelaxationTime, $IsPrint) = @_;
	my $vdrift = $this->ElectronDriftVelocity($V, $Thickness, $ElectronMobility, $IsPrint);
	my $ldrift = $vdrift * $MomentumRelaxationTime;
	$this->printf("  Mean free path(vdrift): %g nm\n", $ldrift * 1.0e9) if($IsPrint);
	return $ldrift;
}

sub ElectronThermalMeanFreePath
{
	my ($this, $Temperature, $ElectronEffectiveMass, $MomentumRelaxationTime, $IsPrint) = @_;
	my $vth = $this->ElectronThermalVelocity($Temperature, $ElectronEffectiveMass, $IsPrint);
	my $lth = $vth * $MomentumRelaxationTime;
	$this->printf("  Mean free path(vth): %g nm\n", $lth * 1.0e9) if($IsPrint);
	return $lth;
}

sub ElectronThermalVelocity
{
	my ($this, $Temperature, $ElectronEffectiveMass, $IsPrint) = @_;
$ElectronEffectiveMass = 1.0e-100 if($ElectronEffectiveMass == 0.0);
	my $vth = sqrt(3.0 * $kB * $Temperature / $ElectronEffectiveMass);
	$this->printf("  Thermal velocity(e): %e m/s\n", $vth) if($IsPrint);
	return $vth;
}

sub ElectronDriftConductivity
{
	my ($this, $ElectronMobility, $Ne, $IsPrint) = @_;
	my $sigma = $e * $Ne * $ElectronMobility;
	$this->printf("  Conductivity(e)    : %g S/cm\n",  $sigma * 1.0e2) if($IsPrint);
}

sub ElectronDriftMobility
{
	my ($this, $ElectronEffectiveMass, $MomentumRelaxationTime, $IsPrint) = @_;
	my $ElectronMobility = $e * $MomentumRelaxationTime / $ElectronEffectiveMass;
	$this->printf("  Mobility(e)    : %g cm2/Vs\n",  $ElectronMobility * 1.0e4) if($IsPrint);
	return $ElectronMobility;
}

sub EmissionCurrent 
{
	my ($this, $V, $Thickness, $Temperature, 
		$EffectiveRichardson, $WorkFunction, 
		$DielectricConstant, $ElectrodeSpacing,
		$IsPrint) = @_;
	my $E = $V / $Thickness;
	my $J0 = $EffectiveRichardson * $Temperature * $Temperature
				* exp( -$WorkFunction / $kB / $Temperature );

	my $SchottkySlope   = sqrt($e*$e*$e / 4.0 / $pi / $DielectricConstant);
	my $SchottkySlopekT = $SchottkySlope / $kB / $Temperature;

# Space charge limited current: Child-Langmuir's equation
	my $KSCLC = 4.0/9.0 * $e0 * sqrt(2.0*$e/$me);

# Fowler-Nordheim's equation
	my $KFN  = $e*$e*$e / 8.0 / $pi / $h;
	my $eKFN = 8.0*$pi/3.0/$h/$e*sqrt(2.0*$me);

	my $JSchottky = ($E >= 0.0) ? $J0 * exp($SchottkySlopekT * sqrt($E)) : 0.0;
	my $JSCLC     = ($V >= 0.0) ? $KSCLC * $V**1.5 / $ElectrodeSpacing / $ElectrodeSpacing : 0.0;
	my $JFN       = ($E > 0.0)  ? $KFN * $E*$E/$WorkFunction * exp(-$eKFN * $WorkFunction**1.5 / $E) : 0.0;
		
	my $y = ($E >= 0) ? 3.79e-4 * sqrt($E*1.0e-2) / ($WorkFunction / $e) : 0.0;
	my $v = $vArray->YVal(0, $y);
	my $t = $tArray->YVal(0, $y);
	my $JFNCorrected = ($E > 0.0)  ? $KFN * $E*$E/$WorkFunction /$t/$t * exp(-$eKFN * $WorkFunction**1.5 / $E * $v) : 0.0;

	if($IsPrint) {
		$this->printf("Thermionic current(No Schottky): J0 = %e A/cm2\n", $J0 * 1.0e-4) if($IsPrint);
		$this->printf("Shottky Slope: %12.6g/kT = %12.6g (m/V)^0.5\n", $SchottkySlope, $SchottkySlopekT) if($IsPrint);
		$this->printf("JSchottky=%12.6g A/cm2\n", $JSchottky*1.0e-4);
		$this->printf("JSCLC    =%12.6g A/cm2\n", $JSCLC*1.0e-4);
		$this->printf("JFN      =%12.6g A/cm2\n", $JFN*1.0e-4);
		$this->printf("JFN(Corr)=%12.6g A/cm2\n", $JFNCorrected*1.0e-4);
	}
}

sub CalEffectiveMassFromRichardsonDushman
{
	my ($this, $EffectiveRichardson, $IsPrint) = @_;
	my $mef = $EffectiveRichardson / (4.0*$pi * $e * $kB * $kB / $h / $h / $h);
	return $mef;
}

sub RichardsonDushman
{
	my ($this, $ElectronEffectiveMass, $IsPrint) = @_;

	my $Richardson = 4.0*$pi * $me * $e * $kB * $kB / $h / $h / $h;
	my $mef = $ElectronEffectiveMass / $me;
	my $EffectiveRichardson = $Richardson * $mef;
	if($IsPrint) {
		$this->printf("Richardson Constant: %e  A/m2/K2\n", $Richardson);
		$this->printf("Effective Ricahrdson constant=%6.3fme: %e  A/m2/K2\n", $mef, $EffectiveRichardson);
	}
	return ($Richardson, $EffectiveRichardson);
}

sub SchottkyFrenkelSlope
{
	my ($this, $Temperature, $DielectricConstant, $IsPrint) = @_;
	my $SchottkySlope   = sqrt($e*$e*$e / 4.0 / $pi / $DielectricConstant);
	my $SchottkySlopekT = $SchottkySlope / $kB / $Temperature;
	printf("Shottky Slope: %12.6g/kT = %12.6g (m/V)^0.5\n", $SchottkySlope, $SchottkySlopekT);
	return ($SchottkySlope, $SchottkySlopekT);
}

sub PooleFrenkelCurrent
{
	my ($this, $V, $Thickness, $Temperature, $PFSlopekT, $U, $NIon, $dHopping, $fIon) = @_;
	my $E = $V / $Thickness;
	my $eFactor = exp(-$U / $kB / $Temperature);
	my $K = 2.0 * $NIon  * $e * $dHopping * $fIon * $eFactor;
	my $eK = $e * $dHopping / 2.0 / $kB / $Temperature;
	my $J = $K * Sci::sinh($eK * $E);
	my $JPF = $J * exp( $PFSlopekT * sqrt(abs($E)) );
	return $JPF;
}

sub PooleFrenkelSlope
{
	my ($this, $Temperature, $DielectricConstant, $IsPrint) = @_;
	my $PFSlope       = sqrt($e*$e*$e / $pi / $DielectricConstant);
	my $PFSlopekT = $PFSlope / $kB / $Temperature;
	printf("Poole-Frenkel Slope: %12.6g/kT = %12.6g (m/V)^0.5\n", $PFSlope, $PFSlopekT) if($IsPrint);
}

sub HoppingMobility
{
	my ($this, $Temperature, $NIon, $U, $dHopping, $fIon, $IsPrint) = @_;
	my $sigmaOhmic = $this->HoppingSigma($Temperature, $NIon, $U, $dHopping, $fIon, $IsPrint);
	my $mobilitySigma = $sigmaOhmic / $NIon / $e;
	my $mobilityEinstein = $e * $dHopping * $dHopping * $fIon / $kB / $Temperature;
	if($IsPrint) {
		$this->printf("Mobility(sigma)         : %g cm2/Vs\n", $mobilitySigma * 1.0e4);
		$this->printf("Mobility(Einstein's law): %g cm2/Vs\n", $mobilityEinstein * 1.0e4);
	}
	return ($mobilitySigma, $mobilityEinstein);
}

sub HoppingSigma
{
	my ($this, $Temperature, $NIon, $U, $dHopping, $fIon, $IsPrint) = @_;
	my $eFactor = exp(-$U / $kB / $Temperature);
	my $K = 2.0 * $NIon  * $e * $dHopping * $fIon * $eFactor;
	my $eK = $e * $dHopping / 2.0 / $kB / $Temperature;
	my $sigmaOhmic = $K * $eK;
	$this->printf("sigma(ohmic) = %g S/cm\n", $sigmaOhmic * 1.0e-2) if($IsPrint);
	return $sigmaOhmic;
}

sub CalPlasmaFrequency
{
# Ne: m-3
	my ($this, $Ne, $m) = @_;
#print "e=$e  me=$me  e0=$e0  pi=$pi\n";
#print "Ne=$Ne  m=$me\n";
	return sqrt( $Ne *$e*$e / $me / $m / $e0) / 2.0 / $pi;
}

sub CalMISCV
{
	my ($this, $T, $Vg, $VFB, $Epss, $Epsi, $dgate, $ds, $Cox, $Cs,
			$Ss, $LDe, $LDh, $LD, $n0, $p0, $CorrectByCs, $EPS, $nMaxIter, $Vs, $dVs) = @_;
	$EPS = 1.0e-6    if(!defined $EPS);
	$nMaxIter = 100  if(!defined $nMaxIter);
	$Vs  = $Vg * 0.5 if(!defined $Vs);
	$dVs = 1.0e-4    if(!defined $dVs);

	my ($Vox, $Qs, $EField, $HFCg, $LFCg, $HFCd, $LFCd);
	for(my $i = 0 ; $i < $nMaxIter ; $i++) {
		my $Vgm1;
		($Vgm1, $Vox, $EField, $Qs, $HFCg, $LFCg, $HFCd, $LFCd) 
				= $this->CalMISCVFromVs($T, $Vs-$dVs, $VFB, $Epss, $Epsi, $dgate, $ds, $Cox, $Cs, 
											$Ss, $LDe, $LDh, $LD, $n0, $p0, $CorrectByCs);
		my ($Vgp1, $Vox2, $EField2, $Qs2, $HFCg2, $LFCg2, $HFCd2, $LFCd2) 
				= $this->CalMISCVFromVs($T, $Vs+$dVs, $VFB, $Epss, $Epsi, $dgate, $ds, $Cox, $Cs, 
											$Ss, $LDe, $LDh, $LD, $n0, $p0, $CorrectByCs);
		my $dVgdVs = ($Vgp1 - $Vgm1) / $dVs;
		my $dVs = -($Vgm1+$Vgp1 - 2*$Vg) / $dVgdVs;
#print "  i=$i: Vs=$Vs + $dVs => ";
		$Vs += $dVs;
#print "$Vs\n";
		last if(abs($dVs) < $EPS);
	}

	return ($Vs, $Vox, $EField, $Qs, $HFCg, $LFCg, $HFCd, $LFCd);
}

sub CalMISCVFromVs
{
	my ($this, $T, $Vs, $VFB, $Epss, $Epsi, $dgate, $ds, $Cox, $Cs,
			$Ss, $LDe, $LDh, $LD, $n0, $p0, $CorrectByCs) = @_;

	my $n0p0 = $n0 / $p0;
	my $bVs = $e * $Vs / $kB / $T;
	my $Fs = sqrt( exp(-$bVs) + $bVs - 1.0 + $n0p0 * (exp($bVs) - $bVs - 1.0) );
	$Fs = -$Fs if($Vs < 0.0);
	my $EField = 2.0 * $kB * $T / $e / $LD * $Fs ;

	my $Qs  = -$Epss * $EField; #-2.0 * $Epss / $LD * $Fs * $kB * $T / $e;
	my $Vox = -$dgate * $Qs / $Epsi;
	my $Vg  = $Vs + $Vox + $VFB;

	my ($LFCd, $HFCd);
	if(abs($bVs) < 1.0e-4) {
		$LFCd = 1.414213562 * sqrt(1.0 + $n0p0) * $Epss / $LD;
		$LFCd *= 1.0 - 2.0/3.0 * (-$p0+$n0) / ($p0+$n0) * $bVs;
	}
	else {
		$LFCd = $Epss / $LD * (-exp(-$bVs) + 1.0 + $n0p0 * (exp($bVs) - 1.0)) / $Fs;
	}
	$LFCd = $LFCd * $Ss;
#	my $Vbi = -$Vs;
	my $Vbi = -$Vs - $kB * $T / $e;
	my $CDe = sqrt(0.5 * $e * $n0 * $Epss / abs($Vbi)) * $Ss;
	my $CDh = sqrt(0.5 * $e * $p0 * $Epss / abs($Vbi)) * $Ss;
	my $CDepletion = ($CDe > $CDh)? $CDe : $CDh;
	if($CorrectByCs) {
		if($CDepletion < $Cs) {
			$CDepletion = $Cs;
		}
		if($LFCd < $CDepletion) {
			$LFCd = $CDepletion;
		}
	}
	if($Vbi > 0.0) {
		$HFCd = $CDepletion;
	}
	else {
		$HFCd = $LFCd;
	}
#print "C=$Cox, $LFCd, $HFCd\n";
	my $LFCg = $Cox * $LFCd / ($Cox + $LFCd);
	my $HFCg = $Cox * $HFCd / ($Cox + $HFCd);

	return ($Vg, $Vox, $EField, $Qs, $HFCg, $LFCg, $HFCd, $LFCd);
}

sub CalMISCVDebyeLength
{
	my ($this, $T, $Ncarrier, $er) = @_; # m
#	$Ncarrier = 1e-100 if($Ncarrier <= 0.0);
	return sqrt(2.0*$er*$e0 * $kB * $T / $e / $e / $Ncarrier);
}

sub CalSubthresholdSwing
{
	my ($this, $T, $er, $dgate, $Dit) = @_; # m
	my $C = $this->CalCapacitance($er, 1.0*1.0, $dgate); # F/m2
#print "C=$C F/m2\n";
#print "eDit/C=", $e * $Dit / $C, "\n";	
#print "e0=$e0\n";
	return 0.0 if($C == 0.0);
	return $kB * $T * log(10.0) * (1.0 + $e * $Dit / $C) / $e;
}

sub CalDit
{
	my ($this, $T, $er, $S, $dgate, $SS) = @_; # m
	my $Cox = &CalC($er, $S, $dgate);
	return ($Cox / $e) * ($SS * $e / ($kB * $T * log(10.0)) - 1.0);
}

sub CalCapacitance
{
	my ($this, $er, $S, $dgate) = @_; # m
	return 0.0 if($dgate == 0.0);
	my $C = $er * $e0 * $S / $dgate; # F/m2
#print "C=$C m  S=$S m2  d=$dgate m\n";
	return $C;
}

sub CalCapacitanceThickness
{
	my ($this, $er, $S, $C) = @_; # m
	return $er * $e0 * $S / $C;
}


1;