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

static char *sccsid = "@(#)spline.c	4.1 (Berkeley) 10/1/80";
#include <stdio.h>

#define NP 1000
#define INF 1.e37

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