printf("\n bpatch.txt routines for creating and analyzing Bezier patches");
printf("\nwith(plots): with(plottools): with(linalg):");
with(plots): with(plottools): with(linalg):

printf("Commands in patch.txt");
printf("\n  BPatch(integer m, integer n, matrix pts) - computes a Bezier patch from a matrix of control points");
printf("\n  PatchPlot(vector patch) - plots a Bezier patch with an isoparametric net");
printf("\n  ContNet(integer m, integer n, matrix pts) - plots the control net of a bicubic Bezier patch");
printf("\n           from a matrix of control points");
#printf("\n  Deriv1(real u, real v, vector patch) - plots the first derivative vectors for the patch");
#printf("\n           at the point corresponding to (u,v)");
#printf("\n  Deriv2(real u, real v,vector patch) - plots the second derivative vectors for the patch");
#printf("\n           at the point corresponding to (u,v)");

BPatch:=proc (m,n,pts)
  local i,j,k,XCrd,YCrd,ZCrd,X,Y,Z,Bi,Bj,ans;

  XCrd:= array(0..m,0..n):
  YCrd:= array(0..m,0..n):
  ZCrd:= array(0..m,0..n):
  for i from 0 to m do
    for j from 0 to n do
      k:=1+i+(m+1)*j:
      XCrd[i,j]:=pts[k,1]:
      YCrd[i,j]:=pts[k,2]:
      ZCrd[i,j]:=pts[k,3]:
    od:
  od:
  unassign('i','j','k'):

  X:=simplify(sum(binomial(m,i)*u^i*(1-u)^(m-i)*sum(binomial(n,j)*v^j*(1-v)^(n-j)*XCrd[i,j],j=0..n),i=0..m)):
  Y:=simplify(sum(binomial(m,i)*u^i*(1-u)^(m-i)*sum(binomial(n,j)*v^j*(1-v)^(n-j)*YCrd[i,j],j=0..n),i=0..m)):
  Z:=simplify(sum(binomial(m,i)*u^i*(1-u)^(m-i)*sum(binomial(n,j)*v^j*(1-v)^(n-j)*ZCrd[i,j],j=0..n),i=0..m)):

  ans:=vector([X,Y,Z]):
  RETURN (evalm(ans))
end:
 
PatchPlot:=proc(patch)
   plot3d(patch,u=0..1,v=0..1,scaling=constrained,axes=boxed,tickmarks=[0,0,0]):
end:

ContNet:=proc(m,n,pts)
  local i,j,k,b,p,q,A,B,C:

  b:=array(0..m,0..n):
  for i from 0 to m do
    for j from 0 to n do
      k:=1+i+(m+1)*j:
      b[i,j]:=[pts[k,1],pts[k,2],pts[k,3]]:
    od:
  od:
  unassign('i','j'):
  for i from 0 to m do
    p[i]:=curve([seq(b[i,j],j=0..n)],color=black):
  od:
  unassign('i'):
  for j from 0 to n do
    q[j]:=curve([seq(b[i,j],i=0..m)],color=black):
  od:
  unassign('j'):
  A:=display([seq(p[i],i=0..m)],scaling=constrained,axes=none,thickness=1,tickmarks=[0,0,0]):
  B:=display([seq(q[j],j=0..n)],scaling=constrained,axes=none,thickness=1,tickmarks=[0,0,0]):
  C:=pointplot3d(pts,thickness=4,color=red,tickmarks=[0,0,0]):

  RETURN (display({A,B,C}));
end:

Deriv1:=proc(u0,v0,patch)
  local du,dv,p0,p1,p2,L1,L2:

  du:=map(diff,patch,u):
  dv:=map(diff,patch,v):

  p0:=subs({u=u0,v=v0},evalm(patch)):
  p1:=subs({u=u0,v=v0},evalm(patch+(1/3)*du)):
  p2:=subs({u=u0,v=v0},evalm(patch+(1/3)*dv)):

  L1:=line([p0[1],p0[2],p0[3]],[p1[1],p1[2],p1[3]],color=black,thickness=3):
  L2:=line([p0[1],p0[2],p0[3]],[p2[1],p2[2],p2[3]],color=black,thickness=3):
  RETURN (display(L1,L2));

end:

Deriv2:=proc(u0,v0,patch)
  local d2u,d2v,dudv,p0,p1,p2,p3,L1,L2,L3:

  d2u:=[diff(patch[1],u$2),diff(patch[2],u$2),diff(patch[3],u$2)]:
  d2v:=[diff(patch[1],u$2),diff(patch[2],u$2),diff(patch[3],v$2)]:
  dudv:=[diff(patch[1],u,v),diff(patch[2],u,v),diff(patch[3],u,v)]:

  p0:=subs({u=u0,v=v0},evalm(patch)):
  p1:=subs({u=u0,v=v0},evalm(patch+(1/10)*d2u)):
  p2:=subs({u=u0,v=v0},evalm(patch+(1/10)*d2v)):
  p3:=subs({u=u0,v=v0},evalm(patch+(1/10)*dudv)):

  L1:=line([p0[1],p0[2],p0[3]],[p1[1],p1[2],p1[3]],color=blue,thickness=3):
  L2:=line([p0[1],p0[2],p0[3]],[p2[1],p2[2],p2[3]],color=green,thickness=3):
  L3:=line([p0[1],p0[2],p0[3]],[p3[1],p3[2],p3[3]],color=red,thickness=3):

  RETURN (display(L1,L2,L3));
end: