#===============================================
# Optics
#===============================================
package Optics;
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);

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 ODToAlpha { my($this,$OD,$thickness)=@_; return $this->AbsToAlpha($OD, $thickness); }
sub ABSToAlpha
{
	my ($ABS, $thickness) = @_; # m
	return log(10.0) * $ABS / ($thickness * 1.0e-2);
}

sub EpsToConductivity
{
	my ($E, $e2) = @_;
	my $omega = 2.0 * $pi * eVToHz($E);
	return $omega * $e2;
}

sub EpsToNK
{
	my ($e1, $e2) = @_;
	my $ec = cplx($e1, -$e2);
	my $nc = sqrt($ec);
	return (Re($nc), -Im($nc));

#	my $a = sqrt($e1*$e1 + $e2*$e2);
#	my $n = sqrt(0.5 * ($a + $e1));
#	my $k = sqrt(0.5 * ($a - $e1));
#	return ($n, $k);
}

sub NKToEps
{
	my ($n, $k) = @_;
	my $e1 = $n*$n - $k*$k;
	my $e2 = 2.0 * $n * $k;
	return ($e1, $e2);
}

sub AlphaToK
{
	my ($alpha, $wavelength) = @_; # nm
	return $alpha / (4.0 * $pi / ($wavelength*1.0e-7));
}

sub KToAlpha
{
	my ($wavelength, $k) = @_; # nm
	return $k * (4.0 * $pi / ($wavelength*1.0e-7));
}

sub eVToKiser
{
	my ($energy) = @_; #eV
	return 0.0 if($energy == 0.0);
#$e = 1.60218e-19;
#$c = 299792458.0;
#$h = 6.62608e-34;
#print"$e $c $h\n";	
	my $wl = $c * $h / ($energy * $e); # m
	return 1.0 / ($wl * 100.0); # cm-1
}

sub nmToeV
{
	my ($wl) = @_; # nm
	return 0.0 if($wl == 0.0);
	return $c * $h / ($wl*1.0e-9 * $e); # eV
}

sub eVTonm
{
	my ($energy) = @_; # energy in eV
	return 0.0 if($energy == 0.0);
	return $c * $h / ($energy * $e) * 1.0e9; # nm
}

sub eVToHz
{
	my ($energy) = @_; # energy in eV
#$h = 6.62608e-34;
#$e = 1.60218e-19;
#print "$energy * $e / $h\n";
	return $energy * $e / $h; # Hz
}

sub HzToeV
{
	my ($f) = @_;
	return $f * $h / $e;
}

sub eVToOmega
{
	my ($energy) = @_; # energy in eV
	return $pi2 * $energy * $e / $h; # Hz
}

sub OmegaToeV
{
	my ($f) = @_;
	return $f * $h / $e / $pi2;
}

sub HzTonm
{
	my ($f) = @_;
	return eVTonm(HzToeV($f));
}

sub nmToHz
{
	my ($wl) = @_; # nm
	return $c / ($wl*1.0e-9);
}

sub CalEffectiveMassFromPlasmaEnergy
{
	my ($Ep, $Ne) = @_;
	my $wp = 2.0 * $pi * Optics::eVToHz($Ep);
	my $meeff = $e * $e * $Ne / $e0 / $me / $wp / $wp;
	return $meeff;
}

sub CalPlasmaFrequency
{
# Ne: m-3
	my ($Ne, $m) = @_;
#print "e=$e  me=$me  e0=$e0  pi=$pi\n";
#print "Ne=$Ne  m=$me\n";
	return sqrt( $Ne *$e*$e / $me / $m / $e0) / 2.0 / $pi;
}


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

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub KKConversionFromEps2toEps1
{
	my ($this, $pE, $pEps2, $fconst, $ExtraporateLow, $ExtraporateHigh, $Integrator) = @_;

	my @ec;
	my @func;
	for(my $i = 0 ; $i < @$pE ; $i++) {
		my $E = $pE->[$i];
		$func[$i] = $E * $pEps2->[$i];
	}

	return $this->KKConversionGeneral($pE, \@func, $fconst, 
					$ExtraporateLow, $ExtraporateHigh, $Integrator);
}

