#===============================================
# OpticalMaterial
#===============================================
package OpticalMaterial;
use Common;
@ISA = qw(Common);

#公開したいサブルーチン
#@EXPORT = qw(aa );

use strict;
use Math::Complex;

use Utils;
use JFile;
use CSV;
use GraphData;
use Sci::GeneralFileFormat;
use Sci::Science;
use Sci::Optics;
use Sci::Ellipsometry;
use Sci::Optimize;

my $pi   = Sci::pi();
my $e0   = Sci::e0();
my $e    = Sci::e();
my $c    = Sci::c();
my $h    = Sci::h();
my $hbar = Sci::hbar();
my $me   = Sci::me();
my $kB   = Sci::kB();

sub nDielectricModel { return shift->{nDielectricModel}; }
sub pDielectricModel { return shift->{pDielectricModel}; }
sub Thickness        { return shift->{Thicnkess}; }
sub SetThickness     { my($this,$d)=@_; return $this->{Thicnkess} = $d; }
sub SetCoherency     { my($this,$c)=@_; return $this->{Coherency} = $c; }
sub Coherency        { my($this)=@_; return 1.0 if(!defined $this->{Coherency}); return $this->{Coherency}; }
sub SetName          { my($this,$n)=@_; return $this->{Name} = $n; }
sub Name             { my($this)=@_; return $this->{Name}; }
sub SetVariable
{
	my ($this, $Name, $value) = @_;
	my ($ModelName, $VarName) = ($Name =~ /^(.*?)::(.*)$/);
	if(defined $ModelName) {
		my $n = $this->nDielectricModel();
		my $p = $this->pDielectricModel();
		if($ModelName =~ /^\d+$/) {
			return $p->[$ModelName]->{$VarName} = $value;
		}
		for(my $i = 0 ; $i < $n ; $i++) {
			if($p->[$i]->{name} eq $ModelName) {
				return $p->[$i]->{$VarName} = $value;
			}
		}
		return undef;
	}
#print "Name=$Name\n";
	return $this->{$Name} = $value;
}

sub pVariable
{
	my ($this, $Name) = @_;
	my ($ModelName, $VarName) = ($Name =~ /^(.*?)::(.*)$/);
	if(defined $ModelName) {
		my $n = $this->nDielectricModel();
		my $p = $this->pDielectricModel();
		if($ModelName =~ /^\d+$/) {
			return \$p->[$ModelName]->{$VarName};
		}
		for(my $i = 0 ; $i < $n ; $i++) {
			if($p->[$i]->{name} eq $ModelName) {
#print "Hit: $ModelName=$VarName => ", $p->[$i]->{$VarName}, "\n";
				return \$p->[$i]->{$VarName};
			}
		}
		return undef;
	}
#print "Name=$Name\n";
	return \$this->{$Name};
}

