/* @(#)mutil.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:59 */
/* 
   Cristian Levin - ESO La Silla
   
   mutil.c

   miscellaneous numerical routines
   mainly based on Numerical Recipes

*/

/* system includes */

#include <stdio.h>
#include <math.h>

/* FEROS specific includes */

#include <mutil.h>
#include <proto_nrutil.h>
#include <proto_mutil.h>

/******************************************************

 * lfit(): general linear least squares solution.
 *
 ******************************************************/
void lfit
#ifdef __STDC__
(
 double x[], double y[], int ndata, double a[],
 int mfit, void (*funcs) (double, double *afunc, int mfit)
 )
#else
     (
      x, y, ndata, a, mfit, funcs
      )
     
     double x[], y[], a[];
     int ndata, mfit;
     void (*funcs) ();

#endif

{
  int k, j, i;
  double **covar, **beta, *afunc;

  beta = dmatrix (1, mfit, 1, 1);
  covar = dmatrix (1, mfit, 1, mfit);
  afunc = dvector (1, mfit);

  for (j = 1; j <= mfit; j++)
    {
      for (k = 1; k <= mfit; k++)
	covar[j][k] = 0.0;
      beta[j][1] = 0.0;
    }

  for (i = 0; i < ndata; i++)
    {
      (*funcs) (x[i], afunc, mfit);
      for (j = 1; j <= mfit; j++)
	{
	  for (k = 1; k <= j; k++)
	    covar[j][k] += afunc[j] * afunc[k];
	  beta[j][1] += y[i] * afunc[j];
	}
    }

for (j = 2; j <= mfit; j++)
  for (k = 1; k <= j - 1; k++)
     covar[k][j] = covar[j][k];

     gaussj (covar, mfit, beta, 1);	/* Matrix solution */

     for (j = 1; j <= mfit; j++)	/* solution to coefficients a */
     a[j] = beta[j][1];

     free_dvector (afunc, 1, mfit);
     free_dmatrix (beta, 1, mfit, 1, 1);
     free_dmatrix (covar, 1, mfit, 1, mfit);
     }


/******************************************************************

 * gaussj(): linear equation solution by Gauss-Jordan elimination
 *
 ******************************************************************/

#define SWAP(a, b) 	{float temp = (a); (a) = (b); (b) = temp; }

void gaussj 
#ifdef __STDC__
(
 double **a, int n, double **b, int m
 )
#else
     (
      a, n, b, m
      )
     double **a, **b;
     int m, n;
#endif

     /* IN    : a[1..n][1..n] */
     /* IN/OUT: b[1..n][1..m] */

{
  int *indxc, *indxr, *ipiv;	/* for bookeeping in the pivoting */
  int i, icol, irow, j, k;
  double big, dum, pivinv;

  indxc = ivector (1, n);
  indxr = ivector (1, n);
  ipiv = ivector (1, n);

  for (i = 1; i <= n; i++)
    ipiv[i] = 0.0;

  for (i = 1; i <= n; i++)
    {				/* main loop over columns to be reduced */
      big = 0.0;
      for (j = 1; j <= n; j++)
	if (ipiv[j] != 1)
	  for (k = 1; k <= n; k++)
	    {
	      if (ipiv[k] == 0)
		{
		  if (fabs (a[j][k]) >= big)
		    {
		      big = fabs (a[j][k]);
		      irow = j;
		      icol = k;
		    }
		}
	      else if (ipiv[k] > 1)
		{
		  nrerror ("Tolerance too small for fitting");
		  return;
		}
	    }
      ++(ipiv[icol]);

      if (irow != icol)
	{
	  for (j = 1; j <= n; j++)
	    SWAP (a[irow][j], a[icol][j]);
	  for (j = 1; j <= m; j++)
	    SWAP (b[irow][j], b[icol][j]);
	}

      indxr[i] = irow;
      indxc[i] = icol;
      if (a[icol][icol] == 0.0)
	{
	  nrerror ("Tolerance too small for fitting");
	  return;
	}

      pivinv = 1.0 / a[icol][icol];
      a[icol][icol] = 1.0;
      for (j = 1; j <= n; j++)
	a[icol][j] *= pivinv;
      for (j = 1; j <= m; j++)
	b[icol][j] *= pivinv;

      for (j = 1; j <= n; j++)
	if (j != icol)
	  {
	    dum = a[j][icol];
	    a[j][icol] = 0.0;
	    for (k = 1; k <= n; k++)
	      a[j][k] -= a[icol][k] * dum;
	    for (k = 1; k <= m; k++)
	      b[j][k] -= b[icol][k] * dum;
	  }
    }				/* main loop */

  for (i = n; i >= 1; i--)
    if (indxr[i] != indxc[i])
      for (k = 1; k <= n; k++)
	SWAP (a[k][indxr[i]], a[k][indxc[i]]);

  free_ivector (ipiv, 1, n);
  free_ivector (indxc, 1, n);
  free_ivector (indxr, 1, n);
}

