[unix-history] / usr / src / usr.bin / spline / spline.c

static char *sccsid = "@(#)spline.c	4.2 (Berkeley) %G%";
#include <stdio.h>
#include <math.h>

#define NP 1000
#define INF HUGE

struct proj { int lbf,ubf; float a,b,lb,ub,quant,mult,val[NP]; } x,y;
float *diag, *r;
float dx = 1.;
float ni = 100.;
int n;
int auta;
int periodic;
float konst = 0.0;
float zero = 0.;

/* Spline fit technique
let x,y be vectors of abscissas and ordinates
    h   be vector of differences h\e9i\e8=x\e9i\e8-x\e9i-1\e\e9\e8\e8
    y"  be vector of 2nd derivs of approx function
If the points are numbered 0,1,2,...,n+1 then y" satisfies
(R W Hamming, Numerical Methods for Engineers and Scientists,
2nd Ed, p349ff)
	h\e9i\e8y"\b\e9i-1\e9\e8\e8+2(h\e9i\e8+h\e9i+1\e8)y"\b\e9i\e8+h\e9i+1\e8y"\b\e9i+1\e8
	
	= 6[(y\e9i+1\e8-y\e9i\e8)/h\e9i+1\e8-(y\e9i\e8-y\e9i-1\e8)/h\e9i\e8]   i=1,2,...,n

where y"\b\e90\e8 = y"\b\e9n+1\e8 = 0
This is a symmetric tridiagonal system of the form

	| a\e91\e8 h\e92\e8               |  |y"\b\e91\e8|      |b\e91\e8|
	| h\e92\e8 a\e92\e8 h\e93\e8            |  |y"\b\e92\e8|      |b\e92\e8|
	|    h\e93\e8 a\e93\e8 h\e94\e8         |  |y"\b\e93\e8|  =   |b\e93\e8|
	|         .           |  | .|      | .|
	|            .        |  | .|      | .|
It can be triangularized into
	| d\e91\e8 h\e92\e8               |  |y"\b\e91\e8|      |r\e91\e8|
	|    d\e92\e8 h\e93\e8            |  |y"\b\e92\e8|      |r\e92\e8|
	|       d\e93\e8 h\e94\e8         |  |y"\b\e93\e8|  =   |r\e93\e8|
	|          .          |  | .|      | .|
	|             .       |  | .|      | .|
where
	d\e91\e8 = a\e91\e8

	r\e90\e8 = 0

	d\e9i\e8 = a\e9i\e8 - h\e9i\e8\b\e82\e9/d\e9i-1\e8	1<i<\b_n

	r\e9i\e8 = b\e9i\e8 - h\e9i\e8r\e9i-1\e8/d\e9i-1\ei\e8	1<\b_i<\b_n

the back solution is
	y"\b\e9n\e8 = r\e9n\e8/d\e9n\e8

	y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y"\b\e9i+1\e8)/d\e9i\e8	1<\b_i<n

superficially, d\e9i\e8 and r\e9i\e8 don't have to be stored for they can be
recalculated backward by the formulas

	d\e9i-1\e8 = h\e9i\e8\b\e82\e9/(a\e9i\e8-d\e9i\e8)	1<i<\b_n

	r\e9i-1\e8 = (b\e9i\e8-r\e9i\e8)d\e9i-1\e8/h\e9i\e8	1<i<\b_n

unhappily it turns out that the recursion forward for d
is quite strongly geometrically convergent--and is wildly
unstable going backward.
There's similar trouble with r, so the intermediate
results must be kept.

Note that n-1 in the program below plays the role of n+1 in the theory

Other boundary conditions\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b_________________________

The boundary conditions are easily generalized to handle

	y\e90\e8\b" = ky\e91\e8\b", y\e9n+1\e8\b\b\b"   = ky\e9n\e8\b"

for some constant k.  The above analysis was for k = 0;
k = 1 fits parabolas perfectly as well as stright lines;
k = 1/2 has been recommended as somehow pleasant.

All that is necessary is to add h\e91\e8 to a\e91\e8 and h\e9n+1\e8 to a\e9n\e8.


Periodic case\b\b\b\b\b\b\b\b\b\b\b\b\b_____________

To do this, add 1 more row and column thus

	| a\e91\e8 h\e92\e8            h\e91\e8 |  |y\e91\e8\b"|     |b\e91\e8|
	| h\e92\e8 a\e92\e8 h\e93\e8            |  |y\e92\e8\b"|     |b\e92\e8|
	|    h\e93\e8 a\e94\e8 h\e94\e8         |  |y\e93\e8\b"|     |b\e93\e8|
	|                     |  | .|  =  | .|
	|             .       |  | .|     | .|
	| h\e91\e8            h\e90\e8 a\e90\e8 |  | .|     | .|

where h\e90\e8=\b_ h\e9n+1\e8

The same diagonalization procedure works, except for
the effect of the 2 corner elements.  Let s\e9i\e8 be the part
of the last element in the i\e8th\e9 "diagonalized" row that
arises from the extra top corner element.

		s\e91\e8 = h\e91\e8

		s\e9i\e8 = -s\e9i-1\e8h\e9i\e8/d\e9i-1\e8	2<\b_i<\b_n+1

After "diagonalizing", the lower corner element remains.
Call t\e9i\e8 the bottom element that appears in the i\e8th\e9 colomn
as the bottom element to its left is eliminated

		t\e91\e8 = h\e91\e8

		t\e9i\e8 = -t\e9i-1\e8h\e9i\e8/d\e9i-1\e8

Evidently t\e9i\e8 = s\e9i\e8.
Elimination along the bottom row
introduces further corrections to the bottom right element
and to the last element of the right hand side.
Call these corrections u and v.

	u\e91\e8 = v\e91\e8 = 0

	u\e9i\e8 = u\e9i-1\e8-s\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8

	v\e9i\e8 = v\e9i-1\e8-r\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8	2<\b_i<\b_n+1

The back solution is now obtained as follows

	y"\b\e9n+1\e8 = (r\e9n+1\e8+v\e9n+1\e8)/(d\e9n+1\e8+s\e9n+1\e8+t\e9n+1\e8+u\e9n+1\e8)

	y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8*y\e9i+1\e8-s\e9i\e8*y\e9n+1\e8)/d\e9i\e8	1<\b_i<\b_n

Interpolation in the interval x\e9i\e8<\b_x<\b_x\e9i+1\e8 is by the formula

	y = y\e9i\e8x\e9+\e8 + y\e9i+1\e8x\e9-\e8 -(h\e82\e9\b\e9i+1\e8/6)[y"\b\e9i\e8(x\e9+\e8-x\e9+\e8\e8\b3\e9)+y"\b\e9i+1\e8(x\e9-\e8-x\e9-\e8\b\e83\e9)]
where
	x\e9+\e8 = x\e9i+1\e8-x

	x\e9-\e8 = x-x\e9i\e8
*/

