#!/usr/bin/perl

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

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

use strict;
#use warnings;

use Math::Matrix;
use Math::MatrixReal;
use PDL::Opt::Simplex;

use Deps;
use Utils;
use JFile;

use MyApplication;
use CSV;
use Sci::Science;

#===============================================
# デバッグ関係変数
#===============================================
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# 大域変数
#===============================================
my $kB = Sci::kB();
my $e  = Sci::e();

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

my $LF        = Deps::LF();
my $DirSep    = Deps::DirSep();
my $RegDirSep = Deps::RegDirSep();

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

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

#===============================================
# スクリプト大域変数
#===============================================
my $InitialDirectory = Deps::GetWorkingDirectory();

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",
		"--Action=[FDFit|Deconvolution|Convolution]",       '');
$App->AddArgument("--Temperature",   "--Temperature=val [Def: 0]", 0);
$App->AddArgument("--AppFunction",   "--AppFunction=[Lorentizan|Gaussian] [Def: Lorentzian]", 'Lorentzian');
$App->AddArgument("--Wa",         "--Wa=val [Def: 0]", 0);
$App->AddArgument("--EF",         "--EF=val [Def: 0]", 0);
$App->AddArgument("--C0",         "--C0=val [Def: 0]", 0);
$App->AddArgument("--BG",         "--BG=val [Def: 0]", 0);
$App->AddArgument("--BGLeft",     "--BGLeft=val [Def: '']", '');
$App->AddArgument("--EBG0",       "--EBG0=val [Def: -0.1]", -0.1);
$App->AddArgument("--WBG",        "--WBG=val [Def: 0.1]", 0.1);
$App->AddArgument("--DOS",        "--DOS=[CB|Metal] [Def: CB]", 'CB');
$App->AddArgument("--ECBM",       "--ECBM=val [Def: -0.4]", -0.4);
$App->AddArgument("--D0CB",       "--D0CB=val [Def: 0.0]", 0.0);
$App->AddArgument("--EVBM",       "--EVBM=val [Def: -3]", -3.0);
$App->AddArgument("--D0VB",       "--D0VB=val [Def: 0.0]", 0.0);
$App->AddArgument("--Dumping",    "--Dumping=val [Def: 1.0]", 1.0);
$App->AddArgument("--nMaxIter",   "--nMaxIter=integer [Def: 10]", 10);
$App->AddArgument("--EPS",        "--EPS=val [Def: 1.0e-6]", 1.0e-6);
$App->AddArgument("--XMin",       "--XMin=val [Def: 0]", 0);
$App->AddArgument("--XMax",       "--XMin=val [Def: 0]", 0);
$App->AddArgument("--IgnoreZeroValue", "--IgnoreZeroValue=[0|1] [Def: 1]", '');
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my $Args = $App->Args();
#my $form = new CGI;
#$Args->SetCGIForm($form);
#$Args->parseInput($WebCharCode);

#==========================================
# メイン関数スタート
#==========================================

#Utils::InitHTML("Research", $WebCharSet, "_self");

my $Debug = $Args->GetGetArg("DebugMode");
$App->SetDebug($Debug);
my $Action = $Args->GetGetArg("Action");

my $ret = 0;
if($Action =~ /FDFit/i) {
	&FDFit();
}
elsif($Action =~ /Deconvolution/i) {
	&Deconvolution();
}
elsif($Action =~ /Convolution/i) {
	&Convolution();
}
else {
	$App->print("Error: Invald Action: $Action\n");
}

#Utils::EndHTML();

exit $ret;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================

my ($pXObs, $pYObs);
my ($XMin, $XMax);
my $T;
#my ($C0, $EF, $Wa);

sub BG
{
	my ($E, $BG, $BGLeft, $EBG0, $WBG) = @_;
	return $BG + ($BGLeft-$BG) / (1.0 + exp(($E-$EBG0)/$WBG));
}

sub DOS
{
	my ($E, $D0CB, $ECBM, $D0VB, $EVBM, $DOSType) = @_;
	my $DE = 0.0;
	if($DOSType =~ /CB/i) {
		$DE = $D0CB * sqrt($E - $ECBM) if($E > $ECBM);
	}
	if($DOSType =~ /VB/i) {
		$DE = $D0VB * sqrt($EVBM - $E) if($E < $EVBM);
	}
	if($DOSType =~ /Metal/i) {
		$DE = $D0CB if($E <= $ECBM);
	}
	return $DE;
}

