#!/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 Math::Complex;

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

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

#=========================================================
# 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;

my $SubstrateEpsPath = "Glass-Substrate-DielectricFunction.csv";
my $TRDataPath = "Zn035-AsDeposited.csv";
my $InitialCSV = "FilmInitial.csv";
my $FinalCSV   = "FilmFinal.csv";

my $nSkip = 10;
my $WLMin = 310e-9; # nm
my $WLMax = 2400e-9; # nm
$WLMin =  600e-9; # nm
$WLMax = 2500e-9; # nm

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

print "Optimization method: $Method\n";
print "Measured WL range  : $WL0 - $WL1, $WLStep ($nData)\n";
print "Fitting WL range   : ", $WLMin*1e9, " - ", $WLMax*1e9, " nm\n";

#=========================================================
# 積層モデルの作製
#=========================================================
my $SubstrateThickness = 0.5e-3; # m, Substrate thickness

#my $FittingMode = "TR";
#my $FittingMode = "T";
#my $FittingMode = "R";
my $FittingMode = "alpha";

my $Ne = 1.0e20 * 1.0e6;      # m3,    carrier density
my $sigma = 145.0e0 * 1.0e-2; # S/m,   conductivity
my $mobility = 9.5 * 1.0e-4;  # m2/Vs, mobility

my $ScalingFactor      = 1.0161;
my $FilmThickness      = 365.878; # nm, Film thickness

# Constant background
my $e1inf =   0.812686;
my $e2inf =   0.0;

# For Tauc-Lorentz
my $A     =  90.1;
my $Eg    =   2.786; # eV, Tauc gap
my $En0   =   6.290; # eV, Lorentz oscillation energy
my $C     =   6.786; # eV, Damping factor

# For Urbach
my @Em      = (3.1,      3.1);
my @E0      = (1.905,    0.5); # eV, Urbach energy
my @AUrbach = (0.0,      0.0000);
#my @AUrbach = (0.1034,   0.0000);

# For Lorentz
my $LorentzNe    = 0.4;    # 1e26 m-3
my $LorentzE0    = 2.805;  # eV
my $Lorentzgamma = 0.3091; # eV

# For Drude
my $DrudeEp      = 0.764814; # eV
my $Drudetau     = 9.81047;  # fs

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

my $Layers = new MultiLayer;
my $Air       = new OpticalMaterial("air");
my $Substrate = new OpticalMaterial();
my $Film      = new OpticalMaterial();
$Layers->AddLayer($Air,       -1.0,                  0.0);
$Layers->AddLayer($Film,      $FilmThickness*1.0e-9, 1.0);
$Layers->AddLayer($Substrate, $SubstrateThickness,   0.0);
$Layers->AddLayer($Air,       -1.0,                  0.0);

$Substrate->AddDielectricModel(
	"${SubstrateEpsPath}-1",
	"File",
	"Path", $SubstrateEpsPath,
	"pWL",  $pWL,
	);

SetDielectricParameters(
		$Film,  $FilmThickness,
		$e1inf, $e2inf, $A, $Eg, $En0, $C,
		$Em[0], $E0[0], $AUrbach[0],
		$Em[1], $E0[1], $AUrbach[1],
		$LorentzNe, $LorentzE0, $Lorentzgamma,
		$DrudeEp, $Drudetau,
		\@CaucyA, \@Caucyn,
		);

print "File e\n";
for(my $E = 0.4 ; $E < 4.0 ; $E += 0.2) {
	my ($e1, $e2) = $Film->CalEps($E);
	printf("%6.4g eV: e=(%8.4g, %8.4g)\n", $E, $e1, $e2);
}
#exit;

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

#=========================================================
# フィッティング設定
#=========================================================
my ($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"ScalingFactor", \$ScalingFactor, 1, $ScalingFactor * 0.1,
			"FilmThickenss", \$FilmThickness, 1, $FilmThickness * 0.1,
			);
($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"e1inf",         \$e1inf,         1, 1,
			"e2inf",         \$e2inf,         0, 1,
			);
($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"A",             \$A,             0, $A,
			"Eg",            \$Eg,            0, $Eg  * 0.1,
			"En0",           \$En0,           0, $En0 * 0.1,
			"C",             \$C,             0, $C,
			);
($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"Em0",           \$Em[0],         0, $Em[0] * 0.1,
			"E00",           \$E0[0],         0, $E0[0] * 0.3,
			"AUrbach0",      \$AUrbach[0],    0, $AUrbach[0],
			);
($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"Em1",           \$Em[1],         0, $Em[1] * 0.1,
			"E01",           \$E0[1],         0, $E0[1] * 0.1,
			"AUrbach1",      \$AUrbach[1],    0, $AUrbach[1],
			);
