/*
 * precompute bernstein factors and their derivates
 */
#declare subdiv = 16;	// subdivision level 16 means that 16 * 16 * 2 triangles will be created per patch

#declare B0 = array[subdiv + 1];
#declare B1 = array[subdiv + 1];
#declare B2 = array[subdiv + 1];
#declare B3 = array[subdiv + 1];
#declare dB0 = array[subdiv + 1];
#declare dB1 = array[subdiv + 1];
#declare dB2 = array[subdiv + 1];
#declare dB3 = array[subdiv + 1];
#declare i = 0;
#while ( i <= subdiv)
	#declare s = i / subdiv;
	#declare s1 = 1 - s;
	#declare s_2 = s * s;
	#declare s_3 = s_2 * s;
	#declare s1_2 = s1 * s1;
	#declare s1_3 = s1_2 * s1;
	#declare B0[i] = s1_3;
	#declare B1[i] = 3 * s * s1_2;
	#declare B2[i] = 3 * s_2 * s1;
	#declare B3[i] = s_3;
	#declare dB0[i] = -3 * s1_2;
	#declare dB1[i] = -6 * s * s1 + 3 * s1_2;
	#declare dB2[i] = -3 * s_2 + 6 * s * s1;
	#declare dB3[i] = 3 * s_2;
	#declare i = i + 1;
#end

/*
 * evaluate the bezier patch and build a mesh2 object
 */
#macro evalPatch(b)
	#local p = array[subdiv + 1][subdiv + 1];
	#local n = array[subdiv + 1][subdiv + 1];
	#local S = 0;
	#while (S <= subdiv)
		#local T = 0;
		#while (T <= subdiv)
			#local p0 = b[0] * B0[T] + b[1] * B1[T] + b[2] * B2[T] + b[3] * B3[T];
			#local p1 = b[4] * B0[T] + b[5] * B1[T] + b[6] * B2[T] + b[7] * B3[T];
			#local p2 = b[8] * B0[T] + b[9] * B1[T] + b[10] * B2[T] + b[11] * B3[T];
			#local p3 = b[12] * B0[T] + b[13] * B1[T] + b[14] * B2[T] + b[15] * B3[T];
			#local p[S][T] = p0 * B0[S] + p1 * B1[S] + p2 * B2[S] + p3 * B3[S];
			#local ds = p0 * dB0[S] + p1 * dB1[S] + p2 * dB2[S] + p3 * dB3[S];
			#local dt = B0[S] * (b[0] * dB0[T] + b[1] * dB1[T] + b[2] * dB2[T] + b[3] * dB3[T]) +
				    B1[S] * (b[4] * dB0[T] + b[5] * dB1[T] + b[6] * dB2[T] + b[7] * dB3[T]) +
				    B2[S] * (b[8] * dB0[T] + b[9] * dB1[T] + b[10] * dB2[T] + b[11] * dB3[T]) +
				    B3[S] * (b[12] * dB0[T] + b[13] * dB1[T] + b[14] * dB2[T] + b[15] * dB3[T]);
			#local n[S][T] = vnormalize(vcross(ds, dt));
			#local T = T + 1;
		#end
		#local S = S + 1;
	#end
	#local N = (subdiv + 1) * (subdiv + 1);
	mesh2 {
		vertex_vectors {
			N,
	#local S = 0;
	#while (S <= subdiv)
		#local T = 0;
		#while (T <= subdiv)
			p[S][T],
			#local T = T + 1;
		#end
		#local S = S + 1;
	#end
		}
		normal_vectors {
			N,
	#local S = 0;
	#while (S <= subdiv)
		#local T = 0;
		#while (T <= subdiv)
			n[S][T],
			#local T = T + 1;
		#end
		#local S = S + 1;
	#end
	#local F = subdiv * subdiv * 2;
		}
		face_indices {
			F,
	#local S = 0;
	#local D = subdiv + 1;
	#while (S < subdiv)
		#local T = 0;
		#while (T < subdiv)
			#local i1 = T * D + S;
			#local i2 = (T + 1) * D + S;
			#local i3 = (T + 1) * D + S + 1;
			#local i4 = T * D + S + 1;
			<i1, i2, i3>,
			<i3, i4, i1>,
			#local T = T + 1;
		#end
		#local S = S + 1;
	#end
		}
		normal_indices {
			F,
	#local S = 0;
	#local D = subdiv + 1;
	#while (S < subdiv)
		#local T = 0;
		#while (T < subdiv)
			#local i1 = T * D + S;
			#local i2 = (T + 1) * D + S;
			#local i3 = (T + 1) * D + S + 1;
			#local i4 = T * D + S + 1;
			<i1, i2, i3>,
			<i3, i4, i1>,
			#local T = T + 1;
		#end
		#local S = S + 1;
	#end
		}
	}
#end

/* ================= START OF TEST SCENE ================= */

/* cam and light */
camera {
	location <1.5,8,-2>
	look_at <2,0,2>
}

light_source {
	<-100,200,0>
	color rgb 1
}

/* the 16 controlpoints of a bezier patch */
#declare patch = array[16] {
	<0 1,0>,<1, 0,0>,<2, 0,0>,<3, 1,0>,
	<0, 0,1>,<1, 1,1>,<2, 1,1>,<3, 0,1>,
	<0, 0,2>,<1, 2,2>,<2, 0,2>,<3, -2,2>,
	<0,-1,3>,<1, 1,3>,<2,-1,3>,<3,-1,3>
};

/* render the test patch */
object {
	evalPatch(patch)
	pigment { color rgb <1,0,0> }
	finish { ambient 0.5 diffuse 0.5 specular 1 roughness 0.01 }
}