#!/usr/bin/perl

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "d:/Programs/Perl/lib", "c:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;
#use Math::Complex;

use CSV;
use JFile;
use MyHTMLApplication;
use Sci qw(nmToeV eVTonm $e);
use MultiColumnData;
use Crystal::VASP;
#use Sci::Color;
#use Sci::SolarCell;
use Sci::Science;
use Sci::Optics;
use Sci::OpticalMaterial;
use Sci::MultiLayer;
#use Sci::Optimize;
#use Optimize;

#===============================================
# Debug variables
#===============================================
my $PrintLevel = 0;

#===============================================
# Global variables
#===============================================
my $DataPath = "D:\\Programs\\Perl\\Color\\data";
   $DataPath = "C:\\Programs\\Perl\\Color\\data" if(!-d $DataPath);
   $DataPath = Deps::MakePath($ENV{'PerlDir'}, ["Color", "data"], 0) if(!-d $DataPath);
#my $StandardLightPath         = Deps::MakePath($DataPath, "CIEStandardLight.txt", 0);
#my $DayLightPath              = Deps::MakePath($DataPath, "DayLightS.txt", 0);
#my $HalogenLampPath           = Deps::MakePath($DataPath, "HalogenLamp.lsc", 0);
my $SolarSpectraPath          = "SolarSpectra.txt";
   $SolarSpectraPath          = Deps::MakePath($DataPath, "SolarSpectra.txt", 0) if(!-f $SolarSpectraPath);

my ($dModelE, $ModelE1)      = (0.02, 6.0); # eV
my $SubstrateThickness       = 0.5e-3; # m
my $Substrateer              = 1.5;
my $Substrateei              = 0.0;
my $TopElectrodeThickness    = 100.0; # nm
my $TopElectrodeer           = 4.0;
my $TopElectrodeei           = 0.0;
my $BottomElectrodeThickness = 100.0; # nm
my $BottomElectrodeer        = 4.0;
my $BottomElectrodeei        = 0.0;

my $OutLightSourceCSV     = "NormalizedLightSource.csv";
my $OutOpticalSpectraCSV  = "OpticalSpectra.csv";
my $OutFilmSpectraCSV     = "FilmDielectricFunction.csv";
my $OutMultilayerCSV      = "MultilayerSpectra.csv";

#===============================================
# Character code variables
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

#===============================================
# Create Application object
#===============================================
my $App = new MyHTMLApplication;
my $ret = $App->Initialize();
exit(-1) if($ret < 0);

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);
$App->SetOutputMode("console");

#==========================================
# Read command-line arguments
#==========================================
$App->AddArgument("--Eg",        "--Eg=val (Def:)" ,                 undef);
$App->AddArgument("--Voc",       "--Voc=val (Def:)",                 undef);
$App->AddArgument("--VocRatio",  "--VocRatio=val (Def:0.8)",           0.7);
$App->AddArgument("--FF",        "--FF=val (Def: 0.8)",                0.8);
$App->AddArgument("--Power",     "--Power=val (Def: 1000.0 W/m2)" , 1000.0);
$App->AddArgument("--Thickness", "--Thickness=val (Def: 100.0 nm)" , 100.0);
exit 1 if($App->ReadArgs(1, "sjis", 0) != 1);

my $Args = $App->Args();
my %ArgHash = $Args->GetArgHash();
#foreach my $key (keys %ArgHash) {
#print "key: $key: $ArgHash{$key}\n";
#	if($key =~ /^Param:(.*?)$/i) {
#		$ParamHash{$1} = $ArgHash{$key};
#print "$1:  $ArgHash{$key}\n";
#	}
#}

my $RootDir          = $Args->GetGetArg(0);
   $RootDir          = '.' if($RootDir eq '');
my $Eg               = $Args->GetGetArg("Eg",        1, 1);
my $Voc              = $Args->GetGetArg("Voc",       1, 1);
my $VocRatio         = $Args->GetGetArg("VocRatio",  1, 1);
#   $VocRatio         = 0.8 if($VocRatio eq '');
my $FillFactor       = $Args->GetGetArg("FF",        1, 1);
#   $FillFactor       = 0.8 if($FillFactor eq '');
my $LightSourcePower = $Args->GetGetArg("Power",     1, 1);
#   $LightSourcePower = 1000.0 if($LightSourcePower eq '');
my $Thickness        = $Args->GetGetArg("Thickness", 1, 1);

