This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Data Structures | |
| struct | polynom3 |
Defines | |
| #define | POLYMAX 100 |
| #define | SP_NATURAL 1 |
| #define | SP_DERIVATIVE 2 |
| #define | SP_PERIODIC 3 |
| #define | splint(poly, x) (poly.a + (x)*(poly.b + (x)*(poly.c + (x)*poly.d))) |
| #define | splint1(poly, x) (poly.b + (x)*(2.*poly.c + (x)*3.*poly.d)) |
| #define | splint2(poly, x) (2.*poly.c + (x)*6.*poly.d) |
Functions | |
| real | polynomial (real *xp, real *yp, real *z, real x, real y) |
| void | CubicPolynom (real x0, real y0, real x1, real y1, real x2, real y2, real x3, real y3, real *ex, real *fx, real *gx, real *ey, real *fy, real *gy, real *rd2) |
| void | spline (real1D x, real1D y, int n, real yp1, real ypn, real1D d2y, polynom3 *poly, int flag) |
| void | tridag (real1D a, real1D b, real1D c, real1D r, real1D u, int n) |
| void | Pspline (real1D x, real1D y, real1D d2y, real lp, int n, polynom3 *poly) |
| real | splint_find (real1D xa, polynom3 *poly, int n, real x) |
| int | distmin (polynom3 px, polynom3 py, real xp, real yp, real *tmin, real t1, real t2, real prec) |
| int | polyroot (polynom3 poly, real x1, real x2, real *x, real prec) |
| int | polyroots (polynom3 poly, real x1, real x2, real *x, real y, real prec) |
| void | bcucof (real2D f, int i, int j, real *c, int nx, int ny) |
| real | bcuint (real t, real u, real *c) |
| real | polydistmin (polynom3 px1, polynom3 px2, polynom3 py1, polynom3 py2, real *t1, real *t2, real ftol) |
| real | trapzd (real(*func)(polynom3, polynom3, real), polynom3 px, polynom3 py, real a, real b, int n) |
| real | qtrap (real(*func)(polynom3, polynom3, real), polynom3 px, polynom3 py, real a, real b, real eps) |
|
|
Definition at line 2 of file interpolation.h. Referenced by distmin(), and polyroot(). |
|
|
Definition at line 4 of file interpolation.h. Referenced by interface_splines(), and spline(). |
|
|
Definition at line 3 of file interpolation.h. Referenced by interface_splines(), and spline(). |
|
|
Definition at line 5 of file interpolation.h. |
|
|
|
Definition at line 8 of file interpolation.h. Referenced by advect(), avg_kappa_normals(), closest_normal(), curvature_write(), extra_poly(), interfaces_xmoment(), interfaces_ymoment(), kappa_normals(), and polyroot(). |
|
|
Definition at line 9 of file interpolation.h. Referenced by avg_kappa_normals(), curvature_write(), and kappa_normals(). |
|
||||||||||||||||||||||||||||
|
Definition at line 376 of file interpolation.c.
00377 {
00378 static int wt[16][16]= {
00379 {1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00380 {0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0},
00381 {-3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0},
00382 {2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0},
00383 {0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0},
00384 {0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0},
00385 {0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1},
00386 {0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1},
00387 {-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0},
00388 {0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0},
00389 {9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2},
00390 {-6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2},
00391 {2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0},
00392 {0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0},
00393 {-6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1},
00394 {4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1}
00395 };
00396 int k;
00397 real xx, x[16];
00398
00399 /* function */
00400 x[0] = f[i][j]; x[1] = f[i+1][j]; x[2] = f[i+1][j+1]; x[3] = f[i][j+1];
00401
00402 /* x component of the gradient */
00403 x[4] = 0.5*(f[i+1][j] - f[i-1][j]);
00404 if (i+2 <= nx) {
00405 x[5] = 0.5*(f[i+2][j] - f[i][j]);
00406 x[6] = 0.5*(f[i+2][j+1] - f[i][j+1]);
00407 }
00408 else { /* periodic boundary conditions */
00409 x[5] = 0.5*(f[i+2 - nx + 2][j] - f[i][j]);
00410 x[6] = 0.5*(f[i+2 - nx + 2][j+1] - f[i][j+1]);
00411 }
00412 x[7] = 0.5*(f[i+1][j+1] - f[i-1][j+1]);
00413
00414 /* y component of the gradient */
00415 x[8] = 0.5*(f[i][j+1] - f[i][j-1]);
00416 x[9] = 0.5*(f[i+1][j+1] - f[i+1][j-1]);
00417 if (j+2 <= ny) {
00418 x[10] = 0.5*(f[i+1][j+2] - f[i+1][j]);
00419 x[11] = 0.5*(f[i][j+2] - f[i][j]);
00420 }
00421 else /* zero gradient */
00422 x[10] = x[11] = 0.0;
00423
00424 /* cross derivatives */
00425 x[12] = 0.25*(f[i+1][j+1] - f[i+1][j-1] - f[i-1][j+1] + f[i-1][j-1]);
00426 if (i+2 <= nx) {
00427 x[13] = 0.25*(f[i+2][j+1] - f[i+2][j-1] - f[i][j+1] + f[i][j-1]);
00428 if (j+2 <= ny) {
00429 x[14] = 0.25*(f[i+2][j+2] - f[i+2][j] - f[i][j+2] + f[i][j]);
00430 x[15] = 0.25*(f[i+1][j+2] - f[i+1][j] - f[i-1][j+2] + f[i-1][j]);
00431 }
00432 else /* zero gradient */
00433 x[14] = x[15] = 0.0;
00434 }
00435 else { /* periodic boundary conditions */
00436 x[13] = 0.25*(f[i+2 - nx + 2][j+1] - f[i+2 - nx + 2][j-1]
00437 - f[i][j+1] + f[i][j-1]);
00438 if (j+2 <= ny) {
00439 x[14] = 0.25*(f[i+2-nx+2][j+2] - f[i+2-nx+2][j] - f[i][j+2] + f[i][j]);
00440 x[15] = 0.25*(f[i+1][j+2] - f[i+1][j] - f[i-1][j+2] + f[i-1][j]);
00441 }
00442 else /* zero gradient */
00443 x[14] = x[15] = 0.0;
00444 }
00445
00446 for (i = 0; i <= 15; i++) {
00447 xx = 0.0;
00448 for (k = 0; k <= 15; k++) xx += wt[i][k]*x[k];
00449 c[i] = xx;
00450 }
00451 }
|
|
||||||||||||||||
|
Definition at line 454 of file interpolation.c. References real.
00455 {
00456 int i;
00457 real ans = 0.0;
00458
00459 for (i = 3; i >= 0; i--)
00460 ans = t*(ans) + ((c[4*i + 3]*u + c[4*i + 2])*u + c[4*i + 1])*u + c[4*i];
00461 return ans;
00462 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 473 of file interpolation.c.
00478 {
00479 real d0, d1, d2, det;
00480 real a1, a2, a3;
00481
00482 d0 = sqrt(sq(x1 - x0) + sq(y1 - y0));
00483 d1 = d0 + sqrt(sq(x2 - x1) + sq(y2 - y1));
00484 d2 = d1 + sqrt(sq(x3 - x2) + sq(y3 - y2));
00485
00486 det = (d0*d1*(-d0 + d1)*d2*(-d0 + d2)*(-d1 + d2));
00487
00488 a1 = sq(d0)*(x0*(d1 - d2) + d2*x2 - d1*x3);
00489 a2 = sq(d1)*(x0*(d2 - d0) - d2*x1 + d0*x3);
00490 a3 = sq(d2)*(x0*(d0 - d1) + d1*x1 - d0*x2);
00491
00492 *ex = (a1 + a2 + a3) / det;
00493 *fx = - (d0*a1 + d1*a2 + d2*a3) / det;
00494 *gx = (sq(d0)*sq(d2)*(d2 - d0)*
00495 (d0*d0*d0*(x0 - x2) + d1*d1*d1*(x1 - x0)) +
00496 sq(d0)*sq(d1)*(d0 - d1)*
00497 (d0*d0*d0*(x0 - x3) + d2*d2*d2*(x1 - x0)))
00498 / (d0*d0*d0*det);
00499
00500 a1 = sq(d0)*(y0*(d1 - d2) + d2*y2 - d1*y3);
00501 a2 = sq(d1)*(y0*(d2 - d0) - d2*y1 + d0*y3);
00502 a3 = sq(d2)*(y0*(d0 - d1) + d1*y1 - d0*y2);
00503
00504 *ey = (a1 + a2 + a3) / det;
00505 *fy = - (d0*a1 + d1*a2 + d2*a3) / det;
00506 *gy = (sq(d0)*sq(d2)*(d2 - d0)*
00507 (d0*d0*d0*(y0 - y2) + d1*d1*d1*(y1 - y0)) +
00508 sq(d0)*sq(d1)*(d0 - d1)*
00509 (d0*d0*d0*(y0 - y3) + d2*d2*d2*(y1 - y0)))
00510 / (d0*d0*d0*det);
00511
00512 *rd2 = d2;
00513 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 194 of file interpolation.c. References polynom3::a, polynom3::b, polynom3::c, polynom3::d, POLYMAX, real, and sq. Referenced by closest_normal().
00196 {
00197 real a, b, c, d, e, f;
00198 real guess = 0.5*(t1 + t2), xh, xl, fl, fh, dxold, dx, f1, df, temp;
00199 int niter = 0;
00200
00201 a = px.a*px.b + py.a*py.b - px.b*xp - py.b*yp;
00202 b = sq(px.b) + sq(py.b) + 2.*(px.a*px.c + py.a*py.c - px.c*xp - py.c*yp);
00203 c = 3.*(px.b*px.c + py.b*py.c + px.a*px.d + py.a*py.d - px.d*xp - py.d*yp);
00204 d = 2.*(sq(px.c) + sq(py.c) + 2.*px.b*px.d + 2.*py.b*py.d);
00205 e = 5.*(px.c*px.d + py.c*py.d);
00206 f = 3.*(sq(px.d) + sq(py.d));
00207 #define dist(x) (a+(x)*(b+(x)*(c+(x)*(d+(x)*(e+(x)*f)))))
00208 #define dist1(x) (b+(x)*(2.*c+(x)*(3.*d+(x)*(4.*e+(x)*5.*f))))
00209
00210 fl = dist(t1);
00211 fh = dist(t2);
00212 if (fabs(fl) < prec) { *tmin = t1; return 1; }
00213 if (fabs(fh) < prec) { *tmin = t2; return 1; }
00214 if (fl < 0.0) {
00215 xl = t1;
00216 xh = t2;
00217 }
00218 else {
00219 xl = t2;
00220 xh = t1;
00221 }
00222 dxold = fabs(t2 - t1);
00223 dx = dxold;
00224 f1 = dist(guess);
00225 df = dist1(guess);
00226 while (niter++ < POLYMAX) {
00227 if ((((guess - xh)*df - f1)*((guess - xl)*df - f1) >= 0.0) ||
00228 (fabs(2.0*f1) > fabs(dxold*df))) {
00229 dxold = dx;
00230 dx = 0.5*(xh - xl);
00231 guess = xl + dx;
00232 }
00233 else {
00234 dxold = dx;
00235 dx = f1 / df;
00236 temp = guess;
00237 guess -= dx;
00238 }
00239 f1 = dist(guess);
00240 if (fabs(f1) < prec) { if (guess > t2 || guess < t1) return 0;
00241 *tmin = guess; return 1; }
00242 df = dist1(guess);
00243 if (f1 < 0.0)
00244 xl = guess;
00245 else
00246 xh = guess;
00247 }
00248 return 0;
00249
00250 #undef dist
00251 #undef dist1
00252 }
|
|
||||||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||
|
Definition at line 524 of file interpolation.c.
00525 {
00526 int i, j, k, i1;
00527 real a[NPTS][NPTS], lambda[NPTS], max, temp, xp2, yp2;
00528
00529 for (i = 0; i < NPTS; i++) {
00530 xp2 = xp[i]*xp[i]; yp2 = yp[i]*yp[i];
00531 a[i][0] = 1.0;
00532 a[i][1] = yp[i];
00533 a[i][2] = yp2;
00534 a[i][3] = xp[i];
00535 a[i][4] = xp[i]*yp[i];
00536 a[i][5] = xp[i]*yp2;
00537 a[i][6] = xp2;
00538 a[i][7] = xp2*yp[i];
00539 a[i][8] = xp2*yp2;
00540 }
00541
00542 /* solve the system of equations using gauss pivoting */
00543 for (j = 0; j <= NPTS - 2; j++) {
00544 i1 = j;
00545 max = fabs(a[j][j]);
00546 for (i = j + 1; i < NPTS; i++) {
00547 if (fabs(a[i][j]) >= max) {
00548 i1 = i;
00549 max = fabs(a[i][j]);
00550 }
00551 }
00552 for (k = j; k < NPTS; k++) {
00553 max = a[j][k]; a[j][k] = a[i1][k]; a[i1][k] = max;
00554 }
00555 max = z[j]; z[j] = z[i1]; z[i1] = max;
00556 for (i = j + 1; i < NPTS; i++) {
00557 if (a[j][j] == 0.) {
00558 printf("polynomial(%g, %g): no solution\n", x, y);
00559 for (i = 0; i < NPTS; i++)
00560 printf("%g %g\n", xp[i], yp[i]);
00561 return 0.0;
00562 }
00563 else
00564 temp = a[i][j] / a[j][j];
00565 for (k = j; k < NPTS; k++)
00566 a[i][k] -= temp * a[j][k];
00567 z[i] -= temp * z[j];
00568 }
00569 }
00570 if (a[NPTS - 1][NPTS - 1] == 0.) {
00571 printf("polynomial(%g, %g): no solution\n", x, y);
00572 for (i = 0; i < NPTS; i++)
00573 printf("%g %g\n", xp[i], yp[i]);
00574 return 0.0;
00575 }
00576 else
00577 lambda[NPTS - 1] = z[NPTS - 1] / a[NPTS - 1][NPTS - 1];
00578 for (i = NPTS - 2; i >= 0; i--) {
00579 lambda[i] = z[i];
00580 for (k = i + 1; k < NPTS; k++)
00581 lambda[i] -= a[i][k] * lambda[k];
00582 lambda[i] /= a[i][i];
00583 }
00584
00585 /* compute the value for coordinate (x,y)
00586 using the approximating function */
00587 xp2 = x*x;
00588 yp2 = y*y;
00589 return lambda[0] + lambda[1]*y + lambda[2]*yp2 +
00590 x*(lambda[3] + lambda[4]*y + lambda[5]*yp2) +
00591 xp2*(lambda[6] + lambda[7]*y + lambda[8]*yp2);
00592 }
|
|
||||||||||||||||||||||||
|
Definition at line 255 of file interpolation.c. References POLYMAX, real, splint, and splint1. Referenced by polyroots().
00256 {
00257 real guess = 0.5*(x1 + x2), xh, xl, fl, fh, dxold, dx, f, df, temp;
00258 int niter = 0;
00259
00260 /* Solving for the polynomial root using a combination of
00261 bissection and Newton's solver */
00262 fl = splint(poly, x1);
00263 fh = splint(poly, x2);
00264 if (fl == 0.0) {
00265 *x = x1;
00266 return 0;
00267 }
00268 if (fh == 0.0) {
00269 *x = x2;
00270 return 0;
00271 }
00272 if (fl < 0.0) {
00273 xl = x1;
00274 xh = x2;
00275 }
00276 else {
00277 xl = x2;
00278 xh = x1;
00279 }
00280 dxold = fabs(x2 - x1);
00281 dx = dxold;
00282 f = splint(poly, guess);
00283 df = splint1(poly, guess);
00284 while (niter++ < POLYMAX)
00285 {
00286 if ((((guess - xh)*df - f)*((guess - xl)*df - f) >= 0.0) ||
00287 (fabs(2.0*f) > fabs(dxold*df))) {
00288 dxold = dx;
00289 dx = 0.5*(xh - xl);
00290 guess = xl + dx;
00291 }
00292 else {
00293 dxold = dx;
00294 dx = f / df;
00295 temp = guess;
00296 guess -= dx;
00297 }
00298 f = splint(poly, guess);
00299 if (fabs(f) < prec) {
00300 *x = guess;
00301 return niter;
00302 }
00303 df = splint1(poly, guess);
00304 if (f < 0.0)
00305 xl = guess;
00306 else
00307 xh = guess;
00308 }
00309 fprintf(stderr, "WARNING polyroot(): Can not find a root in [%g,%g]\n", x1, x2);
00310 return POLYMAX + 1;
00311 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 314 of file interpolation.c. References polynom3::a, polynom3::b, polynom3::c, polynom3::d, polyroot(), real, splint, and sq. Referenced by cut_interface(), interface_fluxes(), interface_hfrac(), and interface_vfrac().
00315 {
00316 real delta, a, b, c, t1, t2, y1, y2, yt1, yt2;
00317 int n = 0;
00318
00319 poly.a -= y;
00320
00321 if (x1 > x2) { a = x1; x1 = x2; x2 = a; }
00322 y1 = splint(poly, x1);
00323 y2 = splint(poly, x2);
00324 /* find the extrema contained in {x1, x2} (if any) */
00325 a = 3.*poly.d;
00326 b = 2.*poly.c;
00327 c = poly.b;
00328 delta = sq(b) - 4.*a*c;
00329 if (delta < 0.0) /* no extremum anywhere -> one root maximum */
00330 t1 = t2 = x2 + 1.;
00331 else if (delta == 0.0) { /* one extremum -> two roots maximum */
00332 t1 = -b / (2.*a);
00333 t2 = x2 + 1.;
00334 if (t1 > t2) {a = t1; t1 = t2; t2 = a; }
00335 }
00336 else { /* two extrema -> three roots maximum */
00337 delta = sqrt(delta);
00338 t1 = (-b + delta) / (2.*a);
00339 t2 = (-b - delta) / (2.*a);
00340 if (t1 > t2) {a = t1; t1 = t2; t2 = a; }
00341 }
00342
00343 /* no extremum in [x1,x2] -> one root maximum in [x1,x2] */
00344 if ((t1 < x1 && t2 > x2) || t1 > x2 || t2 < x1) {
00345 if (y1*y2 <= 0.0) {
00346 polyroot(poly, x1, x2, x, prec);
00347 return 1;
00348 }
00349 return 0;
00350 }
00351
00352 yt1 = splint(poly, t1);
00353 yt2 = splint(poly, t2);
00354
00355 /* one extremum near x1 */
00356 if (t1 >= x1 && y1*yt1 <= 0.0)
00357 polyroot(poly, x1, t1, &(x[n++]), prec);
00358 if (t1 < x1 && t2 <= x2 && y1*yt2 <= 0.0)
00359 polyroot(poly, x1, t2, &(x[n++]), prec);
00360
00361 /* one extremum near x2 */
00362 if (t2 <= x2 && y2*yt2 <= 0.0)
00363 polyroot(poly, t2, x2, &(x[n++]), prec);
00364 if (t2 > x2 && t1 <= x2 && y2*yt1 <= 0.0)
00365 polyroot(poly, t1, x2, &(x[n++]), prec);
00366
00367 /* two extrema in [x1,x2] */
00368 if (t1 >= x1 && t2 <= x2 && yt1*yt2 <= 0.0)
00369 polyroot(poly, t1, t2, &(x[n++]), prec);
00370
00371 return n;
00372 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 105 of file interpolation.c. References polynom3::a, polynom3::b, polynom3::c, polynom3::d, real, real1D, and tridag(). Referenced by interface_splines().
00106 {
00107 real1D a, b, c, r, z, u;
00108 real alpha, gamma, fact, x1, x2, y1, y2, z1, z2, h;
00109 int i;
00110
00111 a = (real1D) malloc(n*sizeof(real));
00112 b = (real1D) malloc(n*sizeof(real));
00113 c = (real1D) malloc(n*sizeof(real));
00114 r = (real1D) malloc(n*sizeof(real));
00115 z = (real1D) malloc(n*sizeof(real));
00116 u = (real1D) malloc(n*sizeof(real));
00117
00118 b[0] = x[1] - x[n-1] + lp;
00119 c[0] = 0.5*(x[1] - x[0]);
00120 r[0] = 3.0*((y[1] - y[0])/(x[1] - x[0]) -
00121 (y[0] - y[n-1])/(x[0] - x[n-1] + lp));
00122 for (i = 1; i < n - 1; i++) {
00123 a[i] = 0.5*(x[i] - x[i-1]);
00124 b[i] = x[i+1] - x[i-1];
00125 c[i] = 0.5*(x[i+1] - x[i]);
00126
00127 r[i] = 3.0*((y[i+1] - y[i])/(x[i+1] - x[i]) -
00128 (y[i] - y[i-1])/(x[i] - x[i-1]));
00129 }
00130 a[n-1] = 0.5*(x[n-1] - x[n-2]);
00131 b[n-1] = x[0] - x[n-2] + lp;
00132 r[n-1] = 3.0*((y[0] - y[n-1])/(x[0] - x[n-1] + lp) -
00133 (y[n-1] - y[n-2])/(x[n-1] - x[n-2]));
00134 alpha = 0.5*(x[0] - x[n-1] + lp);
00135
00136 gamma = -b[0];
00137 b[0] -= gamma;
00138 b[n-1] -= alpha*alpha/gamma;
00139 tridag(a, b, c, r, d2y, n);
00140
00141 u[0] = gamma;
00142 u[n-1] = alpha;
00143 for (i = 1; i < n - 1; i++)
00144 u[i] = 0.0;
00145 tridag(a, b, c, u, z, n);
00146
00147 fact = (d2y[0] + alpha*d2y[n-1]/gamma)/(1.0 + z[0] + alpha*z[n-1]/gamma);
00148 for (i = 0; i < n; i++)
00149 d2y[i] -= fact*z[i];
00150
00151 /* Fill the coefficients of the interpolating polynomial for each segment */
00152 for (i = 0; i < n - 1; i++)
00153 {
00154 x1 = x[i]; y1 = y[i]; z1 = d2y[i];
00155 x2 = x[i+1]; y2 = y[i+1]; z2 = d2y[i+1];
00156 h = x2 - x1;
00157 poly[i].a = (x2*y1 - x1*y2)/h + h*(x1*z2 - x2*z1)/6.
00158 + (x2*x2*x2*z1 - x1*x1*x1*z2)/(6.*h);
00159 poly[i].b = (y2 - y1)/h + h*(z1 - z2)/6.
00160 + (x1*x1*z2 - x2*x2*z1)/(2.*h);
00161 poly[i].c = (x2*z1 - x1*z2)/(2.*h);
00162 poly[i].d = (z2 - z1)/(6.*h);
00163 }
00164 x1 = x[n - 1]; y1 = y[n - 1]; z1 = d2y[n - 1];
00165 x2 = x[0] + lp; y2 = y[0]; z2 = d2y[0];
00166 h = x2 - x1;
00167 poly[n - 1].a = (x2*y1 - x1*y2)/h + h*(x1*z2 - x2*z1)/6.
00168 + (x2*x2*x2*z1 - x1*x1*x1*z2)/(6.*h);
00169 poly[n - 1].b = (y2 - y1)/h + h*(z1 - z2)/6.
00170 + (x1*x1*z2 - x2*x2*z1)/(2.*h);
00171 poly[n - 1].c = (x2*z1 - x1*z2)/(2.*h);
00172 poly[n - 1].d = (z2 - z1)/(6.*h);
00173
00174 free(a); free(b); free(c); free(r); free(z); free(u);
00175 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 618 of file interpolation.c. References real, and trapzd().
00621 {
00622 int j;
00623 real s,olds;
00624
00625 olds = -1.0e30;
00626 for (j = 1; j <= 20; j++) {
00627 s = trapzd(func,px,py,a,b,j);
00628 if (fabs(s - olds) < eps*fabs(olds)) {
00629 printf("niter: %d\n", j);
00630 return s;
00631 }
00632 olds = s;
00633 }
00634 fprintf(stderr, "WARNING qtrap(): too many steps\n");
00635 return s;
00636 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 26 of file interpolation.c. References polynom3::a, polynom3::b, polynom3::c, polynom3::d, real, real1D, SP_DERIVATIVE, and SP_NATURAL. Referenced by interface_splines().
00028 {
00029 int i, k;
00030 real p, qn, sig, un, *u;
00031 real x1, y1, z1, x2, y2, z2, h;
00032
00033 u = (real1D) malloc(sizeof(real) * (n - 1));
00034
00035 switch(flag) {
00036 case SP_NATURAL: d2y[0] = u[0] = 0.0; break;
00037 case SP_DERIVATIVE:
00038 d2y[0] = -0.5;
00039 u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
00040 break;
00041 default:
00042 fprintf(stderr, "spline(): Bad flag parameter\n");
00043 exit(1);
00044 }
00045
00046 for (i = 1; i < n - 1; i++) {
00047 sig = (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
00048 p = sig * d2y[i-1] + 2.0;
00049 d2y[i] = (sig - 1.0) / p;
00050 u[i] = (y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1])
00051 / (x[i] - x[i-1]);
00052 u[i] = (6.0 * u[i] / (x[i+1] - x[i-1]) - sig * u[i-1]) / p;
00053 }
00054
00055 switch(flag) {
00056 case SP_NATURAL: d2y[n - 1] = 0.0;
00057 break;
00058 case SP_DERIVATIVE:
00059 qn = 0.5;
00060 un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])
00061 / (x[n-1] - x[n-2]));
00062 d2y[n - 1] = (un - qn * u[n-2]) / (qn * d2y[n-2] + 1.0);
00063 break;
00064 }
00065
00066 for (k = n - 2; k >= 0; k--)
00067 d2y[k] = d2y[k] * d2y[k+1] + u[k];
00068 free(u);
00069
00070 /* Fill the coefficients of the interpolating polynomial for each segment */
00071 for (i = 0; i < n - 1; i++)
00072 {
00073 x1 = x[i]; y1 = y[i]; z1 = d2y[i];
00074 x2 = x[i+1]; y2 = y[i+1]; z2 = d2y[i+1];
00075 h = x2 - x1;
00076 poly[i].a = (x2*y1 - x1*y2)/h + h*(x1*z2 - x2*z1)/6.
00077 + (x2*x2*x2*z1 - x1*x1*x1*z2)/(6.*h);
00078 poly[i].b = (y2 - y1)/h + h*(z1 - z2)/6.
00079 + (x1*x1*z2 - x2*x2*z1)/(2.*h);
00080 poly[i].c = (x2*z1 - x1*z2)/(2.*h);
00081 poly[i].d = (z2 - z1)/(6.*h);
00082 }
00083 }
|
|
||||||||||||||||||||
|
Definition at line 178 of file interpolation.c. References polynom3::a, polynom3::b, polynom3::c, polynom3::d, real, and real1D. Referenced by interface_write_splines().
00179 {
00180 int klo, khi, k;
00181
00182 /* then search for the nearest lower point */
00183 klo = 0;
00184 khi = n - 1;
00185 while (khi - klo > 1) {
00186 k = (khi + klo) >> 1;
00187 if (xa[k] > x) khi = k;
00188 else klo = k;
00189 }
00190 return poly[klo].a + x*(poly[klo].b + x*(poly[klo].c + x*poly[klo].d));
00191 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 596 of file interpolation.c. References real. Referenced by qtrap().
00599 {
00600 real x,tnm,sum,del;
00601 static real s;
00602 int it,j;
00603
00604 if (n == 1) {
00605 return (s=0.5*(b-a)*((*func)(px,py,a)+(*func)(px,py,b)));
00606 } else {
00607 for (it=1,j=1;j<n-1;j++) it <<= 1;
00608 tnm=it;
00609 del=(b-a)/tnm;
00610 x=a+0.5*del;
00611 for (sum=0.0,j=1;j<=it;j++,x+=del) sum += (*func)(px,py,x);
00612 s=0.5*(s+(b-a)*sum/tnm);
00613 return s;
00614 }
00615 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 86 of file interpolation.c. Referenced by Pspline().
00087 {
00088 int j;
00089 real bet, *gam;
00090
00091 gam = (real *)malloc(n*sizeof(real));
00092 bet = b[0];
00093 u[0] = r[0]/bet;
00094 for (j = 1; j < n; j++) {
00095 gam[j] = c[j-1]/bet;
00096 bet = b[j] - a[j]*gam[j];
00097 u[j] = (r[j] - a[j]*u[j-1])/bet;
00098 }
00099 for (j = n - 2; j >= 0; j--)
00100 u[j] -= gam[j+1]*u[j+1];
00101 free(gam);
00102 }
|
1.2.18