sub new
{
	my ($module, $material, $name) = @_;
	my $this = {};
	bless $this;
	$this->Initialize();
	$this->MakeMaterial($material, $name);
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub Initialize
{
	my ($this) = @_;
	$this->{nDielectricModel}   = 0;
	$this->{pDielectricModel}   = [];
	$this->{idxDielectricModel} = {};
}

sub MakeMaterial
{
	my ($this, $material, $name) = @_;
	$this->SetName($material);
	if($material =~ /^specify$/i) {
		$this->SetName($name);
	}
	elsif($material =~ /^air$/i) {
		$this->AddDielectricModel(
			"Air",
			"Constant",
			"e1inf", 1.0,
			"e2inf", 0.0,
			);
	}
	elsif($material =~ /^a-Si:H(-\d+)?$/i) {
		my $idx = $1;
		if(!defined $idx or $idx eq '-1') {
#APL 69 (1996) 371; 2137(ERRATUM)
			$this->AddDielectricModel(
				"TaucLorentz",
				"TaucLorentz",
				"e1inf", 1.15,
				"e2inf", 0.0,
				"Eg",    1.20,
				"En0",   3.45,
				"A",     122,
				"C",     2.54
				);
		}
		elsif($idx eq '-2') {
#APL 69 (1996) 371; 2137(ERRATUM)
			$this->AddDielectricModel(
				"TaucLorentz",
				"TaucLorentz",
				"e1inf", 0.74,
				"e2inf", 0.0,
				"Eg",    1.39,
				"En0",   3.58,
				"A",     169,
				"C",     2.15
				);
		}
	}
	elsif($material =~ /^SiO$/i) {
#APL 69 (1996) 371; 2137(ERRATUM)
		$this->AddDielectricModel(
			"TaucLorentz",
			"TaucLorentz",
			"e1inf", 2.06,
			"e2inf", 0.0,
			"Eg",    1.59,
			"En0",   3.93,
			"A",     58.0,
			"C",     4.23
			);
	}
	elsif($material =~ /^As2S3$/i) {
#APL 69 (1996) 371; 2137(ERRATUM)
		$this->AddDielectricModel(
			"TaucLorentz",
			"TaucLorentz",
			"e1inf", 2.50,
			"e2inf", 0.0,
			"Eg",    2.37,
			"En0",   3.75,
			"A",    161,
			"C",     4.60
			);
	}
	elsif($material =~ /^Si3N4$/i) {
#APL 69 (1996) 371; 2137(ERRATUM)
		$this->AddDielectricModel(
			"TaucLorentz",
			"TaucLorentz",
			"e1inf", 3.10,
			"e2inf", 0.0,
			"Eg",    4.50,
			"En0",   6.78,
			"A",     59.2,
			"C",     0.49
			);
	}
	elsif($material =~ /^DrudeDefault$/i) {
		$this->AddDielectricModel(
			"Drude",
			"Drude",
			"e1inf", 1.0,
			"e2inf", 0.0,
			"Ep",    0.764814, # eV
			"tau",   9.8e-15,  # fs
			);
	}
	elsif($material =~ /^LorentzDefault$/i) {
		$this->AddDielectricModel(
			"Lorentz",
			"Lorentz",
			"e1inf", 1.0,
			"e2inf", 0.0,
			"Ne",    0.4e21 * 1e6, # m-3
			"E0",    2.805,      # eV
			"gamma", 0.3091,     # eV
			"fa",    1.0,
			);
	}
}

sub SaveDielectricFunctionsCSV
{
	my ($this, $path, $Emin, $EStep, $nE, $pObs1, $pObs2) = @_;

	my $out = new JFile();
	if(!$out->Open($path, "w")) {
		return undef;
	}
	$out->print("E(eV)");
	if($pObs1) {
		$out->print(",e1(obs),e2(obs)");
	}
	$out->print(",e1(cal),e2(cal)\n");

	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Emin + $i * $EStep;
		my ($e1c, $e2c)   = $this->CalEps($E);
		$out->print("$E");
		if($pObs1) {
			$out->print(",$pObs1->[$i],$pObs2->[$i]");
		}
		$out->print(",$e1c,$e2c\n");
	}
	$out->Close();
	return 1;
}

sub ChangeDielectricParameters
{
	my ($this, $idx, @param) = @_;
#print "idx=$idx\n";
	my $p = $this->pDielectricModel();
	my $pdm = $p->[$idx];
	for(my $i = 0 ; $i < @param ; $i += 2) {
		$pdm->{$param[$i]} = $param[$i+1];
	}
}