sub FermiDiracFunction
{
	my ($E, $T, $EF, $C0, $Wa) = @_; # eV
#print "E=$E EF=$EF kB=$kB T=$T e=$e Wa=$Wa\n";
	return $C0 / (1.0 + exp(($E-$EF)/($kB*$T/$e+$Wa)));
}

sub S2
{
	my ($pVariables) = @_;
	my $EF = $pVariables->[0];
	my $C0 = $pVariables->[1];
	my $BG = $pVariables->[2];
	my $Wa = $pVariables->[3];
	my $nData = @$pXObs;
#print "Var: $EF $C0 $BG $Wa\n";

	my $c = 0;
	my $S2 = 0.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x    = $pXObs->[$i];
		my $yobs = $pYObs->[$i];
		next if($x < $XMin or $x > $XMax);
		my $ycal = FermiDiracFunction($x, $T, $EF, $C0, $Wa) + $BG;
		my $d = $ycal - $yobs;
#print "[$i]: ($x, $yobs, $ycal)\n";
		$S2 += $d * $d;
		$c++;
	}
	return sqrt($S2 / $c);
}

sub Differentiate1
{
	my ($pFunc, $pVariables, $idx) = @_;
	my $dv = abs($pVariables->[$idx]) * 0.001;
	$dv = 0.001 if($dv < 0.001);

	my @v;
	for(my $i = 0 ; $i < @$pVariables ; $i++) {
		$v[$i] = $pVariables->[$i];
	}
	$v[$idx] = $pVariables->[$idx] + $dv;

	my $S0 = &$pFunc($pVariables);
	my $S1 = &$pFunc(\@v);
	return ($S1 - $S0) / $dv;
}

sub Differentiate2
{
	my ($pFunc, $pVariables, $idx, $jdx) = @_;
	my $dv = abs($pVariables->[$jdx]) * 0.001;
	$dv = 0.001 if($dv < 0.001);

	my @v;
	for(my $i = 0 ; $i < @$pVariables ; $i++) {
		$v[$i] = $pVariables->[$i];
	}
	$v[$jdx] = $pVariables->[$jdx] + $dv;

	my $S0 = Differentiate1($pFunc, $pVariables, $idx);
	my $S1 = Differentiate1($pFunc, \@v, $idx);

	return ($S1 - $S0) / $dv;
}

sub Optimize
{
	my ($pFunc, $pVariables, $nMaxIter, $EPS) = @_;
	$nMaxIter = 10 if(!defined $nMaxIter);
	$EPS = 1.0e-6 if(!defined $EPS);

	my $n = @$pVariables;
	for(my $iter = 0 ; $iter < $nMaxIter ; $iter++) {
		my $S = &$pFunc($pVariables);
		print "Iter[$iter]: S=$S\n";

#一次微分の計算
		my @v;
		my $Si = new Math::MatrixReal($n, 1);
		for(my $i = 0 ; $i < $n ; $i++) {
			my $val = Differentiate1($pFunc, $pVariables, $i);
#print "dFdx[$i]=$val\n";
			$Si->assign($i+1, 1, $val);
		}
#print "Si:\n";
#print $Si;

#二次微分の計算
		my $Sij = new Math::MatrixReal($n, $n);
		for(my $i = 0 ; $i < $n ; $i++) {
			for(my $j = $i ; $j < $n ; $j++) {
				my $val = Differentiate2($pFunc, $pVariables, $i, $j);
#print "dF2dx2[$i][$j]=$val\n";
				$Sij->assign($i+1, $j+1, $val);
				$Sij->assign($j+1, $i+1, $val) if($i != $j);
			}
		}
		my @kj;
		for(my $j = 0 ; $j < $n ; $j++) {
			$kj[$j] = 1.0 / $Sij->element($j+1, $j+1);
		}

#Dumping factorを変えながら最小点を探す
		my $MinS = $S;
		my $PrevS = 1.0e100;
		my $MinDump = -1;
		my @MinV;
		for(my $id = 0 ; $id <= 10 ; $id++) {
			my $Dump = 0.01 * (2**abs($id));
			my $Sij2 = new Math::MatrixReal($n, $n);
			for(my $i = 0 ; $i < $n ; $i++) {
				for(my $j = 0 ; $j < $n ; $j++) {
					$Sij2->assign($i+1, $j+1, $kj[$j] * $Sij->element($i+1, $j+1));
				}
				$Sij2->assign($i+1, $i+1, $Sij2->element($i+1, $i+1) + $Dump);
			}
#print "Sij:\n";
#print $Sij;
			my $RSij2 = $Sij2->inverse();
#print "Inverse Sij:\n";
#print $RSij;
			my $V = -$RSij2 * $Si;
			$V = -$V if($id < 0);
#print "V:\n";
#print $V;
			my @NewV;
			for(my $i = 0 ; $i < $n ; $i++) {
				$NewV[$i] = $pVariables->[$i] + $V->element($i+1, 1) * $kj[$i];
#print "  v[$i](d=$Dump): $pVariables->[$i] => $NewV[$i]\n";
			}
			my $NewS = &$pFunc(\@NewV);
#print "  NewS(d=$Dump): $NewS\n";
			if($NewS < $MinS) {
				$MinS = $NewS;
				$MinDump = $Dump;
				for(my $i = 0 ; $i < $n ; $i++) {
					$MinV[$i] = $NewV[$i];
				}
			}
			last if($NewS > $PrevS);
			$PrevS = $NewS;
		}

		if(defined $MinV[0]) {
			my $MaxD = 0;
			for(my $i = 0 ; $i < $n ; $i++) {
#収束判定条件のチェック
				my $k = abs($MinV[$i]);
				if($k < abs($pVariables->[$i])) {
					$k = abs($pVariables->[$i]);
				}
				$k = 1.0 if($k == 0.0);
				my $d = ($MinV[$i] - $pVariables->[$i]) / $k;
				if($MaxD < abs($d)) {
					$MaxD = abs($d);
				}

print "  v[$i](d=$MinDump): $pVariables->[$i] => ";
				$pVariables->[$i] = $MinV[$i];
print "$MinV[$i]\n";
			}
#収束判定
			if($MaxD < $EPS) {
				print "  Convergence reached.: Stop Optimization\n";
				last;
			}
		}
		else {
#最小点がみつからなかったら止める
			print "  Warning: Better solution could not be found.\n";
			print "    Stop Optimization\n";
			last;
		}
	}

	return 0.0;
}

