(*Display Proceedures*)

uvPPO3[surfaces_, a_, b_, c_] := 
  param[surfaces, {u, 0, 1}, {v, 0, 1}, AspectRatio -> 1, 
      ViewPoint -> {a, b, c}, Axes -> True, PlotPoints -> 40] /. 
    param -> ParametricPlot3D


dispobj[obj_, a_, b_, c_] := 
  Show[Graphics3D[obj], Lighting -> False, AspectRatio -> Automatic, 
    ViewPoint -> {a, b, c}, PlotRange -> All, Axes -> True]


(*-----4 masses adding to 1 for 0 ² u ² 1 and 0 ² v ² 1---- -*)

XX = {(1 - v)(1 - u), (1 - v) u, v u, v(1 - u)}

(*-----4 points - 
    and the surface generated by the center of gravity----*)(*-----of the 4 \
masses above--*)

P1 = {1, 0, 0}

P2 = {1, 1, 0}

P3 = {0, 1, 1}

P4 = {0, 0, 1}

uvPPO3[XX[[1]] P1 + XX[[2]]P2 + XX[[3]]P3 + XX[[4]]P4, 2, 2, 3]

(* ----    Random Number generator  ---- -*)

ra[a_] := 1 + Floor[a Random[]]

(* ----    constructing a generic  n x n matrix ---- -*)

mkmat[n_] := Table[a[i, j], {i, 1, n}, {j, 1, n}]

A4 = mkmat[4]

MF[M_] := MatrixForm[M]

MF[A4]

(* ----    constructing an (n - 1)x(n - 1) , 
  matrix from an n x n matrix ---- -*)

quadra[M_, u_, v_, w_, z_] := Block[{i, j, ut, n}, (n = Length[M];
      ut = 
        Table[u M[[i]][[j]] + v M[[i]][[j + 1]] + w M[[i + 1]][[j + 1]] + 
            z M[[i + 1]][[j]], {i, 1, n - 1}, {j, 1, n - 1}];
      Return[ut])]

B3 = quadra[A4, u, v, w, z]

B3[[1]][[1]]

B3[[3]][[2]]

B3[[2]][[3]]

C2 = quadra[B3, u, v, w, z]

C2[[2]][[2]]

Expand[C2[[1]][[1]]]

DD = Expand[quadra[C2, u, v, w, z]]

(* ---- Computing the coefficients of the control points in the final point ----*)


givecoes[PT_, n_] := Block[{i, M},
    (M = Table[Coefficient[PT, a[i, j]], {i, 1, n}, {j, 1, n}];
      Return[MatrixForm[M]])]

givecoes[DD[[1]][[1]], 4]

(* ---- The final point gives the equyation of the bezier surface ----*)

C2 =  {{c[1, 1], c[1, 2]}, {c[2, 1], c[2, 2]}}

FinPt[C_, u_, v_, w_, z_] := 
  u C[[1]][[1]] + v C[[1]][[2]] + w C[[2]][[2]] + z C[[2]][[1]]

FinPt [C2, u, v, w, z]

(* ---- a typical set of 16 control points  ----*)

boo = {{{1, 1, 5}, {1, 2, 3}, {1, 3, 3}, {1, 4, 5}}, {{2, 1, 3}, {2, 2, 
        1}, {2, 3, 1}, {2, 4, 3}}, {{3, 1, 3}, {3, 2, 1}, {3, 3, 1}, {3, 4, 
        3}}, {{4, 1, 5}, {4, 2, 3}, {4, 3, 3}, {4, 4, 5}}}

MF[%]

(* ----  Constructing the parametric equation of the bezier surface  ----*)
   
(* ----         with Displays of the intermediate matrices  ----*)

DBezSurf[M_, u_, v_, w_, z_] := Block[{i, prev, new}, (n = Length[M];
      S[1] = M; Print[MF[S[1]]];
      prev = M;
      Do[(new = quadra[prev, u, v, w, z];
          prev = new;
          S[i] = Expand[new];
          Print[i]; Print[MF[new]]), {i, 1, n - 2}];
      S[n] = Expand[FinPt[prev, u, v, w, z]];
      Print[MF[S[n]]];
      Return[Expand[S[n]]])]

DBezSurf[boo, .3, .1, .4, .2]

(* ----  Constructing the parametric equation of the bezier surface  ----*)

(* ----         without Displays of the intermediate matrices  ----*)

(* ---- this is the proceedure you use to construct your desired surface ----*)

BezSurf[M_, u_, v_, w_, z_] := Block[{i, prv, new}, (n = Length[M];
      S[1] = M;
      prev = M;
      Do[(new = quadra[prev, u, v, w, z];
          prev = new;
          S[i] = new;), {i, 1, n - 2}];
      S[n] = FinPt[prev, u, v, w, z];
      Return[Expand[S[n]]])]

now = BezSurf[boo, v u, v(1 - u), (1 - v)(1 - u), (1 - v)u]

see = uvPPO3[now, 3, 3, 4]

soo = {{{1, 1, 8}, {1, 2, 3}, {1, 3, 3}, {1, 4, 8}}, {{2, 1, 3}, {2, 2, 
        10}, {2, 3, 10}, {2, 4, 3}}, {{3, 1, 3}, {3, 2, 10}, {3, 3, 10}, {3, 
        4, 3}}, {{4, 1, 8}, {4, 2, 3}, {4, 3, 3}, {4, 4, 8}}}


wow = BezSurf[soo, v u, v(1 - u), (1 - v)(1 - u), (1 - v)u]

wee = uvPPO3[wow, 3, 3, 4]

coo = {{{1, 1, 5}, {1, 2, 5}, {1, 3, 5}, {1, 4, 5}},
    {{2, 1, 3}, {2, 2, 8}, {2, 3, 8}, {2, 4, 3}},
    {{3, 1, 3}, {3, 2, 8}, {3, 3, 8}, {3, 4, 3}},
    {{4, 1, 5}, {4, 2, 5}, {4, 3, 5}, {4, 4, 5}}}

cow = BezSurf[coo, v u, v(1 - u), (1 - v)(1 - u), (1 - v)u]

cee = uvPPO3[cow, 3, 3, 4]

(*--- translation Proceedure  -- -*)

trans[surf_, a_, b_, c_] := surf + {a, b, c}

trans[{2, 3, 4}, a, b, c]

(*--- experimenting  -- -*)

tcow = trans[cow, 0, 2, 0]


tcee = uvPPO3[tcow, 3, 3, 4]

Show[{cee, tcee}]

(* --- showing the control squares  -- -*)

squares[M_] := Block[{i, j, ut, n}, (n = Length[M];
      ut = 
        Table[Line[{M[[i]][[j]], M[[i]][[j + 1]], M[[i + 1]][[j + 1]], 
              M[[i + 1]][[j]], M[[i]][[j]]}], {i, 1, n - 1}, {j, 1, n - 1}];
      Return[ut])]


orsq = squares[boo]

pua = dispobj[orsq, 2, 2, 4]

Show[see, pua]


(* --- procedures to illustrate the motion of the point when we vary the weights ---*)


(*--- making random control points  ----*)

mkrasurf[n_, a_] := Table[Table[{ra[a], ra[a], ra[a]}, {n}], {n}]


rapts = mkrasurf[4, 8]


(*--- making random surface from the random control points  ----*)


raeq = BezSurf[rapts, v u, v(1 - u), (1 - v)(1 - u), (1 - v)u]


uvPPO3[raeq, -9, 3, 13]

(* --- procedures to illustrate the motion of the point when we vary the weights -- -*)


