#!/usr/bin/perl

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

use strict;

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

my $SourceCSV = "Lee-SC-IV.csv";
my $UseDark = 1;
my $Ilabel  = ($UseDark)? "Id.*" : "Ip.*";
my $F       = 100.0; # mW/cm2
my $T       = 300.0; # K
my $nSkip   = 5;
my ($FitVMin, $FitVMax) = (-0.5, 1.0);
my $nSimplexIter = 2;

my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($SourceCSV);
my $OutTXT    = Deps::MakePath("$drive$directory", "${filebody}-Out.txt", 0);
my $OutCSV    = Deps::MakePath("$drive$directory", "${filebody}-Out.csv", 0);
my $FitCSV    = Deps::MakePath("$drive$directory", "${filebody}-Fit.csv", 0);

my $Method = "Amoeba::Simplex";
#my $Method = "PDL::Simplex";
#my $Method = "ModifiedNewton";
#my $Method = "LinearOptimization"; # 線形パラメータのみ。
my $EPS = 1.0e-7;
my $nMaxIter = 3000;
my $iPrintLevel = 2;

my $txt = new JFile;
$txt->Open($OutTXT, "w") or die "$!: Can not write to [$OutTXT].\n";
$txt->print("Source   : $SourceCSV\n");
$txt->print("OutCSV   : $OutCSV\n");
$txt->print("FitCSV   : $FitCSV\n");
$txt->print("For dark?: $UseDark\n");
$txt->print("Current data taken from: $Ilabel\n");
$txt->print("Light intensity: $F mW/cm2\n");
$txt->print("Temperature    : $T K\n");
$txt->print("Fitting range  : $FitVMin - $FitVMax V, every $nSkip points\n");
$txt->print("nSimplexIter   : $nSimplexIter\n");

my $csv = new CSV;
if($csv->Read($SourceCSV)) {
	print("Read [$SourceCSV]\n");
}
else {
	print("  Error!!: Can not read [$SourceCSV]\n");
	exit;
}
my $pLabelArray = $csv->LabelArray();
print "Label: ", join(',', @$pLabelArray), "\n";
my $pDataArray  = $csv->DataArray();
my $pV = $csv->GetXData("V.*");
my $pI = $csv->GetYData($Ilabel);
my $pJ = $csv->GetYData($Ilabel);
if(!$pV) {
	print "Error: Can not find the data column V\n";
	exit;
}
if(!$pI) {
	print "Error: Can not find the data column I\n";
	exit;
}
if(!$pJ) {
	print "Error: Can not find the data column J\n";
	exit;
}
my $nData = @$pV;
if($nData == 0) {
	print "Error: nData = 0\n";
	exit;
}
my $Area = $pI->[0] / $pJ->[0];
my ($Vmin, $Vmax) = Utils::CalMinMax($pV);
my ($Imin, $Imax) = Utils::CalMinMax($pI);
print "pV=$pV\n";
print "pI=$pI\n";
print "nData = $nData\n";
print "Area  = $Area cm2\n";
print "\n";
$txt->print("\n");
$txt->print("Source data\n");
$txt->print("  nData=$nData\n");
$txt->print("  Electrode area: $Area cm2\n");

my $pW    = [];
my $pdJdV = [];
my ($Vop, $Jop);
my $Wmax = 0.0;
my ($Vsat, $Jsat);
my $dJdVmin = 1.0e10;
my $dJdVmaxAtForward = -1.0e10;
my $VforR;
my $dJdVminAtReverse = 1.0e10;
my $VrevR;
for(my $i = 0 ; $i < $nData ; $i++) {
	my $V = $pV->[$i];
	$pW->[$i] = -$V * $pJ->[$i] * 1.0e3;
	$pdJdV->[$i] = Algorism::DifferentiateByCubicSpline($pV, $pJ, $V);

	next if($i == $nData-1 or $i == 0);
	if($V >= 0.0) {
		if($dJdVmaxAtForward < $pdJdV->[$i]) {
			$dJdVmaxAtForward = $pdJdV->[$i];
			$VforR = $V;
		}
	}
	if($V >= 0.0 and $pJ->[$i] <= 0.0) {
		if($Wmax < $pW->[$i]) {
			$Wmax = $pW->[$i];
			$Vop  = $V;
			$Jop  = -$pJ->[$i] * 1.0e3; #mA/cm2
		}
	}
	if($V <= 0.0) {
		if($dJdVmin > $pdJdV->[$i]) {
			$dJdVmin = $pdJdV->[$i];
			$Vsat = $V;
		}
		if($dJdVminAtReverse > $pdJdV->[$i] and $pdJdV->[$i] >= 0.0) {
			$dJdVminAtReverse = $pdJdV->[$i];
			$VrevR = $V;
		}
	}
}
$Jsat = -$csv->YVal($pV, $pJ, $Vsat) * 1.0e3; #mA/cm2

$txt->print("\n");
$txt->print("Electric parameters\n");
my $Jsc = -$csv->YVal($pV, $pJ, 0.0) * 1.0e3; #mA/cm2
print "Jsc  = $Jsc A/cm2\n";
my $Voc = $csv->XValByBinaryMethod($pV, $pJ, 0.0);#, 1.0e-6, 300, 0.0, $Vmax, 0.01);
print "Voc  = $Voc mA/cm2\n";
print "Jsat = $Jsat mA/cm2 (at V=$Vsat V)\n";

print "Vop  = $Vop V\n";
print "Jop  = $Jop mA/cm2\n";
print "Wmax = $Wmax mW/cm2\n";
my $FF  = $Wmax / ($Jsc * $Voc);
print "FF   = $FF\n";
my $Eff = $Voc * $Jsc * $FF * 100.0 / $F;
print "Eff  = $Eff %\n";

my $RsOn  = 1.0 / $dJdVmaxAtForward; #$pdJdV->[$nData-2];
my $RsVoc = 1.0 / $csv->YVal($pV, $pdJdV, $Voc);
my $Rsh0  = 1.0 / $csv->YVal($pV, $pdJdV, 0.0);
my $RsRev = 1.0 / $dJdVminAtReverse;
print "Rs(on)   = $RsOn ohm cm2 (V=$VforR V)\n";
print "Rs(Voc)  = $RsVoc ohm cm2\n";
print "Rsh(0)   = $Rsh0 ohm cm2\n";
print "Rsh(rev) = $RsRev ohm cm2 (V=$VrevR V)\n";
print "\n";

$txt->print("Jsc  = $Jsc A/cm2\n");
$txt->print("Voc  = $Voc mA/cm2\n");
$txt->print("Jsat = $Jsat mA/cm2 (at V=$Vsat V)\n");
$txt->print("Vop  = $Vop V\n");
$txt->print("Jop  = $Jop mA/cm2\n");
$txt->print("Wmax = $Wmax mW/cm2\n");
$txt->print("FF   = $FF\n");
$txt->print("Eff  = $Eff %\n");
$txt->print("Rs(on)   = $RsOn ohm cm2 (V=$VforR V)\n");
$txt->print("Rs(Voc)  = $RsVoc ohm cm2\n");
$txt->print("Rsh(0)   = $Rsh0 ohm cm2\n");
$txt->print("Rsh(rev) = $RsRev ohm cm2 (V=$VrevR V)\n");
$txt->print("\n");

my $log10        = log(10.0);
my $Joffset      = 1.0e-10;
my $plogJ        = [];
my $pdlogJdV     = [];
my $pDiodeFactor = [];
my $dlogJdVmin   = 1.0e10;
my $Vdfactor;
my $dlogJdVdfactor;
my $ekT = $e / $kB / $T;
for(my $i = 0 ; $i < $nData ; $i++) {
	my $V = $pV->[$i];
	$plogJ->[$i] = log(abs($pJ->[$i] + $Jsc*1.0e-3 + $Joffset));# / $log10;
}
for(my $i = 0 ; $i < $nData ; $i++) {
	my $V = $pV->[$i];
	$pdlogJdV->[$i] = Algorism::DifferentiateByCubicSpline($pV, $plogJ, $V);
	if($V >= 0.0) {
		if($dlogJdVmin > $pdlogJdV->[$i]) {
			$dlogJdVmin = $pdlogJdV->[$i];
			$Vdfactor   = $V;
			$dlogJdVdfactor = $pdlogJdV->[$i];
		}
	}
	if($V > 0.0 and $pdlogJdV->[$i] > 0.0) {
		$pDiodeFactor->[$i] = $ekT * $V / $pdlogJdV->[$i];
	}
}
print "dlog(|J|)/dV(min)=$dlogJdVdfactor (V=$Vdfactor V)\n";

my $out = new JFile;
print("Write to [$FitCSV]\n");
$out->Open($OutCSV, "w") or die "$!: Can not write to [$FitCSV]\n";
$out->print("V(V),I(A),J(A/cm2),|J|(A/cm2),W(mW/cm2),dJ/dV(S/cm),dlog(|J|)/dV,n\n");
for(my $i = 0 ; $i < $nData ; $i++) {
	my $absJ = abs($pJ->[$i]);
	$out->print("$pV->[$i],$pI->[$i],$pJ->[$i],$absJ,$pW->[$i],$pdJdV->[$i],$pdlogJdV->[$i],$pDiodeFactor->[$i]\n");
}
$out->Close();

my $pV1    = [];
my $plogJ1 = [];
my $c = 0;
for(my $i= 0 ; $i < $nData ; $i++) {
	next if($pV->[$i] <= 0.0);
	$pV1->[$c]    = $pV->[$i];
	$plogJ1->[$c] = $plogJ->[$i];
#print "${i}[$c]: $pV1->[$c], $plogJ1->[$c]\n";
	$c++;
}

my $iPrintLevel = 0;
my $m = 2;
my ($pCi, $S2) = Optimize->new()->Optimize("mlsq", $pV1, $plogJ1, $m, $iPrintLevel);
print "C0=$pCi->[0]\n";
print "C1=$pCi->[1]\n";
print "C2=$pCi->[2]\n";
print "S2=$S2\n";

my $J01  = exp($pCi->[0]);
my $n1   = $ekT / $pCi->[1]; #2.0;
my $Rs  = $RsOn;
my $Rsh = $RsRev;
#$Rsh = 1e3;
my $Jpv = $Jsat*1.0e-3 - $J01;
$Jpv = 0 if($UseDark);
print "\n";
print "Intial values:\n";
print "J01 = $J01 A/cm2\n";
print "n1  = $n1\n";
print "Jpv = $Jpv A/cm2\n";
print "Rs  = $Rs ohm cm2\n";
print "Rsh = $Rsh ohm cm2\n";

$txt->print("\n");
$txt->print("Intial values:\n");
$txt->print("J01 = $J01 A/cm2\n");
$txt->print("n1  = $n1\n");
$txt->print("Jpv = $Jpv A/cm2\n");
$txt->print("Rs  = $Rs ohm cm2\n");
$txt->print("Rsh = $Rsh ohm cm2\n");

my ($OptVars, $MinVal);
my @VarLabel  = (  "J01",     "n1",    "Jpv",    "Rsh",    "Rs");
for(my $i = 0 ; $i < $nSimplexIter ; $i++) {
	my @Guess = (   $J01,      $n1,     $Jpv,     $Rsh,     $Rs);
	my @OptId = (      1,        1,        1,        1,       1);
$OptId[2] = 0 if($UseDark);
	my @Scale = ($J01*0.1, $n1*0.1, $Jpv*0.1, $Rsh*0.1, $Rs*0.1);
	my $Opt = new Optimize;
	($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, 2,#$iPrintLevel,
				sub { &MinimizationFunc(@_); }, \&PDLDisplayFunc,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
	($J01, $n1, $Jpv, $Rsh, $Rs) = @$OptVars;
	print "Optimized at (",join(',',@{$OptVars}),")=$MinVal\n";
}

$txt->print("\n");
$txt->print("Optimized parameters\n");
for(my $i = 0 ; $i < @VarLabel ; $i++) {
	$txt->print("$VarLabel[$i]: $OptVars->[$i]\n");
}
$txt->print("MinS: $MinVal\n");

$out->Open($FitCSV, "w") or die "$!: Can not write to [$FitCSV]\n";
$out->print("V(V),log(|Jcal|),log(|Jobs|),|Jcal|(A/cm2),|Jobs|(A/cm2),Jcal(A/cm2),Jobs(A/cm2),Vdiode(V)\n");
for(my $i = 0 ; $i < $nData ; $i++) {
	my $Vtarget = $pV->[$i];
	my ($Vdiode, $Jtotal) = &Cal($Vtarget, $J01, $n1, $Jpv, $Rsh, $Rs);
#print "$Vtarget, $Vdiode, $Jtotal\n";
	my $absJcal = abs($Jtotal);
	my $absJobs = abs($pJ->[$i]);
	my $logJcal = log($absJcal) / log(10.0);
	my $logJobs = log($absJobs) / log(10.0);
	if($i < @VarLabel) {
		$out->print("$Vtarget,$logJcal,$logJobs,$absJcal,$absJobs,$Jtotal,$pJ->[$i],$Vdiode,"
					."$VarLabel[$i]=,$OptVars->[$i]\n");
	}
	elsif($i == @VarLabel) {
		$out->print("$Vtarget,$logJcal,$logJobs,$absJcal,$absJobs,$Jtotal,$pJ->[$i],$Vdiode,"
					."MinS=,$MinVal\n");
	}
	else {
		$out->print("$Vtarget,$logJcal,$logJobs,$absJcal,$absJobs,$Jtotal,$pJ->[$i],$Vdiode\n");
	}
}
$out->Close();
$txt->Close();

exit;

#==============================================
# 最小化関数のサブルーチン
# MinimizationFunc    : 
# PDLDisplaySimplex: PDL::Opt::Simplex用のDisplay関数
#==============================================
sub MinimizationFunc {
	my ($pVars, $iPrintLevel) = @_;
	
	my ($J0, $n, $Jpv, $Rsh, $Rs) = @$pVars;
	my $S = 0.0;
	my $c = 0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		next if($i % $nSkip != 0);
		my $Vtarget = $pV->[$i];
		next if($Vtarget < $FitVMin or $FitVMax < $Vtarget);

		my ($Vdiode, $Jtotal) = &Cal($Vtarget, $J0, $n, $Jpv, $Rsh, $Rs);
		$S += (log(abs($Jtotal)) - log(abs($pJ->[$i])))**2 if($Jtotal != 0.0 and $pJ->[$i] != 0.0);
#print "$i: $Vtarget, $Jtotal, $pJ->[$i]\n";
		$c++;
	}
	$S += 1.0e5 if($J0  <= 0.0);
	$S += 1.0e5 if($n   <= 0.0);
	$S += 1.0e5 if($Rs  < 0.0);
	$S += 1.0e5 if($Rsh < 0.0);
	$S += 1.0e5 if($Jpv < 0.0);
printf("J0,n,Jpv,Rsh,Rs=%6.3g\t%6.3f\t%6.3g\t%6.3g\t%6.3g\t(%6.3g)\n", 
		$J0, $n, $Jpv, $Rsh, $Rs, $S) if($iPrintLevel >= 2);

	return $S;
}

sub Cal
{
	my ($Vtarget, $J0, $n, $Jpv, $Rsh, $Rs) = @_;

	my $Opt    = new Optimize;
	my $Method = "ModifiedNewton";
	my @Guess  = ($Vtarget);
	my @OptId  = (     1);
	my @Scale  = (   0.1);

	my ($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				sub { &IVMinimizationFunc($Vtarget, $J0, $n, $Jpv, $Rsh, $Rs, @_); }, \&PDLDisplayFunc,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
	my ($Vdiode) = @$OptVars;
if($MinVal > 1.0e-4) {
	print "Error i Cal: Could not get convergence for V=$Vtarget (MinVal=$MinVal)\n";
#	exit;
}
	my ($Vtotal, $Jtotal) = &CalIV($Vdiode, $J0, $n, $Jpv, $Rsh, $Rs);
	return ($Vdiode, $Jtotal);
}

sub CalIV
{
	my ($Vdiode, $J0, $n, $Jpv, $Rsh, $Rs) = @_;

	my $Jdiode = $J0 * (exp($ekT * $Vdiode / $n) - 1.0) - $Jpv;
	my $Jtotal = $Jdiode + $Vdiode / $Rsh;
	my $Vtotal = $Vdiode + $Rs * $Jtotal;
	return ($Vtotal, $Jtotal);
}

sub IVMinimizationFunc {
	my ($Vtarget, $J0, $n, $Jpv, $Rsh, $Rs, $pVars, $iPrintLevel) = @_;

	my $S = 0.0;
	my ($Vdiode) = @$pVars;
	my ($Vtotal, $Jtotal) = &CalIV($Vdiode, $J0, $n, $Jpv, $Rsh, $Rs);
	return ($Vtotal - $Vtarget)**2;
}

sub PDLDisplayFunc {
	my ($piddle) = @_;
#print "Display: piddle=$piddle\n";
}