sub FDFit
{
#print "300 K = ", 300.0*$kB / $e, "\n";

	$App->print("\n\n<b>Fitting XPS Data to Fermi-Dirac distribution by Simplex method:</b>\n");
	my $InputFile = $Args->GetGetArg(0);
	$App->print("  Input File: $InputFile\n");
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($InputFile);
	my $OutCSVFile = Deps::MakePath("$drive$dir", "${filebody}-out.csv", 0);
	$App->print("  Output CSV File: $OutCSVFile\n");

	$XMin = $Args->GetGetArg("XMin") + 0.0;
	$XMax = $Args->GetGetArg("XMax") + 0.0;
	$App->print("  X range: $XMin - $XMax eV\n");
	$T  = $Args->GetGetArg("Temperature") + 0.0;
	my $Wa = $Args->GetGetArg("Wa") + 0.0;
	my $C0 = $Args->GetGetArg("C0") + 0.0;
	my $BG = $Args->GetGetArg("BG") + 0.0;
	my $EF = $Args->GetGetArg("EF") + 0.0;
	$App->print("  Temperature: $T K\n");
	$App->print("  Initial parameters:\n");
	$App->print("    Apparatus resolution: $Wa eV\n");
	$App->print("    Scaling factor: $C0\n");
	$App->print("    Background: $BG\n");
	$App->print("    Fermi energy: $EF eV\n");
	my $nMaxIter = $Args->GetGetArg("nMaxIter");
	$App->print("  nMaxIter: $nMaxIter\n");
	my $EPS = $Args->GetGetArg("EPS");
	$App->print("  EPS: $EPS\n");

	$App->print("  Read [$InputFile]\n");
	my $csv = new CSV;
	if($csv->Read($InputFile)) {
	}
	else {
		$App->print("  Error!!: Can not read [$InputFile]\n");
		return -1;
	}
	my $pLabelArray = $csv->LabelArray();
	my $pDataArray  = $csv->DataArray();
	$pXObs = $pDataArray->[0];
	$pYObs = $pDataArray->[1];
	my $nData = @$pXObs;
$App->print("    nData   : $nData\n");
	for(my $i = 0 ; $i < $nData ; $i++) {
		$pXObs->[$i] *= -1;
	}
	$csv->Close();

	my $pVariables = [$EF, $C0, $BG, $Wa];
	my $S2 = Optimize(\&S2, $pVariables, $nMaxIter, $EPS);

	$App->print("  Write to [$OutCSVFile]\n");
	my $out = new JFile();
	if(!$out->Open($OutCSVFile, "w")) {
		$App->print("  Error:: Can not write to [$OutCSVFile].\n");
		return;
	}
	$out->print("E(eV),Obs,Cal(ini),Cal(Fin)\n");
	my $EF1 = $pVariables->[0];
	my $C01 = $pVariables->[1];
	my $BG1 = $pVariables->[2];
	my $Wa1 = $pVariables->[3];
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x    = $pXObs->[$i];
		my $yobs = $pYObs->[$i];
		next if($x < $XMin or $x > $XMax);
		my $ycal0 = FermiDiracFunction($x, $T, $EF, $C0, $Wa) + $BG;
		my $ycal1 = FermiDiracFunction($x, $T, $EF1, $C01, $Wa1) + $BG1;
		$out->print("$x,$yobs,$ycal0,$ycal1\n");
	}
	$out->Close();

	$App->print("\n\n<b>Fitting XPS Data to Fermi-Dirac distribution by Simplex method: Finished.</b>\n");
}

