#===============================================
# Ellipsometry
#===============================================
package Ellipsometry;
use MultiColumnData;
use Sci::Optics;
@ISA = qw(MultiColumnData Optics);

use strict;
use Math::Complex;

use Sci qw($pi $e0 $e $c $h $hbar $me $kB $todeg $torad);
my $pi2  = $pi + $pi;

#============================================================
#　静的関数
#============================================================

#============================================================
#　変数等取得関数
#============================================================

#============================================================
#　コンストラクタ、デストラクタ
#============================================================
BEGIN { }

sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;

	$this->Initialize();
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　一般関数
#============================================================
sub Read
{
	my ($this, $filename, $IsNumerical, $YLabel) = @_;

	my $ext = Deps::ExtractExtention($filename);
	if($ext =~ /\.ref$/i or $ext =~ /\.asp$/i or $ext =~ /\.isa$/i) {
		my ($nData, $pLabelArray, @DataArray) = $this->ReadASP($filename, $YLabel);
#		return ($nData, $pLabelArray, @DataArray);

		$this->GetXData("E.*");
		$this->GetYData("T.*");
		my $pSorted = $this->SortXYByX(1);
		$nData = $this->nData(0);
		my ($x0, $xend, $xstep) = (0.0, 100.0, 0.05);
		my ($xmin, $xmax, $ymin, $ymax) = $pSorted->CalMinMax();
		my $func = \&Algorism::Interpolate;
#		my $func = \&Algorism::InterpolateByCubicSpline;
		my $pX = $pSorted->pX();
		my $pY = $pSorted->pY();
		my (@E, @T);
		my $c = 0;
		for(my $i = 0 ; $i < $nData ; $i++) {
			my $xx = $x0 + $xstep * $i;
			last if($xx > $xend);
			next if($xx <$xmin or $xmax < $xx);

			$E[$c] = $xx;
			$T[$c] = &$func($pX, $pY, $xx);
			$c++;
		}

		$this->SetTitle(Deps::ExtractFileBody($filename));
		$this->SetFileType("TEXT");
		$this->SetLabelArray(["E(eV)", "T"]);
		$this->SetDataArray([\@E, \@T]);
		return ($c, ["E(eV)", "T"], [\@E, \@T]);
	}
	if($ext =~ /\.txt$/i) {
		my ($nData, $pLabelArray, @DataArray) = $this->ReadSPring8NIMSHXPES($filename, $YLabel);
		my $pX = $this->GetXData("E.*");
		my $pY = $this->GetYData("I.*");
print "pX=$pX  /  pY=$pY\n";
#		return ($nData, $pLabelArray, @DataArray);
		return ($nData, ['E(eV)', 'I'], [$pX, $pY]);
	}
	if($ext =~ /\.txt$/i) {
		my $csv = new CSV;
		if(!$csv->ReadFreeFormat($filename, "[\\s,]+", ["wl(nm)", "T"], 1, 1, 2)) {
print "Error in Ellipsometry::Read [Free Format Text]: Can not read [$filename].\n";
			$this->SetTitle('');
			$this->SetFileType('');
			return undef;
		}
		my $pwl = $csv->GetXData("wl.*");
		$csv->AddDataByConversion("E(eV)", $pwl, sub { my ($x)=@_; return Optics::nmToeV($x); });
#		my $nDataArray = $csv->nDataArray();
#print "nDataArray=$nDataArray\n";
#		my $pLabelArray = $csv->LabelArray();
#print "LabelArray=", join(', ', @$pLabelArray), "\n";
		$csv->GetXData("E.*");
		$csv->GetYData("T.*");
		my $pSorted = $csv->SortXYByX(1);
		my $nData = $csv->nData(0);
		my ($x0, $xend, $xstep) = (0.0, 100.0, 0.05);
		my ($xmin, $xmax, $ymin, $ymax) = $pSorted->CalMinMax();
		my $func = \&Algorism::Interpolate;
#		my $func = \&Algorism::InterpolateByCubicSpline;
		my $pX = $pSorted->pX();
		my $pY = $pSorted->pY();
		my (@E, @T);
		my $c = 0;
		for(my $i = 0 ; $i < $nData ; $i++) {
			my $xx = $x0 + $xstep * $i;
			last if($xx > $xend);
			next if($xx <$xmin or $xmax < $xx);

			$E[$c] = $xx;
			$T[$c] = &$func($pX, $pY, $xx);
			$c++;
		}

		$this->SetTitle(Deps::ExtractFileBody($filename));
		$this->SetFileType("TEXT");
		$this->SetLabelArray(["E(eV)", "T"]);
		$this->SetDataArray([\@E, \@T]);
		return ($c, ["E(eV)", "T"], [\@E, \@T]);
	}
	if($ext =~ /\.csv$/i) {
		$this->SetIncidentAngle(0.0);
		my ($nData, $pLabelArray, @DataArray) = CSV::GetArraysFromFileSimple($filename, undef, 
								$this->SkipBlankData());
#		my ($nData, $pLabelArray, @DataArray) = CSV::GetArraysFromFile($filename);
		if(!defined $nData) {
print "Error in Ellipsometry::Read [CSV]: Can not read [$filename].\n";
			$this->SetTitle('');
			$this->SetFileType('');
			return undef;
		}
		$this->SetTitle(Deps::ExtractFileBody($filename));
		$this->SetFileType("SPE");
		$this->SetLabelArray($pLabelArray);
		$this->SetDataArray(\@DataArray);
		return ($nData, $pLabelArray, @DataArray);
	}

	my ($nData, $pLabelArray, @DataArray) = $this->SUPER::Read($filename, $IsNumerical, 'T');
	return ($nData, $pLabelArray, @DataArray);
}

