// Last edited on 2011-09-12 02:17:20 by stolfilocal

#macro bzpatch (P, txp, rad, txg)
  // Bezier patch.
  // P = control points (a 4x4 array)
  // txp = patch texture.
  // rad = radius of grid rods (0 = no grid).
  // txg = texture of grid.

  #macro grid_rod(p,q,rrad)
    // Rod between {p} and {q}, unless too close to be visible:
    #if (vlength(p - q) > 1.5*rrad) 
       cylinder{ p,q,rrad texture{ txg } }
    #end
  #end
  
  #macro grid_ball(p,srad)
    sphere{ P[i][j],srad texture{ txg } }
  #end

  union{

    bicubic_patch{
      type 1 u_steps 3 v_steps 3
      P[3][0],P[3][1],P[3][2],P[3][3]
      P[2][0],P[2][1],P[2][2],P[2][3]
      P[1][0],P[1][1],P[1][2],P[1][3]
      P[0][0],P[0][1],P[0][2],P[0][3] 
      texture {txp}
    }

    #local srad = 2*rad;  // Normal ball radius.
    #local rrad = rad;    // Normal rod radius.
    
    #if (rad > 0)
      #local i = 0;
      #while (i < 4)
        #local j = 0;
        #while (j < 4)
          #local smag = 1.0;
          #local jmag = 1.0;
          #local imag = 1.0;
          #if (((i = 1) | (i = 2)) & ((j = 1) | (j = 2)))
            #local smag = 1.5;
          #end
          grid_ball(P[i][j],srad*smag)
          #if (i < 3)
            #if ((i = 1) & ((j = 1) | (j = 2)))
              #local imag = 2.0;
            #end
            grid_rod(P[i][j],P[i+1][j],rrad*imag)
          #end
          #if (j < 3)
            #if ((j = 1) & ((i = 1) | (i = 2)))
              #local jmag = 2.0;
            #end
            grid_rod(P[i][j],P[i][j+1],rrad*jmag)
          #end
          #local j = j+1;
        #end
        #local i = i+1;
      #end
    #end
  }

#end

#macro bzassign(A,B)
  // Assigns {A[i][j] = B[i][j]} for all {i,j} in {0..3}.
  #local i = 0;
  #while (i < 4)
    #local j = 0;
    #while (j < 4)
      #declare A[i][j] = B[i][j];
      #local j = j + 1;
    #end
    #local i = i + 1;
  #end
#end

#macro bzpoint(P,uu,vv)
  // Point on Bézier patch with local coordinates {uu,vv}.
  // If both {uu,vv} are 0 or 1, the result is {P[uu][vv]}.
  #local Q = array[4][4];
  bzassign(Q,P)
  // DeCasteljau on each row od {Q}:
  #local k = 3;
  #while (k > 0)
    #local i = 0;
    #while (i < 4)
      #local j = 0;
      #while (j < k)
        #local Q[i][j] = (1-vv)*Q[i][j] + vv*Q[i][j+1];
        #local j = j + 1;
      #end
      #local i = i+ 1;
    #end
    #local k = k - 1;
  #end
  // DeCasteljau on first column of {Q}:
  #local k = 3;
  #while (k > 0)
    #local i = 0;
    #while (i < k)
      #local Q[i][0] = (1-uu)*Q[i][0] + uu*Q[i+1][0];
      #local i = i+ 1;
    #end
    #local k = k - 1;
  #end
  // Return the result:
  #local res = Q[0][0];
  res
#end