#===============================================
# MultiLayer
#===============================================
package MultiLayer;
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 qw($pi $pi2 $e0 $e $c $h $hbar $me $kB $torad $todeg);
use Sci::Optics;

sub nLayer           { return shift->{nLayer}; }
sub pLayerArray      { return shift->{pLayerArray}; }
sub pLayer           { my($this,$idx)=@_; return $this->{pLayerArray}->[$idx]; }
sub SetIncidentAngle { my($this,$a)=@_; return $this->{IncidentAngle} = $a; }
sub IncidentAngle    { my($this)=@_; return $this->{IncidentAngle}; }
sub SetVariable
{
	my ($this, $idx, $VarName, $val) = @_;
	my $pLayer = $this->pLayer($idx);
	return $pLayer->SetVariable($VarName, $val);
}

sub pVariable
{
	my ($this, $idx, $VarName) = @_;
	my $pLayer = $this->pLayer($idx);
	if(!defined $pLayer) {
		print "MultiLayer::pVariable: Undefined pLayer for index=$idx.\n";
		return undef;
	}
	my $p = $pLayer->pVariable($VarName);
	if(!defined $p) {
		print "MultiLayer::pVariable: Undefined Variable [$VarName] for index=$idx.\n";
		return undef;
	}
	return $p;
}

sub SetCoherency
{
	my ($this, $idx, $Coherency) = @_;
	my $pLayer = $this->pLayer($idx);
#	$pLayer->SetCoherency($Coherency);
	return $this->{"Coherency$idx"} = $Coherency;
}

sub Coherency
{
	my ($this, $idx) = @_;
	return $this->{"Coherency$idx"};
}

sub SetLayerThickness
{
	my ($this, $idx, $Thickness) = @_;
	my $pLayer = $this->pLayer($idx);
#	$pLayer->SetThickness($Thickness);
	return $this->{"Thickness$idx"} = $Thickness;
}

sub LayerThickness
{
	my ($this, $idx) = @_;
	return $this->{"Thickness$idx"};
}

sub pLayerThickness
{
	my ($this, $idx) = @_;
	my $p = \$this->{"Thickness$idx"};
#my $d = $this->LayerThickness($idx);
#print "p=$p, $$p ($d) [idx=$idx]\n";
	return $p;
}

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

sub DESTROY
{
	my $this = shift;
}

sub Initialize
{
	my ($this) = @_;
	$this->{nLayer} = 0;
	$this->{pLayerArray} = [];
}

sub LightAngle
{
	my ($this, $idx, $E, $Angle0) = @_;
	$Angle0 = $this->IncidentAngle() if(!defined $Angle0);
	if(!defined $Angle0) {
print "MultiLayer: Error!!: IncidentAngle not specified.\n";
		return undef;
	}

	my $pLayer0 = $this->pLayer(0);
	my ($n0, $k0) = $pLayer0->CalNK($E);
	my $pLayer = $this->pLayer($idx);
	my ($n, $k) = $pLayer->CalNK($E);
	my $Angle = Ellipsometry::SnellsRawComplex(undef, 1.0, 0.0, $n, $k, $Angle0);
	return $Angle;
}

sub AddLayer
{
# $pMaterial = Object: OpticalMaterial
# $Thickness: m
# $Coherency: 0.0 ~ 1.0 (valid for 0 and non-zero at present)
	my ($this, $pMaterial, $Thickness, $Coherency) = @_;
	$Coherency = 1.0 if(!defined $Coherency);

	my $n = $this->{nLayer};
	$this->{pLayerArray}->[$n] = $pMaterial;
#print "n=$n c=$Coherency\n";
	if(defined $Thickness) {
		$this->SetLayerThickness($n, $Thickness);
#		$pMaterial->SetThickness($Thickness);
#		$this->{"Thickness$n"} = $Thickness;
	}
	$this->SetCoherency($n, $Coherency);
#	$pMaterial->SetCoherency($Coherency);
#	$this->{"Coherency$n"} = $Coherency;
	return $this->{nLayer}++;
}

