508 lines
16 KiB
C++
Executable File
508 lines
16 KiB
C++
Executable File
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: TBA_XXZ.cc
|
|
|
|
Purpose: TBA for the gapless XXZ antiferromagnet.
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
// First, define some useful kernels:
|
|
|
|
DP XXZ_phi1_kernel (DP zeta, DP lambda)
|
|
{
|
|
return(2.0 * atan(tanh(lambda)/tan(0.5*zeta)));
|
|
}
|
|
|
|
DP XXZ_phi2_kernel (DP zeta, DP lambda)
|
|
{
|
|
return(2.0 * atan(tanh(lambda)/tan(zeta)));
|
|
}
|
|
|
|
DP XXZ_a1_kernel (DP sinzeta, DP coszeta, DP lambda)
|
|
{
|
|
DP expmin2abslambda = exp(-2.0 * fabs(lambda));
|
|
|
|
DP answer = sinzeta * expmin2abslambda
|
|
/(PI * (0.5 * (1.0 + expmin2abslambda * expmin2abslambda) - coszeta * expmin2abslambda));
|
|
|
|
return(answer);
|
|
}
|
|
|
|
DP XXZ_da1dlambda_kernel (DP sinzeta, DP coszeta, DP lambda)
|
|
{
|
|
DP expmin2abslambda = exp(-2.0 * fabs(lambda));
|
|
DP signlambda = lambda >= 0.0 ? 1.0 : -1.0;
|
|
DP answer = -sinzeta * signlambda * expmin2abslambda * (1.0 - expmin2abslambda * expmin2abslambda)
|
|
/(PI * pow((0.5 * (1.0 + expmin2abslambda * expmin2abslambda)
|
|
- coszeta * expmin2abslambda), 2.0));
|
|
|
|
return(answer);
|
|
}
|
|
|
|
DP XXZ_a2_kernel (DP sin2zeta, DP cos2zeta, DP lambda)
|
|
{
|
|
DP expmin2abslambda = exp(-2.0 * fabs(lambda));
|
|
|
|
DP answer = sin2zeta * expmin2abslambda
|
|
/(PI * (0.5 * (1.0 + expmin2abslambda * expmin2abslambda) - cos2zeta * expmin2abslambda));
|
|
|
|
return(answer);
|
|
}
|
|
|
|
void Iterate_XXZ_rhotot_GS (DP zeta, DP B, Root_Density& rhotot_GS)
|
|
{
|
|
// Performs an iteration of the TBA equation for rhotot_GS:
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < rhotot_GS.Npts; ++i) rhotot_GS.prev_value[i] = rhotot_GS.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sinzeta = sin(zeta);
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP coszeta = cos(zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < rhotot_GS.Npts; ++i) {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < rhotot_GS.Npts; ++j)
|
|
conv += fabs(rhotot_GS.lambda[j]) > B ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, rhotot_GS.lambda[i] - rhotot_GS.lambda[j])
|
|
* rhotot_GS.prev_value[j] * rhotot_GS.dlambda[j];
|
|
rhotot_GS.value[i] = XXZ_a1_kernel(sinzeta, coszeta, rhotot_GS.lambda[i]) - conv;
|
|
}
|
|
|
|
// Calculate the sum of differences:
|
|
rhotot_GS.diff = 0.0;
|
|
for (int i = 0; i < rhotot_GS.Npts; ++i)
|
|
rhotot_GS.diff += fabs(rhotot_GS.value[i] - rhotot_GS.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_rhotot_GS (DP Delta, DP B, DP lambdamax, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the ground state
|
|
// root distribution of the XXZ model.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_rhotot_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density rhotot_GS(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
int niter = 0;
|
|
do {
|
|
Iterate_XXZ_rhotot_GS (zeta, B, rhotot_GS);
|
|
niter ++;
|
|
} while (rhotot_GS.diff > req_prec);
|
|
|
|
return(rhotot_GS);
|
|
}
|
|
|
|
DP Return_GS_Sz_tot_value (DP B, Root_Density& rhotot_GS)
|
|
{
|
|
DP integral = 0.0;
|
|
for (int i = 0; i < rhotot_GS.Npts; ++i)
|
|
integral += fabs(rhotot_GS.lambda[i]) > B ? 0.0 : rhotot_GS.value[i] * rhotot_GS.dlambda[i];
|
|
|
|
return(0.5 - integral);
|
|
}
|
|
|
|
void Iterate_XXZ_eps_GS (DP zeta, DP Hz, Root_Density& eps_GS)
|
|
{
|
|
// Performs an iteration of the TBA equation for eps_GS:
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < eps_GS.Npts; ++i) eps_GS.prev_value[i] = eps_GS.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sinzeta = sin(zeta);
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP coszeta = cos(zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < eps_GS.Npts; ++i) {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < eps_GS.Npts; ++j)
|
|
conv += eps_GS.prev_value[j] > 0.0 ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, eps_GS.lambda[i] - eps_GS.lambda[j])
|
|
* eps_GS.prev_value[j] * eps_GS.dlambda[j];
|
|
eps_GS.value[i] = Hz - PI * sinzeta * XXZ_a1_kernel(sinzeta, coszeta, eps_GS.lambda[i]) - conv;
|
|
}
|
|
|
|
// Calculate the sum of differences:
|
|
eps_GS.diff = 0.0;
|
|
for (int i = 0; i < eps_GS.Npts; ++i)
|
|
eps_GS.diff += fabs(eps_GS.value[i] - eps_GS.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_eps_GS (DP Delta, DP Hz, DP lambdamax, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the ground state
|
|
// epsilon function of the XXZ model.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_eps_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density eps_GS(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
int niter = 0;
|
|
do {
|
|
Iterate_XXZ_eps_GS (zeta, Hz, eps_GS);
|
|
niter++;
|
|
} while (eps_GS.diff > req_prec);
|
|
|
|
return(eps_GS);
|
|
}
|
|
|
|
void Iterate_XXZ_depsdlambda_GS (DP zeta, DP B, Root_Density& depsdlambda_GS)
|
|
{
|
|
// Performs an iteration of the TBA equation for depsdlambda_GS:
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < depsdlambda_GS.Npts; ++i) depsdlambda_GS.prev_value[i] = depsdlambda_GS.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sinzeta = sin(zeta);
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP coszeta = cos(zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < depsdlambda_GS.Npts; ++i) {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < depsdlambda_GS.Npts; ++j)
|
|
conv += fabs(depsdlambda_GS.lambda[j]) > B ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, depsdlambda_GS.lambda[i] - depsdlambda_GS.lambda[j])
|
|
* depsdlambda_GS.prev_value[j] * depsdlambda_GS.dlambda[j];
|
|
depsdlambda_GS.value[i] = -PI * sinzeta * XXZ_da1dlambda_kernel(sinzeta, coszeta, depsdlambda_GS.lambda[i]) - conv;
|
|
}
|
|
|
|
// Calculate the sum of differences:
|
|
depsdlambda_GS.diff = 0.0;
|
|
for (int i = 0; i < depsdlambda_GS.Npts; ++i)
|
|
depsdlambda_GS.diff += fabs(depsdlambda_GS.value[i] - depsdlambda_GS.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_depsdlambda_GS (DP Delta, DP B, DP lambdamax, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the ground state
|
|
// (d/d\lambda) epsilon function of the XXZ model.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_rho_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density depsdlambda_GS(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
int niter = 0;
|
|
do {
|
|
Iterate_XXZ_depsdlambda_GS (zeta, B, depsdlambda_GS);
|
|
niter++;
|
|
} while (depsdlambda_GS.diff > req_prec);
|
|
|
|
return(depsdlambda_GS);
|
|
}
|
|
|
|
void Iterate_XXZ_b2BB_lambda_B (DP zeta, DP B, Root_Density& b2BB)
|
|
{
|
|
// Calculates the vector corresponding to b2BB (lambda,B)
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < b2BB.Npts; ++i) b2BB.prev_value[i] = b2BB.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < b2BB.Npts; ++i) {
|
|
if (fabs(b2BB.lambda[i]) > B) b2BB.value[i] = 0.0;
|
|
else {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < b2BB.Npts; ++j)
|
|
conv += fabs(b2BB.lambda[j]) > B ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, b2BB.lambda[i] - b2BB.lambda[j])
|
|
* b2BB.prev_value[j] * b2BB.dlambda[j];
|
|
b2BB.value[i] = -XXZ_a2_kernel(sin2zeta, cos2zeta, b2BB.lambda[i] - B)
|
|
- conv;
|
|
}
|
|
}
|
|
// Calculate the sum of differences:
|
|
b2BB.diff = 0.0;
|
|
for (int i = 0; i < b2BB.Npts; ++i)
|
|
b2BB.diff += fabs(b2BB.value[i] - b2BB.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_b2BB_lambda_B (DP Delta, DP B, DP lambdamax, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the inverse
|
|
// kernel b2BB (lambda, B) used in Kbar function.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_Kbackflow_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density b2BB_lambda_B(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
int niter = 0;
|
|
do {
|
|
Iterate_XXZ_b2BB_lambda_B (zeta, B, b2BB_lambda_B);
|
|
niter++;
|
|
} while (b2BB_lambda_B.diff > req_prec);
|
|
|
|
return(b2BB_lambda_B);
|
|
}
|
|
|
|
void Iterate_XXZ_b2BB_lambda_lambdap (DP zeta, DP B, DP lambdap, Root_Density& b2BB)
|
|
{
|
|
// Calculates the vector corresponding to b2BB (lambda,lambdap)
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < b2BB.Npts; ++i) b2BB.prev_value[i] = b2BB.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < b2BB.Npts; ++i) {
|
|
if (fabs(b2BB.lambda[i]) > B) b2BB.value[i] = 0.0;
|
|
else {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < b2BB.Npts; ++j)
|
|
conv += fabs(b2BB.lambda[j]) > B ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, b2BB.lambda[i] - b2BB.lambda[j])
|
|
* b2BB.prev_value[j] * b2BB.dlambda[j];
|
|
b2BB.value[i] = -XXZ_a2_kernel(sin2zeta, cos2zeta, b2BB.lambda[i] - lambdap)
|
|
- conv;
|
|
}
|
|
}
|
|
// Calculate the sum of differences:
|
|
b2BB.diff = 0.0;
|
|
for (int i = 0; i < b2BB.Npts; ++i)
|
|
b2BB.diff += fabs(b2BB.value[i] - b2BB.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_b2BB_lambda_lambdap (DP Delta, DP B, DP lambdap, DP lambdamax, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the inverse
|
|
// kernel b2BB (lambda, lambdap) used in Kbar function.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_Kbackflow_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density b2BB_lambda_lambdap(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
int niter = 0;
|
|
do {
|
|
Iterate_XXZ_b2BB_lambda_lambdap (zeta, B, lambdap, b2BB_lambda_lambdap);
|
|
niter++;
|
|
} while (b2BB_lambda_lambdap.diff > req_prec);
|
|
|
|
return(b2BB_lambda_lambdap);
|
|
}
|
|
|
|
void Iterate_XXZ_Kbackflow_GS (DP zeta, DP B, DP lambda_p, DP lambda_h, Root_Density& Kbackflow_GS)
|
|
{
|
|
// Performs an iteration of the TBA equation for Kbackflow_GS:
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < Kbackflow_GS.Npts; ++i) Kbackflow_GS.prev_value[i] = Kbackflow_GS.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < Kbackflow_GS.Npts; ++i) {
|
|
if (false && fabs(Kbackflow_GS.lambda[i]) > B) Kbackflow_GS.value[i] = 0.0;
|
|
else {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < Kbackflow_GS.Npts; ++j)
|
|
conv += fabs(Kbackflow_GS.lambda[j]) > B ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, Kbackflow_GS.lambda[i] - Kbackflow_GS.lambda[j])
|
|
* Kbackflow_GS.prev_value[j] * Kbackflow_GS.dlambda[j];
|
|
Kbackflow_GS.value[i] = -XXZ_a2_kernel(sin2zeta, cos2zeta, Kbackflow_GS.lambda[i] - lambda_p)
|
|
+ XXZ_a2_kernel(sin2zeta, cos2zeta, Kbackflow_GS.lambda[i] - lambda_h) - conv;
|
|
}
|
|
}
|
|
// Calculate the sum of differences:
|
|
Kbackflow_GS.diff = 0.0;
|
|
for (int i = 0; i < Kbackflow_GS.Npts; ++i)
|
|
Kbackflow_GS.diff += fabs(Kbackflow_GS.value[i] - Kbackflow_GS.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_Kbackflow_GS (DP Delta, DP B, DP lambdamax, DP lambda_p, DP lambda_h, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the ground state
|
|
// backflow function of the XXZ model.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_Kbackflow_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density Kbackflow_GS(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
int niter = 0;
|
|
do {
|
|
Iterate_XXZ_Kbackflow_GS (zeta, B, lambda_p, lambda_h, Kbackflow_GS);
|
|
niter++;
|
|
} while (Kbackflow_GS.diff > req_prec);
|
|
|
|
return(Kbackflow_GS);
|
|
}
|
|
|
|
void Iterate_XXZ_Fbackflow_GS (DP zeta, DP B, DP lambda_p, DP lambda_h, Root_Density& Fbackflow_GS)
|
|
{
|
|
// Performs an iteration of the TBA equation for Fbackflow_GS:
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < Fbackflow_GS.Npts; ++i) Fbackflow_GS.prev_value[i] = Fbackflow_GS.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < Fbackflow_GS.Npts; ++i) {
|
|
if (fabs(Fbackflow_GS.lambda[i]) > B) Fbackflow_GS.value[i] = 0.0;
|
|
else {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < Fbackflow_GS.Npts; ++j)
|
|
conv += fabs(Fbackflow_GS.lambda[j]) > B ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, Fbackflow_GS.lambda[i] - Fbackflow_GS.lambda[j])
|
|
* Fbackflow_GS.prev_value[j] * Fbackflow_GS.dlambda[j];
|
|
Fbackflow_GS.value[i] = (-XXZ_phi2_kernel(zeta, Fbackflow_GS.lambda[i] - lambda_p)
|
|
+ XXZ_phi2_kernel(zeta, Fbackflow_GS.lambda[i] - lambda_h))/twoPI - conv;
|
|
}
|
|
}
|
|
// Calculate the sum of differences:
|
|
Fbackflow_GS.diff = 0.0;
|
|
for (int i = 0; i < Fbackflow_GS.Npts; ++i)
|
|
Fbackflow_GS.diff += fabs(Fbackflow_GS.value[i] - Fbackflow_GS.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_Fbackflow_GS (DP Delta, DP B, DP lambdamax, DP lambda_p, DP lambda_h, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the ground state
|
|
// F backflow function of the XXZ model.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_Fbackflow_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density Fbackflow_GS(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
int niter = 0;
|
|
do {
|
|
Iterate_XXZ_Fbackflow_GS (zeta, B, lambda_p, lambda_h, Fbackflow_GS);
|
|
niter++;
|
|
} while (Fbackflow_GS.diff > req_prec);
|
|
|
|
return(Fbackflow_GS);
|
|
}
|
|
|
|
void Iterate_XXZ_Z_GS (DP zeta, DP B, Root_Density& Z_GS)
|
|
{
|
|
// Performs an iteration of the TBA equation for the dressed charge:
|
|
|
|
// First, copy the actual values into the previous ones:
|
|
for (int i = 0; i < Z_GS.Npts; ++i) Z_GS.prev_value[i] = Z_GS.value[i];
|
|
|
|
// Calculate the new values:
|
|
DP conv = 0.0;
|
|
DP sin2zeta = sin(2.0*zeta);
|
|
DP cos2zeta = cos(2.0*zeta);
|
|
for (int i = 0; i < Z_GS.Npts; ++i) {
|
|
// First, calculate the convolution
|
|
conv = 0.0;
|
|
for (int j = 0; j < Z_GS.Npts; ++j)
|
|
conv += fabs(Z_GS.lambda[j]) > B ? 0.0 :
|
|
XXZ_a2_kernel (sin2zeta, cos2zeta, Z_GS.lambda[i] - Z_GS.lambda[j])
|
|
* Z_GS.prev_value[j] * Z_GS.dlambda[j];
|
|
Z_GS.value[i] = 1.0 - conv;
|
|
}
|
|
|
|
// Calculate the sum of differences:
|
|
Z_GS.diff = 0.0;
|
|
for (int i = 0; i < Z_GS.Npts; ++i)
|
|
Z_GS.diff += fabs(Z_GS.value[i] - Z_GS.prev_value[i]);
|
|
|
|
return;
|
|
}
|
|
|
|
Root_Density XXZ_Z_GS (DP Delta, DP B, DP lambdamax, int Npts, DP req_prec)
|
|
{
|
|
|
|
// This function returns a Root_Density object corresponding to the ground state
|
|
// dressed charge function of the XXZ model.
|
|
|
|
if (Delta >= 1.0 || Delta <= 0.0) ABACUSerror("Delta out of bounds in XXZ_Z_GS.");
|
|
|
|
DP zeta = acos(Delta);
|
|
|
|
Root_Density Z_GS(Npts, lambdamax);
|
|
|
|
// We now attempt to find a solution...
|
|
|
|
do {
|
|
Iterate_XXZ_Z_GS (zeta, B, Z_GS);
|
|
} while (Z_GS.diff > req_prec);
|
|
|
|
return(Z_GS);
|
|
}
|
|
|
|
|
|
} // namespace ABACUS
|