# bezcurve3d.txt  - commands to compute bezier curves given matrix of
# control pts
printf("\nwith(plots): with(plottools): with(linalg):");
with(plots): with(plottools): with(linalg):

printf("\nCommands in bezcurve3d.txt");
printf("\n  BezCurve(matrix pts,var) - computes Bezier curve from matrix of control points in var");
printf("\n  BezPlot(curve,var,color) - plots curve in parameter var in the color given with axes");

# bezcurv - computes bezier curve given matrix of control pts
BezCurve:=proc(pts,var)
  local d,i,m,xp,yp,zp,x,y,z,ans:
  m:=rowdim(pts):
  d:=m-1:
  for i from 0 to d do
    xp[i]:=pts[i+1,1]:
    yp[i]:=pts[i+1,2]:
    zp[i]:=pts[i+1,3]:
  od:
  unassign('i'):
  x:=evalf(sum(binomial(d,i)*var^i*(1-var)^(d-i)*xp[i],i=0..d)):
  y:=evalf(sum(binomial(d,i)*var^i*(1-var)^(d-i)*yp[i],i=0..d)):
  z:=evalf(sum(binomial(d,i)*var^i*(1-var)^(d-i)*zp[i],i=0..d)):
  RETURN (vector([x,y,z]));
  end:
    
BezPlot:=proc(curve,var,col)
 spacecurve([curve[1],curve[2],curve[3]],var=0..1,scaling=constrained,axes=boxed,color=col,tickmarks=[0,0,0]);
 end: 

Basis:=proc(U,K)
  local m, N, i, k, j, dU1, dU2, d1, d2, P;
   
  m:=vectdim(U)-1:
  N:=array(0..m,0..m):
  for i from 0 to m-1 do
    N[i,0]:=Heaviside(t-U[i+1])-Heaviside(t-U[i+2]):
  od:
  unassign('i'):

  for k from 1 to K+1 do
    for i from 0 to m-k-1 do
      j:=i+1:
      dU1:=U[j+k]-U[j]:
      dU2:=U[j+k+1]-U[j+1]:
      d1:=0:
      d2:=0:
      if dU1 > 0 then d1:=(t-U[j])/dU1 fi:
      if dU2 > 0 then d2:=(U[j+k+1]-t)/dU2 fi:
      N[i,k]:=d1*N[i,k-1]+d2*N[i+1,k-1]:
    od:
  od:
  unassign('i','j','k'):
 
  RETURN (N);
  end:

BSpline3d:=proc(pts,U)
  local n, m, i, k, x, y, z, N, A;
 
  n:=rowdim(pts)-1:
  m:=vectdim(U)-1:
  k:=n+1-m:
  N:=Basis(U,k):
  x:=sum(N[i,k]*pts[i+1,1],i=0..n):
  y:=sum(N[i,k]*pts[i+1,2],i=0..n):
  z:=sum(N[i,k]*pts[i+1,3],i=0..n):
  RETURN ([x,y,z]);
  end: