#!/usr/bin/perl

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

use strict;
use Sci qw($kB $e);
use Sci::Material;
use Sci::TFT;
use Sci::Optimize;

#================================
# Global parameters
#================================
my $ElementsCSV = "ElementsIV.csv";
my $CircuitCSV  = "CircuitIV.csv";

my $material = new Material;

#================================
# Diode parameters
#================================
my $T  = 300.0;
my $S  = 60e-6 * 20e-6;
my $I0 = 1.0e-30 * $S;
my $nd = 1.0;
my $Rs = 1.0e5;
my $kd = $e / ($nd*$kB*$T);

#================================
# TFT parameters
#================================
my $TFTModel = 'GCA';
#my $TFTModel = 'Exact';
my $TFTW     = 200e-6;
my $TFTL     =  50e-6;
my $TFTNes   = 1.0e15; # cm-3
my $TFTEpss  = 13.0;
my $TFTuFE   = 10.0;   # cm2/Vs
my $TFTEpsg  = 4.0;
my $TFTdg    = 150.0;  # nm
my $Vth      = 0.0;    # V
my $VFB      = 0.0;    # V
my $VB       = 0.0;    # V
my $TFTCox   = $material->CalCapacitance($TFTEpsg, 0.01*0.01, $TFTdg*1.0e-9); # F/cm2
print "Cox: $TFTCox F/cm2\n";

my $tft = new TFT;
$tft->SetParameters(
	T     => $T,
	W     => $TFTW,
	L     => $TFTL,
	Nes   => $TFTNes,
	Epss  => $TFTEpss,
	uFE   => $TFTuFE,
	Cox   => $TFTCox,
	Vth   => $Vth,
	VFB   => $VFB,
	VB    => $VB,
	Model => $TFTModel,
	);
$tft->PrepareForCalculation();

#================================
# Circuit parameters
#================================
#my $Circuit = "pTFTNormal";
my $Circuit = "pTFTNormal2";
#my $Circuit = "nTFTNormal";
#my $Circuit = "nTFTInverted";

my $VGND = 0.0;
my $Vdd  = 10.0;

#================================
# Optimization parameters
#================================
#my $Method = "PDL::Simplex";
my $Method = "Amoeba::Simplex";
#my $Method = "ModifiedNewton";
#my $Method = "LinearOptimization"; # 線形パラメータのみ。

my $EPS               = 1.0e-4;
my $IEPS              = 1.0e-12;
my $NewtonDump        = 0.0;
my $NewtonDiffV1      = 0.001;
my $nMaxIter          = 300;
my $BisectionnMaxIter = 5;
my $BisectionV1EPS    = 0.01;
my $iPrintLevel       = -1;

#================================
# Main
#================================
print "Optimization method: $Method\n";
print "Circuit: $Circuit\n";
print "\n";

print "Open $ElementsCSV\n";
my $out  = JFile->new($ElementsCSV, "w") or die "$!: Can not write to [$ElementsCSV].\n";
print "Open $CircuitCSV\n";
my $cout = JFile->new($CircuitCSV, "w") or die "$!: Can not write to [$CircuitCSV].\n";

$out->print("\n");
print("Diode I-V from V\n");
$out->print("Diode I-V from V\n");
$out->print("V,I,I(Rs=0)\n");
my $Iprev;
for(my $V = 0.0 ; $V <= 5.0 ; $V += 0.2) {
	my ($I  , $Ir,   $V1O)   = &DiodeRsI($V, 0);
	my ($I00, $Ir00, $V1O00) = &DiodeRsI($V, 1);
	$out->print("$V,$I,$I00\n");
	$Iprev = $I;
}

$out->print("\n");
print("TFT I-V\n");
$out->print("TFT I-V\n");
$out->print("Vds");
for(my $Vgs = 2.0 ; $Vgs <= 20.0 ; $Vgs += 2.0) {
	$out->print(",Ids(Vgs=$Vgs)");
}
$out->print("\n");
for(my $Vds = 0.0 ; $Vds <= 20.0 ; $Vds += 1.0) {
	$out->printf("$Vds");
	for(my $Vgs = 2.0 ; $Vgs <= 20.0 ; $Vgs += 2.0) {
		my $Ids = &nTFTI($Vgs, $Vds);
		$out->print(",$Ids");
	}
	$out->print("\n");
}

