package PES;
use Common;
@ISA = qw(Common);

#公開したいサブルーチン
@EXPORT = qw();

use strict;
#use warnings;

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

use Deps;
use Utils;
use JFile;

use CSV;
use Sci::Science;

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

#============================================================
#　変数等取得関数
#============================================================

#============================================================
#　コンストラクタ、デストラクタ
#============================================================
sub new
{
	my ($module, $app) = @_;
	my $this = {};
	bless $this;
	$this->InitializeArgs();
	$this->SetApplication($app) if($app);
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　一般関数
#============================================================
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 BG
{
	my ($E, $BG, $BGLeft, $EBG0, $WBG) = @_;
	return $BG + ($BGLeft-$BG) / (1.0 + exp(($E-$EBG0)/$WBG));
}

sub CBDOS
{
	my ($E, $D0CB, $ECBM) = @_;
	my $DE = 0.0;
	$DE = $D0CB * sqrt($E - $ECBM) if($E > $ECBM);
	return $DE;
}

sub VBDOS
{
	my ($E, $D0VB, $EVBM) = @_;
	my $DE = 0.0;
	$DE = $D0VB * sqrt($EVBM - $E) if($E < $EVBM);
	return $DE;
}

sub MetalDOS
{
	my ($E, $D0CB, $EF, $k0Metal) = @_;
	my $DE = 0.0;
#print "slop: [ECBM = $EF]: $k0Metal * ($E-$EF) = ", $k0Metal * ($E-$EF), "\n";
	if(defined $EF) {
		$DE = $D0CB + $k0Metal * ($E-$EF) if($E <= $EF);
	}
	else {
		$DE = $D0CB + $k0Metal * $E;
	}
	return $DE;
}

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 TotalDOS
{
	my ($E, $pVars, $IsMetal) = @_;

	my $y = 0.0;
	my $CB = 0.0;
	my $VB = 0.0;
	if($IsMetal) {
		$CB = MetalDOS($E, $pVars->{D0CB}, $pVars->{EF}, $pVars->{k0Metal}) if($pVars->{D0CB} != 0.0);
#		$CB = MetalDOS($E, $pVars->{D0CB}, undef, $pVars->{k0Metal}) if($pVars->{D0CB} != 0.0);
		$y += $CB;
	}
	else {
		$CB = CBDOS($E, $pVars->{D0CB}, $pVars->{ECBM}) if($pVars->{D0CB} != 0.0);
		$VB = VBDOS($E, $pVars->{D0VB}, $pVars->{EVBM}) if($pVars->{D0VB} != 0.0);
		$y += $CB + $VB;
	}
	my @GL;
	for(my $i = 0 ; $i < $pVars->{nGL} ; $i++) {
		if($pVars->{"GLC0$i"}) {
			$GL[$i] = Sci::GaussLorentz($E, $pVars->{"GLE0$i"}, $pVars->{"GLC0$i"},
						$pVars->{"GLWL$i"}, $pVars->{"GLGFraction$i"}, $pVars->{"GLGWRatio$i"});
		}
		else {
			$GL[$i] = 0.0;
		}
		$y += $GL[$i];
	}

	my $FD = FermiDiracFunction($E, $pVars->{T}, $pVars->{EF}, 1.0, 0.0);
	$y *= $FD;
	return ($y, $CB, $VB, \@GL, $FD);
}

sub EnergyCut
{
	my ($pX, $pY, $pVars) = @_;

	my $K0 = $pVars->{K0Ecut};
	return if($K0 <= 0.0);

	my $Ecut  = $pVars->{Ecut};
	my $BGcut = $pVars->{BGcut};
	my $C0    = $pVars->{C0Exp};
	my $W     = $pVars->{WExp};
	my $beta  = $pVars->{BetaExp};
	if($C0 > 0.0) {
		for(my $i = 0 ; $i < @$pY ; $i++) {
			my $x = $pX->[$i];
			$pY->[$i] += $C0 * exp(-(($x - $Ecut) / $W)**$beta);
		}
	}

	for(my $i = 0 ; $i < @$pY ; $i++) {
		my $x = $pX->[$i];
		if($x <= $Ecut) {
			$pY->[$i] = $BGcut;
		}
		else {
			my $y = $K0 * ($x - $Ecut);
			if($y < $pY->[$i]) {
				$pY->[$i] = $y;
			}
			else {
				last;
			}
		}
	}

	return $pY;
}

sub CalTotalDOSSpectrum
{
	my ($pX, $pVars, $IsMetal, $pObsX, $pObsY) = @_;

#print "n=$pVars->{nGL}: CL0: $pVars->{GLC00},$pVars->{GLE00},$pVars->{GLWL0},$pVars->{Ecut}\n";
	my $nData = @$pX;
	my @CalY;
	my (@CBFD, @VBFD, @CB, @VB, @GLArray, @FD, @BG);
	for(my $i = 0 ; $i < $pVars->{nGL} ; $i++) {
		$GLArray[$i] = [];
	}
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x = $pX->[$i];
		my $pGL;
		($CalY[$i], $CB[$i], $VB[$i], $pGL, $FD[$i]) = TotalDOS($x, $pVars, $IsMetal);
		$CBFD[$i] = $CB[$i] * $FD[$i];
		$VBFD[$i] = $VB[$i] * $FD[$i];
		for(my $j = 0 ; $j < $pVars->{nGL} ; $j++) {
			$GLArray[$j][$i] = $pGL->[$j];
		}
	}


#	if(!defined $pObsX or !defined $pObsY or $pVars->{BGSL0} == 0.0 or ($pVars->{EBGSL1} == $pVars->{EBGSL0})) {
	if($pVars->{BGSL0} == 0.0 or ($pVars->{EBGSL1} == $pVars->{EBGSL0})) {
##print "Use step BG\n";
		for(my $i = 0 ; $i < $nData ; $i++) {
			my $x = $pX->[$i];
			$BG[$i] = $pVars->{BG};
			for(my $j = 0 ; $j < $pVars->{nStepBG} ; $j++) {
				$BG[$i] += BG($x, 0.0, $pVars->{"BGLeft$j"}, $pVars->{"EBG0$j"}, $pVars->{"WBG$j"})
								 if($pVars->{"BGLeft$j"} != 0.0);
			}
			$CalY[$i] = $CalY[$i] + $BG[$i];
		}
	}
	else {
##print "Use Shirley BG\n";
		my $nXData = @$pObsX;
#my $nYData = @$pObsY;
#print "Obs: n=$nData, $nXData, $nYData; p=$pObsX, $pObsY\n";
		my $pInteg = Algorism::IntegrateByTrapezoid($pObsX, $pObsY);
		my $IE1    = Algorism::Interpolate($pObsX, $pInteg, $pVars->{EBGSL1});
		my $IE0    = Algorism::Interpolate($pObsX, $pInteg, $pVars->{EBGSL0});
		for(my $i = 0 ; $i < $nData ; $i++) {
			my $x = $pX->[$i];
			my $IE = Algorism::Interpolate($pObsX, $pInteg, $x);
			$BG[$i] = $pVars->{BG} + $pVars->{BGSL0} * ($IE - $IE1) / ($IE0 - $IE1);
		}

		my @NetY;
		for(my $i = 0 ; $i < $nXData ; $i++) {
			if(!defined $BG[$i]) {
				$NetY[$i] = $pObsY->[$i] - 0.0;
			}
			else {
				$NetY[$i] = $pObsY->[$i] - $BG[$i];
			}
		}
#print "p2=$pObsX\n";
#my $nXData = @$pObsX;
#my $nYData = @NetY;
#print "Net: n=$nXData, $nYData; p=$pObsX, $pObsY\n";
		$pInteg = Algorism::IntegrateByTrapezoid($pObsX, \@NetY);
		$IE1    = Algorism::Interpolate($pObsX, $pInteg, $pVars->{EBGSL1});
		$IE0    = Algorism::Interpolate($pObsX, $pInteg, $pVars->{EBGSL0});
		for(my $i = 0 ; $i < $nData ; $i++) {
			my $x = $pX->[$i];
			my $IE = Algorism::Interpolate($pObsX, $pInteg, $x);
			$BG[$i] = $pVars->{BG} + $pVars->{BGSL0} * ($IE - $IE1) / ($IE0 - $IE1);
		}

		for(my $i = 0 ; $i < $nData ; $i++) {
			$CalY[$i] = $CalY[$i] + $BG[$i];
		}
	}

	EnergyCut($pX, \@CalY, $pVars);

	my $pYConvoluted = Convolution($pX, \@CalY, $pVars->{AppFunction}, $pVars->{Wa});

	return ($pYConvoluted, \@CB, \@VB, \@CBFD, \@VBFD, \@GLArray, \@FD, \@BG); #\@CalY;
}