print "Arguments:\n";
print "  Eg              : $Eg eV\n"                 if($Eg ne '');
print "  Voc             : $Voc V\n"                 if($Voc ne '');
print "  VocRatio        : $VocRatio\n"              if($VocRatio ne '');
print "  FillFactor      : $FillFactor\n"            if($FillFactor ne '');
print "  LightSourcePower: $LightSourcePower W/m2\n" if($LightSourcePower ne '');
print "  Thickness       : $Thickness nm\n"          if($Thickness ne '');
print "\n";

#=========================================================
# Main program
#=========================================================

#=============================================================
# Read and save solar spectra
#=============================================================
print "Read solar spectra from [$SolarSpectraPath]\n";
my ($pWL, $pMCDPAll, $pMCDPDirect) = &ReadSolarSpectra($SolarSpectraPath);
my $nPData         = $pMCDPAll->nData();
my $pMCDIntPAll    = $pMCDPAll->IntegrateByTrapezoid();
my $pMCDIntPDirect = $pMCDPDirect->IntegrateByTrapezoid();
my $TotalPAll      = $pMCDIntPAll->YValByIndex($nPData-1);
my $TotalPDirect   = $pMCDIntPDirect->YValByIndex($nPData-1);
printf "  Int(P(All))    = %10.4g W/m2\n", $TotalPAll;
printf "  Int(P(Direct)) = %10.4g W/m2\n", $TotalPDirect;
print  "    Normalize these spectra to $LightSourcePower W/m2...\n";
$pMCDPAll->MultiplyY($LightSourcePower / $TotalPAll);
$pMCDPDirect->MultiplyY($LightSourcePower / $TotalPDirect);
$pMCDIntPAll->MultiplyY($LightSourcePower / $TotalPAll);
$pMCDIntPDirect->MultiplyY($LightSourcePower / $TotalPDirect);
my $pPAll       = $pMCDPAll->GetYData();
my $pPDirect    = $pMCDPDirect->GetYData();
my $pIntPAll    = $pMCDIntPAll->GetYData();
my $pIntPDirect = $pMCDIntPDirect->GetYData();

print  "    Calculate photon flux...\n";
my $pMCDFAll    = $pMCDPAll->Duplicate();
my $pMCDFDirect = $pMCDPDirect->Duplicate();
my $pFAll       = $pMCDFAll->GetYData();
my $pFDirect    = $pMCDFDirect->GetYData();
for(my $i = 0 ; $i < $nPData ; $i++) {
	my $wl          = $pWL->[$i];
	my $E           = Sci::nmToeV($wl);
	$pFAll->[$i]    = $pPAll->[$i] / $E / $e;
	$pFDirect->[$i] = $pPDirect->[$i] / $E / $e;
#print "wl=$wl\tE=$E\t$FAll[$i]\t$FDirect[$i]\n";
}
my $pMCDFDirect    = new MultiColumnData($pWL, $pFDirect);
my $pMCDIntFAll    = $pMCDFAll->IntegrateByTrapezoid();
my $pMCDIntFDirect = $pMCDFDirect->IntegrateByTrapezoid();
my $pIntFAll       = $pMCDIntFAll->GetYData();
my $pIntFDirect    = $pMCDIntFDirect->GetYData();
print "    Save the normlized spectra to [$OutLightSourceCSV]\n";
my $out = JFile->new($OutLightSourceCSV, 'w') or die "$!: Can not write to [$OutLightSourceCSV].\n";
$out->print("WL(nm),E(eV),P(All) (W/m2nm),P(Direct) (W/m2nm),Int(All) (W/m2),Int(Direct) (W/m2),"
		   ."F(All) (/m2snm),F(Direct) (/m2snm),Int(All) (/m2s),Int(Direct) (/m2s)\n");
for(my $i = 0 ; $i < $nPData ; $i++) {
	my $wl = $pWL->[$i];
	my $E  = nmToeV($wl);
	$out->print("$wl,$E,$pPAll->[$i],$pPDirect->[$i],$pIntPAll->[$i],$pIntPDirect->[$i],"
			   ."$pFAll->[$i],$pFDirect->[$i],$pIntFAll->[$i],$pIntFDirect->[$i]\n");
}
$out->Close();
print "\n";