#$pfには、KK変換(積分範囲０〜∞)の被積分関数の分子を与える
#つまり、ε2(E)のKK変換をする場合、Eε2(E)の配列のレファレンスを与える
sub KKConversionGeneral
{
	my ($this, $pE, $pf, $fconst, $ExtraporateLow, $ExtraporateHigh, $Integrator) = @_;

	my $nData = @$pE;
	my $MinE = $pE->[0];
	my $dE = $pE->[1] - $pE->[0];
	my @Q;
	my $Coeff = 2.0/$pi * $dE * 2.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $Eg = $pE->[$i];
		my $yg = $pf->[$i];
		my $idx = int(($Eg - $MinE) / $dE + 0.001);
		my $j0 = 0;
		$j0 = 1  if($idx % 2 == 0);

		my $v = 0.0;
		if($Integrator eq "Rectangular") {
			for(my $j = $j0 ; $j < $nData ; $j += 2) {
				my $Ej = $pE->[$j];
				$v += ($pf->[$j] - $yg) / ($Ej*$Ej - $Eg*$Eg);
			}
		}
		elsif($Integrator eq "Simpson") {
			for(my $j = $j0 ; $j <= $nData-5 ; $j += 2) {
				next if($i-3 <= $j and $j <= $i);
				my $Ej  = $pE->[$j];
				my $Ej2 = $pE->[$j+2];
				my $Ej4 = $pE->[$j+4];
				$v += 1.0/3.0 * (
					    0.5 * ($pf->[$j]   - $yg) / ($Ej*$Ej   - $Eg*$Eg)
					  + 2.0 * ($pf->[$j+2] - $yg) / ($Ej2*$Ej2 - $Eg*$Eg)
					  + 0.5 * ($pf->[$j+4] - $yg) / ($Ej4*$Ej4 - $Eg*$Eg)
					);
			}
		}
		else {
			$this->print("Error!!: Invalid Integrator [$Integrator]\n");
			return -5;
		}
#		$v *= $Coeff * $Eg;
		$v *= $Coeff;

		if($ExtraporateLow  eq "Drude") {
		}
		elsif($ExtraporateLow  eq "None") {
		}
		else {
			$this->print("Error!!: Invalid Extraporation to low f [$ExtraporateLow]\n");
			return -3;
		}
		if($ExtraporateHigh eq "Constant") {
			my $jmax = $nData - 1;
			$jmax-- if(($j0+$jmax) % 2 != 0);
			$j0 = 0 if($j0 == 1);
			my $Emax    = $pE->[$jmax] + $dE;
			$Emax += $dE if(abs($Emax-$Eg) < 1.0);
			my $ymax = $pf->[$jmax] - $yg;
			$v += $ymax / $pi * log( ($Emax+$Eg)/($Emax-$Eg) );
		}
		elsif($ExtraporateHigh eq "e2~1/E^3") {
			my $jmax = $nData - 1;
			$jmax-- if(($j0+$jmax) % 2 != 0);
			$j0 = 0 if($j0 == 1);
			my $Emax    = $pE->[$jmax] + $dE;
			$Emax += $dE if(abs($Emax-$Eg) < 1.0);
			my $ymax = $pf->[$jmax] - $yg;
			my $K = $ymax * $Emax * $Emax * 2.0/$pi;
			$v += $K / $Eg / $Eg * ( -1.0/$Emax + 0.5/$Eg * log( ($Emax+$Eg) / abs($Emax-$Eg) ) );
		}
		elsif($ExtraporateHigh eq "e2~1/E^4") {
			my $jmax = $nData - 1;
			$jmax-- if(($j0+$jmax) % 2 != 0);
			$j0 = 0 if($j0 == 1);
			my $Emax    = $pE->[$jmax] + $dE;
			$Emax += $dE if(abs($Emax-$Eg) < 1.0);
			my $ymax = $pf->[$jmax] - $yg;
			my $K = $ymax * $Emax * $Emax * $Emax * 2.0/$pi;
			$v += $K / 2.0 / $Eg / $Eg * ( 1.0/$Eg/$Eg * log( $Emax*$Emax / abs($Emax*$Emax-$Eg*$Eg) ) - 1.0 / $Emax/$Emax);
		}
		elsif($ExtraporateHigh  eq "None") {
		}
		else {
			$this->print("Error!!: Invalid Extraporation to high f [$ExtraporateHigh]\n");
			return -3;
		}

		$Q[$i] = $fconst + $v;
	}

	return (\@Q);
}

