/* interp - spline.c */

/*	This module handles spline interpolation of Toby
 *	Orloff's interpolation program.
 *
 *	JSL (Jim Larson) 6/26/90
 */

/*
 * Copyright (c) 1990, Geometry Supercomputer Project
 *                     University of Minnesota
 *                     1200 Washington Ave. S
 *                     Minneapolis, MN  55415
 *
 * email address: software@geom.umn.edu
 *
 * This software is copyrighted as noted above.  It may be freely copied,
 * modified, and redistributed, provided that the copyright notice is
 * preserved on all copies.
 *
 * There is no warranty or other guarantee of fitness for this software,
 * it is provided solely "as is".  Bug reports or fixes may be sent
 * to the authors, who may or may not act on them as they desire.
 *
 * You may not include this software in a program or other software product
 * without supplying the source, or without informing the end-user that the
 * source is available for no extra charge.
 *
 * If you modify this software, you should include a notice giving the
 * name of the person performing the modification, the date of modification,
 * and the reason for such modification.
 */


#include <stdio.h>
#include <math.h>
#include "types.h"


/*	See Bartels, Beatty, and Barsky, p. 426.
 *
 *	To interpolate a given segment, we need to know the values
 *	and slopes of the endpoints.  These determine a unique cubic
 *	polynomial function that meets those endpoint conditions.
 *	An efficient way of calculating the interpolated points is to
 *	make use of Hermitian interpolation polynomials:
 *		h00(x) = 2x^3 - 3x^2 + 1,
 *		h10(x) = -2x^3 + 3x^2,
 *		h01(x) = x^3 - 2x^2 + x,
 *		h11(x) = x^3 - x^2.
 *	These polynomials have the property that the a-th derivative of
 *	h<i><j>(b) = ((i == b) and (j == a)) where a,b,i,j = 0 or 1.
 */

#define h00	(2*x3 - 3*x2 + 1)
#define h10	(-2*x3 + 3*x2)
#define h01	(x3 - 2*x2 + x)
#define h11	(x3 - x2)

/*	Suppose we have a pair of points (x0, y0) and (x1, y1), and we
 *	want a cubic spline with slopes d0 and d1 at the respective end-
 *	points.  The unique cubic function can be calculated as
 *		len = x1 - x0;
 *		g(x) = (x - x0)/len;
 *		f(x) = y0 * h00(g(x)) + y1 * h10(g(x)) +
 *		       len * (d0 * h10(g(x)) + h11(g(x)));
 *
 *	If a number of segments are being computed with the same number
 *	of divisions, the values hnm(x) can be saved for the different
 *	values of x, and then merely multiplied with the values and
 *	derivatives to get the interpolated values.
 */

/*	There are several ways to determine the slopes of the data
 *	points.  One method is to set the slopes so that the whole
 *	interpolated curve has a continuous second derivative, but
 *	this has the undesired effect that there is no local control:
 *	changes to the value of one data point affect the entire curve.
 *	A second method employed here is to make the slope at point
 *	x_i equal to the slope of the chord from x_(i-1) to x_(i+1), or
 *		d_i = (y_(i+1) - y(i-1))/(x_(i+1) - x(i-1)).
 *	This method satisfies the criterion of local control.
 */

/*	To find a derivative for the endpoints, we will choose a value
 *	that will make the second derivative of the spline zero at the
 *	endpoint.  To do this, just use the simple formula
 *		d_unknown = (3 * (y_1 - y_0)/(len) - d_other)/2.
 *	This formula works for both left and right endpoints.
 *	Thus, to use cubic spline interpolation, we need at least
 *	three data points.
 */

#define ERR_MESG "Spline interpolation requires at least three data points.\n"

/*	To generate new data outside of the data points given, we need
 *	to extrapolate rather than interpolate.  The most reasonable
 *	way to proceed is to use the same function as the last spline.
 */

/*	Here is the main routine called by the interp.c module.
 */

