/* spline.c - curve smoothing ported from edition 7 source */
#include <u.h>
#include <libc.h>
#include <bio.h>

#define NP 1000		/* max #  of points read */
#define INF 1.e37	/* default upper and lower limits */

struct proj { 
	int	lbf, ubf; 
	float	a, b, lb, ub, quant, mult, val[NP]; 
} x, y;

Biobuf	bout;
float	*diag, *r;
float	dx = 1.0;
float	ni = 100.0;
int	n;
int	auta;
int	periodic;
float	konst = 0.0;
float	zero = 0.0;

float rhs(int);
int spline(void);
void readin(Biobuf *bp);
void getlim(struct proj *);
void usage(void);

float
rhs(int i)
{
	int	i_;
	double	zz;
	i_ = i == n - 1 ? 0 : i;
	zz = (y.val[i] - y.val[i-1]) / (x.val[i] - x.val[i-1]);
	return(6 * ((y.val[i_+1] - y.val[i_]) / (x.val[i+1] - x.val[i]) - zz));
}


int
spline(void)
{
	float	d, s, u, v, hi, hi1;
	float	h;
	float	D2yi, D2yi1, D2yn1, x0, x1, yy, a;
	int	end;
	float	corr;
	int	i, j, m;
	if (n < 3) 
		return(0);
	if (periodic) 
		konst = 0;
	d = 1;
	r[0] = 0;
	s = periodic ? -1 : 0;
	for (i = 0; ++i < n - !periodic; ) {	/* triangularize */
		hi = x.val[i] - x.val[i-1];
		hi1 = i == n - 1 ? x.val[1] - x.val[0] : 
		    x.val[i+1] - x.val[i];
		if (hi1 * hi <= 0) 
			return(0);
		u = i == 1 ? zero : u - s * s / d;
		v = i == 1 ? zero : v - s * r[i-1] / d;
		r[i] = rhs(i) - hi * r[i-1] / d;
		s = -hi * s / d;
		a = 2 * (hi + hi1);
		if (i == 1) 
			a += konst * hi;
		if (i == n - 2) 
			a += konst * hi1;
		diag[i] = d = i == 1 ? a : 
		    a - hi * hi / d;
	}
	D2yi = D2yn1 = 0;
	for (i = n - !periodic; --i >= 0; ) {	/* back substitute */
		end = i == n - 1;
		hi1 = end ? x.val[1] - x.val[0] : 
		    x.val[i+1] - x.val[i];
		D2yi1 = D2yi;
		if (i > 0) {
			hi = x.val[i] - x.val[i-1];
			corr = end ? 2 * s + u : zero;
			D2yi = (r[i] - hi1 * D2yi1 - s * D2yn1 + end * v) / 
			    (diag[i] + corr);
			if (end) 
				D2yn1 = D2yi;
			if (i > 1) {
				a = 2 * (hi + hi1);
				if (i == 1) 
					a += konst * hi;
				if (i == n - 2) 
					a += konst * hi1;
				d = diag[i-1];
				s = -s * d / hi;
			}
		} else 
			D2yi = D2yn1;

		if (!periodic) {
			if (i == 0) 
				D2yi = konst * D2yi1;
			if (i == n - 2) 
				D2yi1 = konst * D2yi;
		}
		if (end) 
			continue;
		m = (hi1 > 0) ? ni : -ni;
		m = 1.001 * m * hi1 / (x.ub - x.lb);
		if (m <= 0) 
			m = 1;
		h = hi1 / m;
		for (j = m; j > 0 || i == 0 && j == 0; j--) {	/* interpolate */
			x0 = (m - j) * h / hi1;
			x1 = j * h / hi1;
			yy = D2yi * (x0 - x0 * x0 * x0) + D2yi1 * (x1 - x1 * x1 * x1);
			yy = y.val[i] * x0 + y.val[i+1] * x1 - hi1 * hi1 * yy / 6;
			Bprint(&bout, "%f ", x.val[i] + j * h);
			Bprint(&bout, "%f\n", yy);
		}
	}
	return(1);
}


void
readin(Biobuf *bp) 
{
	int num;
	char *p, *a[4];

	for (n = 0; n < NP && (p = Brdline(bp, '\n')) != nil; n++){
		p[Blinelen(bp)-1] = 0;
		if ((num = tokenize(p, a, nelem(a))) == 0)
			continue;
		if (num == 1 || auta){
			x.val[n] = n * dx * x.lb;
			y.val[n] = atof(a[0]);
		}
		if(num >= 2){
			x.val[n] = atof(a[0]);
			y.val[n] = atof(a[1]);
		}
	}
}


void
getlim(struct proj *p)
{
	int	i;
	for (i = 0; i < n; i++) {
		if (!p->lbf && p->lb > (p->val[i])) 
			p->lb = p->val[i];
		if (!p->ubf && p->ub < (p->val[i])) 
			p->ub = p->val[i]; 
	}
}


void
main(int argc, char *argv[])
{
	int i;
	char *p;
	Biobuf *bp, bin;

	x.lbf = x.ubf = y.lbf = y.ubf = 0;
	x.lb = INF;
	x.ub = -INF;
	y.lb = INF;
	y.ub = -INF;

	Binit(&bin, 0, OREAD);
	Binit(&bout, 1, OWRITE);

	ARGBEGIN{
		case 'a': 
			auta = 1;
			dx = 1;
			if ((p = ARGF()) != nil)
				dx = atof(p);
			break;
		case 'k': 
			konst = atof(EARGF(usage()));
			break;
		case 'n': 
			ni = atof(EARGF(usage()));
			break;
		case 'p':
			periodic = 1;
			break;
		case 'x':
			x.lbf = atof(EARGF(usage()));
			if ((p = ARGF()) != nil)
				x.ubf = atof(p);
			break;
	}ARGEND;
	if (auta && !x.lbf) 
		x.lb = 1;

	if (argc == 0)
		readin(&bin);
	else
		for(; argc--; argv++){
			if ((bp = Bopen(*argv, OREAD)) == nil){
				fprint(2, "%s: %s cannot open, %r", argv0, *argv);
				continue;
			}
			readin(bp);
			Bterm(bp);
		}

	getlim(&x);
	getlim(&y);

	if ((diag = malloc((n+1) * sizeof(dx))) == nil)
		sysfatal("no memory");
	if ((r = malloc((n+1) * sizeof(dx))) == nil)
		sysfatal("no memory");

	if (spline() == 0) 
		for (i = 0; i < n; i++) {
			Bprint(&bout, "%f ", x.val[i]);
			Bprint(&bout, "%f\n", y.val[i]); 
		}
	Bflush(&bout);
	Bterm(&bout);
	exits(0);
}

void
usage(void)
{
	fprint(2, "usage: %s [-a n] [-k n] [-n n] [-p] [-x lo [hi]]\n", argv0);
	exits("usage");
}