Code:
int knot(int n,int c,int x[60])
{
int nplusc,nplus2,i;
nplusc = n + c;
nplus2 = n + 2;
x[0] = 0;
for (i = 0; i < nplusc; i++){
if ( (i > c) && (i < nplus2) )
x[i+1] = x[i] + 1;
else
x[i] = x[i+1];
}
return 0;
}
float rbasis(int c,float u,int npts,int x[60],float h[28],float r[30])
{
int nplusc;
int i,k;
float d,e;
float sum;
float temp[36];
nplusc = npts + c;
/* printf("knot vector is \n");
for (i = 1; i <= nplusc; i++){
printf(" %d %d \n", i,x[i]);
}
printf("u is %lf \n", u);
*/
/* calculate the first order nonrational basis functions n[i] */
for (i = 0; i < nplusc; i++){
if (( u >= x[i]) && (u < x[i+1]))
temp[i] = 1;
else
temp[i] = 0;
}
/* calculate the higher order nonrational basis functions */
for (k = 0; k < c+1; k++){
for (i = 0; i < nplusc-k+1; i++){
if (temp[i] != 0) { /* if the lower order basis function is zero skip the calculation */
d = ((u-x[i])*temp[i])/(x[i+k]-x[i]);}
else{
d = 0;}
if (temp[i+1] != 0) { /* if the lower order basis function is zero skip the calculation */
e = ((x[i+k]-u)*temp[i+1])/(x[i+k]-x[i]);}
else {
e = 0;} /*line for check */
temp[i] = d + e;
}
}
if (u == (float)x[nplusc]){ /* pick up last point */
temp[npts] = 1;
}
/*
printf("Nonrational basis functions are \n");
for (i=1; i<= npts; i++){
printf("%lf ", temp[i]);
}
printf("\n");
*/
/* calculate sum for denominator of rational basis functions */
sum = 0.;
for (i = 0; i < npts; i++){
sum = sum + temp[i]*h[i];
}
/* form rational basis functions and put in r vector */
for (i = 0; i < npts; i++){
if (sum != 0){
r[i] = (temp[i]*h[i])/(sum);}
else
r[i] = 0;
}
return 0;
}
float rbspline()
{
int i,j,icount,jcount;
int i1;
int x[60]; /* data points store vector */
int nplusc;
int ispline;
int npts,k;
int p1;
float step;
float u;
float nbasis[30];
float temp;
float b[60]; /* control vertices matrix*/
float p[203]; /* allows for up to 200 points on curve */
float h[28];
FILE *splineinput;
FILE *airfoilcoordinates;
splineinput = fopen("controlpoints2.txt","r");
airfoilcoordinates = fopen("airfoilcoordinates.dat","w");
npts = 28;
k = 6; /* third ordech=getc();r, change for other orders */
p1 = 100; /* eleven points on curve */
for (ispline = 0; ispline < 2*npts; ispline++)
{
fscanf(splineinput,"%f\n ",&b[ispline]);
}
printf("\nPolygon points\n");
for (ispline = 0; ispline < 2*npts; ispline=ispline+2)
{
printf(" %1.5f %1.5f \n",b[ispline],b[ispline+1]);
}
for (ispline=0; ispline < npts; ispline++)
{
h[ispline] = 1.0;
}
for (ispline = 0; ispline < 2*p1; ispline++)
{
p[ispline] = 0.;
}
nplusc = npts + k;
/* zero and redimension the knot vector and the basis array */
for(i = 0; i < npts; i++){
nbasis[i] = 0.;
}
for(i = 0; i < nplusc; i++){
x[i] = 0.;
}
/* generate the uniform open knot vector */
knot(npts,k,x);
/*
printf("The knot vector is ");
for (i = 1; i <= nplusc; i++){
printf(" %d ", x[i]);
}
printf("\n");
*/
icount = 0;
/* calculate the points on the rational B-spline curve */
u = 0;
step = ((float)x[nplusc])/((float)(p1-1));
for (i1 = 0; i1 < p1; i1++){
if ((float)x[nplusc] - u < 5e-6){
u = (float)x[nplusc];
}
rbasis(k,u,npts,x,h,nbasis); /* generate the basis function for this value of t */
/*
printf("u = %f \n",u);
printf("nbasis = ");
for (i = 1; i <= npts; i++){
printf("%f ",nbasis[i]);
}
printf("\n");
*/
for (j = 0; j < 2; j++){ /* generate a point on the curve */
jcount = j;
p[icount+j] = 0.;
for (i = 0; i < npts; i++){ /* Do local matrix multiplication */
temp = nbasis[i]*b[jcount];
p[icount + j] = p[icount + j] + temp;
/*
printf("jcount,nbasis,b,nbasis*b,p = %d %f %f %f %f\n",jcount,nbasis[i],b[jcount],temp,p[icount+j]);
*/
jcount = jcount + 2;
}
}
/*
printf("icount, p %d %f %f %l \n",icount,p[icount+1],p[icount+2],p[icount+3]);
*/
icount = icount + 2;
u = u + step;
}
printf("\nHomogeneous weighting vector is \n");
for (ispline = 0; ispline < npts; ispline++){
printf(" %1.4f ", h[ispline]);
}
printf("\n");
/*printf("\nCurve points\n"); */
for (ispline = 0; ispline < 2*p1; ispline=ispline+2)
{
fprintf(airfoilcoordinates,"%1.5f %1.5f \n",p[ispline],p[ispline+1]);
}
fclose(splineinput);
fclose(airfoilcoordinates);
return 0;
}