#!/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::Material;
use Sci::Optimize;
use CSV;

my $SourceCSV = "VBM-DOS.csv";
my $FitCSV    = "VBM-DOS-Fit.csv";

my $delta         = 5.5; # nm
my $FilmThickness = 30.0; # nm

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

my $EPS = 1.0e-7;
my $nMaxIter = 1000;
my $iPrintLevel = 2;

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 $pQ = $csv->GetXData("ANGLE.*");
#my $pI = $csv->GetYData("N\\(as\\).*");
my $pI = $csv->GetYData("N\\(dry\\).*");
#my $pI = $csv->GetYData("N\\(wet\\).*");
my $nData = @$pQ;
print "nData = $nData\n";
print "pQ=$pQ\n";
print "pI=$pI\n";
print "Q=", join(', ', @$pQ), "\n";
print "I=", join(', ', @$pI), "\n";

my $loga1 = log(9e20);
my $l1    = 0.5; # nm
my $loga2 = log(1.0);
my $l2    = 2.5; # nm
my $logbg = log(1.0);

my @Guess = ( $loga1, $l1, $loga2, $l2, $logbg);
my @OptId = (      1,   1,      1,   1,      1);
my @Scale = (    1.0, 1.0,    1.0, 1.0,    1.0);

print "use $Method\n";

my $Opt = new Optimize;

#必要なければ、PDLDisplayFunc以降は指定する必要はない
my ($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				\&MinimizationFunc, \&PDLDisplayFunc,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
( $loga1, $l1, $loga2, $l2, $logbg) = @$OptVars;
print "Optimized at (",join(',',@{$OptVars}),")=$MinVal\n";

my $out = new JFile;
$out->Open($FitCSV, "w") or die "$!: Can not write to [$FitCSV]\n";
$out->print("Q(degree),deff(nm),Iobs,Ical\n");
for(my $i = 0 ; $i < $nData ; $i++) {
	my $Q = $pQ->[$i];
	my $sinQ = sin($Q*$torad);
	my $deff = $delta * $sinQ;
	
	my $Ical = &Cal($Q, $loga1, $l1, $loga2, $l2, $logbg);
	$out->print("$Q,$deff,$pI->[$i],$Ical\n");
}

my $a1 = exp($loga1) * $delta;
my $a2 = exp($loga2) * $delta;
my $bg = exp($logbg);
$out->print("loga1=,$loga1,,DOS1=,$a1,cm-3\n");
$out->print("l1=,$l1,nm\n");
$out->print("loga2=,$loga2,,DOS2=,$a2,cm-3\n");
$out->print("l2=,$l2,nm\n");
$out->print("logbg=,$logbg,,bg=,$bg,cm-3\n");
$out->print("\n");

my $n = 100;
my $dstep = $FilmThickness / $n;
$out->print("d(nm),DOS(cm-3)\n");
for(my $i = 0 ; $i <= $n ; $i++) {
	my $d = $dstep * $i;
	my $DOS = $a1 * exp(-$d / $l1) + $a2 * exp(-$d / $l2) + $bg;
	$out->print("$d,$DOS\n");
}
$out->Close();

exit;

#==============================================
# 最小化関数のサブルーチン
# MinimizationFunc    : 
# PDLDisplaySimplex: PDL::Opt::Simplex用のDisplay関数
#==============================================
sub MinimizationFunc {
	my ($pVars, $iPrintLevel) = @_;

	my $S = 0.0;
	my ($loga1, $l1, $loga2, $l2, $logb) = @$pVars;
	for(my $i = 0 ; $i < $nData ; $i++) {
		next if(!defined $pI->[$i] or $pI->[$i] eq '');
		my $Q = $pQ->[$i];
		my $Ical = &Cal($Q, $loga1, $l1, $loga2, $l2, $logbg);
		$S += log($Ical / $pI->[$i])**2;
	}
	$S += 1.0e5 if($l1 < 0.0);
	$S += 1.0e5 if($l2 < 0.0);
printf("loga1,l1,loga2,l2,logbg=%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t(%g)\n", 
		$loga1, $l1, $loga2, $l2, $logbg, $S) if($iPrintLevel >= 2);
	return $S;
}

sub Cal
{
	my ($Q, $loga1, $l1, $loga2, $l2, $logbg) = @_;
	
	my $sinQ = sin($Q*$torad);
	my $a1 = exp($loga1);
	my $a2 = exp($loga2);
	my $bg = exp($logbg);
	my $invleff1 = $sinQ / $l1 + 1.0 / $delta;
	my $invleff2 = $sinQ / $l2 + 1.0 / $delta;

	return $bg 
			+ $a1 / $invleff1 * (1.0 - exp(-$invleff1 * $FilmThickness / $sinQ))
			+ $a2 / $invleff2 * (1.0 - exp(-$invleff2 * $FilmThickness / $sinQ));

}

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