#!/usr/bin/perl

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

use strict;

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

my $SourceCSV    = "c 1 Lsr.csv";
$SourceCSV = $ARGV[0] if($ARGV[0] ne '');
my $UseDark      = 1;
$UseDark = $ARGV[1] if($ARGV[1] ne '');
my $DoLSQ        = 1;
$DoLSQ = $ARGV[2] if($ARGV[2] ne '');

my $F            = 100.0; # mW/cm2
my $T            = 300.0; # K
my $nSkip        = 10;
my ($ILabel, $JLabel, $FitVMin, $FitVMax);
if($UseDark) {
	$ILabel  = "I.*";
#	$ILabel  = "Id.*";
	$JLabel  = "J.*";
	$FitVMin = -0.5;
	$FitVMax =  3.0;
}
else {
	$ILabel  = "I.*";
#	$ILabel  = "Ip.*";
	$JLabel  = "J.*";
	$FitVMin = -3.0;
	$FitVMax =  5.0;
	$ILabel       = ($UseDark)? "Id.*" : "Ip.*";
}
my $RsId         = 0; # Flag to choose whether total Rs is optimized or not
my $nSimplexIter = 1;

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 $ekT = $e / ($kB*$T);

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->mprint("Source   : $SourceCSV\n");
$txt->mprint("OutCSV   : $OutCSV\n");
$txt->mprint("FitCSV   : $FitCSV\n");
$txt->mprint("For dark?: $UseDark\n");
$txt->mprint("Current data taken from: $ILabel\n");
$txt->mprint("Light intensity: $F mW/cm2\n");
$txt->mprint("Temperature    : $T K\n");
$txt->mprint("Fitting range  : $FitVMin - $FitVMax V, every $nSkip points\n");
$txt->mprint("nSimplexIter   : $nSimplexIter\n");

print("Read [$SourceCSV]\n");

my $SC = new SolarCell;
$SC->SetLightIntensity($F);
$SC->SetTemplature($T);
my $iData = 0;
my ($prevV, $VIndex);
my $VLabel = "V.*";
my ($pV, $pI, $pJ, $nData) = $SC->ReadIVData($SourceCSV, $VLabel, $ILabel, $JLabel, 3, 
	sub {
		my ($csv, $pLabelList, $pDataArray, $pCur) = @_;
		if(!defined $VIndex) {
			$VIndex = $csv->FindLabelIndex([$VLabel], $pLabelList);
		}
		if(!defined $VIndex) {
			print "Error: Can not find V column data.\n";
			exit;
		}
		my $V = $pCur->[$VIndex];
#print "$iData for Vidx=$VIndex: $V\n";
		$prevV = $V - 1.0e-5 if(!defined $prevV);
		if($V <= $prevV) {
			return -1;
		}
		$prevV = $V;
		$iData++;
		return 1;
	}
	);
print("nData=$nData\n");
if(!$pJ) {
	print "Error: Can not find the data column I\n";
	exit;
}
if($nData == 0) {
	print "Error: nData = 0\n";
	exit;
}
my $prevV = $pV->[0] - 1.0e-5;
for(my $i = 0 ; $i < $nData ; $i++) {
	if($pV->[$i] <= $prevV) {
		$SC->{nData} = $i;
		last;
	}
	print "$i\t$pV->[$i]\t$pJ->[$i]\n";
	$prevV = $pV->[$i];
}

$SC->Analyze(1);
if(!$SC->WriteAnalyzeCSV($OutCSV)) {
	die "$!: Can not write to [$OutCSV]\n";
}

$txt->mprint("\n");
$txt->mprint("Source data\n");
$txt->mprint("  nData=$SC->{nData}\n");
$txt->mprint("  Electrode area: $SC->{Area} cm2\n");

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

my $Rs  = $SC->{RsOn};
my $Rsh = $SC->{RsRev};
my $Jpv = ($UseDark)? 0 : $SC->{Jpv};

my $J01 = $SC->{J0} / 2.0;
my $n1  = $SC->{n};
my $Rs1 = $Rs * 2.0;
my ($J02, $n2, $Rs2) = ($J01, $n1+1.0, $Rs1);
my $Rs  = 0.0;

$txt->mprint("\n");
$txt->mprint("Intial values:\n");
$txt->mprint("J01 = $J01 A/cm2\n");
$txt->mprint("n1  = $n1\n");
$txt->mprint("Rs1 = $Rs1 ohm cm2\n");
$txt->mprint("J02 = $J02 A/cm2\n");
$txt->mprint("n2  = $n2\n");
$txt->mprint("Rs2 = $Rs2 ohm cm2\n");
$txt->mprint("Rs  = $Rs ohm cm2\n");
$txt->mprint("Rsh = $Rsh ohm cm2\n");
$txt->mprint("Jpv = $Jpv A/cm2\n");
exit if(!$DoLSQ);

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

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

my $out = new JFile;
$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) = &CalSCFTotalIV($Vtarget, $J01, $n1, $Rs1, $J02, $n2, $Rs2, $Rs, $Rsh, $Jpv);
#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 ($J01, $n1, $Rs1, $J02, $n2, $Rs2, $Rs, $Rsh, $Jpv) = @$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) = &CalSCFTotalIV($Vtarget, $J01, $n1, $Rs1, $J02, $n2, $Rs2, $Rs, $Rsh, $Jpv);
		if($Jtotal != 0.0 and $pJ->[$i] != 0.0) {
			if($Jtotal * $pJ->[$i] < 0.0) {
				my $absd = abs(($Jtotal-$pJ->[$i]) / $pJ->[$i]);
				$S += (log($absd))**2;
			}
			else {
				my $absd = abs($Jtotal / $pJ->[$i]);
				$S += (log($absd))**2;
			}
		}
#print "$i: $Vtarget, $Jtotal, $pJ->[$i]\n";
		$c++;
	}
	$S += 1.0e5 if($J01 <= 0.0);
	$S += 1.0e5 if($n1  <= 0.0);
	$S += 1.0e5 if($Rs1 < 0.0);
	$S += 1.0e5 if($J02 <= 0.0);
	$S += 1.0e5 if($n2  <= 0.0);
	$S += 1.0e5 if($Rs2 < 0.0);
	$S += 1.0e5 if($Rs  < 0.0);
	$S += 1.0e5 if($Rsh < 0.0);
	$S += 1.0e5 if($Jpv < 0.0);
if($iPrintLevel >= 2) {
	printf("(J0,n,Rs)2,Rs,Rsh,Jpv\n");
	printf("  (%5.3g, %4.3f, %5.3g)\n", $J01, $n1, $Rs1);
	printf("  (%5.3g, %4.3f, %5.3g), %6.3f, %6.3f, %6.3g (%5.3g)\n", $J02, $n2, $Rs2, $Rs, $Rsh, $Jpv, $S);
}

	return $S;
}

sub CalSCFTotalIV
{
	my ($Vtarget, $J01, $n1, $Rs1, $J02, $n2, $Rs2, $Rs, $Rsh, $Jpv) = @_;

	my $dd = new DiodeRsParallelRsRshJpv;
	$dd->SetTemperature($T);
	return $dd->CalSCFTotalIV($Vtarget, $J01, $n1, $Rs1, $J02, $n2, $Rs2, $Rs, $Rsh, $Jpv);
}