sub Deconvolution
{
	$App->print("\n<b>Deconvolute CSV File:</b>\n");

	my $InputFile = $Args->GetGetArg(0);
	$App->print("  Input File: $InputFile\n");
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($InputFile);
	my $OutCSVFile = Deps::MakePath("$drive$dir", "${filebody}-out.csv", 0);
	$App->print("  Output CSV File: $OutCSVFile\n");

	my $AppFunction = $Args->GetGetArg("AppFunction");
	$AppFunction = 'Lorentian' if(!defined $AppFunction);
	$App->print("  AppFunction: $AppFunction\n");

	$XMin = $Args->GetGetArg("XMin") + 0.0;
	$XMax = $Args->GetGetArg("XMax") + 0.0;
	$App->print("  X range: $XMin - $XMax eV\n");
	my $Wa = $Args->GetGetArg("Wa") + 0.0;
	my $EF = $Args->GetGetArg("EF") + 0.0;
	my $Dumping = $Args->GetGetArg("Dumping");
	$App->print("  Apparatus resolution: $Wa eV\n");
	$App->print("  Fermi energy: $EF eV\n");
	$App->print("  Dumping Factor: $Dumping\n");
	my $nMaxIter = $Args->GetGetArg("nMaxIter");
	$App->print("  nMaxIter: $nMaxIter\n");
	my $EPS = $Args->GetGetArg("EPS");
	$App->print("  EPS: $EPS\n");

	my $T  = $Args->GetGetArg("Temperature") + 0.0;
	my $C0 = $Args->GetGetArg("C0") + 0.0;
	my $BG = $Args->GetGetArg("BG") + 0.0;
	$App->print("  Temperature: $T\n");
	$App->print("  C0: $C0\n");
	$App->print("  BG: $BG\n");

	$App->print("  Read [$InputFile]\n");
	my $csv = new CSV;
	if($csv->Read($InputFile)) {
	}
	else {
		$App->print("  Error!!: Can not read [$InputFile]\n");
		return -1;
	}
	my $pLabelArray = $csv->LabelArray();
	my $pDataArray  = $csv->DataArray();
	$pXObs = $pDataArray->[0];
	$pYObs = $pDataArray->[1];
	my $nData = @$pXObs;
$App->print("    nData   : $nData\n");
	for(my $i = 0 ; $i < $nData ; $i++) {
		$pXObs->[$i] *= -1;
	}
	$csv->Close();

	my $CorrectZero = 1;

	my @Y;
	$Y[0] = [];
	for(my $i = 0 ; $i < $nData ; $i++) {
		$Y[0]->[$i] = $pYObs->[$i];
		$Y[0]->[$i] = 0.0 if($CorrectZero and $Y[0]->[$i] < 0.0);
	}

	my $AppFunc = sub { Sci::Lorentzian(@_); };
	if($AppFunction =~ /^g/i) {
		print "Use Gaussian\n";
		$AppFunc = sub { Sci::Gaussian(@_); };
	}
	else {
		print "Use Lorentian\n";
	}

	my $dX = $pXObs->[1] - $pXObs->[0];
	my $hii = &$AppFunc(0.0, 0.0, $Wa) * $dX;
	my $K = $Dumping / $hii;
print "hii=$hii K=$K\n";
	for(my $iter = 1 ; $iter <= $nMaxIter ; $iter++) {
		my $iprev = $iter-1;
		$Y[$iter] = [];
		for(my $i = 0 ; $i < $nData ; $i++) {
			my $xi = $pXObs->[$i];
			my $S0 = 0.0;
			if($i > 0) {
				for(my $j = 0 ; $j <= $i-1 ; $j++) {
					my $xj = $pXObs->[$j];
					my $hij = &$AppFunc($xi, $xj, $Wa) * $dX;
					$S0 += $hij * $Y[$iter]->[$j];
				}
			}
			for(my $j = $i ; $j < $nData ; $j++) {
				my $xj = $pXObs->[$j];
				my $hij = &$AppFunc($xi, $xj, $Wa) * $dX;
				$S0 += $hij * $Y[$iprev]->[$j];
			}
			$Y[$iter]->[$i] = $Y[$iprev]->[$i] + $K * ($pYObs->[$i] - $S0);
			$Y[$iter]->[$i] = 0.0 if($CorrectZero and $Y[$iter]->[$i] < 0.0);
		}
	}

	$App->print("  Write to [$OutCSVFile]\n");
	my $out = new JFile();
	if(!$out->Open($OutCSVFile, "w")) {
		$App->print("  Error:: Can not write to [$OutCSVFile].\n");
		return;
	}
	my $nInterval = int($nMaxIter / 10 + 0.1);
	$out->print("E(eV),");
	for(my $i = 0 ; $i <= $nMaxIter ; $i += $nInterval) {
		my $c = $i;
		$out->print("Cal($c),");
	}
	$out->print("AppFunc,FD(E)\n");
	my $xcenter = ($XMin + $XMax) / 2.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x    = $pXObs->[$i];
		my $yobs = $pYObs->[$i];
		next if($x < $XMin or $x > $XMax);
		$out->print("$x,");
		for(my $j = 0 ; $j <= $nMaxIter ; $j += $nInterval) {
			$out->print("$Y[$j]->[$i],");
		}
		my $yapp = &$AppFunc($x, $xcenter, $Wa);
		$yapp = 0 if(abs($yapp) < 1.0e-100);
		my $FD = FermiDiracFunction($x, $T, $EF, $C0, 0.0) + $BG;
		$out->print("$yapp,$FD\n");
	}
	$out->Close();


	$App->print("\n\n<b>Deconvolute CSV File: Finished</b>\n");
}

