double interpol1d(double x,int n,double *xbuf,double *buf);
double interpolAngular(double x,int n,double *xbuf,double *buf);
double interpol2d(double x,double y,int nx,int ny,double *xbuf,double *ybuf,double *buf);

double interpol1d(double x,int n,double *xbuf,double *buf)
/*----------------------------------------------------------------------
 * Parabolic blending interpolation of a function of one variable.
 *
 * Input parameters:
 *    x    abscissa
 *    n    number of data points of abscissa and ordinate
 *    xbuf array with pre-defined values of x, ascending
 *    buf  array with function values, where buf[i] is function of xbuf[i]
 *
 * Terms of use:
 *    You are not allowed to claim that this code is the work of anybody else
 *    than Henning Weede.
 *    The author is not liable for any damage caused by this code.
 *    The author has no duty to waist his time answering questions related
 *    to this code.
 *    The GNU Public License applies.
 *
 * Author: Henning Weede
 *----------------------------------------------------------------------
 */{
   int ix;
   double vor,nach,tmp;

   /* check if ascending order (overflow protection) */
   tmp=fabs(xbuf[n-1]-xbuf[0])/10000.;
   for(ix=1; ix<n; ++ix)
      if(xbuf[ix]-xbuf[ix-1]<tmp)
         {
         fprintf(stderr,"interpol1d(): abscissa not in ascending order.\n");
         exit(1);
         }

   /* if less than 4 data points: other methods */
   switch(n)
      {
      case 0:
         fprintf(stderr,"interpol1d(): abscissa array is empty.\n");
         exit(1);
      case 1:
         return(buf[0]);
      case 2:
         tmp=(x-xbuf[0])/(xbuf[1]-xbuf[0]);
         return((1.-tmp)*buf[0]+tmp*buf[1]);
      case 3:
         tmp=buf[0]*(x-xbuf[1])*(x-xbuf[2])/((xbuf[0]-xbuf[1])*(xbuf[0]-xbuf[2]))
            +buf[1]*(x-xbuf[0])*(x-xbuf[2])/((xbuf[1]-xbuf[0])*(xbuf[1]-xbuf[2]))
            +buf[2]*(x-xbuf[0])*(x-xbuf[1])/((xbuf[2]-xbuf[0])*(xbuf[2]-xbuf[1]));
         return(tmp);
      }

   /* in which interval are we? */
   for(ix=1; ix<n-1; ++ix)
      if(xbuf[ix]>x) break;

   /* interpolate in pre-parabola, return only this if no post-parabola */
   if(ix>1)
      vor=buf[ix-2]*(x-xbuf[ix-1])*(x-xbuf[ix  ])/((xbuf[ix-2]-xbuf[ix-1])*(xbuf[ix-2]-xbuf[ix  ]))
         +buf[ix-1]*(x-xbuf[ix-2])*(x-xbuf[ix  ])/((xbuf[ix-1]-xbuf[ix-2])*(xbuf[ix-1]-xbuf[ix  ]))
         +buf[ix  ]*(x-xbuf[ix-2])*(x-xbuf[ix-1])/((xbuf[ix  ]-xbuf[ix-2])*(xbuf[ix  ]-xbuf[ix-1]));
   if(ix>=n-1) return(vor);

   /* interpolate in post-parabola, return only this if no pre-parabola */
   if(ix<n-1)
      nach=buf[ix-1]*(x-xbuf[ix  ])*(x-xbuf[ix+1])/((xbuf[ix-1]-xbuf[ix  ])*(xbuf[ix-1]-xbuf[ix+1]))
          +buf[ix  ]*(x-xbuf[ix-1])*(x-xbuf[ix+1])/((xbuf[ix  ]-xbuf[ix-1])*(xbuf[ix  ]-xbuf[ix+1]))
          +buf[ix+1]*(x-xbuf[ix-1])*(x-xbuf[ix  ])/((xbuf[ix+1]-xbuf[ix-1])*(xbuf[ix+1]-xbuf[ix  ]));
   if(ix<=1) return(nach);

   /* linearly blend post- and pre-parabola */
   tmp=(x-xbuf[ix-1])/(xbuf[ix]-xbuf[ix-1]);
   return(vor*(1.-tmp)+nach*tmp);
   } /* interpol1d() */

double interpolAngular(double x,int n,double *xbuf,double *buf)
/*----------------------------------------------------------------------
 * Parabolic blending interpolation of a function of an angle.
 *
 * Input parameters:
 *    x    angular abscissa [rad]
 *    n    number of data points of abscissa and ordinate
 *    xbuf array with pre-defined values of x [rad], ascending
 *    buf  array with function values, where buf[i] is function of xbuf[i]
 *
 * Terms of use:
 *    You are not allowed to claim that this code is the work of anybody else
 *    than Henning Weede.
 *    The author is not liable for any damage caused by this code.
 *    The author has no duty to waist his time answering questions related
 *    to this code.
 *    The GNU Public License applies.
 *
 * Author: Henning Weede
 *----------------------------------------------------------------------
 */{
   int i0,i1,i2,i3,ix;
   double vor,nach,tmp,x0,x1,x2,x3,y0,y1,y2,y3;

   for(ix=1; ix<n; ++ix)
      if(xbuf[ix]-xbuf[ix-1]<PI/18000.)
         {
         fprintf(stderr,"interpolAngular(): angular abscissa not in ascending order.\n");
         exit(1);
         }
   if(xbuf[0]-xbuf[n-1]+2*PI>PI/2.)
       {
       fprintf(stderr,"interpolAngular(): huge gap between beginning and end of angular abscissa.\n");
       exit(1);
       }

   while(x<xbuf[0])   x+=2*PI;
   while(x>xbuf[n-1]) x-=2*PI;

   for(i2=0; i2<n; ++i2)
      if(xbuf[i2]>x) break;
   if(i2==n) i2=0;
   i1=i2-1; if(i1<0) i1+=n;
   i0=i1-1; if(i0<0) i0+=n;
   i3=i2+1; if(i3>n-1) i3-=n;

   x0=xbuf[i0];
   x1=xbuf[i1];
   x2=xbuf[i2];
   x3=xbuf[i3];

   while(x1>x0+2*PI) x1-=2*PI;
   while(x1<x0) x1+=2*PI;
   while(x2>x1+2*PI) x2-=2*PI;
   while(x2<x1) x2+=2*PI;
   while(x3>x2+2*PI) x3-=2*PI;
   while(x3<x2) x3+=2*PI;
   while(x<x1) x+=2*PI;
   while(x>x2) x-=2*PI;

   y0=buf[i0];
   y1=buf[i1];
   y2=buf[i2];
   y3=buf[i3];

   vor=y0*(x-x1)*(x-x2)/((x0-x1)*(x0-x2))
      +y1*(x-x0)*(x-x2)/((x1-x0)*(x1-x2))
      +y2*(x-x0)*(x-x1)/((x2-x0)*(x2-x1));
   nach=y1*(x-x2)*(x-x3)/((x1-x2)*(x1-x3))
       +y2*(x-x1)*(x-x3)/((x2-x1)*(x2-x3))
       +y3*(x-x1)*(x-x2)/((x3-x1)*(x3-x2));
   tmp=(x-x1)/(x2-x1);
   return(vor*(1.-tmp)+nach*tmp);

   } /* interpolAngular() */