sub Convolution
{
	my ($pX, $pY, $AppFunction, $Wa) = @_;
	return $pY if($Wa <= 0.0);

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

	my $nData = @$pX;
	my $dX    = $pX->[1] - $pX->[0];
	my $absdX = abs($dX);
	my @Y;
	if($Wa <= 0.0) {
		for(my $i = 0 ; $i < $nData ; $i++) {
			$Y[$i] = $pY->[$i];
		}
		return \@Y;
	}

	for(my $i = 0 ; $i < $nData ; $i++) {
		$Y[$i] = 0.0;
	}

	my $nExt = int($Wa * 10.0 / $absdX + 0.01);
	$nExt = 10 if($nExt < 10);
	for(my $i = -$nExt ; $i <= $nData+$nExt ; $i++) {
		my $x0 = $pX->[0] + $dX * $i;
		my $y  = 0.0;
		if($i < 0) {
			$y = $pY->[0];
		}
		elsif($i >= $nData) {
			$y = $pY->[$nData-1];
		}
		else {
			$y = $pY->[$i];
		}
		for(my $j = 0 ; $j < $nData ; $j++) {
			my $xj  = $pX->[$j];
			my $hij = &$AppFunc($x0, $xj, $Wa) * $absdX;
#			$Y[$j] += $hij * $pY->[$j];
			$Y[$j] += $hij * $y;
		}
	}
	return \@Y;
}

