adjust for double precision floating point

This commit is contained in:
Axel Kohlmeyer
2022-01-04 23:02:17 -05:00
parent d0f203127d
commit 40abc0886c
2 changed files with 4 additions and 4 deletions

View File

@ -578,7 +578,7 @@ void AngleTable::spline(double *x, double *y, int n,
double p,qn,sig,un; double p,qn,sig,un;
double *u = new double[n]; double *u = new double[n];
if (yp1 > 0.99e30) y2[0] = u[0] = 0.0; if (yp1 > 0.99e300) y2[0] = u[0] = 0.0;
else { else {
y2[0] = -0.5; y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1); u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
@ -590,7 +590,7 @@ void AngleTable::spline(double *x, double *y, int n,
u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]); u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p; u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
} }
if (ypn > 0.99e30) qn = un = 0.0; if (ypn > 0.99e300) qn = un = 0.0;
else { else {
qn = 0.5; qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2])); un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));

View File

@ -541,7 +541,7 @@ void BondTable::spline(double *x, double *y, int n,
double p,qn,sig,un; double p,qn,sig,un;
double *u = new double[n]; double *u = new double[n];
if (yp1 > 0.99e30) y2[0] = u[0] = 0.0; if (yp1 > 0.99e300) y2[0] = u[0] = 0.0;
else { else {
y2[0] = -0.5; y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1); u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
@ -553,7 +553,7 @@ void BondTable::spline(double *x, double *y, int n,
u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]); u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p; u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
} }
if (ypn > 0.99e30) qn = un = 0.0; if (ypn > 0.99e300) qn = un = 0.0;
else { else {
qn = 0.5; qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2])); un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));