#すべての層のe1,e2をリストで返す
sub CalEps
{
	my ($this, $E) = @_;
	my $n = $this->nLayer();
	my $pLayerArray = $this->pLayerArray();
	my @e;
	for(my $i = 0 ; $i < $n*2 ; $i += 2) {
		my $pLayer = $pLayerArray->[$i/2];
		($e[$i], $e[$i+1]) = $pLayer->CalEps($E);
	}
	return @e;
}

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

	my $pLayer  = $this->pLayer($idx);
	my ($n, $k) = $pLayer->CalNK($E);
	my $N       = cplx($n, -$k);
	my $Angle   = $this->LightAngle($idx, $E);
	my $cosQ    = cos($Angle*$torad);
	my $etas    = $N * $cosQ;
	my $etap    = $N / $cosQ;
	return ($etas, $etap, $N, $n, $k, $Angle, $cosQ);
}

sub CalPsiDelta
{
	my ($this, $E) = @_;
	my ($rs, $rp)     = $this->CalFresnelCoefficient($E);
	my ($Psi, $Delta) = Ellipsometry::ReflectionFresnelCoefficientToPsiDelta(undef, $rs, $rp);
	return ($Psi, $Delta);
}

sub DecomposeLayersByCoherency
{
	my ($this) = @_;
	
	my $Angle = $this->IncidentAngle();
	my @LayersBlock;
	my $nl = 0;
#最初の層は特別
	my $pl = $this->pLayer(0);
	my $d  = $this->LayerThickness(0);
	my $c  = $this->Coherency(0);
	$LayersBlock[$nl] = new MultiLayer();
	$LayersBlock[$nl]->SetIncidentAngle($Angle);
	$LayersBlock[$nl]->AddLayer($pl, $d, $c);

	my $nlayer = $this->nLayer();
	my $prevc;
	for(my $i = 1 ; $i < $nlayer - 1 ; $i++) {
		my $pl = $this->pLayer($i);
		my $d  = $this->LayerThickness($i);
		my $c  = $this->Coherency($i);
		if(!defined $prevc) {
			$prevc = $c;
		}
		if($c != $prevc) {
			$LayersBlock[$nl]->AddLayer($pl, $d, $c);
			$nl++;
			$prevc = $c;
			$LayersBlock[$nl] = new MultiLayer();
			$LayersBlock[$nl]->SetIncidentAngle($Angle);
		}
		$LayersBlock[$nl]->AddLayer($pl, $d, $c);
	}

	$pl = $this->pLayer($nlayer-1);
	$d  = $this->LayerThickness($nlayer-1);
	$c  = $this->Coherency($nlayer-1);
	$LayersBlock[$nl]->AddLayer($pl, $d, $c);

if(0) {
	for(my $i = 0 ; $i < @LayersBlock ; $i++) {
#		my $pl = $LayersBlock[$i]->pLayer(0);
		my $c  = $LayersBlock[$i]->Coherency(0);
		print "***Layers $i: $c\n";
		$LayersBlock[$i]->PrintLayers();
	}
	exit;
}
	return \@LayersBlock;
}

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

	my ($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp);
	my $pLayersBlock = $this->DecomposeLayersByCoherency();
	my $nBlock = @$pLayersBlock;
#print "nBlock=$nBlock\n";
	if($nBlock == 1) {
		($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp)
			= $this->CalReflectionFresnelCoefficient_FromBottom($E);
	}
	elsif($nBlock == 2) {
		($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp)
			= $pLayersBlock->[0]->CalReflectionFresnelCoefficient_FromBottom($E);
		($Rs, $Rp, $Ts, $Tp) = $this->CalIncoherntRT($E, $Rs, $Rp, $Ts, $Tp, $pLayersBlock->[1]);
	}
	else {
		print "Error in MultiLayer::CalFresnelCoefficient: # of Layers blocks [$nBlock] is not supported.\n";
		exit;
	}

#上部から計算し、下部反射率を正確に計算していないので間違えている。
#	($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp);
#			= $this->CalReflectionFresnelCoefficient_FromTop($E);
#下部から計算して正しい値を出す
#	($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp);
#			= $this->CalReflectionFresnelCoefficient_FromBottom($E);
#CalReflectionFresnelCoefficient_FromBottomと同じ正しい値を出す
#	($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp);
#			= $this->CalFresnelCoefficient_ByCharacteristicMatrix($E);

	$rp23 = -$rp23;
	return ($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp);
}

