/* This program, spline, interpolates scalar or vector-valued input data
   using splines with tension, including piecewise cubic (zero-tension)
   splines.  When acting as a real-time filter, it uses cubic Bessel
   interpolation instead.  Written by Robert S. Maier
   <rsm@math.arizona.edu>, based on earlier work by Rich Murphey.
   Copyright (C) 1989-1999 Free Software Foundation, Inc.

   References:

   D. Kincaid and [E.] W. Cheney, Numerical Analysis, Brooks/Cole,
   2nd. ed., 1996, Section 6.4.

   C. de Boor, A Practical Guide to Splines, Springer-Verlag, 1978, 
   Chapter 4.

   A. K. Cline, "Scalar and Planar-Valued Curve Fitting Using Splines under
   Tension", Communications of the ACM 17 (1974), 218-223.

   The tension in a spline is set with the -T (i.e., --tension) option.  By
   definition, a one-dimensional spline with tension satisfies the
   differential equation y''''=sgn(tension)*(tension**2)y''.  The default
   value for the tension is zero.  If tension=0 then a spline with tension
   reduces to a conventional piecewise cubic spline.  In the limits
   tension->+infinity and tension->-infinity, a spline with tension reduces
   to a piecewise linear (`broken line') interpolation.

   To oversimplify a bit, 1.0/tension is the maximum abscissa range over
   which the spline likes to curve, at least when tension>0.  So increasing
   the tension far above zero tends to make the spline contain short curved
   sections, separated by sections that are almost straight.  The curved
   sections will be centered on the user-specified data points.  The
   behavior of the spline when tension<0 is altogether different: it will
   tend to oscillate, though as tension->-infinity the oscillations are
   damped out.

   Tension is a `dimensionful' quantity.  If tension=0 (the cubic spline
   case), then the computation of the spline is scale-invariant.  But if
   the tension is nonzero, then when the abscissa values are multiplied by
   some common positive factor, the tension should be divided by the same
   factor to obtain a scaled version of the original spline.

   The algorithms of Kincaid and Cheney have been extended to include
   support for periodicity.  To obtain a periodic spline, with or without
   tension, the user uses the -p (i.e., --periodic) option and supplies
   input data satisfying y[n]=y[0].  Also, in the non-periodic case the
   algorithms have been extended to include support for a parameter k,
   which appears in the two boundary conditions y''[0]=ky''[1] and
   y''[n]=ky''[n-1].  The default value of k is 1.0.  The parameter k,
   which is specified with the -k (i.e. --boundary-condition) option, is
   ignored for periodic splines (using the -k option with the -p option
   will elicit a warning).

   If the -f option is specified, then an altogether different (real-time)
   algorithm for generating interpolating points will be used, so that this
   program can be used as a real-time filter.  If -f is specified then the
   -t option, otherwise optional, must also be used.  (I.e., the minimum
   and maximum abscissa values for the interpolating points must be
   specified, and optionally the spacing between them as well.  If the
   spacing is not specified on the command line, then the interval
   [tmin,tmax] will be subdivided into a default number of intervals [100],
   unless the default number of intervals is overridden with the -n option.

   The real-time algorithm that is used when the -f option is specified is
   cubic Bessel interpolation.  (The -T, -p, and -k options are ignored
   when -f is specified; using them will elicit a warning.)  Interpolation
   in this case is piecewise cubic, and the slopes at either end of each
   sub-interval are found by fitting a parabola through each successive
   triple of points.  That is, the slope at t=t_n is found by fitting a
   parabola through the points at t_(n-1), t_n, and t_(n+1).  This
   interpolation scheme yields a spline that is only once, rather than
   twice, continuously differentiable.  However, it has the feature that
   all computations are local rather than global, so it is suitable for
   real-time work.

   Since the above was written, the -d option has been added, to permit the
   splining of multidimensional data.  All components of a d-dimensional
   data set (a d-dimensional vector y is specified at each t) are splined
   in the same way, as if they were one-dimensional functions of t.  All
   options that apply to 1-dimensional datasets, such as -T, -p, -k, -f,
   etc., apply to d-dimensional ones also. */

#include "sys-defines.h"
#include "getopt.h"

#define	ARG_NONE	0
#define	ARG_REQUIRED	1
#define	ARG_OPTIONAL	2

/* states for cubic Bessel DFA; occupancy of data point queue */
enum { STATE_ZERO, STATE_ONE, STATE_TWO, STATE_THREE };

/* types of auto-abscissa */
enum { AUTO_NONE, AUTO_INCREMENT, AUTO_BY_DISTANCE };

#define FUZZ 0.0000001		/* potential roundoff error */

/* Minimum value for magnitude of x, for such functions as x-sinh(x),
   x-tanh(x), x-sin(x), and x-tan(x) to have acceptable accuracy.  If the
   magnitude of x is smaller than this value, these functions of x will be
   computed via power series to accuracy O(x**6). */
#define TRIG_ARG_MIN 0.001

/* Maximum value for magnitude of x, beyond which we approximate
   x/sinh(x) and x/tanh(x) by |x|exp(-|x|). */
#define TRIG_ARG_MAX 50.0

struct option long_options[] =
{
  {"no-of-intervals",	ARG_REQUIRED,	NULL, 'n'},
  {"periodic",		ARG_NONE,	NULL, 'p'},
  {"y-dimension",	ARG_REQUIRED,	NULL, 'd'},
  {"t-limits",		ARG_REQUIRED,	NULL, 't'}, /* 1 or 2 or 3 */
  {"t-limits",		ARG_REQUIRED,	NULL, 'x'}, /* obsolescent; hidden */
  {"tension",		ARG_REQUIRED, 	NULL, 'T'},
  {"boundary-condition",ARG_REQUIRED,	NULL, 'k'},
  {"auto-abscissa",	ARG_OPTIONAL,	NULL, 'a'}, /* 0 or 1 or 2 */
  {"auto-dist-abscissa",ARG_NONE,	NULL, 'A'},
  {"filter",		ARG_NONE,	NULL, 'f'},
  {"precision",		ARG_REQUIRED,	NULL, 'P'},
  {"suppress-abscissa",	ARG_NONE,	NULL, 's'},
  /* ascii or double */
  {"input-type",	ARG_REQUIRED,	NULL, 'I'},
  {"output-type",	ARG_REQUIRED,	NULL, 'O'},
  /* Long options with no equivalent short option alias */
  {"version",		ARG_NONE,	NULL, 'V' << 8},
  {"help",		ARG_NONE,	NULL, 'h' << 8},
  {NULL, 		0, 		0,     0}
};

/* null-terminated list of options that we don't show to the user */
int hidden_options[] = { (int)'x', 0 };

/* type of data in input and output streams */
typedef enum
{
  T_ASCII, T_SINGLE, T_DOUBLE, T_INTEGER
}
data_type;

data_type input_type = T_ASCII;
data_type output_type = T_ASCII;

const char *progname = "spline"; /* name of this program */

const char *usage_appendage = " [FILE]...\n\
With no FILE, or when FILE is -, read standard input.\n";

/* forward references */
bool do_bessel ____P ((FILE *input, int ydimension, int auto_abscissa, double auto_t, double auto_delta, double first_t, double last_t, double spacing_t, int precision, bool suppress_abscissa));
bool is_monotonic ____P ((int n, double *t));
bool read_data ____P ((FILE *input, int *len, int *used, int auto_abscissa, double auto_t, double auto_delta, double **t, int ydimension, double **y, double **z));
bool read_float ____P((FILE *input, double *dptr));
bool skip_whitespace ____P ((FILE *stream));
bool write_point ____P((double t, double *y, int ydimension, int precision, bool suppress_abscissa));
double interpolate ____P ((int n, double *t, double *y, double *z, double x, double tension, bool periodic));
double quotient_sin_func ____P((double x, double y));
double quotient_sinh_func ____P((double x, double y));
double sin_func ____P((double x));
double sinh_func ____P((double x));
double tan_func ____P((double x));
double tanh_func ____P((double x));
int read_point ____P ((FILE *input, double *t, double *y, int ydimension, bool *first_point, int auto_abscissa, double *auto_t, double auto_delta, double *stored));
void do_bessel_range ____P ((double abscissa0, double abscissa1, double *value0, double *value1, double *slope0, double *slope1, double first_t, double last_t, double spacing_t, int ydimension, int precision, bool endit, bool suppress_abscissa));
void do_spline ____P ((int used, int len, double **t, int ydimension, double **y, double **z, double tension, bool periodic, bool spec_boundary_condition, double boundary_condition, int precision, double first_t, double last_t, double spacing_t, int no_of_intervals, bool spec_first_t, bool spec_last_t, bool spec_spacing_t, bool spec_no_of_intervals, bool suppress_abscissa));
void fit ____P ((int n, double *t, double *y, double *z, double k, double tension, bool periodic));
void maybe_emit_oob_warning ____P ((void));
void non_monotonic_error ____P((void));
void output_dataset_separator ____P ((void));
void set_format_type ____P ((char *s, data_type *typep));
/* from libcommon */
extern void display_usage ____P((const char *progname, const int *omit_vals, const char *appendage, bool fonts));
extern void display_version ____P((const char *progname)); 
extern voidptr_t xcalloc ____P ((size_t nmemb, size_t size));
extern voidptr_t xmalloc ____P ((size_t size));
extern voidptr_t xrealloc ____P ((voidptr_t p, size_t length));
extern char *xstrdup ____P ((const char *s));


