#===============================================
# Quantum
#===============================================
package Quantum;
use Exporter;
use Common;
@ISA = qw(Exporter Common);

@EXPORT = qw();

use strict;
use File::Basename;
use Math::Complex;

use JFile;
use Sci qw($pi $a0 $h $hbar $me $kB $torad);
use Sci::Algorism;
use Crystal::Atom;

#============================================================
#@ÓIo[֐
#============================================================
my $a0InA = $a0 * 1.0e10;

#============================================================
#@ϐ擾֐ȂǍ
#============================================================

#============================================================
#@RXgN^AfXgN^
#============================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#============================================================
#@pNXŒ苠Ȃ֐
#============================================================
sub GetOrbName
{
	my ($this, $n, $l, $m) = @_;
	my @names1   = ('s'); 										# for l=0
	my @names2   = ('py', 'pz', 'px');							# for l=1, m=-1, 0, 1
	my @names3   = ('dxy', 'dyz', 'd3z2-r2', 'dzx', 'dx2-y2');	# for l=2, m=-2,-1,0,1,2
	my @names4   = ('fy(3x2-y2)', 'fxyz', 'fy(5z2-r2)', 'fz(5z2-3r2)',
		        'fx(5z2-r2)', 'fz(x2-y2)', 'fx(x2-3y2)');			# for l=3, m=-3,-2,-1,0,1,2,3
	my @names    = (\@names1, \@names2, \@names3, \@names4);
	my @alphabet = qw(s p d f g h i j k l m n);

	my $pNames = $names[$l];
	my $name   = ($pNames)? $pNames->[$m+$l] : $alphabet[$l];
	return "$n$name";
}