sub EpsToNK
{
	my ($this, $e1, $e2) = @_;
	return Optics::EpsToNK($e1, $e2);
}

sub NKToEps
{
	my ($this, $n, $k) = @_;
	return Optics::NKToEps($n, $k);
}

# P-M=45, A=+-45
# Configuration II: M = +-90
# Configuration III: +-45
sub PsiDeltaTolIsIc
{
	my ($this, $Psi, $Delta, $PmM, $A, $M) = @_;
	my $Psi2 = $Psi + $Psi;
	my $A2   = $A + $A;
	my $M2   = $M + $M;
	my $sin2PmM  = sin(($PmM+$PmM)*$torad);
	my $sin2A    = sin($A2*$torad);
	my $cos2A    = cos($A2*$torad);
	my $sin2M    = sin($M2*$torad);
	my $cos2M    = cos($M2*$torad);
	my $sin2Psi  = sin($Psi2*$torad);
	my $cos2Psi  = cos($Psi2*$torad);
	my $sinDelta = sin($Delta*$torad);
	my $cosDelta = cos($Delta*$torad);
# Mの符号が変
print "2A=$sin2A, $cos2A (A2=$A2)\n";
print "2M=$sin2M, $cos2M (M2=$M2)\n";

	my $Is = $sin2PmM * $sin2A * $sin2Psi * $sinDelta;
	my $Ic = $sin2PmM * ($sin2M * ($cos2Psi - $cos2A) + $sin2A * $cos2M * $sin2Psi * $cosDelta);
	return ($Is, $Ic);
}

# A: +1 for A=+45, -1 for A=-45
# M: +1 for M=+90, -1 for M=-90
sub PsiDeltaToIsIcConfigrationII
{
	my ($this, $Psi, $Delta, $signA, $signM) = @_;
	$signA = 1 if(!defined $signA);
	$signM = 1 if(!defined $signM);

	my $Psi2 = $Psi + $Psi;
	my $sin2Psi  = sin($Psi2*$torad);
	my $sinDelta = sin($Delta*$torad);
	my $cosDelta = cos($Delta*$torad);
	
	my $Is = $signA * $sin2Psi * $sinDelta;
	my $Ic = $signM * $signA * $sin2Psi * $cosDelta;
	return ($Is, $Ic);
}

sub IsIcToPsiDeltaConfigrationII
{
	my ($this, $Is, $Ic, $signA, $signM) = @_;
	$signA = 1 if(!defined $signA);
	$signM = 1 if(!defined $signM);

	my $Delta = atan2($Is, $Ic);
	my $Psi   = asin($signA * $Is / sin($Delta)) * 0.5;
	return ($Psi*$todeg, $Delta*$todeg);
}