int
#ifdef _HAVE_PROTOS
main (int argc, char *argv[])
#else
main (argc, argv)
     int argc;
     char *argv[];
#endif
{
  int option;
  int opt_index;
  int errcnt = 0;		/* errors encountered */
  bool show_version = false;	/* remember to show version message */
  bool show_usage = false;	/* remember to output usage message */
  bool dataset_follows;

  /* parameters controlled by command line options: */
  bool filter = false;		/* act as a filter (cubic Bessel)? */
  bool periodic = false;	/* spline should be periodic? */
  bool spec_boundary_condition = false;	/* user-specified boundary cond'n? */
  bool spec_first_t = false, spec_last_t = false, spec_spacing_t = false;
  bool spec_no_of_intervals = false; /* user-specified number of intervals? */
  bool suppress_abscissa = false; /* for each point, print ordinate only? */
  double boundary_condition = 1.0; /* force  y''_1 = k * y''_0, etc. */
  double delta_t = 1.0;		/* increment of auto abscissa */
  double first_t = 0.0, last_t = 0.0, spacing_t = 0.0; /* values of limits */
  double tension = 0.0;		/* `tension' parameter */
  double t_start = 0.0;		/* start of auto abscissa */
  int auto_abscissa = AUTO_NONE; /* automatic generation of abscissa? */
  int no_of_intervals = 100;	/* no. of intervals to divide abs. range */
  int precision = 6;		/* default no. of significant digits printed */
  int ydimension = 1;		/* dimension of each point's ordinate */

  /* used in argument parsing */
  double local_first_t, local_last_t, local_spacing_t;
  double local_t_start, local_delta_t;
  int local_precision;

  for ( ; ; )
    {
      option = getopt_long (argc, argv, "fpsAd:I:O:P:k:n:t:x:T:a::", long_options, &opt_index);
      if (option == 0)
	option = long_options[opt_index].val;
      
      switch (option)
	{
	  /* ----------- options with no argument --------------*/

	case 'p':		/* construct periodic, i.e., closed spline */
	  periodic = true;
	  break;
	case 'f':		/* act as filter */
	  filter = true;
	  break;
	case 's':		/* don't output t values */
	  suppress_abscissa = true;
	  break;
	case 'A':		/* delta t = inter-y distance */
	  auto_abscissa = AUTO_BY_DISTANCE;
	  t_start = 0.0;
	  break;
	case 'V' << 8:		/* Version */
	  show_version = true;
	  break;
	case 'h' << 8:		/* Help */
	  show_usage = true;
	  break;

	  /*--------------options with a single argument--------*/

	case 'I':
	  set_format_type (optarg, &input_type);
	  break;
	case 'O':
	  set_format_type (optarg, &output_type);
	  break;
	case 'd':		/* dimensionality of ordinate variable */
	  if (sscanf (optarg, "%d", &ydimension) <= 0 || ydimension < 1)
	    {
	      fprintf (stderr, 
		       "%s: error: bad ordinate dimension `%s' (must be positive integer)\n", 
		       progname, optarg);
	      errcnt++;
	    }
	  break;
	case 'k':
	  if (sscanf (optarg, "%lf", &boundary_condition) <= 0)
	    {
	      fprintf (stderr, 
		       "%s: error: bad boundary condition argument `%s'\n",
		       progname, optarg);
	      errcnt++;
	    }
	  else
	    spec_boundary_condition = true;
	  break;
	case 'T':
	  if (sscanf (optarg, "%lf", &tension) <= 0)
	    {
	      fprintf (stderr, 
		       "%s: error: bad tension argument `%s'\n",
		       progname, optarg);
	      errcnt++;
	    }
	  break;
	case 'n':		/* number of intervals */
	  if (sscanf (optarg, "%d", &no_of_intervals) <= 0)
	    {
	      fprintf (stderr, 
		       "%s: error: bad requested number of intervals `%s'\n", 
		       progname, optarg);
	      errcnt++;
	    }
	  else
	    spec_no_of_intervals = true;
	  break;
	case 'P':		/* precision */
	  if (sscanf (optarg, "%d", &local_precision) <= 0)
	    {
	      fprintf (stderr, "%s: error: bad requested precision `%s' (must be positive integer)\n", 
		       progname, optarg);
	      errcnt++;
	    }
	  else
	    {
	      if (local_precision <= 0)
		fprintf (stderr, 
			 "%s: ignoring bad precision value `%s' (must be positive integer)\n",
			 progname, optarg);
	      else
		precision = local_precision;
	    }
	  break;

	  /*------------options with 0 or more args ----------*/

	case 'a':		/* Auto-abscissa, ARG OPTIONAL [0,1,2] */
	  auto_abscissa = AUTO_INCREMENT;
	  if (optind >= argc)
	    break;
	  if (sscanf (argv[optind], "%lf", &local_delta_t) <= 0)
	    break;
	  delta_t = local_delta_t;
	  optind++;	/* tell getopt we recognized delta_t */
	  if (optind >= argc)
	    break;
	  if (sscanf (argv [optind], "%lf", &local_t_start) <= 0)
	    break;
	  t_start = local_t_start;
	  optind++;	/* tell getopt we recognized t_start */
	  break;

	  /*--------------options with 1 or more arguments------*/

	case 't':		/* t axis limits, ARG REQUIRED [1,2,3] */
	case 'x':		/* obsolescent variant */
	  if (sscanf (optarg, "%lf", &local_first_t) <= 0)
	    break;
	  first_t = local_first_t;
	  spec_first_t = true;
	  if (optind >= argc)
	    break;
	  if (sscanf (argv [optind], "%lf", &local_last_t) <= 0)
	    break;
	  last_t = local_last_t;
	  spec_last_t = true;
	  optind++;	/* tell getopt we recognized last_t */
	  if (optind >= argc)
	    break;
	  if (sscanf (argv [optind], "%lf", &local_spacing_t) <= 0)
	    break;
	  spacing_t = local_spacing_t;
	  spec_spacing_t = true;
	  optind++;	/* tell getopt we recognized spacing_t */
	  break;

	  /*---------------- End of options ----------------*/

	default:		/* Default, unknown option */
	  errcnt++;
	  break;
	}			/* endswitch */

      if ((option == EOF))
	{
	  errcnt--;
	  break;		/* break out of option processing */
	}
    }
				/* endwhile */
  if (errcnt > 0)
    {
      fprintf (stderr, "Try `%s --help' for more information\n", progname);
      return EXIT_FAILURE;
    }
  if (show_version)
    {
      display_version (progname);
      return EXIT_SUCCESS;
    }
  if (show_usage)
    {
      display_usage (progname, hidden_options, usage_appendage, false);
      return EXIT_SUCCESS;
    }

  /* Some sanity checks on user-supplied options. */

  if (no_of_intervals < 1)
    {
      fprintf (stderr, 
	       "%s: error: cannot subdivide abscissa range into %d intervals\n", 
	       progname, no_of_intervals);
      return EXIT_FAILURE;
    }

  if (periodic)
    {
      if (spec_boundary_condition)
	fprintf (stderr, 
		 "%s: for periodic splines, setting of boundary condition not supported\n", 
		 progname);
      boundary_condition = 0.0;
    }

  if (filter)
    /* acting as a filter, so use cubic Bessel interpolation */
    {
      if (!spec_first_t || !spec_last_t)
	{
	  fprintf (stderr,
		   "%s: error: acting as filter, must specify abscissa range with -t option\n",
		   progname);
	  return EXIT_FAILURE;
	}

      if (!spec_spacing_t) 
	spacing_t = (last_t - first_t) / no_of_intervals;
      else			/* user specified spacing */
	{
	  if (spec_no_of_intervals)
	    fprintf (stderr, "%s: ignoring specified number of intervals\n",
		     progname);
	  if ((last_t - first_t) * spacing_t < 0.0)
	    {
	      fprintf (stderr, "%s: specified spacing is of wrong sign, corrected\n",
		       progname);
	      spacing_t = -spacing_t;
	    }

	  /* N.B. if spacing specified, should optionally contract first_t and
	     last_t to make them integer multiples of spacing; cf. graph */
	}
      
      if (spec_boundary_condition)
	fprintf (stderr, 
		 "%s: acting as filter, so setting of boundary condition not supported\n",
		 progname);
      if (tension != 0.0)
	fprintf (stderr, 
		 "%s: acting as filter, so nonzero tension not supported\n",
		 progname);
      if (periodic)
	fprintf (stderr, 
		 "%s: acting as filter, so periodicity not supported\n",
		 progname);

      if (optind < argc)
	{
	  /* call do_bessel() on each file specified on the command line,
	     generating a spline from each dataset in the file */
	  for (; optind < argc; optind++)
	    {
	      FILE *data_file;
	      
	      /* open file, treating "-" as stdin */
	      if (strcmp (argv[optind], "-") == 0)
		data_file = stdin;
	      else
		{
		  data_file = fopen (argv[optind], "r");
		  if (data_file == NULL)
		    {
		      fprintf (stderr, "%s: %s: %s\n", progname, argv[optind], strerror(errno));
		      return EXIT_FAILURE;
		    }
		}		

	      /* loop through datasets in file (may be more than one) */
	      do
		{
		  dataset_follows = do_bessel (data_file, ydimension,
					       auto_abscissa, t_start, delta_t,
					       first_t, last_t, spacing_t, 
					       precision, suppress_abscissa);

		  /* output a separator between successive datasets */
		  if (dataset_follows || (optind + 1 != argc))
		    output_dataset_separator();
		  
		} while (dataset_follows);

	      /* close file */
	      if (data_file != stdin) /* don't close stdin */
		{
		  if (fclose (data_file) < 0)
		    {
		      fprintf (stderr, 
			       "%s: error: couldn't close input file `%s'\n",
			       progname, argv[optind]);
		      return EXIT_FAILURE;
		    }
		}
	    }
	}
      else			/* no files spec'd, read stdin instead */
	/* loop through datasets read from stdin (may be more than one) */
	do
	  {
	    dataset_follows = do_bessel (stdin, ydimension,
					 auto_abscissa, t_start, delta_t,
					 first_t, last_t, spacing_t, 
					 precision, suppress_abscissa);
	    
	    /* output a separator between successive datasets */
	    if (dataset_follows)
	      output_dataset_separator();
	  }
	while (dataset_follows);	/* keep going if no EOF yet */
    }

  else
    /* not acting as filter, so use spline interpolation (w/ tension) */
    {
      double *t, **y, **z;	/* ordinate, abscissa, 2nd derivative arrays */
      int i, len, used;

      if (optind < argc)	/* files spec'd on command line */
	{

	  /* call do_spline() on each file specified on the command line,
	     generating a spline from each dataset contained in the file */
	  for (; optind < argc; optind++)
	    {
	      FILE *data_file;
	      
	      /* open file, treat "-" as meaning stdin */
	      if (strcmp (argv[optind], "-") == 0)
		data_file = stdin;
	      else
		{
		  data_file = fopen (argv[optind], "r");
		  if (data_file == NULL)
		    {
		      fprintf (stderr, "%s: error: couldn't open file `%s'\n",
			       progname, argv[optind]);
		      return EXIT_FAILURE;
		    }
		}
	      
	      /* loop through datasets in file (may be more than one) */
	      do
		{
		  len = 16;	/* initial value of storage length */
		  used = -1;	/* initial value of array size, minus 1 */
	      
		  t = (double *)xmalloc (sizeof(double) * len);
		  y = (double **)xmalloc (sizeof(double *) * ydimension);
		  z = (double **)xmalloc (sizeof(double *) * ydimension);
		  for (i = 0; i < ydimension; i++)
		    {
		      y[i] = (double *)xmalloc (sizeof(double) * len);
		      z[i] = (double *)xmalloc (sizeof(double) * len);
		    }
		  
		  dataset_follows = read_data (data_file, &len, &used, 
					       auto_abscissa, t_start, delta_t,
					       &t, ydimension, y, z);
		  /* read_data() may reallocate t,y[*],z[*], and update
		     len, used; on exit, used + 1 is number of data points */
		  
		  /* spline the dataset and output interpolating points */
		  do_spline (used, len, 
			     &t, ydimension, y, z, tension, periodic,
			     spec_boundary_condition, boundary_condition, 
			     precision,
			     first_t, last_t, spacing_t, no_of_intervals,
			     spec_first_t, spec_last_t, spec_spacing_t, 
			     spec_no_of_intervals, suppress_abscissa);

		  /* output a separator between successive datasets */
		  if (dataset_follows || (optind + 1 != argc))
		    output_dataset_separator();
		  
		  free (z);
		  free (y);
		  free (t);
		}
	      while (dataset_follows);	/* keep going if no EOF yet */
	      
	      /* close file */
	      if (data_file != stdin) /* don't close stdin */
		{
		  if (fclose (data_file) < 0)
		    {
		      fprintf (stderr, 
			       "%s: error: couldn't close input file `%s'\n",
			       progname, argv[optind]);
		      return EXIT_FAILURE;
		    }
		}
	    }
	}
      else			/* no files spec'd, read stdin instead */
	/* loop through datasets read from stdin (may be more than one) */
	do
	  {
	    len = 16;		/* initial value for array size */
	    used = -1;	/* initial number of stored points, minus 1 */
	    
	    t = (double *)xmalloc (sizeof(double) * len);
	    y = (double **)xmalloc (sizeof(double *) * ydimension);
	    z = (double **)xmalloc (sizeof(double *) * ydimension);
	    for (i = 0; i < ydimension; i++)
	      {
		y[i] = (double *)xmalloc (sizeof(double) * len);
		z[i] = (double *)xmalloc (sizeof(double) * len);
	      }
	    
	    dataset_follows = read_data (stdin, &len, &used, 
				     auto_abscissa, t_start, delta_t, 
				     &t, ydimension, y, z);
	    /* read_data() may reallocate t,y[*],z[*], and update len,
	       used; on exit, used + 1 is number of data points */
	    
	    /* spline the dataset and output interpolating points */
	    do_spline (used, len, 
		       &t, ydimension, y, z, tension, periodic,
		       spec_boundary_condition, boundary_condition, precision,
		       first_t, last_t, spacing_t, no_of_intervals,
		       spec_first_t, spec_last_t, spec_spacing_t, 
		       spec_no_of_intervals, suppress_abscissa);
	    
	    /* output a separator between successive datasets */
	    if (dataset_follows)
	      output_dataset_separator();
	    
	    for (i = 0; i < ydimension; i++)
	      {
		free (z[i]);
		free (y[i]);
	      }
	    free (z);
	    free (y);
	    free (t);
	  }
	while (dataset_follows);	/* keep going if no EOF yet */
      
    }

  return EXIT_SUCCESS;
}


