We Have Mail

Listing 8 ratint.c

      Function RATINT} for the RPFT Code
/******************************************************
 * RATINT - Diagonal rational function interpolation in
 * the arrays xa[1..n] and ya[1..n].
 *****************************************************/
void ratint(double xa[], double ya[], double *c,
       double *d, int n, double x, double *y)
{ int m,i,ns=1;
  double w,t,hh,h,dd;
  static double miny=1.e99;

  if (miny>1.e90) for (i=1; i<=n; ++i)
                     if (ya[i]<miny) miny=ya[i];
  hh=fabs(x-xa[1]);
  for (i=1;i<=n;i++)
  { h=fabs(x-xa[i]);
    if (h == 0.0) {*y=ya[i]; return; }
    else if (h < hh) { ns=i; hh=h; }
    c[i]=ya[i]-miny; d[i]=ya[i]-miny+1.e-50;
}
*y=ya[ns--] - miny;
for (m=1;m<n;m++)
{ for (i=1;i<=n-m;i++)
  { w=c[i+1]-d[i] ; h=xa[i+m]-x; t=(xa[i]-x)*d[i]/h;
    dd=t-c[i+1];
    if (fabs(t)>1.e15)
      fprintf(stderr,"Probable loss of accuracy in"
       "RATINT. fabs(t) > 1.e15 for X = %.8G\n",x);
    if (dd == 0.0)
    { fprintf(stderr,"Error in routine ratint. The"
         "function may have a pole at x=%.8G\n",x);
      exit(1);
    }
    dd=w/dd; d[i]=c[i+1]*dd; c[i]=t*dd;
  }
  *y += (2*ns < (n-m) ? c[ns+l] : d[ns--]);
 }
 *y += miny; return;
}

/* End of File */