print("Circuit I-V (Vdd=$Vdd)\n");

my $V1  = $Vdd / 2.0;
if($Circuit eq "pTFTNormal") {
	$cout->print("Vg,I(OLED),I(TFT),V1,V(OLED),Vgs(TFT),Vds(TFT)\n");
	for(my $Vg = 0.0 ; $Vg <= 20.0 ; $Vg += 1.0) {
		my ($IOLED, $Itft, $V1, $VOLED, $Vgs, $Vds) = &pTFTCircuitI(\&GetPotentialsForpTFTNormal, $Vdd, $Vg, $V1);
		$cout->printf("$Vg,$IOLED,$Itft,$V1, $VOLED, $Vgs, $Vds\n");
	}
}
elsif($Circuit eq "pTFTNormal2") {
	$cout->print("Vg,I(OLED),I(TFT),V1,V(OLED),Vgs(TFT),Vds(TFT)\n");
	for(my $Vg = 0.0 ; $Vg <= 20.0 ; $Vg += 1.0) {
		my ($IOLED, $Itft, $V1, $VOLED, $Vgs, $Vds) = &pTFTCircuitI(\&GetPotentialsForpTFTNormal2, -$Vdd, -$Vg, -$V1);
		$cout->printf("$Vg,$IOLED,$Itft,$V1, $VOLED, $Vgs, $Vds\n");
	}
}
elsif($Circuit eq "nTFTNormal") {
	$cout->print("Vg,I(OLED),I(TFT),V1,V(OLED),Vgs(TFT),Vds(TFT)\n");
	for(my $Vg = 0.0 ; $Vg <= 20.0 ; $Vg += 1.0) {
		my ($IOLED, $Itft, $V1, $VOLED, $Vgs, $Vds) = &nTFTCircuitI(\&GetPotentialsFornTFTNormal, $Vdd, $Vg, $V1);
		$cout->printf("$Vg,$IOLED,$Itft,$V1, $VOLED, $Vgs, $Vds\n");
	}
}
elsif($Circuit eq "nTFTInverted") {
	$cout->print("Vg,I(OLED),I(TFT),V1,V(OLED),Vgs(TFT),Vds(TFT)\n");
	for(my $Vg = 0.0 ; $Vg <= 20.0 ; $Vg += 1.0) {
		my ($IOLED, $Itft, $V1, $VOLED, $Vgs, $Vds) = &nTFTCircuitI(\&GetPotentialsFornTFTInverted, $Vdd, $Vg, $V1);
		$cout->printf("$Vg,$IOLED,$Itft,$V1, $VOLED, $Vgs, $Vds\n");
	}
}

$cout->Close();
$out->Close();

exit;

#==============================================
# Subroutines
#==============================================
sub GetOptimize
{
	my $optimize = new Optimize;
#	$optimize->SetDumpingFactor($DumpingFactor, $DumpingFactor, 1);
#	$optimize->SetDiffV($DiffV);

	return $optimize;
}

sub GetPotentialsFornTFTInverted
{
	my ($VGND, $Vdd, $Vg, $V1) = @_;

	my $Vgs   = $Vg  - $VGND;
	my $Vds   = $V1  - $VGND;
	my $VOLED = $Vdd - $V1;

	return ($VOLED, $Vds, $Vgs);
}

sub GetPotentialsForpTFTNormal
{
	my ($VGND, $Vdd, $Vg, $V1) = @_;

	my $Vgs   = $Vg  - $Vdd;
	my $Vds   = $Vdd - $V1;
	my $VOLED = $V1  - $VGND;

	return ($VOLED, $Vds, $Vgs);
}

sub GetPotentialsForpTFTNormal2
{
	my ($VGND, $Vdd, $Vg, $V1) = @_;

	my $Vgs   = $Vg   - $VGND;
	my $Vds   = $VGND - $V1;
	my $VOLED = $V1   - $Vdd;

	return ($VOLED, $Vds, $Vgs);
}

sub GetPotentialsFornTFTNormal
{
	my ($VGND, $Vdd, $Vg, $V1) = @_;

	my $Vgs   = $Vg  - $V1;
	my $Vds   = $Vdd - $V1;
	my $VOLED = $V1  - $VGND;

	return ($VOLED, $Vds, $Vgs);
}