sub CalInterfaceFresnelCoefficient
{
# 光の入射は $idx1 =>  $idx2 の方向
	my ($this, $E, $idx1, $idx2) = @_;

	my ($e1s, $e1p, $N1, $n1, $k1, $Angle1, $cosQ1) = $this->CalEta($E, $idx1);
	$e1s = 1.0e-3 if(Re($e1s) == 0.0);
	$e1p = 1.0e-3 if(Re($e1p) == 0.0);
	my ($e2s, $e2p, $N2, $n2, $k2, $Angle2, $cosQ2) = $this->CalEta($E, $idx2);

	my $ks = 1.0 / ($e1s + $e2s);
	my $kp = 1.0 / ($e1p + $e2p);
	my $r12s = ($e1s - $e2s) * $ks;
	my $r12p = ($e1p - $e2p) * $kp;
	my $t12s = ($e1s + $e1s) * $ks;
	my $t12p = ($e1p + $e1p) * $kp;

	my $R12s = abs($r12s);
	$R12s = $R12s * $R12s;
	my $T12s = abs($t12s);
	$T12s = Re($e2s) / Re($e1s) * $T12s * $T12s;

	my $R12p = abs($r12p);
	$R12p = $R12p * $R12p;
	my $T12p = abs($t12p);
	$T12p = Re($e2p) / Re($e1p) * $T12p * $T12p;

	return ($r12s, $r12p, $t12s, $t12p, $R12s, $R12p, $T12s, $T12p);
}

#idx層中を端から端まで片道走行する際の光学距離から求めた位相変化
sub PhaseShift
{
	my ($this, $E, $idx) = @_; 
	my $wl = Optics::eVTonm($E) * 1.0e-9; # m
	my $pLayer = $this->pLayer($idx);
	my $d      = $this->LayerThickness($idx); # m
#print "idx=$idx d=$d  name=", $pLayer->Name(), "\n";
if($d < 0) {
	print "Error: MultiLayer::PhaseShift: Negative thickness ($d) for idx=$idx\n";
#	return undef;
}
	my $N     = $pLayer->CalComplexN($E);
	my $Angle = $this->LightAngle($idx, $E);
	my $cosQ  = cos($Angle * $torad);
# 膜には干渉を考慮した位相変化を、基板の位相変化には干渉を考慮せず実際の光路長を使った位相変化を返す
	my $k = $pi2 / $wl * $d * $N;
	my $PSfilm      = $k * $cosQ;
	my $PSsubstrate = $k / $cosQ;
#print "idx=$idx d=$d  $PSfilm\n";
	return ($PSfilm, $PSsubstrate);
}

sub CalIncoherntRT
{
	my ($this, $E, $Rs01, $Rp01, $Ts01, $Tp01, $pLayersBlock) = @_;

	my $n =  $pLayersBlock->nLayer();
#print "n=$n: Rs=$Rs01 => ";
	my ($Rs012, $Rp012, $Ts012, $Tp012);
	for(my $i = 0 ; $i < $n - 1 ; $i++) {
#print "i=$i\n";
		my ($rs12, $rp12, $ts12, $tp12, $Rs12, $Rp12, $Ts12, $Tp12)
				= $pLayersBlock->CalInterfaceFresnelCoefficient($E, $i, $i+1);
		my ($b2, $B2) = $pLayersBlock->PhaseShift($E, $i);
# 振幅の位相変化B2から、エネルギー透過係数Tiをだすため、B2は2倍する
		$B2 = $B2 + $B2;
		my $Ti = abs(exp(-$B2 * i));
		my $Ti2 = $Ti * $Ti;
		my $Ks = 1.0 / (1.0 - $Rs01 * $Rs12 * $Ti2);
		my $Kp = 1.0 / (1.0 - $Rp01 * $Rp12 * $Ti2);
		$Rs012 = $Rs01 + $Ks * $Ts01 * $Rs12 * $Ts01 * $Ti2;
		$Rp012 = $Rp01 + $Kp * $Tp01 * $Rp12 * $Tp01 * $Ti2;
		$Ts012 = $Ks * $Ts01 * $Ts12 * $Ti;
		$Tp012 = $Kp * $Tp01 * $Tp12 * $Ti;
		$Rs01 = $Rs012;
		$Rp01 = $Rp012;
		$Ts01 = $Ts012;
		$Tp01 = $Tp012;
	}
#print " $Rs012\n";
	return ($Rs012, $Rp012, $Ts012, $Tp012);
}