float
rhs(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));
}

spline(){
	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 = (end*v+r[i]-hi1*D2yi1-s*D2yn1)/
				(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;
			printf("%f ",x.val[i]+j*h);
			printf("%f\n",yy);
			}
		}
	return(1);
	}
readin() {
	for(n=0;n<NP;n++){
		if(auta) x.val[n] = n*dx+x.lb;
		else if(!getfloat(&x.val[n])) break;
		if(!getfloat(&y.val[n])) break; } }

getfloat(p)
	float *p;{
	char buf[30];
	register c;
	int i;
	extern double atof();
	for(;;){
		c = getchar();
		if (c==EOF) {
			*buf = '\0';
			return(0);
		}
		*buf = c;
		switch(*buf){
			case ' ':
			case '\t':
			case '\n':
				continue;}
		break;}
	for(i=1;i<30;i++){
		c = getchar();
		if (c==EOF) {
			buf[i] = '\0';
			break;
		}
		buf[i] = c;
		if('0'<=c && c<='9') continue;
		switch(c) {
			case '.':
			case '+':
			case '-':
			case 'E':
			case 'e':
				continue;}
		break; }
	buf[i] = ' ';
	*p = atof(buf);
	return(1); }

getlim(p)
	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]; }
	}


main(argc,argv)
	char *argv[];{
	extern char *malloc();
	int i;
	x.lbf = x.ubf = y.lbf = y.ubf = 0;
	x.lb = INF;
	x.ub = -INF;
	y.lb = INF;
	y.ub = -INF;
	while(--argc > 0) {
		argv++;
again:		switch(argv[0][0]) {
		case '-':
			argv[0]++;
			goto again;
		case 'a':
			auta = 1;
			numb(&dx,&argc,&argv);
			break;
		case 'k':
			numb(&konst,&argc,&argv);
			break;
		case 'n':
			numb(&ni,&argc,&argv);
			break;
		case 'p':
			periodic = 1;
			break;
		case 'x':
			if(!numb(&x.lb,&argc,&argv)) break;
			x.lbf = 1;
			if(!numb(&x.ub,&argc,&argv)) break;
			x.ubf = 1;
			break;
		default:
			fprintf(stderr, "Bad agrument\n");
			exit(1);
		}
	}
	if(auta&&!x.lbf) x.lb = 0;
	readin();
	getlim(&x);
	getlim(&y);
	i = (n+1)*sizeof(dx);
	diag = (float *)malloc((unsigned)i);
	r = (float *)malloc((unsigned)i);
	if(r==NULL||!spline()) for(i=0;i<n;i++){
		printf("%f ",x.val[i]);
		printf("%f\n",y.val[i]); }
}
numb(np,argcp,argvp)
	int *argcp;
	float *np;
	char ***argvp;{
	double atof();
	char c;
	if(*argcp<=1) return(0);
	c = (*argvp)[1][0];
	if(!('0'<=c&&c<='9' || c=='-' || c== '.' )) return(0);
	*np = atof((*argvp)[1]);
	(*argcp)--;
	(*argvp)++; 
	return(1); }