($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"LorentzNe1",    \$LorentzNe,     0, $LorentzNe * 0.5,
			"LorentzE01",    \$LorentzE0,     0, $LorentzE0 * 0.1,
			"LorentzG1",     \$Lorentzgamma,  0, $Lorentzgamma * 0.3,
			);
($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"DrudeEp1",      \$DrudeEp,       0, $DrudeEp * 0.5,
			"Drudetau1",     \$Drudetau,      0, $Drudetau * 0.5,
			);
($pName, $ppVariable, $pGuess, $pOptId, $pScale) = $optimize->AddParameters(
			"CaucyA0",       \$CaucyA[0],     0, $CaucyA[0],
			"CaucyA1",       \$CaucyA[1],     0, $CaucyA[1],
			"CaucyA2",       \$CaucyA[2],     0, $CaucyA[2],
			"CaucyA3",       \$CaucyA[3],     0, $CaucyA[3],
			"CaucyA4",       \$CaucyA[4],     0, $CaucyA[4],
			"CaucyA5",       \$CaucyA[5],     0, $CaucyA[5],
		);

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

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

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

my $e  = Sci::e();
my $me = Sci::me();
my $meeff           = Optics::CalEffectiveMassFromPlasmaEnergy($DrudeEp, $Ne);
my $mobilityFCA     = $e * $Drudetau * 1.0e-15 / $meeff / $me * 1.0e4; # cm2/Vs
my $conductivityFCA = $e * $Ne * 1.0e-6 * $mobilityFCA; # S/cm
print "Ne(obs)     : ", $Ne*1.0e-6, " cm-3\n";
print "me*         : $meeff\n";
print "mobility    : (FCA)$mobilityFCA  (Meas)", $mobility*1.0e4, " cm2/Vs\n";
print "conducitivty: (FCA)$conductivityFCA  (Meas)", $sigma*1.0e2, " S/cm\n";

exit;

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


#==================================================================
# 関数
#==================================================================
# AddDielectricModelは、初回はモデルを追加する。
# $nameと同じものが存在したら、パラメータの更新を行う
sub SetDielectricParameters
{
	my ($Material,$FilmThickness,
		$e1inf, $e2inf, $A, $Eg, $En0, $C,
		$Em0, $E00, $AUrbach0,
		$Em1, $E01, $AUrbach1,
		$LorentzNe, $LorentzE0, $Lorentzgamma,
		$DrudeEp, $Drudetau,
		$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";

	$FilmThickness = 1.0e-10 if($FilmThickness < 1.0e-10);
	$e1inf = 1.0e-15 if($e1inf < 1.0e-15);
	$e2inf = 1.0e-15 if($e2inf < 1.0e-15);
	$Eg    = 1.0e-15 if($Eg    < 1.0e-15);
	$En0   = 1.0e-15 if($En0   < 1.0e-15);
	$A     = 1.0e-15 if($A     < 1.0e-15);
	$C     = 1.0e-15 if($C     < 1.0e-15);

	$Material->SetThickness($FilmThickness*1.0e-9);
	$Material->AddDielectricModel(
		"Constant1",
		"Constant",
		"e1inf", $e1inf,
		"e2inf", $e2inf,
		);

	$Material->AddDielectricModel(
		"TaucLorentz1",
		"TaucLorentz",
		"A",     $A,
		"Eg",    $Eg,
		"En0",   $En0,
		"C",     $C
		);
	$Material->AddDielectricModel(
		"UrbachTail1",
		"Urbach",
		"Em", $Em0,
		"E0", $E00,
		"A",  $AUrbach0,
		);
	$Material->AddDielectricModel(
		"UrbachTail2",
		"Urbach",
		"Em", $Em1,
		"E0", $E01,
		"A",  $AUrbach1,
		);
	$Material->AddDielectricModel(
		"Lorentz1",
		"Lorentz",
		"Ne", $LorentzNe*1.0e26, # m-3
		"E0", $LorentzE0,
		"G",  $Lorentzgamma,
		);
	$Material->AddDielectricModel(
		"Drude1",
		"Drude",
		"Ep",  $DrudeEp,
		"tau", $Drudetau*1.0e-15, # s
		);
	$Material->AddDielectricModel(
		"Caucy1",
		"Caucy",
		"pA", $pCaucyA,
		"pn", $pCaucyn,
		);
}

sub CalS2
{
	my ($pVars, $iPrintLevel) = @_;
	my ($lScalingFactor, $lFilmThickness,
		$le1inf, $le2inf, $lA, $lEg, $lEn0, $lC, 
		$lEm0, $lE00, $lAUrbach0,
		$lEm1, $lE01, $lAUrbach1,
		$lLorentzNe, $lLorentzE0, $lLorentzgamma,
		$lDrudeEp, $lDrudetau,
		$lCaucyA0, $lCaucyA1, $lCaucyA2,
		$lCaucyA3, $lCaucyA4, $lCaucyA5,
		) = @$pVars;

	my @lCaucyA = ($lCaucyA0, $lCaucyA1, $lCaucyA2, $lCaucyA3, $lCaucyA4, $lCaucyA5);
	my $Penalty = 0.0;
	if($lFilmThickness < 1.0e-10 or $lEg < 1.0e-15 or $lEn0 < 1.0e-15 or 
	   $lA < 1.0e-15 or $lC < 1.0e-15 or $lLorentzNe <= 0.0 or
	   $lDrudeEp < 1.0e-15 or $lDrudetau < 1.0e-20) {
		$Penalty = 10.0;
	}
	SetDielectricParameters(
		$Film, $lFilmThickness,
		$le1inf, $le2inf, $lA, $lEg, $lEn0, $lC,
		$lEm0, $lE00, $lAUrbach0,
		$lEm1, $lE01, $lAUrbach1,
		$lLorentzNe, $lLorentzE0, $lLorentzgamma,
		$lDrudeEp,   $lDrudetau,
		\@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] * $lScalingFactor;
		my $R  = $pR->[$i] * $lScalingFactor;
		my ($e1cal, $e2cal, $Rcal, $Tcal) = CalEpsTR($E);
		my $alphaobs = -log($T   /(1.0-$R   )) / ($lFilmThickness*1.0e-7);
		my $alphacal = -log($Tcal/(1.0-$Rcal)) / ($lFilmThickness*1.0e-7);
#print "E=$E  T=$Tcal  R=$Rcal  a=$alphacal d=$lFilmThickness\n";

		if($FittingMode =~ /T/i) {
			$S2 += ($T - $Tcal)**2;
		}
		if($FittingMode =~ /R/i) {
			$S2 += ($R - $Rcal)**2;
		}
		if($FittingMode =~ /alpha/i) {
			$S2 += ($alphaobs - $alphacal)**2;
		}
		$count++;
	}

	$S2 = sqrt($S2 / $count) + $Penalty;
	PrintParameters($S2, 0, @$pVars);
	return $S2;
}

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

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

