#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';
#use lib 'd:\Working/lib';

BEGIN {
my $BaseDir = $ENV{'PerlDir'};
$BaseDir = $ENV{'TkPerlDir'} unless(defined $BaseDir);
@INC = ("$BaseDir/lib", "$BaseDir/VNL", @INC) if($BaseDir);
}

use strict;
#use warnings;

use Sci::Science;
use Sci::Optimize;
use Sci::OpticalMaterial;
use Sci::MultiLayer;
use Sci::Ellipsometry;

my $e  = Sci::e();
my $me = Sci::me();

my $SamplePath = "NT080.BEF";
my $SiPath     = "SI111025.SPE";
my $InitialSpectraCSV = "DielectricFunction-Initial.csv";
my $FinalSpectraCSV   = "DielectricFunction-Final.csv";

#=========================================================
# 大域変数
#=========================================================
my (@Tinitial, @Rinitial);
my $Angle = 70.0;

#=========================================================
# Object作製
# 解析条件設定
#=========================================================
my $optics = new Optics;

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

print "Optimization method: $Method\n";

#=========================================================
# 入力データ読み込み
#=========================================================
my $el = new Ellipsometry();
$el->SetIncidentAngle($Angle);

print "Read sample datafrom [$SamplePath].\n";
if(!$el->Read($SamplePath)) {
	print "Error: Can not read [$SamplePath].\n";
	exit;
}
my $pObsE    = $el->GetData('eV');
my $pObsEps1 = $el->GetData('e1');
my $pObsEps2 = $el->GetData('e2');
my $nData = @$pObsE;
my $Emin  = $pObsE->[0];
my $Emax  = $pObsE->[$nData-1];
my $Estep = $pObsE->[1] - $pObsE->[0];
my $nE    = $nData;
print "E range: $Emin - $Emax, $Estep ($nE)\n";
if(0) {
for(my $i = 0 ; $i < $nData ; $i++) {
	my $E  = $pObsE->[$i];
	my $e1 = $pObsEps1->[$i];
	my $e2 = $pObsEps2->[$i];
	print "$E: $e1, $e2\n";
}
}

#=========================================================
# 積層モデルの作製
#=========================================================
my $Layers = new MultiLayer;
my $Air    = new OpticalMaterial("air");
my $SiO2   = new OpticalMaterial("Specify", "SiO2");
my $aSi    = new OpticalMaterial("a-Si:H");
my $Si     = new OpticalMaterial("Specify", "Si");

my $FilmThickness      = 0.3;    # nm, Film thickness
my $SubstrateThickness = 0.5e-3; # m, Substrate thickness

$Layers->AddLayer($Air,  -1.0,                  0.0);
$Layers->AddLayer($SiO2, $FilmThickness*1.0e-9, 1.0);
$Layers->AddLayer($Si,   -1.0,                  0.0);

$Layers->SetIncidentAngle($Angle);

print "Layers:\n";
my $nLayer = $Layers->nLayer();
for(my $i = 0 ; $i < $nLayer ; $i++) {
	my $pLayer = $Layers->pLayer($i);
	my $name   = $pLayer->Name();
	my $d      = $Layers->LayerThickness($i) * 1.0e9; # nm
	print "$i: $name ($d nm)\n";
}

$SiO2->AddDielectricModel(
	"Constant",
	"Constant",
	"e1inf", 3.9,
	"e2inf", 0.0,
	);
$Si->AddDielectricModel(
	"Si",
	"EFile",
	"Path", $SiPath,
	);

#=========================================================
# 初期モデルをファイルに保存
#=========================================================
print "Save initial spectrum to [$InitialSpectraCSV].\n";
if(!$Layers->SaveDielectricFunctionsCSV($InitialSpectraCSV, $Emin, $Estep, $nE, $pObsEps1, $pObsEps2)) {
	print "Error: Can not write to [$InitialSpectraCSV].\n";
	exit;
}

#=========================================================
# フィッティング設定
#=========================================================
my ($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"FilmThickness", $Layers->pLayerThickness(1),    1, 1.0e-10,
			"SiO2-e1",       $Layers->pVariable(1, 'e1inf'), 1, 1,
			);

#=========================================================
# 最適化の実行
#=========================================================
my ($OptVars, $MinVal) = $optimize->Optimize(
				$Method, undef, undef, undef,
				$EPS, $nMaxIter, $iPrintLevel,
				sub { &CalS2($pObsE, $pObsEps1, $pObsEps2, $Layers, @_); },
				undef,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
print "\nOptimized at:\n";
$optimize->PrintParameters(1);

#=========================================================
# 元の変数に最適化した結果を入力
#=========================================================
my ($lFilmThickness, $SiO2e1) = @$OptVars;
$Layers->SetLayerThickness(1, $lFilmThickness);
$Layers->SetVariable(1, 'e1inf', $SiO2e1);

#=========================================================
# 最適化モデルをファイルに保存
#=========================================================
print "Save final spectrum to [$FinalSpectraCSV].\n";
if(!$Layers->SaveDielectricFunctionsCSV($FinalSpectraCSV, $Emin, $Estep, $nE, $pObsEps1, $pObsEps2)) {
	print "Error: Can not write to [$FinalSpectraCSV].\n";
	exit;
}

exit;

#==================================================================
# メインルーチン終了
#==================================================================


#==================================================================
# 関数
#==================================================================
sub CalS2
{
	my ($pObsE, $pObsEps1, $pObsEps2, $Layers, $pVars, $iPrintLevel) = @_;
	my ($lFilmThickness, $SiO2e1) = @$pVars;

	my $Penalty = 0.0;
	if($lFilmThickness < 1.0e-10 or $SiO2e1 < 1.0) {
		$Penalty = 1.0e10;
	}

	$Layers->SetLayerThickness(1, $lFilmThickness);
	$Layers->SetVariable(1, 'e1inf', $SiO2e1);

	my $nData = @$pObsE;
	my $S2 = 0.0;
	my $count = 0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $E  = $pObsE->[$i];
		my $ObsE1 = $pObsEps1->[$i];
		my $ObsE2 = $pObsEps2->[$i];

		my ($Psi, $Delta) = $Layers->CalPsiDelta($E);
		my ($e1c, $e2c)   = $el->PsiDeltaToEps($Psi, $Delta, $Angle);
		
		$S2 += ($ObsE1 - $e1c)**2 + ($ObsE2 - $e2c)**2;
		$count++;
	}

	$S2 = sqrt($S2 / $count) + $Penalty;

print "S2=$S2: Thickness: $lFilmThickness\n";
print "        SiO2-e1: $SiO2e1\n";

	return $S2;
}