void
#ifdef _HAVE_PROTOS
set_format_type (char *s, data_type *typep)
#else
set_format_type (s, typep)
     char *s;
     data_type *typep;
#endif
{
  switch (s[0])
    {
    case 'a':
    case 'A':
      *typep = T_ASCII;
      break;
    case 'f':
    case 'F':
      *typep = T_SINGLE;
      break;
    case 'd':
    case 'D':
      *typep = T_DOUBLE;
      break;
    case 'i':
    case 'I':
      *typep = T_INTEGER;
      break;
    default:
      {
	fprintf (stderr, "%s: error: invalid data format type `%s'\n",
		 progname, s);
	exit (EXIT_FAILURE);
      }
      break;
    }
}


/* fit() computes the array z[] of second derivatives at the knots, i.e.,
   internal data points.  The abscissa array t[] and the ordinate array y[]
   are specified.  On entry, have n+1 >= 2 points in the t, y, z arrays,
   numbered 0..n.  The knots are numbered 1..n-1 as in Kincaid and Cheney.
   In the periodic case, the final knot, i.e., (t[n-1],y[n-1]), has the
   property that y[n-1]=y[0]; moreover, y[n]=y[1].  The number of points
   supplied by the user was n+1 in the non-periodic case, and n in the
   periodic case.  When this function is called, n>=1 in the non-periodic
   case, and n>=2 in the periodic case. */

/* Algorithm: the n-1 by n-1 tridiagonal matrix equation for the vector of
   2nd derivatives at the knots is reduced to upper diagonal form.  At that
   point the diagonal entries (pivots) of the upper diagonal matrix are in
   the vector u[], and the vector on the right-hand side is v[].  That is,
   the equation is of the form Ay'' = v, where a_(ii) = u[i], and a_(i,i+1)
   = alpha[i].  Here i=1..n-1 indexes the set of knots.  The matrix
   equation is solved by back-substitution for y''[], i.e., for z[]. */

void
#ifdef _HAVE_PROTOS
fit (int n, double *t, double *y, double *z, double k, double tension,
     bool periodic)
#else
fit (n, t, y, z, k, tension, periodic)
     int n;
     double *t, *y, *z;
     double k;			/* y''_1 = k y''_0, etc. */
     double tension;
     bool periodic;
