// cubic interpolation function for equidistant 4 points
// 0 <= mu <=1


float cubic_interpolate( float y0, float y1, float y2, float y3, float mu ) {

   float a0, a1, a2, a3, mu2;

   mu2 = mu*mu;
   a0 = y3 - y2 - y0 + y1; //p
   a1 = y0 - y1 - a0;
   a2 = y2 - y0;
   a3 = y1;

   return ( a0*mu*mu2 + a1*mu2 + a2*mu + a3 );
}