#=============================================================
# Read VASP data if Eg is not given by command-line argument
#=============================================================
my ($EF, $VBM, $CBM);
my ($nOpt, $pOptE, $pOptWL, $pMCDaxx, $pMCDayy, $pMCDazz, 
	$pMCDerxx, $pMCDeixx, $pMCDeryy, $pMCDeiyy, $pMCDerzz, $pMCDeizz);

if($Eg eq '') {
	print "Read VASP data\n";

	my $vasp = new VASP;
	my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($RootDir);
	my $VCRelaxDir = Deps::MakePath($RootDir, "VCRelax", 0);
	my $SCFDir     = Deps::MakePath($RootDir, "SCF", 0);
	my $DOSDir     = Deps::MakePath($RootDir, "DOS", 0);
	my $OpticsPath = Deps::MakePath($DOSDir,  "Optics.csv", 0);

	my $INCARHash          = $vasp->ReadINCARtoHash  (Utils::MakePath($SCFDir,     'INCAR',   '/', 0));
#	my $KPOINTSHash        = $vasp->ReadKPOINTStoHash(Utils::MakePath($SCFDir,     'KPOINTS', '/', 0));
	my $VCRelaxKPOINTSHash = $vasp->ReadKPOINTStoHash(Utils::MakePath($VCRelaxDir, 'KPOINTS', '/', 0));
	my $POTCARHash         = $vasp->ReadPOTCARtoHash (Utils::MakePath($SCFDir,     'POTCAR',  '/', 0));
	my $POSCARHash         = $vasp->ReadPOSCARtoHash (Utils::MakePath($SCFDir,     'POSCAR',  '/', 0));
	my $SCFOUTCARHash      = $vasp->ReadOUTCARtoHash (Utils::MakePath($SCFDir,     'OUTCAR',  '/', 0));
	my $VCROUTCARHash      = $vasp->ReadOUTCARtoHash (Utils::MakePath($VCRelaxDir, 'OUTCAR',  '/', 0));

	my $Status = '';
	if(-d $VCRelaxDir and $VCROUTCARHash->{RequiredAccuracyMessage} =~ /^\s*$/) {
		$Status = 'Relaxation not converged(1)';
	}
	elsif(-d $VCRelaxDir and $VCROUTCARHash->{AbortMessage} !~ /is\s+reached/) {
		$Status = 'Relaxation not converged(2)';
	}
	elsif(-d $VCRelaxDir and $VCROUTCARHash->{TotCPUTime} eq '') {
		$Status = 'VCR aborted(1)';
	}
	elsif($SCFOUTCARHash->{AbortMessage} =~ /^\s*$/) {
		$Status = 'SCF not converged(1)';
	}
	elsif($SCFOUTCARHash->{TotCPUTime} eq '') {
		$Status = 'SCF aborted(1)';
	}
	($EF, $VBM, $CBM, $Eg) = VASP->new()->ReadEgFromDOSOUTCAR(Deps::MakePath($DOSDir, 'OUTCAR', 0));

	$App->print("  ChemicalComposition=$POSCARHash->{ChemicalComposition}\n");
	$App->print("  SumChemicalComposition=$POSCARHash->{SumChemicalComposition}\n");
	$App->print("  FormulaUnit=$POSCARHash->{FormulaUnit}\n");

	$POTCARHash->{PseudoPotential} =~ s/ .*$//s;
	$App->print("  PseudoPotential=$POTCARHash->{PseudoPotential}\n");
	$App->print("  ISPIN=$INCARHash->{ISPIN}\n");
	$App->print("  PREC=$INCARHash->{PREC}\n");
	$App->print("  EDIFF=$INCARHash->{EDIFF}\n");
	$App->print("  EDIFFG=$INCARHash->{EDIFFG}\n");
	$App->print("  LDAU=$INCARHash->{LDAU}\n")         if($INCARHash->{LDAU});
	$App->print("  LDAUTYPE=$INCARHash->{LDAUTYPE}\n") if($INCARHash->{LDAUTYPE});
	$App->print("  LDAUU=$INCARHash->{LDAUU}\n")       if($INCARHash->{LDAUU});
	$App->print("  LDAUJ=$INCARHash->{LDAUJ}\n")       if($INCARHash->{LDAUJ});
	$App->print("  ISMEAR=$INCARHash->{ISMEAR}\n");
	$App->print("  SIGMA=$INCARHash->{SIGMA}\n");
#	$App->print("  KPOINTS=$KPOINTSHash->{KPOINTS}\n");
	print "\n";
	$App->print("  Etot=$SCFOUTCARHash->{TOTEN}\n");
	$App->print("  EF(SCF)=$SCFOUTCARHash->{EF}\n");
	$App->print("  EF(DOS)=$EF\n");
	$App->print("  EVBM(DOS)=$VBM\n");
	$App->print("  ECBM(DOS)=$CBM\n");
	$App->print("  Eg(DOS)=$Eg\n");
	print "\n";

	print "  Read optical spectra from [$OpticsPath]...\n";
	($nOpt, $pOptE, $pOptWL, $pMCDaxx, $pMCDayy, $pMCDazz, $pMCDerxx, $pMCDeixx, $pMCDeryy, $pMCDeiyy, $pMCDerzz, $pMCDeizz) 
			= &ReadOpticalSpectra($OpticsPath);
	print "\n";
}