#endif
{
  double *h, *b, *u, *v, *alpha, *beta;
  double *uu = NULL, *vv = NULL, *s = NULL;
  int i;

  if (n == 1)			/* exactly 2 points, use straight line */
    {
      z[0] = z[1] = 0.0;
      return;
    }

  h = (double *)xmalloc (sizeof(double) * n);
  b = (double *)xmalloc (sizeof(double) * n);
  u = (double *)xmalloc (sizeof(double) * n);
  v = (double *)xmalloc (sizeof(double) * n);
  alpha = (double *)xmalloc (sizeof(double) * n);
  beta = (double *)xmalloc (sizeof(double) * n);

  if (periodic)
    {
      s = (double *)xmalloc (sizeof(double) * n); 
      uu = (double *)xmalloc (sizeof(double) * n); 
      vv = (double *)xmalloc (sizeof(double) * n); 
    }

  for (i = 0; i <= n - 1 ; ++i)
    {
      h[i] = t[i + 1] - t[i];
      b[i] = 6.0 * (y[i + 1] - y[i]) / h[i]; /* for computing RHS */
    }

  if (tension < 0.0)		/* must rule out sin(tension * h[i]) = 0 */
    {
      for (i = 0; i <= n - 1 ; ++i)
	if (sin (tension * h[i]) == 0.0)
	  {
	    fprintf (stderr, "%s: error: specified negative tension value is singular\n", progname);
	    exit (EXIT_FAILURE);
	  }
    }
  if (tension == 0.0)
    {
      for (i = 0; i <= n - 1 ; ++i)
	{
	  alpha[i] = h[i];	/* off-diagonal = alpha[i] to right */
	  beta[i] = 2.0 * h[i];	/* diagonal = beta[i-1] + beta[i] */
	}
    }
  else
    if (tension > 0.0)
      /* `positive' (really real) tension, use hyperbolic trig funcs */
      {
	for (i = 0; i <= n - 1 ; ++i)
	  {
	    double x = tension * h[i];
	    double xabs = (x < 0.0 ? -x : x);

	    if (xabs < TRIG_ARG_MIN)
	      /* hand-compute (6/x^2)(1-x/sinh(x)) and (3/x^2)(x/tanh(x)-1)
                 to improve accuracy; here `x' is tension * h[i] */
	      {
		alpha[i] = h[i] * sinh_func(x);
		beta[i] = 2.0 * h[i] * tanh_func(x);
	      }
	    else if (xabs > TRIG_ARG_MAX)
	      /* in (6/x^2)(1-x/sinh(x)) and (3/x^2)(x/tanh(x)-1),
		 approximate x/sinh(x) and x/tanh(x) by 2|x|exp(-|x|)
		 and |x|, respectively */
	      {
		int sign = (x < 0.0 ? -1 : 1);

		alpha[i] = ((6.0 / (tension * tension))
			   * ((1.0 / h[i]) - tension * 2 * sign * exp(-xabs)));
		beta[i] = ((6.0 / (tension * tension))
			   * (tension - (1.0 / h[i])));
	      }
	    else
	      {
		alpha[i] = ((6.0 / (tension * tension))
			    * ((1.0 / h[i]) - tension / sinh(x)));
		beta[i] = ((6.0 / (tension * tension))
			   * (tension / tanh(x) - (1.0 / h[i])));
	      }
	  }
      }
    else				/* tension < 0 */
      /* `negative' (really imaginary) tension,  use circular trig funcs */
      {
	for (i = 0; i <= n - 1 ; ++i)
	  {
	    double x = tension * h[i];
	    double xabs = (x < 0.0 ? -x : x);

	    if (xabs < TRIG_ARG_MIN)
	      /* hand-compute (6/x^2)(1-x/sin(x)) and (3/x^2)(x/tan(x)-1)
                 to improve accuracy; here `x' is tension * h[i] */
	      {
		alpha[i] = h[i] * sin_func(x);
		beta[i] = 2.0 * h[i] * tan_func(x);
	      }
	    else
	      {
		alpha[i] = ((6.0 / (tension * tension))
		           * ((1.0 / h[i]) - tension / sin(x)));
		beta[i] = ((6.0 / (tension * tension))
			   * (tension / tan(x) - (1.0 / h[i])));
	      }
	  }
      }
  
  if (!periodic && n == 2)
      u[1] = beta[0] + beta[1] + 2 * k * alpha[0];
  else
    u[1] = beta[0] + beta[1] + k * alpha[0];

  v[1] = b[1] - b[0];
  
  if (u[1] == 0.0)
    {
      fprintf (stderr, 
	       "%s: error: as posed, problem of computing spline is singular\n",
	       progname);
      exit (EXIT_FAILURE);
    }

  if (periodic)
    {
      s[1] = alpha[0];
      uu[1] = 0.0;
      vv[1] = 0.0;
    }

  for (i = 2; i <= n - 1 ; ++i)
    {
      u[i] = (beta[i] + beta[i - 1]
	      - alpha[i - 1] * alpha[i - 1] / u[i - 1]
	      + (i == n - 1 ? k * alpha[n - 1] : 0.0));

      if (u[i] == 0.0)
	{
	  fprintf (stderr, 
		   "%s: error: as posed, problem of computing spline is singular\n",
		   progname);
	  exit (EXIT_FAILURE);
	}


      v[i] = b[i] - b[i - 1] - alpha[i - 1] * v[i - 1] / u[i - 1];

      if (periodic)
	{
	  s[i] = - s[i-1] * alpha[i-1] / u[i-1];
	  uu[i] = uu[i-1] - s[i-1] * s[i-1] / u[i-1];
	  vv[i] = vv[i-1] - v[i-1] * s[i-1] / u[i-1];
	}
    }
      
  if (!periodic)
    {
      /* fill in 2nd derivative array */
      z[n] = 0.0;
      for (i = n - 1; i >= 1; --i)
	z[i] = (v[i] - alpha[i] * z[i + 1]) / u[i];
      z[0] = 0.0;
      
      /* modify to include boundary condition */
      z[0] = k * z[1];
      z[n] = k * z[n - 1];
    }
  else		/* periodic */
    {
      z[n-1] = (v[n-1] + vv[n-1]) / (u[n-1] + uu[n-1] + 2 * s[n-1]);
      for (i = n - 2; i >= 1; --i)
	z[i] = ((v[i] - alpha[i] * z[i + 1]) - s[i] * z[n-1]) / u[i];

      z[0] = z[n-1];
      z[n] = z[1];
    }

  if (periodic)
    {
      free (vv);
      free (uu);
      free (s);
    }
  free (beta);
  free (alpha);
  free (v);
  free (u);
  free (b);
  free (h);
}


/* interpolate() computes an approximate ordinate value for a given
   abscissa value, given an array of data points (stored in t[] and y[],
   containing abscissa and ordinate values respectively), and z[], the
   array of 2nd derivatives at the knots (i.e. internal data points).
   
   On entry, have n+1 >= 2 points in the t, y, z arrays, numbered 0..n.
   The number of knots (i.e. internal data points) is n-1; they are
   numbered 1..n-1 as in Kincaid and Cheney.  In the periodic case, the
   final knot, i.e., (t[n-1],y[n-1]), has the property that y[n-1]=y[0];
   also, y[n]=y[1].  The number of data points supplied by the user was n+1
   in the non-periodic case, and n in the periodic case.  When this
   function is called, n>=1 in the non-periodic case, and n>=2 in the
   periodic case. */

double
#ifdef _HAVE_PROTOS
interpolate (int n, double *t, double *y, double *z, double x, 
	     double tension, bool periodic)
#else
interpolate (n, t, y, z, x, tension, periodic)
     int n;
     double *t, *y, *z, x;
     double tension;
     bool periodic;