# A: +1 for A=+45, -1 for A=-45
# M: +1 for M=+45, -1 for M=-45
sub PsiDeltaToIsIcConfigrationIII
{
	my ($this, $Psi, $Delta, $signA, $signM) = @_;
	$signA = 1 if(!defined $signA);
	$signM = 1 if(!defined $signM);

	my $Psi2 = $Psi + $Psi;
	my $sin2Psi  = sin($Psi2*$torad);
	my $cos2Psi  = cos($Psi2*$torad);
	my $sinDelta = sin($Delta*$torad);
	
	my $Is = $signA * $sin2Psi * $sinDelta;
	my $Ic = $signM * $cos2Psi;
	return ($Is, $Ic);
}

sub IsIcToPsiDeltaConfigrationIII
{
	my ($this, $Is, $Ic, $signA, $signM) = @_;
	$signA = 1 if(!defined $signA);
	$signM = 1 if(!defined $signM);

	my $Psi   = acos($signM * $Ic) * 0.5;
	my $Delta = asin($signA * $Is / sin($Psi+$Psi));
	return ($Psi*$todeg, $Delta*$todeg);
}

sub MultMatrix
{
	my ($this, $pM1, $pM2) = @_;
#print "M1=($pM1->{11}, $pM1->{12})\n";
#print "   ($pM1->{21}, $pM1->{22})\n";
#print "M2=($pM2->{11}, $pM2->{12})\n";
#print "   ($pM2->{21}, $pM2->{22})\n";
	my %M;
	$M{11} = $pM1->{11}*$pM2->{11} + $pM1->{12}*$pM2->{21};
	$M{12} = $pM1->{11}*$pM2->{12} + $pM1->{12}*$pM2->{22};
	$M{21} = $pM1->{21}*$pM2->{11} + $pM1->{22}*$pM2->{21};
	$M{22} = $pM1->{21}*$pM2->{12} + $pM1->{22}*$pM2->{22};
#print "Mt=($M{11}, $M{12})\n";
#print "   ($M{21}, $M{22})\n";
	return %M;
}

sub CharacteristicMatrixToFresnelCoefficient
{
	my ($this, $Polarization, $pM, $nIncident, $kIncident, $nTransmission, $kTransmission, $QIncident, $QTransmission) = @_;
	my $NIncident        = cplx($nIncident, -$kIncident);
	my $NTransmission    = cplx($nTransmission, -$kTransmission);
#print "N=$NIncident,$NTransmission\n";
#print "Q=$QIncident, $QTransmission\n";
	my $cosQIncident     = cos($QIncident    *$torad);
	my $cosQTransmission = cos($QTransmission*$torad);
	my $eIncident;
	my $eTransmission;
	if($Polarization eq 's') {
		$eIncident     = $NIncident     * $cosQIncident;
		$eTransmission = $NTransmission * $cosQTransmission;
	}
	elsif($Polarization eq 'p') {
		$eIncident     = $NIncident     / $cosQIncident;
		$eTransmission = $NTransmission / $cosQTransmission;
	}
	else {
print "Error in Ellipsometry::CharacteristicMatrixToFresnelCoefficient: Invalid Polarization [$Polarization].\n";
		return undef;
	}
#print "eta=$eIncident,$eTransmission\n";

	my $m11 = $pM->{11};
	my $m12 = $pM->{12};
	my $m21 = $pM->{21};
	my $m22 = $pM->{22};
#print "M=($m11, $m12)\n";
#print "  ($m21, $m22)\n";

#	my $rhoUpper = $eIncident*$m11 - $eTransmission*$m22 + i * ($eIncident*$eTransmission*$m12 - $m21);
#	my $Lower    = $eIncident*$m11 + $eTransmission*$m22 + i * ($eIncident*$eTransmission*$m12 + $m21);
	my $rhoUpper = $eIncident*$m11 + $eIncident*$eTransmission*$m12 - $m21 - $eTransmission*$m22;
	my $Lower    = $eIncident*$m11 + $eIncident*$eTransmission*$m12 + $m21 + $eTransmission*$m22;
	my $rho = $rhoUpper / $Lower;
	my $tauUpper = 2.0 * $eIncident;
	my $tau = $tauUpper / $Lower;
	my $R = abs($rho);
	$R = $R * $R;
	my $T = abs($tau);
	$T = Re($eTransmission) / Re($eIncident) * $T * $T;
	return ($rho, $tau, $R, $T);
}

