1  
   2  float cubic_interpolate( float y0, float y1, float y2, float y3, float mu ) {
   3  
   4     float a0, a1, a2, a3, mu2;
   5  
   6     mu2 = mu*mu;
   7     a0 = y3 - y2 - y0 + y1; //p
   8     a1 = y0 - y1 - a0;
   9     a2 = y2 - y0;
  10     a3 = y1;
  11  
  12     return ( a0*mu*mu2 + a1*mu2 + a2*mu + a3 );
  13  }