#仮想界面を使って計算
sub CalReflectionFresnelCoefficient_FromBottom
{
	my ($this, $E) = @_;
	my $wl = Optics::eVTonm($E) * 1.0e-9; # m
	my $nLayer = $this->nLayer();
	if($nLayer <= 1) {
print "Note: MultiLayer::CalReflectionFresnelCoefficient: "
	 ."Use Ellipsometry::CalFresnelCoefficients for nLayer=$nLayer\n";
		my ($nsub, $ksub) = $this->pLayer(0)->CalNK($E);
		my ($rs, $ts, $rp, $tp) 
			= Ellipsometry::CalFresnelCoefficients(undef, 1.0, 0.0, $nsub, $ksub, $this->IncidentAngle());
		return ($rs, $rp, $rp, $tp);
	}

# 光の入射は $idx1=$idx0-1  =>  $idx0 の方向ｓ
# 一番下の部層
	my $idx0 = $nLayer-1;
# 下の界面でのFresnel係数
	my ($rs23, $rp23, $ts23, $tp23, $Rs23, $Rp23, $Ts23, $Tp23)
			= $this->CalInterfaceFresnelCoefficient($E, $idx0-1, $idx0);
	my ($esBottom, $epBottom) = $this->CalEta($E, $idx0);

# idx: 考えている、光が透過する層
	for(my $idx = $nLayer-2 ; $idx >= 1 ; $idx--) {
#print "i2=$idx\n";
# 一つ上の界面でのFresnel係数
		my ($rs12, $rp12, $ts12, $tp12, $Rs12, $Rp12, $Ts12, $Tp12)
				= $this->CalInterfaceFresnelCoefficient($E, $idx-1, $idx);

		my ($esTop, $epTop) = $this->CalEta($E, $idx-1);

#干渉性
# 真ん中(idx)での片道位相差
		my ($b2, $B2) = $this->PhaseShift($E, $idx);
		my $expb2  = exp(-$b2 * i);
		my $exp2b2 = $expb2 * $expb2;
		my $ks = 1.0 / (1.0 + $rs12 * $rs23 * $exp2b2);
		my $kp = 1.0 / (1.0 + $rp12 * $rp23 * $exp2b2);
		my $rs123 = $ks * ($rs12 + $rs23 * $exp2b2);
		my $rp123 = $kp * ($rp12 + $rp23 * $exp2b2);
		my $ts123 = $ks * ($ts12 * $ts23 * $expb2 );
		my $tp123 = $kp * ($tp12 * $tp23 * $expb2 );
#print "$E: [$idx0-$idx1-$idx2] rs123=($rs123)  rp123=($rp123)\n";

#非干渉性
		my $C1 = $this->Coherency($idx); #pLayer($idx)->Coherency();
#my $name = $this->pLayer($idx)->Name();
#print "idx=$idx: c=$C1 name=$name\n";
		if($C1 < 0.999999) {
			my $Cr = 1.0 - $C1;	
			my $ks = Re($esBottom) / Re($esTop);
			my $kp = Re($epBottom) / Re($epTop);
# 振幅の位相変化B2から、エネルギー透過係数Tiをだすため、B2は2倍する
			$B2 = $B2 + $B2;
			my $Ti = abs(exp(-$B2 * i));
			my $Ti2 = $Ti * $Ti;
#print "idx=$idx: Rs23=$Rs23, Rs12=$Rs12\n"; # Rp23=$Rp23, Rp12=$Rp12\n";
			my $Ks = 1.0 / (1.0 - $Rs12 * $Rs23 * $Ti2);
			my $Kp = 1.0 / (1.0 - $Rp12 * $Rp23 * $Ti2);
			my $Rs123 = $Rs12 + $Ks * $Ts12 * $Rs23 * $Ts12 * $Ti2;
			my $Rp123 = $Rp12 + $Kp * $Tp12 * $Rp23 * $Tp12 * $Ti2;
#			my $Rs123 = $Rs12 + $Ks * $Ts12 * $Rs12 * $Ts12 * $Ti2;
#			my $Rp123 = $Rp12 + $Kp * $Tp12 * $Rp12 * $Tp12 * $Ti2;
			my $Ts123 = $Ks * $Ts12 * $Ts23 * $Ti;
			my $Tp123 = $Kp * $Tp12 * $Tp23 * $Ti;
#print "idx=$idx: Rs123+Ts123=", $Ts123+$Rs123, ", Rs12+Ts12=", $Rs12+$Ts12, ", Rs23+Ts23=", $Rs23+$Ts23, "\n";
			my $rRs123 = sqrt($Rs123);
			my $rRp123 = sqrt($Rp123);
#print "idx=$idx: ks=$ks (", Re($esBottom), ",", Re($esTop), ")\n";
			my $tTs123 = sqrt($Ts123 / $ks);
			my $tTp123 = sqrt($Tp123 / $kp);

			$rs123 = $C1 * $rs123 + $Cr * $rRs123;
			$rp123 = $C1 * $rp123 + $Cr * $rRp123;
			$ts123 = $C1 * $ts123 + $Cr * $tTs123;
			$tp123 = $C1 * $tp123 + $Cr * $tTp123;
		}

		$rp23 = $rp123;
		$rs23 = $rs123;
		$tp23 = $tp123;
		$ts23 = $ts123;

#ここまでの界面の反射率、透過率を、次の計算の反射率、透過率として保存
		$Rs23 = abs($rs23);
		$Rs23 = $Rs23 * $Rs23;
		$Rp23 = abs($rp23);
		$Rp23 = $Rp23 * $Rp23;
		$Ts23 = abs($ts23);
		$Ts23 = Re($esBottom) / Re($esTop) * $Ts23 * $Ts23;
		$Tp23 = abs($tp23);
		$Tp23 = Re($epBottom) / Re($epTop) * $Tp23 * $Tp23;
	}

	my ($esTop, $epTop) = $this->CalEta($E, 0);
	my $Rs = abs($rs23);
	$Rs = $Rs * $Rs;
	my $Ts = abs($ts23);
	$Ts = Re($esBottom) / Re($esTop) * $Ts * $Ts;
	my $Rp = abs($rp23);
	$Rp = $Rp * $Rp;
	my $Tp = abs($tp23);
	$Tp = Re($epBottom) / Re($epTop) * $Tp * $Tp;

	return ($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp);
}