/********************************************************************

 * polint(): polynomial interpolation or extrapolation from N input
 *           points, using the Neville's algorithm.
 *
 ********************************************************************/
void  polint 
#ifdef __STDC__
(
 float xa[], float ya[], int n, float x, float *y, float *dy
 )
#else
     (
      xa, ya, n, x, y, dy
      )
     float xa[], ya[], x, *y, *dy;
     int n;
#endif

{
  int i, m, ns = 1;
  float den, dif, dift, ho, hp, w;
  float *c, *d;

  dif = fabs (x - xa[1]);
  c = fvector (1, n);
  d = fvector (1, n);
  for (i = 1; i <= n; i++)
    {				/* find the index ns of the closest */
      if ((dift = fabs (x - xa[i])) < dif)
	{			/* table entry. */
	  ns = i;
	  dif = dift;
	}
      c[i] = ya[i];		/* initialize the tableau of c's and d's */
      d[i] = ya[i];
    }
  *y = ya[ns--];		/* initial aprox. of y */
  for (m = 1; m < n; m++)
    {
      for (i = 1; i <= n - m; i++)
	{
	  ho = xa[i] - x;
	  hp = xa[i + m] - x;
	  w = c[i + 1] - d[i];
	  if ((den = ho - hp) == 0.0)
	    {
	      nrerror ("Error in polynomial interpolation");
	      return;
	    }
	  den = w / den;
	  d[i] = hp * den;
	  c[i] = ho * den;
	}
      *y += (*dy = (2 * ns < (n - m) ? c[ns + 1] : d[ns--]));
    }
  free_fvector (c, 1, n);
  free_fvector (d, 1, n);
}

/************************************************************

 * eval_fpoly(), eval_dpoly(): fast polynomial evaluation.
 *
 ************************************************************/
float eval_fpoly 
#ifdef __STDC__
(
 float x, float *a, int n
 )
#else
     (
      x, a, n
      )
     float x, *a;
     int n;
#endif

     /*  *a: coefficients (range: 1..n) */

{
  int i;
  float y = 0.0;

  for (i = n; i >= 1; i--)
    y = x * y + a[i];

  return (y);
}

double eval_dpoly 
#ifdef __STDC__
(
 double x, double *a, int n
 )
#else
     (x, a, n
      )
     double x, *a;
     int n;
#endif

     /* *a: coefficients (range: 1..n) */

{
  int i;
  double y = 0.0;

  for (i = n; i >= 1; i--)
    y = x * y + a[i];

  return (y);
}


void dpoly 
#ifdef __STDC__
(
 double x, double p[], int np
 )
#else
     (
      x, p, np
      )
     double x, p[];
     int np;
#endif
{
  int j;

  p[1] = 1.0;
  for (j = 2; j <= np; j++)
    p[j] = p[j - 1] * x;
}

/************************************************************

 * median(): efficient search of the median.
 *
 ************************************************************/

float median 
#ifdef __STDC__
(
 float x[], int n
 )
#else
     (
      x, n
      )
     float x[];
     int n;
#endif
{
  int n2, n2p;
#ifdef __STDC__
  void piksrt (int, float[]);
#else
  void piksrt ();
#endif

  piksrt (n, x);
  n2p = (n2 = n / 2) + 1;
  return ((n % 2) ? x[n2p] : 0.5 * (x[n2] + x[n2p]));
}

/************************************************************

 * piksrt(): sorting by straight insertion.
 *
 ************************************************************/

void piksrt
#ifdef __STDC__
(
 int n, float arr[]
 )
#else
     (
      n, arr
      )
     float arr[];
     int n;
#endif

{
  int i, j;
  float a;

  for (j = 2; j <= n; j++)
    {
      a = arr[j];
      i = j - 1;
      while (i > 0 && arr[i] > a)
	{
	  arr[i + 1] = arr[i];
	  i--;
	}
      arr[i + 1] = a;
    }
}

