home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Reference Guide / C-C++ Interactive Reference Guide.iso / c_ref / csource3 / 167_01 / spline.c < prev    next >
Encoding:
Text File  |  1985-11-15  |  17.6 KB  |  735 lines
  1. /* 
  2. HEADER:     CUG
  3. TITLE:        SPLINE.C - Interpolate Smooth Curve;
  4. VERSION:    3.00;
  5. DATE:        09/26/86;
  6. DESCRIPTION:    "SPLINE takes pairs of numbers from the standard
  7.         input as abscissae and ordinates of a function.
  8.         It produces a similar set, which is approximately
  9.         equally spaced and includes the input set, on the
  10.         standard output. The cubic spline output has two
  11.         continuous derivatives and sufficiently many
  12.         points to look smooth when plotted.";
  13. KEYWORDS:    spline, graphics, plot, filter, UNIX;
  14. SYSTEM:        ANY;
  15. FILENAME:    SPLINE.DOC;
  16. WARNINGS:    NONE;
  17. CRC:        xxxx
  18. SEE-ALSO:    SPLINE.C;
  19. AUTHORS:    Ian Ashdown - byHeart Software;
  20. COMPILERS:    ANY;
  21. REFERENCES:    AUTHORS: Bell Laboratories;
  22.         TITLE:     UNIX Programmer's Manual, Volume 1
  23.              p. 145, Holt, Rinehart and Winston;
  24.         AUTHORS: R.W. Hamming;
  25.         TITLE:     Numerical Methods for Scientists and
  26.              Engineers, 2nd Edition
  27.              pp. 349 - 355, McGraw-Hill (1973);
  28.         AUTHORS: Forsythe, Malcolm & Moler;
  29.         TITLE:     Computer Methods for Mathematical
  30.              Computations
  31.              pp. 70 - 83, Prentice-Hall;
  32. ENDREF
  33. */
  34.  
  35. /*-------------------------------------------------------------*/
  36.  
  37. /* SPLINE.C - Interpolate Smooth Curve
  38.  *
  39.  * Version 3.00        September 26th, 1986
  40.  *
  41.  * Modifications:
  42.  *
  43.  *   V1.00 (85/11/01)    - beta test release
  44.  *   V2.00 (85/12/25)    - general revision
  45.  *   V2.01 (86/09/13)    - corrected command-line option bug for
  46.  *              Microsoft C compiler
  47.  *   V3.00 (86/09/26)    - added non-uniform abscissa spacing
  48.  *              capability 
  49.  *
  50.  * Copyright 1985,1986:    Ian Ashdown
  51.  *            byHeart Software
  52.  *            620 Ballantree Road
  53.  *            West Vancouver, B.C.
  54.  *            Canada V7S 1W3
  55.  *
  56.  * This program may be copied for personal, non-commercial use
  57.  * only, provided that the above copyright notice is included in
  58.  * all copies of the source code. Copying for any other use
  59.  * without previously obtaining the written permission of the
  60.  * author is prohibited.
  61.  *
  62.  * Synopsis:    SPLINE [option] ...
  63.  *
  64.  * Description:    SPLINE takes pairs of numbers from the standard
  65.  *        input as abscissae and ordinates of a function.
  66.  *        (A minimum of four pairs is required.) It
  67.  *        produces a similar set, which is approximately
  68.  *        equally spaced and includes the input set, on the
  69.  *        standard output. The cubic spline output (R.W.
  70.  *        Hamming, "Numerical Methods for Scientists and
  71.  *        Engineers", 2nd ed. 349ff) has two continuous
  72.  *        derivatives and sufficiently many points to look
  73.  *        smooth when plotted.
  74.  *
  75.  *        The following options are recognized, each as a
  76.  *        separate argument:
  77.  *
  78.  *        -a  Supply abscissae automatically (they are
  79.  *            missing from the input); spacing is given by
  80.  *            the next argument or is assumed to be 1 if
  81.  *            next argument is not a number.
  82.  *
  83.  *        -k  The constant "k" is used in the boundary
  84.  *            value computation
  85.  *
  86.  *            y " = ky " , y " = ky "
  87.  *             0    1     n         n-1
  88.  *
  89.  *            is set by the next argument. By default,
  90.  *            k = 0. A value of k = 0.5 often results in a
  91.  *            smoother curve at the endpoints than the
  92.  *            default value. Negative values for k are not
  93.  *            allowed. Cannot be used with -p option.
  94.  *
  95.  *        -n  Next argument (which must be an integer)
  96.  *            specifies the number of intervals that are to
  97.  *            occur between the lower and upper "x" limits.
  98.  *            If -n option is not given, default spacing is
  99.  *            100 intervals.
  100.  *
  101.  *        -p  Make output periodic, i.e. match derivatives
  102.  *            at ends. First and last input values must
  103.  *            agree. Cannot be used with -k option.
  104.  *
  105.  *        -x  Next 1 (or 2) arguments are lower (and upper)
  106.  *            "x" limits. Normally these limits are
  107.  *            calculated from the data. Automatic abscissae
  108.  *            start at lower limit (default 0). If either
  109.  *            argument is outside of the range of
  110.  *            abscissae, it is ignored.
  111.  *
  112.  * Diagnostics:    When data is not strictly monotone in "x", SPLINE
  113.  *        reproduces the input without interpolating extra
  114.  *        points.
  115.  *
  116.  * Bugs:    A limit of 1000 input points is silently
  117.  *        enforced.
  118.  *
  119.  *        The -n option has not been implemented in
  120.  *        accordance with the "UNIX Programmer's Manual"
  121.  *        specification. This was done to avoid ambiguities
  122.  *        when the -n option follows the -x option with one
  123.  *        argument.
  124.  *
  125.  *        At certain negative values for the -k option (for
  126.  *        example, k equals -4.0), the curve becomes
  127.  *        discontinuous. The -k option value has thus been
  128.  *        arbitrarily constrained to be greater than or
  129.  *        equal to zero.
  130.  *
  131.  * Credits:    The above description is a reworded and expanded
  132.  *        version of that appearing in the "UNIX Programmer's
  133.  *        Manual", copyright 1979, 1983 Bell Laboratories.
  134.  */
  135.  
  136. /*** Definitions ***/
  137.  
  138. #define FALSE         0
  139. #define TRUE         1
  140. #define MAX_SIZE  1000    /* Input point array limit */
  141.  
  142. #define ILL_ARG         0    /* Error codes */
  143. #define ILL_CMB      1
  144. #define ILL_KVL         2
  145. #define ILL_NVL         3
  146. #define ILL_OPT      4
  147. #define ILL_XVL         5
  148. #define INS_INP      6
  149. #define MIS_KVL      7
  150. #define MIS_NVL      8
  151. #define MIS_XVL      9
  152. #define MIS_YVL        10
  153. #define NMT_ORD        11
  154.  
  155. #define SQUARE(a) a*a
  156. #define CUBE(a) a*a*a
  157.  
  158. /*** Typedefs ***/
  159.  
  160. typedef int BOOL;    /* Boolean flag */
  161.  
  162. /*** Include Files ***/
  163.  
  164. #include <stdio.h>
  165. #include <ctype.h>
  166. #include <math.h>
  167.  
  168. /*** Main Body of Program ***/
  169.  
  170. int main(argc,argv)
  171. int argc;
  172. char **argv;
  173. {
  174.   int n = 0,
  175.       i,
  176.       j,
  177.       n_val = 0,
  178.       atoi();
  179.   float    x[MAX_SIZE],
  180.     y[MAX_SIZE],
  181.     h[MAX_SIZE-1];
  182.   double a_val = 1.0,
  183.      k_val = 0.0,
  184.      x1_val = 0.0,
  185.      x2_val = 0.0,
  186.      x_intvl,
  187.      ix,
  188.      iy,
  189.      d2y[MAX_SIZE],
  190.      atof(),
  191.      spl_int();
  192.   char buffer[257],
  193.        *temp,
  194.        *gets();
  195.   BOOL aflag = FALSE,    /* Command-line option flags */
  196.        kflag = FALSE,
  197.        pflag = FALSE,
  198.        x1flag = FALSE,
  199.        x2flag = FALSE,
  200.        is_float();
  201.   void spl_coeff(),
  202.        pspl_coeff(),
  203.        error();
  204.  
  205.   /* Parse the command line for user-selected options */
  206.  
  207.   while(--argc)
  208.   {
  209.     temp = *++argv;
  210.     if(*temp != '-')    /* Check for legal option flag */
  211.       error(ILL_OPT,*argv);
  212.     else
  213.     {
  214.       temp++;
  215.       switch(toupper(*temp))
  216.       {
  217.     case 'A':    /* "-a" option */
  218.       aflag = TRUE;
  219.       if(argc > 1 && is_float(*(argv+1)))
  220.       {
  221.         argc--;
  222.         argv++;
  223.          if((a_val = atof(*argv)) <= 0.00)
  224.           error(ILL_ARG,*argv);
  225.       }
  226.       break;
  227.     case 'K':    /* "-k" option */
  228.       if(pflag == TRUE)
  229.         error(ILL_CMB,NULL);
  230.       kflag = TRUE;
  231.       if(argc > 1 && is_float(*(argv+1)))
  232.       {
  233.         argc--;
  234.         argv++;
  235.         k_val = atof(*argv);
  236.         if(k_val < 0.00)
  237.           error(ILL_KVL,*argv);
  238.         break;
  239.       }
  240.       else
  241.         error(MIS_KVL,NULL);
  242.     case 'N':    /* "-n" option */
  243.       if(argc > 1)
  244.       {
  245.         argc--;
  246.         argv++;
  247.         if((n_val = atoi(*argv)) < 1)
  248.           error(ILL_NVL,*argv);
  249.         else
  250.           break;
  251.       }
  252.       else
  253.         error(MIS_NVL,NULL);
  254.     case 'P':    /* "-p" option */
  255.       if(kflag == TRUE)
  256.         error(ILL_CMB,NULL);
  257.       pflag = TRUE;
  258.       break;
  259.     case 'X':    /* "-x" option */
  260.       x1flag = TRUE;
  261.       if(argc > 1 && is_float(*(argv+1)))
  262.       {
  263.         argc--;
  264.         argv++;
  265.         x1_val = atof(*argv);
  266.       }
  267.       else
  268.         error(MIS_XVL,NULL);
  269.       if(argc > 1 && is_float(*(argv+1)))
  270.       {
  271.         x2flag = TRUE;
  272.         argc--;
  273.         argv++;
  274.         x2_val = atof(*argv);
  275.         if(x2_val <= x1_val)
  276.           error(ILL_XVL,x2_val);
  277.       }
  278.       break;
  279.     default:    /* "-n" option */
  280.       error(ILL_OPT,*argv);
  281.       }
  282.     }
  283.   }
  284.   if(n_val == 0)    /* Set "n_val" if not given */
  285.     n_val = 100;
  286.  
  287.   /* Get the input data */
  288.  
  289.   while(1)        /* ... while there is more input data */
  290.   {
  291.     if(aflag == TRUE)    /* Automatic abscissae were called for */
  292.     {
  293.       if(n == 0)
  294.     x[0] = x1_val;
  295.       else
  296.     x[n] = x[n-1] + a_val;
  297.     }
  298.     else        /* Abscissae supplied with input data */
  299.     {
  300.       if(gets(buffer))
  301.     x[n] = atof(buffer);
  302.       else
  303.     break;
  304.     }
  305.     if(gets(buffer))    /* Read in the corresponding ordinate */
  306.       y[n] = atof(buffer);
  307.     else
  308.     {
  309.       if(aflag == TRUE)
  310.     break;
  311.       else
  312.     error(MIS_YVL,NULL);
  313.     }
  314.     if(++n == MAX_SIZE)    /* Maximum amount of input data? */
  315.       break;
  316.   }
  317.   if(n < 4)        /* Check for insufficient input data */
  318.     error(INS_INP,NULL);
  319.  
  320.   /* Check for non-monotonic abscissae. Output input data set
  321.    * without interpolation if true.
  322.    */
  323.  
  324.   for(i = 0; i < n-1; i++)
  325.     if((h[i] = (x[i+1] - x[i])) <= 0.0)
  326.     {
  327.       for(i = 0; i < n; i++)
  328.     printf("%g\n%g\n",x[i],y[i]);
  329.       exit();
  330.     }
  331.  
  332.   /* Calculate abscissa interval. Use "-x" option values if
  333.    * they were given unless they fall outside the range of
  334.    * given (or calculated) abscissae.
  335.    */ 
  336.  
  337.   if(x1flag == FALSE || x1_val < x[0])
  338.     x1_val = x[0];
  339.   if(x2flag == FALSE || x2_val > x[n-1])
  340.     x2_val = x[n-1];
  341.   x_intvl = (x2_val - x1_val)/n_val;
  342.  
  343.   /* Find the coefficients */
  344.  
  345.   if(pflag == FALSE)
  346.     spl_coeff(y,d2y,h,n,k_val);
  347.   else
  348.     pspl_coeff(y,d2y,h,n);
  349.  
  350.   /* Interpolate and output results */
  351.  
  352.   ix = x1_val;
  353.   i = 0;
  354.   for(j = 0; j <= n_val; j++)
  355.   {
  356.     while(ix >= x[i+1] && i < n - 2)
  357.       i++;
  358.     iy = spl_int(ix,x,y,d2y,h,i);
  359.     printf("%g\n%g\n",ix,iy);
  360.     ix += x_intvl;
  361.   }
  362. }
  363.  
  364. /*** Functions ***/
  365.  
  366. /* SPL_COEFF() - Calculate spline coefficients and return in
  367.  *         vector "d2y[]". Matrix to be solved has the
  368.  *         typical form:
  369.  *
  370.  * +-                         -++-       -+   +-      -+
  371.  * | c[1]  h[1]  0     ....    0       0      || y"[1]   | = | d[1]   |
  372.  * | h[1]  c[2]  h[2]  ....    0       0      || y"[2]   |   | d[2]   |
  373.  * | 0     h[2]  c[3]  ....    0       0      || y"[3]   |   | d[3]   |
  374.  * | ....  ....  ....  ....    ....    ....   || .....   |   | ....   |
  375.  * | 0     0     ....  h[n-4]  c[n-3]  h[n-3] || y"[n-3] |   | d[n-3] |
  376.  * | 0     0     ....  0       h[n-3]  c[n-2] || y"[n-2] |   | d[n-2] |
  377.  * +-                         -++-       -+   +-      -+
  378.  *
  379.  * where c[1] = (2 + k_val) * h[0] + 2 * h[1]
  380.  *       c[i] = 2 * (h[i-1] + h[i]), i = 2,...,n-3
  381.  *       c[n-2] = 2 * h[n-3] + (2 + k_val) * h[n-2]
  382.  *
  383.  *            +-                   -+
  384.  *                | y[i+1] - y[i]   y[i] - y[i-1] |
  385.  *       d[i] = 6 * | ------------- - ------------- | , i = 1,...,n-2
  386.  *            |     h[i]             h[i-1]        |
  387.  *            +-                   -+
  388.  */
  389.  
  390. void spl_coeff(y,d2y,h,n,k_val)
  391. float y[],
  392.       h[];
  393. double d2y[],
  394.        k_val;
  395. int n;
  396. {
  397.   double c[MAX_SIZE-1];
  398.   int i;
  399.  
  400.   /* Set up the (symmetric tridiagonal) matrix. Array "d2y[]"
  401.    * initially holds the constants vector, then is overlaid with
  402.    * the calculated variables vector.
  403.    */
  404.  
  405.   c[1] = (2.0 + k_val) * h[0] + 2.0 * h[1];
  406.   c[n-2] = 2.0 * h[n-3] + (2.0 + k_val) * h[n-2];
  407.   for(i = 2; i < n-2; i++)
  408.     c[i] = 2.0 * (h[i-1] + h[i]);
  409.  
  410.   for(i = 1; i < n-1; i++)
  411.     d2y[i] = 6.0 * ((y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1]);
  412.  
  413.   /* Reduce the matrix to upper triangular form */ 
  414.  
  415.   for(i = 1; i < n-2; i++)
  416.   {
  417.     c[i+1] -= h[i] * h[i]/c[i];
  418.     d2y[i+1] -= d2y[i] * h[i]/c[i];
  419.   }
  420.  
  421.   /* Solve using back substitution */
  422.  
  423.   d2y[n-2] /= c[n-2];
  424.   for(i = n-3; i > 0; i--)
  425.     d2y[i] = (d2y[i] - h[i] * d2y[i+1])/c[i];
  426.  
  427.   /* Solve for end conditions */
  428.  
  429.   d2y[0] = d2y[1] * k_val;
  430.   d2y[n-1] = d2y[n-2] * k_val;
  431. }
  432.  
  433. /* PSPL_COEFF() - Calculate periodic spline coefficients and
  434.  *          return in vector "d2y[]". Matrix to be solved
  435.  *          has the typical form:
  436.  *
  437.  * +-                         -++-       -+   +-      -+
  438.  * | c[1]  h[1]  0     ....    0       h[0]   || y"[1]   | = | d[1]   |
  439.  * | h[1]  c[2]  h[2]  ....    0       0      || y"[2]   |   | d[2]   |
  440.  * | 0     h[2]  c[3]  ....    0       0      || y"[3]   |   | d[3]   |
  441.  * | ....  ....  ....  ....    ....    ....   || .....   |   | ....   |
  442.  * | 0     0     ....  h[n-3]  c[n-2]  h[n-2] || y"[n-2] |   | d[n-2] |
  443.  * | h[0]  0     ....  0       h[n-2]  c[n-1] || y"[n-1] |   | d[n-1] |
  444.  * +-                         -++-       -+   +-      -+
  445.  *
  446.  * where c[i] = 2 * (h[i-1] + h[i]), i = 1,...,n-2
  447.  *       c[n-1] = 2 * (h[0] + h[n-2])
  448.  *
  449.  *            +-                   -+
  450.  *            | y[i+1] - y[i]   y[i] - y[i-1] |
  451.  *       d[i] = 6 * | ------------- - ------------- | , i = 1,...,n-2
  452.  *            |     h[i]             h[i-1]        |
  453.  *            +-                   -+
  454.  *
  455.  *              +-                 -+
  456.  *              | y[1] - y[0]   y[n-1] - y[n-2] |
  457.  *       d[n-1] = 6 * | ----------- - --------------- |
  458.  *              |    h[0]               h[n-2]     |
  459.  *              +-                 -+
  460.  *
  461.  */
  462.  
  463. void pspl_coeff(y,d2y,h,n)
  464. float y[],
  465.       h[];
  466. double d2y[];
  467. int n;
  468. {
  469.   double a,
  470.      b,
  471.      f,
  472.      fabs();
  473.   static double c[MAX_SIZE-1],
  474.         e[MAX_SIZE-1];
  475.   int i;
  476.  
  477.   /* Check for matching end ordinates */
  478.  
  479.   if(fabs(y[n-1] - y[0]) > 0.0)
  480.     error(NMT_ORD,NULL);
  481.  
  482.   /* Set up the matrix. Array "d2y[]" initially holds the constants
  483.    * vector, then is overlaid with the calculated variables vector.
  484.    */
  485.  
  486.   for(i = 1; i < n-1; i++)
  487.   {
  488.     c[i] = 2.0 * (h[i-1] + h[i]);
  489.     d2y[i] = 6.0 * ((y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1]);
  490.   }
  491.   c[n-1] = 2.0 * (h[0] + h[n-2]);
  492.   d2y[n-1] = 6.0 * ((y[1] - y[0])/h[0] - (y[n-1] - y[n-2])/h[n-2]);
  493.  
  494.   /* Initialize array e[] as nth column of matrix M[][]. */
  495.  
  496.   e[1] = h[0];
  497.   for(i = 2; i < n-2; i++)
  498.     e[i] = 0.0;
  499.   e[n-2] = h[n-2];
  500.   e[n-1] = c[n-1];
  501.  
  502.   /* Initialize variable f as matrix element M[n][1] */
  503.  
  504.   f = h[0];
  505.  
  506.   /* Reduce matrix to upper triangular form */
  507.  
  508.   for(i = 1; i < n-2; i++)
  509.   {
  510.     a = h[i]/c[i];
  511.     b = f/c[i];
  512.     c[i+1] -= h[i] * a;
  513.     d2y[i+1] -= d2y[i] * a;
  514.     e[i+1] -= e[i] * a;
  515.     d2y[n-1] -= d2y[i] * b;
  516.     e[n-1] -= e[i] * b;
  517.     f = -f * a;        /* Now matrix element M[n-1][i] */
  518.   }
  519.   f = f + h[n-2];    /* Now matrix element M[n-1][n-2] */
  520.   d2y[n-1] -= d2y[n-2] * f/c[n-2];
  521.   e[n-1] -= e[n-2] * f/c[n-2];
  522.  
  523. /* Solve using back substitution */
  524.  
  525.   d2y[n-1] /= e[n-1];
  526.   d2y[n-2] = (d2y[n-2] - e[n-2] * d2y[n-1])/c[n-2];
  527.   for(i = n-3; i > 0; i--)
  528.     d2y[i] = (d2y[i] - h[i] * d2y[i+1] - e[i] * d2y[n-1])/c[i];
  529.  
  530.   d2y[0] = d2y[n-1]; /* End condition */
  531. }
  532.  
  533. /* SPL_INT - Interpolate points using spline function */
  534.  
  535. double spl_int(ix,x,y,d2y,h,i)
  536. float x[],
  537.       y[],
  538.       h[];
  539. double ix,
  540.        d2y[];
  541. int i;
  542. {
  543.   double iy,
  544.      t1,
  545.      t2;
  546.  
  547.   t1 = (x[i+1] - ix)/h[i];
  548.   t2 = (ix - x[i])/h[i];
  549.   iy = y[i] * t1 + y[i+1] * t2 - SQUARE(h[i]) * (d2y[i] * (t1 -
  550.       CUBE(t1)) + d2y[i+1] * (t2 - CUBE(t2)))/6.0;
  551.   return iy;
  552. }
  553.  
  554. /* IS_FLOAT() - Check that character string is in correct floating
  555.  *        point format. Return TRUE if correct, FALSE
  556.  *        otherwise. The algorithm used is a deterministic
  557.  *        finite state machine. Using the regular
  558.  *        expression terminology of Unix's "lex", the
  559.  *        character string must be of the form:
  560.  *
  561.  *            -?d*.?d*(e|E(\+|-)?d+)?
  562.  *
  563.  *        where d = 0|1|2|3|4|5|6|7|8|9
  564.  */
  565.  
  566. BOOL is_float(str)
  567. char *str;
  568. {
  569.   int c,        /* Next FSM input character */
  570.       state = 0;    /* Current FSM state */
  571.  
  572.   while(c = *str++)
  573.   {
  574.     switch(state)
  575.     {
  576.       case 0:        /* FSM State 0 */
  577.     switch(c)
  578.     {
  579.       case '-':    
  580.         state = 1;
  581.         break;
  582.       case '.':
  583.         state = 3;
  584.         break;
  585.       default:
  586.         if(isdigit(c))
  587.           state = 2;
  588.         else
  589.           return FALSE;
  590.         break;
  591.     }
  592.     break;
  593.       case 1:        /* FSM State 1 */
  594.     switch(c)
  595.     {
  596.       case '.':
  597.         state = 2;
  598.         break;
  599.       default:
  600.         if(isdigit(c))
  601.           state = 2;
  602.         else
  603.           return FALSE;
  604.         break;
  605.     }
  606.     break;
  607.       case 2:        /* FSM State 2 */
  608.     switch(c)
  609.     {
  610.       case '.':
  611.         state = 4;
  612.         break;
  613.       case 'e':
  614.       case 'E':
  615.         state = 5;
  616.         break;
  617.       default:
  618.         if(isdigit(c))
  619.           state = 2;
  620.         else
  621.           return FALSE;
  622.         break;
  623.     }
  624.     break;
  625.       case 3:        /* FSM State 3 */
  626.     if(isdigit(c))
  627.       state = 4;
  628.     else
  629.       return FALSE;
  630.     break;
  631.       case 4:        /* FSM State 4 */
  632.     switch(c)
  633.     {
  634.       case 'e':
  635.       case 'E':
  636.         state = 5;
  637.         break;
  638.       default:
  639.         if(isdigit(c))
  640.           state = 4;
  641.         else
  642.           return FALSE;
  643.         break;
  644.     }
  645.     break;
  646.       case 5:        /* FSM State 5 */
  647.     switch(c)
  648.     {
  649.       case '+':
  650.       case '-':
  651.         state = 6;
  652.         break;
  653.       default:
  654.         if(isdigit(c))
  655.           state = 7;
  656.         else
  657.           return FALSE;
  658.         break;
  659.     }
  660.     break;
  661.       case 6:        /* FSM State 6 */
  662.     if(isdigit(c))
  663.       state = 7;
  664.     else
  665.       return FALSE;
  666.     break;
  667.       case 7:        /* FSM State 7 */
  668.     if(isdigit(c))
  669.       state = 7;
  670.     else
  671.       return FALSE;
  672.     break;
  673.     }
  674.   }
  675.   return TRUE;
  676. }
  677.  
  678. /* ERROR() - Error reporting. Returns an exit status of 2 to the
  679.  *         parent process.
  680.  */
  681.  
  682. void error(n,str)
  683. int n;
  684. char *str;
  685. {
  686.   fprintf(stderr,"\007\n*** ERROR - ");
  687.   switch(n)
  688.   {
  689.     case ILL_ARG:
  690.       fprintf(stderr,"Argument must be greater than zero: %s",
  691.       str);
  692.       break;
  693.     case ILL_CMB:
  694.       fprintf(stderr,"Cannot use -k option with -p option");
  695.       break;
  696.     case ILL_KVL:
  697.       fprintf(stderr,"Illegal value for -k option: %s",str);
  698.       break;
  699.     case ILL_NVL:
  700.       fprintf(stderr,"Illegal value for -n option: %s",str);
  701.       break;
  702.     case ILL_OPT:
  703.       fprintf(stderr,"Illegal command line option: %s",str);
  704.       break; 
  705.     case ILL_XVL:
  706.       fprintf(stderr,"Illegal value for -x option: %s",str);
  707.       break;
  708.     case INS_INP:
  709.       fprintf(stderr,"Insufficient input data");
  710.       break; 
  711.     case MIS_KVL:
  712.       fprintf(stderr,"Missing value for -k option");
  713.       break;
  714.     case MIS_NVL:
  715.       fprintf(stderr,"Missing value for -n option");
  716.       break;
  717.     case MIS_XVL:
  718.       fprintf(stderr,"Missing value for -x option");
  719.       break;
  720.     case MIS_YVL:
  721.       fprintf(stderr,"Missing ordinate value");
  722.       break;
  723.     case NMT_ORD:
  724.       fprintf(stderr,"End ordinates do not match");
  725.       break;
  726.     default:
  727.       break;
  728.   }
  729.   fprintf(stderr," ***\n\nUsage: spline [-aknpx]\n");
  730.   exit(2);
  731. }
  732.  
  733. /*** End of SPLINE.C ***/
  734.  case ILL_KVL:
  735.       fprintf(stderr,"Illegal value