sub AddDielectricModel
{
	my ($this, $name, $model, @param) = @_;
	my $idx = $this->{idxDielectricModel}->{$name};
#print "name=[$name] idx=$idx\n";
	return $this->ChangeDielectricParameters($idx, @param) if(defined $idx);

	my %dm;
	$dm{e1inf} = 0.0;
	$dm{e2inf} = 0.0;
	$dm{fa}    = 0.0;
	$dm{name}  = $name;
	$dm{model} = $model;
	if($model =~ /^Constant$/i) {
		$dm{model} = 'Constant';
	}
	elsif($model =~ /^Tauc.*Lorentz$/i) {
		$dm{model} = 'TaucLorentz';
	}
	elsif($model =~ /^Urbach/i) {
		$dm{model} = 'Urbach';
	}
	elsif($model =~ /^Lorentz/i) {
		$dm{model} = 'Lorentz';
	}
	elsif($model =~ /^Drude/i) {
		$dm{model} = 'Drude';
	}
	elsif($model =~ /^Caucy$/i) {
		$dm{model} = 'Caucy';
	}
	elsif($model =~ /^Sellmeier$/i) {
		$dm{model} = 'Sellmeier';
	}
	elsif($model =~ /^EpsArray$/i) {
		$dm{model} = 'EpsArray';
	}
	elsif($model =~ /^File$/i) {
		$dm{model} = 'File';
	}
	elsif($model =~ /^EFile$/i) {
		$dm{model} = 'EFile';
	}
	elsif($model =~ /^EFileForTranparent$/i) {
		$dm{model} = 'EFileForTranparent';
	}
	elsif($model =~ /^EMA$/i) {
		$dm{model} = 'EMA';
	}
	else {
print "OpticalMaterial::AddDielectricModel: Invalid model [$model].\n";
		exit;
	}

	for(my $i = 0 ; $i < @param ; $i += 2) {
		$dm{$param[$i]} = $param[$i+1];
#print "dm{$param[$i]} = ", $param[$i+1], "\n";
	}

	if($model =~ /^EpsArray$/i) {
		my $pE  = $dm{pE};
#		my $pe1 = $dm{pe1};
#		my $pe2 = $dm{pe2};
		$dm{EStart} = $pE->[0];
		$dm{EEnd}   = $pE->[@$pE - 1];
		$dm{EStep}  = $pE->[1] - $pE->[0];
		$dm{nEData} = @$pE;
	}
	elsif($model =~ /^File$/i) {
		my $path = $dm{Path};
		my $pWL  = $dm{pWL};
#print "p=$path / $pWL\n";
		my ($pe1, $pe2) = CSV::RemakeDataArrayFromFile($path, $pWL, 1);
		$dm{pe1}     = $pe1;
		$dm{pe2}     = $pe2;
		$dm{WLStart} = $pWL->[0];
		$dm{WLEnd}   = $pWL->[@$pWL - 1];
		$dm{WLStep}  = $pWL->[1] - $pWL->[0];
		$dm{nWLData} = @$pWL;
	}
	elsif($model =~ /^EFile$/i or $model =~ /^EFileForTranparent$/i) {
		my $path = $dm{Path};
		my $el = new Ellipsometry;
		my ($nData, $pLabelArray, @DataArray) = $el->Read($path);
#print "Labels: ", join(', ', @$pLabelArray), "\n";
		if(!defined $nData) {
print "Error in OpticalMaterials::AddDielectricModel: Can not read [$path].\n";
			return undef;
		}
		my $pE      = $el->GetData(".*eV.*", ".*E.*");
		my $pe1     = $el->GetData("e1.*");
		my $pe2     = $el->GetData("e2.*");
		if($model =~ /^EFileForTranparent$/i) {
			my $nData = @$pE;
			$pe2 = [];
			for(my $i = 0 ; $i < $nData ; $i++) {
				$pe2->[$i] = 0.0;
			}
		}
if(!defined $pE or !defined $pe1 or !defined $pe2) {
	print "Error in OpticalMaterials::AddDielectricModel: E, e1 and/or e2 can not found in [$path].\n";
	return undef;
}
#print "p=$pE, $pe1, $pe2\n";
		$dm{pe1}    = $pe1;
		$dm{pe2}    = $pe2;
		$dm{EStart} = $pE->[0];
		$dm{EEnd}   = $pE->[@$pE - 1];
		$dm{EStep}  = $pE->[1] - $pE->[0];
		$dm{nEData} = @$pE;
	}
	elsif($model =~ /^EMA$/i) {
		$dm{ScreeningFactor} = 2.0;
		if($dm{Model} eq 'Maxwell-Garnett' or $dm{Model} eq 'MG') {
			$dm{Model}    = 'MG';
			$dm{pEMAFunc} = \&EMAZeroMaxwellGarnett;
		}
		elsif($dm{Model} eq 'Bruggeman') {
			$dm{pEMAFunc} = \&EMAZeroBruggemann;
		}
		elsif($dm{Model} eq 'Average') {
			$dm{pEMAFunc} = \&CalEpsByEMAAverage;
		}
	}

	my $n = $this->nDielectricModel();
	my $p = $this->pDielectricModel();
	$p->[$n] = \%dm;
# $nameのindexを格納
	$this->{idxDielectricModel}->{$name} = $n;
	my $ret = ++$this->{nDielectricModel};
	return $ret;
}