void dsvdcmp 
#ifdef __STDC__
(
 double **a, int m, int n, double w[], double **v
 )
#else
     (
      a, m, n, w, v
      )
     double **a, w[], **v;
     int m, n;
#endif

{

  int flag, i, its, j, jj, k, l, nm;
  double anorm, c, f, g, h, s, scale, x, y, z, *rv1;

  rv1 = dvector (1, n);
  g = scale = anorm = 0.0;
  for (i = 1; i <= n; i++)
    {
      l = i + 1;
      rv1[i] = scale * g;
      g = s = scale = 0.0;
      if (i <= m)
	{
	  for (k = i; k <= m; k++)
	    scale += fabs (a[k][i]);
	  /*	  if (scale)   */
	  if ((int) scale)   
	    {
	      for (k = i; k <= m; k++)
		{
		  a[k][i] /= scale;
		  s += a[k][i] * a[k][i];
		}
	      f = a[i][i];
	      g = -SIGN (sqrt (s), f);
	      h = f * g - s;
	      a[i][i] = f - g;
	      for (j = l; j <= n; j++)
		{
		  for (s = 0.0, k = i; k <= m; k++)
		    s += a[k][i] * a[k][j];
		  f = s / h;
		  for (k = i; k <= m; k++)
		    a[k][j] += f * a[k][i];
		}
	      for (k = i; k <= m; k++)
		a[k][i] *= scale;
	    }
	}
      w[i] = scale * g;
      g = s = scale = 0.0;
      if (i <= m && i != n)
	{
	  for (k = l; k <= n; k++)
	    scale += fabs (a[i][k]);
	  /*	  if (scale) */
	  if ((int) scale)   
	    {
	      for (k = l; k <= n; k++)
		{
		  a[i][k] /= scale;
		  s += a[i][k] * a[i][k];
		}
	      f = a[i][l];
	      g = -SIGN (sqrt (s), f);
	      h = f * g - s;
	      a[i][l] = f - g;
	      for (k = l; k <= n; k++)
		rv1[k] = a[i][k] / h;
	      for (j = l; j <= m; j++)
		{
		  for (s = 0.0, k = l; k <= n; k++)
		    s += a[j][k] * a[i][k];
		  for (k = l; k <= n; k++)
		    a[j][k] += s * rv1[k];
		}
	      for (k = l; k <= n; k++)
		a[i][k] *= scale;
	    }
	}
      anorm = DMAX (anorm, (fabs (w[i]) + fabs (rv1[i])));
    }
  for (i = n; i >= 1; i--)
    {
      if (i < n)
	{
	  /*	  if (g) */
	  if ((int) g)
	    {
	      for (j = l; j <= n; j++)
		v[j][i] = (a[i][j] / a[i][l]) / g;
	      for (j = l; j <= n; j++)
		{
		  for (s = 0.0, k = l; k <= n; k++)
		    s += a[i][k] * v[k][j];
		  for (k = l; k <= n; k++)
		    v[k][j] += s * v[k][i];
		}
	    }
	  for (j = l; j <= n; j++)
	    v[i][j] = v[j][i] = 0.0;
	}
      v[i][i] = 1.0;
      g = rv1[i];
      l = i;
    }
  for (i = IMIN (m, n); i >= 1; i--)
    {
      l = i + 1;
      g = w[i];
      for (j = l; j <= n; j++)
	a[i][j] = 0.0;
      if ((int) g)
	{
	  g = 1.0 / g;
	  for (j = l; j <= n; j++)
	    {
	      for (s = 0.0, k = l; k <= m; k++)
		s += a[k][i] * a[k][j];
	      f = (s / a[i][i]) * g;
	      for (k = i; k <= m; k++)
		a[k][j] += f * a[k][i];
	    }
	  for (j = i; j <= m; j++)
	    a[j][i] *= g;
	}
      else
	for (j = i; j <= m; j++)
	  a[j][i] = 0.0;
      /*      ++a[i][i];  */
      a[i][i] = a[i][i]+1;
    }
  for (k = n; k >= 1; k--)
    {
      for (its = 1; its <= 30; its++)
	{
	  flag = 1;
	  for (l = k; l >= 1; l--)
	    {
	      nm = l - 1;
	      if ((double) (fabs (rv1[l]) + anorm) == anorm)
		{
		  flag = 0;
		  break;
		}
	      if ((double) (fabs (w[nm]) + anorm) == anorm)
		break;
	    }
	  if (flag)
	    {
	      c = 0.0;
	      s = 1.0;
	      for (i = l; i <= k; i++)
		{
		  f = s * rv1[i];
		  rv1[i] = c * rv1[i];
		  if ((double) (fabs (f) + anorm) == anorm)
		    break;
		  g = w[i];
		  h = dpythag (f, g);
		  w[i] = h;
		  h = 1.0 / h;
		  c = g * h;
		  s = -f * h;
		  for (j = 1; j <= m; j++)
		    {
		      y = a[j][nm];
		      z = a[j][i];
		      a[j][nm] = y * c + z * s;
		      a[j][i] = z * c - y * s;
		    }
		}
	    }
	  z = w[k];
	  if (l == k)
	    {
	      if (z < 0.0)
		{
		  w[k] = -z;
		  for (j = 1; j <= n; j++)
		    v[j][k] = -v[j][k];
		}
	      break;
	    }
	  if (its == 30)
	    nrerror ("no convergence in 30 dsvdcmp iterations");
	  x = w[l];
	  nm = k - 1;
	  y = w[nm];
	  g = rv1[nm];
	  h = rv1[k];
	  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
	  g = dpythag (f, 1.0);
	  f = ((x - z) * (x + z) + h * ((y / (f + SIGN (g, f))) - h)) / x;
	  c = s = 1.0;
	  for (j = l; j <= nm; j++)
	    {
	      i = j + 1;
	      g = rv1[i];
	      y = w[i];
	      h = s * g;
	      g = c * g;
	      z = dpythag (f, h);
	      rv1[j] = z;
	      c = f / z;
	      s = h / z;
	      f = x * c + g * s;
	      g = g * c - x * s;
	      h = y * s;
	      y *= c;
	      for (jj = 1; jj <= n; jj++)
		{
		  x = v[jj][j];
		  z = v[jj][i];
		  v[jj][j] = x * c + z * s;
		  v[jj][i] = z * c - x * s;
		}
	      z = dpythag (f, h);
	      w[j] = z;
	      if ((int) z)
		{
		  z = 1.0 / z;
		  c = f * z;
		  s = h * z;
		}
	      f = c * g + s * y;
	      x = c * y - s * g;
	      for (jj = 1; jj <= m; jj++)
		{
		  y = a[jj][j];
		  z = a[jj][i];
		  a[jj][j] = y * c + z * s;
		  a[jj][i] = z * c - y * s;
		}
	    }
	  rv1[l] = 0.0;
	  rv1[k] = f;
	  w[k] = x;
	}
    }
  free_dvector (rv1, 1, n);
}