sub CalCharacteristicMatrix
{
# thickness: m
	my ($this, $Polarization, $E, $n, $k, $thickness, $Angle) = @_;

	my %M;
	if($Polarization eq 'unit') {
		$M{11} = 1.0;
		$M{12} = 0.0;
		$M{21} = 0.0;
		$M{22} = 1.0;
		return %M;
	}

	my $wl = Optics::eVTonm($E) * 1.0e-9; # m
#print "Angle=$Angle\n";
	my $cosAngle  = cos($Angle*$torad);
	my $N         = cplx($n, -$k);
	my $Delta;
	if($thickness > 0.0) {
		$Delta = $pi2 / $wl * $N * $thickness * $cosAngle;
	}
	else {
		$Delta = 0.0;
	}
#print "Delta=$Delta\n";
	my $cosDelta  = cos($Delta);
	my $sinDelta  = sin($Delta);

	my $eta;
	if($Polarization eq 's') {
		$eta = $N * $cosAngle;
	}
	elsif($Polarization eq 'p') {
		$eta = $N / $cosAngle;
	}
	else {
print "Error in Ellipsometry::CalCharacteristicMatrix: Invalid Polarization [$Polarization].\n";
		return undef;
	}

#print "$Polarization: d=$Delta  N=($N)\n";
	$M{11} = $cosDelta;
	$M{12} = i * $sinDelta / $eta;
	$M{21} = i * $sinDelta * $eta;
	$M{22} = $cosDelta;
#print "Mc=($M{11}, $M{12})\n";
#print "   ($M{21}, $M{22})\n";

	return %M;
}

sub ReflectionFresnelCoefficientToPsiDelta
{
	my ($this, $rs, $rp) = @_;

	if($rs == 0.0) {
		return (90.0, 0.0);
	}
	if($rp == 0.0) {
		return (0.0, 0.0);
	}

	my $rho    = $rp / $rs;
	my $tanPsi = abs($rho);
	my $Psi    = atan2($tanPsi, 1.0) * $todeg;
	my $Delta  = arg($rho) * $todeg;
#	my $Delta  = -arg($rho) * $todeg;
#	my $rr     = Re($rho);
#	my $ri     = Im($rho);
#	my $Delta  = -atan2($ri, $rr) * $todeg;
#print "Psi=$Psi, $Delta\n";
	$Delta += 360.0 if($Delta < 0.0);
	return ($Psi, $Delta);
}

sub SnellsRawComplex
{
	my ($this, $nIncident, $kIncident, $nSubstrate, $kSubstrate, $AngleIncident) = @_;
	$AngleIncident = $this->IncidentAngle() if(!defined $AngleIncident);
	my $NIncident      = cplx($nIncident,  -$kIncident);
	my $NSubstrate     = cplx($nSubstrate, -$kSubstrate);
	$NSubstrate = 1.0e-3 if($NSubstrate == 0.0);
#print "N=$NIncident / $NSubstrate\n";
	my $sinQ = sin($AngleIncident*$torad) * $NIncident / $NSubstrate;
	my $AngleSubstrate = Sci::asin($sinQ) * $todeg;
	return $AngleSubstrate;
}

#sub SnellsRaw
#{
#	my ($this, $nIncident, $nSubstrate, $AngleIncident) = @_;
#	$AngleIncident = $this->IncidentAngle() if(!defined $AngleIncident);
#	my $sinQ = sin($AngleIncident*$torad) * $nIncident / $nSubstrate;
#	my $AngleSubstrate = Sci::asin($sinQ) * $todeg;
#	return $AngleSubstrate;
#}

# rpの符号が逆？
sub CalFresnelCoefficients
{
	my ($this, $nIncident, $kIncident, $nSubstrate, $kSubstrate, $AngleIncident) = @_;
	$AngleIncident = $this->IncidentAngle() if(!defined $AngleIncident);
	my $NIncident      = cplx($nIncident, -$kIncident);
	my $cosQIncident   = cos($AngleIncident*$torad);
	my $NSubstrate     = cplx($nSubstrate, -$kSubstrate);
	my $AngleSubstrate = SnellsRawComplex($this, $nIncident, $kIncident, $nSubstrate, $kSubstrate, $AngleIncident);
	my $cosQSubstrate  = cos($AngleSubstrate*$torad);
#print "AngleIncident=$AngleIncident\n";
#print "AngleSubstrate=$AngleSubstrate\n";

	my $ap = ($NSubstrate * $cosQIncident + $NIncident * $cosQSubstrate);
	my $rp = ($NIncident * $cosQSubstrate - $NSubstrate * $cosQIncident) / $ap;
	my $tp = 2.0 * $NIncident * $cosQSubstrate / $ap;
	my $as = ($NIncident * $cosQIncident + $NSubstrate * $cosQSubstrate);
	my $rs = ($NIncident * $cosQIncident - $NSubstrate * $cosQSubstrate) / $as;
	my $ts = 2.0 * $NIncident * $cosQIncident / $as;
	return ($rs, $ts, -$rp, $tp);
}