#================================================
# Maximam absorption case (infinite thickness)
#================================================
my $TotalPhotonFlux = $pIntFDirect->[$nPData-1];
my $nPhotonsEg      = Algorism::Interpolate($pWL, $pIntFDirect, eVTonm($Eg));
my $PFRatio         = $nPhotonsEg / $TotalPhotonFlux * 100.0; # %
my $JscMax          = $nPhotonsEg * $e;
my $Voc             = $Eg  * $VocRatio;
my $EffMax          = $Eg  * $JscMax               / $LightSourcePower * 100.0; # %
my $SupposedEffMax  = $Voc * $JscMax * $FillFactor / $LightSourcePower * 100.0; # %
print "Solar(Direct)\n";
printf "  Total photon flux                  : %10.4g m^-2\n", $TotalPhotonFlux;
printf "  Total photon flux from Eg=%6.4f eV: %10.4g m^-2 (%6.4f %)\n",  $Eg, $nPhotonsEg, $PFRatio;
printf "  Maximum Jsc from Eg=%6.4f eV: %6.4f mA/cm2\n", $Eg, $JscMax * 0.1;
printf "  Maximum Voc : %6.4f V\n", $Eg;
printf "    Maximum Eff =Eg*Jsc    : %6.2f %\n", $EffMax;
printf "  Supposed Voc: %6.4f V\n", $Voc;
printf "  Supposed FF : %6.4f\n",   $FillFactor;
printf "    Supposed Eff=Voc*Jsc*FF: %6.2f %\n", $SupposedEffMax;
#my $ApproxFFByEqCircuit  = ($Voc - log($Voc + 0.72)) / ($Voc + 1.0);
#my $ApproxEffByEqCircuit = $Voc * $JscMax * $ApproxFFByEqCircuit / $LightSourcePower * 100.0; # %
#printf "  Approx. FF(eq. circuit): %6.4f\n",   $ApproxFFByEqCircuit;
#printf "    Approx. Eff=Voc*Jsc*FF(eq. circuit): %6.2f %\n", $ApproxEffByEqCircuit;
print "\n";