sub CalComplexN
{
	my ($this, $E) = @_;
	my ($e1, $e2) = $this->CalEps($E);
	my ($n, $k)   = Optics::EpsToNK($e1, $e2);
	return cplx($n, -$k);
}

sub CalNK
{
	my ($this, $E) = @_;
	my ($e1, $e2) = $this->CalEps($E);
	return Optics::EpsToNK($e1, $e2);
}

sub CalEps
{
	my ($this, $E) = @_;

	my $n = $this->nDielectricModel();
	my $p = $this->pDielectricModel();
	my ($e1, $e2) = (0.0, 0.0);
	my $S2;
#my $name = $this->Name();
#print "$name\n";
	for(my $i = 0 ; $i < $n ; $i++) {
		my $pModel = $p->[$i];
#print "$name [$i]: $pModel->{model}\n";
		my ($ee1, $ee2);
		if($pModel->{model} eq 'Constant') {
			$ee1 = $pModel->{e1inf};
			$ee2 = $pModel->{e2inf};
		}
		if($pModel->{model} eq 'Constant') {
			($ee1, $ee2) = ($pModel->{e1inf}, $pModel->{e2inf});
		}
		elsif($pModel->{model} =~ /^Tauc.*Lorentz$/i) {
			($ee1, $ee2) = Optics::TaucLorentz(undef, $E, 
								$pModel->{e1inf}, $pModel->{e2inf}, $pModel->{A},
								$pModel->{Eg}, $pModel->{En0}, $pModel->{C});
#print "$E, e=$ee1,$ee2\n";
		}
		elsif($pModel->{model} eq 'Urbach') {
			($ee1, $ee2) = Optics::Urbach(undef, $E, $pModel->{Em}, $pModel->{E0}, $pModel->{A});
		}
		elsif($pModel->{model} eq 'Lorentz') {
			($ee1, $ee2) = Optics::Lorentz(undef, $E, $pModel->{e1inf}, $pModel->{e2inf},
							1.0,
							$pModel->{Ne}, $pModel->{E0}, $pModel->{G});
		}
		elsif($pModel->{model} eq 'Drude') {
			($ee1, $ee2) = Optics::Drude(undef, $E, $pModel->{e1inf}, $pModel->{e2inf},
							$pModel->{Ep}, $pModel->{tau});
		}
		elsif($pModel->{model} eq 'Caucy') {
			my $wl = Optics::eVTonm($E); # nm
#my $pA = $pModel->{pA};
#print "A=$pA->[0],$pA->[1]\n";
			($ee1, $ee2) = Optics::Caucy(undef, $wl, $pModel->{pA}, $pModel->{pn});
		}
		elsif($pModel->{model} eq 'Sellmeier') {
			my $wl = Optics::eVTonm($E); # nm
			($ee1, $ee2) = Optics::Sellmeier(undef, $wl, $pModel->{A}, $pModel->{pB}, $pModel->{pwl0});
		}
		elsif($pModel->{model} eq 'File') {
			my $wl0    = $pModel->{WLStart};
			my $wlstep = $pModel->{WLStep};
			my $wl = Optics::eVTonm($E);
			my $idx = int( ($wl - $wl0) / $wlstep + 0.001 );
			$ee1 = $pModel->{pe1}->[$idx];
			$ee2 = $pModel->{pe2}->[$idx];
		}
		elsif($pModel->{model} eq 'EFile' or $pModel->{model} eq 'EFileForTranparent' or $pModel->{model} eq 'EpsArray') {
			my $E0    = $pModel->{EStart};
			my $Emax  = $pModel->{EEnd};
			my $EStep = $pModel->{EStep};
			my $idx = int( ($E - $E0) / $EStep + 0.001 );
			if($idx < 0 or $idx >= $pModel->{nEData}) {
print "Warning in OpticalMaterial::CalEPS: idx [$idx] out of range [0, $pModel->{nEData}] for E=$E [$E0, $Emax]\n";
				$idx = 0 if($idx < 0);
				$idx = $pModel->{nEData} - 1 if($idx >= $pModel->{nEData});
			}
			$ee1 = $pModel->{pe1}->[$idx];
			$ee2 = $pModel->{pe2}->[$idx];
#print "E=$E ($E0,$EStep): $idx ($ee1,$ee2)\n";
		}
		elsif($pModel->{model} eq 'EMA') {
#			my $Method = "Amoeba::Simplex";
#			my $Method = "PDL::Simplex";
			my $Method = "ModifiedNewton";
			my $EPS         = 1.0e-5;
			my $nMaxIter    = 100;
			my $iPrintLevel = 0;
#print "Optimization method: $Method\n";

			my $EMATheory  = $pModel->{Model};
			my $pMaterials = $pModel->{pMaterials};
			my $pVf        = $pModel->{pVolumeFractions};
			my $nPhase     = @$pMaterials;
#print "nPhase=$nPhase, pVf=", join(', ', @$pVf), "\n";
			if($EMATheory eq 'Average') {
				return $this->CalEpsByEMAAverage($E, $pMaterials, $pVf);
				return $this->CalEpsByEMAAverage($E, $pMaterials, $pVf);
			}

			my @e;
			for(my $j = 0 ; $j < $nPhase ; $j++) {
				my ($e1, $e2) = $pMaterials->[$j]->CalEps($E);
				$e[$j] = cplx($e1, -$e2);
			}
			if(!defined $this->{InitialEps1}) {
				$this->SetInitialEps($pModel->{name}, $E);
			}
			my ($e1h, $e2h) = ($this->{InitialEps1}, $this->{InitialEps2});
			my $optimize = new Optimize;
			$optimize->AddParameters(
					"e1h", \$e1h, 1, 10.0,
					"e2h", \$e2h, 1, 10.0,
					);
			$optimize->SetDumpingFactor(0.0, 0.0, 1);
			my ($OptVars, $MinVal) = $optimize->Optimize(
					$Method, undef, undef, undef,
					$EPS, $nMaxIter, $iPrintLevel,
					sub { $this->CalEMAS2($pModel->{pEMAFunc}, $pModel->{ScreeningFactor}, \@e, $pVf, @_); },
					undef,
					sub { Optimize::BuildDifferentialMatrixes(@_); },
					);
			$ee1 = $OptVars->[0];
			$ee2 = $OptVars->[1];
			if($MinVal > 1.0e-5) {
print "Warning in OpticalMaterial::CalEps: EMA Eps ($ee1,$ee2) at E=$E did not converge (S2=$MinVal).\n";
			}
			$S2 = $MinVal;
			($this->{InitialEps1}, $this->{InitialEps2}) = ($ee1, $ee2);
		}
		else {
			print("Error in OpticalMaterial::CalEps: Invalid model [$pModel->{model}].\n");
			exit;
#			return undef;
		}
		$e1 += $ee1;
		$e2 += $ee2;
#print "i=$i, $ee1, $ee2, $e1, $e2\n";
	}
	return ($e1, $e2, $S2);
}