#特性行列を使って計算
sub CalFresnelCoefficient_ByCharacteristicMatrix
{
	my ($this, $E) = @_;
	my $nLayer = $this->nLayer();
	if($nLayer <= 2) {
		my ($rs, $ts, $rp, $tp) = $this->CalFresnelCoefficientForBilayer($E);
		return ($rs, $rp, $rp, $tp);
	}

	my $wl = Optics::eVTonm($E) * 1.0e-9; # m
	my %Mstotal = Ellipsometry::CalCharacteristicMatrix(undef, 'unit');
	my %Mptotal = Ellipsometry::CalCharacteristicMatrix(undef, 'unit');

# 光の入射は $idx1=$idx0-1  =>  $idx0 の方向ｓ
	for(my $i = 0 ; $i < $nLayer ; $i++) {
		my $idx0        = $i;
		my $pLayer0     = $this->pLayer($idx0);
		my $thickness0  = $this->LayerThickness($idx0);
		my ($e10, $e20) = $pLayer0->CalEps($E);
		my ($n0, $k0)   = Optics::EpsToNK($e10, $e20);
		my $Angle0      = $this->LightAngle($idx0, $E);
#print "Angle0=$Angle0\n";
		my %Ms = Ellipsometry::CalCharacteristicMatrix(undef, 's', $E, $n0, $k0, $thickness0, $Angle0);
		my %Mp = Ellipsometry::CalCharacteristicMatrix(undef, 'p', $E, $n0, $k0, $thickness0, $Angle0);
		%Mstotal = Ellipsometry::MultMatrix(undef, \%Mstotal, \%Ms);
		%Mptotal = Ellipsometry::MultMatrix(undef, \%Mptotal, \%Mp);
	}

	my $pLayer0         = $this->pLayer(0);
	my ($e10, $e20)     = $pLayer0->CalEps($E);
	my ($n0, $k0)       = Optics::EpsToNK($e10, $e20);
	my $pLayersub       = $this->pLayer($nLayer-1);
	my ($e1sub, $e2sub) = $pLayersub->CalEps($E);
	my ($nsub, $ksub)   = Optics::EpsToNK($e1sub, $e2sub);
	my $Angle0      = $this->LightAngle(0,         $E);
	my $Anglesub    = $this->LightAngle($nLayer-1, $E);

	my ($rs, $ts, $Rs, $Ts) = Ellipsometry::CharacteristicMatrixToFresnelCoefficient(undef, 's', \%Mstotal, $n0, $k0, $nsub, $ksub, $Angle0, $Anglesub);
	my ($rp, $tp, $Rp, $Tp) = Ellipsometry::CharacteristicMatrixToFresnelCoefficient(undef, 'p', \%Mptotal, $n0, $k0, $nsub, $ksub, $Angle0, $Anglesub);
#print "Last: rs123=($rs)  rp123=($rp)\n";
	return ($rs, $rp, $ts, $tp, $Rs, $Rp, $Ts, $Tp);
}

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

	my $Angle = $this->IncidentAngle();
	my $el = new Ellipsometry();
	$el->SetIncidentAngle($Angle);

	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(pseudo),e2(pseudo)");

	my $nLayer = $this->nLayer();
	for(my $i = 0 ; $i < $nLayer ; $i++) {
		my $pLayer = $this->pLayer($i);
		my $name   = $pLayer->Name();
		$out->print(",e1($name),e2($name)");
	}
	$out->print("\n");

	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Emin + $i * $EStep;
		my ($Psi, $Delta) = $this->CalPsiDelta($E);
		my ($e1c, $e2c)   = $el->PsiDeltaToEps($Psi, $Delta, $Angle);
		$out->print("$E");
		if($pObs1) {
			$out->print(",$pObs1->[$i],$pObs2->[$i]");
		}
		$out->print(",$e1c,$e2c");
		for(my $i = 0 ; $i < $nLayer ; $i++) {
			my $pLayer = $this->pLayer($i);
			my ($e1, $e2) = $pLayer->CalEps($E);
			$out->print(",$e1,$e2");
		}
		$out->print("\n");
	}
	$out->Close();
	return 1;
}

