#!/usr/bin/perl

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

BEGIN {
my $BaseDir = $ENV{'PerlDir'};
$BaseDir = $ENV{'TkPerlDir'} unless(defined $BaseDir);
@INC = ("$BaseDir/lib", "$BaseDir/VNL", @INC);
}

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

use CSV;
use JFile;
use MyHTMLApplication;
use Sci qw(nmToeV eVTonm $e);
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 = "$ENV{'PerlDir'}/Color/data" 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          = Deps::MakePath($DataPath, "SolarSpectra.txt", 0);

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, $pPAll, $pPDirect) = &ReadSolarSpectra($SolarSpectraPath);
my $nPData       = @$pWL;
#IntegrateBySimpson
my $pIntPAll     = Algorism::IntegrateByTrapezoid($pWL, $pPAll);
my $pIntPDirect  = Algorism::IntegrateByTrapezoid($pWL, $pPDirect);
my $TotalPAll    = $pIntPAll->[$nPData-1];
my $TotalPDirect = $pIntPDirect->[$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";
for(my $i = 0 ; $i < $nPData ; $i++) {
	$pPAll->[$i]       *= $LightSourcePower / $TotalPAll;
	$pPDirect->[$i]    *= $LightSourcePower / $TotalPDirect;
	$pIntPAll->[$i]    *= $LightSourcePower / $TotalPAll;
	$pIntPDirect->[$i] *= $LightSourcePower / $TotalPDirect;
}
print  "    Calculate photon flux...\n";
my (@FAll, @FDirect);
for(my $i = 0 ; $i < $nPData ; $i++) {
	my $wl = $pWL->[$i];
	my $E  = Sci::nmToeV($wl);
	$FAll[$i]    = $pPAll->[$i] / $E / $e;
	$FDirect[$i] = $pPDirect->[$i] / $E / $e;
#print "wl=$wl\tE=$E\t$FAll[$i]\t$FDirect[$i]\n";
}
my $pIntFAll     = Algorism::IntegrateByTrapezoid($pWL, \@FAll);
my $pIntFDirect  = Algorism::IntegrateByTrapezoid($pWL, \@FDirect);

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),P(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],"
			   ."$FAll[$i],$FDirect[$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, $paxx, $payy, $pazz, $perxx,$peixx);
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, $paxx, $payy, $pazz, 
		$perxx, $peixx) = &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  = Algorism::Interpolate($pOptE, $paxx, $E);
		my $ayy  = Algorism::Interpolate($pOptE, $payy, $E);
		my $azz  = Algorism::Interpolate($pOptE, $pazz, $E);