sub KKConversionFromEps1toEps2
{
	my ($this, $pE, $pEps2, $fconst, $ExtraporateLow, $ExtraporateHigh, $Integrator) = @_;

	my @ec;
	my @func;
	for(my $i = 0 ; $i < @$pE ; $i++) {
		my $E = $pE->[$i];
		$func[$i] = $pEps2->[$i] - $fconst;
	}

	return $this->KKReverseConversionGeneral($pE, \@func, $fconst, 
					$ExtraporateLow, $ExtraporateHigh, $Integrator);
}


#$pfには、KK変換(積分範囲０〜∞)の被積分関数の分子を与える
#つまり、ε2(E)のKK変換をする場合、Eε2(E)の配列のレファレンスを与える
sub KKReverseConversionGeneral
{
	my ($this, $pE, $pf, $fconst, $ExtraporateLow, $ExtraporateHigh, $Integrator) = @_;

	my $nData = @$pE;
	my $MinE = $pE->[0];
	my $dE = $pE->[1] - $pE->[0];
	my @Q;
	my $Coeff = 2.0/$pi * $dE * 2.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $Eg = $pE->[$i];
		my $yg = $pf->[$i];
		my $idx = int(($Eg - $MinE) / $dE + 0.001);
		my $j0 = 0;
		$j0 = 1  if($idx % 2 == 0);

		my $v = 0.0;
		if($Integrator eq "Rectangular") {
			for(my $j = $j0 ; $j < $nData ; $j += 2) {
				my $Ej = $pE->[$j];
				$v += ($pf->[$j] - $yg) / ($Ej*$Ej - $Eg*$Eg);
			}
		}
		elsif($Integrator eq "Simpson") {
			for(my $j = $j0 ; $j <= $nData-5 ; $j += 2) {
				next if($i-3 <= $j and $j <= $i);
				my $Ej  = $pE->[$j];
				my $Ej2 = $pE->[$j+2];
				my $Ej4 = $pE->[$j+4];
				$v += 1.0/3.0 * (
					    0.5 * ($pf->[$j]   - $yg) / ($Ej*$Ej   - $Eg*$Eg)
					  + 2.0 * ($pf->[$j+2] - $yg) / ($Ej2*$Ej2 - $Eg*$Eg)
					  + 0.5 * ($pf->[$j+4] - $yg) / ($Ej4*$Ej4 - $Eg*$Eg)
					);
			}
		}
		else {
			$this->print("Error!!: Invalid Integrator [$Integrator]\n");
			return -5;
		}
		$v *= -$Coeff * $Eg;

		if($ExtraporateLow  eq "Drude") {
		}
		elsif($ExtraporateLow  eq "None") {
		}
		else {
			$this->print("Error!!: Invalid Extraporation to low f [$ExtraporateLow]\n");
			return -3;
		}
		if($ExtraporateHigh eq "Constant") {
			my $jmax = $nData - 1;
			$jmax-- if(($j0+$jmax) % 2 != 0);
			$j0 = 0 if($j0 == 1);
			my $Emax    = $pE->[$jmax] + $dE;
			$Emax += $dE if(abs($Emax-$Eg) < 1.0);
			my $ymax = $pf->[$jmax] - $yg;
			$v += -$ymax / $pi * log( ($Emax+$Eg)/($Emax-$Eg) );
		}
		elsif($ExtraporateHigh eq "e2~1/E^3") {
			my $jmax = $nData - 1;
			$jmax-- if(($j0+$jmax) % 2 != 0);
			$j0 = 0 if($j0 == 1);
			my $Emax    = $pE->[$jmax] + $dE;
			$Emax += $dE if(abs($Emax-$Eg) < 1.0);
			my $ymax = $pf->[$jmax] - $yg;
			my $K = $ymax * $Emax * $Emax * 2.0/$pi;
			$v += -$K / $Eg * ( -1.0/$Emax + 0.5/$Eg * log( ($Emax+$Eg) / abs($Emax-$Eg) ) );
		}
		elsif($ExtraporateHigh eq "e2~1/E^4") {
			my $jmax = $nData - 1;
			$jmax-- if(($j0+$jmax) % 2 != 0);
			$j0 = 0 if($j0 == 1);
			my $Emax    = $pE->[$jmax] + $dE;
			$Emax += $dE if(abs($Emax-$Eg) < 1.0);
			my $ymax = $pf->[$jmax] - $yg;
			my $K = $ymax * $Emax * $Emax * $Emax * 2.0/$pi;
			$v += -$K / 2.0 / $Eg * ( 1.0/$Eg/$Eg * log( $Emax*$Emax / abs($Emax*$Emax-$Eg*$Eg) ) - 1.0 / $Emax/$Emax);
		}
		elsif($ExtraporateHigh  eq "None") {
		}
		else {
			$this->print("Error!!: Invalid Extraporation to high f [$ExtraporateHigh]\n");
			return -3;
		}

		$Q[$i] = $v;
	}

	return (\@Q);
}

