#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';

use strict;
use Math::Complex;

use Sci::Science;

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();

my $e1inf = 1.15;
my $A     = 122;
my $Eg    = 1.20; # eV
my $En0   = 3.45; # eV
my $C     = 2.54; # eV

my $OutputFile = "Optics.csv";
my $EStart = 1.0; # V
my $EEnd   = 5.0; # V
my $nE     = 301;
my $EStep  = ($EEnd-$EStart) / ($nE-1);

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

my $alpha    = sqrt(4.0 * $En02 - $C2);
my $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 );

open(OUT,">$OutputFile") or die "$!\n";
print OUT "hv(eV),e1(Tauc),e2(Tauc),e1(Tauc-Lorentz),e2(Tauc-Lorentz),n,k\n";
for(my $i = 0 ; $i < $nE ; $i++) {
	my $E  = $EStart + $i * $EStep;
	my $E2 = $E*$E;

	my $e2Tauc = 0.0;
	my $e2TL   = 0.0;

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

	my $e1Tauc = $e1inf;

#E == Egの場合、log(|E-Eg|)の2つの項が打ち消しあうので、|E-Eg|を1と置き換えても結果は変わらない
	my $EminusEg = 1.0;
	if($E != $Eg) {
		$EminusEg = abs($E - $Eg);
	}
	my $EplusEg  = $E + $Eg;
	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 ) );

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

	my $ec = cplx($e1TL, -$e2TL);
	my $nc = sqrt($ec);
	my $n = Re($nc);
	my $k = -Im($nc);
	print OUT "$E,$e1Tauc,$e2Tauc,$e1TL,$e2TL,$n,$k\n";
}

close(OUT);
exit;