#endif
{
  double diff, updiff, reldiff, relupdiff, h;
  double value;
  int is_ascending = (t[n-1] < t[n]);
  int i = 0, k;

  /* in periodic case, map x to t[0] <= x < t[n] */
  if (periodic && (x - t[0]) * (x - t[n]) > 0.0)
    x -= ((int)(floor( (x - t[0]) / (t[n] - t[0]) )) * (t[n] - t[0]));

  /* do binary search to find interval */
  for (k = n - i; k > 1;)
    {
      if (is_ascending ? x >= t[i + (k>>1)] : x <= t[i + (k>>1)])
	{
	  i = i + (k>>1);
	  k = k - (k>>1);
	}
      else
	k = k>>1;
    }

  /* at this point, x is between t[i] and t[i+1] */
  h = t[i + 1] - t[i];
  diff = x - t[i];
  updiff = t[i+1] - x;
  reldiff = diff / h;
  relupdiff = updiff / h;

  if (tension == 0.0)
  /* evaluate cubic polynomial in nested form */
    value =  y[i] 
      + diff
	* ((y[i + 1] - y[i]) / h - h * (z[i + 1] + z[i] * 2.0) / 6.0
	   + diff * (0.5 * z[i] + diff * (z[i + 1] - z[i]) / (6.0 * h)));
  
  else if (tension > 0.0)
    /* `positive' (really real) tension, use sinh's */
    {
      if (fabs(tension * h) < TRIG_ARG_MIN)
	/* hand-compute (6/y^2)(sinh(xy)/sinh(y) - x) to improve accuracy;
	   here `x' means reldiff or relupdiff and `y' means tension*h */
	value = (y[i] * relupdiff + y[i+1] * reldiff
		 + ((z[i] * h * h / 6.0) 
		    * quotient_sinh_func (relupdiff, tension * h))
		 + ((z[i+1] * h * h / 6.0) 
		    * quotient_sinh_func (reldiff, tension * h)));
      else if (fabs(tension * h) > TRIG_ARG_MAX)
	/* approximate 1/sinh(y) by 2 sgn(y) exp(-|y|) */
	{
	  int sign = (h < 0.0 ? -1 : 1);

	  value = (((z[i] * (exp (tension * updiff - sign * tension * h) 
			     + exp (-tension * updiff - sign * tension * h))
		     + z[i + 1] * (exp (tension * diff - sign * tension * h) 
				   + exp (-tension * diff - sign * tension*h)))
		    * (sign / (tension * tension)))
		   + (y[i] - z[i] / (tension * tension)) * (updiff / h)
		   + (y[i + 1] - z[i + 1] / (tension * tension)) * (diff / h));
	}
      else
	value = (((z[i] * sinh (tension * updiff) 
		   + z[i + 1] * sinh (tension * diff))
		  / (tension * tension * sinh (tension * h)))
		 + (y[i] - z[i] / (tension * tension)) * (updiff / h)
		 + (y[i + 1] - z[i + 1] / (tension * tension)) * (diff / h));
    }
  else
    /* `negative' (really imaginary) tension, use sin's */
    {
      if (fabs(tension * h) < TRIG_ARG_MIN)
	/* hand-compute (6/y^2)(sin(xy)/sin(y) - x) to improve accuracy;
	   here `x' means reldiff or relupdiff and `y' means tension*h */
	value = (y[i] * relupdiff + y[i+1] * reldiff
		 + ((z[i] * h * h / 6.0) 
		    * quotient_sin_func (relupdiff, tension * h))
		 + ((z[i+1] * h * h / 6.0) 
		    * quotient_sin_func (reldiff, tension * h)));
      else
	value = (((z[i] * sin (tension * updiff) 
		   + z[i + 1] * sin (tension * diff))
		  / (tension * tension * sin (tension * h)))
		 + (y[i] - z[i] / (tension * tension)) * (updiff / h)
		 + (y[i + 1] - z[i + 1] / (tension * tension)) * (diff / h));
    }
  
  return value;
}


/* is_monotonic() check whether an array of data points, read in by
   read_data(), has monotonic abscissa values. */
bool
#ifdef _HAVE_PROTOS
is_monotonic (int n, double *t)
#else
is_monotonic (n, t)
     int n;			/* array size n+1, n>=1 */
     double *t;
#endif
{
  bool is_ascending;

  if (t[n-1] < t[n])
    is_ascending = true;
  else if (t[n-1] > t[n])
    is_ascending = false;
  else				/* equality */
    return false;

  while (n>0)
    {
      n--;
      if (is_ascending == true ? t[n] >= t[n+1] : t[n] <= t[n+1])
	return false;
    };
  return true;
}


/* read_float reads a single floating point quantity from an input file
   (in either ascii or double format).  Return value indicates whether it
   was read successfully. */
bool 
#ifdef _HAVE_PROTOS
read_float (FILE *input, double *dptr)
#else
read_float (input, dptr)
     FILE *input;
     double *dptr;
#endif
{
  int num_read;
  double dval;
  float fval;
  int ival;

  switch (input_type)
    {
    case T_ASCII:
    default:
      num_read = fscanf (input, "%lf", &dval);
      break;
    case T_SINGLE:
      num_read = fread ((voidptr_t) &fval, sizeof (fval), 1, input);
      dval = fval;
      break;
    case T_DOUBLE:
      num_read = fread ((voidptr_t) &dval, sizeof (dval), 1, input);
      break;
    case T_INTEGER:
      num_read = fread ((voidptr_t) &ival, sizeof (ival), 1, input);
      dval = ival;
      break;
    }
  if (num_read <= 0)
    return false;
  if (dval != dval)
    {
      fprintf (stderr, "%s: encountered a NaN (not-a-number) in binary input file, treating as EOF\n",
	       progname);
      return false;		/* effectively eof */
    }
  else
    {
      *dptr = dval;
      return true;
    }
}

/* Emit a pair of doubles, in specified output representation.  Inform user
   if any of the emitted values was out-of-bounds for single-precision or
   integer format. */
bool 
#ifdef _HAVE_PROTOS
write_point (double t, double *y, int ydimension, int precision, bool suppress_abscissa)
#else
write_point (t, y, ydimension, precision, suppress_abscissa)
     double t, *y;
     int ydimension, precision;
     bool suppress_abscissa;
#endif
{
  int i, num_written = 0;
  float ft, fy;
  int it, iy;

  switch (output_type)
    {
    case T_ASCII:
    default:
      if (suppress_abscissa == false)
	num_written += printf ("%.*g ", precision, t);
      for (i = 0; i < ydimension - 1; i++)
	num_written += printf ("%.*g ", precision, y[i]);
      num_written += printf ("%.*g\n", precision, y[ydimension - 1]);
      break;
    case T_SINGLE:
      if (suppress_abscissa == false)
	{
	  ft = FROUND(t);
	  if (ft == FLT_MAX || ft == -(FLT_MAX))
	    {
	      maybe_emit_oob_warning();
	      if (ft == FLT_MAX)
		ft *= 0.99999;	/* kludge */
	    }
	  num_written += fwrite ((voidptr_t) &ft, sizeof (ft), 1, stdout);
	}
      for (i = 0; i < ydimension; i++)
	{
	  fy = y[i];
	  if (fy == FLT_MAX || fy == -(FLT_MAX))
	    {
	      maybe_emit_oob_warning();
	      if (fy == FLT_MAX)
		fy *= 0.99999;	/* kludge */
	    }
	  num_written += fwrite ((voidptr_t) &fy, sizeof (fy), 1, stdout);
	}
      break;
    case T_DOUBLE:
      if (suppress_abscissa == false)
	num_written += fwrite ((voidptr_t) &t, sizeof (t), 1, stdout);
      for (i = 0; i < ydimension; i++)
	num_written += fwrite ((voidptr_t) &(y[i]), sizeof (double), 1, stdout);
      break;
    case T_INTEGER:
      if (suppress_abscissa == false)
	{
	  it = IROUND(t);
	  if (it == INT_MAX || it == -(INT_MAX))
	    {
	      maybe_emit_oob_warning();
	      if (it == INT_MAX)
		it--;
	    }
	  num_written += fwrite ((voidptr_t) &it, sizeof (it), 1, stdout);
	}
      for (i = 0; i < ydimension; i++)
	{
	  iy = IROUND(y[i]);
	  if (iy == INT_MAX || iy == -(INT_MAX))
	    {
	      maybe_emit_oob_warning();
	      if (iy == INT_MAX)
		iy--;
	    }
	  num_written += fwrite ((voidptr_t) &iy, sizeof (iy), 1, stdout);
	}
      break;
    }
  
  return (num_written > 0 ? true : false); /* i.e. return successp */
}

/* read_point() attempts to read a data point from an input file
   (auto-abscissa is supported, as are both ascii and double formats).
   Return value is 0 if a data point was read, 1 if no data point could be
   read (i.e. EOF or garbage in file).  A return value of 2 is special: it
   indicates that an explicit end-of-dataset indicator was seen in the input
   stream.  For an ascii stream this is two newlines in succession; for a
   double stream this is a DBL_MAX, etc. */
int
#ifdef _HAVE_PROTOS
read_point (FILE *input, double *t, double *y, int ydimension, 
	    bool *first_point,
	    int auto_abscissa, double *auto_t, double auto_delta, 
	    double *stored)
#else
read_point (input, t, y, ydimension, first_point, auto_abscissa, auto_t, auto_delta, stored)
     FILE *input;
     double *t, *y;
     int ydimension;
     bool *first_point;
     int auto_abscissa;
     double *auto_t, auto_delta;
     double *stored;