spline_interp (in, out, attr)
     FileValueList *in, *out;
     InterpolationAttributes *attr;
{
  FileValueList *infile, *outfile;
  int i, num;
  static int files_ok = 0;

  num = attr->num;

  if (!files_ok)
    for (i = 0, infile = in; i <3; i++, infile = infile->next)
      if (!infile)
	{
	  fprintf (stderr, ERR_MESG);
	  return;
	}

  /* Interpolate values for each requested output time. */

  for (outfile = out; outfile; outfile = outfile->next)
    {
      double lderiv, rderiv;
      ValueList *lvalue, *rvalue, *llvalue, *rrvalue, *outvalue;
      double time, lltime, ltime, rtime, rrtime;
      double len, x, x2, x3;
      int llexists, rrexists;	/* Do surrounding data points exist? */

      /* Initialize values for this output file. */

      time = outfile->time;
      infile = in;
      outvalue = outfile->value.next;
      llexists = 0;
      rrexists = 1;
      ltime = infile->time;
      lvalue = infile->value.next;

      /* Check to see if we hit first data point. */

      if (ltime == time)
	{
	  for (i = 0; i < num; i++)
	    {
	      outvalue->value = lvalue->value;
	      lvalue = lvalue->next;
	      outvalue = outvalue->next;
	    }
	  continue;		/* Go to next outfile. */
	}

      infile = infile->next;
      rtime = infile->time;
      rvalue = infile->value.next;

      do
	{
	  if (rtime == time)	/* Check to see if we hit this data point. */
	    {
	      for (i = 0; i < num; i++)
		{
		  outvalue->value = rvalue->value;
		  rvalue = rvalue->next;
		  outvalue = outvalue->next;
		}
	      break;
	    }

	  if (rtime > time	/* We are at the correct segment. */
	      || !(infile->next)) /* We need to extrapolate the right side. */
	    {
	      len = rtime - ltime;
	      x = (time - ltime)/len;
	      x2 = x * x;
	      x3 = x2 * x;

	      if (!llexists)	/* We are at the first segment. */
		{
		  rrexists = 1;
		  rrtime = infile->next->time;
		  rrvalue = infile->next->value.next;
		}
	      else if (infile->next) /* Check for values to the right. */
		{
		  rrexists = 1;
		  rrtime = infile->next->time;
		  rrvalue = infile->next->value.next;
		}
	      else		/* No values to the right. */
		rrexists = 0;

	      for (i = 0; i < num; i++)
		{
		  /* Calculate derivatives at the endpoints. */

		  if (!llexists)	/* We are at the first segment. */
		    {
		      rderiv = (rrvalue->value - lvalue->value)
			/(rrtime - ltime);
		      lderiv = (3 * (rvalue->value - lvalue->value)/len
				- rderiv)/2.0;
		    }
		  else if (!rrexists) /* We are at the last segment. */
		    {
		      lderiv = (rvalue->value - llvalue->value)
			/(rtime - lltime);
		      rderiv = (3 * (rvalue->value - lvalue->value)/len
				- lderiv)/2.0;
		    }
		  else		/* We are in a middle segment. */
		    {
		      rderiv = (rrvalue->value - lvalue->value)
			/(rrtime - ltime);
		      lderiv = (rvalue->value - llvalue->value)
			/(rtime - lltime);
		    }

		  outvalue->value = lvalue->value * h00 + rvalue->value * h10
		    + lderiv * h01 + rderiv * h11;

		  outvalue = outvalue->next;
		  if (llexists)
		    llvalue = llvalue->next;
		  lvalue = lvalue->next;
		  rvalue = rvalue->next;
		  if (rrexists)
		    rrvalue = rrvalue->next;
		}
	      break;
	    }

	  /* This is the wrong segment; continue down the infile list. */

	  lltime = ltime;
	  llvalue = lvalue;
	  llexists = 1;
	  ltime = rtime;
	  lvalue = rvalue;
	  infile = infile->next;
	  rtime = infile->time;
	  rvalue = infile->value.next;
	}
      while (1);
    }				/* Now go to next outfile. */
}

