+++ /dev/null
-/* spline.c -- Do spline interpolation. */
-
-/******************************************************************************
- David L. Fox
- 2140 Braun Dr.
- Golden, CO 80401
-
-This program has been placed in the public domain by its author.
-
-Version Date Change
-1.1 17 Dec 1985 Modify getdata() to realloc() more memory if needed.
-1.0 14 May 1985
-
-spline [options] [file]
-
-Spline reads pairs of numbers from file (or the standard input, if file is missing).
-The pairs are interpreted as abscissas and ordinates of a function. The output of
-spline consists of similar pairs generated from the input by interpolation with
-cubic splines. (See R. W. Hamming, Numerical Methods for Scientists and Engineers,
-2nd ed., pp. 349ff.) There are sufficiently many points in the output to appear
-smooth when plotted (e.g., by graph(1)). The output points are approximately evenly
-spaced and include the input points.
-
-There may be one or more options of the form: -o [argument [arg2] ].
-
-The available options are:
-
--a Generate abscissa values automatically. The input consists of a list of
- ordinates. The generated abscissas are evenly spaced with spacing given by
- the argument, or 1 if the next argument is not a number.
-
--f Set the first derivative of the spline function at the left and right end
- points to the first and second arguments following -f. If only one numerical
- argument follows -f then that value is used for both left and right end points.
-
--k The argument of the k option is the value of the constant used to calculate
- the boundary values by y''[0] = ky''[1] and y''[n] = ky''[n-1]. The default
- value is k = 0.
-
--m The argument gives the maximum number of input data points. The default is 1000.
- If the input contains more than this many points additional memory will be
- allocated. This option may be used to slightly increase efficiency for large
- data sets or reduce the amount of memory required for small data sets.
-
--n The number of output points is given by the argument. The default value is 100.
-
--p The splines used for interpolation are forced to be periodic, i.e. y'[0] = y'[n].
- The first and last ordinates should be equal.
-
--s Set the second derivative of the spline function at the left and right end
- points to the first and second arguments following -s. If only one numerical
- argument follows -s then that value is used for both left and right end points.
-
--x The argument (and arg2, if present) are the lower (and upper) limits on the
- abscissa values in the output. If the x option is not specified these limits
- are calculated from the data. Automatic abscissas start at the lower limit
- (default 0).
-
-The data need not be monotonic in x but all the x values must be distinct.
-Non-numeric data in the input is ignored.
-******************************************************************************/
-
-#include <a:stdio.h>
-#include <ctype.h>
-
-/* The constant DOUBLE is a compile time switch.
- If #define DOUBLE appears here double pecision variables are
- used to store all data and parameters.
- Otherwise, float variables are used for the data.
- For most smoothing and and interpolation applications single
- precision is more than adequate.
- Double precision is used to solve the system of linear equations
- in either case.
-*/
-/* #define DOUBLE */
-
-#ifdef DOUBLE
-# define double real; /* Type used for data storage. */
-# define IFMT "%F" /* Input format. */
-# define OFMT "%18.14g %18.14g\n" /* Output format. */
-#else
-# define float real; /* Type used for data storage. */
-# define IFMT "%f" /* Input format. */
-# define OFMT "%8.5g %8.5g\n" /* Output format. */
-#endif
-
-/* Numerical constants: These may be machine and/or precision dependent. */
-#define HUGE (1.e38) /* Largest floating point number. */
-#define EPS 1.e-5 /* Test for zero when solving linear equations. */
-
-/* Default parameters */
-#define NPTS 1000
-#define NINT 100
-#define DEFSTEP 1. /* Default step size for automatic abcissas. */
-#define DEFBVK 0. /* Default boundary value constant. */
-
-/* Boundary condition types. */
-#define EXTRAP 0 /* Extrapolate second derivative:
- y''(0) = k * y''(1), y''(n) = k * y''(n-1) */
-#define FDERIV 1 /* Fixed first derivatives y'(0) and y'(n). */
-#define SDERIV 2 /* Fixed second derivatives y''(0) and y''(n). */
-#define PERIOD 3 /* Periodic: derivatives equal at end points. */
-
-/* Token types for command line processing. */
-#define OPTION 1
-#define NUMBER 2
-#define OTHER 3
-
-/* Define error numbers. */
-#define MEMERR 1
-#define NODATA 2
-#define BADOPT 4
-#define BADFILE 5
-#define XTRAARG 6
-#define XTRAPTS 7
-#define SINGERR 8
-#define DUPX 9
-#define BADBC 10
-#define RANGERR 11
-
-/* Constants and flags are global. */
-int aflag = FALSE; /* Automatic abscissa flag. */
-real step = DEFSTEP; /* Spacing for automatic abscissas. */
-int bdcode = EXTRAP; /* Type of boundary conditions:
- 0 extrapolate f'' with coeficient k.
- 1 first derivatives given
- 2 second derivatives given
- 3 periodic */
-real leftd = 0., /* Boundary values of derivatives. */
- rightd = 0.;
-real k = DEFBVK; /* Boundary value constant. */
-int rflag = 0; /* 0: take range from data, 1: min. given, 2: min. & max. given. */
-real xmin = HUGE, xmax; /* Output range. */
-unsigned nint = NINT; /* Number of intervals in output. */
-unsigned maxpts = NPTS; /* Maximun number of points. */
-unsigned nknots = 0; /* Number of input points. */
-int sflag = FALSE; /* If TRUE data must be sorted. */
-char *datafile; /* Name of data file */
-
-int xcompare(); /* Function to compare x values of two data points. */
-
-struct pair {
- real x, y;
- } *data; /* Points to list of data points. */
-
-struct bandm {
- double diag, upper, rhs;
- } *eqns; /* Points to elements of band matrix used to calculate
- coefficients of the splines. */
-
-main(argc, argv)
-int argc;
-char **argv;
-{
- setup(argc, argv);
- if (nknots == 1) { /* Cannot interpolate with one data point. */
- printf(OFMT,data->x, data->y);
- exit(0);
- }
- if (sflag) /* Sort data if needed. */
- qsort(data, nknots, sizeof(struct pair), xcompare);
- calcspline();
- interpolate();
-}
-
-/* xcompare -- Compare abcissas of two data points (for qsort). */
-xcompare(arg1, arg2)
-struct pair *arg1, *arg2;
-{
- if (arg1->x > arg2->x)
- return(1);
- else if (arg1->x < arg2->x)
- return(-1);
- else
- return(0);
-}
-
-/* error -- Print error message and abort fatal errors. */
-error(errno, auxmsg)
-int errno;
-char *auxmsg;
-{ char *msg;
- int fatal, usemsg;
- static char *usage =
- "usage: spline [options] [file]\noptions:\n",
- *options = "-a spacing\n-k const\n-n intervals\n-m points\n-p\n-x xmin xmax\n";
-
- fatal = TRUE; /* Default is fatal error. */
- usemsg = FALSE; /* Default no usage message. */
- fprintf(stderr, "spline: ");
- switch(errno) {
- case MEMERR:
- msg = "not enough memory for %u data points\n";
- break;
- case NODATA:
- msg = "data file %s empty\n";
- break;
- case BADOPT:
- msg = "unknown option: %c\n";
- usemsg = TRUE;
- break;
- case BADFILE:
- msg = "cannot open file: %s\n";
- break;
- case XTRAARG:
- msg = "extra argument ignored: %s\n";
- fatal = FALSE;
- usemsg = TRUE;
- break;
- case XTRAPTS:
- fatal = FALSE;
- msg = "%s";
- break;
- case DUPX:
- msg = "duplicate abcissa value: %s\n";
- break;
- case RANGERR:
- msg = "xmax < xmin not allowed %s\n";
- break;
- /* The following errors "can't happen." */
- /* If they occur some sort of bug is indicated. */
- case SINGERR:
- msg = "singular matrix encountered %s\n";
- break;
- case BADBC:
- msg = "internal error: bad boundary value code %d\n";
- break;
- default:
- fprintf(stderr, "unknown error number: %d\n", errno);
- exit(1);
- }
- fprintf(stderr, msg, auxmsg);
- if (usemsg)
- fprintf(stderr,"%s%s", usage, options);
- if (fatal)
- exit(1);
-}
-
-/* setup -- Initalize all constants and read data. */
-setup(argc, argv)
-int argc;
-char **argv;
-{ char *malloc();
- FILE fdinp, /* Source of input. */
- doarg();
-
- fdinp = doarg(argc, argv); /* Process command line arguments. */
-
- /* Allocate memory for data and band matrix of coefficients. */
- if ((data = malloc(maxpts*sizeof(struct pair))) == NULL) {
- error(MEMERR, (char *)maxpts);
- }
-
- getdata(fdinp); /* Read data from fdinp. */
- if (fdinp != stdin)
- fclose(fdinp); /* Close input data file. */
- if (nknots == 0) {
- error(NODATA, datafile);
- }
- /* Allocate memory for calculation of spline coefficients. */
- if ((eqns = malloc((nknots+1)*sizeof(struct bandm))) == NULL) {
- error(MEMERR, (char *)nknots);
- }
-}
-
-/* doarg -- Process arguments. */
-FILE
-doarg(argc, argv)
-int argc;
-char **argv;
-{ int i, type;
- double atof();
- char *s, str[15];
- FILE fdinp;
-
- s = argv[i=1];
- type = gettok(&i, &s, argv, str);
- do {
- if (type == OPTION) {
- switch(*str) {
- case 'a': /* Automatic abscissa. */
- aflag = TRUE;
- rflag = rflag < 2 ? 1 : rflag;
- if (xmin == HUGE) /* Initialize xmin, if needed. */
- xmin = 0.;
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- if ((step = atof(str)) <= 0.)
- error(RANGERR,"");
- type = gettok(&i, &s, argv, str);
- }
- break;
- case 'f': /* Fix first derivative at boundaries. */
- bdcode = FDERIV;
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- leftd = atof(str);
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- rightd = atof(str);
- type = gettok(&i, &s, argv, str);
- }
- else {
- rightd = leftd;
- }
- }
- break;
- case 'k': /* Set boundary value constant. */
- bdcode = EXTRAP;
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- k = atof(str);
- type = gettok(&i, &s, argv, str);
- }
- break;
- case 'm': /* Set number of intervals in output. */
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- maxpts = (unsigned)atoi(str);
- type = gettok(&i, &s, argv, str);
- }
- break;
- case 'n': /* Set number of intervals in output. */
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- nint = (unsigned)atoi(str);
- type = gettok(&i, &s, argv, str);
- }
- break;
- case 'p': /* Require periodic interpolation function. */
- bdcode = PERIOD;
- type = gettok(&i, &s, argv, str);
- break;
- case 's': /* Fix second derivative at boundaries. */
- bdcode = SDERIV;
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- leftd = atof(str);
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- rightd = atof(str);
- type = gettok(&i, &s, argv, str);
- }
- else {
- rightd = leftd;
- }
- }
- break;
- case 'x': /* Set range of x values. */
- rflag = 1;
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- xmin = atof(str);
- if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
- xmax = atof(str);
- rflag = 2;
- type = gettok(&i, &s, argv, str);
- if (xmax <= xmin)
- error(RANGERR, "");
- }
- }
- break;
- default:
- error(BADOPT, (char *)argv[i][1]);
- }
- }
- else {
- if (argc > i) {
- datafile = argv[i];
- if ((fdinp = fopen(argv[i++], "r")) == NULL) {
- error(BADFILE, argv[i-1]);
- }
- if (argc > i)
- error(XTRAARG, argv[i]);
- }
- else
- fdinp = stdin;
- break;
- }
- } while (i < argc);
- return fdinp;
-}
-
-/* gettok -- Get one token from command line, return type. */
-gettok(indexp, locp, argv, str)
-int *indexp; /* Pointer to index in argv array. */
-char **locp; /* Pointer to current location in argv[*indexp]. */
-char **argv;
-char *str;
-{ char *s;
- char *strcpy(), *strchr();
- int type;
-
- s = *locp;
- while (isspace(*s) || *s == ',')
- ++s;
- if (*s == '\0') /* Look at next element in argv. */
- s = argv[++*indexp];
- if (*s == '-' && isalpha(s[1])) {
- /* Found an option. */
- *str = *++s;
- str[1] = '\0';
- ++s;
- type = OPTION;
- }
- else if (is_number(s)) {
- while (is_number(s))
- *str++ = *s++;
- *str = '\0';
- type = NUMBER;
- }
- else {
- strcpy(str, s);
- s = strchr(s, '\0');
- type = OTHER;
- }
- *locp = s;
- return(type);
-}
-
-/* is_number -- Return TRUE if argument is the ASCII representation of a number. */
-is_number(string)
-char *string;
-{
- if (isdigit(*string) ||
- *string == '.' ||
- *string == '+' ||
- (*string == '-' && (isdigit(string[1]) || string[1] == '.')))
- return(TRUE);
- else
- return(FALSE);
-}
-
-/* getdata -- Read data points from fdinp. */
-getdata(fdinp)
-FILE fdinp;
-{ int n, i;
- real lastx, min, max;
- struct pair *dp;
- char msg[60], *realloc();
-
- nknots = 0;
- lastx = -HUGE;
- dp = data; /* Point to head of list of data. */
- min = HUGE;
- max = -HUGE;
- do {
- if (aflag) { /* Generate abcissa. */
- dp->x = xmin + nknots*step;
- n = 0;
- }
- else { /* Read abcissa. */
- while ((n = (fscanf(fdinp,IFMT,&dp->x))) == 0)
- ; /* Skip non-numeric input. */
- }
- if (n == 1) {
- if (min > dp->x)
- min = dp->x;
- if (max < dp->x)
- max = dp->x;
- if (lastx > dp->x) { /* Check for monotonicity. */
- sflag = TRUE; /* Data must be sorted. */
- }
- lastx = dp->x;
- }
- /* Read ordinate. */
- while ((n = (fscanf(fdinp, IFMT, &dp->y))) == 0)
- ; /* Skip non-numeric input. */
- if (++nknots >= maxpts) { /* Too many points, allocate more memory. */
- if ((data = realloc(data, (maxpts *= 2)*sizeof(struct pair))) == NULL) {
- error(MEMERR, (char *)maxpts);
- }
- dp = data + nknots;
- sprintf(msg, "more than %d input points, more memory allocated\n",
- maxpts/2);
- error(XTRAPTS,msg);
- }
- else
- ++dp;
- } while (n == 1);
- --nknots;
- if (aflag) { /* Compute maximum ordinate. */
- max = xmin + (nknots-1)*step;
- }
- if (rflag < 2) {
- xmax = max;
- if (rflag < 1)
- xmin = min;
- }
-}
-
-/* calcspline -- Calculate coeficients of spline functions. */
-calcspline()
-{
- calccoef();
- solvband();
-}
-
-/* calccoef -- Calculate coefficients of linear equations determining spline functions. */
-calccoef()
-{ int i, j;
- struct bandm *ptr, tmp1, tmp2;
- double h, h1;
- char str[10];
-
- ptr = eqns;
- /* Initialize first row of matrix. */
- if ((h1 = data[1].x - data[0].x) == 0.) {
- sprintf(str, "%8.5g", data[0].x);
- error(DUPX, str);
- }
- switch(bdcode) { /* First equation depends on boundary conditions. */
- case EXTRAP:
- ptr->upper = -k;
- ptr->diag = 1.;
- ptr->rhs = 0.;
- break;
- case FDERIV: /* Fixed first derivatives at ends. */
- ptr->upper = 1.;
- ptr->diag = 2.;
- h = data[1].x - data[0].x;
- ptr->rhs = 6.*((data[1].y - data[0].y)/(h*h) - leftd/h);
- break;
- case SDERIV:
- ptr->upper = 0.;
- ptr->diag = 1.;
- ptr->rhs = leftd;
- break;
- case PERIOD: /* Periodic splines. */
- ptr->upper = 1.;
- ptr->diag = 2.;
- break;
- default:
- error(BADBC, (char *) bdcode);
- }
- ++ptr;
-
- /* Initialize rows 1 to n-1, sub-diagonal band is assumed to be 1. */
- for (i=1; i < nknots-1; ++i, ++ptr) {
- h = h1;
- if ((h1 = data[i+1].x - data[i].x) == 0.) {
- sprintf(str, "%8.5g", data[i].x);
- error(DUPX, str);
- }
- ptr->diag = 2.*(h + h1)/h;
- ptr->upper = h1/h;
- ptr->rhs = 6.*((data[i+1].y-data[i].y)/(h*h1) -
- (data[i].y - data[i-1].y)/(h*h));
- }
- /* Initialize last row. */
- switch(bdcode) {
- case EXTRAP: /* Standard case. */
- ptr->diag = 1.;
- ptr->upper = -k; /* Upper holds actual sub-diagonal value. */
- ptr->rhs = 0.;
- break;
- case FDERIV: /* Fixed first derivatives at ends. */
- ptr->upper = 1.;
- ptr->diag = 2.;
- h = data[nknots-1].x - data[nknots-2].x;
- ptr->rhs = 6.*((data[nknots-2].y - data[nknots-1].y)/(h*h) + rightd/h);
- break;
- case SDERIV:
- ptr->upper = 0.;
- ptr->diag = 1.;
- ptr->rhs = rightd;
- break;
- case PERIOD: /* Use periodic boundary conditions. */
- /* First and last row are not in tridiagonal form. */
- h = data[1].x - data[0].x;
- h1 = data[nknots-1].x - data[nknots-2].x;
- ptr->diag = 1.;
- ptr->upper = 0.;
- tmp1.diag = -1.;
- tmp1.upper = 0;
- tmp1.rhs = 0.;
- tmp2.diag = 2.*h1/h;
- tmp2.upper = h1/h;
- tmp2.rhs = 6.*((data[1].y - data[0].y)/(h*h) -
- (data[nknots-1].y - data[nknots-2].y)/(h1*h));
- /* Transform periodic boundary equations to tri-diagonal form. */
- for (i = 1; i < nknots - 1; ++i) {
- tmp1.upper /= tmp1.diag;
- tmp1.rhs /= tmp1.diag;
- ptr->diag /= tmp1.diag;
- ptr->upper /= tmp1.diag;
- tmp1.diag = tmp1.upper - eqns[i].diag;
- tmp1.upper = -eqns[i].upper;
- tmp1.rhs -= eqns[i].rhs;
- tmp2.upper /= tmp2.diag;
- tmp2.rhs /= tmp2.diag;
- eqns->diag /= tmp2.diag;
- eqns->upper /= tmp2.diag;
- tmp2.diag = tmp2.upper - eqns[nknots-1-i].diag/eqns[nknots-1-i].upper;
- tmp2.upper = -1./eqns[nknots-1-i].upper;
- tmp2.rhs -= eqns[nknots-1-i].rhs/eqns[nknots-1-i].upper;
- }
- /* Add in remaining terms of boundary condition equation. */
- ptr->upper += tmp1.diag;
- ptr->diag += tmp1.upper;
- ptr->rhs = tmp1.rhs;
- eqns->diag += tmp2.upper;
- eqns->upper += tmp2.diag;
- eqns->rhs = tmp2.rhs;
- break;
- default:
- error(BADBC, (char *) bdcode);
- }
-}
-
-/* solvband -- Solve band matrix for spline functions. */
-solvband()
-{ int i, flag;
- struct bandm *ptr;
- double k1;
- double fabs();
- int fcompare();
-
- ptr = eqns;
- flag = FALSE;
- /* Make a pass to triangularize matrix. */
- for (i=1; i < nknots - 1; ++i, ++ptr) {
- if (fabs(ptr->diag) < EPS) {
- /* Near zero on diagonal, pivot. */
- if (fabs(ptr->upper) < EPS)
- error(SINGERR, "");
- flag = TRUE;
- ptr->diag = i; /* Keep row index in diag.
- Actual value of diag is always 1. */
- if (i == nknots - 2) {
- flag = 2;
- /* Exchange next to last and last rows. */
- k1 = ptr->rhs/ptr->upper;
- if (fabs((ptr+1)->upper) < EPS)
- error(SINGERR, "");
- ptr->rhs = (ptr+1)->rhs/(ptr+1)->upper;
- ptr->upper = (ptr+1)->diag/(ptr+1)->upper;
- (ptr+1)->upper = 0.;
- (ptr+1)->rhs = k1;
- }
- else {
- ptr->rhs = (ptr+1)->rhs - (k1 = ptr->rhs/ptr->upper)*(ptr+1)->diag;
- ptr->upper = (ptr+1)->upper; /* This isn't super-diagonal element
- but rather one to its right.
- Must watch for this when
- back substituting. */
- (++ptr)->diag = i-1;
- ++i;
- ptr->upper = 0;
- ptr->rhs = k1;
- (ptr+1)->rhs -= ptr->rhs;
- }
- }
- else {
- ptr->upper /= ptr->diag;
- ptr->rhs /= ptr->diag;
- ptr->diag = i-1; /* Used to reorder solution if needed. */
- (ptr+1)->diag -= ptr->upper;
- (ptr+1)->rhs -= ptr->rhs;
- }
- }
- /* Last row is a special case. */
- if (flag != 2) {
- /* If flag == 2 last row is already computed. */
- ptr->upper /= ptr->diag;
- ptr->rhs /= ptr->diag;
- ptr->diag = ptr - eqns;
- ++ptr;
- ptr->diag -= (ptr-1)->upper * ptr->upper;
- ptr->rhs -= (ptr-1)->rhs * ptr->upper;
- ptr->rhs /= ptr->diag;
- ptr->diag = ptr - eqns;
- }
-
- /* Now make a pass back substituting for solution. */
- --ptr;
- for ( ; ptr >= eqns; --ptr) {
- if ((int)ptr->diag != ptr - eqns) {
- /* This row and one above have been exchanged in a pivot. */
- --ptr;
- ptr->rhs -= ptr->upper * (ptr+2)->rhs;
- }
- else
- ptr->rhs -= ptr->upper * (ptr+1)->rhs;
- }
- if (flag) { /* Undo reordering done by pivoting. */
- qsort(eqns, nknots, sizeof(struct bandm), fcompare);
- }
-}
-
-/* fcompare -- Compare two floating point numbers. */
-fcompare(arg1, arg2)
-real *arg1, *arg2;
-{
- if (arg1 > arg2)
- return(1);
- else if (arg1 < arg2)
- return(-1);
- else
- return(0);
-}
-
-/* interpolate -- Do spline interpolation. */
-interpolate()
-{ int i;
- struct pair *dp;
- struct bandm *ep;
- real h, xi, yi, hi, xu, xl, limit;
-
- h = (xmax - xmin)/nint;
- ep = eqns;
- dp = data + 1;
- for (xi = xmin; xi < xmax + 0.25*h; xi += h) {
- while (dp->x < xi && dp < data + nknots - 1) {
- ++dp; /* Skip to correct interval. */
- ++ep;
- }
- if (dp < data + nknots - 1 && dp->x < xmax)
- limit = dp->x;
- else
- limit = xmax;
- for ( ; xi < limit - 0.25*h; xi += h) {
- /* Do interpolation. */
- hi = dp->x - (dp-1)->x;
- xu = dp->x - xi;
- xl = xi - (dp-1)->x;
- yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
- (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
- printf(OFMT, xi, yi);
- }
- if (limit != dp->x) { /* Interpolate. */
- hi = dp->x - (dp-1)->x;
- xu = dp->x - xmax;
- xl = xmax - (dp-1)->x;
- yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
- (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
- printf(OFMT, xmax, yi);
- }
- else { /* Print knot. */
- printf(OFMT, dp->x, dp->y);
- xi = dp->x;
- }
- }
-}
projects / thirdparty / coreutils.git / commitdiff
