/* @(#)interpolate.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:59 */
/*  
    
    Otmar Stahl

    interpolate.c
    
    spline interpolation of array 

*/

/* FEROS specific includes */

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

int interpolate
#ifdef __STDC__
(
 float in[], int howmany_in, float by_howmuch, float out[], int howmany_out
 )
#else
     (
      in, howmany_in, by_howmuch, out, howmany_out
      )
     float in[], out[], by_howmuch;
     int howmany_in, howmany_out;
#endif

{

  /* simple spline interpolation on equidistant grid */

  float diffx, sum, *xpos, *y2, *tmp_in, x, startpos, stepsize;
  int i, in_entries, out_entries;

  diffx = -by_howmuch;
  in_entries = howmany_in;
  out_entries = howmany_out;
  stepsize = (float) out_entries / in_entries;

  startpos = 0.;
  while (startpos >= -0.5)
    startpos -= stepsize;
  startpos += stepsize;
  sum = 0.0;

  xpos = vector(1, in_entries);
  y2 = vector(1, in_entries);
  tmp_in = vector(1, in_entries);
 
  for (i = 1; i <= in_entries; i++)
    {
      tmp_in[i] = in[i - 1];
      xpos[i] = startpos + (((float) i) - 1.) * stepsize;
    }
  spline(xpos, tmp_in, in_entries, 0., -0., y2);
  for (i = 0; i < out_entries; i++)
    {
      x = ((float) i)  + diffx;
      splint(xpos, tmp_in, y2, in_entries, x, &out[i]);
    }
  for (i = 0; i < out_entries; i++)
    {
      sum += out[i];
    }
  for (i = 0; i < out_entries; i++)
    {
      out[i] /= sum;
    }
  free_vector(xpos, 1, in_entries); 
  free_vector(y2, 1, in_entries);
  free_vector(tmp_in, 1, in_entries);
  return 0;
}