
#include "spline.h"


/******************************************************************************/
/* CALCUL DES POINTS A TRACER POUR AFFICHER UNE SPLINE                        */
/*                                                                            */
/* Calcule une liste de points pour tracer une courbe spline fermee de points */
/* de controle donnes.                                                        */
/*                                                                            */
/* ENTREES :                                                                  */
/*   -> la liste des points de controle                                       */
/*   -> n : nombre de points a calculer entre deux points de controle         */
/* SORTIES :                                                                  */
/*   -> points [x,y] de la courbe calcules                                    */
/******************************************************************************/

void mexFunction(int nb_sorties,      mxArray *sorties[],
                 int nb_entrees,const mxArray *entrees[]) {

  /* variables locales */

  long nb_points,n,i,j;
  double *sx,*sy,h,*x,*y,s;
  polynome_t *polynomes_x,*polynomes_y,px,py;


  /* lecture des entrees */

  if (nb_entrees!=2) mexErrMsgTxt("Il faut deux entrees !");
  if (nb_sorties>1) mexErrMsgTxt("Il faut au maximum une sortie !");
  if (mxGetN(entrees[0])!=2 || mxGetM(entrees[0])<3)
    mexErrMsgTxt("La matrice des points de la forme est invalide !");
  nb_points=mxGetM(entrees[0]);
  sx=mxGetPr(entrees[0]);
  sy=mxGetPr(entrees[0])+nb_points;
  polynomes_x=(polynome_t *)malloc((nb_points-1)*sizeof(polynome_t));
  polynomes_y=(polynome_t *)malloc((nb_points-1)*sizeof(polynome_t));
  coeffs_spline(sx,polynomes_x,nb_points);
  coeffs_spline(sy,polynomes_y,nb_points);
  if (*mxGetPr(entrees[1])>1.) n=(long)*mxGetPr(entrees[1]);
  else                         n=1;
  h=1./n;


  /* allocation memoire pour les sorties */

  sorties[0]=mxCreateDoubleMatrix(n*(nb_points-1)+1,2,mxREAL);
  x=mxGetPr(sorties[0]);
  y=mxGetPr(sorties[0])+n*(nb_points-1)+1;


  /* construction de x et y */

  for (i=0;i<nb_points-1;i++) {
    px=polynomes_x[i];
    py=polynomes_y[i];
    for (j=0;j<n;j++) {
      s=h*j;
      x[i*n+j]=sx[i]+s*(px.c1+s*(px.c2+s*px.c3));
      y[i*n+j]=sy[i]+s*(py.c1+s*(py.c2+s*py.c3));
    }
  }
  x[n*(nb_points-1)]=x[0];
  y[n*(nb_points-1)]=y[0];


  /* liberation de la memoire */

  free(polynomes_x);
  free(polynomes_y);
}

