#!/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);
}

use strict;
#use warnings;
#use Math::Complex;

use CSV;
use Sci::Science;
#use Sci::Optimize;
use Optimize;
#use Sci::Optics;
use Optics;
use GraphData;
use OpticalMaterial;
use MultiLayer;

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

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

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

my $TRDataPath = "glassT.csv";
my $InitialCSV = "GlassInitial.csv";
my $FinalCSV   = "GlassFinal.csv";

my $Model = "TaucLorentz";
#my $Model = "Caucy";

my $nSkip = 10;
my $WLMin = 267e-9; # nm
my $WLMax = 2400e-9; # nm
if($Model eq 'Caucy') {
	$WLMin = 267e-9; #390e-9; # nm
	$WLMax = 2000e-9; # nm
}

#=========================================================
# 入力データ読み込み
#=========================================================
print "Read transmission spectrum from [$TRDataPath].\n";
my ($nData, $pTRLabelArray, $pWL, $pT) = CSV::GetArraysFromFile($TRDataPath);
if(!defined $nData) {
	print "Can not read[$TRDataPath].\n";
	exit;
}

print "Optimization method: $Method\n";
print "Optical model      : $Model\n";
print "Fitting WL range   : ", $WLMin*1e9, " - ", $WLMax*1e9, " nm\n";
print "nData              : $nData\n";

#=========================================================
# 積層モデルの作製
#=========================================================
my $thickness  = 0.5e-3; # m, Thickness
my $Layers = new MultiLayer;
my $Air       = new OpticalMaterial("air");
my $Substrate = new OpticalMaterial();
$Layers->AddLayer($Air,       -1.0,       0.0);
$Layers->AddLayer($Substrate, $thickness, 0.0);
$Layers->AddLayer($Air,       -1.0,       0.0);

# For Tauc-Lorentz
my $e1inf = 0.6385;
my $e2inf = 0.0;
my $A     = 252.1;
my $Eg    =  8.5; # eV, Tauc gap
my $En0   = 12.0; # eV, Lorentz oscillation energy
my $C     =  8.542; # eV, Damping factor

# For Urbach
my @Em      = (5.0,       5.0);
my @E0      = (0.3668,    0.8); # eV, Urbach energy
my @AUrbach = (0.0005036, 0.0005);

# For Caucy
my @CaucyA = (1.152302, 3364.8616, -1541.0569, 0, 0, 0);
my @Caucyn = (0,        2,        4          , 6, 8, 10);

SetDielectricParameters(
		$Substrate, 
		$e1inf, $e2inf, $A, $Eg, $En0, $C,
		$Em[0], $E0[0], $AUrbach[0],
		$Em[1], $E0[1], $AUrbach[1],
		\@CaucyA, \@Caucyn,
		);

#=========================================================
# フィッティング設定
#=========================================================
my ($pName, $ppVariable, $pGuess, $pOptId, $pScale);
if($Model eq 'TaucLorentz') {
	($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"e1inf",    \$e1inf,      1, $e1inf,
			"e2inf",    \$e2inf,      0, 1,
			"A",        \$A,          1, $A,
			"Eg",       \$Eg,         0, $Eg,
			"En0",      \$En0,        0, $En0,
			"C",        \$C,          1, $C,
			"Em0",      \$Em[0],      0, $Em[0],
			"E00",      \$E0[0],      0, $E0[0],
			"AUrbach0", \$AUrbach[0], 1, $AUrbach[0],
			"Em1",      \$Em[1],      0, $Em[1],
			"E01",      \$E0[1],      1, $E0[1],
			"AUrbach1", \$AUrbach[1], 1, $AUrbach[1],
		);
}
elsif($Model eq 'Caucy') {
	($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"CaucyA0",  \$CaucyA[0],    1, $CaucyA[0],
			"CaucyA1",  \$CaucyA[1],    1, $CaucyA[1],
			"CaucyA2",  \$CaucyA[2],    1, $CaucyA[2],
			"CaucyA3",  \$CaucyA[3],    1, $CaucyA[3],
			"CaucyA4",  \$CaucyA[4],    1, $CaucyA[4],
			"CaucyA5",  \$CaucyA[5],    1, $CaucyA[5],
		);
}

#=========================================================
# 初期モデルをファイルに保存
#=========================================================
print "Save initial spectrum to [$InitialCSV].\n";
&SaveSpectrum($InitialCSV);

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

#=========================================================
# 元の変数に最適化した結果を入力
#=========================================================
$optimize->UpdateParameters();
SetDielectricParameters(
	$Substrate, 
	$e1inf, $e2inf, $A, $Eg, $En0, $C,
	$Em[0], $E0[0], $AUrbach[0],
	$Em[1], $E0[1], $AUrbach[1],
	@CaucyA, \@Caucyn,
	);

#=========================================================
# 最適化モデルをファイルに保存
#=========================================================
print "Save optimized result to [$FinalCSV].\n";
&SaveSpectrum($FinalCSV, \@Tinitial, \@Rinitial);

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