#endif
{
  bool success;
  int i, items_read, lookahead;

 head:

  if (input_type == T_ASCII)
    {
      bool two_newlines;

      /* skip whitespace, up to but not including 2nd newline */
      two_newlines = skip_whitespace (input);
      if (two_newlines)
	/* end-of-dataset indicator */
	return 2;
    }
  if (feof (input))
    return 1;

  if (input_type == T_ASCII)
    {
      lookahead = getc (input);
      ungetc (lookahead, input);
      if (lookahead == (int)'#')	/* comment line */
	{
	  char c;
	  
	  do 
	    {
	      items_read = fread (&c, sizeof (c), 1, input);
	      if (items_read <= 0)
		return 1;	/* EOF */
	    }
	  while (c != '\n');
	  ungetc ((int)'\n', input); /* push back \n at the end of # line */
	  goto head;
	}
    }

  if (auto_abscissa != AUTO_NONE) /* i.e. AUTO_INCREMENT or AUTO_BY_DISTANCE */
    {
      /* read 1st component of y */
      success = read_float (input, &(y[0]));
      if (!success)		/* e.g., EOF */
	return 1;
      if ((input_type == T_DOUBLE && y[0] == DBL_MAX)
	  || (input_type == T_SINGLE && y[0] == (double)FLT_MAX)
	  || (input_type == T_INTEGER && y[0] == (double)INT_MAX))
	/* end-of-dataset indicator */
	return 2;

      /* read other components of y */
      for (i = 1; i < ydimension; i++)
	{
	  success = read_float (input, &(y[i]));
	  if (!success)		/* effectively EOF (could be garbage) */
	    {
	      fprintf (stderr, "%s: input file terminated prematurely\n",
		       progname);
	      return 1;
	    }
	}

      /* t is kept track of, not read from file; two different methods */
      if (auto_abscissa == AUTO_INCREMENT)
	{
	  *t = *auto_t;
	  *auto_t += auto_delta;	/* update */
	}
      else			/* AUTO_BY_DISTANCE */
	{
	  if (*first_point == true)
	    {
	      *t = *auto_t;
	      *first_point = false;
	    }
	  else		/* compute distance to previous point */
	    {
	      double distsq = 0.0;

	      for (i = 0; i < ydimension; i++)
		distsq += (y[i] - stored[i])*(y[i] - stored[i]);
	      *auto_t += sqrt (distsq);
	      *t = *auto_t;
	    }
	  for (i = 0; i < ydimension; i++)	  
	    stored[i] = y[i];	/* store current point */
	}

      /* successfully read all components of y */
      return 0;
    }
  else
    {
      /* read t */
      success = read_float (input, t);
      if (!success)		/* e.g., EOF */
	return 1;
      if ((input_type == T_DOUBLE && *t == DBL_MAX)
	  || (input_type == T_SINGLE && *t == (double)FLT_MAX)
	  || (input_type == T_INTEGER && *t == (double)INT_MAX))
	/* end-of-dataset indicator */
	return 2;

      /* read components of y */
      for (i = 0; i < ydimension; i++)
	{
	  success = read_float (input, &(y[i]));
	  if (!success)		/* effectively EOF (could be garbage) */
	    {
	      fprintf (stderr, "%s: input file terminated prematurely\n",
		       progname);
	      return 1;
	    }
	}

      /* successfully read both t and all components of y */
      return 0;
    }
}

/* read_data() reads a single dataset from an input file, and stores it.
   If the stream is in ascii format, end-of-dataset is signalled by two
   newlines in succession.  If the stream is in double format,
   end-of-dataset is signalled by the occurrence of a DBL_MAX, etc.

   Return value is true if the dataset is ended by an explicit
   end-of-dataset, and false if the dataset is terminated by EOF.  That is,
   return value indicates whether another dataset is expected to follow. */
bool
#ifdef _HAVE_PROTOS
read_data (FILE *input, int *len, int *used, int auto_abscissa,
	   double auto_t, double auto_delta, 
	   double **t, int ydimension, double **y, double **z)
#else
read_data (input, len, used, auto_abscissa, auto_t, auto_delta, t, ydimension, y, z)
     FILE *input;
     int *len, *used;
     int auto_abscissa;
     double auto_t, auto_delta;
     double **t;
     int ydimension;
     double **y, **z;
#endif
{
  bool first = true;
  int i, success;
  double tt, *yy, *stored;

  yy = (double *)xmalloc (sizeof(double) * ydimension);
  stored = (double *)xmalloc (sizeof(double) * ydimension);
  for ( ; ; )
    {
      if ((++ *used) >= *len)
	{
	  *len *= 2;
	  *t = (double *)xrealloc (*t, sizeof(double) * *len);
	  for (i = 0; i < ydimension; i++)
	    {
	      y[i] = (double *)xrealloc (y[i], sizeof(double) * *len);
	      z[i] = (double *)xrealloc (z[i], sizeof(double) * *len);
	    }
	}

      success = read_point (input, &tt, yy, ydimension, &first,
			    auto_abscissa, &auto_t, auto_delta, stored);

      switch (success)
	{
	case 0:			/* good data point */
	  (*t)[*used] = tt;
	  for (i = 0; i < ydimension; i++)
	    y[i][*used] = yy[i];
	  break;
	case 1:			/* end of dataset, EOF seen */
	  (*used)--;
	  free (stored);
	  free (yy);
	  return false;
	case 2:			/* end of dataset, but input continues */
	  (*used)--;
	  free (stored);
	  free (yy);
	  return true;
	}
    }
}


/* do_spline() is the main routine for piecewise cubic spline
   interpolation, supporting both periodicity and a user-specified boundary
   condition parameter.  Nonzero tension may be specified, in which case
   the interpolate() routine, which this calls, will use not cubic
   polynomials but rather expressions involving hyperbolic sines.

   t[] and y[] are the arrays in which the abscissa and ordinate values of
   the user-specified data points are stored, and z[] is the array in which
   the 2nd derivatives at the knots (data points in the interior of the
   interval) will be stored.  used+1 is the effective size of each of these
   arrays.  The number of points supplied by the user was used+1 in the
   non-periodic case.  It was used+0 in the periodic case.  

   The reason that the number of elements is greater by one in the periodic
   case is that the first user-supplied data point occurs also at the end.
   In fact, in the periodic case this function will increment the size of
   the array once more, since the periodic interpolation algorithm requires
   the first two data points, not just the first, to appear at the end. */
void
#ifdef _HAVE_PROTOS
do_spline (int used, int len, double **t, int ydimension, double **y, double **z, 
	   double tension, bool periodic, bool spec_boundary_condition,
	   double k, int precision, double first_t, double last_t, 
	   double spacing_t, int no_of_intervals, bool spec_first_t, 
	   bool spec_last_t, bool spec_spacing_t, 
	   bool spec_no_of_intervals, bool suppress_abscissa)
#else
do_spline (used, len, t, ydimension, y, z, tension, periodic, spec_boundary_condition, k, precision, first_t, last_t, spacing_t, no_of_intervals, spec_first_t, spec_last_t, spec_spacing_t, spec_no_of_intervals, suppress_abscissa)
     int used;			/* used+1 elements stored in (*t)[] etc. */
     int len;			/* length of each array */
     double **t;
     int ydimension;
     double **y, **z;		/* we use ** because may have to realloc */
     double tension;
     bool periodic;
     bool spec_boundary_condition;
     double k;			/* boundary condition: y''_1 = k y''_0, etc. */
     int precision;
     double first_t, last_t, spacing_t; 
     int no_of_intervals;
     bool spec_first_t, spec_last_t, spec_spacing_t, spec_no_of_intervals;
     bool suppress_abscissa;
#endif
{
  int range_count = 0;		/* number of req'd datapoints out of range */
  int lastval = 0;		/* last req'd point = 1st/last data point? */
  int i;

  if (used + 1 == 0)		/* zero data points in array */
    /* don't output anything (i.e. effectively output a null dataset) */
    return;

  if (used+1 == 1)		/* a single data point in array */
    {
      fprintf (stderr, 
	       "%s: cannot construct a spline from a single data point\n", 
	       progname);
      /* don't output anything (i.e. effectively output a null dataset) */
      return;
    }

  if (!periodic && used+1 <= 2)
    {
      if (spec_boundary_condition)
	fprintf (stderr, 
		 "%s: only 2 data points, so ignoring specified boundary condition\n", 
		 progname);
      k = 0.0;
    }

  if (!is_monotonic (used, *t))
    non_monotonic_error();	/* self-explanatory */

  if (periodic)
    {
      bool print_warning = false;
      
      for (i = 0; i < ydimension; i++)
	{
	  if (y[i][used] != y[i][0])
	    print_warning = true;
	  y[i][used] = y[i][0];
	}
      if (print_warning == true)
	fprintf (stderr, "%s: setting final y value equal to initial to ensure periodicity\n", 
		 progname); 

      /* add pseudo-point at end (to accord with periodicity) */
      if (used + 1 >= len)
	{
	  len++;
	  *t = (double *)xrealloc (*t, sizeof(double) * len);
	  for (i = 0; i < ydimension; i++)
	    {
	      y[i] = (double *)xrealloc (y[i], sizeof(double) * len);
	      z[i] = (double *)xrealloc (z[i], sizeof(double) * len);
	    }
	}
      (*t)[used + 1] = (*t)[used] + ((*t)[1] - (*t)[0]);
      for (i = 0; i < ydimension; i++)
	y[i][used + 1] = y[i][1];
    }

  /* compute z[], array of 2nd derivatives at each knot */
  for (i = 0; i < ydimension; i++)
    fit (used + (periodic ? 1 : 0), /* include pseudo-point if any */
	 *t, y[i], z[i], k, tension, periodic);

  if (!spec_first_t) 
    first_t = (*t)[0];
  if (!spec_last_t) 
    last_t = (*t)[used];	/* used+1 data points in all */
  if (!spec_spacing_t) 
    {
      if (no_of_intervals > 0)
	spacing_t = (last_t - first_t) / no_of_intervals;
      else
	spacing_t = 0;		/* won't happen */
    }
  else				/* user specified spacing */
    {
      if ((last_t - first_t) * spacing_t < 0.0)
	{
	  fprintf (stderr, "%s: specified spacing is of wrong sign, corrected\n",
		   progname);
	  spacing_t = -spacing_t;
	}
      if (spec_no_of_intervals)
	fprintf (stderr, "%s: ignoring specified number of intervals\n",
		 progname);
      no_of_intervals = (int)(fabs((last_t - first_t) / spacing_t) + FUZZ);
    }

  if (last_t == (*t)[0])
    lastval = 1;
  else if (last_t == (*t)[used])
    lastval = 2;

  for (i = 0; i <= no_of_intervals; ++i)
    {
      double x;

      x = first_t + spacing_t * i;

      if (i == no_of_intervals)
	{
	  /* avoid numerical fuzz */
	  if (lastval == 1)	/* left end of input */
	    x = (*t)[0];
	  else if (lastval == 2) /* right end of input */
	    x = (*t)[used];
	}

      if (periodic || (x - (*t)[0]) * (x - (*t)[used]) <= 0)
	{
	  int j;
	  double *yy;

	  yy = (double *)xmalloc (sizeof(double) * ydimension); 
	  for (j = 0; j < ydimension; j++)
	    yy[j] = interpolate (used, *t, y[j], z[j], x, 
				 tension, periodic);
	  write_point (x, yy, ydimension, precision, suppress_abscissa);
	  free (yy);
	}
      else
	range_count++;
    }

  switch (range_count)
    {
    case 0:
      break;
    case 1:
      fprintf (stderr, 
	       "%s: one requested point could not be computed (out of data range)\n", 
	       progname);
      break;
    default:
      fprintf (stderr, 
	       "%s: %d requested points could not be computed (out of data range)\n", 
	       progname, range_count);
      break;
    }
}