sub SetInitialEps
{
	my ($this, $ModelName, $E, $e1, $e2) = @_;
	if(!defined $e1) {
		my $n = $this->nDielectricModel();
		my $p = $this->pDielectricModel();
		my $pModel;
		my $pMaterials;
		for(my $i = 0 ; $i < $n ; $i++) {
			$pModel = $p->[$i];
			if($pModel->{name} eq $ModelName) {
				$pMaterials = $pModel->{pMaterials};
				last;
			}
		}
		return undef if(!defined $pModel);

		$pMaterials = $pModel->{pMaterials};
		my $pVf        = $pModel->{pVolumeFractions};
		my $nPhase = @$pMaterials;
		my @e;
		for(my $j = 0 ; $j < $nPhase ; $j++) {
			my ($e1, $e2) = $pMaterials->[$j]->CalEps($E);
			$e[$j] = cplx($e1, -$e2);
		}
		($e1, $e2) = $this->CalInitialEMAEps(\@e, $pVf);
	}
	$this->{InitialEps1} = $e1;
	$this->{InitialEps2} = $e2;
	return ($e1, $e2);
}

sub GetnEMAComponent
{
	my ($this, $ModelName) = @_;

	my $n = $this->nDielectricModel();
	my $p = $this->pDielectricModel();
	for(my $i = 0 ; $i < $n ; $i++) {
		my $pModel = $p->[$i];
		if($pModel->{name} eq $ModelName) {
			my $pMaterials = $pModel->{pMaterials};
			return scalar @$pMaterials;
		}
	}
	return 0;
}