sub PrintLayers
{
	my ($this) = @_;

	print("Layers:\n");
	my $nLayer = $this->nLayer();
	for(my $i = 0 ; $i < $nLayer ; $i++) {
		my $pLayer = $this->pLayer($i);
		my $name   = $pLayer->Name();
		my $d      = $this->LayerThickness($i) * 1.0e9; # nm
		my $c      = $this->Coherency($i);
		if($d < 0.0) {
			print("$i: $name (infinite)\n");
		}
		elsif($d < 500.0) {
			print("$i: $name ($d nm) [Coherency: $c]\n");
		}
		elsif($d < 500.0e3) {
			$d *= 1.0e-3;
			print("$i: $name ($d um) [Coherency: $c]\n");
		}
		elsif($d < 500.0e6) {
			$d *= 1.0e-6;
			print("$i: $name ($d mm) [Coherency: $c]\n");
		}
		elsif($d < 500.0e9) {
			$d *= 1.0e-9;
			print("$i: $name ($d m) [Coherency: $c]\n");
		}
	}
}

#=================================================
# 代替があるので不要なサブルーチン
#=================================================
sub CalEpsRTForSubstrateAtNormalncidence
{
	my ($this, $E) = @_;
	my ($Rtots, $Rtotp, $Ttots, $Ttotp) = $this->CalEpsRTForSubstrate($E, 0.0);
	return ($Rtots, $Rtotp);
}