/* do_bessel() is the main routine for doing real-time cubic Bessel
   interpolation of a dataset.  If the input stream is in ascii format,
   end-of-dataset is signalled by two newlines in succession.  If the
   stream is in double format, end-of-dataset is signalled by the
   occurrence of a DBL_MAX, etc.

   Return value is true if the dataset is ended by an explicit
   end-of-dataset, and false if the dataset is terminated by EOF.  That is,
   return value indicates whether another dataset is expected to follow. */
bool
#ifdef _HAVE_PROTOS
do_bessel (FILE *input, int ydimension, int auto_abscissa, double auto_t, 
	   double auto_delta, double first_t, double last_t, 
	   double spacing_t, int precision, bool suppress_abscissa)
#else
do_bessel (input, ydimension, auto_abscissa, auto_t, auto_delta, first_t, last_t, spacing_t, precision, suppress_abscissa)
     FILE *input;
     int ydimension;
     int auto_abscissa;
     double auto_t, auto_delta;     
     double first_t, last_t, spacing_t;
     int precision;
     bool suppress_abscissa;
#endif
{
  bool first = true;
  double t, *y, *s0, *s1, *s2, *stored;
  double tt[4], **yy;
  int direction = (last_t > first_t ? 1 : -1);
  int state = STATE_ZERO;
  int i, success;

  y = (double *)xmalloc (sizeof(double) * ydimension); 
  s0 = (double *)xmalloc (sizeof(double) * ydimension); 
  s1 = (double *)xmalloc (sizeof(double) * ydimension); 
  s2 = (double *)xmalloc (sizeof(double) * ydimension); 
  yy = (double **)xmalloc (4 * sizeof(double *));
  stored = (double *)xmalloc (sizeof(double) * ydimension);
  for (i = 0; i < 4; i++)
    yy[i] = (double *)xmalloc (ydimension * sizeof(double));

  for ( ; ; )
    {
      success = read_point (input, &t, y, ydimension, &first,
			    auto_abscissa, &auto_t, auto_delta, stored);
      
      if (success == 0)		/* got a new data point */
	{
	  /* use our DFA to process the new data point */
	  switch (state)
	    {
	    case STATE_ZERO:	/* just store point */
	      tt[0] = t;
	      for (i = 0; i < ydimension; i++)
		yy[0][i] = y[i];
	      state = STATE_ONE;
	      break;
	    case STATE_ONE:	/* just store point */
	      tt[1] = t;
	      if (direction * (tt[1] - tt[0]) <= 0)
		non_monotonic_error();
	      for (i = 0; i < ydimension; i++)
		yy[1][i] = y[i];
	      state = STATE_TWO;
	      break;
	    case STATE_TWO:	/* store point, and process */
	      tt[2] = t;
	      if (direction * (tt[2] - tt[1]) <= 0)
		non_monotonic_error();
	      for (i = 0; i < ydimension; i++)
		{
		  yy[2][i] = y[i];
		  
		  /* fit parabola through 0,1,2 to compute slopes at 0,1*/
		  s0[i] = (((tt[1]-tt[0]) * ((yy[0][i]-yy[2][i]) / (tt[0]-tt[2]))
			 + (tt[0]-tt[2]) * ((yy[1][i]-yy[0][i]) / (tt[1]-tt[0])))
			/ (tt[1]-tt[2]));
		  s1[i] = (((tt[2]-tt[1]) * ((yy[1][i]-yy[0][i]) / (tt[1]-tt[0]))
			 + (tt[1]-tt[0]) * ((yy[2][i]-yy[1][i]) / (tt[2]-tt[1])))
			/ (tt[2]-tt[0]));
		}

	      /* output spline points in range between points 0, 1 */
	      do_bessel_range (tt[0], tt[1], yy[0], yy[1], s0, s1,
			       first_t, last_t, spacing_t, 
			       ydimension, precision, false,
			       suppress_abscissa);
	      
	      state = STATE_THREE;
	      break;
	    case STATE_THREE:	/* store point, and process */
	      tt[3] = t;
	      if (direction * (tt[3] - tt[2]) <= 0)
		non_monotonic_error();
	      for (i = 0; i < ydimension; i++)
		{
		  yy[3][i] = y[i];
		  
		  /* fit parabola through points 1,2,3 to compute slope at 2 */
		  s2[i] = (((tt[3]-tt[2]) * ((yy[2][i]-yy[1][i]) / (tt[2]-tt[1]))
			 + (tt[2]-tt[1]) * ((yy[3][i]-yy[2][i]) / (tt[3]-tt[2])))
			/ (tt[3]-tt[1]));
		}
	      
	      /* output spline points in range between points 1, 2 */
	      do_bessel_range (tt[1], tt[2], yy[1], yy[2], s1, s2, 
			       first_t, last_t, spacing_t, 
			       ydimension, precision, false,
			       suppress_abscissa);
	      
	      /* shift points down */
	      tt[0] = tt[1];
	      tt[1] = tt[2];
	      tt[2] = tt[3];
	      for (i = 0; i < ydimension; i++)
		{
		  yy[0][i] = yy[1][i];
		  yy[1][i] = yy[2][i];
		  yy[2][i] = yy[3][i];
		  /* shift down the only knot slope worth keeping */
		  s1[i] = s2[i];
		}

	      break;
	    }
	}
      else		/* didn't get a point, so wind things up */
	{
	  switch (state)
	    {
	    case STATE_ZERO:
	      /* silently output a null dataset (i.e., don't output anything) */
	      break;
	    case STATE_ONE:
	      fprintf (stderr, "%s: cannot construct a spline from a single data point\n", 
		       progname);
	      /* output a null dataset (i.e., don't output anything) */
	      break;
	    case STATE_TWO:
	      /* have two points: do linear interp between points 0, 1 */
	      for (i = 0; i < ydimension; i++)
		s0[i] = s1[i] = (yy[1][i] - yy[0][i])/(tt[1]-tt[0]);
	      do_bessel_range (tt[0], tt[1], yy[0], yy[1], s0, s1, 
			       first_t, last_t, spacing_t, 
			       ydimension, precision, true,
			       suppress_abscissa);
	      break;
	    case STATE_THREE:
	      /* already did 1st of 2 intervals, so do 2nd one too */

	      /* fit parabola through points 0,1,2 to compute slope at 2 */
	      for (i = 0; i < ydimension; i++)
		s2[i] = (((tt[0]-tt[2]) * ((yy[2][i]-yy[1][i]) / (tt[2]-tt[1]))
		       + (tt[2]-tt[1]) * ((yy[0][i]-yy[2][i]) / (tt[0]-tt[2])))
		      / (tt[0]-tt[1]));

	      /* output spline points in range between points 1, 2 */
	      do_bessel_range (tt[1], tt[2], yy[1], yy[2], s1, s2, 
			       first_t, last_t, spacing_t, 
			       ydimension, precision, true,
			       suppress_abscissa);
	      break;
	    }

	  /* free storage before return */
	  for (i = 0; i < 4; i++)
	    free (yy[i]);
	  free (stored);
	  free (yy);
	  free (s2);
	  free (s1);
	  free (s0);
	  free (y);

	  /* return indication of whether end-of-dataset was seen in stream */
	  return (success == 2 ? true : false);
	}
    }
}