#		my $axx  = Algorism::Interpolate($pOptWL, $paxx, $wl);
#		my $ayy  = Algorism::Interpolate($pOptWL, $payy, $wl);
#		my $azz  = Algorism::Interpolate($pOptWL, $pazz, $wl);
#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, \@FDirect, $wl);
		$FDirectAxx[$i] = $FDirect[$i] * $Axx[$i];
		$FDirectAyy[$i] = $FDirect[$i] * $Ayy[$i];
		$FDirectAzz[$i] = $FDirect[$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  = Algorism::Interpolate($pOptE, $paxx, $E);
		my $ayy  = Algorism::Interpolate($pOptE, $payy, $E);
		my $azz  = Algorism::Interpolate($pOptE, $pazz, $E);
#		my $axx = Algorism::Interpolate($pOptWL, $paxx, $wl);
#		my $ayy = Algorism::Interpolate($pOptWL, $payy, $wl);
#		my $azz = Algorism::Interpolate($pOptWL, $pazz, $wl);
		my $FD  = Algorism::Interpolate($pWL, \@FDirect, $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: Not implemented
#===========================================
if(defined $nOpt) {
	printf "Calculated from optical absorption for %6.2f nm mulatilayer interference model.\n", $Thickness;
	my (@ModelE, @erxx, @eixx);
	for(my $i = 0 ; ; $i++) {
		my $E = $dModelE * ($i+1);
		last if($E > $ModelE1);

		$ModelE[$i] = $E;
		$erxx[$i]   = Algorism::Interpolate($pOptE, $perxx, $E);
		$eixx[$i]   = Algorism::Interpolate($pOptE, $peixx, $E);
	}
	my $nModelData = @ModelE;

	my $optics = new Optics;
	my $Layers = new MultiLayer;
	my $Air             = new OpticalMaterial("air");
	my $TopElectrode    = new OpticalMaterial();
	my $Film            = new OpticalMaterial();
	my $BottomElectrode = new OpticalMaterial();
#	my $Substrate       = new OpticalMaterial("air");
	my $Substrate       = new OpticalMaterial();
	$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);

	$Substrate->SetThickness($SubstrateThickness);
	my $ret = $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,
		);
	$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,
		);

	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,ei\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]\n");
	}
	$out->Close();

	my (@Axx, @Txx, @Rxx, @FDirectAxx);
	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 $FD   = Algorism::Interpolate($pWL, \@FDirect, $wl);
		$FDirectAxx[$i] = $FDirect[$i] * $Axx[$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; # %
	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;

	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\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 = Algorism::Interpolate($pOptE, $paxx, $E);
		my $FD  = Algorism::Interpolate($pWL, \@FDirect, $wl);
		$out->print("$pWL->[$i],$E,$e1xx,$e2xx,$axx,$Ts,$Rs,$Axx[$i],$FD,$FDirectAxx[$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)
	return (\@X, \@Y1, \@Y2);
}

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

	my $csv = new CSV;
	$csv->Read($path) or die "$!: Can not read [$path].\n";
	my $pLabelArray = $csv->LabelArray();
	my $nData = $csv->nData();
print "    nData: $nData\n";

	my $pDataArray  = $csv->DataArray();
	print "    ";
for(my $il = 0 ; $il < @$pLabelArray ; $il++) {
	print "," if($il > 0);
	print "$pLabelArray->[$il]";
}
print "\n";
#exit;

	my $iE    = $csv->FindLabelIndex(['E.*']);
	my $iWL   = $csv->FindLabelIndex(['wl.*']);
	my $iaxx  = $csv->FindLabelIndex(['a(xx)']);
	my $iayy  = $csv->FindLabelIndex(['a(yy)']);
	my $iazz  = $csv->FindLabelIndex(['a(zz)']);
	my $ierxx = $csv->FindLabelIndex(['er(xx)']);
	my $ieixx = $csv->FindLabelIndex(['ei(xx)']);
print "    idexes: iE=$iE  iWL=$iWL  iaxx=$iaxx  iayy=$iayy  iazz=$iazz\n";

	my $pE    = $pDataArray->[$iE];
	my $pWL   = $pDataArray->[$iWL];
	my $paxx  = $pDataArray->[$iaxx];
	my $payy  = $pDataArray->[$iayy];
	my $pazz  = $pDataArray->[$iazz];
	my $perxx = $pDataArray->[$ierxx];
	my $peixx = $pDataArray->[$ieixx];
#	for(my $i = 0 ; $i < $nData / 2 ; $i++) {
#		($pE->[$i],   $pE->[$nData-1-$i])   = Utils::Swap($pE->[$i],   $pE->[$nData-1-$i]);
#		($pWL->[$i],  $pWL->[$nData-1-$i])  = Utils::Swap($pWL->[$i],  $pWL->[$nData-1-$i]);
#		($paxx->[$i], $paxx->[$nData-1-$i]) = Utils::Swap($paxx->[$i], $paxx->[$nData-1-$i]);
#		($payy->[$i], $payy->[$nData-1-$i]) = Utils::Swap($payy->[$i], $payy->[$nData-1-$i]);
#		($pazz->[$i], $pazz->[$nData-1-$i]) = Utils::Swap($pazz->[$i], $pazz->[$nData-1-$i]);
#	}
#for(my $i = 0 ; $i < $nData ; $i++) {
#print "$pWL->[$i]\t$paxx->[$i]\n";
#}

	return ($nData, $pE, $pWL, $paxx, $payy, $pazz, $perxx, $peixx);
}