#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';

use strict;

use GraphData;
use Sci::Science;
use Sci::Material;
my $Material = new Material;

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 $OutputFile = "IV.csv";
my $Temperature = 300.0; # K

my $AppliedVoltage = 10.0; # V
my $Thickness      = 1.0 * 1.0e-6; # m Thickness

my $ElectronEffectiveMass = 0.067 * $me;     # kg
my $DielectricConstant    = 1.0 * $e0;
my $MomentumRelaxationTime = 2.0 * 1.0e-15; # s
my $Ne = 1.0e10 * 1.0e6;                    # m-3  Electron density
my $ElectronLifeTime       = 1.0 * 1.0e-6;  # s    Life time of electron

my $VStart = -1000.0; # V
my $VEnd   =  1000.0; # V
my $nV     = 101;
my $VStep  = ($VEnd-$VStart) / ($nV-1);

#For IonHopping1D
my $NIon = 1.0e20 * 1.0e6; # m-3  Ion density
my $fIon = 1.0e12;         # Hz   Vibration frequency
my $dHopping = 3.0e-10;    # m    Hopping distance
my $U = 1.0 * $e;          # J    Activation energy for hopping
#IonHopping1D();

#For CalSemiconductorParameters
#CalSemiconductorParameters();

#For ElectronEmission
my $MetalNe            = 5.0e22 * 1.0e6; # m-3  Electron density
my $WorkFunction       = 2.5 * $e; # J
my $ElectrodeSpacing   = 1.0e-6; # m
#ElectronEmission();

#For PNJunction
my $DonorDensity    = 1.0e17 * 1.0e6; # m-3
my $AcceptorDensity = 1.0e17 * 1.0e6; # m-3
my $nDielectricConstant = 11.9 * $e0;
my $pDielectricConstant = 11.9 * $e0;
my $nFermiLevel = -3.5 * $e; # J
my $pFermiLevel = -4.5 * $e; # J
PNJunction();

#For SchottkyDiode
my $Vbi = 0.5; # eV
my $SchottkyBarrierHeight = 1.0; # eV
my $ElectronMobility = 10.0 * 1.0e4; # m2/Vs
$VStart = -3.0; # V
$VEnd   =  1.0; # V
$nV     = 101;
$VStep  = ($VEnd-$VStart) / ($nV-1);
#SchottkyDiode();

exit;

sub SchottkyDiode()
{
	my $RichardsonDushman = 4.0*$pi * $me * $e * $kB * $kB / $h / $h / $h;
	printf("Richardson Constant: %e  A/m2/K2\n", $RichardsonDushman);
	my $mef = $ElectronEffectiveMass / $me;
	my $RichardsonDushmanef = $RichardsonDushman * $mef;
	printf("Effective Richardson constant=%6.3fme: %e  A/m2/K2\n", 
			$mef, $RichardsonDushmanef);

	printf("DielectricConstant(n): $nDielectricConstant\n");
	printf("ND: %12.6g cm-3\n", $DonorDensity * 1e-6);
	my $Ne = $DonorDensity;
	printf("Ne: %12.6g cm-3\n", $DonorDensity * 1e-6);
	printf("$ElectronMobility: %12.6g cm2/Vs\n", $ElectronMobility * 1.0e-4);
	printf("Vbi: %12.6g V\n", $Vbi);
	printf("$SchottkyBarrierHeight: %12.6g V\n", $SchottkyBarrierHeight);

	my $InvkT  = 1.0 / $kB / $Temperature;
	my $eInvkT = $e * $InvkT;
	my $ed = ($e / $nDielectricConstant)**1.5;
	my $K1 = $e * $InvkT * sqrt($ed * sqrt(2.0)/4.0/$pi * sqrt($DonorDensity));
	my $K2 = exp(-$SchottkyBarrierHeight * $eInvkT);
	my $KD1 = $e * $ElectronMobility * $Ne * sqrt(2.0 * $e * $DonorDensity / $nDielectricConstant);
	my $KD2 = exp(-$eInvkT * $Vbi);
	
#return;
	open(OUT,">$OutputFile") or die "$!\n";

	print "V(V),J(TE)(A/cm2),J(w/ Schottky eff)(A/cm2),J(Diff)(A/cm2)\n";
	print OUT "V(V),J(TE)(A/cm2),J(TE w/ Schottky eff)(A/cm2),J(Diff)(A/cm2),|J(TE)|(A/cm2),|J(TE w/ Schottky eff)|(A/cm2),|J(Diff)|(A/cm2),\n";
	for(my $i = 0 ; $i < $nV ; $i++) {
		my $V = $VStart + $i * $VStep;

		my $JTE1 = $RichardsonDushmanef * $Temperature * $Temperature * $K2 * ( exp($eInvkT * $V) - 1.0 );
		my $JTE2 = ($Vbi > $V) ? $JTE1 * exp($K1 * sqrt($Vbi - $V)) : $JTE1;
		my $JDiff = ($Vbi > $V) ? $KD1 * sqrt($Vbi - $V) * $KD2 * ( exp($eInvkT * $V) - 1.0 ) : 0.0;

		printf("%12.6g,%12.6g,%12.6g,%12.6g\n", 
			$V, $JTE1*1.0e-4, $JTE2*1.0e-4, $JDiff*1.0e-4);
		printf(OUT "%12.6g,%12.6g,%12.6g,%12.6g,%12.6g,%12.6g,%12.6g\n", 
			$V, $JTE1*1.0e-4, $JTE2*1.0e-4, $JDiff*1.0e-4,
			abs($JTE1)*1.0e-4, abs($JTE2)*1.0e-4, abs($JDiff)*1.0e-4);
	}

	close(OUT);
}

