/************************************************************************/
/*                                                                      */
/* INTERPOLATION                                                        */
/*                                                                      */
/* Fonction qui calcule une fonction interpolante "y" a partir des      */
/* valeurs tabulees ys=f(xs) pour les points stockés dans le vecteur    */
/* "x".                                                                 */
/* On utilise un algorithme d'interpolation cubique par des splines     */
/* naturelles.                                                          */
/*                                                                      */
/************************************************************************/

void interpolation(double y[],double xs[],double ys[],int n,double x[]) {
  int i,j,cour,prec;
  double p,sig,*ys2,*temp,a,b,c,d,e,f,g,a0,a1,a2,a3,xc;

  /* Calcul des derivees secondes au points xs[i] */
  ys2=(double *)malloc(n*sizeof(double));
  temp=(double *)malloc(n*sizeof(double));
  ys2[0]=temp[0]=0.0;
  for (i=1;i<n-1;i++) {
    sig=(xs[i]-xs[i-1])/(xs[i+1]-xs[i-1]);
    p=sig*ys2[i-1]+2.0;
    ys2[i]=(sig-1.0)/p;
    temp[i]=(ys[i+1]-ys[i])/(xs[i+1]-xs[i])-(ys[i]-ys[i-1])/(xs[i]-xs[i-1]);
    temp[i]=(6.0*temp[i]/(xs[i+1]-xs[i-1])-sig*temp[i-1])/p;
  }
  ys2[n-1]=0.0;
  for (j=n-2;j>=0;j--) ys2[j]=ys2[j]*ys2[j+1]+temp[j];
  free(temp);

  /* Calcul de la fonction interpolante */
  cour=0;
  for (j=0;j<(n-1);j++) {

    /* calcul des coefficients du polynome */
    a=xs[j];
    b=xs[j+1];
    c=b-a;
    d=ys[j];
    e=ys[j+1];
    f=ys2[j];
    g=ys2[j+1];
    a0=(b*d-a*e+CUBE(b)*f/6-CUBE(a)*g/6)/c+c*(a*g-b*f)/6;
    a1=(e-d-CARRE(b)*f/2+CARRE(a)*g/2)/c+c*(f-g)/6;
    a2=(b*f-a*g)/(2*c);
    a3=(g-f)/(6*c);

    prec=cour;
    while (x[cour]!=xs[j+1]) cour++;

    /* calcul de la valeur du polynome au point x[i] */
    for (i=prec;i<cour;i++) {
      xc=x[i];
      y[i]=a0+a1*xc+a2*CARRE(xc)+a3*CUBE(xc);
    }

  }
  y[cour]=ys[n-1];

  free(ys2);
}
