#!/usr/bin/perl

BEGIN {
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "c:/Programs/Perl/lib", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;

use Sci qw($pi $e $e0 erfc erf);
use Sci::Algorism;

#===============================
# Fundamental parameters
#===============================
my $KC = $e / 4.0 / $pi / $e0 * 1.0e10;

#===============================
# File configuration
#===============================
my $DiffFile = 'diff.csv';

# Maximum connection radius allowed
my $rcmax = 1.4;

# R mesh for output
my $N = 1500;
my ($rmin, $rstep) = (0.1, 0.01);
#my ($rmin, $rstep) = (0.01, 0.01);
my $rmax = $rmin + ($N - 1) * $rstep;

# Numrical differentiation
my $h          = 0.0001; # delta r
my $nDiffOrder = 7;      # 2, 3, 5, 7, -1=Analytical

# Binary search criteria
my $nMaxIter = 100;
my $f2EPS = 1.0e-6;

#===============================
# Potential definition
#===============================
my $type = 'buckingham';

# Atoms to be output
my @atoms = ('O', 'Zn', 'Ga', 'In');
my $atom1 = ($ARGV[0])? $ARGV[0] : $atoms[0];
my $atom2 = ($ARGV[1])? $ARGV[1] : $atoms[0];

my $PotentialFile = "Potential-$atom1-$atom2.table";


# Coulomb potential parameters

# $SplitEwaldSum: 1にするとPotential.csvにはCoulomb項を入れず、
#                 LAMMPS .inにewald項を入れる ($fCoulomb = 1)
my $SplitEwaldSum = 1;
my $alpha = 0.3;

# $fCoulomb: 最終的なポテンシャル中のCoulomb項の比率
my $fCoulomb = 1.0;

my %Zi;
$Zi{In} =  3.0; #1.8; #3.0;
$Zi{Ga} =  3.0; #1.8; #3.0;
$Zi{Zn} =  2.0; #1.2; #2.0;
$Zi{O}  = -2.0; #-1.2; #-2.0;

# Born repulsion parameters
my (%Aij, %rhoij);

$Aij{In}{In} = $Aij{In}{In} =    0.0;
$Aij{In}{Ga} = $Aij{Ga}{In} =    0.0;
$Aij{In}{Zn} = $Aij{Zn}{In} =    0.0;
$Aij{In}{O}  = $Aij{O}{In}  = 1293.600;
$Aij{Ga}{Ga} = $Aij{Ga}{Ga} =    0.0;
$Aij{Ga}{Zn} = $Aij{Zn}{Ga} =    0.0;
$Aij{Ga}{O}  = $Aij{O}{Ga}  = 2339.766;
$Aij{Zn}{Zn} = $Aij{Zn}{Zn} =    0.0;
$Aij{Zn}{O}  = $Aij{O}{Zn}  =  600.300;
$Aij{O}{O}   = $Aij{O}{O}   =   25.410;

$rhoij{In}{In} = $rhoij{In}{In} = 1.0;
$rhoij{In}{Ga} = $rhoij{Ga}{In} = 1.0;
$rhoij{In}{Zn} = $rhoij{Zn}{In} = 1.0;
$rhoij{In}{O}  = $rhoij{O}{In}  = 0.3310;
$rhoij{Ga}{Ga} = $rhoij{Ga}{Ga} = 1.0;
$rhoij{Ga}{Zn} = $rhoij{Zn}{Ga} = 1.0;
$rhoij{Ga}{O}  = $rhoij{O}{Ga}  = 0.2740;
$rhoij{Zn}{Zn} = $rhoij{Zn}{Zn} = 1.0;
$rhoij{Zn}{O}  = $rhoij{O}{Zn}  = 0.3370;
$rhoij{O}{O}   = $rhoij{O}{O}   = 0.6940;

# Dispersion parameters
my %Cij;

$Cij{In}{In} = $Cij{In}{In} = 0.0;
$Cij{In}{Ga} = $Cij{Ga}{In} = 0.0;
$Cij{In}{Zn} = $Cij{Zn}{In} = 0.0;
$Cij{In}{O}  = $Cij{O}{In}  = 4.325;
$Cij{Ga}{Ga} = $Cij{Ga}{Ga} = 0.0;
$Cij{Ga}{Zn} = $Cij{Zn}{Ga} = 0.0;
$Cij{Ga}{O}  = $Cij{O}{Ga}  = 0.0;
$Cij{Zn}{Zn} = $Cij{Zn}{Zn} = 0.0;
$Cij{Zn}{O}  = $Cij{O}{Zn}  = 0.0;
$Cij{O}{O}   = $Cij{O}{O}   = 32.32;

# Morse potential parameters
my %CR12;
$CR12{In}{In} = 0.0;

my (%Deij, %a0ij, %r0ij);

$Deij{In}{In} = $Deij{In}{In} = 0.0;
$Deij{In}{Ga} = $Deij{Ga}{In} = 0.0;
$Deij{In}{Zn} = $Deij{Zn}{In} = 0.0;
$Deij{In}{O}  = $Deij{O}{In}  = 0.072975;
$Deij{Ga}{Ga} = $Deij{Ga}{Ga} = 0.0;
$Deij{Ga}{Zn} = $Deij{Zn}{Ga} = 0.0;
$Deij{Ga}{O}  = $Deij{O}{Ga}  = 0.038263;
$Deij{Zn}{Zn} = $Deij{Zn}{Zn} = 0.0;
$Deij{Zn}{O}  = $Deij{O}{Zn}  = 0.001221;
$Deij{O}{O}   = $Deij{O}{O}   = 0.042395;

$a0ij{In}{In} = $a0ij{In}{In} = 0.0;
$a0ij{In}{Ga} = $a0ij{Ga}{In} = 0.0;
$a0ij{In}{Zn} = $a0ij{Zn}{In} = 0.0;
$a0ij{In}{O}  = $a0ij{O}{In}  = 1.9342;
$a0ij{Ga}{Ga} = $a0ij{Ga}{Ga} = 0.0;
$a0ij{Ga}{Zn} = $a0ij{Zn}{Ga} = 0.0;
$a0ij{Ga}{O}  = $a0ij{O}{Ga}  = 2.1052;
$a0ij{Zn}{Zn} = $a0ij{Zn}{Zn} = 0.0;
$a0ij{Zn}{O}  = $a0ij{O}{Zn}  = 2.347100;
$a0ij{O}{O}   = $a0ij{O}{O}   = 1.379316;

$r0ij{In}{In} = $r0ij{In}{In} = 0.0;
$r0ij{In}{Ga} = $r0ij{Ga}{In} = 0.0;
$r0ij{In}{Zn} = $r0ij{Zn}{In} = 0.0;
$r0ij{In}{O}  = $r0ij{O}{In}  = 2.82499;
$r0ij{Ga}{Ga} = $r0ij{Ga}{Ga} = 0.0;
$r0ij{Ga}{Zn} = $r0ij{Zn}{Ga} = 0.0;
$r0ij{Ga}{O}  = $r0ij{O}{Ga}  = 2.668100;
$r0ij{Zn}{Zn} = $r0ij{Zn}{Zn} = 0.0;
$r0ij{Zn}{O}  = $r0ij{O}{Zn}  = 3.167030;
$r0ij{O}{O}   = $r0ij{O}{O}   = 3.618701;

# 1/r^n potential parameters
my $nRn = 12;
my %CRnij;

$CRnij{In}{In} = $CRnij{In}{In} = 0.0;
$CRnij{In}{Ga} = $CRnij{Ga}{In} = 0.0;
$CRnij{In}{Zn} = $CRnij{Zn}{In} = 0.0;
$CRnij{In}{O}  = $CRnij{O}{In}  = 0.0;
$CRnij{Ga}{Ga} = $CRnij{Ga}{Ga} = 0.0;
$CRnij{Ga}{Zn} = $CRnij{Zn}{Ga} = 0.0;
$CRnij{Ga}{O}  = $CRnij{O}{Ga}  = 0.0;
$CRnij{Zn}{Zn} = $CRnij{Zn}{Zn} = 0.0;
$CRnij{Zn}{O}  = $CRnij{O}{Zn}  = 0.0;
$CRnij{O}{O}   = $CRnij{O}{O}   = 0.0;


#===============================
# Main routine
#===============================

print "\n";
print "Make LAMMPS connected potential [$DiffFile] for $type potential\n";
print "\n";
print "r range: $rmin - $rmax, $rstep step for $N mesh\n";
print "\n";

print     "# Potential type: $type\n";
print     "#  fCoulomb=$fCoulomb\n";

print "# Differentiation: h=$h  nOrder=$nDiffOrder\n";
print "\n";

my $Z1 = $Zi{$atom1};
my $Z2 = $Zi{$atom2};
my $De  = $Deij{$atom1}{$atom2};
my $a0  = $a0ij{$atom1}{$atom2};
my $r0  = $r0ij{$atom1}{$atom2};
my $CRn = $CRnij{$atom1}{$atom2};
my $A   = $Aij{$atom1}{$atom2};
my $rho = $rhoij{$atom1}{$atom2};
my $C   = $Cij{$atom1}{$atom2};
#my @args = ($fCoulomb, $Z1, $Z2, $De, $a0, $r0, $nRn, $CRn, $SplitEwaldSum, 0.0);
my @args = ($fCoulomb, $Z1, $Z2, $A, $rho, $C, $SplitEwaldSum, 0.0);
print "for $atom1-$atom2: ", join(', ', @args), "\n";

#===
# U, 微分(1階〜3階)の計算
#===
my ($pr, $pf, $pf1, $pf2, $pf3) = &CalDifferentials($nDiffOrder, $N, $rmin, $rstep, @args);

#===
# U, 微分(1階〜3階)のをdiff.csvに保存
#===
open(OUT, ">$DiffFile") or die "Can not write to [$DiffFile].\n";
print OUT "i,r,f,f',f'',f'''\n";
for(my $k = 0 ; $k < $N ; $k++) {
	printf OUT "%6d,%g,%g,%g,%g,%g\n", $k+1, $pr->[$k], $pf->[$k], $pf1->[$k], $pf2->[$k], $pf3->[$k];
}
close(OUT);

#===
# f'''=0となるrcを探す
#===
my ($rc, $fc, $f1c, $f2c, $f3c) = &Find3rdDerivativeZero($nDiffOrder, $pr, $pf, $pf1, $pf2, $pf2);

#===
# rc > rcmaxのときはrc=rcmaxとし、rcmaxにおけるU,微分を計算しなおす
#===
if(defined $rc and $rc > $rcmax) {
	print "Warning: Calculated rc=$rc exceeds the allowed rcmax=$rcmax\n";
	$rc = $rcmax;
	$fc  = Algorism::Interpolate($pr, $pf, $rc);
	$f1c = Algorism::Interpolate($pr, $pf1, $rc);
	$f2c = Algorism::Interpolate($pr, $pf2, $rc);
	$f3c = Algorism::Interpolate($pr, $pf3, $rc);
	print "         Use rc=$rcmax (f=$fc  f'=$f1c  f''=$f2c  f'''=$f3c\n";
	print "\n";
}

#===
# 解析解で計算。うまくいかないので使わない
#===
#USR(r) = B / r^n + Dr^2
#my ($Bp, $Dp, $np, $Bm, $Dm, $nm) = &CalBDnAnalytical($fc, $f1c, $f2c);

#===
# nを変えてrcにおけるf''の差を計算
# df''=0になるnを探すための補助だが、後のbinary searchの初期値を広くとればいいので
# 使わない
#===
if(0) {
for(my $n = 2.0 ; $n > 0.2 ; $n -= 0.2) {
	my ($B, $D, $df2) = &CalBD($rc, $fc, $f1c, $f2c, $n);

	my $USR2 = ($n+1)*$n*$B / $rc**($n+2) + 2.0*$D;
	printf "n=%3.2f  B=%12.6f  D=%12.6f: f''=%12.4g  df''=%12.6g\n", $n, $B, $D, $USR2, $df2;
#	my $USR  = $B / $rc**$n + $D * $rc*$rc;
#	my $USR1 = -$n*$B / $rc**($n+1) + 2.0*$D * $rc;
#	print "   USR(rc): f = $USR  f'=$USR1  f''=$USR2\n";
}
print "\n";
}

#===
# rcが見つかった場合のみ、Binary searchでdf''=0になるnを探す
#===
# Binary searchのnの初期範囲
my ($n0, $n1) = (0.01, 20.0);

my ($nm, $Bm, $Dm, $df2m);
if($rc) {
print "Binary search for df'' = 0\n";
my ($B0, $D0, $df20) = &CalBD($rc, $fc, $f1c, $f2c, $n0, 0);
my ($B1, $D1, $df21) = &CalBD($rc, $fc, $f1c, $f2c, $n1, 0);
print "Initial range:\n";
printf "  df''=%12.4g for n=$n0\n", $df20;
printf "  df''=%12.4g for n=$n1\n", $df21;
if($df20 * $df21 > 0.0) {
	print "Warning: Initial values must have different signs:\n";
	print "         rc is not found.\n";
	print "\n";
}

print "Binary search... (nMaxIter=$nMaxIter  EPS(f'')=$f2EPS)\n";
for(my $i = 0 ; $i < $nMaxIter ; $i++) {
	$nm = ($n0 + $n1) / 2.0;
	($Bm, $Dm, $df2m) = &CalBD($rc, $fc, $f1c, $f2c, $nm, 0);
	printf "  df''=%12.4g for n=(n0+n1)/2=$nm\n", $df2m;

	if(abs($df2m) < $f2EPS) {
		print "** Final solution: n=$nm  B=$Bm  D=$Dm: df''=$df2m\n";
		last;
	}
	elsif($df2m * $df20 < 0.0) {
		$n1 = $nm;
		next;
	}
	elsif($df2m * $df21 < 0.0) {
		$n0 = $nm;
		next;
	}
}
print "\n";
}

#===
# Uを計算。$SplitEwaldSum=1の場合、逆空間和 $KC*Zi*Zj/r*erf(ar)を引く
#===
my @U;
for(my $k = 0 ; $k < $N ; $k++) {
	my $r = $pr->[$k];
	if($r < $rc) {
		$U[$k] = $Bm / $r**$nm + $Dm * $r*$r;
	}
	else {
		$U[$k] = &Potential($r, @args);
	}
	$U[$k] -= $fCoulomb * &Coulomb($r, $Z1, $Z2, 0.0) * erf($alpha * $r) if($SplitEwaldSum);
}

#===
# Potential.tableへ保存
# 微分で7点公式を使った場合、最初と最後の3点ずつは計算できないので保存しない
#===
print "Save potential table to [$PotentialFile]\n";
open(OUT, ">$PotentialFile") or die "Can not write to [$PotentialFile].\n";

print OUT "#\n";
print OUT "# Potential type: $type\n";
print     "# Potential type: $type\n";
print OUT "#  fCoulomb=$fCoulomb\n";
print     "#  fCoulomb=$fCoulomb\n";
if($SplitEwaldSum) {
	print OUT "#  SplitEwaldSum=$SplitEwaldSum (Short range Coulomb term included in this table)\n";
	print     "#  SplitEwaldSum=$SplitEwaldSum (Short range Coulomb term included in this table)\n";
	print OUT "#  Ewald parameter (G vector)=$alpha\n";
	print     "#  Ewald parameter (G vector)=$alpha\n";
}
else {
	print OUT "#  SplitEwaldSum=$SplitEwaldSum (Full Coulomb term included in this table)\n";
	print     "#  SplitEwaldSum=$SplitEwaldSum (Full Coulomb term included in this table)\n";
}
print OUT "#\n";
print OUT "#  connected at rc = $rc (df''=$df2m)\n";
print OUT "#  r > rc: $atom1 - $atom2: A=$A  rho=$rho  C=$C\n";
#print OUT "#  r > rc: $atom1 - $atom2: De=$De  a0=$a0  r0=$r0\n";
print OUT "#  r < rc: n=$nm  B=$Bm  D=$Dm\n";
print OUT "# Differentiation: h=$h  nOrder=$nDiffOrder\n";
print     "# Differentiation: h=$h  nOrder=$nDiffOrder\n";
print OUT "#\n";
print OUT "\n";
print "for $atom1-$atom2: ", join(', ', @args), "\n";
print OUT "$atom1-$atom2\n";
printf OUT "N\t\t%d\tR\t%g\t%g\n", $N - 6, $rmin + $rstep*3, $rmax - $rstep*3;
printf     "N\t\t%d\tR\t%g\t%g\n", $N - 6, $rmin + $rstep*3, $rmax - $rstep*3;
print OUT "\n";
my $c = 1;
for(my $k = 3 ; $k < $N-3 ; $k++) {
	my $r = $pr->[$k];
	my $F;
	if($k < 3 or $k >= $N-3) {
		$F = 0.0;
	}
	else {
		$F = -Algorism::Differentiate7WithIndex(\@U, $k, $rstep);
	}

	printf OUT "%6d\t\t%8.6f\t\t%12.6g\t\t%12.6g\n", $c, $r, $U[$k], $F;
#	printf     "%6d\t\t%8.6f\t\t%12.6g\t\t%12.6g\n", $c, $r, $U[$k], $F;
	$c++;
}
print OUT "\n";
close(OUT);

exit;


#===============================
# Functions
#===============================
sub Potential
{
#	my ($r, $fCoulomb, $Z1, $Z2, $De, $a0,  $r0, $nRn, $CRn, $SplitEwaldSum, $alpha) = @_;
	my ($r, $fCoulomb, $Z1, $Z2, $A,  $rho, $C,  $SplitEwaldSum, $alpha) = @_;

	$alpha = 0.0 if(!$SplitEwaldSum);
	return $fCoulomb * &Coulomb($r, $Z1, $Z2, $alpha) 
#		 + &Morse($r, $De, $a0, $r0)
#		 + &Rn($r, $nRn, $CRn);
		 + &Born($r, $A, $rho)
		 + &Dispersion($r, $C);
}

sub DifferentiatePotential
{
	my ($nDiffOrder, $r, @args) = @_;

		if($nDiffOrder == 2) {
			return Algorism::Differentiate2(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}
		elsif($nDiffOrder == 3) {
			return Algorism::Differentiate3(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}
		elsif($nDiffOrder == 5) {
			return Algorism::Differentiate5(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}
		elsif($nDiffOrder == 7) {
			return Algorism::Differentiate7(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}

		print "Error: Invalid nOrder=$nDiffOrder for differentiation "
			 ."(Choose from 2, 3, 5, or 7)\n";
		exit;
}

sub Coulomb
{
	my ($r, $Z1, $Z2, $alpha) = @_;
	my $U = $KC * $Z1 * $Z2 / $r * erfc($alpha * $r);
#print "UC=$U (KC=$KC)\n";
	return $U;
}

sub Born
{
	my ($r, $A, $rho) = @_;
	return $A * exp(-$r / $rho);
}

sub Morse
{
	my ($r, $De, $a0, $r0) = @_;

	my $v = 1.0 - exp(-$a0 * ($r - $r0));
	my $U = $De * ($v*$v - 1.0);
#print "UM=$U\n";
	return $U;
}

sub Dispersion
{
	my ($r, $C) = @_;

	return -&Rn($r, 6, $C);
}

sub Rn
{
	my ($r, $nRn, $CR) = @_;

	return $CR / $r**$nRn;
}

sub CalDifferentials
{
	my ($nDiffOrder, $N, $rmin, $rstep, @args) = @_;

	my (@r, @f, @f1, @f2, @f3);
	for(my $k = 0 ; $k < $N ; $k++) {
		$r[$k]  = $rmin + $k * $rstep;
		$f[$k]  = &Potential($r[$k], @args);
		$f1[$k] = &DifferentiatePotential($nDiffOrder, $r[$k], @args);
	}

	for(my $k = 3 ; $k < $N-3 ; $k++) {
		$f2[$k] = Algorism::Differentiate7WithIndex(\@f1, $k, $rstep);
	}

	for(my $k = 6 ; $k < $N-6 ; $k++) {
		$f3[$k] = Algorism::Differentiate7WithIndex(\@f2, $k, $rstep);
	}

	return (\@r, \@f, \@f1, \@f2, \@f3);
}

sub Find3rdDerivativeZero
{
	my ($nDiffOrder, $pr, $pf, $pf1, $pf2, $pf2) = @_;

	my $rc;
	print "Find zero 3rd derivative radius rc\n";
	for(my $k = 7 ; $k < $N-6 ; $k++) {
		my $f1 = $pf3->[$k-1];
		my $f2 = $pf3->[$k];
		if($f1 * $f2 <= 0.0) {
			my $r1 = $rmin + ($k-1) * $rstep;
			my $r2 = $rmin + $k * $rstep;
			$rc = $r1 + (0.0 - $f1) / ($f2 - $f1) * $rstep;

			$fc  = Algorism::InterpolateBySimpson($pr, $pf,  $rc);
			$f1c = Algorism::InterpolateBySimpson($pr, $pf1, $rc);
			$f2c = Algorism::InterpolateBySimpson($pr, $pf2, $rc);
			$f3c = Algorism::InterpolateBySimpson($pr, $pf3, $rc);

			print "\n";
			print "  3rd derivative becomes 0 at rc = $rc\n";
			print "     cf: index = ", $k-1, " - $k: r = $r1 - $r2: f''' = $f1 - $f2\n";
			print "    f(rc)    = $fc\n";
			print "    f'(rc)   = $f1c\n";
			print "    f''(rc)  = $f2c\n";
			print "    f'''(rc) = $f3c\n";

			last;
		}
	}
	if(!defined $rc) {
		print "Error in Find3rdDerivativeZero: rc could not be found\n";
	}
	print "\n";

	return ($rc, $fc, $f1c, $f2c, $f3c);
}

sub CalBD
{
	my ($rc, $fc, $f1c, $f2c, $n, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

	print "Calculate B and D from given n = $n\n" if($IsPrint);

#  1. B * (1 / rc^n)      + D * (rc^2)   = f(rc)
#  2. B * (-n / rc^(n+1)) + D * (2 * rc) = f1(rc)
	my $a11 = 1.0 / $rc**$n;
	my $a12 = $rc * $rc;
	my $a21 = -$n / $rc**($n+1);
	my $a22 = 2.0 * $rc;

	my $det = $a11 * $a22 - $a12 * $a21;
	my $B  = ($fc  * $a22 - $f1c * $a12) / $det;
	my $D  = ($f1c * $a11 - $fc  * $a21) / $det;

##  3.  n(n+1)B / rc^(n+2) + 2*D = f2(rc)
	my $df2 = $n*($n+1) * $B / $rc**($n+2) + 2.0 * $D - $f2c;

	return ($B, $D, $df2);
}

sub CalBDnAnalytical
{
	my ($fc, $f1c, $f2c) = @_;

	print "Calculate B, D, and n from analytical solution\n";
	
# USR(r) = B / r^n + Dr^2
# Conditions:
#  1. B / rc^n + D * rc^2 = f(rc)
#     => B = (f - D * rc^2) * rc^n
#  2. -nB / rc^(n+1) + 2*D * rc = f1(rc)
#     => D = (f1 + n(f - D * rc^2) / rc) / (2*rc) 
#  3.  n(n+1)B / rc^(n+2) + 2*D = f2(rc)
#
# SQR = sqrt(20 f1^2 - 24 f1*f2 + 4f2^2 - 28f1 * f + 20 f2*f + 9 f^2)
# D = (1/8) * (6f1 - 2f2 - 3f +- SQR)
# n = (-2f1 - 2f2 + 3f -+ SQR) / [4(f1 - f)]
# B = (f - Drc^2) * rc^n

	my $SQR = sqrt(20.0 * $f1c * $f1c - 24.0 * $f1c * $f2c + 4.0 * $f2c * $f2c 
				 - 28.0 * $f1c * $fc  + 20.0 * $f2c * $fc  + 9.0 * $fc  * $fc);
	my $Dp = (1.0/8.0) * (6.0 * $f1c - 2.0 * $f2c - 3.0 * $fc + $SQR);
	my $np = (-2.0 * $f1c - 2.0 * $f2c + 3.0 * $fc - $SQR) / 4.0 / ($f1c - $fc);
	my $Bp = ($fc - $Dp * $rc * $rc) * $rc**$np;
	my $Dm = (1.0/8.0) * (6.0 * $f1c - 2.0 * $f2c - 3.0 * $fc - $SQR);
	my $nm = (-2.0 * $f1c - 2.0 * $f2c + 3.0 * $fc + $SQR) / 4.0 / ($f1c - $fc);
	my $Bm = ($fc - $Dm * $rc * $rc) * $rc**$nm;

	my $USRp  = $Bp / $rc**$np + $Dp * $rc*$rc;
	my $USR1p = -$np*$Bp / $rc**($np+1) + 2.0*$Dp * $rc;
	my $USR2p = ($np+1)*$np*$Bp / $rc**($np+2) + 2.0*$Dp;
	my $USRm  = $Bm / $rc**$nm + $Dm * $rc*$rc;
	my $USR1m = -$nm*$Bm / $rc**($nm+1) + 2.0*$Dm * $rc;
	my $USR2m = ($nm+1)*$nm*$Bm / $rc**($nm+2) + 2.0*$Dm;

	print "+: B = $Bp  D = $Dp  n = $np\n";
	print "   USR(rc): f = $USRp  f'=$USR1p  f''=$USR2p\n";
	print "-: B = $Bm  D = $Dm  n = $nm\n";
	print "   USR(rc): f = $USRm  f'=$USR1m  f''=$USR2m\n";
	print "\n";

	return ($Bp, $Dp, $np, $Bm, $Dm, $nm);
}