#declare surface =
  mesh
{
#declare iu      = (umax-umin)/uiter      //increment per step
#declare iv      = (vmax-vmin)/viter
#declare counter = 0
#declare vv      = vmin
  
#while (vv<=vmax)		//outer loop
  #declare uu = umin
  #while (uu<=umax)               //inner loop

//POINTS
    #declare x1=#include "uuvv.x"   //Calculate the
    #declare y1=#include "uuvv.y"   //four points in
    #declare z1=#include "uuvv.z"   //the rectangle defined
    #declare uu=uu+iu               //by uu,uu+iu,vv,vv+iv
    #declare x2=#include "uuvv.x"   //Uses uu and vv because
    #declare y2=#include "uuvv.y"   //u and v are built in
    #declare z2=#include "uuvv.z"   //constants.
    #declare vv=vv+iv
    #declare x3=#include "uuvv.x"
    #declare y3=#include "uuvv.y"
    #declare z3=#include "uuvv.z"
    #declare uu=uu-iu
    #declare x4=#include "uuvv.x"
    #declare y4=#include "uuvv.y"
    #declare z4=#include "uuvv.z"
    #declare vv=vv-iv
//NORMALS   used for smoothing triangles
    #declare vv=vv-iv               //Calculate surface normals
    #declare nx1=#include "uuvv.x"  //to the nine rectangles in
    #declare ny1=#include "uuvv.y"  //the immediate neighbourhood
    #declare nz1=#include "uuvv.z"
    #declare uu=uu+iu
    #declare nx2=#include "uuvv.x"
    #declare ny2=#include "uuvv.y"
    #declare nz2=#include "uuvv.z"
    #declare uu=uu+iu
    #declare vv=vv+iv
    #declare nx3=#include "uuvv.x"
    #declare ny3=#include "uuvv.y"
    #declare nz3=#include "uuvv.z"
    #declare vv=vv+iv
    #declare nx4=#include "uuvv.x"
    #declare ny4=#include "uuvv.y"
    #declare nz4=#include "uuvv.z"
    #declare vv=vv+iv
    #declare uu=uu-iu
    #declare nx5=#include "uuvv.x"
    #declare ny5=#include "uuvv.y"
    #declare nz5=#include "uuvv.z"
    #declare uu=uu-iu
    #declare nx6=#include "uuvv.x"
    #declare ny6=#include "uuvv.y"
    #declare nz6=#include "uuvv.z"
    #declare uu=uu-iu
    #declare vv=vv-iv
    #declare nx7=#include "uuvv.x"
    #declare ny7=#include "uuvv.y"
    #declare nz7=#include "uuvv.z"
    #declare vv=vv-iv
    #declare nx8=#include "uuvv.x"
    #declare ny8=#include "uuvv.y"
    #declare nz8=#include "uuvv.z"
    #declare uu=uu+iu
//Here is where normals are calculated.
    #declare n1=vcross(<x1-nx1,y1-ny1,z1-nz1>,<x1-nx8,y1-ny8,z1-nz8>)
    #declare n2=vcross(<x1- x2,y1- y2,z1- z2>,<x1-nx1,y1-ny1,z1-nz1>)
    #declare n3=vcross(<x1- x4,y1- y4,z1- z4>,<x1- x2,y1- y2,z1- z2>)
    #declare n4=vcross(<x1-nx8,y1-ny8,z1-nz8>,<x1- x4,y1- y4,z1- z4>)
    #declare normal1=n1+n2+n3+n4
    #declare n1=vcross(<x2-nx2,y2-ny2,z2-nz2>,<x2- x1,y2- y1,z2- z1>)
    #declare n2=vcross(<x2-nx3,y2-ny3,z2-nz3>,<x2-nx2,y2-ny2,z2-nz2>)
    #declare n3=vcross(<x2- x3,y2- y3,z2- z3>,<x2-nx3,y2-ny3,z2-nz3>)
    #declare n4=vcross(<x2- x1,y2- y1,z2- z1>,<x2- x3,y2- y3,z2- z3>)
    #declare normal2=n1+n2+n3+n4
    #declare n1=vcross(<x3- x2,y3- y2,z3- z2>,<x3- x4,y3- y4,z3- z4>)
    #declare n2=vcross(<x3-nx4,y3-ny4,z3-nz4>,<x3- x2,y3- y2,z3- z2>)
    #declare n3=vcross(<x3-nx5,y3-ny5,z3-nz5>,<x3-nx4,y3-ny4,z3-nz4>)
    #declare n4=vcross(<x3- x4,y3- y4,z3- z4>,<x3-nx5,y3-ny5,z3-nz5>)
    #declare normal3=n1+n2+n3+n4
    #declare n1=vcross(<x4- x1,y4- y1,z4- z1>,<x4-nx7,y4-ny7,z4-nz7>)
    #declare n2=vcross(<x4- x3,y4- y3,z4- z3>,<x4- x1,y4- y1,z4- z1>)
    #declare n3=vcross(<x4-nx6,y4-ny6,z4-nz6>,<x4- x3,y4- y3,z4- z3>)
    #declare n4=vcross(<x4-nx7,y4-ny7,z4-nz7>,<x4-nx6,y4-ny6,z4-nz6>)
    #declare normal4=n1+n2+n3+n4
//TRIANGLES
    smooth_triangle  {
      <x1, y1, z1 >, normal1,
      <x2, y2, z2 >, normal2,
      <x3, y3, z3 >, normal3
      }

    smooth_triangle  {
      <x1, y1, z1 >, normal1,
      <x3, y3, z3 >, normal3,
      <x4, y4, z4 >, normal4
      }

#if (tofile > 0)

    #debug concat("smooth_triangle{<",
                               str(x1,0,-1),",",str(y1,0,-1),",",str(z1,0,-1),">,<",
                               str(normal1.x,0,-1),",",str(normal1.y,0,-1),",",str(normal1.z,0,-1),">,<",
                               str(x2,0,-1),",",str(y2,0,-1),",",str(z2,0,-1),">,<",
                               str(normal2.x,0,-1),",",str(normal2.y,0,-1),",",str(normal2.z,0,-1),">,<",
                               str(x3,0,-1),",",str(y3,0,-1),",",str(z3,0,-1),">,<",
                               str(normal3.x,0,-1),",",str(normal3.y,0,-1),",",str(normal3.z,0,-1),
                               ">}","\n")
    #debug concat("smooth_triangle{<",
                               str(x1,0,-1),",",str(y1,0,-1),",",str(z1,0,-1),">,<",
                               str(normal1.x,0,-1),",",str(normal1.y,0,-1),",",str(normal1.z,0,-1),">,<",
                               str(x2,0,-1),",",str(y2,0,-1),",",str(z2,0,-1),">,<",
                               str(normal3.x,0,-1),",",str(normal3.y,0,-1),",",str(normal3.z,0,-1),">,<",
                               str(x3,0,-1),",",str(y3,0,-1),",",str(z3,0,-1),">,<",
                               str(normal4.x,0,-1),",",str(normal4.y,0,-1),",",str(normal4.z,0,-1),
                               ">}","\n")
#end




    #declare uu = uu+iu
  #end                            //inner loop
#if (tofile = 0)
  #debug concat("Looping on  ",str(counter,0,-1)," of ",str(viter,0,-1),"\n")
#end
  #declare counter = counter + 1
  #declare vv = vv+iv
#end				//outer loop
}				//mesh