function f(u,p){return (u**3)*(10-p+(2*p-15)*u+(6-p)*u**2); } 
function g(u,q,p){ u=scaleu(u); return q*u+2*q*u**2 + (10-12*q-p)*u**3 + (2*p+14*q-15)*u**4 + (6-5*q-p)*u**5; }
function h(u,q,p){ u=scaleu(u); return q*u + 2*q*u**2 - 2*q*u**4 - q*u**5; }

function a0(t, tb){ tb=int(t)-1; return t > bigtp[tb]  ?  0 : h( (t-bigtp[tb])/(tb-bigtp[tb]) , q[tb], p[-1] ); } 
function a1(t, tb){ tb=int(t)-1; return t > bigtp[tb+1] ? 0 : g( (t-bigtp[tb+1])/(tb+1-bigtp[tb+1]) ,q[tb+1], p[0] ); }
function a2(t, tb){ tb=int(t)-1; return t < bigtm[tb+2] ? 0 : g( (t-bigtm[tb+2])/(tb+2-bigtm[tb+2]) ,q[tb+2], p[1] ); }
function a3(t, tb){ tb=int(t)-1; return t < bigtm[tb+3] ? 0 : h( (t-bigtm[tb+3])/(tb+3-bigtm[tb+3]) ,q[tb+3], p[2] ); }

function scaleu(u)
 {
   #return u;  # required if s==0
   if( u<0.5 )
      u=2*u-1;
   else 
      u=2*(u-0.5);
  return u;
 }

# this is test code from X-splines: a spline model designed for the end user.
# this version implements the third stage (section 5), with the control parameter s[] -1 to 1
# except that s[] must be -1, ie, must always go through the control points. Cannot
# figure out the maths to do what the author did, make -1 to 1 possible on the fly.
# Rather than continue debugging that, modified the test code to read the coards
# from stdin and output interpolated coards to stdout

{
  x[n]=$1;  y[n]=$2;  z[n]=$3;  s[n]=-1;
  n++;
}

BEGIN {
  n=1;
 }

END { 

   s[n]=-1;

# def=-1;
# x[1]=1;  y[1]=5;   z[1]=5;   s[1]=def;   # ==0 and curve goes through that point,  ==1 and curves goes past it
# x[2]=3;  y[2]=10;  z[2]=5;   s[2]=def;   
# x[3]=10; y[3]=10;  z[3]=1;   s[3]=def;
# x[4]=13; y[4]=5;   z[4]=1;   s[4]=def;
# x[5]=10; y[5]=0;   z[5]=10;  s[5]=-1;
# x[6]=20; y[6]=0;   z[6]=1;   s[6]=def;
# x[7]=15; y[7]=5;   z[7]=5;   s[7]=def;
#                              s[8]=def;

   delta=1;   # t_k-1 - t_k , which is assumed to be a constant for all k 

   for(t=1;t<=n+1;t++)
    {
      if( s[t]<0 )
       { q[t]=-s[t]/2;  s[t]=-s[t]; }   # she never says what to do with s[] when operating with q[]
      else 
         q[t]=s[t]/2;   # nor does she say what to do with q[] when s[]>0  . xfig code says to use f()
    }

   for(t=1;t<=n;t++)
    {
      bigtp[t]=t+1+s[t+1]*delta;
      bigtm[t]=t-1-s[t-1]*delta;
    }

   #for(u=0;u<=1;u+=0.01)
   #   printf "%f\t%f\t%f\t%f\n",u,f(u,8), g(2*(u-0.5),0.5,8), h(2*u-1,0.5,8);   
   #
   #exit; 

   #n=8;

   for(t=1;t<=n-1;t+=0.01)
    {
      bt=int(t-1);
  
      k=bt;


      p[-1]=(2/(delta**2)) * (k-bigtp[k])**2;
      p[0] =(2/(delta**2)) * (k+1-bigtp[k+1])**2;
      p[1] =(2/(delta**2)) * (k+2-bigtm[k+2])**2;
      p[2] =(2/(delta**2)) * (k+3-bigtm[k+3])**2;

      cx[t]= ( a0(t)*x[bt] + a1(t)*x[bt+1] + a2(t)*x[bt+2] + a3(t)*x[bt+3] ) / (a0(t) + a1(t) + a2(t) + a3(t) );
      cy[t]= ( a0(t)*y[bt] + a1(t)*y[bt+1] + a2(t)*y[bt+2] + a3(t)*y[bt+3] ) / (a0(t) + a1(t) + a2(t) + a3(t) );
      cz[t]= ( a0(t)*z[bt] + a1(t)*z[bt+1] + a2(t)*z[bt+2] + a3(t)*z[bt+3] ) / (a0(t) + a1(t) + a2(t) + a3(t) );

      printf "%f %f %f %f\n",t, cx[t],cy[t],cz[t];
    }
 }  # END