#=================================
# Single-path absorption case
#=================================
my $d = $Thickness * 1.0e-7; # cm
if(defined $nOpt) {
#my ($nOpt, $pOptE, $pOptWL, $paxx, $payy, $pazz);
#my $out=JFile->new("a.csv","w") or die $!;
#$out->print("wl,E,axx,ayy,azz\n");
#for(my $i = 0 ; $i < $nOpt ; $i++) {
#	$out->print("$pOptWL->[$i],$pOptE->[$i],$paxx->[$i],$payy->[$i],$pazz->[$i]\n");
#}
#$out->Close();

	printf "Calculated from optical absorption for %6.2f nm single-path tranmission.\n", $Thickness;
	my (@Axx, @Ayy, @Azz, @FDirectAxx, @FDirectAyy, @FDirectAzz);
	for(my $i = 0 ; $i < $nPData ; $i++) {
		my $wl   = $pWL->[$i];
		my $E    = Sci::nmToeV($wl);
		my $axx  = $pMCDaxx->YValLinear($E);
		my $ayy  = $pMCDayy->YValLinear($E);
		my $azz  = $pMCDazz->YValLinear($E);
#print "$i: $wl\t$E\t$axx\n";
		$Axx[$i] = 1.0 - exp(-$axx * $d);
		$Ayy[$i] = 1.0 - exp(-$ayy * $d);
		$Azz[$i] = 1.0 - exp(-$azz * $d);
		my $FD   = Algorism::Interpolate($pWL, $pFDirect, $wl);
		$FDirectAxx[$i] = $pFDirect->[$i] * $Axx[$i];
		$FDirectAyy[$i] = $pFDirect->[$i] * $Ayy[$i];
		$FDirectAzz[$i] = $pFDirect->[$i] * $Azz[$i];
	}
	my $pIntFDirectAxx  = Algorism::IntegrateByTrapezoid($pWL, \@FDirectAxx);
	my $pIntFDirectAyy  = Algorism::IntegrateByTrapezoid($pWL, \@FDirectAyy);
	my $pIntFDirectAzz  = Algorism::IntegrateByTrapezoid($pWL, \@FDirectAzz);
	my $nPhotonsAxx     = Algorism::Interpolate($pWL, $pIntFDirectAxx, eVTonm($Eg));
	my $nPhotonsAyy     = Algorism::Interpolate($pWL, $pIntFDirectAyy, eVTonm($Eg));
	my $nPhotonsAzz     = Algorism::Interpolate($pWL, $pIntFDirectAzz, eVTonm($Eg));
	my $JscAxx          = $nPhotonsAxx * $e;
	my $JscAyy          = $nPhotonsAyy * $e;
	my $JscAzz          = $nPhotonsAzz * $e;
	my $EffAxx          = $Voc * $JscAxx * $FillFactor / $LightSourcePower * 100.0; # %
	my $EffAyy          = $Voc * $JscAyy * $FillFactor / $LightSourcePower * 100.0; # %
	my $EffAzz          = $Voc * $JscAzz * $FillFactor / $LightSourcePower * 100.0; # %
	printf "  Photon flux absorbed by axx: %10.4g m^-2 (%4.2f mA/cm2)\n", $nPhotonsAxx, $JscAxx*0.1;
	printf "  Photon flux absorbed by ayy: %10.4g m^-2 (%4.2f mA/cm2)\n", $nPhotonsAyy, $JscAyy*0.1;
	printf "  Photon flux absorbed by azz: %10.4g m^-2 (%4.2f mA/cm2)\n", $nPhotonsAzz, $JscAzz*0.1;
	printf "    Supposed Effxx=Voc*JscAxx*FF: %6.2f %\n", $EffAxx;
	printf "    Supposed Effyy=Voc*JscAyy*FF: %6.2f %\n", $EffAyy;
	printf "    Supposed Effzz=Voc*JscAzz*FF: %6.2f %\n", $EffAzz;

	print "  Save the optical spectra to [$OutOpticalSpectraCSV]\n";
	my $out = JFile->new($OutOpticalSpectraCSV, 'w') or die "$!: Can not write to [$OutOpticalSpectraCSV].\n";
	$out->print("WL(nm),E(eV),F(Direct),axx,ayy,azz,Axx,Ayy,Azz,F*Axx,F*Ayy,F*Azz\n");
	for(my $i = 0 ; $i < $nPData ; $i++) {
		my $wl  = $pWL->[$i];
		my $E   = Sci::nmToeV($wl);
		my $axx  = $pMCDaxx->YValLinear($E);
		my $ayy  = $pMCDayy->YValLinear($E);
		my $azz  = $pMCDazz->YValLinear($E);
		my $FD  = $pMCDFDirect->YValLinear($wl);
		$out->print("$pWL->[$i],$E,$FD,$axx,$ayy,$azz,$Axx[$i],$Ayy[$i],$Azz[$i],"
				   ."$FDirectAxx[$i],$FDirectAyy[$i],$FDirectAzz[$i]\n");
	}
	$out->Close();
	print "\n";
}

