1*3e12c5d1SDavid du Colombier /* 2*3e12c5d1SDavid du Colombier Produce spline (uniform knots, second order) 3*3e12c5d1SDavid du Colombier from guiding points 4*3e12c5d1SDavid du Colombier */ 5*3e12c5d1SDavid du Colombier #include "mplot.h" 6*3e12c5d1SDavid du Colombier void splin(int mode, int num[], double *ff[]){ 7*3e12c5d1SDavid du Colombier int i, *np, n; 8*3e12c5d1SDavid du Colombier double xa, ya, xc, yc, *xp, *yp, *xp0, *yp0, *xpe, *ype; 9*3e12c5d1SDavid du Colombier double **fp; 10*3e12c5d1SDavid du Colombier np = num; 11*3e12c5d1SDavid du Colombier fp = ff; 12*3e12c5d1SDavid du Colombier while((n = *np++)){ 13*3e12c5d1SDavid du Colombier xp = *fp++; 14*3e12c5d1SDavid du Colombier yp = xp + 1; 15*3e12c5d1SDavid du Colombier xp0 = xp; 16*3e12c5d1SDavid du Colombier yp0 = yp; 17*3e12c5d1SDavid du Colombier xpe = xp0 + 2 * (n - 1); 18*3e12c5d1SDavid du Colombier ype = yp0 + 2 * (n - 1); 19*3e12c5d1SDavid du Colombier if (n < 3) { 20*3e12c5d1SDavid du Colombier line(*xp, *yp, *(xp + 2), *(yp + 2)); 21*3e12c5d1SDavid du Colombier continue; 22*3e12c5d1SDavid du Colombier } 23*3e12c5d1SDavid du Colombier if (mode == 4) { /*closed curve*/ 24*3e12c5d1SDavid du Colombier xa = 0.5 * (*xpe + *(xpe - 2)); 25*3e12c5d1SDavid du Colombier xc = 0.5 * (*xpe + *xp0); 26*3e12c5d1SDavid du Colombier ya = 0.5 * (*ype + *(ype - 2)); 27*3e12c5d1SDavid du Colombier yc = 0.5 * (*ype + *yp0); 28*3e12c5d1SDavid du Colombier parabola(xa, ya, xc, yc, *xpe, *ype); 29*3e12c5d1SDavid du Colombier xa = 0.5 * (*xpe + *xp0); 30*3e12c5d1SDavid du Colombier xc = 0.5 * (*(xp0 + 2) + *xp0); 31*3e12c5d1SDavid du Colombier ya = 0.5 * (*ype + *yp0); 32*3e12c5d1SDavid du Colombier yc = 0.5 * (*(yp0 + 2) + *yp0); 33*3e12c5d1SDavid du Colombier parabola(xa, ya, xc, yc, *xp0, *yp0); 34*3e12c5d1SDavid du Colombier } 35*3e12c5d1SDavid du Colombier else { /*open curve with multiple endpoints*/ 36*3e12c5d1SDavid du Colombier if (mode % 2) /*odd mode makes first point double*/ 37*3e12c5d1SDavid du Colombier line(*xp0,*yp0,0.5*(*xp0+*(xp0+2)),0.5*(*yp0+*(yp0+2))); 38*3e12c5d1SDavid du Colombier } 39*3e12c5d1SDavid du Colombier xp += 2; 40*3e12c5d1SDavid du Colombier yp += 2; 41*3e12c5d1SDavid du Colombier for (i = 1; i < (n - 1); i++, xp += 2, yp += 2) { 42*3e12c5d1SDavid du Colombier xa = 0.5 * (*(xp - 2) + *xp); 43*3e12c5d1SDavid du Colombier xc = 0.5 * ( *xp + *(xp + 2)); 44*3e12c5d1SDavid du Colombier ya = 0.5 * (*(yp - 2) + *yp); 45*3e12c5d1SDavid du Colombier yc = 0.5 * ( *yp + *(yp + 2)); 46*3e12c5d1SDavid du Colombier parabola(xa, ya, xc, yc, *xp, *yp); 47*3e12c5d1SDavid du Colombier } 48*3e12c5d1SDavid du Colombier if(mode >= 2 && mode != 4) 49*3e12c5d1SDavid du Colombier line(0.5*(*(xpe-2)+*xpe),0.5*(*(ype-2)+*ype),*xpe,*ype); 50*3e12c5d1SDavid du Colombier } 51*3e12c5d1SDavid du Colombier } 52