#!/usr/bin/perl

use strict;

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

use JFile;
use Crystal::Atom;
use Crystal::Quantum;

use Sci::Algorism;
use Sci qw($pi);

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

my $AtomName = 'B';
my $atom     = new Atom($AtomName);
my $Z        =  $atom->AtomicNumber();

my $nr     = 3001;
my $rstep  = 0.01;
my $r0     = 0.0;

my $RConnect  = 1.0;
#my $RConnect  = 1.5;
my $nApproximateForForwardSolver = 5;
my $nApproximateForReverseSolver = 2;

my $nMaxIter = 100;
my $PEPS     = 1.0e-6;
my (@r, @Pf, @Pr);

my ($B0, $B1) = (10.0, 150.0);
my $B = 10.0;

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

#my $S2 = $Quantum->CalPForwardReverse($E, \@r, \@Pf, \@Pr, $nr, $r0, $rstep, $RConnect, $BForward, $BReverse);

$B = $B0;
$Quantum->CalTFPReverse(\@r, \@Pr, $nr, $r0, $rstep, $RConnect, $B, $nApproximateForReverseSolver);
my $S0 = $Pr[0] - 1.0;
$B = $B1;
$Quantum->CalTFPReverse(\@r, \@Pr, $nr, $r0, $rstep, $RConnect, $B, $nApproximateForReverseSolver);

my $S1 = $Pr[0] - 1.0;
for(my $i = 0 ; $i < $nMaxIter ; $i++) {
	$B = ($B0 + $B1) / 2.0;
	$Quantum->CalTFPReverse(\@r, \@Pr, $nr, $r0, $rstep, $RConnect, $B, $nApproximateForReverseSolver);
	my $Sm = $Pr[0] - 1.0;
printf "Cycle %2d: B=%8.4f - %8.4f - %8.4f, S=%12.4g - %12.4g - %12.4g\n", $i, $B0, $B, $B1, $S0, $Sm, $S1;
	if(abs($Sm) < $PEPS) {
		last;
	}
	if($S0*$Sm < 0.0) {
		$B1 = $B;
		$S1 = $Sm;
	}
	else {
		$B0 = $B;
		$S0 = $Sm;
	}
}
print "Converged at B = $B, P(0)=$Pr[0]\n";

$Quantum->SaveTFP($OutFile, \@r, \@Pf, \@Pr, $nr, $r0, $rstep, $RConnect);

exit;

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