#!/usr/bin/perl

use strict;

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

use JFile;
use Sci qw($a0);
use Crystal::Atom;
use Crystal::Quantum;

my $a0InA = $a0 * 1.0e10;

#===============================================
# Parameters
#===============================================
my $OutFile = "HSForwardReverse.csv";

my $AtomName = 'B';
my $atom     = new Atom($AtomName);

my $n      =  1;
my $l      =  0;
my $m      =  0;
#my $Z      =  5.0;
my $Z      =  $atom->AtomicNumber();

my $nr     = 4001;
my $rstep  = 0.001;
my $r0     = 0.0;

my $E0        = -30;
my $RConnect  = 1.0;
#my $RConnect  = 1.5;
my $nApproximateForForwardSolver = 2;
my $nApproximateForReverseSolver = 2;

# For global search
my ($N0, $N1, $NStep) = (0.9, 5.0, 0.2);

# For Newton-Raphson method
my $dE       = 0.001;
my $nMaxIter = 100;
my $DF       = 0.1;
my $EEPS     = 1.0e-6;

my $GlobalSearch  = 1;
my $NewtonRaphson = 1;

my (@r, @Pf, @Pr);
my @Pcal;

my $E = -24.1052;
#my $S2 = $Quantum->CalPForwardReverse($l, $E, \@r, \@Pf, \@Pr, $Z, $nr, $r0, $rstep, $RConnect,
#				$nApproximateForForwardSolver, $nApproximateForReverseSolver);
#$Quantum->SaveP($OutFile, $l, $E, \@r, \@Pf, \@Pr, $Z, $nr, $r0, $rstep, $RConnect, @Pcal);
#exit;

#===============================================
# Main
#===============================================
my $Quantum = new Quantum;

#============================
# Global search
#============================
if($GlobalSearch) {
print "Global search from N = $N0 to $N1 at $NStep step\n";
my ($PrevE, $PrevS2);
my $N = $N0;
while(1) {
	$E = -$Z*$Z / $N/$N;

	my $S2 = $Quantum->CalPForwardReverse($l, $E, \@r, \@Pf, \@Pr,$Z, $nr, $r0, $rstep, $RConnect,
				$nApproximateForForwardSolver, $nApproximateForReverseSolver);
printf "  N = %6.2f  E = %6.4f: S2 = %g\n", $N, $E, $S2;
	if(defined $PrevS2) {
		if($PrevS2 * $S2 < 0.0) {
			$E0 = $PrevE + (0 - $PrevS2) / ($S2 - $PrevS2) * ($E - $PrevE);
			last;
		}
		else {
			$PrevE  = $E;
			$PrevS2 = $S2;
		}
	}

	$PrevE = $E;
	$PrevS2 = $S2;

	$N += $NStep;
	last if($N > $N1);
}

for(my $i = 0 ; $i < @r ; $i++) {
#	my $Rr = $Quantum->CalRadialWaveFunction2($Z, $n, $l, $m, $r[$i]);
#	$Pcal[$i] = $r[$i] * $Rr;
	my $Rr = $Quantum->CalRadialWaveFunction2($Z, $n, $l, $m, $r[$i] * $a0InA);
	$Pcal[$i] = $r[$i] * $a0InA**(1.5) * $Rr;
#print "$i: r=$r[$i]: P=$Pcal[$i]\n";
}

print "  Transition E = $E0\n";
my $S2 = $Quantum->CalPForwardReverse($l, $E0, \@r, \@Pf, \@Pr, $Z, $nr, $r0, $rstep, $RConnect,
			$nApproximateForForwardSolver, $nApproximateForReverseSolver);
$Quantum->SaveP($OutFile, $l, $E, \@r, \@Pf, \@Pr, $Z, $nr, $r0, $rstep, $RConnect, \@Pcal);
print "\n";
#exit;
}

#============================
# Newton-Raphson
#============================
if($NewtonRaphson) {
print "Newton-Rphson method from E = $E0 for max $nMaxIter cycles\n";
$E = $E0;
for(my $c = 1 ; $c <= $nMaxIter ; $c++) {
print "Newton cycle $c:\n";
	my $S2  = $Quantum->CalPForwardReverse($l, $E,     \@r, \@Pf, \@Pr,
				$Z, $nr, $r0, $rstep, $RConnect,
				$nApproximateForForwardSolver, $nApproximateForReverseSolver);
	my $S21 = $Quantum->CalPForwardReverse($l, $E+$dE, \@r, \@Pf, \@Pr,
				$Z, $nr, $r0, $rstep, $RConnect,
				$nApproximateForForwardSolver, $nApproximateForReverseSolver);
print "  Ei(E)    = $E: S2 = $S2\n";
print "  Ei(E+dE) = ", $E+$dE, ": S2 = $S21\n";
	my $dS2dE = ($S21 - $S2) / $dE;
print "    dS2/dE = $dS2dE\n";

	my $sign = ($dS2dE < 0.0)? -1.0 : 1.0;
	my $deltaE = - ($S2+$S21)/2.0 / ($sign * (abs($dS2dE) + $DF));
print "    E: $E => ";
	$E += $deltaE;
print "$E\n\n";

	if(abs($deltaE) < $EEPS) {
print "   **** Reaches EPS = $EEPS (> dE = $deltaE) at E = $E\n";
		last;
	}
}

$E0 = $E;
print "\n";
}


$Quantum->CalPForwardReverse($l, $E0, \@r, \@Pf, \@Pr,
			$Z, $nr, $r0, $rstep, $RConnect,
			$nApproximateForForwardSolver, $nApproximateForReverseSolver);
#$Quantum->Normalize($nr, \@Pr, $rstep);
my $f = $Quantum->CalRrNormalizationFactor($nr, \@Pr, $rstep);
for(my $i = 0 ; $i < $nr ; $i++) {
	$Pf[$i] /= $f;
	$Pr[$i] /= $f;
}

my ($S, $Scal) = (0.0, 0.0);
for(my $i = 0 ; $i < $nr ; $i++) {
	$S    += $Pr[$i]   * $Pr[$i]   * $rstep;
	$Scal += $Pcal[$i] * $Pcal[$i] * $rstep;
}
print "S      = $S\n";
print "S(cal) = $Scal\n";

$Quantum->SaveP($OutFile, $l, $E, \@r, \@Pf, \@Pr, $Z, $nr, $r0, $rstep, $RConnect, \@Pcal);

print "Converged E = $E0\n";

exit;

#===============================================
# Subroutines
#===============================================