sub PsiDeltaToNK
{
	my ($this, $Psi, $Delta) = @_;
	my ($e1, $e2) = $this->PsiDeltaToEps($Psi, $Delta);
	return Optics::EpsToNK($e1, $e2);
}

sub PsiDeltaToEps
{
	my ($this, $Psi, $Delta, $Angle) = @_;
	$Angle = $this->IncidentAngle() if(!defined $Angle);

#	my ($n, $k) = $this->PsiDeltaToNK_OldVersion($Psi, $Delta, $Angle);
#	my ($e1, $e2) = $this->NKToEps($n, $k);
#	return ($e1, $e2);

	my $tanPsi = Sci::tan($Psi*$torad);
	my $expiDelta = exp(cplx(0.0, $Delta*$torad));
#	my $expiDelta = exp(cplx(0.0, -$Delta*$torad));
	my $rho = $tanPsi * $expiDelta;

	my $sinQ     =      sin($Angle*$torad);
	my $tanQ     = Sci::tan($Angle*$torad);
	my $a = (1.0 - $rho) / (1.0 + $rho);
	my $eps = $sinQ*$sinQ * (1.0 + $tanQ*$tanQ * $a*$a);

	return (Re($eps), -Im($eps));
#	return (Re($eps), Im($eps));
}

#実数を使った方法。PsiDeltaToEpsと同じ結果がでるはず
sub PsiDeltaToNK_OldVersion
{
	my ($this, $Psi, $Delta, $Angle) = @_;
	$Angle = $this->IncidentAngle() if(!defined $Angle);
	my $sinQ     =      sin(  $Angle*$torad);
	my $tanQ     = Sci::tan(  $Angle*$torad);
	my $sin2Psi  =      sin(2.0*$Psi*$torad);
	my $cos2Psi  =      cos(2.0*$Psi*$torad);
	my $tan2Psi  = Sci::tan(2.0*$Psi*$torad);
	my $sinDelta =      sin(  $Delta*$torad);
	my $cosDelta =      cos(  $Delta*$torad);

	my $A = (1.0 + $sin2Psi*$cosDelta);
	my $G2 = $sinQ*$sinQ * $tanQ*$tanQ * $cos2Psi*$cos2Psi / $A / $A;
	my $n2mk2 =  $sinQ*$sinQ + (1.0 + $tan2Psi*$tan2Psi * $sinDelta*$sinDelta) * $G2;
	my $nk2   = 2.0 * $G2 * $tan2Psi * $sinDelta;

	my $n2 = ($n2mk2 + sqrt($n2mk2*$n2mk2 + $nk2*$nk2)) / 2.0;
	my $n = sqrt($n2);
	my $k = $nk2 / 2.0 / $n;

#my $n2mk2c = $n*$n - $k*$k;
#my $nk2c = 2.0 * $n*$k;
#print "n2-k2: $n2mk2, $n2mk2c  2nk: $nk2, $nk2c\n";

	return ($n, $k);
}

#近似解
sub ApproximatePsiDeltaToNK
{
	my ($this, $Psi, $Delta, $Angle) = @_;
	$Angle = $this->IncidentAngle() if(!defined $Angle);
	my $sinQ     = sin($Angle*$torad);
	my $tanQ     = Sci::tan($Angle*$torad);
	my $sin2Psi  = sin(2.0*$Psi*$torad);
	my $cos2Psi  = cos(2.0*$Psi*$torad);
	my $tan2Psi  = Sci::tan(2.0*$Psi*$torad);
	my $sinDelta = sin($Delta*$torad);
	my $cosDelta = cos($Delta*$torad);
	my $tanDelta = $sinDelta / $cosDelta;

	my $n = $sinQ * $tanQ * $cos2Psi / (1.0 + $cosDelta * $sin2Psi);
	my $k = $n * $tanDelta * $tan2Psi;
	return ($n, $k);
}

1;