sub SpeculateElectronicConfiguration
{
	my ($this, $AtomicNumber, $Charge) = @_;

	my $nElectron = $AtomicNumber - $Charge;

	my $ResCharge = $nElectron;
	my $n1s = 0.0;
	if($ResCharge > 2.0) {
		$n1s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n1s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n2s = 0.0;
	if($ResCharge > 2.0) {
		$n2s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n2s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n2p = 0.0;
	if($ResCharge > 6.0) {
		$n2p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n2p = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n3s = 0.0;
	if($ResCharge > 2.0) {
		$n3s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n3s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n3p = 0.0;
	if($ResCharge > 6.0) {
		$n3p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n3p = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n3d = 0.0;
	my $n4s = 0.0;
	if($Charge == 0.0) {
		if($ResCharge > 2.0) {
			$n4s = 2.0;
			$ResCharge -= 2.0;
		}
		else {
			$n4s = $ResCharge;
			$ResCharge = 0.0;
		}
		if($ResCharge > 10.0) {
			$n3d = 10.0;
			$ResCharge -= 10.0;
		}
		else {
			$n3d = $ResCharge;
			$ResCharge = 0.0;
		}
	}
	else {
		if($ResCharge > 10.0) {
			$n3d = 10.0;
			$ResCharge -= 10.0;
		}
		else {
			$n3d = $ResCharge;
			$ResCharge = 0.0;
		}
		if($ResCharge > 2.0) {
			$n4s = 2.0;
			$ResCharge -= 2.0;
		}
		else {
			$n4s = $ResCharge;
			$ResCharge = 0.0;
		}
	}
	my $n4p = 0.0;
	if($ResCharge > 6.0) {
		$n4p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n4p = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n4d = 0.0;
	if($ResCharge > 10.0) {
		$n4d = 10.0;
		$ResCharge -= 10.0;
	}
	else {
		$n4d = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n5s = 0.0;
	if($ResCharge > 2.0) {
		$n5s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n5s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n5p = 0.0;
	if($ResCharge > 6.0) {
		$n5p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n5p = $ResCharge;
		$ResCharge = 0.0;
	}

	my $nTotal = $n1s + $n2s + $n2p + $n3s + $n3p + $n3d + $n4s + $n4p + $n4d + $n5s + $n5p;
	my @n = (
		[$n1s],
		[$n2s, $n2p],
		[$n3s, $n3p, $n3d],
		[$n4s, $n4p, $n4d],
		[$n5s, $n5p],
		);
	my %ret = (
		pnArray      => \@n,
		AtomicNumber => $AtomicNumber,
		nElectron    => $nElectron,
		ResCharge    => $ResCharge,
		nTotal  => $nTotal,
		n1s     => $n1s,
		n2s     => $n2s,
		n2p     => $n2p,
		n3s     => $n3s,
		n3p     => $n3p,
		n3d     => $n3d,
		n4s     => $n4s,
		n4p     => $n4p,
		n4d     => $n4d,
		n5s     => $n5s,
		n5p     => $n5p,
		);
	return \%ret;
}

sub CalSTO
{
	my ($this, $IonName) = @_;

	my ($AtomicNumber, $AtomName, $Charge) = AtomType::GetAtomInformation($IonName);
	my $pHash = $this->SpeculateElectronicConfiguration($AtomicNumber, $Charge);
	if($pHash->{ResCharge }> 0.0) {
		print("Error: Ne=$AtomicNumber - $Charge is not supported.\n");
		exit;
	}
	if($pHash->{nTotal} ne $pHash->{nElectron}) {
		print("Error: Inconsistent Ne (Ntotal=$pHash->{nTotal} != Ne=$pHash->{nElectron}).\n");
		exit;
	}

	my $nTotal    = $pHash->{nTotal};
	my $nElectron = $pHash->{nElectron};
	my $ResCharge = $pHash->{ResCharge};
	my $n1s       = $pHash->{n1s};
	my $n2s       = $pHash->{n2s};
	my $n2p       = $pHash->{n2p};
	my $n3s       = $pHash->{n3s};
	my $n3p       = $pHash->{n3p};
	my $n3d       = $pHash->{n3d};
	my $n4s       = $pHash->{n4s};
	my $n4p       = $pHash->{n4p};
	my $n4d       = $pHash->{n4d};
	my $n5s       = $pHash->{n5s};
	my $n5p       = $pHash->{n5p};

	my @ne = ($n1s, $n2s+$n2p, $n3s+$n3p+$n3d, $n4s+$n4p+$n4d, $n5s+$n5p);
	my $MaxN = 0.0;
	for(my $i = 0 ; $i < @ne ; $i++) {
		if($ne[$i] > 0.0) {
			$MaxN = $i+1;
		}
	}
	print("Max n: $MaxN\n");
	my @nStar = (undef, 1, 2, 3, 3.7, 4.0, 4.2);
	my $nStar = $nStar[$MaxN];
	print("n*   : $nStar\n");

	my $Is_d = 0;
	if(($n3d > 0.0 and $n4s == 0.0) or ($n4d > 0.0 and $n5s == 0.0)) {
		$Is_d = 1;
	}

#(i)  z*=z-s�Ƃ���B������z�͊j�d�ׁAs�͎Օ��萔
#(ii) �d�q��(1s)(2s,2p)(3s,3p)(3d)������Ƃ����Q�ɂ킯��B
#	my @nGroup = ($n1s, $n2s+$n2p, $n3s+$n3p, $n3d, $n4s+$n4p, $n4d, $n5s+$n5p);
	my @nGroup = ($n1s, $n2s+$n2p, $n3s+$n3p+$n3d, $n4s+$n4p+$n4d, $n5s+$n5p);
	my $nMax = 0;
	for(my $i = 0 ; $i < @nGroup ; $i++) {
		if($nGroup[$i] > 0.0) {
			$nMax = $i;
		}
	}
	print("n(max) for groups: $nMax\n");

	my $stotal = 0.0;
#(iii) �l���Ă���d�q����ʎq��n�̑傫���d�q�����s�ւ̊�^�̓[���Ƃ���B
#(iv)  �l���Ă���d�q�̑�����Q�̑��̓d�q����͂��ꂼ��0.35�Ƃ���B�������A1s�Q�̏ꍇ��0.30�Ƃ���B
	my $s;
	my $n = $nGroup[$nMax]-1.0;
	my $k;
	if($nGroup[$nMax] > 1.0) {
		if($nMax == 0) { # 1sQ
			$k = 0.30;
		}
		else {
			$k = 0.35;
		}
		$s = $n * $k;
	}
print "s(outmost)  = $n * $k = $s\n";
	$stotal += $s;

#(v) �l���Ă���d�q��s�܂���p�d�q�Ȃ�A������n�̈�������d�q�����0.85����
#�������Ad�d�q�Ȃ�΁A����������̓d�q����̊�^�͂��ׂ�1.00�Ƃ���B
	$n = $nGroup[$nMax-1];
	if($nMax >= 1 and $nGroup[$nMax-1] > 1.0) {
#		if($nMax == 3) { # dQ
		if($Is_d) { # dQ
			$k = 1.0;
		}
		else {
			$k = 0.85;
		}
		$s = $n * $k;
	}
print "s(n(max)-1) = $n * $k = $s\n";
	$stotal += $s;

#�����Ɠ����̓d�q�����1.00���Š�^����B�������Ad�d�q�Ȃ�΁A����������̓d�q����̊�^�͂��ׂ�1.00�Ƃ���B
	$k = 1.0;
	$n = 0.0;
	for(my $i = $nMax - 2 ; $i >= 0 ; $i--) {
			$n += $nGroup[$i] * 1.0;
	}
	$s = $n * $k;
	$stotal += $s;
print "s(inner)    = $n * $k = $s\n";
print "s(total)    = $stotal\n";

	my $zStar = $AtomicNumber - $stotal;
	my $Rmax  = $nStar*$nStar / $zStar * $a0InA;
	my $n1    = $nStar - 1;

	return ($zStar, $nStar, $n1, $Rmax);
}

sub HPG
{
	my ($this, $n, $x, $x0, $b) = @_;

	my $s = ($x - $x0) / $b;
	my $G = &Gaussian($s, 0.0, 1.0);
	my $H = Sci::HermitePolynomial2($n, $s);
	my $k = sqrt(2.0**$n * sqrt($pi) * Sci::Factorial($n));
	return $k * $H * $G;
}

my $kPhi = sqrt(1.0 / $pi);
sub CalPhi
{
	my ($this, $m, $phi) = @_; # phi in rad
	return 1.0 if($m == 0);
	return $kPhi * cos($m * $phi) if($m > 0);
	return $kPhi * sin($m * $phi);
}

my $kPhi2 = sqrt(0.5 / $pi);
sub CalPhi2
{
	my ($this, $m, $phi) = @_; # phi in rad
	return $kPhi2 * exp(i * $m * $phi);
}

sub CalTheta
{
	my ($this, $l, $m, $Q) = @_; # Q in rad

	return sqrt(1.0/2.0) 					if($l == 0);			# s
	my $cosQ = cos($Q);
	my $sinQ = sin($Q);
	return  sqrt(3.0/4.0)    * $cosQ 			if($l == 1 and $m ==  0);	# pz
	return  sqrt(3.0/4.0)    * $sinQ 			if($l == 1 and $m == -1);	# px
	return -sqrt(3.0/4.0)    * $sinQ 			if($l == 1 and $m == +1);	# py
	return  sqrt(15.0/16.0)  * $sinQ*$sinQ 			if($l == 2 and $m == +2);	# dx2-y2
	return -sqrt(15.0/4.0)   * $sinQ*$cosQ 			if($l == 2 and $m == +1);	# dzx
	return  sqrt(5.0/16.0)   * (3.0*$cosQ*$cosQ-1.0) 	if($l == 2 and $m ==  0);	# dz2
	return  sqrt(15.0/4.0)   * $sinQ*$cosQ 			if($l == 2 and $m == -1);	# dyz
	return  sqrt(15.0/4.0)   * $sinQ*$sinQ 			if($l == 2 and $m == -2);	# dxy
	return -sqrt(35.0/32.0)  * $sinQ*$sinQ*$sinQ 		if($l == 3 and $m == +3);	# fx(x2-3y2)
	return  sqrt(105.0/16.0) * $sinQ*$sinQ*$cosQ 		if($l == 3 and $m == +2);	# fz(x2-y2)
	return -sqrt(21.0/32.0)  * $sinQ*(5.0*$cosQ*$cosQ-1.0) 	if($l == 3 and $m == +1);	# fx(5z2-r2)
	return  sqrt(7.0/16.0)   * $cosQ*(5.0*$cosQ*$cosQ-3.0) 	if($l == 3 and $m ==  0);	# fz(5z2-3r2)
	return  sqrt(21.0/31.0)  * $sinQ*(5.0*$cosQ*$cosQ-1.0) 	if($l == 3 and $m == -1);	# fy(5z2-r2)
	return  sqrt(105.0/1.0)  * $sinQ*$sinQ*$cosQ 		if($l == 3 and $m == -2);	# fxyz
	return  sqrt(35.0/8.0)   * $sinQ*$sinQ*$sinQ 		if($l == 3 and $m == -3);	# fy(3x2-y2)
}

sub CalTheta2
{
	my ($this, $l, $m, $Q) = @_; # Q in rad
	
	my $ma   = ($m < 0)? -$m : $m;
	my $sign = ($m <= 0)? 1 : (-1)**$m;
	my $lmm    = Sci::Factorial($l - $ma);
	my $lpm    = Sci::Factorial($l + $ma);
	my $k    = sqrt( ($l+0.5) * $lmm / $lpm);

	return $sign * $k * Sci::LegendrePolynomial2($l, $ma, cos($Q));
}

# Normalized with respect to r from 0 to inf for r^2*Rr^2
sub CalRadialWaveFunction2
{
	my ($this, $Z, $n, $l, $m, $r) = @_;

	my $rho    = 2.0 * $r / $a0InA * $Z / $n;
	my $L      = Sci::LaguerrePolynomial2($n+$l, 2.0*$l+1, $rho);
	my $nl1    = Sci::Factorial($n - $l - 1);
	my $nl     = Sci::Factorial($n + $l);
	my $c      = -(2.0 * $Z / $n / $a0InA)**1.5 * sqrt($nl1 / 2.0 / $n / $nl);
	my $Rr     = $c * $rho**$l * $L * exp(-$rho / 2.0);
	return $Rr / $nl;
}

# Normalized with respect to r from 0 to inf for r^2*Rr^2
sub CalRadialWaveFunction
{
	my ($this, $Z, $n, $l, $m, $r) = @_;

	my $Z32  = ($Z / $n / $a0InA)**1.5;
	my $rho  = 2.0 * $Z / $n / $a0InA * $r;

	my $Rr;
	if($n == 1 and $l == 0) {	# 1s
		$Rr = 2.0 * $Z32 * exp(-$rho/2.0);
#		$Rr = $Z32 * 2.0 * exp(-$rho/2.0);
#k=1
	}
	elsif($n == 2 and $l == 0) {	# 2s
		$Rr = 2.0 * $Z32 * (2.0-$rho/2.0) * exp(-$rho/2.0);
#		$Rr = 0.5 * sqrt(2.0) * $Z32 * (2.0-$rho) * exp(-$rho/2.0);
#k=4/sqrt(2)
	}
	elsif($n == 2 and $l == 1) {	# 2p
		$Rr = 1.0/sqrt(3.0) * $Z32 * $rho * exp(-$rho/2.0);
#		$Rr = 0.5 * sqrt(6.0) * $Z32 * $rho * exp(-$rho/2.0);
#k=2/sqrt(2)
	}
	elsif($n == 3 and $l == 0) {	# 3s
		$Rr = 2.0/27.0 * $Z32 * (27.0-27.0*$rho+9.0/2.0*$rho*$rho) * exp(-$rho/2.0);
#		$Rr = 1.0/9.0 * sqrt(3.0) * $Z32 * (6.0-6.0*$rho+$rho*$rho) * exp(-$rho/2.0);
#k=2/3/sqrt(3)
	}
	elsif($n == 3 and $l == 1) {	# 3p
		$Rr = 1.0/9.0*sqrt(2.0) * $Z32 * (6.0-1.5*$rho)*$rho * exp(-$rho/2.0);
#		$Rr = 1.0/9.0 * sqrt(6.0) * $Z32 * (4.0-$rho)*$rho * exp(-$rho/2.0);
#k=2/3/sqrt(3)
	}
	elsif($n == 3 and $l == 2) {	# 3d
		$Rr = 1.0/6.0*sqrt(2.0/5.0) * $Z32 * $rho*$rho * exp(-$rho/2.0);
#		$Rr = 1.0/9.0 * sqrt(30.0) * $Z32 * $rho*$rho * exp(-$rho/2.0);
#k=9/6*sqrt(2/5/30)=3/2/5/sqrt(3)
	}
	elsif($n == 4 and $l == 0) {	# 4s
		$Rr = 8.0 * 1.0/96.0 * $Z32 * (24.0-36.0*$rho+12.0*$rho*$rho-$rho*$rho*$rho) * exp(-$rho/2.0);
#		$Rr = 1.0/96.0 * $Z32 * (24.0-36.0*$rho+12.0*$rho*$rho-$rho*$rho*$rho) * exp(-$rho/2.0);
	}
	elsif($n == 4 and $l == 1) {	# 4p
		$Rr = 0.53333 * 1.0/32.0 * sqrt(15.0) * $Z32 * (20.0-10.0*$rho+$rho*$rho)*$rho * exp(-$rho/2.0);
#		$Rr = 1.0/32.0 * sqrt(15.0) * $Z32 * (20.0-10.0*$rho+$rho*$rho)*$rho * exp(-$rho/2.0);
	}
	elsif($n == 4 and $l == 2) {	# 4d
		$Rr = 1.6 * 1.0/96.0 * sqrt(5.0) * $Z32 * (6.0-$rho)*$rho*$rho * exp(-$rho/2.0);
#		$Rr = 1.0/96.0 * sqrt(5.0) * $Z32 * (6.0-$rho)*$rho*$rho * exp(-$rho/2.0);
	}
	elsif($n == 4 and $l == 3) {	# 4f
		$Rr = 0.228571 * 1.0/96.0 * sqrt(35.0) * $Z32 * $rho*$rho*$rho * exp(-$rho/2.0);
#		$Rr = 1.0/96.0 * sqrt(35.0) * $Z32 * $rho*$rho*$rho * exp(-$rho/2.0);
	}
	else {
		return undef;
	}
#print "x=$x  r=$r rho=$rho R=$Rr\n";
	return $Rr;
}

sub CalTFPForwardReverse
{
	my ($this, $pr, $pPf, $pPr, $nr, $r0, $rstep, $RConnect, $BForward, $BReverse) = @_;
	my $irc = ($RConnect - $r0) / $rstep;

	$this->CalTFPForward($pr, $pPf, $BForward);
	$this->CalTFPReverse($pr, $pPr, $BReverse);

	my $dPdrf = ($pPf->[$irc+1] - $pPf->[$irc-1]) / 2.0 / $rstep;
	my $dPdrr = ($pPr->[$irc+1] - $pPr->[$irc-1]) / 2.0 / $rstep;
#print "Forward (rc=$RConnect): P = rR = $Pf[$irc]  dP/dF = $dPdrf\n";
#print "Reverse (rc=$RConnect): P = rR = $Pr[$irc]  dP/dF = $dPdrr\n";

#	my $Pc = ($pPf->[$irc] > $pPr->[$irc])? $pPr->[$irc] : $pPf->[$irc];
	my $Pc = $pPf->[$irc];
	my $kf  = $Pc / $pPf->[$irc];
	my $kr  = $Pc / $pPr->[$irc];
#print "\nAdjust to P = $Pc at r = $RConnect\n";

	for(my $i = 0 ; $i < $nr ; $i++) {
		$pPf->[$i] *= $kf;
		$pPr->[$i] *= $kr;
	}
	my $dPdrf = ($pPf->[$irc+1] - $pPr->[$irc-1]) / 2.0 / $rstep;
	my $dPdrr = ($pPr->[$irc+1] - $pPr->[$irc-1]) / 2.0 / $rstep;
#print "Forward (rc=$RConnect): P = rR = $Pf[$irc]  dP/dF = $dPdrf\n";
#print "Reverse (rc=$RConnect): P = rR = $Pr[$irc]  dP/dF = $dPdrr\n";

	my $S2 = ($dPdrf - $dPdrr) / (abs($dPdrf) + abs($dPdrf)) * 2.0;
#	my $S2 = (($dPdrf - $dPdrr) / $Pc)**2;
#print "  S2 = $S2 at E = $E\n";

	return $S2;
}

sub CalTFPForward
{
	my ($this, $pr, $pP, $nr, $r0, $rstep, $RConnect, $B, $nApproximateForForwardSolver) = @_;
	my $irc = ($RConnect - $r0) / $rstep;

	for(my $i = 0 ; $i < $nr ; $i++) {
		$pr->[$i] = $r0 + $i * $rstep;
	}

	for(my $i = 0 ; $i <= $nApproximateForForwardSolver ; $i++) {
		$pP->[$i] = $this->ApproximateTFRZero($pr->[$i], undef, 1.0);
#print "f$i r=$pr->[$i] P=", $pP->[$i], "\n";
	}
	for(my $i = $nApproximateForForwardSolver ; $i < $nr-1 ; $i++) {
		my $r = $pr->[$i];
		$pP->[$i+1] = 2.0*$pP->[$i] - $pP->[$i-1] + $rstep*$rstep * $pP->[$i]**1.5 / sqrt($r);
		$pP->[$i] = 0.0 if($pP->[$i] < 0.0);
#print "f$i r=$r P=", $pP->[$i+1], "\n";
	}
#	my $nNode = &Normalize($nr, $pP, $rstep, "f");

	return;
}

sub CalTFPReverse
{
	my ($this, $pr, $pP, $nr, $r0, $rstep, $RConnect, $B, $nApproximateForReverseSolver) = @_;
	my $irc = ($RConnect - $r0) / $rstep;

	for(my $i = 0 ; $i < $nr ; $i++) {
		$pr->[$i] = $r0 + $i * $rstep;
	}

	for(my $i = 1 ; $i <= $nApproximateForReverseSolver ; $i++) {
		$pP->[$nr-$i] = $this->ApproximateTFRInf($pr->[$nr - $i], $B, 1.0);
#print "r$i r=", $pr->[$nr-$i], " P=", $pP->[$nr-$i], "\n";
	}
	for(my $i = $nr-$nApproximateForReverseSolver ; $i >= 1 ; $i--) {
		my $r = $pr->[$i];
		$pP->[$i-1] = 2.0*$pP->[$i] - $pP->[$i+1] + $rstep*$rstep * $pP->[$i]**1.5 / sqrt($r);
		$pP->[$i-1] = 0.0 if($pP->[$i-1] < 0.0);
#print "r$i r=$r P=", $pP->[$i+1], "\n";
	}

#	my $nNode = &Normalize($nr, $pP, $rstep, "r");

	return;
}

sub ApproximateTFRZero
{
	my ($this, $r, $B, $k) = @_;
	$B = -1.588076779 if(!defined $B);
	my $f = 1.0 + $B * $r + 4.0/3.0 * $r**1.5 + 0.4 * $B * $r**2.5 + 1.0/3.0 * $r*$r*$r;
	return $k * $f;
}

sub ApproximateTFRInf
{
	my ($this, $r, $B, $k) = @_;
	return undef if($r <= 0.0);
	$B = 140.0 if(!defined $B);
	my $f = $B / ($r*$r*$r);
	return $k*$f;

	my $lambda = 3.886;
	my $r0     = 5.2415;
	my $f = 1.0 / (1.0 + ($r/$r0)**(3.0/$lambda))**$lambda;
	return $k * $f;
}

sub SaveTFP
{
	my ($this, $OutFile, $pr, $pPf, $pPr, $nr, $r0, $rstep, $RConnect) = @_;
	my $irc = ($RConnect - $r0) / $rstep;

	my @dPdr;
	for(my $i = 0 ; $i < $nr ; $i++) {
		if($i < $irc) {
			$dPdr[$i] = ($pPf->[$i+1] - $pPf->[$i]) / $rstep;
		}
		else {
			$dPdr[$i] = ($pPr->[$i+1] - $pPr->[$i]) / $rstep;
		}
	}

	my @rl1;
	my $kf = $pPf->[1] / $this->ApproximateTFRZero($pr->[1], undef, 1.0);
	for(my $i = 0 ; $i < $nr ; $i++) {
		$rl1[$i] = &ApproximateRZero($pr->[$i], undef, $kf);
	}
	my @expE;
	my $kr = $pPr->[$nr-1] / $this->ApproximateTFRInf($pr->[$nr-1], undef, 1.0);
	for(my $i = 0 ; $i < $nr ; $i++) {
		$expE[$i] = &ApproximateRInf($pr->[$i], undef, $kr);
	}

	my $out = JFile->new($OutFile, "w") or die "$!: Can not write to [$OutFile].\n";
	$out->print("r,Pf,Pr,dP/dr,r**(l+1),exp(-sqrt(E)r)\n");
	for(my $i = 0 ; $i < @$pr ; $i++) {
		$out->print("$pr->[$i],$pPf->[$i],$pPr->[$i],$dPdr[$i],$rl1[$i],$expE[$i]\n");
	}
	$out->Close();
}

sub CalPForwardReverse
{
	my ($this, $l, $E, $pr, $pPf, $pPr, $Z, $nr, $r0, $rstep, $RConnect,
		$nApproximateForForwardSolver, $nApproximateForReverseSolver) = @_;

	my $irc = ($RConnect - $r0) / $rstep;

	$this->CalPForward($l, $E, $pr, $pPf, $Z, $nr, $r0, $rstep, $nApproximateForForwardSolver);
	$this->CalPReverse($l, $E, $pr, $pPr, $Z, $nr, $r0, $rstep, $nApproximateForReverseSolver);

	my $dPdrf = ($pPf->[$irc+1] - $pPf->[$irc-1]) / 2.0 / $rstep;
	my $dPdrr = ($pPr->[$irc+1] - $pPr->[$irc-1]) / 2.0 / $rstep;
#print "Forward (rc=$RConnect): P = rR = $Pf[$irc]  dP/dF = $dPdrf\n";
#print "Reverse (rc=$RConnect): P = rR = $Pr[$irc]  dP/dF = $dPdrr\n";

#	my $Pc = ($pPf->[$irc] > $pPr->[$irc])? $pPr->[$irc] : $pPf->[$irc];
	my $Pc = $pPf->[$irc];
#	my $Pc = abs($pPr->[$irc]);
	my $kf = $Pc / $pPf->[$irc];
	my $kr = $Pc / $pPr->[$irc];
#print "\nAdjust to P = $Pc at r = $RConnect\n";

	for(my $i = 0 ; $i < $nr ; $i++) {
		$pPf->[$i] *= $kf;
		$pPr->[$i] *= $kr;
	}
	my $dPdrf = ($pPf->[$irc+1] - $pPr->[$irc-1]) / 2.0 / $rstep;
	my $dPdrr = ($pPr->[$irc+1] - $pPr->[$irc-1]) / 2.0 / $rstep;
#print "Forward (rc=$RConnect): P = rR = $Pf[$irc]  dP/dF = $dPdrf\n";
#print "Reverse (rc=$RConnect): P = rR = $Pr[$irc]  dP/dF = $dPdrr\n";

	my $S2 = ($dPdrf - $dPdrr) / (abs($dPdrf) + abs($dPdrf)) * 2.0;
#	my $S2 = (($dPdrf - $dPdrr) / $Pc)**2;
#print "  S2 = $S2 at E = $E\n";

	return $S2;
}

sub CalPForward
{
	my ($this, $l, $E, $pr, $pP, $Z, $nr, $r0, $rstep, $nApproximateForForwardSolver) = @_;

	for(my $i = 0 ; $i < $nr ; $i++) {
		$pr->[$i] = $r0 + $i * $rstep;
	}

	for(my $i = 0 ; $i <= $nApproximateForForwardSolver ; $i++) {
		$pP->[$i] = $this->ApproximateRZero($pr->[$i], $l, $E, 1.0);
	}

	for(my $i = $nApproximateForForwardSolver ; $i < $nr-1 ; $i++) {
		my $r = $pr->[$i];
		my $g = -$E - 2.0 * $Z / $r + $l*($l+1) / $r / $r;
		$pP->[$i+1] = (2.0 + $rstep*$rstep * $g) * $pP->[$i] - $pP->[$i-1];
	}

	my $nNode = $this->Normalize($nr, $pP, $rstep);

	return ($nNode);
}


sub CalPReverse
{
	my ($this, $l, $E, $pr, $pP, $Z, $nr, $r0, $rstep, $nApproximateForReverseSolver) = @_;

	for(my $i = 0 ; $i < $nr ; $i++) {
		$pr->[$i] = $r0 + $i * $rstep;
	}

	for(my $i = 1 ; $i <= $nApproximateForReverseSolver ; $i++) {
		$pP->[$nr-$i] = $this->ApproximateRInf($pr->[$nr - $i], $l, $E, 1.0);
	}
	for(my $i = $nr-$nApproximateForReverseSolver ; $i >= 1 ; $i--) {
		my $i1 = $i + 1;
		my $r = $pr->[$i];
		my $g = -$E - 2.0 * $Z / $r + $l*($l+1) / $r / $r;
		$pP->[$i-1] = (2.0 + $rstep*$rstep * $g) * $pP->[$i] - $pP->[$i+1];
	}

	my $nNode = $this->Normalize($nr, $pP, $rstep);

	return ($nNode);
}

sub CalRrNormalizationFactor
{
	my ($this, $nr, $pP, $rstep) = @_;

	my $f = 0.0;
	for(my $i = 0 ; $i < $nr ; $i++) {
		$f += $pP->[$i]*$pP->[$i] * $rstep;
	}
	$f = sqrt($f);
#print "f = $f\n";

	return $f;
}

sub Normalize
{
	my ($this, $nr, $pP, $rstep) = @_;

	my $f = $this->CalRrNormalizationFactor($nr, $pP, $rstep);
	my $nNode = 0;
	for(my $i = 0 ; $i < $nr ; $i++) {
		$pP->[$i] /= $f;
		if($i > 0 and $pP->[$i] * $pP->[$i-1] < 0.0) {
			$nNode++;
		}
	}
	return $nNode;
}

sub ApproximateRZero
{
	my ($this, $r, $l, $E, $k) = @_;
	return $k * $r**($l+1);
}

sub ApproximateRInf
{
	my ($this, $r, $l, $E, $k) = @_;
	return $k * $r * exp(-$r * sqrt(-$E));
}

sub SaveP
{
	my ($this, $OutFile, $l, $E, $pr, $pPf, $pPr, $Z, $nr, $r0, $rstep, $RConnect, $pPcal) = @_;
	my $irc = ($RConnect - $r0) / $rstep;
	$pPcal = () if(!defined $pPcal);

	my @dPdr;
	for(my $i = 0 ; $i < $nr ; $i++) {
		if($i < $irc) {
			$dPdr[$i] = ($pPf->[$i+1] - $pPf->[$i]) / $rstep;
		}
		else {
			$dPdr[$i] = ($pPr->[$i+1] - $pPr->[$i]) / $rstep;
		}
	}

	my @rl1;
	my $kf = $pPf->[1] / $this->ApproximateRZero($pr->[1], $l, $E, 1.0);
	for(my $i = 0 ; $i < $nr ; $i++) {
		$rl1[$i] = $this->ApproximateRZero($pr->[$i], $l, $E, $kf);
	}
	my @expE;
	my $kr = $pPr->[$nr-1] / $this->ApproximateRInf($pr->[$nr-1], $l, $E, 1.0);
	for(my $i = 0 ; $i < $nr ; $i++) {
		$expE[$i] = $this->ApproximateRInf($pr->[$i], $l, $E, $kr);
	}

	my $out = JFile->new($OutFile, "w") or die "$!: Can not write to [$OutFile].\n";
	$out->print("r,Pf,Pr,dP/dr,r**(l+1),rexp(-sqrt(E)r),P(cal)\n");
	for(my $i = 0 ; $i < @$pr ; $i++) {
		$out->print("$pr->[$i],$pPf->[$i],$pPr->[$i],$dPdr[$i],$rl1[$i],$expE[$i],$pPcal->[$i]\n");
	}
	$out->Close();
}

1;