double interpol2d(double x,double y,int nx,int ny,double *xbuf,double *ybuf,double *buf)
/*----------------------------------------------------------------------
 * Parabolic blending interpolation of a function of two variables.
 *
 * Input parameters:
 *    x,y    abscissae
 *    nx,ny  number of data points of abscissae
 *    xbuf   array with pre-defined values of x, ascending
 *    ybuf   array with pre-defined values of y, ascending
 *    buf    array with function values, where buf[i+j*nx] is function of xbuf[i] and ybuf[j]
 *
 * Return value: function value at x,y
 *
 * Terms of use:
 *    You are not allowed to claim that this code is the work of anybody else
 *    than Henning Weede.
 *    The author is not liable for any damage caused by this code.
 *    The author has no duty to waist his time answering questions related
 *    to this code.
 *    The GNU Public License applies.
 *
 * Author: Henning Weede
 *----------------------------------------------------------------------
 */{
   int i,ix,iy,j;
   double linDecx,linDecy,linIncx,linIncy,
      quadAftx1,quadAftx2,quadAftx3,quadAfty1,quadAfty2,quadAfty3,
      quadPrex0,quadPrex1,quadPrex2,quadPrey0,quadPrey1,quadPrey2,
      result,tmp,
      shapFuncx[4],shapFuncy[4];

   tmp=fabs(xbuf[nx-1]-xbuf[0])/1000.;
   for(i=1; i<nx; ++i)
      if(xbuf[i]<=xbuf[i-1]+tmp)
         {
         fprintf(stderr,"interpol2d(): 1st abscissa not in ascending order.\n");
         exit(1);
         }
   tmp=fabs(ybuf[ny-1]-ybuf[0])/1000.;
   for(i=1; i<ny; ++i)
      if(ybuf[i]<=ybuf[i-1]+tmp)
         {
         fprintf(stderr,"interpol2d(): 2nd abscissa not in ascending order.\n");
         exit(1);
         }

   for(ix=0; ix<nx-2; ++ix)
      if(xbuf[ix+1]>x) break;

   if(ix>0)
      {
      quadPrex0=(x-xbuf[ix  ])*(x-xbuf[ix+1])/((xbuf[ix-1]-xbuf[ix  ])*(xbuf[ix-1]-xbuf[ix+1]));
      quadPrex1=(x-xbuf[ix-1])*(x-xbuf[ix+1])/((xbuf[ix  ]-xbuf[ix-1])*(xbuf[ix  ]-xbuf[ix+1]));
      quadPrex2=(x-xbuf[ix-1])*(x-xbuf[ix  ])/((xbuf[ix+1]-xbuf[ix-1])*(xbuf[ix+1]-xbuf[ix  ]));
      }
   if(ix<nx-2)
      {
      quadAftx1=(x-xbuf[ix+1])*(x-xbuf[ix+2])/((xbuf[ix  ]-xbuf[ix+1])*(xbuf[ix  ]-xbuf[ix+2]));
      quadAftx2=(x-xbuf[ix  ])*(x-xbuf[ix+2])/((xbuf[ix+1]-xbuf[ix  ])*(xbuf[ix+1]-xbuf[ix+2]));
      quadAftx3=(x-xbuf[ix  ])*(x-xbuf[ix+1])/((xbuf[ix+2]-xbuf[ix  ])*(xbuf[ix+2]-xbuf[ix+1]));
      }

   linIncx=(x-xbuf[ix  ])/(xbuf[ix+1]-xbuf[ix]);
   linDecx=(xbuf[ix+1]-x)/(xbuf[ix+1]-xbuf[ix]);

   if(ix==0)
      {
      shapFuncx[1] = quadAftx1;
      shapFuncx[2] = quadAftx2;
      shapFuncx[3] = quadAftx3;
      }
   else if(ix==nx-2)
      {
      shapFuncx[0] = quadPrex0;
      shapFuncx[1] = quadPrex1;
      shapFuncx[2] = quadPrex2;
      }
   else
      {
      shapFuncx[0] = linDecx * quadPrex0;
      shapFuncx[1] = linDecx * quadPrex1 + linIncx * quadAftx1;
      shapFuncx[2] = linDecx * quadPrex2 + linIncx * quadAftx2;
      shapFuncx[3] = linIncx * quadAftx3;
      }

   for(iy=0; iy<ny-2; ++iy)
      if(ybuf[iy+1]>y) break;

   if(iy>0)
      {
      quadPrey0=(y-ybuf[iy  ])*(y-ybuf[iy+1])/((ybuf[iy-1]-ybuf[iy  ])*(ybuf[iy-1]-ybuf[iy+1]));
      quadPrey1=(y-ybuf[iy-1])*(y-ybuf[iy+1])/((ybuf[iy  ]-ybuf[iy-1])*(ybuf[iy  ]-ybuf[iy+1]));
      quadPrey2=(y-ybuf[iy-1])*(y-ybuf[iy  ])/((ybuf[iy+1]-ybuf[iy-1])*(ybuf[iy+1]-ybuf[iy  ]));
      }
   if(iy<ny-2)
      {
      quadAfty1=(y-ybuf[iy+1])*(y-ybuf[iy+2])/((ybuf[iy  ]-ybuf[iy+1])*(ybuf[iy  ]-ybuf[iy+2]));
      quadAfty2=(y-ybuf[iy  ])*(y-ybuf[iy+2])/((ybuf[iy+1]-ybuf[iy  ])*(ybuf[iy+1]-ybuf[iy+2]));
      quadAfty3=(y-ybuf[iy  ])*(y-ybuf[iy+1])/((ybuf[iy+2]-ybuf[iy  ])*(ybuf[iy+2]-ybuf[iy+1]));
      }

   linIncy=(y-ybuf[iy  ])/(ybuf[iy+1]-ybuf[iy]);
   linDecy=(ybuf[iy+1]-y)/(ybuf[iy+1]-ybuf[iy]);

   if(iy==0)
      {
      shapFuncy[1] = quadAfty1;
      shapFuncy[2] = quadAfty2;
      shapFuncy[3] = quadAfty3;
      }
   else if(iy==ny-2)
      {
      shapFuncy[0] = quadPrey0;
      shapFuncy[1] = quadPrey1;
      shapFuncy[2] = quadPrey2;
      }
   else
      {
      shapFuncy[0] = linDecy * quadPrey0;
      shapFuncy[1] = linDecy * quadPrey1 + linIncy * quadAfty1;
      shapFuncy[2] = linDecy * quadPrey2 + linIncy * quadAfty2;
      shapFuncy[3] = linIncy * quadAfty3;
      }

   result=0.;
   for(i=0; i<4; ++i)
   for(j=0; j<4; ++j)
      if((ix-1+i>=0)&&(ix-1+i<=nx-1)&&(iy-1+j>=0)&&(iy-1+j<=ny-1))
         result += shapFuncx[i]*shapFuncy[j]*buf[ix-1+i + (iy-1+j)*nx];

   return(result);
   } /* interpol2d() */