sub GetEMAComponentMaterial
{
	my ($this, $ModelName, $idx) = @_;

	my $n = $this->nDielectricModel();
	my $p = $this->pDielectricModel();
	for(my $i = 0 ; $i < $n ; $i++) {
		my $pModel = $p->[$i];
		if($pModel->{name} eq $ModelName) {
			my $pMaterials = $pModel->{pMaterials};
			return $pMaterials->[$idx];
		}
	}
	return undef;
}

sub GetEMAComponentVolumeFraction
{
	my ($this, $ModelName, $idx) = @_;

	my $n = $this->nDielectricModel();
	my $p = $this->pDielectricModel();
	for(my $i = 0 ; $i < $n ; $i++) {
		my $pModel = $p->[$i];
		if($pModel->{name} eq $ModelName) {
			my $pVf = $pModel->{pVolumeFractions};
			return $pVf->[$idx];
		}
	}
	return undef;
}

sub CalEMAS2
{
	my ($this, $pEMAFunc, $ScreeningFactor, $pe, $pVf, $pVars, $iPrintLevel) = @_;
	my ($e1h, $e2h) = @$pVars;

	my $eh = cplx($e1h, -$e2h);
	my $S = &$pEMAFunc($this, $ScreeningFactor, $eh, $pe, $pVf);
	$S = abs($S);
	my $S2 = $S * $S;
#print "S2=$S2: eh=($eh)\n";
	return $S2;
}

sub CalEpsByEMAAverage
{
	my ($this, $E, $pMaterial, $pVf) = @_;

	my ($e1, $e2) = (0.0, 0.0);
	for(my $i = 0 ; $i < @$pMaterial ; $i++) {
#print "$i: $pVf,$pMaterial: $pVf->[$i],$pMaterial->[$i]\n";
		my ($e11, $e21) = $pMaterial->[$i]->CalEps($E);
#print "    e1=$e1\n";
		$e1 += $pVf->[$i] * $e11;
		$e2 += $pVf->[$i] * $e21;
	}
	return ($e1, $e2);
}

sub EMAZeroBruggemann
{
	my ($this, $ScreeningFactor, $eh, $pe, $pVf) = @_;
	my $nPhase = @$pe;
	my $S = 0.0;
	for(my $i = 0 ; $i < $nPhase ; $i++) {
		$S += $pVf->[$i] * ($pe->[$i] - $eh) / ($pe->[$i] + $ScreeningFactor*$eh);
	}
	return $S;
}

sub EMAZeroMaxwellGarnett
{
	my ($this, $ScreeningFactor, $eh, $pe, $pVf) = @_;
	my $nPhase = @$pe;

	my $S = ($eh - $pe->[0]) / ($eh + $ScreeningFactor*$pe->[0]);
	for(my $i = 1 ; $i < $nPhase ; $i++) {
		$S -= $pVf->[$i] * ($pe->[$i] - $pe->[0]) / ($pe->[$i] + $ScreeningFactor*$pe->[0]);
	}
	return $S;

#	if($nPhase != 2) {
#		print "Warning in EMAZeroMaxwellGarnett: nPhase [$nPhase] should be 2.\n"
#	}
#	my $S = ($eh - $pe->[0]) / ($eh + $ScreeningFactor*$pe->[0])
#		  - $pVf->[1] * ($pe->[1] - $pe->[0]) / ($pe->[1] + $ScreeningFactor*$pe->[0]);
#	return $S;
}

sub CalInitialEMAEps
{
	my ($this, $pe, $pVf) = @_;
	my $nPhase = @$pe;
	my ($e1h, $e2h) = (0.0, 0.0);
	for(my $i = 0 ; $i < $nPhase ; $i++) {
		$e1h +=  Re($pe->[$i]) * $pVf->[$i];
		$e2h += -Im($pe->[$i]) * $pVf->[$i];
	}
if(0) {
	for(my $i = 0 ; $i < $nPhase ; $i++) {
		$e1h +=  Re($pe->[$i]);
		$e2h += -Im($pe->[$i]);
	}
	$e1h /= $nPhase;
	$e2h /= $nPhase;
}
	return ($e1h, $e2h);
}
1;