sub nTFTCircuitI
{
	my ($pGetPotentialFunc, $Vdd, $Vg, $V1) = @_;

	my $optimize = &GetOptimize();
	$optimize->AddParameters(
#name, var, ID, scale, min, max
		"V1",  \$V1, 1, 0.1, 0.0, undef,
		sub {
			$V1 = $_[0];
		},
		);

	my ($OptVars, $MinVal) = $optimize->Optimize(
		$Method, undef, undef, undef,
		$EPS, $nMaxIter, $iPrintLevel,
		sub { &CalSFornTFTCurcuitI($pGetPotentialFunc, $Vdd, $Vg, @_); },
		undef,
		sub { Optimize::BuildDifferentialMatrixes(@_); },
		);
	$optimize->RecoverParameters($OptVars);
if($EPS < $MinVal or $MinVal < 0) {
	print "Error in CircuitI: Not converged for Vdd = $Vdd, Vg = $Vg [V1=$V1, S=$MinVal (EPS=$EPS)\n";
	exit;
}
#	print "\nOptimized at I = $I for V = $V (S = $MinVal):\n";
#	$optimize->PrintParameters(1, $OptVars, $MinVal, 1);

	my ($VOLED, $Vds, $Vgs) = &$pGetPotentialFunc($VGND, $Vdd, $Vg, $V1);

	my $Itft  = &nTFTI($Vgs, $Vds);
	my ($IOLED, $Ir, $V1O) = &DiodeRsI($VOLED, 0);
printf "Vg=%5.2f VOLED=%5.2f Vds=%5.2f Vgs=%5.2f (%5.2f - %5.2f - %5.2f) Itft=%8.3e IOLED=%8.3e\n",
		$Vg, $VOLED, $Vds, $Vgs, $VGND, $V1, $Vdd, $Itft, $IOLED;

	return ($IOLED, $Itft, $V1, $VOLED, $Vgs, $Vds);
}

sub CalSFornTFTCurcuitI {
	my ($pGetPotentialFunc, $Vdd, $Vg, $pVars, $iPrintLevel) = @_;
	my ($V1) = @$pVars;
#print "V,I=$V,$I\n";

	my ($VOLED, $Vds, $Vgs) = &$pGetPotentialFunc($VGND, $Vdd, $Vg, $V1);

	my $Itft  = &nTFTI($Vgs, $Vds);
	my ($IOLED, $Ir, $V1O) = &DiodeRsI($VOLED, 0);
#print "VOLED=$VOLED\tVg=$Vg\tVds=$Vds\tVgs=$Vgs ($VGND - $V1 - $Vdd) Itft=$Itft\tIOLED=$IOLED\n";

	my $S   = ($Itft - $IOLED)**2;
#printf("a,b,c=%6.3f\t%6.3f\t%6.3f\t(%g)\n", $a, $b, $c, $S) if($iPrintLevel >= 2);

	return $S;
}

sub pTFTCircuitI
{
	my ($pGetPotentialFunc, $Vdd, $Vg, $V1) = @_;

	my $optimize = &GetOptimize();
	$optimize->AddParameters(
#name, var, ID, scale, min, max
		"V1",  \$V1, 1, 0.1, 0.0, undef,
		sub {
			$V1 = $_[0];
		},
		);

	my ($OptVars, $MinVal) = $optimize->Optimize(
		$Method, undef, undef, undef,
		$EPS, $nMaxIter, $iPrintLevel,
		sub { &CalSForpTFTCurcuitI($pGetPotentialFunc, $Vdd, $Vg, @_); },
		undef,
		sub { Optimize::BuildDifferentialMatrixes(@_); },
		);
	$optimize->RecoverParameters($OptVars);
if($EPS < $MinVal or $MinVal < 0) {
	print "Error in CircuitI: Not converged for Vdd = $Vdd, Vg = $Vg [V1=$V1, S=$MinVal (EPS=$EPS)\n";
	exit;
}
#	print "\nOptimized at I = $I for V = $V (S = $MinVal):\n";
#	$optimize->PrintParameters(1, $OptVars, $MinVal, 1);

	my ($VOLED, $Vds, $Vgs) = &$pGetPotentialFunc($VGND, $Vdd, $Vg, $V1);

	my $Itft  = &pTFTI($Vgs, $Vds);
	my ($IOLED, $Ir, $V1O) = &DiodeRsI($VOLED, 0);
printf "Vg=%5.2f VOLED=%5.2f Vds=%5.2f Vgs=%5.2f (%5.2f - %5.2f - %5.2f) Itft=%8.3e IOLED=%8.3e\n",
		$Vg, $VOLED, $Vds, $Vgs, $VGND, $V1, $Vdd, $Itft, $IOLED;

	return ($IOLED, $Itft, $V1, $VOLED, -$Vgs, $Vds);
}

