diff options
| -rw-r--r-- | src/misc.c | 18 | ||||
| -rw-r--r-- | src/vector.c | 4 |
2 files changed, 15 insertions, 7 deletions
@@ -352,6 +352,7 @@ double interp2d(double x, double y, double *xp, double *yp, double *zp, size_t n * See https://en.wikipedia.org/wiki/Bilinear_interpolation. */ size_t i, j; double q11, q12, q21, q22; + double idx, idy; if (x <= xp[0]) { fprintf(stderr, "x < xp[0]!\n"); @@ -363,14 +364,18 @@ double interp2d(double x, double y, double *xp, double *yp, double *zp, size_t n exit(1); } - i = (x-xp[0])/(xp[1]-xp[0]); + idx = 1.0/(xp[1]-xp[0]); + + i = (x-xp[0])*idx; if (i > n1-2) { fprintf(stderr, "i > n1 - 2!\n"); exit(1); } - j = (y-yp[0])/(yp[1]-yp[0]); + idy = 1.0/(yp[1]-yp[0]); + + j = (y-yp[0])*idy; if (j > n2-2) { fprintf(stderr, "j > n2 - 2!\n"); @@ -382,7 +387,7 @@ double interp2d(double x, double y, double *xp, double *yp, double *zp, size_t n q21 = zp[(i+1)*n2+j]; q22 = zp[(i+1)*n2+j+1]; - return (q11*(xp[i+1]-x)*(yp[j+1]-y) + q21*(x-xp[i])*(yp[j+1]-y) + q12*(xp[i+1]-x)*(y-yp[j]) + q22*(x-xp[i])*(y-yp[j]))/((xp[i+1]-xp[i])*(yp[i+1]-yp[i])); + return (q11*(xp[i+1]-x)*(yp[j+1]-y) + q21*(x-xp[i])*(yp[j+1]-y) + q12*(xp[i+1]-x)*(y-yp[j]) + q22*(x-xp[i])*(y-yp[j]))*idx*idy; } double interp1d(double x, double *xp, double *yp, size_t n) @@ -392,14 +397,17 @@ double interp1d(double x, double *xp, double *yp, size_t n) * * If x < xp[0] returns yp[0] and if x > xp[n-1] returns yp[n-1]. */ size_t i; + double idx; if (x <= xp[0]) return yp[0]; - i = (x-xp[0])/(xp[1]-xp[0]); + idx = 1.0/(xp[1]-xp[0]); + + i = (x-xp[0])*idx; if (i > n-2) return yp[n-1]; - return yp[i] + (yp[i+1]-yp[i])*(x-xp[i])/(xp[i+1]-xp[i]); + return yp[i] + (yp[i+1]-yp[i])*(x-xp[i])*idx; } int isclose(double a, double b, double rel_tol, double abs_tol) diff --git a/src/vector.c b/src/vector.c index d0e304a..061c3d6 100644 --- a/src/vector.c +++ b/src/vector.c @@ -4,8 +4,8 @@ void normalize(double *a) { double c; - c = NORM(a); - DIV(a,c); + c = 1.0/NORM(a); + MUL(a,c); } void rotate(double *dest, double *v, double *n, double phi) |