13e12c5d1SDavid du Colombier /*
23e12c5d1SDavid du Colombier Produce spline (uniform knots, second order)
33e12c5d1SDavid du Colombier from guiding points
43e12c5d1SDavid du Colombier */
53e12c5d1SDavid du Colombier #include "mplot.h"
splin(int mode,int num[],double * ff[])63e12c5d1SDavid du Colombier void splin(int mode, int num[], double *ff[]){
73e12c5d1SDavid du Colombier int i, *np, n;
83e12c5d1SDavid du Colombier double xa, ya, xc, yc, *xp, *yp, *xp0, *yp0, *xpe, *ype;
93e12c5d1SDavid du Colombier double **fp;
103e12c5d1SDavid du Colombier np = num;
113e12c5d1SDavid du Colombier fp = ff;
123e12c5d1SDavid du Colombier while((n = *np++)){
133e12c5d1SDavid du Colombier xp = *fp++;
143e12c5d1SDavid du Colombier yp = xp + 1;
153e12c5d1SDavid du Colombier xp0 = xp;
163e12c5d1SDavid du Colombier yp0 = yp;
173e12c5d1SDavid du Colombier xpe = xp0 + 2 * (n - 1);
183e12c5d1SDavid du Colombier ype = yp0 + 2 * (n - 1);
193e12c5d1SDavid du Colombier if (n < 3) {
20*7dd7cddfSDavid du Colombier plotline(*xp, *yp, *(xp + 2), *(yp + 2));
213e12c5d1SDavid du Colombier continue;
223e12c5d1SDavid du Colombier }
233e12c5d1SDavid du Colombier if (mode == 4) { /*closed curve*/
243e12c5d1SDavid du Colombier xa = 0.5 * (*xpe + *(xpe - 2));
253e12c5d1SDavid du Colombier xc = 0.5 * (*xpe + *xp0);
263e12c5d1SDavid du Colombier ya = 0.5 * (*ype + *(ype - 2));
273e12c5d1SDavid du Colombier yc = 0.5 * (*ype + *yp0);
283e12c5d1SDavid du Colombier parabola(xa, ya, xc, yc, *xpe, *ype);
293e12c5d1SDavid du Colombier xa = 0.5 * (*xpe + *xp0);
303e12c5d1SDavid du Colombier xc = 0.5 * (*(xp0 + 2) + *xp0);
313e12c5d1SDavid du Colombier ya = 0.5 * (*ype + *yp0);
323e12c5d1SDavid du Colombier yc = 0.5 * (*(yp0 + 2) + *yp0);
333e12c5d1SDavid du Colombier parabola(xa, ya, xc, yc, *xp0, *yp0);
343e12c5d1SDavid du Colombier }
353e12c5d1SDavid du Colombier else { /*open curve with multiple endpoints*/
363e12c5d1SDavid du Colombier if (mode % 2) /*odd mode makes first point double*/
37*7dd7cddfSDavid du Colombier plotline(*xp0,*yp0,0.5*(*xp0+*(xp0+2)),0.5*(*yp0+*(yp0+2)));
383e12c5d1SDavid du Colombier }
393e12c5d1SDavid du Colombier xp += 2;
403e12c5d1SDavid du Colombier yp += 2;
413e12c5d1SDavid du Colombier for (i = 1; i < (n - 1); i++, xp += 2, yp += 2) {
423e12c5d1SDavid du Colombier xa = 0.5 * (*(xp - 2) + *xp);
433e12c5d1SDavid du Colombier xc = 0.5 * ( *xp + *(xp + 2));
443e12c5d1SDavid du Colombier ya = 0.5 * (*(yp - 2) + *yp);
453e12c5d1SDavid du Colombier yc = 0.5 * ( *yp + *(yp + 2));
463e12c5d1SDavid du Colombier parabola(xa, ya, xc, yc, *xp, *yp);
473e12c5d1SDavid du Colombier }
483e12c5d1SDavid du Colombier if(mode >= 2 && mode != 4)
49*7dd7cddfSDavid du Colombier plotline(0.5*(*(xpe-2)+*xpe),0.5*(*(ype-2)+*ype),*xpe,*ype);
503e12c5d1SDavid du Colombier }
513e12c5d1SDavid du Colombier }
52