sub Convolution
{
	$App->print("\n<b>Convolute CSV File:</b>\n");

	my $InputFile = $Args->GetGetArg(0);
	$App->print("  Input File: $InputFile\n");
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($InputFile);
	my $OutCSVFile = Deps::MakePath("$drive$dir", "${filebody}-out.csv", 0);
	$App->print("  Output CSV File: $OutCSVFile\n");

	my $AppFunction = $Args->GetGetArg("AppFunction");
	$AppFunction = 'Lorentian' if(!defined $AppFunction);
	$App->print("  AppFunction: $AppFunction\n");

	$XMin = $Args->GetGetArg("XMin") + 0.0;
	$XMax = $Args->GetGetArg("XMax") + 0.0;
	$App->print("  X range: $XMin - $XMax eV\n");
	my $Wa = $Args->GetGetArg("Wa") + 0.0;
	my $EF = $Args->GetGetArg("EF") + 0.0;
	$App->print("  Apparatus resolution: $Wa eV\n");
	$App->print("  Fermi energy: $EF eV\n");

	my $Dumping = $Args->GetGetArg("Dumping");
#	$App->print("  Dumping Factor: $Dumping\n");
	my $nMaxIter = $Args->GetGetArg("nMaxIter");
#	$App->print("  nMaxIter: $nMaxIter\n");
	my $EPS = $Args->GetGetArg("EPS");
#	$App->print("  EPS: $EPS\n");

	my $T  = $Args->GetGetArg("Temperature") + 0.0;
	my $C0 = $Args->GetGetArg("C0") + 0.0;
	$App->print("  Temperature: $T\n");
	$App->print("  C0: $C0\n");

	my $BG     = $Args->GetGetArg("BG") + 0.0;
	my $BGLeft = $Args->GetGetArg("BGLeft") + 0.0;
	$BGLeft = $BG if($BGLeft eq '');
	my $EBG0   = $Args->GetGetArg("EBG0") + 0.0;
	my $WBG    = $Args->GetGetArg("WBG") + 0.0;
	$App->print("  BG    : $BG\n");
	$App->print("  BGLeft: $BGLeft\n");
	$App->print("  EBG0  : $EBG0 eV\n");
	$App->print("  WBG   : $BG eV\n");

	my $DOS = $Args->GetGetArg("DOS");
	$DOS = 'Metal' if(!defined $DOS);
	my $ECBM = $Args->GetGetArg("ECBM") + 0.0;
	$ECBM = -0.4 if(!defined $ECBM);
	my $D0CB = $Args->GetGetArg("D0CB") + 0.0;
	$D0CB = 0 if(!defined $D0CB);
	my $EVBM = $Args->GetGetArg("EVBM") + 0.0;
	$EVBM = -3 if(!defined $EVBM);
	my $D0VB = $Args->GetGetArg("D0VB") + 0.0;
	$D0VB = 0 if(!defined $D0VB);
	$App->print("  DOS: $DOS\n");
	$App->print("  ECBM: $ECBM eV\n");
	$App->print("  D0CB: $D0CB\n");
	$App->print("  EVBM: $EVBM eV\n");
	$App->print("  D0VB: $D0VB\n");

	$App->print("  Read [$InputFile]\n");
	my $csv = new CSV;
	if($csv->Read($InputFile)) {
	}
	else {
		$App->print("  Error!!: Can not read [$InputFile]\n");
		return -1;
	}
	my $pLabelArray = $csv->LabelArray();
	my $pDataArray  = $csv->DataArray();
	$pXObs = $pDataArray->[0];
	$pYObs = $pDataArray->[1];
	my $nData = @$pXObs;
$App->print("    nData   : $nData\n");
	for(my $i = 0 ; $i < $nData ; $i++) {
		$pXObs->[$i] *= -1;
	}
	$csv->Close();

	my $AppFunc = sub { Sci::Lorentzian(@_); };
	if($AppFunction =~ /^g/i) {
		print "Use Gaussian\n";
		$AppFunc = sub { Sci::Gaussian(@_); };
	}
	else {
		print "Use Lorentian\n";
	}

	my $dX = $pXObs->[1] - $pXObs->[0];
	my @Y;
	for(my $i = 0 ; $i < $nData ; $i++) {
		$Y[$i] = 0.0;
	}
	my $nExt = int($Wa * 10.0 / $dX + 0.01);
	for(my $i = -$nExt ; $i <= $nData+$nExt ; $i++) {
		my $x0 = $pXObs->[0] + $dX * $i;
		my $FD = FermiDiracFunction($x0, $T, $EF, 1.0, 0.0);
		my $DE = DOS($x0, $D0CB, $ECBM, $D0VB, $EVBM, $DOS);
		for(my $j = 0 ; $j < $nData ; $j++) {
			my $xj = $pXObs->[$j];
			my $hij = &$AppFunc($x0, $xj, $Wa) * $dX;
			$Y[$j] += $hij * $FD * $DE;
		}
	}
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x0 = $pXObs->[0] + $dX * $i;
		$Y[$i] += BG($x0, $BG, $BGLeft, $EBG0, $WBG);
	}

	$App->print("  Write to [$OutCSVFile]\n");
	my $out = new JFile();
	if(!$out->Open($OutCSVFile, "w")) {
		$App->print("  Error:: Can not write to [$OutCSVFile].\n");
		return;
	}
	my $nInterval = int($nMaxIter / 10 + 0.1);
	$out->print("E(eV),yobs,D(E)FD(E),yconv,D(E),BG(E),FD(E),AppFunc\n");
	my $xcenter = ($XMin + $XMax) / 2.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x    = $pXObs->[$i];
		my $yobs = $pYObs->[$i];
		next if($x < $XMin or $x > $XMax);
		my $FD = FermiDiracFunction($x, $T, $EF, 1.0, 0.0);
		my $DE = DOS($x, $D0CB, $ECBM, $D0VB, $EVBM, $DOS);
		my $BGE  = BG($x, $BG, $BGLeft, $EBG0, $WBG);
		my $DEFD = $FD * $DE + $BGE;
		my $yapp = &$AppFunc($x, $xcenter, $Wa);
		$yapp = 0 if(abs($yapp) < 1.0e-100);
		$out->print("$x,$yobs,$DEFD,$Y[$i],$DE,$BGE,$FD,$yapp\n");
	}
	$out->Close();


	$App->print("\n\n<b>Convolute CSV File: Finished</b>\n");
}