sub CalFresnelCoefficientForBilayer
{
	my ($this, $E, $AngleIncident) = @_;
	$AngleIncident = $this->IncidentAngle() if(!defined $AngleIncident);
	my $nLayer = $this->nLayer();
	my $wl = Optics::eVTonm($E) * 1.0e-9; # m
	my ($ninc, $kinc);
	my ($nsub, $ksub);
	if($nLayer == 1) {
		($ninc, $kinc) = (1.0, 0.0);
		($nsub, $ksub) = $this->pLayer(0)->CalNK($E);
	}
	elsif($nLayer == 2) {
		($ninc, $kinc) = $this->pLayer(0)->CalNK($E);
		($nsub, $ksub) = $this->pLayer(1)->CalNK($E);
	}
	else {
print "Error in Multilayer::CalFresnelCoefficientForBilayer: Invalid nLayer [$nLayer].\n";
		return undef;
	}
	my ($rs, $ts, $rp, $tp) 
		= Ellipsometry::CalFresnelCoefficients(undef, $ninc, $kinc, $nsub, $ksub, $AngleIncident);
	return ($rs, $ts, $rp, $tp);
}

sub CalEpsRTForSubstrate
{
	my ($this, $E, $Angle) = @_;
	my $n = $this->nLayer();
	if($n != 3) {
		print "MultiLayer::CalEpsRTForSubstrateAtNormalncidence: "
			 ."nLayer [$n] must be 3 including top and bottom air layers.\n";
		return undef;
	}

	my ($rho1s, $rho1p, $tau1s, $tau1p, $Rcal1s, $Rcal1p, $Tcal1s, $Tcal1p)
				 = $this->CalInterfaceFresnelCoefficient($E, 0, 1);
	my ($rho2s, $rho2p, $tau2s, $tau2p, $Rcal2s, $Rcal2p, $Tcal2s, $Tcal2p)
				 = $this->CalInterfaceFresnelCoefficient($E, 0, 1);
	my ($PSFilm, $PSSub) = $this->PhaseShift($E, 1);
	my $Tint  = abs(exp(-2.0*$PSSub * i));
	my $Tint2 = $Tint * $Tint;
	my $Rtots = $Rcal1s + $Tcal1s*$Tcal1s*$Rcal2s*$Tint2 / (1.0 - $Rcal1s*$Rcal2s*$Tint2);
	my $Rtotp = $Rcal1p + $Tcal1p*$Tcal1p*$Rcal2p*$Tint2 / (1.0 - $Rcal1p*$Rcal2p*$Tint2);
	my $Ttots = $Tcal1s*$Tcal2s*$Tint / (1.0 - $Rcal1s*$Rcal2s*$Tint2);
	my $Ttotp = $Tcal1p*$Tcal2p*$Tint / (1.0 - $Rcal1p*$Rcal2p*$Tint2);

	return ($Rtots, $Rtotp, $Ttots, $Ttotp);
}

