Code for obtaining the basis function by using the NLE

// Nurbsbasis. cpp: defines the entry point of the console application.
//

# Include "stdafx. H"
# Include <stdio. h>
# Include <iostream>
Using namespace STD;

Int findspan (int n, int P, double U, double * U)
{
If (u = U [n + 1]) return N;
Int low = P; int high = n + 1;
Int mid = (low + high)/2;
While (U <u [Mid] | u> = U [Mid + 1]) {
If (U <u [Mid]) High = mid;
Else low = mid;
Mid = (low + high)/2;
}
Return mid;
}

Void basisfuns (int I, double U, int P, double * u, double * n)
{
Double saved = 0.0, temp = 0.0;
Double * Left = new double [p + 1];
Double * Right = new double [p + 1];
N [0] = 1.0;
For (Int J = 1; j <= P; j ++)
{
Left [J] = u-u [I + 1-J];
Right [J] = U [I + J]-U;
Saves = 0.0;
For (INT r = 0; r <j; r ++)
{
// Lower triangle
Temp = right [R + 1] + left [J-R];
Temp = N [R]/temp;

// Upper triangle
N [R] = saved + right [R + 1] * temp;
Saved = left [J-R] * temp;
}
N [J] = saved;
}
}

Void dersbasisfuns (int I, double U, int P, int N, double * u, double ** & ders)
{
// Compute nonzero basis functions and their
// Derivatives. First Section is a2.2 modified
// To store functions and knot differences
// Input: I, U, P, N, u
// Output: ders
Double saved, temp;
Double ** NDU;
Double **;

// Set up a [p + 1] * [p + 1] matrix NDU
NDU = new double * [p + 1];
If (NDU = NULL)
{
Cout <"error" <Endl;
}

For (Int J = 0; j <= P; j ++)
{
NDU [J] = new double [p + 1];
If (NDU [J] = NULL)
{
Cout <"error" <Endl;
For (int K = 0; k <= J; k ++)
{
Delete [] NDU [k];
}
Delete [] NDU;
}
}

// Set up a [2] * [p + 1] matrix

A = new double * [2];
If (A = NULL)
{
Cout <"error" <Endl;
}

For (Int J = 0; j <= 1; j ++)
{
A [J] = new double [p + 1];
If (A [J] = NULL)
{
Cout <"error" <Endl;
For (int K = 0; k <= J; k ++)
{
Delete [] a [k];
}
Delete [];
}
}

// Set up a [n + 1] * [p + 1] matrix ders

/* Ders = new double * [n + 1];
If (Ders = NULL)
{
Cout <"error" <Endl;
}

For (Int J = 0; j <= 1; j ++)
{
Ders [J] = new double [p + 1];
If (Ders [J] = NULL)
{
Cout <"error" <Endl;
For (int K = 0; k <= J; k ++)
{
Delete [] Ders [k];
}
Delete [] Ders;
}
}*/

NDU [0] [0] = 1.0;
Double * Left = new double [p + 1];
Double * Right = new double [p + 1];

For (Int J = 1; j <= P; j ++)
{
Left [J] = u-u [I + 1-J];
Right [J] = U [I + J]-U;
Saves = 0.0;
For (INT r = 0; r <j; r ++)
{
// Lower triangle
NDU [J] [r] = right [R + 1] + left [J-R];
Temp = NDU [r] [J-1]/NDU [J] [r];

// Upper triangle
NDU [r] [J] = saved + right [R + 1] * temp;
Saved = left [J-R] * temp;
}
NDU [J] [J] = saved;
}

For (Int J = 0; j <= P; j ++)
{
// Load the basis functions
Ders [0] [J] = NDU [J] [p];
}

Int S1, S2, rk, PK, J1, J2;
Double D;
// This section computes the derivatives (Eq. [2.9])
For (INT r = 0; r <= P; r ++)
{
S1 = 0; S2 = 1; // alternate rows in array
A [0] [0] = 1.0;
// Loop to compute K' th Derivative
For (int K = 1; k <= N; k ++)
{
D = 0.0;
Rk = r-k;
PK = p-K;
If (r> = K)
{
A [s2] [0] = A [S1] [0]/NDU [PK + 1] [rk];
D = A [s2] [0] * NDU [rk] [PK];
}
If (rk> =-1)
J1 = 1;
Else
J1 =-rk;
If (R-1 <= PK)
J2 = k-1;
Else
J2 = P-r;
For (Int J = J1; j <= J2; j ++)
{
A [s2] [J] = (a [S1] [J]-A [S1] [J-1])/NDU [PK + 1] [rk + J];
D + = A [s2] [J] * NDU [rk + J] [PK];
}
If (r <= PK)
{
A [s2] [k] =-A [S1] [k-1]/NDU [PK + 1] [r];
D + = A [s2] [k] * NDU [r] [PK];
}
Ders [k] [r] = D;
// Switch rows
Int J = S1;
S1 = S2;
S2 = J;
}
}
// Multiply through by the correct factors
// (Eq.2.9)
Int r = P;
For (int K = 1; k <= N; k ++)
{
For (Int J = 0; j <= P; j ++)
{
Ders [k] [J] * = R;
}
R * = (p-k );
}
}
Int _ tmain (INT argc, _ tchar * argv [])
{
Int I0 = 0;
Int p0 = 2;
Double U0 [] = {0, 0, 1, 2, 3, 4, 5, 5 };
Double U0 = 5.0/2.0;
Int num0 = sizeof (U0)/sizeof (double)-P0-1;
I0 = findspan (num0, P0, U0, U0 );
Double * N0 = new double [P0 + 1];
Basisfuns (I0, U0, P0, U0, N0 );
Int nders = 2;
Double ** ders0;

Ders0 = new double * [nders + 1];
If (ders0 = NULL)
{
Cout <"error" <Endl;
}

For (Int J = 0; j <nders + 1; j ++)
{
Ders0 [J] = new double [P0 + 1];
If (ders0 [J] = NULL)
{
Cout <"error" <Endl;
For (int K = 0; k <= J; k ++)
{
Delete [] ders0 [k];
}
Delete [] ders0;
}
}

Dersbasisfuns (I0, U0, P0, nders, U0, ders0 );
For (INT I = 0; I <nders + 1; I ++)
{
For (Int J = 0; j <P0 + 1; j ++)
{
Double derivative = ders0 [I] [J];
Cout <ders0 [I] [J] <"\ t ";
}
Cout <Endl;
}
Return 0;
}