void dsvbksb 
#ifdef __STDC__
(
 double **u, double w[], double **v, int m,
 int n, double b[], double x[]
 )
#else
     (
      u, w, v, m, n, b, x
      )
     double **u, w[], **v, b[], x[];
     int m, n;
#endif

{

  int jj, j, i;
  double s, *tmp;

  tmp = dvector (1, n);
  for (j = 1; j <= n; j++)
    {
      s = 0.0;
      if ((int) w[j])
	{
	  for (i = 1; i <= m; i++)
	    s += u[i][j] * b[i];
	  s /= w[j];
	}
      tmp[j] = s;
    }
  for (j = 1; j <= n; j++)
    {
      s = 0.0;
      for (jj = 1; jj <= n; jj++)
	s += v[j][jj] * tmp[jj];
      x[j] = s;
    }
  free_dvector (tmp, 1, n);
}

double dpythag 
#ifdef __STDC__
(
 double a, double b
 )
#else
     (
      a, b
      )
     double a, b;
#endif
{
  double absa, absb;
  absa = fabs (a);
  absb = fabs (b);
  if (absa > absb)
    return absa * sqrt (1.0 + DSQR (absb / absa));
  else
    return (absb == 0.0 ? 0.0 : absb * sqrt (1.0 + DSQR (absa / absb)));
}

float select_pos 
#ifdef __STDC__
(
 int k, int n, float arr[]
 )
#else
     (
      k, n, arr
      )
     int k, n;
     float arr[];
#endif

     /*      select the kth smallest value in arr[1...n], NR, pg. 342 */

#define SWAPD(a, b) 	{float temp = (a); (a) = (b); (b) = temp; }

{

  int i, ir, j, l, mid;
  /*  double a, temp; */
  float a; 

  l = 1;
  ir = n;
  for (;;)
    {
      if (ir <= l + 1)
	{
	  if (ir == l + 1 && arr[ir] < arr[l])
	    {
	      SWAPD (arr[l], arr[ir])
		}
	  return arr[k];
	}
      else
	{
	  mid = (l + ir) >> 1;
	  SWAPD (arr[mid], arr[l + 1])
	    if (arr[l + 1] > arr[ir])
	      {
		SWAPD (arr[l + 1], arr[ir])
		  }
	  if (arr[l] > arr[ir])
	    {
	      SWAPD (arr[l], arr[ir])
		}
	  if (arr[l + 1] > arr[l])
	    {
	      SWAPD (arr[l + 1], arr[l])
		}
	  i = l + 1;
	  j = ir;
	  a = arr[l];
	  for (;;)
	    {
	      do
		i++;
	      while (arr[i] < a);
	      do
		j--;
	      while (arr[j] > a);
	      if (j < i)
		break;
	      SWAPD (arr[i], arr[j])
		}
	  arr[l] = arr[j];
	  arr[j] = a;
	  if (j >= k)
	    ir = j - 1;
	  if (j <= k)
	    l = i;
	}
    }
} /* end of select_pos */
#undef SWAPD
#undef SWAP
/*---------------------------------------------------------------------------*/
float pik_median 
#ifdef __STDC__
(
 int n, float arr[]
 )
#else
     ( 
      n, arr 
      )
     float arr[];
     int   n;
#endif

     /*              find median of distribution (use piksort)       */

{
  int i, j, index;
  float a, b[101];


  for (j = 0; j < n; j++)
    b[j] = arr[j];

  for (j = 1; j < n; j++)
    {
      a = b[j];
      i = j - 1;
      while (i >= 0 && b[i] > a)
	{
	  b[i + 1] = b[i];
	  i--;
	}
      b[i + 1] = a;
    }
  index = (n - 1) / 2;
  return(b[index]);
}

float heap_median 
#ifdef __STDC__
(
 int n, float arr[]
 )
#else
     (
      n, arr
      )
     float arr[];
     int   n;
#endif
     /*              find median of distribution (heapsort)           */
{
  int i, j, index;
  float a, b[101];

  for (j = 0; j < n; j++)
    b[j] = arr[j];

  for (j = 1; j < n; j++)
    {
      a = b[j];
      i = j - 1;
      while (i >= 0 && b[i] > a)
	{
	  b[i + 1] = b[i];
	  i--;
	}
      b[i + 1] = a;
    }
  index = (n - 1) / 2;
  return(b[index]);
}

/*---------------------------------------------------------------------------*/
void fit_poly 
#ifdef __STDC__
(
 float inimage[], float ximage[], int npixi, 
 float outimage[], float ximageo[], int npixo, 
 int fit_deg, double ad[]
 )
#else
     ( 
      inimage, ximage, npixi, outimage, ximageo, npixo, fit_deg, ad 
      )
     float inimage[], ximage[], outimage[], ximageo[];
     int npixi, npixo, fit_deg;
     double ad[];
#endif

     /*      fit a polynomial (Numerical Recipes), modified for echelle!  */

{
  float xout;
  double *xin, *yin;
  int i;

  xin = dvector (0, npixi);
  yin = dvector (0, npixi);

  for (i = 0; i < npixi; i++)
    {
      xin[i] = (double) ximage[i];
      yin[i] = (double) inimage[i];
    }

  lfit (xin, yin, npixi, ad, fit_deg, dpoly);

  for (i = 0; i < npixo; i++)
    {
      xout = ximageo[i];
      outimage[i] = (float) eval_dpoly (xout, ad, fit_deg);
    }

  free_dvector(xin, 0, npixi), free_dvector(yin, 0, npixi);

}

void spline
#ifdef __STDC__
(
 float x[], float y[], int n, float yp1, float ypn, float y2[]
 )
#else
     (
      x, y, n, yp1, ypn, y2
      )
     float x[], y[], yp1, ypn, y2[];
     int n;
#endif
{
  int i,k;
  float p,qn,sig,un,*u;

  u=vector(1,n-1);
  if (yp1 > 0.99e30)
    y2[1]=u[1]=0.0;
  else {
    y2[1] = -0.5;
    u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
  }
  for (i=2;i<=n-1;i++) {
    sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
    p=sig*y2[i-1]+2.0;
    y2[i]=(sig-1.0)/p;
    u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
    u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
  }
  if (ypn > 0.99e30)
    qn=un=0.0;
  else {
    qn=0.5;
    un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
  }
  y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  for (k=n-1;k>=1;k--)
    y2[k]=y2[k]*y2[k+1]+u[k];
  free_vector(u,1,n-1);
}

void splint
#ifdef __STDC__
(
 float xa[], float ya[], float y2a[], int n, float x, float *y
 )
#else
     (
      xa, ya, y2a, n, x, y
      )
     float xa[], ya[], y2a[], x, *y;
     int n;
#endif
{
  int klo,khi,k;
  float h,b,a;

  klo=1;
  khi=n;
  while (khi-klo > 1) {
    k=(khi+klo) >> 1;
    if (xa[k] > x) khi=k;
    else klo=k;
  }
  h=xa[khi]-xa[klo];
  if (h == 0.0) nrerror("Bad xa input to routine splint");
  a=(xa[khi]-x)/h;
  b=(x-xa[klo])/h;
  *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}


int splie2 
#ifdef __STDC__
(
 float x1a[], float x2a[], float **ya, int m, int n, float **y2a
 )
#else
     (
      x1a, x2a, ya, m, n, y2a
      )
     float x1a[], x2a[], **ya, **y2a;
     int m, n; 
#endif
{
  int j;
  for (j=1; j<m; j++)
    {
      spline(x2a, ya[j], n, 1.0e30, 1.0e30, y2a[j]); 
    }
  return 0;
}

int splin2
#ifdef __STDC__
(
 float x1a[], float x2a[], float **ya, float **y2a, int m, int n,
 int num, float x1[], float x2, float y[]
 )
#else
     (
      x1a, x2a, ya, y2a, m, n, num, x1, x2, y
      )
     float x1a[], x2a[], **ya, **y2a, x1[], x2, *y;
     int m, n, num;
#endif
{
  int j;
  float *ytmp,*yytmp;
  
  ytmp=vector(1,n);
  yytmp=vector(1,n);
  
  for (j=1;j<=m;j++)
    splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]);
  
  spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp); 

  for (j=0;j<num;j++)
    {
      splint(x1a,yytmp,ytmp,m,x1[j],&y[j]);
    }
  free_vector(yytmp,1,n);
  free_vector(ytmp,1,n);
  return 0;
}

int splin3
#ifdef __STDC__
(
 float x1a[], float x2a[], float **ya, float **y2a, int m, int n,
 int num, float x1[], float x2, float y[]
 )
#else
     (
      x1a, x2a, ya, y2a, m, n, num, x1, x2, y
      )
     float x1a[], x2a[], **ya, **y2a, x1[], x2, *y;
     int m, n, num;
#endif
{

  int j;
  float *ytmp,*yytmp;

  int klo,khi,k;
  float h,b,a;

  klo=1;
  khi=m;

  ytmp=vector(1,m);
  yytmp=vector(1,m);
  for (j=1;j<=m;j++)
    splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]);

  spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp); 

  while (khi-klo > 1) 
    {
      k=(khi+klo) >> 1;
      if (x1a[k] > x1[0]) khi=k;
      else klo=k;
    }
  h=x1a[khi]-x1a[klo];
  if (h == 0.0) nrerror("Bad xa input to routine splint");

  for (j=0;j<num;j++)
    {
      if ( (x1[j] > x1a[1]) && (x1[j] < x1a[m]) )
	{
	  if (x1[j] > x1a[khi])
	    {
	      while (x1[j] > x1a[khi])
		{
		  klo++;
		  khi++;
		}
	      h=x1a[khi]-x1a[klo];
	      if (h == 0.0) nrerror("Bad xa input to routine splint");
	    }
	}
      
      a=(x1a[khi]-x1[j])/h;
      b=(x1[j]-x1a[klo])/h;
      y[j] = a*yytmp[klo] + b*yytmp[khi] + 
	((a*a*a-a)*ytmp[klo] + (b*b*b-b)*ytmp[khi])*(h*h)/6.0;
    }
  free_vector(yytmp,1,m);
  free_vector(ytmp,1,m);
  return 0;
}