#仮想界面を使って計算(rpの符号が反転する => Δの定義などがひっくり返っている?)
#上部から計算し、下部反射率を正確に計算していないので間違えている。
sub CalReflectionFresnelCoefficient_FromTop
{
	my ($this, $E) = @_;
	my $wl = Optics::eVTonm($E) * 1.0e-9; # m
	my $nLayer = $this->nLayer();
	if($nLayer <= 1) {
print "Note: MultiLayer::CalReflectionFresnelCoefficient: "
	 ."Use Ellipsometry::CalFresnelCoefficients for nLayer=$nLayer\n";
		my ($nsub, $ksub) = $this->pLayer(0)->CalNK($E);
		my ($rs, $ts, $rp, $tp) 
			= Ellipsometry::CalFresnelCoefficients(undef, 1.0, 0.0, $nsub, $ksub, $this->IncidentAngle());
#		my ($rs, $ts, $rp, $tp) = $this->CalFresnelCoefficientForBilayer($E);
		return ($rs, $rp, $rp, $tp);
	}

# 光の入射は $idx=0 => $nLayer-1 の方向ｓ
	my ($rs01, $rp01, $ts01, $tp01, $Rs01, $Rp01, $Ts01, $Tp01)
			 = $this->CalInterfaceFresnelCoefficient($E, 0, 1);
	my ($esTop, $epTop) = $this->CalEta($E, 0);

# idx: 考えている、光が透過する層
	for(my $idx = 1 ; $idx <= $nLayer-2 ; $idx++) {
#print "i1=$idx\n";
# 一つ下の層でのFresnel係数
		my ($rs12, $rp12, $ts12, $tp12, $Rs12, $Rp12, $Ts12, $Tp12)
					= $this->CalInterfaceFresnelCoefficient($E, $idx, $idx+1);

		my ($esBottom, $epBottom) = $this->CalEta($E, $idx+1);

#干渉性
		my ($b2, $B2) = $this->PhaseShift($E, $idx);
		my $expb2     = exp(-$b2 * i);
		my $exp2b2    = $expb2 * $expb2;
		my $ks = 1.0 / (1.0 + $rs01 * $rs12 * $exp2b2);
		my $kp = 1.0 / (1.0 + $rp01 * $rp12 * $exp2b2);
		my $rs012 = $ks * ($rs01 + $rs12 * $exp2b2);
		my $rp012 = $kp * ($rp01 + $rp12 * $exp2b2);
		my $ts012 = $ks * ($ts01 * $ts12 * $expb2 );
		my $tp012 = $kp * ($tp01 * $tp12 * $expb2 );

#非干渉性
		my $C1 = $this->Coherency($idx); #pLayer($idx)->Coherency();
		if($C1 < 0.999999) {
			my $Cr = 1.0 - $C1;	
			my $ks = Re($esBottom) / Re($esTop);
			my $kp = Re($epBottom) / Re($epTop);
# 振幅の位相変化B2から、エネルギー透過係数Tiをだすため、B2は2倍する
			$B2 = $B2 + $B2;
			my $Ti = abs(exp(-$B2 * i));
			my $Ti2 = $Ti * $Ti;
			my $Ks = 1.0 / (1.0 - $Rs01 * $Rs12 * $Ti2);
			my $Kp = 1.0 / (1.0 - $Rp01 * $Rp12 * $Ti2);
			my $Rs012 = $Rs01 + $Ks * $Ts01 * $Rs12 * $Ts01 * $Ti2;
			my $Rp012 = $Rp01 + $Kp * $Tp01 * $Rp12 * $Tp01 * $Ti2;
#			my $Rs012 = $Rs01 + $Ks * $Ts01 * $Rs01 * $Ts01 * $Ti2;
#			my $Rp012 = $Rp01 + $Kp * $Tp01 * $Rp01 * $Tp01 * $Ti2;
			my $Ts012 = $Ks * $Ts01 * $Ts12 * $Ti;
			my $Tp012 = $Kp * $Tp01 * $Tp12 * $Ti;
			my $rRs012 = sqrt($Rs012);
			my $rRp012 = sqrt($Rp012);
			my $tTs012 = sqrt($Ts012 / $ks);
			my $tTp012 = sqrt($Tp012 / $kp);

			$rs012 = $C1 * $rs012 + $Cr * $rRs012;
			$rp012 = $C1 * $rp012 + $Cr * $rRp012;
			$ts012 = $C1 * $ts012 + $Cr * $tTs012;
			$tp012 = $C1 * $tp012 + $Cr * $tTp012;
		}

#ここまでの界面のフレネル係数を、次の計算のフレネル係数として保存
		$rs01 = $rs012;
		$rp01 = $rp012;
		$ts01 = $ts012;
		$tp01 = $tp012;

#ここまでの界面の反射率、透過率を、次の計算の反射率、透過率として保存
		$Rs01 = abs($rs01);
		$Rs01 = $Rs01 * $Rs01;
		$Rp01 = abs($rp01);
		$Rp01 = $Rp01 * $Rp01;
		$Ts01 = abs($ts01);
		$Ts01 = Re($esBottom) / Re($esTop) * $Ts01 * $Ts01;
		$Tp01 = abs($tp01);
		$Tp01 = Re($epBottom) / Re($epTop) * $Tp01 * $Tp01;
	}

	my ($esBottom, $epBottom) = $this->CalEta($E, $nLayer-1);
	my $Rs = abs($rs01);
	$Rs = $Rs * $Rs;
	my $Ts = abs($ts01);
	$Ts = Re($esBottom) / Re($esTop) * $Ts * $Ts;
	my $Rp = abs($rp01);
	$Rp = $Rp * $Rp;
	my $Tp = abs($tp01);
	$Tp = Re($epBottom) / Re($epTop) * $Tp * $Tp;

	return ($rs01, $rp01, $ts01, $tp01, $Rs, $Rp, $Ts, $Tp);
}
1;