void
#ifdef _HAVE_PROTOS
non_monotonic_error (void)
#else
non_monotonic_error ()
#endif
{
  fprintf (stderr, "%s: error: abscissa values not monotonic\n",
	   progname);
  exit (EXIT_FAILURE);
}


/* do_bessel_range() computes spline points separated by spacing_t, within
   the abscissa interval abscissa0 <= t < abscissa1, that happen to lie in
   the desired range first_t <= t <= last_t.  It writes them to standard
   output.  The ordinate values value0 and value1, and endpoint slopes
   slope0 and slope1, are specified.  If `endit' is set, then the intervals
   stretch slightly farther than abscissa1 and last_t, to compensate for
   roundoff error. */

void
#ifdef _HAVE_PROTOS
do_bessel_range (double abscissa0, double abscissa1, double *value0, 
		 double *value1, double *slope0, double *slope1, 
		 double first_t, double last_t, double spacing_t, 
		 int ydimension, int precision, bool endit,
		 bool suppress_abscissa)
#else
do_bessel_range (abscissa0, abscissa1, value0, value1, slope0, slope1,
		 first_t, last_t, spacing_t, ydimension, precision, 
		 endit, suppress_abscissa)
     double abscissa0, abscissa1, *value0, *value1, *slope0, *slope1;
     double first_t, last_t, spacing_t;
     int ydimension;
     int precision;
     bool endit;		/* last interval to be treated */
     bool suppress_abscissa;
#endif
{
  int direction = ((last_t > first_t) ? 1 : -1); /* sign of spacing_t */
  int i, j;
  int imin1 = (int)((abscissa0 - first_t) / spacing_t - 1);
  int imax1 = (int)((abscissa1 - first_t) / spacing_t + 1);
  int imin2 = 0;
  int imax2 = (int)((last_t - first_t) / spacing_t + 1);
  int imin, imax;
  
  /* compute maximum interval over which i must range */
  imin = IMAX (imin1, imin2);
  imax = IMIN (imax1, imax2);
  for (i = imin; i <= imax; i++)
    {
      double t;

      t = first_t + i * spacing_t;

      if ((direction * t >= direction * abscissa0)
	  && (direction * t >= direction * first_t)
	  /* stretch slightly if `endit' is set */
	  && ((direction * t < (direction 
				* (abscissa1 
				   + (endit ? 
				      FUZZ * (abscissa1 - abscissa0) : 0.)))))
	  && (direction * t <= (direction
			       * (last_t
				  + (endit ? FUZZ * (last_t - first_t) : 0.)))))
	{
	  double diff = t - abscissa0;
	  double updiff = abscissa1 - t;
	  double h = abscissa1 - abscissa0;
	  double *y;
	  bool success;

	  y = (double *)xmalloc (sizeof(double) * ydimension); 
	  for (j = 0; j < ydimension; j++)
	    {
	      /* should use a nested form */
	      y[j] = (value1[j] * (-2.0 * diff * diff * diff / (h * h * h)
				   + 3.0 * diff * diff / (h * h))
		+ value0[j] * (-2.0 * updiff * updiff * updiff / (h * h * h)
			           + 3.0 * updiff * updiff / (h * h)))
		+ ((slope1[j] * (diff * diff * diff / (h * h) 
			      - diff * diff / h)
		- (slope0[j] * (updiff * updiff * updiff / (h * h) 
				 - updiff * updiff / h))));
	    }
	  
	  success = write_point (t, y, 
				 ydimension, precision, suppress_abscissa);
	  if (!success)
	    {
	      fprintf (stderr, 
		       "%s: error: unable to write to standard output\n",
		       progname);
	      exit (EXIT_FAILURE);
	    }
	  free (y);
	}	  
    }
}


/* Output a separator between datasets.  For ascii-format output streams
   this is an extra newline (after the one that the spline ended with,
   yielding two newlines in succession).  For double-format output streams
   this is a DBL_MAX, etc. */

void
#ifdef _HAVE_PROTOS
output_dataset_separator (void)
#else
output_dataset_separator ()
#endif
{
  double ddummy;
  float fdummy;
  int idummy;

  switch (output_type)
    {
    case T_ASCII:
    default:
      printf ("\n");
      break;
    case T_DOUBLE:
      ddummy = DBL_MAX;
      fwrite ((voidptr_t) &ddummy, sizeof(ddummy), 1, stdout);
      break;
    case T_SINGLE:
      fdummy = FLT_MAX;
      fwrite ((voidptr_t) &fdummy, sizeof(fdummy), 1, stdout);
      break;
    case T_INTEGER:
      idummy = INT_MAX;
      fwrite ((voidptr_t) &idummy, sizeof(idummy), 1, stdout);
      break;
    }
}

/* skip_whitespace() skips whitespace in an ascii-format input file,
   up to but not including a second newline.  Return value indicates
   whether or not two newlines were in fact seen.  (For ascii-format
   input files, two newlines signals an end-of-dataset.) */

bool
#ifdef _HAVE_PROTOS
skip_whitespace (FILE *stream)
#else
skip_whitespace (stream)
     FILE *stream;
#endif
{
  int lookahead;
  int nlcount = 0;
  
  do 
    {
      lookahead = getc (stream);
      if (lookahead == (int)'\n')
	  nlcount++;
    }
  while (lookahead != EOF 
	 && isspace((unsigned char)lookahead)
	 && nlcount < 2);

  if (lookahead == EOF)
    return false;
  
  ungetc (lookahead, stream);
  return (nlcount == 2 ? true : false);
}

void
#ifdef _HAVE_PROTOS
maybe_emit_oob_warning (void)
#else
maybe_emit_oob_warning ()
#endif
{
  static bool warning_written = false;

  if (!warning_written)
    {
      fprintf (stderr, "%s: approximating one or more out-of-bounds output values\n", progname);
      warning_written = true;
    }
}


/* Following four functions compute (6/x^2)(1-x/sinh(x)),
   (3/x^2)(x/tanh(x)-1), (6/x^2)(1-x/sin(x)), and (3/x^2)(x/tan(x)-1) via
   the first three terms of the appropriate power series.  They are used
   when |x|<TRIG_ARG_MIN, to avoid loss of significance.  Errors are
   O(x**6). */
double
#ifdef _HAVE_PROTOS
sinh_func (double x) 
#else
sinh_func (x)
     double x;
#endif
{
  /* use 1-(7/60)x**2+(31/2520)x**4 */
  return 1.0 - (7.0/60.0)*x*x + (31.0/2520.0)*x*x*x*x;
}

double
#ifdef _HAVE_PROTOS
tanh_func (double x) 
#else
tanh_func (x)
     double x;
#endif
{
  /* use 1-(1/15)x**2+(2/315)x**4 */
  return 1.0 - (1.0/15.0)*x*x + (2.0/315.0)*x*x*x*x;
}

double
#ifdef _HAVE_PROTOS
sin_func (double x) 
#else
sin_func (x)
     double x;
#endif
{
  /* use -1-(7/60)x**2-(31/2520)x**4 */
  return -1.0 - (7.0/60.0)*x*x - (31.0/2520.0)*x*x*x*x;
}

double
#ifdef _HAVE_PROTOS
tan_func (double x) 
#else
tan_func (x)
     double x;
#endif
{
  /* use -1-(1/15)x**2-(2/315)x**4 */
  return -1.0 - (1.0/15.0)*x*x - (2.0/315.0)*x*x*x*x;
}


/* Following two functions compute (6/y^2)(sinh(xy)/sinh(y)-x) and
   (6/y^2)(sin(xy)/sin(y)-x), via the first three terms of the appropriate
   power series in y.  They are used when |y|<TRIG_ARG_MIN, to avoid loss
   of significance.  Errors are O(y**6). */
double
#ifdef _HAVE_PROTOS
quotient_sinh_func (double x, double y) 
#else
quotient_sinh_func (x, y)
     double x, y;
#endif
{
  return ((x*x*x-x) + (x*x*x*x*x/20.0 - x*x*x/6.0 + 7.0*x/60.0)*(y*y)
	  + (x*x*x*x*x*x*x/840.0 - x*x*x*x*x/120.0 + 7.0*x*x*x/360.0
	     -31.0*x/2520.0)*(y*y*y*y));
}

double
#ifdef _HAVE_PROTOS
quotient_sin_func (double x, double y) 
#else
quotient_sin_func (x, y)
     double x, y;
#endif
{
  return (- (x*x*x-x) + (x*x*x*x*x/20.0 - x*x*x/6.0 + 7.0*x/60.0)*(y*y)
	  - (x*x*x*x*x*x*x/840.0 - x*x*x*x*x/120.0 + 7.0*x*x*x/360.0
	     -31.0*x/2520.0)*(y*y*y*y));
}