sub CalSForpTFTCurcuitI 
{
	my ($pGetPotentialFunc, $Vdd, $Vg, $pVars, $iPrintLevel) = @_;
	my ($V1) = @$pVars;
#print "V,I=$V,$I\n";

	my ($VOLED, $Vds, $Vgs) = &$pGetPotentialFunc($VGND, $Vdd, $Vg, $V1);

	my $Itft               = &pTFTI($Vgs, $Vds);
	my ($IOLED, $Ir, $V1O) = &DiodeRsI($VOLED, 0);
#print "VOLED=$VOLED\tVg=$Vg\tVds=$Vds\tVgs=$Vgs ($VGND - $V1 - $Vdd)\tItft=$Itft\tIOLED=$IOLED\n";

	my $S   = ($Itft - $IOLED)**2;
#printf("a,b,c=%6.3f\t%6.3f\t%6.3f\t(%g)\n", $a, $b, $c, $S) if($iPrintLevel >= 2);

	return $S;
}

sub nTFTI
{
	my ($Vgs, $Vds) = @_;
	return $tft->IdsByModel($Vgs, $Vds);
}

sub pTFTI
{
	my ($Vgs, $Vds) = @_;
	return $tft->IdsByModel(-$Vgs, $Vds);
}

sub DiodeI
{
	my ($V) = @_;
	return $I0 * (exp($kd * $V) - 1.0);
}

sub DiodeV
{
	my ($I) = @_;
	return -1e10 if($I <= -$I0);
	return log($I / $I0 + 1.0) / $kd;
}


sub DiodeRsV
{
	my ($I, $IgnoreRs) = @_;

	my $Vd = &DiodeV($I);

	return $Vd if($IgnoreRs);
	return $Rs * $I + $Vd;
}

sub GetPotentialsForDiodeRs
{
	my ($VGND, $V, $V1) = @_;

	my $Vdiode = $V1 - $VGND;
	my $Vr     = $V  - $V1;

	return ($Vdiode, $Vr);
}

sub DiodeRsI
{
	my ($V, $IgnoreRs, $V1) = @_;

#	return 0.0 if($V == 0.0);

	my $I = &DiodeI($V);
	return $I  if($IgnoreRs);
#	return $I  if($IgnoreRs or $V <= 0.0);

	my ($pVars, $Smin, $nIter) = Optimize->BisectionMethod(0.0, $V, 
					sub { &CalDiffForDiodeRsI($V, @_); }, 
					$BisectionnMaxIter, $BisectionV1EPS, $IEPS, $iPrintLevel);
	$V1 = $pVars->[0];
#print "Bisection: V1=$V1   S=$Smin  nIter=$nIter\n";
	($pVars, $Smin, $nIter) = Optimize->SolveByModifiedNewton1D($V1, $NewtonDiffV1, $NewtonDump,
					sub { &CalDiffForDiodeRsI($V, @_); }, 
					$nMaxIter, $EPS, $IEPS, $iPrintLevel);
	$V1 = $pVars->[0];
#print "ModifiedNewton: V1=$V1   S=$Smin  nIter=$nIter\n";

	my ($Vdiode, $Vr) = &GetPotentialsForDiodeRs($VGND, $V, $V1);

	my $Idiode = &DiodeI($Vdiode);
	my $Ir     = $Vr / $Rs;

	return ($Idiode, $Ir, $V1);
}

sub CalDiffForDiodeRsI
{
	my ($V, $V1, $iPrintLevel) = @_;

#	return 0.0 if($V == 0.0);

	my ($Vdiode, $Vr) = &GetPotentialsForDiodeRs($VGND, $V, $V1);

	my $Idiode = &DiodeI($Vdiode);
	my $Ir     = $Vr / $Rs;

	my $S = $Idiode - $Ir;
#print "V=$V  V1=$V1  S=$S = $Idiode - $Ir\n";
	return $S;
}