sub PNJunction()
{
	printf("DielectricConstant(n): $nDielectricConstant\n");
	printf("DielectricConstant(p): $pDielectricConstant\n");
	printf("ND: %12.6g cm-3\n", $DonorDensity * 1e-6);
	printf("NA: %12.6g cm-3\n", $AcceptorDensity * 1e-6);

	my $Vbi = ($nFermiLevel - $pFermiLevel) / $e; # V
	printf("Built-in potential: %12.6g V\n", $Vbi);

	$Material->Nc($ElectronEffectiveMass, 1.0, $Temperature, 1);
	$Material->CalEFFromNe($Ne, $Temperature, 1);

	my $ND = $DonorDensity;
	my $NA = $DonorDensity;
	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 );
	printf("Depletion layer thickness: %12.6g nm\n", $W * 1e9);
	my $xn = $W / $Ntot * $NA;
	my $xp = $W / $Ntot * $ND;
	printf("  n: %12.6g nm\n", $xn * 1e9);
	printf("  p: %12.6g nm\n", $xp * 1e9);
	my $Enmax = $e * $ND / $nDielectricConstant * $xn;
	my $Epmax = $e * $NA / $pDielectricConstant * $xp;
	printf("Maxium electric field:\n");
	printf("  n: %12.6g V/cm\n", $Enmax * 1.0e-2);
	printf("  p: %12.6g V/cm\n", $Epmax * 1.0e-2);
}

sub ElectronEmission
{
	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();

	my $RichardsonDushman = 4.0*$pi * $me * $e * $kB * $kB / $h / $h / $h;
	printf("Richardson Constant: %e  A/m2/K2\n", $RichardsonDushman);
	my $mef = $ElectronEffectiveMass / $me;
	my $RichardsonDushmanef = $RichardsonDushman * $mef;
	printf("Effective Ricahrdson constant=%6.3fme: %e  A/m2/K2\n", 
			$mef, $RichardsonDushmanef);

	my $D0 = 1.0 / 2.0 / $pi / $pi * (2.0 * $ElectronEffectiveMass / $hbar / $hbar)**(1.5);
			# D(E) = $D * sqrt(E-EC)
	my $EF = ($MetalNe / $D0 * 1.5)**(2.0/3.0);
	printf("D(E) = %12.6g*sqrt(E-EC)\n", $D0);
	printf("EF(metal): %6.3f eV\n", $EF / $e);

	printf("Work function: %8.4g eV\n", $WorkFunction / $e);
	my $J0 = $RichardsonDushmanef * $Temperature * $Temperature
				* exp( -$WorkFunction / $kB / $Temperature );
	printf("Thermionic current(No Schottky): J0 = %e A/cm2\n", $J0 * 1.0e-4);

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

# 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);
	
	open(OUT,">$OutputFile") or die "$!\n";

	print "V(V),E(V/cm),JSchottky(A/cm2),JSCLC(A/cm2),JFN(A/cm2)\n";
	print OUT "V(V),E(V/cm),JSchottky(A/cm2),JSCLC(A/cm2),JFN(A/cm2),JFNCorrected(A/cm2),y,v,t\n";
	for(my $i = 0 ; $i < $nV ; $i++) {
		my $V = $VStart + $i * $VStep;
		my $E = $V / $Thickness;

		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;

		printf("%12.6g,%12.6g,%12.6g,%12.6g,%12.6g\n", 
			$V, $E * 1.0e-2, $JSchottky*1.0e-4, $JSCLC*1.0e-4, $JFN*1.0e-4);
		printf(OUT "%12.6g,%12.6g,%12.6g,%12.6g,%12.6g,%g,%g,%g\n",
			$V, $E * 1.0e-2, $JSchottky*1.0e-4, $JSCLC*1.0e-4, $JFN*1.0e-4, $JFNCorrected*1.0e-4, $y, $v, $t);
	}

	close(OUT);
}