#===========================================
# Film / substrate model
#===========================================
if(defined $nOpt) {
	printf "Calculated from optical absorption for %6.2f nm mulatilayer interference model.\n", $Thickness;
	my (@ModelE, @erxx, @eixx, @eryy, @eiyy, @erzz, @eizz);
	for(my $i = 0 ; ; $i++) {
		my $E = $dModelE * ($i+1);
		last if($E > $ModelE1);

		$ModelE[$i] = $E;
		$erxx[$i]   = $pMCDerxx->YValLinear($E);
		$eixx[$i]   = $pMCDeixx->YValLinear($E);
		$eryy[$i]   = $pMCDeryy->YValLinear($E);
		$eiyy[$i]   = $pMCDeiyy->YValLinear($E);
		$erzz[$i]   = $pMCDerzz->YValLinear($E);
		$eizz[$i]   = $pMCDeizz->YValLinear($E);
	}
	my $nModelData = @ModelE;

	my $optics = new Optics;
	my $Layers   = new MultiLayer;
	my $Layersyy = new MultiLayer;
	my $Layerszz = new MultiLayer;
	my $Air             = new OpticalMaterial("air");
	my $TopElectrode    = new OpticalMaterial();
	my $Film            = new OpticalMaterial();
	my $Filmyy          = new OpticalMaterial();
	my $Filmzz          = new OpticalMaterial();
	my $BottomElectrode = new OpticalMaterial();
#	my $Substrate       = new OpticalMaterial("air");
	my $Substrate       = new OpticalMaterial();

	$Substrate->SetThickness($SubstrateThickness);
	$Substrate->AddDielectricModel(
		"Substrate",
		"Constant",
		pE    => \@ModelE,
		e1inf => $Substrateer,
		e2inf => $Substrateei,
		);
#print "ret=$ret\n";
	$Film->SetThickness($Thickness * 1.0e-9);
	$Film->AddDielectricModel(
		"Film",
		"EpsArray",
		pE    => \@ModelE,
		pe1 => \@erxx,
		pe2 => \@eixx,
		);
	$Filmyy->SetThickness($Thickness * 1.0e-9);
	$Filmyy->AddDielectricModel(
		"Film",
		"EpsArray",
		pE    => \@ModelE,
		pe1 => \@eryy,
		pe2 => \@eiyy,
		);
	$Filmzz->SetThickness($Thickness * 1.0e-9);
	$Filmzz->AddDielectricModel(
		"Film",
		"EpsArray",
		pE    => \@ModelE,
		pe1 => \@erzz,
		pe2 => \@eizz,
		);
	$TopElectrode->SetThickness($TopElectrodeThickness * 1.0e-9);
	$TopElectrode->AddDielectricModel(
		"TopElectrode",
		"Constant",
		pE    => \@ModelE,
		e1inf => $TopElectrodeer,
		e2inf => $TopElectrodeei,
		);
	$BottomElectrode->SetThickness($BottomElectrodeThickness * 1.0e-9);
	$BottomElectrode->AddDielectricModel(
		"BottomElectrode",
		"Constant",
		pE    => \@ModelE,
		e1inf => $BottomElectrodeer,
		e2inf => $BottomElectrodeei,
		);

	$Layers->AddLayer($Air,             -1.0,                0.0);
	$Layers->AddLayer($TopElectrode,    $TopElectrodeThickness * 1.0e-9, 1.0);
	$Layers->AddLayer($Film,            $Thickness * 1.0e-9, 1.0);
	$Layers->AddLayer($BottomElectrode, $BottomElectrodeThickness * 1.0e-9, 1.0);
	$Layers->AddLayer($Substrate,       $SubstrateThickness, 0.0);
	$Layers->AddLayer($Air,             -1.0,                0.0);
	$Layers->SetIncidentAngle(0.0);

	$Layersyy->AddLayer($Air,             -1.0,                0.0);
	$Layersyy->AddLayer($TopElectrode,    $TopElectrodeThickness * 1.0e-9, 1.0);
	$Layersyy->AddLayer($Filmyy,          $Thickness * 1.0e-9, 1.0);
	$Layersyy->AddLayer($BottomElectrode, $BottomElectrodeThickness * 1.0e-9, 1.0);
	$Layersyy->AddLayer($Substrate,       $SubstrateThickness, 0.0);
	$Layersyy->AddLayer($Air,             -1.0,                0.0);
	$Layersyy->SetIncidentAngle(0.0);

	$Layerszz->AddLayer($Air,             -1.0,                0.0);
	$Layerszz->AddLayer($TopElectrode,    $TopElectrodeThickness * 1.0e-9, 1.0);
	$Layerszz->AddLayer($Filmzz,          $Thickness * 1.0e-9, 1.0);
	$Layerszz->AddLayer($BottomElectrode, $BottomElectrodeThickness * 1.0e-9, 1.0);
	$Layerszz->AddLayer($Substrate,       $SubstrateThickness, 0.0);
	$Layerszz->AddLayer($Air,             -1.0,                0.0);
	$Layerszz->SetIncidentAngle(0.0);

	print "  Save the film dielectric function to [$OutFilmSpectraCSV]\n";
	my $out = JFile->new($OutFilmSpectraCSV, 'w') or die "$!: Can not write to [$OutFilmSpectraCSV].\n";
	$out->print("WL(nm),E(eV),er(xx),ei(xx),er(yy),ei(yy),er(zz),ei(zz)\n");
	for(my $i = 0 ; $i < $nModelData ; $i++) {
		my $E  = $ModelE[$i];
		my $wl = Sci::eVTonm($E);
		$out->print("$wl,$E,$erxx[$i],$eixx[$i],$eryy[$i],$eiyy[$i],$erzz[$i],$eizz[$i]\n");
	}
	$out->Close();

	my (@Axx, @Txx, @Rxx, @Ayy, @Tyy, @Ryy, @Azz, @Tzz, @Rzz, @FDirectAxx, @FDirectAyy, @FDirectAzz);
	for(my $i = 0 ; $i < $nPData ; $i++) {
		my $wl   = $pWL->[$i];
		my $E    = Sci::nmToeV($wl);

		my ($e1, $e2) = $Film->CalEps($E);
		my ($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp)         = $Layers->CalFresnelCoefficient($E);
		$Txx[$i] = $Ts;
		$Rxx[$i] = $Rs;
		$Axx[$i] = 1.0 - $Ts - $Rs;
		my ($rs23, $rp23, $ts23, $tp23, $Rsyy, $Rpyy, $Tsyy, $Tpyy) = $Layersyy->CalFresnelCoefficient($E);
		$Tyy[$i] = $Tsyy;
		$Ryy[$i] = $Rsyy;
		$Ayy[$i] = 1.0 - $Tsyy - $Rsyy;
		my ($rs23, $rp23, $ts23, $tp23, $Rszz, $Rpzz, $Tszz, $Tpzz) = $Layerszz->CalFresnelCoefficient($E);
		$Tzz[$i] = $Tszz;
		$Rzz[$i] = $Rszz;
		$Azz[$i] = 1.0 - $Tszz - $Rszz;
		my $FD   = Algorism::Interpolate($pWL, $pFDirect, $wl);
		$FDirectAxx[$i] = $pFDirect->[$i] * $Axx[$i];
		$FDirectAyy[$i] = $pFDirect->[$i] * $Ayy[$i];
		$FDirectAzz[$i] = $pFDirect->[$i] * $Azz[$i];
	}
	my $pIntFDirectAxx  = Algorism::IntegrateByTrapezoid($pWL, \@FDirectAxx);
	my $nPhotonsAxx     = Algorism::Interpolate($pWL, $pIntFDirectAxx, eVTonm($Eg));
	my $JscAxx          = $nPhotonsAxx * $e;
	my $EffAxx          = $Voc * $JscAxx * $FillFactor / $LightSourcePower * 100.0; # %
	my $pIntFDirectAyy  = Algorism::IntegrateByTrapezoid($pWL, \@FDirectAyy);
	my $nPhotonsAyy     = Algorism::Interpolate($pWL, $pIntFDirectAyy, eVTonm($Eg));
	my $JscAyy          = $nPhotonsAyy * $e;
	my $EffAyy          = $Voc * $JscAyy * $FillFactor / $LightSourcePower * 100.0; # %
	my $pIntFDirectAzz  = Algorism::IntegrateByTrapezoid($pWL, \@FDirectAzz);
	my $nPhotonsAzz     = Algorism::Interpolate($pWL, $pIntFDirectAzz, eVTonm($Eg));
	my $JscAzz          = $nPhotonsAzz * $e;
	my $EffAzz          = $Voc * $JscAzz * $FillFactor / $LightSourcePower * 100.0; # %
	printf "  Photon flux absorbed by axx: %10.4g m^-2 (%4.2f mA/cm2)\n", $nPhotonsAxx, $JscAxx*0.1;
	printf "    Supposed Effxx=Voc*JscAxx*FF: %6.2f %\n", $EffAxx;
	printf "  Photon flux absorbed by ayy: %10.4g m^-2 (%4.2f mA/cm2)\n", $nPhotonsAyy, $JscAyy*0.1;
	printf "    Supposed Effyy=Voc*JscAyy*FF: %6.2f %\n", $EffAyy;
	printf "  Photon flux absorbed by azz: %10.4g m^-2 (%4.2f mA/cm2)\n", $nPhotonsAzz, $JscAzz*0.1;
	printf "    Supposed Effzz=Voc*JscAzz*FF: %6.2f %\n", $EffAzz;

	print "  Save the optical spectra to [$OutMultilayerCSV]\n";
	my $out = JFile->new($OutMultilayerCSV, 'w') or die "$!: Can not write to [$OutMultilayerCSV].\n";
	$out->print("WL(nm),E(eV),erxx,eixx,axx,Txx,Rxx,Axx,F(Direct),F*Axx,eryy,eiyy,ayy,Tyy,Ryy,Ayy,F*Ayy,erzz,eizz,azz,Tzz,Rzz,Azz,F*Azz\n");
	for(my $i = 0 ; $i < $nPData ; $i++) {
		my $wl  = $pWL->[$i];
		my $E   = Sci::nmToeV($wl);
		my ($e1xx, $e2xx) = $Film->CalEps($E);
		my ($rs23, $rp23, $ts23, $tp23, $Rs, $Rp, $Ts, $Tp) = $Layers->CalFresnelCoefficient($E);
		my $axx = $pMCDaxx->YValLinear($E);
		my ($e1yy, $e2yy) = $Filmyy->CalEps($E);
		my ($rs23, $rp23, $ts23, $tp23, $Rsyy, $Rpyy, $Tsyy, $Tpy) = $Layersyy->CalFresnelCoefficient($E);
		my $ayy = $pMCDayy->YValLinear($E);
		my ($e1zz, $e2zz) = $Filmzz->CalEps($E);
		my ($rs23, $rp23, $ts23, $tp23, $Rszz, $Rpzz, $Tszz, $Tpzz) = $Layerszz->CalFresnelCoefficient($E);
		my $azz = $pMCDazz->YValLinear($E);
		my $FD  = $pMCDFDirect->YValLinear($wl);
		$out->print("$pWL->[$i],$E,$e1xx,$e2xx,$axx,$Ts,$Rs,$Axx[$i],$FD,$FDirectAxx[$i],$e1yy,$e2yy,$ayy,$Tsyy,$Rsyy,$Ayy[$i],$FDirectAyy[$i],$e1zz,$e2zz,$azz,$Tszz,$Rszz,$Azz[$i],$FDirectAzz[$i]\n");
	}
	$out->Close();
	print "\n";
}