#==================================================================
# 関数
#==================================================================
# AddDielectricModelは、初回はモデルを追加する。
# $nameと同じものが存在したら、パラメータの更新を行う
sub SetDielectricParameters
{
	my ($Substrate,
		$e1inf, $e2inf, $A, $Eg, $En0, $C,
		$Em0, $E00, $AUrbach0,
		$Em1, $E01, $AUrbach1,
		$pCaucyA, $pCaucyn,
		) = @_;

#print "e=($e1inf, $e2inf), A=$A, Eg=$Eg, En0=$En0, C=$C\n";
#print "Em=$Em0, E=$E00, A=$AUrbach0\n";
#print "Em=$Em1, E=$E01, A=$AUrbach1\n";

	if($Model eq 'TaucLorentz') {
		$Substrate->AddDielectricModel(
			"TaucLorentz1",
			"TaucLorentz",
			"e1inf", $e1inf,
			"e2inf", $e2inf,
			"A",     $A,
			"Eg",    $Eg,
			"En0",   $En0,
			"C",     $C
			);
		$Substrate->AddDielectricModel(
			"UrbachTail1",
			"Urbach",
			"Em", $Em0,
			"E0", $E00,
			"A",  $AUrbach0,
			);
		$Substrate->AddDielectricModel(
			"UrbachTail2",
			"Urbach",
			"Em", $Em1,
			"E0", $E01,
			"A",  $AUrbach1,
			);
	}
	elsif($Model eq 'Caucy') {
		$Substrate->AddDielectricModel(
			"Caucy1",
			"Caucy",
			"pA", $pCaucyA,
			"pn", $pCaucyn,
			);
	}
}

sub CalS2
{
	my ($pVars, $iPrintLevel) = @_;
	my ($le1inf, $le2inf, $lA, $lEg, $lEn0, $lC, 
		$lEm0, $lE00, $lAUrbach0,
		$lEm1, $lE01, $lAUrbach1,
		);
	my @lCaucyA;

	if($Model eq 'TaucLorentz') {
		($le1inf, $le2inf, $lA, $lEg, $lEn0, $lC, 
		 $lEm0, $lE00, $lAUrbach0,
		 $lEm1, $lE01, $lAUrbach1,
		 ) = @$pVars;
	}
	elsif($Model eq 'Caucy') {
		@lCaucyA = @$pVars;
	}
	SetDielectricParameters(
		$Substrate, 
		$le1inf, $le2inf, $lA, $lEg, $lEn0, $lC,
		$lEm0, $lE00, $lAUrbach0,
		$lEm1, $lE01, $lAUrbach1,
		\@lCaucyA, \@Caucyn,
		);

	my $S2 = 0.0;
	my $count = 0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		next if($i % $nSkip != 0);
		my $WL = $pWL->[$i]; # nm
		next if($WL*1e-9 < $WLMin or $WL*1e-9 > $WLMax);
#		my $nu = Optics::nmToHz($WL); # Hz
		my $E  = Optics::nmToeV($WL);

		my $T  = $pT->[$i];
		my $Tnet = $T;
		my ($e1cal, $e2cal, $Rcal, $Tcal) = CalEpsTR($WL, $E);
		$S2 += ($Tnet - $Tcal)**2;
		$count++;
	}
	$S2 = sqrt($S2 / $count);
	PrintParameters($S2, @$pVars);
	return $S2;
}

sub CalEpsTR
{
	my ($WL, $E) = @_;

	my ($e1air1, $e2air1, $e1sub, $e2sub, $e1air2, $e1air3, $Tcal, $Rcal) = $Layers->CalEpsTRForSubstrateAtNormalncidence($E);
	return  ($e1sub, $e2sub, $Rcal, $Tcal);
}

sub SaveSpectrum
{
	my ($path, $pTinitial, $pRinitial) = @_;

	open(OUT, ">$path") or die "$!\n";
	if($pTinitial) {
		print OUT "WL(nm),nu(Hz),E(eV),T(obs),T(initial),T(final),R(initial),R(final),alpha(cm-1),e1(cal),e2(cal),n(cal),k(cal)\n";
	}
	else {
		print OUT "WL(nm),nu(Hz),E(eV),T(obs),T(initial),R(initial),alpha(cm-1),e1(cal),e2(cal),n(cal),k(cal)\n";
	}
	for(my $i = 0 ; $i < $nData ; $i++) {
#		next if(($i+1) % $nSkip == 0);
		my $WL = $pWL->[$i]; # nm
		next if($WL*1e-9 < $WLMin or $WL*1e-9 > $WLMax);
		my $nu = Optics::nmToHz($WL); # Hz
		my $E  = Optics::nmToeV($WL);

		my $T  = $pT->[$i];
		my $Tnet = $T;

		my ($e1cal, $e2cal, $Rcal, $Tcal) = CalEpsTR($WL, $E);
		my $alpha = -log($Tnet) / ($thickness * 1.0e2); # cm-1, alpha
		my ($ncal, $kcal) = Optics::EpsToNK($e1cal, $e2cal);
		if($pTinitial) {
			print OUT "$WL,$nu,$E,$T,$pTinitial->[$i],$Tcal,$pRinitial->[$i],$Rcal,$alpha,$e1cal,$e2cal,$ncal,$kcal\n";
		}
		else {
			$Tinitial[$i] = $Tcal if(!$pTinitial);
			$Rinitial[$i] = $Rcal if(!$pTinitial);
			print OUT "$WL,$nu,$E,$T,$Tcal,$Rcal,$alpha,$e1cal,$e2cal,$ncal,$kcal\n";
		}
	}
	close(OUT);
	return 1;
}

sub PrintParameters
{
	my ($val, @Vars) = @_;
	printf("f=%12.6g: ", $val);
	for(my $i = 0 ; $i < @Vars ; $i++) {
		if($i == 0) {
			printf("%6.4g", $Vars[$i]);
		}
		else {
			printf(", %6.4g", $Vars[$i]);
		}
	}
	print("\n");

}