#!/usr/bin/perl

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

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'PerlDir'};
$BaseDir = $ENV{'TkPerlDir'} unless(defined $BaseDir);
#print "\n\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", @INC);
}

use strict;
#use warnings;

use Deps;
use Utils;
use JFile;

use MyApplication;
use Sci::Science;
use Sci::Optimize;

my $e0 = Sci::e0();
my $e  = Sci::e();

my $IsPrint = 0;
my $IsSave  = 1;

my $Ne         = 1.0e12; # cm-3
my $er         = 3.0;
my $mobility   = 1.0e-4;   # cm2/Vs
my $thickness  = 200.0 * 1e-7; # cm
my $ndMesh     = 1001;
my $dthickness = $thickness / ($ndMesh - 1) * 0.01; # m

my $I          = 2.0;          # A/m2
my $Istart     =  0.0 * 1.0e4; # A/m2
my $Iend       =  3.0 * 1.0e4;
my $nI         = 101;
my $Istep      = ($Iend - $Istart) / ($nI - 1);

print "I: $Istart - $Iend, $Istep ($nI)\n";
print "x: 0 - $thickness, $dthickness ($ndMesh)\n";

my $Eohmic = $I / $e / ($Ne*1e6) / ($mobility*1e-4);
my $Vohmic = $Eohmic * ($thickness * 1.0e-2);

print "d=", $thickness * 1.0e-2, " m\n";
print "I0=$I A/m2\n";
print "Eohmic=$Eohmic V/m\n";
print "Vohmic=$Vohmic V\n";

my $E0 = $Eohmic;
my ($pX, $pE, $pV, $prho, $pV2) = CalSCLC($E0, $I, $IsPrint);
my $Vohmic = $pV->[$ndMesh-1];
&Save("OhmicLimit.csv", $pX, $pE, $pV, $prho, $pV2) if($IsSave);

($pX, $pE, $pV, $prho, $pV2) = CalSCLC(0.0, $I, $IsPrint);
my $VSCLC = $pV->[$ndMesh-1];
&Save("SCLCLimit.csv", $pX, $pE, $pV, $prho, $pV2) if($IsSave);

print "\n";
printf "Vohmic=%12.6g V\n", $Vohmic;
printf "VSCLC =%12.6g V\n", $VSCLC;

my $sigma = $e * ($Ne * 1.0e6) * ($mobility * 1.0e-4);
my $KSCLC = 9.0 / 8.0 * ($mobility * 1.0e-4) * ($e0 * $er) / ($thickness * 1.0e-2);

my $Method = "Amoeba::Simplex";
#my $Method = "ModifiedNewton";
my $EPS = 1.0e-9;
my $nMaxIter = 300;
my @Guess = ($Eohmic);
my @OptId = (1);
my @Scale = ($Eohmic);
my $iPrintLevel = 0;
my $Opt = new Optimize;

open(OUT, ">IV.csv") or die "$!\n";
print OUT "V(V),I(A/m2),Iohm(A/m2),ISCLC(A/m2),E0(V/m),Pohmic(%)\n";
for(my $i = 0 ; $i < $nI ; $i++) {
	my $I = $Istart + $i * $Istep;
	my $Eohmic = $I / $e / ($Ne*1e6) / ($mobility*1e-4);
#print "I=$I A  Eohmic=$Eohmic V/m\n";
	@Guess = ($Eohmic);

	my ($OptVars, $MinVal);
	if($I == 0.0) {
		$OptVars = [0.0];
		$MinVal  = 0.0;
	}
	else {
		($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				sub { 
					my ($pVars, $iPrintLevel) = @_;
					my ($pX, $pE, $pV, $prho, $pV2) = &CalSCLC($pVars->[0], $I, 0);
					my $V = $pV->[$ndMesh-1];
					if($V < 0.0) {
						return 1.0e100;
					}
					else {
						return $V;
					}
				},
				sub {},
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
	}
	my $Pohmic;
	if($I == 0.0) {
		$Pohmic = 100.0;
	}
	else {
		$Pohmic = $OptVars->[0] / $Eohmic * 100.0;
	}
	
	my $E0   = $OptVars->[0];
	my $V    = $MinVal;
	my $E    = $V / ($thickness * 1.0e-2);
	my $Iohm = $sigma * $E;
	my $ISCLC = $KSCLC * $E * $E;	
	printf "I=%8.3g A/m2  Emin=%12.6g V/m  Vmin=%12.6g V  Pohmic=%5.1f %%\n", 
			$I, $OptVars->[0], $MinVal, $Pohmic;
	print OUT "$V,$I,$Iohm,$ISCLC,$E0,$Pohmic\n";
}
close(OUT);
exit;

sub CalSCLC
{
	my ($E0, $I, $IsPrint) = @_;

	my $Eohmic = $I / $e / ($Ne*1e6) / ($mobility*1e-4);
	my $eu     = $er * $e0 * ($mobility * 1.0e-4);
	my $eNee   = $e * ($Ne * 1.0e6) / $er / $e0;

	my @x   = (0.0);
	my @V   = (0.0); # Potential
	my @V2  = (0.0); # d2V/dx2
	my @E   = ($E0); # dV/dx
	my @rho = (0.0); # Charge density: d2V/dx2 = rho / (er*e0)

if($IsPrint) {
	print "x(m)\tE(V/m)\tV(V)\trho(cm-3)\n";
}
	for(my $i = 1 ; $i < $ndMesh ; $i++) {
		my $ip = $i - 1;
		$x[$i] = $i * $dthickness;
		my $Eavg;
		if($E[$ip] == 0.0) {
			$Eavg = sqrt(2.0 * $I / $eu * ($x[$i] + $x[$ip]) / 2.0);
		}
		else {
			$Eavg = $E[$ip];
			if($i > 1) {
				$Eavg = $E[$ip] + ($E[$ip] - $E[$ip-1]) / 2.0;
			}
		}
		$V2[$i]  = -$eNee * (-$Eohmic / $Eavg + 1.0);
		$E[$i]   = $E[$ip] + ($V2[$i]+$V2[$ip])/2.0 * $dthickness;
		$V[$i]   = $V[$ip] + ($E[$i]+$E[$ip])/2.0 * $dthickness;
		$rho[$i] = $V2[$i] * $er * $e0 / $e * 1.0e-6;

if($IsPrint) {
		printf("%8.4g\t%12.6g\t%12.6g\t%12.6g\n", $x[$i], $E[$i], $V[$i], $rho[$i]);
}
	}
	return (\@x, \@E, \@V, \@rho, \@V2);
}

sub Save
{
	my ($OutFile, $pX, $pE, $pV, $prho, $pV2) = @_;

	open(OUT, ">$OutFile") or die "$!\n";
	print OUT "x(m),E(V/m),V(V),rho(cm-3),d2V/dx2\n";
	printf(OUT "%8.4g,%12.6g,%12.6g,%12.6g,%12.6g\n", 
			$pX->[0], $pE->[0], $pV->[0], $prho->[0], $pV2->[0]);
	for(my $i = 1 ; $i < $ndMesh ; $i++) {
		printf(OUT "%8.4g,%12.6g,%12.6g,%12.6g,%12.6g\n", 
			$pX->[$i], $pE->[$i], $pV->[$i], $prho->[$i], $pV2->[$i]);
	}
	close(OUT);
}