sub SaveSpectrum
{
	my ($path, $thickness, $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(obs),R(initial),R(final),"
				 ."aa(obs),aa(cal),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(obs),R(initial),"
				 ."aa(obs),aa(cal),alpha(cm-1),e1(cal),e2(cal),n(cal),k(cal)\n";
	}
	for(my $i = 0 ; $i < $nData ; $i++) {
		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] * $ScalingFactor;
		$T = 1.0e-12 if($T <= 0.0);
		my $R  = $pR->[$i] * $ScalingFactor;

		my ($e1cal, $e2cal, $Rcal, $Tcal) = CalEpsTR($E);
		my $alpha = -log($T) / ($thickness * 1.0e2); # cm-1, alpha
		my ($ncal, $kcal) = Optics::EpsToNK($e1cal, $e2cal);
		my $alphaobs = -log($T/(1.0-$R)) / ($thickness * 1.0e2);
		my $alphacal = -log($Tcal/(1.0-$Rcal)) / ($thickness * 1.0e2);

		if($pTinitial) {
			print OUT "$WL,$nu,$E,$T,$pTinitial->[$i],$Tcal,"
					 ."$R,$pRinitial->[$i],$Rcal,$alphaobs,$alphacal,$alpha,"
					 ."$e1cal,$e2cal,$ncal,$kcal\n";
		}
		else {
			$Tinitial[$i] = $Tcal if(!$pTinitial);
			$Rinitial[$i] = $Rcal if(!$pTinitial);
			print OUT "$WL,$nu,$E,$T,$Tcal,$R,$Rcal,$alphaobs,$alphacal,$alpha,"
					 ."$e1cal,$e2cal,$ncal,$kcal\n";
		}
	}
	close(OUT);
	return 1;
}

sub PrintParameters
{
	my ($val, $IsDetail, @Vars) = @_;

	if($IsDetail) {
		print("  Return value: $val\n");
		for(my $i = 0 ; $i < @Vars ; $i++) {
			printf("    %2d: %15s: %g\n", $i, $pName->[$i], $Vars[$i]);
		}
		return;
	}

	printf("f=%5.3g: ", $val);
	for(my $i = 0 ; $i < @Vars ; $i++) {
		my $v = $Vars[$i];
		if($i == 0) {
			printf("%6.4g", $v);
		}
		else {
			printf(",%6.4g", $v);
		}
		if(($i+1) % 9 == 0) {
			print("\n         ");
		}
	}
	print("\n");

}