squares[M_] := Block[{i, j, ut, n}, (n = Length[M];
      ut = 
        Table[Line[{M[[i]][[j]], M[[i]][[j + 1]], M[[i + 1]][[j + 1]], 
              M[[i + 1]][[j]], M[[i]][[j]]}], {i, 1, n - 1}, {j, 1, n - 1}];
      Return[ut])]


 polygs[M_, c_] := Block[{i, j, ut, n}, (n = Length[M];
      ut = 
        Table[{Hue[c], 
            Polygon[{M[[i]][[j]], M[[i]][[j + 1]], M[[i + 1]][[j + 1]], 
                M[[i + 1]][[j]], M[[i]][[j]]}]}, {i, 1, n - 1}, {j, 1, 
            n - 1}];
      Return[ut])]


pol = polygs[boo, .3]


voo = dispobj[pol, 2, 2, 4]

Show[{see, voo}]

(*---- constructing a set of weights from four arbitrary positive  numbers  ----*)


gtweights[A_, B_, C_, D_] := {A, B, C, D}/(A + B + C + D)

(*---- proceedure for constructing a graphic object consisting of a colored point *)

pt[P_, c_, s_] := {Hue[c], PointSize[s], Point[P]}

(*---- proceedure for constructing a table of colored points *)


gpts[M_, c_, s_] := 
  Table[pt[M[[i]][[j]], c, s], {i, 1, Length[M]}, {j, 1, Length[M]}]


(*---- making the successive control matrices ----*)

mkpolos[Pts_, u_, v_] := Block[{n, ut, i, M}, (n = Length[Pts];
      prev = Pts;
      M[1] = Pts;
      Do[(new = quadra[prev, v u, v(1 - u), (1 - v)(1 - u), (1 - v)u];
          M[i] = new;
          prev = new), {i, 2, n}];
      M[n] = prev;
      Return[Table[M[i], {i, 1, n}]])]


baba = mkpolos[boo, .2, .5]

Map[MF, %]


(*--- converting the control matrices into control Polygons -- -*)


cvrt[seq_] := Table[polygs[seq[[i]], i/10], {i, 1, Length[seq] - 1}]


cvrt[baba]

dispobj[%, 2, 3, 7]

(*--- getting the corners of the control polygons ----*)

P4=gpts[baba[[1]],.8,.04]

cpp = dispobj[P4, 1, 2, 8]


pol = polygs[boo, .3]
 
voo=dispobj[pol,2,2,10]

Show[{voo, cpp}]

(* getting the displays of the various stages ----*)

stage[M_, c1_, c2_, s_] := {polygs[M, c1], gpts[M, c2, s]}

R1 = stage[baba[[1]], .2, .8, .04];


R2 = stage[baba[[2]], .4, .9, .03];


R3 = stage[baba[[3]], .5, .99, .02];


dispobj[{R1, R2, R3}, 1, 2, 12]


(*---- Showing the mechanism that produces a finalpoint on the surface ---- *)


mkshow[PTS_, a_, b_, c_, d_, l_, m_, n_] := 
  Block[{u, v, w, z, R1, R2, R3, R4, PP, BB, CC, 
      DD}, ({u, v, w, z} = gtweights[a, b, c, d];
      BB = quadra[PTS, u, v, w, z];
      CC = quadra[BB, u, v, w, z];
      PP = quadra[CC, u, v, w, z];
      R1 = stage[PTS, .1, .6, .04];
      R2 = stage[BB, .4, .8, .03];
      R3 = stage[CC, .5, .1, .03];
      R4 = pt[Flatten[PP], 0, .04];
      dispobj[{R1, R2, R3, R4}, l, m, n])]

mkshow[boo, ra[5], ra[5], ra[5], ra[5], 1, 2, 9]

mkshow[boo, ra[5], ra[5], ra[5], ra[5], 1, 2, 9]

mkshow[boo, ra[5], ra[5], ra[5], ra[5], 1, 2, 9]

mkshow[boo, ra[5], ra[5], ra[5], ra[5], 1, 2, 9]

mkshow[boo, ra[5], ra[5], ra[5], ra[5], 1, 2, 9]