#include "all.h"
#include <math.h>
#define PI 3.1415926
int polar_xy(double r, double angle, double *dx, double *dy);
int xy_polar(double dx,double dy,double *radius,double *angle);

int glefitcf_(long *mode,float *x, float *y, long *l, long *m, float *u, float *v, long *n);
fitbez(double **pxv, double **pyv, int **pmv,int *pnp, int multivalued)
{
	double *xv, *yv;
	int *m;
	int n,i;
	float *xin,*yin,*xout,*yout,*xx,*yy,*xxout,*yyout;
	long mode,nsub,outp,ninp;


	if (*pnp>138  || *pnp<=2) return;
	xin = myallocz(*pnp*sizeof(*xin));
	yin = myallocz(*pnp*sizeof(*xin));
	xout = myallocz(300*sizeof(*xin));
	yout = myallocz(300*sizeof(*xin));
	xxout = xout; yyout = yout;
	/* myfree(*pmv); */
	*pmv = myallocz(300*sizeof(i));

	xv = *pxv;
	yv = *pyv;
	xx = xin; yy = yin;
	for (i=0;i<*pnp;i++) {
		*xx++ = *xv++;
		*yy++ = *yv++;
	}
	mode = 1;
	if (multivalued) mode=2;
	nsub = floor(280/(*pnp-1));
	outp = 300;
	outp = (*pnp-1)*nsub+1;
	ninp = *pnp;
	glefitcf_(&mode,xin,yin,&ninp,&nsub,xout,yout,&outp);
	xv = myallocz(300*sizeof(*xv));
	yv = myallocz(300*sizeof(*xv));
	/* myfree(*pxv); 	myfree(*pyv); */
	*pxv = xv;
	*pyv = yv;
	*pnp = outp;
	for (i=0;i<*pnp;i++) {
		*xv++ = *xout++;
		*yv++ = *yout++;
	}
	/* should call FREE to deal old points */
	m = *pmv;
	for (i=0;i<=*pnp;i++) *m++ = 0;	/* no missing values */

	myfree(xin);
	myfree(yin);
	myfree(xxout);
	myfree(yyout);
}