Commit	Line	Data
762ee6b6	1	static char *sccsid = "@(#)spline.c 4.2 (Berkeley) %G%";
365384c2	2	#include <stdio.h>
762ee6b6	3	#include <math.h>
365384c2 BJ	4
365384c2 BJ	5	#define NP 1000
762ee6b6	6	#define INF HUGE
365384c2 BJ	7
	8	struct proj { int lbf,ubf; float a,b,lb,ub,quant,mult,val[NP]; } x,y;
	9	float diag, r;
	10	float dx = 1.;
	11	float ni = 100.;
	12	int n;
	13	int auta;
	14	int periodic;
	15	float konst = 0.0;
	16	float zero = 0.;
	17
	18	/* Spline fit technique
	19	let x,y be vectors of abscissas and ordinates
	20	h be vector of differences h\e9i\e8=x\e9i\e8-x\e9i-1\e\e9\e8\e8
	21	y" be vector of 2nd derivs of approx function
	22	If the points are numbered 0,1,2,...,n+1 then y" satisfies
	23	(R W Hamming, Numerical Methods for Engineers and Scientists,
	24	2nd Ed, p349ff)
	25	h\e9i\e8y"\b\e9i-1\e9\e8\e8+2(h\e9i\e8+h\e9i+1\e8)y"\b\e9i\e8+h\e9i+1\e8y"\b\e9i+1\e8
	26
	27	= 6[(y\e9i+1\e8-y\e9i\e8)/h\e9i+1\e8-(y\e9i\e8-y\e9i-1\e8)/h\e9i\e8] i=1,2,...,n
	28
	29	where y"\b\e90\e8 = y"\b\e9n+1\e8 = 0
	30	This is a symmetric tridiagonal system of the form
	31
	32	\| a\e91\e8 h\e92\e8 \| \|y"\b\e91\e8\| \|b\e91\e8\|
	33	\| h\e92\e8 a\e92\e8 h\e93\e8 \| \|y"\b\e92\e8\| \|b\e92\e8\|
	34	\| h\e93\e8 a\e93\e8 h\e94\e8 \| \|y"\b\e93\e8\| = \|b\e93\e8\|
	35	\| . \| \| .\| \| .\|
	36	\| . \| \| .\| \| .\|
	37	It can be triangularized into
	38	\| d\e91\e8 h\e92\e8 \| \|y"\b\e91\e8\| \|r\e91\e8\|
	39	\| d\e92\e8 h\e93\e8 \| \|y"\b\e92\e8\| \|r\e92\e8\|
	40	\| d\e93\e8 h\e94\e8 \| \|y"\b\e93\e8\| = \|r\e93\e8\|
	41	\| . \| \| .\| \| .\|
	42	\| . \| \| .\| \| .\|
	43	where
	44	d\e91\e8 = a\e91\e8
	45
	46	r\e90\e8 = 0
	47
	48	d\e9i\e8 = a\e9i\e8 - h\e9i\e8\b\e82\e9/d\e9i-1\e8 1<i<\b_n
	49
	50	r\e9i\e8 = b\e9i\e8 - h\e9i\e8r\e9i-1\e8/d\e9i-1\ei\e8 1<\b_i<\b_n
	51
	52	the back solution is
	53	y"\b\e9n\e8 = r\e9n\e8/d\e9n\e8
	54
	55	y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y"\b\e9i+1\e8)/d\e9i\e8 1<\b_i<n
	56
	57	superficially, d\e9i\e8 and r\e9i\e8 don't have to be stored for they can be
	58	recalculated backward by the formulas
	59
	60	d\e9i-1\e8 = h\e9i\e8\b\e82\e9/(a\e9i\e8-d\e9i\e8) 1<i<\b_n
	61
	62	r\e9i-1\e8 = (b\e9i\e8-r\e9i\e8)d\e9i-1\e8/h\e9i\e8 1<i<\b_n
	63
	64	unhappily it turns out that the recursion forward for d
	65	is quite strongly geometrically convergent--and is wildly
	66	unstable going backward.
	67	There's similar trouble with r, so the intermediate
	68	results must be kept.
	69
	70	Note that n-1 in the program below plays the role of n+1 in the theory
71
72	Other boundary conditions\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b_________________________
73
74	The boundary conditions are easily generalized to handle
75
76	y\e90\e8\b" = ky\e91\e8\b", y\e9n+1\e8\b\b\b" = ky\e9n\e8\b"
77
78	for some constant k. The above analysis was for k = 0;
79	k = 1 fits parabolas perfectly as well as stright lines;
80	k = 1/2 has been recommended as somehow pleasant.
81
82	All that is necessary is to add h\e91\e8 to a\e91\e8 and h\e9n+1\e8 to a\e9n\e8.
83
84
85	Periodic case\b\b\b\b\b\b\b\b\b\b\b\b\b_____________
86
87	To do this, add 1 more row and column thus
88
89	\| a\e91\e8 h\e92\e8 h\e91\e8 \| \|y\e91\e8\b"\| \|b\e91\e8\|
90	\| h\e92\e8 a\e92\e8 h\e93\e8 \| \|y\e92\e8\b"\| \|b\e92\e8\|
91	\| h\e93\e8 a\e94\e8 h\e94\e8 \| \|y\e93\e8\b"\| \|b\e93\e8\|
92	\| \| \| .\| = \| .\|
93	\| . \| \| .\| \| .\|
94	\| h\e91\e8 h\e90\e8 a\e90\e8 \| \| .\| \| .\|
95
96	where h\e90\e8=\b_ h\e9n+1\e8
97
98	The same diagonalization procedure works, except for
99	the effect of the 2 corner elements. Let s\e9i\e8 be the part
100	of the last element in the i\e8th\e9 "diagonalized" row that
101	arises from the extra top corner element.
102
103	s\e91\e8 = h\e91\e8
104
105	s\e9i\e8 = -s\e9i-1\e8h\e9i\e8/d\e9i-1\e8 2<\b_i<\b_n+1
106
107	After "diagonalizing", the lower corner element remains.
108	Call t\e9i\e8 the bottom element that appears in the i\e8th\e9 colomn
109	as the bottom element to its left is eliminated
110
111	t\e91\e8 = h\e91\e8
112
113	t\e9i\e8 = -t\e9i-1\e8h\e9i\e8/d\e9i-1\e8
114
115	Evidently t\e9i\e8 = s\e9i\e8.
116	Elimination along the bottom row
117	introduces further corrections to the bottom right element
118	and to the last element of the right hand side.
119	Call these corrections u and v.
120
121	u\e91\e8 = v\e91\e8 = 0
122
123	u\e9i\e8 = u\e9i-1\e8-s\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8
124
125	v\e9i\e8 = v\e9i-1\e8-r\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8 2<\b_i<\b_n+1
126
127	The back solution is now obtained as follows
128
129	y"\b\e9n+1\e8 = (r\e9n+1\e8+v\e9n+1\e8)/(d\e9n+1\e8+s\e9n+1\e8+t\e9n+1\e8+u\e9n+1\e8)
130
131	y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y\e9i+1\e8-s\e9i\e8y\e9n+1\e8)/d\e9i\e8 1<\b_i<\b_n
132
133	Interpolation in the interval x\e9i\e8<\b_x<\b_x\e9i+1\e8 is by the formula
134
135	y = y\e9i\e8x\e9+\e8 + y\e9i+1\e8x\e9-\e8 -(h\e82\e9\b\e9i+1\e8/6)[y"\b\e9i\e8(x\e9+\e8-x\e9+\e8\e8\b3\e9)+y"\b\e9i+1\e8(x\e9-\e8-x\e9-\e8\b\e83\e9)]
136	where
137	x\e9+\e8 = x\e9i+1\e8-x
138
139	x\e9-\e8 = x-x\e9i\e8
140	*/
141
142	float
143	rhs(i){
144	int i_;
145	double zz;
146	i_ = i==n-1?0:i;
147	zz = (y.val[i]-y.val[i-1])/(x.val[i]-x.val[i-1]);
148	return(6*((y.val[i_+1]-y.val[i_])/(x.val[i+1]-x.val[i]) - zz));
149	}
150
151	spline(){
152	float d,s,u,v,hi,hi1;
153	float h;
154	float D2yi,D2yi1,D2yn1,x0,x1,yy,a;
155	int end;
156	float corr;
157	int i,j,m;
158	if(n<3) return(0);
159	if(periodic) konst = 0;
160	d = 1;
161	r[0] = 0;
162	s = periodic?-1:0;
163	for(i=0;++i<n-!periodic;){ /* triangularize */
164	hi = x.val[i]-x.val[i-1];
165	hi1 = i==n-1?x.val[1]-x.val[0]:
166	x.val[i+1]-x.val[i];
167	if(hi1*hi<=0) return(0);
168	u = i==1?zero:u-s*s/d;
169	v = i==1?zero:v-s*r[i-1]/d;
170	r[i] = rhs(i)-hi*r[i-1]/d;
171	s = -hi*s/d;
172	a = 2*(hi+hi1);
173	if(i==1) a += konst*hi;
174	if(i==n-2) a += konst*hi1;
175	diag[i] = d = i==1? a:
176	a - hi*hi/d;
177	}
178	D2yi = D2yn1 = 0;
179	for(i=n-!periodic;--i>=0;){ /* back substitute */
180	end = i==n-1;
181	hi1 = end?x.val[1]-x.val[0]:
182	x.val[i+1]-x.val[i];
183	D2yi1 = D2yi;
184	if(i>0){
185	hi = x.val[i]-x.val[i-1];
186	corr = end?2*s+u:zero;
187	D2yi = (endv+r[i]-hi1D2yi1-s*D2yn1)/
188	(diag[i]+corr);
189	if(end) D2yn1 = D2yi;
190	if(i>1){
191	a = 2*(hi+hi1);
192	if(i==1) a += konst*hi;
193	if(i==n-2) a += konst*hi1;
194	d = diag[i-1];
195	s = -s*d/hi;
196	}}
197	else D2yi = D2yn1;
198	if(!periodic) {
199	if(i==0) D2yi = konst*D2yi1;
200	if(i==n-2) D2yi1 = konst*D2yi;
201	}
202	if(end) continue;
203	m = hi1>0?ni:-ni;
204	m = 1.001mhi1/(x.ub-x.lb);
205	if(m<=0) m = 1;
206	h = hi1/m;
207	for(j=m;j>0\|\|i==0&&j==0;j--){ /* interpolate */
208	x0 = (m-j)*h/hi1;
209	x1 = j*h/hi1;
210	yy = D2yi(x0-x0x0x0)+D2yi1(x1-x1x1x1);
211	yy = y.val[i]x0+y.val[i+1]x1 -hi1hi1yy/6;
212	printf("%f ",x.val[i]+j*h);
213	printf("%f\n",yy);
214	}
215	}
216	return(1);
217	}
218	readin() {
219	for(n=0;n<NP;n++){
220	if(auta) x.val[n] = n*dx+x.lb;
221	else if(!getfloat(&x.val[n])) break;
222	if(!getfloat(&y.val[n])) break; } }
223
224	getfloat(p)
225	float *p;{
226	char buf[30];
227	register c;
228	int i;
229	extern double atof();
230	for(;;){
231	c = getchar();
232	if (c==EOF) {
233	*buf = '\0';
234	return(0);
235	}
236	*buf = c;
237	switch(*buf){
238	case ' ':
239	case '\t':
240	case '\n':
241	continue;}
242	break;}
243	for(i=1;i<30;i++){
244	c = getchar();
245	if (c==EOF) {
246	buf[i] = '\0';
247	break;
248	}
249	buf[i] = c;
250	if('0'<=c && c<='9') continue;
251	switch(c) {
252	case '.':
253	case '+':
254	case '-':
255	case 'E':
256	case 'e':
257	continue;}
258	break; }
259	buf[i] = ' ';
260	*p = atof(buf);
261	return(1); }
262
263	getlim(p)
264	struct proj *p; {
265	int i;
266	for(i=0;i<n;i++) {
267	if(!p->lbf && p->lb>(p->val[i])) p->lb = p->val[i];
268	if(!p->ubf && p->ub<(p->val[i])) p->ub = p->val[i]; }
269	}
270
271
272	main(argc,argv)
273	char *argv[];{
274	extern char *malloc();
275	int i;
276	x.lbf = x.ubf = y.lbf = y.ubf = 0;
277	x.lb = INF;
278	x.ub = -INF;
279	y.lb = INF;
280	y.ub = -INF;
281	while(--argc > 0) {
282	argv++;
283	again: switch(argv[0][0]) {
284	case '-':
285	argv[0]++;
286	goto again;
287	case 'a':
288	auta = 1;
289	numb(&dx,&argc,&argv);
290	break;
291	case 'k':
292	numb(&konst,&argc,&argv);
293	break;
294	case 'n':
295	numb(&ni,&argc,&argv);
296	break;
297	case 'p':
298	periodic = 1;
299	break;
300	case 'x':
301	if(!numb(&x.lb,&argc,&argv)) break;
302	x.lbf = 1;
303	if(!numb(&x.ub,&argc,&argv)) break;
304	x.ubf = 1;
305	break;
306	default:
307	fprintf(stderr, "Bad agrument\n");
308	exit(1);
309	}
310	}
311	if(auta&&!x.lbf) x.lb = 0;
312	readin();
313	getlim(&x);
314	getlim(&y);
315	i = (n+1)*sizeof(dx);
316	diag = (float *)malloc((unsigned)i);
317	r = (float *)malloc((unsigned)i);
318	if(r==NULL\|\|!spline()) for(i=0;i<n;i++){
319	printf("%f ",x.val[i]);
320	printf("%f\n",y.val[i]); }
321	}
322	numb(np,argcp,argvp)
323	int *argcp;
324	float *np;
325	char ***argvp;{
326	double atof();
327	char c;
328	if(*argcp<=1) return(0);
329	c = (*argvp)[1][0];
330	if(!('0'<=c&&c<='9' \|\| c=='-' \|\| c== '.' )) return(0);
331	np = atof((argvp)[1]);
332	(*argcp)--;
333	(*argvp)++;
334	return(1); }
335