sub Deconvolution
{
	my ($pXObs, $pYObs, $AppFunction, $Wa, $CorrectZero, $Dumping, $nMaxIter, $EPS) = @_;
	$CorrectZero = 1      if(!defined $CorrectZero);
	$Dumping     = 0.2    if(!defined $Dumping);
	$nMaxIter    = 30     if(!defined $nMaxIter);
	$EPS         = 1.0e-3 if(!defined $EPS);

	my $nData = @$pXObs;
	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 $absdX = abs($dX);
	my $hii = &$AppFunc(0.0, 0.0, $Wa) * $absdX;
	my $K = $Dumping / $hii;
print "hii=$hii K=$K\n";
	my $iter = 1;
	for($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) * $absdX;
					$S0 += $hij * $Y[$iter]->[$j];
				}
			}
			for(my $j = $i ; $j < $nData ; $j++) {
				my $xj = $pXObs->[$j];
				my $hij = &$AppFunc($xi, $xj, $Wa) * $absdX;
				$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);
		}
		my $diff =DeconvolutionDiff($Y[$iter-1], $Y[$iter]);
print "iter: $iter (sigma2=$diff)\n";
		if($diff < $EPS) {
print "Converged at iter=$iter\n";
			last;
		}
	}
	$iter = $nMaxIter-1 if($iter >= $nMaxIter);
	return ($Y[$iter], @Y);
}

sub DeconvolutionDiff
{
	my ($Y1, $Y2) = @_;
	
	my $Y2sum = 0.0;
	my $Y2diff = 0.0;
	for(my $i = 0 ; $i < @$Y1 ; $i++) {
		my $y1 = $Y1->[$i];
		my $y2 = $Y2->[$i];
		$Y2sum += $y1 * $y1;
		my $d = $y1 - $y2;
		$Y2diff += $d * $d;
	}
	$Y2diff /= $Y2sum if($Y2sum > 0.0);
	return $Y2diff;
}

1;