exit;

#===========================
# Subroutines
#===========================
sub ReadSolarSpectra
{
	my ($FileName) = @_;

	my $in = new JFile;
	return undef if(!$in->Open($FileName, "r"));

	my (@X, @Y1, @Y2);
	my $FirstLine = $in->ReadLine();
	my $count = 0;
	while(!$in->eof()) {
		my $line = $in->ReadLine();
		last if($in->eof());
		my @a = Utils::Split("\\s+", $line);
#print "$a[0]\t$a[1]\t$a[2]\n";
		$X[$count]  = $a[0];
		$Y1[$count] = $a[1];
		$Y2[$count] = $a[2];
		$count++;
	}
	$in->Close();

# lambda (nm), Solar(All), Solar(Direct)
	my $pMCDY1  = new MultiColumnData(\@X, \@Y1);
	my $pMCDY2  = new MultiColumnData(\@X, \@Y2);

	return (\@X, $pMCDY1, $pMCDY2);
}

sub ReadOpticalSpectra
{
	my ($path) = @_;

	my $in    = new MultiColumnData($path, 1) or die "$!: Can not read [$path].\n";
#	my $nDataArray = $in->nDataArray();
	my $nData = $in->nData();
	my $pE    = $in->GetData('E.*');
	my $pWL   = $in->GetData('wl.*');
	$pMCDaxx  = $in->Duplicate('E.*', 'a(xx)');
	$pMCDayy  = $in->Duplicate('E.*', 'a(yy)');
	$pMCDazz  = $in->Duplicate('E.*', 'a(zz)');
	$pMCDerxx = $in->Duplicate('E.*', 'er(xx)');
	$pMCDeixx = $in->Duplicate('E.*', 'ei(xx)');
	$pMCDeryy = $in->Duplicate('E.*', 'er(yy)');
	$pMCDeiyy = $in->Duplicate('E.*', 'ei(yy)');
	$pMCDerzz = $in->Duplicate('E.*', 'er(zz)');
	$pMCDeizz = $in->Duplicate('E.*', 'ei(zz)');
	return ($nData, $pE, $pWL, $pMCDaxx, $pMCDayy, $pMCDazz, $pMCDerxx, $pMCDeixx, $pMCDeryy, $pMCDeiyy, $pMCDerzz, $pMCDeizz);
}