package DiodeRs;
use Exporter;
@ISA = qw(Exporter);

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

use strict;
use File::Basename;

use Utils;
use Sci qw($pi $kB $R $NA $c $e $e0 $u0 $me $mp $mn $h $hbar $torad $todeg);
use Sci::Optimize;

#============================================================
#　変数等取得、再定義関数
#============================================================
sub SetTemperature
{
	my ($this, $T) = @_;
	$this->{T} = $T;
	$this->{ekT} = $e / ($kB * $T);
}

#============================================================
#　コンストラクタ、デストラクタ
#============================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	$this->SetTemperature(300.0);

	return $this;
}

sub DESTROY { my $this = shift; }

#============================================================
#　一般関数
#============================================================
sub SetParameters
{
	my ($this, %args) = @_;
	
	Utils::MergeHash($this, \%args);
}

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

	$this->{I0} = 1.0   if(!defined $this->{I0});
	$this->{nd} = 1.0   if(!defined $this->{nd});
	$this->{T}  = 300.0 if(!defined $this->{T});

	$this->{Rs} = 0.0   if(!defined $this->{Rs});

	$this->{BisectionV1EPS}    = 0.01    if(!defined $this->{BisectionV1EPS});
	$this->{IEPS}              = 1.0e-10 if(!defined $this->{IEPS});
	$this->{BisectionnMaxIter} = 5       if(!defined $this->{BisectionnMaxIter});

	$this->{NewtonDiffV1}      = 0.001   if(!defined $this->{NewtonDiffV1});
	$this->{NewtonDump}        = 0.0     if(!defined $this->{NewtonDump});
	$this->{V1EPS}             = 1.0e-4  if(!defined $this->{V1EPS});
	$this->{nMaxIter}          = 30      if(!defined $this->{nMaxIter});

	$this->{iPrintLevel}       = -2      if(!defined $this->{iPrintLevel});

	$this->{kd} = $e / ($this->{nd} * $kB * $this->{T});
}

sub RsI
{
	my ($this, $V) = @_;
	return $V / $this->{Rs};
}

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

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

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

	my $I = $this->DiodeI($V);
	return $I if($IgnoreRs or $this->{Rs} == 0.0);

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

	my ($Vdiode, $Vr) = $this->GetPotentialsForDiodeRs(0.0, $V, $V1);

	my $Idiode = $this->DiodeI($Vdiode);
	my $Ir     = $this->RsI($Vr);

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

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

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

	return ($Vdiode, $Vr);
}

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

	my ($Vdiode, $Vr) = $this->GetPotentialsForDiodeRs(0.0, $V, $V1);
	my $Idiode = $this->DiodeI($Vdiode);
	my $Ir     = $this->RsI($Vr);
	my $S      = $Idiode - $Ir;

	return $S;
}

#================================
# obsolete
#================================
sub CalSCFTotalIV
{
	my ($this, $Vtarget, $J0, $n, $Rs, $LSQMethod, $EPS, $nMaxIter, $iPrintLevel) = @_;
	$LSQMethod   = "ModifiedNewton" if(!defined $LSQMethod);
	$EPS         = 1.0e-4           if(!defined $EPS);
	$nMaxIter    = 100              if(!defined $nMaxIter);
	$iPrintLevel = 0                if(!defined $iPrintLevel);

	if($Rs == 0.0) {
		my ($Vtotal, $Jtotal) = $this->CalDiodeIV($Vtarget, $J0, $n, $Rs);
		return ($Vtotal, $Jtotal);
	}

	my $Opt    = new Optimize;
	my @Guess  = ($Vtarget);
	my @OptId  = (     1);
	my @Scale  = (   0.1);
	my ($OptVars, $MinVal) = $Opt->Optimize($LSQMethod, \@Guess, \@OptId, \@Scale, 
		$EPS, $nMaxIter, $iPrintLevel,
		sub { $this->MinimizationFuncIVDiode($Vtarget, $J0, $n, $Rs, @_); }, 
		sub { },
		sub { Optimize::BuildDifferentialMatrixes(@_); },
		);
	my ($Vdiode) = @$OptVars;
if($MinVal > $EPS) {
	print "Error i Cal: Could not get convergence for V=$Vtarget (MinVal=$MinVal)\n";
#	exit;
}
	my ($Vtotal, $Jtotal) = $this->CalDiodeIV($Vdiode, $J0, $n, $Rs);
	return ($Vdiode, $Jtotal);
}

sub CalDiodeIV
{
	my ($this, $Vdiode, $J0, $n, $Rs) = @_;
	my $Jdiode = $J0 * (exp($this->{ekT} * $Vdiode / $n) - 1.0);
	my $Vtotal = $Vdiode + $Rs * $Jdiode;
	return ($Vtotal, $Jdiode);
}

sub MinimizationFuncIVDiode {
	my ($this, $Vtarget, $J0, $n, $Rs, $pVars, $iPrintLevel) = @_;

	my $S = 0.0;
	my ($Vdiode) = @$pVars;
	my ($Vtotal, $Jtotal) = $this->CalDiodeIV($Vdiode, $J0, $n, $Rs);

	return ($Vtotal - $Vtarget)**2;
}

1;