no Math::Complex;

sub TaucLorentz
{
	my ($this, $E, $e1inf, $e2inf, $A, $Eg, $En0, $C) = @_;
#print "E=$E, e=($e1inf,$e2inf) A=$A Eg=$Eg En0=$En0 C=$C\n";

	return ($e1inf, $e2inf) if($A <= 0.0 or $C <= 0.0);

my $Eg2  = $Eg*$Eg;
my $En02 = $En0*$En0;
my $C2   = $C*$C;

my $alpha = 0.1;
my $gamma = 0.1;
if(2.0*$En02 > $C2) {
	$alpha = sqrt(4.0 * $En02 - $C2);
	$gamma = sqrt($En02 - $C2/2.0);
}
my $alpha2   = $alpha*$alpha;
my $gamma2   = $gamma*$gamma;

my $lnEn0g   = log( ($En02+$Eg2+$alpha*$Eg) / ($En02+$Eg2-$alpha*$Eg) );
my $atanEgaC = $pi - Sci::atan( (2.0*$Eg+$alpha)/$C ) + Sci::atan( (-2.0*$Eg+$alpha)/$C );
my $atangEg  = $pi + 2.0*Sci::atan( 2.0 * ($gamma2-$Eg2) / $alpha / $C );

	my $E2 = $E*$E;

	my $e2TL = $e2inf;
	if($E > $Eg) {
#print "E=$E($Eg)\n";
		my $e2Tauc = $A * ($E - $Eg)**2 / $E / $E;
		my $L = $En0 * $C * $E / ( ($E2-$En02)*($E2-$En02) + $C2*$E2 );
		$e2TL = $e2inf + $e2Tauc * $L;
	}

#E == Egの場合、log(|E-Eg|)の2つの項が打ち消しあうので、|E-Eg|を1と置き換えても結果は変わらない
	my $EminusEg = 1.0;
	if($E != $Eg) {
		$EminusEg = abs($E - $Eg);
	}
	my $EplusEg  = $E + $Eg;
#print "Eg=$Eg, $EminusEg, $EplusEg\n";# if($EminusEg < 0);
	my $aln      = ($Eg2-$En02) * $E2 + $Eg2*$C2 - $En02 * ($En02 + 3.0*$Eg2);
	my $atan     = ($E2-$En02) * ($En02+$Eg2) + $Eg2*$C2;
	my $gsi4     = ($E2-$gamma2)*($E2-$gamma2) + $alpha2*$C2/4.0;
	my $lnE = log( $EminusEg / $EplusEg );
	my $lnEngrt  = log( $EminusEg * $EplusEg / sqrt( ($En02-$Eg2)*($En02-$Eg2) + $Eg2*$C2 ) );

#print "En0=$En0, alpha=$alpha, gsi4=$gsi4, E=$E\n";
	my $e1TL = $e1inf;
	$e1TL += + ($A*$C/$pi/$gsi4) * ($aln/2.0/$alpha/$En0) * $lnEn0g;
	$e1TL += - ($A/$pi/$gsi4) * ($atan/$En0) * $atanEgaC;
	$e1TL += + 2.0 * ($A*$En0/$pi/$gsi4/$alpha) * $Eg * ($E2-$gamma2) * $atangEg;
	$e1TL += - ($A*$En0*$C/$pi/$gsi4) * ($E2+$Eg2)/$E * $lnE;
	$e1TL += + (2.0*$A*$En0*$C/$pi/$gsi4) * $Eg * $lnEngrt;

	return ($e1TL, $e2TL);
}

