backend/nsl/nsl_interp.c

0001 /*
0002     File                 : nsl_interp.c
0003     Project              : LabPlot
0004     Description          : NSL interpolation functions
0005     --------------------------------------------------------------------
0006     SPDX-FileCopyrightText: 2016 Stefan Gerlach <stefan.gerlach@uni.kn>
0007
0008     SPDX-License-Identifier: GPL-2.0-or-later
0009 */
0010
0011 #include "nsl_interp.h"
0012 #include "nsl_common.h"
0013
0014 const char* nsl_interp_type_name[] = {i18n("Linear"),
0015                                       i18n("Polynomial"),
0016                                       i18n("Cubic Spline (Natural)"),
0017                                       i18n("Cubic Spline (Periodic)"),
0018                                       i18n("Akima-Spline (Natural)"),
0019                                       i18n("Akima-Spline (Periodic)"),
0020                                       i18n("Steffen-Spline"),
0021                                       i18n("Cosine"),
0022                                       i18n("Exponential"),
0023                                       i18n("Piecewise Cubic Hermite (PCH)"),
0024                                       i18n("Rational Functions")};
0025 const char* nsl_interp_pch_variant_name[] = {i18n("Finite Differences"), i18n("Catmull-Rom"), i18n("Cardinal"), i18n("Kochanek-Bartels (TCB)")};
0026 const char* nsl_interp_evaluate_name[] = {i18n("Function"), i18n("Derivative"), i18n("Second Derivative"), i18n("Integral")};
0027
0028 int nsl_interp_ratint(const double* x, const double* y, int n, double xn, double* v, double* dv) {
0029     int i, a = 0, b = n - 1;
0030     while (b - a > 1) { /* find interval using bisection */
0031         int j = (int)floor((a + b) / 2.);
0032         if (x[j] > xn)
0033             b = j;
0034         else
0035             a = j;
0036     }
0037
0038     int ns = a; /* nearest index */
0039     if (fabs(xn - x[a]) > fabs(xn - x[b]))
0040         ns = b;
0041
0042     if (xn == x[ns]) { /* exact point */
0043         *v = y[ns];
0044         *dv = 0;
0045         return 1;
0046     }
0047
0048     double* c = (double*)malloc(n * sizeof(double));
0049     double* d = (double*)malloc(n * sizeof(double));
0050     for (i = 0; i < n; i++)
0051         c[i] = d[i] = y[i];
0052
0053     *v = y[ns--];
0054
0055     double t, dd;
0056     int m;
0057     for (m = 1; m < n; m++) {
0058         for (i = 0; i < n - m; i++) {
0059             t = (x[i] - xn) * d[i] / (x[i + m] - xn);
0060             dd = t - c[i + 1];
0061             if (dd == 0.0) /* pole */
0062                 dd += DBL_MIN;
0063             dd = (c[i + 1] - d[i]) / dd;
0064             d[i] = c[i + 1] * dd;
0065             c[i] = t * dd;
0066         }
0067
0068         *dv = (2 * (ns + 1) < (n - m) ? c[ns + 1] : d[ns--]);
0069         *v += *dv;
0070     }
0071
0072     free(c);
0073     free(d);
0074
0075     return 0;
0076 }