#ifndef lint static char *sccsid = "@(#)spline.c 4.3 (Berkeley) 9/21/85"; #endif #include #include #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 hi=xi-xi-1 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) hiy"i-1+2(hi+hi+1)y"i+hi+1y"i+1 = 6[(yi+1-yi)/hi+1-(yi-yi-1)/hi] i=1,2,...,n where y"0 = y"n+1 = 0 This is a symmetric tridiagonal system of the form | a1 h2 | |y"1| |b1| | h2 a2 h3 | |y"2| |b2| | h3 a3 h4 | |y"3| = |b3| | . | | .| | .| | . | | .| | .| It can be triangularized into | d1 h2 | |y"1| |r1| | d2 h3 | |y"2| |r2| | d3 h4 | |y"3| = |r3| | . | | .| | .| | . | | .| | .| where d1 = a1 r0 = 0 di = ai - hi2/di-1 1=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;nlbf && 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