use Math::Complex;

sub CalNefromLorentzA
{
	my ($this, $A, $fa) = @_;
	my $K = 2.0 * $pi * $e / $h;
	return $A / ($fa * $e*$e / $me / $e0 / $K / $K);
}

sub CalLorentzAfromNe
{
	my ($this, $Ne, $fa) = @_;
	my $K = 2.0 * $pi * $e / $h;
	return $fa * $Ne * $e*$e / $me / $e0 / $K / $K;
}

sub Lorentz
{
# Ne: m-3
	my ($this, $E, $e1inf, $e2inf, $fa, $Ne, $E0, $gamma) = @_;
#print "einf=($e1inf,$e2inf)  fa=$fa  Ne=$Ne  E0=$E0   g=$gamma\n";
	my $K = 2.0 * $pi * $e / $h;
	my $ac = cplx($fa * $Ne * $e*$e / $me / $e0 / $K / $K, 0.0);
#print "ac=", Re($ac), " + i ", -Im($ac), "\n";

	my $ec = $ac / ($E0*$E0 - $E*$E + i * $gamma * $E0);
	my $e1 = $e1inf + Re($ec);
	my $e2 = $e2inf - Im($ec);
#print "Lorentz: E=$E  e=($e1,$e2)\n";
	return ($e1, $e2);
}

sub Drude
{
	my ($this, $E, $e1inf, $e2inf, $Ep, $tau) = @_;
#print "e=$e1inf,$e2inf, Ep=$Ep  tau=$tau\n";
	my $omega = 2.0 * $pi * eVToHz($E);
	my $ot = $omega * $tau;
	my $wpe   = 2.0 * $pi * $Ep * $e / $h;
#print "wp=$wpe ot=$ot  ";
	my $K = $wpe*$wpe*$tau / (1.0 + $ot*$ot);
	my $e1 = $e1inf - $K * $tau;
	my $e2 = $e2inf + $K / $omega;
#print "E=$E  e=($e1,$e2) [$Ep,$tau]\n";
	return ($e1, $e2);
}

# n = sum(Ai / wl^ni)
sub Caucy
{
	my ($this, $wl, $pA, $pn) = @_;  # wl (nm)
	my $n = 0.0;
	for(my $i = 0 ; $i < @$pA ; $i++) {
#print "$i: $pn->[$i]: $pA->[$i], wl=$wl\n";
		if($pn->[$i] == 0) {
			$n += $pA->[$i];
		}
		else {
			$n += $pA->[$i] / ($wl)**($pn->[$i]);
		}
	}
	return ($n*$n, 0.0);
}

# eps = A + sum(Bi*WL^2/(WL^2-WL0i^2)
sub Sellmeier
{
	my ($this, $wl, $A, $pB, $pwl0) = @_;
	my $e1 = $A;
	my $wl2  = $wl*$wl;
	for(my $i = 0 ; $i < @$pB ; $i++) {
		my $wl0 = $pwl0->[$i];
		$e1 += $pB->[$i] * $wl2 / ($wl2 - $wl0*$wl0);
	}
	return ($e1, 0.0);
}

sub Urbach
{
	my ($this, $E, $Em, $E0, $A) = @_;
	my $e2 = 0.0;
	if($E > $Em) {
		$e2 += $A;
	}
	else {
		$e2 += $A * exp(-($Em - $E) / $E0);
	}
	return (0.0, $e2);
}

sub Urbach2
{
	my ($this, $E, $pEm, $pE0, $pA) = @_;
	my $e2 = 0.0;
	for(my $i = 0 ; $i < @$pEm ; $i++) {
		if($E > $pEm->[$i]) {
			$e2 += $pA->[$i];
		}
		else {
			$e2 += $pA->[$i] * exp(-($pEm->[$i] - $E) / $pE0->[$i]);
		}
	}
	return (0.0, $e2);
}