sub CalSemiconductorParameters
{
	my $ElectronMobility = $e * $MomentumRelaxationTime / $ElectronEffectiveMass;
	my $sigma = $e * $Ne * $ElectronMobility;
	printf("Electron:\n");
	printf("  Effective mass : %g me\n", $ElectronEffectiveMass / $me);
	printf("  Relaxation time: %g s\n",  $MomentumRelaxationTime);
	printf("  Density     : %g cm-3\n",    $Ne * 1.0e-6);
	printf("  Mobility    : %g cm2/Vs\n",  $ElectronMobility * 1.0e4);
	printf("  Conductivity: %g S/cm\n",    $sigma * 1.0e-2);
	
	my $vth = sqrt(3.0 * $kB * $Temperature / $ElectronEffectiveMass);
	printf("  Thermal velocity: %e m/s\n", $vth);

	my $E = $AppliedVoltage / $Thickness;
	my $vdrift = $ElectronMobility * $E;
	printf("  Drift velocity  : %e m/s at E=%g V/m\n", $vdrift, $E);
	
	my $lth = $vth * $MomentumRelaxationTime;
	printf("  Mean free path(vth): %g nm\n", $lth * 1.0e9);

	my $ldrift = $vdrift * $MomentumRelaxationTime;
	printf("  Mean free path(vdrift): %g nm\n", $ldrift * 1.0e9);

	my $ElectronDiffusionConstant = $ElectronMobility / $e * $kB * $Temperature;
	printf("  Diffusion constant: %g cm2/s\n", $ElectronDiffusionConstant * 1.0e4);
	my $ElectronDiffusionLength = sqrt($ElectronDiffusionConstant * $ElectronLifeTime);
	printf("  Diffusion length  : %g nm for lifetime of $ElectronLifeTime s\n", 
				$ElectronDiffusionLength * 1.0e9, $ElectronLifeTime);

	my $D0 = 1.0 / 2.0 / $pi / $pi * (2.0 * $ElectronEffectiveMass / $hbar / $hbar)**(1.5);
			# D(E) = $D * sqrt(E-EC)
#$Ne=1e21 * 1.0e6;
	my $EF = ($Ne / $D0 * 1.5)**(2.0/3.0);
	printf("  Dc(E-EC) = %12.6g*sqrt(E-EC)\n", $D0);
	printf("  EF-EC: %12.6g eV\n", $EF / $e);
}

sub IonHopping1D
{
#	my $eFactor0 =$U / $kB / $Temperature;
	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 = $NIon * $e * $e * $dHopping * $dHopping * $fIon / $kB / $Temperature * $eFactor;
	my $sigmaOhmic = $K * $eK;
	printf("sigma(ohmic) = %g S/cm\n", $sigmaOhmic * 1.0e-2);

	my $mobilitySigma = $sigmaOhmic / $NIon / $e;
	printf("Mobility(sigma)         : %g cm2/Vs\n", $mobilitySigma * 1.0e4);
	my $mobilityEinstein = $e * $dHopping * $dHopping * $fIon / $kB / $Temperature;
	printf("Mobility(Einstein's law): %g cm2/Vs\n", $mobilityEinstein * 1.0e4);
	
	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);

	open(OUT,">$OutputFile") or die "$!\n";

	print "V(V),E(V/cm),JHopping(A/cm2),JPF(A/cm2)\n";
	print OUT "V(V),E(V/cm),JHopping(A/cm2),JPF(A/cm2)\n";
	for(my $i = 0 ; $i < $nV ; $i++) {
		my $V = $VStart + $i * $VStep;
		my $E = $V / $Thickness;

		my $J = $K * Sci::sinh($eK * $E);
		my $JPF = $J * exp( $PFSlopekT * sqrt(abs($E)) );

		printf("%12.6g,%12.6g,%12.6g,%12.6g\n", $V, $E * 1.0e-2, $J*1.0e-4, $JPF*1.0e-4);
		printf(OUT "%12.6g,%12.6g,%12.6g,%12.6g\n", $V, $E * 1.0e-2, $J*1.0e-4, $JPF*1.0e-4);
#print "e=$eFactor0,$eFactor,$K,$eK\n";
	}

	close(OUT);
}